Analyzing the monotonicity of numerical schemes with examples

Tianhao Zhao

Abstract

The important paper by Barles and Souganidis (1991) provides a powerful framework to analyze a class of 2nd-order non-linear PDEs which are commonly used in economics. Among all the required conditions, the monotonicity of the numerical scheme is critical and often violated in practice. This blog post explains how to analyze the monotonicity of a wide class of numerical schemes with multiple examples that share similar structures with some commonly used economics models. One can follow the examples to design their own numerical schemes.

Notation

A class of 2nd-order non-linear PDE

Before discussing the monotonicity of numerical schemes, let’s set up the problem and introduce the notations. Let’s consider a generic non-linear PDE that is abstracted from an infinite-horizon stochastic optimal control problem with discounting:

where:

• $x=(x^1,\dots ,x^k{)}^{\prime }\in {\mathbb{R}}^k$$x = (x^1,\dots,x^k)' \in \R^k$ is the vector of state variables where $x^j$$x^j$ represents the $j$$j$-th dimension element of the vector

• ${\mu }_j:{\mathbb{R}}^k\times {\mathbb{R}}^{1\times k}\times {\mathbb{R}}^{k\times k}\to \mathbb{R}$$\mu_j:\R^k\times \R^{1\times k}\times \R^{k\times k}\to \R$ is the drift function

• ${\sigma }_{mj}\mapsto \mathbb{R}$$\sigma_{mj}\mapsto \R$ is the volatility function that is an element of the volatility matrix

• $Dv,D^{2}v$$Dv,D^2 v$ are the Jacobian and Hessian of function $v$$v$

• $f\mapsto \mathbb{R}$$f\mapsto \R$ is the flow payoff/utility function

• $\rho >0$$\rho>0$ is the discounting rate

• $\mathcal{W}$$\mathcal{W}$ is a 1-dimensional Brownian motion

• $\mathcal{B}$$\mathscr{B}$ is a collection of boundary conditions as functionals

Remark:

• A generic volatility matrix is of size $k\times m$$k\times m$ in which there are $k$$k$ states and a $m$$m$-dimensional Brownian motion. WLOG, the above formulation assumes a full-size volatility matrix such that the dimension of the Brownian motion is $k$$k$ which is exactly the dimension of the state variable $x$$x$. Neither “wide” ($m>k$$m>k$) or “long” ($m<k$$m<k$) shape of the volatility matrix does not change the conclusion.

• If $v$$v$ is a functional and $x$$x$ includes infinite-dimensional objects (e.g. distribution function), then the conclusion of this post still holds. The above formulation can be seen as a finite-moment approximation of such problems.

In the left part of this post, the long function arguments are ignored if they don’t affect the interpretation. Meanwhile, for convenience, let’s define the following clamping functions:

Discretization & grid

To numerically solve this PDE, we need to discretize it over a hyper rectangle $\mathcal{X}:={\otimes }_{j=1}^k[{\underset{―}{x}}^j,{\overline{x}}^j]$$\mathcal{X} := \otimes_{j=1}^k [\underline{x}^j,\bar{x}^j]$. In the case of linear state constraints (such as borrowing constraints), one can always rectanglize the computational domain by redefining $x$$x$. We don’t discuss the case of state constraints that cannot be rectanglized.

A standard even-spaced dense grid can be defined by the number of grid points along each dimension: $\mathcal{N}:=(N^1,\dots ,N^k)$$\mathcal{N} := (N^1,\dots,N^k)$. The total number of grid points is $N:=\prod _{j=1}^k{N}^j$$N:=\prod_{j=1}^k N^j$. To index a specific grid point, we use multi-index $i:=(i^1,\dots ,i^k)$$i := (i^1,\dots,i^k)$ where $i^j=1,\dots ,N^j$$i^j = 1,\dots, N^j$ for each $j=1,\dots ,k$$j=1,\dots,k$. Let $e^j\in {\mathbb{Z}}^k$$e^j\in\Z^{k}$ be a unit vector in which all elements are zero except the $j$$j$-th element being 1. And let ${\mathrm{\Delta }}^j:=\frac{{\overline{v}}^j-{\underset{―}{v}}^j}{N^j-1}$$\Delta^j := \frac{\bar{v}^j-\underline{v}^j}{N^j-1}$ be the grid mesh size of dimension $j$$j$. Thus, the first right neighbor along dimension $j$$j$ of the current grid point $x[i]$$x[i]$ can be denoted as $x[i+{\mathrm{\Delta }}^j{e}^j]$$x[i+\Delta^j e^j]$.

For a more generic type of grid representation, we use a set of finite $N$$N$ discrete supporting nodes ${\hat{\mathcal{X}}}_N\subset \mathcal{X}$$\hat{\mathcal{X}}_N \subset \mathcal{X}$. We will discuss the generic interpolation supported by ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$.

