-
Notifications
You must be signed in to change notification settings - Fork 49
Expand file tree
/
Copy pathchaisemartin_dhaultfoeuille.py
More file actions
8636 lines (7990 loc) · 392 KB
/
Copy pathchaisemartin_dhaultfoeuille.py
File metadata and controls
8636 lines (7990 loc) · 392 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
de Chaisemartin-D'Haultfoeuille (dCDH) estimator for reversible-treatment DiD.
The dCDH estimator is the only modern DiD estimator in the diff-diff library
that handles **non-absorbing (reversible) treatments** — treatment can switch
on AND off over time. All other staggered estimators in the library
(``CallawaySantAnna``, ``SunAbraham``, ``ImputationDiD``, ``TwoStageDiD``,
``EfficientDiD``, ``WooldridgeDiD``) assume treatment is absorbing.
Phase 1 ships the contemporaneous-switch case ``DID_M`` (= ``DID_1`` at
horizon ``l = 1`` of the dynamic companion paper). Phases 2 and 3 add
dynamic horizons and covariates respectively, on the *same* class — see
``ROADMAP.md`` for the full progression. The forward-compatibility
parameters in :meth:`ChaisemartinDHaultfoeuille.fit` raise
``NotImplementedError`` with phase pointers until later phases land.
References
----------
- de Chaisemartin, C. & D'Haultfoeuille, X. (2020). Two-Way Fixed Effects
Estimators with Heterogeneous Treatment Effects. *American Economic
Review*, 110(9), 2964-2996.
- de Chaisemartin, C. & D'Haultfoeuille, X. (2022, revised 2023).
Difference-in-Differences Estimators of Intertemporal Treatment Effects.
NBER Working Paper 29873. Web Appendix Section 3.7.3 contains the
cohort-recentered plug-in variance formula implemented here.
"""
import warnings
from typing import Any, Dict, List, Optional, Sequence, Set, Tuple
import numpy as np
import pandas as pd
from diff_diff.chaisemartin_dhaultfoeuille_bootstrap import (
ChaisemartinDHaultfoeuilleBootstrapMixin,
)
from diff_diff.chaisemartin_dhaultfoeuille_results import (
ChaisemartinDHaultfoeuilleResults,
DCDHBootstrapResults,
)
from diff_diff.linalg import solve_ols
from diff_diff.utils import safe_inference
__all__ = [
"ChaisemartinDHaultfoeuille",
"chaisemartin_dhaultfoeuille",
"twowayfeweights",
"TWFEWeightsResult",
]
# =============================================================================
# Public dataclass for the standalone TWFE diagnostic helper
# =============================================================================
class TWFEWeightsResult:
"""
Lightweight container for the standalone ``twowayfeweights`` helper.
Returned by :func:`twowayfeweights`. Mirrors the per-cell decomposition
information that the dCDH estimator stores on its results object when
``twfe_diagnostic=True``, but available as a standalone function for
users who only want the diagnostic without fitting the full estimator.
"""
__slots__ = ("weights", "fraction_negative", "sigma_fe", "beta_fe")
def __init__(
self,
weights: pd.DataFrame,
fraction_negative: float,
sigma_fe: float,
beta_fe: float,
) -> None:
self.weights = weights
self.fraction_negative = fraction_negative
self.sigma_fe = sigma_fe
self.beta_fe = beta_fe
def __repr__(self) -> str:
return (
f"TWFEWeightsResult(beta_fe={self.beta_fe:.4f}, "
f"fraction_negative={self.fraction_negative:.4f}, "
f"sigma_fe={self.sigma_fe:.4f}, n_cells={len(self.weights)})"
)
# =============================================================================
# Shared validation + cell aggregation helper
# =============================================================================
def _validate_and_aggregate_to_cells(
data: pd.DataFrame,
outcome: str,
group: str,
time: str,
treatment: str,
weights: Optional[np.ndarray] = None,
) -> pd.DataFrame:
"""
Validate input data and aggregate to ``(g, t)`` cells per the dCDH contract.
Used by both :meth:`ChaisemartinDHaultfoeuille.fit` and
:func:`twowayfeweights` so the validation rules and aggregation
behavior are identical across the two public entry points.
The contract (matching ``REGISTRY.md`` ``## ChaisemartinDHaultfoeuille``):
1. **Required columns** ``outcome``, ``group``, ``time``, ``treatment``
must all be present in ``data`` (raises ``ValueError`` listing
any missing).
2. **Treatment** must coerce to numeric and contain no ``NaN``
(raises ``ValueError`` — silent dropping would change cell counts
without informing the user).
3. **Outcome** must coerce to numeric and contain no ``NaN`` (same
reasoning).
4. **Treatment** must be numeric. Both binary ``{0, 1}`` and
non-binary (ordinal or continuous) treatment are supported.
Non-binary treatment requires ``L_max >= 1`` in ``fit()`` because
the per-period DID path uses binary joiner/leaver categorization.
5. **Cell aggregation** via ``groupby([group, time]).agg(...)``
producing ``y_gt`` (cell mean of ``outcome``), ``d_gt`` (cell
mean of ``treatment``), and ``n_gt`` (count of original
observations in the cell).
6. **Within-cell-varying treatment** (any cell where ``d_min !=
d_max``) raises ``ValueError``. Treatment must be constant
within each ``(group, time)`` cell; fuzzy DiD is deferred to a
separate dCdH 2018 paper. Pre-aggregate your data to constant
cell-level treatment before calling ``fit()`` or
``twowayfeweights()``.
When ``weights`` is provided (survey pweights), cell means use weighted
averages: ``y_gt = sum(w_i * y_i) / sum(w_i)``. An additional column
``w_gt`` (total weight per cell) is included in the output for downstream
IF expansion.
Returns the aggregated cell DataFrame with columns
``[group, time, y_gt, d_gt, n_gt]`` (plus ``w_gt`` when weighted),
sorted by ``[group, time]`` with a fresh index.
Raises
------
ValueError
On missing columns; NaN values in any of the ``group``, ``time``,
``treatment``, or ``outcome`` columns (``group`` and ``time`` are
rejected pre-``groupby`` because ``groupby`` silently drops NaN
keys, which would change the estimation sample without warning);
non-numeric treatment / outcome that cannot be coerced via
``pd.to_numeric``; or within-cell-varying treatment (any
``(group, time)`` cell where ``d_min != d_max``, since fuzzy DiD
is out of scope and deferred to a separate dCdH 2018 paper).
Integer-coded non-binary treatment (the ``by_path`` /
``paths_of_interest`` requirement) is enforced separately at
``fit()`` time, not here at aggregation time — this helper
accepts continuous ``d_gt`` cell means and lets ``fit()`` decide
whether the integer-only contract applies.
Under the survey-weighted path (``weights`` is not ``None``),
zero-weight rows are pre-filtered before any NaN / coercion /
within-cell validation per the ``SurveyDesign.subpopulation()``
out-of-sample contract — invalid values in zero-weight rows
therefore do NOT raise. NaN / coercion / within-cell checks
still apply to all positive-weight rows.
"""
# 1. Required columns
missing = [c for c in (outcome, group, time, treatment) if c not in data.columns]
if missing:
raise ValueError(
f"ChaisemartinDHaultfoeuille / twowayfeweights: column(s) {missing!r} "
f"not found in data. Required columns: outcome, group, time, treatment."
)
df = data.copy()
# 1a. SurveyDesign.subpopulation() contract: zero-weight rows are
# out-of-sample. Pre-filter them *before* any NaN/coercion validation
# so that invalid values in excluded rows do not abort the fit.
if weights is not None:
weights_arr = np.asarray(weights, dtype=np.float64)
pos_mask = weights_arr > 0
if not pos_mask.all():
df = df.loc[pos_mask].reset_index(drop=True)
weights = weights_arr[pos_mask]
# 1b. Group and time NaN checks (before groupby, which silently drops NaN keys)
n_nan_group = int(df[group].isna().sum())
if n_nan_group > 0:
raise ValueError(
f"Group column {group!r} contains {n_nan_group} NaN value(s). "
"groupby silently drops NaN keys, which would change the "
"estimation sample without warning. Drop or impute NaN group "
"values before calling fit() or twowayfeweights()."
)
n_nan_time = int(df[time].isna().sum())
if n_nan_time > 0:
raise ValueError(
f"Time column {time!r} contains {n_nan_time} NaN value(s). "
"groupby silently drops NaN keys, which would change the "
"estimation sample without warning. Drop or impute NaN time "
"values before calling fit() or twowayfeweights()."
)
# 2. Treatment numeric coercion + NaN check
try:
df[treatment] = pd.to_numeric(df[treatment])
except (ValueError, TypeError) as exc:
raise ValueError(
f"Could not coerce treatment column {treatment!r} to numeric: {exc}"
) from exc
n_nan_treat = int(df[treatment].isna().sum())
if n_nan_treat > 0:
raise ValueError(
f"Treatment column {treatment!r} contains {n_nan_treat} NaN value(s). "
"ChaisemartinDHaultfoeuille requires non-missing treatment indicators "
"on every observation; impute or drop NaN treatment rows before fitting "
"so the dropped count is explicit."
)
# 3. Outcome numeric coercion + NaN check
try:
df[outcome] = pd.to_numeric(df[outcome])
except (ValueError, TypeError) as exc:
raise ValueError(f"Could not coerce outcome column {outcome!r} to numeric: {exc}") from exc
n_nan_outcome = int(df[outcome].isna().sum())
if n_nan_outcome > 0:
raise ValueError(
f"Outcome column {outcome!r} contains {n_nan_outcome} NaN value(s). "
"Drop or impute missing outcomes before calling fit() so the "
"exclusion is explicit (silently averaging over present values "
"would distort per-cell means)."
)
# 4. Treatment must be numeric (binary or non-binary both accepted)
# No longer enforces {0, 1} - non-binary and continuous treatment supported.
# 5. Cell aggregation (compute min/max for within-cell check)
if weights is not None:
# Survey-weighted cell aggregation (zero-weight rows already
# filtered upstream at step 1a).
# y_gt = sum(w_i * y_i) / sum(w_i) within each (g, t) cell.
# Treatment is constant within cells (checked below), so weighted
# and unweighted means are identical for d_gt.
df["_w_"] = weights
df["_wy_"] = weights * df[outcome].values
g_obj = df.groupby([group, time], as_index=False)
cell = g_obj.agg(
_wy_sum=("_wy_", "sum"),
w_gt=("_w_", "sum"),
d_gt=(treatment, "mean"),
d_min=(treatment, "min"),
d_max=(treatment, "max"),
n_gt=(treatment, "count"),
)
cell["y_gt"] = cell["_wy_sum"] / cell["w_gt"]
cell = cell.drop(columns=["_wy_sum"])
# Zero-weight cells: drop entirely so downstream validators
# (ragged-panel, baseline requirement) don't see them.
zero_w_mask = cell["w_gt"] <= 0
if zero_w_mask.any():
cell = cell[~zero_w_mask].reset_index(drop=True)
df.drop(columns=["_w_", "_wy_"], inplace=True)
else:
cell = df.groupby([group, time], as_index=False).agg(
y_gt=(outcome, "mean"),
d_gt=(treatment, "mean"),
d_min=(treatment, "min"),
d_max=(treatment, "max"),
n_gt=(treatment, "count"),
)
# 6. Within-cell-varying treatment rejection.
# All observations in a cell must have the same treatment value
# (for both binary and non-binary treatment). Detect by checking
# that cell min equals cell max.
non_constant_mask = cell["d_min"] != cell["d_max"]
if non_constant_mask.any():
n_non_constant = int(non_constant_mask.sum())
example_cells = cell.loc[non_constant_mask, [group, time, "d_gt", "d_min", "d_max"]].head(5)
raise ValueError(
f"Within-cell-varying treatment detected in {n_non_constant} "
f"(group, time) cell(s). dCDH requires treatment to be "
f"constant within each (group, time) cell. Cells where "
f"d_min != d_max indicate that some units have different "
f"treatment values. Pre-aggregate your data to constant "
f"cell-level treatment before calling fit() or "
f"twowayfeweights(). Fuzzy DiD is deferred to a separate "
f"dCDH paper (see ROADMAP.md out-of-scope). Affected cells "
f"(first 5):\n{example_cells}"
)
# Drop the min/max columns; keep d_gt as float (no int cast - supports
# ordinal and continuous treatment). w_gt retained when weighted.
drop_cols = ["d_min", "d_max"]
cell = cell.drop(columns=drop_cols)
# Sort to ensure deterministic order in downstream operations
cell = cell.sort_values([group, time]).reset_index(drop=True)
return cell
def _validate_paths_of_interest(
paths_of_interest: Any,
) -> List[Tuple[int, ...]]:
"""Validate and canonicalize ``paths_of_interest`` to ``List[Tuple[int, ...]]``.
Rejects non-sequence inputs, empty lists, non-tuple/list path entries,
empty path entries, non-int elements (including ``bool`` and
``np.bool_``), and entries with mixed lengths. Numpy integer types
(``np.integer``) are accepted and canonicalized to Python ``int``
so the resulting tuples are usable as dict keys interchangeably with
paths emitted by ``_enumerate_treatment_paths`` (which casts via
``int(round(float(v)))``).
"""
if not isinstance(paths_of_interest, (list, tuple)):
raise ValueError(
f"paths_of_interest must be a list/tuple of int tuples, "
f"got {type(paths_of_interest).__name__}."
)
if len(paths_of_interest) == 0:
raise ValueError("paths_of_interest must be non-empty.")
canonical: List[Tuple[int, ...]] = []
for i, p in enumerate(paths_of_interest):
if not isinstance(p, (list, tuple)):
raise ValueError(
f"paths_of_interest[{i}] must be a tuple/list of ints, " f"got {type(p).__name__}."
)
if len(p) == 0:
raise ValueError(f"paths_of_interest[{i}] must be non-empty.")
canonical_path: List[int] = []
for j, v in enumerate(p):
if isinstance(v, (bool, np.bool_)) or not isinstance(v, (int, np.integer)):
raise ValueError(
f"paths_of_interest[{i}][{j}] must be an int, got "
f"{v!r} of type {type(v).__name__}."
)
canonical_path.append(int(v))
canonical.append(tuple(canonical_path))
lens = {len(p) for p in canonical}
if len(lens) > 1:
raise ValueError(
f"paths_of_interest entries must all have the same length "
f"(L_max+1); got mixed lengths {sorted(lens)}."
)
return canonical
# =============================================================================
# Main estimator class
# =============================================================================
class ChaisemartinDHaultfoeuille(ChaisemartinDHaultfoeuilleBootstrapMixin):
"""
de Chaisemartin-D'Haultfoeuille (dCDH) estimator.
The only modern DiD estimator in the library that handles **reversible
(non-absorbing) treatments** - treatment may switch on AND off over
time. Computes the contemporaneous-switch DiD ``DID_M`` from the
AER 2020 paper (equivalently ``DID_1`` at horizon ``l = 1`` of the
dynamic companion paper, NBER WP 29873) plus the full multi-horizon
event study ``DID_l`` for ``l = 1..L_max`` via the ``L_max`` parameter
on :meth:`fit`.
Supported:
- Headline ``DID_M`` plus multi-horizon ``DID_l`` event study
- Joiners-only ``DID_+`` and leavers-only ``DID_-`` decompositions
- Single-lag placebo ``DID_M^pl`` and dynamic placebos ``DID^{pl}_l``
(computed automatically by default; gate via ``placebo=False``)
- Analytical SE via the cohort-recentered plug-in formula from Web
Appendix Section 3.7.3; multiplier bootstrap clustered at the group
level by default via ``n_bootstrap``; under ``survey_design`` with
strictly-coarser PSUs the bootstrap automatically upgrades to
PSU-level Hall-Mammen wild clustering (see REGISTRY.md
``ChaisemartinDHaultfoeuille`` Note on survey + bootstrap)
- Normalized estimator ``DID^n_l``, cost-benefit aggregate ``delta``,
and sup-t simultaneous confidence bands
- Residualization-style covariate adjustment (``DID^X``) via
``controls=``, group-specific linear trends (``DID^{fd}``) via
``trends_linear=True``, state-set-specific trends via
``trends_nonparam=``, heterogeneity testing, non-binary treatment,
HonestDiD sensitivity integration on placebos via ``honest_did=True``
- Per-path event-study disaggregation via ``by_path=k`` (top-k most
common observed treatment paths within the window
``[F_g-1, F_g-1+L_max]``; requires ``drop_larger_lower=False``;
supports binary or integer-coded discrete treatment) or via
``paths_of_interest=[(...), ...]`` for an explicit user-specified
path subset (Python-only API; mutex with ``by_path=k``)
- Survey support via ``survey_design=``: pweight with strata/PSU/FPC
via Taylor Series Linearization (analytical) or replicate-weight
variance (BRR/Fay/JK1/JKn/SDR)
- TWFE decomposition diagnostic from Theorem 1 of AER 2020
Only ``aggregate`` on :meth:`fit` still raises ``NotImplementedError``.
Parameters
----------
alpha : float, default=0.05
Significance level for confidence intervals.
cluster : str, optional, default=None
Must be ``None`` (the default). User-specified clustering via
this kwarg is not supported — passing any non-``None`` value
raises ``NotImplementedError`` at construction time (and the
same gate fires from ``set_params``). The effective clustering
depends on how you call ``fit()``:
- **Default (no survey_design)**: clustered at the group level
via the cohort-recentered influence-function plug-in
(analytical SEs) and the multiplier bootstrap.
- **Under ``survey_design`` with auto-inject or explicit
``psu=group``**: PSU coincides with the group and the
group-level and PSU-level paths are bit-identical.
- **Under ``survey_design`` with strictly-coarser PSUs**: the
multiplier bootstrap automatically upgrades to PSU-level
Hall-Mammen wild clustering.
So dCDH does NOT always cluster at the group level — see
REGISTRY.md ``ChaisemartinDHaultfoeuille`` Notes on cluster
contract and survey + bootstrap for the full matrix. Custom
user-specified clustering at a coarser or finer level than the
group is a planned extension.
n_bootstrap : int, default=0
Number of multiplier-bootstrap iterations. ``0`` (default) uses
only the analytical SE. Set to ``999`` or higher for stable
bootstrap inference.
bootstrap_weights : str, default="rademacher"
Type of multiplier-bootstrap weights: ``"rademacher"``,
``"mammen"``, or ``"webb"``. Ignored unless ``n_bootstrap > 0``.
seed : int, optional
Random seed for the multiplier bootstrap.
placebo : bool, default=True
If ``True`` (default), automatically compute the single-lag
placebo ``DID_M^pl`` (AER 2020 placebo specification) on the same data.
Set to ``False`` to skip the placebo computation for speed; the
results object will still expose ``placebo_*`` fields, but with
NaN values and ``placebo_available=False``.
twfe_diagnostic : bool, default=True
If ``True`` (default), compute the TWFE decomposition diagnostic
from Theorem 1 of AER 2020: per-``(g, t)`` weights, fraction of
treated cells with negative weights, and ``sigma_fe`` (the
smallest cell-effect standard deviation that could flip the sign
of the plain TWFE coefficient). The diagnostic answers "what
would the plain TWFE estimator say on the data you passed in?",
so it runs on the **FULL pre-filter cell sample** (the same
input as the standalone :func:`twowayfeweights` function), NOT
on the post-filter estimation sample used by ``DID_M``. When
the ragged-panel filter or ``drop_larger_lower`` drops groups,
the fitted ``results.twfe_*`` values describe a LARGER sample
(pre-filter) than ``results.overall_att`` and a ``UserWarning``
is emitted to make the divergence explicit. See REGISTRY.md
``ChaisemartinDHaultfoeuille`` ``Note (TWFE diagnostic sample
contract)`` for the full rationale.
drop_larger_lower : bool, default=True
If ``True`` (default, matches R ``DIDmultiplegtDYN``), drops
groups whose treatment switches more than once (multi-switch
groups) before estimation. This is required for the analytical
variance formula to be consistent with the AER 2020 Theorem 3
point estimate — both formulas operate on the same post-drop
dataset. Setting to ``False`` is supported for diagnostic
comparison but produces an inconsistent estimator-variance
pairing for multi-switch groups; a warning is emitted.
by_path : int, optional, default=None
If set to a positive integer ``k``, disaggregate the per-horizon
event study by the observed treatment trajectory in the window
``[F_g - 1, F_g, ..., F_g - 1 + L_max]``, reporting ATT + SE +
inference for the ``k`` most common observed paths (ties broken
lexicographically on the path tuple). If ``k`` exceeds the number
of observed paths, all paths are returned and a ``UserWarning``
is emitted. ``None`` (the default) disables the disaggregation.
Requires ``drop_larger_lower=False`` (multi-switch groups are
the object of interest) and ``L_max >= 1`` (the path window
depends on ``L_max``). Compatible with non-binary integer-coded
treatment (D in Z); path tuples become integer-state tuples
like ``(0, 2, 2, 2)``. D values must be integer-valued
(``D == round(D)``); a ``ValueError`` is raised at fit-time on
continuous D. Compatible with ``survey_design`` for analytical
Binder TSL SE and replicate-weight bootstrap; per-path SE
routes through the cell-period allocator, with non-path
switcher-side contributions skipped (control contributions
remain unchanged, matching the joiners/leavers IF convention).
``n_bootstrap > 0`` (multiplier bootstrap) under
``survey_design`` is not yet supported and raises
``NotImplementedError``. Top-k path ranking under
``survey_design`` remains group-cardinality-based (unweighted),
not population-weight-based — survey weights do not affect
which paths are selected as "top-k".
Compatible with ``heterogeneity="<col>"`` — per-path
heterogeneity coefficient is computed by re-running the
Lemma 7 regression on each path-restricted switcher
subsample. Cohort dummies absorb baseline (no R-divergence
warning needed). Surfaces on
``results.path_heterogeneity_effects`` keyed
``{path: {l: {beta, se, t_stat, p_value, conf_int, n_obs}}}``
and on ``to_dataframe(level="by_path")`` via ``het_*``
columns. Mirrors R ``did_multiplegt_dyn(..., by_path,
predict_het)`` per-by_level. Composes with ``survey_design``
(analytical Binder TSL + replicate-weight) via the existing
cell-period IF allocator path. Incompatible with
``design2`` and ``honest_did`` (each combination raises
``NotImplementedError`` in the current release).
Mutually exclusive with ``paths_of_interest`` — use
``by_path=k`` for top-k automatic ranking by frequency, or
``paths_of_interest=[(...), ...]`` for an explicit user-
specified path list. Setting both raises ``ValueError``.
Compatible with ``controls`` (DID^X residualization) -- the
per-baseline OLS residualization runs once on first-differenced
``Y`` BEFORE path enumeration, so per-path point estimates,
bootstrap SE, per-path placebos, and per-path sup-t bands all
consume the residualized ``Y_mat`` automatically (Frisch-
Waugh-Lovell). Per-period effects remain unadjusted, consistent
with the existing ``controls`` + per-period DID contract.
**Deviation from R on multi-baseline switcher panels:** R
``did_multiplegt_dyn(..., by_path, controls)`` re-runs the
per-baseline residualization on each path's restricted
subsample (path's switchers + same-baseline not-yet-treated
controls), so its residualization coefficients vary per path
when switchers have different baseline values. Our global-
residualization architecture coincides with R on single-
baseline panels (every switcher shares the same ``D_{g,1}``)
and per-path point estimates match exactly on the one-
observation-per-``(g, t)`` regime; on multi-observation-per-
cell panels the existing DID^X cell-weighting deviation from
R applies (see ``docs/methodology/REGISTRY.md`` "Note (Phase
3 DID^X covariate adjustment)"; independent of the by_path
lift). On multi-baseline switcher panels, point estimates can
diverge — a ``UserWarning`` is emitted at fit-time when this
configuration is detected. SE inherits the cross-path cohort-
sharing deviation from R documented for ``path_effects``.
Compatible with ``trends_linear`` (DID^{fd} group-specific
linear trends) -- first-differencing replaces ``Y`` with
``Z = Y_t - Y_{t-1}`` once globally before path enumeration,
so per-path raw second-differences DID^{fd}_{path, l} surface
on ``path_effects[path]["horizons"][l]`` automatically. Per-path
cumulated level effects ``delta_{path, l} = sum_{l'=1..l}
DID^{fd}_{path, l'}`` are surfaced on the new
``results.path_cumulated_event_study[path][l]`` field
(mirroring the global ``linear_trends_effects`` cumulation;
inner dict keyed by horizon directly, no ``"horizons"`` wrapper).
SE on the cumulated layer is the conservative upper bound
(sum of per-horizon component SEs, NaN-consistent), matching
the global ``linear_trends_effects`` SE convention. Path
enumeration runs on the post-first-differenced ``N_mat_fd``:
switchers with ``F_g==2`` fail the window-eligibility check
and are dropped from path enumeration entirely, so a path
whose switchers all have ``F_g < 3`` is silently absent from
``path_effects`` (the existing global ``F_g < 3`` warning
still fires). Per-path R parity matches R
``did_multiplegt_dyn(..., by_path, trends_lin)`` on per-path
cumulated point estimates under single-baseline panels with
sufficient pre-window depth (``F_g >= 4`` for every selected-
path switcher). R re-runs the per-path full pipeline on each
path's restricted subsample; same multi-baseline divergence
pattern as ``controls`` (a ``UserWarning`` fires when switcher
baselines take multiple values). **F_g=3 boundary-case
divergence:** `F_g=3` switchers have only 1 valid pre-window
Z value after first-differencing and the ``time==1`` filter,
which causes Python's global-then-disaggregate architecture
to diverge from R's per-path full-pipeline call (30%+ on
point estimates observed empirically). A separate
``UserWarning`` fires at fit-time when the panel includes any
`F_g=3` switchers and `by_path + trends_linear` is set, so
practitioners hitting this boundary regime see the divergence
flag explicitly. **Placebo under trends_linear returns RAW
per-horizon values, not cumulated** -- there is no per-path
placebo cumulation surface (verified empirically against R
via the existing ``joiners_only_trends_lin`` parity scenario).
Compatible with ``trends_nonparam`` (state-set trends) -- the
set membership column is validated and stored once globally
(time-invariance, NaN rejection, partition coarseness checks
unchanged); per-path analytical SE, bootstrap SE, per-path
placebos, and per-path sup-t bands all inherit the
set-restricted control pool automatically through the
``set_ids`` parameter threaded through the per-path IF
helpers. Per-path R parity matches R
``did_multiplegt_dyn(..., by_path, trends_nonparam)`` on
per-path point estimates under single-baseline panels.
Compatible with ``n_bootstrap > 0`` -- the top-k paths are
enumerated once on the observed data (paths held fixed across
bootstrap draws, matching R ``did_multiplegt_dyn(..., by_path,
bootstrap=B)``) and bootstrap SE / percentile CI / percentile
p-value are written to ``path_effects[path]["horizons"][l]``
in place of the analytical fields. See REGISTRY.md for the
full bootstrap contract.
Compatible with ``placebo=True`` -- when both are active,
per-path backward-horizon placebos ``DID^{pl}_{path, l}`` for
``l = 1..L_max`` are surfaced on
``results.path_placebo_event_study[path][-l]`` (negative-int
keys mirroring ``placebo_event_study``). The same per-path SE
convention is applied backward (joiners/leavers IF precedent;
cohort-recentered plug-in with path-specific divisor); the
cross-path cohort-sharing deviation from R is inherited from
the analytical event-study path.
With ``n_bootstrap > 0``, per-path joint sup-t simultaneous
confidence bands are also computed across horizons
``1..L_max`` within each path. A path-specific critical value
``c_p`` (constructed from a fresh shared-weights multiplier-
bootstrap draw per path) is surfaced at top level as
``results.path_sup_t_bands[path] = {"crit_value", "alpha",
"n_bootstrap", "method", "n_valid_horizons"}``, applied
per-horizon as ``cband_conf_int`` on
``path_effects[path]["horizons"][l]``, and rendered as
``cband_lower`` / ``cband_upper`` columns on
``results.to_dataframe(level="by_path")`` (mirroring the
OVERALL ``level="event_study"`` schema). Bands cover joint
inference WITHIN a single path across horizons; they do NOT
provide simultaneous coverage across paths. Python-only
library extension; R ``did_multiplegt_dyn`` provides no joint
bands at any surface. See REGISTRY.md ``Note (Phase 3 by_path
per-path joint sup-t bands)``.
SE convention: per-path IF parallels the joiners / leavers
construction — the switcher-side contribution is zeroed for
groups not in the selected path, and the cohort structure and
control pool are unchanged. Plug-in SE uses the path-specific
divisor ``N_l_path`` (count of path switchers eligible at horizon
``l``), matching how ``joiners_se`` / ``leavers_se`` use their
respective counts as divisors. See REGISTRY.md
``ChaisemartinDHaultfoeuille`` ``Note`` on ``by_path`` for the
full contract.
Results are exposed on ``results.path_effects`` as a dict keyed
by the path tuple, with nested ``"horizons"`` dicts per
horizon ``l``. Also available via
``results.to_dataframe(level="by_path")``.
paths_of_interest : list of tuple of int, optional, default=None
Explicit user-specified treatment paths to disaggregate by, as
an alternative to ``by_path=k``'s top-k automatic ranking.
Each path tuple must have length ``L_max + 1`` and represents
the treatment trajectory in the window
``[F_g - 1, F_g, ..., F_g - 1 + L_max]``, e.g.
``[(0, 1, 1, 1), (0, 1, 0, 0)]`` for two paths under
``L_max=3``. Mutually exclusive with ``by_path``; setting both
raises ``ValueError``.
Validation:
- Each path element must be an ``int`` (``bool`` and
``np.bool_`` rejected; ``np.integer`` accepted and
canonicalized to Python ``int``).
- All paths must have the same length (uniformity validated
at ``__init__``; length match against ``L_max + 1``
validated at fit-time).
- Empty list raises ``ValueError``.
- Duplicate paths are deduplicated with a ``UserWarning``.
- A path with zero observed groups in the panel emits a
``UserWarning`` and is omitted from ``path_effects``.
Compatible with non-binary integer treatment (paths can
contain integer states like ``(0, 2, 2)``).
Compatible with all downstream surfaces inherited by
``by_path``: bootstrap, per-path placebos, per-path joint
sup-t bands, ``controls``, ``trends_linear``,
``trends_nonparam``, ``survey_design`` (analytical Binder
TSL + replicate-weight; multiplier bootstrap under survey
remains gated, same as ``by_path=k``), and ``heterogeneity``
(per-path heterogeneity coefficient surfaces on
``results.path_heterogeneity_effects``). Mechanical
extension to path enumeration; no methodology change.
**Order semantics**: paths appear in
``results.path_effects`` in the user-specified order, modulo
deduplication and unobserved-path filtering.
**Python-only API extension; no R equivalent.** R's
``did_multiplegt_dyn(..., by_path=k)`` only accepts a positive
int (top-k) or ``-1`` (all paths); there is no list-based
path selection in R.
Results expose the same surfaces as ``by_path``:
``results.path_effects`` (dict keyed by path tuple),
``results.path_placebo_event_study``,
``results.path_sup_t_bands``,
``results.path_cumulated_event_study`` (under
``trends_linear``), and the ``level="by_path"`` DataFrame.
rank_deficient_action : str, default="warn"
Action when the TWFE decomposition diagnostic OLS encounters a
rank-deficient design matrix: ``"warn"``, ``"error"``, or
``"silent"``. Only used when ``twfe_diagnostic=True``.
Attributes
----------
results_ : ChaisemartinDHaultfoeuilleResults
Estimation results after calling :meth:`fit`.
is_fitted_ : bool
Whether the model has been fitted.
Notes
-----
The analytical CI is **conservative** under Assumption 8 (independent
groups) of the dynamic companion paper, and exact only under iid
sampling. This is documented as a deliberate deviation from "default
nominal coverage" in ``REGISTRY.md``.
Examples
--------
Basic single-switch panel:
>>> from diff_diff import ChaisemartinDHaultfoeuille
>>> from diff_diff.prep_dgp import generate_reversible_did_data
>>> data = generate_reversible_did_data(n_groups=80, n_periods=6, seed=42)
>>> est = ChaisemartinDHaultfoeuille()
>>> results = est.fit(
... data, outcome="outcome", group="group",
... time="period", treatment="treatment",
... )
>>> abs(results.overall_att - 2.0) < 1.0 # close to the true effect
True
"""
def __init__(
self,
alpha: float = 0.05,
cluster: Optional[str] = None,
n_bootstrap: int = 0,
bootstrap_weights: str = "rademacher",
seed: Optional[int] = None,
placebo: bool = True,
twfe_diagnostic: bool = True,
drop_larger_lower: bool = True,
by_path: Optional[int] = None,
paths_of_interest: Optional[Sequence[Sequence[int]]] = None,
rank_deficient_action: str = "warn",
) -> None:
# Parameter validation
if rank_deficient_action not in ("warn", "error", "silent"):
raise ValueError(
f"rank_deficient_action must be 'warn', 'error', or 'silent', "
f"got '{rank_deficient_action}'"
)
if bootstrap_weights not in ("rademacher", "mammen", "webb"):
raise ValueError(
f"bootstrap_weights must be 'rademacher', 'mammen', or 'webb', "
f"got '{bootstrap_weights}'"
)
if not 0.0 < alpha < 1.0:
raise ValueError(f"alpha must be in (0, 1), got {alpha}")
if n_bootstrap < 0:
raise ValueError(f"n_bootstrap must be non-negative, got {n_bootstrap}")
if by_path is not None:
if isinstance(by_path, bool) or not isinstance(by_path, int):
raise ValueError(
f"by_path must be None or a positive int, got "
f"{by_path!r} of type {type(by_path).__name__}."
)
if by_path <= 0:
raise ValueError(
f"by_path must be a positive int (top-k most common paths), "
f"got {by_path}. Use by_path=None to disable, or "
f"paths_of_interest for explicit path selection."
)
if paths_of_interest is not None:
paths_of_interest = _validate_paths_of_interest(paths_of_interest)
if by_path is not None and paths_of_interest is not None:
raise ValueError(
"by_path and paths_of_interest are mutually exclusive. "
"Use by_path=k for top-k automatic ranking, OR "
"paths_of_interest=[(...), ...] for explicit user-"
"specified paths. Set one and leave the other as None."
)
if cluster is not None:
raise NotImplementedError(
f"cluster={cluster!r}: user-specified clustering is not "
f"supported in ChaisemartinDHaultfoeuille. dCDH clusters at "
f"the group level by default via the cohort-recentered "
f"influence-function plug-in (analytical SEs) and the "
f"multiplier bootstrap. Under survey_design with strictly-"
f"coarser PSUs, bootstrap clustering automatically upgrades "
f"to PSU-level Hall-Mammen wild. To use the default path, "
f"pass cluster=None (the "
f"default). Custom clustering is reserved for a future "
f"phase. See REGISTRY.md ChaisemartinDHaultfoeuille section "
f"for the full contract."
)
self.alpha = alpha
self.cluster = cluster
self.n_bootstrap = n_bootstrap
self.bootstrap_weights = bootstrap_weights
self.seed = seed
self.placebo = placebo
self.twfe_diagnostic = twfe_diagnostic
self.drop_larger_lower = drop_larger_lower
self.by_path = by_path
self.paths_of_interest = paths_of_interest
self.rank_deficient_action = rank_deficient_action
self.is_fitted_ = False
self.results_: Optional[ChaisemartinDHaultfoeuilleResults] = None
# ------------------------------------------------------------------
# sklearn-style parameter introspection
# ------------------------------------------------------------------
def get_params(self) -> Dict[str, Any]:
"""Return all ``__init__`` parameters as a dictionary."""
return {
"alpha": self.alpha,
"cluster": self.cluster,
"n_bootstrap": self.n_bootstrap,
"bootstrap_weights": self.bootstrap_weights,
"seed": self.seed,
"placebo": self.placebo,
"twfe_diagnostic": self.twfe_diagnostic,
"drop_larger_lower": self.drop_larger_lower,
"by_path": self.by_path,
"paths_of_interest": self.paths_of_interest,
"rank_deficient_action": self.rank_deficient_action,
}
def set_params(self, **params: Any) -> "ChaisemartinDHaultfoeuille":
"""
Set estimator parameters (sklearn-compatible).
**Transactional**: validation runs after the candidate mutations,
and if any rule fails the estimator state is rolled back to its
pre-call values before the exception is re-raised. Callers can
therefore retry with corrected params on the same instance
without repairing inconsistent intermediate state.
"""
# Snapshot current values for the keys we are about to set so
# we can roll back on validation failure (transactional semantics).
for key in params:
if not hasattr(self, key):
raise ValueError(f"Unknown parameter: {key}")
snapshot = {key: getattr(self, key) for key in params}
try:
for key, value in params.items():
setattr(self, key, value)
self._validate_invariants()
except Exception:
for key, value in snapshot.items():
setattr(self, key, value)
raise
return self
def _validate_invariants(self) -> None:
"""Run the post-mutation validation rules. Mirrors `__init__`."""
# Re-run __init__ validation rules so the post-set state is valid.
if self.rank_deficient_action not in ("warn", "error", "silent"):
raise ValueError(
f"rank_deficient_action must be 'warn', 'error', or 'silent', "
f"got '{self.rank_deficient_action}'"
)
if self.bootstrap_weights not in ("rademacher", "mammen", "webb"):
raise ValueError(
f"bootstrap_weights must be 'rademacher', 'mammen', or 'webb', "
f"got '{self.bootstrap_weights}'"
)
if not 0.0 < self.alpha < 1.0:
raise ValueError(f"alpha must be in (0, 1), got {self.alpha}")
if self.n_bootstrap < 0:
raise ValueError(f"n_bootstrap must be non-negative, got {self.n_bootstrap}")
if self.by_path is not None:
if isinstance(self.by_path, bool) or not isinstance(self.by_path, int):
raise ValueError(
f"by_path must be None or a positive int, got "
f"{self.by_path!r} of type {type(self.by_path).__name__}."
)
if self.by_path <= 0:
raise ValueError(
f"by_path must be a positive int (top-k most common paths), "
f"got {self.by_path}. Use by_path=None to disable, or "
f"paths_of_interest for explicit path selection."
)
if self.paths_of_interest is not None:
self.paths_of_interest = _validate_paths_of_interest(self.paths_of_interest)
if self.by_path is not None and self.paths_of_interest is not None:
raise ValueError(
"by_path and paths_of_interest are mutually exclusive. "
"Use by_path=k for top-k automatic ranking, OR "
"paths_of_interest=[(...), ...] for explicit user-"
"specified paths. Set one and leave the other as None."
)
if self.cluster is not None:
raise NotImplementedError(
f"cluster={self.cluster!r}: user-specified clustering is "
f"not supported in ChaisemartinDHaultfoeuille. dCDH clusters "
f"at the group level by default; under survey_design with "
f"strictly-coarser PSUs the bootstrap automatically upgrades "
f"to PSU-level Hall-Mammen wild clustering. Pass cluster=None "
f"(the default) to use this path. User-specified custom "
f"clustering is reserved for a future phase. See REGISTRY.md "
f"ChaisemartinDHaultfoeuille section for the full contract."
)
# ------------------------------------------------------------------
# fit
# ------------------------------------------------------------------
def fit(
self,
data: pd.DataFrame,
outcome: str,
group: str,
time: str,
treatment: str,
# ---------- forward-compat parameters ----------
aggregate: Optional[str] = None,
L_max: Optional[int] = None,
controls: Optional[List[str]] = None,
trends_linear: Optional[bool] = None,
trends_nonparam: Optional[Any] = None,
honest_did: bool = False,
# ---------- Phase 3 extensions ----------
heterogeneity: Optional[str] = None,
design2: bool = False,
# ---------- deferred (separate effort) ----------
survey_design: Any = None,
) -> ChaisemartinDHaultfoeuilleResults:
"""
Fit the dCDH estimator on individual-level panel data.
Parameters
----------
data : pd.DataFrame
Individual-level panel. Must contain columns for ``outcome``,
``group``, ``time``, and ``treatment``. The estimator
internally aggregates to ``(group, time)`` cells.
outcome : str
Outcome variable column name.
group : str
Group identifier column name. Treatment must be constant
within each ``(group, time)`` cell after aggregation;
``ValueError`` is raised if any cell has fractional
treatment after grouping (within-cell-varying treatment
indicates a fuzzy design not supported in Phase 1).
time : str
Time period column name. Must be sortable.
treatment : str
Per-observation treatment column. Must be numeric and constant
within each ``(group, time)`` cell. Both binary ``{0, 1}`` and
non-binary (ordinal or continuous) treatment are supported.
Non-binary treatment requires ``L_max >= 1``.
aggregate : str, optional
**Reserved for Phase 3.** Must be ``None``; any other value
raises ``NotImplementedError``.
L_max : int, optional
Maximum event-study horizon. When set, computes ``DID_l``
for ``l = 1, ..., L_max`` using the per-group building block
from Equation 3 of the dynamic companion paper. When
``None`` (default), only the ``l = 1`` contemporaneous-
switch estimator ``DID_M`` is computed (Phase 1 behavior).
Must be a positive integer not exceeding the number of
post-baseline periods in the panel.
controls : list of str, optional
Column names for covariate adjustment via residualization-style
``DID^X`` (Web Appendix Section 1.2). Requires ``L_max >= 1``.
One ``theta_hat`` per baseline treatment value, estimated by
OLS on not-yet-treated observations. NOT doubly-robust.
trends_linear : bool, optional
If ``True``, estimate group-specific linear trends via
``DID^{fd}`` (Web Appendix Section 1.3, Lemma 6). Requires
``L_max >= 1`` and at least 3 time periods.
trends_nonparam : str, optional
Column name for state-set membership. Restricts the control
pool to groups in the same set (Web Appendix Section 1.4).
Requires ``L_max >= 1`` and time-invariant values per group.
honest_did : bool, default=False
Run HonestDiD sensitivity analysis (Rambachan & Roth 2023) on
the placebo + event study surface. Requires ``L_max >= 1``.
Default: relative magnitudes (DeltaRM, Mbar=1.0), targeting
the equal-weight average over all post-treatment horizons
(``l_vec=None``). Results stored on
``results.honest_did_results``; ``None`` with a warning if
the solver fails. For custom parameters (e.g., targeting
the on-impact effect only via ``l_vec``), call
``compute_honest_did(results, ...)`` post-hoc instead.
heterogeneity : str, optional
Column name for a time-invariant covariate to test for
heterogeneous effects (Web Appendix Section 1.5, Lemma 7).
Partial implementation: post-treatment regressions only
(no placebo regressions or joint null test). Cannot be
combined with ``controls``, ``trends_linear``, or
``trends_nonparam``. Requires ``L_max >= 1``. Under
``by_path`` / ``paths_of_interest``, per-path
heterogeneity coefficients also surface on
``results.path_heterogeneity_effects`` and on
``to_dataframe(level="by_path")`` via ``het_*`` columns.
design2 : bool, default=False
If ``True``, identify and report switch-in/switch-out
(Design-2) groups. Convenience wrapper (descriptive summary,
not full paper re-estimation). Requires
``drop_larger_lower=False`` to retain 2-switch groups.
survey_design : SurveyDesign, optional
Survey design specification for design-based inference.
Supports ``weight_type='pweight'`` with two variance paths:
(1) Taylor Series Linearization using strata / PSU / FPC
(analytical) via the **cell-period IF allocator** that