set seed 547
# This generates the dataset for the Heckit Example in slides
# Change working directory first

# *************** Before the loop *********************
nulldata 5000
matrix SIG = {1,0.9;0.9,1}
matrix A = cholesky(SIG)

#****** open a loop, to be repeated R=500 times **********
#loop 500 --progressive --quiet

#************* variables creation ***********************
 genr kids = uniform()<0.6
 summary kids
 genr u=uniform()
 genr educ=8*(u<=.15)+12*(u>.15)*(u<=.6)+16*(u>.6)*(u<=.9)+21*(u>.9)
 genr u1 = normal()
 genr u2 = normal()
 genr u = A[1,1]*u1+A[1,2]*u2
 genr v = A[2,1]*u1+A[2,2]*u2
 genr Net_U = -0.5 + 0.03* educ -1.5*kids + v
 genr work=Net_U>=0
 genr wage_offer=5+0.07*educ+u
 genr wage=NA*(work=0)+wage_offer*(work=1)
 genr wage2=0*(work=0)+wage_offer*(work=1)
# estimation, testing, and fit 
 ols wage const educ --robust
 ols wage2 const educ --robust
 heckit wage const educ ; work const educ kids --two-step
 heckit wage const educ ; work const educ kids
 
restrict --quiet
     b[lambda]=0
 end restrict
 
 restrict --quiet
     b[5]=0.03
     b[6]=-1.5
 end restrict
 

 # marginal effects on wage offers of another year of education
 genr coeff=$coeff
 genr beta=coeff[1:2]
 genr gamma=coeff[4:6]
 series educ0=educ
 matrix x0={const,educ0}
 series educ1=educ+1
 matrix x1={const,educ1}
 series x1b = x1*beta
 series x0b = x0*beta
 genr Mg_educ = mean(x1b-x0b)

 # marginal effects on observed wages of another year of education
 genr coeff=$coeff
 genr beta=coeff[1:2]
 genr gamma=coeff[4:6]
 series educ0=educ
 matrix x0={const,educ0}
 matrix z0={const,educ0,kids}
 series w_hat0 = x0*beta+coeff[3]*pdf(N,z0*gamma)./cdf(N,z0*gamma)
 series educ1=educ+1
 matrix x1={const,educ1}
 matrix z1={const,educ1,kids}
 series w_hat1 = x1*beta+coeff[3]*pdf(N,z1*gamma)./cdf(N,z1*gamma)
 series Mg_educ2 = w_hat1-w_hat0
 summary Mg_educ2 --by=kids --simple 
 summary Mg_educ2 --simple

quit