Eurostat’s carbon footprint with FIGARO in R

Author

Pablo Piñero

Published

June 20, 2024

Loading EUROSTAT Carbon Footprint results

First, we will upload the EUROSTAT Carbon Footprint results for year 2021 to R (the most recent year available). The results are in format .csv. In this example, we will paste the path and name of the file using the paste0() function. This way, if we store the data in a different location, we just need to change the path in the pathToData object. We will call the results CF.

pathToData <- "/home/NET1/b7_training/"

CF <- readr::read_csv(paste0(pathToData, "CO2_footprints_2021.csv"))

By using the head() function, we can observe the data structure. The CF object created is a data.frame with 12 columns, containing information on the year, country where the emissions are occurring (ref_area), industry, and the final consumer country (counterpart_area) and category, and other attributes:

head(CF)
# A tibble: 6 × 12
  time_period ref_area industry counterpart_area sto    obs_value decimals
        <dbl> <chr>    <chr>    <chr>            <chr>      <dbl>    <dbl>
1        2021 AR       A01      AR               P3_S13     694.         3
2        2021 AR       A02      AR               P3_S13      31.1        3
3        2021 AR       A03      AR               P3_S13      58.0        3
4        2021 AR       B        AR               P3_S13     853.         3
5        2021 AR       C10T12   AR               P3_S13      24.0        3
6        2021 AR       C13T15   AR               P3_S13      10.3        3
# ℹ 5 more variables: unit_measure <chr>, unit_mult <dbl>, obs_status <chr>,
#   conf_status <lgl>, last_update <dttm>

Obtaning the Carbon Footprint of the EU

The Carbon Footprint, also known as Consumption-based Responsibility, is an estimate of the greenhouse gases (GHGs) induced by a country’s domestic final demand. This includes: i) direct household emissions (from housing and vehicles), ii) emissions from domestic production (excluding exports), and iii) emissions from foreign economic activities whose production is destined for imports into the country. If we want to calculate the Carbon Footprint for the EU in 2021 by Member State, we can process the data to obtain this figure in different ways. One approach is to first filter by the final consumer country (counterpart_area) and then sum up all industries. It is also recommended to create a specific vector with all the EU countries’ labels to facilitate data processing.

The following chunk of code can be used to carry out this task:

EUcountries <- c("AT", "BE", "BG", "CY", "CZ", "DE", "DK", 
                 "EE", "ES", "FI", "FR", "GR", "HR", "HU", 
                "IE", "IT", "LT", "LU", "LV", "MT", "NL", 
                  "PL", "PT", "RO", "SE", "SI", "SK")

CFms <- CF %>% 
  filter(counterpart_area %in% EUcountries) %>% 
  group_by(time_period, counterpart_area) %>% 
  summarise(obs_value = sum(obs_value)/1000)

head(CFms)
# A tibble: 6 × 3
# Groups:   time_period [1]
  time_period counterpart_area obs_value
        <dbl> <chr>                <dbl>
1        2021 AT                    85.4
2        2021 BE                   113. 
3        2021 BG                    36.6
4        2021 CY                    11.1
5        2021 CZ                    90.0
6        2021 DE                   829. 

The object CFms is an aggregated version of the original results file. It has 27 rows, one for each member state, showing the total emissions in million tons, as we have also divided by 1000. To obtain the total CO2 carbon footprint emissions for that year in the EU, we can simply sum all entries in the file:

sum(CFms$obs_value)
[1] 3468.402

Therefore, the emissions due to EU consumption in 2021 account for approximately 3.5 billion tons CO2e.

Obtaning the Territorial emissions or Production-based emissions of the EU for comparison

The territorial emissions, also known as National inventories or Production-based Responsibility, account for the quantities of greenhouse gases physically emitted within the country. This includes emissions from households (cars and dwellings) and economic activities (fossil energy consumption, industrial processes, and emissions from agriculture). To process the data and obtain such emissions, we can follow a similar approach as above, which involves filtering first and then aggregating the data conveniently. In this case, we need to use the group_by() function to group the emitting country, which corresponds to ref_area in FIGARO terminology.

The following chunk of code can be used to carry out this task:

TEms <- CF %>% 
  filter(ref_area %in% EUcountries) %>% 
  group_by(time_period, ref_area) %>% 
  summarise(obs_value = sum(obs_value)/1000)

