###################
### INTODUCTION ###
###################
# This script plots figures with the results and produces other key statistics mentioned in the main body of the paper
# It also produces some annex figures whenever these use nearly the exact same code as a figure in the main paper.
# It requires the following packages: dplyr, ggplot, ggalluvial, tidyr, joyn, forcats, stringr, scales
library(dplyr)
library(ggalluvial)
library(ggplot2)
rm(list=ls())

######################
### FIGURE 1A & 1B ###
######################

for (group in c("region","incgroup")) {   
  # Class
  load("02-intermediatedata/class.Rda")  
  class <- class |> rename("group" = group)
  
  # By emitter
  load("03-outputdata/welfare_cf_country.Rda")
  emitter <- welfare_cf_country |>
    tidyr::pivot_longer(-c(code,year,pop),values_to="value",names_to="type")|>
    tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
    filter(discount=="df0" & year==2024 & type=="ea" & outcome=="poor") |>
    select(-c(outcome,pop,type,discount,year)) |>
    mutate(measure = recode(measure, "300" = "$3.00", "830" = "$8.30","npl"="National poverty line","gap"="Prosperity gap","mea"="Income loss")) |>
    rename("emitter" = "value") |>
    joyn::joyn(class,by="code",match_type="m:1",reportvar=FALSE) |>
    group_by(group,measure) |>
    summarize(emitter = sum(emitter,na.rm=TRUE)/10^6) |>
    rename("group_emitter" = "group") |>
    group_by(measure) |>
    mutate(damageshare = emitter/sum(emitter)) |>
    ungroup()
  
  # By subject
  load("03-outputdata/damage_country.Rda")
  subject <- damage_country |>
    tidyr::pivot_longer(-c(code,year,pop),values_to="value",names_to="type")|>
    tidyr::separate(type,c("outcome","measure","type","discount"),sep="_") |>
    filter(year==2024 & discount=="df0" & outcome=="poor") |>
    rename(subject = value) |>
    select(code,year,subject,measure) |>
    mutate(measure = recode(measure, "300" = "$3.00", "830" = "$8.30","npl"="National poverty line","gap"="Prosperity gap","mea"="Income loss")) |>
    joyn::joyn(class,by="code",match_type="m:1",reportvar=FALSE) |>
    group_by(group,measure) |>
    summarize(subject = sum(subject,na.rm=TRUE)/10^6) |>
    rename("group_subject" = "group") |>
    ungroup()
  
  rm(welfare_cf_country,damage_country)  
  
  alluvial <- expand.grid(group_emitter=unique(emitter$group_emitter),group_subject=unique(subject$group_subject),measure=unique(emitter$measure)) |>
    joyn::joyn(emitter,by=c("measure","group_emitter"),match_type="m:1",reportvar=FALSE) |>
    joyn::joyn(subject,by=c("measure","group_subject"),match_type="m:1",reportvar=FALSE) |>
    filter(measure=="$3.00") |>
    mutate(value=damageshare*subject) |>
    select(-c(damageshare,subject,emitter)) |>
    arrange(group_emitter,group_subject) |>
    mutate(id=row_number())
  
  alluvial1 <- alluvial |> 
    select(-group_subject) |>
    rename("group" = "group_emitter") |>
    mutate(type="Caused in 2024")
  
  alluvial2 <- alluvial |> 
    select(-group_emitter) |>
    rename("group" = "group_subject") |>
    mutate(type="To be experienced")
  
  alluvial <- rbind(alluvial1,alluvial2)
  
  if (group=="incgroup") {
    alluvial$group <- factor(alluvial$group, levels = c("High income", "Upper middle income", "Lower middle income","Low income"))
  }
  if (group=="region") {
    alluvial <- alluvial |> mutate(group = if_else(group=="Middle East, North Africa, Afghanistan & Pakistan","Middle East & North Africa",group))
    alluvial$group <- factor(alluvial$group, levels = c("East Asia & Pacific", "Europe & Central Asia", "North America","South Asia",
                                                        "Latin America & Caribbean","Middle East & North Africa","Sub-Saharan Africa"))
    
  }
  
  if (group=="incgroup") {
    ggplot(alluvial,
           aes(x = type, stratum = group, alluvium = id, y=value, fill = group)) +
      scale_x_discrete(expand = c(.1, .1)) +
      geom_flow() +
      geom_stratum(alpha = .5) +
      geom_text(stat = "stratum",
                aes(label = after_stat(round(count))))  + theme_minimal() +
      ylab("Million years in extreme poverty") + xlab("") +
      theme(legend.title=element_blank(),axis.text.x = element_text(size = 10)) +
      scale_fill_manual(values = c("#56B4E9","darkgreen","#E69F00","darkred")) +
      guides(fill=guide_legend(nrow=4,byrow=TRUE))
    
    ggsave("05-Figures/Fig1a.jpg", width = 5, height = 4.5)
  }
  if (group=="region") {
    ggplot(alluvial,
           aes(x = type, stratum = group, alluvium = id, y=value, fill = group)) +
      scale_x_discrete(expand = c(.1, .1)) +
      geom_flow() +
      geom_stratum(alpha = .5) +
      geom_text(stat = "stratum",
                aes(label = after_stat(round(count))))  + theme_minimal() +
      ylab("Million years in extreme poverty") + xlab("") +
      theme(legend.title=element_blank(),axis.text.x = element_text(size = 10)) +
      scale_fill_manual(values = c("black","darkblue","#56B4E9","darkgreen","#E69F00","grey","darkred")) +
      guides(fill=guide_legend(nrow=7,byrow=TRUE))
    ggsave("05-Figures/Fig1b.jpg", width = 5, height =4.5)
  }
}

