clear;

% Scalars
T = 100; % number of years in simulation
Delta = 1/1000; % inverse of number of steps per year
S = T/Delta; % number of steps in simulation
alpha = 1/3; % cap share of output
n = .01;
g = .02;
delta = .03;

% Vector with saving rates
s=zeros(S,1);
s_0=.1;
s_1=.2;
t_hat=S/5;
s(1:t_hat-1,1)=s_0;
s(t_hat:S,1)=s_1;

% Other vectors
time=zeros(S,1);
k=zeros(S,1);

% Steady state levels of k
k_star_0=(s_0/(n+g+delta))^(1/(1-alpha));
k_star_1=(s_1/(n+g+delta))^(1/(1-alpha));

%Initial values

time(1) = 0;
k(1) = k_star_0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Loop on the dynamical system

for t=1:S-1

time(t+1)=T*t/S;
k(t+1)=k(t)+Delta*( s(t)*(k(t)^alpha)-(n+g+delta)*k(t) );

end

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

% Output per effective worker
y=k.^alpha;
y_star_0=k_star_0^alpha;
y_star_1=k_star_1^alpha;

% Productivity
A=exp(g*time);

% Output per worker
x=A.*y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solutions to Problem 12 start here:

% Part (i): pre-shock path of ln(x)
x_preshock=A.*y_star_0;

% Part (ii): cons per eff worker
c=(1-s).*y;

% Part (iii): normalized output per eff worker
y_norm=y/y_star_0;

% Part (iv): line that is tangent to y at t_hat
slope=(s_1-s_0)*alpha*(k_star_0)^(2*alpha-1);
tang_line=y_star_0*ones(S,1)+slope*(time-time(t_hat,1));

% Continued with plots below

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

% Setting current directory (where I want to save my pdf files)
cd C:\Users\lagerlof\Dropbox\OnlineLectures4020\Matlab;

% Plot over income per eff worker
plot(time,y,'b-','linewidth',1.5);
hold on;
plot(time,ones(S,1)*y_star_0,'r--','linewidth',1);
plot(time,ones(S,1)*y_star_1,'g-.','linewidth',1);
hold off;
Legend=legend('$y(t)$','$y_{0}^{*}$','$y_{1}^{*}$');
ylim([.8*y_star_0,1.1*y_star_1]);
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','west');
saveas(gcf, 'Figure_1.pdf');

% Plot over time path for log income per worker
plot(time,log(x),'b-','linewidth',1.5);
Legend=legend('$\ln(x(t))$');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','northwest');
saveas(gcf, 'Figure_2.pdf');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots to solve Problem 12:

% Plot for part (i)
plot(time,log(x),'b-','linewidth',1.5);
hold on;
plot(time,log(x_preshock),'r:','linewidth',.75);
hold off;
Legend=legend('$\ln(x(t))$','pre-shock path');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','northwest');
saveas(gcf, 'Figure_3.pdf');

% Plot for part (ii)
plot(time,c,'b-','linewidth',1.5);
Legend=legend('$c(t)$');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','east');
saveas(gcf, 'Figure_4.pdf');


% Plot for part (iii)
plot(time,y_norm,'b-','linewidth',1.5);
ylim([.8,1.1*max(y_norm)]);
Legend=legend('$\tilde y(t)$');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','east');
saveas(gcf, 'Figure_5.pdf');


% Plot for part (iv)
plot(time,y,'b-','linewidth',1.5);
hold on;
plot(time,ones(S,1)*y_star_0,'r--','linewidth',1);
plot(time,ones(S,1)*y_star_1,'g-.','linewidth',1);
plot(time,tang_line,'k:','linewidth',.75);
hold off;
Legend=legend('$y(t)$','$y_{0}^{*}$','$y_{1}^{*}$','$y_{0}^{*}+\dot y(\hat t)(t-\hat t)$');
ylim([.8*y_star_0,1.1*y_star_1])
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','east');
saveas(gcf, 'Figure_6.pdf');









