Solving a PDE which comes from an HJB in DifferentialEquations.jl

I was wondering whether one could use the Finite Difference Method developed in the recent GSoC project DiffEqOperators.jl to solve the following PDE, which comes from an HJB that describes a household consumption/saving problem from economics:

where a is wealth and y is income, see section 4 in this paper on page 43 for details (boundary conditions, etc).

I thought about discretizing the y-variable and solving only in the a-direction.

The regular way to solve those PDEs is writing your own FDM or use a specialized FDM solver package.

1 Like

In principle there’s a way to rewrite it as just a derivative in a, and time step a system of ODEs in the a direction. You need to make sure that the operators you’re choosing give a stable discretization of the PDE though. I am not an expert in HJB itself so maybe @jlperla could chime in.

While there is some work to be done to get the BCs to all line up to the same convention:

the Neuman0 BC is exactly what you need here. So a lot of those operators should be pretty easy to write down with the lazy operators of DiffEqOperators.jl. The tricky thing is that da isn’t isolated. You’ll probably need to chain rule on the da term to isolate dv/da, but then you’ll get something that discretizes nicely. This source

says it’s possible to get there though I don’t know the exactly how your formulation lines up with their discussion there.

(Though note I am looking for collaborators to help build an HJB solver for Julia so this can be put to rest since it comes up quite often. Or if you’re just willing to throw papers at me to let me know what economists are interested in, I’ll read them.)

@MaximilianJHuber As you have figured out, these HJBE have their own solution methods. At this point I do not believe you can use “standard” solvers could work and that the DiffEqOperators.jl is exactly the right way to help solve them (as you have figured out). There are a few tricks that would need to be added to it, in particular upwind procedures for multiple dimensions, etc.

I would scan and if you haen’t already. They go through the code for that paper you referenced in detail.

Also, I have had grad students working on notes for the discretization of the operators which you can look at and in particular

I have been talking to @ChrisRackauckas about trying to get the differential operators used for solving upwind HJBE, but I am overwhlemed with work and have failed him. Send me an email and maybe we can help coordinate a group of people interested in Julia continuous time methods for economics/finance.


Only a relatively small number of economists use continuous time methods, and an even smaller fraction of them solves these models numerically. It is not taught in most graduate programs.

When the functions are sufficiently smooth, a common choice is collocation, eg with complete polynomials or splines. You need a good initial guess for convergence, linear or higher order perturbations are usually used for this. “Applied Computational Economics and Finance” by Miranda and Fackler has the most detailed treatment of these problems I have seen in economics, could be a good starting point.

Discretizations are the last resort, because they converge slower and it is tricky to get them to be accurate. A very robust set of methods were developed by Kushner and Dupuis (see their textbook by Springer), this basically converts into a discrete time problem with sparse & local transition laws for Ito processes, and you get the usual benefits (contraction, etc).

Value functions can have high curvature in some regions (my understanding is that this is what ODE people call “stiff”, but I may be mistaken). Tricks and transformations can cope with this, and also with non-smooth regions and similar. More of an art than science, as @jlperla said, there are no “standard” methods.

1 Like

This is changing rapidly, thankfully! With a little push, I am pretty sure economists are going to have idiot-proof toolboxes to solve them. There is a large group of us eager for this to happen.

Yes, the upwind methods that are referenced in the Moll notes are very related to this.

I disagree a little here. I think what I meant to say is that you can’t use standard PDE solvers. But there are standard methods (i.e. the upwind finite differences). Once you battle through the math, they turn out to be relatively straightforward to implement, and have great convergence properties due to the connection between upwind methods and viscosity solutions (i.e. barles souganidis montone sequences).

Also, optimal stopping problems are WAY easier in continuous time since you can write things as LCP methods (see the links on that github page I sent or on Ben Moll’s notes).

I previously had done all my work with collocation solutions to these sorts of problems but you run into all sorts of convergence issues (even without kinks) and you have the perpetual issue of dealing with transervality conditions. Time spent on upwind finite differences is well spent, and I have faith that there will be some great libraries to solve these things automagically (e.g. is a starting point)

1 Like

This is interesting, I will read the references there.

The most tricky part of these problems for me is not the partial equilibrium, but getting a GE (but not necessarily HA) model to converge. Do the benefits you mention still apply? Are there any examples you would recommend reading/working through (or implementing in Julia :smile:)?

Check out and mostly
the ben moll links. This is a HUGE area to work on for macro and
macro-finance, and I really hope you get interested. Julia is the perfect
language for this.

The information on this stuff has not entirely diffused yet, but there are
a lot of people interested in it… mostly concentrated at Chicago,
Columbia, Princeton, Penn, NYU, WashU, and (now) UBC.

I am very busy for the next few weeks, but if you are interested in this I
would love to talk about it.


We do have arbitrary order Upwind operators in DiffEqOperators.jl, though as listed above we do need to go back through and check consistency of the BC conventions.

I’d be careful there. Depends on what you mean by standard PDE solver. You definitely can put a discretization of this into a DAE solver like IDA, though it will look really odd since the “timestepping” operator is reversed in some parts of space for a due to upwinding. And with the right massaging this seems like it can fit into an ODE solver with the right mass matrix (which is likely time dependent). It’s not simple and will take a lot of thought, but I am pretty sure there’s a way to get a 4th-5th adaptive discretization out by combining the right tools with high order discretization operators.

Yeah, let’s get in on this early and build the right software. We’re already close.

Great. I would love it if a DAE solver could work to solve this stuff, but I am not sure I see it.

Check out in the meantime to see someone who did a first take on doing this sort of library.

The other big thing is that the discretized operators can be used for solving optimal stopping problems, which ends up mapping into complementarity problems. For me in the medium-term, this is going to be very important. Check out Option Pricing and Linear Complementarity I am not sure how a DAE can work there, but it could be

Finally, for me I will need to get the diffeqoperators to support jump diffusions in the medium-term for a paper. I am hiring an RA who is very competent and may be able to work to add that kind of thing given some support. More later.

Yeah, I think finishing up the operators is the best short-term goal. A lot of people seem interested in them for many different purposes, and it would be a good service to get those done, stabilized, efficient, and benchmarked.

It can be done in EconPDEs. For reference, I have uploaded this PDE in the examples folder.


Hi, I am having trouble plotting graphs in Julia. I tried to plot the DiTella model found in EconPDEs examples folder using:

plot results

using Plots
surface(grid[:x], grid[:ν], result[:p])

and got the following error message:

MethodError: no method matching getindex(::typeof(grid), ::Symbol)

[1] top-level scope at In[28]:3

Any help would be greatly appreciated. Thanks

Please provide a mwe (minimal working example) - what is grid here? In any case the error message is implying that that is not the right way of extracting fields from the grid object.

Avoid calling any object grid in your userspace, because grid is a function in Plots. So replace the name of your object grid by stategrid for instance.

Thank you both for the responses. The comment from 2LxtkNuVTZof](/u/2LxtkNuVTZof worked!

Since this popped up again, I would like to point people to our large dimensional neural network based PDE solver that is in NeuralNetDiffEq.jl. A 100-dimensional HJB is one of the examples.