OpenFin/backend/services/strategy_engine.py

"""
Generic N-leg (1-4) option strategy pricer: entry cost with real broker spread,
scenario repricing at a future horizon on a shocked vol surface, payoff curves,
greeks, and bounded-risk / non-directional checks.

A "leg" dict: {expiry_date, days_to_expiry, strike, option_type ("call"/"put"),
               position ("long"/"short"), quantity}
"""
import math
from typing import Any, Dict, List, Optional

import numpy as np
from scipy.optimize import minimize_scalar

from services.options_pricer import black_scholes
from services.option_chain import find_quote
from services.vol_surface import Surface, ScenarioSurface

DEFAULT_SPREAD_PCT = 0.05  # fallback relative bid/ask spread when no live quote is found
DEFAULT_CONTRACT_SIZE = 100_000  # notional per 1 contract/lot (e.g. a standard FX lot); "quantity" on a leg is the number of these


def to_native(obj: Any) -> Any:
    """Recursively converts numpy scalars (bool_, int64, float64, ...) to native Python
    types. Comparisons/aggregations over numpy-typed inputs (e.g. bid/ask sourced from a
    DB row that came back as a numpy type, or scipy's norm.cdf) can leave a stray
    numpy.bool_/numpy.float64 buried in a nested result — FastAPI's default JSON encoder
    doesn't know those types and fails ('X object is not iterable', then a secondary
    'vars() argument must have __dict__ attribute' from its own fallback). Applied once
    at the API boundary (price_combo/payoff_curves/optimizer results) rather than chasing
    the exact field through the whole pricing pipeline."""
    if isinstance(obj, dict):
        return {k: to_native(v) for k, v in obj.items()}
    if isinstance(obj, (list, tuple)):
        return [to_native(v) for v in obj]
    if isinstance(obj, np.generic):
        return obj.item()
    return obj


def _sign(leg: Dict[str, Any]) -> int:
    return 1 if leg.get("position", "long") == "long" else -1


def _intrinsic(S: float, K: float, option_type: str) -> float:
    return max(0.0, S - K) if option_type == "call" else max(0.0, K - S)


def _quote_spread_pct(quote: Optional[Dict[str, Any]]) -> float:
    if not quote or quote["bid"] <= 0 or quote["ask"] <= 0:
        return DEFAULT_SPREAD_PCT
    mid = (quote["bid"] + quote["ask"]) / 2
    if mid <= 0:
        return DEFAULT_SPREAD_PCT
    return (quote["ask"] - quote["bid"]) / mid


def entry_price(leg: Dict[str, Any], chain_slice: Dict[str, Any], surface_now: Surface, r: float) -> Dict[str, float]:
    """Real execution price (crossing the spread) + theoretical mid, for one leg today.
    A "stock" leg (Covered Call/Protective Put/Collar's underlying position, not an
    option) has no strike/vol at all — its price IS the spot, no spread modeled (the
    underlying's own spread is typically far tighter than any option on it, and this
    tool doesn't have a live quote for it the way find_quote does for options)."""
    if leg["option_type"] == "stock":
        S = chain_slice["spot"]
        return {"exec_price": S, "mid": S, "spread_pct": 0.0}
    quote = find_quote(chain_slice, leg["expiry_date"], leg["strike"], leg["option_type"])
    T = max(leg["days_to_expiry"], 0.001) / 365
    sigma = surface_now.iv_at(leg["strike"], leg["days_to_expiry"])
    theo_mid = black_scholes(chain_slice["spot"], leg["strike"], T, r, sigma, leg["option_type"])["price"]

    if quote and quote["bid"] > 0 and quote["ask"] > 0:
        exec_price = quote["ask"] if leg.get("position", "long") == "long" else quote["bid"]
        mid = quote["mid"] or theo_mid
    else:
        spread = theo_mid * DEFAULT_SPREAD_PCT
        exec_price = theo_mid + spread / 2 if leg.get("position", "long") == "long" else max(0.0, theo_mid - spread / 2)
        mid = theo_mid

    return {"exec_price": exec_price, "mid": mid, "spread_pct": _quote_spread_pct(quote)}