# Total future poverty caused by 2024 emissions in millions
(emitter |> filter(measure=="$3.00") |> summarize(emitter=sum(emitter)))

# Millions of extreme poor in 2024
load("02-intermediatedata/welfare_global.Rda")
(welfare_global <- welfare_global |> filter(year==2024 & region=="World") |> select(poor_300_ac) |> mutate(poor_300_ac = poor_300_ac/10^6))

# Damage experience stats with total dollar loss
(subject <- subject |>
    filter(measure=="Income loss") |>
    mutate(share = subject/sum(subject)*100))

# Share from avarage growth shocks and share from the distribution of these shocks
# To generate these numbers, do the following 7 steps
# 1. Run files 1a-2e in order, and everything of 3a-main.R above line 128. 
# 2. In 2a-damages_historical.R, replace line 55 
      # from
      # for(i in 1990:2024) {
      # to  
      # for(i in 2024:2024) {
# 3. In 2a-damages_historical.R, add an "#" before lines 78, 80, 81, and 83
# 4. In 2a-damages_historical.R, add an "#" before lines 173, 208, 209, 210 and 211
# 5. Run 2a-damages_historical.R. It will return an error, that's ok.
# 6. When the code has run, execute the line below, without the "#" 
# sum(damage_country_temp$poor_300_da_df0)/10^6
# 7. Undo steps 2, 3, and 4
 
######################
### FIGURE 1C & S1 ###
######################

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

gamma <- gamma |>  
  filter(discount=="0") |>
  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",
                                             "Years in\nextreme poverty\nper ktCH4",
                                             "Years in\nmoderate poverty\nper ktCH4",
                                             "Years in\nnationally-defined\npoverty per ktCH4",
                                             "Increase in\nprosperity gap\nper ktCH4",
                                             "Income loss ($)\nper tCH4",
                                             "Years in\nextreme poverty\nper ktN2O",
                                             "Years in\nmoderate poverty\nper ktN2O",
                                             "Years in\nnationally-defined\npoverty per ktN2O",
                                             "Increase in\nprosperity gap\nper ktN2O",
                                             "Income loss ($)\nper tN2O"))))


ggplot(gamma[gamma$gas=="CO2",], aes(x=year,y=gamma)) + 
  geom_line() + facet_wrap(~combine, scales = "free",nrow=1) + ylim(0,NA) +  theme_minimal() + 
  labs(x="",y="", fill="") +
  theme(legend.position="bottom") +    theme(plot.title = element_text(hjust = 0.5))
