%% Run this script after DIS_new_bridge.m
% Database: EXIOBASE 3.7

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script is meant to do the following:

% 1.  Harmonise the HBS expenditure and EXIOBASE
% 1.1 Bridge the HBS and EXIOBASE on a country level
% 1.2 Import and bridge to EXIOBASE classification on a person level
% 1.3 Account for under- and over-reporting of expenditure in the HBS
% 1.4 Transform from PPs to BPs

% 2.  Quantifying energy
% 2.1 Prepare EXIOBASE's indirect energy multipliers on a country level
% 2.2 Prepare EXIOBASE's direct energy multipliers on a country level
% 2.3 Compute household-level indirect energy footprints
% 2.4 Compute household-level direct energy footprints
% 2.5 Compute and export total energy footprints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

%% 1. Harmonise the HBS expenditure and EXIOBASE
% 1.1 Bridge the HBS and EXIOBASE on a country level
load('new_bridge.mat')

% import the survey expenditure in purchaser prices and remove spending
% from the survey where there is none 
su_ybp(:,sum(CES_pp)==0)=0; 
su_ypp(:,sum(CES_pp)==0)=0; 
su_ytax(:,sum(CES_pp)==0)=0; 
su_ytrade(:,sum(CES_pp)==0)=0; 
su_ytrans(:,sum(CES_pp)==0)=0; 
su_diff=CES_pp./su_ypp;
su_diff(isnan(su_diff))=0;
su_diff(isinf(su_diff))=0;
CES_pp (su_diff==0) = 0;

% Take care of under- and over-reporting for each country on a product level
% some sectors - e.g. 160 Sea and water travel - require a different procedure as it is only a small % of the sample that consumes it 
% this is especially important for the sectors with the highest intensities
% first, we rescale the pp-expenditure vector for under- and over-reporting
% to match the mrio data
% second, we account for the leftover difference between the new exp vector and the
% mrio, which is caused by zero-exp in survey classification (for some
% countries). Diff could be both positive or negative.

YHH_pp_cap(:,sum(CES_pp)==0)=0; 
diff_ratio = YHH_pp_cap./CES_pp;
diff_ratio(isnan(diff_ratio)) = 0; 
diff_ratio(isinf(diff_ratio)) = 1; 

%outlier treatment
diff_ratio_adj_out = (YHH_pp_cap-CES_pp)./((YHH_pp_cap+CES_pp)./2);
diff_ratio_adj_out(isnan(diff_ratio_adj_out)) = 0; 
mean =mean(diff_ratio_adj_out(:)) %close to zero
std = std(diff_ratio_adj_out(:))
diff_ratio_adj_noout = diff_ratio_adj_out;
diff_ratio_adj_noout(diff_ratio_adj_noout>mean+1.96*std)=mean+1.96*std; %corresponding to 95% probability
diff_ratio_adj_noout(diff_ratio_adj_noout<mean-1.96*std)=mean-1.96*std;

%replace the outliers in diff_ratio
CES_pp_noout = (YHH_pp_cap.*(2-diff_ratio_adj_noout))./(2+diff_ratio_adj_noout);
diff_ratio_noout = YHH_pp_cap./CES_pp_noout;
diff_ratio_noout(isnan(diff_ratio_noout)) = 0; 

CES_pp_adj =CES_pp.*diff_ratio_noout; %% rescale depending on exio estimate

% Identify margins sectors:
trdMsectors = 152:155;
tspMsectors = 157:163;

%% 1.2 Import and bridge to EXIOBASE classification on a person level
% We import the data for all EXIOBASE countries except for Germany and
% Sweden. These countries had different product classification in their HBS
% and since they did not provide any data on disability, we have to exclude
% them from the analysis anyway. In any case, don't be alarmed about their
% zero-consumption.

% The microdata table (275247x66 table) contains household expenditure in several COICOP categories
% It also contains the variables HA04 (household ID), COUNTRY (country ID)
% and HB05 (household size). For a list of all variables, see the
% background files: "DIS_background_files","microdata_structure"
microdata = readtable('micro_data.csv');

% single out the expenditure variables
ces_meur_mic = (table2array(microdata(:,1:63)))/1000000;
ces_meur_mic = ces_meur_mic';
ces_meur_mic(isnan(ces_meur_mic))=0;

% create an index based on the country variable depending on the number of
% observations from that country
% e.g. for ex: 1 is for BE, 7178 is for BG, etc
names =table2array(microdata(:,65));
[abbr, Ia]= unique(names,'stable');
hhsize = table2array(microdata(:,66));

