Vanilla Options Pricing — Black-Scholes, Binomial Tree, and Greeks:

A Python module that prices European and American vanilla options under two models built from scratch — closed-form Black-Scholes (1973) and the Cox-Ross-Rubinstein (1979) binomial tree — and computes the full set of sensitivities (delta, gamma, vega, theta, and rho) analytically and by finite difference. Pricing and Greek surfaces are generated over user-defined moneyness and time-to-expiry grids. A scenario engine quantifies P&L attribution for common structures (call spreads, straddles, strangles) under parallel implied-vol shocks, spot moves, and theta decay. Output is exposed through a typer CLI and a Streamlit dashboard with real-time sliders for spot, strike, maturity, rate, and volatility. Market price histories are sourced from yfinance; implied-vol proxies are inferred from historical return distributions and stored in DuckDB. Built on Python 3.11+ using pandas, numpy, scipy, plotly, streamlit, and duckdb; packaged with hatchling and tested with pytest against known Black-Scholes reference values.

%==========%

I. Streamlit Dashboard (app.py):

The dashboard exposes the full pricing, Greeks, surface, and scenario engine via a browser UI. The sidebar controls the underlying ticker, spot, strike, maturity, rate, volatility, and option type. All four tabs read from a single shared @st.cache_data call that re-runs only when the sidebar inputs change. Tab 3 offers a toggle between a 3-D surface and a 2-D heatmap for each Greek. Tab 4 lets the user select a pre-built structure or define custom legs, then apply a spot move and a vol shock simultaneously.

# app.py (structure)
st.set_page_config(page_title="Options Pricing", layout="wide")

with st.sidebar:
    ticker  = st.text_input("Underlying (yfinance)", "SPY")
    S       = st.slider("Spot (S)",              50.0,  200.0, 100.0, step=0.5)
    K       = st.slider("Strike (K)",            50.0,  200.0, 100.0, step=0.5)
    T       = st.slider("Time to Expiry (years)", 0.05,   3.0,   1.0, step=0.05)
    r       = st.slider("Risk-Free Rate (%)",     0.0,   10.0,   5.0, step=0.25) / 100.0
    sigma   = st.slider("Implied Vol (%)",        1.0,  100.0,  20.0, step=0.5)  / 100.0
    otype   = st.radio("Option Type", ["call", "put"], horizontal=True)
    n_steps = st.select_slider("Binomial Steps (N)", [50, 100, 200, 500], value=200)

tab1, tab2, tab3, tab4 = st.tabs(
    ["💲 Pricing", "🔢 Greeks", "📈 Surfaces", "⚡ Scenarios"]
)
# tab1: BS and CRR prices side-by-side; implied vol inversion from a market price input;
#       put-call parity check card showing |C - P - (S - K·e^{-rT})|
# tab2: Greek table comparing analytical BS vs finite-diff CRR; bar chart for each Greek
# tab3: surface selector (price / delta / gamma / vega / theta), toggle 3-D surface ↔ 2-D heatmap
# tab4: structure builder (call_spread / put_spread / straddle / strangle / custom legs),
#       dS + dσ + dt sliders, full P&L decomposition table + bar chart of attribution
  

# 1. Navigate to the project root
cd assets/projects/options_pricing

# 2. Create and activate a virtual environment (Python 3.11+ required)
python -m venv .venv
# Windows:
.venv\Scripts\Activate.ps1
# macOS / Linux:
source .venv/bin/activate

# 3. Install the package and all dependencies
pip install -e ".[dev]"

# 4. (Optional) copy the env template and add an Alpha Vantage API key
cp .env.example .env

# 5. Pre-fetch price data and build the DuckDB snapshot table
python scripts/download_data.py

# 6. Launch
streamlit run src/options_pricing/app.py
  

The dashboard below runs entirely in the browser via stlite (Streamlit on WebAssembly — no server required). First load downloads Pyodide and may take 20–40 seconds; subsequent loads are cached. The live-data tab is not available in the browser version.

%==========%

II. Project Layout:

