capt program drop lomackinlay
*! version 1.0.4  23may2006    C F Baum
*  from kpss 1.2.1, ipshin 1.0.8,
*  Lo, MacKinlay, A Non-Random Walk Down Wall Street (1999);
*  Campbell, Lo, MacKinlay, Econometrics of Financial Markets (1997); 
*  Lo and MacKinlay, "Stock market prices do not follow random walks: evidence
*  from a simple specification test", Rev. Fin. Studies 1:1, 1988, 41-66; and
* "A small-sample overlapping variance-ratio test", Tse, Ng and Zhang, 2002
*  www.mysmu.edu/faculty/yktse/JTSA_R.pdf 
*
*  1.0.2: output formatting, handle numlist of q values
*  1.0.3: deal with gaps, add gaps option
*  1.0.4: rewrite using LMK notation

* Allin Cottrell, Wake Forest University, 2006/08/02: fix several
* serious problems with version 1.0.4.  The modified program could 
* be formulated more elegantly.

program define lomackinlay, rclass
	version 9.2

	syntax varname(ts) [if] [in] [ , Q(numlist int >0) Gaps Robust]  

   	marksample touse
			/* get time variables */
    _ts tvar panelvar `if' `in', sort onepanel
	markout `touse' `timevar'
	tsreport if `touse' //, report
	if r(N_gaps) & "`gaps'"=="" {
		di in red "sample may not contain gaps"
		exit
	}

* require data to be tsset
	qui tsset
	local rts = r(tsfmt)
	local rtvar = r(timevar)
* if gaps option used, create new tsset variable
	if "`gaps'" ~= "" {
		tempvar timevar
		qui {
			gen `timevar' = _n  if `touse'
			tsset `timevar'
			local rts = r(tsfmt)
		}
	}
 			
	su `r(timevar)' if `touse', meanonly
	local n1 = r(min)
	local n2 = r(max)

        // Cottrell: add variable "cd1x2" below
	tempvar cd1x cd1x2 dx cdx dxl cdxl dprod
	tempname enn enn0 n vara varb mu qmu sumcdiff mult m v rs p vr 
	// Cottrell: add macro "deltaden" below
	tempname deltanum deltaden vhatr type
	qui {
		gen double `dx' = .
		gen double `cdx' = .
	}	
* get q: span of differencing
	local nq: word count `q'
	local sq: list sort q
	if !`nq' {
		local sq 2 4 8 16
	}
* ensure that all calcs done on smallest indicated sample
	local nsq: word count `sq'
	local maxlag: word `nsq' of `sq'
	markout `touse' L(1/`maxlag').`varlist'
	
* output header
	di " "
	di as text "Lo-MacKinlay modified overlapping Variance Ratio statistic for `varlist'"
	di as text "["`rts' `n1' " - " `rts' `n2' " ]" _n
	di as text "q {col 10} N {col 20} VR {col 32} R_s {col 42} p>|z|" _n "{hline 50}"
	
* generate centered first-differenced series
	qui {
	        // Cottrell: add "cd1x2" to drop list
		capt drop `cd1x' `cd1x2'
		gen double `cd1x' = D.`varlist' if `touse'
		sum `cd1x', meanonly
                // Cottrell: save mu, the mean of the first differences		
		scalar `mu' = r(mean)
		replace `cd1x' = `cd1x' - `mu'
* and its variance (LMK 3.2.7a)
		sum `cd1x'
		scalar `vara' = r(Var)
		// Cottrell: save cd1x2 and deltaden, next 3 lines --
                // these will be needed in calculating the robust S.E.
		gen double `cd1x2' = `cd1x'^2
		sum `cd1x2'
		scalar `deltaden' = r(sum) 
	}

