---
title: "Credit markets, property rights, and the commons"
subtitle: "TWFE, Sun–Abraham estimates, and a joint test of cohort deviations"
date: 2026-10-08
format:
  html:
    theme: cosmo
    toc: true
    toc-depth: 2
    embed-resources: true
    self-contained-math: true
    code-fold: true
    code-summary: "Show Python"
    code-copy: true
    html-math-method: katex
    include-in-header:
      text: |
        <style>
        body {font-family: system-ui, sans-serif; font-size: 17px; line-height: 1.65;}
        h1,h2,h3 {letter-spacing: -.025em;}
        .subtitle {max-width: 55rem; color: #526477;}
        main, .cell, .cell-output-display {min-width: 0; max-width: 100%;}
        .cell-output-display {overflow-x: auto;}
        figure img, .cell-output-display img {max-width: 100%; height: auto !important;}
        table {display: block; max-width: 100%; overflow-x: auto;}
        p code, li code {overflow-wrap: anywhere; white-space: normal;}
        .cell-output-display table {font-size: .87rem;}
        .callout {border-radius: 8px;}
        .quarto-title-block {border-bottom: 3px solid #247a83; padding-bottom: 1rem;}
        pre {font-size: .8rem;}
        .katex-display {overflow-x:auto;overflow-y:hidden;max-width:100%;}
        </style>
jupyter: python3
execute:
  warning: false
  message: false
  cache: false
  error: false
  fig-dpi: 160
  fig-width: 10
  fig-height: 4.5
---

[Noack and Costello (2024)](https://doi.org/10.1086/729065) study credit markets and property rights in fisheries. This example uses their prepared stock–country panel to estimate the effect of property-rights adoption on fishing mortality. It compares static TWFE, dynamic TWFE, and Sun–Abraham estimates, then tests whether post-treatment effects differ across adoption cohorts.

## 1. Source data and variable map

```{python}
#| label: setup
#| include: false
from pathlib import Path
import hashlib
import json
import re
import warnings
import os
import subprocess
import tempfile
from importlib.metadata import version

import duckreg
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import pyfixest as pf
from IPython.display import Markdown, display
from scipy import stats

PROJECT = next(p for p in [Path.cwd(), *Path.cwd().parents]
               if (p / "repl-archives").is_dir())
HERE = PROJECT / "analyses/credit-markets"
RESULTS = HERE / "results"
RESULTS.mkdir(exist_ok=True)
SOURCE = (PROJECT / "repl-archives/twfe-microdata-replication-private/"
          "outputs/estimation_2026_10_06/group_c/CMPRC-balanced-original.csv")
SOURCE_HASH = hashlib.sha256(SOURCE.read_bytes()).hexdigest()
VERSIONS = {p: version(p) for p in ["polars", "pyfixest", "duckreg", "numpy", "pandas"]}
plt.rcParams.update({"axes.spines.top": False, "axes.spines.right": False,
                     "font.size": 11, "axes.titleweight": "bold"})

def table(frame, digits=5):
    display(Markdown(frame.to_markdown(index=False, floatfmt=f".{digits}g")))
```

The source is the prepared panel `CMPRC-balanced-original.csv`. The archive builder, `scripts/estimation/papers_group_c.py::commons()`, defines its variables.

```{python}
#| label: source-head
df = pl.read_csv(SOURCE)
df.head()
```

| Role | Archived variable | Meaning and handling |
|---|---|---|
| Outcome, **Y** | `y` | One outcome: `asinh(ffmsy)`, already computed by the archive builder as `log(ffmsy + sqrt(1 + ffmsy**2))`. `ffmsy` measures fishing mortality relative to its maximum-sustainable-yield benchmark. No further rescaling is applied. |
| Treatment, **W** | `d` | Binary, cardinality two: 0 and 1. The builder assigns EEZ declaration time to units with `max(pr) > 0.5`, then sets `d = 1` from that time onward. The underlying continuous `pr` measure is not in this prepared file. |
| Unit, **i** | `unit` | Integer code for a stock–country pair; `assignment` duplicates this identifier. |
| Period, **t** | `time` | Calendar year. There are no gaps in the annual time grid. |
| Covariates, **X_it** | None | The prepared specification is unconditional apart from fixed effects.  |
| Inference | `cluster` | Country identifier. Country effects are absorbed by stock–country unit effects, but uncertainty is clustered at country level. |
| Weights | `weight` | Always 1: estimation is unweighted. |
| Timing/provenance | `cohort`, `obs_id` | Supplied adoption year (`cohort = 0` for never-treated) and source-row identifier. Neither is an additional treatment variant. |

```{python}
#| label: data-checks
assert df.select(pl.struct(["unit", "time"]).n_unique()).item() == df.height
assert df.select(pl.all().null_count()).row(0) == (0,) * df.width
assert set(df["d"].unique()) == {0, 1}
assert df["weight"].unique().to_list() == [1.0]
assert (df["assignment"] == df["unit"]).all()
assert df.group_by("unit").len()["len"].unique().to_list() == [53]
assert df.group_by("unit").agg(pl.col("cluster").n_unique())["cluster"].max() == 1
assert df.sort(["unit", "time"]).select(
    pl.col("d").diff().over("unit").min()).item() == 0
assert df["time"].unique().sort().to_list() == list(range(1960, 2013))

table(pd.DataFrame({"Quantity": ["Observations", "Stock–country units", "Countries",
                                 "Years", "Outcome variables", "Treatment values"],
                    "Value": [df.height, df["unit"].n_unique(), df["cluster"].n_unique(),
                              "1960–2012 (53)", "1: y", "2: {0, 1}"]}))
```

### Construct adoption time and relative time

For each unit, adoption time `g` is the minimum observed year with `d = 1`. Never-treated units receive PyFixest's sentinel `g = 0`; their event time remains missing. On the annual grid, relative time is simply `time - g`.

```{python}
#| label: treatment-timing
df = df.with_columns(
    pl.when(pl.col("d") == 1).then(pl.col("time")).otherwise(None)
      .min().over("unit").fill_null(0).cast(pl.Int64).alias("g")
).with_columns(
    pl.when(pl.col("g") > 0).then(pl.col("time") - pl.col("g"))
      .otherwise(None).alias("event_time")
)
first_year = df["time"].min()
always = df.filter((pl.col("g") > 0) & (pl.col("g") <= first_year))
eligible = df.filter((pl.col("g") == 0) | (pl.col("g") > first_year))
assert (eligible["g"] == eligible["cohort"]).all()
assert eligible.select(
    (((pl.col("g") > 0) & (pl.col("time") >= pl.col("g"))).cast(pl.Int64)
     == pl.col("d")).all()).item()

data = eligible.to_pandas()
data["ever_treated"] = data["g"].gt(0).astype(int)
# The never-treated sentinel is used only with an ever-treated interaction.
data["relative_time"] = data["event_time"].fillna(-1).astype(int)
unit_info = data[["unit", "cluster", "g"]].drop_duplicates()
cohort_sizes = unit_info.loc[unit_info.g.gt(0)].groupby("g").size()
cohorts = unit_info.groupby("g").agg(units=("unit", "size"),
                                    countries=("cluster", "nunique")).reset_index()
table(cohorts)
```

```{python}
#| echo: false
display(Markdown(
    f"**Common estimation sample:** {len(data):,} observations, "
    f"{data.unit.nunique():,} stock–country units, {data.cluster.nunique()} countries; "
    f"{int((unit_info.g > 0).sum())} eventually treated units in {len(cohort_sizes)} adoption cohorts "
    f"and {int((unit_info.g == 0).sum())} never-treated units. "
    f"The {always['unit'].n_unique()} always-treated units ({always.height:,} observations) "
    "are excluded from **all** models. Their supplied adoption years are 1950/1952, before "
    "the panel begins; the within-panel minimum would incorrectly suggest 1960. "
    "There are no treatment reversals or additional missing-value exclusions."
))
```

## 2. Visualize adoption timing

```{python}
#| label: fig-adoption
#| fig-cap: "Purple = untreated; yellow = treated. A seeded sample of 100 eligible stock–country units, sorted by adoption timing. The regressions use all 723 eligible units."
rng = np.random.default_rng(20261008)
plot_units = rng.choice(np.sort(data.unit.unique()), size=min(100, data.unit.nunique()),
                        replace=False)
fig, ax = plt.subplots(figsize=(10, 5))
pf.panelview(data[data.unit.isin(plot_units)], unit="unit", time="time", treat="d",
             sort_by_timing=True, xlab="Calendar year", ylab="Sampled stock–country unit",
             title="Treatment adoption", noticks=True, ax=ax)
tick_years = [1960, 1970, 1980, 1990, 2000, 2012]
ax.set_xticks(np.array(tick_years) - first_year, labels=tick_years)
plt.tight_layout()
plt.show()
```

**Subsampling applies only to this plot:** 100 units without replacement, seed 20261008. It does not change the estimation sample, cohort weights, or covariance calculations.

## 3. Static TWFE and dynamic specifications

### Static two-way fixed effects

This benchmark regresses `y` on `d` with stock–country and calendar-year effects. CRV1 standard errors are clustered by country. Under heterogeneous effects with staggered adoption, the coefficient need not equal a convex average ATT.

```{python}
#| label: static-twfe
static = pf.feols("y ~ d | unit + time", data=data, vcov={"CRV1": "cluster"})
table(static.tidy().reset_index())
```

### Dynamic TWFE

The dynamic specification includes relative-time indicators interacted with ever-treated status, with year −1 omitted.

```{python}
#| label: dynamic-twfe
dynamic_twfe = pf.feols(
    "y ~ i(relative_time, ever_treated, ref=-1) | unit + time",
    data=data, vcov={"CRV1": "cluster"},
)
twfe_period = dynamic_twfe.tidy().copy()
twfe_period.index = twfe_period.index.str.extract(r"relative_time::([^:]+)",
                                                 expand=False).astype(float)
twfe_period.index.name = "period"
twfe_period = twfe_period.sort_index()
assert len(twfe_period) > 1
```

All observed event times, −50 to +45, enter the regression. Never-treated units have zero exposure to every event-time indicator.

### Saturated / Sun–Abraham event study

```{python}
#| label: saturated
with warnings.catch_warnings(record=True) as sa_warnings:
    warnings.simplefilter("always")
    saturated = pf.event_study(
        data=data, yname="y", idname="unit", tname="time", gname="g",
        cluster="cluster", estimator="saturated", att=False,
    )
sa_period = saturated.aggregate(agg="period", weighting="shares").astype(float)
assert np.isfinite(sa_period.to_numpy()).all()
assert saturated._N == dynamic_twfe._N == static._N == len(data)
```

The saturated model estimates cohort-by-relative-time effects with year −1 omitted for each cohort. `.aggregate(agg="period", weighting="shares")` returns event-time effects weighted by each contributing cohort's observation share. Never-treated observations anchor the comparison.

```{python}
#| label: fig-event-study
#| fig-cap: "Dynamic TWFE and share-weighted Sun–Abraham estimates. Error bars are pointwise 95% normal intervals using country-clustered covariance; not simultaneous bands. The displayed window is −20 to +20, but both models were fitted on all observed event times."
fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
for ax, frame, label, color in [
    (axes[0], twfe_period, "Dynamic TWFE", "#916440"),
    (axes[1], sa_period, "Sun–Abraham: cohort shares", "#147d87"),
]:
    shown = frame.loc[(frame.index >= -20) & (frame.index <= 20)]
    ax.errorbar(shown.index, shown["Estimate"], yerr=1.96*shown["Std. Error"],
                fmt="o", ms=3, capsize=2, color=color, elinewidth=1)
    ax.scatter([-1], [0], marker="s", facecolors="white", edgecolors=color, zorder=4)
    ax.axhline(0, color="#536477", lw=.8)
    ax.axvline(-.5, color="#88939c", lw=.8, ls="--")
    ax.set(title=label, xlabel="Years relative to adoption")
axes[0].set_ylabel("Effect on asinh(ffmsy)")
plt.tight_layout()
plt.show()
```

### Aggregate post-treatment effects

The scalar ATT averages post-treatment effects using treated-unit-year shares. Its standard error uses the full coefficient covariance matrix. Dynamic TWFE uses the same post-period weights.

```{python}
#| label: scalar-aggregation
sa_names = saturated.coef().index
parts = sa_names.str.extract(r"rel_time::(?P<e>[^:]+):first_treated_period::(?P<g>[^:]+)")
parts = parts.astype({"e": float, "g": int})
assert parts.notna().all().all()
periods = sa_period.index.to_numpy(dtype=float)
sa_weights = np.zeros((len(periods), len(sa_names)))
for row, event in enumerate(periods):
    hit = parts["e"].eq(event).to_numpy()
    sa_weights[row, hit] = parts.loc[hit, "g"].map(cohort_sizes)
    sa_weights[row] /= sa_weights[row].sum()
np.testing.assert_allclose(sa_weights @ saturated._beta_hat,
                           sa_period["Estimate"], atol=1e-9)
np.testing.assert_allclose(np.sqrt(np.diag(sa_weights @ saturated._vcov @ sa_weights.T)),
                           sa_period["Std. Error"], atol=1e-9)

post_mass = data.loc[data.d.eq(1)].groupby("event_time").size()
post_weights = post_mass.reindex(periods, fill_value=0).to_numpy(dtype=float)
post_weights /= post_weights.sum()
r_sa = post_weights @ sa_weights
twfe_events = dynamic_twfe.coef().index.str.extract(r"relative_time::([^:]+)",
                                                  expand=False).astype(float)
r_twfe = post_mass.reindex(twfe_events, fill_value=0).to_numpy(dtype=float)
r_twfe /= r_twfe.sum()

def linear_summary(label, fit, weights):
    estimate = float(weights @ fit._beta_hat)
    se = float(np.sqrt(weights @ fit._vcov @ weights))
    return {"Summary": label, "Estimate": estimate, "Country-cluster SE": se,
            "95% lower": estimate-1.96*se, "95% upper": estimate+1.96*se}

pooled = pd.DataFrame([
    linear_summary("Static TWFE coefficient", static, np.ones(1)),
    linear_summary("Dynamic TWFE: post unit-year average", dynamic_twfe, r_twfe),
    linear_summary("Sun–Abraham: post unit-year ATT", saturated, r_sa),
])
table(pooled)
```

These scalar intervals use a normal critical value for comparability with the event-time aggregation; the earlier native static output uses PyFixest's finite-cluster t reference. The Sun–Abraham scalar weights cohort/event-time cells in proportion to treated observations. Earlier cohorts contribute more post-treatment years; this is not an equal-cohort or common-horizon average.

## 4. Joint F-test of cohort deviations

The [TestingInEventStudies implementation](https://github.com/apoorvalal/TestingInEventStudies/blob/67fcce76736e78863a53da9696b405ab31037672/code/saturated.py) constructs common event-time indicators and interactions for all treated cohorts except the first. The interactions are deviations from the reference cohort's effect path.

$$
Y_{it}=\alpha_i+\lambda_t+\sum_{e\ne-1}\beta_e D_{it}^{e}
+\sum_{g\ne g_0}\sum_{e\ne-1}\delta_{ge}D_{it}^{e}1\{G_i=g\}+\varepsilon_{it}.
$$

The joint null is $\delta_{ge}=0$ for every estimated deviation with $e>0$. Adoption-year and pre-treatment deviations remain unrestricted.

```{python}
#| label: reference-deviations-f-test
# Run the unchanged design builder in its pinned reference environment.
with tempfile.TemporaryDirectory() as temporary:
    input_csv = Path(temporary) / "sample.csv"
    data.to_csv(input_csv, index=False)
    spec = {
        "reference_repo": str(PROJECT / "TestingInEventStudies"),
        "data": str(input_csv), "outcome": "y", "treatment": "d",
        "unit": "unit", "time": "time", "cluster": "cluster",
        "test_mode": "classical", "test_horizon": None,
        "output": str(RESULTS / "deviations-f-test.json"),
    }
    completed = subprocess.run(
        [str(PROJECT / "TestingInEventStudies/.venv/bin/python"),
         str(HERE / "credit-markets-reference-f-test.py"), json.dumps(spec)],
        env=dict(os.environ, OPENBLAS_NUM_THREADS="1", OMP_NUM_THREADS="1",
                 NUMBA_NUM_THREADS="2", MPLBACKEND="Agg", PYTHONDONTWRITEBYTECODE="1"),
        capture_output=True, text=True,
    )
    if completed.returncode:
        raise RuntimeError(completed.stderr)
f_test = json.loads((RESULTS / "deviations-f-test.json").read_text())
table(pd.DataFrame([{
    "F": f_test["F"], "Numerator df": f_test["df_num"],
    "Denominator df": f_test["df_den"], "p-value": f_test["pvalue"],
}]))
```

```{python}
#| echo: false
display(Markdown(
    f"The imported design has **{f_test['design_columns']} columns** after absorbing "
    f"unit and year effects. The restricted regression removes **{f_test['df_num']} "
    f"post-treatment deviations**. **F({f_test['df_num']}, {f_test['df_den']}) = "
    f"{f_test['F']:.3f}, p = {f_test['pvalue']:.3g}.**"
))
```

This is the homoskedastic nested-model F-test, computed from the two residual sums of squares. The coefficient tables above retain country-clustered standard errors.

$$
F=\frac{(RSS_R-RSS_U)/q}{RSS_U/(n-\operatorname{rank}(X_U))}.
$$

[Download the F-test runner](credit-markets-reference-f-test.py). It imports `test_treatment_heterogeneity(..., retmod=True)` and uses the resulting deviations design; the source repository is unchanged.

## 5. Results and reproducibility record

```{python}
#| echo: false
display(Markdown(
    f"The static TWFE coefficient is **{pooled.iloc[0]['Estimate']:.3f}** "
    f"(cluster SE {pooled.iloc[0]['Country-cluster SE']:.3f}); the saturated model's "
    f"treated-unit-year-weighted post ATT is **{pooled.iloc[2]['Estimate']:.3f}** "
    f"(cluster SE {pooled.iloc[2]['Country-cluster SE']:.3f}). "
    "The models use the same observations, but do not impose the same effect structure "
    "or assign the same implicit treatment-effect weights."
))
```

[Download the executable Quarto source](credit-markets-event-study.qmd). The source expects the existing `big-did-replication` directory layout and access to its prepared CSV. Environment manifests: [pyproject.toml](credit-markets-pyproject.toml), [uv.lock](credit-markets-uv.lock).

From the project root:

```bash
uv sync --locked
QUARTO_PYTHON="$PWD/.venv/bin/python" \
  OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 \
  quarto render analyses/credit-markets/credit-markets-event-study.qmd
```

```{python}
#| label: save-results
#| include: false
static.tidy().to_csv(RESULTS / "static-twfe.csv")
twfe_period.to_csv(RESULTS / "dynamic-twfe.csv")
sa_period.to_csv(RESULTS / "sun-abraham-period.csv")
pooled.to_csv(RESULTS / "pooled-estimates.csv", index=False)
cohorts.to_csv(RESULTS / "cohorts.csv", index=False)
record = {
    "versions": VERSIONS, "source_sha256": SOURCE_HASH,
    "raw_observations": df.height, "eligible_observations": len(data),
    "eligible_units": data.unit.nunique(), "country_clusters": data.cluster.nunique(),
    "always_treated_units_excluded": always["unit"].n_unique(),
    "never_treated_units": int((unit_info.g == 0).sum()),
    "adoption_cohorts": len(cohort_sizes), "plot_units": len(plot_units),
    "plot_seed": 20261008,
    "sa_coefficients": len(sa_names), "dynamic_twfe_coefficients": len(twfe_period),
    "deviations_f_test": f_test,
    "pooled": pooled.to_dict("records"),
    "sa_warnings": [str(w.message) for w in sa_warnings],
}
(RESULTS / "analysis.json").write_text(json.dumps(record, indent=2) + "\n")
```

```{python}
#| echo: false
table(pd.DataFrame(list(VERSIONS.items()), columns=["Package", "Version"]))
display(Markdown(f"Source SHA-256: `{SOURCE_HASH}`. Quarto 1.8.27; Python 3.12.9."))
```

### References and implementation sources

- Noack and Costello (2024), [published article](https://doi.org/10.1086/729065); [working paper](https://www.nber.org/papers/w29889).
- PyFixest [LLM documentation index](https://pyfixest.org/llms.txt), [DiD tutorial and treatment-timing plots](https://pyfixest.org/tutorials/difference-in-differences.html#examining-treatment-timing), [event-study API](https://pyfixest.org/reference/did.estimation.event_study.html), and [saturated event-study API](https://pyfixest.org/reference/did.saturated_twfe.SaturatedEventStudy.html).
- Sun and Abraham (2021), [Estimating dynamic treatment effects in event studies with heterogeneous treatment effects](https://doi.org/10.1016/j.jeconom.2020.09.006); Lal (2025), [Testing in Event Studies](https://arxiv.org/abs/2503.05125).