*
* Example 9.1 from pp 365-368
* ?? Unconditional variance calcs on MVGARCH
* ?? Spec for 2nd model (clearly wrong in part)
*
open data hkja.dat
data(format=free,org=columns) 1 469 hk ja
spgraph(vfields=2,header='Figure 9.1 Time Plots of Daily Log Returns')
graph(header='(a) Hong Kong')
# hk
graph(header='(b) Japan')
# ja
spgraph(done)
*
* Univariate GARCH models
*
garch(p=1,q=1,reg,resids=hkres,hseries=hkvar) / hk
# hk{1}
*
* Diagnostics on the standardized residuals
*
set ustd   = hkres/sqrt(hkvar)
set ustdsq = ustd**2
corr(qstats,dfc=1,number=4) ustd
corr(qstats,dfc=2,number=4) ustdsq
*
disp 'Unconditional Variance for Hong Kong' %beta(2)/(1-%beta(3)-%beta(4))
*
garch(p=1,q=1,nomean,resids=jares,hseries=javar) / ja
set ustd   = jares/sqrt(javar)
set ustdsq = ustd**2
corr(qstats,dfc=1,number=4) ustd
corr(qstats,dfc=2,number=4) ustdsq
disp 'Unconditional Variance for Japan' %beta(1)/(1-%beta(2)-%beta(3))
*
spgraph(vfields=2,header='Estimated Volatilities')
graph(header='(a) Hong Kong')
# hkvar
graph(header='(b) Japan')
# javar
spgraph(done)
*
* Bivariate GARCH models
*
* Do the mean equations. We need to define a model (using GROUP) which
* includes the two mean equations. The equation instruction for Japan
* uses the "empty" option, because the equation has no explanatory variables.
*
equation hkmean hk
# hk{6}
equation(empty) jamean ja
group bimean hkmean jamean
*
* Constant Correlation model estimates
*
garch(p=1,q=1,iters=200,mv=cc,model=bimean,method=bhhh,hmatrices=hbi,rvector=rbi)
compute gbeta=%beta
*
* Generate the standardized residuals and do diagnostics on them
*
set hkstd = rbi(t)(1)/sqrt(hbi(t)(1,1))
set jpstd = rbi(t)(2)/sqrt(hbi(t)(2,2))
@MVQStat(lags=4,dfc=1)
# hkstd jpstd
@MVQStat(lags=8,dfc=1)
# hkstd jpstd
*
set hkstdsq = hkstd**2
set jpstdsq = jpstd**2
@MVQStat(lags=4,dfc=4)
# hkstdsq jpstdsq
@MVQStat(lags=8,dfc=4)
# hkstdsq jpstdsq
*
* Unconditional variances
*
compute ucvarhk=gbeta(2)/(1-gbeta(4)-gbeta(6))
compute ucvarjp=gbeta(3)/(1-gbeta(5)-gbeta(7))
disp 'Unconditional Variance for Hong Kong' ucvarhk
disp 'Unconditional Variance for Japan' ucvarjp
*
* Restricted MVGARCH. Because this is rather non-standard, it's
* necessary to use MAXIMIZE rather than GARCH.
*
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) p6
nonlin(parmset=garchparms) c1 c2 a11 a22 b11 b21 b22 rho
*
* Do a separate estimate of the mean model by non-linear
* systems estimation. (It's actually a linear system, but
* it's easier to use nlsystem since everything is already set
* up as FRMLs). This gives us an initial estimate for the
* p6 parameter and also the fixed variance estimate of the
* covariance matrix, which is used in initializing the pre-sample
* variance.
*
frml resid(1) = hk-p6*hk{6}
frml resid(2) = ja
nlsystem(parmset=meanparms,resids=u) / resid
gset uu = %sigma
gset h  = %sigma
*
frml hd(1) = c1+a11*uu{1}(1,1)+b11*h{1}(1,1)
frml hd(2) = c2+a22*uu{1}(2,2)+b21*h{1}(1,1)+b22*h{1}(2,2)
frml hf = hx(1,1)=hd(1),hx(2,2)=hd(2),hx(1,2)=rho*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
compute a22 = 0.1 , b21 = 0.2 , b22 = 0.5
compute rho = 0.1
*
* The estimates here end up being quite different from those shown in the book.
* GARCH models with fairly large coefficients on the lagged variance terms tend
* to be quite sensitive to the choice from the pre-sample variances. The estimates
* here use the fixed matrix estimate from the nlsystem. Tsay's estimates use
* a zero matrix. His results can be (almost exactly) duplicated by switching the
* gset h above to = %zeros(2,2) and changing the starting point for estimation on
* the maximize to 8. (The only other difference is with the pre-sample values for
* uu, which has relatively little effect on the estimates).
*
maximize(parmset=meanparms+garchparms,pmethod=simplex,piters=20,method=bfgs,iters=200) logl 7 *
*
set hkstd = resid(1)/sqrt(h(t)(1,1))
set jpstd = resid(2)/sqrt(h(t)(2,2))
*
@MVQStat(lags=4,dfc=1)
# hkstd jpstd
@MVQStat(lags=8,dfc=1)
# hkstd jpstd
set hkstdsq = hkstd**2
set jpstdsq = jpstd**2
@MVQStat(lags=4,dfc=5)
# hkstdsq jpstdsq
@MVQStat(lags=8,dfc=5)
# hkstdsq jpstdsq
*
* Unconditional variances
*
compute lrmat=||1.0-a11-b11,0.0|0.0-b21,1.0-a22-b22||
compute ucvar=inv(lrmat)*||c1|c2||