options-pricing/
├── pyproject.toml                          # Build config, deps, ruff + pytest settings
├── .env.example                            # API key template (Alpha Vantage)
├── data/                                   # Populated by scripts/download_data.py
├── notebooks/
│   └── analytical_vs_numerical.ipynb       # BS vs CRR convergence, stability, runtimes
├── scripts/
│   └── download_data.py                    # Fetches prices via yfinance; writes vol-proxy table to DuckDB
├── src/options_pricing/
│   ├── data/
│   │   ├── schemas.py                      # Pydantic v2: OptionContract, MarketSnapshot
│   │   ├── fetchers.py                     # yfinance + Alpha Vantage price fetchers
│   │   └── loaders.py                      # DuckDB loader / CSV fallback
│   ├── models/
│   │   ├── black_scholes.py                # Closed-form BS for European options + IV inversion
│   │   └── binomial.py                     # CRR binomial tree: European and American exercise
│   ├── greeks/
│   │   └── analytical.py                   # Analytical BS Greeks + model-agnostic finite-difference Greeks
│   ├── surfaces/
│   │   └── builder.py                      # Pricing and Greek surfaces over (S/K × T) grid
│   ├── scenarios/
│   │   └── analysis.py                     # Scenario P&L and Greek attribution for spreads / straddles
│   ├── report/
│   │   └── plots.py                        # Plotly: 3-D surfaces, Greek heatmaps, payoff diagrams
│   ├── cli.py                              # Typer CLI: price | greeks | surface | scenario
│   └── app.py                              # Streamlit dashboard (4 tabs)
└── tests/
    ├── conftest.py                         # Fixtures: ATM call/put params, known BS reference values
    ├── test_black_scholes.py               # Put-call parity, boundary conditions, IV roundtrip
    ├── test_binomial.py                    # CRR convergence to BS, American exercise premium
    ├── test_greeks.py                      # FD vs analytical Greeks, Greek sign invariants
    └── test_scenarios.py                   # Scenario attribution vs full revaluation
  

%==========%

III. Data Layer — schemas.py & loaders.py:

OptionContract captures the parameters needed to price a single leg; MarketSnapshot stores the spot and historical-vol proxies fetched for a given underlying and date. loaders.py first queries DuckDB (built by download_data.py); it falls back to a CSV if the database is absent. The 21-day and 63-day historical-vol columns serve as proxy inputs to the pricer when no exchange-quoted implied vol is available.

# schemas.py
from __future__ import annotations
from datetime import date
from typing import Literal
from pydantic import BaseModel, field_validator

OptionType    = Literal["call", "put"]
ExerciseStyle = Literal["european", "american"]

class OptionContract(BaseModel):
    underlying:     str
    option_type:    OptionType
    exercise_style: ExerciseStyle = "european"
    K:     float          # strike price
    T:     float          # time to expiry (years)
    r:     float = 0.05   # risk-free rate, continuously compounded
    sigma: float          # implied or proxy vol (annualised, decimal)

    @field_validator("K", "T", "sigma")
    @classmethod
    def must_be_positive(cls, v: float) -> float:
        if v <= 0:
            raise ValueError(f"must be positive, got {v}")
        return v

class MarketSnapshot(BaseModel):
    underlying:    str
    spot:          float
    date:          date
    hist_vol_21d:  float   # 21-day rolling std proxy (annualised)
    hist_vol_63d:  float   # 63-day proxy
  

# loaders.py
from __future__ import annotations
from pathlib import Path
import duckdb
import pandas as pd

_DB = Path("data/options_pricing.duckdb")

def load_snapshot(underlying: str, as_of: str | None = None) -> pd.DataFrame:
    """Return most recent MarketSnapshot row for a given underlying from DuckDB."""
    con = duckdb.connect(str(_DB), read_only=True)
    if as_of:
        q = ("SELECT * FROM market_snapshots WHERE underlying = ? AND date <= ?"
             " ORDER BY date DESC LIMIT 1")
        return con.execute(q, [underlying, as_of]).df()
    q = "SELECT * FROM market_snapshots WHERE underlying = ? ORDER BY date DESC LIMIT 1"
    return con.execute(q, [underlying]).df()

def load_prices(path: str | Path) -> pd.DataFrame:
    path = Path(path)
    if not path.exists():
        raise FileNotFoundError(f"Price file not found: {path}")
    df = pd.read_csv(path, index_col=0, parse_dates=True)
    return df.sort_index().astype(float)
  

%==========%

IV. Model I — Black-Scholes (black_scholes.py):

The Black-Scholes (1973) formula gives a closed-form price for a European option under geometric Brownian motion with constant drift \(\mu\) and volatility \(\sigma\). The risk-neutral measure replaces \(\mu\) with the risk-free rate \(r\), yielding the standard \(d_1\) and \(d_2\) statistics: $$d_1 = \frac{\ln(S/K) + (r + \tfrac{1}{2}\sigma^2)\,T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T}$$ from which the call and put prices follow immediately: $$C = S\,\Phi(d_1) - K e^{-rT}\,\Phi(d_2), \qquad P = K e^{-rT}\,\Phi(-d_2) - S\,\Phi(-d_1)$$ where \(\Phi\) is the standard normal CDF. Implied volatility inversion uses Brent’s method, which is globally convergent on a bounded bracket and avoids the numerical instability of Newton–Raphson near the strike wings where vega approaches zero.

