diff-diff 2.5.0
A library for Difference-in-Differences causal inference analysis
Requires Python: <3.14,>=3.9
Keywords: causal-inference, difference-in-differences, econometrics, statistics, treatment-effects

# diff-diff

A Python library for Difference-in-Differences (DiD) causal inference analysis with an sklearn-like API and statsmodels-style outputs.

## Installation

```bash
pip install diff-diff
```

Or install from source:

```bash
git clone https://github.com/igerber/diff-diff.git
cd diff-diff
pip install -e .
```

## Quick Start

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

# Create sample data
data = pd.DataFrame({
    'outcome': [10, 11, 15, 18, 9, 10, 12, 13],
    'treated': [1, 1, 1, 1, 0, 0, 0, 0],
    'post': [0, 0, 1, 1, 0, 0, 1, 1]
})

# Fit the model
did = DifferenceInDifferences()
results = did.fit(data, outcome='outcome', treatment='treated', time='post')

# View results
print(results)  # DiDResults(ATT=3.0000, SE=1.7321, p=0.1583)
results.print_summary()
```

Output:
```
======================================================================
             Difference-in-Differences Estimation Results
======================================================================

Observations:                      8
Treated units:                     4
Control units:                     4
R-squared:                    0.9055

----------------------------------------------------------------------
Parameter           Estimate    Std. Err.     t-stat      P>|t|
----------------------------------------------------------------------
ATT                   3.0000       1.7321      1.732     0.1583
----------------------------------------------------------------------

95% Confidence Interval: [-1.8089, 7.8089]

Signif. codes: '***' 0.001, '**' 0.01, '*' 0.05, '.' 0.1
======================================================================
```

## Features

- **sklearn-like API**: Familiar `fit()` interface with `get_params()` and `set_params()`
- **Pythonic results**: Easy access to coefficients, standard errors, and confidence intervals
- **Multiple interfaces**: Column names or R-style formulas
- **Robust inference**: Heteroskedasticity-robust (HC1) and cluster-robust standard errors
- **Wild cluster bootstrap**: Valid inference with few clusters (<50) using Rademacher, Webb, or Mammen weights
- **Panel data support**: Two-way fixed effects estimator for panel designs
- **Multi-period analysis**: Event-study style DiD with period-specific treatment effects
- **Staggered adoption**: Callaway-Sant'Anna (2021), Sun-Abraham (2021), Borusyak-Jaravel-Spiess (2024) imputation, Two-Stage DiD (Gardner 2022), and Stacked DiD (Wing, Freedman & Hollingsworth 2024) estimators for heterogeneous treatment timing
- **Triple Difference (DDD)**: Ortiz-Villavicencio & Sant'Anna (2025) estimators with proper covariate handling
- **Synthetic DiD**: Combined DiD with synthetic control for improved robustness
- **Triply Robust Panel (TROP)**: Factor-adjusted DiD with synthetic weights (Athey et al. 2025)
- **Event study plots**: Publication-ready visualization of treatment effects
- **Parallel trends testing**: Multiple methods including equivalence tests
- **Goodman-Bacon decomposition**: Diagnose TWFE bias by decomposing into 2x2 comparisons
- **Placebo tests**: Comprehensive diagnostics including fake timing, fake group, permutation, and leave-one-out tests
- **Honest DiD sensitivity analysis**: Rambachan-Roth (2023) bounds and breakdown analysis for parallel trends violations
- **Pre-trends power analysis**: Roth (2022) minimum detectable violation (MDV) and power curves for pre-trends tests
- **Power analysis**: MDE, sample size, and power calculations for study design; simulation-based power for any estimator
- **Data prep utilities**: Helper functions for common data preparation tasks
- **Validated against R**: Benchmarked against `did`, `synthdid`, and `fixest` packages (see [benchmarks](docs/benchmarks.rst))

## Tutorials

We provide Jupyter notebook tutorials in `docs/tutorials/`:

| Notebook | Description |
|----------|-------------|
| `01_basic_did.ipynb` | Basic 2x2 DiD, formula interface, covariates, fixed effects, cluster-robust SE, wild bootstrap |
| `02_staggered_did.ipynb` | Staggered adoption with Callaway-Sant'Anna and Sun-Abraham, group-time effects, aggregation methods, Bacon decomposition |
| `03_synthetic_did.ipynb` | Synthetic DiD, unit/time weights, inference methods, regularization |
| `04_parallel_trends.ipynb` | Testing parallel trends, equivalence tests, placebo tests, diagnostics |
| `05_honest_did.ipynb` | Honest DiD sensitivity analysis, bounds, breakdown values, visualization |
| `06_power_analysis.ipynb` | Power analysis, MDE, sample size calculations, simulation-based power |
| `07_pretrends_power.ipynb` | Pre-trends power analysis (Roth 2022), MDV, power curves |
| `08_triple_diff.ipynb` | Triple Difference (DDD) estimation with proper covariate handling |
| `09_real_world_examples.ipynb` | Real-world data examples (Card-Krueger, Castle Doctrine, Divorce Laws) |
| `10_trop.ipynb` | Triply Robust Panel (TROP) estimation with factor model adjustment |

## Data Preparation

diff-diff provides utility functions to help prepare your data for DiD analysis. These functions handle common data transformation tasks like creating treatment indicators, reshaping panel data, and validating data formats.

### Generate Sample Data

Create synthetic data with a known treatment effect for testing and learning:

```python
from diff_diff import generate_did_data, DifferenceInDifferences

