C++20 + Python● Live

Greeks and Risk Attribution Engine

Full analytic Greeks — first and second-order — for European options under Black-Scholes. Bump-and-reprice validation. Portfolio dollar Greeks with industry-standard conventions. P&L Taylor decomposition (delta/gamma/theta/vega/ vanna/volga/residual). Scenario grid via full reprice. Five Python Plotly visualisations.

C++20Analytic GreeksBump-revalDollar GreeksP&L AttributionScenario AnalysisPlotly

Interactive Explorer

Greeks & Risk Attribution Engine

Black-Scholes analytic Greeks · P&L decomposition · scenario analysis

Option type

Spot (S)$100

Strike (K)$100

Vol (σ)20%

Expiry (T)1.00y

Rate (r)5.00%

Div yield (q)2.0%

Price: $9.2270

Δ 0.5869Γ 0.0190ν 0.3790θ 0.0139Λ -0.0948Ω 0.0002

Delta∂V/∂S

0.5869

Delta

0.5869

Gamma

1.8951

Vega

0.3790

Theory and Derivation

1. Model Assumptions

All Greeks in this engine are derived from the Black-Scholes model. The underlying follows geometric Brownian motion under the risk-neutral measure $\mathbb{Q}$ :

$dS = (r - q)\, S\, dt + \sigma\, S\, dW^Q$

with constant risk-free rate $r$ , continuous dividend yield $q$ , and volatility $\sigma$ . European exercise; no transaction costs. All rates are continuously compounded and annualised.

The auxiliary quantities used throughout:

$d_1 = \frac{\ln(S/K) + (r - q + \frac{1}{2}\sigma^2)\,T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T}$

2. First-Order Greeks

The standard Greek set for a European call ( $C$ ) or put ( $P$ ):

$\Delta_C = e^{-qT} N(d_1), \qquad \Delta_P = -e^{-qT} N(-d_1)$

$\Gamma = \frac{e^{-qT}\,\phi(d_1)}{S\,\sigma\,\sqrt{T}} \quad\text{(identical for call and put)}$

$\mathcal{V} = S\,e^{-qT}\,\phi(d_1)\,\sqrt{T} \quad\text{per unit vol (identical for call and put)}$

$\Theta_C = -\frac{S\,e^{-qT}\,\phi(d_1)\,\sigma}{2\sqrt{T}} - r\,K\,e^{-rT}\,N(d_2) + q\,S\,e^{-qT}\,N(d_1)$

$\rho_C = K\,T\,e^{-rT}\,N(d_2)$

Put-call parity implies $\Delta_C - \Delta_P = e^{-qT}$ , $\Gamma_C = \Gamma_P$ , and $\mathcal{V}_C = \mathcal{V}_P$ . These identities provide exact validation benchmarks.

Every Greek satisfies the Black-Scholes PDE as a consequence of the pricing formula:

$\Delta\,(r-q)\,S + \tfrac{1}{2}\,\Gamma\,\sigma^2\,S^2 + \Theta_{\text{annual}} = r\,V$

This identity is enforced as a unit test, verifying internal consistency of the Greek implementations.

3. Second-Order and Cross Greeks

The first step is the key partial derivative $\partial d_1/\partial\sigma$ . Differentiating $d_1 = [\ln(S/K) + (r-q+\sigma^2/2)T]/(\sigma\sqrt{T})$ by the quotient rule:

$\frac{\partial d_1}{\partial \sigma} = -\frac{d_2}{\sigma}, \qquad \frac{\partial d_2}{\partial \sigma} = -\frac{d_1}{\sigma}$

Vanna — $\partial^2 V/\partial S\,\partial\sigma = \partial\Delta/\partial\sigma = \partial\mathcal{V}/\partial S$ :

$\Lambda = \frac{\partial \Delta}{\partial \sigma} = e^{-qT}\,\phi(d_1)\cdot\frac{\partial d_1}{\partial \sigma} = -\frac{e^{-qT}\,\phi(d_1)\,d_2}{\sigma}$

Vanna is the same for calls and puts — confirmed by Schwarz's theorem: $\partial\Delta/\partial\sigma = \partial^2 V/\partial S\,\partial\sigma = \partial\mathcal{V}/\partial S$ . It is zero at ATM ( $d_2 = 0$ ), positive for OTM calls ( $d_2 < 0$ ), and negative for ITM calls.

Volga / Vomma — $\partial^2 V/\partial\sigma^2$ :

$\Omega = \frac{\partial \mathcal{V}}{\partial \sigma} = \mathcal{V}\cdot\frac{d_1\,d_2}{\sigma}$

Since $d_1 > 0 > d_2$ at ATM (for positive rates and short-to-medium tenors), volga is negative at ATM — the BS price function is locally concave in vol near the money. It turns positive for far OTM/ITM options where d₁ and d₂ share the same sign. The vanna-volga approximation exploits this structure to construct smile-consistent prices.

Charm — $\partial\Delta/\partial t$ , delta bleed per calendar day:

$\text{Charm}_C = \frac{1}{365}\left[q\,e^{-qT}\,N(d_1) - \frac{e^{-qT}\,\phi(d_1)\,\bigl[2(r-q)T - d_2\,\sigma\sqrt{T}\bigr]}{2T\,\sigma\sqrt{T}}\right]$

Charm captures delta's daily time decay. A long call has negative charm near expiry as delta collapses from 0.5 toward 0 or 1. Dealers hedge dynamically by monitoring the combined theta and charm exposure.

