*
* @VARFromDLM(options)
* Generates a VAR representation for a selection of variables from a
* stationary DLM.
*
* The underlying state space model takes the form
*
*  X(t)=AX(t-1)+Fw(t)
*
* The observables take the form
*
*  Y(t)=C'X(t)
*
* The procedure generates the implied VAR coefficients and covariance
* matrix of residuals given this structure.
*
* Options:
*   A =A matrix
*   F =F matrix [identity]
*   SW=Covariance matrix of W
*   C =loading matrix [identity]
*
*   LAGS=number of lags in the generated VAR [1]
*
*   CVOUT=Implied VAR covariance matrix
*   PHI=RECT[RECT] of implied VAR coefficients
*    This is a LAGS x LAGS array of m x m matrices. PHI(k,l) is an mxm
*    matrix of lag coefficients at lag <<l>> for a <<k>> lag VAR. Thus, the
*    final VAR will have coefficients in PHI(LAGS,1),...,PHI(LAGS,LAGS).
*
* Revision Schedule:
*  01/2010 Written by Tom Doan, Estima
*
procedure VARFromDLM
*
option rect    a
option symm    sw
option rect    f
option rect    c
option integer lags   1
option rect[rect] *phi
option symm       *cvout
*
local integer    nx ny j k
local rect[rect] phil
local vect[symm] vlfin
local vect[rect] cxfin cyfin
local rect       cx sfin
local symm       swx s0fin
*
if .not.%defined(a) {
   disp "#### Usage: @VARFromDLM(A=transition matrix,other options)"
   return
}
compute nx=%rows(a)
*
if %defined(f)
   compute swx=%mqform(sw,tr(f))
else
   compute swx=sw
*
if %defined(c)
   compute cx=c
else
   compute cx=%identity(nx)
compute ny=%cols(cx)
*
* Compute the ergodic covariance matrix
*
compute s0fin=%psdinit(a,swx)
*
* Build up full system covariogram. Because 1-based indices are used,
* cxx(j) is in entry j+1.
*
dim cxfin(lags+1) cyfin(lags+1)
compute cxfin(1)=s0fin
*
do i=2,lags+1
   compute cxfin(i)=a*cxfin(i-1)
end do i
*
* Generate the covariogram for the extracted variables.
*
do i=1,lags+1
   compute cyfin(i)=tr(cx)*cxfin(i)*cx
end do i
*
* Apply the Durbin-Levinson recursion to generate the VAR coefficients.
*
dim phil(lags,lags)
dim vlfin(lags+1)
*
compute vlfin(1)=cyfin(1)
*
do k=1,lags
   compute sfin=cyfin(k+1)
   do j=1,k-1
      compute sfin=sfin-phil(k-1,j)*cyfin(k+1-j)
   end do j
   compute phil(k,k)=inv(vlfin(k))*sfin
   do j=1,k-1
      compute phil(k,j)=phil(k-1,j)-phil(k,k)*phil(k-1,k-j)
   end do j
   do j=k+1,lags
      compute phil(k,j)=%zeros(ny,ny)
   end do j
   compute vlfin(k+1)=vlfin(k)-phil(k,k)*vlfin(k)*tr(phil(k,k))
end do k
*
if %defined(phi)
   compute phi=phil
if %defined(cvout)
   compute cvout=vlfin(k+1)
end