def value_at(
    legs: List[Dict[str, Any]],
    S: float,
    eval_days_from_now: float,
    surface: Any,
    r: float,
    contract_size: float = DEFAULT_CONTRACT_SIZE,
) -> float:
    """Signed portfolio value (BS reprice for unexpired legs, intrinsic for expired ones).
    check_bounded_risk calls this ~700 times per candidate it evaluates — second-order
    Greeks are never read here, so they're skipped (include_second_order=False) rather than
    computed and discarded on every one of those calls."""
    total = 0.0
    for leg in legs:
        qty = leg.get("quantity", 1)
        sign = _sign(leg)
        if leg["option_type"] == "stock":
            price = S  # a share/lot of the underlying is worth exactly the spot, always
        else:
            remaining = leg["days_to_expiry"] - eval_days_from_now
            if remaining <= 0:
                price = _intrinsic(S, leg["strike"], leg["option_type"])
            else:
                sigma = surface.iv_at(leg["strike"], remaining)
                price = black_scholes(S, leg["strike"], remaining / 365, r, sigma, leg["option_type"], include_second_order=False)["price"]
        total += sign * price * qty * contract_size
    return total


def greeks_at(legs: List[Dict[str, Any]], S: float, eval_days_from_now: float, surface: Any, r: float) -> Dict[str, float]:
    net = {
        "delta": 0.0, "gamma": 0.0, "theta": 0.0, "vega": 0.0, "rho": 0.0,
        "vanna": 0.0, "charm": 0.0, "vomma": 0.0, "veta": 0.0, "speed": 0.0, "color": 0.0, "zomma": 0.0,
    }
    for leg in legs:
        qty = leg.get("quantity", 1)
        sign = _sign(leg)
        if leg["option_type"] == "stock":
            # d(spot)/d(spot) = 1, and every other Greek (gamma, theta, vega, rho, the
            # second-order ones) is exactly zero for a position in the underlying itself.
            net["delta"] += sign * qty
            continue
        remaining = max(leg["days_to_expiry"] - eval_days_from_now, 0.001)
        sigma = surface.iv_at(leg["strike"], remaining)
        g = black_scholes(S, leg["strike"], remaining / 365, r, sigma, leg["option_type"])
        for k in net:
            net[k] += g[k] * qty * sign
    return {k: round(v, 6) for k, v in net.items()}


def vanna_simulation(
    legs: List[Dict[str, Any]], S: float, eval_days_from_now: float, surface: Any, r: float,
    spot_shock_pct: float = -5.0, iv_shock_pts: float = 8.0,
) -> Dict[str, float]:
    """A concrete joint spot+IV shock reprice — "if spot drops 5% and IV jumps 8pts, what
    actually happens to my net delta" — rather than a bare "vanna is positive/negative"
    label. Uses a real Black-Scholes reprice (not the linear vanna approximation) so it's
    accurate for shocks this large, matching the same "show a simulation, not a sign"
    principle the payoff diagram already uses elsewhere in Strategy Builder."""
    delta_before = greeks_at(legs, S, eval_days_from_now, surface, r)["delta"]

    class _ShockedSurface:
        def iv_at(self, strike: float, days: float) -> float:
            return max(0.01, surface.iv_at(strike, days) + iv_shock_pts / 100.0)

    S_shocked = S * (1 + spot_shock_pct / 100.0)
    delta_after = greeks_at(legs, S_shocked, eval_days_from_now, _ShockedSurface(), r)["delta"]

    return {
        "spot_shock_pct": spot_shock_pct,
        "iv_shock_pts": iv_shock_pts,
        "delta_before": delta_before,
        "delta_after": delta_after,
        "delta_change": round(delta_after - delta_before, 6),
    }


