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.

[Torsten]

Could you include a more detailed mathematical problem formulation similar to what is done in the section "Finite Time" of the Wikipedia article about optimal control

https://en.wikipedia.org/wiki/Optimal_control

?

Can you explain why there was no solution and that fsolve stopped?

Either the boundary value problem is numerically hard to solve or your problem statement is incorrect or both.

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 complex numbers come from expressions like k^alpha or k^(alpha-1) if k becomes negative or r^beta if r becomes negative or (lambda_s / (lambda_k * A * beta * k^alpha)).^(1/(beta-1)) if lambda_s / (lambda_k * A * beta * k^alpha becomes negative.

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?

"fsolve" solves problems of the type res(x) = 0 and you have to supply res. kT and sT are set to zero in the params vector - thus you try to find lambda(1) and lambda(2) such that k(T) - kT = 0 and s(T) - sT = 0. And vars_out(end,1) = k(T) and vars_out(end,2) = s(T).

For the modified terminal states should res = [vars_out(end, 1)*vars_out(end, 3) -0, vars_out(end, 2) - sT] ?

Yes. But a*b = 0 means a = 0 or b = 0. Thus you better should test whether res = [vars_out(end, 3) -0, vars_out(end, 2) - sT] works.

The vars_out or vars is not stored as any vector/matrix when I run the code.

As far as I know, only variables computed in the script part are accessible in the workspace. But why do you need them ? Either the iteration converges or not - I don't think that intermediate results can help to get a better understand about what goes wrong. The final results are stored and plotted after the call to "find_initial_lambda" in the code.

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)?

I don't know what you mean. For given values of lambda(1) and lambda(2), the trajectories for k, s, lambda_k and lambda_s over time are obtained by calling

[~, t, vars] = find_initial_lambda(lambda,params);

as shown in the code.

Can you please explain the steps like this? You have defined the function earlier and then d_vars at the end? Why so?

I don't know what you mean. In the script part, you call "fsolve" to adjust the lambda values such that k(T) - kT = 0 and s(T) - sT = 0. In order to compute k(T) and s(T) within the iteration process for "fsolve", you have to call the ode integrator "ode45" within the function "find_initial_lambda", and the ode function is "ode_func" which returns the time derivatives in the array "d_vars".

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)?

If you mean you want the trajectories at t = 0, t=10,..., t=100, change t_span = [0 T]; in function "find_initial_lambda" to tspan = 0:10:100. Further, isn't it necessary for the utility function to appear somewhere in the system of differential equations to be solved ? I don't see a log(c) expression somewhere in the code.