# 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(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 to count absolute decoupling flips
gdp_abs <- function(a,b,d) {
  twoyears <- master %>% select(uq(a),uq(b),uq(d),year,isocode)
  twoyears <-twoyears[complete.cases(twoyears),]
  names(twoyears)<-c("x","y","z","year","isocode")
  twoyears <- twoyears %>% 
    mutate(abs = ifelse(x>0 & z<=0,1,0),
           abs_on = ifelse(x<=0 & y>0 & z <=0,1,0),
           abs_off = ifelse(x>0 & y<=0 & z <=0,1,0)) %>% 
  group_by(year) %>% 
    summarize(total = sum(abs),
        count = n(),
        total_on = sum(abs_on),
        total_off = sum(abs_off),
        total = sum(abs)+total_off,
        share_on = total_on/total,
        share_off = total_off/total,
        share_abs = total/count
   ) %>% 
  return(twoyears)
}

#~GDP_gr56,~GDP_gr56, ~ener_gr_10
twoyears <- master %>% select(GDP_gr56,GDP_gr100,DMC_gr_10,year,isocode) %>% 
  filter(year == 1991) 
twoyears <-twoyears[complete.cases(twoyears),]
names(twoyears)<-c("x","y","z","year","isocode")
twoyears <- twoyears %>% 
  mutate(abs = ifelse(x>0 & z<=0,1,0),
         abs_on = ifelse(x<=0 & y>0 & z <=0,1,0),
         abs_off = ifelse(x>0 & y<=0 & z <=0,1,0),
         sign_up = ifelse(x<=0 & y>0,1,0),
         sign_down = ifelse(x>0 & y<=0,1,0))

gdp_abs_en <- function(a,b,d) {
  twoyears <- master_en %>% select(uq(a),uq(b),uq(d),year,isocode)
  twoyears <-twoyears[complete.cases(twoyears),]
  names(twoyears)<-c("x","y","z","year","isocode")
  twoyears <- twoyears %>% 
    mutate(abs = ifelse(x>0 & z<=0,1,0),
           abs_on = ifelse(x<=0 & y>0 & z <=0,1,0),
           abs_off = ifelse(x>0 & y<=0 & z <=0,1,0)) %>% 
    group_by(year) %>% 
    summarize(total = sum(abs),
              count = n(),
              total_on = sum(abs_on),
              total_off = sum(abs_off),
              total = sum(abs)+total_off,
              share_on = total_on/total,
              share_off = total_off/total,
              share_abs = total/count
    ) %>% 
    return(twoyears)
}

# plot parameter
labelx = 0.705

pdf("Zenodo/6.2_plot.pdf",10,13)
par(mar=c(2.6,4.3,1,2.7),mfcol= c(4,2),oma = c(0,0,1.3,0))
absolute <- gdp_abs_en(~GDP_gr56,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21,31)], labels = absolute$year[c(1,11,21,31)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
#text(coord[1]-2,1.35,"PWT 5.6 vs 10.0",pos=4)
mtext(side = 3, adj = -.12, "a)", font = 2, line = 0.1)
mtext(side = 3, line = -0.5, "PWT 5.6 vs 10.0",cex = labelx, adj = 0.03)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)+40),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)
legend("topleft",legend=c("Absolute decoupling in late but not early vintage","Absolute decoupling in early but not late vintage","Share of countries absolutely decoupling in early vintage","Country total (right y-axis)"),col=c(viridis(6)[c(5,3)],1,1),pch=c(15,15,NA,NA),lty = c(NA,NA,1,2),pt.cex = c(1.5,1.5,1,1),bty="n",cex = 0.97,title = "Legend for all plots")


absolute <- gdp_abs_en(~GDP_gr61,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21,31,41)], labels = absolute$year[c(1,11,21,31,41)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 6.1 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "b)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

absolute <- gdp_abs_en(~GDP_gr70,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21,31,41)], labels = absolute$year[c(1,11,21,31,41)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 7.0 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "c)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

absolute <- gdp_abs_en(~GDP_gr90,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21,31,41,51)], labels = absolute$year[c(1,11,21,31,41,51)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 9.0 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "d)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

########################## MATERIALS ##

absolute <- gdp_abs(~GDP_gr56,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11)], labels = absolute$year[c(1,11)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 5.6 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "e)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = .9, "Country count",cex = labelx)

absolute <- gdp_abs(~GDP_gr61,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21)], labels = absolute$year[c(1,11,21)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 6.1 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "f)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

absolute <- gdp_abs(~GDP_gr70,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F)#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.6) 
axis(side = 1, at = coord[c(1,11,21)], labels = absolute$year[c(1,11,21)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -0.5, "PWT 7.0 vs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "g)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

absolute <- gdp_abs(~GDP_gr90,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10)
coord = barplot(as.matrix(t(absolute[,7:6])), beside=F,plot=F)
barplot(as.matrix(t(absolute[,7:6])), beside=F,col=viridis(6)[c(3,5)], ann = F, axes = F,ylim = c(0,.333))#names.arg = absolute$year,
lines(coord,absolute$share_abs)
axis(side = 2, at = pretty(range(c(0,max(absolute$share_on+absolute$share_off,na.rm=T)))),tcl = 0.25,las=2,hadj=0.55) 
axis(side = 1, at = coord[c(1,11,21,31)], labels = absolute$year[c(1,11,21,31)],padj=-1.4,lwd=0)
mtext(side = 1, line = 1.3, "Year",cex = labelx)
mtext(side = 2, line = 2.1, "Share of absolutely decoupling countries  (bars)\nof all countries (line)",cex = labelx)
mtext(side = 3, line = -1.5, "PWT 9.0\nvs 10.0",cex = labelx, adj = 0.03)
mtext(side = 3, adj = -.12, "h)", font = 2, line = 0.1)
par(new = TRUE)                        
plot(coord, absolute$count, type = "l", lty = 2, axes = F, ann = F,ylim = c(0,max(absolute$count)),xlim = c(coord[1]-.5,rev(coord)[1]+0.5))
axis(side = 4, at = pretty(range(c(0,max(absolute$count)))),tcl = 0.25,padj = -1.5) 
mtext(side = 4, line = 0.9, "Country count",cex = labelx)

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()

# generate datapoints file
fig6_data <- gdp_abs_en(~GDP_gr56,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6a") %>% select(-total,-total_on,-total_off)
fig_add <-gdp_abs_en(~GDP_gr61,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6b") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs_en(~GDP_gr70,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6c") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs_en(~GDP_gr90,~GDP_gr100, ~ener_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6d") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs(~GDP_gr56,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6e") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs(~GDP_gr61,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6f") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs(~GDP_gr70,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6g") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
fig_add <-gdp_abs(~GDP_gr90,~GDP_gr100, ~DMC_gr_10) %>% 
  filter(count>10) %>% mutate(figure_panel = "Figure 6h") %>% select(-total,-total_on,-total_off)
fig6_data <- bind_rows(fig6_data, fig_add)
write_csv(fig6_data,"Zenodo/6.3_fig6_datapoints.csv")