Speed ( $\partial\Gamma/\partial S$ ), Zomma ( $\partial\Gamma/\partial\sigma$ ), Color ( $\partial\Gamma/\partial t$ ), Veta ( $\partial\mathcal{V}/\partial t$ ):

$\text{Speed} = -\frac{\Gamma}{S}\left(1 + \frac{d_1}{\sigma\sqrt{T}}\right)$

$\text{Zomma} = \frac{\Gamma\,(d_1\,d_2 - 1)}{\sigma}$

Speed and zomma are the sensitivity of gamma to spot and vol respectively — critical for gamma scalping: as gamma itself moves with the market, a gamma-neutral book is only momentarily neutral.

4. Bump-and-Reprice: Finite Difference Validation

Analytic Greeks are validated against central finite differences. For a bump size $h$ , the standard stencils are:

$\Delta \approx \frac{V(S+h) - V(S-h)}{2h}, \qquad \Gamma \approx \frac{V(S+h) - 2V(S) + V(S-h)}{h^2}$

$\Lambda \approx \frac{V(S+h_S, \sigma+h_\sigma) - V(S+h_S,\sigma-h_\sigma) - V(S-h_S,\sigma+h_\sigma) + V(S-h_S,\sigma-h_\sigma)}{4\,h_S\,h_\sigma}$

All central stencils have O(h²) truncation error. Default bumps are chosen to balance truncation error against floating-point cancellation: 1% of spot for delta/gamma, 0.1 vol-point for vega/vanna/volga, 1 bp for rho, 1 calendar hour for theta. In production systems, bump sizes are calibrated per instrument type and risk factor; a universal bump size is a common source of Greek hedging error.

5. Dollar Greeks and Portfolio Aggregation

Raw Black-Scholes Greeks are dimensionally inconsistent across positions (delta is dimensionless; gamma is per dollar of spot; vega is per unit vol). Desks convert to dollar Greeks for portfolio-level aggregation:

$\$\Delta = \Delta \cdot S \cdot \text{qty} \cdot \text{mult}$

$\$\Gamma = \Gamma \cdot (0.01\,S)^2 \cdot \text{qty} \cdot \text{mult}$

$\$\mathcal{V} = \mathcal{V}_{\text{unit}} \cdot 0.01 \cdot \text{qty} \cdot \text{mult}$

$\$\Delta$ is the dollar P&L per $1 absolute spot move. $\$\Gamma$ is the dollar P&L per 1% spot move (from the quadratic term). $\$\mathcal{V}$ is the dollar P&L per 1 vol-point. A hedged book has near-zero aggregate dollar Greeks; residual exposure after hedging is the unhedged risk.

6. P&L Attribution: Taylor Decomposition

For a position with start-of-day market state $(S_0, \sigma_0, T_0)$ and end-of-day state $(S_1, \sigma_1, T_1 = T_0 - \Delta t)$ , the Taylor expansion of the P&L gives:

$\Delta V \approx \underbrace{\Delta\cdot\Delta S}_{\text{delta}} + \underbrace{\tfrac{1}{2}\Gamma\cdot(\Delta S)^2}_{\text{gamma}} + \underbrace{\theta\cdot\Delta t}_{\text{theta}} + \underbrace{\mathcal{V}\cdot\Delta\sigma}_{\text{vega}} + \underbrace{\Lambda\cdot\Delta S\cdot\Delta\sigma}_{\text{vanna}} + \underbrace{\tfrac{1}{2}\Omega\cdot(\Delta\sigma)^2}_{\text{volga}} + \varepsilon$

The residual $\varepsilon = \Delta V_{\text{actual}} - \Delta V_{\text{explained}}$ captures higher-order Taylor terms ( $\mathcal{O}(\Delta S^3,\,\Delta\sigma^3)$ ), model error (departure from GBM), and discretisation error in snapshot timing. For a delta-neutral book with overnight moves of 1–2%, the residual should be under 1% of total explained. Persistent large residuals signal model breakdown or data quality issues.

On structured desks, the P&L explain is run as a formal daily process. Risk management and finance require agreement between risk-sensitivities-based explain and accounting P&L to within a defined tolerance. Unexplained P&L is an audit trigger.

7. Model Limitations

—Constant volatility. All Greeks assume σ is a fixed model parameter. In reality, implied vol moves with spot and time. Vanna captures the joint ΔS−Δσ exposure, but the vol dynamics model is needed to determine Δσ itself.
—Near-expiry singularities. Gamma, charm, color, speed, zomma, veta all diverge as T → 0 or σ → 0. For near-expiry options, numerical stability requires minimum T of ~1 hour.
—GBM dynamics. No jumps, no stochastic vol. Jump risk (gap risk) is completely unmodelled. For equities, the vol skew embeds expected jump size — BS Greeks are wrong in a jump world.
—Single-factor assumption. Scenario grid applies the same shock to all positions' spots. Cross-asset or multi-leg books require a correlation model for joint scenarios.
—Discretisation in P&L explain. Taylor decomposition uses SOD Greeks as coefficients. For large intraday moves, the approximation degrades. Production systems recompute Greeks at intermediate states or use higher-order stencils.
—American options. This engine covers European exercise only. For American options, use the LSMC or finite-difference American pricer; analytic Greeks are not available in closed form.

Interview Angle

Junior (L1). Define delta, gamma, vega, theta, rho. Explain the sign conventions. What does it mean for gamma to be positive? Why is theta negative for long options? Reproduce the Black-Scholes PDE from the Greeks. What is put-call parity for delta?

