####################
### INTRODUCTION ###
####################
# This script creates the figures and statitics used in materials and methods
# It requires the following packages: dplyr,ggplot2, tidyr, joyn
library(dplyr)
library(ggplot2)
rm(list=ls())

#################
### FIGURE S5 ###
#################

# To generate these numbers, do the following 14 steps
# 1. Run files 1a-2e in order.
# 2. In 1g-welfare.R change line 44 
     # from
     # growth <- growth |> filter(scenario=="med") |> select(-scenario)
     # to 
     # growth <- growth |> filter(scenario=="low") |> select(-scenario)
# 3. Run 1g-welfare.R
# 4. In 2a-damages_historical.R, change line 211
     # from
     # save(ea_gas,         file="03-outputdata/ea_gas.Rda")
     # to
     # save(ea_gas,         file="03-outputdata/ea_gas_low.Rda")
# 5. Run 2a-damages_historical.R
# 6. In 1g-welfare.R change line 44 
     # from
     # growth <- growth |> filter(scenario=="low") |> select(-scenario)
     # to 
     # growth <- growth |> filter(scenario=="hig") |> select(-scenario)
# 7. Run 1g-welfare.R
# 8. In 2a-damages_historical.R, change line 211
     # from
     # save(ea_gas,         file="03-outputdata/ea_gas_low.Rda")
     # to
     # save(ea_gas,         file="03-outputdata/ea_gas_hig.Rda")
# 9. Run 2a-damages_historical.R
# 10. In 1g-welfare.R change line 44 
     # from
     # growth <- growth |> filter(scenario=="hig") |> select(-scenario)
     # to 
     # growth <- growth |> filter(scenario=="med") |> select(-scenario)
# 11. Run 1g-welfare.R
# 12. In 2a-damages_historical.R, change line 211
     # from
     # save(ea_gas,         file="03-outputdata/ea_gas_hig.Rda")
     # to
     # save(ea_gas,         file="03-outputdata/ea_gas.Rda")
# 13. Run 2a-damages_historical.R
# 14. Run lines 54-108 of this script

rm(list=ls())
load("02-intermediatedata/ghg.Rda")
co2 <- ghg |>
  group_by(year) |>
  summarize(co2 = sum(co2,na.rm=TRUE)) |>
  ungroup() |>
  filter(year<=2024)

load("03-outputdata/ea_gas_low.Rda")
ea_gas_low <- ea_gas |>
              tidyr::pivot_longer(-c(gas,year),values_to="value",names_to="type")|>
              tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
              filter(outcome=="poor" & type=="ea" & gas=="co2" & discount=="df0") |>
              mutate(growth = "low")


load("03-outputdata/ea_gas_hig.Rda")
ea_gas_hig <- ea_gas |>
               tidyr::pivot_longer(-c(gas,year),values_to="value",names_to="type")|>
               tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
               filter(outcome=="poor" & type=="ea" & gas=="co2" & discount=="df0") |>
               mutate(growth = "high")

load("03-outputdata/ea_gas.Rda")
ea_gas      <- ea_gas |>
  tidyr::pivot_longer(-c(gas,year),values_to="value",names_to="type")|>
  tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
  filter(outcome=="poor" & type=="ea" & gas=="co2" & discount=="df0") |>
  mutate(growth = "baseline")

gamma <- rbind(ea_gas,ea_gas_low,ea_gas_hig) |>
  joyn::joyn(co2,match_type="m:1",by=c("year"),reportvar=FALSE,verbose=FALSE) |>
  mutate(gamma=value/co2) |>
  select(year,measure,gamma,growth) |>
  mutate(measure = recode(measure, 
                          "300" = "Years in\nextreme poverty\nper ktCO2", 
                          "830" = "Years in\nmoderate poverty\nper ktCO2",
                          "npl"="Years in\nnationally-defined\npoverty per ktCO2",
                          "gap"="Increase in\nprosperity gap\nper ktCO2",
                          "mea"="Income loss ($)\nper tCO2"),
         gamma = if_else(measure=="Income loss ($)\nper tCO2",-gamma*365/1000,gamma)) |>
        mutate(measure = factor(measure, levels=c("Years in\nextreme poverty\nper ktCO2",
                                             "Years in\nmoderate poverty\nper ktCO2",
                                              "Years in\nnationally-defined\npoverty per ktCO2",
                                               "Increase in\nprosperity gap\nper ktCO2",
                                               "Income loss ($)\nper tCO2")))



