*
* @vratio(options) series start end
*
* Performs a variance ratio test as described in Campbell, Lo and
* MacKinlay, The Econometrics of Financial Markets, section 2.4.
*
* Parameters:
*   series       The series to be tested. This is the *undifferenced* series, such
*                as the log price series.
*   start end    The range of entries to use [defaults to full range]
*
* Options:
*     [PRINT]/NOPRINT
*     TITLE=title of report ["Variance Ratio Test, Series xxx"]
*
*     LAGS=number of lags examined [2]
*     ROBUST/[NOROBUST]
*        If ROBUST, computes heteroscedasticity robust standard errors
*     [CORRECT]/NOCORRECT
*        If CORRECT, corrects for bias in the variance.
*     [CENTER]/NOCENTER
*        If CENTER, subtracts a fixed mean from the differenced data.
*     DIFF/[NODIFF]
*        If DIFF, difference the data, which means that the null is that
*        the input series is a (possibly drifting) random walk rather than
*        a series of independent increments.

* Variables Defined:
*   %%VRATIO         Variance Ratio Statistic
*   %CDSTAT          Test statistic
*   %SIGNIF          Asymptotic p-value
*
* Revision Schedule:
*   03/2005 Substantially rewritten by Tom Doan of Estima from original procedure
*      submitted by Christopher J. Zorn, Department of Political Science,
*      Ohio State University. Added small sample corrections to variance calculations,
*      report of asymptotic p-values added, NOPRINT and ROBUST options
*   05/2006 Document DIFFS option. Correct coding for NODIFFS branch. Add output of
*      VR statistic (not just test statistic). Fix integer overflow for very large
*      number of observations in the CORRECTED case.
*   09/2019 Switch to use of REPORT instead of DISPLAY
*
procedure vratio series start end
type series series
type integer start end
*
option integer lags    2
option switch  print   1
option string  title
option switch  robust  0
option switch  correct 1
option switch  center  1
option switch  diff    0
*
local integer  startp endp nobs i
local series   xdiff xseries devsqr
local real     meandif varone varlag
local real     ssqlag ssqone
local vector   lwform
local report   vreport
local string   ltitle
*
declare real %%vratio
*
if .not.%defined(series) {
   display "##Syntax: @VRATIO(options) <<series>> start end"
   return
}
*
* xseries is the series which should behave as a RW
* under the null. xdiff is its increments
*
inquire(series=series)  startp<<start endp<<end
if diff {
   set xdiff startp+1 endp = series - series{1}
   set xseries startp endp = series
}
else {
   set xdiff startp endp = series
   acc series startp endp xseries
}
*
statistics(noprint) xdiff startp+1 endp
compute nobs = %nobs
*
* Check for missing values
*
if %nmiss {
   display "##@VRATIO cannot be applied because the series contains missing values."
   return
}

if center
   compute meandif=%mean,varone=%variance,ssqone=varone*(nobs-1)
else {
   compute meandif=0.0
   sstats startp+1 endp xdiff^2>>ssqone
   compute varone=ssqone/nobs
}

sstats startp+lags endp (xseries-xseries{lags}-(lags*meandif))^2>>ssqlag

*
* This does the bias correction for the variance as in Campbell, Lo and
* MacKinlay, 2.4.36 with nobs=nq (in their notation).
*

if center.and.correct
   compute varlag  = ssqlag*nobs/(lags*(nobs-lags+1.0)*(nobs-lags))
else
   compute varone  = ssqone/nobs , varlag = ssqlag/(nobs*lags)

compute %%vratio = varlag / (lags * varone)
*
* If ROBUST, this does the heteroscedasticity robust asymptotic
* distribution This uses MCOV with a non-standard window applied to the
* square of the deviation between xdiff and its sample mean.
*
if robust {
   set devsqr startp+1 endp = (xdiff-meandif)^2
   dim lwform(2*lags-1)
   ewise lwform(i)=%if(i==lags,0.0,2.0*(%min(i,2.0*lags-i)/lags)^2)
   mcov(lwform=lwform,lags=lags-1)
   # devsqr
   compute %variance=%cmom(1,1)/ssqone^2
}
else
   compute %variance=2.0*(2*lags-1)*(lags-1)/(nobs*(3.0*lags))
compute %cdstat = (%%vratio-1)/sqrt(%variance)
compute %signif = %ztest(%cdstat)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Variance-Ratio Test on series "+%l(series)
*
report(use=vreport,action=define,title=ltitle)
report(use=vreport,atrow=1,atcol=1,span) ltitle
report(use=vreport,atrow=2,atcol=1,span) "From "+%datelabel(startp)+" to "+%datelabel(endp)
if robust
   report(use=vreport,atrow=3,atcol=1,span) "Computed with heteroscedasticity robust standard errors"
report(use=vreport,atrow=5,atcol=1,span) "Null is Variance Ratio=1. Two-Tailed Test."
report(use=vreport,atrow=6,atcol=1) "VR("+lags+")"       %%vratio
report(use=vreport,atrow=7,atcol=1) "Test Statistic"     %cdstat
report(use=vreport,atrow=8,atcol=1) "Asymptotic P-value" %signif
if print
   report(use=vreport,action=show)
end
