1 · Estimand — which effect are we trying to learn?
This tab shows why the estimand has to come before the estimator. The same DiD design can answer different economic questions depending on whether the target is the typical unit, the population total, or a level effect. With heavy-tailed outcomes, these estimands can differ sharply, and sometimes even imply effects with opposite signs. Use the sliders below to change the data-generating process and see how the true estimands and the estimated coefficients move.
Typical-unit %
ΔΔ E[log Y]
Percent change in outcome for a representative unit, where each unit contributes equally, regardless of its size.
“Did per-unit usage rise by ~5%?”
Population-total %
ΔΔ log E[Y]
Percent change in the mean of Y — the economically relevant margin under heavy-tailed outcomes.
“Did total revenue rise by ~5%?”
Level effect
ΔΔ E[Y]
Absolute change in units of Y; the average change per unit equals the change in the mean.
“Did this add ~$2M of revenue?”
2 · DGP — Data Generating Process
The simulation creates a matched treated–control panel that simplifies the setting for educational purposes, isolating three ingredients: outcome concentration (the share of outcomes in the top decile), heterogeneous treatment effects (separate effects for the viral head and the long tail), and treatment-induced changes in variance (the proportional change in treated-unit volatility). Treated and control units have the same pre-period baseline levels by construction, so differences across estimators come from the estimand they target and how they weight the outcome distribution. For the case where treated and control units start at different levels, see the Levels-OLS bias map tab →
Presets:
Where the population total comes from
(pre-treatment distribution)
(pre-treatment distribution)
True estimand values
True typical-unit % effect
—
equal-weighted percent change across units
True population-total % effect
—
percent change in total treated outcomes
Realized top-decile share
pre-treatment top-decile share in this simulation draw
3 · Estimator — how to estimate it
Each estimator targets a different estimand and applies a different implicit or explicit weighting scheme. Log(1+Y) OLS weights observations equally, so the many low-volume units dominate the estimate. Weighted log OLS uses explicit pre-period outcome shares. PPML weights through the predicted mean and is the natural estimator for the population-total percent effect. Levels OLS targets an effect in units of the outcome. When treatment also changes outcome variance, log-transformed specifications can introduce a bias and even reverse sign. Setting Δ variance to zero isolates the role of weighting and estimand choice.
Estimated parameters relative to the true DGP estimands
Red lines show true DGP effects: solid for the population-total and dashed for the typical-unit. The black line is the point estimate from the corresponding regression, the box spans ±1 standard error (SE), and the whiskers span the 95% confidence interval (CI), clustered by matched pair. Levels OLS is reported as a percent of the treated pre-period mean, with its SE and CI scaled by the delta method.
1 · What changes when baseline levels differ?
Tab 1 assumed matched treated and control units, so baseline level differences were zero by construction. This tab asks what happens when the matched-baseline assumption is relaxed. We keep the Tab 1 DGP fixed, including the treatment effects, variance shift, outcome concentration, and panel structure. We then add two choices, defined by the sliders below: a treated-control baseline gap (b) and common proportional growth (g). Each maps to a factor shown below, and the two multiply to give the levels-OLS bias under the multiplicative trend assumed here. Under an additive trend, the bias would instead appear in PPML and log-based estimators.
Good news: with matched baseline levels between treated and control outcomes, this bias channel largely disappears; there is no level gap for proportional growth to amplify, so researchers do not need to worry about additive vs. multiplicative trend structure for credible inference.
b — baseline levels gap
treated units start higher or lower than controls; matching on levels is switched off
g — common growth
all units share the same proportional growth rate per period
b × g
—
gap factor b/(1+b)
×
—
trend-scale factor
=
—
predicted levels-OLS bias
Bias component of levels OLS (% of counterfactual mean)
Click or drag anywhere on the map to move ✕ — the sliders follow, and the simulated check below re-runs when you release. The white bands are zero bias: the entire b = 0 column and g = 0 row.
2 · Building intuition — why the two parameters multiply
Intuition for how the bias arises. Treatment effect set = 0 for instructional purposes.
Move the b and g sliders (or drag the ✕) to see how the bias changes.
Move the b and g sliders (or drag the ✕) to see how the bias changes.
Notes: The figure sets the treatment effect to zero to isolate the baseline-gap channel. The left panel shows outcomes in levels; the right panel shows the corresponding log-scale implication under common proportional growth. The shaded area is the component that levels DiD would attribute to treatment. The simulation below restores the Tab 1 treatment effects.
Under common proportional growth g, absolute changes scale with baseline levels: a treated unit (1+b) times larger mechanically adds (1+b) times more under the no-treatment counterfactual, and a levels DiD attributes the difference to treatment. With balanced windows the bias is exactly [b/(1+b)] × [(Ḡpost − Ḡpre)/Ḡpost], with Ḡ the window-average growth factor; in this 10+10 design the trend-scale factor is 1 − (1+g)−10. Either parameter at zero closes the channel: matching on baseline levels shrinks b toward zero, and no common growth means no trend-scale factor. The bias is separate from the treatment effect and adds to the scaled levels estimate. Proportional estimators (PPML, log specifications) impose proportional counterfactual trends and carry no bias from this channel. See Roth & Sant’Anna (2023) for more details.
3 · Estimator behavior
Simulated check at ✕, with Tab 1’s treatment effects
The solid red line is the true population-total estimand from Tab 1’s DGP; the dark dashed line is the levels-OLS estimate from this simulated dataset. The red shaded band shows the realized bias, the distance between the truth and the estimate. The equation card above gives the predicted bias, so any difference between the two reflects sampling noise. The black line is the point estimate from the corresponding regression, the box spans ±1 standard error (SE), and the whiskers show the 95% confidence interval (CI), clustered by matched pair. Levels OLS is reported as a percent of the counterfactual mean, with the SE and CI scaled by the delta method.
3 · Synthetic DiD in heavy tails
Synthetic difference-in-differences (Arkhangelsky et al., 2021) reweights control units (ω) and pre-periods (λ) so the control group matches the treated group’s pre-trend, then estimates the effect from a weighted two-way fixed-effects DiD. The reweighting relaxes parallel trends and absorbs level differences between the groups. We apply it to the same simulated panel as Tab 1 and Tab 2, with the treated units as the treated group and the matched control units as the donor pool. Synthetic DiD can be run on either scale: in levels it targets the population-total effect, in logs the typical unit. Levels is the more common default, though Arkhangelsky et al. themselves apply it to log GDP. In levels its period weights reduce the baseline-gap bias of Tab 2 but do not remove it, and with heavy-tailed outcomes the levels estimate is the least precise of the estimators here.
The sliders below set the outcome concentration, the baseline gap and growth that generate the levels bias, and the treatment-induced variance change that drives the log bias; every other feature of the data-generating process follows Tab 1.
The design: treated units, donor pool, synthetic control
Treated units (black) and control units (grey) begin at different levels — the baseline gap b. The reweighted donors form a synthetic control (teal) that matches the treated group through the pre-period; the gap that opens after treatment is the estimate. Unit fixed effects absorb the level difference, so the estimate comes from trends rather than levels.
Estimators against the two true estimands (with 95% CIs)
Point estimates with ±1 SE boxes and 95% confidence intervals against the two true estimands. The SDID intervals use the jackknife, which refits the weights — the default in the synthdid package; Levels OLS and PPML use the usual clustered SE, since they carry no estimated weights. PPML lands on the population total; Levels OLS is pulled toward zero by the baseline-gap × growth interaction (the Tab 2 bias); SDID levels removes part of that bias; SDID logs sits at the typical-unit estimand. SDID levels uses the same percent basis as Levels OLS, so it is unbiased when the baseline gap is zero, and the residual otherwise is the genuine b × g bias.
A confidence interval needs a standard error, and for synthetic DiD the choice matters: the estimator reweights the data before it estimates, so a standard error that treats those weights as fixed is not estimating the right object. The convention in applied work is to use one of the resampling estimators that refit the weights; the synthdid package defaults to the jackknife, which is what the SDID intervals above use. Four estimators are worth comparing:
- ω,λ fixed (clustered). The cluster-robust standard error computed once with the ω and λ weights held fixed (the “ω,λ fixed” column below). This is the baseline to be wary of for SDID: it ignores that the weights were themselves estimated, so it is not the right object. Shown for comparison, not as the default.
- Jackknife. Drop groups of units and re-run the whole procedure — refitting the weights — on each subset, then read the variance off the leave-out estimates. The package default, and the SDID interval reported in the chart above.
- Bootstrap. Resample the matched pairs with replacement, re-run on each sample, and take the spread of the estimates.
- Placebo. Reassign treatment among the control units to trace out a null distribution — built for a small treated group against a larger donor pool (Arkhangelsky et al., 2021). This design has equally many treated and control units, so to run it at the full sample we mint extra control-only units — possible only because we know the data-generating process — and draw the placebo assignments from that enlarged donor pool.
The two charts below plot the interval under all four methods — for synthetic DiD in levels and in logs — against the empirical spread of the estimate over fresh draws of the panel, the quantity a standard error is meant to recover. Both share the y-axis of the chart above.
Standard errors for the current draw.
Inference for SDID · levels
Inference for SDID · logs
Each chart fixes the SDID point estimate (levels left, logs right) and shows its 95% interval under each variance estimator, on the same axis as the chart above. The ω,λ-fixed, jackknife and bootstrap SEs all track the empirical re-draw spread closely here — on these panels the ω,λ-fixed SE is even slightly conservative, so the case against it is one of principle (it stops estimating the right object once the weights are estimated), not a large numerical gap in this example. Placebo is the exception: re-randomizing which units are “treated” makes every draw a different head-versus-tail level imbalance, so in a heavy tail it over-covers — running several times wider than the re-draw spread — and settles toward the others only as you lower the top-decile share. The real estimate uses all of the treated units at once and never incurs that imbalance, so the placebo standard error overstates its uncertainty.
Synthetic DiD improves on a plain levels DiD, but it remains a levels estimator: under the multiplicative data-generating process of Tab 1 and Tab 2 it carries a residual baseline-gap bias and, with heavy-tailed outcomes, a wider confidence interval than PPML. PPML reaches the same population-total effect more precisely. The estimand is set by the choice of outcome scale, not by the estimator; the typical-unit effect is the subject of Tab 1.
Related literature
American Economic Review · “Synthetic Difference-in-Differences”
Reweights control units (ω) and pre-periods (λ) to align pre-trends, then runs a weighted two-way-FE DiD. The estimator behind Tab 3 — and, being a levels estimator, heir to the bias and variance lessons of the other tabs.
Econometrica · “When Is Parallel Trends Sensitive to Functional Form?”
Parallel trends can hold in levels but fail in logs, or the reverse — so the estimate depends on the chosen transformation. The premise behind the Levels-OLS bias map.
Review of Economics and Statistics · “The Log of Gravity”
When the error variance depends on the mean, log-linear OLS is biased for percent effects on the mean; PPML is not. Why PPML anchors the population-total % here.
Quarterly Journal of Economics · “Logs with Zeros? Some Problems and Solutions”
With zeros, log(1+Y) effects depend on the units of Y and have no clean percent interpretation. Here zeros are rare and treatment barely moves them, so log(1+Y) remains a practical typical-unit summary.
Journal of Econometric Methods · “Dif-in-Dif Estimators of Multiplicative Treatment Effects”
If treatment changes the outcome’s variance, log-transformed DiD picks up the variance change on top of the mean effect — the Δ-variance slider channel on Tab 1.
Working paper · “On the Perils of Log Dependent Variables and Difference-in-Differences”
With baseline level differences between groups and a common trend, log DiD can flip the sign of a levels DiD; gives the condition for the flip. The b × g logic of Tab 2, in log form.
The Econometrics Journal · “Simple Approaches to Nonlinear Difference-in-Differences with Panel Data”
Poisson regression with two-way fixed effects identifies the treatment effect under an exponential-mean model. The PPML specification used throughout this companion.
1 · What this code replicates
The applet generates the data; the code estimates the specifications. The buttons below export the simulated panel for the current slider settings as a CSV file. The code below (in Python, R, and Julia) reads that file, runs the specifications from Tabs 1 and 2 — unchanged from v5. A separate, clearly-fenced block at the end of each script adds the Tab 3 synthetic-DiD estimator, so it is easy to see exactly what is new. The zip also includes a reference panel for the default settings.
2 · Download
↓ Panel (CSV) — Tab 1 settings
↓ Panel (CSV) — with Tab 2’s b and g
↓ DGP script (R)
↓ Full package (.zip)
zip: DGP generator + Python, R, and Julia scripts + reference panel + README
3 · Code
did_lab_dgp.R (base R — generates the panel CSV from the slider settings + a seed)
#!/usr/bin/env Rscript
# ============================================================================
# did_lab_dgp.R — data-generating process for the DiD Estimand Lab.
#
# Reproduces the applet's simulated panel from the slider parameters + a seed
# and writes did_lab_panel.csv, the file the estimation scripts
# (did_lab_replication.{py,R,jl}) read. It is a faithful port of the applet's
# generator: the same mulberry32 PRNG, the same Acklam inverse-normal CDF, and
# the same draw order, so for given settings it reproduces the applet's
# "Panel (CSV)" download (the RNG stream is bit-identical; y values match up to
# at most a 1-unit rounding difference from last-bit transcendental rounding
# across math libraries).
#
# Set the parameters below, then: Rscript did_lab_dgp.R
#
# The model (matching Tabs 1-2): 2*n_units units in matched pairs. Unit sizes
# are lognormal, placed at stratified quantiles so the concentration slider maps
# exactly into the realized top-decile share. Outcomes are multiplicative,
# Y = base * (1+g)^t * (1+effect*D) * eps * 100, with a pair-common shock and an
# idiosyncratic shock whose size follows a size-volatility scaling law; the
# treatment scales every treated unit's idiosyncratic volatility by (1+dvar).
# ============================================================================
## ---------------------- parameters (match the applet sliders) -------------
seed <- 4 # RNG seed
n_units <- 800 # matched pairs; the panel has 2 * n_units units
top_share_pct <- 76 # top-decile concentration, % (applet: concentration)
beta_head_pct <- -5 # treatment effect on the top-decile (head) units, % (betaHead)
beta_tail_pct <- 0 # treatment effect on the remaining (tail) units, % (betaTail)
gap_pct <- 0 # baseline gap b, treated vs control, % (Tab 2 slider; 0 for Tab 1)
growth_pct <- 0 # common growth g per period, % (Tab 2 slider; 0 for Tab 1)
dvar <- -0.35 # treatment-induced change in idiosyncratic variance (Δ variance), a fraction
out_file <- "did_lab_panel.csv"
## fixed structural constants (not exposed as sliders)
Tpre <- 10L # pre-treatment periods
Tpost <- 10L # post-treatment periods
d_cv <- 0.30 # baseline coefficient-of-variation anchor
gamma <- 0.20 # size-volatility scaling exponent
## ---------------------- exact ports of the applet's primitives ------------
# Acklam inverse-normal CDF (vectorized) — matches the applet's qnorm()
qnorm_ak <- function(p) {
a <- c(-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2,
1.383577518672690e2, -3.066479806614716e1, 2.506628277459239e0)
b <- c(-5.447609879822406e1, 1.615858368580409e2, -1.556989798598866e2,
6.680131188771972e1, -1.328068155288572e1)
cc <- c(-7.784894002430293e-3, -3.223964580411365e-1, -2.400758277161838e0,
-2.549732539343734e0, 4.374664141464968e0, 2.938163982698783e0)
d <- c( 7.784695709041462e-3, 3.224671290700398e-1, 2.445134137142996e0,
3.754408661907416e0)
pl <- 0.02425; out <- numeric(length(p))
lo <- p < pl; hi <- p > 1 - pl; mid <- !lo & !hi
if (any(lo)) { q <- sqrt(-2*log(p[lo])); out[lo] <- (((((cc[1]*q+cc[2])*q+cc[3])*q+cc[4])*q+cc[5])*q+cc[6]) / ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1) }
if (any(hi)) { q <- sqrt(-2*log(1-p[hi])); out[hi] <- -(((((cc[1]*q+cc[2])*q+cc[3])*q+cc[4])*q+cc[5])*q+cc[6]) / ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1) }
if (any(mid)) { q <- p[mid]-0.5; r <- q*q; out[mid] <- (((((a[1]*r+a[2])*r+a[3])*r+a[4])*r+a[5])*r+a[6])*q / (((((b[1]*r+b[2])*r+b[3])*r+b[4])*r+b[5])*r+1) }
out
}
# sigma such that the top size-decile holds `share` of the lognormal total (finite-N bisection)
sigma_for_share <- function(share, N) {
q <- qnorm_ak((seq_len(N) - 0.5) / N); nTop <- ceiling(N / 10)
f <- function(s) { bb <- exp(s * q); sum(bb[(N - nTop + 1):N]) / sum(bb) }
lo <- 0.05; hi <- 4
for (i in 1:40) { m <- (lo + hi) / 2; if (f(m) < share) lo <- m else hi <- m }
(lo + hi) / 2
}
# mulberry32 uniforms (vectorized). The only state update is a += C each draw,
# so the k-th draw uses state (seed + k*C) mod 2^32 — a closed form, no loop.
mulberry_uniforms <- function(seed, n) {
C <- 1831565813; M <- 4294967296 # 0x6D2B79F5, 2^32
u2i <- function(u) as.integer(ifelse(u >= 2147483648, u - M, u)) # uint32 -> signed int bit pattern
i2u <- function(i) ifelse(i < 0, as.numeric(i) + M, as.numeric(i)) # back to uint32 double
xr <- function(x, y) i2u(bitwXor(u2i(x), u2i(y)))
orr <- function(x, y) i2u(bitwOr (u2i(x), u2i(y)))
shr <- function(x, k) floor(x / (2^k)) # logical (zero-fill) right shift
imul<- function(x, y) { xh <- x %/% 65536; xl <- x %% 65536; yh <- y %/% 65536; yl <- y %% 65536
(((xh*yl + xl*yh) %% 65536) * 65536 + xl*yl) %% M } # Math.imul (low 32 bits)
a <- (seed + seq_len(n) * C) %% M
t <- imul(xr(a, shr(a, 15)), orr(1, a))
t <- xr((t + imul(xr(t, shr(t, 7)), orr(61, t))) %% M, t)
xr(t, shr(t, 14)) / M
}
## ---------------------- the data-generating process -----------------------
simulate_panel <- function() {
N <- n_units; Tt <- Tpre + Tpost
sigmaPop <- sigma_for_share(top_share_pct / 100, N)
betaHead <- beta_head_pct / 100; betaTail <- beta_tail_pct / 100
gap <- gap_pct / 100; growth <- growth_pct / 100
dVol <- sqrt(max(1 + dvar, 0.01)) - 1
sigC <- 0.12; s2C <- log(1 + sigC * sigC)
# unit sizes at stratified lognormal quantiles; head = top 10% by size
baseC <- exp(sigmaPop * qnorm_ak((seq_len(N) - 0.5) / N))
baseT <- baseC * (1 + gap)
nHead <- round(0.1 * N)
head <- logical(N); head[order(baseC, decreasing = TRUE)[seq_len(nHead)]] <- TRUE
medBase <- sort(baseC)[floor(N / 2) + 1]
theta <- ifelse(head, betaHead, betaTail)
# idiosyncratic volatility per unit: sigma_i = d * (b_i/b_med)^(-gamma), clipped
sigBase <- pmin(0.6, pmax(0.12, d_cv * (baseC / medBase)^(-gamma)))
idio0 <- sqrt(pmax(sigBase^2 - sigC^2, 4e-4)) # length N
idio1 <- idio0 * max(0.05, 1 + dVol) # treated post-period volatility
s2_0 <- log(1 + idio0^2); s2_1 <- log(1 + idio1^2)
# normal draws in the applet's consumption order: per unit, T common then
# T control-arm then T treated-arm; each normal = 2 uniforms (Box-Muller).
nU <- N * 6L * Tt
U <- mulberry_uniforms(seed, nU)
Z <- sqrt(-2 * log(U[seq(1, nU, 2)] + 1e-12)) * cos(2 * pi * U[seq(2, nU, 2)]) # N*3T normals
Zm <- matrix(Z, nrow = N, ncol = 3L * Tt, byrow = TRUE) # row i = unit i's draws
Zc <- Zm[, 1:Tt, drop = FALSE] # common shock
Za0 <- Zm[, (Tt + 1):(2 * Tt), drop = FALSE] # control arm
Za1 <- Zm[, (2 * Tt + 1):(3 * Tt), drop = FALSE] # treated arm
ln <- function(Zmat, s2) exp(sqrt(s2) * Zmat - s2 / 2) # mean-1 lognormal multiplier (s2 recycles per row)
common <- ln(Zc, s2C) # N x T, shared by both arms
eps0 <- common * ln(Za0, s2_0) # control: idio0 in every period
s2t <- cbind(matrix(s2_0, N, Tpre), matrix(s2_1, N, Tpost)) # treated: idio0 pre, idio1 post
eps1 <- common * exp(sqrt(s2t) * Za1 - s2t / 2)
gfac <- (1 + growth)^(0:(Tt - 1)) # length T
eff <- cbind(matrix(1, N, Tpre), matrix(1 + theta, N, Tpost)) # treated effect kicks in post
yc <- floor((baseC %o% gfac) * eps0 * 100 + 0.5) # control outcomes (Math.round)
yt <- floor((baseT %o% gfac) * eff * eps1 * 100 + 0.5) # treated outcomes
# assemble in the CSV order: unit ascending (controls 0..N-1, treated N..2N-1), period ascending
post <- c(rep(0L, Tpre), rep(1L, Tpost))
ctrl <- data.frame(unit = rep(0:(N - 1), each = Tt), period = rep(0:(Tt - 1), N),
treat = 0L, post = rep(post, N), y = as.integer(as.vector(t(yc))))
trt <- data.frame(unit = rep(N:(2 * N - 1), each = Tt), period = rep(0:(Tt - 1), N),
treat = 1L, post = rep(post, N), y = as.integer(as.vector(t(yt))))
df <- rbind(ctrl, trt)
df$pair <- df$unit %% N
df$D <- df$treat * df$post
df[, c("unit", "period", "pair", "treat", "post", "D", "y")]
}
df <- simulate_panel()
write.csv(df, out_file, row.names = FALSE, quote = FALSE)
cat(sprintf("wrote %s : %d rows, %d units x %d periods (seed %d, top-decile %d%%, b=%g%%, g=%g%%)\n",
out_file, nrow(df), length(unique(df$unit)), length(unique(df$period)),
seed, top_share_pct, gap_pct, growth_pct))
did_lab_replication.py (requires numpy, pandas, scipy, pyfixest, synthdid — click to view)
"""
DiD estimation script (Python), companion to Winkler et al. (2026):
https://papers.ssrn.com/sol3/papers.cfm?abstract_id=6143552
Estimates the companion's four specifications on the exported panel, all with unit
and period fixed effects and one Treat x Post regressor (D):
levels OLS reported as % of the counterfactual mean
(treated pre-mean x control post/pre ratio); delta-method CI
log(1+Y) OLS
weighted log(1+Y) OLS (pre-period mean outcome weights)
PPML reported as exp(b) - 1
Input: did_lab_panel.csv. Download it from the applet ("Panel (CSV)" buttons
on the Replication tab; any slider settings) or use the bundled copy
(default settings, seed 4). Reference values for the bundled panel:
n = 32,000 sum(y) = 23,551,023
levels -0.03861113 [-4.4695%, -3.2527%] log(1+Y) +0.00760246 [-0.3962%, +1.9167%]
w-log -0.03704146 [-4.3714%, -3.0369%] PPML -0.03859994 [-4.4669%, -3.2492%]
(package small-sample conventions can shift CI ends in the 3rd decimal)
True estimand values are properties of the DGP, not recoverable from the CSV;
read them off the applet's Tab 1 cards.
Requires: numpy, pandas, scipy, pyfixest (pip install numpy pandas scipy pyfixest)
"""
import numpy as np
import pandas as pd
import pyfixest as pf
from pyfixest.estimation import demean
from scipy.stats import t as t_dist
df = pd.read_csv("did_lab_panel.csv")
print(f"rows = {len(df)} sum(y) = {int(df.y.sum())}")
# scale for the levels estimate: treated pre-mean x control post/pre ratio
mean_treat_pre = df.loc[(df.treat == 1) & (df.post == 0), "y"].mean()
mean_ctrl_pre = df.loc[(df.treat == 0) & (df.post == 0), "y"].mean()
mean_ctrl_post = df.loc[(df.treat == 0) & (df.post == 1), "y"].mean()
scale = mean_treat_pre * mean_ctrl_post / mean_ctrl_pre
n_pairs = df["pair"].nunique()
tcrit = t_dist.ppf(0.975, df=n_pairs - 1)
df["log_y"] = np.log1p(df.y)
w_pre = df[df.post == 0].groupby("unit")["y"].mean().rename("w_pre")
df = df.merge(w_pre, on="unit")
# SEs: CRV1 clustered by matched pair. With 1:1 matching WITHOUT replacement
# (this panel), each unit belongs to exactly one pair, so two-way unit-and-pair
# clustering is identical to pair clustering. If you match WITH replacement
# (controls reused across pairs), units no longer nest in pairs, so cluster
# two-way instead: vcov={"CRV1": "unit+pair"}.
V = {"CRV1": "pair"}
m_levels = pf.feols("y ~ D | unit + period", data=df, vcov=V)
m_log = pf.feols("log_y ~ D | unit + period", data=df, vcov=V)
m_wlog = pf.feols("log_y ~ D | unit + period", data=df[df.w_pre > 0],
weights="w_pre", vcov=V) # pyfixest requires strictly
# positive weights, so drop zero-weight units (they carry
# no information either way)
m_ppml = pf.fepois("y ~ D | unit + period", data=df, vcov=V)
# (fepois may note that all-zero units were dropped for separation; this is expected)
b = lambda m: m.coef()["D"]
se = lambda m: m.tidy()["Std. Error"]["D"]
# delta-method CI for the scaled levels effect theta = b/scale: the scale is
# estimated from the same sample and strongly negatively correlated with b,
# so a fixed-scale CI would be ~2.5x too wide.
D_dem, _ = demean(df[["D"]].to_numpy(float),
df[["unit", "period"]].to_numpy(), np.ones(len(df)))
D_dem = D_dem[:, 0]
denom = D_dem @ D_dem
resid_lv = m_levels.resid()
score_pair = np.bincount(df["pair"], weights=D_dem * resid_lv,
minlength=n_pairs)
pair_mean = lambda rows: df[rows].groupby("pair")["y"].mean().to_numpy()
pm_treat_pre = pair_mean((df.treat == 1) & (df.post == 0))
pm_ctrl_pre = pair_mean((df.treat == 0) & (df.post == 0))
pm_ctrl_post = pair_mean((df.treat == 0) & (df.post == 1))
theta = b(m_levels) / scale
psi = (score_pair / (denom * scale)
- (theta / n_pairs) * (
(pm_treat_pre - pm_treat_pre.mean()) / pm_treat_pre.mean()
- (pm_ctrl_pre - pm_ctrl_pre.mean()) / pm_ctrl_pre.mean()
+ (pm_ctrl_post - pm_ctrl_post.mean()) / pm_ctrl_post.mean()))
se_theta = np.sqrt((psi**2).sum()) * n_pairs / (n_pairs - 1)
ci = lambda est, s: f"[{100*(est - tcrit*s):+.4f}%, {100*(est + tcrit*s):+.4f}%]"
print(f"levels {theta:+.8f} {ci(theta, se_theta)} (delta method)")
print(f"log1y {b(m_log):+.8f} {ci(b(m_log), se(m_log))}")
print(f"wlog {b(m_wlog):+.8f} {ci(b(m_wlog), se(m_wlog))}")
d = b(m_ppml)
print(f"ppml {np.exp(d)-1:+.8f} [{100*(np.exp(d - tcrit*se(m_ppml))-1):+.4f}%, "
f"{100*(np.exp(d + tcrit*se(m_ppml))-1):+.4f}%]")
# Standard regression table (raw coefficients).
# NOTE: levels is in outcome units here (not the scaled % above), PPML is the
# log coefficient (not exp(b)-1), and SEs use pyfixest's default convention.
print()
pf.etable([m_levels, m_log, m_wlog, m_ppml], type="md")
# ============================================================================
# Synthetic difference-in-differences (Tab 3) — via the synthdid package.
# Everything above is the v5 four-estimator script, unchanged. SDID is shown the
# way researchers run it: the published package, not a hand-rolled solver.
# pip install synthdid (Python port of Arkhangelsky et al. 2021)
# Point estimate AND jackknife SE come from the package; the jackknife refits the
# omega/lambda weights on each leave-out sample, so it accounts for the estimated
# weights — the package default, and the right object for a reweighting estimator
# (a clustered SE with the weights held fixed is not). Run in BOTH levels and logs,
# the same levels-vs-logs split as the four estimators above: levels reuses the
# counterfactual-mean `scale` already computed once for Levels OLS; logs is exp(att)-1.
# Bundled-panel reference (jackknife):
# sdid-lv -0.03784 [-5.4839%, -2.0847%] sdid-lg +0.00700 [-0.5842%, +2.0017%]
# ============================================================================
from synthdid.synthdid import Synthdid
m_sdid_lv = Synthdid(df, unit="unit", time="period", treatment="D", outcome="y").fit()
m_sdid_lv.vcov(method="jackknife") # refits the weights on each leave-out
a_lv, s_lv = m_sdid_lv.att, m_sdid_lv.se # level effect; reuse `scale` from Levels OLS
print(f"sdid-lv {a_lv/scale:+.8f} {ci(a_lv/scale, s_lv/scale)} (synthdid levels, jackknife)")
m_sdid_lg = Synthdid(df, unit="unit", time="period", treatment="D", outcome="log_y").fit()
m_sdid_lg.vcov(method="jackknife")
a_lg, s_lg = m_sdid_lg.att, m_sdid_lg.se # log points; report as exp(att)-1
print(f"sdid-lg {np.expm1(a_lg):+.8f} "
f"[{100*np.expm1(a_lg - tcrit*s_lg):+.4f}%, {100*np.expm1(a_lg + tcrit*s_lg):+.4f}%] (synthdid logs, jackknife)")
did_lab_replication.R (requires fixest, synthdid — click to view)
# ==============================================================================
# DiD estimation script (R), companion to Winkler et al. (2026):
# https://papers.ssrn.com/sol3/papers.cfm?abstract_id=6143552
#
# Estimates the companion's four specifications on the exported panel, all with unit
# and period fixed effects and one Treat x Post regressor (D):
# levels OLS reported as % of the counterfactual mean
# (treated pre-mean x control post/pre ratio); delta-method CI
# log(1+Y) OLS
# weighted log(1+Y) OLS (pre-period mean outcome weights)
# PPML reported as exp(b) - 1
#
# Input: did_lab_panel.csv. Download it from the applet ("Panel (CSV)" buttons
# on the Replication tab; any slider settings) or use the bundled copy
# (default settings, seed 4). Reference values for the bundled panel:
# n = 32,000 sum(y) = 23,551,023
# levels -0.03861113 [-4.4695%, -3.2527%] log(1+Y) +0.00760246 [-0.3962%, +1.9167%]
# w-log -0.03704146 [-4.3714%, -3.0369%] PPML -0.03859994 [-4.4669%, -3.2492%]
# (package small-sample conventions can shift CI ends in the 3rd decimal)
# True estimand values are properties of the DGP, not recoverable from the CSV;
# read them off the applet's Tab 1 cards.
#
# Requires: fixest (install.packages("fixest"))
# ==============================================================================
library(fixest)
df <- read.csv("did_lab_panel.csv")
cat(sprintf("rows = %d sum(y) = %.0f\n", nrow(df), sum(as.numeric(df$y))))
# scale for the levels estimate: treated pre-mean x control post/pre ratio
mean_treat_pre <- mean(df$y[df$treat == 1 & df$post == 0])
mean_ctrl_pre <- mean(df$y[df$treat == 0 & df$post == 0])
mean_ctrl_post <- mean(df$y[df$treat == 0 & df$post == 1])
scale <- mean_treat_pre * mean_ctrl_post / mean_ctrl_pre
n_pairs <- length(unique(df$pair))
tcrit <- qt(0.975, df = n_pairs - 1)
df$log_y <- log1p(df$y)
w_pre <- tapply(df$y[df$post == 0], df$unit[df$post == 0], mean)
df$w_pre <- as.numeric(w_pre[as.character(df$unit)])
# SEs: CRV1 clustered by matched pair. With 1:1 matching WITHOUT replacement
# (this panel), each unit belongs to exactly one pair, so two-way unit-and-pair
# clustering is identical to pair clustering. If you match WITH replacement
# (controls reused across pairs), units no longer nest in pairs, so cluster
# two-way instead: cluster = ~unit + pair.
m_levels <- feols(y ~ D | unit + period, data = df, cluster = ~pair)
m_log <- feols(log_y ~ D | unit + period, data = df, cluster = ~pair)
m_wlog <- feols(log_y ~ D | unit + period, data = df[df$w_pre > 0, ],
weights = ~w_pre, cluster = ~pair) # zero-weight units drop out anyway
m_ppml <- fepois(y ~ D | unit + period, data = df, cluster = ~pair)
# delta-method CI for the scaled levels effect theta = b/scale: the scale is
# estimated from the same sample and strongly negatively correlated with b,
# so a fixed-scale CI would be ~2.5x too wide.
D_dem <- demean(X = data.frame(D = df$D), f = df[, c("unit", "period")], tol = 1e-8)[, 1]
resid_lv <- resid(m_levels)
denom <- sum(D_dem^2)
score_pair <- tapply(D_dem * resid_lv, df$pair, sum)
pm_treat_pre <- tapply(df$y[df$treat == 1 & df$post == 0], df$pair[df$treat == 1 & df$post == 0], mean)
pm_ctrl_pre <- tapply(df$y[df$treat == 0 & df$post == 0], df$pair[df$treat == 0 & df$post == 0], mean)
pm_ctrl_post <- tapply(df$y[df$treat == 0 & df$post == 1], df$pair[df$treat == 0 & df$post == 1], mean)
theta <- unname(coef(m_levels)["D"]) / scale
psi <- score_pair / (denom * scale) -
(theta / n_pairs) * ((pm_treat_pre - mean(pm_treat_pre)) / mean(pm_treat_pre) -
(pm_ctrl_pre - mean(pm_ctrl_pre)) / mean(pm_ctrl_pre) +
(pm_ctrl_post - mean(pm_ctrl_post)) / mean(pm_ctrl_post))
se_theta <- sqrt(sum(psi^2)) * n_pairs / (n_pairs - 1)
ci <- function(est, s) sprintf("[%+.4f%%, %+.4f%%]", 100*(est - tcrit*s), 100*(est + tcrit*s))
cat(sprintf("levels %+.8f %s (delta method)\n", theta, ci(theta, se_theta)))
cat(sprintf("log1y %+.8f %s\n", coef(m_log)["D"], ci(coef(m_log)["D"], se(m_log)["D"])))
cat(sprintf("wlog %+.8f %s\n", coef(m_wlog)["D"], ci(coef(m_wlog)["D"], se(m_wlog)["D"])))
d <- coef(m_ppml)["D"]
sd_ <- se(m_ppml)["D"]
cat(sprintf("ppml %+.8f [%+.4f%%, %+.4f%%]\n", exp(d) - 1,
100*(exp(d - tcrit*sd_) - 1), 100*(exp(d + tcrit*sd_) - 1)))
# Standard regression table (raw coefficients).
# NOTE: levels is in outcome units here (not the scaled % above), PPML is the
# log coefficient (not exp(b)-1), and SEs use fixest's default convention.
cat("\n")
etable(m_levels, m_log, m_wlog, m_ppml)
# ==============================================================================
# Synthetic difference-in-differences (Tab 3) — via the synthdid package.
# Everything above is the v5 four-estimator script, unchanged. SDID is shown the
# way researchers run it: the published package, not a hand-rolled solver.
# install.packages("remotes"); remotes::install_github("synth-inference/synthdid")
# Point estimate AND jackknife SE come from the package; the jackknife refits the
# omega/lambda weights, so it accounts for the estimated weights — the package
# default, and the right object for a reweighting estimator (a clustered SE with the
# weights held fixed is not). Run in BOTH levels and logs, the same levels-vs-logs
# split as the four estimators above: levels reuses the counterfactual-mean `scale`
# already computed once for Levels OLS; logs is exp(att)-1.
# Bundled-panel reference (jackknife):
# sdid-lv -0.03784 [-5.4839%, -2.0847%] sdid-lg +0.00700 [-0.5842%, +2.0017%]
# ==============================================================================
library(synthdid)
sdid_fit <- function(outcome) { # fit + jackknife SE for one outcome column
s <- panel.matrices(df, unit = "unit", time = "period", outcome = outcome, treatment = "D")
est <- synthdid_estimate(s$Y, s$N0, s$T0)
c(att = as.numeric(est), se = sqrt(vcov(est, method = "jackknife")))
}
lv <- sdid_fit("y") # level effect; reuse `scale` from Levels OLS
cat(sprintf("sdid-lv %+.8f %s (synthdid levels, jackknife)\n",
lv[["att"]]/scale, ci(lv[["att"]]/scale, lv[["se"]]/scale)))
lg <- sdid_fit("log_y") # log points; report as exp(att)-1
cat(sprintf("sdid-lg %+.8f [%+.4f%%, %+.4f%%] (synthdid logs, jackknife)\n",
expm1(lg[["att"]]), 100*expm1(lg[["att"]] - tcrit*lg[["se"]]), 100*expm1(lg[["att"]] + tcrit*lg[["se"]])))
did_lab_replication.jl (requires DataFrames, FixedEffectModels, GLFixedEffectModels — click to view)
# ==============================================================================
# DiD estimation script (Julia), companion to Winkler et al. (2026):
# https://papers.ssrn.com/sol3/papers.cfm?abstract_id=6143552
#
# Estimates the companion's four specifications on the exported panel, all with unit
# and period fixed effects and one Treat x Post regressor (D):
# levels OLS reported as % of the counterfactual mean
# (treated pre-mean x control post/pre ratio); delta-method CI
# log(1+Y) OLS
# weighted log(1+Y) OLS (pre-period mean outcome weights)
# PPML reported as exp(b) - 1
#
# Input: did_lab_panel.csv. Download it from the applet ("Panel (CSV)" buttons
# on the Replication tab; any slider settings) or use the bundled copy
# (default settings, seed 4). Reference values for the bundled panel:
# n = 32,000 sum(y) = 23,551,023
# levels -0.03861113 [-4.4695%, -3.2527%] log(1+Y) +0.00760246 [-0.3962%, +1.9167%]
# w-log -0.03704146 [-4.3714%, -3.0369%] PPML -0.03859994 [-4.4669%, -3.2492%]
# (package small-sample conventions can shift CI ends in the 3rd decimal)
# True estimand values are properties of the DGP, not recoverable from the CSV;
# read them off the applet's Tab 1 cards.
#
# Requires: import Pkg; Pkg.add(["DelimitedFiles","DataFrames","Distributions",
# "FixedEffects","FixedEffectModels","GLFixedEffectModels","RegressionTables"])
# ==============================================================================
using DelimitedFiles, DataFrames, Statistics, Printf
using FixedEffectModels, GLFixedEffectModels
using FixedEffects: FixedEffect, solve_residuals!
using RegressionTables
using Distributions: TDist
M, hdr = readdlm("did_lab_panel.csv", ',', Int; header = true)
df = DataFrame(M, vec(Symbol.(hdr)))
@printf("rows = %d sum(y) = %d\n", nrow(df), sum(df.y))
# scale for the levels estimate: treated pre-mean x control post/pre ratio
mean_treat_pre = mean(df.y[(df.treat .== 1) .& (df.post .== 0)])
mean_ctrl_pre = mean(df.y[(df.treat .== 0) .& (df.post .== 0)])
mean_ctrl_post = mean(df.y[(df.treat .== 0) .& (df.post .== 1)])
scale = mean_treat_pre * mean_ctrl_post / mean_ctrl_pre
n_pairs = length(unique(df.pair))
tcrit = quantile(TDist(n_pairs - 1), 0.975)
df.log_y = log1p.(df.y)
w_pre = combine(groupby(df[df.post .== 0, :], :unit), :y => mean => :w_pre)
df = innerjoin(df, w_pre, on = :unit)
# SEs: CRV1 clustered by matched pair. With 1:1 matching WITHOUT replacement
# (this panel), each unit belongs to exactly one pair, so two-way unit-and-pair
# clustering is identical to pair clustering. If you match WITH replacement
# (controls reused across pairs), units no longer nest in pairs, so cluster
# two-way instead: Vcov.cluster(:unit, :pair).
m_levels = reg(df, @formula(y ~ D + fe(unit) + fe(period)), Vcov.cluster(:pair),
save = :residuals)
m_log = reg(df, @formula(log_y ~ D + fe(unit) + fe(period)), Vcov.cluster(:pair))
m_wlog = reg(df[df.w_pre .> 0, :], @formula(log_y ~ D + fe(unit) + fe(period)),
Vcov.cluster(:pair), weights = :w_pre) # zero-weight units drop out anyway
m_ppml = nlreg(df, @formula(y ~ D + fe(unit) + fe(period)),
Poisson(), LogLink(), Vcov.cluster(:pair),
separation = [:fe]) # drop all-zero units, like the other languages
bof(m) = coef(m)[findfirst(==("D"), coefnames(m))]
sof(m) = stderror(m)[findfirst(==("D"), coefnames(m))]
# delta-method CI for the scaled levels effect theta = b/scale: the scale is
# estimated from the same sample and strongly negatively correlated with b,
# so a fixed-scale CI would be ~2.5x too wide.
D_dem = solve_residuals!(Float64.(df.D),
[FixedEffect(df.unit), FixedEffect(df.period)])[1]
resid_lv = Float64.(residuals(m_levels))
denom = sum(abs2, D_dem)
df.score = D_dem .* resid_lv
score_pair = sort!(combine(groupby(df, :pair), :score => sum => :s), :pair).s
pair_mean(rows) = sort!(combine(groupby(df[rows, :], :pair), :y => mean => :m), :pair).m
pm_treat_pre = pair_mean((df.treat .== 1) .& (df.post .== 0))
pm_ctrl_pre = pair_mean((df.treat .== 0) .& (df.post .== 0))
pm_ctrl_post = pair_mean((df.treat .== 0) .& (df.post .== 1))
theta = bof(m_levels) / scale
psi = score_pair ./ (denom * scale) .-
(theta / n_pairs) .* ((pm_treat_pre .- mean(pm_treat_pre)) ./ mean(pm_treat_pre) .-
(pm_ctrl_pre .- mean(pm_ctrl_pre)) ./ mean(pm_ctrl_pre) .+
(pm_ctrl_post .- mean(pm_ctrl_post)) ./ mean(pm_ctrl_post))
se_theta = sqrt(sum(abs2, psi)) * n_pairs / (n_pairs - 1)
ci(est, s) = @sprintf("[%+.4f%%, %+.4f%%]", 100*(est - tcrit*s), 100*(est + tcrit*s))
@printf("levels %+.8f %s (delta method)\n", theta, ci(theta, se_theta))
@printf("log1y %+.8f %s\n", bof(m_log), ci(bof(m_log), sof(m_log)))
@printf("wlog %+.8f %s\n", bof(m_wlog), ci(bof(m_wlog), sof(m_wlog)))
d = bof(m_ppml)
@printf("ppml %+.8f [%+.4f%%, %+.4f%%]\n", exp(d) - 1,
100*(exp(d - tcrit*sof(m_ppml)) - 1), 100*(exp(d + tcrit*sof(m_ppml)) - 1))
# Standard regression table (raw coefficients).
# NOTE: levels is in outcome units here (not the scaled % above), PPML is the
# log coefficient (not exp(b)-1), and SEs use the package default convention.
println(regtable(m_levels, m_log, m_wlog, m_ppml))
# ============================================================================
# Synthetic difference-in-differences (Tab 3).
# No widely-used Julia SDID package exists, so researchers typically run SDID via
# the R or Python `synthdid` package (see those tabs) — that is the implementation
# to cite. The manual solver below reproduces the package point estimate (validated
# against it). Note: the SE it prints is the pair-clustered SE with the weights held
# FIXED — the baseline to be wary of for SDID; the package reports the jackknife,
# which refits the weights. Self-contained: needs just `df` and `tcrit` from above.
# Bundled-panel reference: sdid -0.03784 (point; package jackknife CI [-5.48%, -2.09%])
# ============================================================================
using LinearAlgebra: dot
function sdid_weights(A, Bv, ridge; iters = 400)
nVar = size(A, 2); w = fill(1.0 / nVar, nVar); wn2 = 1.0 / nVar
for _ in 1:iters
r = A * w .- Bv; grad = 2.0 .* (A' * r .+ ridge .* w)
km = argmin(grad); gd = grad[km] - dot(grad, w)
gd >= -1e-13 && break
Ad = A[:, km] .- (r .+ Bv); denom = 2.0 * (dot(Ad, Ad) + ridge * (1.0 - 2.0 * w[km] + wn2))
g = clamp(denom > 0 ? -gd / denom : 1.0, 0.0, 1.0)
wn2 = (1 - g)^2 * wn2 + 2 * g * (1 - g) * w[km] + g * g; w .*= (1 - g); w[km] += g
end
return w
end
units = sort(unique(df.unit)); periods = sort(unique(df.period))
Tn = length(periods); Tpre = length(unique(df.period[df.post .== 0])); Tpost = Tn - Tpre
uidx = Dict(u => i for (i, u) in enumerate(units)); pidx = Dict(p => i for (i, p) in enumerate(periods))
tdict = Dict{Int,Int}(); for row in eachrow(df); tdict[row.unit] = row.treat; end
Y = zeros(Float64, length(units), Tn)
for row in eachrow(df); Y[uidx[row.unit], pidx[row.period]] = row.y; end
treat_u = [tdict[u] for u in units]
Ytr = Y[treat_u .== 1, :]; Yco = Y[treat_u .== 0, :]; nTr = size(Ytr, 1); nCo = size(Yco, 1)
agg = vec(mean(Ytr, dims = 1)); cPreM = vec(mean(Yco[:, 1:Tpre], dims = 2))
dd = vec(diff(Yco, dims = 2)); sdc = sqrt(mean((dd .- mean(dd)).^2))
omega = sdid_weights(permutedims(Yco[:, 1:Tpre] .- cPreM), agg[1:Tpre] .- mean(agg[1:Tpre]), sqrt(nTr*Tpost)*sdc^2*Tpre)
coPostM = vec(mean(Yco[:, Tpre+1:Tn], dims = 2)); perPreM = vec(mean(Yco[:, 1:Tpre], dims = 1))
lam = sdid_weights(Yco[:, 1:Tpre] .- permutedims(perPreM), coPostM .- mean(coPostM), 1e-6*sdc^2*nCo)
co_units = units[treat_u .== 0]
unit_w = Dict(u => 1.0 / nTr for u in units); for (k, u) in enumerate(co_units); unit_w[u] = omega[k]; end
period_w = Dict(periods[t] => (t <= Tpre ? lam[t] : 1.0 / Tpost) for t in 1:Tn)
df.w_sdid = [unit_w[u] * period_w[p] for (u, p) in zip(df.unit, df.period)]
m_sdid = reg(df[df.w_sdid .> 0, :], @formula(y ~ D + fe(unit) + fe(period)), weights = :w_sdid, Vcov.cluster(:pair))
b_sdid = coef(m_sdid)[1]; se_sdid = stderror(m_sdid)[1]
# reuse the counterfactual-mean `scale` computed once above for Levels OLS
@printf("sdid %+.8f [%+.4f%%, %+.4f%%] (levels; weights-fixed SE — see R/Python for synthdid jackknife, levels & logs)\n",
b_sdid / scale, 100*(b_sdid - tcrit*se_sdid)/scale, 100*(b_sdid + tcrit*se_sdid)/scale)
DGP and specification details
Balanced matched-pair panel. Each treated unit i has baseline b0i(1+b); its matched control has b0i, where b is the Tab 2 baseline-gap slider (b = 0 on Tab 1). The b0i sit at the midpoint quantiles of lognormal(0, σ²pop), and σpop is solved numerically on that quantile grid so the realized top-decile share of the total equals the “% of outcome in top decile” slider. Panel: 10 pre + 10 post periods. Outcomes: Yit = b0i·(1+g)t·(1+θi)Dit·εit·100, rounded to integer counts, with g the Tab 2 common-growth slider (g = 0 on Tab 1). θi equals the head slider for the top size decile and the tail slider for the rest. εit is mean-one multiplicative lognormal noise, eσZ−σ²/2, composed of a pair-common component (cv 0.12, shared by the two units of a pair in each period) and a unit-idiosyncratic component. Baseline total cv: σi = 0.3·(b0i/b0,med)−0.2, clipped to [0.12, 0.60]. Treatment multiplies the idiosyncratic component of treated-unit volatility by √(1+Δvar) in post periods, for all treated units; the pair-common component is unchanged. True estimands: typical-unit = equal-weighted mean of θi; population-total = baseline-weighted mean of θi. Specifications: every estimator uses unit and period fixed effects with a single Treat×Post regressor — OLS in levels, OLS on log(1+Y), weighted OLS on log(1+Y) with pre-period mean outcome weights, and Poisson pseudo-maximum-likelihood (log link, reported as eδ−1). Standard errors are CRV1, clustered by matched pair, with t critical values on pairs−1 degrees of freedom; two-way clustering by unit and pair is identical here because units are nested within pairs. The levels coefficient is reported as a percent of the treated pre-period mean multiplied by the control post/pre growth ratio, with a delta-method confidence interval. The Tab 2 bias surface is the exact probability limit [b/(1+b)] · [(Ḡpost−Ḡpre)/Ḡpost], with Ḡ the within-window mean of (1+g)t.