Numerical Replication of the Sequeira Application • compdid

This guide is for numerical replication of Sant’Anna and Xu (2026), Table 4. It is intentionally different from the applied user guide:

vignette("sequeira-application", package = "compdid")

The applied guide shows how external researchers should use compdid with current best-practice defaults. This replication guide instead recreates the paper-final application output. To make the published numbers reproducible, it uses the archived paper-final local-polynomial nuisance fitted values included in sequeira_application_reference, then recomputes the ATT estimators, stationarity tests, and clustered bootstrap with exported package functions.

This distinction matters. The current package deliberately uses more guarded defaults than the historical application scripts. In particular, modern GPS fits normalize rows by default, and the local multinomial-logit optimizer has more stable controls. Those choices are better for new analyses, but they can change the Sequeira application estimates slightly.

Data And Reference Objects

library(compdid)

data("sequeira_bribes", package = "compdid")
data("sequeira_application_reference", package = "compdid")

ref <- sequeira_application_reference
ref$table4
ref$tests

The reference object contains:

the paper-final Table 4 values;
the paper-final local-polynomial GPS and OR fitted values;
paper-final local-polynomial bandwidths;
the seed and number of bootstrap draws used in the application;
source hashes used in the package audit.

Preprocessing

The replication uses the same deterministic transformations as the original application scripts.

normalize01 <- function(data, cols) {
  out <- data
  for (col in cols) {
    x <- as.numeric(out[[col]])
    rng <- range(x, na.rm = TRUE)
    denom <- rng[[2L]] - rng[[1L]]
    out[[col]] <- if (is.finite(denom) && denom > 0) {
      (x - rng[[1L]]) / denom
    } else {
      0
    }
  }
  out
}

prepare_locpol_data <- function(data, ref) {
  data$lba_value[2783L] <- 0
  data$monitor <- data$monitor - 1
  data$psi <- 2 - data$psi
  data$ltonnage <- log(data$tonnage + 1)

  keep <- c(
    ref$locpol$ynames,
    ref$locpol$tname,
    ref$locpol$dname,
    "lvalue_tonnage",
    "ltonnage",
    "lvalue_shipment_metical",
    "tariff2007",
    ref$locpol$xunames,
    ref$locpol$gname,
    "clear_agent",
    "hc_group"
  )
  data[complete.cases(data[, keep]), , drop = FALSE]
}

make_sequeira_dp <- function(data, ref, outcome_index) {
  x_cont <- ref$locpol$xcnames[[outcome_index]]
  data_norm <- normalize01(data, x_cont)
  d <- as.numeric(data_norm[[ref$locpol$dname]])
  post <- as.numeric(data_norm[[ref$locpol$tname]])
  covariates <- cbind(
    as.matrix(data_norm[, x_cont, drop = FALSE]),
    as.matrix(data_norm[, ref$locpol$xunames, drop = FALSE])
  )
  storage.mode(covariates) <- "double"

  list(
    y = as.numeric(data_norm[[ref$locpol$ynames[[outcome_index]]]]),
    d = d,
    post = post,
    dpost = cbind(
      d * post,
      d * (1 - post),
      (1 - d) * post,
      (1 - d) * (1 - post)
    ),
    covariates = covariates,
    dim_covariates = c(length(x_cont), length(ref$locpol$xunames), 0L),
    group = data_norm[[ref$locpol$gname]]
  )
}

app <- prepare_locpol_data(sequeira_bribes, ref)
stopifnot(nrow(app) == ref$locpol$n)

DR Rows

The DR rows in Table 4 are reproduced by passing the paper-final nuisance fitted values through the exported second-stage functions.

rows <- vector("list", length(ref$locpol$ynames))
tests <- vector("list", length(ref$locpol$ynames))

