# Figure 1 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(viridis)
a <- read_csv("Zenodo/1.2_fig1a-sourcedata.csv") #set your working directory here
a <- a %>% mutate(excess = (gdp18-gdp70)/gdp70)

pdf("Zenodo/1.4_fig1a_plot.pdf",5.7,5.2)
par(mar=c(4,4,0.3,4.5))
plot(a$yr,a$gdp70,type="o",pch=19,ylim=c(min(a$gdp18),max(a$gdp18)),xlab = "Year",ylab="Index, 1929 = 1",col=viridis(4)[1]) # frame.plot=F,
lines(a$yr,a$gdp18,type="o",pch=17,col=viridis(4)[2])
par(new = TRUE)
plot(a$yr, a$excess, type="o",pch=15,col=viridis(4)[3], axes = FALSE, bty = "n", xlab = "", ylab = "")
axis(side=4, at = pretty(c(0,0.4)), labels = paste0(pretty(c(0,0.4))*100,"%"),las = 2,hadj=0.3,tcl=-0.2)
mtext(side = 4, "Excess of 2018 GDP over 1975 GDP",line=2.6)
mtext(side = 3, "a)",font =2,adj=-.2,line = -0.7)
legend("topleft",legend = c("US GDP index (1980 revision)","US GDP index (2018 revision)","Excess 2018 over 1980 GDP\nin per cent (right hand y-axis)"),bty="n",lty=1,pch=c(19,17,15),col=viridis(4)[1:3])
dev.off()

b <- read_csv("Zenodo/1.3_fig1b-sourcedata.csv")
b <- b %>% mutate(across(2:7, ~ .x/.x[1]),
                  across(2:7, ~ (log(.x)-log(lag(.x,10)))/10,.names = "gr_{.col}")
)

pdf("Zenodo/1.5_fig1b_plot.pdf",5.7,5.2)
par(mar=c(4,4,0.3,.3))
barplot(t(b[c(12,22,32,42,52,62,72,82),8:13]), names.arg = c("1930s","1940s","1950s","1960s","1970s","1980s","1990s","2000s"), beside = T,col=viridis(6),xlab = "Decade",cex.names = 0.9)
mtext(side = 2, line = 2.5, "Average annual GNP growth rate")
mtext(side = 3, "b)",font =2,adj=-.16,line = -0.7)
legend("topright",legend=c("1965","1980","1987","1999","2008","2018"), title="GNP revision year",ncol=2,pch=c(rep(15,6)),col=viridis(6),lwd=1,lty=0,bty="n",pt.cex=c(2))
dev.off()

#generate datapoints files
write_csv(a,"Zenodo/1.6_fig1a_datapoints.csv")

bout <- b[c(12,22,32,42,52,62,72,82),8:13]
bout <-cbind(c("1930s","1940s","1950s","1960s","1970s","1980s","1990s","2000s"),bout)
names(bout)[1]<-"Decade"
write_csv(bout,"Zenodo/1.7_fig1b_datapoints.csv")