Solving differential equations with state variables and co-states

The problem is quite straightforward in terms of max present value of the the sum of discounted utilities. Social discount rate is \rho. Utility function is given by log (c).

There are two equations of motion:

capital accumulation or investment = production - consumption

resource stock accumulation = - extraction of resources (there is a negative sign)

so there are two associated co-states.

production function is standard Cobb-Douglas with capital and non-renewable resource.

There are two controls, consumption and flow of resource extraction (call c and r) and two states, capital and stock of the non-renewable resource ( call k and s). The initial and final values of the states k and s are given. And the initial values of the co-states lambdas need to be found optimally. The graphs of the states, the co-states and the controls would be helpful.

It gives a system of differential equations and we need to find the optimal values of the initial values of lambdas.

%% Title: Solving Hamiltonian Objective Function using system of Differential Equations
% 1. Parameter Setup
clear; clc;
%%
% Preferences
rho = 0.05;         % Discount rate
%sigma = 2;          % Inverse of elasticity of intertemporal substitution (for CRRA)
u = @(c) log (c); % Utility function
%% 
% Technology
%delta = 0.05;       % Capital depreciation rate
alpha = 0.90;       % Capital share in production
beta = 1-alpha;        % Resource share in production
A = 3.968;              % TFP
f = @(k, r) A * k.^alpha .* r.^beta; % Cobb-Douglas production function
% Initial conditions
%K0 = 100;           % Initial capital stock
%S0 = 1000;          % Initial resource stock
K0 = 4.913;           % Initial capital stock
S0 = 11.5;          % Initial resource stock
KT = 0;              % final capital stock
ST = 0;               % final resource stock
% Time setup
T = 100;            % Time horizon
% dt = 0.1;           % Time step
dt = 10;           % Time step
time_periods = 0:dt:T;
N = length(time_periods);

%%
% Upon solving the Hamiltonian maximization with constraints, we get the
% following foc's
% 2. Define the System of ODEs (based on First-Order Conditions)
function d_vars = ode_system(t, vars, params)
% Unpack variables and parameters
k = vars(1);
s = vars(2);
lambda_k = vars(3);
lambda_s = vars(4);
params.rho = rho;
params.alpha = alpha;
params.beta = beta;
params.A = A;

% Optimal controls (consumption and resource extraction)
c = (lambda_k).^(-1);
r = (lambda_s / (lambda_k * A * beta * k^alpha)).^(1/(beta-1));
% System of ODEs
dk_dt = A * k^alpha * r^beta - c;
ds_dt = -r;
d_lambda_k_dt = lambda_k * (rho - A * alpha * k^(alpha-1) * r^beta);
d_lambda_s_dt = lambda_s * rho;

d_vars = [dk_dt; ds_dt; d_lambda_k_dt; d_lambda_s_dt];
end

%%
% 3. Solve using a Shooting Method
% Define terminal conditions
% S(T) = 0 and K(T) = 0 (transversality condition)
final_cond = @(vars_T) [vars_T(1); vars_T(2)];

% Define the function to find the initial shadow prices
find_initial_lambda = @(lambda0) (final_cond(ode_solver(lambda0)))';

% ODE solver setup
ode_solver = @(lambda0) solve_ode(lambda0, K0, S0, T, dt, @ode_system, params);

% Use fsolve to find the initial shadow prices (lambda_k0, lambda_s0)
x0 = [1; 1]; % Initial guess for lambda_k0 and lambda_s0
initial_lambdas = fsolve(find_initial_lambda, x0);

Unrecognized function or variable 'ode_solver'.

