|
| 1 | +#!/usr/bin/env Rscript |
| 2 | +# Golden generator: wild_bootstrap_se vs R `fwildclusterboot::boottest` |
| 3 | +# (Roodman, MacKinnon, Nielsen & Webb 2019; Cameron-Gelbach-Miller 2008). |
| 4 | +# |
| 5 | +# Writes a fixed few-cluster 2x2 DiD design and the R reference so that |
| 6 | +# tests/test_wild_bootstrap.py::TestWildBootstrapParityR can pin Python output |
| 7 | +# against R without requiring R at test time. |
| 8 | +# |
| 9 | +# Outputs (checked into the repo): |
| 10 | +# benchmarks/data/wild_cluster_boot_test_data.csv (cluster, treated, post, y) |
| 11 | +# benchmarks/data/wild_cluster_boot_golden.json |
| 12 | +# |
| 13 | +# Usage: |
| 14 | +# Rscript benchmarks/R/generate_wild_cluster_boot_golden.R |
| 15 | +# |
| 16 | +# Notes: |
| 17 | +# - boottest defaults are the WCR (restricted) bootstrap: impose_null=TRUE, |
| 18 | +# type="rademacher". boottest switches to full enumeration once B reaches the |
| 19 | +# number of possible draws (2^G); with G=6 and B=99999 it enumerates all |
| 20 | +# 2^6=64 raw sign-vectors (deterministic, no RNG), of which 2^(G-1)=32 have |
| 21 | +# distinct |t|. The library's wild_bootstrap_se uses the same trigger |
| 22 | +# (2^n_clusters <= n_bootstrap) and also materializes all 2^G raw vectors, |
| 23 | +# reporting n_bootstrap = 2^G; both yield the same p-value and CI. |
| 24 | +# - The DGP uses a deliberately weak effect + heavy noise so the bootstrap |
| 25 | +# p-value is INTERIOR (not 0/1), letting the test pin the exact p-value. |
| 26 | +# boottest counts strict exceedances |t*| > |t0|; the library matches this |
| 27 | +# (it floors the reported p at 1/(B+1), inactive for an interior p). |
| 28 | +# - se is the analytical CR1 cluster-robust SE = (G/(G-1))((N-1)/(N-k)) sandwich, |
| 29 | +# which the library reports as `se` and uses to studentize the test. boottest |
| 30 | +# studentizes with the same CR1 SE, so teststat == coef/se. |
| 31 | +# - The CI is obtained by inverting the bootstrap test (boottest's grid search |
| 32 | +# vs the library's bisection); they agree to ~1e-4 on this design. |
| 33 | + |
| 34 | +suppressMessages({ |
| 35 | + library(fwildclusterboot) |
| 36 | + library(fixest) |
| 37 | + library(jsonlite) |
| 38 | +}) |
| 39 | + |
| 40 | +set.seed(20240624) |
| 41 | +G <- 6 |
| 42 | +obs_per_cluster <- 8 |
| 43 | +rows <- list() |
| 44 | +i <- 1 |
| 45 | +for (c in 0:(G - 1)) { |
| 46 | + is_treated <- as.integer(c < G / 2) |
| 47 | + cluster_effect <- rnorm(1, 0, 1.5) |
| 48 | + for (o in 1:obs_per_cluster) { |
| 49 | + for (period in c(0, 1)) { |
| 50 | + y <- 4.0 + cluster_effect + 1.0 * period |
| 51 | + if (is_treated == 1 && period == 1) y <- y + 0.3 # weak effect |
| 52 | + y <- y + rnorm(1, 0, 4.0) # heavy noise -> interior p |
| 53 | + rows[[i]] <- data.frame(cluster = c, treated = is_treated, post = period, y = y) |
| 54 | + i <- i + 1 |
| 55 | + } |
| 56 | + } |
| 57 | +} |
| 58 | +df <- do.call(rbind, rows) |
| 59 | +df$inter <- df$treated * df$post |
| 60 | + |
| 61 | +data_path <- "benchmarks/data/wild_cluster_boot_test_data.csv" |
| 62 | +write.csv(df[, c("cluster", "treated", "post", "y")], data_path, row.names = FALSE) |
| 63 | + |
| 64 | +m <- feols(y ~ treated + post + inter, data = df, cluster = ~cluster) |
| 65 | +coef_est <- as.numeric(coef(m)["inter"]) |
| 66 | +se_cr1 <- as.numeric(se(m)["inter"]) # CR1 clustered SE |
| 67 | + |
| 68 | +run_bt <- function(ptype) { |
| 69 | + set.seed(1) |
| 70 | + bt <- boottest(m, param = "inter", clustid = ~cluster, B = 99999, |
| 71 | + type = "rademacher", impose_null = TRUE, |
| 72 | + p_val_type = ptype, conf_int = TRUE, sign_level = 0.05) |
| 73 | + list(p_val = as.numeric(bt$p_val), |
| 74 | + t_stat = as.numeric(bt$t_stat), |
| 75 | + conf_int = as.numeric(bt$conf_int), |
| 76 | + n_clusters = as.integer(bt$N_G)) |
| 77 | +} |
| 78 | + |
| 79 | +golden <- list( |
| 80 | + n_clusters = G, |
| 81 | + coef = coef_est, |
| 82 | + se_cr1 = se_cr1, |
| 83 | + two_tailed = run_bt("two-tailed"), |
| 84 | + equal_tailed = run_bt("equal-tailed") |
| 85 | +) |
| 86 | + |
| 87 | +json_path <- "benchmarks/data/wild_cluster_boot_golden.json" |
| 88 | +write(toJSON(golden, auto_unbox = TRUE, digits = 12, pretty = TRUE), json_path) |
| 89 | +cat("Wrote", data_path, "and", json_path, "\n") |
| 90 | +cat(sprintf("coef=%.6f se=%.6f two-tailed p=%.6f CI=[%.6f, %.6f]\n", |
| 91 | + coef_est, se_cr1, golden$two_tailed$p_val, |
| 92 | + golden$two_tailed$conf_int[1], golden$two_tailed$conf_int[2])) |
0 commit comments