generate_dcdh_dynr_test_values.R

#!/usr/bin/env Rscript
# Generate golden values for de Chaisemartin-D'Haultfoeuille (dCDH)
# parity tests at horizon l=1.
#
# This script fits the R `DIDmultiplegtDYN` package's `did_multiplegt_dyn`
# function (the official R implementation of the dCDH dynamic-effects
# companion paper, NBER WP 29873) on a set of canonical reversible-treatment
# scenarios. At horizon l=1, did_multiplegt_dyn computes DID_1, which is
# numerically identical to DID_M of the AER 2020 paper. Phase 1 of the
# Python diff-diff dCDH implementation tests for parity against these
# golden values.
#
# Usage:
#   Rscript benchmarks/R/generate_dcdh_dynr_test_values.R
#
# Prerequisites:
#   install.packages("DIDmultiplegtDYN")  # CRAN v2.3.3+
#   install.packages("jsonlite")
#
# Output:
#   benchmarks/data/dcdh_dynr_golden_values.json
#
# Each scenario exports:
#   - data: the simulated dataset (so Python tests use identical data)
#   - params: the dCDH options used
#   - results: DID_1 point estimate, SE, CI, placebo, switcher counts

library(DIDmultiplegtDYN)
library(jsonlite)
suppressMessages(library(polars))  # required by DIDmultiplegtDYN >= 2.x

# Pin DIDmultiplegtDYN to an exact version because the `by_path` output
# slots (res$by_levels, res$by_level_i) were introduced in v2.3.3 and
# the structure is not version-stable per the R package's own docs. A
# floor constraint (`>= 2.3.3`) could silently drift the fixture schema
# when regenerated against a future release. Update this pin *and* re-
# run TestDCDHDynRParityByPath when bumping to a newer known-compatible
# release; extend to an explicit allowlist (e.g. `%in% c("2.3.3",
# "2.3.4")`) once a second version is verified.
stopifnot(packageVersion("DIDmultiplegtDYN") == "2.3.3")

cat("Generating dCDH golden values via DIDmultiplegtDYN at l=1...\n")

output_path <- file.path("benchmarks", "data", "dcdh_dynr_golden_values.json")

