*
* @BJDIFF( options )  series   start   end
*
* Generates a table of Schwarz criteria for various combinations of differencings
* and mean extractions, as a preliminary step in Box-Jenkins modeling.
*
* Box-Jenkins model identification step
*
* Parameters
*   series      Series to analyze
*   start end   Range to analyze [range of series]
*
* Options
*   DIFF=maximum regular differencings[0]
*   SDIFFS=maximum seasonal differencings[0]
*   TRANS=[NONE]/LOG/ROOT
*      Transformation to apply to data
*   SPAN=seasonal span
*   WIDTH=window width for spectral smoothing
*
* Variables Defined:
*   %%AUTOD     Recommended number of regular differencings
*   %%AUTODS    Recommended number of seasonal differencings
*   %%AUTOCONST Recommended choice for CONST option
*
* Revision schedule:
*   05/1999 Written by Tom Doan, Estima
*   07/2005 Revised to use non-parametric estimate of the variance
*   01/2009 Add WIDTH option to control spectral smoothing
*   01/2009 Correct %NOBS in NPWNVariance.
*   02/2009 Define %%AUTOD, %%AUTODS and %%AUTOCONST
*   02/2009 Flip sign on criterion
*

* @NPWNVariance( options )  series  start end
*
* Computes a non-parametric (frequency domain) estimate of the white noise
* variance for the moving average representation of a process
*
* Option:
*   [DEMEAN]/NODEMEAN   extract the mean before estimating the spectral density
*   WIDTH=window width for spectral smoothing
*
procedure  NPWNVariance  series  start  end
type series    series
type integer   start end
*
option switch  demean 1
option integer width
*
local   int    startl endl widthl ordinate
local   series workx
local   real   integral
*
declare real   %wnvariance
*
inquire(series=series) startl<<start endl<<end

compute %nobs    = endl-startl+1
compute ordinate = 2*%freqsize(%nobs)
if %defined(width)
   compute widthl = width
else
   compute widthl = ordinate/4
*
*  Remove the mean
*

if demean
   diff(center) series startl endl workx
else
   set workx startl endl = series

freq 4 ordinate
rtoc startl endl
#  workx
#    1

*
*  Compute the log spectral density
*

fft 1
cmult(scale=1.0/%nobs) 1 1 / 2
cset 4 = 1.0
window(width=widthl,mask=4) 2 / 3
cln  3
*
*  Compute its integral over [-pi,pi]
*
sstats 1 ordinate %real(%z(t,3))>>integral
compute %wnvariance=exp(integral/(ordinate))
compute %nobs    = endl-startl+1
end
*
* BJDIFF procedure itself
*
proc bjdiff series start end
type series    series
type integer   start end
*
option integer diff   0
option integer sdiff  0
option choice  trans  1  none  log  root
option integer span
option integer width
*
local integer  nbeg nend spanl i j
local series   transfrm  diffed
local real     sbc bestsbc
*
dec   integer  %%autod %%autods %%autoconst
*
* Get the series range and seasonal, to be overridden by the start and end
* parameters and span option if the user specifies them.
*
inquire(series=series) nbeg<<start nend<<end
inquire(seasonal) spanl<<span
*
* Transform the series
*
if trans == 1
   set transfrm nbeg nend = series
else if trans == 2
   set transfrm nbeg nend = log(series)
else
   set transfrm nbeg nend = sqrt(series)

report(action=define)
report(atrow=1,atcol=1,span) "BJDiff Table, Series "+%l(series)
report(atrow=2,atcol=1) "Reg Diff" "Seas Diff" "Intercept" "Crit"
compute %%autod=-1
do j=0,diff
   do i=0,sdiff
      diff(diffs=j,sdiffs=i,span=spanl) transfrm nbeg+j+i*spanl nend diffed
      @NPWNVariance(width=width,nodemean) diffed nbeg+j+i*spanl nend
      compute sbc=log(%wnvariance)
      if %%autod==-1.or.sbc<bestsbc
         compute %%autod=j,%%autods=i,%%autoconst=0,bestsbc=sbc
      report(row=new,atcol=1,align=decimal) j i "No" sbc
      @NPWNVariance(width=width,demean) diffed nbeg+j+i*spanl nend
      compute sbc=log(%wnvariance)+log(%nobs)/%nobs
      if sbc<bestsbc
         compute %%autod=j,%%autods=i,%%autoconst=1,bestsbc=sbc
      report(row=new,atcol=1,align=decimal) j i "Yes" sbc
   end do i
end do j
*
report(action=format,atcol=4,width=9)
report(action=format,tag=minimum,special=onestar,align=decimal)
report(action=show)
end bjdiff