# Generate panel data with 100 units, 4 periods, and a treatment effect of 5
data = generate_did_data(
    n_units=100,
    n_periods=4,
    treatment_effect=5.0,
    treatment_fraction=0.5,  # 50% of units are treated
    treatment_period=2,       # Treatment starts at period 2
    seed=42
)

# Verify the estimator recovers the treatment effect
did = DifferenceInDifferences()
results = did.fit(data, outcome='outcome', treatment='treated', time='post')
print(f"Estimated ATT: {results.att:.2f} (true: 5.0)")
```

### Create Treatment Indicators

Convert categorical variables or numeric thresholds to binary treatment indicators:

```python
from diff_diff import make_treatment_indicator

# From categorical variable
df = make_treatment_indicator(
    data,
    column='state',
    treated_values=['CA', 'NY', 'TX']  # These states are treated
)

# From numeric threshold (e.g., firms above median size)
df = make_treatment_indicator(
    data,
    column='firm_size',
    threshold=data['firm_size'].median()
)

# Treat units below threshold
df = make_treatment_indicator(
    data,
    column='income',
    threshold=50000,
    above_threshold=False  # Units with income <= 50000 are treated
)
```

### Create Post-Treatment Indicators

Convert time/date columns to binary post-treatment indicators:

```python
from diff_diff import make_post_indicator

# From specific post-treatment periods
df = make_post_indicator(
    data,
    time_column='year',
    post_periods=[2020, 2021, 2022]
)

# From treatment start date
df = make_post_indicator(
    data,
    time_column='year',
    treatment_start=2020  # All years >= 2020 are post-treatment
)

# Works with datetime columns
df = make_post_indicator(
    data,
    time_column='date',
    treatment_start='2020-01-01'
)
```

### Reshape Wide to Long Format

Convert wide-format data (one row per unit, multiple time columns) to long format:

```python
from diff_diff import wide_to_long

# Wide format: columns like sales_2019, sales_2020, sales_2021
wide_df = pd.DataFrame({
    'firm_id': [1, 2, 3],
    'industry': ['tech', 'retail', 'tech'],
    'sales_2019': [100, 150, 200],
    'sales_2020': [110, 160, 210],
    'sales_2021': [120, 170, 220]
})

# Convert to long format for DiD
long_df = wide_to_long(
    wide_df,
    value_columns=['sales_2019', 'sales_2020', 'sales_2021'],
    id_column='firm_id',
    time_name='year',
    value_name='sales',
    time_values=[2019, 2020, 2021]
)
# Result: 9 rows (3 firms × 3 years), columns: firm_id, year, sales, industry
```

### Balance Panel Data

Ensure all units have observations for all time periods:

```python
from diff_diff import balance_panel

# Keep only units with complete data (drop incomplete units)
balanced = balance_panel(
    data,
    unit_column='firm_id',
    time_column='year',
    method='inner'
)

# Include all unit-period combinations (creates NaN for missing)
balanced = balance_panel(
    data,
    unit_column='firm_id',
    time_column='year',
    method='outer'
)

# Fill missing values
balanced = balance_panel(
    data,
    unit_column='firm_id',
    time_column='year',
    method='fill',
    fill_value=0  # Or None for forward/backward fill
)
```

### Validate Data

Check that your data meets DiD requirements before fitting:

```python
from diff_diff import validate_did_data

# Validate and get informative error messages
result = validate_did_data(
    data,
    outcome='sales',
    treatment='treated',
    time='post',
    unit='firm_id',      # Optional: for panel-specific validation
    raise_on_error=False  # Return dict instead of raising
)

if result['valid']:
    print("Data is ready for DiD analysis!")
    print(f"Summary: {result['summary']}")
else:
    print("Issues found:")
    for error in result['errors']:
        print(f"  - {error}")

for warning in result['warnings']:
    print(f"Warning: {warning}")
```

### Summarize Data by Groups

Get summary statistics for each treatment-time cell:

```python
from diff_diff import summarize_did_data

summary = summarize_did_data(
    data,
    outcome='sales',
    treatment='treated',
    time='post'
)
print(summary)
```

Output:
```
                        n      mean       std       min       max