Error in solution>@(lambda0)(final_cond(ode_solver(lambda0)))' (line 64)
find_initial_lambda = @(lambda0) (final_cond(ode_solver(lambda0)))';
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fsolve (line 167)
fuser = feval(funfcn{3},x,varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% Final run with optimal initial values
[~, vars] = ode_solver(initial_lambdas);

% Extract optimal paths
k_optimal = vars(:, 1);
s_optimal = vars(:, 2);
lambda_k_optimal = vars(:, 3);
lambda_s_optimal = vars(:, 4);
c_optimal = (lambda_k_optimal).^(-1/sigma);
r_optimal = (lambda_s_optimal ./ (lambda_k_optimal .* A .* k_optimal.^alpha)).^(1/(beta-1));

%%
% 4. Auxiliary function to solve ODE
function [t_out, vars_out] = solve_ode(lambda0, k0, s0, T, dt, ode_func, params)
y0 = [k0; s0; lambda0(1); lambda0(2)];
t_span = [0 T];
[t_out, vars_out] = ode45(@(t, y) ode_func(t, y, params), t_span, y0);
end

Answer 1

[Torsten]

This code works, but gives no solution for the parameters given.

%% Title: Solving Hamiltonian Objective Function using system of Differential Equations
% 1. Parameter Setup
% Preferences
rho = 0.05;         % Discount rate
%% 
% Technology
alpha = 0.90;       % Capital share in production
beta = 1-alpha;        % Resource share in production
A = 3.968;              % TFP
% Initial conditions
k0 = 4.913;           % Initial capital stock
s0 = 11.5;          % Initial resource stock
% Final conditions
kT = 0;              % final capital stock
sT = 0;               % final resource stock
% Time setup
T = 100;            % Time horizon

params.rho = rho;
params.alpha = alpha;
params.beta = beta;
params.A = A;
params.k0 = k0;
params.s0 = s0;
params.kT = kT;
params.sT = sT;
params.T = T;

% Use fsolve to find the initial shadow prices (lambda_k0, lambda_s0)
lambda0 = [1; 1]; % Initial guess for lambda_k0 and lambda_s0
initial_lambdas = fsolve(@(lambda)find_initial_lambda(lambda,params), lambda0)

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

initial_lambdas = 
   0.9845 - 0.0000i
   1.0150 - 0.0000i

% Final run with optimal initial values
[~, t, vars] = find_initial_lambda(initial_lambdas,params);
% Extract optimal paths
k_optimal = vars(:, 1);
s_optimal = vars(:, 2);
lambda_k_optimal = vars(:, 3);
lambda_s_optimal = vars(:, 4);
c_optimal = (lambda_k_optimal).^(-1);
r_optimal = (lambda_s_optimal ./ (lambda_k_optimal .* A .* k_optimal.^alpha)).^(1/(beta-1));
% Now plot what you like
figure(1)
plot(t,k_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

figure(2)
plot(t,s_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

●

figure(3)
plot(t,lambda_k_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

figure(4)
plot(t,lambda_s_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

figure(5)
plot(t,c_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

figure(6)
plot(t,r_optimal)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

function [res,t_out,vars_out] = find_initial_lambda(lambda,params)
  k0 = params.k0;
  s0 = params.s0;
  kT = params.kT;
  sT = params.sT;
  T = params.T;
  y0 = [k0; s0; lambda(1); lambda(2)];
  t_span = [0 T];
  [t_out, vars_out] = ode45(@(t, y) ode_func(t, y, params), t_span, y0);
  res = [vars_out(end,1)-kT ; vars_out(end,2)-sT];
end

function d_vars = ode_func(t,vars,params)
  k = vars(1);
  s = vars(2);
  lambda_k = vars(3);
  lambda_s = vars(4);
  rho = params.rho;
  alpha = params.alpha;
  beta = params.beta;
  A = params.A;
  % Optimal controls (consumption and resource extraction)
  c = (lambda_k).^(-1);
  r = (lambda_s / (lambda_k * A * beta * k^alpha)).^(1/(beta-1));
  % System of ODEs
  dk_dt = A * k^alpha * r^beta - c;
  ds_dt = -r;
  d_lambda_k_dt = lambda_k * (rho - A * alpha * k^(alpha-1) * r^beta);
  d_lambda_s_dt = lambda_s * rho;
  d_vars = [dk_dt; ds_dt; d_lambda_k_dt; d_lambda_s_dt];
end

[SUPRATIM]

Hello @Torsten,

Thank you for your answer. Can you explain why there was no solution and that fsolve stopped?

I think there is some problem with the terminal states, they should be actually sT=0 and lambda_kT*kT = 0. Maybe this would eliminate the complex number problem? The time paths plotted are showing the right shape. And I also do not understand the res = [vars_out(end,1)-kT ; vars_out(end,2)-sT] should equal zero or be as close to zero? For the modified terminal states should res = [vars_out(end, 1)*vars_out(end, 3) -0, vars_out(end, 2) - sT] ?

Few things:

1. The vars_out or vars is not stored as any vector/matrix when I run the code.
2. After defining the ode system with function d_vars = ode_system(t, vars, params) and then defining d_vars= [], can the function be solved as a function of initial lambdas (or guesses)?
3. Can you please explain the steps like this? You have defined the function earlier and then d_vars at the end? Why so?
4. The vars_out or vars is not stored as any vector/matrix when I run the code. Please have the vars_out as a matrix storing all the variables, c_optimal, r_optimal, k, s, lambda_k, lambda_s in each column. It is easier the output that way.
5. Can the problem be please solved in time steps of 10 years with T= 100. And I also need the present value of consumption c which is Cmax = sum0 to sum T U(c_t)/(1+ rho)^t where U(c_t) = log(c) or log(c_optimal)?

Thank you for your help.