Aprendizado de máquina duplo em séries temporais e refutação por blocos

# causal.py


```py
"""Shared causal inference utilities for Ch15 notebooks and case study DML.

Provides:
- block_permute(): Block permutation preserving autocorrelation
- manual_dml_timeseries(): Walk-forward DML with embargo
- run_dml_analysis(): Full DML pipeline (naive + DML + refutation)

Used by teaching notebooks (02-04, 07) and case study 09_causal_dml.py.
"""

from __future__ import annotations

import os

# HistGradientBoostingRegressor uses OMP-parallel histogram construction whose
# floating-point reduction order is non-deterministic across threads, so the placebo
# loop is only bit-reproducible at a fixed thread count.
#
# Setting OMP_NUM_THREADS here does NOT achieve that on its own. Every notebook imports
# case_studies.research (and through it sklearn) before this module, so the OpenMP and
# OpenBLAS runtimes have already read the variable by the time this line runs: measured
# with threadpoolctl, the notebook import order leaves the openmp pool at 16 and openblas
# at 24 despite os.environ reporting 1. The env pin is kept because it does work when this
# module is imported first, but DML_THREAD_LIMIT below is the mechanism that holds, and the
# limit is recorded in the resolved specification so two runs at different counts cannot
# share an identity.
os.environ.setdefault("OMP_NUM_THREADS", "1")

# Fixed rather than derived from the host: a value like -1 varies with the machine, so a
# result would not be identity-stable across the readers' hardware.
DML_THREAD_LIMIT = 1
# 1 -> 2 on 2026-09-10: the block-permutation refutation moved from comparing raw effects
# to comparing HAC t-statistics. That changes a registered
# value, so it has to move the identity - and causal rows have no migration path, so every
# causal row in every case study refits rather than being re-keyed. That is the intended
# cost: a stored refutation_p computed on raw effects is anti-conservative, always in the
# direction of "Passes", and there is no arithmetic that converts one into the other
# because the t-scale draws were never recorded.
CAUSAL_RUNNER_VERSION = 2

import importlib.metadata
import json
import platform
import re
import time
from dataclasses import dataclass
from datetime import UTC, datetime
from typing import TYPE_CHECKING, Any

import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
from sklearn.ensemble import HistGradientBoostingRegressor
from statsmodels.regression.linear_model import OLS
from threadpoolctl import threadpool_limits

from utils.modeling import RANDOM_SEED, seed_everything

if TYPE_CHECKING:
    from case_studies.research.workspace import Study


from case_studies.utils.preview_fields import DML_PREVIEW_FIELDS as _DML_PREVIEW_FIELDS
from case_studies.utils.warning_policy import warn_the_reader


@dataclass(frozen=True)
class DMLResearchContext:
    analysis: pd.DataFrame
    treatment_col: str
    outcome_col: str
    confounder_cols: tuple[str, ...]
    time_col: str
    entity_col: str
    n_folds: int
    embargo: int
    n_placebo: int
    block_size: int
    seed: int
    horizon: int
    expected_step: pd.Timedelta
    nuisance_params: dict[str, Any]
    runtime_provenance: dict[str, Any]


def observation_step(frame: Any, date_col: str = "timestamp") -> pd.Timedelta:
    """Measure the spacing between consecutive decision times in *frame*.

    Returns the most common gap between distinct sorted values of ``date_col``,
    which is the observation grid the panel is actually recorded on. Sessions,
    weekends and holidays introduce larger gaps; taking the mode rather than the
    minimum or the mean makes those irrelevant.

    Accepts a Polars or pandas frame, or anything exposing the column through
    ``__getitem__``.
    """
    column = frame[date_col]
    values = pd.Series(column.to_list() if hasattr(column, "to_list") else list(column))
    stamps = pd.to_datetime(values).drop_duplicates().sort_values()
    if len(stamps) < 2:
        raise ValueError(f"{date_col!r} has fewer than two distinct values; no grid to measure")
    gaps = stamps.diff().dropna()
    return pd.Timedelta(gaps.mode().iloc[0])


def embargo_from_buffer(
    label_buffer: str,
    *,
    periods_per_year: int | None = None,
    observed_step: str | pd.Timedelta | None = None,
) -> int:
    """Convert a label buffer string to an integer embargo period count.

    The embargo is counted in *observation periods*, so converting a duration
    into one requires knowing how long an observation period is. Pass
    ``observed_step`` — from :func:`observation_step` on the frame being
    analysed — and the conversion is exact: the number of periods spanning the
    buffer, rounded up, at least one.

    Without ``observed_step`` the conversion falls back to a fixed assumption
    about bar size per unit, which is correct only when that assumption holds:

    - D: one period per ``value`` days, correct on a daily grid
    - H/h: the number of ``value``-hour bars in one day, so "8H" gives a
      one-day embargo on 8-hour bars
    - M: ``value`` monthly groups when ``periods_per_year=12``, else
      ``value`` months x 21 trading days
    - T/min: the number of ``value``-minute spans in 15 minutes, which assumes
      the panel is recorded in 15-minute bars

    That last assumption is the one that bites, because it is wrong by exactly
    the ratio between the assumed bar and the real one, and nothing in the result
    shows it. A "15min" buffer resolves to a single period whatever the panel is
    recorded at, so on a one-minute grid it yields a one-minute embargo against a
    fifteen-minute label. Pass ``observed_step`` on any sub-daily panel rather
    than relying on the declared bar size, which can disagree with the artifacts.

    A month buffer has no fixed length and is rejected when ``observed_step`` is
    supplied; use the ``periods_per_year`` branch for it.
    """
    import math
    import re

    if observed_step is not None:
        step = pd.Timedelta(observed_step)
        if step <= pd.Timedelta(0):
            raise ValueError(f"observed_step must be positive, got {observed_step!r}")
        # pandas deprecated the uppercase hour alias; the buffers are authored by
        # hand in setup.yaml and still use it.
        normalized = re.sub(r"(?<=\d)H\b", "h", label_buffer.strip())
        if re.match(r"\d+\s*M\b", normalized):
            raise ValueError(
                f"A month buffer ({label_buffer!r}) has no fixed length, so it cannot be "
                f"divided by an observation step. Use the periods_per_year branch by "
                f"omitting observed_step."
            )
        span = pd.Timedelta(normalized)
        # The floor of one period is there so a buffer shorter than a bar still
        # embargoes the bar it resolves inside. A zero-length buffer is not that
        # case: it declares that the label resolves on its own row's timestamp, so
        # there is nothing to embargo, and the floor would invent a gap the case
        # study did not ask for. The `periods_per_year` branch below answers zero
        # for the same declaration, and the two must not disagree.
        if span == pd.Timedelta(0):
            return 0
        return max(1, math.ceil(span / step))

    match = re.match(r"(\d+)(D|H|h|M|T|min)", label_buffer.strip())
    if not match:
        raise ValueError(f"Cannot parse label_buffer: {label_buffer}")
    value, unit = int(match.group(1)), match.group(2)
    # The per-unit conversions below divide by the value and are built eagerly, so
    # every unit raised on a zero, including the D branch that would have returned it
    # unchanged had it been reached. A case study whose outcome resolves on the row's
    # own timestamp declares exactly that, for example `labels.horizons` in
    # us_firm_characteristics.
    if value == 0:
        return 0
    return {
        "D": value,
        "H": max(1, 24 // value),
        "h": max(1, 24 // value),
        "M": value if periods_per_year == 12 else value * 21,
        "T": max(1, value // 15),
        "min": max(1, value // 15),
    }[unit]


# A segment boundary is a hole in the observation series, not a step the calendar
# always takes. Splitting wherever the gap differs from the cadence at all cut a
# daily series at every weekend, leaving five-row segments that could not hold two
# blocks of any useful size; the short-segment path then shuffled them. Four
# cadences clears a weekend (three) and a long weekend (four) on a daily series and
# still catches a real hole.
GAP_TOLERANCE_CADENCES = 4


def block_permute(
    arr: np.ndarray,
    block_size: int,
    rng: np.random.Generator | None = None,
    groups: np.ndarray | None = None,
    units: np.ndarray | None = None,
    expected_step: str | pd.Timedelta | None = None,
    gap_tolerance: str | pd.Timedelta | None = None,
) -> np.ndarray:
    """Permute array in blocks to preserve autocorrelation structure.

    Essential for refutation tests on time series data. Random permutation
    destroys autocorrelation, making placebo tests too easy to pass.

    Parameters
    ----------
    arr : array-like
        Array to permute.
    block_size : int
        Size of blocks to preserve.
    rng : np.random.Generator, optional
        Random number generator for reproducibility.
    groups : array-like, optional
        Ordered decision time for each row. For a single time series, this
        validates that each row is one decision time.
    units : array-like, optional
        Panel entity for each row. When supplied with ``groups``, treatment is
        block-permuted within each entity, so ``block_size`` counts that
        entity's ordered decision times rather than flattened panel rows.

    Returns
    -------
    np.ndarray
        Block-permuted array.
    """
    arr = np.asarray(arr)
    if rng is None:
        rng = np.random.default_rng()

    segments = _permutation_segments(len(arr), groups, units, expected_step, gap_tolerance)
    result = np.array(arr, copy=True)
    for idx in segments:
        result[idx] = _permute_one_segment(arr[idx], block_size, rng)
    return result


def _permutation_segments(
    n: int,
    groups: np.ndarray | None,
    units: np.ndarray | None,
    expected_step: str | pd.Timedelta | None,
    gap_tolerance: str | pd.Timedelta | None,
) -> list[np.ndarray]:
    """The maximal stretches of rows a block permutation may reorder within.

    One entity's uninterrupted run of decision times is one segment. Blocks never cross
    a segment boundary, because the rows on either side are not adjacent in time - they
    belong to different entities, or to the same entity either side of a gap.

    This is the single definition of where the series is cut. ``block_permute`` permutes
    within these segments and ``_immobile_masks`` counts what they leave standing, so
    the two can never disagree about the segmentation.
    """
    if units is not None:
        if groups is None:
            raise ValueError("groups are required when units are supplied")
        group_arr = np.asarray(groups)
        unit_arr = np.asarray(units)
        if len(group_arr) != n or len(unit_arr) != n:
            raise ValueError("groups and units must have the same length as arr")
        segments: list[np.ndarray] = []
        for unit in pd.unique(unit_arr):
            idx = np.flatnonzero(unit_arr == unit)
            unit_groups = group_arr[idx]
            if len(unit_groups) > 1 and np.any(unit_groups[1:] <= unit_groups[:-1]):
                raise ValueError("groups must be strictly increasing within each unit")
            segments.extend(
                idx[bounds]
                for bounds in _contiguous_runs(unit_groups, expected_step, gap_tolerance)
            )
        return segments

    if groups is not None:
        group_arr = np.asarray(groups)
        if len(group_arr) != n:
            raise ValueError("groups must have the same length as arr")
        if len(np.unique(group_arr)) != n:
            raise ValueError("units are required when decision times contain multiple rows")
        if n > 1 and np.any(group_arr[1:] <= group_arr[:-1]):
            raise ValueError("groups must be strictly increasing")
        positions = np.arange(n)
        return [
            positions[bounds]
            for bounds in _contiguous_runs(group_arr, expected_step, gap_tolerance)
        ]

    return [np.arange(n)]


def _contiguous_runs(
    group_arr: np.ndarray,
    expected_step: str | pd.Timedelta | None,
    gap_tolerance: str | pd.Timedelta | None,
) -> list[slice]:
    """Split one entity's ordered decision times wherever the gap exceeds the tolerance."""
    n = len(group_arr)
    if expected_step is None or n <= 1:
        return [slice(0, n)]
    cadence = pd.Timedelta(expected_step)
    tolerance = (
        pd.Timedelta(gap_tolerance)
        if gap_tolerance is not None
        else GAP_TOLERANCE_CADENCES * cadence
    )
    timestamps = pd.to_datetime(group_arr, utc=True)
    steps = np.asarray(timestamps[1:] - timestamps[:-1])
    boundaries = np.flatnonzero(steps > tolerance) + 1
    if not boundaries.size:
        return [slice(0, n)]
    starts = np.r_[0, boundaries]
    stops = np.r_[boundaries, n]
    return [slice(int(a), int(b)) for a, b in zip(starts, stops, strict=True)]


def _permute_one_segment(
    values: np.ndarray, block_size: int, rng: np.random.Generator
) -> np.ndarray:
    """Reorder whole blocks within one uninterrupted stretch."""
    n = len(values)
    n_blocks = n // block_size
    if n_blocks < 2:
        # Not enough room for two blocks, so there is no permutation to make at this
        # block size. Returning the segment intact preserves the dependence the
        # caller asked to keep; shuffling it would destroy exactly that, which is
        # what the old `rng.permutation(arr)` did to every weekend-bounded segment
        # of a daily series. A caller that permutes nothing at all is caught by
        # `_assert_placebo_permutation_possible`, not here: inside a panel, one short unit staying
        # put while the others move is correct.
        return np.array(values, copy=True)

    pieces = [
        values[idx * block_size : (idx + 1) * block_size] for idx in rng.permutation(n_blocks)
    ]

    # The trailing rows that do not fill a whole block stay where they are. This is a
    # property of block permutation itself, not of the data: it happens to any segment
    # whose length is not a multiple of the block size, however long the segment is.
    remainder_start = n_blocks * block_size
    if remainder_start < n:
        pieces.append(values[remainder_start:])

    return np.concatenate(pieces)


def _immobile_masks(
    n: int, block_size: int, segments: list[np.ndarray]
) -> tuple[np.ndarray, np.ndarray]:
    """Which rows no draw can move, split by the two unrelated reasons.

    The first mask is rows in a segment too short to hold two blocks. Those are frozen
    because of how the data is shaped against the block size, they carry their observed
    values into every placebo, and a large share of them makes the refutation
    uninformative - which is worth telling the caller about.

    The second is the trailing remainder of a segment whose length is not a multiple of
    the block size. Those rows also never move, but for a reason intrinsic to block
    permutation that no choice of gap tolerance can remove, and shrinking the block size
    to chase them is the wrong response. They are reported, not warned about.
    """
    short = np.zeros(n, dtype=bool)
    remainder = np.zeros(n, dtype=bool)
    for idx in segments:
        n_blocks = len(idx) // block_size
        if n_blocks < 2:
            short[idx] = True
        elif n_blocks * block_size < len(idx):
            remainder[idx[n_blocks * block_size :]] = True
    return short, remainder


def _walk_forward_indices(
    n_rows: int,
    n_folds: int,
    embargo: int,
    groups: np.ndarray | None = None,
) -> list[tuple[np.ndarray, np.ndarray]]:
    """Build expanding-window folds in rows or complete decision-time groups."""
    if groups is None:
        fold_size = n_rows // (n_folds + 1)
        if fold_size == 0:
            raise ValueError(
                f"{n_folds}-fold walk-forward needs at least {n_folds + 1} rows, got {n_rows}."
            )
        folds = []
        for fold in range(n_folds):
            train_end = (fold + 1) * fold_size
            test_start = train_end + embargo
            test_end = min(test_start + fold_size, n_rows)
            folds.append((np.arange(0, train_end), np.arange(test_start, test_end)))
        return folds

    group_arr = np.asarray(groups)
    if len(group_arr) != n_rows:
        raise ValueError("groups must have the same length as the input arrays")
    group_starts = np.flatnonzero(np.r_[True, group_arr[1:] != group_arr[:-1]])
    ordered_groups = group_arr[group_starts]
    if len(np.unique(ordered_groups)) != len(ordered_groups) or (
        len(ordered_groups) > 1 and np.any(ordered_groups[1:] < ordered_groups[:-1])
    ):
        raise ValueError("groups must be sorted and contiguous")

    # A fold is `len(ordered_groups) // (n_folds + 1)` decision times wide, and integer
    # division makes that 0 whenever the panel holds fewer complete decision times than
    # folds. Every train and test slice is then empty, every fold is skipped downstream,
    # and the estimate comes back NaN over zero observations with nothing raised - the
    # caller's guard counts ROWS (`(n_folds + 1) * 50 + n_folds * embargo`), which a wide
    # panel clears on a handful of dates. Measured on a 10,000-row cap of
    # us_firm_characteristics/09_causal_dml: 8,180 rows, 4 complete decision months,
    # 5 folds, `4 // 6 == 0`, full summary printed, DML effect nan.
    #
    # The count that has to clear the geometry is decision times, so it is checked here,
    # at the one place the geometry is built, rather than in each caller.
    fold_size = len(ordered_groups) // (n_folds + 1)
    if fold_size == 0:
        raise ValueError(
            f"{n_folds}-fold walk-forward needs at least {n_folds + 1} complete decision "
            f"times, got {len(ordered_groups)}. On a panel the fold geometry is sized in "
            f"decision times, not rows, so a row-count minimum does not constrain it."
        )
    folds = []
    for fold in range(n_folds):
        train_end = (fold + 1) * fold_size
        test_start = train_end + embargo
        test_end = min(test_start + fold_size, len(ordered_groups))
        train_groups = ordered_groups[:train_end]
        test_groups = ordered_groups[test_start:test_end]
        folds.append(
            (
                np.flatnonzero(np.isin(group_arr, train_groups)),
                np.flatnonzero(np.isin(group_arr, test_groups)),
            )
        )
    return folds


