-
Notifications
You must be signed in to change notification settings - Fork 19
Expand file tree
/
Copy pathimputation.py
More file actions
2009 lines (1767 loc) · 75.7 KB
/
imputation.py
File metadata and controls
2009 lines (1767 loc) · 75.7 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
"""
Borusyak-Jaravel-Spiess (2024) Imputation DiD Estimator.
Implements the efficient imputation estimator for staggered
Difference-in-Differences from Borusyak, Jaravel & Spiess (2024),
"Revisiting Event-Study Designs: Robust and Efficient Estimation",
Review of Economic Studies.
The estimator:
1. Runs OLS on untreated observations to estimate unit + time fixed effects
2. Imputes counterfactual Y(0) for treated observations
3. Aggregates imputed treatment effects with researcher-chosen weights
Inference uses the conservative clustered variance estimator (Theorem 3).
"""
import warnings
from typing import Any, Dict, List, Optional, Set, Tuple
import numpy as np
import pandas as pd
from scipy import sparse, stats
from scipy.sparse.linalg import spsolve
from diff_diff.imputation_bootstrap import ImputationDiDBootstrapMixin, _compute_target_weights
from diff_diff.imputation_results import ( # noqa: F401 (re-export)
ImputationBootstrapResults,
ImputationDiDResults,
)
from diff_diff.linalg import solve_ols
from diff_diff.utils import safe_inference
# =============================================================================
# Main Estimator
# =============================================================================
class ImputationDiD(ImputationDiDBootstrapMixin):
"""
Borusyak-Jaravel-Spiess (2024) imputation DiD estimator.
This is the efficient estimator for staggered Difference-in-Differences
under parallel trends. It produces shorter confidence intervals than
Callaway-Sant'Anna (~50% shorter) and Sun-Abraham (2-3.5x shorter)
under homogeneous treatment effects.
The estimation procedure:
1. Run OLS on untreated observations to estimate unit + time fixed effects
2. Impute counterfactual Y(0) for treated observations
3. Aggregate imputed treatment effects with researcher-chosen weights
Inference uses the conservative clustered variance estimator from Theorem 3
of the paper.
Parameters
----------
anticipation : int, default=0
Number of periods before treatment where effects may occur.
alpha : float, default=0.05
Significance level for confidence intervals.
cluster : str, optional
Column name for cluster-robust standard errors.
If None, clusters at the unit level by default.
n_bootstrap : int, default=0
Number of bootstrap iterations. If 0, uses analytical inference
(conservative variance from Theorem 3).
bootstrap_weights : str, default="rademacher"
Type of bootstrap weights: "rademacher", "mammen", or "webb".
seed : int, optional
Random seed for reproducibility.
rank_deficient_action : str, default="warn"
Action when design matrix is rank-deficient:
- "warn": Issue warning and drop linearly dependent columns
- "error": Raise ValueError
- "silent": Drop columns silently
horizon_max : int, optional
Maximum event-study horizon. If set, event study effects are only
computed for |h| <= horizon_max.
aux_partition : str, default="cohort_horizon"
Controls the auxiliary model partition for Theorem 3 variance:
- "cohort_horizon": Groups by cohort x relative time (tightest SEs)
- "cohort": Groups by cohort only (more conservative)
- "horizon": Groups by relative time only (more conservative)
Attributes
----------
results_ : ImputationDiDResults
Estimation results after calling fit().
is_fitted_ : bool
Whether the model has been fitted.
Examples
--------
Basic usage:
>>> from diff_diff import ImputationDiD, generate_staggered_data
>>> data = generate_staggered_data(n_units=200, seed=42)
>>> est = ImputationDiD()
>>> results = est.fit(data, outcome='outcome', unit='unit',
... time='time', first_treat='first_treat')
>>> results.print_summary()
With event study:
>>> est = ImputationDiD()
>>> results = est.fit(data, outcome='outcome', unit='unit',
... time='time', first_treat='first_treat',
... aggregate='event_study')
>>> from diff_diff import plot_event_study
>>> plot_event_study(results)
Notes
-----
The imputation estimator uses ALL untreated observations (never-treated +
not-yet-treated periods of eventually-treated units) to estimate the
counterfactual model. There is no ``control_group`` parameter because this
is fundamental to the method's efficiency.
References
----------
Borusyak, K., Jaravel, X., & Spiess, J. (2024). Revisiting Event-Study
Designs: Robust and Efficient Estimation. Review of Economic Studies,
91(6), 3253-3285.
"""
def __init__(
self,
anticipation: int = 0,
alpha: float = 0.05,
cluster: Optional[str] = None,
n_bootstrap: int = 0,
bootstrap_weights: str = "rademacher",
seed: Optional[int] = None,
rank_deficient_action: str = "warn",
horizon_max: Optional[int] = None,
aux_partition: str = "cohort_horizon",
):
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 aux_partition not in ("cohort_horizon", "cohort", "horizon"):
raise ValueError(
f"aux_partition must be 'cohort_horizon', 'cohort', or 'horizon', "
f"got '{aux_partition}'"
)
self.anticipation = anticipation
self.alpha = alpha
self.cluster = cluster
self.n_bootstrap = n_bootstrap
self.bootstrap_weights = bootstrap_weights
self.seed = seed
self.rank_deficient_action = rank_deficient_action
self.horizon_max = horizon_max
self.aux_partition = aux_partition
self.is_fitted_ = False
self.results_: Optional[ImputationDiDResults] = None
# Internal state preserved for pretrend_test()
self._fit_data: Optional[Dict[str, Any]] = None
def fit(
self,
data: pd.DataFrame,
outcome: str,
unit: str,
time: str,
first_treat: str,
covariates: Optional[List[str]] = None,
aggregate: Optional[str] = None,
balance_e: Optional[int] = None,
survey_design: object = None,
) -> ImputationDiDResults:
"""
Fit the imputation DiD estimator.
Parameters
----------
data : pd.DataFrame
Panel data with unit and time identifiers.
outcome : str
Name of outcome variable column.
unit : str
Name of unit identifier column.
time : str
Name of time period column.
first_treat : str
Name of column indicating when unit was first treated.
Use 0 (or np.inf) for never-treated units.
covariates : list of str, optional
List of covariate column names.
aggregate : str, optional
Aggregation mode: None/"simple" (overall ATT only),
"event_study", "group", or "all".
balance_e : int, optional
When computing event study, restrict to cohorts observed at all
relative times in [-balance_e, max_h].
survey_design : SurveyDesign, optional
Survey design specification for design-based inference. Supports
pweight only (aweight/fweight raise ValueError). FPC raises
NotImplementedError. PSU is used as cluster variable for Theorem 3
variance. Strata enters survey df for t-distribution inference.
Both analytical (n_bootstrap=0) and bootstrap inference are supported.
Returns
-------
ImputationDiDResults
Object containing all estimation results.
Raises
------
ValueError
If required columns are missing or data validation fails.
"""
# Validate inputs
required_cols = [outcome, unit, time, first_treat]
if covariates:
required_cols.extend(covariates)
missing = [c for c in required_cols if c not in data.columns]
if missing:
raise ValueError(f"Missing columns: {missing}")
# Create working copy
df = data.copy()
# Resolve survey design if provided
from diff_diff.survey import (
_inject_cluster_as_psu,
_resolve_effective_cluster,
_resolve_survey_for_fit,
_validate_unit_constant_survey,
)
resolved_survey, survey_weights, _, survey_metadata = _resolve_survey_for_fit(
survey_design, data, "analytical"
)
# Validate within-unit constancy for panel survey designs
if resolved_survey is not None:
if resolved_survey.uses_replicate_variance:
raise NotImplementedError(
"ImputationDiD does not yet support replicate-weight survey "
"designs. Use a TSL-based survey design (strata/psu/fpc)."
)
_validate_unit_constant_survey(data, unit, survey_design)
if resolved_survey.weight_type != "pweight":
raise ValueError(
f"ImputationDiD survey support requires weight_type='pweight', "
f"got '{resolved_survey.weight_type}'. The survey variance math "
f"assumes probability weights (pweight)."
)
if resolved_survey.fpc is not None:
raise NotImplementedError(
"ImputationDiD does not yet support FPC (finite population "
"correction) in SurveyDesign. Weights, strata (for survey df), "
"and PSU (for cluster-robust variance) are supported."
)
# Bootstrap + survey supported via PSU-level multiplier bootstrap.
# Ensure numeric types
df[time] = pd.to_numeric(df[time])
df[first_treat] = pd.to_numeric(df[first_treat])
# Validate absorbing treatment: first_treat must be constant within each unit
ft_nunique = df.groupby(unit)[first_treat].nunique()
non_constant = ft_nunique[ft_nunique > 1]
if len(non_constant) > 0:
example_unit = non_constant.index[0]
example_vals = sorted(df.loc[df[unit] == example_unit, first_treat].unique())
warnings.warn(
f"{len(non_constant)} unit(s) have non-constant '{first_treat}' "
f"values (e.g., unit '{example_unit}' has values {example_vals}). "
f"ImputationDiD assumes treatment is an absorbing state "
f"(once treated, always treated) with a single treatment onset "
f"time per unit. Non-constant first_treat violates this assumption "
f"and may produce unreliable estimates.",
UserWarning,
stacklevel=2,
)
# Coerce to per-unit value so downstream code
# (_never_treated, _treated, _rel_time) uses a single
# consistent first_treat per unit.
df[first_treat] = df.groupby(unit)[first_treat].transform("first")
# Identify treatment status
df["_never_treated"] = (df[first_treat] == 0) | (df[first_treat] == np.inf)
# Check for always-treated units (treated in all observed periods)
min_time = df[time].min()
always_treated_mask = (~df["_never_treated"]) & (df[first_treat] <= min_time)
n_always_treated = df.loc[always_treated_mask, unit].nunique()
if n_always_treated > 0:
warnings.warn(
f"{n_always_treated} unit(s) are treated in all observed periods "
f"(first_treat <= {min_time}). These units have no untreated "
"observations and cannot contribute to the counterfactual model. "
"Their treatment effects will be imputed but may be unreliable.",
UserWarning,
stacklevel=2,
)
# Create treatment indicator D_it
# D_it = 1 if t >= first_treat and first_treat > 0
# With anticipation: D_it = 1 if t >= first_treat - anticipation
effective_treat = df[first_treat] - self.anticipation
df["_treated"] = (~df["_never_treated"]) & (df[time] >= effective_treat)
# Identify Omega_0 (untreated) and Omega_1 (treated)
omega_0_mask = ~df["_treated"]
omega_1_mask = df["_treated"]
n_omega_0 = int(omega_0_mask.sum())
n_omega_1 = int(omega_1_mask.sum())
if n_omega_0 == 0:
raise ValueError(
"No untreated observations found. Cannot estimate counterfactual model."
)
if n_omega_1 == 0:
raise ValueError("No treated observations found. Nothing to estimate.")
# Identify groups and time periods
time_periods = sorted(df[time].unique())
treatment_groups = sorted([g for g in df[first_treat].unique() if g > 0 and g != np.inf])
if len(treatment_groups) == 0:
raise ValueError("No treated units found. Check 'first_treat' column.")
# Unit info
unit_info = (
df.groupby(unit).agg({first_treat: "first", "_never_treated": "first"}).reset_index()
)
n_treated_units = int((~unit_info["_never_treated"]).sum())
# Control units = units with at least one untreated observation
units_in_omega_0 = df.loc[omega_0_mask, unit].unique()
n_control_units = len(units_in_omega_0)
# Cluster variable
cluster_var = self.cluster if self.cluster is not None else unit
if self.cluster is not None and self.cluster not in df.columns:
raise ValueError(
f"Cluster column '{self.cluster}' not found in data. "
f"Available columns: {list(df.columns)}"
)
# Resolve effective cluster and inject cluster-as-PSU for survey variance
if resolved_survey is not None:
cluster_ids_raw = df[cluster_var].values if cluster_var in df.columns else None
effective_cluster_ids = _resolve_effective_cluster(
resolved_survey,
cluster_ids_raw,
cluster_var if self.cluster is not None else None,
)
resolved_survey = _inject_cluster_as_psu(resolved_survey, effective_cluster_ids)
# When survey PSU is present, use it as the effective cluster for
# Theorem 3 variance (PSU overrides unit-level clustering)
if resolved_survey.psu is not None:
# Create a temporary column with PSU IDs for cluster_var
df["_survey_cluster"] = resolved_survey.psu
cluster_var = "_survey_cluster"
# Recompute metadata after PSU injection
if resolved_survey.psu is not None and survey_metadata is not None:
from diff_diff.survey import compute_survey_metadata
raw_w = (
data[survey_design.weights].values.astype(np.float64)
if survey_design.weights
else np.ones(len(data), dtype=np.float64)
)
survey_metadata = compute_survey_metadata(resolved_survey, raw_w)
# Compute relative time
df["_rel_time"] = np.where(
~df["_never_treated"],
df[time] - df[first_treat],
np.nan,
)
# ---- Step 1: OLS on untreated observations ----
unit_fe, time_fe, grand_mean, delta_hat, kept_cov_mask = self._fit_untreated_model(
df, outcome, unit, time, covariates, omega_0_mask, weights=survey_weights
)
# ---- Rank condition checks ----
# Check: every treated unit should have >= 1 untreated period (for unit FE)
treated_unit_ids = df.loc[omega_1_mask, unit].unique()
units_with_fe = set(unit_fe.keys())
units_missing_fe = set(treated_unit_ids) - units_with_fe
# Check: every post-treatment period should have >= 1 untreated unit (for time FE)
post_period_ids = df.loc[omega_1_mask, time].unique()
periods_with_fe = set(time_fe.keys())
periods_missing_fe = set(post_period_ids) - periods_with_fe
if units_missing_fe or periods_missing_fe:
parts = []
if units_missing_fe:
sorted_missing = sorted(units_missing_fe)
parts.append(
f"{len(units_missing_fe)} treated unit(s) have no untreated "
f"periods (units: {sorted_missing[:5]}"
f"{'...' if len(units_missing_fe) > 5 else ''})"
)
if periods_missing_fe:
sorted_missing = sorted(periods_missing_fe)
parts.append(
f"{len(periods_missing_fe)} post-treatment period(s) have no "
f"untreated units (periods: {sorted_missing[:5]}"
f"{'...' if len(periods_missing_fe) > 5 else ''})"
)
msg = (
"Rank condition violated: "
+ "; ".join(parts)
+ ". Affected treatment effects will be NaN."
)
if self.rank_deficient_action == "error":
raise ValueError(msg)
elif self.rank_deficient_action == "warn":
warnings.warn(msg, UserWarning, stacklevel=2)
# "silent": continue without warning
# ---- Step 2: Impute treatment effects ----
tau_hat, y_hat_0 = self._impute_treatment_effects(
df,
outcome,
unit,
time,
covariates,
omega_1_mask,
unit_fe,
time_fe,
grand_mean,
delta_hat,
)
# Store tau_hat in dataframe
df["_tau_hat"] = np.nan
df.loc[omega_1_mask, "_tau_hat"] = tau_hat
# ---- Step 3: Aggregate ----
# Always compute overall ATT (simple aggregation)
finite_mask = np.isfinite(tau_hat)
valid_tau = tau_hat[finite_mask]
if len(valid_tau) == 0:
overall_att = np.nan
elif survey_weights is not None:
# Survey-weighted ATT: use treated obs' survey weights
treated_survey_w = survey_weights[omega_1_mask.values]
w_finite = treated_survey_w[finite_mask]
overall_att = float(np.average(valid_tau, weights=w_finite))
else:
overall_att = float(np.mean(valid_tau))
# ---- Conservative variance (Theorem 3) ----
# Build weights matching the ATT: proportional to survey weights for
# finite tau_hat, uniform when no survey
overall_weights = np.zeros(n_omega_1)
n_valid = int(finite_mask.sum())
if n_valid > 0:
if survey_weights is not None:
treated_sw = survey_weights[omega_1_mask.values]
sw_finite = treated_sw[finite_mask]
overall_weights[finite_mask] = sw_finite / sw_finite.sum()
else:
overall_weights[finite_mask] = 1.0 / n_valid
if n_valid == 0:
overall_se = np.nan
else:
overall_se = self._compute_conservative_variance(
df=df,
outcome=outcome,
unit=unit,
time=time,
first_treat=first_treat,
covariates=covariates,
omega_0_mask=omega_0_mask,
omega_1_mask=omega_1_mask,
unit_fe=unit_fe,
time_fe=time_fe,
grand_mean=grand_mean,
delta_hat=delta_hat,
weights=overall_weights,
cluster_var=cluster_var,
kept_cov_mask=kept_cov_mask,
survey_weights=survey_weights,
)
# Survey degrees of freedom for t-distribution inference
_survey_df = resolved_survey.df_survey if resolved_survey is not None else None
overall_t, overall_p, overall_ci = safe_inference(
overall_att, overall_se, alpha=self.alpha, df=_survey_df
)
# Event study and group aggregation
event_study_effects = None
group_effects = None
if aggregate in ("event_study", "all"):
event_study_effects = self._aggregate_event_study(
df=df,
outcome=outcome,
unit=unit,
time=time,
first_treat=first_treat,
covariates=covariates,
omega_0_mask=omega_0_mask,
omega_1_mask=omega_1_mask,
unit_fe=unit_fe,
time_fe=time_fe,
grand_mean=grand_mean,
delta_hat=delta_hat,
cluster_var=cluster_var,
treatment_groups=treatment_groups,
balance_e=balance_e,
kept_cov_mask=kept_cov_mask,
survey_weights=survey_weights,
survey_df=_survey_df,
)
if aggregate in ("group", "all"):
group_effects = self._aggregate_group(
df=df,
outcome=outcome,
unit=unit,
time=time,
first_treat=first_treat,
covariates=covariates,
omega_0_mask=omega_0_mask,
omega_1_mask=omega_1_mask,
unit_fe=unit_fe,
time_fe=time_fe,
grand_mean=grand_mean,
delta_hat=delta_hat,
cluster_var=cluster_var,
treatment_groups=treatment_groups,
kept_cov_mask=kept_cov_mask,
survey_weights=survey_weights,
survey_df=_survey_df,
)
# Build treatment effects dataframe
treated_df = df.loc[omega_1_mask, [unit, time, "_tau_hat", "_rel_time"]].copy()
treated_df = treated_df.rename(columns={"_tau_hat": "tau_hat", "_rel_time": "rel_time"})
# Weights consistent with actual ATT: zero for NaN tau_hat
tau_finite = treated_df["tau_hat"].notna()
n_valid_te = int(tau_finite.sum())
if n_valid_te > 0:
if survey_weights is not None:
# Survey-weighted: use normalized survey weights for treated obs
treated_sw = survey_weights[omega_1_mask.values]
sw_finite = np.where(tau_finite, treated_sw, 0.0)
sw_sum = sw_finite.sum()
treated_df["weight"] = sw_finite / sw_sum if sw_sum > 0 else 0.0
else:
treated_df["weight"] = np.where(tau_finite, 1.0 / n_valid_te, 0.0)
else:
treated_df["weight"] = 0.0
# Store fit data for pretrend_test
self._fit_data = {
"df": df,
"outcome": outcome,
"unit": unit,
"time": time,
"first_treat": first_treat,
"covariates": covariates,
"omega_0_mask": omega_0_mask,
"omega_1_mask": omega_1_mask,
"cluster_var": cluster_var,
"unit_fe": unit_fe,
"time_fe": time_fe,
"grand_mean": grand_mean,
"delta_hat": delta_hat,
"kept_cov_mask": kept_cov_mask,
"survey_design": survey_design,
}
# Pre-compute cluster psi sums for bootstrap
psi_data = None
if self.n_bootstrap > 0 and n_valid > 0:
try:
# Extract survey weights for untreated obs (same as analytical path)
_sw_0 = survey_weights[omega_0_mask.values] if survey_weights is not None else None
# Extract survey weights for treated obs (event-study/group bootstrap paths)
_sw_1 = survey_weights[omega_1_mask.values] if survey_weights is not None else None
psi_data = self._precompute_bootstrap_psi(
df=df,
outcome=outcome,
unit=unit,
time=time,
first_treat=first_treat,
covariates=covariates,
omega_0_mask=omega_0_mask,
omega_1_mask=omega_1_mask,
unit_fe=unit_fe,
time_fe=time_fe,
grand_mean=grand_mean,
delta_hat=delta_hat,
cluster_var=cluster_var,
kept_cov_mask=kept_cov_mask,
overall_weights=overall_weights,
event_study_effects=event_study_effects,
group_effects=group_effects,
treatment_groups=treatment_groups,
tau_hat=tau_hat,
balance_e=balance_e,
survey_weights_0=_sw_0,
survey_weights_1=_sw_1,
)
except Exception as e:
warnings.warn(
f"Bootstrap pre-computation failed: {e}. " "Skipping bootstrap inference.",
UserWarning,
stacklevel=2,
)
psi_data = None
# Bootstrap
bootstrap_results = None
if self.n_bootstrap > 0 and psi_data is not None:
bootstrap_results = self._run_bootstrap(
original_att=overall_att,
original_event_study=event_study_effects,
original_group=group_effects,
psi_data=psi_data,
resolved_survey=resolved_survey,
)
# Update inference with bootstrap results
overall_se = bootstrap_results.overall_att_se
overall_t = (
overall_att / overall_se if np.isfinite(overall_se) and overall_se > 0 else np.nan
)
overall_p = bootstrap_results.overall_att_p_value
overall_ci = bootstrap_results.overall_att_ci
# Update event study
if event_study_effects and bootstrap_results.event_study_ses:
for h in event_study_effects:
if (
h in bootstrap_results.event_study_ses
and event_study_effects[h].get("n_obs", 1) > 0
):
event_study_effects[h]["se"] = bootstrap_results.event_study_ses[h]
assert bootstrap_results.event_study_cis is not None
event_study_effects[h]["conf_int"] = bootstrap_results.event_study_cis[h]
assert bootstrap_results.event_study_p_values is not None
event_study_effects[h]["p_value"] = bootstrap_results.event_study_p_values[
h
]
eff_val = event_study_effects[h]["effect"]
se_val = event_study_effects[h]["se"]
event_study_effects[h]["t_stat"] = safe_inference(
eff_val, se_val, alpha=self.alpha
)[0]
# Update group effects
if group_effects and bootstrap_results.group_ses:
for g in group_effects:
if g in bootstrap_results.group_ses:
group_effects[g]["se"] = bootstrap_results.group_ses[g]
assert bootstrap_results.group_cis is not None
group_effects[g]["conf_int"] = bootstrap_results.group_cis[g]
assert bootstrap_results.group_p_values is not None
group_effects[g]["p_value"] = bootstrap_results.group_p_values[g]
eff_val = group_effects[g]["effect"]
se_val = group_effects[g]["se"]
group_effects[g]["t_stat"] = safe_inference(
eff_val, se_val, alpha=self.alpha
)[0]
# Construct results
self.results_ = ImputationDiDResults(
treatment_effects=treated_df,
overall_att=overall_att,
overall_se=overall_se,
overall_t_stat=overall_t,
overall_p_value=overall_p,
overall_conf_int=overall_ci,
event_study_effects=event_study_effects,
group_effects=group_effects,
groups=treatment_groups,
time_periods=time_periods,
n_obs=len(df),
n_treated_obs=n_omega_1,
n_untreated_obs=n_omega_0,
n_treated_units=n_treated_units,
n_control_units=n_control_units,
alpha=self.alpha,
bootstrap_results=bootstrap_results,
_estimator_ref=self,
survey_metadata=survey_metadata,
)
self.is_fitted_ = True
return self.results_
# =========================================================================
# Step 1: OLS on untreated observations
# =========================================================================
def _iterative_fe(
self,
y: np.ndarray,
unit_vals: np.ndarray,
time_vals: np.ndarray,
idx: pd.Index,
max_iter: int = 100,
tol: float = 1e-10,
weights: Optional[np.ndarray] = None,
) -> Tuple[Dict[Any, float], Dict[Any, float]]:
"""
Estimate unit and time FE via iterative alternating projection (Gauss-Seidel).
Converges to the exact OLS solution for both balanced and unbalanced panels.
For balanced panels, converges in 1-2 iterations (identical to one-pass).
For unbalanced panels, typically 5-20 iterations.
Parameters
----------
weights : np.ndarray, optional
Survey weights. When provided, uses weighted group means
(sum(w*x)/sum(w)) instead of unweighted means.
Returns
-------
unit_fe : dict
Mapping from unit -> unit fixed effect.
time_fe : dict
Mapping from time -> time fixed effect.
"""
n = len(y)
alpha = np.zeros(n) # unit FE broadcast to obs level
beta = np.zeros(n) # time FE broadcast to obs level
# Precompute per-group weight sums (invariant across iterations)
if weights is not None:
w_series = pd.Series(weights, index=idx)
wsum_t = w_series.groupby(time_vals).transform("sum").values
wsum_u = w_series.groupby(unit_vals).transform("sum").values
with np.errstate(invalid="ignore", divide="ignore"):
for iteration in range(max_iter):
resid_after_alpha = y - alpha
if weights is not None:
wr_t = pd.Series(resid_after_alpha * weights, index=idx)
beta_new = wr_t.groupby(time_vals).transform("sum").values / wsum_t
else:
beta_new = (
pd.Series(resid_after_alpha, index=idx)
.groupby(time_vals)
.transform("mean")
.values
)
resid_after_beta = y - beta_new
if weights is not None:
wr_u = pd.Series(resid_after_beta * weights, index=idx)
alpha_new = wr_u.groupby(unit_vals).transform("sum").values / wsum_u
else:
alpha_new = (
pd.Series(resid_after_beta, index=idx)
.groupby(unit_vals)
.transform("mean")
.values
)
# Check convergence on FE changes
max_change = max(
np.max(np.abs(alpha_new - alpha)),
np.max(np.abs(beta_new - beta)),
)
alpha = alpha_new
beta = beta_new
if max_change < tol:
break
unit_fe = pd.Series(alpha, index=idx).groupby(unit_vals).first().to_dict()
time_fe = pd.Series(beta, index=idx).groupby(time_vals).first().to_dict()
return unit_fe, time_fe
@staticmethod
def _iterative_demean(
vals: np.ndarray,
unit_vals: np.ndarray,
time_vals: np.ndarray,
idx: pd.Index,
max_iter: int = 100,
tol: float = 1e-10,
weights: Optional[np.ndarray] = None,
) -> np.ndarray:
"""Demean a vector by iterative alternating projection (unit + time FE removal).
Converges to the exact within-transformation for both balanced and
unbalanced panels. For balanced panels, converges in 1-2 iterations.
Parameters
----------
weights : np.ndarray, optional
Survey weights. When provided, uses weighted group means
(sum(w*x)/sum(w)) instead of unweighted means.
"""
result = vals.copy()
# Precompute per-group weight sums (invariant across iterations)
if weights is not None:
w_series = pd.Series(weights, index=idx)
wsum_t = w_series.groupby(time_vals).transform("sum").values
wsum_u = w_series.groupby(unit_vals).transform("sum").values
with np.errstate(invalid="ignore", divide="ignore"):
for _ in range(max_iter):
if weights is not None:
wr_t = pd.Series(result * weights, index=idx)
time_means = wr_t.groupby(time_vals).transform("sum").values / wsum_t
else:
time_means = (
pd.Series(result, index=idx).groupby(time_vals).transform("mean").values
)
result_after_time = result - time_means
if weights is not None:
wr_u = pd.Series(result_after_time * weights, index=idx)
unit_means = wr_u.groupby(unit_vals).transform("sum").values / wsum_u
else:
unit_means = (
pd.Series(result_after_time, index=idx)
.groupby(unit_vals)
.transform("mean")
.values
)
result_new = result_after_time - unit_means
if np.max(np.abs(result_new - result)) < tol:
result = result_new
break
result = result_new
return result
@staticmethod
def _compute_balanced_cohort_mask(
df_treated: pd.DataFrame,
first_treat: str,
all_horizons: List[int],
balance_e: int,
cohort_rel_times: Dict[Any, Set[int]],
) -> np.ndarray:
"""Compute boolean mask selecting treated obs from balanced cohorts.
A cohort is 'balanced' if it has observations at every relative time
in [-balance_e, max(all_horizons)].
Parameters
----------
df_treated : pd.DataFrame
Post-treatment observations (Omega_1).
first_treat : str
Column name for cohort identifier.
all_horizons : list of int
Post-treatment horizons in the event study.
balance_e : int
Number of pre-treatment periods to require.
cohort_rel_times : dict
Maps each cohort value to the set of all observed relative times
(including pre-treatment) from the full panel. Built by
_build_cohort_rel_times().
"""
if not all_horizons:
return np.ones(len(df_treated), dtype=bool)
max_h = max(all_horizons)
required_range = set(range(-balance_e, max_h + 1))
balanced_cohorts = set()
for g, horizons in cohort_rel_times.items():
if required_range.issubset(horizons):
balanced_cohorts.add(g)
return df_treated[first_treat].isin(balanced_cohorts).values
@staticmethod
def _build_cohort_rel_times(
df: pd.DataFrame,
first_treat: str,
) -> Dict[Any, Set[int]]:
"""Build mapping of cohort -> set of observed relative times from full panel.
Precondition: df must have '_never_treated' and '_rel_time' columns
(set by fit() before any aggregation calls).
"""
treated_mask = ~df["_never_treated"]
treated_df = df.loc[treated_mask]
result: Dict[Any, Set[int]] = {}
ft_vals = treated_df[first_treat].values
rt_vals = treated_df["_rel_time"].values
for i in range(len(treated_df)):
h = rt_vals[i]
if np.isfinite(h):
result.setdefault(ft_vals[i], set()).add(int(h))
return result
def _fit_untreated_model(
self,
df: pd.DataFrame,
outcome: str,
unit: str,
time: str,
covariates: Optional[List[str]],
omega_0_mask: pd.Series,
weights: Optional[np.ndarray] = None,
) -> Tuple[
Dict[Any, float], Dict[Any, float], float, Optional[np.ndarray], Optional[np.ndarray]
]:
"""
Step 1: Estimate unit + time FE on untreated observations.
Uses iterative alternating projection (Gauss-Seidel) to compute exact
OLS fixed effects for both balanced and unbalanced panels. For balanced
panels, converges in 1-2 iterations (identical to one-pass demeaning).
Parameters
----------
weights : np.ndarray, optional
Full-panel survey weights (same length as df). The untreated subset
is extracted internally via omega_0_mask. When None, unweighted.
Returns
-------
unit_fe : dict
Unit fixed effects {unit_id: alpha_i}.
time_fe : dict
Time fixed effects {time_period: beta_t}.
grand_mean : float
Grand mean (0.0 — absorbed into iterative FE).
delta_hat : np.ndarray or None
Covariate coefficients (if covariates provided).
kept_cov_mask : np.ndarray or None
Boolean mask of shape (n_covariates,) indicating which covariates
have finite coefficients. None if no covariates.
"""
df_0 = df.loc[omega_0_mask]
w_0 = weights[omega_0_mask.values] if weights is not None else None
if covariates is None or len(covariates) == 0:
# No covariates: estimate FE via iterative alternating projection
# (exact OLS for both balanced and unbalanced panels)
y = df_0[outcome].values.copy()
unit_fe, time_fe = self._iterative_fe(
y, df_0[unit].values, df_0[time].values, df_0.index, weights=w_0
)
# grand_mean = 0: iterative FE absorb the intercept
return unit_fe, time_fe, 0.0, None, None
else:
# With covariates: iteratively demean Y and X, OLS for delta,
# then recover FE from covariate-adjusted outcome
y = df_0[outcome].values.copy()
X_raw = df_0[covariates].values.copy()
units = df_0[unit].values
times = df_0[time].values
n_cov = len(covariates)
# Step A: Iteratively demean Y and all X columns to remove unit+time FE
y_dm = self._iterative_demean(y, units, times, df_0.index, weights=w_0)
X_dm = np.column_stack(
[
self._iterative_demean(X_raw[:, j], units, times, df_0.index, weights=w_0)
for j in range(n_cov)
]
)
# Step B: OLS for covariate coefficients on demeaned data
result = solve_ols(
X_dm,
y_dm,
return_vcov=False,
rank_deficient_action=self.rank_deficient_action,
column_names=covariates,
weights=w_0,
)
delta_hat = result[0]
# Mask of covariates with finite coefficients (before cleaning)
# Used to exclude rank-deficient covariates from variance design matrices
kept_cov_mask = np.isfinite(delta_hat)
# Replace NaN coefficients with 0 for adjustment
# (rank-deficient covariates are dropped)