%Lina I. Brand Correa and Anne Owen, 2017

% For more information on the context, rationale and methodology, please
% refer to the paper: Brand-Correa, L.I., Broackway, P.E., Copeland, C.L.,
% Foxon, T.J., Owen, A., and Taylor, P.G. (2017) "Developing and Input-Output 
% Based method to Estimate a National-Level Energy Return on Investment (EROI)"
% Energies, 10, 520.

% Variable names correspond to those used in the paper, unless stated otherwise.

clear all

%% Load the data

load Z.mat; % transaction matrix 
load Y.mat; % final demand (two columns, one for the UK one for RoW)
load F.mat; % energy 'production', i.e. our energy extension (f in the paper)

% See Figure A1 for the basic structure of the UK's MRIO

%% EROI EiE Generate and preallocate variables 

year = 17; % years (1997 to 2013)
i = 106; % sectors
s = 6; % regions (UK, Rest of Europe, Middle East, China and RoW)
c = 4; % position of UK's energy sector
g = 5; % position of UK's energy sector
ensec = g-c+1; % amount of energy sectors

% preallocating matrices with zeros for speed of calculations
A = zeros(i*s,i*s,year); 
L = zeros(i*s,i*s,year);  
y = zeros(i*s,1,year); 
x = zeros(1,i*s,year);
e = zeros(1,i*s,year);
Foot = zeros(i*s,i*s,year);

Z0 = zeros(i*s,i*s,year);
A0 = zeros(i*s,i*s,year);
L0 = zeros(i*s,i*s,year);
F0 = zeros(i*s,i*s,year);

% matrix of zeros to be able to create Z0, A0 and F0. See Equation (A11)
boxmask = zeros(ensec,ensec);  

%% EROI EiE Calculation
% these steps are described in detail in Appendix A

% Step 1: Generate basic variables 
for yr = 1:year; %increases in each yr-loop by one year
    x(1,1:i*s,yr) = (sum(Z(:,:,yr),2)+sum(Y(:,:,yr),2))'; % total output vector (see Equation (A1))
    x(x==0) = 0.000000001; % to avoid errors when dividing by x
    A(:,:,yr) = Z(:,:,yr)./repmat(x(:,:,yr),i*s,1); % direct requirements matrix (see Equation (A2))
    L(:,:,yr) = inv(eye(i*s)-A(:,:,yr)); % Leontief matrix (see Equation (A5))
    e(:,:,yr) = F(:,:,yr)./(x(:,:,yr)); % energy intensity coefficient vector (see Equation (A7)) 
    y(:,:,yr)= sum(Y(:,:,yr),2); % total final demand vector (UK + RoW)
    Foot(:,:,yr) = diag(e(:,:,yr))*L(:,:,yr)*diag(y(:,:,yr)); % UK footprint (F in the paper) for each year (see Equation (A10))
        
% Step 2: Generate box varaibles    
    Z0(:,:,yr) = Z(:,:,yr); 
    for rw = 1:s; 
        row = (rw*i)-(i); % increases in each rw-loop by the number of sectors, starting at 0
        Z0(row+c:row+g,c:g,yr) = Z(row+c:row+g,c:g,yr).*boxmask; % Z0 matrix (see Equation (A11))
        A0(:,:,yr) = Z0(:,:,yr)./repmat(x(:,:,yr),i*s,1); % A0 matrix (see Equation (A12))
        L0(:,:,yr) = inv(eye(i*s) - A0(:,:,yr)); % L0 matrix (see Equation (A12))
        F0(:,:,yr) = diag(e(:,:,yr))*L0(:,:,yr)*diag(y(:,:,yr)); % F0 matrix (see Equation (A12))
       
% Step 3: Calculate EiE (see Equation (A13))
        Ei1(rw,:,yr) = Foot(row+c,:,yr) - F0(row+c,:,yr);
        Ei2(rw,:,yr) = Foot(row+g,:,yr) - F0(row+g,:,yr); 
        Ei(rw,yr) = (sum(Ei1(rw,:,yr))+sum(Ei2(rw,:,yr)));  
    end
    EiE(yr) = sum(Ei(:,yr));    
end

