pathToData <- "/home/NET1/b7_training/"
CF <- readr::read_csv(paste0(pathToData, "CO2_footprints_2021.csv"))Eurostat’s carbon footprint with FIGARO in R
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.
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"))
figWe 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.