*
* Example of a multivariate GARCH-M model using square roots of the
* variances in the mean equation. The mean model is a VAR(1)+square
* roots of variances.
*
all 6237
open data g10xrate.xls
data(format=xls,org=columns) / usxjpn usxfra usxsui
*
set xjpn = 100.0*log(usxjpn/usxjpn{1})
set xfra = 100.0*log(usxfra/usxfra{1})
set xsui = 100.0*log(usxsui/usxsui{1})
*
compute n=3
compute gstart=3,gend=6237
*
dec series[vect] yv
dec frml[vect] residv
dec vect[series] u(n)
*
* The paths of the covariance matrices and uu' are saved in the
* SERIES[SYMM] named H and UU. UX and HX are used for the current values
* of the residual vector and H matrices
*
declare series[symm] h uu
*
* ux is used when extracting a u vector
*
declare symm hx(n,n)
declare vect ux(n)
*
* These will be the parameters for the mean equations. These are
* adjusted to add variance or covariance terms as needed.
*
dec vect b(n)
dec rect ar(n,n) gm(n,n)
nonlin(parmset=meanparms) b ar
*
* Mean model = VAR(1) with sqrt(h) "M" term
*
frml residv = yv-b-ar*yv{1}-gm*%sqrt(%xdiag(h))
*
gset h  = %zeros(n,n)
gset yv = ||xjpn,xfra,xsui||
*
* Get guess values for the mean parameters ignoring the "M" term by
* running an NLSYSTEM. (The system is actually linear, so this converges
* right away).
*
compute gm=%const(0.0)
nlsystem(parmset=meanparms,frml=residv) gstart gend
compute rr=%sigma
*
* Reset the parameter set for the means to include gm
*
nonlin(parmset=meanparms) b ar gm
*************************************************************************
*************************************************************************
*
* From this point, everything is the same as it would be for any GARCH
* model done using MAXIMIZE.
*
* These are used to initialize pre-sample variances.
*
gset h  * gend = rr
gset uu * gend = rr
*
declare frml[symm] hf
*
* Common formula for the log likelihood
*
frml logl = $
    hx    = hf(t) , $
    h(t)  = hx, $
    ux    = residv , $
    uu(t) = %outerxx(ux), $
    %pt(u,t,ux),$
    %logdensity(hx,ux)
*************************************************************
*
* Simple DVEC GARCH(1,1) model for the variance.
*
dec symm vcs(n,n) vas(n,n) vbs(n,n)
compute vcs=rr,vbs=%const(0.50),vas=%const(0.05)
nonlin(parmset=garchparms) vcs vas vbs
frml hf = vcs+vbs.*h{1}+vas.*uu{1}
maximize(parmset=meanparms+garchparms,$
   pmethod=simplex,piters=5,method=bfgs,iters=400) logl gstart gend
**************************************************************
*
* BEKK model for the variance
*
dec rect vab(n,n) vbb(n,n)
compute vab=.05*%identity(n),vbb=.50*%identity(n)
nonlin(parmset=garchparms) vcs vab vbb
compute vcs=%decomp(rr)
frml hf = vcs*tr(vcs)+vbb*h{1}*tr(vbb)+vab*uu{1}*tr(vab)
maximize(parmset=meanparms+garchparms,$
  pmethod=simplex,piters=5,method=bfgs,iters=400) logl gstart gend
**************************************************************
*
* CCC
* The correlations are parameterized using an (n-1)x(n-1) matrix for
* the subdiagonal. The (i,j) element of this will actually be the
* correlation between i+1 and j.
*
dec symm qc(n-1,n-1)
dec vect vcv(n) vbv(n) vav(n)
*
function hfcccgarch time
type symm hfcccgarch
type integer time
do i=1,n
   compute hx(i,i)=vcv(i)+vav(i)*h(time-1)(i,i)+vbv(i)*uu(time-1)(i,i)
   do j=1,i-1
     compute hx(i,j)=qc(i-1,j)*sqrt(hx(j,j)*hx(i,i))
   end do j
end do i
compute hfcccgarch=hx
end
*
frml hf = hfcccgarch(t)
nonlin(parmset=garchparms) vcv vbv vav qc
compute vcv=%xdiag(rr),vbv=%const(0.05),vav=%const(0.05),qc=%const(0.0)
maximize(parmset=meanparms+garchparms,$
  pmethod=simplex,piters=5,method=bfgs,iters=500) logl gstart gend

