Computer Lab 9

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 FRM
source("functions10A.r")
source("usastocks.r")

### usastocks.r has a combined dataset for SP500 and Nasdaq. A key issue
# in compiling this file was to ensure that all of the trading days coincide.


########## Exercise 1:
# Here we will show two ways to simulate a bivariate dataset with
# Normal(0,1) marginals but different dependency structures.

# (a) (Bivariate normal distribution)

 rho=0.8
 z1=rnorm(10000)
 z2=rho*z1+sqrt(1-rho^2)*rnorm(10000)
	# the last two lines mean that both z1 and z2 are N(0,1)
	# (z1,z2) is bivariate normal with cor(z1,z2)=rho
 setaxes(-4,4,-4,4)
 points(z1,z2,pch=20,cex=0.5)	# pch=20 means that we get solid dots
				# cex=0.5 means character expansion 0.5 times normal

# next do a QQnorm plot to verify graphically that the marginals are normal

 qqnorm(z1)
 qqnorm(z2)

# (b) Normal marginals using the t-3 copula
#     Keep z1 and z2 as above, but now use these in a different way
#     in order to create a standard bivariate t

 nu=3
 y=rchisq(10000,df=nu)	# generates 10000 i.i.d. chi-sqaured r.v.'s
			# with nu degrees of freedom
 x1=z1/sqrt(y/nu)	# this defines a standard t_nu random variable 
 x2=z2/sqrt(y/nu)
 u1=pt(x1,df=nu)	# u1 will be i.i.d. ~ U[0,1]
 u2=pt(x2,df=nu)

# first plot the empirical t-3 copula 

 setaxes(0,1,0,1)
 points(u1,u2,pch=20,cex=0.5)

# now generate the normal random variables

 w1=qnorm(u1)
 w2=qnorm(u2)
 setaxes(-4,4,-4,4)
 points(w1,w2,pch=20,cex=0.5)

# It helps to be able to compare the results of (a) and (b) side
# by side as follows:

 par(mfrow=c(1,2))	# mfrow=c(1,2) means set up a plotting region that
			# will have one row and two plots per row
 setaxes(-4,4,-4,4)
 points(w1,w2,pch=20,cex=0.5)
 setaxes(-4,4,-4,4)
 points(z1,z2,pch=20,cex=0.5)

# and then maximise the plotting window to view the plots on the full screen
# Alternatively, distort the graphics window until the two plots are reasonably square.

# Do these two scatter plots look different?
# Compare and contrast the two plots.
# If they do look different check that the marginals of (w1,w2) are normal


# What about w1-w2 ?
# What about z1-z2 ? (and what, theoretically is the distribution of z1-z2?)
# e.g. try

qqnorm(w1-w2)
qqnorm(z1-z2)

# (c) To see the impact of the dependency structure on risk measurement,
#     now calculate the mean excess gain of w1-w2 in excess of a threshold of 1

 dw=w1-w2
 dw=dw[order(-dw)]	# put dw into decreasing order
 umin=1			# set the minimum threshold
 n1=sum(dw > umin)	# count how many of w1-w2 exceed the threshold
 excess.gains=dw[1:n1]-umin	# calculate the conditional excess gains
 mean(excess.gains)	# calculate the mean

# Alternatively the following will achieve the same result:

 mean(dw[dw > umin]-umin)

# Compare this result with what you get for z1-z2




Exercise 2:

Investigate the distinguishing characteristics of the Clayton, Gumbel,
Gaussian and t copulas by using the following functions:

# Contour plots can be generated as follows

 plot.clayton.copula(theta=1.5,pl=1)
 plot.gumbel.copula(theta=2.3,pl=1)
 plot.gauss.copula(rho=0.5,pl=1)
 plot.t.copula(df=4.5,rho=0.7,pl=1)

# Perspective plots can be generated by changing the argument pl=1
# above to pl=2

 plot.clayton.copula(theta=1.5,pl=2)
 plot.gumbel.copula(theta=2.3,pl=2)
 plot.gauss.copula(rho=0.5,pl=2)
 plot.t.copula(df=4.5,rho=0.7,pl=2)

# Simulations can be obtained using the following functions

 copula.data=simulate.clayton.copula(theta=1.5,ns=10000)	
	# This function takes a little bit of time to run, so a counter is
	# output on the screen (which goes up to 1000) so that you can monitor
	# progress
	# ns is the number of simulations
	# the output from the function is an array that has two rows
	# each row of length ns. 
	# Row 1 is the simulated u1, Row 2 is the simulated u2.
	# We store this output in the array called copula.data
 setaxes(0,1,0,1)	# now do a scatterplot of (u1,u2)
 points(copula.data[1,],copula.data[2,],pch=20,cex=0.5)

 copula.data=simulate.gumbel.copula(theta=3,ns=10000)	
 setaxes(0,1,0,1)
 points(copula.data[1,],copula.data[2,],pch=20,cex=0.5)

 copula.data=simulate.gauss.copula(rho=0.7,ns=10000)	
 setaxes(0,1,0,1)
 points(copula.data[1,],copula.data[2,],pch=20,cex=0.5)

 copula.data=simulate.t.copula(df=4,rho=0.4,ns=10000)	
 setaxes(0,1,0,1)
 points(copula.data[1,],copula.data[2,],pch=20,cex=0.5)

