clear all cls local N = 10000 // Sample size local R = 5000 // Simulation repetitions local L = 20000000 // Sample size for approximate true value calculations local S = 10 // Repetitions for approximate true value calculations set seed 111 program define mkdata syntax, [n(integer 1000)] clear quietly set obs `n' // 1. Generating data from probit, logit, and misspecified generate x1 = rchi2(2)-2 generate x2 = rbeta(4,2)>.2 generate u = runiform() generate e = ln(u) -ln(1-u) generate ep = rnormal() generate xb = .5*(1 - x1 + x2) generate y = xb + e > 0 generate yp = xb + ep > 0 // 1. Computing probit & logit marginal and treatment effects generate m1 = exp(xb)*(-.5)/(1+exp(xb))^2 generate m2 = exp(1 -.5*x1)/(1+ exp(1 -.5*x1 )) - /// exp(.5 -.5*x1)/(1+ exp(.5 -.5*x1 )) generate m1p = normalden(xb)*(-.5) generate m2p = normal(1 -.5*x1 ) - normal(.5 -.5*x1) // 4. Computing probit & logit marginal and treatment effects at means quietly mean x1 x2 matrix A = r(table) scalar a = .5 -.5*A[1,1] + .5*A[1,2] scalar b1 = 1 -.5*A[1,1] scalar b0 = .5 -.5*A[1,1] generate mean1 = exp(a)*(-.5)/(1+exp(a))^2 generate mean2 = exp(b1)/(1+ exp(b1)) - exp(b0)/(1+ exp(b0)) generate mean1p = normalden(a)*(-.5) generate mean2p = normal(b1) - normal(b0) end mkdata, n(`L') local values "m1 m2 mean1 mean2 m1p m2p mean1p mean2p" local means "mx1 mx2 meanx1 meanx2 mx1p mx2p meanx1p meanx2p" local n : word count `values' forvalues i= 1/`n' { local a: word `i' of `values' local b: word `i' of `means' sum `a', meanonly local `b' = r(mean) } postfile lpm y1l y1l_r y1lp y1lp_r y2l y2l_r y2lp y2lp_r /// using simslpm, replace forvalues i=1/`R' { quietly { mkdata, n(`N') logit y x1 i.x2, vce(robust) margins, dydx(*) atmeans post vce(unconditional) local y1l = _b[x1] test _b[x1] = `meanx1' local y1l_r = (r(p)<.05) local y2l = _b[1.x2] test _b[1.x2] = `meanx2' local y2l_r = (r(p)<.05) regress y x1 i.x2, vce(robust) margins, dydx(*) atmeans post vce(unconditional) local y1lp = _b[x1] test _b[x1] = `meanx1' local y1lp_r = (r(p)<.05) local y2lp = _b[1.x2] test _b[1.x2] = `meanx2' local y2lp_r = (r(p)<.05) post lpm (`y1l') (`y1l_r') (`y1lp') (`y1lp_r') /// (`y2l') (`y2l_r') (`y2lp') (`y2lp_r') } if (`i'/50) == int(`i'/50) { di ". `i'" } else { di _c "." } } postclose lpm use simslpm, clear sum postfile lp y1l y1l_r y1lp y1lp_r y2l y2l_r y2lp y2lp_r /// using simslp, replace forvalues i=1/`R' { quietly { mkdata, n(`N') logit y x1 i.x2, vce(robust) margins, dydx(*) post vce(unconditional) local y1l = _b[x1] test _b[x1] = `mx1' local y1l_r = (r(p)<.05) local y2l = _b[1.x2] test _b[1.x2] = `mx2' local y2l_r = (r(p)<.05) regress y x1 i.x2, vce(robust) margins, dydx(*) post vce(unconditional) local y1lp = _b[x1] test _b[x1] = `mx1' local y1lp_r = (r(p)<.05) local y2lp = _b[1.x2] test _b[1.x2] = `mx2' local y2lp_r = (r(p)<.05) post lp (`y1l') (`y1l_r') (`y1lp') (`y1lp_r') /// (`y2l') (`y2l_r') (`y2lp') (`y2lp_r') } if (`i'/50) == int(`i'/50) { di ". `i'" } else { di _c "." } } postclose lp use simslp, clear sum postfile lppm y1l y1l_r y1lp y1lp_r y2l y2l_r y2lp y2lp_r /// using simslppm, replace forvalues i=1/`R' { quietly { mkdata, n(`N') probit yp x1 i.x2, vce(robust) margins, dydx(*) atmeans post vce(unconditional) local y1l = _b[x1] test _b[x1] = `meanx1p' local y1l_r = (r(p)<.05) local y2l = _b[1.x2] test _b[1.x2] = `meanx2p' local y2l_r = (r(p)<.05) regress yp x1 i.x2, vce(robust) margins, dydx(*) atmeans post vce(unconditional) local y1lp = _b[x1] test _b[x1] = `meanx1p' local y1lp_r = (r(p)<.05) local y2lp = _b[1.x2] test _b[1.x2] = `meanx2p' local y2lp_r = (r(p)<.05) post lppm (`y1l') (`y1l_r') (`y1lp') (`y1lp_r') /// (`y2l') (`y2l_r') (`y2lp') (`y2lp_r') } if (`i'/50) == int(`i'/50) { di ". `i'" } else { di _c "." } } postclose lppm use simslppm, clear sum postfile lpp y1l y1l_r y1lp y1lp_r y2l y2l_r y2lp y2lp_r /// using simslpp, replace forvalues i=1/`R' { quietly { mkdata, n(`N') probit yp x1 i.x2, vce(robust) margins, dydx(*) post vce(unconditional) local y1l = _b[x1] test _b[x1] = `mx1p' local y1l_r = (r(p)<.05) local y2l = _b[1.x2] test _b[1.x2] = `mx2p' local y2l_r = (r(p)<.05) regress yp x1 i.x2, vce(robust) margins, dydx(*) post vce(unconditional) local y1lp = _b[x1] test _b[x1] = `mx1p' local y1lp_r = (r(p)<.05) local y2lp = _b[1.x2] test _b[1.x2] = `mx2p' local y2lp_r = (r(p)<.05) post lpp (`y1l') (`y1l_r') (`y1lp') (`y1lp_r') /// (`y2l') (`y2l_r') (`y2lp') (`y2lp_r') } if (`i'/50) == int(`i'/50) { di ". `i'" } else { di _c "." } } postclose lpp use simslpp, clear sum