Panelgerechtes Double Machine Learning für ETF-Momentum-Effekte

# Double Machine Learning for Momentum Causal Effect


# Double Machine Learning for Momentum Causal Effect

**Chapter 15: Causal Estimation**
**Docker image**: `ml4t`
**Section Reference**: See Section 15.4 for DML theory and the ETF factor application

## Purpose
**Double machine learning** (DML) estimates the effect of a momentum signal on forward
returns while adjusting for volatility and regime, on a panel of ETFs. The estimate is one
of two things this notebook is about. The other is that a panel breaks every correction
that has a notion of "nearby": the standard error, the cross-validation folds, the
permutation blocks and the temporal subsets all count in rows unless they are told not to,
and on a stacked panel a row is a different ETF rather than the next day.

## Learning Objectives
- LO1: Explain what confounding does to a factor regression, and what adjustment can and
  cannot recover
- LO2: Fit a DML estimator with EconML for a continuous treatment
- LO3: Compare the raw and adjusted estimates and read the direction of the bias
- LO4: Build folds, standard errors and permutations that count in decision times
- LO5: Read a refutation test, including one that cannot refute anything

## Cross-References
- **Upstream**: Chapter 8 (ETF momentum features)
- **Downstream**: Chapter 16 (strategy simulation), Chapter 19 (risk management)
- **Related**: [`02_dowhy_causal_graph`](02_dowhy_causal_graph.ipynb) (graphical approach), [`04_dml_crypto_regime`](04_dml_crypto_regime.ipynb) (regime effects)

## Methodological Notes
Following Chernozhukov et al. (2017) and de Prado (2018):
- **WalkForwardCV** from ml4t-diagnostics, with purging and embargo, built over the
  panel's decision times
- **Driscoll-Kraay standard errors**, which aggregate by decision time before applying the
  Newey-West kernel
- **Block permutation within entity**, so the placebo treatment keeps its persistence

**Prerequisites**: [`02_dowhy_causal_graph`](02_dowhy_causal_graph.ipynb) for the DoWhy
workflow, and an ETF modeling dataset built by the features pipeline

## Causal Design Contract

| Element                   | This notebook                                                                                |
|---------------------------|----------------------------------------------------------------------------------------------|
| Unit                      | ETF-date row from the ETF modeling panel                                                     |
| Treatment                 | `skip_recent_6_1` - 6/1 momentum factor (continuous)                                          |
| Outcome                   | `fwd_ret_21d` - 21-day forward return                                                         |
| Controls (W in EconML)    | `vol_21d`, `vol_126d`, `regime`, `yield_curve_slope` - all backward-looking, pre-treatment    |
| Effect modifiers (X)      | None in this notebook; constant ATE target. See `04_dml_crypto_regime` for the X-slot example |
| Identification assumption | Selection on observables given the four controls; sufficient pre-treatment information       |
| Main failure mode         | Unobserved confounding (sentiment / macro shocks), nuisance-model misspecification, panel-time leakage if CV is not date-grouped |
| Estimand                  | Marginal effect of a one-unit change in 6/1 momentum on 21-day forward return after adjustment |

## Setup

```python
"""Double Machine Learning for Momentum Causal Effect - estimate causal effect of momentum on forward returns."""

import warnings

import numpy as np
import pandas as pd
from ml4t.diagnostic.splitters import WalkForwardCV
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import Ridge

from case_studies.utils.causal import (
    block_permute,
    empirical_permutation_p,
    manual_dml_timeseries,
)
from utils.modeling import load_modeling_dataset
from utils.reproducibility import set_global_seeds

# scikit-learn repeats a notice, once per nuisance fit, that a frame carrying feature names
# was fitted and a bare array predicted; EconML does that internally. Convergence and
# numerical warnings stay visible.
warnings.filterwarnings("ignore", category=UserWarning, module="sklearn.utils.validation")

# Statsmodels for HAC standard errors
import statsmodels.api as sm
from econml.dml import LinearDML
from statsmodels.regression.linear_model import OLS
```

## Configuration

```python
# Configuration - readers can modify these
CASE_STUDY_ID = "etfs"
PRIMARY_LABEL = "fwd_ret_21d"
MAX_SYMBOLS = 0
CV_FOLDS = 5
MAX_SAMPLES = 50000  # Temporal subsample if dataset too large
SEED = 42

# Cross-validation parameters for WalkForwardCV
FORWARD_HORIZON = 21  # 21-day forward returns
LABEL_HORIZON = FORWARD_HORIZON  # Purge overlapping samples
EMBARGO_PCT = 0.01  # 1% embargo after test set

# Refutation test parameters
N_PLACEBO_PERMUTATIONS = 100
BLOCK_SIZE = 21  # Block size for permutation (match forward horizon)
```

