capture program drop treatoprobit program define treatoprobit, eclass sortpreserve version 15 *!version 1.0.0 cagregory if replay() { if ("`e(cmd)'" != "treatobprobit") error 301 Replay `0' } else Estimate `0' end program define Estimate, eclass syntax anything(equalok) [if] [in] [pweight iweight fweight], treat(varlist) /// [vce(string) level(integer 95) ] di "Hi!" *parse gettoken depvar2 indvar2: anything gettoken depvar1 indvar1: treat local indvar2 "`indvar2' `depvar1'" *options mlopts mlopts, `options' //ml options if "`vce'"!= "" { gettoken firstvce restvce: vce if "`firstvce'" == "`cluster'" { local clopt "vce(cluster `restvce')" } if "`firstvce'" == "`robust'" { local robustopt "vce(robust)" } } if "`weight'" != "" { tempvar wvar quietly gen double `wvar' `exp' local wgt "[`weight'=`wvar']" local awgt "[aw=`wvar']" } *mark sample marksample touse markout `touse' `depvar1' markout `touse' `depvar2' markout `touse' `indvar1' markout `touse' `indvar2' markout `touse' `cluster' *get startvals qui levelsof(`depvar2') if `touse' local t: word count `r(levels)' global ncut = `t'-1 *treatment eq: probit tempname b_trt startmat ystar di in gr _newline "Estimating Treatment equation" probit `depvar1' `indvar1' `wgt' if `touse', nocoef iter(30) qui predict `ystar', xb mat `b_trt' = e(b) *outcome eq: ordered probit di in gr _newline "Estimating Ordered Outcome Equation" oprobit `depvar2' `indvar2' `wgt' if `touse', nocoef tempname b b_out cuts zstar qui predict `zstar', xb mat `b' = e(b) local j = colsof(`b') mat `b_out' =`b'[1,1..`j'-$ncut] *cutpoints mat `cuts'= `b'[1,(`j'-$ncut)+1..`j'] *rho qui corr `zstar' `ystar' local rho_init = min(r(rho),.5) *starting values mat `startmat' = `b_trt', `b_out', `rho_init', `cuts' local cutpts forvalues i = 1/$ncut { local cutpts "`cutpts' /cut`i'" } qui levelsof(`depvar2') global nchoices: word count `r(levels)' *displaying ancilliary parameters local athrho diparm(atanh_rho) local rho diparm(atanh_rho, tanh label("rho")) local di_cuts forv i = 1/$ncut { local di_cuts "`di_cuts' diparm(cut`i') " } local di_cuts "`di_cuts' diparm(__sep__)" *full model constraint 1 _b[cut1:_cons] = 0 di in gr _newline "Estimating Full Model" ml model lf2 treatoprobit_work ("`depvar1'": `depvar1' = `indvar1') /// ("`depvar2'": `depvar2' = `indvar2', noconstant) /// /atanh_rho `cutpts' if `touse' /// `wgt' /// , /// title("Treatment Effects Ordered Probit Regression") /// init(`startmat', copy) /// search(off) /// maximize /// `clopt' /// `mlopts' /// `robustopt' /// `di_cuts' /// `athrho' /// technique(nr) /// nonrtol /// `rho' *ml init `startmat', copy //ml check *ml maximize, nonrtol difficult tempname lrtest pval ereturn local marginsnotok _ALL ereturn local weight `wvar' ereturn scalar k_aux = $ncut+2 ereturn local cmd "treatoprobit" ereturn local depvar "`depvar1' `depvar2'" ereturn local treatx "`indvar1'" ereturn scalar cuts = $ncut ereturn local predict "treatoprobit_p" ereturn local indvarout `indvar2' treatoprobit_replay , level(95) end program define treatoprobit_replay syntax [,level(integer 95)] ml display, level(`level') neq(2) noomitted tempname lrtest pval qui test _b[atanh_rho:_cons]=0 local `lrtest' = round(r(chi2),.01) local `pval' = round(r(p),.01) di in gr `"{p} Test of independent equations = ``lrtest''. Probability of independent equations ``pval''.{p_end}"' end //end parsing and calling command