# 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:** 16\
**Page:** 2

<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 28, 2022, 11:18pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/24 "2022-08-28T23:18:21Z")

</div>

It turns out the firm investment model used in the example above can be solved very precisely in all 4 quadrants using e.g. `MatrixEquations.jl`.  
Hat tip to @baggepinnen & @andreasvarga for their help  
Here are infinite horizon versions of all 4 quadrants:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/a/3acb46d5f5d708c8b8579faec5849eb2424d635b.png)  
Here is the math to formulate the firm investment problem as LQ/LQG:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/4/94f85a293d912a93d720dcaf20370e40e7bb4013.png)  
I’ll post the code later…

1. Above I solve infinite horizon versions in all 4 quadrants.  
`MatrixEquations.jl` currently does not solve finite horizon versions (that I could find…)  
[QuantEcon](https://github.com/QuantEcon/QuantEcon.jl/blob/master/src/lqcontrol.jl).jl solves the finite horizon version, only in discrete time (not continuous time)  
Would be nice if the Riccatti solvers in both packages were benchmarked/compared…
2. Symbolics.jl works insanely well (at least for what I needed so far).  
Though it took me a lot of time to figure out how to use it. I’m sure the docs will improve…

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [September 10, 2022, 10:49am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/25 "2022-09-10T10:49:21Z")

</div>

Sorry for interferring in this discusssion, but I have the feeling that the function `ared` available in [MatrixEquations](https://github.com/andreasvarga/MatrixEquations.jl) is improperly used in this example. My background is control engineering and the points discussed below reflect my understanding of the formulation of LQG problems in control theory.

1. For a given pair (A,B) (I will discuss its properties at point 2), `ared` essentially computes an “optimal” stabilizing state feedback F such that the eigenvalues of A-B\*F are stable (i.e., have moduli less than one for a discrete-time setting). For this, it is necessary that the pair (A,B) is _stabilizable_ (i.e., such a feedback must exist). Moreover, it is standard when solving LQG problems to choose the control and state weighting matrices R and Q _positive semidefinite_. For this setting, `arec` computes an “optimal” stabilizing solution as expected.

2. In this example the pair (A,B) is not stabilizable, because the unstable eigenvalue in 1 remains fixed for any feedback. Moreover, R is negative definite and Q is indefinite (just a symmetric matrix with both positive and negative eigenvalues). With this setting, `ared` will **guarantedly** fail.

3. By replacing (A,B) with the stabilizable pair `(sqrt(β)*A,sqrt(β)*B)`, and with the above choice of R and Q, a “solution” F is computed by `ared`, which has however no practical effect on the resulting closed-loop dynamics, because the closed-loop eigenvalues computed in CLSEIG are the same as the eigenvalues of `sqrt(β)*A`. Therefore, the usefulness of this result seems to me highly questionable.

4. If R is replaced by -R and Q = I is used, then at least one of the closed loop eigenvalue will be improved.

---

<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 10, 2022, 9:16pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/26 "2022-09-10T21:16:51Z")

</div>

1. Thank you very much for taking the time to write your feedback.  
The whole point of this post is to bring together people from different fields who solve intertemporal optimization problems.  
Nobel laureate [Tom Sargent](https://www.youtube.com/watch?v=_7zvFkM0NSQ) brought many tools from engineering (LQR, Kalman filter) to [1st year graduate economics](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjBrt3Q7Yr6AhXHGVkFHe01DdMQFnoECDwQAQ&url=http%3A%2F%2Fwww.sfu.ca%2F~kkasa%2FRecursive_Macroeconomic_Theory_Ljungqvist_Sargent_2018.pdf&usg=AOvVaw08_9a0sGfBJvUPeNi6lVkn).  
I wish we all collaborated on solving control problems.  
For those who don’t know, @andreasvarga is an author of the Systems and Control Library [SLICOT](http://slicot.org/) implemented in Fortran. He has shown that comparable performance can be achieved in Julia.
2. The deterministic discrete time problem I seek to solve (the simplest I could think of):  
V(x\_{t}) = \max\_{u} \sum\_{i=0}^{i=\infty} \beta^{i} \left( z \times x\_{t+i} - u\_{t+i} -\frac{1}{2} u\_{t+i}^{2} \right)   
subject to: x\_{t+1} = (1-\delta) x\_{t} + u\_{t}  
Initial condition: x\_{0} given  
Terminal condition: \lim\_{t \to \infty} \beta^{t} \lambda\_{t} x\_{t+1} =0  
Where:  
State variable: x\_{t} {my current stock of capital}  
Control variable: u\_{t} {capital investment}  
Parameters: \delta \in (0,1), \rho \>0 \Rightarrow \beta \equiv\frac{1}{1+\rho} \in (0,1), z\>\rho +\delta  
Reward function: z \times x\_{t} - u\_{t} -\frac{1}{2} u\_{t}^{2} {profit = revenue - cost of investment}  
Note: investment has two costs, purchase cost u\_{t} and quadratic installation cost \frac{1}{2} u\_{t}^{2}  
Note: maximizing f is equivalent to minimizing -f, we can rewrite the sequential problem:  
V(x\_{t}) = \min\_{u} \sum\_{i=0}^{i=\infty} \beta^{i} (-1)\times \left( z \times x\_{t+i} - u\_{t+i} -\frac{1}{2} u\_{t+i}^{2} \right)   
Unlike the previous note: I will proceed by formulating this as a minimization problem  
Bellman equation:  
V(x\_{t}) = \min\_{u} \left\{ -1 \times \left( z \times x\_{t} - u\_{t} -\frac{1}{2} u\_{t}^{2}\right) + \beta V(x\_{t+1}) \right\}   
x\_{t+1} = (1-\delta) x\_{t} + u\_{t}  
This simple problem has a closed form solution:  
Policy function: u(x\_{t}) = \frac{z}{\rho + \delta} -1  
Value function: V(x\_{t})= -\left( V\_{ss} - Q\_{ss} K\_{ss} \right) - z \frac{1 + \rho }{\delta + \rho} x\_{t} {find V\_{ss}, Q\_{ss}, K\_{ss} in posts above}
3. Since this problem has a quadratic reward/cost function and a linear transition function, it can be formulated as a discrete time LQR.  
Define the augmented state vector \hat{x}\_{t}' \equiv \begin{bmatrix} 1 & x\_{t}\end{bmatrix}  
f(x\_{t}, u\_{t}) = \hat{x}\_{t}' Q \hat{x}\_{t} + u\_{t}' R u\_{t} + 2 \hat{x}\_{t}' S u\_{t} = -1 \times \left( z \times x\_{t} - u\_{t} -\frac{1}{2} u\_{t}^{2}\right)  
where:  
Q= -1 \times \begin{bmatrix} 0 & z/2\\ z/2 & 0 \end{bmatrix} \Rightarrow \hat{x}\_{t}' Q \hat{x}\_{t} = -1 \times\begin{bmatrix} 1 & x\_{t}\end{bmatrix} \begin{bmatrix} 0 & z/2\\ z/2 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ x\_{t}\end{bmatrix} = -z \times x\_{t}  
R= -1 \times \begin{bmatrix} -\frac{1}{2} \end{bmatrix} \Rightarrow u\_{t}' R u\_{t} =\frac{1}{2} u\_{t}^{2}  
S= -1 \times\begin{bmatrix} -1/2 \\ 0 \end{bmatrix} \Rightarrow 2 \hat{x}\_{t}' S u\_{t} = -1 \times2 \begin{bmatrix} 1 & x\_{t}\end{bmatrix} \begin{bmatrix} -1/2 \\ 0 \end{bmatrix} \begin{bmatrix} u\_{t}\end{bmatrix} = u\_{t}   
x\_{t+1} = g(x\_{t}, u\_{t}) = A \hat{x}\_{t} + B u\_{t} = (1-\delta) x\_{t} + u\_{t}  
where  
A = \begin{bmatrix} 1 & 0 \\ 0 & (1-\delta)\end{bmatrix}, B = \begin{bmatrix} 0 \\ 1\end{bmatrix} \Rightarrow \begin{bmatrix} 1 \\ x\_{t+1} \end{bmatrix} = A \hat{x}\_{t} + B u\_{t}= \begin{bmatrix} 1 \\ (1-\delta) x\_{t} + u\_{t} \end{bmatrix}  
Some LQR solvers (like QuantEcon.jl) take the discount factor \beta as input, MatrixEquations.jl does not.  
Following your advice, define A=\sqrt{\beta}A, B=\sqrt{\beta}B.

```julia
using MatrixEquations, LinearAlgebra; 

δ=0.1; ρ=0.05; z= ρ + δ +.02; β=1/(1+ρ); # parameters

A= √β*[1.0 0.0; 0.0 (1.0 - δ)];
B= √β*[0.0; 1.0];
Q=-1*[0.0 z/2; z/2 0.0];
R=-1*[-0.5;;];  
S=-1*[-1/2; 0.0];

P, CLSEIG, F = ared(A,B,R,Q,S);

# Issue 1: both eigenvalues of A-B*F are less than 1. 
# I don't understand the problem???
julia> CLSEIG
2-element Vector{ComplexF64}:
 0.8783100656536799 - 0.0im
 0.9759000729485324 - 0.0im

julia> eigvals(A-B*F)
2-element Vector{Float64}:
 0.8783100656536799
 0.9759000729485332

julia> eigvals(A)
2-element Vector{Float64}:
 0.8783100656536799
 0.9759000729485332

# Issue 2: R is positive definite & Q is indefinite.
# is that a problem? 
julia> eigvals(R)
1-element Vector{Float64}:
 0.5
#
julia> eigvals(Q)
2-element Vector{Float64}:
 -0.085
  0.085

# Issue 3: if ared() gives the wrong answer, 
# why is it exactly equal to the closed form solution?
I_SS = (z-ρ-δ)/(ρ+δ);
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0+ρ)/ρ)*D_SS;
Q_SS = (1+ρ)*(1+I_SS);

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

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

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

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

```

I would love to get to the bottom of this.

Update: while the matrix `Q` above is indeed indefinite, recall the augmented state vector is \hat{x}\_{t}' \equiv \begin{bmatrix} 1 & x\_{t}\end{bmatrix}, hence \hat{x}\_{t}' Q \hat{x}\_{t} = -z \times x\_{t} for any x\_{t}.  
Economists are only interested in the case where capital is strictly positive x\_{t}\>0, hence the quadratic form will always be strictly negative \hat{x}\_{t}' Q \hat{x}\_{t} = -z \times x\_{t}\<0.  
Not sure if this is relevant.

Also: the condition `Q >= 0` is sufficient, but not necessary.

Also: I’m away from computer but should compute eigvals(A-B\*F) dividing A, B by sqrt(beta) to undo normalization.  
I think, Will get one Eigenvalue \<1, one =1.  
I think due to dummy state variable =1

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 11, 2022, 11:42am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/27 "2022-09-11T11:42:18Z")

</div>

I confess I did not go through all the technical details of this (and the previous) posts, but I suspect that one issue here might be related to what is the goal of the optimization. In (optimal) control theory, it is conventional to formulate the problem of optimal control as **minimization of cost**. See the classical monographs by Kwakernaak & Sivan, Athans & Falb, Bryson & Ho, Kirk, Lewis, …

Indeed, this is one of the conceptual clashes between the optimal control theory as developed within control theory and reinforcement learning as originally developed within computer science, because in RL they are **maximizing the value**.

As a consequence, when the sum-of-squares type of a cost function is considered in control theory texts and software, namely,

J(x\_0, u\_0,u\_1,\ldots,u\_{N-1}) = \frac{1}{2}x\_N^\top S x\_N + \frac{1}{2} \sum\_{k=0}^{N-1} \left [x\_k^\top Q x\_k + u\_k^\top R u\_k\right]

with x\_{k+1} = f(x\_k,u\_k), \; x\_0 \; \mathrm{given}, it is necessary to impose nonnegativity conditions on the weigting matrices S, Q, R, that is,

S\succeq 0, \; Q \succeq 0, \; R \succeq 0.

The condition is strenghtened to R\succ 0 for some scenarios such as the ARE-based solution to LQ-optimal control problem (the cost function is as above but the system dynamics is linear) when no constraint on the magnitude of u\_k is imposed.

To summarize, while applying classical results originating in optimal control theory such as LQ-optimal control, you should be aware that the optimal control problem has surely been formulated as a minimization of some cost. Reformulation to maximation is straightforward, extension to saddle point optimization has also been done ( **LQ differential games** ).

---

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

---

<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 12, 2022, 9:04pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/29 "2022-09-12T21:04:18Z")

</div>

I solved the Discrete time, Deterministic, infinite horizon version w/ SDDP.jl.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/2/d274d4974489941c30fa71364e177110cbb96210.png) ![image](https://global.discourse-cdn.com/julialang/original/3X/b/1/b1050b373c9f700df7dde39590752090aedbf022.png)

```julia
using SDDP, Plots
import Ipopt

###################################################
δ, ρ = 0.1, 0.05; β = 1 / (1 + ρ);
z = ρ + δ + 0.02;
u_SS = (z - ρ - δ) / (ρ + δ);
x_SS = u_SS / δ;
x0 = 0.5 * x_SS;
x0 = 1.5 * x_SS;

# f(x,u;δ=δ) = (1 - δ)*x + u; # transition fcn 
# F(x,u;z=z) = z*x -u -0.5*u^2; # payoff fcn
opt = Ipopt.Optimizer;
T_sim=100;
#####################################################
graph = SDDP.LinearGraph(1)
SDDP.add_edge(graph, 1 => 1, β)

model = SDDP.PolicyGraph(
    graph; sense = :Max, upper_bound = 1000, optimizer = opt,
) do sp, _
    @variable(sp, x, SDDP.State, initial_value = x0) # state variable
    @variable(sp, u) # control variable
    @constraint(sp, x.out == (1 - δ) * x.in + u) # transition fcn 
    @stageobjective(sp, z * x.in - u - 0.5 * u^2) # payoff fcn 
    return
end

SDDP.train(model; iteration_limit = 10)

sims = SDDP.simulate(model, 1, [:x, :u];
    sampling_scheme = SDDP.Historical([(1, nothing) for t in 1:T_sim])
);

sim_K = map(data -> data[:x].out, sims[1])
sim_i = map(data -> data[:u], sims[1])

plot(legend=:topright, ylims = (0,2))
plot!(sim_K, c=1, lab="K")
plot!(1:T_sim, tt -> x_SS, c=1, l=:dash, lw=3, lab = "K_SS")
plot!(sim_i, c=2, lab="i")
plot!(1:T_sim, tt -> u_SS, c=2, l=:dash, lw=3, lab = "i_SS")

```

---

<div class="post-metadata">

**Author:** ![joaquimg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joaquimg/32/223_2.png) [@joaquimg](https://discourse.julialang.org/u/joaquimg)\
**Post date:** [September 13, 2022, 4:18am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/30 "2022-09-13T04:18:31Z")

</div>

SDDP was designed to be discrete+stochastic. The S in SDDP stands for Stochastic.  
In order to converge smoothly, it requires convexity and stage-independence of random variables.  
There are generalizations to handle more general cases in exchange for performance.  
Note that many practical non-convex problems with stagewise dependent random variables are solved approximately with the SDDP algorithm every day to operate multiple power systems in many countries.

---

<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 13, 2022, 6:19pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/31 "2022-09-13T18:19:24Z")

</div>

Here is the solution using [InfiniteOpt](https://github.com/infiniteopt/InfiniteOpt.jl/discussions/194).jl.  
Even though this package is intended to solve finite horizon models, we can set the Salvage value equal to the value function in the infinite horizon case (which I solved in closed form above).

```julia
using InfiniteOpt, Ipopt, Plots
δ=0.10; # capital depreciation
ρ=0.05; # discount rate
T = 90.0 # life horizon
z=ρ+δ +0.02    
# z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0)/(ρ))*D_SS;
Q_SS = (1+I_SS);
k0 = 0.5K_SS # endowment

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
S(k) = exp(-ρ*T)*(V_SS - Q_SS*K_SS + Q_SS*k)
discount(t; ρ=ρ) = exp(-ρ*t) # discount function
BC(k, i; δ=δ) = i - δ*k # LOM 

opt = Ipopt.Optimizer # desired solver
ns = 10_000; # number of gridpoints
m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount) + S(k(T)))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, c2, k >= 0) # Can't sell more than all your capital.
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, i(T) >=-1)
@constraint(m, k(T) >=0)

optimize!(m)
termination_status(m)

i_opt = value(i)
k_opt = value(k)
ts = supports(t)
opt_obj = objective_value(m) # V(k0, 0)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:topright, ylims=(0,2));
plot!(ts[ix], k_opt[ix], color = 1, lab = "k: InfiniteOpt")
plot!(ts[ix], tt -> K_SS, color = 1, lab = "K_SS", l=:dash, lw=3)
plot!(ts[ix], i_opt[ix], color = 2, lab = "i: InfiniteOpt")
plot!(ts[ix], tt -> I_SS, color = 2, lab = "I_SS", l=:dash, lw=3)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/2/0/20b543d65d52f8a42a98b3be3c7d1ddc50924de0.png)

---

<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:** [November 13, 2022, 9:45pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/32 "2022-11-13T21:45:49Z")

</div>

It would be awesome to see if the discrete-time version of this problem can be solved with @Norman package [SequenceJacobians](https://github.com/junyuan-chen/SequenceJacobians.jl).jl.

It’s a wip but I was able to figure out how to solve the simple RBC model w/ his pkg in the examples:

```julia
using SequenceJacobians, NLsolve;

# Define RBC Model ################################################################
@simple function firm(K, L, Z, α, δ)
    r = α * Z * (lag(K) / L)^(α-1) - δ
    w = (1-α) * Z * (lag(K) / L)^α
    Y = Z * lag(K)^α * L^(1-α)
    return r, w, Y
end

@simple function household(K, L, w, eis, frisch, φ, δ)
    C = (w / (φ * L^(1/frisch)))^eis
    I = K - (1 - δ) * lag(K)
    return C, I
end

@simple function mkt_clearing(r, C, Y, I, K, L, w, eis, β)
    goods_mkt = Y - C - I
    euler = C^(-1/eis) - β*(1+lead(r))*lead(C)^(-1/eis)
    walras = C + K - (1+r)*lag(K) - w*L
    return goods_mkt, euler, walras
end

rbcblocks() = (firm_blk(), household_blk(), mkt_clearing_blk())
################################################################
m = model(rbcblocks()) # vsrcs(m) == [:K, :L, :Z, :α, :δ, :eis, :frisch, :φ, :β]
calis = [:L=>1, :eis=>1, :frisch=>1, :δ=>0.025, :α=>0.11]
tars = [:goods_mkt=>0, :r=>0.01, :euler=>0, :Y=>1]
inits = [:φ=>0.9, :β=>0.99, :K=>2, :Z=>1]

# [:K init 2, :L=1, :Z init 1.0, :α=0.11, :δ=0.025, :eis=1, :frisch=1, :φ init 0.9, :β init 0.99]

ss = SteadyState(m, calis, inits, tars)
f!(y, x) = residuals!(y, ss, x)
r = solve!(NLsolve_Solver, f!, ss.inits) # r[1] # getval(ss, :K)

J = TotalJacobian(model(ss), [:Z,:K,:L], [:euler, :goods_mkt], getvarvals(ss), 300)
GJ = GEJacobian(J, :Z)
dZ = zeros(300)
dZ[11:end] .= getvarvals(ss)[:Z] .* 0.01 .* 0.8.^(0:289)
irfs = linirf(GJ, :Z=>dZ)
dCdZ = irfs[:Z][:C]
dwdZ = irfs[:Z][:w]
dYdZ = irfs[:Z][:Y]

ε = randn(10_000)
# simulate(GJ,exovar, endovar, ε, ρ, σ=1.0; kwargs...)
s = simulate(GJ, :Z, :K, ε, 0.9)
s = simulate(GJ, :Z, :K, ε, 0.9, 1.00)
#
s = simulate(GJ, :Z, :K, ε, 0.9, 0.02)
s = simulate(GJ, :Z, :K, ε, 0.9, 0.00)
s = simulate(GJ, :Z, :K, ε, 0.9, 0.10)
#
s = simulate(GJ, :Z, :K, ε, 0.2, 0.02)
s = simulate(GJ, :Z, :K, ε, 0.2, 0.00)
s = simulate(GJ, :Z, :K, ε, 0.2, 0.10)

```

---

<div class="post-metadata">

**Author:** ![Norman](https://avatars.discourse-cdn.com/v4/letter/n/97f17d/32.png) [@Norman](https://discourse.julialang.org/u/Norman)\
**Post date:** [November 13, 2022, 10:09pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/33 "2022-11-13T22:09:57Z")

</div>

Hi, @Albert_Zevelev! Thanks for your interest.

What I was experimenting with [SequenceJacobians.jl](https://github.com/junyuan-chen/SequenceJacobians.jl) is essentially a Julia implementation of the Python package [shade-econ/sequence-jacobian](https://github.com/shade-econ/sequence-jacobian) developed by the original authors who wrote the paper. You might want to read their [Econometrica paper](https://doi.org/10.3982/ECTA17434) to figure out what these packages are about. I wrote the Julia code in order to help me understand exactly how their methods work in practice (and also because I am less proficient in Python for more complicated work).

---

<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:** [November 13, 2022, 10:33pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/34 "2022-11-13T22:33:06Z")

</div>

Thanks. I’m not having luck applying it to the following very simple model (w/o equilibrium prices).

State variables: K\_{t}, Z\_{t} {capital stock & tech shock}  
Control variables: I\_{t} {firm investment}

LOM K: K\_{t+1} = I\_{t} + (1-\delta ) \* K\_{t}  
LOM Z: Z\_{t+1} = (1- \rho\_{z} )\mu\_{z} + (\rho\_{z} ) \* Z\_{t} +\sigma\_{z} \varepsilon\_{t}, \varepsilon\_{t} \sim\_{iid} N(0,1)  
EE: I\_{t+1} = \frac{\rho + \delta - Z\_{t}}{1 - \delta} + \frac{1 +\rho }{1-\delta} I\_{t}  
note: this model has a unique steady-state

---

<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:** [November 14, 2022, 7:44pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/35 "2022-11-14T19:44:38Z")

</div>

@Norman I got your pkg to work, but needed to push investment 1-period forward.  
(to be clear, the pkg no longer gives an error & solves for the correct steady state, I’m sure more things need to be corrected…)

```julia
using SequenceJacobians, NLsolve;
@simple function firm(K, δ)
    I = K - (1 - δ) * lag(K)
    return I
end
@simple function mkt_clearing(Z, I, δ, ρ)
    euler = lead(I) - ((ρ+δ-Z)/(1-δ)) - ((1+ρ)/(1-δ)) * I
    return euler 
end
rbcblocks() = (firm_blk(), mkt_clearing_blk())
m = model(rbcblocks()) # vsrcs(m) == [:K, :δ, :Z, :ρ]
calis = [:δ=>.10, :ρ=>0.05, :Z=>.17]
tars = [:euler=>0] # 
inits = [:K=>1]
ss = SteadyState(m, calis, inits, tars)
f!(y, x) = residuals!(y, ss, x)
r = solve!(NLsolve_Solver, f!, ss.inits) # r[1] # getval(ss, :K)

# julia> getvarvals(ss)
# (K = 1.33, ρ = 0.05, Z = 0.17, δ = 0.1, I = 0.1333, euler = 1.345e-7)

n_irf = 300; #n_irf = 2_000;

J = TotalJacobian(model(ss), [:K,:Z], [:euler], getvarvals(ss), n_irf)
GJ = GEJacobian(J, :Z)
dZ = zeros(n_irf)
dZ[11:end] .= getvarvals(ss)[:Z] .* 0.01 .* 0.8.^(0:(n_irf-11))
irfs = linirf(GJ, :Z=>dZ)

dKdZ = irfs[:Z][:K]
dIdZ = irfs[:Z][:I]

using Plots; 
plot(dKdZ); plot!(dIdZ)

##############################
#SIM 
##############################
n_sim = 10_000;
ε = randn(n_sim)
#ε = zeros(n_sim)
# simulate(GJ,exovar, endovar, ε, ρ, σ=1.0; kwargs...)
plot()
simulate(GJ, :Z, :K, ε, 0.9) |> plot! 
simulate(GJ, :Z, :K, ε, 0.9, 1.00) |> plot! #verify default value σ=1
#
ρ_Z = 0.01; σ_Z = 1e-4;
plot()
simulate(GJ, :Z, :K, ε, ρ_Z, 9*σ_Z) |> plot!
simulate(GJ, :Z, :K, ε, ρ_Z, 2*σ_Z) |> plot!
simulate(GJ, :Z, :K, ε, ρ_Z, 0*σ_Z) |> plot!
#
ρ_Z = 0.99; σ_Z = 1e-4;
plot()
simulate(GJ, :Z, :K, ε, ρ_Z, 9*σ_Z) |> plot!
simulate(GJ, :Z, :K, ε, ρ_Z, 2*σ_Z) |> plot!
simulate(GJ, :Z, :K, ε, ρ_Z, 0*σ_Z) |> plot!

```

Some questions:  
Q1: Is it possible to use your package to simulate starting from K\_{0} = 1.5\ \times K\_{ss}?  
So I get something like the left two quadrants:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/f/1f555777cbe32925b05d580627a5ed98a4d24d88.png)

Q2: is it possible to use your package to plot the policy function for the firm’s optimal investment decision?  
Q3: plot the value function?

PS: I found other SSJ Julia implementations [here](https://github.com/jprodriguesumn/SequenceSpaceJacobians)

---

<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 2, 2022, 7:16pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/36 "2022-12-02T19:16:19Z")

</div>

I forgot to include:  
[GDSGE](https://www.gdsge.com/) (in C++ & Matlab) seems like an impressive package for global solutions (using PFI).

---

<div class="post-metadata">

**Author:** ![thorek1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thorek1/32/20989_2.png) [@thorek1](https://discourse.julialang.org/u/thorek1)\
**Post date:** [March 5, 2023, 9:24pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/37 "2023-03-05T21:24:23Z")

</div>

Hi,

I have another package to add which solves discrete time dynamic optimisation problems with perturbation methods: [MacroModelling.jl](https://github.com/thorek1/MacroModelling.jl)

Solving the firm investment problem in it’s deterministic version would look like this:

```julia
using MacroModelling

@model firm_investment_problem begin
    K[0] = (1 - δ) * K[-1] + I[0]
    Z[0] = (1 - ρ) * μ + ρ * Z[-1] 
    I[1] = ((ρ + δ - Z[0])/(1 - δ)) + ((1 + ρ)/(1 - δ)) * I[0]
end

@parameters firm_investment_problem begin
    ρ = 0.05
    δ = 0.10
    μ = .17
    σ = .2
end

steady_state = get_steady_state(firm_investment_problem)

plot(firm_investment_problem, initial_state = collect(steady_state(:,:Steady_state)) .* [1, 1.5, 1])

```

 ![irf __firm_investment_problem__ no_shock__1](https://global.discourse-cdn.com/julialang/original/3X/b/2/b2b67d7b580cb3cc59dd98f016ab2aa0d7c07a4f.png)  
The `@model` macro is used to define the model equations and parameters are declared in the `@parameters` macro. The package then takes care of variable and parameter declaration and solves the steady state symbolically (if possible and numerically otherwise) for you.

The `get_steady_state` function returns the steady state and its derivatives wrt model parameters and the `plot` function with the `initial_state` for `K` set 50% higher delivers the simulations seen above.

There are about 13 DSGE models (some as large as 100 equations) already implemented using this package. The package has a lot of standard output implemented (autocorrelation, forecast error variance decomposition, variance, …) and can solve models up to third order using perturbation methods. Check out the [documentation](https://thorek1.github.io/MacroModelling.jl/stable) for more details.

---

<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:** [April 17, 2024, 6:47am UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/38 "2024-04-17T06:47:13Z")

</div>

There is a new package:

> **[GitHub - control-toolbox/OptimalControl.jl: Solvers of optimal control problems](https://github.com/control-toolbox/OptimalControl.jl)**
>
> Solvers of optimal control problems. Contribute to control-toolbox/OptimalControl.jl development by creating an account on GitHub.

Presented at last year’s JuliaCon:

[![](https://global.discourse-cdn.com/julialang/original/3X/d/d/dd6f11f65c23431d33fcc20126c73a4ac27e4ffa.jpeg "On solving optimal control problems with Julia | Caillau, Cots, Gergaud, Martinon | JuliaCon 2023") ](https://www.youtube.com/watch?v=RYUtVnzLj5k)

---

<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:** [April 22, 2024, 10:54pm UTC](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285/39 "2024-04-22T22:54:24Z")

</div>

New package for MDPs in Julia:

> **[GitHub - Zinoex/IntervalMDP.jl: GPU-accelerated value iteration for Interval...](https://github.com/Zinoex/IntervalMDP.jl)**
>
> GPU-accelerated value iteration for Interval Markov Decision Processes - Zinoex/IntervalMDP.jl

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