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…
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> 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?
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
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.
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.
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. Crandall and Lions (1983): if P \Rightarrow HJB equation has unique viscosity solution
Lions won the fields medal for this & other work in 1994. Finite difference 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 methods are the most fruitful for these types of HJBs.
A new generation of solvers use deep-learning methods & claim to solve HJBs of >100 dimensions. I have not seen them successfully solve this type of Lifecycle economics problem yet…
One way NeuralPDE.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.
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.
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.
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.
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)
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.
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.