# ---------------------------------------------------------------------------
# Helper: Python-mirror reversible-treatment generator.
# Mirrors generate_reversible_did_data() in diff_diff/prep_dgp.py at the
# STRUCTURAL level — the two implementations apply the same pattern logic
# (single_switch / joiners_only / leavers_only / mixed_single_switch) and
# the same fixed-effect / treatment-effect / time-trend / noise model. They
# do NOT produce bit-identical draws even with the same seed: R's set.seed
# and NumPy's default_rng use different RNGs and the parity tests don't
# rely on RNG identity. Instead, the parity tests load THIS R script's
# golden-value JSON output and pass the SAME data (group/period/treatment/
# outcome columns) to the Python estimator, so both sides operate on
# byte-identical input regardless of how it was originally generated.
# ---------------------------------------------------------------------------
gen_reversible <- function(n_groups, n_periods, pattern, seed,
                           p_switch = 0.2, initial_treat_frac = 0.3,
                           cycle_length = 2, treatment_effect = 2.0,
                           heterogeneous_effects = FALSE, effect_sd = 0.5,
                           group_fe_sd = 2.0, time_trend = 0.1, noise_sd = 0.5,
                           n_never_treated = 20, n_always_treated = 20,
                           L_max = 3) {
  # n_never_treated and n_always_treated add stable control cohorts so
  # both Python (AER 2020 zero-retention) and R DIDmultiplegtDYN (dynamic
  # paper, drop-cohort) implementations have controls available at every
  # period — eliminating the methodology divergence and giving a clean
  # parity comparison. The total returned panel has
  # n_groups + n_never_treated + n_always_treated groups.
  set.seed(seed)

  # --- Build the (n_groups, n_periods) treatment matrix ---
  D <- matrix(0L, nrow = n_groups, ncol = n_periods)

  if (pattern == "single_switch") {
    initial_treated <- runif(n_groups) < initial_treat_frac
    switch_times <- sample.int(n_periods - 1, n_groups, replace = TRUE)  # 1..(n_periods-1)
    for (g in seq_len(n_groups)) {
      st <- switch_times[g] + 1L  # convert to 1-indexed switch period
      if (initial_treated[g]) {
        D[g, seq_len(st - 1L)] <- 1L
        D[g, st:n_periods] <- 0L
      } else {
        D[g, seq_len(st - 1L)] <- 0L
        D[g, st:n_periods] <- 1L
      }
    }
  } else if (pattern == "joiners_only") {
    switch_times <- sample.int(n_periods - 1, n_groups, replace = TRUE)
    for (g in seq_len(n_groups)) {
      st <- switch_times[g] + 1L
      D[g, st:n_periods] <- 1L
    }
  } else if (pattern == "leavers_only") {
    switch_times <- sample.int(n_periods - 1, n_groups, replace = TRUE)
    for (g in seq_len(n_groups)) {
      st <- switch_times[g] + 1L
      D[g, seq_len(st - 1L)] <- 1L
    }
  } else if (pattern == "mixed_single_switch") {
    switch_times <- sample.int(n_periods - 1, n_groups, replace = TRUE)
    n_joiners <- n_groups %/% 2
    for (g in seq_len(n_groups)) {
      st <- switch_times[g] + 1L
      if (g <= n_joiners) {
        D[g, st:n_periods] <- 1L
      } else {
        D[g, seq_len(st - 1L)] <- 1L
      }
    }
  } else if (pattern == "multi_path_reversible") {
    # Deterministic multi-path DGP designed for by_path R-parity:
    #   - 4 distinct joiner-style target paths with unequal frequencies
    #     (so top-k ranking produces unique ranks with no ties)
    #   - path assignment is a DETERMINISTIC FUNCTION OF F_g, so each
    #     cohort (D_{g,1}, F_g, S_g) contains switchers from a single
    #     path. This avoids cross-path cohort sharing in the
    #     cohort-recentered influence function, which otherwise blows
    #     out SE parity with R's re-run-per-path convention.
    #   - post-window treatment is stable at path[L_max+1] (clean
    #     control-pool eligibility — no post-window contamination)
    #
    # Each group:
    #   - F_g in [2, n_periods - L_max] so the length-(L_max+1) window
    #     [F_g-1, F_g-1+L_max] fits the panel
    #   - path is determined by F_g: two F_g values per path (groups in
    #     {F_g in {2,3}} share path 1, {F_g in {4,5}} share path 2,
    #     {F_g in {6}} = path 3, {F_g in {7}} = path 4); within a path,
    #     F_g distribution yields n_groups * path_prop total groups
    max_switch <- n_periods - L_max - 1L
    stopifnot(max_switch >= 1L)
    # With n_periods=10, L_max=3: max_switch=6, F_g in [2,7] (6 values).

    target_paths <- list(
      c(0L, 1L, 1L, 1L),  # sustained on (rank 1)
      c(0L, 1L, 1L, 0L),  # on then off  (rank 2)
      c(0L, 1L, 0L, 0L),  # on briefly   (rank 3)
      c(0L, 1L, 0L, 1L)   # on-off-on    (rank 4, truncated under by_path=3)
    )
    stopifnot(length(target_paths[[1]]) == L_max + 1L)

    # Per-F_g path assignment: 2 F_g values for paths 1 & 2, 1 for paths
    # 3 & 4. This keeps each (D_{g,1}, F_g, S_g) cohort single-path.
    f_g_to_path <- c(1L, 1L, 2L, 2L, 3L, 4L)
    stopifnot(length(f_g_to_path) == max_switch)
    # Group counts per F_g (rank 1 > rank 2 > rank 3 > rank 4 with unique
    # ranks — no frequency ties, robust to R's undocumented tiebreak):
    #   F_g=2: 20 groups (path 1)
    #   F_g=3: 20 groups (path 1) → rank 1 has 40 switchers total
    #   F_g=4: 15 groups (path 2)
    #   F_g=5: 10 groups (path 2) → rank 2 has 25 switchers total
    #   F_g=6: 10 groups (path 3) → rank 3 has 10 switchers
    #   F_g=7:  5 groups (path 4) → rank 4 has  5 switchers (excluded by by_path=3)
    counts_per_F_g <- c(20L, 20L, 15L, 10L, 10L, 5L)
    stopifnot(sum(counts_per_F_g) == n_groups)

    # Build the group-to-(F_g, path) assignment, deterministic with seed
    g_idx <- 1L
    for (f_idx in seq_along(counts_per_F_g)) {
      F_g <- f_idx + 1L  # F_g in [2, 7] for f_idx in [1, 6]
      path_idx <- f_g_to_path[f_idx]
      target <- target_paths[[path_idx]]
      n_here <- counts_per_F_g[f_idx]
      for (k in seq_len(n_here)) {
        g <- g_idx
        # Pre-baseline [1 .. F_g-2]: initial state
        if (F_g >= 3L) {
          D[g, 1:(F_g - 2L)] <- target[1]
        }
        # Window [F_g-1 .. F_g-1+L_max]: exactly the target path
        for (j in 0:L_max) {
          D[g, F_g - 1L + j] <- target[j + 1L]
        }
        # Post-window [F_g+L_max .. n_periods]: stable at path[L_max+1]
        if (F_g + L_max <= n_periods) {
          D[g, (F_g + L_max):n_periods] <- target[L_max + 1L]
        }
        g_idx <- g_idx + 1L
      }
    }
  } else {
    stop(sprintf("Unknown pattern: %s", pattern))
  }

  # --- Append stable control cohorts ---
  if (n_never_treated > 0) {
    D <- rbind(D, matrix(0L, nrow = n_never_treated, ncol = n_periods))
  }
  if (n_always_treated > 0) {
    D <- rbind(D, matrix(1L, nrow = n_always_treated, ncol = n_periods))
  }
  n_groups <- nrow(D)

  # --- Generate fixed effects, true effects, outcomes ---
  group_fe <- rnorm(n_groups, mean = 0, sd = group_fe_sd)
  if (heterogeneous_effects) {
    true_effects <- matrix(rnorm(n_groups * n_periods, mean = treatment_effect, sd = effect_sd),
                           nrow = n_groups, ncol = n_periods)
  } else {
    true_effects <- matrix(treatment_effect, nrow = n_groups, ncol = n_periods)
  }
  true_effects[D == 0] <- 0.0

  period_arr <- 0:(n_periods - 1)
  noise <- matrix(rnorm(n_groups * n_periods, mean = 0, sd = noise_sd),
                  nrow = n_groups, ncol = n_periods)
  Y <- 10.0 + matrix(group_fe, nrow = n_groups, ncol = n_periods) +
       matrix(time_trend * period_arr, nrow = n_groups, ncol = n_periods, byrow = TRUE) +
       true_effects + noise

  # --- Build long-format data frame ---
  group_col <- rep(seq_len(n_groups) - 1L, each = n_periods)  # 0-indexed
  period_col <- rep(period_arr, n_groups)
  treatment_col <- as.vector(t(D))
  outcome_col <- as.vector(t(Y))

  data.frame(
    group = group_col,
    period = period_col,
    treatment = treatment_col,
    outcome = outcome_col
  )
}