Control - Pre        250  100.5000   15.2340   65.0000  145.0000
Control - Post       250  105.2000   16.1230   68.0000  152.0000
Treated - Pre        250  101.2000   14.8900   67.0000  143.0000
Treated - Post       250  115.8000   17.5600   72.0000  165.0000
DiD Estimate           -    9.9000         -         -         -
```

### Create Event Time for Staggered Designs

For designs where treatment occurs at different times:

```python
from diff_diff import create_event_time

# Add event-time column relative to treatment timing
df = create_event_time(
    data,
    time_column='year',
    treatment_time_column='treatment_year'
)
# Result: event_time = -2, -1, 0, 1, 2 relative to treatment
```

### Aggregate to Cohort Means

Aggregate unit-level data for visualization:

```python
from diff_diff import aggregate_to_cohorts

cohort_data = aggregate_to_cohorts(
    data,
    unit_column='firm_id',
    time_column='year',
    treatment_column='treated',
    outcome='sales'
)
# Result: mean outcome by treatment group and period
```

### Rank Control Units

Select the best control units for DiD or Synthetic DiD analysis by ranking them based on pre-treatment outcome similarity:

```python
from diff_diff import rank_control_units, generate_did_data

# Generate sample data
data = generate_did_data(n_units=50, n_periods=6, seed=42)

# Rank control units by their similarity to treated units
ranking = rank_control_units(
    data,
    unit_column='unit',
    time_column='period',
    outcome_column='outcome',
    treatment_column='treated',
    n_top=10  # Return top 10 controls
)

print(ranking[['unit', 'quality_score', 'pre_trend_rmse']])
```

Output:
```
   unit  quality_score  pre_trend_rmse
0    35         1.0000          0.4521
1    42         0.9234          0.5123
2    28         0.8876          0.5892
...
```

With covariates for matching:

```python
# Add covariate-based matching
ranking = rank_control_units(
    data,
    unit_column='unit',
    time_column='period',
    outcome_column='outcome',
    treatment_column='treated',
    covariates=['size', 'age'],  # Match on these too
    outcome_weight=0.7,          # 70% weight on outcome trends
    covariate_weight=0.3         # 30% weight on covariate similarity
)
```

Filter data for SyntheticDiD using top controls:

```python
from diff_diff import SyntheticDiD

# Get top control units
top_controls = ranking['unit'].tolist()

# Filter data to treated + top controls
filtered_data = data[
    (data['treated'] == 1) | (data['unit'].isin(top_controls))
]

# Fit SyntheticDiD with selected controls
sdid = SyntheticDiD()
results = sdid.fit(
    filtered_data,
    outcome='outcome',
    treatment='treated',
    unit='unit',
    time='period',
    post_periods=[3, 4, 5]
)
```

## Usage

### Basic DiD with Column Names

```python
from diff_diff import DifferenceInDifferences

did = DifferenceInDifferences(robust=True, alpha=0.05)
results = did.fit(
    data,
    outcome='sales',
    treatment='treated',
    time='post_policy'
)

# Access results
print(f"ATT: {results.att:.4f}")
print(f"Standard Error: {results.se:.4f}")
print(f"P-value: {results.p_value:.4f}")
print(f"95% CI: {results.conf_int}")
print(f"Significant: {results.is_significant}")
```

### Using Formula Interface

```python
# R-style formula syntax
results = did.fit(data, formula='outcome ~ treated * post')

# Explicit interaction syntax
results = did.fit(data, formula='outcome ~ treated + post + treated:post')

# With covariates
results = did.fit(data, formula='outcome ~ treated * post + age + income')
```

### Including Covariates

```python
results = did.fit(
    data,
    outcome='outcome',
    treatment='treated',
    time='post',
    covariates=['age', 'income', 'education']
)
```

### Fixed Effects

Use `fixed_effects` for low-dimensional categorical controls (creates dummy variables):

```python
# State and industry fixed effects
results = did.fit(
    data,
    outcome='sales',
    treatment='treated',
    time='post',
    fixed_effects=['state', 'industry']
)

# Access fixed effect coefficients
state_coefs = {k: v for k, v in results.coefficients.items() if k.startswith('state_')}
```

Use `absorb` for high-dimensional fixed effects (more efficient, uses within-transformation):

```python
# Absorb firm-level fixed effects (efficient for many firms)
results = did.fit(
    data,
    outcome='sales',
    treatment='treated',
    time='post',
    absorb=['firm_id']
)
```

Combine covariates with fixed effects:

```python
results = did.fit(
    data,
    outcome='sales',
    treatment='treated',
    time='post',
    covariates=['size', 'age'],           # Linear controls
    fixed_effects=['industry'],            # Low-dimensional FE (dummies)
    absorb=['firm_id']                     # High-dimensional FE (absorbed)
)
```

### Cluster-Robust Standard Errors

```python
did = DifferenceInDifferences(cluster='state')
results = did.fit(
    data,
    outcome='outcome',
    treatment='treated',
    time='post'
)
```

### Wild Cluster Bootstrap

When you have few clusters (<50), standard cluster-robust SEs are biased. Wild cluster bootstrap provides valid inference even with 5-10 clusters.

```python
# Use wild bootstrap for inference
did = DifferenceInDifferences(
    cluster='state',
    inference='wild_bootstrap',
    n_bootstrap=999,
    bootstrap_weights='rademacher',  # or 'webb' for <10 clusters, 'mammen'
    seed=42
)
results = did.fit(data, outcome='y', treatment='treated', time='post')