Mid-level (L2). Derive vanna from scratch (quotient rule on d₁, then chain rule). Why is volga negative ATM and positive far OTM? Explain the P&L attribution formula — what does each term represent, and what does a large residual signal? How do you size the bump for a finite-difference Greek to minimise error?

Senior (L3). Describe the full daily P&L explain process on a structured desk. Why does a BS delta hedge fail to explain P&L in a stochastic vol world? How does a vanna-volga pricer use Λ and Ω to add a smile-consistent correction to a flat-vol price? What is the connection between charm and the term structure of delta? When would color matter more than gamma for a risk manager?

C++20 Implementation

Production-grade. C++20, CMake, Catch2. Strict warnings. No UB.

Greeks.hpp — interface and type definitions

// Greeks.hpp — analytic and finite-difference Greeks for European options.
// Model: dS = (r−q)S dt + σS dW  under Q.
// Conventions: r, q continuously compounded annualised decimal; σ annualised decimal;
//              T in years; theta/charm/color/veta per calendar day.

#pragma once
#include "common.hpp"
#include <vector>

namespace bb {

// -------------------------------------------------------------------------
// Full analytic Greek set — first and second order
// -------------------------------------------------------------------------
struct GreeksResult {
    double price;
    // First order
    double delta;   ///< ∂V/∂S
    double gamma;   ///< ∂²V/∂S²
    double vega;    ///< ∂V/∂σ  (per unit vol; ÷100 for per vol-point)
    double theta;   ///< ∂V/∂t  per calendar day (negative for long options)
    double rho;     ///< ∂V/∂r  (per unit rate; ÷100 for per bp)
    // Second order / cross
    double vanna;   ///< ∂²V/∂S∂σ = ∂delta/∂σ = ∂vega/∂S
    double volga;   ///< ∂²V/∂σ²  (Vomma)
    double charm;   ///< ∂²V/∂S∂t per calendar day (delta bleed)
    double speed;   ///< ∂³V/∂S³  = ∂gamma/∂S
    double zomma;   ///< ∂³V/∂S²∂σ = ∂gamma/∂σ
    double color;   ///< ∂³V/∂S²∂t per calendar day (gamma bleed)
    double veta;    ///< ∂²V/∂σ∂t per calendar day (vega bleed)
};

GreeksResult bs_greeks_full(double S, double K, double r, double q,
                             double sigma, double T, OptionType type);

// -------------------------------------------------------------------------
// Bump-and-reprice finite-difference Greeks (validation)
// -------------------------------------------------------------------------
struct BumpGreeksResult {
    double delta, gamma, vega, theta, rho, vanna, volga;
};

BumpGreeksResult bs_greeks_bump(double S, double K, double r, double q,
                                 double sigma, double T, OptionType type,
                                 double h_spot  = 0.0,    // 0 = auto (1% of S)
                                 double h_vol   = 0.001,
                                 double h_rate  = 0.0001,
                                 double h_time  = 1.0/365.0);

// -------------------------------------------------------------------------
// Portfolio dollar Greeks
// -------------------------------------------------------------------------
struct PortfolioPosition {
    double S, K, r, q, sigma, T;
    OptionType type;
    double quantity;    // signed (long > 0, short < 0)
    double multiplier;  // shares per contract (100 for equity options)
};

struct DollarGreeks {
    double dollar_delta;   // $ per $1 spot move
    double dollar_gamma;   // $ per 1% spot move²
    double dollar_vega;    // $ per 1 vol-point
    double dollar_theta;   // $ per calendar day
    double dollar_vanna;   // $ per 1% spot × 1 vol-pt
    double dollar_volga;   // $ per (vol-pt)²
    double total_value;
};

DollarGreeks portfolio_dollar_greeks(const std::vector<PortfolioPosition>& positions);

// -------------------------------------------------------------------------
// P&L attribution
// -------------------------------------------------------------------------
struct PnlAttribution {
    double delta_pnl, gamma_pnl, theta_pnl, vega_pnl;
    double vanna_pnl, volga_pnl;
    double total_explained, actual_pnl, unexplained;
};

PnlAttribution pnl_explain(double S_old, double S_new,
                             double sigma_old, double sigma_new,
                             double K, double r, double q,
                             double T_old, OptionType type,
                             double quantity, double multiplier,
                             double delta_t);

// -------------------------------------------------------------------------
// Scenario analysis
// -------------------------------------------------------------------------
struct ScenarioPoint {
    double spot_shock_pct, vol_shock_abs, pnl, portfolio_value;
};

std::vector<ScenarioPoint> scenario_grid(
    const std::vector<PortfolioPosition>& positions,
    const std::vector<double>& spot_shocks,
    const std::vector<double>& vol_shocks);

} // namespace bb

Greeks.cpp — analytic Greek derivations (excerpts)

// Greeks.cpp — selected implementation excerpts

// Core quantities for any Black-Scholes Greek.
// d1 = [ln(S/K) + (r − q + σ²/2)T] / (σ√T)
// d2 = d1 − σ√T
struct D1D2 { double d1, d2; };

static inline D1D2 compute_d1d2(double S, double K, double r, double q,
                                  double sigma, double T) noexcept {
    const double sqrtT = std::sqrt(T);
    const double d1 = (std::log(S/K) + (r - q + 0.5*sigma*sigma)*T) / (sigma*sqrtT);
    return {d1, d1 - sigma*sqrtT};
}