head(TEms)
# A tibble: 6 × 3
# Groups:   time_period [1]
  time_period ref_area obs_value
        <dbl> <chr>        <dbl>
1        2021 AT           62.3 
2        2021 BE           95.3 
3        2021 BG           44.2 
4        2021 CY            7.04
5        2021 CZ           88.1 
6        2021 DE          718.  

We can calculate the total for the EU once again and compare it with the carbon footprint by simply adding up the values for all member states:

sum(TEms$obs_value)
[1] 3009.012

The territorial or production-based emissions amount to 3.0 billion tons. Therefore, the EU consumption-based CO2 emissions are higher than production-based, approximately 15% higher in 2021. We can now also determine the net trade balance of CO2 emissions:

sum(CFms$obs_value) - sum(TEms$obs_value)
[1] 459.3903

The net trade balance of CO2 is 459 million tons, meaning that the EU is a net importer of CO2 emissions from the world in 2021.

Analyzing the EU carbon Footprint by source industry

We could process the data to analyze the carbon footprint by industry source. To do this, we need to use the group_by() function to group the data by industry, and then aggregate the results. Additionally, we could add a column to display the shares of each industry in the total EU carbon footprint.

Remember that in FIGARO industries are classified according to the Classification of Economic Activities in the European Community NACE Rev. 2.

The following code can be used to perform this task:

IndCF <- CF %>% 
  filter(counterpart_area %in% EUcountries) %>% 
  group_by(time_period, industry) %>% 
  summarise(obs_value = sum(obs_value)/1000) %>% 
  ungroup() %>% 
  mutate(share = (obs_value/sum(obs_value))*100) %>% 
  arrange(-share)

IndCF
# A tibble: 65 × 4
   time_period industry obs_value share
         <dbl> <chr>        <dbl> <dbl>
 1        2021 D35         1024.  29.5 
 2        2021 HH           697.  20.1 
 3        2021 C24          236.   6.79
 4        2021 C23          211.   6.09
 5        2021 C20          177.   5.10
 6        2021 H49          163.   4.71
 7        2021 C19          113.   3.25
 8        2021 B            102.   2.95
 9        2021 A01           96.7  2.79
10        2021 H51           77.2  2.23
# ℹ 55 more rows

We can observe that 29.5 % of the EU carbon footprint is produced in industry D35 Electricity, gas, steam and air conditioning supply, 20 % directly emitted by households (HH), 6.8 % in C24 Manufacture of basic metals, and 6.1% in C23 Manufacture of other non-metallic mineral products.

Studying the foreign fraction of the EU carbon footprint

Let’s assume that we now want to analyze where in the world, that is outside the EU, CO2 emissions are generated to maintain the EU final consumption. We need to exclude (!) EU emissions and aggregate by non-EU ref_area. We can also add here a column indicating the percentages of emissions for each country.

The following code chunk can be used to accomplish this task:

originCF <- CF %>% 
  filter(counterpart_area %in% EUcountries, 
         !ref_area %in% EUcountries) %>% 
  group_by(time_period, ref_area) %>% 
  summarise(obs_value = sum(obs_value)/1000) %>% 
  ungroup() %>% 
  mutate(share = (obs_value/sum(obs_value))*100) %>% 
  arrange(-share)

originCF
# A tibble: 19 × 4
   time_period ref_area obs_value  share
         <dbl> <chr>        <dbl>  <dbl>
 1        2021 FIGW1       317.   29.8  
 2        2021 CN          293.   27.6  
 3        2021 RU          165.   15.5  
 4        2021 IN           54.8   5.16 
 5        2021 US           54.6   5.14 
 6        2021 TR           28.4   2.67 
 7        2021 GB           19.9   1.87 
 8        2021 KR           19.2   1.81 
 9        2021 ZA           18.4   1.73 
10        2021 JP           16.7   1.57 
11        2021 CA           12.5   1.18 
12        2021 NO           10.8   1.02 
13        2021 SA           10.6   0.994
14        2021 BR           10.3   0.971
15        2021 MX            9.09  0.856
16        2021 ID            8.43  0.794
17        2021 AU            6.65  0.626
18        2021 CH            3.72  0.350
19        2021 AR            3.09  0.291

