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;
rho=.02;
theta=1;

% Vectors
time=zeros(S,1);
k=zeros(S,1);
c=zeros(S,1);

% Steady state levels of k, c
k_star=(alpha/(rho+theta*g))^(1/(1-alpha));
c_star=k_star^alpha-(n+g)*k_star;

%Initial values

time(1) = 0;
k(1)=.25*k_star;

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

% First guess given in template file; too high
c(1)=.5*c_star;

% Try unsilencing the commands below, one by one

% Lower guess; too low
% c(1)=.4*c_star;

% Better guess; looks okay (for T=100)
% c(1)=.45*c_star;

% Should work for T=200:
% c(1)=.45001316415*c_star;

% See further code at the bottom

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

for t=1:S-1

time(t+1)=T*t/S;

k(t+1)=max(0,k(t)+Delta*( k(t)^alpha-c(t)-(n+g)*k(t) ));

c(t+1)=c(t)+(c(t)*Delta/theta)*(alpha*k(t)^(alpha-1) - rho-theta*g);
if k(t+1)==0 
    c(t+1)=0;
end

end

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

% Output per effective worker
y=k.^alpha;
y_star_0=k_star^alpha;

% Productivity
A=exp(g*time);

% Output per worker
x=A.*y;

% Normalized output per eff worker
y_norm=y/y_star_0;


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

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

% Time path for cap per eff worker
plot(time,k,'b-','linewidth',1.5);
hold on;
plot(time,ones(S,1)*k_star,'r--','linewidth',1);
hold off;
Legend=legend('$k(t)$','$k^{*}$');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','northwest');
saveas(gcf, 'Figure_R1.pdf');


% Time path for cons per eff worker
plot(time,c,'b-','linewidth',1.5);
hold on;
plot(time,ones(S,1)*c_star,'r--','linewidth',1);
hold off;
Legend=legend('$c(t)$','$c^{*}$');
Label=xlabel('Year');
set(Label,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','northwest');
saveas(gcf, 'Figure_R2.pdf');


% Creating some vectors to draw loci in phase diagram
k_max=(1/(n+g))^(1/(1-alpha)); % Point where (k-dot=0)-locus intersects k axis
k_span=linspace(0,.5*k_max,S)'; % Column vector with S elements going from 0 to halfway towards k_max
c_span=linspace(0,max(c),S)'; % Column vector with S elements going from 0 to max value of c in simulation above
k_dot_0_locus=k_span.^alpha - (n+g)*k_span; % Column vector to plot (k-dot=0)-locus
c_dot_0_locus=k_star*ones(S,1); % Column vector to plot (k-dot=0)-locus

% Phase digram
plot(k,c,'-b*','MarkerSize',4,'linewidth',.5);
hold on;
plot(k_span,k_dot_0_locus,'r:','linewidth',1.5);
plot(c_dot_0_locus,c_span,'g-','linewidth',1.5); % Reverse order since c is on vertical axis
hold off;
Legend=legend('Trajectory','$\dot k=0$','$\dot c=0$');
Label_x=xlabel('k');
Label_y=ylabel('c');
set(Label_x,'Interpreter','latex','FontSize',15);
set(Label_y,'Interpreter','latex','FontSize',15);
set(Legend,'Interpreter','latex','FontSize',15,'location','east');
saveas(gcf, 'Figure_R3.pdf');

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

% To check if approching steady state through full simulation:
% see if gaps stay negative

gap_k=k-k_star;
gap_c=c-c_star;

% Check values of last 5 elements in gap vectors
% Skip ; to make result show in command window
test_k=gap_k(S-4:S)
test_c=gap_c(S-4:S)





