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

Replication of empirical results in Chapter 10: Functional forms beyond default

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)

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


*** Higher-degree polynomials test

** 2nd degree polynomial
* Baseline
use "Chapter_10_functional.dta", replace
reg  emissions_pc gdppc i.year, cluster(country)
outreg2 using table1, replace excel level(95)  e(rmse)
margins, exp((_b[gdppc]*gdppc)) at(gdppc=(1(1)88)) saving(margins_base, replace)
capture drop coeff_b
capture drop se_b
predictnl coeff_b=((_b[gdppc]*gdppc)), se(se_b)  

* Robustness test model
reg emissions_pc c.gdppc##c.gdppc i.year , cluster(country)
capture drop coeff_r
capture drop se_r
predictnl coeff_r=((_b[gdppc]*gdppc+_b[c.gdppc*c.gdppc]*gdppc*gdppc)), se(se_r)   
outreg2 using table1, append excel level(95)  e(rmse)
margins, exp(_b[gdppc]*gdppc+_b[c.gdppc*c.gdppc]*gdppc*gdppc) at(gdppc=(1(1)88)) saving(margins_rob, replace)
combomarginsplot margins_base margins_rob, labels("Baseline" "Robustness test") title(Predicted CO2 per capita emissions) ytitle("") xtitle("GDP per capita in thousand USD") graphregion(color(white)) bgcolor(white)
graph save "Higher-degree polynomial test 1", replace
graph export "Higher-degree polynomial test 1.eps", as(eps) replace



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


*** Higher-order polynomials test

** 3rd degree polynomial
* Baseline
use "Chapter_10_functional.dta", replace
reg  emissions_pc gdppc  i.year, cluster(country)
margins, exp(_b[gdppc]*gdppc) at(gdppc=(1(1)88)) saving(margins_base, replace)
capture drop coeff_b
capture drop se_b
predictnl coeff_b=((_b[gdppc]*gdppc)), se(se_b)  

* Robustness test model
reg emissions_pc c.gdppc##c.gdppc##c.gdppc  i.year, cluster(country)
outreg2 using table1, append excel level(95)  e(rmse)
margins, exp(_b[gdppc]*gdppc+_b[c.gdppc*c.gdppc]*gdppc*gdppc+_b[c.gdppc*c.gdppc*c.gdppc]*gdppc*gdppc*gdppc) at(gdppc=(1(1)88)) saving(margins_rob, replace)
combomarginsplot margins_base margins_rob, labels("Baseline" "Robustness test") title(Predicted CO2 per capita emissions) ytitle("") xtitle("GDP per capita in thousand USD") graphregion(color(white)) bgcolor(white)
graph save "Higher-degree polynomial test 2", replace
graph export "Higher-degree polynomial test 2.eps", as(eps) replace

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

** 4th degree polynomial
* Baseline
use "Chapter_10_functional.dta", replace
reg  emissions_pc gdppc  i.year, cluster(country) 
margins, exp(_b[gdppc]*gdppc) at(gdppc=(1(1)88)) saving(margins_base, replace)
capture drop coeff_b
capture drop se_b
predictnl coeff_b=((_b[gdppc]*gdppc)), se(se_b)  

* Robustness test model
reg emissions_pc c.gdppc##c.gdppc##c.gdppc##c.gdppc  i.year, cluster(country)
outreg2 using table1, append excel level(95)  e(rmse)
margins, exp(_b[gdppc]*gdppc+_b[c.gdppc*c.gdppc]*gdppc*gdppc+_b[c.gdppc*c.gdppc*c.gdppc]*gdppc*gdppc*gdppc+_b[c.gdppc*c.gdppc*c.gdppc*c.gdppc]*gdppc*gdppc*gdppc*gdppc) at(gdppc=(1(1)88)) saving(margins_rob, replace)
combomarginsplot margins_base margins_rob, labels("Baseline" "Robustness test") title(Predicted CO2 per capita emissions) ytitle("") xtitle("GDP per capita in thousand USD") graphregion(color(white)) bgcolor(white)



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


*** Semi-parametric test
** Percentiles

use "Chapter_10_functional.dta", replace

* Baseline
reg co2 gdppc i.year, cluster(country)
outreg2 using table2, replace excel level(95)  e(rmse)
predictnl coeff_b=((_b[gdppc]*gdppc)), se(se_b)  

* Robustness test model
capture dropvars gdppc_percentiles
quietly reg co2 gdppc
xtile gdppc_percentiles=gdppc, nq(10)
tab gdppc_percentiles
capture dropvars percentile*
gen percentile1=gdppc_percentiles==1
gen percentile2=gdppc_percentiles==2
gen percentile3=gdppc_percentiles==3
gen percentile4=gdppc_percentiles==4
gen percentile5=gdppc_percentiles==5
gen percentile6=gdppc_percentiles==6
gen percentile7=gdppc_percentiles==7
gen percentile8=gdppc_percentiles==8
gen percentile9=gdppc_percentiles==9
gen percentile10=gdppc_percentiles==10
reg co2 percentile* i.year, nocons cluster(countryid)
outreg2 using table2, append excel level(95)   e(rmse)
predictnl coeff_r=(xb() ), se(se_r)   

* Estimating degree of effect robustness
gen lb_b= coeff_b -invnormal(0.975)*se_b
gen ub_b= coeff_b +invnormal(0.975)*se_b
gen lb1=lb_b/se_r
gen  ub1=ub_b/se_r
gen   lb2=lb1-(coeff_r/se_r)
gen   ub2=ub1-(coeff_r/se_r)
gen    ubd=normal(ub2)
gen   lbd=normal(lb2)
gen   robustness=ubd-lbd
collapse robustness gdppc, by(percentile*)
twoway scatter robustness gdppc, title("") ytitle("") xtitle("GDP per capita in thousand USD") graphregion(color(white)) bgcolor(white)




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



*** Semi-parametric test
** Equal bin width 

use "Chapter_10_functional.dta", replace

* Baseline
reg co2 gdppc i.year, cluster(country)
predictnl coeff_b=((_b[gdppc]*gdppc)), se(se_b)  

* Robustness test model
capture dropvars below*
gen below10=0
replace below10=1 if gdppc<=10
gen below20=0
replace below20=1 if gdppc<=20 & gdppc>10
gen below30=0
replace below30=1 if gdppc<=30 & gdppc>20
gen below40=0
replace below40=1 if gdppc<=40 & gdppc>30
gen below50=0
replace below50=1 if gdppc<=50 & gdppc>40
gen below60=0
replace below60=1 if gdppc<=60 & gdppc>50
gen below70=0
replace below70=1 if gdppc<=70 & gdppc>60
gen below80=0
replace below80=1 if gdppc<=80 & gdppc>70
gen below90=0
replace below90=1 if gdppc<=90 & gdppc>80
reg co2 below* i.year, noconst cluster(countryid)
outreg2 using table2, replace excel level(95)  e(rmse)
predictnl coeff_r=(xb() ), se(se_r)   

* Estimating degree of effect robustness
gen lb_b= coeff_b -invnormal(0.975)*se_b
gen ub_b= coeff_b +invnormal(0.975)*se_b
gen lb1=lb_b/se_r
gen  ub1=ub_b/se_r
gen   lb2=lb1-(coeff_r/se_r)
gen   ub2=ub1-(coeff_r/se_r)
gen    ubd=normal(ub2)
gen   lbd=normal(lb2)
gen   robustness=ubd-lbd
collapse robustness gdppc, by(below*)
twoway scatter robustness gdppc, title("") ytitle("") xtitle("GDP per capita in thousand USD") graphregion(color(white)) bgcolor(white)