ggplot(gamma, aes(x=year,y=gamma,group=growth,color=as.factor(growth))) + 
  geom_line() + facet_wrap(~measure, scales = "free",nrow=1) + ylim(0,NA) +  theme_minimal() + 
  labs(x="",y="", fill="") +
  theme(legend.position="bottom", plot.title = element_text(hjust = 0.5), legend.title=element_blank(), legend.margin=margin(-10, 0, 0, 0)) +
  scale_color_manual(values=c("#56B4E9","#E69F00","darkgreen"),
                     labels = c("baseline", "high", "low"))
ggsave("05-Figures/FigS5.jpg", width = 7, height = 3)

# Change in gamma in 2024 with low growth
(gamma <- gamma |> 
         filter(year==2024 & growth!="high") |>
         filter(measure %in% c("Income loss ($)\nper tCO2","Years in\nextreme poverty\nper ktCO2")) |>
         group_by(measure) |>
         arrange(growth) |>
         summarize(impact = 100*(gamma[2]/gamma[1]-1)))

#################
### FIGURE S6 ###
#################
load("02-Intermediatedata/tempincrease_global.Rda")
figure <- tempincrease_global |> filter(year==2024) |> select(-tempincrease_ghg,-year) |>
          rename(CO2 = tempincrease_co2, CH4 = tempincrease_ch4, N2O = tempincrease_n2o) |>
          tidyr::pivot_longer(-c(impactyear,scenario), names_to="gas",values_to="tempincrease") |>
  mutate(gas = factor(gas, levels = c("CO2", "CH4", "N2O"), # The internal factor levels 
               labels = c(expression(CO[2]),expression(CH[4]), expression(N[2]*O))))
         
ggplot(figure) + geom_line(aes(x=impactyear,y=tempincrease)) + facet_wrap(~gas,labeller = label_parsed) + theme_minimal() +
  ylab("Global temperature increase (\u00B0C)") + xlab("") +  scale_x_continuous(
    breaks = c(2024, 2040, 2060, 2080))
ggsave("05-Figures/FigS6.jpg", width = 8, height = 3.5)

###################################
### INEQUALITY B10/T10 IMPACSTS ###
###################################

# Load relevant files
load("02-Intermediatedata/welfare.Rda") 
welfare <- welfare |> filter(year==2024) |>
  group_by(code) |> 
  mutate(welf_share = welf/sum(welf),
         damage_share = welf_share^0.64,
         damage_share = damage_share/sum(damage_share),
         welfmean_cf  = mean(welf),
         welf_cf     =  welf-1000*welfmean_cf*0.01*damage_share) |>
  ungroup() |>
  group_by(code) |>
  summarise(
    threshold_10 = quantile(welf, 0.10, na.rm = TRUE),
    threshold_90 = quantile(welf, 0.90, na.rm = TRUE),
    welf_b10 = mean(welf[welf <= threshold_10], na.rm = TRUE),
    welf_t10 = mean(welf[welf >= threshold_90], na.rm = TRUE),
    welf_cf_b10 = mean(welf_cf[welf <= threshold_10], na.rm = TRUE),
    welf_cf_t10 = mean(welf_cf[welf >= threshold_90], na.rm = TRUE)) |>
  ungroup() |>
  mutate(damage_ratio_t10 = (1-welf_cf_t10/welf_t10)*100,
         damage_ratio_b10 = (1-welf_cf_b10/welf_b10)*100,
         ratio_t10_b10 = damage_ratio_b10/damage_ratio_t10)

summary(welfare$ratio_t10_b10)  

#################
### FIGURE S7 ###
#################
rm(list=ls())
load("03-outputdata/damage_annual.Rda")
damage_annual <- damage_annual |>
  select(-year,-responsibilityshare_ea,-scenario_ea) |>
  rename(year=impactyear_ea) |>
  tidyr::pivot_longer(-c(gas,year),values_to="value",names_to="type") |>
  tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
  filter(discount=="df0") |>
  select(year,gas,measure,value) |>
  mutate(value = if_else(measure=="mea",-365*value/10^9,value/10^6),
         measure = if_else(measure=="300","Million years in\nextreme poverty",measure),
         measure = if_else(measure=="830","Million years in\nmoderate poverty",measure),
         measure = if_else(measure=="npl","Million years in\nnationally-defined\npoverty",measure),
         measure = if_else(measure=="gap","Million decrease in\nprosperity gap",measure),
         measure = if_else(measure=="mea","Billion dollars\nlost",measure),
         measure = factor(measure, levels=c("Million years in\nextreme poverty",
                                            "Million years in\nmoderate poverty",
                                            "Million years in\nnationally-defined\npoverty",
                                            "Million decrease in\nprosperity gap",
                                            "Billion dollars\nlost")))