# Results include bootstrap-based SE and p-value
print(f"ATT: {results.att:.3f} (SE: {results.se:.3f})")
print(f"P-value: {results.p_value:.4f}")
print(f"95% CI: {results.conf_int}")
print(f"Inference method: {results.inference_method}")
print(f"Number of clusters: {results.n_clusters}")
```

**Weight types:**
- `'rademacher'` - Default, ±1 with p=0.5, good for most cases
- `'webb'` - 6-point distribution, recommended for <10 clusters
- `'mammen'` - Two-point distribution, alternative to Rademacher

Works with `DifferenceInDifferences` and `TwoWayFixedEffects` estimators.

### Two-Way Fixed Effects (Panel Data)

```python
from diff_diff import TwoWayFixedEffects

twfe = TwoWayFixedEffects()
results = twfe.fit(
    panel_data,
    outcome='outcome',
    treatment='treated',
    time='year',
    unit='firm_id'
)
```

### Multi-Period DiD (Event Study)

For settings with multiple pre- and post-treatment periods. Estimates treatment × period
interactions for ALL periods (pre and post), enabling parallel trends assessment:

```python
from diff_diff import MultiPeriodDiD

# Fit full event study with pre and post period effects
did = MultiPeriodDiD()
results = did.fit(
    panel_data,
    outcome='sales',
    treatment='treated',
    time='period',
    post_periods=[3, 4, 5],      # Periods 3-5 are post-treatment
    reference_period=2,          # Last pre-period (e=-1 convention)
    unit='unit_id',              # Optional: warns if staggered adoption detected
)

# Pre-period effects test parallel trends (should be ≈ 0)
for period, effect in results.pre_period_effects.items():
    print(f"Pre {period}: {effect.effect:.3f} (SE: {effect.se:.3f})")

# Post-period effects estimate dynamic treatment effects
for period, effect in results.post_period_effects.items():
    print(f"Post {period}: {effect.effect:.3f} (SE: {effect.se:.3f})")

# View average treatment effect across post-periods
print(f"Average ATT: {results.avg_att:.3f}")
print(f"Average SE: {results.avg_se:.3f}")

# Full summary with pre and post period effects
results.print_summary()
```

Output:
```
================================================================================
            Multi-Period Difference-in-Differences Estimation Results
================================================================================

Observations:                      600
Pre-treatment periods:             3
Post-treatment periods:            3

--------------------------------------------------------------------------------
Average Treatment Effect
--------------------------------------------------------------------------------
Average ATT       5.2000       0.8234      6.315      0.0000
--------------------------------------------------------------------------------
95% Confidence Interval: [3.5862, 6.8138]

Period-Specific Effects:
--------------------------------------------------------------------------------
Period            Effect     Std. Err.     t-stat      P>|t|
--------------------------------------------------------------------------------
3                 4.5000       0.9512      4.731      0.0000***
4                 5.2000       0.8876      5.858      0.0000***
5                 5.9000       0.9123      6.468      0.0000***
--------------------------------------------------------------------------------

Signif. codes: '***' 0.001, '**' 0.01, '*' 0.05, '.' 0.1
================================================================================
```

### Staggered Difference-in-Differences (Callaway-Sant'Anna)

When treatment is adopted at different times by different units, traditional TWFE estimators can be biased. The Callaway-Sant'Anna estimator provides unbiased estimates with staggered adoption.

```python
from diff_diff import CallawaySantAnna

# Panel data with staggered treatment
# 'first_treat' = period when unit was first treated (0 if never treated)
cs = CallawaySantAnna()
results = cs.fit(
    panel_data,
    outcome='sales',
    unit='firm_id',
    time='year',
    first_treat='first_treat',  # 0 for never-treated, else first treatment year
    aggregate='event_study'      # Compute event study effects
)

# View results
results.print_summary()

# Access group-time effects ATT(g,t)
for (group, time), effect in results.group_time_effects.items():
    print(f"Cohort {group}, Period {time}: {effect['effect']:.3f}")

# Event study effects (averaged by relative time)
for rel_time, effect in results.event_study_effects.items():
    print(f"e={rel_time}: {effect['effect']:.3f} (SE: {effect['se']:.3f})")

# Convert to DataFrame
df = results.to_dataframe(level='event_study')
```

Output:
```
=====================================================================================
          Callaway-Sant'Anna Staggered Difference-in-Differences Results
