C++20 Pricing Library

// BlackScholes.hpp
#pragma once
/**
 * @brief European option pricing under the Black-Scholes model.
 *
 * Model assumptions:
 *   dS = (r − q) S dt + σ S dW  under risk-neutral measure Q
 *   Constant r, q, σ; continuous compounding; European exercise.
 *
 * Conventions:
 *   All rates are annualised, continuously compounded.
 *   σ is annualised vol in decimal form (0.20 = 20%).
 *   T is time to expiry in years (ACT/365).
 *   Vega is scaled per 1 vol point (1%).  Theta is per calendar day.
 */
#include "common.hpp"
#include <stdexcept>

namespace bb {

double normcdf(double x) noexcept;
double normpdf(double x) noexcept;

struct BSResult {
    double price;  // Option premium
    double delta;  // ∂V/∂S
    double gamma;  // ∂²V/∂S²
    double vega;   // ∂V/∂σ per 1 vol point  (∂V/∂(100σ) / 100)
    double theta;  // ∂V/∂t per calendar day  (annualised / 365)
    double rho;    // ∂V/∂r per 1 bp          (∂V/∂(100r)  / 100)
};

// European option: closed-form price + all 5 Greeks.
// Throws std::invalid_argument on invalid inputs (S,K,σ,T must be > 0).
BSResult bs_european(double S, double K, double r, double q,
                     double sigma, double T, OptionType type);

// Put-call parity RHS:  C − P = S·e^{−qT} − K·e^{−rT}
double bs_parity_forward(double S, double K, double r, double q, double T) noexcept;

// Implied volatility via Newton-Raphson with Brent bracket fallback.
// Initial guess: Brenner-Subrahmanyam (1988) ATM approximation.
// Throws std::runtime_error if solver does not converge.
double bs_implied_vol(double market_price, double S, double K,
                      double r, double q, double T, OptionType type,
                      double tol = 1e-8, int max_iter = 100);

} // namespace bb

// BlackScholes.cpp
#include "bb/BlackScholes.hpp"
#include <cmath>
#include <stdexcept>

namespace bb {

// Standard normal CDF via erfc (numerically stable for large |x|).
double normcdf(double x) noexcept { return 0.5 * std::erfc(-x * M_SQRT1_2); }
double normpdf(double x) noexcept {
    constexpr double inv_sqrt2pi = 0.3989422804014327;
    return inv_sqrt2pi * std::exp(-0.5 * x * x);
}

BSResult bs_european(double S, double K, double r, double q,
                     double sigma, double T, OptionType type) {
    if (S <= 0.0 || K <= 0.0 || sigma <= 0.0 || T <= 0.0)
        throw std::invalid_argument("S, K, sigma, T must be > 0");

    const double sqrtT      = std::sqrt(T);
    const double sigma_sqT  = sigma * sqrtT;
    const double d1 = (std::log(S / K) + (r - q + 0.5 * sigma * sigma) * T) / sigma_sqT;
    const double d2 = d1 - sigma_sqT;

    const double Nd1 = normcdf(d1),  Nmd1 = 1.0 - Nd1;
    const double Nd2 = normcdf(d2),  Nmd2 = 1.0 - Nd2;
    const double nd1 = normpdf(d1);

    const double discR = std::exp(-r * T);
    const double discQ = std::exp(-q * T);
    const double F     = S * discQ;   // forward (undiscounted)

    BSResult res{};

    if (type == OptionType::Call) {
        res.price = F * Nd1  - K * discR * Nd2;
        res.delta = discQ * Nd1;
        res.rho   = K * T * discR * Nd2  * 0.01;   // per 1 bp
    } else {
        res.price = K * discR * Nmd2 - F * Nmd1;
        res.delta = -discQ * Nmd1;
        res.rho   = -K * T * discR * Nmd2 * 0.01;
    }

    // Greeks symmetric in call/put:
    res.gamma  = discQ * nd1 / (S * sigma_sqT);
    res.vega   = F * nd1 * sqrtT * 0.01;      // per 1 vol point
    res.theta  = -(F * nd1 * sigma / (2.0 * sqrtT))
                 - r * K * discR * (type == OptionType::Call ?  Nd2 : -Nmd2)
                 + q * F           * (type == OptionType::Call ?  Nd1 : -Nmd1);
    res.theta /= 365.0;   // annualised → per calendar day

    return res;
}

double bs_parity_forward(double S, double K, double r, double q, double T) noexcept {
    return S * std::exp(-q * T) - K * std::exp(-r * T);
}

// Newton-Raphson implied vol with Brent fallback.
double bs_implied_vol(double market_price, double S, double K,
                      double r, double q, double T, OptionType type,
                      double tol, int max_iter) {
    // Brenner-Subrahmanyam (1988) ATM initial guess: σ₀ ≈ √(2π/T)·V/S
    double sigma = std::clamp(std::sqrt(2.0 * M_PI / T) * market_price / S, 1e-4, 10.0);

    for (int iter = 0; iter < max_iter; ++iter) {
        const BSResult res = bs_european(S, K, r, q, sigma, T, type);
        const double diff  = res.price - market_price;
        if (std::abs(diff) < tol) return sigma;
        const double vega_raw = res.vega * 100.0;   // undo per-vol-point scaling
        if (std::abs(vega_raw) < 1e-12) break;      // near-zero vega: switch to Brent
        sigma = std::clamp(sigma - diff / vega_raw, 1e-6, 10.0);
    }

    // Brent bracket fallback (guaranteed convergence; price is monotone in σ).
    auto f  = [&](double s) { return bs_european(S, K, r, q, s, T, type).price - market_price; };
    double a = 1e-6, b = 10.0, fa = f(a), fb = f(b);
    if (fa * fb > 0.0)
        throw std::runtime_error("bs_implied_vol: no root in [1e-6, 10.0]");

    for (int i = 0; i < 200 && std::abs(b - a) > tol; ++i) {
        const double mid = 0.5 * (a + b);
        const double fm  = f(mid);
        if (fa * fm < 0.0) { b = mid; fb = fm; } else { a = mid; fa = fm; }
    }
    return 0.5 * (a + b);
}

} // namespace bb