GreeksResult bs_greeks_full(double S, double K, double r, double q,
                             double sigma, double T, OptionType type) {
    const auto [d1, d2] = compute_d1d2(S, K, r, q, sigma, T);
    const double sqrtT  = std::sqrt(T);
    const double phi    = normpdf(d1);
    const double eqT    = std::exp(-q*T), erT = std::exp(-r*T);
    const double Nd1    = normcdf(d1),  Nd2  = normcdf(d2);
    const double Nm1    = normcdf(-d1), Nm2  = normcdf(-d2);
    const bool   call   = (type == OptionType::Call);

    // -----------------------------------------------------------------------
    // First-order Greeks
    // -----------------------------------------------------------------------
    const double price = call ? S*eqT*Nd1 - K*erT*Nd2
                              : K*erT*Nm2 - S*eqT*Nm1;
    const double delta = call ?  eqT*Nd1 : -eqT*Nm1;
    const double gamma = eqT*phi / (S*sigma*sqrtT);     // same call/put
    const double vega  = S*eqT*phi*sqrtT;               // same call/put, per unit vol

    // Θ:  ∂V/∂t per calendar day.  We compute ∂V/∂T then negate and scale.
    // ∂V_call/∂T = −(S·eqT·φ·σ)/(2√T) − r·K·e^{−rT}·N(d2) + q·S·eqT·N(d1)
    const double dVdT = -(S*eqT*phi*sigma)/(2*sqrtT)
                        + (call ? -r*K*erT*Nd2 + q*S*eqT*Nd1
                                : +r*K*erT*Nm2 - q*S*eqT*Nm1);
    const double theta = -dVdT / 365.0;

    const double rho = call ?  K*T*erT*Nd2 : -K*T*erT*Nm2;

    // -----------------------------------------------------------------------
    // Second-order / cross Greeks
    // -----------------------------------------------------------------------

    // Vanna: ∂delta/∂σ.
    //   Δ_call = eqT·N(d1);  ∂d1/∂σ = −d2/σ  (derived via quotient rule on d1).
    //   ∂Δ/∂σ = eqT·φ(d1)·(−d2/σ) = −eqT·φ(d1)·d2/σ  — same for calls and puts.
    const double vanna = -eqT*phi*d2/sigma;

    // Volga (Vomma): ∂vega/∂σ.
    //   vega = S·eqT·φ(d1)·√T;  ∂/∂σ = S·eqT·√T·(−d1·φ)·(−d2/σ) = vega·d1·d2/σ.
    //   Note: volga < 0 ATM (where d1·d2 < 0); > 0 for far OTM/ITM.
    const double volga = vega*d1*d2/sigma;

    // Charm: ∂delta/∂t per calendar day.  Convention: positive charm = delta rises.
    //   Standard Haug form for ∂Δ/∂T then negated and scaled:
    const double charm_common = eqT*phi*(2*(r-q)*T - d2*sigma*sqrtT)
                                / (2*T*sigma*sqrtT);
    const double charm = (call ? q*eqT*Nd1 - charm_common
                               : -q*eqT*Nm1 - charm_common) / 365.0;

    // Speed: ∂gamma/∂S.
    //   Γ = eqT·φ/(S·σ·√T);  ∂Γ/∂S = −(Γ/S)·(1 + d1/(σ√T)).
    const double speed = -gamma*(1 + d1/(sigma*sqrtT))/S;

    // Zomma: ∂gamma/∂σ.
    //   ∂Γ/∂σ = Γ·(d1·d2 − 1)/σ.
    const double zomma = gamma*(d1*d2 - 1)/sigma;

    // Color: ∂gamma/∂t per calendar day.
    //   ∂Γ/∂T = −Γ·[2qT + 1 + d1·(2(r−q)T − d2·σ·√T)/(σ·√T)] / (2T)
    const double dGdT = gamma*(2*q*T + 1 + d1*(2*(r-q)*T - d2*sigma*sqrtT)
                                           /(sigma*sqrtT)) / (-2*T);
    const double color = -dGdT / 365.0;

    // Veta: ∂vega/∂t per calendar day.
    //   ∂V/∂T = vega·[−q − (1 + d1²)/(2T) + d1·(r−q+σ²/2)/(σ√T)]
    const double dVegadT = vega*(-q - (1+d1*d1)/(2*T)
                                    + d1*(r-q+0.5*sigma*sigma)/(sigma*sqrtT));
    const double veta = -dVegadT / 365.0;

    return {price, delta, gamma, vega, theta, rho,
            vanna, volga, charm, speed, zomma, color, veta};
}

Greeks.cpp — P&L attribution

// P&L attribution — Taylor decomposition vs. actual full reprice.
//
//   ΔV ≈ Δ·ΔS + ½Γ·(ΔS)² + θ·Δt + ν·Δσ + Λ·ΔS·Δσ + ½Ω·(Δσ)² + residual
//
// All components are in dollar terms (position size × multiplier included).

