cap program drop small ciplot mysuest

program small, eclass
	args df N
	mat b=e(b)
	mat V=e(V)*((`N'-1)/`df') // this gives exactly the same standard error as single-equation "reg ..., cluster(...)"
	ereturn post b V, obs(`N') dof(`=min(`df',e(N_clust)-1)') // this may give a wider confidence interval because "reg ..., cluster(...)" seems to always use #(clusters)-1 as dof for t-distribution
	ereturn local cmd "regress"
	end 

program ciplot // this is overkill because I use it only once now, but rewriting it would certainly have created errors
	args c p sigma2_cc sigma2_pp sigma2_cp df
	local error = "invt(`df',0.975)*sqrt((x^2)*(`sigma2_pp') + ((1+x)^2)*(`sigma2_cc') - x*(1+x)*(`sigma2_cp'))"
	twoway (function y=(1+x)*`c' - x*`p' + `error', range(-1 0.5) lcolor(black) lpattern(dash)) ///
	       (function y=(1+x)*`c' - x*`p' - `error', range(-1 0.5) lcolor(black) lpattern(dash)) ///
		   (function y=(1+x)*`c' - x*`p'          , range(-1 0.5) lcolor(black) lpattern(solid)), ///
	       xtitle("assumed elasticity") ytitle("implied residual latent crime") xline(0, lstyle(grid)) yline(0, lstyle(grid)) legend(off) subtitle("crime proxied by overall victimization") ///
		   name(eta_k_plot_victim, replace)
	end

program mysuest, eclass properties(mi) // code from I. Canette & Y. Marchenko, StataCorp
	args m1 m2
	qui `m1'
	local df=e(df_r)
	estimates store est1
	qui `m2'
	estimates store est2
	qui suest est1 est2
	small `df' `=e(N)'
	ereturn local title "Seemingly unrelated estimation"
	end

	
*** 1. ICVS & ICPS (left panel of figure 6 & associated tests)

use victim9 $imputable $notimputable I002A I003 numISO USA using MI_ICVS_dataset.dta, clear
sort numISO I002A
collapse (last) victim9 $imputable $notimputable I002A I003, by(numISO USA) // retain only most recent instead of all obs. (and pw) -- to sync with ICPS, and avoid problems of inconsistent weights betweeen equations
tempfile ICVS
save `ICVS'

use i_ln_pris $imputable $notimputable I002A I003 year numISO USA if !mi(i_ln_pris) & !mi(`notimputable') & year==2005 using paneldata.dta, clear
drop year
merge 1:1 numISO using `ICVS', nogenerate update

foreach var of varlist $imputable { // this generates the missing dummies
	quietly misstable sum `var'
	if r(N_eq_dot)!=. {
		gen md_`var'=mi(`var')
		recode `var' (.=0)
		}
	}

* basic estimates
reg i_ln_pris USA $notimputable $imputable md*
estimates store pri
reg victim9 USA $notimputable $imputable i.I003 i.I002A md* // not quietly -- need to check consistency with main results
estimates store vic
local N=e(N)
local df=e(df_r)

* replace normal standard errors by robust standard errors, adjusted for small sample
suest pri vic, vce(robust)
small `df' `N'
ereturn display
ereturn list

* tests
nlcom (eta: [vic_mean]USA/([pri_mean]USA-[vic_mean]USA)), df(`=e(df_r)') // calculate the implied elasticity of crime w.r.t. incarceration rate (not: punishment per crime)
forvalues eta=-.3(.01)-.2 { // search for eta where lower bound of epsilon_k^USA is zero; includes Abrams (2013) elasticity (-.25)
	lincom (1+`eta')*[vic_mean]USA-`eta'*[pri_mean]USA
	}
lincom 0.26*[vic_mean]USA+0.74*[pri_mean]USA	// testing if Drago et al. (2009) elasticity (-.74) is compatible with epsilon_k^USA==0 when using victimization data
lincom [pri_mean]USA-[vic_mean]USA 				// this calculates epsilon_pi^USA as implied by US victimization data (cf. footnote 34 in section 5.2)

* graph
mat V=e(V)
local c=rownumb(V,"vic_mean:USA")
local p=rownumb(V,"pri_mean:USA")
ciplot `=_b["vic_mean:USA"]' `=_b["pri_mean:USA"]' `=V[`c',`c']' `=V[`p',`p']' `=V[`c',`p']' `df'

estimates drop _all


* 2. IHME & ICPS (right panel of figure 6 & associated tests)

use MIdataset.dta, clear
tempfile estimates
mi estimate, saving(`estimates', replace) imputations(`=1+2*$imputations'/`=3*$imputations'): ///
	mysuest "reg i_ln_pris USA $notimputable $imputable" "reg ln_IHMEviol USA $notimputable $imputable"
	
* tests
mi estimate ([est1_mean]USA-[est2_mean]USA)				using `estimates', nocoef // this calculates epsilon_pi^USA as implied by US homicide data (cf. footnote 34 in section 5.2)
mi estimate (0.75*[est2_mean]USA+0.25*[est1_mean]USA)	using `estimates', nocoef // testing if Abrams (2013) elasticity (-.25) is compatible with epsilon_k^USA==0 when using homicide data

* graph
forvalues x=-1(0.02)0.5 {
	mi estimate ((1+`x')*[est2_mean]USA-`x'*[est1_mean]USA) using `estimates', nocoef
	mat lincom=(nullmat(lincom) \ `x',e(b_Q_mi),e(V_Q_mi),e(df_Q_mi))
	}
matrix colnames lincom = eta epsilon_k sigma2 df
svmat lincom, names(col)
gen ub=epsilon_k+invt(df,0.975)*sqrt(sigma2)
gen lb=epsilon_k-invt(df,0.975)*sqrt(sigma2)
twoway (rline ub lb eta, lcolor(black) lpattern(dash)) ///
	   (line epsilon_k eta, lcolor(black) lpattern(solid)), ///
	   xtitle("assumed elasticity") ytitle("implied residual latent crime") yscale(range(-2(1)2)) xline(0, lstyle(grid)) yline(0, lstyle(grid)) legend(off) subtitle("crime proxied by homicide") ///
	   name(eta_k_plot_homicide, replace)
	   

**************************************************
***** combining both graphs into figure 6 ********
graph combine eta_k_plot_victim eta_k_plot_homicide, title("Confidence intervals for residual US latent crime by crime-punishment elasticity", size(msmall)) ///
	note("Point estimates and 95% confidence intervals from linear combinations of regression coefficients (see formula in section 5.2)." ///
	"The right plot has been drawn using MI estimates. The left plot derives from regressions with dummies indicating missing values.", size(vsmall) placement(s)) ///
	saving(results/figure6.gph, replace)
discard