```python
set_global_seeds(SEED)
print(f"Seed: {SEED}")
```

## 1. Load ETF Features from Modeling Pipeline

We use `load_modeling_dataset()` to load pre-computed features (Ch8),
temporal features (Ch9), and labels, joined and ready for analysis.
Real-data only - no synthetic fallback. If the modeling dataset is missing,
the notebook fails loudly with a clear error rather than silently switching
to a synthetic substitute that would publish indistinguishable numbers.

```python
# A load failure is fatal rather than a fallback to synthetic data.
mds = load_modeling_dataset(CASE_STUDY_ID, PRIMARY_LABEL, max_symbols=MAX_SYMBOLS)

treatment_col = "skip_recent_6_1"
outcome_col = mds.label_col
confounder_cols = ["vol_21d", "vol_126d", "regime", "yield_curve_slope"]
entity_col = mds.entity_cols[0]

available = set(mds.dataset.columns)
missing = [c for c in [treatment_col, outcome_col] + confounder_cols if c not in available]
if missing:
    raise RuntimeError(
        f"Required columns missing from modeling dataset "
        f"{CASE_STUDY_ID}/{PRIMARY_LABEL}: {missing}. "
        f"Available features: {mds.feature_names[:20]}... "
        f"Set ML4T_DATA_PATH and rebuild the Ch8 features pipeline for case "
        f"study '{CASE_STUDY_ID}'."
    )

# Convert to pandas for sklearn/econml, sorted by date
analysis_cols = [mds.date_col] + mds.entity_cols + [treatment_col, outcome_col] + confounder_cols
df = (
    mds.dataset.select([c for c in analysis_cols if c in available])
    .drop_nulls()
    .sort(mds.date_col)
    .to_pandas()
)
# Temporal subsample if too large: take the most recent N unique dates so the
# subsample never cuts through a cross-section. `df.iloc[-MAX_SAMPLES:]` would
# slice at row level on a stacked panel and leave a fragmented final date.
if len(df) > MAX_SAMPLES:
    rows_per_date = df.groupby(mds.date_col).size().median()
    n_dates = int(np.ceil(MAX_SAMPLES / max(rows_per_date, 1)))
    keep_dates = df[mds.date_col].drop_duplicates().iloc[-n_dates:]
    df = df[df[mds.date_col].isin(keep_dates)].reset_index(drop=True)
    print(f"Taking most recent {n_dates} dates ({len(df):,} rows) from {len(keep_dates):,} dates")

print(f"Analysis data: {df.shape[0]:,} rows × {df.shape[1]} columns")
print(f"Date range: {df[mds.date_col].min()} to {df[mds.date_col].max()}")
print(f"Treatment: {treatment_col}, Outcome: {outcome_col}")
```

### The Panel Keys, and Why Everything Below Needs Them

A row here is an ETF *and* a date, and the frame is sorted by date, so consecutive rows
are usually different ETFs on the same day rather than the same ETF on consecutive days.
Every correction in this notebook that has a notion of "nearby" - the standard error, the
cross-validation folds, the permutation blocks, the temporal subsets - has to count in
decision times, not in rows. Counting in rows on this frame measures a slice of one day's
cross-section and calls it three weeks.

The two arrays below carry that information, and each of those four places takes them.

```python
decision_times = df[mds.date_col].to_numpy()
entities = df[entity_col].to_numpy()

unique_dates = np.sort(pd.unique(decision_times))
date_position = pd.Series(np.arange(len(unique_dates)), index=unique_dates)
row_date_position = date_position.reindex(decision_times).to_numpy()

print(
    f"{len(df):,} rows over {len(unique_dates):,} decision times "
    f"and {df[entity_col].nunique()} entities "
    f"({len(df) / len(unique_dates):.0f} rows per decision time)"
)


def panel_folds(n_splits):
    """Walk-forward folds built over decision times, then expanded to panel rows.

    WalkForwardCV counts `label_horizon` and the embargo in the positions it is handed. Fed
    the panel's rows it would purge a fraction of one date; fed the ordered unique dates it
    purges the 21 trading days the forward return actually spans. The row indices come back
    by membership, so no fold boundary cuts through a cross-section.
    """
    splitter = WalkForwardCV(
        n_splits=n_splits,
        label_horizon=LABEL_HORIZON,
        embargo_pct=EMBARGO_PCT,
        expanding=True,
    )
    folds = []
    for train_dates, test_dates in splitter.split(np.arange(len(unique_dates)).reshape(-1, 1)):
        folds.append(
            (
                np.flatnonzero(np.isin(row_date_position, train_dates)),
                np.flatnonzero(np.isin(row_date_position, test_dates)),
            )
        )
    return folds
```