ggsave("05-Figures/Fig1c.jpg", width = 7, height =3)

ggplot(gamma[gamma$gas!="CO2",], aes(x=year,y=gamma)) + 
  geom_line() + facet_wrap(~combine, scales = "free",nrow=2) + ylim(0,NA) +  theme_minimal() + 
  labs(x="",y="", fill="") +
  theme(legend.position="bottom") +    theme(plot.title = element_text(hjust = 0.5))
ggsave("05-Figures/FigS1.jpg", width = 8, height = 6)

# Gamma in 2024
load("03-outputdata/gamma.Rda")
(gamma <- gamma |>
         filter(gas=="co2" & year=="2024" & discount=="0"))

#################
### FIGURE 2A ###
#################
load("03-outputdata/welfare_cf_global.Rda")

welfare_cf_global <- welfare_cf_global |>
  tidyr::pivot_longer(-c(region,year,pop),values_to="value",names_to="type") |>
  tidyr::separate(type,c("outcome","discount","measure","type"),sep="_")  |>
  filter(discount %in% c("ndf","df0") & type!="da" & measure=="300" & outcome=="rate") |> 
  mutate(measure = recode(measure, "300" = "$3.00")) |>
  filter(region %in% c("World","East Asia & Pacific","North America","Sub-Saharan Africa")) |>
  mutate(type = factor(type, levels = c("ea","ac","cf"))) 

ggplot() + 
  geom_area(welfare_cf_global[welfare_cf_global$type!="cf",],mapping = aes(x=year,y=value,fill=type)) + facet_wrap(~forcats::fct_rev(region),ncol=4) + 
  labs(x="",y="%", fill="") + theme_minimal() +   
  scale_fill_manual(values=c("darkgreen","#56B4E9"),labels=c("Future poverty caused by emissions","Domestic poverty rate")) + 
  theme(legend.position="bottom",legend.title=element_blank(),plot.title = element_text(hjust = 0.5), legend.margin = margin(t = -10, r = 0, b = 0, l = 0, unit = "pt")) +
  geom_line(welfare_cf_global[welfare_cf_global$type=="cf",],mapping = aes(x=year,y=value,color=type),linetype="dashed")  +
  scale_color_manual(values=c("#E69F00"), labels=c("Emissions-adjusted poverty rate")) + 
  guides(fill=guide_legend(reverse = TRUE)) 
ggsave("05-Figures/Fig2a.jpg", width = 7, height = 3)

#################
### FIGURE 2B ###
#################

load("03-outputdata/welfare_cf_country.Rda")
welfare_cf_country <- welfare_cf_country |>
  tidyr::pivot_longer(-c(code,year,pop),values_to="value",names_to="type") |>
  tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
  filter(discount %in% c("ndf","df0")) |>
  filter(measure==300 & year==2024 & outcome=="rate") |>
  select(-discount,-year,-measure,-outcome,-pop) |>
  filter(code %in% c("QAT","GAB","BGD")) |>
  mutate(code = if_else(code=="BGD","Bangladesh",if_else(code=="GAB","Gabon","Qatar"))) |>
  mutate(type = factor(type, levels = c("ea","ac","cf"))) 

ggplot() + 
  geom_bar(welfare_cf_country[welfare_cf_country$type!="cf",], mapping=aes(fill=type, y=value, x=code),position="stack", stat="identity") + 
  scale_fill_manual(values=c("darkgreen","#56B4E9"), 
                    labels=c("Future poverty caused by emissions","Domestic poverty rate")) + 
  theme_minimal() +
  geom_text(welfare_cf_country[welfare_cf_country$type!="cf",], mapping=aes(y=value, x=code,label = paste0(round(value,1), "%")), 
            position = position_stack(vjust = 0.5), 
            color = "white", size = 3) + ylab("%")+ xlab("") + 
  theme(legend.position="bottom",legend.box = "vertical",rows=3,legend.title=element_blank(),plot.title = element_text(hjust = 0.5),legend.margin = margin(t = -10, r = 0, b = 0, l = 0, unit = "pt"))  + 
  guides(fill=guide_legend(reverse = TRUE,nrow=2,byrow=TRUE)) +
  geom_point(welfare_cf_country[welfare_cf_country$type=="cf",],mapping = aes(x=code,y=value,color=type),size=3) +
  scale_color_manual(values=c("#E69F00"), labels=c("Emissions-adjusted poverty rate")) +
  geom_text(welfare_cf_country[welfare_cf_country$type=="cf",], mapping=aes(y=value, x=code,label = paste0(round(value,1), "%")), 
            vjust = -0.5,color = "#E69F00", size = 3) +
  coord_cartesian(clip = "off")
