*
* @MVGARCHVarMats(variances=VARIANCE model from GARCH)
*
* This pulls the matrices for a variance recursion out of the output
* from the GARCH model. Do this immediately after the GARCH instruction
* and repeat the VARIANCES option.
*
* Options:
*   VARIANCES=[SIMPLE]/VARMA/EXPONENTIAL/SPILLOVER/KOUTMOS
*    (note that EXPONENTIAL and KOUTMOS are provided merely to give an
*     informative error message. They don't admit a simple recursion formula
*     for the variance).
*
* For GARCH(1,1) models with CC or DCC, and VARIANCES=SIMPLE, VARMA or
* SPILLOVER, the variances themselves (not the covariances) can be
* extrapolated using a matrix recursion:
*
*   var(t) = C + A eps(t-1).^2 + B var(t-1)
*
* where var(t) is the n-vector of variances, C is an n-vector, A and B
* are n x n matrices, and eps(t-1).^2 is the vector of squares of
* residuals (.^2 meaning the elementwise squares). With
* VARIANCES=SIMPLE, A and B are diagonal, for VARIANCES=SPILLOVER, A is
* a full n x n matrix and B is diagonal, and for VARIANCES=VARMA, both A
* and B are full matrices.
*
* Variables Defined:
*   %%VAR_C   The C vector above
*   %%VAR_A   The A matrix above
*   %%VAR_B   The B matrix above
*
procedure MVGARCHVarMats
*
option choice variances 1 simple varma exponential spillover koutmos
option switch asymmetric 0
*
dec vect %%var_c
dec rect %%var_a %%var_b %%var_d
local integer cpos apos bpos dpos
*
if variances==3.or.variances==5 {
   disp "@MVGARCHVarMats doesn't apply to exponential or Koutmos variances"
   return
}
if variances==1 {
   compute cpos=%nregmean+1
   compute apos=cpos+%nvar
   compute bpos=apos+%nvar
   compute dpos=bpos+%nvar
   compute %%var_c=%xsubvec(%beta,cpos,cpos+%nvar-1)
   compute %%var_a=%diag(%xsubvec(%beta,apos,apos+%nvar-1))
   compute %%var_b=%diag(%xsubvec(%beta,bpos,bpos+%nvar-1))
   compute %%var_d=%if(asymmetric,%diag(%xsubvec(%beta,dpos,dpos+%nvar-1)),%zeros(%nvar,%nvar))
}
else
if variances==2 {
   compute cpos=%nregmean+1
   compute apos=cpos+%nvar
   compute bpos=apos+%nvar^2
   compute dpos=bpos+%nvar^2
   compute %%var_c=%xsubvec(%beta,cpos,cpos+%nvar-1)
   compute %%var_a=tr(%vectorect(%xsubvec(%beta,apos,apos+%nvar*%nvar-1),%nvar))
   compute %%var_b=tr(%vectorect(%xsubvec(%beta,bpos,bpos+%nvar*%nvar-1),%nvar))
   compute %%var_d=%if(asymmetric,%diag(%xsubvec(%beta,dpos,dpos+%nvar-1)),%zeros(%nvar,%nvar))
}
else
if variances==4 {
   compute cpos=%nregmean+1
   compute apos=cpos+%nvar
   compute bpos=apos+%nvar^2
   compute dpos=bpos+%nvar
   compute %%var_c=%xsubvec(%beta,cpos,cpos+%nvar-1)
   compute %%var_a=tr(%vectorect(%xsubvec(%beta,apos,apos+%nvar*%nvar-1),%nvar))
   compute %%var_b=%diag(%xsubvec(%beta,bpos,bpos+%nvar-1))
   compute %%var_d=%if(asymmetric,%diag(%xsubvec(%beta,dpos,dpos+%nvar-1)),%zeros(%nvar,%nvar))
}
end
