*
* 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.
*
* Syntax:
*
*    @yulelags( options )  series  start   end
*
* 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]
*    type=[aic]/bic
*
proc yulelags series start end
type series series
type integer start end
*
option integer maxlags 25
option choice  type    1  aic  bic
*
local series xc
local symm cxx
local vect cxy ic z
local rect s
local real rss fiddle
local integer i maxl
*
* Compute the needed covariances
*
corr(number=maxlags+1,covariances,noprint) series / xc
*
* Use T/2 as the maximum lags if user provides nothing
*
compute maxl=maxlags
if %nobs/2<maxl
   compute maxl=%nobs/2
*
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.
*
compute fiddle=%if(type==1,2.0/%nobs,log(%nobs)/%nobs)
dim ic(maxl+1)
compute ic(1)=log(rss)+fiddle
do i=1,maxl
   compute rss=rss-z(i)**2
   compute ic(i+1)=log(rss)+fiddle*(i+1)
end do i
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,special=onestar,tag=minimum,align=decimal)
report(action=format,atcol=2,width=9)
report(action=show,window='IC Lag Analysis for '+%l(series))
end


