capture estimates drop

***********************************************************
***********************************************************
************ Figure 4: raw data ***************************

use numISO year BJShom ln_BJShom homrate ln_homrate numICD victim9 I002A I003 UNCTSprisrate ln_UNCTSprisrate if !mi(BJShom) | !mi(homrate) | !mi(victim9) using paneldata.dta, clear

* generate demeaned values
gen ln_victim9=ln(victim9)
quietly xtreg ln_homrate i.numICD i.year, fe
predict demeaned_hom, ue
quietly xtreg ln_victim9 i.I002A i.I003, fe
predict demeaned_victim, ue
quietly xtreg ln_UNCTSprisrate i.year, fe
predict demeaned_pris, ue
label var demeaned_hom "Homicides per 100,000 (WHO)"
label var demeaned_victim "Victimization rate (ICVS)"
label var demeaned_pris "Incarceration per 100,000 (UNCTS)"

* create US dataset only
keep if numISO=="USA":numISO3code
drop numISO
merge 1:1 year using USprisoners.dta, nogenerate keepusing(incarcerationrate) // this is US inc. rate from BJS NCJ 230113, 243936, 173414 (in that order of preference)
drop if year>2010

* levels
replace BJShom=BJShom/10
replace incarceration=incarceration/1000
sum year if !mi(incarceration) | !mi(BJShom) | !mi(victim9)
twoway (line incarceration year, lpattern(solid) lcolor(blue) mcolor(blue)) ///
	(line BJShom year, lpattern(longdash) lcolor(red) mcolor(red)) ///
	(line victim9 year, lpattern(shortdash) lcolor(green) mcolor(green)), ///
	yscale(range(0)) ylabel(0(0.5)1) ytitle("") ///
	xlabel(`=10*floor(`r(min)'/10)'(10)`=10*ceil(`r(max)'/10)') ///
	subtitle("Levels") ///
	legend(cols(3) label(1 "Incarceration per 100") label(2 "Homicides per 10,000") label(3 "Victimization per 1") size(vsmall) symxsize(small))  ///
	scheme(s1color) name(figure4a, replace)

* demeaned
sum year if !mi(demeaned_pris) | !mi(demeaned_hom) | !mi(demeaned_victim)
twoway (lpoly demeaned_pris year, lpattern(solid) lcolor(blue) mcolor(blue)) ///
	(lpoly demeaned_hom year, lpattern(longdash) lcolor(red) mcolor(red)) ///
	(lpoly demeaned_victim year, lpattern(shortdash) lcolor(green) mcolor(green)), ///
	yscale(range(0)) ylabel(0(.5)1.5) ///
	xtitle("year") xlabel(`=10*floor(`r(min)'/10)'(10)`=10*ceil(`r(max)'/10)') ///
	subtitle("Logs net of annual constant-sample world mean") legend(off) ///
	scheme(s1color) name(figure4b, replace)
	
* combine
grc1leg figure4a figure4b, ///
	scheme(s1color) cols(1) ysize(10.5) xsize(7.5) title("US crime and punishment over time") legendfrom(figure4a) ///
	note("Data sources: Upper graph: BJS, UCR, and ICVS. Lower graph: UNCTS, WHO, and ICVS." "The lower graph smooths US fixed effects plus annual residuals from a country fixed effect regression of the respective variable," ///
	"log-transformed, on year dummies and, in the case of WHO and ICVS, ICD version and sweep dummies, respectively.", size(tiny))
graph display, ysize(10.5) xsize(7.5)
graph save results/figure4, replace
discard

***********************************************************
***********************************************************
************ Figure 5: residuals **************************

*********************************************
****** 5a: Extensions of main regressions ***

tempfile graphs // this is a bin for accumulating the results. It also stores the raw US data, just in case a comparison is desired.

*** ICVS
use MI_ICVS_dataset.dta, clear
set emptycells drop

mi estimate: reg victim9 ibn.year#1.USA $notimputable $imputable i.I003 i.I002A [pw=1/freq], cluster(numISO)
mi extract 0
keep if USA==1
keep year victim9
save `graphs', replace
clear

matrix B=e(b_mi)
matrix V=vecdiag(e(V_mi))
matrix R=(1989,1992,1996,2000,2004\B[1,1..5]\V[1,1..5])'
matrix colnames R = year ICVSresidual ICVSresidual_se
svmat R, names(col)
replace ICVSresidual_se=sqrt(ICVSresidual_se)
merge 1:1 year using `graphs', nogenerate
save `graphs', replace


*** IHME and ICPS
tempfile IHME_ICPS
use year numISO ln_IHMEviolentdeaths i_ln_pris USA $notimputable $imputable if (!mi(ln_IHME) | !mi(i_ln_pris)) & !mi(legor_uk, fedccp, rgdp) using paneldata.dta, clear

* making country coverage consistent with main regression table for year 2005
foreach compcontrol of global notimputable {
drop if mi(`compcontrol')
}

* adding missing dummies
local imputed ""
foreach var of varlist $imputable {
	quietly misstable sum `var'
	if r(N_eq_dot)!=. local imputed "`imputed' `var'"
	}
foreach indepvar in `imputed' {
gen md_`indepvar'=mi(`indepvar')
recode `indepvar' (.=0)
}

save `IHME_ICPS', replace

* merge raw data into ICVS data
keep if USA==1
keep year ln_IHME i_ln_pris
merge 1:1 year using `graphs', nogenerate
save `graphs', replace

