*
* @BJIDENT  series   start   end
* Box-Jenkins model identification step
*
* Parameters:
*   series      dependent variable
*    start  end  range over which to estimate (by default, range of series)
*
* Options
*   DIFF=maximum regular differencings[0]
*   SDIFFS=maximum seasonal differencings[0]
*   TRANS=[NONE]/LOG/ROOT
*     Transformation to apply to data
*
*   [GRAPH]/NOGRAPH
*     Do graphs?
*
*   SPAN=seasonal span
*   NUMBER=number of autocorrelations to compute [default depends upon data]
*
*   METHOD=YULE/[BURG]
*     Selects method used to compute correlations (see CORR instruction for details)
*
*   REPORT/[NOREPORT]
*   QSTATS/[NOQSTATS]
*     REPORT shows a table with the ACF's and PACF's
*     When used with REPORT, QSTATS includes the Q-statistics in the report.
*
*   SPIKE/[NOSPIKE]
*     SPIKE selects a "spike" graph, rather than a bar graph, with the
*     regular and partial autocorrelation graphs separated vertically.
*
*   SEPARATE/NOSEPARATE  [default is NOSEPARATE except when using SPIKE]
*     Use SEPARATE if you want the correlations and partial correlations
*     plotted as separate graphs (one a single page).
*
*   [ZERO]/NOZERO
*     Display the first correlation
*
*   MAX=maximum/minimum value on vertical axis scale [1]
*
*  Revision schedule (written by Tom Doan, Estima):
*    05/1999 Adjustments to improve appearance of graphs, adjust numbers
*            of correlations to match input.
*    01/2003 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/2004 Updated to use REPORT
*    03/2007 Added MAX and ZERO options
*    04/2007 Added SPIKE, SEPARATE options.
*    09/2010 Add TITLE to report, use local report
*    06/2014 Add header over group of statistics in report
*    04/2017 Fix formatting in report (could misalign numbers with -)
*
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
option real    max    1.0
option switch  zero   1
option switch  spike  0
option switch  separate spike
*
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
local integer 	startper numberarg styleopt
local integer  leftgroup rightgroup
local string   ltitle
local	report   creport
*
*  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


* Set argument for NUMBER option and the graph start period
* based on whether or not the ZERO option is on:

if zero
   compute numberarg=0, startper=1
  else
   compute numberarg=1, startper=2

*
* Set STYLE option depending on whether or not user selected SPIKE.
* (STYLE = 3 for bar, STYLE = 12 for spike):
*
if spike
   compute styleopt = 12
else
   compute styleopt = 3

*
*  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.
*
compute ltitle="Correlations for "+%l(series)

report(use=creport,action=define,title=ltitle)
report(use=creport,atcol=1,atrow=2,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.and.separate) {
            spgraph(vfields=2,header=dispstr)
             graph(max=max,min=-max,style=styleopt,header="Correlations",$
                number=numberarg,nodates,window=dispstr,overlay=step,ovcount=2,noovkey,ovsame) 3
             # corrs    startper numc+1
             # upper    2        numc+1 1
             # lower    2        numc+1 1
             graph(max=max,min=-max,style=styleopt,header="Partial Correlations",$
                number=numberarg,nodates,window=dispstr,overlay=step,ovcount=2,noovkey,ovsame) 3
             # partials startper numc+1
             # upper    2        numc+1 1
             # lower    2        numc+1 1
            spgraph(done)
      }
      if (graph.and.(.not.separate)) {
             graph(max=max,min=-max,header=dispstr,style=styleopt,$
                number=numberarg,nodates,window=dispstr,overlay=step,ovcount=2,key=loright,noovkey,ovsame) 4
             # corrs    startper numc+1
             # partials startper numc+1
             # upper    2 numc+1 1
             # lower    2 numc+1 1
      }

      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))
      compute leftgroup=%reportcol+1
      report(use=creport,col=new,atrow=2,align=center)    "Corr"
      report(use=creport,col=current,atrow=3,fillby=cols) %xsubvec(corrs,2,numc+1)
      report(use=creport,col=new,atrow=2,align=center)    "Partial"
      report(use=creport,col=current,atrow=3,fillby=cols) %xsubvec(partials,2,numc+1)
      if (qstats) {
         report(use=creport,col=new,atrow=2,align=center)    "LB Q"
         report(use=creport,col=current,atrow=3,fillby=cols) %xsubvec(qstat,2,numc+1)
         report(use=creport,col=new,atrow=2,align=center)    "Q Signif"
         report(use=creport,col=current,atrow=3,fillby=cols) %xsubvec(qsignif,2,numc+1)
      }
      compute rightgroup=%reportcol
      report(use=creport,atrow=1,atcol=leftgroup,tocol=rightgroup,span,align=center) dispstr
   enddo j
enddo i
report(use=creport,action=format,atcol=2,width=9)
if report
   report(use=creport,action=show)
end
