# Solving the 4 quadrants of dynamic optimization problems in Julia. Help Wanted!

Solving dynamic optimization problems is at the heart of economics (and many other fields).
The goal of this post is to explore various tools for solving dynamic optimization problems in the Julia Eco-system and to hopefully see if there are opportunities for synergy between the many different communities that solve these types of problems.

This post will consider a firm that chooses capital investment to maximize its value.
Here are the 4 quadrants of modeling:

We will study the case F(K)=zK, p_t =1, and C(I,K)=I + \frac{1}{2} I^2 because it has simple closed form solutions (for comparison).
Here is the corresponding closed form solution:

More detail for the deterministic cases.

Once a user (economist/engineer/etc) chooses a modeling framework (e.g. deterministic/discrete-time) there are often many different ways to solve the same problem.

Hence, a modeler has 3 steps

1. Decide how to formulate the problem (e.g. Deterministic/Continuous-Time)
2. Decide which mathematical problem solve:
-Sometimes it’s easier to solve a system of three BV-DAEs
-Sometimes it’s easier to solve a system of two first-order BV-ODEs
-Sometimes it’s easier to solve one second-order BV-ODE
-Sometimes it’s easier to solve the HJB-pde (in this case it’s an ODE)
-Note: not all approaches involve solving equations (e.g. see InfiniteOpt.jl)
3. Decide which method(s) to use to solve the mathematical problem specified above
-E.g. there are many methods a modeler can use to solve a system of two first-order difference equations (perturbation, projection, and hybrids).

There are ENORMOUS trade-offs users need to make.
For example, perturbation methods (SolveDSGE.jl) are a powerful approach to solve large systems quickly, but generally work for models w/ well defined steady-states.
They also need to be extended to solve models w/ occasionally binding constraints (OccBin is not yet implemented in SolveDSGE.jl, see Julia examples here).

Dynamic programming approaches (such as VFI) can solve a very wide variety of problems, but they tend to be very-very slow. Currently, Julia does not have a generic equivalent to Matlab’s VFIToolkit. There was a Julia GSoC project that solved the optimal growth model w/ parallel VFI, but it was not generic.

Package/Sample Code Description Note
SolveDSGE.jl -Code link Solves difference equations w/ Perturbation/Projection Only for discrete time models w/ well defined steady-states @RJDennis
POMDPs.jl -Code link Solves Bellman equations w/ various methods discussion @zsunberg
JuMP.jl SDDP.jl Description Not sure it solves infinite horizon models. @ odow & co
Package Description Note

Package/Sample Code Description Note
SolveDSGE.jl -Code link Solves stochastic difference equations w/ Perturbation/Projection Only for discrete time models w/ well defined steady-states @RJDennis
POMDPs.jl Solves Bellman equations w/ various methods Please post code below using e.g. a 2-state Markov Chain for z. @zsunberg
Package Description Note

Package/Sample Code Description Note
-Code link Solves HJB equations using Ben Moll’s method My custom solver. Works pretty nicely.
EconPDEs.jl -Code link Solves HJB equations @matthieu
InfiniteOpt.jl Solves Hamiltonian Only finite horizon currently. @pulsipher
SciML.jl Solve DEs TO DO. Only finite horizon currently. Try @MainBatDealer & @ptoche trick to solve a DE w/ an infinite horizon boundary constraint. @ChrisRackauckas
Flux.jl -Code link Solve HJB w/ NN @Goosee
Package Description Note

Package/Sample Code Description Note
EconPDEs.jl Solves HJB equations To Do. @matthieu
Flux.jl Solves HJB equations w/ Deep Learning @Goosee
Package Description Note

There are many other dynamic programming methods (VFI/PFI/ECM…) that currently don’t have generic/robust implementations in Julia.
@sglyon has 7+ methods
@floswald has others
Unfortunately @sglyon’s code is not generic like VFIToolkit.m

To Do:
check out QuantEcon.jl
solve w/ SDDP.jl|
solve 2-state Markov Chain version w/ POMDPs.jl
solve CTS-time w/ SciML.jl
solve stochastic w/ EconPDEs.jl
solve infinite horizon w/ InfiniteOpt.jl
check out Dynare.jl (not maintained?)
Dolo.jl
DifferentiableStateSpaceModels.jl
MFGNet.jl

Oren Levintal has a state-of-the-art global package in Matlab
Robert Kirkby has a phenomenal VFI package in Matlab

Please let me know what I missed!

22 Likes

