백테스트 불확실성, 선택 편향과 전략 짝비교

# uncertainty.py


```py
"""Backtest performance uncertainty: block bootstrap, HAC SE, PSR/DSR, paired comparisons.

Operates on the daily out-of-sample strategy return series persisted at
``run_log/backtest/{hash}/daily_returns.parquet``. Series-level CIs are computed
via stationary block bootstrap; Sharpe SE uses the López de Prado (2025) closed
form with autocorrelation, skew, and kurtosis corrections; mean-return SE uses
Newey-West HAC. Selection bias across K variants uses the library's
:func:`deflated_sharpe_ratio`. Challenger-vs-baseline uncertainty uses paired
stationary block bootstrap on daily-return differences.

All bootstraps share a single seed so the same call is reproducible. Default
``n_boot=2000``; tune down for sweeps if needed.

Block-length policy
-------------------
``resolve_block_length(case_study, label, returns)`` picks the block length:

1. ``setup.yaml.labels.{label}.rebalance_step`` if present (canonical).
2. Falls back to :func:`ml4t.diagnostic.evaluation.stats._optimal_block_size`.
3. Floored at the label's forward-return horizon (``ret_5d`` → ≥5,
   ``ret_to_expiry`` → keeps the floor at 1, since horizon is variable).

Baseline registry
-----------------
:data:`STAGE_BASELINE` declares the natural baseline for each backtest stage,
used by the Ch20 paired-bootstrap synthesis:

- ``signal``  → equal-weight benchmark (per case study, registered separately)
- ``allocation``    → ``signal`` leader of the same (label, family)
- ``risk_overlay``  → ``allocation`` leader (sized, no overlay)
- ``cost_sensitivity`` → ``risk_overlay`` leader (sized and overlaid, frictionless)

Each stage is benchmarked against the leader of the stage before it, so the
chain follows the order the backtest sequence runs: size positions, apply risk
controls, then measure what realistic costs take off the winner. A stage that a
case study has not run is skipped, and the benchmark falls back to the nearest
earlier stage that has rows.

Per-case-study baselines for the signal stage live in
:data:`SIGNAL_BASELINE_BY_CASE_STUDY`; populate this when the equal-weight
benchmark name in the registry is non-default.
"""

from __future__ import annotations

import hashlib
import json
import re
import warnings
from collections.abc import Iterable, Mapping, Sequence
from dataclasses import dataclass
from itertools import combinations
from typing import Any, Final, Literal, cast

import numpy as np
import polars as pl


class EntireRegistry:
    """The population scope that reads every registered row, retired generations included."""

    __slots__ = ()

    def __repr__(self) -> str:
        return "ENTIRE_REGISTRY"


class NoCarrier:
    """The carrier scope that re-ranks the candidates on raw Sharpe instead of naming one."""

    __slots__ = ()

    def __repr__(self) -> str:
        return "NO_CARRIER"


ENTIRE_REGISTRY: Final = EntireRegistry()
"""Ask the paired and cohort producers to compute over the whole registry.

The scope arguments that decide what a published number covers carry no default, so a
caller that wants the widest reading has to write it down. That is the whole point: the
wide reading is almost never what a strategy-analysis notebook means, and when it was a
default, omitting the argument type-checked, ran, and produced plausible numbers that
were wrong exactly when the registry held a generation the notebook does not report.
Naming it here also makes the wide callers greppable, which they were not.
"""

NO_CARRIER: Final = NoCarrier()
"""Ask ``populate_paired_metrics`` to rank its own candidates on raw Sharpe.

The legacy ranking, kept because the rung-pinned case studies restrict on a dimension
``resolve_canonical_rank1_lineage`` does not know. It is not the canonical selection: it
orders on raw Sharpe and applies neither the common-support re-ranking nor the
restrictions the resolver holds, so a caller that can resolve the lineage should pass it
rather than this.
"""

PredictionScope = Iterable[str] | EntireRegistry
"""The population a published number is computed over: a list of prediction hashes,
or :data:`ENTIRE_REGISTRY`."""

CarrierScope = Mapping[str, Any] | NoCarrier
"""The lineage pairs #2-6 are pinned to: a ``resolve_canonical_rank1_lineage`` result,
or :data:`NO_CARRIER`."""


def periods_per_year_from_setup(case_study: str) -> int:
    """Resolve periods-per-year from the case study's setup.yaml.

    Reads ``evaluation.periods_per_year`` from
    ``case_studies/{case_study}/config/setup.yaml``. This is the
    annualization convention of the **daily_returns** grid that the
    backtester actually writes (NYSE-like 5d/wk → 252, 7d/wk crypto →
    365, genuinely monthly us_firm → 12). It is NOT the trade cadence
    or rebalance frequency.

    Raises FileNotFoundError if the case study has no setup.yaml, or
    KeyError if `evaluation.periods_per_year` is not declared — both
    are programming errors that callers should fix upstream rather
    than silently absorb.
    """
    import yaml

    from utils.paths import get_case_study_dir

    setup_path = get_case_study_dir(case_study) / "config" / "setup.yaml"
    with setup_path.open() as f:
        setup = yaml.safe_load(f)
    evaluation = setup.get("evaluation", {}) if isinstance(setup, dict) else {}
    if "periods_per_year" not in evaluation:
        raise KeyError(
            f"{setup_path} is missing evaluation.periods_per_year; add it "
            "(252 for daily 5d/wk markets, 365 for 7d/wk crypto, 12 for "
            "genuinely monthly us_firm)."
        )
    return int(evaluation["periods_per_year"])


# ---------------------------------------------------------------------------
# Block-length resolver
# ---------------------------------------------------------------------------


_LABEL_HORIZON_RE = re.compile(r"(\d+)\s*d\b")


def _label_horizon_floor(label: str | None) -> int:
    """Best-effort horizon (days) implied by a label name; 1 if unknown."""
    if not label:
        return 1
    match = _LABEL_HORIZON_RE.search(label)
    if match:
        return max(int(match.group(1)), 1)
    return 1


def resolve_block_length(
    case_study: str | None,
    label: str | None,
    returns: np.ndarray,
    *,
    explicit: int | None = None,
) -> int:
    """Resolve block length: rebalance_step → optimal → floor at label horizon."""
    if explicit is not None and explicit > 0:
        return int(explicit)

    rebalance_step: int | None = None
    if case_study and label:
        try:
            from case_studies.utils.backtest_loaders import get_rebalance_step

            rebalance_step = int(get_rebalance_step(case_study, label))
        except Exception:
            rebalance_step = None

    floor = _label_horizon_floor(label)

    if rebalance_step and rebalance_step > 0:
        return max(rebalance_step, floor)

    scale = float(np.std(returns))
    scale_floor = np.finfo(float).eps * max(1.0, abs(float(np.mean(returns))))
    if returns.size >= 10 and scale <= scale_floor:
        return max(int(returns.size ** (1 / 3)), floor, 1)

    from ml4t.diagnostic.evaluation.stats import _optimal_block_size

    optimal = int(round(float(_optimal_block_size(returns))))
    return max(optimal, floor, 1)


# ---------------------------------------------------------------------------
# Baseline registry (Ch20 paired-bootstrap synthesis)
# ---------------------------------------------------------------------------


#: Stage order of the backtest sequence. Each stage's benchmark is the leader of
#: the nearest preceding stage that has rows.
STAGE_SEQUENCE: tuple[str, ...] = (
    "signal",
    "allocation",
    "risk_overlay",
    "cost_sensitivity",
)


#: The block of a strategy spec each stage introduces. ``cost_sensitivity`` has no entry
#: because it is terminal - nothing is ever built on top of a cost sweep.
STAGE_CARRIER_BLOCK: dict[str, str] = {
    "signal": "signal",
    "allocation": "allocation",
    "risk_overlay": "risk",
}


def carried_blocks(stage: str) -> tuple[str, ...]:
    """Every strategy block a backtest at ``stage`` has inherited or introduced."""
    if stage not in STAGE_SEQUENCE:
        return ()
    upto = STAGE_SEQUENCE[: STAGE_SEQUENCE.index(stage) + 1]
    return tuple(STAGE_CARRIER_BLOCK[s] for s in upto if s in STAGE_CARRIER_BLOCK)


def descends_from(challenger: dict, baseline: dict, baseline_stage: str) -> bool:
    """Is ``challenger`` a strategy built on top of ``baseline``?

    `champion_lineage` takes the best backtest at each stage independently, so its
    entries can be siblings rather than parent and child - two strategies that branch
    off the same allocation carrier, say, one adding a risk overlay and one sweeping
    costs. Comparing those two attributes the whole difference between two unrelated
    strategies to whichever stage happens to come second in the chain.

    Descent requires the challenger to match the baseline on the *whole prefix* the
    baseline carries, not only on the block its own stage introduced. A shared
    prediction hash fixes the predictions and nothing else: signal-stage runs vary
    the signal method and ``top_k``, so an allocation leader can differ from the
    signal leader in the one place the comparison is meant to hold fixed. Checking a
    single block would pass it.
    """
    return all(challenger.get(b) == baseline.get(b) for b in carried_blocks(baseline_stage))


STAGE_BASELINE: dict[str, str] = {
    "signal": "equal_weight",
    "allocation": "signal_leader",
    "risk_overlay": "allocation_leader",
    "cost_sensitivity": "risk_overlay_leader",
}


SIGNAL_BASELINE_BY_CASE_STUDY: dict[str, str] = {
    "etfs": "equal_weight",
    "nasdaq100_microstructure": "equal_weight",
    "sp500_equity_option_analytics": "equal_weight",
    "sp500_options": "equal_weight",
    "us_firm_characteristics": "equal_weight",
    "us_equities_panel": "equal_weight",
    "fx_pairs": "equal_weight",
    "crypto_perps_funding": "equal_weight",
    "cme_futures": "equal_weight",
}


# ---------------------------------------------------------------------------
# Series-level uncertainty
# ---------------------------------------------------------------------------


@dataclass
class _Stats:
    sharpe: float
    sortino: float
    ann_return: float
    volatility: float
    max_drawdown: float
    calmar: float


def _sample_stats(returns: np.ndarray, periods_per_year: int) -> _Stats:
    """Point-estimate statistics on a return series."""
    if len(returns) < 2:
        return _Stats(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
    mu = float(np.mean(returns))
    sd = float(np.std(returns, ddof=1))
    sharpe = (mu / sd * np.sqrt(periods_per_year)) if sd > 0 else 0.0
    # Downside deviation averages the squared shortfall over EVERY period, not over the
    # periods that fell. Dividing by the count of negative returns instead inflates the
    # ratio by sqrt(n / n_negative), and since `backtest_metrics.sortino` is written by
    # the engine's standard definition, that made the point estimate and the interval
    # around it two different estimators: on us_firm_characteristics' validation rank-1,
    # 99 periods with 20 negative, a stored 13.876 against a bootstrap CI of
    # [4.22, 9.65] - the point outside its own interval, and a forest plot that could
    # not be drawn.
    shortfall = np.minimum(returns, 0.0)
    dsd = float(np.sqrt(np.mean(shortfall**2)))
    sortino = (mu / dsd * np.sqrt(periods_per_year)) if dsd > 0 else 0.0
    cum = np.cumprod(1.0 + returns)
    total_return = float(cum[-1] - 1.0)
    n_years = len(returns) / periods_per_year
    # Guard the negative-cumulative-return case: if the strategy lost more
    # than 100 % (cumulative growth ≤ 0), the geometric mean would be
    # complex. Report `ann_return = -1` (total wipeout) instead.
    base = 1.0 + total_return
    if n_years <= 0:
        ann_return = 0.0
    elif base <= 0.0:
        ann_return = -1.0
    else:
        ann_return = float(base ** (1.0 / n_years) - 1.0)
    vol = sd * np.sqrt(periods_per_year)
    running_max = np.maximum.accumulate(cum)
    # Avoid div-by-zero / negative running_max (post-bankruptcy paths).
    safe_max = np.where(running_max > 0, running_max, np.nan)
    dd = (cum - running_max) / safe_max
    max_dd_raw = float(np.nanmin(dd)) if np.any(np.isfinite(dd)) else 0.0
    max_dd = max_dd_raw if np.isfinite(max_dd_raw) else 0.0
    calmar = (ann_return / abs(max_dd)) if max_dd < 0 else 0.0
    return _Stats(sharpe, sortino, ann_return, vol, max_dd, calmar)


def _newey_west_mean_se(returns: np.ndarray, lag: int) -> float:
    """Newey-West HAC standard error of the sample mean."""
    n = len(returns)
    if n < 3:
        return float("nan")
    r = returns - np.mean(returns)
    gamma0 = float(np.dot(r, r) / n)
    s = gamma0
    for h in range(1, min(lag, n - 1) + 1):
        gamma_h = float(np.dot(r[:-h], r[h:]) / n)
        w = 1.0 - h / (lag + 1.0)
        s += 2.0 * w * gamma_h
    s = max(s, 0.0)
    return float(np.sqrt(s / n))


def _sharpe_se_lo(returns: np.ndarray, periods_per_year: int) -> float:
    """LdP-2025 Sharpe SE with autocorrelation, skewness, kurtosis."""
    from ml4t.diagnostic.evaluation.stats import compute_sharpe_variance

    n = len(returns)
    if n < 4:
        return float("nan")
    mu = float(np.mean(returns))
    sd = float(np.std(returns, ddof=1))
    scale_floor = np.finfo(float).eps * max(1.0, abs(mu))
    if sd <= scale_floor:
        return float("nan")
    sr = mu / sd  # native frequency
    centered = returns - mu
    m2 = float(np.mean(centered**2))
    if m2 <= scale_floor**2:
        return float("nan")
    skew = float(np.mean(centered**3) / m2**1.5)
    kurt = float(np.mean(centered**4) / m2**2)  # Pearson convention (normal=3)
    previous = returns[:-1]
    following = returns[1:]
    if float(np.std(previous)) == 0.0 or float(np.std(following)) == 0.0:
        rho = 0.0
    else:
        rho = float(np.corrcoef(previous, following)[0, 1])
    if not np.isfinite(rho) or abs(rho) >= 0.999:
        rho = 0.0
    var = compute_sharpe_variance(
        sharpe=sr,
        n_samples=n,
        skewness=skew,
        kurtosis=kurt,
        autocorrelation=rho,
        n_trials=1,
    )
    if var <= 0 or not np.isfinite(var):
        return float("nan")
    se_native = float(np.sqrt(var))
    return se_native * np.sqrt(periods_per_year)


def _stationary_bootstrap_metrics(
    returns: np.ndarray,
    *,
    periods_per_year: int,
    block_length: int,
    n_boot: int,
    seed: int,
) -> dict[str, np.ndarray]:
    """Run a stationary bootstrap and return arrays of resampled metrics."""
    from ml4t.diagnostic.evaluation.stats import _stationary_bootstrap_indices

    rng = np.random.default_rng(seed)
    sharpes = np.empty(n_boot)
    sortinos = np.empty(n_boot)
    ann_rets = np.empty(n_boot)
    vols = np.empty(n_boot)
    max_dds = np.empty(n_boot)
    calmars = np.empty(n_boot)

    for i in range(n_boot):
        idx = _stationary_bootstrap_indices(len(returns), float(block_length), rng)
        sample = returns[idx]
        stats = _sample_stats(sample, periods_per_year)
        sharpes[i] = stats.sharpe
        sortinos[i] = stats.sortino
        ann_rets[i] = stats.ann_return
        vols[i] = stats.volatility
        max_dds[i] = stats.max_drawdown
        calmars[i] = stats.calmar

    return {
        "sharpe": sharpes,
        "sortino": sortinos,
        "ann_return": ann_rets,
        "volatility": vols,
        "max_drawdown": max_dds,
        "calmar": calmars,
    }


def _percentile_ci(arr: np.ndarray, alpha: float = 0.05) -> tuple[float, float]:
    arr = arr[np.isfinite(arr)]
    if arr.size == 0:
        return float("nan"), float("nan")
    return (
        float(np.percentile(arr, 100 * alpha / 2)),
        float(np.percentile(arr, 100 * (1 - alpha / 2))),
    )


def compute_backtest_uncertainty(
    daily_returns: np.ndarray | pl.Series | pl.DataFrame,
    *,
    periods_per_year: int = 252,
    block_length: int | None = None,
    case_study: str | None = None,
    label: str | None = None,
    n_boot: int = 2000,
    seed: int = 0,
) -> dict[str, float]:
    """Compute series-level uncertainty for one backtest.

    Returns a flat dict suitable for upsert into ``backtest_metrics``:

    - ``sharpe_se_lo``                — Lo / LdP-2025 closed-form SE
    - ``sharpe_ci95_lo`` / ``_hi``    — block-bootstrap percentile CI
    - ``sortino_ci95_lo`` / ``_hi``
    - ``ann_return_hac_se``           — Newey-West HAC SE of mean return (annualized)
    - ``ann_return_ci95_lo`` / ``_hi``
    - ``max_dd_ci95_lo`` / ``_hi``
    - ``calmar_ci95_lo`` / ``_hi``
    - ``psr_pvalue``                  — 1 − P(true SR > 0) under PSR
    - ``bootstrap_block_length``
    - ``bootstrap_n``
    """
    arr = _coerce_returns(daily_returns)
    if arr.size < 4:
        return {}

    block = resolve_block_length(case_study, label, arr, explicit=block_length)
    boot = _stationary_bootstrap_metrics(
        arr,
        periods_per_year=periods_per_year,
        block_length=block,
        n_boot=n_boot,
        seed=seed,
    )

    point = _sample_stats(arr, periods_per_year)
    sharpe_se = _sharpe_se_lo(arr, periods_per_year)

    # NW lag: at least the bootstrap block (≈ rebalance step)
    nw_lag = max(block - 1, int(np.floor(4 * (len(arr) / 100) ** (2 / 9))))
    mean_se_native = _newey_west_mean_se(arr, lag=nw_lag)
    ann_return_hac_se = (
        mean_se_native * periods_per_year if np.isfinite(mean_se_native) else float("nan")
    )

    # PSR
    psr_pvalue = float("nan")
    try:
        from ml4t.diagnostic.evaluation.stats import deflated_sharpe_ratio

        psr = deflated_sharpe_ratio(arr, periods_per_year=periods_per_year)
        psr_pvalue = float(psr.p_value)
    except Exception:
        psr_pvalue = float("nan")

    sh_lo, sh_hi = _percentile_ci(boot["sharpe"])
    so_lo, so_hi = _percentile_ci(boot["sortino"])
    ar_lo, ar_hi = _percentile_ci(boot["ann_return"])
    md_lo, md_hi = _percentile_ci(boot["max_drawdown"])
    cl_lo, cl_hi = _percentile_ci(boot["calmar"])

    return {
        "sharpe_se_lo": _to_float(sharpe_se),
        "sharpe_ci95_lo": _to_float(sh_lo),
        "sharpe_ci95_hi": _to_float(sh_hi),
        "sortino_ci95_lo": _to_float(so_lo),
        "sortino_ci95_hi": _to_float(so_hi),
        "ann_return_hac_se": _to_float(ann_return_hac_se),
        "ann_return_ci95_lo": _to_float(ar_lo),
        "ann_return_ci95_hi": _to_float(ar_hi),
        "max_dd_ci95_lo": _to_float(md_lo),
        "max_dd_ci95_hi": _to_float(md_hi),
        "calmar_ci95_lo": _to_float(cl_lo),
        "calmar_ci95_hi": _to_float(cl_hi),
        "psr_pvalue": _to_float(psr_pvalue),
        "bootstrap_block_length": float(block),
        "bootstrap_n": float(n_boot),
    }


# ---------------------------------------------------------------------------
# Paired uncertainty: challenger vs baseline
# ---------------------------------------------------------------------------


def joint_returns(
    challenger: np.ndarray | pl.Series,
    baseline: np.ndarray | pl.Series,
    *,
    challenger_overlays_baseline: bool = False,
) -> tuple[np.ndarray, np.ndarray]:
    """Coerce a paired return series to the precondition of a paired bootstrap.

    :func:`compute_paired_uncertainty` requires two arrays of the same length whose position
    ``i`` is the same session on both sides, and it refuses the pair rather than bootstrap a
    misaligned one. Coercing each side on its own does not deliver that: ``_coerce_returns``
    drops non-finite values and the leading run of zeros per series, so two series with
    different amounts of leading inactivity part company. Joining on the timestamp beforehand
    does not save it either, because the per-side trim happens after.

    So both decisions are taken once, over both series: keep a session only where both sides
    are finite, then start where the comparison becomes defined.

    Where it becomes defined depends on what the pair is, which is why the caller has to say.
    A leading flat run on the challenger has two possible meanings and they are
    indistinguishable in the numbers:

    ``challenger_overlays_baseline=False`` (the default, and the strategy-versus-benchmark
        case): the two series are independent, each live from its own first traded session.
        A strategy has a warmup prefix before its first signal while an equal-weight
        benchmark is invested from the first joined session, and those rows are pre-sample
        for the strategy rather than a result. The sample starts where **both** are trading,
        which the code below reads as the first session on which both returns are non-zero;
        see the comment there for the difference and why it is preserved.

    ``challenger_overlays_baseline=True`` (the risk-overlay case): the challenger runs on top
        of the baseline, so both are live from the same session and a flat challenger there
        is a position it chose to hold - the largest instance of the effect the comparison
        exists to measure. Starting where both traded would delete exactly those rows and pull
        the measured difference toward zero in the direction the overlay is under test. The
        sample starts where **either** has traded.

    Returns two empty arrays when no session qualifies.
    """
    c = _as_return_array(challenger)
    b = _as_return_array(baseline)
    if c.size != b.size:
        raise ValueError(
            f"a paired series must arrive aligned; got {c.size} and {b.size} observations"
        )
    finite = np.isfinite(c) & np.isfinite(b)
    c, b = c[finite], b[finite]
    if c.size == 0:
        return c, b
    if challenger_overlays_baseline:
        # Either side having traded starts the sample, so the first index where anything is
        # non-zero: the earlier of the two firsts, or nothing if neither ever traded.
        first_c = np.flatnonzero(c != 0.0)
        first_b = np.flatnonzero(b != 0.0)
        starts = [int(x[0]) for x in (first_c, first_b) if x.size]
        if not starts:
            return c[:0], b[:0]
        start = min(starts)
    else:
        # The first session on which both are simultaneously non-zero, which is the rule the
        # per-case-study producer and `20_strategy_synthesis/01_aggregate_synthesis.py` have
        # both applied since they were split apart. It is not quite the rule the paragraph
        # above states: the later starter's own first session is skipped when the other side
        # happens to post an exactly zero return on it, and those observations are live on
        # both series. Correcting that moves every default pair in the registry and obliges a
        # re-execution of the Chapter 20 synthesis, so it is left as it stands here rather
        # than changed underneath a comparison this function was only asked to make paired.
        both = np.flatnonzero((c != 0.0) & (b != 0.0))
        if not both.size:
            return c[:0], b[:0]
        start = int(both[0])
    return c[start:], b[start:]


def compute_paired_uncertainty(
    challenger: np.ndarray | pl.Series,
    baseline: np.ndarray | pl.Series,
    *,
    periods_per_year: int = 252,
    block_length: int | None = None,
    case_study: str | None = None,
    label: str | None = None,
    n_boot: int = 2000,
    seed: int = 0,
    challenger_overlays_baseline: bool = False,
) -> dict[str, float]:
    """Paired stationary bootstrap on daily-return differences.

    The two series must arrive the same length and aligned by date, so that position ``i``
    is the same session on both sides; a pair that does not is refused with an empty mapping
    rather than truncated. Which rows to drop is then decided over both series at once by
    :func:`joint_returns`, so a caller does not have to coerce them beforehand.
    ``challenger_overlays_baseline`` is passed straight through and says which pair this is;
    read that function before choosing it, because the default is right for a strategy
    against a benchmark and wrong for a risk overlay against its carrier.

    Returns a flat dict for upsert into ``backtest_paired_metrics``, and an empty mapping
    when fewer than four sessions survive.
    """
    from ml4t.diagnostic.evaluation.stats import _stationary_bootstrap_indices

    c_raw = _as_return_array(challenger)
    b_raw = _as_return_array(baseline)
    # Caller's contract: the two series arrive pre-aligned by timestamp, so position i is
    # the same session on both sides. Nothing here can recover that if they do not, because
    # truncating to the shorter one would compare different sessions. Refuse instead.
    if c_raw.size != b_raw.size:
        return {}
    # The coercion is taken once over both series rather than per side. `_coerce_returns`
    # trims each series' own leading run of zeros, which is exactly what a risk overlay
    # produces - it sits out sessions its carrier trades - and the two arrays then part
    # company, so the size check above refused every overlay in `17_risk_management`.
    c, b = joint_returns(c_raw, b_raw, challenger_overlays_baseline=challenger_overlays_baseline)
    if c.size < 4:
        return {}

    diff = c - b
    block = resolve_block_length(case_study, label, diff, explicit=block_length)

    point_c = _sample_stats(c, periods_per_year)
    point_b = _sample_stats(b, periods_per_year)
    sharpe_diff = point_c.sharpe - point_b.sharpe
    ret_diff = point_c.ann_return - point_b.ann_return
    max_dd_diff = point_c.max_drawdown - point_b.max_drawdown

    # Information ratio on the diff series; require sd above a real-world floor
    # (1bp/day) — degenerate near-constant differences are not informative
    sd_diff = float(np.std(diff, ddof=1))
    info_ratio = (
        float(np.mean(diff) / sd_diff * np.sqrt(periods_per_year))
        if sd_diff > 1e-6
        else float("nan")
    )

    # Paired bootstrap: same indices applied to both series
    rng = np.random.default_rng(seed)
    sharpe_diffs = np.empty(n_boot)
    ret_diffs = np.empty(n_boot)
    max_dd_diffs = np.empty(n_boot)
    irs = np.empty(n_boot)
    wins = 0

    for i in range(n_boot):
        idx = _stationary_bootstrap_indices(c.size, float(block), rng)
        cs = _sample_stats(c[idx], periods_per_year)
        bs = _sample_stats(b[idx], periods_per_year)
        sharpe_diffs[i] = cs.sharpe - bs.sharpe
        ret_diffs[i] = cs.ann_return - bs.ann_return
        max_dd_diffs[i] = cs.max_drawdown - bs.max_drawdown
        d = c[idx] - b[idx]
        sd = float(np.std(d, ddof=1))
        irs[i] = float(np.mean(d) / sd * np.sqrt(periods_per_year)) if sd > 1e-6 else float("nan")
        if cs.sharpe > bs.sharpe:
            wins += 1

    sd_lo, sd_hi = _percentile_ci(sharpe_diffs)
    rd_lo, rd_hi = _percentile_ci(ret_diffs)
    mdd_lo, mdd_hi = _percentile_ci(max_dd_diffs)
    ir_lo, ir_hi = _percentile_ci(irs)

    # Two-sided bootstrap p-value for sharpe_diff != 0 (centered around the bootstrap mean)
    centered = sharpe_diffs - np.mean(sharpe_diffs)
    p_value = float(np.mean(np.abs(centered) >= abs(sharpe_diff)))

    return {
        "sharpe_diff": _to_float(sharpe_diff),
        "sharpe_diff_ci95_lo": _to_float(sd_lo),
        "sharpe_diff_ci95_hi": _to_float(sd_hi),
        "ret_diff": _to_float(ret_diff),
        "ret_diff_ci95_lo": _to_float(rd_lo),
        "ret_diff_ci95_hi": _to_float(rd_hi),
        "max_dd_diff": _to_float(max_dd_diff),
        "max_dd_diff_ci95_lo": _to_float(mdd_lo),
        "max_dd_diff_ci95_hi": _to_float(mdd_hi),
        "info_ratio": _to_float(info_ratio),
        "info_ratio_ci95_lo": _to_float(ir_lo),
        "info_ratio_ci95_hi": _to_float(ir_hi),
        "prob_challenger_wins": float(wins) / float(n_boot),
        "p_value": p_value,
        "bootstrap_block_length": float(block),
        "bootstrap_n": float(n_boot),
    }


def compute_independent_diff_uncertainty(
    challenger: np.ndarray | pl.Series,
    baseline: np.ndarray | pl.Series,
    *,
    periods_per_year: int = 252,
    block_length: int | None = None,
    case_study: str | None = None,
    label: str | None = None,
    n_boot: int = 2000,
    seed: int = 0,
) -> dict[str, float]:
    """Difference CI for two return series that share no timestamps.

    Use when challenger and baseline come from non-overlapping windows (e.g.
    holdout vs validation of the same lineage). Each side is bootstrapped over its
    own window and the difference distribution is formed from those draws.

    What is independent here is the two *resampling* draws, and that is forced by
    the windows sharing no observations: there is no difference series to resample,
    so there is nothing to pair on. It is not a claim that the two Sharpes are
    independent. They are not — the same strategy, the same market factor and a
    volatility regime spanning the boundary make them dependent whether or not
    their windows touch — and nothing below needs them to be.

    What the interval covers, and what it does not
    ----------------------------------------------
    Resampling inside a window conditions on that window's realized returns, so
    each side's spread is the sampling noise *given* those returns. The interval is
    calibrated for the difference between the two windows' own population Sharpes,
    and the dependence between the windows does not disturb that: in a Gaussian
    simulation (n=126 a side, 400 replications, 400 draws, a mean shift common to
    both windows) coverage of the nominal 95% interval was 0.945 to 0.948 at every
    correlation between the two windows' shifts from +1 through 0 to -1.

    It is not calibrated for the question the validation-to-holdout decay is
    usually read as asking — the strategy has one edge, is this gap noise? A regime
    that lands differently on the two windows moves their Sharpes apart, and a
    resampler that never looks outside either window cannot see it. Coverage of
    that target in the same simulation: 0.945 when the two windows carry an
    identical shift (nothing to miss), 0.873 when the shifts are independent, 0.850
    when they are opposed. Mean interval width was 7.82 in all three, because the
    interval cannot widen for something it cannot see.

    So the error runs toward under-coverage and never toward over-coverage, and the
    reading that this interval is conservative because it ignores a positive
    covariance is wrong: the same conditioning that drops the covariance term drops
    it from both marginals too, and what is left uncovered is the part of the
    regime the two windows do not share. There is no fix inside a resampling
    scheme — estimating that term needs a model of how a regime carries across the
    boundary, not a different resample — so this is a limit on the reading rather
    than a defect in the arithmetic. A decay this interval leaves unresolved is
    weaker evidence of stability than the same interval over one window would be.

    Returns the same dict shape as :func:`compute_paired_uncertainty` so registry
    callers are interchangeable. ``info_ratio`` columns are NaN — there is no diff
    *series* to ratio when the windows are disjoint. Block length is resolved per
    side from ``(case_study, label)``.
    """
    from ml4t.diagnostic.evaluation.stats import _stationary_bootstrap_indices

    c = _coerce_returns(challenger)
    b = _coerce_returns(baseline)
    if c.size < 4 or b.size < 4:
        return {}

    # Resolve block length per side: disjoint windows can have different
    # autocorrelation structure (different volatility regimes / sample
    # sizes), so a block tuned to one side would under- or over-state
    # bootstrap variance on the other.
    block_c = resolve_block_length(case_study, label, c, explicit=block_length)
    block_b = resolve_block_length(case_study, label, b, explicit=block_length)

    point_c = _sample_stats(c, periods_per_year)
    point_b = _sample_stats(b, periods_per_year)
    sharpe_diff = point_c.sharpe - point_b.sharpe
    ret_diff = point_c.ann_return - point_b.ann_return
    max_dd_diff = point_c.max_drawdown - point_b.max_drawdown

    rng = np.random.default_rng(seed)
    sharpes_c = np.empty(n_boot)
    sharpes_b = np.empty(n_boot)
    rets_c = np.empty(n_boot)
    rets_b = np.empty(n_boot)
    mdds_c = np.empty(n_boot)
    mdds_b = np.empty(n_boot)

    for i in range(n_boot):
        idx_c = _stationary_bootstrap_indices(c.size, float(block_c), rng)
        idx_b = _stationary_bootstrap_indices(b.size, float(block_b), rng)
        cs = _sample_stats(c[idx_c], periods_per_year)
        bs = _sample_stats(b[idx_b], periods_per_year)
        sharpes_c[i] = cs.sharpe
        sharpes_b[i] = bs.sharpe
        rets_c[i] = cs.ann_return
        rets_b[i] = bs.ann_return
        mdds_c[i] = cs.max_drawdown
        mdds_b[i] = bs.max_drawdown

    sharpe_diffs = sharpes_c - sharpes_b
    ret_diffs = rets_c - rets_b
    max_dd_diffs = mdds_c - mdds_b
    wins = float(np.sum(sharpes_c > sharpes_b))

    sd_lo, sd_hi = _percentile_ci(sharpe_diffs)
    rd_lo, rd_hi = _percentile_ci(ret_diffs)
    mdd_lo, mdd_hi = _percentile_ci(max_dd_diffs)

    centered = sharpe_diffs - np.mean(sharpe_diffs)
    p_value = float(np.mean(np.abs(centered) >= abs(sharpe_diff)))

    return {
        "sharpe_diff": _to_float(sharpe_diff),
        "sharpe_diff_ci95_lo": _to_float(sd_lo),
        "sharpe_diff_ci95_hi": _to_float(sd_hi),
        "ret_diff": _to_float(ret_diff),
        "ret_diff_ci95_lo": _to_float(rd_lo),
        "ret_diff_ci95_hi": _to_float(rd_hi),
        "max_dd_diff": _to_float(max_dd_diff),
        "max_dd_diff_ci95_lo": _to_float(mdd_lo),
        "max_dd_diff_ci95_hi": _to_float(mdd_hi),
        "info_ratio": float("nan"),
        "info_ratio_ci95_lo": float("nan"),
        "info_ratio_ci95_hi": float("nan"),
        "prob_challenger_wins": wins / float(n_boot),
        "p_value": p_value,
        # Schema column is single-valued; report the larger block so the
        # value is conservative w.r.t. autocorrelation. n_c / n_b expose
        # the actual per-side post-coerce sample sizes for callers that
        # need an accurate "n_overlap"-equivalent on the disjoint path.
        "bootstrap_block_length": float(max(block_c, block_b)),
        "bootstrap_block_length_c": float(block_c),
        "bootstrap_block_length_b": float(block_b),
        "bootstrap_n": float(n_boot),
        "n_c": float(c.size),
        "n_b": float(b.size),
    }


# ---------------------------------------------------------------------------
# Selection adjustment across K variants (DSR + reality check + PBO + MinTRL)
# ---------------------------------------------------------------------------


def compute_selection_adjustment(
    returns_by_variant: dict[str, np.ndarray | pl.Series],
    *,
    periods_per_year: int = 252,
) -> dict[str, Any]:
    """Selection-bias adjustment across K candidate strategies — **raw-K only**.

    .. deprecated::
        Returns ``dsr / dsr_pvalue / expected_max_sharpe / min_trl_periods``
        using raw trial counts (no Marchenko-Pastur or effective-rank
        correlation correction), which overcounts trials when variants are
        correlated. The recommended replacement is
        :func:`compute_cohort_metrics`, which surfaces raw / MP / ER DSR
        alongside RAS and is persisted to the ``cohort_metrics`` table.
        Consumers should read ``cohort_metrics`` (e.g. via
        :func:`case_studies.utils.notebook_render.selection_adjusted_leader_table`
        or :meth:`BacktestExplorer.deflated_sharpe`) rather than
        recomputing from this helper.

    Combines:

    - **DSR** for the best-of-K leader (Sharpe haircut for selection bias)
    - **Expected max Sharpe** under the null
    - **MinTRL** — periods needed for leader to reach significance at α=0.05
    - **White's reality check** — bootstrap p-value against the best benchmark
    - **PBO** — Probability of Backtest Overfitting across folds (caller must
      pass per-fold returns by variant via ``returns_by_variant`` keyed
      ``"{variant}__fold{i}"``); skipped if no fold-level keys are present.
    """
    from ml4t.diagnostic.evaluation.stats import (
        compute_min_trl,
        deflated_sharpe_ratio,
    )

    # Series prep
    arrays = {}
    for name, ret in returns_by_variant.items():
        a = _coerce_returns(ret)
        if a.size >= 4 and float(np.std(a, ddof=1)) > 1e-10:
            arrays[name] = a
    if not arrays:
        return {}

    names = list(arrays.keys())
    arr_list = [arrays[n] for n in names]
    sharpes = {n: _sample_stats(arrays[n], periods_per_year).sharpe for n in names}
    leader = max(sharpes, key=lambda name: sharpes[name])

    out: dict[str, Any] = {
        "leader": leader,
        "leader_sharpe": float(sharpes[leader]),
        "k_variants": float(len(arr_list)),
    }

    # DSR
    try:
        dsr = deflated_sharpe_ratio(arr_list, periods_per_year=periods_per_year)
        out["dsr"] = float(dsr.deflated_sharpe)
        out["dsr_pvalue"] = float(dsr.p_value)
        out["expected_max_sharpe"] = float(dsr.expected_max_sharpe)
        out["min_trl_periods"] = float(dsr.min_trl)
        out["dsr_significant"] = bool(dsr.is_significant)
    except Exception:
        pass

    # MinTRL standalone (for the leader, against SR=0 benchmark)
    try:
        leader_arr = arrays[leader]
        mtrl = compute_min_trl(
            leader_arr,
            periods_per_year=periods_per_year,
        )
        out["leader_min_trl"] = float(mtrl.min_trl)
    except Exception:
        pass

    return out


def compute_reality_check(
    challenger_returns: dict[str, np.ndarray | pl.Series],
    benchmark_returns: np.ndarray | pl.Series | pl.DataFrame,
    *,
    block_size: int | None = None,
    n_bootstrap: int = 2000,
    seed: int = 0,
) -> dict[str, Any]:
    """White's reality check: do any of K challengers beat the benchmark?

    Returns ``{p_value, test_statistic, best_strategy, k_strategies}``.
    """
    from ml4t.diagnostic.evaluation.stats import whites_reality_check

    bench = _coerce_returns(benchmark_returns)
    names = list(challenger_returns.keys())
    arrs: list[np.ndarray] = []
    keep_names: list[str] = []
    for n in names:
        a = _coerce_returns(challenger_returns[n])
        if a.size == bench.size and float(np.std(a, ddof=1)) > 1e-10:
            arrs.append(a)
            keep_names.append(n)
    if not arrs:
        return {}
    strategies = np.column_stack(arrs)
    rc = whites_reality_check(
        returns_benchmark=bench,
        returns_strategies=strategies,
        bootstrap_samples=n_bootstrap,
        block_size=block_size,
        random_state=seed,
    )
    best_idx = _leader_index(np.mean(strategies - bench.reshape(-1, 1), axis=0), keep_names)
    return {
        "reality_check_pvalue": float(rc.get("p_value", float("nan"))),
        "reality_check_statistic": float(rc.get("test_statistic", float("nan"))),
        "reality_check_best": keep_names[best_idx],
        "k_strategies": float(len(keep_names)),
    }


# ---------------------------------------------------------------------------
# Internal helpers
# ---------------------------------------------------------------------------


def _leader_index(scores: np.ndarray, names: Sequence[str]) -> int:
    """The index of the highest score, and of the lowest hash among those that tie it.

    ``argmax`` returns the FIRST maximum, and "first" here is the order the caller's
    mapping happened to be built in - which is SQLite's row order for a cohort listing,
    a property of one database file rather than of the results. Two readers with the
    same registry content and different insert histories get different leaders, and so
    do the same rows before and after a rebuild.

    Exact ties are not exotic in this data; they are produced by construction. Two
    backtests of one model that differ only in a specification field the returns do not
    depend on - an overlay that never triggers, a re-run that reproduces its inputs -
    book identical return series and therefore identical Sharpes to the last bit.
    Measured 2026-09-07 across the nine production registries: 17 of 253 cohorts have a
    tied leader, and in 12 of them the tie decides which hash is reported.

    The hash is the tie-break because it is the only key that is a property of the result
    rather than of how it was stored, so the same registry content answers the same way
    everywhere. NaN never compares equal, so a NaN score cannot enter the tied set.
    """
    best = np.nanmax(scores)
    tied = [i for i in range(len(names)) if scores[i] == best]
    return min(tied, key=lambda i: names[i])


def _as_return_array(x: np.ndarray | pl.Series | pl.DataFrame) -> np.ndarray:
    """The return series as a float array, with no row dropped.

    Separate from :func:`_coerce_returns` because a paired comparison has to decide which
    rows to drop over both series at once; see :func:`joint_returns`.
    """
    if isinstance(x, pl.DataFrame):
        for col in ("daily_return", "ret", "return", "value"):
            if col in x.columns:
                arr = x[col].to_numpy()
                break
        else:
            arr = x[x.columns[-1]].to_numpy()
    elif isinstance(x, pl.Series):
        arr = x.to_numpy()
    else:
        arr = np.asarray(x).flatten()
    return arr.astype(np.float64, copy=False)


def _coerce_returns(x: np.ndarray | pl.Series | pl.DataFrame) -> np.ndarray:
    arr = _as_return_array(x)
    arr = arr[np.isfinite(arr)]
    # Engine-mode parquets often carry leading zero rows from bars before the
    # first signal. Including them dilates uncertainty by underestimating
    # variance and overstating effective sample size.
    if arr.size > 0:
        nonzero = np.flatnonzero(arr != 0.0)
        if nonzero.size > 0:
            arr = arr[nonzero[0] :]
    return arr


def _to_float(v: Any) -> float:
    try:
        f = float(v)
    except (TypeError, ValueError):
        return float("nan")
    if not np.isfinite(f):
        return float("nan")
    return f


def load_daily_returns(case_study: str, backtest_hash: str) -> np.ndarray | None:
    """Load persisted daily returns for a backtest hash; None if missing."""
    from utils.paths import get_case_study_dir

    path = (
        get_case_study_dir(case_study)
        / "run_log"
        / "backtest"
        / backtest_hash
        / "daily_returns.parquet"
    )
    if not path.exists():
        return None
    df = pl.read_parquet(path)
    return _coerce_returns(df)


def _normalized_timestamp(dtype) -> pl.Expr:
    """Return an expression casting a daily-returns ``timestamp`` to ``Datetime("us")``, tz-naive.

    Daily-returns parquets across stages and case studies write this column with inconsistent
    dtypes - ``Date`` for monthly-rebalance aggregations, ``Datetime[ms]`` for some engine paths,
    ``Datetime[us]`` for others. Polars refuses to join or compare across them, so anything that
    puts two of these frames side by side has to normalize first.

    This lived only inside :func:`_align_variants_on_timestamp`, which meant a caller that joined
    two frames itself got the raw dtypes and a comparison error. Defining it once and applying it
    where the frame is loaded is what makes every caller safe rather than only that one.
    """
    expr = pl.col("timestamp")
    if dtype == pl.Date:
        return expr.cast(pl.Datetime("us"))
    if isinstance(dtype, pl.Datetime):
        if getattr(dtype, "time_zone", None) is not None:
            expr = expr.dt.replace_time_zone(None)
        if getattr(dtype, "time_unit", "us") != "us":
            expr = expr.cast(pl.Datetime("us"))
    return expr


def load_daily_returns_with_timestamp(case_study: str, backtest_hash: str) -> pl.DataFrame | None:
    """Load persisted daily returns as a (timestamp, ret) frame.

    Cohort selection statistics that pass ``correlation_method`` to
    :func:`ml4t.diagnostic.evaluation.stats.deflated_sharpe_ratio` need
    an equal-length N×K matrix across variants; this helper preserves the
    timestamp so the caller can inner-join on it before stacking.

    Returns ``None`` if the parquet is missing. Unlike :func:`load_daily_returns`,
    leading zero rows are NOT stripped here — that strip happens after
    cross-variant alignment (otherwise variants land on different windows).
    """
    from utils.paths import get_case_study_dir

    path = (
        get_case_study_dir(case_study)
        / "run_log"
        / "backtest"
        / backtest_hash
        / "daily_returns.parquet"
    )
    if not path.exists():
        return None
    df = pl.read_parquet(path)
    ret_col = next(
        (c for c in ("daily_return", "ret", "return", "value") if c in df.columns),
        df.columns[-1],
    )
    if "timestamp" not in df.columns:
        return None
    return df.select(
        _normalized_timestamp(df.schema["timestamp"]).alias("timestamp"),
        pl.col(ret_col).cast(pl.Float64).alias("ret"),
    ).drop_nulls()


def stored_cohort_members(members_json: str | None) -> list[str] | None:
    """The members a stored ``cohort_metrics`` row says its correction covers.

    ``None`` for a row written before the members were persisted, which is a different
    state from a cohort that is empty: the first cannot be verified at all, the second
    verifies trivially. Every ``cohort_metrics`` row in the fleet on 2026-09-07 - 109 of
    them across five registries - was in the first state.
    """
    if members_json is None:
        return None
    members = json.loads(members_json)
    return sorted(str(member) for member in members)


def cohort_membership_diff(
    stored: Iterable[str], in_hand: Iterable[str]
) -> tuple[list[str], list[str]]:
    """``(missing, extra)`` between the members a stored row covers and a cohort in hand.

    ``missing`` is stored and absent from the cohort the reader assembled; ``extra`` is
    assembled and not covered by the correction. Empty lists mean the stored correction is
    about exactly this cohort.

    This is what a digest comparison cannot do. Two cohorts of the same size with one
    member swapped have different digests and the same ``k_variants``, so the digest
    establishes that they differ and says nothing about how - a swapped member and a
    changed selection rule look identical through it.
    """
    stored_set = {str(member) for member in stored}
    in_hand_set = {str(member) for member in in_hand}
    return sorted(stored_set - in_hand_set), sorted(in_hand_set - stored_set)


def _align_variants_on_timestamp(
    returns_by_hash: dict[str, pl.DataFrame],
) -> tuple[np.ndarray, list[str]] | None:
    """Inner-join per-hash return frames on timestamp; return (T×K matrix, hashes).

    Variants with too-few observations after alignment are dropped. Returns
    ``None`` if fewer than 2 variants survive or fewer than 4 timestamps remain.

    Every frame is normalized through :func:`_normalized_timestamp` before joining, because
    the polars inner-join refuses to match across the dtypes these parquets carry.
    """
    frames: dict[str, pl.DataFrame] = {}
    for name, frame in returns_by_hash.items():
        if frame is None or frame.is_empty():
            continue
        if "timestamp" not in frame.columns or "ret" not in frame.columns:
            continue
        # Still normalized here as well as in the loader: `returns_by_hash` may hold frames a
        # caller assembled itself. One helper, so the two cannot drift apart.
        ts_expr = _normalized_timestamp(frame.schema["timestamp"])
        frames[name] = frame.select(
            ts_expr.alias("timestamp"),
            pl.col("ret").cast(pl.Float64).alias(name),
        )
    if len(frames) < 2:
        return None
    names = list(frames.keys())
    joined = frames[names[0]]
    for name in names[1:]:
        joined = joined.join(frames[name], on="timestamp", how="inner")
    if joined.height < 4:
        return None
    matrix = joined.select(names).to_numpy().astype(np.float64, copy=False)
    finite_rows = np.isfinite(matrix).all(axis=1)
    matrix = matrix[finite_rows]
    if matrix.shape[0] < 4:
        return None
    return matrix, names


def _distinct_trials(
    matrix: np.ndarray, names: list[str], *, keep: int
) -> tuple[np.ndarray, list[str]]:
    """Collapse columns carrying the same result series into one trial each.

    A regularisation grid that runs past the point where the penalty stops binding
    submits several configurations and produces one result. On
    ``nasdaq100_microstructure``/``fwd_dir_15m``, L1 logistic at C=10 and C=100 agree to
    six decimals in log loss and to the coefficient in sparsity (195 of 198 non-zero),
    under liblinear and saga independently: at C >= 10 the penalty is barely binding and
    a further tenfold weakening changes the fitted model not at all. Each configuration
    still occupies a training row, a prediction set and a backtest.

    The multiple-testing adjustment asks how many chances the selection had to find a
    high Sharpe by luck. A configuration that reproduces another one's series supplies no
    chance, so counting it inflates K, inflates the expected maximum a zero-skill cohort
    would reach, and understates the deflated Sharpe underneath.

    Equality is exact and needs no tolerance: two configurations that converged to the
    same fitted model emit the same predictions and so the same returns, bit for bit. A
    read-only scan of the fleet on 2026-09-07 found 375 of 15152 registered backtests
    reproducing another one in the same ``(stage, label)``, in runs of up to six (``ols``
    through ``ridge_a10.0`` on sp500_equity_option_analytics and us_firm_characteristics);
    rounding to twelve decimals first found exactly the same 375, so a tolerance would
    widen the rule without finding anything. Two series that differ only in the sign of a
    zero are left as two, which errs toward counting a trial that is not one.

    ``keep`` is the column that must survive as its group's representative — the cohort
    leader, so every adjusted statistic still refers to the row the caller reports.
    Returns the collapsed matrix and its names, in first-seen group order.
    """
    groups: dict[bytes, list[int]] = {}
    for column in range(matrix.shape[1]):
        key = np.ascontiguousarray(matrix[:, column]).tobytes()
        groups.setdefault(key, []).append(column)
    representatives = [keep if keep in members else members[0] for members in groups.values()]
    return matrix[:, representatives], [names[column] for column in representatives]


def cohort_member_digest(hashes: Iterable[str]) -> str:
    """Identify a cohort by its members rather than by how many it has.

    Order-independent, so the digest does not depend on how the caller happened to
    assemble the cohort, and duplicates collapse - a hash is in the cohort or it is not.
    """
    unique = sorted(set(str(h) for h in hashes))
    return hashlib.sha256("\n".join(unique).encode()).hexdigest()


def compute_cohort_metrics(
    returns_by_hash: dict[str, pl.DataFrame],
    *,
    periods_per_year: int,
    baseline_returns: pl.DataFrame | np.ndarray | None = None,
    fold_returns_by_hash: dict[str, np.ndarray] | None = None,
    rademacher_n_simulations: int = 2000,
    rademacher_seed: int = 0,
) -> dict[str, Any]:
    """Compute the full cohort selection-bias bundle for a set of variants.

    Returns a flat dict matching the ``cohort_metrics`` table schema (minus
    identity columns ``cohort_type / stage / label / family``, which the
    caller adds). Empty dict if alignment fails or too few variants survive.

    The ``leader_hash`` value emitted in the result is one of the dict keys
    of ``returns_by_hash`` — the contract is that those keys ARE the
    ``backtest_runs.backtest_hash`` strings used as the natural identifier
    everywhere downstream. The ``cohort_metrics`` table has a
    ``leader_hash REFERENCES backtest_runs(backtest_hash) NOT NULL`` FK,
    so passing non-hash dict keys (synthetic variant names, family labels,
    …) will fail at insert with a foreign-key violation. Callers compose
    the dict from ``load_daily_returns_with_timestamp(case_study, hash)``
    keyed on the backtest hash — do not key on family/method names.

    Every estimator below runs on the cohort's *distinct results*, not on the
    configurations submitted: a regularisation grid that runs past the point where the
    penalty stops binding hands in several configurations and produces one series, and
    each of those occupies a training row, a prediction set and a backtest. The
    correction is for how many chances the selection had to find a high Sharpe by luck,
    and a repeat supplies none. ``k_variants`` is that count, because it is the K a
    notebook prints beside a deflated Sharpe and it has to be the K that deflated it;
    ``k_variants_submitted`` records the configurations, and ``member_digest`` still
    covers every aligned member. See :func:`_distinct_trials`.

    Estimators
    ----------
    - Raw-K DSR (no correlation adjustment) via
      :func:`ml4t.diagnostic.evaluation.stats.deflated_sharpe_ratio`.
    - MP-K DSR (``correlation_method="marchenko_pastur"``).
    - ER-K DSR (``correlation_method="effective_rank"``).
    - Rademacher Adjusted Sharpe (RAS) — lower bound on leader Sharpe.
    - White's Reality Check vs ``baseline_returns`` (optional).
    - Probability of Backtest Overfitting (CSCV) on
      ``fold_returns_by_hash`` (optional). The per-fold Sharpe matrix is
      partitioned into all C(S, S/2) IS/OOS half-fold combinations; the
      library's :func:`compute_pbo` then operates on the IS / OOS mean
      Sharpe matrices.

    Parameters
    ----------
    returns_by_hash
        Mapping ``backtest_hash → pl.DataFrame[timestamp, ret]`` (use
        :func:`load_daily_returns_with_timestamp`). Dict keys MUST be
        registered ``backtest_runs.backtest_hash`` values — see contract
        note above.
    periods_per_year
        Annualization factor (use :func:`periods_per_year_from_setup`).
    baseline_returns
        If provided, used as Reality Check benchmark. Same frame layout
        as a variant, or a numpy array of returns aligned to the variant
        intersection.
    fold_returns_by_hash
        Mapping ``backtest_hash → per-fold Sharpe ratios (1D)``. All
        variants must share fold cardinality. Skipped if not provided.
    """
    from ml4t.diagnostic.evaluation.stats import (
        RASResult,
        compute_min_trl,
        deflated_sharpe_ratio,
        effective_number_of_trials,
        rademacher_complexity,
        ras_sharpe_adjustment,
    )

    aligned = _align_variants_on_timestamp(returns_by_hash)
    if aligned is None:
        return {}
    matrix, names = aligned
    n_periods, k_variants = matrix.shape
    if k_variants < 2:
        return {}

    sharpes = _sharpe_per_column(matrix, periods_per_year)
    if np.all(np.isnan(sharpes)):
        return {}
    leader_idx = _leader_index(sharpes, names)
    leader_hash = names[leader_idx]
    leader_arr = matrix[:, leader_idx]

    # Every selection adjustment below runs on the distinct results rather than on the
    # submitted configurations: a grid that saturates hands in several configurations
    # and produces one series, and the correction is for how many chances the selection
    # had. `k_variants` and `member_digest` still describe the cohort's membership, which
    # is what a reader matches a stored correction against.
    trial_matrix, trial_names = _distinct_trials(matrix, names, keep=leader_idx)
    k_trials = trial_matrix.shape[1]
    if k_trials < 2:
        # The same rule as the `k_variants < 2` guard above, applied to what the cohort
        # actually tried. A cohort whose configurations all produced one series ran no
        # selection, so there is nothing to correct for and no correction to report. Left
        # unguarded it writes a half-row: `dsr_raw` computed at K=1, which is the
        # undeflated Sharpe, beside NULL MP and ER columns whose estimators refuse a
        # single strategy - and a reader takes the first for a corrected figure.
        return {}
    trial_leader_idx = trial_names.index(leader_hash)
    trial_sharpes = _sharpe_per_column(trial_matrix, periods_per_year)

    out: dict[str, Any] = {
        "leader_hash": leader_hash,
        # The trials every adjustment below faced, which is what a K printed beside a
        # deflated Sharpe has to be. Equal to the membership unless the cohort holds
        # configurations that produced the same series; see `_distinct_trials`.
        "k_variants": int(k_trials),
        # The configurations the cohort holds. Larger than `k_variants` exactly when a
        # grid saturated, and the two together are what says by how much.
        "k_variants_submitted": int(k_variants),
        # `names` is the cohort the correction below is actually computed over, after
        # alignment has dropped whatever could not be aligned. The digest identifies it
        # cheaply and the members are the fact a reader checks against: the digest is
        # one-way, so on its own it turns verification into a replay of every selection
        # rule that was in force when the row was written. Both come from `names` here so
        # they cannot describe different cohorts.
        "member_digest": cohort_member_digest(names),
        "members_json": json.dumps(sorted(names)),
        "periods_per_year": float(periods_per_year),
        "leader_sharpe": float(sharpes[leader_idx]),
    }

    # Per-estimator failures are caught narrowly and surfaced via warnings
    # so a regression in the library API or a degenerate input shape shows
    # up as a one-line emission rather than a silent NULL in the registry.
    _ESTIMATOR_ERRORS = (ValueError, TypeError, np.linalg.LinAlgError, ZeroDivisionError)

    # Leader Sortino + MinTRL
    try:
        out["leader_sortino"] = float(_sortino(leader_arr, periods_per_year))
    except _ESTIMATOR_ERRORS as exc:
        warnings.warn(f"leader_sortino failed for {leader_hash}: {exc}", stacklevel=2)
        out["leader_sortino"] = None
    try:
        mtrl = compute_min_trl(leader_arr, periods_per_year=periods_per_year)
        out["leader_min_trl"] = float(mtrl.min_trl)
    except _ESTIMATOR_ERRORS as exc:
        warnings.warn(f"leader_min_trl failed for {leader_hash}: {exc}", stacklevel=2)
        out["leader_min_trl"] = None

    # Effective trials — MP and ER
    try:
        et_mp = effective_number_of_trials(trial_matrix, method="marchenko_pastur")
        out["n_trials_effective_mp"] = float(et_mp.k_eff)
    except _ESTIMATOR_ERRORS as exc:
        warnings.warn(f"n_trials_effective_mp failed for {leader_hash}: {exc}", stacklevel=2)
        out["n_trials_effective_mp"] = None
    try:
        et_er = effective_number_of_trials(trial_matrix, method="effective_rank")
        out["n_trials_effective_er"] = float(et_er.k_eff)
    except _ESTIMATOR_ERRORS as exc:
        warnings.warn(f"n_trials_effective_er failed for {leader_hash}: {exc}", stacklevel=2)
        out["n_trials_effective_er"] = None

    # DSR — raw, MP, ER (three calls; library handles K correctly per method)
    arr_list = [trial_matrix[:, i] for i in range(k_trials)]
    methods: tuple[
        tuple[str, Literal["marchenko_pastur", "effective_rank"] | None],
        ...,
    ] = (
        ("raw", None),
        ("mp", "marchenko_pastur"),
        ("er", "effective_rank"),
    )
    for suffix, method in methods:
        try:
            if method is not None:
                dsr = deflated_sharpe_ratio(