*
* @VARLagSelect( options )  start  end
* # list of variables
*
* VARLagSelect chooses the lag length which minimizes one of the
* information criteria.
*
* Options:
*   LAGS=(maximum) number of lags [1]
*   DET=NONE/[CONSTANT] deterministic variables to include
*   CRIT=AIC/SBC/BIC/HQ/GTOS/[NONE] Criterion to use in selecting the lag
*     length (SBC and BIC are synonyms). TEST does a general to specific
*     pruning using LR tests. If NONE, it just estimates the model with the
*     specified number of lags.
*   MODEL=model to define from chosen lag length [not used]
*   [PRINT]/NOPRINT
*
* Variables Defined:
*   %%AUTOP = chosen number of lags
*
* Revision Schedule:
*   06/2005 Written by Tom Doan, Estima.
*   03/2009 Add PRINT/NOPRINT option and define %%AUTOP
*   10/2009 Add CRIT=GTOS option
*
procedure VARLagSelect start end
type integer start end
*
option integer lags 1
option model   *model
option choice  criterion 6 aic sbc bic hq gtos none
option choice  det       2 none constant
option switch  print     1
*
declare integer %%autop
*
local integer n i j l ntotal ndet
local vect[int] list
*
local rect      cmomx cmombase
local rect      pullcoeff modelcoeff
local vect[int] reorder
local real      lrbase penalty cvalue lastlogl lrcrit
local vect[int] reglist
local integer   bestlag
local real      bestvalue
local string    critlab
*
* We need N(L+1)+1 variable/lag combinations, for lags 0 to L of the
* endogenous variables, plus the intercept
*
enter(varying) list
compute n=%rows(list)
compute ntotal=n*(lags+1)
*
* Generate list of required variables. Start with the dependent variables
*
dim reglist(0)
do i=1,n
   compute reglist=%rladdone(reglist,list(i))
end do i
*
* Now add any required lags. These are grouped by lags to make it easier
* to use the sweep operator to do the regressions.
*
do j=1,lags
   do i=1,n
      compute reglist=%rladdlag(reglist,list(i),j)
   end do i
end do j
*
* Add in the deterministic components
*
if det==2 {
   compute reglist=%rladdone(reglist,constant)
   compute ndet=1
}
else
   compute ndet=0
compute ntotal=ntotal+ndet
*
* Compute the cross product matrix
*
cmom start end
# reglist
*
* Scale it to a covariance matrix
*
compute cmomx=%cmom/%nobs
*
if criterion==6 {
   compute cmomx=%sweeplist(cmomx,%seq(n+1,ntotal))
   compute bestlag=lags
}
if criterion<>6 {
   if criterion==1
      compute critlab="AICC"
   else
   if criterion==2.or.criterion==3
      compute critlab="SBC/BIC"
   else
   if criterion==4
      compute critlab="HQ"
   else
      compute critlab="F-Test"

	if print
      report(action=define,hlabels=||"Lags",critlab||)
   *
   * Sweep out deterministic variables. Save the resulting matrix so we
   * can easily go back to the optimum number of lags.
   *
   if ndet>0
      compute cmomx=%sweep(cmomx,ntotal)
   compute cmombase=cmomx
   compute lrbase  =-2.0*%logconcdensity(%xsubmat(cmomx,1,n,1,n),%nobs)
   compute lastlogl=lrbase
   compute %nregsystem=n*ndet
   compute penalty=%if(criterion==1,%nregsystem*(2.0*n*%nobs/(n*%nobs-%nregsystem-1)),$
                   %if(criterion<>4,%nregsystem*log(%nobs),$
                                    %nregsystem*(2.0*log(log(%nobs)))))
   compute cvalue=lrbase+penalty
   if criterion==5
      compute cvalue=0.0,lrcrit=%invchisqr(.05,n^2)
   compute bestlag=0,bestvalue=cvalue
   if print
      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=l*n**2+n*ndet
      compute penalty=%if(criterion==1,%nregsystem*(2.0*n*%nobs/(n*%nobs-%nregsystem-1)),$
                      %if(criterion<>4,%nregsystem*log(%nobs),$
                                       %nregsystem*(2.0*log(log(%nobs)))))
      compute cvalue=lrbase+penalty
      if criterion==5 {
         compute cvalue=lastlogl-lrbase
         if cvalue>lrcrit
            compute bestlag=l
      }
      else
      if cvalue<bestvalue
         compute bestlag=l,bestvalue=cvalue
      if print
         report(row=new,atcol=1) l cvalue
      compute lastlogl=lrbase
   end do l
   *
   if print {
      if criterion==5
         report(action=format,atcol=2,atrow=bestlag+1,special=onestar)
      else
         report(action=format,atcol=2,tag=minimum,special=onestar)
      report(action=format,atcol=2,width=10,align=decimal)
      report(action=show)
   }
   compute %%autop=bestlag
}
if .not.%defined(model)
   return
*
* Set up a VAR model in standard order
*
system(model=model)
variables list
if bestlag>0
   lags 1 to bestlag
if ndet>0
   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. This reorders them.
*
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))
end



