# 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:** 1
**Showing post:** 5

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

---

_[View the full topic](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285)._