EventStudyInteracts.py

EventStudyInteracts.py estimates the Sun and Abraham interaction-weighted event-study model in Python. The regressions are run by pyfixest, so the package keeps pyfixest's high-dimensional fixed effects, clustered covariance matrices, demeaning backends, solvers, weights, and speed.

This repository is maintained with Codex-assisted development. The estimator logic is intentionally narrow: eventreg() parses the user formula, builds the cohort-by-relative-time saturated regression, lets pyfixest.feols() do the high-dimensional fixed-effect estimation and covariance calculation, computes the IW count-share weights, and then forms the IW coefficient vector and delta-method covariance matrix. To guard correctness, the test suite checks direct pyfixest/TWFE equivalence in single-cohort designs, pyfixest saturated parity in the id + time case, manual binning and heterogeneity behavior, ATE/lincom covariance formulas, multi-FE plus continuous-control examples, and external parity against Stata eventstudyinteract and Julia EventStudyInteracts.jl when those tools are available.

Agent Skill

A usage-oriented skill for Codex, Claude Code, and similar coding agents is included at skills/eventstudyinteracts/SKILL.md. The skill tells an agent how to prepare data, write eventreg() formulas, use controls and multiple fixed effects, create binning and heterogeneity columns, plot results, and compute ATE/lincom.

To install it for Codex in this local setup:

mkdir -p ~/.codex/skills
cp -R skills/eventstudyinteracts ~/.codex/skills/

To install it for Claude Code as a user skill:

mkdir -p ~/.claude/skills
cp -R skills/eventstudyinteracts ~/.claude/skills/

For a project-local Claude Code skill, copy the same folder to .claude/skills/eventstudyinteracts/ inside the project.

Why This Package Exists

pyfixest already has a fast saturated event-study estimator:

import pyfixest as pf

fit = pf.event_study(
    data,
    yname="y",
    idname="unit",
    tname="time",
    gname="first_treated_period",
    estimator="saturated",
)

That is the right tool when the model is the built-in unit and time fixed-effect event-study. EventStudyInteracts.py is for the cases where the Stata eventstudyinteract and Julia EventStudyInteracts.jl interface is more natural:

• relative-time dummy columns are written directly in the formula, for example ln_wage ~ g_18 + g_17 + ... + g18 | idcode + year;
• continuous controls are needed, for example covariates="south + tenure + age";
• absorbed fixed effects go beyond unit and time, for example idcode + year + ind_code^year or idcode + province^year;
• binning and heterogeneity are created manually as columns and then estimated;
• ate() and lincom() are computed from the IW covariance matrix.

Internally, eventreg() calls pyfixest.feols() for the saturated cohort-by-relative-time regression and for the share covariance regressions. The wrapper adds the Stata/Julia IW aggregation logic and result methods.

Installation

The package is not on PyPI yet, so installing by package name is not available until the first PyPI release.

Current GitHub install for users:

python -m pip install "git+https://github.com/xiaobaaaa/EventStudyInteracts.py.git"

This also installs pyfixest. fit.iplot_aggregate() uses the same Matplotlib-based plotting style as pyfixest's saturated event-study object, so there is no separate plotting extra to install.

After a PyPI release, the intended command is:

python -m pip install eventstudyinteracts

Requires Python 3.10+ and pyfixest>=0.50.1.

Basic nlswork Example

This follows the nlswork example from the Stata help file. The first example uses only idcode and year fixed effects so it can be compared directly with pyfixest.event_study(..., estimator="saturated"), Stata eventstudyinteract, and Julia EventStudyInteracts.jl.

import numpy as np
import pandas as pd
import eventstudyinteracts as esi

# Load the same NLS worker data used by the Stata help example.
df = pd.read_stata("dataset/nlswork.dta")

# first_union is the first calendar year in which a worker is unionized.
df["union_year"] = np.where(df["union"].fillna(0).eq(1), df["year"], np.nan)
df["first_union"] = df.groupby("idcode")["union_year"].transform("min")
df = df.drop(columns=["union_year"])

