Black and Scholes's d2¶

Black and Scholes's d2 = d1 - s, with s = sigma sqrt(T / P) the volatility over the option's life (lusine.finance.d1, whose parameters it takes): Hull's d1 - sigma sqrt(T).

N(d2) is the risk-neutral chance the option finishes in the money. Source: Black & Scholes, The Pricing of Options and Corporate Liabilities, J. Political Economy 81 (1973); Hull, Options, Futures, and Other Derivatives, 9th ed. (Pearson, 2014), ch. 15. Checked against a reference in numpy, which has its own tests.

Source

Black & Scholes, J. Political Economy 81 (1973); Hull, Options, Futures, and Other Derivatives, 9th ed. (2014), ch. 15.

Proof

reference, invariant.

d2.trait

namespace lusine.finance;

float d2(float<_> spot, float<_> strike, float<year^-1> rate, float volatility, float<year> period, float<year> years)
{
  return d1(spot, strike, rate, volatility, period, years) - volatility * sqrt(years / period);
}

Uses, without proving it here: finance.d1.

D1 - sigma sqrt(T).

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ndtr
from scipy.stats import norm
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32
YEAR = 31557600.0                        # lusine's year row, the Julian year in seconds: the host speaks base units

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from scipy.special import ndtr from scipy.stats import norm from lusine import trait_run EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32 YEAR = 31557600.0 # lusine's year row, the Julian year in seconds: the host speaks base units

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('finance/d2.trait'))

print(trait_run.source('finance/d2.trait'))

// lusine.finance.d2: Black and Scholes's d2 = d1 - s, with s = sigma sqrt(T / P) the volatility over the option's life
// (lusine.finance.d1, whose parameters it takes): Hull's d1 - sigma sqrt(T). N(d2) is the risk-neutral chance the
// option finishes in the money.
// Source: Black & Scholes, The Pricing of Options and Corporate Liabilities, J. Political Economy 81 (1973); Hull,
// Options, Futures, and Other Derivatives, 9th ed. (Pearson, 2014), ch. 15.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float d2(float<_> spot, float<_> strike, float<year^-1> rate, float volatility, float<year> period, float<year> years)
{
  return d1(spot, strike, rate, volatility, period, years) - volatility * sqrt(years / period);
}

The program¶

