* use "C:\data\STATA\participation.dta", clear
use "H:\Documents and Settings\Lee\My Documents\Document\stata\participation.dta", clear

gen A2 = A*A

quietly probit lfp I A A2 E F no nu, technique(nr)
predict index, xb
predict prb, pr

generate Ft = prb
generate ft = normalden(index)

/* variance */
gen vt = Ft*(1-Ft)

/* transform variables for OPG regression */
gen ystar=(lfp-Ft)/sqrt(vt)
gen Istar=ft*I/sqrt(vt)
gen Astar=ft*A/sqrt(vt)
gen A2star=ft*A2/sqrt(vt)
gen Estar=ft*E/sqrt(vt)
gen nustar=ft*nu/sqrt(vt)
gen nostar=ft*no/sqrt(vt)
gen Fstar=ft*F/sqrt(vt)
gen cons = ft/sqrt(vt)

gen Istaro=index*ft*I/sqrt(vt)
gen Astaro=index*ft*A/sqrt(vt)
gen A2staro=index*ft*A2/sqrt(vt)
gen Estaro=index*ft*E/sqrt(vt)
gen nustaro=index*ft*nu/sqrt(vt)
gen nostaro=index*ft*no/sqrt(vt)
gen Fstaro=index*ft*F/sqrt(vt)
gen conso =index*ft/sqrt(vt)

/* BRMR for heteroskedasticity related to all variables in x */
reg ystar cons Istar Astar A2star Estar nustar nostar Fstar Istaro Astaro A2staro Estaro nustaro nostaro Fstaro, noconst 
scalar ess = e(mss)
scalar pval = chi2tail(7,ess)
scalar list ess pval

/* BRMR 'RESET' -- though the results don't quite match those of DM? */
gen index2=index*index*ft/sqrt(vt)
reg ystar cons Istar Astar A2star Estar Fstar nustar nostar index2, noconst 
scalar ess = e(mss)
scalar pval = chi2tail(1,ess)
scalar list ess pval