# ---------------------------------------------------------------------------
# Helper: extract DID_1 (l=1) results from did_multiplegt_dyn output
# ---------------------------------------------------------------------------
extract_dcdh_l1 <- function(res) {
  # did_multiplegt_dyn returns a results object with $results$Effects matrix.
  # The Effects matrix has one row per effect (Effect_1, Effect_2, ...) and
  # columns: Estimate, SE, LB CI, UB CI, N, Switchers, N.w, Switchers.w. We
  # pull the row for "Effect_1" (the l=1 effect, == DID_M of AER 2020).
  effects <- res$results$Effects
  if (is.null(effects)) {
    stop("did_multiplegt_dyn returned no Effects; check the input data")
  }

  out <- list(
    overall_att = as.numeric(effects[1, "Estimate"]),
    overall_se = as.numeric(effects[1, "SE"]),
    overall_ci_lo = as.numeric(effects[1, "LB CI"]),
    overall_ci_hi = as.numeric(effects[1, "UB CI"]),
    n_switchers = as.numeric(effects[1, "N"])
  )

  # Placebo at lag 1 if available (Placebos has the same column layout)
  placebos <- res$results$Placebos
  if (!is.null(placebos) && nrow(placebos) >= 1) {
    out$placebo_effect <- as.numeric(placebos[1, "Estimate"])
    out$placebo_se <- as.numeric(placebos[1, "SE"])
    out$placebo_ci_lo <- as.numeric(placebos[1, "LB CI"])
    out$placebo_ci_hi <- as.numeric(placebos[1, "UB CI"])
  }

  out
}

