####################
### INTRODUCTION ###
####################
# This script calculates the welfare damages by emissions from 1990-2024. It takes several hours to run
# It requires the following packages: dplyr, joyn, stringr
library(dplyr)
rm(list=ls())

# Load relevant files
load("02-intermediatedata/welfare.Rda") 
load("02-intermediatedata/pop.Rda")
load("02-intermediatedata/welfare_country.Rda")
load("02-intermediatedata/precipitation_inequality.Rda")
load("02-intermediatedata/npl.Rda")
load("02-intermediatedata/responsibilityshare_country.Rda")
load("02-intermediatedata/responsibilityshare_gas.Rda")

#######################################
### CALCULATE WELFARE DAMAGE SHARES ###
#######################################

# Calculate welfare shares and mean
welfare <- welfare |>
group_by(code,year) |> 
mutate(welf_share = welf/sum(welf),
       damage_share = welf_share^0.64,
       damage_share = damage_share/sum(damage_share)) |> 
       ungroup() |>
select(-welf_share)

########################
### BOTTOM 50 SHARES ###
########################
b50 <- welfare |> 
       group_by(code,year) |>
       arrange(welf) |>
       mutate(b50 = row_number()<=500) |>
       ungroup() |>
       group_by(code,year,b50) |>
       summarize(welf = sum(welf)) |>
       group_by(code,year) |>
       mutate(b50share = welf/sum(welf)*100) |>
       filter(b50==TRUE) |>
       select(code,year,b50share) |>
       ungroup()
  
load("02-Intermediatedata/growthshock.Rda")
growthshock_full <- growthshock |>
                filter(scenario=="historical") |>
                select(-scenario)