=====================================================================================

Total observations:                     600
Treated units:                           35
Control units:                           15
Treatment cohorts:                        3
Time periods:                             8
Control group:                never_treated

-------------------------------------------------------------------------------------
                  Overall Average Treatment Effect on the Treated
-------------------------------------------------------------------------------------
Parameter         Estimate     Std. Err.     t-stat      P>|t|   Sig.
-------------------------------------------------------------------------------------
ATT                 2.5000       0.3521       7.101     0.0000   ***
-------------------------------------------------------------------------------------

95% Confidence Interval: [1.8099, 3.1901]

-------------------------------------------------------------------------------------
                          Event Study (Dynamic) Effects
-------------------------------------------------------------------------------------
Rel. Period       Estimate     Std. Err.     t-stat      P>|t|   Sig.
-------------------------------------------------------------------------------------
0                   2.1000       0.4521       4.645     0.0000   ***
1                   2.5000       0.4123       6.064     0.0000   ***
2                   2.8000       0.5234       5.349     0.0000   ***
-------------------------------------------------------------------------------------

Signif. codes: '***' 0.001, '**' 0.01, '*' 0.05, '.' 0.1
=====================================================================================
```

**When to use Callaway-Sant'Anna vs TWFE:**

| Scenario | Use TWFE | Use Callaway-Sant'Anna |
|----------|----------|------------------------|
| All units treated at same time | ✓ | ✓ |
| Staggered adoption, homogeneous effects | ✓ | ✓ |
| Staggered adoption, heterogeneous effects | ✗ | ✓ |
| Need event study with staggered timing | ✗ | ✓ |
| Fewer than ~20 treated units | ✓ | Depends on design |

**Parameters:**

```python
CallawaySantAnna(
    control_group='never_treated',  # or 'not_yet_treated'
    anticipation=0,                  # Periods before treatment with effects
    estimation_method='dr',          # 'dr', 'ipw', or 'reg'
    alpha=0.05,                      # Significance level
    cluster=None,                    # Column for cluster SEs
    n_bootstrap=0,                   # Bootstrap iterations (0 = analytical SEs)
    bootstrap_weights='rademacher',  # 'rademacher', 'mammen', or 'webb'
    seed=None                        # Random seed
)
```

**Multiplier bootstrap for inference:**

With few clusters or when analytical standard errors may be unreliable, use the multiplier bootstrap for valid inference. This implements the approach from Callaway & Sant'Anna (2021).

```python
# Bootstrap inference with 999 iterations
cs = CallawaySantAnna(
    n_bootstrap=999,
    bootstrap_weights='rademacher',  # or 'mammen', 'webb'
    seed=42
)
results = cs.fit(
    data,
    outcome='sales',
    unit='firm_id',
    time='year',
    first_treat='first_treat',
    aggregate='event_study'
)

# Access bootstrap results
print(f"Overall ATT: {results.overall_att:.3f}")
print(f"Bootstrap SE: {results.bootstrap_results.overall_att_se:.3f}")
print(f"Bootstrap 95% CI: {results.bootstrap_results.overall_att_ci}")
print(f"Bootstrap p-value: {results.bootstrap_results.overall_att_p_value:.4f}")

# Event study bootstrap inference
for rel_time, se in results.bootstrap_results.event_study_ses.items():
    ci = results.bootstrap_results.event_study_cis[rel_time]
    print(f"e={rel_time}: SE={se:.3f}, 95% CI=[{ci[0]:.3f}, {ci[1]:.3f}]")
```

**Bootstrap weight types:**
- `'rademacher'` - Default, ±1 with p=0.5, good for most cases
- `'mammen'` - Two-point distribution matching first 3 moments
- `'webb'` - Six-point distribution, recommended for very few clusters (<10)

**Covariate adjustment for conditional parallel trends:**

When parallel trends only holds conditional on covariates, use the `covariates` parameter:

```python
# Doubly robust estimation with covariates
cs = CallawaySantAnna(estimation_method='dr')  # 'dr', 'ipw', or 'reg'
results = cs.fit(
    data,
    outcome='sales',
    unit='firm_id',
    time='year',
    first_treat='first_treat',
    covariates=['size', 'age', 'industry'],  # Covariates for conditional PT
    aggregate='event_study'
)
```

### Sun-Abraham Interaction-Weighted Estimator

The Sun-Abraham (2021) estimator provides an alternative to Callaway-Sant'Anna using an interaction-weighted (IW) regression approach. Running both estimators serves as a useful robustness check—when they agree, results are more credible.

```python
from diff_diff import SunAbraham

# Basic usage
sa = SunAbraham()
results = sa.fit(
    panel_data,
    outcome='sales',
    unit='firm_id',
    time='year',
    first_treat='first_treat'  # 0 for never-treated, else first treatment year
)

# View results
results.print_summary()