# ---------------------------------------------------------------------------
# Helper: convert data frame to exportable list
# ---------------------------------------------------------------------------
export_data <- function(df) {
  list(
    group = as.numeric(df$group),
    period = as.numeric(df$period),
    treatment = as.numeric(df$treatment),
    outcome = as.numeric(df$outcome)
  )
}

scenarios <- list()

# Golden value datasets use n_groups=80 to keep the JSON file small (~50KB)
# while still being large enough to exercise inference.
N_GOLDEN <- 80

# ---------------------------------------------------------------------------
# Scenario 1: single_switch — mix of joiners and leavers
# ---------------------------------------------------------------------------
cat("  Scenario 1: single_switch_mixed\n")
d1 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 6,
                     pattern = "single_switch", seed = 101)
res1 <- did_multiplegt_dyn(
  df = d1, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 1, placebo = 1, ci_level = 95
)
scenarios$single_switch_mixed <- list(
  data = export_data(d1),
  params = list(pattern = "single_switch", n_groups = N_GOLDEN, n_periods = 6,
                seed = 101, effects = 1, placebo = 1, ci_level = 95),
  results = extract_dcdh_l1(res1)
)

# ---------------------------------------------------------------------------
# Scenario 2: joiners_only — pure staggered adoption
# ---------------------------------------------------------------------------
cat("  Scenario 2: joiners_only\n")
d2 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 6,
                     pattern = "joiners_only", seed = 102)
res2 <- did_multiplegt_dyn(
  df = d2, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 1, placebo = 1, ci_level = 95
)
scenarios$joiners_only <- list(
  data = export_data(d2),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 6,
                seed = 102, effects = 1, placebo = 1, ci_level = 95),
  results = extract_dcdh_l1(res2)
)

# ---------------------------------------------------------------------------
# Scenario 3: leavers_only — pure staggered removal
# ---------------------------------------------------------------------------
cat("  Scenario 3: leavers_only\n")
d3 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 6,
                     pattern = "leavers_only", seed = 103)
res3 <- did_multiplegt_dyn(
  df = d3, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 1, placebo = 1, ci_level = 95
)
scenarios$leavers_only <- list(
  data = export_data(d3),
  params = list(pattern = "leavers_only", n_groups = N_GOLDEN, n_periods = 6,
                seed = 103, effects = 1, placebo = 1, ci_level = 95),
  results = extract_dcdh_l1(res3)
)

# ---------------------------------------------------------------------------
# Scenario 4: mixed_single_switch — deterministic 50/50 joiners/leavers
# ---------------------------------------------------------------------------
cat("  Scenario 4: mixed_single_switch\n")
d4 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 6,
                     pattern = "mixed_single_switch", seed = 104)
res4 <- did_multiplegt_dyn(
  df = d4, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 1, placebo = 1, ci_level = 95
)
scenarios$mixed_single_switch <- list(
  data = export_data(d4),
  params = list(pattern = "mixed_single_switch", n_groups = N_GOLDEN, n_periods = 6,
                seed = 104, effects = 1, placebo = 1, ci_level = 95),
  results = extract_dcdh_l1(res4)
)

# ---------------------------------------------------------------------------
# Scenario 5: hand-calculable 4-group panel from the worked example.
# This is the panel used by tests/test_methodology_chaisemartin_dhaultfoeuille.py
# in test_hand_calculable_4group_3period_joiners_and_leavers. We capture R's
# answer here so the Python test can assert exact agreement.
# ---------------------------------------------------------------------------
cat("  Scenario 5: hand_calculable_worked_example\n")
d5 <- data.frame(
  group =     c(1, 1, 1,  2, 2, 2,  3, 3, 3,  4, 4, 4),
  period =    c(0, 1, 2,  0, 1, 2,  0, 1, 2,  0, 1, 2),
  treatment = c(0, 1, 1,  1, 1, 0,  0, 0, 0,  1, 1, 1),
  outcome   = c(10, 13, 14,  10, 11, 9,  10, 11, 12,  10, 11, 12)
)
res5 <- did_multiplegt_dyn(
  df = d5, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 1, placebo = 0, ci_level = 95
)
scenarios$hand_calculable_worked_example <- list(
  data = export_data(d5),
  params = list(description = "4-group hand-calculable panel from plan worked example",
                effects = 1, placebo = 0, ci_level = 95,
                expected_did_m = 2.5, expected_did_plus = 2.0, expected_did_minus = 3.0),
  results = extract_dcdh_l1(res5)
)

