Computer Lab 5

Web page: http://www.ma.hw.ac.uk/~andrewc/erm

Go to the module web page and, if you have not already done so
copy the following files to your directory H:\FRM
functions10A.r
sp500.r
nasdaq.r

Load each one into your R workspace, if you have not already done so:
for example 

source("functions10A.r")
source("sp500.r")
source("nasdaq.r")

### Note some graphical results with discussion will appear later
### on the module web page.

x=sp500$adjclose
tv=sp500$t
n=length(x)
tv=tv[2:n]
d=diff(log(x))
nd=length(d)

# Note that tv is the same length as d

nv=(1:nd)[tv > 1980] 	# nv is a set of index numbers
tv=tv[nv]			# tv[nv] extracts the dates which are all after 1980
d=d[nv]			# This line extracts the daily log returns after 1980
nd=length(d)		# recalculate the length of d

plot(tv,d,pch=20)

# The next few lines take a while to run!
# They output the current log-likelihood on the screen so you can watch progress
# Don't worry about warning messages that appear at the end. These occur if the 
# optimisation algorithm steps outside the admissible region for parameters
# at any stage (e.g. it tries a value for nu that gives infinite variance.

res1=fit.garch11.normal(d)	# Fit GARCH(1,1) assuming Z(t)~N(0,1)
				# This is the same as the quasi maximum likelihood
				# which assumes the Z(t) are normal
res2=fit.garch11.t(d)		# Fit GARCH(1,1); Z(t)~t_nu with E[Z]=0, Var[Z]=1
res3=fit.garch11.nct(d)		# Fit GARCH(1,1); Z(t)~NCT with E[Z]=0, Var[Z]=1

# Now type 
names(res1)	# This gives a list of the different components of the
		# output stored in res1 etc.

# res1$sigma is the vector of volatilities (i.e. conditional standard deviations)
# res1$Z     is the vector of standardised residuals
# res1$ll    is the value of the maximum log-likelihood (useful for AIC, BIC)
# res1$x     is the original process (e.g. daily log returns) being modelled
# res1$par   is the vector of parameter values

n1=length(res1$x)
plot(1:n1,res1$x,pch=20,cex=0.5)	# plot the daily log returns
lines(1:n1,res1$sigma,col=2)		# plot the fitted volatility sigma(t)

names(res1$par)	# Lists the names of the parameter values stored in res1$par
names(res2$par)	# Lists the names of the parameter values stored in res2$par
names(res3$par)	# Lists the names of the parameter values stored in res3$par

pv1=res1$par	# return the parameter values for the Normal GARCH(1,1)
			# This is also the Quasi Maximum Likelihood output
			# The vector has 5 values 
			# sigma(1), alpha0, alpha1, beta1, mu

ZQ=res1$Z	# QML volatility standardised residuals
		# Also the true volatility standardised residuals if the 
		# Normal hypothesis is correct

####### Exercise 1 ######

# Now test for the following questions:
# Plot the ZQ vector
# Are the ZQ i.i.d.?
# Are the ZQ Normally distributed?
# Work out how to simulate the next day's daily log return as follows:
# Suppose that you are using the GARCH(1,1)-Normal model
# First work out sigma(n+1) given the values of sigma(n), alpha0,
# alpha1, beta1, and Z(n) from the res1 output.
# Then note that Z(n+1) is N(0,1).
# Hence delta(n+1), conditional of the history up to time n,
# is normal with mean mu and variance sigma(n+1)^2.
# Next do the same for the t distribution. But you will need to
# work out the number of degrees of freedom, the location and the
# scaling parameter for the three parameter t distribution
# for the conditional distribution of delta(n+1).
# Finally, do the same for the NCT. This is much more difficult!
# You will have to first work out the location and the scaling parameter
# parameters are for the NCT given that the Z(t) have mean 0 and  
# variance 1. You will need to (a) look up Chapter 4
# to get the formulae for E[X] and E[X^2], where X has a standard
# NCT distribution given the degrees of freedom and the non-centrality
# parameter NCP. Then work out constants a=a(nu,NCP) and b=b(nu,NCP)
# that result if E[a+bX]=0 and Var(a+bX)=1.
# Define Z(t)=a+b.X(t) where the X(t) are i.i.d. standard NCT(nu,NCP)'s.
# Finally, delta(n+1) = mu + sigma(n+1) Z(n+1)
#                     = mu + sigma(n+1) [ a + b.Z(n+1) ]
#                     = [mu + sigma(n+1).a] + [sigma(n+1).b].X(n+1)
# i.e. NCT with nu degrees of freedom, non-centrality parameter NCP
# location parameter [mu + sigma(n+1).a] and scaling parameter 
# [sigma(n+1).b]

###########################################################

pv2=res2$par	# return the parameter values for the t-GARCH(1,1)
			# The vector has 6 values 
			# sigma(1), alpha0, alpha1, beta1, mu, nu
pv3=res3$par	# return the parameter values for the t-GARCH(1,1)
			# The vector has 7 values 
			# sigma(1), alpha0, alpha1, beta1, mu, nu, ncp

###### Exercise 2 #######

# Compare the estimated sigma(t) volatility functions resulting from
# the three MLE functions (normal, t, and NCT assumptions)
# and stored in results files res1, res2, res3

setaxes(1980,2006,-0.08,0.08)
points(tv,d,pch=20,cex=0.5)

# Now overlay res1$sigma, res2$sigma, res3$sigma using different colours
# Are they quite similar? Where do the differences occur?

# Also try:

plot(res1$sigma,res2$sigma,pch=20,cex=0.5)
abline(0,1,col=2)
plot(res1$sigma,res3$sigma,pch=20,cex=0.5)
abline(0,1,col=2)
plot(res2$sigma,res3$sigma,pch=20,cex=0.5)
abline(0,1,col=2)

# In each case try to explain what you see by referring to the parameter 
# values pv1, pv2, pv3

# What do you conclude?

########## Exercise 2B ################

# Now compare the Z(t)'s

cdfplot(res1$Z)
cdfplot(res2$Z,a=1,co=2)
cdfplot(res3$Z,a=1,co=3)

# What do you conclude?

############ Exercise 2C #################

# Now compare the Z(t)'s by individual date

plot(res1$Z,res2$Z,pch=20,cex=0.5)
abline(0,1)

plot(res1$Z,res3$Z,pch=20,cex=0.5)
abline(0,1)

plot(res2$Z,res3$Z,pch=20,cex=0.5)
abline(0,1)

# What do you conclude?

######## Exercise 3 - Repeat all of the above using the Nasdaq data

x=nasdaq$adjclose
tv=nasdaq$t
n=length(x)
tv=tv[2:n]
d=diff(log(x))
	nd=length(d)
plot(tv,d,pch=20)

####### Exercise 4 ######

# Go back to the SP500 data
# The data go up to and including the daily log return delta(nd)

# Simulate 1000 independent realisations of delta(nd+1)

# HINT: Work out how to simulate the next day's daily log return.
# (i.e. use the values in res1$par and res1$sigma[n1]
# to calculate res1$sigma[n1+1].