originCF indicates that the primary producers are the ‘Rest of the world’ (31%), China (27.5%), and Russia (13.5%). We could also analyze which industries in the extra-EU countries are generating the emissions. To estimate these figures, we need to slightly adjust the code and again use the group_by() function to group the data by industry.

The following code chunk can be used to achieve this:

originIndCF <- CF %>% 
  filter(counterpart_area %in% EUcountries, 
         !ref_area %in% EUcountries) %>% 
  group_by(time_period, industry) %>% 
  summarise(obs_value = sum(obs_value)/1000) %>% 
  ungroup() %>% 
  mutate(share = (obs_value/sum(obs_value))*100) %>% 
  arrange(-share)

originIndCF
# A tibble: 64 × 4
   time_period industry obs_value share
         <dbl> <chr>        <dbl> <dbl>
 1        2021 D35          410.  38.6 
 2        2021 C24          140.  13.1 
 3        2021 C20           97.5  9.18
 4        2021 B             84.9  8.00
 5        2021 C23           73.3  6.90
 6        2021 H49           49.0  4.61
 7        2021 H50           42.0  3.95
 8        2021 H51           33.9  3.19
 9        2021 C19           30.0  2.82
10        2021 A01           22.3  2.10
# ℹ 54 more rows

We can observe in object originIndCF that 37.1 % of the foreign fraction of the EU carbon footprint is produced in industry D35 Electricity, gas, steam and air conditioning supply, 13.5 % in C24 Manufacture of basic metals, 9.2% in C20 Manufacture of chemicals and chemical products, and 8.4% in B Mining and quarrying.

EU carbon footprint for all years: looping

Eurostat data is typically published on a yearly basis, but in many cases, we are interested in studying trends over time. In this tutorial, we will show how to loop through the different years and calculate the carbon footprint of the European Union for the period 2010-2021.

Additionally, we will create a basic line chart using the R package plotly to visualize the trend.

The following code chunk will accomplish this task:

processResults <- function(time_period){
  
CF <- readr::read_csv(paste0(pathToData, "CO2_footprints_", time_period,".csv")) %>% 
  filter(counterpart_area %in% EUcountries) %>% 
  group_by(time_period) %>% 
  summarise(obs_value = sum(obs_value)/1000)

return(CF)

}

CFyear <- do.call(rbind, lapply(2010:2021, processResults))
  
CFyear
# A tibble: 12 × 2
   time_period obs_value
         <dbl>     <dbl>
 1        2010     4210.
 2        2011     4129.
 3        2012     3888.
 4        2013     3793.
 5        2014     3672.
 6        2015     3621.
 7        2016     3635.
 8        2017     3659.
 9        2018     3699.
10        2019     3572.
11        2020     3178.
12        2021     3468.
fig <- plot_ly(CFyear, x = ~time_period, y = ~obs_value, type = 'scatter', mode = 'lines+markers') %>% 
  layout( yaxis = list(rangemode = "tozero"))

fig

We can observe a steady decline in the carbon footprint until 2015, followed by an increase from 2016 to 2018, and then a sharp decrease in 2020.

Loading the neccesary data for replicating EUROSTAT estimation

Now, we will calculate the required input-output components needed to replicate the EUROSTAT estimation. Remember that the Leontief model is defined by:

\[x=g'(I-A)^{-1}Y =g'LY\] where \(g'=b'x^-1\), that is, \(g\) refers to industry emissions by unit of output \(x\), \(L\) is the Leontief multipliers matrix, and \(Y\) is the final use by country. First, we upload the FIGARO tables in matrix (matrixFigaro) and long format (dtFigaro) to R:

matrixFigaro <- readr::read_csv(paste0(pathToData, "matrix_eu-ic-io_ind-by-ind_23ed_2021.csv"))

dtFigaro <- readr::read_csv(paste0(
  pathToData,"flatfile_eu-ic-io_ind-by-ind_23ed_2021.csv"))

We could use dim() to know a bit more about our input-output matrix, and unique() to have a clearer idea of the classifications used in the data.

dim(matrixFigaro)
[1] 2950 3175

The matrix has 2950 rows and 3175 columns. Therefore, we can infer that the intermediate and final use matrices are provided together. This is not the case in other input-output databases, for instance, the OECD’s input-output database provides two different files for intermediate and final uses.

We could extract the counterpartArea countries available in the data:

