The RATS Software Forum

This quickly scans a set of moving average models looking for the one which minimizes one of four information criteria. This can be used in ARMA model identification, and also for diagnostics (best lag should be zero if the series is serially uncorrelated).

Code: Select all

*
* @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
disp theta
*
* 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