clear;

cd 'C:\Users\lagerlof\Dropbox\MultipleStStHood\Replication_February_2021\'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameters
alpha=.5;
sigma=2;
rho=(alpha/sigma)+1-alpha;
beta=.95;
phi_init=.01; % initial level of phi
eta=1;
z_unbar=.01;
z_upbar=.99;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hor=100; % Number of periods for simulation
soc=200; % Number of societies (or countries)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kappa=(  (  beta*sigma*(1-alpha) +  sigma  + alpha*beta  )  /...
    ( sigma *(1-alpha*beta) + alpha*beta )  )^rho;

omega = beta*sigma*(1-alpha) / ( beta*sigma*(1-alpha)  + sigma+ alpha*beta );
% omega is shorthand to calc tau when constr. not binding

tau_unbar=1- (beta*sigma*(1-alpha) / ( beta*sigma*(1-alpha)  + sigma*(1-beta)+ alpha*beta ));

gama = (100*(( eta*sigma*(1-alpha)/alpha )^(alpha/sigma)) * ((phi_init/omega)^rho ) )^(1/(1-alpha));
% Set gamma to make initial D=100

D_init=(( alpha/( eta*sigma*(1-alpha) ) )^(alpha/sigma) )*(gama^(1-alpha)) * ( (omega/phi_init)^rho );
E=   sigma*beta   /   (  sigma*beta*(1-alpha)   + sigma  +  alpha*beta   );
% E is shorthand to calc z(t+1) when constr. not binding

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculating thresholds for B (given phi=phi_init, and associated levels of X_unbar, X_upbar)
X_unbar=(z_unbar/phi_init)*( sigma + (1-sigma)*alpha*beta )/(beta*sigma);
X_upbar=(1/phi_init)*( z_upbar*( sigma*beta*(1-alpha)  + sigma + alpha*beta )/(beta*sigma)  - z_unbar );
B_hat=(1/z_unbar)*((X_unbar^(1-rho)  /   ( kappa*D_init* (phi_init^rho) )  )^(1/alpha) );
B_2hat=(1/z_upbar)*((X_upbar  /   ( kappa*D_init* (phi_init*X_upbar+z_unbar-z_upbar)^rho)  )^(1/alpha) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Vectors
Time_x=zeros(hor,1);
z_x=zeros(hor,soc);
Y_x=zeros(hor,soc);
tau_x=zeros(hor,soc);
z_mean=z_unbar*ones(hor,1);
phi_x=zeros(hor,1);

ind=zeros(hor,1); % Indicator for which constraint binds

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Levels of B for all soc societies
weight_x=0:1:soc-1; % Vector with soc elements
weight_x=weight_x/(soc-1); % Norm to fall btw 0 and 1
B_x=(1-weight_x)*B_2hat+weight_x*B_hat;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Start values for vectors

Y_unbar_x=(kappa*D_init*(B_x.^alpha)*(z_unbar^(alpha-1))*((phi_init*z_unbar)^rho)).^(1/(1-rho));
z_x(1,:)=z_unbar;
Y_x(1,:)=Y_unbar_x(1,:);
tau_x(1,:)=tau_unbar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The loop

for t=1:1:hor-1
    
% Things that change over time
Time_x(t+1,1)=t; % time runs from 0

% Shock to extractive capacity
if t==41
z_x(t,1:10:end)=z_upbar;
end

% Setting phi
z_mean(t,1)=mean(z_x(t,:),2);
w_n=(z_mean(t,1)-z_unbar)/(z_upbar-z_unbar); % weight
phi=(1-w_n)*phi_init+w_n*phi_init*30;
phi_x(t,1)=phi;

% Updating things that depend on phi (but not B)
D=(( alpha/( eta*sigma*(1-alpha) ) )^(alpha/sigma) )*(gama^(1-alpha)) * ( (omega/phi)^rho );
X_unbar=(z_unbar/phi)*( sigma + (1-sigma)*alpha*beta )/(beta*sigma);
X_upbar=(1/phi)*( z_upbar*( sigma*beta*(1-alpha)  + sigma + alpha*beta )/(beta*sigma)  - z_unbar );
    
    
for i=1:1:soc
    
% Things that vary across societies
    
% B = land prod for soc i
B=B_x(1,i);

% If no constraint binds
Y_x(t+1,i)=D*(B^alpha)*(z_x(t,i)^(alpha-1))*(phi*Y_x(t,i)*z_x(t,i)+z_unbar)^rho;
z_x(t+1,i)=E*(phi*Y_x(t,i)*z_x(t,i)+z_unbar);
tau_x(t,i)=(1-omega)-omega*(z_unbar/(phi*z_x(t,i)*Y_x(t,i)));
 
 
% If one of the constraints binds
if z_x(t,i)*Y_x(t,i)<X_unbar
 Y_x(t+1,i)=kappa*D*(B^alpha)*(z_x(t,i)^(alpha-1))*((phi*z_x(t,i)*Y_x(t,i))^rho);
 z_x(t+1,i)=z_unbar;
 tau_x(t,i)=tau_unbar;
 ind(t,1)=1;
elseif z_x(t,i)*Y_x(t,i)>X_upbar    
 Y_x(t+1,i)=kappa*D*(B^alpha)*(z_x(t,i)^(alpha-1))*(phi*Y_x(t,i)*z_x(t,i)+z_unbar-z_upbar)^rho;
 z_x(t+1,i)=z_upbar;
 tau_x(t,i)=(tau_unbar)+(1-tau_unbar)*(z_upbar-z_unbar)/((phi*z_x(t,i)*Y_x(t,i)));
 ind(t,1)=2;
end
 
 
end    
end
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fixing vectors for last period (t=hor)
for i=1:1:soc
tau_x(hor,i)=(1-omega)-omega*(z_unbar/(phi_x(hor,1)*z_x(hor,i)*Y_x(hor,i)));
if z_x(hor,i)*Y_x(hor,i)<X_unbar
tau_x(hor,i)=tau_unbar;
elseif z_x(hor,i)*Y_x(hor,i)>X_upbar    
tau_x(hor,i)=tau_unbar+(1-tau_unbar)*(z_upbar-z_unbar)/((phi_x(hor,1)*z_x(hor,i)*Y_x(hor,i)));
end
end

z_mean(hor,1)=mean(z_x(hor,:),2);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Figure 2

Line_width=1.5;

plot(Time_x,log(z_x(:,2)),'b-*','LineWidth',Line_width)
hold on;
plot(Time_x,log(z_x(:,soc/2)),'r-o','LineWidth',Line_width)
plot(Time_x,log(z_x(:,soc)),'g-+','LineWidth',Line_width)
plot(Time_x,log(z_mean),'k:','LineWidth',2)
hold off;

A1subh1=legend('Low-$B$ society','Intermediate-$B$ society','High-$B$ society','Log of mean $z_{t}$ across societies','Location','northwest');
legend('boxoff');
ylabel=ylabel('Log of extractive capacity, ln($z_{t}$)');
xlabel=xlabel('Time ($t$)');
hold off;


set(ylabel,'Interpreter','latex','Fontsize',14);
set(xlabel,'Interpreter','latex','Fontsize',14);
set(A1subh1,'Interpreter','latex','Fontsize',14);
set(gcf, 'units','normalized','outerposition',[0 0 .8 .99]);
set(gcf,'PaperOrientation','landscape');
set(gca,'TickLabelInterpreter','latex','Fontsize',14)

print('TablesFigs\Figure_2','-dpdf');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Figure 3

state_hist=mean(z_x,1);
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
clear xlabel ylabel;
s=scatter(B_x,state_hist,'o','filled');
s.LineWidth = 0.6;
s.MarkerEdgeColor = 'b';
s.MarkerFaceColor = 'r';
ylabel_2=ylabel('Mean extractive capacity over 100 periods');
xlabel_2=xlabel('Land productivity, B');
set(ylabel_2,'Interpreter','latex','Fontsize',14);
set(xlabel_2,'Interpreter','latex','Fontsize',14);
set(gcf, 'units','normalized','outerposition',[0 0 .8 .99]);
set(gcf,'PaperOrientation','landscape');
set(gca,'TickLabelInterpreter','latex','Fontsize',14);
box on;

print('TablesFigs\Figure_3','-dpdf');


