* How should we handle rho under twostep when the estimate of rho is outside * the range {-1, 1}? * * Bottom line: The best coverage probs for the nulls are for option rhosigma * version 6 global N 50 global rho .8 global selrate .4 global reps 1000 set seed 123456789 cap program drop simit program define simit if "`1'" == "?" { global S_1 p1 pt1 pl1 pfs1 pml1 /* */ p2 pt2 pl2 pfs2 pml2 /* */ p3 pt3 pl3 pfs3 pml3 /* */ pc ptc plc pfsc pmlc /* */ b1 b2 b3 bc /* */ bml1 bml2 bml3 bmlc /* */ rho rhoml exit } /* simulate the data */ drop _all set obs \$N gennorm u v, corr(\$rho) replace u = 2 * u gen double x1 = uniform() gen double x2 = uniform() gen double x3 = x1 + uniform() gen double y = 0.5 + x1 - 2*x2 + 3*x3 + u gen double z1 = x3 + uniform() gen double z2 = z1 + x1 + uniform() gen double zg = -2 - 2*z1 + 2*z2 + v /* this is mean 0 */ sum zg gen double zprb = normprob(zg) gen byte yseen = zprb > (1-\$selrate) noi heckman y x1-x3, sel(yseen = z1 z2) twostep rhoforce test x1 = 1 local p1 = r(p) test x2 = -2 local p2 = r(p) test x3 = 3 local p3 = r(p) test _cons = .5 local pc = r(p) local rho = e(rho) local b1 = _b[x1] local b2 = _b[x2] local b3 = _b[x3] local bc = _b[_cons] heckman y x1-x3, sel(yseen = z1 z2) twostep rhotrunc test x1 = 1 local pt1 = r(p) test x2 = -2 local pt2 = r(p) test x3 = 3 local pt3 = r(p) test _cons = .5 local ptc = r(p) heckman y x1-x3, sel(yseen = z1 z2) twostep rholim test x1 = 1 local pl1 = r(p) test x2 = -2 local pl2 = r(p) test x3 = 3 local pl3 = r(p) test _cons = .5 local plc = r(p) heckman y x1-x3, sel(yseen = z1 z2) twostep rhosigma test x1 = 1 local pfs1 = r(p) test x2 = -2 local pfs2 = r(p) test x3 = 3 local pfs3 = r(p) test _cons = .5 local pfsc = r(p) heckman y x1-x3, sel(yseen = z1 z2) nonrtolerance difficult test x1 = 1 local pml1 = r(p) test x2 = -2 local pml2 = r(p) test x3 = 3 local pml3 = r(p) test _cons = .5 local pmlc = r(p) local rhoml = e(rho) local bml1 = _b[x1] local bml2 = _b[x2] local bml3 = _b[x3] local bmlc = _b[_cons] post `1' `p1' `pt1' `pl1' `pfs1' `pml1' /* */ `p2' `pt2' `pl2' `pfs2' `pml2' /* */ `p3' `pt3' `pl3' `pfs3' `pml3' /* */ `pc' `ptc' `plc' `pfsc' `pmlc' /* */ (`b1') (`b2') (`b3') (`bc') /* */ (`bml1') (`bml2') (`bml3') (`bmlc') /* */ (`rho') (`rhoml') end simul simit, reps(\$reps) dots * How many times was the covariance matrix non-pd with the -rhoforce- method? count if p1 == . local nonpd = r(N) * How many times is the consistent estimate of rho outside the [-1, 1] * interval? count if rho < -1 | rho > 1 local outside = r(N) * How does the empirical variance of the ML estimate compare with that of the * twostep estimate sum b1 bml1 b2 bml2 b3 bml3 bc bmlc * Coverage probabilities * Note for above, p# is the p-value with no adjustment, even if rho * outside [-1,1], this is what -heckman- did * previously * pt# is with rho truncated to -1 or 1 if outside range * pfs# is with rho truncated and sigma adjusted to be * consistent with the new rho. * pl# is with rho truncated only where used as a weight, * (see expression for D in computing V_twostep in * [R] heckman, Methods and Formulas) * pmlX is with estimation by ML cap prog drop Report program define Report args pval local prelist p pt pl pfs pml local names rhoforce rhotrunc rholimited rhosigma MLE local blist 1 2 3 c local trate = 100 * `pval' di _n "Rejection rates for `trate'% tests against true nulls" di "(nominal rates are `pval') di "Rho method" _col(12) "| B_1 B_2 B_3 B_cons" di _dup(11) "-" "+" _dup(36) "-" tokenize `prelist' local i 1 while "``i''" != "" { local name : word `i' of `names' di "`name'" _col(12) "|" _c local j 1 local Nj : word count `blist' while `j' <= `Nj' { local coef : word `j' of `blist' qui count if ``i''`coef' < `pval' di %9.3f r(N)/\$reps _c local j = `j' + 1 } di "" local i = `i' + 1 } end cap prog drop Header program define Header args outside nonpd di _n di " Data generating process: N = \$N " di " rho = \$rho" di " selection rate = \$selrate" di " repetitions = \$reps" di _n "Number of times estimate of rho outside [-1,1] = `outside'" di "Number of times rhoforce produced non-PD VCE = `nonpd'" end Header `outside' `nonpd' Report .01 Report .05 Report .10