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');