clear
use "participation.dta", clear

gen A2=A*A

/* Probit */
probit lfp I A A2 E F 
/* Predicted values for index */
predict index, xb

/* Use normal() and normalden() to generate cdf and pdf */
generate Ft = normal(index)
generate ft = normalden(index)

/* variance */
gen vt = Ft*(1-Ft)
/* you may need a constant /*
gen one = 1

/* transform variables for BRMR *?
gen ystar=(lfp-Ft)

   [etc.]

/* Weighted least squares */
reg y x2 x3 [aweight=1/vt], noconst 
/* model sum of squares for test */
scalar ess = e(mss)
scalar pval = chi2tail(2,ess)
scalar list ess pval