Add Event Study (aka Dynamic Difference in Differences) functionality

AGENTS.md

@@ -36,6 +36,12 @@ See the [python-environment skill](.agents/skills/python-environment/SKILL.md) f
 
 ## Documentation
 
+- **Dual README files**: The project has two files that must be kept in sync:
+  - `README.md` (GitHub landing page)
+  - `docs/source/index.md` (documentation website homepage)
+
+  When adding major new features to the Features table or making other content changes, update **both files**. The Features table lists all quasi-experimental methods supported by CausalPy.
+- **Reporting statistics**: When adding new experiment types, update the "Experiment Support" table in `docs/source/knowledgebase/reporting_statistics.md` to document `effect_summary()` support status for the new experiment.
 - **Structure**: Notebooks (how-to examples) go in `docs/source/notebooks/`, knowledgebase (educational content) goes in `docs/source/knowledgebase/`
 - **Notebook naming**: Use lowercase hyphen-separated words with the full method name spelled out, then the dataset/variant token (if any), then the backend (`pymc` or `sklearn`). Pattern: `{method}[-{variant}]-{backend}.ipynb` (e.g., `difference-in-differences-pymc.ipynb`, `regression-discontinuity-drinking-sklearn.ipynb`, `synthetic-control-brexit-pymc.ipynb`). Prefer descriptive words over acronyms for SEO: `interrupted-time-series` over `its`, `regression-discontinuity` over `rd`, `synthetic-control` over `sc`, `instrumental-variables` over `iv`. When renaming an existing notebook, update `docs/source/notebooks/gallery.yaml`, run `make gallery`, add an entry to `rediraffe_redirects` in `docs/source/conf.py`, and never remove older redirect keys.
 - **MyST directives**: Use `:::{note}` and other MyST features for callouts and formatting

ARCHITECTURE.md