def manual_dml_timeseries(
    Y: np.ndarray,
    T: np.ndarray,
    X: np.ndarray,
    n_folds: int = 5,
    embargo: int = 21,
    model_y=None,
    model_t=None,
    return_residuals: bool = False,
    hac_maxlags: int | None = None,
    horizon: int | None = None,
    groups: np.ndarray | None = None,
    thread_limit: int = DML_THREAD_LIMIT,
) -> dict:
    """Walk-forward DML with embargo for temporal data.

    Follows Chernozhukov et al. (2017) and de Prado (2018):
    1. Split data into K temporal folds (not random)
    2. For each fold, train on earlier data, predict on later
    3. Embargo gap between train and test prevents autocorrelation leakage
    4. HAC standard errors account for residual autocorrelation

    Parameters
    ----------
    Y : array
        Outcome variable.
    T : array
        Treatment variable.
    X : array
        Confounder matrix.
    n_folds : int
        Number of temporal folds.
    embargo : int
        Gap periods between train and test sets.
    model_y, model_t : sklearn estimator, optional
        Nuisance models for E[Y|X] and E[T|X].
    return_residuals : bool
        If True, include residual arrays in result dict.
    hac_maxlags : int or None
        HAC (Newey-West) bandwidth. If given, used verbatim. If None, resolved
        from `horizon` (see below).
    horizon : int or None
        Label horizon in observation periods. Overlapping h-period forward
        returns induce MA(h-1) structure, so the HAC bandwidth must satisfy
        L >= h - 1. When `hac_maxlags` is None, the bandwidth is
        `max(horizon - 1, cube-root-of-n)` if `horizon` is given, else the
        cube-root rule alone. Pass this for any overlapping label of horizon
        >= ~10 periods, or the standard error is understated and the
        t-statistic overstated.
    groups : array-like or None
        Ordered decision-time group for each observation. For panel data,
        supply the timestamp column so folds and embargoes operate on complete
        decision times rather than arbitrary rows.

    Returns
    -------
    dict
        Keys: theta, se_iid, se_hac, t_stat_iid, t_stat_hac, p_value_hac,
        n_obs, n_periods, hac_maxlags, covariance_type. If return_residuals:
        also Y_res, T_res.

        `covariance_type` is "driscoll_kraay" when groups were supplied,
        "newey_west" when they were not, and "failed" when the robust
        covariance did not produce a usable one. On "failed", `se_hac`,
        `t_stat_hac` and `p_value_hac` are NaN and `hac_maxlags` is 0;
        `se_iid` still carries the HC0 standard error under its own name.
        The NaN is deliberate: a fallback that returns HC0 under the name
        `se_hac` reports a smaller standard error and a smaller p-value than
        the correct one, and every caller that does not read
        `covariance_type` takes it for a robust result.
    """
    # Pinned here, where the nuisance models are actually fitted, so that every caller is
    # covered: run_dml_analysis, the six case-study DML stages that call it directly, and
    # the chapter-15 notebooks that call this function themselves and rely on the default
    # HistGradientBoostingRegressor. The module-level OMP_NUM_THREADS setdefault does not
    # bind for any of them, because sklearn is imported first.
    with threadpool_limits(limits=thread_limit):
        seed_everything(RANDOM_SEED)

        n = len(Y)

        # Initialize residual arrays
        Y_res = np.full(n, np.nan)
        T_res = np.full(n, np.nan)

        folds = _walk_forward_indices(n, n_folds, embargo, groups=groups)

        for train_idx, test_idx in folds:
            if len(test_idx) == 0:
                continue

            if len(train_idx) < 50 or len(test_idx) < 10:
                continue

            # Fit nuisance models on training data (clone to avoid mutation)
            from sklearn.base import clone

            _default_y = HistGradientBoostingRegressor(max_iter=50, max_depth=3, random_state=42)
            _default_t = HistGradientBoostingRegressor(max_iter=50, max_depth=3, random_state=42)
            my = clone(model_y) if model_y is not None else _default_y
            mt = clone(model_t) if model_t is not None else _default_t

            my.fit(X[train_idx], Y[train_idx])
            mt.fit(X[train_idx], T[train_idx])

            Y_res[test_idx] = Y[test_idx] - my.predict(X[test_idx])
            T_res[test_idx] = T[test_idx] - mt.predict(X[test_idx])

        # Drop observations without residuals
        valid = ~np.isnan(Y_res) & ~np.isnan(T_res)
        Y_v = Y_res[valid]
        T_v = T_res[valid]
        n_valid = len(Y_v)
        valid_groups = np.asarray(groups)[valid] if groups is not None else None
        n_periods = len(np.unique(valid_groups)) if valid_groups is not None else n_valid

        empty = {
            "theta": np.nan,
            "se_iid": np.nan,
            "se_hac": np.nan,
            "t_stat_iid": np.nan,
            "t_stat_hac": np.nan,
            "p_value_hac": np.nan,
            "n_obs": n_valid,
            "n_periods": n_periods,
            "hac_maxlags": 0,
            # Nothing was estimated on this path, so naming the estimator that would
            # have run reports a success that did not happen. The invariant the whole
            # dict now holds: `covariance_type == "failed"` exactly when `se_hac` is
            # not a robust standard error.
            "covariance_type": "failed",
        }
        if n_valid < 50:
            if return_residuals:
                empty["Y_res"] = Y_res
                empty["T_res"] = T_res
            return empty

        # Final stage: Y_res = alpha + theta * T_res + epsilon
        # Must include intercept: cross-fitting residuals may have non-zero mean
        # when training data varies across folds (expanding window).
        if hac_maxlags is None:
            auto = max(1, int(n_periods ** (1 / 3)))
            # Overlapping h-period labels need L >= h-1; the cube-root rule is
            # horizon-blind and under-lags long-horizon overlapping returns.
            hac_maxlags = max(horizon - 1, auto) if horizon else auto
            hac_maxlags = min(hac_maxlags, max(1, n_periods // 2))

        T_const = sm.add_constant(T_v)
        ols_iid = OLS(Y_v, T_const).fit()
        theta = ols_iid.params[1]

        # HC0 standard error
        se_iid = np.sqrt(ols_iid.cov_HC0[1, 1])

        # Serial-correlation-robust standard error with frequency-adaptive bandwidth.
        # Panel rows share decision times, so ordinary row-wise Newey-West treats
        # cross-sectional observations as extra time periods and understates risk.
        # Driscoll-Kraay aggregates the score by decision time and remains robust to
        # general cross-sectional dependence.
        #
        # `covariance_type` is the field a reader and a registry query use to decide what
        # `se_hac` is, so it is assigned from what happened rather than from which branch
        # was intended. It used to be set unconditionally at the bottom of this function,
        # which meant a run that fell back to HC0 still reported a successful
        # Driscoll-Kraay - and `se_hac` still carried a number, seeded from `se_iid`.
        # HC0 on an overlapping panel label is the estimator this chain exists to
        # replace, so that fallback understated the standard error in the direction of
        # significance under a name that said it was robust.
        covariance_type = "driscoll_kraay" if groups is not None else "newey_west"
        se_hac = np.nan
        if valid_groups is not None and n_periods < 2:
            # A groupsum HAC over a single decision time is not an estimator: one group
            # score, no lag structure to estimate. statsmodels does not raise on it - it
            # returns a variance at the rounding floor. Measured end to end against
            # this function on `origin/main`, at n_periods = 1 over 60 rows:
            # se_iid = 0.1237, se_hac = 3.69e-16, t_stat_hac = 2.28e15,
            # p_value_hac = 2.80e-16, covariance_type = "driscoll_kraay", and
            # `register_causal_run`'s finiteness check accepted all of it. The
            # `n_valid >= 50` guard above does not bound this, because sixty entities
            # sharing one timestamp satisfies it. Two periods is also where the
            # bandwidth cap becomes self-consistent: `max(1, n_periods // 2)` floors
            # the bandwidth at one lag even when the sample holds none.
            covariance_type = "failed"
            warn_the_reader(
                f"Driscoll-Kraay needs at least two decision times; got {n_periods} over "
                f"{n_valid} rows. Reporting se_hac, t_stat_hac and p_value_hac as NaN; "
                f"se_iid carries the HC0 standard error under its own name.",
                source="manual_dml_timeseries",
                key=("dk_periods", n_periods, n_valid),
                category=RuntimeWarning,
            )
        else:
            try:
                if valid_groups is not None:
                    time_codes = pd.factorize(valid_groups, sort=False)[0]
                    robust = ols_iid.get_robustcov_results(
                        cov_type="hac-groupsum",
                        time=time_codes,
                        maxlags=hac_maxlags,
                        use_correction="hac",
                        df_correction=False,
                    )
                else:
                    robust = ols_iid.get_robustcov_results(
                        cov_type="HAC",
                        maxlags=hac_maxlags,
                        use_correction=True,
                    )
                cov = robust.cov_params()
                variance = cov.iloc[1, 1] if hasattr(cov, "iloc") else cov[1, 1]
                if not np.isfinite(variance) or variance <= 0:
                    # Raised into the handler below rather than returned, so that a
                    # covariance which comes back unusable without raising is reported
                    # the same way as one that raises. Without this the estimator name
                    # would say driscoll_kraay while se_hac was NaN.
                    raise ValueError(
                        f"robust covariance produced a non-positive variance: {variance!r}"
                    )
                se_hac = np.sqrt(variance)
            except (np.linalg.LinAlgError, ValueError) as exc:
                # Narrow, because `except Exception` also swallowed every programming
                # error in these lines and degraded it silently to HC0. These two are
                # what a singular or ill-conditioned fit raises; anything else here is
                # a defect in this function and propagates.
                covariance_type = "failed"
                warn_the_reader(
                    f"robust covariance failed ({type(exc).__name__}: {exc}). Reporting "
                    f"se_hac, t_stat_hac and p_value_hac as NaN; se_iid carries the HC0 "
                    f"standard error under its own name.",
                    source="manual_dml_timeseries",
                    key=("dk_failed", type(exc).__name__, str(exc)),
                    category=RuntimeWarning,
                )

        if covariance_type == "failed":
            # The bandwidth that was computed is not a bandwidth that was applied.
            hac_maxlags = 0

        t_stat_hac = theta / se_hac if se_hac > 0 else np.nan
        p_value_hac = (
            2 * stats.t.sf(abs(t_stat_hac), df=max(n_periods - 2, 1))
            if not np.isnan(t_stat_hac)
            else np.nan
        )

        result = {
            "theta": theta,
            "se_iid": se_iid,
            "se_hac": se_hac,
            "t_stat_iid": theta / se_iid if se_iid > 0 else np.nan,
            "t_stat_hac": t_stat_hac,
            "p_value_hac": p_value_hac,
            "n_obs": n_valid,
            "n_periods": n_periods,
            "hac_maxlags": hac_maxlags,
            "covariance_type": covariance_type,
        }

        if return_residuals:
            result["Y_res"] = Y_res
            result["T_res"] = T_res

        return result


REFUTATION_ALPHA = 0.05


def _placebo_is_unchanged(original: np.ndarray, permuted: np.ndarray) -> bool:
    """Whether a placebo draw returned the observed treatment.

    ``np.array_equal`` calls two arrays different wherever either holds a NaN, so a
    frame the resolver's ``drop_nulls()`` never touched - which is every frame the
    case-study notebooks pass to ``run_dml_analysis`` directly - would report every
    identity draw as a real permutation. Compare the non-null positions and require
    the null positions to agree.
    """
    return not _placebo_moved_mask(original, permuted).any()


def _placebo_moved_fraction(original: np.ndarray, permuted: np.ndarray) -> float:
    """The share of comparable rows this draw actually moved.

    Separate from `_placebo_moved_mask` on purpose. That one answers "are these the same
    series at all", so a disagreement in null positions means every row counts as moved;
    here that convention would report a frozen panel as fully permuted the moment one
    draw shifted a missing value. Rows that are null on either side are simply not
    comparable, so they are left out of both the count and the denominator.
    """
    original = np.asarray(original)
    permuted = np.asarray(permuted)
    if original.shape != permuted.shape:
        return 0.0
    if not np.issubdtype(original.dtype, np.floating):
        comparable = np.ones(original.shape, dtype=bool)
    else:
        comparable = ~np.isnan(original) & ~np.isnan(permuted)
    if not comparable.any():
        return 0.0
    return float(np.mean(original[comparable] != permuted[comparable]))


def _placebo_moved_mask(original: np.ndarray, permuted: np.ndarray) -> np.ndarray:
    """Which comparable rows the permutation actually moved.

    ``True`` where the row differs, ``False`` where it sits still. Rows whose original
    value is missing are never comparable, so they are reported as unmoved rather than
    counted as evidence either way. A shape or null-position mismatch means the two are
    not the same series at all, so every row counts as moved.
    """
    original = np.asarray(original)
    permuted = np.asarray(permuted)
    if original.shape != permuted.shape:
        return np.ones(original.shape, dtype=bool)
    if not np.issubdtype(original.dtype, np.floating):
        return original != permuted
    missing = np.isnan(original)
    if not np.array_equal(missing, np.isnan(permuted)):
        return np.ones(original.shape, dtype=bool)
    moved = np.zeros(original.shape, dtype=bool)
    moved[~missing] = original[~missing] != permuted[~missing]
    return moved


# Below this many successful placebo draws the permutation test is not computed at all:
# the plus-one correction floors the empirical p at 1/(n+1), so under ten draws no data
# could produce a pass and a number would be reported that no test earned. When it is not
# computed, `refutation` stays empty and `refutation_p` is registered NULL - a missing
# measurement, which is what it is.
MIN_PLACEBO_DRAWS = 10


def placebo_request_is_on_the_boundary(n_placebo: int) -> bool:
    """Would one failed draw take the whole refutation with it?

    This is not enforced at run time, and the attempt to do so is worth recording. Nine
    tests in tests/test_causal_adapter.py request two to six draws on purpose, to
    exercise the block-span and permutation-feasibility logic without paying for a
    refutation they never read; refusing every small request turned all nine red for a
    property they are not about. The boundary is a property of a *declared reduction* -
    a config that says how a real run should be made cheap - not of every call.

    Zero means "do not refute", which is a different statement from "refute with too
    few draws to say anything".
    """
    return 0 < n_placebo < MIN_PLACEBO_DRAWS + PLACEBO_REQUEST_MARGIN


# Enough that one draw failing does not take the whole test with it. Small on purpose:
# a placebo draw is a full nuisance refit, so this is the least that makes the boundary
# unreachable by a single failure rather than a comfortable cushion.
PLACEBO_REQUEST_MARGIN = 5


def _assert_placebo_permutation_possible(
    unchanged_draws: int, n_draws: int, block_size: int, short_segment_fraction: float = 0.0
) -> None:
    """A refutation whose every placebo equals the observed treatment measures nothing.

    ``block_permute`` leaves a segment intact when it cannot hold two blocks of the
    requested size, which is the right thing to do to one short unit in a panel. If it
    happens to *every* segment - the block is larger than the longest uninterrupted
    stretch of observations - then the "permuted" treatment is the observed treatment,
    every placebo effect equals the observed effect, and the refutation reports p = 1
    while looking like it ran. Fail instead of publishing that.

    The test is over the whole set of draws, not each one. ``rng.permutation(n_blocks)``
    can return the identity by chance - one time in two at two blocks, one in six at
    three - so failing on a single unchanged draw aborts runs that are structurally
    fine, with a message asserting something false about the data. Every draw coming
    back unchanged is the structural condition; the chance of that happening to a
    series that can be permuted falls off as the draws multiply.
    """
    if n_draws and unchanged_draws == n_draws:
        raise ValueError(
            f"block permutation with block_size={block_size} left the treatment "
            f"unchanged on all {n_draws} placebo draws: no uninterrupted segment of "
            "the series holds two blocks of that size. Either the block size exceeds "
            "the data's contiguous runs or the gap tolerance is splitting the series "
            "too finely."
        )
    if n_draws and short_segment_fraction > 0:
        # `.1%` rounds any share below 0.05% to "0.0%", so the warning read "cannot move
        # 0.0% of the treatment rows" - a sentence that tells the reader nothing is wrong
        # while warning that something is, which is how a warning gets trained out. Two
        # significant figures on the share keeps a genuinely tiny frozen fraction legible
        # as tiny rather than as zero.
        share = f"{short_segment_fraction:.2%} of"
        if short_segment_fraction < 0.0001:
            share = f"a {short_segment_fraction:.2e} fraction of"
        warn_the_reader(
            f"block permutation with block_size={block_size} cannot move "
            f"{share} the treatment rows: they sit in segments "
            "too short to hold two blocks, so the placebo distribution holds them at "
            "their observed values and the refutation p-value is biased toward 1. Read "
            "placebo_frozen_fraction alongside the p-value, and lower block_size or "
            "widen gap_tolerance if the frozen share is large.",
            source="block_permutation",
            key=("frozen_blocks", block_size, share),
            stacklevel=4,
        )


def empirical_permutation_p(placebo_effects: np.ndarray, observed_effect: float) -> float:
    """Two-sided Monte Carlo p-value for the block-permutation refutation.

    The observed statistic is itself one draw the permutation distribution can
    produce, so both the count and the denominator take the plus-one correction
    (Davison and Hinkley 1997; Phipson and Smyth 2010). Without it, a run in
    which no placebo reaches the observed effect reports ``p = 0.000`` - a claim
    no finite number of permutations can support. With ``n`` placebo draws the
    smallest p-value the test can report is ``1 / (n + 1)``.

    The statistic is the caller's choice and this function does not care which, but
    `run_dml_analysis` passes t-statistics rather than raw effects, and the difference is
    not cosmetic: a permuted treatment is no longer predictable from the controls, so
    ``var(T_res)`` - the second stage's whole denominator - inflates, and every placebo
    theta is shrunk toward zero by arithmetic. Comparing thetas therefore measures a
    narrower distribution than the null it stands for.

    Parameters
    ----------
    placebo_effects : np.ndarray
        The statistic from each successful placebo permutation.
    observed_effect : float
        The same statistic computed on the unpermuted data.

    Returns
    -------
    float
        The fraction of the permutation distribution at least as extreme as the
        observed value in absolute value, in ``(0, 1]``.
    """
    placebo = np.asarray(placebo_effects, dtype=float)
    # Both sides have to be finite, and neither check is defensive - a NaN on either side
    # silently biases the answer in the same direction, toward significance.
    #
    # `np.abs(x) >= abs(nan)` is False for every x, so a non-finite observed statistic
    # scores zero placebos as extreme and this returns 1/(n+1), the smallest p-value the
    # test can produce. Measured through `run_dml_analysis` with its own guard removed, 24
    # draws, only the observed fit's covariance failing: empirical_p = 0.04, and
    # `classify_refutation` published "Passes" against an undefined observed statistic.
    # That guard covers `run_dml_analysis`; `15_causal_estimation/03_econml_dml.py` and
    # `04_dml_crypto_regime.py` call this function directly and would not be covered by it.
    #
    # A NaN inside `placebo` fails its own comparison the same way and is counted as "not
    # extreme", which shrinks the numerator. Both notebooks already append only finite
    # draws, so this raises for no caller that exists; it is here because the filtering
    # belongs to whoever builds the array and the failure is silent if they forget.
    if not np.isfinite(observed_effect):
        raise ValueError(
            f"empirical_permutation_p needs a finite observed statistic, got "
            f"{observed_effect!r}. Every comparison against it is False, so the p-value "
            f"would be 1/(n+1) - the most significant value the test can report - with "
            f"nothing to say it was not measured."
        )
    if placebo.size and not np.isfinite(placebo).all():
        raise ValueError(
            f"empirical_permutation_p needs finite placebo draws, got "
            f"{int((~np.isfinite(placebo)).sum())} non-finite of {placebo.size}. They "
            f"count as not extreme and shrink the p-value. Drop the failed draws and "
            f"report how many there were."
        )
    at_least_as_extreme = int(np.sum(np.abs(placebo) >= abs(observed_effect)))
    return (1.0 + at_least_as_extreme) / (1.0 + placebo.size)


def classify_refutation(empirical_p: float, n_successful: int | None = None) -> str:
    """Pass, fail, or too few draws to tell, at the 5 % level.

    Returns "Passes" if the empirical placebo p-value is below 5 % (the observed effect
    cannot be reproduced by permutation in most placebo runs), and "Fails" otherwise.

    "Underpowered" is the third answer, and it is the honest one whenever the number of
    successful draws puts "Passes" out of reach. The plus-one correction floors the
    reported p-value at ``1 / (n + 1)``, so at 19 successful draws or fewer the smallest
    value the test can produce is already at or above 5 % and "Fails" would be published
    whatever the data said - untrue by construction in the same way ``p = 0.000`` was at
    the other end. The preview tier runs ten draws, so this is the ordinary case there,
    not an exotic one.

    ``n_successful`` is optional so that a caller holding only a p-value keeps the
    two-way answer; pass it wherever the draw count is known.
    """
    if n_successful is not None and 1.0 / (n_successful + 1) >= REFUTATION_ALPHA:
        return "Underpowered"
    return "Passes" if empirical_p < REFUTATION_ALPHA else "Fails"


def _resolve_panel_columns(
    df: pd.DataFrame,
    time_col: str | None,
    entity_col: str | None,
) -> tuple[str | None, str | None]:
    """Resolve canonical panel columns while preserving single-series inputs."""
    if time_col is None and "timestamp" in df.columns:
        time_col = "timestamp"
    if entity_col is None:
        entity_col = next((name for name in ("symbol", "product") if name in df.columns), None)
    if entity_col is not None and time_col is None:
        raise ValueError("time_col is required when an entity column is present")
    if time_col is not None and df[time_col].duplicated().any() and entity_col is None:
        raise ValueError("entity_col is required when decision times contain multiple rows")
    return time_col, entity_col


def run_dml_analysis(
    df: pd.DataFrame,
    treatment_col: str,
    outcome_col: str,
    confounder_cols: list[str],
    n_folds: int = 5,
    embargo: int = 21,
    n_placebo: int = 100,
    block_size: int = 21,
    seed: int = 42,
    hac_maxlags: int | None = None,
    horizon: int | None = None,
    time_col: str | None = None,
    entity_col: str | None = None,
    model_y=None,
    model_t=None,
    expected_step: str | pd.Timedelta | None = None,
    thread_limit: int = DML_THREAD_LIMIT,
) -> dict:
    """Full DML analysis pipeline: naive OLS, DML, and refutation tests.

    Parameters
    ----------
    df : pd.DataFrame
        Analysis dataset sorted by time.
    treatment_col : str
        Treatment variable column name.
    outcome_col : str
        Outcome variable column name.
    confounder_cols : list[str]
        Confounder column names.
    n_folds : int
        Number of walk-forward CV folds.
    embargo : int
        Gap periods between train and test.
    n_placebo : int
        Number of block permutation replications.
    block_size : int
        Block size for permutation test.
    seed : int
        Random seed.
    hac_maxlags : int or None
        HAC bandwidth passed through to the second stage. If None, resolved
        from `horizon`.
    horizon : int or None
        Label horizon in observation periods, forwarded to the second-stage
        HAC regression so the Newey-West bandwidth satisfies L >= horizon - 1.
        Pass it for overlapping labels (horizon >= ~10). It is the outcome
        horizon, read with
        `resolve_label_horizon(case_study_id, label, setup)`, and not the CV
        buffer, which bounds a different quantity and can be longer. When both
        `horizon` and `hac_maxlags` are None, the bandwidth falls back to the
        horizon-blind cube-root rule and a warning is emitted.
    time_col : str or None
        Ordered decision-time column. Inferred from canonical ``timestamp``
        when omitted. Required for panel data so cross-fitting, embargoes, and
        placebo blocks keep each decision time intact.
    entity_col : str or None
        Panel entity column. Inferred from canonical ``symbol`` or ``product``
        when omitted. Supply with ``time_col`` for non-canonical panels so
        placebo blocks are permuted within entity histories.

    Returns
    -------
    dict
        Comprehensive results with keys: naive_effect, naive_n_obs, dml_result,
        confounding_bias, confounding_bias_pct, refutation (z_score,
        empirical_p, placebo_mean, placebo_std, placebo_effects,
        refutation_class), p_value_hac, hac_maxlags, and n_obs.
    """
    # Also wrapped here, not only inside manual_dml_timeseries. The naive-OLS comparison
    # below runs np.linalg.lstsq after that call returns, and LAPACK's dgelsd reaches
    # threaded BLAS on a tall design - so naive_effect, and confounding_bias_pct which is a
    # difference between it and the pinned theta, would otherwise vary with the ambient
    # pool while the spec records deterministic_reduction: True. The inner context nests
    # harmlessly and still covers callers that use manual_dml_timeseries directly.
    with threadpool_limits(limits=thread_limit):
        # Input validation
        time_col, entity_col = _resolve_panel_columns(df, time_col, entity_col)
        n = len(df)
        min_rows = (n_folds + 1) * 50 + n_folds * embargo
        if n < min_rows:
            raise ValueError(
                f"Need at least {min_rows} rows for {n_folds}-fold CV with embargo={embargo}, got {n}"
            )
        if df[treatment_col].std() < 1e-10:
            raise ValueError(f"Treatment '{treatment_col}' has near-zero variance")
        if df[outcome_col].std() < 1e-10:
            raise ValueError(f"Outcome '{outcome_col}' has near-zero variance")

        if hac_maxlags is None and horizon is None:
            warn_the_reader(
                "no horizon or hac_maxlags given; the second-stage HAC bandwidth falls "
                "back to the horizon-blind cube-root rule, which under-lags overlapping "
                "labels of horizon >= ~10 and overstates the t-statistic. Pass the outcome "
                "horizon in observation periods: read the horizon with "
                "resolve_label_horizon(case_study_id, label, setup) - not the CV buffer, "
                "which can be longer - and convert it against the panel's own cadence. "
                "embargo_from_buffer without observed_step applies per-unit defaults "
                "instead, which read 24H as one period on an eight-hour panel.",
                source="run_dml_analysis",
                key=("hac_fallback", treatment_col, outcome_col),
            )

        _dml_started_at = datetime.now(UTC).isoformat()
        _dml_t0 = time.perf_counter()

        rng = np.random.default_rng(seed)

        T = df[treatment_col].values
        Y = df[outcome_col].values
        X = df[confounder_cols].values
        groups = df[time_col].values if time_col is not None else None
        units = df[entity_col].values if entity_col is not None else None

        # DML estimate
        dml = manual_dml_timeseries(
            Y,
            T,
            X,
            n_folds=n_folds,
            embargo=embargo,
            return_residuals=True,
            hac_maxlags=hac_maxlags,
            horizon=horizon,
            groups=groups,
            model_y=model_y,
            model_t=model_t,
            thread_limit=thread_limit,
        )

        # Compare the adjusted estimate with naive OLS on the exact second-stage
        # population. Earlier walk-forward dates have no out-of-fold residuals and
        # cannot enter only one side of the comparison.
        valid = np.isfinite(dml["Y_res"]) & np.isfinite(dml["T_res"])
        naive_n_obs = int(valid.sum())
        T_const = np.column_stack([np.ones(naive_n_obs), T[valid]])
        naive_coef = np.linalg.lstsq(T_const, Y[valid], rcond=None)[0]
        naive_effect = naive_coef[1]

        # Confounding bias
        dml_effect = dml["theta"]
        bias = naive_effect - dml_effect
        bias_pct = bias / abs(dml_effect) * 100 if dml_effect != 0 else 0.0

        # Block permutation refutation
        placebo_effects = []
        placebo_t_stats = []
        placebo_n_obs = []
        unchanged_draws = 0
        moved_fractions: list[float] = []
        # What no draw can move is a property of the segmentation and the block size, so
        # it is computed from them rather than inferred from what the draws happened to
        # do. The two reasons a row never moves are unrelated and only one of them is
        # worth a warning: see `_immobile_masks`.
        short_frozen, remainder_frozen = _immobile_masks(
            len(T),
            block_size,
            _permutation_segments(len(T), groups, units, expected_step, None),
        )
        for _ in range(n_placebo):
            T_perm = block_permute(
                T,
                block_size,
                rng=rng,
                groups=groups,
                units=units,
                expected_step=expected_step,
            )
            moved_fractions.append(_placebo_moved_fraction(T, T_perm))
            unchanged_draws += _placebo_is_unchanged(T, T_perm)
            perm_result = manual_dml_timeseries(
                Y,
                T_perm,
                X,
                n_folds=n_folds,
                embargo=embargo,
                hac_maxlags=hac_maxlags,
                horizon=horizon,
                groups=groups,
                model_y=model_y,
                model_t=model_t,
                thread_limit=thread_limit,
            )
            if not np.isnan(perm_result["theta"]) and np.isfinite(perm_result["t_stat_hac"]):
                if perm_result["n_obs"] != dml["n_obs"]:
                    raise RuntimeError(
                        "Observed and placebo DML statistics use different second-stage samples"
                    )
                placebo_effects.append(perm_result["theta"])
                placebo_t_stats.append(float(perm_result["t_stat_hac"]))
                placebo_n_obs.append(int(perm_result["n_obs"]))

        frozen_fraction = float(short_frozen.mean()) if short_frozen.size else 0.0
        remainder_fraction = float(remainder_frozen.mean()) if remainder_frozen.size else 0.0
        _assert_placebo_permutation_possible(
            unchanged_draws, n_placebo, block_size, frozen_fraction
        )

        refutation = {}
        observed_t = float(dml["t_stat_hac"])
        if not np.isfinite(observed_t):
            # The draws are fine; the statistic they would be compared against is not.
            # `empirical_permutation_p` counts placebos at least as extreme as the
            # observed one, and every `>=` against NaN is False, so the comparison would
            # return the smallest p-value the test can produce - 1/(n+1) - and publish
            # "Passes" on an undefined observed statistic. The permutation p-value is the
            # one number in this dict that is not a diagnostic, so the answer is to
            # withhold the verdict rather than to qualify it: `refutation` stays empty,
            # which is the shape a caller already handles for too few draws, and
            # `covariance_type` on the fit says which of the two happened.
            #
            # `run_resolved_causal_request` refuses the run before this matters. Direct
            # callers of `run_dml_analysis` - the chapter-15 notebooks - do not, and they
            # are the ones who would have read the verdict.
            warn_the_reader(
                f"the observed t-statistic is not finite "
                f"(covariance_type={dml.get('covariance_type')!r}), so the "
                f"{len(placebo_t_stats)} placebo draws have nothing to be compared "
                f"against. Reporting no refutation rather than a verdict computed "
                f"against NaN.",
                source="run_dml_analysis",
                key=("nonfinite_t", dml.get("covariance_type"), len(placebo_t_stats)),
                category=RuntimeWarning,
            )
        elif len(placebo_effects) >= MIN_PLACEBO_DRAWS:
            # THE TEST IS ON THE T-STATISTIC, NOT ON THETA, and the difference is not
            # cosmetic: comparing thetas made this refutation anti-conservative on every
            # run ever recorded.
            #
            # DML's second stage regresses the residualized outcome on the residualized
            # treatment, so var(T_res) is the estimator's whole denominator. Permuting the
            # treatment also frees it from the controls: the first stage can no longer
            # predict it, its residual keeps essentially all of its variance, and the
            # placebo estimator therefore divides by a much larger number than the observed
            # one does. A placebo theta comes out smaller than the observed theta by
            # arithmetic, whether or not there is any alignment to find - so the permutation
            # distribution is narrower than the null it is supposed to represent, and the
            # observed effect looks extreme against it more often than it should. The bias
            # runs one way, toward "Passes", so no stored value is safe to read as evidence.
            #
            # Measured on a synthetic panel with theta EXACTLY ZERO by construction, so
            # every rejection is a false positive that needs no interpretation:
            #
            #   var(T_res)/var(T)     observed 0.100   placebo mean 1.132   -> 11.4x
            #   placebo theta sd / observed se_hac                             0.167
            #   empirical p on raw effects                                    0.0164
            #   empirical p on t-statistics                                   0.5902
            #   placebo t distribution                        mean 0.017, sd 1.069
            #
            # The same ratio measured 10.9x on the crypto_perps_funding panel, which is a
            # different dataset entirely. The t-statistic is what cancels the denominator:
            # each draw divides by its own standard error, and what is left is the
            # alignment the refutation is about. Its null comes out centred on zero with
            # unit spread without anyone tuning for it, which is the calibration a
            # permutation test is supposed to have.
            #
            # `placebo_effects` stays in this dict, and no notebook plots it any more - the three
            # that draw the permutation distribution all read `placebo_t_stats`, because that is
            # the scale the verdict is decided on. It is kept because the effect scale
            # is the one a reader can interpret against the estimate, so a row carries both and
            # the registry schema says the same. `placebo_t_stats` is what `empirical_p`,
            # `z_score`, `placebo_mean` and `placebo_std` describe. `refutation_statistic` names
            # the scale in this dict, for a caller holding the fit; it is not registered and
            # `CausalResult.metrics` does not expose it, so what tells a reader of a stored row
            # which comparison produced its p-value is `refutation_placebo_t_json` being
            # non-NULL.
            placebo_arr = np.array(placebo_effects)
            placebo_t_arr = np.array(placebo_t_stats)
            p_mean = float(np.mean(placebo_t_arr))
            p_std = float(np.std(placebo_t_arr))
            z = (observed_t - p_mean) / p_std if p_std > 0 else np.inf
            emp_p = empirical_permutation_p(placebo_t_arr, observed_t)
            ref_class = classify_refutation(emp_p, len(placebo_t_stats))
            refutation = {
                "refutation_statistic": "t_stat_hac",
                "z_score": z,
                "empirical_p": emp_p,
                "placebo_mean": p_mean,
                "placebo_std": p_std,
                "observed_t_stat": observed_t,
                "placebo_effect_mean": float(np.mean(placebo_arr)),
                "placebo_effect_std": float(np.std(placebo_arr)),
                "n_successful": len(placebo_t_stats),
                "placebo_frozen_fraction": frozen_fraction,
                "placebo_remainder_fraction": remainder_fraction,
                "placebo_moved_fraction": float(np.mean(moved_fractions))
                if moved_fractions
                else 0.0,
                "n_folds": n_folds,
                "placebo_n_obs": placebo_n_obs,
                "placebo_effects": placebo_effects,
                "placebo_t_stats": placebo_t_stats,
                "refutation_class": ref_class,
            }

        return {
            "naive_effect": naive_effect,
            "naive_n_obs": naive_n_obs,
            "dml_result": dml,
            "confounding_bias": bias,
            "confounding_bias_pct": bias_pct,
            "refutation": refutation,
            "p_value_hac": dml.get("p_value_hac", np.nan),
            "hac_maxlags": dml.get("hac_maxlags", 0),
            "n_obs": len(df),
            "started_at": _dml_started_at,
            "elapsed_s": time.perf_counter() - _dml_t0,
        }


def format_dml_summary(results: dict) -> str:
    """Format DML analysis results for display."""
    dml = results["dml_result"]
    p_hac = results.get("p_value_hac", dml.get("p_value_hac", np.nan))
    hac_lags = results.get("hac_maxlags", dml.get("hac_maxlags", "?"))
    lines = [
        "=" * 60,
        "DML ANALYSIS SUMMARY",
        "=" * 60,
        f"Analysis rows: {results['n_obs']:,}",
        f"Second-stage rows: {dml.get('n_obs', results['n_obs']):,}",
        f"Second-stage decision times: {dml.get('n_periods', results['n_obs']):,}",
        f"Covariance: {dml.get('covariance_type', 'newey_west').replace('_', '-').title()}",
        (
            "HAC bandwidth: not applied; the robust covariance did not produce a "
            "usable estimate and SE (HAC) below is NaN"
            if dml