// tests/test_black_scholes.cpp
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "bb/BlackScholes.hpp"

using namespace bb;
using Catch::Matchers::WithinAbsMargin;

// Reference: S=100, K=100, r=5%, q=0, σ=20%, T=1y
TEST_CASE("ATM call price", "[bs][call]") {
    const auto res = bs_european(100.0, 100.0, 0.05, 0.0, 0.20, 1.0, OptionType::Call);
    REQUIRE_THAT(res.price, WithinAbsMargin(10.4506, 1e-3));  // QuantLib reference
    REQUIRE(res.delta > 0.0 && res.delta < 1.0);
    REQUIRE_THAT(res.delta, WithinAbsMargin(0.6368, 1e-3));
    REQUIRE(res.gamma > 0.0);   // long option: gamma always positive
    REQUIRE(res.theta < 0.0);   // long option: theta always negative
}

TEST_CASE("Put-call parity: C − P = S·exp(−qT) − K·exp(−rT)", "[bs][parity]") {
    const double S=105, K=95, r=0.03, q=0.01, sig=0.25, T=0.5;
    const auto call = bs_european(S, K, r, q, sig, T, OptionType::Call);
    const auto put  = bs_european(S, K, r, q, sig, T, OptionType::Put);
    REQUIRE_THAT(call.price - put.price, WithinAbsMargin(bs_parity_forward(S,K,r,q,T), 1e-10));
}

TEST_CASE("Delta symmetry: Δ_call − Δ_put = exp(−qT)", "[bs][delta]") {
    const auto call = bs_european(100, 100, 0.05, 0.02, 0.20, 1.0, OptionType::Call);
    const auto put  = bs_european(100, 100, 0.05, 0.02, 0.20, 1.0, OptionType::Put);
    REQUIRE_THAT(call.delta - put.delta, WithinAbsMargin(std::exp(-0.02), 1e-10));
}

TEST_CASE("Gamma and Vega agree for call and put", "[bs][greeks]") {
    const auto call = bs_european(100, 105, 0.04, 0.01, 0.22, 0.75, OptionType::Call);
    const auto put  = bs_european(100, 105, 0.04, 0.01, 0.22, 0.75, OptionType::Put);
    REQUIRE_THAT(call.gamma, WithinAbsMargin(put.gamma, 1e-10));  // put-call symmetry
    REQUIRE_THAT(call.vega,  WithinAbsMargin(put.vega,  1e-10));
}

