* This example is in Cameron and Trivedi, ch 5
clear
* Generated data for heteroskedasticity
set seed 10101
quietly set obs 500
generate double x2 = 5 * rnormal(0)
generate double x3 = 5 * rnormal(0)
generate double e = 5 * rnormal(0)
generate double u = sqrt(exp(-1.+0.2*x2))*e
generate double y = 1 + 1*x2+1*x3+u

summarize

* OLS with default (inconsistent) std errors
reg y x2 x3

* OLS with robust std errors
reg y x2 x3, vce(hc3)

* Diagnostic plots
quietly reg y x2 x3
predict double uhat, resid
generate double absu = abs(uhat)
quietly twoway (scatter absu x2) (lowess absu x2, bw(0.4) lw(thick)),scale(1.2) xscale(titleg(*5)) yscale(titleg(*5)) plotr(style(none)) name(gls1)
quietly twoway (scatter absu x3) (lowess absu x3, bw(0.4) lw(thick)),scale(1.2) xscale(titleg(*5)) yscale(titleg(*5)) plotr(style(none)) name(gls2)
graph combine gls1 gls2

* Breusch-Pagan test
reg y x2 x3
hettest x2, iid

* Manual Version
quietly regress y x2 x3
predict ehat, res
gen ehat2=ehat*ehat
quietly reg ehat2 x2
scalar LM = e(N)*e(r2)
scalar pvalue = chi2tail(1,LM)
scalar list LM pvalue

* Whites Test
quietly reg y x2 x3
estat imtest, white 

* FGLS estimation
drop uhat
quietly reg y x2 x3
predict double uhat, resid
generate double uhatsq = uhat*uhat
generate double one = 1
 nl (uhatsq = exp({xb: x2 one})), nolog
 predict double varu, yhat
 
 regress y x2 x3 [aweight=1/varu]