PnlAttribution pnl_explain(double S_old, double S_new,
                             double sigma_old, double sigma_new,
                             double K, double r, double q,
                             double T_old, OptionType type,
                             double quantity, double multiplier,
                             double delta_t) {
    const double T_new = T_old - delta_t;
    if (T_new <= 0.0) throw std::invalid_argument("pnl_explain: option expires mid-period");

    // Greeks at SOD state — Taylor coefficients.
    const GreeksResult g = bs_greeks_full(S_old, K, r, q, sigma_old, T_old, type);
    const double scale = quantity * multiplier;
    const double dS    = S_new  - S_old;
    const double dSig  = sigma_new - sigma_old;
    const double dt_days = delta_t * 365.0;

    const double delta_pnl = g.delta * dS          * scale;
    const double gamma_pnl = 0.5 * g.gamma * dS*dS * scale;
    const double theta_pnl = g.theta * dt_days      * scale;
    const double vega_pnl  = g.vega  * dSig         * scale;
    const double vanna_pnl = g.vanna * dS * dSig    * scale;
    const double volga_pnl = 0.5 * g.volga * dSig*dSig * scale;

    const double explained = delta_pnl + gamma_pnl + theta_pnl
                           + vega_pnl  + vanna_pnl + volga_pnl;

    // Full mark-to-market change.
    const double V_old   = g.price * scale;
    const double V_new   = bs_european(S_new, K, r, q, sigma_new, T_new, type).price * scale;
    const double actual  = V_new - V_old;

    return {delta_pnl, gamma_pnl, theta_pnl, vega_pnl,
            vanna_pnl, volga_pnl, explained, actual, actual - explained};
}

Greeks.cpp — scenario grid

// Scenario grid — full reprice over (spot shock, vol shock) pairs.
// Spot is shocked relative (ss=0.05 → S×1.05); vol absolute (vs=0.05 → σ+0.05).
// Vol is clamped to [1e-4, 5.0] to prevent pricing failures under extreme shocks.

std::vector<ScenarioPoint> scenario_grid(
    const std::vector<PortfolioPosition>& positions,
    const std::vector<double>& spot_shocks,
    const std::vector<double>& vol_shocks) {

    // Base portfolio value (zero shocks).
    double base = 0.0;
    for (const auto& pos : positions)
        base += bs_european(pos.S, pos.K, pos.r, pos.q,
                             pos.sigma, pos.T, pos.type).price
                * pos.quantity * pos.multiplier;

    std::vector<ScenarioPoint> grid;
    grid.reserve(spot_shocks.size() * vol_shocks.size());

    for (double ss : spot_shocks) {
        for (double vs : vol_shocks) {
            double shocked = 0.0;
            for (const auto& pos : positions) {
                const double S_sh  = pos.S * (1.0 + ss);
                const double sv_sh = std::clamp(pos.sigma + vs, 1e-4, 5.0);
                shocked += bs_european(S_sh, pos.K, pos.r, pos.q,
                                        sv_sh, pos.T, pos.type).price
                           * pos.quantity * pos.multiplier;
            }
            grid.push_back({ss, vs, shocked - base, shocked});
        }
    }
    return grid;
}

test_greeks.cpp — selected Catch2 tests (16 total)

// test_greeks.cpp — selected Catch2 v3 tests.

// 1. Put-call parity for delta: Δ_call − Δ_put = e^{−qT}
TEST_CASE("Greeks: put-call parity for delta", "[greeks][pcp]") {
    auto gc = bs_greeks_full(100, 100, 0.05, 0.02, 0.20, 1.0, OptionType::Call);
    auto gp = bs_greeks_full(100, 100, 0.05, 0.02, 0.20, 1.0, OptionType::Put);
    REQUIRE_THAT(gc.delta - gp.delta, WithinAbs(std::exp(-0.02*1.0), 1e-12));
}

// 2. Black-Scholes PDE: Δ·(r−q)·S + ½·Γ·σ²·S² + θ_annual = r·V
TEST_CASE("Greeks: satisfy Black-Scholes PDE", "[greeks][pde]") {
    for (double K : {85.0, 100.0, 115.0}) {
        auto g = bs_greeks_full(100, K, 0.05, 0.02, 0.20, 1.0, OptionType::Call);
        double lhs = g.delta*(0.05-0.02)*100 + 0.5*g.gamma*0.04*10000 + g.theta*365;
        REQUIRE_THAT(lhs, WithinAbs(0.05*g.price, 1e-8));
    }
}

// 3. Vanna symmetry: ∂delta/∂σ = ∂vega/∂S (Schwarz's theorem)
TEST_CASE("Greeks: vanna symmetry", "[greeks][symmetry]") {
    auto g = bs_greeks_full(100, 100, 0.05, 0.02, 0.20, 1.0, OptionType::Call);
    const double h = 1e-4;
    double dDelta_dSig = (bs_greeks_full(100,100,0.05,0.02,0.20+h,1.0,OptionType::Call).delta
                        - bs_greeks_full(100,100,0.05,0.02,0.20-h,1.0,OptionType::Call).delta)/(2*h);
    double dVega_dS    = (bs_greeks_full(100+0.01,100,0.05,0.02,0.20,1.0,OptionType::Call).vega
                        - bs_greeks_full(100-0.01,100,0.05,0.02,0.20,1.0,OptionType::Call).vega)/0.02;
    REQUIRE_THAT(g.vanna, WithinAbs(dDelta_dSig, 5e-5));
    REQUIRE_THAT(g.vanna, WithinAbs(dVega_dS,    5e-5));
}

// 4. Analytic vs. bump delta (central FD, h=0.01)
TEST_CASE("Greeks: analytic vs. bump delta", "[greeks][bump]") {
    for (double K : {85.0, 100.0, 115.0}) {
        auto ga = bs_greeks_full(100, K, 0.05, 0.02, 0.20, 1.0, OptionType::Call);
        auto gb = bs_greeks_bump(100, K, 0.05, 0.02, 0.20, 1.0, OptionType::Call, 0.01);
        REQUIRE_THAT(ga.delta, WithinAbs(gb.delta, 1e-6));
    }
}