# models/black_scholes.py
from __future__ import annotations
import numpy as np
from scipy.stats import norm
from scipy.optimize import brentq
from dataclasses import dataclass
from typing import Literal

OptionType = Literal["call", "put"]

@dataclass
class BSResult:
    option_type: OptionType
    price: float
    d1: float; d2: float
    S: float; K: float; T: float; r: float; sigma: float

def _d1d2(S: float, K: float, T: float, r: float, sigma: float) -> tuple[float, float]:
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    return d1, d1 - sigma * np.sqrt(T)

def bs_price(S: float, K: float, T: float, r: float, sigma: float,
             option_type: OptionType = "call") -> float:
    """Black-Scholes (1973) closed-form price for a European option."""
    if T <= 0.0:
        return max(S - K, 0.0) if option_type == "call" else max(K - S, 0.0)
    d1, d2 = _d1d2(S, K, T, r, sigma)
    if option_type == "call":
        return float(S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2))
    return float(K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1))

def bs_result(S: float, K: float, T: float, r: float, sigma: float,
              option_type: OptionType = "call") -> BSResult:
    d1, d2 = _d1d2(S, K, T, r, sigma)
    return BSResult(option_type=option_type, price=bs_price(S, K, T, r, sigma, option_type),
                    d1=d1, d2=d2, S=S, K=K, T=T, r=r, sigma=sigma)

def implied_vol(market_price: float, S: float, K: float, T: float, r: float,
                option_type: OptionType = "call") -> float:
    """Brent's-method implied-vol inversion. Returns NaN if the price violates no-arbitrage bounds."""
    disc      = np.exp(-r * T)
    lb        = max(S - K * disc, 0.0) if option_type == "call" else max(K * disc - S, 0.0)
    if market_price <= lb or market_price >= S:
        return np.nan
    try:
        return brentq(
            lambda s: bs_price(S, K, T, r, s, option_type) - market_price,
            1e-6, 10.0, xtol=1e-7,
        )
    except ValueError:
        return np.nan
  

%==========%

V. Model II — CRR Binomial Tree (binomial.py):

The Cox-Ross-Rubinstein (1979) binomial model discretises the asset path into \(N\) time steps of length \(\Delta t = T/N\). Each step the stock moves up by factor \(u\) or down by factor \(d = 1/u\), chosen to match the GBM volatility: $$u = e^{\sigma\sqrt{\Delta t}}, \qquad d = \frac{1}{u}, \qquad p = \frac{e^{r\Delta t} - d}{u - d}$$ where \(p\) is the risk-neutral up-probability. The option value is computed by backward induction from the \(N+1\) terminal payoffs: $$V_i^{(n)} = e^{-r\Delta t}\!\left[p\,V_{i}^{(n+1)} + (1-p)\,V_{i+1}^{(n+1)}\right]$$ For American options the early-exercise constraint is imposed at each node: $$V_i^{(n)} \leftarrow \max\!\left(V_i^{(n)},\;\text{intrinsic}_i^{(n)}\right)$$ As \(N \to \infty\) the CRR price converges to the Black-Scholes price for European options, providing a direct numerical check on the analytical formula. The implementation vectorises the terminal payoff array and loops backward through the tree without allocating a full \((N+1)\times(N+1)\) matrix.

# models/binomial.py
from __future__ import annotations
import numpy as np
from dataclasses import dataclass
from typing import Literal

OptionType    = Literal["call", "put"]
ExerciseStyle = Literal["european", "american"]

@dataclass
class BinomialResult:
    option_type:    OptionType
    exercise_style: ExerciseStyle
    price: float
    S: float; K: float; T: float; r: float; sigma: float; N: int

def crr_price(S: float, K: float, T: float, r: float, sigma: float,
              N: int = 200, option_type: OptionType = "call",
              exercise_style: ExerciseStyle = "european") -> float:
    """Cox-Ross-Rubinstein (1979) binomial tree. O(N) space via single-vector backward induction."""
    dt   = T / N
    u    = np.exp(sigma * np.sqrt(dt))
    d    = 1.0 / u
    p    = (np.exp(r * dt) - d) / (u - d)
    disc = np.exp(-r * dt)

    # Terminal stock prices  (j = 0 … N: j up-moves, N-j down-moves)
    j  = np.arange(N + 1, dtype=float)
    ST = S * u ** (N - j) * d ** j

    # Terminal payoffs
    V = np.maximum(ST - K, 0.0) if option_type == "call" else np.maximum(K - ST, 0.0)

    # Backward induction
    for step in range(N - 1, -1, -1):
        V = disc * (p * V[:-1] + (1.0 - p) * V[1:])
        if exercise_style == "american":
            j_s = np.arange(step + 1, dtype=float)
            St  = S * u ** (step - j_s) * d ** j_s
            ev  = np.maximum(St - K, 0.0) if option_type == "call" else np.maximum(K - St, 0.0)
            V   = np.maximum(V, ev)

    return float(V[0])