ggsave("05-Figures/Fig2b.jpg", width = 3, height = 4)


######################
### FIGURE 2C & S2 ###
######################
load("02-intermediatedata/gdp.Rda")
gdp <- gdp |> filter(year==2024)

load("02-intermediatedata/class.Rda")
load("03-outputdata/welfare_cf_country.Rda")
welfare_cf_country <- welfare_cf_country |>
  tidyr::pivot_longer(-c(code,year,pop),values_to="value",names_to="type")|>
  tidyr::separate(type,c("outcome","discount","measure","type"),sep="_") |>
  filter(discount %in% c("ndf","df0") & type!="da") |>
  select(-discount) |> 
  tidyr::pivot_wider(names_from=c(outcome,type)) |>
  mutate(measure = recode(measure, "300" = "$3.00", "830" = "Moderate povety","npl"="Poverty at national lines","gap"="Prosperity gap","mea" = "Income loss")) |>
  filter(year==2024) |>
  joyn::joyn(gdp,by="code",match_type="m:1",reportvar=FALSE,verbose=FALSE) |>
  joyn::joyn(class,by="code",match_type="m:1",reportvar=FALSE,verbose=FALSE)|>
  select(gdppc,region,measure,rate_cf,rate_ac,rate_ea,code,incgroup) |>
  mutate(rate_ea = if_else(measure=="Income loss",-rate_ea,rate_ea),
         ratio = rate_ea/rate_ac, 
         ratio = if_else(ratio>100,NA,ratio)) |>
  mutate(incgroup = stringr::str_remove(incgroup," income"))

ggplot(data=welfare_cf_country[welfare_cf_country$measure=="$3.00",],aes(x=gdppc,y=ratio)) + 
  geom_point(size=0) +
  geom_text(aes(color=incgroup),label=welfare_cf_country$code[welfare_cf_country$measure=="$3.00"],size=3) +
  geom_smooth(color="black",se=F,linewidth=1.2) + 
  xlab("GDP/capita (2021 USD PPP)") + ylab("Future global poverty caused by emissions\nrelative to domestic poverty") + 
  scale_x_continuous(trans='log10', labels = scales::comma,lim=c(800,82000)) + theme_minimal() +
  scale_y_continuous(trans='log10', labels = scales::comma) + 
  scale_color_manual(values = c("darkred","#E69F00","darkgreen","#56B4E9"), breaks = c("Low","Lower middle","Upper middle","High")) + 
  theme(legend.position="bottom",plot.title = element_text(hjust = 0.5),legend.spacing = unit(0, "pt"),plot.margin = unit(c(10, 20, 10, 10), "pt")) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", size=1) +
  guides(color=guide_legend(nrow=1,title="Income group:",byrow=TRUE))
ggsave("05-Figures/Fig2c.jpg", width = 5.5, height = 4)

ggplot(data=welfare_cf_country[welfare_cf_country$measure!="$3.00",],aes(x=gdppc,y=ratio)) + 
  geom_point(size=0) + facet_wrap(~measure) +
  geom_text(aes(color=incgroup),label=welfare_cf_country$code[welfare_cf_country$measure!="$3.00"],size=3) +
  geom_smooth(color="black",se=F,linewidth=1.2) + 
  xlab("GDP/capita (2021 USD PPP)") + ylab("Future global deprivation by emissions\nrelative to domestic 2024 deprivation") + 
  scale_x_continuous(trans='log10', labels = scales::comma,lim=c(800,82000)) + theme_minimal() +
  scale_y_continuous(trans='log10', labels = scales::comma) + 
  scale_color_manual(values = c("darkred","#E69F00","darkgreen","#56B4E9"), breaks = c("Low","Lower middle","Upper middle","High")) + 
  theme(legend.position="bottom",plot.title = element_text(hjust = 0.5),legend.spacing = unit(0, "pt"),plot.margin = unit(c(10, 20, 10, 10), "pt")) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", size=1) +
  guides(color=guide_legend(nrow=1,title="Income group:",byrow=TRUE))
