*! 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