*
* @YuleVAR(options)  start end
* # list of dependent variables
*
* Estimates a VAR on stationary data using the Yule-Walker equations.
*
* Parameters:
*    start  end  range to use for estimation. Note that because it uses the
*                stationary covariance calculation, there are no observations
*                lost to lags.
*
* Options:
*   LAGS=(maximum) number of lags [1]
*   MODEL=(output) model computed by procedure
*   CRITERION=AIC/SBC/HQ/[NONE]
*     If you choose one of the criteria (AIC, SBC or HQ), the lag length is
*     chosen by minimizing the chosen criterion. By default, the lag length
*     is what you choose with the LAGS option.
*   PRINT/[NOPRINT]
*   DET=NONE/[CONSTANT]
*     Deterministic components, either none or the constant.
*
* Revision Schedule
*   05/2005 Written by Tom Doan, Estima.
*   10/2009 Added CRIT=HQ option
*
procedure YuleVAR start end
type integer start end
*
option integer lags 1
option model   *model
option choice  criterion 4 aic sbc hq none
option switch  print  0
option choice  det    2 none constant
*
local integer n i j k l ntotal
local integer startl endl
*
local symm cmom
local rect cmomx cmombase
local vect means
local rect pullcoeff modelcoeff
local vect[int] reorder
local real lrbase penalty cvalue
local series covars
local vect[int] list
local integer ndet
local real bestvalue
local integer bestlag
*
* We need N(L+1)+1 variable/lag combinations, for lags 0 to L of the
* endogenous variables, plus (optionally) the intercept
*
enter(varying) list
compute n=%rows(list)
compute ndet  =(det==2)
compute ntotal=n*(lags+1)+ndet
*
dim cmom(ntotal,ntotal) means(n)
*
* Figure out the desired range
*
inquire(reglist) startl<<start endl<<end
# list
*
* Generate the cross product matrix but with Cxy(s) used when the
* difference between the x and y lags is s.
*
do i=1,n
   stats(noprint) list(i) startl endl
   compute means(i)=%mean
   do j=1,n
      cross(noprint,cov,center=(det==2),corrs=covars) list(i) list(j) startl endl 0 lags
      do l=0,lags
         do k=0,lags-l
            compute cmom(k*n+i,(k+l)*n+j)=covars(l+1)
         end do k
      end do l
   end do j
end do i
*
* Copy over from symmetric to rectangular.
*
compute cmomx=cmom
*
* The covariances already have the means pulled out. Border the array
* with the effect of sweeping out the intercept (means in the bottom
* row, -means in the right column, 1 in the lower right corner).
*
if ndet==1 {
   do i=1,n
      do l=0,lags
         compute cmomx(ntotal,l*n+i)=means(i),cmomx(l*n+i,ntotal)=-means(i)
      end do l
   end do i
   compute cmomx(ntotal,ntotal)=1.0
}
*
if criterion==4 {
   compute cmomx=%sweeplist(cmomx,%seq(n+1,ntotal-ndet))
   compute bestlag=lags
}
else {
   compute cmombase=cmomx
   if criterion==1
      report(action=define,hlabels=||"Lags","AICC"||)
   else
   if criterion==2
      report(action=define,hlabels=||"Lags","SBC"||)
   else
      report(action=define,hlabels=||"Lags","HQ"||)

   compute lrbase=-2.0*%logconcdensity(%xsubmat(cmomx,1,n,1,n),%nobs)
   compute %nregsystem=n*ndet
   compute penalty=$
      %if(criterion==1,%nregsystem*2.0*%nobs*n/(%nobs*n-%nregsystem-1),$
      %if(criterion==2,%nregsystem*log(%nobs),$
                       %nregsystem*2.0*log(log(%nobs))))
   compute cvalue=lrbase+penalty
   compute bestvalue=cvalue,bestlag=0
   report(row=new,atcol=1) 0 cvalue
   do l=1,lags
      compute cmomx =%sweeplist(cmomx,%seq(n*l+1,n*(l+1)))
      compute lrbase=-2.0*%logconcdensity(%xsubmat(cmomx,1,n,1,n),%nobs)
      compute %nregsystem=n*ndet+l*n**2
      compute penalty=$
       %if(criterion==1,%nregsystem*2.0*%nobs*n/(%nobs*n-%nregsystem-1),$
       %if(criterion==2,%nregsystem*log(%nobs),$
                        %nregsystem*2.0*log(log(%nobs))))
      compute cvalue=lrbase+penalty
      if cvalue<bestvalue
         compute bestvalue=cvalue,bestlag=l
      report(row=new,atcol=1) l cvalue
   end do l
   *
   report(action=format,tag=minimum,atcol=2,align=decimal,special=onestar)
   report(action=show)
   if bestlag>0
      compute cmomx=%sweeplist(cmombase,%seq(n+1,n*(bestlag+1)))
   else
      compute cmomx=cmombase
}
*
* Set up a VAR model in standard order
*
if %defined(model) {
   system(model=model)
   variables list
   if bestlag>0
      lags 1 to bestlag
   if ndet
      det constant
   end(system)
   *
   * The coefficients are in the first n columns, below the top nxn
   * corner (which has the covariance matrix of residuals).
   *
   compute pullcoeff=%xsubmat(cmomx,n+1,ntotal,1,n)
   *
   * However, the coefficients are grouped by lag, rather than by
   * variable, and include possibly unnecessary lags
   *
   dim reorder(bestlag*n+ndet) modelcoeff(bestlag*n+ndet,n)
   if bestlag>0
      ewise reorder(i)=(%clock(i,bestlag)-1)*n+(i-1)/bestlag+1
   if ndet
      compute reorder(bestlag*n+ndet)=n*lags+1
   ewise modelcoeff(i,j)=pullcoeff(reorder(i),j)
   compute %modelsetcoeffs(model,modelcoeff)
   compute %modelsetvcv(model,%xsubmat(cmomx,1,n,1,n))
   *
   if print
      disp
      disp "Estimated VAR Equations"
      disp
      do i=1,n
         modify(print) %modeleqn(model,i)
      end do i
}
end

