# 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:** 28

<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: [September 11, 2022, 8:56pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/28 "2022-09-11T20:56:42Z")

</div>

I just did a test drive w/ @davidanthoff’s dynamic programming solver [Judyp](https://github.com/anthofflab/Judyp.jl).jl.  
WOW!  
Here are two simulations: w/ capital starting below (left) & above (right) the steady state  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/a/ca9a728f9d57ae16219ca278327a319e9da129b1.png) ![image](https://global.discourse-cdn.com/julialang/original/3X/c/9/c9994be0aa8bd6a503a7d9b17930d80e4f454056.png)

```julia
using Judyp, Plots
par=(δ=0.1, ρ=0.05, z= 0.05 + 0.1 +.02, β=1/(1+0.05));

I_SS = (par.z-par.ρ-par.δ)/(par.ρ+par.δ);
K_SS = I_SS/par.δ;
K0 = 0.5*K_SS; # start below SS
K0 = 1.5*K_SS; # start above SS

p1 = DynProgProblem(); 

set_transition_function!(p1) do k, k_new, x, up, p
    k_new[1] = (1-p.δ)*k[1] + x[1] 
end

set_payoff_function!(p1) do s, x, p
    p.z*s[1] - x[1] - 0.5*(x[1])^(2.0)
end

set_discountfactor!(p1, par.β)

set_exogenous_parameters!(p1, par)

add_state_variable!(p1,Symbol("K_1"),K0,0.0,3.0,10) #10 nodes [0,3]

add_choice_variable!(p1,Symbol("I_1"),0., 1.) # bounds [0,1] initial 0.5

T_sim = 150; N_sim=0; # number Monte Carlo runs

res1 = solve(p1) # res = solve(problem, solvers=[solver], print_level=1)
simres1 = simulate(res1, T_sim, N_sim);

plot(legend=:topright, ylims = (0,2))
plot!(simres1.state_vars[:], c=1, lab="k")
plot!(1:T_sim, tt -> K_SS, c=1, l=:dash, lw=3, lab = "K_SS")
plot!(simres1.choice_vars[:], c=2, lab="i")
plot!(1:T_sim, tt -> I_SS, c=2, l=:dash, lw=3, lab = "K_SS")

```

It would be awesome to benchmark w/ [VFIToolkit](https://github.com/vfitoolkit/VFIToolkit-matlab).m and @zsunberg’s [POMDPs](https://github.com/JuliaPOMDP/POMDPs.jl).jl and @odow’s [SDDP](https://github.com/odow/SDDP.jl).jl.

---

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