# diff-diff

> A Python library for Difference-in-Differences causal inference analysis. Provides sklearn-like estimators with statsmodels-style output for econometric analysis.

- Version: 3.9.0
- Repository: https://github.com/igerber/diff-diff
- License: MIT
- Dependencies: numpy, pandas, scipy (no statsmodels dependency)
- Optional: Rust backend for performance (via maturin)

## Quick Start

```python
import pandas as pd
from diff_diff import DifferenceInDifferences, generate_did_data

# Generate synthetic data with known treatment effect
data = generate_did_data(n_units=200, treatment_effect=5.0, seed=42)

# Fit basic 2x2 DiD
did = DifferenceInDifferences()
results = did.fit(data, outcome='outcome', treatment='treated', post='post')
print(results.summary())
print(f"ATT: {results.att:.3f} (SE: {results.se:.3f})")
```

## Design Patterns

- **sklearn-like API**: All estimators use `fit()` method, `get_params()`/`set_params()` for configuration.
- **Formula interface**: Supports R-style formulas like `"outcome ~ treated * post"`.
- **Results objects**: Rich dataclass containers with `summary()`, `to_dict()`, `to_dataframe()`.
- **Estimator aliases**: Short names available (e.g., `DiD`, `CS`, `SA`, `BJS`, `SDiD`, `TWFE`, `DDD`, `EDiD`, `SCM`, `Bacon`). The former `CDiD`/`Gardner`/`Stacked` aliases are deprecated (3.9, removed in 4.0) - use `ContinuousDiD`/`TwoStageDiD`/`StackedDiD`.

## Practitioner Workflow (based on Baker et al. 2025)

For rigorous DiD analysis, follow the 8-step framework (call `diff_diff.get_llm_guide("practitioner")`).
After estimation, call:

```python
from diff_diff import practitioner_next_steps
guidance = practitioner_next_steps(results)
```

Returns context-aware guidance on remaining diagnostic steps (parallel trends
testing, sensitivity analysis, heterogeneity checks, robustness comparisons).

## Estimators

### DifferenceInDifferences

Basic 2x2 Difference-in-Differences estimator.

```python
DifferenceInDifferences(
    vcov_type: str | None = None,           # Variance family: "hc1" (default), "classical", "hc2", "hc2_bm", "conley"
    cluster: str | None = None,             # Column for cluster-robust SEs
    alpha: float = 0.05,                    # Significance level
    inference: str = "analytical",          # "analytical" or "wild_bootstrap" (wild_bootstrap requires cluster=)
    n_bootstrap: int = 999,                 # Bootstrap replications; >= 2 under inference="wild_bootstrap"
    bootstrap_weights: str = "rademacher",  # "rademacher", "webb", or "mammen"
    seed: int | None = None,                # Random seed
    rank_deficient_action: str = "warn",    # "warn", "error", or "silent"
    df_convention: str = "residual",        # Analytical t/p/CI df: "residual" (n-K, default), "cluster" (Stata/fixest G-1, clustered fits), or "normal" (deliberate z); survey df + hc2_bm BM-DOF keep precedence; default flips at v4
)
```

**Alias:** `DiD`

**fit() parameters:**

```python
did.fit(
    data: pd.DataFrame,
    outcome: str = None,                   # Outcome variable column
    treatment: str = None,                 # Treatment indicator column (0/1)
    post: str = None,                      # Post-treatment indicator column (0/1)
    formula: str = None,                   # R-style formula (e.g., "y ~ treated * post")
    covariates: list[str] = None,          # Linear control variables
    fixed_effects: list[str] = None,       # Low-dimensional FE (dummy variables)
    absorb: list[str] = None,             # High-dimensional FE (within-transformation)
) -> DiDResults
```

`time=` remains accepted as a deprecated alias for `post=` (FutureWarning; removed in 4.0), and the `robust=` constructor flag as a deprecated alias for `vcov_type=`.

**Usage:**

```python
from diff_diff import DifferenceInDifferences

did = DifferenceInDifferences()
results = did.fit(data, outcome='y', treatment='treated', post='post')
results.print_summary()

# Formula interface
results = did.fit(data, formula='y ~ treated * post')

# With covariates and fixed effects
results = did.fit(data, outcome='y', treatment='treated', post='post',
                  covariates=['age', 'income'], absorb=['firm_id'])
```

### TwoWayFixedEffects

Two-Way Fixed Effects estimator for panel data. Inherits from DifferenceInDifferences.

```python
TwoWayFixedEffects(
    robust: bool = True,
    cluster: str | None = None,   # Auto-clusters at unit level if None
    alpha: float = 0.05,
    df_convention: str = "residual",  # Analytical t/p/CI df: "residual" (default), "cluster" (G-1), or "normal" (z); flips at v4
)
```

**Alias:** `TWFE`

**fit() parameters:**

```python
twfe.fit(
    data: pd.DataFrame,
    outcome: str,
    treatment: str,
    post: str,                    # 0/1 post dummy (renamed from time=, which warns through 3.9)
    unit: str,
    covariates: list[str] = None,
    event_study: bool = False,    # per-period event study returning EventStudyResults
    spec: str = "within",         # "within" (unit-FE event study) | "pooled" (the MultiPeriodDiD design)
    reference_period: Any = None, # ES mode: omitted period (default: last pre-period, e=-1)
    post_periods: list = None,    # ES mode: post-treatment period values
    time: str = ...,              # ES mode: the CALENDAR column (keyword); static: deprecated alias for post
) -> DiDResults | EventStudyResults
```

**Usage:**

```python
from diff_diff import TwoWayFixedEffects

twfe = TwoWayFixedEffects()
results = twfe.fit(data, outcome='y', treatment='treated', post='post', unit='unit_id')
results.print_summary()

# Event-study mode (absorbs the deprecated MultiPeriodDiD; spec="pooled"
# reproduces its design; wild bootstrap raises in this mode; the unit
# auto-cluster applies with the static carve-outs)
es = twfe.fit(data, outcome='y', treatment='treated', unit='unit_id',
              event_study=True, time='period', post_periods=[4, 5, 6, 7])
es.print_summary()
```

**Note:** TWFE can be biased with staggered treatment timing and heterogeneous effects. Consider CallawaySantAnna, SunAbraham, or ImputationDiD for staggered designs.

### MultiPeriodDiD

DEPRECATED (3.9, removed in 4.0; ledger row M-010): use
`TwoWayFixedEffects(...).fit(..., event_study=True)` - `spec="pooled"`
reproduces this design exactly. Constructing MultiPeriodDiD (or its
EventStudy alias) emits a FutureWarning. Event-study style DiD with
period-specific treatment effects. Inherits from DifferenceInDifferences.

```python
MultiPeriodDiD(
    robust: bool = True,
    cluster: str | None = None,
    alpha: float = 0.05,
    df_convention: str = "residual",  # Analytical t/p/CI df: "residual" (default), "cluster" (G-1), or "normal" (z); flips at v4
)
```

**Alias:** `EventStudy`

**fit() parameters:**

```python
mp_did.fit(
    data: pd.DataFrame,
    outcome: str,
    treatment: str,
    time: str,
    post_periods: list = None,             # Which periods are post-treatment
    covariates: list[str] = None,
    fixed_effects: list[str] = None,
    absorb: list[str] = None,
    reference_period: Any = None,          # Reference period (default: last pre-period)
) -> MultiPeriodDiDResults
```

**Usage:**

```python
from diff_diff import MultiPeriodDiD, plot_event_study

did = MultiPeriodDiD()
results = did.fit(data, outcome='sales', treatment='treated',
                  time='period', post_periods=[4, 5, 6, 7])
results.print_summary()
plot_event_study(results)
```

### CallawaySantAnna

Callaway-Sant'Anna (2021) estimator for staggered DiD with heterogeneous treatment effects.

```python
CallawaySantAnna(
    control_group: str = "never_treated",        # "never_treated" or "not_yet_treated"
    anticipation: int = 0,                       # Anticipation periods
    estimation_method: str = "dr",               # "dr", "ipw", or "reg"
    alpha: float = 0.05,
    cluster: str | None = None,                  # Cluster col; activates CR1 on the IF via synthesized SurveyDesign(psu=col). None → per-unit IF.
    n_bootstrap: int = 0,                        # 0 = analytical SEs, 999+ recommended
    bootstrap_weights: str | None = None,        # "rademacher", "mammen", or "webb"
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    base_period: str = "varying",                # "varying" or "universal"
    cband: bool = True,                          # Simultaneous confidence bands
    pscore_trim: float = 0.01,                   # Propensity score trimming bound
    vcov_type: str = "hc1",                      # {"hc1"} only — IF-based variance per Callaway & Sant'Anna (2021). Analytical-sandwich {classical, hc2, hc2_bm} and conley REJECTED at __init__ (see REGISTRY.md IF-vs-sandwich subsection).
)
```

**Alias:** `CS`

**fit() parameters:**

```python
cs.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,              # Column: first treatment period (0 or inf for never-treated)
    covariates: list[str] = None,
    aggregate: str = None,         # None, "simple", "event_study", "group", or "all"
    balance_e: int = None,         # Balance event study at this relative period
) -> CallawaySantAnnaResults
```

**Usage:**

```python
from diff_diff import CallawaySantAnna, plot_event_study

# NOTE: no n_bootstrap here - post-fit aggregate() is ANALYTICAL-ONLY and
# recompute levels raise NotImplementedError on a bootstrapped fit
# ('simple' relays the stored bootstrap inference - see below).
cs = CallawaySantAnna(estimation_method="dr")
results = cs.fit(data, outcome='outcome', unit='unit', time='period',
                 first_treat='first_treat')
results.print_summary()

# Aggregate post-fit - no refit. Returns a NEW object; `results` is unchanged.
event_study = results.aggregate('event_study')   # -> EventStudyResults
by_cohort   = results.aggregate('group')         # -> AggregationResult
```

Fit-time `aggregate=` / `balance_e=` are DEPRECATED since 3.9 (removed in 4.0,
ledger rows M-020 / M-117) and emit a `FutureWarning`. The downstream consumers
all accept the post-fit container directly - `plot_event_study`,
`compute_honest_did` and `compute_pretrends_power` each take
`results.aggregate('event_study')` - so the only case still requiring the
fit-time path is **bootstrap inference**: CallawaySantAnna's `aggregate()`
fails closed on a bootstrapped fit's RECOMPUTE levels ('event_study'/'group') rather than substituting analytical inference for
percentile-bootstrap statistics.

```python
# Post-fit route (recommended): aggregate once, feed any consumer.
es = results.aggregate('event_study')
plot_event_study(es)
honest = compute_honest_did(es, M=1.0)

# Bootstrap fits only: aggregate at fit time.
boot = CallawaySantAnna(estimation_method="dr", n_bootstrap=999, seed=42)
plotted = boot.fit(data, outcome='outcome', unit='unit', time='period',
                   first_treat='first_treat', aggregate='event_study')
plot_event_study(plotted)
```

### ChaisemartinDHaultfoeuille

de Chaisemartin & D'Haultfœuille (2020/2022) estimator for **non-absorbing (reversible) treatments**. The most general library estimator for treatments that switch on AND off over time (allows dynamic/carryover effects + joiner/leaver decomposition); `LPDiD` (`non_absorbing="first_entry"`/`"effect_stabilization"`) and `TROP` (`non_absorbing=True`, no-dynamic-effects) also handle non-absorbing treatment under stronger assumptions. Ships `DID_M` (= `DID_1` at horizon `l = 1`) plus the full multi-horizon event study `DID_l` for `l = 1..L_max` from the dynamic companion paper (NBER WP 29873). Includes normalized estimator `DID^n_l`, cost-benefit aggregate `delta`, dynamic placebos `DID^{pl}_l`, and sup-t simultaneous confidence bands.

```python
ChaisemartinDHaultfoeuille(
    alpha: float = 0.05,
    cluster: str | None = None,                  # Must be None; non-None raises NotImplementedError
    n_bootstrap: int = 0,                        # 0 = analytical SE only
    bootstrap_weights: str = "rademacher",       # "rademacher", "mammen", or "webb"
    seed: int | None = None,
    placebo: bool = True,                        # Auto-compute single-lag placebo
    twfe_diagnostic: bool = True,                # Auto-compute Theorem 1 TWFE decomposition
    drop_larger_lower: bool = True,              # Drop multi-switch groups (matches R DIDmultiplegtDYN)
    by_path: int | None = None,                  # Top-k per-path event study; requires drop_larger_lower=False, L_max>=1; supports binary or integer-coded discrete D (D in Z); composes with survey_design (analytical TSL + replicate-weight; multiplier bootstrap n_bootstrap>0 still gated under survey) and heterogeneity (per-path predict_het, mirrors R did_multiplegt_dyn(..., by_path, predict_het); composes with placebo=True for per-path placebo predict_het on backward horizons — survey_design + placebo + heterogeneity warns and falls back to forward-horizon-only heterogeneity until the pre-period cell allocator is derived); mutex with paths_of_interest
    paths_of_interest: list[tuple[int, ...]] | None = None,  # User-specified path subset, alternative to by_path=k (Python-only API; mutex with by_path; composes with survey_design, heterogeneity, and placebo same as by_path=k)
    rank_deficient_action: str = "warn",         # Used by TWFE diagnostic OLS
)
```

**Alias:** `DCDH`

**fit() parameters:**

```python
est.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,                                   # Unit (group) identifier
    time: str,
    treatment: str,                              # Per-observation binary treatment or non-binary intensity
    # ---- multi-horizon ----
    L_max: int | None = None,                    # Max horizon; None = l=1 only
    # ---- covariates and extensions ----
    aggregate: Any = NOT_SUPPLIED,               # DEPRECATED (M-026): warns; non-None raises - use results.aggregate() post-fit
    covariates: list[str] | None = None,         # DID^X residualization-style covariates
    trends_linear: bool | None = None,           # DID^{fd} group-specific linear trends
    trends_nonparam: Any | None = None,          # DID^s state-set-specific trends
    honest_did: bool = False,                    # HonestDiD sensitivity on placebos
    # ---- survey support ----
    survey_design: SurveyDesign | None = None,   # pweight (TSL or replicate BRR/Fay/JK1/JKn/SDR); n_bootstrap > 0 uses Hall-Mammen PSU multiplier
) -> ChaisemartinDHaultfoeuilleResults
```

`L_max` controls multi-horizon computation. `covariates`, `trends_linear`, `trends_nonparam`, `honest_did`, `heterogeneity`, `design2`, and `survey_design` are all supported. Fit-time `aggregate=` is DEPRECATED (row M-026: it never computed anything here; supplying it warns and any non-None value raises) - aggregation is post-fit: `results.aggregate('event_study')` returns the unified `EventStudyResults` view (Phase-1 fits the 2-row l=1 view, `L_max>=1` the multi-horizon surface) and `results.aggregate('simple')` a one-row `AggregationResult` relaying the overall estimand (DID_M / DID_1 / delta) bit-exactly; both are pure views, so bootstrap fits are permitted. The dCDH container is deliberately NOT accepted by `compute_honest_did`/`compute_pretrends_power` (their l1 placebo semantics need the native dCDH branch - pass the results object itself). The 3.8 names `group=` and `controls=` remain accepted as deprecated aliases (FutureWarning; removed in 4.0).

**Usage:**

```python
from diff_diff import ChaisemartinDHaultfoeuille
from diff_diff.prep import generate_reversible_did_data

data = generate_reversible_did_data(
    n_groups=80, n_periods=6, pattern="single_switch", seed=42,
)

est = ChaisemartinDHaultfoeuille()
results = est.fit(
    data, outcome="outcome", unit="group",
    time="period", treatment="treatment",
)
results.print_summary()

# Decomposition
print(f"DID_M (overall):  {results.overall_att:.3f}")
print(f"DID_+ (joiners):  {results.joiners_att:.3f}")
print(f"DID_- (leavers):  {results.leavers_att:.3f}")
print(f"Placebo (DID^pl): {results.placebo_effect:.3f}")

# Multi-horizon event study
results = est.fit(data, outcome="outcome", unit="group",
                  time="period", treatment="treatment", L_max=3)
for h in sorted(results.event_study_effects):
    e = results.event_study_effects[h]
    print(f"  DID_{h} = {e['effect']:.3f} (SE={e['se']:.3f})")
print(f"Cost-benefit delta: {results.cost_benefit_delta['delta']:.3f}")
df = results.to_dataframe("event_study")  # includes placebos as negative horizons
```

**Standalone TWFE diagnostic** (without fitting the full estimator):

```python
from diff_diff import twowayfeweights

diagnostic = twowayfeweights(
    data, outcome="outcome", unit="group", time="period", treatment="treatment",
)
print(f"Plain TWFE coefficient: {diagnostic.beta_fe:.3f}")
print(f"Fraction of negative weights: {diagnostic.fraction_negative:.3f}")
print(f"sigma_fe (sign-flipping threshold): {diagnostic.sigma_fe:.3f}")
```

**Notes:**
- Validated against R `DIDmultiplegtDYN` v2.3.3 at horizon `l = 1` via `tests/test_chaisemartin_dhaultfoeuille_parity.py`
- Placebo SE contract: single-period placebo `DID_M^pl` (`L_max=None`) has `NaN` SE because the per-period aggregation path has no influence-function derivation; inference fields stay NaN-consistent **even when `n_bootstrap > 0`** for the single-period path (single-period bootstrap covers only `DID_M`, `DID_+`, and `DID_-`). Multi-horizon dynamic placebos `DID^{pl}_l` (`L_max >= 1`) have valid analytical SE via the placebo influence function (same cohort-recentered structure as positive horizons, applied to backward outcome differences), with bootstrap SE override when `n_bootstrap > 0` — bootstrap at `L_max >= 1` covers `DID_M`, `DID_+`, `DID_-`, per-horizon event-study effects (`event_study_effects[l]`), and placebo horizons (`placebo_event_study[-l]`), with shared weights across horizons for valid joint (sup-t) bands. This is a library extension beyond the dynamic companion paper, which states Theorem 1 variance for `DID_l` only.
- The analytical CI is conservative under Assumption 8 (independent groups) of the dynamic companion paper, exact only under iid sampling
- Survey design supported: pweight with strata/PSU/FPC via Taylor Series Linearization (analytical) or replicate-weight variance (BRR/Fay/JK1/JKn/SDR). Opt-in PSU-level Hall-Mammen wild bootstrap via `n_bootstrap > 0`. Replicate + `n_bootstrap > 0` rejected with `NotImplementedError` (replicate variance is closed-form)

### SunAbraham

Sun-Abraham (2021) interaction-weighted estimator for staggered DiD.

```python
SunAbraham(
    control_group: str = "never_treated",        # "never_treated" or "not_yet_treated"
    anticipation: int = 0,
    alpha: float = 0.05,
    cluster: str | None = None,                  # Defaults to unit-level clustering (dropped on explicit vcov_type='hc2' / 'classical')
    n_bootstrap: int = 0,                        # 0 = analytical cluster-robust SEs
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    vcov_type: str = "hc1",                       # {"classical","hc1","hc2","hc2_bm","conley"}; classical/hc2/hc2_bm via full-dummy saturated design; conley (Conley-1999 spatial-HAC) via within-transform — also pass conley_coords=(lat_col,lon_col), conley_cutoff_km, conley_lag_cutoff (unit auto-cluster dropped; explicit cluster= -> spatial+cluster product kernel). survey_design=/weights/n_bootstrap rejected for non-hc1 (use hc1 default for surveys)
    conley_coords: tuple[str, str] | None = None, # (lat_col, lon_col); required for vcov_type="conley"
    conley_cutoff_km: float | None = None,        # km bandwidth; required for vcov_type="conley"
    conley_metric: str = "haversine",             # "haversine" | "euclidean"
    conley_kernel: str = "bartlett",              # "bartlett" | "uniform"
    conley_lag_cutoff: int | None = None,         # within-unit Bartlett max lag (0 = spatial-only)
    df_convention: str = "residual",              # Analytical t/p/CI df for cells AND aggregates (3.9: aggregates share the cells' df — previously z): "residual" (default), "cluster" (G-1), "normal" (z); survey/BM DOF keep precedence; flips at v4
)
```

**Alias:** `SA`

**fit() parameters:**

```python
sa.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    covariates: list[str] = None,
) -> SunAbrahamResults
```

**Usage:**

```python
from diff_diff import SunAbraham

sa = SunAbraham()
results = sa.fit(data, outcome='outcome', unit='unit',
                 time='period', first_treat='first_treat')
results.print_summary()
```

