Case study: from Black–Scholes to the heat equation

The equation, and four obstacles

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

Here $r$ is the continuous interest rate and $\sigma$ the volatility, both annual decimals. A call with strike $K$ and expiry $T$ has the terminal condition $V(S, T) = \max(S - K, 0)$. The inputs are those of case study 2.4.6: $K = 100$, $T = 0.5$, $r = 0.03$, $\sigma = 0.2$, and today $t = 0$ with $S = 100$.

Next to $u_s = u_{xx}$, the equation has four differences, and each change below removes one.

1. Its condition is at the end, $t = T$, not at the start.
2. Its coefficients $S^2$ and $S$ change with the stock price.
3. It has terms in $V_S$ and $V$ that the heat equation lacks.
4. Its diffusion term has a factor $\frac{1}{2}\sigma^2$, not 1.

Step 1: time to expiry

Case study 6.3.4 found that the equation runs backwards in time, so it needs a condition at expiry. Measure time backwards from expiry instead: $\tau = T - t$, the time to expiry in years. Then $V_t = -V_\tau$ (lesson 6.3.2), and

This runs forwards in $\tau$, as heat does (lesson 6.4.1). The payoff becomes a starting profile at $\tau = 0$, and today is $\tau = 0.5$.

Step 2: the log of the price

Measure the price on a log scale relative to the strike: $x = \ln(S/K)$, so $S = Ke^x$. A price of 100 is $x = 0$, the strike itself; 110 is $x = 0.09531$ and 90 is $x = -0.105361$. Lesson 6.2.2 used the same logarithm on the equidimensional equation, whose coefficients also grow with powers of the variable.

Write $V(S, t) = W(x, \tau)$. The chain rule (lessons 2.3.3 and 4.2.2) gives the partials in $S$, with $\frac{dx}{dS} = \frac{1}{S}$:

1. 1
   $$V_S = \frac{W_x}{S}$$

   Why?
2. 2
   $$V_{SS} = \frac{W_{xx} - W_x}{S^2}$$

   Why?
3. 3
   $$\begin{gathered} S\,V_S = W_x \\[4pt] S^2\,V_{SS} = W_{xx} - W_x \end{gathered}$$

   Why?

Substituting these into the equation of step 1 gives an equation with constant coefficients:

Prices from 0 to infinity become all of $x$, from $-\infty$ to $\infty$: the endless rod of lesson 6.4.2. The payoff becomes $W(x, 0) = K\max(e^x - 1, 0)$.

Step 3: a factor that removes two terms

The terms in $W_x$ and $W$ remain. Write $W = K e^{ax + b\tau}\,u(x, \tau)$ for two numbers $a$ and $b$ still to choose. The product rule gives each partial as $Ke^{ax + b\tau}$ times a bracket:

Put these into the equation of step 2 and divide by $Ke^{ax + b\tau}$. The $u_x$ terms add up to $\left(\sigma^2 a + r - \frac{1}{2}\sigma^2\right)u_x$, and the right $a$ makes that 0. The right $b$ then cancels every term in $u$, leaving only $u_\tau = \frac{1}{2}\sigma^2 u_{xx}$.

Step 4: a clock for heat

$u_\tau = \frac{1}{2}\sigma^2 u_{xx}$ is a heat equation with diffusivity $D = \frac{1}{2}\sigma^2 = 0.02$ per year. As in lesson 6.4.1, change the clock to $s = \frac{1}{2}\sigma^2\tau$; then $u_s = u_{xx}$. Today, with half a year to go, $s = 0.02 \times 0.5 = 0.01$.

$s$ has no units, because $\sigma^2$ is per year and $\tau$ is in years. The whole half-year is only 0.01 of heat-equation time: a short spell of smoothing.

The payoff becomes a kink

At $\tau = 0$, $u = W e^{-ax}/K$, so the starting profile is

using $1 - a = 1.25$ and $-a = 0.25$. It is 0 for $x \le 0$, and its slope just right of 0 is $1.25 - 0.25 = 1$. So near the strike it is lesson 6.4.2's kink $\max(x, 0)$, bending up a little more further out.

The slope is 1 for every input, since $(1 - a) - (-a) = 1$. So every call becomes a kink-like profile; only how much it bends up depends on $r$ and $\sigma$.

Solving with the heat kernel

Now lesson 6.4.2 applies. With $s = 0.01$ today,

and the price is $V = Ke^{ax + b\tau}u(x, s)$. At $S = 100$, $x = 0$, so $V = 100\,e^{-0.015625}\,u(0, 0.01)$.

The bell $\Phi(\cdot, 0.01)$ has width $\sqrt{2s} = \sqrt{\sigma^2\tau} = 0.141421$. That is $\sigma\sqrt{T}$, the standard deviation of the six-month log return in case study 5.6.5. In lesson 6.4.2's reading, $u(x, s)$ is $u_0$ averaged over a random shift of $x$: here, the stock's random log return to expiry.

quad from the corner at 0 to ten widths above $x$ gives $u(0, 0.01) \approx 0.0647136$, and the price is 6.371028. The kink alone would give $\sqrt{s/\pi} \approx 0.0564$: the extra upward bend of $u_0$ makes $u$ 1.147 times as large.

The two rows agree to within $10^{-12}$ at every price. Three months later, at $t = 0.25$, $s$ is 0.005 and the call is worth 4.3576 at $S = 100$ and 11.5715 at $S = 110$, again matching the formula.

Back to the Black–Scholes formula

