# Figure 5 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(tidyverse)
library(readxl)
library(lazyeval)
library(countrycode)
library(viridis)

load(file = "Zenodo/5.2_sourcedata_energy.RData")
load(file = "Zenodo/5.3_sourcedata_materials.RData")

master <- masterdata_materials 
master_en <- masterdata_energy %>% 
  arrange(isocode, year)

rm(masterdata_materials)
rm(masterdata_energy)

# functions that selected two GDP vintages for materials and energy (suffix _en)
gdp_select <- function(a,b) {
  twoyears <- master %>% select(uq(a),uq(b),year)
  twoyears <-twoyears[complete.cases(twoyears),]
  names(twoyears)<-c("x","y","year")#,"iso")
  twoyears <- twoyears %>% mutate(flip = case_when(
    x<=0 & y>0 ~ 1,  # retroactive change to rel decoupling
    x>0 & y<=0 ~ -1,  # retroactive change to recoupling
    TRUE ~ 0
  ))
  return(twoyears)
}

gdp_select_en <- function(a,b) {
  twoyears <- master_en %>% select(uq(a),uq(b),year)#,country)
  twoyears <-twoyears[complete.cases(twoyears),]
  #twoyears <- twoyears
  names(twoyears)<-c("x","y","year")#,"iso")
  twoyears <- twoyears %>% mutate(flip = case_when(
    x<=0 & y>0 ~ 1,  # retroactive change to rel decoupling
    x>0 & y<=0 ~ -1,  # retroactive change to recoupling
    TRUE ~ 0
  ))
  return(twoyears)
}

# parameters for plot
labelx = 0.705
plcex <- 0.75

pdf("Zenodo/5.4_fig5_plot.pdf",10,13)
par(mar=c(.2,.2,0.2,0.2),oma=c(3.9,3.9,1.3,0),mfcol=c(4,2))

dd<-gdp_select_en(~ey_gr56, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 5.6 & 10.0", adj=0.06,cex = plcex)
mtext(side =3, line = -1.3,"a)", adj=0.01, font = 2)

dd<-gdp_select_en(~ey_gr61, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 6.1 & 10.0", adj=0.06,cex = plcex)
mtext(side =3, line = -1.3,"b)", adj=0.01, font = 2)

dd<-gdp_select_en(~ey_gr70, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 7.0 & 10.0", adj=0.06,cex=plcex)
mtext(side =3, line = -1.3,"c)", adj=0.01, font = 2)

dd<-gdp_select_en(~ey_gr90, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 9.0 & 10.0", adj=0.06,cex=plcex)
mtext(side =3, line = -1.3,"d)", adj=0.01, font = 2)

dd<-gdp_select(~ey_gr56, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n",yaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 5.6 & 10.0", adj=0.06,cex=plcex)
mtext(side =3, line = -1.3,"e)", adj=0.01, font = 2)

dd<-gdp_select(~ey_gr61, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n",yaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 6.1 & 10.0", adj=0.06,cex= plcex)
mtext(side =3, line = -1.3,"f)", adj=0.01, font = 2)

dd<-gdp_select(~ey_gr70, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",xaxt="n",yaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 7.0 & 10.0", adj=0.06,cex=plcex)
mtext(side =3, line = -1.3,"g)", adj=0.01, font = 2)

dd<-gdp_select(~ey_gr90, ~ey_gr100)
plot(0,0, xlim = c(-0.08,0.1), ylim = c(-0.08,0.11), ann=F,pch=20,col="white",yaxt="n")
abline(v=0,lwd=0.5);abline(h=0,lwd=0.5)
abline(-.02,1,lwd=0.5);abline(.02,1,lwd=0.5)
points(filter(dd,flip==0)[1:2],col=alpha(viridis(4)[2],0.3),pch=20)
points(filter(dd,flip==-1)[1:2],col=alpha(viridis(4)[3],0.4),pch=20)
points(filter(dd,flip==1)[1:2],col=alpha(viridis(4)[1],0.4),pch=20)
text(-0.045,0.065,"Retroactive\nrecoupling",col=viridis(4)[1])
text(0.075,-0.07,"Retroactive\nrelative\ndecoupling",col=viridis(4)[3])
text(0,0.11,paste("Corr.:", round(cor(dd[1:2])[2],3)),pos=4)
mtext(side =3, line = -1.3,"PWT 9.0 & 10.0", adj=0.06,cex=plcex)
mtext(side =3, line = -1.3,"h)", adj=0.01, font = 2)

mtext("10 year average energy intensity rate of change, older vintage", side = 1, line =2.6,
      outer = TRUE)

mtext("10 year average energy intensity rate of change, newer vintage", side = 2.5, line =2, 
      outer = TRUE)

mtext("Energy, 10 year rates of change", side = 3, line = 0, outer = T,adj = 0.16)
mtext("Materials, 10 year rates of change", side = 3, line = 0, outer = T,adj = 0.85)
dev.off()

# source data by panel
dd<-gdp_select_en(~ey_gr56, ~ey_gr100)
write_csv(dd,"Zenodo/5.5a_fig5a_datapoints.csv")
dd<-gdp_select_en(~ey_gr61, ~ey_gr100)
write_csv(dd,"Zenodo/5.5b_fig5b_datapoints.csv")
dd<-gdp_select_en(~ey_gr70, ~ey_gr100)
write_csv(dd,"Zenodo/5.5c_fig5c_datapoints.csv")
dd<-gdp_select_en(~ey_gr90, ~ey_gr100)
write_csv(dd,"Zenodo/5.5d_fig5d_datapoints.csv")
dd<-gdp_select(~ey_gr56, ~ey_gr100)
write_csv(dd,"Zenodo/5.5e_fig5e_datapoints.csv")
dd<-gdp_select(~ey_gr61, ~ey_gr100)
write_csv(dd,"Zenodo/5.5f_fig5f_datapoints.csv")
dd<-gdp_select(~ey_gr70, ~ey_gr100)
write_csv(dd,"Zenodo/5.5g_fig5g_datapoints.csv")
dd<-gdp_select(~ey_gr90, ~ey_gr100)
write_csv(dd,"Zenodo/5.5h_fig5h_datapoints.csv")