## 2. The Unadjusted Slope

An OLS regression of forward returns on momentum, with nothing else in it. It is the
benchmark the adjusted estimates are read against, and it is confounded by construction:
volatility and the yield-curve regime both move the momentum signal and move forward
returns, so its slope carries their contribution as well as momentum's. Which way that
pushes the slope is an empirical question the comparison below answers.

### Which Robust Standard Error a Panel Takes

Overlapping 21-day returns make consecutive observations of one ETF correlated, and the
usual answer is a Newey-West standard error with the bandwidth set to the label horizon.
On this frame that answer is applied to the wrong axis. `cov_type="HAC"` runs its kernel
down the rows, and 21 rows here are a fifth of one day's cross-section, so the correction
treats different ETFs on the same day as if they were successive days.

**Driscoll-Kraay** is the version that fits a panel. It aggregates the regression score by
decision time first and applies the Newey-West kernel to that time series, which makes it
robust both to the serial correlation the overlap creates and to whatever the ETFs share
on a given day. statsmodels reaches it through `cov_type="hac-groupsum"` with a `time`
argument, which is the same call `case_studies/utils/causal.py` makes for every case study.

```python
y = df[outcome_col].values
X_with_const = sm.add_constant(df[[treatment_col]])

ols_iid = OLS(y, X_with_const).fit()
naive_estimate = float(ols_iid.params.iloc[1])
se_iid = float(ols_iid.bse.iloc[1])

# Driscoll-Kraay: the kernel runs over decision times, with the bandwidth at the label
# horizon because a 21-day forward return overlaps for 20 of every 21 days.
HAC_LAGS = FORWARD_HORIZON
time_codes = pd.factorize(decision_times, sort=False)[0]
ols_dk = ols_iid.get_robustcov_results(
    cov_type="hac-groupsum",
    time=time_codes,
    maxlags=HAC_LAGS,
    use_correction="hac",
    df_correction=False,
)
se_hac = float(np.sqrt(np.asarray(ols_dk.cov_params())[1, 1]))
t_stat_hac = naive_estimate / se_hac

naive_ci = (naive_estimate - 1.96 * se_hac, naive_estimate + 1.96 * se_hac)

print("=" * 60)
print("NAIVE ESTIMATE (ignoring confounders)")
print("=" * 60)
print(f"Coefficient: {naive_estimate:.6f}")
print(f"Standard Error (IID): {se_iid:.6f}")
print(f"Standard Error (Driscoll-Kraay): {se_hac:.6f}")
print(f"SE Inflation (DK/IID): {se_hac / se_iid:.2f}x")
print(f"95% CI (Driscoll-Kraay): [{naive_ci[0]:.6f}, {naive_ci[1]:.6f}]")
print(f"t-statistic (Driscoll-Kraay): {t_stat_hac:.2f}")
```

## 3. Double Machine Learning Estimate under Observed-Confounder Adjustment

DML uses a three-step orthogonalization process:

1. **Predict outcome from confounders**: $\hat{Y} = g(X)$ → residual $\tilde{Y} = Y - \hat{Y}$
2. **Predict treatment from confounders**: $\hat{T} = m(X)$ → residual $\tilde{T} = T - \hat{T}$
3. **Regress residual outcome on residual treatment**: $\tilde{Y} \sim \theta \tilde{T}$

The coefficient $\theta$ is an *orthogonalized* estimate of the treatment
effect under the maintained assumption that the specified controls capture
the relevant pre-treatment confounding variation. DML reduces sensitivity
to nuisance-model errors but does not, by itself, solve unobserved
confounding, simultaneity, interference, or bad-control bias.

The folds come from `WalkForwardCV` rather than `KFold` or `TimeSeriesSplit`, because
cross-fitting a forward-looking label needs two things neither of those provides. **Purging**
drops training rows whose 21-day label window overlaps the test window; the **embargo**
leaves a gap after the test window so autocorrelation does not carry back into the next
training set (Chernozhukov et al. 2017, de Prado 2018).

The confounders go into EconML's `W` slot, which holds controls used for residualization.
The `X` slot is for effect modifiers, the variables along which the treatment effect is
allowed to vary. Putting plain confounders in `X` still produces a single ATE, so nothing
visibly breaks, but it teaches the wrong habit for `04_dml_crypto_regime`, where the two
slots carry different variables and the distinction decides what the model estimates.

