
pacman::p_load(pacman, rio, raster, rgdal, stars, 
               ggplot2, maps,  sf, RColorBrewer,
               #rasterVis, 
               ggthemes, viridis, tidyverse,
               RColorBrewer,
               extrafont, dplyr, mapdeck, echarts4r, DT,
               data.table,readr, plotly, htmlwidgets, mapview,
               leafpop, mapedit, reshape2, rasterVis, tidync, pastecs,
               foreign, haven, Rfast, forcats, gdata,
               biscale, cowplot, gtools, PerformanceAnalytics,
               fastDummies, ipumsr, stringr, rgeos, maptools, tmap,
               scales, sjlabelled, jtools, huxtable, texreg,
               robustbase, readxl, miceadds , estimatr) 



memory.limit(10000000)

########################################################
########################################################

# CSV file is from HH survey excel sheet of CPAT excel file  

HH <- read_csv("data/HH-EECs-V3.csv")

hh <- HH[HH$type == "Deciles" & HH$stat_type == "mean" , ]

hh <- hh %>% 
  mutate(biomass_share = ccl_share  + fwd_share) %>% 
  mutate(fossil_elec_share =  coa_share + die_share + ethanol_share + gso_share + ker_share + lpg_share + nga_share + oil_share + ely_share) %>% 
  mutate(fossil_share =  coa_share + die_share + ethanol_share + gso_share + ker_share + lpg_share + nga_share + oil_share)  


# biomass includes charcoal, crop waste, dung, and wood based on the presentation
# in CPAT ethanol is considered as biomass see the distribution tab of the excel sheet

###

# country classification by income 

hh <- hh %>% 
  mutate(income_class = ifelse(iso3 == "MDG" | iso3 == "MLI" | iso3 == "RWA" , "LICs" , 
                               ifelse( iso3 == "BGD" | iso3 == "BOL" | iso3 == "CIV" | iso3 == "EGY" | iso3 == "GHA" | iso3 == "HND" |
                                         iso3 == "IDN" | iso3 == "IND" | iso3 == "LKA" | iso3 == "MAR" | iso3 == "NPL" | iso3 == "PAK" | 
                                         iso3 == "PHL" | iso3 == "UKR" | iso3 == "VNM"  , "LMICs" ,
                                       ifelse( iso3 == "ARG" | iso3 == "BGR" | iso3 == "BRA" | iso3 == "CHN" | iso3 == "COL" | iso3 == "CRI" |
                                                 iso3 == "DOM" | iso3 == "ECU" | iso3 == "KAZ" | iso3 == "MEX"  | iso3 == "MKD"  | iso3 == "MYS"  | 
                                                 iso3 == "PER" | iso3 == "SRB" | iso3 == "THA" | iso3 == "TUR"  , "UMICs" ,
                                               ifelse( iso3 == "AUT" | iso3 == "BEL" | iso3 == "CAN" | iso3 == "CHL" | iso3 == "CYP" | iso3 == "CZE" |
                                                         iso3 == "DEU" | iso3 == "DNK" | iso3 == "ESP" | iso3 == "EST" | iso3 == "FIN" | iso3 == "FRA" |
                                                         iso3 == "GBR" | iso3 == "GRC" | iso3 == "HRV" | iso3 == "HUN" | iso3 == "IRL" | iso3 == "ITA" |
                                                         iso3 == "LTU" | iso3 == "LUX" | iso3 == "LVA" | iso3 == "MLT" | iso3 == "NLD" | iso3 == "POL" |
                                                         iso3 == "PRT" | iso3 == "ROU" | iso3 == "SVK" | iso3 == "SVN" | iso3 == "SWE" | iso3 == "URY" |
                                                         
                                                         iso3 == "USA", "HICs", NA)))))




# Budget shares

