si.Rmd

title: "The energy and carbon inequality corridor for a 1.5°C compatible and just Europe: supplementary information"
author:
  - Ingram S. Jaccard:
      email: jaccard@pik-potsdam.de
      institute: [PIK]
      correspondence: false
  - Peter-Paul Pichler:
      email: pichler@pik-potsdam.de
      institute: [PIK]
      correspondence: false
  - Johannes Többen:
      email: toebben@pik-potsdam.de
      institute: [PIK, GWS]
      correspondence: false
  - Helga Weisz:
      email: weisz@pik-potsdam.de
      institute: [PIK, HU]
      correspondence: false
institute:
  - PIK: Social Metabolism and Impacts, Potsdam Institute for Climate Impact Research, Member of the Leibniz Association, PO Box 60 12 03, Potsdam, 14412, Germany
  - HU: Department of Cultural History & Theory and Department of Social Sciences, Humboldt-University Berlin, Unter den Linden 6, Berlin, 10117, Germany
  - GWS: Gesellschaft für Wirtschaftliche Strukturforschung (GWS) mbH, Heinrichstraße 30, 49080 Osnabrück, Germany
output:
  word_document:
    fig_caption: yes
    reference_docx: ../templates/template.docx
    pandoc_args:
    - --lua-filter=../templates/scholarly-metadata.lua
    - --lua-filter=../templates/author-info-blocks.lua
    - --lua-filter=../templates/pagebreak.lua
bibliography: references.bib
csl: ../templates/vancouver.csl
content: |
  35 pages, 13 tables, 7 figures
always_allow_html: yes

Content: r rmarkdown::metadata$content

knitr::opts_chunk$set(echo = TRUE)
devtools::source_gist("c83e078bf8c81b035e32c3fc0cf04ee8", 
                      filename = 'render_toc.R')

render_toc("si.Rmd")

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = FALSE,
  comment = "#>",
  fig.path = "../figures/",
  dpi = 300
)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,
               janitor,
               here,
               wbstats,
               ISOcodes,
               viridis,
               hrbrthemes,
               wesanderson,
               glue,
               ggridges,
               patchwork,
               kableExtra,
               readxl,
               flextable)

pal <- wes_palette("Cavalcanti1", 5, type = "discrete")
extrafont::loadfonts()

options(scipen=999)

library(here)

Supplementary materials and methods

In this paper we first stratified EXIOBASE national household energy and carbon footprints by income quintile, using data from the EUROSTAT Household Budget Survey (HBS). We then ranked all of these national income quintiles (140 quintiles in total: 28 countries x 5 national income quintiles) according to their mean expenditure in purchasing power standards (PPS), and aggregated them into 10 European expenditure deciles. The resultant household energy and carbon footprints per European expenditure decile in the year 2015 are then our level of analysis through the whole main paper. In this SI document, we explain these steps in more detail, discuss limitations, and provide some supplementary results at the end. All data, relevant tables, and R code are available at: https://gitlab.pik-potsdam.de/pichler/europe.inequality.

Income-stratified national household energy and carbon footprints

The first step was stratifying an Environmentally-Extended Multi-Regional Input-Output (EE-MRIO) model's household final demand expenditure by income quantile, and then multiplying this income-stratified expenditure by 'total' (direct and indirect supply chain) environmental intensities calculated in the EE-MRIO to estimate income-stratified national household environmental footprints.

An EE-MRIO can calculate national environmental footprints from final demand expenditure. Final demand expenditure on different production sectors, in monetary units, is multiplied by 'direct and indirect supply chain' environmental intensities per production sector to calculate the environmental footprint (for one year). The 'direct and indirect supply chain' intensities per production sector estimate the total physical amount of environmental extension, whether emissions, energy use, etc., anywhere along the supply chain per monetary unit of final demand expenditure.

In this way, the global amount of environmental pressure in that year is allocated to each country based on their final demand expenditure (this 'footprint' approach is also called consumption-based accounting as opposed to territorial-based accounting). The 'direct and indirect supply chain' intensities per production sector are calculated using standard EE-MRIO calculations [@miller_input-output_1985].

The final demand expenditure, and thus the national footprint, is typically disaggregated into several final demand categories, including households, non-profits serving households, government, gross-fixed capital formation, change in inventories and valuables. Current publicly available EE-MRIOs do not have final demand expenditure stratified by income quantile. The underlying distribution of final demand expenditure and footprint across income groups is therefore hidden.

In order to stratify the EE-MRIO final demand expenditure and footprint by income quantile, a second data input is needed with information on the income/expenditure distribution. This could be a household budget survey (HBS) stratified by income quantile or, at a more aggregate level, the national income distribution. In this paper we use:

1. EE-MRIO: EXIOBASE version3.7 industry-by-industry, from: https://zenodo.org/record/3583071#.XjC7kSN4wpY [accessed on 12.03.2020] [@stadler_exiobase_2018; @stadler_exiobase_2019]

2. HBS: European household budget survey from EUROSTAT, macro-data, from : https://ec.europa.eu/eurostat/web/household-budget-surveys/database [accessed on 22.05.2020] [@eurostat_database_2021]

We discuss each of these in turn, including additional data inputs needed to complement the HBS. We then present our methodology for our results from the main paper and some limitations. The final sub-section in 'supplementary materials and methods' presents an alternative methodology, and the final section of this supplementary information document presents supplementary results from the methodology used in the main paper and the alternative methodology. The whole analysis was performed in RStudio [@r_core_team_r:_2020], using a variety of R packages [@wickham_welcome_2019; @firke_janitor:_2021; @muller_here:_2020; @piburn_wbstats:_2020; @buchta_isocodes:_2020; @garnier_viridis:_2018; @rudis_hrbrthemes:_2020; @ram_wesanderson:_2018; @hester_glue:_2020; @wilke_ggridges:_2021; @pedersen_patchwork:_2020; @zhu_kableextra:_2020; @wickham_readxl:_2019; @gohel_flextable:_2021; @south_rworldmap:_2016; @arnold_ggthemes:_2021].

EXIOBASE

We use standard input-output calculations to calculate 'direct and indirect supply chain' intensity vectors in EXIOBASE [@miller_input-output_1985]. EXIOBASE publishes the A matrix, the final demand matrix, the satellite extensions matrix, and satellite extensions direct from final demand matrix. We use the industry-by-industry (ixi) EXIOBASE data tables from EXIOBASE version3.7. This means 163 industry production sectors and 6 final demand categories for 49 regions worldwide (44 countries and 5 rest-of-world regions), from 1995 - 2016. All monetary units are in million current Euros. Stadler et al. (2018) [@stadler_exiobase_2018] describe the EXIOBASE version3.7 compilation procedure in detail, including nine supporting information documents with further detailed information on the compilation of the monetary tables (S1), energy (S2), emissions (S3), and others.

For each year, we first load the A matrix and calculate the Leontief inverse (the inverse of the A matrix). We load the final demand matrix and calculate total output (x) by pre-multiplying the Leontief inverse (L) by the row sums of the final demand matrix (Y):

x = L ° sum(Y)

We then load the satellite extensions matrix and extract the relevant extensions. To calculate direct intensity vectors (DIV) we divide the satellite extension vectors (f) by total output (x):

DIV = f/x

The 'direct and indirect supply chain' intensity vectors (TIV) are calculated by pre-multiplying the direct intensity vectors (DIV) by the Leontief inverse (L):

TIV = DIV ° L

The footprint is then calculated by row-wise multiplying the TIV by final demand:

fp = TIV * Y

Before that final step, however, we stratify national household final demand expenditure by income quintile according to the structure of the EUROSTAT HBS, explained in the 'EUROSTAT HBS' section below.

The results in the main paper also present the footprint broken down by its domestic, other European, and non-European parts. To calculate these domestic and foreign parts of the footprint, we row-wise multiply the direct intensity vectors (DIV) by the Leontief inverse (L):

TIV breakdown = DIV * L

Environmental extensions

The environmental extensions we use are emissions of CO2-equivalence (in kilograms) and 'energy carrier net: total' (in terajoules). Our final energy use extension is created by subtracting 'energy carrier net: loss' from 'energy carrier net: total'. We create the CO2-equivalence extension by summing together the greenhouse gases CO2, CH4, N2O, SF6, HFCs, and PFCs, from combustion, noncombustion, agriculture and waste. We use Global Warming Potential (GWP) values for a 100-year time horizon taken from the IPCC Fifth Assessment Report [@myhre_g._anthropogenic_2013 (p.73-79)]: 28 for CH4, 265 for N2O and 23500 for SF6 (HFCs and PFCs are in CO2-equivalence already in the EXIOBASE environmental extensions).

The creation of the net energy use extensions (final energy use extensions and losses) in EXIOBASE is explained in the supplementary information documents of [@usubiagaliano_energy_2020], [@mastrucci_framework_2020] and [@vita_durable_2021]. The extended energy balances of the International Energy Agency (IEA) are first split into gross energy supply and use tables, and converted from the territorial to residence principle following the System of Environmental Economic Accounting (SEEA). This first step was initially explained by Stadler et al. (2018) [@stadler_exiobase_2018] in their Supporting Information 2, where they describe the compilation of the EXIOBASE version 3 original energy extensions: primary energy supply, gross energy supply, gross energy use, and emission-relevant energy use.