% transform to 200 sectors and per capita expenditure
% the loop takes care of all the countries but the last one
% abbr_index is the index for countries with the survey order
CEShh_pp_mic = zeros(size(ces_meur_mic,2),200);
CES_pp_mic = zeros(size(ces_meur_mic,2),200);

for i =1:length(abbr)-1
    
    for j = 1:length(countries)
        if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i))))&&i~=5&&i~=23)
            country_code = cell2mat(meta.raw_countries(j));
            CEShh_pp_mic(Ia(i):Ia(i+1)-1,:)=ces_meur_mic(:,Ia(i):Ia(i+1)-1)'*new_bridge.(country_code);
            CES_pp_mic(Ia(i):Ia(i+1)-1,:)=CEShh_pp_mic(Ia(i):Ia(i+1)-1,:)./hhsize(Ia(i):Ia(i+1)-1,:);
        end
    end
end

for j = 1:length(countries)
    if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end)))))
        country_code = cell2mat(meta.raw_countries(j));
        CEShh_pp_mic(Ia(end):end,:)=ces_meur_mic(:,Ia(end):end)'*new_bridge.(country_code);
        CES_pp_mic(Ia(end):end,:)=CEShh_pp_mic(Ia(end):end,:)./hhsize(Ia(end):end,:);
    end
end

% account for zero exp sectors in the MRIO
for i =1:length(abbr)-1
    for j = 1:length(countries)
        for k=1:200
            if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i))))&&su_diff(k,j)==0)
                CEShh_pp_mic(Ia(i):Ia(i+1)-1,k) = 0;
                CES_pp_mic(Ia(i):Ia(i+1)-1,k) = 0;
            end
        end
    end
end
for j = 1:length(countries)
       for k=1:200
    if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end))))&&su_diff(k,j)==0)
        CEShh_pp_mic(Ia(end):end,k) = 0;
            CES_pp_mic(Ia(end):end,k) = 0;
        end
    end
end

%% 2.2 Account for under- and over-reporting of expenditure in the HBS
diff_ratio_noout = diff_ratio_noout';

CES_pp_adj_mic = zeros(size(CES_pp_mic,1),200);
for i =1:length(abbr)-1
    for j = 1:length(countries)
        if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i)))))
            CES_pp_adj_mic(Ia(i):Ia(i+1)-1,:) =CES_pp_mic(Ia(i):Ia(i+1)-1,:).*repmat(diff_ratio_noout(j,:),Ia(i+1)-Ia(i),1);
        end
    end
end
for j = 1:length(countries)
    if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end)))))
        CES_pp_adj_mic(Ia(end):end,:)=CES_pp_mic(Ia(end):end,:).*repmat(diff_ratio_noout(j,:),size(CES_pp_mic,1)-Ia(end)+1,1);
    end
end

% If there are some inconsistencies in terms of averages - they come from
% zero expenditure in the survey and positive MRIO entry. Even though the
% bridges were optimised, they may still be such cases. No further
% additions or substractions are done here - in order to avoid negative
% footprints. 

%% 2.3 Transform from purchaser prices (PP) to basic prices (BP)
% We do that in consideration of taxes, transport and trade margins
su_ypp=su_ypp';

su_diff_mic = zeros(size(CES_pp_mic,1),200);
CES_tax_mic = zeros(size(CES_pp_mic,1),200);
CES_trade_mic = zeros(size(CES_pp_mic,1),200);
CES_trans_mic = zeros(size(CES_pp_mic,1),200);

for i =1:length(abbr)-1
    for j = 1:length(countries)
        if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i)))))
            su_diff_mic(Ia(i):Ia(i+1)-1,:)=CES_pp_adj_mic(Ia(i):Ia(i+1)-1,:)./repmat(su_ypp(j,:),Ia(i+1)-Ia(i),1);
            su_diff_mic(isnan(su_diff_mic))=0;
            CES_tax_mic(Ia(i):Ia(i+1)-1,:) = su_diff_mic(Ia(i):Ia(i+1)-1,:).*su_ytax(:,j)';
            CES_trade_mic(Ia(i):Ia(i+1)-1,:) = su_diff_mic(Ia(i):Ia(i+1)-1,:).*su_ytrade(:,j)';
            CES_trans_mic(Ia(i):Ia(i+1)-1,:) = su_diff_mic(Ia(i):Ia(i+1)-1,:).*su_ytrans(:,j)';
        end
    end