// 5. P&L attribution: residual < 1% for small overnight move
TEST_CASE("Greeks: P&L attribution accuracy", "[greeks][pnl]") {
    auto attr = pnl_explain(100, 101, 0.20, 0.205, 100, 0.05, 0.02,
                             1.0, OptionType::Call, 1.0, 100.0, 1.0/365.0);
    double tol = std::abs(attr.actual_pnl) * 0.01;
    REQUIRE_THAT(attr.unexplained, WithinAbs(0.0, tol + 1e-6));
}

// 6. Scenario grid: zero shock = zero P&L
TEST_CASE("Greeks: scenario zero-shock", "[greeks][scenario]") {
    std::vector<PortfolioPosition> pos = {{100,100,0.05,0.02,0.20,1.0,OptionType::Call,10,100}};
    auto grid = scenario_grid(pos, {-0.05, 0.0, 0.05}, {-0.05, 0.0, 0.05});
    REQUIRE_THAT(grid[4].pnl, WithinAbs(0.0, 1e-10));  // centre = (0,0)
}

Python Reporting — Plotly Visualisations

Five production-style charts. Run end-to-end in a clean Python 3.10+ environment. Dependency: scipy, numpy, plotly. No other quant libraries required.

1. Greeks Surface — Delta, Gamma, Vega over Strike × Expiry

3-D surface plots of the three primary Greeks as a function of strike and time to expiry. Delta surface is smooth; Gamma surface peaks sharply ATM at short tenors (pin risk). Vega surface peaks ATM at intermediate tenors.

"""
Greeks surface — Delta, Gamma, Vega as functions of moneyness and time to expiry.
Plotly interactive 3-D surface; no external quant library required.
"""
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.stats import norm

def bs_greeks(S, K, r, q, sigma, T):
    sqrtT = np.sqrt(T)
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*sqrtT)
    d2 = d1 - sigma*sqrtT
    phi = norm.pdf(d1)
    eqT = np.exp(-q*T)

    delta = eqT * norm.cdf(d1)
    gamma = eqT * phi / (S * sigma * sqrtT)
    vega  = S * eqT * phi * sqrtT  # per unit vol
    return delta, gamma, vega

# Grid: moneyness (S/K) × time to expiry
S, r, q, sigma = 100.0, 0.05, 0.02, 0.20
strikes = np.linspace(70, 130, 60)
tenors  = np.linspace(0.05, 2.0, 40)

K_grid, T_grid  = np.meshgrid(strikes, tenors)
m_grid = S / K_grid  # moneyness = S/K

delta_surf, gamma_surf, vega_surf = bs_greeks(S, K_grid, r, q, sigma, T_grid)

fig = make_subplots(
    rows=1, cols=3,
    subplot_titles=["Delta (Δ)", "Gamma × 100 (Γ)", "Vega per vol-pt (ν)"],
    specs=[[{"type": "surface"}] * 3],
)

for col, (surf, title) in enumerate(
    [(delta_surf, "Delta"), (gamma_surf * 100, "Gamma×100"),
     (vega_surf * 0.01, "Vega/pp")], start=1
):
    fig.add_trace(
        go.Surface(x=strikes, y=tenors, z=surf,
                   colorscale="RdBu", showscale=(col == 3),
                   name=title),
        row=1, col=col,
    )

fig.update_layout(
    title="Black-Scholes Greeks: Strike × Expiry surface",
    height=500,
    scene=dict(xaxis_title="Strike K", yaxis_title="Expiry T (yr)", zaxis_title="Δ"),
    scene2=dict(xaxis_title="Strike K", yaxis_title="Expiry T (yr)", zaxis_title="Γ×100"),
    scene3=dict(xaxis_title="Strike K", yaxis_title="Expiry T (yr)", zaxis_title="ν/pp"),
)
fig.show()

2. Vanna-Volga Surface — Cross-Greek Exposure

Vanna and Volga as joint functions of strike and implied vol. Vanna is antisymmetric around ATM. Volga has a saddle structure: negative ATM, positive at the wings. These two Greeks are the core inputs to the vanna-volga smile construction method.

"""
Vanna-Volga surface — second-order cross Greeks across (moneyness, vol) space.
Desks use vanna and volga to quantify exposure to joint spot-vol and vol² moves.
"""
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.stats import norm

def vanna_volga(S, K, r, q, sigma, T):
    sqrtT = np.sqrt(T)
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*sqrtT)
    d2 = d1 - sigma*sqrtT
    phi = norm.pdf(d1)
    eqT = np.exp(-q*T)
    vega = S * eqT * phi * sqrtT

    vanna = -eqT * phi * d2 / sigma         # ∂²V/∂S∂σ
    volga = vega * d1 * d2 / sigma          # ∂²V/∂σ²
    return vanna, volga

S, r, q, T = 100.0, 0.05, 0.02, 1.0
strikes = np.linspace(60, 140, 80)
vols    = np.linspace(0.10, 0.60, 60)

K_grid, sig_grid = np.meshgrid(strikes, vols)
vanna_surf, volga_surf = vanna_volga(S, K_grid, r, q, sig_grid, T)

fig = make_subplots(
    rows=1, cols=2,
    subplot_titles=["Vanna (Λ) = ∂²V/∂S∂σ", "Volga (Ω) = ∂²V/∂σ²"],
    specs=[[{"type": "surface"}, {"type": "surface"}]],
)
fig.add_trace(go.Surface(x=strikes, y=vols, z=vanna_surf,
                          colorscale="RdBu", showscale=True), row=1, col=1)