Stacking

To well-define the linear algebra we will use later, we need an order of the elements in ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$. For even-spaced dense grid $(\mathcal{X},\mathcal{N})$$(\mathcal{X},\mathcal{N})$, the stacking is done by Cartesian product. This leads to a lexicographical sorting of the multi-indices $\{i\}$$\{i\}$. A vector of the sorted multi-indices is denoted by ${\mathcal{I}}_N\in {\mathbb{R}}^N$$\mathcal{I}_N\in\R^N$. For the more generic type of grid, one can always determine an order of the stacking. We use the same notation ${\mathcal{I}}_N$$\mathcal{I}_N$.

Remark: an alternative way of stacking representation is to use a vector of supporting nodes rather than “set+sorting”. The reason of the split here is to: 1. represent neighbor nodes in the space; 2. accommodate the fact that users have to manually define the sorting and keep it consistent in their computer programs.

In general, a grid (structure) can be represented by a 2-tuple $({\hat{\mathcal{X}}}_N,{\mathcal{I}}_N)$$(\hat{\mathcal{X}}_N,\mathcal{I}_N)$ of supporting node set and the stacking order.

Similarly, the stacking of function value $v(x[i])$$v(x[i])$ at grid nodes $x[i]\in {\hat{\mathcal{X}}}_N$$x[i]\in\hat{\mathcal{X}}_N$ sorted by ${\mathcal{I}}_N$$\mathcal{I}_N$ can be represented by $({\mathcal{V}}_N,{\mathcal{I}}_N)$$(\mathcal{V}_N,\mathcal{I}_N)$, or a vector $\mathbf{V}\in {\mathbb{R}}^N$$\mathbf{V}\in\R^N$.

Interpolation & stencil

An interpolation of function $v(x):{\mathbb{R}}^k\to \mathbb{R}$$v(x):\R^k\to\R$ is the linear combination of basis functions ${\phi }^p(x)\mapsto \mathbb{R}$$\varphi^p(x)\mapsto \R$

where $P$$P$ is the degree of the interpolation; $\hat{v}(x)$$\hat{v}(x)$ is called an interpolant; and $\overrightarrow{\theta }:=({\theta }^1,\dots ,{\theta }^P{)}^{\prime }\in {\mathbb{R}}^{P\times 1}$$\vec\theta := (\theta^1,\dots,\theta^P)' \in \R^{P\times 1}$ is the interpolation coefficient. In this post, we particularly discuss a class of linear interpolations in which ${\theta }^p$$\theta^p$ is the function value evaluated at supporting nodes from set ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$. Readers can quickly verify that most conclusions also hold for more generic types of interpolations in which $\overrightarrow{\theta }$$\vec\theta$ is not necessarily to be the value of $v(x)$$v(x)$ at the supporting nodes (e.g. orthogonal polynomials).

A stencil is a “template” of linear operations. Consider the following $N+1$$N+1$ degree interpolation:

where $C(z,({\hat{\mathcal{X}}}_N,{\mathcal{I}}_N))\in \mathbb{R}$$C(z,(\hat{\mathcal{X}}_N,\mathcal{I}_N)) \in \R$ is a constant which can be seen as $1\cdot C(\cdot )$$1\cdot C(\cdot)$; $\hat{x}[i]\in {\hat{\mathcal{X}}}_N$$\hat{x}[i]\in\hat{\mathcal{X}}_N$ is supporting node. We call the stacking basis vector $\overrightarrow{\phi }(z)$$\vec\varphi(z)$ a basis stencil which can fully determine the interpolated value at point $z\in \mathcal{X}$$z\in\mathcal{X}$ given the grid $({\hat{\mathcal{X}}}_N,\mathcal{I}{}_N)$$(\hat{\mathcal{X}}_N,\mathcal{I}{_N})$. The introduction of constant term $C$$C$ is to accommodate the strategy of hypothetical (ghost) outside nodes.

For interpolations where the coefficient is not the function value at supporting nodes, we can define a similar structure and also the basis stencil:

Remark: There are two common methods to handle boundary conditions:

1. Assume the boundary conditions hold exactly on the boundary nodes

2. Assume there exist hypothetical nodes outside but close to the boundary nodes. Assume boundary conditions hold on these hypothetical nodes.

The 2nd method is usually employed to generate solutions with smoother boundary behaviors.

The defined stencil representation allows for standard arithmetic such as addition, deduction, and scalars multiplication.

Monotonicity of numerical schemes

For second-order elliptic and parabolic PDEs, the equation's ellipticity determines both the system's tendency to reach a stationary state and its smoothness. In finite difference methods, a numerical scheme for a PDE consists of difference equations applied to discretized nodes and function values, approximating the original PDE. Within this framework, the monotonicity of a scheme is crucial for ensuring numerical stability and smoothness, and these properties are closely related.

