*
* @RSStatistic( options )  x  start  end
* Computes either the classical R/S statistic, or Lo's modified version.
*
* Note that because the R/S statistic is usually computed for many
* different samples, this procedure doesn't display any output. Instead,
* use it to do the number-crunching, and retrieve the result with %CDSTAT.
*
* Parameters:
*   x          Series to analyze
*   start end  Range to analyze [defined range of x]
*
* Options:
*   CLASSICAL/[NOCLASSICAL]
*     If CLASSICAL, it computes the classical form (range divided by sample
*     standard deviation). If NOCLASSICAL (the default) it computes Lo's
*     modified form, where the scale is the square root of the long-run
*     variance. Lo's statistic is also divided by the square root of the
*     number of observations.
*
*   LAGS=Bartlett window width for computing Lo's statistic [sqrt(observations)]
*
* Variables defined:
*  %CDSTAT is the computed R/S statistic
*
* References:
*   Mandelbrot and Wallis(1969), "Computer Experiments with Fractional
*    Gaussian Noise", Water Resources Res., vol 5, 228-267.
*
*   Lo(1991) "Long-term Memory in Stock Market Prices", Econometrica,
*    vol 59, 1279-1313.
*
* Revision Schedule:
*   02/2007 Written by Tom Doan, Estima
*
procedure RSStatistic x start end
type series 	x
type integer 	start end
*
option switch 	classical 0
option integer lags
*
local integer 	startl endl lagsl
local series 	dx sumdx
local real 		range stdvar lrvar
*
inquire(series=x) startl<<start endl<<end
diff(center) x startl endl dx
*
* Accumulate the changes and find the range.
*
acc dx startl endl sumdx
ext(noprint) sumdx startl endl
compute range=%maximum-%minimum
*
* Compute the two options for the scale factor.
*
stats(noprint) dx startl endl
compute stdvar=%variance*(%nobs-1)/%nobs

if .not.%defined(lags)
   compute lagsl=fix(sqrt(%nobs))
else
   compute lagsl=lags
*
* Because the matrix produced by MCOV is a sum (not an average) we can
* just use its square root directly to get Lo's statistic.
*
mcov(lags=lagsl,lwindow=bartlett) / dx
# constant
compute lrvar=%cmom(1,1)
*
if classical
   compute %cdstat=range/sqrt(stdvar)
else
   compute %cdstat=range/sqrt(lrvar)
end
