Bump-and-Reval — Numerical Greeks and Convergence Order

Companion notebook to the Bump-and-Reval article (Greeks & Sensitivities, Module 5). Standard library only.

§1 — Forward difference delta: O(h) convergence demonstration
§2 — Central difference delta: O(h²) convergence and V-shape
§3 — Central difference gamma: O(h²) convergence
§4 — Mixed partial (vanna) by 4-point stencil
§5 — Validation summary: convergence orders verified

💡 Run cells in order. §2–§5 depend on §1 helpers.

§0 — Analytic Greeks (Reference Implementations)

import math

EPS_MACH = 2.220446049250313e-16  # IEEE 754 double machine epsilon

def norm_cdf(x):
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))

def norm_pdf(x):
    return math.exp(-0.5 * x * x) / math.sqrt(2.0 * math.pi)

def bs_d1_d2(S, K, r, sigma, T):
    d1 = (math.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * math.sqrt(T))
    return d1, d1 - sigma * math.sqrt(T)

def bs_price(S, K, r, sigma, T, opt="call"):
    """Reference Black-Scholes pricer — the function we bump and reval."""
    d1, d2 = bs_d1_d2(S, K, r, sigma, T)
    if opt == "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)

# Analytic Greeks (reference / ground truth)
def analytic_delta(S, K, r, sigma, T, opt="call"):
    d1, _ = bs_d1_d2(S, K, r, sigma, T)
    return norm_cdf(d1) if opt == "call" else norm_cdf(d1) - 1.0

def analytic_gamma(S, K, r, sigma, T):
    d1, _ = bs_d1_d2(S, K, r, sigma, T)
    return norm_pdf(d1) / (S * sigma * math.sqrt(T))

def analytic_vega(S, K, r, sigma, T):
    d1, _ = bs_d1_d2(S, K, r, sigma, T)
    return S * math.sqrt(T) * norm_pdf(d1)

def analytic_vanna(S, K, r, sigma, T):
    d1, d2 = bs_d1_d2(S, K, r, sigma, T)
    return -norm_pdf(d1) * d2 / sigma

# Numerical Greek approximations
def fwd_diff_delta(S, K, r, sigma, T, h):
    """Forward difference: O(h) truncation error."""
    return (bs_price(S + h, K, r, sigma, T) - bs_price(S, K, r, sigma, T)) / h

def central_diff_delta(S, K, r, sigma, T, h):
    """Central difference: O(h²) truncation error."""
    return (bs_price(S + h, K, r, sigma, T) - bs_price(S - h, K, r, sigma, T)) / (2 * h)

def central_diff_gamma(S, K, r, sigma, T, h):
    """Second-order central difference for gamma: O(h²)."""
    v0 = bs_price(S, K, r, sigma, T)
    vp = bs_price(S + h, K, r, sigma, T)
    vm = bs_price(S - h, K, r, sigma, T)
    return (vp - 2 * v0 + vm) / h**2

def mixed_diff_vanna(S, K, r, sigma, T, h_s, h_v):
    """4-point stencil for mixed partial d²V/(dS dsigma): O(h_s², h_v²)."""
    vpp = bs_price(S + h_s, K, r, sigma + h_v, T)
    vpm = bs_price(S + h_s, K, r, sigma - h_v, T)
    vmp = bs_price(S - h_s, K, r, sigma + h_v, T)
    vmm = bs_price(S - h_s, K, r, sigma - h_v, T)
    return (vpp - vpm - vmp + vmm) / (4 * h_s * h_v)

print("✓ Analytic and numerical Greek functions loaded.")
S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
print(f"\nReference: S=K=100, r=0.05, sigma=0.20, T=1.0")
print(f"  Analytic delta = {analytic_delta(S,K,r,sigma,T):.8f}")
print(f"  Analytic gamma = {analytic_gamma(S,K,r,sigma,T):.8f}")
print(f"  Analytic vega  = {analytic_vega(S,K,r,sigma,T):.8f}")
print(f"  Analytic vanna = {analytic_vanna(S,K,r,sigma,T):.8f}")
print(f"  Machine epsilon: {EPS_MACH:.3e}")

§1 — Forward Difference Delta: O(h) Convergence

Forward difference: Δ_FD = [V(S+h) − V(S)] / h. Truncation error = ½Γh + O(h²). Doubling h doubles the error — first-order convergence.

S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
delta_true = analytic_delta(S, K, r, sigma, T)
gamma_true = analytic_gamma(S, K, r, sigma, T)

print(f"True delta = {delta_true:.10f}")
print()
print("Forward difference convergence (O(h)):")
print(f"{'h':>10}  {'Delta_FD':>12}  {'|Error|':>12}  {'log10(err)':>12}  {'Theoretical ½Γh':>18}")
print("-" * 72)

