*
* Example 10.7 from pp 478-480
*
open data m-ibmspln.dat
calendar(m) 1926
data(format=free,org=columns) 1926:1 1999:12 ibm sp500
*
* Estimate the simplified AR(3) model
*
system(model=simplified)
variables ibm sp500
det constant sp500{1 3}
end(system)
*
estimate(sigma,resids=a)
*
* Do the principal components on the %sigma matrix produced by ESTIMATE, Note
* that due to a slight difference in the scaling of the covariance matrix, the
* eigenvalues themselves are different. However, the proportions and eigenvectors
* are the same.
*
@prinfactors(print,vectors=ev) %sigma
*
* Principal components of original (centered) data
*
vcv(center)
# ibm sp500
@prinfactors(print,vectors=ev) %sigma
set x = ev(1,1)*ibm+ev(2,1)*sp500
*
garch(reg,p=1,q=1,hseries=hx,resids=ax) / x
# constant x{1}
*
set ustd = ax/sqrt(hx)
set ustdsq = ustd**2
@regcorrs(number=10) ustd
@regcorrs(number=10) ustdsq
*
spgraph(vfields=2,footer="Figure 10.15 GARCH Model applied to principal components")
graph(header="(a) The first principal component")
# x
graph(header="(b) A fitted volatility process")
# hx
spgraph(done)
*
nonlin(parmset=meanparms) b10 b11 b12 b13 b20
frml a1f = ibm-b10-b11*ibm{1}-b12*ibm{2}-b13*sp500{2}
frml a2f = sp500-b20
*
nlsystem(parmset=meanparms) / a1f a2f
*
set a1t = a1f
set a2t = a2f
*
set g11t = %sigma(1,1)
set g22t = %sigma(2,2)
set rho  = %sigma(1,2)/sqrt(%sigma(1,1)*%sigma(2,2))
*
nonlin(parmset=garchparms) c10 c11 c12 c20 c22 q0 q1 q2
compute c10=%sigma(1,1),c11=0.0,c12=0.0
compute c20=%sigma(2,2),c22=0.0
compute q0=0.0,q1=1.0,q2=0.0
*
frml g11f = c10+c11*a1t{1}**2+c12*hx
frml g22f = c20+c22*hx
frml qf   = q0+q1*rho{1}+q2*a1t{1}*a2t{1}/sqrt(g11t{1}*g22t{1})
frml logl = a1t=a1f,a2t=a2f,g11t=g11f,g22t=g22f,rho=%logistic(qf,1.0),$
   %logdensity(||g11t|rho*sqrt(g11t*g22t),g22t||,||a1t,a2t||)
maximize(parmset=meanparms+garchparms) logl 4 *
*
set a1std = a1t/sqrt(g11t)
set a2std = a2t/sqrt(g22t)
*
@mvqstat(lags=8)
# a1std a2std
set a1stdsq = a1std**2
set a2stdsq = a2std**2
@mvqstat(lags=8)
# a1stdsq a2stdsq