*
* Example 10.6 from pp 472-477
*
open data d-spcscointc.txt
data(format=free,org=columns) 1 2275 sp500 cisco intel
*
* Sample statistics
*
table / sp500 cisco intel
vcv(center)
# sp500 cisco intel
*
* Multivariate Q statistic
*
@mvqstat(lags=1)
# sp500 cisco intel
@mvqstat(lags=4)
# sp500 cisco intel
@mvqstat(lags=8)
# sp500 cisco intel
*
@mvident(lags=6)
# sp500 cisco intel
*
garch(reg,p=1,q=1,hseries=h,resids=u) / sp500
# constant sp500{1 3 4 5}
*
set ustd = u/sqrt(h)
@regcorrs(dfc=4,number=10) ustd
set ustdsq = ustd**2
@regcorrs(dfc=2,number=10) ustdsq
*
* Bivariate GARCH model
*
equation sp500eq sp500
# constant sp500{3}
equation ciscoeq cisco
# constant sp500{2} cisco{2}
group bimodel sp500eq ciscoeq
*
garch(model=bimodel,p=1,q=1,pmethod=simplex,piters=10)
*
* Cholesky factor three variable model
*
nonlin(parmset=meanparms) c1 p3 c2 p21 p22 c3 p31 p33
*
* Basic residuals
*
frml a1f = sp500-c1-p3*sp500{3}
frml a2f = cisco-c2-p21*sp500{2}-p22*cisco{2}
frml a3f = intel-c3-p31*sp500{1}-p33*intel{1}
*
nlsystem(parmset=meanparms) / a1f a2f a3f
*
set a1t = a1f
set a2t = a2f
set a3t = a3f
set b2t = 0.0
set b3t = 0.0
*
compute cholesky=%decomp(%sigma)
set g11t = cholesky(1,1)**2
set g22t = cholesky(2,2)**2
set g33t = cholesky(3,3)**2
set q21t = cholesky(2,1)
set q31t = cholesky(3,1)
set q32t = cholesky(3,2)
*
nonlin(parmset=garchparms) a0 a1 a2 b0 b1 b2 d0 d1 d2 d5 t0 t1 t2 u0 u1 u2 w0 w1 w2
*
frml g11f = a0+a1*a1t{1}**2+a2*g11t{1}
frml q21f = t0+t1*q21t{1}  +t2*a2t{1}
frml b2f  = a2f-q21t*a1f
frml g22f = b0+b1*b2t{1}**2+b2*g22t{1}
frml q31f = u0+u1*q31t{1}  +u2*a3t{1}
frml q32f = w0+w1*q32t{1}  +w2*a2t{1}
frml b3f  = a3f-q31t*a1f-q32t*a2f
frml g33f = d0+d1*b3t{1}**2+d2*g33t{1}+d5*g22t{1}
*
frml garchln = a1t=a1f,g11t=g11f,q21t=q21f,b2t=b2f,g22t=g22f,q31t=q31f,q32t=q32f,b3t=b3f,g33t=g33f,$
   %logdensity(g11t,a1t)+%logdensity(g22t,b2t)+%logdensity(g33t,b3t)
*
compute a0 =  0.01  , a1 = 0.05 , a2 = 0.94
compute t0 =  0.28  , t1 = 0.82 , t2 = -0.035
compute b0 =  0.17  , b1 = 0.08 , b2 = 0.89
compute u0 =  0.04  , u1 = 0.97 , u2 = 0.01
compute w0 =  0.006 , w1 = 0.98 , w2 = 0.004
compute d0 =  1.38  , d1 = 0.06 , d2 = 0.64   , d5 = -0.027
*
maximize(parmset=meanparms+garchparms,method=bhhh,iterations=250) garchln 5 *
*
set fv1   = g11t
set resi1 = a1t/sqrt(fv1)
set residsq = resi1**2
cor(qstats,number=12,span=4) resi1
cor(qstats,number=12,span=4) residsq
*
set fv2   = g22t+q21t**2*g11t
set resi2 = a2t/sqrt(fv2)
set residsq = resi2**2
cor(qstats,number=12,span=4) resi2
cor(qstats,number=12,span=4) residsq
*
set fv3   = g33t+q31t**2*g11t+q32t**2*g22t
set resi3 = a3t/sqrt(fv3)
set residsq = resi3**2
cor(qstats,number=12,span=4) resi3
cor(qstats,number=12,span=4) residsq
*
* Compute and graph correlations
*
set rho21 = q21t*sqrt(g11t/fv2)
set rho31 = q31t*sqrt(g11t/fv3)
set rho32 = (q31t*q21t*g11t+q32t*g22t)/sqrt(fv2*fv3)
spgraph(vfields=3,footer="Figure 10.14 Time-varying correlation coefficients")
graph(header="(a) S&P 500 index and Cisco")
# rho21
graph(header="(b) S&P 500 index and Intel")
# rho31
graph(header="(a) Cisco and Intel")
# rho32
spgraph(done)