%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ExtractMultipliers for Norway from EXIOBASE.
%
% The end purpose of this script is to produce multipliers per unit of
% imports for Norway
% A number of steps are undertaken, and the final results of this script are
% the estimation of emissions embodied in imports, and imports.
% Multipliers can then be calculated as the division between the emissions
% embodied in imports and imports.
% The calculation is done for both total supply chain emissions
% (multipliers) and for direct emissions (direct emissions associated with 
% imported goods). The latter is useful for extracting scope 2 emissions.
%
%
% Richard Wood
% 26.11.2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1 (Manual, may be automated later through API on Zenodo
% EXIOBASE data must be downloaded from the web.
% EXIOBASE v3.8.2 provides estimates to 2019 available at:
% https://ntnu.box.com/s/e9yr097f02aihj5l1fcgxn31g5bbd5uk - these are matlab data files
% https://zenodo.org/record/5589597 - these are txt files that need converting to matlab structures
% This script is set up to use the data in matlab structures, and using the
% txt will require conversion.

% You need these files:
% IOT_2018_ixi.mat
% IOT_ZFdom_2018_ixi.mat

% The capital use matrices are also needed,  Source: https://zenodo.org/record/5787624 
% Kbar_exio_v3_8_2_2018ixi.mat

mrio_data_path='C:\Users\ALEC\Documents\DFØ\working\IOT_2018_ixi\'

% Source: https://zenodo.org/record/5787624
capital_data_path = 'C:\Users\ALEC\Documents\DFØ\working\capital\'

%set years to process
extractyears=[2018]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2 - set up meta data - work out which emission categories we are
% interested in:

load([mrio_data_path,'IOT_2018_ixi.mat'],'meta')

%specify the emission name
EmissNames = {'GHG emissions AR5 (GWP100)|kg CO2 eq.||GWP100 (IPCC, 2010)'}
% if the above does not work, try uncommenting the line below, as I believe
% different matlab verssions have different string conventions
% EmissNames = {'"GHG emissions AR5 (GWP100)|kg CO2 eq.||GWP100 (IPCC, 2010)"'} 


%%
% find the index of the specified emissions in EXIOBASE meta's
% characterisation labels:
n_emis=length(EmissNames);
for i=1:n_emis
    EmissIndx(i)=find(strcmp(EmissNames{i},meta.labsC));
    tmp=split(meta.labsC(EmissIndx(i)),'|');
    EmissUnit_raw(i)=tmp(2);
    %     % if matching F (the basic emission rows) rather than C (the
    %     charachterised emission rows - e.g. with GWP applied)
    %     if isempty(meta.labsF{EmissIndx(i),3})
    %         EmissUnit_raw(i)=meta.labsF(EmissIndx(i),2);
    %     else
    %         EmissUnit_raw(i)=meta.labsF(EmissIndx(i),3);
    %     end
end

%%
% specify some generic meta information for variables later (using this
% standard arrangment:
meta_generic_dim3.dim1=meta.secLabsZ(1:meta.NSECTORS,1:4);
meta_generic_dim3.dim2=EmissNames;
meta_generic_dim3.dim3=[meta.countries,meta.raw_countries,meta.countrynames];

meta_generic_dim3_full=meta_generic_dim3;
meta_generic_dim3_full.dim1=meta.secLabsZ(:,1:4);

meta_generic_dim2.dim1=meta.secLabsZ(1:meta.NSECTORS,1:4);
meta_generic_dim2.dim2=[meta.countries,meta.raw_countries,meta.countrynames];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Extract the relevant data and calculate footprints


yrcount=0;
for yr=extractyears
    yrcount=yrcount+1;
    
    %initialise the variables:
    footprint_trade=meta_generic_dim3_full;
    footprint_fd=meta_generic_dim3_full;
    footprint_trade_aggregated=meta_generic_dim3;
    footprint_fd_aggregated=meta_generic_dim3;
    trade_aggregated=meta_generic_dim2;
    fd_aggregated=meta_generic_dim2;
    
    %load the main MRIO file, and also a file that contains the bilateral
    %trade data:
    load([mrio_data_path,'IOT_',num2str(yr),'_ixi.mat'])
    load([mrio_data_path,'IOT_ZFdom_',num2str(yr),'_ixi.mat'],'trade_flat')
    load([capital_data_path,'Kbar_exio_v3_8_2_',num2str(yr),'ixi.mat'])
    % calculate Leontief multipliers for all Emissions at once:
    multipliers=((IO.char(EmissIndx,:)*IO.S)/(eye(meta.Zdim)-full(IO.A)))';
    
    KCoef=KbarCfc_ixi/sparse(diag(sum(KbarCfc_ixi,1)));
    KCoef(isnan(KCoef))=0;
    KCoef=KCoef*diag(IO.S(9,:));
    
    multipliers_with_K=((IO.char(EmissIndx,:)*IO.S)/(eye(meta.Zdim)-full(IO.A+KCoef)))';
    direct_intensities=(IO.char(EmissIndx,:)*IO.S)';
    direct_emissions_full=(IO.char(EmissIndx,:)*IO.F)';
    
    
    %specify final demand and footprints, all in full 163*49 industry
    %classification (needs to be aggregated later) - rows show the industry
    %and country of last stage of production:
    for ncnt=1:49
        finaldemand.val(:,ncnt)=full(sum(IO.Y(:,(ncnt-1)*7+1:ncnt*7),2));
        direct_emissions.val(:,ncnt)=direct_emissions_full((ncnt-1)*163+1:ncnt*163);
    end
    for ncnt=1:49
        footprint_fd.val(:,:,ncnt)=multipliers.*repmat(finaldemand.val(:,ncnt),1,n_emis);
        footprint_fd_w_k.val(:,:,ncnt)=multipliers_with_K.*repmat(finaldemand.val(:,ncnt),1,n_emis);
        direct_footprint_fd.val(:,:,ncnt)=direct_intensities.*repmat(finaldemand.val(:,ncnt),1,n_emis);
    end
    for ncnt=1:49
        footprint_trade.val(:,:,ncnt)=multipliers.*repmat(trade_flat(:,ncnt),1,n_emis);
        footprint_trade_w_k.val(:,:,ncnt)=multipliers_with_K.*repmat(trade_flat(:,ncnt),1,n_emis);
        direct_footprint_trade.val(:,:,ncnt)=direct_intensities.*repmat(trade_flat(:,ncnt),1,n_emis);
    end
    
    % now lets create an aggregation matrix that will (second line) aggregate all the
    % imports from different countries into total imports into the importing country:
    agg_source_country=kron(ones(49,1),eye(meta.NSECTORS));
    
    % now aggregate the trade and final demand:
    trade_aggregated.val=agg_source_country'*trade_flat;
    fd_aggregated.val=agg_source_country'*finaldemand.val;
    
    
    for ncnt=1:meta.NCOUNTRIES
        % aggregate final demand footprints
        footprint_fd_aggregated.val(:,:,ncnt)=agg_source_country'*footprint_fd.val(:,:,ncnt);
        footprint_fd_aggregated_w_k.val(:,:,ncnt)=agg_source_country'*footprint_fd_w_k.val(:,:,ncnt);
        direct_footprint_fd_aggregated.val(:,:,ncnt)=agg_source_country'*direct_footprint_fd.val(:,:,ncnt);
        
        % do the same aggregation for footprints of trade
        footprint_trade_aggregated.val(:,:,ncnt)=agg_source_country'*footprint_trade.val(:,:,ncnt);
        footprint_trade_aggregated_w_k.val(:,:,ncnt)=agg_source_country'*footprint_trade_w_k.val(:,:,ncnt);
        direct_footprint_trade_aggregated.val(:,:,ncnt)=agg_source_country'*direct_footprint_trade.val(:,:,ncnt);
    end
    
    % switch on to save detailed matlab data if desired:
    %     save(['extracted_exiobase_data_3_8_',num2str(yr),'.mat'],'multipliers', 'finaldemand','footprint_fd','footprint_fd_aggregated','footprint_trade','trade_flat',...
    %         'trade_aggregated','footprint_trade_aggregated','meta')
    
    
    %%
    % now for the dump of the output for the coupled model - the aggregated
    % trade and footprint data which will be remapped to NO classification
    % in the coupled model.
    
    %  do 3 regions:
    %  Nor
    %  EU
    %  GLobal
    cnt_index(1,:)=strcmp(meta.countries,'NOR'); regname(1)={'NO'};
    cnt_index(2,:)=logical([repmat(1,28,1);repmat(0,21,1)])'; regname(2)={'EU'};
    cnt_index(3,:)=logical([repmat(1,49,1)])'; regname(3)={'WO'};
    k_suffix=[{''},{'_w_Capital'}];
    
    for inkl_cap=1:2
        for nreg=1%:3
            if inkl_cap==1
                % simply extract the data for region of interest:
                trade_aggregated.reg=sum(trade_aggregated.val(:, cnt_index(nreg,:)),2);
                footprint_trade_aggregated.reg=sum(footprint_trade_aggregated.val(:,:, cnt_index(nreg,:)),3);
                fd_aggregated.reg=sum(fd_aggregated.val(:, cnt_index(nreg,:)),2);
                footprint_fd_aggregated.reg=sum(footprint_fd_aggregated.val(:,:, cnt_index(nreg,:)),3);
            else
                % simply extract the data for region of interest:
                trade_aggregated.reg=sum(trade_aggregated.val(:, cnt_index(nreg,:)),2);
                footprint_trade_aggregated.reg=sum(footprint_trade_aggregated_w_k.val(:,:, cnt_index(nreg,:)),3);
                fd_aggregated.reg=sum(fd_aggregated.val(:, cnt_index(nreg,:)),2);
                footprint_fd_aggregated.reg=sum(footprint_fd_aggregated_w_k.val(:,:, cnt_index(nreg,:)),3);
            end
            
            direct_footprint_trade_aggregated.reg=sum(direct_footprint_trade_aggregated.val(:,:, cnt_index(nreg,:)),3);
            direct_footprint_fd_aggregated.reg=sum(direct_footprint_fd_aggregated.val(:,:, cnt_index(nreg,:)),3);
            
            direct_emissions.reg=direct_emissions.val(:,cnt_index(nreg,:));
            
            % and write data out into a table in the "Multiplier" Spreadsheet
            
            % Firstly, the import data
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],yr,['EXIOrawImp_',num2str(yr)],'a1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],{'Imports'},['EXIOrawImp_',num2str(yr)],'B1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],{'million Euro'},['EXIOrawImp_',num2str(yr)],'B2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],trade_aggregated.dim1.CodeNr,['EXIOrawImp_',num2str(yr)],'A3')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],trade_aggregated.reg,['EXIOrawImp_',num2str(yr)],'B3')
            % Secondly the footprints of imports:
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissNames,['EXIOrawImp_',num2str(yr)],'C1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissUnit_raw,['EXIOrawImp_',num2str(yr)],'C2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],footprint_trade_aggregated.reg,['EXIOrawImp_',num2str(yr)],'C3')
            % third the direct footprint of imports
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissNames,['EXIOrawImp_',num2str(yr)],'D1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissUnit_raw,['EXIOrawImp_',num2str(yr)],'D2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],direct_footprint_trade_aggregated.reg,['EXIOrawImp_',num2str(yr)],'D3')
            
            
            % Next, the final demand data:
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],yr,['EXIOrawFd_',num2str(yr)],'A1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],{'Final Demand'},['EXIOrawFd_',num2str(yr)],'B1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],{'million Euro'},['EXIOrawFd_',num2str(yr)],'B2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],fd_aggregated.dim1.CodeNr,['EXIOrawFd_',num2str(yr)],'A3')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],fd_aggregated.reg,['EXIOrawFd_',num2str(yr)],'B3')
            % Secondly the footprints of fd:
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissNames,['EXIOrawFd_',num2str(yr)],'C1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissUnit_raw,['EXIOrawFd_',num2str(yr)],'C2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],footprint_fd_aggregated.reg,['EXIOrawFd_',num2str(yr)],'C3')
            % third the direct footprint of fd
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissNames,['EXIOrawFd_',num2str(yr)],'D1')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],EmissUnit_raw,['EXIOrawFd_',num2str(yr)],'D2')
            xlswrite(['Multipliers_',regname{nreg},'_bp_EXIOBASE',k_suffix{inkl_cap},'.xlsx'],direct_footprint_fd_aggregated.reg,['EXIOrawFd_',num2str(yr)],'D3')
            
            
            % Direct (production based) emissions
            xlswrite(['Direct_Emissions_',regname{nreg},'_EXIOBASE.xlsx'],yr,['EXIO_prod_',num2str(yr)],'A1')
            xlswrite(['Direct_Emissions_',regname{nreg},'_EXIOBASE.xlsx'],trade_aggregated.dim1.CodeNr,['EXIO_prod_',num2str(yr)],'A3')
            xlswrite(['Direct_Emissions_',regname{nreg},'_EXIOBASE.xlsx'],EmissNames,['EXIO_prod_',num2str(yr)],'B1')
            xlswrite(['Direct_Emissions_',regname{nreg},'_EXIOBASE.xlsx'],EmissUnit_raw,['EXIO_prod_',num2str(yr)],'B2')
            xlswrite(['Direct_Emissions_',regname{nreg},'_EXIOBASE.xlsx'],direct_emissions.reg,['EXIO_prod_',num2str(yr)],'B3')
            
        end %regional loop
    end % include capital loop
end % year loop


