EconPDEs.jl

EconPDEs.jl solves nonlinear ODEs and PDEs that arise in economic models, especially Hamilton-Jacobi-Bellman equations.

You specify three things: (i) the state space, (ii) an initial guess for the solution, and (iii) a function that returns the PDE residual at each grid point. The package handles the rest — it builds the finite-difference derivatives, applies upwind stencils, assembles the sparse Jacobian, and solves the system by pseudo-transient Newton iteration.

Use it for stationary or time-dependent HJBs on one or more state variables, with one or more value functions.

EconPDEs.jl is robust, fast, and solves virtually all textbook continuous-time models; see the Examples for a list.

Installation

The package is registered in the Julia General registry:

] add EconPDEs

Current versions require Julia 1.10 or later.

Quickstart

Here is a small linear equation on one state, with a mean-reverting drift:

$$\rho v(x) = x + \mu(x) v_x(x), \qquad \mu(x) = -\kappa (x - 0.5).$$

pdesolve finds the stationary solution by embedding the equation in a time-dependent problem, solved backward in time from an initial guess until the time derivative is zero.

$$\rho v(x, t) = x + \mu(x) v_x(x, t) + v_t(x, t),$$

Here is how to solve this equation using this package:

using EconPDEs

const rho = 0.05
const kappa = 0.2

# 1. Define a `NamedTuple` representing the state space as a grid:
grid = (; x = range(0.0, 1.0, length = 50))

# 2. Define a `NamedTuple` representing an initial guess for the solution:
guess = (; v = zeros(length(grid.x)))

# 3. Define a function  encoding the PDE at each grid point:
function pde(state::NamedTuple, u::NamedTuple)
    # `state` holds the coordinates of the current grid point
    (; x) = state
    # `u` holds the solution and its finite-difference derivatives at that point.
    (; v, vx_up, vx_down, vxx) = u
    # The function must return the time derivative implied by rho*v = x + mu*v_x + v_t.
    mu = -kappa * (x - 0.5)
    vt = -(x + mu * (mu >= 0 ? vx_up : vx_down) - rho * v)
    return (; vt)::NamedTuple
end

# 4. Pass these three objects to `pdesolve`
result = pdesolve(pde, grid, guess)
# Solving for 1 unknown (v) on a 50-point grid
# Iter TimeStep Residual
# ---- -------- --------
#    1 1.00e+00 7.51e-02
#    2 1.09e+01 4.45e-02
#    3 1.85e+03 4.67e-04
#    4 1.76e+08 1.22e-10
# Converged after 4 time steps (63ms)
# EconPDEResult
#   solution:      v (50)
#   residual_norm: 1.22e-10
#   converged:     true (tolerance 1.49e-08)

Note that, in encoding the pde, we used the forward difference when the drift is positive, the backward one otherwise (upwinding).

The same interface extends to multiple states...

grid = (; x = range(0.0, 1.0, length = 50),
          y = collect(range(0.5, 1.5, length = 10)))
guess = (; v = zeros(length(grid.x), length(grid.y)))

function pde(state, u)
    (; x, y) = state
    (; v, vx_up, vx_down, vy_up, vy_down, vxx, vyy, vxy_up, vxy_down) = u
    ...
    return (; vt)
end

... as well as multiple value functions (i.e., coupled PDES).

grid = (; x = range(0.0, 1.0, length = 50))
guess = (; v = zeros(length(grid.x)), w = zeros(length(grid.x)))

function pde(state, u)
    (; x) = state
    (; v, vx_up, vx_down, vxx, w, wx_up, wx_down, wxx) = u
    ...
    return (; vt, wt)
end

Documentation

Read the documentation for further details. The documentation also includes a gallery of runnable examples (see examples/ for the original scripts):

• consumption-saving models, including borrowing constraints, income risk, and risky assets;
• asset-pricing models, including habit, long-run risk, rare disasters, and intermediary finance;
• corporate-finance models, including optimal default and investment with financing constraints.

Citation

Please cite EconPDEs.jl if it is useful in your work. You can use the "Cite this repository" button on the GitHub page, or the following BibTeX entry:

@software{EconPDEs.jl,
  author = {Gomez, Matthieu},
  title  = {{EconPDEs.jl}: {Solving PDEs in economics}},
  url    = {https://github.com/matthieugomez/EconPDEs.jl},
}