impact_audit.py

"""Audit the SIMULATED fill-impact term against executable REAL L5 depth.

For each simulated child timestamp, the module walks the contemporaneous five-level AAPL
book. Kappa is fitted on the first execution window and audited on the second. Orders larger
than cumulative L5 depth are labelled unsupported rather than extrapolated as calibration.
"""

from pathlib import Path
from typing import Any, Iterable, Mapping, Tuple

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


_PRICE_SCALE = 10_000.0


def adjusted_participation(values, over_displayed_multiplier: float):
    """Return the exact piecewise participation transform used by the C# fill model."""

    values = np.asarray(values, dtype="float64")
    if over_displayed_multiplier < 1.0:
        raise ValueError("over_displayed_multiplier must be at least one")
    return np.where(
        values <= 1.0,
        values,
        1.0 + over_displayed_multiplier * (values - 1.0),
    )


def build_depth_walk_sample(
    message_path: Path,
    orderbook_path: Path,
    fills: pd.DataFrame,
    *,
    levels: int,
    tick_size: float,
    configured_kappa: float,
    over_displayed_multiplier: float,
) -> pd.DataFrame:
    """Walk REAL book depth at every SIMULATED fill timestamp."""

    if levels < 2 or tick_size <= 0 or configured_kappa < 0:
        raise ValueError("Depth-audit parameters are invalid")
    required = {
        "fill_id",
        "parent_id",
        "data_label",
        "side",
        "size",
        "timestamp",
        "fill_price",
        "mid_at_fill",
    }
    missing = required.difference(fills.columns)
    if missing:
        raise ValueError(f"Fills missing columns: {sorted(missing)}")
    if set(fills["data_label"].unique()) != {"SIMULATED"}:
        raise ValueError("Depth audit accepts SIMULATED fills only")

    message_times = pd.read_csv(message_path, header=None, usecols=[0]).iloc[:, 0]
    columns = _book_columns(levels)
    books = pd.read_csv(
        orderbook_path,
        header=None,
        usecols=range(4 * levels),
        names=columns,
    )
    if len(message_times) != len(books):
        raise ValueError("L5 message and order-book row counts differ")
    timestamps = message_times.to_numpy(dtype="float64")
    if np.any(np.diff(timestamps) < 0):
        raise ValueError("L5 messages must be chronological")

    result = fills.sort_values("timestamp").reset_index(drop=True).copy()
    positions = np.searchsorted(
        timestamps, result["timestamp"].to_numpy(dtype="float64"), side="right"
    ) - 1
    if (positions < 0).any():
        raise ValueError("A fill precedes the first L5 book state")
    states = books.iloc[positions].reset_index(drop=True).astype("float64")
    for level in range(1, levels + 1):
        states[f"ask_price_{level}"] /= _PRICE_SCALE
        states[f"bid_price_{level}"] /= _PRICE_SCALE

    rows = []
    for fill, state in zip(result.itertuples(index=False), states.itertuples(index=False)):
        side = int(fill.side)
        prices = np.array(
            [getattr(state, f"ask_price_{level}") for level in range(1, levels + 1)]
            if side == 1
            else [getattr(state, f"bid_price_{level}") for level in range(1, levels + 1)]
        )
        sizes = np.array(
            [getattr(state, f"ask_size_{level}") for level in range(1, levels + 1)]
            if side == 1
            else [getattr(state, f"bid_size_{level}") for level in range(1, levels + 1)],
            dtype="int64",
        )
        walk = _walk_book(prices, sizes, int(fill.size))
        best_ask = float(state.ask_price_1)
        best_bid = float(state.bid_price_1)
        mid = (best_ask + best_bid) / 2.0
        touch = best_ask if side == 1 else best_bid
        level_one_size = int(state.ask_size_1 if side == 1 else state.bid_size_1)
        participation = int(fill.size) / level_one_size
        transformed = float(
            adjusted_participation([participation], over_displayed_multiplier)[0]
        )
        tick_bps = 10_000.0 * tick_size / mid
        configured_beyond_ticks = configured_kappa * transformed
        walk_beyond_ticks = (
            side * (walk["vwap"] - touch) / tick_size if walk["supported"] else np.nan
        )
        walk_immediate_cost_bps = (
            10_000.0 * side * (walk["vwap"] - mid) / mid
            if walk["supported"]
            else np.nan
        )
        rows.append(
            {
                "market_data_label": "REAL",
                "data_label": "SIMULATED",
                "fill_id": fill.fill_id,
                "parent_id": fill.parent_id,
                "timestamp": float(fill.timestamp),
                "side": side,
                "size": int(fill.size),
                "mid_l5": mid,
                "touch_l5": touch,
                "level_one_size": level_one_size,
                "level_one_participation": participation,
                "adjusted_participation": transformed,
                "cumulative_l5_depth": int(sizes.sum()),
                "depth_coverage_fraction": min(1.0, float(sizes.sum() / int(fill.size))),
                "full_l5_supported": bool(walk["supported"]),
                "walk_vwap_l5": walk["vwap"],
                "walk_levels_used": walk["levels_used"],
                "walk_beyond_touch_ticks": walk_beyond_ticks,
                "walk_immediate_cost_bps": walk_immediate_cost_bps,
                "configured_beyond_touch_ticks": configured_beyond_ticks,
                "configured_immediate_cost_bps": 10_000.0
                * side
                * (touch - mid)
                / mid
                + configured_beyond_ticks * tick_bps,
                "recorded_fill_cost_bps": 10_000.0
                * side
                * (float(fill.fill_price) - float(fill.mid_at_fill))
                / float(fill.mid_at_fill),
                "top_mid_difference_ticks": (mid - float(fill.mid_at_fill)) / tick_size,
            }
        )
    return pd.DataFrame(rows)


