// Code to generate results in Roodman Lancet letter of September 15, 2002 mata mata set matafavor speed // procedure to run an mi estimate, then rescue some test statistics by aggregating over imputations cap program drop MI program define MI, eclass tempfile estimates mi estimate (longrun: _b[govdisgdp_up_net]/(1-_b[L.ghegdp_s])), post saving(`estimates'): `0' // estimate est store miEstimate preserve clear set obs $Nimputations gen double pcaR2 = . gen double hansenp = . forvalues i=1/$Nimputations { est use `estimates', number(`i') qui replace pcaR2 = e(pcaR2) in `i' qui replace hansenp = e(hansenp) in `i' } est restore miEstimate est drop miEstimate collapse pcaR2 hansenpmean=hansenp (p50) hansenpmedian=hansenp, fast foreach stat of varlist * { ereturn scalar `stat' = `stat' } restore mat b = e(b_Q_mi) ereturn scalar longrun = b[1,1] mi testtransform longrun ereturn scalar longrunp = r(p) end postfile PCA str10(MI source collapse) hansenp j alpha alpha_CI beta beta_CI longrunbeta longrunbetap components pcaR2 hansenpmedian hansenpmean /// using PCA, replace foreach MI in "" MI { // MI and non-MI runs use "Lu et al..dta", clear ren mghegdp_who_s ghegdp_who_s ren mghegdp_john_s ghegdp_imf_s ren pridisgdp_up_net pridisbgdp_up_net ren mgdppc_usd06_imf gdppc_usd_imf ren mggegdpwb ggegdp_wb encode country, gen(countrynumber) tsset countrynumber year drop if country=="Angola" | country=="Eritrea" // Lu et al. Table 2 xtabond2 ghegdp_who_s l.ghegdp_who_s govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, /// gmm(l.ghegdp_who_s l.govdisgdp_up_net) iv(govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, eq(diff)) robust h(2) di "Long-run effect: " _b[govdisgdp_up_net]/(1-_b[L.ghegdp_who_s]) testnl _b[govdisgdp_up_net]/(1-_b[L.ghegdp_who_s]) = 0 xtabond2 ghegdp_imf_s l.ghegdp_imf_s govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, /// gmm(l.ghegdp_imf_s l.govdisgdp_up_net) iv(govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, eq(diff)) robust h(2) di "Long-run effect: " _b[govdisgdp_up_net]/(1-_b[L.ghegdp_imf_s]) testnl _b[govdisgdp_up_net]/(1-_b[L.ghegdp_imf_s]) = 0 gen byte sample = e(sample) if "`MI'"!="" { // Put multiply imputed file from Lu in form needed for Stata's -mi- tempfile temp keep countrynumber year sample save `temp', replace use "Lu mi data set 2012.dta", clear keep _mj - countrynumber joinby countrynumber year using `temp', unmatched(master) recode sample (. = 0) save `temp', replace collapse (first) _mj hiv_prevalence-sample (sd) sd_who=ghegdp_who_s sd_imf=ghegdp_imf_s, by(_mi country iso3 year) order _mi _mj replace ghegdp_who_s =. if sd_who replace ghegdp_imf_s =. if sd_imf drop sd_* replace _mj=0 append using `temp' save `temp', replace mi import ice, automatic drop if country=="Angola" | country=="Eritrea" mi set global Nimputations = r(M) mi xtset countrynumber year sort countrynumber year // MI versions of published regressions foreach source in who imf { ren ghegdp_`source'_s ghegdp_s // MI version of published regressions MI xtabond2 ghegdp_s l.ghegdp_s govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence if sample, /// gmm(l.ghegdp_s l.govdisgdp_up_net) iv(govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, eq(diff)) robust h(2) // fixing inadvertent treatment of health aid as exogenous--basis for "Fixes to these two problems..." paragraph in letter MI xtabond2 ghegdp_s l.ghegdp_s govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence if sample, /// gmm(l.ghegdp_s l.govdisgdp_up_net) iv( pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, eq(diff)) robust h(2) di "Long-run effect: " e(longrun) ren ghegdp_s ghegdp_`source'_s } } // generate data for figure (plus non-MI and collapsed variants) foreach collapse in "" collapse { forvalues components=2/`=cond("`collapse'"=="", 130, 22)' { foreach source in who imf { ren ghegdp_`source'_s ghegdp_s `MI' xtabond2 ghegdp_s l.ghegdp_s govdisgdp_up_net pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence if sample, /// gmm(l.ghegdp_s l.govdisgdp_up_net, `collapse') /// iv(pridisbgdp_up_net drdisgdp gdppc_usd_imf ggegdp_wb hiv_prevalence, eq(diff)) /// robust pca components(`components') small if "`MI'" == "" { scalar longrun = _b[govdisgdp_up_net]/(1-_b[L.ghegdp_s]) testnl _b[govdisgdp_up_net]/(1-_b[L.ghegdp_s]) = 0 scalar longrunp = r(p) mat beta_df = e(df_r) mat alpha_df = e(df_r) } else { scalar longrun = e(longrun) scalar longrunp = e(longrunp) mat df = e(df_mi) mat beta_df = df[1, "govdisgdp_up_net"] // HT Joe Dieleman mat alpha_df = df[1, "l.ghegdp_s"] } scalar beta_CI = _se[govdisgdp_up_net] * invttail(beta_df [1,1], .05/2) // half width of 95% CI. scalar alpha_CI = _se[l.ghegdp_s ] * invttail(alpha_df[1,1], .05/2) post PCA ("`MI'") ("`source'") ("`collapse'") (0`e(hansenp)') (`e(j)') (_b[l.ghegdp_s]) (alpha_CI) (_b[govdisgdp_up_net]) (beta_CI) /// (longrun) (longrunp) (`components') (`e(pcaR2)') (0`e(hansenpmedian)') (0`e(hansenpmean)') ren ghegdp_s ghegdp_`source'_s } } } } postclose PCA use PCA, clear gen beta_hi = beta + beta_CI gen beta_lo = beta - beta_CI gen alpha_hi = alpha + alpha_CI gen alpha_lo = alpha - alpha_CI replace components = components*1.02 if source=="imf" * Lancet letter graph, with non-MI and collapsed variants foreach MI in "" MI { foreach collapse in collapse "" { sum beta_hi if components>=4 & MI=="`MI'" & collapse=="`collapse'" twoway rcap beta_hi beta_lo components if source=="imf", yaxis(2) lcolor(maroon) lwidth(medthin) || /// rcap beta_hi beta_lo components if source=="who" , yaxis(2) lcolor(navy) lwidth(medthin) || /// line hansenp`=cond("`MI'"=="","","median")' components if source=="who", lcolor(navy) lwidth(medthin) || /// line hansenp`=cond("`MI'"=="","","median")' components if source=="imf", lcolor(maroon) lwidth(medthin) || /// if components>=4 & MI=="`MI'" & collapse=="`collapse'", xscale(log) /// legend(order(3 4) margin(zero) bmargin(zero) region(style(none)) label(3 "WHO") label(4 "IMF")) /// xtitle("Number of retained principal components of instrument set") /// ytitle("", axis(1)) ytitle("", axis(2)) xlabel(5 10 20 50 100) /// yscale(range(0 `=(`r(max)'+3)/4')) /// ylabel(0 .25 .5 .75 1, grid angle(horiz)) /// ylabel(-3 "–3" -2 "–2" -1 "–1" 0 1 "+1", axis(2)) ylabel(, axis(2) angle(horiz)) /// scheme(s1color) plotregion(style(none) margin(zero)) graphregion(margin(zero)) xsize(7.5) ysize(5) /// name(gr`MI'`collapse', replace) } }