% growth_pop_decline.m
% matlab file for the simulation of transitional dynamics in
% Strulik (2013): Long-Run Economic Growth Despite Population Decline
% the code produces Figure 1 and 2 of the paper

clear all; clc

% select iso-elastic (iso=1) or logistic (iso=0) growth rate of ideas
iso=0

% select initial initial value of human capital and line color in Figure
init_val=3

h(1)=1;

if init_val==2
h(1)=0.1;
end

if init_val==3
h(1)=10;
end
% **********************************************************

% ****** calibration and steady state **********************+

% parameters
sig=1;
beta=2;
lam=0.75;
tet=0.8
nu=1.79

eta=0.07
alpha=0.39
gamma=0.26
phi=0.137


%nu=10
mu=30; %length of a generation
gA_max=(1+0.0038).^(mu)-1;
omega=100;


% household decisions:

s_star=(eta*gamma*phi)/(alpha-gamma*(eta+tet));
n_star=(alpha-gamma*(eta+tet))/(1+alpha)/phi;
q_star=gamma*tet*phi/(alpha-gamma*(eta+tet));

% implied steady state growth rates

g_h=(1/n_star)^(lam*sig*tet/(lam*sig*tet-beta*(1-tet)))-1
g_A=(1+g_h)^((1-tet)/sig/tet)-1
g_y=(1+g_A)^sig*(1+g_h)-1

% generational growth rates converted to annual ones
g_A=(1+g_A)^(1/mu)-1
g_h=(1+g_h)^(1/mu)-1
g_y_star=(1+g_y)^(1/mu)-1;
g_L=n_star^(1/mu)-1


%************ adjustment dynamics *****************

t(1)=1;
L(1)=1;
A(1)=100;
w(1)=A(1)^sig;

for i=1:1:160
w(i)=A(i)^sig;
H(i)=(1-phi*n_star-s_star*n_star)*h(i)*L(i); 
q(i)=q_star*w(i)*h(i);
h(i+1)=nu*q(i)^tet*s_star^eta;
L(i+1)=n_star*L(i);

if iso==1
gA(i)=max(0,A(i)^(-beta)*H(i)^lam);
else
gA(i)=gA_max./(1+omega.*exp((-H(i)/A(i))));
end

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


t(i+1)=t(i)+1;

g_A(i+1)=(A(i+1)-A(i))/A(i);
g_L(i+1)=(L(i+1)-L(i))/L(i);
g_h(i+1)=(h(i+1)-h(i))/h(i);

end

h(end)=[];L(end)=[];A(end)=[];t(end)=[];g_A(end)=[];g_L(end)=[];g_h(end)=[];
g_y=(1+g_A).^sig.*(1+g_h)-1;
g_A=(1+g_A).^(1/mu)-1;
g_L=(1+g_L).^(1/mu)-1;
g_h=(1+g_h).^(1/mu)-1;
g_y=(1+g_y).^(1/mu)-1;

t=t*mu+1700;

from=1900;
to=3000 %max(t);

if init_val==1
    close
end

set(gcf,'PaperUnits', 'centimeters','PaperType','a4letter','Units','centimeters','PaperPosition', [00.5 1.5 22 6],'Position',[0.5 1.5 22 6],'Color','white');

subplot(1,2,1)
if init_val==1
plot(t,g_A,'b-','LineWidth',1)
hold on
end
if init_val==2
plot(t,g_A,'r--','LineWidth',1.5)
end
if init_val==3
plot(t,g_A,'g-.','LineWidth',1.5', 'Color',[0 0.5 0])
end
axis([from to 0 0.01])
%axis([from to 0.002 0.003])
set(gca,'YTick',[0 0.002 0.004 0.006 0.008 0.01 0.012])
%set(gca,'YTickLabel',[0 0.002 0.004 0.006 0.008 0.010 0.012])
set(gca,'FontSize',10);
ylabel('growth ideas (g_A)')
xlabel('year')
grid on
box on
hold on


subplot(1,2,2)
if init_val==1
plot(t,g_h,'b-','LineWidth',1)
hold on
end
if init_val==2
plot(t,g_h,'r--','LineWidth',1.5)
end
if init_val==3
plot(t,g_h,'g-.','LineWidth',1.5', 'Color',[0 0.5 0])
end
axis([from to 0.0 0.04])
%axis([from to 0.01 0.02])
%set(gca,'YTickLabel',[0 0.002 0.004 0.006 0.008 0.010 0.012])
set(gca,'FontSize',10);
ylabel('growth human capital(g_h)')
xlabel('year')
grid on
box on
hold on