# Event study effects (by relative time to treatment)
for rel_time, effect in results.event_study_effects.items():
    print(f"e={rel_time}: {effect['effect']:.3f} (SE: {effect['se']:.3f})")

# Overall ATT
print(f"Overall ATT: {results.overall_att:.3f} (SE: {results.overall_se:.3f})")

# Cohort weights (how each cohort contributes to each event-time estimate)
for rel_time, weights in results.cohort_weights.items():
    print(f"e={rel_time}: {weights}")
```

**Parameters:**

```python
SunAbraham(
    control_group='never_treated',  # or 'not_yet_treated'
    anticipation=0,                  # Periods before treatment with effects
    alpha=0.05,                      # Significance level
    cluster=None,                    # Column for cluster SEs
    n_bootstrap=0,                   # Bootstrap iterations (0 = analytical SEs)
    bootstrap_weights='rademacher',  # 'rademacher', 'mammen', or 'webb'
    seed=None                        # Random seed
)
```

**Bootstrap inference:**

```python
# Bootstrap inference with 999 iterations
sa = SunAbraham(
    n_bootstrap=999,
    bootstrap_weights='rademacher',
    seed=42
)
results = sa.fit(
    data,
    outcome='sales',
    unit='firm_id',
    time='year',
    first_treat='first_treat'
)

# Access bootstrap results
print(f"Overall ATT: {results.overall_att:.3f}")
print(f"Bootstrap SE: {results.bootstrap_results.overall_att_se:.3f}")
print(f"Bootstrap 95% CI: {results.bootstrap_results.overall_att_ci}")
print(f"Bootstrap p-value: {results.bootstrap_results.overall_att_p_value:.4f}")
```

**When to use Sun-Abraham vs Callaway-Sant'Anna:**

| Aspect | Sun-Abraham | Callaway-Sant'Anna |
|--------|-------------|-------------------|
| Approach | Interaction-weighted regression | 2x2 DiD aggregation |
| Efficiency | More efficient under homogeneous effects | More robust to heterogeneity |
| Weighting | Weights by cohort share at each relative time | Weights by sample size |
| Use case | Robustness check, regression-based inference | Primary staggered DiD estimator |

**Both estimators should give similar results when:**
- Treatment effects are relatively homogeneous across cohorts
- Parallel trends holds

**Running both as robustness check:**

```python
from diff_diff import CallawaySantAnna, SunAbraham

# Callaway-Sant'Anna
cs = CallawaySantAnna()
cs_results = cs.fit(data, outcome='y', unit='unit', time='time', first_treat='first_treat')

# Sun-Abraham
sa = SunAbraham()
sa_results = sa.fit(data, outcome='y', unit='unit', time='time', first_treat='first_treat')

# Compare
print(f"Callaway-Sant'Anna ATT: {cs_results.overall_att:.3f}")
print(f"Sun-Abraham ATT: {sa_results.overall_att:.3f}")

# If results differ substantially, investigate heterogeneity
```

### Borusyak-Jaravel-Spiess Imputation Estimator

The Borusyak et al. (2024) imputation estimator is the **efficient** estimator for staggered DiD under parallel trends, producing ~50% shorter confidence intervals than Callaway-Sant'Anna and 2-3.5x shorter than Sun-Abraham under homogeneous treatment effects.

```python
from diff_diff import ImputationDiD, imputation_did

# Basic usage
est = ImputationDiD()
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat')
results.print_summary()

# Event study
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat',
                  aggregate='event_study')

# Pre-trend test (Equation 9)
pt = results.pretrend_test(n_leads=3)
print(f"F-stat: {pt['f_stat']:.3f}, p-value: {pt['p_value']:.4f}")

# Convenience function
results = imputation_did(data, 'outcome', 'unit', 'period', 'first_treat',
                         aggregate='all')
```

```python
ImputationDiD(
    anticipation=0,         # Number of anticipation periods
    alpha=0.05,             # Significance level
    cluster=None,           # Cluster variable (defaults to unit)
    n_bootstrap=0,          # Bootstrap iterations (0=analytical inference)
    seed=None,              # Random seed
    horizon_max=None,       # Max event-study horizon
    aux_partition="cohort_horizon",  # Variance partition: "cohort_horizon", "cohort", "horizon"
)
```

**When to use Imputation DiD vs Callaway-Sant'Anna:**

| Aspect | Imputation DiD | Callaway-Sant'Anna |
|--------|---------------|-------------------|
| Efficiency | Most efficient under homogeneous effects | Less efficient but more robust to heterogeneity |
| Control group | Always uses all untreated obs | Choice of never-treated or not-yet-treated |
| Inference | Conservative variance (Theorem 3) | Multiplier bootstrap |
| Pre-trends | Built-in F-test (Equation 9) | Separate testing |

### Two-Stage DiD (Gardner 2022)

Two-Stage DiD addresses TWFE bias in staggered adoption designs by estimating unit and time fixed effects on untreated observations only, then regressing the residualized outcomes on treatment indicators. Point estimates match the Imputation DiD estimator (Borusyak et al. 2024); the key difference is that Two-Stage DiD uses a GMM sandwich variance estimator that accounts for first-stage estimation error, while Imputation DiD uses a conservative variance (Theorem 3).

```python
from diff_diff import TwoStageDiD