ggplot(damage_annual, aes(x = year, y = value, fill = gas)) +
  stat_smooth(
    geom = "area",
    method = "loess",  # Choose smoothing method, e.g., "loess" or "gam"
    span = 0.7,        # Adjust the span (e.g., 0.5) to control smoothing (for loess)
    se = FALSE,        # Remove the confidence band
    position = "stack" # Use "stack" if you want a stacked area plot for the smoothed values
  ) +
  facet_wrap(~measure, scale = "free",nrow=1) +
  labs(x = "", y = "",fill="") +
  theme_minimal() +
  scale_fill_manual(
    values = c("darkgreen", "#56B4E9", "#E69F00"),
    labels = c(expression(CH[4]), expression(CO[2]), expression(N[2]*O))
  ) + scale_x_continuous(
    breaks = c(2025,2050,2075,2100)) +  theme(legend.position="bottom",legend.margin=margin(-10, 0, 0, 0)) 
ggsave("05-Figures/FigS7.jpg", width = 7, height = 3)

#################
### FIGURE S8 ###
#################

load("03-outputdata/gamma.Rda")

gamma <- gamma |>  
  filter(gas =="co2") |>
  mutate(measure = recode(measure, "300" = "Years in\nextreme poverty\nper kt", 
                          "830" = "Years in\nmoderate poverty\nper kt",
                          "npl"="Years in\nnationally-defined\npoverty per kt",
                          "gap"="Increase in\nprosperity gap\nper kt",
                          "mea"="Income loss ($)\nper t"),
         gas = toupper(gas)) |>
  mutate(gamma = if_else(measure=="Income loss ($)\nper t",-gamma*365/1000,gamma)) |>
  tidyr::unite(combine,measure,gas, sep = "",remove=FALSE) |>
  mutate(across(combine, ~factor(., levels=c("Years in\nextreme poverty\nper ktCO2",
                                             "Years in\nmoderate poverty\nper ktCO2",
                                             "Years in\nnationally-defined\npoverty per ktCO2",
                                             "Increase in\nprosperity gap\nper ktCO2",
                                             "Income loss ($)\nper tCO2")))) 

ggplot(gamma, aes(x=year,y=gamma,group=discount,color=as.factor(discount))) + 
  geom_line() + facet_wrap(~combine, scales = "free",nrow=1) + ylim(0,NA) +  theme_minimal() + 
  labs(x="",y="", fill="") +
  theme(legend.position="bottom", plot.title = element_text(hjust = 0.5), legend.title=element_blank(),legend.margin=margin(-10, 0, 0, 0)) +
  scale_color_manual(values=c("#56B4E9","#E69F00","darkgreen"), 
                     labels=c("0%","2%","4%")) 
ggsave("05-Figures/FigS8.jpg", width = 7, height = 3)

###############################
### CUMULATIVE LOSS BY 2030 ###
###############################
load("03-outputdata/damage_2030.Rda")
sum(damage_2030$added_poor)/10^6

#################
### FIGURE S3 ###
#################
load("02-intermediatedata/elasticities.Rda")

# Winsorization 
load("02-intermediatedata/elasticities.Rda")
(elasticities |>
  group_by(gas) |>
  summarize(min = round(min(elasticity, na.rm = TRUE),2),
            max = round(max(elasticity, na.rm = TRUE),2)))

ggplot(elasticities) + geom_density(aes(x=elasticity_unwinsorized)) + facet_wrap(~gas,scales="free_y",labeller = label_parsed) + xlim(-0.4,1.4) + theme_minimal() + 
  geom_vline(xintercept=0,linetype="dashed") + xlab("Elasticity") +   geom_vline(xintercept=1,linetype="dashed") 
ggsave("05-Figures/FigS3.jpg", width = 8, height = 3.5)

# Table S1 regression output
# See 1k-elasticities.R
