From 100cc7e420afeb8c41a0f46a968cd07c9cde3089 Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 29 Jun 2026 09:15:06 -0400 Subject: [PATCH] test(lpdid): add R-parity validation harness (Dube et al. 2025), Phase B2 Pin the absorbing LPDiD estimator against the method authors' own R recipes (danielegirardi/lpdid) with an alexCardazzi/lpdid cross-check gate: - benchmarks/R/generate_lpdid_golden.R: in-R panel (+ interior-gap unit) and 6 variants (variance-weighted, reweight, pmd, direct-covariate, pooled, RA-point); writes committed lpdid_test_panel.csv + lpdid_golden.json. - tests/test_methodology_lpdid.py: skip-guarded parity (att/se to ~1e-12, cross-platform asserted at 1e-6/1e-7). - benchmarks/python/coverage_lpdid_ra.py + lpdid_ra_coverage.json: ungated Monte-Carlo study validating the RA influence-function SE calibration (~0.95). Resolves the two provisional REGISTRY deviation notes in the library's favour with no estimator change: the RA SE matches the Stata teffects convention (point-anchored, SE pinned + coverage-validated; no R-package analogue), and the pooled estimand matches the authors' fixed-composition recipe (correcting the prior "horizon-stacked" wording). no_composition documented as more paper-faithful than the R packages (B1-tested). Ticks the REGISTRY B2 checklist box. Co-Authored-By: Claude Opus 4.8 (1M context) --- CHANGELOG.md | 15 + TODO.md | 1 + benchmarks/R/generate_lpdid_golden.R | 314 +++++++++++ benchmarks/data/lpdid_golden.json | 76 +++ benchmarks/data/lpdid_ra_coverage.json | 76 +++ benchmarks/data/lpdid_test_panel.csv | 720 +++++++++++++++++++++++++ benchmarks/python/coverage_lpdid_ra.py | 216 ++++++++ docs/methodology/REGISTRY.md | 7 +- tests/test_methodology_lpdid.py | 228 ++++++++ 9 files changed, 1650 insertions(+), 3 deletions(-) create mode 100644 benchmarks/R/generate_lpdid_golden.R create mode 100644 benchmarks/data/lpdid_golden.json create mode 100644 benchmarks/data/lpdid_ra_coverage.json create mode 100644 benchmarks/data/lpdid_test_panel.csv create mode 100644 benchmarks/python/coverage_lpdid_ra.py create mode 100644 tests/test_methodology_lpdid.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 381d64f89..874e7e4fa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **`LPDiD` R-parity validation (absorbing).** `tests/test_methodology_lpdid.py` pins the + estimator against the method authors' own R recipes (`danielegirardi/lpdid` event-study / + reweight / premean / pooled `fixest::feols` specifications) with an `alexCardazzi/lpdid` + cross-check gate on the event-study variants, generated by + `benchmarks/R/generate_lpdid_golden.R`. Variance-weighted, reweighted, premean, + direct-covariate, pooled, and regression-adjustment-point estimates match to ~1e-12 + (pooled is anchored to the authors' fixed-composition recipe only — `alexCardazzi`'s pooled + uses a laxer clean-control window and is recorded for transparency, not gated). The + regression-adjustment standard error (canonically Stata `teffects ra` only — no R-package + analogue) is pinned and its calibration validated by an ungated Monte-Carlo coverage study + (`benchmarks/python/coverage_lpdid_ra.py`, ≈0.95 coverage across cluster counts). Confirms + the two previously-provisional REGISTRY notes (RA influence-function variance convention; + pooled fixed-composition estimand) resolve in the library's favour with no estimator change. + ## [3.6.0] - 2026-06-29 ### Added diff --git a/TODO.md b/TODO.md index bcaea1561..a4b0835f5 100644 --- a/TODO.md +++ b/TODO.md @@ -116,6 +116,7 @@ exists but parity can't be verified without a local toolchain. | `CallawaySantAnna` bootstrap: align p-value computation with R `did`'s symmetric-percentile method (former "CallawaySantAnna Bootstrap Improvements" section). | `staggered.py` | — | Low | | **`bias_corrected_local_linear` (lprobust) Phase-1c follow-ups:** extend golden parity to `kernel ∈ {triangular, uniform}` (epa-only today); expose `vce ∈ {hc0,hc1,hc2,hc3}` on the public wrapper once R goldens exist (port supports all four; needs a per-mode generator + a hc2/hc3 q-fit-leverage decision); clustered-DGP auto-bandwidth parity is **blocked upstream** on an nprobust singleton-cluster bug in `lpbwselect.mse.dpi` (Phase-1c DGP 4 uses manual `h=b=0.3`). | `_nprobust_port.py`, `local_linear.py`, `generate_nprobust_lprobust_golden.R` | Phase 1c | Low-Med | | `HeterogeneousAdoptionDiD` Stute-family Stata-bridge parity: no public R `Stutetest` package exists; would add `benchmarks/stata/generate_stute_golden.do` + a Stata dependency. | `benchmarks/stata/`, `tests/test_stute_test_parity.py` | follow-up | Low | +| **`LPDiD` regression-adjustment SE — no runnable R reference.** The RA influence-function cluster SE is canonically Stata `teffects ra ... atet vce(cluster)` only; no R package computes it (`alexCardazzi/lpdid` does direct covariate inclusion, not RA). Today the RA *point* is R-anchored (~1e-12), the SE is pinned + MC-coverage-validated (`coverage_lpdid_ra.py`). Follow-up: contribute the RA path to `alexCardazzi/lpdid` so a runnable R RA reference exists — only a *trusted* anchor once cross-checked vs Stata `teffects` (else circular). | `tests/test_methodology_lpdid.py`, `benchmarks/python/coverage_lpdid_ra.py` | #B2 follow-up | Low | | `HeterogeneousAdoptionDiD` Phase-3 R-parity: ships coverage-rate validation on synthetic DGPs, not tight point parity vs `chaisemartin::stute_test` / `yatchew_test` (needs bootstrap-seed-semantics + `B` alignment across numpy/R). | `tests/test_had_pretests.py` | Phase 3 | Low | ### Parked — pending user demand / out of scope diff --git a/benchmarks/R/generate_lpdid_golden.R b/benchmarks/R/generate_lpdid_golden.R new file mode 100644 index 000000000..25417f0ae --- /dev/null +++ b/benchmarks/R/generate_lpdid_golden.R @@ -0,0 +1,314 @@ +#!/usr/bin/env Rscript +# Golden generator: LPDiD vs the method authors' reference recipes +# (Dube, Girardi, Jorda & Taylor 2025), with an alexCardazzi/lpdid cross-check. +# +# The paper specifies NO standard-error formula (Section 1 defers to "standard +# techniques"), so every estimate is validated against the reference *software*: +# - PRIMARY: the authors' own R recipe (danielegirardi/lpdid example scripts) - +# fixest::feols(long_diff ~ treat_switch | time, vcov = ~unit) with the +# reghdfe-style small-sample correction setFixest_ssc(adj, cluster.adj), plus +# the reweighted / premean / pooled / regression-adjustment variants. +# - CROSS-CHECK GATE: alexCardazzi::lpdid() (third-party, pinned commit) must +# agree where conditions match. The POOLED cross-check is intentionally skipped +# (alexCardazzi uses a laxer pooled clean-control window than the authors' +# clean-through-window-end recipe); alex's pooled value is recorded in `meta`. +# +# Regression-adjustment (RA, variant ra_cov): the canonical RA standard error is +# Stata `teffects ra ... atet vce(cluster)` only - no R package implements it +# (alexCardazzi uses direct covariate inclusion, not RA). We therefore anchor the +# RA POINT estimate here (full-interaction route == teffects point) and pin the +# library's influence-function SE on the Python side as a documented regression +# value; its calibration is validated out-of-band by the ungated Monte-Carlo +# coverage study benchmarks/python/coverage_lpdid_ra.py. See REGISTRY.md "## LPDiD". +# +# Outputs (checked into the repo): +# benchmarks/data/lpdid_test_panel.csv +# benchmarks/data/lpdid_golden.json +# +# Usage: +# Rscript benchmarks/R/generate_lpdid_golden.R + +suppressMessages({ + library(dplyr); library(tsibble); library(slider); library(fixest) + library(sandwich); library(data.table); library(jsonlite); library(lpdid) +}) +setFixest_ssc(ssc(adj = TRUE, cluster.adj = TRUE)) # authors' reghdfe-style convention + +SEED <- 20260629L +set.seed(SEED) +PRE <- 3L; POST <- 4L +ALEX_SHA <- "ba64563983861be5e616f6020842c1a1cdf17a27" +TOL_XCHECK <- 1e-8 # authors' recipe vs alexCardazzi agreement gate + +# Verify the installed alexCardazzi/lpdid is exactly the pinned commit, so the +# cross-check gate is meaningful (fail-closed, not just metadata). +.alex_sha <- packageDescription("lpdid")$RemoteSha +if (is.null(.alex_sha) || is.na(.alex_sha)) { + stop(sprintf("installed `lpdid` has no RemoteSha; reinstall via remotes::install_github(\"alexCardazzi/lpdid\", ref = \"%s\")", + ALEX_SHA)) +} +if (.alex_sha != ALEX_SHA) { + stop(sprintf("installed `lpdid` commit %s != pinned ALEX_SHA %s; reinstall the pinned commit", + .alex_sha, ALEX_SHA)) +} + +# ============================================================ +# 1. Panel: staggered absorbing, ~63% never-treated, heterogeneous + persistent, +# exogenous confounder x, plus ONE interior-gap unit (pins the reindex / Dev 5). +# ============================================================ +Tt <- 12L +cohort <- c(rep(4L, 12), rep(8L, 10), rep(0L, 38)) # 1-12 @ t=4, 13-22 @ t=8, 23-60 never +units <- seq_along(cohort) +unit_fe <- rnorm(length(units), 0, 2) +time_fe <- rnorm(Tt, 0, 1) +rows <- vector("list", length(units) * Tt); idx <- 0L +for (u in units) { + g <- cohort[u]; eps_prev <- 0 + for (t in 1:Tt) { + eps <- 0.4 * eps_prev + rnorm(1); eps_prev <- eps + treated <- as.integer(g > 0 && t >= g) + eff <- if (treated == 1) { + k <- t - g; base <- if (g == 4) 2.0 else 3.5; slope <- if (g == 4) 0.5 else 0.2 + base + slope * k + } else 0 + x <- rnorm(1) + y <- unit_fe[u] + time_fe[t] + eff + 0.7 * x + eps + idx <- idx + 1L + rows[[idx]] <- list(unit = u, time = t, treat = treated, y = y, x = x) + } +} +panel <- rbindlist(rows) +panel <- panel[!(unit == 60L & time == 7L)] # interior gap on a never-treated unit +setorder(panel, unit, time) +panel_path <- file.path("benchmarks", "data", "lpdid_test_panel.csv") +dir.create(dirname(panel_path), recursive = TRUE, showWarnings = FALSE) +fwrite(panel, panel_path) +message(sprintf("Wrote panel: %s (%d rows, %d units, 1 interior gap)", + panel_path, nrow(panel), uniqueN(panel$unit))) + +# ============================================================ +# 2. Prep on the interior-calendar grid (fill_gaps == our reindex), absorbing-fill +# treatment on the grid, features computed on the grid, restrict back to observed. +# ============================================================ +td <- panel[treat == 1, .(treat_date = min(time)), by = unit] +maxT <- max(panel$time) +prep <- panel %>% + left_join(td, by = "unit") %>% + as_tsibble(key = unit, index = time) %>% + fill_gaps() %>% # interior gaps -> explicit NA rows + group_by(unit) %>% arrange(time, .by_group = TRUE) %>% + mutate(treat_date = if (all(is.na(treat_date))) NA_real_ else suppressWarnings(min(treat_date, na.rm = TRUE))) %>% + mutate( + treat = as.integer(!is.na(treat_date) & time >= treat_date), # absorbing-fill on grid + tdiff = treat - lag(treat, 1L), + Ly = lag(y, 1L), + premean2 = (lag(y, 1L) + lag(y, 2L)) / 2, # fixed-k=2 premean (pmd variant) + obs = !is.na(y) + ) %>% ungroup() %>% as_tibble() +prep$tdiff[is.na(prep$tdiff) | prep$tdiff < 0] <- 0 + +# clean sample for one horizon; restrict to OBSERVED rows with a finite long diff. +clean_h <- function(h, base, kind) { + p <- prep %>% group_by(unit) %>% arrange(time, .by_group = TRUE) + if (kind == "post") { + p <- p %>% mutate(Dy = lead(y, h) - .data[[base]], Fh = lead(treat, h)) %>% ungroup() %>% + filter(obs & !is.na(Dy) & !is.na(tdiff) & !is.na(Fh) & (tdiff == 1 | Fh == 0)) + } else { + p <- p %>% mutate(Dy = lag(y, h) - .data[[base]]) %>% ungroup() %>% + filter(obs & !is.na(Dy) & !is.na(tdiff) & !is.na(treat) & (tdiff == 1 | treat == 0)) + } + p +} +# identified == both treated-switch and clean-control present at >1 time level +identified <- function(d) nrow(d) > 0 && length(unique(d$time[d$tdiff == 1])) > 0 && + any(d$tdiff == 0) + +feols_es <- function(formula, weights = NULL) function(d) { + if (!identified(d)) return(c(NA_real_, NA_real_)) + m <- if (is.null(weights)) feols(formula, data = d, vcov = ~unit, warn = FALSE, notes = FALSE) + else feols(formula, data = d, weights = weights, vcov = ~unit, warn = FALSE, notes = FALSE) + c(coef(m)[["tdiff"]], se(m)[["tdiff"]]) +} + +post_h <- 0:POST +pre_h <- 2:PRE # h = -2 .. -PRE placebos (h = -1 is the fixed reference) + +run_es <- function(base, fml) { + out <- list() + for (h in post_h) out[[as.character(h)]] <- feols_es(fml)(clean_h(h, base, "post")) + for (h in pre_h) out[[as.character(-h)]] <- feols_es(fml)(clean_h(h, base, "pre")) + out +} + +# ---- girardi reweighting (FWL inverse weights) for the equally-weighted ATT ---- +reweighted_es <- function() { + get_w <- function(d) { + mw <- feols(tdiff ~ 1 | time, data = d, vcov = "iid", warn = FALSE, notes = FALSE) + d$nw <- residuals(mw); d$nw[d$tdiff != 1] <- NA_real_ + den <- sum(d$nw, na.rm = TRUE) + d <- d %>% group_by(time) %>% + mutate(w = nw / den, gw = suppressWarnings(max(w, na.rm = TRUE)), + gw = ifelse(is.infinite(gw), NA_real_, gw), + w = ifelse(is.na(w), gw, w)) %>% ungroup() + 1 / d$w + } + out <- list() + w0 <- NULL + for (h in post_h) { + d <- clean_h(h, "Ly", "post") + if (!identified(d)) { out[[as.character(h)]] <- c(NA_real_, NA_real_); next } + d$rw <- get_w(d); if (h == 0) w0 <- d %>% select(unit, time, rw0 = rw) + d <- d %>% filter(rw > 0) + out[[as.character(h)]] <- feols_es(Dy ~ tdiff | time, weights = ~rw)(d) + } + for (h in pre_h) { # pre uses the h=0 weights (reweight_0) + d <- clean_h(h, "Ly", "pre") %>% left_join(w0, by = c("unit", "time")) + d <- d %>% filter(!is.na(rw0) & rw0 > 0) + out[[as.character(-h)]] <- feols_es(Dy ~ tdiff | time, weights = ~rw0)(d) + } + out +} + +# ---- RA via full-interaction route (teffects ra ATET POINT; SE pinned in Python) ---- +ra_point_es <- function() { + atet_h <- function(d) { + if (!identified(d)) return(c(NA_real_, NA_real_, NA_real_)) + d2 <- d %>% transmute(Dy, dtr = factor(tdiff), x, tf = factor(time), unit) + fi <- lm(Dy ~ dtr * (tf + x), data = d2) + b <- coef(fi); keep <- !is.na(b); nm <- names(b)[keep]; bk <- b[keep] + mm <- model.matrix(fi)[, nm, drop = FALSE] + V0 <- sandwich::vcovCL(fi, cluster = d2$unit, cadjust = FALSE, type = "HC0") + g <- colMeans(mm[d2$dtr == "1", , drop = FALSE]); g[!grepl("^dtr1", nm)] <- 0 + c(sum(g * bk), sqrt(as.numeric(t(g) %*% V0 %*% g)), NA_real_) # [att, conditional CR0 se (ref only)] + } + out <- list() + for (h in post_h) out[[as.character(h)]] <- atet_h(clean_h(h, "Ly", "post")) + for (h in pre_h) out[[as.character(-h)]] <- atet_h(clean_h(h, "Ly", "pre")) + out +} + +vw_es <- run_es("Ly", Dy ~ tdiff | time) +direct_es <- run_es("Ly", Dy ~ tdiff + x | time) +pmd_es <- run_es("premean2", Dy ~ tdiff | time) +ew_es <- reweighted_es() +ra_es <- ra_point_es() + +# ---- pooled (authors' recipe): mean(y over [t, t+POST]) - Ly, clean-through-F.POST ---- +pooled_tbl <- prep %>% group_by(unit) %>% arrange(time, .by_group = TRUE) %>% + mutate(ypost = slide_dbl(y, ~ mean(.x, na.rm = FALSE), .before = 0, .after = POST, .complete = TRUE), + ypre = slide_dbl(y, ~ mean(.x, na.rm = FALSE), .before = PRE, .after = -2, .complete = TRUE), + pooly_post = ypost - Ly, pooly_pre = ypre - Ly, + Fpost = lead(treat, POST)) %>% ungroup() +pp <- pooled_tbl %>% filter(obs & !is.na(pooly_post) & !is.na(tdiff) & !is.na(Fpost) & (tdiff == 1 | Fpost == 0)) +mpost <- feols(pooly_post ~ tdiff | time, data = pp, vcov = ~unit, warn = FALSE, notes = FALSE) +ppre <- pooled_tbl %>% filter(obs & !is.na(pooly_pre) & !is.na(tdiff) & !is.na(treat) & (tdiff == 1 | treat == 0)) +mpre <- feols(pooly_pre ~ tdiff | time, data = ppre, vcov = ~unit, warn = FALSE, notes = FALSE) +pooled <- list(post = c(coef(mpost)[["tdiff"]], se(mpost)[["tdiff"]]), + pre = c(coef(mpre)[["tdiff"]], se(mpre)[["tdiff"]])) + +# ============================================================ +# 3. alexCardazzi::lpdid cross-check (VW / direct / pmd / nocomp event-study + pooled-in-meta) +# ============================================================ +alex_df <- as.data.frame(panel) +alex_fit <- function(...) tryCatch( + lpdid(alex_df, window = c(-PRE, POST), y = "y", unit_index = "unit", + time_index = "time", treat_status = "treat", ...), + error = function(e) list(error = conditionMessage(e))) +alex_es <- function(fit) { + if (!is.null(fit$error)) return(list(error = fit$error)) + setNames(lapply(seq_along(fit$window), + function(i) c(fit$coeftable$Estimate[i], fit$coeftable$`Std. Error`[i])), + as.character(fit$window)) +} +alex_vw <- alex_es(alex_fit()) +alex_direct <- alex_es(alex_fit(controls = ~ x)) +alex_pmd <- alex_es(alex_fit(pmd = TRUE, pmd_lag = 2)) +alex_nocomp <- alex_es(alex_fit(nocomp = TRUE)) +alex_pooled <- tryCatch({ + f <- lpdid(alex_df, window = c(-PRE, POST), y = "y", unit_index = "unit", + time_index = "time", treat_status = "treat", pooled = TRUE) + c(f$coeftable$Estimate[1], f$coeftable$`Std. Error`[1]) +}, error = function(e) c(NA_real_, NA_real_)) + +# cross-check gate: authors' recipe must agree with alexCardazzi where comparable. +# Fail-closed - a cross-check that cannot run (alex error, no overlapping horizons, +# all-NA comparisons, or non-finite difference) is a FAILURE, not a silent pass, so +# the gate can never be bypassed into writing ungated goldens. +xcheck <- function(name, ours, theirs) { + if (!is.null(theirs$error)) { + stop(sprintf("alexCardazzi cross-check UNAVAILABLE for %s: %s", name, theirs$error)) + } + hs <- intersect(names(ours), names(theirs)) + if (length(hs) == 0L) { + stop(sprintf("alexCardazzi cross-check for %s: no overlapping horizons", name)) + } + diffs <- sapply(hs, function(h) { + a <- ours[[h]]; b <- theirs[[h]] + if (any(is.na(a)) || any(is.na(b))) NA_real_ else max(abs(a[1:2] - b[1:2])) + }) + if (all(is.na(diffs))) { + stop(sprintf("alexCardazzi cross-check for %s: all overlapping horizons non-comparable (NA)", name)) + } + d <- max(diffs, na.rm = TRUE) + message(sprintf(" [xcheck %s] max|authors - alex| = %.2e over h={%s}", + name, d, paste(hs, collapse = ","))) + if (!is.finite(d) || d > TOL_XCHECK) { + stop(sprintf("alexCardazzi cross-check FAILED for %s (%.2e > %.0e)", name, d, TOL_XCHECK)) + } +} +message("alexCardazzi cross-check (authors' feols recipe vs the package):") +xcheck("vw", vw_es, alex_vw) +xcheck("direct", direct_es, alex_direct) +xcheck("pmd", pmd_es, alex_pmd) +# nocomp is intentionally NOT a golden parity variant: the library's no_composition +# fixes the realized sample across ALL post horizons (the paper's fixed-composition +# intent, Section 3.6) and excludes cohorts with p_g > T-H, whereas alexCardazzi uses a +# looser per-horizon sample and a stricter treat_date < maxT-POST cutoff - so our path is +# MORE faithful to the paper but has no matching R-package reference. It is validated by +# the pure-Python B1 tests in tests/test_lpdid.py. alexCardazzi's looser-semantics value +# is recorded in meta for transparency only. +message(" [nocomp] no golden parity variant (library is more paper-faithful than any R pkg;") +message(" B1-tested; alexCardazzi looser-semantics value recorded in meta).") + +# ============================================================ +# 4. Write golden JSON +# ============================================================ +golden <- list( + meta = list( + estimator = "LPDiD (Dube, Girardi, Jorda & Taylor 2025) - absorbing", + r_version = R.version.string, + fixest_version = as.character(packageVersion("fixest")), + lpdid_alexcardazzi_version = as.character(packageVersion("lpdid")), + lpdid_alexcardazzi_commit = ALEX_SHA, + seed = SEED, pre_window = PRE, post_window = POST, + se_convention = paste( + "Default/weighted/direct/pmd/pooled SEs: feols cluster-robust at unit with the", + "reghdfe-style finite-sample factor (G/(G-1))*((n-1)/(n-k)) via", + "setFixest_ssc(adj=TRUE, cluster.adj=TRUE) - the authors' convention."), + ra_se_note = paste( + "ra_cov records the RA POINT estimate (full-interaction == teffects ra atet).", + "The canonical RA SE is Stata teffects (unconditional IF, NO finite-sample", + "factor) - no R package computes it; the library influence-function SE is pinned", + "on the Python side and calibration-validated by benchmarks/python/coverage_lpdid_ra.py.", + "ra_cov[h] = c(att, conditional_CR0_se_for_reference_only)."), + pooled_note = paste( + "Pooled uses the authors' fixed-composition window-mean recipe", + "(mean(y over [t,t+POST]) - y_{t-1}, clean through F.POST). alexCardazzi's pooled", + "uses a laxer clean-control window, so its pooled value is recorded here for", + "transparency but NOT used as the cross-check gate."), + nocomp_note = paste( + "no_composition is NOT a golden parity variant: the library fixes the realized", + "sample across all post horizons (paper Section 3.6) and excludes cohorts p_g > T-H,", + "more faithful to the paper than alexCardazzi's looser per-horizon version (recorded", + "below). The library path is validated by the pure-Python B1 tests."), + alex_pooled_post = alex_pooled, + alex_nocomp_looser = alex_nocomp + ), + vw_es = vw_es, ew_es = ew_es, pmd_es = pmd_es, + direct_es = direct_es, ra_cov = ra_es, pooled = pooled +) +golden_path <- file.path("benchmarks", "data", "lpdid_golden.json") +# na = "null" so R NA serializes as JSON null (not the string "NA"), which the +# Python loader reads as None and handles in its missing-value branch. +write_json(golden, golden_path, auto_unbox = TRUE, pretty = TRUE, digits = 12, na = "null") +message(sprintf("Wrote golden: %s", golden_path)) diff --git a/benchmarks/data/lpdid_golden.json b/benchmarks/data/lpdid_golden.json new file mode 100644 index 000000000..35db69fdf --- /dev/null +++ b/benchmarks/data/lpdid_golden.json @@ -0,0 +1,76 @@ +{ + "meta": { + "estimator": "LPDiD (Dube, Girardi, Jorda & Taylor 2025) - absorbing", + "r_version": "R version 4.5.2 (2025-10-31)", + "fixest_version": "0.13.2", + "lpdid_alexcardazzi_version": "0.6.1", + "lpdid_alexcardazzi_commit": "ba64563983861be5e616f6020842c1a1cdf17a27", + "seed": 20260629, + "pre_window": 3, + "post_window": 4, + "se_convention": "Default/weighted/direct/pmd/pooled SEs: feols cluster-robust at unit with the reghdfe-style finite-sample factor (G/(G-1))*((n-1)/(n-k)) via setFixest_ssc(adj=TRUE, cluster.adj=TRUE) - the authors' convention.", + "ra_se_note": "ra_cov records the RA POINT estimate (full-interaction == teffects ra atet). The canonical RA SE is Stata teffects (unconditional IF, NO finite-sample factor) - no R package computes it; the library influence-function SE is pinned on the Python side and calibration-validated by benchmarks/python/coverage_lpdid_ra.py. ra_cov[h] = c(att, conditional_CR0_se_for_reference_only).", + "pooled_note": "Pooled uses the authors' fixed-composition window-mean recipe (mean(y over [t,t+POST]) - y_{t-1}, clean through F.POST). alexCardazzi's pooled uses a laxer clean-control window, so its pooled value is recorded here for transparency but NOT used as the cross-check gate.", + "nocomp_note": "no_composition is NOT a golden parity variant: the library fixes the realized sample across all post horizons (paper Section 3.6) and excludes cohorts p_g > T-H, more faithful to the paper than alexCardazzi's looser per-horizon version (recorded below). The library path is validated by the pure-Python B1 tests.", + "alex_pooled_post": [3.487463576882, 0.3495904202394], + "alex_nocomp_looser": { + "-3": [-0.7532927692809, 0.5003193796022], + "-2": [-0.250198606317, 0.6160670067758], + "-1": [0, 0], + "0": [1.46387943226, 0.545454853302], + "1": [2.335727820716, 0.5511865769279], + "2": [2.865671038762, 0.5566693486612], + "3": [3.245223911323, 0.4728816808867], + "4": [4.100734546643, 0.4370902112454] + } + }, + "vw_es": { + "0": [2.891454184758, 0.4972838592697], + "1": [3.215287698688, 0.4242105863868], + "2": [3.570283061763, 0.4502488079078], + "3": [3.968564586288, 0.3904079437881], + "4": [4.269176444588, 0.396524220757], + "-2": [0.1183626967809, 0.4413440724331], + "-3": [-0.1254433057345, 0.4292726043484] + }, + "ew_es": { + "0": [2.903628272087, 0.4973282590126], + "1": [3.221973963653, 0.4235473199649], + "2": [3.57527972142, 0.4499773773075], + "3": [3.971884855509, 0.390089371568], + "4": [4.265997983221, 0.3934421813065], + "-2": [0.1195811966608, 0.4402860720602], + "-3": [-0.1219233071656, 0.4298940947632] + }, + "pmd_es": { + "0": [2.832272836367, 0.4337549369753], + "1": [3.156106350297, 0.3522180898809], + "2": [3.511101713373, 0.394076563941], + "3": [3.902081106111, 0.3654320104818], + "4": [4.269986848274, 0.3537196151693], + "-2": [0.05918134839047, 0.2206720362165], + "-3": [-0.184624654125, 0.3713701578016] + }, + "direct_es": { + "0": [2.876086144279, 0.4110153204559], + "1": [3.217671119032, 0.433940624349], + "2": [3.571778263524, 0.4571207641365], + "3": [3.969558831006, 0.3912243590417], + "4": [4.268878161094, 0.3995261701574], + "-2": [0.1185563879358, 0.442300060411], + "-3": [-0.1239894706087, 0.4328879644821] + }, + "ra_cov": { + "0": [2.888871838355, 0.2774049784919, null], + "1": [3.22526891068, 0.3671167840304, null], + "2": [3.577950221628, 0.3920693801281, null], + "3": [3.973278128237, 0.363703531293, null], + "4": [4.265103781445, 0.3614251956151, null], + "-2": [0.1204444745696, 0.4103497220604, null], + "-3": [-0.1197046094519, 0.4026535247538, null] + }, + "pooled": { + "post": [3.532265685012, 0.3521082333958], + "pre": [-0.003540304476795, 0.377308359478] + } +} diff --git a/benchmarks/data/lpdid_ra_coverage.json b/benchmarks/data/lpdid_ra_coverage.json new file mode 100644 index 000000000..a1a5bbb85 --- /dev/null +++ b/benchmarks/data/lpdid_ra_coverage.json @@ -0,0 +1,76 @@ +{ + "study": "LPDiD RA influence-function SE coverage", + "nominal": 0.95, + "tau": { + "0": 1.0, + "1": 1.5, + "2": 2.0, + "3": 2.5 + }, + "note": "RA IF variance is asymptotic (teffects convention, no finite-sample factor); empirical coverage of the true effect holds near nominal across G for both the event-study horizons and the headline pooled-row RA CI - the t(G-1) reference keeps it well-calibrated even at modest G. Underwrites REGISTRY ## LPDiD Deviation 2.", + "sweep": [ + { + "n_units": 30, + "reps": 500, + "per_horizon": { + "0": 0.954, + "1": 0.952, + "2": 0.94, + "3": 0.968 + }, + "mean_event_coverage": 0.9535, + "pooled_att_coverage": 0.968, + "n_valid": { + "0": 500, + "1": 500, + "2": 500, + "3": 500 + }, + "n_pooled_valid": 500, + "n_fit_errors": 0, + "min_valid_share": 1.0 + }, + { + "n_units": 100, + "reps": 500, + "per_horizon": { + "0": 0.942, + "1": 0.932, + "2": 0.944, + "3": 0.942 + }, + "mean_event_coverage": 0.94, + "pooled_att_coverage": 0.946, + "n_valid": { + "0": 500, + "1": 500, + "2": 500, + "3": 500 + }, + "n_pooled_valid": 500, + "n_fit_errors": 0, + "min_valid_share": 1.0 + }, + { + "n_units": 300, + "reps": 500, + "per_horizon": { + "0": 0.95, + "1": 0.948, + "2": 0.946, + "3": 0.936 + }, + "mean_event_coverage": 0.945, + "pooled_att_coverage": 0.958, + "n_valid": { + "0": 500, + "1": 500, + "2": 500, + "3": 500 + }, + "n_pooled_valid": 500, + "n_fit_errors": 0, + "min_valid_share": 1.0 + } + ] +} \ No newline at end of file diff --git a/benchmarks/data/lpdid_test_panel.csv b/benchmarks/data/lpdid_test_panel.csv new file mode 100644 index 000000000..5827255ac --- /dev/null +++ b/benchmarks/data/lpdid_test_panel.csv @@ -0,0 +1,720 @@ +unit,time,treat,y,x +1,1,0,-0.737494152804906,0.294019556493925 +1,2,0,-1.45401392275265,-0.598844687188877 +1,3,0,0.743166315353444,1.14556035897094 +1,4,1,2.17379174056721,0.598931517876722 +1,5,1,1.11629494796406,0.853216548263565 +1,6,1,2.54915684708201,0.563026964359905 +1,7,1,2.81725203396961,0.479281608158206 +1,8,1,5.83581008102388,0.531293506020565 +1,9,1,5.13263653497541,-1.36160831678921 +1,10,1,3.14137523098187,-1.44382449655213 +1,11,1,3.5865855789882,-0.266844857953363 +1,12,1,6.17202200834693,1.51519276722892 +2,1,0,1.646001749588,0.210739579039542 +2,2,0,-0.933936349951549,-0.236565802121641 +2,3,0,-1.92831242499949,0.758552727506621 +2,4,1,2.30614944775482,0.700518102709616 +2,5,1,3.48524326239861,0.817886128251562 +2,6,1,2.97658900835797,-0.0311394346005515 +2,7,1,1.20268590865461,0.580882232704267 +2,8,1,4.33148010630701,-0.20180514731093 +2,9,1,6.32606812640205,-0.138705206869484 +2,10,1,3.52934491322057,0.599149463754405 +2,11,1,4.53552985929723,1.8001746854719 +2,12,1,4.92919284991123,-0.608106762873418 +3,1,0,0.0874535430131846,1.47708614464285 +3,2,0,0.46840946706061,1.11785752974202 +3,3,0,-2.63004947914542,-0.160723787604591 +3,4,1,1.66111545005929,-0.34512418076391 +3,5,1,3.02030038287294,0.945258834774553 +3,6,1,1.87774535096995,-0.818197787538497 +3,7,1,2.58304471042042,-1.0822868377944 +3,8,1,5.44387303908824,0.604401838989833 +3,9,1,7.8886093365895,0.936695558871581 +3,10,1,2.99090064911511,-0.464802307628442 +3,11,1,5.64731031668647,-0.671578507831782 +3,12,1,6.19883543130713,0.481976946108813 +4,1,0,-0.576472510381969,0.316856746783563 +4,2,0,3.0008479152642,0.670010315012874 +4,3,0,0.301515168123848,1.27802585444943 +4,4,1,2.25480963956575,-0.508633973286516 +4,5,1,3.9569418953648,0.97695801737316 +4,6,1,1.18877781128471,-1.21729409971739 +4,7,1,3.59208377980335,-0.638228775534505 +4,8,1,6.3792430335011,-0.645061518389099 +4,9,1,9.10984266624578,-0.329606260599108 +4,10,1,3.68678179266504,0.0577247424260506 +4,11,1,6.83099268297783,1.64640330460313 +4,12,1,4.39023973386204,-2.32003690347998 +5,1,0,-0.823475154934965,0.0570316367313002 +5,2,0,-0.00650422661640793,0.900898442159465 +5,3,0,-0.760541384984146,0.125724796933208 +5,4,1,2.43040437513097,0.836337372561178 +5,5,1,1.70995155773565,-0.677272565628402 +5,6,1,4.90021102807239,-0.979908968506341 +5,7,1,2.89928643915708,-1.29643182574397 +5,8,1,6.38616692608323,0.140614892597938 +5,9,1,8.40033721845525,-0.250741602837939 +5,10,1,4.48585537338913,0.80398786810283 +5,11,1,3.35113592135563,-0.497273499253265 +5,12,1,6.51357847220368,0.336411487689359 +6,1,0,0.104678625125971,-0.235406691588465 +6,2,0,0.324734334034314,0.16777575363705 +6,3,0,-1.25348276253401,0.432399220399878 +6,4,1,-1.33681316524322,-2.07526717426322 +6,5,1,1.00436444877079,-0.55053828691214 +6,6,1,0.677816837932778,-1.86228234564614 +6,7,1,1.85863531380113,-0.402069246608021 +6,8,1,4.85030573266757,1.82341730811489 +6,9,1,8.92737088567792,1.82994904984382 +6,10,1,3.921619183976,2.02623287400165 +6,11,1,3.02696134752673,-1.68710636859447 +6,12,1,6.87570526191578,-0.143529888435072 +7,1,0,-2.97918771671275,-0.975418459563763 +7,2,0,0.53267443161121,0.585834418570614 +7,3,0,-3.68116031465737,-1.39963613396742 +7,4,1,3.03106928103086,1.8040940567051 +7,5,1,1.95641452936406,-1.12067221926158 +7,6,1,1.16220423440467,-0.964930370911195 +7,7,1,1.89047797440029,0.346643442950692 +7,8,1,3.44469348155771,-0.414304552906162 +7,9,1,5.86979416555203,-0.0145154169533943 +7,10,1,1.92196991741356,0.599383868514409 +7,11,1,1.39885892352642,-0.802159634619603 +7,12,1,3.06493967286531,-1.92354411164103 +8,1,0,-0.47832060263948,1.82462169101574 +8,2,0,-1.50823921524675,-0.463087408357031 +8,3,0,-1.87616692146639,0.334542395752953 +8,4,1,0.359915556678139,-0.772864882492901 +8,5,1,3.30996505479371,2.16273846105896 +8,6,1,1.18890122745909,-0.73683293075437 +8,7,1,3.15371121272238,-0.336158251666822 +8,8,1,4.16864405631247,-0.41367584711935 +8,9,1,8.06078644691527,1.30163307550858 +8,10,1,2.81309076186481,1.20328054717193 +8,11,1,5.21865845134671,0.0667070309563468 +8,12,1,2.46023023645355,-3.35600410968268 +9,1,0,1.62268451942774,-0.10684551501785 +9,2,0,3.03744732994524,-0.398433987538361 +9,3,0,2.68687706828423,0.711228015881105 +9,4,1,3.82364065666178,-0.763059370239301 +9,5,1,3.94680848333765,-1.41190061560896 +9,6,1,5.22878678350309,-1.64429769441713 +9,7,1,5.14490066209261,0.399327113745293 +9,8,1,6.57639196813046,-3.24442955656723 +9,9,1,9.82043653541356,-1.53636694254794 +9,10,1,6.57470352058585,0.0345399451872066 +9,11,1,7.48034625490587,-1.16835051101502 +9,12,1,9.69583379955807,0.454607514866113 +10,1,0,2.26360423170703,1.34946795288848 +10,2,0,5.44564282039709,2.55251417617279 +10,3,0,0.179081004513205,-1.0163574700282 +10,4,1,2.71931418914351,0.0208916887576389 +10,5,1,4.09038990367592,0.655193159215659 +10,6,1,6.55566990747321,0.697740911944428 +10,7,1,6.19838657168851,0.650098834932255 +10,8,1,8.58434381850559,2.41841960958219 +10,9,1,10.3386890166897,-0.184816285806871 +10,10,1,8.28329836947991,-0.249629572274525 +10,11,1,8.51177285101523,1.35302124143609 +10,12,1,9.61389988928289,-1.3457799357315 +11,1,0,2.35434327562484,-0.597117206305814 +11,2,0,2.51365358173191,-0.629482347099027 +11,3,0,1.45045125349503,0.439797205122496 +11,4,1,3.25392580767134,-0.697529014915505 +11,5,1,4.93808550992731,0.206586368375649 +11,6,1,6.06206660391056,-0.255628778087342 +11,7,1,4.01328997136615,1.09631992171615 +11,8,1,7.7273082339794,0.947420766922036 +11,9,1,8.78353527709928,0.748077890119802 +11,10,1,6.14559190090774,-0.576096154322538 +11,11,1,7.65354827949617,-0.404718700335898 +11,12,1,8.92050045455535,0.0533667003477756 +12,1,0,-0.515334669295465,0.207322536291339 +12,2,0,-0.938779067913301,-0.391019334214117 +12,3,0,-1.22712757378703,-0.730644352585501 +12,4,1,0.84144508053606,0.188802081061093 +12,5,1,0.98690667549782,-1.55646240655012 +12,6,1,3.09395909690716,-1.14803742894224 +12,7,1,3.78812832184841,-0.529335578450324 +12,8,1,5.88233262363563,0.183674804793185 +12,9,1,10.1172777897868,1.27174627225927 +12,10,1,4.41299621316629,-0.272063950419708 +12,11,1,6.45157696896915,0.531165871047484 +12,12,1,8.9161337137643,0.999659167218455 +13,1,0,-1.65101102511835,1.70425554089418 +13,2,0,0.338313411263848,1.69871969577489 +13,3,0,-2.08001511137864,1.32026088559086 +13,4,0,-0.525263988459146,0.572769792613532 +13,5,0,-2.44395186587404,1.18681835624052 +13,6,0,-1.51130099592879,0.606773726409144 +13,7,0,-3.63192872140446,-0.576319512582858 +13,8,1,4.69825093496256,0.114128686249434 +13,9,1,3.54080519762415,0.782478366154311 +13,10,1,2.30882686549308,0.522820407932903 +13,11,1,3.22480072666722,1.69535369929462 +13,12,1,2.58893902688746,-0.541654317056403 +14,1,0,-0.197247330738393,-0.822975039591858 +14,2,0,0.959213275613963,0.241133688127109 +14,3,0,-0.867191511512427,-0.526864475541956 +14,4,0,-0.825430205424154,-1.09651124697817 +14,5,0,-2.70389180186096,-1.29269806254471 +14,6,0,-2.13552637182237,-0.526705295096638 +14,7,0,0.292895346921166,0.310960163354357 +14,8,1,4.79516464634862,-0.678127123592227 +14,9,1,6.07275434306472,-1.68828061910892 +14,10,1,3.03095127564554,-0.268286951211828 +14,11,1,2.78815689230502,-1.38846266226257 +14,12,1,5.16714561634034,-0.789170989459405 +15,1,0,-3.00278893270643,-0.788919076608423 +15,2,0,1.004304845345,1.72565692963065 +15,3,0,-0.672548212343948,0.231431158527356 +15,4,0,-0.635707403008096,-0.217297372116625 +15,5,0,-0.26239838875756,1.36142947120505 +15,6,0,-1.24757421041228,0.331594251419177 +15,7,0,0.471654904573003,0.382472089934457 +15,8,1,1.93975826180716,-2.04739717678117 +15,9,1,5.84576828154384,0.267329197841454 +15,10,1,0.948980200118755,-0.549392116608249 +15,11,1,5.26814289964413,1.64789906810577 +15,12,1,1.8903280756713,-1.92779886746096 +16,1,0,2.20707768806675,-0.466549048858403 +16,2,0,1.52341937643323,-0.0484310874078651 +16,3,0,1.04508675397628,0.022571428856729 +16,4,0,2.60402281337196,0.593741785343849 +16,5,0,0.619824652600638,-1.02403923476629 +16,6,0,3.17983187497731,-0.0349604961243257 +16,7,0,0.631473677838844,-0.3318641815444 +16,8,1,6.61358927982589,1.03995574944946 +16,9,1,9.01837636388115,-0.544529108314221 +16,10,1,4.65162486994155,-0.113111775852241 +16,11,1,3.83237853095296,-0.886473995935048 +16,12,1,5.50611532309779,-1.31753626887725 +17,1,0,0.846154482449593,0.955125936266205 +17,2,0,-0.97403124274424,-0.867956589796371 +17,3,0,-0.720078138125112,-1.26775296446946 +17,4,0,3.66142587251981,1.04873067297599 +17,5,0,4.30826871539927,-0.0103805875559447 +17,6,0,2.46240660941483,-1.25398437596292 +17,7,0,1.32675025458811,0.924237089183059 +17,8,1,8.11980858616472,0.228440276206889 +17,9,1,7.54404126770492,-0.425596110737993 +17,10,1,5.8107648180251,1.32799039638161 +17,11,1,3.92014028011573,-1.41295887557999 +17,12,1,6.58828802847078,-1.0007557402778 +18,1,0,2.52663222351609,0.707395560182477 +18,2,0,2.30541163099368,0.497535261978263 +18,3,0,2.23709457212529,-0.251643460322685 +18,4,0,1.99505207366139,-1.39495735362309 +18,5,0,2.78970785101564,1.53576229323095 +18,6,0,0.7720091187181,-2.02214454951655 +18,7,0,0.931588163683373,-0.0760332712788952 +18,8,1,8.50282885163737,0.138322610576226 +18,9,1,9.639240182544,-0.350910136946422 +18,10,1,5.11165633965826,-0.325523185829146 +18,11,1,6.76912278267685,0.549960854483448 +18,12,1,7.36280560795315,0.270478955381992 +19,1,0,4.08497538927326,2.51576275212678 +19,2,0,2.63808400495744,-0.122007947568827 +19,3,0,0.946918291219366,-0.718412507995877 +19,4,0,1.44876338238382,-0.617977673396205 +19,5,0,2.90213444537096,1.76066999629702 +19,6,0,-1.21628897604624,-0.553511145870931 +19,7,0,-2.31392096551496,-2.13378120608274 +19,8,1,6.13626477154039,1.50095439052937 +19,9,1,6.93865411286563,-1.22192491065537 +19,10,1,3.98233152713351,-0.945183530098006 +19,11,1,3.63285842657419,-0.319426862161486 +19,12,1,4.78387087199889,-0.0703091001332301 +20,1,0,-3.6639388096845,-1.3255271633008 +20,2,0,-1.84737062184273,-1.26906406909434 +20,3,0,-2.00724672707923,1.42329108101332 +20,4,0,-3.18713493522943,-1.78120480039952 +20,5,0,-3.20815151056298,-2.04036747701499 +20,6,0,-0.972907978661607,-0.498084815680567 +20,7,0,-1.30360987954144,1.79460203515538 +20,8,1,5.58945991361791,0.508097416926942 +20,9,1,8.50298895041954,0.842829163852509 +20,10,1,2.35785720743897,-0.0712697790990908 +20,11,1,3.29820009456482,-0.332057403379186 +20,12,1,7.05046852962555,1.6015154045222 +21,1,0,0.555819329795658,0.489661836638027 +21,2,0,-0.621197991536694,-1.13387126378313 +21,3,0,-1.07091561285194,1.56827623091664 +21,4,0,1.33445069777761,0.861193811926297 +21,5,0,-1.4888047731579,-0.755283228723595 +21,6,0,-1.7051964266187,0.465287063591231 +21,7,0,-1.73249998890916,-0.104860824539837 +21,8,1,4.3133145603143,-1.01724483250107 +21,9,1,6.15747520554308,-0.465698975781445 +21,10,1,5.83636694775812,1.25959525483508 +21,11,1,5.8804943652373,1.20423298219188 +21,12,1,6.47689397780052,1.78370032406465 +22,1,0,1.1772503663919,-0.133132565087424 +22,2,0,-0.423252821729625,-1.63465116748968 +22,3,0,1.27028074249367,1.95302212089924 +22,4,0,1.37563663029062,-0.985684377269498 +22,5,0,2.87447383946282,0.0876386650818677 +22,6,0,2.97441363764663,0.048832480281275 +22,7,0,0.262197120192222,-2.21292125998977 +22,8,1,7.4028009101668,-0.0167699080054312 +22,9,1,9.99833345771663,1.37790891285919 +22,10,1,6.18063870640222,0.0376236958050915 +22,11,1,5.52191521784527,-0.608828884398675 +22,12,1,6.51014777088796,0.244244438545098 +23,1,0,5.18743592455813,-0.551876533690783 +23,2,0,5.7621539446331,-0.134386626365369 +23,3,0,4.65022334145646,-0.734772169804066 +23,4,0,6.43348618968253,-0.488358648642302 +23,5,0,3.28782932023952,0.296907504885126 +23,6,0,5.53040580943089,1.01434302794517 +23,7,0,3.24429577857551,0.0900548684849309 +23,8,0,6.29080399480756,-0.376436759897836 +23,9,0,8.24162447982853,0.0225412005310621 +23,10,0,3.61945322757036,-0.0504369222470352 +23,11,0,3.31476405662982,-0.382029466759738 +23,12,0,5.25052787162893,0.641580616945033 +24,1,0,3.20483735684682,-1.17305951105671 +24,2,0,3.4338983487765,1.19046616835329 +24,3,0,1.5399242220226,1.31988849934392 +24,4,0,1.88074380065049,0.513683597661488 +24,5,0,4.37408304206491,1.09168547455265 +24,6,0,-0.381432046663667,-0.525725662389131 +24,7,0,0.787363451551131,0.236343131770071 +24,8,0,2.90736895080204,0.536046416679769 +24,9,0,1.62849406300933,-0.758307954420094 +24,10,0,1.84068636443256,-0.164896903990545 +24,11,0,1.99235011668306,0.203345963536017 +24,12,0,1.98756961918718,0.131147382531175 +25,1,0,4.54402184238126,0.433680855937419 +25,2,0,4.02126832197123,-0.0223413140107889 +25,3,0,1.73892999012968,0.167773259523592 +25,4,0,4.54094797894667,1.39822999479551 +25,5,0,4.46936515412985,0.369947648855894 +25,6,0,2.99454138430091,-0.728920574997309 +25,7,0,3.82247904177426,1.111675573869 +25,8,0,4.91975567394294,-0.227241658206269 +25,9,0,7.43141985493347,0.479514765057957 +25,10,0,3.99803734299001,-0.279226036482864 +25,11,0,3.43238233994097,0.0236322220969042 +25,12,0,3.95439797740879,0.86516714344901 +26,1,0,2.40787362970815,-0.540829014730045 +26,2,0,2.5018930037711,-0.15538354853331 +26,3,0,2.19482510878346,0.724767100300454 +26,4,0,1.34459968293508,-0.791582984379063 +26,5,0,2.25872420073315,0.0487380582080548 +26,6,0,0.117709593577951,-0.718571362760114 +26,7,0,1.28894465392044,0.770151772729173 +26,8,0,3.95689763856931,-0.902138982971895 +26,9,0,4.29779136563814,0.831318942049728 +26,10,0,0.547819418929679,-0.420157632208387 +26,11,0,2.68709119955431,0.147051183355696 +26,12,0,5.25357531902167,2.49300048675213 +27,1,0,1.38895724738085,-0.535366516276821 +27,2,0,0.947317668635352,-0.138212934474133 +27,3,0,-0.809163975319325,-0.465831892691777 +27,4,0,0.509255860047552,-0.189121431556514 +27,5,0,0.425551435827195,-1.32625944617564 +27,6,0,0.972314021720786,0.684057341793356 +27,7,0,1.99215560513198,1.29575806154921 +27,8,0,1.66914038457447,0.441664900713668 +27,9,0,4.56646216754484,-0.528124245358283 +27,10,0,2.61367506369179,1.45030891374817 +27,11,0,1.30251577533073,1.82006279724695 +27,12,0,0.00911831343448533,-1.2689769227096 +28,1,0,7.0573731413536,0.476394254010928 +28,2,0,6.53180569828399,-0.946393662084651 +28,3,0,7.10509979823759,1.21692488072893 +28,4,0,5.18378982012065,0.362530580178856 +28,5,0,5.71839470304202,0.942129392137432 +28,6,0,3.3280137739762,-0.589312041287479 +28,7,0,4.41931219425191,-0.111175250828361 +28,8,0,4.76213565837385,-0.513700277953363 +28,9,0,7.97924838803743,2.07904580780684 +28,10,0,2.29869668004013,-0.2655825700308 +28,11,0,5.37065384185892,0.825289584110518 +28,12,0,8.0343831762633,0.964462454503418 +29,1,0,-2.47323525610785,0.0224951161580127 +29,2,0,-0.93481926307591,0.990431064567552 +29,3,0,-4.19879336360184,-1.64594248664167 +29,4,0,-3.65122036987343,0.159429265361569 +29,5,0,-2.2495340314537,1.74536397866834 +29,6,0,-3.90500098936637,-0.607898525604142 +29,7,0,-5.65107987645139,-1.57979996383504 +29,8,0,-1.485917867964,0.125935244690911 +29,9,0,-1.49574541745671,-1.71005097704522 +29,10,0,-2.11902840390673,0.913194085998252 +29,11,0,-5.44566795577166,-0.282849711598633 +29,12,0,-2.87924426874221,0.0763290964825009 +30,1,0,-1.41659612189996,0.027065097570904 +30,2,0,-0.0431038987367287,-0.113046107259987 +30,3,0,-0.526873780850899,1.45185712042365 +30,4,0,0.588902329566442,0.558678398252781 +30,5,0,1.27852641508301,0.112264269144144 +30,6,0,-0.196511979019153,0.919591128198704 +30,7,0,-1.25304644236438,-1.42891246363454 +30,8,0,0.0651261166794273,0.330307771066111 +30,9,0,2.08414152542685,-0.720215680106555 +30,10,0,-0.0214842878065051,-0.730789443809481 +30,11,0,2.48245754414828,1.61539967224588 +30,12,0,1.10052293786728,0.578542683597726 +31,1,0,-0.540231406877571,-0.876547153983497 +31,2,0,-0.101148828837663,0.692882830849928 +31,3,0,-0.906265244139917,0.0359723709047931 +31,4,0,-0.504967399401651,0.451897465734737 +31,5,0,0.507725156058648,1.25778112084516 +31,6,0,0.0892688005456319,-0.182245250963657 +31,7,0,-0.904249912082254,0.54018482647941 +31,8,0,-0.120592312220213,-0.683007710662242 +31,9,0,4.43043433201751,1.30452873878101 +31,10,0,1.81013895938113,2.19527360768442 +31,11,0,-0.901718444052005,-0.0652648351919655 +31,12,0,-1.26426466580684,-1.18238250861937 +32,1,0,0.633210533425812,0.596668116453535 +32,2,0,1.031704322743,0.829625220898127 +32,3,0,-0.819154016693427,0.0897407946047844 +32,4,0,1.7306665461513,0.214567394687701 +32,5,0,2.11285603046468,0.588678510092288 +32,6,0,0.396313287272259,-0.976261114913195 +32,7,0,2.89640902929524,0.646729062190554 +32,8,0,3.0984446318974,0.997873705288833 +32,9,0,3.47747666554884,-1.17278310197572 +32,10,0,-2.35993829336561,-1.3897079873787 +32,11,0,0.943194827645753,0.959497498216224 +32,12,0,0.963570016025643,-0.263954223343817 +33,1,0,3.3280046575675,-1.5132855808917 +33,2,0,5.87046766964973,-0.106085654103736 +33,3,0,5.1951097204621,0.707526158648266 +33,4,0,4.44955421208068,-1.69490253833151 +33,5,0,4.60987942189189,-1.23070440706148 +33,6,0,5.59335298698762,1.07390299892405 +33,7,0,3.63419851119118,0.901980678631844 +33,8,0,5.56768893736595,-1.22531546534879 +33,9,0,8.07891022217767,0.675204900213813 +33,10,0,2.69411142224171,-0.591831183115575 +33,11,0,3.61045467449367,-0.346771251355396 +33,12,0,7.74598921443928,1.03959066520478 +34,1,0,0.355053174636498,-0.316728835996765 +34,2,0,-0.575700324508086,-0.989742235535738 +34,3,0,-3.09444765757537,-0.953400119670862 +34,4,0,-0.671275735785913,-0.256137247929617 +34,5,0,0.0668627804404278,-0.373840583348533 +34,6,0,-1.11884361416075,-1.27880950296398 +34,7,0,-0.869514831096626,0.780044183199852 +34,8,0,1.04344677206794,-1.36516789637856 +34,9,0,3.96834487061594,0.72702583270107 +34,10,0,1.75299595615307,1.06790506637333 +34,11,0,-0.655464079311993,-0.78776067042406 +34,12,0,-0.13466554515034,0.777147199172931 +35,1,0,-3.53836196208955,-0.741044952830363 +35,2,0,-2.81600375423442,-1.37311386251735 +35,3,0,-2.68490120652665,1.71630631762094 +35,4,0,-1.98501076115915,0.359611963067763 +35,5,0,-1.8310131588552,0.740473892669465 +35,6,0,-3.19343301444019,-0.346922332809428 +35,7,0,-2.7525915892214,0.202074439852311 +35,8,0,-0.077236440338722,0.206129290133392 +35,9,0,-1.9875822068138,-2.67149835200588 +35,10,0,-3.0892603480014,-1.95697596435864 +35,11,0,-2.55173839063754,1.13469613214802 +35,12,0,-2.10859901665625,-0.214919940130488 +36,1,0,2.40817327244319,-0.26913060492782 +36,2,0,2.88019405027981,-0.721352860332151 +36,3,0,2.27632229172175,0.332997672840165 +36,4,0,3.29097299439555,-0.967749811031793 +36,5,0,3.40695157113828,0.264722379800319 +36,6,0,2.03599246910829,-1.89662228609729 +36,7,0,2.82442478580712,0.906925038667774 +36,8,0,4.16559521605119,1.16158489065585 +36,9,0,2.62046698227912,-0.961291572842877 +36,10,0,2.04500151731009,0.530607351399573 +36,11,0,1.84153735850372,0.182165286248023 +36,12,0,3.76063167875658,1.38298492622236 +37,1,0,-1.82503693162603,-0.961068745111079 +37,2,0,-1.91373360764556,-0.00554726738275648 +37,3,0,-2.44279736686784,1.33609765940011 +37,4,0,-0.911337380114975,0.528384216666723 +37,5,0,-2.76084280850644,0.588247456740233 +37,6,0,-0.782305799187594,-2.38505054605771 +37,7,0,0.538489197845528,1.46701811212371 +37,8,0,0.976845035589607,-0.43986029149716 +37,9,0,3.6061428260028,-0.258194893643888 +37,10,0,0.411619112629319,0.609425393773574 +37,11,0,-0.279658830119308,1.28417251827404 +37,12,0,-2.23949526742441,0.388355688634707 +38,1,0,1.4917748522338,0.765427616543368 +38,2,0,2.9107046776598,0.576852197822651 +38,3,0,-0.0761724795611929,-0.0742504606686619 +38,4,0,1.93313149793402,0.744716528019745 +38,5,0,0.606798776043749,-0.906954665301575 +38,6,0,-0.162929296319538,0.729798268985491 +38,7,0,-1.78443757755426,-1.3960848148549 +38,8,0,2.23110431838977,-1.24818184538643 +38,9,0,6.20661716704552,-0.0421316807324575 +38,10,0,2.45698298076082,1.0978444965581 +38,11,0,1.58231632644836,-0.706132547579953 +38,12,0,2.45683711856108,-0.378590870967114 +39,1,0,0.342609071657654,1.07547920230413 +39,2,0,0.664103829608586,1.43730525474925 +39,3,0,-1.86783805455493,-1.51535174927019 +39,4,0,-1.15624418283872,0.97217858191987 +39,5,0,-3.14309284463586,-1.2042791802108 +39,6,0,0.229232276185235,0.578626551319532 +39,7,0,-2.21714115755208,1.11142145767804 +39,8,0,0.906478888291919,-0.50606111485816 +39,9,0,3.08739137936308,-0.216639620793998 +39,10,0,-0.66714082243591,0.0607215845641283 +39,11,0,-0.38916404285819,0.358488144876812 +39,12,0,2.55950686852103,2.25608884605937 +40,1,0,-1.40505995123184,0.468134459463131 +40,2,0,-2.27628737115565,-0.0223233741272032 +40,3,0,-3.18712931373241,-0.394789877129397 +40,4,0,-2.12548865731953,0.0859814930637629 +40,5,0,-2.51744899505794,-1.08678351280393 +40,6,0,-0.877678318484733,0.955463189501957 +40,7,0,-2.84884505639338,-0.137092651915791 +40,8,0,-0.850615986806339,0.70496451913069 +40,9,0,0.733013941979194,0.0727075968425066 +40,10,0,-3.85714646679222,-1.06756363903849 +40,11,0,-3.29536545142536,-1.48933761772441 +40,12,0,-0.727906532160943,-0.259959138299631 +41,1,0,1.69988777325805,-1.7161415526547 +41,2,0,2.63490934306202,-0.176290868407023 +41,3,0,-1.07524506760411,-1.11786095085709 +41,4,0,1.19543778759141,-0.86073277887357 +41,5,0,0.255169327931451,-0.770090997365728 +41,6,0,-0.502475504542747,-0.289355341152884 +41,7,0,0.133415451078003,-0.254970431527807 +41,8,0,4.10806665625316,0.681579168411048 +41,9,0,5.10063879572885,-0.351002668314993 +41,10,0,-0.407503131456817,-1.30159835527013 +41,11,0,2.48625620984401,1.19225640066255 +41,12,0,3.50624776284883,2.10285179266958 +42,1,0,1.75996455214612,0.682403463750592 +42,2,0,2.32982762334151,1.20543554813838 +42,3,0,-1.65570747936954,-0.92417812193615 +42,4,0,0.164228466673098,-0.0797698164233437 +42,5,0,-0.646434939254295,-0.925258829273138 +42,6,0,1.35013028545223,1.71754572183517 +42,7,0,-1.28995209293625,-1.18994336308786 +42,8,0,2.50516980762985,1.72108224548895 +42,9,0,5.31497768274573,2.13373948117311 +42,10,0,0.0294305424694456,0.477262538293024 +42,11,0,-0.0779095931942182,-1.58836526991872 +42,12,0,2.24591287455693,0.904883115949513 +43,1,0,3.77462383051829,0.472791600960528 +43,2,0,3.31990098735607,-1.26755639754727 +43,3,0,1.77901555228324,-0.137899055181795 +43,4,0,1.67238122775318,-1.48185206321787 +43,5,0,1.61942767433817,-0.692867939394719 +43,6,0,1.7217799782557,-1.39079275115377 +43,7,0,3.41376291317851,0.0579499939575364 +43,8,0,3.68132342015881,-1.63838004892575 +43,9,0,7.51624480070154,-0.0374546395224234 +43,10,0,3.98073958635326,1.81911976409549 +43,11,0,3.68120569504312,-0.34810441352505 +43,12,0,4.85451063593129,0.948675858582723 +44,1,0,1.47602528227136,0.417667012488565 +44,2,0,1.70688457090465,0.374653784634821 +44,3,0,2.1842888850979,0.997429924446094 +44,4,0,4.08629273074079,0.399906196038096 +44,5,0,2.2671329975573,0.987925003954454 +44,6,0,3.36889485617273,0.472038324760354 +44,7,0,2.95321942102128,-0.728940639204687 +44,8,0,4.23769352190405,0.436448015473034 +44,9,0,5.55423145214422,-1.15807165759875 +44,10,0,4.59475983150857,0.00427274639777722 +44,11,0,3.95440195050926,1.05840209151928 +44,12,0,3.94443713987822,-0.861992506198017 +45,1,0,3.78884986864662,0.767988411058982 +45,2,0,2.57906865758606,0.129910368746091 +45,3,0,-1.1768547393589,0.784581967337071 +45,4,0,1.47287973447031,-0.38163258664992 +45,5,0,2.10885948756905,-0.407863868314852 +45,6,0,3.38256404469586,-0.556523770763695 +45,7,0,1.19773404812737,0.499562694403946 +45,8,0,3.8589696725146,0.402350119087397 +45,9,0,6.26145467762668,-0.294770688705694 +45,10,0,2.32176698794748,-0.352824173057148 +45,11,0,-1.05463203826736,-0.656997753752131 +45,12,0,3.51152005618872,0.0823020143542185 +46,1,0,-2.91920510288547,0.783330735338723 +46,2,0,-3.46355265685011,0.368745341311848 +46,3,0,-4.29657213321202,-1.08920895206655 +46,4,0,0.403246085235683,0.763914317395046 +46,5,0,-1.68696168164031,-0.992338932179358 +46,6,0,-2.19445461147953,0.857236859389004 +46,7,0,-4.8889926626569,-0.574862590117261 +46,8,0,-0.575611163971309,-0.920657066460542 +46,9,0,-0.00677576042379358,-2.18272636567385 +46,10,0,-2.99494600065924,-0.787955949697475 +46,11,0,-2.30469418285391,-0.525786409882226 +46,12,0,-0.378778044009238,-1.04429774791425 +47,1,0,0.579905765486989,-0.198267385953029 +47,2,0,1.56507423512865,2.12978956232106 +47,3,0,0.769373744736294,0.412641630748787 +47,4,0,0.211719374852828,0.970358695853843 +47,5,0,1.62776507969353,0.237513570764907 +47,6,0,0.591238853105474,-0.814813548131481 +47,7,0,0.798374976809924,0.267653626666534 +47,8,0,3.14680921508134,1.1160641799409 +47,9,0,2.45489946120529,0.62267607494546 +47,10,0,-0.59978033659247,0.256798309828087 +47,11,0,-0.912999656187098,-0.190456768139995 +47,12,0,1.58535077078276,1.91979762838446 +48,1,0,-0.884790662756267,-1.20684006820841 +48,2,0,-0.38743165082115,0.983933774044088 +48,3,0,-3.4859926019302,-1.20811171260519 +48,4,0,-1.39715065787311,0.268001187406735 +48,5,0,-2.77005451654606,0.220873337272988 +48,6,0,-0.672896874753092,-0.341685993903828 +48,7,0,-1.88331562055561,0.178960030909073 +48,8,0,-1.03806789881738,0.299118858595076 +48,9,0,4.84418093294561,-0.276117140791951 +48,10,0,0.79121790979812,-1.80384916742251 +48,11,0,1.41030394503036,1.7057494038789 +48,12,0,-0.719594325629335,-0.677392816185464 +49,1,0,4.21233451706684,-0.494640270047562 +49,2,0,1.96232249389761,-1.78589831580548 +49,3,0,2.63576881272911,0.379165873424303 +49,4,0,4.27819779624658,-0.285625399269694 +49,5,0,3.63629173449065,-0.0303076635965998 +49,6,0,4.12080413771879,-0.392824035302597 +49,7,0,4.38122484195532,-1.36185230819394 +49,8,0,7.03682908225946,0.906072059771075 +49,9,0,8.29332758205941,-1.0910232270343 +49,10,0,3.21184716579145,-1.38684869640319 +49,11,0,2.84730496622726,-0.360097716537831 +49,12,0,3.6597548104162,-1.04712924982014 +50,1,0,-1.79754977485735,-0.404123669432513 +50,2,0,-1.03076422646484,0.0402003408299969 +50,3,0,-3.52773156004194,-0.88532373931668 +50,4,0,-2.02037030407731,-1.67746997455102 +50,5,0,-3.82321787604754,-1.64574568396371 +50,6,0,-0.687033617600416,1.02476538461847 +50,7,0,-0.848014392750634,0.852702042161652 +50,8,0,0.740872706964419,1.4469218235384 +50,9,0,2.03000448403283,-0.836297871793999 +50,10,0,1.36038800811183,2.4815759406954 +50,11,0,-0.799366783980821,-0.230293634932176 +50,12,0,1.30377016319762,-0.546493250953308 +51,1,0,0.204068547305652,2.06502081293007 +51,2,0,-1.53231430493351,0.100055579334191 +51,3,0,-1.29419755352033,-0.995849176274067 +51,4,0,-2.20541945971837,-0.433481812452031 +51,5,0,-0.342290285515447,-0.164351222325515 +51,6,0,-2.37513789513402,-0.82100097151649 +51,7,0,-0.639545653826987,0.627428167449279 +51,8,0,0.772066744827248,-1.16884654067355 +51,9,0,3.62156164018377,1.42405081529065 +51,10,0,-0.0769974233575679,1.18993228113453 +51,11,0,-1.5631278747557,0.371715667592506 +51,12,0,-1.26482946734126,-1.1371856995142 +52,1,0,-1.22067935033525,0.479717931768143 +52,2,0,-2.425361444624,0.325328736826614 +52,3,0,-3.67957566446919,-1.22491796206433 +52,4,0,-2.78537729236243,0.787747450146284 +52,5,0,-1.26872932434369,1.16863866933012 +52,6,0,-1.65991054401408,0.453370393188987 +52,7,0,-2.76784543351631,0.534921291821843 +52,8,0,-1.6831201948102,-0.758387398691424 +52,9,0,-2.59252763747817,-0.497751923081385 +52,10,0,-4.03447955968344,-0.137052421482375 +52,11,0,-2.59998270848158,-0.239828780220526 +52,12,0,-2.78612070742159,-1.38436184817713 +53,1,0,1.90250951412155,1.77294766575095 +53,2,0,2.90196314093133,-0.580505778033514 +53,3,0,-0.760872970845981,0.0387864235563768 +53,4,0,0.710640837873169,-1.30505910377026 +53,5,0,0.900264913835888,-1.94784345603364 +53,6,0,1.65553591956421,0.633728425127798 +53,7,0,2.03572400333897,1.3537283384036 +53,8,0,1.39819327354254,-0.183224610176147 +53,9,0,4.30370820489215,0.691583655272762 +53,10,0,-0.47868548558802,-0.145869543362797 +53,11,0,-0.550874634451685,-0.921683648928497 +53,12,0,-0.554250919936976,-1.66026089182096 +54,1,0,-2.99730186117659,1.14847318561438 +54,2,0,-2.55301026535854,0.942334968186566 +54,3,0,-4.66156133857617,0.233222776709222 +54,4,0,-4.88204073507546,-0.571178972802094 +54,5,0,-2.53146914272851,1.54015309779897 +54,6,0,-6.50141377783043,-0.507367468828441 +54,7,0,-6.33859897621678,-0.527152107694665 +54,8,0,-1.81151288324156,-0.874343694356267 +54,9,0,-0.0268271480869565,-0.39015866374651 +54,10,0,-2.88142879807442,0.692019424394767 +54,11,0,-4.03399033985865,-0.864288795477511 +54,12,0,-1.85846679864857,0.807191049657966 +55,1,0,2.76017557754694,1.0339399617585 +55,2,0,2.35747078099768,0.300368175694152 +55,3,0,-0.0635153891598355,1.11281323044467 +55,4,0,0.871831593494461,0.340911889689798 +55,5,0,0.361629393187786,-0.037784337862638 +55,6,0,0.211487805602075,-1.18525693272821 +55,7,0,0.720721390747349,-0.187602455256423 +55,8,0,3.09799938246643,-0.452513259817792 +55,9,0,5.50777010404874,1.00349453254821 +55,10,0,-0.839585474762627,-0.422064525167276 +55,11,0,1.53454902336008,1.30692566372475 +55,12,0,1.86198875365478,-0.212972178931566 +56,1,0,0.247894288218761,0.636952640920173 +56,2,0,-2.53312049728091,-0.32773404268582 +56,3,0,-5.06926903032977,-0.872347093859596 +56,4,0,-2.63255742000783,-0.223945964731547 +56,5,0,-2.29413094950663,-0.496582483155375 +56,6,0,-3.0331001633746,-0.51656653905686 +56,7,0,-0.813405150719412,1.51363820318432 +56,8,0,-1.52557639883955,0.36135682153405 +56,9,0,1.9303994384096,-0.779747875710288 +56,10,0,-3.25187202782774,-0.185590838639731 +56,11,0,-0.839021407138039,0.64091911686072 +56,12,0,0.445858472513891,0.109913586762345 +57,1,0,1.22158206817226,-0.167386892852535 +57,2,0,3.42905979347097,0.591468087360757 +57,3,0,-0.524615163385438,0.313888104082703 +57,4,0,-0.141008936510904,0.689862080632296 +57,5,0,3.38313071374649,2.139457802089 +57,6,0,2.40887701319438,1.10655786650145 +57,7,0,0.266083505511773,0.684951427749022 +57,8,0,-0.803765372814441,-1.20774825942069 +57,9,0,3.46308709648723,-0.437460128198174 +57,10,0,0.039150221524857,-0.723615616521702 +57,11,0,1.26883329938231,0.485156008623437 +57,12,0,1.13724023902607,-1.99810668787088 +58,1,0,1.0151274054329,0.451081233147967 +58,2,0,3.96749094667702,3.22357099777109 +58,3,0,-1.33155006453001,-0.748603382174008 +58,4,0,-0.689106363823263,0.0156023217283191 +58,5,0,-1.07075536044487,0.184579264971641 +58,6,0,-0.478159519121169,-0.386987401657154 +58,7,0,-0.867116075335817,-0.901042182426348 +58,8,0,2.28883099177806,0.674660564981887 +58,9,0,2.91510615278916,0.589747808040576 +58,10,0,-1.12473075741351,0.104251913070669 +58,11,0,-2.88652904484522,-1.62351787464913 +58,12,0,1.03076980024476,0.553877472634272 +59,1,0,4.15463490740598,1.52897153403284 +59,2,0,3.19334878989643,-0.999598917291693 +59,3,0,1.82749289946515,-0.392909406444298 +59,4,0,4.01785645723831,-1.0711655343625 +59,5,0,3.56040052677102,-0.0607296063419901 +59,6,0,4.93800429893885,0.602025890356261 +59,7,0,3.21668706702915,1.75218778470294 +59,8,0,1.15462225929893,-0.40802224850385 +59,9,0,7.23268317390191,0.0306359308030901 +59,10,0,2.3000251089639,-0.560522928574665 +59,11,0,3.59318036333106,0.0396951634163755 +59,12,0,3.7844201625516,-0.0342236654339522 +60,1,0,1.15254412082406,0.466559108899187 +60,2,0,3.20662254958859,1.509853539609 +60,3,0,0.423303189910291,-1.49237616594706 +60,4,0,2.05924858596634,1.09947032178921 +60,5,0,-0.0906012265639716,-1.65002096995697 +60,6,0,-0.157186329645417,-1.06155803593584 +60,8,0,0.412125119799941,0.564313376505113 +60,9,0,2.05979971568282,-2.02100178014148 +60,10,0,-0.946028786077311,-0.279184819077689 +60,11,0,2.69091956546499,1.32159479029395 +60,12,0,2.6560775199568,1.17220075989494 diff --git a/benchmarks/python/coverage_lpdid_ra.py b/benchmarks/python/coverage_lpdid_ra.py new file mode 100644 index 000000000..4273ff22b --- /dev/null +++ b/benchmarks/python/coverage_lpdid_ra.py @@ -0,0 +1,216 @@ +#!/usr/bin/env python3 +"""Monte-Carlo coverage study for the LP-DiD regression-adjustment (RA) SE. + +The canonical RA standard error is Stata ``teffects ra ... atet vce(cluster)`` +only - no R package computes it (``alexCardazzi`` uses direct covariate +inclusion, not RA), so the R-parity harness +(``tests/test_methodology_lpdid.py``) can anchor the RA *point* estimate but only +*pins* the library influence-function SE as a documented regression value. This +study validates that SE the way coverage is the ultimate test of a standard +error: simulate panels with a KNOWN dynamic treatment effect, fit the RA path +(``reweight=True`` + covariate), and check that the reported 95% CI covers the +true effect at ~95%. + +The RA IF variance is asymptotic (no finite-sample factor - the ``teffects`` +convention); the study sweeps the cluster count G and checks that empirical +coverage holds near nominal. In practice it is well-calibrated (~0.95) even at +modest G because the reported CI uses a ``t(G-1)`` reference distribution, which +widens to compensate at small G. Coverage at/near ~0.95 across G is the +validation. Mirrors +``benchmarks/python/coverage_sdid.py``: it lives under ``benchmarks/`` (never run +in gated CI - CI's isolated-install copies only ``tests/``) and writes an +artifact that underwrites ``docs/methodology/REGISTRY.md`` ## LPDiD Deviation 2. + +Usage:: + + python benchmarks/python/coverage_lpdid_ra.py # default sweep + python benchmarks/python/coverage_lpdid_ra.py --reps 200 # quick +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + +_REPO_ROOT = Path(__file__).resolve().parent.parent.parent +if str(_REPO_ROOT) not in sys.path: + sys.path.insert(0, str(_REPO_ROOT)) + +from diff_diff.lpdid import LPDiD # noqa: E402 + +PRE, POST = 2, 3 +# Homogeneous dynamic effect (so the true RA-ATT at horizon h is exactly TAU[h]). +TAU = {0: 1.0, 1: 1.5, 2: 2.0, 3: 2.5} +BETA_X = 0.7 +# True pooled-post effect: mean long-difference over [0, POST] == mean of TAU (the +# post window [0, POST] spans exactly TAU's horizons), so the pooled-row RA estimand is: +TRUE_POOLED = sum(TAU.values()) / len(TAU) + + +def simulate_panel(n_units: int, rng: np.random.Generator, n_periods: int = 10) -> pd.DataFrame: + """Staggered absorbing panel, ~40% never-treated, homogeneous effects + AR(1) noise.""" + cohorts = rng.choice([4, 6, 0], size=n_units, p=[0.3, 0.3, 0.4]) + unit_fe = rng.normal(0, 2.0, n_units) + time_fe = rng.normal(0, 1.0, n_periods + 1) + rows = [] + for u in range(n_units): + g = int(cohorts[u]) + eps_prev = 0.0 + for t in range(1, n_periods + 1): + eps = 0.3 * eps_prev + rng.normal(0, 1.0) + eps_prev = eps + treated = int(g > 0 and t >= g) + eff = TAU.get(t - g, 0.0) if treated else 0.0 + x = rng.normal(0, 1.0) + y = unit_fe[u] + time_fe[t] + eff + BETA_X * x + eps + rows.append((u + 1, t, treated, y, x)) + return pd.DataFrame(rows, columns=["unit", "time", "treat", "y", "x"]) + + +def run_coverage(n_units: int, reps: int, seed: int) -> dict: + rng = np.random.default_rng(seed) + hits = {h: 0 for h in TAU} + valid = {h: 0 for h in TAU} + pooled_hits = 0 + pooled_valid = 0 + n_fit_errors = 0 + for _ in range(reps): + panel = simulate_panel(n_units, rng) + try: + res = LPDiD(pre_window=PRE, post_window=POST, reweight=True, cluster="unit").fit( + panel, + outcome="y", + unit="unit", + time="time", + treatment="treat", + covariates=["x"], + ) + except Exception: + n_fit_errors += 1 + continue + es = res.event_study.set_index("horizon") + for h, tau in TAU.items(): + if h not in es.index: + continue + lo, hi = es.loc[h, "conf_low"], es.loc[h, "conf_high"] + if not (np.isfinite(lo) and np.isfinite(hi)): + continue + valid[h] += 1 + if lo <= tau <= hi: + hits[h] += 1 + # Headline pooled-row RA CI coverage of the true pooled effect - validates the + # SE that backs results.att/results.se (not just the event-study horizons). + plo, phi = res.conf_int + if np.isfinite(plo) and np.isfinite(phi): + pooled_valid += 1 + if plo <= TRUE_POOLED <= phi: + pooled_hits += 1 + coverage = {h: (hits[h] / valid[h] if valid[h] else float("nan")) for h in TAU} + mean_event_coverage = float(np.nanmean(list(coverage.values()))) + pooled_att_coverage = (pooled_hits / pooled_valid) if pooled_valid else float("nan") + min_valid_share = (min(min(valid.values()), pooled_valid) / reps) if reps else 0.0 + return { + "n_units": n_units, + "reps": reps, + "per_horizon": coverage, + "mean_event_coverage": mean_event_coverage, + "pooled_att_coverage": pooled_att_coverage, + "n_valid": valid, + "n_pooled_valid": pooled_valid, + "n_fit_errors": n_fit_errors, + "min_valid_share": min_valid_share, + } + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--reps", type=int, default=500, help="Monte-Carlo reps per G (default 500)") + ap.add_argument( + "--clusters", + type=int, + nargs="+", + default=[30, 100, 300], + help="cluster counts G to sweep (default 30 100 300)", + ) + ap.add_argument("--seed", type=int, default=20260629) + ap.add_argument( + "--out", + type=str, + default=str(_REPO_ROOT / "benchmarks" / "data" / "lpdid_ra_coverage.json"), + ) + args = ap.parse_args() + + results = [] + print(f"LP-DiD RA-SE Monte-Carlo coverage (nominal 0.95), reps={args.reps}") + print(f"{'G':>5} {'mean_ev':>8} {'pooled':>8} " + " ".join(f"h={h}" for h in TAU) + " errs") + for i, g in enumerate(args.clusters): + r = run_coverage(g, args.reps, args.seed + i) + results.append(r) + ph = " ".join(f"{r['per_horizon'][h]:.3f}" for h in TAU) + print( + f"{g:>5} {r['mean_event_coverage']:>8.3f} {r['pooled_att_coverage']:>8.3f} " + f"{ph} {r['n_fit_errors']}" + ) + + # Surface a broken or miscalibrated regeneration loudly rather than writing a + # misleading artifact (it underwrites REGISTRY ## LPDiD Deviation 2). Mass fit + # failures or NaN coverage are always hard errors. Out-of-band coverage is a hard + # error for substantive runs (reps >= 200, incl. the committed default 500), a warning + # for small noisy diagnostic runs. The aggregates (headline pooled-row RA coverage and + # the event-horizon mean) use a tight band; individual horizons are noisier, so a + # per-horizon miscalibration is caught with a slightly wider band (so a bad horizon + # cannot hide behind a good average). + SUBSTANTIVE_REPS = 200 + AGG_BAND = (0.90, 0.98) + PER_HORIZON_BAND = (0.88, 0.99) + + def _check(r: dict, name: str, cov: float, band: tuple) -> None: + if not np.isfinite(cov): + raise RuntimeError(f"G={r['n_units']}: {name} coverage is NaN") + if band[0] <= cov <= band[1]: + return + msg = f"G={r['n_units']} {name} coverage {cov:.3f} outside {list(band)}" + if args.reps >= SUBSTANTIVE_REPS: + raise RuntimeError( + msg + f" (reps={args.reps}) - the committed artifact underwrites a REGISTRY " + "claim, so a substantive run must not write an out-of-band result" + ) + print(f" WARNING: {msg} (reps={args.reps} < {SUBSTANTIVE_REPS}; noisy)") + + for r in results: + if r["min_valid_share"] < 0.5: + raise RuntimeError( + f"G={r['n_units']}: only {r['min_valid_share']:.0%} of reps produced a valid CI " + f"({r['n_fit_errors']} fit errors) - regeneration looks broken" + ) + _check(r, "pooled-att", r["pooled_att_coverage"], AGG_BAND) + _check(r, "mean-event", r["mean_event_coverage"], AGG_BAND) + for h, cov in r["per_horizon"].items(): + _check(r, f"h={h}", cov, PER_HORIZON_BAND) + + artifact = { + "study": "LPDiD RA influence-function SE coverage", + "nominal": 0.95, + "tau": {str(k): v for k, v in TAU.items()}, + "note": ( + "RA IF variance is asymptotic (teffects convention, no finite-sample " + "factor); empirical coverage of the true effect holds near nominal across G " + "for both the event-study horizons and the headline pooled-row RA CI - the " + "t(G-1) reference keeps it well-calibrated even at modest G. Underwrites " + "REGISTRY ## LPDiD Deviation 2." + ), + "sweep": results, + } + Path(args.out).parent.mkdir(parents=True, exist_ok=True) + with open(args.out, "w") as f: + json.dump(artifact, f, indent=2) + print(f"wrote {args.out}") + + +if __name__ == "__main__": + main() diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index b801bddb4..b25c950a3 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1844,11 +1844,12 @@ Weights are **always non-negative** (the central result). Via Frisch-Waugh-Lovel The paper specifies no standard-error formula (Section 1 defers to "standard, well-understood techniques"); the reference Stata `lpdid` uses `vce(cluster unit)`. The entries below document diff-diff's inference and scope choices. 1. **Note:** Standard errors are **cluster-robust at the unit level by default** - `cluster=None` auto-clusters at the unit identifier and the results record `cluster_name`/`n_clusters` - with a `t(G-1)` reference distribution (G = realized clusters in each horizon's clean-control sample). Matches Stata `lpdid` `vce(cluster unit)`; the paper prescribes no SE. -2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance currently carries **no finite-sample factor** (the ImputationDiD convention), while both paths share the `t(G-1)` reference. This RA-vs-default small-sample-scaling asymmetry is documented here and will be validated against the reference R packages and reconciled in the R-parity follow-up (PR-B2). +2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)` (the authors' `lpdid_regression_adjustment.do` `margins`/`kmatch` degrees-of-freedom comments prove `teffects` applies neither factor), while the default path matches `feols`/`reghdfe`. The RA *point* estimate is R-anchored to ~1e-13 (full-interaction `i.dtreat##(i.time c.x)` == `teffects` point; `tests/test_methodology_lpdid.py::test_ra_covariate_point`). The RA *standard error* itself has **no runnable R reference** (no R package computes the RA IF variance - `alexCardazzi` uses direct covariate inclusion, not RA; the canonical RA SE is Stata `teffects` only), so it is **pinned** as a documented regression value (`test_ra_covariate_se_regression_pin`) and its calibration is validated by the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). 3. **Note:** Direct covariate inclusion (`reweight=False` with covariates/absorb) emits a `UserWarning`: per online Appendix B.2.2 it preserves the non-negative LP-DiD weighting result only under linear and homogeneous covariate effects, so the regression-adjustment path (`reweight=True`) is preferred. 4. **Deviation from R:** Scope - this release implements the **absorbing-treatment main path only** (the estimator raises on non-absorbing input). Non-absorbing treatment (Section 4.2) and survey-design support are deferred to later PRs; the reference `lpdid` packages support non-absorbing treatment. 5. **Note:** LP-DiD's per-unit quantities (outcome lags `ylags`, first-difference lags `dylags`, integer-`pmd` premean baselines, treatment-entry detection) are **calendar** quantities (`t-1`, `t-k`), so the estimator requires integer-valued, globally consecutive `time` labels. A unit with an **interior time gap** is handled by reindexing that unit to its complete interior calendar grid `[min_t, max_t]`, computing the features on the grid, then **restricting back to the observed rows** - so a lag/first-difference spanning a gap is NaN and the observation fails closed (never the previous-*observed* row), and no synthetic gap row enters a regression. A gap-free panel skips this entirely and is bit-identical. **Entry = first OBSERVED treated period** (`min(t | D_it=1)`): an unobserved pre-onset gap cannot move a cohort earlier, the only well-defined convention when the true switch falls in an unobserved period. -6. **Note (pooled estimand):** The pooled pre/post ATT (the headline `results.att` is the pooled-post row) is the **unit-equal-weighted average of each unit-event-time's mean long difference** over the window - `mean_h(y_{i,t+h}) - baseline_{i,t}`, one observation per (unit, event-time), regressed on the treatment-switch indicator with event-time fixed effects on the **fixed-composition** sample (only units observing *every* pooled target, with clean controls required through `max(h)`). This equals the mean of the per-horizon event-study coefficients on a balanced panel; under cross-horizon composition changes it differs from the authors' horizon-**stacked** pooled regression (each (unit, event-time, horizon) long difference as a separate observation). It is therefore reported as a fixed-composition pooled ATT, and exact parity with the reference `lpdid` pooled specification is validated/reconciled in the R-parity follow-up (PR-B2). +6. **Note (pooled estimand):** The pooled pre/post ATT (the headline `results.att` is the pooled-post row) is the **unit-equal-weighted average of each unit-event-time's mean long difference** over the window - `mean_h(y_{i,t+h}) - baseline_{i,t}`, one observation per (unit, event-time), regressed on the treatment-switch indicator with event-time fixed effects on the **fixed-composition** sample (only units observing *every* pooled target, with clean controls required through `max(h)`). This equals the mean of the per-horizon event-study coefficients on a balanced panel. **PR-B2 validated it against the authors' runnable R reference**: the pooled estimand matches the authors' own R pooled recipe (`danielegirardi/lpdid`: a `slider` window-mean minus `y_{t-1}` on the clean-through-window-end sample) to ~1e-13 (`tests/test_methodology_lpdid.py::test_pooled`). A prior version of this note speculated the authors used a horizon-**stacked** pooled regression; the authors' R reference in fact uses this same fixed-composition mean-long-difference, so that speculation was incorrect. Unlike the event-study variants (where `alexCardazzi` is a cross-check gate), pooled is anchored to the authors' recipe **only**: `alexCardazzi`'s pooled uses a **laxer** clean-control window, so it differs and is recorded in the golden `meta` for transparency, not as a parity target. +7. **Deviation from R:** `no_composition` is intentionally more faithful to the paper's fixed-composition intent (Section 3.6) than the R packages: it fixes the realized sample across *all* post horizons (every post coefficient shares one sample, even on unbalanced panels) and excludes cohorts with `p_g > T-H`, whereas `alexCardazzi/lpdid` uses a looser per-horizon sample and a stricter `treat_date < T-H` cutoff. It therefore has **no exact R-package anchor** and is validated by the pure-Python tests in `tests/test_lpdid.py` (the R-parity golden omits it; `alexCardazzi`'s looser-semantics value is recorded in the golden `meta`). ### Implementation Checklist @@ -1862,7 +1863,7 @@ The paper specifies no standard-error formula (Section 1 defers to "standard, we - [x] `LPDiDResults` with `summary()` / `to_dict()` / cluster metadata (PR-B1) - [x] doc-deps.yaml mapping for `diff_diff/lpdid.py` + `lpdid_results.py`; llms.txt / llms-full.txt catalog entries (PR-B1, test-enforced) - [x] B1 pure-Python tests: analytical DGPs + cross-estimator equivalence (CS / BJS / DiD; Cengiz-stacked dropped, documented) + unbalanced / interior-gap / RA-overlap / pmd-missing edge cases (PR-B1) -- [ ] B2: self-generated R-parity (authors' `danielegirardi/lpdid` + `alexCardazzi/lpdid` cross-check) +- [x] B2: self-generated R-parity (authors' `danielegirardi/lpdid` recipes + `alexCardazzi/lpdid` cross-check; VW / reweight / pmd / direct / pooled / RA-point to ~1e-12; RA SE pinned + MC-coverage-validated; `no_composition` more paper-faithful than R, B1-tested) (PR-B2) - [ ] Non-absorbing extension (Section 4.2) - deferred to a later PR - [ ] Survey-design support - deferred to a later PR diff --git a/tests/test_methodology_lpdid.py b/tests/test_methodology_lpdid.py new file mode 100644 index 000000000..ac3b3bb3f --- /dev/null +++ b/tests/test_methodology_lpdid.py @@ -0,0 +1,228 @@ +"""R-parity tests for ``LPDiD`` (Dube, Girardi, Jorda & Taylor 2025), absorbing. + +Pins our ``LPDiD`` against goldens generated by +``benchmarks/R/generate_lpdid_golden.R`` from the method authors' own reference +recipes (the ``danielegirardi/lpdid`` ``fixest::feols`` event-study / reweight / +premean / pooled specifications) with an ``alexCardazzi/lpdid`` cross-check gate. +The paper specifies no SE formula, so validation is against the reference +software, not the paper. Point estimates and cluster-robust SEs match to ~1e-12 +on the reference platform; the tests assert ATT ``abs=1e-6`` / SE ``abs=1e-7`` for +cross-platform robustness (matching the ``didimputation`` parity precedent). + +Two paths get special handling, both documented in ``docs/methodology/REGISTRY.md`` +``## LPDiD`` Deviations: + +* **Regression adjustment (RA) standard error** (``ra_cov``): the canonical RA SE + is Stata ``teffects ra ... atet vce(cluster)`` only - no R package computes it + (``alexCardazzi`` uses direct covariate inclusion, not RA). We therefore anchor + the RA *point* estimate against R (full-interaction == ``teffects`` point) and + pin the library influence-function SE as a documented regression guard + (``RA_SE_PIN``); its calibration is validated out-of-band by the ungated + ``benchmarks/python/coverage_lpdid_ra.py`` Monte-Carlo coverage study. +* **no_composition** is intentionally not a parity variant: the library fixes the + realized sample across all post horizons (paper Section 3.6), more faithful than + any R package's per-horizon version, so it has no external anchor. It is + validated by the pure-Python tests in ``tests/test_lpdid.py``. + +Regenerate the goldens with:: + + Rscript benchmarks/R/generate_lpdid_golden.R +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from diff_diff.lpdid import LPDiD + +_DATA = Path(__file__).parent.parent / "benchmarks" / "data" +GOLDEN_PATH = _DATA / "lpdid_golden.json" +PANEL_PATH = _DATA / "lpdid_test_panel.csv" +_R_FIXTURE_AVAILABLE = GOLDEN_PATH.is_file() and PANEL_PATH.is_file() + +# Library regression-adjustment influence-function SE per horizon (reweight=True, +# covariates=["x"]). NOT R-parity: the canonical RA SE is Stata teffects only, with +# no runnable R analogue (see module docstring / REGISTRY ## LPDiD Deviation 2). +# Pinned regression guard. If the committed panel changes, refresh these from the +# library: fit LPDiD(pre_window=3, post_window=4, reweight=True, cluster="unit") on +# lpdid_test_panel.csv with covariates=["x"], then read event_study["se"] (RA_SE_PIN) +# and results.se / the pooled "pre" row se (POOLED_RA_SE_PIN). +RA_SE_PIN = { + -3: 0.4329626086028361, + -2: 0.4359293779913195, + 0: 0.45516589747218855, + 1: 0.4452290481291023, + 2: 0.4617939686659572, + 3: 0.39295354484873873, + 4: 0.39398602601018545, +} +# Headline (pooled) RA influence-function SE, same convention as RA_SE_PIN. +POOLED_RA_SE_PIN = {"post": 0.35917407532121975, "pre": 0.3785290141761256} + + +@pytest.fixture(scope="module") +def golden() -> dict: + if not _R_FIXTURE_AVAILABLE: + pytest.skip( + "R LPDiD parity fixture not present. Run " + "`Rscript benchmarks/R/generate_lpdid_golden.R` to regenerate " + "`benchmarks/data/lpdid_golden.json`." + ) + with GOLDEN_PATH.open("r") as f: + return json.loads(f.read()) + + +@pytest.fixture(scope="module") +def panel() -> pd.DataFrame: + if not _R_FIXTURE_AVAILABLE: + pytest.skip("R LPDiD parity fixture not present.") + return pd.read_csv(PANEL_PATH) + + +def _es_map(res) -> dict: + return { + int(r.horizon): (float(r.coefficient), float(r.se)) for r in res.event_study.itertuples() + } + + +def _num(v) -> float: + """Coerce a golden cell (number, JSON null, or a stray ``"NA"`` string) to float; + any non-numeric value becomes NaN so the missing-value branch fires instead of a + ``TypeError`` from calling ``np.isfinite`` on a non-number.""" + try: + return float(v) + except (TypeError, ValueError): + return float("nan") + + +def _assert_es(ours: dict, golden_es: dict, *, check_se: bool = True) -> None: + """Every golden horizon must be present and finite in ours, no silent skips.""" + assert golden_es, "golden variant has no horizons" + for hk, gv in golden_es.items(): + h = int(hk) + assert h in ours, f"missing horizon {h}" + o_att, o_se = ours[h] + g_att = _num(gv[0]) + if not np.isfinite(g_att): + # NaN-consistency: a missing golden horizon requires BOTH our att and se NaN + # (the library's all-inference-fields-NaN-together contract), not just att. + assert not np.isfinite(o_att) and not np.isfinite( + o_se + ), f"h={h}: golden att NA but ours not fully NaN (att={o_att}, se={o_se})" + continue + assert np.isfinite(o_att), f"h={h}: non-finite att" + assert o_att == pytest.approx(g_att, abs=1e-6), f"att h={h}" + if check_se: + g_se = _num(gv[1]) + if not np.isfinite(g_se): + assert not np.isfinite(o_se), f"h={h}: golden se NA but ours finite" + else: + assert np.isfinite(o_se), f"h={h}: non-finite se" + assert o_se == pytest.approx(g_se, abs=1e-7), f"se h={h}" + + +def _window(golden: dict): + return golden["meta"]["pre_window"], golden["meta"]["post_window"] + + +def test_assert_es_tolerates_missing_golden_horizons() -> None: + """A missing golden cell (JSON null, stray 'NA' string, or NaN) is handled as an + empty-result horizon - consistent when ours is also NaN, a failure when ours is + finite - never a TypeError.""" + for na in (None, "NA", float("nan")): + _assert_es( + {0: (1.0, 0.1), 1: (float("nan"), float("nan"))}, + {"0": [1.0, 0.1], "1": [na, na]}, + ) + # ours finite where golden is NA -> inconsistent, must fail + with pytest.raises(AssertionError): + _assert_es({1: (5.0, 0.1)}, {"1": [na, na]}) + # NaN-consistency: att NaN but se still finite must also fail (both must be NaN) + with pytest.raises(AssertionError): + _assert_es({1: (float("nan"), 0.1)}, {"1": [na, na]}) + + +class TestLPDiDParityR: + """Pin ``LPDiD`` against the authors' reference recipes + alexCardazzi cross-check.""" + + def test_vw_event_study(self, golden: dict, panel: pd.DataFrame) -> None: + pre, post = _window(golden) + res = LPDiD(pre_window=pre, post_window=post, reweight=False, cluster="unit").fit( + panel, outcome="y", unit="unit", time="time", treatment="treat" + ) + _assert_es(_es_map(res), golden["vw_es"]) + + def test_ew_reweight_event_study(self, golden: dict, panel: pd.DataFrame) -> None: + pre, post = _window(golden) + res = LPDiD(pre_window=pre, post_window=post, reweight=True, cluster="unit").fit( + panel, outcome="y", unit="unit", time="time", treatment="treat" + ) + _assert_es(_es_map(res), golden["ew_es"]) + + def test_pmd_event_study(self, golden: dict, panel: pd.DataFrame) -> None: + pre, post = _window(golden) + res = LPDiD(pre_window=pre, post_window=post, reweight=False, pmd=2, cluster="unit").fit( + panel, outcome="y", unit="unit", time="time", treatment="treat" + ) + _assert_es(_es_map(res), golden["pmd_es"]) + + def test_direct_covariate_event_study(self, golden: dict, panel: pd.DataFrame) -> None: + pre, post = _window(golden) + # Direct inclusion (reweight=False + covariates) must emit the Appendix B.2.2 + # homogeneity warning (REGISTRY ## LPDiD Deviation 3). + with pytest.warns(UserWarning, match="direct inclusion"): + res = LPDiD(pre_window=pre, post_window=post, reweight=False, cluster="unit").fit( + panel, covariates=["x"], outcome="y", unit="unit", time="time", treatment="treat" + ) + _assert_es(_es_map(res), golden["direct_es"]) + + def test_pooled(self, golden: dict, panel: pd.DataFrame) -> None: + pre, post = _window(golden) + res = LPDiD(pre_window=pre, post_window=post, reweight=False, cluster="unit").fit( + panel, outcome="y", unit="unit", time="time", treatment="treat" + ) + g_post = golden["pooled"]["post"] + assert np.isfinite(res.att) and np.isfinite(res.se) + assert res.att == pytest.approx(g_post[0], abs=1e-6) + assert res.se == pytest.approx(g_post[1], abs=1e-7) + # pre pooled row + pre_row = res.pooled.loc[res.pooled["window"] == "pre"] + assert len(pre_row) == 1 + g_pre = golden["pooled"]["pre"] + assert float(pre_row["coefficient"].iloc[0]) == pytest.approx(g_pre[0], abs=1e-6) + assert float(pre_row["se"].iloc[0]) == pytest.approx(g_pre[1], abs=1e-7) + + def test_ra_covariate_point(self, golden: dict, panel: pd.DataFrame) -> None: + """RA point estimate == authors' teffects/full-interaction point (R-anchored).""" + pre, post = _window(golden) + res = LPDiD(pre_window=pre, post_window=post, reweight=True, cluster="unit").fit( + panel, covariates=["x"], outcome="y", unit="unit", time="time", treatment="treat" + ) + _assert_es(_es_map(res), golden["ra_cov"], check_se=False) + + def test_ra_covariate_se_regression_pin(self, panel: pd.DataFrame) -> None: + """RA influence-function SE pinned (no runnable R reference; see RA_SE_PIN). + + Calibration validated separately by benchmarks/python/coverage_lpdid_ra.py. + """ + res = LPDiD(pre_window=3, post_window=4, reweight=True, cluster="unit").fit( + panel, covariates=["x"], outcome="y", unit="unit", time="time", treatment="treat" + ) + es = {int(r.horizon): float(r.se) for r in res.event_study.itertuples()} + for h, pinned in RA_SE_PIN.items(): + assert h in es, f"missing horizon {h}" + assert np.isfinite(es[h]) and es[h] > 0, f"h={h}: RA SE not finite/positive" + assert es[h] == pytest.approx(pinned, abs=1e-6), f"RA SE pin drift at h={h}" + # Headline (pooled) RA SE is the same convention - pin it too so pooled RA + # inference drift is also caught (not just the event-study horizons). + assert np.isfinite(res.se) and res.se > 0 + assert res.se == pytest.approx( + POOLED_RA_SE_PIN["post"], abs=1e-6 + ), "pooled-post RA SE drift" + pre_se = float(res.pooled.loc[res.pooled["window"] == "pre", "se"].iloc[0]) + assert pre_se == pytest.approx(POOLED_RA_SE_PIN["pre"], abs=1e-6), "pooled-pre RA SE drift"