% THIS CODE REPRODUCES FIG 10 IN LAGERLÖF (RED 2006)

clear;

T = 100; % number of periods in simulation

% Numbers from Table  1 and/or 2

gamma = 0.355;

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

tau = 0.28;

alpha = 0.6;

theta = 1;

nstar = 1;


% Endogenous variables being targeted

gstar = 2.362;

estar = 0.075;


% Setting rho and astar to match targets abive; see Footnote 13

rho=-(gstar+2*estar)/(2*tau) + ( ((gstar+2*estar)/(2*tau))^2 - (estar^2-tau*gstar)/(tau^2)  )^.5;

astar = gstar/(estar+rho*tau);
 

Lhat = rho/(theta*(1-rho)); % 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.05*Lhat; % See p. 20

A(1,1) = 0.951;

e(1,1) = 0;

g(1,1) = 0.048;

n(1,1) = 1;

 
% Loop on the dynamical system

for i=1:T
    
    time(i+1)=i;
    
    % Eq. (22)
    z(i) = ((e(i)+rho*tau)/(e(i)+rho*tau+g(i)))^alpha*(A(i)/L(i))^(1-alpha);
    
    % Eq. (19)
    a(i) = min(theta*L(i), astar);
    
    % Eq. (18), or (24)
    g(i+1) = (e(i)+rho*tau)*a(i);
    
    % Eq. (20), or (24)
    e(i+1) = max(0, sqrt(g(i+1)*tau*(1-rho))-rho*tau);
    
    % Eq. (23):
    if z(i)<zbar & z(i)>cbar
        n(i) = (1-cbar/z(i))/(tau+e(i+1));
    elseif z(i)<=cbar       
        n(i) = 0;        
    else
        n(i) = gamma/(tau+e(i+1));
    end

    % Eq. (24)
    L(i+1) = L(i)*n(i);
    
    % Eq. (24)
    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

% Alt way to compute annual gn:
% gn =100*(n.^(1/20)-1);

% Getting all vectors of same size, plus rescaling e

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


% Figure

set(gcf,'Position',[100, 80, 1100, 600]);

plot(time,e,'r:',time,gn,'b-',time,gA,'g.-',time,gz,'k--','linewidth',3);
Legend=legend('Education ($e_{t}$) times 100','Growth rate of $L_{t}$',...
    'Growth rate of $A_{t}$','Growth rate of $z_{t}$');
ylim([-.2,8.2]);
Label=xlabel('Period');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','east');












