Case study: the normal integrals behind Black–Scholes

The future stock price

The Black–Scholes model is derived in the Stochastic Calculus course and used throughout the Derivatives Modelling course. Here we take one claim of it as given: the stock price at expiry $T$ is

Here $z$ is a standard normal number: any real value, weighted by the density $\varphi(z)$ of lesson 2.6.3. Each value of $z$ is one possible future, and $\varphi(z)$ says how much weight it gets.

The inputs follow the desk's conventions. $S_0 = 100$ is today's price and $K = 100$ the strike, both in dollars. $T = 0.5$ is in years, and $r = 0.03$ and $\sigma = 0.2$ are annual decimals, with $r$ continuously compounded (case study 1.4.6).

Where the option starts to pay

A call pays $\max(S_T - K, 0)$ at expiry: the ramp of lesson 1.1.2. It pays nothing unless $S_T > K$. Setting $S_T = K$ and solving for $z$ with logs (lesson 1.4.4) gives the corner $z^* = -d_2$, where

These $d_1$ and $d_2$ are the ones in the formula of case study 2.4.6. For our inputs $z^* \approx -0.035355$, so the call pays off in every future with $z$ above that: a little more than half of them.

The price as an integral

The model's price is the payoff averaged over all futures, each weighted by $\varphi(z)$, then discounted by $e^{-rT}$ to today.

Where $N(d_2)$ and $N(d_1)$ come from

Above $z^*$ the payoff is $S_T - K$, so the price splits into two integrals over $z > -d_2$: a stock part and a strike part.

The strike part

$\int_{-d_2}^{\infty}\varphi(z)\,dz$ is the area under the bell above $-d_2$. The bell is symmetric (lesson 2.6.3), so that equals the area below $d_2$, which is $N(d_2)$.

So the strike part is $K e^{-rT}N(d_2) \approx 50.644782$. Here $N(d_2) \approx 0.514102$ is the model's chance that the call ends in the money, above the strike.

The stock part, by substitution

1. 1
   $$\begin{aligned} S_T(z)\,\varphi(z) &= \frac{S_0\,e^{(r - \sigma^2/2)T}}{\sqrt{2\pi}} \\ &\quad \times e^{\sigma\sqrt{T}z - z^2/2} \end{aligned}$$

   Why?
2. 2
   $$\begin{aligned} &\sigma\sqrt{T}\,z - \frac{z^2}{2} \\ &= -\frac{(z - \sigma\sqrt{T})^2}{2} + \frac{\sigma^2 T}{2} \end{aligned}$$

   Why?
3. 3
   $$\begin{aligned} &S_T(z)\,\varphi(z) \\ &= S_0\,e^{rT}\,\varphi\big(z - \sigma\sqrt{T}\big) \end{aligned}$$

   Why?
4. 4
   $$\begin{aligned} &e^{-rT}\int_{-d_2}^{\infty} S_T\,\varphi\,dz \\ &= S_0\int_{-d_1}^{\infty}\varphi(w)\,dw \end{aligned}$$

   Why?
5. 5
   $$\begin{aligned} &e^{-rT}\int_{-d_2}^{\infty} S_T\,\varphi\,dz \\ &= S_0\,N(d_1) \end{aligned}$$

   Why?

The stock part is $S_0 N(d_1) \approx 57.01581$. Together, $c = S_0 N(d_1) - K e^{-rT} N(d_2) \approx 57.01581 - 50.644782 = 6.371028$: the formula of case study 2.4.6.

Step 3 over the whole line gives one more check. The shifted bell has area 1, so the model's average future price is $S_0 e^{rT} \approx 101.5113$: the stock grows on average at the rate $r$.

On the desk

$N$ has no formula, so it is computed

$N(x)$ is an area under $\varphi$, and lesson 2.6.3 showed it has no formula in familiar functions. Every pricing library computes it numerically. SciPy's norm.cdf uses tuned approximations accurate to about 16 significant digits, close to the limit of a float.

So "the Black–Scholes formula" is exact algebra wrapped around two numerical integrals. Pricing by quad directly, as above, is the same calculation done more slowly; it is useful for payoffs with no formula.

Cutting the infinite integral off

Any numerical method must stop somewhere. Cutting the call integral at $z = L$ instead of $\infty$ leaves out the tail beyond $L$.