def calibrate_depth_model(
    depth_walks: pd.DataFrame,
    *,
    calibration_end_seconds: float,
    configured_kappa: float,
    over_displayed_multiplier: float,
    tick_size: float,
    participation_edges: Iterable[float],
    minimum_bin_observations: int = 5,
) -> Mapping[str, Any]:
    """Fit kappa on supported early children and audit supported later children."""

    supported = depth_walks.loc[depth_walks["full_l5_supported"]].copy()
    calibration = supported.loc[
        supported["timestamp"].lt(calibration_end_seconds)
    ].copy()
    validation = supported.loc[
        supported["timestamp"].ge(calibration_end_seconds)
    ].copy()
    if calibration.empty or validation.empty:
        raise ValueError("L5 audit requires supported children on both sides of the split")
    x = calibration["adjusted_participation"].to_numpy(dtype="float64")
    y = calibration["walk_beyond_touch_ticks"].to_numpy(dtype="float64")
    weights = calibration["size"].to_numpy(dtype="float64")
    denominator = np.sum(weights * np.square(x))
    candidate_kappa = max(0.0, float(np.sum(weights * x * y) / denominator))

    tick_bps = 10_000.0 * tick_size / validation["mid_l5"]
    validation["candidate_beyond_touch_ticks"] = (
        candidate_kappa * validation["adjusted_participation"]
    )
    validation["candidate_immediate_cost_bps"] = (
        10_000.0
        * validation["side"]
        * (validation["touch_l5"] - validation["mid_l5"])
        / validation["mid_l5"]
        + validation["candidate_beyond_touch_ticks"] * tick_bps
    )
    configured_mae = _weighted_mean(
        (
            validation["configured_immediate_cost_bps"]
            - validation["walk_immediate_cost_bps"]
        ).abs(),
        validation["size"],
    )
    candidate_mae = _weighted_mean(
        (
            validation["candidate_immediate_cost_bps"]
            - validation["walk_immediate_cost_bps"]
        ).abs(),
        validation["size"],
    )
    performance_pass = candidate_mae < configured_mae
    accepted_kappa = candidate_kappa if performance_pass else configured_kappa
    curve = _depth_curve(validation, participation_edges, minimum_bin_observations)
    curve["candidate_accepted"] = performance_pass
    support = _support_curve(depth_walks, participation_edges)
    summary = pd.DataFrame(
        [
            {
                "real_data_label": "REAL_L5_BOOK",
                "model_data_label": "SIMULATED_MODEL",
                "calibration_end_seconds": calibration_end_seconds,
                "calibration_supported_children": len(calibration),
                "validation_supported_children": len(validation),
                "configured_kappa": configured_kappa,
                "l5_candidate_kappa": candidate_kappa,
                "accepted_kappa": accepted_kappa,
                "holdout_performance_pass": performance_pass,
                "simulated_share_fraction_fully_supported_l5": float(
                    depth_walks.loc[depth_walks["full_l5_supported"], "size"].sum()
                    / depth_walks["size"].sum()
                ),
                "mean_l5_walk_cost_bps": _weighted_mean(
                    validation["walk_immediate_cost_bps"], validation["size"]
                ),
                "mean_configured_cost_bps": _weighted_mean(
                    validation["configured_immediate_cost_bps"], validation["size"]
                ),
                "mean_candidate_cost_bps": _weighted_mean(
                    validation["candidate_immediate_cost_bps"], validation["size"]
                ),
                "configured_bias_bps": _weighted_mean(
                    validation["configured_immediate_cost_bps"]
                    - validation["walk_immediate_cost_bps"],
                    validation["size"],
                ),
                "candidate_bias_bps": _weighted_mean(
                    validation["candidate_immediate_cost_bps"]
                    - validation["walk_immediate_cost_bps"],
                    validation["size"],
                ),
                "configured_mae_bps": configured_mae,
                "candidate_mae_bps": candidate_mae,
                "top_mid_alignment_mae_ticks": float(
                    depth_walks["top_mid_difference_ticks"].abs().mean()
                ),
            }
        ]
    )
    return {
        "calibration": calibration,
        "validation": validation,
        "curve": curve,
        "support": support,
        "summary": summary,
        "candidate_kappa": candidate_kappa,
        "accepted_kappa": accepted_kappa,
        "performance_pass": performance_pass,
    }