* IHME residuals
use `IHME_ICPS', clear
reg ln_IHME ibn.year#1.USA i.year $notimputable $imputable md_*, cluster(numISO) // since IHME only has 3 data points (1990, 2005, 2010), year dummies are equivalent to a quadratic in year
clear

matrix B=e(b)
matrix V=vecdiag(e(V))
matrix R=(1990,2005,2010\B[1,1..3]\V[1,1..3])'
matrix colnames R = year IHMEresidual IHMEresidual_se
svmat R, names(col)
replace IHMEresidual_se=sqrt(IHMEresidual_se)
merge 1:1 year using `graphs', nogenerate
save `graphs', replace

* ICPS
use `IHME_ICPS', clear
gen year2=year^2
reg i_ln_pris ibn.year#1.USA year year2 $notimputable $imputable md_*, cluster(numISO)
clear

matrix R=(e(b)\vecdiag(e(V)))'
matrix colnames R = ICPSresidual ICPSresidual_se
local names : rowfullnames R
svmat R, names(col)
gen year=""
forvalues i=1/`: word count `names'' {
  replace year=`"`: word `i' of `names''"' in `i'
}
drop if !regexm(year,".year#1.USA")
replace year=substr(year,1,4)
destring year, replace

replace ICPSresidual_se=sqrt(ICPSresidual_se)

merge 1:1 year using `graphs', nogenerate
save `graphs', replace


*** Graph

sort year

label var IHMEresidual "residual homicide rate (IHME-GBD)"
label var ICPSresidual "residual incarceration rate (ICPS)"
label var ICVSresidual "residual victimization rate (ICVS)"

gen ICPShigh=ICPSresidual+2*ICPSresidual_se	// these are the confidence bounds
gen ICPSlow=ICPSresidual-2*ICPSresidual_se
gen IHMEhigh= IHMEresidual+2* IHMEresidual_se
gen IHMElow= IHMEresidual-2* IHMEresidual_se
gen ICVShigh= ICVSresidual+2* ICVSresidual_se
gen ICVSlow= ICVSresidual-2* ICVSresidual_se
	
twoway (connected ICPSresidual IHMEresidual ICVSresidual year, lpattern(solid longdash shortdash) lcolor(blue red green) mcolor(blue red green)) ///
	(rspike ICPShigh ICPSlow year, lcolor(blue)) (rspike IHMEhigh IHMElow year, lcolor(red)) (rspike ICVShigh ICVSlow year, lcolor(green)), ///
	ylabel(0(1)2) legend(order(1 2 3) label(1 "residual incarceration") label(2 "residual homicide") label(3 "residual victimization")) ///
	subtitle("US dummies from panel regressions using ICPS/IHME/ICVS data") scheme(s1color) name(figure5a, replace)


************************************************************
****** 5b: longer time series with UNCTS and WHOSIS data ***

use year numISO ln_UNCTSprisrate ln_homrate numICD USA $historical $notimputable if (!mi(ln_UNCTSprisrate) | !mi(ln_homrate)) using paneldata.dta, clear
gen year2 = year^2
foreach indepvar in $historical $notimputable { // missing dummies (NB: no democracy and propdem data before 1975, which as dummies are usually not imputed; also no free data before 1972, so don't go further back)
	gen md_`indepvar'=mi(`indepvar')
	recode `indepvar' (.=0)
}

preserve
reg ln_homrate ibn.year#1.USA year* $historical $notimputable md_* i.numICD if year>=1972, cluster(numISO)
parmest, norestore
keep parm estimate min max
drop if substr(reverse(parm),1,3)!="ASU"
drop if substr(parm,5,1)=="o"
gen year=substr(parm,1,4)
drop parm
destring year, replace
rename (estimate min max) hom_=
tempfile UNresults
save `UNresults', replace

restore
reg ln_UNCTSprisrate ibn.year#1.USA year* $historical $notimputable md_* if year>=1972, cluster(numISO)
parmest, norestore
keep parm estimate min max
drop if substr(reverse(parm),1,3)!="ASU"
drop if substr(parm,5,1)=="o"
gen year=substr(parm,1,4)
drop parm
destring year, replace
rename (estimate min max) pris_=

merge 1:1 year using `UNresults', nogenerate
sort year
gen hyear=year+0.5 // prevents the confidence spikes from overlapping

twoway (line pris_estimate hom_estimate year, lpattern(solid longdash) lcolor(blue red) mcolor(blue red)) ///
	(rspike pris_max pris_min hyear, lcolor(blue)) (rspike hom_max hom_min year, lcolor(red)), ///
	ylabel(0(1)2)  legend(off) subtitle("US dummies from panel regressions using UNCTS/WHO data") ///
	scheme(s1color) name(figure5b, replace)


	
*********************************************************
******* creating figure 5 by combining graphs ***********

grc1leg figure5a figure5b, ///
	scheme(s1color) cols(1) ysize(10.5) xsize(7.5) title("Log-residual US crime and punishment over time", size(medium)) legendfrom(figure5a) ///
	note("Estimated US dummies and 95% confidence intervals from panel regressions with missing value dummies (OLS+) and a quadratic trend." ///
	"Standard errors are clustered at the country level. All scales are in logs. The independent variables are the same as in tables 3 and 4 except" ///
	"that graph 2 omits labor laws, unemployment, and teen birth rates, and uses Gini data from UTIP instead of the OECD or World Bank." ///
	, size(tiny))
graph display, ysize(10.5) xsize(7.5)
graph save results/figure5, replace
discard
