*! actest9 2.0.13 22Jan2015
*! based on actest 2.0.12 and livreg2 Mata library 1.1.07
*! see end of file for version comments

program define actest9, rclass sortpreserve
version 9.2
	syntax [varlist(default=none ts)] [if] [in] [aw fw pw iw]					///
		[, 	LAGs(numlist integer>0 sort) STRICTexog q0 bp small					///
			Robust CLuster(varlist) BW(string) kernel(string) psd(string) sw ]

************* START SETUP BLOCKS FOR UNIVARIATE AND POST-REGRESSION *************

	local univariate=("`varlist'"~="")

	if `univariate' {			// Start setup block for univariate test

* univariate implies strict exogeneity and (in effect) OLS
		local strictexog "strictexog"
		local ols=1

		local varcount :  word count `varlist'
		if `varcount' > 1 {
di as err "error - can specify at most one variable as argument to actest"
		exit 198
		}

		marksample touse

* Weight support
		tempvar wvar
		if "`weight'" ~= "" {
* exp has "= " as well as weight expression
			qui gen double `wvar' `exp' if `touse'
		}
		else {
			local weight ""
			qui gen byte `wvar'=1 if `touse'
		}

* For the univariate case, the "residuals" are just the demeaned variable
		tempvar uhat
		sum `varlist' if `touse' [`weight'`exp'], meanonly
		qui gen double `uhat'=`varlist'-r(mean) if `touse'
* Local varlists x and z are just a vector of 1s
		tempvar one
		qui gen byte `one'=1 if `touse'
		local x "`one'"
		local z "`one'"
* Local cons=1 (since it's the only "regressor"
		local cons=1

* dofminus used by ivreg2, needed by m_omega()
		local dofminus=0


	}							// End setup block for univariate test
	else {						// Start setup block for post-regression test

		if "`e(cmd)'" == "" {
			error 301			// error 301 = "last estimates not found"
		}

* Partialling-out not (yet) supported
		if "`e(fwlcons)'" != ""  | (e(partial_ct)>0 & e(partial_ct)<.) {
di as err "actest not allowed after ivreg2 with partialling-out option"
			error 499
		}

		tempvar touse
		gen byte `touse'=e(sample)

* Is previous estimation command supported?
		if									///
			"`e(cmd)'" ~= "ivreg2" &		///
			"`e(cmd)'" ~= "ivreg28" &		///
			"`e(cmd)'" ~= "ivreg29" &		///
			"`e(cmd)'" ~= "ivreg210" &		///
			"`e(cmd)'" ~= "ivreg" &			///
			"`e(cmd)'" ~= "regress" &		///
			"`e(cmd)'" ~= "ivregress" &		///
			"`e(cmd)'" ~= "newey" &			///
			"`e(cmd)'" ~= "newey2" {
di as err "`e(cmd)' not supported by actest"
			exit 198
		}
		local cmd `e(cmd)'

* Weights support
		if "`weight'" ~= "" {
di as err "[`weight'=`exp'] not allowed - cannot change weights used by original estimation"
			exit 198
		}
		local weight "`e(wtype)'"
		local exp    "`e(wexp)'"
		tempvar wvar
		if "`weight'" ~= "" {
* exp has "= " as well as weight expression
			qui gen double `wvar' `exp' if `touse'
		}
		else {
			qui gen byte `wvar'=1 if `touse'
		}

* define the x matrix, and check if constant appears
		local x : colnames(e(b))
		local x : subinstr local x "_cons" "", count(local cons)

* define the z matrix
		local z "`e(insts)'"
		local ols=("`z'"=="")
		if `ols' {
			local z "`x'"
		}

* add constant to both x and z if present in reglist
		if `cons'==1 {
			tempvar one
			qui gen byte `one'=1 if `touse'
			local x "`x' `one'"
			local z "`z' `one'"
		}
		else {
			local noconstant "noconstant"
		}
/* cfb2.0.10
* Replace time-series ops with temp vars
		tsrevar `x'
		local x "`r(varlist)'"
		tsrevar `z'
		local z "`r(varlist)'"
*/

* Generate residuals
		tempvar uhat
		qui predict double `uhat' if `touse', resid

* dofminus used by ivreg2, needed by m_omega()
		if "`e(dofminus)'"=="" {
			local dofminus=0
		}
		else {
			local dofminus=`e(dofminus)'
		}
	}						// End setup block for post-regression test

//////// Start setup block shared by univariate and post-regression tests //////////

* Incompatible options
	if "`strictexog'"~="" & "`cluster'"~="" {
		di as err "Incompatible options: strictexog and cluster-robust"
		exit 198
	}
	if "`bp'"~="" & "`robust'`cluster'"~="" {
		di as err "Incompatible options: bp not allowed with robust/cluster-robust"
		exit 198
	}
* Other incompatibilities...?

* Box-Pierce (diagonal V_r) implies q0 (no s.c. under H0)
	if "`bp'"~="" {
		local q0 "q0"
	}

* Check tsset and create time and panel vars
	capture tsset
	local tvar "`r(timevar)'"
* If no time var, exit with error
	if "`tvar'"=="" {
		di as err "Data must be tsset"
		error 198
	}
	if "`r(tdelta)'" != "" {
		local tdelta = `r(tdelta)'
	}
	else {
		local tdelta=1
	}
* If not panel data, ivar will be empty
	local ivar "`r(panelvar)'"

* Default lags
	if "`lags'"=="" {
		local lags=1
	}
* Process lags
	local numlags : list sizeof lags
* tokenize and then take first and last element from lags
* If only one lag supplied, then 1st lag=1 and last lag=supplied
* If multiple lags supplied use lowest and highest
	if `numlags'==1 {
		local llo = 1
		local lhi = `lags'
	}
	else {
		tokenize `lags'
		local llo = `1'
		local lhi = ``numlags''
	}

* If bw specified, must provide kernel as well.
* If kernel specified, must provide bw as well.
	if	("`kernel'"~="" & "`bw'"=="") |	///
		("`kernel'"=="" & "`bw'"~="") {
		di as err "error: must supply both bandwidth and kernel with kernel-robust option"
		exit 198
	}

* Check valid bw and kernel and replace with unabbreviated kernel name
* and bw converted to scalar
	if "`kernel'" ~= "" {
* bw argument is string
* s_vkernel is in livreg2 mlib.
		mata: s_actest9_vkernel("`kernel'", "`bw'", "`ivar'")
		local kernel `r(kernel)'
		local bw = `r(bw)'
	}
* Otherwise convert bw string "0" to scalar 0
	else if "`bw'"=="0" | "`bw'"=="" {
		local bw = 0
	}
* Check bandwidth is legal
	sum `tvar' if `touse', meanonly
	local T = r(max)-r(min)+1
	if `bw' > (`T'-1)/`tdelta' {
di as err "invalid bandwidth in option bw() - cannot exceed timespan of data"
		exit 198
	}

* Default behavior: truncated kernel with bw=q=(min lag - 1).
* When min lag = 1, q=0 and no serial correlation assumed.
* Indicated by kernel="default".
* Overridden if kernel or cluster-robust specified
	if "`kernel'`cluster'"=="" {
		local kernel "default"
	}
