Black-Scholes-Merton, Black-76 and the Greeks

"""M08 exercise: Black-Scholes-Merton, Black-76 and the Greeks.

Read notes §8.6-§8.9 first. Then run
    cd exercises && .venv/bin/python -m pytest modules/M08-option-pricing/tests/test_m08_black_scholes.py

Conventions used everywhere in M08:
  * r, q are continuously compounded rates per year; T is in years (ACT/365); sigma per sqrt(year).
  * vega is per 1.00 change in sigma (divide by 100 for "per vol point").
  * theta is dV/dt per year of calendar time (divide by 365 for "per day"). It is minus dV/dT.
  * rho is per 1.00 change in r.
  * The learning goal is Black-Scholes, not the error function: for N(x) use either
    scipy.stats.norm.cdf (imported below) or math.erfc, as you prefer.
  * The tests compare your functions with toolkit.options on a grid. Write your own
    first; reading toolkit/options.py before you have is reading the solution.
"""

from __future__ import annotations

import math

from scipy.stats import norm  # noqa: F401  (norm.cdf / norm.pdf are allowed)


def norm_cdf(x: float) -> float:
    """N(x) = P(Z <= x) for Z ~ N(0,1): norm.cdf(x), or 0.5 * erfc(-x / sqrt(2))."""
    raise NotImplementedError


def norm_pdf(x: float) -> float:
    """phi(x) = exp(-x^2/2) / sqrt(2 pi)."""
    raise NotImplementedError


def d1_d2(S: float, K: float, T: float, r: float, sigma: float, q: float = 0.0) -> tuple[float, float]:
    """Return (d1, d2) of notes eq. (8.22):

        d1 = [ln(S/K) + (r - q + sigma^2/2) T] / (sigma sqrt(T)),   d2 = d1 - sigma sqrt(T).

    Raise ValueError unless S, K, T, sigma are all > 0.
    Check: d1_d2(100, 100, 1, 0.05, 0.2) == (0.35, 0.15).
    """
    raise NotImplementedError


def bs_price(S: float, K: float, T: float, r: float, sigma: float, kind: str = "call", q: float = 0.0) -> float:
    """Black-Scholes-Merton price of a European option on an asset with continuous yield q (notes §8.7).

        call = S e^{-qT} N(d1) - K e^{-rT} N(d2)
        put  = K e^{-rT} N(-d2) - S e^{-qT} N(-d1)

    Limits you must handle (notes §8.7, Remark 8.7.3):
      * T <= 0: return intrinsic value max(S-K,0) or max(K-S,0).
      * sigma <= 0 (and T > 0): return max(+-(S e^{-qT} - K e^{-rT}), 0), the discounted
        intrinsic value of the forward.
    kind is 'call' or 'put'; anything else raises ValueError.
    Reference: bs_price(100, 100, 1, 0.05, 0.2) = 10.4505835722; put 5.5735260223.
    """
    raise NotImplementedError


def black76(F: float, K: float, T: float, r: float, sigma: float, kind: str = "call") -> float:
    """Black (1976) price of a European option on a futures/forward price F (notes §8.8):

        call = e^{-rT} [F N(d1) - K N(d2)],  put = e^{-rT} [K N(-d2) - F N(-d1)],
        d1 = [ln(F/K) + sigma^2 T / 2] / (sigma sqrt(T)),  d2 = d1 - sigma sqrt(T).

    If T <= 0 or sigma <= 0 return e^{-rT} times the intrinsic value in F.
    Reference: black76(20, 20, 4/12, 0.09, 0.25, 'put') = 1.116641 (Hull).
    """
    raise NotImplementedError


def bs_greeks(S: float, K: float, T: float, r: float, sigma: float, kind: str = "call", q: float = 0.0) -> dict:
    """Analytic Greeks of bs_price (notes §8.9, Table 8.2). Return a dict with keys
    'delta', 'gamma', 'vega', 'theta', 'rho' (units as in the module docstring).

        delta_call = e^{-qT} N(d1),           delta_put = -e^{-qT} N(-d1)
        gamma      = e^{-qT} phi(d1) / (S sigma sqrt T)
        vega       = S e^{-qT} phi(d1) sqrt T
        theta_call = -S e^{-qT} phi(d1) sigma / (2 sqrt T) - r K e^{-rT} N(d2) + q S e^{-qT} N(d1)
        theta_put  = -S e^{-qT} phi(d1) sigma / (2 sqrt T) + r K e^{-rT} N(-d2) - q S e^{-qT} N(-d1)
        rho_call   = K T e^{-rT} N(d2),       rho_put = -K T e^{-rT} N(-d2)

    Reference (S=K=100, T=1, r=5%, sigma=20%, call): delta 0.636831, gamma 0.018762,
    vega 37.524035, theta -6.414028, rho 53.232482.
    """
    raise NotImplementedError


def parity_residual(C: float, P: float, S: float, K: float, T: float, r: float, q: float = 0.0) -> float:
    """Put-call parity residual C - P - (S e^{-qT} - K e^{-rT}) (notes Proposition 8.2.1). Zero under no arbitrage."""
    raise NotImplementedError


def parity_residual_futures(C: float, P: float, F: float, K: float, T: float, r: float) -> float:
    """Parity on the forward: C - P - e^{-rT} (F - K) (notes eq. (8.5)). The practical Indian form:
    take F from the futures market and no dividend estimate is needed."""
    raise NotImplementedError