ggsave("05-Figures/FigS2.jpg", width = 8, height = 6)

#####################
### FIGURE 3 & S4 ###
#####################
rm(list=ls())

## Median and interquartile range  of elasticities
load("02-intermediatedata/elasticities.Rda")
(elasticities <- elasticities |>
  group_by(gas) |>
  summarize(p25 = round(quantile(elasticity_unwinsorized, probs = 0.25, na.rm = TRUE),2),
            p50 = round(quantile(elasticity_unwinsorized, probs = 0.50, na.rm = TRUE),2),
            p75 = round(quantile(elasticity_unwinsorized, probs = 0.75, na.rm = TRUE),2)))

load("03-outputdata/cost_benefit.Rda")
cost_benefit <- cost_benefit |>
  mutate(incgroup = stringr::str_remove(incgroup," income"),
         incgroup = factor(incgroup, levels = c("Low", "Lower middle", "Upper middle","High")),
         measure = recode(measure, "300" = "$3.00", "830" = "$8.30","npl"="National poverty line","gap"="Prosperity gap","mea" = "Income loss")) 

# Figure 3a
ggplot(cost_benefit[cost_benefit$measure=="Prosperity gap" & cost_benefit$type=="ebre",],aes(x=gdppc)) + 
  geom_point(aes(y=factor),size=0) +
  geom_text(aes(y=factor,color=incgroup),label=cost_benefit$code[cost_benefit$measure=="Prosperity gap" & cost_benefit$type=="ebre"],size=3) +
  geom_smooth(aes(y=factor),color="black",linewidth=1.2) +
  xlab("GDP/capita") + ylab("Break-even emission-growth elasticity") + 
  scale_x_continuous(trans='log10', labels = scales::comma,lim=c(800,150000)) + 
  scale_y_continuous(trans='log10',labels = scales::comma) + theme_minimal() +
  scale_color_manual(values = c("#56B4E9","#E69F00","darkgreen","darkred")) +
  theme(legend.position="bottom",legend.title=element_blank(),plot.title = element_text(hjust = 0.5)) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", size=1)
ggsave("05-Figures/Fig3a.jpg", width = 4.5, height = 4.5)

# Figure 3b
table <- cost_benefit |>
  filter(measure %in% c("Prosperity gap","Income loss") & type!="ebre") |>
  mutate(lessthan1 = factor<1,
         lessthan10 = factor<10) |>
  filter(!is.na(gdppc) & factor!=-Inf) |>
  group_by(incgroup,measure,type) |>
  summarize(lessthan1  = 100*mean(lessthan1, na.rm=TRUE),
            lessthan10 = 100*mean(lessthan10, na.rm=TRUE)) |>
  ungroup() |>
  tidyr::pivot_longer(-c(incgroup,measure,type),names_to="factor",values_to="value") 

ggplot(table[table$measure=="Prosperity gap" & table$type=="ehis",],aes(x=incgroup,y=value,fill=factor)) +
  geom_bar(stat="identity",position="dodge") +
  scale_fill_manual(values=c("darkgreen","#56B4E9"),
                    labels=c(expression(elasticity > tilde(epsilon)),expression(elasticity > frac(tilde(epsilon),10)))) +
  theme_minimal() + ylab("Share of countries (%)") + ylim(0,100) +
  geom_text(table[table$measure=="Prosperity gap" & table$type=="ehis",],
            mapping=aes(y=value, x=incgroup,color=factor,label = paste0(round(value,1), "%")),
            position = position_dodge(width = .9),vjust=-0.3,show.legend = FALSE,size = 3) +
  scale_color_manual(values=c("darkgreen","#56B4E9")) +
  xlab("Income group") + theme(legend.position="bottom",legend.spacing = unit(0, "pt"),legend.title=element_blank(),plot.title = element_text(hjust = 0.5))  +
  guides(fill=guide_legend(nrow=1,byrow=TRUE))