unique(dtFigaro$counterpartArea)
 [1] "AR"    "AT"    "AU"    "BE"    "BG"    "BR"    "CA"    "CH"    "CN"   
[10] "CY"    "CZ"    "DE"    "DK"    "EE"    "ES"    "FI"    "FIGW1" "FR"   
[19] "GB"    "GR"    "HR"    "HU"    "ID"    "IE"    "IN"    "IT"    "JP"   
[28] "KR"    "LT"    "LU"    "LV"    "MT"    "MX"    "NL"    "NO"    "PL"   
[37] "PT"    "RO"    "RU"    "SA"    "SE"    "SI"    "SK"    "TR"    "US"   
[46] "ZA"   

That function shows that 46 countries are included in the data (27 member states + 18 main trade partners + Rest of the World region).

Countries explicitly covered are the 27 EU Member States (Belgium, Bulgaria, Czechia, Denmark, Germany, Estonia, Ireland, Greece, Spain, France, Croatia, Italy, Cyprus, Latvia, Lithuania, Luxembourg, Hungary, Malta, the Netherlands, Austria, Poland, Portugal, Romania, Slovenia, Slovakia, Finland, and Sweden), and 18 main EU trading partners (Argentina, Australia, Brazil, Canada, China, India, Indonesia, Japan, Norway, Mexico, Russia, Saudi Arabia, South Africa, Switzerland, Türkiye, the United Kingdom, and the United States).

unique(dtFigaro$refArea)
 [1] "AR"    "AT"    "AU"    "BE"    "BG"    "BR"    "CA"    "CH"    "CN"   
[10] "CY"    "CZ"    "DE"    "DK"    "EE"    "ES"    "FI"    "FIGW1" "FR"   
[19] "GB"    "GR"    "HR"    "HU"    "ID"    "IE"    "IN"    "IT"    "JP"   
[28] "KR"    "LT"    "LU"    "LV"    "MT"    "MX"    "NL"    "NO"    "PL"   
[37] "PT"    "RO"    "RU"    "SA"    "SE"    "SI"    "SK"    "TR"    "US"   
[46] "ZA"    "W2"   

When doing the same by row, that is, by ref_area we see that there is one extra region called W2. This is a specific code that EUROSTAT uses for value added components, and could be useful if we want to filter for those variables.

unique(dtFigaro$rowIi)
 [1] "A01"     "A02"     "A03"     "B"       "C10T12"  "C13T15"  "C16"    
 [8] "C17"     "C18"     "C19"     "C20"     "C21"     "C22"     "C23"    
[15] "C24"     "C25"     "C26"     "C27"     "C28"     "C29"     "C30"    
[22] "C31_32"  "C33"     "D35"     "E36"     "E37T39"  "F"       "G45"    
[29] "G46"     "G47"     "H49"     "H50"     "H51"     "H52"     "H53"    
[36] "I"       "J58"     "J59_60"  "J61"     "J62_63"  "K64"     "K65"    
[43] "K66"     "L"       "M69_70"  "M71"     "M72"     "M73"     "M74_75" 
[50] "N77"     "N78"     "N79"     "N80T82"  "O84"     "P85"     "Q86"    
[57] "Q87_88"  "R90T92"  "R93"     "S94"     "S95"     "S96"     "T"      
[64] "U"       "D21X31"  "OP_RES"  "OP_NRES" "D1"      "D29X39"  "B2A3G"  

Similarly, we could explore industry classifications by row, that is, using rowIi. We can observe that there are 64 industries, along with the value added components, and other variables.

The FIGARO tables adhere to the statistical concepts and definitions of the European System of Accounts ESA 2010. Accordingly, D21X31 is Taxes less subsidies on products, D1 is Compensation of employees, D29X39 is Other taxes less subsidies on production, B2A3G is Gross operating surplus, OP_RES is Direct purchases abroad by residents, and OP_NRES is Purchases on the domestic territory by non-residents.

unique(dtFigaro$colIi)
 [1] "A01"    "A02"    "A03"    "B"      "C10T12" "C13T15" "C16"    "C17"   
 [9] "C18"    "C19"    "C20"    "C21"    "C22"    "C23"    "C24"    "C25"   