```python
Y = df[outcome_col].to_numpy()
T = df[treatment_col].to_numpy()
W = df[confounder_cols].to_numpy()  # controls used for residualization

# Walk-forward folds over decision times, expanded to rows (see `panel_folds` above).
cv = panel_folds(CV_FOLDS)
print(
    f"Using WalkForwardCV over decision times (label_horizon={LABEL_HORIZON} trading days, "
    f"embargo={EMBARGO_PCT:.1%})"
)
for i, (train_idx, test_idx) in enumerate(cv):
    print(f"  fold {i}: train {len(train_idx):,} rows, test {len(test_idx):,} rows")

dml = LinearDML(
    model_y=GradientBoostingRegressor(n_estimators=50, max_depth=3, random_state=SEED),
    model_t=GradientBoostingRegressor(n_estimators=50, max_depth=3, random_state=SEED),
    cv=cv,
    random_state=SEED,
)

dml.fit(Y, T, W=W)

dml_estimate = float(dml.ate())
dml_ci_lower, dml_ci_upper = (float(v) for v in dml.ate_interval(alpha=0.05))
bias = naive_estimate - dml_estimate
bias_pct = 100 * bias / abs(dml_estimate) if dml_estimate != 0 else float("inf")
```

```python
print("\n" + "=" * 60)
print("DOUBLE ML ESTIMATE (controlling for confounders)")
print("=" * 60)
print(f"Average Treatment Effect: {dml_estimate:.6f}")
if not np.isnan(dml_ci_lower):
    print(f"95% CI: [{dml_ci_lower:.6f}, {dml_ci_upper:.6f}]")
else:
    print(
        "95% CI: unavailable (EconML inference returned NaN - "
        "likely numerically degenerate first-stage residualization)"
    )

print("\n" + "=" * 60)
print("COMPARISON: Naive vs DML")
print("=" * 60)
print(f"Naive estimate: {naive_estimate:.6f}")
print(f"DML estimate:   {dml_estimate:.6f}")
print(f"Bias (Naive - DML): {bias:.6f}")
print(f"Bias percentage: {bias_pct:.1f}%")

if abs(naive_estimate) > abs(dml_estimate):
    print("\n-> The adjusted effect is smaller: the controls were inflating the raw slope")
else:
    print("\n-> The adjusted effect is larger: the controls were masking part of the slope")
```

**Interpretation**: The naive OLS estimate captures both the causal momentum effect and
spurious correlation induced by shared confounders (volatility, regime). DML's
orthogonalization removes this confounding by residualizing both treatment and outcome
against the confounders before estimating the final coefficient.

The direction of bias reveals the confounding structure. If the naive effect is
smaller in magnitude than the DML estimate, confounders *mask* the true effect -
for example, high volatility reduces both momentum and returns simultaneously.
If the naive effect is larger, confounders *inflate* the apparent predictive power.

**What it changes downstream**: a position size derived from the raw slope inherits
whatever the confounders contributed to it. The adjusted estimate is the one to size on
when volatility and regime are not separately hedged, and the gap between the two is how
much the raw slope was carrying.

## 4. The Same Estimate, Assembled by Hand

`manual_dml_timeseries` from `case_studies/utils/causal.py` runs the three steps in the
open: cross-fit the two nuisance models over walk-forward folds with an embargo, take the
residuals, regress one on the other. Handed `groups`, it builds those folds over decision
times and reports a Driscoll-Kraay standard error; without them it would do both by row.

```python
# Use shared manual DML implementation
dml_result = manual_dml_timeseries(
    df[outcome_col].values,
    df[treatment_col].values,
    df[confounder_cols].values,
    n_folds=CV_FOLDS,
    embargo=LABEL_HORIZON,
    groups=decision_times,
    horizon=FORWARD_HORIZON,
)

manual_ate = dml_result["theta"]
manual_se_iid = dml_result["se_iid"]
manual_se_hac = dml_result["se_hac"]
manual_t_hac = dml_result["t_stat_hac"]
manual_ci = (manual_ate - 1.96 * manual_se_hac, manual_ate + 1.96 * manual_se_hac)

print("\n" + "=" * 60)
print("MANUAL DML WITH WALK-FORWARD CV + EMBARGO")
print("=" * 60)
print(f"ATE estimate: {manual_ate:.6f}")
print(f"Standard Error (IID): {manual_se_iid:.6f}")
print(f"Standard Error (HAC): {manual_se_hac:.6f}")
print(f"SE Inflation (HAC/IID): {manual_se_hac / manual_se_iid:.2f}x")
print(f"95% CI (HAC): [{manual_ci[0]:.6f}, {manual_ci[1]:.6f}]")
print(f"t-statistic (Driscoll-Kraay): {manual_t_hac:.2f}")
```

## 5. Refutation Tests

**Critical for causal validity**: We validate the DML estimate using refutation tests.