end
for j = 1:length(countries)
    if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end)))))
        su_diff_mic(Ia(end):end,:)=CES_pp_adj_mic(Ia(end):end,:)./repmat(su_ypp(j,:),size(su_diff_mic,1)-Ia(end)+1,1);
        su_diff_mic(isnan(su_diff_mic))=0;
        CES_tax_mic(Ia(end):end,:) = su_diff_mic(Ia(end):end,:).*su_ytax(:,j)';
        CES_trade_mic(Ia(end):end,:) = su_diff_mic(Ia(end):end,:).*su_ytrade(:,j)';
        CES_trans_mic(Ia(end):end,:) = su_diff_mic(Ia(end):end,:).*su_ytrans(:,j)';
    end
end

su_diff_ybp_prel_mic = CES_pp_adj_mic-CES_tax_mic-CES_trade_mic-CES_trans_mic;

%Calculate total margins:
totaltrdM_CREEA_mic = -sum(CES_trade_mic(:,trdMsectors),2);
totaltspM_CREEA_mic = -sum(CES_trans_mic(:,tspMsectors),2);

%Calculate relative distribution among margins sectors:
trdMsectorweights_mic = -CES_trade_mic(:,trdMsectors)./totaltrdM_CREEA_mic;
tspMsectorweights_mic = -CES_trans_mic(:,tspMsectors)./totaltspM_CREEA_mic;
trdMsectorweights_mic(isnan(trdMsectorweights_mic))=0;
tspMsectorweights_mic(isnan(tspMsectorweights_mic))=0;

% Redistribute margins to the margins sectors:
% Calculate total margins:
totaltrdM_CES_mic = sum(CES_trade_mic,2);
totaltspM_CES_mic = sum(CES_trans_mic,2);

% Distribute totals on individual margins sectors:
CES_bp_adj_mic = su_diff_ybp_prel_mic;

for i=1:meta.NCOUNTRIES
    CES_bp_adj_mic(:,trdMsectors) = CES_bp_adj_mic(:,trdMsectors)+trdMsectorweights_mic.*totaltrdM_CES_mic;
    CES_bp_adj_mic(:,tspMsectors) = CES_bp_adj_mic(:,tspMsectors)+tspMsectorweights_mic.*totaltspM_CES_mic;
end

clear totaltrdM_CES_mic totaltspM_CES_mic tspMsectorweights_mic 
clear su_diff_tax su_diff_trade CES_bp_prelim su_diff_trans su_diff_ybp_prel_mic su_ybp su_ypp su_ytax su_ytrade su_ytrans totaltrdM_CES totaltrdM_CREEA totaltrdM_CREEA_mic tspMsectors trdMsectorweights_mic trdMsectorweights trdMsectors totaltspM_CES totaltspM_CREEA totaltspM_CREEA_mic  i j an
clear CES_pp_adj_mic CES_pp_mic CES_tax_mic CES_trade_mic CES_trans_mic CEShh_pp_mic hhsize

%% 2.  Quantifying energy
% 2.1 Prepare EXIOBASE's indirect energy multipliers on a country level
% Load EXIOBASE data: Import A.txt, Y.txt, F.txt and F_hh.txt from EXIOBASE
% Import the files as numeric matrices, pxp, 2010
% https://zenodo.org/record/4277368

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

tic
Yhh = Y37(:,1:7:end);
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 = zeros(size(char_energy));

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;

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

for jj =1:size(S,1)
    M_ind_200.(S_name{jj})= FP_ind_200.(S_name{jj})'./Yhh_200;
    M_ind_200.(S_name{jj})(isnan(M_ind_200.(S_name{jj})))=0;
end
toc

% calculate an EU-average for the hot water sector which otherwise gives
% huge multipliers for some countries (PL, BG, LT)

for i=1:28
    for jj =1:size(S,1)
        M_ind_200.(S_name{jj})(148,i) = sum(FP_ind_200.(S_name{jj})(1:28,148))./sum(Yhh_200(148,1:28));
    end
end

clear temp temp2 X A37 Y37 Sfull L 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.2 Prepare EXIOBASE's direct energy multipliers on a country level
% Calculate direct emissions, (a) proportionate to expenditure or (b)
% fixed per capita for the few products with no expenditure 

%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;

FP_dir_200.Carbon = zeros(size(FP_dir_200.Energy));

%direct emission multipliers - with direct emissions varying by expenditure
for jj =1:length(S_name)
M_d_200.(S_name{jj}) = FP_dir_200.(S_name{jj})./Yhh_200';
M_d_200.(S_name{jj})(isnan(M_d_200.(S_name{jj})))=0;
M_d_200.(S_name{jj})(isinf(M_d_200.(S_name{jj})))=0;
end

% calculate an EU-average for the hot water sector which otherwise gives
% huge multipliers for some countries (DK, NL and LV)
for i=1:28
    for jj =1:size(S,1)
        M_d_200.(S_name{jj})(i,148) = sum(FP_dir_200.(S_name{jj})(1:28,148))./sum(Yhh_200(148,1:28));
    end
