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

**URL:** <https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285>\
**Category:** Finance and Economics\
**Tags:** optimization\
**Created:** [December 17, 2021, 10:38pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285 "2021-12-17T22:38:40Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [December 17, 2021, 10:38pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/1 "2021-12-17T22:38:40Z")

</div>

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:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/b/dbc4f12246916003d6ba02876a468aa8df59e718.png)

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:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/5/55b82ebfe900dd900b6ffb60c0c78e0c5a0ffd88.png)

More detail for the deterministic cases.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/9/b9253e91502d1fb0988fc2849a4e45fd4d3fecd8.jpeg)

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.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/3/933b30d8a53cbf24b06d592aa1eeee786ab3bf3a.jpeg)

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](https://github.com/pulsipher/InfiniteOpt.jl).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](https://github.com/RJDennis/SolveDSGE.jl).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](https://github.com/RJDennis/SolveDSGE.jl).jl, see [Julia examples here](https://github.com/pcuba/RED-PLC-Estimation)).

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](https://github.com/vfitoolkit/VFIToolkit-matlab). There was a Julia GSoC project that solved the optimal growth model w/ parallel VFI, but it was not generic.

**Quadrant: Discrete-time/deterministic**

| Package/Sample Code | Description | Note |
| --- | --- | --- |
| [SolveDSGE](https://github.com/RJDennis/SolveDSGE.jl).jl -[Code link](https://github.com/azev77/4Quadrants/blob/main/DiscreteTime_Deterministic_SolveDSGE.jl) | Solves difference equations w/ Perturbation/Projection | Only for discrete time models w/ well defined steady-states @RJDennis |
| [POMDPs](https://github.com/JuliaPOMDP/POMDPs.jl).jl -[Code link](https://github.com/azev77/4Quadrants/blob/main/DiscreteTime_Deterministic_POMDPs.jl) | Solves Bellman equations w/ various methods | [discussion](https://github.com/JuliaPOMDP/POMDPs.jl/discussions/351) @zsunberg |
| [JuMP](https://github.com/jump-dev/JuMP.jl).jl [SDDP](https://github.com/odow/SDDP.jl).jl | Description | Not sure it solves infinite horizon models. @ odow & co |
| Package | Description | Note |

**Quadrant: Discrete-time/Stochastic**

| Package/Sample Code | Description | Note |
| --- | --- | --- |
| [SolveDSGE](https://github.com/RJDennis/SolveDSGE.jl).jl -[Code link](https://github.com/azev77/4Quadrants/blob/main/DiscreteTime_Stochastic_SolveDSGE.jl) | Solves stochastic difference equations w/ Perturbation/Projection | Only for discrete time models w/ well defined steady-states @RJDennis |
| [POMDPs](https://github.com/JuliaPOMDP/POMDPs.jl).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 |

**Quadrant: Continuous -time/deterministic**

| Package/Sample Code | Description | Note |
| --- | --- | --- |
| -[Code link](https://github.com/azev77/4Quadrants/blob/main/ContinuousTime_Deterministic_Zevelev.jl) | Solves HJB equations using Ben Moll’s method | My custom solver. Works pretty nicely. |
| [EconPDEs](https://github.com/matthieugomez/EconPDEs.jl).jl -[Code link](https://github.com/azev77/4Quadrants/blob/main/ContinuousTime_Deterministic_EconPDEs.jl) | Solves HJB equations | @matthieu |
| [InfiniteOpt](https://github.com/pulsipher/InfiniteOpt.jl).jl | Solves Hamiltonian | Only finite horizon currently. @pulsipher |
| [SciML](https://github.com/SciML).jl | Solve DEs | TO DO. Only finite horizon currently. Try @MainBatDealer & @ptoche [trick](https://discourse.julialang.org/t/solving-boundary-value-differential-equation-problems-for-economics/72871) to solve a DE w/ an infinite horizon boundary constraint. @ChrisRackauckas |
| [Flux](https://discourse.julialang.org/t/solving-hjb-pde-using-deep-learning/60639/12).jl -[Code link](https://github.com/azev77/4Quadrants/blob/main/ContinuousTime_Deterministic_NeuralNetwork.jl) | Solve HJB w/ NN | @Goosee |
| Package | Description | Note |

**Quadrant: Continuous -time/Stochastic**

| Package/Sample Code | Description | Note |
| --- | --- | --- |
| [EconPDEs](https://github.com/matthieugomez/EconPDEs.jl).jl | Solves HJB equations | To Do. @matthieu |
| [Flux](https://discourse.julialang.org/t/solving-hjb-pde-using-deep-learning/60639/12).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](https://github.com/sglyon/CLMMJuliaPythonMatlab)  
@floswald has others  
Unfortunately @sglyon’s code is not generic like [VFIToolkit](https://github.com/vfitoolkit/VFIToolkit-matlab).m

To Do:  
check out [QuantEcon](https://github.com/QuantEcon/QuantEcon.jl).jl  
solve w/ [SDDP](https://github.com/odow/SDDP.jl).jl|  
solve 2-state Markov Chain version w/ [POMDPs](https://github.com/JuliaPOMDP/POMDPs.jl).jl  
solve CTS-time w/ [SciML](https://github.com/SciML).jl  
solve stochastic w/ [EconPDEs](https://github.com/matthieugomez/EconPDEs.jl).jl  
solve infinite horizon w/ [InfiniteOpt](https://github.com/pulsipher/InfiniteOpt.jl).jl  
check out [Dynare](https://git.dynare.org/camilomrch/Dynare.jl).jl (not maintained?)  
[Dolo](https://github.com/EconForge/Dolo.jl).jl  
[DifferentiableStateSpaceModels](https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl).jl  
[MFGNet](https://github.com/EmoryMLIP/MFGnet.jl).jl

[Oren Levintal](https://sites.google.com/site/orenlevintal/taylor-projection) has a state-of-the-art global package in Matlab  
[Robert Kirkby](https://github.com/vfitoolkit/VFIToolkit-matlab) has a phenomenal VFI package in Matlab

Please let me know what I missed!

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [December 30, 2021, 11:29pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/4 "2021-12-30T23:29:43Z")

</div>

Here is a discussion on how to solve the finite-horizon version of this problem w/ [InfiniteOpt](https://github.com/pulsipher/InfiniteOpt.jl).jl  
[https://github.com/pulsipher/InfiniteOpt.jl/discussions/194](https://github.com/pulsipher/InfiniteOpt.jl/discussions/194)

---

<div class="post-metadata">

**Author:** ![ptoche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ptoche/32/23554_2.png) [@ptoche](https://discourse.julialang.org/u/ptoche)\
**Post date:** [January 11, 2022, 4:05am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/5 "2022-01-11T04:05:15Z")

</div>

The infinite-horizon continuous-time problems with a saddle-point structure could be solved using the time reversal technique outlined [here](https://discourse.julialang.org/t/solving-boundary-value-differential-equation-problems-for-economics/72871/17). 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.

```julia
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)

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

savefig(plt, "investment-saddle.png")

# 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

```

 ![investment-streamplot-black](https://global.discourse-cdn.com/julialang/original/3X/e/4/e4a39bae0280ed6be145eebb2f61f29ce5a664c8.png)

![investment-saddle](https://global.discourse-cdn.com/julialang/original/3X/e/2/e2fe65e7425b20200fae2c2b885cdbb734f9f9ad.png)

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [January 11, 2022, 10:05am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/6 "2022-01-11T10:05:53Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [January 11, 2022, 4:45pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/7 "2022-01-11T16:45:23Z")

</div>

@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](https://github.com/azev77/4Quadrants/blob/main/DiscreteTime_Stochastic_SolveDSGE.jl) 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](https://github.com/JuliaStats/Distributions.jl).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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 11, 2022, 4:54pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/8 "2022-01-11T16:54:41Z")

</div>

Note that this fills in a separate quadrant as well:

[https://github.com/SciML/DifferenceEquations.jl](https://github.com/SciML/DifferenceEquations.jl)

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

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [January 11, 2022, 5:28pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/9 "2022-01-11T17:28:16Z")

</div>

Got it!

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [January 11, 2022, 8:50pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/11 "2022-01-11T20:50:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 11, 2022, 9:49pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/12 "2022-01-11T21:49:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [January 12, 2022, 9:24am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/13 "2022-01-12T09:24:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 5, 2022, 9:54pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/14 "2022-02-05T21:54:11Z")

</div>

In the same spirit as this post, @sdwfrost, solves the SIR model in [all 4 quadrants here](https://github.com/epirecipes/sir-julia):  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/3/d376d755ee7a5757515ff9e73ea50ac8bbfa2783.png)  
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](https://arxiv.org/abs/2007.16069).

Thanks to @ChrisRackauckas for the link.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [April 27, 2022, 1:02am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/15 "2022-04-27T01:02:14Z")

</div>

Your recent post made me remember about this thread. This book was recently released and it seems to address all possible solution approaches to discrete dynamic problems:  
[https://www.wiley.com/en-au/Reinforcement+Learning+and+Stochastic+Optimization%3A+A+Unified+Framework+for+Sequential+Decisions-p-9781119815044](https://www.wiley.com/en-au/Reinforcement+Learning+and+Stochastic+Optimization%3A+A+Unified+Framework+for+Sequential+Decisions-p-9781119815044)

The description says several things that could be potentially useful for economics:

> _Reinforcement Learning and Stochastic Optimization_ offers a single canonical framework that can model any sequential decision problem using five core components: state variables, decision variables, exogenous information variables, transition function, and objective function. This book highlights twelve types of uncertainty that might enter any model and pulls together the diverse set of methods for making decisions, known as policies, into four fundamental classes that span every method suggested in the academic literature or used in practice.
> 
> _Reinforcement Learning and Stochastic Optimization_ is the first book to provide a balanced treatment of the different methods for modeling and solving sequential decision problems, following the style used by most books on machine learning, optimization, and simulation.

Could it be useful for economists?

---

<div class="post-metadata">

**Author:** ![youainti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/youainti/32/30451_2.png) [@youainti](https://discourse.julialang.org/u/youainti)\
**Post date:** [May 27, 2022, 6:37pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/16 "2022-05-27T18:37:06Z")

</div>

Currently the link to the book is broken: Here is an updated link

> **[Reinforcement Learning and Stochastic Optimization: A Unified Framework for...](https://www.wiley.com/en-us/Reinforcement+Learning+and+Stochastic+Optimization%3A+A+Unified+Framework+for+Sequential+Decisions-p-9781119815051)**
>
> REINFORCEMENT LEARNING AND STOCHASTIC OPTIMIZATION Clearing the jungle of stochastic optimization Sequential decision problems, which consist of “decision, information, decision, information,” are ubiquitous, spanning virtually every human...

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [June 27, 2022, 6:43pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/17 "2022-06-27T18:43:50Z")

</div>

I have written the package [SpeedMapping.jl](https://github.com/NicolasL-S/SpeedMapping.jl) to accelerate the convergence of fixed-point mappings. You may want to try it to speed up applications such as VFIs. Feel free to contact me for details.

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [June 27, 2022, 7:24pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/18 "2022-06-27T19:24:55Z")

</div>

@nicolas It would be greeat if this sort of thing could be built into the sciml wrapper ecosystem - if it isn’t already. Then it could be easily benchmarked reatlive to anderson acceleration/etc. @ChrisRackauckas might have other plans for fixed point agorithms/noninear solvers but I think that [NonlinearSolve.jl: High-Performance Unified Nonlinear Solvers · NonlinearSolve.jl](https://nonlinearsolve.sciml.ai/dev/) might be the plan?

Also, the wrappers implement (or will eventually) forward and reverse-mode AD rules, so you get differentiable programming for free.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [June 27, 2022, 11:11pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/19 "2022-06-27T23:11:09Z")

</div>

@nicolas I’m very curious to apply your package. Can you please help me.  
Let’s consider the discrete-time, deterministic example.  
Here is my boring/slow VFI code:

```julia
δ=0.10; r=0.15; z= r + δ + .02;
u(K,I;z=z)= z*K -I -(.5)*(I^(2.0))
#
K_ss = (z - (r+δ))/(δ*(r+δ))
D_ss = u(K_ss, δ*K_ss)
V_ss = ((1+r)/r)*D_ss
Q_ss = (1+r)*(1+δ*K_ss)
V_cf(k) = (V_ss - Q_ss*K_ss) + Q_ss*k
Grid_K = 0.01*K_ss : .01 : 3.00*K_ss
#
using Interpolations, Plots, LinearAlgebra
approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc = Line())
struct MDP γ; S; A; T; R end
RHS(P,V,s,a) = P.R(s,a) + P.γ*V(T(s,a))
#
γ = (1.0 / (1.0 + r))
S = [Grid_K;]
A(s;δ=δ) = [range( -(1.0 -δ)*s +.01, maximum(S), length=length(S));]  
T(s,a) = (1.0 -δ)*s + a
R(s,a) = u(s,a)
P1 = MDP(γ, S, A, T, R)
V0 = [u(s, δ*s) for s in P1.S]
#
function MDP_VFI(P::MDP, k_max::Int64, V)
    pol = similar(V)
    for k = 1:k_max 
        v̂ = approx(P.S, V)
        rhs = [[RHS(P,v̂,s,a) for a ∈ P.A(s)] for s ∈ P.S ]
        rhs = hcat(rhs...)
        Vp = maximum(rhs, dims=1) |> vec
        V = Vp 
        #plot!(Vp) |> display
        #pol = (k==k_max) ? argmax(rhs, dims=1) |> vec : nothing 
    end
    return V #, pol
end
V = MDP_VFI(P1, 500, V0)

plot(S, V)
plot!(S, V_cf.(S))

```

Start with a guess V\_0(k) over `Grid_K `.  
Update rule:  
V\_{t+1}(K) =\underset{I}{\max} \{ u(K, I) +\frac{1}{1+r}V\_{t}( I +(1-\delta)K ) \}  
Stop when: d(V\_{t+1}(K), V\_{t}(K))\<\varepsilon

How do I solve VFI problems w/ your package?!?!

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [June 28, 2022, 2:41am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/20 "2022-06-28T02:41:54Z")

</div>

@Albert_Zevelev I have modified your function to slightly reduce memory allocation (just to make the speed comparison as accurate as possible) and make the loop stop sooner (when norm(V - Vp) \< 1e-10).

To use the package, just write a map V → Vp and feed it to speedmapping with a starting point. On my computer, I get 3.3s running the accelerated loop and 1.7s using SpeedMapping, a rather modest improvement. On top of my head, I would say your problem converges relatively quickly already. Perhaps also the discretization (taking the maximum over array elements) may reduce the accuracy of the extrapolation compared to an ideal case with an infinite number of support points.

You also mentioned Anderson Acceleration. I have adapted to Julia an improved version of AA based on a paper by Henderson and Varadhan (2019) which I also tested out of curiosity. The speed is the same as SpeedMapping so either algorithm would probably do.

Here is the code:

```julia
δ = 0.10;
r = 0.15;
z = r + δ + 0.02;
u(K, I; z=z) = z * K - I - (0.5) * (I^(2.0))
#
K_ss = (z - (r + δ)) / (δ * (r + δ))
D_ss = u(K_ss, δ * K_ss)
V_ss = ((1 + r) / r) * D_ss
Q_ss = (1 + r) * (1 + δ * K_ss)
V_cf(k) = (V_ss - Q_ss * K_ss) + Q_ss * k
Grid_K = 0.01*K_ss:0.01:3.00*K_ss
#
using Interpolations, Plots, LinearAlgebra
approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc=Line())
struct MDP
    γ
    S
    A
    T
    R
end
RHS(P, V, s, a) = P.R(s, a) + P.γ * V(T(s, a))
#
γ = (1.0 / (1.0 + r))
S = [Grid_K;]
A(s; δ=δ) = [range(-(1.0 - δ) * s + 0.01, maximum(S), length=length(S));]
T(s, a) = (1.0 - δ) * s + a
R(s, a) = u(s, a)
P1 = MDP(γ, S, A, T, R)
V0 = [u(s, δ * s) for s in P1.S]

function MDP_VFI(P::MDP, k_max::Int64, Vin; tol=1e-10)
    # pol = similar(V)
    V = copy(Vin)
    Vp = similar(V)
    for k = 1:k_max
        v̂ = approx(P.S, V)
        rhs = [[RHS(P, v̂, s, a) for a ∈ P.A(s)] for s ∈ P.S]
        rhs = hcat(rhs...)
        Vp .= maximum(rhs, dims=1) |> vec
        if norm(V - Vp) < tol
            break
        end
        V .= Vp
        #plot!(Vp) |> display
        #pol = (k==k_max) ? argmax(rhs, dims=1) |> vec : nothing 
    end
    return V #, pol
end
V_VFI = MDP_VFI(P1, 500, V0)

plot(S, V)
plot!(S, V_cf.(S))

# To use SpeedMapping, write the loop as a map
function map!(Vout, Vin, P)
    v̂ = approx(P.S, Vin)
    rhs = [[RHS(P, v̂, s, a) for a ∈ P.A(s)] for s ∈ P.S]
    rhs = hcat(rhs...)
    Vout .= maximum(rhs, dims=1) |> vec
end

using SpeedMapping
V_SM = speedmapping(V0; (m!)=(Vout, Vin) -> map!(Vout, Vin, P1), tol=1e-10).minimizer

# Solutions are similar
norm(V_VFI - V_SM)

using BenchmarkTools
@benchmark V = MDP_VFI(P1, 500, V0)
@benchmark speedmapping(V0; (m!)=(Vout, Vin) -> map!(Vout, Vin, P1), tol=1e-10)

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [June 28, 2022, 3:53am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/21 "2022-06-28T03:53:28Z")

</div>

Thanks! It works well! The bigger the grid, the more advantage to using your package.

1. Both solvers still could be much faster  
My guess is I could use a faster `approx()` among other things…
2. Would it be hard to generalize this solver to n-dimensions?  
(similar “generic” VFI solvers like [VFIToolkit.m](https://github.com/vfitoolkit/VFIToolkit-matlab) have been floating around, but nothing in Julia)  
I have some code implementing other algorithms including [7+ methods](https://github.com/sglyon/CLMMJuliaPythonMatlab) …

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [June 28, 2022, 4:31am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/22 "2022-06-28T04:31:35Z")

</div>

Sure. It does work with n dimensions already, as well as other number types like Float32 or BigFloat (on Github, I also have a new version for Complex numbers, but I’m waiting a bit to release it officially in case bugs are found). If you have some application that doesn’t work, just ask and I’ll see what I can do.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 12, 2022, 6:42pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/23 "2022-08-12T18:42:15Z")

</div>

Update:  
I was able to formulate this as a Linear Quadratic Regulator optimal control problem & solved it w/ @andreasvarga’s amazing package [MatrixEquations](https://github.com/andreasvarga/MatrixEquations.jl).jl.

I did the infinite horizon, discrete time, deterministic version.  
LQR/LQG can be used in all 4 quadrants for both finite/infinite horizon problems!

Formulation (note: there is more than 1 way to formulate this problem as an LQR. be careful!)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/2/a2f6032a1fa80eea7d42cce01386cf20bb5b33bc.jpeg)

Code:

```julia
using MatrixEquations, LinearAlgebra; 
δ=0.1; r=0.05; z= r + δ +.02; β=1.0/(1.0 + r); # parameters

I_SS = (z-r-δ)/(r+δ);
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0+r)/r)*D_SS;
Q_SS = (1+r)*(1+I_SS);

P_sol = [(V_SS - Q_SS*K_SS) 0.0; Q_SS 0.0]; # closed form P. Value = x'*P*x
F_sol = [-I_SS 0.0]; # closed form u. Policy = -F_sol*x

A=[1.0 0.0; 0.0 (1.0 - δ)];
B=[0.0; 1.0];
Q=[0.0 0.0; z 0.0];
R=[-0.5];
S=[-1.0; 0.0];

S = 0.5*S; # my Riccatti is different from MatrixEquations.jl
#Q = -1.0*Q; R = -1.0*R; S = -1.0*S; # max to min. DO NOT NEED!
A = sqrt(β)*A; B = sqrt(β)*B; # scale by discount factor
Q = (Q + Q') ./ 2.0; R = (R + R') ./ 2.0; P_sol = (P_sol + P_sol') ./ 2.0; # make symmetric  
                       
P, CLSEIG, F = ared(A,B,R,Q,S); # solve for numerical value function = x'*P*x

# Compare Value function P. CORRECT!
julia> P
2×2 Matrix{Float64}:
 0.186667 0.595
 0.595 0.0

julia> P_sol
2×2 Matrix{Float64}:
 0.186667 0.595
 0.595 0.0

# Compare Policy function F. u = -F*x. CORRECT!
julia> F
1×2 Matrix{Float64}:
 -0.133333 -0.0

julia> F_sol
1×2 Matrix{Float64}:
 -0.133333 0.0

Fn = pinv(R .+ (B')*P*B) * ((B')*P*A +S'); # numerical policy = -F_n*x. SAME as F from ared().
julia> Fn
1×2 Matrix{Float64}:
 -0.133333 0.0

```

LQR problems that are used regularly by engineers in robotics, aeronautics & other fields are also common in economics.  
The world would be a much better place if we coordinated on (one or few) high quality Julia implementations, where engineers/economists/nurses/plumbers all submitted bug reports, feature requests, pull requests…

[Next page](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285.md?page=2)