The net energy use accounts (final energy use accounts and losses) were developed by [@usubiagaliano_energy_2020] and [@mastrucci_framework_2020] for later EXIOBASE 3 versions. They refer to the final use of energy products plus all losses of energy, and minus exports of energy products [@mastrucci_framework_2020]. Net energy use is isolated from the gross energy use tables, and allocated to the EXIOBASE industries, products and final demand categories as in the original energy extensions [@usubiagaliano_energy_2020; @mastrucci_framework_2020]. Total net energy is disaggregated into final use categories with a separate losses extension. The energy losses extension subtracted from total net energy gives final energy use.

The conversion to the residence principle means that the EXIOBASE energy extensions refer to the functional border of a country's economy. In this case, the system border is defined by the 'residence' of the agent. This means that energy supply and use from international transport by ships, airplanes, fishing vessels, cars and trucks is allocated to the resident units of a country, independent from where these activities take place. In EXIOBASE version3.7, because emissions from these transport activities are estimated from the energy extensions via emission factors, the emissions extensions follow the residence principle as well.

Environmental extensions direct from households

For CO2-equivalence emissions and energy use direct from households, the EXIOBASE extensions provide one number of physical emissions and energy use direct from households per country. This number is made up of some mix between, primarily, direct household vehicle fuel use and direct household fuel use for housing heating and cooling, along with some other, smaller miscellaneous uses. To estimate the split between these three activities, we use two further EUROSTAT emissions and energy tables to split 'Total activities by households' between heating/cooling (HH_HEAT), transport activities (HH_TRA), and other (HH_OTH). Full definitions of what is included in these can be found in EUROSTAT's 'manual for air emissions accounts' (2015, p.66) [@eurostat_manual_2015]. While this disaggregation exists for nearly all EUROSTAT HBS countries, 2015 is the earliest year in our sample with complete coverage (except for Turkey in energy, which we impute using the Bulgarian splits). Therefore, we also use the 2015 splits between these 3 categories for our 2005 and 2010 estimates (the 2005 and 2010 results are shown only in this SI document). The two data tables are:

1. For energy, the EUROSTAT data table 'Energy supply and use by NACE Rev. 2 activity' [env_ac_pefasu] at: http://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=env_ac_pefasu [accessed on 03.06.2020]. [@eurostat_eurostat_2021]

2. For emissions to air, we download the EUROSTAT data table 'Air emissions accounts by NACE Rev. 2 activity' [env_ac_ainah_r2] at: https://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=env_ac_ainah_r2&lang=en [accessed on 03.06.2020]. [@eurostat_eurostat_2021-1]

EUROSTAT HBS

The EUROSTAT HBS compiles national household budget surveys from European countries. Countries provide EUROSTAT with their HBSs in micro-data form, and EUROSTAT has aggregated and published these survey data in macro-data form every five years since 1988. Of the publicly available data tables, we use two:

1. 'Mean consumption expenditure by income quintile (hbs_exp_t133)'

2. 'Structure of consumption expenditure by income quintile and COICOP consumption purpose (hbs_str_t223)'

'Mean consumption expenditure by income quintile' is expressed in purchasing power standards (PPS) per year, country and income quintile, in two possible units; per household and per adult equivalent. Purchasing power standard (PPS) is a common currency measure that eliminates differences in national price-levels. The income quintiles are determined by household, ie. all households in the sample are ranked by income and evenly distributed into five quintiles, thus there are the same number of households in each income quintile. For per adult equivalent, the mean consumption expenditure by income quintile is adjusted for household size between countries, years, and income quintiles, using the modified OECD scale: the first adult in the household is given a weight of 1.0, each adult thereafter 0.5, and each child 0.3 [@eurostat_description_2016].

The 'Structure of consumption expenditure by income quintile and COICOP consumption purpose' table gives the distribution of consumption expenditure across COICOP consumption categories in each income quintile, expressed in 'parts per mille (pm)'. There are three 'levels' of COICOP breakdown. All countries have level 1 (12 consumption categories) and 2, but there are only a few with level 3. For the current analysis, we select our own mix of COICOP level 1 and 2. We use COICOP level 1, except for the categories of food (CP01), housing (CP04) and mobility (CP07), where we use level 2, due to their importance for energy use and emissions.

This is primarily because of the diversity of level 2 categories, especially within those three aggregate level 1 categories. This is most clearly seen in the housing category, where, at level 2, housing is broken down into: 'Actual rentals for housing' (CP041), 'Imputed rentals for housing' (CP042), 'Maintenance and repair of the dwelling' (CP043), 'Water supply and miscellaneous services relating to the dwelling' (CP044), and 'Electricity, gas and other fuels' (CP045). Corresponding all housing-related EXIOBASE production sectors only to the HBS level 1 'housing' consumption category (CP04) would obscure the difference between level 2 categories with extremely different effects on the footprint (electricity consumption vs. rental payments). The COICOP consumption categories we thus use are: CP011 ('Food'), CP012 ('Non-alcoholic beverages'), CP02 ('Alcoholic beverages, tobacco and narcotics'), CP03 ('Clothing and footwear'), 'rent' (we collapse CP041 with CP042), CP043, CP044, CP045, CP05 ('Furnishings, household equipment and routine household maintenance'), CP06 ('Health'), CP071 ('Purchase of vehicles'), CP072 ('Operation of personal transport equipment'), CP073 ('Transport services'), CP08 ('Communications'), CP09 ('Recreation and culture'), CP10 ('Education'), CP11 ('Restaurants and hotels'), CP12 ('Miscellaneous goods and services').

As mentioned above, the EUROSTAT HBS 'mean consumption expenditure by income quintile' macro-data table is expressed per household and per adult equivalent. We use households per adult equivalent as our unit of analysis throughout the main paper, and in order to express our footprint estimates per adult equivalent, we needed the total number of private households in each country and year, and for this, two additional data tables were needed. We used:

1. The 'lfst_hhnhtych' table from EUROSTAT (Number of private households), selecting all available years and total households [accessed on 04.05.2020]. [@eurostat_eurostat_2021-2]

2. Norway is missing from the 'lfst_hhnhtych' EUROSTAT table. We download Norwegian data from the Norwegian statistical office: https://www.ssb.no/en/statbank/table/10986/. We select 'Private households', 'the whole country', no household type selection, all years (2005-2019), and continue with 'Table - Layout 1', then save the table as a 'Tab delimited without heading (csv)' file [accessed on 04.05.2020]. [@statistics_norway_10986:_2021]

Method for main paper results

For our results in the main paper, we first stratify the EXIOBASE national household final demand expenditure (before calculating the footprint) by income quintile, using the relative expenditure shares of each income quintile on the COICOP consumption categories in the EUROSTAT HBS. We do this for 2005, 2010 and 2015, but only present 2015 in the main paper.

The initial step was creating a correspondence table, assigning each production sector in EXIOBASE to one COICOP consumption category in the HBS. Because the share of each consumption category/production sector in the total amount of expenditure is not identical between the HBS and EXIOBASE, there are two possible methods for then stratifying the EXIOBASE household final demand expenditure: one that keeps the EXIOBASE production sector shares of total EXIOBASE expenditure intact (breaking the HBS consumption category shares), and one that keeps the HBS consumption category shares of total HBS expenditure intact (breaking the EXIOBASE production sector shares).

Our results in the main paper use the first method, keeping the EXIOBASE production sector shares of total EXIOBASE expenditure intact. This means that the total footprint is identical to when it is calculated in EXIOBASE without any stratification by income quintile. The alternative method, on the other hand, results in a different total footprint because a different amount of final demand expenditure in each production sector is now multiplied by the same original 'direct and indirect supply chain' intensity for that production sector. We present the alternative method and some results in the last sections of this document.

In this current section we use a two sector toy model to illustrate the method behind the main paper results. In the two sector model, the term 'sector' represents EXIOBASE production sectors that correspond to one HBS consumption category.

Step 1) Multiply 'mean consumption expenditure by income quintile' in purchasing power standards per household (pps hh) by the 'structure of consumption expenditure by income quintile and COICOP consumption purpose' (in 'parts per mille' or pm) to calculate the consumption expenditure structure in 'pps hh'. Then calculate the shares of each income quintile within each sector. Table S1 shows the two sector example, where 'pps hh' is multiplied by the shares of each sector ('parts per mille' divided by 1000) to calculate the expenditure on each sector per income quintile in 'pps hh' ('s1 (pps hh)' and 's2 (pps hh)'). 's1 (q share)' and 's2 (q share)' are each income quintile's share of the total amount of expenditure (in 'pps hh') on that sector, according to the HBS.

#"HBS structure with calculations of quintile shares per sector."

quintile = c("q1","q2","q3","q4","q5")
pps_hh = c(10,13,20,33,40)

pm_sector1 = c(200,300,400,500,600)
pm_sector2 = c(800,700,600,500,400)

hbs = data.frame(quintile, pps_hh, pm_sector1, pm_sector2) %>%
  mutate(pps_hh_sector1 = pps_hh*(pm_sector1/1000),
         sector_1_shares = pps_hh_sector1/sum(pps_hh_sector1),
         pps_hh_sector2 = pps_hh*(pm_sector2/1000),
         sector_2_shares = pps_hh_sector2/sum(pps_hh_sector2)) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile", 
                    pps_hh = "pps hh",
                    pm_sector1 = "s1 (pm)",
                    pm_sector2 = "s2 (pm)",
                    pps_hh_sector1 = "s1 (pps hh)",
                    sector_1_shares = "s1 (q share)",
                    pps_hh_sector2 = "s2 (pps hh)",
                    sector_2_shares = "s2 (q share)") %>%
  fit_to_width(max_width = 6.8)

Step 2) Join the HBS income quintile shares per sector to the EE-MRIO by sector, and multiply these shares by the EE-MRIO household final demand expenditure per sector ('eemrio hh fd') to stratify it by income quintile (Table S2).