# ---------------------------------------------------------------------------
# Phase 2: Multi-horizon scenarios (effects > 1)
# ---------------------------------------------------------------------------

# Helper: extract multi-horizon results from did_multiplegt_dyn output
extract_dcdh_multi <- function(res, n_effects, n_placebos = 0) {
  effects <- res$results$Effects
  if (is.null(effects)) {
    stop("did_multiplegt_dyn returned no Effects; check the input data")
  }

  out <- list(effects = list(), placebos = list())

  for (i in seq_len(min(n_effects, nrow(effects)))) {
    out$effects[[as.character(i)]] <- list(
      overall_att = as.numeric(effects[i, "Estimate"]),
      overall_se = as.numeric(effects[i, "SE"]),
      overall_ci_lo = as.numeric(effects[i, "LB CI"]),
      overall_ci_hi = as.numeric(effects[i, "UB CI"]),
      n_switchers = as.numeric(effects[i, "N"])
    )
  }

  placebos <- res$results$Placebos
  if (!is.null(placebos) && n_placebos > 0) {
    for (i in seq_len(min(n_placebos, nrow(placebos)))) {
      out$placebos[[as.character(i)]] <- list(
        effect = as.numeric(placebos[i, "Estimate"]),
        se = as.numeric(placebos[i, "SE"]),
        ci_lo = as.numeric(placebos[i, "LB CI"]),
        ci_hi = as.numeric(placebos[i, "UB CI"])
      )
    }
  }

  out
}

# Scenario 6: joiners_only multi-horizon (L_max=3, placebo=3)
# Uses n_periods=8 to give enough room for 3 positive + 3 placebo horizons
cat("  Scenario 6: joiners_only_multi_horizon\n")
d6 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                     pattern = "joiners_only", seed = 106)
res6 <- did_multiplegt_dyn(
  df = d6, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 3, placebo = 3, ci_level = 95
)
scenarios$joiners_only_multi_horizon <- list(
  data = export_data(d6),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 8,
                seed = 106, effects = 3, placebo = 3, ci_level = 95),
  results = extract_dcdh_multi(res6, n_effects = 3, n_placebos = 3)
)

# Scenario 7: leavers_only multi-horizon (L_max=3, placebo=3)
cat("  Scenario 7: leavers_only_multi_horizon\n")
d7 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                     pattern = "leavers_only", seed = 107)
res7 <- did_multiplegt_dyn(
  df = d7, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 3, placebo = 3, ci_level = 95
)
scenarios$leavers_only_multi_horizon <- list(
  data = export_data(d7),
  params = list(pattern = "leavers_only", n_groups = N_GOLDEN, n_periods = 8,
                seed = 107, effects = 3, placebo = 3, ci_level = 95),
  results = extract_dcdh_multi(res7, n_effects = 3, n_placebos = 3)
)

# Scenario 8: mixed_single_switch multi-horizon (L_max=5, placebo=4)
# Uses n_periods=10 for far horizons
cat("  Scenario 8: mixed_single_switch_multi_horizon\n")
d8 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 10,
                     pattern = "mixed_single_switch", seed = 108)
res8 <- did_multiplegt_dyn(
  df = d8, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 5, placebo = 4, ci_level = 95
)
scenarios$mixed_single_switch_multi_horizon <- list(
  data = export_data(d8),
  params = list(pattern = "mixed_single_switch", n_groups = N_GOLDEN, n_periods = 10,
                seed = 108, effects = 5, placebo = 4, ci_level = 95),
  results = extract_dcdh_multi(res8, n_effects = 5, n_placebos = 4)
)

# Scenario 9: joiners_only long panel multi-horizon (L_max=5, placebo=5)
# Uses n_periods=12 and n_groups=80 for thorough coverage
cat("  Scenario 9: joiners_only_long_multi_horizon\n")
d9 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 12,
                     pattern = "joiners_only", seed = 109)
res9 <- did_multiplegt_dyn(
  df = d9, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 5, placebo = 5, ci_level = 95
)
scenarios$joiners_only_long_multi_horizon <- list(
  data = export_data(d9),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 12,
                seed = 109, effects = 5, placebo = 5, ci_level = 95),
  results = extract_dcdh_multi(res9, n_effects = 5, n_placebos = 5)
)