# Basic usage
est = TwoStageDiD()
results = est.fit(data, outcome='outcome', unit='unit', time='period', first_treat='first_treat')
results.print_summary()
```

**Event study:**

```python
# Event study aggregation with visualization
results = est.fit(data, outcome='outcome', unit='unit', time='period',
                  first_treat='first_treat', aggregate='event_study')
plot_event_study(results)
```

**Parameters:**

```python
TwoStageDiD(
    anticipation=0,                   # Periods of anticipation effects
    alpha=0.05,                       # Significance level for CIs
    cluster=None,                     # Column for cluster-robust SEs (defaults to unit)
    n_bootstrap=0,                    # Bootstrap iterations (0 = analytical GMM SEs)
    seed=None,                        # Random seed
    rank_deficient_action='warn',     # 'warn', 'error', or 'silent'
    horizon_max=None,                 # Max event-study horizon
)
```

**When to use Two-Stage DiD vs Imputation DiD:**

| Aspect | Two-Stage DiD | Imputation DiD |
|--------|--------------|---------------|
| Point estimates | Identical | Identical |
| Variance | GMM sandwich (accounts for first-stage error) | Conservative (Theorem 3, may overcover) |
| Intuition | Residualize then regress | Impute counterfactuals then aggregate |
| Reference impl. | R `did2s` package | R `didimputation` package |

Both estimators are the efficient estimator under homogeneous treatment effects, producing shorter confidence intervals than Callaway-Sant'Anna or Sun-Abraham.

### Stacked DiD (Wing, Freedman & Hollingsworth 2024)

Stacked DiD addresses TWFE bias in staggered adoption settings by constructing a "clean" comparison dataset for each treatment cohort and stacking them together. Each cohort's sub-experiment compares units treated at that cohort's timing against units that are not yet treated (or never treated) within a symmetric event-study window. This avoids the "bad comparisons" problem in TWFE while retaining a regression-based framework that practitioners familiar with event studies will find intuitive.

```python
from diff_diff import StackedDiD, generate_staggered_data

# Generate sample data
data = generate_staggered_data(n_units=200, n_periods=12,
                                cohort_periods=[4, 6, 8], seed=42)

# Fit stacked DiD with event study
est = StackedDiD(kappa_pre=2, kappa_post=2)
results = est.fit(data, outcome='outcome', unit='unit',
                  time='period', first_treat='first_treat',
                  aggregate='event_study')
results.print_summary()

# Access stacked data for custom analysis
stacked = results.stacked_data

# Convenience function
from diff_diff import stacked_did
results = stacked_did(data, 'outcome', 'unit', 'period', 'first_treat',
                      kappa_pre=2, kappa_post=2, aggregate='event_study')
```

**Parameters:**

```python
StackedDiD(
    kappa_pre=1,                          # Pre-treatment event-study periods
    kappa_post=1,                         # Post-treatment event-study periods
    weighting='aggregate',                # 'aggregate', 'population', or 'sample_share'
    clean_control='not_yet_treated',      # 'not_yet_treated', 'strict', or 'never_treated'
    cluster='unit',                       # 'unit' or 'unit_subexp'
    alpha=0.05,                           # Significance level
    anticipation=0,                       # Anticipation periods
    rank_deficient_action='warn',         # 'warn', 'error', or 'silent'
)
```

> **Note:** Group aggregation (`aggregate='group'`) is not supported because the pooled
> stacked regression cannot produce cohort-specific effects. Use `CallawaySantAnna` or
> `ImputationDiD` for cohort-level estimates.

**When to use Stacked DiD vs Callaway-Sant'Anna:**

| Aspect | Stacked DiD | Callaway-Sant'Anna |
|--------|-------------|-------------------|
| Approach | Stack cohort sub-experiments, run pooled TWFE | 2x2 DiD aggregation |
| Symmetric windows | Enforced via kappa_pre / kappa_post | Not required |
| Control group | Not-yet-treated (default) or never-treated | Never-treated or not-yet-treated |
| Covariates | Passed to pooled regression | Doubly robust / IPW |
| Intuition | Familiar event-study regression | Nonparametric aggregation |

**Convenience function:**

```python
# One-liner estimation
results = stacked_did(
    data,
    outcome='outcome',
    unit='unit',
    time='period',
    first_treat='first_treat',
    kappa_pre=3,
    kappa_post=3,
    aggregate='event_study'
)
```

### Triple Difference (DDD)

Triple Difference (DDD) is used when treatment requires satisfying two criteria: belonging to a treated **group** AND being in an eligible **partition**. The `TripleDifference` class implements the methodology from Ortiz-Villavicencio & Sant'Anna (2025), which correctly handles covariate adjustment (unlike naive implementations).

```python
from diff_diff import TripleDifference, triple_difference

