*
* BJAutofit chooses the minimum AICC or BIC model for the series y
* It uses maximum likelihood estimation to ensure that the estimates
* are done over a consistent time interval.
*
* Syntax:
*   @BJAutoFit(pmax=maximum AR,qmax=maximum MA, other options ) y  start end
*
* Options:
*   PMAX=maximum number of AR lags to consider [0]
*   QMAX=maximum number of MA lags to consider [0]
*   CRIT=[AIC]/BIC/SBC criterion to use (BIC and SBC are synonyms)
*   DEMEAN/[NODEMEAN] extract mean before estimating ARMA models
*   CONST/[NOCONST] estimate mean as part of model
*   DIFFS=number of preliminary differences [0]
*   SDIFFS=number of preliminary seasonal differences[0]
*
procedure bjautofit y start end
type series y
type integer start end
*
option integer pmax   0
option integer qmax   0
option choice  crit   1  aic bic sbc
option switch  demean 0
option switch  const  0
option integer diffs  0
option integer sdiffs 0
*
local integer p q
local string clabel
local real   penalty
*
if .not.%defined(y) {
   display 'Syntax: @BJAutoFit(pmax=max AR,qmax =max MA,other options) y start end'
   return
}

report(action=define)
if crit==1
   compute clabel='AIC'
else
   compute clabel='BIC'

report(row=new,atcol=1,span) clabel+' analysis of models for series '+%l(y)
report(row=new,atcol=2,span) 'MA'
report(row=new,atcol=1) 'AR' %seq(0,qmax)
do p=0,pmax
   report(row=new,atcol=1) p
   do q=0,qmax
      boxjenk(noprint,ar=p,ma=q,maxl,const=const,demean=demean,diffs=diffs,sdiffs=sdiffs,cvcrit=.001) y start end
*      compute penalty=%if(crit==1,2*(p+q+1)*%nobs/(%nobs-p-q-2),(p+q)*log(%nobs))
      compute penalty=%if(crit==1,2*(p+q),(p+q)*log(%nobs))
      report(row=current,atcol=q+2) -2*%logl+penalty
   end do q
end do p
report(action=format,atcol=2,tocol=qmax+2,picture='*.####')
report(action=format,tag=minimum,atcol=2,tocol=qmax+2,align=decimal,special=one)
report(action=show)
end bjautofit