fig.add_trace(go.Surface(x=strikes, y=vols, z=volga_surf,
                          colorscale="Viridis", showscale=True), row=1, col=2)

fig.update_layout(
    title=f"Vanna-Volga surface (S={S}, r={r:.0%}, q={q:.0%}, T={T}y)",
    height=500,
    scene=dict( xaxis_title="Strike K", yaxis_title="Implied Vol σ", zaxis_title="Vanna"),
    scene2=dict(xaxis_title="Strike K", yaxis_title="Implied Vol σ", zaxis_title="Volga"),
)
fig.show()

3. P&L Attribution Waterfall

Annotated bar chart decomposing overnight P&L into delta, gamma, theta, vega, vanna, volga, and unexplained residual. Run for any market move scenario. The residual should be near zero for small moves; large residuals signal model error or data issues.

"""
P&L attribution waterfall — Taylor decomposition vs. actual full reprice.
Annotated bar chart showing delta / gamma / theta / vega / vanna / volga / residual.
"""
import numpy as np
import plotly.graph_objects as go
from scipy.stats import norm

def bs_greeks_full(S, K, r, q, sigma, T, call=True):
    sqrtT = np.sqrt(T)
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*sqrtT)
    d2 = d1 - sigma*sqrtT
    phi = norm.pdf(d1)
    eqT, erT = np.exp(-q*T), np.exp(-r*T)
    Nd1, Nd2 = norm.cdf(d1), norm.cdf(d2)

    price = (S*eqT*Nd1 - K*erT*Nd2) if call else (K*erT*norm.cdf(-d2) - S*eqT*norm.cdf(-d1))
    delta = eqT*Nd1 if call else -eqT*norm.cdf(-d1)
    gamma = eqT*phi / (S*sigma*sqrtT)
    vega  = S*eqT*phi*sqrtT
    dVdT  = -(S*eqT*phi*sigma)/(2*sqrtT) + (-r*K*erT*Nd2 + q*S*eqT*Nd1 if call
                                              else r*K*erT*norm.cdf(-d2) - q*S*eqT*norm.cdf(-d1))
    theta = -dVdT / 365
    vanna = -eqT*phi*d2/sigma
    volga = vega*d1*d2/sigma
    return dict(price=price, delta=delta, gamma=gamma, vega=vega,
                theta=theta, vanna=vanna, volga=volga)

# Parameters
S0, K, r, q, sigma0, T0, qty, mult = 100.0, 100.0, 0.05, 0.02, 0.20, 1.0, 10, 100

# Market move
dS_pct, dSig_pp, dt_days = 2.0, 1.5, 1.0
S1     = S0 * (1 + dS_pct/100)
sigma1 = sigma0 + dSig_pp/100
T1     = T0 - dt_days/365

# Compute
g   = bs_greeks_full(S0, K, r, q, sigma0, T0)
g1  = bs_greeks_full(S1, K, r, q, sigma1, T1)
scl = qty * mult
dS, dSig = S1 - S0, sigma1 - sigma0

components = {
    "Delta":      g["delta"] * dS          * scl,
    "Gamma":      0.5 * g["gamma"] * dS**2 * scl,
    "Theta":      g["theta"] * dt_days      * scl,
    "Vega":       g["vega"]  * dSig         * scl,
    "Vanna":      g["vanna"] * dS * dSig   * scl,
    "Volga":      0.5 * g["volga"] * dSig**2 * scl,
}
actual     = (g1["price"] - g["price"]) * scl
explained  = sum(components.values())
unexplained = actual - explained
components["Unexplained"] = unexplained

labels = list(components.keys())
values = list(components.values())
colors = ["#3b82f6","#10b981","#ef4444","#f59e0b","#8b5cf6","#ec4899","#94a3b8"]

fig = go.Figure(go.Bar(
    x=labels, y=values,
    marker_color=[c if v >= 0 else c for c, v in zip(colors, values)],
    marker_line_color="white",
    marker_line_width=1,
    text=[f"${v:+.2f}" for v in values],
    textposition="outside",
))
fig.add_hline(y=actual, line_dash="dash", line_color="black",
              annotation_text=f"Actual PnL: ${actual:+.2f}")
fig.update_layout(
    title=f"P&L Attribution — ΔS={dS_pct:+.1f}%, Δσ={dSig_pp:+.1f}pp, Δt={dt_days}d",
    yaxis_title="P&L ($)",
    xaxis_title="Attribution Component",
    height=420,
)
fig.show()

4. Scenario Heatmap — Full-Reprice Grid

Colour-coded grid of portfolio P&L over a 9×9 spot-shock × vol-shock grid. Green = profit, red = loss, intensity proportional to magnitude. Here shown for a long straddle: the gamma convexity profile (profits for large spot moves) and vega profile (profits for vol spikes) are both visible.

"""
Scenario heatmap — full-reprice P&L grid over (spot shock, vol shock).
Annotated colour-coded grid: green profit, red loss; intensity ∝ magnitude.
"""
import numpy as np
import plotly.graph_objects as go
from scipy.stats import norm

def bs_price(S, K, r, q, sigma, T, call=True):
    sqrtT = np.sqrt(T)
    d1 = (np.log(S/K) + (r-q+0.5*sigma**2)*T)/(sigma*sqrtT)
    d2 = d1 - sigma*sqrtT
    eqT, erT = np.exp(-q*T), np.exp(-r*T)
    return (S*eqT*norm.cdf(d1) - K*erT*norm.cdf(d2) if call
            else K*erT*norm.cdf(-d2) - S*eqT*norm.cdf(-d1))

