*
* MVGARCHFore computes out-of-sample variance forecasts for many of the
* standard multivariate GARCH models. This applies to the model most
* recently estimated with GARCH.
*
* Syntax:
*
* @MVGARCHFore( options)  h u
*
* Parameters:
*   h   (input and output) is the SERIES[SYMM] which has the in-sample
*       estimates of the covariances (produced with the HMATRICES option
*       on GARCH), and will be extended to hold the forecasts of the
*       covariances
*   u   (input) is the SERIES[VECT] of the residuals. This is produced by
*       the RVECTORS option on GARCH.
*
* Options:
*
*  MV=[standard]/bekk/diagonal/cc/dcc/vech/dvech
*    This should be the same choice that you made on the GARCH instruction.
*    Note that DCC has no simple extrapolation formula. You also can't
*    apply this if you used the VARIANCES=EXP or VARIANCES=VARMA with CC or
*    DIAG.
*  STEPS=number of out-of-sample forecast steps
*
*  ACCUMULATE/[NOACCUMULATE]
*    If ACCUMULATE, returns the partial sums of the forecasts (which would be
*    the forecast of the variance of multiple step returns if the raw data
*    are period-to-period returns).
*
*  Revision Schedule:
*   02/2006 Written by Tom Doan, Estima.
*   05/2006 Revised, adding comments and correcting error for BEKK branch
*   08/2010 Added check for required inputs
*   11/2025 Added ACCUMULATE option
*
* Utility function which converts a correlation matrix whose subdiagonal
* values are packed into an N-1 x N-1 "symmetric" (the coding used by
* the GARCH instruction for the free parameters in the CC model) and a
* vector of variances into a full covariance matrix.
*
function %corrsdtocv corrsd var
type symm %corrsdtocv corrsd var
type vector var
*
local integer i j n
*
compute n = %rows(corrsd)
dim %corrsdtocv(n+1,n+1)
ewise %corrsdtocv(i,j)=%if(i==j,var(i),corrsd(i-1,j)*sqrt(var(i)*var(j)))
end
*
procedure MVGARCHFore h u
type series[symm] *h
type series[vect] u
*
option choice mv 1 standard bekk diagonal cc dcc vech dvech
option integer steps 12
option switch accum 0
*
local integer ncomp
local vect cx ax bx
local series[symm] uu
local rect af bf
local rect cb
local symm cs
local symm rr
local rect ab bb
local integer i j k l m n
*
gset uu %regstart() %regend() = %outerxx(u)
compute ncomp=%nvar*(%nvar+1)/2
if .not.%defined(h).or..not.%defined(u) {
   disp "Usage: @MVGARCHFore(options)  H_Matrices  U_Vectors"
   return
}
if mv==5 {
   disp "No Simple Extrapolative Formula Exists for DCC Model"
   return
}
else
if mv==1.or.mv==7 {
   *
   * Standard/DVEC - convert to equivalent VECH representation
   *
   dim cx(ncomp) ax(ncomp) bx(ncomp)
   ewise cx(i)=%beta(%nregmean+i)
   ewise ax(i)=%beta(%nregmean+ncomp+i)
   ewise bx(i)=%beta(%nregmean+2*ncomp+i)
   compute af=%diag(ax)
   compute bf=%diag(bx)
}
else
if mv==2 {
   *
   * BEKK - convert to equivalent VECH representation. This is a bit
   * clumsy because it requires packing an NxN Kroneker product into the
   * shape needed for calculations on packed symmetrics. Because only
   * the (i,j) element shows in the packed matrix (i>j), the coefficient
   * on (j,i) needs to be combined into that.
   *
   * Note that the A and B coefficient matrices are stored in the order
   * A(1,1),A(1,2),..., which is actually the way A' and B' are stored
   * internally. Thus ab and bb below (which are just reshaped
   * subvectors) will be A' and B'.
   *
   dim cx(ncomp)
   dim cb(%nvar,%nvar) cs(%nvar,%nvar)
   compute cs=%vectosymm(%xsubvec(%beta,%nregmean+1,%nregmean+ncomp),%nvar)
   ewise cb(i,j)=%if(j<=i,cs(i,j),0.0)
   compute cs=cb*tr(cb)
   compute cx=%vec(cs)
   *
   dim ab(%nvar,%nvar) bb(%nvar,%nvar)
   compute ab=%vectorect(%xsubvec(%beta,%nregmean+ncomp+1         ,$
                                        %nregmean+ncomp+  %nvar^2),%nvar)
   compute bb=%vectorect(%xsubvec(%beta,%nregmean+ncomp+%nvar^2+1,$
                                        %nregmean+ncomp+2*%nvar^2),%nvar)
   dim af(ncomp,ncomp) bf(ncomp,ncomp)
   ewise af(m,n)=i=%symmrow(m),j=%symmcol(m),k=%symmrow(n),l=%symmcol(n),$
            %if(k==l,ab(i,k)*ab(j,l),ab(i,k)*ab(j,l)+ab(i,l)*ab(j,k))
   ewise bf(m,n)=i=%symmrow(m),j=%symmcol(m),k=%symmrow(n),l=%symmcol(n),$
            %if(k==l,bb(i,k)*bb(j,l),bb(i,k)*bb(j,l)+bb(i,l)*bb(j,k))
}
else
if mv==6 {
   *
   * VECH - just copy and reshape matrices
   *
   compute cx=%xsubvec(%beta,%nregmean+1,%nregmean+ncomp)
   dim af(ncomp,ncomp) bf(ncomp,ncomp)
   compute af=tr(%vectorect(%xsubvec(%beta,%nregmean+ncomp+1         ,$
          %nregmean+ncomp+ncomp^2)  ,ncomp))
   compute bf=tr(%vectorect(%xsubvec(%beta,%nregmean+ncomp+ncomp^2+1,$
          %nregmean+ncomp+2*ncomp^2),ncomp))
}
else
if mv==3.or.mv==4 {
   compute cx=%xsubvec(%beta,%nregmean+1        ,%nregmean+%nvar)
   compute ax=%xsubvec(%beta,%nregmean+1+%nvar  ,%nregmean+%nvar*2)
   compute bx=%xsubvec(%beta,%nregmean+1+%nvar*2,%nregmean+%nvar*3)
   if mv==3
      compute rr=%zeros(%nvar-1,%nvar-1)
   else
      compute rr=%vectosymm(%xsubvec(%beta,%nregmean+1+%nvar*3,$
                         %nregmean+%nvar*3+%nvar*(%nvar-1)/2),%nvar-1)
}
gset uu %regend()+1 %regend()+steps = %zeros(%nvar,%nvar)
if mv==1.or.mv==2.or.mv==6
   gset h  %regend()+1 %regend()+steps = $
              uu=%vectosymm(cx+bf*%vec(h{1})+af*%vec(uu{1}),%nvar)
else
   gset h  %regend()+1 %regend()+steps = $
              uu=%corrsdtocv(rr,cx+bx.*%xdiag(h{1})+ax.*%xdiag(uu{1}))

if accum
   gset h %regend()+1 %regend()+steps = %if(t==%regend()+1,h,h{1}+h)
end