# ---------------------------------------------------------------------------
# Phase 3: Covariate and linear-trends scenarios
# ---------------------------------------------------------------------------

# Helper: add a covariate column to a panel. The covariate is correlated with
# switch timing (confounding) but the true effect is constant.
add_covariate <- function(df, seed = 42, x_effect = 1.5) {
  set.seed(seed)
  n <- nrow(df)
  groups <- unique(df$group)
  # Group-level base value (correlated with which groups switch)
  x_base <- setNames(rnorm(length(groups), 0, 1), groups)
  # Time-varying component
  df$X1 <- x_base[as.character(df$group)] + 0.3 * df$period + rnorm(n, 0, 0.2)
  # Add covariate effect to outcome
  df$outcome <- df$outcome + x_effect * df$X1
  df
}

# Scenario 10: joiners_only with controls (L_max=2)
cat("  Scenario 10: joiners_only_controls\n")
d10 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                      pattern = "joiners_only", seed = 110)
d10 <- add_covariate(d10, seed = 210, x_effect = 1.5)
res10 <- did_multiplegt_dyn(
  df = d10, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 2, placebo = 1, ci_level = 95,
  controls = "X1"
)
scenarios$joiners_only_controls <- list(
  data = list(
    group = as.numeric(d10$group),
    period = as.numeric(d10$period),
    treatment = as.numeric(d10$treatment),
    outcome = as.numeric(d10$outcome),
    X1 = as.numeric(d10$X1)
  ),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 8,
                seed = 110, effects = 2, placebo = 1, ci_level = 95,
                controls = "X1"),
  results = extract_dcdh_multi(res10, n_effects = 2, n_placebos = 1)
)

# Scenario 11: joiners_only with trends_lin (L_max=2)
cat("  Scenario 11: joiners_only_trends_lin\n")
d11 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                      pattern = "joiners_only", seed = 111)
# Add group-specific linear trends to outcome
set.seed(311)
groups11 <- unique(d11$group)
g_trends <- setNames(rnorm(length(groups11), 0, 0.5), groups11)
d11$outcome <- d11$outcome + g_trends[as.character(d11$group)] * d11$period
res11 <- did_multiplegt_dyn(
  df = d11, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 2, placebo = 1, ci_level = 95,
  trends_lin = TRUE
)
scenarios$joiners_only_trends_lin <- list(
  data = export_data(d11),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 8,
                seed = 111, effects = 2, placebo = 1, ci_level = 95,
                trends_lin = TRUE),
  results = extract_dcdh_multi(res11, n_effects = 2, n_placebos = 1)
)

# Scenario 12: joiners_only with both controls and trends_lin (L_max=2)
cat("  Scenario 12: joiners_only_controls_trends_lin\n")
d12 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                      pattern = "joiners_only", seed = 112)
d12 <- add_covariate(d12, seed = 212, x_effect = 1.5)
# Add group-specific linear trends
set.seed(312)
groups12 <- unique(d12$group)
g_trends12 <- setNames(rnorm(length(groups12), 0, 0.5), groups12)
d12$outcome <- d12$outcome + g_trends12[as.character(d12$group)] * d12$period
res12 <- did_multiplegt_dyn(
  df = d12, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 2, placebo = 1, ci_level = 95,
  controls = "X1", trends_lin = TRUE
)
scenarios$joiners_only_controls_trends_lin <- list(
  data = list(
    group = as.numeric(d12$group),
    period = as.numeric(d12$period),
    treatment = as.numeric(d12$treatment),
    outcome = as.numeric(d12$outcome),
    X1 = as.numeric(d12$X1)
  ),
  params = list(pattern = "joiners_only", n_groups = N_GOLDEN, n_periods = 8,
                seed = 112, effects = 2, placebo = 1, ci_level = 95,
                controls = "X1", trends_lin = TRUE),
  results = extract_dcdh_multi(res12, n_effects = 2, n_placebos = 1)
)

# ---------------------------------------------------------------------------
# Phase 3 extension: by_path per-path event-study disaggregation
# ---------------------------------------------------------------------------