[17] "C26"    "C27"    "C28"    "C29"    "C30"    "C31_32" "C33"    "D35"   
[25] "E36"    "E37T39" "F"      "G45"    "G46"    "G47"    "H49"    "H50"   
[33] "H51"    "H52"    "H53"    "I"      "J58"    "J59_60" "J61"    "J62_63"
[41] "K64"    "K65"    "K66"    "L"      "M69_70" "M71"    "M72"    "M73"   
[49] "M74_75" "N77"    "N78"    "N79"    "N80T82" "O84"    "P85"    "Q86"   
[57] "Q87_88" "R90T92" "R93"    "S94"    "S95"    "S96"    "T"      "U"     
[65] "P3_S13" "P3_S14" "P3_S15" "P51G"   "P5M"   

Finally, by column we identify the 64 NACE FIGARO industries, from A01 agriculture to U international organizations, along with the final use categories.

P3_S13 is Final consumption expenditure by general government. P3_S14 is Final consumption expenditure by households. P3_S15 is Final consumption expenditure by NPISH (Non-Profit Institutions Serving Households), P51G is Gross fixed capital formation, and P5M is Changes in inventories and valuables.

Next, we can calculate \(x\), \(L\) and \(Y\). For making easier the data processing, we need also certain label vectors:

IndustriesFigaro <- c("A01", "A02", "A03", "B", "C10T12","C13T15","C16", 
                       "C17", "C18", "C19", "C20", "C21", "C22", "C23", "C24",
                       "C25", "C26", "C27", "C28", "C29", "C30", "C31_32","C33", 
                       "D35", "E36", "E37T39","F", "G45", "G46", "G47", "H49",
                       "H50", "H51", "H52", "H53", "I", "J58", "J59_60","J61", "J62_63",
                       "K64", "K65", "K66", "L", "M69_70","M71", "M72", "M73", "M74_75",
                       "N77", "N78", "N79", "N80T82","O84", "P85", "Q86", "Q87_88","R90T92",
                       "R93", "S94", "S95", "S96", "T", "U")

FinalDemandFigaro <- c("P3_S13","P3_S14","P3_S15","P51G","P5M")

CountriesFigaro <- unique(dtFigaro$counterpartArea)

LabelsZFigaro <- paste0(rep(CountriesFigaro, each= 64),"_", IndustriesFigaro)

LabelsYFigaro <- paste0(rep(CountriesFigaro, each= 64),"_", FinalDemandFigaro)

LabelsYFigaroEU <- paste0(rep(EUcountries, each= 64),"_", FinalDemandFigaro)

Then we could simply filter, and proceed using R functions for matrix algebra, such as solve(), and diag():

x <- colSums(matrixFigaro %>% 
    select(all_of(LabelsZFigaro)))

Z <- matrixFigaro %>% 
    filter(rowLabels %in% LabelsZFigaro) %>% 
    select(all_of(LabelsZFigaro))

A <- t(t(Z)/x)
A[is.na(A)] <- 0
rownames(A) <- colnames(A)
  
I <- diag(x = 1, nrow=length(x), ncol=length(x))

L <- solve(I-A, sparse=TRUE, tol = 1e-19)
colnames(L) <- rownames(L)

Y <- matrixFigaro %>% 
    filter(rowLabels %in% LabelsZFigaro) %>% 
    select(all_of(LabelsYFigaroEU)) %>% 
  as.matrix()

We need to process also the CO2 emissions data. The CO2 emissions data for the model is extracted from the results file provided by EUROSTAT.

b <- readr::read_csv(paste0(pathToData, "CO2_footprints_2021.csv")) %>% 
    group_by(ref_area, industry) %>% 
    summarise(obsValue = sum(obs_value)) %>% ungroup() %>% 
    mutate(rowLabels = paste0(ref_area, "_", industry)) %>% 
    arrange(match(rowLabels, LabelsZFigaro)) %>% 
    mutate(counterpartArea = ref_area, rowIi = "CO2") %>% 
    rename(colIi = industry,
           refArea = ref_area) %>% 
    select(rowIi, counterpartArea, colIi, obsValue)

g <- b %>%
  filter(colIi!= "HH") %>%
  pull(obsValue)/x

g[!is.finite(g)] <- 0

\(b\) is the territorial emissions by industry, and \(g\) is the territorial emissions per unit of industry output. We now can apply the Leontief model for obtaining the EU carbon footprint:

gLY <- (L*g) %*% Y 