1. **Temporal Placebo**: Regress Y on *lead* of T (should be ~0 if no reverse causality)
2. **Block Permutation**: Shuffle treatment in blocks to preserve autocorrelation
3. **Subset Stability**: Check if effect is stable across temporal subsets

```python
print("\n" + "=" * 60)
print("REFUTATION TESTS")
print("=" * 60)

# The lead is taken within symbol; a row-level shift would mix ETFs at the boundary.
print("\n1. TEMPORAL PLACEBO TEST (lead of treatment)")
df_placebo = df.sort_values([entity_col, mds.date_col]).copy()
df_placebo["treatment_lead"] = df_placebo.groupby(entity_col)[treatment_col].shift(-FORWARD_HORIZON)
# Back to date order before the fit. The shift needed entity-major rows; the folds need
# every row of a decision time adjacent, and `manual_dml_timeseries` rejects groups that
# are not sorted and contiguous rather than silently splitting one date across two folds.
df_placebo = (
    df_placebo.dropna(subset=["treatment_lead", outcome_col])
    .sort_values(mds.date_col, kind="stable")
    .reset_index(drop=True)
)

if len(df_placebo) > 100:
    placebo_result = manual_dml_timeseries(
        df_placebo[outcome_col].values,
        df_placebo["treatment_lead"].values,
        df_placebo[confounder_cols].values,
        n_folds=CV_FOLDS,
        embargo=LABEL_HORIZON,
        groups=df_placebo[mds.date_col].to_numpy(),
        horizon=FORWARD_HORIZON,
    )
    placebo_effect = placebo_result["theta"]

    print(f"   Lead treatment effect (DML): {placebo_effect:.6f}")
    print(f"   Original DML effect:         {manual_ate:.6f}")
    print(f"   Ratio |lead / original|:     {abs(placebo_effect) / abs(manual_ate):.3f}")
else:
    placebo_effect = None
    print("   Insufficient data for placebo test")
```

**Reading the temporal placebo.** A ratio near one is what this test returns here, and it
is not evidence of reverse causality. The treatment is a momentum signal that barely moves
over 21 trading days, so the horizon-shifted placebo is nearly the same variable as the
real treatment and reproduces a similar effect. For a treatment this persistent, the
shifted-signal placebo cannot refute anything.

The block permutation below is the refuter that can, because it breaks the alignment
between treatment and outcome while keeping the treatment's own persistence intact - but
only if the blocks are built along time.

### Blocks Along Time, Not Along Rows

`block_permute` given a bare array counts `BLOCK_SIZE` in the positions it is handed. On
this frame that is about a fifth of one day's cross-section, so a "block" is a handful of
ETFs on a single day. Permuting those destroys the serial dependence the test is supposed
to preserve, which is the iid shuffle `block_permute`'s own docstring warns makes placebo
tests too easy to pass.

Passing `groups` and `units` makes it permute within each ETF, where `BLOCK_SIZE` counts
that ETF's own ordered trading days and a block is the three weeks the name promises.

Each placebo also runs the estimator being tested, with the same fold count and embargo. A
placebo fitted on three folds against an estimate fitted on five is a different estimator
on a different number of cross-fitted rows, and the null would then be centred wherever
that difference puts it rather than where the absence of an effect does. The observed
statistic itself moves between the two settings.

```python
# Test 2: Block Permutation Test (uses shared block_permute)
print(f"\n2. BLOCK PERMUTATION TEST ({N_PLACEBO_PERMUTATIONS} permutations)")
placebo_effects = []
placebo_t_stats = []
permutation_failures = 0
T_original = df[treatment_col].values
rng = np.random.default_rng(SEED)

for i in range(N_PLACEBO_PERMUTATIONS):
    T_permuted = block_permute(
        T_original,
        BLOCK_SIZE,
        rng=rng,
        groups=decision_times,
        units=entities,
    )

    df_perm = df.copy()
    df_perm[treatment_col] = T_permuted

    try:
        perm_result = manual_dml_timeseries(
            df_perm[outcome_col].values,
            df_perm[treatment_col].values,
            df_perm[confounder_cols].values,
            n_folds=CV_FOLDS,  # the estimate's own setting; see the markdown above
            embargo=LABEL_HORIZON,
            groups=decision_times,
            horizon=FORWARD_HORIZON,
        )
        if np.isfinite(perm_result["t_stat_hac"]):
            placebo_t_stats.append(perm_result["t_stat_hac"])
            placebo_effects.append(perm_result["theta"])
        else:
            permutation_failures += 1
    except Exception as exc:
        permutation_failures += 1
        print(f"   Permutation {i} failed: {type(exc).__name__}: {exc}")

PERMUTATION_MIN_SUCCESS = max(10, int(0.5 * N_PLACEBO_PERMUTATIONS))
print(f"   Permutations: {len(placebo_effects)} successful, {permutation_failures} failed")
if len(placebo_effects) < PERMUTATION_MIN_SUCCESS:
    raise RuntimeError(
        f"Block-permutation placebo: only {len(placebo_effects)} successful "
        f"runs (need ≥{PERMUTATION_MIN_SUCCESS}); refuter cannot be trusted."
    )
```