hh <- hh %>% 
  rename(npg_share = nga_share) %>% 
  group_by(income_class , quant_cons) %>%  mutate(ccl_cs = mean(ccl_share)) %>%  mutate(coa_cs = mean(coa_share)) %>%  mutate(die_cs = mean(die_share)) %>% 
  mutate(ethanol_cs = mean(ethanol_share)) %>%  mutate(ely_cs = mean(ely_share)) %>% mutate(fwd_cs = mean(fwd_share)) %>%  mutate(gso_cs = mean(gso_share)) %>% 
  mutate(ker_cs = mean(ker_share)) %>% mutate(lpg_cs = mean(lpg_share)) %>%  mutate(npg_cs = mean(npg_share)) %>% mutate(oil_cs = mean(oil_share)) %>% 
  mutate(biomass_cs = mean(biomass_share, na.rm = TRUE)) %>% mutate(fossil_elec_cs = mean(fossil_elec_share, na.rm = TRUE)) %>% mutate(fossil_cs = mean(fossil_share, na.rm = TRUE))



hh_budget_share <- hh %>% 
  group_by(income_class , quant_cons) %>%
  summarise(across(c(ccl_share, coa_share, die_share, ethanol_share , ely_share ,fwd_share , gso_share ,  ker_share , lpg_share , npg_share , oil_share ), list(mean = mean)))



hh_budget_share <- hh_budget_share %>% 
  mutate(High_pollutants = (ccl_share_mean + fwd_share_mean +coa_share_mean + ker_share_mean)) %>% 
  mutate(Medium_pollutants = (die_share_mean + ethanol_share_mean + gso_share_mean + oil_share_mean)) %>% 
  mutate(Low_pollutants = (ely_share_mean + lpg_share_mean + npg_share_mean ))

hh_share_budget <- melt(setDT(hh_budget_share), id.vars = c("income_class" , "quant_cons" ), variable.name = "fuels")


# Budgets

hh <- hh %>% 
  rename(NPG = "Natural gas")

fuels <- c("Charcoal"  ,  "Coal"  , "Diesel"  , "Ethanol"  ,  "Electricity" , "Firewood"  ,  "Gasoline"   ,  "Kerosene" , "LPG" ,    "NPG"   ,  "Oil"  )

##########################################################

# Graphs

########

# Budget shares

# High_pollutants

ggplot(hh_share_budget[hh_share_budget$fuels == "High_pollutants"  , ],                         # Draw barplot with grouping & stacking
       aes(x = as.factor(quant_cons),
           y = value,
           fill = fuels)) + 
  geom_bar(stat = "identity",
           position = "dodge") +
  facet_grid(~ factor(income_class, levels=c('LICs', 'LMICs', 'UMICs', 'HICs'))) +
  labs(x = "Income decile", y = "Budget shares", group = "Country Groups") +
  scale_fill_manual("Fuels", values = c("High_pollutants" = "red3", "Medium_pollutants" = "blue3", "Low_pollutants" = "green3"))


# Medium_pollutants


ggplot(hh_share_budget[hh_share_budget$fuels == "Medium_pollutants"  , ],                         # Draw barplot with grouping & stacking
       aes(x = as.factor(quant_cons),
           y = value,
           fill = fuels)) + 
  geom_bar(stat = "identity",
           position = "dodge") +
  facet_grid(~ factor(income_class, levels=c('LICs', 'LMICs', 'UMICs', 'HICs'))) +
  labs(x = "Income decile", y = "Budget shares", group = "Country Groups") +
  scale_fill_manual("Fuels", values = c("High_pollutants" = "red3", "Medium_pollutants" = "blue3", "Low_pollutants" = "green3"))


# Low_pollutants

ggplot(hh_share_budget[hh_share_budget$fuels == "Low_pollutants"  , ],                         # Draw barplot with grouping & stacking
       aes(x = as.factor(quant_cons),
           y = value,
           fill = fuels)) + 
  geom_bar(stat = "identity",
           position = "dodge") +
  facet_grid(~ factor(income_class, levels=c('LICs', 'LMICs', 'UMICs', 'HICs'))) +
  labs(x = "Income decile", y = "Budget shares", group = "Country Groups") +
  scale_fill_manual("Fuels", values = c("High_pollutants" = "red3", "Medium_pollutants" = "blue3", "Low_pollutants" = "green3"))






