si.Rmd

As a numerical example, we adjust a total final energy of 53 GJ per capita from the LED scenario (Grubler et al. (2018) [@grubler_low_2018]), first by the household share of the total European energy footprint in 2015 (around 0.62, calculated in EXIOBASE), and then the share of total adult equivalents in the total European population in 2015 (also around `r ae_share_of_pop`, calculated using the EUROSTAT HBS, number of households per country, and population data per country). A total final energy of 53 GJ/capita is therefore adjusted to a household final energy of 53 GJ/adult equivalent in Europe ((53 total GJ/capita * 0.62 household share of total footprint)/0.62 adult equivalent share of total population = 53 household GJ/adult equivalent).

The decarbonisation scenario final energy numbers in 2050, presented in Table 1 of the main paper, were originally in total GJ per capita: 94 GJ/capita (SSP2-1.9), 87 GJ/capita (SSP1-1.9), 84 GJ/capita (IEA ETP B2DS), 64 GJ/capita (GEA-efficiency), 53 GJ/capita (LED), and 15.3 GJ/capita (DLE). Because of the similar relative shares of the household part of the total European energy footprint (~0.62), and the adult equivalent share of the total population in our sample (also around 0.62), these final energy numbers end up close to the same when adjusted to household per adult equivalent. The original numbers for the SSP and GEA-efficiency scenarios are from the International Institute for Applied Systems Analysis (IIASA) scenario database [@riahi_shared_2017 @gea_gea_nodate]. The SSP total GJ/capita numbers are for the OECD region, while the GEA-efficiency total GJ/capita number is for their 'West EU' region. The LED total GJ/capita number is from Grubler et al. (2018) [@grubler_low_2018], and the IEA ETP B2DS total GJ/capita number is from the Supplementary Table 27 in the supplementary information document of Grubler et al. (2018) [@grubler_low_2018]. The LED and IEA ETP B2DS total GJ/capita numbers are both for the Global North region. We also refer in the main paper to the LED numbers for the Global South (20 total GJ/capita) and the world (27 total GJ/capita). Finally, the DLE number is one number for the world, and while they give a range of 13-18.4 total GJ/capita, we take their average of 15.3 total GJ/capita [@millward-hopkins_providing_2020].

Our European expenditure deciles were constructed having the exact same number of adult equivalents per decile. When comparing with external per capita numbers, however, there are not the same number of population per decile because of differences in non-adult-equivalent-normalized people per household between income quintiles per country, and between countries. In Figure Sx we show an estimate of population per European expenditure decile. We use this to re-estimate our energy footprint per European expenditure decile in per capita terms, and then re-create Figure 5 from the main paper (Figure Sxx below) in per capita terms. 

Because we have the number of adult equivalents per country, and the total population per country, we could use both to calculate the total population per adult equivalent ratio for each country. We applied these ratio to the adult equivalents from different countries making up each European expenditure decile, as so we could estimate total population per European expenditure decile taking into account differences in household per capita between countries, but not between income quintiles within each country. Figure Sx shows these population estimates per European expenditure decile. 

```{r figureSx, fig.cap="**Figure S5:**"}

ggplot(a, aes(x=factor(eu_ntile_name), y=value)) +
  geom_col() +
  ylab("number (in millions)") +
  xlab("") +
  facet_wrap(~indicator, scales="free_y") +
  theme_ipsum() +
  theme(axis.text.x = element_text(angle = 90)) 

ggsave(here("analysis", "figures", "figureSx.pdf"), device=cairo_pdf)

```

## Figure 5 from manuscript in household final energy per capita

Here we re-create Figure 5 from the main paper after estimating our energy footprint per European expenditure decile in per capita terms (using the population per decile estimates from above), instead of per adult equivalent terms. Now we only adjust the decarbonisation scenario final energy numbers from total GJ per capita to household GJ per capita, so, for example, 53 total GJ/capita becomes: 53 * 0.62 = 33 household GJ/capita. 

```{r }

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(eu_q_rank,total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)) %>%
  left_join(a_tmp, by = "eu_q_rank") %>%
  mutate(pc_energy_use_tj = total_energy_use_tj_new/total_population,
         pc_energy_use_gj = pc_energy_use_tj*1000)

df_energy_deciles = pdat_final_demand %>%
  select(eu_q_rank, pc_energy_use_gj)

ineq_curr = df_energy_deciles$pc_energy_use_gj[10]/df_energy_deciles$pc_energy_use_gj[1]

```

```{r }

library(readxl)

df_scenario_info = read_excel(here("analysis/data/raw/scenarios_additions.xlsx"), sheet="overview") %>%
  select(scenario, fe_gj_pc = final_energy_gj_per_capita_2050,
         ccs_required = primary_energy_fossil_w_ccs2050_ej,
         description) %>%
  arrange(fe_gj_pc) %>%
  mutate(fe_gj_pc = fe_gj_pc*0.62,
         fe_gj_pc = round(fe_gj_pc),
         ccs_required = round(ccs_required)) 

mea = c(min(df_scenario_info$fe_gj_pc),max(df_scenario_info$fe_gj_pc))
mer = c(min(df_scenario_info$fe_gj_pc),53) #before was 33: multiplied 53 by 0.62 and rounded up

c_mean_mea = round((17+mea[2])/2) #multiplied 27 by 0.62 and rounded up
c_mean_mer = round((mer[1]+mer[2])/2)

```

