% Matlab program for the deterministic continuous-time infinite-horizon Ramsey-Cass-Koopmans growth model.
% Patrick Toche©  created July 2002, modified 22/04/2003
% Interactive parameter selection.

function ramsey1

clear all;
disp('The deterministic continuous-time infinite-horizon neoclassical growth model with CES utility and CES production. Patrick Toche©.')

global A alpha beta g d n s b p;

parameters=input('Baseline Model Parameters? Yes=1, No=0  ');

if parameters==1;
     A=1;
     alpha=1/3;
     beta=1.5;
     g=0.01;
     d=0.1;
     n=0;
     s=1;
     b=0.1;
     p=0.9;
else

A=input('Scale Parameter:  ');
if A<0
     error('Invalid Scale Parameter')
end

alpha=input('Capital Income Share:  ');
if alpha<0
     error('Invalid Capital Income Share')
end
if alpha>1
     error('Invalid Capital Income Share')
end

warning('For an elasticity of substitution between capital and labour close to 1/2, convergence will obtain for horizons T=10 or less, while for beta close to 2, it will require T=500 or more')
beta=input('Elasticity of Substitution between Capital and Labour:  ');
if beta<0
     error('Invalid Elasticity')
end

g=input('Labour-Augmenting Growth Rate:  ');
if g<0
    error('Invalid Growth Rate')
end
if g>1
    error('Invalid Growth Rate')
end

d=input('Depreciation Rate:  ');
if d<0
    error('Invalid Depreciation Rate')
end
if d>1
    error('Invalid Depreciation Rate')
end

s=input('Elasticity of Intertemporal Substitution (1 is log):  ');
if s<0
    error('Invalid Elasticity')
end

b=input('Discount Rate (Time Preference):  ');
if b<0
    error('Invalid Discount Rate')
end
if b>1
    error('Invalid Discount Rate')
end

n=input('Population Growth Rate:  ');
if n<0
    error('Invalid Population Growth Rate')
end
if n>1
    error('Invalid Population Growth Rate')
end

if b<n+(1-1/s)*g
    error('Transversality Condition Violated')
end

p=input('Percentage Deviation from Steady State (e.g. p=0.5):  ');
if p<0
    error('Invalid Percentage Deviation')
end

end


global T0 T N NMax RelTol AbTol;

     T0=10;
     T=T0-1+input('Horizon Length:  ');
     N=10;
     NMax=50;
     RelTol=1e-6;
     AbsTol=1e-3;

global c0 css k0 kss s0 sss;

tic;
fprintf('Parameters    A = %7.5f, alpha = %7.5f, beta = %7.5f, g = %7.5f, d = %7.5f, n = %7.5f, s = %7.5f, b = %7.5f\n',A,alpha,beta,g,d,n,s,b)

%%% Compute the initial steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kss=feval(@z,d+b+g/s); css=feval(@Y,kss)-(n+d+g)*kss;
ess=[css;kss]; sss=feval(@S,css,kss);
fprintf('Steady State css = %7.5f, kss =%7.5f, sss =%7.5f\n',css,kss,sss)

%%%  Computing eigenvalues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
odess=feval(@ODE,1,ess); %fprintf('odess = %7.5f\n',odess);
jss=zeros(2,2); jss=feval(@J,1,ess); %fprintf('jss = %7.5f\n',jss);
Jss=eig(jss); disp('The eigenvalues are:'); disp(Jss);
tr=trace(jss); dt=det(jss); fprintf('tr = %7.9f, dt = %7.9f\n',tr,dt);
if T<T0+1
    error('The horizon is too short to compute the dynamics')
end

%%% Define new parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c0=css; k0=p*kss; s0=sss;

%%% Computing the dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
guess=[c0 k0];
options = bvpset('FJacobian',@J,'Vectorized','on','Stats','off','RelTol',RelTol,'AbsTol',AbsTol','NMax',NMax);
solinit = bvpinit(linspace(0,T0,N),guess);
sol = bvp4c(@ODE,@BC,solinit,options);
t = linspace(0,T0);
e = deval(sol,t);
e(3,:) = feval(@S,e(1,:),e(2,:));
cTss=e(1,end)/css; kTss=e(2,end)/kss;
fprintf('With T = %g, c(T) = %7.5f, k(T) = %7.5f, c(T)/css = %7.5f, k(T)/kss = %7.5f\n',T0,e(1,end),e(2,end),cTss,kTss)

%%% Drawing figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);
clf;
plot(t,e(1,:),'-',t(end),e(1,end),'o');
title('Consumption');
xlabel('t');
ylabel('c(t)');
hold on;
drawnow;