ggsave("05-Figures/Fig3b.jpg", width = 4, height = 4)

# Figure S4a
ggplot(cost_benefit[cost_benefit$measure=="Income loss" & cost_benefit$type=="ebre",],aes(x=gdppc)) + 
  geom_point(aes(y=factor),size=0) +
  geom_text(aes(y=factor,color=incgroup),label=cost_benefit$code[cost_benefit$measure=="Income loss" & cost_benefit$type=="ebre"],size=3) +
  geom_smooth(aes(y=factor),color="black",linewidth=1.2) +
  xlab("GDP/capita") + ylab("Break-even emission-growth elasticity") + 
  scale_x_continuous(trans='log10', labels = scales::comma,lim=c(800,150000)) + 
  scale_y_continuous(trans='log10',labels = scales::comma) + theme_minimal() +
  scale_color_manual(values = c("#56B4E9","#E69F00","darkgreen","darkred")) +
  theme(legend.position="bottom",legend.title=element_blank(),plot.title = element_text(hjust = 0.5)) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", size=1)
ggsave("05-Figures/FigS4a.jpg", width = 4.5, height = 4.5)

# Figure S4b
ggplot(table[table$measure=="Income loss" & table$type=="ehis",],aes(x=incgroup,y=value,fill=factor)) + geom_bar(stat="identity",position="dodge") + 
  scale_fill_manual(values=c("darkgreen","#56B4E9"), 
                    labels=c(expression(elasticity > tilde(epsilon)),expression(elasticity > frac(tilde(epsilon),10)))) + 
  theme_minimal() + ylab("Share of countries (%)") + ylim(0,100) +
  geom_text(table[table$measure=="Income loss" & table$type=="ehis",], mapping=aes(y=value, x=incgroup,color=factor,label = paste0(round(value,1), "%")),position = position_dodge(width = .9),vjust=-0.3,show.legend = FALSE,size = 3) +
  scale_color_manual(values=c("darkgreen","#56B4E9")) +
  xlab("Income group") +  theme(legend.position="bottom",legend.spacing = unit(0, "pt"),legend.title=element_blank(),plot.title = element_text(hjust = 0.5))  + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE))
ggsave("05-Figures/FigS4b.jpg", width = 4, height = 4)

# Growth spillovers
(spillovers <- cost_benefit |>
    filter(measure %in% c("Prosperity gap","gap_spillover")) |>
    filter(code=="USA" & type=="ehis") |>
    tidyr::pivot_wider(names_from="measure",values_from="factor")) 

# Figure S4c
ggplot(table[table$measure=="Prosperity gap" & table$type=="emin",],aes(x=incgroup,y=value,fill=factor)) + 
  geom_bar(stat="identity",position="dodge") +
  scale_fill_manual(values=c("darkgreen","#56B4E9"), 
                    labels=c(expression(elasticity > tilde(epsilon)),expression(elasticity > frac(tilde(epsilon),10)))) + 
  theme_minimal() + ylab("Share of countries (%)") + ylim(0,100) +
  geom_text(table[table$measure=="Prosperity gap" & table$type=="emin",], mapping=aes(y=value, x=incgroup,color=factor,label = paste0(round(value,1), "%")),position = position_dodge(width = .9),vjust=-0.3,show.legend = FALSE,size = 3) +
  scale_color_manual(values=c("darkgreen","#56B4E9")) +
  xlab("Income group") +  theme(legend.position="bottom",legend.spacing = unit(0, "pt"),legend.title=element_blank(),plot.title = element_text(hjust = 0.5))  + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE))
ggsave("05-Figures/FigS4c.jpg", width = 4, height = 4)

# Damages under various RCP pathways
load("03-outputdata/damage_future.Rda")
(results <- damage_future |>
           group_by(scenario) |>
           summarize(poor = sum(poor_300_da_df0)/10^9) |>
           ungroup())