q_share_of_sector = hbs %>%
  select(quintile,pps_hh_sector1,pps_hh_sector2) %>%
  gather(sector,value,-quintile) %>%
  mutate(sector = dplyr::recode(sector,
                                "pps_hh_sector1" = 1,
                                "pps_hh_sector2" = 2)) %>%
  group_by(sector) %>%
  mutate(q_share_of_sector = value/sum(value)) %>%
  select(-value) %>%
  spread(quintile,q_share_of_sector) 

eemrio_hh_fd = c(300,800)

eemrio = data.frame(q_share_of_sector, eemrio_hh_fd) %>%
  mutate(q1_eemrio = q1*eemrio_hh_fd,
         q2_eemrio = q2*eemrio_hh_fd,
         q3_eemrio = q3*eemrio_hh_fd,
         q4_eemrio = q4*eemrio_hh_fd,
         q5_eemrio = q5*eemrio_hh_fd) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(eemrio) %>%
  autofit(add_w=0.5) %>%
  set_header_labels(sector = "sector",
                    q1 = "q1 share",
                    q2 = "q2 share",
                    q3 = "q3 share",
                    q4 = "q4 share",
                    q5 = "q5 share",
                    eemrio_hh_fd = "eemrio hh fd",
                    q1_eemrio = "q1 fd",
                    q2_eemrio = "q2 fd",
                    q3_eemrio = "q3 fd",
                    q4_eemrio = "q4 fd",
                    q5_eemrio = "q5 fd") %>%
  fit_to_width(max_width = 6.8)

Step 3) Once the EE-MRIO household final demand expenditure is stratified by income quintile, it is multiplied by the EE-MRIO 'direct and indirect supply chain' intensities per sector to calculate the EE-MRIO household footprint per sector stratified by income quintile (Table S3). The total footprint per sector ('total fp') is the same as if we had simply multiplied the non-stratified EE-MRIO household final demand expenditure per sector ('eemrio hh fd') by the 'direct and indirect supply chain' intensity vector ('TIV').

TIV = c(0.2,0.4)

footprint = data.frame(eemrio,TIV) %>%
  select(sector,
         eemrio_hh_fd,
         q1_eemrio,
         q2_eemrio,
         q3_eemrio,
         q4_eemrio,
         q5_eemrio,
         TIV) %>%
  mutate(q1_footprint = q1_eemrio*TIV,
         q2_footprint = q2_eemrio*TIV,
         q3_footprint = q3_eemrio*TIV,
         q4_footprint = q4_eemrio*TIV,
         q5_footprint = q5_eemrio*TIV,
         total_footprint = eemrio_hh_fd*TIV) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(footprint) %>%
  autofit() %>%
  set_header_labels(sector = "sector",
                    eemrio_hh_fd = "eemrio hh fd",
                    q1_eemrio = "q1 fd",
                    q2_eemrio = "q2 fd",
                    q3_eemrio = "q3 fd",
                    q4_eemrio = "q4 fd",
                    q5_eemrio = "q5 fd",
                    TIV = "TIV",
                    q1_footprint = "q1 fp",
                    q2_footprint = "q2 fp",
                    q3_footprint = "q3 fp",
                    q4_footprint = "q4 fp",
                    q5_footprint = "q5 fp",
                    total_footprint = "total fp") %>%
  fit_to_width(max_width = 6.8)

Correspondence between EXIOBASE and EUROSTAT HBS

Table S4 shows which EUROSTAT HBS consumption category each EXIOBASE production sector corresponds to, along with our grouping of the EXIOBASE production sectors into the five aggregated final consumption categories we present in the main paper. The aggregated final consumption categories of 'goods' and 'services' refer to other material goods and other services that are not included in the three aggregated categories of food, housing and mobility. The direct energy use and carbon from households correspond to: heating/cooling (HH_HEAT) to CP045, mobility (HH_TRA) to CP072, and other (HH_OTHER) to CP05, respectively.

labels = read.csv(here("analysis/preprocessing/income-stratified-footprints/Exiobase_T_labels_ixi_w_coicop_mapping.csv")) %>%
  select(V2,coicop,five_sectors) %>%
  unique()

flextable(labels) %>%
  autofit() %>%
  set_header_labels(V2 = "exiobase industry production sector",
                    coicop = "coicop consumption category",
                    five_sectors = "aggregate final consumption category") %>%
  width(j=1, width = 4) %>%
  width(j=2, width = 1.5) %>%
  width(j=3, width = 1.5) 

Country and year coverage

The EXIOBASE version3.7 industry-by-industry is available for the years 1995 to 2016, albeit with the caveat that the original data series ends in 2011, the 2012-2016 estimates are based on trade and macro-economic data, and caution is advised when analysing trends over time. The EUROSTAT HBS is available for the years: 1988, 1994, 1999, 2005, 2010 and 2015, although not all countries are available for all years, and the years 1988, 1994, and 1999, in particular, have many missing countries. Table S5 shows the country and year coverage between EXIOBASE and the EUROSTAT HBS for the years 2005, 2010 and 2015. Rows with black text show countries that are represented in EXIOBASE, and an 'x' for those years where the necessary EUROSTAT HBS data also exists for that country. Rows with red text show countries where EUROSTAT HBS data exists, but who are not represented individually in EXIOBASE (they are in 'rest-of-world' categories in EXIOBASE). In Figure 3 of the paper, the 'non-Europe' distinction refers to all countries outside our sample (except Italy and Luxembourg), but is imperfect as there are some smaller countries on the European continent included there: Switzerland, Iceland, several Balkan states, and all of the European microstates (except Malta). Each of these countries either has no EUROSTAT HBS data at all, or is not included as a stand-alone country in EXIOBASE. Italy and Luxembourg are included in the 'Europe' distinction (even though they are outside of our sample for our empirical results), as they have EUROSTAT HBS data, are stand-alone countries in EXIOBASE, but did not have EUROSTAT HBS data disaggregated by income quintile for all three years of 2005, 2010 and 2015, although this could change in future EUROSTAT HBS waves. The 'Europe', 'non-Europe' distinction was used for simplicity in the figure.

geo = c("Austria",
        "Belgium",
        "Bulgaria",
        "Cyprus",
        "Czech Rep.",
        "Germany",
        "Denmark",
        "Estonia",
        "Greece",
        "Spain",
        "Finland",
        "France",
        "Croatia",
        "Hungary",
        "Ireland",
        "Italy",
        "Lithuania",
        "Luxembourg",
        "Latvia",
        "Montenegro",
        "North Macedonia",
        "Malta",
        "Netherlands",
        "Norway",
        "Poland",
        "Portugal",
        "Romania",
        "Serbia",
        "Sweden",
        "Slovenia",
        "Slovakia",
        "Turkey",
        "UK",
        "Kosovo")

year_2015 = c("x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x")

year_2010 = c("x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "",
              "x",
              "",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "",
              "x",
              "x",
              "x",
              "x",
              "x",
              "")
  
year_2005 = c("x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "x",
              "",
              "x",
              "x",
              "x",
              "x",
              "x",
              "")

country_year_coverage = data.frame(geo,
                                   year_2015,
                                   year_2010,
                                   year_2005)

flextable(country_year_coverage) %>%
  autofit() %>%
  set_header_labels(geo = "Country",
                    year_2015 = "2015",
                    year_2010 = "2010",
                    year_2005 = "2005") %>%
  color(i=c(20,21,28,34),color='red') %>%
  fit_to_width(max_width = 7.5) 

Purchaser price versus base price

The EUROSTAT HBS is reported in purchaser price (consumers report their expenditure in the prices they paid) while we use EXIOBASE version3.7 in base price. Because we decided to keep the EXIOBASE production sector shares the same (and just stratify by income quintile within each production sector), which also keeps the total footprint the same, this distinction between purchaser price and base price does not matter. The 'alternative method' section of this document shows that, in the alternative methodology where we break the EXIOBASE production sector shares and keep the consumption category shares of the HBS, the distinction between purchaser price and base price does matter.

In the methodology used for the main paper, once we have multiplied HBS 'pps hh' per income quintile by the sectoral shares in 'pm', to get sectoral expenditure per income quintile in 'pps hh', we include the fact that sector 1 and sector 2 have base price to purchaser price ratios that are different. The base price to purchaser price ratio for sector 1 ('s1 bp') is 0.5, whereas for sector 2 ('s2 bp') it is 1 (ie. there are no trade, transport, tax or subsidy margins in sector 2). This results in sectoral expenditure per income quintile in 'pps hh' now expressed in base price ('s1 pps hh bp' and 's2 pps hh bp') (Table S6).

bp_share_in_pp_sector_1 = c(0.5,0.5,0.5,0.5,0.5)
bp_share_in_pp_sector_2 = c(1,1,1,1,1)

hbs_bp = data.frame(quintile, 
                    pps_hh, 
                    pm_sector1, 
                    pm_sector2, 
                    bp_share_in_pp_sector_1, 
                    bp_share_in_pp_sector_2) %>%
  mutate(pps_hh_sector1 = pps_hh*(pm_sector1/1000),
         pps_hh_sector1_bp = pps_hh_sector1*bp_share_in_pp_sector_1,
         sector_1_shares = pps_hh_sector1/sum(pps_hh_sector1),
         sector_1_shares_bp = pps_hh_sector1_bp/sum(pps_hh_sector1_bp),
         pps_hh_sector2 = pps_hh*(pm_sector2/1000),
         pps_hh_sector2_bp = pps_hh_sector2*bp_share_in_pp_sector_2,
         sector_2_shares = pps_hh_sector2/sum(pps_hh_sector2),
         sector_2_shares_bp = pps_hh_sector2_bp/sum(pps_hh_sector2_bp),
         pps_hh_bp = pps_hh_sector1_bp + pps_hh_sector2_bp)
  
