*
* Example 22.8 from page 786
*
open data tablef4-1[1].txt
data(format=prn,org=columns,skip=36) 1 753 lfp whrs kl6 k618 wa we ww rpwg hhrs ha he hw faminc mtr wmed wfed un cit ax
set agesq = wa**2
set axsq  = ax**2
set kids  = (kl6+k618)>0
*
* Probit for labor force participation. Compute the "lambda" statistics using the
* Mills option, and (minus the) "delta" statistics using DMills. The deltas are
* used in recomputing the covariance matrix.
*
ddv lfp
# constant wa agesq faminc we kids
prj(mills=lambda,dmills=delta,xvx=prbfitvar) fit
compute probitxx=%xx
set delta = -delta
*
* OLS on those in the labor force
*
linreg(smpl=lfp) ww
# constant ax axsq we cit
*
* Second step of Heckit procedure
*
linreg(smpl=lfp,title="Heckit w/o Covariance Correction") ww
# constant ax axsq we cit lambda
compute xxheckit=%xx
*
* Compute the corrected estimate of the variance and the corresponding estimate
* of rho**2.
*
sstats(mean,smpl=lfp) / delta>>deltabar
compute sigsq=%rss/%nobs+deltabar*%beta(6)**2
compute rhosq=%beta(6)**2/sigsq
*
* Covariance matrix recalculation assuming deltas are known.
*
set fiddle = (1-rhosq*delta)
mcov(nosquare,smpl=lfp,lastreg) / fiddle
linreg(create,lastreg,form=chisquared,$
   covmat=xxheckit*%cmom*xxheckit*sigsq,$
   title="Heckit with Partially Corrected Covariances")
*
* Covariance matrix recalculation allowing for the estimation error in the probit
* model. The calculation of the "Q" is done differently than shown, because the W
* var(gamma) W' can be (and was) computed using the XVX option on the PRJ
* immediately after the probit.
*
set fiddle2 = (1-rhosq*delta+rhosq*prbfitvar*delta**2)
mcov(nosquare,smpl=lfp,lastreg) / fiddle2
linreg(create,lastreg,form=chisquared,$
   covmat=xxheckit*%cmom*xxheckit*sigsq,$
   title="Heckit with Fully Corrected Covariances")
*
* ML by EM algorithm
*
frml(regressors,vector=b) wagefrml
# constant ax axsq we cit
frml(regressors,vector=g) lfpfrml
# constant wa agesq faminc we kids
nonlin sigsq b g
compute siguv=-.131
ddv lfp
# constant wa agesq faminc we kids
compute g=%beta
linreg(smpl=lfp) ww
# constant ax axsq we cit lambda
compute b=%xsubvec(%beta,1,5)
*
function emstep time
type real emstep
type integer time
*
local real lindex v lem sigu yem
compute lindex = lfpfrml(time)
compute v      = ww(time)-wagefrml(time)
compute lem    = lindex+siguv/sigsq*v
compute sigu   = sqrt(1-siguv**2/sigsq)
compute yem    = sigu*%mills(lem/sigu)
compute emstep = yem
end

frml emlikely = %if(lfp,%logdensity(||1.0|siguv,sigsq||,||emstep(t),ww-wagefrml||),log(%cdf(-lfpfrml)))
maximize(pmethod=simplex,piters=10) emlikely