# ry is event time: calendar year minus first union year.
df["ry"] = df["year"] - df["first_union"]

# Controls are never-union workers. eventreg uses this indicator as the
# authoritative control definition, so first_union may be missing for controls.
df["never_union"] = df["first_union"].isna()

# For package comparisons, keep treated workers observed in the omitted period.
# The omitted period is ry == -1, matching ref = -1 in saturated event studies.
ids_with_ref = df.loc[df["ry"].eq(-1), "idcode"].unique()
df = df.loc[df["never_union"] | df["idcode"].isin(ids_with_ref)].copy()

# Create all observed relative-time columns except the omitted period g_1.
relmin = int(abs(np.nanmin(df["ry"])))
relmax = int(abs(np.nanmax(df["ry"])))

event_vars = []
for k in range(relmin, 1, -1):
    name = f"g_{k}"
    df[name] = df["ry"].eq(-k).fillna(False).astype(int)
    event_vars.append(name)

for k in range(0, relmax + 1):
    name = f"g{k}"
    df[name] = df["ry"].eq(k).fillna(False).astype(int)
    event_vars.append(name)

fit = esi.eventreg(
    f"ln_wage ~ {' + '.join(event_vars)} | idcode + year",
    data=df,
    cohort="first_union",
    control_cohort="never_union",
    vcov={"CRV1": "idcode"},
    demeaner_backend="numba",
)

eventreg() already returns the IW relative-time estimates. Use tidy() for the event-study table:

fit.tidy()
fit.coef()
fit.se()
fit.vcov()
fit.summary()

fit.aggregate() is kept only as a pyfixest-compatible alias. It is not a required second estimation step. In this package:

fit.tidy()       # main event-study table
fit.aggregate()  # same IW estimates, with pyfixest-style aggregate() naming

Plotting

For a pyfixest-style event-study plot:

fit.iplot_aggregate()

To save the plot:

ax = fit.iplot_aggregate(show=False)
ax.figure.savefig("nlswork_eventstudy.png", dpi=200, bbox_inches="tight")

Manual plotting is also straightforward:

import matplotlib.pyplot as plt

plot_df = fit.aggregate(period_index=True)
err_low = plot_df["Estimate"] - plot_df["2.5%"]
err_high = plot_df["97.5%"] - plot_df["Estimate"]

fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(
    plot_df.index,
    plot_df["Estimate"],
    yerr=[err_low, err_high],
    fmt="o-",
    capsize=3,
)
ax.axhline(0, color="black", linewidth=1, linestyle="--")
ax.axvline(-1, color="gray", linewidth=1, linestyle=":")
ax.set_xlabel("Relative time")
ax.set_ylabel("IW estimate")
fig.tight_layout()

Basic Comparison Results

The comparison sample has 20,318 rows before singleton removal, 19,813 rows used by the fixed-effect regressions, and 3,541 workers. The omitted period is g_1. The first table reports selected coefficients. The second table reports the corresponding cluster-idcode standard errors.

term	eventreg	pyfixest saturated	Stata eventstudyinteract	Julia EventStudyInteracts.jl
g_20	0.111521	0.111521	0.111521	0.111521
g_19	0.118273	0.118273	0.118273	0.118273
g_18	-0.065599	-0.065599	-0.065599	-0.065599
g_17	0.045408	0.045408	0.045408	0.045408
g_16	-0.023468	-0.023468	-0.023468	-0.023468
g_15	-0.208510	-0.208510	-0.208510	-0.208510
g_14	-0.239402	-0.239402	-0.239402	-0.239402
g_3	-0.024659	-0.024659	-0.024659	-0.024659
g_2	0.001476	0.001476	0.001476	0.001476
g0	0.131608	0.131608	0.131608	0.131608
g1	0.177707	0.177707	0.177707	0.177707
g2	0.132955	0.132955	0.132955	0.132955
g3	0.149590	0.149590	0.149590	0.149590
g4	0.173574	0.173574	0.173574	0.173574
g18	-0.120038	-0.120038	-0.120038	-0.120038

