capture program drop _all program define treatoprobitsim_p, eclass *! version 1.0.0 cagregory version 11 tempname cuts local cutpts = e(cuts) mat `cuts' = J(1,`cutpts',0) forv i = 1/`cutpts' { mat `cuts'[1,`i']=_b[cut`i':_cons] } local nchoices=`cutpts'+1 local PFX forv i = 1/`nchoices'{ local PFX "`PFX' P1`i' P0`i' TE`i' TT`i' SETE`i' SETT`i' " } local myopts PTR XBOUT `PFX' LF //di "`myopts'" _pred_se "`myopts'" `0' if (`s(done)') { exit } local vtyp `s(typ)' local varn `s(varn)' local 0 `"`s(rest)'"' syntax [if] [in] [, `myopts' noOFFset ] marksample touse local pfx forv i = 1/`nchoices'{ local pfx "`pfx'`p1`i''`p0`i''`te`i''`tt`i''`sete`i''`sett`i''" } local type "`ptr'`xbout'`pfx'`lf'" //di "`type'" local ytr : word 1 of `e(depvar)' local yout : word 2 of `e(depvar)' tempvar zbp xbt psel qui _predict double `zbp' if `touse', xb eq(#1) qui _predict double `xbt' if `touse', xb eq(#2) qui g double `psel' = normal(`zbp'/(sqrt(2))) if `touse' tempname cutmat_tr cutmat_out lambda betatreat delta local facskew = e(facskew) local neginf = minfloat() local posinf = maxfloat() mat `cutmat_out'= `neginf',`cuts', `posinf' local `delta' = _b[lambda:_cons] local `betatreat' = _b[`yout':`ytr'] tempname xbtzero xbtone qui gen double `xbtzero' = `xbt'-``betatreat''*`ytr' qui gen double `xbtone' = `xbtzero'+``betatreat'' qui gen double `lambda' = `ytr'*``delta''+(1-`ytr')*(-``delta'') qui count if `touse' scalar nobs = `r(N)' mata : simobs = rows(_treatoprobit_rnd) mata : nobs = st_numscalar("nobs") mata : simxN = nobs-simobs local facmean = e(facmean) local facscale = e(facscale) local facdensity = e(facdensity) local startpoint = e(startpoint) local nsim = e(sesim) mata: st_numscalar("simxtra", simxN) //scalar list simxtra if (scalar(simxtra))>0 { if "`facdensity'"=="normal" { mata: _treatoprobit_rnd_x = `facscale'* /// invnormal(halton(simxN*_treatoprobit_S,1,`startpoint',0)) } if "`facdensity'"=="uniform" { mata: _treatoprobit_rnd_x = `facscale'*sqrt(12)* /// ((halton(simxN*_treatoprobit_S,1,`startpoint',0)):-0.5) } if "`facdensity'"=="chi2" { local k = `=8/(`facskew'*`facskew')' local sgn = `=sign(`facskew')' di `sgn' mata: _treatoprobit_rnd_x = `facscale'/sqrt(2*`k')*`sgn'* /// (invchi2(`k',halton(simxN*_treatoprobit_S, 1,`startpoint',0)):-`k') } if "`facdensity'"=="gamma"{ mata: _treatoprobit_rnd_x = invgammap(`facmean', halton(simxN*_treatoprobit_S,1, `startpoint',0)):-`facmean' } if "`facdensity'"=="logit" { mata: _treatoprobit_rnd_x= `facscale'*logit(halton(simxN*_treatoprobit_S,1,`startpoint',0)) } if "`facdensity'"=="lognormal" { mata: _treatoprobit_rnd_x = `facscale'*exp(invnormal(halton(simxN*_treatoprobit_S,1,`startpoint',0))) mata: _treatoprobit_mean = mean(_treatoprobit_rnd_x) mata: _treatoprobit_rnd_x = _treatoprobit_rnd_x:-_treatoprobit_mean } mata: _treatoprobit_rnd_x=colshape(_treatoprobit_rnd_x,_treatoprobit_S) mata: _treatoprobit_rnd = _treatoprobit_rnd\_treatoprobit_rnd_x } if (missing("`type'") | "`type'" == "p11") { if missing("`type'") noisily display as text "(option p11 assumed; Pr(`ytr'=1, `yout'=1)" local trtv=1 local outv=1 tempvar vtemp qui g double `vtemp'=. //sum `zbp' qui generate `vtyp' `varn' = . if `touse' mata: mata_treatoprobitsim_predict("`vtemp'", "`zbp'", "`xbt'", "`cutmat_out'", "`lambda'", /// _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`trtv'", "`outv'" ) qui replace `varn'=`vtemp' if `touse' label variable `varn' "Pr(`ytr'=1,`yout'=1)" } if "`type'"=="lf" { tempvar zbeta likefun qui g double `zbeta' = (2*(`ytr')-1)*`zbp' tempvar xbj xbjM1 qui g double `xbj'=0 qui g double `xbjM1'=0 forvalues i = 1/`nchoices' { local k = `i'+1 local cutjM1 = `cutmat_out'[1,`i'] local cutj = `cutmat_out'[1,`k'] tempvar m`i' xb`i'j xb`i'jM1 qui g double `m`i'' = `yout'==`i' qui g double `xb`i'jM1' = `m`i''*(`cutjM1'-`xbt') qui g double `xb`i'j' = `m`i''*(`cutj'-`xbt') qui replace `xbj'=`xbj'+`xb`i'j' qui replace `xbjM1' = `xbjM1'+`xb`i'jM1' } qui g double `likefun' = . mata: mata_treatoprobitsim_predict_lf("`likefun'", "`zbp'", "`xbj'", "`xbjM1'", "`lambda'", /// _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`ytr'" ) qui gen `vtyp' `varn' = `likefun' if `touse' label variable `varn' "Likelihood Contribution" } forv i = 2/`nchoices' { if "`type'"=="p1`i'" { local trtv=1 local outv=`i' tempvar vtemp qui g double `vtemp'=. //di `cut_hout' qui generate `vtyp' `varn' = . if `touse' mata: mata_treatoprobitsim_predict("`vtemp'", "`zbp'", "`xbt'", "`cutmat_out'", "`lambda'", /// _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`trtv'", "`outv'" ) qui replace `varn'=`vtemp' if `touse' label variable `varn' "Pr(`ytr'=1,`yout'=`i')" } } forv i = 1/`nchoices' { if "`type'"=="p0`i'" { local trtv=0 local outv=`i' tempvar vtemp qui g double `vtemp'=. qui generate `vtyp' `varn' = . if `touse' mata: mata_treatoprobitsim_predict("`vtemp'", "`zbp'", "`xbt'", "`cutmat_out'", "`lambda'", /// _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`trtv'", "`outv'" ) qui replace `varn'=`vtemp' if `touse' label variable `varn' "Pr(`ytr'=0,`yout'=`i')" } } forv i=1/`nchoices' { if "`type'"=="te`i'" { local outv = `i' tempvar atexb qui g double `atexb' = . if `touse' mata: mata_treatoprobitsim_predict_te("`atexb'", "`zbp'", "`xbtone'", "`xbtzero'", /// "`cutmat_out'", _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`outv'", "`lambda'" ) qui generate `vtype' `varn'= `atexb' if `touse' label variable `varn' "Marginal Effect of `ytr' on `yout1' outcome `i'" } /*if "`type'"=="sete`i'" { tempname W t sehat mat `W' = e(V) local `t'=`cutpts' local xvars `e(xvarout)' local tvars `e(treatrhs)' local wcount: word count `tvars' local tcount = `wcount'+1 local outv = `i' tempvar sehat qui g double `sehat' = . mata: treatoprobitsim_predict_sete("`sehat'", "`xbtone'", "`cutmat_out'", "`xbtzero'", /// _treatoprobit_rnd, _treatoprobit_S, "`cutpts'", "`outv'", "`W'", "`xvars'", /// "`tcount'" ) qui generate `vtype' `varn'= `sehat' if `touse' label variable `varn' "Standard Error of the Marginal Effect of `ytr' on `yout1' outcome `i'" }*/ if "`type'"=="sete`i'" { local outv = `i' tempname beta vcv mat `beta' = e(b) mat `vcv' = e(V) local xvars `e(xvarout)' local xbeta foreach x of local xvars { local xbeta "`xbeta' beta_`x'" } local treatx `e(treatrhs)' local trtx foreach x of local treatx { local trtx "`trtx' btreat_`x'" } local trtx "`trtx' bconstant" local loading lambda local cuts forv i = 1/`cutpts' { local cuts "`cuts' betacut_`i'" } local simnames `trtx' `xbeta' `loading' `cuts' preserve capture drop _all scalar betaC = colsof(`beta') qui drawnorm `simnames', n(`nsim') means(`beta') cov(`vcv') mata: simbetas = st_data(.,"`simnames'") restore tempname constant qui g double `constant' = 1 local treatx "`treatx' `constant'" tempvar atexb qui g double `atexb' = . if `touse' mata: treatoprobitsim_predict_sete("`atexb'", _treatoprobit_rnd, _treatoprobit_S, /// simbetas, "`xvars'", "`treatx'", "`touse'", "`outv'", "`nsim'","`cutpts'", "`lambda'" ) qui generate `vtype' `varn'= `atexb' if `touse' label variable `varn' "Standard Error of the Effect of `ytr' on `yout1' outcome `i'" } } forv i=1/`nchoices' { if "`type'"=="tt`i'" { local outv = `i' tempvar zxb1 qui g double `zxb1' = . if `touse' mata: mata_treatoprobitsim_predict_tt("`zxb1'", "`zbp'", "`xbtone'", "`xbtzero'", /// "`cutmat_out'","`lambda'", _treatoprobit_rnd, _treatoprobit_S, "`touse'", "`outv'" ) qui generate `vtype' `varn' = (`zxb1')/`psel' if `touse' & `ytr'==1 label variable `varn' "Marginal Effect of `ytr' on `yout1' if treated" } if "`type'"=="sett`i'" { local outv = `i' tempname beta vcv mat `beta' = e(b) mat `vcv' = e(V) local xvars `e(xvarout)' local xbeta foreach x of local xvars { local xbeta "`xbeta' beta_`x'" } local treatx `e(treatrhs)' local trtx foreach x of local treatx { local trtx "`trtx' btreat_`x'" } local trtx "`trtx' bconstant" local loading lambda local cuts forv i = 1/`cutpts' { local cuts "`cuts' betacut_`i'" } local simnames `trtx' `xbeta' `loading' `cuts' preserve capture drop _all scalar betaC = colsof(`beta') qui drawnorm `simnames', n(`nsim') means(`beta') cov(`vcv') mata: simbetas = st_data(.,"`simnames'") restore tempname constant qui g double `constant' = 1 local treatx "`treatx' `constant'" tempvar atexb qui g double `atexb' = . if `touse' mata: treatoprobitsim_predict_sett("`atexb'", _treatoprobit_rnd, _treatoprobit_S, /// simbetas, "`xvars'", "`treatx'", "`touse'", "`outv'", "`nsim'","`cutpts'", "`lambda'" ) qui generate `vtype' `varn'= `atexb' if `touse' & `ytr'==1 label variable `varn' "Standard Error of the Marginal Effect of `ytr' on `yout1' outcome `i' if treated" } } if "`type'" =="ptr" { qui generate `vtyp' `varn' = normal(`zbp') if `touse' label variable `varn' "Probability of selection into `ytr'" } if "`type'" =="xbout" { qui generate `vtyp' `varn' = `xbt' if `touse' label variable `varn' "Outcome index" } end mata: function mata_treatoprobitsim_predict(string scalar newv, string scalar zb, string scalar xb, string scalar cutmat, string /// scalar lmbd, real matrix rnd, real scalar S, string scalar touse, string scalar treat, string scalar out) { st_view(zbeta,.,zb, touse) st_view(xbeta,.,xb, touse) st_view(lam,.,lmbd, touse) cutpts = st_matrix(cutmat) c = cols(cutpts) trt = strtoreal(treat) outr = strtoreal(out) k = 2*trt-1 eta = k:*zbeta etarnd = eta:+rnd cutj = cutpts[1,outr+1] cutjM1 = cutpts[1,outr] xbj = cutj:-xbeta xbjrnd = xbj:-lam:*rnd xbjM1 = cutjM1:-xbeta xbjM1rnd = xbjM1:-lam:*rnd newvar = rowsum((normal(etarnd)):*((normal(xbjrnd)):-(normal(xbjM1rnd)))):/S st_store(.,newv,newvar) } end mata: function mata_treatoprobitsim_predict_lf(string scalar newv, string scalar zb, string scalar xbj, string scalar /// xbjM1, string scalar lmbd, real matrix rnd, real scalar S, string scalar touse, string scalar y) { st_view(y1,.,y) st_view(zbetai,., zb) st_view(xbji,., xbj) st_view(xbjM1i,., xbjM1) st_view(lam,.,lmbd) sign = ((2:*y1):-1) zbetarnd = (sign:*zbetai):+rnd normzb = normal(zbetarnd) xbj = xbji:-lam:*rnd xbjM1 = xbjM1i:-lam:*rnd normxb0r = normal(xbj):-normal(xbjM1) treatoprobitsim = normzb:*normxb0r L = (rowsum(treatoprobitsim)):/S L=rowmax((L , J(rows(L),1,smallestdouble()))) st_store(.,newv,L) } mata: function mata_treatoprobitsim_predict_tt(string scalar newv, string scalar zb, string scalar xb1, string scalar xb0, /// string scalar cutmat, string scalar lmbd, real matrix rnd, real scalar S, string scalar touse, string scalar out) { st_view(zbeta,.,zb, touse) st_view(xbeta1,.,xb1, touse) st_view(xbeta0,.,xb0, touse) st_view(lam,.,lmbd,touse) etarnd =zbeta:+rnd cutpts = st_matrix(cutmat) c = cols(cutpts1) outr = strtoreal(out) cutj = cutpts[1,outr+1] cutjM1 = cutpts[1,outr] xb1j = cutj:-xbeta1 xb0j = cutj:-xbeta0 xb1jrnd = xb1j:-lam:*rnd xb0jrnd = xb0j:-lam:*rnd xb1jM1 = cutjM1:-xbeta1 xb0jM1 = cutjM1:-xbeta0 xb1jM1rnd = xb1jM1:-lam:*rnd xb0jM1rnd = xb0jM1:-lam:*rnd newvar = rowsum((normal(etarnd)):*((normal(xb1jrnd):-normal(xb1jM1rnd)):-(normal(xb0jrnd) /// :-normal(xb0jM1rnd)))):/S st_store(.,newv,newvar) } end mata: function treatoprobitsim_predict_sett(string scalar newv, real matrix rnd, real scalar S, /// real matrix simbetas, string scalar xvar, string scalar trtx, /// string scalar touse, string scalar out, string scalar numsim, /// string scalar cuts, string scalar lam) { st_view(X,.,tokens(xvar),touse) st_view(tX,.,tokens(trtx),touse) st_view(lambda,.,lam,touse) n = rows(X) cts = strtoreal(cuts) outr = strtoreal(out) nsim = strtoreal(numsim) csim = cols(simbetas) rsim = rows(simbetas) tCols = cols(tX) xCols = cols(X) neginf = st_numscalar("c(minfloat)") posinf = st_numscalar("c(maxfloat)") ones = J(n,1,1) zeros = J(n,1,0) lambda_t = lambda:*ones /*indexes*/ cbeg = csim-cts+1 cend = csim lambdaidx = cbeg-1 xend = lambdaidx-1 xbeg = lambdaidx-xCols trend = xbeg-1 sim_te = J(n,nsim,.) for(i=1;i<=nsim;i++){ cutpts = neginf, simbetas[i,cbeg..cend], posinf beta = simbetas[i,xbeg..xend]' cutj = cutpts[1,outr+1] cutjM1 = cutpts[1,outr] x1 = (X[.,1..xCols-1],ones) x0 = (X[.,1..xCols-1],zeros) beta_t = simbetas[i,1..trend]' xbeta1 = x1*beta xbeta0 = x0*beta eta = (tX*beta_t) psel = normal(eta:/(sqrt(2))) etarnd = (eta:+rnd) xb1j = cutj:-xbeta1 xb0j = cutj:-xbeta0 xb1jrnd = xb1j:-lambda_t:*rnd xb0jrnd = xb0j:-lambda_t:*rnd xb1jM1 = cutjM1:-xbeta1 xb0jM1 = cutjM1:-xbeta0 xb1jM1rnd = xb1jM1:-lambda_t:*rnd xb0jM1rnd = xb0jM1:-lambda_t:*rnd sim_te[.,i] = (rowsum((normal(etarnd)):*((normal(xb1jrnd):-normal(xb1jM1rnd)):-(normal(xb0jrnd) /// :-normal(xb0jM1rnd)))):/S):/psel } mu_sim = rowsum(sim_te):/nsim newvar = rowsum(sqrt((sim_te:-mu_sim):^2)):/nsim st_store(.,newv,newvar) } end mata: function mata_treatoprobitsim_predict_te(string scalar newv, string scalar zb, string scalar xb1, string scalar xb0, /// string scalar cutmat, real matrix rnd, real scalar S, string scalar touse, string scalar out, string scalar lambda) { st_view(zbeta,.,zb, touse) st_view(xbeta1,.,xb1, touse) st_view(xbeta0,.,xb0, touse) st_view(lam,.,lambda,touse) cutpts = st_matrix(cutmat) c = cols(cutpts1) outr = strtoreal(out) cutj = cutpts[1,outr+1] cutjM1 = cutpts[1,outr] xb1j = cutj:-xbeta1 xb0j = cutj:-xbeta0 xb1jrnd = xb1j:-lam:*rnd xb0jrnd = xb0j:-lam:*rnd xb1jM1 = cutjM1:-xbeta1 xb0jM1 = cutjM1:-xbeta0 xb1jM1rnd = xb1jM1:-lam:*rnd xb0jM1rnd = xb0jM1:-lam:*rnd newvar = rowsum(((normal(xb1jrnd):-normal(xb1jM1rnd)):-(normal(xb0jrnd) /// :-normal(xb0jM1rnd)))):/S st_store(.,newv,newvar) } end mata: function treatoprobitsim_predict_sete(string scalar newv, real matrix rnd, real scalar S, /// real matrix simbetas, string scalar xvar, string scalar trtx, /// string scalar touse, string scalar out, string scalar numsim, /// string scalar cuts, string scalar lam) { st_view(X,.,tokens(xvar),touse) st_view(tX,.,tokens(trtx),touse) st_view(lambda,.,lam,touse) n = rows(X) cts = strtoreal(cuts) outr = strtoreal(out) nsim = strtoreal(numsim) csim = cols(simbetas) rsim = rows(simbetas) tCols = cols(tX) xCols = cols(X) neginf = st_numscalar("c(minfloat)") posinf = st_numscalar("c(maxfloat)") ones = J(n,1,1) zeros = J(n,1,0) lambda_t = lambda /*indexes*/ cbeg = csim-cts+1 cend = csim lambdaidx = cbeg-1 xend = lambdaidx-1 xbeg = lambdaidx-xCols trend = xbeg-1 sim_te = J(n,nsim,.) for(i=1;i<=nsim;i++){ cutpts = neginf, simbetas[i,cbeg..cend], posinf beta = simbetas[i,xbeg..xend]' cutj = cutpts[1,outr+1] cutjM1 = cutpts[1,outr] x1 = (X[.,1..xCols-1],ones) x0 = (X[.,1..xCols-1],zeros) xbeta1 = x1*beta xbeta0 = x0*beta xb1j = cutj:-xbeta1 xb0j = cutj:-xbeta0 xb1jrnd = xb1j:-lambda_t:*rnd xb0jrnd = xb0j:-lambda_t:*rnd xb1jM1 = cutjM1:-xbeta1 xb0jM1 = cutjM1:-xbeta0 xb1jM1rnd = xb1jM1:-lambda_t:*rnd xb0jM1rnd = xb0jM1:-lambda_t:*rnd sim_te[.,i] = rowsum(((normal(xb1jrnd):-normal(xb1jM1rnd)):-(normal(xb0jrnd) /// :-normal(xb0jM1rnd)))):/S } mu_sim = rowsum(sim_te):/nsim newvar = rowsum(sqrt((sim_te:-mu_sim):^2)):/nsim st_store(.,newv,newvar) } end /* mata: function treatoprobitsim_predict_sete(string scalar newv, string scalar xb_1, string scalar cutmat, /// string scalar xb_0, real matrix rnd, real scalar S, string scalar cuts, string scalar out, /// string scalar Mat, string scalar xvar, string scalar tvec) { n=st_nobs() st_view(X,.,tokens(xvar)) st_view(xbeta1,.,xb_1) st_view(xbeta0,.,xb_0) JMat = st_matrix(Mat) cutptsmat = st_matrix(cutmat) cutpts = strtoreal(cuts) //cols(cutpts1) tcols = strtoreal(tvec) nchoices = cutpts+1 kx = cols(X) outr = strtoreal(out) cutj1 = cutptsmat[1,outr+1] cutjM11 = cutptsmat[1,outr] cutj0 = cutptsmat[1,outr+1] cutjM10 = cutptsmat[1,outr] xb1j = cutj1:-xbeta1 xb0j = cutj0:-xbeta0 xb1jrnd = xb1j:-rnd xb0jrnd = xb0j:-rnd xb1jM1 = cutjM11:-xbeta1 xb0jM1 = cutjM10:-xbeta0 xb1jM1rnd = xb1jM1:-rnd xb0jM1rnd = xb0jM1:-rnd newvar = J(n,1,.) tzero = J(n,tcols,0) kappa1 = J(n,cutpts, 0) delta1 = J(n, kx, 0) kappa0 = J(n, cutpts, 0) delta0 = J(n, kx, 0) rho = J(n,1,0) zeta = tzero, delta0, delta1, rho, kappa0, kappa1 for (i=1;i<=kx;i++) { delta0[.,i] = (rowsum((normalden(xb0jrnd):-normalden(xb0jM1rnd))):/S):*X[.,i] delta1[.,i] = (rowsum((normalden(xb1jrnd):-normalden(xb1jM1rnd))):/S):*X[.,i] } delta = -1:*(delta1:-delta0) if (outr==1){ kappa0[.,1] = (rowsum(normalden(xb0jrnd)):/S) kappa1[.,1] = (rowsum(normalden(xb1jrnd)):/S) kappamat = (kappa1:-kappa0) } if (outr==nchoices){ kappa0[.,cutpts] = -1:*(rowsum(normalden(xb0jM1rnd)):/S) kappa1[.,cutpts] = -1:*(rowsum(normalden(xb1jM1rnd)):/S) kappamat = (kappa1:-kappa0) } if (outr>1 & outr