hbs_bp_pps = hbs_bp %>%
  select(quintile, 
         pps_hh,
         pm_sector1,
         pm_sector2,
         pps_hh_sector1,
         bp_share_in_pp_sector_1,
         pps_hh_sector1_bp,
         pps_hh_sector2,
         bp_share_in_pp_sector_2,
         pps_hh_sector2_bp) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs_bp_pps) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile",
                    pps_hh = "pps hh",
                    pm_sector1 = "s1 (pm)",
                    pm_sector2 = "s2 (pm)",
                    pps_hh_sector1 = "s1 (pps hh)",
                    bp_share_in_pp_sector_1 = "s1 bp",
                    pps_hh_sector1_bp = "s1 (pps hh bp)",
                    pps_hh_sector2 = "s2 (pps hh)",
                    bp_share_in_pp_sector_2 = "s2 bp",
                    pps_hh_sector2_bp = "s2 (pps hh bp)") %>%
  fit_to_width(max_width = 6.8)

Table S7 shows that even as the sectoral expenditure for sector 1 is 50% lower in base price than in purchaser price, because we only have one base price to purchaser price ratio per sector (ie. we assume that every income quintile in a country faces the same base price to purchaser price ratio when purchasing a good or service from a given sector), the sectoral shares per income quintile in base price are the same as they were in Table S1.

hbs_bp_shares = hbs_bp %>%
  select(quintile, 
         pps_hh_bp,
         pps_hh_sector1_bp,
         sector_1_shares_bp,
         pps_hh_sector2_bp,
         sector_2_shares_bp) %>%
  mutate(pm_sector1_bp = (pps_hh_sector1_bp/(pps_hh_sector1_bp + pps_hh_sector2_bp))*1000,
         pm_sector2_bp = (pps_hh_sector2_bp/(pps_hh_sector1_bp + pps_hh_sector2_bp))*1000) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs_bp_shares) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile",
                    pps_hh_bp = "pps hh bp",
                    pps_hh_sector1_bp = "s1 (pps hh bp)",
                    sector_1_shares_bp = "s1 (q share bp)",
                    pps_hh_sector2_bp = "s2 (pps hh bp)",
                    sector_2_shares_bp = "s2 (q share bp)",
                    pm_sector1_bp = "s1 (pm bp)",
                    pm_sector2_bp = "s2 (pm bp)") %>%
  fit_to_width(max_width = 6.8)

If we now do the same as in Table S3, joining these income quintile shares per sector to the EE-MRIO by sector, and multiplying them by the EE-MRIO household final demand expenditure per sector to stratify it by income quintile, we end up with the same stratification which will result in the same stratified footprint as before.

q_share_of_sector_bp = hbs_bp_shares %>%
  select(quintile,pps_hh_sector1_bp,pps_hh_sector2_bp) %>%
  gather(sector,value,-quintile) %>%
  mutate(sector = dplyr::recode(sector,
                                "pps_hh_sector1_bp" = 1,
                                "pps_hh_sector2_bp" = 2)) %>%
  group_by(sector) %>%
  mutate(q_share_of_sector_bp = value/sum(value)) %>%
  select(-value) %>%
  spread(quintile,q_share_of_sector_bp) 

eemrio_hh_fd = c(300,800)

eemrio_bp = data.frame(q_share_of_sector_bp, eemrio_hh_fd) %>%
  mutate(q1_eemrio = q1*eemrio_hh_fd,
         q2_eemrio = q2*eemrio_hh_fd,
         q3_eemrio = q3*eemrio_hh_fd,
         q4_eemrio = q4*eemrio_hh_fd,
         q5_eemrio = q5*eemrio_hh_fd) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(eemrio_bp) %>%
  autofit() %>%
  set_header_labels(sector = "sector",
                    q1 = "q1 share bp",
                    q2 = "q2 share bp",
                    q3 = "q3 share bp",
                    q4 = "q4 share bp",
                    q5 = "q5 share bp",
                    eemrio_hh_fd = "eemrio hh fd",
                    q1_eemrio = "q1 fd",
                    q2_eemrio = "q2 fd",
                    q3_eemrio = "q3 fd",
                    q4_eemrio = "q4 fd",
                    q5_eemrio = "q5 fd") %>%
  fit_to_width(max_width = 6.8)

European household expenditure deciles

Each national income quintile has a household final demand expenditure, household energy footprint and household carbon footprint estimate allocated to it, which represents the average in the quintile, after these initial steps. Then, to calculate European household expenditure deciles, we first ranked all these national income quintiles (140 in total: 28 European countries x 5 national income quintiles each) according to their mean household expenditure in PPS and aggregated the result to 10 European expenditure groups. This distribution allowed us to analyze the total European household energy and carbon footprints per these European expenditure deciles. We included only those countries with EUROSTAT HBS and EXIOBASE data in 2015, 2010, and 2005, which excludes Italy (no 2010 or 2015 necessary EUROSTAT HBS data, i.e. no data per income quintile) and Luxembourg (no 2010 EUROSTAT HBS data), but includes Norway and Turkey.

Each national income quintile is thus allocated to one of the 10 European expenditure deciles (some national income quintiles at the boundaries between deciles are split between two deciles). Figure S1 shows the population share of each country in our bottom 4 European expenditure deciles in 2015. A 100% share thus means that all 5 national income quintiles of that country fall within the bottom 4 European expenditure deciles. This does not imply that there are no high-income households in those countries, but because this method uses average expenditure data from the national income quintiles, the aggregation cuts off the lower and higher tails of the respective national expenditure distributions.

# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

tmp3 = pdat_country_summary_by_cquint_and_euntile %>%
  filter(eu_q_rank < 5, year== 2015,  indicator=="pae_fd_ke") %>%
  group_by(iso3, quint) %>%
  summarise(value = sum(value)) %>%
  group_by(iso3) %>%
  summarise(pop = n()/5*100) 

tmp4 = pdat_country_summary_by_cquint_and_euntile %>%
  ungroup() %>%
  filter(quint == 1, year== 2015,  indicator=="pae_fd_ke") %>%
  transmute(iso3, pop2 = 0) %>%
  left_join(tmp3, by="iso3") %>%
  transmute(iso3, popshare = if_else(is.na(pop), pop2, pop+pop2))

library(rworldmap)
library(ggthemes)

quantile_rank_map = joinCountryData2Map(tmp4, joinCode = "ISO_A3", nameJoinColumn = "iso3")
quantile_rank_map_poly = fortify(quantile_rank_map) #extract polygons 

quantile_rank_map_poly = merge(quantile_rank_map_poly, quantile_rank_map@data, by.x="id", by.y="ADMIN", all.x=T)
quantile_rank_map_poly = quantile_rank_map_poly %>% arrange(id, order)

map1 = ggplot() + 
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group), color="black") +
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group, 
                                        fill=factor(popshare))) + 
  #scale_fill_gradient(name="Mean decile \nrank", low="skyblue", high = "skyblue4", na.value = "#EEEEEE") +
  scale_fill_viridis(direction = -1, discrete = T, na.value = "white", name="Population\n share in\n bottom 4\n European\n deciles (%)") +
  theme_map() +
  coord_map("bonne", lat0 = 50,xlim = c(-9, 44), ylim = c(35, 70), clip="on") +
  theme(legend.position = c(0.75, 0.28))
  

map1

Limitations

While the EUROSTAT HBS is compiled for cross-country comparison purposes and aims for harmonization, there remains imperfect harmonization in the frequency of surveys, timing, content and structure between countries and years [@eurostat_eu_2020]. Some types of households may also be excluded from the samples, including super-rich households, for example in Germany, which excludes households with over €18,000 monthly net income [@eurostat_eu_2020]. Sensitive goods and services, such as alcohol, may be under-reported in household budget surveys, while expenditure on infrequent purchases such as a vehicle may create artificially large expenditure differences between households depending on the timing of the survey [@eurostat_eu_2020]. The EUROSTAT HBS macro-data also does not report direct foreign purchases, and we assumed that the expenditure shares between income quintiles of direct final demand for foreign goods and services was the same as direct final demand for domestic goods and services.

We used EXIOBASE as the EE-MRIO for this study because of its European focus, with nearly all countries in the EUROSTAT HBS also found as stand-alone countries in EXIOBASE, its detailed satellite extension data, and its year coverage (specifically version3.7, industry-by-industry), but there are well known limitations when using and selecting an EE-MRIO [@moran_convergence_2014]. The production sectors in EXIOBASE are harmonized across countries and years, but needing to map the EUROSTAT HBS to EXIOBASE meant that the most recent year of 2015 could only use the industry-by-industry version of EXIOBASE version3.7. This version assumes fixed product sales. Furthermore, because EXIOBASE version3.7 is extrapolated beyond 2011, caution should be used when comparing results over time. This, and the fact that harmonization guidelines in the EUROSTAT HBS have changed over time, were the justification for presenting only the 2015 results in the main paper, and presenting 2005 and 2010 results only here in the SI. We also show 2010 results using the product-by-product version of EXIOBASE version3.4 in the final section of this SI document.

