clear 
clc
%single intervention
n=10;
n_max=10000;
while n<n_max;
%random matrix generation for technical coefficient matrix
%A with size n * n
c = 0;
d = 0.001; %max of the tech coefficients
A = (d-c).*rand(n,n) + c;
b=rand(1,n);
y=rand(n,1 );
n
L= (eye(n)-A);

%The initial Leontief inverse
inv_L=inv(L); 
%select a single sector k which will be tested 
%for the circular intervention
k=3;
%select the sector j whose input to sector k to be decreased, 
%and the shift amount c_j
j= 1;
c_j=0.5;
%select the sector i whose input to sector k to be increased, 
%and the shift amount c_i
i= 4;
c_i=0.3;

%modified  L after the interventions
modified_L=L;
modified_L(i,k)= L(i,k) - c_i*A(j,k); 
modified_L(j,k)= L(j,k) + c_j*A(j,k); 

%emissions before the intervention
m=b*inv_L*y;

%change in the emissions after the intervention 
%(default calculation)
tic
delta_m_default=b * inv(modified_L) *y-m;
time_d= toc;

%change in the emissions after the intervention using
%Sherman- Morrison formula
tic
delta_m_SM = (A(j,k)*(c_i* b * inv_L(:,i ) - c_j* b* inv_L(:,j) ) * inv_L(k,:)*y) / (1-A(j,k) * (c_i*inv_L(k,i)-c_j * inv_L(k,j)));
time_w=toc;

%store the ratios of the calculation times
ratio(n)=time_d/time_w;
n=n+10;
end

%plotting the performance plots
figure;
scatter(10:10:9990,ratio(10:10:9990))
xlabel('n, the size of the Leontief inverse')
ylabel('Ratio of calculation times')