@@ -55,6 +55,7 @@ Instantiation fits eagerly in `__init__`: `_build_design_matrices()` → `_prepa
 | `InterruptedTimeSeries` | ITS | OLS + Bayes | 3-period design via `treatment_end_time` |
 | `PiecewiseITS` | Segmented ITS | OLS + Bayes | Fits full series; `step()`/`ramp()` transforms |
 | `DifferenceInDifferences` | DiD | OLS + Bayes | Effect from interaction coefficient |
+| `EventStudy` | Event Study / Dynamic DiD | OLS + Bayes | Simultaneous treatment only; event-time dummies |
 | `StaggeredDifferenceInDifferences` | Staggered DiD | OLS + Bayes | Fits untreated obs only |
 | `SyntheticControl` | SC | OLS + Bayes | Multi-unit; control/treated unit lists, no formula |
 | `SyntheticDifferenceInDifferences` | SDiD | OLS + Bayes | Tau computed analytically from weight posteriors |

CONTRIBUTING.md

@@ -122,6 +122,8 @@ For more instructions see the [Pull request checklist](#pull-request-checklist)
 
     To update an existing environment after changes to `environment.yml` (e.g. after pulling or after regenerating it from `pyproject.toml`):
 
+    **Note:** Use `python -m pip` instead of just `pip` to ensure you're using the conda environment's pip, not a system pip.
+
     ```bash
     conda env update --file environment.yml --prune
     ```

README.md

@@ -133,6 +133,7 @@ CausalPy provides methods for common causal inference decision contexts:
 | Geographical lift | Measures the impact of an intervention in a specific geographic area by comparing it to similar areas without the intervention. Commonly used in marketing to assess regional campaigns. |
 | ANCOVA | Analysis of Covariance combines ANOVA and regression to control for the effects of one or more quantitative covariates. Used when comparing group means while controlling for other variables. |
 | Differences in Differences | Compares the changes in outcomes over time between a treatment group and a control group. Used in observational studies to estimate causal effects by accounting for time trends. |
+| Event Study | Estimates dynamic treatment effects over event time (time relative to treatment). Extends difference-in-differences by estimating separate effects for each time period, enabling pre-trend validation and analysis of how causal effects evolve. |
 | Staggered Difference-in-Differences | Estimates event-time treatment effects when different units adopt treatment at different times, using an imputation approach that models untreated outcomes and compares observed outcomes to counterfactual predictions. |
 | Regression discontinuity | Identifies causal effects by exploiting a cutoff or threshold in an assignment variable. Used when treatment is assigned based on a threshold value of an observed variable, allowing comparison just above and below the cutoff. |
 | Regression kink designs | Focuses on changes in the slope (kinks) of the relationship between variables rather than jumps at cutoff points. Used to identify causal effects when treatment intensity changes at a threshold. |

causalpy/__init__.py

@@ -24,6 +24,7 @@
 
 from .data import load_data
 from .experiments.diff_in_diff import DifferenceInDifferences
+from .experiments.event_study import EventStudy
 from .experiments.instrumental_variable import InstrumentalVariable
 from .experiments.interrupted_time_series import InterruptedTimeSeries
 from .experiments.inverse_propensity_weighting import InversePropensityWeighting
@@ -52,6 +53,7 @@
     "create_causalpy_compatible_class",
     "DifferenceInDifferences",
     "EstimateEffect",
+    "EventStudy",
     "extract_lift_for_mmm",
     "GenerateReport",
     "InstrumentalVariable",

causalpy/data/simulate_data.py

@@ -76,7 +76,7 @@ def generate_synthetic_control_data(
     :param seed:
         Random seed for reproducibility
 
-    Example
+    Examples
     --------
     >>> from causalpy.data.simulate_data import generate_synthetic_control_data
     >>> df, weightings_true = generate_synthetic_control_data(
@@ -211,7 +211,7 @@ def generate_did(seed: int | None = None) -> pd.DataFrame:
     :param seed:
         Random seed for reproducibility
 
-    Example
+    Examples
     --------
     >>> from causalpy.data.simulate_data import generate_did
     >>> df = generate_did(seed=42)
@@ -277,7 +277,7 @@ def generate_regression_discontinuity_data(
     :param seed:
         Random seed for reproducibility
 
-    Example
+    Examples
     --------
     >>> import pathlib
     >>> from causalpy.data.simulate_data import generate_regression_discontinuity_data
@@ -316,7 +316,7 @@ def generate_ancova_data(
     :param seed:
         Random seed for reproducibility
 
-    Example
+    Examples
     --------
     >>> import pathlib
     >>> from causalpy.data.simulate_data import generate_ancova_data
@@ -482,11 +482,229 @@ def generate_multicell_geolift_data(
     return df
 
 
+def generate_event_study_data(
+    n_units: int = 20,
+    n_time: int = 20,
+    treatment_time: int = 10,
+    treated_fraction: float = 0.5,
+    event_window: tuple[int, int] = (-5, 5),
+    treatment_effects: dict[int, float] | None = None,
+    unit_fe_sigma: float = 1.0,
+    time_fe_sigma: float = 0.5,
+    noise_sigma: float = 0.2,
+    predictor_effects: dict[str, float] | None = None,
+    ar_phi: float = 0.9,
+    ar_scale: float = 1.0,
+    seed: int | None = None,
+) -> pd.DataFrame:
+    """
+    Generate synthetic panel data for event study / dynamic DiD analysis.
+
+    Creates panel data with unit and time fixed effects, where a fraction of units
+    receive treatment at a common treatment time. Treatment effects can vary by
+    event time (time relative to treatment). Optionally includes time-varying
+    predictor variables generated via AR(1) processes.
+
+    Parameters
+    ----------
+    n_units : int
+        Total number of units (treated + control). Default 20.
+    n_time : int
+        Number of time periods. Default 20.
+    treatment_time : int
+        Time period when treatment occurs (0-indexed). Default 10.
+    treated_fraction : float
+        Fraction of units that are treated. Default 0.5.
+    event_window : tuple[int, int]
+        Range of event times (K_min, K_max) for which treatment effects are defined.
+        Default (-5, 5).
+    treatment_effects : dict[int, float], optional
+        Dictionary mapping event time k to treatment effect beta_k.
+        Default creates effects that are 0 for k < 0 (pre-treatment)
+        and gradually increase post-treatment.
+    unit_fe_sigma : float
+        Standard deviation for unit fixed effects. Default 1.0.
+    time_fe_sigma : float
+        Standard deviation for time fixed effects. Default 0.5.
+    noise_sigma : float
+        Standard deviation for observation noise. Default 0.2.
+    predictor_effects : dict[str, float], optional
+        Dictionary mapping predictor names to their true coefficients.
+        Each predictor is generated as an AR(1) time series that varies over time
+        but is the same for all units at a given time. For example,
+        ``{'temperature': 0.3, 'humidity': -0.1}`` creates two predictors.
+        Default None (no predictors).
+    ar_phi : float
+        AR(1) autoregressive coefficient controlling persistence of predictors.
+        Values closer to 1 produce smoother, more persistent series.
+        Default 0.9.
+    ar_scale : float
+        Standard deviation of the AR(1) innovation noise for predictors.
+        Default 1.0.
+    seed : int, optional
+        Random seed for reproducibility.
+
+    Returns
+    -------
+    pd.DataFrame
+        Panel data with columns:
+
+        - unit: Unit identifier
+        - time: Time period
+        - y: Outcome variable
+        - treat_time: Treatment time for unit (NaN if never treated)
+        - treated: Whether unit is in treated group (0 or 1)
+        - <predictor_name>: One column per predictor (if predictor_effects provided)
+
+    Examples
+    --------
+    >>> from causalpy.data.simulate_data import generate_event_study_data
+    >>> df = generate_event_study_data(
+    ...     n_units=20, n_time=20, treatment_time=10, seed=42
+    ... )
+    >>> df.shape
+    (400, 5)
+    >>> df.columns.tolist()
+    ['unit', 'time', 'y', 'treat_time', 'treated']
+
+    With predictors:
+
+    >>> df = generate_event_study_data(
+    ...     n_units=10,
+    ...     n_time=10,
+    ...     treatment_time=5,
+    ...     seed=42,
+    ...     predictor_effects={"temperature": 0.3, "humidity": -0.1},
+    ... )
+    >>> df.shape
+    (100, 7)
+    >>> "temperature" in df.columns and "humidity" in df.columns
+    True
+    """
+    if seed is not None:
+        np.random.seed(seed)
+
+    # Default treatment effects: zero pre-treatment, gradual increase post-treatment
+    if treatment_effects is None:
+        treatment_effects = {}
+        for k in range(event_window[0], event_window[1] + 1):
+            if k < 0:
+                treatment_effects[k] = 0.0  # No anticipation
+            else:
+                # Gradual treatment effect that increases post-treatment
+                treatment_effects[k] = 0.5 + 0.1 * k
+
+    # Determine treated units
+    n_treated = int(n_units * treated_fraction)
+    treated_units = set(range(n_treated))
+
+    # Generate unit fixed effects
+    unit_fe = np.random.normal(0, unit_fe_sigma, n_units)
+
+    # Generate time fixed effects
+    time_fe = np.random.normal(0, time_fe_sigma, n_time)
+
+    # Generate predictor time series (if any)
+    # Each predictor is an AR(1) series that varies over time but is the same
+    # for all units at a given time
+    predictors: dict[str, np.ndarray] = {}
+    if predictor_effects is not None:
+        for predictor_name in predictor_effects:
+            predictors[predictor_name] = generate_ar1_series(
+                n=n_time, phi=ar_phi, scale=ar_scale
+            )
+
+    # Build panel data
+    data = []
+    for unit in range(n_units):
+        is_treated = unit in treated_units
+        unit_treat_time = treatment_time if is_treated else np.nan
+
+        for t in range(n_time):
+            # Base outcome: unit FE + time FE + noise
+            y = unit_fe[unit] + time_fe[t] + np.random.normal(0, noise_sigma)
+
+            # Add predictor contributions to outcome
+            if predictor_effects is not None:
+                for predictor_name, coef in predictor_effects.items():
+                    y += coef * predictors[predictor_name][t]
+
+            # Add treatment effect for treated units in event window
+            if is_treated:
+                event_time = t - treatment_time
+                if (
+                    event_window[0] <= event_time <= event_window[1]
+                    and event_time in treatment_effects
+                ):
+                    y += treatment_effects[event_time]
+
+            row = {
+                "unit": unit,
+                "time": t,
+                "y": y,
+                "treat_time": unit_treat_time,
+                "treated": 1 if is_treated else 0,
+            }
+            # Add predictor values to the row
+            for predictor_name, series in predictors.items():
+                row[predictor_name] = series[t]
+
+            data.append(row)
+
+    return pd.DataFrame(data)
+
+
 # -----------------
 # UTILITY FUNCTIONS
 # -----------------
 
 
+def generate_ar1_series(
+    n: int,
+    phi: float = 0.9,
+    scale: float = 1.0,
+    initial: float = 0.0,
+) -> np.ndarray:
+    """
+    Generate an AR(1) autoregressive time series.
+
+    The AR(1) process is defined as:
+        x_{t+1} = phi * x_t + eta_t, where eta_t ~ N(0, scale^2)
+
+    Parameters
+    ----------
+    n : int
+        Length of the time series to generate.
+    phi : float
+        Autoregressive coefficient controlling persistence. Values closer to 1
+        produce smoother, more persistent series. Must be in (-1, 1) for
+        stationarity. Default 0.9.
+    scale : float
+        Standard deviation of the innovation noise. Default 1.0.
+    initial : float
+        Initial value of the series. Default 0.0.
+
+    Returns
+    -------
+    np.ndarray
+        Array of length n containing the AR(1) time series.
+
+    Example
+    -------
+    >>> from causalpy.data.simulate_data import generate_ar1_series
+    >>> np.random.seed(42)
+    >>> series = generate_ar1_series(n=10, phi=0.9, scale=0.5)
+    >>> len(series)
+    10
+    """
+    series = np.zeros(n)
+    series[0] = initial
+    innovations = np.random.normal(0, scale, n - 1)
+    for t in range(1, n):
+        series[t] = phi * series[t - 1] + innovations[t - 1]
+    return series
+
+
 def _generate_seasonality(
     n: int,
     amplitude: int,
@@ -509,6 +727,18 @@ def _generate_seasonality(
     return rng.multivariate_normal(np.zeros(n), cov)
 
 
+def generate_seasonality(
+    n: int = 12, amplitude: int = 1, length_scale: float = 0.5
+) -> np.ndarray:
+    """Generate monthly seasonality by sampling from a Gaussian process."""
+    return _generate_seasonality(
+        n=n,
+        amplitude=amplitude,
+        length_scale=length_scale,
+        rng=np.random.default_rng(),
+    )
+
+
 def _periodic_kernel(
     x1: np.ndarray,
     x2: np.ndarray,

causalpy/experiments/__init__.py

@@ -14,6 +14,7 @@
 """CausalPy experiment module."""
 
 from .diff_in_diff import DifferenceInDifferences
+from .event_study import EventStudy
 from .instrumental_variable import InstrumentalVariable
 from .interrupted_time_series import InterruptedTimeSeries
 from .inverse_propensity_weighting import InversePropensityWeighting
@@ -28,6 +29,7 @@
 
 __all__ = [
     "DifferenceInDifferences",
+    "EventStudy",
     "InstrumentalVariable",
     "InterruptedTimeSeries",
     "InversePropensityWeighting",

causalpy/experiments/base.py

@@ -449,7 +449,7 @@ def effect_summary(
         Parameters
         ----------
         window : str, tuple, or slice, default="post"
-            Time window for analysis (ITS/SC only, ignored for DiD/RD):
+            Time window for analysis (ITS/SC only, ignored for DiD/RD/EventStudy):
 
             - "post": All post-treatment time points (default)
             - (start, end): Tuple of start and end times (handles both datetime and integer indices)
@@ -470,10 +470,10 @@ def effect_summary(
             a 95% HDI in :meth:`effect_summary` and a 94% HDI in :meth:`plot`
             with default settings.
         cumulative : bool, default=True
-            Whether to include cumulative effect statistics (ITS/SC only, ignored for DiD/RD)
+            Whether to include cumulative effect statistics (ITS/SC only, ignored for DiD/RD/EventStudy)
         relative : bool, default=True
             Whether to include relative effect statistics (% change vs counterfactual)
-            (ITS/SC only, ignored for DiD/RD)
+            (ITS/SC only, ignored for DiD/RD/EventStudy)
         min_effect : float, optional
             Region of Practical Equivalence (ROPE) threshold (PyMC only, ignored for OLS).
             If provided, reports ``P(|effect| > min_effect)`` for two-sided or

causalpy/experiments/diff_in_diff.py

@@ -138,9 +138,11 @@ def __init__(
 
     def _build_design_matrices(self) -> None:
         """Build design matrices from formula and data using patsy."""
-        y, X = build_formula_matrices(self.formula, self.data)
+        y, X = build_formula_matrices(self.formula, self.data, return_type="dataframe")
         self._y_design_info = y.design_info
         self._x_design_info = X.design_info
+        # Filter data to rows that patsy kept (in case NaN values were dropped)
+        self.data = self.data.loc[X.index]
         self.labels = X.design_info.column_names
         self._y_raw, self._X_raw = np.asarray(y), np.asarray(X)
         self.outcome_variable_name = y.design_info.column_names[0]