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

Replication of empirical results in Chapter 13: Effect Dynamics

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)

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



*** Effect Onset test
* Baseline original model
use "Chapter_13_dynamics.dta", clear
xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
nlcom((_b[left])/(1-_b[lspend]))
scalar  coeff_b=_b[left]
scalar  se_b=_se[left]
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 table1, replace excel 

* Robustness test model
*optimize left and cdem lags per country using r2 
encode country, gen(countryfac) label(country)
drop country
gen country = countryfac
drop countryfac

gen leftlag = left
forval Country = 1/18 {
scalar r2best = 0
scalar leftlag`Country' = 0
forval lag = 0/3 {
drop leftlag
gen leftlag = left
replace leftlag = L`lag'.left if country == `Country'
replace leftlag = left if missing(leftlag)
quietly xtpcse spend unem growthpc depratio leftlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar leftlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
gen cdemlag = cdem
forval Country = 1/18 {
scalar r2best = 0
scalar cdemlag`Country' = 0
forval lag = 0/3 {
drop cdemlag
gen cdemlag = cdem
replace cdemlag = L`lag'.cdem if country == `Country'
replace cdemlag = cdem if missing(cdemlag)
quietly xtpcse spend unem growthpc depratio cdemlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar cdemlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
forval Country = 1/18 {
*display `Country'
scalar bestlag = leftlagoptim`Country'
display bestlag
replace leftlag = L`=bestlag'.left if country==`Country'
replace leftlag = left if missing(leftlag)
}

drop cdemlag
gen cdemlag = cdem
forval Country = 1/18 {
*display `Country'
scalar bestlag = cdemlagoptim`Country'
display bestlag
replace cdemlag = L`=bestlag'.cdem if country==`Country'
replace cdemlag = cdem if missing(cdemlag)
}




xi: xtpcse spend lspend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, 
outreg2 using table1, append excel 

* Estimating degree of effect robustness
scalar coeff_r=_b[leftlag]
scalar se_r=_se[leftlag]
nlcom((_b[leftlag])/(1-_b[lspend]))
nlcom((_b[cdemlag])/(1-_b[lspend]))

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 LEFT is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[cdemlag]
scalar se_r=_se[cdemlag]
quietly xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
scalar  coeff_b=_b[cdem]
scalar  se_b=_se[cdem]
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

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 RIGHT is: " robustness


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


*** Effect Onset Test with Alternative Dynamic Specification Test 
** Baseline original model

* Baseline original model
use "Chapter_13_dynamics.dta", clear
xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
nlcom((_b[left])/(1-_b[lspend]))
scalar  coeff_b=_b[left]
scalar  se_b=_se[left]
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 table1, replace excel 

* Robustness test model
*optimize left and cdem lags per country using r2 
encode country, gen(countryfac) label(country)
drop country
gen country = countryfac
drop countryfac

gen leftlag = left
forval Country = 1/18 {
scalar r2best = 0
scalar leftlag`Country' = 0
forval lag = 0/3 {
drop leftlag
gen leftlag = left
replace leftlag = L`lag'.left if country == `Country'
replace leftlag = left if missing(leftlag)
quietly xtpcse spend unem growthpc depratio leftlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar leftlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
gen cdemlag = cdem
forval Country = 1/18 {
scalar r2best = 0
scalar cdemlag`Country' = 0
forval lag = 0/3 {
drop cdemlag
gen cdemlag = cdem
replace cdemlag = L`lag'.cdem if country == `Country'
replace cdemlag = cdem if missing(cdemlag)
quietly xtpcse spend unem growthpc depratio cdemlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar cdemlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
forval Country = 1/18 {
*display `Country'
scalar bestlag = leftlagoptim`Country'
display bestlag
replace leftlag = L`=bestlag'.left if country==`Country'
replace leftlag = left if missing(leftlag)
}

drop cdemlag
gen cdemlag = cdem
forval Country = 1/18 {
*display `Country'
scalar bestlag = cdemlagoptim`Country'
display bestlag
replace cdemlag = L`=bestlag'.cdem if country==`Country'
replace cdemlag = cdem if missing(cdemlag)
}



xtpcse spend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, corr(ar1)
outreg2 using table1, append excel 
scalar coeff_r=_b[leftlag]
scalar se_r=_se[leftlag]

* 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 LEFT is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[cdemlag]
scalar se_r=_se[cdemlag]
quietly xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
scalar  coeff_b=_b[cdem]
scalar  se_b=_se[cdem]
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

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 RIGHT is: " robustness



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



*** Effect Onset Test with Alternative Dynamic Specification Test 
** Baseline LDV and optimized onset model

* Baseline LDV and optimized onset model
use "Chapter_13_dynamics.dta", clear

*optimize left and cdem lags per country using r2 
encode country, gen(countryfac) label(country)
drop country
gen country = countryfac
drop countryfac

gen leftlag = left
forval Country = 1/18 {
scalar r2best = 0
scalar leftlag`Country' = 0
forval lag = 0/3 {
drop leftlag
gen leftlag = left
replace leftlag = L`lag'.left if country == `Country'
replace leftlag = left if missing(leftlag)
quietly xtpcse spend unem growthpc depratio leftlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar leftlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
gen cdemlag = cdem
forval Country = 1/18 {
scalar r2best = 0
scalar cdemlag`Country' = 0
forval lag = 0/3 {
drop cdemlag
gen cdemlag = cdem
replace cdemlag = L`lag'.cdem if country == `Country'
replace cdemlag = cdem if missing(cdemlag)
quietly xtpcse spend unem growthpc depratio cdemlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar cdemlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
forval Country = 1/18 {
*display `Country'
scalar bestlag = leftlagoptim`Country'
display bestlag
replace leftlag = L`=bestlag'.left if country==`Country'
replace leftlag = left if missing(leftlag)
}

drop cdemlag
gen cdemlag = cdem
forval Country = 1/18 {
*display `Country'
scalar bestlag = cdemlagoptim`Country'
display bestlag
replace cdemlag = L`=bestlag'.cdem if country==`Country'
replace cdemlag = cdem if missing(cdemlag)
}



quietly xtpcse spend lspend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, 
scalar  coeff_b=_b[leftlag]
scalar  se_b=_se[leftlag]


* Robustness test model
xi: xtpcse spend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, corr(ar1)
scalar coeff_r=_b[leftlag]
scalar se_r=_se[leftlag]

* Estimating degree of effect robustness
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
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 LEFT is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[cdemlag]
scalar se_r=_se[cdemlag]
quietly xtpcse spend lspend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, 
scalar  coeff_b=_b[cdemlag]
scalar  se_b=_se[cdemlag]
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

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 RIGHT is: " robustness



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


*** Effect Onset and Random Coefficients Dynamics Tests
** Baseline original model

* Baseline original model
use "Chapter_13_dynamics.dta", clear
xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
nlcom((_b[left])/(1-_b[lspend]))
scalar  coeff_b=_b[left]
scalar  se_b=_se[left]
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 table1, replace excel 

* Robustness test model
*optimize left and cdem lags per country using r2 
encode country, gen(countryfac) label(country)
drop country
gen country = countryfac
drop countryfac

gen leftlag = left
forval Country = 1/18 {
scalar r2best = 0
scalar leftlag`Country' = 0
forval lag = 0/3 {
drop leftlag
gen leftlag = left
replace leftlag = L`lag'.left if country == `Country'
replace leftlag = left if missing(leftlag)
quietly xtpcse spend unem growthpc depratio leftlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar leftlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
gen cdemlag = cdem
forval Country = 1/18 {
scalar r2best = 0
scalar cdemlag`Country' = 0
forval lag = 0/3 {
drop cdemlag
gen cdemlag = cdem
replace cdemlag = L`lag'.cdem if country == `Country'
replace cdemlag = cdem if missing(cdemlag)
quietly xtpcse spend unem growthpc depratio cdemlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar cdemlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
forval Country = 1/18 {
*display `Country'
scalar bestlag = leftlagoptim`Country'
display bestlag
replace leftlag = L`=bestlag'.left if country==`Country'
replace leftlag = left if missing(leftlag)
}

drop cdemlag
gen cdemlag = cdem
forval Country = 1/18 {
*display `Country'
scalar bestlag = cdemlagoptim`Country'
display bestlag
replace cdemlag = L`=bestlag'.cdem if country==`Country'
replace cdemlag = cdem if missing(cdemlag)
}


xtmixed spend lspend unem growthpc depratio trade lowwage fdi leftlag cdemlag i.country i.year country#c.year, cluster(cc) || country: lspend unem growthpc depratio trade lowwage fdi leftlag cdemlag
scalar coeff_r=_b[leftlag]
scalar se_r=_se[leftlag]

* 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 LEFT is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[cdemlag]
scalar se_r=_se[cdemlag]
quietly xi: xtpcse spend lspend  unem growthpc depratio left cdem trade lowwage fdi i.country i.year
scalar  coeff_b=_b[cdem]
scalar  se_b=_se[cdem]
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

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 RIGHT is: " robustness




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


*** Effect Onset and Random Coefficients Dynamics Tests
** Baseline LDV and optimized onset model

use "Chapter_13_dynamics.dta", clear

*optimize left and cdem lags per country using r2 
encode country, gen(countryfac) label(country)
drop country
gen country = countryfac
drop countryfac

gen leftlag = left
forval Country = 1/18 {
scalar r2best = 0
scalar leftlag`Country' = 0
forval lag = 0/3 {
drop leftlag
gen leftlag = left
replace leftlag = L`lag'.left if country == `Country'
replace leftlag = left if missing(leftlag)
quietly xtpcse spend unem growthpc depratio leftlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar leftlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
gen cdemlag = cdem
forval Country = 1/18 {
scalar r2best = 0
scalar cdemlag`Country' = 0
forval lag = 0/3 {
drop cdemlag
gen cdemlag = cdem
replace cdemlag = L`lag'.cdem if country == `Country'
replace cdemlag = cdem if missing(cdemlag)
quietly xtpcse spend unem growthpc depratio cdemlag cdem trade lowwage fdi i.country, correlation(psar1)
scalar r2 = e(r2)
if r2 > r2best scalar cdemlagoptim`Country' = `lag'
if r2 > r2best scalar r2best = r2
}
}

drop leftlag
gen leftlag = left
forval Country = 1/18 {
*display `Country'
scalar bestlag = leftlagoptim`Country'
display bestlag
replace leftlag = L`=bestlag'.left if country==`Country'
replace leftlag = left if missing(leftlag)
}

drop cdemlag
gen cdemlag = cdem
forval Country = 1/18 {
*display `Country'
scalar bestlag = cdemlagoptim`Country'
display bestlag
replace cdemlag = L`=bestlag'.cdem if country==`Country'
replace cdemlag = cdem if missing(cdemlag)
}



quietly xtpcse spend lspend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, 
scalar  coeff_b=_b[leftlag]
scalar  se_b=_se[leftlag]


* Robustness test model
xtmixed spend lspend unem growthpc depratio trade lowwage fdi leftlag cdemlag i.country i.year country#c.year, cluster(cc) || country: lspend unem growthpc depratio trade lowwage fdi leftlag cdemlag
outreg2 using table1, append excel 
nlcom((_b[leftlag])/(1-_b[lspend]))
nlcom((_b[cdemlag])/(1-_b[lspend]))
scalar coeff_r=_b[leftlag]
scalar se_r=_se[leftlag]

* Estimating degree of effect robustness (short-run)
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 LEFT is: " robustness

* Estimating degree of effect robustness
scalar coeff_r=_b[cdemlag]
scalar se_r=_se[cdemlag]
quietly xtpcse spend lspend  unem growthpc depratio trade lowwage fdi leftlag cdemlag  i.country i.year country#c.year, 
scalar  coeff_b=_b[cdemlag]
scalar  se_b=_se[cdemlag]
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

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 RIGHT is: " robustness