In their seminal paper, Barles and Souganidis (1991) outline sufficient conditions for a numerical scheme to converge locally to the true solution, including monotonicity, consistency, and uniform stability. While consistency and uniform stability are typically easier to achieve in economic models, maintaining monotonicity requires special attention, as its absence is a common cause of failure in many proposed methods. In this post, we assume the other conditions are met and focus on the issue of monotonicity.

Consider a numerical scheme of for the boundary problem Eq ($\text{1}$$\ref{eq:01}$)-($\text{3}$$\ref{eq:03}$):

where $t$$t$ is the explicitly written time dimension which is used to update the guess across iterations; $D{\mathcal{V}}_N$$D\mathcal{V}_N$ and $D^2{\mathcal{V}}_N$$D^2 \mathcal{V}_N$ represent the stacked discretized Jacobian and Hessian, respectively which can be approximated by ${\mathcal{V}}_N$$\mathcal{V}_N$ with finite difference.

This equation of $S$$S$ consolidates all components of the discretized boundary problem equations onto one side. Monotonicity then requires that $S$$S$ be monotonically increasing with respect to $v[i+{\mathrm{\Delta }}^t]$$v[i+\Delta^t]$ but monotonically decreasing with respect to $v[i]$$v[i]$ and all other points in the stack ${\mathcal{V}}_N$$\mathcal{V}_N$, i.e., the neighbors of $v[i]$$v[i]$.

Example (Neo-classical growth model)

The neo-classical growth model:

is stationary (thus the time derivatives term $v_t=\partial v/\partial t=0$$v_t=\partial v/\partial t=0$). By plugging in the policy function $c=v_k^{-1/\gamma }$$c=v_k^{-1/\gamma}$, we get a non-linear elliptic PDE:

By adding back the always-zero time derivative term, we make it a non-linear parabolic PDE:

Suppose now we use forward difference along the $k$$k$ dimension to approximate $v_k$$v_k$ and use forward difference along the time $t$$t$ dimension to approximate $v_t$$v_t$ at an inner point $(k[i],v[i])$$(k[i],v[i])$ given discretization $({\hat{\mathcal{X}}}_N,{\mathcal{I}}_N,{\mathcal{V}}_N)$$(\hat{\mathcal{X}}_N,\mathcal{I}_N,\mathcal{V}_N)$, there is:

which is a difference equation about $v[i]$$v[i]$ and its neighbor points along $(t,k)$$(t,k)$ dimensions. Rearranging the RHS to the LHS, we then get the same format as $S(\cdot )=0$$S(\cdot)=0$. The example scheme in Eq ($\text{12}$$\ref{eq:11}$) is a non-linear scheme.

Remarks:

• The monotonicity is required also on the boundary condition Eq ($\text{3}$$\ref{eq:03}$), which is often ignored in practice.

• A commonly used alternative way of writing this (if applicable) is to split the terms about $v[i+{\mathrm{\Delta }}^t{e}^t]$$v[i+\Delta^t e^t]$ on one side of the equation, and put the left terms about ${\mathcal{V}}_N\setminus v[i+{\mathrm{\Delta }}^t{e}^t]$$\mathcal{V}_N\setminus v[i+\Delta^t e^t]$ to the other side. Then, both sides of the equation require monotonically increasing. We will use this style of organizing terms in the left part of this post.

• The example scheme in Eq ($\text{12}$$\ref{eq:11}$) is a one (time) layer scheme because by standing at a $t$$t$ moment, we only use $t$$t$ and $t+{\mathrm{\Delta }}^t$$t+\Delta^t$ information. Some more complex scheme may use more moments information. For example, Crank-Nicolson (two layer implicit method), "leapfrog" method (three layer explicit method), and Alternating direction implicit (ADI) method (multi layer method implicit method, fraction layers).

Explicit (forward Euler) method & CFL condition

Naive case

The explicit (forward Euler) method is one of the most intuitive schemes, relying solely on the previous iteration's guess to update the solution. It is called "forward" because it applies a forward difference along the time dimension. In the context of HJB equation, this approach typically avoids the need to solve any linear or nonlinear systems, requiring only a loop over each grid node. So it is easy to do parallelization. However, this method usually demands a small and adaptively adjusted time step to maintain stability, resulting in the so-called Courant–Friedrichs–Lewy (CFL) condition.

We illustrate the forward Euler method using the neo-classical growth model in Eq ($\text{11}$$\ref{eq:10}$). For convenience, we use $k[i]$$k[i]$ solely to index the spatial dimension $k$$k$, while a prime mark indicates the forward direction in time. This notation is permissible since the scheme involves only one time layer. For instance, $v^{\prime }[i]$$v'[i]$ refers to the right/forward neighbor at time $t+{\mathrm{\Delta }}^t$$t + \Delta^t$ of the current node $v[i]$$v[i]$ at time $t$$t$, with $v_t\approx (v^{\prime }[i]-v[i])/{\mathrm{\Delta }}^t$$v_t \approx (v'[i] - v[i])/\Delta^t$.

