Finite difference PDE methods

Definition Assume the Black-Scholes setting with constant $r$, $q$ and $\sigma > 0$, and let $V(t, S)$ be the value of a European claim with payoff $h(S_T)$, solving

for $t < T$ and $S > 0$, with $V(T, S) = h(S)$. In the remaining time $\tau = T - t$ and the log price $x = \ln S$, the chain rule gives $S V_S = V_x$ and $S^2 V_{SS} = V_{xx} - V_x$, so

Writing $V(\tau, x) = e^{-r\tau} u(\tau, y)$ with $y = x + b\tau$ gives $V_\tau = e^{-r\tau}(u_\tau + b u_y - r u)$, $V_x = e^{-r\tau} u_y$ and $V_{xx} = e^{-r\tau} u_{yy}$. The discount and drift terms cancel and leave the heat equation

which runs forward in $\tau$, starting from the payoff. With $a > 0$ it is well posed forward in time and smooths its initial data, while running it backward would be the ill-posed backward heat problem. When $h(e^y)$ grows at most exponentially in $y$, the solution given by Feynman-Kac is the Gaussian average $u(\tau, y) = \mathbb{E}\bigl[h(e^{y + \sigma\sqrt{\tau}\, Z})\bigr]$ with $Z$ standard normal, since $2a\tau = \sigma^2 \tau$. With coefficients that depend on the state this reduction is not available, and a scheme has to discretize the full operator.

Definition The finite difference method replaces derivatives by difference quotients on a grid and solves the resulting algebraic equations. A uniform spatial grid cuts the line down to an interval $[y_{\min}, y_{\max}]$ with the nodes $y_j = y_{\min} + j\,\Delta y$ for $j = 0, \dots, J$, together with the time levels $\tau_n = n\,\Delta\tau$, and $u_j^n$ denotes the approximation of $u(\tau_n, y_j)$. For a function with four continuous derivatives, Taylor expansion gives

for some $\xi$ between $y - h$ and $y + h$, and a one-sided difference in $\tau$ has an error of order $\Delta\tau$. The values at $y_0$ and $y_J$ come from boundary conditions that approximate the contract far from the strike. The computed answer therefore carries a domain truncation error, which shrinks as the interval widens, besides the errors in $\Delta y$ and $\Delta\tau$, which shrink as the grid is refined. A scheme is consistent when its difference equations approach the PDE as $\Delta\tau, \Delta y \to 0$, and stable when, over a fixed horizon, it amplifies perturbations of its data by at most a fixed factor that does not grow as the grid is refined. For a well-posed linear initial value problem, the Lax equivalence theorem, stated here without proof, says that a consistent scheme converges exactly when it is stable.

Definition The explicit scheme, forward Euler in $\tau$, evaluates the second difference at the old level,

With the ratio $\lambda = a\,\Delta\tau/(\Delta y)^2$ this becomes

for $1 \le j \le J - 1$, with boundary values at $j = 0$ and $j = J$. Each new value comes directly from three old ones, with no linear system, and the local truncation error is of order $\Delta\tau + (\Delta y)^2$ for smooth solutions. The stability condition of this scheme, in this heat scaling, is

equivalently $\Delta\tau \le (\Delta y)^2/\sigma^2$. Von Neumann analysis on an infinite grid inserts the Fourier mode $u_j^n = g^n e^{\mathrm{i} j\theta}$, with $\mathrm{i}$ the imaginary unit, and finds the amplification factor

which satisfies $\lvert g(\theta) \rvert \le 1$ for every $\theta$ exactly when $\lambda \le \tfrac12$. Above one half the alternating mode $\theta = \pi$, the pattern $1, -1, 1, \dots$, is multiplied by $1 - 4\lambda < -1$ at every step and grows without bound. On a fixed bounded grid with fixed boundary values the largest stable ratio is slightly above one half, and it tends to one half as $\Delta y \to 0$, so $\lambda \le \tfrac12$ is the condition that works for every grid.

Definition A scheme satisfies the discrete maximum principle when every new interior value lies between the smallest and the largest of the values it is computed from, which are the old interior values and the boundary values used in the step. For the explicit scheme with $0 \le \lambda \le \tfrac12$, the weights $\lambda$, $1 - 2\lambda$ and $\lambda$ are nonnegative and add to $1$, so each new value is a convex combination of old values and

where $u_0^n$ and $u_J^n$ are the boundary values. As consequences, nonnegative data stay nonnegative, data between two bounds stay between them when the boundary values do too, and with zero boundary values the largest absolute value never increases, which proves stability in the maximum norm. For $\lambda > \tfrac12$ the center weight is negative and the principle fails. The implicit scheme satisfies the principle for every $\lambda > 0$, and for Crank-Nicolson $\lambda \le 1$ is a sufficient condition. The property is especially useful for comparison arguments and for American obstacle problems, whose algorithms lie beyond this lesson.

Definition The implicit scheme, backward Euler in $\tau$, evaluates the second difference at the new level. With the discrete operator $(L_h u)_j = a\,(u_{j-1} - 2u_j + u_{j+1})/(\Delta y)^2$ at the interior nodes it reads

where the boundary values $\lambda u_0^{n+1}$ and $\lambda u_J^{n+1}$ move to the right-hand side of the first and last rows. Each interior row is

The matrix is tridiagonal with diagonal $1 + 2\lambda$ and off-diagonal entries $-\lambda$, so it is strictly diagonally dominant with nonpositive off-diagonal entries. It is therefore an M-matrix, whose inverse has only nonnegative entries, and the system is solved in a number of operations proportional to $J$ by the tridiagonal (Thomas) algorithm. The amplification factor

lies in $(0, 1]$ for every $\lambda > 0$, so the scheme is stable with no restriction on $\Delta\tau$, and it satisfies the discrete maximum principle for every $\lambda$. Its error is of order $\Delta\tau + (\Delta y)^2$ for smooth solutions, only first order in time, so a large stable step can still be inaccurate. Unconditional stability is not unconditional accuracy.

Definition The Crank-Nicolson scheme averages the explicit and implicit spatial differences,

with the boundary values moved to the right-hand side. Each interior row is

Its error is of order $(\Delta\tau)^2 + (\Delta y)^2$ for sufficiently smooth solutions. The amplification factor

satisfies $\lvert g(\theta) \rvert \le 1$ for every $\lambda$, so the scheme is stable in the discrete $\ell^2$ norm with no restriction on $\Delta\tau$. Stability does not give the maximum principle. The left matrix is an M-matrix for every $\lambda$, and the right-hand weights $\tfrac12\lambda$, $1 - \lambda$ and $\tfrac12\lambda$ are nonnegative when $\lambda \le 1$, so $\lambda \le 1$ is a sufficient condition for the discrete maximum principle. For large $\lambda$ the alternating mode has $g(\pi) = (1 - 2\lambda)/(1 + 2\lambda)$, close to $-1$, so it flips sign at every step and decays slowly, while the exact heat flow removes it almost at once. Near the kink of a call payoff this appears as oscillations in the value and especially in delta and gamma. A common remedy, Rannacher smoothing, takes a few implicit steps, often with half the step size, before switching to Crank-Nicolson.