def price_combo(
    legs: List[Dict[str, Any]],
    chain_slice: Dict[str, Any],
    surface_now: Surface,
    surface_scenario: ScenarioSurface,
    horizon_days: int,
    r: float = 0.05,
    contract_size: float = DEFAULT_CONTRACT_SIZE,
    precise: bool = True,
) -> Dict[str, Any]:
    spot_now = chain_slice["spot"]
    spot_scenario = surface_scenario.spot

    entry_ref = 0.0
    entry_ref_mid = 0.0
    for leg in legs:
        ep = entry_price(leg, chain_slice, surface_now, r)
        sign = _sign(leg)
        qty = leg.get("quantity", 1)
        entry_ref += sign * ep["exec_price"] * qty * contract_size
        entry_ref_mid += sign * ep["mid"] * qty * contract_size

    # Scenario exit: apply each leg's own bid/ask spread (est. from entry quote) to the
    # theoretical scenario value, since we don't have a live quote for the future date.
    scenario_mid = 0.0
    scenario_exec = 0.0
    for leg in legs:
        qty = leg.get("quantity", 1)
        sign = _sign(leg)
        if leg["option_type"] == "stock":
            theo = exec_price = spot_scenario  # no spread modeled for the underlying itself
        else:
            remaining = max(leg["days_to_expiry"] - horizon_days, 0.001)
            sigma = surface_scenario.iv_at(leg["strike"], remaining)
            theo = black_scholes(spot_scenario, leg["strike"], remaining / 365, r, sigma, leg["option_type"])["price"]
            quote = find_quote(chain_slice, leg["expiry_date"], leg["strike"], leg["option_type"])
            spread_pct = _quote_spread_pct(quote)
            exec_price = theo * (1 - spread_pct / 2) if leg.get("position", "long") == "long" else theo * (1 + spread_pct / 2)
        scenario_mid += sign * theo * qty * contract_size
        scenario_exec += sign * exec_price * qty * contract_size

    net_pnl = scenario_exec - entry_ref
    broker_cost = (entry_ref - entry_ref_mid) + (scenario_mid - scenario_exec)

    # Uses surface_scenario (not surface_now): for a single-expiry combo (condor,
    # butterfly, straddle...) all legs are simultaneously intrinsic at eval_days, so the
    # surface is never actually consulted there and this changes nothing. For a
    # calendar/diagonal spread, the far leg is still alive at the near leg's expiry and
    # DOES need a vol assumption to be priced — using surface_now there silently mixed
    # "today's vol" into a number sitting next to net_pnl (which uses the scenario's
    # shocked vol), producing a max_gain that could be below net_pnl. Pricing both with
    # the same scenario vol view keeps them consistent.
    bounded = check_bounded_risk(legs, entry_ref, surface_scenario, spot_now, r, contract_size, precise)
    delta_now = greeks_at(legs, spot_now, 0, surface_now, r)["delta"]
    delta_scenario = greeks_at(legs, spot_scenario, horizon_days, surface_scenario, r)["delta"]

    return to_native({
        "entry_cost": round(entry_ref, 2),
        "entry_cost_mid": round(entry_ref_mid, 2),
        "scenario_value": round(scenario_exec, 2),
        "scenario_value_mid": round(scenario_mid, 2),
        "net_pnl": round(net_pnl, 2),
        "broker_spread_cost": round(broker_cost, 2),
        "max_gain": bounded["max_gain"],
        "max_loss": bounded["max_loss"],
        "bounded_risk": bounded["bounded"],
        "greeks_now": greeks_at(legs, spot_now, 0, surface_now, r),
        "greeks_scenario": greeks_at(legs, spot_scenario, horizon_days, surface_scenario, r),
        "net_delta_now": delta_now,
        "net_delta_scenario": delta_scenario,
        # Skipped during the optimizer's bulk scan (precise=False, hundreds of candidates
        # per request) — only computed for the single position actually loaded/priced.
        "vanna_simulation": vanna_simulation(legs, spot_now, 0, surface_now, r) if precise else None,
    })


