 program define heck
 version 6.0
 syntax [, obs(integer 1000) vcm(string)]
 cap log close
 set log l 95
 set d l 95
 set more off
 log using heck_xmp, replace
 clear
 set obs `obs'
 di _n(2) _dup(80) "*" _n(2) "Simulate the heckman/lee selection model"
 di "Step 1." _n(2) _col(10) "Generate the x variables and pick the coefficient vectors. Note exclusion restrictions" _n(2) _dup(80) "*" _n(2)
 gen x2 = invnorm(uniform())
 gen x3 = invnorm(uniform())
 gen x4 = invnorm(uniform())
 summ x2 x3 x4
 matrix b1 = (2.0,  0.8, -1.5, 0.0)
 mat coln b1 = _cons x2 x3 x4
 matrix b2 = (1.0, 0.0, -2.0, 0.05)
 mat coln b2 = _cons x2 x3 x4
 matrix l b1 , title("True coefficients in the regression equation")
 matrix l b2,  title("True coefficients in the selection equation")
 di _n(2) _dup(80) "*" _n(2) "Step 2. " _n(2) _col(10) "Generate independent normal random numbers (pseudo)" _n(2) _dup(80) "*" _n(2)
 gen u1 = invnorm(uniform())
 gen u2 = invnorm(uniform())
 corr u1 u2, covar means
 cap matrix l `vcm'
 if _rc {
   tempname v
   matrix `v'= ( 5.0 , 3 \ 3 , 2.0)
   loc vcm "`v'"
   }
 di _n(2) _dup(80) "*" _n(2) "Step 3. " _n(2) _col(10) " Choose a variance matrix for the epsilons"
 di  _col(10) " this matrix must satisfy the normalization that the REDUCED-FORM"
 di  _col(10) " probit has variance=1, which requires  abs(v11 + v22 - 2v12) = 1" _n(2) _dup(80) "*" _n(2)
 matrix l `vcm' , title("True variance matrix for structural errors")
 di in w _dup(80) "*" _n(2) "Step 4. " _n(2)  _col(10) " Generate the correlated errors using Cholesky decomposition of v" _n(2) _dup(80) "*" _n(2)
 tempname chol
 matrix `chol' = cholesky(`vcm')
 mat l `chol' , title("Cholesky Decomp. of True variance matrix")
 gen eps1 = u1*`chol'[1,1] + u2*`chol'[1,2]
 gen eps2 = u1*`chol'[2,1] + u2*`chol'[2,2]
 corr eps1 eps2, mean covar
 graph eps1 eps2, title(Correlated Standard Normals) saving(corr,replace)
 di _n(2) _dup(80) "*" _n(2) "Step 5. " _n(2)  _col(10) "  Generate x\beta and the y variables" _n(2) _dup(80) "*" _n(2)
 set trace on
 mat score xb1 = b1
 mat score xb2 = b2
 gen yst1 = xb1 + eps1
 gen yst2 = xb2 + eps2
 set trace off
 summ xb1 xb2 yst1 yst2
 di _n(2) _dup(80) "*" _n(2) "Step 6. " _n(2) _col(10) "  Check the data: OLS should return the b1 and b2 vectors (nearly)" _n(2) _dup(80) "*" _n(2)
 regress yst1 x2 x3 x4
 regress yst2 x2 x3 x4
 di _n(2) _dup(80) "*" _n(2) "Step 7. " _n(2)  _col(10) " Carry out the sample selection" _n(2) _dup(80) "*" _n(2)
 gen byte wking = cond(yst1>yst2,1,0)
 gen y = yst1 if wking==1
 summ y wking
 count if yst1>yst2
 tab wking
 save heck , replace
di _n(2) _dup(80) "*" _n(2) "Step 8. " _n(2) _col(10) " Bias in OLS if selection is ignored " _n(2) _dup(80) "*" _n(2)
regress y x2 x3 x4
di _n(2) _dup(80) "*" _n(2) "Step 9. " _n(2) _col(10) " Heckman 2-step: reduced form probit, inverse mills, corrected OLS"
di _col(10) "compare results to structural coefficients - why are they different" _n(2) _dup(80) "*" _n(2)
probit wking x2 x3 x4
predict wkhat, xb
mat pihat = get(_b)
gen lamhat = -normd(wkhat) / normprob(wkhat)
regress y x2 x3 lamhat
mat bhat = get(_b)
di _n(2) _dup(80) "*" _n(2) _col(10) "Estimated Coeff on x2 in first equation would be"
di _col(10) bhat[1,1] "-" pihat[1,1] "=" bhat[1,1]-pihat[1,1]
di _n(2) _col(10) "Step 10. " _n(2) _col(10) " FIML estimates of the selection model"  _n(2) _dup(80) "*" _n(2)
heckman y x2 x3,  select(x3 x4)
end