*
*  GARCHMV.PRG
*  Updated, March 2003 to include dynamic conditional correlation
*  Also, all models have been rewritten to stuff the uu' into a
*  matrix of series. This allows a bit more flexibility in handling
*  both initial conditions, and parameterizing the functions.
*
cal 1986 1 12
all 1996:12
*
open data returns.xls
data(format=xls,org=cols)
COMPUTE GSTART=1986:2 , GEND=1996:12
*
*  Parameters for the regression function
*
dec vect[series] y(2) u(2)
dec vect[frml] resid(2)
set y(1) = sp500
set y(2) = spmidcap
*
NONLIN(parmset=MEANPARMS) B11 B21
FRML RESID(1) = (Y(1)-B11)
FRML RESID(2) = (Y(2)-B21)
*
*  Do initial regression. Copy initial values for regression parameters
*
LINREG Y(1) / u(1)
# CONSTANT
COMPUTE B11 = %BETA(1)
LINREG Y(2) / u(2)
# CONSTANT
COMPUTE B21 = %BETA(1)
*
*  Get the covariance matrix of the residuals.
*
VCV(MATRIX=RR,NOPRINT)
# U
*
* h will have the sequence of variance estimates
* uu will have the sequence of uu' matrices
*
declare symm[series] h(2,2)
declare symm[series] uu(2,2)
*
* hx and uux are used when extracting elements from h and uu.
* ux is used when extracting a u vector
*
declare symm hx(2,2) uux(2,2)
declare vect ux(2)
*
* This is used to initialize pre-sample variances.
* If you want the pre-sample uu' to be the unconditional variance,
* change the right side of the set uu(i,j) to rr(i,j) (same as h).
*
do i=1,2
   do j=1,i
      set h(i,j)  = rr(i,j)
      set uu(i,j) = 0.0
   end do j
end do i
*
*
*  This is a standard log likelihood formula for any bivariate
*  ARCH, GARCH, ARCH-M,... The difference among these will be in
*  the definitions of HF and RESID. The function %XT pulls information
*  out of a matrix of SERIES, while %PT puts information into one.
*
declare frml[symm] hf
*
FRML LOGL = $
    U(1) = RESID(1) , U(2) = RESID(2) ,$
    HX = HF(T) , $
    UX = %XT(U,T) , UUX = %OUTERXX(UX),$
    %PT(H,T,HX),%PT(UU,T,%OUTERXX(UX)),$
    %LOGDENSITY(HX,UX)
*
*  Simple GARCH(1,1)
*
dec symm vc(2,2) va(2,2) vb(2,2)
nonlin(parmset=garchparms) vc va vb
frml hf = ||vc(1,1)+va(1,1)*h(1,1){1}+vb(1,1)*uu(1,1){1}|$
            vc(1,2)+va(1,2)*h(1,2){1}+vb(1,2)*uu(1,2){1},$
            vc(2,2)+va(2,2)*h(2,2){1}+vb(2,2)*uu(2,2){1}||
*
*  Initialize GARCH parameters
*
compute vc = rr , vb = %mscalar(0.05) , va = %mscalar(0.05)
*
*  Use simplex for a few iterations to get initial conditions
*  straightened out
*
NLPAR(SUBITS=50)
defaults maximize(trace)
MAXIMIZE(parmset=meanparms+garchparms,METHOD=SIMPLEX,ITERS=5) LOGL GSTART GEND
MAXIMIZE(parmset=meanparms+garchparms,METHOD=Bfgs,ITERS=100) LOGL GSTART GEND
*
*  Constant correlation
*
dec vect vbv(2) vav(2)
nonlin(parmset=garchparms) vc vbv vav
frml hf = (h11=vc(1,1)+vav(1)*h(1,1){1}+vbv(1)*uu(1,1){1}),$
          (h22=vc(2,2)+vav(2)*h(2,2){1}+vbv(2)*uu(2,2){1}),$
          ||h11|vc(1,2)*sqrt(h11*h22),h22||