### Compared on t-Statistics, Not on Effect Sizes

Permuting the treatment also frees it from the controls. The second stage regresses the
residualized outcome on the residualized treatment, so whatever the controls explain leaves
the denominator; a permuted treatment is no longer explained by them, its residual variance
is larger, and a placebo effect is mechanically smaller than the observed one whether or not
there is anything to find. Comparing raw effects against that distribution reports
significance the standard error does not support. Each permutation's t-statistic divides by
its own standard error, so the scale cancels and only the alignment between treatment and
outcome is left.

The p-value is the fraction of the placebo distribution at least as extreme as the observed
statistic, and it carries a plus-one correction because the observed statistic is itself one
draw that distribution could produce. Without the correction, a run in which no placebo
reaches it reports zero, which no finite number of permutations can establish. With n draws
the smallest reportable value is 1 / (n + 1), printed beside it as the floor. It is not a
false discovery rate, which is what this quantity used to be called.

**What cancelling the scale does not buy.** It removes one known bias and does not make the
test calibrated. A permutation test is valid when the permuted labels are exchangeable under
the null, and a treatment that persistent confounders predict is not: the shuffle breaks the
confounding along with the effect, so an estimate that is itself biased sits far from its own
permutation null and the test reports that distance. Measured on twelve synthetic panels with
a true effect of exactly zero and persistent AR(1) confounders, the studentized comparison
still rejects at the conventional five percent level on five of the twelve, against eleven for
the raw-effect comparison. `10_case_study_insights` carries that measurement. Read the count below as a comparison against a shuffle
rather than as a significance level.

```python
if len(placebo_t_stats) > 10:
    placebo_mean = np.mean(placebo_t_stats)
    placebo_std = np.std(placebo_t_stats)
    z_score = (manual_t_hac - placebo_mean) / placebo_std if placebo_std > 0 else np.inf
    # The plus-one correction is why the floor below is 1 / (n + 1); see the markdown above.
    permutation_p = empirical_permutation_p(np.asarray(placebo_t_stats), manual_t_hac)

    print(f"   Placebo t mean: {placebo_mean:.4f}")
    print(f"   Placebo t std:  {placebo_std:.4f}")
    print(f"   Observed t (Driscoll-Kraay): {manual_t_hac:.4f}")
    print(f"   Z-score vs placebo: {z_score:.2f}")
    print(
        f"   Permutation p-value: {permutation_p:.4f} (floor {1 / (len(placebo_t_stats) + 1):.4f})"
    )

    print(
        f"   Placebo draws at least as extreme: "
        f"{int(round(permutation_p * (len(placebo_t_stats) + 1))) - 1} of {len(placebo_t_stats)}"
    )
    print(
        f"   Placebo effect spread {np.std(placebo_effects):.6f} against a Driscoll-Kraay "
        f"standard error of {manual_se_hac:.6f}; the comparison above is on t-statistics, "
        f"which holds whatever the ratio of those two turns out to be"
    )
else:
    print("   Insufficient successful permutations")
    z_score = None
    permutation_p = None
```

**The z-score and the p-value can disagree, and the count is the one that holds.** The
z-score measures how far the observed t-statistic sits from the placebo *mean* in placebo
standard deviations, which is a statement about a normal distribution centred where the
placebos are. The permutation p-value counts how many placebo draws reach its magnitude.
When the placebo distribution is not centred near zero the two answer different questions,
and only the count is a statement about the null the test actually built. Read the count
printed above, and the mean and standard deviation beside it, before reading the z-score.

**This null is not centred at zero, and that is worth looking at rather than reporting.**
A permutation that implemented "no association" would put the placebo statistics around
zero. Two obvious explanations do not survive a check: demeaning the fitted residuals by
symbol leaves the placebo estimates positive, and so does demeaning them by date, so
neither a preserved between-symbol nor a preserved between-date component accounts for it.
What does is not settled here. Until it is, the count is the statistic to read - it asks
how often a permuted treatment reaches the observed statistic, which stays a fair question
whatever the distribution's centre - and the z-score, which measures distance from that
centre in placebo standard deviations, is a statement about a distribution the test has not
explained.