def crr_result(S: float, K: float, T: float, r: float, sigma: float,
               N: int = 200, option_type: OptionType = "call",
               exercise_style: ExerciseStyle = "european") -> BinomialResult:
    return BinomialResult(
        option_type=option_type, exercise_style=exercise_style,
        price=crr_price(S, K, T, r, sigma, N, option_type, exercise_style),
        S=S, K=K, T=T, r=r, sigma=sigma, N=N,
    )
  

%==========%

VI. Greeks — Analytical Sensitivities (greeks/analytical.py):

All five first-order Greeks are derived analytically from the Black-Scholes formula. Gamma is identical for calls and puts by put-call parity. Theta and rho have opposite signs across the two option types. Vega is scaled to per 1 vol-point (\(\partial C / \partial\sigma / 100\)) and theta to per calendar day (\(\partial C / \partial t / 365\)) so that the numbers are directly interpretable in a dashboard context.

# greeks/analytical.py
from __future__ import annotations
import numpy as np
from scipy.stats import norm
from dataclasses import dataclass
from typing import Callable, Literal

OptionType = Literal["call", "put"]

@dataclass
class Greeks:
    delta: float
    gamma: float
    vega:  float   # per 1 vol-point  (∂price/∂σ / 100)
    theta: float   # per calendar day  (∂price/∂t / 365)
    rho:   float   # per 1 bp          (∂price/∂r / 100)

def bs_greeks(S: float, K: float, T: float, r: float, sigma: float,
              option_type: OptionType = "call") -> Greeks:
    """Analytical Black-Scholes Greeks for European vanilla options."""
    sqrt_T = np.sqrt(T)
    d1     = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * sqrt_T)
    d2     = d1 - sigma * sqrt_T
    phi_d1 = norm.pdf(d1)
    disc   = np.exp(-r * T)

    gamma = phi_d1 / (S * sigma * sqrt_T)
    vega  = S * phi_d1 * sqrt_T / 100.0

    if option_type == "call":
        delta = norm.cdf(d1)
        theta = (-(S * phi_d1 * sigma) / (2.0 * sqrt_T) - r * K * disc * norm.cdf(d2)) / 365.0
        rho   =   K * T * disc * norm.cdf(d2)  / 100.0
    else:
        delta = norm.cdf(d1) - 1.0
        theta = (-(S * phi_d1 * sigma) / (2.0 * sqrt_T) + r * K * disc * norm.cdf(-d2)) / 365.0
        rho   = -K * T * disc * norm.cdf(-d2) / 100.0

    return Greeks(delta=delta, gamma=gamma, vega=vega, theta=theta, rho=rho)


def fd_greeks(price_fn: Callable, S: float, K: float, T: float, r: float,
              sigma: float, option_type: OptionType = "call",
              dS: float = 0.01, dsigma: float = 0.001,
              dr: float = 0.0001, dt: float = 1.0 / 365) -> Greeks:
    """Central finite-difference Greeks — model-agnostic; works with CRR or any pricer."""
    p0    = price_fn(S,      K, T,              r,      sigma,          option_type=option_type)
    p_up  = price_fn(S + dS, K, T,              r,      sigma,          option_type=option_type)
    p_dn  = price_fn(S - dS, K, T,              r,      sigma,          option_type=option_type)
    delta = (p_up - p_dn) / (2.0 * dS)
    gamma = (p_up - 2.0 * p0 + p_dn) / dS**2
    vega  = (price_fn(S, K, T, r, sigma + dsigma, option_type=option_type)
           - price_fn(S, K, T, r, sigma - dsigma, option_type=option_type)) / (2.0 * dsigma * 100.0)
    theta = -(price_fn(S, K, max(T - dt, 1e-5), r, sigma, option_type=option_type) - p0)
    rho   = (price_fn(S, K, T, r + dr, sigma, option_type=option_type)
           - price_fn(S, K, T, r - dr, sigma, option_type=option_type)) / (2.0 * dr * 100.0)
    return Greeks(delta=delta, gamma=gamma, vega=vega, theta=theta, rho=rho)
  

%==========%

VII. Pricing & Greek Surfaces (surfaces/builder.py):

Surfaces are evaluated over a two-dimensional grid: moneyness \(S/K\) on the vertical axis and time to expiry \(T\) on the horizontal axis. Moneyness is used rather than raw spot so that the surface is strike-independent and directly comparable across underlyings. build_price_surface and build_greek_surface both return a typed Surface dataclass consumed by plots.py, keeping the computation layer decoupled from the visualisation layer. The 3-D surface and 2-D heatmap views are generated from the same underlying Surface.values array.