term	eventreg SE	pyfixest saturated SE	Stata eventstudyinteract SE	Julia EventStudyInteracts.jl SE
g_20	0.117627	0.117627	0.117627	0.117630
g_19	0.123875	0.123875	0.123875	0.123878
g_18	0.143241	0.143241	0.143241	0.143245
g_17	0.091193	0.091193	0.091193	0.091196
g_16	0.134788	0.134788	0.134788	0.134791
g_15	0.098071	0.097923	0.098069	0.098073
g_14	0.218520	0.218520	0.218520	0.218526
g_3	0.023623	0.023185	0.023618	0.023624
g_2	0.020685	0.020575	0.020683	0.020685
g0	0.014452	0.013996	0.014447	0.014452
g1	0.021273	0.020693	0.021267	0.021273
g2	0.023436	0.023104	0.023432	0.023436
g3	0.030111	0.029412	0.030104	0.030112
g4	0.028911	0.028690	0.028909	0.028912
g18	0.067475	0.067475	0.067474	0.067476

The exact comparison code is in scripts/check_external_parity.py. Current maximum absolute differences against eventreg are:

comparison	coefficient	standard error	reason
pyfixest saturated, share_vcov="none"	<= 1e-12	<= 1e-12	same fixed-share covariance convention
pyfixest saturated, default share_vcov="regression"	0.0	0.003700	default adds share-step covariance
Stata eventstudyinteract	1.12e-08	0.000038	external numerical/output tolerance
Julia EventStudyInteracts.jl	2.52e-08	0.000006	external numerical/output tolerance

Point estimates use count-share IW weights. For valid relative-time columns that are 0/1 and mutually exclusive within each row, count-share weights and share-regression weights are the same mathematical weights; the regression form can introduce tiny floating-point differences. Therefore this package uses count shares for point estimates.

There are two inference conventions:

fit = esi.eventreg(..., share_vcov="regression")  # default
fit_fixed_share = esi.eventreg(..., share_vcov="none")

• share_vcov="regression" adds the share-step covariance term delta' V_share delta, estimated using pyfixest share regressions. This is the default because it follows the IW variance calculation used by the Stata/Julia-style estimator.
• share_vcov="none" treats the shares as fixed. This matches pyfixest saturated aggregate() exactly in the unit/time fixed-effect comparison.

Controls And Multiple Fixed Effects

Controls are passed through covariates=. Fixed effects are written after | with pyfixest/fixest syntax.

fit_fe = esi.eventreg(
    f"ln_wage ~ {' + '.join(event_vars)} | idcode + year + ind_code^year",
    data=df,
    cohort="first_union",
    control_cohort="never_union",
    covariates="south + tenure + age",
    vcov={"CRV1": "idcode"},
    demeaner_backend="rust-cg",
)

ind_code^year is an absorbed industry-by-year fixed effect. A province-by-year fixed effect is written the same way, for example province^year.

The deterministic multi-FE test uses:

esi.eventreg(
    "y ~ rel0 + rel1 | unit + time + market",
    data=df,
    cohort="cohort",
    control_cohort="control",
    covariates="x",
    vcov="iid",
)

and checks the result against Stata eventstudyinteract and Julia EventStudyInteracts.jl.

Binning

Binning is manual. Create the exact binned column you want, then put that column in the formula.

# One lead bin for all periods <= -4.
df["g_l4"] = df["ry"].le(-4).fillna(False).astype(int)

bin_vars = ["g_l4", "g_3", "g_2"]
bin_vars.extend([f"g{k}" for k in range(0, 19)])

fit_bin = esi.eventreg(
    f"ln_wage ~ {' + '.join(bin_vars)} | idcode + year",
    data=df,
    cohort="first_union",
    control_cohort="never_union",
    covariates="south",
    vcov={"CRV1": "idcode"},
)

The output terms are exactly the columns in bin_vars.

Last-Treated Controls

If the last-treated cohort is the control cohort, mark it and restrict the sample to pre-treatment periods for that last-treated cohort.

df["last_union"] = df["first_union"].eq(88).fillna(False)
sample = df.loc[df["first_union"].notna() & df["year"].lt(88)].copy()

