vault-dash/app/core/pricing/american_pricing.py

from __future__ import annotations

from dataclasses import dataclass
from datetime import date, timedelta
from typing import Literal

import QuantLib as ql

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

DEFAULT_RISK_FREE_RATE: float = 0.045
DEFAULT_VOLATILITY: float = 0.16
DEFAULT_DIVIDEND_YIELD: float = 0.0
DEFAULT_GLD_PRICE: float = 460.0


@dataclass(frozen=True)
class AmericanOptionInputs:
    """Inputs for American option pricing via a binomial tree.

    This module is intended primarily for GLD protective puts, where early
    exercise can matter in stressed scenarios.

    Example:
        >>> params = AmericanOptionInputs(
        ...     spot=460.0,
        ...     strike=420.0,
        ...     time_to_expiry=0.5,
        ...     option_type="put",
        ... )
        >>> params.steps
        500
    """

    spot: float = DEFAULT_GLD_PRICE
    strike: float = DEFAULT_GLD_PRICE
    time_to_expiry: float = 0.5
    risk_free_rate: float = DEFAULT_RISK_FREE_RATE
    volatility: float = DEFAULT_VOLATILITY
    option_type: OptionType = "put"
    dividend_yield: float = DEFAULT_DIVIDEND_YIELD
    steps: int = 500
    valuation_date: date | None = None
    tree: str = "crr"


@dataclass(frozen=True)
class AmericanPricingResult:
    """American option price and finite-difference Greeks."""

    price: float
    delta: float
    gamma: float
    theta: float
    vega: float
    rho: float


def _validate_option_type(option_type: str) -> OptionType:
    option = option_type.lower()
    if option not in {"call", "put"}:
        raise ValueError("option_type must be either 'call' or 'put'")
    return option  # type: ignore[return-value]


def _to_quantlib_option_type(option_type: OptionType) -> ql.Option.Type:
    return ql.Option.Call if option_type == "call" else ql.Option.Put


def _build_dates(time_to_expiry: float, valuation_date: date | None) -> tuple[ql.Date, ql.Date]:
    if time_to_expiry <= 0.0:
        raise ValueError("time_to_expiry must be positive")

    valuation = valuation_date or date.today()
    maturity = valuation + timedelta(days=max(1, round(time_to_expiry * 365)))
    return (
        ql.Date(valuation.day, valuation.month, valuation.year),
        ql.Date(maturity.day, maturity.month, maturity.year),
    )


def _american_price(
    params: AmericanOptionInputs,
    *,
    spot: float | None = None,
    risk_free_rate: float | None = None,
    volatility: float | None = None,
    time_to_expiry: float | None = None,
) -> float:
    option_type = _validate_option_type(params.option_type)
    used_spot = params.spot if spot is None else spot
    used_rate = params.risk_free_rate if risk_free_rate is None else risk_free_rate
    used_vol = params.volatility if volatility is None else volatility
    used_time = params.time_to_expiry if time_to_expiry is None else time_to_expiry

    if used_spot <= 0 or used_vol <= 0 or used_time <= 0:
        raise ValueError("spot, volatility, and time_to_expiry must be positive")
    if params.steps < 10:
        raise ValueError("steps must be at least 10 for binomial pricing")

    valuation_ql, maturity_ql = _build_dates(used_time, params.valuation_date)
    ql.Settings.instance().evaluationDate = valuation_ql

    day_count = ql.Actual365Fixed()
    calendar = ql.NullCalendar()

    spot_handle = ql.QuoteHandle(ql.SimpleQuote(used_spot))
    dividend_curve = ql.YieldTermStructureHandle(ql.FlatForward(valuation_ql, params.dividend_yield, day_count))
    risk_free_curve = ql.YieldTermStructureHandle(ql.FlatForward(valuation_ql, used_rate, day_count))
    volatility_curve = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(valuation_ql, calendar, used_vol, day_count))

    process = ql.BlackScholesMertonProcess(
        spot_handle,
        dividend_curve,
        risk_free_curve,
        volatility_curve,
    )

    payoff = ql.PlainVanillaPayoff(_to_quantlib_option_type(option_type), params.strike)
    exercise = ql.AmericanExercise(valuation_ql, maturity_ql)
    option = ql.VanillaOption(payoff, exercise)
    option.setPricingEngine(ql.BinomialVanillaEngine(process, params.tree, params.steps))
    return float(option.NPV())


def american_option_price_and_greeks(
    params: AmericanOptionInputs,
) -> AmericanPricingResult:
    """Price an American option and estimate Greeks with finite differences.

    Notes:
        - The price uses a QuantLib binomial tree engine.
        - Greeks are finite-difference approximations because closed-form
          American Greeks are not available in general.
        - Theta is annualized and approximated by rolling one calendar day forward.

    Args:
        params: American option inputs.

    Returns:
        A price and finite-difference Greeks.

    Example:
        >>> params = AmericanOptionInputs(
        ...     spot=460.0,
        ...     strike=400.0,
        ...     time_to_expiry=0.5,
        ...     risk_free_rate=0.045,
        ...     volatility=0.16,
        ...     option_type="put",
        ... )
        >>> result = american_option_price_and_greeks(params)
        >>> result.price > 0
        True
    """

    base_price = _american_price(params)

    spot_bump = max(0.01, params.spot * 0.01)
    vol_bump = 0.01
    rate_bump = 0.0001
    dt = 1.0 / 365.0

    price_up = _american_price(params, spot=params.spot + spot_bump)
    price_down = _american_price(params, spot=max(1e-8, params.spot - spot_bump))
    delta = (price_up - price_down) / (2.0 * spot_bump)
    gamma = (price_up - 2.0 * base_price + price_down) / (spot_bump**2)

    vega_up = _american_price(params, volatility=params.volatility + vol_bump)
    vega_down = _american_price(params, volatility=max(1e-6, params.volatility - vol_bump))
    vega = (vega_up - vega_down) / (2.0 * vol_bump)

    rho_up = _american_price(params, risk_free_rate=params.risk_free_rate + rate_bump)
    rho_down = _american_price(params, risk_free_rate=params.risk_free_rate - rate_bump)
    rho = (rho_up - rho_down) / (2.0 * rate_bump)

    if params.time_to_expiry <= dt:
        theta = 0.0
    else:
        shorter_price = _american_price(params, time_to_expiry=params.time_to_expiry - dt)
        theta = (shorter_price - base_price) / dt

    return AmericanPricingResult(
        price=base_price,
        delta=delta,
        gamma=gamma,
        theta=theta,
        vega=vega,
        rho=rho,
    )