compute vc = rr, vc(1,2)=vc(1,2)/sqrt(vc(1,1)*vc(2,2))
compute vbv=%const(0.05),vav=%const(0.05)
MAXIMIZE(parmset=meanparms+garchparms,METHOD=SIMPLEX,ITERS=5) LOGL GSTART GEND
MAXIMIZE(parmset=meanparms+garchparms,METHOD=Bfgs,ITERS=100) LOGL GSTART GEND
*
*  Full vech parameterization
*
dec vect vcv(3)
dec rect var(3,3) vbr(3,3)
nonlin(parmset=garchparms) vcv var vbr
*
dec frml[vect] hfv
frml hfv = $
  ([VECTOR] VECHH=%VEC(%XT(H,T-1))),$
  ([VECTOR] VECHU=%VEC(%XT(UU,T-1))),$
  VCV + VAR * VECHH + VBR * VECHU
frml hf = (vechh=hfv(t)),||vechh(1)|vechh(2),vechh(3)||
*
compute vcv = %vec(rr)
compute var = %mscalar(0.05), vbr = %mscalar(0.05)
*
MAXIMIZE(parmset=meanparms+garchparms,METHOD=SIMPLEX,ITERS=5) LOGL GSTART GEND
MAXIMIZE(parmset=meanparms+garchparms,METHOD=BFGS,ITERS=100) LOGL GSTART GEND
*
*  Positive definite parameterization (BEKK,EK)
*  This enforces a positive definite covariance matrix by writing the
*  covariance matrix evolution as
*
*   V(t) = C'C + B'u(t)u(t)'B + A'V(t-1)A
*
*  Note that the parameters are not globally identified: changing the signs
*  of all members of C,B or A will have no effect on the function value.
*  Using METHOD=SIMPLEX to begin is quite important with this setup, to
*  pull the estimates away from zero before starting the derivative-based
*  methods.
*
dec rect var(2,2) vbr(2,2)
dec rect vcr(2,2)
nonlin(parmset=garchparms) var vbr vcr vcr(1,2)=0.0

FRML HF = $
  (HX=%XT(H,T-1)),(UUX=%XT(UU,T-1)),$
  %INNERXX(VCR)+%MQFORM(HX,VAR)+%MQFORM(UUX,VBR)
*
*  Initialize c's from the decomp of the covariance matrix
*
COMPUTE vcr = %decomp(rr)
compute var = %mscalar(.05) , vbr = %mscalar(.05)
*
MAXIMIZE(parmset=meanparms+garchparms,METHOD=SIMPLEX,ITERS=5) LOGL GSTART GEND
MAXIMIZE(parmset=meanparms+garchparms,METHOD=Bfgs,ITERS=100) LOGL GSTART GEND
*
* Dynamic Conditional Constant Correlation
* Engle, JBES 2002, pp 339-350
*
* This uses a bivariate GARCH(1,1) model with fixed parameters to generate
* a sequence of p.d. matrices (Q's) which are used only for their implied correlation.
* The main model is otherwise structured like the constant correlation variety.
*
declare symm[series] q(2,2)
declare symm qx(2,2)
*
declare frml[symm] qf
*
* Initialize the q sequence
*
do i=1,2
   do j=1,i
      set q(i,j)  = rr(i,j)
   end do j
end do i

dec vect vbv(2) vav(2) vcv(2)
dec real a b
*
* a and b are the parameters governing the "GARCH" process of the Q sequence
*
nonlin(parmset=garchparms) vcv vbv vav a b
*
frml qf = (qx=(1-a-b)*rr+a*%xt(uu,t-1)+b*%xt(q,t-1)),%pt(q,t,qx),qx
frml hf = qf(t),rho=%if(a<1.and.b<1,qx(1,2)/sqrt(qx(1,1)*qx(2,2)),%na),$
          (h11=vcv(1)+vav(1)*h(1,1){1}+vbv(1)*uu(1,1){1}),$
          (h22=vcv(2)+vav(2)*h(2,2){1}+vbv(2)*uu(2,2){1}),$
          ||h11|rho*sqrt(h11*h22),h22||
*
* Initialize the c's to the diagonal elements of rr, others with typical values
*
compute vcv=%xdiag(rr),vbv=%const(0.05),vav=%const(0.05), a=0.05, b=0.05
MAXIMIZE(parmset=meanparms+garchparms,METHOD=SIMPLEX,ITERS=10) LOGL GSTART GEND
MAXIMIZE(parmset=meanparms+garchparms,METHOD=Bfgs,ITERS=100) LOGL GSTART GEND