*
* VARLAG.RPF
* RATS Version 8, User's Guide, Example 7.1.
*
open data haversample.rat
calendar(q) 1959
data(format=rats) 1959:1 2006:4 dgnp fm2 gdph ih lr
*
set loggdp = log(gdph)
set logm2  = log(fm2)
set logp   = log(dgnp)
set logi   = log(ih)
set unemp  = lr
*
system(model=usamodel)
variables loggdp unemp logi logp logm2
lags 1 to 8
det constant
end(system)
estimate(noprint,resids=resids8) 1961:1 *
compute logdetu=%logdet
compute mcorr=%nreg,ntotal=%nregsystem
*
system(model=usamodel)
variables loggdp unemp logi logp logm2
lags 1 to 4
det constant
end(system)
estimate(noprint,resids=resids4) 1961:1 *
compute logdetr=%logdet
*
* Using RATIO
*
compute nrestr=ntotal-%nregsystem
ratio(degrees=nrestr,mcorr=mcorr,title="Test of 8 vs 4 Lags")
# resids8
# resids4
*
* Using calculated statistic
*
cdf(title="Test of 8 vs 4 Lags") chisqr $
     (%nobs-mcorr)*(logdetr-logdetu) nrestr
*
* Using VARLagSelect procedure
*
@varlagselect(lags=8,crit=aic)
# loggdp unemp logi logp logm2
*
* This shows the workings of VARLagSelect, showing columns for both AIC
* and SBC, and also shows sequential LR tests. Note that this uses a
* version of AIC with small-sample correction.
*
report(action=define,$
  hlabel=||"Lags","AIC","SBC","LR Test","P-Value"||)
dec real lastll
do lags=1,8
   system(model=usamodel)
   variables loggdp unemp logi logp logm2
   lags 1 to lags
   det constant
   end(system)
   estimate(noprint) 1961:1 *
   compute ll =%logl
   compute sbc=-2.0*ll/%nobs+%nregsystem*log(%nobs)/%nobs
   compute aic=-2.0*ll/%nobs+$
     %nregsystem*2.0*%nvar/(%nvar*%nobs-%nregsystem-1)
   report(row=new,atcol=1,align=decimal) lags aic sbc
   if lags>1
      report(row=current,atcol=4) 2*(ll-lastll) $
        %chisqr(2*(ll-lastll),%nvar*%nvar)
   compute lastll=ll
end do lags
report(action=format,atcol=2,tocol=2,special=onestar,$
  tag=min,align=decimal)
report(action=format,atcol=3,tocol=3,special=onestar,$
  tag=min,align=decimal)
report(action=format,atcol=2,tocol=3,width=8)
report(action=format,atcol=4,tocol=5,picture="*.####")
report(action=show)