Here is a discussion on how to solve the finite-horizon version of this problem w/ InfiniteOpt .jl

2 Likes

The infinite-horizon continuous-time problems with a saddle-point structure could be solved using the time reversal technique outlined here. At some point some more general technique should be put together though.

Edit Adding some code for the simple deterministic saddle-path. As it turns out, the optimal trajectory features a constant level of investment (if I copied the equations above correctly); for this reason I have drawn the optimal trajectory with dashes, so we can see the black nullcline below it. Solving this problem is just like solving the Ramsey-Cass-Koopmans model (see the link above). Both models should be part of a general setup, as discussed below by Chris and Jesse.

using DifferentialEquations
using Plots

# Parameters
p = Dict(:δ => 0.1, :ρ => 0.05, :z => 0.3)

"""firm investment problem
K: state variable, with given initial value K0
I: control variable, free to jump"""
function inv!(du,u,p::Dict,t)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
K, I = u
du[1] = I - δ * K
du[2] = (ρ + δ - z) + (ρ + δ)*I
nothing
end

"""compute the stationary state"""
function ss(p::Dict)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
Iss = z / (ρ+δ) - 1
Kss = Iss/δ
return Dict(:Kss => Kss, :Iss => Iss)
end

# The stationary state
SS = ss(p)  # returns 10.0 and 1.0 for reference parameter values
Kss, Iss = ss(p)[:Kss], ss(p)[:Iss]

# event handler to prevent variables going out of bounds
function set_callback()
c1(u, t, integrator) = u[1]  # stay in positive range
c2(u, t, integrator) = u[2]  # stay in positive range
c3(u, t, integrator) = 2*Kss - u[1]  # don't stray too far from stationary state
c4(u, t, integrator) = 2*Iss - u[2]  # don't stray too far from stationary state
function terminate_affect!(integrator)
terminate!(integrator)
end
# set up affect upon callback
cb1 = ContinuousCallback(c1, terminate_affect!)
cb2 = ContinuousCallback(c2, terminate_affect!)
cb3 = ContinuousCallback(c3, terminate_affect!)
cb4 = ContinuousCallback(c4, terminate_affect!)
cb = CallbackSet(cb1, cb2, cb3, cb4)
return cb
end

cb = set_callback()

# Run time "backwards" | start below the stationary state
tspan = (0.0, -100.0)
tweak = 1.01  # I0 starts above Iss
u0 = [0.999999 * Kss; 0.999999 * tweak * Iss]
ode = ODEProblem(inv!, u0, tspan, p)
sol1 = solve(ode, abstol=1e-8, reltol=1e-8, callback=cb)

# Run time "backwards" | start above the stationary state
tspan = (0.0, -100.0)
tweak = 0.99 # I0 starts below Iss
u0 = [1.00001 * Kss; 1.00001 * tweak * Iss]
ode = ODEProblem(inv!, u0, tspan, p)
sol2 = solve(ode, abstol=1e-8, reltol=1e-8, callback=cb)

# initialize plot
plt = plot(title="Policy Function in the Neoclassical Investment Model", xlims=(0,2*Kss), ylims=(0,2*Iss), xlabel="K(t)", ylabel="I(t)", legend=false)

# plot nullclines
plot!(plt, x-> p[:δ] * x, 0, 2*Kss, color=:black, linewidth=3, legend=false)
plot!(plt, x-> Iss, 0, 2*Kss, color=:black, linewidth=3, legend=false)
plot!(plt, sol1, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Iss), color=:blue, linewidth=2, linestyle=:dash)
plot!(plt, sol2, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Iss), color=:blue, linewidth=2, linestyle=:dash)

scatter!((Kss,Iss), label=["Stationary State" "."], color=:red)

# make the streamplot
using CairoMakie

"""firm investment model equations as points in (K,I) space"""
function invEq(K, I, p)
δ, ρ, z = p[:δ], p[:ρ], p[:z]
dK = I - δ * K
dI = (ρ + δ - z) + (ρ + δ)*I
return Point(dK, dI)
end

let
fig = Figure(fontsize = 20)
ax = Axis(fig, xlabel="K", ylabel="I", title="Streamplot of Neoclassical Investment Model", backgroundcolor=:black)
streamplot!(ax, (x,y)->invEq(x,y,p), 0..2*Kss, 0..2*Iss, colormap = Reverse(:plasma),
gridsize= (32,32), arrow_size = 10)
fig[1, 1] = ax
display(fig)
save("investment-streamplot-black.png", fig)
end


1 Like

