*
* Example 10.5 from pp 462-464
*
open data m-ibmspln.dat
calendar(m) 1926
data(format=free,org=columns) 1926:1 1999:12 ibm sp500
*
equation ibmeq ibm
# constant ibm{1 2} sp500{2}
equation speq sp500
# constant
group bimean ibmeq speq
*
* The volatility model in (10.25) is almost what you would get by using the
* combination of MV=CC and VARIANCES=VARMA. The VARIANCES option specifies the
* formula used for computing the variances; VARMA is, in effect, a vector ARMA
* specification for those. The only difference between this and (9.22) is the 0
* restriction on the 1,2 element on the lagged squared residuals matrix.
*
garch(mv=cc,variances=varma,p=1,q=1,model=bimean,hmatrices=hbi,rvector=rbi,$
  pmethod=simplex,piters=10,iters=200)
set ibmstd = rbi(t)(1)/sqrt(hbi(t)(1,1))
set spstd  = rbi(t)(2)/sqrt(hbi(t)(2,2))
@MVQStat(lags=4,dfc=3)
# ibmstd spstd
@MVQStat(lags=8,dfc=3)
# ibmstd spstd
*
set ibmstdsq = ibmstd**2
set spstdsq  = spstd**2
@MVQStat(lags=4,dfc=4)
# ibmstdsq spstdsq
@MVQStat(lags=8,dfc=4)
# ibmstdsq spstdsq
*
* Continued (page 466)
*
* Compute the moving window correlations, creating the series "movingcorr"
*
do time=1935:12,1999:12
   cmom(corr) time-199 time
   # ibm sp500
   set movingcorr time time = %cmom(1,2)
end do time
graph(footer="Figure 10.9 Sample Correlation with Moving Window")
# movingcorr
*
* Time varying correlation model
*
compute n=2
dec vect[series] u(n)
dec vect[frml] resid(n)
dec vect[frml] hd(n)
dec frml[symm] hf
dec series[symm] uu h
dec symm hx(n,n) uux(n,n)
dec vect ux(n)
*
nonlin(parmset=meanparms) p10 p11 p12 p20
nonlin(parmset=garchparms) c1 c2 a11 a21 a22 b11 b12 b21 b22 q0 q1 q2
*
frml resid(1) = ibm-p10-p11*ibm{1}-p12*sp500{2}
frml resid(2) = sp500-p20
nlsystem(parmset=meanparms,resids=u) / resid
gset uu   = %sigma
gset h    = %sigma
set rhotv = .60
*
frml hd(1) = c1+a11*uu{1}(1,1)+b11*h{1}(1,1)+b12*h{1}(2,2)
frml hd(2) = c2+a21*uu{1}(1,1)+a22*uu{1}(2,2)+b21*h{1}(1,1)+b22*h{1}(2,2)
frml qf    = q0+q1*rhotv{1}+q2*u(1){1}*u(2){1}/sqrt(h{1}(1,1)*h{1}(2,2))
frml hf    = hx(1,1)=hd(1),hx(2,2)=hd(2),rhotv=%logistic(qf,1.0),hx(1,2)=rhotv*sqrt(hx(1,1)*hx(2,2)),hx
frml logl = $
    hx = hf(t) , $
    %do(i,1,n,u(i)=resid(i)) , $
    ux = %xt(u,t), $
    h(t)=hx, uu(t)=%outerxx(ux), $
    %logdensity(hx,ux)
compute c1  = 0.1 , a11 = 0.1 , b11 = 0.4 , b12 = 0.0
compute c2  = 0.1 , a21 = 0.1 , a22 = 0.1 , b21 = 0.2 , b22 = 0.5
compute q0  = 0.5 , q1=0.0, q2=0.0
*
maximize(parmset=meanparms+garchparms,pmethod=simplex,piters=20,method=bfgs,iters=200) logl 7 *
set ibmstd = resid(1)/sqrt(h(t)(1,1))
set spstd  = resid(2)/sqrt(h(t)(2,2))
*
@MVQStat(lags=4,dfc=2)
# ibmstd spstd
@MVQStat(lags=8,dfc=2)
# ibmstd spstd
*
set ibmstdsq = ibmstd**2
set spstdsq  = spstd**2
@MVQStat(lags=4,dfc=7)
# ibmstdsq spstdsq
@MVQStat(lags=8,dfc=7)
# ibmstdsq spstdsq
*
graph(footer="Figure 10.10 Fitted Conditional Correlation Coefficient",vgrid=.612)
# rhotv