The correspondence between the EUROSTAT HBS COICOP consumption categories and the industry production sectors in EXIOBASE is not one-to-one. Both the EUROSTAT HBS and EXIOBASE are limited in their consumption category/production sector level of detail. The share of each consumption category/production sector in the total amount of expenditure is also not identical between the HBS and EXIOBASE. As discussed in the 'methods for main paper results' section of this SI document, there are alternative methods for stratifying the EXIOBASE household final demand expenditure: one that keeps the EXIOBASE production sector shares of total expenditure intact, and one that keeps the HBS consumption category shares of total expenditure intact. Our results in the main paper use the first method, keeping the EXIOBASE sectoral shares of total expenditure intact, which means that the total footprint is identical to when it is calculated in EXIOBASE without any stratification by income quantile. The alternative method, on the other hand, results in a different total footprint because a different amount of final demand expenditure in each sector is multiplied by the same original 'direct and indirect supply chain' intensities, but stays faithful to the original HBS consumption category shares of total HBS expenditure. We show the alternative method below, and some results in the last sections of this SI document. We also stay at a relatively aggregated level for our energy and carbon footprints, as our primary goal was to connect aspects of the aggregate European household environmental footprint distribution with the decarbonisation and minimum energy use scenarios. Work specifically investigating distributional aspects of the European carbon footprint, sometimes at a finer sectoral level than we do here, can be found in the refs. [@ivanova_unequal_2020; @gore_t._confronting_2020; @sommer_carbon_2017] and [@bianco_understanding_2019].

Finally, the main limitation of the European expenditure deciles is the fact that they use average expenditure data from the national income quintiles. This aggregation cuts off the lower and higher tails of the respective national expenditure distributions. Our constructed European expenditure deciles thus likely underestimate European expenditure inequality. Also, when calculating energy and carbon intensities per harmonized European expenditure decile, and comparing them across deciles, one aspect not harmonized is the possible effect of different national energy subsidy levels. Varying energy subsidy levels between countries could play a role in the different energy and carbon intensities we observe per European expenditure deciles [@sovacool_reviewing_2017].

Alternative method for income-stratified national household energy and carbon footprints

Our methodology used for the main paper, and explained in the sections above, keeps the production sector shares of EE-MRIO household final demand expenditure (and subsequently the footprint) the same as they are found in the original EE-MRIO household final demand expenditure when not stratified by income quantile. The alternative method is to keep the consumption category shares of total HBS expenditure the same as they are found in the HBS. This means taking the total sum of household final demand expenditure from the EE-MRIO and stratifying it first based on the share of each income quantile's consumption expenditure of the total consumption expenditure as found in the HBS, before stratifying into sectors as well using the HBS 'parts per mille' per sector. This leads to a different total footprint than the original EE-MRIO footprint (when not stratified by income quantile), because a different amount of final demand expenditure in each sector is now multiplied by the same 'direct and indirect supply chain' intensities per sector that are originally calculated in the EE-MRIO.

Here the distinction between purchaser price and base price matters, as we directly take the HBS expenditure structure (which is reported in purchaser price) and apply it to the aggregated sum of EE-MRIO household final demand expenditure which is originally in base price. To correct for this, we can either convert the household budget survey to base price or take the total EE-MRIO household final demand expenditure in purchaser price, stratify it according to the HBS structure, then convert to base price. Here we show how we transform the household budget survey from purchaser price to base price using EUROSTAT margin data tables.

We use two data tables from the OECD statistics database, which are provided to the OECD by EUROSTAT (for European countries). The units here are: 'Trade and transport margins in percentage of final consumption expenditure by households' and 'Taxes less subsidies on product in percentage final consumption expenditure by households':

• Trade and Transport margin tables, Taxes less Subsidies tables (National Accounts > Annual National Accounts > Supply and Use Tables > SUT Indicators): https://stats.oecd.org/index.aspx?queryid=84864# [@oecd_sut_2021]

  1. Trade and transport margins in percentage of final consumption expenditure by households: Export as text file (csv) - customize to include years 2005 to 2015 [accessed on: 15-04-2020]

  2. Taxes less subsidies on product in percentage of final consumption expenditure by households: Export as text file (csv) - customize to include years 2005 to 2015 [accessed on: 15-04-2020]

We first map the HBS COICOP consumption categories to the margin data 'PRODUCT' sectors. Most countries have margin data here starting only in 2010. We take the 2010 margin data and apply this to all three years in our sample (the margins do not change significantly over time). This was possible for 20 countries. There are then two countries with only 'rent' margins missing: the Czech Republic and Slovakia. Here we apply zero to both trade and transport, and taxes less subsidies margins. There are then five countries with some missing margin data (Estonia, Ireland, Lithuania, Norway and Sweden), and three countries with no margin data at all (Germany, Spain, Turkey). We impute their margin percentages using the margin percentages of a similar, neighboring country. For Estonia and Lithuania, we use the Latvian margin percentages. For Ireland we use the UK margin percentages. For Norway and Sweden we use the Finnish margin percentages. For Germany we use the Austrian margin percentages. For Spain we use the Portuguese margin percentages, and finally, for Turkey we use the Bulgarian margin percentages. These countries appear to be close approximations based on examination of the 2010 values of another table; 'trade and transport margins in percentage of total supply at purchaser's prices.' The need to apply margin percentage data is a potential limitation of this method. If we relied on this method to produce main paper results, acquiring and applying more precise margin percentage data would arguably be a requirement.

bp_share_in_pp_sector_1 = c(0.5,0.5,0.5,0.5,0.5)
bp_share_in_pp_sector_2 = c(1,1,1,1,1)

hbs_bp = data.frame(quintile, 
                    pps_hh, 
                    pm_sector1, 
                    pm_sector2, 
                    bp_share_in_pp_sector_1, 
                    bp_share_in_pp_sector_2) %>%
  mutate(pps_hh_sector1 = pps_hh*(pm_sector1/1000),
         pps_hh_sector1_bp = pps_hh_sector1*bp_share_in_pp_sector_1,
         sector_1_shares = pps_hh_sector1/sum(pps_hh_sector1),
         sector_1_shares_bp = pps_hh_sector1_bp/sum(pps_hh_sector1_bp),
         pps_hh_sector2 = pps_hh*(pm_sector2/1000),
         pps_hh_sector2_bp = pps_hh_sector2*bp_share_in_pp_sector_2,
         sector_2_shares = pps_hh_sector2/sum(pps_hh_sector2),
         sector_2_shares_bp = pps_hh_sector2_bp/sum(pps_hh_sector2_bp),
         pps_hh_bp = pps_hh_sector1_bp + pps_hh_sector2_bp)
  
hbs_bp_pps = hbs_bp %>%
  select(quintile, 
         pps_hh,
         pm_sector1,
         pm_sector2,
         pps_hh_sector1,
         bp_share_in_pp_sector_1,
         pps_hh_sector1_bp,
         pps_hh_sector2,
         bp_share_in_pp_sector_2,
         pps_hh_sector2_bp) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs_bp_pps) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile",
                    pps_hh = "pps hh",
                    pm_sector1 = "s1 (pm)",
                    pm_sector2 = "s2 (pm)",
                    pps_hh_sector1 = "s1 (pps hh)",
                    bp_share_in_pp_sector_1 = "s1 bp",
                    pps_hh_sector1_bp = "s1 (pps hh bp)",
                    pps_hh_sector2 = "s2 (pps hh)",
                    bp_share_in_pp_sector_2 = "s2 bp",
                    pps_hh_sector2_bp = "s2 (pps hh bp)") %>%                 
  fit_to_width(max_width = 7)

hbs_bp_shares = hbs_bp %>%
  select(quintile, 
         pps_hh_bp,
         pps_hh_sector1_bp,
         pps_hh_sector2_bp) %>%
  mutate(pm_sector1_bp = (pps_hh_sector1_bp/(pps_hh_sector1_bp + pps_hh_sector2_bp))*1000,
         pm_sector2_bp = (pps_hh_sector2_bp/(pps_hh_sector1_bp + pps_hh_sector2_bp))*1000) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs_bp_shares) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile",
                     pps_hh_bp = "pps hh bp",
                     pps_hh_sector1_bp = "s1 (pps hh bp)",
                     pps_hh_sector2_bp = "s2 (pps hh bp)",
                     pm_sector1_bp = "s1 (pm bp)",
                     pm_sector2_bp = "s2 (pm bp)") %>%
  fit_to_width(max_width = 7.5)

Table S11 shows the HBS data on the left-hand side ('pps hh', 's1 (pm)' and 's2 (pm)') as before, but now the income quintile shares of mean consumption expenditure (the sum of 'pps hh') are calculated and shown as 'mean exp'. These shares are multiplied by total EE-MRIO household final demand expenditure, which in this case is 1100 (300 in sector 1 plus 800 in sector 2) (see Table S2), to calculate the amount of EE-MRIO household final demand expenditure per income quintile. These are then multiplied by the sectoral parts per mille shares to calculate EE-MRIO household final demand expenditure per sector and income quintile ('hh fd s1' and 'hh fd s2').

quintile = c("q1","q2","q3","q4","q5")
pps_hh = c(9,11,16,25,28)

pm_sector1 = c(111,176,250,333,429)
pm_sector2 = c(889,824,750,667,571)

hbs_alt_method = data.frame(quintile, 
                            pps_hh, 
                            pm_sector1, 
                            pm_sector2,
                            bp_share_in_pp_sector_1,
                            bp_share_in_pp_sector_2) %>%
  mutate(mean_expenditure_share = pps_hh/sum(pps_hh),
         eemrio_hh_fd = 1100*mean_expenditure_share,
         hh_sector1 = eemrio_hh_fd*(pm_sector1/1000),
         hh_sector1_bp = hh_sector1*bp_share_in_pp_sector_1,
         hh_sector2 = eemrio_hh_fd*(pm_sector2/1000),
         hh_sector2_bp = hh_sector2*bp_share_in_pp_sector_2)