### ImputationDiD

Borusyak-Jaravel-Spiess (2024) imputation DiD estimator. Efficient estimator producing shorter CIs than CS/SA under homogeneous effects.

```python
ImputationDiD(
    anticipation: int = 0,
    alpha: float = 0.05,
    cluster: str | None = None,                  # Defaults to unit-level clustering
    vcov_type: str = "hc1",                      # {"hc1"} only — IF-based variance per Borusyak et al. (2024) Theorem 3
    n_bootstrap: int = 0,                        # 0 = analytical (Theorem 3 variance)
    bootstrap_weights: str = "rademacher",       # "rademacher", "mammen", or "webb"
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    horizon_max: int | None = None,              # Max event-study horizon
    aux_partition: str = "cohort_horizon",        # "cohort_horizon", "cohort", or "horizon"
    pretrends: bool = False,                     # Include pre-treatment horizons in event study
    leave_one_out: bool = False,                 # BJS 2024 App. A.9 leave-one-out finite-sample variance (larger, less-downward-biased SE)
    df_convention: str = "residual",             # Pretrends lead-regression per-lead t/p/CI only (3.9: leads moved z -> t(residual)): "residual" (default), "cluster" (G-1), "normal" (z); BJS overall/post inference + Wald F knob-independent; flips at v4
)
```

**Alias:** `BJS`

**fit() parameters:**

```python
imp.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    covariates: list[str] = None,
    aggregate: Any = NOT_SUPPLIED,  # DEPRECATED (M-021): warns, still works - use results.aggregate() post-fit
    balance_e: Any = NOT_SUPPLIED,  # DEPRECATED (M-118): moves onto results.aggregate(balance_e=)
    survey_design: SurveyDesign = None,  # Optional design-based inference (pweight + analytical strata/PSU/FPC or replicate BRR/Fay/JK1/JKn/SDR)
) -> ImputationDiDResults
```

**Usage:**

```python
from diff_diff import ImputationDiD, plot_event_study

est = ImputationDiD()
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat')
results.print_summary()
es = results.aggregate('event_study')   # post-fit (M-021); balance_e= lives here too
plot_event_study(es)
```

### TwoStageDiD

Gardner (2022) two-stage DiD estimator. Point estimates match ImputationDiD; uses GMM sandwich variance.

```python
TwoStageDiD(
    anticipation: int = 0,
    alpha: float = 0.05,
    cluster: str | None = None,
    n_bootstrap: int = 0,
    bootstrap_weights: str = "rademacher",
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    horizon_max: int | None = None,
    vcov_type: str = "hc1",          # {"hc1"} only — Gardner (2022) two-stage GMM cluster-sandwich; analytical-sandwich {classical, hc2, hc2_bm} and conley REJECTED at __init__/fit (see REGISTRY.md IF-vs-sandwich subsection)
)
```

**Alias:** `Gardner` (deprecated 3.9, removed 4.0 - emits FutureWarning; use `TwoStageDiD`)

**fit() parameters:**

```python
ts.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    covariates: list[str] = None,
    aggregate: Any = NOT_SUPPLIED,  # DEPRECATED (M-022): warns, still works - use results.aggregate() post-fit
    balance_e: Any = NOT_SUPPLIED,  # DEPRECATED (M-119): moves onto results.aggregate(balance_e=)
) -> TwoStageDiDResults
```

**Usage:**

```python
from diff_diff import TwoStageDiD

est = TwoStageDiD()
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat')
results.print_summary()
es = results.aggregate('event_study')   # post-fit (M-022); balance_e= lives here too
```

### SpilloverDiD

Butts (2021) ring-indicator spillover-aware DiD. Augments two-stage Gardner with ring-indicator covariates that identify direct effect on treated (`tau_total`) and per-ring spillover effects on near-control units (`delta_j`). Handles non-staggered and staggered timing in one estimator. Recommends `vcov_type="conley"` with cutoff = `d_bar` (paper Section 3.1).

```python
SpilloverDiD(
    rings: list[float],                  # K+1 sorted breakpoints; K rings
    d_bar: float | None = None,          # Far-away cutoff (defaults to max(rings))
    vcov_type: str = "hc1",              # "hc1", "conley", or default cluster
    conley_coords: tuple[str, str] | None = None,  # (lat_col, lon_col), required
    conley_metric: str = "haversine",    # or "euclidean" / callable
    conley_cutoff_km: float | None = None,
    conley_lag_cutoff: int | None = None,
    cluster: str | None = None,
    alpha: float = 0.05,
    anticipation: int = 0,
    event_study: bool = False,           # Wave C: per-event-time × ring decomposition (Butts Table 2)
    horizon_max: int | None = None,      # Bin event-times outside [-H,+H] into endpoint pools (event-study mode); H>=1 or None — H=0 rejected (use event_study=False for aggregate spec)
    rank_deficient_action: str = "warn",
)
```

**fit() parameters:**

```python
sp.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    treatment: str | None = None,        # binary D_it; auto-converted to first_treat
    first_treat: str | None = None,      # OR onset time per unit (Gardner)
    covariates: list[str] | None = None, # Deferred: NotImplementedError if non-None
    survey_design: object = None,        # Supported (Wave E): analytical/TSL HC1 / CR1 / Conley; replicate weights not yet
) -> SpilloverDiDResults
```

**Restrictions and Wave C/D status:**

- `covariates=` raises `NotImplementedError` (planned follow-up). Gardner two-stage requires covariate effects estimated on the untreated-and-unexposed Omega_0 subsample at stage 1; appending raw covariates only at stage 2 silently biases `tau_total` / `delta_j` on panels with time-varying covariates.
- `survey_design=` is **supported** (Wave E): design-based variance via Binder (1983) Taylor-series linearization for `vcov_type="hc1"` (HC1) and `cluster=<col>` (CR1), and for `vcov_type="conley"` via a panel-aware stratified-Conley sandwich on per-period PSU totals (`conley_lag_cutoff > 0` adds a within-PSU serial Bartlett HAC and requires an effective PSU — supply `survey_design.psu` or inject one via `cluster=<col>`). Accepts `pweight` with strata / PSU / FPC. Replicate-weight variance (BRR / Fay / JK1 / JKn / SDR) is NOT yet supported (`NotImplementedError`): per Gerber (2026) Appendix A the IF-reweighting shortcut does not apply because `gamma_hat` is weight-sensitive, so correct support requires a per-replicate full re-fit.
- `vcov_type="classical"` raises `NotImplementedError` (Wave D restriction). Wave D GMM first-stage correction has not been derived for the homoskedastic meat structure `sigma_hat^2 * (X_10' X_10)`. Use `vcov_type="hc1"`, `vcov_type="conley"`, or pair with `cluster=<col>` for CR1 — all three apply the Wave D GMM correction.
- `event_study=True` SHIPPED (Wave C): emits per-event-time `tau_k` and per-(ring, event-time) `delta_jk` as `att_dynamic: pd.DataFrame` (indexed by event-time `k`) plus MultiIndex `spillover_effects: pd.DataFrame` (levels `(ring_label, event_time)`). TwoStageDiD-compatible `event_study_effects: Dict[int, Dict]` alias also emitted for `plot_event_study` consumption — `_extract_plot_data` prefers the new `reference_period` attribute over the legacy `n_obs==0` heuristic. (DiagnosticReport integration: WIRED — `SpilloverDiDResults` routes to parallel_trends (event-study on these direct-effect dynamics), design_effect, and heterogeneity; bacon excluded since spillover identifies off far-away controls, not TWFE 2x2 comparisons.) (schema: `{k: {"effect", "se", "n_obs", "t_stat", "p_value", "conf_int": (low, high)}}` mirroring `two_stage.py:1355-1389`). Reference period `ref_period = -1 - anticipation` (TwoStageDiD `two_stage.py:486` convention); reference row uses `coef=0.0, se=0.0, n_obs=0, conf_int=(0.0, 0.0)`. Scalar `att` field becomes a sample-share-weighted average of post-treatment `tau_k` (`att = sum_{k>=0} w_k * tau_k` with `w_k = n_treated_at_k / total`) with SE from linear-combination inference `Var(att) = w' V_subset w` on the post-treatment vcov block — no separate fit. **Two-clock K_it:** direct-effect clock is `K_direct = t - effective_first_treat(i)` for ever-treated rows; spillover clock is `K_spill = t - earliest-in-range-cohort-onset(i)` (running min across activated cohorts, NaN pre-trigger). `K_spill >= 0` structurally; negative-k spillover cells are rectangularly emitted with `coef = NaN, n_obs = 0`. **`horizon_max` semantics:** bins event-times outside `[-H, +H]` into endpoint pools (no observations dropped — divergence from TwoStageDiD which filters; intentional, per `feedback_no_silent_failures`). With `horizon_max=None`, auto-detects bin set from observed K. **Validation:** `horizon_max < 0` raises `ValueError`; `ref_period < -horizon_max` (i.e., `anticipation > horizon_max - 1`) raises `ValueError` — silently floor-shifting the reference would change identification. **Reduce-to-aggregate:** under constant-tau DGP with `horizon_max=None`, the share-weighted scalar `att` reproduces Wave B's aggregate bit-identically. **Note:** `horizon_max=0` does NOT reduce to Wave B (binning collapses pre-treatment K values to `k=0`, making `D^0 = D_i` ever-treated indicator rather than `D_it`). Per-event-time SEs include the Wave D Gardner GMM first-stage correction (see next bullet).
- Stage-2 variance applies the Gardner GMM first-stage uncertainty correction across HC1 / Conley / cluster (Wave D, SHIPPED). The IF outer-product formula `psi_i = gamma_hat' X_{10,i} eps_{10,i} - X_{2,i} eps_{2,i}` is used unconditionally; kernel `K` is path-dependent (identity for HC1, block-indicator for cluster, spatial kernel for Conley). Documented synthesis of Butts (2021) §3.1 + Gardner (2022) §4 + Conley (1999); no reference software combines all three. Point estimates unchanged from Wave B/C; SE values shift upward by 1-few percent.
- Only nearest-treated rings supported; `ring_method="count"` (count of treated neighbors in ring) not yet exposed

**Usage:**

```python
from diff_diff import SpilloverDiD

est = SpilloverDiD(
    rings=[0, 50, 100, 200],
    conley_coords=("lat", "lon"),
    vcov_type="conley",
    conley_cutoff_km=200.0,
    conley_lag_cutoff=0,
)
results = est.fit(data, outcome="y", unit="unit", time="time", treatment="D")
print(f"tau_total = {results.att:.4f}")
print(results.spillover_effects)         # per-ring DataFrame
```

### SyntheticDiD

Synthetic Difference-in-Differences (Arkhangelsky et al. 2021). Combines DiD with synthetic control by re-weighting control units.

```python
SyntheticDiD(
    zeta_omega: float | None = None,        # Unit weight regularization (auto-computed if None)
    zeta_lambda: float | None = None,       # Time weight regularization (auto-computed if None)
    alpha: float = 0.05,
    variance_method: str = "placebo",       # "placebo", "bootstrap", or "jackknife"
    n_bootstrap: int = 200,                 # Replications for variance estimation
    seed: int | None = None,
)
```

**Alias:** `SDiD`

**fit() parameters:**

```python
sdid.fit(
    data: pd.DataFrame,
    outcome: str,
    treatment: str,
    unit: str,
    time: str,
    post_periods: list,
) -> SyntheticDiDResults
```

**Usage:**

```python
from diff_diff import SyntheticDiD

sdid = SyntheticDiD(seed=42)
results = sdid.fit(data, outcome='outcome', treatment='treated',
                   unit='unit', time='period', post_periods=[5, 6, 7, 8])
results.print_summary()
weights_df = results.get_unit_weights_df()
```

### SyntheticControl

Classic Synthetic Control Method (Abadie, Diamond & Hainmueller 2010; Abadie & Gardeazabal 2003). Builds a single treated unit's counterfactual as a convex combination of never-treated "donor" units. **Donor weights only** (no time weights, no ridge) — distinct from `SyntheticDiD`.

**Alias:** `SCM`

```python
SyntheticControl(
    v_method: str = "nested",               # "nested" (pre-MSPE) | "cv" (out-of-sample, ADH 2015; needs predictors spanning BOTH train/val windows — default single-period lags are REJECTED) | "inverse_variance" (1/Var(X) on RAW predictors, bypasses standardize) | "custom"
    custom_v=None,                          # diagonal V (len = #predictors); required iff v_method="custom"; forbidden otherwise
    optimizer_options: dict | None = None,  # merged into scipy.optimize.minimize (outer V search; nested/cv)
    n_starts: int = 4,                      # multistart count for the outer V search (nested/cv)
    inner_max_iter: int = 10000,            # Frank-Wolfe inner-solve cap
    inner_min_decrease: float = 1e-5,       # inner-solve convergence scale (scale-aware)
    standardize: str = "std",               # "std" (per-row SD, ddof=1) or "none" (deviation from R)
    alpha: float = 0.05,
    seed: int | None = None,                # seeds the multistart Dirichlet draws
    v_cv_t0: int | None = None,             # cv train/val split index; default len(pre)//2 (Abadie 2021 t0=T0/2); None unless v_method="cv"
)
```

**fit() parameters:**

```python
scm.fit(
    data: pd.DataFrame,
    outcome: str,
    treatment: str,                         # ABSORBING 0/1 indicator (treated unit, post periods)
    unit: str,
    time: str,
    *,
    post_periods: list | None = None,       # inferred from D==1 periods if None
    treated_unit=None,                      # inferred as the single ever-treated unit if None
    predictors: list[str] | None = None,    # columns averaged over predictor_window
    predictors_op: str = "mean",            # "mean" | "sum" (linear combinations only, ADH 2010 §2.3)
    predictor_window: list | None = None,   # pre-periods to average over (default: all pre)
    special_predictors=None,                # list of (var, periods, op) triples
    pre_period_outcomes=None,               # "all" or list; per-period outcome lags
    donor_pool: list | None = None,         # explicit never-treated donors (default: all)
) -> SyntheticControlResults
```

**Inference:** NONE analytical — `se`/`t_stat`/`p_value`/`conf_int` are always NaN. `att` is the mean post-period gap. Significance via in-space placebo permutation inference: `results.in_space_placebo()` reassigns treatment to each donor, refits against the other J-1 donors (the real treated unit is excluded from every placebo pool), and sets `placebo_p_value = rank/(n_placebos+1)` from the post/pre RMSPE-ratio. The permutation `placebo_p_value` is a SEPARATE field from the (NaN) `p_value`; `is_significant` stays bound to `p_value`. **ADH-2015 §4 robustness (opt-in, analytical inference unchanged):** `results.leave_one_out()` drops each reportably-weighted donor (weight > 1e-6) and re-fits (per-drop ATT/`delta_att` — large `delta_att` ⇒ single-donor dependence); `results.in_time_placebo()` backdates the intervention and checks for a spurious pre-period gap (TRUNCATE windowing — predictor windows in the held-out region are dropped); `results.regression_weights()` computes the implied regression-counterfactual weights `W^reg = X0a'(X0a X0a')^{-1}X1a` (intercept-augmented, sums to 1 at full row rank) and flags donors outside `[0,1]` — the extrapolation an OLS counterfactual incurs but the simplex-constrained SC cannot (pure linear algebra, no refit; min-norm + `_regw_rank_deficient` when `k+1>J`); `results.sparse_synthetic_control(sizes=None, max_subsets=50000)` exhaustively searches `C(J,l)` size-`l` donor subsets holding `V` FIXED at the baseline (default sweep `[1,2,3]` skips over-cap sizes with a warning; an explicitly-requested oversize `l` raises), reporting how fit/ATT degrade as the synthetic is forced sparse (+ `get_sparse_synthetic_control_gaps()`). **Confidence sets by test inversion (Firpo-Possebom 2018 §4, opt-in, also non-analytical):** `results.test_sharp_null(effect, gamma=0.1)` tests `H_0: α_1t = f(t)` by re-ranking the in-space placebo gaps (no refits; `test_sharp_null(0)` == `placebo_p_value`), and `results.confidence_set(family="constant"|"linear", gamma=0.1)` inverts it into a confidence set for the effect path (constant-effect interval / linear-slope set, strict `p>gamma`), surfaced on `effect_confidence_set` / `get_confidence_set_df()` — the analytical `conf_int` still stays NaN. **Conformal inference (Chernozhukov-Wüthrich-Zhu 2021, opt-in, also non-analytical):** unlike the Firpo path (which re-ranks the cross-unit placebo gaps), the conformal layer fits its OWN time-permutation-invariant constrained-LS proxy (eqs 3–4, no V-matrix) under the null on ALL periods and permutes residuals over time for the single treated unit. `results.conformal_test(effect, q=1, scheme="moving_block")` gives a joint sharp-null p-value (eqs 1–2; `q∈{1,2,∞}`); `results.conformal_confidence_intervals(alpha=0.1)` gives pointwise per-period CIs (Algorithm 1 — each period `t` uses `Z=(pre, t)`, the other post-periods dropped); `results.conformal_average_effect(alpha=0.1)` gives an average-effect CI by collapsing into `T*`-blocks (Appendix A.1). All three accept `alternative={"two-sided","greater","less"}` (one-sided uses the SIGNED average-effect statistic per CWZ Remark 1, q fixed at 1; CIs become half-lines `[lower,+inf)`/`(-inf,upper]`) and `covariates=[...]` (pivoted-variable matching rows stacked RAW into the constrained-LS proxy per the note after eq 6; residuals/p-value stay outcome-only). `scheme="moving_block"` (default, serial-dependence-robust) or `"iid"` (finer). Surfaced on `conformal_inference` / `get_conformal_grid_df()`; `conf_int` stays NaN. Predictor periods must lie within the pre window; `post_periods` must be a contiguous suffix cross-checked against `D` (no anticipation).

**Usage:**

```python
from diff_diff import SyntheticControl

scm = SyntheticControl(v_method='nested', seed=0)
# Set predictor_window explicitly when a covariate is observed on only a subset of
# the pre periods — the default averages over ALL pre periods and fails closed on
# non-finite cells.
results = scm.fit(data, outcome='gdpcap', treatment='treated', unit='region', time='year',
                  predictors=['invest', 'school.high'],
                  predictor_window=[1964, 1965, 1966, 1967, 1968, 1969],
                  special_predictors=[('gdpcap', [1960, 1965, 1969], 'mean')])
results.print_summary()
gap_df = results.get_gap_df()       # period, gap, phase
weights_df = results.get_weights_df()  # unit, weight (descending)
```

### TripleDifference

Triple Difference (DDD) estimator following Ortiz-Villavicencio & Sant'Anna (2025).

```python
TripleDifference(
    estimation_method: str = "dr",            # "dr", "reg", or "ipw"
    cluster: str | None = None,
    vcov_type: str = "hc1",                   # {"hc1"} only — IF-based variance per Ortiz-Villavicencio & Sant'Anna (2025). Analytical-sandwich {classical, hc2, hc2_bm} and conley REJECTED at __init__ (see REGISTRY.md IF-vs-sandwich subsection).
    alpha: float = 0.05,
    pscore_trim: float = 0.01,
    rank_deficient_action: str = "warn",
)
```

**Alias:** `DDD`

**fit() parameters:**

```python
ddd.fit(
    data: pd.DataFrame,
    outcome: str,
    group: str,                # Treated group indicator (0/1)
    partition: str,            # Eligible partition indicator (0/1)
    post: str,                 # Post-treatment indicator (0/1)
    covariates: list[str] = None,
) -> TripleDifferenceResults
```

**Usage:**

```python
from diff_diff import TripleDifference

ddd = TripleDifference(estimation_method="dr")
results = ddd.fit(data, outcome='outcome', group='group',
                  partition='partition', post='post')
results.print_summary()
```

### ContinuousDiD

