igerber · igerber · Jul 12, 2026 · Jul 12, 2026 · Jul 12, 2026 · Jul 12, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -30,6 +30,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   main — no `ready-for-ci` label needed; the aggregate `Lint Gate` job is the single
   required-check name. Plus `.git-blame-ignore-revs` covering the 2026-07 bulk
   normalization commits, and lint sections in CONTRIBUTING.md / CLAUDE.md.
+- **ChangesInChanges (Athey & Imbens 2006) + QDiD estimators.** Distributional
+  DiD for the canonical 2x2 design with continuous outcomes: `ChangesInChanges`
+  (alias `CiC`) recovers the treated group's full counterfactual outcome
+  distribution and quantile treatment effects via the CDF transformation
+  `F_10(F_00^{-1}(F_01(y)))`, invariant to monotone outcome rescaling; `QDiD`
+  is the paper's additive quantile-DiD comparison estimator. Point estimation
+  matches R `qte` 1.3.1 exactly (CiC via an exact port of R's type-1 quantile
+  arithmetic incl. its fuzz handling - bit-exact golden parity; QDiD via qte's
+  type-7 form, documented in REGISTRY.md as a finite-sample deviation from the
+  paper's transformation, population-equivalent). Bootstrap-only inference
+  (panel unit-block or pooled repeated-cross-section resampling per qte's
+  schemes; seedable, `n_bootstrap=200` default; replicate-SD SEs, normal
+  intervals, sup-t uniform bands at qte's fixed 95% level). Diagnostics:
+  interior-range guard (eq. 17) with NaN inference outside, support-violation /
+  heavy-ties / QDiD-non-monotonicity warnings. Deferred and documented:
+  covariates, discrete-outcome bounds, analytical SEs, staggered designs.
+  Golden fixtures: `benchmarks/R/generate_qte_golden.R` (qte pinned 1.3.1) ->
+  `benchmarks/data/qte_golden.json`.
 
 ### Changed
 - **Internal: repo-wide lint normalization + pinned tooling.** black/ruff/mypy are now

diff --git a/README.md b/README.md
@@ -118,6 +118,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`.
 - [StaggeredTripleDifference](https://diff-diff.readthedocs.io/en/stable/api/staggered.html#staggeredtripledifference) - Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD with group-time ATT
 - [WooldridgeDiD](https://diff-diff.readthedocs.io/en/stable/api/wooldridge_etwfe.html) - Wooldridge (2023, 2025) ETWFE: saturated OLS, logit/Poisson QMLE (ASF-based ATT). Alias `ETWFE`.
 - [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html) - Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting), variance- or equally-weighted ATT, for absorbing or non-absorbing (reversible) treatment
+- [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html) - Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: full counterfactual distribution and quantile treatment effects via CDF transformation, plus the QDiD comparison estimator; bootstrap inference; R qte parity. Alias `CiC`
 - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings
 
 ## Diagnostics & Sensitivity

diff --git a/TODO.md b/TODO.md
@@ -45,6 +45,11 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 ### Testing / docs
 
+| Issue | Location | Origin | Effort | Priority |
+|-------|----------|--------|--------|----------|
+| ChangesInChanges/QDiD tutorial notebook (2x2 distributional walkthrough: QTE grid, interior range, uniform bands, CiC-vs-QDiD comparison) - deferred from the implementation PR as a documented decision. | `docs/tutorials/` | #682 | Mid | Low |
+| `practitioner_next_steps()` dedicated handler for `ChangesInChangesResults` (currently falls back to `_handle_generic`, which is safe; a dedicated handler is the established full-integration step, cf. HAD Phase 5). | `diff_diff/practitioner.py` | #682 | Quick | Low |
+
 ---
 
 ## Deferred / Documented
@@ -105,6 +110,11 @@ Doable in principle, but no current caller and/or explicitly out of paper scope.
 
 | Issue | Location | PR | Priority |
 |-------|----------|----|----------|
+| ChangesInChanges covariates (Melly-Santangelo 2015 QR pipeline: per-cell quantile-regression conditional CDFs -> conditional CiC -> integrate over treated-group covariates). No R parity target exists (the MS Stata code is the only implementation; distinct from Kranker's `cic`); would need simulation-based validation. Reviewed: `docs/methodology/papers/melly-santangelo-2015-review.md`. | `diff_diff/changes_in_changes.py` | #682 | Low |
+| ChangesInChanges discrete-outcome bounds + DCIC point identification (Athey-Imbens Sections 4/5.2 incl. Imbens-Manski intervals; Kranker's Stata `cic` is the reference). The shipped ties warning marks the boundary of the continuous scope. | `diff_diff/changes_in_changes.py` | #682 | Low |
+| ChangesInChanges analytical SEs (Athey-Imbens Theorems 5.1-5.3 influence functions, panel 5.5-5.7, Appendix B covariances; needs the footnote-31 boundary density estimator - note the review's suspected half-range/midpoint typo). Bootstrap is the shipped inference. | `diff_diff/changes_in_changes.py` | #682 | Low |
+| Staggered/multi-period distributional DiD (Athey-Imbens Section 6 / Ciaccio arXiv:2408.01208v2; `ecic` is the staggered event-study CiC lineage - a distinct method from Ciaccio's copula approach, do not conflate). Reviewed: `docs/methodology/papers/ciaccio-2024-review.md`; ROADMAP row is reviewed-deferred pending demand. | `diff_diff/changes_in_changes.py` | #682 | Low |
+| ChangesInChanges treatment-on-the-controls (Athey-Imbens Theorem 3.2: group-label exchange + negation; no qte equivalent to anchor parity). | `diff_diff/changes_in_changes.py` | #682 | Low |
 | Rust-backend CR2 Bell-McCaffrey port (`return_dof` in the Rust vcov dispatch + CR2 algebra) — **premise re-scoped 2026-07-09**: the scores-based DOF + low-rank factored `A_g` changes made the NumPy CR2-BM path BLAS-bound (`O(n_g k²)` per cluster; 4.1s→38ms at n=100k/k=40), so a Rust port buys ~nothing and adds a parity surface. Revisit only if profiling shows CR2-BM hot again. | `rust/src/linalg.rs` | — | Low |
 | Clustered-CR1 inference df **default flip to `"cluster"` (G−1) at v4** — the opt-in `df_convention=` knob landed 2026-07 (DiD/TWFE/MPD + LinearRegression; REGISTRY §TwoWayFixedEffects deviation note); the remaining work is the major-version default change (moves every clustered p-value/CI) + migration note + flipping `TestDfConvention`/`test_moderate_t_pins_residual_df_convention` expectations. Also evaluate extending the knob to standalone estimators with CR1-t inference at that time. | `diff_diff/linalg.py::LinearRegression`, `diff_diff/estimators.py`, `diff_diff/twfe.py` | — | Medium |
 | CallawaySantAnna **unbalanced-panel R parity — LANDED** via `allow_unbalanced_panel=True` (matches R `did::att_gt(allow_unbalanced_panel=TRUE)` / `DRDID::reg_did_rc`: ATT bit-exact on cells AND dynamic aggregation via fixed unit-cohort-mass `pg` + a per-unit WIF; SE up to the documented CR1 `sqrt(G/(G-1))` factor). The earlier "weighting" framing was a mis-diagnosis — on unbalanced panels the dominant divergence from R is the *estimator* (within-cell differencing vs RC-on-pooled-obs), not only the weighting; both are resolved by the flag. The DEFAULT path keeps within-cell differencing as a documented design choice and now emits a `UserWarning` on unbalanced input (no-silent-failures). **Remaining deferred:** `survey_design=` × `allow_unbalanced_panel=` (per-obs vs per-unit weight resolution — currently fail-closed `NotImplementedError`); and covariate / ipw / dr × the flag R-parity verification (the RC path supports them; the committed golden covers `reg` no-cov). | `staggered.py`, `staggered_aggregation.py` | SE-audit D3 | Low |

diff --git a/benchmarks/R/generate_qte_golden.R b/benchmarks/R/generate_qte_golden.R
@@ -0,0 +1,214 @@
+#!/usr/bin/env Rscript
+# Golden fixtures for the diff-diff ChangesInChanges / QDiD estimators vs qte 1.3.1.
+#
+# Regenerate with:
+#   Rscript benchmarks/R/generate_qte_golden.R
+#
+# Output: benchmarks/data/qte_golden.json, consumed by
+# tests/test_changes_in_changes_parity.py (which pytest.skips when absent).
+#
+# Version contract: qte is pinned to 1.3.1 here, in benchmarks/R/requirements.R
+# (pinned_versions), and in the parity test's test_metadata_versions_match.
+# Bump all three in lockstep when re-anchoring.
+#
+# Point fixtures are generated with se=FALSE so they are fully deterministic
+# (qte's bootstrap is not seedable through its public API). The SE block runs
+# seeded (set.seed before each call) with iters=999; Python compares those SEs
+# statistically, not bit-exactly.
+
+suppressMessages({
+  library(qte)
+  library(jsonlite)
+})
+
+stopifnot(packageVersion("qte") == "1.3.1")
+
+probs <- seq(0.05, 0.95, 0.05)
+
+# ---------------------------------------------------------------------------
+# Synthetic DGPs (long format: id, period 0/1, treat 0/1 group indicator, y)
+# ---------------------------------------------------------------------------
+
+make_normal_2x2 <- function(n_treated, n_control, seed) {
+  # Additive time shift + heterogeneous treatment effect monotone in the unit
+  # unobservable; full support overlap so the whole grid is interior.
+  set.seed(seed)
+  n <- n_treated + n_control
+  treat <- c(rep(1L, n_treated), rep(0L, n_control))
+  u <- rnorm(n, mean = 0, sd = 1)
+  y_pre <- u + rnorm(n, sd = 0.3)
+  effect <- 0.8 + 0.4 * pnorm(u)
+  y_post <- u + 0.5 + rnorm(n, sd = 0.3) + treat * effect
+  data.frame(
+    id = rep(seq_len(n), times = 2),
+    period = rep(c(0L, 1L), each = n),
+    treat = rep(treat, times = 2),
+    y = c(y_pre, y_post)
+  )
+}
+
+make_lognormal_2x2 <- function(n_treated, n_control, seed) {
+  # Skewed outcome with a scale-type change over time: CiC and QDiD are
+  # numerically far apart, exercising the inverse-CDF composition away from
+  # the additive special case.
+  set.seed(seed)
+  n <- n_treated + n_control
+  treat <- c(rep(1L, n_treated), rep(0L, n_control))
+  u <- rnorm(n, mean = 0, sd = 0.6)
+  y_pre <- exp(u + rnorm(n, sd = 0.2))
+  y_post <- exp(u * 1.3 + 0.3 + rnorm(n, sd = 0.2)) + treat * (0.5 + 0.5 * exp(u / 2))
+  data.frame(
+    id = rep(seq_len(n), times = 2),
+    period = rep(c(0L, 1L), each = n),
+    treat = rep(treat, times = 2),
+    y = c(y_pre, y_post)
+  )
+}
+
+dgps <- list(
+  normal_2x2_n500 = make_normal_2x2(250, 250, seed = 20260712),
+  lognormal_2x2_n300 = make_lognormal_2x2(150, 150, seed = 20260713),
+  # Small unbalanced-cell design stressing near-integer n*p type-1 boundaries.
+  smalln_2x2_n60 = make_normal_2x2(35, 25, seed = 20260714)
+)
+
+# lalonde (qte-shipped): 1975/1978 two-period panel; re78 zeros give heavy ties.
+data(lalonde, package = "qte")
+lal <- lalonde.psid.panel[lalonde.psid.panel$year %in% c(1975, 1978), ]
+dgps$lalonde_psid <- data.frame(
+  id = lal$id,
+  period = ifelse(lal$year == 1978, 1L, 0L),
+  treat = as.integer(lal$treat),
+  y = as.numeric(lal$re)
+)
+
+# ---------------------------------------------------------------------------
+# Point fixtures: {CiC, QDiD} x {panel, repeated cross-section}, se = FALSE
+# ---------------------------------------------------------------------------
+
+run_point <- function(df, method, panel) {
+  fn <- if (method == "cic") CiC else QDiD
+  res <- fn(
+    formla = y ~ treat,
+    t = 1, tmin1 = 0, tname = "period",
+    data = df,
+    panel = panel,
+    idname = if (panel) "id" else NULL,
+    se = FALSE,
+    probs = probs
+  )
+  list(ate = res$ate, qte = as.numeric(res$qte))
+}
+
+scenarios <- list()
+for (dgp_name in names(dgps)) {
+  df <- dgps[[dgp_name]]
+  entry <- list(
+    data = list(
+      id = df$id,
+      period = df$period,
+      treat = df$treat,
+      y = df$y
+    ),
+    results = list()
+  )
+  for (method in c("cic", "qdid")) {
+    for (panel in c(TRUE, FALSE)) {
+      key <- sprintf("%s_%s", method, if (panel) "panel" else "rcs")
+      entry$results[[key]] <- run_point(df, method, panel)
+      message(sprintf(
+        "%s / %s: ate = %.10g", dgp_name, key, entry$results[[key]]$ate
+      ))
+    }
+  }
+  scenarios[[dgp_name]] <- entry
+}
+
+# ---------------------------------------------------------------------------
+# SE block (statistical parity): normal_2x2_n500 only, seeded, iters = 999
+# ---------------------------------------------------------------------------
+
+run_se <- function(df, method, panel, seed) {
+  fn <- if (method == "cic") CiC else QDiD
+  set.seed(seed)
+  res <- fn(
+    formla = y ~ treat,
+    t = 1, tmin1 = 0, tname = "period",
+    data = df,
+    panel = panel,
+    idname = if (panel) "id" else NULL,
+    se = TRUE, iters = 999,
+    probs = probs
+  )
+  list(
+    ate = res$ate, qte = as.numeric(res$qte),
+    ate_se = res$ate.se, qte_se = as.numeric(res$qte.se),
+    sup_t_crit = as.numeric(res$c)
+  )
+}
+
+se_block <- list(
+  cic_panel = run_se(dgps$normal_2x2_n500, "cic", TRUE, seed = 42),
+  cic_rcs = run_se(dgps$normal_2x2_n500, "cic", FALSE, seed = 42),
+  qdid_panel = run_se(dgps$normal_2x2_n500, "qdid", TRUE, seed = 42),
+  qdid_rcs = run_se(dgps$normal_2x2_n500, "qdid", FALSE, seed = 42)
+)
+message("SE block done.")
+
+# ---------------------------------------------------------------------------
+# Micro-fixtures: raw R type-1 / type-7 quantile outputs on adversarial probs
+# ---------------------------------------------------------------------------
+
+micro_case <- function(x, p) {
+  list(
+    x = x,
+    probs = p,
+    type1 = as.numeric(quantile(x, probs = p, type = 1, names = FALSE)),
+    type7 = as.numeric(quantile(x, probs = p, names = FALSE))
+  )
+}
+
+set.seed(20260715)
+x7 <- sort(rnorm(7))
+x20 <- sort(rnorm(20))
+x60 <- sort(rlnorm(60))
+x101 <- sort(rnorm(101))
+# ECDF-composed probabilities (the parity-critical case): ranks k/n00 from one
+# sample evaluated in another - float products can land an ulp off an integer.
+p_comp20 <- ecdf(x7)(x20)
+p_comp60 <- ecdf(x20)(x60)
+grid_edges <- c(0, 1e-12, 1 / 7, 2 / 7, 0.25, 1 / 3, 0.5, 19 / 20, 0.95, 0.999, 1)
+p_near_int101 <- c((0:101) / 101, (1:100) / 101 + 1e-16, (1:100) / 101 - 1e-16)
+
+micro <- list(
+  n7_grid = micro_case(x7, grid_edges),
+  n20_composed = micro_case(x20, p_comp20),
+  n60_composed = micro_case(x60, p_comp60),
+  n101_near_integer = micro_case(x101, pmin(pmax(p_near_int101, 0), 1))
+)
+
+# ---------------------------------------------------------------------------
+# Write JSON
+# ---------------------------------------------------------------------------
+
+out <- list(
+  metadata = list(
+    description = paste(
+      "Golden fixtures for diff-diff ChangesInChanges/QDiD parity vs qte.",
+      "Point fixtures use se=FALSE (deterministic); the se_block is seeded",
+      "(set.seed before each call, iters=999) and compared statistically."
+    ),
+    qte_version = as.character(packageVersion("qte")),
+    r_version = as.character(getRversion()),
+    probs = probs,
+    point_atol = 1e-10,
+    n_scenarios = length(scenarios)
+  ),
+  scenarios = scenarios,
+  se_block = se_block,
+  quantile_cases = micro
+)
+
+out_path <- file.path("benchmarks", "data", "qte_golden.json")
+write_json(out, path = out_path, digits = 17, auto_unbox = TRUE, null = "null")
+message(sprintf("Wrote %s (%.1f KB)", out_path, file.size(out_path) / 1024))
diff --git a/benchmarks/R/requirements.R b/benchmarks/R/requirements.R
@@ -18,6 +18,7 @@ required_packages <- c(
   "YatchewTest",   # Yatchew (1997) linearity test (HAD yatchew R-parity)
   "nprobust",      # Calonico-Cattaneo-Farrell local-linear (DIDHAD dependency)
   "Synth",         # Abadie-Diamond-Hainmueller (2010) synthetic control (SyntheticControl R-parity; ships data(basque))
+  "qte",           # Callaway qte package (Athey-Imbens CiC + QDiD R-parity; ships data(lalonde))
 
   # Utilities
   "jsonlite",      # JSON output for Python interop
@@ -71,7 +72,9 @@ install_pinned_version <- function(pkg, version) {
 pinned_versions <- list(
   DIDHAD = "2.0.0",
   YatchewTest = "1.1.1",
-  nprobust = "0.5.0"
+  nprobust = "0.5.0",
+  # CiC/QDiD R-parity (generate_qte_golden.R + tests/test_changes_in_changes_parity.py)
+  qte = "1.3.1"
 )
 
 install_github_if_missing <- function(pkg, repo) {

diff --git a/benchmarks/data/qte_golden.json b/benchmarks/data/qte_golden.json
diff --git a/diff_diff/__init__.py b/diff_diff/__init__.py
@@ -51,6 +51,14 @@
     ChaisemartinDHaultfoeuilleResults,
     DCDHBootstrapResults,
 )
+from diff_diff.changes_in_changes import (
+    ChangesInChanges,
+    QDiD,
+)
+from diff_diff.changes_in_changes_results import (
+    ChangesInChangesResults,
+    QDiDResults,
+)
 from diff_diff.continuous_did import (
     ContinuousDiD,
     ContinuousDiDResults,
@@ -299,6 +307,7 @@
 ETWFE = WooldridgeDiD
 DCDH = ChaisemartinDHaultfoeuille
 HAD = HeterogeneousAdoptionDiD
+CiC = ChangesInChanges
 
 __version__ = "3.7.0"
 __all__ = [
@@ -309,6 +318,8 @@
     "SyntheticDiD",
     "CallawaySantAnna",
     "ChaisemartinDHaultfoeuille",
+    "ChangesInChanges",
+    "QDiD",
     "ContinuousDiD",
     "SunAbraham",
     "ImputationDiD",
@@ -326,6 +337,7 @@
     "CS",
     "CDiD",
     "DCDH",
+    "CiC",
     "SA",
     "BJS",
     "Gardner",
@@ -346,6 +358,8 @@
     "CallawaySantAnnaResults",
     "CSBootstrapResults",
     "GroupTimeEffect",
+    "ChangesInChangesResults",
+    "QDiDResults",
     "ContinuousDiDResults",
     "DoseResponseCurve",
     "SunAbrahamResults",