hbs_alt_method_fd = hbs_alt_method %>%
  select(quintile,
         pps_hh,
         pm_sector1,
         pm_sector2,
         mean_expenditure_share,
         eemrio_hh_fd,
         hh_sector1,
         hh_sector2) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(hbs_alt_method_fd) %>%
  autofit() %>%
  set_header_labels(quintile = "quintile",
                    pps_hh = "pps hh",
                    pm_sector1 = "s1 (pm)",
                    pm_sector2 = "s2 (pm)",
                    mean_expenditure_share = "mean exp",
                    eemrio_hh_fd = "eemrio hh fd",
                    hh_sector1 = "hh fd s1",
                    hh_sector2 = "hh fd s2") %>%
  fit_to_width(max_width = 7.5)

Table S12 shows that we then have our EE-MRIO household final demand expenditure stratified by income quintile.

eemrio_alt_method = hbs_alt_method %>%
  select(quintile,hh_sector1,hh_sector2) %>%
  gather(sector,value,-quintile) %>%
  spread(quintile,value) %>%
  mutate(sector = dplyr::recode(sector,
                                "hh_sector1" = "1",
                                "hh_sector2" = "2")) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(eemrio_alt_method) %>%
  autofit() %>%
  set_header_labels(sector = "sector",
                    q1 = "q1 fd",
                    q2 = "q2 fd",
                    q3 = "q3 fd",
                    q4 = "q4 fd",
                    q5 = "q5 fd") %>%
  width(j=c(1:6), width = 0.8) 

If we now multiply this by the 'direct and indirect supply chain' intensity vector ('TIV'), we calculate the footprint stratified by income quintile, but our total footprint is different than in the other methodology because now different amounts of EE-MRIO household final demand expenditure per sector are being multiplied by the same 'direct and indirect supply chain' intensity vector as previously calculated inside the EE-MRIO (compare Table S13 with Table S3).

footprint_alt_method = data.frame(eemrio_alt_method, TIV) %>%
  mutate(q1_footprint = q1*TIV,
         q2_footprint = q2*TIV,
         q3_footprint = q3*TIV,
         q4_footprint = q4*TIV,
         q5_footprint = q5*TIV,
         total_footprint = q1_footprint +
           q2_footprint +
           q3_footprint +
           q4_footprint +
           q5_footprint) %>%
  mutate_if(is.numeric, round, digits = 2)

flextable(footprint_alt_method) %>%
  autofit() %>%
  set_header_labels(sector = "sector",
                    q1 = "q1 fd",
                    q2 = "q2 fd",
                    q3 = "q3 fd",
                    q4 = "q4 fd",
                    q5 = "q5 fd",
                    TIV = "TIV",
                    q1_footprint = "q1 fp",
                    q2_footprint = "q2 fp",
                    q3_footprint = "q3 fp",
                    q4_footprint = "q4 fp",
                    q5_footprint = "q5 fp",
                    total_footprint = "total fp") %>%
  fit_to_width(max_width = 7.5)

Supplementary Results

We ran the same analysis for the years 2015, 2010 and 2005. We presented only 2015 in the main paper as we focused more on the distributional aspect of the environmental footprints and their relationship to decarbonisation scenarios and minimum energy use, rather than the time component. As mentioned in the 'limitations' section above, caution is needed when comparing the EUROSTAT HBS data over time, as well as the version of EXIOBASE we used (version3.7, industry-by-industry). We also estimated the environmental footprints per European expenditure decile using the product-by-product version of EXIOBASE for the year 2010, and then the year 2015 again using the alternative methodology that was explained in the supplementary materials and methods sections in this document above.

We show each of these in turn in the format of Figure 1 from the main paper: 1) the year 2005 using the main paper method and EXIOBASE industry-by-industry version, then 2) the year 2010 using the main paper method and EXIOBASE industry-by-industry version, 3) the year 2010 using the main paper method but EXIOBASE product-by-product version, and finally 4) the year 2015 using the alternative methodology and EXIOBASE industry-by-industry version.

Year 2005, using main method, EXIOBASE industry-by-industry

# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator %in% c("total_energy_use_tj", "total_energy_use_losses_tj")) %>%
  group_by(eu_ntile_name, indicator) %>%
  summarise(value = sum(value)*0.000001) %>%
  pivot_wider(id_cols = eu_ntile_name, names_from = indicator, values_from = value) %>% 
  mutate(total_energy_use_final_tj = total_energy_use_tj-total_energy_use_losses_tj) %>% 
  select(-total_energy_use_tj) %>% 
  pivot_longer(cols=-eu_ntile_name, names_to = "indicator", values_to = "value") %>% 
  mutate(component = if_else(indicator=="total_energy_use_losses_tj", "Loss", "Final energy")) %>% 
  mutate(component = factor(component, levels=c("Loss", "Final energy"))) %>% 
  ggplot(aes(x=eu_ntile_name, y=value, fill=component)) +
    geom_col(position = position_stack()) +
    scale_fill_manual(name = "", values = pal[c(3,1)]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Net energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  theme(legend.position = c(0.2,0.91),
        legend.text = element_text(size=11))

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2005) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2005) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2005) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2005) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2005) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2005) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (tCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

# values in text

## inequality
exp = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

exp_total = round((exp %>% summarise(value = sum(value)))$value,digits=1)

exp_10_10 = round((exp %>% filter(eu_q_rank == 10))$value/(exp %>% filter(eu_q_rank == 1))$value,digits = 1)

energy = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

energy_total = round((energy %>% summarise(value = sum(value)))$value,digits=1)

energy_10_10 = round((energy %>% filter(eu_q_rank == 10))$value/(energy %>% filter(eu_q_rank == 1))$value,digits = 1)

co2eq = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name))

co2eq_total = round((co2eq %>% summarise(value = sum(value)))$value,digits=1)

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2005,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2005) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2005) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

In 2005, total household final demand expenditure was r exp_total trn€, the energy footprint r energy_total EJ, and the carbon footprint r co2eq_total MtCO2eq. We estimated the 90:10 ratios at r exp_10_10 for expenditure, r energy_10_10 for the net energy footprint, and r co2eq_10_10 for the carbon footprint. Total expenditure ranged from r exp_bottom_decile trn€ to r exp_top_decile trn€ (or r fd_pae_bottom_decile€ to r fd_pae_top_decile€ per adult equivalent) across bottom and top decile, the net energy footprint from r energy_bottom_decile EJ to r energy_top_decile EJ (or r energy_pae_bottom_decile GJ/ae to r energy_pae_top_decile GJ/ae), and the carbon footprint from r co2eq_bottom_decile MtCO2eq to r co2eq_top_decile MtCO2eq (or r co2eq_pae_bottom_decile tCO2eq/ae to r co2eq_pae_top_decile tCO2eq/ae). The average net energy intensity of consumption was r mean_energy_intens_bottom_decile MJ/€ in the bottom decile and r mean_energy_intens_top_decile MJ/€ in the top decile. The carbon intensity of net energy was r mean_co2eq_of_energy_intens_bottom_decile tCO2eq/TJ in the bottom decile, and r mean_co2eq_of_energy_intens_top_decile tCO2eq/TJ in the top decile.

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') &
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))
a

Year 2010, using main method, EXIOBASE industry-by-industry

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator %in% c("total_energy_use_tj", "total_energy_use_losses_tj")) %>%
  group_by(eu_ntile_name, indicator) %>%
  summarise(value = sum(value)*0.000001) %>%
  pivot_wider(id_cols = eu_ntile_name, names_from = indicator, values_from = value) %>% 
  mutate(total_energy_use_final_tj = total_energy_use_tj-total_energy_use_losses_tj) %>% 
  select(-total_energy_use_tj) %>% 
  pivot_longer(cols=-eu_ntile_name, names_to = "indicator", values_to = "value") %>% 
  mutate(component = if_else(indicator=="total_energy_use_losses_tj", "Loss", "Final energy")) %>% 
  mutate(component = factor(component, levels=c("Loss", "Final energy"))) %>% 
  ggplot(aes(x=eu_ntile_name, y=value, fill=component)) +
    geom_col(position = position_stack()) +
    scale_fill_manual(name = "", values = pal[c(3,1)]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Net energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  theme(legend.position = c(0.2,0.91),
        legend.text = element_text(size=11))

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (tCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

# values in text

## inequality
exp = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

exp_total = round((exp %>% summarise(value = sum(value)))$value,digits=1)

exp_10_10 = round((exp %>% filter(eu_q_rank == 10))$value/(exp %>% filter(eu_q_rank == 1))$value,digits = 1)

energy = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

energy_total = round((energy %>% summarise(value = sum(value)))$value,digits=1)

energy_10_10 = round((energy %>% filter(eu_q_rank == 10))$value/(energy %>% filter(eu_q_rank == 1))$value,digits = 1)

co2eq = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name))

co2eq_total = round((co2eq %>% summarise(value = sum(value)))$value,digits=1)

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