```{r , eval = FALSE}

# run once to save file. If changing scenarios included or input data, need to re-run and save 'scenarios_fine.rds'

# vectorized function that returns scaled quantiles given 
#quantile index column and column with quantile averages
#qidx: quantile index
#qavg_to_scale: column to scale
#first_target: target value of first quantile
#mean_target: target mean of scaled quantiles
scaled_quantiles <- function(.data, 
                          qidx, 
                          qavg_to_scale, 
                          first_target, 
                          mean_target) {
  
  # cumbersomely extract current quantile mean
  mean_current = .data %>%
    ungroup() %>%
    summarise(mean_cur = first(mean({{qavg_to_scale}}))) %>%
    pull(mean_cur)
    
  # cumbersomely extract current first wuantile value
  first_current = .data %>%
    ungroup() %>%
    arrange({{qidx}}) %>%
    summarise(first_cur = first({{qavg_to_scale}})) %>%
    pull(first_cur)
  
  df_tmp = .data %>%
    mutate(tmp = {{qavg_to_scale}}*mean_target/mean_current)
  
  first_tmp = df_tmp$tmp[1]
  
  df_tmp = df_tmp %>%
    mutate(scaled = mean_target-(mean_target-tmp) * (mean_target-first_target)/(mean_target-first_tmp)) %>%
    select({{qidx}}, scaled) %>%
    mutate(v_mean = mean_target, v_first = first_target)
}

## run once to save file
df_all = NULL
for (min_energy in seq(from=mer[1], to=mer[2], by=0.25)) {
  for (mean_energy in seq(from=mea[1], to=mea[2], by=0.25)) {
    if (min_energy <= mean_energy) {
      df_all = df_all %>%
      bind_rows(df_energy_deciles %>%
                  scaled_quantiles(eu_q_rank, pc_energy_use_gj, min_energy, mean_energy))
    }
  }
}

saveRDS(df_all, here("analysis/data/derived/si/scenarios_fine.rds"))
```

```{r , fig.width=7, fig.height=5.5, fig.cap="**Figure S6:**"}

round_by = 10

df_all = readRDS(here("analysis/data/derived/si/scenarios_fine.rds")) %>%
  filter(eu_q_rank %in% c(1,10)) %>%
  group_by(v_mean, v_first) %>%
  summarise(ratio = last(scaled)/first(scaled)) %>%
  mutate(bin_ratio = if_else((ratio*100)%%round_by > round_by*0.5, 
                             ratio*100+(round_by-(ratio*100)%%round_by),
                             ratio*100-(ratio*100)%%round_by)) %>%
  group_by(bin_ratio, v_first) %>%
  summarise(v_mean = mean(v_mean)) %>%
  mutate(bin_ratio = bin_ratio*0.01)

df_scenario = df_all %>%
  filter(v_mean %in% df_scenario_info$fe_gj_pc) %>%
  left_join(df_scenario_info, by=c("v_mean"="fe_gj_pc")) %>%
  filter(!(scenario == "DLE")) %>%
  mutate(scenario = dplyr::recode(scenario,
              "LED" = "LED (33)",
              "IEA ETP B2DS" = "IEA ETP B2DS (52)",
              "GEA efficiency" = "GEA efficiency (40)",
              "SSP1-1.9" = "SSP1-1.9 (54)",
              "SSP2-1.9" = "SSP2-1.9 (58)",
              "SSP4-1.9" = "SSP4-1.9 (64)",
              "SSP3-3.4" = "SSP3-3.4 (90)",
              "SSP5-6.0" = "SSP5-6.0 (94)"))

library(wesanderson)

a = df_all %>%
  ggplot(aes(x=v_first, y=bin_ratio, fill=v_mean)) +
  #geom_tile() +
  geom_hline(yintercept = ineq_curr, alpha=0.8, color="grey", linetype=2) +
  geom_line(data=df_scenario, aes(color=scenario, group=scenario)) +
  annotate(geom="text", x=max(df_all$v_first)-7.5,y=ineq_curr+0.6,label = "Current (2015) 10:10 ratio") +
  #scale_fill_viridis("Mean energy\navailable") +
  #scale_fill_gradient("Mean energy\navailable (GJ/cap)",
  #                    low=wes_palette("Chevalier1")[3], 
  #                    high = wes_palette("Rushmore1")[4]) +
  #scale_color_manual(values=wes_palette("Darjeeling1")) +
  theme_minimal() +
  labs(x="Minimum energy requirement (GJ/cap)", y="Maximum energy inequality (10:10 ratio)", color = "Scenario")+
  theme(text=element_text(family="Liberation Sans Narrow"),
        axis.text.x = element_text(size = 13),
        axis.text.y = element_text(size = 13)) +
  #scale_color_discrete(name = "Scenario (mean GJ/cap)") +
  xlim(10,53) +
  #ylim(0,15) +
  scale_y_continuous(breaks = c(2.5,5,7.5,10,12.5,15), limits = c(1,15)) +
  scale_color_discrete(name = "Scenario\n(mean GJ/cap)", breaks=c("LED (33)",
                               "GEA efficiency (40)",
                               "IEA ETP B2DS (52)",
                               "SSP1-1.9 (54)",
                               "SSP2-1.9 (58)",
                               "SSP4-1.9 (64)",
                               "SSP3-3.4 (90)",
                               "SSP5-6.0 (94)")) #+
  #theme_ipsum()
a
ggsave(here("analysis", "figures", "figureSxx.pdf"))
```

\newpage

# References 

<!-- The following line ensures the references appear here for the MS Word or HTML output files, rather than right at the end of the document (this will not work for PDF files):  -->

<div id="refs"></div>