fit_last = esi.eventreg(
    f"ln_wage ~ {' + '.join(bin_vars)} | idcode + year",
    data=sample,
    cohort="first_union",
    control_cohort="last_union",
    covariates="south",
    vcov={"CRV1": "idcode"},
)

For eventreg(), the control_cohort indicator defines the controls. The cohort value for controls may be missing, 0, or the last-treated cohort value, as long as control_cohort is correct.

Binning And Heterogeneity

Binning combines periods. Heterogeneity splits periods by group. They are separate operations, and both are created as ordinary 0/1 columns.

het_vars = []

for group in [0, 1]:
    # Binned pre-period for each college-graduate group.
    name = f"g_l4_collgrad{group}"
    df[name] = (
        df["ry"].le(-4).fillna(False) & df["collgrad"].eq(group)
    ).astype(int)
    het_vars.append(name)

    # Separate leads -3 and -2 for each group.
    for k in range(3, 1, -1):
        name = f"g_{k}_collgrad{group}"
        df[name] = (
            df["ry"].eq(-k).fillna(False) & df["collgrad"].eq(group)
        ).astype(int)
        het_vars.append(name)

    # Separate post-treatment periods 0 through 18 for each group.
    for k in range(0, 19):
        name = f"g{k}_collgrad{group}"
        df[name] = (
            df["ry"].eq(k).fillna(False) & df["collgrad"].eq(group)
        ).astype(int)
        het_vars.append(name)

fit_het = esi.eventreg(
    f"ln_wage ~ {' + '.join(het_vars)} | idcode + year",
    data=df,
    cohort="first_union",
    control_cohort="never_union",
    covariates="south",
    vcov={"CRV1": "idcode"},
)

Compare group-specific effects with lincom():

fit_het.lincom(
    {"g1_collgrad1": 1.0, "g1_collgrad0": -1.0},
    name="collgrad_diff_g1",
)

ATE And Lincom

ate() and lincom() do not re-estimate the model. They use the IW coefficient vector and IW covariance matrix:

