/***********************************************************************************
***********************************************************************************

Replication of empirical results in Chapter 8: Concept Validyt and Measurement

of 

Eric Neumayer and Thomas Plümper: Robustness Tests for Quantitative Research, 
   Cambridge University Press 2017

***********************************************************************************
***********************************************************************************/

* Stata version 12 or higher needed
version 12

***********************************************************************************
***********************************************************************************

* change relative path to the directory where the data are located               
local DIR =  "C:\Users\Eric\Documents\Book\Robustness\"
cd "`DIR'"

***********************************************************************************
***********************************************************************************

* Install user-written ado-files
capture net install outreg2, from(http://fmwww.bc.edu/RePEc/bocode/o)
capture net install mimrgns, from(http://fmwww.bc.edu/RePEc/bocode/m)


***********************************************************************************
***********************************************************************************


*** Randomized Components Test

scalar drop _all
matrix drop _all
*reference model
use "Chapter_8_measurement_data_2.dta", clear
set level 95
logit onset i.warl  gdpenl  lpopl1  lmtnest i.ncontig i.Oil i.nwstate  i.instab ethfrac relfrac i.deml i.anocl, level(90)
margins, dydx(_all) post
scalar anoc_lpolity2 =  _b[1.anocl]
scalar low_anoc_lpolity2 =  _b[1.anocl] - invnormal(1-(100-c(level))/200)*_se[1.anocl]
scalar up_anoc_lpolity2 = _b[1.anocl] + invnormal(1-(100-c(level))/200)*_se[1.anocl]

capture  program drop polity2sim_dum 
use "Chapter_8_measurement_data_3.dta", clear
matrix Weights_default = J(1,5,1)
matrix Range = (Weights_default-Weights_default/2 \ Weights_default+Weights_default/2)

program polity2sim_dum
version 12
drop _all
matrix Weights = ((Range[2,1]-Range[1,1])*runiform()+Range[1,1] \ (Range[2,2]-Range[1,2])*runiform()+Range[1,2] \ (Range[2,3]-Range[1,3])*runiform()+Range[1,3] \ (Range[2,4]-Range[1,4])*runiform()+Range[1,4] \ (Range[2,5]-Range[1,5])*runiform()+Range[1,5])'
use "Chapter_8_measurement_data_2.dta", clear
tsset ccode year
gen polity2lsim = lxrcomppolity*Weights[1,1]+lxropenpolity*Weights[1,2]+lxconstpolity*Weights[1,3]+lparregpolity*Weights[1,4]+lparcomppolity*Weights[1,5]
capture dropvars polity2lsim_anoc polity2lsim_dem
gen polity2lsim_anoc=0
replace polity2lsim_anoc=1 if polity2lsim>=-5 & polity2lsim<=5
replace polity2lsim_anoc=. if polity2lsim==.
gen polity2lsim_dem=0
replace polity2lsim_dem=1 if polity2lsim>5 & polity2lsim<.
replace polity2lsim_dem=. if polity2lsim==.
logit onset i.warl  gdpenl  lpopl1  lmtnest i.ncontig i.Oil i.nwstate  i.instab ethfrac relfrac i.polity2lsim_anoc i.polity2lsim_dem, level(90)
estat ic
matrix IC = r(S)
scalar BIC = IC[1,6]
margins, dydx(_all) post
scalar anoc_se_lpolity2sim = _se[1.polity2lsim_anoc]
scalar dem_se_lpolity2sim = _se[1.polity2lsim_dem]
scalar anoc_lpolity2sim = _b[1.polity2lsim_anoc]
scalar dem_lpolity2sim =  _b[1.polity2lsim_dem]
scalar low_anoc_lpolity2sim = _b[1.polity2lsim_anoc] - invnormal(1-(100-c(level))/200)*_se[1.polity2lsim_anoc]
scalar up_anoc_lpolity2sim = _b[1.polity2lsim_anoc] + invnormal(1-(100-c(level))/200)*_se[1.polity2lsim_anoc]
scalar low_dem_lpolity2sim =  _b[1.polity2lsim_dem] - invnormal(1-(100-c(level))/200)*_se[1.polity2lsim_dem]
scalar up_dem_lpolity2sim = _b[1.polity2lsim_dem] + invnormal(1-(100-c(level))/200)*_se[1.polity2lsim_dem]
scalar  lb1=low_anoc_lpolity2/anoc_se_lpolity2sim
scalar  ub1=up_anoc_lpolity2/anoc_se_lpolity2sim
scalar  lb2=lb1-(anoc_lpolity2sim/anoc_se_lpolity2sim)
scalar  ub2=ub1-(anoc_lpolity2sim/anoc_se_lpolity2sim)
scalar  ubd=normal(ub2)
scalar  lbd=normal(lb2)
scalar  robustness = ubd-lbd
scalar z_anoc_lpolity2 = normal(-anoc_lpolity2sim/anoc_se_lpolity2sim)
end

simulate rho = robustness  anoc_lpolity2sim = anoc_lpolity2sim low_anoc_lpolity2sim = low_anoc_lpolity2sim up_anoc_lpolity2sim = up_anoc_lpolity2sim dem_lpolity2sim = dem_lpolity2sim low_dem_lpolity2sim = low_dem_lpolity2sim up_dem_lpolity2sim = up_dem_lpolity2sim anoc_se_lpolity2sim = anoc_se_lpolity2sim dem_se_lpolity2sim = dem_se_lpolity2sim bic = BIC zero_anoc_lpolity2 = z_anoc_lpolity2, reps(1000) seed(1234): polity2sim_dum
gen countzero_dem_lpolity2 = (low_dem_lpolity2sim < 0) & (up_dem_lpolity2sim > 0)
gen countzero_anoc_lpolity2 = (low_anoc_lpolity2sim < 0) & (up_anoc_lpolity2sim > 0)
describe

** without weighting
summarize
_pctile rho, p(2.5, 97.5)
matrix range_robustness = (r(r1),r(r2))
matrix rownames range_robustness = rho
matrix colnames range_robustness = lower_025 upper_975
quietly summarize rho
display "Estimated degree of effect robustness is: " r(mean)
display  "95% of unweighted simulated effect robustness ranges between: "
matrix list range_robustness


***********************************************************************************
***********************************************************************************



*** Principal Component Jitter Test

scalar drop _all
matrix drop _all
*reference model
use "Chapter_8_measurement_data_2.dta", clear
set level 95
logit onset i.warl  gdpenl  lpopl1  lmtnest i.ncontig i.Oil i.nwstate  i.instab ethfrac relfrac i.deml i.anocl, level(90)
margins, dydx(_all) post
scalar anoc_lpolity2 =  _b[1.anocl]
scalar low_anoc_lpolity2 =  _b[1.anocl] - invnormal(1-(100-c(level))/200)*_se[1.anocl]
scalar up_anoc_lpolity2 = _b[1.anocl] + invnormal(1-(100-c(level))/200)*_se[1.anocl]
use "Chapter_8_measurement_data_3.dta", clear
pca xrcomppolity xropenpolity xconstpolity parregpolity parcomppolity, components (1) cov vce(normal)
matrix Weights = e(L)
scalar theoreticmin = (-2*abs(Weights[1,1])-1*abs(Weights[2,1])-3*abs(Weights[3,1])-2*abs(Weights[4,1])-2*abs(Weights[5,1]))
scalar theoreticmax = (2*abs(Weights[1,1])+1*abs(Weights[2,1])+4*abs(Weights[3,1])+0*abs(Weights[4,1])+3*abs(Weights[5,1]))
matrix Range = (Weights-Weights/2 , Weights+Weights/2)'
capture  program drop polity2jitter_dum

program polity2jitter_dum
version 12
drop _all
matrix Weightsjitter = ((Range[2,1]-Range[1,1])*runiform()+Range[1,1] \ (Range[2,2]-Range[1,2])*runiform()+Range[1,2] \ (Range[2,3]-Range[1,3])*runiform()+Range[1,3] \ (Range[2,4]-Range[1,4])*runiform()+Range[1,4] \ (Range[2,5]-Range[1,5])*runiform()+Range[1,5])'
use "Chapter_8_measurement_data_2.dta", clear
tsset ccode year
gen score = lxrcomppolity*Weightsjitter[1,1]+lxropenpolity*Weightsjitter[1,2]+lxconstpolity*Weightsjitter[1,3]+lparregpolity*Weightsjitter[1,4]+lparcomppolity*Weightsjitter[1,5]
gen  polity2ljitter = (score - theoreticmin)/(theoreticmax-theoreticmin)*20-10
gen polity2ljitter_anoc=0
replace polity2ljitter_anoc=1 if polity2ljitter>=-5 & polity2ljitter<=5
replace polity2ljitter_anoc=. if polity2ljitter==.
gen polity2ljitter_dem=0
replace polity2ljitter_dem=1 if polity2ljitter>5 & polity2ljitter<.
replace polity2ljitter_dem=. if polity2ljitter==.
logit onset i.warl  gdpenl  lpopl1  lmtnest i.ncontig i.Oil i.nwstate  i.instab ethfrac relfrac i.polity2ljitter_anoc i.polity2ljitter_dem
estat ic
matrix IC = r(S)
scalar BIC = IC[1,6]
margins, dydx(_all) post
scalar anoc_se_lpolity2jitter = _se[1.polity2ljitter_anoc]
scalar dem_se_lpolity2jitter = _se[1.polity2ljitter_dem]
scalar anoc_lpolity2jitter = _b[1.polity2ljitter_anoc]
scalar dem_lpolity2jitter =  _b[1.polity2ljitter_dem]
set level 90
scalar low_anoc_lpolity2jitter = _b[1.polity2ljitter_anoc] - invnormal(1-(100-c(level))/200)*_se[1.polity2ljitter_anoc]
scalar up_anoc_lpolity2jitter = _b[1.polity2ljitter_anoc] + invnormal(1-(100-c(level))/200)*_se[1.polity2ljitter_anoc]
scalar low_dem_lpolity2jitter =  _b[1.polity2ljitter_dem] - invnormal(1-(100-c(level))/200)*_se[1.polity2ljitter_dem]
scalar up_dem_lpolity2jitter = _b[1.polity2ljitter_dem] + invnormal(1-(100-c(level))/200)*_se[1.polity2ljitter_dem]
scalar  lb1=low_anoc_lpolity2/anoc_se_lpolity2jitter
scalar  ub1=up_anoc_lpolity2/anoc_se_lpolity2jitter
scalar  lb2=lb1-(anoc_lpolity2jitter/anoc_se_lpolity2jitter)
scalar  ub2=ub1-(anoc_lpolity2jitter/anoc_se_lpolity2jitter)
scalar  ubd=normal(ub2)
scalar  lbd=normal(lb2)
scalar  robustness = ubd-lbd
scalar z_anoc_lpolity2 = normal(-anoc_lpolity2jitter/anoc_se_lpolity2jitter)
end

simulate rho = robustness  anoc_lpolity2jitter = anoc_lpolity2jitter low_anoc_lpolity2jitter = low_anoc_lpolity2jitter up_anoc_lpolity2jitter = up_anoc_lpolity2jitter dem_lpolity2jitter = dem_lpolity2jitter low_dem_lpolity2jitter = low_dem_lpolity2jitter up_dem_lpolity2jitter = up_dem_lpolity2jitter anoc_se_lpolity2jitter = anoc_se_lpolity2jitter dem_se_lpolity2jitter = dem_se_lpolity2jitter bic = BIC zero_anoc_lpolity2 = z_anoc_lpolity2, reps(1000) seed(1234): polity2jitter_dum
gen countzero_dem_lpolity2 = (low_dem_lpolity2jitter < 0) & (up_dem_lpolity2jitter > 0)
gen countzero_anoc_lpolity2 = (low_anoc_lpolity2jitter < 0) & (up_anoc_lpolity2jitter > 0)
describe

**without weighting
summarize
_pctile rho, p(2.5, 97.5)
matrix range_robustness = (r(r1),r(r2))
matrix rownames range_robustness = rho
matrix colnames range_robustness = lower_025 upper_975
quietly summarize rho
display "Estimated degree of effect robustness is: " r(mean)
display  "95% of unweighted simulated effect robustness ranges between: "
matrix list range_robustness








***********************************************************************************
***********************************************************************************


 

*** Measurement Error Injection Test

discard
clear
set more off

* create the test data for the monte carlo experiment



	postutil clear
	capture program drop quakesim	
	program quakesim
	version 8.0
	
	set seed 10010

		
	tempname quake5
	
	postfile `quake5' b_mag5above1015 se_mag5above1015 b_depth_min se_depth_min b_actualmaxmag10_adj se_actualmaxmag10_adj /*
			*/ b_0ncsumactualmaxmag se_0ncsumactualmaxmag /*
			*/ b_1ncsumactualmaxmag se_1ncsumactualmaxmag /*
			*/ b_dem5_polity se_dem5_polity b_corr_icrg se_corr_icrg b_lngdppc se_lngdppc b_lnpop se_lnpop using quake5, replace


			
	qui {	

	forvalues i = 1/1000 {
	
	drop _all
	
	use "Chapter_7_population_data.dta", clear
	
	gen choice=runiform()
	
	gen measurelow=0.5+runiform()

	gen measurehigh=(0.5+runiform())*2
	replace measurehigh=(0.5+runiform())/2 if runiform()<.5

	gen quakekill_hlp= killed_earthquake

	replace quakekill_hlp= killed_earthquake*measurelow 
	replace quakekill_hlp= killed_earthquake*measurehigh if (inc_developed==0 | dem5_polity==0)
	gen killed_earthquake_sim=round(quakekill_hlp) 
	drop quakekill_hlp	

	nbreg killed_earthquake_sim  mag5above1015 depth_min   actualmaxmag10_adj    notcorr_icrg_dum#c.sumactualmaxmagabove610 dem5_polity  lngdppc  corr_icrg  lnpop    if  maxmag_annual_incloffcoast>=5 &  maxmag_annual_incloffcoast<.,  cluster(country) nolrtest


	post `quake5' (_b[mag5above1015]) (_se[mag5above1015]) (_b[depth_min]) (_se[depth_min]) (_b[actualmaxmag10_adj]) (_se[actualmaxmag10_adj]) /*
			*/ (_b[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]) (_se[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]) /*
			*/ (_b[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]) (_se[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]) /*
			*/ (_b[dem5_polity]) (_se[dem5_polity]) (_b[corr_icrg]) (_se[corr_icrg]) (_b[lngdppc]) (_se[lngdppc]) (_b[lnpop]) (_se[lnpop]) 
			
		
	drop choice measure* killed_earthquake_sim

		}
	
	}


	postclose `quake5'

	end
	
di "NOTE: The simulations will take a few minutes. Please wait..."
di ""
	
	quakesim

use "quake5.dta", clear
su  b_0ncsumactualmaxmag se_0ncsumactualmaxmag b_1ncsumactualmaxmag se_1ncsumactualmaxmag	



* Estimate degrees of effect robustness
capture  program drop robdegree_sim

program robdegree_sim
version 12
drop _all

use "Chapter_7_population_data.dta", clear

set level 95
quietly xi: nbreg killed_earthquake  actualmaxmag10_adj mag5above1015 depth_min        notcorr_icrg_dum#c.sumactualmaxmagabove610   dem5_polity  lngdppc  corr_icrg  lnpop    if   maxmag_annual_incloffcoast>=5 &  maxmag_annual_incloffcoast<.,  cluster(country)
scalar  coeff_b=_b[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]
scalar  se_b=_se[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]
scalar lb_b_corr= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b_corr= coeff_b +invnormal(1-(100-c(level))/200)*se_b
scalar  coeff_b=_b[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]
scalar  se_b=_se[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]
scalar lb_b_notcorr= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b_notcorr= coeff_b +invnormal(1-(100-c(level))/200)*se_b


	*forvalues i = 1/1000 {
	
		drop _all

		use "Chapter_7_population_data.dta", clear
		
		gen choice=runiform()
	
		gen measurelow=0.5+runiform()

		gen measurehigh=(0.5+runiform())*2
		replace measurehigh=(0.5+runiform())/2 if runiform()<.5

		gen quakekill_hlp= killed_earthquake

		replace quakekill_hlp= killed_earthquake*measurelow 
		replace quakekill_hlp= killed_earthquake*measurehigh if (inc_developed==0 | dem5_polity==0)
		gen killed_earthquake_sim=round(quakekill_hlp) 
		drop quakekill_hlp	

		nbreg killed_earthquake_sim  mag5above1015 depth_min   actualmaxmag10_adj    notcorr_icrg_dum#c.sumactualmaxmagabove610 dem5_polity  lngdppc  corr_icrg  lnpop    if  maxmag_annual_incloffcoast>=5 &  maxmag_annual_incloffcoast<.,  cluster(country) nolrtest
				scalar coeff_r=_b[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]
				scalar se_r=_se[1.notcorr_icrg_dum#c.sumactualmaxmagabove610]
				scalar  lb1=lb_b_notcorr/se_r
				scalar  ub1=ub_b_notcorr/se_r
				scalar  lb2=lb1-(coeff_r/se_r)
				scalar  ub2=ub1-(coeff_r/se_r)
				scalar  ubd=normal(ub2)
				scalar  lbd=normal(lb2)
				scalar  robustness_notcorrupt=ubd-lbd
				scalar coeff_r=_b[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]
				scalar se_r=_se[0.notcorr_icrg_dum#c.sumactualmaxmagabove610]
				scalar  lb1=lb_b_corr/se_r
				scalar  ub1=ub_b_corr/se_r
				scalar  lb2=lb1-(coeff_r/se_r)
				scalar  ub2=ub1-(coeff_r/se_r)
				scalar  ubd=normal(ub2)
				scalar  lbd=normal(lb2)
				scalar  robustness_corrupt=ubd-lbd
			

	*}



end

simulate rho_notcorrupt = robustness_notcorrupt  rho_corrupt = robustness_corrupt, reps(1000) seed(1234): robdegree_sim

_pctile rho_notcorrupt, p(2.5, 97.5)
matrix range_robustness = (r(r1),r(r2))
matrix rownames range_robustness = rho_notcorrupt
matrix colnames range_robustness = lower_025 upper_975
quietly summarize rho_notcorrupt
di "Estimated degree of effect robustness in non-corrupt countries is: " r(mean)
display  "95% of unweighted simulated effect robustness ranges between: "
matrix list range_robustness

_pctile rho_corrupt, p(2.5, 97.5)
matrix range_robustness = (r(r1),r(r2))
matrix rownames range_robustness = rho_corrupt
matrix colnames range_robustness = lower_025 upper_975
quietly summarize rho_corrupt
di "Estimated degree of effect robustness in corrupt countries is: " r(mean)
display  "95% of unweighted simulated effect robustness ranges between: "
matrix list range_robustness