market_data/analytics.py

from __future__ import annotations

import math
from statistics import NormalDist

import pandas as pd

from .schemas import OptionChain


NORMAL = NormalDist()


def realized_volatility(
    prices: pd.Series,
    windows: tuple[int, ...] = (5, 10, 20, 30, 60),
    trading_days: int = 252,
) -> dict[str, float | None]:
    close = prices.dropna().astype(float)
    returns = close.pct_change().dropna()
    output: dict[str, float | None] = {}

    for window in windows:
        key = f"{window}d"
        if len(returns) < window:
            output[key] = None
            continue
        output[key] = float(returns.tail(window).std(ddof=1) * math.sqrt(trading_days))

    return output


def _norm_pdf(value: float) -> float:
    return math.exp(-0.5 * value * value) / math.sqrt(2 * math.pi)


def black_scholes_greeks(
    spot: float,
    strike: float,
    time_to_expiry: float,
    volatility: float,
    risk_free_rate: float = 0.0,
    dividend_yield: float = 0.0,
    option_type: str = "call",
) -> dict[str, float | None]:
    if spot <= 0 or strike <= 0 or time_to_expiry <= 0 or volatility <= 0:
        return {
            "delta": None,
            "gamma": None,
            "vega": None,
            "theta": None,
            "rho": None,
        }

    sqrt_t = math.sqrt(time_to_expiry)
    d1 = (
        math.log(spot / strike)
        + (risk_free_rate - dividend_yield + 0.5 * volatility * volatility) * time_to_expiry
    ) / (volatility * sqrt_t)
    d2 = d1 - volatility * sqrt_t
    discount_dividend = math.exp(-dividend_yield * time_to_expiry)
    discount_rate = math.exp(-risk_free_rate * time_to_expiry)
    option_type = option_type.lower()

    if option_type == "put":
        delta = discount_dividend * (NORMAL.cdf(d1) - 1)
        theta = (
            -spot * discount_dividend * _norm_pdf(d1) * volatility / (2 * sqrt_t)
            + dividend_yield * spot * discount_dividend * NORMAL.cdf(-d1)
            - risk_free_rate * strike * discount_rate * NORMAL.cdf(-d2)
        ) / 365
        rho = -strike * time_to_expiry * discount_rate * NORMAL.cdf(-d2) / 100
    else:
        delta = discount_dividend * NORMAL.cdf(d1)
        theta = (
            -spot * discount_dividend * _norm_pdf(d1) * volatility / (2 * sqrt_t)
            - dividend_yield * spot * discount_dividend * NORMAL.cdf(d1)
            + risk_free_rate * strike * discount_rate * NORMAL.cdf(d2)
        ) / 365
        rho = strike * time_to_expiry * discount_rate * NORMAL.cdf(d2) / 100

    gamma = discount_dividend * _norm_pdf(d1) / (spot * volatility * sqrt_t)
    vega = spot * discount_dividend * _norm_pdf(d1) * sqrt_t / 100

    return {
        "delta": float(delta),
        "gamma": float(gamma),
        "vega": float(vega),
        "theta": float(theta),
        "rho": float(rho),
    }


def nearest_atm_iv(chain: OptionChain) -> float | None:
    if chain.underlying_price is None:
        return None
    contracts = chain.calls + chain.puts
    valid = [
        contract
        for contract in contracts
        if contract.implied_volatility is not None and contract.implied_volatility > 0
    ]
    if not valid:
        return None
    nearest = min(valid, key=lambda contract: abs(contract.strike - chain.underlying_price))
    return nearest.implied_volatility


def simple_skew(chain: OptionChain) -> float | None:
    if chain.underlying_price is None:
        return None
    otm_puts = [
        contract
        for contract in chain.puts
        if contract.strike < chain.underlying_price and contract.implied_volatility
    ]
    otm_calls = [
        contract
        for contract in chain.calls
        if contract.strike > chain.underlying_price and contract.implied_volatility
    ]
    if not otm_puts or not otm_calls:
        return None
    put = max(otm_puts, key=lambda contract: contract.strike)
    call = min(otm_calls, key=lambda contract: contract.strike)
    return float((put.implied_volatility or 0) - (call.implied_volatility or 0))


def summarize_option_chain(chain: OptionChain, realized_vol_20d: float | None = None) -> dict:
    atm_iv = nearest_atm_iv(chain)
    return {
        "symbol": chain.symbol,
        "expiration": chain.expiration,
        "underlying_price": chain.underlying_price,
        "atm_iv": atm_iv,
        "iv_rv_spread_20d": (
            float(atm_iv - realized_vol_20d)
            if atm_iv is not None and realized_vol_20d is not None
            else None
        ),
        "skew_put_minus_call": simple_skew(chain),
        "call_count": len(chain.calls),
        "put_count": len(chain.puts),
    }


def rank_current_iv_against_rv(
    current_iv: float | None,
    realized_vols: dict[str, float | None],
) -> float | None:
    if current_iv is None:
        return None
    rv_values = [value for value in realized_vols.values() if value is not None]
    if len(rv_values) < 2:
        return None
    low = min(rv_values)
    high = max(rv_values)
    if high <= low:
        return None
    return max(0.0, min(1.0, (current_iv - low) / (high - low)))


def classify_volatility_regime(
    current_iv: float | None,
    realized_vol_20d: float | None,
    term_structure_slope: float | None,
    skew: float | None,
) -> dict:
    if current_iv is None or realized_vol_20d is None:
        return {
            "regime": "unknown",
            "vol_signal": "insufficient_iv_or_rv",
            "confidence": "low",
            "notes": ["Need both option implied volatility and realized volatility."],
        }

    iv_rv_spread = current_iv - realized_vol_20d
    notes = []
    if iv_rv_spread > 0.08:
        regime = "high_implied_vol_premium"
        vol_signal = "short_vol_candidate"
        notes.append("Current ATM IV is materially above 20D realized volatility.")
    elif iv_rv_spread < -0.04:
        regime = "low_implied_vol_discount"
        vol_signal = "long_vol_candidate"
        notes.append("Current ATM IV is below 20D realized volatility.")
    else:
        regime = "balanced_iv_vs_rv"
        vol_signal = "neutral_vol"
        notes.append("Current ATM IV is close to 20D realized volatility.")

    if term_structure_slope is not None:
        if term_structure_slope > 0.04:
            notes.append("Term structure is upward sloping.")
        elif term_structure_slope < -0.04:
            notes.append("Term structure is inverted or front-loaded.")
    if skew is not None and abs(skew) > 0.05:
        notes.append("Put-call skew is elevated in the sampled expiration.")

    confidence = "medium" if len(notes) >= 2 else "low"
    return {
        "regime": regime,
        "vol_signal": vol_signal,
        "confidence": confidence,
        "notes": notes,
    }