In Eq ($\text{12}$$\ref{eq:11}$), we provided an example of a fully non-linear scheme which is very hard to characterize its property. However, this is not the only approach. Due to the structure of the HJB equation, a common alternative is to treat the flow utility and drift terms as scalar functions without discretizing them. This renders the PDE linear with respect to the partial derivatives of the unknown function. Consistent with the directional choice in Eq ($\text{12}$$\ref{eq:11}$), we continue to use forward differences along the $k$$k$ dimension. Thus,

Rearranging the scheme:

The monotonicity requires:

while there is no special restriction on $u[i]$$u[i]$ if we do not discretize this term with $v[i]$$v[i]$ and other nodes. The first inequality always holds because ${\mathrm{\Delta }}^t>0$$\Delta^t>0$ for sure. The CFL condition comes from the 2nd inequality which implies $\mathrm{\forall }i\in {\mathcal{I}}_N$$\forall i\in \mathcal{I}_N$:

Eq ($\text{19}$$\ref{eq:18}$) is in the typical form of CFL condition that we could see from the textbook. It is easy to see that ${\mathrm{\Delta }}^t\to 0$$\Delta^t \to 0$ quickly as ${\mu }_k[i]$$\mu_k[i]$ increasing, which leads to slow convergence of the iteration.

The monotonicity of Eq ($\text{14}$$\ref{eq:13}$) scheme depends on the computation domain $\mathcal{X}$$\mathcal{X}$: larger the space, smaller time step required, then slower convergence of the algorithm, $\dots$$\dots$ so on until being infeasible. Meanwhile, The inequality of Eq ($\text{17}$$\ref{eq:16}$) is irrelevant to ${\mathrm{\Delta }}^t$$\Delta^t$ such that the scheme is monotonic only in the region of net saving. These bad properties show that Eq ($\text{14}$$\ref{eq:13}$) is not a working scheme. We need to find a strategy to improve Eq ($\text{14}$$\ref{eq:13}$).

Upwind scheme example: flux terms

In the context of second-order elliptic, parabolic, or hyperbolic PDEs, there is a strategy to improve a scheme’s monotonicity called the upwind scheme. This terminology reflects the physical intuition behind the approach, particularly for heat or wave equations. Here, the direction of differencing is set against the direction of wave or heat propagation, ensuring that information from the "origin of the wave/heat" (considered more crucial) is accurately incorporated. Mathematically, this strategy involves clamping the flux (drift) terms, which act as coefficients for the partial derivative terms, to ensure they remain positive or negative as required.

Let’s stick to the growth model example of Eq ($\text{11}$$\ref{eq:10}$), but this time we discretize the flux term $v_k\cdot {\mu }_k$$v_k \cdot \mu_k$ as:

For an arbitrary value of the drift ${\mu }_k[i]$$\mu_k[i]$, only one direction of differencing is acceptable, or else $v_k[i]\cdot {\mu }_k[i]=0$$v_k[i] \cdot \mu_k[i] = 0$, indicating a stationary state. But why Eq ($\text{20}$$\ref{eq:19}$) improves monotonicity? If we simplify Eq ($\text{20}$$\ref{eq:19}$), there is:

Noticing that all the terms of Eq ($\text{21}$$\ref{eq:20}$) appear on the RHS of the discretization. Specifically, let’s apply upwind scheme to the growth model example to get a full picture of this:

The monotonicity conditions for $v^{\prime }[i]$$v'[i]$, $v[i+{\mathrm{\Delta }}^k]$$v[i+\Delta^k]$, and $v[i-{\mathrm{\Delta }}^k]$$v[i-\Delta^k]$ are automatically satisfied. However, to ensure the monotonicity condition holds for the $v[i]$$v[i]$ term, we need to solve for the CFL condition with respect to ${\mathrm{\Delta }}^t$$\Delta^t$.

Compared to Eq ($\text{14}$$\ref{eq:13}$), applying the upwind scheme enables us to ensure monotonicity simply by adjusting the time step ${\mathrm{\Delta }}^t$$\Delta^t$. The ${\mathrm{\Delta }}^t$$\Delta^t$ now can handle both net saving and de-saving cases. Soon, we will demonstrate that the upwind scheme is still effective in the multi-dimensional case.