CFrep <- reshape2::melt(gLY, varnames = c("icioiCol", "icioiRow"), 
                        value.name = "obsValue") %>% 
  mutate(refArea = sub("_.*", "", icioiRow),
         counterpartArea = sub("_.*", "", icioiCol),
         rowIi = sub(".*?_", "", icioiRow),
         colIi = sub(".*?_", "", icioiCol)) %>%  
  select(refArea, rowIi, counterpartArea, colIi, obsValue) %>% 
  bind_rows(b %>%
  filter(counterpartArea %in% EUcountries,
    colIi== "HH"))

sum(CFrep$obsValue)/1000
[1] 3465.424

We obtain 3.5 billion tons of CO2 emissions in 2021. The number does not match exactly the one provided above (3465 vs. 3468 thousand tons) because the results file is calculated with the confidential data (‘Reference data’ in FIGARO terminology), while this replication uses the public data (‘Dissemination data’).

However, performing the calculations ourselves allow for new possibilities in the analysis. For instance, we could break down the EU carbon footprint by final industry, that is, by the industry producing the final product consumed in the EU.

For doing that, we need to slightly modify the model expression to: \[x=<g>(I-A)^{-1}y =<g>L<y>\] where \(<g>\) is the diagonalized vector of \(g\), and \(y=Yi\), being \(i\) the summation vector.

gLy <- t(t(L*g)*rowSums(Y))

CFprod <- reshape2::melt(gLy, varnames = c("icioiRow", "icioiCol"), 
                        value.name = "obsValue") %>% 
  mutate(refArea = sub("_.*", "", icioiRow),
         counterpartArea = sub("_.*", "", icioiCol),
         rowIi = sub(".*?_", "", icioiRow),
         colIi = sub(".*?_", "", icioiCol)) %>%  
  select(refArea, rowIi, counterpartArea, colIi, obsValue) %>% 
  bind_rows(b %>%
  filter(counterpartArea %in% EUcountries,
    colIi== "HH"))

CFprodAg <- CFprod %>% 
  group_by(colIi) %>% 
  summarise(obsValue = sum(obsValue)/1000) %>% 
  ungroup() %>% 
  mutate(share = (obsValue/sum(obsValue))*100) %>% 
  arrange(-share)

CFprodAg
# A tibble: 65 × 3
   colIi  obsValue share
   <chr>     <dbl> <dbl>
 1 HH        697.  20.1 
 2 D35       335.   9.67
 3 F         323.   9.32
 4 C10T12    198.   5.73
 5 C29       145.   4.18
 6 C19       104.   3.00
 7 O84       102.   2.93
 8 Q86        96.4  2.78
 9 C28        91.3  2.63
10 G47        87.0  2.51
# ℹ 55 more rows

We can observe that according to CFprodAg 20% of the emissions are produced directly in the households (HH), 9.7% due to the consumption of electricity, 9.3 by the consumption of products form the construction sector (F), 5.7% due to products from C10T12 Manufacture of food products; beverages and tobacco products, 4.2% due to products from C29 Manufacture of motor vehicles, trailers and semi-trailers, 3% due to products from C19 Manufacture of coke and refined petroleum products, and 2.9% due to services provided by O84 Public administration and defence; compulsory social security.

Increasing granularity: FIGARO-E3

The official FIGARO tables have a resolution of 64 industries and products. However, for certain analyses higher granularity might be needed. In that case, we could use the FIGARO-E3 database. This is an inter-country supply and use, input-output experimental statistics for the year 2015.

For avoiding memory issues, we have prepared a ‘toy’ FIGARO-E3 database with only two regions EU and extra-EU, instead of the original 46 countries. The data is stored using R storing format .rds.

First, we upload the input-output components:

Z <- readRDS(paste0(pathToData, "Z.rds"))

Y <- readRDS(paste0(pathToData, "Y.rds"))

X <- readRDS(paste0(pathToData, "X.rds"))

CO2e <- readRDS(paste0(pathToData, "CO2e.rds"))

As we did earlier, we could use function unique() to have a closer look to the data:

