*
* @yulelags( options )  series  start   end
*
* YuleLags computes an information criterion for various lags of AR
* processes using the Yule-Walker estimates based upon the sample
* covariances. The Yule-Walker algorithm is particularly good in this
* setting, because all the sums of squares can be derived from a single
* matrix operation.
*
* The newer procedure ARAutoLags is recommended instead.
*
* Parameters:
*    series       series to analyze
*    start end    range of series to use (by default, the maximum possible)
*
* Options:
*    maxlags=maximum number of lags to consider [25]
*    crit=[aic]/bic/caic/hq
*       Criterion to use. CAIC is AIC corrected for degrees of freedom
*         (usually called AICC).
*    table/[notable]   Show full table of results (not just best)
*
* Variables Defined:
*    %%autop = best number of lags
*
* Revision Schedule:
*   05/2005 Written by Tom Doan, Estima
*   02/2009 Changed TYPE option to CRIT option and added corrected AIC and
*           HQ. Defined %%autop. Improved output. Newer procedure ARAUTOLAGS
*           written.
*
proc yulelags series start end
type series series
type integer start end
*
option integer maxlags 25
option choice  crit    1  aic  bic  caic  hq
option switch  table   0
*
local series   xc
local symm     cxx
local vect     cxy ic z
local rect     s
local real     rss fiddle
local integer  i j maxl
*
* Don't go above T/2 for the maximum lags
*
compute maxl=maxlags
stats(noprint) series start end
if %nobs/2<maxl
   compute maxl=%nobs/2
*
* Compute the needed covariances
*
corr(number=maxl+1,covariances,method=yule,noprint) series start end xc
dim cxx(maxl,maxl)
ewise cxx(i,j)=xc(i-j+1)
dim cxy(maxl)
ewise cxy(i)=xc(i+1)
*
compute s=inv(%decomp(cxx))
compute z=s*cxy
*
compute rss=xc(1)
*
* Figure out the IC penalty (per free parameter). Note that, because the
* covariances are already scaled through by the number of observations,
* the "rss" computed here is actually the sigma squared.
*
* We use a /T variation on the criteria.
*
if crit==1
   compute fiddle=2.0/%nobs,critlabel="AIC"
else
if crit==2
   compute fiddle=log(%nobs)/%nobs,critlabel="Schwarz/Bayesian IC"
else
if crit==3
   compute fiddle=2.0*log(log(%nobs))/%nobs,critlabel="Hannan-Quinn"
else
   compute fiddle=2.0/(%nobs-1.0),critlabel="AIC-Corrected"
*
dim ic(maxl+1)
compute ic(1)=log(rss)+fiddle
do i=1,maxl
   if crit==4 ; compute fiddle=fiddle*(%nobs-i)/(%nobs-(i+1))
   compute rss=rss-z(i)**2
   compute ic(i+1)=log(rss)+fiddle*(i+1)
end do i
compute %%autop=%minindex(ic)-1
if table {
   report(action=define,hlabels=||"Lags","IC"||)
   report(atcol=1,atrow=1,fillby=cols) %seq(0,maxl)
   report(atcol=2,atrow=1,fillby=cols,align=decimal) ic
   *
   * Tag the minimum IC value
   *
   report(action=format,atcol=2,picture="*.###",width=10,special=onestar,tag=minimum,align=decimal)
   report(action=show,window=critlabel+" Lag Analysis for "+%l(series))
}
else
   display "Minimum "+critlabel+" Lags for "+%l(series)+" =" %%autop
end