Remark: The illustrative examples by Benjamin Moll suggests an estimate of ${\mathrm{\Delta }}^t\approx 0.9\cdot \frac{{\mathrm{\Delta }}^k}{|{\mu }_k[i]|}$$\Delta^t \approx 0.9 \cdot \frac{\Delta^k}{|\mu_k[i]|}$, which works well for small $\rho$$\rho$. If you expect your discounting rate is significantly larger than the usual choice (e.g. $0.05$$0.05$), then a safer estimate is ${\mathrm{\Delta }}^t\approx 0.9\cdot \frac{{\mathrm{\Delta }}^k}{\rho {\mathrm{\Delta }}^k+|{\mu }_k[i]|}$$\Delta^t \approx 0.9 \cdot \frac{\Delta^k}{\rho \Delta^k + |\mu_k[i]|}$.

Upwind scheme example: diffusion terms

HJB equation is a 2nd-order PDE where there could be diffusion terms as Eq ($\text{1}$$\ref{eq:01}$). We have shown how to apply upwind scheme onto the flux terms. Now, let’s discuss how to apply upwind scheme to the diffusion terms. WLOG, let’s consider a stochastic neo-classical growth model:

where $v_{zz}\cdot \frac{1}{2}{\sigma }_z^2$$v_{zz} \cdot \frac{1}{2} \sigma^2_z$ represents the diffusion term in this equation. As before, we do not discretize ${\sigma }_z^2$$\sigma^2_z$ but instead discretize $v_{zz}$$v_{zz}$. Similarly, we use a three-point finite difference, given by $v_{zz}\approx (v[i+{\mathrm{\Delta }}^z{e}^z]-2v[i]+v[i-{\mathrm{\Delta }}^z{e}^z])/({\mathrm{\Delta }}^z{)}^2$$v_{zz} \approx ( v[i+\Delta^z e^z] - 2 v[i] + v[i-\Delta^z e^z] )/(\Delta^z)^2$. By clamping the (co)variance ${\sigma }_z^2$$\sigma^2_z$, we have:

Readers can instantly figure out the unconditional monotonicity of $v[i+{\mathrm{\Delta }}^z{e}^z]$$v[i+\Delta^z e^z]$ and $v[i-{\mathrm{\Delta }}^z{e}^z]$$v[i-\Delta^z e^z]$ terms. Then, collecting the $v[i]$$v[i]$ terms here and the $v[i]$$v[i]$ terms from the flux terms of state $k$$k$ and $z$$z$, we have the following CFL condition:

One may have realized that the CFL condition, in the case of explicit scheme + upwind scheme, has a formula for generic models as Eq ($\text{1}$$\ref{eq:01}$):

As we noted at the beginning of this section, the explicit method requires progressively smaller time steps as the dimensionality of the problem and the size of the computational space increase. For a medium-sized problem, this typically means millions of iterations are needed to achieve convergence.

Implicit (backward Euler) method

Compared to the explicit method, the implicit method, or backward Euler method, is often preferred in practice. For the HJB equation with an evenly spaced dense grid, this scheme is unconditionally monotonic for any size of ${\mathrm{\Delta }}^t$$\Delta^t$ when the upwind scheme is applied. However, the trade-off with the implicit method is that it requires solving a linear or nonlinear system at each iteration, which can be large but sparse in most cases.

The explicit method is “forward” in the meaning of forwardly difference along the time dimension (standing at moment $t$$t$). The implicit method is “backward” such that we backwardly difference along the time dimension, but we are no longer standing at moment $t$$t$ but standing at moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$ and looking back at moment $t$$t$.

A fully implicit scheme is one in which only $v_t$$v_t$ uses information from moment $t$$t$, while all other discretizations along spatial dimensions use information from moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$. This type of scheme typically results in nonlinear schemes (with respect to partial derivatives), which are challenging to analyze. Let’s use the deterministic growth model example again, but this time we have:

Observe that in the fully implicit scheme, both $u^{\prime }[i]$$u'[i]$ and ${\mu }_k^{\prime }[i]$$\mu'_k[i]$ are evaluated at moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$. These two terms are highly nonlinear in $D{v}^{\prime }$$Dv'$. We cannot take them as scalar functions as before because we cannot evaluate them without the knowledge about ${\mathcal{V}}^{\prime }$$\mathcal{V}'$. Their partial derivatives with respect to $v^{\prime }[i]$$v'[i]$ and its neighboring values being extremely complex. In most economic models, analyzing these derivatives is impractical.

Remark: We did not apply upwind scheme to Eq ($\text{32}$$\ref{eq:31}$). However, readers can quickly figure out that upwind scheme does not work here due to the high non-linearity of $u$$u$ and ${\mu }_k$$\mu_k$.

Thus, we focus on a class of quasi-implicit methods, where we evaluate the flow utility term, flux terms, and diffusion terms using the estimate at moment $t$$t$, while leaving the partial derivatives at moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$. This scheme corresponds precisely to the “implicit method” used by Benjamin Moll and other economists. Under such a quasi-implicit scheme and applying upwind scheme, the growth model is discretized as:

All coefficients of $v^{\prime }[i]$$v'[i]$, $v[i]$$v[i]$, and its spatial neighbors satisfy the monotonicity condition unconditionally, regardless of the size of the time step ${\mathrm{\Delta }}^t$$\Delta^t$. When collecting all $N+2k$$N+2k$ difference equations (including the $2k$$2k$ boundary conditions), one can obtain a linear system about the $N$$N$ nodes at moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$. Solving this linear system updates the guess of $\mathcal{V}$$\mathcal{V}$ for time-invariant HJB (Markovian), or update the guess of moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$ for finite-horizon HJB equations.

Remark: 1. In later sections of this post, we briefly introduce how to check the monotonicity when a transition matrix $L$$L$ is available. The check procedure suggests how to adaptively adjust the time step size ${\mathrm{\Delta }}^t$$\Delta^t$ per iteration. 2. A transition matrix is a linear approximation of the infinitesimal generator operator such that $L\cdot \mathbf{V}\approx \mathcal{L}\mathbf{V}$$L\cdot \mathbf{V} \approx \mathcal{L}\mathbf{V}$; it is basically a stacking of the flux & covariance terms in our previous examples.

Non-monotonicity of a class of approximation-based methods

The standard finite difference method discussed above is non-parametric, relying on a dense grid with no dependencies between any two nodes in ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$. However, this dense grid approach is highly susceptible to the curse of dimensionality, becoming quickly infeasible for three-dimensional and higher-dimensional problems. One may consider function interpolation or approximation methods (to approximate $v$$v$), which offer a wide range of choices. However, it is important to recognize that when these approximations of $v$$v$ are combined with finite difference methods, the monotonicity condition often fails to hold for a large class of these function approximation methods, even when the implicit method and upwind scheme are applied. This section is an intuitive illustration of the proof idea in Garcke and Ruttscheidt (2019).

We consider the two types of interpolation as outlined in the notation section.

First type

The first type uses interpolation coefficients $\overrightarrow{\theta }$$\vec\theta$ that represent the function values at supporting nodes from the set ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$. This type of interpolation can be abstracted as Eq ($\text{6}$$\ref{eq:06}$). Some examples: local methods such as piecewise polynomial, spline, etc. The multi-dimensional piecewise linear (aka multi-linear) interpolation exactly fits into the discussion of this sub-section.

The key idea here is that any approximated value $\hat{v}(z),z\in \mathcal{X}$$\hat{v}(z),z\in\mathcal{X}$ can be represented as a linear combination of all supporting $v[i]\in \mathcal{V}$$v[i]\in\mathcal{V}$. Taking this a step further, any linear operation on the approximated function values (e.g., finite difference) among arbitrary points $z_1,z_2,\cdots \in \mathcal{X}$$z_1, z_2, \dots \in \mathcal{X}$ can be represented by a basis stencil. For example:

where $\mathbf{V}$$\mathbf{V}$, as previously arranged, is the stack of function values at the supporting nodes, formatted as a vector. Notably, all finite difference approximations of partial derivatives fall within this category.

To illustrate the non-monotonicity issue when combining the above interpolation with the finite difference method, we use the stochastic growth model in Eq ($\text{27}$$\ref{eq:26}$) as an example. Here, we apply finite differences and the quasi-implicit scheme described in the previous section, but with a key modification: the finite difference is now computed between a supporting node and an interpolated neighboring point, rather than between two supporting nodes, for a given spatial mesh step size $({\mathrm{\Delta }}^k,{\mathrm{\Delta }}^z)$$(\Delta^k,\Delta^z)$. For instance, $v_k[i]\approx [\hat{v}(k[i]+{\mathrm{\Delta }}^k,z[i])-v[i]]/{\mathrm{\Delta }}^k$$v_k[i] \approx [\hat{v}(k[i] + \Delta^k, z[i]) - v[i]] / \Delta^k$.

Remark: The reason why not directly differencing two neighbor supporting nodes is that:

1. In some interpolations, there may not exist a “neighboring” node for supporting nodes in ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$

2. Even neighboring supporting nodes can be defined, there would be inconsistency issue which diverges the solution from the true solution. (Garcke & Ruttscheidt, 2019)

Remark: In this illustration, we set up a difference equation system over the $N$$N$ supporting nodes from ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$ such that the finite difference happens between a supporting node and an interpolated neighboring point. However, the conclusion of non-monotonicity trivially holds for a more generic case. This implies the failure of the following tries:

1. Sampling or designing an even-spaced “virtual” grid over $\mathcal{X}$$\mathcal{X}$

2. Evaluate the interpolants at these points

3. Apply standard finite difference + quasi implicit scheme

Discretizing Eq ($\text{27}$$\ref{eq:26}$) and apply quasi implicit scheme and upwind scheme:

