Case study: pricing a European option on a grid

From the equation to a grid

Case study 6.3.4 stated the Black–Scholes equation for the value $V(S, t)$ of an option on a stock with price $S$ at time $t$:

At expiry $T$ a call is worth its payoff, $V(S, T) = \max(S - K, 0)$: a terminal condition (lesson 6.3.2). The inputs are case study 2.4.6's: strike $K = 100$, $T = 0.5$ years, rate $r = 0.03$ and volatility $\sigma = 0.2$, both yearly decimals, with the stock at 100 today.

Case study 6.4.4 counted time backwards from expiry, $\tau = T - t$, so that $V_t = -V_\tau$. Then the payoff becomes a starting profile:

This runs forward in $\tau$ from the payoff, as lesson 6.5.1's rod ran forward from its first profile. A grid can work with $S$ directly, without 6.4.4's change of variables: it only needs the coefficients at each grid point.

Take grid prices $S_i = i\,\Delta S$ for $i = 0, \ldots, M$, up to $S_{\max} = M\,\Delta S$, and $N$ time steps of $\Delta\tau = \frac{T}{N}$. Replace $V_{SS}$ by the second difference and $V_S$ by the central difference (lesson 2.8.1). Because $S_i = i\,\Delta S$, the factors of $\Delta S$ cancel: $\frac{S_i^2}{\Delta S^2} = i^2$ and $\frac{S_i}{2\Delta S} = \frac{i}{2}$.

The weights change from row to row, but the pattern is lesson 6.5.2's. At the strike, $i = 100$, the second-difference weight times $\Delta\tau$ is $\frac{1}{2}\sigma^2 i^2\,\Delta\tau$, which plays the part of lesson 6.5.1's $r$.

Boundary conditions

Each step needs the two end values, as in lesson 6.5.1. Case study 6.3.4 found them, written here with $\tau$:

• At $S = 0$. A stock worth 0 stays at 0, so the call never pays: $V(0, \tau) = 0$. The equation agrees: at $S = 0$ it reads $V_\tau = -rV$.
• At $S_{\max}$. Far above the strike, the call is almost sure to be exercised, so it is worth the stock less the discounted strike: $V(S_{\max}, \tau) = S_{\max} - Ke^{-r\tau}$. By put–call parity (case study 3.3.4), the error is the put's value there.
• A put has the mirror conditions: $V(0, \tau) = Ke^{-r\tau}$ and $V(S_{\max}, \tau) = 0$.

Case study 6.2.4 found that a wrong value at a far boundary drags the whole curve, and fixed it by putting the true value there. The grid makes the same choice: its far value is the one the call nearly has, not 0 or the payoff.

$S_{\max}$ must be far enough above the strike for the second condition to hold. Over six months the log of the stock price spreads by about $\sigma\sqrt{T} = 0.1414$ (case study 5.6.5). Four times that above the strike is $100e^{4 \times 0.1414} \approx 176$.

The last row is the true call price at $S_{\max}$ minus the boundary value. At 120 it is 57 cents off, and the price at 100 is 2.4 cents low. From 150 up, $S_{\max}$ no longer changes the price at 100 to 6 decimals.

Pricing the call

Take $S_{\max} = 200$ with $M = 200$, so $\Delta S = 1$, and $N = 100$ steps of $\Delta\tau = 0.005$. The strike and the spot, both 100, are grid prices.

The pricer starts with Rannacher's start (lesson 6.5.2): two implicit steps of half the length, then Crank–Nicolson. The next section shows why.

The grid gives 6.367513 at $S = 100$, against the formula's 6.371028: 0.0035 low, within half a cent. The remaining error comes from the grid spacing, and the refinement section below shows it shrinking.

The grid holds a value only at its nodes. For a spot between them, step from the nearest node $S_j$ with case study 2.7.4's delta–gamma estimate, using the grid's own delta and gamma (next section):

At 103.3 this gives 8.396155, against the formula's 8.399483: off by 0.0033, the grid's own error.

Delta and gamma from the grid

Case study 2.8.4 found delta and gamma by bumping the stock price and repricing. On a grid the neighbours are already there: the values at 99 and 101 are a $1 bump each way, at no extra cost.

At $S = 100$ the grid gives delta 0.570080 and gamma 0.027791, against the formula's 0.570158 and 0.027772. The figure shows gamma at every node beside the formula.

Too few time steps

Lesson 6.5.2 warned that Crank–Nicolson with large steps leaves a zig-zag at a kink, and a payoff has one at the strike. Take only $N = 10$ steps, so the ratio $\frac{1}{2}\sigma^2 i^2\,\Delta\tau$ at the strike is 10:

Plain Crank–Nicolson's gamma swings between $-0.0728$ and 0.1846, exactly at the strike where the spot sits. With 5 steps its price is 6.5925, 22 cents high. Rannacher's start brings the 5-step price to 6.3580, within 1.3 cents.

Grid refinement and its cost

Halve $\Delta S$ and $\Delta\tau$ together, keeping $S_{\max} = 200$ and $\Delta\tau = \frac{\Delta S}{200}$:

Each halving divides the error by 4, as Crank–Nicolson's $\Delta\tau^2 + \Delta S^2$ promises, and multiplies the work by 4. The error is within a cent from $\Delta S = 1$ on. Delta and gamma settle too, to 0.570153 and 0.027773 on the finest grid.

Because the error falls by 4 each time, lesson 2.8.2's Richardson extrapolation applies: $\frac{4 \times 6.367513 - 6.356929}{3} = 6.371041$, only $1.3 \times 10^{-5}$ from the formula.

Set that against the other methods. Case study 5.11.4's tree with 1,000 steps has 501,501 nodes and an error of 0.0014; the grid with 80,000 values does better, at 0.00088.