figure(2);
clf;
plot(t,e(2,:),'-',t(end),e(2,end),'o');
title('Capital');
xlabel('t');
ylabel('k(t)');
hold on;
drawnow;

%%% Adjusting end point and udating existing figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for TT = T0+1:T

  solinit = bvpinit(sol,[0 TT]);
  sol = bvp4c(@ODE,@BC,solinit,options);
  t = linspace(0,TT);
  e = deval(sol,t);
  e(3,:) = feval(@S,e(1,:),e(2,:));

  cTss=e(1,end)/css; kTss=e(2,end)/kss;
  fprintf('With T = %g, c(T) = %7.5f, k(T) = %7.5f, c(T)/css = %7.5f, k(T)/kss = %7.5f\n',TT,e(1,end),e(2,end),cTss,kTss)


  figure(1);
  plot(t,e(1,:),'-',t(end),e(1,end),'o');
  figure(2);
  plot(t,e(2,:),'-',t(end),e(2,end),'o');
  drawnow;

end

hold off;

diary off;

save('ramsey.mat');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ode = ODE(t,e)
global g d n s b;
ode = [(feval(@Yk,e(2,:))-d-b-g/s).*s.*e(1,:);
       feval(@Y,e(2,:))-(n+d+g).*e(2,:)-e(1,:)];

function bc = BC(e0,eT)
global css kss c0 k0;
bc = zeros(2,1);
bc(1) = e0(2)-k0;
bc(2) = eT(2)-kss;

function j = J(t,e)
global g d n s b;
j = zeros(2,2);
j(1,1) = 0;
j(1,2) = feval(@Ykk,e(2,:)).*s.*e(1,:);
j(2,1) = -1;
j(2,2) = feval(@Yk,e(2,:))-(n+d+g);

function rr = Y(k) %Production
global A alpha beta;
if beta==1
rr = A*k.^(alpha);  %Cobb-Douglas Production
else
rr = A*(alpha*k.^((beta-1)/beta)+1-alpha).^(beta/(beta-1));  %CES Production
end

function rr = Yk(k) %Marginal Product
global A alpha beta;
if beta==1
rr = A*alpha*k.^(alpha-1);
else
rr = A*alpha*(alpha +(1-alpha)*k.^(-(beta-1)/beta)).^(1/(beta-1));
end

function rr = Ykk(k) %Derivative of the Marginal Product
global A alpha beta;
if beta==1
rr = A*alpha*(alpha-1)*k.^(alpha-2);
else
rr = A*alpha/beta*(alpha-1).*k^(-(2*beta-1)/beta)*(alpha+(1-alpha).*k^(-(beta-1)/beta))^(-(-2+beta)/(beta-1));
end

function rr = z(x) %Inverse of Yk
global A alpha beta;
if beta==1
rr = (x/A/alpha)^(1/(alpha-1));
else
rr = (alpha*(-(x/A/alpha)^beta*A+x)/(-1+alpha)/x)^(-beta/(beta-1));
end

function rr = S(c,k) %Saving Rate
rr = 1-c./feval(@Y,k);

% Matlab program for the deterministic continuous-time infinite-horizon Ramsey-Cass-Koopmans growth model.
% Patrick Toche©
% Please report any problems to patrick.toche@free.fr
% Program uses Matlab built-in boundary value problem routine bvp4c
% Inelastic labor supply.
% CES production function and CES utility function
% Definition of the parameters of the model
% scale production parameter                                  A
% capital income share                                        alpha
% elasticity of substitution capital-labour                   beta
% growth of labour-augmenting technological progress          g
% depreciation rate of capital                                d
% growth rate of population                                   n
% elaticity of intertemporal substitution in consumption      s
% discount rate (rate of time preference)                     b