def audit_parent_costs(
    parent_metrics: pd.DataFrame,
    depth_walks: pd.DataFrame,
    *,
    accepted_kappa: float,
    candidate_kappa: float,
    tick_size: float,
) -> Tuple[pd.DataFrame, pd.DataFrame]:
    """Replace impact only for children fully covered by REAL L5 depth."""

    view = depth_walks.copy()
    view["accepted_impact_ticks"] = (
        accepted_kappa * view["adjusted_participation"]
    )
    view["candidate_impact_ticks"] = (
        candidate_kappa * view["adjusted_participation"]
    )
    view["configured_impact_notional"] = (
        view["size"] * view["configured_beyond_touch_ticks"] * tick_size
    )
    view["candidate_impact_notional"] = (
        view["size"] * view["candidate_impact_ticks"] * tick_size
    )
    view["support_adjusted_impact_notional"] = np.where(
        view["full_l5_supported"],
        view["size"] * view["accepted_impact_ticks"] * tick_size,
        view["configured_impact_notional"],
    )
    view["unsupported_size"] = np.where(view["full_l5_supported"], 0, view["size"])
    impact = view.groupby("parent_id", as_index=False).agg(
        configured_impact_notional=("configured_impact_notional", "sum"),
        candidate_impact_notional=("candidate_impact_notional", "sum"),
        support_adjusted_impact_notional=("support_adjusted_impact_notional", "sum"),
        unsupported_shares=("unsupported_size", "sum"),
    )
    parent = parent_metrics.merge(impact, on="parent_id", how="left", validate="one_to_one")
    denominator = parent["arrival_mid"] * parent["requested_size"]
    for prefix in ("configured", "candidate", "support_adjusted"):
        parent[f"{prefix}_impact_component_bps"] = (
            10_000.0 * parent[f"{prefix}_impact_notional"] / denominator
        )
    parent["supported_model_artifact_bps"] = (
        parent["configured_impact_component_bps"]
        - parent["support_adjusted_impact_component_bps"]
    )
    parent["full_extrapolation_artifact_bps"] = (
        parent["configured_impact_component_bps"]
        - parent["candidate_impact_component_bps"]
    )
    parent["support_adjusted_is_bps"] = (
        parent["implementation_shortfall_bps"] - parent["supported_model_artifact_bps"]
    )
    parent["full_extrapolation_is_sensitivity_bps"] = (
        parent["implementation_shortfall_bps"]
        - parent["full_extrapolation_artifact_bps"]
    )
    parent["unsupported_share_fraction"] = (
        parent["unsupported_shares"] / parent["executed_size"]
    )
    notional = parent["arrival_mid"] * parent["requested_size"]
    portfolio = pd.DataFrame(
        [
            {
                "data_label": "SIMULATED",
                "calibration_source": "REAL_L5_BOOK",
                "parents": len(parent),
                "portfolio_raw_is_bps": _weighted_mean(
                    parent["implementation_shortfall_bps"], notional
                ),
                "unsupported_simulated_share_fraction": float(
                    view.loc[~view["full_l5_supported"], "size"].sum()
                    / view["size"].sum()
                ),
                "portfolio_support_adjusted_is_bps": _weighted_mean(
                    parent["support_adjusted_is_bps"], notional
                ),
                "portfolio_supported_model_artifact_bps": _weighted_mean(
                    parent["supported_model_artifact_bps"], notional
                ),
                "portfolio_rejected_candidate_extrapolation_is_sensitivity_bps": _weighted_mean(
                    parent["full_extrapolation_is_sensitivity_bps"], notional
                ),
                "portfolio_rejected_candidate_extrapolation_artifact_bps": _weighted_mean(
                    parent["full_extrapolation_artifact_bps"], notional
                ),
            }
        ]
    )
    return parent, portfolio


