%% This is a sample code for structural decomposition analysis, taking 5 decomposed factors as an example.

%% INPUT variables
% E_1995 and E_1996: sectoral resource uses or emissions in 1995 and 1996, respectively
% X_1995 and X_1996: sectoral total outputs in 1995 and 1996, respectively
% Z_1995 and Z_1996: inter-sectoral transations matrices in 1995 and 1996, respectively
% FS_1995 and FS_1996: final demand structure matrices indicating each
% sector's share in each final demand category in 1995 and 1996, respectively
% FL_1995 and FL_1996: per capita final demand level in 1995 and 1996, respectively
% P_1995 and P_1996: the population in 1995 and 1996, respectively

%% OUTPUT variables
% vector Q: relative contributions of each decomposed factor

%% codes
n = length(Z_1995); % n is the number of sectors
I = eye(n); % I is an identify matrix

% calculating sectoral emission intensity vectors in 1995 and 1996, respectively
E0 = zeros(n,1);
E1 = zeros(n,1);
for i = 1:n    
    if X_1995(i,1) ~= 0      
        E0(i,1) = E_1995(i,1)/X_1995(i,1);
    else
        E0(i,1) = 0;
    end    
end
for i = 1:n    
    if X_1996(i,1) ~= 0      
        E1(i,1) = E_1996(i,1)/X_1996(i,1);
    else
        E1(i,1) = 0;
    end    
end
E0=(E0)';
E1=(E1)';
DE=E1-E0; % changes in sectoral emission intensity

% calculating direct input coefficient matrices in 1995 and 1996, respectively
A_1995 = zeros(n);
A_1996 = zeros(n);
for i = 1:n
   for j = 1:n
       if X_1995(j,1) ~= 0
           A_1995(i,j) = Z_1995(i,j)/X_1995(j,1);              
       else
           A_1995(i,j) = 0;
       end 
   end
end
for i = 1:n
   for j = 1:n
       if X_1996(j,1) ~= 0
           A_1996(i,j) = Z_1996(i,j)/X_1996(j,1);              
       else
           A_1996(i,j) = 0;
       end 
   end
end

% calculating Leontief inverse matrices in 1995 and 1996, respectively
L0 = inv(I-A_1995);
L1 = inv(I-A_1996);
DL=L1-L0; % changes in Leontief inverse matrices indicating production structure

FS0=FS_1995;
FS1=FS_1996;
DFS=FS1-FS0; % changes in final demand structure

FL0=FL_1995;
FL1=FL_1996;
DFL=FL1-FL0; % changes in per capita final demand level

P0=P_1995;
P1=P_1996;
DP=P1-P0; % changes in the population

% calculating relative contributions of each decomposed factor, Q
QE=1/120*(24*P0*DE*L0*FS0*FL0+6*P0*DE*L0*FS0*FL1+6*P0*DE*L0*FS1*FL0+6*P0*DE*L1*FS0*FL0+6*P1*DE*L0*FS0*FL0+4*P0*DE*L0*FS1*FL1+4*P0*DE*L1*FS0*FL1+4*P1*DE*L0*FS0*FL1+4*P1*DE*L0*FS1*FL0+4*P1*DE*L1*FS0*FL0+4*P0*DE*L1*FS1*FL0+6*P0*DE*L1*FS1*FL1+6*P1*DE*L0*FS1*FL1+6*P1*DE*L1*FS0*FL1+6*P1*DE*L1*FS1*FL0+24*P1*DE*L1*FS1*FL1);
QL=1/120*(24*P0*E0*DL*FS0*FL0+6*P0*E0*DL*FS0*FL1+6*P0*E0*DL*FS1*FL0+6*P0*E1*DL*FS0*FL0+6*P1*E0*DL*FS0*FL0+4*P0*E0*DL*FS1*FL1+4*P0*E1*DL*FS0*FL1+4*P1*E0*DL*FS0*FL1+4*P1*E0*DL*FS1*FL0+4*P1*E1*DL*FS0*FL0+4*P0*E1*DL*FS1*FL0+6*P0*E1*DL*FS1*FL1+6*P1*E0*DL*FS1*FL1+6*P1*E1*DL*FS0*FL1+6*P1*E1*DL*FS1*FL0+24*P1*E1*DL*FS1*FL1);
QFS=1/120*(24*P0*E0*L0*DFS*FL0+6*P0*E0*L0*DFS*FL1+6*P0*E0*L1*DFS*FL0+6*P0*E1*L0*DFS*FL0+6*P1*E0*L0*DFS*FL0+4*P0*E0*L1*DFS*FL1+4*P0*E1*L0*DFS*FL1+4*P1*E0*L0*DFS*FL1+4*P1*E0*L1*DFS*FL0+4*P1*E1*L0*DFS*FL0+4*P0*E1*L1*DFS*FL0+6*P0*E1*L1*DFS*FL1+6*P1*E0*L1*DFS*FL1+6*P1*E1*L0*DFS*FL1+6*P1*E1*L1*DFS*FL0+24*P1*E1*L1*DFS*FL1);
QFL=1/120*(24*P0*E0*L0*FS0*DFL+6*P0*E0*L0*FS1*DFL+6*P0*E0*L1*FS0*DFL+6*P0*E1*L0*FS0*DFL+6*P1*E0*L0*FS0*DFL+4*P0*E0*L1*FS1*DFL+4*P0*E1*L0*FS1*DFL+4*P1*E0*L0*FS1*DFL+4*P1*E0*L1*FS0*DFL+4*P1*E1*L0*FS0*DFL+4*P0*E1*L1*FS0*DFL+6*P0*E1*L1*FS1*DFL+6*P1*E0*L1*FS1*DFL+6*P1*E1*L0*FS1*DFL+6*P1*E1*L1*FS0*DFL+24*P1*E1*L1*FS1*DFL);
QP=1/120*(24*DP*E0*L0*FS0*FL0+6*DP*E0*L0*FS0*FL1+6*DP*E0*L0*FS1*FL0+6*DP*E0*L1*FS0*FL0+6*DP*E1*L0*FS0*FL0+4*DP*E0*L0*FS1*FL1+4*DP*E0*L1*FS0*FL1+4*DP*E1*L0*FS0*FL1+4*DP*E1*L0*FS1*FL0+4*DP*E1*L1*FS0*FL0+4*DP*E0*L1*FS1*FL0+6*DP*E0*L1*FS1*FL1+6*DP*E1*L0*FS1*FL1+6*DP*E1*L1*FS0*FL1+6*DP*E1*L1*FS1*FL0+24*DP*E1*L1*FS1*FL1);
Q=[QE;QL;QFS;QFL;QP];