which appears very similar to Eq ($\text{33}$$\ref{eq:33}$). One might mistakenly assume that Eq ($\text{36}$$\ref{eq:40}$) shares the same unconditional monotonicity as the conclusion from the previous section. However, it is crucial to recognize that monotonicity is required for all supporting nodes, not the interpolated points. Therefore, an additional step is needed: substituting the basis stencil and rearranging terms as in Eq ($\text{35}$$\ref{eq:35}$). A trick here is to represent the supporting node as an interpolated value but using unit vector $e^i\in {\mathbb{R}}^N$$e^i \in \R^N$ as basis stencil.

After simplification, Eq ($\text{37}$$\ref{eq:37}$) becomes a linear equation that simultaneously about all supporting nodes in ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$:

where ${\mathcal{A}}^{\prime }\in {\mathbb{R}}^{1\times N}$$\mathscr{A}'\in\R^{1\times N}$ is a complex linear combination of all basis stencils; $\mathcal{C}\in \mathbb{R}$$\mathscr{C}\in\R$ is the summation of all constant terms. The monotonicity condition, in fact, requires Eq ($\text{38}$$\ref{eq:38}$) to satisfy:

1. All elements in the basis stencil $\overrightarrow{\phi }[i]$$\vec\varphi[i]$ must be non-negative such that the coefficient of every supporting node’s function value $v[q],q\in {\mathcal{I}}_N$$v[q],q\in\mathcal{I}_N$ of moment $t$$t$ is non-negative.

2. All elements in the basis stencil ${\mathcal{A}}^{\prime }$$\mathscr{A}'$, except the current $i$$i$ node’s coefficient of moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$, must be non-positive.

3. Meanwhile, the $i$$i$ node’s coefficient of moment $t+{\mathrm{\Delta }}^t$$t+\Delta^t$ must be non-negative.