def plot_impact_audit(
    curve: pd.DataFrame,
    support: pd.DataFrame,
) -> Tuple[plt.Figure, np.ndarray]:
    """Overlay REAL L5 walk cost and show supported/unsupported child shares."""

    labels = list(support["participation_bin"])
    reportable = curve.loc[curve["reportable"]].set_index("participation_bin")
    aligned = reportable.reindex(labels)
    figure, axes = plt.subplots(
        2, 1, figsize=(7.5, 6.5), sharex=True, gridspec_kw={"height_ratios": [2, 1]}
    )
    x = np.arange(len(labels))
    axes[0].plot(x, aligned["l5_walk_cost_bps"], marker="o", label="REAL L5 walk VWAP")
    candidate_accepted = bool(curve["candidate_accepted"].iloc[0])
    axes[0].plot(
        x,
        aligned["configured_cost_bps"],
        marker="o",
        label="SIMULATED model (configured κ)",
    )
    axes[0].plot(
        x,
        aligned["candidate_cost_bps"],
        marker="o",
        label=(
            "SIMULATED model (L5 candidate κ)"
            if candidate_accepted
            else "SIMULATED model (L5 candidate κ — REJECTED)"
        ),
    )
    axes[0].set(
        title="Impact realism audit — second execution window holdout",
        ylabel="Immediate one-way cost (bps)",
    )
    axes[0].legend()
    axes[1].bar(
        x,
        100.0 * support["supported_share_fraction"],
        color="#2e7d32",
        label="Covered by L5",
    )
    axes[1].bar(
        x,
        100.0 * support["unsupported_share_fraction"],
        bottom=100.0 * support["supported_share_fraction"],
        color="#b23a48",
        label="Beyond L5",
    )
    axes[1].set_xticks(x)
    axes[1].set_xticklabels(labels, rotation=30, ha="right")
    axes[1].set(
        xlabel="Child size / displayed level-one size",
        ylabel="All child shares (%)",
        title="REAL five-level depth coverage",
    )
    axes[1].legend()
    return figure, axes