# surfaces/builder.py
from __future__ import annotations
import numpy as np
from dataclasses import dataclass
from typing import Callable

@dataclass
class Surface:
    values:    np.ndarray   # shape (n_S, n_T)
    moneyness: np.ndarray   # S/K grid
    expiries:  np.ndarray   # T grid (years)
    label:     str

def build_price_surface(
    price_fn: Callable,
    K: float = 100.0, r: float = 0.05, sigma: float = 0.20,
    S_range: tuple[float, float] = (0.70, 1.30),
    T_range: tuple[float, float] = (0.05, 2.0),
    n_S: int = 50, n_T: int = 40,
    option_type: str = "call",
) -> Surface:
    """Price surface over a (moneyness × time-to-expiry) grid."""
    moneyness = np.linspace(S_range[0], S_range[1], n_S)
    expiries  = np.linspace(T_range[0], T_range[1], n_T)
    Z = np.array([
        [price_fn(m * K, K, t, r, sigma, option_type=option_type) for t in expiries]
        for m in moneyness
    ])
    return Surface(values=Z, moneyness=moneyness, expiries=expiries,
                   label=f"BS {option_type.capitalize()} Price (σ={sigma:.0%})")

def build_greek_surface(
    greek_fn: Callable,
    greek_name: str,
    K: float = 100.0, r: float = 0.05, sigma: float = 0.20,
    S_range: tuple[float, float] = (0.70, 1.30),
    T_range: tuple[float, float] = (0.05, 2.0),
    n_S: int = 50, n_T: int = 40,
    option_type: str = "call",
) -> Surface:
    """Greek surface over the same (moneyness × time-to-expiry) grid."""
    moneyness = np.linspace(S_range[0], S_range[1], n_S)
    expiries  = np.linspace(T_range[0], T_range[1], n_T)
    attr = greek_name.lower()
    Z = np.array([
        [getattr(greek_fn(m * K, K, t, r, sigma, option_type=option_type), attr)
         for t in expiries]
        for m in moneyness
    ])
    return Surface(values=Z, moneyness=moneyness, expiries=expiries,
                   label=f"{greek_name.capitalize()} surface ({option_type})")
  

# report/plots.py (surface and heatmap sections)
import numpy as np
import plotly.graph_objects as go
from .surfaces.builder import Surface

def plot_surface(surface: Surface, colorscale: str = "RdBu_r") -> go.Figure:
    fig = go.Figure(data=go.Surface(
        x=surface.expiries, y=surface.moneyness, z=surface.values,
        colorscale=colorscale, showscale=True,
    ))
    fig.update_layout(
        title=surface.label,
        scene=dict(
            xaxis_title="Time to Expiry (years)",
            yaxis_title="Moneyness (S/K)",
            zaxis_title=surface.label.split(" surface")[0],
        ),
        template="plotly_white", margin=dict(l=0, r=0, b=0, t=40),
    )
    return fig

def plot_greek_heatmap(surface: Surface) -> go.Figure:
    fig = go.Figure(data=go.Heatmap(
        x=surface.expiries, y=surface.moneyness, z=surface.values,
        colorscale="RdBu", zmid=0.0,
    ))
    fig.update_layout(
        title=surface.label,
        xaxis_title="Time to Expiry (years)",
        yaxis_title="Moneyness (S/K)",
        template="plotly_white",
    )
    return fig

def plot_payoff(legs: list, S_range: tuple[float, float] = (60.0, 140.0),
                n: int = 300) -> go.Figure:
    from .models.black_scholes import bs_price
    spots  = np.linspace(S_range[0], S_range[1], n)
    payoff = np.array([
        sum(leg.quantity * bs_price(s, leg.K, 0.0, 0.0, 1e-5, leg.option_type)
            for leg in legs)
        for s in spots
    ])
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=spots, y=payoff, name="Payoff at expiry",
                             line=dict(color="#636EFA", width=2)))
    fig.add_hline(y=0, line_dash="dash", line_color="grey", line_width=1)
    fig.update_layout(xaxis_title="Spot (S)", yaxis_title="P&L",
                      title="Payoff Diagram at Expiry", template="plotly_white")
    return fig
  

The interactive 3-D surfaces below are generated from the same build_price_surface and build_greek_surface functions used by the Streamlit dashboard. Each surface spans a 60×50 moneyness–expiry grid (moneyness 0.70–1.30, T 0.05–2.0 yr) at a fixed 20% implied vol. Rotate or zoom with the mouse; hover for exact values.

%==========%

VIII. Scenario Analysis (scenarios/analysis.py):