* Cluster => robust
* Do here, after options are processed, since specifying just "robust" as option means no AC or HAC.
	if "`cluster'"~="" {
		local robust "robust"
	}

* Every time a weight is used, must multiply by scalar wf ("weight factor")
* wf=1 for no weights, fw and iw, wf = scalar that normalizes sum to be N if aw or pw
	sum `wvar' if `touse' [`weight'`exp'], meanonly
	if "`weight'"=="" | "`weight'"=="fweight" | "`weight'"=="iweight" {
* Effective number of observations is sum of weight variable.
* If weight is "", weight var must be column of ones and N is number of rows
		local wf=1
		local N=r(sum_w)
	}
	else {
		local wf=r(N)/r(sum_w)
		local N=r(N)
	}

// Replace TS operators with temp vars; tweak lists to deal with omitted regressors,
// as tsrevar does not cope with o. prefix
/*
	if `c(stata_version)'>=11 {
		foreach var in `x' `z' {
			_ms_parse_parts `var'

	di in red "`var'"
	ret li

			if `r(omit)' {
				local remove "`var'"
				local x : list x - remove
				local z : list z - remove
			}
		}
	}
*/
// tsrevar -> fvrevar in V11+ so that factor variables won't choke
	loc fix ts
	if `c(stata_version)'>=11 {
		loc fix fv
	}
	`fix'revar `x'
	local x "`r(varlist)'"
	`fix'revar `z'
	local z "`r(varlist)'"

* tindex variable required by m_omega()
	tempvar tindex
	qui gen `tindex'=1 if `touse'
	qui replace `tindex'=sum(`tindex') if `touse'


************* END SETUP BLOCK *************

* Calculate ACs up to lhi for specified sample
	tempname ac regest
* If starting with a lag>1, fill out list of ACs with zeros
	if `llo'>1 {
		mat `ac' = J(1,`llo'-1,0)
	}
	forvalues lag = `llo'/`lhi' {
		tempvar luhat
* Generate required lags of residual vector.
* Lag=0 if uhat available but lag is not because it's a leading ob (e.g., before data start)
		qui gen double `luhat'=L`lag'.`uhat' if `touse'
		qui replace `luhat'=0 if `luhat'==. & `uhat'<. & `touse'
		local lagresids "`lagresids' `luhat'"
* capture in case this is a univariate version and there is nothing to hold/restore
		capture _estimates hold `regest', restore
		qui reg `luhat' `uhat' if `touse' [`weight'`exp'], `noconstant'
		mat `ac' = nullmat(`ac'), _b[`uhat']
		capture _estimates unhold `regest'
	}

* Code lifted from ivreg2. Last block of code before test stats.
********** CLUSTER SETUP **********************************************

* Mata code requires data are sorted on (1) the first var cluster if there
* is only one cluster var; (2) on the 3rd and then 1st if two-way clustering,
* unless (3) two-way clustering is combined with kernel option, in which case
* the data are tsset and sorted on panel id (first cluster variable) and time
* id (second cluster variable).
* Second cluster var is optional and requires an identifier numbered 1..N_clust2,
* unless combined with kernel option, in which case it's the time variable.
* Third cluster var is the intersection of 1 and 2, unless combined with kernel
* opt, in which case it's unnecessary.
* Sorting on "cluster3 cluster1" means that in Mata, panelsetup works for
* both, since cluster1 nests cluster3.
* Note that it is possible to cluster on time but not panel, in which case
* cluster1 is time, cluster2 is empty and data are sorted on panel-time.
* Note also that if data are sorted here but happen to be tsset, will need
* to be re-tsset after estimation code concludes.

		if "`cluster'"!="" {
			local clopt "cluster(`cluster')"
			tokenize `cluster'
			local cluster1 "`1'"
			local cluster2 "`2'"
			if "`kernel'"~="" {
* kernel requires either that cluster1 is time var and cluster2 is empty
* or that cluster1 is panel var and cluster2 is time var.
* Either way, data must be tsset and sorted for panel data.
				if "`cluster2'"~="" {
* Allow backwards order
					if "`cluster1'"=="`tvar'" & "`cluster2'"=="`ivar'" {
						local cluster1 "`2'"
						local cluster2 "`1'"
					}
					if "`cluster1'"~="`ivar'" | "`cluster2'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset panel & time vars."
di as err "       tsset panel var=`ivar'; tsset time var=`tvar'; cluster vars=`cluster1',`cluster2'"
						exit 198
					}
				}
				else {
					if "`cluster1'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset time variable."
di as err "       tsset time var=`tvar'; cluster var=`cluster1'"
						exit 198
					}
				}
			}
* Simple way to get quick count of 1st cluster variable without disrupting sort
* clusterid1 is numbered 1.._Nclust1.
			tempvar clusterid1
			qui egen `clusterid1'=group(`cluster1') if `touse'
			sum `clusterid1' if `touse', meanonly
			if "`cluster2'"=="" {
				local N_clust=r(max)
				local N_clust1=.
				local N_clust2=.
				if "`kernel'"=="" {
* Single level of clustering and no kernel-robust, so sort on single cluster var.
* kernel-robust already sorted via tsset.
					sort `cluster1'
				}
			}
			else {
				local N_clust1=r(max)
				if "`kernel'"=="" {
					tempvar clusterid2 clusterid3
* New cluster id vars are numbered 1..N_clust2 and 1..N_clust3
					qui egen `clusterid2'=group(`cluster2') if `touse'
					qui egen `clusterid3'=group(`cluster1' `cluster2') if `touse'
* Two levels of clustering and no kernel-robust, so sort on cluster3/nested in/cluster1
* kernel-robust already sorted via tsset.
					sort `clusterid3' `cluster1'
					sum `clusterid2' if `touse', meanonly
					local N_clust2=r(max)
				}
				else {
* Need to create this only to count the number of clusters
					tempvar clusterid2
					qui egen `clusterid2'=group(`cluster2') if `touse'
					sum `clusterid2' if `touse', meanonly
					local N_clust2=r(max)
* Now replace with original variable
					local clusterid2 `cluster2'
				}
				local N_clust=min(`N_clust1',`N_clust2')
			}
		}
		else {
* No cluster options but for Mata purposes, set N_clust=0
				local N_clust=0
		}

************************************************************************************************