unique(X$colPi)
  [1] "A01_A"    "A01_B"    "A01_C"    "A01_D"    "A01_E"    "A01_F"   
  [7] "A01_G"    "A01_H"    "A01_I"    "A01_J"    "A01_K"    "A01_L"   
 [13] "A01_M"    "A01_N"    "A01_O"    "A01_P"    "A01_Q"    "A02"     
 [19] "A03"      "B05"      "B06_A"    "B06_B"    "B06_C"    "B07_A"   
 [25] "B07_B"    "B07_C"    "B07_D"    "B07_E"    "B07_F"    "B07_G"   
 [31] "B07_H"    "B08_A"    "B08_B"    "B08_C"    "C10_A"    "C10_B"   
 [37] "C10_C"    "C10_D"    "C10_E"    "C10_F"    "C10_G"    "C10_H"   
 [43] "C10_I"    "C10_J"    "C11"      "C12"      "C13"      "C14"     
 [49] "C15"      "C16_A"    "C16_B"    "C17_A"    "C17_B"    "C17_C"   
 [55] "C18"      "C191"     "C192"     "C20_A"    "C20_B"    "C20_C"   
 [61] "C20_D"    "C21"      "C22"      "C23_A"    "C23_B"    "C23_C"   
 [67] "C23_D"    "C23_E"    "C23_F"    "C23_G"    "C24_A"    "C24_B"   
 [73] "C24_C"    "C24_D"    "C24_E"    "C24_F"    "C24_G"    "C24_H"   
 [79] "C24_I"    "C24_J"    "C24_K"    "C24_L"    "C24_M"    "C24_N"   
 [85] "C25"      "C26_A"    "C26_B"    "C27"      "C28_A"    "C28_B"   
 [91] "C29"      "C30"      "C31_32"   "C33"      "D3511_A"  "D3511_B" 
 [97] "D3511_C"  "D3511_D"  "D3511_E"  "D3511_F"  "D3511_G"  "D3511_H" 
[103] "D3511_I"  "D3511_J"  "D3511_K"  "D3511_L"  "D3512"    "D3513"   
[109] "D352"     "D353"     "E36"      "E37T39_A" "E37T39_B" "E37T39_C"
[115] "E37T39_D" "E37T39_E" "E37T39_F" "E37T39_G" "E37T39_H" "E37T39_I"
[121] "E37T39_J" "E37T39_K" "E37T39_L" "E37T39_M" "E37T39_N" "E37T39_O"
[127] "E37T39_P" "E37T39_Q" "E37T39_R" "E37T39_S" "E37T39_T" "E37T39_U"
[133] "E37T39_V" "F_A"      "F_B"      "G45"      "G46"      "G47_A"   
[139] "G47_B"    "H49_A"    "H49_B"    "H49_C"    "H50_A"    "H50_B"   
[145] "H51"      "H52"      "H53"      "I"        "J58"      "J59_60"  
[151] "J61"      "J62_63"   "K64"      "K65"      "K66"      "L"       
[157] "M69_70"   "M71"      "M72"      "M73"      "M74_75"   "N77"     
[163] "N78"      "N79"      "N80T82"   "O84"      "P85"      "Q86"     
[169] "Q87_88"   "R90T92"   "R93"      "S94"      "S95"      "S96"     
[175] "T"        "U"       

We can observe that this new database has 176 industries, that is, 112 industries more than the official FIGARO tables.

We could obtain the components of the Leontief model, as we did above:

x <- X$obsValue
A <- t(t(Z)/x)
A[is.na(A)] <- 0
rownames(A) <- colnames(A)

I <- diag(x = 1, nrow=length(x), ncol=length(x))

L <- solve(I-A, sparse=TRUE, tol = 1e-19)
colnames(L) <- rownames(L)

Now, we could have a closer look to the emission data. In this case, we have data not only on CO2 but also on CH4, N2O and F-gases.

co2 <- CO2e %>% 
  filter(!colIi %in% FinalDemandFigaro) %>% 
  select(codeIndicator, icioCol, obsValue) %>% 
  tidyr::pivot_wider(names_from = icioCol, values_from = obsValue) %>% 
  select(-codeIndicator) %>% 
  as.matrix()

rownames(co2) <- c("CH4", "CO2", "F-gases", "N2O")

g <- t(t(co2)/x)
g[!is.finite(g)] <- 0

co2hh <- CO2e %>% 
  filter(colIi %in% FinalDemandFigaro,
         counterpartArea == "EU") %>% 
  group_by(codeIndicator) %>% 
  summarise(obsValue = sum(obsValue)/1000) %>% 
  ungroup() %>% 
  mutate(refArea = "EU",
         counterpartArea = "EU",
         colIi = "HH",
         rowIi = "HH")

