This package provides the function pdesolve that solves (system of) nonlinear ODEs/PDEs arising in economic models (i.e. parabolic/elliptic PDEs arising from HJB equations). It is:

• robust: upwinding + fully implicit time stepping (see here)
• fast: sparse matrices + Newton acceleration
• simple-to-use

Examples

The examples folder solves a variety of macro-finance models:

• Habit Model (Campbell Cochrane (1999) and Wachter (2005))
• Long Run Risk Model (Bansal Yaron (2004))
• Disaster Model (Wachter (2013))
• Heterogeneous Agent Models (Garleanu Panageas (2015), Di Tella (2017), Haddad (2015))
• Consumption with Borrowing Constraint (Wang Wang Yang (2016), Achdou Han Lasry Lions Moll (2018))
• Investment with Borrowing Constraint (Bolton Chen Wang (2009))

A Simple Example

For instance, to solve the PDE giving the price-dividend ratio in the Long Run Risk model with time-varying drift:

using EconPDEs

# Define a discretized state space
# An OrderedDict in which each key corresponds to a dimension of the state space.
stategrid = OrderedDict(:μ => range(-0.05, stop = 0.1, length = 500))

# Define an initial guess for the value functions
# An OrderedDict in which each key corresponds to a value function to solve for, 
# specified as an array with the same dimension as the state space
y0 = OrderedDict(:V => ones(500))

# Define a function that encodes the PDE. 
# The function takes three arguments:
# 1. A named tuple giving the current value of the state. 
# 2. A named tuple giving the value function(s) (as well as its derivatives)
# at the current value of the state. 
# 3. (Optional) Current time t
# It returns two tuples:
# 1. a tuple with the opposite of the time derivative of each value function
# 2. a tuple with the drift of each state variable (internally used for upwinding)
function f(state::NamedTuple, sol::NamedTuple)
	μbar = 0.018 ; ϑ = 0.00073 ; θμ = 0.252 ; νμ = 0.528 ; ρ = 0.025 ; ψ = 1.5 ; γ = 7.5
	Vt = 1  - ρ * sol.V + (1 - 1 / ψ) * (state.μ - 0.5 * γ * ϑ) * sol.V + θμ * (μbar - state.μ) * sol.Vμ +
	0.5 * νμ^2 * ϑ * sol.Vμμ  + 0.5 * (1 / ψ - γ) / (1- 1 / ψ) * νμ^2 *  ϑ * sol.Vμ^2/sol.V
	(Vt,), (θμ * (μbar - state.μ),)
end

# The function `pdesolve` takes four arguments:
# 1. the function encoding the PDE
# 2. the discretized state space
# 3. the initial guess for the value functions
# 4. a time grid with decreasing values 
@time pdesolve(f, stategrid, y0, range(1000, stop = 0, length = 100))
#> 0.220390 seconds (3.07 M allocations: 219.883 MiB, 18.28% gc time)

# To solve directly for the stationary solution, 
# i.e. the solution of the PDE with ∂tV = 0,
# simply omit the time grid
@time pdesolve(f, stategrid, y0)
#>  0.018544 seconds (301.91 k allocations: 20.860 MiB)

More complicated ODEs / PDES (including PDE with two state variables or systems of multiple PDEs) can be found in the examples folder.

Boundary Conditions

When solving a PDE using a finite scheme approach, one needs to specify the value of the solution outside the grid ("ghost node") to construct the second derivative and, in some cases, the first derivative at the boundary.

By default, the values at the ghost node is assumed to equal the value at the boundary node (reflecting boundaries). Specify different values for values at the ghost node using the option bc (see BoltonChenWang.jl for an example).

Installation

The package is registered in the General registry and so can be installed at the REPL with ] add EconPDEs.