% RAMSEY.M
% Finds the policy function for the Ramsey model for k < steady state
% c = per capita consumption
% k = per capita capital
% ks = steady state capital

% Initialize parameters
n = .014;           % population growth rate
d = .10;            % depreciation rate of aggregate capital
delta = n + d;      % depreciation rate of k
p = .065;           % rate of time preference
theta = 2;          % (inverse of) intertemporal elas. of substitution
a = .5;             % capital share (production per capita is k^a)
rho = p - n;        % exponent in dynastic utility function

% Calculate steady state values
ks = (a/(rho+delta))^(1/(1-a));
cs = (ks^a - delta*ks);

% Find the policy function
k0 = ks/80; kf = ks;
% MATLAB needs initial conditions, so run backwards from steady-state
% define ssd = ks - k
ssd0 = ks - k0;
ssdf = ks - kf;
global theta rho delta a ks cs
[ssd,c] = ode23('cprime',ssdf,ssd0,cs,.00001);
k = ks - ssd;

% Calculate some interesting functions of k
y = k.^a;                     % Output
s = (y-c)./y;                 % saving rate
ys = ks.^a;                   % steady state output
ss = (ys - cs)/ys;            % steady state saving rate
kdotis0 = y - delta*k;        % those (c,k) for which kdot = 0
gk = k.^(a-1) - c./k - delta; % capital growth rate
plot(k,c,k,kdotis0)           % Plots policy function c(k) & kdot = 0 schedule


---------------------
% CPRIME.M
function cdot = cprime(ssd,c);

k = ks - ssd;
if k == ks
     cdot = -.5*(rho + sqrt(rho*rho+4*cs*(1-a)*a*(ks^(a-2))/theta));
else
     cdot = c.*(rho + delta - a*(k.^(a-1))) / (theta*((k.^a) - c - delta*k));
end     % if k == ks