diff --git a/README.md b/README.md index 3c5acfc07..98155765d 100644 --- a/README.md +++ b/README.md @@ -122,6 +122,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`. - [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 via `method="qdid"`; bootstrap inference; R qte parity. Alias `CiC` - [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html) - Lee & Wooldridge (2025, 2026) rolling-transformation DiD: unit-specific demean/detrend converts panel to cross-section, staggered adoption, `estimation_method` in `reg`/`ipw`/`dr`/`psm` (the papers' RA/IPW/IPWRA plus propensity-score matching), exact small-N inference on the classical collapsed regression - [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html) - Chang (2020) double/debiased machine learning DiD: staggered ATT(g,t) with cross-fitted ML nuisance learners (DML2) and Neyman-orthogonal scores, for flexible/high-dimensional covariate adjustment under conditional parallel trends; panel or declared repeated cross sections (`panel=False`); survey/cluster support on both lanes (bad-control lane: panel only, `cluster=` only); Caetano, Callaway, Payne & Sant'Anna (2026) bad-control score via `fit(bad_control=, bad_control_covariates=)` +- [DIDOVBSensitivity](https://diff-diff.readthedocs.io/en/stable/api/did_ovb.html) - Wang, Sant'Anna, Chernozhukov & Cinelli (2026) canonical two-period DiD omitted-variable-bias sensitivity analysis with cross-fitted nuisance learners, restricted bias bounds, robustness values, and extreme robustness values - [DurationDiD](https://diff-diff.readthedocs.io/en/stable/api/duration_did.html) - Deaner & Ku (2026) causal duration DiD for a binary absorbing outcome (spell ended) in a two-group common-timing design: restricts the groups' untreated hazards (`method="cd"` additive gap or `method="ph"` ratio) instead of outcome levels, imputes the treated counterfactual survival, reports the per-date absorption ATT with whole-individual bootstrap pointwise and simultaneous bands plus a fixed-anchor pre-treatment specification test - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings diff --git a/benchmarks/R/benchmark_did_ovb_minwage.R b/benchmarks/R/benchmark_did_ovb_minwage.R new file mode 100644 index 000000000..b6b4921ed --- /dev/null +++ b/benchmarks/R/benchmark_did_ovb_minwage.R @@ -0,0 +1,134 @@ +#!/usr/bin/env Rscript + +# R application oracle for Wang et al.'s minimum-wage example. +# This file uses ranger only as an executed reference implementation; no +# ranger source is copied into diff-diff. + +suppressPackageStartupMessages({ + library(jsonlite) + library(ranger) +}) + +args <- commandArgs(trailingOnly = TRUE) +input <- if (length(args) >= 1) args[[1]] else "/private/tmp/CS_RR/data/min_wage_CS.rds" +output <- if (length(args) >= 2) args[[2]] else "/private/tmp/did_ovb_minwage_r.json" +fold_output <- if (length(args) >= 3) args[[3]] else "/private/tmp/did_ovb_minwage_folds.csv" +seed <- if (length(args) >= 4) as.integer(args[[4]]) else 42L +n_folds <- if (length(args) >= 5) as.integer(args[[5]]) else 5L +num_trees <- if (length(args) >= 6) as.integer(args[[6]]) else 1000L +mtry <- if (length(args) >= 7) as.integer(args[[7]]) else 2L +min_node_size <- if (length(args) >= 8) as.integer(args[[8]]) else 10L +splitrule <- if (length(args) >= 9) args[[9]] else "variance" + +if (grepl("\\.rds$", tolower(input))) { + raw <- readRDS(input) +} else { + raw <- read.csv(input, stringsAsFactors = TRUE, check.names = FALSE) +} + +dat <- raw[raw$year %in% c(2006, 2007) & raw$first.treat %in% c(0, 2007), , drop = FALSE] +dat$treated <- as.integer(dat$first.treat == 2007) +dat <- dat[order(dat$countyreal, dat$year), , drop = FALSE] + +units <- unique(dat$countyreal) +if (length(units) != 1961L || sum(dat$treated[dat$year == 2006]) != 584L) { + stop("unexpected minimum-wage application sample") +} + +wide <- dat[dat$year == 2006, c("countyreal", "treated", "region", "white", "hs", "pov", "lpop", "lmedinc"), drop = FALSE] +post <- dat[dat$year == 2007, c("countyreal", "lemp"), drop = FALSE] +pre <- dat[dat$year == 2006, c("countyreal", "lemp"), drop = FALSE] +wide <- merge(wide, pre, by = "countyreal", suffixes = c("", "_pre"), sort = FALSE) +wide <- merge(wide, post, by = "countyreal", suffixes = c("", "_post"), sort = FALSE) +wide$delta_y <- wide$lemp_post - wide$lemp + +set.seed(seed) +folds <- integer(nrow(wide)) +for (d in c(0L, 1L)) { + idx <- which(wide$treated == d) + folds[idx] <- sample(rep(0:(n_folds - 1L), length.out = length(idx))) +} +write.csv(data.frame(countyreal = wide$countyreal, fold_id = folds), fold_output, row.names = FALSE) + +x_names <- c("region", "white", "hs", "pov", "lpop", "lmedinc") +x <- wide[, x_names, drop = FALSE] +# The Python public API currently accepts numeric covariates. Encode region by +# its public integer levels so both application runners consume identical X. +x$region <- as.numeric(x$region) +x <- as.matrix(x) +d <- wide$treated +y <- wide$delta_y +ps <- numeric(nrow(wide)) +m <- numeric(nrow(wide)) + +for (fold in 0:(n_folds - 1L)) { + test <- folds == fold + train <- !test + ps_fit <- ranger( + x = as.data.frame(x[train, , drop = FALSE]), + y = d[train], + num.trees = num_trees, + mtry = mtry, + min.node.size = min_node_size, + splitrule = splitrule, + seed = seed + fold + ) + out_fit <- ranger( + x = as.data.frame(x[train & d == 0, , drop = FALSE]), + y = y[train & d == 0], + num.trees = num_trees, + mtry = mtry, + min.node.size = min_node_size, + splitrule = splitrule, + seed = seed + fold + ) + ps[test] <- predict(ps_fit, data = as.data.frame(x[test, , drop = FALSE]))$predictions + m[test] <- predict(out_fit, data = as.data.frame(x[test, , drop = FALSE]))$predictions +} + +ps <- pmin(pmax(ps, 0.01), 0.99) +p <- mean(d) +omega <- (ps / (1 - ps)) / (p / (1 - p)) +residual <- y - m +score <- (d / p - (1 - d) / (1 - p) * omega) * residual +short_att <- mean(score) +sigma2 <- mean(residual[d == 0]^2) +nu2 <- mean(omega[d == 0]^2) +scale <- sqrt(sigma2 * nu2) +theta_if <- score - short_att +short_se <- sqrt(mean(theta_if^2) / length(y)) + +# This is the same plug-in/IF root search used by the Python result object. +sigma_if <- (1 - d) / (1 - p) * (residual^2 - sigma2) +nu_if <- (1 - d) / (1 - p) * (omega^2 - nu2) +scale_if <- (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale) +contains <- function(strength, kind, alpha = 0.05) { + multiplier <- if (kind == "rv") strength / sqrt(1 - strength) else sqrt(strength / (1 - strength)) + radius <- multiplier * scale + lo_if <- theta_if - multiplier * scale_if + hi_if <- theta_if + multiplier * scale_if + z <- qnorm(1 - alpha / 2) + lo_se <- sqrt(mean(lo_if^2) / length(y)) + hi_se <- sqrt(mean(hi_if^2) / length(y)) + short_att - radius - z * lo_se <= 0 && 0 <= short_att + radius + z * hi_se +} +find_rv <- function(kind) { + if (!contains(1 - 1e-12, kind)) return(NaN) + lo <- 0 + hi <- 1 - 1e-12 + for (i in seq_len(70)) { + mid <- (lo + hi) / 2 + if (contains(mid, kind)) hi <- mid else lo <- mid + } + hi +} + +result <- list( + n = nrow(wide), treated = sum(d), control = sum(d == 0), n_folds = n_folds, + seed = seed, num_trees = num_trees, mtry = mtry, + min_node_size = min_node_size, splitrule = splitrule, + short_att = short_att, short_se = short_se, sigma2_control = sigma2, + nu2_selection = nu2, scale = scale, rv = find_rv("rv"), xrv = find_rv("xrv") +) +write_json(result, output, auto_unbox = TRUE, digits = 17, pretty = TRUE) +cat(toJSON(result, auto_unbox = TRUE, digits = 8, pretty = TRUE), "\n") diff --git a/benchmarks/R/generate_did_ovb_parity.R b/benchmarks/R/generate_did_ovb_parity.R new file mode 100644 index 000000000..840dcec24 --- /dev/null +++ b/benchmarks/R/generate_did_ovb_parity.R @@ -0,0 +1,129 @@ +#!/usr/bin/env Rscript + +# Independent parity oracle for the canonical DiD OVB implementation. +# The calculations below are written from the paper's displayed formulas; +# this script does not import or call dml.sensemakr. + +suppressPackageStartupMessages(library(jsonlite)) + +args <- commandArgs(trailingOnly = TRUE) +output_path <- if (length(args) >= 1) args[[1]] else "benchmarks/data/did_ovb_r_results.json" +panel_path <- if (length(args) >= 2) args[[2]] else "benchmarks/data/real/mpdta.csv" + +zcrit <- qnorm(0.975) +data <- read.csv(panel_path, check.names = FALSE) +data <- data[data$year %in% c(2006, 2007) & data$`first.treat` %in% c(0, 2007), ] +data <- data[order(data$countyreal, data$year), ] + +pre <- data[data$year == 2006, ] +post <- data[data$year == 2007, ] +stopifnot(nrow(pre) == nrow(post), all(pre$countyreal == post$countyreal)) + +d <- as.numeric(pre$`first.treat` == 2007) +dy <- post$lemp - pre$lemp +x <- pre$lpop +n <- length(d) +n_folds <- 2L +fold_ids <- (seq_len(n) - 1L) %% n_folds +p <- mean(d) + +ps <- numeric(n) +m <- numeric(n) +for (fold in 0:(n_folds - 1L)) { + test <- fold_ids == fold + train <- !test + ps_fit <- glm(d ~ x, family = binomial(), subset = train) + m_fit <- lm(dy ~ x, subset = train & d == 0) + ps[test] <- predict(ps_fit, newdata = data.frame(x = x[test]), type = "response") + m[test] <- predict(m_fit, newdata = data.frame(x = x[test])) +} + +ps <- pmin(pmax(ps, 0.01), 0.99) +omega <- (ps / (1 - ps)) / (p / (1 - p)) +residual <- dy - m +score <- (d / p - (1 - d) / (1 - p) * omega) * residual +short_att <- mean(score) +sigma2 <- mean(residual[d == 0]^2) +nu2 <- mean(omega[d == 0]^2) +scale <- sqrt(sigma2 * nu2) +theta_if <- score - short_att +sigma_if <- (1 - d) / (1 - p) * (residual^2 - sigma2) +nu_if <- (1 - d) / (1 - p) * (omega^2 - nu2) +scale_if <- (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale) +short_se <- sqrt(mean(theta_if^2) / n) + +bounds <- function(trend_r2, selection_r2, rho_max = 1, alpha = 0.05) { + multiplier <- rho_max * sqrt(trend_r2) * sqrt(selection_r2 / (1 - selection_r2)) + radius <- multiplier * scale + lower_if <- theta_if - multiplier * scale_if + upper_if <- theta_if + multiplier * scale_if + lower_se <- sqrt(mean(lower_if^2) / n) + upper_se <- sqrt(mean(upper_if^2) / n) + z <- qnorm(1 - alpha / 2) + list( + lower = short_att - radius, + upper = short_att + radius, + radius = radius, + lower_se = lower_se, + upper_se = upper_se, + lower_ci = short_att - radius - z * lower_se, + upper_ci = short_att + radius + z * upper_se, + trend_r2 = trend_r2, + selection_r2 = selection_r2, + rho_max = rho_max, + alpha = alpha + ) +} + +robustness <- function(null_value = 0, alpha = 0.05) { + base <- bounds(0, 0, alpha = alpha) + if (base$lower_ci <= null_value && null_value <= base$upper_ci) { + return(list(rv = 0, xrv = 0, rv_bounds = base, xrv_bounds = base)) + } + search <- function(kind) { + contains <- function(s) { + b <- if (kind == "rv") bounds(s, s, alpha = alpha) else bounds(1, s, alpha = alpha) + b$lower_ci <= null_value && null_value <= b$upper_ci + } + lo <- 0 + hi <- 1 - 1e-12 + if (!contains(hi)) { + b <- if (kind == "rv") bounds(hi, hi, alpha = alpha) else bounds(1, hi, alpha = alpha) + return(list(value = NaN, bounds = b)) + } + for (i in seq_len(80)) { + mid <- (lo + hi) / 2 + if (contains(mid)) hi <- mid else lo <- mid + } + b <- if (kind == "rv") bounds(hi, hi, alpha = alpha) else bounds(1, hi, alpha = alpha) + list(value = hi, bounds = b) + } + rv <- search("rv") + xrv <- search("xrv") + list(rv = rv$value, xrv = xrv$value, rv_bounds = rv$bounds, xrv_bounds = xrv$bounds) +} + +result <- list( + settings = list( + panel = panel_path, + periods = c(2006, 2007), + treatment_cohort = 2007, + n_folds = n_folds, + fold_ids = fold_ids, + pscore_trim = 0.01, + alpha = 0.05, + null_value = -0.1 + ), + n_obs = n, + n_treated = sum(d == 1), + n_control = sum(d == 0), + short_att = short_att, + short_se = short_se, + sigma2_control = sigma2, + nu2_selection = nu2, + scale = scale, + bounds = bounds(1, 0.5), + robustness = robustness(-0.1, 0.05) +) + +write_json(result, output_path, auto_unbox = TRUE, digits = 17, pretty = TRUE) diff --git a/benchmarks/R/generate_did_ovb_simulation.R b/benchmarks/R/generate_did_ovb_simulation.R new file mode 100644 index 000000000..c9f2d0e0c --- /dev/null +++ b/benchmarks/R/generate_did_ovb_simulation.R @@ -0,0 +1,58 @@ +#!/usr/bin/env Rscript + +# Generate one draw from Appendix E.1 of Wang et al. and estimate the +# short DiD OVB components using the correctly specified parametric nuisances. +# The CSV is intentionally shared with the Python runner for cross-language +# comparison; it contains the latent U only because this is a simulation. + +suppressPackageStartupMessages(library(jsonlite)) + +args <- commandArgs(trailingOnly = TRUE) +csv_path <- if (length(args) >= 1) args[[1]] else "benchmarks/data/did_ovb_simulation.csv" +json_path <- if (length(args) >= 2) args[[2]] else "benchmarks/data/did_ovb_simulation_r.json" +n <- if (length(args) >= 3) as.integer(args[[3]]) else 500L +p <- if (length(args) >= 4) as.numeric(args[[4]]) else 0.5 +seed <- if (length(args) >= 5) as.integer(args[[5]]) else 20260928L + +set.seed(seed) +d <- rbinom(n, 1L, p) +x <- rnorm(n, mean = ifelse(d == 0, 0.3, 0), sd = ifelse(d == 0, sqrt(6), sqrt(3))) +u <- rnorm(n, mean = ifelse(d == 0, 0.3, 0), sd = ifelse(d == 0, sqrt(6), sqrt(3))) +delta_y <- 1 + x + u + 2 * d + rnorm(n, sd = sqrt(2)) + +fold_ids <- (seq_len(n) - 1L) %% 10L +ps <- numeric(n) +m <- numeric(n) +for (fold in 0:9) { + test <- fold_ids == fold + train <- !test + ps_fit <- glm(d ~ x + I(x^2), family = binomial(), subset = train) + m_fit <- lm(delta_y ~ x, subset = train & d == 0) + ps[test] <- predict(ps_fit, newdata = data.frame(x = x[test]), type = "response") + m[test] <- predict(m_fit, newdata = data.frame(x = x[test])) +} + +ps <- pmin(pmax(ps, 0.01), 0.99) +p_hat <- mean(d) +omega <- (ps / (1 - ps)) / (p_hat / (1 - p_hat)) +residual <- delta_y - m +score <- (d / p_hat - (1 - d) / (1 - p_hat) * omega) * residual +short_att <- mean(score) +sigma2 <- mean(residual[d == 0]^2) +nu2 <- mean(omega[d == 0]^2) +scale <- sqrt(sigma2 * nu2) + +write.csv( + data.frame(unit = seq_len(n), pre = 0, post = delta_y, treated = d, x = x, u = u), + csv_path, + row.names = FALSE +) +write_json( + list( + n = n, p = p, seed = seed, n_folds = 10L, + short_att = short_att, sigma2_control = sigma2, + nu2_selection = nu2, scale = scale, + treatment_effect = 2, omitted_bias = -mean(u[d == 1]) + mean(u[d == 0]) + ), + json_path, auto_unbox = TRUE, digits = 17, pretty = TRUE +) diff --git a/benchmarks/R/run_did_ovb_simulation_mc.R b/benchmarks/R/run_did_ovb_simulation_mc.R new file mode 100644 index 000000000..90a6e1572 --- /dev/null +++ b/benchmarks/R/run_did_ovb_simulation_mc.R @@ -0,0 +1,82 @@ +#!/usr/bin/env Rscript + +# Monte Carlo runner for the Appendix E.1 correctly specified parametric DGP. +# This reports coverage for the short estimand theta_s = 1.7. It is a compact +# reproducibility diagnostic; the paper's full tables additionally use true +# bias factors and several sensitivity-statistic surfaces. + +suppressPackageStartupMessages(library(jsonlite)) +args <- commandArgs(trailingOnly = TRUE) +output <- if (length(args) >= 1) args[[1]] else "benchmarks/data/did_ovb_simulation_mc_r.json" +n <- if (length(args) >= 2) as.integer(args[[2]]) else 500L +reps <- if (length(args) >= 3) as.integer(args[[3]]) else 20L +p <- if (length(args) >= 4) as.numeric(args[[4]]) else 0.5 +seed <- if (length(args) >= 5) as.integer(args[[5]]) else 20260928L +alpha <- if (length(args) >= 6) as.numeric(args[[6]]) else 0.05 +misspecified <- if (length(args) >= 7) as.logical(as.integer(args[[7]])) else FALSE + +one_rep <- function(n, p, seed) { + set.seed(seed) + d <- rbinom(n, 1L, p) + x <- rnorm(n, ifelse(d == 0, 0.3, 0), ifelse(d == 0, sqrt(6), sqrt(3))) + u <- rnorm(n, ifelse(d == 0, 0.3, 0), ifelse(d == 0, sqrt(6), sqrt(3))) + dy <- 1 + x + u + 2 * d + rnorm(n, sd = sqrt(2)) + folds <- (seq_len(n) - 1L) %% 10L + x_fit <- if (misspecified) exp(x / 2) else x + x_prop <- cbind(1, x_fit, x_fit^2) + x_out <- cbind(1, x_fit) + ps <- numeric(n); m <- numeric(n) + for (fold in 0:9) { + test <- folds == fold; train <- !test + ps_fit <- glm.fit(x_prop[train, , drop = FALSE], d[train], family = binomial()) + m_fit <- lm.fit(x_out[train & d == 0, , drop = FALSE], dy[train & d == 0]) + ps[test] <- plogis(x_prop[test, , drop = FALSE] %*% ps_fit$coefficients) + m[test] <- x_out[test, , drop = FALSE] %*% m_fit$coefficients + } + ps <- pmin(pmax(ps, 0.01), 0.99) + phat <- mean(d); omega <- (ps / (1 - ps)) / (phat / (1 - phat)) + residual <- dy - m + score <- (d / phat - (1 - d) / (1 - phat) * omega) * residual + att <- mean(score); se <- sqrt(mean((score - att)^2) / n) + sigma2 <- mean(residual[d == 0]^2); nu2 <- mean(omega[d == 0]^2) + scale <- sqrt(sigma2 * nu2) + theta_if <- score - att + sigma_if <- (1 - d) / (1 - phat) * (residual^2 - sigma2) + nu_if <- (1 - d) / (1 - phat) * (omega^2 - nu2) + scale_if <- (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale) + z <- qnorm(1 - alpha / 2) + contains <- function(strength, kind) { + multiplier <- if (kind == "rv") strength / sqrt(1 - strength) else sqrt(strength / (1 - strength)) + radius <- multiplier * scale + lower_if <- theta_if - multiplier * scale_if + upper_if <- theta_if + multiplier * scale_if + lower_se <- sqrt(mean(lower_if^2) / n); upper_se <- sqrt(mean(upper_if^2) / n) + att - radius - z * lower_se <= 0 && 0 <= att + radius + z * upper_se + } + find_rv <- function(kind) { + if (!contains(1 - 1e-12, kind)) return(NaN) + lo <- 0; hi <- 1 - 1e-12 + for (i in seq_len(70)) { + mid <- (lo + hi) / 2 + if (contains(mid, kind)) hi <- mid else lo <- mid + } + hi + } + c(att = att, se = se, cover = as.numeric(1.7 >= att - z * se && 1.7 <= att + z * se), + rv = find_rv("rv"), xrv = find_rv("xrv")) +} + +draws <- t(vapply(seq_len(reps), function(i) one_rep(n, p, seed + i - 1L), numeric(5))) +result <- list( + n = n, reps = reps, p = p, alpha = alpha, misspecified = misspecified, + seed = seed, true_short_att = 1.7, + mean_att = mean(draws[, "att"]), sd_att = sd(draws[, "att"]), + mean_se = mean(draws[, "se"]), bias = mean(draws[, "att"]) - 1.7, + coverage = mean(draws[, "cover"]), + mean_rv = mean(draws[, "rv"], na.rm = TRUE), + mean_xrv = mean(draws[, "xrv"], na.rm = TRUE), + rv_sd = sd(draws[, "rv"], na.rm = TRUE), + xrv_sd = sd(draws[, "xrv"], na.rm = TRUE) +) +write_json(result, output, auto_unbox = TRUE, digits = 17, pretty = TRUE) +cat(toJSON(result, auto_unbox = TRUE, digits = 8, pretty = TRUE), "\n") diff --git a/benchmarks/README.md b/benchmarks/README.md index e25f0b956..f4c1bf546 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -110,6 +110,95 @@ benchmarks/ ## Estimator Comparisons +## DiD OVB sensitivity parity + +`R/generate_did_ovb_parity.R` is an independent R oracle for the canonical +two-period omitted-variable-bias sensitivity implementation. It uses the +2006--2007 treated/never-treated slice of `data/real/mpdta.csv`, fixed folds, +and the paper's displayed plug-in and influence-function formulas. Regenerate +the JSON fixture with: + +```bash +Rscript benchmarks/R/generate_did_ovb_parity.R \ + benchmarks/data/did_ovb_r_results.json \ + benchmarks/data/real/mpdta.csv +``` + +The corresponding Python test is `tests/test_did_ovb_r_parity.py`; it compares +the short ATT, scale components, bounds, RV, and XRV. This is a clean-room +parity harness and does not import the GPL-3 `dml.sensemakr` source. + +The paper's minimum-wage application uses a separate external data file from +the authors' `CS_RR` repository. After exporting `data/min_wage_CS.rds` to CSV +in R, run the learner-level application check with: + +```bash +Rscript -e 'write.csv(readRDS("data/min_wage_CS.rds"), "min_wage_CS.csv", row.names=FALSE)' +python benchmarks/python/benchmark_did_ovb_minwage.py min_wage_CS.csv +``` + +This runner uses sklearn random forests as an explicit approximation to the +paper's tuned `ranger` specification. Its output is a diagnostic, not a paper +replication claim; exact application parity requires matching `ranger`, its +cross-validation choices, fold assignments, and multiplier-bootstrap settings. + +An R oracle using the actual `ranger` implementation is also provided. It +writes the fold assignment used by the R run so the Python estimator can use +the same units and folds: + +```bash +R_LIBS_USER=/path/to/r-library Rscript benchmarks/R/benchmark_did_ovb_minwage.R \ + /path/to/min_wage_CS.rds /tmp/did_ovb_minwage_r.json \ + /tmp/did_ovb_minwage_folds.csv 42 5 1000 2 10 variance +PYTHONPATH=. python benchmarks/python/benchmark_did_ovb_minwage.py \ + /path/to/min_wage_CS.csv --folds-csv /tmp/did_ovb_minwage_folds.csv \ + --n-estimators 1000 --max-features 2 --min-samples-leaf 10 \ + --seed 42 --n-folds 5 --output /tmp/did_ovb_minwage_python.json +``` + +The R and Python forests are deliberately reported as separate estimates: +`ranger` and sklearn do not implement identical tree-growth and probability +prediction rules. The shared fold file isolates that learner implementation +difference from sample construction and cross-fitting differences. The +published application estimate is approximately `-0.0366`; the Python +diagnostic with the fixed seed and default learner settings is approximately +`-0.0362` on the shared data. + +The Appendix E.1 simulation DGP has a shared R/Python fixture runner: + +```bash +Rscript benchmarks/R/generate_did_ovb_simulation.R \ + benchmarks/data/did_ovb_simulation.csv \ + benchmarks/data/did_ovb_simulation_r.json +PYTHONPATH=. python benchmarks/python/benchmark_did_ovb_simulation.py \ + benchmarks/data/did_ovb_simulation.csv \ + --r-json benchmarks/data/did_ovb_simulation_r.json +``` + +This is the first single-draw cross-language check. The full 5,000-repetition +coverage and sensitivity-statistics Monte Carlo tables remain a separate task. + +The compact Monte Carlo diagnostics are available as: + +```bash +Rscript benchmarks/R/run_did_ovb_simulation_mc.R /tmp/ovb_mc_r.json 500 20 +PYTHONPATH=. python benchmarks/python/run_did_ovb_simulation_mc.py \ + --n 500 --reps 20 --output /tmp/ovb_mc_python.json +``` + +Both runners report the mean short ATT, Monte Carlo standard deviation, mean +estimated standard error, bias relative to the true short ATT, and 95% CI +coverage. The paper-scale 5,000-repetition tables and sensitivity-statistic +coverage surfaces are not silently substituted by this compact diagnostic. + +For the Appendix E.2 misspecification design, pass `--misspecified` to Python +and `1` as the seventh R argument. Both then fit the nuisance models using +`X*=exp(X/2)` while generating outcomes from the original `X`. + +The Python-only forest diagnostic is +`benchmarks/python/run_did_ovb_simulation_forest.py`. It is not an R parity +test because the current R environment does not have `ranger` installed. + | diff-diff | Reference Package | Reference | Status | |-----------|-----------|-----------|--------| | `CallawaySantAnna` | `did::att_gt` | Callaway & Sant'Anna (2021) | ✓ Integrated | diff --git a/benchmarks/data/did_ovb_r_results.json b/benchmarks/data/did_ovb_r_results.json new file mode 100644 index 000000000..7d9b87090 --- /dev/null +++ b/benchmarks/data/did_ovb_r_results.json @@ -0,0 +1,63 @@ +{ + "settings": { + "panel": "benchmarks/data/real/mpdta.csv", + "periods": [2006, 2007], + "treatment_cohort": 2007, + "n_folds": 2, + "fold_ids": [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1], + "pscore_trim": 0.01, + "alpha": 0.05, + "null_value": -0.1 + }, + "n_obs": 440, + "n_treated": 131, + "n_control": 309, + "short_att": -0.0256949187178233, + "short_se": 0.0163569461630401, + "sigma2_control": 0.02696266378721, + "nu2_selection": 1.04962434592, + "scale": 0.168228024841021, + "bounds": { + "lower": -0.193922943558844, + "upper": 0.142533106123198, + "radius": 0.168228024841021, + "lower_se": 0.0226380821635521, + "upper_se": 0.0222126833838563, + "lower_ci": -0.238292769278465, + "upper_ci": 0.186069165555547, + "trend_r2": 1, + "selection_r2": 0.5, + "rho_max": 1, + "alpha": 0.05 + }, + "robustness": { + "rv": 0.217077391238398, + "xrv": 0.0567711175755117, + "rv_bounds": { + "lower": -0.0669667293686721, + "upper": 0.0155768919330254, + "radius": 0.0412718106508488, + "lower_se": 0.0168540192023375, + "upper_se": 0.0167145789785925, + "lower_ci": -0.1, + "upper_ci": 0.048336864747817, + "trend_r2": 0.217077391238398, + "selection_r2": 0.217077391238398, + "rho_max": 1, + "alpha": 0.05 + }, + "xrv_bounds": { + "lower": -0.0669667293686721, + "upper": 0.0155768919330254, + "radius": 0.0412718106508488, + "lower_se": 0.0168540192023375, + "upper_se": 0.0167145789785925, + "lower_ci": -0.1, + "upper_ci": 0.048336864747817, + "trend_r2": 1, + "selection_r2": 0.0567711175755117, + "rho_max": 1, + "alpha": 0.05 + } + } +} diff --git a/benchmarks/python/benchmark_did_ovb_minwage.py b/benchmarks/python/benchmark_did_ovb_minwage.py new file mode 100644 index 000000000..620407aab --- /dev/null +++ b/benchmarks/python/benchmark_did_ovb_minwage.py @@ -0,0 +1,224 @@ +#!/usr/bin/env python3 +"""Run the Wang et al. minimum-wage OVB application on exported RDS data. + +The authors distribute the source data as ``min_wage_CS.rds``. This runner +expects a user-created CSV export so that the repository does not redistribute +the external replication data. It intentionally exposes the learner settings +instead of claiming that sklearn's random forest is identical to R's ranger. +""" + +from __future__ import annotations + +import argparse +import json + +import numpy as np +import pandas as pd + +from diff_diff import DIDOVBSensitivity + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("csv", help="CSV export of CS_RR/data/min_wage_CS.rds") + parser.add_argument("--n-estimators", type=int, default=1000) + parser.add_argument("--max-features", type=int, default=2) + parser.add_argument("--min-samples-leaf", type=int, default=10) + parser.add_argument("--seed", type=int, default=42) + parser.add_argument("--n-folds", type=int, default=5) + parser.add_argument( + "--folds-csv", + help="Optional countyreal/fold_id CSV produced by the R application oracle", + ) + parser.add_argument("--output", help="Optional JSON output path") + parser.add_argument( + "--tune", + action="store_true", + help="tune RF/ExtraTrees hyperparameters within each outer fold", + ) + args = parser.parse_args() + + try: + from sklearn.ensemble import RandomForestRegressor + except ImportError as exc: # pragma: no cover - exercised by environments + raise SystemExit("Install the optional scikit-learn dependency first") from exc + + data = pd.read_csv(args.csv) + data = data.loc[data["year"].isin([2006, 2007]) & data["first.treat"].isin([0, 2007])].copy() + data["treated"] = (data["first.treat"] == 2007).astype(int) + + forest_kwargs = { + "n_estimators": args.n_estimators, + "max_features": args.max_features, + "min_samples_leaf": args.min_samples_leaf, + "random_state": args.seed, + "n_jobs": -1, + } + + from sklearn.base import BaseEstimator + from sklearn.ensemble import ExtraTreesRegressor, RandomForestRegressor + + class _RegressionPropensity(BaseEstimator): + def __init__( + self, + n_estimators=1000, + max_features=2, + min_samples_leaf=10, + random_state=42, + n_jobs=-1, + ): + self.n_estimators = n_estimators + self.max_features = max_features + self.min_samples_leaf = min_samples_leaf + self.random_state = random_state + self.n_jobs = n_jobs + + def fit(self, X, y): + self.model_ = RandomForestRegressor( + n_estimators=self.n_estimators, + max_features=self.max_features, + min_samples_leaf=self.min_samples_leaf, + random_state=self.random_state, + n_jobs=self.n_jobs, + ).fit(X, y) + return self + + def predict_proba(self, X): + probability = self.model_.predict(X) + return np.column_stack([1 - probability, probability]) + + propensity = _RegressionPropensity(**forest_kwargs) + outcome = RandomForestRegressor(**forest_kwargs) + if args.tune: + from sklearn.model_selection import GridSearchCV + + # This is the sklearn analogue of the paper's ranger grid. The + # estimator is cloned independently inside every outer fold by the + # existing cross-fitting machinery. + prop_grid = [ + { + "model": [_RegressionPropensity(**forest_kwargs)], + "model__max_features": [2, 4, 6], + "model__min_samples_leaf": [10, 15, 25, 50, 75, 100, 125, 150], + }, + { + "model": [ + _RegressionPropensity( + n_estimators=forest_kwargs["n_estimators"], + min_samples_leaf=forest_kwargs["min_samples_leaf"], + random_state=forest_kwargs["random_state"], + n_jobs=forest_kwargs["n_jobs"], + ) + ], + "model__max_features": [2, 4, 6], + "model__min_samples_leaf": [10, 15, 25, 50, 75, 100, 125, 150], + }, + ] + out_grid = [ + { + "model": [RandomForestRegressor(**forest_kwargs)], + "model__max_features": [2, 4, 6], + "model__min_samples_leaf": [10, 15, 25, 50, 75, 100, 125, 150], + }, + { + "model": [ExtraTreesRegressor(**forest_kwargs)], + "model__max_features": [2, 4, 6], + "model__min_samples_leaf": [10, 15, 25, 50, 75, 100, 125, 150], + }, + ] + # A one-step wrapper keeps the public estimator API unchanged while + # allowing GridSearchCV to switch between ranger-like split rules. + from sklearn.base import BaseEstimator + from sklearn.pipeline import Pipeline + + class _GridLearner(BaseEstimator): + def __init__(self, estimator=None, grid=None, scoring=None, cv=3): + self.estimator = estimator + self.grid = grid + self.scoring = scoring + self.cv = cv + + def fit(self, X, y): + params = self.grid or {} + self.search_ = GridSearchCV( + self.estimator, + params, + scoring=self.scoring, + cv=self.cv, + n_jobs=-1, + refit=True, + ).fit(X, y) + return self + + def predict(self, X): + return self.search_.predict(X) + + def predict_proba(self, X): + return self.search_.predict_proba(X) + + propensity = _GridLearner( + estimator=Pipeline([("model", _RegressionPropensity(**forest_kwargs))]), + grid=prop_grid, + scoring="neg_root_mean_squared_error", + cv=3, + ) + outcome = _GridLearner( + estimator=Pipeline([("model", RandomForestRegressor(**forest_kwargs))]), + grid=out_grid, + scoring="neg_root_mean_squared_error", + cv=3, + ) + + fit_kwargs = {} + if args.folds_csv: + folds = pd.read_csv(args.folds_csv) + if not {"countyreal", "fold_id"}.issubset(folds.columns): + raise SystemExit("--folds-csv must contain countyreal and fold_id columns") + unit_order = data.loc[data["year"] == 2006, "countyreal"].drop_duplicates() + fold_map = dict(zip(folds["countyreal"], folds["fold_id"], strict=True)) + try: + fit_kwargs["fold_ids"] = unit_order.map(fold_map).to_numpy(dtype=int) + except (TypeError, ValueError) as exc: + raise SystemExit("--folds-csv does not cover every application unit") from exc + + result = DIDOVBSensitivity( + n_folds=args.n_folds, + seed=args.seed, + propensity_learner=propensity, + outcome_learner=outcome, + ).fit( + data, + outcome="lemp", + treatment="treated", + time="year", + unit="countyreal", + covariates=["region", "white", "hs", "pov", "lpop", "lmedinc"], + **fit_kwargs, + ) + + print(result.summary()) + print("RV:", result.robustness_value().rv) + print("XRV:", result.robustness_value().xrv) + payload = { + "n": result.n_obs, + "treated": result.n_treated, + "control": result.n_control, + "n_folds": result.n_folds, + "seed": result.seed, + "short_att": result.short_att, + "short_se": result.short_se, + "sigma2_control": result.sigma2_control, + "nu2_selection": result.nu2_selection, + "scale": result.scale, + "rv": result.robustness_value().rv, + "xrv": result.robustness_value().xrv, + } + print(json.dumps(payload, indent=2)) + if args.output: + with open(args.output, "w", encoding="utf-8") as handle: + json.dump(payload, handle, indent=2) + handle.write("\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/python/benchmark_did_ovb_simulation.py b/benchmarks/python/benchmark_did_ovb_simulation.py new file mode 100644 index 000000000..2cd196218 --- /dev/null +++ b/benchmarks/python/benchmark_did_ovb_simulation.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +"""Estimate the shared Appendix E.1 simulation fixture in Python.""" + +from __future__ import annotations + +import argparse +import json + +import pandas as pd +from sklearn.linear_model import LinearRegression, LogisticRegression +from sklearn.pipeline import make_pipeline +from sklearn.preprocessing import PolynomialFeatures + +from diff_diff import DIDOVBSensitivity + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("csv") + parser.add_argument("--r-json") + args = parser.parse_args() + + wide = pd.read_csv(args.csv) + data = pd.concat( + [ + wide[["unit", "treated", "x"]].assign(time=0, outcome=wide["pre"]), + wide[["unit", "treated", "x"]].assign(time=1, outcome=wide["post"]), + ], + ignore_index=True, + ) + propensity = make_pipeline( + PolynomialFeatures(degree=2, include_bias=False), + LogisticRegression(max_iter=2000, solver="lbfgs"), + ) + outcome = LinearRegression() + result = DIDOVBSensitivity( + n_folds=10, + seed=0, + propensity_learner=propensity, + outcome_learner=outcome, + pscore_trim=0.01, + ).fit( + data, + outcome="outcome", + treatment="treated", + time="time", + unit="unit", + covariates=["x"], + fold_ids=(wide["unit"].to_numpy() - 1) % 10, + ) + output = result.to_dict() + print(json.dumps(output, indent=2)) + if args.r_json: + reference = json.loads(open(args.r_json).read()) + print("R short_att:", reference["short_att"]) + print("R scale:", reference["scale"]) + print("Python-R short_att:", result.short_att - reference["short_att"]) + print("Python-R scale:", result.scale - reference["scale"]) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/python/run_did_ovb_simulation_forest.py b/benchmarks/python/run_did_ovb_simulation_forest.py new file mode 100644 index 000000000..aad0f7b88 --- /dev/null +++ b/benchmarks/python/run_did_ovb_simulation_forest.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +"""Random-forest diagnostic for the Appendix E.2 misspecification DGP. + +This is intentionally a Python-only diagnostic because the current R runtime +does not have the paper's ``ranger`` package installed. It uses the same +cross-fitting score as the parametric runner and reports whether a flexible +learner reduces the misspecification bias. +""" + +from __future__ import annotations + +import argparse +import json + +import numpy as np +from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor + + +def one_rep(n: int, p: float, seed: int, alpha: float, trees: int) -> tuple[float, float, float]: + rng = np.random.default_rng(seed) + d = rng.binomial(1, p, n) + x = rng.normal(np.where(d == 0, 0.3, 0.0), np.where(d == 0, np.sqrt(6), np.sqrt(3))) + u = rng.normal(np.where(d == 0, 0.3, 0.0), np.where(d == 0, np.sqrt(6), np.sqrt(3))) + dy = 1 + x + u + 2 * d + rng.normal(0, np.sqrt(2), n) + z = np.exp(x / 2).reshape(-1, 1) + ps = np.empty(n) + m = np.empty(n) + folds = np.arange(n) % 10 + for fold in range(10): + test = folds == fold + train = ~test + prop = RandomForestClassifier( + n_estimators=trees, + max_features=1, + min_samples_leaf=10, + random_state=seed + fold, + n_jobs=-1, + ).fit(z[train], d[train]) + outcome = RandomForestRegressor( + n_estimators=trees, + max_features=1, + min_samples_leaf=10, + random_state=seed + fold, + n_jobs=-1, + ).fit(z[train & (d == 0)], dy[train & (d == 0)]) + ps[test] = prop.predict_proba(z[test])[:, 1] + m[test] = outcome.predict(z[test]) + ps = np.clip(ps, 0.01, 0.99) + phat = float(d.mean()) + omega = (ps / (1 - ps)) / (phat / (1 - phat)) + score = (d / phat - (1 - d) / (1 - phat) * omega) * (dy - m) + att = float(score.mean()) + se = float(np.sqrt(np.mean((score - att) ** 2) / n)) + from scipy.stats import norm + + zcrit = norm.ppf(1 - alpha / 2) + return att, se, float(att - zcrit * se <= 1.7 <= att + zcrit * se) + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--n", type=int, default=500) + parser.add_argument("--reps", type=int, default=100) + parser.add_argument("--p", type=float, default=0.5) + parser.add_argument("--alpha", type=float, default=0.05) + parser.add_argument("--trees", type=int, default=100) + parser.add_argument("--seed", type=int, default=20260928) + parser.add_argument("--output") + args = parser.parse_args() + draws = np.asarray( + [ + one_rep(args.n, args.p, args.seed + i, args.alpha, args.trees) + for i in range(args.reps) + ] + ) + result = { + "n": args.n, + "reps": args.reps, + "p": args.p, + "alpha": args.alpha, + "trees": args.trees, + "true_short_att": 1.7, + "mean_att": float(draws[:, 0].mean()), + "sd_att": float(draws[:, 0].std(ddof=1)), + "mean_se": float(draws[:, 1].mean()), + "bias": float(draws[:, 0].mean() - 1.7), + "coverage": float(draws[:, 2].mean()), + } + payload = json.dumps(result, indent=2) + print(payload) + if args.output: + with open(args.output, "w", encoding="utf-8") as handle: + handle.write(payload + "\n") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/python/run_did_ovb_simulation_mc.py b/benchmarks/python/run_did_ovb_simulation_mc.py new file mode 100644 index 000000000..951d9453a --- /dev/null +++ b/benchmarks/python/run_did_ovb_simulation_mc.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +"""Monte Carlo diagnostic for the Appendix E.1 parametric DGP.""" + +from __future__ import annotations + +import argparse +import json + +import numpy as np +from scipy.special import expit + +from diff_diff.linalg import solve_logit, solve_ols + + +def one_rep(n: int, p: float, seed: int, alpha: float, misspecified: bool) -> tuple[float, float, float]: + rng = np.random.default_rng(seed) + d = rng.binomial(1, p, n) + x = rng.normal(np.where(d == 0, 0.3, 0.0), np.where(d == 0, np.sqrt(6), np.sqrt(3))) + u = rng.normal(np.where(d == 0, 0.3, 0.0), np.where(d == 0, np.sqrt(6), np.sqrt(3))) + dy = 1 + x + u + 2 * d + rng.normal(0, np.sqrt(2), n) + x_fit = np.exp(x / 2) if misspecified else x + x_short = np.column_stack([x_fit]) + x_prop = np.column_stack([x_fit, x_fit**2]) + ps = np.empty(n) + m = np.empty(n) + folds = np.arange(n) % 10 + for fold in range(10): + test = folds == fold + train = ~test + beta_p, _ = solve_logit(x_prop[train], d[train], max_iter=25, tol=1e-8) + eta = np.column_stack([np.ones(test.sum()), x_prop[test]]) @ beta_p + ps[test] = expit(eta) + beta_m, _, _ = solve_ols( + np.column_stack([np.ones((np.sum(train & (d == 0)), 1)), x_short[train & (d == 0)]]), + dy[train & (d == 0)], + return_vcov=False, + ) + m[test] = np.column_stack([np.ones(test.sum()), x_short[test]]) @ beta_m + ps = np.clip(ps, 0.01, 0.99) + phat = float(d.mean()) + omega = (ps / (1 - ps)) / (phat / (1 - phat)) + residual = dy - m + score = (d / phat - (1 - d) / (1 - phat) * omega) * residual + att = float(score.mean()) + se = float(np.sqrt(np.mean((score - att) ** 2) / n)) + sigma2 = float(np.mean(residual[d == 0] ** 2)) + nu2 = float(np.mean(omega[d == 0] ** 2)) + scale = float(np.sqrt(sigma2 * nu2)) + theta_if = score - att + sigma_if = (1 - d) / (1 - phat) * (residual**2 - sigma2) + nu_if = (1 - d) / (1 - phat) * (omega**2 - nu2) + scale_if = (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale) + from scipy.stats import norm + + z = norm.ppf(1 - alpha / 2) + lo = att - z * se + hi = att + z * se + + def contains(strength: float, kind: str) -> bool: + if kind == "rv": + multiplier = strength / np.sqrt(1 - strength) + else: + multiplier = np.sqrt(strength / (1 - strength)) + radius = multiplier * scale + lower_if = theta_if - multiplier * scale_if + upper_if = theta_if + multiplier * scale_if + lower_se = np.sqrt(np.mean(lower_if**2) / n) + upper_se = np.sqrt(np.mean(upper_if**2) / n) + return att - radius - z * lower_se <= 0 <= att + radius + z * upper_se + + def rv(kind: str) -> float: + if contains(1 - 1e-12, kind): + lo_s, hi_s = 0.0, 1 - 1e-12 + for _ in range(70): + mid = (lo_s + hi_s) / 2 + if contains(mid, kind): + hi_s = mid + else: + lo_s = mid + return hi_s + return float("nan") + + return att, se, float(lo <= 1.7 <= hi), rv("rv"), rv("xrv") + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--n", type=int, default=500) + parser.add_argument("--reps", type=int, default=20) + parser.add_argument("--p", type=float, default=0.5) + parser.add_argument("--alpha", type=float, default=0.05) + parser.add_argument( + "--misspecified", + action="store_true", + help="fit nuisance models using X*=exp(X/2), as in Appendix E.2", + ) + parser.add_argument("--seed", type=int, default=20260928) + parser.add_argument("--output") + args = parser.parse_args() + if not 0 < args.alpha < 1: + raise SystemExit("--alpha must be in (0, 1)") + draws = np.asarray( + [ + one_rep(args.n, args.p, args.seed + i, args.alpha, args.misspecified) + for i in range(args.reps) + ] + ) + result = { + "n": args.n, + "reps": args.reps, + "p": args.p, + "alpha": args.alpha, + "misspecified": args.misspecified, + "seed": args.seed, + "true_short_att": 1.7, + "mean_att": float(draws[:, 0].mean()), + "sd_att": float(draws[:, 0].std(ddof=1)), + "mean_se": float(draws[:, 1].mean()), + "bias": float(draws[:, 0].mean() - 1.7), + "coverage": float(draws[:, 2].mean()), + "mean_rv": float(np.nanmean(draws[:, 3])), + "mean_xrv": float(np.nanmean(draws[:, 4])), + "rv_sd": float(np.nanstd(draws[:, 3], ddof=1)), + "xrv_sd": float(np.nanstd(draws[:, 4], ddof=1)), + } + payload = json.dumps(result, indent=2) + print(payload) + if args.output: + with open(args.output, "w", encoding="utf-8") as handle: + handle.write(payload + "\n") + + +if __name__ == "__main__": + main() diff --git a/changelog.d/20260928-did-ovb-application.md b/changelog.d/20260928-did-ovb-application.md new file mode 100644 index 000000000..8053ed163 --- /dev/null +++ b/changelog.d/20260928-did-ovb-application.md @@ -0,0 +1,5 @@ +### Added +- **DiD OVB application benchmarks**: add an R `ranger` oracle and a paired Python runner for the minimum-wage application, including shared cross-fitting folds and documented learner differences. + +### Documentation +- **Document the application replication path**: record the external data provenance, paper tuning grid, and the limits of R/Python random-forest numerical parity. diff --git a/diff_diff/__init__.py b/diff_diff/__init__.py index dcbb3ea77..b646903ed 100644 --- a/diff_diff/__init__.py +++ b/diff_diff/__init__.py @@ -96,6 +96,12 @@ run_all_placebo_tests, run_placebo_test, ) +from diff_diff.did_ovb import ( + DIDOVBBounds, + DIDOVBRobustness, + DIDOVBSensitivity, + DIDOVBSensitivityResults, +) from diff_diff.dml_did import DMLDiD from diff_diff.dml_did_results import DMLDiDResults from diff_diff.duration_did import DurationDiD diff --git a/diff_diff/did_ovb.py b/diff_diff/did_ovb.py new file mode 100644 index 000000000..3db8635fd --- /dev/null +++ b/diff_diff/did_ovb.py @@ -0,0 +1,425 @@ +"""Clean-room omitted-variable-bias sensitivity analysis for canonical DiD. + +This module implements the estimable components and sensitivity bounds described +by Wang, Sant'Anna, Chernozhukov, and Cinelli, "Omitted Variable Bias in +Difference-in-Differences Designs". It is an original Python implementation; +the GPL-3 ``dml.sensemakr`` source is not imported or translated. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any, Dict, Optional, Sequence + +import numpy as np +import pandas as pd +from scipy.stats import norm + +from diff_diff._crossfit import FoldAssignment, assign_folds, cross_fit_predict +from diff_diff._learners import make_learner + +__all__ = [ + "DIDOVBSensitivity", + "DIDOVBSensitivityResults", + "DIDOVBBounds", + "DIDOVBRobustness", +] + + +def _validate_probability(value: float, name: str, *, allow_one: bool = True) -> float: + value = float(value) + upper = 1.0 if allow_one else np.nextafter(1.0, 0.0) + if not np.isfinite(value) or value < 0.0 or value > upper: + right = "1" if allow_one else "1 (exclusive)" + raise ValueError(f"{name} must be in [0, {right}], got {value!r}") + return value + + +def _validate_alpha(alpha: float) -> float: + alpha = float(alpha) + if not np.isfinite(alpha) or not 0.0 < alpha < 0.5: + raise ValueError(f"alpha must be in (0, 0.5), got {alpha!r}") + return alpha + + +@dataclass(frozen=True) +class DIDOVBBounds: + """Bias bounds and delta-method confidence intervals.""" + + lower: float + upper: float + radius: float + lower_se: float + upper_se: float + lower_ci: float + upper_ci: float + rho_max: float + trend_r2: float + selection_r2: float + alpha: float + + def to_dict(self) -> Dict[str, Any]: + return { + "lower": self.lower, + "upper": self.upper, + "radius": self.radius, + "lower_se": self.lower_se, + "upper_se": self.upper_se, + "lower_ci": self.lower_ci, + "upper_ci": self.upper_ci, + "rho_max": self.rho_max, + "trend_r2": self.trend_r2, + "selection_r2": self.selection_r2, + "alpha": self.alpha, + } + + +@dataclass(frozen=True) +class DIDOVBRobustness: + """RV/XRV search results and the intervals at the reported strengths.""" + + null_value: float + alpha: float + rv: float + xrv: float + rv_bounds: DIDOVBBounds + xrv_bounds: DIDOVBBounds + + def to_dict(self) -> Dict[str, Any]: + return { + "null_value": self.null_value, + "alpha": self.alpha, + "rv": self.rv, + "xrv": self.xrv, + "rv_bounds": self.rv_bounds.to_dict(), + "xrv_bounds": self.xrv_bounds.to_dict(), + } + + +@dataclass +class DIDOVBSensitivityResults: + """Fitted canonical DiD OVB sensitivity results.""" + + short_att: float + short_se: float + sigma2_control: float + nu2_selection: float + scale: float + n_obs: int + n_treated: int + n_control: int + n_folds: int + seed: Optional[int] + propensity_learner: Any + outcome_learner: Any + _theta_if: np.ndarray + _scale_if: np.ndarray + + def _radius(self, trend_r2: float, selection_r2: float, rho_max: float) -> float: + _validate_probability(trend_r2, "trend_r2") + _validate_probability(selection_r2, "selection_r2", allow_one=False) + _validate_probability(rho_max, "rho_max") + selection_factor = np.sqrt(selection_r2 / (1.0 - selection_r2)) + return float(rho_max * np.sqrt(trend_r2) * selection_factor * self.scale) + + def bounds( + self, + *, + trend_r2: float = 1.0, + selection_r2: float = 0.5, + rho_max: float = 1.0, + alpha: float = 0.05, + ) -> DIDOVBBounds: + """Return OVB bounds under paper-scale trend/selection restrictions. + + ``trend_r2`` is the share of residual control-trend variation explained + by the omitted confounder. ``selection_r2`` is the residual share of + treatment-odds variation; its implied selection factor is + ``sqrt(selection_r2 / (1-selection_r2))``. + """ + alpha = _validate_alpha(alpha) + radius = self._radius(trend_r2, selection_r2, rho_max) + selection_factor = np.sqrt(selection_r2 / (1.0 - selection_r2)) + multiplier = float(rho_max * np.sqrt(trend_r2) * selection_factor) + lower_if = self._theta_if - multiplier * self._scale_if + upper_if = self._theta_if + multiplier * self._scale_if + lower_se = float(np.sqrt(np.mean(lower_if**2) / self.n_obs)) + upper_se = float(np.sqrt(np.mean(upper_if**2) / self.n_obs)) + critical = float(norm.ppf(1.0 - alpha / 2.0)) + return DIDOVBBounds( + lower=self.short_att - radius, + upper=self.short_att + radius, + radius=radius, + lower_se=lower_se, + upper_se=upper_se, + lower_ci=self.short_att - radius - critical * lower_se, + upper_ci=self.short_att + radius + critical * upper_se, + rho_max=float(rho_max), + trend_r2=float(trend_r2), + selection_r2=float(selection_r2), + alpha=alpha, + ) + + def robustness_value( + self, *, null_value: float = 0.0, alpha: float = 0.05, tolerance: float = 1e-10 + ) -> DIDOVBRobustness: + """Find the paper's common-strength RV and selection-only XRV.""" + alpha = _validate_alpha(alpha) + null_value = float(null_value) + base = self.bounds(trend_r2=0.0, selection_r2=0.0, alpha=alpha) + if base.lower_ci <= null_value <= base.upper_ci: + zero = self.bounds(trend_r2=0.0, selection_r2=0.0, alpha=alpha) + return DIDOVBRobustness(null_value, alpha, 0.0, 0.0, zero, zero) + + def search(kind: str) -> tuple[float, DIDOVBBounds]: + def interval_contains(strength: float) -> bool: + if kind == "rv": + b = self.bounds(trend_r2=strength, selection_r2=strength, alpha=alpha) + else: + b = self.bounds(trend_r2=1.0, selection_r2=strength, alpha=alpha) + return b.lower_ci <= null_value <= b.upper_ci + + high = np.nextafter(1.0, 0.0) + if not interval_contains(high): + if kind == "rv": + b = self.bounds(trend_r2=high, selection_r2=high, alpha=alpha) + else: + b = self.bounds(trend_r2=1.0, selection_r2=high, alpha=alpha) + return float("nan"), b + low = 0.0 + for _ in range(80): + mid = (low + high) / 2.0 + if interval_contains(mid): + high = mid + else: + low = mid + if high - low <= tolerance: + break + if kind == "rv": + b = self.bounds(trend_r2=high, selection_r2=high, alpha=alpha) + else: + b = self.bounds(trend_r2=1.0, selection_r2=high, alpha=alpha) + return float(high), b + + rv, rv_bounds = search("rv") + xrv, xrv_bounds = search("xrv") + return DIDOVBRobustness(null_value, alpha, rv, xrv, rv_bounds, xrv_bounds) + + def contour( + self, *, null_value: float = 0.0, alpha: float = 0.05, n_grid: int = 101 + ) -> pd.DataFrame: + """Return confidence-bound endpoints over a common-strength grid.""" + alpha = _validate_alpha(alpha) + if isinstance(n_grid, bool) or int(n_grid) != n_grid or n_grid < 2: + raise ValueError(f"n_grid must be an integer >= 2, got {n_grid!r}") + rows = [] + for strength in np.linspace(0.0, np.nextafter(1.0, 0.0), int(n_grid)): + rv = self.bounds(trend_r2=float(strength), selection_r2=float(strength), alpha=alpha) + xrv = self.bounds(trend_r2=1.0, selection_r2=float(strength), alpha=alpha) + rows.append( + { + "strength": float(strength), + "rv_lower": rv.lower_ci, + "rv_upper": rv.upper_ci, + "xrv_lower": xrv.lower_ci, + "xrv_upper": xrv.upper_ci, + } + ) + return pd.DataFrame(rows) + + def to_dict(self) -> Dict[str, Any]: + return { + "short_att": self.short_att, + "short_se": self.short_se, + "sigma2_control": self.sigma2_control, + "nu2_selection": self.nu2_selection, + "scale": self.scale, + "n_obs": self.n_obs, + "n_treated": self.n_treated, + "n_control": self.n_control, + "n_folds": self.n_folds, + "seed": self.seed, + "propensity_learner": repr(self.propensity_learner), + "outcome_learner": repr(self.outcome_learner), + } + + def to_dataframe(self) -> pd.DataFrame: + return pd.DataFrame([self.to_dict()]) + + def summary(self) -> str: + return "\n".join( + [ + "OVB Sensitivity Analysis for Difference-in-Differences", + f"Short ATT: {self.short_att:.6g}", + f"Short ATT SE: {self.short_se:.6g}", + f"Control trend var: {self.sigma2_control:.6g}", + f"Selection scale: {self.nu2_selection:.6g}", + f"OVB scale S0: {self.scale:.6g}", + f"Observations: {self.n_obs}", + ] + ) + + +class DIDOVBSensitivity: + """Cross-fitted canonical two-period DiD OVB sensitivity estimator.""" + + def __init__( + self, + *, + n_folds: int = 5, + seed: Optional[int] = None, + propensity_learner: Any = "logit", + outcome_learner: Any = "linear", + pscore_trim: float = 0.01, + ) -> None: + if isinstance(n_folds, bool) or int(n_folds) != n_folds or n_folds < 2: + raise ValueError(f"n_folds must be an integer >= 2, got {n_folds!r}") + if not 0.0 < float(pscore_trim) < 0.5: + raise ValueError("pscore_trim must be in (0, 0.5)") + self.n_folds = int(n_folds) + self.seed = seed + self.propensity_learner = propensity_learner + self.outcome_learner = outcome_learner + self.pscore_trim = float(pscore_trim) + + def fit( + self, + data: pd.DataFrame, + *, + outcome: str, + treatment: str, + time: str, + unit: str, + covariates: Sequence[str], + outcome_change: Optional[str] = None, + fold_ids: Optional[Sequence[int]] = None, + ) -> DIDOVBSensitivityResults: + """Fit on a balanced two-period panel.""" + if not isinstance(data, pd.DataFrame): + raise TypeError("data must be a pandas DataFrame") + covariates = list(covariates) + required = [outcome, treatment, time, unit, *covariates] + missing = [name for name in required if name not in data.columns] + if missing: + raise KeyError(f"missing columns: {missing}") + periods = list(pd.unique(data[time])) + if len(periods) != 2: + raise ValueError("DIDOVBSensitivity requires exactly two periods") + try: + periods = sorted(periods) + except TypeError as exc: + raise ValueError("time values must be sortable") from exc + + grouped = data.groupby(unit, sort=False, dropna=False) + if not bool((grouped.size() == 2).all()): + raise ValueError("each unit must have exactly one observation in each period") + pre = data.loc[data[time] == periods[0]].set_index(unit) + post = data.loc[data[time] == periods[1]].set_index(unit) + if not pre.index.is_unique or not post.index.is_unique: + raise ValueError("unit must identify one row per period") + common = pre.index.intersection(post.index) + if len(common) != len(pre) or len(common) != len(post): + raise ValueError("the two periods must contain the same units") + pre = pre.loc[common] + post = post.loc[common] + + d_pre = pd.to_numeric(pre[treatment], errors="coerce").to_numpy(dtype=float) + d_post = pd.to_numeric(post[treatment], errors="coerce").to_numpy(dtype=float) + if not np.all(np.isfinite(d_pre)) or not np.all(np.isfinite(d_post)): + raise ValueError("treatment must be finite and binary") + if not np.array_equal(d_pre, d_post) or not np.all(np.isin(d_pre, [0.0, 1.0])): + raise ValueError("treatment must be time-invariant and binary") + y_pre = pd.to_numeric(pre[outcome], errors="coerce").to_numpy(dtype=float) + y_post = pd.to_numeric(post[outcome], errors="coerce").to_numpy(dtype=float) + x = pre[covariates].apply(pd.to_numeric, errors="coerce").to_numpy(dtype=float) + if outcome_change is not None: + if outcome_change not in data.columns: + raise KeyError(f"missing columns: [{outcome_change!r}]") + dy_pre = pd.to_numeric(pre[outcome_change], errors="coerce").to_numpy(dtype=float) + dy_post = pd.to_numeric(post[outcome_change], errors="coerce").to_numpy(dtype=float) + delta_y = dy_post - dy_pre + else: + delta_y = y_post - y_pre + valid = np.isfinite(delta_y) & np.all(np.isfinite(x), axis=1) + if not np.all(valid): + d_pre, delta_y, x = d_pre[valid], delta_y[valid], x[valid] + if len(delta_y) < 4: + raise ValueError("at least four complete units are required") + n_treated = int(np.sum(d_pre == 1.0)) + n_control = int(np.sum(d_pre == 0.0)) + if n_treated < self.n_folds or n_control < self.n_folds: + raise ValueError("n_folds cannot exceed the number of treated or control units") + + if fold_ids is None: + rng = np.random.default_rng(self.seed) + folds = assign_folds(len(delta_y), self.n_folds, rng=rng, stratify=d_pre) + else: + provided_folds = np.asarray(fold_ids) + if provided_folds.shape != (len(delta_y),) or not np.issubdtype( + provided_folds.dtype, np.integer + ): + raise ValueError("fold_ids must be an integer vector with one entry per unit") + provided_folds = provided_folds.astype(np.int64, copy=False) + if np.any(provided_folds < 0) or np.any(provided_folds >= self.n_folds): + raise ValueError("fold_ids must lie in [0, n_folds)") + if len(np.unique(provided_folds)) != self.n_folds: + raise ValueError("fold_ids must contain every fold") + folds = FoldAssignment( + n_folds=self.n_folds, + n_units=len(delta_y), + fold_ids=provided_folds, + bitgen_state={}, + bitgen_name="provided", + stratify_labels=d_pre.copy(), + ) + ps_result = cross_fit_predict( + make_learner(self.propensity_learner, kind="classifier"), + x, + d_pre, + folds, + predict_method="predict_proba", + context_label="DIDOVBSensitivity propensity learner", + ) + outcome_result = cross_fit_predict( + make_learner(self.outcome_learner, kind="regressor"), + x, + delta_y, + folds, + predict_method="predict", + fit_mask=d_pre == 0.0, + context_label="DIDOVBSensitivity outcome learner", + ) + p = float(np.mean(d_pre)) + ps = np.clip(ps_result.oof_predictions, self.pscore_trim, 1.0 - self.pscore_trim) + omega = (ps / (1.0 - ps)) / (p / (1.0 - p)) + residual = delta_y - outcome_result.oof_predictions + score = (d_pre / p - (1.0 - d_pre) / (1.0 - p) * omega) * residual + short_att = float(np.mean(score)) + sigma2 = float(np.mean(residual[d_pre == 0.0] ** 2)) + nu2 = float(np.mean(omega[d_pre == 0.0] ** 2)) + if not np.isfinite(short_att) or not np.isfinite(sigma2) or not np.isfinite(nu2): + raise ValueError("OVB components are non-finite") + scale = float(np.sqrt(max(0.0, sigma2 * nu2))) + theta_if = score - short_att + sigma_if = (1.0 - d_pre) / (1.0 - p) * (residual**2 - sigma2) + nu_if = (1.0 - d_pre) / (1.0 - p) * (omega**2 - nu2) + scale_if = np.zeros_like(theta_if) + if scale > 0.0: + scale_if = (nu2 * sigma_if + sigma2 * nu_if) / (2.0 * scale) + short_se = float(np.sqrt(np.mean(theta_if**2) / len(theta_if))) + return DIDOVBSensitivityResults( + short_att=short_att, + short_se=short_se, + sigma2_control=sigma2, + nu2_selection=nu2, + scale=scale, + n_obs=len(delta_y), + n_treated=n_treated, + n_control=n_control, + n_folds=self.n_folds, + seed=self.seed, + propensity_learner=self.propensity_learner, + outcome_learner=self.outcome_learner, + _theta_if=theta_if, + _scale_if=scale_if, + ) diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index f7125471e..a081bd90f 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -82,6 +82,7 @@ The site is organized into 5 sections, each with a landing page: - [QDiD](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html): **Deprecated 3.9, removed 4.0 - use `ChangesInChanges(method="qdid")`.** Athey & Imbens (2006) quantile DiD comparison estimator (additive quantile-by-quantile DiD, matching R `qte::QDiD()` including its covariate branch via `covariates=`); same bootstrap machinery as ChangesInChanges. The paper recommends CiC over QDiD (scale-dependent model with testable restrictions; a non-monotonicity warning fires when violated - unconditional fits only, the covariate-path counterfactual quantile curve is monotone by construction). - [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html): Lee & Wooldridge (2025, 2026) rolling-transformation DiD — unit-specific demean/detrend converts panel to cross-section, supports staggered adoption with flexible control groups. Signature: `LWDiD(rolling='demean', estimation_method='reg', vcov_type='hc1', cluster=None, control_group='not_yet_treated', alpha=0.05, n_bootstrap=0, seed=None, pscore_trim=0.01, n_neighbors=1, caliper=None, with_replacement=True, n_jobs=1).fit(data, outcome, unit, time, treatment, first_treat=None, covariates=None)`. `estimation_method` values: `reg` (papers' RA), `ipw`, `dr` (papers' IPWRA, doubly robust), `psm`; `vcov_type` values: `classical`/`hc1`/`hc2`/`hc3` for `reg`; `ipw`/`dr` accept `hc1` only (influence-function variance); `psm` accepts `hc1` as configuration only - PSM inference is unavailable (NaN) pending an Abadie-Imbens matching variance; cluster-robust inference via the constructor's `cluster=` column (hc1/CR1 only, not a `vcov_type` value; rejected for `psm`). Per-period effects: post-fit `results.aggregate('event_study')`. - [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html): Chang (2020) double/debiased machine learning DiD — staggered ATT(g,t) with cross-fitted ML nuisances (DML2) and Neyman-orthogonal scores; covariates REQUIRED (conditional parallel trends). Signature: `DMLDiD(propensity_learner='logit', outcome_learner='linear', n_folds=5, control_group='never_treated', anticipation=0, alpha=0.05, n_bootstrap=0, bootstrap_weights=None, seed=None, base_period='varying', cband=True, pscore_trim=0.01, panel=True, cluster=None).fit(data, outcome, unit, time, first_treat, covariates, survey_design=None, bad_control=None, bad_control_covariates=None)`. `panel=False` = declared repeated cross sections (Chang Case 2: level outcomes, row-unique unit IDs, lambda-corrected variance). Survey/cluster support on BOTH designs (bad-control lane: panel only, `cluster=` only): `survey_design=` (pweight full-design TSL — weighted moments, PSU-cohesive folds, design-based variance with t-inference; a library extension of Chang's i.i.d. theory), replicate-weight designs (BRR/Fay/JK1/JKn/SDR — per-cell AND aggregate IF-reweighting variance, df = rank-1; cluster= and bootstrap combinations rejected), and coarser-than-unit `cluster=` (variance/folds only, kernels stay unweighted). Bad controls (Caetano, Callaway, Payne & Sant'Anna 2026): `fit(..., bad_control='x', bad_control_covariates=[outcome])` runs the paper's Neyman-orthogonal doubly-robust score (Eq. 10) - parallel trends conditional on the bad control's untreated path plus covariate unconfoundedness given its base-period value, W and Z - and stores the per-cell `ATT_X(g,t)` pre-test (`results.bad_control_summary()`); a bad control must NOT be in `covariates`. Learners: string names (`linear`/`ridge`/`sieve` regressors, `logit` classifier) or any object with fit/predict(_proba) (sklearn-compatible); `SieveLearner(k_max, criterion)` is exported for adaptive polynomial nuisances. Aggregation is POST-FIT: `results.aggregate('event_study'/'group'/'simple')`, plus `'total'` on panel non-survey fits (RCS and declared-survey fits fail 'total' closed); sup-t bands via bootstrap replay. With seed=None point estimates vary across fits (random folds); set seed for reproducibility. +- [DIDOVBSensitivity](https://diff-diff.readthedocs.io/en/stable/api/did_ovb.html): Wang, Sant'Anna, Chernozhukov & Cinelli (2026) canonical two-period DiD omitted-variable-bias sensitivity estimator. Signature: `DIDOVBSensitivity(n_folds=5, seed=None, propensity_learner='logit', outcome_learner='linear', pscore_trim=0.01).fit(data, outcome, treatment, time, unit, covariates, outcome_change=None, fold_ids=None)`. Requires a balanced two-period panel and time-invariant binary treatment; cross-fits the treatment propensity and untreated outcome change, reports the short ATT, control-trend scale, selection scale, restricted bias bounds via `bounds(trend_r2, selection_r2)`, and `robustness_value(null_value=0.0)` with RV/XRV. Use observed covariates or pre-trend benchmarks only when their identifying interpretation is substantively justified; this diagnostic does not repair a violated parallel-trends assumption. - [DurationDiD](https://diff-diff.readthedocs.io/en/stable/api/duration_did.html): Deaner & Ku (2026) causal duration DiD — two-group, common-timing design for a binary ABSORBING outcome (spell ended; once 1, always 1). Identifies off a restriction on the groups' UNTREATED HAZARDS, never outcome levels: `method='cd'` (constant additive hazard gap) or `method='ph'` (constant hazard ratio, mean-of-ratios estimator), fitted on the pre-treatment cumulative hazards (default equal weights over every eligible pre-treatment date; a window via `pre_periods=`/`pre_period_weights=`). Reports the absorption ATT at every post-treatment date and its uniform average as `att`, with the paper's whole-individual pooled bootstrap (centered pointwise intervals + simultaneous max-|t| band) and a fixed-anchor pre-treatment specification test in `results.pretest`. Signature: `DurationDiD(method='cd', n_bootstrap=1000, alpha=0.05, seed=None).fit(data, outcome, unit, time, treatment, last_pre_period=..., pre_periods=None, pre_period_weights=None)`; `last_pre_period` (the last untreated date) is REQUIRED and never inferred. Balanced equally spaced numeric grid, fixed 0/1 group indicator; no covariates, staggering, censoring, survey or cluster support. Every inference family is fully available or fully withheld (`inference_status`). Not admitted by DiagnosticReport/BusinessReport; `results.aggregate('event_study')` returns the unified container. - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html): Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings diff --git a/docs/api/did_ovb.rst b/docs/api/did_ovb.rst new file mode 100644 index 000000000..ff7fcbe4e --- /dev/null +++ b/docs/api/did_ovb.rst @@ -0,0 +1,68 @@ +DiD OVB Sensitivity +=================== + +``DIDOVBSensitivity`` implements the canonical two-period DiD omitted-variable +bias decomposition from Wang, Sant'Anna, Chernozhukov, and Cinelli. It is a +clean-room Python implementation for `diff-diff`; it is not a translation of +the GPL-3 ``dml.sensemakr`` source. + +The estimator accepts a balanced two-period panel, cross-fits a control trend +regression and treatment propensity, and reports the short ATT together with +the paper's estimable scale components. The result object then evaluates +restricted bias bounds and the RV/XRV sensitivity statistics. + +Basic usage +----------- + +.. code-block:: python + + from diff_diff import DIDOVBSensitivity + + result = DIDOVBSensitivity(n_folds=5, seed=42).fit( + data, + outcome="lemp", + treatment="treated", + time="year", + unit="county", + covariates=["lpop", "white", "pov"], + ) + + print(result.summary()) + bounds = result.bounds(trend_r2=1.0, selection_r2=0.5) + robustness = result.robustness_value(null_value=0.0) + contour = result.contour(null_value=0.0) + +``trend_r2`` is the residual trend variation explained by an omitted +confounder. ``selection_r2`` is the residual share of treatment-odds variation; +the corresponding selection factor is +``sqrt(selection_r2 / (1 - selection_r2))``. This parameterization follows the +paper's :math:`C_{0D}^2` decomposition and is intentionally different from a +generic regression ``R^2`` argument. + +R parity +-------- + +The independent parity oracle is generated with: + +.. code-block:: bash + + Rscript benchmarks/R/generate_did_ovb_parity.R \ + benchmarks/data/did_ovb_r_results.json \ + benchmarks/data/real/mpdta.csv + +The committed Python parity test compares the short ATT, scale components, +bounds, RV, and XRV. The fixture uses the 2006--2007 minimum-wage slice of the +public ``mpdta`` panel with the same fixed folds and learner specification in +both languages. + +Attribution and licensing +------------------------- + +The method is based on Wang et al., *Omitted Variable Bias in +Difference-in-Differences Designs*. The related R package ``dml.sensemakr`` is +GPL-3; its source is not included here. Users should cite the paper and the R +package when comparing results. + +.. automodule:: diff_diff.did_ovb + :members: + :undoc-members: diff --git a/docs/api/index.rst b/docs/api/index.rst index 5c02bb7ec..611dbd8be 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -36,6 +36,7 @@ regression discontinuity, and the Goodman-Bacon decomposition diagnostic: diff_diff.QDiD diff_diff.LWDiD diff_diff.DMLDiD + diff_diff.DIDOVBSensitivity diff_diff.DurationDiD diff_diff.BaconDecomposition diff_diff.StaggeredTripleDifference @@ -84,6 +85,9 @@ Result containers returned by estimators: diff_diff.changes_in_changes_results.ChangesInChangesResults diff_diff.lwdid_results.LWDiDResults diff_diff.dml_did_results.DMLDiDResults + diff_diff.DIDOVBSensitivityResults + diff_diff.DIDOVBBounds + diff_diff.DIDOVBRobustness diff_diff.duration_did_results.DurationDiDResults diff_diff.duration_did_results.DurationDiDPretestResults diff_diff.Comparison2x2 @@ -380,6 +384,7 @@ Estimators changes_in_changes lwdid dml_did + did_ovb duration_did bacon diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 16ae58353..d6d6c908e 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -1085,6 +1085,27 @@ sources: - path: docs/tutorials/33_bad_controls.ipynb type: tutorial + # ── DiD OVB sensitivity (did_ovb group) ──────────────────────────── + + diff_diff/did_ovb.py: + drift_risk: medium + docs: + - path: docs/api/did_ovb.rst + type: api_reference + - path: README.md + section: "Estimators (one-line catalog entry)" + type: user_guide + - path: docs/references.rst + type: user_guide + - path: diff_diff/guides/llms.txt + section: "Estimators" + type: user_guide + - path: docs/tutorials/34_did_ovb.ipynb + type: tutorial + - path: benchmarks/README.md + section: "DiD OVB sensitivity parity" + type: testing + # ── TROP (trop group) ────────────────────────────────────────────── diff_diff/trop.py: diff --git a/docs/index.rst b/docs/index.rst index 08f7380af..15394a0f8 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -187,6 +187,8 @@ Supported Estimators - Lee & Wooldridge (2025, 2026) rolling-transformation DiD; ``rolling='detrend'`` handles heterogeneous linear trends * - :class:`~diff_diff.DMLDiD` - Chang (2020) double/debiased ML DiD; staggered ATT(g,t) with cross-fitted nuisance learners (panel or declared repeated cross sections; survey/cluster support); Caetano et al. (2026) bad-control score via ``fit(bad_control=)`` (bad-control lane: panel only, ``cluster=`` only) + * - :class:`~diff_diff.DIDOVBSensitivity` + - Wang et al. (2026) two-period DiD omitted-variable-bias sensitivity analysis with cross-fitted nuisances, restricted bounds, RV, and XRV * - :class:`~diff_diff.DurationDiD` - Deaner & Ku (2026) causal duration DiD for binary absorbing outcomes; untreated-hazard restriction (common dynamics or proportional hazards), whole-individual bootstrap bands * - :class:`~diff_diff.QDiD` diff --git a/docs/references.rst b/docs/references.rst index e9f76999d..667805508 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -6,6 +6,8 @@ This library implements methods from the following scholarly works. Difference-in-Differences ------------------------- +- **Wang, J., Sant'Anna, P. H. C., Chernozhukov, V., & Cinelli, C. (2026).** "Omitted Variable Bias in Difference-in-Differences Designs." *arXiv preprint arXiv:2609.19386*. https://arxiv.org/abs/2609.19386 + - **Ashenfelter, O., & Card, D. (1985).** "Using the Longitudinal Structure of Earnings to Estimate the Effect of Training Programs." *The Review of Economics and Statistics*, 67(4), 648-660. https://doi.org/10.2307/1924810 - **Card, D., & Krueger, A. B. (1994).** "Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania." *The American Economic Review*, 84(4), 772-793. https://www.jstor.org/stable/2118030 diff --git a/docs/superpowers/plans/2026-09-28-did-ovb-sensitivity.md b/docs/superpowers/plans/2026-09-28-did-ovb-sensitivity.md new file mode 100644 index 000000000..2155559a0 --- /dev/null +++ b/docs/superpowers/plans/2026-09-28-did-ovb-sensitivity.md @@ -0,0 +1,161 @@ +# DiD OVB Sensitivity Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Implement a clean-room canonical DiD OVB sensitivity estimator in Python and verify its applied results against an independent R implementation. + +**Architecture:** Add a focused estimator and results module that consumes a balanced two-period panel, uses existing cross-fitting/learner helpers, and exposes bias bounds plus RV/XRV and contour data. Keep the current `DMLDiD` contract unchanged; parity is a separate R script and committed fixture. + +**Tech Stack:** Python, NumPy, pandas, SciPy, pytest, existing `diff_diff._crossfit` learners, Rscript for parity. + +--- + +### Task 1: Add the failing public API contract tests + +**Files:** +- Create: `tests/test_did_ovb_sensitivity.py` + +- [ ] **Step 1: Write tests for two-period input normalization and missing-input errors.** + +```python +def test_fit_accepts_balanced_two_period_panel_and_reports_short_att(): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit( + panel, outcome="y", treatment="treated", time="period", + unit="unit", covariates=["x"], + ) + assert np.isfinite(result.short_att) + assert result.n_obs == len(panel) // 2 + +def test_fit_rejects_more_than_two_periods(): + with pytest.raises(ValueError, match="exactly two periods"): + DIDOVBSensitivity().fit(panel_3_periods, outcome="y", treatment="treated", + time="period", unit="unit", covariates=["x"]) +``` + +- [ ] **Step 2: Run the focused tests and confirm they fail because the API is absent.** + +Run: `pytest -q tests/test_did_ovb_sensitivity.py` + +Expected: collection or import failure naming `DIDOVBSensitivity`. + +### Task 2: Implement input normalization and cross-fitted estimable components + +**Files:** +- Create: `diff_diff/did_ovb.py` +- Modify: `diff_diff/__init__.py` +- Test: `tests/test_did_ovb_sensitivity.py` + +- [ ] **Step 1: Add the smallest estimator skeleton and normalize a balanced panel.** + +Implement `fit()` with explicit validation for unit uniqueness, two periods, +finite outcome/covariates, binary treatment, and at least one treated and one +control unit. Construct `delta_y` as post minus pre and retain one row per unit. + +- [ ] **Step 2: Add failing tests for deterministic cross-fitting and component identities.** + +```python +def test_components_obey_scale_identity(): + result = DIDOVBSensitivity(n_folds=2, seed=11).fit(...) + assert result.scale == pytest.approx( + np.sqrt(result.sigma2_control * result.nu2_selection) + ) + assert result.short_att == pytest.approx(result.short_att_against_control) +``` + +- [ ] **Step 3: Implement cross-fitted nuisance predictions and orthogonal plug-ins.** + +Use `assign_folds` and `cross_fit_predict` with the existing learner factory; +estimate the control trend regression and treatment propensity, compute the +paper's short ATT, residual control variance, selection scale, and retained +diagnostics. Keep all arrays in the result only when needed for delta-method +variance and parity diagnostics. + +- [ ] **Step 4: Run focused tests and make them pass.** + +Run: `pytest -q tests/test_did_ovb_sensitivity.py -k 'normaliz or component'` + +### Task 3: Add bounds, inference, RV/XRV, and structured result output + +**Files:** +- Modify: `diff_diff/did_ovb.py` +- Modify: `diff_diff/__init__.py` +- Test: `tests/test_did_ovb_sensitivity.py` + +- [ ] **Step 1: Add failing tests for bounds and monotone robustness searches.** + +```python +def test_bounds_are_symmetric_and_expand_with_strength(): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit(...) + narrow = result.bounds(trend_strength=.2, selection_strength=.2) + wide = result.bounds(trend_strength=.5, selection_strength=.5) + assert narrow.lower < result.short_att < narrow.upper + assert wide.radius > narrow.radius + +def test_null_inside_original_interval_has_zero_rv_and_xrv(): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit(...) + out = result.robustness_value(null_value=result.short_att, alpha=.05) + assert out.rv == 0.0 + assert out.xrv == 0.0 +``` + +- [ ] **Step 2: Implement result dataclasses and normal/delta-method intervals.** + +Expose `bounds()`, `robustness_value()`, `contour()`, `summary()`, +`to_dict()`, and `to_dataframe()`. Validate strengths in `[0, 1]`, alpha in +`(0, .5)`, and return explicit NaNs/errors for non-identifiable fits. + +- [ ] **Step 3: Implement RV/XRV grid search against the paper's widest-CI definition.** + +Use a deterministic one-dimensional search over common strength for RV and +selection-only strength for XRV, retaining the grid and CI endpoints used to +make the decision. Do not silently replace a failed endpoint with a point +estimate. + +- [ ] **Step 4: Run focused tests and then the full relevant Python suite.** + +Run: `pytest -q tests/test_did_ovb_sensitivity.py tests/test_dml_did.py tests/test_methodology_dml_did.py` + +### Task 4: Add R application parity harness + +**Files:** +- Create: `benchmarks/R/generate_did_ovb_parity.R` +- Create: `benchmarks/data/did_ovb_parity_panel.csv` +- Create: `benchmarks/data/did_ovb_r_results.json` +- Create: `tests/test_did_ovb_r_parity.py` +- Modify: `benchmarks/R/requirements.R` +- Modify: `benchmarks/README.md` + +- [ ] **Step 1: Write the parity test before the R oracle exists.** + +The test must load the committed panel and R JSON, run the Python estimator with +the exact recorded settings, and compare short ATT, sigma/nu, scale, bounds, +RV, and XRV with named tolerances. Skip only when `Rscript` or the required R +package is unavailable; fail if the fixture is missing or malformed. + +- [ ] **Step 2: Implement the R script from public formulas and fixed settings.** + +The script writes JSON containing settings, sample counts, intermediate +components, and final sensitivity outputs. It must not import Python or call +the GPL-3 source through an undocumented bridge. + +- [ ] **Step 3: Generate the fixture and run Python/R parity.** + +Run: `Rscript benchmarks/R/generate_did_ovb_parity.R` + +Then: `pytest -q tests/test_did_ovb_r_parity.py` + +Expected: parity passes with all component-level comparisons inside the declared tolerance. + +### Task 5: Document, audit, and prepare the branch + +**Files:** +- Modify: `docs/api/dml_did.rst` or create `docs/api/did_ovb.rst` +- Modify: `docs/methodology/REGISTRY.md` +- Modify: `README.md` only if the public API is promoted there +- Modify: `LICENSE` notices only if reused MIT material is included + +- [ ] **Step 1: Document clean-room provenance, citations, limitations, and R parity command.** +- [ ] **Step 2: Run formatting/type checks and the complete targeted test set.** +- [ ] **Step 3: Inspect `git diff`, verify no GPL-3 source was copied, and confirm the working tree contains only intended files.** +- [ ] **Step 4: Commit locally with a focused message.** +- [ ] **Step 5: Do not push until the user explicitly asks for the reviewed commit/PR.** diff --git a/docs/superpowers/specs/2026-03-18-wooldridge-did-design.md b/docs/superpowers/specs/2026-03-18-wooldridge-did-design.md new file mode 100644 index 000000000..954771cf8 --- /dev/null +++ b/docs/superpowers/specs/2026-03-18-wooldridge-did-design.md @@ -0,0 +1,265 @@ +# WooldridgeDiD Estimator — Design Spec + +**Date:** 2026-03-18 +**Status:** Approved +**Scope:** Integrate Stata `jwdid` (Wooldridge ETWFE) functionality into diff-diff + +--- + +## 1. Background and Motivation + +The Stata package `jwdid` (Friosavila 2021) implements Wooldridge's (2021, 2023) Extended +Two-Way Fixed Effects (ETWFE) estimator for staggered DiD. Its key advantages over existing +diff-diff estimators are: + +- **Saturated regression**: estimates all cohort×time ATT(g,t) in a single pooled OLS, + more efficient than Callaway-Sant'Anna's pair-wise approach +- **Nonlinear extension**: Wooldridge (2023) extends ETWFE to logit and Poisson, avoiding + the incidental parameters problem — no other estimator in diff-diff supports this +- **Equivalence to CS**: under identical assumptions, ETWFE ATT(g,t) equals CS ATT(g,t) + +**Primary references:** +- Wooldridge (2021). "Two-Way Fixed Effects, the Two-Way Mundlak Regression, and + Difference-in-Differences Estimators." SSRN 3906345. +- Wooldridge (2023). "Simple approaches to nonlinear difference-in-differences with panel + data." *The Econometrics Journal*, 26(3), C31–C66. +- Friosavila (2021). `jwdid`: Stata module. SSC s459114. + +--- + +## 2. Architecture Overview + +### New files +| File | Purpose | +|------|---------| +| `diff_diff/wooldridge.py` | `WooldridgeDiD` estimator class | +| `diff_diff/wooldridge_results.py` | `WooldridgeDiDResults` dataclass | +| `tests/test_wooldridge.py` | Full test suite | + +### Modified files +| File | Change | +|------|--------| +| `diff_diff/__init__.py` | Export `WooldridgeDiD`, `WooldridgeDiDResults` | +| `docs/methodology/REGISTRY.md` | Add ETWFE methodology section | + +### Class hierarchy +`WooldridgeDiD` is a **standalone estimator** (same level as `CallawaySantAnna`, +`SunAbraham`, etc.), not inheriting from `DifferenceInDifferences`. It implements its own +`get_params` / `set_params`. + +--- + +## 3. Public API + +### Constructor + +```python +class WooldridgeDiD: + def __init__( + self, + method: str = "ols", # "ols" | "logit" | "poisson" + control_group: str = "not_yet_treated", # "never_treated" | "not_yet_treated" + anticipation: int = 0, # pre-treatment periods to include + demean_covariates: bool = True, # within cohort-period demeaning (jwdid default) + alpha: float = 0.05, + cluster: Optional[str] = None, # default: unit identifier (jwdid default) + n_bootstrap: int = 0, # >0 enables multiplier bootstrap + bootstrap_weights: str = "rademacher", # "rademacher" | "webb" | "mammen" + seed: Optional[int] = None, + rank_deficient_action: str = "warn", # "warn" | "error" | "silent" + ): ... +``` + +### fit() + +```python +def fit( + self, + data: pd.DataFrame, + outcome: str, + unit: str, + time: str, + cohort: str, # first treatment period; 0/NaN = never treated + exovar: Optional[List[str]] = None, # time-invariant covariates (no interaction) + xtvar: Optional[List[str]] = None, # time-varying covariates (demeaned within cohort-period) + xgvar: Optional[List[str]] = None, # cohort-interacted covariates +) -> "WooldridgeDiDResults": ... +``` + +**Notes:** +- `cohort` column convention: integer = first treatment period, 0 or NaN = never treated. + Consistent with `CallawaySantAnna`'s `cohort` parameter. +- Default clustering is at the `unit` level (matches `jwdid` default of `vce(cluster ivar)`). +- `demean_covariates=True` corresponds to `jwdid` default; `False` corresponds to `xasis` option. + +### get_params / set_params + +```python +def get_params(self) -> Dict[str, Any]: ... # returns all constructor params +def set_params(self, **params) -> "WooldridgeDiD": ... # sklearn-compatible +``` + +--- + +## 4. Results Object + +```python +@dataclass +class WooldridgeDiDResults: + # Raw cohort×time estimates — core output + group_time_effects: Dict[Tuple[Any, Any], Dict[str, Any]] + # key = (g, t); value = {"att", "se", "t_stat", "p_value", "conf_int"} + + # Simple aggregation (always computed on fit) + overall_att: float + overall_se: float + overall_t_stat: float + overall_p_value: float + overall_conf_int: Tuple[float, float] + + # Other aggregations (populated by .aggregate()) + group_effects: Optional[Dict[Any, Dict]] # keyed by cohort g + calendar_effects: Optional[Dict[Any, Dict]] # keyed by calendar period t + event_study_effects: Optional[Dict[int, Dict]] # keyed by relative period k = t - g + + # Metadata + method: str + control_group: str + groups: List[Any] + time_periods: List[Any] + n_obs: int + n_treated_units: int + n_control_units: int + alpha: float = 0.05 + + # Methods + def aggregate(self, type: str) -> "WooldridgeDiDResults": ... + # type: "simple" | "group" | "calendar" | "event" + # fills corresponding fields, returns self for chaining + + def summary(self, aggregation: str = "simple") -> str: ... + def to_dataframe(self, aggregation: str = "event") -> pd.DataFrame: ... + def plot_event_study(self, **kwargs) -> None: ... + def __repr__(self) -> str: ... +``` + +**Inference rule:** ALL inference fields (t_stat, p_value, conf_int) computed together +via `safe_inference()` from `diff_diff.utils`. Never computed individually. + +--- + +## 5. Internal Computation + +### 5a. Linear ETWFE (`method="ols"`) + +Faithful port of `jwdid` + `reghdfe`: + +1. **Filter observations**: keep control group (never- or not-yet-treated at time t) plus + all treated units. Drop observations where `t < g - anticipation`. + +2. **Build interaction matrix**: for each (g, t) with `t >= g - anticipation`, create + column `1(G_i = g) * 1(T = t)`. These are the β_{g,t} regressors. + +3. **Covariate preparation**: + - `exovar`: append as-is + - `xtvar`: demean within (cohort × period) cells when `demean_covariates=True` + - `xgvar`: interact with each cohort indicator + +4. **Absorb unit + time FE**: within-transformation (existing `absorb` mechanism in + `linalg.py`), not explicit dummies. + +5. **Solve**: `linalg.solve_ols()` → extract β_{g,t} coefficients and vcov submatrix. + +6. **Inference**: `linalg.compute_robust_vcov()` with unit-level clustering by default, + then `safe_inference()` for each (g, t) cell. + +7. **Bootstrap**: multiplier bootstrap supported for all inference; + wild cluster bootstrap supported for linear only (same as `DifferenceInDifferences`). + +### 5b. Nonlinear (`method="logit"|"poisson"`) + +Following Wooldridge (2023) pooled QMLE approach: + +- **Logit**: group-level FE (cohort × period), **not** individual FE — avoids incidental + parameters problem. Log-likelihood: Bernoulli QLL. +- **Poisson**: individual FE absorbed via PPML (iterative within-transformation). + Log-likelihood: Poisson QLL. + +Optimization: `scipy.optimize.minimize` (L-BFGS-B). Vcov from numerical Hessian +(`scipy.optimize.approx_fprime` second differences). + +**ATT computation via Average Structural Function (ASF):** +Coefficients on treatment interactions are not directly ATTs. Must compute: +``` +ATT(g,t) = mean[ g(X_i'β̂ + δ̂_{g,t}) - g(X_i'β̂) ] over treated units in (g,t) +``` +where `g(·)` = logistic or exp. Delta method for SE propagation. + +Bootstrap: multiplier bootstrap only (no wild cluster bootstrap for nonlinear). + +### 5c. Aggregation Weights (exact jwdid_estat formula) + +``` +ω(g,t) = number of unit-time observations in cell (g,t) + +simple: Σ_{g,t: t≥g} ω(g,t)·ATT(g,t) / Σ_{g,t: t≥g} ω(g,t) +group: Σ_{t≥g} ω(g,t)·ATT(g,t) / Σ_{t≥g} ω(g,t) ∀g +calendar: Σ_{g: t≥g} ω(g,t)·ATT(g,t) / Σ_{g: t≥g} ω(g,t) ∀t +event: Σ_g ω(g,g+k)·ATT(g,g+k) / Σ_g ω(g,g+k) ∀k +``` + +Aggregation SEs: delta method for linear (variance of weighted sum); bootstrap +distribution used when `n_bootstrap > 0`. + +--- + +## 6. Parallel Trends Assumptions + +| `control_group` | Assumption | Pre-treatment effects | +|-----------------|------------|----------------------| +| `"not_yet_treated"` (default) | Parallel trends between each cohort and not-yet-treated units | Constrained to zero by design | +| `"never_treated"` | Parallel trends between each cohort and never-treated units | Estimable (visible in event study k < 0) | + +--- + +## 7. Testing Strategy + +### test_wooldridge.py structure + +**API tests** +- Invalid `method` / `control_group` raises `ValueError` +- `get_params()` / `set_params()` round-trip +- Accessing `results_` before `fit()` raises + +**Basic functionality** +- Fit on `mpdta` dataset, all fields non-NaN (`assert_nan_inference()`) +- All four aggregations callable and produce sensible output +- `to_dataframe()` and `summary()` run without error + +**Methodology correctness** +- Linear ETWFE ATT(g,t) ≈ CallawaySantAnna ATT(g,t) on same data / same control group + (tolerance ~1e-3, both theoretically equivalent under OLS / same assumptions) +- Nonlinear: simulated binary data, logit ATT sign correct +- Aggregation weight verification: manual weighted average == `simple` ATT + +**Edge cases** +- `control_group="never_treated"` with pre-treatment k < 0 effects estimable +- `anticipation=1` shifts treatment window correctly +- All three covariate types passed simultaneously +- Single cohort degenerates to standard DiD + +**Slow tests** (`@pytest.mark.slow`) +- Bootstrap SE convergence (`ci_params.bootstrap(n, min_n=199)`, threshold 0.40/0.15) +- Nonlinear bootstrap + +--- + +## 8. Documentation + +- `docs/methodology/REGISTRY.md`: add "WooldridgeDiD / ETWFE" section with: + - Academic sources (Wooldridge 2021, 2023; Friosavila 2021) + - Estimator equation (saturated model) + - SE methods (unit-level cluster, multiplier bootstrap, wild cluster bootstrap for OLS) + - Edge cases: nonlinear ASF computation, covariate demeaning + - Note: `**Deviation from Stata:** nonlinear bootstrap uses multiplier (jwdid uses delta method)` +- Export as `WooldridgeDiD` and alias `ETWFE` in `__init__.py` diff --git a/docs/superpowers/specs/2026-09-28-did-ovb-sensitivity-design.md b/docs/superpowers/specs/2026-09-28-did-ovb-sensitivity-design.md new file mode 100644 index 000000000..244336509 --- /dev/null +++ b/docs/superpowers/specs/2026-09-28-did-ovb-sensitivity-design.md @@ -0,0 +1,101 @@ +# DiD OVB Sensitivity Design + +## Goal + +Add a clean-room Python implementation of the omitted-variable-bias sensitivity +framework for canonical two-period DiD, integrated with `diff-diff`'s existing +learner and result conventions and validated against an independent R oracle. + +## Scope + +The first release targets the paper's canonical DiD application rather than a +general replacement for the GPL-3 `dml.sensemakr` package. It will provide: + +- cross-fitted estimates of the short ATT, residual trend variance, and observed + selection scale; +- the OVB decomposition into scale, trend, selection, and alignment factors; +- bias bounds and delta-method confidence intervals under user-supplied + restrictions; +- point-estimate, confidence-interval, and contour-grid RV/XRV calculations; +- a structured results object with `summary()`, `to_dict()`, and + `to_dataframe()`; +- a deterministic R/Python parity fixture using the same canonical two-period + data and nuisance specifications. + +The implementation will not copy R source, function bodies, tests, or package +layout from `dml.sensemakr`. It will cite Wang et al. and the R package as an +independent reference implementation. The existing `DMLDiD` estimator remains +unchanged in the first slice; integration is through a new estimator/module so +existing behavior stays bit-stable. + +## Mathematical contract + +For observations with outcome change `delta_y`, treatment `d`, and covariates +`x`, the module estimates the paper's short parameter and scale components: + +\[ + \theta_s = E[\Delta Y-g_s(X)\mid D=1],\quad + \sigma^2_{0s}=E[(\Delta Y-g_s(X))^2\mid D=0],\quad + \nu^2_{0s}=E[(O_X/O)^2\mid D=0]. +\] + +The reported scale is `S0 = sqrt(sigma2_0s * nu2_0s)`. Given restrictions +`abs(rho) <= rho_max`, `C_delta_y <= trend_max`, and `C_d <= selection_max`, +the bias radius is `rho_max * trend_max * selection_max * S0` and the bound +estimates are `theta_s +/- radius`. Benchmarking against an observed covariate +and pre-trend extrapolation are separate, explicit result methods and are not +silently mixed into the baseline bound. + +RV/XRV follow the paper's definitions: they are the smallest common selection +and trend strength compatible with the requested null entering the widest +confidence interval; if the original confidence interval already contains the +null, both values are zero. The implementation exposes the search grid and the +resulting endpoint calculations so parity failures are diagnosable. + +## API and data flow + +The public surface will be: + +```python +from diff_diff import DIDOVBSensitivity + +result = DIDOVBSensitivity(n_folds=5, seed=42).fit( + data, outcome="y", treatment="treated", time="post", + covariates=["x1", "x2"], unit="unit", +) +result.bounds(trend_strength=1.0, selection_strength=0.25) +result.robustness_value(null_value=0.0, alpha=0.05) +result.contour(null_value=0.0, alpha=0.05) +``` + +The estimator accepts either a supplied `outcome_change` column or a balanced +two-period panel from which the change is constructed. Nuisance learners use +the repository's existing learner factory and cross-fitting helpers. No hidden +network download is performed by the estimator. + +## R parity + +The parity harness will use a committed, small canonical fixture and an R script +under `benchmarks/R/`. It will pin fold assignments, seed, clipping, learner +configuration, and null/alpha settings. The Python test compares the short ATT, +sigma/nu scale components, decomposition radius, bounds, and RV/XRV against +machine-readable R output with tolerances documented next to the fixture. +Tests requiring R will skip explicitly when `Rscript` or the required R package +is unavailable; they will never replace missing R output with Python-generated +goldens. + +## License and attribution + +The Python implementation is original MIT-licensed work in this repository. +The paper, `dml.sensemakr`, and `CS_RR` are cited in module and documentation +references. Any reused MIT-licensed `CS_RR` material or data receives its own +copyright/license notice. GPL-3 source from `dml.sensemakr` is not copied or +translated mechanically. + +## Non-goals for the first release + +- porting every general-purpose ATE/ATT feature in `dml.sensemakr`; +- reproducing R plotting objects or R-specific S3 classes; +- adding automatic data downloads; +- changing the existing `DMLDiD` result schema; +- pushing to `upstream` before local tests, R parity, and review are complete. diff --git a/docs/tutorials/34_did_ovb.ipynb b/docs/tutorials/34_did_ovb.ipynb new file mode 100644 index 000000000..68b0c9cb2 --- /dev/null +++ b/docs/tutorials/34_did_ovb.ipynb @@ -0,0 +1,132 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# DiD omitted-variable-bias sensitivity\n", + "\n", + "This tutorial applies `DIDOVBSensitivity` to a balanced two-period panel. The estimator follows Wang, Sant'Anna, Chernozhukov, and Cinelli (2026): it cross-fits the treatment propensity and the untreated outcome change, then reports restricted omitted-variable-bias bounds, the robustness value (RV), and the extreme robustness value (XRV).\n", + "\n", + "The sensitivity parameters describe a hypothetical omitted confounder. They do not prove that parallel trends holds and should be justified using the study context." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "from diff_diff import DIDOVBSensitivity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Create a balanced two-period panel\n", + "\n", + "The outcome is observed before and after treatment. The treatment is fixed at the unit level, while the covariates are measured before treatment." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.default_rng(20260929)\n", + "n = 400\n", + "unit = np.arange(n)\n", + "x1 = rng.normal(size=n)\n", + "x2 = rng.normal(size=n)\n", + "propensity = 1 / (1 + np.exp(-(0.35 * x1 - 0.25 * x2)))\n", + "treated = rng.binomial(1, propensity)\n", + "y_pre = 0.5 * x1 - 0.25 * x2 + rng.normal(scale=1.0, size=n)\n", + "y_post = y_pre + 0.75 + 0.40 * treated + 0.20 * x1 + rng.normal(scale=1.0, size=n)\n", + "\n", + "data = pd.DataFrame({\n", + " 'unit': np.tile(unit, 2),\n", + " 'time': np.repeat([0, 1], n),\n", + " 'treated': np.tile(treated, 2),\n", + " 'x1': np.tile(x1, 2),\n", + " 'x2': np.tile(x2, 2),\n", + " 'outcome': np.concatenate([y_pre, y_post]),\n", + "})\n", + "data.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Fit the sensitivity estimator\n", + "\n", + "Set a seed when you want the cross-fitting partition to be reproducible. The default learners are a logistic propensity model and a linear untreated-trend model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "result = DIDOVBSensitivity(n_folds=5, seed=17).fit(\n", + " data,\n", + " outcome='outcome',\n", + " treatment='treated',\n", + " time='time',\n", + " unit='unit',\n", + " covariates=['x1', 'x2'],\n", + ")\n", + "print(result.summary())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Translate sensitivity assumptions into ATT bounds\n", + "\n", + "`trend_r2` is the fraction of residual untreated-trend variation explained by the omitted confounder. `selection_r2` is the corresponding fraction of residual treatment-odds variation. The bound uses the paper's DiD OVB scale decomposition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bounds = result.bounds(trend_r2=0.25, selection_r2=0.25)\n", + "print(bounds)\n", + "\n", + "robustness = result.robustness_value(null_value=0.0, alpha=0.05)\n", + "print('RV:', robustness.rv)\n", + "print('XRV:', robustness.xrv)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The RV is the minimum equal strength of outcome-trend and treatment-selection confounding required to move the two-sided confidence interval to the null value. The XRV allows the omitted confounder to explain all remaining outcome-trend variation and therefore isolates the selection-strength requirement." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 5e6596ee3..0951fa2f8 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -260,6 +260,13 @@ Modern estimators for designs the basic toolkit cannot handle. Approach 1 via base-period covariates, the DMLDiD bad-control lane, and reading the ATT_X pre-test. + .. grid-item-card:: DiD OVB Sensitivity + :link: 34_did_ovb + :link-type: doc + + Quantify how strong an omitted confounder must be to overturn a + canonical two-period DiD estimate. + .. toctree:: :maxdepth: 1 @@ -279,6 +286,7 @@ Modern estimators for designs the basic toolkit cannot handle. LWDiD Rolling Transformations <31_lwdid> Double ML DiD (Chang 2020) <32_dml_did> Bad Controls (Caetano et al.) <33_bad_controls> + DiD OVB Sensitivity <34_did_ovb> Study Design ------------ diff --git a/tests/test_did_ovb_r_parity.py b/tests/test_did_ovb_r_parity.py new file mode 100644 index 000000000..8501bf56e --- /dev/null +++ b/tests/test_did_ovb_r_parity.py @@ -0,0 +1,46 @@ +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from diff_diff import DIDOVBSensitivity + + +@pytest.mark.skipif( + __import__("shutil").which("Rscript") is None, + reason="Rscript is required for the DiD OVB parity fixture", +) +def test_python_matches_independent_r_application_fixture(): + root = Path(__file__).parents[1] + fixture = json.loads((root / "benchmarks/data/did_ovb_r_results.json").read_text()) + data = pd.read_csv(root / "benchmarks/data/real/mpdta.csv") + data = data[data["year"].isin([2006, 2007]) & data["first.treat"].isin([0, 2007])].copy() + data["treated"] = (data["first.treat"] == 2007).astype(int) + data = data.sort_values(["countyreal", "year"], kind="stable").reset_index(drop=True) + fold_ids = np.arange(fixture["n_obs"], dtype=int) % fixture["settings"]["n_folds"] + + result = DIDOVBSensitivity(n_folds=2, seed=0).fit( + data, + outcome="lemp", + treatment="treated", + time="year", + unit="countyreal", + covariates=["lpop"], + fold_ids=fold_ids, + ) + expected = fixture + + for name in ("short_att", "short_se", "sigma2_control", "nu2_selection", "scale"): + assert getattr(result, name) == pytest.approx(expected[name], rel=5e-5, abs=5e-7) + + bounds = result.bounds(trend_r2=1.0, selection_r2=0.5) + for name in ("lower", "upper", "radius", "lower_se", "upper_se", "lower_ci", "upper_ci"): + assert getattr(bounds, name) == pytest.approx(expected["bounds"][name], rel=5e-5, abs=5e-7) + + robustness = result.robustness_value( + null_value=expected["settings"]["null_value"], alpha=expected["settings"]["alpha"] + ) + assert robustness.rv == pytest.approx(expected["robustness"]["rv"], rel=5e-5, abs=5e-7) + assert robustness.xrv == pytest.approx(expected["robustness"]["xrv"], rel=5e-5, abs=5e-7) diff --git a/tests/test_did_ovb_sensitivity.py b/tests/test_did_ovb_sensitivity.py new file mode 100644 index 000000000..a527ffad3 --- /dev/null +++ b/tests/test_did_ovb_sensitivity.py @@ -0,0 +1,131 @@ +import numpy as np +import pandas as pd +import pytest + +from diff_diff import DIDOVBSensitivity + + +@pytest.fixture +def panel(): + rows = [] + rng = np.random.default_rng(123) + for unit in range(80): + treated = int(unit < 40) + x = (unit - 39.5) / 40.0 + baseline = 0.4 * x + rng.normal(scale=0.1) + trend = 0.25 * x + rng.normal(scale=0.08) + effect = 0.8 if treated else 0.0 + rows.extend( + [ + {"unit": unit, "period": 0, "treated": treated, "x": x, "y": baseline}, + { + "unit": unit, + "period": 1, + "treated": treated, + "x": x, + "y": baseline + trend + effect, + }, + ] + ) + return pd.DataFrame(rows) + + +def test_fit_accepts_balanced_two_period_panel_and_reports_short_att(panel): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit( + panel, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + + assert np.isfinite(result.short_att) + assert result.n_obs == len(panel) // 2 + assert result.short_se > 0 + + +def test_fit_rejects_more_than_two_periods(panel): + three_periods = pd.concat( + [panel, panel.assign(period=2, y=panel["y"] + 0.1)], ignore_index=True + ) + with pytest.raises(ValueError, match="exactly two periods"): + DIDOVBSensitivity().fit( + three_periods, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + + +def test_components_obey_scale_identity(panel): + result = DIDOVBSensitivity(n_folds=2, seed=11).fit( + panel, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + + assert result.scale == pytest.approx(np.sqrt(result.sigma2_control * result.nu2_selection)) + assert result.sigma2_control > 0 + assert result.nu2_selection >= 1 + + +def test_bounds_are_symmetric_and_expand_with_strength(panel): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit( + panel, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + narrow = result.bounds(trend_r2=0.2, selection_r2=0.2) + wide = result.bounds(trend_r2=0.5, selection_r2=0.5) + + assert narrow.lower < result.short_att < narrow.upper + assert narrow.lower == pytest.approx(result.short_att - narrow.radius) + assert narrow.upper == pytest.approx(result.short_att + narrow.radius) + assert wide.radius > narrow.radius + + +def test_null_inside_original_interval_has_zero_rv_and_xrv(panel): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit( + panel, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + robustness = result.robustness_value(null_value=result.short_att, alpha=0.05) + + assert robustness.rv == 0.0 + assert robustness.xrv == 0.0 + + +def test_result_exports_and_contour_are_structured(panel): + result = DIDOVBSensitivity(n_folds=2, seed=7).fit( + panel, + outcome="y", + treatment="treated", + time="period", + unit="unit", + covariates=["x"], + ) + + contour = result.contour(null_value=0.0, n_grid=11) + assert list(contour.columns) == [ + "strength", + "rv_lower", + "rv_upper", + "xrv_lower", + "xrv_upper", + ] + assert len(contour) == 11 + assert result.to_dict()["short_att"] == pytest.approx(result.short_att) + assert "OVB Sensitivity" in result.summary()