# Helper: extract per-path by_path results. When did_multiplegt_dyn is
# called with by_path=k, the result object has no $results slot; instead
# per-path results live at res$by_level_1, res$by_level_2, ... in rank
# order (1 = most frequent observed path). res$by_levels is a character
# vector of comma-joined path labels (e.g. "0,1,1,1") in the same order.
extract_dcdh_by_path <- function(res, n_effects) {
  by_levels <- res$by_levels
  out <- list()
  for (i in seq_along(by_levels)) {
    slot <- res[[paste0("by_level_", i)]]
    effects <- slot$results$Effects
    horizons <- list()
    for (h in seq_len(min(n_effects, nrow(effects)))) {
      horizons[[as.character(h)]] <- list(
        effect = as.numeric(effects[h, "Estimate"]),
        se = as.numeric(effects[h, "SE"]),
        ci_lo = as.numeric(effects[h, "LB CI"]),
        ci_hi = as.numeric(effects[h, "UB CI"]),
        n_switchers = as.numeric(effects[h, "Switchers"]),
        n_obs = as.numeric(effects[h, "N"])
      )
    }
    out[[i]] <- list(
      path = by_levels[i],
      frequency_rank = i,
      horizons = horizons
    )
  }
  list(by_path = out)
}

# Scenario 13: mixed_single_switch + by_path=2 (basic 2-path case).
# The mixed_single_switch DGP produces joiners (path 0,1,1,1) and
# leavers (path 1,0,0,0) as its only two observed paths at L_max=3, so
# by_path=2 captures both and tests core per-path parity.
cat("  Scenario 13: mixed_single_switch_by_path\n")
d13 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 8,
                      pattern = "mixed_single_switch", seed = 113)
res13 <- did_multiplegt_dyn(
  df = d13, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 3, by_path = 2, ci_level = 95
)
scenarios$mixed_single_switch_by_path <- list(
  data = export_data(d13),
  params = list(pattern = "mixed_single_switch", n_groups = N_GOLDEN,
                n_periods = 8, seed = 113, effects = 3, by_path = 2,
                ci_level = 95),
  results = extract_dcdh_by_path(res13, n_effects = 3)
)

# Scenario 14: multi_path_reversible + by_path=3 (top-k ranking case).
# The `multi_path_reversible` pattern is a DETERMINISTIC multi-path DGP:
# path assignment is a fixed function of F_g (so every (D_{g,1}, F_g,
# S_g) cohort contains switchers from a single path), path proportions
# are fixed at 20/20/15/10/10/5 across the 6 F_g values, and
# post-window treatment is stable at path[L_max+1]. by_path=3 exercises
# top-k selection when observed paths exceed k (4 observed paths, top-3
# selected). n_periods=10 gives every switch_time a complete length-
# (L_max+1) window. The old `p_switch`-driven random-toggle variant
# (pre-PR) blew out SE parity with R via cross-path cohort mixing;
# see the REGISTRY.md `Note (Phase 3 by_path ...)` Deviation bullet.
cat("  Scenario 14: multi_path_reversible_by_path\n")
d14 <- gen_reversible(n_groups = N_GOLDEN, n_periods = 10,
                      pattern = "multi_path_reversible", seed = 114,
                      L_max = 3)
res14 <- did_multiplegt_dyn(
  df = d14, outcome = "outcome", group = "group", time = "period",
  treatment = "treatment", effects = 3, by_path = 3, ci_level = 95
)
scenarios$multi_path_reversible_by_path <- list(
  data = export_data(d14),
  params = list(pattern = "multi_path_reversible", n_groups = N_GOLDEN,
                n_periods = 10, seed = 114, effects = 3, by_path = 3,
                ci_level = 95),
  results = extract_dcdh_by_path(res14, n_effects = 3)
)

# ---------------------------------------------------------------------------
# Write output
# ---------------------------------------------------------------------------
dir.create(dirname(output_path), showWarnings = FALSE, recursive = TRUE)
writeLines(
  toJSON(
    list(scenarios = scenarios,
         generator = "generate_reversible_did_data v1",
         dcdh_package = paste0("DIDmultiplegtDYN ", utils::packageVersion("DIDmultiplegtDYN"))),
    auto_unbox = TRUE,
    digits = 10,
    pretty = TRUE
  ),
  output_path
)
cat(sprintf("dCDH golden values written to %s\n", output_path))
cat(sprintf("File size: %.1f KB\n", file.info(output_path)$size / 1024))
cat(sprintf("Scenarios: %d\n", length(scenarios)))