Continuous Difference-in-Differences estimator (Callaway, Goodman-Bacon & Sant'Anna 2024). Estimates dose-response curves ATT(d) and ACRT(d).

```python
ContinuousDiD(
    degree: int = 3,                          # B-spline degree (3 = cubic)
    num_knots: int = 0,                       # Interior knots
    dvals: np.ndarray | None = None,          # Custom dose evaluation grid
    control_group: str = "never_treated",     # "never_treated", "not_yet_treated", or "lowest_dose" (Remark 3.1, P(D=0)=0)
    anticipation: int = 0,
    base_period: str = "varying",             # "varying" or "universal"
    alpha: float = 0.05,
    n_bootstrap: int = 0,
    bootstrap_weights: str = "rademacher",
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    covariates: list[str] | None = None,      # DEPRECATED here (M-084; removed in 4.0) — pass covariates= to fit() instead
    estimation_method: str = "dr",            # "reg" or "dr" (used only with covariates); "ipw" not
                                              #   supported on the dose curve (raises NotImplementedError)
    pscore_trim: float = 0.01,                # dr propensity trimming
    epv_threshold: float = 10.0,              # dr propensity events-per-variable
    pscore_fallback: str = "error",           # dr propensity-failure action: "error" (fail-closed
                                              #   default) or "unconditional" (opt-in reg-like fallback)
    treatment_type: str = "continuous",       # "continuous" (B-spline curve) or "discrete"
                                              #   (saturated per-dose-level regression, CGBS Eq 4.1)
)
```

Covariates give conditional parallel trends: each (g,t) cell's control counterfactual becomes a
covariate-adjusted prediction. `reg`/`dr` share the ATT(d) *shape* and ACRT(d); `dr` (doubly-robust,
default) differs only in the overall_att/ATT(d) level and SE. `overall_att`+SE match DRDID
reg_did_panel/drdid_panel. `covariates=` + `survey_design=` is not yet supported.

`treatment_type="discrete"` fits a **saturated regression** for a multi-valued dose: one indicator per
distinct dose level, so `ATT(d_j) = mean_{D=d_j}(ΔY) − control` (a per-level 2×2 DiD) and `ACRT(d_j)`
is the paper's backward finite difference on `{0, d_1, …, d_J}` (`ACRT(d_1) = ATT(d_1)/d_1`, so binary
`D ∈ {0,1}` gives `ACRT = ATT`). It reuses the full inference stack (analytical / bootstrap /
covariate / survey) and reduces to the per-level 2×2 DiD SE. Multi-cohort fits need a shared dose
support across cohorts (else `NotImplementedError`); an off-support `dvals` value raises `ValueError`.

**Alias:** `CDiD` (deprecated 3.9, removed 4.0 - emits FutureWarning; use `ContinuousDiD`)

**fit() parameters:**

```python
cdid.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    dose: str,                     # Column with continuous treatment dose
    aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-025) - removed in 4.0; use post-fit results.aggregate();
                                   #   legacy values None/"dose"/"eventstudy" still route until 4.0 (any
                                   #   supplied value warns; "dose" was always a no-op - the curves are
                                   #   always computed)
    covariates: list[str] = None,  # Conditional parallel trends (X). None = unconditional
) -> ContinuousDiDResults
```

**Usage:**

```python
from diff_diff import ContinuousDiD

est = ContinuousDiD(n_bootstrap=199, seed=42)
results = est.fit(data, outcome='outcome', unit='unit', time='period',
                  first_treat='first_treat', dose='dose')
results.print_summary()
dose_table = results.aggregate('dose')      # ATT(d) + ACRT(d) rows (view; any fit)
overall = results.aggregate('simple')       # att + acrt rows (view; any fit)
# Binarized event study: analytical fits only (bootstrapped fits raise -
# re-fit with n_bootstrap=0, or the deprecated fit-time aggregate='eventstudy'
# until 4.0):
es = ContinuousDiD(seed=42).fit(data, outcome='outcome', unit='unit',
                                time='period', first_treat='first_treat',
                                dose='dose').aggregate('event_study')
```

### HeterogeneousAdoptionDiD

HeterogeneousAdoption DiD estimator (de Chaisemartin, Ciccia, D'Haultfœuille & Knau 2026). Targets a Weighted Average Slope (WAS) at the dose support boundary on **Heterogeneous Adoption Designs** — designs where treatment varies in dose intensity across units. Comparison comes from dose variation across units. The estimator does NOT require dropping never-treated units: a small share of never-treated units is fully compatible (paper edge case — Garrett et al. 2020 retained 12 untreated counties out of 2,954), and on staggered event-study panels never-treated units are explicitly retained as the untreated-group comparison (paper Appendix B.2). Uses a bias-corrected local-linear estimator at the dose support boundary on continuous-dose designs (Design 1' / Design 1) and a 2SLS Wald-IV estimator on the mass-point design.

```python
HeterogeneousAdoptionDiD(
    design: str = "auto",              # "auto" / "continuous_at_zero" / "continuous_near_d_lower" / "mass_point"
    d_lower: float | None = None,      # Support infimum; auto-detected when None
    kernel: str = "epanechnikov",      # Local-linear kernel
    alpha: float = 0.05,
    vcov_type: str | None = None,      # Mass-point only: "classical" (default) or "hc1"
    robust: bool | None = None,        # DEPRECATED alias for vcov_type= (M-047; FutureWarning, removed in 4.0)
    cluster: str | None = None,        # Cluster column for cluster-robust SEs (all designs; event-study adds a cluster-robust sup-t band). cluster= + survey_design= rejected
    n_bootstrap: int = 999,            # Multiplier-bootstrap iterations for sup-t bands (event-study, when weighted/survey OR clustered)
    seed: int | None = None,
)
```

**Alias:** `HAD`

**fit() parameters:**

```python
had.fit(
    data: pd.DataFrame,
    outcome: str,
    dose: str,
    time: str,
    unit: str,
    first_treat: str | None = None,            # Required on staggered panels (last-cohort auto-filter trigger)
    aggregate: Any = NOT_SUPPLIED,             # DEPRECATED (M-027) - removed in 4.0; fit() selects the mode from the panel shape (2 distinct periods -> overall scalar WAS; more -> per-horizon event-study WAS) and post-fit results.aggregate() is the aggregation surface; supplying a value warns then runs the legacy routing
    *,
    cband: bool = True,                        # Simultaneous (sup-t) confidence bands on event-study fits that are survey-weighted OR clustered (keyword-only)
    survey_design: SurveyDesign | None = None, # Survey-design kwarg (weights, strata, PSU, FPC) — the sole weighting entry
    trends_lin: bool = False,                  # Eq 17 linear-trend detrending. Requires a multi-period panel (the event-study mode); needs F>=3 (pre-period depth) for the regression; rejects survey_design= (raises NotImplementedError under trends_lin).
    covariates: Any | None = None,             # NOT IMPLEMENTED — non-None raises NotImplementedError (deferred Appendix B.1 / Theorem 6 covariate-adjusted extension; pre-residualize the outcome on covariates as a workaround)
    outcome_col: str = ..., dose_col: str = ..., time_col: str = ..., unit_col: str = ..., first_treat_col: str = ...,  # deprecated aliases for the five column params (FutureWarning; removed in 4.0)
) -> HeterogeneousAdoptionDiDResults | HeterogeneousAdoptionDiDEventStudyResults
```

**Usage:**

```python
from diff_diff import HeterogeneousAdoptionDiD, did_had_pretest_workflow

# Both surfaces select their mode from the panel shape (M-027/M-139): a
# two-period panel runs the overall (single-period WAS) estimator / pretest
# battery, a multi-period panel runs the event-study ones. Use distinct data
# objects for the two regimes.

# Vet the testable identifying assumptions on the two-period panel first:
report = did_had_pretest_workflow(
    data_2p, outcome='y', unit='unit', time='t',
    dose='d', first_treat='first_treat')
print(report.summary())

# Single-period scalar WAS (the overall mode) on the two-period panel:
est = HeterogeneousAdoptionDiD()
results = est.fit(data_2p, outcome='y', unit='unit',
                  time='t', dose='d',
                  first_treat='first_treat')
print(results.summary())

# Multi-period per-horizon WAS (the event-study mode) on the multi-period panel:
es = est.fit(data_mp, outcome='y', unit='unit',
             time='t', dose='d',
             first_treat='first_treat')

# Post-fit aggregation (rows M-027/M-122) - pure views, no kit:
print(results.aggregate('simple').to_dataframe())   # 1-row WAS relay (target = the estimand label)
print(es.aggregate('event_study').to_dataframe())   # unified EventStudyResults container
```

**Staggered panels.** On multi-cohort multi-period panels (the event-study mode), `fit()` auto-filters to the last treatment cohort plus never-treated units (paper Appendix B.2) and emits a `UserWarning` naming kept/dropped counts. The estimand is then a **last-cohort-only WAS**, not a multi-cohort average. For full multi-cohort staggered support, see `ChaisemartinDHaultfoeuille`.

**Mass-point + survey constraint.** When fitting `design="mass_point"` with `survey_design=SurveyDesign(...)`, `vcov_type="hc1"` is required: the survey path composes the standard error via Binder-TSL on the HC1-scale influence function, so the default classical sandwich path raises `NotImplementedError` — on both the static path and the event-study path (the event-study rejection fires regardless of `cband`, since the Binder-TSL analytical SE consumes the HC1-scaled IF either way). The one exception is `cluster=`: a clustered mass-point fit resolves to the CR1 sandwich regardless of `vcov_type` (and uses the clustered sup-t band), so there is no classical-vs-HC1 mismatch. Passing `vcov_type="hc1"` is a safe default on weighted survey + sup-t examples since `vcov_type` is unused on the continuous designs (CCT-2014 robust SE is the only formula there).

### RegressionDiscontinuity

