diff --git a/analysis/preprocessing/full_code_europe.inequality.Rmd b/analysis/preprocessing/full_code_europe.inequality.Rmd
index f427a5fbabb233926db36165195f0f281b529207..e2120d87668ce2703252f9252296674a2b3bfd89 100644
--- a/analysis/preprocessing/full_code_europe.inequality.Rmd
+++ b/analysis/preprocessing/full_code_europe.inequality.Rmd
@@ -43,15 +43,31 @@ pacman::p_load(tidyverse,
 
 ```
 
+# EXIOBASE
 
-# Exiobase
+EE-MRIO: EXIOBASE version3 ixi, from: https://zenodo.org/record/3583071#.XjC7kSN4wpY [accessed on 12.03.2020] @stadler_exiobase_2018
+
+We use standard input-output calculations to calculate total intensity vectors in EXIOBASE. 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. 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 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$): 
+
+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$): 
+
+The total intensity vectors ($TIV$) are calculated by pre-multiplying the direct intensity vectors ($DIV$) by the Leontief inverse ($L$):
+
+The footprint is then calculated by row-wise multiplying the TIV by final demand:
+
+Before that final step, however, we decompose national household final demand by income quintile according to the structure of the household budget survey, 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$):
 
 ```{r exiobase, eval = FALSE}
 
 # EXIOBASE_cluster_ixi_version
 
 # data directories (on cluster)
-data_dir_exiobase = paste("/",file.path("data","metab","Exiobase", fsep=.Platform$file.sep),sep="")
+#data_dir_exiobase = paste("/",file.path("data","metab","Exiobase", fsep=.Platform$file.sep),sep="")
+data_dir_exiobase = here("analysis", "preprocessing", "EXIOBASE")
 
 years_exiobase_ixi = c(2005,2010,2015)
 
@@ -125,11 +141,9 @@ write.csv(satellite_FD, paste0(data_dir_exiobase, "/IOT_", year_current, "_ixi/s
 
 # CO2 combustion air
 CO2_combustion_air = satellite[24,]
-#write.csv(CO2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/co2_combustion_air_",year_current,"_ixi.csv"))
 DIV_co2_combustion_air = CO2_combustion_air/total_output
 DIV_co2_combustion_air[is.na(DIV_co2_combustion_air)]=0
 DIV_co2_combustion_air[DIV_co2_combustion_air == Inf]<-0
-#write.csv(DIV_co2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_co2_combustion_air_",year_current,"_ixi.csv"))
 TIV_co2_combustion_air = as.vector(DIV_co2_combustion_air) %*% L
 
 write.csv(TIV_co2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_co2_combustion_air_",year_current,"_ixi.csv"))
@@ -240,13 +254,10 @@ write.csv(TIV_country_breakdown_co2_waste_fossil_air_w_labels, paste0(data_dir_e
 
 # CH4 combustion air
 CH4_combustion_air = satellite[25,]
-#write.csv(CH4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/ch4_combustion_air_",year_current,"_ixi.csv"))
 CH4_combustion_air = CH4_combustion_air*28
-#write.csv(CH4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/ch4_CO2eq_combustion_air_",year_current,"_ixi.csv"))
 DIV_ch4_combustion_air = CH4_combustion_air/total_output
 DIV_ch4_combustion_air[is.na(DIV_ch4_combustion_air)]=0
 DIV_ch4_combustion_air[DIV_ch4_combustion_air == Inf]<-0
-#write.csv(DIV_ch4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_ch4_CO2eq_combustion_air_",year_current,"_ixi.csv"))
 TIV_ch4_combustion_air = as.vector(DIV_ch4_combustion_air) %*% L
 
 write.csv(TIV_ch4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_ch4_CO2eq_combustion_air_",year_current,"_ixi.csv"))
@@ -457,13 +468,10 @@ write.csv(TIV_country_breakdown_ch4_waste_air_w_labels, paste0(data_dir_exiobase
 
 # N2O combustion air
 N2O_combustion_air = satellite[26,]
-#write.csv(N2O_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/n2o_combustion_air_",year_current,"_ixi.csv"))
 N2O_combustion_air = N2O_combustion_air*265
-#write.csv(N2O_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/n2o_CO2eq_combustion_air_",year_current,"_ixi.csv"))
 DIV_n2o_combustion_air = N2O_combustion_air/total_output
 DIV_n2o_combustion_air[is.na(DIV_n2o_combustion_air)]=0
 DIV_n2o_combustion_air[DIV_n2o_combustion_air == Inf]<-0
-#write.csv(DIV_n2o_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_n2o_CO2eq_combustion_air_",year_current,"_ixi.csv"))
 TIV_n2o_combustion_air = as.vector(DIV_n2o_combustion_air) %*% L
 
 write.csv(TIV_n2o_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_n2o_CO2eq_combustion_air_",year_current,"_ixi.csv"))
@@ -500,13 +508,10 @@ write.csv(TIV_country_breakdown_n2o_agriculture_air_w_labels, paste0(data_dir_ex
 
 # SF6 air
 SF6_air = satellite[424,]
-#write.csv(SF6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/sf6_air_",year_current,"_ixi.csv"))
 SF6_air = SF6_air*23500
-#write.csv(SF6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/sf6_CO2eq_air_",year_current,"_ixi.csv"))
 DIV_sf6_air = SF6_air/total_output
 DIV_sf6_air[is.na(DIV_sf6_air)]=0
 DIV_sf6_air[DIV_sf6_air == Inf]<-0
-#write.csv(DIV_sf6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_sf6_CO2eq_air_",year_current,"_ixi.csv"))
 TIV_sf6_air = as.vector(DIV_sf6_air) %*% L
 
 write.csv(TIV_sf6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_sf6_CO2eq_air_",year_current,"_ixi.csv"))
@@ -523,11 +528,9 @@ write.csv(TIV_country_breakdown_sf6_air_w_labels, paste0(data_dir_exiobase, "/IO
 
 # HFC air
 HFC_air = satellite[425,]
-#write.csv(HFC_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/hfc_CO2eq_air_",year_current,"_ixi.csv"))
 DIV_hfc_air = HFC_air/total_output
 DIV_hfc_air[is.na(DIV_hfc_air)]=0
 DIV_hfc_air[DIV_hfc_air == Inf]<-0
-#write.csv(DIV_hfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_hfc_CO2eq_air_",year_current,"_ixi.csv"))
 TIV_hfc_air = as.vector(DIV_hfc_air) %*% L
 
 write.csv(TIV_hfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_hfc_CO2eq_air_",year_current,"_ixi.csv"))
@@ -544,11 +547,9 @@ write.csv(TIV_country_breakdown_hfc_air_w_labels, paste0(data_dir_exiobase, "/IO
 
 # PFC air
 PFC_air = satellite[426,]
-#write.csv(PFC_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/satellite/pfc_CO2eq_air_",year_current,"_ixi.csv"))
 DIV_pfc_air = PFC_air/total_output
 DIV_pfc_air[is.na(DIV_pfc_air)]=0
 DIV_pfc_air[DIV_pfc_air == Inf]<-0
-#write.csv(DIV_pfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/DIV_pfc_CO2eq_air_",year_current,"_ixi.csv"))
 TIV_pfc_air = as.vector(DIV_pfc_air) %*% L
 
 write.csv(TIV_pfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_ixi/TIV_pfc_CO2eq_air_",year_current,"_ixi.csv"))
@@ -658,11 +659,9 @@ write.csv(satellite_FD, paste0(data_dir_exiobase, "/IOT_", year_current, "_pxp/s
 
 # CO2 combustion air
 CO2_combustion_air = satellite[24,]
-#write.csv(CO2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_combustion_air_",year_current,"_pxp.csv"))
 DIV_co2_combustion_air = CO2_combustion_air/total_output
 DIV_co2_combustion_air[is.na(DIV_co2_combustion_air)]=0
 DIV_co2_combustion_air[DIV_co2_combustion_air == Inf]<-0
-#write.csv(DIV_co2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_combustion_air_",year_current,"_pxp.csv"))
 TIV_co2_combustion_air = as.vector(DIV_co2_combustion_air) %*% L
 
 write.csv(TIV_co2_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_combustion_air_",year_current,"_pxp.csv"))
@@ -680,11 +679,9 @@ write.csv(TIV_country_breakdown_co2_combustion_air_w_labels, paste0(data_dir_exi
 # CO2 non-combustion air
 ## cement
 CO2_noncombustion_cement_air = satellite[93,]
-#saveRDS(CO2_noncombustion_cement_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_noncombustion_cement_air_",year_current,".rds"))
 DIV_co2_noncombustion_cement_air = CO2_noncombustion_cement_air/total_output
 DIV_co2_noncombustion_cement_air[is.na(DIV_co2_noncombustion_cement_air)]=0
 DIV_co2_noncombustion_cement_air[DIV_co2_noncombustion_cement_air == Inf]<-0
-#saveRDS(DIV_co2_noncombustion_cement_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_noncombustion_cement_air_",year_current,".rds"))
 TIV_co2_noncombustion_cement_air = as.vector(DIV_co2_noncombustion_cement_air) %*% L
 
 write.csv(TIV_co2_noncombustion_cement_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_noncombustion_cement_air_",year_current,"_pxp.csv"))
@@ -700,11 +697,9 @@ write.csv(TIV_country_breakdown_co2_noncombustion_cement_air_w_labels, paste0(da
 
 ## lime
 CO2_noncombustion_lime_air = satellite[94,]
-#saveRDS(CO2_noncombustion_lime_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_noncombustion_lime_air_",year_current,".rds"))
 DIV_co2_noncombustion_lime_air = CO2_noncombustion_lime_air/total_output
 DIV_co2_noncombustion_lime_air[is.na(DIV_co2_noncombustion_lime_air)]=0
 DIV_co2_noncombustion_lime_air[DIV_co2_noncombustion_lime_air == Inf]<-0
-#saveRDS(DIV_co2_noncombustion_lime_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_noncombustion_lime_air_",year_current,".rds"))
 TIV_co2_noncombustion_lime_air = as.vector(DIV_co2_noncombustion_lime_air) %*% L
 
 write.csv(TIV_co2_noncombustion_lime_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_noncombustion_lime_air_",year_current,"_pxp.csv"))
@@ -721,11 +716,9 @@ write.csv(TIV_country_breakdown_co2_noncombustion_lime_air_w_labels, paste0(data
 
 # CO2 agriculture peat decay air
 CO2_agriculture_peatdecay_air = satellite[428,]
-#saveRDS(CO2_agriculture_peatdecay_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_agriculture_peatdecay_air_",year_current,".rds"))
 DIV_co2_agriculture_peatdecay_air = CO2_agriculture_peatdecay_air/total_output
 DIV_co2_agriculture_peatdecay_air[is.na(DIV_co2_agriculture_peatdecay_air)]=0
 DIV_co2_agriculture_peatdecay_air[DIV_co2_agriculture_peatdecay_air == Inf]<-0
-#saveRDS(DIV_co2_agriculture_peatdecay_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_agriculture_peatdecay_air_",year_current,".rds"))
 TIV_co2_agriculture_peatdecay_air = as.vector(DIV_co2_agriculture_peatdecay_air) %*% L
 
 write.csv(TIV_co2_agriculture_peatdecay_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_agriculture_peatdecay_air_",year_current,"_pxp.csv"))
@@ -743,11 +736,9 @@ write.csv(TIV_country_breakdown_co2_agriculture_peatdecay_air_w_labels, paste0(d
 # CO2 waste air
 ## biogenic
 CO2_waste_biogenic_air = satellite[438,]
-#saveRDS(CO2_waste_biogenic_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_waste_biogenic_air_",year_current,".rds"))
 DIV_co2_waste_biogenic_air = CO2_waste_biogenic_air/total_output
 DIV_co2_waste_biogenic_air[is.na(DIV_co2_waste_biogenic_air)]=0
 DIV_co2_waste_biogenic_air[DIV_co2_waste_biogenic_air == Inf]<-0
-#saveRDS(DIV_co2_waste_biogenic_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_waste_biogenic_air_",year_current,".rds"))
 TIV_co2_biogenic_air = as.vector(DIV_co2_waste_biogenic_air) %*% L
 
 write.csv(TIV_co2_biogenic_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_biogenic_air_",year_current,"_pxp.csv"))
@@ -763,11 +754,9 @@ write.csv(TIV_country_breakdown_co2_biogenic_air_w_labels, paste0(data_dir_exiob
 
 ## fossil
 CO2_waste_fossil_air = satellite[439,]
-#saveRDS(CO2_waste_fossil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/co2_waste_fossil_air_",year_current,".rds"))
 DIV_co2_waste_fossil_air = CO2_waste_fossil_air/total_output
 DIV_co2_waste_fossil_air[is.na(DIV_co2_waste_fossil_air)]=0
 DIV_co2_waste_fossil_air[DIV_co2_waste_fossil_air == Inf]<-0
-#saveRDS(DIV_co2_waste_fossil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_co2_waste_fossil_air_",year_current,".rds"))
 TIV_co2_waste_fossil_air = as.vector(DIV_co2_waste_fossil_air) %*% L
 
 write.csv(TIV_co2_waste_fossil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_co2_waste_fossil_air_",year_current,"_pxp.csv"))
@@ -784,16 +773,12 @@ write.csv(TIV_country_breakdown_co2_waste_fossil_air_w_labels, paste0(data_dir_e
 
 # CH4 combustion air
 CH4_combustion_air = satellite[25,]
-#saveRDS(CH4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_combustion_air_",year_current,".rds"))
 CH4_combustion_air = CH4_combustion_air*28
-#saveRDS(CH4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_combustion_air_",year_current,".rds"))
 DIV_ch4_combustion_air = CH4_combustion_air/total_output
 DIV_ch4_combustion_air[is.na(DIV_ch4_combustion_air)]=0
 DIV_ch4_combustion_air[DIV_ch4_combustion_air == Inf]<-0
-#saveRDS(DIV_ch4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_combustion_air_",year_current,".rds"))
 TIV_ch4_combustion_air = as.vector(DIV_ch4_combustion_air) %*% L
 
-#saveRDS(TIV_ch4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_combustion_air_",year_current,".rds"))
 write.csv(TIV_ch4_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_combustion_air_",year_current,"_pxp.csv"))
 
 TIV_breakdown_ch4_combustion_air = as.vector(DIV_ch4_combustion_air) * L
@@ -809,13 +794,10 @@ write.csv(TIV_country_breakdown_ch4_combustion_air_w_labels, paste0(data_dir_exi
 # CH4 non-combustion air
 ## gas
 CH4_noncombustion_gas_air = satellite[68,]
-#saveRDS(CH4_noncombustion_gas_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_gas_air_",year_current,".rds"))
 CH4_noncombustion_gas_air = CH4_noncombustion_gas_air*28
-#saveRDS(CH4_noncombustion_gas_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_gas_air_",year_current,".rds"))
 DIV_ch4_noncombustion_gas_air = CH4_noncombustion_gas_air/total_output
 DIV_ch4_noncombustion_gas_air[is.na(DIV_ch4_noncombustion_gas_air)]=0
 DIV_ch4_noncombustion_gas_air[DIV_ch4_noncombustion_gas_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_gas_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_gas_air_",year_current,".rds"))
 TIV_ch4_noncombustion_gas_air = as.vector(DIV_ch4_noncombustion_gas_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_gas_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_gas_air_",year_current,"_pxp.csv"))
@@ -831,13 +813,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_gas_air_w_labels, paste0(data_
 
 ## oil
 CH4_noncombustion_oil_air = satellite[69,]
-#saveRDS(CH4_noncombustion_oil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_oil_air_",year_current,".rds"))
 CH4_noncombustion_oil_air = CH4_noncombustion_oil_air*28
-#saveRDS(CH4_noncombustion_oil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_oil_air_",year_current,".rds"))
 DIV_ch4_noncombustion_oil_air = CH4_noncombustion_oil_air/total_output
 DIV_ch4_noncombustion_oil_air[is.na(DIV_ch4_noncombustion_oil_air)]=0
 DIV_ch4_noncombustion_oil_air[DIV_ch4_noncombustion_oil_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_oil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_oil_air_",year_current,".rds"))
 TIV_ch4_noncombustion_oil_air = as.vector(DIV_ch4_noncombustion_oil_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_oil_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_oil_air_",year_current,"_pxp.csv"))
@@ -853,13 +832,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_oil_air_w_labels, paste0(data_
 
 ## anthracite
 CH4_noncombustion_anthracite_air = satellite[70,]
-#saveRDS(CH4_noncombustion_anthracite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_anthracite_air_",year_current,".rds"))
 CH4_noncombustion_anthracite_air = CH4_noncombustion_anthracite_air*28
-#saveRDS(CH4_noncombustion_anthracite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_anthracite_air_",year_current,".rds"))
 DIV_ch4_noncombustion_anthracite_air = CH4_noncombustion_anthracite_air/total_output
 DIV_ch4_noncombustion_anthracite_air[is.na(DIV_ch4_noncombustion_anthracite_air)]=0
 DIV_ch4_noncombustion_anthracite_air[DIV_ch4_noncombustion_anthracite_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_anthracite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_anthracite_air_",year_current,".rds"))
 TIV_ch4_noncombustion_anthracite_air = as.vector(DIV_ch4_noncombustion_anthracite_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_anthracite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_anthracite_air_",year_current,"_pxp.csv"))
@@ -875,13 +851,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_anthracite_air_w_labels, paste
 
 ## bituminous coal
 CH4_noncombustion_bituminouscoal_air = satellite[71,]
-#saveRDS(CH4_noncombustion_bituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_bituminouscoal_air_",year_current,".rds"))
 CH4_noncombustion_bituminouscoal_air = CH4_noncombustion_bituminouscoal_air*28
-#saveRDS(CH4_noncombustion_bituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_bituminouscoal_air_",year_current,".rds"))
 DIV_ch4_noncombustion_bituminouscoal_air = CH4_noncombustion_bituminouscoal_air/total_output
 DIV_ch4_noncombustion_bituminouscoal_air[is.na(DIV_ch4_noncombustion_bituminouscoal_air)]=0
 DIV_ch4_noncombustion_bituminouscoal_air[DIV_ch4_noncombustion_bituminouscoal_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_bituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_bituminouscoal_air_",year_current,".rds"))
 TIV_ch4_noncombustion_bituminouscoal_air = as.vector(DIV_ch4_noncombustion_bituminouscoal_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_bituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_bituminouscoal_air_",year_current,"_pxp.csv"))
@@ -897,13 +870,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_bituminouscoal_air_w_labels, p
 
 ## coking coal
 CH4_noncombustion_cokingcoal_air = satellite[72,]
-#saveRDS(CH4_noncombustion_cokingcoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_cokingcoal_air_",year_current,".rds"))
 CH4_noncombustion_cokingcoal_air = CH4_noncombustion_cokingcoal_air*28
-#saveRDS(CH4_noncombustion_cokingcoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_cokingcoal_air_",year_current,".rds"))
 DIV_ch4_noncombustion_cokingcoal_air = CH4_noncombustion_cokingcoal_air/total_output
 DIV_ch4_noncombustion_cokingcoal_air[is.na(DIV_ch4_noncombustion_cokingcoal_air)]=0
 DIV_ch4_noncombustion_cokingcoal_air[DIV_ch4_noncombustion_cokingcoal_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_cokingcoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_cokingcoal_air_",year_current,".rds"))
 TIV_ch4_noncombustion_cokingcoal_air = as.vector(DIV_ch4_noncombustion_cokingcoal_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_cokingcoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_cokingcoal_air_",year_current,"_pxp.csv"))
@@ -919,13 +889,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_cokingcoal_air_w_labels, paste
 
 ## lignite
 CH4_noncombustion_lignite_air = satellite[73,]
-#saveRDS(CH4_noncombustion_lignite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_lignite_air_",year_current,".rds"))
 CH4_noncombustion_lignite_air = CH4_noncombustion_lignite_air*28
-#saveRDS(CH4_noncombustion_lignite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_lignite_air_",year_current,".rds"))
 DIV_ch4_noncombustion_lignite_air = CH4_noncombustion_lignite_air/total_output
 DIV_ch4_noncombustion_lignite_air[is.na(DIV_ch4_noncombustion_lignite_air)]=0
 DIV_ch4_noncombustion_lignite_air[DIV_ch4_noncombustion_lignite_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_lignite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_lignite_air_",year_current,".rds"))
 TIV_ch4_noncombustion_lignite_air = as.vector(DIV_ch4_noncombustion_lignite_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_lignite_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_lignite_air_",year_current,"_pxp.csv"))
@@ -941,13 +908,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_lignite_air_w_labels, paste0(d
 
 ## subbituminous coal
 CH4_noncombustion_subbituminouscoal_air = satellite[74,]
-#saveRDS(CH4_noncombustion_subbituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_subbituminouscoal_air_",year_current,".rds"))
 CH4_noncombustion_subbituminouscoal_air = CH4_noncombustion_subbituminouscoal_air*28
-#saveRDS(CH4_noncombustion_subbituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_subbituminouscoal_air_",year_current,".rds"))
 DIV_ch4_noncombustion_subbituminouscoal_air = CH4_noncombustion_subbituminouscoal_air/total_output
 DIV_ch4_noncombustion_subbituminouscoal_air[is.na(DIV_ch4_noncombustion_subbituminouscoal_air)]=0
 DIV_ch4_noncombustion_subbituminouscoal_air[DIV_ch4_noncombustion_subbituminouscoal_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_subbituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_subbituminouscoal_air_",year_current,".rds"))
 TIV_ch4_noncombustion_subbituminouscoal_air = as.vector(DIV_ch4_noncombustion_subbituminouscoal_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_subbituminouscoal_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_subbituminouscoal_air_",year_current,"_pxp.csv"))
@@ -963,13 +927,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_subbituminouscoal_air_w_labels
 
 ## oil refinery
 CH4_noncombustion_oilrefinery_air = satellite[75,]
-#saveRDS(CH4_noncombustion_oilrefinery_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_noncombustion_oilrefinery_air_",year_current,".rds"))
 CH4_noncombustion_oilrefinery_air = CH4_noncombustion_oilrefinery_air*28
-#saveRDS(CH4_noncombustion_oilrefinery_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_noncombustion_oilrefinery_air_",year_current,".rds"))
 DIV_ch4_noncombustion_oilrefinery_air = CH4_noncombustion_oilrefinery_air/total_output
 DIV_ch4_noncombustion_oilrefinery_air[is.na(DIV_ch4_noncombustion_oilrefinery_air)]=0
 DIV_ch4_noncombustion_oilrefinery_air[DIV_ch4_noncombustion_oilrefinery_air == Inf]<-0
-#saveRDS(DIV_ch4_noncombustion_oilrefinery_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_noncombustion_oilrefinery_air_",year_current,".rds"))
 TIV_ch4_noncombustion_oilrefinery_air = as.vector(DIV_ch4_noncombustion_oilrefinery_air) %*% L
 
 write.csv(TIV_ch4_noncombustion_oilrefinery_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_noncombustion_oilrefinery_air_",year_current,"_pxp.csv"))
@@ -986,13 +947,10 @@ write.csv(TIV_country_breakdown_ch4_noncombustion_oilrefinery_air_w_labels, past
 
 # CH4 agriculture air
 CH4_agriculture_air = satellite[427,]
-#saveRDS(CH4_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_agriculture_air_",year_current,".rds"))
 CH4_agriculture_air = CH4_agriculture_air*28
-#saveRDS(CH4_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_agriculture_air_",year_current,".rds"))
 DIV_ch4_agriculture_air = CH4_agriculture_air/total_output
 DIV_ch4_agriculture_air[is.na(DIV_ch4_agriculture_air)]=0
 DIV_ch4_agriculture_air[DIV_ch4_agriculture_air == Inf]<-0
-#saveRDS(DIV_ch4_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_agriculture_air_",year_current,".rds"))
 TIV_ch4_agriculture_air = as.vector(DIV_ch4_agriculture_air) %*% L
 
 write.csv(TIV_ch4_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_agriculture_air_",year_current,"_pxp.csv"))
@@ -1009,13 +967,10 @@ write.csv(TIV_country_breakdown_ch4_agriculture_air_w_labels, paste0(data_dir_ex
 
 # CH4 waste air
 CH4_waste_air = satellite[436,]
-#saveRDS(CH4_waste_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_waste_air_",year_current,".rds"))
 CH4_waste_air = CH4_waste_air*28
-#saveRDS(CH4_waste_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/ch4_CO2eq_waste_air_",year_current,".rds"))
 DIV_ch4_waste_air = CH4_waste_air/total_output
 DIV_ch4_waste_air[is.na(DIV_ch4_waste_air)]=0
 DIV_ch4_waste_air[DIV_ch4_waste_air == Inf]<-0
-#saveRDS(DIV_ch4_waste_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_ch4_CO2eq_waste_air_",year_current,".rds"))
 TIV_ch4_waste_air = as.vector(DIV_ch4_waste_air) %*% L
 
 write.csv(TIV_ch4_waste_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_ch4_CO2eq_waste_air_",year_current,"_pxp.csv"))
@@ -1032,16 +987,12 @@ write.csv(TIV_country_breakdown_ch4_waste_air_w_labels, paste0(data_dir_exiobase
 
 # N2O combustion air
 N2O_combustion_air = satellite[26,]
-#saveRDS(N2O_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/n2o_combustion_air_",year_current,".rds"))
 N2O_combustion_air = N2O_combustion_air*265
-#saveRDS(N2O_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/n2o_CO2eq_combustion_air_",year_current,".rds"))
 DIV_n2o_combustion_air = N2O_combustion_air/total_output
 DIV_n2o_combustion_air[is.na(DIV_n2o_combustion_air)]=0
 DIV_n2o_combustion_air[DIV_n2o_combustion_air == Inf]<-0
-#saveRDS(DIV_n2o_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_n2o_CO2eq_combustion_air_",year_current,".rds"))
 TIV_n2o_combustion_air = as.vector(DIV_n2o_combustion_air) %*% L
 
-#saveRDS(TIV_n2o_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_n2o_CO2eq_combustion_air_",year_current,".rds"))
 write.csv(TIV_n2o_combustion_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_n2o_CO2eq_combustion_air_",year_current,"_pxp.csv"))
 
 TIV_breakdown_n2o_combustion_air = as.vector(DIV_n2o_combustion_air) * L
@@ -1056,13 +1007,10 @@ write.csv(TIV_country_breakdown_n2o_combustion_air_w_labels, paste0(data_dir_exi
 
 # N2O agriculture air
 N2O_agriculture_air = satellite[430,]
-#saveRDS(N2O_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/n2o_agriculture_air_",year_current,".rds"))
 N2O_agriculture_air = N2O_agriculture_air*265
-#saveRDS(N2O_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/n2o_CO2eq_agriculture_air_",year_current,".rds"))
 DIV_n2o_agriculture_air = N2O_agriculture_air/total_output
 DIV_n2o_agriculture_air[is.na(DIV_n2o_agriculture_air)]=0
 DIV_n2o_agriculture_air[DIV_n2o_agriculture_air == Inf]<-0
-#saveRDS(DIV_n2o_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_n2o_CO2eq_agriculture_air_",year_current,".rds"))
 TIV_n2o_agriculture_air = as.vector(DIV_n2o_agriculture_air) %*% L
 
 write.csv(TIV_n2o_agriculture_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_n2o_CO2eq_agriculture_air_",year_current,"_pxp.csv"))
@@ -1079,13 +1027,10 @@ write.csv(TIV_country_breakdown_n2o_agriculture_air_w_labels, paste0(data_dir_ex
 
 # SF6 air
 SF6_air = satellite[424,]
-#write.csv(SF6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/sf6_air_",year_current,"_pxp.csv"))
 SF6_air = SF6_air*23500
-#write.csv(SF6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/sf6_CO2eq_air_",year_current,"_pxp.csv"))
 DIV_sf6_air = SF6_air/total_output
 DIV_sf6_air[is.na(DIV_sf6_air)]=0
 DIV_sf6_air[DIV_sf6_air == Inf]<-0
-#write.csv(DIV_sf6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_sf6_CO2eq_air_",year_current,"_pxp.csv"))
 TIV_sf6_air = as.vector(DIV_sf6_air) %*% L
 
 write.csv(TIV_sf6_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_sf6_CO2eq_air_",year_current,"_pxp.csv"))
@@ -1102,11 +1047,9 @@ write.csv(TIV_country_breakdown_sf6_air_w_labels, paste0(data_dir_exiobase, "/IO
 
 # HFC air
 HFC_air = satellite[425,]
-#write.csv(HFC_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/hfc_CO2eq_air_",year_current,"_pxp.csv"))
 DIV_hfc_air = HFC_air/total_output
 DIV_hfc_air[is.na(DIV_hfc_air)]=0
 DIV_hfc_air[DIV_hfc_air == Inf]<-0
-#write.csv(DIV_hfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_hfc_CO2eq_air_",year_current,"_pxp.csv"))
 TIV_hfc_air = as.vector(DIV_hfc_air) %*% L
 
 write.csv(TIV_hfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_hfc_CO2eq_air_",year_current,"_pxp.csv"))
@@ -1123,11 +1066,9 @@ write.csv(TIV_country_breakdown_hfc_air_w_labels, paste0(data_dir_exiobase, "/IO
 
 # PFC air
 PFC_air = satellite[426,]
-#write.csv(PFC_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/satellite/pfc_CO2eq_air_",year_current,"_pxp.csv"))
 DIV_pfc_air = PFC_air/total_output
 DIV_pfc_air[is.na(DIV_pfc_air)]=0
 DIV_pfc_air[DIV_pfc_air == Inf]<-0
-#write.csv(DIV_pfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/DIV_pfc_CO2eq_air_",year_current,"_pxp.csv"))
 TIV_pfc_air = as.vector(DIV_pfc_air) %*% L
 
 write.csv(TIV_pfc_air, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_pfc_CO2eq_air_",year_current,"_pxp.csv"))
@@ -1149,7 +1090,6 @@ DIV_e_u[is.na(DIV_e_u)]=0
 DIV_e_u[DIV_e_u == Inf]<-0
 TIV_e_u = as.vector(DIV_e_u) %*% L
 
-#saveRDS(TIV_e_u, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_energy_carrier_use_",year_current,".rds"))
 write.csv(TIV_e_u, paste0(data_dir_exiobase, "/IOT_",year_current,"_pxp/TIV_energy_carrier_use_",year_current,"_pxp.csv"))
 
 TIV_breakdown_e_u = as.vector(DIV_e_u) * L
@@ -1166,6 +1106,18 @@ write.csv(TIV_country_breakdown_e_u_w_labels, paste0(data_dir_exiobase, "/IOT_",
 
 # isf
 
+HBS: European household budget survey from EUROSTAT, macro-data, from : https://ec.europa.eu/eurostat/web/household-budget-surveys/database [acessed on 22.05.2020]
+
+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)'
+
+
+1) The 'lfst_hhnhtych' table from EUROSTAT, selecting all available years and total households [accessed on 04.05.2020]. 
+
+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].
+
+
 ```{r income-stratified-footprints, eval = FALSE}
 
 # income-stratified-footprints directory