close all;
clear all;
pack;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Self-organized criticality in a model of production and inventory dynamics
%% A raw and simple program to numerically simulate the model developed by:
%% P. Bak, K. Chen, J. Scheinkman and M. Woodford (1993),
%% Aggregate Fluctuations from Independent Sectoral Shocks: Self-Organized Criticality
%% in a Model of Production and Inventory Dynamics, Ricerche Economiche, 47 (1), 3-30.
%% Author: Marco Sandri - Verona - Italy - Email: sandri.marco@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters definition
%%%%%%%%%%%%%%%%%%%%%%%%%
N = 500;
n=5;
T = 2^n;
gamma = 5/6;
p = N^(-gamma);
% Declare matrices
X1 = unidrnd(2,N,N)-1;
X2 = zeros(N);
Y = zeros(N);
S = zeros(N);
Ytot = zeros(T,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time evolution of the economy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for t = 1:T
% Exogenous shocks S(t)
S(1,:) = (rand(1,N)< p);
% Calculate the state X(t) of the system
for i = 1:N
if i > 1
S(i,2:N) = (Y(i-1,1:N-1)+Y(i-1,2:N))/2;
end
x = X1(i,:); s = S(i,:);
tmp = 1-((x==0 & s==0) | (x==1 & s==0) | (x==1 & s==1));
y= tmp*2;
X2(i,:) = x+y-s;
Y(i,:) = y;
end
X2(1,:)= S(1,:);
X1=X2;
Ytot(t) = sum(sum(Y(:,2:N)));
disp(t);
end
%%%%%%%%%%%%%%%%%%%%%%
% Draw variuos graphs
%%%%%%%%%%%%%%%%%%%%%%
figure(1);
Ytot = Ytot/N^(3*(1-gamma));
bins=logspace(-2.2,-0.8,100);
[f,x]=hist(Ytot,bins);
sel = find(f>0);
P=polyfit(log10(x(sel(2:end-1))),log10(f(sel(2:end-1))),1);
log10fhat=polyval(P,log10(x(sel(2:end-1))));
% Time evolution
subplot(2,2,1);
plot(Ytot(T/2:T/2+199));
xlabel('Time'); ylabel('Y(t)');
title('Time evolution of Y(t)');
% Spectral analysis of Y(t)
subplot(2,2,2);
Y = fft(Ytot-mean(Ytot),T);
Pyy = Y.* conj(Y) / T;
f = 1000*(0:T/2)/T;
plot(f,Pyy(1:T/2+1));
title('Frequency content of y');
xlabel('frequency (Hz)');
axis tight;
% Histogram
subplot(2,2,3);
hist(Ytot,500);
xlabel('Y'); ylabel('f(Y)');
title('Frequency distribution of aggregate production Y');
axis([0 50 0 1000])
% Log-log plot
subplot(2,2,4);
plot(log10(x(sel(2:end-1))),log10(f(sel(2:end-1))),'.',log10(x(sel(2:end-1))),log10fhat);
%loglog(x(sel(2:end-1)),f(sel(2:end-1)),'.',x(sel(2:end-1)),10.^log10fhat);
xlabel('log_{10}(Y)'); ylabel('log_{10}(f(Y))');
title('Frequency distribution of aggregate production Y');