Case study 5.7.5's Monte Carlo price, 6.3765 with a standard error of 0.0095 from a million draws, is far rougher still.

The explicit scheme of lesson 6.5.1 is no shortcut here. Its ratio $\frac{1}{2}\sigma^2 i^2\,\Delta\tau$ grows like $i^2$, so the top of the grid sets the limit, about $\Delta\tau \le \frac{1}{\sigma^2M^2}$: at least $\sigma^2M^2T = 800$ steps for $M = 200$. With 400 steps it blows up; with 800 it gives 6.3684.

An American put

Case study 5.11.4 priced an American put on a tree by taking, at each node, the larger of holding and exercising. A grid does the same at each time step: after each Crank–Nicolson solve, replace every value by the larger of it and the payoff $K - S_i$.

The European column reaches the Black–Scholes put, 4.8822, as the call did. The American column creeps up on the tree's 5.0098, more slowly. Its distance from the tree falls from 0.0059 to 0.0020, 0.0007 and 0.0003: about 2 to 3 times per halving, not 4.

The reason is a new corner. Where the put switches from holding to exercising, its value meets the payoff, and that point moves as time passes. On the 800-price grid the highest price where the put is exercised today is 81.25, matching the tree's critical price of about 81.

import math

from scipy.stats import norm


def bs(S, K, T, r, sigma, kind="call"):
    d1 = (math.log(S / K) + (r + sigma**2 / 2) * T) / (sigma * math.sqrt(T))
    call = S * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d1 - sigma * math.sqrt(T))
    return call if kind == "call" else call - S + K * math.exp(-r * T)


def expect_value_error(call, why):
    try:
        call()
    except ValueError:
        return
    raise AssertionError(f"expected ValueError: {why}")


def test_shared_call():
    """The 6-month call prices within 1 cent of 6.3710, with delta within 0.001 of 0.570158 and gamma within 0.0001 of 0.027772"""
    price, delta, gamma = grid_price(100, 100, 0.5, 0.03, 0.2)
    expect_close(price, bs(100, 100, 0.5, 0.03, 0.2), what="price", abs_tol=0.01, rel=0)
    expect_close(delta, 0.570158, what="delta", abs_tol=1e-3, rel=0)
    expect_close(gamma, 0.027772, what="gamma", abs_tol=1e-4, rel=0)


def test_case_study_grid():
    """With the default grid (S up to 200 in steps of 1, 100 time steps) the price is the case study's 6.367513"""
    expect_close(grid_price(100, 100, 0.5, 0.03, 0.2)[0], 6.367513, what="price on the default grid", abs_tol=1e-4, rel=0)


def test_few_time_steps():
    """With only 5 time steps, Rannacher's start keeps the price within 2 cents and gamma within 0.005"""
    price, delta, gamma = grid_price(100, 100, 0.5, 0.03, 0.2, N=5)
    expect_close(price, bs(100, 100, 0.5, 0.03, 0.2), what="price with N = 5", abs_tol=0.02, rel=0)
    expect_close(gamma, 0.027772, what="gamma with N = 5", abs_tol=0.005, rel=0)


def test_put_call_parity():
    """On the same grid, call minus put equals S0 − K e^(−rT) to within 1e-5"""
    call = grid_price(100, 100, 0.5, 0.03, 0.2, "call")[0]
    put = grid_price(100, 100, 0.5, 0.03, 0.2, "put")[0]
    expect_close(call - put, 100 - 100 * math.exp(-0.015), what="call − put", abs_tol=1e-5, rel=0)


def test_spot_between_nodes():
    """A spot of 103.3, between nodes, prices within 1 cent of the formula"""
    expect_close(grid_price(103.3, 100, 0.5, 0.03, 0.2)[0], bs(103.3, 100, 0.5, 0.03, 0.2), what="price at S0 = 103.3", abs_tol=0.01, rel=0)


def test_other_options():
    """A put with strike 90, a call deep in the money and a one-year put at 40% volatility each price within 1 cent"""
    expect_close(grid_price(100, 90, 0.5, 0.03, 0.2, "put")[0], bs(100, 90, 0.5, 0.03, 0.2, "put"), what="put, K = 90", abs_tol=0.01, rel=0)
    expect_close(grid_price(150, 100, 0.5, 0.03, 0.2)[0], bs(150, 100, 0.5, 0.03, 0.2), what="call, S0 = 150", abs_tol=0.01, rel=0)
    expect_close(grid_price(100, 110, 1, 0.05, 0.4, "put", S_max=400, M=400)[0], bs(100, 110, 1, 0.05, 0.4, "put"),
                 what="put, K = 110, T = 1, σ = 0.4", abs_tol=0.01, rel=0)


def test_bad_inputs():
    """Non-positive inputs, an unknown kind, too small a grid, a fractional step count or a spot off the grid raise ValueError"""
    expect_value_error(lambda: grid_price(0, 100, 0.5, 0.03, 0.2), "S0 = 0")
    expect_value_error(lambda: grid_price(100, 100, 0, 0.03, 0.2), "T = 0")
    expect_value_error(lambda: grid_price(100, 100, 0.5, 0.03, -0.2), "a negative volatility")
    expect_value_error(lambda: grid_price(100, 100, 0.5, 0.03, 0.2, "digital"), "'digital' is not call or put")
    expect_value_error(lambda: grid_price(100, 100, 0.5, 0.03, 0.2, M=2), "M = 2 leaves one inside point")
    expect_value_error(lambda: grid_price(100, 100, 0.5, 0.03, 0.2, N=10.5), "10.5 time steps")
    expect_value_error(lambda: grid_price(250, 100, 0.5, 0.03, 0.2), "S0 = 250 lies beyond S_max = 200")