*! version 1: 10 April 2021, Gueorgui I. Kolev
*! Computes robust and/or clustered variance matrix of estimates post xtgls.
*! Methods and Formulas are in Kolev (2014) "Robust variance estimation in panel data
*! generalized least squares regression."
program define xtglsr
version 11
preserve
syntax [, CLUSTER(varname) MINUS(integer 0)]
if "`e(cmd)'" != "xtgls" {
display as error "Can be used only after xtgls"
error 301
}
qui xtset
local panelis `r(panelvar)'
local timeis `r(timevar)'
if "`r(panelvar)'"=="" | "`r(timevar)'"=="" {
display as error "You need to xtset your data, and specify both panel and time identifiers"
error 459
}
* If the option cluster is specified, restrict sample to non missing cluster var.
if !missing("`cluster'") qui drop if missing(`cluster')
* Tempnames for the residuals, Inverse Sigma, and the weighted residuals.
tempname residual InvSigma weightedresidual
tempfile thedata
* Replay in case somebody is calling xtgls repeatedly for the correct homoskedastic variance.
quietly {
`e(cmdline)'
mat `InvSigma' = invsym(e(Sigma)) // Sigma^-1, Sigma is the covariance of errors across equations.
egen Individual_Observations = group(`timeis' `panelis') if e(sample)
* Generate the residuals
predict double `residual' if e(sample), xb // only the linear prediction is available post xtgls
replace `residual' = `e(depvar)' - `residual'
* Generate the weighted residual
save `thedata', replace
keep `residual' `panelis' `timeis'
levelsof `panelis', local(levelspanelis) // Extract the levels: panelvar might be something bizarre, like -3, 0, 5, 300, etc.
* We reshape the data to wide. matrix score operates on variables.
reshape wide `residual', i(`timeis') j(`panelis')
ds `timeis', not // at the point in the data there are only `residual'-stub variables, and `timeis'.
local theresiduals `r(varlist)'
mat colnames `InvSigma' = `theresiduals' // assign the names of the residuals as column names of Sigma^-1.
forvalues i = 1/`=rowsof(`InvSigma')' { // Loop over the rows of Sigma^-1.
matrix tempvec = `InvSigma'[`i',1...] // Extract the given row.
matrix score `weightedresidual'`: word `i' of `levelspanelis''= tempvec // Form a linear combination
// of the unweighted residual vectors u1 u2 .. uM, where the weights are given by the row of Sigma^-1.
} // Closes the Forvalues loop.
drop `theresiduals'
reshape long `weightedresidual', i(`timeis') j(`panelis')
merge 1:1 `timeis' `panelis' using `thedata', nogenerate
* Once we have calculated the set of weighted residuals, we call _robust.
if "`cluster'" != "" { // If the user has taken control of the clustering, we let the user decide what the clustring is on.
_robust `weightedresidual', minus(`minus') cluster(`cluster')
}
else if e(vt) == "heteroskedastic with cross-sectional correlation" { // If the user has not taken control of clustering,
// and the estimator is cross sectionally correlated, we cluster by time.
_robust `weightedresidual', minus(`minus') cluster(`timeis')
}
else _robust `weightedresidual', minus(`minus') cluster(Individual_Observations) // If the user has not specified clustering,
// and the estimator is NOT cross sectionally correlated, we cluster by observation, that is (Time X PanelId) identifier.
// I am doind this because _robust seems to be too smart for its own good: When I specify robust only, and the data is xtset
// _robust computes clustered by PanelID standard errors. In the context of -xtgls- this is incorrect and we do NOT want it.
} // Closes the Quietly brace.
xtgls, // We replay xtgls, but Robust Variance = Bread*Meat*Bread has been substituted
// by _robust for what was only Bread when we fit the xtgls initially.
end