The scenario engine aggregates Greeks across a multi-leg structure and decomposes the P&L of any joint shock into first- and second-order Taylor contributions, then compares the sum against a full revaluation. For a call spread facing a +20% implied-vol shock with spot up $5, the delta, gamma, vega, and theta contributions can be reported separately, making it clear which Greeks are driving P&L. Four pre-built structures are available; arbitrary legs can be passed at the call site.

# scenarios/analysis.py
from __future__ import annotations
import numpy as np
from dataclasses import dataclass
from typing import Callable

@dataclass
class Leg:
    option_type: str   # "call" or "put"
    K:           float
    quantity:    float  # +1 long, −1 short; fractional for notional scaling

@dataclass
class ScenarioResult:
    structure:   str
    full_pnl:    float   # full revaluation P&L
    delta_pnl:   float   # Δ · dS
    gamma_pnl:   float   # ½ · Γ · dS²
    vega_pnl:    float   # ν · d(σ × 100)
    theta_pnl:   float   # Θ · dt_days
    greek_total: float   # sum of Greek contributions (Taylor approximation)

STRUCTURES: dict[str, list[Leg]] = {
    "call_spread": [Leg("call", 100.0, +1), Leg("call", 110.0, -1)],
    "put_spread":  [Leg("put",   90.0, -1), Leg("put",  100.0, +1)],
    "straddle":    [Leg("call", 100.0, +1), Leg("put",  100.0, +1)],
    "strangle":    [Leg("call", 110.0, +1), Leg("put",   90.0, +1)],
}

def _portfolio_value(legs: list[Leg], price_fn: Callable,
                     S: float, T: float, r: float, sigma: float) -> float:
    return sum(
        leg.quantity * price_fn(S, leg.K, T, r, sigma, option_type=leg.option_type)
        for leg in legs
    )

def run_scenario(
    legs: list[Leg],
    price_fn: Callable,
    greek_fn: Callable,
    S: float, T: float, r: float, sigma: float,
    dS: float = 0.0, dsigma: float = 0.0, dt_days: float = 0.0,
    structure_name: str = "custom",
) -> ScenarioResult:
    """Compute scenario P&L with Taylor attribution (Δ, Γ, ν, Θ) vs full revaluation."""
    V0 = _portfolio_value(legs, price_fn, S, T, r, sigma)

    # Aggregate Greeks across legs
    delta = vega = theta = gamma = 0.0
    for leg in legs:
        g      = greek_fn(S, leg.K, T, r, sigma, option_type=leg.option_type)
        delta += leg.quantity * g.delta
        vega  += leg.quantity * g.vega
        theta += leg.quantity * g.theta
        gamma += leg.quantity * g.gamma

    # Taylor contributions
    delta_pnl = delta * dS
    gamma_pnl = 0.5 * gamma * dS**2
    vega_pnl  = vega  * (dsigma * 100.0)   # ν is per vol-point; dsigma is decimal
    theta_pnl = theta * dt_days

    # Full revaluation
    S1     = S + dS
    sigma1 = sigma + dsigma
    T1     = max(T - dt_days / 365.0, 1e-5)
    V1     = _portfolio_value(legs, price_fn, S1, T1, r, sigma1)

    return ScenarioResult(
        structure=structure_name, full_pnl=V1 - V0,
        delta_pnl=delta_pnl, gamma_pnl=gamma_pnl,
        vega_pnl=vega_pnl,   theta_pnl=theta_pnl,
        greek_total=delta_pnl + gamma_pnl + vega_pnl + theta_pnl,
    )
  

Example: call spread (\(K_1=100, K_2=110\)), \(S=100\), \(T=0.5\text{y}\), \(\sigma=20\%\), implied vol shock \(+20\%\) (\(\Delta\sigma=+0.20\)), spot up $5:

from options_pricing.scenarios.analysis import run_scenario, STRUCTURES
from options_pricing.models.black_scholes import bs_price
from options_pricing.greeks.analytical import bs_greeks

result = run_scenario(
    legs=STRUCTURES["call_spread"],
    price_fn=bs_price, greek_fn=bs_greeks,
    S=100.0, T=0.5, r=0.05, sigma=0.20,
    dS=5.0, dsigma=0.20, dt_days=0.0,
    structure_name="call_spread",
)
# result.full_pnl    →  full revaluation
# result.delta_pnl   →  Δ · 5
# result.vega_pnl    →  ν · 20   (20 vol-points)
# result.gamma_pnl   →  ½ · Γ · 25
  

%==========%

IX. CLI — cli.py:

Four subcommands share common parameter flags. When --model crr is selected, finite-difference Greeks are computed against the binomial pricer rather than the analytical BS formulas.