#################################
### DISCOUNTED WELFARE DAMAGE ###
#################################
for(i in 1990:2024) {
print(i)
# Merge on relevant welfare shock
x     <- paste0("growthshock_tot",i)
x_tmu <- paste0("growthshock_tmu",i)
x_oth <- paste0("growthshock_oth",i)
growthshock <- growthshock_full |> 
                rename("growthshock"=x,"growthshock_tmu"=x_tmu,"growthshock_oth"=x_oth) |>
                select(code,year,growthshock,growthshock_tmu,growthshock_oth) |> 
                joyn::joyn(elasticity,by="code",          match_type="m:1",reportvar=FALSE,verbose=FALSE) |>
                joyn::joyn(b50,       by=c("code","year"),match_type="1:1",reportvar=FALSE,verbose=FALSE) |>
                mutate(b50c  = elasticity*growthshock_oth,
                       alpha = (b50c*b50share)/(50-b50share)/100) |>
                select(code,year,growthshock,growthshock_tmu,alpha)
# Merge data
temp <- welfare |>
        # Removes years not used, speeds up the code a bit
        filter(between(year,i,i+76)) |>
        joyn::joyn(growthshock, by=c("code","year"),match_type="m:1",reportvar=FALSE,verbose=FALSE,keep="inner") |>
        group_by(code,year) |>           
        # Calculate counterfactual welfare with no change to inequality
        mutate(welf_cf     = welf/(1+growthshock/100),
               # Calculate change in inequality from warming
                welfmean_cf = mean(welf_cf),
               # 1000*welfmean_cf*(growthshock) = total dollar loss 
                welf_cf     = if_else(growthshock_tmu<0,welf-1000*welfmean_cf*(growthshock_tmu/100)*damage_share,welf_cf),
                welf_cf     = if_else(welf_cf<0.28,0.28,welf_cf),
               # Calculate change in inequality from other channels 
                welf_cf     = (1+alpha)*welf_cf-alpha*welfmean_cf,
               welf_cf     = if_else(welf_cf<0.28,0.28,welf_cf)) |>
        select(welf_cf,code,year) |>
        ungroup() |>
        # Add on national poverty lines
        joyn::joyn(npl,match_type="m:1",by="code",reportvar=FALSE,keep="left",verbose=FALSE) |>
        # Calculate country-year level poverty and prosperity gap
        group_by(code,year) |>
        summarize(rate_300_cf = mean(welf_cf<3.00)*100,
                  rate_830_cf = mean(welf_cf<8.30)*100,
                  rate_npl_cf = mean(welf_cf<npl)*100,
                  rate_mea_cf = mean(welf_cf),
                  rate_gap_cf = mean(if_else(28/welf_cf<100,28/welf_cf,100))) |>
        ungroup() |>
        # Country-year annual counterfactual welfare statistics               
        joyn::joyn(pop,by=c("year","code"),match_type="1:1",reportvar=FALSE,verbose=FALSE,keep="inner") |>
        # Merge with actual welfare estimates |>
        joyn::joyn(welfare_country,by=c("code","year"),match_type="1:1",reportvar=FALSE,verbose=FALSE,keep="inner") |>
        mutate(add_300     = poor_300_ac-rate_300_cf/100*pop,
               add_830     = poor_830_ac-rate_830_cf/100*pop,
               add_npl     = poor_npl_ac-rate_npl_cf/100*pop,
               add_mea     = poor_mea_ac-rate_mea_cf*pop,
               add_gap     = poor_gap_ac-rate_gap_cf*pop,
               npv_df0_300 = add_300/(1.00)^(year-i),
               npv_df0_830 = add_830/(1.00)^(year-i),
               npv_df0_npl = add_npl/(1.00)^(year-i),
               npv_df0_mea = add_mea/(1.00)^(year-i),
               npv_df0_gap = add_gap/(1.00)^(year-i),
               npv_df2_300 = add_300/(1.02)^(year-i),
               npv_df2_830 = add_830/(1.02)^(year-i),
               npv_df2_npl = add_npl/(1.02)^(year-i),
               npv_df2_mea = add_mea/(1.02)^(year-i),
               npv_df2_gap = add_gap/(1.02)^(year-i),
               npv_df4_300 = add_300/(1.04)^(year-i),
               npv_df4_830 = add_830/(1.04)^(year-i),
               npv_df4_npl = add_npl/(1.04)^(year-i),
               npv_df4_mea = add_mea/(1.04)^(year-i),
               npv_df4_gap = add_gap/(1.04)^(year-i)) |>
        select(year,code,starts_with("npv")) 

# Aggregated to the country committing the damage
ea_country_temp <- temp |>
                   select(-contains("df2"),-contains("df4")) |>
                   group_by(year) |>
                   summarize(across(starts_with("npv"),~sum(.),.names="{col}")) |>
                   ungroup() |>
                   rename("impactyear" = "year") |>
                   mutate(year=i) |>
                   joyn::joyn(responsibilityshare_country,by=c("year","impactyear"),match_type="1:m",reportvar=FALSE,keep="inner",verbose=FALSE) |>
                   mutate(across(starts_with("npv"),~responsibilityshare*.,.names="{col}")) |>
                   group_by(code,year) |>
                   summarize(across(starts_with("npv"),~sum(.),.names="{col}")) |>
                   ungroup() |>
                   rename_with(~ stringr::str_replace(.x,pattern = "npv",replacement = "poor"))  |>
                   rename_at(vars(-code,-year), ~ paste0(., "_ea")) 

# Aggregated to the gas causing the damage
ea_gas_temp <- temp |>
               group_by(year) |>
               summarize(across(starts_with("npv"),~sum(.),.names="{col}")) |>
               ungroup() |>
               rename("impactyear" = "year") |>
               mutate(year=i) |>
               joyn::joyn(responsibilityshare_gas,by=c("year","impactyear"),match_type="1:m",reportvar=FALSE,keep="inner",verbose=FALSE) |>
               mutate(across(starts_with("npv"),~responsibilityshare*.,.names="{col}")) |>
               group_by(gas,year) |>
               summarize(across(starts_with("npv"),~sum(.),.names="{col}")) |>
               ungroup() |>
               rename_with(~ stringr::str_replace(.x,pattern = "npv",replacement = "poor"))  |>
               rename_at(vars(-gas,-year), ~ paste0(., "_ea")) 

# Aggregated to total damage in 2030
  damage_2030_temp <- temp |>
    filter(year==2030) |>
    summarize(added_poor = sum (npv_df0_300)) |>
    mutate(year=i)  
    
# Aggregated 2024 damage in future year
if (i==2024) {
damage_annual <- temp |>
                 select(-contains("df2"),-contains("df4")) |>
                 group_by(year) |>
                 summarize(across(starts_with("npv"),~sum(.),.names="{col}")) |>
                 ungroup() |>
                 rename("impactyear" = "year") |>
                 mutate(year=i) |>
                 joyn::joyn(responsibilityshare_gas,by=c("year","impactyear"),match_type="1:m",reportvar=FALSE,keep="inner",verbose=FALSE) |>
                 mutate(across(starts_with("npv"),~responsibilityshare*.,.names="{col}")) |>
                 rename_with(~ stringr::str_replace(.x,pattern = "npv",replacement = "poor"))  |>
                 rename_at(vars(-gas,-year), ~ paste0(., "_ea")) 
 save(damage_annual, file="03-outputdata/damage_annual.Rda")
}

# Aggregated to the country experiencing the damage
damage_country_temp <- temp |>
                       group_by(code) |>
                       summarize(poor_300_da_df0 = sum(npv_df0_300),
                                 poor_830_da_df0 = sum(npv_df0_830),
                                 poor_npl_da_df0 = sum(npv_df0_npl),
                                 poor_mea_da_df0 = sum(npv_df0_mea),
                                 poor_gap_da_df0 = sum(npv_df0_gap)) |>
                       ungroup() |>
                       mutate(year=i) |>
                       joyn::joyn(pop,by=c("year","code"),match_type="1:1",reportvar=FALSE,verbose=FALSE,keep="inner") |>
                       mutate(rate_300_da_df0 = poor_300_da_df0/pop*100,
                              rate_830_da_df0 = poor_830_da_df0/pop*100,
                              rate_npl_da_df0 = poor_npl_da_df0/pop*100,
                              rate_mea_da_df0 = poor_mea_da_df0/pop,
                              rate_gap_da_df0 = poor_gap_da_df0/pop) 

# Append with results from other emission years
if (i == 1990) {
  damage_2030 <- damage_2030_temp  
  ea_country     <- ea_country_temp  
  ea_gas         <- ea_gas_temp 
  damage_country <- damage_country_temp  
}
if (i != 1990) {
  ea_country     <- rbind(ea_country,ea_country_temp)
  ea_gas         <- rbind(ea_gas,ea_gas_temp)
  damage_2030 <- rbind(damage_2030,damage_2030_temp)  
  damage_country <- rbind(damage_country,damage_country_temp)  
}
}

 save(damage_country, file="03-outputdata/damage_country.Rda")
 save(damage_2030,    file="03-outputdata/damage_2030.Rda")
 save(ea_country,     file="03-outputdata/ea_country.Rda")
 save(ea_gas,         file="03-outputdata/ea_gas.Rda")