estimate = a' beta_iw
standard error = sqrt(a' V_iw a)

For simple event-time names, fit.ate() averages all nonnegative periods, for example g0 through g18 in the nlswork example:

fit.ate(name="ATE_g0_g18")

Pass periods explicitly when you want a different window or group-specific columns:

fit.ate(["g0", "g1", "g2", "g3", "g4"], name="ATE_g0_g4")
fit_het.ate(
    [f"g{k}_collgrad1" for k in range(0, 5)],
    name="ATE_collgrad1_g0_g4",
)

Current nlswork result for the full post-treatment window:

term	estimate	standard error
ATE_g0_g18	0.111872	0.023011

lincom(weights) is the general linear-combination method. The dictionary maps coefficient names to weights, so {"g3": 1.0, "g0": -1.0} means beta_g3 - beta_g0.

fit.lincom({"g3": 1.0, "g0": -1.0}, name="g3_minus_g0")

Current nlswork result:

term	estimate	standard error
g3_minus_g0	0.017982	0.027433

Two-Way Clustered Standard Error Correctness Check

The purpose of this check is to verify the standard errors, not only the point estimates.

In a general staggered-adoption design, the IW estimator combines many cohort-specific event-study coefficients. It is hard to see from the final table whether the covariance matrix has been aggregated correctly. So the test uses a simpler design where the correct answer is known exactly: a non-staggered DID panel with only one treated cohort.

• 96 units and 10 periods;
• half the units are treated in period 5 and half are never treated;
• the true dynamic effects are fixed in the data-generating process;
• the fitted model is y ~ g_3 + g_2 + g0 + g1 + g2 + g3 | id + t;
• the covariance is clustered by both id and t.

With one treated cohort, all IW weights are one and the share covariance is zero. The Sun-Abraham IW estimator therefore collapses to the ordinary TWFE event-study regression with the same formula. That makes direct pyfixest.feols(..., vcov={"CRV1": "id + t"}) the reference result.

The required conclusion is:

• eventreg() coefficients must equal direct pyfixest/TWFE coefficients;
• eventreg().vcov() must equal the direct pyfixest/TWFE covariance matrix;
• eventreg().ate() and eventreg().lincom() standard errors must equal sqrt(a' V a) from that same covariance matrix.

Current Python check, at 1e-12 tolerance:

term	direct pyfixest estimate	eventreg estimate	direct standard error	eventreg standard error
g_3	0.480228	0.480228	0.119375	0.119375
g_2	-0.023928	-0.023928	0.119225	0.119225
g0	-0.711354	-0.711354	0.122300	0.122300
g1	1.595773	1.595773	0.122213	0.122213
g2	2.817989	2.817989	0.120473	0.120473
g3	1.957438	1.957438	0.119352	0.119352
ATE_g0_g3	1.414962	1.414962	0.120690	0.120690

Conclusion: this package passes the two-way clustered standard-error check. For the non-staggered single-cohort design, eventreg() is numerically equal to direct pyfixest/TWFE for coefficients, the full covariance matrix, and the ATE standard error.

The external parity script also compares Stata eventstudyinteract's e(V_iw) in this single-cohort family of checks. It reports Stata as an expected difference:

Stata eventstudyinteract vcov comparison: EXPECTED_DIFF

EXPECTED_DIFF means the Stata eventstudyinteract coefficient vector matches, but its reported e(V_iw) covariance does not match the direct TWFE/saturated full-FE covariance target in this diagnostic. It is not a failure of this package. It is the reason the correctness target for standard errors is direct pyfixest/TWFE and the fixed Julia implementation, not Stata's e(V_iw) in this edge case.

Exporting Tables

eventreg() returns an IW result object, not a pyfixest model object. Use the result's table methods:

fit.etable(type="df")
fit.etable(type="tex")
fit.etable(type="html", file_name="eventstudy.html")

fit.tidy().to_csv("eventstudy.csv")
fit.ate(name="ATE_g0_g18").to_csv("ate.csv")
fit.lincom({"g3": 1.0, "g0": -1.0}).to_csv("lincom.csv")

For Word output, export CSV/HTML or pass the returned pandas DataFrame to the table package you already use.

pyfixest Options

The expensive regression work is handled by pyfixest.feols(). Important pyfixest options are forwarded:

fit_fast = esi.eventreg(
    f"ln_wage ~ {' + '.join(event_vars)} | idcode + year + ind_code^year",
    data=df,
    cohort="first_union",
    control_cohort="never_union",
    covariates="south + tenure + age",
    vcov={"CRV1": "idcode + year"},
    weights="person_weight",
    weights_type="aweights",
    fixef_rm="singleton",
    fixef_tol=1e-8,
    fixef_maxiter=100_000,
    collin_tol=1e-6,
    solver="scipy.linalg.solve",
    demeaner_backend="numba",   # also rust, rust-cg, jax, cupy, scipy
    share_vcov="regression",    # "none" matches pyfixest saturated SE
    use_compression=False,
)

By default, first-stage covariance scaling follows pyfixest. For the explicit saturated/full absorbed-FE small-sample correction used in the Julia parity check, pass fixef_vcov_correction="full".

Verification

To reproduce the usage examples in this README, run:

PYTHONPATH=src python examples/nlswork_julia_replication.py

That script loads dataset/nlswork.dta, creates the first_union, ry, relative-time, binning, and heterogeneity columns, then estimates the same examples shown above.

To reproduce the cross-package comparison tables and the standard-error correctness checks, run:

PYTHONPATH=src python scripts/check_external_parity.py

The parity script checks:

• nlswork same-sample parity against pyfixest saturated, Stata, and Julia;
• multi-FE plus continuous-control parity against Stata and Julia;
• single-cohort equality to direct pyfixest TWFE, including two-way clustering;
• ate() and lincom() standard errors from the IW covariance matrix.

The Stata, R, and Julia parts run only when those tools are installed locally; otherwise the script prints SKIP for the missing external package. The Python and pyfixest checks always run.