The two halves are cut at a decision time rather than at a row, so neither holds a
fragment of a cross-section. That is the same reason the subsample near the top keeps
whole dates.

```python
print("\n3. SUBSET STABILITY TEST (temporal halves)")
midpoint_date = unique_dates[len(unique_dates) // 2]
df_first_half = df[decision_times < midpoint_date]
df_second_half = df[decision_times >= midpoint_date]

if len(df_first_half) > 100 and len(df_second_half) > 100:
    result_first = manual_dml_timeseries(
        df_first_half[outcome_col].values,
        df_first_half[treatment_col].values,
        df_first_half[confounder_cols].values,
        n_folds=3,
        groups=df_first_half[mds.date_col].to_numpy(),
        horizon=FORWARD_HORIZON,
    )
    result_second = manual_dml_timeseries(
        df_second_half[outcome_col].values,
        df_second_half[treatment_col].values,
        df_second_half[confounder_cols].values,
        n_folds=3,
        groups=df_second_half[mds.date_col].to_numpy(),
        horizon=FORWARD_HORIZON,
    )

    effect_first = result_first["theta"]
    se_first = result_first["se_hac"]
    effect_second = result_second["theta"]
    se_second = result_second["se_hac"]

    print(f"   First half effect:  {effect_first:.6f} (SE: {se_first:.6f})")
    print(f"   Second half effect: {effect_second:.6f} (SE: {se_second:.6f})")

    diff = abs(effect_first - effect_second)
    diff_se = np.sqrt(se_first**2 + se_second**2)
    diff_z = diff / diff_se if diff_se > 0 else 0

    print(f"   Difference: {diff:.6f} (z = {diff_z:.2f})")

    print(f"   First half: {len(df_first_half):,} rows before {midpoint_date}")
    print(f"   Second half: {len(df_second_half):,} rows from {midpoint_date}")
else:
    print("   Insufficient data for subset test")
```

## 6. Nuisance Model Sensitivity

DML results depend on the quality of nuisance models. Let's compare different choices.

```python
results = []

nuisance_models = [
    ("Linear", Ridge(alpha=1.0), Ridge(alpha=1.0)),
    (
        "GBM (shallow)",
        GradientBoostingRegressor(n_estimators=50, max_depth=2, random_state=SEED),
        GradientBoostingRegressor(n_estimators=50, max_depth=2, random_state=SEED),
    ),
    (
        "GBM (deep)",
        GradientBoostingRegressor(n_estimators=100, max_depth=5, random_state=SEED),
        GradientBoostingRegressor(n_estimators=100, max_depth=5, random_state=SEED),
    ),
]

cv_sensitivity = panel_folds(3)

for name, model_y, model_t in nuisance_models:
    dml_test = LinearDML(model_y=model_y, model_t=model_t, cv=cv_sensitivity, random_state=SEED)
    dml_test.fit(Y, T, W=W)
    ate = float(dml_test.ate())
    results.append({"Nuisance Model": name, "ATE Estimate": ate})

sensitivity_df = pd.DataFrame(results)
print("\n" + "=" * 60)
print("NUISANCE MODEL SENSITIVITY")
print("=" * 60)
print(sensitivity_df.to_string(index=False))
# A ratio of magnitudes reports 1.00x, which reads as no spread at all, when the
# specifications straddle zero symmetrically - the one case where the nuisance choice
# changes the conclusion rather than its size. The signed range cannot do that.
lo = float(sensitivity_df["ATE Estimate"].min())
hi = float(sensitivity_df["ATE Estimate"].max())
print(
    f"\nAcross {len(sensitivity_df)} nuisance specifications the estimate spans "
    f"{lo:.6f} to {hi:.6f}, a range of {hi - lo:.6f} in {outcome_col} units"
)
print(f"  Specifications agreeing in sign: {'all' if lo * hi > 0 else 'not all'}")
```

## 7. Results Summary

```python
print("\n" + "=" * 60)
print("CHAPTER 15 NOTEBOOK RESULTS: 03_econml_dml.py")
print("=" * 60)

results_summary = {
    "naive_estimate": naive_estimate,
    "naive_se_iid": se_iid,
    "naive_se_hac": se_hac,
    "naive_ci_lower": naive_ci[0],
    "naive_ci_upper": naive_ci[1],
    "manual_dml_estimate": manual_ate,
    "manual_dml_se_iid": manual_se_iid,
    "manual_dml_se_hac": manual_se_hac,
}

if dml_estimate is not None:
    results_summary["econml_dml_estimate"] = dml_estimate
    results_summary["econml_dml_ci_lower"] = dml_ci_lower
    results_summary["econml_dml_ci_upper"] = dml_ci_upper
    results_summary["bias_from_confounding"] = bias
    results_summary["bias_pct"] = bias_pct

if z_score is not None:
    results_summary["placebo_z_score"] = z_score
if permutation_p is not None:
    results_summary["permutation_p"] = permutation_p

for key, value in results_summary.items():
    if isinstance(value, float):
        print(f"{key}: {value:.6f}")
    else:
        print(f"{key}: {value}")
```

