*
* ACF performs an autocorrelation analysis on a series. This is a simpler
* routine than the related BJIDENT and REGCORRS.
*
* Syntax:
*
*  @ACF( options )  series start end
*
* Parameters:
*
*    series     series to analyze
*    start end  range to end (by default, the full range of series)
*
*  Options:
*     NUMBER=number of autocorrelations to compute [depends upon number of data points]
*     [GRAPH]/NOGRAPH Graph the correlations
*     REPORT/[NOREPORT]
*     METHOD=YULE/[BURG]
*     TITLE=Title of graph window and report window
*     HEADER=header string for graph
*     FOOTER=footer string for graph
*     PRINT/[NOPRINT] Include printed output from the correlations
*     [QSTATS]/NOQSTATS Computes and displays Ljung-Box Q Statistics
*     DFC=degrees of freedom correction for Q statistics
*
*  Revision Schedule:
*    05/2005 Stripped out of the more involved REGCORRS procedure
*    06/2007 Add FOOTER option
*    05/2014 Modernize to use report title, put header and footer on
*            the SPGRAPH rather than the GRAPH.
*
procedure ACF series start end
type series    series
type integer   start end
*
option switch  graph     1
option switch  qstats    1
option integer number
option choice  method    2   yule burg
option string  title
option string  header
option string  footer
option switch  report    0
option integer dfc       0
option switch  print     0
*
local string   ltitle
local report   areport
*
local string   qstring
local integer  ncorr i
local series   corrs pcorrs upper lower
local real     qstat
*
if .not.%defined(series) {
   disp "Syntax: @ACF( options )  series  start  end"
   return
}

clear corrs
corr(number=number,qstats,print=print,partial=pcorrs,method=method,dfc=dfc) series start end corrs
disp(store=qstring) "Q=" *.## %cdstat "P-value" #.##### %signif
inquire(series=corrs) * ncorr

if %defined(title)
   compute ltitle=title
else
if %defined(header)
   compute ltitle=header
else
   compute ltitle="Autocorrelation Function of "+%l(series)

if graph {
   set upper 1 ncorr = 2.0/sqrt(%nobs)
   set lower 1 ncorr = -2.0/sqrt(%nobs)
   *
   * We use SPGRAPH as a wrapper because we might be using GRTEXT to add the
   * final Q and the criteria to the graph.
   *
   spgraph(header=ltitle,footer=footer)
   graph(max=1.0,min=-1.0,number=1,style=bargraph,$
      overlay=line,ovcount=2,ovsame) 3
   # corrs 2 ncorr
   # upper 2 ncorr 1
   # lower 2 ncorr 1
   if qstats
      grtext(entry=2,y=-.9,align=left,size=16) qstring
   spgraph(done)
}
*
* Do the report of the correlations, partial correlations, running LB Q
* statistics and their significance levels.
*
if report {
   report(use=areport,action=define,hlabels=||"Lag","Corr","Partial","LB Q","Q Signif"||)
   compute qstat=0.0
   do i=2,ncorr
      compute qstat=qstat+(%nobs*(%nobs+2.0)/(%nobs-i+1))*corrs(i)**2
      report(use=areport,row=new,atcol=1) i-1 corrs(i) pcorrs(i) qstat
      if i-1-dfc>0
         report(use=areport,row=current,atcol=5) %chisqr(qstat,i-1-dfc)
   end do i
   report(use=areport,action=format,atcol=2,tocol=3,width=6)
   report(use=areport,action=format,atcol=4,width=9)
   report(use=areport,action=format,atcol=5,width=6)
   report(use=areport,action=show,window=ltitle)
}
end acf