set.seed(ref$locpol$seed)
for (i in seq_along(ref$locpol$ynames)) {
  dp <- make_sequeira_dp(app, ref, i)
  ps <- ref$locpol$ps_fit[[i]]$fitted.values
  or <- ref$locpol$or_fit[[i]]$fitted.values

  stationary <- drdid_stationary(
    dp$y, dp$d, dp$post, ps, or,
    stabilized = TRUE,
    boot = FALSE,
    inffunc = TRUE
  )
  nonstationary <- drdid_nonstationary(
    dp$y, dp$d, dp$post, ps, or,
    stabilized = TRUE,
    boot = FALSE,
    inffunc = TRUE
  )
  clustered <- drdid_cluster_bootstrap(
    nonstationary_fit = nonstationary,
    stationary_fit = stationary,
    group = dp$group,
    nboot = ref$locpol$nboot,
    seed = NULL,
    alpha = 0.05,
    keep_boots = FALSE
  )

  rows[[i]] <- rbind(
    data.frame(
      estimator = "DR DiD tau_sz (no-compositional changes)",
      outcome = ref$table4$outcome[[i]],
      estimate = as.numeric(stationary$ATT),
      se = as.numeric(stationary$se),
      cluster_se = clustered$stationary$se
    ),
    data.frame(
      estimator = "DR DiD tau_dr (robust to compositional changes)",
      outcome = ref$table4$outcome[[i]],
      estimate = as.numeric(nonstationary$ATT),
      se = as.numeric(nonstationary$se),
      cluster_se = clustered$nonstationary$se
    )
  )
  tests[[i]] <- data.frame(
    outcome = ref$tests$outcome[[i]],
    wald_stat = as.numeric(clustered$stationarity_test$statistic),
    unclustered_p_value = as.numeric(clustered$stationarity_test$p.value),
    clustered_p_value = clustered$stationarity_test$clustered.p.value
  )
}

dr_table <- do.call(rbind, rows)
dr_tests <- do.call(rbind, tests)
dr_table
dr_tests

The corresponding installed check script runs the same logic and asserts that the package output matches the stored paper-final reference:

script <- system.file("replication", "verify_sequeira_saved_nuisance.R",
                      package = "compdid")
system2(file.path(R.home("bin"), "Rscript"), script)

TWFE Rows

The TWFE rows use twfe_did_rc(). The installed application script contains the full fixed-effect design used in the paper and writes the corresponding CSV/RData files:

script <- system.file("replication", "sequeira_application.R", package = "compdid")
system2(
  file.path(R.home("bin"), "Rscript"),
  c(script, "--saved-nuisance", "--out-dir=/tmp/compdid_sequeira_exact")
)

This writes:

result_sequeira_fe_compdid.RData;
result_sequeira_fe_compdid.csv;
result_sequeira_ml_compdid.RData;
result_sequeira_ml_compdid.csv;
table4_paper_check_estimates.csv;
table4_paper_check_tests.csv.

Why Not Refit Everything?

The paper-final GPS rows are clipped but not row-normalized; their row sums range from 1.00 to 1.02. Current package defaults normalize GPS rows because they are probabilities over the four treatment-period cells. That is the right default for new users.

The historical local-polynomial GPS cross-validation is also toolchain-sensitive. On the package development machine, a full from-scratch historical CV audit completed in about 2 hours 34 minutes. It matched the TWFE rows and the outcome-regression local-polynomial fits, but the live GPS CV optimizer selected different bandwidths for several outcomes. The DR rows therefore did not numerically match the paper table. This is why exact replication uses the archived paper-final nuisance fitted values.

Advanced Historical Sensitivity

Researchers who want to understand the role of the old GPS convention can use real estimator controls directly. The example below shows the relevant controls at fixed paper-final bandwidths. This is a sensitivity exercise, not the recommended analysis path.

i <- 1L
dp <- make_sequeira_dp(app, ref, i)
paper_bws <- c(ref$locpol$ps_fit[[i]]$bws[, 1], -1)

ps_audit <- locpol_ps_fit(
  dp,
  bws = paper_bws,
  cv_method = "cv.ml",
  lp_type = "logit",
  list_control = list(
    ps_min = 1e-2,
    normalize_ps = FALSE,
    mnl_start = "zero",
    mnl_stable = FALSE,
    mnl_gradtol = NULL,
    mnl_maxit = 10000L,
    mnl_reltol = sqrt(.Machine$double.eps)
  )
)

range(rowSums(ps_audit$fitted.values))

For new applications, use normalize_ps = TRUE and the default guarded optimizer settings. Use the best-practice Sequeira guide for that workflow.