Finite Difference Methods for Pricing Derivatives

Worked Solution

How to Think About It: A pricing PDE with terminal condition is solved numerically by discretizing time and the state variable, replacing derivatives by finite differences, and stepping from the payoff at expiry back to today. The only design decision is *when* to evaluate the space derivative: at the known time level (explicit), at the unknown one (implicit), or halfway (Crank-Nicolson). That choice trades computational simplicity against stability.

Approach: Transform Black-Scholes to the heat equation, set up the grid, derive the three schemes from Taylor expansions, and state their properties. A short Python implementation of the explicit scheme is included and checked against the Black-Scholes formula.

Formal Solution:

Part (a): From Black-Scholes to the heat equation

*Step 1 -- The PDE.* A European derivative $V(S, t)$ on a non-dividend stock satisfies $$\frac{\partial V}{\partial t} + \tfrac12\sigma^2 S^2\frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0, \qquad V(S, T) = \text{payoff}(S).$$

*Step 2 -- Change of variables.* Put $S = Ke^{x}$, $t = T - 2\tau/\sigma^2$, and $V = K\,v(x, \tau)$. Then $v_\tau = v_{xx} + (k - 1)v_x - kv$ with $k = 2r/\sigma^2$, a constant-coefficient equation. Removing the first-order and zero-order terms with $v = e^{\alpha x + \beta\tau}u(x, \tau)$, $\alpha = -\tfrac12(k - 1)$, $\beta = -\tfrac14(k + 1)^2$, leaves the heat equation $$\frac{\partial u}{\partial \tau} = \frac{\partial^2 u}{\partial x^2},$$ with initial condition (at $\tau = 0$, i.e. expiry) $u(x, 0) = e^{-\alpha x}\,\text{payoff}(Ke^x)/K$; for a call, $u(x, 0) = \max\!\left(e^{(k+1)x/2} - e^{(k-1)x/2}, 0\right)$. The terminal-value problem has become an initial-value problem marching forward in $\tau$.

Part (b): Grid and difference approximations

*Step 3 -- Grid.* Truncate $x$ to $[x_{\min}, x_{\max}]$ (far enough into and out of the money) and set $x_j = x_{\min} + j\,\delta x$, $j = 0, \dots, J$; $\tau_n = n\,\delta\tau$, $n = 0, \dots, N$, with $N\delta\tau = \sigma^2 T/2$. Write $u_j^n \approx u(x_j, \tau_n)$. Boundary values $u_0^n$ and $u_J^n$ come from the asymptotic behaviour of the option (a call is worthless as $x \to -\infty$ and behaves like the forward minus strike as $x \to \infty$).

*Step 4 -- Differences.* Taylor expansion gives the central second difference $$\frac{\partial^2 u}{\partial x^2}(x_j, \tau_n) = \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{\delta x^2} + O(\delta x^2),$$ and for the time derivative either the forward difference $(u_j^{n+1} - u_j^n)/\delta\tau + O(\delta\tau)$ or the backward one at level $n + 1$.

Part (c): The three schemes

*Step 5 -- Explicit scheme.* Evaluate $u_{xx}$ at the known level $n$: $$\frac{u_j^{n+1} - u_j^n}{\delta\tau} = \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{\delta x^2} \quad\Longrightarrow\quad u_j^{n+1} = a\,u_{j-1}^n + (1 - 2a)\,u_j^n + a\,u_{j+1}^n, \qquad a = \frac{\delta\tau}{\delta x^2}.$$ Each new value is a weighted average of three old neighbours, $O(J)$ work per step, no linear algebra. Truncation error $O(\delta\tau) + O(\delta x^2)$. It is stable only if $a \le \tfrac12$ (all weights nonnegative; von Neumann analysis gives amplification factor $1 - 4a\sin^2(\theta/2)$), which forces $\delta\tau \le \delta x^2/2$, i.e. many time steps when the space grid is fine.

