---
title: "Credit Markets: alternative counterfactual models"
subtitle: "CS-DID, SDID, matrix completion, horizontal ridge and partially pooled staggered synthetic control"
date: 2026-10-09
format:
  html:
    theme: cosmo
    toc: true
    embed-resources: true
    self-contained-math: true
    html-math-method: katex
    code-fold: true
    code-summary: "Code"
    code-copy: true
    include-in-header:
      text: |
        <style>
        body{font-family:system-ui,sans-serif;line-height:1.6}
        main,.cell,.cell-output-display{min-width:0;max-width:100%}
        .cell-output-display,table{overflow-x:auto}
        table{display:block;font-size:.86rem}
        img{max-width:100%;height:auto!important}
        .katex-display{overflow-x:auto;overflow-y:hidden}
        p code,li code{white-space:normal;overflow-wrap:anywhere}
        </style>
jupyter: python3
execute:
  warning: false
  error: false
  cache: false
  fig-dpi: 160
---

[Noack and Costello (2024)](https://doi.org/10.1086/729065) study property rights and fishing mortality. This comparison extends the [Credit Markets replication](study-20.html) with Callaway–Sant’Anna DID and four synthetic-control or outcome-model estimators, on the same sample. CS-DID reports both the dynamic/event-study aggregate and the simple treated-unit-year aggregate. The latter matches the averaging weights used by the other counterfactual summaries.

[Results JSON](credit-markets-counterfactuals-results.json) · [Quarto source](credit-markets-counterfactuals.qmd) · [Code and environment](credit-markets-counterfactuals-code.zip) · [Archive index](index.html)

```{python}
#| label: setup
import sys
from pathlib import Path

import crabbymetrics as cm
import matplotlib.pyplot as plt
```

```{python}
#| label: project-paths
ROOT = next(
    p
    for p in [Path.cwd(), *Path.cwd().parents]
    if (p / "templates/full_archive").exists()
)
OUT = ROOT / "analyses/credit-markets/crabbymetrics"
sys.path.insert(0, str(OUT / "templates"))
```

```{python}
#| label: prepare-panel
from fit_crabby import fit_one
from report_support import (
    adoption_plot,
    event_plot,
    prepare_panel,
    publish_results,
    results,
    scalar_plot,
    show_table,
)

data, Y, W, baseline = prepare_panel()
```

## Data and sample

The outcome is `asinh(ffmsy)`: fishing mortality relative to the maximum-sustainable-yield benchmark. Treatment is EEZ/property-rights adoption, coded 0 or 1. Units are stock–country pairs; time is annual, 1960–2012. There are no additional covariates in this specification.

All models use the existing **723-unit × 53-year panel**, comprising 529 treated units in 15 adoption cohorts and 194 never-treated units. The 21 units already treated in 1960 are excluded. No further observations are dropped. The panel contains **38,319 observations and 17,811 treated unit-years**, across 47 countries.

The outcome is divided by **1.3291711871762193**, the original full-source sample SD before exclusions. It is not demeaned. Adoption is the first treated year within a unit; event time is year minus adoption year.

```{python}
#| label: data-head
preview = data[["time", "d", "g", "y"]].copy()
preview["time"] += 1959
preview["g"] = preview["g"].where(preview["g"].eq(0), preview["g"] + 1959)
preview = preview.rename(
    columns={
        "time": "Year",
        "d": "Treatment",
        "g": "Adoption year",
        "y": "Outcome / SD",
    }
)
preview.head()
```

```{python}
#| label: adoption
adoption_plot(W)
plt.show()
```

The adoption plot shows a seeded sample of 100 units. Estimation uses all 723.

## Estimators

| Estimator | Counterfactual model and settings |
|:--|:--|
| CS-DID | Python `csdid.ATTgt`, doubly robust estimation with no covariates, not-yet-treated controls and universal cohort-specific year −1 reference. `aggte(typec="dynamic")` supplies the event path and equal-event-time scalar; `aggte(typec="simple")` supplies the matched treated-unit-year scalar. |
| SDID | `crabbymetrics.SyntheticDID`: donor and pre-period weights fitted separately for each adoption cohort; default data-dependent penalties. |
| Matrix completion | `crabbymetrics.MatrixCompletion`: unit effects, time effects and a nuclear-norm-regularized low-rank component, fitted to all untreated cells. `lambda_fraction=0.25`; exact SVD; objective tolerance $10^{-12}$. |
| Horizontal ridge | `crabbymetrics.HorizontalPanelRidge`: regress each cohort’s pre-treatment mean on the never-treated outcome histories, with an intercept and ridge penalty 1. Predict using contemporaneous donor outcomes. |
| Multisynth (individual) | `time_cohort=FALSE`: 529 unit-level synthetic controls fitted jointly, with partial pooling. `nu=NULL`, `lambda=0`, fixed-effect adjustment and nonnegative donor weights. |
| Multisynth (cohort-averaged) | `time_cohort=TRUE`: 15 cohort-level synthetic controls fitted jointly, with partial pooling and the same other settings. |

CS-DID uses never-treated units plus other cohorts untreated in both comparison periods. SDID, horizontal ridge and multisynth use the 194 never-treated donors. Matrix completion also uses untreated observations from units that adopt later. The three crabbymetrics estimators follow the [staggered-panel vignette](https://apoorvalal.github.io/crabbymetrics/examples/staggered-panel-event-study.html), with a tighter matrix-completion stopping tolerance.

### Callaway–Sant’Anna DID

[CS-DID](https://github.com/d2cml-ai/csdid) estimates one ATT for each adoption cohort and calendar year. This fit uses `est_method="dr"`, `control_group="notyettreated"`, `xformla="y ~ 1"`, `anticipation=0` and `base_period="universal"`. With no covariates, each group-time point estimate is the cohort’s mean outcome change from its year −1 minus the corresponding change among never-treated and eligible not-yet-treated units. The focal cohort is excluded from its own controls, and eligibility requires remaining untreated through both the baseline and comparison year.

```{python}
#| label: csdid
from fit_csdid import fit_csdid

csdid = fit_csdid(data)
from csdid_aggregations import aggregation_table

show_table(aggregation_table(csdid))
```

The same not-yet-treated group-time estimates support three distinct scalar aggregations, as defined in the [aggregation documentation](https://bcallaway11.github.io/did/reference/aggte.html):

- **Dynamic (`typec="dynamic"`):** average cohort effects using cohort-size weights within each event time, then average the 46 post-event estimates equally over years **0–45**. `balance_e=None` leaves cohort composition free to change with the horizon.
- **Simple (`typec="simple"`):** weight every available post-treatment cohort–year ATT by that cohort's size. On this complete balanced panel, this is exactly the **17,811 treated-unit-year average** used by Sun–Abraham and the other counterfactual summaries. It is the matched-aggregation CS-DID row in the main table and figure.
- **Group (`typec="group"`):** first average post-treatment years within each cohort, then average those cohort means using cohort-size weights. This gives each treated unit equal total weight regardless of its number of observed post-treatment years. It does **not** match the treated-unit-year target here.

Let $n_g$ be cohort size, $T_g$ its observed post-period count, and $\bar\tau_g=T_g^{-1}\sum_{t\geq g}\widehat{ATT}(g,t)$. Then

$$
\widehat\theta_{simple}=
\frac{\sum_g n_g T_g\bar\tau_g}{\sum_g n_g T_g},
\qquad
\widehat\theta_{group}=
\frac{\sum_g n_g\bar\tau_g}{\sum_g n_g}.
$$

For event time $e$, let $G_e$ contain cohorts observed at that exposure. The dynamic scalar is

$$
\widehat\theta_{dynamic}=\frac{1}{46}\sum_{e=0}^{45}
\frac{\sum_{g\in G_e}n_g\widehat{ATT}(g,g+e)}{\sum_{g\in G_e}n_g}.
$$

The group and simple weights differ because $T_g$ ranges from **3 to 46 years**. The two main CS-DID rows therefore are not separate model fits: they use identical controls, observations and group-time effects, with different aggregation weights. The group summary is retained in this section and JSON as an additional diagnostic, not substituted for the matched-cell row.

All three standard errors use analytical influence functions, summed within the 47 country clusters (`clustervar="cluster"`, `bstrap=False`, `cband=False`). They include uncertainty in estimated cohort shares. Independent calculations reproduce each scalar, its influence function and its standard error. Intervals are pointwise 95% normal intervals, not simultaneous bands.

### SDID

This implementation handles staggered adoption by **fitting a separate block-adoption SDID for each of the 15 cohorts**. Each fit contains that cohort’s treated units and the 194 never-treated donors, split into pre/post periods at the cohort’s adoption date. Other adoption cohorts are excluded from its donor pool, including units not yet treated at that date.

Donor weights match the cohort’s mean pre-treatment outcome history. Pre-period weights match the donors’ average post-treatment outcomes using their pre-treatment observations. Both sets of weights are nonnegative and sum to one, with intercepts and the regularization settings below. The cohort ATT is the mean post-treatment treated–synthetic gap minus the time-weighted pre-treatment gap. Cohort ATTs are then weighted by treated-unit count times observed post-period count, yielding the same 17,811-cell target. **There is no cross-cohort partial-pooling parameter.** See the [pinned implementation](https://github.com/apoorvalal/crabbymetrics/blob/896ceaf9fa885524409ad7ddd15874034ef27033/src/estimators/synthetic.rs#L557).

```{python}
sdid = fit_one("SDID", cm.SyntheticDID(), {}, Y, W)
```

### Matrix completion

```{python}
mc_settings = {
    "lambda_fraction": 0.25,
    "max_iterations": 1000,
    "tolerance": 1e-12,
    "fit_unit_effects": True,
    "fit_time_effects": True,
    "svd_method": "exact",
}
mc = fit_one("Matrix completion", cm.MatrixCompletion(**mc_settings), mc_settings, Y, W)
```

### Horizontal ridge

```{python}
ridge = fit_one(
    "Horizontal ridge", cm.HorizontalPanelRidge(1.0), {"penalty": 1.0}, Y, W
)
```

### Partially pooled staggered synthetic control

Two R **augsynth** specifications are reported on the same sample:

- **Individual-unit fits (`time_cohort=FALSE`):** one donor-weight vector for each of the 529 treated units, fitted jointly to balance individual and pooled pre-treatment outcomes.
- **Cohort-averaged fits (`time_cohort=TRUE`):** one donor-weight vector for each of the 15 adoption cohorts, fitted jointly to balance cohort-mean and pooled pre-treatment outcomes. This uses cohort-mean histories, like SDID and horizontal ridge, while retaining cross-cohort partial pooling.

These are different fitting objectives, not merely two ways to average a single set of fitted effects. Both reported ATTs use the **same 17,811 treated-unit-year weights**, with `fixedeff=TRUE`, `scm=TRUE` and `lambda=0`. Each fit selects its own shared pooling parameter using `nu=NULL`.

The horizon is set to `n_leads=46`, covering event years 0–45 wherever observed, with `n_lags=50` using all available pre-treatment history. This preserves all **17,811 treated unit-years**. The default three-year post horizon would change the target. At the specified full horizon, the package's donor-eligibility rule leaves exactly the same **194 never-treated donors** for every adoption cohort; the fitted donor masks are checked directly.

```{python}
#| label: multisynth
from report_support import multisynth_comparison, run_multisynth

multisynth = run_multisynth()
multisynth_individual = run_multisynth(individual=True)
show_table(multisynth_comparison([multisynth_individual, multisynth]), digits=5)
```

The cohort-averaged R call is shown below. The individual-unit call changes `time_cohort` to `FALSE`; the source bundle includes both executable scripts.

```r
fit <- multisynth(
  y ~ d, stock_id, time, panel,
  n_leads = 46, n_lags = 50,
  nu = NULL, lambda = 0,
  fixedeff = TRUE, scm = TRUE, time_cohort = TRUE,
  eps_abs = 1e-6, eps_rel = 1e-6
)
effect <- predict(fit, att = TRUE)
```

## Pooling and tuning parameters

**Each multisynth specification uses one shared partial-pooling parameter $\nu$:** 0.15163794 for the 529 individual-unit fits and 0.60511294 for the 15 cohort-averaged fits. Neither selects $\nu$ separately by cohort. In the implemented objective, $\nu=0$ weights only separate fit imbalance and $\nu=1$ weights only pooled imbalance, with package-specific normalization of both terms.

With `nu=NULL`, the package first fits separate synthetic controls and sets $\nu$ to its ratio of pooled to average separate-fit imbalance. This is a data-derived default, not cross-validation. Both realized values are reported above and repeated in the cohort table below. The weight-regularization parameter is $\lambda=0$; this is distinct from $\nu$, and no ridge augmentation is requested.

SDID and horizontal ridge still fit cohort-mean histories independently, with no cross-cohort partial-pooling parameter. Matrix completion estimates a shared low-rank structure jointly, but does not have a multisynth-style $\nu$.

| Estimator | Pooling used in these fits | Pooling parameter | Regularization |
|:--|:--|:--|:--|
| CS-DID | Group-time DID contrasts | N/A | None in this intercept-only specification |
| SDID | Separate cohort-mean fits; no cross-cohort coupling | N/A | Cohort-specific $\zeta_\omega$ and $\zeta_\lambda$, shown below |
| Horizontal ridge | Separate cohort-mean fits; no cross-cohort coupling | N/A | Ridge penalty 1 for every cohort |
| Multisynth (individual) | Jointly fitted unit-level weights; individual and pooled balance | Shared $\nu=0.15163794$ | Weight penalty $\lambda=0$ |
| Multisynth (cohort-averaged) | Jointly fitted cohort-level weights; cohort and pooled balance | Shared $\nu=0.60511294$ | Weight penalty $\lambda=0$ |
| Matrix completion | One shared low-rank matrix fitted jointly to all untreated cells | No separate $\nu$ | Nuclear-norm penalty `lambda_l=0.0010985526453507334`, from `lambda_fraction=0.25` |

### Cohort-specific settings

All synthetic-control fits use 194 never-treated donors. The table reports actual tuning values in outcome-SD units. SDID uses its data-dependent defaults: $\zeta_{\omega,g}=(N_g T_{\mathrm{post},g})^{1/4}\widehat\sigma_g$ and $\zeta_{\lambda,g}=10^{-6}\widehat\sigma_g$, where $\widehat\sigma_g$ is computed from donor pre-treatment first differences. Both multisynth specifications’ shared $\nu$ values and common $\lambda$ are repeated by cohort for reference; horizontal ridge's penalty is 1 throughout.

```{python}
#| label: cohort-tuning
from report_support import show_cohort_tuning

tuning_models, _, _ = results()
show_cohort_tuning(tuning_models)
```

The results JSON includes explicit `partial_pooling` metadata for every estimator and a `cohort_tuning` record for each adoption cohort. A null parameter for the other methods means **not applicable**, not an unreported fitted value.

## Controls, estimands and the CS-DID difference

The calendar panel is balanced: every one of the 723 units is observed in all 53 years. **No additional covariates** enter this specification (`xformla="y ~ 1"` for CS-DID). Unit/time fixed effects, lagged-outcome histories, learned donor weights and low-rank factors are parts of the respective estimators, not an added common covariate vector. The outcome scaling and eligible observations are identical across rows.

| Method | Comparison information |
|:--|:--|
| CS-DID | Never-treated units plus eligible not-yet-treated units; controls are untreated in both comparison periods |
| Sun–Abraham | Never-treated reference group in this saturated fit |
| SDID, horizontal ridge, both multisynth fits | Never-treated donor units |
| BJS / did2s and matrix completion | All untreated observations, including pre-adoption observations of eventually treated units |
| Static and dynamic TWFE | Regression comparisons can involve already-treated units; neither is simply a never-treated/not-yet-treated design |

The following audit reruns CS-DID on the identical sample, changing only the control group and aggregation. Both use the universal period −1 baseline, no covariates, and the same country clustering. Intervals are 95% analytical influence-function intervals.

```{python}
#| echo: false
import json

audit = json.loads((OUT / "control-aggregation-audit/audit.json").read_text())
show_table(audit["comparisons"])
```

With treated-cell weights, changing from never-treated to not-yet-treated controls moves CS-DID from **−0.5138 to −0.4901**. With equal-event-time weights, the same control change moves it from **−0.4130 to −0.3930**. Holding not-yet-treated controls fixed, changing aggregation instead moves the estimate by **0.0971 SD**. Thus most of the earlier plotted separation came from weighting, not the choice of controls. The revised main comparison displays both CS-DID weighting rules. The never-treated, treated-cell CS-DID point estimate matches Sun–Abraham here, although their reported covariance calculations differ.

The **event-time panel is not composition-balanced**, despite the balanced calendar panel: event year 0 includes 529 treated units, whereas event year 45 includes only 11. Equal-event-time averaging gives these two event years equal weight; treated-cell averaging does not. No horizon-balancing restriction is applied.

These rows therefore do **not** all estimate the same scalar target. In particular, static TWFE is a coefficient benchmark, and dynamic TWFE's aggregated coefficients need not equal a clean treated-cell ATT under effect heterogeneity. Identical final averaging weights do not remove event-study contamination. The alternative estimators also use different identifying restrictions on untreated outcomes.

## Aggregate estimates

**Two CS-DID summaries are shown: equal event years 0–45 and equal treated unit-years.** The other dynamic and counterfactual summaries average effects over **17,811 treated unit-years**, giving earlier cohorts more weight because they contribute more post-treatment observations. Static TWFE is its regression coefficient. These scalar weighting rules differ; none is Peter’s original equal-weight event-years 0–20 summary.

```{python}
models, spec, estimates = results()
from report_support import estimate_rows

show_table(estimate_rows(estimates))
scalar_plot(estimates)
plt.show()
```

SDID and matrix completion are close to BJS and Sun–Abraham. Horizontal ridge gives a small positive estimate. All estimates use the same sample and outcome scale.

For both multisynth fits, the reported scalar weights each observed treated unit-year equally. The package's final `predict()` row instead averages individual post-treatment means equally, or cohort post-treatment means weighted by cohort size, without the additional post-period-count weighting; the JSON records those alternative summaries separately. The plotted event path and the reported scalar use the existing archive target weights.

Both multisynth specifications remain in this comparison: **individual −0.3574 SD; cohort-averaged −0.4733 SD**. The earlier −0.4495 SD came from separate ridge-augmented cohort fits, a different estimator not included in these two multisynth rows.

The scalar figure displays 95% intervals for all eleven summaries. Regression rows retain country-clustered CRV1 intervals; both CS-DID summaries use country-clustered influence-function intervals. SDID, matrix completion and horizontal ridge use **999 whole-country pairs-bootstrap refits**, recomputing weights, data-dependent tuning and cohort shares while holding the source outcome-SD denominator fixed; their intervals are bootstrap percentiles. Both multisynth fits use **9,999 country-clustered Mammen wild-bootstrap draws** from the package’s fixed-weight prediction procedure, with no jackknife or QP refits. A common multiplier is applied to all units within each of the 47 countries, and every draw is aggregated over the same 17,811 treated cells. These are basic bootstrap intervals conditional on fitted weights, pooling/tuning and the observed design—not full-refit pairs-bootstrap intervals. Their 95% ranges are **[−0.7055, −0.2514]** for cohort fits and **[−0.6508, −0.0838]** for individual fits.

## Event studies

### Counterfactual models

Each curve is centered on its own event-time −1 value, as in the crabbymetrics example. This is a plotting normalization only: the aggregate estimates above are not shifted. Curves show years −10 through 20; the JSON includes the complete paths and the number of contributing units at each event time.

For SDID, the package’s event summary is the raw synthetic-control gap, whereas its scalar ATT also subtracts a time-weighted pre-treatment gap. The SDID path below applies that subtraction within each cohort before aggregation. Its uncentered post-period average therefore equals the package’s ATT. Both raw and corrected gaps are retained in the JSON.

```{python}
event_plot(models, spec)
plt.show()
```

### Individual versus cohort-averaged multisynth

This isolates the two multisynth event paths. Both are centered at their own event-time −1 value for plotting only; their reported scalar ATTs are not recentered. The JSON includes both full event paths and both fits’ cohort-level summaries.

```{python}
event_plot(models, spec, only_multisynth=True)
plt.show()
```

### Regression benchmarks

These are dynamic TWFE, BJS, Sun–Abraham and CS-DID on the identical sample. Shaded bands are country-clustered 95% pointwise intervals for TWFE, Sun–Abraham and CS-DID, using each estimator’s reported covariance calculation. CS-DID now uses not-yet-treated controls, while Sun–Abraham retains its earlier specification. BJS is centered on its estimated period −1 gap and shown as a point-estimate path. CS-DID uses dashed lines and triangle markers.

```{python}
event_plot(models, spec, matrix=False)
plt.show()
```

## Numerical checks

```{python}
from validate_models import validate_models

checks = validate_models(Y, W, models)
```

- An independent closed-form calculation reproduces horizontal ridge to within $9\times10^{-11}$ SD.
- All 795 CS-DID group-time estimates match direct two-period DID contrasts to within $7\times10^{-16}$ SD. All 17,811 treated cells remain available; the dynamic scalar averages post-event times equally, while the simple scalar weights all treated cells equally. Independent influence functions, including cohort-share uncertainty, reproduce the dynamic, simple and group aggregates and country-clustered standard errors.
- Increasing SDID’s optimization budget from 1,000 to 5,000 iterations leaves the ATT unchanged.
- Matrix completion converges after 44 iterations at the reported tolerance. Tightening it from $10^{-12}$ to $10^{-14}$ changes the ATT by less than $1.6\times10^{-6}$ SD.
- Aggregation checks recover each scalar from its event path: equal post-event weights for CS-DID dynamic; treated-unit-year weights for CS-DID simple and the other new estimators. No treated cells are omitted.
- Multisynth's event-time estimates agree with `predict(fit, att=TRUE)`. A separate calendar-time reconstruction from the fitted weights and fixed-effect adjustment reproduces the reported treated-cell average. Solver statuses and balance diagnostics are included in the JSON.
- The 15-cohort multisynth fit uses the unmodified package’s native matrix construction. The 529-unit fit uses a call-scoped sparse construction of the identical optimization matrix; a four-unit check matches the unmodified package’s weights and predictions within $10^{-12}$. All solver calls converge. Both scalar estimates independently match their calendar-time reconstructions.

```{python}
#| label: export
exported = publish_results(models, spec, estimates)
print("Saved credit-markets-counterfactuals-results.json")
print("crabbymetrics:", sdid["crabbymetrics"])
print("augsynth:", multisynth["augsynth"])
print("csdid:", csdid["package"]["csdid"])
```

## Reproduction and downloads

The source uses Python **csdid 0.4.2** with **drdid 1.1.6**, crabbymetrics **0.9.1** from [commit 896ceaf](https://github.com/apoorvalal/crabbymetrics/tree/896ceaf9fa885524409ad7ddd15874034ef27033), and R augsynth **0.2.0** from [commit 7e70072](https://github.com/ebenmichael/augsynth/tree/7e70072232fe75057aa8b5d7af08c7e1aef53ae8). `uv add csdid` updates the project environment and lockfile. All code and environment files are in the [reproduction bundle](credit-markets-counterfactuals-code.zip). The source data stay in the existing local archive.

```bash
QUARTO_PYTHON=.venv/bin/python quarto render \
  analyses/credit-markets/crabbymetrics/credit-markets-counterfactuals.qmd
```

- Scalar estimates: [PNG](credit-markets-scalar.png) · [PDF](credit-markets-scalar.pdf).
- Counterfactual event studies: [PNG](credit-markets-matrix-events.png) · [PDF](credit-markets-matrix-events.pdf).
- Multisynth comparison: [PNG](credit-markets-multisynth-events.png) · [PDF](credit-markets-multisynth-events.pdf).
- Regression event studies: [PNG](credit-markets-regression-events.png) · [PDF](credit-markets-regression-events.pdf).
- [Machine-readable results](credit-markets-counterfactuals-results.json), including cohort diagnostics, settings, event paths and validation results.