In 2010, total household final demand expenditure was r exp_total trn€, the net energy footprint r energy_total EJ, and the carbon footprint r co2eq_total MtCO2eq. We estimated the 90:10 ratios at r exp_10_10 for expenditure, r energy_10_10 for the net energy footprint, and r co2eq_10_10 for the carbon footprint. Total expenditure ranged from r exp_bottom_decile trn€ to r exp_top_decile trn€ (or r fd_pae_bottom_decile€ to r fd_pae_top_decile€ per adult equivalent) across bottom and top decile, the net energy footprint from r energy_bottom_decile EJ to r energy_top_decile EJ (or r energy_pae_bottom_decile GJ/ae to r energy_pae_top_decile GJ/ae), and the carbon footprint from r co2eq_bottom_decile MtCO2eq to r co2eq_top_decile MtCO2eq (or r co2eq_pae_bottom_decile tCO2eq/ae to r co2eq_pae_top_decile tCO2eq/ae). The average net energy intensity of consumption was r mean_energy_intens_bottom_decile MJ/€ in the bottom decile and r mean_energy_intens_top_decile MJ/€ in the top decile. The carbon intensity of net energy was r mean_co2eq_of_energy_intens_bottom_decile tCO2eq/TJ in the bottom decile, and r mean_co2eq_of_energy_intens_top_decile tCO2eq/TJ in the top decile.

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') &
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

Year 2010, using main method, EXIOBASE product-by-product

# load data wrangling functions
source(here("analysis", "R", "si", "wrangler_functions_pxp.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

p_exp_and_co2eq_totals = p1 + p3

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (tCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

p_co2eq_intensity_of_consumption = p1

# values in text

## inequality
exp = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

exp_total = round((exp %>% summarise(value = sum(value)))$value,digits=1)

exp_10_10 = round((exp %>% filter(eu_q_rank == 10))$value/(exp %>% filter(eu_q_rank == 1))$value,digits = 1)

energy = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name))

energy_total = round((energy %>% summarise(value = sum(value)))$value,digits=1)

energy_10_10 = round((energy %>% filter(eu_q_rank == 10))$value/(energy %>% filter(eu_q_rank == 1))$value,digits = 1)

co2eq = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name))

co2eq_total = round((co2eq %>% summarise(value = sum(value)))$value,digits=1)

co2eq_10_10 = round((co2eq %>% filter(eu_q_rank == 10))$value/(co2eq %>% filter(eu_q_rank == 1))$value,digits = 1)

## total per decile

exp_bottom_decile = round((exp %>% filter(eu_q_rank == 1))$value, digits = 1)

exp_top_decile = round((exp %>% filter(eu_q_rank == 10))$value, digits = 1)

energy_bottom_decile = round((energy %>% filter(eu_q_rank == 1))$value, digits = 1)

energy_top_decile = round((energy %>% filter(eu_q_rank == 10))$value, digits = 1)

co2eq_bottom_decile = round((co2eq %>% filter(eu_q_rank == 1))$value, digits = 1)

co2eq_top_decile = round((co2eq %>% filter(eu_q_rank == 10))$value, digits = 1)


## per adult equivalent per decile

aeu = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2010,
         indicator == "total_adult_eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value),
            eu_ntile_name = first(eu_ntile_name)) 

exp_pae = exp %>%
  rename(total_fd_me = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_fd_e = (total_fd_me/value)*1000000000000)

fd_pae_bottom_decile = round((exp_pae %>% filter(eu_q_rank == 1))$pae_fd_e, digits = 0)

fd_pae_top_decile = round((exp_pae %>% filter(eu_q_rank == 10))$pae_fd_e, digits = 0)

energy_pae = energy %>%
  rename(total_energy_use_tj = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_energy_use_gj = (total_energy_use_tj/value)*1000000000)

energy_pae_bottom_decile = round((energy_pae %>% filter(eu_q_rank == 1))$pae_energy_use_gj, digits = 1)

energy_pae_top_decile = round((energy_pae %>% filter(eu_q_rank == 10))$pae_energy_use_gj, digits = 1)

co2eq_pae = co2eq %>%
  rename(total_co2eq_kg = value) %>%
  left_join(aeu, by = c("eu_q_rank", "eu_ntile_name")) %>%
  mutate(pae_co2eq_t = (total_co2eq_kg/value)*1000000)

co2eq_pae_bottom_decile = round((co2eq_pae %>% filter(eu_q_rank == 1))$pae_co2eq_t, digits = 1)

co2eq_pae_top_decile = round((co2eq_pae %>% filter(eu_q_rank == 10))$pae_co2eq_t, digits = 1)


## intensities

mean_energy_intens = dat_country_summary_by_eu_ntile %>%
                 filter(year == 2010) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me))

mean_energy_intens_bottom_decile = round((mean_energy_intens %>% filter(eu_q_rank == 1))$pe_energy_use_mj, digits = 1)

mean_energy_intens_top_decile = round((mean_energy_intens %>% filter(eu_q_rank == 10))$pe_energy_use_mj, digits = 1)

mean_co2eq_of_energy_intens = dat_country_summary_by_eu_ntile %>%
  filter(year == 2010) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj) %>%
  group_by(eu_q_rank) %>%
  summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj))

mean_co2eq_of_energy_intens_bottom_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 1))$intensity_e_c, digits = 1)

mean_co2eq_of_energy_intens_top_decile = round((mean_co2eq_of_energy_intens %>% filter(eu_q_rank == 10))$intensity_e_c, digits = 1)

In 2010 but using the product-by-product version of EXIOBASE, total household final demand expenditure was r exp_total trn€, and the carbon footprint r co2eq_total MtCO2eq. We estimated the 90:10 ratios at r exp_10_10 for expenditure, and r co2eq_10_10 for the carbon footprint. Total expenditure ranged from r exp_bottom_decile trn€ to r exp_top_decile trn€ (or r fd_pae_bottom_decile€ to r fd_pae_top_decile€ per adult equivalent) across bottom and top decile, and the carbon footprint from r co2eq_bottom_decile MtCO2eq to r co2eq_top_decile MtCO2eq (or r co2eq_pae_bottom_decile tCO2eq/ae to r co2eq_pae_top_decile tCO2eq/ae). The 'net energy use' energy extensions used for calculating the energy footprints in the main paper were not available for the product-by-product version used here, EXIOBASE version3.4.

#a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
#  plot_layout(guides = 'collect')  & 
#  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
#        legend.position = 'bottom',
#        axis.title.y = element_text(size=13, hjust = 0.5),
#        axis.text.x = element_text(size = 12),
#        axis.text.y = element_text(size = 12),
#        legend.text = element_text(size=12),
#        legend.title = element_text(size=13))

a = p_exp_and_co2eq_totals + p_co2eq_intensity_of_consumption + plot_annotation(tag_levels = 'a') +
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

Year 2015, using alternative method, EXIOBASE industry-by-industry

There is relatively good agreement between the alternative method footprints and the EXIOBASE footprints, except in Bulgaria, where the alternative method footprint was around 3 times larger. Including Bulgaria, the alternative method footprints were on average around 20% larger than the EXIOBASE footprints. With Bulgaria removed they were around 10% larger on average. Eastern European countries especially had larger alternative method footprints due to high intensities in electricity production and hot water supply, and then more expenditure in CP045 multiplied by these intensities than the expenditure in EXIOBASE. This is why the figure shows such high footprints in the bottom European deciles compared to the method keeping EXIOBASE footprints the same (ms results).

# load data wrangling functions
source(here("analysis", "R", "si", "wrangler_functions_method2_ixi.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "coicop")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "coicop")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_top = p1 + p2 + p3

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (tCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p_bottom = p2 + p1 + p3

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') &
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

a

Units of analysis

# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

# ae vs. per capita

dat_adult_eq = read_csv(here("analysis", "data", "derived", "si", "adult_eq_per_household.csv")) %>%
  filter(year==2015)

dat_tmp = dat_country_summary_by_cquint_and_euntile %>%
  filter(year == 2015) %>%
  left_join(dat_adult_eq %>%
              select(iso2, iso3, quint, adult_e_p_hh), by=c("iso2", "quint")) 

europe28 = c("AUT",
             "BEL",
             "BGR",
             "CYP",
             "CZE",
             "DEU",
             "DNK",
             "EST",
             "GRC",
             "ESP",
             "FIN",
             "FRA",
             "HRV",
             "HUN",
             "IRL",
             "LTU",
             "LVA",
             "MLT",
             "NLD",
             "NOR",
             "POL",
             "PRT",
             "ROU",
             "SWE",
             "SVN",
             "SVK",
             "TUR",
             "GBR")

population = wbstats::wb(country = europe28, 
                         indicator = "SP.POP.TOTL", 
                         startdate = 1990, enddate = 2017) %>%
  select(iso3 = iso3c, population = value, year = date) %>%
  mutate(year = parse_number(year))

dat_people_per_adult_eq = dat_tmp %>%
  ungroup() %>% 
  select(year, iso2, iso3, quint, total_adult_eq) %>%
  group_by(year, iso2, iso3) %>%
  summarise(total_adult_eq = sum(total_adult_eq)) %>%
  left_join(population, by = c("year","iso3")) %>%
  mutate(pop_per_ae = population/total_adult_eq) %>%
  ungroup() %>%
  select(iso2, pop_per_ae)

dat_tmp2 = dat_tmp %>%
  left_join(dat_people_per_adult_eq, by = "iso2") %>%
  mutate(pop_estimate_per_quintile = total_adult_eq*pop_per_ae) %>%
  ungroup() %>%
  select(iso2,iso3,quint,eu_q_rank,
                              total_adult_eq,
                              pop_estimate_per_quintile) %>%
  rename(total_population = pop_estimate_per_quintile)

eu_ntile_name = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")

a_tmp = dat_tmp2 %>%
  group_by(eu_q_rank) %>%
  summarise(total_population = sum(total_population)) %>%
  cbind(eu_ntile_name)

a = dat_tmp2 %>%
  group_by(eu_q_rank) %>%
  summarise(total_adult_eq = sum(total_adult_eq),
            total_population = sum(total_population)) %>%
  cbind(eu_ntile_name) %>%
  pivot_longer(cols = c(-eu_q_rank,-eu_ntile_name), names_to = "indicator", values_to="value") %>%
  mutate(indicator = dplyr::recode(indicator,
                                   "total_adult_eq" = "Adult equivalents",
                                   "total_population" = "Total population")) %>%
  mutate(value = value/1000000)

ae_vs_pop = dat_tmp2 %>%
  summarise(total_adult_eq = sum(total_adult_eq),
            total_population = sum(total_population)) 

ae_share_of_pop = round(ae_vs_pop$total_adult_eq/ae_vs_pop$total_population,digits=2)

Through the main paper we use household per adult equivalent as our unit of analysis, following the EUROSTAT HBS. This meant that we converted decarbonisation scenario final energy numbers from total per capita to household per adult equivalent to better compare them with our environmental footprint estimates. We converted them using 1) the household share of the total footprint, and 2) the adult equivalent share of the total population.

As a numerical example, we convert a total final energy of 53 GJ per capita from the Low Energy Demand (LED) scenario (Grubler et al. (2018) [@grubler_low_2018]), first by the household share of the total European final energy footprint in 2015 (around 0.65, calculated in EXIOBASE), and then the share of total adult equivalents in the total European population in 2015 (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 converted to a household final energy of 55 GJ/adult equivalent in Europe ((53 total GJ/capita * 0.65 household share of total footprint)/0.63 adult equivalent share of total population = 55 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 - Global North), and 15.3 GJ/capita (DLE) (115 GJ/capita - SSP2-2.6). Because of the similar relative shares of the household part of the total European energy footprint (around 0.65), and the adult equivalent share of the total population in our sample (around 0.63), these final energy numbers end up close to the same when converted 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; @iiasa_gea_2012]. 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 SSP numbers are from the MESSAGE-GLOBIOM model, and the GEA-efficiency number is from the MESSAGE model. The LED total GJ/capita numbers are 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 global (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 is 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.

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 ratios to the adult equivalents from different countries making up each European expenditure decile, 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 S6 shows these population estimates per European expenditure decile. We use this to re-estimate our household energy footprint per European expenditure decile in per capita terms, and then re-create Figure 5 from the main paper (Figure S7 below) in per capita terms.

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)) 