# Basic usage
ddd = TripleDifference(estimation_method='dr')  # doubly robust (recommended)
results = ddd.fit(
    data,
    outcome='wages',
    group='policy_state',       # 1=state enacted policy, 0=control state
    partition='female',         # 1=women (affected by policy), 0=men
    time='post'                 # 1=post-policy, 0=pre-policy
)

# View results
results.print_summary()
print(f"ATT: {results.att:.3f} (SE: {results.se:.3f})")

# With covariates (properly incorporated, unlike naive DDD)
results = ddd.fit(
    data,
    outcome='wages',
    group='policy_state',
    partition='female',
    time='post',
    covariates=['age', 'education', 'experience']
)
```

**Estimation methods:**

| Method | Description | When to use |
|--------|-------------|-------------|
| `"dr"` | Doubly robust | Recommended. Consistent if either outcome or propensity model is correct |
| `"reg"` | Regression adjustment | Simple outcome regression with full interactions |
| `"ipw"` | Inverse probability weighting | When propensity score model is well-specified |

```python
# Compare estimation methods
for method in ['reg', 'ipw', 'dr']:
    est = TripleDifference(estimation_method=method)
    res = est.fit(data, outcome='y', group='g', partition='p', time='t')
    print(f"{method}: ATT={res.att:.3f} (SE={res.se:.3f})")
```

**Convenience function:**

```python
# One-liner estimation
results = triple_difference(
    data,
    outcome='wages',
    group='policy_state',
    partition='female',
    time='post',
    covariates=['age', 'education'],
    estimation_method='dr'
)
```

**Why use DDD instead of DiD?**

DDD allows for violations of parallel trends that are:
- Group-specific (e.g., economic shocks in treatment states)
- Partition-specific (e.g., trends affecting women everywhere)

As long as these biases are additive, DDD differences them out. The key assumption is that the *differential* trend between eligible and ineligible units would be the same across groups.

### Event Study Visualization

Create publication-ready event study plots:

```python
from diff_diff import plot_event_study, MultiPeriodDiD, CallawaySantAnna, SunAbraham

# From MultiPeriodDiD (full event study with pre and post period effects)
did = MultiPeriodDiD()
results = did.fit(data, outcome='y', treatment='treated',
                  time='period', post_periods=[3, 4, 5], reference_period=2)
plot_event_study(results, title="Treatment Effects Over Time")

# From CallawaySantAnna (with event study aggregation)
cs = CallawaySantAnna()
results = cs.fit(data, outcome='y', unit='unit', time='period',
                 first_treat='first_treat', aggregate='event_study')
plot_event_study(results, title="Staggered DiD Event Study (CS)")

# From SunAbraham
sa = SunAbraham()
results = sa.fit(data, outcome='y', unit='unit', time='period',
                 first_treat='first_treat')
plot_event_study(results, title="Staggered DiD Event Study (SA)")

# From a DataFrame
df = pd.DataFrame({
    'period': [-2, -1, 0, 1, 2],
    'effect': [0.1, 0.05, 0.0, 2.5, 2.8],
    'se': [0.3, 0.25, 0.0, 0.4, 0.45]
})
plot_event_study(df, reference_period=0)

# With customization
ax = plot_event_study(
    results,
    title="Dynamic Treatment Effects",
    xlabel="Years Relative to Treatment",
    ylabel="Effect on Sales ($1000s)",
    color="#2563eb",
    marker="o",
    shade_pre=True,           # Shade pre-treatment region
    show_zero_line=True,      # Horizontal line at y=0
    show_reference_line=True, # Vertical line at reference period
    figsize=(10, 6),
    show=False                # Don't call plt.show(), return axes
)
```

### Synthetic Difference-in-Differences

Synthetic DiD combines the strengths of Difference-in-Differences and Synthetic Control methods by re-weighting control units to better match treated units' pre-treatment outcomes.

```python
from diff_diff import SyntheticDiD

# Fit Synthetic DiD model
sdid = SyntheticDiD()
results = sdid.fit(
    panel_data,
    outcome='gdp_growth',
    treatment='treated',
    unit='state',
    time='year',
    post_periods=[2015, 2016, 2017, 2018]
)

# View results
results.print_summary()
print(f"ATT: {results.att:.3f} (SE: {results.se:.3f})")

# Examine unit weights (which control units matter most)
weights_df = results.get_unit_weights_df()
print(weights_df.head(10))

# Examine time weights
time_weights_df = results.get_time_weights_df()
print(time_weights_df)
```

Output:
```
===========================================================================
         Synthetic Difference-in-Differences Estimation Results
=======================================================