*
* @MVGARCHtoVECH( options )
* Extracts a VECH representation out of the most recent GARCH
* instruction.
*
* Options:
*
* MV=[standard]/bekk/diagonal/vech
* ASYMMETRIC/[NOASYMMETRIC]
*   Repeat the options used on the GARCH instruction. These are the only
*   multivariate types which can be converted to VECH.
*
* The VECH representation has
*
* vech H(t) = C + B vech H(t-1) + A vech e(t-1)e(t-1)' + D vect v(t-1)v(t-1)'
*
* where v(t)=e(t) with positive elements zeroed out. The D is only
* present if you add the ASYMMETRIC option. vech takes a symmetric
* matrix and strings it out into a vector according to the rows of the
* lower triangle, that is, the components of a 3 variable model will be
* in the order
*
*  1
*  2 3
*  4 5 6
*
* Variables Defined:
*   %%VECH_C  C vector in the formula above
*   %%VECH_A  A matrix in the formula above
*   %%VECH_B  B matrix in the formula above
*   %%VECH_D  D matrix in the formula above
*
dec vector %%VECH_C
dec rect   %%VECH_A %%VECH_B

*  Revision Schedule:
*   04/2012 Written by Tom Doan, Estima.
*   08/2013 Correct MV=VECH,ASYMMETRIC case
*
procedure MVGARCHtoVECH
*
option choice mv 1 standard bekk diagonal vech
option switch asymmetric 0
*
local integer ncomp
local vect cx ax bx dx
local rect cb
local symm cs
local rect ab bb db
local integer i j k l m n
*
dec vector %%VECH_C
dec rect   %%VECH_A %%VECH_B
*
compute ncomp=%nvar*(%nvar+1)/2
dim %%VECH_C(ncomp) %%VECH_A(ncomp,ncomp) %%VECH_B(ncomp,ncomp)
compute %%VECH_D=%zeros(ncomp,ncomp)
*
if mv==1 {
   *
   * Standard/DVEC - convert to equivalent VECH representation
   *
   dim cx(ncomp) ax(ncomp) bx(ncomp) dx(ncomp)
   ewise cx(i)=%beta(%nregmean+i)
   ewise ax(i)=%beta(%nregmean+ncomp+i)
   ewise bx(i)=%beta(%nregmean+2*ncomp+i)
   compute %%VECH_A=%diag(ax)
   compute %%VECH_B=%diag(bx)
   compute %%VECH_C=cx
   if asymmetric {
      ewise dx(i)=%beta(%nregmean+3*ncomp+i)
      compute %%VECH_D=dx
	}
}
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, B and D coefficient matrices are stored in the
   * order A(1,1),A(1,2),..., which is actually the way A', B' and D' are
   * stored internally. Thus ab, bb and db below (which are just reshaped
   * subvectors) will be A', B' and D'.
   *
   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) db(%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)
   ewise %%VECH_A(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 %%VECH_B(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))
   compute %%VECH_C=CX
   if asymmetric {
      compute db=%vectorect(%xsubvec(%beta,%nregmean+ncomp+2*%nvar^2+1,$
                                        %nregmean+ncomp+3*%nvar^2),%nvar)
      ewise %%VECH_D(m,n)=i=%symmrow(m),j=%symmcol(m),k=%symmrow(n),l=%symmcol(n),$
            %if(k==l,db(i,k)*db(j,l),db(i,k)*db(j,l)+db(i,l)*db(j,k))
   }
}
else
if mv==4 {
   *
   * VECH - just copy and reshape matrices
   *
   compute %%VECH_C=%xsubvec(%beta,%nregmean+1,%nregmean+ncomp)
   compute %%VECH_A=tr(%vectorect(%xsubvec(%beta,%nregmean+ncomp+1         ,$
          %nregmean+ncomp+ncomp^2)  ,ncomp))
   compute %%VECH_B=tr(%vectorect(%xsubvec(%beta,%nregmean+ncomp+ncomp^2+1,$
          %nregmean+ncomp+2*ncomp^2),ncomp))
   if asymmetric {
      compute %%VECH_D=tr(%vectorect(%xsubvec(%beta,%nregmean+ncomp+2*ncomp^2+1,$
          %nregmean+ncomp+3*ncomp^2),ncomp))
   }
}
else
if mv==3 {
   dim cx(%nvar) ax(%nvar) bx(%nvar) dx(%nvar)
   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)
   compute %%VECH_C=%vec([symm] %diag(cx))
   compute %%VECH_A=%diag(%vec([symm] %diag(ax)))
   compute %%VECH_B=%diag(%vec([symm] %diag(bx)))
   if asymmetric {
      compute dx=%xsubvec(%beta,%nregmean+1+%nvar*3,%nregmean+%nvar*4)
      compute %%VECH_D=%diag(%vec([symm] %diag(dx)))
   }
}
*
end

