* Born and Breitung (2016) HR-test for serial correlation
*! Version 1.0.2 14sep2016
* Contact jesse.wursten@kuleuven.be for bug reports/inquiries.

/*
cap mata mata drop hr_statistic()
cap mata mata drop hr_statistic_unbalanced()
cap program drop xthrtest
*/

program define xthrtest, rclass
	version 12

	** Technicalities
	syntax [varlist(default=none)] [if] [in] [, force]

	*** Postestimation?
	tempvar residuals
	if "`varlist'" == "" {
		predict `residuals', ue
		local varlist = "`residuals'"
		local postEstimation = "1"
	}

	marksample toUse, novarlist

	*** Obtain time and panel variables
	qui xtset
	local panelvar = r(panelvar)
	local timevar = r(timevar)


	** Print results header
	if "`postEstimation'" == "" di as result _newline "Heteroskedasticity-robust Born and Breitung (2016) HR-test on `varlist'"
	else 						di as result _newline "Heteroskedasticity-robust Born and Breitung (2016) HR-test as postestimation"
	di as text "Panelvar: `panelvar'"
	di as text "Timevar: `timevar'"
	di as text "{hline 30}{c TT}{hline 23}{c TT}{hline 16}{c TT}{hline 14}{c TRC}"
	di as text _col(2) %~28s = "Variable"  _col(30) " {c |}" _skip(2) "HR-stat" _skip(4) "p-value" _skip(2) " {c |}" _skip(6) "N" _skip(4) "maxT" " {c |}" %~14s = "balance?" "{c |}"
	di as text "{hline 30}{c +}{hline 23}{c +}{hline 16}{c +}{hline 14}{c RT}"

	** Calculate statistic
	tempname hr pvalue
	tempvar t id mean_resid toUse2
	local j = 1
	foreach var of local varlist {
		*** Obtain number of time and panel units
		qui egen `t' = group(`timevar') if `toUse' == 1 & ~missing(`var')
		sum `t', meanonly
		local timelength = r(max)

		qui egen `id' = group(`panelvar') if `toUse' == 1 & ~missing(`var')
		sum `id', meanonly
		local panelunits = r(max)
		drop `t' `id'

		*** Check balancedness
		**** Unbalanced
		qui count if ~missing(`var') & `toUse' == 1
		local totalObsUsed = r(N)
		local requiredObs = `timelength'*`panelunits'

		if `totalObsUsed' != `requiredObs' 	local balance = "unbalanced"

		**** With gaps
		tempvar gap
		qui bysort `panelvar' (`timevar'): gen `gap' = 1 if missing(`var') & ~missing(L.`var', F.`var') & `toUse' == 1
		qui count if `gap' == 1
		if r(N) != 0 local balance = "gaps (error)"

		**** Balanced
		if "`balance'" == "" local balance = "balanced"

		*** Unless force is specified ...
		if "`force'" == "" {
			* Test if residuals include the fixed effect
			tempvar mean_resid
			qui bysort `panelvar' (`timevar'): egen `mean_resid' = mean(`var') if `toUse' == 1
			qui sum `mean_resid' if `toUse' == 1
			drop `mean_resid'
			if r(sd) < 0.001 | abs(r(max)) < 0.001 | abs(r(min)) < 0.001 {
				*noisily di as error _col(2) %~28s = abbrev("`var'",28)
				noisily di as error "`var': " _continue
				noisily di as error _col(4) "Residuals do not appear to include the fixed effect."
				noisily di as error _col(4) "This test is made to function with ue = c_i + e_it."
				noisily di as error _col(4) "If you are sure that your residuals do indeed include" _newline _col(4) "the fixed effect (programming bugs happen),"
				noisily di as error _col(4) "specify 'force' to skip this test."
				continue
			}
		}

		*** Exclude missings (balanced case only)
		qui gen `toUse2' = `toUse'
		qui replace `toUse2' = 0 if missing(`var')

		*** Run calculations
		if "`balance'" == "balanced" mata: hr_statistic("`var'", "`toUse2'", `timelength', `panelunits')
		else if "`balance'" == "unbalanced" {
			tempvar originalObs
			gen byte `originalObs' = 1
			tsfill, full
			mata: hr_statistic_unbalanced("`var'", "`toUse'", `panelunits')
			qui drop if `originalObs' != 1
		}
		scalar `hr' = round(HR, 0.001)
		scalar `pvalue' = round(pvalue, 0.001)

		else if "`balance'" == "gaps (error)" {
			scalar `hr' = .
			scalar `pvalue' = .
		}

		if "`postEstimation'" == "1" local var = "Post Estimation"
		di _col(2) %~28s = abbrev("`var'",28) _col(30)  " {c +}" _skip(3) %4.2f = `hr' _col(46) %4.3f = `pvalue' _col(55) "{c +}" %7.0f = `panelunits' %8.0f = `timelength' " {c +}" %~14s = "`balance'" "{c RT}"

		** Prep for return
		drop `toUse2'
		mat hrs = (nullmat(hrs), HR)
		mat ps = (nullmat(ps), pvalue)

		** Return as scalar (more useful in some cases)
		return scalar hr`j' = HR
		return scalar pvalue`j' = pvalue
		local j = `j' + 1
	}

	** Notes
	di as text "{hline 30}{c BT}{hline 23}{c BT}{hline 16}{c BT}{hline 14}{c BRC}"
	di _col(2) "Notes: Under H0, HR ~ N(0,1)"
	di as txt _col(5) "H0: No first-order serial correlation."
	di as txt _col(5) "Ha: Some first order serial correlation."

	** Return
	return matrix HR = hrs
	return matrix p = ps
end

mata:
	void hr_statistic(string scalar varname, string scalar toUse, real scalar T, real scalar N) {
		real matrix UE, V0, V1, diagMatrix, zeroMaker; real scalar HR, pvalue

		// ue = fe + e (currently we just use e, as it also seems to work)
		// me = e
		UE = rowshape(st_data(., varname, toUse), N)'

		// V0
		V0 = J(T-3, T, 0) :- 1:/range(2, T-2, 1)

		diagMatrix = J(T-3, 1, 0), I(T-3), J(T-3, 2, 0)

		zeroMaker = J(T-3, 1, .), lowertriangle(J(T-3, T-3, 1)), J(T-3, 2, 0)
		_editvalue(zeroMaker, 1, .)
		zeroMaker = zeroMaker :+ 1:/range(2, T-2, 1)
		_editmissing(zeroMaker, 0)

		V0 = V0 + diagMatrix + zeroMaker


		// V1
		V1 = J(T-3, T, 0) :- 1:/range(T-2, 2, -1)

		diagMatrix = J(T-3, 2, 0), I(T-3), J(T-3, 1, 0)

		zeroMaker = J(T-3, 2, 0), uppertriangle(J(T-3, T-3, 1)), J(T-3, 1, .)
		_editvalue(zeroMaker, 1, .)
		zeroMaker = zeroMaker :+ 1:/range(T-2, 2, -1)
		_editmissing(zeroMaker, 0)

		V1 = V1 + diagMatrix + zeroMaker

		// At
		At = V0'V1

		// Z
		Z = J(1, N, .)
		for(n=1; n<=N; n++) {
			Z[1, n] = UE[., n]'*At*UE[.,n]
		}

		Z2 = (diagonal(UE'*At*UE))'

		// Qp-tilde
		HR = sum(Z)/sqrt(sum(Z:^2) - 1/N*(sum(Z))^2)
		pvalue = 2*(1-normal(abs(HR)))

		// Store results
		st_numscalar("HR", HR)
		st_numscalar("pvalue", pvalue)
	}
end

mata:
	void hr_statistic_unbalanced(string scalar varname, string scalar toUse, real scalar N) {
		real matrix UE, V0, V1, diagMatrix, zeroMaker; real scalar HR, pvalue

		UE = rowshape(st_data(., varname, toUse), N)'
		Ti = colnonmissing(UE)

		Z = J(1, N, .)
		for(n=1; n<=N; n++) {
			V0 = J(Ti[n]-3, Ti[n], 0) :- 1:/range(2, Ti[n]-2, 1)
			diagMatrix = J(Ti[n]-3, 1, 0), I(Ti[n]-3), J(Ti[n]-3, 2, 0)
			zeroMaker = J(Ti[n]-3, 1, .), lowertriangle(J(Ti[n]-3, Ti[n]-3, 1)), J(Ti[n]-3, 2, 0)
			_editvalue(zeroMaker, 1, .)
			zeroMaker = zeroMaker :+ 1:/range(2, Ti[n]-2, 1)
			_editmissing(zeroMaker, 0)

			V0 = V0 + diagMatrix + zeroMaker


			V1 = J(Ti[n]-3, Ti[n], 0) :- 1:/range(Ti[n]-2, 2, -1)
			diagMatrix = J(Ti[n]-3, 2, 0), I(Ti[n]-3), J(Ti[n]-3, 1, 0)
			zeroMaker = J(Ti[n]-3, 2, 0), uppertriangle(J(Ti[n]-3, Ti[n]-3, 1)), J(Ti[n]-3, 1, .)
			_editvalue(zeroMaker, 1, .)
			zeroMaker = zeroMaker :+ 1:/range(Ti[n]-2, 2, -1)
			_editmissing(zeroMaker, 0)

			V1 = V1 + diagMatrix + zeroMaker

			Ati = V0'V1
			UEi = select(UE[.,n], UE[.,n]:!=.)
			Z[1, n] = UEi'*Ati*UEi
		}

		// Qp-tilde
		HR = sum(Z)/sqrt(sum(Z:^2) - 1/N*(sum(Z))^2)
		pvalue = 2*(1-normal(abs(HR)))

		// Store results
		st_numscalar("HR", HR)
		st_numscalar("pvalue", pvalue)
	}
end