*
* @BJAutoFit(pmax=maximum AR,qmax=maximum MA, other options ) y  start end
*
* BJAutofit chooses the model for the series <<y>> which minimizes one
* of the information criteria. It uses maximum likelihood estimation to
* ensure that the estimates are done over a consistent time interval.
*
* Options:
*   PMAX=maximum number of AR lags to consider [0]
*   QMAX=maximum number of MA lags to consider [0]
*   CRIT=[AIC]/BIC/SBC/HQ/CAIC criterion to use
*     AIC is (uncorrected) Akaike information criterion.
*     BIC and SBC are synonyms for Bayesian or Schwarz criterion.
*     HQ is Hannan-Quinn.
*     CAIC is the AIC-corrected for degrees of freedom
*   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]
*   [PRINT]/NOPRINT controls whether the table is displayed
*   TITLE=title for table [title=title for table ["'criterion' analysis of 'series'"]
*
* Variables Defined:
*   %%AUTOP  number of AR parameters selected
*   %%AUTOQ  number of MA parameters selected
*
* Revision Schedule:
*   04/2005 Written by Tom Doan, Estima
*   06/2008 Added PRINT/NOPRINT option. Defined %%autop and %%autoq
*   01/2009 Add CRIT=CAIC (for corrected AIC)
*   07/2010 Switch to use of LOCAL REPORT. Add TITLE option.
*   12/2020 Add QUIET option to suppress use of any REPORT features
*
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 hq caic
option switch  demean 0
option switch  const  0
option integer diffs  0
option integer sdiffs 0
option switch  print  1
option string  title
option switch  quiet  0
*
local integer  p q
local string   clabel
local real     ic bestic penalty
local string   ltitle
local report   ireport
*
declare integer %%autop %%autoq
*
if .not.%defined(y) {
   display "##Syntax: @BJAutoFit(pmax=max AR,qmax =max MA,other options) <<y>> start end"
   return
}

if crit==1
   compute clabel="AIC"
else
   compute clabel="BIC"

if %defined(title)
   compute ltitle=title
else
   compute ltitle=clabel+" analysis of models for series "+%l(y)

if (.not.quiet) {
   report(use=ireport,action=define,title=ltitle)
   report(use=ireport,row=new,atcol=1,span) ltitle
   report(use=ireport,row=new,atcol=2,span) "MA"
   report(use=ireport,row=new,atcol=1) "AR" %seq(0,qmax)
}
compute %%autop=%%autoq=0
do p=0,pmax
   if (.not.quiet)
      report(use=ireport,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
      if crit==5
         compute penalty=2*(p+q+1)*%nobs/(%nobs-p-q-2)
      else
      if crit==1
         compute penalty=(p+q)*2.0
      else
      if crit==4
         compute penalty=(p+q)*2.0*log(log(%nobs))
      else
         compute penalty=(p+q)*log(%nobs)
      compute ic=-2*%logl+penalty
      if p==0.and.q==0.or.ic<bestic
         compute %%autop=p,%%autoq=q,bestic=ic
      if .not.quiet
         report(use=ireport,row=current,atcol=q+2) ic
   end do q
end do p
if .not.quiet {
   report(use=ireport,action=format,atcol=2,tocol=qmax+2,picture="*.####")
   report(use=ireport,action=format,tag=minimum,atcol=2,tocol=qmax+2,align=decimal,special=one)
   if print
      report(use=ireport,action=show)
}
end bjautofit

