clear;

clc;

 % To Graph Dr. Lagerlof Figure 5

T = 100; % number of periods in simulation



% Setting parameters -- numbers from Table  1 and/or 2

gamma = 0.255;

cbar = 1; %this is ctilda in GW+Lagerlof

tou = 0.15; % From Table 1

alpha = 0.6;

theta = 1;

nstar = 1;

% Endogenous variables

gstar = 2.362;

estar = 0.075;

% Computing ro and astar as described in Footnote 13

ro=-(gstar+2*estar)/(2*tou) + ( ((gstar+2*estar)/(2*tou))^2 - (estar^2-tou*gstar)/(tou^2)  )^.5;

astar = gstar/(estar+ro*tou);


Lhat = ro/(theta*(1-ro)); % See eq. (21)

zbar = cbar/(1-gamma); %this is ztilda in GW+Lagerlof


% Matrices

time=zeros(T+1,1);

L=zeros(T,1);

n=zeros(T,1);

A=zeros(T,1);

z=zeros(T+1,1);

g=zeros(T+1,1);

e=zeros(T,1);

a=zeros(T,1);

gz=zeros(T+1,1);

gA=zeros(T+1,1);

gn=zeros(T+1,1);



%Initial values; see Table 2

L(1,1) = 0.364; % See p. 20

A(1,1) = 0.87;

e(1,1) = 0;

g(1,1) = 0.048;

n(1,1) = 1;





%Dynamic System

for i=1:T

    time(i+1)=i;

    z(i) = ((e(i)+ro*tou)/(e(i)+ro*tou+g(i)))^alpha*(A(i)/L(i))^(1-alpha);

    a(i) = min(theta*L(i), astar);

    g(i+1) = (e(i)+ro*tou)*a(i);

    e(i+1) = max(0, sqrt(g(i+1)*tou*(1-ro))-ro*tou);

    if z(i)<zbar

        n(i) = (1-cbar/z(i))/(tou+e(i+1));

    else

        n(i) = gamma/(tou+e(i+1));

    end

    L(i+1) = L(i)*n(i);

    A(i+1) = A(i)*(1+g(i+1));

end



for j=1:T

    % To make rates annual; each period=20 years
    gz(j+1) =100*( [1+(z(j+1)-z(j))/z(j)]^(1/20)-1);
    gA(j+1) = 100*([1+(A(j+1)-A(j))/A(j)]^(1/20)-1);
    gn(j+1) =100*( [n(j)]^(1/20)-1);

end

% Getting all vectors of same size

e=100*e(2:100);
gn=gn(1:99);
gA=gA(2:100);
gz=gz(2:100);
time=time(1:99);


% Figure

plot(time,e,':',time,gn,'-',time,gA,'.-',time,gz,'--');
h2=legend('education ($e_{t}$) times 100','growth rate of $L_{t}$','growth rate of $A_{t}$','growth rate of $z_{t}$');
h3=title('Figure5');
ylim([-.2,8.2]);
xlabel('Period (t)');
set(h2,'Interpreter','latex');
% set(h3,'Interpreter','latex');