Credit markets, property rights, and the commons

TWFE, Sun–Abraham estimates, and a joint test of cohort deviations

Published

October 8, 2026

Noack and Costello (2024) 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

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

Show Python 
df = pl.read_csv(SOURCE)
df.head()
shape: (5, 9)
obs_id unit time assignment cluster y d weight cohort
i64 i64 i64 i64 i64 f64 i64 f64 i64
1327 1686 1960 1686 5 1.137982 0 1.0 1967
1328 1686 1961 1686 5 1.137982 0 1.0 1967
1329 1686 1962 1686 5 1.719052 0 1.0 1967
1330 1686 1963 1686 5 1.078451 0 1.0 1967
1331 1686 1964 1686 5 3.163037 0 1.0 1967
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.
Show Python 
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}"]}))
Quantity Value
Observations 39432
Stock–country units 744
Countries 47
Years 1960–2012 (53)
Outcome variables 1: y
Treatment values 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.

Show Python 
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)
g units countries
0 194 40
1967 11 2
1969 2 1
1971 1 1
1976 10 3
1977 255 10
1978 81 4
1979 46 3
1980 7 2
1982 4 1
1983 82 1
1984 1 1
1995 15 1
2000 1 1
2005 12 1
2010 1 1

Common estimation sample: 38,319 observations, 723 stock–country units, 47 countries; 529 eventually treated units in 15 adoption cohorts and 194 never-treated units. The 21 always-treated units (1,113 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

Show Python 
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()
Figure 1: Purple = untreated; yellow = treated. A seeded sample of 100 eligible stock–country units, sorted by adoption timing. The regressions use all 723 eligible units.

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.

Show Python 
static = pf.feols("y ~ d | unit + time", data=data, vcov={"CRV1": "cluster"})
table(static.tidy().reset_index())
Coefficient Estimate Std. Error t value Pr(> t )
d -0.60473 0.092302 -6.5516 4.2968e-08 -0.79052 -0.41893

Dynamic TWFE

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

Show Python 
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

Show Python 
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.

Show Python 
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()
Figure 2: 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.

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.

Show Python 
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)
Summary Estimate Country-cluster SE 95% lower 95% upper
Static TWFE coefficient -0.60473 0.092302 -0.78564 -0.42382
Dynamic TWFE: post unit-year average -0.67764 0.11325 -0.89961 -0.45568
Sun–Abraham: post unit-year ATT -0.68295 0.085152 -0.84985 -0.51605

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

Show Python 
# 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"],
}]))
F Numerator df Denominator df p-value
2.9413 379 36764 5.8121e-72

The imported design has 780 columns after absorbing unit and year effects. The restricted regression removes 379 post-treatment deviations. F(379, 36764) = 2.941, p = 5.81e-72.

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. It imports test_treatment_heterogeneity(..., retmod=True) and uses the resulting deviations design; the source repository is unchanged.

5. Results and reproducibility record

The static TWFE coefficient is -0.605 (cluster SE 0.092); the saturated model’s treated-unit-year-weighted post ATT is -0.683 (cluster SE 0.085). 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. The source expects the existing big-did-replication directory layout and access to its prepared CSV. Environment manifests: pyproject.toml, uv.lock.

From the project root:

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
Package Version
polars 2.0.0
pyfixest 0.60.0
duckreg 0.4.1
numpy 2.5.3
pandas 3.0.6

Source SHA-256: 9b62b770cc25e997c5c5c2f01f8b8a3c7307c240dc98ee455bae57ac2bc6529b. Quarto 1.8.27; Python 3.12.9.

References and implementation sources