The integral can also be done by hand. Above 0, $u_0$ is two exponentials, $e^{c_1y} - e^{c_2y}$ with $c_1 = 1 - a = 1.25$ and $c_2 = -a = 0.25$. Lesson 6.4.2 gave each one cut at 0:

Multiply by $Ke^{ax + b\tau}$ and add up the exponents. Since $2s = \sigma^2\tau$, they simplify for every $r$ and $\sigma$:

• First piece: $a + c_1 = 1$ and $b\tau + c_1^2 s = 0$, so it is $Ke^x N(\cdot) = S\,N(d_1)$.
• Second piece: $a + c_2 = 0$ and $b\tau + c_2^2 s = -r\tau$, so it is $Ke^{-r\tau}N(d_2)$.
• Inside $N$: $2c_1 s = \left(r + \frac{1}{2}\sigma^2\right)\tau$ and $2c_2 s = \left(r - \frac{1}{2}\sigma^2\right)\tau$.

The last line makes the arguments $d_1 = \frac{\ln(S/K) + (r + \sigma^2/2)\tau}{\sigma\sqrt{\tau}}$ and $d_2 = d_1 - \sigma\sqrt{\tau}$. So the heat solution is

the formula of case study 2.4.6. Today $d_1 = 0.176777$ and $d_2 = 0.035355$, and the two pieces are 57.01581 and 50.644782: the stock part and the strike part of case study 2.6.5.

The heat route and case study 2.6.5's route are the same calculation. That case study averaged the payoff over the normal log return. Here the heat kernel does the averaging, and the factor $e^{ax + b\tau}$ supplies the discounting and the drift.

On the desk

Black and Scholes found their formula this way, by reducing their equation to the heat equation. Using the change of variables in code brings its own traps.

The put, from the same equation

The put pays $\max(K - S, 0)$. Steps 1 to 4 are unchanged; only the starting profile changes, to $u_0 = \max\left(e^{0.25x} - e^{1.25x}, 0\right)$, which is positive for $x < 0$. The heat integral gives 4.882222, case study 5.11.4's European put. Put–call parity holds: $6.371028 - 4.882222 = 1.488806 = S - Ke^{-rT}$.

A factor that can overflow

The factor is mild here, but $a = \frac{1}{2} - \frac{r}{\sigma^2}$ explodes as $\sigma$ shrinks. At $\sigma = 0.5\%$, $a = -1199.5$, so $u_0(y)$ contains $e^{1200.5\,y}$. That overflows a float for $y$ above about 0.59, a stock price 1.81 times the strike.

The cure is to multiply before exponentiating. Inside the integral, $Ke^{ax + b\tau}u_0(y)$ is $Ke^{b\tau}e^{a(x - y)}\max(e^y - 1, 0)$, and on the bell's window $x - y$ is tiny. Written that way, a call with $S = 200$ at 0.5% volatility prices at 101.488806 instead of crashing.

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))
    d2 = d1 - sigma * math.sqrt(T)
    if kind == "call":
        return S * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d2)
    return K * math.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)


def test_shared_call():
    """The 6-month at-the-money call is worth 6.3710"""
    expect_close(price_by_heat(100, 100, 0.5, 0.03, 0.2), _bs(100, 100, 0.5, 0.03, 0.2), what="price_by_heat(100, 100, 0.5, 0.03, 0.2)", rel=1e-8)


def test_shared_put():
    """The matching put is worth 4.8822"""
    expect_close(price_by_heat(100, 100, 0.5, 0.03, 0.2, "put"), _bs(100, 100, 0.5, 0.03, 0.2, "put"), what="price_by_heat(..., 'put')", rel=1e-8)


def test_out_of_the_money_call():
    """A strike-120 call is worth 0.8914"""
    expect_close(price_by_heat(100, 120, 0.5, 0.03, 0.2), _bs(100, 120, 0.5, 0.03, 0.2), what="price_by_heat(100, 120, 0.5, 0.03, 0.2)", rel=1e-8)


def test_other_inputs():
    """A 1-year put with S = 80, K = 100, r = 5% and 35% volatility is worth 21.1650"""
    expect_close(price_by_heat(80, 100, 1, 0.05, 0.35, "put"), _bs(80, 100, 1, 0.05, 0.35, "put"), what="price_by_heat(80, 100, 1, 0.05, 0.35, 'put')", rel=1e-8)


def test_low_volatility():
    """At 0.5% volatility a call with S = 200, K = 100 is worth 101.4888, with no overflow"""
    try:
        value = price_by_heat(200, 100, 0.5, 0.03, 0.005)
    except OverflowError:
        raise AssertionError("OverflowError: with a = 1/2 - r/sigma^2 = -1199.5, combine e^(a(x - y)) inside the integral")
    expect_close(value, _bs(200, 100, 0.5, 0.03, 0.005), what="price_by_heat(200, 100, 0.5, 0.03, 0.005)", rel=1e-8)


def test_at_expiry():
    """At T = 0 the price is the payoff: 5 for a call at 105 and for a put at 95"""
    expect_close(price_by_heat(105, 100, 0, 0.03, 0.2), 5.0, what="call at expiry")
    expect_close(price_by_heat(95, 100, 0, 0.03, 0.2, "put"), 5.0, what="put at expiry")


def test_bad_inputs():
    """A zero volatility or an unknown kind raises ValueError"""
    for args, kind in [((100, 100, 0.5, 0.03, 0.0), "call"), ((100, 100, 0.5, 0.03, 0.2), "straddle")]:
        try:
            price_by_heat(*args, kind=kind)
        except ValueError:
            continue
        raise AssertionError(f"expected ValueError for sigma={args[4]}, kind={kind!r}")