# Solve a routine economics PDE from an HJB w/ ModelingToolkit
**URL:** https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718
**Category:** Modelling & Simulations
**Tags:** economics, pde, finance
**Created:** [November 7, 2020, 2:46am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718 "2020-11-07T02:46:03Z")
**Posts on this page:** 14
**Page:** 2
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [November 12, 2020, 10:21pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/21 "2020-11-12T22:21:01Z")
I can further simplify the PDE: \rho V(t,a\_{t}) = \frac{\gamma \times V\_{a}^{1-\frac{1}{\gamma} } }{1-\gamma }+ V\_a \times (y + ra\_{t}) + V\_{t}
W/ \gamma=2: \rho V(t,a\_{t}) = -2 V\_{a}^{\frac{1}{2}}+ V\_a \times (y + ra\_{t}) + V\_{t}
Issue: V\_{a}^{\frac{1}{2} } is not real when V\_a \<0.
Trick: define the piecewise function: f(V\_a) = 1\left\{ V\_a \>0 \right\} \left( -2 V\_{a}^{\frac{1}{2}} \right)
Problem: ModelingToolkit.jl doesn’t like these types of piecewise functions
-I prob don’t know the correct syntax…
```julia
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux, Plots
@parameters t a θ
@variables V(..)
@derivatives Da'~a
@derivatives Dt'~t
γ = 2.00; ρ = 0.01; r = 0.01; y = 0.00; T = 5.0; ω = r - (r-ρ)/γ;
u(c,γ)= (c^(1.0 -γ))/(1.0 -γ);
VS(t,a;y=y,r=r,γ=γ,ω=ω,T=T) = u(a+(y/r), γ) * ( (1.0 - exp(-ω*(T-t))) / ω )^γ;
v = V(t,a,θ);
va = Da(V(t,a,θ));
vt = Dt(V(t,a,θ));
f(va) = va > 0.0 ? γ*(va^(1.0 - 1.0/γ) )/(1.0 - γ) : 0.0;
eq = ρ*v ~ f(va) + va*(y + r*a) + vt
```
Error message
```julia
julia> f(va)
ERROR: TypeError: non-boolean (Operation) used in boolean context
Stacktrace:
[1] f(::Operation) at .\untitled-7e85d02f47feb4021c0dd3ed76113a34:22
[2] top-level scope at none:1
```
@ChrisRackauckas is there a way to include this type of piecewise function of a derivative?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 12, 2020, 10:30pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/22 "2020-11-12T22:30:32Z")
Use `IfElse.ifelse(x,y,z)`.
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [November 12, 2020, 10:51pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/23 "2020-11-12T22:51:15Z")
I get an error:
```julia
julia> using IfElse
julia> f(x)=IfElse.ifelse(x > 0.0, x^(.5), -1)
f (generic function with 1 method)
julia> f(10)
3.1622776601683795
julia> f(-10)
ERROR: DomainError with -10.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(::Float64) at .\math.jl:37
[2] ^ at .\math.jl:888 [inlined]
[3] ^ at .\promotion.jl:343 [inlined]
[4] f(::Int64) at .\none:1
[5] top-level scope at none:1
```
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [November 12, 2020, 11:13pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/24 "2020-11-12T23:13:50Z")
Kenneth Judd has some slides online about Chebyshev colocation methods that has a continuous time consumption savings problem as one of the examples, at least that’s what my verrrry distant memory tells me.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 13, 2020, 6:16am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/25 "2020-11-13T06:16:51Z")
It’s a function so it evaluates before the call.
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [November 13, 2020, 7:02am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/26 "2020-11-13T07:02:00Z")
Update:
1. Recall the value function to the optimization problem:
V(t,a\_{t}) := \max \int\_{\tau=t}^{\tau = T} e^{-\rho (\tau -t)} u(c\_{\tau})d\tau \text{ s.t. }
\dot{a}\_{t} = y + ra\_{t} - c\_{t},
a\_{0} \text{ given, } a\_{T}=0.
2. It turns out there are multiple solutions to the HJB-PDE:
\left[\begin{array}{l} \rho V(t,a\_{t}) = \frac{c\_{t}^{1-\gamma}}{1-\gamma} + V\_{a}(t,a\_{t})\times \left(y + ra\_{t} - c\_{t} \right) + V\_{t}(t,a\_{t}) \\ c(t,a\_{t}) = (V\_{a}(t,a\_{t}))^{-\frac{1}{\gamma}} \\ V(T,a\_{T}) = 0 \end{array} \right]
Example: if \gamma \>1, V(t,a\_{t})=0 satisfies the PDE & boundary.
Only one of the solutions to this PDE is the same solution to the optimization problem.
This unique “nice” solution is the [viscosity solution](http://www.princeton.edu/~moll/viscosity_slides.pdf).
**Crandall and Lions** (1983): if P \Rightarrow HJB equation has unique viscosity solution
Lions won the [fields medal](https://en.wikipedia.org/wiki/Pierre-Louis_Lions) for this & other work in 1994.
[Finite difference](http://www.princeton.edu/~moll/ECO521_2016/Lecture3_ECO521.pdf) scheme “converges” to the unique viscosity solution under three conditions: monotonicity, consistency, stability.
Bottom line: V(T,a\_{T}) = 0 is not enough to pin down the “nice” solution, we also need an appropriate boundary inequality.
I’ve looked around & it appears that [finite difference](https://benjaminmoll.com/codes/) methods are the most fruitful for these types of HJBs.
A [new generation of solvers use deep-learning methods](https://github.com/frankhan91/DeepBSDE) & claim to solve HJBs of \>100 dimensions. I have not seen them successfully solve this type of Lifecycle economics problem yet…
One way [NeuralPDE](https://github.com/SciML/NeuralPDE.jl).jl can find the “nice” solution is if they allowed the user to restrict the solution V(t,a\_{t}) w/: V\_a \>0 and V\_{aa}\<0. That is monotonic and concave neural networks.
---
### Author:  [@jlperla](https://discourse.julialang.org/u/jlperla)
#### Post date: [November 13, 2020, 5:09pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/27 "2020-11-13T17:09:12Z")
> [@Albert\_Zevelev](#):
>
> A [new generation of solvers use deep-learning methods](https://github.com/frankhan91/DeepBSDE) & claim to solve HJBs of \>100 dimensions. I have not seen them successfully solve this type of Lifecycle economics problem yet…
It isn’t even close. There is no implementation that I am currently aware of that uses these techniques for solving even a low dimensional HJBE that actually has an “H” in it induced by a standard continuous control. There is some work on problems with the free-boundary induced by options which can implement optimal stopping problems, but it is very complicated. I think that before jumping into this sort of research you probably should replicate the state-of-the-art algorithms there (and look to see if there is any progress there - probably in the control literature). Regardless, it is going to be very involved and has a low likelihood of beating the highly-tuned and carefully studied upwind finite-difference methods that people use in control for problems of this size.
I think the modelingtoolkit is a red-herring here. Its purpose is to do tracing to generate jacobians/etc. for complicated functions, but the thing you are posing is simple enough that you don’t need it and could handcode the jacobians by hand if you chose to. Not to say you should, but I don’t think that is the core of the issue here.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 14, 2020, 12:16am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/28 "2020-11-14T00:16:27Z")
> [@jlperla](#):
>
> It isn’t even close. There is no implementation that I am currently aware of that uses these techniques for solving even a low dimensional HJBE that actually has an “H” in it induced by a standard continuous control. There is some work on problems with the free-boundary induced by options which can implement optimal stopping problems, but it is very complicated
Indeed, and there are improvements to these methods in the Julia NeuralPDE.jl package
[https://neuralpde.sciml.ai/dev/examples/100\_HJB/](https://neuralpde.sciml.ai/dev/examples/100_HJB/)
But they don’t quite solve the macroeconomics problems yet. While still HJBEs, it doesn’t have all of the add-ons that are usually tagged on in those applications. But I’m working (on the side, not as a main project, so no timeline) with @jlperla to fix that. It’ll take some time though. I suspect they will be efficient here in the end (unlike general PINNs which don’t specialize enough), but there’s work to be done to get there.
---
### Author:  [@Honza9723](https://discourse.julialang.org/u/Honza9723)
#### Post date: [November 14, 2020, 12:56am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/29 "2020-11-14T00:56:51Z")
What about this?
> **[Deep\_Learning.pdf](https://maximilianvogler.github.io/My_Website/Deep_Learning.pdf)**
>
> 980.32 KB
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [November 14, 2020, 2:13am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/30 "2020-11-14T02:13:19Z")
@Honza9723 yep.
Related Matlab code is here: [https://github.com/jesusfv/financial-frictions](https://github.com/jesusfv/financial-frictions)
And [MLEcon](https://github.com/vduarte/MLEcon/blob/master/examples/DSGE%20-%20Caldara%20et%20al%202012.py).py claims to solve continuous time DSGE models w/ TensorFlow…
It would be great if the FD methods in Moll et al were automated in Julia. (possibly one day in ModelingToolkit?) This would allow side-by-side comparisons w/ Deep NN methods (when they work for these problems on Julia)…
EconPDEs.jl does this for up to two state variables.
DSGE.jl has `solve_hjb()` & `solve_kfe()` but they don’t work well currently. (They even have an example of a continuous time Krusell Smith model…)
There are other Julia packages for solving HJBs (HJBSolver.jl) which are currently not maintained.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 14, 2020, 3:47am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/31 "2020-11-14T03:47:46Z")
> [@Albert\_Zevelev](#):
>
> It would be great if the FD methods in Moll et al were automated in Julia. (possibly one day in ModelingToolkit?) This would allow side-by-side comparisons w/ Deep NN methods (when they work for these problems on Julia)…
That’s what we’re doing with DiffEqOperators.jl. It’s not there yet though.
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [November 17, 2020, 7:44am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/32 "2020-11-17T07:44:53Z")
It sounds like SciML is gonna be a an incredible continuous-time modeling package.
Many things I’ve ungraciously asked for seem to be in progress:
1. boundary value DAEs (like Fortran’s Coldae)
btw neither Matlab nor Mathematica have this…
2. [optimal control](https://github.com/SciML/OptimalControl.jl/issues/4)
3. now FD methods to solve HJBs (DiffEqOperators.jl)
This is Julian greed at its best.
I really hope SciML has the funding & labor to make these dreams a reality!
(yes I realize the things I listed are prob not urgent priorities at the moment)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 17, 2020, 12:34pm UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/33 "2020-11-17T12:34:16Z")
Optimal control is coming really soon. It’s being tackled through MTK. Boundary value DAEs and more HJB stuff is not as high of a probability, but BVP DAEs I think will probably come first and hopefully as one of our summers.
---
### Author:  [@ptoche](https://discourse.julialang.org/u/ptoche)
#### Post date: [April 21, 2021, 5:29am UTC](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718/34 "2021-04-21T05:29:10Z")
Hi Albert, just came across this thread, a little late to the party! Have you made progress with this? I have always found it more intuitive to solve for the marginal utility/value, that is differentiate the PDE once to get rid of V. I’m not sure if it makes a difference to the solver, but your V\_a will always be positive, so you can force that condition (unlike with V which changes sign). If you’re looking for solutions, rather than trying to be general about the PDE, you can also transform and solve for the inverse A(c) (wealth as a function of consumption) or solve for A(V\_a) (wealth as a function of the marginal utility of wealth). It’s all in this [1986 article](https://pubsonline.informs.org/doi/abs/10.1287/moor.11.2.261).
[Previous page](https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718.md?page=1)