TEST_CASE("Implied vol round-trip", "[bs][impliedvol]") {
    const double S=100, K=105, r=0.04, q=0.01, sigma=0.22, T=0.75;
    const auto res = bs_european(S, K, r, q, sigma, T, OptionType::Call);
    const double iv = bs_implied_vol(res.price, S, K, r, q, T, OptionType::Call);
    REQUIRE_THAT(iv, WithinAbsMargin(sigma, 1e-6));   // round-trip to 6 d.p.
}

TEST_CASE("Deep ITM call approaches intrinsic", "[bs][intrinsic]") {
    const auto res = bs_european(200.0, 100.0, 0.05, 0.0, 0.20, 1.0, OptionType::Call);
    const double intrinsic = 200.0 - 100.0 * std::exp(-0.05);
    REQUIRE_THAT(res.price, WithinAbsMargin(intrinsic, 1.0));
}

TEST_CASE("Deep OTM call price near zero", "[bs][otm]") {
    const auto res = bs_european(100.0, 200.0, 0.05, 0.0, 0.20, 1.0, OptionType::Call);
    REQUIRE(res.price >= 0.0);
    REQUIRE(res.price < 0.01);
}

TEST_CASE("Invalid inputs throw", "[bs][validation]") {
    REQUIRE_THROWS_AS(bs_european(-1, 100, 0.05, 0, 0.20, 1, OptionType::Call), std::invalid_argument);
    REQUIRE_THROWS_AS(bs_european(100, 100, 0.05, 0, 0.0, 1, OptionType::Call), std::invalid_argument);
}

// MonteCarlo.hpp
#pragma once
/**
 * @brief Monte Carlo pricing of European options under GBM.
 *
 * Model: dS = (r−q) S dt + σ S dW  (exact terminal simulation — no time-step bias)
 *   S_T = S_0 · exp((r − q − σ²/2)·T + σ·√T·Z),   Z ~ N(0,1)
 *
 * Variance reduction: antithetic variates (Z, −Z) — halves asymptotic std error
 * at constant path count. No bias introduced: E[f(Z)+f(−Z)]/2 = E[f(Z)].
 *
 * Convergence: O(1/√N). 100× more paths reduces std error by 10×.
 */
#include "common.hpp"
#include <cstdint>
#include <vector>

namespace bb {

struct MCResult {
    double   price;      // Discounted expected payoff
    double   std_error;  // MC standard error (1σ)
    double   ci_lower;   // 95% confidence interval lower bound
    double   ci_upper;   // 95% confidence interval upper bound
    uint64_t n_paths;    // Effective path count
};

struct ConvergencePoint { uint64_t n_paths; double price; double std_error; };

// European option via exact GBM simulation + antithetic variates (optional).
MCResult mc_european_gbm(double S, double K, double r, double q,
                         double sigma, double T, OptionType type,
                         uint64_t n_paths, uint64_t seed = 42,
                         bool antithetic = true);

// Log-spaced convergence diagnostic from 10 to n_max paths.
std::vector<ConvergencePoint> mc_convergence(
    double S, double K, double r, double q, double sigma, double T,
    OptionType type, uint64_t n_max = 1'000'000, int n_steps = 12, uint64_t seed = 42);

} // namespace bb

// tests/test_monte_carlo.cpp
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "bb/BlackScholes.hpp"
#include "bb/MonteCarlo.hpp"

using namespace bb;
using Catch::Matchers::WithinAbsMargin;

