*
* @GMAutoFit(options)  series start end
* does an automatic fit of a multiplicative seasonal ARMA model to a series
* following the procedure described in Gomez and Maravall, "Automatic Modeling
* Methods for Univariate Series", in Peña, Tiao and Tsay, eds., "A Course in Time
* Series Analysis", New York: Wiley, 2001.
*
* Parameters:
*   series      Series to analyze
*   start end   Range to analyze [range of series]
*
* Options
*   DIFF=number of regular differencings[0]
*   SDIFFS=number of seasonal differencings[0]
*   CONST/[NOCONST]  Include constant in all models
*   REGULAR=maximum number of regular AR or MA parameters (each polynomial) [2]
*   SEASONAL=maximum number of seasonal AR or MA parameters (each polynomial) [1]
*   FULL/[NOFULL] FULL does an exhaustive search rather than the GM search
*   REPORT/[NOREPORT] shows a table with all the models examined
*
* Description:
*   @GMAUTOFIT searches for an approximate minimum BIC model with no more than
*   the
*
* Defines:
*   %%AUTOP,%%AUTOQ,%%AUTOPS,%%AUTOQS are the chosen number of regular AR,
*    regular MA,seasonal AR and seasonal MA lags
*
procedure GMAutoFit series  start end
type series    series
type integer   start end
*
option integer diffs    0
option integer sdiffs   0
option integer regular  2
option integer seasonal 1
option switch  const    0
option switch  report   0
option switch  full     0
*
local real      ic bestic
local integer   p q ps qs
local integer   blockrow
local integer   plower pupper qlower qupper
local report    gmreport
*
declare integer %%autop %%autoq %%autops %%autoqs
if report {
   report(use=gmreport,action=define)
   report(use=gmreport,atrow=1,atcol=1,span) "Search for Minimum BIC Model"
   report(use=gmreport,atrow=2,atcol=1,span) "Series "+%label(series)
   report(use=gmreport,atrow=3,atcol=1,span) "with "+diffs+" regular and "+sdiffs+" seasonal differences"
   report(use=gmreport,atrow=5,atcol=1,align=center) "AR" "MA" "AR(s)" "MA(s)" "LogL" "BIC"
}
*
* For the GM search, fix the regular ARMA model at (3,0) and do an exhaustive
* search over the seasonal parameters.
*
if full
   compute plower=0,pupper=regular,qlower=0,qupper=regular,blockrow=%reportrow+1
else
   compute plower=3,pupper=3,qlower=0,qupper=0
*
compute bestic=%na
do p=plower,pupper
   do q=qlower,qupper
      do ps=0,seasonal
         do qs=0,seasonal
            boxjenk(ar=p,ma=q,sar=ps,sma=qs,$
               diffs=diffs,sdiffs=sdiffs,const=const,cvcrit=.001,noprint,maxl,method=bfgs) series start end
            compute ic=-2.0*%logl+%nreg*log(%nobs)
            if report ; report(use=gmreport,row=new,atcol=1) p q ps qs %logl ic
            if .not.%valid(bestic).or.ic<bestic
               compute %%autop=p,%%autoq=q,%%autops=ps,%%autoqs=qs,bestic=ic
         end do qs
      end do ps
   end do q
end do p
if full ; branch final
*
* Skip a line between blocks
*
if report ; report(use=gmreport,row=new)
*
* Fix the seasonal ARMA model at the optimal values from the previous step.
* Do an exhaustive search over the regular parameters
*
compute ps=%%autops,qs=%%autoqs
do p=0,regular
   do q=0,regular
      boxjenk(ar=p,ma=q,sar=ps,sma=qs,$
        diffs=diffs,sdiffs=sdiffs,const=const,cvcrit=.001,noprint,maxl,method=bfgs) series start end
      compute ic=-2.0*%logl+%nreg*log(%nobs)
      if report ; report(use=gmreport,row=new,atcol=1) p q ps qs %logl ic
      if p==0.and.q==0.or.ic<bestic
         compute %%autop=p,%%autoq=q,bestic=ic
   end do q
end do p
*
* Skip a line between blocks
*
if report ; report(use=gmreport,row=new)
*
* Fix the regular ARMA model at the optimal values from the previous step and
* search again over the seasonal parameters
*
compute blockrow=%reportrow+1
compute p=%%autop,q=%%autoq
do ps=0,seasonal
   do qs=0,seasonal
      boxjenk(ar=p,ma=q,sar=ps,sma=qs,$
        diffs=diffs,sdiffs=sdiffs,const=const,cvcrit=.001,noprint,maxl,method=bfgs) series start end
      compute ic=-2.0*%logl+%nreg*log(%nobs)
      if report ; report(use=gmreport,row=new,atcol=1) p q ps qs %logl ic
      if ps==0.and.qs==0.or.ic<bestic
         compute %%autops=ps,%%autoqs=qs,bestic=ic
   end do qs
end do ps
*
:final
if report {
   report(use=gmreport,action=format,atcol=1,tocol=4,width=6)
   report(use=gmreport,action=format,atcol=5,tocol=5,width=11)
   report(use=gmreport,action=format,atcol=6,tocol=6,atrow=blockrow,tag=minimum,special=onestar)
   report(use=gmreport,action=format,atcol=6,tocol=6,width=11,align=decimal)
   report(use=gmreport,action=show)
}
end




