clear;

len=500;

%exogenous parameters

cbar=1;
alpha=.88;
gama=.72;
D=4.6;
theta=0.5;
beta=0.18;
q=2.7;
Abar=17;

% Checking that Assumption 2 holds

must_be_negative_2=beta*cbar/[(1-beta)*q]-1-alpha

% Checking that Assumption 1 holds

must_be_negative_1=alpha*(1+gama)^(1-alpha)-1

numm=cbar*(1+gama)^(1-alpha);
denn=1-alpha*(1+gama)^(1-alpha);
psiconst=(numm/denn)^(1/alpha);

%matrices

PX=ones(len,1);
PXs1=ones(len,1);
PXs2=ones(len,1);
PXs3=ones(len,1);

PsiX=ones(len,1);
PhiX=ones(len,1);
OmegaX=ones(len,1);

ALX=ones(len,1);

PLXfree=ones(len,1);
PLXslave=ones(len,1);
PLXslave2=ones(len,1);

%functions

Psi=@(P) psiconst*P;

OmegaNum=@(P) cbar*((1+gama)^(1-alpha))*(P.^(1+alpha));
OmegaDen=@(P) P-(1+gama)^(1-alpha);

Omega=@(P) (OmegaNum(P)/OmegaDen(P))^(1/alpha);

Phi=@(P) (alpha)^(-1/(1-alpha))*(cbar/(1-alpha))^(1/alpha)*P^(-alpha/(1-alpha));

AL=@(P) Abar+(D^(1/theta))*P;

PLfree=@(P) ((q/beta)^(1/alpha))*P;

PLslave=@(P) ((1+gama)^((1-alpha)/alpha))*((q/beta)+cbar)^(1/alpha)*P;

% Drawing the loci and the institutional borders in the phase diagram

for Pn=1:len;

% Adjusting the scales

PX(Pn)=0+11*(Pn-1)/len;

PXs1(Pn)=(1/alpha)+10*(Pn-1)/len;

PXs2(Pn)=1/alpha-.06*(Pn-1)/len;

PXs3(Pn)=1/alpha+5*(Pn-1)/len;

% The institutional borders

PsiX(Pn)=Psi(PXs1(Pn));

PhiX(Pn)=Phi(PXs1(Pn));

OmegaX(Pn)=Omega(PXs2(Pn));

ALX(Pn)=AL(PX(Pn));

PLXfree(Pn)=PLfree(PX(Pn));

% The loci

PLXslave(Pn)=PLslave(PX(Pn));
PLXslave2(Pn)=PLslave(PXs3(Pn));

end;

%Path

hor=100;

% Functions for updated pop and technology

PnewFree=@(A,P) (beta/q)*(A^alpha)*(P^(1-alpha));

PnewSlave=@(A,P) (beta/q)*(A^alpha*(P/(1+gama))^(1-alpha)-cbar*P);

Anew=@(A,P) Abar+D*(A-Abar)^(1-theta)*(P^theta);

%initial values

A=Abar*1.1;
P=0.5;
Apath(1)=A;
Ppath(1)=P;

% The "for" algorithm to create the paths

for t=2:hor;

Time(t)=t-1;

if (A > Omega(P)) & (A > Psi(P)) & (P>(1+gama)^(1-alpha));
% Updating pop under slavery
Ppath(t)=PnewSlave(A,P);
Switch1=t-1.5;
else;
% Updating pop under free labor or egalitarianism
Ppath(t)=PnewFree(A,P);
end

if (P<1/alpha) | (A<Phi(P));
Switch2=t-1.5;
end;

%Updating technology
Apath(t)=Anew(A,P);

w(t)=(1-alpha)*(A/P)^alpha;

P=Ppath(t);
A=Apath(t);

end;

% Creating vertical and horizontal lines for the diagram

hort=(1/alpha)*ones(2,1);
vert=ones(2,1);
vert(1)=Phi(1/alpha);
vert(2)=Psi(1/alpha);

% Creating 2x2 plots

%First plot: institutional borders and the path

subplot(2,2,1);
plot(PXs1,PsiX,'b',PXs2,OmegaX,'b',PXs1,PhiX,'b',...
hort,vert,'b',Ppath,Apath,'ko-')
xlim([0 5])
ylim([0 100])
str1 = {'Slavery'};
str2 = {'Egalitarianism'};
str3 = {'Free labor'};
text(1.35,70,str1);
text(0.2,5,str2);
text(2.5,30,str3);
title('Institutional regions and the path');
xlabel('\itP_{t}');
ylabel('\itA_{t}');

%Second plot: the loci and the path

subplot(2,2,2);
plot(PX,ALX,'b:',PX,PLXfree,'--',PX,PLXslave,'r',Ppath,Apath,'ko-')
xlim([0 5])
ylim([0 100])
title('Loci and the path');
legend('{\bf\itL^{A}(P)}','{\bf\itL^{E/F}(P)}','{\bf\itL^S(P)}','Path')
xlabel('\itP_{t}');
ylabel('\itA_{t}');

% Third plot: Loci, institutional regions, and the path

subplot(2,2,3);
plot(PXs1,PsiX,'b',PXs2,OmegaX,'b',PXs1,PhiX,'b',...
hort,vert,'b',PX,ALX,'b:',PX,PLXfree,'--',PX,PLXslave,'r',Ppath,Apath,'ko-')
xlim([0 5])
ylim([0 100])
title('Loci, institutional regions, and the path');
xlabel('\itP_{t}');
ylabel('\itA_{t}');

% Indexing the paths for P and A

Ppathind=Ppath./Ppath(1);
Apathind=Apath./Apath(1);

% Fourth plot: the paths of A, P over time

subplot(2,2,4);
plot(Time,Ppathind,'k*-',Time,Apathind,'ko-')
xlim([0 30])
ylim([0 10])
title('The paths over time (initial level=1)');
line([Switch2,Switch2],[0,10])
line([Switch1,Switch1],[0,10])
xlabel('Period');
legend('Population','Technology')

str1 = {'Slavery'};
str2 = {'Egalitarianism'};
str3 = {'Free labor'};
text(0,8,str1);
text(0.2,5,str2);
text(8,1,str3);

%Test=(1-beta)*w';

%Checking that average product > subsistence
% (only needs to hold in egal regime)
must_be_positive_1=(Apath'./Ppath').^(alpha)-cbar
%Checking that wage > subsistence
% (only needs to hold in free labor regime)
must_be_positive_2=(1-alpha)*(Apath'./Ppath').^alpha-cbar