# cli.py (key subcommands)
@app.command("price")
def cmd_price(
    S: float = 100.0, K: float = 100.0, T: float = 1.0,
    r: float = 0.05,  sigma: float = 0.20,
    model: str = "bs", style: str = "european", option_type: str = "call",
) -> None:
    """Price a European or American option under BS or CRR binomial tree."""
    if model == "bs":
        res = bs_result(S, K, T, r, sigma, option_type)
        console.print(f"[green]BS {option_type.upper()} price:[/] {res.price:.6f}")
        console.print(f"  d1={res.d1:.4f}  d2={res.d2:.4f}")
    else:
        res = crr_result(S, K, T, r, sigma, N=500, option_type=option_type, exercise_style=style)
        console.print(f"[green]CRR {option_type.upper()} ({style}) price:[/] {res.price:.6f}")

@app.command("greeks")
def cmd_greeks(
    S: float = 100.0, K: float = 100.0, T: float = 1.0,
    r: float = 0.05, sigma: float = 0.20, option_type: str = "call",
    mode: str = "analytical",
) -> None:
    """Print analytical or finite-difference Black-Scholes Greeks."""
    g = bs_greeks(S, K, T, r, sigma, option_type) if mode == "analytical" \
        else fd_greeks(bs_price, S, K, T, r, sigma, option_type)
    table = rich.table.Table(title=f"Greeks ({mode})")
    for field, val in [("Delta", g.delta), ("Gamma", g.gamma),
                        ("Vega",  g.vega),  ("Theta", g.theta), ("Rho", g.rho)]:
        table.add_row(field, f"{val:.6f}")
    console.print(table)

@app.command("surface")
def cmd_surface(
    K: float = 100.0, r: float = 0.05, sigma: float = 0.20,
    greek: str = "price", option_type: str = "call", output: str = "surface.html",
) -> None:
    """Export a Plotly pricing or Greek surface to an HTML file."""
    if greek == "price":
        surf = build_price_surface(bs_price, K=K, r=r, sigma=sigma, option_type=option_type)
        fig  = plot_surface(surf)
    else:
        surf = build_greek_surface(bs_greeks, greek, K=K, r=r, sigma=sigma, option_type=option_type)
        fig  = plot_greek_heatmap(surf)
    fig.write_html(output)
    console.print(f"[green]Surface written to {output}[/]")

@app.command("scenario")
def cmd_scenario(
    structure: str = "call_spread",
    S: float = 100.0, T: float = 1.0, r: float = 0.05, sigma: float = 0.20,
    dS: float = 0.0, dsigma: float = 0.0, dt_days: float = 0.0,
) -> None:
    """Run P&L scenario on a pre-built structure and print Greek attribution."""
    legs   = STRUCTURES[structure]
    result = run_scenario(legs, bs_price, bs_greeks, S, T, r, sigma,
                          dS=dS, dsigma=dsigma, dt_days=dt_days,
                          structure_name=structure)
    # prints: full_pnl | delta_pnl | gamma_pnl | vega_pnl | theta_pnl | greek_total
  

# Install
pip install -e ".[dev]"

# Price an ATM call under Black-Scholes
options-price price --S 100 --K 100 --T 1 --r 0.05 --sigma 0.20 --model bs

# Price an American put under CRR (500 steps)
options-price price --S 100 --K 105 --T 0.5 --sigma 0.25 --model crr --style american --option-type put

# Print analytical Greeks for an OTM call
options-price greeks --S 100 --K 110 --T 0.5 --sigma 0.20

# Export a gamma surface to HTML
options-price surface --K 100 --sigma 0.20 --greek gamma --output gamma_surface.html

# Scenario: call spread, spot +$5, implied vol +20%
options-price scenario --structure call_spread --S 100 --T 0.5 --sigma 0.20 --dS 5 --dsigma 0.20
  

%==========%

X. Notebook: Analytical vs Numerical (analytical_vs_numerical.ipynb):

The notebook benchmarks the Black-Scholes closed form against the CRR binomial tree across three dimensions: convergence (error as a function of \(N\)), numerical stability (behaviour near zero time-to-expiry and deep in/out of the money), and runtime (wall-clock seconds for \(N \in \{50, 100, 200, 500, 1000, 2000\}\)). Key findings reproduced in the notebook:

The convergence plot shows the familiar odd-even oscillation of CRR error as \(N\) increases, which is mitigated by using the average of the \(N\)-step and \((N+1)\)-step prices (Richardson extrapolation). The runtime plot demonstrates that BS is three to four orders of magnitude faster, justifying the analytical formula for interactive dashboards and Monte Carlo Greeks, while the binomial tree remains the reference implementation for American exercise and barrier-like path constraints.

%==========%

XI. Test Suite:

All tests are fully offline. The atm fixture in conftest.py provides a canonical ATM European call parameter set \((S=K=100,\;T=1\text{y},\;r=5\%,\;\sigma=20\%)\) against which known BS reference values can be checked. Tests verify mathematical invariants — put-call parity to machine precision, delta in \((0,1)\), gamma and vega positivity, theta negativity for long options, American put early-exercise premium ≥ 0, CRR convergence to BS, and implied-vol inversion roundtrip accuracy.

# conftest.py
@pytest.fixture
def atm() -> dict:
    """ATM European params: S=K=100, T=1y, r=5%, σ=20%."""
    return dict(S=100.0, K=100.0, T=1.0, r=0.05, sigma=0.20)
  

# test_black_scholes.py
class TestBlackScholes:
    def test_put_call_parity(self, atm):
        """C − P = S − K·e^{−rT} holds to machine precision."""
        C   = bs_price(**atm, option_type="call")
        P   = bs_price(**atm, option_type="put")
        fwd = atm["S"] - atm["K"] * np.exp(-atm["r"] * atm["T"])
        assert abs((C - P) - fwd) < 1e-10

    def test_at_expiry_call_intrinsic(self):
        assert abs(bs_price(110.0, 100.0, 0.0, 0.05, 0.20, "call") - 10.0) < 1e-10

    def test_deep_otm_call_near_zero(self):
        assert bs_price(50.0, 100.0, 0.1, 0.05, 0.20, "call") < 1e-5

    def test_implied_vol_roundtrip(self, atm):
        price = bs_price(**atm, option_type="call")
        iv    = implied_vol(price, **atm, option_type="call")
        assert abs(iv - atm["sigma"]) < 1e-6

# test_binomial.py
class TestBinomial:
    def test_convergence_to_bs(self, atm):
        """CRR with N=1000 agrees with BS to within 0.01."""
        bs_val  = bs_price(**atm, option_type="call")
        crr_val = crr_price(**atm, N=1000, option_type="call", exercise_style="european")
        assert abs(crr_val - bs_val) < 0.01

    def test_american_put_ge_european_put(self, atm):
        eur = crr_price(**atm, N=200, option_type="put", exercise_style="european")
        ame = crr_price(**atm, N=200, option_type="put", exercise_style="american")
        assert ame >= eur - 1e-10

    def test_american_call_no_dividend_equals_european(self, atm):
        eur = crr_price(**atm, N=200, option_type="call", exercise_style="european")
        ame = crr_price(**atm, N=200, option_type="call", exercise_style="american")
        assert abs(ame - eur) < 1e-8

# test_greeks.py
class TestGreeks:
    def test_delta_call_in_range(self, atm):
        assert 0.0 < bs_greeks(**atm, option_type="call").delta < 1.0

    def test_delta_put_in_range(self, atm):
        assert -1.0 < bs_greeks(**atm, option_type="put").delta < 0.0

    def test_gamma_positive(self, atm):
        assert bs_greeks(**atm, option_type="call").gamma > 0.0

    def test_vega_positive(self, atm):
        assert bs_greeks(**atm, option_type="call").vega > 0.0

    def test_theta_negative_long_call(self, atm):
        assert bs_greeks(**atm, option_type="call").theta < 0.0

    def test_fd_delta_matches_analytical(self, atm):
        analytical = bs_greeks(**atm, option_type="call").delta
        fd         = fd_greeks(bs_price, **atm, option_type="call").delta
        assert abs(analytical - fd) < 1e-4
  

%==========%

XII. Notebook — Analytical vs Numerical:

The notebook below benchmarks the Black-Scholes closed form against the CRR binomial tree across four experiments: (1) convergence of CRR price to the BS exact value as N grows; (2) wall-clock runtime comparison; (3) numerical stability at extreme parameters (short-dated ATM, deep OTM, deep ITM); and (4) Richardson extrapolation to suppress the CRR odd-even oscillation. Section 5 quantifies the American early-exercise premium for puts over a range of spot prices.

%==========%

XIII. Configuration & Data Sources:

Team:

Theodosios Dimitrasopoulos, personal project.

Tools & methods:

Python 3.11, pandas, NumPy, SciPy (norm CDF/PDF, Brent’s method), Plotly (interactive 3-D surfaces and heatmaps), Streamlit (4-tab dashboard), yfinance / Alpha Vantage (market data), DuckDB (local snapshot store), Pydantic v2 (schema validation), Typer (CLI), rich (console output), pytest, ruff, hatchling. Pricing methodology: Black-Scholes (1973) closed form for European options; Cox-Ross-Rubinstein (1979) binomial tree for European and American exercise; analytical first- and second-order Greeks (delta, gamma, vega, theta, rho) and central finite-difference Greeks; Brent’s-method implied-vol inversion; pricing and Greek surfaces over a moneyness × expiry grid; scenario P&L attribution decomposed into delta, gamma, vega, and theta contributions vs full revaluation.