# Figure 6 code for
# Inconsistent definitions of GDP: Implications for estimates of decoupling
# created by Gregor Semieniuk, October 2022

# Article: https://doi.org/10.1016/j.ecolecon.2023.108000
# Zenodo repository: https://doi.org/10.5281/zenodo.8381368

library(haven)
library(tidyverse)
library(viridis)

#Implicit price deflator Table 1.1.9 BEA: 54.536 (1985), 107.747 (2017) (2012=100)) for horizontal adjustment
defl <- 107.747/54.536
defl71 <- 107.747/87.504
msiz <- 0.9

pdf("Zenodo/7.4_fig7-plot.pdf",8,4)
par(mfrow=c(1,2),mar=c(4,4.2,1,.1),lwd=1.5)
a<-read_dta("Zenodo/7.2a_fig7a_sourcedata_zgraph_smo.dta")
b<-read_dta("Zenodo/7.2b_fig7a_sourcedata_zgraph_smo_r70.dta")
d<-read_dta("Zenodo/7.2c_fig7a_sourcedata_zgraph_smo_r.dta")
e<-read_dta("Zenodo/7.2d_fig7a_sourcedata_zgraph_smo_ro.dta")
f<-read_dta("Zenodo/7.2e_fig7a_sourcedata_zgraph_smo_rna.dta")

plot(a$gdpch*defl,a$ypred,type="l",xlim=c(0,35),ylim=c(5,140),col=viridis(6)[1],xlab="2017 int. kUSD/person",ylab = "",lwd=2)
lines(b$gdpch*defl71,b$ypred,type="l",col=viridis(6)[2],lwd=2)
lines(d$gdpch,d$ypred,type="l",col=viridis(6)[3],lwd=2)
lines(e$gdpch,e$ypred,type="l",col=viridis(6)[4],lwd=2)
lines(f$gdpch,f$ypred,type="l",col=viridis(6)[5],lwd=2)
text(9.6,95,"mark 5",pos = 2, col = viridis(5)[1], cex = msiz)
text(18,97,"7.0",pos = 4, col = viridis(6)[2], cex = msiz)
text(29,56.3,"10.0\nrgdpe", col = viridis(6)[3], pos = 1, cex = msiz)
text(26.5,74.5,"10.0\nrgdpo", col = viridis(6)[4], cex = msiz)
text(36,133,"10.0 rgdpna", pos = 2, col = viridis(6)[5], cex = msiz)
mtext(side = 3, line=0, adj=0,"Smoke concentration in cities")
#legend("bottomright",legend= c("mark 5, rgdpch", "7.0, rgdpch", "10.0 rdpe", "10.0 rgdpo", "10.0 rdgpna"), lty=1,lwd=3,col=viridis(6)[1:5], title = "PWT version & series", ncol=2, bty="n", cex=0.8)
text(37,4,"PWT mark 5 & 7.0 use rgdpch", pos = 2, cex = msiz)
mtext(side=2,line =2.7, text = expression(paste(mu,"g/",m^3)))
mtext(side = 3, line =-0.3, adj=-0.23, "a)" ,font =2)

a<-read_dta("Zenodo/7.3a_fig7b_sourcedata_zgraph_mer.dta")
b<-read_dta("Zenodo/7.3b_fig7b_sourcedata_zgraph_mer_r70.dta")
d<-read_dta("Zenodo/7.3c_fig7b_sourcedata_zgraph_mer_r.dta")
e<-read_dta("Zenodo/7.3d_fig7b_sourcedata_zgraph_mer_ro.dta")
f<-read_dta("Zenodo/7.3e_fig7b_sourcedata_zgraph_mer_rna.dta")

plot(a$gdpch*defl,a$ypred,type="l",xlim=c(0,36),col=viridis(6)[1],ylim=c(-0.3,.85),xlab="2017 int. kUSD/person",ylab = "",lwd=2)
lines(b$gdpch*defl71,b$ypred,type="l",col=viridis(6)[2],lwd=2)
lines(d$gdpch,d$ypred,type="l",col=viridis(6)[3],lwd=2)
lines(e$gdpch,e$ypred,type="l",col=viridis(6)[4],lwd=2)
lines(f$gdpch,f$ypred,type="l",col=viridis(6)[5],lwd=2)
text(33.4,0.79,"mark 5",pos = 2, col = viridis(5)[1], cex = msiz)
text(34.8,0.6,"7.0", pos = 3, col = viridis(6)[2], cex = msiz)
text(4,-.14,"10.0\nrgdpe", col = viridis(6)[3], cex = msiz)
text(7,0.075,"10.0 rgdpo", pos = 4, col = viridis(6)[4], cex = msiz)
text(3.3,.3,"10.0 rgdpna", pos = 4, col = viridis(6)[5], cex = msiz)
mtext(side = 3, line=0, adj=0,"Mercury concentration in rivers")
mtext(side=2,line =2.7, text = expression(paste(mu,"g/l")))
mtext(side = 3, line =-0.3, adj=-0.23, "b)" ,font =2)

dev.off()

# datapoints file
a<-read_dta("Zenodo/7.2a_fig7a_sourcedata_zgraph_smo.dta") %>% 
  mutate(gdpch = gdpch*defl, version = "mark5") %>% 
  select(gdpch,ypred, version)
b<-read_dta("Zenodo/7.2b_fig7a_sourcedata_zgraph_smo_r70.dta") %>% 
  mutate(gdpch = gdpch*defl71, version = "7.0" ) %>% 
  select(gdpch,ypred, version) 
d<-read_dta("Zenodo/7.2c_fig7a_sourcedata_zgraph_smo_r.dta")  %>% 
  mutate(version = "10.0 rgdpe") %>% 
  select(gdpch,ypred, version)
e<-read_dta("Zenodo/7.2d_fig7a_sourcedata_zgraph_smo_ro.dta")  %>% 
  mutate(version = "10.0 rgdpo") %>% 
  select(gdpch,ypred, version)
f<-read_dta("Zenodo/7.2e_fig7a_sourcedata_zgraph_smo_rna.dta")  %>% 
  mutate(version = "10.0 rgdpna") %>% 
  select(gdpch,ypred, version)
smoke <-bind_rows(a,b,d,e,f)

a<-read_dta("Zenodo/7.3a_fig7b_sourcedata_zgraph_mer.dta") %>% 
  mutate(gdpch = gdpch*defl, version = "mark5") %>% 
  select(gdpch,ypred, version)
b<-read_dta("Zenodo/7.3b_fig7b_sourcedata_zgraph_mer_r70.dta") %>% 
  mutate(gdpch = gdpch*defl71, version = "7.0" ) %>% 
  select(gdpch,ypred, version) 
d<-read_dta("Zenodo/7.3c_fig7b_sourcedata_zgraph_mer_r.dta")  %>% 
  mutate(version = "10.0 rgdpe") %>% 
  select(gdpch,ypred, version)
e<-read_dta("Zenodo/7.3d_fig7b_sourcedata_zgraph_mer_ro.dta")  %>% 
  mutate(version = "10.0 rgdpo") %>% 
  select(gdpch,ypred, version)
f<-read_dta("Zenodo/7.3e_fig7b_sourcedata_zgraph_mer_rna.dta")  %>% 
  mutate(version = "10.0 rgdpna") %>% 
  select(gdpch,ypred, version)
mercury <-bind_rows(a,b,d,e,f)

write_csv(smoke,"Zenodo/7.5_fig7a_datapoints.csv")
write_csv(mercury,"Zenodo/7.6_fig7b_datapoints.csv")