Step 3 shows where the weight sits. The stock part's integrand is a bell centred at $\sigma\sqrt{T}$, not at 0, so the tail starts later when $\sigma\sqrt{T}$ is large.

For a 5-year option with $\sigma = 0.8$, $\sigma\sqrt{T} \approx 1.79$ and the price is $65.6312. Cutting at 4 then leaves out $1.3, about 2%. Measure the cut from $\sigma\sqrt{T}$ instead: 6 beyond it leaves out less than $0.0000001 for both options.

Tails and precision

A float keeps about 16 significant digits. So $1 - N(x)$, a tiny number found by subtracting from 1, loses digits as $x$ grows.

At $x = 9$, $N(9)$ is $1 - 1.1 \times 10^{-19}$, which rounds to exactly 1, and 1 - norm.cdf(9) gives 0. By symmetry the same tail is $N(-9) \approx 1.1 \times 10^{-19}$, computed with no subtraction; SciPy also offers norm.sf(x) for $1 - N(x)$. Worse, 1 - quad(norm.pdf, -math.inf, 9) returns a small negative number: a negative probability.

Desks meet this in far out-of-the-money options. The put integral, done the same way, gives $p = K e^{-rT}N(-d_2) - S_0 N(-d_1)$.

For a strike of 40, writing $N(-d)$ as $1 - N(d)$ gets the fifth digit of the price wrong. Such tiny prices feed ratios, such as implied volatilities, where every digit counts.

Bad inputs and units

• At expiry, $T = 0$, the formulas for $d_1$ and $d_2$ divide by zero. There is nothing left to average: the price is the payoff, $\max(S_0 - K, 0)$.
• With $\sigma \to 0$ there is no randomness, and the price tends to $\max(S_0 - Ke^{-rT}, 0)$, about $1.488806 here. $\sigma = 0$ itself divides by zero, so code must treat it separately.
• Volatility as a decimal. Typing $\sigma = 20$ for 20% prices the call at $100.0, the whole stock.
• Time in years. Typing $T = 6$ for six months gives $27.1509 instead of $6.3710.
• Positive prices. $\ln(S_0/K)$ needs $S_0 > 0$ and $K > 0$; a zero or negative input raises an error.

import math

from scipy.stats import norm


def _bs_call(S0, K, T, r, sigma):
    d1 = (math.log(S0 / K) + (r + sigma**2 / 2) * T) / (sigma * math.sqrt(T))
    return S0 * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d1 - sigma * math.sqrt(T))


def test_call_matches_formula():
    """The 6-month at-the-money call is worth about 6.3710, as the formula says"""
    expect_close(call_by_integral(100, 100, 0.5, 0.03, 0.2), _bs_call(100, 100, 0.5, 0.03, 0.2), what="call_by_integral", rel=1e-9)


def test_out_of_the_money_call():
    """A strike-120 call matches the formula, about 0.8914"""
    expect_close(call_by_integral(100, 120, 0.5, 0.03, 0.2), _bs_call(100, 120, 0.5, 0.03, 0.2), what="call, K = 120", rel=1e-8)


def test_put_call_parity():
    """Call minus put equals S0 − K e^(−rT)"""
    c = call_by_integral(80, 75, 1.0, 0.04, 0.3)
    p = put_by_integral(80, 75, 1.0, 0.04, 0.3)
    expect_close(c - p, 80 - 75 * math.exp(-0.04), what="c − p", rel=1e-8)


def test_put_positive():
    """A put is worth something: the at-the-money put is about 4.8822"""
    expect_close(put_by_integral(100, 100, 0.5, 0.03, 0.2), 4.882221902473731, what="put_by_integral", rel=1e-8)


def test_at_expiry():
    """At expiry (T = 0) each option is worth its payoff: a call 10 when S0 = 110, a put 10 when S0 = 90"""
    expect_close(call_by_integral(110, 100, 0.0, 0.03, 0.2), 10.0, what="call_by_integral at T = 0, S0 = 110")
    expect_close(call_by_integral(90, 100, 0.0, 0.03, 0.2), 0.0, what="call_by_integral at T = 0, S0 = 90")
    expect_close(put_by_integral(90, 100, 0.0, 0.03, 0.2), 10.0, what="put_by_integral at T = 0, S0 = 90")