rndseed 123456
!n=100         '  this is is the sample size
!nrep=5000    ' this is the number of MC replications
                        ' assign "true" values for the paramters
smpl 1 !n
scalar b1=1
scalar b2=2
series ystar
series y
series xbeta=b1+b2*x             ' use this to see how many xbeta values lie in unit interval

vector(!nrep) coefols
vector(!nrep) coefp
vector(!nrep) meff

                              ' next, calculate the TRUE marginal effect of X (at the sample mean)
scalar meff_true=@cnorm(b1+b2*@mean(x))*(1-@cnorm(b1+b2*@mean(x)))*b2
                               ' start the Monte Carlo loop
for !i=1 to !nrep
y=0
ystar=b1+b2*x+@rnorm
smpl 1 !n if ystar>0
y=1
smpl 1 !n
equation eq1.ls y c x
coefols(!i)=c(2)        ' this is also the marginal effect for x in the LPM
equation eq2.binary y c x
coefp(!i)=c(2)     ' next, calculate the estimated marginal effect of X (at the sample mean)
meff(!i)=@cnorm(c(1)+c(2)*@mean(x))*(1-@cnorm(c(1)+c(2)*@mean(x)))*c(2)
next                      ' end of the Monte Carlo loop
' convert vectors to series, simply to facilitate plots
smpl 1 !nrep
mtos(coefols, coefolss)
mtos(coefp, coefps)
mtos(meff,meffs)