*
* @CUMPDGM( options )  series  start  end
*
* Performs a Durbin's Cumulated Periodogram for serial correlation. The
* test statistic is a Kolmogorov-Smirnov statistic.
*
* Parameters:
*   series      series to test
*   start  end  range to test (by default, the defined range of <<series>>)
*
* Options:
*   [PRINT]/NOPRINT
*   [GRAPH]/NOGRAPH
*
* Reference:
*   J. Durbin(1969), "Tests for Serial Correlation in Regression Analysis Based
*     on the Periodogram of Least Squares Residuals", Biometrika, 1-16.
*
* Revision Schedule:
*   01/1992 Converted to RATS 4.0, Rejection Limits Added
*   01/1995 Added symbol font Xlabels, GRTEXT
*   10/2007 Added GRPARM save and recall
*   08/2010 Switch to use of REPORT. Add PRINT and GRAPH options.
*
procedure cumpdgm series start end
type series 	series
type integer 	start end
*
option switch	print	1
option switch	graph 1
*
local integer 	nobs startl endl ordinate half
local series  	white_noise actual freqs
local string  	gaplabel
local vector[integer] grparms
local report	creport
local string	ltitle

if .not.%defined(series) {
   display "Syntax: @CUMPDGM   series  start  end"
   return

}

inquire(series=series)  startl<<start   endl<<end

compute nobs=endl-startl+1
compute ordinate=%freqsize(nobs)

compute half = (ordinate+1)/2

frequency 3 ordinate
rtoc(pad)  startl  endl  1
# series
#   1
fft 1
cmult 1 1
cacc(scale) 1 1 half
cunits 2
cacc(scale) 2 1 half
cset 3 1 half = %z(t,1)-%z(t,2)
cxt(noprint,part=absval) 3 1 half

compute ltitle="Cumulated Periodogram Test for Series "+%l(series)
report(use=creport,action=define,title=ltitle)
report(use=creport,atrow=1,atcol=1,span) ltitle
report(use=creport,atrow=2,atcol=1) "Maximum Gap" %maximum
report(use=creport,atrow=3,atcol=1) "At Frequency"  2*%pi*float(%maxent)/ordinate
report(use=creport,atrow=5,atcol=1,span) "Approximate Rejection Limits"
report(use=creport,atrow=6,atcol=1) "1%" 1.63/SQRT(HALF)
report(use=creport,atrow=7,atcol=1) "5%" 1.36/SQRT(HALF)
report(use=creport,atrow=8,atcol=1) "10%" 1.22/SQRT(HALF)
report(use=creport,action=format,picture="#.####")
if print
   report(use=creport,action=show)
if .not.graph
   return
*
*
*
ctor 1 half
# 1 2
# actual white_noise
label actual white_noise
# "Actual" "White Noise"
set freqs 1 half = t-1
*
compute grparms = %grparm()
grparm(font="Symbol") axislabels *
spgraph(header="Cumulated Periodogram Test")
scatter(style=line,key=upleft,xlabels=||"0","p/4","p/2","3p/4","p"||) 2
# freqs white_noise 1 half
# freqs actual      1 half
display(store=gaplabel) "gap = " #.#### %maximum
if actual(%maxent)<white_noise(%maxent)
   grtext(alignment=left,x=%maxent-1,y=actual(%maxent)-.01) gaplabel
else
   grtext(alignment=right,x=%maxent-1,y=actual(%maxent)+.01) gaplabel
spgraph(done)
grparm(recall=grparms)
end