## Key Methodological Points

### The four choices that decide whether the estimate means anything
1. **Cross-fitting with purging and embargo**, over decision times rather than rows.
   Purging drops training rows whose labels overlap the test window; the embargo adds a
   buffer after it (Chernozhukov et al. 2017, de Prado 2018).
2. **Driscoll-Kraay standard errors**, so the serial-correlation correction is applied
   along time and the cross-sectional dependence within a date is absorbed rather than
   counted as extra observations.
3. **Block permutation within entity**, so the placebo distribution is built from a
   treatment that keeps its persistence.
4. **A nuisance-model sweep**, because the point estimate depends on the first stage and
   the spread across specifications is the honest width of the finding.

### Three Intervals on One Effect

The notebook produces three statements of uncertainty about the same quantity, and they do
not agree. EconML's `ate_interval` treats the residualized observations as independent
draws. The manual DML path reports a Driscoll-Kraay standard error, which aggregates by
decision time and absorbs whatever the ETFs share on a day. The permutation test compares
the estimate against a null built by shuffling the treatment within each ETF.

The spread between them is not a defect in any one of them. It is the price of the panel:
an interval is only as good as its account of what is independent, and 52,000 ETF-days are
not 52,000 independent observations.

```python
se_inflation = se_hac / se_iid

print("Quantitative Findings")
print("-" * 40)
print(f"SE inflation, Driscoll-Kraay over iid: {se_inflation:.2f}x")
if dml_estimate is not None:
    direction = "larger" if abs(naive_estimate) > abs(dml_estimate) else "smaller"
    print(
        f"Adjustment moves the slope by {abs(bias_pct):.1f}% - the unadjusted estimate is "
        f"{direction} in magnitude"
    )
    if not np.isnan(dml_ci_lower):
        print(
            f"EconML DML: {dml_estimate:.6f}, iid 95% CI "
            f"[{dml_ci_lower:.6f}, {dml_ci_upper:.6f}], width {dml_ci_upper - dml_ci_lower:.6f}"
        )
    else:
        print(f"EconML DML: {dml_estimate:.6f} (interval unavailable)")
print(
    f"Manual DML:  {manual_ate:.6f}, Driscoll-Kraay 95% CI "
    f"[{manual_ci[0]:.6f}, {manual_ci[1]:.6f}], width {manual_ci[1] - manual_ci[0]:.6f}"
)
print(
    f"Naive OLS:   {naive_estimate:.6f}, Driscoll-Kraay 95% CI [{naive_ci[0]:.6f}, {naive_ci[1]:.6f}]"
)
if z_score is not None:
    print(f"Placebo z-score (on t-statistics): {z_score:.2f}")
if permutation_p is not None:
    print(f"Permutation p-value: {permutation_p:.4f}")
```

## Key Takeaways

1. **DML is an adjusted estimate, not a free lunch**. The orthogonalized
   coefficient is interpretable as a causal effect only under the maintained
   assumptions: pre-treatment controls are sufficient, positivity holds, no
   interference, and the specified DAG is correct.

2. **The robust standard error has to match the data's axis.** Overlapping 21-day labels
   make the iid standard error too small, and on a stacked panel the ordinary Newey-West
   correction fixes that along the wrong axis: its lags run across the cross-section.
   Driscoll-Kraay aggregates by decision time first, which is what the SE inflation
   printed above measures.

3. **Manual DML matches EconML conceptually and differs numerically.** The point estimate
   moves with nuisance-model flexibility, and the sweep in section 6 shows by how much.
   Neyman orthogonality is a guarantee about first-order sensitivity to nuisance error,
   not about agreement across nuisance choices.

4. **A refutation counts in the units its name claims.** Block permutation preserves
   autocorrelation only when the blocks run along an entity's own trading days; on a
   flattened panel the same call permutes a slice of one day's cross-section, which is
   the iid shuffle the test exists to avoid. The same is true of the walk-forward folds
   and of the temporal halves.

5. **The sign of the confounding bias is a result, not a rule.** Here the controls change
   the momentum slope in one direction; on the crypto funding premium in
   `02_dowhy_causal_graph` they change it in the other. The comparison printed above says
   which happened on this data, and it is the printed numbers, not this sentence, that
   settle it.