A lattice of options, one an item, the kernel bound over them: the prices in euros, the rate a year, the volatility over a period (a year here, Hull's), the life in years. The host speaks the base units, so a time is stored in seconds and a rate per second, and d2 is dimensionless. kernel below is that call; closed is the independent account, Black-Scholes in float64 with scipy's exact normal distribution, at the parameters as the device holds them. held rounds them so.

In [3]:

CALL = False
BS = """namespace bs;

lusine.version version = { { number = 1 } };

struct option {
  host float<EUR> spot; host float<EUR> strike; host float<year^-1> rate; host float volatility;
  host float<year> period; host float<year> years;
  float d2;
};
lattice<option, 1> options(spacing = 1 m, boundary = wall);

void main()
{
  options.d2 = lusine.finance.d2(options.spot, options.strike, options.rate, options.volatility,
                                    options.period, options.years);
}
"""
print(BS)


def held(call, spot, strike, rate, volatility, years, period=1.0):
    """The parameters as the device holds them, in the host's units (euros, a year, years): each rounded to binary32 in
    the base units it is stored in."""
    f = lambda v, scale=1.0: np.float32(np.asarray(v, dtype=float) * scale).astype(float) / scale
    return (np.asarray(call, dtype=bool), f(spot), f(strike), f(rate, 1 / YEAR), f(volatility), f(years, YEAR),
            f(period, YEAR))


def kernel(call, spot, strike, rate, volatility, years, period=1.0):
    """The kernel's d2 for each option: prices in euros, the rate a year, the volatility over `period` years and the
    life in years, stored in the base units (seconds); the result back in dimensionless, one an option (a
    single option's alone)."""
    n = max(np.size(a) for a in (call, spot, strike, rate, volatility, years, period))
    run = trait_run.Program(['finance/d2.trait'], BS, {'bs.options': [n]})
    for field, v in (('spot', spot), ('strike', strike), ('rate', np.asarray(rate) / YEAR),
                     ('volatility', volatility), ('period', np.asarray(period) * YEAR),
                     ('years', np.asarray(years) * YEAR)):
        run.store('bs.options', field, np.broadcast_to(np.asarray(v), (n,)))
    out = run.steps('bs.main', 1, [('bs.options', 'd2')])[0][0]
    return out if n > 1 else float(out[0])


def closed(call, spot, strike, rate, volatility, years, period=1.0):
    """The independent account: Black-Scholes in float64 with scipy's exact normal distribution, at the parameters
    the device holds."""
    c, S, K, r, v, T, P = held(call, spot, strike, rate, volatility, years, period)
    s = v * np.sqrt(T / P)
    a = (np.log(S / K) + (r + 0.5 * (v * v / P)) * T) / s
    b = a - s
    df = K * np.exp(-r * T)
    return b


rng = np.random.default_rng(391)
N = 1000
SPOT, STRIKE = rng.uniform(20, 200, N), rng.uniform(20, 200, N)
RATE, VOL, LIFE = rng.uniform(-0.02, 0.15, N), rng.uniform(0.05, 0.9, N), rng.uniform(0.02, 5.0, N)
CALLS = rng.random(N) < 0.5

CALL = False BS = """namespace bs; lusine.version version = { { number = 1 } }; struct option { host float<EUR> spot; host float<EUR> strike; host float<year^-1> rate; host float volatility; host float<year> period; host float<year> years; float d2; }; lattice<option, 1> options(spacing = 1 m, boundary = wall); void main() { options.d2 = lusine.finance.d2(options.spot, options.strike, options.rate, options.volatility, options.period, options.years); } """ print(BS) def held(call, spot, strike, rate, volatility, years, period=1.0): """The parameters as the device holds them, in the host's units (euros, a year, years): each rounded to binary32 in the base units it is stored in.""" f = lambda v, scale=1.0: np.float32(np.asarray(v, dtype=float) * scale).astype(float) / scale return (np.asarray(call, dtype=bool), f(spot), f(strike), f(rate, 1 / YEAR), f(volatility), f(years, YEAR), f(period, YEAR)) def kernel(call, spot, strike, rate, volatility, years, period=1.0): """The kernel's d2 for each option: prices in euros, the rate a year, the volatility over `period` years and the life in years, stored in the base units (seconds); the result back in dimensionless, one an option (a single option's alone).""" n = max(np.size(a) for a in (call, spot, strike, rate, volatility, years, period)) run = trait_run.Program(['finance/d2.trait'], BS, {'bs.options': [n]}) for field, v in (('spot', spot), ('strike', strike), ('rate', np.asarray(rate) / YEAR), ('volatility', volatility), ('period', np.asarray(period) * YEAR), ('years', np.asarray(years) * YEAR)): run.store('bs.options', field, np.broadcast_to(np.asarray(v), (n,))) out = run.steps('bs.main', 1, [('bs.options', 'd2')])[0][0] return out if n > 1 else float(out[0]) def closed(call, spot, strike, rate, volatility, years, period=1.0): """The independent account: Black-Scholes in float64 with scipy's exact normal distribution, at the parameters the device holds.""" c, S, K, r, v, T, P = held(call, spot, strike, rate, volatility, years, period) s = v * np.sqrt(T / P) a = (np.log(S / K) + (r + 0.5 * (v * v / P)) * T) / s b = a - s df = K * np.exp(-r * T) return b rng = np.random.default_rng(391) N = 1000 SPOT, STRIKE = rng.uniform(20, 200, N), rng.uniform(20, 200, N) RATE, VOL, LIFE = rng.uniform(-0.02, 0.15, N), rng.uniform(0.05, 0.9, N), rng.uniform(0.02, 5.0, N) CALLS = rng.random(N) < 0.5

namespace bs;

lusine.version version = { { number = 1 } };

struct option {
  host float<EUR> spot; host float<EUR> strike; host float<year^-1> rate; host float volatility;
  host float<year> period; host float<year> years;
  float d2;
};
lattice<option, 1> options(spacing = 1 m, boundary = wall);

void main()
{
  options.d2 = lusine.finance.d2(options.spot, options.strike, options.rate, options.volatility,
                                    options.period, options.years);
}

Reference: scipy's Black-Scholes¶

A thousand options: each d2 within 1e-6 (1 + |d2|) of the float64 account.

In [4]:

b = kernel(None, SPOT, STRIKE, RATE, VOL, LIFE)
want = closed(True, SPOT, STRIKE, RATE, VOL, LIFE)
err = np.abs(b - want) / (1 + np.abs(want))
assert np.max(err) <= 1e-6, np.max(err)
print(f'largest error {np.max(err):.2e} of 1 + |d2|')

b = kernel(None, SPOT, STRIKE, RATE, VOL, LIFE) want = closed(True, SPOT, STRIKE, RATE, VOL, LIFE) err = np.abs(b - want) / (1 + np.abs(want)) assert np.max(err) <= 1e-6, np.max(err) print(f'largest error {np.max(err):.2e} of 1 + |d2|')

largest error 3.44e-07 of 1 + |d2|

Invariant: the chance of finishing in the money¶

N(d2) is the risk-neutral chance a call finishes in the money: a stock drawn at expiry from ln S_T ~ N(ln S + (r - sigma^2 / 2) T, sigma^2 T), 400 000 times, finishes above the strike that often, within three standard errors. Hull's Example 15.6 (S0 = 42, K = 40, r = 10%, sigma = 20%, six months).

In [5]:

b = kernel(None, 42.0, 40.0, 0.10, 0.20, 0.5)
z = np.random.default_rng(3912).standard_normal(400_000)
end = 42.0 * np.exp((0.10 - 0.02) * 0.5 + 0.2 * np.sqrt(0.5) * z)
seen = np.mean(end > 40.0)
se = np.sqrt(seen * (1 - seen) / len(z))
assert abs(seen - ndtr(b)) <= 3 * se, (seen, ndtr(b))
print(f'N(d2) = {ndtr(b):.4f}; the draws finish in the money {seen:.4f} of the time')

b = kernel(None, 42.0, 40.0, 0.10, 0.20, 0.5) z = np.random.default_rng(3912).standard_normal(400_000) end = 42.0 * np.exp((0.10 - 0.02) * 0.5 + 0.2 * np.sqrt(0.5) * z) seen = np.mean(end > 40.0) se = np.sqrt(seen * (1 - seen) / len(z)) assert abs(seen - ndtr(b)) <= 3 * se, (seen, ndtr(b)) print(f'N(d2) = {ndtr(b):.4f}; the draws finish in the money {seen:.4f} of the time')

N(d2) = 0.7349; the draws finish in the money 0.7351 of the time