feat(rdd): covariate-adjusted RD - covariates= on fit() with rdrobust 4.0.0 parity (CCFT 2019)

Adds the CCFT 2019 additive common-coefficient covariate adjustment to
RegressionDiscontinuity (sharp AND fuzzy), parity-targeting R rdrobust
4.0.0 end-to-end. The estimand is unchanged (precision only, unlike the
DiD estimators' conditional-parallel-trends covariates role); bandwidths
are covariate-aware (Z threads into every pilot of all three selector
chains with a per-pilot partialled gamma); collinear covariates are
dropped with a warning naming them under covs_drop=True (R's exact
dqrdc2 rank/pivot semantics incl. the name-length column sort, ported
directly; covs_drop=False is a deterministic strict error).

Degenerate covariate adjustment is GUARDED, not reproduced (documented
Deviation from R): R's ginv(tol=1e-20) inverts a float-noise singular
value on constant covariates / full dummy sets, silently returning
platform-dependent estimates; diff-diff excludes per-column degeneracies
(a constant covariate reproduces the fit without it bit-for-bit),
applies a scale-invariant stabilized cut for rank-deficient sets (a full
dummy set reproduces the drop-one-category fit - span invariance), and
warns naming the columns. Well-posed systems match R at machine
precision.

Surfaces: covariates= (fit-time, library-wide name), covs_drop=
constructor knob, name-keyed covariate_coefficients /
first_stage_covariate_coefficients + covariates/covariates_dropped/
covs_drop result echoes, covariate-adjusted summary() banner. Goldens
23 -> 32 configs (9 covariate configs incl. msetwo/cercomb2 chain pins,
collinear-drop, ties, fuzzy, fuzzy-sharpbw; all 23 pre-existing configs
reproduced exactly on regeneration; worst covariate-config disagreement
6.2e-12). R-free anchors: the CCFT 2019 partial-out identity
tau_adj = tau_unadj - gamma' tau_Z exact at common manual bandwidths
(both conventional and bias-corrected rows), span/order invariance,
CI shrinkage. Docs: REGISTRY covariate subsection + deviation notes,
CCFT-2019/2017 review checkbox flips, api rst, llms guides, README,
CHANGELOG, references, choosing-estimator.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QGca52n6H8oDDXALjjrsp4

Author: igerber

CHANGELOG.md

@@ -8,8 +8,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-- **`RegressionDiscontinuity` - sharp AND fuzzy regression discontinuity estimation
-  with robust bias-corrected inference (alias `RDD`).** Local-polynomial RD per
+- **`RegressionDiscontinuity` - sharp, fuzzy, AND covariate-adjusted regression
+  discontinuity estimation with robust bias-corrected inference (alias `RDD`).** Local-polynomial RD per
   Calonico, Cattaneo & Titiunik (2014), parity-targeting R `rdrobust` 4.0.0
   end-to-end: all 10 data-driven bandwidth selectors (`mserd` default,
   `msetwo`/`msesum`/comb and the CER-optimal variants),
@@ -29,6 +29,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   robust CI contains zero (documented deviation - R is silent; CCT 2014 Theorem 3
   "guard and warn", Feir-Lemieux-Marmer weak-IV inference a documented seam); R's
   exact no-variation-no-jump identification error is raised on both entry points.
+  **Covariate adjustment** via `fit(..., covariates=[...])` (R's `covs=`; CCFT 2019):
+  additive common-coefficient adjustment with a pooled-across-sides gamma that
+  leaves the estimand UNCHANGED (precision only - explicitly unlike the DiD
+  estimators' conditional-parallel-trends `covariates` role; the operative,
+  testable requirement is covariate balance at the cutoff, and the placebo recipe
+  - fit each covariate as the outcome - is documented); bandwidths are
+  covariate-aware (Z stacked into every pilot with a per-pilot gamma); collinear
+  columns are dropped with a warning NAMING them under `covs_drop=True` (R's
+  default and its exact dqrdc2 rank/pivot semantics incl. the name-length column
+  sort, ported directly; `covs_drop=False` = deterministic strict error); fitted
+  gammas are exposed name-keyed (`covariate_coefficients`, fuzzy
+  `first_stage_covariate_coefficients`) and `covariates`/`covariates_dropped`/
+  `covs_drop` echo on the results. Degenerate adjustments are GUARDED, not
+  reproduced (documented deviation): R's `ginv(tol=1e-20)` inverts a float-noise
+  singular value on constant covariates / full dummy sets, silently returning
+  platform-dependent estimates; diff-diff excludes degenerate columns, applies a
+  scale-invariant stabilized cut for rank-deficient sets (a full dummy set
+  reproduces the drop-one-category fit exactly), and warns naming the columns.
   **Canonical binding:** `att`/`se`/`t_stat`/`p_value`/`conf_int` are ONE
   coherent row - the robust bias-corrected row (`att = tau_bc`, CI centered on it,
   `t_stat == att/se`), preserving the library-wide field identities; the new
@@ -37,20 +55,26 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   row, and `summary()` prints the familiar three-row table. Estimation-path port
   (`rdrobust_fit`: Q_q bias-correction score matrix, conventional/robust NN
   sandwiches, fuzzy ratio/first-stage variances) validated against a new estimates
-  golden (`benchmarks/data/rdrobust_estimates_golden.json`, 23 configurations incl.
-  the Senate anchors and 7 fuzzy configs - default/sharpbw/manual-h/epa/msetwo/
-  one-sided-compliance/ties) at rtol=1e-9 in `tests/test_rdd_parity.py` (+ port-level
-  linearized-bias pins in `tests/test_rdrobust_port.py`); R-free methodology
-  anchors (CCT 2014 Remark 7 bias-corrected == local-quadratic equivalence at rel
-  1e-10 across all kernels, perfect-compliance == sharp reproduction, bandwidth
-  auto-switch locks, invariances, joint-NaN degenerate contracts) in
+  golden (`benchmarks/data/rdrobust_estimates_golden.json`, 32 configurations incl.
+  the Senate anchors, 7 fuzzy configs - default/sharpbw/manual-h/epa/msetwo/
+  one-sided-compliance/ties - and 9 covariate configs - default/manual-h/msetwo/
+  cercomb2/epa/collinear-drop/ties/fuzzy/fuzzy-sharpbw with `coef_covs` gamma
+  pins) at rtol=1e-9 in `tests/test_rdd_parity.py` (+ port-level linearized-bias,
+  gamma-matrix, and dqrdc2 rank/pivot pins in `tests/test_rdrobust_port.py`);
+  R-free methodology anchors (CCT 2014 Remark 7 bias-corrected == local-quadratic
+  equivalence at rel 1e-10 across all kernels, perfect-compliance == sharp
+  reproduction, bandwidth auto-switch locks, invariances, joint-NaN degenerate
+  contracts, and the CCFT 2019 partial-out identity - exact at common manual
+  bandwidths - plus covariate span-/order-invariance and CI-shrinkage anchors) in
   `tests/test_rdd_methodology.py`; API/validation suite in `tests/test_rdd.py`.
   Deviations from R (each labeled in the REGISTRY section): warn-instead-of-silent
   NaN drops, warn-and-ignore `b`-without-`h`, the weak-first-stage warning,
-  warn-and-ignore `sharpbw` on sharp fits, fail-closed targeted errors on
-  degenerate designs, and the canonical-binding note above. Covariates (CCFT 2019 -
-  review on file), cluster-robust variance, weights, kink estimands, weak-IV-robust
-  fuzzy inference, and rdplot/density diagnostics are documented follow-ups.
+  warn-and-ignore `sharpbw` on sharp fits (and `covs_drop=False` without
+  covariates), fail-closed targeted errors on degenerate designs, the guarded
+  degenerate covariate adjustment, and the canonical-binding note above.
+  Cluster-robust variance, weights, kink estimands, weak-IV-robust fuzzy
+  inference, a packaged covariate-balance helper, and rdplot/density diagnostics
+  are documented follow-ups.
 - **Internal: mypy enforced at zero errors.** Triaged the 184 pre-existing
   `mypy diff_diff` errors to an enforceable zero and added a blocking Mypy job to
   the Lint CI workflow (pinned `mypy==2.1.0` + pinned numpy/pandas/scipy for stub
@@ -134,9 +158,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - `diff_diff/guides/llms-autonomous.txt` no longer lists regression discontinuity as
-  out of scope: sharp AND fuzzy RD route to `RegressionDiscontinuity`; only kink
-  designs and the covariate-adjusted / cluster-robust RD variants are referred to
-  external tooling.
+  out of scope: sharp, fuzzy, AND covariate-adjusted RD route to
+  `RegressionDiscontinuity`; only kink designs and the cluster-robust RD variant
+  are referred to external tooling.
 - **Internal: repo-wide lint normalization + pinned tooling.** black/ruff/mypy are now
   pinned exactly in the `dev` extra (`black==26.3.1`, `ruff==0.15.13`, `mypy==2.1.0`;
   the tools require Python >= 3.10 — the library floor stays 3.9); full `black` +

README.md

@@ -112,7 +112,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`.
 - [TripleDifference](https://diff-diff.readthedocs.io/en/stable/api/triple_diff.html) - triple difference (DDD) estimator for designs requiring two criteria for treatment eligibility
 - [ContinuousDiD](https://diff-diff.readthedocs.io/en/stable/api/continuous_did.html) - Callaway, Goodman-Bacon & Sant'Anna (2024) continuous treatment DiD with dose-response curves
 - [HeterogeneousAdoptionDiD](https://diff-diff.readthedocs.io/en/stable/api/had.html) - de Chaisemartin, Ciccia, D'Haultfœuille & Knau (2026) for designs where **no unit remains untreated**; local-linear estimator at the dose support boundary returning Weighted Average Slope (WAS) on Design 1' (`d̲ = 0` / QUG) or `WAS_{d̲}` on Design 1 (`d̲ > 0`, continuous-near-d̲ or mass-point), with a multi-period event-study extension (last-treatment cohort, pointwise CIs). **Panel-only** in this release - repeated cross-sections rejected by the validator. Alias `HAD`.
-- [RegressionDiscontinuity](https://diff-diff.readthedocs.io/en/stable/api/regression_discontinuity.html) - Calonico, Cattaneo & Titiunik (2014) sharp AND fuzzy regression discontinuity with robust bias-corrected inference and rdrobust-parity bandwidth selection (all 10 selectors, mass-point handling; fuzzy via `treatment_col=` with a first-stage block and weak-identification warning). Canonical `att` is the bias-corrected estimate with a coherent robust CI (rdrobust's printed headline is `att_conventional`). Alias `RDD`.
+- [RegressionDiscontinuity](https://diff-diff.readthedocs.io/en/stable/api/regression_discontinuity.html) - Calonico, Cattaneo & Titiunik (2014) sharp, fuzzy, AND covariate-adjusted regression discontinuity with robust bias-corrected inference and rdrobust-parity bandwidth selection (all 10 selectors, mass-point handling; fuzzy via `treatment_col=` with a first-stage block and weak-identification warning; covariates via `covariates=` - CCFT 2019, same estimand, covariate-aware bandwidths). Canonical `att` is the bias-corrected estimate with a coherent robust CI (rdrobust's printed headline is `att_conventional`). Alias `RDD`.
 - [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html) - Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments; optional covariate balancing (Ustyuzhanin 2026)
 - [EfficientDiD](https://diff-diff.readthedocs.io/en/stable/api/efficient_did.html) - Chen, Sant'Anna & Xie (2025) efficient DiD with optimal weighting for tighter SEs
 - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html) - Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment

benchmarks/R/generate_rdrobust_estimates_golden.R

@@ -1,8 +1,11 @@
-# Golden-value generator for the diff-diff RD ESTIMATION port - sharp AND
-# fuzzy (diff_diff/_rdrobust_port.py::rdrobust_fit and the public
-# RegressionDiscontinuity estimator; 23 configs across four synthetic DGPs
-# + the Senate data, incl. 7 fuzzy configs with full first-stage
-# tau_T/se_T/z_T/pv_T/ci_T blocks and per-side take-up coefficients).
+# Golden-value generator for the diff-diff RD ESTIMATION port - sharp,
+# fuzzy, AND covariate-adjusted (diff_diff/_rdrobust_port.py::rdrobust_fit
+# and the public RegressionDiscontinuity estimator; 32 configs across five
+# synthetic DGPs + the Senate data, incl. 7 fuzzy configs with full
+# first-stage tau_T/se_T/z_T/pv_T/ci_T blocks and 9 covariate configs with
+# coef_covs (gamma) pins; covariate names deliberately differ in length
+# and are passed UNSORTED so every covariate config also pins rdrobust's
+# order(nchar(colnames)) column sort, rdrobust.R:131).
 #
 # Deliberately a SEPARATE file/JSON from generate_rdrobust_golden.R so the
 # bandwidth fixtures reviewed in the machinery PR are never regenerated.
@@ -27,19 +30,24 @@ TARBALL_SHA256 <- "78f0d6b4bdec4091cc8f42f6f1598704747f95926446d3aaee381ea1d613a
 run_estimate <- function(y, x, c = 0, masspoints = "adjust", kernel = "tri",
                          p = 1, q = 2, h = NULL, b = NULL, rho = NULL,
                          level = 95, bwselect = "mserd",
-                         fuzzy = NULL, sharpbw = FALSE) {
+                         fuzzy = NULL, sharpbw = FALSE,
+                         covs = NULL, covs_drop = TRUE) {
   args <- list(y = y, x = x, c = c, masspoints = masspoints, kernel = kernel,
                p = p, q = q, level = level, bwselect = bwselect,
-               sharpbw = sharpbw)
+               sharpbw = sharpbw, covs_drop = covs_drop)
   if (!is.null(h)) args$h <- h
   if (!is.null(b)) args$b <- b
   if (!is.null(rho)) args$rho <- rho
   if (!is.null(fuzzy)) args$fuzzy <- fuzzy
+  if (!is.null(covs)) args$covs <- covs
   r <- suppressWarnings(do.call(rdrobust, args))
   out <- list(
     c = c, masspoints = masspoints, kernel = kernel, p = p, q = q,
     bwselect = bwselect,
     fuzzy_in = !is.null(fuzzy), sharpbw = sharpbw,
+    covs_in = !is.null(covs),
+    covs_names = if (is.null(covs)) NA else colnames(covs),
+    covs_drop = covs_drop,
     h_in = if (is.null(h)) NA else h,
     b_in = if (is.null(b)) NA else b,
     rho_in = if (is.null(rho)) NA else rho,
@@ -66,6 +74,12 @@ run_estimate <- function(y, x, c = 0, masspoints = "adjust", kernel = "tri",
     out$beta_t_p_l <- unname(as.vector(r$beta_T_p_l))
     out$beta_t_p_r <- unname(as.vector(r$beta_T_p_r))
   }
+  if (!is.null(covs)) {
+    # Common projection coefficients gamma (dZ_kept x 1 sharp, x 2 fuzzy)
+    # over the covariates KEPT after covs_drop, in R's nchar-sorted column
+    # order; the row count pins WHICH columns survived the drop.
+    out$coef_covs <- unname(as.matrix(r$coef_covs))
+  }
   out
 }
 
@@ -75,14 +89,15 @@ golden$metadata <- list(
   rdrobust_version = as.character(packageVersion("rdrobust")),
   rdrobust_tarball_sha256 = TARBALL_SHA256,
   seeds = list(dgp_lee_smooth = 42L, dgp_ties_moderate = 123L,
-               dgp_asymmetric_scaled = 777L, dgp_fuzzy = 314L),
+               dgp_asymmetric_scaled = 777L, dgp_fuzzy = 314L,
+               dgp_covs = 2718L),
   generator = "benchmarks/R/generate_rdrobust_estimates_golden.R",
   algorithm = paste(
-    "rdrobust() sharp AND fuzzy estimation blocks (three-row coef/se/z/pv/ci,",
-    "counts, per-side beta_p; fuzzy configs add the first-stage",
-    "tau_T/se_T/z_T/pv_T/ci_T rows and per-side beta_T_p) for the vce='nn'",
-    "no-covariate path, complementing the bandwidth fixtures in",
-    "rdrobust_golden.json."
+    "rdrobust() sharp, fuzzy, AND covariate-adjusted estimation blocks",
+    "(three-row coef/se/z/pv/ci, counts, per-side beta_p; fuzzy configs add",
+    "the first-stage tau_T/se_T/z_T/pv_T/ci_T rows and per-side beta_T_p;",
+    "covariate configs add the coef_covs gamma matrix) for the vce='nn'",
+    "path, complementing the bandwidth fixtures in rdrobust_golden.json."
   ),
   r_version = R.version.string
 )
@@ -162,6 +177,42 @@ golden$dgp_fuzzy <- list(
   )
 )
 
+# Covariate DGP: two informative covariates with NAME LENGTHS that differ
+# and are passed UNSORTED (c("zlong", "zb")) so R's order(nchar) column
+# sort (rdrobust.R:131) is exercised by every config; zdup is an EXACT
+# linear combination for the covs_drop config. covs_ties reuses the
+# 2dp-rounded running variable (masspoints machinery x covariates).
+set.seed(2718)
+n5 <- 1200
+x5 <- 2 * rbeta(n5, 2, 4) - 1
+zlong <- 0.5 * x5 + rnorm(n5, sd = 0.8)
+zb <- rbinom(n5, 1, 0.4)
+y5 <- 0.4 * x5 + 0.9 * (x5 >= 0) + 0.7 * zlong + 0.3 * zb + rnorm(n5, sd = 0.3)
+t5 <- rbinom(n5, 1, ifelse(x5 >= 0, 0.75, 0.2))
+zdup <- 1.5 * zlong - 0.5 * zb
+x5_ties <- round(x5, 2)
+covs2 <- cbind(zlong = zlong, zb = zb)
+covs3 <- cbind(zlong = zlong, zb = zb, zdup = zdup)
+
+golden$dgp_covs <- list(
+  x = x5, y = y5, t = t5, zlong = zlong, zb = zb, zdup = zdup,
+  x_ties = x5_ties,
+  configs = list(
+    covs_default        = run_estimate(y5, x5, covs = covs2),
+    covs_manual_h       = run_estimate(y5, x5, covs = covs2, h = 0.2),
+    covs_msetwo         = run_estimate(y5, x5, covs = covs2,
+                                       bwselect = "msetwo"),
+    covs_cercomb2       = run_estimate(y5, x5, covs = covs2,
+                                       bwselect = "cercomb2"),
+    covs_epa            = run_estimate(y5, x5, covs = covs2, kernel = "epa"),
+    covs_drop_collinear = run_estimate(y5, x5, covs = covs3),
+    covs_ties           = run_estimate(y5, x5_ties, covs = covs2),
+    fuzzy_covs          = run_estimate(y5, x5, covs = covs2, fuzzy = t5),
+    fuzzy_covs_sharpbw  = run_estimate(y5, x5, covs = covs2, fuzzy = t5,
+                                       sharpbw = TRUE)
+  )
+)
+
 senate_path <- "benchmarks/data/rdrobust_senate.csv"
 stopifnot(file.exists(senate_path))
 senate <- read.csv(senate_path)