*
* @MAAutoLags( options) series  start   end
*
* MAAutoLags computes an information criterion for various lags of MA
* processes using the innovations algorithm, finding the "optimal"
* choice under that criterion and (optionally) producing a table with
* the values for all lags.
*
* 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:
*    %%autoq = best number of lags
*
* Revision Schedule:
*   02/2009 Written by Tom Doan, Estima
*
proc MAAutoLags series start end
type series    series
type integer   start end
*
option integer maxlags 25
option choice  crit    1  aic bic hq caic
option switch  table   0
*
type series    series
type integer   start end
*
local symm     kappa
local rect     theta
local vect     v
local real     s
local series   covarx
local integer  i j n k
local integer  startl endl
local vect     ic
local real     fiddle
local integer  i maxl
local string   critlabel
*
declare integer %%autoq
*
* 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
inquire(series=series) startl<<start endl<<end
*
dim kappa(maxl+1,maxl+1) theta(maxl,maxl) v(maxl+1)
corr(covariances,number=maxl+1,noprint,method=yule) series startl endl covarx
*
* Create the Toeplitz matrix kappa
*
ewise kappa(i,j)=covarx(i-j+1)
*
compute v(1)=kappa(1,1)
compute theta=%const(0.0)
*
* 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"
*
do n=1,maxl
   do k=0,n-1
      compute s=kappa(n+1,k+1)
      do j=0,k-1
         compute s=s-theta(k,k-j)*theta(n,n-j)*v(j+1)
      end do j
      compute theta(n,n-k)=s/v(k+1)
   end do k
   compute s=kappa(n+1,n+1)
   do j=0,n-1
      compute s=s-theta(n,n-j)**2*v(j+1)
   end do j
   compute v(n+1)=s
end do n
*
* 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"
*
* These count the extracted mean as a parameter
*
dim ic(maxl+1)
do j=1,maxl+1
   compute ic(j)=log(v(j))+fiddle*j
   *
   * Adjust the multiplier for the corrected AIC
   *
   if crit==4 ; compute fiddle=fiddle*(%nobs-i)/(%nobs-(j+1))
end do j
*
compute %%autoq=%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)+" =" %%autoq
end

