Black-Scholes: a European option's value¶

The Black-Scholes value of a European option on a stock paying no dividend, a call's c = S N(d1) - K e^(-rT) N(d2) and a put's p = K e^(-rT) N(-d2) - S N(-d1), in the money unit of S and K (the kernel is generic over it).

The parameters are lusine.finance.d1's; call chooses the call's value, else the put's. An option with no price, volatility or time left has no value here (d1 divides by nothing): its intrinsic value is the game's. Source: Black & Scholes, J. Political Economy 81 (1973); Hull, Options, Futures, and Other Derivatives, 9th ed. (Pearson, 2014), ch. 15, Example 15.6. 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, Example 15.6.

Proof

reference, invariant, limit.

option_value.trait

namespace lusine.finance;

float<_> option_value(bool call, float<_> spot, float<_> strike, float<year^-1> rate, float volatility,
                      float<year> period, float<year> years)
{
  let a = d1(spot, strike, rate, volatility, period, years);
  let b = d2(spot, strike, rate, volatility, period, years);
  let discounted = strike * exp(-rate * years);
  if (call) return spot * normal_cdf(x = a) - discounted * normal_cdf(x = b);
  return discounted * normal_cdf(x = -b) - spot * normal_cdf(x = -a);
}

Uses, without proving it here: finance.d1, finance.d2, finance.normal_cdf.

A European call's value S N(d1) - K e^(-rT) N(d2) and a put's K e^(-rT) N(-d2) - S N(-d1), N finance.normal_cdf's.

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/option_value.trait'))

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