*Step 6 -- Implicit scheme.* Evaluate $u_{xx}$ at the new level $n + 1$: $$-a\,u_{j-1}^{n+1} + (1 + 2a)\,u_j^{n+1} - a\,u_{j+1}^{n+1} = u_j^n.$$ Each step is a tridiagonal linear system in the $J - 1$ unknowns, solved in $O(J)$ by the Thomas algorithm. Unconditionally stable for every $a > 0$; truncation error still $O(\delta\tau) + O(\delta x^2)$.

*Step 7 -- Crank-Nicolson.* Average the explicit and implicit right-hand sides (centred in time at $\tau_{n + 1/2}$): $$-\tfrac{a}{2}u_{j-1}^{n+1} + (1 + a)u_j^{n+1} - \tfrac{a}{2}u_{j+1}^{n+1} = \tfrac{a}{2}u_{j-1}^{n} + (1 - a)u_j^{n} + \tfrac{a}{2}u_{j+1}^{n}.$$ Also tridiagonal and unconditionally stable, and now $O(\delta\tau^2) + O(\delta x^2)$ accurate. It is the workhorse for vanilla and barrier options; for non-smooth payoffs it can show oscillations near the kink, often damped by a few fully implicit start-up steps (Rannacher smoothing).

*Step 8 -- Back to prices.* After stepping to $\tau_N$, $V(S, 0) = K e^{\alpha x + \beta\tau_N}u(x, \tau_N)$ at $x = \ln(S/K)$; interpolate between grid points. American options are handled by applying $u \ge$ intrinsic value after every step (explicit) or via a projected SOR / penalty method (implicit).

```python import numpy as np def bs_call_explicit_fd(S0, K, r, sigma, T, J=600, xmin=-3.0, xmax=3.0, a=0.45): k = 2 * r / sigma**2 alpha, beta = -(k - 1) / 2, -(k + 1)**2 / 4 dx = (xmax - xmin) / J tau_max = sigma**2 * T / 2 N = int(np.ceil(tau_max / (a * dx2))); dt = tau_max / N; a = dt / dx2 x = np.linspace(xmin, xmax, J + 1) u = np.maximum(np.exp((k + 1) / 2 * x) - np.exp((k - 1) / 2 * x), 0.0) # payoff at tau = 0 tau = 0.0 for _ in range(N): tau += dt u[1:-1] = a * u[:-2] + (1 - 2 * a) * u[1:-1] + a * u[2:] u[0] = 0.0 u[-1] = np.exp((k + 1) / 2 * xmax + (k + 1)**2 / 4 * tau) - np.exp((k - 1) / 2 * xmax + (k - 1)**2 / 4 * tau) V = K * np.exp(alpha * x + beta * tau) * u return np.interp(np.log(S0 / K), x, V) # bs_call_explicit_fd(100, 100, 0.05, 0.2, 1.0) -> 10.4499 ; Black-Scholes: 10.4506 ```

Answer:

(a) With $S = Ke^x$, $\tau = \sigma^2(T - t)/2$ and $V = Ke^{\alpha x + \beta\tau}u$ ($\alpha = -(k-1)/2$, $\beta = -(k+1)^2/4$, $k = 2r/\sigma^2$), the Black-Scholes PDE becomes the heat equation $u_\tau = u_{xx}$ with the payoff as initial condition.

(b) Discretize $x$ and $\tau$ on a grid, approximate $u_{xx}$ by the central second difference $(u_{j+1} - 2u_j + u_{j-1})/\delta x^2$ and $u_\tau$ by a one-step difference, and march from the payoff toward today with boundary values from the option's asymptotics.

(c) Explicit: $u_j^{n+1} = a u_{j-1}^n + (1 - 2a)u_j^n + a u_{j+1}^n$, $a = \delta\tau/\delta x^2$, trivial to compute but stable only for $a \le 1/2$. Implicit: tridiagonal system $-a u_{j-1}^{n+1} + (1 + 2a)u_j^{n+1} - a u_{j+1}^{n+1} = u_j^n$, unconditionally stable, first order in time. Crank-Nicolson: average of the two, unconditionally stable and second order in time.