end

% calculate fixed emissions for the zero- and low-expenditure products
FP_d_fixed.Energy = zeros(49,200);
FP_d_fixed.Carbon = zeros(49,200);
for jj =1:length(S_name)
    for i=[16,17,20,22:27,29,59:61,64,83]
        for j=1:49
            M_d_200.(S_name{jj})(j,i)= 0;
            FP_d_fixed.(S_name{jj})(j,i) = FP_dir_200.(S_name{jj})(j,i)./pop(j,:);
        end
    end
end

%% 2.3 Compute household-level indirect energy footprints
% Import the consumption category concordance bridging the 200 EXIOBASE products to 18 consumption categories
cat_bridge_disability = xlsread('DIS_background_files.xlsx','bridge','C3:S202');

% Calculate indirect energy footprints for households by multiplying the household expenditure by product with the product-level indirect energy multiplier. 
% Then aggregate the indirect energy footprints by category and in total
for jj = 1:length(S_name)
for i =1:length(abbr)-1
    for j = 1:length(countries)
                if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i)))))
            FP_ind_mic.(S_name{jj})(Ia(i):Ia(i+1)-1,:) = CES_bp_adj_mic(Ia(i):Ia(i+1)-1,:)*M_ind_200.(S_name{jj})(:,j);
            FP_ind_200_mic.(S_name{jj})(Ia(i):Ia(i+1)-1,:) = CES_bp_adj_mic(Ia(i):Ia(i+1)-1,:)*diag(M_ind_200.(S_name{jj})(:,j));
        end
    end
end
end 

end_end = size(CES_bp_adj_mic,1);

for jj = 1:length(S_name)
for j = 1:length(countries)
    if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end)))))
        FP_ind_mic.(S_name{jj})(Ia(end):end_end,:) = CES_bp_adj_mic(Ia(end):end_end,:)*M_ind_200.(S_name{jj})(:,j);
        FP_ind_200_mic.(S_name{jj})(Ia(end):end_end,:) = CES_bp_adj_mic(Ia(end):end_end,:)*diag(M_ind_200.(S_name{jj})(:,j));
    end
end
end

for jj =1:length(S_name)
     FP_ind_cat_mic.(S_name{jj}) = FP_ind_200_mic.(S_name{jj})*cat_bridge_disability;
     FP_ind_mic.(S_name{jj}) = sum(FP_ind_200_mic.(S_name{jj}),2);  
end

%% 2.4 Compute household-level direct energy footprints
% (a) variable: proportionate to expenditure (M_d_200) or (b) % fixed per capita: for the few products with no expenditure (FP_d_fixed)
for jj = 1:length(S_name)
    for i =1:length(abbr)-1
        for j = 1:length(countries)
            if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(i)))))
           
                FP_d_200_mic.(S_name{jj})(Ia(i):Ia(i+1)-1,:) = CES_bp_adj_mic(Ia(i):Ia(i+1)-1,:)*diag(M_d_200.(S_name{jj})(j,:))+repmat(FP_d_fixed.(S_name{jj})(j,:),Ia(i+1)-Ia(i),1);
            end
        end
    end
end

for jj = 1:length(S_name)
    for j = 1:length(countries)
        if(strcmp(cell2mat(countries(j)),cell2mat(cellstr(abbr(end)))))
            FP_d_200_mic.(S_name{jj})(Ia(end):end_end,:) = CES_bp_adj_mic(Ia(end):end_end,:)*diag(M_d_200.(S_name{jj})(j,:))+repmat(FP_d_fixed.(S_name{jj})(j,:),size(CES_bp_adj_mic,1)-Ia(end)+1,1);
            
        end
    end
end

for jj =1:length(S_name)
     FP_d_cat_mic.(S_name{jj}) = FP_d_200_mic.(S_name{jj})*cat_bridge_disability;
     FP_d_mic.(S_name{jj}) = sum(FP_d_200_mic.(S_name{jj}),2);  
end

%% 2.5 Compute and export total energy footprints
for jj =1:length(S_name)
FP_cat_mic.(S_name{jj}) = FP_d_cat_mic.(S_name{jj}) + FP_ind_cat_mic.(S_name{jj});
FP_mic.(S_name{jj}) = FP_d_mic.(S_name{jj}) + FP_ind_mic.(S_name{jj});
end 

csvwrite('DIS_energy_footprints_cat.csv',FP_cat_mic.Energy)
csvwrite('DIS_energy_footprints_total.csv',FP_mic.Energy)