These three conditions actually impose around $(2_{\text{moments}}\cdot 2_{\text{dim}}\cdot 2_{\text{neighbor}})N$$(2_{\text{moments}} \cdot 2_{\text{dim}} \cdot 2_{\text{neighbor}})N$ sign restrictions on $\{\overrightarrow{\phi },{\overrightarrow{\phi }}^{\prime }\}$$\{\vec\varphi, \vec\varphi'\}$ and require these restrictions to hold in each iteration. If elements of $\overrightarrow{\phi }\in {\mathbb{R}}^N$$\vec\varphi \in \mathbb{R}^N$ can be negative, then satisfying these numerous restrictions becomes nearly impossible. This can be verified by manually summarizing the ${\mathcal{A}}^{\prime }$$\mathscr{A}'$ term and examining it. A more straightforward example illustrates this vulnerability: if any element of the basis stencil $\overrightarrow{\phi }[i]$$\vec\varphi[i]$ is negative, then monotonicity with respect to $\mathbf{V}$$\mathbf{V}$ (or $\mathcal{V}$$\mathcal{V}$) is violated, leaving no feasible solution for economists to restore it.

Even though we are discussing a specific example, it contains all the ingredients needed to trivially generalize it to an arbitrary $k$$k$ dimension HJB equation driven by an $m$$m$-dimensional Brownian motion, as shown in Eq ($\text{1}$$\ref{eq:01}$). I leave the formal statement for those who are interested in further exploration.

The issue is evident, but what could be a potential solution? The answer lies in using interpolation methods that represent any point $z\in \mathcal{X}$$z \in \mathcal{X}$ as a convex combination of the function values of supporting nodes in ${\hat{\mathcal{X}}}_N$$\hat{\mathcal{X}}_N$. Examining Eq ($\text{38}$$\ref{eq:38}$), if all elements of the basis stencils used in this equation are non-negative, then applying the upwind scheme to Eq ($\text{38}$$\ref{eq:38}$) preserves unconditional monotonicity. There are specific types of interpolation that meet this convexity requirement.

Remark: The non-negativity of basis stencil or function elements is NOT equivalent to the monotonicity of the interpolant. Some readers may recognize methods like monotonic interpolations (e.g., monotonic cubic splines). These interpolations, while designed to ensure monotonicity of the interpolant, do not necessarily guarantee monotonicity of numerical schemes—two forms of "monotonicity" that are fundamentally different.

More specifically, monotonic interpolation ensures that within a single iteration, if the guessed solution is monotonic along certain dimensions, then evaluating the interpolant will preserve this property. In contrast, the monotonicity of a numerical scheme ensures that, as the guessed solution is updated from one iteration to the next, the monotonicity of the previous iteration’s guess is maintained in the next iteration, thereby preserving the solution’s shape over successive updates. This is sometimes referred as “shape-preserving” property of the scheme.

Second type

The second type of interpolation we consider does not use function values at supporting nodes as interpolation coefficients but instead utilizes a generic ${\overrightarrow{\theta }}_P$$\vec\theta_P$. Here, the degree of interpolation $P$$P$ is not necessarily equal to the number of supporting nodes $N$$N$. Examples of this type of interpolation include global methods such as orthogonal polynomials (spectral), neural networks, and Gaussian processes, among others.

In short, this type of interpolation also suffers from the same issue as the first type. The key to understand this is to recognize that, after substituting the interpolation into Eq ($\text{38}$$\ref{eq:38}$), we are simply projecting the original space $\mathcal{V}$$\mathcal{V}$ onto the space of interpolation coefficients ${\overrightarrow{\theta }}_P$$\vec\theta_P$, while everything else remains unchanged. For example (note that now even $v[i]$$v[i]$ is the evaluated value),

which results in a collection of $P+2k$$P+2k$ difference equations. If the stencil elements cannot be guaranteed to be non-negative, then the same issue of irreparable monotonicity arises once again.

The solution in this case is the same as for the first type of interpolation: employ interpolation methods that ensure non-negative stencils for all $z\in \mathcal{X}$$z \in \mathcal{X}$. For example, one might consider using ReLU, Sigmoid, or other non-negative activation functions in neural network implementations.

Bonus section: determine which scheme to use

The combined strategy of the quasi-implicit method and the upwind scheme is sufficient for handling a broad class of economic models in research. However, certain models with unique structures may require customization beyond the standard approach. Learning how to design tailored linear or nonlinear schemes for specific problems can be highly beneficial (even though most customized schemes works less good than the standard approach when the latter is available).

A snack example, the absolute value of a partial derivative $|v_k|$$|v_k|$ appears in some applications (e.g. adjustment costs). The trick to monotonically discretize it is:

i.e. use backward difference when $v_k\ge 0$$v_k\geq 0$, and use forward difference when $v_k<0$$v_k<0$.

The general steps to determine what numerical scheme to use is:

1. Define the economic model as a boundary problem and solve policy functions in the interior and on the boundary.

2. Plugging the policy functions back and observe:

   1. How tractable the problem is? If maximum principle can gives the analytical solution, then solve it. If not, then consider numerically solve it.

   2. Is finite difference the best choice? If yes, move on. If not, try generalized residual methods (I’ll discuss this in another post)

   3. How non-linear the problem is? If linear enough or easy to linearize, then consider linearize the PDE and solve it analytically or numerically. e.g. some portfolio or risk-neutral problems can be approximated as a Linear Quadratic Regulator (LQR) problem which is tractable

   4. Dimensionality of the problem? If the dimensionality is low, then primarily use even-spaced grid + standard strategy. If the dimensionality is high, then choose proper approximation methods that satisfy the above non-negativity conditions, and potentially combine it with other techniques such as sparse grids.

3. If finite difference is your final choice, then

   1. Design your grid: even-spaced dense grid? uneven-spaced dense grid (e.g. Chebyshev)? sparse grid (e.g. Smolyak)? adaptive sparse grid (e.g. multi-linear ASG)?

   2. Does your model has mixed frequency of noises that require multi-grid methods? $\to$$\to$ If yes, then things go more complex and you may need an advanced course of numerical PDE.

   3. Based on your grid, decide if to use approximation methods

      1. Even-spaced dense grid: No need to use approximation, everything is standard

      2. Uneven-spaced dense grid: No need to use approximation, but needs to be careful about the mesh step size for difference

      3. (Regular or adaptive) sparse grid: Must use approximation. Carefully choose approximation methods that satisfy the above conditions

   4. Design the discretization strategy, i.e. the scheme. This is the heart part.

      1. How many layers do you need? $\to$$\to$ This is needed in some special cases such as doing dimension reduction using methods like ADI, or you are handling some special time-dependent (parabolic) problems.

      2. Is fully implicit scheme possible to be unconditionally monotonic? $\to$$\to$ If yes, then use it. This is faster convergent than the quasi implicit scheme. If not, one may need to compare it (after applying CFL condition) with a quasi implicit scheme.

      3. Any non-linearity that must be discretized? $\to$$\to$ If no, then great and just fit your PDE into the standard strategy. If yes, you need more tricks to handle these non-linearities before moving forward.

      4. Do you get an excellent idea which beats the standard strategy even though the latter is applicable? $\to$$\to$ If yes, then go ahead. If no, stay with the standard strategy.

   5. Carefully double check if all Barles-Souganidis sufficient conditions are satisfied.

   6. Carefully verify the chosen scheme on the discretized boundary conditions. The PDE, esp. elliptic/parabolic/hyperbolic PDEs are pinned down greatly by boundary conditions. For multi-dimensional problems, pay extra attention to “corners” i.e. where multiple boundary conditions need to simultaneously hold.

   7. Write your code, debug, test and run.

   Reference

   1. Garcke, J., & Ruttscheidt, S. (2019). Finite differences on sparse grids for continuous time heterogeneous agent models.

   2. Barles, G., & Souganidis, P. E. (1991). Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic analysis, 4(3), 271-283.