###################
### INTODUCTION ###
###################
# This script prepares the welfare data
# It requires the following packages: dplyr, joyn, pipr
library(dplyr)
rm(list=ls())

#######################
### HISTORICAL DATA ###
#######################
load("02-intermediatedata/class.Rda")

# Distributions loaded below come from this paper: https://doi.org/10.1596/43363
load("01-inputdata/welfare/comparable_pref.Rda")

# Create welfare type
lineup <- pipr::get_stats(fill_gaps = TRUE,nowcast=TRUE) |>
          filter(year==2025) |>
          select(country_code,welfare_type) |>
          rename("code"="country_code") |>
          distinct() |>
          mutate(welfare_type = if_else(code=="PHL","income",welfare_type),
                 welfare_type = if_else(code=="XKX","consumption",welfare_type),
                 welfare_type = if_else(code=="ALB","income",welfare_type))

# Convert consumption vectors to income
inccon <- expand.grid(year=1990:2025,code=unique(class$code)) |>
          joyn::joyn(lineup,by="code",match_type="m:1", reportvar=FALSE)
welfare <- comparable_pref  |>
           # Convert to income
           mutate(welf=((welf-0.42)/0.96)^(1/0.94)) |>
           # Bottom-code welfare at 1 cent
           mutate(welf=replace(welf, welf <0.28 | is.na(welf),0.28)) |>
           # Keep relevant years
           filter(year>=1990) |>
           joyn::joyn(inccon,by=c("code","year"),match_type="m:1",reportvar=FALSE) 
rm(lineup,comparable_pref)

###################
### PROJECTIONS ###
###################
load("02-intermediatedata/growth.Rda")
growth <- growth |> filter(scenario=="med") |> select(-scenario)
inccon <- inccon |> filter(year==2025) |> select(code,welfare_type)

welfare_growth <- joyn::joyn(growth,inccon,by="code",match_type="m:1",reportvar=FALSE) |>
                  group_by(code) |>
                  arrange(year) |>
                  mutate(# Calculate growth
                         welfmultiplier = 1+growth,
                         welfmultiplier = if_else(year==2025,1,welfmultiplier),
                         # Calculate cumulative growth
                         cumwelfmultiplier = lag(cumprod(welfmultiplier),default=1),
                         cumwelfmultiplier = cumwelfmultiplier*welfmultiplier[row_number()==2],
                         cumwelfmultiplier = if_else(year==2025,1,cumwelfmultiplier)) |>
                  ungroup() |>
                  select(code,year,cumwelfmultiplier)
rm(growth,inccon)

# Fill 2025 with actual welfare vector
welfare2025 <- welfare |> 
               filter(year==2025) |> 
               select(code,welf) |>
               group_by(code) |>
               mutate(obs =row_number()) |>
               rename(welf2025 = welf)

# Prepare empty file to store results
welfare_future <- expand.grid(code=unique(welfare$code),year=2025:2100,obs=1:1000) |>
                  joyn::joyn(welfare2025,by=c("code","obs"),match_type="m:1",reportvar=FALSE) |>   
                  joyn::joyn(welfare_growth,by=c("code","year"),match_type="m:1",reportvar=FALSE) |>
                  select(-obs) |>
                  # Project welfare
                  mutate(welf = if_else(year!=2025,welf2025*cumwelfmultiplier,welf2025),
                         welf = replace(welf, welf <0.01,0.01)) |>
                  select(code,year,welf)

# Append historical data
welfare <- welfare |> 
           filter(year!=2025) |> 
           select(-welfare_type) |>
           bind_rows(welfare_future) 

rm(welfare_growth,welfare2025,welfare_future)

save(welfare,file="02-intermediatedata/welfare.Rda")

#############################
### NATIONAL POVERTY LINE ###
#############################
# Construct national poverty lines in order to derive povety rates at these liens
load("02-intermediatedata/welfare.Rda")

npl <- welfare |>
  filter(year==2025) |>
  group_by(code) |>
  summarize(npl = max(1.30+median(welf)/2,3)) |>
  ungroup()
save(npl,file="02-intermediatedata/npl.Rda")

#############################################
### COUNTRY-YEAR-LEVEL POVERTY STATISTICS ###
#############################################
load("02-Intermediatedata/pop.Rda")

# Calculate country-year level poverty and prosperity gap
welfare_country <- welfare |> 
                   joyn::joyn(npl,match_type="m:1",by="code",reportvar=FALSE) |>
                   group_by(code,year) |>
                   summarize(rate_300_ac = mean(welf<3.00)*100,
                             rate_830_ac = mean(welf<8.30)*100,
                             rate_npl_ac = mean(welf<npl)*100,
                             rate_mea_ac = mean(welf),
                             rate_gap_ac = mean(if_else(28/welf<100,28/welf,100))) |>
                   ungroup() |>
                   joyn::joyn(pop,by=c("year","code"),match_type="1:1",reportvar=FALSE, keep="left") |>
                   mutate(poor_300_ac = rate_300_ac/100*pop,
                          poor_830_ac = rate_830_ac/100*pop,
                          poor_npl_ac = rate_npl_ac/100*pop,
                          poor_mea_ac = rate_mea_ac*pop,
                          poor_gap_ac = rate_gap_ac*pop)

save(welfare_country,file="02-intermediatedata/welfare_country.Rda")

##########################################
### REGIONAL/GLOBAL POVERTY STATISTICS ###
##########################################
welfare_region <- welfare_country |>
                  joyn::joyn(class,by="code",match_type="m:1",reportvar=FALSE) |>
                  group_by(year,region) |> 
                  summarize(rate_300_ac = weighted.mean(rate_300_ac,pop),
                            rate_830_ac = weighted.mean(rate_830_ac,pop),
                            rate_npl_ac = weighted.mean(rate_npl_ac,pop),
                            rate_mea_ac = weighted.mean(rate_mea_ac,pop),
                            rate_gap_ac = weighted.mean(rate_gap_ac,pop),
                            pop         = sum(pop),
                            poor_300_ac = sum(poor_300_ac),
                            poor_830_ac = sum(poor_830_ac),
                            poor_npl_ac = sum(poor_npl_ac),
                            poor_mea_ac = sum(poor_mea_ac),
                            poor_gap_ac = sum(poor_gap_ac))

welfare_global <- welfare_country |>
                  joyn::joyn(class,by="code",match_type="m:1",reportvar=FALSE) |>
                  group_by(year) |> 
                  summarize(rate_300_ac = weighted.mean(rate_300_ac,pop),
                            rate_830_ac = weighted.mean(rate_830_ac,pop),
                            rate_npl_ac = weighted.mean(rate_npl_ac,pop),
                            rate_mea_ac = weighted.mean(rate_mea_ac,pop),
                            rate_gap_ac = weighted.mean(rate_gap_ac,pop),
                            pop         = sum(pop),
                            poor_300_ac = sum(poor_300_ac),
                            poor_830_ac = sum(poor_830_ac),
                            poor_npl_ac = sum(poor_npl_ac),
                            poor_mea_ac = sum(poor_mea_ac),
                            poor_gap_ac = sum(poor_gap_ac)) |>
                  mutate(region="World") |>
                  rbind(welfare_region)

save(welfare_global,file="02-intermediatedata/welfare_global.Rda")