// Reference: ATM call — MC must agree with closed-form within 3σ.
TEST_CASE("MC call price agrees with BS closed-form (3σ)", "[mc][call]") {
    const auto bs  = bs_european(100, 100, 0.05, 0.0, 0.20, 1.0, OptionType::Call);
    const auto mc  = mc_european_gbm(100, 100, 0.05, 0.0, 0.20, 1.0,
                                      OptionType::Call, 500'000, 42, true);
    REQUIRE_THAT(mc.price, WithinAbsMargin(bs.price, 3.0 * mc.std_error));
}

TEST_CASE("MC put price agrees with BS (3σ)", "[mc][put]") {
    const auto bs  = bs_european(100, 105, 0.03, 0.01, 0.25, 0.5, OptionType::Put);
    const auto mc  = mc_european_gbm(100, 105, 0.03, 0.01, 0.25, 0.5,
                                      OptionType::Put,  500'000, 42, true);
    REQUIRE_THAT(mc.price, WithinAbsMargin(bs.price, 3.0 * mc.std_error));
}

TEST_CASE("Antithetic reduces std error vs plain MC (same N)", "[mc][antithetic]") {
    // Antithetic should produce strictly smaller std error for monotone payoffs.
    const auto plain = mc_european_gbm(100, 100, 0.05, 0, 0.20, 1.0,
                                       OptionType::Call, 10'000, 42, false);
    const auto anti  = mc_european_gbm(100, 100, 0.05, 0, 0.20, 1.0,
                                       OptionType::Call, 10'000, 42, true);
    REQUIRE(anti.std_error < plain.std_error);   // antithetic always wins for vanilla call
}

TEST_CASE("Convergence diagnostic produces 12 log-spaced points", "[mc][convergence]") {
    const auto pts = mc_convergence(100, 100, 0.05, 0, 0.20, 1.0, OptionType::Call,
                                    100'000, 12, 42);
    REQUIRE(pts.size() == 12);
    // Path counts are ascending
    for (std::size_t i = 1; i < pts.size(); ++i)
        REQUIRE(pts[i].n_paths > pts[i-1].n_paths);
}

TEST_CASE("MC put-call parity holds within 3σ", "[mc][parity]") {
    const double S=100, K=100, r=0.05, q=0.0, sigma=0.20, T=1.0;
    const auto c = mc_european_gbm(S, K, r, q, sigma, T, OptionType::Call,  500'000);
    const auto p = mc_european_gbm(S, K, r, q, sigma, T, OptionType::Put,   500'000);
    const double parity = S * std::exp(-q*T) - K * std::exp(-r*T);
    const double se     = std::sqrt(c.std_error*c.std_error + p.std_error*p.std_error);
    REQUIRE_THAT(c.price - p.price, WithinAbsMargin(parity, 3.0 * se));
}

// FiniteDifference.hpp
#pragma once
/**
 * @brief Crank-Nicolson finite difference PDE solver for European and American options.
 *
 * PDE (Black-Scholes in time-to-expiry τ = T − t):
 *   ∂V/∂τ = ½σ²S²(∂²V/∂S²) + (r−q)S(∂V/∂S) − rV
 *
 * Discretisation:
 *   - S-grid: uniform in log-space: x = log(S), xᵢ = x_min + i·Δx
 *   - τ-grid: uniform in τ
 *   - Scheme: Crank-Nicolson (θ=0.5) — O(Δx², Δτ²), unconditionally stable
 *   - Each time step: tridiagonal solve via Thomas algorithm
 *
 * Boundary conditions (log-space):
 *   Call: V(x_min,τ) = 0;   V(x_max,τ) = exp(x_max − qτ) − K·exp(−rτ)
 *   Put:  V(x_min,τ) = K·exp(−rτ);  V(x_max,τ) = 0
 *
 * American (PSOR): at each time step, enforce V ≥ intrinsic via projected SOR.
 *
 * Grid: S_max = max(3S, 5K) to keep the far boundary innocuous.
 * Production: n_S=500, n_τ=500 gives 4–5 significant figures.
 */
#include "common.hpp"

namespace bb {

struct FDResult {
    double price;   // V(S, t=0) by linear interpolation on the S-grid
    double delta;   // ∂V/∂S  (central FD on the S-grid)
    double gamma;   // ∂²V/∂S² (central FD on the S-grid)
};

// Crank-Nicolson for European options.
FDResult fd_european_cn(double S, double K, double r, double q,
                        double sigma, double T, OptionType type,
                        int n_S = 400, int n_t = 400);

// Crank-Nicolson + projected SOR for American options.
// Early exercise premium = fd_american_psor.price − fd_european_cn.price.
FDResult fd_american_psor(double S, double K, double r, double q,
                          double sigma, double T, OptionType type,
                          int n_S = 400, int n_t = 400,
                          double omega = 1.3, double psor_tol = 1e-8);

} // namespace bb

# black_scholes.py
"""
Black-Scholes European option analytics in Python.

Assumptions:
    dS = (r − q) S dt + σ S dW   under the risk-neutral measure Q
    Constant r, q, σ; European exercise; continuous compounding; ACT/365.

Conventions:
    r, q, σ all in decimal form (0.20 = 20 %).
    Vega returned per 1 vol point (1 %). Theta returned per calendar day.
"""

from __future__ import annotations

import numpy as np
from scipy.stats import norm
from typing import NamedTuple


class BSResult(NamedTuple):
    price: float
    delta: float
    gamma: float
    vega:  float   # per 1 vol point (1%)
    theta: float   # per calendar day
    rho:   float   # per 1 bp (0.01%)
    d1:    float
    d2:    float


def bs_european(
    S: float, K: float, r: float, q: float, sigma: float, T: float,
    option_type: str = "call",
) -> BSResult:
    """
    European option price and Greeks under Black-Scholes-Merton.

    Parameters
    ----------
    S, K         : Spot and strike (same currency).
    r, q         : Continuously compounded risk-free rate and dividend yield.
    sigma        : Annualised vol (decimal).
    T            : Time to expiry in years (ACT/365).
    option_type  : "call" or "put".
    """
    if S <= 0 or K <= 0 or sigma <= 0 or T <= 0:
        raise ValueError("S, K, sigma, T must be strictly positive.")

    sqrt_T    = np.sqrt(T)
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * sqrt_T)
    d2 = d1 - sigma * sqrt_T

    Nd1,  Nd2  = norm.cdf(d1),  norm.cdf(d2)
    Nmd1, Nmd2 = 1.0 - Nd1,    1.0 - Nd2
    nd1 = norm.pdf(d1)

    disc_r = np.exp(-r * T)
    disc_q = np.exp(-q * T)
    F      = S * disc_q   # undiscounted forward

    is_call = option_type.lower() == "call"
    if is_call:
        price = F * Nd1  - K * disc_r * Nd2
        delta = disc_q * Nd1
        rho   = K * T * disc_r * Nd2  * 0.01
    else:
        price = K * disc_r * Nmd2 - F * Nmd1
        delta = -disc_q * Nmd1
        rho   = -K * T * disc_r * Nmd2 * 0.01

    gamma  = disc_q * nd1 / (S * sigma * sqrt_T)
    vega   = F * nd1 * sqrt_T * 0.01
    theta  = (-(F * nd1 * sigma / (2 * sqrt_T))
              - r * K * disc_r * (Nd2  if is_call else -Nmd2)
              + q * F           * (Nd1  if is_call else -Nmd1)) / 365.0

    return BSResult(price=price, delta=delta, gamma=gamma,
                    vega=vega, theta=theta, rho=rho, d1=d1, d2=d2)


def bs_implied_vol(
    market_price: float,
    S: float, K: float, r: float, q: float, T: float,
    option_type: str = "call",
    tol: float = 1e-8, max_iter: int = 100,
) -> float:
    """
    Newton-Raphson implied vol inversion with scipy.brentq fallback.

    Initial guess: Brenner-Subrahmanyam (1988) ATM approximation.
    """
    from scipy.optimize import brentq

    sigma = np.clip(np.sqrt(2 * np.pi / T) * market_price / S, 1e-4, 10.0)
    for _ in range(max_iter):
        res  = bs_european(S, K, r, q, sigma, T, option_type)
        diff = res.price - market_price
        if abs(diff) < tol:
            return sigma
        vega_raw = res.vega * 100.0
        if abs(vega_raw) < 1e-12:
            break
        sigma = np.clip(sigma - diff / vega_raw, 1e-6, 10.0)

    # Brent fallback
    f = lambda s: bs_european(S, K, r, q, s, T, option_type).price - market_price
    return brentq(f, 1e-6, 10.0, xtol=tol)

# visualization.py
"""
Five Plotly charts for the Black-Scholes Pricing Library.

All figures use plotly_dark template and the platform colour palette:
    Indigo  #6366f1  — primary series
    Emerald #10b981  — reference / benchmark
    Amber   #f59e0b  — secondary series
    Rose    #f43f5e  — put / warning
"""

from __future__ import annotations

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.stats   import norm as scipy_norm

from black_scholes import bs_european, bs_implied_vol
from monte_carlo   import mc_european_gbm, mc_convergence

INDIGO  = "#6366f1"
EMERALD = "#10b981"
AMBER   = "#f59e0b"
ROSE    = "#f43f5e"
SLATE   = "#94a3b8"

_LAYOUT = dict(
    template="plotly_dark",
    paper_bgcolor="rgba(0,0,0,0)",
    plot_bgcolor="rgba(0,0,0,0)",
    font=dict(color=SLATE, size=11),
    margin=dict(l=56, r=20, t=48, b=40),
)


# ── Chart 1: Price surface over (S, σ) ────────────────────────────────────────

def plot_price_surface(
    K: float = 100.0, r: float = 0.05, q: float = 0.0, T: float = 1.0,
    n_S: int = 40, n_sigma: int = 40,
    option_type: str = "call",
) -> go.Figure:
    """
    3-D surface of BS option price as a function of spot S and vol σ.

    Fixing K, r, q, T and sweeping S ∈ [60, 140] and σ ∈ [5%, 60%] shows:
    - How price increases monotonically with σ (vega > 0 always).
    - The curvature in S for a given σ (gamma effect).
    - The asymptotic approaches to intrinsic and zero at the boundary.
    """
    S_grid     = np.linspace(60, 140, n_S)
    sigma_grid = np.linspace(0.05, 0.60, n_sigma)
    Z = np.zeros((n_sigma, n_S))

    for i, sigma in enumerate(sigma_grid):
        for j, S in enumerate(S_grid):
            Z[i, j] = bs_european(S, K, r, q, sigma, T, option_type).price

    fig = go.Figure(data=go.Surface(
        x=S_grid, y=sigma_grid * 100, z=Z,
        colorscale="Plasma",
        colorbar=dict(title="Price", tickfont=dict(color="white")),
    ))
    fig.update_layout(
        **_LAYOUT,
        title=f"BS {option_type.capitalize()} Price — Spot × Volatility Surface",
        scene=dict(
            xaxis_title="Spot S",
            yaxis_title="Volatility σ (%)",
            zaxis_title="Option Price",
        ),
        height=500,
    )
    return fig


# ── Chart 2: Greeks profile vs spot ───────────────────────────────────────────

def plot_greeks_profile(
    K: float = 100.0, r: float = 0.05, q: float = 0.0,
    sigma: float = 0.20, T: float = 1.0,
    option_type: str = "call",
) -> go.Figure:
    """
    Delta, Gamma, and Vega as a function of spot S ∈ [50, 150].

    Panel 1 — Delta (Δ = ∂V/∂S):
        Call delta ∈ (0, e^{−qT}), increasing in S.
        Put  delta ∈ (−e^{−qT}, 0), increasing toward 0 as S rises.

    Panel 2 — Gamma (Γ = ∂²V/∂S²):
        Peaked ATM, symmetric in S for the log-normal model.
        Identical for calls and puts (put-call symmetry on second-order Greeks).

    Panel 3 — Vega (ν = ∂V/∂σ per 1 vol point):
        Peaked ATM; positive for all long options.
        At-the-money vega ≈ S·√T / (100·√(2π)).
    """
    S_grid = np.linspace(50, 150, 200)
    delta, gamma, vega = [], [], []

    for S in S_grid:
        res = bs_european(S, K, r, q, sigma, T, option_type)
        delta.append(res.delta)
        gamma.append(res.gamma)
        vega.append(res.vega)

    fig = make_subplots(
        rows=1, cols=3,
        subplot_titles=["Delta (Δ)", "Gamma (Γ)", "Vega (ν) per vol%"],
        horizontal_spacing=0.10,
    )
    color = INDIGO if option_type == "call" else ROSE

    fig.add_trace(go.Scatter(x=S_grid, y=delta, mode="lines",
        line=dict(color=color, width=2), name="Δ"), row=1, col=1)
    fig.add_trace(go.Scatter(x=S_grid, y=gamma, mode="lines",
        line=dict(color=EMERALD, width=2), name="Γ"), row=1, col=2)
    fig.add_trace(go.Scatter(x=S_grid, y=vega,  mode="lines",
        line=dict(color=AMBER, width=2), name="ν"), row=1, col=3)

    # ATM strike marker on each panel
    for col in range(1, 4):
        fig.add_vline(x=K, line=dict(color=SLATE, dash="dash", width=1),
                      row=1, col=col)

    fig.update_layout(**_LAYOUT, title="Black-Scholes Greeks Profile", height=320)
    fig.update_xaxes(title_text="Spot S")
    return fig


# ── Chart 3: Monte Carlo convergence ─────────────────────────────────────────

def plot_mc_convergence(
    S: float = 100.0, K: float = 100.0, r: float = 0.05,
    q: float = 0.0, sigma: float = 0.20, T: float = 1.0,
    n_max: int = 500_000, seed: int = 42,
) -> go.Figure:
    """
    MC call price and ±2σ confidence band as a function of path count (log scale).

    Shows O(1/√N) convergence. With antithetic variates, the effective variance is
    halved: σ²_anti = (σ²_plain + Cov[f(Z), f(−Z)]) / 2. For monotone payoffs the
    covariance term is negative, giving strictly faster convergence.

    Benchmark: BS closed-form.
    """
    pts    = mc_convergence(S, K, r, q, sigma, T, "call", n_max=n_max, seed=seed)
    ns     = np.array([p.n_paths   for p in pts])
    prices = np.array([p.price     for p in pts])
    stdes  = np.array([p.std_error for p in pts])

    analytic = bs_european(S, K, r, q, sigma, T, "call").price

    fig = go.Figure()

    # ±2σ confidence band
    fig.add_trace(go.Scatter(
        x=np.concatenate([ns, ns[::-1]]),
        y=np.concatenate([prices + 2 * stdes, (prices - 2 * stdes)[::-1]]),
        fill="toself", fillcolor="rgba(99,102,241,0.15)",
        line=dict(color="rgba(0,0,0,0)"), name="±2σ band",
    ))
    fig.add_trace(go.Scatter(
        x=ns, y=prices, mode="lines+markers",
        line=dict(color=INDIGO, width=2), marker=dict(size=4), name="MC price",
    ))
    fig.add_hline(
        y=analytic, line=dict(color=EMERALD, dash="dash", width=2),
        annotation_text=f"Analytic = {analytic:.4f}",
        annotation_font_color=EMERALD,
    )
    fig.update_layout(
        **_LAYOUT,
        title="Monte Carlo Convergence — BS Call Price (antithetic variates)",
        xaxis_title="Number of paths", xaxis_type="log",
        yaxis_title="Call Price",
        height=400,
    )
    return fig


# ── Chart 4: Crank-Nicolson grid convergence ──────────────────────────────────

def plot_cn_convergence(
    S: float = 100.0, K: float = 100.0, r: float = 0.05,
    q: float = 0.0, sigma: float = 0.20, T: float = 1.0,
    option_type: str = "call",
) -> go.Figure:
    """
    Absolute error of the Crank-Nicolson price vs the BS closed-form as a function
    of the number of grid points N (with n_S = n_τ = N).

    The Crank-Nicolson scheme is O(Δx², Δτ²) = O(1/N²) in both spatial and time
    steps. Doubling N reduces the error by a factor of ≈ 4 — visible as slope −2
    on the log-log plot.

    Limitation: near the payoff discontinuity (digital options), spurious
    Gibbs-like oscillations degrade convergence to O(1/N). Rannacher smoothing
    (two fully-implicit steps at τ=0) restores second-order accuracy.
    """
    from finite_difference import fd_european_cn

    grid_sizes = [25, 50, 100, 200, 400, 800]
    analytic   = bs_european(S, K, r, q, sigma, T, option_type).price
    errors     = []

    for n in grid_sizes:
        fd_price = fd_european_cn(S, K, r, q, sigma, T, option_type, n_S=n, n_t=n).price
        errors.append(abs(fd_price - analytic))

    # Reference: O(1/N²) slope anchored at the coarsest grid
    ref = [errors[0] * (grid_sizes[0] / n) ** 2 for n in grid_sizes]

    fig = go.Figure()
    fig.add_trace(go.Scatter(
        x=grid_sizes, y=errors, mode="lines+markers",
        line=dict(color=INDIGO, width=2), marker=dict(size=6), name="CN error",
    ))
    fig.add_trace(go.Scatter(
        x=grid_sizes, y=ref, mode="lines",
        line=dict(color=EMERALD, dash="dash", width=1.5),
        name="O(1/N²) reference",
    ))
    fig.update_layout(
        **_LAYOUT,
        title="Crank-Nicolson Grid Convergence — Absolute Error vs BS Analytic",
        xaxis_title="Grid size N  (n_S = n_τ = N)", xaxis_type="log",
        yaxis_title="Absolute error", yaxis_type="log",
        height=400,
    )
    return fig


# ── Chart 5: Discrete delta-hedging P&L distribution ─────────────────────────

def plot_hedging_pnl(
    S0: float = 100.0, K: float = 100.0, r: float = 0.05,
    q: float = 0.0, sigma: float = 0.20, T: float = 1.0,
    n_paths: int = 20_000, n_steps: int = 52,
    seed: int = 42,
) -> go.Figure:
    """
    P&L distribution of a dynamically delta-hedged call under discrete rebalancing.

    Under continuous hedging (Black-Scholes), the hedged portfolio is riskless.
    Under weekly rebalancing (n_steps=52), residual P&L arises from two effects:

        ΔP&L ≈ ½ Γ (ΔS)² − ½ σ² S² Γ Δt

    This is the gamma P&L: the hedger gains when realised variance exceeds the
    implied variance used to price the option. The distribution is approximately
    normal by CLT (many small gamma P&L contributions), centred near zero.

    The standard deviation of the terminal P&L scales as:
        std(P&L) ≈ (option_vega) × σ / √(2 n_steps)

    which quantifies the irreducible hedging noise from discrete rebalancing.
    """
    from scipy.stats import norm as scipy_norm

    rng  = np.random.default_rng(seed)
    dt   = T / n_steps
    disc = np.exp(-r * T)

    pnl = np.zeros(n_paths)

    for path in range(n_paths):
        S = S0
        t = 0.0
        # Short call + delta hedge
        call0     = bs_european(S, K, r, q, sigma, T, "call")
        portfolio = call0.price               # cash received for the call
        delta     = call0.delta
        bond      = portfolio - delta * S     # bond account (money market)

        for step in range(n_steps):
            tau = T - t
            dW  = rng.standard_normal() * np.sqrt(dt)
            S  *= np.exp((r - q - 0.5 * sigma**2) * dt + sigma * dW)
            t  += dt

            # Roll bond at risk-free rate
            bond  = bond * np.exp(r * dt)

            if step < n_steps - 1:
                tau_new   = T - t
                new_res   = bs_european(S, K, r, q, sigma, tau_new, "call")
                new_delta = new_res.delta
                # Rebalance: buy/sell (new_delta - delta) shares
                bond -= (new_delta - delta) * S
                delta = new_delta

        # At expiry: close hedge
        payoff = max(0.0, S - K)
        pnl[path] = delta * S + bond - payoff  # short call + delta → residual

    mu, std = np.mean(pnl), np.std(pnl)
    x_ref   = np.linspace(mu - 4 * std, mu + 4 * std, 300)

    fig = go.Figure()
    fig.add_trace(go.Histogram(
        x=pnl, histnorm="probability density", nbinsx=80,
        marker_color=INDIGO, opacity=0.75,
        name=f"P&L  (μ={mu:.3f}, σ={std:.3f})",
    ))
    fig.add_trace(go.Scatter(
        x=x_ref, y=scipy_norm.pdf(x_ref, mu, std), mode="lines",
        line=dict(color=EMERALD, width=2, dash="dash"),
        name="Normal reference",
    ))
    fig.add_vline(x=0, line=dict(color=SLATE, dash="dot", width=1))
    fig.update_layout(
        **_LAYOUT,
        title=f"Discrete Delta-Hedging P&L Distribution — {n_steps} rebalancing steps",
        xaxis_title="Terminal hedging P&L",
        yaxis_title="Probability density",
        height=400,
    )
    return fig


# ── Entry point ────────────────────────────────────────────────────────────────

if __name__ == "__main__":
    S, K, r, q, sigma, T = 100.0, 100.0, 0.05, 0.0, 0.20, 1.0
    res = bs_european(S, K, r, q, sigma, T, "call")
    print(f"Call price: {res.price:.4f}  Delta: {res.delta:.4f}  "
          f"Gamma: {res.gamma:.6f}  Vega: {res.vega:.4f}/vol%")

    plot_price_surface(K=K, r=r, q=q, T=T).show()
    plot_greeks_profile(K=K, r=r, q=q, sigma=sigma, T=T, option_type="call").show()
    plot_mc_convergence(S=S, K=K, r=r, q=q, sigma=sigma, T=T).show()
    plot_cn_convergence(S=S, K=K, r=r, q=q, sigma=sigma, T=T).show()
    plot_hedging_pnl(S0=S, K=K, r=r, q=q, sigma=sigma, T=T, n_paths=20_000, n_steps=52).show()