*
*  BJIDENT
*  Box-Jenkins model identification step
*
*  Syntax:
*  @BJIDENT  series   start   end
*
*  Options
*     DIFF=maximum regular differencings[0]
*     SDIFFS=maximum seasonal differencings[0]
*     TRANS=[NONE]/LOG/ROOT
*        Transformation to apply to data
*     [GRAPH]/NOGRAPH
*        Do High-resolution graphs?
*     SPAN=seasonal span
*     METHOD=YULE/[BURG]
*     NUMBER=number of autocorrelations to compute [depends upon data]
*     REPORT/[NOREPORT] REPORT shows a table with the ACF's and PACF's
*     QSTATS/[NOQSTATS] When used with REPORT, QSTATS includes the Q-statistics
*         in the report.
*
*  Written by Tom Doan, Estima.
*  Revision schedule:
*    05/99 Adjustments to improve appearance of graphs, adjust numbers
*          of correlations to match input.
*    01/03 Add METHOD option to allow either Burg (default) or Yule-Walker
*          correlation calculations, NUMBER option to control the number
*          of autocorrelations and PRINT to allow both graphing and printint.
*    01/04 Updated to use REPORT
*
proc bjident 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 switch  graph  1
option switch  print  0
option integer span
option choice  method 2  yule burg
option integer number
option switch  report 0
option switch  qstats 0
*
local integer nbeg nend spanl i j
local series  corrs     partials
local series  transfrm  diffed     upper lower
local series  qstat     qsignif
local integer numc
local string  dispstr

*
*  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)

*
*  If doing seasonal differences, make sure we do two years worth of
*  lags.
*

if %defined(number)
   compute numc   = number
else if sdiff>0
   compute numc   = %imax(3*spanl,25)
else
   compute numc   = 25

*
*  If we have a short data set, don't compute quite so many lags
*

if numc>(nend-nbeg+1)/2
   compute numc = (nend-nbeg+1)/2

*
*  Loop over differences and seasonal differences. Difference data
*  and compute correlations and partial autocorrelations. Graph them
*  if requested, otherwise print.
*

if report {
   report(action=define)
   report(atcol=1,atrow=1,fillby=cols) 'Lag' %seq(1,numc)
}
do i=0,diff
   do j=0,sdiff
      diff(diff=i,sdiff=j) transfrm nbeg+i+spanl*j nend diffed
      if sdiff==0
         disp(store=dispstr) i 'Differences of' %l(series)
      else
         disp(store=dispstr) i 'Regular' j 'Seasonal of' %l(series)
      corr(print=.not.graph.or.print,number=numc,partial=partials,stderrs=upper,method=method) $
           diffed nbeg+i+spanl*j nend corrs
      set upper 2 number+1 = 2*upper
      set lower 2 number+1 = -upper
      if (graph) {
         graph(max=1.0,min=-1.0,header=dispstr,style=bargraph,$
            number=0,nodates,window=dispstr,overlay=step,ovcount=2,key=loright,noovkey,ovsame) 4
         # corrs    1 numc+1
         # partials 1 numc+1
         # upper    2 numc+1 1
         # lower    2 numc+1 1
      }
      if report {
         set(first=0.0) qstat   1 numc+1 = qstat{1}+%nobs*(%nobs+2.0)/(%nobs-(t-1))*corrs(t)**2
         set(first=0.0) qsignif 1 numc+1 = %chisqr(qstat,(t-1))
         report(col=new,atrow=1,fillby=cols,align=center) 'Corr'     %xsubvec(corrs,2,numc+1)
         report(col=new,atrow=1,fillby=cols,align=center) 'Partial'  %xsubvec(partials,2,numc+1)
         if (qstats) {
         report(col=new,atrow=1,fillby=cols,align=center) 'LB Q'     %xsubvec(qstat,2,numc+1)
         report(col=new,atrow=1,fillby=cols,align=center) 'Q Signif' %xsubvec(qsignif,2,numc+1)
         }
      }
   enddo j
enddo i
if report {
   report(action=format,atcol=2,width=9)
   report(action=show,window='Correlations for '+%l(series))
}
end




