*
* @ARAutoLags( options )  series  start   end
*
* ARAutoLags computes an information criterion for various lags of AR processes
* using the Burg or Yule estimates of the partial autocorrelations.
*
* Parameters:
*    series       series to analyze
*    start end    range of series to use (by default, the maximum possible)
*
* Options:
*    method=[yule]/burg
*    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:
*   02/2009 Written by Tom Doan, Estima
*
proc ARAutoLags series start end
type series    series
type integer   start end
*
option integer maxlags 25
option choice  method  1  yule burg
option choice  crit    1  aic bic hq caic
option switch  table   0
*
local series   xc pc
local vect     ic
local real     sigmasq fiddle
local integer  i maxl
local string   critlabel
*
declare integer %%autop
*
* 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 partial autocorrelations and covariances
*
corr(number=maxl,covariances,method=method,noprint,partial=pc) series start end xc
compute sigmasq=xc(1)
*
* Figure out the IC penalty (per free parameter).
* 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(sigmasq)+fiddle
*
* These count the extracted mean as a parameter
*
do i=1,maxl
   *
   * Adjust the multiplier for the corrected AIC
   *
   if crit==4 ; compute fiddle=fiddle*(%nobs-i)/(%nobs-(i+1))
   compute sigmasq=sigmasq*(1-pc(i+1)**2)
   compute ic(i+1)=log(sigmasq)+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



