Explicit Finite Differences for a Call: Stability, Cost, Convergence
You price a European call with $S_0 = K = 100$, $r = 0.05$, $\sigma = 0.2$, $T = 1$ (no dividends) by explicit finite differences. With $k = 2r/\sigma^2$, the change of variables $S = Ke^x$, $\tau = \sigma^2(T - t)/2$, $V = K e^{\alpha x + \beta\tau}u(x, \tau)$, $\alpha = -(k - 1)/2$, $\beta = -(k + 1)^2/4$ turns the Black-Scholes PDE into the heat equation $u_\tau = u_{xx}$ for $\tau \in [0, \sigma^2 T/2]$, with the payoff as initial condition. You truncate to $x \in [-3, 3]$ with $J$ equal space steps, use $N$ equal time steps, and march $$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}.$$
(a) What Black-Scholes price should the grid reproduce?
(b) With $J = 600$, what is the smallest number of time steps $N$ for which the scheme is stable?
(c) Keeping $a$ fixed, you refine to $J = 1200$. By about what factor does the total number of grid-point updates grow?
(d) At $a \approx 0.45$ the grid gives $P_{600} = 10.44992$ and $P_{1200} = 10.45042$. What do these numbers say about the error, and what should a run with $J = 2400$ give?
Hints
- At the money, $d_1 = (r + \sigma^2/2)\sqrt{T}/\sigma$ and $d_2 = d_1 - \sigma\sqrt{T}$; the call is $S_0 N(d_1) - K e^{-rT} N(d_2)$.
- The weight $1 - 2a$ must not go negative: stability needs $a \le 1/2$. Here $\delta x = 6/J$ and the total heat-equation time is $\tau_{\max} = \sigma^2 T/2$, so $\delta\tau = \tau_{\max}/N$.
- At fixed $a$, $\delta\tau = a\,\delta x^2$, so the truncation error $O(\delta\tau) + O(\delta x^2)$ is $O(\delta x^2)$ overall. Compare $P_{600}$ and $P_{1200}$ with the price from (a).
Worked Solution
How to Think About It: A finite difference price is only useful if you know how far it is from the truth and what it costs to get closer. Three facts settle that for the explicit scheme: the stability bound $a \le 1/2$ fixes the minimum number of time steps, refining at fixed $a$ costs $2 \times 4 = 8$ times the work, and the error at fixed $a$ is $O(\delta x^2)$. The Black-Scholes formula gives the exact answer to measure against.
Quick Estimate: At the money the call is roughly $0.4\,\sigma S\sqrt{T}$ plus about half the rate carry, $0.4 \times 20 + 0.5 \times 100(1 - e^{-0.05}) \approx 8 + 2.4 \approx 10.4$. With $\delta x = 0.01$, $a \le 1/2$ means $\delta\tau \le 5 \times 10^{-5}$ over a total time of $0.02$: about $400$ steps.
Formal Solution:
*Step 1 -- The benchmark (a).* $d_1 = (0.05 + 0.02)/0.2 = 0.35$, $d_2 = 0.15$, $N(d_1) = 0.63683$, $N(d_2) = 0.55962$: $$C = 100 \times 0.63683 - 100e^{-0.05} \times 0.55962 = 63.683 - 53.232 = 10.4506.$$
*Step 2 -- Stability (b).* The explicit update multiplies a Fourier mode by $1 - 4a\sin^2(\theta/2)$, which stays in $[-1, 1]$ for all $\theta$ only if $a \le 1/2$ (equivalently, all three weights are nonnegative). Here $\delta x = 6/600 = 0.01$ and $\tau_{\max} = \sigma^2 T/2 = 0.02$, so $$a = \frac{0.02/N}{0.01^2} = \frac{200}{N} \le \frac12 \iff N \ge 400.$$ At $N = 400$ the run is stable ($10.4454$, less accurate than at $a = 0.45$); at $N = 300$ ($a = 0.67$) it explodes to about $-2 \times 10^{64}$.
*Step 3 -- Cost of refinement (c).* Doubling $J$ halves $\delta x$. Keeping $a = \delta\tau/\delta x^2$ fixed then needs $\delta\tau$ cut by $4$, so $N$ grows by $4$ (from $445$ to $1778$ at $a = 0.45$). The number of node updates is about $J \times N$, so it grows by $2 \times 4 = 8$ (the exact ratio of $(J - 1)N$ is $7.998$).
*Step 4 -- Error and extrapolation (d).* Against $10.45058$, the errors are $P_{600} - C = -0.00066$ and $P_{1200} - C = -0.00016$, a ratio of about $4$. That is what $O(\delta x^2)$ predicts: at fixed $a$, $\delta\tau = a\,\delta x^2$, so the first-order time error is also $O(\delta x^2)$. The next halving should cut the error by another $4$: $$P_{2400} \approx 10.45058 - \frac{0.00016}{4} \approx 10.45054,$$ and the run gives exactly $10.45054$. Writing $P_J \approx C + c\,\delta x^2$ and eliminating $c$ gives Richardson extrapolation, $$\frac{4P_{1200} - P_{600}}{3} = \frac{41.80168 - 10.44992}{3} = 10.45059,$$ within $10^{-5}$ of Black-Scholes at no extra cost. Truncating $x$ to $[-3, 3]$ is not the issue: widening to $[-4, 4]$ at the same $\delta x$ changes the price by less than $10^{-12}$.
The scheme in code:
```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, J=600) -> 10.44992 (N = 445) # bs_call_explicit_fd(100, 100, 0.05, 0.2, 1.0, J=1200) -> 10.45042 (N = 1778) # bs_call_explicit_fd(100, 100, 0.05, 0.2, 1.0, J=2400) -> 10.45054 (N = 7112); Black-Scholes: 10.45058 ```
Answer: (a) $10.4506$. (b) $N = 400$ (need $a = 200/N \le 1/2$). (c) About $8$ times the work ($2\times$ the nodes, $4\times$ the steps). (d) The error falls by about $4$ per halving of $\delta x$ ($-0.00066$ to $-0.00016$), as $O(\delta x^2)$ predicts at fixed $a$, so $J = 2400$ should give about $10.45054$; Richardson extrapolation $(4P_{1200} - P_{600})/3 = 10.45059$ already matches Black-Scholes.
Intuition
The explicit scheme is a three-point weighted average, cheap per step but only stable while every weight is nonnegative, $a = \delta\tau/\delta x^2 \le 1/2$. That ties the time step to the square of the space step: here $\delta x = 0.01$ forces at least $400$ time steps, and each halving of $\delta x$ doubles the nodes and quadruples the steps, so the work grows eightfold. The reward is that the error is $O(\delta x^2)$ at fixed $a$: it falls by $4$ per refinement, which you can see in the two runs and exploit with Richardson extrapolation, $(4P_{1200} - P_{600})/3$, to land on the Black-Scholes price to five decimals. Implicit and Crank-Nicolson schemes remove the stability constraint, which is why they are the default for production grids.