foreach q of local sq {
* calculate n = int(N/q)
	sum `varlist' if `touse', meanonly
* total N available before differencing
	scalar `enn0' = r(N)
	scalar `n' = int(`enn0'/`q')

* generate q-differenced series
	qui	replace `dx' = `varlist' - L`q'.`varlist' if `touse'

* N, mean of q-differenced series
	sum `dx', meanonly
	scalar `enn' = r(N)
	// Cottrell: the next line is bad, mu shouldn't be redefined
	// scalar `mu' = r(mean)
	scalar `qmu' = `mu'*`q'
* center q-differenced series and square
	qui replace `cdx' = (`dx' - `qmu')^2
* sum of squares of q-differenced series 
	qui sum `cdx'
	scalar `sumcdiff' = r(sum) 
// di as err `sumcdiff'
* calculate multiplier with d.f. adjustment (3.2.7b)
	scalar `mult' = 1/(`q'*(`n'*`q'-`q'+1)*(1-`q'/(`n'*`q')))
// di as err `mult'
	scalar `varb' = `mult'*`sumcdiff'
* calculate VR_q and MR_q
	scalar `vr' = `varb'/`vara' 
	scalar `m' = `vr' -1
// di as err `m'
* and its variance V (3.2.6b)
	scalar `v' = (2*(2*`q'-1)*(`q'-1))/(3*`q')
	if "`robust'" == "" {
* standardized R statistic for homoskedastic case
                // Cottrell: fix omission of factor of sqrt(nq) in the next line
                // scalar `rs' = `m'/sqrt(`v')
		scalar `rs' = sqrt(`n'*`q')*`m'/sqrt(`v')
		local type homoskedastic
}
	else {
* for robust statistic, calculate delta_j sequence (3.2.9c)
		qui {
		capt drop `dxl' `cdxl' `dprod'
		gen double `dxl' = .
		gen double `cdxl' = .
		gen double `dprod' = .
		}
* set trace on
		scalar `vhatr' = 0
		local q1 = `q'-1
		forv j = 1/`q1' {
			local j1 = `j'-1
			qui	{
                                // Cottrell: scrap the following line
                                // replace `dxl' = L`j'.`varlist' - L`j1'.`varlist' if `touse'
                                // Cottrell: second term in numerator of delta is just the
                                // jth lag of the first term, which is just the squared, centered
                                // first differences of the original series
                                replace `cdxl' = L`j'.`cd1x2' if `touse'
				// Cottrell: scrap 2 wrong lines (`qmu' doesn't belong here)
				// replace `cdxl' =  (`dxl' - `qmu')^2
				// replace `dprod' = `cdx'*`cdxl'
				replace `dprod' = `cd1x2'*`cdxl'
			}
			sum `dprod', meanonly
			// Cottrell: remove extraneous factor of nq in next line: you
                        // need this only if you multiply m by sqrt(nq), but multiplying 
                        // on both sides is clearly redundant.
			// scalar `deltanum' = `n'*`q'*r(sum)
			scalar `deltanum' = r(sum)
			// Cottrell: use correct denominator element `deltaden' below
			// local delta`j' = `deltanum' / (`sumcdiff'^2) 
			local delta`j' = `deltanum' / (`deltaden'^2) 
* calculate vhat_q for robust case (3.2.9b)
			scalar `vhatr' = `vhatr' + (2*(`q'-`j')/`q')^2 * `delta`j''
		}
* standardized R statistic for heteroskedastic case
	scalar `rs' = `m'/sqrt(`vhatr')
	local type heteroskedastic
}
* calculate normal tail prob
	scalar `p' = 2 * normprob(-abs(`rs'))
* display statistic
// di in r `vara' " " `varb' _n
	di "`q' {col 10}" `enn' "{col 20}" %5.3f `vr' "{col 30}" %8.4f `rs' "{col 42}" %5.4f `p' 
	
* set up returns
	return scalar N_`q' = `enn'
	return scalar q_`q' = `q'
	return scalar R_`q' = `rs'
	return scalar p_`q' = `p'
	
} // end q loop
	if "`robust'" != "" {
			di _n "Test statistics robust to heteroskedasticity" _n
	}
	return local type = "`type'"
	
* restore orig tsset if needed
	if "`gaps'" ~= "" {
		qui tsset `rtvar'
	}
	
	end