% Matlab program for the continuous-time infinite-horizon Ramsey-Poisson growth model.
% Patrick TocheŠ
% December 2002, modified 26th/12/2002
% Please report any problems to patrick.toche@free.fr
% Program uses Matlab built-in boundary value problem routine bvp4c
% Inelastic labor supply.
% Cobb-Douglas production function and CES utility function
% Definition of the parameters of the model
% elasticity of intertemporal substitution s
% discount rate (rate of time preference) rho
% rate of interest r
% growth of labour income g
% poisson jump probability q
% labour income level in current state y
% labour income level in other state b
function rp1
clear all; delete thediary.out; diary thediary.out;
disp('The continuous-time infinite-horizon partial-equilibrium Ramsey-Poisson model with CES utility. Patrick TocheŠ.')
global s rho r g q y b T0 T N NMax RelTol AbTol;
global css kss c0 k0 c00;
s=1/3;
rho=0.05;
r=0.04;
g=0.02;
q=0.025;
y=1;
b=0;
T0=25;
T=input('Horizon Length: ');
N=10;
NMax=100;
RelTol=1e-6;
AbsTol=1e-3;
tic;
fprintf('Parameters sigma = %7.5f, rho = %7.5f, r = %7.5f, g = %7.5f, q = %7.5f\n',s,rho,r,g,q);
%%% Compute the initial steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=s*(rho-r)+r;
fprintf('m = %7.5f\n',m)
kss=y/((s*(rho-r)+r)*((g/s+q+rho-r)/q)^(s)-(r-g));
css=feval(@C,kss);
ess=[css; kss];
fprintf('css = %7.5f, kss =%7.5f\n',css,kss)
%%% Define new parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c0=css;
%q=1.1*q; fprintf('Parameters q = %7.4f\n',q)
%c00=0.18; k0=0.0007*kss; %for s=1 and T0=5
%c00=2*css; k0=2*kss;
c00=0.2696; k0=0.09*kss; %for s=1/3 and T0=25
%%% Compute the new steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=s*(rho-r)+r;
fprintf('m = %7.5f\n',m)
kss=y/((s*(rho-r)+r)*((g/s+q+rho-r)/q)^(s)-(r-g));
css=feval(@C,kss);
ess=[css; kss];
fprintf('css = %7.5f, kss =%7.5f\n',css,kss)
%%% Computing eigenvalues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
odess=feval(@ODE,1,ess); %fprintf('odess = %7.5f\n',odess);
jss=feval(@J,1,ess); %fprintf('jss = %7.5f\n',jss);
Jac=eig(jss); disp('The eigenvalues are:'); disp(Jac);
tr=trace(jss); dt=det(jss); fprintf('tr = %7.9f, dt = %7.9f\n',tr,dt);
if T