Simulation Study

/simulation-study — Monte Carlo Simulation Study

Design and run a Monte Carlo experiment that characterizes an estimator's finite-sample behavior, then review it for the bugs that quietly invalidate simulation evidence.

Input: $ARGUMENTS — a description of the estimator(s) and DGP to study (e.g., "compare 2SLS vs LIML under weak instruments with heteroskedasticity"), or a pointer to an existing script/paper whose simulation you want to reproduce or extend.

Constraints

• Follow .claude/rules/simulation-conventions.md — the simulation contract (DGP, truth, estimand, MCSE, assumption regime) is non-negotiable.
• Declare the assumption regime in the script header and respect the firewall — an out-of-assumption run never supports a within-assumption claim (simulation-conventions.md §2).
• Follow .claude/rules/r-code-conventions.md for general R standards (header, library() at top, relative paths, numerical discipline).
• Save the script to scripts/R/ with a numbered, descriptive name (e.g., scripts/R/sim_2sls_vs_liml.R).
• Save outputs (per-rep raw tibble, summary table, figures) to output/.
• saveRDS() the per-replication raw results, not just the summary — re-aggregation and the review pass need them.
• Run the sim-reviewer agent on the generated script before presenting results, then address Critical/High findings.

Workflow Phases

Phase 0: Pre-Flight Report

Before writing any code, produce a Pre-Flight Report showing you have pinned down the experiment. This prevents the most common failure mode — a beautiful results table built on a mismatched estimand or a coverage-against-the-estimate bug.

## Pre-Flight Report — Simulation Design

**Research question:** [what finite-sample property is being demonstrated]
**Target estimand:** [ATT / ATE / coefficient θ — and how its TRUE value is computed from the DGP params]
**Maintained assumptions:** [the FULL list the estimator(s) under study require — A1 … An, every one]
**Regime:** [IN-ASSUMPTION — all hold | OUT-OF-ASSUMPTION — relaxes A[k] only, severity grid {…}, targeting pseudo-estimand …]
**Verification:** [per assumption, the checkable property of the DGP that establishes it — by construction or by an assertion]
**DGP:** [structure + the parameters that define it; what is held fixed vs. varied]
**Estimator grid:** [list each estimator + which estimand it targets + how it returns est/se/CI]
**Design grid:** [sample sizes, parameter values, scenarios to sweep]
**Replications R:** [value] → implied MCSE on coverage ≈ sqrt(0.95·0.05/R) = [value]
**Metrics:** bias, empirical SE, RMSE, coverage, size/power — each with MCSE
**Conventions read:** simulation-conventions.md, r-code-conventions.md

If the estimand or its true value is ambiguous, stop and ask before writing code.

If an assumption cannot be verified — you cannot name the property of the DGP that establishes it — the run is not IN-ASSUMPTION, and per the firewall no within-assumption claim may rest on it. Say so in the Pre-Flight Report rather than letting the header assert what was never checked.

Phase 1: The DGP

Write one parameterized function that returns a dataset. Compute and return (or store) the true target value from the parameters.

generate_data <- function(n, params) {
  # ... generate covariates, treatment, outcome from params ...
  list(data = df, truth = compute_truth(params))   # truth from params, never from an estimate
}

The header's Verified lines are earned here: every assumption the regime block claims holds by a check gets that check written into the script (a large-draw assertion, a condition number, a stopifnot() on the parameter bounds), run once at setup. An assumption whose verification exists only in the comment is asserted, not verified.

Phase 2: Estimator Grid

Each estimator is a function data -> list(est, se, ci_lo, ci_hi, converged). State the estimand each one targets; an estimator scored against a mismatched truth is a bug, not a finding.

Phase 3: Replication Engine

• set.seed(YYYYMMDD) once. For parallel reps use RNGkind("L'Ecuyer-CMRG") and furrr::furrr_options(seed = TRUE).
• One run = generate data → run every estimator → record a row per estimator with est, se, ci_lo, ci_hi, converged.
• Pre-allocate / bind results into a tibble of R × (#estimators) rows. Track non-convergence; never silently drop.

Phase 4: Metrics & Summary

Per estimator × scenario, against truth:

• Bias = mean(est) - truth (+ MCSE = sd(est)/sqrt(R))
• Empirical SE = sd(est); RMSE = sqrt(mean((est - truth)^2))
• Coverage = mean(ci_lo <= truth & truth <= ci_hi) (+ MCSE = sqrt(p(1-p)/R))
• Size / power = rejection rate under the null / alternative DGP
• Failures = count of non-converged reps

Build a tidy summary table; report MCSE next to every headline metric.

Phase 5: Figures

Use ggplot2 with the project theme: bias / coverage vs. sample size (or scenario), with reference lines (0 bias, nominal coverage). Transparent background, explicit dimensions (per r-code-conventions.md §4).