* Test description
	di _n as text "Cumby-Huizinga test for autocorrelation" _c
	if "`ivar'"=="" {											//  Time-series large-T tests
		if `univariate' & `ols' & "`bp'"~="" & "`small'"=="" {
			di as text " (Box-Pierce)"
		}
		else if `univariate' & `ols' & "`bp'"~="" & "`small'"~="" {
			di as text " (Ljung-Box)"
		}
		else if `ols' & "`bp'"~="" & "`small'"=="" {
			di as text " (Modified Box-Pierce)"
		}
		else if `ols' & "`bp'"~="" & "`small'"~="" {
			di as text " (Modified Ljung-Box)"
		}
		else if `ols' & "`strictexog'`robust'`cluster'"=="" & `bw'==0 {
			di as text " (Breusch-Godfrey)"
		}
		else {
			di
		}
	}
	else {														//  Panel large-N tests
		if "`cluster'"~="" {
			di as text " (Arellano-Bond)"
		}
		else {
			di
		}
	}

	if `univariate' {
		local testobj "disturbance"
	}
	else {
		local testobj "variable"
	}

	di as text "  H0: `testobj' is MA process up to order q"
	di as text "  HA: serial correlation present at specified lags >q"

	di as text "{hline 41}{c TT}{hline 35}"
	if "`q0'"~="" {
		di	as text "  H0: q=0 (serially uncorrelated)"			///
			as text _col(42) "{c |}"							///
			as text "  H0: q=0 (serially uncorrelated)"
		di	as text "  HA: s.c. present at range specified"		///
			as text _col(42) "{c |}"							///
			as text "  HA: s.c. present at lag specified"
	}
	else if `llo'==1 {
		di	as text "  H0: q=0 (serially uncorrelated)"			///
			as text _col(42) "{c |}"							///
			as text "  H0: q=specified lag-1"
		di	as text "  HA: s.c. present at range specified"		///
			as text _col(42) "{c |}"							///
			as text "  HA: s.c. present at lag specified"
	}
	else {
		di	as text "  H0: `testobj' is MA(q), q=" `llo'-1		///
			as text _col(42) "{c |}"							///
			as text "  H0: q=specified lag-1"
		di	as text "  HA: s.c. present at range specified"		///
			as text _col(42) "{c |}"							///
			as text "  HA: s.c. present at lag specified"
	}

	di	as text "{hline 11}{c TT}{hline 29}{c +}{hline 5}{c TT}{hline 29}"
	di	as text _col(5) "lags"													///
		as text _col(12) "{c |}"												///
		as text _col(19) "chi2"													///
		as text _col(29) "df"													///
		as text _col(36) "p-val"												///
		as text _col(42) "{c |}"												///
		as text _col(44) "lag"													///
		as text _col(48) "{c |}"												///
		as text _col(55) "chi2"													///
		as text _col(65) "df"													///
		as text _col(72) "p-val"
	di	as text "{hline 11}{c +}{hline 29}{c +}{hline 5}{c +}{hline 29}"

* C-H notation:
* q=last lag BEFORE lags to be tested
* NB: q=0 if testing lags starting with lag 1
* Start testing at lag q+1
* s=number of ACs to be tested
* q+s=last lag to be tested

* Loop through lags
* local macro lag is highest lag being tested in that iteration
* Matrix m stores results
* Matrix r in loop is one row (lag set) of results
* Tokenize stored lagged resids; each time through loop,
* add to list (capuhat) and mac shift to next resid.
	tempname m vrhat vr
	tokenize `lagresids'
	forvalues lag = `llo'/`lhi' {
		tempname lqs lqsp r
		local uq "`1'"
		mac shift
		local capuhat "`capuhat' `uq'"

* Call to test set of lags starting at q+1
		local q=`llo'-1
		local s=`lag'-`q'
		if ("`kernel'"=="default" & ((`q'==0) | "`q0'"~="")) {
* q=0, serially uncorrelated under null => not AC or HAC
			local loopkernel ""
			local loopbw=0
		}
		else if ("`kernel'"=="default" & `q'>0) {
* q>0, serially correlated up to q => HAC, truncated kernel w/bw=q
			local loopkernel "Truncated"
			local loopbw=`q'
		}
		else {
* user-defined, including cluster-robust
			local loopkernel "`kernel'"
			local loopbw=`bw'
		}
		mata: s_acstat(						///
						"`x'",				///
						"`z'",				///
						`q',				///
						`s',				///
						"`uhat'",			///
						"`capuhat'",		///
						"`ac'",				///
						"`strictexog'",		///
						"`bp'",				///
						"`small'",			///
						"`touse'",			///
						"`weight'",			///
						"`wvar'",			///
						`wf',				///
						"`tvar'",			///
						`tdelta',			///
						"`ivar'",			///
						"`tindex'",			///
						"`robust'",			///
						"`clusterid1'",		///
						"`clusterid2'",		///
						"`clusterid3'",		///
						`loopbw',			///
						"`loopkernel'",		///
						"`psd'",			///
						`dofminus',			///
						"`sw'"				///
						)
		scalar `lqs'=r(lqs)
		scalar `lqsp'=chiprob(`s',`lqs')
		mat `r'=`llo',`lag',`lqs',`lqsp'
		mat `vrhat'=r(vrhat)
		mat `vr'=r(vr)

		di										///
			as result _col(2) %3.0f `llo'		///
			as text _col(6) "-"					///
			as result _col(7) %3.0f `lag'		///
			as text _col(12) "{c |}"			///
			as result _col(17) %7.3f `lqs' _c   ///
			as result _col(30) "`s'"
		if `r(npsd)' {
			di as text "*" _c
			local npsdflag = "npsdflag"
		}
		di	as result _col(35) %5.4f `lqsp' _c

* Call to test at specified lag only
		local q=`lag'-1
		local s=1
		if ("`kernel'"=="default" & ((`q'==0) | "`q0'"~="")) {
* q=0, serially uncorrelated under null => not AC or HAC
			local loopkernel ""
			local loopbw=0
		}
		else if ("`kernel'"=="default" & `q'>0) {
* q>0, serially correlated up to q => truncated kernel w/bw=q
			local loopkernel "Truncated"
			local loopbw=`q'
		}
		else {
* user-defined, including cluster-robust
			local loopkernel "`kernel'"
			local loopbw=`bw'
		}
		mata: s_acstat(						///
						"`x'",				///
						"`z'",				///
						`q',				///
						`s',				///
						"`uhat'",			///
						"`uq'",				///
						"`ac'",				///
						"`strictexog'",		///
						"`bp'",				///
						"`small'",			///
						"`touse'",			///
						"`weight'",			///
						"`wvar'",			///
						`wf',				///
						"`tvar'",			///
						`tdelta',			///
						"`ivar'",			///
						"`tindex'",			///
						"`robust'",			///
						"`clusterid1'",		///
						"`clusterid2'",		///
						"`clusterid3'",		///
						`loopbw',			///
						"`loopkernel'",		///
						"`psd'",			///
						`dofminus',			///
						"`sw'"				///
						)
		scalar `lqs'=r(lqs)
		scalar `lqsp'=chiprob(1,`lqs')
		mat `r'=`r',`lag',`lqs',`lqsp'

		di										///
			as text _col(42) "{c |}"			///
			as result _col(44) %3.0f `lag'		///
			as text _col(48) "{c |}"			///
			as result _col(53) %7.3f `lqs' _c   ///
			as result _col(66) "1"
		if `r(npsd)' {
			di as text "*" _c
			local npsdflag = "npsdflag"
		}
		di	as result _col(71) %5.4f `lqsp'

		matrix rownames `r' = lags_`llo'_`lag'
		mat `m'=nullmat(`m') \ `r'

	}

	noi di as text "{hline 11}{c BT}{hline 29}{c BT}{hline 5}{c BT}{hline 29}"
	matrix colnames `m' = lag_l lag_u chi2_lu p_lu lag chi2 p

* Footnotes
	if ~`univariate' {
		if "`strictexog'"=="" {
			di as text "  Test allows predetermined regressors/instruments"
		}
		else {
			di as text "  Test requires strictly exogenous regressors/instruments"
		}
	}
	if "`cluster'"~="" {
		di as text "  Test robust to heteroskedasticity and within-cluster autocorrelation"
	}
	else if "`robust'"~="" {
		di as text "  Test robust to heteroskedasticity"
	}
	else {
		di as text "  Test requires conditional homoskedasticity"
	}
	if "`npsdflag'" ~= "" {
		di as text "  * Eigenvalues adjusted to make matrix positive semidefinite"
	}

	tempname lqs lqsp
	local lastrow=rowsof(`m')
	return scalar df=(`m'[`lastrow',2]-`m'[`lastrow',1]+1)
	return scalar p=`m'[`lastrow',4]
	return scalar chi2=`m'[`lastrow',3]

	return scalar maxlag=`lhi'
	return scalar minlag=`llo'
	return local kernel = "`kernel'"
	return scalar bw=`bw'
	return matrix results=`m'

	return matrix vrhat=`vrhat'
	return matrix vr=`vr'

