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

Replication of empirical results in Chapter 14: Spatial Correlation and Dependence

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)

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



*** Unit fixed effects and lagged dependent variable test

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
tsset
capture drop lsumtot
gen lsumtot=l.sumtot
xi: xtpoisson sumtot lsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, fe
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel



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



*** Year fixed effects test

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots i.year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel



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



*** Unit and year fixed effects and lagged dependent variable test

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
tsset
capture drop lsumtot
gen lsumtot=l.sumtot
xi: xtpoisson sumtot lsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots i.year if year>2000, fe
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel



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


*** Spatial error tests
** connectivity variable: contiguity

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots  serr_lnsumtot_contig_tots year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, replace excel


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



*** Spatial error tests
** connectivity variable: inverse distance

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots  serr_lnsumtot_inv_dist_tots year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel


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



*** Spatial error tests
** connectivity variable: common region

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots  serr_lnsumtot_comreg_tots year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel


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



*** Spatial error tests
** connectivity variable: common language

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots  serr_lnsumtot_comlang_tots year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel



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




*** Alternative connectivity test 

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots sl_lnsumtot_w_cult_grp2_peers sl_lnsumtot_cinc_tots year if year>2000, cluster(country)
scalar coeff_r=_b[sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[sl_pi_lnsumtot_wpeers_tots]

* Estimating degree of effect robustness
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness is: " robustness
outreg2 using table, append excel


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



*** Heterogeneous exposure test

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
replace sl_pi_lnsumtot_wi_k_low_tots=0 if sl_pi_lnsumtot_wi_k_low_tots==.
replace sl_pi_lnsumtot_wi_k_mid_tots=0 if sl_pi_lnsumtot_wi_k_mid_tots==.
replace sl_pi_lnsumtot_wi_k_high_tots=0 if sl_pi_lnsumtot_wi_k_high_tots==.
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wi_k_low_tots sl_pi_lnsumtot_wi_k_mid_tots sl_pi_lnsumtot_wi_k_high_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
outreg2 using table, replace excel
test sl_pi_lnsumtot_wi_k_low_tots= sl_pi_lnsumtot_wi_k_mid_tots= sl_pi_lnsumtot_wi_k_high_tots , mtest

* Estimating degree of effect robustness
scalar coeff_r=_b[sl_pi_lnsumtot_wi_k_high_tots]
scalar se_r=_se[sl_pi_lnsumtot_wi_k_high_tots]
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness in peers of high exposure is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[sl_pi_lnsumtot_wi_k_mid_tots]
scalar se_r=_se[sl_pi_lnsumtot_wi_k_mid_tots]
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness in peers of medium exposure is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[sl_pi_lnsumtot_wi_k_low_tots]
scalar se_r=_se[sl_pi_lnsumtot_wi_k_low_tots]
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness in peers of low exposure is: " robustness


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



*** Heterogeneous responsiveness test

use "Chapter_14_spatial.dta", clear
* Baseline
set level 95
xi: poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
margins, at((mean) _all (p25) sl_pi_lnsumtot_wpeers_tots )
margins, at((mean) _all (p75) sl_pi_lnsumtot_wpeers_tots )
scalar  coeff_b=_b[sl_pi_lnsumtot_wpeers_tots]
scalar  se_b=_se[sl_pi_lnsumtot_wpeers_tots]
scalar lb_b= coeff_b -invnormal(1-(100-c(level))/200)*se_b
scalar ub_b= coeff_b +invnormal(1-(100-c(level))/200)*se_b
outreg2 using table, replace excel 

* Robustness test model
capture drop sumtotasobservedann_3grps
egen sumtotasobservedann_3grps=mean(sumtot), by(i_pi_totsum_rank_3grps year)
capture drop belowasobservedann_3grps
gen belowasobservedann_3grps=0
replace belowasobservedann_3grps=. if sumtot==.
replace belowasobservedann_3grps=1 if sumtot<sumtotasobservedann_3grps
capture drop sl_pi_wpeers_above sl_pi_wpeers_below
gen sl_pi_wpeers_above=sl_pi_lnsumtot_wpeers_tots*(1-belowasobservedann_3grps)
gen sl_pi_wpeers_below=sl_pi_lnsumtot_wpeers_tots*(belowasobservedann_3grps)
poisson sumtot icsumtot gdppcconst gov_left1 gov_right1 predinc_total belowasobservedann_3grps belowasobservedann_3grps#c.sl_pi_lnsumtot_wpeers_tots sl_pi_lnsumtot_wnonpeers_tots year if year>2000, cluster(country)
outreg2 using table, append excel

* Estimating degree of effect robustness
scalar coeff_r=_b[0.belowasobservedann_3grps#c.sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[0.belowasobservedann_3grps#c.sl_pi_lnsumtot_wpeers_tots]
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness in countries above peer group mean is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[1.belowasobservedann_3grps#c.sl_pi_lnsumtot_wpeers_tots]
scalar se_r=_se[1.belowasobservedann_3grps#c.sl_pi_lnsumtot_wpeers_tots]
scalar lb1=lb_b/se_r
scalar ub1=ub_b/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=ubd-lbd
di "Estimated degree of effect robustness in countries below peer group mean is: " robustness