*
* Performs the battery of independence tests used extensively in
* Brockwell and Davis (2002), "Introduction to Time Series and
* Forecasting" 2nd ed, Springer.
*
* Syntax:
*
*  @bdindtests( options ) x start end
*
* Parameters:
*
*   x           series to analyze
*   start end   range to analyze [defaults to full range of x]
*
* Option:
*   NUMBER=number of autocorrelations in the Ljung-Box and McLeod-Li tests
*
* Revision Schedule:
*   05/05 Written by Tom Doan, Estima
*
proc BDIndTests x start end
type series x
type integer start end
*
option integer number
*
local real 		tp ds rt
local real 		etp eds ert vtp vds vrt ztp zds zrt
local real 		nobs
local integer 	s ncorr
local integer 	startl endl
local real 		lbstat lbsignif mlstat mlsignif
local series 	xsq
local string 	title
*
inquire(series=x) startl<<start endl<<end
*
compute nobs=endl-startl+1
*
* Count the sign changes
*
sstats startl+1 endl (x>x{1})>>ds
*
* Count the turning points
*
sstats startl+2 endl (x>x{1}.and.x{1}<x{2}).or.(x<x{1}.and.x{1}>x{2})>>tp
*
* Do the pair count for the rank test
*
compute rt=0
do t=startl,endl
   do s=t+1,endl
      compute rt=rt+(x(s)>x(t))
   end do s
end do t
*
* Compute the expected values and asymptotic variances of
* the test statistics
*
compute etp=2.0*(nobs-2)/3
compute vtp=(16.0*nobs-29)/90
compute eds=0.5*(nobs-1)
compute vds=(nobs+1)/12.0
compute ert=.25*nobs*(nobs-1)
compute vrt=nobs*(nobs-1)*(2.0*nobs+5)/72.0
*
corr(number=number,noprint,qstats) x startl endl
compute lbstat=%qstat,lbsignif=%qsignif
set xsq startl endl = x^2
corr(number=number,noprint,qstats) xsq startl endl
compute mlstat=%qstat,mlsignif=%qsignif
*
compute title="Independence Tests for Series "+%l(x)
report(action=define,title=title)
report(row=new,atcol=1,span) title
report(row=new,atcol=1) "Test" "Statistic" "P-Value"
report(row=new,atcol=1) "Ljung-Box Q("+%ndfq+")" lbstat lbsignif
report(row=new,atcol=1) "McLeod-Li("+%ndfq+")" mlstat mlsignif
report(row=new,atcol=1) "Turning Points" ztp=(tp-etp)/sqrt(vtp) %ztest(ztp)
report(row=new,atcol=1) "Difference Sign" zds=(ds-eds)/sqrt(vds) %ztest(zds)
report(row=new,atcol=1) "Rank Test" zrt=(rt-ert)/sqrt(vrt) %ztest(zrt)
report(action=format,atcol=2,width=10)
report(action=format,atcol=3,picture="#.####")
report(action=show)
end