def check_bounded_risk(
    legs: List[Dict[str, Any]], entry_ref: float, surface: Any, spot: float, r: float = 0.05,
    contract_size: float = DEFAULT_CONTRACT_SIZE, precise: bool = True,
) -> Dict[str, Any]:
    """
    Scan a wide log-spaced spot range at expiry and inspect both tails independently for LOSS
    vs GAIN direction. "Bounded risk" only requires the loss side to be capped — a long
    straddle (loss capped at the premium, gain uncapped) is a textbook risk-bounded,
    non-directional trade and must not be excluded just because its gain is open-ended.

    A tail is loss-bounded if moving further to that extreme does not make the P&L any worse
    than a point already well into that tail; symmetrically for gain-bounded.

    "At expiry" means the earliest expiry among the legs. For a single-expiry combo that's
    every leg simultaneously (pure intrinsic, no vol assumption involved at all). For a
    calendar/diagonal spread it's only the near leg — the far leg is still alive and needs
    `surface` to be priced, so max_gain/max_loss there is only as good as that vol input.

    `precise=False` skips only the 1-D refinement below (the expensive part — two
    minimize_scalar calls). The dense near-strike grid stays even in the fast path: a long
    straddle's max loss (or a butterfly's max gain) sits exactly AT a strike too, same as a
    calendar spread's peak, and the wide log-spaced tail grid alone is coarse enough near
    the money to miss it by an order of magnitude, not just a rounding difference (confirmed:
    on a real long straddle, tail-grid-only reported max_loss=-117 when the true V-bottom at
    the strike was -1140 — nearly 10x off, and wrong enough to silently let a candidate past
    a max_loss_cap filter it should have failed). The optimizer's bulk scan (hundreds of
    candidates) uses this fast path for speed; the single-candidate /price call refines
    further with minimize_scalar for pixel-perfect precision.
    """
    eval_days = min(l["days_to_expiry"] for l in legs)

    def f(s: float) -> float:
        return value_at(legs, s, eval_days, surface, r, contract_size) - entry_ref

    # Wide, log-spaced tail grid — used to detect whether the payoff flattens out (bounded)
    # toward either extreme, or keeps moving further away.
    tail_grid = np.geomspace(spot * 0.05, spot * 20, 300)
    tail_values = [f(float(p)) for p in tail_grid]

    tail_n = max(3, len(tail_values) // 30)
    tol = max(abs(entry_ref), 1.0) * 0.01
    lo_edge, lo_in = tail_values[0], tail_values[tail_n]
    hi_edge, hi_in = tail_values[-1], tail_values[-1 - tail_n]

    loss_bounded = (lo_edge >= lo_in - tol) and (hi_edge >= hi_in - tol)
    gain_bounded = (lo_edge <= lo_in + tol) and (hi_edge <= hi_in + tol)

    # Dense linear sweep across the legs' own strikes — needed even in the fast path (see
    # docstring above), only the final 1-D refinement is reserved for precise=True.
    strikes = [l["strike"] for l in legs]
    lo_k, hi_k = min(strikes) * 0.7, max(strikes) * 1.3
    near_grid = np.linspace(max(lo_k, spot * 0.05), min(hi_k, spot * 20), 400)
    near_values = [f(float(p)) for p in near_grid]
    combined_values = tail_values + near_values

    if not precise:
        return {
            "bounded": loss_bounded,
            "max_loss": round(min(combined_values), 2) if loss_bounded else None,
            "max_gain": round(max(combined_values), 2) if gain_bounded else None,
        }

    # Any FIXED grid — however dense — is a different finite sampling of the same
    # continuous curve than whatever grid a chart or another caller uses, so two "close but
    # not identical" readings of the same peak are pretty much guaranteed (confirmed: this
    # grid vs the payoff chart's own grid gave two different peak heights for the same
    # trade). Refine with a bounded 1-D optimizer around the grid's own best point instead.
    grid_all = np.concatenate([tail_grid, near_grid])
    values_all = np.concatenate([tail_values, near_values])
    order = np.argsort(grid_all)
    grid_sorted, values_sorted = grid_all[order], values_all[order]

    # ...then refine with a bounded 1-D optimizer in a NARROW bracket around that grid
    # point (the payoff at a fixed date/vol is smooth in spot, built from Black-Scholes).
    # A wide bracket is actively harmful here: tested directly on a real calendar spread,
    # minimize_scalar given a wide bound (±50% of the strike) converged on a flat false
    # optimum far from the true spike (-65 instead of +1170) because Brent's method isn't
    # guaranteed to find the global optimum over a non-unimodal interval. Anchoring the
    # bracket tightly around the grid's own best point keeps the search unimodal, and
    # taking max()/min() against the grid value means refinement can never do worse.
    def refine(is_max: bool) -> float:
        idx = int(np.argmax(values_sorted)) if is_max else int(np.argmin(values_sorted))
        grid_value = float(values_sorted[idx])
        step = grid_sorted[min(idx + 1, len(grid_sorted) - 1)] - grid_sorted[max(idx - 1, 0)]
        step = float(step) if step > 0 else spot * 0.001
        lo_b, hi_b = grid_sorted[idx] - 3 * step, grid_sorted[idx] + 3 * step
        if lo_b >= hi_b:
            return grid_value
        res = minimize_scalar((lambda s: -f(s)) if is_max else f, bounds=(float(lo_b), float(hi_b)), method="bounded")
        refined = -res.fun if is_max else res.fun
        return max(refined, grid_value) if is_max else min(refined, grid_value)

    return {
        "bounded": loss_bounded,
        "max_loss": round(refine(False), 2) if loss_bounded else None,
        "max_gain": round(refine(True), 2) if gain_bounded else None,
    }


def payoff_curve_expiry(legs: List[Dict[str, Any]], surface_now: Surface, spot: float, n: int = 100, contract_size: float = DEFAULT_CONTRACT_SIZE) -> List[Dict[str, float]]:
    eval_days = min(l["days_to_expiry"] for l in legs)
    prices = np.linspace(spot * 0.5, spot * 1.5, n)
    return [
        {"underlying": round(float(p), 2), "pnl": round(float(value_at(legs, float(p), eval_days, surface_now, 0.05, contract_size)), 2)}
        for p in prices
    ]


def expected_pnl_scenario(
    legs: List[Dict[str, Any]],
    surface_scenario: ScenarioSurface,
    horizon_days: int,
    r: float,
    entry_ref: float,
    n: int = 200,
    contract_size: float = DEFAULT_CONTRACT_SIZE,
) -> float:
    """
    Probability-weighted expected P&L at the scenario date: integrates the payoff over a
    risk-neutral lognormal density for the underlying, centered on the scenario spot with a
    variance driven by the scenario surface's ATM IV over the residual time to the nearest
    leg's expiry (the remaining uncertainty once the scenario date is reached).
    """
    spot_scenario = surface_scenario.spot
    remaining = max(min(l["days_to_expiry"] for l in legs) - horizon_days, 1)
    sigma = surface_scenario.iv_at(spot_scenario, remaining)
    T = remaining / 365
    sd = sigma * math.sqrt(T)
    if sd <= 1e-6:
        return value_at(legs, spot_scenario, horizon_days, surface_scenario, r, contract_size) - entry_ref

    mu = math.log(spot_scenario) + (r - 0.5 * sigma ** 2) * T
    grid = np.geomspace(spot_scenario * 0.15, spot_scenario * 4, n)
    log_grid = np.log(grid)
    density = np.exp(-0.5 * ((log_grid - mu) / sd) ** 2) / (grid * sd * math.sqrt(2 * math.pi))
    pnl = np.array([value_at(legs, float(s), horizon_days, surface_scenario, r, contract_size) - entry_ref for s in grid])

    numerator = np.trapz(pnl * density, grid)
    denominator = np.trapz(density, grid)
    return float(numerator / denominator) if denominator > 1e-12 else float(pnl.mean())


def time_decay_slices(
    legs: List[Dict[str, Any]], prices: np.ndarray, surface: Any, eval_days_expiry: float, r: float,
    entry_ref: float, contract_size: float = DEFAULT_CONTRACT_SIZE, n_slices: int = 4,
) -> List[Dict[str, Any]]:
    """Payoff curve at n_slices evenly-spaced elapsed-day checkpoints between today (0) and
    the nearest leg's expiry — the "T+0/T+10/T+20..." view that shows how the curve morphs
    from today's time-value-laden shape into the kinked expiry payoff, instead of only the
    two endpoints at_expiry/at_scenario give. Uses the same scenario vol view as those two
    curves (see payoff_curves' own comment) so the only thing that varies between slices is
    time decay, not the vol assumption."""
    day_points = np.linspace(0, eval_days_expiry, n_slices)
    slices = []
    for d in day_points:
        d = float(d)
        label = "Aujourd'hui" if d < 0.5 else ("Échéance" if d >= eval_days_expiry - 0.5 else f"J+{round(d)}")
        points = [
            {"underlying": round(float(p), 4), "pnl": round(float(value_at(legs, float(p), d, surface, r, contract_size) - entry_ref), 2)}
            for p in prices
        ]
        slices.append({"days_from_now": round(d, 1), "label": label, "points": points})
    return slices


def payoff_heatmap(
    legs: List[Dict[str, Any]], surface: Any, eval_days_expiry: float, r: float, spot: float,
    entry_ref: float, contract_size: float = DEFAULT_CONTRACT_SIZE, n_prices: int = 9, n_days: int = 7,
) -> Dict[str, Any]:
    """Price x days-to-expiry grid of P&L — rows are elapsed-day checkpoints from today down
    to expiry (top-to-bottom reading matches watching the position age), columns are
    underlying prices zoomed closer to spot than the line chart (a heatmap only reads well
    over the range where the color actually varies)."""
    lo, hi = spot * 0.85, spot * 1.15
    price_points = np.linspace(lo, hi, n_prices)
    day_points = np.linspace(0, eval_days_expiry, n_days)
    rows = [
        {
            "days_from_now": round(float(d), 1),
            "pnl": [round(float(value_at(legs, float(p), float(d), surface, r, contract_size) - entry_ref), 2) for p in price_points],
        }
        for d in day_points
    ]
    return {"prices": [round(float(p), 4) for p in price_points], "rows": rows}


def payoff_curves(
    legs: List[Dict[str, Any]],
    chain_slice: Dict[str, Any],
    surface_now: Surface,
    surface_scenario: ScenarioSurface,
    horizon_days: int,
    r: float = 0.05,
    n: int = 100,
    contract_size: float = DEFAULT_CONTRACT_SIZE,
) -> Dict[str, Any]:
    spot = chain_slice["spot"]
    priced = price_combo(legs, chain_slice, surface_now, surface_scenario, horizon_days, r, contract_size)
    entry_ref = priced["entry_cost"]

    lo, hi = spot * 0.6, spot * 1.4
    # A calendar/ratio spread's payoff can spike sharply right at a strike — a plain
    # uniform sweep over the full 0.6x-1.4x range (n points) can straddle right over that
    # peak without ever sampling it (same issue fixed in check_bounded_risk). Blend a
    # coarse baseline (overall shape) with a dense window around the legs' own strikes.
    # The exact spot/scenario spot are forced in as sample points too — otherwise hovering
    # right on the "Spot"/"Scénario" reference line reads the nearest grid point, which can
    # sit meaningfully off the true value on a peak this steep, and disagree with the
    # entry_cost/net_pnl tiles (computed at the exact spot, not a grid sample).
    baseline = np.linspace(lo, hi, n)
    strikes = [l["strike"] for l in legs]
    lo_k, hi_k = max(min(strikes) * 0.9, lo), min(max(strikes) * 1.1, hi)
    near_strikes = np.linspace(lo_k, hi_k, n * 3)
    exact_points = np.array([spot, surface_scenario.spot] + strikes)
    prices = np.unique(np.concatenate([baseline, near_strikes, exact_points]))
    prices.sort()
    eval_days_expiry = min(l["days_to_expiry"] for l in legs)

    # Both curves share the scenario's volatility view (surface_scenario) — only the
    # evaluation DATE differs: "à échéance" prices the day the near leg expires (matching
    # the Max gain/Max perte tile, itself computed with surface_scenario for the same
    # reason — see check_bounded_risk), "à J+8" prices your chosen scenario horizon. Using
    # surface_now here instead would silently mix in today's un-shocked vol, producing a
    # curve whose peak doesn't match the Max gain tile right next to it.
    at_expiry = [
        {"underlying": round(float(p), 4), "pnl": round(float(value_at(legs, float(p), eval_days_expiry, surface_scenario, r, contract_size) - entry_ref), 2)}
        for p in prices
    ]
    at_scenario = [
        {"underlying": round(float(p), 4), "pnl": round(float(value_at(legs, float(p), horizon_days, surface_scenario, r, contract_size) - entry_ref), 2)}
        for p in prices
    ]

    # Coarser price grid than the two headline curves above — this trades some precision
    # for keeping a single /price request's added cost bounded (n_slices/heatmap cells x
    # their own price points, on top of the ~1400 value_at calls at_expiry/at_scenario
    # above already need).
    slice_prices = np.linspace(lo, hi, 80)
    time_slices = time_decay_slices(legs, slice_prices, surface_scenario, eval_days_expiry, r, entry_ref, contract_size)
    heatmap = payoff_heatmap(legs, surface_scenario, eval_days_expiry, r, spot, entry_ref, contract_size)

    return {"at_expiry": at_expiry, "at_scenario": at_scenario, "time_slices": time_slices, "heatmap": heatmap, **priced}