Figure 5 from main paper in household final energy per capita

Here we re-create Figure 5 from the main paper after estimating our household 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 convert 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.65 = 34 household GJ/capita.

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_energy_use_losses_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_losses_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]

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.65,
         fe_gj_pc = round(fe_gj_pc),
         ccs_required = round(ccs_required)) 

mea_seq = c(df_scenario_info$fe_gj_pc, seq(from=50, to=300, by=50))
#mea = c(16,300)
mer = c(10,60)

# 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.1)) {
  for (mean_energy in mea_seq) {
    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_extrafine.rds"))

round_by = 10

df_all = readRDS(here("analysis/data/derived/si/scenarios_extrafine.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_grid = df_all %>%
  filter(!(v_mean %in% df_scenario_info$fe_gj_pc))

df_scenario_info = df_scenario_info %>% 
  mutate(labl = paste0(scenario, " (", fe_gj_pc, ")"))

df_scenario = df_all %>%
  filter((v_mean %in% df_scenario_info$fe_gj_pc)) %>%
  left_join(df_scenario_info %>% select(labl, v_mean = fe_gj_pc))

df_scenario$labl = factor(df_scenario$labl, levels = c("DLE (10)",
                                                       "LED global (18)",
                                                       "LED North (34)",
                                                       "GEA-efficiency (42)",
                                                       "IEA ETP B2DS (55)",
                                                       "SSP1-1.9 (56)",
                                                       "SSP2-1.9 (61)",
                                                       "SSP2-2.6 (75)",
                                                       "SSP2-3.4 (85)",
                                                       "SSP2-4.5 (91)"))

#df_dle = tibble(bin_ratio = 1, v_first = 10, v_mean = 10)

#scenario_label = df_scenario_info %>%
#  mutate(labl = paste0(scenario, " (", fe_gj_pc, ")")) %>%
#  pull(labl)

df_hm = df_scenario %>% 
  group_by(v_mean) %>% 
  slice_min(v_first) %>% 
  mutate(dle = if_else(v_mean == 10, "yes", "no"))

library(ggsci)

a = df_grid %>%
  ggplot(aes(x=v_first, y=bin_ratio, group=v_mean)) +
  geom_smooth(aes(linetype="Maximum energy\nsupply (GJ/ae)"), se=FALSE, color="grey", size=0.5) +
  scale_linetype_manual(name="", values = c(2)) +
  scale_size_manual(values = c("yes"=1, "no"=-1))+
  geom_point(data=df_hm, aes(color=labl, size=dle)) +
  geom_smooth(data=df_scenario, aes(color=labl), se=FALSE) +
  #geom_point(data=df_scenario %>% filter(v_mean==16), aes(color=labl)) +
  geom_hline(yintercept = ineq_curr, color="grey") +
  annotate(geom="text",x=56,y=12,label="300", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=46,y=9.5,label="200", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=33.5,y=6,label="100", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=60,y=ineq_curr-0.28,
           label="energy use inequality 2015", size=3.75, hjust=1, color="grey40") +
  guides(size=FALSE) +
  scale_color_npg(name = "Scenario") +
  lims(x=c(9.5,60), y=c(1,13)) +
  labs(x= "Minimum energy use (GJ/ae)", y="Energy use inequality (90:10 ratio)") +
  theme_minimal() +
  theme(text=element_text(family="Liberation Sans Narrow"),
      axis.title.y = element_text(size=13, hjust = 0.5),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      legend.text = element_text(size=12),
      legend.title = element_text(size=13))

a

Formula for adjusting distribution

We first scale the vector \boldsymbol{e} = (e_1,...e_{10}) of mean final energy use values per decile to obtain a vector \boldsymbol{e_s} that conforms to the new target mean e_{mean} (e.g. from the energy supply scenario):

\boldsymbol{e_s} = \boldsymbol{e} \cdot e_{mean}/\hat{e}

where \hat{e} is the current mean final energy use across deciles. Next, we transform this vector \boldsymbol{e_s} to conform to the minimum final energy use constraint e_{min}, i.e. we set the first decile mean to the minimum final energy use value and scale all other decile means in a way that preserves their relative distance, the target mean e_{mean}, and the minimum energy use constraint in the first decile. We calculate the target decile values \boldsymbol{e_t} like this:

\boldsymbol{e_t} = e_{mean} - (e_{mean}-\boldsymbol{e_s}) \cdot \frac{\left(e_{mean}-e_{min}\right)}{\left(e_{mean}-e_{s,1}\right)}

library(viridis)
mer_const = 28
mea_const = 90


pdat_intensity = 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(final_energy_intensity = (total_energy_use_tj-total_energy_use_losses_tj)/(total_fd_me)) %>%  ## changed to subtract losses to get 'final' energy
  select(five_sectors, eu_q_rank, 
         final_energy_intensity) %>%
  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_intensity, by="five_sectors") %>%
  mutate(total_final_energy_use_tj_new = (total_fd_me)*final_energy_intensity) %>%
  #mutate(total_energy_use_tj_diff = (total_energy_use_tj-total_energy_use_losses_tj)-total_energy_use_tj_new) %>%  ## changed here too
  select(eu_q_rank,total_final_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_final_energy_use_tj_new = sum(total_final_energy_use_tj_new)) %>%
  mutate(pae_final_energy_use_tj = total_final_energy_use_tj_new/33417583,
         pae_final_energy_use_gj = pae_final_energy_use_tj*1000)

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

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



df_mea = NULL
for (min_energy in seq(from=16, to=112, by=10)) {
    if (min_energy <= mea_const) {
      df_mea = df_mea %>%
      bind_rows(df_energy_deciles %>%
                  scaled_quantiles(eu_q_rank, pae_final_energy_use_gj, min_energy, mea_const))
  }
}


p1 = df_mea %>% 
  mutate(exp_decile = if_else(eu_q_rank < 10, paste0("D0", eu_q_rank), paste0("D",eu_q_rank))) %>% 
  ggplot(aes(x=exp_decile, y=scaled)) +
  geom_point(aes(color=v_first)) +
  labs(x="Expenditure decile", y="Final energy footprint (GJ/ae)") +
  scale_color_viridis("Minimum\nenergy use\n(GJ/ae)")+
  ggtitle(paste0("Maximum inequality at ", mea_const, " GJ/ae\nmaximum energy supply")) +
  theme_minimal()


df_mer = NULL
for (mean_energy in seq(from=16, to=119, by=10)) {
    if (mean_energy >= mer_const) {
      df_mer = df_mer %>%
      bind_rows(df_energy_deciles %>%
                  scaled_quantiles(eu_q_rank, pae_final_energy_use_gj, mer_const, mean_energy))
    }
}

p2 = df_mer %>% 
  mutate(exp_decile = if_else(eu_q_rank < 10, paste0("D0", eu_q_rank), paste0("D",eu_q_rank))) %>% 
  ggplot(aes(x=exp_decile, y=scaled)) +
  geom_point(aes(color=v_mean)) +
  labs(x="Expenditure decile", y="Final energy footprint (GJ/ae)") +
  scale_color_viridis("Maximum energy\nsupply (GJ/ae)")+
  ggtitle(paste0("Maximum inequality at ", mer_const, " GJ/ae\nminimum energy use")) +
  theme_minimal() 

p1 + p2 + plot_annotation(tag_levels = 'a') &
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12, angle=90),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

\newpage

Supplementary references