Metadata-Version: 2.4
Name: eventstudyinteracts
Version: 0.1.0
Summary: A pyfixest-backed Python implementation of interaction-weighted event-study estimators.
Author: EventStudyInteracts.py maintainers
License-Expression: MIT
Project-URL: Homepage, https://github.com/xiaobaaaa/EventStudyInteracts.py
Project-URL: Source, https://github.com/xiaobaaaa/EventStudyInteracts.py
Project-URL: Issues, https://github.com/xiaobaaaa/EventStudyInteracts.py/issues
Keywords: event study,difference-in-differences,Sun Abraham,pyfixest,fixed effects
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.25.2
Requires-Dist: pandas>=1.1.0
Requires-Dist: pyfixest>=0.50.1
Requires-Dist: scipy<1.16,>=1.6
Provides-Extra: dev
Requires-Dist: pytest>=8; extra == "dev"
Requires-Dist: ruff>=0.8; extra == "dev"
Dynamic: license-file

# 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:

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

To install it for Claude Code as a user skill:

```bash
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:

```python
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:

```bash
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:

```bash
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`.

```python
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:

```python
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:

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

## Plotting

For a pyfixest-style event-study plot:

```python
fit.iplot_aggregate()
```

To save the plot:

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

Manual plotting is also straightforward:

```python
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:

```python
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.

```python
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:

```python
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.

```python
# 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.

```python
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.

```python
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()`:

```python
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:

```text
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:

```python
fit.ate(name="ATE_g0_g18")
```

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

```python
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`.

```python
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:

```text
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:

```python
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:

```python
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:

```bash
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:

```bash
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.