end

version 9.2
mata:

// ********* struct ms_actest9_vcvorthog ******************* //
struct ms_actest9_vcvorthog {
	string scalar	ename, Znames, touse, weight, wvarname
	string scalar	robust, clustvarname, clustvarname2, clustvarname3, kernel
	string scalar	sw, psd, ivarname, tvarname, tindexname
	real scalar		wf, N, bw, tdelta, dofminus
	real matrix		ZZ
	pointer matrix	e
	pointer matrix	Z
	pointer matrix	wvar
}

void s_acstat(string scalar xvars,
				string scalar zvars,
				real scalar q,
				real scalar s,
				string scalar uvars,
				string scalar capuvars,
				string scalar ac,
				string scalar strictexog,
				string scalar bp,
				string scalar small,
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				string scalar tvarname,
				scalar tdelta,
				string scalar ivarname,
				string scalar tindexname,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				scalar bw,
				string scalar kernel,
				string scalar psd,
				scalar dofminus,
				string scalar sw)
{

	xv=tokens(xvars)
	vx= xv[|1,.|]
	zv=tokens(zvars)
	vz= zv[|1,.|]
	uv=tokens(uvars)
	vu= uv[|1,.|]
	capuv=tokens(capuvars)
	vg= capuv[|1,.|]
	wv=tokens(wvarname)
	vw= wv[|1,.|]

	st_view(X=.,.,vx,touse)
	st_view(Z=.,.,vz,touse)
	st_view(uhat=.,.,vu,touse)
	st_view(capuhat=.,.,vg,touse)

	T=rows(uhat)

	h=cols(Z)

	st_view(eta,.,(vz,vg),touse)
	st_view(wvar,   ., vw, touse)
	ZZ = quadcross(eta, wf*wvar, eta)
	XZ = quadcross(X, wf*wvar, Z)

	struct ms_actest9_vcvorthog scalar vcvo
	vcvo.ename			= uv[|1,1|]
	vcvo.Znames			= (zv,capuv)
	vcvo.touse			= touse
	vcvo.weight			= weight
	vcvo.wvarname		= wvarname
	vcvo.robust			= robust
	vcvo.clustvarname	= clustvarname
	vcvo.clustvarname2	= clustvarname2
	vcvo.clustvarname3	= clustvarname3
	vcvo.kernel			= kernel
	vcvo.sw				= sw
	vcvo.psd			= psd
	vcvo.ivarname		= ivarname
	vcvo.tvarname		= tvarname
	vcvo.tindexname		= tindexname
	vcvo.wf				= wf
	vcvo.N				= T
	vcvo.bw				= bw
	vcvo.tdelta			= tdelta
	vcvo.dofminus		= dofminus
	vcvo.ZZ				= ZZ

	vcvo.e		= &uhat
	vcvo.Z		= &eta
	vcvo.wvar	= &wvar

	psi=m_actest9_omega(vcvo)

// psda option: Stock-Watson 2008 Econometrica, Remark 8, say replace neg EVs with abs(EVs).
	if (det(psi) < 0) {
// Use S-W approach to make PSD
		vcvo.psd = "psda"
		psi=m_actest9_omega(vcvo)
		st_numscalar("r(npsd)", 1)
	}
	else if (vcvo.psd ~= "" ) {
// User had already supplied either psda or psd0 option to make PSD
		st_numscalar("r(npsd)", 1)
	}
	else {
		st_numscalar("r(npsd)",0)
	}

// Construct estimate of V_r (eqn 13)

	allr = st_matrix(ac)
	qmax = q + s
// (9)
	r = allr[|(q+1)\qmax|]'
// (22)
	sigma2 = quadcross(uhat, wf*wvar, uhat) / T

// s x s submatrix of psi is vr * sigma^4
	vr = psi[|(h+1,h+1)\(h+s,h+s)|] / (sigma2^2)

// Force vr to be identity matrix (Box-Pierce)
	if (bp~="") {
		vr = I(rows(vr))
	}


// Prop 1 - variance
	if (strictexog=="strictexog") {
// If strictly exogenous, Prop 1 simplifies: V_rhat = V_r
		vrhat = vr
	}
	else {
// Not strictly exogenous: Prop 1 is: V_rhat = V_r + BV_dB' + CD'B' + BDC'

// (5) construct A and Ainv
		A = psi[|(1,1)\(h,h)|]
		Ainv = invsym(A)

// (24)
		dhat = T * invsym(XZ * Ainv * XZ') * XZ * Ainv

// (23) depends on capuhat
		bhat = - quadcross(capuhat, wf*wvar, X) / T / sigma2

		ul = bhat * dhat
		ur = J(s,s,0)
		ll = J(s,h,0)
		lr = I(s) / sigma2
		phi = (ul, ur \ ll, lr)

// (20)
		ppp = phi * psi * phi'

		bvb = ppp[|(1,1)\(s,s)|]
		bdc = ppp[|(1,s+1)\(s,2*s)|]
		cdb = bdc'
// Proposition 1, not strictly exogenous
		if (bp=="") {
// General case
			vrhat = vr + bvb + cdb + bdc
		}
		else {
// C-H case (iii) aka Hayashi "Modified Box-Pierce"
			vrhat = I(rows(vr)) - bvb
		}
	}

	vrhatinv = invsym(vrhat)

// Ljung-Box small-sample correction
	if (small~="") {
		for (i=1; i<=rows(r); i++) {
			r[i,1] = sqrt((T+2)/(T-(q+i))) * r[i,1]
		}
	}

// Proposition 2 test stat
	lqs = T * r' * vrhatinv * r
	st_numscalar("r(lqs)",lqs)
	st_matrix("r(vrhat)",vrhat)
	st_matrix("r(vr)",vr)

}

// *********************************************************************** //
// **************** SUPPORT CODE (prev in livreg2.mlib ******************* //
// *********************************************************************** //


// ************************* s_actest9_vkernel ***************************** //
// Program checks whether kernel and bw choices are valid.
// s_actest9_vkernel is called from Stata.
// Arguments are the kernel name (req), bandwidth (req) and ivar name (opt).
// All 3 are strings.
// Returns results in r() macros.
// r(kernel) - name of kernel (string)
// r(bw) - bandwidth (scalar)

void s_actest9_vkernel(	string scalar kernel,
						string scalar bwstring,
						string scalar ivar
				)
{

// Check bandwidth
	if (bwstring=="auto") {
		bw=-1
	}
	else {
		bw=strtoreal(bwstring)
		if (bw==.) {
			printf("{err}bandwidth option bw() required for HAC-robust estimation\n")
			exit(102)
		}
		if (bw<=0) {
			printf("{err}invalid bandwidth in option bw() - must be real > 0\n")
			exit(198)
		}
	}

// Check ivar
	if (bwstring=="auto" & ivar~="") {
			printf("{err}Automatic bandwidth selection not available for panel data\n")
			exit(198)
	}

// Check kernel
// Valid kernel list is abbrev, full name, whether special case if bw=1
// First in list is default kernel = Barlett
	vklist = 	(	("", "bartlett", "0")
				\	("bar", "bartlett", "0")
				\	("bartlett", "bartlett", "0")
				\	("par", "parzen", "0")
				\	("parzen", "parzen", "0")
				\	("tru", "truncated", "1")
				\	("truncated", "truncated", "1")
				\	("thann", "tukey-hanning", "0")
				\	("tukey-hanning", "tukey-hanning", "0")
				\	("thamm", "tukey-hamming", "0")
				\	("tukey-hamming", "tukey-hamming", "0")
				\	("qua", "quadratic spectral", "1")
				\	("qs", "quadratic spectral", "1")
				\	("quadratic-spectral", "quadratic spectral", "1")
				\	("quadratic spectral", "quadratic spectral", "1")
				\	("dan", "danielle", "1")
				\	("danielle", "danielle", "1")
				\	("ten", "tent", "1")
				\	("tent", "tent", "1")
			)
	kname=strltrim(strlower(kernel))
	pos = (vklist[.,1] :== kname)

// Exit with error if not in list
	if (sum(pos)==0) {
		printf("{err}invalid kernel\n")
		exit(198)
		}

	vkname=strproper(select(vklist[.,2],pos))
	st_global("r(kernel)", vkname)
	st_numscalar("r(bw)",bw)

// Warn if kernel is type where bw=1 means no lags are used
	if (bw==1 & select(vklist[.,3],pos)=="0") {
		printf("{result}Note: kernel=%s", vkname)
		printf("{result} and bw=1 implies zero lags used.  Standard errors and\n")
		printf("{result}      test statistics are not autocorrelation-consistent.\n")
	}
}  // end of program s_actest9_vkernel


// ************************ m_actest9_omega ************************************** //

real matrix m_actest9_omega(struct ms_actest9_vcvorthog scalar vcvo)
{
	if (vcvo.clustvarname~="") {
		st_view(clustvar, ., vcvo.clustvarname, vcvo.touse)
		info = panelsetup(clustvar, 1)
		N_clust=rows(info)
		if (vcvo.clustvarname2~="") {
			st_view(clustvar2, ., vcvo.clustvarname2, vcvo.touse)
			if (vcvo.kernel=="") {
				st_view(clustvar3, ., vcvo.clustvarname3, vcvo.touse) // needed only if not panel tsset
			}
		}
	}

	if (vcvo.kernel~="") {
		st_view(t,    ., st_tsrevar(vcvo.tvarname),  vcvo.touse)
		T=max(t)-min(t)+1
	}

	if ((vcvo.kernel=="Bartlett") | (vcvo.kernel=="Parzen") | (vcvo.kernel=="Truncated") ///
		 | (vcvo.kernel=="Tukey-Hanning")| (vcvo.kernel=="Tukey-Hamming")) {
		window="lag"
	}
	else if ((vcvo.kernel=="Quadratic Spectral") | (vcvo.kernel=="Danielle") | (vcvo.kernel=="Tent")) {
		window="spectral"
	}
	else if (vcvo.kernel~="") {
// Should never reach this point
printf("\n{error:Error: invalid kernel}\n")
		exit(error(3351))
	}

	L=cols(*vcvo.Z)
	K=cols(*vcvo.e)		// ivreg2 always calls with K=1; ranktest may call with K>=1.

// Covariance matrices
// shat * 1/N is same as estimated S matrix of orthog conditions

// Block for homoskedastic and AC.  dof correction if any incorporated into sigma estimates.
	if ((vcvo.robust=="") & (vcvo.clustvarname=="")) {
// ZZ is already calculated as an external
		ee = quadcross(*vcvo.e, vcvo.wf*(*vcvo.wvar), *vcvo.e)
		sigma2=ee/(vcvo.N-vcvo.dofminus)
		shat=sigma2#vcvo.ZZ
		if (vcvo.kernel~="") {
			if (window=="spectral") {
				TAU=T/vcvo.tdelta-1
			}
			else {
				TAU=vcvo.bw
			}
			tnow=st_data(., vcvo.tindexname)
			for (tau=1; tau<=TAU; tau++) {
				kw = m_actest9_calckw(tau, vcvo.bw, vcvo.kernel)
				if (kw~=0) {						// zero weight possible with some kernels
													// save an unnecessary loop if kw=0
													// remember, kw<0 possible with some kernels!
					lstau = "L"+strofreal(tau)
					tlag=st_data(., lstau+"."+vcvo.tindexname)
					tmatrix = tnow, tlag
					svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
					tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix
// if no lags exist, tmatrix has zero rows.
					if (rows(tmatrix)>0) {
// col 1 of tmatrix has row numbers of all rows of data with this time period that have a corresponding lag
// col 2 of tmatrix has row numbers of all rows of data with lag tau that have a corresponding ob this time period.
// Should never happen that fweights or iweights make it here,
// but if they did the next line would be sqrt(wvari)*sqrt(wvari1) [with no wf since not needed for fw or iw]
						wv = (*vcvo.wvar)[tmatrix[.,1]]		///
									:* (*vcvo.wvar)[tmatrix[.,2]]*(vcvo.wf^2)	// inner weighting matrix for quadcross
						sigmahat = quadcross((*vcvo.e)[tmatrix[.,1],.],   wv ,(*vcvo.e)[tmatrix[.,2],.])	///
									/ (vcvo.N-vcvo.dofminus)					// large dof correction
						ZZhat    = quadcross((*vcvo.Z)[tmatrix[.,1],.], wv, (*vcvo.Z)[tmatrix[.,2],.])
						ghat = sigmahat#ZZhat
						shat=shat+kw*(ghat+ghat')
					}
				}	// end non-zero kernel weight block
			}	// end tau loop
		}  // end kernel code
// Note large dof correction (if there is one) has already been incorporated
	shat=shat/vcvo.N
	}  // end homoskedastic, AC code

// Block for robust HC and HAC but not Stock-Watson and single clustering.
// Need to enter for double-clustering if one cluster is time.
	if ( (vcvo.robust~="") & (vcvo.sw=="") & ((vcvo.clustvarname=="")		///
			| ((vcvo.clustvarname2~="") & (vcvo.kernel~="")))  ) {
		if (K==1) {										// simple/fast where e is a column vector
			if ((vcvo.weight=="fweight") | (vcvo.weight=="iweight")) {
				wv = (*vcvo.e:^2) :* *vcvo.wvar
			}
			else {
				wv = (*vcvo.e :* *vcvo.wvar * vcvo.wf):^2		// wf needed for aweights and pweights
			}
			shat=quadcross(*vcvo.Z, wv, *vcvo.Z)		// basic Eicker-Huber-White-sandwich-robust vce
		}
		else {											// e is a matrix so must loop
			shat=J(L*K,L*K,0)
			for (i=1; i<=rows(*vcvo.e); i++) {
				eZi=((*vcvo.e)[i,.])#((*vcvo.Z)[i,.])
				if ((vcvo.weight=="fweight") | (vcvo.weight=="iweight")) {
// wvar is a column vector. wf not needed for fw and iw (=1 by dfn so redundant).
					shat=shat+quadcross(eZi,eZi)*((*vcvo.wvar)[i])
				}
				else {
					shat=shat+quadcross(eZi,eZi)*((*vcvo.wvar)[i] * vcvo.wf)^2	//  **** ADDED *vcvo.wf
				}
			}
		}
		if (vcvo.kernel~="") {
// Spectral windows require looping through all T-1 autocovariances
			if (window=="spectral") {
				TAU=T/vcvo.tdelta-1
			}
			else {
				TAU=vcvo.bw
			}
			tnow=st_data(., vcvo.tindexname)
			for (tau=1; tau<=TAU; tau++) {
				kw = m_actest9_calckw(tau, vcvo.bw, vcvo.kernel)
				if (kw~=0) {						// zero weight possible with some kernels
													// save an unnecessary loop if kw=0
													// remember, kw<0 possible with some kernels!
					lstau = "L"+strofreal(tau)
					tlag=st_data(., lstau+"."+vcvo.tindexname)
					tmatrix = tnow, tlag
					svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
					tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix

// col 1 of tmatrix has row numbers of all rows of data with this time period that have a corresponding lag
// col 2 of tmatrix has row numbers of all rows of data with lag tau that have a corresponding ob this time period.
// if no lags exist, tmatrix has zero rows
					if (rows(tmatrix)>0) {
						if (K==1) {										// simple/fast where e is a column vector
// wv is inner weighting matrix for quadcross
							wv   = (*vcvo.e)[tmatrix[.,1]] :* (*vcvo.e)[tmatrix[.,2]]		///
								:* (*vcvo.wvar)[tmatrix[.,1]] :* (*vcvo.wvar)[tmatrix[.,2]] * (vcvo.wf^2)
							ghat = quadcross((*vcvo.Z)[tmatrix[.,1],.], wv, (*vcvo.Z)[tmatrix[.,2],.])
						}
						else {										// e is a matrix so must loop
							ghat=J(L*K,L*K,0)
							for (i=1; i<=rows(tmatrix); i++) {
								wvari =(*vcvo.wvar)[tmatrix[i,1]]
								wvari1=(*vcvo.wvar)[tmatrix[i,2]]
								ei    =(*vcvo.e)[tmatrix[i,1],.]
								ei1   =(*vcvo.e)[tmatrix[i,2],.]
								Zi    =(*vcvo.Z)[tmatrix[i,1],.]
								Zi1   =(*vcvo.Z)[tmatrix[i,2],.]
								eZi =ei#Zi
								eZi1=ei1#Zi1
// Should never happen that fweights or iweights make it here, but if they did
// the next line would be ghat=ghat+eZi'*eZi1*sqrt(wvari)*sqrt(wvari1)
// [without *vcvo.wf since wf=1 for fw and iw]
								ghat=ghat+quadcross(eZi,eZi1)*wvari*wvari1 * (vcvo.wf^2)	// ADDED * (vcvo.wf^2)
							}
						}
						shat=shat+kw*(ghat+ghat')
					}	// end non-zero-obs accumulation block
				}	// end non-zero kernel weight block
			}	// end tau loop
		}  // end kernel code
// Incorporate large dof correction if there is one
	shat=shat/(vcvo.N-vcvo.dofminus)
	}  // end HC/HAC code

	if (vcvo.clustvarname~="") {
// Block for cluster-robust
// 2-level clustering: S = S(level 1) + S(level 2) - S(level 3 = intersection of levels 1 & 2)
// Prepare shat3 if 2-level clustering
		if (vcvo.clustvarname2~="") {
			if (vcvo.kernel~="") {	// second cluster variable is time
									// shat3 was already calculated above as shat
				shat3=shat*(vcvo.N-vcvo.dofminus)
			}
			else {					// calculate shat3
									// data were sorted on clustvar3-clustvar1 so
									// clustvar3 is nested in clustvar1 and Mata panel functions
									// work for both.
				info3 = panelsetup(clustvar3, 1)
				if (rows(info3)==rows(*vcvo.e)) {	// intersection of levels 1 & 2 are all single obs
													// so no need to loop through row by row
					if (K==1) {										// simple/fast where e is a column vector
						wv = (*vcvo.e :* *vcvo.wvar * vcvo.wf):^2
						shat3=quadcross(*vcvo.Z, wv, *vcvo.Z)		// basic Eicker-Huber-White-sandwich-robust vce
					}
					else {											// e is a matrix so must loop
						shat3=J(L*K,L*K,0)
						for (i=1; i<=rows(*vcvo.e); i++) {
							eZi=((*vcvo.e)[i,.])#((*vcvo.Z)[i,.])
							shat3=shat3+quadcross(eZi,eZi)*((*vcvo.wvar)[i] * vcvo.wf)^2	//  **** ADDED *vcvo.wf
							}
						}
				}
				else {								// intersection of levels 1 & 2 includes some groups of obs
					N_clust3=rows(info3)
					shat3=J(L*K,L*K,0)
					for (i=1; i<=N_clust3; i++) {
						esub=panelsubmatrix(*vcvo.e,i,info3)
						Zsub=panelsubmatrix(*vcvo.Z,i,info3)
						wsub=panelsubmatrix(*vcvo.wvar,i,info3)
						wv = esub :* wsub * vcvo.wf
						if (K==1) {							// simple/fast where e is a column vector
							eZ = quadcross(1, wv, Zsub)		// equivalent to colsum(wv :* Zsub)
						}
						else {
							eZ = J(1,L*K,0)
							for (j=1; j<=rows(esub); j++) {
								eZ=eZ+(esub[j,.]#Zsub[j,.])*wsub[j,.] * vcvo.wf	//  **** ADDED *vcvo.wf
							}
						}
						shat3=shat3+quadcross(eZ,eZ)
					}
				}
			}
		}

// 1st level of clustering, no kernel-robust
// Entered unless 1-level clustering and kernel-robust
		if (!((vcvo.kernel~="") & (vcvo.clustvarname2==""))) {
			shat=J(L*K,L*K,0)
			for (i=1; i<=N_clust; i++) {		// loop through clusters, adding Z'ee'Z
												// for indiv cluster in each loop
				esub=panelsubmatrix(*vcvo.e,i,info)
				Zsub=panelsubmatrix(*vcvo.Z,i,info)
				wsub=panelsubmatrix(*vcvo.wvar,i,info)
				if (K==1) {						// simple/fast if e is a column vector
					wv = esub :* wsub * vcvo.wf
					eZ = quadcross(1, wv, Zsub)		// equivalent to colsum(wv :* Zsub)
				}
				else {
					eZ=J(1,L*K,0)
					for (j=1; j<=rows(esub); j++) {
						eZ=eZ+(esub[j,.]#Zsub[j,.])*wsub[j,.]*vcvo.wf	//  **** ADDED *vcvo.wf
					}
				}
				shat=shat+quadcross(eZ,eZ)
			}	// end loop through clusters
		}

// 2-level clustering, no kernel-robust
		if ((vcvo.clustvarname2~="") & (vcvo.kernel=="")) {
			imax=max(clustvar2)					// clustvar2 is numbered 1..N_clust2
			shat2=J(L*K,L*K,0)
			for (i=1; i<=imax; i++) {			// loop through clusters, adding Z'ee'Z
												// for indiv cluster in each loop
				svar=(clustvar2:==i)			// mimics panelsubmatrix but doesn't require sorted data
				esub=select(*vcvo.e,svar)		// it is, however, noticably slower.
				Zsub=select(*vcvo.Z,svar)
				wsub=select(*vcvo.wvar,svar)
				if (K==1) {						// simple/fast if e is a column vector
					wv = esub :* wsub * vcvo.wf
					eZ = quadcross(1, wv, Zsub)		// equivalent to colsum(wv :* Zsub)
				}
				else {
					eZ=J(1,L*K,0)
					for (j=1; j<=rows(esub); j++) {
						eZ=eZ+(esub[j,.]#Zsub[j,.])*wsub[j,.]*vcvo.wf	//  **** ADDED *vcvo.wf
					}
				}
				shat2=shat2+quadcross(eZ,eZ)
			}
		}

// 1st level of cluster, kernel-robust OR
// 2-level clustering, kernel-robust and time is 2nd cluster variable
		if (vcvo.kernel~="") {
			shat2=J(L*K,L*K,0)
// First, standard cluster-robust, i.e., no lags.
			i=min(t)
			while (i<=max(t)) {  				// loop through all T clusters, adding Z'ee'Z
												// for indiv cluster in each loop
				eZ=J(1,L*K,0)
				svar=(t:==i)					// select obs with t=i
				if (colsum(svar)>0) {			// there are obs with t=i
					esub=select(*vcvo.e,svar)
					Zsub=select(*vcvo.Z,svar)
					wsub=select(*vcvo.wvar,svar)
					if (K==1) {						// simple/fast if e is a column vector
						wv = esub :* wsub * vcvo.wf
						eZ = quadcross(1, wv, Zsub)		// equivalent to colsum(wv :* Zsub)
					}
					else {
// MISSING LINE IS NEXT
						eZ=J(1,L*K,0)
						for (j=1; j<=rows(esub); j++) {
							eZ=eZ+(esub[j,.]#Zsub[j,.])*wsub[j,.]*vcvo.wf	//  **** ADDED *vcvo.wf
						}
					}
					shat2=shat2+quadcross(eZ,eZ)
				}
				i=i+vcvo.tdelta
			} // end i loop through all T clusters

// Spectral windows require looping through all T-1 autocovariances
			if (window=="spectral") {
				TAU=T/vcvo.tdelta-1
			}
			else {
				TAU=vcvo.bw
			}

			for (tau=1; tau<=TAU; tau++) {
				kw = m_actest9_calckw(tau, vcvo.bw, vcvo.kernel)	// zero weight possible with some kernels
															// save an unnecessary loop if kw=0
															// remember, kw<0 possible with some kernels!
				if (kw~=0) {
					i=min(t)+tau*vcvo.tdelta				// Loop through all possible ts (time clusters)
					while (i<=max(t)) {						// Start at earliest possible t
						svar=t[.,]:==i						// svar is current, svar1 is tau-th lag
						svar1=t[.,]:==(i-tau*vcvo.tdelta)	// tau*vcvo.tdelta is usually just tau
						if ((colsum(svar)>0)				// there are current & lagged obs
								& (colsum(svar1)>0))	 {
							wv  = select((*vcvo.e),svar)  :* select((*vcvo.wvar),svar)  * vcvo.wf
							wv1 = select((*vcvo.e),svar1) :* select((*vcvo.wvar),svar1) * vcvo.wf
							Zsub =select((*vcvo.Z),svar)
							Zsub1=select((*vcvo.Z),svar1)
							if (K==1) {						// simple/fast, e is column vector
								eZsub = quadcross(1, wv, Zsub)		// equivalent to colsum(wv :* Zsub)
								eZsub1= quadcross(1, wv1, Zsub1)	// equivalent to colsum(wv :* Zsub)
							}
							else {
								eZsub=J(1,L*K,0)
								for (j=1; j<=rows(Zsub); j++) {
									wvj =wv[j,.]
									Zj  =Zsub[j,.]
									eZsub=eZsub+(wvj#Zj)
								}
								eZsub1=J(1,L*K,0)
								for (j=1; j<=rows(Zsub1); j++) {
									wv1j =wv1[j,.]
									Z1j  =Zsub1[j,.]
									eZsub1=eZsub1+(wv1j#Z1j)
								}
							}
							ghat=quadcross(eZsub,eZsub1)
							shat2=shat2+kw*(ghat+ghat')
						}
						i=i+vcvo.tdelta
					}
				}	// end non-zero kernel weight block
			}	// end tau loop

// If 1-level clustering, shat2 just calculated above is actually the desired shat
			if (vcvo.clustvarname2=="") {
				shat=shat2
			}
		}

// 2-level clustering, completion
// Cameron-Gelbach-Miller/Thompson method:
// Add 2 cluster variance matrices and subtract 3rd
		if (vcvo.clustvarname2~="") {
			shat = shat+shat2-shat3
		}

// Note no dof correction required for cluster-robust
	shat=shat/vcvo.N
	} // end cluster-robust code

	if (vcvo.sw~="") {
// Stock-Watson adjustment.  Calculate Bhat in their equation (6).  Also need T=panel length.
// They define for balanced panels.  Since T is not constant for unbalanced panels, need
// to incorporate panel-varying 1/T, 1/(T-1) and 1/(T-2) as weights in summation.

		st_view(ivar, ., st_tsrevar(vcvo.ivarname), vcvo.touse)
		info_ivar = panelsetup(ivar, 1)

		shat=J(L*K,L*K,0)
		bhat=J(L*K,L*K,0)
		N_panels=0
		for (i=1; i<=rows(info_ivar); i++) {
			esub=panelsubmatrix(*vcvo.e,i,info_ivar)
			Zsub=panelsubmatrix(*vcvo.Z,i,info_ivar)
			wsub=panelsubmatrix(*vcvo.wvar,i,info_ivar)
			Tsub=rows(esub)
			if (Tsub>2) {			// SW cov estimator defined only for T>2
				N_panels=N_panels+1
				sigmahatsub=J(K,K,0)
				ZZsub=J(L*K,L*K,0)
				shatsub=J(L*K,L*K,0)
				for (j=1; j<=rows(esub); j++) {
					eZi=esub[j,1]#Zsub[j,.]
					if ((vcvo.weight=="fweight") | (vcvo.weight=="iweight")) {
						shatsub=shatsub+quadcross(eZi,eZi)*wsub[j]*vcvo.wf
						sigmahatsub=sigmahatsub + quadcross(esub[j,1],esub[j,1])*wsub[j]*vcvo.wf
						ZZsub=ZZsub+quadcross(Zsub[j,.],Zsub[j,.])*wsub[j]*vcvo.wf
					}
					else {
						shatsub=shatsub+quadcross(eZi,eZi)*((wsub[j]*vcvo.wf)^2)
						sigmahatsub=sigmahatsub + quadcross(esub[j,1],esub[j,1])*((wsub[j]*vcvo.wf)^2)
						ZZsub=ZZsub+quadcross(Zsub[j,.],Zsub[j,.])*((wsub[j]*vcvo.wf)^2)
					}
				} // end loop through j obs of panel i
				shat=shat + shatsub*(Tsub-1)/(Tsub-2)
				bhat=bhat + ZZsub/Tsub#sigmahatsub/(Tsub-1)/(Tsub-2)
			}
		} // end loop through i panels

// Note that Stock-Watson incorporate an N-n-k degrees of freedom correction in their eqn 4
// for what we call shat.  We use only an N-n degrees of freedom correction, i.e., we ignore
// the k regressors.  This is because this is an estimate of S, the VCV of orthogonality conditions,
// independently of its use to obtain an estimate of the variance of beta.  Makes no diff aysmptotically.
// Ignore dofminus correction since this is explicitly handled here.
// Use number of valid panels in denominator (SW cov estimator defined only for panels with T>2).
		shat=shat/(vcvo.N-N_panels)
		bhat=bhat/N_panels
		shat=shat-bhat
	} // end Stock-Watson block

	_makesymmetric(shat)

// shat may not be positive-definite.  Use spectral decomposition to obtain an invertable version.
// Extract Eigenvector and Eigenvalues, replace EVs, and reassemble shat.
// psda option: Stock-Watson 2008 Econometrica, Remark 8, say replace neg EVs with abs(EVs).
// psd0 option: Politis (2007) says replace neg EVs with zeros.
	if (vcvo.psd~="") {
		symeigensystem(shat,Evec,Eval)
		if (vcvo.psd=="psda") {
			Eval = abs(Eval)
		}
		else {
			Eval = Eval + (abs(Eval) - Eval)/2
		}
		shat = Evec*diag(Eval)*Evec'
		_makesymmetric(shat)
	}

	return(shat)

} // end of program m_actest9_omega

// *********************************************************************** //
// *********************************************************************** //

real scalar m_actest9_calckw(	real scalar tau,
							real scalar bw,
							string scalar kernel)
	{
				karg = tau / bw
				if (kernel=="Truncated") {
					kw=1
				}
				if (kernel=="Bartlett") {
					kw=(1-karg)
				}
				if (kernel=="Parzen") {
					if (karg <= 0.5) {
						kw = 1-6*karg^2+6*karg^3
					}
					else {
						kw = 2*(1-karg)^3
					}
				}
				if (kernel=="Tukey-Hanning") {
					kw=0.5+0.5*cos(pi()*karg)
				}
				if (kernel=="Tukey-Hamming") {
					kw=0.54+0.46*cos(pi()*karg)
				}
				if (kernel=="Tent") {
					kw=2*(1-cos(tau*karg)) / (karg^2)
				}
				if (kernel=="Danielle") {
					kw=sin(pi()*karg) / (pi()*karg)
				}
				if (kernel=="Quadratic Spectral") {
					kw=25/(12*pi()^2*karg^2) /*
						*/ * ( sin(6*pi()*karg/5)/(6*pi()*karg/5) /*
						*/     - cos(6*pi()*karg/5) )
				}
				return(kw)
	}  // end kw

// *********************************************************************** //
// *********************************************************************** //


end

* Version notes
* Cumby-Huizinga-Arellano-Bond general specification test of serial correlation after OLS/IV/GMM estimation
* Cumby-Huizinga, NBER Technical Working Paper 90, 1990.
* Cumby-Huizinga, Econometrica, 60:1, 1992, 185-195
* Arellano-Bond, Review of Economic Studies, 58:2, 1991, 277-297
* 1.0.0: 27Jan2007
* 1.0.1: 31Jan2007 switch to Bartlett weights if Psi not p.d.;
*        allow for N=0 (regress without robust vce)
* 1.0.2: 02Feb2007 Fold in ivsc logic for i.i.d. errors
* 1.0.3: 04Feb2007 correct handling of A, disallow cluster
* 1.0.4: 14jul2007 guard against multiple panels, corrgram under if
* 1.0.5: 19Aug2007 Added trap for partialling-out by ivreg2
* 2.0.0: 22Dec2012 Complete rewrite
* 2.0.1: 30Dec2012 Explicit bw(0) allowed. Report if S-W eigenvalue adjustment used.
*                  2-way clustering support added.
* 2.0.2: 01Mar2013 Bug fix for capturing det(psi)<0 => use S-W method to force PSD
* 2.0.3: 24Jun2013 Redubbed actest. Support for univariate test added.
* 2.0.04 27Jun2013 Renaming/reconfigation of user-specified options and options inherited from estimation
* 2.0.05 03Jul2013 Added strictexog, q0, bp, small. Removed cmd arg to m_actest.  Restructured reporting, modified tests.
* 2.0.06 04Jul2013 Output indicates which test statistics required npsd adjustment
* 2.0.07 23Jul2013 Output cleanup
* 2.0.08 23Jul2013 Catch "last estimates not found" error
* 2.0.09 25Jul2013 Modify tsrevar calls to fvrevar calls
* 2.0.10 21Jan2014 Fixed reporting bug with 2-way clustering and kernel-robust
*                  that would give wrong count for 2nd cluster variable.
*                  Fixed bug that would cause crash if TS operators used with cluster.
*                  Tweaked output header so that BP, LB, BG, AB tests cited.
*                  Changed lags option to take lag as minimal abbreviation
* 2.0.11 22jan2014 Allow TS ops in depvar
* 2.0.12 24jan2014 Fixed minor bug in footnotes (wouldn't announce strict/weak exogeneity)
* 2.0.13 22Jan2013 First version of actest9. Mata library now internal with names incorporating "_actest9_".
*                  Changed subroutine name from acstat to s_acstat.

* Version notes for imported version of Mata library
* 1.1.01     First version of library.
*            Contains struct ms_vcvorthog, m_omega, m_calckw, s_vkernel.
*            Compiled in Stata 9.2 for compatibility with ranktest 1.3.01 (a 9.2 program).
* 1.1.02     Add routine cdsy. Standardized spelling/caps/etc. of QS as "Quadratic Spectral"
* 1.1.03     Corrected spelling of "Danielle" kernel in m_omega()
* 1.1.04     Fixed weighting bugs in robust and cluster code of m_omega where K>1
* 1.1.05     Added whichlivreg2(.) to aid in version control
* 1.1.06     Fixed remaining weighting bug (see 1.1.04) in 2-way clustering when interection
*            of clustering levels is groups
* 1.1.07     Fixed HAC bug that crashed m_omega(.) when there were no obs for a particular lag