' example subroutine for Lo-MacKinlay variance ratio testing (v1.0: 2009-04-30 - gts)
'
' *** note: this routine has undergone some testing and matches a variety
' of known results but is strictly use-at-your-own-risk ***
'
subroutine local vrtest(group g1, scalar lags, matrix vars, matrix vr, matrix ses, matrix probs)

' arguments
'	g1 - group containing series on which to perform the test
'	lags - highest lag to perform test (will compute from 2 to number of lags)
'	vars, vr, ses, probs - matrices to contain the output as described below
'
' output matrices  (columns are for series, rows are lags from 2 to number of lags):
'	vars - estimates of the centered variances (allow for drift)
'	vr - variance ratios
'	ses - asymptotic normal ses (assuming homoskedasticity)
'	probs - p-values for the vrtest
'	
' notes:
' 1. innovations are assumed to be differences of the original series
' 2. estimates variances assuming a constant (allows for drift)
' 3. variances estimated with no bias correction
' 4. ses and probs assume homoskedasticity
' 5. results do not account for joint testing
' 6. output matrices are not initialized to NAs
'
' example of code to call the routine on data in an EViews workfile
' --------------------------------------------------------------------------------------
' include vrtest.sub
' group group01 x y z
' matrix vars
' matrix vr
' matrix ses
' matrix probs
' call vrtest(group01, 3, vars, vr, ses, probs)
' --------------------------------------------------------------------------------------

' get the sizes
!n = g1.count
!l = lags

' redefine output matrix sizes
'  (rows are lags from 2 to number of lags, columns are series)

matrix(!l, !n) vars
matrix(!l-1, !n) vr
matrix(!l-1, !n) ses
matrix(!l-1, !n) probs

' loop over series
for !i=1 to !n
	' get the first difference of the series and the mean of the difference
	frml diff!i = d(g1(!i))
	scalar mu = @mean(diff!i)
	scalar obs = @obs(diff!i)

	' form the centered difference
	frml cdiff!i = diff!i - mu

	' get the base variance estimate
	scalar cvardiff = @sumsq(cdiff!i) / obs
	vars(1, !i) = cvardiff

	' loop over lags
	for !j=2 to !l
		' form the centered j-th difference
		frml ldiff!j = g1(!i) - @lag(g1(!i), !j) - !j*mu

		' get the variance estimate
		scalar vardiff!j = @sumsq(ldiff!j) / (!j*obs)
		vars(!j, !i) = vardiff!j

		' form the ratio
		vr(!j-1, !i) = vardiff!j/cvardiff

		' get the ses
		ses(!j-1, !i) = @sqrt(2*(2*!j - 1)*(!j-1)/(3*!j*obs))

		' get the z-stat prob
		scalar zstat!i!j = @abs(vr(!j-1, !i) - 1)/ses(!j-1, !i)
		probs(!j-1, !i) = 2*(1-@cnorm(zstat!i!j))

	next

next

endsub