# Portfolio: long straddle (long call + long put, same strike)
S0, K, r, q, sigma0, T, qty, mult = 100.0, 100.0, 0.05, 0.02, 0.20, 0.5, 10, 100

spot_shocks = np.array([-0.10,-0.075,-0.05,-0.025, 0, 0.025, 0.05, 0.075, 0.10])
vol_shocks  = np.array([-0.10,-0.075,-0.05,-0.025, 0, 0.025, 0.05, 0.075, 0.10])

base = (bs_price(S0, K, r, q, sigma0, T, call=True) +
        bs_price(S0, K, r, q, sigma0, T, call=False)) * qty * mult

pnl = np.zeros((len(spot_shocks), len(vol_shocks)))
for i, ss in enumerate(spot_shocks):
    for j, vs in enumerate(vol_shocks):
        S_sh  = S0*(1+ss)
        sig_sh = max(sigma0+vs, 0.0001)
        shocked = (bs_price(S_sh, K, r, q, sig_sh, T, call=True) +
                   bs_price(S_sh, K, r, q, sig_sh, T, call=False)) * qty * mult
        pnl[i, j] = shocked - base

# Annotation: formatted P&L strings
text = [[f"${pnl[i,j]:+,.0f}" for j in range(len(vol_shocks))]
        for i in range(len(spot_shocks))]

fig = go.Figure(go.Heatmap(
    z=pnl,
    x=[f"{v*100:+.1f}pp" for v in vol_shocks],
    y=[f"{s*100:+.1f}%" for s in spot_shocks],
    colorscale="RdYlGn",
    zmid=0,
    text=text, texttemplate="%{text}",
    colorbar_title="P&L ($)",
))
fig.update_layout(
    title="Scenario P&L Heatmap — Long Straddle (rows: ΔS, cols: Δσ)",
    xaxis_title="Vol Shock (pp)",
    yaxis_title="Spot Shock (%)",
    height=500,
)
fig.show()

5. Dollar Greeks Ladder — Portfolio Risk by Position

Grouped bar chart showing DollarDelta, DollarGamma, and DollarVega contribution per position. Long and short positions net against each other. The net portfolio exposure determines the hedge trade required to flatten each risk factor.

"""
Dollar Greeks ladder — portfolio risk aggregation by instrument.
Sorted bar chart showing DollarDelta, DollarGamma, DollarVega contribution per position.
Desk convention:
  DollarDelta = Δ × S × qty × mult          ($ per $1 spot move)
  DollarGamma = Γ × (0.01S)² × qty × mult  ($ per 1% spot move²)
  DollarVega  = Vega_unit × 0.01 × qty × mult ($ per vol-point)
"""
import numpy as np
import plotly.graph_objects as go
from scipy.stats import norm

def bs_full(S, K, r, q, sigma, T, call):
    sqrtT = np.sqrt(T)
    d1 = (np.log(S/K)+(r-q+0.5*sigma**2)*T)/(sigma*sqrtT)
    d2 = d1-sigma*sqrtT
    phi, eqT = norm.pdf(d1), np.exp(-q*T)
    delta = eqT*norm.cdf(d1) if call else -eqT*norm.cdf(-d1)
    gamma = eqT*phi/(S*sigma*sqrtT)
    vega  = S*eqT*phi*sqrtT
    return delta, gamma, vega

# Example portfolio: mixed calls and puts on same underlying
positions = [
    dict(label="Long ATM Call",  S=100,K=100,r=0.05,q=0.02,sig=0.20,T=1.0,call=True, qty=+50,mult=100),
    dict(label="Short OTM Call", S=100,K=110,r=0.05,q=0.02,sig=0.22,T=1.0,call=True, qty=-30,mult=100),
    dict(label="Long OTM Put",   S=100,K= 90,r=0.05,q=0.02,sig=0.25,T=0.5,call=False,qty=+20,mult=100),
    dict(label="Short ATM Put",  S=100,K=100,r=0.05,q=0.02,sig=0.20,T=0.5,call=False,qty=-10,mult=100),
]

dollar_deltas, dollar_gammas, dollar_vegas, labels = [], [], [], []
for pos in positions:
    d, g, v = bs_full(pos["S"],pos["K"],pos["r"],pos["q"],pos["sig"],pos["T"],pos["call"])
    s = pos["qty"]*pos["mult"]
    dollar_deltas.append(d * pos["S"] * s)
    dollar_gammas.append(g * (0.01*pos["S"])**2 * s)
    dollar_vegas .append(v * 0.01 * s)
    labels.append(pos["label"])

fig = go.Figure()
for name, vals, color in [
    ("DollarDelta ($/$1 spot)", dollar_deltas, "#3b82f6"),
    ("DollarGamma ($/1% spot²)", dollar_gammas, "#10b981"),
    ("DollarVega ($/vol-pt)",  dollar_vegas,  "#f59e0b"),
]:
    fig.add_trace(go.Bar(name=name, x=labels, y=vals,
                          marker_color=[color if v>=0 else "#ef4444" for v in vals],
                          text=[f"${v:+,.0f}" for v in vals], textposition="outside"))

fig.update_layout(
    barmode="group",
    title="Portfolio Dollar Greeks Ladder",
    yaxis_title="Dollar Risk ($)",
    xaxis_title="Position",
    height=420,
    legend=dict(orientation="h", yanchor="bottom", y=1.02),
)
fig.show()

← All Projects LSMC Pricer →SVI Calibrator →