|
| 1 | +#!/usr/bin/env Rscript |
| 2 | +# Generate R `pretrends` parity goldens for diff-diff PreTrendsPower (PR-C). |
| 3 | +# |
| 4 | +# This script is committed in PR-B (PreTrendsPower implementation audit, |
| 5 | +# Roth 2022); the JSON goldens at ../data/r_pretrends_golden.json are |
| 6 | +# DEFERRED to PR-C. Running this script writes the JSON to that path; PR-C |
| 7 | +# pins the R `pretrends` package commit / release, runs this script, and |
| 8 | +# commits the resulting JSON to land the parity tests. |
| 9 | +# |
| 10 | +# Requires: |
| 11 | +# - R 4.4+ (tested on 4.5.2) |
| 12 | +# - install.packages("remotes") |
| 13 | +# - remotes::install_github("jonathandroth/pretrends", ref = "<PR-C-PIN>") |
| 14 | +# - install.packages("jsonlite") |
| 15 | +# |
| 16 | +# **R `pretrends` commit pin (TODO — PR-C):** the audited revision MUST be |
| 17 | +# recorded here before parity assertions are committed. As of 2026-05-18 |
| 18 | +# (PR-B implementation date) the script targets the default `main` branch |
| 19 | +# at https://github.com/jonathandroth/pretrends with no pin. PR-C will |
| 20 | +# replace `<PR-C-PIN>` with the exact commit hash AND verify the surface |
| 21 | +# claims documented in REGISTRY.md `## PreTrendsPower` and the paper |
| 22 | +# review's "R `pretrends` package version pin (provisional)" Gaps bullet. |
| 23 | +# |
| 24 | +# Output: ../data/r_pretrends_golden.json |
| 25 | +# |
| 26 | +# diff-diff PreTrendsPower with `pretest_form='nis'` (the new default per |
| 27 | +# PR-B Step 2) is expected to match the values in this JSON at atol=1e-6 |
| 28 | +# along a three-tier contract: |
| 29 | +# (1) NIS box probability `P(β̂_pre ∈ B_NIS(Σ))` at fixed M values on |
| 30 | +# all 3 fixtures; |
| 31 | +# (2) MDV / gamma_p (slope at target power 0.5 and 0.8) on regular and |
| 32 | +# irregular pre-period grids; |
| 33 | +# (3) γ-unit MDV invariance: PR-B's "skip L2 norm for linear with |
| 34 | +# relative_times" path produces MDV in Roth's γ units exactly, |
| 35 | +# matching R's `slope_for_power()` which also reports γ. |
| 36 | +# |
| 37 | +# Three fixtures (matched to test_methodology_pretrends.py expectations): |
| 38 | +# 1. uniform_3_pre_periods_no_anticipation — K=3 regular grid (t ∈ {-3, -2, -1}), |
| 39 | +# never-treated control. Default-case parity baseline. |
| 40 | +# 2. irregular_pre_periods — K=3 with relative_times = [-5, -3, -1]. |
| 41 | +# Exercises the PR-B γ-unit linear-pattern fix. |
| 42 | +# 3. anticipation_shifted — K=4 with anticipation=1 (pre-cutoff at t<-1, |
| 43 | +# so pre-periods are {-5, -4, -3, -2}). Verifies the pre-period filter |
| 44 | +# logic in `_extract_pre_period_params`. |
| 45 | +# |
| 46 | +# Run: |
| 47 | +# cd benchmarks/R && Rscript generate_pretrends_golden.R |
| 48 | + |
| 49 | +suppressPackageStartupMessages({ |
| 50 | + library(pretrends) |
| 51 | + library(jsonlite) |
| 52 | +}) |
| 53 | + |
| 54 | +stopifnot(packageVersion("pretrends") >= "0.1.0") |
| 55 | + |
| 56 | +# --------------------------------------------------------------------------- |
| 57 | +# DGP helper: build a synthetic event-study coefficient vector + VCV under a |
| 58 | +# stylized null DGP (β = 0, Σ_22 ~ correlated). Mirrors the simulation |
| 59 | +# fixtures in test_methodology_pretrends.py. |
| 60 | +# --------------------------------------------------------------------------- |
| 61 | + |
| 62 | +build_event_study_fixture <- function( |
| 63 | + pre_periods, |
| 64 | + post_periods, |
| 65 | + sigma2 = 0.04, |
| 66 | + rho = 0.3, |
| 67 | + seed = 42L |
| 68 | +) { |
| 69 | + # Generate a correlated equicorrelation Σ across all (pre + post) periods. |
| 70 | + # Realized β̂ drawn from N(0, Σ) — null DGP, no real treatment effect. |
| 71 | + set.seed(seed) |
| 72 | + all_periods <- c(pre_periods, post_periods) |
| 73 | + K_total <- length(all_periods) |
| 74 | + Sigma <- sigma2 * (rho * matrix(1, K_total, K_total) + (1 - rho) * diag(K_total)) |
| 75 | + beta_hat <- MASS::mvrnorm(1, mu = rep(0, K_total), Sigma = Sigma) |
| 76 | + |
| 77 | + list( |
| 78 | + beta_hat = beta_hat, |
| 79 | + Sigma = Sigma, |
| 80 | + all_periods = all_periods, |
| 81 | + pre_periods = pre_periods, |
| 82 | + post_periods = post_periods |
| 83 | + ) |
| 84 | +} |
| 85 | + |
| 86 | +# --------------------------------------------------------------------------- |
| 87 | +# Extract R pretrends() output into a fixture-shaped list. |
| 88 | +# --------------------------------------------------------------------------- |
| 89 | + |
| 90 | +extract_pretrends <- function(fixture_data, fixture_name) { |
| 91 | + beta_hat <- fixture_data$beta_hat |
| 92 | + Sigma <- fixture_data$Sigma |
| 93 | + pre_periods <- fixture_data$pre_periods |
| 94 | + post_periods <- fixture_data$post_periods |
| 95 | + all_periods <- fixture_data$all_periods |
| 96 | + |
| 97 | + # R `pretrends` expects: betahat (coefficient vector), sigma (VCV matrix), |
| 98 | + # tVec (relative-time labels including the reference period 0, omitted |
| 99 | + # from betahat / sigma per convention), referencePeriod = 0, alpha = 0.05. |
| 100 | + |
| 101 | + # The `slopes_for_power` helper returns gamma values at target power. |
| 102 | + # For the three-tier parity contract, we capture both NIS power at a fixed |
| 103 | + # slope and the inverse (γ_p MDV) at target power 0.5 and 0.8. |
| 104 | + |
| 105 | + # NIS power at fixed gamma values (for tier-1 parity): |
| 106 | + gamma_test_values <- c(0.0, 0.2, 0.5, 1.0) |
| 107 | + power_values <- sapply(gamma_test_values, function(g) { |
| 108 | + # Build δ = γ * |t| for pre-periods (Roth's δ_t = γ·t convention, |
| 109 | + # using |t| since pre-period t < 0). |
| 110 | + delta_pre <- g * abs(pre_periods) |
| 111 | + # `pretrends` package: pretrends() with explicit delta vector. |
| 112 | + # The exact R API: pretrends(betahat, sigma, tVec, referencePeriod, |
| 113 | + # deltahypothesis, ...). |
| 114 | + # PR-C: replace this stub with the actual R pretrends() call and |
| 115 | + # extract the rejection probability. |
| 116 | + NA_real_ # PR-C will populate |
| 117 | + }) |
| 118 | + |
| 119 | + # γ_p MDV: solve for γ such that NIS rejection probability = target power. |
| 120 | + # R `slope_for_power(betahat, sigma, tVec, referencePeriod, power)`. |
| 121 | + gamma_p_values <- sapply(c(0.5, 0.8), function(p) { |
| 122 | + # PR-C: replace with actual R slope_for_power() call. |
| 123 | + NA_real_ |
| 124 | + }) |
| 125 | + |
| 126 | + list( |
| 127 | + panel = list( |
| 128 | + pre_periods = as.integer(pre_periods), |
| 129 | + post_periods = as.integer(post_periods), |
| 130 | + all_periods = as.integer(all_periods), |
| 131 | + beta_hat = as.numeric(beta_hat), |
| 132 | + Sigma = Sigma |
| 133 | + ), |
| 134 | + r_power_at_gamma = list( |
| 135 | + gamma_test_values = as.numeric(gamma_test_values), |
| 136 | + power_values = as.numeric(power_values) |
| 137 | + ), |
| 138 | + r_gamma_p = list( |
| 139 | + target_power = c(0.5, 0.8), |
| 140 | + gamma_p_values = as.numeric(gamma_p_values) |
| 141 | + ), |
| 142 | + fixture_name = fixture_name |
| 143 | + ) |
| 144 | +} |
| 145 | + |
| 146 | +# --------------------------------------------------------------------------- |
| 147 | +# Fixtures |
| 148 | +# --------------------------------------------------------------------------- |
| 149 | + |
| 150 | +cat("Building fixture 1: uniform_3_pre_periods_no_anticipation...\n") |
| 151 | +f1 <- build_event_study_fixture( |
| 152 | + pre_periods = c(-3L, -2L, -1L), |
| 153 | + post_periods = c(1L, 2L, 3L), |
| 154 | + seed = 101L |
| 155 | +) |
| 156 | +fixture_1 <- extract_pretrends(f1, "uniform_3_pre_periods_no_anticipation") |
| 157 | + |
| 158 | +cat("Building fixture 2: irregular_pre_periods...\n") |
| 159 | +# K=3 with t ∈ {-5, -3, -1}. Tests PR-B's γ-unit linear-pattern fix: |
| 160 | +# pre-PR-B Python with normalized count-based weights would silently report |
| 161 | +# MDV in [0.45, 0.30, 0.15] / sqrt(0.3) units, not γ. R `slope_for_power()` |
| 162 | +# always reports γ; Python's PR-B Step 4 makes the two match at atol=1e-6. |
| 163 | +f2 <- build_event_study_fixture( |
| 164 | + pre_periods = c(-5L, -3L, -1L), |
| 165 | + post_periods = c(1L, 2L, 3L), |
| 166 | + seed = 202L |
| 167 | +) |
| 168 | +fixture_2 <- extract_pretrends(f2, "irregular_pre_periods") |
| 169 | + |
| 170 | +cat("Building fixture 3: anticipation_shifted...\n") |
| 171 | +# K=4 pre-periods with anticipation=1. Real pre-treatment cutoff is t < -1, |
| 172 | +# so the {-5, -4, -3, -2} cells are the genuine pre-periods; t=-1 is the |
| 173 | +# anticipation window. Tests the pre-period filtering logic. |
| 174 | +f3 <- build_event_study_fixture( |
| 175 | + pre_periods = c(-5L, -4L, -3L, -2L), # genuine pre-periods (cutoff = -1) |
| 176 | + post_periods = c(1L, 2L, 3L), |
| 177 | + seed = 303L |
| 178 | +) |
| 179 | +fixture_3 <- extract_pretrends(f3, "anticipation_shifted") |
| 180 | + |
| 181 | +# --------------------------------------------------------------------------- |
| 182 | +# Write JSON |
| 183 | +# --------------------------------------------------------------------------- |
| 184 | + |
| 185 | +out <- list( |
| 186 | + meta = list( |
| 187 | + generated_at = format(Sys.Date()), |
| 188 | + pretrends_version = as.character(packageVersion("pretrends")), |
| 189 | + pretrends_commit = "<PR-C-PIN>", # TODO PR-C: replace with actual git SHA |
| 190 | + r_version = R.version.string, |
| 191 | + description = paste( |
| 192 | + "Roth (2022) PreTrendsPower parity goldens for diff-diff", |
| 193 | + "compute_pretrends_power / PreTrendsPower (PR-C parity target).", |
| 194 | + "Parity at atol=1e-6 along a three-tier contract:", |
| 195 | + "(1) NIS box probability at fixed γ values on all 3 fixtures;", |
| 196 | + "(2) γ_p MDV (slope at target power 0.5 and 0.8) on regular and", |
| 197 | + "irregular grids;", |
| 198 | + "(3) γ-unit MDV invariance: PR-B's skip-L2-norm path produces MDV", |
| 199 | + "in Roth's γ units exactly, matching R's slope_for_power().", |
| 200 | + "See diff-diff/docs/methodology/papers/roth-2022-review.md for", |
| 201 | + "the full derivation." |
| 202 | + ) |
| 203 | + ), |
| 204 | + uniform_3_pre_periods_no_anticipation = fixture_1, |
| 205 | + irregular_pre_periods = fixture_2, |
| 206 | + anticipation_shifted = fixture_3 |
| 207 | +) |
| 208 | + |
| 209 | +out_path <- "../data/r_pretrends_golden.json" |
| 210 | +write_json(out, out_path, pretty = TRUE, digits = NA, auto_unbox = TRUE) |
| 211 | +cat(sprintf("Wrote %s\n", out_path)) |
| 212 | +cat("\n") |
| 213 | +cat("PR-C TODO checklist:\n") |
| 214 | +cat(" [ ] Replace <PR-C-PIN> commit-hash placeholder above with actual\n") |
| 215 | +cat(" git SHA from https://github.com/jonathandroth/pretrends.\n") |
| 216 | +cat(" [ ] Replace the NA_real_ stubs in extract_pretrends() with the\n") |
| 217 | +cat(" actual pretrends::pretrends() / slope_for_power() calls.\n") |
| 218 | +cat(" [ ] Verify the surface claims in REGISTRY.md PreTrendsPower\n") |
| 219 | +cat(" Reference implementations section against the pinned revision.\n") |
| 220 | +cat(" [ ] Activate tests/test_methodology_pretrends.py::TestPretrendsParityR\n") |
| 221 | +cat(" (currently skips via @pytest.mark.skipif when the JSON is missing).\n") |
| 222 | +cat(" [ ] Flip METHODOLOGY_REVIEW.md PreTrendsPower row from\n") |
| 223 | +cat(" **Complete** (R parity pending) → **Complete**.\n") |
0 commit comments