bump_sizes = [1.0, 0.5, 0.1, 0.05, 0.01, 0.005, 0.001, 1e-4, 1e-6, 1e-8]
prev_err = None
for h in bump_sizes:
    d_fd = fwd_diff_delta(S, K, r, sigma, T, h)
    err = abs(d_fd - delta_true)
    log_err = math.log10(err) if err > 0 else -999
    theory = 0.5 * gamma_true * h  # leading O(h) term
    ratio_str = ""
    if prev_err is not None and err > 0:
        ratio = prev_err / err
        ratio_str = f" (ratio {ratio:.1f})"
    print(f"{h:>10.2e}  {d_fd:>12.8f}  {err:>12.2e}  {log_err:>12.2f}  {theory:>18.2e}{ratio_str}")
    prev_err = err

print()
print("At large h: error ≈ ½Γh — first-order (ratio ≈ 10 for each 10x h reduction).")
print("At very small h (< 1e-7): floating-point cancellation dominates — error increases.")

§2 — Central Difference Delta: O(h²) Convergence and the V-Shape

Central difference: Δ_CD = [V(S+h) − V(S−h)] / (2h). Truncation error = (h²/6)·V''' + O(h⁴). Error decreases as h² for medium h, then increases due to floating-point cancellation. The optimal h is near √(ε_mach) ≈ 1.5×10⁻⁸ for forward difference, and (ε_mach)^(1/3) ≈ 6×10⁻⁶ for central difference.

S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
delta_true = analytic_delta(S, K, r, sigma, T)

print(f"True delta = {delta_true:.10f}")
print()
print("Central difference convergence (O(h²)):")
print(f"{'h':>10}  {'Delta_CD':>12}  {'|Error|':>12}  {'log10(err)':>12}")
print("-" * 54)

bump_sizes_cd = [1.0, 0.5, 0.1, 0.05, 0.01, 0.005, 0.001, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-10, 1e-12]
prev_err_cd = None
for h in bump_sizes_cd:
    d_cd = central_diff_delta(S, K, r, sigma, T, h)
    err = abs(d_cd - delta_true)
    log_err = math.log10(err) if err > 1e-16 else -16
    ratio_str = ""
    if prev_err_cd is not None and prev_err_cd > 0 and err > 0:
        ratio = prev_err_cd / err
        ratio_str = f" (ratio {ratio:.1f})"
    print(f"{h:>10.2e}  {d_cd:>12.8f}  {err:>12.2e}  {log_err:>12.2f}{ratio_str}")
    prev_err_cd = err

print()
print("Truncation regime: error ∝ h² (ratio ≈ 100 per 10x h reduction)")
print("Optimal h ≈ ε_mach^(1/3) ≈ 6e-6: minimum error before cancellation dominates.")
print("Cancellation regime (h < 1e-7): error grows as ε_mach/h.")

§3 — Numerical Gamma: O(h²) Convergence

Gamma: [V(S+h) − 2V(S) + V(S−h)] / h². Error = (h²/12)·V'''' + O(h⁴). The optimal h for gamma is larger than for delta: h_opt ≈ ε_mach^(1/4) ≈ 1.2×10⁻⁴.

S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
gamma_true = analytic_gamma(S, K, r, sigma, T)

print(f"True gamma = {gamma_true:.10f}")
print()
print("Central difference GAMMA convergence (O(h²)):")
print(f"{'h':>10}  {'Gamma_FD':>12}  {'|Error|':>12}  {'log10(err)':>12}")
print("-" * 54)

bump_sizes_g = [1.0, 0.5, 0.1, 0.01, 0.001, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8]
prev_err_g = None
for h in bump_sizes_g:
    g_fd = central_diff_gamma(S, K, r, sigma, T, h)
    err = abs(g_fd - gamma_true)
    log_err = math.log10(err) if err > 1e-16 else -16
    ratio_str = ""
    if prev_err_g is not None and prev_err_g > 0 and err > 0:
        ratio = prev_err_g / err
        ratio_str = f" (ratio {ratio:.1f})"
    print(f"{h:>10.2e}  {g_fd:>12.8f}  {err:>12.2e}  {log_err:>12.2f}{ratio_str}")
    prev_err_g = err

# Theoretical optimal h for gamma: eps_mach^(1/4)
h_opt_gamma = EPS_MACH**(1.0/4.0)
print(f"\nTheoretical optimal h for gamma: eps_mach^(1/4) = {h_opt_gamma:.3e}")
g_at_opt = central_diff_gamma(S, K, r, sigma, T, h_opt_gamma)
print(f"Error at h_opt = {abs(g_at_opt - gamma_true):.3e}")

§4 — Mixed Partial: Vanna by 4-Point Stencil

Vanna = ∂²V/(∂S∂σ) via the 4-point stencil:
[V(S+h,σ+k) − V(S+h,σ−k) − V(S−h,σ+k) + V(S−h,σ−k)] / (4hk)
Accuracy: O(h², k²). Requires 4 extra pricings (plus base if needed).

S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
vanna_true = analytic_vanna(S, K, r, sigma, T)