# Exercise 2a:
# For the Clayton copula try out all three functions for theta=1.5.
# Can you see how the scatterplot of the simulated data matches the contour
# and perspective plots?
# Try different values for theta and see how the shape of the density changes.
# Remember that we must have theta > 0!
# 
# Exercise 2b-d:
# Repeat exercise 1a for the Gumbel (theta > 1), Gaussian and t copulas.
# 




########## Exercise 3:
###### Load up the file "lab9data.r" (copy it from the module web page first!)

source("lab9data.r")

# This file contains 4 simulated copulas named
# lab9data$va1
# lab9data$va2
# lab9data$va3
# lab9data$va4

To view all four copulas at the same time type:

 par1()
 par(mfrow=c(2,2))	# this sets us up to do 4 plots on one screen
 par(mar=c(4,4,1,1))	# This reduces the margin widths around the individual plots

 setaxes(0,1,0,1)
 points(lab9data$va1[1,],lab9data$va1[2,],pch=20,cex=0.3)

 setaxes(0,1,0,1)
 points(lab9data$va2[1,],lab9data$va2[2,],pch=20,cex=0.3)

 setaxes(0,1,0,1)
 points(lab9data$va3[1,],lab9data$va3[2,],pch=20,cex=0.3)

 setaxes(0,1,0,1)
 points(lab9data$va4[1,],lab9data$va4[2,],pch=20,cex=0.3)

# Exercise 3a:
# Confirm that these are copulas by verifying that the marginals 
# are approximately uniform[0,1]
# for example:

 hist(lab9data$va1[1,])

# Exercise 3b:
# The four copulas are were generated using the Clayton, Gumbel, 
# Gaussian and t copulas but they are in a different order now.
# Use what you have learned about the characteristics of the 
# four copulas in Exercise 1 to determine which simulated dataset
# was generated by which copula.



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

Now start to analyse the S&P500 and Nasdaq data

tv=usa.stocks$t1		# dates corredponding to the daily log returns d1 and d2 below
x1=usa.stocks$sp500
x2=usa.stocks$nasdaq
x1=x1/x1[1]*1000		# normalised indices that start at 1000 at the beginning of day 1
x2=x2/x2[1]*1000
d1=diff(log(x1))
d2=diff(log(x2))

res1=fit.garch11.normal(d1)   # Fit the GARCH(1,1) model using QML to SP500
res2=fit.garch11.normal(d2)   # Fit the GARCH(1,1) model using QML to Nasdaq

Z1=res1$Z
Z2=res2$Z

nz=length(Z1)
za=array(0,c(2,nz))
za[1,]=Z1
za[2,]=Z2

# Fit the NCT to the marginals for Z1 and Z2 independently

pv31=mle.nct(Z1)	# pv31=(mu,sigma,df,ncp) format
pv32=mle.nct(Z2)

# Now calculate the probability transformed random variables
u1=pt((Z1-pv31[1])/pv31[2],df=pv31[3],ncp=pv31[4])
u2=pt((Z2-pv32[1])/pv32[2],df=pv32[3],ncp=pv32[4])

# Verify that the NCT gives a good fit using a chi-squared test
chi2test.unif(u1,4)	# 4 means 4 constraints, corresponding to 4 parameters estimated
chi2test.unif(u2,4)

# Put u1 and u2 into an array that will be used as input for the copula fitting
# functions
ua=array(0,c(2,nz))
ua[1,]=u1
ua[2,]=u2

# Use MLE to estimate the copula parameters
cv1=mle.clayton.copula(ua)
cv2=mle.gumbel.copula(ua)
cv3=mle.gauss.copula(ua)
cv4=mle.t.copula(ua)

# In each case note down the log-likelihood and decide which model is best

###### Section 6.7 non-parametric alternative ######
# Rerun the above estimation routine using the non-parametric approach
# based on rank statistics

u1A=(rank(Z1)-0.5)/nz	# A for "Alternative" is added to the name of the vector etc.
u2A=(rank(Z2)-0.5)/nz

uaA=array(0,c(2,nz))
uaA[1,]=u1A
uaA[2,]=u2A

cv1A=mle.clayton.copula(uaA)
cv2A=mle.gumbel.copula(uaA)
cv3A=mle.gauss.copula(uaA)
cv4A=mle.t.copula(uaA)

# Are the parameter estimates about the same

###### Exercise 4A:

# Use the following commands to simulate (U1,U2) from the fitted
# copulas, then plot each one and compare with the market copula
# From above cv1 has the Clayton parameter estimate, cv2=Gumbel,
# cv3=Gaussian, cv4=t-copula rho and df.

# Market copula
 	setaxes(0,1,0,1)
 	points(ua[1,],ua[2,],pch=20,cex=0.5)

# Simulated copulas - you must fill in the blanks!

 copula.data=simulate.clayton.copula(theta=....,ns=nz)	
 	setaxes(0,1,0,1)
 	points(copula.data[1,],copula.data[2,],pch=20,cex=0.5)
 copula.data=simulate.gumbel.copula(theta=...,ns=nz)	
 copula.data=simulate.gauss.copula(rho=...,ns=nz)	
 copula.data=simulate.t.copula(df=...,rho=...,ns=nz)	

# In each case try to plot the simulated copula alongside the market copula