%% Creating 2015 energy and carbon footprints for the EU HBS analysis
%For the 2010 version please contact Diana Ivanova
%Diana Ivanova, amended by Milena Buchs to create more aggregate
%consumption categories with 2015 data



%% 2.  Quantifying energy & emissions

% 2.1 Prepare EXIOBASE's INDIRECT energy and carbon footprints on a country
% level - based on 200 product categories
% Load EXIOBASE data: Import A.txt, Y.txt, F.txt and F_hh.txt from EXIOBASE
% into Matlab
% Import the files as numeric matrices, pxp, 2015
% https://zenodo.org/record/4277368#.YLktD_lKjIW; https://zenodo.org/record/4588235#.YLksZflKjIU


cd 'enter your file path'
load('IO_start.mat')

number_countries = length(meta.raw_countries);
countries=meta.raw_countries;

load('A37.mat')
load('Y37.mat')
load('F37.mat')
load('Fhh37.mat')

tic
Yhh = Y37(:,1:7:end); %only import households here
Yhh_200 = aggarr(Yhh, 200,49,1);
L= inv(eye(size(A37))-A37);
X = L*Y37;

% construct own characterisation matrix for net energy - taking the rows between 1339 and 1412
char_energy= xlsread('DIS_background_files.xlsx','char_energy','B5:BMR5'); 
char_ghg = xlsread('DIS_background_files.xlsx','char_energy','B6:BMR6');

Sfull = F37./sum(X,2)';
Sfull(isnan(Sfull))=0;

S_name = ["Energy", "Carbon"]
S(1,:) = char_energy*Sfull; %compare with meta.labsC, net energy in EXIOBASE
S(2,:) = char_ghg*Sfull;
S(isnan(S))=0;

%% create INDIRECT energy and carbon footprints per 200 consumption category by country

for j = 1:length(countries)
    for jj =1:size(S,1)
        temp = S(jj,:) * L * diag(Yhh(:,j));
        temp2=aggarr(temp,1,200,1);
        FP_ind_200.(S_name{jj})(j,:) =temp2;
    end
end

%write indirect emissions and energy to an excel file (will need to
%separate out energy and emission tabs)
writetable(struct2table(FP_ind_200), 'FP_ind_200.xlsx')
clear temp temp2 X A37 Y37 Sfull L 









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.2 Prepare EXIOBASE's DIRECT energy and carbon footrpints on a country level
% - for 200 product categories

%direct energy use
ext_bridge= xlsread('DIS_background_files.xlsx','ext_bridge','C2:GT75'); 
Fhh_char = diag(char_energy)*Fhh37;
FP_dir_200.Energy = Fhh_char(1339:1412,1:7:end)'*ext_bridge;
writetable(struct2table(FP_dir_200), 'FP_dir_200.xlsx')




%%
%direct carbon footprint (first total footprint per country)
%1707 are energy/carbon extensions
%343 are countries times final demand sectors, the 7th are households, only
%need these

Fc=Fhh37'; %this is the carbon footprint file
b=char_ghg'; %this is the file for the GHG emission selectors (which env indicators are included)
FP_dir.Carbon=Fc*b; %create a structure by multiplying the carbon footprints with the GHG indicators
C_dir=FP_dir.Carbon(1:7:end,:); %create a table with total direct GHG per country
%use electricity emissions from indirect carbon emissions (already there
%for each country)
%convert gj to emissions from direct energy
writematrix(C_dir, 'Carbon_dir_tot.xlsx')




%direct carbon footprint for consumption categories with Richard Wood's
%file
hhld_emis =IO.hhld_emis;
for i = 1: size(hhld_emis,3)
    temp(:,:,i) =char_ghg*hhld_emis(:,:,i); 
end

hhld_emis_ghg=squeeze(temp(1,:,:));
FP_cdir_200.Carbon = hhld_emis_ghg(:,1:7:end)';







%% 2.3 Compute energy footprints per country for aggregate consumption categories
% Import the consumption category concordance bridging the 200 EXIOBASE products to 11 consumption categories
cat_bridge = xlsread('DIS_background_files.xlsx','bridge_mb','C3:M202');

%% Multiply the INDIRECT footprint matrix with the conversion bridge matrix
% to get indirect footprints by category and country
% this is for both energy and carbon
for jj = 1:length(S_name)
        FP_ind_cat.(S_name{jj}) = FP_ind_200.(S_name{jj})*cat_bridge;
end 

%write this to an excel file, will need to split up tabs on there
writetable(struct2table(FP_ind_cat), 'FP_ind_cat.xlsx');





%% DIRECT energy footprints per category (calculate direct emissions from
%%here for gas, other home energy fuels and motor fuel

%Multiply the DIRECT energy footprint matrix with the conversion bridge matrix
% to get indirect footprints by category and country
E_name = ["Energy"]
for jj = 1:length(E_name)
        FP_endir_cat.(E_name{jj}) = FP_dir_200.(E_name{jj})*cat_bridge;
end 

%write this to an excel file, this is for direct energy only
writetable(struct2table(FP_endir_cat), 'FP_endir_cat.xlsx');









% get DIRECT GHG footprints 2015 per consumption category
C_name = ["Carbon"]
for jj = 1:length(C_name)
        FP_cdir_cat.(C_name{jj}) = FP_cdir_200.(C_name{jj})*cat_bridge;
end 


%write this to an excel file, this is for direct energy only
writetable(struct2table(FP_cdir_cat), 'FP_cdir_cat.xlsx');







% DIRECT GHG footprints 2010 per consumption category (49x200 file)
load('FPdirGHG2010.mat')
C_name = ["Carbon"]
D.Carbon=FPdirGHG492002010;
for jj = 1:length(C_name)
        FP_cdir2010_cat.(C_name{jj}) = D.(C_name{jj})*cat_bridge;
end 

%write this to an excel file, this is for direct energy only
writetable(struct2table(FP_cdir2010_cat), 'FP_cdir2010_cat.xlsx');