I find this interesting but the problem description is too economical (half-pun) for me to follow. Any chance that we can have more details on the mathematical problem, in particular the stochastic continuous quadrant?

1 Like

@mschauer
1: A firm begins period “t” w/ capital K_t
2: it uses the capital to produce output according to technology F(K_t)=z\times K_t, where z is a productivity parameter (constant in deterministic models, stochastic in stochastic models)
3: \delta\times K_t of a firms capital depreciates after production (\delta is a constant depreciation parameter)
4: the firm invests (buys capital) I_t, at purchase price I_t and pays quadratic installation cost .5*I_t^2
5: the firm pays the owners a dividend D_{t} \equiv zK_{t} - I_t - \frac{1}{2} I_t^2
-The firm faces a trade-off when investing: lower dividend today vs higher dividend in the future
6: capital evolves according to the Law of Motion (LOM)
K_{t+1}=(1-\delta)K_{t} +I_{t} (discrete-time)
\dot{K}=I_t -\delta K_{t} (continuous time)
-the firm begins the next period w/ K_{t+1} or K_{t+dt}

I fully specified the deterministic (discrete time & continuous time) models & solved them in closed form above (& compare them w/ the various numerical solutions).

The only missing component for the stochastic models is the productivity process.
I omitted it bc I haven’t (yet) solved the stochastic model in closed form.
It’s a bit of an art & science to specify a stochastic process that is interesting & tractable.
Here is code where I solve the stochastic-discrete time model w/ an AR-1 tech shock:
z_{t+1} = 0.5 + 0.5*z_{t} + ε_{t+1}

For the stochastic, continuous-time case, economists typically use a diffusion process, a Jump/Poisson process, or a Jump Diffusion.
I left that part to readers like you to specify your favorite process (maybe GBM?) that makes the problem not too hard to solve.

Warning: sometimes the process needs to be bounded for the solution to be well-defined
@matthieu sometimes uses a diffusion w/ reflecting boundaries that stays between [\underline{z},\bar{z}]

Researchers from economics/bio/physics/etc often use Distributions.jl when working w/ random variables.
The world would be a much better place if researchers from different fields came together to solve intertemporal optimization problems.

3 Likes

Note that this fills in a separate quadrant as well:

That’s stochastic discrete time, separating out the mathematical part from the DSGE modeling.

3 Likes

Got it!

1 Like

What @ChrisRackauckas is saying is important for all of these when we want to make composable packages (especially for AD). This cuts across some of the quadrants you have outlined.

What many people think of as a “solution” is really just spitting out a new dynamic process (potentially stochastic) as an IVP. Given those solutions, it would be great if we all used the same software rather than having those orthogonal parts built directly into the control-theoretic solvers themselves.

DifferenceEquations (when it is ready) and all of the existing SDE/DE/etc. stuff already in SciML for continuous time fulfill those parts of the puzzle.

5 Likes

Exactly. Essentially, there’s a few methods used here, along with their adjoints etc. to speedup over what Zygote would generate, so we want to implement it once in a package just about the mathematical problem and the solution process. We believe modeling packages should be kept separate from this, so say a DSGE modeling package should be separate from the solver, and reuse these solver components.

DiffEq has been used for a lot of this, where we have the jump process and the jump diffusion ones (along with the differential equations of course), but there’s been a noticeable gap in the interface for discrete time stochastic models.

3 Likes

Any discrete time stochastic dynamical system and any static hierarchical PPL model can be described as composition of Markov kernels, I think that is the next ecosystem-integrating abstraction we should pin down.

1 Like

In the same spirit as this post, @sdwfrost, solves the SIR model in all 4 quadrants here:

A few observations:

1. @sdwfrost solves this model 26+ ways, doing a great job of illustrating the Julia ecosystem. Bravo!
2. These models are generally models are generally initial value problems & there (currently) exist more packages & other resources for solving these types of problems.
Problems that economists are interested in tend to be boundary value problems (often w/ infinite horizon terminal conditions) which are less easy to solve w/ some “standard solvers” not designed for econ.
(Note: I’ve had some success w/ Mathematica, POMDPs.jl. I don’t see why discrete time Dynamic Programming solvers can’t be used for econ problems…)
3. The standard SIR model has no closed form solution.
My motto is “trust by verify”, hence I strongly prefer packages verify their solvers work w/ as many closed form solutions as possible. (This is why I focus on econ models w/ known closed form solutions here.)
A slightly modified SIR does have a closed form solution.

Thanks to @ChrisRackauckas for the link.

1 Like