Once we have all necessary components, we could make the calculations:

gLY <- (g %*% L %*% Y [,1:5])/1000

sum(gLY)+ sum(co2hh$obsValue)
[1] 4655.177
gLY
         EU_P3_S13  EU_P3_S14  EU_P3_S15   EU_P51G    EU_P5M
CH4      51.263793  363.02359  2.3959318  73.19746  8.425231
CO2     348.287524 1891.08470 23.2165226 780.46521 48.209174
F-gases   9.420244   58.52193  0.5461871  14.24468  1.704623
N2O      16.902380  171.97599  0.9262335  19.15557  6.136410

In this case, we have applied the following equation: \[x=G(I-A)^{-1}Y =GLY\] where \(G\) is the emission intensity matrix, describing in each row the emission intensities by industry for the greenhouse gases.

As we did above, we could obtain the carbon footprint and shares of each final product:

gLy <- t(t(L*colSums(g))*rowSums(Y [,1:5]))

CFprod <- reshape2::melt(gLy, varnames = c("icioiRow", "icioiCol"),
                        value.name = "obsValue") %>%
  mutate(refArea = sub("_.*", "", icioiRow),
         counterpartArea = sub("_.*", "", icioiCol),
         rowIi = sub(".*?_", "", icioiRow),
         colIi = sub(".*?_", "", icioiCol)) %>%
  select(refArea, rowIi, counterpartArea, colIi, obsValue)

CFprodAg <- CFprod %>%
  group_by(colIi) %>%
  summarise(obsValue = sum(obsValue)/1000) %>%
  ungroup() %>%
  mutate(share = (obsValue/sum(obsValue))*100) %>%
  arrange(-share)

CFprodAg
# A tibble: 176 × 3
   colIi   obsValue share
   <chr>      <dbl> <dbl>
 1 F_A        360.   9.25
 2 C10_J      168.   4.32
 3 D3511_A    154.   3.96
 4 C192       151.   3.88
 5 I          147.   3.78
 6 D353       141.   3.63
 7 O84        137.   3.52
 8 C29        127.   3.27
 9 Q86        116.   2.98
10 G47_A       96.6  2.48
# ℹ 166 more rows

We can observe that 9.25 % are emissions due to F_A Contruction, 4.3% C10_J Processing of Food products nec, 4% D3511_A Production of electricity by coal, C192 Petroleum Refinery, 3.8% I Hotels and restaurants, 3.6% D353 Steam and hot water supply, etc.

If you wish to focus on the subcomponents of D35, you could use the function substr()

CFd35 <- CFprod %>%
  filter(substr(colIi, 1,3) == "D35") %>% 
  group_by(colIi) %>%
  summarise(obsValue = sum(obsValue)/1000) %>%
  ungroup() %>%
  mutate(share = (obsValue/sum(obsValue))*100) %>%
  arrange(-share)

CFd35
# A tibble: 16 × 3
   colIi   obsValue   share
   <chr>      <dbl>   <dbl>
 1 D3511_A  1.54e+2 3.13e+1
 2 D353     1.41e+2 2.87e+1
 3 D3513    9.18e+1 1.86e+1
 4 D3511_B  3.22e+1 6.54e+0
 5 D352     2.86e+1 5.80e+0
 6 D3512    1.83e+1 3.72e+0
 7 D3511_F  1.07e+1 2.17e+0
 8 D3511_C  7.16e+0 1.45e+0
 9 D3511_G  4.22e+0 8.58e-1
10 D3511_D  2.73e+0 5.54e-1
11 D3511_E  8.49e-1 1.72e-1
12 D3511_L  4.37e-1 8.88e-2
13 D3511_J  1.11e-1 2.25e-2
14 D3511_K  7.81e-2 1.58e-2
15 D3511_H  6.38e-3 1.30e-3
16 D3511_I  2.92e-8 5.92e-9

In this way, we could have a better understanding on how the emissions in D35 occur: 31% D3511_A Production of electricity by coal, 28.7% in D353 Steam and hot water supply, 18.6% in D3513 Distribution and trade of electricity, 6.5% in D3511_B Production of electricity by gas, 5.8% in D352 Manufacture of gas, distribution of gaseous fuels through mains, etc.