// lusine.finance.option_value: the Black-Scholes value of a European option on a stock paying no dividend, a call's c =
// S N(d1) - K e^(-rT) N(d2) and a put's p = K e^(-rT) N(-d2) - S N(-d1), in the money unit of S and K (the kernel is
// generic over it). The parameters are lusine.finance.d1's; `call` chooses the call's value, else the put's. An option
// with no price, volatility or time left has no value here (d1 divides by nothing): its intrinsic value is the game's.
// Source: Black & Scholes, J. Political Economy 81 (1973); Hull, Options, Futures, and Other Derivatives, 9th ed.
// (Pearson, 2014), ch. 15, Example 15.6.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float<_> option_value(bool call, float<_> spot, float<_> strike, float<year^-1> rate, float volatility,
                      float<year> period, float<year> years)
{
  let a = d1(spot, strike, rate, volatility, period, years);
  let b = d2(spot, strike, rate, volatility, period, years);
  let discounted = strike * exp(-rate * years);
  if (call) return spot * normal_cdf(x = a) - discounted * normal_cdf(x = b);
  return discounted * normal_cdf(x = -b) - spot * normal_cdf(x = -a);
}

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 the value is read back in euros. 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 = True
BS = """namespace bs;

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

struct option {
  host bool call; 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<EUR> option_value;
};
lattice<option, 1> options(spacing = 1 m, boundary = wall);

void main()
{
  options.option_value = lusine.finance.option_value(options.call, 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 option_value 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 euros, 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/option_value.trait'], BS, {'bs.options': [n]})
    for field, v in (('call', call), ('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', 'option_value')])[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 np.where(c, S * ndtr(a) - df * ndtr(b), df * ndtr(-b) - S * ndtr(-a))


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 = True BS = """namespace bs; lusine.version version = { { number = 1 } }; struct option { host bool call; 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<EUR> option_value; }; lattice<option, 1> options(spacing = 1 m, boundary = wall); void main() { options.option_value = lusine.finance.option_value(options.call, 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 option_value 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 euros, 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/option_value.trait'], BS, {'bs.options': [n]}) for field, v in (('call', call), ('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', 'option_value')])[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 np.where(c, S * ndtr(a) - df * ndtr(b), df * ndtr(-b) - S * ndtr(-a)) 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 bool call; 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<EUR> option_value;
};
lattice<option, 1> options(spacing = 1 m, boundary = wall);

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

Reference: Hull's examples¶

Hull's printed values (Options, Futures, and Other Derivatives, 9th ed., 2014) are quoted from sources that were opened, the book itself not: Example 15.6 (S0 = 42, K = 40, r = 10%, sigma = 20%, six months) as QFRM's BS_Simple documentation cites it, and ch. 19's running example (S0 = 49, K = 50, r = 5%, sigma = 20%, 20 weeks) as the black-scholes-greeks README quotes it. Each is held to half a unit in the last digit printed.

In [4]:

c, p = kernel([True, False], 42.0, 40.0, 0.10, 0.20, 0.5)
g = kernel(True, 49.0, 50.0, 0.05, 0.20, 20 / 52)
assert abs(c - 4.76) <= 0.005 and abs(p - 0.81) <= 0.005, (c, p)
assert abs(g - 2.40) <= 0.005, g
print(f'Example 15.6: call {c:.4f} (4.76), put {p:.4f} (0.81); ch. 19: call {g:.4f} (2.40)')

c, p = kernel([True, False], 42.0, 40.0, 0.10, 0.20, 0.5) g = kernel(True, 49.0, 50.0, 0.05, 0.20, 20 / 52) assert abs(c - 4.76) <= 0.005 and abs(p - 0.81) <= 0.005, (c, p) assert abs(g - 2.40) <= 0.005, g print(f'Example 15.6: call {c:.4f} (4.76), put {p:.4f} (0.81); ch. 19: call {g:.4f} (2.40)')

Example 15.6: call 4.7594 (4.76), put 0.8086 (0.81); ch. 19: call 2.4005 (2.40)

Reference: scipy's Black-Scholes¶

A thousand options across moneyness, rates from -2% to 15%, volatilities from 5% to 90% and lives from a week to five years, calls and puts: each value within 2e-6 of its spot from the float64 account. Where the value is small beside the spot, a deep out-of-the-money option, its relative error is A&S 26.2.17's in the tail of N (finance.normal_cdf).

In [5]:

v = kernel(CALLS, SPOT, STRIKE, RATE, VOL, LIFE)
want = closed(CALLS, SPOT, STRIKE, RATE, VOL, LIFE)
err = np.abs(v - want) / SPOT
assert np.max(err) <= 2e-6, np.max(err)
big = want > 1e-3 * SPOT
print(f'largest error {np.max(err):.2e} of the spot; relative, where the value is above 1e-3 of it, '
      f'{np.max(np.abs(v[big] / want[big] - 1)):.2e}')

v = kernel(CALLS, SPOT, STRIKE, RATE, VOL, LIFE) want = closed(CALLS, SPOT, STRIKE, RATE, VOL, LIFE) err = np.abs(v - want) / SPOT assert np.max(err) <= 2e-6, np.max(err) big = want > 1e-3 * SPOT print(f'largest error {np.max(err):.2e} of the spot; relative, where the value is above 1e-3 of it, ' f'{np.max(np.abs(v[big] / want[big] - 1)):.2e}')

largest error 7.23e-07 of the spot; relative, where the value is above 1e-3 of it, 2.47e-04

Invariant: put-call parity, and in any money¶

c - p = S - K e^(-rT) for each option, to binary32's rounding of the larger of the two (a few EPS of the spot). The kernel is generic over the money unit: bound to an integer price in yen and a strike in yen, the same program gives the same value.

In [6]:

c = kernel(True, SPOT, STRIKE, RATE, VOL, LIFE)
p = kernel(False, SPOT, STRIKE, RATE, VOL, LIFE)
_, S, K, r, _, T, _ = held(True, SPOT, STRIKE, RATE, VOL, LIFE)
gap = np.abs(c - p - (S - K * np.exp(-r * T))) / SPOT
assert np.max(gap) <= 16 * EPS, np.max(gap)
YEN = BS.replace('host float<EUR> spot;', 'host integer<JPY> spot;').replace('float<EUR>', 'float<JPY>')
run = trait_run.Program(['finance/option_value.trait'], YEN, {'bs.options': [2]})
for field, value in (('call', [True, False]), ('spot', [4200, 4200]), ('strike', [4000.0, 4000.0]),
                     ('rate', 0.10 / YEAR), ('volatility', 0.20), ('period', YEAR), ('years', 0.5 * YEAR)):
    run.store('bs.options', field, np.broadcast_to(np.asarray(value), (2,)))
yen = run.steps('bs.main', 1, [('bs.options', 'option_value')])[0][0]
assert abs(yen[0] - 476.0) <= 0.5 and abs(yen[1] - 81.0) <= 0.5, yen
print(f'parity within {np.max(gap) / EPS:.1f} EPS of the spot; in yen, a call {yen[0]:.2f} and a put {yen[1]:.2f}')

c = kernel(True, SPOT, STRIKE, RATE, VOL, LIFE) p = kernel(False, SPOT, STRIKE, RATE, VOL, LIFE) _, S, K, r, _, T, _ = held(True, SPOT, STRIKE, RATE, VOL, LIFE) gap = np.abs(c - p - (S - K * np.exp(-r * T))) / SPOT assert np.max(gap) <= 16 * EPS, np.max(gap) YEN = BS.replace('host float<EUR> spot;', 'host integer<JPY> spot;').replace('float<EUR>', 'float<JPY>') run = trait_run.Program(['finance/option_value.trait'], YEN, {'bs.options': [2]}) for field, value in (('call', [True, False]), ('spot', [4200, 4200]), ('strike', [4000.0, 4000.0]), ('rate', 0.10 / YEAR), ('volatility', 0.20), ('period', YEAR), ('years', 0.5 * YEAR)): run.store('bs.options', field, np.broadcast_to(np.asarray(value), (2,))) yen = run.steps('bs.main', 1, [('bs.options', 'option_value')])[0][0] assert abs(yen[0] - 476.0) <= 0.5 and abs(yen[1] - 81.0) <= 0.5, yen print(f'parity within {np.max(gap) / EPS:.1f} EPS of the spot; in yen, a call {yen[0]:.2f} and a put {yen[1]:.2f}')

parity within 5.3 EPS of the spot; in yen, a call 475.94 and a put 80.86

Limit: the value against the spot¶

As the spot falls far below the strike a call is worth nothing and a put its discounted strike less the spot; far above, a call is worth the spot less the discounted strike and a put nothing; between, both lie above those bounds.

In [7]:

S = np.linspace(5.0, 150.0, 300)
c, p = kernel(True, S, 50.0, 0.05, 0.2, 1.0), kernel(False, S, 50.0, 0.05, 0.2, 1.0)
df = 50.0 * np.exp(-0.05)
assert np.all(c >= np.maximum(S - df, 0) - 1e-4) and np.all(p >= np.maximum(df - S, 0) - 1e-4)
assert c[0] < 1e-6 and abs(p[0] - (df - S[0])) < 1e-4 and abs(c[-1] - (S[-1] - df)) < 1e-4 and p[-1] < 1e-4
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(S, c, label='call')
ax.plot(S, p, label='put')
ax.plot(S, np.maximum(S - df, 0), ':', color='grey', label='bounds')
ax.plot(S, np.maximum(df - S, 0), ':', color='grey')
ax.set_xlabel('spot (EUR)')
ax.set_ylabel('value (EUR)')
ax.legend()
plt.show()

S = np.linspace(5.0, 150.0, 300) c, p = kernel(True, S, 50.0, 0.05, 0.2, 1.0), kernel(False, S, 50.0, 0.05, 0.2, 1.0) df = 50.0 * np.exp(-0.05) assert np.all(c >= np.maximum(S - df, 0) - 1e-4) and np.all(p >= np.maximum(df - S, 0) - 1e-4) assert c[0] < 1e-6 and abs(p[0] - (df - S[0])) < 1e-4 and abs(c[-1] - (S[-1] - df)) < 1e-4 and p[-1] < 1e-4 fig, ax = plt.subplots(figsize=(7, 3.5)) ax.plot(S, c, label='call') ax.plot(S, p, label='put') ax.plot(S, np.maximum(S - df, 0), ':', color='grey', label='bounds') ax.plot(S, np.maximum(df - S, 0), ':', color='grey') ax.set_xlabel('spot (EUR)') ax.set_ylabel('value (EUR)') ax.legend() plt.show()