/****************************************************************************
* Donohue and Wolfers (2006) 
* Figure 7
* Dezhbakhsh, Rubin, and Shepherd (DRS) Distribution of Estimates Under the
* Null of No Deterrent Effect
****************************************************************************/


clear
cd  "..."
program drop _all
set more off
set mem 500m
set mat 500

/*Most of this file contains the code creating the DRSRandomization dataset, which takes a 
very long time to execute, and thus is commented out*/

/*program define setlags
* Program creates the variables required for models 4-6
	drop if YEAR==76
	drop if COUNTY==51570 & YEAR==80 & MURDO1==0
 	tsset COUNTY YEAR
	* But I'll butcher it slightly
 	egen yr=rank(YEAR), by(COUNTY)
 	tsset COUNTY yr
 	sort COUNTY yr
	* Prepare Table 5 variables
	for Q in any PAMURC PEXECF PEXECP PEXECA PSENTP PSENTF PSENTA: gen Qb=Q
	replace PAMURCb=. if MURDO1==0 & AMUR~=.
	for Q in any PAMURC PEXECF PEXECP PEXECA: for Z in num 1/6: gen QZb=lZ.Qb if Qb==.
	for Z in num 1/6: replace PEXECFZb=0 if PEXECFZb==. & YEAR-Z<77
	for Z in num 1/6: replace PAMURCb=PAMURCZb if PAMURCb==. & MURDO1==0 & AMUR~=.
	for Z in num 1/6: for Q in any PEXECF PEXECP PEXECA: replace Qb=QZb if Qb==.
end

program define firststage
	* First Stage, Table 4
	areg  PAMURC POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PAMURCP2 , xbd

	areg  PSENTP POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PSENTPP2 , xbd
	areg  PSENTF POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PSENTFP2 , xbd
	areg  PSENTA POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC] if (FIPSSTAT==4 & (YEAR==82 | YEAR==83))~=1, absorb(COUNTY) 
	predict PSENTAP2, xbd

	areg  PEXECP POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECPP2, xbd
	areg  PEXECF POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECFP2, xbd
	areg  PEXECA POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECAP2, xbd

	* First stage, Table 5
	areg  PAMURCb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PAMURCP2b , xbd

	areg  PSENTPb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PSENTPP2b , xbd
	areg  PSENTFb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PSENTFP2b , xbd
	areg  PSENTAb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC] if (FIPSSTAT==4 & (YEAR==82 | YEAR==83))~=1, absorb(COUNTY) 
	predict PSENTAP2b, xbd

	areg  PEXECPb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECPP2b, xbd
	areg  PEXECFb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECFP2b, xbd
	areg  PEXECAb POLEXPD JLEXPD RATAGA RATROB NRARATE PRISADM vote1 vote2 vote3 vote4 vote5 vote6 DENSITYM RPCIM RPCUI RPCPI PP1019 PP2029 PPB PPN PPM  IYEAR* [aweight=POPC], absorb(COUNTY) 
	predict double PEXECAP2b, xbd
end


** Randomization tests
use CPcounty5_send.dta, clear
drop if YEAR==76
setlags
firststage
drop if COUNTY==. | YEAR==.
drop if COUNTY==51570 & YEAR==80 /* Drop repeated observations*/
gen sample=RATMUR+RATAGA+RATROB+PAMURCP2b+PSENTPP2b+PP1019+PP2029+PPB+PPN+PPM+NRARATE+RPCPI+RPCIM+RPCUI+DENSITYM+YEAR~=. 
egen obs=sum(sample), by(COUNTY)
drop if obs==0
fillin COUNTY YEAR
drop _fillin

sort COUNTY YEAR
areg RATMUR RATAGA RATROB PAMURCP2b PSENTPP2b PP1019 PP2029 PPB PPN PPM NRARATE RPCPI RPCIM RPCUI DENSITYM IYEAR_78-IYEAR_96 [aweight=POPC], absorb(COUNTY)
predict RATMURhat, xbd /* This is necessary for randomizing residuals*/
predict RATMURres, res
preserve
egen countyid=group(COUNTY)
replace COUNTY=countyid
keep COUNTY YEAR RATMUR POPC RATMURres
compress
sort COUNTY YEAR
save randomization, replace
restore

for any RATMUR RATMURres POPC COUNTY: rename X real_X

program drop _all
program define shuffle
	gen double rand=uniform()
	egen double avrand=mean(rand), by(real_COUNTY)
	egen countyid=rank(avrand), by(YEAR)
	gen int COUNTY=countyid
	sort COUNTY YEAR
	merge COUNTY YEAR using randomization, update replace
	drop _merge rand avrand countyid COUNTY
end

gen n=_n
for any betaA seA betaB seB nA nB: gen X=.
gen RATMUR_rand=.

set seed 1
#delimit ;
for Z in num 1/1000:
	shuffle if Z==Z \ /*Randomize y*/
	areg RATMUR RATAGA RATROB PAMURCP2b PSENTPP2b PEXECPP2b PP1019 PP2029 PPB PPN PPM NRARATE RPCPI RPCIM RPCUI DENSITYM IYEAR_78-IYEAR_96 [aweight=POPC], absorb(real_COUNTY) \ 
	replace betaA=_b[PEXECPP2b] if n==Z \
	replace seA=_se[PEXECPP2b] if n==Z \ 
	replace nA=e(N) if n==Z \ /*Randomize residuals*/
	replace RATMUR_rand=RATMURhat+RATMURres \
	areg RATMUR_rand RATAGA RATROB PAMURCP2b PSENTPP2b PEXECPP2b PP1019 PP2029 PPB PPN PPM NRARATE RPCPI RPCIM RPCUI DENSITYM IYEAR_78-IYEAR_96 [aweight=POPC], absorb(real_COUNTY) \
	replace betaB=_b[PEXECPP2b] if n==Z \
	replace seB=_se[PEXECPP2b] if n==Z \ 
	replace nB=e(N) if n==Z 
;
#delimit cr

for any A B: gen effectX=-(2293*(betaX/320)+1) \ gen seeffectX=sqrt(2293*2293*(seX*seX/(320*320))) \ gen tX=betaX/seX

preserve
keep beta* se* effect* n*
save DRSRandomization, replace*/


use DRSRandomization, clear
#delimit;
kdensity effectA if abs(effectA)<150,
	xline(18.9, lpattern(dash)) 
	/*title("Distribution of Estimates Under the Null of No Deterrent Effect")*/
	xtitle("Estimated Life-Life Tradeoff" "Lives saved by executing one more death row inmate") 
	xlabel(-100(25)100)
	text(.014 60 "DRS (2003) preferred estimate") 
	n(1000)
;
#delimit cr