print(f"True vanna = {vanna_true:.8f}")
print()
print("Mixed partial (vanna) convergence — 4-point stencil:")
print(f"{'h_s':>8} {'h_v':>8} {'Vanna_FD':>12} {'|Error|':>12}")
print("-" * 46)

stencil_params = [
    (1.0, 0.01),
    (0.5, 0.005),
    (0.1, 0.001),
    (0.05, 0.0005),
    (0.01, 0.0001),
    (0.001, 0.00001),
]
for h_s, h_v in stencil_params:
    vanna_fd = mixed_diff_vanna(S, K, r, sigma, T, h_s, h_v)
    err = abs(vanna_fd - vanna_true)
    print(f"{h_s:>8.3f} {h_v:>8.5f} {vanna_fd:>12.6f} {err:>12.2e}")

print()
print("4-point stencil achieves O(h_s², h_v²) — same accuracy as central differences.")
print("Cost: 4 extra pricings (vs 2 for single-variable central difference).")

§5 — Validation Summary: Convergence Orders Verified

S, K, r, sigma, T = 100.0, 100.0, 0.05, 0.20, 1.0
results = []

delta_true = analytic_delta(S, K, r, sigma, T)
gamma_true = analytic_gamma(S, K, r, sigma, T)
vanna_true = analytic_vanna(S, K, r, sigma, T)

# 1. Forward difference is O(h): error ratio ~10 when h decreases 10x
h1, h2 = 0.01, 0.001
err1_fd = abs(fwd_diff_delta(S, K, r, sigma, T, h1) - delta_true)
err2_fd = abs(fwd_diff_delta(S, K, r, sigma, T, h2) - delta_true)
ratio_fd = err1_fd / err2_fd if err2_fd > 0 else 0
results.append((f"FD delta O(h): error ratio h1/h2=10 gives ~10 (got {ratio_fd:.1f})", 5 < ratio_fd < 15))

# 2. Central difference is O(h²): error ratio ~100 when h decreases 10x
err1_cd = abs(central_diff_delta(S, K, r, sigma, T, h1) - delta_true)
err2_cd = abs(central_diff_delta(S, K, r, sigma, T, h2) - delta_true)
ratio_cd = err1_cd / err2_cd if err2_cd > 0 else 0
results.append((f"CD delta O(h²): error ratio h1/h2=10 gives ~100 (got {ratio_cd:.1f})", 50 < ratio_cd < 200))

# 3. Central difference is more accurate than forward difference at same h
h = 0.01
err_fd = abs(fwd_diff_delta(S, K, r, sigma, T, h) - delta_true)
err_cd = abs(central_diff_delta(S, K, r, sigma, T, h) - delta_true)
results.append(f"CD more accurate than FD at h={h}: err_FD={err_fd:.2e} > err_CD={err_cd:.2e}")
results[-1] = (results[-1], err_cd < err_fd)

# 4. Gamma FD is O(h²)
err1_g = abs(central_diff_gamma(S, K, r, sigma, T, 0.1) - gamma_true)
err2_g = abs(central_diff_gamma(S, K, r, sigma, T, 0.01) - gamma_true)
ratio_g = err1_g / err2_g if err2_g > 0 else 0
results.append((f"CD gamma O(h²): error ratio ~100 (got {ratio_g:.1f})", 50 < ratio_g < 300))

# 5. Vanna 4-point stencil converges to analytic at h=0.01, k=0.001
vanna_fd = mixed_diff_vanna(S, K, r, sigma, T, 0.01, 0.001)
err_v = abs(vanna_fd - vanna_true)
results.append((f"Vanna 4-point stencil: error={err_v:.2e} < 1e-4", err_v < 1e-4))

# 6. At very small h (1e-12), forward difference error is larger than at h=0.001 (cancellation)
err_small_h = abs(fwd_diff_delta(S, K, r, sigma, T, 1e-12) - delta_true)
err_medium_h = abs(fwd_diff_delta(S, K, r, sigma, T, 0.001) - delta_true)
results.append((f"Cancellation: err at h=1e-12 ({err_small_h:.2e}) > err at h=0.001 ({err_medium_h:.2e})",
                err_small_h > err_medium_h))

print("=" * 68)
print("VALIDATION SUMMARY — Bump-and-Reval Numerical Greeks Module")
print("=" * 68)
for label, passed in results:
    print(f"  {'✓' if passed else '✗'}  {label}")
n_pass = sum(1 for _, p in results if p)
print("=" * 68)
print(f"  {n_pass}/{len(results)} checks passed.")
print()
print("Key conclusions:")
print("  1. Central difference is O(h²); forward difference is O(h).")
print("  2. Gamma FD uses V(S+h) − 2V(S) + V(S−h) with the same 2 extra pricings as delta.")
print("  3. Mixed partials require 4 extra pricings — O(h², k²) accuracy.")
print("  4. Floating-point cancellation sets a floor on accuracy for very small h.")