def _walk_book(prices: np.ndarray, sizes: np.ndarray, quantity: int) -> Mapping[str, Any]:
    remaining = quantity
    notional = 0.0
    levels_used = 0
    for price, available in zip(prices, sizes):
        take = min(remaining, int(available))
        if take > 0:
            notional += take * float(price)
            remaining -= take
            levels_used += 1
        if remaining == 0:
            break
    supported = remaining == 0
    vwap = notional / quantity if supported else np.nan
    return {
        "supported": supported,
        "vwap": vwap,
        "levels_used": levels_used,
    }


def _depth_curve(
    validation: pd.DataFrame,
    participation_edges: Iterable[float],
    minimum_observations: int,
) -> pd.DataFrame:
    view, labels = _with_bins(validation, participation_edges)
    rows = []
    for label in labels:
        group = view.loc[view["participation_bin"].eq(label)]
        if group.empty:
            continue
        weights = group["size"]
        rows.append(
            {
                "real_data_label": "REAL_L5_BOOK",
                "model_data_label": "SIMULATED_MODEL",
                "participation_bin": label,
                "observations": len(group),
                "shares": int(group["size"].sum()),
                "l5_walk_cost_bps": _weighted_mean(
                    group["walk_immediate_cost_bps"], weights
                ),
                "configured_cost_bps": _weighted_mean(
                    group["configured_immediate_cost_bps"], weights
                ),
                "candidate_cost_bps": _weighted_mean(
                    group["candidate_immediate_cost_bps"], weights
                ),
                "configured_bias_bps": _weighted_mean(
                    group["configured_immediate_cost_bps"]
                    - group["walk_immediate_cost_bps"],
                    weights,
                ),
                "candidate_bias_bps": _weighted_mean(
                    group["candidate_immediate_cost_bps"]
                    - group["walk_immediate_cost_bps"],
                    weights,
                ),
                "reportable": len(group) >= minimum_observations,
            }
        )
    return pd.DataFrame(rows)


def _support_curve(
    depth_walks: pd.DataFrame,
    participation_edges: Iterable[float],
) -> pd.DataFrame:
    view, labels = _with_bins(depth_walks, participation_edges)
    total = view["size"].sum()
    rows = []
    for label in labels:
        group = view.loc[view["participation_bin"].eq(label)]
        supported = int(group.loc[group["full_l5_supported"], "size"].sum())
        unsupported = int(group.loc[~group["full_l5_supported"], "size"].sum())
        rows.append(
            {
                "data_label": "SIMULATED",
                "market_data_label": "REAL_L5_BOOK",
                "participation_bin": label,
                "fills": len(group),
                "supported_shares": supported,
                "unsupported_shares": unsupported,
                "supported_share_fraction": supported / total,
                "unsupported_share_fraction": unsupported / total,
            }
        )
    return pd.DataFrame(rows)


def _with_bins(
    frame: pd.DataFrame,
    participation_edges: Iterable[float],
) -> Tuple[pd.DataFrame, list]:
    edges = sorted(set(float(edge) for edge in participation_edges))
    if not edges or edges[0] != 0.0:
        raise ValueError("Participation edges must start at zero")
    if not np.isinf(edges[-1]):
        edges.append(float("inf"))
    labels = [
        f"{left:g}-{right:g}x" if np.isfinite(right) else f">{left:g}x"
        for left, right in zip(edges, edges[1:])
    ]
    view = frame.copy()
    view["participation_bin"] = pd.cut(
        view["level_one_participation"],
        edges,
        labels=labels,
        include_lowest=True,
        right=True,
    )
    return view, labels


def _book_columns(levels: int) -> list:
    columns = []
    for level in range(1, levels + 1):
        columns.extend(
            [
                f"ask_price_{level}",
                f"ask_size_{level}",
                f"bid_price_{level}",
                f"bid_size_{level}",
            ]
        )
    return columns


def _weighted_mean(values, weights) -> float:
    values = np.asarray(values, dtype="float64")
    weights = np.asarray(weights, dtype="float64")
    valid = np.isfinite(values) & np.isfinite(weights) & (weights > 0)
    if not valid.any():
        return float("nan")
    return float(np.average(values[valid], weights=weights[valid]))