Regression discontinuity estimator - sharp and fuzzy, with optional covariate adjustment (Calonico, Cattaneo & Titiunik 2014; covariates per Calonico, Cattaneo, Farrell & Titiunik 2019), parity-targeting R rdrobust 4.0.0. SHARP (default): treatment is assigned by a known threshold of an observed running variable (`running >= cutoff`; units exactly at the cutoff are treated); no treatment column. FUZZY: pass the OBSERVED take-up column via `fit(..., takeup=...)` (R's `fuzzy=`) - the estimand becomes the local Wald ratio (complier LATE at the cutoff for BINARY take-up under monotonicity; ratio-of-jumps otherwise - the `estimand` field says which) with a linearized bias correction, and the results gain a full `first_stage*` three-row block. Point estimation via kernel-weighted local polynomials on each side; data-driven MSE/CER-optimal bandwidths (all 10 rdrobust selectors; fuzzy selects on the ratio objective by default, with a sharp-on-Y switch under one-sided perfect compliance or `sharpbw=True`); robust bias-corrected inference. COVARIATE ADJUSTMENT: pass `fit(..., covariates=[...])` (R's `covs=`) for the CCFT 2019 additive common-coefficient adjustment - the estimand is UNCHANGED (precision only, unlike the DiD estimators' conditional-parallel-trends role); requires covariate BALANCE at the cutoff (testable: fit each covariate as the outcome and inspect its RD p-value); bandwidths are covariate-aware. Cross-sectional - no panel/time dimension.

```python
RegressionDiscontinuity(
    cutoff: float = 0.0,               # Known threshold c of the running variable
    p: int = 1,                        # Local-polynomial order, 0..20 (1 = local linear, 0 = local constant)
    q: int | None = None,              # Bias-regression order; explicit q needs p < q <= 20; None -> p + 1 unvalidated (R-exact: p=20 yields q=21)
    kernel: str = "triangular",        # triangular / epanechnikov / uniform ("tri"/"epa"/"uni" accepted)
    bwselect: str = "mserd",           # 10-option rdrobust menu: mserd/msetwo/msesum/msecomb1/msecomb2 + cer* variants
    h: float | None = None,            # Manual main bandwidth (both sides); h alone -> b = h
    b: float | None = None,            # Manual bias bandwidth (requires h; ignored with a warning otherwise)
    rho: float | None = None,          # h/b ratio; with h -> b = h/rho; without h applies to SELECTED bandwidths
    vcov_type: str = "nn",             # Only "nn" in this release (same-side nearest-neighbor variance)
    nnmatch: int = 3,                  # Minimum NN matches for the variance
    masspoints: str = "adjust",        # "adjust" (rdrobust default) / "check" / "off"
    bwcheck: int | None = None,        # Force >= this many unique support points into the window
    bwrestrict: bool = True,           # Clamp bandwidths to the observed running-variable range
    scaleregul: float = 1.0,           # IK-style regularization scale (0 removes)
    sharpbw: bool = False,             # Fuzzy only: select bandwidths on the sharp reduced form (R's sharpbw); auto under one-sided perfect compliance
    covs_drop: bool = True,            # Covariate fits only: drop collinear covariates with a warning naming them (R default); False = strict error
    alpha: float = 0.05,               # rdrobust level = 100*(1-alpha)
)
```

**Alias:** `RDD`

**fit() parameters:**

```python
rd.fit(
    data: pd.DataFrame,
    outcome: str,
    running: str,
    takeup: str | None = None,         # None = sharp; a column name = fuzzy (observed take-up; any numeric, typically binary)
    covariates: list[str] | None = None,  # Pre-determined covariate columns (R's covs=); additive common-coefficient adjustment, SAME estimand
) -> RegressionDiscontinuityResults
```

The 3.8 names `outcome_col=`, `running_col=`, and `treatment_col=` remain accepted as deprecated aliases (FutureWarning; removed in 4.0).

```python
from diff_diff import RegressionDiscontinuity

rd = RegressionDiscontinuity(cutoff=0.0)
results = rd.fit(df, outcome="y", running="score")
results.att        # ROBUST bias-corrected estimate (canonical binding: t_stat == att/se, conf_int centered on att)
results.att_conventional  # rdrobust's printed headline coefficient (conventional local-polynomial estimate)
results.conf_int   # robust bias-corrected CI
results.h_left, results.b_left  # selected bandwidths
print(results.summary())  # three-row Conventional / Bias-Corrected / Robust table, as in rdrobust

fuzzy = rd.fit(df, "y", "score", takeup="takeup")  # fuzzy RD
fuzzy.att          # linearized bias-corrected local Wald ratio, robust row (complier LATE for binary take-up)
fuzzy.first_stage  # take-up jump (bias-corrected; full three-row first_stage* mirror available)
fuzzy.estimand     # "fuzzy (LATE for compliers at the cutoff)" (binary take-up) / "fuzzy (local Wald ratio at the cutoff; non-binary take-up)" / "sharp (ATE at the cutoff)"

adj = rd.fit(df, "y", "score", covariates=["age", "income"])  # covariate-adjusted (estimand unchanged, shorter CIs)
adj.covariate_coefficients  # {"age": ..., "income": ...} - nuisance projection gammas, NOT causal effects
balance = rd.fit(df, "age", "score")  # balance placebo: covariate as outcome; small p-value = imbalance, do not adjust
```

Canonical fields are ONE coherent row (the robust row): att = bias-corrected
estimate, se = robust SE. rdrobust prints the conventional estimate as its
headline - that is att_conventional here, with a full inference row of its own.
Fuzzy fits warn when the first-stage robust CI contains zero (weak
identification; R is silent) and raise R's exact error when the take-up
variable has no variation and no jump. Covariate-adjusted fits drop collinear
covariates with a warning naming them (covs_drop=True, R's default; False =
strict error) and guard degenerate adjustments (constant covariates are
excluded, full dummy sets take a stabilized cut - both warned; R silently
returns platform-dependent noise there). Cluster-robust variance, weights,
kink estimands, weak-IV-robust fuzzy inference, and a packaged
covariate-balance helper are documented follow-ups; missing rows are dropped
WITH a warning (R drops silently); N < 20 falls back to full-range bandwidths
exactly as rdrobust does (overriding manual h).

### RDPlot

RD plot diagnostic (Calonico, Cattaneo & Titiunik 2015), parity-targeting R
rdrobust 4.0.0's `rdplot()`. NOT a treatment-effect estimator: builds the
standard exploratory RD plot - per-side global order-p polynomial fits
(default p=4, uniform kernel over each side's full range; deliberately
different defaults from estimation) plus binned local sample means with a
DATA-DRIVEN number of bins per side. All 8 `binselect` selectors =
{evenly (es) / quantile (qs) spaced} x {IMSE-optimal / mimicking variance
(mv)} x {spacings / polynomial-regression (pr) variance estimators}; the
result reports `J`, `J_IMSE`, `J_MV`, the implied scale `rscale = J/J_IMSE`,
and the CCT 2015 Supplement S.1 WIMSE weights (variance `1/(1+rscale^3)`,
bias `rscale^3/(1+rscale^3)`), matching R's `summary.rdplot`.

```python
RDPlot(
    cutoff: float = 0.0,               # Known threshold of the running variable
    p: int = 4,                        # Global polynomial order (R's rdplot default)
    nbins: int | tuple | None = None,  # Manual bins per side (overrides the selector)
    binselect: str = "esmv",           # es/espr/esmv/esmvpr/qs/qspr/qsmv/qsmvpr
    scale: float | tuple | None = None,  # Multiplies the selected J (fractional products take CCT 2015 Eq 2's ceiling; R crashes there)
    kernel: str = "uniform",           # Weighting for the global fit within h
    h: float | tuple | None = None,    # Global-fit window; default = full per-side range
    support: tuple | None = None,      # (lower, upper); only WIDENS the observed range
    masspoints: str = "adjust",        # >= 20% ties: warn + switch spacings binselect to its *pr sibling ("check" warns only; "off" disables)
    ci: float | None = None,           # Per-bin CI level for plotted error bars (CI columns are ALWAYS computed, at 95 when None)
    covs_drop: bool = True,            # Drop collinear covariates with a warning; False = strict error
)
```

**fit() parameters:**

```python
rp.fit(
    data: pd.DataFrame,
    outcome: str,
    running: str,
    covariates: list[str] | None = None,  # Covariate-adjusted plot (R's covs=, covs_eval="mean"); gammas exposed as covariate_coefficients
) -> RDPlotResult
```

The 3.8 names `outcome_col=` and `running_col=` remain accepted as deprecated aliases (FutureWarning; removed in 4.0).

```python
from diff_diff import RDPlot

rp = RDPlot(cutoff=0.0)
res = rp.fit(df, outcome="y", running="score")
res.J, res.J_IMSE, res.J_MV     # selected / IMSE-optimal / mimicking-variance bins per side
res.rscale                       # implied scale; WIMSE weights via res.wimse_variance_weight / res.wimse_bias_weight
res.to_dataframe()               # vars_bins: R-named per-bin columns (rdplot_mean_bin/mean_x/mean_y/min_bin/max_bin/se_y/N/ci_l/ci_r)
res.vars_poly                    # 500-point global-fit curve per side
print(res.summary())             # R summary.rdplot layout incl. implied scale + WIMSE weights
res.plot()                       # optional matplotlib rendering (matplotlib NOT a dependency; ImportError if absent)
```

Empty bins are dropped from vars_bins (R's tapply) while bin_avg/bin_med
cover all partition cells; single-observation bins get se = 0 (R's NA -> 0
fixup); n < 20 raises outright (R behavior - no fallback); a side with
zero outcome variance pins its J to 1 with a warning and echoes R's Inf
J_IMSE. With empty LEFT-side bins, R's own rdplot_min_bin/rdplot_max_bin
come from mirrored grid slots and can disagree with rdplot_mean_bin - an
R quirk replicated for parity (rdplot_mean_bin is the reliable bin
identity). Sampling weights and R's subset= are documented seams (filter
the DataFrame instead).

### StackedDiD

Stacked DiD estimator (Wing, Freedman & Hollingsworth 2024). Addresses TWFE bias with corrective Q-weights.

```python
StackedDiD(
    kappa_pre: int = 1,                       # Pre-treatment event-time periods
    kappa_post: int = 1,                      # Post-treatment event-time periods
    weighting: str = "aggregate",             # "aggregate", "population", or "sample_share"
    control_group: str = "not_yet_treated",   # "not_yet_treated", "strict", or "never_treated"
    clean_control: str | None = None,         # DEPRECATED alias for control_group= (M-043; FutureWarning, removed in 4.0)
    cluster: str = "unit",                    # "unit" or "unit_subexp"
    alpha: float = 0.05,
    anticipation: int = 0,
    rank_deficient_action: str = "warn",
    vcov_type: str = "hc1",                   # {"hc1","hc2_bm"}; classical/hc2 rejected (intrinsically clustered), conley deferred. survey_design=... requires hc1
    balance: str = "none",                    # {"none","entropy"}; "entropy" = CBWSDID covariate balancing (Ustyuzhanin 2026), requires fit(covariates=[...]) + weighting="aggregate", no survey_design
    df_convention: str = "residual",          # Analytical t/p/CI df (3.9: hc1 lane moved z -> t(pooled residual df)): "residual" (default), "cluster" (G-1, positive-weight clusters), "normal" (z); BM/survey df keep precedence; flips at v4
)
```

**Alias:** `Stacked` (deprecated 3.9, removed 4.0 - emits FutureWarning; use `StackedDiD`)

**fit() parameters:**

```python
stacked.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-024): warns; surface always computed - use results.aggregate() post-fit
    population: str = None,        # Required when weighting="population"
    covariates: list[str] = None,  # Columns to balance (requires balance="entropy"); values read at t=a-1-anticipation; balanced windows only
) -> StackedDiDResults
```

**Usage:**

```python
from diff_diff import StackedDiD, plot_event_study

est = StackedDiD(kappa_pre=2, kappa_post=2)
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat')
results.print_summary()
plot_event_study(results)
es = results.aggregate('event_study')  # unified container (M-024)
simple = results.aggregate('simple')   # one-row overall relay
```

### EfficientDiD

Efficient DiD estimator (Chen, Sant'Anna & Xie 2025). Achieves the semiparametric efficiency bound for ATT(g,t) on the **no-covariate path**. Also supports an optional doubly-robust covariate path with all nuisances estimated nonparametrically (sieve-based propensity score ratios + a sieve outcome regression with AIC/BIC order selection + kernel-smoothed conditional covariance): the DR property gives consistency if either the outcome regression or the PS is correctly specified, and the covariate path attains the efficiency bound asymptotically under the paper's regularity conditions (a *growing* polynomial sieve — no fixed order ceiling — whose basis dimension satisfies Assumption C.1's rate for low-dimensional covariates; degree 1 reproduces a linear working model, and `sieve_k_max=1` forces all covariate-path sieves to degree 1). Pass column names to the `covariates` parameter on `fit()`.

```python
EfficientDiD(
    pt_assumption: str = "all",              # "all" (overidentified) or "post" (just-identified)
    alpha: float = 0.05,
    cluster: str | None = None,              # Column name for cluster-robust SEs (Liang-Zeger on EIF values); cluster-level multiplier bootstrap when n_bootstrap > 0
    vcov_type: str = "hc1",                  # {"hc1"} only — IF-based variance per Chen-Sant'Anna-Xie (2025); analytical-sandwich families {classical, hc2, hc2_bm} and conley rejected at __init__/set_params
    n_bootstrap: int = 0,                    # Multiplier bootstrap iterations
    bootstrap_weights: str = "rademacher",   # "rademacher", "mammen", or "webb"
    seed: int | None = None,
    anticipation: int = 0,
    omega_ridge: float = 1e-6,               # Ridge for the Omega* inversion behind the efficient weights (PT-All's overidentified moments make sample Omega* numerically singular); 0 = legacy exact-inverse/pseudoinverse path. See the Omega* ridge Note in the methodology registry.
)
```

**Alias:** `EDiD`

**fit() parameters:**

```python
edid.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    covariates: list[str] = None,  # Time-invariant unit-level covariates; uses doubly-robust sieve path when non-None
    aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-023): warns, still works - use results.aggregate() post-fit
    balance_e: Any = NOT_SUPPLIED, # DEPRECATED (M-120): moves onto results.aggregate(balance_e=)
) -> EfficientDiDResults
```

**Usage:**

```python
from diff_diff import EfficientDiD

edid = EfficientDiD(pt_assumption="all")
results = edid.fit(data, outcome='y', unit='id', time='t',
                   first_treat='first_treat')
results.print_summary()
# Aggregate post-fit (recomputed from retained EIFs; on bootstrapped fits
# the recompute levels raise while 'simple' relays the stored inference):
es = results.aggregate('event_study')       # EventStudyResults container
grp = results.aggregate('group')            # per-cohort AggregationResult
print(grp.to_dataframe())
```

### LPDiD

Local Projections DiD (Dube, Girardi, Jorda & Taylor 2025). Estimates a separate OLS at each event-time horizon of a long difference (`y_{i,t+h} - y_{i,t-1}`) on the treatment-switch indicator plus calendar-time fixed effects (no unit FE), restricted to a flexible "clean control" sample of newly-treated and not-yet-treated units. Excluding already-treated units from the control group removes the negative-weighting bias of naive TWFE, so the default (variance-weighted) estimand has strictly non-negative weights. `reweight=True` yields the equally-weighted ATT (numerically equivalent to Callaway-Sant'Anna); covariates then enter via regression adjustment. Standard errors on the default/weighted path are cluster-robust at the unit level (the paper specifies no SE; matches Stata `lpdid` `vce(cluster unit)`); the regression-adjustment covariate path (`reweight=True`) instead reports an influence-function cluster variance (ImputationDiD/BJS family). Scope: binary treatment; absorbing by default (rejects panels where treatment turns off), with non-absorbing (reversible) treatment available via `non_absorbing` - `"first_entry"` (Dube et al. Eq. 12, the effect of entering for the first time and staying treated) or `"effect_stabilization"` (Eq. 13, requires `stabilization_window=L`; lets units whose treatment has been stable for at least `L` periods act as clean controls, so estimation is feasible with few/no never-treated units). Non-absorbing modes require a gap-free panel within each unit's observed span. Complex-survey designs are supported on the variance-weighted default path via the `survey_design=` argument to `fit()` (probability weights enter the WLS point estimate; the SE is the stratified-PSU Taylor-linearization sandwich with `df = n_PSU - n_strata`, with optional FPC and lonely-PSU handling) — rejected with `reweight=True`, replicate weights, or non-pweight types.

```python
LPDiD(
    pre_window: int = 2,                     # Number of pre-treatment horizons (placebos)
    post_window: int = 0,                    # Number of post-treatment horizons
    control_group: str = "clean",            # "clean" (not-yet-treated) or "never_treated"
    reweight: bool = False,                  # True -> equally-weighted ATT (== Callaway-Sant'Anna); False -> variance-weighted
    no_composition: bool = False,            # Hold the post-treatment composition fixed across post horizons
    pmd: str | int | None = None,            # Base period: None=first-lag (t-1), "max"=premean over all pretreatment periods, int=last-k premean
    alpha: float = 0.05,
    cluster: str | None = None,              # Cluster column for cluster-robust SEs; defaults to the unit identifier
    rank_deficient_action: str = "warn",     # "warn", "error", or "silent"
    non_absorbing: str | None = None,        # None=absorbing; "first_entry" (Eq. 12); "effect_stabilization" (Eq. 13)
    stabilization_window: int | None = None, # The paper's L; required when non_absorbing="effect_stabilization"
    df_convention: str = "cluster",          # DEFAULT "cluster" = the existing Stata lpdid t(G-1) reference (nothing moved in 3.9); "residual" (per-horizon n_eff-k; RA path n_total-k0-1) and "normal" (z) opt in
)
```

**fit() parameters:**

```python
lpdid.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    treatment: str,                    # Binary treatment indicator (0/1); absorbing unless non_absorbing is set
    covariates: list[str] = None,      # Direct inclusion (reweight=False) or regression adjustment (reweight=True)
    ylags: int = 0,                    # Lagged-outcome controls
    dylags: int = 0,                   # Lagged first-difference controls
    absorb: list[str] = None,          # Additional absorbed fixed-effect columns
    post_pooled: int | tuple = None,   # Pooled post-window horizons (int or (start, end))
    pre_pooled: int | tuple = None,    # Pooled pre-window horizons (int or (start, end))
    only_event: bool = False,          # Compute only the event-study table
    only_pooled: bool = False,         # Compute only the pooled pre/post table
    survey_design: SurveyDesign = None,  # Complex-survey design (pweight + optional strata/PSU/FPC); variance-weighted default path only (rejected with reweight=True)
) -> LPDiDResults
```

**Usage:**

```python
from diff_diff import LPDiD

lp = LPDiD(pre_window=5, post_window=10)
results = lp.fit(data, outcome='y', unit='id', time='t', treatment='treated')
results.print_summary()
print(results.event_study)   # per-horizon coefficients
print(results.pooled)        # pooled pre (placebo) / post (ATT) rows
```

### ChangesInChanges

Changes-in-Changes (Athey & Imbens 2006) for the canonical 2x2 design with continuous outcomes. Recovers the treated group's full counterfactual outcome distribution `F_10(F_00^{-1}(F_01(y)))` and reports the ATT plus quantile treatment effects on a grid (default 0.05-0.95 by 0.05, matching R `qte`). The model is invariant to monotone transformations of the outcome (exact for unconditional fits; the covariate QR branch is not equivariant to nonlinear monotone transforms) and, with continuous data, places no testable restrictions. Point estimation matches `qte::CiC()` (v1.3.1) exactly: R type-1 (ceiling-order-statistic) quantiles - the paper's empirical-inverse convention - throughout, no smoothing anywhere. Inference is bootstrap-only in this release: panel mode resamples units (both periods together), repeated cross-section mode draws a pooled row resample; SEs are replicate SDs with symmetric normal-approximation CIs, plus a sup-t critical value for uniform bands at a FIXED 95% level (qte parity - does not follow `alpha`). Quantile effects outside the point-identified interior range keep their point estimates but report NaN inference with a warning (unconditional fits; with covariates the interior bounds do not apply and `q_lower`/`q_upper` are NaN). Covariates (`covariates=[...]` or trailing formula terms) port qte's `xformla` branch exactly: linear quantile regressions in the control cells on qte's fixed internal 0.01-0.99 tau grid, conditional-rank imputation per treated pre-period observation, a conditional-envelope support warning (Melly-Santangelo Assumption 4), and quantile regressions refit inside every bootstrap replicate (~40k small LPs at n_bootstrap=200 - tens of seconds, same cost profile as qte). Covariates must be numeric (dummy-encode categoricals). Deferred (documented in REGISTRY.md): the full Melly-Santangelo covariate estimator, discrete-outcome bounds (a ties warning fires on discrete-looking outcomes), analytical SEs, staggered designs, treatment-on-controls. Note: additive random group-time shocks BIAS CiC (not just its inference) and are undetectable in a 2x2 design (Athey-Imbens p. 476).

```python
ChangesInChanges(
    quantiles: array-like | None = None,  # Quantile grid strictly inside (0,1); None -> 0.05..0.95 step 0.05 (qte default)
    n_bootstrap: int = 200,               # Bootstrap replicates; 0 disables inference (NaN se/t/p/CI)
    alpha: float = 0.05,                  # Pointwise CI level (uniform bands stay fixed at 95%)
    panel: bool = False,                  # True: same units both periods (requires unit=); affects resampling only
    seed: int | None = None,              # Bootstrap RNG seed (numpy default_rng)
    method: str = 'cic',                  # KEYWORD-ONLY. 'cic' | 'qdid' - which 2x2 distributional estimator to fit
)
```

**fit() parameters:**

```python
cic.fit(
    data: pd.DataFrame,
    outcome: str = None,       # Continuous outcome column
    treatment: str = None,     # Binary group indicator (1 = treated group in BOTH periods)
    time: str = None,          # Binary post-period indicator
    formula: str = None,       # Alternative: "y ~ treated * post [+ x1 + x2]" (trailing terms = covariates; formula + explicit kwargs raises)
    covariates: list[str] = None,  # Numeric covariate columns (qte xformla parity; fit-time arg, not in get_params)
    unit: str = None,          # Unit id; required when panel=True, ignored (documented) otherwise
) -> ChangesInChangesResults
```

**Alias:** `CiC`

**Usage:**

```python
from diff_diff import ChangesInChanges

cic = ChangesInChanges(n_bootstrap=200, seed=42)
results = cic.fit(data, outcome='y', treatment='treated', time='post')
results.print_summary()
print(results.quantile_effects)   # per-quantile QTE table
print(results.uniform_bands())    # sup-t simultaneous bands (fixed 95%)
```

### QDiD

Quantile Difference-in-Differences comparison estimator (Athey & Imbens 2006, Section 3.3) for the 2x2 design: `QTE(tau) = Q(y11,tau) - [Q(y10,tau) + Q(y01,tau) - Q(y00,tau)]` with R type-7 linear-interpolation quantiles, matching `qte::QDiD()` (v1.3.1) exactly - including its ATT formula (control-group quantile functions evaluated at treated pre-period own-sample ranks; population-equivalent to the paper's k^QDID transformation but a different finite-sample estimator, see the REGISTRY.md Note). The paper recommends ChangesInChanges over QDiD: QDiD's justifying model is not scale-invariant, forces identical unobservable distributions in all four cells, and places testable restrictions on the data (a warning fires when the implied counterfactual quantile function is non-monotone; unconditional fits only - the covariate-path counterfactual quantile curve is monotone by construction). QDiD's mean effect equals standard DiD's ATT in population. **Deprecated in 3.9, removed in 4.0 - use `ChangesInChanges(method="qdid")` instead** (same engine, identical numbers; only the class spelling is deprecated, the estimator is not). Fit signature (including `covariates=`), bootstrap machinery, and results container are identical to ChangesInChanges (no interior-range guard; `q_lower`/`q_upper` are NaN); the constructor is identical EXCEPT that the merged class additionally carries `method=`, which is what selects this estimator there. Covariate fits use quantile regressions in THREE cells with own-cell conditional ranks and qte's verbatim-ported asymmetric quantile types (type-7 treated-post, type-1 imputed counterfactual).

```python
QDiD(
    quantiles: array-like | None = None,  # Same defaults as ChangesInChanges
    n_bootstrap: int = 200,
    alpha: float = 0.05,
    panel: bool = False,
    seed: int | None = None,
)
```

**Usage:**

```python
from diff_diff import ChangesInChanges

qdid = ChangesInChanges(method='qdid', n_bootstrap=200, seed=42)
results = qdid.fit(data, outcome='y', treatment='treated', time='post')
```

### TROP

Triply Robust Panel estimator (Athey, Imbens, Qu & Viviano 2025). Combines nuclear norm regularization, distance-based unit weights, and time decay weights.

```python
TROP(
    method: str = "local",                     # "local" or "global"
    lambda_time_grid: list[float] = None,     # Time weight decay grid [0, 0.1, 0.5, 1, 2, 5]
    lambda_unit_grid: list[float] = None,     # Unit weight decay grid [0, 0.1, 0.5, 1, 2, 5]
    lambda_nn_grid: list[float] = None,       # Nuclear norm grid [0, 0.01, 0.1, 1, 10]
    max_iter: int = 100,
    tol: float = 1e-6,
    alpha: float = 0.05,
    n_bootstrap: int = 200,
    seed: int | None = None,
    non_absorbing: bool = False,               # False: require absorbing D (reject non-monotonic). True: allow on/off treatment (Eq. 12/Alg. 2), method='local' only; emits a caveat warning (Thm 5.1 is block-only).
)
```

**fit() parameters:**

```python
trop.fit(
    data: pd.DataFrame,
    outcome: str,
    treatment: str,                # Treatment indicator (0/1). Default (non_absorbing=False): absorbing state -- 0 for all pre-treatment periods, 1 for treatment and post-treatment; non-monotonic D raises ValueError. With non_absorbing=True: any on/off pattern (general assignment).
    unit: str,
    time: str,
) -> TROPResults
```

**Usage:**

```python
from diff_diff import TROP

trop = TROP(method='local', seed=42)
results = trop.fit(data, outcome='outcome', treatment='treated',
                   unit='unit', time='period')
results.print_summary()
```

### BaconDecomposition

Goodman-Bacon (2021) decomposition of TWFE into 2x2 DiD comparisons.

```python
BaconDecomposition(
    weights: str = "approximate",   # "approximate" or "exact"
)
```

**Alias:** `Bacon`

**fit() parameters:**

```python
bacon.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
) -> BaconDecompositionResults
```

**Usage:**

```python
from diff_diff import BaconDecomposition, plot_bacon

bacon = BaconDecomposition(weights="exact")
results = bacon.fit(data, outcome='outcome', unit='unit',
                    time='period', first_treat='first_treat')
results.print_summary()
plot_bacon(results)
```

### StaggeredTripleDifference

DEPRECATED in 3.9, removed in 4.0 (ledger row M-013). Use
`TripleDifference().fit(..., unit=, time=, first_treat=, partition=)` - the same
engine, so the numbers are identical. Vocabulary on the merged surface:
`eligibility=` is `partition=`, and `control_group` takes the underscored values
`not_yet_treated`/`never_treated` (this class keeps R's compact `notyettreated`/
`nevertreated` until removal). `cluster=` raises on the merged surface instead of
being accepted-and-ignored. Alias `SDDD` is deprecated with the class.

Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD estimator for designs with two eligibility criteria and staggered treatment timing.

```python
StaggeredTripleDifference(
    estimation_method: str = "dr",          # "dr", "ipw", or "reg"
    control_group: str = "notyettreated",   # "nevertreated" or "notyettreated"
    alpha: float = 0.05,
    anticipation: int = 0,
    base_period: str = "varying",           # "varying" or "universal"
    n_bootstrap: int = 0,
    bootstrap_weights: str = "rademacher",
    seed: int | None = None,
    cband: bool = True,
    pscore_trim: float = 0.01,
    cluster: str | None = None,
    rank_deficient_action: str = "warn",
    epv_threshold: float = 10,              # Min events-per-variable for propensity score
    pscore_fallback: str = "error",         # "error" or "unconditional"
)
```

**Alias:** `SDDD`

**fit() parameters:**

```python
sddd.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,
    eligibility: str,               # Binary eligibility column (0/1)
    covariates: list[str] | None = None,
    aggregate: str | None = None,   # None, "simple", "group", "event_study", "all"
    balance_e: int | None = None,   # Max event time for balanced event study
    survey_design=None,
) -> StaggeredTripleDiffResults
```

**Usage:**

```python
from diff_diff import StaggeredTripleDifference

sddd = StaggeredTripleDifference(estimation_method='dr', n_bootstrap=999)
results = sddd.fit(data, outcome='y', unit='id', time='t',
                   first_treat='ft', eligibility='eligible',
                   aggregate='event_study')
results.print_summary()
```

### WooldridgeDiD

Wooldridge (2023, 2025) Extended Two-Way Fixed Effects (ETWFE) estimator. OLS path uses direct saturated-regression coefficients for ATT(g,t). Logit and Poisson QMLE paths use Average Structural Function (ASF) based ATT with delta-method standard errors.

```python
WooldridgeDiD(
    method: str = "ols",                    # "ols", "logit", or "poisson"
    control_group: str = "not_yet_treated", # "not_yet_treated" or "never_treated"
    anticipation: int = 0,                  # Number of anticipation periods
    demean_covariates: bool = True,         # Demean covariates within cohort*period cells
    alpha: float = 0.05,
    cluster: str | None = None,
    n_bootstrap: int = 0,
    bootstrap_weights: str = "rademacher",
    seed: int | None = None,
    rank_deficient_action: str = "warn",
    vcov_type: str = "hc1",                 # {"classical","hc1","hc2","hc2_bm","conley"}; OLS path only.
                                            # hc1 (default) preserves prior bit-equal within-transform CR1.
                                            # hc2_bm auto-routes to full-dummy + clubSandwich WLS-CR2 algebra.
                                            # classical/hc2 auto-drop the unit auto-cluster (one-way only);
                                            # explicit cluster="X" + one-way raises at the linalg validator.
                                            # conley (Conley-1999 spatial-HAC, OLS only): within-transform via solve_ols — pass conley_coords=(lat,lon), conley_cutoff_km, conley_lag_cutoff; survey/weights/n_bootstrap rejected. method != "ols" requires hc1 (QMLE CR2-BM TBD).
                                            # survey_design= requires hc1 (survey TSL overrides analytical).
    conley_coords: tuple[str, str] | None = None,  # (lat_col, lon_col); required for vcov_type="conley"
    conley_cutoff_km: float | None = None,         # km bandwidth; required for vcov_type="conley"
    conley_metric: str = "haversine",              # "haversine" | "euclidean"
    conley_kernel: str = "bartlett",               # "bartlett" | "uniform"
    conley_lag_cutoff: int | None = None,          # within-unit Bartlett max lag (0 = spatial-only)
    cohort_trends: bool = False,            # Linear dg_i*t cohort trends (W2025 Sec 8 / Eq 8.1); OLS only.
                                            # On all-eventually-treated panels the last cohort's trend
                                            # column is dropped per Sec 5.4, so cohort_trend_coefs
                                            # carries G-1 entries. Rejected with survey_design= and
                                            # with control_group="never_treated".
    df_convention: str = "residual",        # OLS analytical t/p/CI df (3.9: default-hc1 arms moved z -> t(residual)): "residual" (default), "cluster" (G-1, hc1-clustered only), "normal" (z); survey/BM DOF keep precedence; GLM arms knob-independent (explicit non-default warns); flips at v4
)
```

`fit()` additionally accepts `survey_design=` (a `SurveyDesign`) on all three
methods; see the Survey Support section.

**All-eventually-treated panels (no never-treated group).** Use
`control_group="not_yet_treated"` — `"never_treated"` raises when no cohort-0
units exist. Periods at which every unit is treated carry no identified
ATT(g, t), so they are REMOVED from the estimation sample before the solve and
the last cohort becomes the reference (W2025 Section 5.4). The fit warns naming
the dropped periods, the observation count and the cause, and separately names
any cohort left with no estimated cells; `results.groups` excludes those
cohorts. Stata `jwdid` performs the same reduction silently, reporting only a
smaller `N`. SOME covariate specifications on such a panel are still
rank-deficient -- `exovar`, `xgvar`, and `xtvar` with
`demean_covariates=False` -- because `D_{G_max} x X` is not normalized;
`rank_deficient_action="error"` raises on those and coefficients are
unaffected either way. Default `xtvar` (`demean_covariates=True`) is FULL
RANK and fits cleanly.

**Survey designs refuse row-dropping paths.** `survey_design=` combined with
either comparison-support period filtering or unidentified-cohort exclusion
raises `NotImplementedError` rather than deleting rows, because deletion would
remove their PSUs and strata from the TSL variance. Restrict the frame yourself
and re-fit. The refusals are conditional: survey fits that drop nothing are
unaffected.

**Alias:** `ETWFE`

**fit() parameters:**

```python
etwfe.fit(
    data: pd.DataFrame,
    outcome: str,
    unit: str,
    time: str,
    first_treat: str,               # First treatment period (0 or NaN = never treated)
    exovar: list[str] | None = None,  # Time-invariant covariates (no interaction/demeaning)
    xtvar: list[str] | None = None,   # Time-varying covariates (demeaned within cohort*period)
    xgvar: list[str] | None = None,   # Covariates interacted with each cohort indicator
) -> WooldridgeDiDResults
```

The 3.8 name `cohort=` remains accepted as a deprecated alias for `first_treat=` (FutureWarning; removed in 4.0).

**Aggregation types:**

```python
results.aggregate("simple")    # Overall ATT
results.aggregate("group")     # ATT by cohort
results.aggregate("calendar")  # ATT by calendar period
results.aggregate("event_study")  # ATT by event time (relative to treatment)
```

**Usage:**

```python
from diff_diff import WooldridgeDiD

# OLS (linear outcomes). On a binary/count outcome this emits a UserWarning that
# logit/poisson is often the more appropriate specification -- link-scale (not
# level) parallel trends, and less biased/more precise in Wooldridge (2023) sims;
# a different identifying assumption, so a recommended comparison not a free switch.
etwfe = WooldridgeDiD(method='ols')
results = etwfe.fit(data, outcome='y', unit='id', time='t', first_treat='first_treat')
print(results.aggregate("simple"))

# Logit (binary outcomes)
etwfe_logit = WooldridgeDiD(method='logit')
results = etwfe_logit.fit(data, outcome='y_bin', unit='id', time='t', first_treat='first_treat')

# Poisson (count outcomes)
etwfe_pois = WooldridgeDiD(method='poisson')
results = etwfe_pois.fit(data, outcome='y_count', unit='id', time='t', first_treat='first_treat')
```

### Convenience Functions (deprecated 3.9, removed 4.0)

All 8 wrapper functions (`imputation_did`, `two_stage_did`, `stacked_did`, `trop`, `synthetic_control`, `triple_difference`, `bacon_decompose`, `chaisemartin_dhaultfoeuille`) emit a `FutureWarning` since 3.9 and are removed in 4.0 (rows M-070..M-077) - construct the class instead:

```python
# Canonical class form (what the wrappers did internally)
from diff_diff import ImputationDiD, TwoStageDiD, TripleDifference, StackedDiD, TROP, BaconDecomposition

results = ImputationDiD().fit(data, outcome='y', unit='id', time='t', first_treat='ft')
results = TwoStageDiD().fit(data, outcome='y', unit='id', time='t', first_treat='ft')
results = TripleDifference().fit(data, outcome='y', group='g', partition='p', post='t')
results = StackedDiD(kappa_pre=2, kappa_post=2).fit(data, outcome='y', unit='id', time='t',
                                                    first_treat='ft')
results = TROP().fit(data, 'y', 'd', 'id', 't')
results = BaconDecomposition().fit(data, outcome='y', unit='id', time='t', first_treat='ft')
```

## Results Objects

**Flat-alias compatibility note.** Every staggered result class in this
section (those with canonical `overall_*` / `overall_att_*` / `avg_*`
prefixed inference fields) ALSO exposes the unprefixed flat names
`att` / `se` / `conf_int` / `p_value` / `t_stat` as read-only `@property`
aliases over the canonical fields. The canonical prefixed fields remain
the documented and computed surface; the flat aliases are pure
read-throughs for compatibility with external adapters that
`getattr(res, "se", None)`-style query the inference surface (e.g.
`balance.interop.diff_diff.as_balance_diagnostic()`). Tables below list
the canonical names; assume the flat aliases are present on every
staggered class unless explicitly noted otherwise.

**balance interop.** Meta's `balance` package (>= 0.21) ships the
one-way adapter `balance.interop.diff_diff` (`pip install "balance[did]"`,
pins `diff-diff>=3.3,<4`): `to_survey_design(sample)` builds a
`SurveyDesign` from a balance `Sample`'s active weight column plus the
convention columns `stratum`/`psu`/`fpc`; `to_panel_for_did(sample, by=,
outcomes=)` wraps `diff_diff.aggregate_survey` to collapse respondent
microdata into a unit-period panel plus second-stage design;
`fit_did(sample, estimator=, ...)` resolves any exported estimator by
name or short alias (`CS`/`DiD`/`BJS`/`HAD`) and forwards
`survey_design=`, attaching the source Sample to the result as
`_balance_adjustment` for provenance; `as_balance_diagnostic(sample,
res)` joins balance's ASMD/Kish-ESS with `res.survey_metadata`'s
DEFF/effective-n into one flat dict. The diff-diff surface it consumes
is pinned by `tests/test_balance_interop_contract.py`; the workflow is
demonstrated in Tutorial 26 (composition drift & calibration).

### DiDResults

Returned by `DifferenceInDifferences.fit()` and `TwoWayFixedEffects.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | Average Treatment effect on the Treated |
| `se` | `float` | Standard error of ATT |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value (H0: ATT = 0) |
| `conf_int` | `tuple[float, float]` | Confidence interval |
| `n_obs` | `int` | Number of observations |
| `n_treated` | `int` | Number of treated units |
| `n_control` | `int` | Number of control units |
| `alpha` | `float` | Significance level |
| `coefficients` | `dict` | All regression coefficients |
| `vcov` | `np.ndarray` | Variance-covariance matrix |
| `residuals` | `np.ndarray` | Regression residuals |
| `fitted_values` | `np.ndarray` | Fitted values |
| `r_squared` | `float` | R-squared |
| `inference_method` | `str` | "analytical" or "wild_bootstrap" |
| `n_bootstrap` | `int` | Number of bootstrap replications |
| `n_clusters` | `int` | Number of clusters |
| `bootstrap_distribution` | `np.ndarray` | Bootstrap ATT distribution |

**Methods:** `summary(alpha=None)`, `print_summary()`, `to_dict()`, `to_dataframe()`

**Properties:** `is_significant`, `significance_stars`

### MultiPeriodDiDResults

Returned by `MultiPeriodDiD.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `period_effects` | `dict[Any, PeriodEffect]` | Period-specific effects (pre and post) |
| `avg_att` | `float` | Average ATT across post-periods |
| `avg_se` | `float` | SE of average ATT |
| `avg_t_stat` | `float` | T-statistic for average ATT |
| `avg_p_value` | `float` | P-value for average ATT |
| `avg_conf_int` | `tuple[float, float]` | CI for average ATT |
| `n_obs` | `int` | Number of observations |
| `n_treated` | `int` | Number of treated observations |
| `n_control` | `int` | Number of control observations |
| `pre_periods` | `list` | Pre-treatment period identifiers |
| `post_periods` | `list` | Post-treatment period identifiers |
| `reference_period` | `Any` | Reference (omitted) period |
| `r_squared` | `float` | R-squared |
| `vcov` | `np.ndarray` | Variance-covariance matrix |
| `interaction_indices` | `dict` | Period to VCV column index mapping |

**Methods:** `summary()`, `print_summary()`, `get_effect(period)`, `to_dict()`, `to_dataframe()`

**Properties:** `pre_period_effects`, `post_period_effects`, `is_significant`, `significance_stars`

### PeriodEffect

Individual period treatment effect (used in MultiPeriodDiDResults).

| Attribute | Type | Description |
|-----------|------|-------------|
| `period` | `Any` | Time period identifier |
| `effect` | `float` | Treatment effect estimate |
| `se` | `float` | Standard error |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | Confidence interval |

**Properties:** `is_significant`, `significance_stars`

### CallawaySantAnnaResults

Returned by `CallawaySantAnna.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `group_time_effects` | `dict[(g,t), GroupTimeEffect]` | ATT(g,t) for each (group, time) |
| `overall_att` | `float` | Overall ATT |
| `overall_se` | `float` | SE of overall ATT |
| `overall_t_stat` | `float` | T-statistic |
| `overall_p_value` | `float` | P-value |
| `overall_conf_int` | `tuple[float, float]` | CI for overall ATT |
| `groups` | `list` | Treatment cohorts |
| `time_periods` | `list` | All time periods |
| `n_obs` | `int` | Number of observations |
| `event_study_effects` | `dict[int, dict]` | Event study effects by relative time. Populated only by the DEPRECATED fit-time `aggregate=`; `None` after a plain fit - use `aggregate("event_study")` instead |
| `group_effects` | `dict` | Group-level aggregated effects. Same: populated only by fit-time `aggregate=` |

**Methods:** `aggregate(type, weights=None, *, balance_e=None)`, `summary()`, `print_summary()`, `to_dataframe(level="event_study"|"group_time"|"group")`

`aggregate("simple"|"event_study"|"group")` re-aggregates POST-FIT with no
refit, returning a NEW object and leaving the result unchanged.
`"event_study"` returns `EventStudyResults`; the others return
`AggregationResult`. `balance_e=` applies to `"event_study"` only. Raises on
`"calendar"` (CS has no calendar aggregator) and, on a bootstrapped fit, on
the recompute levels (`"event_study"`/`"group"`) rather than substituting
analytical inference for percentile-bootstrap statistics - `"simple"` relays
the stored bootstrap inference with a NaN df column (the per-level rule).

### SunAbrahamResults

Returned by `SunAbraham.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `event_study_effects` | `dict[int, dict]` | Effects by relative time |
| `overall_att` | `float` | Overall ATT |
| `overall_se` | `float` | SE of overall ATT |
| `overall_t_stat` | `float` | T-statistic |
| `overall_p_value` | `float` | P-value |
| `overall_conf_int` | `tuple[float, float]` | CI |
| `cohort_weights` | `dict[int, dict]` | Interaction weights per period |
| `groups` | `list` | Treatment cohorts |
| `n_obs` | `int` | Number of observations |
| `n_treated_units` | `int` | Number of ever-treated units |
| `n_control_units` | `int` | Number of never-treated units |
| `control_group` | `str` | Control group type used |
| `cohort_effects` | `dict` | Cohort-level effects |
| `df_convention` | `str | None` | The configured df convention echoed onto the results ("residual" | "cluster" | "normal"; 3.9) |
| `inference_df` | `float | None` | The df the overall-ATT p/CI actually used (BM contrast df / survey df / knob-resolved fallback; None under bootstrap or "normal") |
| `event_study_df` | `dict[int, float] | None` | Per-event-time df provenance (finite residual df on plain analytic fits since 3.9; None under bootstrap) |

**Methods:** `summary()`, `print_summary()`, `to_dataframe(level="event_study"|"cohort")`

### SyntheticDiDResults

Returned by `SyntheticDiD.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | Average Treatment effect on the Treated |
| `se` | `float` | Standard error (bootstrap or placebo-based) |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | Confidence interval |
| `n_obs` | `int` | Number of observations |
| `n_treated` | `int` | Number of treated units |
| `n_control` | `int` | Number of control units |
| `unit_weights` | `dict` | Control unit synthetic weights |
| `time_weights` | `dict` | Pre-treatment time weights |
| `pre_periods` | `list` | Pre-treatment periods |
| `post_periods` | `list` | Post-treatment periods |
| `variance_method` | `str` | "bootstrap", "jackknife", or "placebo" |
| `variance_effects` | `np.ndarray` | Per-iteration draws (placebo effects, bootstrap ATT draws, or jackknife LOO estimates per `variance_method`); deprecated alias `placebo_effects` (removed v4.0.0) |
| `noise_level` | `float` | Estimated noise level |
| `zeta_omega` | `float` | Unit weight regularization |
| `zeta_lambda` | `float` | Time weight regularization |
| `pre_treatment_fit` | `float` | Pre-treatment RMSE |

**Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `get_unit_weights_df()`, `get_time_weights_df()`

**Validation diagnostics** (call after `fit()`):
- `get_weight_concentration(top_k=5)` - effective N and top-k weight share; flags fragile synthetic controls dominated by a few donor units
- `get_loo_effects_df()` - per-unit leave-one-out influence from the jackknife pass (DataFrame includes both control and treated rows). Requires `variance_method="jackknife"` with unit-level LOO granularity: available on non-survey and pweight-only jackknife fits; raises `NotImplementedError` on full-design survey jackknife (PSU-level LOO, see `result.variance_effects` for raw PSU-level replicates) and `ValueError` when LOO is unavailable (single treated unit, only one control with nonzero effective weight, etc.)
- `in_time_placebo()` - re-estimate on shifted fake treatment dates in the pre-period; near-zero placebo ATTs indicate a credible design
- `sensitivity_to_zeta_omega()` - re-estimate across a grid of unit-weight regularization values; checks ATT robustness to the auto-selected zeta_omega

### SyntheticControlResults

Returned by `SyntheticControl.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | Mean post-period gap (reported point estimate) |
| `se`, `t_stat`, `p_value`, `conf_int` | `float` / tuple | Always NaN — no analytical SE (use `in_space_placebo()`) |
| `placebo_p_value` | `float` | In-space placebo permutation p-value, `rank/(n_placebos+1)` (NaN until `in_space_placebo()` runs) |
| `rmspe_ratio` | `float` | Treated unit's post/pre RMSPE ratio = sqrt(post-MSPE/pre-MSPE) (the placebo test statistic; set at fit) |
| `n_placebos`, `n_failed` | `int` | Placebos in the reference set / excluded for non-convergence |
| `n_obs` | `int` | Treated + donor rows over all periods |
| `n_donors` | `int` | Donor units in the (post-filter) pool |
| `n_pre_periods`, `n_post_periods` | `int` | Period counts |
| `donor_weights` | `dict` | `{donor_id: weight}` on the simplex (near-zeros dropped) |
| `v_weights` | `dict` | `{predictor_label: v}`, trace-normalized |
| `predictor_balance` | `pd.DataFrame` | treated vs synthetic vs donor-mean per predictor |
| `gap_path` | `dict` | `{period: gap}` for all periods |
| `pre_rmspe` | `float` | Pre-treatment fit diagnostic |
| `mspe_v` | `float \| None` | Selected-V objective: pre-period MSPE (nested) or held-out validation MSPE (cv); None for custom / inverse_variance |
| `treated_unit` | `Any` | Treated unit identifier |
| `pre_periods`, `post_periods` | `list` | Calendar-sorted periods |
| `v_method`, `standardize` | `str` | Echoed configuration |
| `v_cv_t0` | `int \| None` | Resolved cv train/validation split index (None unless v_method="cv") |
| `effect_confidence_set` | `dict \| None` | Test-inversion confidence-set summary (Firpo-Possebom 2018 §4): `{family, parameter, gamma, lower, upper, contiguous, status ("ran"/"empty"/"unbounded"), ...}`; None until `confidence_set()` runs. SEPARATE from the always-NaN analytical `conf_int` (a permutation set at level 1−gamma, possibly a set/unbounded). |
| `conformal_inference` | `dict \| None` | CWZ (2021) conformal-inference summary of the most recent run: `{kind ("joint"/"pointwise"/"average"), scheme, status, n_perms, ...}` (joint adds `joint_p_value`/`proxy_converged`; average/pointwise add `alpha` + CI fields); None until a `conformal_*()` method runs. SEPARATE from the always-NaN analytical `conf_int`. |

**Methods:** `in_space_placebo()` (opt-in permutation inference; refits one synthetic control per donor), `get_placebo_df()` (per-unit RMSPE-ratio table incl. the treated row), `leave_one_out()` (ADH-2015 §4 donor robustness; drops each reportably-weighted donor (weight > 1e-6) → per-drop ATT/`delta_att` table) + `get_leave_one_out_df()`/`get_leave_one_out_gaps()`, `in_time_placebo()` (ADH-2015 §4 backdating placebo; reassigns the intervention earlier, TRUNCATE windowing, placebo ATT ~0 if no real pre-effect) + `get_in_time_placebo_df()`/`get_in_time_placebo_gaps()`, `regression_weights()` (ADH-2015 §4 regression-weight extrapolation diagnostic; intercept-augmented `W^reg`, flags donors outside [0,1]; pure linear algebra) + `get_regression_weights_df()`, `sparse_synthetic_control(sizes=None, max_subsets=50000)` (ADH-2015 §4 sparse subset search; exhaustive `C(J,l)` subsets holding V fixed; default-skip vs explicit-raise cap) + `get_sparse_synthetic_control_df()`/`get_sparse_synthetic_control_gaps()`, `test_sharp_null(effect, gamma=0.1)` (Firpo-Possebom 2018 §4: test a sharp null α_1t=f(t) by re-ranking the in-space placebo gaps — `effect` is a scalar or a post-period array; `test_sharp_null(0)` is identically `placebo_p_value`), `confidence_set(family="constant"|"linear", gamma=0.1, bounds=None, n_grid=200)` (invert that test for a confidence set of the effect path — a constant-effect interval / linear-slope set; strict p>gamma membership; exact piecewise-constant breakpoint inversion when bounds=None, else a fixed grid; `conf_int` stays NaN) + `get_confidence_set_df()`, `conformal_test(effect, q=1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, seed=None)` (CWZ 2021 joint sharp-null conformal p-value — fits its own constrained-LS proxy under the null on all periods, permutes residuals over time; `q∈{1,2,∞}`; one-sided `alternative` = signed statistic per Remark 1 with q=1; `covariates=` stacks pivoted-variable matching rows into the proxy), `conformal_confidence_intervals(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", bounds=None, n_grid=100, seed=None)` (pointwise per-period CIs, Algorithm 1 — each period uses Z=(pre, t); one-sided → half-lines), `conformal_average_effect(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", bounds=None, n_grid=200, seed=None)` (average-effect CI by T*-block collapse, Appendix A.1; one-sided → half-line; covariate rows block-collapse identically) + `get_conformal_grid_df()`, `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `get_gap_df()`, `get_weights_df()`

### TripleDifferenceResults

Returned by `TripleDifference.fit()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | ATT estimate |
| `se` | `float` | Standard error |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | Confidence interval |
| `n_obs` | `int` | Total observations |
| `n_treated_eligible` | `int` | Treated + eligible count |
| `n_treated_ineligible` | `int` | Treated + ineligible count |
| `n_control_eligible` | `int` | Control + eligible count |
| `n_control_ineligible` | `int` | Control + ineligible count |
| `estimation_method` | `str` | "dr", "reg", or "ipw" |
| `group_means` | `dict` | Cell means |
| `pscore_stats` | `dict` | Propensity score diagnostics |
| `r_squared` | `float` | R-squared (for "reg") |

**Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`

### BaconDecompositionResults

Returned by `BaconDecomposition.fit()` (and the deprecated `bacon_decompose()` wrapper, removed in 4.0).

| Attribute | Type | Description |
|-----------|------|-------------|
| `twfe_estimate` | `float` | Overall TWFE coefficient |
| `comparisons` | `list[Comparison2x2]` | All 2x2 comparisons |
| `total_weight_treated_vs_never` | `float` | Weight on treated vs never-treated |
| `total_weight_earlier_vs_later` | `float` | Weight on earlier vs later |
| `total_weight_later_vs_earlier` | `float` | Weight on forbidden comparisons |
| `weighted_avg_treated_vs_never` | `float` | Avg effect from clean comparisons |
| `weighted_avg_earlier_vs_later` | `float` | Avg effect from earlier vs later |
| `weighted_avg_later_vs_earlier` | `float` | Avg effect from forbidden comparisons |
| `n_timing_groups` | `int` | Number of treatment timing groups |
| `n_never_treated` | `int` | Number of never-treated units |
| `timing_groups` | `list` | Treatment timing cohorts |
| `n_obs` | `int` | Total observations |
| `decomposition_error` | `float` | Error: TWFE minus weighted sum |

**Methods:** `summary()`, `print_summary()`, `to_dataframe()`

### Comparison2x2

Individual 2x2 DiD comparison (used in BaconDecompositionResults).

| Attribute | Type | Description |
|-----------|------|-------------|
| `treated_group` | `Any` | Timing group used as treated |
| `control_group` | `Any` | Timing group used as control |
| `comparison_type` | `str` | "treated_vs_never", "earlier_vs_later", or "later_vs_earlier" |
| `estimate` | `float` | 2x2 DiD estimate |
| `weight` | `float` | Weight in TWFE average |
| `n_treated` | `int` | Number of treated observations |
| `n_control` | `int` | Number of control observations |
| `time_window` | `tuple[float, float]` | (start, end) time window |

### Common Results Pattern for Staggered Estimators

ImputationDiDResults, TwoStageDiDResults, StackedDiDResults, and EfficientDiDResults share a similar structure:

| Attribute | Type | Description |
|-----------|------|-------------|
| `overall_att` | `float` | Overall ATT |
| `overall_se` | `float` | SE of overall ATT |
| `overall_t_stat` | `float` | T-statistic |
| `overall_p_value` | `float` | P-value |
| `overall_conf_int` | `tuple[float, float]` | CI |
| `event_study_effects` | `dict[int, dict]` | Event study effects (if the DEPRECATED fit-time aggregate included event_study; ALWAYS populated for StackedDiD since 3.9 - row M-024; post-fit `aggregate('event_study')` returns a CONTAINER, it does not populate this field) |
| `group_effects` | `dict` | Group-level effects (if the deprecated fit-time aggregate included group; post-fit `aggregate('group')` returns an AggregationResult container) |
| `groups` | `list` | Treatment cohorts |
| `time_periods` | `list` | All time periods |
| `n_obs` | `int` | Number of observations |
| `n_treated_units` | `int` | Number of treated units |
| `n_control_units` | `int` | Number of control units |

Each event study effect dict contains: `effect`, `se`, `t_stat`, `p_value`, `conf_int`, `n_obs` (or `n_groups`).

**Variance metadata** (`ImputationDiDResults` and `EfficientDiDResults` carry these; other staggered Results may surface a subset):

| Attribute | Type | Description |
|-----------|------|-------------|
| `vcov_type` | `str` | Variance estimator family (`"hc1"` for IF-based estimators; permanently narrow on `ImputationDiD` / `CallawaySantAnna` / `TripleDifference` / `EfficientDiD` per IF-vs-sandwich taxonomy) |
| `cluster_name` | `str | None` | Effective cluster column name (e.g. `"unit"` for default `cluster=None` on `ImputationDiD`; `None` for default on `EfficientDiD` since its per-unit EIF SE is HC1-style not CR1-at-unit); `None` under survey designs (the survey block already names PSU/strata) |
| `n_clusters` | `int | None` | Number of effective clusters; `None` under survey designs and `None` under EfficientDiD's default unclustered fit |
| `df_convention` | `str | None` | On the knob-carrying containers (`StackedDiDResults`, `ImputationDiDResults`, `WooldridgeDiDResults`, `LPDiDResults`, `SunAbrahamResults`): the configured df convention (3.9 / M-127); `StackedDiDResults` additionally carries `inference_df` (the overall-ATT df actually used) |

**Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `to_dict()` (flat dict of headline aliases + `vcov_type` + conditional `cluster_name`/`n_clusters`/`n_bootstrap`/`inference_method`); `aggregate(type, weights=None, *, balance_e=None)` on the shipped post-fit adopters (`StackedDiDResults` views since 3.9/M-024; `EfficientDiDResults` recomputes from retained EIFs since 3.9/M-023; `ImputationDiDResults` and `TwoStageDiDResults` recompute from their PANEL-BACKED kits since 3.9/M-021/M-022 - on bootstrapped fits their recompute levels raise while `'simple'` relays the stored bootstrap inference with a NaN df column (the per-level rule, converged with M-027); `ContinuousDiDResults` is MIXED since 3.9/M-025 - `'simple'`/`'dose'` are views over stored fields that work on any fit incl. bootstrapped, `'event_study'` recomputes from a pruned per-cell IF kit and raises on bootstrapped fits; `HeterogeneousAdoptionDiDResults` and `HeterogeneousAdoptionDiDEventStudyResults` are pure views since 3.9/M-027 - `'simple'` on the overall class, `'event_study'` on the event-study class, no kit, work on pickles from any release)

### ContinuousDiDResults

| Attribute | Type | Description |
|-----------|------|-------------|
| `dose_response_att` | `DoseResponseCurve` | Dose-response curve for ATT |
| `dose_response_acrt` | `DoseResponseCurve` | Dose-response curve for ACRT |
| `overall_att` | `float` | Overall ATT |
| `overall_att_se` | `float` | SE of overall ATT |
| `overall_att_t_stat` | `float` | T-statistic for ATT |
| `overall_att_p_value` | `float` | P-value for ATT |
| `overall_att_conf_int` | `tuple[float, float]` | CI for ATT |
| `overall_acrt` | `float` | Overall ACRT |
| `overall_acrt_se` | `float` | SE of overall ACRT |
| `overall_acrt_t_stat` | `float` | T-statistic for ACRT |
| `overall_acrt_p_value` | `float` | P-value for ACRT |
| `overall_acrt_conf_int` | `tuple[float, float]` | CI for ACRT |
| `group_time_effects` | `dict[tuple, dict]` | Group-time level effects |
| `dose_grid` | `np.ndarray` | Evaluation grid for dose-response |
| `groups` | `list` | Treatment cohorts |
| `time_periods` | `list` | All time periods |
| `n_obs` | `int` | Number of observations |
| `n_treated_units` | `int` | Treated units |
| `n_control_units` | `int` | Control units |
| `event_study_effects` | `dict[int, dict] or None` | Event study effects (populated only by the DEPRECATED fit-time `aggregate="eventstudy"`; the post-fit successor `results.aggregate('event_study')` returns a unified `EventStudyResults` container instead of populating this field) |

**DoseResponseCurve** sub-dataclass:

| Attribute | Type | Description |
|-----------|------|-------------|
| `dose_grid` | `np.ndarray` | Dose values |
| `effects` | `np.ndarray` | Estimated effects at each dose |
| `se` | `np.ndarray` | Standard errors |
| `conf_int_lower` | `np.ndarray` | Lower CI bound |
| `conf_int_upper` | `np.ndarray` | Upper CI bound |
| `target` | `str` | `"att"` or `"acrt"` |

**Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `aggregate(type, ...)` (M-025 MIXED adopter: `'simple'` = 2-row att+acrt view, `'dose'` = 2N-row target-discriminated curve view - both work on any fit; `'event_study'` = kit recompute, analytical fits only)

### HeterogeneousAdoptionDiDResults

Single-period results container for `HeterogeneousAdoptionDiD`. The table below enumerates every public dataclass field; a regression test in `tests/test_guides.py` (`test_llms_full_had_results_class_field_lists_match_real_dataclass`) compares this list against the real `dataclasses.fields()` of the result class.

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | Point estimate of the WAS parameter on the β-scale |
| `se` | `float` | Standard error on the β-scale |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | Confidence interval |
| `alpha` | `float` | CI level used at fit time |
| `design` | `str` | Resolved design: `"continuous_at_zero"`, `"continuous_near_d_lower"`, or `"mass_point"` |
| `target_parameter` | `str` | `"WAS"` (Design 1') or `"WAS_d_lower"` (Design 1 / mass-point) |
| `d_lower` | `float` | Support infimum (`0.0` on Design 1', `min(d)` otherwise) |
| `dose_mean` | `float` | `D_bar = (1/G) * sum(D_{g,2})` |
| `n_obs` | `int` | Units contributing to estimation |
| `n_treated` | `int` | Units with `D > d_lower` |
| `n_control` | `int` | Units at or below `d_lower` |
| `n_mass_point` | `int | None` | Mass-point design only: units exactly at `d_lower`; `None` on continuous designs |
| `n_above_d_lower` | `int | None` | Mass-point design only: units strictly above `d_lower`; `None` on continuous designs |
| `inference_method` | `str` | `"analytical_nonparametric"` or `"analytical_2sls"` |
| `vcov_type` | `str | None` | Mass-point only: `"classical"`, `"hc1"`, or `"cr1"` |
| `cluster_name` | `str | None` | Cluster column name when CR1 cluster-robust SE is requested; `None` otherwise |
| `survey_metadata` | `SurveyMetadata | None` | Repo-standard survey metadata when `survey_design=` is supplied |
| `bandwidth_diagnostics` | `BandwidthResult | None` | MSE-DPI selector output (continuous designs); `None` on `mass_point` |
| `bias_corrected_fit` | `BiasCorrectedFit | None` | Phase 1c bias-corrected local-linear fit object (continuous designs); `None` on `mass_point` |
| `variance_formula` | `str | None` | HAD-specific SE label on weighted fits, populated on BOTH continuous and mass-point designs: `"survey_binder_tsl"` (continuous, Binder 1983 TSL on the `survey_design=` path) or `"survey_binder_tsl_2sls"` (mass-point, Binder 1983 TSL on the `survey_design=` path; requires `vcov_type="hc1"` — the mass-point survey path rejects `vcov_type="classical"`, and `cluster=` + `survey_design=` is rejected, so PSU clustering is expressed via `SurveyDesign(weights='<weight_col>', psu='<cluster_col>')`). `None` on unweighted fits |
| `effective_dose_mean` | `float | None` | Weighted denominator used by the β̂-scale rescaling, populated on weighted fits across all designs: weighted `mean(d)` (`continuous_at_zero`), weighted `mean(d − d_lower)` (`continuous_near_d_lower`), or weighted Wald-IV dose gap `mean(d | Z=1, w) − mean(d | Z=0, w)` (`mass_point`). `None` on unweighted fits |

**Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `aggregate(type, ...)` (M-027 pure view: `'simple'` only - a one-row bit-exact relay whose `target` is the estimand label `WAS`/`WAS_d_lower`; no kit, works on pickles from any release; `'event_study'` needs a multi-period fit)

### HeterogeneousAdoptionDiDEventStudyResults

Per-horizon event-study results container for `HeterogeneousAdoptionDiD`'s event-study mode (multi-period panels; M-027). The anchor horizon `e = -1` is excluded by construction. The table below enumerates every public dataclass field; a regression test (`test_llms_full_had_results_class_field_lists_match_real_dataclass`) compares this list against the real `dataclasses.fields()`.

| Attribute | Type | Description |
|-----------|------|-------------|
| `event_times` | `np.ndarray` | Integer event-time labels `e = t - F`, sorted ascending |
| `att` | `np.ndarray` | Per-horizon WAS point estimates |
| `se` | `np.ndarray` | Per-horizon standard errors |
| `t_stat` | `np.ndarray` | Per-horizon t-statistics |
| `p_value` | `np.ndarray` | Per-horizon p-values |
| `conf_int_low` | `np.ndarray` | Pointwise CI lower bounds |
| `conf_int_high` | `np.ndarray` | Pointwise CI upper bounds |
| `n_obs_per_horizon` | `np.ndarray` | Per-horizon contributing-unit counts |
| `alpha` | `float` | CI level used at fit time |
| `design` | `str` | Shared across horizons (paper Appendix B.2 invariant) |
| `target_parameter` | `str` | Same convention as the single-period result |
| `d_lower` | `float` | Support infimum, shared across horizons |
| `dose_mean` | `float` | `D_bar` on the fit sample |
| `F` | `object` | First-treatment period label |
| `n_units` | `int` | Unique units contributing to the fit (post last-cohort filter) |
| `inference_method` | `str` | `"analytical_nonparametric"` or `"analytical_2sls"` |
| `vcov_type` | `str | None` | Mass-point: `"classical"`, `"hc1"`, or `"cr1"` (with `cluster=`). Continuous: `None`, or `"cr1"` with `cluster=` |
| `cluster_name` | `str | None` | Cluster column name when CR1 is requested; `None` otherwise |
| `survey_metadata` | `SurveyMetadata | None` | Populated on weighted fits |
| `bandwidth_diagnostics` | `list[BandwidthResult | None] | None` | Per-horizon MSE-DPI selector output (continuous designs); `None` on `mass_point`; entries can be `None` on degenerate horizons |
| `bias_corrected_fit` | `list[BiasCorrectedFit | None] | None` | Per-horizon Phase 1c bias-corrected local-linear fit objects; `None` on `mass_point`; entries can be `None` on degenerate horizons |
| `filter_info` | `dict | None` | Staggered last-cohort auto-filter metadata (`F_last`, `n_kept`, `n_dropped`, `dropped_cohorts`); `None` when no filter applied |
| `variance_formula` | `str | None` | HAD-specific SE label applied UNIFORMLY across all horizons, populated on BOTH continuous and mass-point designs: `"survey_binder_tsl"` (continuous, Binder 1983 TSL on the `survey_design=` path) or `"survey_binder_tsl_2sls"` (mass-point, Binder 1983 TSL on the `survey_design=` path; requires `vcov_type="hc1"` — the mass-point survey path rejects `vcov_type="classical"`, and `cluster=` + `survey_design=` is rejected, so PSU clustering is expressed via `SurveyDesign(weights='<weight_col>', psu='<cluster_col>')`). `None` on unweighted fits |
| `effective_dose_mean` | `float | None` | Weighted denominator used by the β̂-scale rescaling, populated on weighted fits across all designs: weighted `sum(w·d)/sum(w)` (`continuous_at_zero`), weighted `sum(w·(d − d_lower))/sum(w)` (`continuous_near_d_lower`), or weighted Wald-IV dose gap (`mass_point`). Scalar (not per-horizon) because the β̂-scale denominator is computed once on the fit sample. `None` on unweighted fits |
| `cband_low` | `np.ndarray | None` | Simultaneous (sup-t) band lower bounds; `None` when `cband=False` or on unweighted, unclustered fits (a clustered fit produces the band even when unweighted) |
| `cband_high` | `np.ndarray | None` | Simultaneous (sup-t) band upper bounds |
| `cband_crit_value` | `float | None` | Sup-t critical value used for the simultaneous band |
| `cband_method` | `str | None` | `"multiplier_bootstrap"` (weighted/survey band) or `"cluster_multiplier_bootstrap"` (clustered band) when populated |
| `cband_n_bootstrap` | `int | None` | Bootstrap iterations used for the band |

**Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `aggregate(type, ...)` (M-027 pure view: `'event_study'` only - the unified `EventStudyResults` container via the `_from_had` passthrough, cband fields included; no kit; `'simple'` needs a two-period fit - no overall WAS is stored here)

### ChangesInChangesResults

Results container shared by `ChangesInChanges` and `QDiD` (the `method` field is `"cic"` or `"qdid"`; the pre-3.9 name `estimator` still reads it with a FutureWarning until 4.0; `QDiDResults` is an alias of this class). Flat-native headline fields `att`, `se`, `t_stat`, `p_value`, `conf_int`; `quantile_effects` is a DataFrame with columns `quantile`, `qte`, `se`, `t_stat`, `p_value`, `conf_low`, `conf_high`. `q_lower`/`q_upper` bound the point-identified interior quantile range for unconditional CiC fits (NaN for QDiD and for covariate fits); `sup_t_crit` is the qte sup-t critical value backing `uniform_bands()` (fixed 95% level). Also carries `n_obs`, `cell_sizes`, `n_bootstrap`, `n_bootstrap_valid`, `panel`, `quantiles`, `alpha`, and `covariates` (the covariate columns of a conditional fit, else None). Methods: `summary()`, `print_summary()`, `to_dict()`, `to_dataframe(level="quantiles"|"att")`, `uniform_bands()`. All inference flows through `safe_inference`/`safe_inference_batch` (joint-NaN contract; `n_bootstrap=0` yields NaN inference everywhere).

### TROPResults

| Attribute | Type | Description |
|-----------|------|-------------|
| `att` | `float` | ATT estimate |
| `se` | `float` | Bootstrap standard error |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | CI |
| `n_obs` | `int` | Number of observations |
| `n_treated` | `int` | Number of treated units |
| `n_control` | `int` | Number of control units |
| `n_treated_obs` | `int` | Number of treated unit-time observations |
| `lambda_time` | `float` | Selected time decay parameter |
| `lambda_unit` | `float` | Selected unit decay parameter |
| `lambda_nn` | `float` | Selected nuclear norm parameter |
| `n_bootstrap` | `int` | Number of bootstrap replications |

**Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`

## Diagnostics

### Placebo Tests

```python
from diff_diff import (
    run_placebo_test,
    placebo_timing_test,
    placebo_group_test,
    permutation_test,
    leave_one_out_test,
    run_all_placebo_tests,
)

# Unified interface
results = run_placebo_test(
    data, outcome='y', treatment='treated', time='period',
    test_type='fake_timing',           # "fake_timing", "fake_group", "permutation", "leave_one_out"
    fake_treatment_period=1,           # For fake_timing
    post_periods=[3, 4, 5],
)

# Run all tests at once
all_results = run_all_placebo_tests(
    data, outcome='y', treatment='treated', time='period', unit='unit_id',
    pre_periods=[0, 1, 2], post_periods=[3, 4, 5],
    n_permutations=500, seed=42,
)
```

**Individual test functions:**

```python
# Fake timing test
placebo_timing_test(data, outcome, treatment, time,
                    fake_treatment_period, post_periods=None, alpha=0.05)

# Fake group test
placebo_group_test(data, outcome, time, unit,
                   fake_treated_units, post_periods=None, alpha=0.05)

# Permutation test
permutation_test(data, outcome, treatment, post, unit,
                 n_permutations=1000, alpha=0.05, seed=None)

# Leave-one-out test
leave_one_out_test(data, outcome, treatment, post, unit, alpha=0.05)
```

`permutation_test` / `leave_one_out_test` accept `time=` as a deprecated alias for `post=` (FutureWarning; removed in 4.0). The `run_placebo_test` / `run_all_placebo_tests` wrappers keep their `time` parameter.

All return `PlaceboTestResults` with attributes: `test_type`, `placebo_effect`, `se`, `t_stat`, `p_value`, `conf_int`, `n_obs`, `is_significant`.

### Parallel Trends Testing

```python
from diff_diff import check_parallel_trends, check_parallel_trends_robust, equivalence_test_trends

# Simple trend comparison
result = check_parallel_trends(
    data, outcome='y', time='period', treatment_group='treated',
    pre_periods=[0, 1, 2],
)

# Distributional comparison (Wasserstein distance + permutation inference)
result = check_parallel_trends_robust(
    data, outcome='y', time='period', treatment_group='treated',
    unit='unit_id', pre_periods=[0, 1, 2],
    n_permutations=1000, seed=42,
)

# TOST equivalence test
result = equivalence_test_trends(
    data, outcome='y', time='period', treatment_group='treated',
    unit='unit_id', pre_periods=[0, 1, 2],
    equivalence_margin=0.5,
)
```

### Wild Cluster Bootstrap

```python
from diff_diff import wild_bootstrap_se, WildBootstrapResults

# Directly via estimator
did = DifferenceInDifferences(inference="wild_bootstrap", n_bootstrap=999,
                              bootstrap_weights="webb", cluster="state")
results = did.fit(data, outcome='y', treatment='treated', post='post')
```

## HAD Pretests

Diagnostic pretests for the `HeterogeneousAdoptionDiD` identifying assumptions (de Chaisemartin, Ciccia, D'Haultfœuille & Knau 2026). The composite workflow `did_had_pretest_workflow` is the recommended entry point — call it before reporting WAS as causal. The workflow follows paper Section 4.2's three-step battery: **step 1** is the QUG support-infimum test (decides whether Design 1' or Design 1 applies); **step 2** is the Assumption 7 pre-trends test (joint Stute on the event-study path; explicitly NOT covered on the overall path because a single-pre-period panel cannot support the joint variant); **step 3** is the Assumption 8 linearity test (`stute_test` or `yatchew_hr_test`). On a two-period panel the workflow runs the overall battery (steps 1 + 3 only) and the returned `verdict` flags the Assumption 7 gap; a multi-period panel selects the event-study battery (M-139 - the mode is panel-inferred, like `fit()`), which closes that gap **when at least one earlier placebo pre-period beyond the base `F-1` exists**. With only the base `F-1` pre-period available (minimal 3-period event-study, or 4-period under `trends_lin=True` where the consumed `F-2` placebo is dropped), the workflow still sets `pretrends_joint=None`, `all_pass=False`, and appends `joint pre-trends skipped (no earlier pre-period)` to the verdict — step 2 stays uncovered.

```python
from diff_diff import (
    did_had_pretest_workflow,
    qug_test, stute_test, yatchew_hr_test,
    stute_joint_pretest, joint_pretrends_test, joint_homogeneity_test,
)

# Composite workflow - the battery is selected from the panel shape (M-139):
#   two-period panel   -> steps 1 + 3 (QUG + Assumption 8 linearity)
#                         step 2 (Assumption 7 pre-trends) NOT covered;
#                         verdict explicitly flags this gap.
#   multi-period panel -> steps 1 + 2 + 3 (QUG + joint Stute pre-trends +
#                         joint homogeneity-linearity Stute).
report = did_had_pretest_workflow(
    data, outcome='y', unit='unit', time='t',
    dose='d', first_treat='first_treat',
    survey_design=None)   # SurveyDesign for survey-aware pretests (Phase 4.5 C)
print(report.summary())
print(report.all_pass, report.verdict)
```

Individual tests:

- `qug_test(d)` — paper Theorem 4 support-infimum test (`H_0: d_lower = 0`; the QUG decides whether Design 1' or Design 1 applies in step 1 of the workflow). Extreme order statistics, Exp(1)/Exp(1) limit law. The QUG itself does NOT test Assumption 5 (which is the Design 1 sign-identification condition and is not testable via pre-trends per registry). **Permanently rejects** non-`None` `survey_design=` (`NotImplementedError`) per Phase 4.5 C0 deferral (the deprecated `survey=`/`weights=` aliases were removed in 3.7.x and now raise `TypeError` on all 7 pretest helpers) — extreme-value functionals are not smooth in the empirical CDF, so standard survey machinery does not yield a calibrated test.
- `stute_test(d, dy)` — Assumption 8 linearity of `E[ΔY|D]` (paper Section 4.2 step 3) via Stute Cramér-von Mises functional with Mammen wild bootstrap. Survey-aware via PSU-level Mammen multiplier bootstrap.
- `yatchew_hr_test(d, dy, *, null="linearity")` — Assumption 8 linearity of `E[ΔY|D]` (alternative test for step 3) via Yatchew (1997) heteroskedasticity-robust variance-ratio test. The `null="mean_independence"` mode (R `YatchewTest::yatchew_test(order=0)`) is also exposed for placebo-style mean-independence testing. Survey-aware via closed-form weighted variance components (no bootstrap).
- `stute_joint_pretest(residuals_dict, d)` — joint Cramér-von Mises across K horizons with shared-η Mammen wild bootstrap (Delgado-Manteiga 2001 / Hlávka-Hušková 2020). Residuals-in core; the two data-in wrappers below construct residuals for the two paper-spelled nulls.
- `joint_pretrends_test(...)` — Assumption 7 joint pre-trends on K pre-periods (paper Section 4.2 step 2 closure on the event-study path).
- `joint_homogeneity_test(...)` — joint linearity-and-homogeneity on K post-periods (event-study step 3 alternative).

The QUG-under-survey deferral is permanent; the linearity-family pretests support `survey_design=` (pweight, PSU, FPC, strata — the stratified Stute calibration shipped with the Phase 4.5 C strata extension) per Phase 4.5 C. Replicate-weight designs and `lonely_psu='adjust'` under singleton strata remain deferred to follow-up PRs.

## Honest DiD Sensitivity Analysis

Rambachan & Roth (2023) robust inference allowing bounded parallel trends violations.

### Delta Restriction Classes

```python
from diff_diff import DeltaSD, DeltaRM, DeltaSDRM

# Smoothness: bounds on second differences
delta_sd = DeltaSD(M=0.5)

# Relative magnitudes: post violations <= Mbar * max pre violation
delta_rm = DeltaRM(Mbar=1.0)

# Combined restriction
delta_sdrm = DeltaSDRM(M=0.5, Mbar=1.0)
```

### HonestDiD Class

```python
from diff_diff import HonestDiD

honest = HonestDiD(
    method="relative_magnitude",     # "smoothness", "relative_magnitude", or "combined"
    M=1.0,                           # Restriction parameter
    alpha=0.05,
    l_vec=None,                      # Weighting vector (None = uniform)
)

# Fit to event study results
bounds = honest.fit(event_study_results)
print(bounds.summary())

# Sensitivity analysis over M grid
sensitivity = honest.sensitivity_analysis(
    event_study_results,
    M_grid=[0, 0.5, 1.0, 1.5, 2.0],
)
sensitivity.plot()
```

### Convenience Functions

```python
from diff_diff import compute_honest_did, sensitivity_plot

bounds = compute_honest_did(results, method="relative_magnitude", M=1.0, alpha=0.05)
sensitivity_plot(results, method="relative_magnitude", M_grid=[0, 0.5, 1, 1.5, 2])
```

### HonestDiDResults

| Attribute | Type | Description |
|-----------|------|-------------|
| `lb` | `float` | Lower bound of identified set |
| `ub` | `float` | Upper bound of identified set |
| `ci_lb` | `float` | Lower bound of robust CI |
| `ci_ub` | `float` | Upper bound of robust CI |
| `M` | `float` | Restriction parameter value |
| `method` | `str` | Restriction type |
| `original_estimate` | `float` | Original point estimate |
| `original_se` | `float` | Original SE |
| `ci_method` | `str` | "FLCI" or "C-LF" |
| `event_study_bounds` | `dict` | Per-period bounds (optional) |

**Properties:** `is_significant` (CI excludes zero)

## Power Analysis

```python
from diff_diff import PowerAnalysis, compute_mde, compute_power, compute_sample_size, simulate_power

# Class-based interface
pa = PowerAnalysis(alpha=0.05, power=0.80, alternative='two-sided')
mde_result = pa.mde(n_treated=50, n_control=50, sigma=1.0)
sample_result = pa.sample_size(effect_size=0.5, sigma=1.0)
power_result = pa.power(effect_size=0.5, n_treated=50, n_control=50, sigma=1.0)

# Convenience functions
mde_result = compute_mde(n_treated=50, n_control=50, sigma=1.0)
power_result = compute_power(effect_size=0.5, n_treated=50, n_control=50, sigma=1.0)
sample_result = compute_sample_size(effect_size=0.5, sigma=1.0)

# Simulation-based power
sim_result = simulate_power(
    effect_sizes=[0.1, 0.5, 1.0, 2.0],
    n_simulations=500, seed=42,
)
```

## Pre-Trends Power Analysis

```python
from diff_diff import PreTrendsPower, compute_pretrends_power, compute_mdv

# Class-based
ptp = PreTrendsPower()
results = ptp.compute(event_study_results, M_grid=[0, 0.5, 1.0, 2.0])

# Convenience functions
results = compute_pretrends_power(event_study_results, M_grid=[0, 0.5, 1.0, 2.0])
mdv = compute_mdv(event_study_results, target_power=0.80)
```

## Visualization

All plotting functions return a matplotlib `Figure` object.

### plot_event_study

```python
from diff_diff import plot_event_study

plot_event_study(
    results,                           # MultiPeriodDiDResults, CS, SA, BJS, TwoStageDiD, StackedDiD, or DataFrame
    effects=None,                      # Manual dict of effects (alternative to results)
    se=None,                           # Manual dict of SEs
    periods=None,
    reference_period=None,
    pre_periods=None,
    post_periods=None,
    alpha=0.05,
    figsize=(10, 6),
    title="Event Study",
    xlabel="Period Relative to Treatment",
    ylabel="Treatment Effect",
    color="#2563eb",
    show_zero_line=True,
    show_reference_line=True,
    shade_pre=True,
    ax=None,
    show=True,
    use_cband=True,                    # Use simultaneous confidence bands if available
)
```

### plot_group_effects

```python
from diff_diff import plot_group_effects

plot_group_effects(
    results,                           # CallawaySantAnnaResults
    groups=None,
    figsize=(10, 6),
    title="Treatment Effects by Cohort",
    alpha=0.05,
    show=True,
    ax=None,
)
```

### plot_sensitivity

```python
from diff_diff import plot_sensitivity

plot_sensitivity(
    sensitivity_results,               # SensitivityResults from HonestDiD
    show_bounds=True,
    show_ci=True,
    breakdown_line=True,
    figsize=(10, 6),
    title="Honest DiD Sensitivity Analysis",
    ax=None,
    show=True,
)
```

### plot_honest_event_study

```python
from diff_diff import plot_honest_event_study

plot_honest_event_study(
    honest_results,                    # HonestDiDResults with event_study_bounds
    periods=None,
    reference_period=None,
    figsize=(10, 6),
    title="Event Study with Honest Confidence Intervals",
    ax=None,
    show=True,
)
```

### plot_bacon

```python
from diff_diff import plot_bacon

plot_bacon(
    results,                           # BaconDecompositionResults
    plot_type="scatter",               # "scatter" or "bar"
    figsize=(10, 6),
    show_weighted_avg=True,
    show_twfe_line=True,
    ax=None,
    show=True,
)
```

### plot_power_curve

```python
from diff_diff import plot_power_curve

plot_power_curve(
    results=None,                      # PowerResults, SimulationPowerResults, or DataFrame
    effect_sizes=None,
    powers=None,
    mde=None,
    target_power=0.80,
    plot_type="effect",                # "effect" or "sample_size"
    figsize=(10, 6),
    show_mde_line=True,
    show_target_line=True,
    ax=None,
    show=True,
)
```

### plot_pretrends_power

```python
from diff_diff import plot_pretrends_power

plot_pretrends_power(
    results=None,                      # PreTrendsPowerResults or PreTrendsPowerCurve
    M_values=None,
    powers=None,
    mdv=None,
    target_power=0.80,
    figsize=(10, 6),
    ax=None,
    show=True,
)
```

## Data Preparation Utilities

### Data Manipulation

```python
from diff_diff import (
    make_treatment_indicator,
    make_post_indicator,
    wide_to_long,
    balance_panel,
    validate_did_data,
    summarize_did_data,
    create_event_time,
    aggregate_to_cohorts,
    rank_control_units,
)

# Create binary treatment indicator
df = make_treatment_indicator(data, column='group', treated_values='A', new_column='treated')
df = make_treatment_indicator(data, column='size', threshold=75, new_column='treated')

# Create binary post indicator
df = make_post_indicator(data, time_column='year', treatment_start=2020, new_column='post')
df = make_post_indicator(data, time_column='year', post_periods=[2020, 2021])

# Reshape wide to long
long_df = wide_to_long(data, value_columns=['y2018', 'y2019', 'y2020'],
                       id_column='unit', time_name='year', value_name='outcome')

# Balance panel (keep only units observed in all periods)
balanced_df = balance_panel(data, unit='unit', time='period')

# Validate DiD data
validation = validate_did_data(data, outcome='y', treatment='treated',
                               time='period', unit='unit')

# Summarize DiD data
summary = summarize_did_data(data, outcome='y', treatment='treated',
                             time='period', unit='unit')

# Create event time column
df = create_event_time(data, time='period', first_treat='first_treat', new_column='event_time')

# Aggregate to cohort level
cohort_df = aggregate_to_cohorts(data, outcome='y', unit='unit', time='period',
                                 first_treat='first_treat')

# Rank control units by similarity to treated
ranking = rank_control_units(data, outcome='y', unit='unit', time='period',
                             treatment='treated')
```

### Data Generation

```python
from diff_diff import (
    generate_did_data,
    generate_staggered_data,
    generate_panel_data,
    generate_event_study_data,
    generate_factor_data,
    generate_ddd_data,
    generate_continuous_did_data,
    generate_synthetic_control_data,
)

# Basic 2x2 DiD data
data = generate_did_data(n_units=100, n_periods=4, treatment_effect=5.0,
                         treatment_fraction=0.5, treatment_period=2, seed=42)

# Staggered adoption data
data = generate_staggered_data(n_units=100, n_periods=10,
                               treatment_effect=2.0, dynamic_effects=True,
                               never_treated_frac=0.3, seed=42)

# Panel data with optional trend violations
data = generate_panel_data(n_units=100, n_periods=8, treatment_period=4,
                           parallel_trends=True, seed=42)

# Event study data
data = generate_event_study_data(n_units=300, n_pre=5, n_post=5,
                                 treatment_effect=5.0, seed=42)

# Factor model data (for TROP)
data = generate_factor_data(n_units=50, n_pre=10, n_post=5,
                            n_treated=10, n_factors=2, seed=42)

# Triple difference data
data = generate_ddd_data(n_per_cell=100, treatment_effect=2.0, seed=42)

# Continuous dose data
data = generate_continuous_did_data(n_units=500, n_periods=4,
                                    att_function="linear", att_slope=2.0, seed=42)

# Single-treated-unit panel for synthetic control (noiseless treated path in the
# donor convex hull; ramping or constant effect). Use with SyntheticControl.
data = generate_synthetic_control_data(n_donors=20, n_pre=60, n_post=5,
                                       effect_type="ramp", seed=0)
```

## Built-in Datasets

```python
from diff_diff import load_card_krueger, load_castle_doctrine, load_divorce_laws, load_mpdta
from diff_diff import load_prop99, load_walmart
from diff_diff import load_dataset, list_datasets, clear_cache

# List available datasets
for name, desc in list_datasets().items():
    print(f"{name}: {desc}")

# Load by name
data = load_dataset("card_krueger")

# Named loaders
ck = load_card_krueger()          # Card & Krueger (1994) minimum wage
castle = load_castle_doctrine()    # Castle Doctrine / Stand Your Ground laws
divorce = load_divorce_laws()      # Synthetic-only unilateral divorce-law fallback
mpdta = load_mpdta()              # County teen-employment panel from R did package
prop99 = load_prop99()            # California Prop 99 smoking (single treated unit)
walmart = load_walmart()          # Walmart entry county panel (staggered, 1,277 counties)

# Force re-download
data = load_card_krueger(force_download=True)

# Clear local cache
clear_cache()
```

## Survey Support

Most estimators accept an optional `survey_design` parameter in `fit()` (`SyntheticControl` rejects it as not yet supported); depth of support varies by estimator - see the compatibility matrix in `docs/choosing_estimator.rst` (Survey Design Support). Pass a `SurveyDesign` object to get design-based variance estimation.

```python
from diff_diff import SurveyDesign, CallawaySantAnna

# Create survey design with strata, PSU, FPC
sd = SurveyDesign(
    weights='weight',           # Sampling weight column
    strata='stratum',           # Stratification variable
    psu='psu',                  # Primary Sampling Unit (cluster)
    fpc='fpc',                  # Finite Population Correction
    weight_type='pweight',      # "pweight", "fweight", or "aweight"
    nest=False,                 # PSU IDs nested within strata?
    lonely_psu='remove',        # "remove", "certainty", or "adjust"
)

# Fit with survey design
cs = CallawaySantAnna(estimation_method='dr')
results = cs.fit(data, outcome='y', unit='id', time='t',
                 first_treat='ft', survey_design=sd)

# Survey metadata on results
print(results.survey_metadata.design_effect)  # Design effect
print(results.survey_metadata.effective_n)  # Effective sample size
print(results.survey_metadata.df_survey)    # Survey degrees of freedom
```

**Replicate weight designs:**

```python
sd = SurveyDesign(
    weights='weight',
    replicate_weights=['rep_0', 'rep_1', ..., 'rep_79'],  # Column names
    replicate_method='SDR',     # "BRR", "Fay", "JK1", "JKn", or "SDR"
)
```

**Subpopulation analysis:**

```python
sd_female, data_female = sd.subpopulation(data, mask=lambda df: df['sex'] == 'F')
```

**Key features:**
- Taylor Series Linearization (TSL) variance with strata + PSU + FPC
- Replicate weight variance: BRR, Fay's BRR, JK1, JKn, SDR (13 of 23 estimators, including dCDH)
- Survey-aware bootstrap: multiplier at PSU (Hall-Mammen wild; dCDH, staggered) or Rao-Wu rescaled (SunAbraham, SyntheticDiD, TROP). SyntheticDiD bootstrap composes Rao-Wu rescaled per-draw weights with the weighted Frank-Wolfe variant of `_sc_weight_fw` (PR #355): each draw solves `min ||A·diag(rw)·ω - b||² + ζ²·Σ rw_i ω_i²` and composes `ω_eff = rw·ω/Σ(rw·ω)` for the SDID estimator. Pweight-only fits use constant `rw = w_control`; full designs use Rao-Wu. SDID's placebo (stratified permutation + weighted FW) and jackknife (PSU-level LOO with stratum aggregation, Rust & Rao 1996) paths also support pweight-only and full strata/PSU/FPC designs
- DEFF diagnostics, subpopulation analysis, weight trimming (`trim_weights`)
- Repeated cross-sections: `CallawaySantAnna(panel=False)`
- Compatibility matrix: see `docs/choosing_estimator.rst` Survey Design Support section

No R or Python package offers design-based variance estimation for modern heterogeneity-robust DiD estimators.

## Linear Algebra Helpers

```python
from diff_diff import LinearRegression, InferenceResult

# Low-level regression helper
reg = LinearRegression(
    include_intercept=True,
    cluster_ids=cluster_array,
)
reg.fit(X, y)
inference = reg.get_inference(coef_index)  # -> InferenceResult
```

### InferenceResult

| Attribute | Type | Description |
|-----------|------|-------------|
| `coefficient` | `float` | Point estimate |
| `se` | `float` | Standard error |
| `t_stat` | `float` | T-statistic |
| `p_value` | `float` | P-value |
| `conf_int` | `tuple[float, float]` | Confidence interval |

### Conley Spatial HAC Standard Errors

Conley (1999) spatial heteroskedasticity-and-autocorrelation-consistent standard
errors. Use when residuals are spatially (and optionally temporally) correlated
(geo experiments, regional shocks, common-supplier effects). Two operating
modes:

- **Cross-sectional:** Direct `compute_robust_vcov` / `LinearRegression` on
  a single-period design.
- **Panel block-decomposed** (matches R `conleyreg` with `lag_cutoff > 0`):
  `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` with
  `vcov_type="conley"` and `conley_lag_cutoff=<int>`. The sandwich sums
  within-period spatial pairs plus within-unit Bartlett serial pairs
  (excluding lag=0 to avoid double-counting); NOT a multiplicative
  product kernel.

`DifferenceInDifferences(vcov_type="conley").fit(..., unit="<col>")` is
supported (Wave A #118). `unit` is a fit-time kwarg (NOT on `__init__`;
unused unless Conley is set; not part of `get_params()` / `set_params()`)
mirroring `MultiPeriodDiD.fit(unit=...)` / `TwoWayFixedEffects.fit(unit=...)`.

```python
import numpy as np
from diff_diff.linalg import LinearRegression
from diff_diff import DifferenceInDifferences, TwoWayFixedEffects

# Cross-sectional design: 1 row per unit, n × 2 lat/lon coords.
reg = LinearRegression(
    vcov_type="conley",
    include_intercept=True,
    conley_coords=coords,                # n × 2 array of (lat, lon) in degrees
    conley_cutoff_km=200.0,              # required; no default
    conley_metric="haversine",           # or "euclidean", or callable
    conley_kernel="bartlett",            # or "uniform"
).fit(X, y)
se = np.sqrt(np.diag(reg.vcov_))

# Panel design: TWFE with within-unit Bartlett serial HAC.
# Static TWFE's `post` column is a binary post indicator (treatment * post
# interaction; `time=` is its deprecated alias through 3.9); only numeric
# encodings are supported on this surface. `conley_lag_cutoff=1` includes
# the cross-period pair under the Bartlett taper. For multi-period panels,
# use the TWFE event-study mode (`event_study=True, time=<calendar col>`;
# MultiPeriodDiD is deprecated).
res = TwoWayFixedEffects(
    vcov_type="conley",
    conley_coords=("lat", "lon"),        # column names on `data`
    conley_cutoff_km=500.0,
    conley_lag_cutoff=1,                 # within-unit Bartlett, lag 1 panel period
).fit(data, outcome="y", treatment="treated", post="post", unit="unit_id")

# Panel design: multi-period event study (TWFE event-study mode; the
# deprecated MultiPeriodDiD behaves identically via the shared core).
# The design builds period dummies (NOT a treated * time product), so
# `time` can be any orderable encoding — int years (2020, 2021, ...),
# YYYYMM (202012, 202101, ...), datetime64, pd.Period, strings.
# `_compute_conley_vcov` normalizes time to dense codes 0..T-1 internally,
# so `conley_lag_cutoff` always counts panel periods regardless of label.
mp_res = TwoWayFixedEffects(
    vcov_type="conley",
    conley_coords=("lat", "lon"),
    conley_cutoff_km=200.0,
    conley_lag_cutoff=2,                 # within-unit Bartlett up to 2 panel periods
).fit(data, outcome="y", treatment="treated", event_study=True,
      spec="pooled", time="period", post_periods=[2, 3], unit="unit_id")

# 2-period DiD on a panel: DiD.fit(unit="<col>") opts into the Conley
# panel block-decomposed sandwich; on a 2-period design the ATT/SE match
# the pooled event-study fit with post_periods=[1], reference_period=0 bit-exactly.
did_res = DifferenceInDifferences(
    vcov_type="conley",
    conley_coords=("lat", "lon"),
    conley_cutoff_km=200.0,
    conley_lag_cutoff=1,
).fit(data, outcome="y", treatment="treated", post="post", unit="unit_id")

# Combined spatial + cluster product kernel: pass cluster=<col> alongside
# Conley to apply K_total[i,j] = K_space(d_ij/h) · 1{c_i = c_j}. On the
# panel path the cluster must be constant within each unit across periods
# (e.g. an above-unit grouping like region); time-varying cluster raises
# ValueError. TWFE's default auto-cluster on the Conley path is silently
# dropped — users opt into the combined kernel explicitly.
combined_res = TwoWayFixedEffects(
    vcov_type="conley",
    cluster="region",                    # above-unit grouping; time-invariant within unit
    conley_coords=("lat", "lon"),
    conley_cutoff_km=500.0,
    conley_lag_cutoff=1,
).fit(data, outcome="y", treatment="treated", post="post", unit="unit_id")
```

**Note on `conley_lag_cutoff` semantics:** the lag is counted in **panel
periods** (number of distinct sorted values in the `time` column), NOT in
raw-label differences. Internally, time labels are normalized to dense codes
`0..T-1` via `np.unique(return_inverse=True)`. For example, MultiPeriodDiD
with `time` values `(2020, 2021, 2022)` and `(202012, 202101, 202102)` both
produce the same codes `(0, 1, 2)` and the same lag matrix. **TWFE caveat:**
the time-label normalization runs inside `_compute_conley_vcov`, but TWFE's
own design step `_treatment_post = treated * time` requires numeric `time`;
non-numeric labels (datetime64, pd.Period, strings) are TWFE-incompatible
end-to-end. Use MultiPeriodDiD if you need datetime/Period/string time
labels. **Deviation from R `conleyreg`:** R uses raw `time` values directly
in the lag computation, which silently mishandles non-dense encodings.
diff-diff is the more robust default; for bit-exact R parity, pass `time`
as a dense integer index.

**Variance estimator — cross-sectional:**

    Var̂(β) = (X'X)^{-1} · ( Σ_{i,j} K(d_ij / h) · X_i ε_i ε_j X_j' ) · (X'X)^{-1}

**Variance estimator — panel block-decomposed (matches R `conleyreg`):**

    XeeX_spatial = Σ_t  Σ_{i,j∈units}    K_space(d_ij/h)             · X_{i,t} ε_{i,t} ε_{j,t} X_{j,t}'
    XeeX_serial  = Σ_u  Σ_{|t-s|≤L,t≠s}  (1 - |t-s|/(L+1))           · X_{u,t} ε_{u,t} ε_{u,s} X_{u,s}'
    Var̂(β)       = (X'X)^{-1} · ( XeeX_spatial + XeeX_serial ) · (X'X)^{-1}

The temporal kernel is hardcoded Bartlett-style regardless of `conley_kernel`
(matches `conleyreg::time_dist.cpp`).

**Kernels:**
- `"bartlett"` (default): `K(u) = max(0, 1 - |u|)` on pairwise distance `d_ij/h`. The radial 1-D form, matching R `conleyreg` / Stata `acreg`. Conley 1999's explicit PSD-guaranteed Bartlett formula (Eq 3.14) is the 2-D **separable product** window on a lattice; the 1-D radial specialization that diff-diff implements is a practitioner convention and is not formally PSD-guaranteed.
- `"uniform"`: `K(u) = 1{|u| ≤ 1}`. Easier to interpret.

Both kernels: `UserWarning` is emitted if the resulting meat has a materially negative eigenvalue (< -1e-12) — neither kernel is formally PSD-guaranteed in the radial 1-D pairwise-distance form.

**Distance metrics:**
- `"haversine"` (default): great-circle in km, Earth's mean radius 6371.01 km (matching R `conleyreg`). Validates `lat ∈ [-90, 90]`, `lon ∈ [-180, 180]`.
- `"euclidean"`: from projected coordinates; user owns the units.
- `callable(coords1, coords2) -> n×n array`: custom distance for non-geographic networks.

**No default bandwidth.** `conley_cutoff_km` is required. Conley (1999) Section 5
recommends a sensitivity grid (e.g., 50, 100, 200, 500 km) and reporting the
SE range.

**Restrictions in this release:**
- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `+ vcov_type="conley"` without `conley_lag_cutoff` → `ValueError` (no defensible default; explicit user choice required).
- `DifferenceInDifferences` / `MultiPeriodDiD` `+ vcov_type="conley"` without `unit=` at fit-time → `ValueError`.
- Combining `vcov_type="conley"` with explicit `cluster=<col>` applies the combined spatial + cluster product kernel (Wave A #119). On the panel path the cluster must be constant within each unit across periods (validator raises `ValueError` otherwise). TWFE's default auto-cluster on the Conley path is silently dropped; users opt into the combined kernel explicitly.
- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley", inference="wild_bootstrap")` → `NotImplementedError` (wild bootstrap does not consume the analytical sandwich).
- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley")` + `survey_design=` → `NotImplementedError` (estimator-level gate; weighted spatial-HAC under probability sampling is an open methodological question).
- `SyntheticDiD(vcov_type="conley")` → `TypeError` (uses bootstrap, not analytical sandwich).
- Generic `LinearRegression(vcov_type="conley", survey_design=...)` → `NotImplementedError`. Generic `LinearRegression` / `compute_robust_vcov` Conley rejects `weights=` for any `weight_type` (`pweight` / `aweight` / `fweight`) → `NotImplementedError` (weighted Conley is not implemented on the generic linalg surface; `compute_robust_vcov` does not accept `survey_design=`, the survey-design surface is `LinearRegression` only). The `pweight` / `survey_design` subset additionally reflects an open methodological question — no canonical extension of Conley (1999) exists for weighted spatial-HAC under probability sampling.
- A sparse k-d-tree fast path auto-activates for `n > 5_000` with `conley_kernel="bartlett"` AND `conley_metric` in `{"haversine", "euclidean"}` (Wave A #120); callable metrics and uniform kernel fall back to the dense path. `n > 20_000` with the dense fallback still emits a memory-OOM `UserWarning`.
- Callable `conley_metric` is validated at the boundary (shape `(n, n)`, finite, non-negative, symmetric to `atol=1e-10`, AND zero on the diagonal `|d(i, i)| ≤ 1e-10` so the kernel reduces to `K(0) = 1` on the HC0 contribution); each failure raises a targeted `ValueError` (Wave A #123).

**Parity:** matches R `conleyreg` (Düsterhöft 2021, CRAN v0.1.9) to ≤ 1e-6
on six benchmark fixtures in
`benchmarks/data/r_conleyreg_conley_golden.json`: three cross-sectional and
three panel fixtures with `lag_cutoff > 0` (`panel_haversine_lag1`,
`panel_haversine_lag2`, `panel_lat_lon_realistic_lag1`). Earth radius
6371.01 km matches conleyreg.

## Rust Backend

diff-diff includes an optional Rust backend for performance-critical operations.

```python
from diff_diff import HAS_RUST_BACKEND

if HAS_RUST_BACKEND:
    print("Rust backend available - computations will be faster")
```

The Rust backend accelerates: OLS solving, robust VCV computation, bootstrap weight generation, synthetic control weights, and simplex projection. It is used transparently when available. Force backend selection via environment variables:

```bash
DIFF_DIFF_BACKEND=python pytest   # Force pure Python
DIFF_DIFF_BACKEND=rust pytest     # Force Rust (fail if unavailable)
```

## Choosing an Estimator

| Scenario | Recommended Estimator |
|----------|----------------------|
| Classic 2x2 design (one treated group, one time split) | `DifferenceInDifferences` |
| Panel data with unit + time FE | `TwoWayFixedEffects` |
| Event study with multiple periods (simultaneous adoption) | `TwoWayFixedEffects` with `event_study=True` (`spec="pooled"` reproduces the deprecated `MultiPeriodDiD`) |
| Staggered treatment timing | `CallawaySantAnna`, `ImputationDiD`, or `SunAbraham` |
| Few treated units / synthetic control | `SyntheticDiD` |
| Interactive fixed effects / factor confounding | `TROP` |
| Continuous treatment intensity, per-dose ATT(d) / ACRT(d) (needs zero-dose controls, or `control_group="lowest_dose"` for Remark 3.1 when P(D=0)=0) | `ContinuousDiD` |
| Continuous treatment intensity, WAS at dose support boundary (compatible with universal rollout or small never-treated share) | `HeterogeneousAdoptionDiD` |
| Two-criterion treatment, simultaneous (2x2x2 DDD) | `TripleDifference` |
| Two-criterion treatment, staggered timing + eligibility | `StaggeredTripleDifference` |
| Nonlinear outcome (binary/count) with staggered timing | `WooldridgeDiD` |
| Diagnosing TWFE bias | `BaconDecomposition` |
| Efficiency-optimal estimation | `EfficientDiD` |
| Corrective weighting for stacked regressions | `StackedDiD` |
| Robustness to parallel trends violations | `HonestDiD` |

## BusinessReport

Plain-English stakeholder narrative from any fitted result type.
Renders `summary()` (short paragraph), `full_report()` (multi-section
markdown), and `to_dict()` (stable AI-legible schema — single source of
truth; prose renders from the dict).

```python
from diff_diff import BusinessReport

report = BusinessReport(
    results,
    outcome_label="Revenue per user",
    outcome_unit="$",  # "$" / "%" / "pp" / "log_points" / "count" recognized
    outcome_direction="higher_is_better",
    business_question="Did the campaign lift revenue?",
    treatment_label="the campaign",
    alpha=0.05,  # single knob: drives both CI level and phrasing threshold
    auto_diagnostics=True,  # default; auto-constructs DiagnosticReport
)

print(report.summary())       # 6-10 sentence paragraph
print(report.full_report())   # structured markdown
report.to_dict()              # AI-legible schema; stable top-level keys
```

Constructor rejects `BaconDecompositionResults` with a helpful TypeError
(Bacon is a diagnostic, not an estimator; wrap the underlying estimator
and pass the Bacon object to `DiagnosticReport(precomputed={'bacon': ...})`).

Schema top-level keys (all always present; missing content uses a
`{"status": "skipped", "reason": "..."}` shape rather than being absent):

- `schema_version`, `estimator`, `context`
- `headline`, `target_parameter`, `assumption`, `pre_trends`, `sensitivity`
- `sample`, `heterogeneity`, `robustness`, `diagnostics`
- `next_steps`, `caveats`, `references`

`target_parameter` (experimental) names what the headline scalar
represents for each estimator — overall ATT, DID_M, DID_1, cost-
benefit delta, dose-response aggregate, ASF-based ETWFE, etc.
Fields: `name` (short stakeholder label), `definition` (full prose),
`aggregation` (machine-readable dispatch tag, e.g., `"simple"`,
`"event_study"`, `"delta"`, `"no_scalar_headline"`),
`headline_attribute` (which raw attribute holds the scalar, or
`None` when no scalar exists by design — e.g., dCDH with
`trends_linear=True` and `L_max>=2`), `reference` (citation).

Status enum values: `ran | skipped | error | not_applicable | not_run | computed`.
`headline.status` also supports `"no_scalar_by_design"` on the dCDH
no-scalar branch.

## DiagnosticReport

Unified diagnostic runner orchestrating `check_parallel_trends`,
`compute_pretrends_power`, `HonestDiD.sensitivity`, `BaconDecomposition`,
plus estimator-native surfaces for SyntheticDiD (`pre_treatment_fit`,
`get_weight_concentration`, `in_time_placebo`, `sensitivity_to_zeta_omega`)
and TROP (factor-model metrics). EfficientDiD PT uses the native
`hausman_pretest`. The `design_effect` section is read-only: it
echoes `survey_metadata.design_effect` / `effective_n` from the
fitted result along with a plain-English band label. The
`epv` section is similarly read-only, reporting from
`results.epv_diagnostics` plus `results.epv_threshold`.

```python
from diff_diff import DiagnosticReport

dr = DiagnosticReport(
    results,
    data=df,  # optional; needed for 2x2 PT, Bacon-from-scratch
    outcome="outcome",
    unit="unit",
    time="period",
    first_treat="first_treat",
    alpha=0.05,
    # Opt-outs (all default True except placebo)
    run_parallel_trends=True,
    run_sensitivity=True,
    run_placebo=False,          # opt-in; not implemented in MVP
    run_bacon=True,
    run_design_effect=True,
    run_heterogeneity=True,
    run_epv=True,
    run_pretrends_power=True,   # drives power-aware PT phrasing
    sensitivity_M_grid=(0.5, 1.0, 1.5, 2.0),
    sensitivity_method="relative_magnitude",
    # Escape hatch for users who ran a diagnostic with custom args:
    precomputed={"sensitivity": my_honest_did_results},
)

dr.run_all()             # triggers compute, caches
print(dr.summary())      # overall-interpretation paragraph
dr.to_dict()             # AI-legible schema
dr.to_dataframe()        # one row per check
dr.applicable_checks     # tuple of checks that will run for this estimator
dr.skipped_checks        # dict of {check: plain-English reason}
```

Schema top-level keys: `schema_version, estimator, headline_metric,
target_parameter, parallel_trends, pretrends_power, sensitivity,
placebo, bacon, design_effect, heterogeneity, epv,
estimator_native_diagnostics, skipped, warnings,
overall_interpretation, next_steps`. The `target_parameter` block
mirrors BR's (same shape and dispatch); see BR's section above for
field semantics including the `headline_attribute=None` /
`aggregation="no_scalar_headline"` case for dCDH
`trends_linear=True, L_max>=2` fits.

### Verdicts and tiers

Pre-trends verdict (three bins, documented in `docs/methodology/REPORTING.md`):

- `joint_p >= 0.30` -> `no_detected_violation`
- `0.05 <= joint_p < 0.30` -> `some_evidence_against`
- `joint_p < 0.05` -> `clear_violation`

Power tier (drives BR phrasing for the `no_detected_violation` verdict):

- `mdv / |att| < 0.25` -> `well_powered`
- `0.25 <= mdv / |att| < 1.0` -> `moderately_powered`
- `mdv / |att| >= 1.0` -> `underpowered`
- power not runnable -> `unknown` (BR falls back to underpowered phrasing)

### Methodology notes

BR and DR do no estimator fitting and do not re-derive variance from
raw data — every effect, SE, p-value, CI, and sensitivity bound is
read from the fitted result or produced by an existing diff-diff
utility (may call `check_parallel_trends`, `BaconDecomposition.fit`, or
`EfficientDiD.hausman_pretest` when the panel + column kwargs are
supplied). The `design_effect` section is read-only: it echoes
`survey_metadata.design_effect` / `effective_n` from the fitted
result rather than calling `compute_deff_diagnostics`. Report-layer
cross-period aggregations are enumerated in
`docs/methodology/REPORTING.md`. Both schemas are experimental in the
current release; see that document for phrasing rules, the
no-traffic-light decision, unit-translation policy, and schema
stability policy.

## MMM Calibration Export

Interop builders converting experiment results into Marketing Mix Model
(MMM) calibration inputs. Pure numpy/pandas - no MMM package is imported,
no result object is introspected (the module is purely additive). Design:
EXPLICIT IN, VALIDATED OUT. Reconciling an experiment estimate to a
calibration input needs the target MMM's row granularity (per-geo vs
national), its time window, and the outcome scale (additive levels vs
log/rate/share) - none of which diff-diff can see - so the CALLER supplies
the already-scoped incremental outcome and its SE (read off summary(),
aggregated to the population/window one MMM row represents). diff-diff
assembles the schema, enforces each consumer's guards, converts to the
lognormal parameterization, pools, and emits snippets. Deriving totals
from a fit is deferred to the post-4.0 results.aggregate() layer, where
the estimator owns its aggregation weights/balance/survey masses.

```python
from diff_diff import (
    SyntheticDiD,
    to_pymc_marketing_lift_test,
    to_meridian_roi_prior,
)

# Single-treated-geo experiment: with exactly ONE treated geo the SDID ATT
# is that geo's lift, so a geo-labelled row is sound. For MULTIPLE treated
# geos the pooled ATT is an average - export each geo's own effect or omit
# the geo dim and match an aggregate lift to aggregate spend.
result = SyntheticDiD().fit(panel, outcome="revenue", treatment="treated",
                            unit="geo", time="week")

# PyMC-Marketing / prophetverse lift-test frame. Pass the scoped effect.
df_lift = to_pymc_marketing_lift_test(
    channel="tv",              # must match the MMM's channel_columns
    x=50_000.0,                # baseline spend for this row's scope
    delta_x=20_000.0,          # spend change (nonzero; negative go-dark,
                               # with x + delta_x >= 0)
    delta_y=result.att,        # measured lift, scoped to this row
    sigma=result.se,           # its SE (finite, > 0)
    dims={"geo": "US-CA"},     # optional model-dim coordinate (one treated geo)
    on_wrong_sign="raise",     # raise|drop|keep - PyMC rejects
                               # sign(delta_y)!=sign(delta_x) AND delta_y==0
                               # (degenerate for its Gamma lift likelihood)
)

# Google Meridian lognormal prior. Caller aggregates the ATT to a total
# incremental outcome over the treated population/window and its SE.
prior = to_meridian_roi_prior(
    incremental_outcome=180_000.0,     # total incremental revenue
    incremental_outcome_se=45_000.0,   # its SE
    spend=200_000.0,                   # channel spend the outcome is over
    parameter="roi_m",                 # roi_m (full-spend/zero-spend return)
                                       # | mroi_m (marginal return)
    se_widening=1.5,                   # >=1 for transferability skepticism
)
prior.roi_mean, prior.roi_sd          # pooled ROI moments
prior.mu, prior.sigma                 # LogNormal params (match Google's
                                      # lognormal_dist_from_mean_std)
prior.to_dict()                       # JSON-ready
print(prior.to_code(                  # ready-to-paste PriorDistribution +
    channel="tv",                     # ModelSpec; roi_m/mroi_m is per-channel
    media_channels=["search", "tv"],  # so channel scope is required (vector
                                      # prior in model channel order; other
                                      # channels keep the Meridian default),
    roi_calibration_period="mask",    # AND time scope (mask expr) or
                                      # full_model_window=True; sets
))                                    # media_prior_type accordingly.
```

Key points:
- Lift frame: one row per experiment, columns [channel, *dims, x, delta_x,
  delta_y, sigma]; guards on sigma>0, delta_x!=0, x>=0, x+delta_x>=0,
  finite delta_y; dims must not collide with reserved columns and must
  share one key set. on_wrong_sign governs wrong-sign AND zero-lift rows
  (drop-all raises; keep warns the frame is not valid PyMC input).
- Meridian: roi = incremental_outcome/spend, spend-weighted pooling,
  roi_sd = sqrt(sum((w_i*sd_i)^2)) (independence assumed - widen via
  se_widening); non-positive pooled ROI raises (lognormal positivity);
  all scaled/pooled outputs re-validated finite-positive.
- parameter="roi_m" vs "mroi_m" picks the Meridian estimand + its default
  for non-experiment channels (LogNormal(0.2,0.9) vs LogNormal(0.0,0.5))
  and the emitted media_prior_type ("roi" vs "mroi").
