diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1d243d81..400da3db 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -38,6 +38,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     change; the machine-precision hetero/cluster lock is deferred — needs an unbalanced-DGP golden).
 
 ### Added
+- **`ImputationDiD` leave-one-out conservative variance** (`leave_one_out`, default `False`) — the
+  Borusyak-Jaravel-Spiess (2024) Supplementary Appendix A.9 finite-sample refinement. The non-LOO
+  auxiliary aggregate `tau_tilde_g` (eq. 8) is built from the fitted `tau_hat_it` and so partially
+  overfits to the noise `epsilon_it`, biasing the conservative variance downward; `leave_one_out=True`
+  recomputes each unit's group aggregate excluding that unit (implemented efficiently by rescaling each
+  treated auxiliary residual by `1/(1 - v_ig**2/sum_j v_jg**2)`, exactly equivalent to the direct
+  leave-one-out at the per-unit cluster sum), yielding a larger, less-downward-biased SE (Prop. A8:
+  unbiased for an upper bound at the default unit clustering). Point estimates are unchanged. A group
+  with a single positive-weight unit (LOO undefined, App. A.9 fn. 51) keeps the non-LOO residual with a
+  `UserWarning`. `leave_one_out` is recorded on `ImputationDiDResults` (and `to_dict()` / `summary()`);
+  replicate-weight survey designs raise `NotImplementedError` (their variance bypasses the
+  influence-function path). Default `False` preserves R `didimputation` parity.
 - **`ContinuousDiD` lowest-dose-as-control** (`control_group="lowest_dose"`, CGBS 2024 Remark 3.1) for
   settings with no untreated group (`P(D=0) = 0`): the lowest-dose group `d_L` becomes the comparison
   and the estimand is `ATT(d) − ATT(d_L)` (with `ATT(d_L) = 0` the omitted reference). It is a
diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md
index 0dc56f17..998adabc 100644
--- a/METHODOLOGY_REVIEW.md
+++ b/METHODOLOGY_REVIEW.md
@@ -561,7 +561,7 @@ and covariate-adjusted specifications.)
 - Multiplier bootstrap on the Theorem-3 influence function (library extension, not in the paper).
 - Survey-design TSL variance on the influence function (library extension).
 - NaN inference for undefined statistics; Proposition-5 refuse-to-estimate (NaN + warning).
-- Leave-one-out variance refinement (Supplementary Appendix A.9) not implemented (finite-sample refinement; tracked as future work).
+- Leave-one-out variance refinement (Supplementary Appendix A.9) implemented as the opt-in `leave_one_out` parameter (default `False`, preserving R `didimputation` parity); efficient residual rescale `1/(1 - v_ig^2/sum_j v_jg^2)` exactly equivalent to direct leave-one-out (see REGISTRY ImputationDiD Note).
 
 **R Comparison Results** (`didimputation` v0.5.0, fixed-seed panel, `benchmarks/data/didimputation_golden.json`):
 
diff --git a/TODO.md b/TODO.md
index 6ca2f22d..17291ace 100644
--- a/TODO.md
+++ b/TODO.md
@@ -30,7 +30,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 |-------|----------|--------|--------|----------|
 | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
 | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0` (`control_group="lowest_dose"`) — **DONE** (discrete + continuous mass-point, single-cohort; estimand `ATT(d)−ATT(d_L)`; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low |
-| `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low |
 | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low |
 | `CallawaySantAnna` no-covariate ipw treats the propensity as unconditional (no correction term); R's `did` fits an intercept-only logit whose estimation effect is identically zero in the IF, so values match — decide whether to document-only (current REGISTRY note) or mirror R's code path for structural parity. | `staggered.py::_ipw_estimation` (no-cov branch) | CS-scaling | Quick | Low |
 | Fold the R `did` 2.5.1 ipw aggregation yardsticks (hardcoded with provenance in `test_golden_ipw_aggregation_se_vs_r_did_251`) into `csdid_golden_values.json` on the next fixture regeneration — the generator already emits the ipw `aggte` blocks; switch the test to read the JSON. | `benchmarks/R/generate_csdid_test_values.R`, `tests/test_csdid_ported.py` | CS-scaling | Quick | Low |
diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt
index 3a09b257..1f042836 100644
--- a/diff_diff/guides/llms-full.txt
+++ b/diff_diff/guides/llms-full.txt
@@ -392,6 +392,7 @@ ImputationDiD(
     horizon_max: int | None = None,              # Max event-study horizon
     aux_partition: str = "cohort_horizon",        # "cohort_horizon", "cohort", or "horizon"
     pretrends: bool = False,                     # Include pre-treatment horizons in event study
+    leave_one_out: bool = False,                 # BJS 2024 App. A.9 leave-one-out finite-sample variance (larger, less-downward-biased SE)
 )
 ```
 
diff --git a/diff_diff/imputation.py b/diff_diff/imputation.py
index 9e703073..9717378a 100644
--- a/diff_diff/imputation.py
+++ b/diff_diff/imputation.py
@@ -127,6 +127,27 @@ class ImputationDiD(ImputationDiDBootstrapMixin):
         pre-trends assessment. Pre-period effects should be ~0 under
         parallel trends. Only affects event_study aggregation; overall
         ATT and group aggregation are unchanged.
+    leave_one_out : bool, default=False
+        If True, apply the Borusyak-Jaravel-Spiess (2024) Supplementary
+        Appendix A.9 leave-one-out finite-sample refinement to the
+        conservative variance. The non-LOO auxiliary aggregate ``tau_tilde_g``
+        is built from the fitted ``tau_hat_it`` and thus partially overfits to
+        the noise ``epsilon_it``, biasing the variance downward. LOO recomputes
+        each unit's group aggregate excluding that unit -- implemented
+        efficiently by rescaling each treated auxiliary residual by
+        ``1 / (1 - v_ig**2 / sum_j v_jg**2)`` (App. A.9), which is exactly
+        equivalent to the direct leave-one-out at the per-unit cluster sum.
+        Yields a larger, less-downward-biased SE (Prop. A8: unbiased for an
+        upper bound). Default False preserves R ``didimputation`` parity; the
+        refinement is an option in the authors' Stata ``did_imputation``. LOO
+        is undefined for a group with a single positive-weight unit (App. A.9
+        footnote 51): such groups fall back to the non-LOO residual with a
+        UserWarning. The Prop. A8 direction (LOO >= non-LOO) is guaranteed at
+        the default unit clustering; coarser ``cluster=`` / analytical
+        ``survey_design=`` / ``n_bootstrap`` compositions apply the same rescale
+        but are a library extension beyond the paper's derivation.
+        Replicate-weight survey designs raise ``NotImplementedError`` (their
+        variance bypasses the influence-function path where the rescale lives).
 
     Attributes
     ----------
@@ -182,6 +203,7 @@ def __init__(
         horizon_max: Optional[int] = None,
         aux_partition: str = "cohort_horizon",
         pretrends: bool = False,
+        leave_one_out: bool = False,
     ):
         if rank_deficient_action not in ("warn", "error", "silent"):
             raise ValueError(
@@ -199,6 +221,7 @@ def __init__(
                 f"got '{aux_partition}'"
             )
         self._validate_vcov_type(vcov_type)
+        self._validate_leave_one_out(leave_one_out)
 
         self.anticipation = anticipation
         self.alpha = alpha
@@ -211,6 +234,7 @@ def __init__(
         self.horizon_max = horizon_max
         self.aux_partition = aux_partition
         self.pretrends = pretrends
+        self.leave_one_out = leave_one_out
 
         self.is_fitted_ = False
         self.results_: Optional[ImputationDiDResults] = None
@@ -275,6 +299,7 @@ def fit(
         # mutations (e.g. set_params(vcov_type="classical")) are re-checked
         # at use rather than silently accepted by the parameter setter.
         self._validate_vcov_type(self.vcov_type)
+        self._validate_leave_one_out(self.leave_one_out)
 
         # Validate inputs
         required_cols = [outcome, unit, time, first_treat]
@@ -345,6 +370,22 @@ def fit(
                 "the design structure implicitly) or use a non-replicate "
                 "survey design (with explicit strata/psu/fpc)."
             )
+        # Reject replicate-weight + leave_one_out=: the BJS 2024 App. A.9
+        # refinement rescales the conservative influence-function auxiliary
+        # residuals, but replicate-weight variance is computed by per-replicate
+        # point-estimate refits (not the IF path), so leave_one_out would
+        # silently have no effect. Fail-closed (no-silent-failures).
+        if _uses_replicate_imp and self.leave_one_out:
+            raise NotImplementedError(
+                "ImputationDiD(leave_one_out=True) is not supported with "
+                "replicate-weight survey designs. The leave-one-out refinement "
+                "(Borusyak, Jaravel & Spiess 2024, Supp. App. A.9) rescales the "
+                "conservative influence-function residuals, but replicate-weight "
+                "variance is computed by per-replicate refits and does not use "
+                "that path — leave_one_out would silently have no effect. Use a "
+                "non-replicate (Taylor-linearization) survey design, or "
+                "leave_one_out=False."
+            )
         # Validate within-unit constancy for panel survey designs
         if resolved_survey is not None:
             _validate_unit_constant_survey(data, unit, survey_design)
@@ -968,6 +1009,7 @@ def _refit_imp(w_r):
             vcov_type=self.vcov_type,
             cluster_name=_cluster_name_for_results,
             n_clusters=_n_clusters_for_results,
+            leave_one_out=self.leave_one_out,
         )
 
         self.is_fitted_ = True
@@ -1683,6 +1725,8 @@ def _compute_auxiliary_residuals_treated(
         # _compute_target_weights) poison its whole group via 0 * NaN = NaN. The
         # previous observation-level pandas sum relied on skipna to drop them.
         contrib = (v_treated != 0.0) & np.isfinite(tau_hat)
+        loo_factor: Optional[pd.Series] = None
+        n_single_loo = 0
         if contrib.any():
             inner = pd.DataFrame(
                 {
@@ -1702,6 +1746,12 @@ def _compute_auxiliary_residuals_treated(
             )
             den_ok = per_group["den"].abs() >= 1e-15
             tau_tilde_map = (per_group["num"] / per_group["den"]).where(den_ok)
+            # BJS 2024 App. A.9 leave-one-out refinement: rescale each treated
+            # residual by 1/(1 - v_ig^2 / sum_j v_jg^2) (== the direct-LOO tau_tilde
+            # exactly, at the per-unit cluster sum). Reuses a_{i,g} = per_unit['a']
+            # and sum_j v_jg^2 = per_group['den']; applied to epsilon_treated below.
+            if self.leave_one_out:
+                loo_factor, n_single_loo = self._leave_one_out_factor(per_unit, per_group)
         else:
             tau_tilde_map = pd.Series(dtype=float)
 
@@ -1720,6 +1770,28 @@ def _compute_auxiliary_residuals_treated(
         # Auxiliary residuals
         epsilon_treated = tau_hat - tau_tilde
 
+        # Leave-one-out rescale (BJS 2024 App. A.9): map each treated obs to its
+        # (group, unit) factor and inflate the residual. Non-contributing rows
+        # (v_it == 0, psi == 0 anyway) and single-positive-weight-unit groups
+        # (LOO undefined, fn. 51) keep factor 1.0.
+        if self.leave_one_out and loo_factor is not None:
+            obs_index = pd.MultiIndex.from_arrays(
+                [group_codes, df_1[unit].values], names=["g", "u"]
+            )
+            factor_per_obs = loo_factor.reindex(obs_index).to_numpy(dtype=float)
+            factor_per_obs = np.where(np.isfinite(factor_per_obs), factor_per_obs, 1.0)
+            epsilon_treated = epsilon_treated * factor_per_obs
+            if n_single_loo > 0:
+                warnings.warn(
+                    f"leave_one_out=True: {n_single_loo} auxiliary group(s) have a single "
+                    f"positive-weight unit, where the leave-one-out variance is undefined "
+                    f"(Borusyak, Jaravel & Spiess 2024, Supp. App. A.9 fn. 51); those groups "
+                    f"keep the non-leave-out residual. A coarser aux_partition reduces "
+                    f"singleton groups.",
+                    UserWarning,
+                    stacklevel=2,
+                )
+
         return epsilon_treated
 
     def _compute_residuals_untreated(
@@ -2496,6 +2568,7 @@ def get_params(self) -> Dict[str, Any]:
             "horizon_max": self.horizon_max,
             "aux_partition": self.aux_partition,
             "pretrends": self.pretrends,
+            "leave_one_out": self.leave_one_out,
         }
 
     def set_params(self, **params) -> "ImputationDiD":
@@ -2507,6 +2580,75 @@ def set_params(self, **params) -> "ImputationDiD":
                 raise ValueError(f"Unknown parameter: {key}")
         return self
 
+    @staticmethod
+    def _validate_leave_one_out(leave_one_out: Any) -> None:
+        """Validate ``leave_one_out`` is a strict bool.
+
+        Called from ``__init__`` AND ``fit()`` so sklearn-style
+        ``set_params(leave_one_out=...)`` mutations are re-checked at use
+        time -- the naive ``set_params`` setter would otherwise accept a
+        truthy string (e.g. "yes") and silently run the LOO refinement.
+        """
+        if not isinstance(leave_one_out, bool):
+            raise TypeError(f"leave_one_out must be a bool, got {type(leave_one_out).__name__}")
+
+    @staticmethod
+    def _leave_one_out_factor(
+        per_unit: pd.DataFrame, per_group: pd.DataFrame
+    ) -> Tuple[pd.Series, int]:
+        """Per-(group, unit) leave-one-out residual-rescale factor (BJS 2024 A.9).
+
+        ``factor_{g,i} = 1 / (1 - v_ig**2 / sum_j v_jg**2)`` with
+        ``v_ig = per_unit['a']`` and ``sum_j v_jg**2 = per_group['den']``. This
+        rescale of ``epsilon_tilde_it`` reproduces the direct leave-one-out
+        aggregate ``tau_tilde_it^LO`` exactly at the per-unit cluster sum
+        ``psi_i = sum_t v_it * epsilon_tilde_it`` (App. A.9). A group with a
+        single positive-weight unit has ``v_ig**2 == sum_j v_jg**2`` so the
+        factor diverges (LOO undefined, App. A.9 fn. 51); those groups fall back
+        to ``1.0`` (non-LOO). A genuinely unit-dominated but >=2-unit group keeps
+        its large finite factor -- that is the paper's intended inflation.
+
+        Returns
+        -------
+        (factor : pd.Series indexed like ``per_unit`` (g, u), n_single_unit_groups : int)
+        """
+        a = per_unit["a"].to_numpy(dtype=float)
+        sq = a**2
+        g_level = per_unit.index.get_level_values("g")
+        u_level = per_unit.index.get_level_values("u")
+        den = per_group["den"].reindex(g_level).to_numpy(dtype=float)  # D_g per (g,u)
+
+        # A group is "singleton" for LOO (App. A.9 fn. 51) when fewer than two
+        # units carry positive squared weight -- covers a true 1-unit group AND
+        # the effective-singleton case (>=2 rows, only one with a_ig != 0).
+        pos_per_group = pd.Series(sq > 0.0, index=per_unit.index).groupby(level="g").sum()
+        single_groups = pos_per_group.index[pos_per_group < 2]
+        is_single = pd.Index(g_level).isin(single_groups)
+
+        den_ok = np.abs(den) >= 1e-15
+        # factor = D_g / (D_g - v_ig^2) = D_g / sum_{j!=i} v_jg^2. Compute the
+        # leave-one-out denominator as the sum of the OTHER units' squared
+        # weights -- NOT as D_g - v_ig^2 after forming the ratio: for a genuinely
+        # dominated (but >=2-unit) group the subtraction loses precision (and can
+        # cancel to 0/negative) in float64 -- a finite-but-wrong or silently
+        # non-LOO factor. The fast subtraction is accurate away from the
+        # near-cancellation boundary; wherever the leave-one-out mass is a tiny
+        # fraction of D (relative loss of >~1e-6), recompute it exactly as the
+        # drop-then-sum of the OTHER units' squared weights. At most one unit per
+        # group can be that dominant, so the recompute stays O(units).
+        other_mass = den - sq
+        suspect = (~is_single) & den_ok & (other_mass <= 1e-6 * den)
+        if suspect.any():
+            sq_series = pd.Series(sq, index=per_unit.index)
+            for pos in np.nonzero(suspect)[0]:
+                grp = sq_series.xs(g_level[pos], level="g")
+                other_mass[pos] = float(grp.drop(u_level[pos]).sum())
+        # Fall back to non-LOO (factor 1.0) only where LOO is genuinely undefined:
+        # a singleton group (fn. 51), a degenerate den, or no other positive mass.
+        fallback = is_single | ~den_ok | (other_mass <= 0.0)
+        factor = np.where(fallback, 1.0, den / np.where(fallback, 1.0, other_mass))
+        return pd.Series(factor, index=per_unit.index), int(len(single_groups))
+
     @staticmethod
     def _validate_vcov_type(vcov_type: str) -> None:
         """Validate ``vcov_type`` membership against ImputationDiD's
diff --git a/diff_diff/imputation_results.py b/diff_diff/imputation_results.py
index 38f2884d..762619d7 100644
--- a/diff_diff/imputation_results.py
+++ b/diff_diff/imputation_results.py
@@ -150,6 +150,9 @@ class ImputationDiDResults:
     vcov_type: str = field(default="hc1")
     cluster_name: Optional[str] = field(default=None)
     n_clusters: Optional[int] = field(default=None)
+    # BJS 2024 Supp. App. A.9 leave-one-out finite-sample variance refinement
+    # (opt-in). Recorded here so reported SEs are self-describing.
+    leave_one_out: bool = field(default=False)
 
     # --- Inference-field aliases (balance/external-adapter compatibility) ---
     @property
@@ -260,6 +263,8 @@ def summary(self, alpha: Optional[float] = None) -> str:
                 lines.append(f"{'Variance estimator:':<30} {vcov_label:>15}")
         if self.n_clusters is not None and self.bootstrap_results is None:
             lines.append(f"{'Number of clusters:':<30} {self.n_clusters:>15}")
+        if self.leave_one_out:
+            lines.append(f"{'Leave-one-out variance:':<30} {'A.9 (BJS 2024)':>15}")
 
         lines.append("")
 
@@ -504,6 +509,7 @@ def to_dict(self) -> Dict[str, Any]:
             "alpha": self.alpha,
             "anticipation": self.anticipation,
             "vcov_type": self.vcov_type,
+            "leave_one_out": self.leave_one_out,
         }
         if self.cluster_name is not None:
             result["cluster_name"] = self.cluster_name
diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
index d216ac96..393cbe01 100644
--- a/docs/methodology/REGISTRY.md
+++ b/docs/methodology/REGISTRY.md
@@ -1578,6 +1578,7 @@ where `W_it(h) = 1[K_it = h]` are lead indicators, estimated on `Omega_0` only.
 - **Bootstrap inference:** Uses multiplier bootstrap on the Theorem 3 influence function: `psi_i = sum_t v_it * epsilon_tilde_it`. Cluster-level psi sums are pre-computed for each aggregation target (overall, per-horizon, per-group), then perturbed with multiplier weights (Rademacher by default; configurable via `bootstrap_weights` parameter to use Mammen or Webb weights, matching CallawaySantAnna). This is a library extension (not in the paper) consistent with CallawaySantAnna/SunAbraham bootstrap patterns.
 - **Auxiliary residuals (Equation 8):** Implements the paper's *unit-clustered* Equation 8 aggregator, `tau_tilde_g = sum_i (sum_{t in G_g,i} v_it)(sum_{t in G_g,i} v_it * tau_hat_it) / sum_i (sum_{t in G_g,i} v_it)^2` (Borusyak-Jaravel-Spiess 2024, eq. 8, p. 3272; minimal-excess-variance derivation in Supplementary Appendix A.8): for each unit form the within-unit weight sum `a_{i,g}` and weighted-effect sum `b_{i,g}` over the unit's observations in group `g`, then combine across units. Groups partition `Omega_1` via `aux_partition` (default `"cohort_horizon"` = cohort × event-time; also `"cohort"` / `"horizon"`). Unimputable (NaN `tau_hat`) and off-target observations carry `v_it = 0` and are excluded from the aggregation — exact for finite `tau_hat` (a zero-weight row adds 0 to both `a` and `b`) and NaN-safe; a group with no contributing observations falls back to the unweighted group mean (a variance no-op, since `psi_g = sum_t v_it * eps_tilde_it = 0` there).
 - **Note (deviation from R):** R `didimputation::did_imputation` computes the auxiliary aggregator as `sum(v_it^2 * tau_hat_it) / sum(v_it^2)` grouped by cohort × event-time only (no partition control is exposed). At that partition each unit contributes at most one observation per group, so the paper's unit-clustered Equation 8 reduces exactly to `sum(v^2 * tau)/sum(v^2)` — i.e. diff-diff matches R at the default `aux_partition="cohort_horizon"` (pinned in `tests/test_methodology_imputation.py::TestImputationDiDParityR`). diff-diff additionally offers the coarser `aux_partition="cohort"` / `"horizon"` groupings (where a unit may contribute several observations to a group), which have no R analogue and are validated by hand-calculation. Both implement the same paper Equation 8; only the available partition granularity differs. The earlier observation-level mean `sum(v*tau)/sum(v)` (pre-3.5.x) coincided with this only under uniform within-group weights; it was corrected to the exact unit-clustered form during the ImputationDiD methodology validation.
+- **Note (leave-one-out variance refinement, Supp. App. A.9):** the opt-in `leave_one_out=True` applies the Borusyak-Jaravel-Spiess (2024) Supplementary Appendix A.9 finite-sample refinement. The non-LOO `tau_tilde_g` (eq. 8) is built from the fitted `tau_hat_it`, which contain the noise `epsilon_it`, so it partially overfits and the auxiliary residuals `epsilon_tilde_it` are too small, biasing the variance **downward**. LOO recomputes each unit's group aggregate excluding that unit, `tau_tilde_it^LO = sum_{j!=i} v_jg^2 T_jg / sum_{j!=i} v_jg^2`, implemented efficiently (A.9) by rescaling each treated residual `epsilon_tilde_it^LO = epsilon_tilde_it / (1 - v_ig^2 / sum_j v_jg^2)` — where `v_ig = sum_{t in G_g,i} v_it` and `sum_j v_jg^2` are already materialized in `_compute_auxiliary_residuals_treated` as `a_{i,g}` and `per_group['den']`. That rescale reproduces the direct-LOO per-unit cluster sum `psi_i = sum_t v_it * epsilon_tilde_it` **exactly** (`psi_i^rescale = a_{i,g}(T_ig - tau_tilde_g) D/(D - a^2) = psi_i^direct-LOO`; a machine-precision-verified identity — **paper-fidelity** to the source `tau_tilde_it^LO`, not merely internal consistency). At the **default unit clustering** LOO gives a larger, less-downward-biased SE (Prop. A8: unbiased for an upper bound on the true variance; an equal-weight K-unit group inflates residuals by exactly `K/(K-1)`). **Default `leave_one_out=False`** preserves R `didimputation` parity (which omits LOO). **Edge (App. A.9 fn. 51):** a group with a single positive-weight unit has an undefined LOO (rescale factor `1/0`); such groups keep the non-LOO residual and the fit emits one consolidated `UserWarning` (a coarser `aux_partition` reduces singletons); a genuinely unit-dominated `>=2`-unit group keeps its large finite factor (intended inflation). **Composition scope:** the rescale operates on `epsilon_tilde`, so it flows through the coarser-`cluster=` CR sum, the analytical survey PSU-TSL variance, and the multiplier bootstrap unchanged — but Prop. A8's upper-bound guarantee and the `LOO >= non-LOO` direction hold at the default unit clustering only (under a coarser `cluster=`, per-unit `psi` inflation can partially cancel), so those compositions are a documented library extension, not paper-derived. **Replicate-weight survey designs** (BRR / Fay / JK1 / JKn / SDR) raise `NotImplementedError` with `leave_one_out=True`: replicate variance is computed by per-replicate point-estimate refits, bypassing the conservative-IF residual path where the LOO rescale lives, so LOO would silently no-op (fail-closed, no-silent-failures). **Effective-singleton guard:** a group's singleton test counts units with *positive* squared weight (not raw rows), and the leave-one-out denominator is the exact sum of the *other* units' squared weights (not `D - v_ig^2`, which can cancel to `<= 0` in float64 for an extremely dominated `>=2`-unit group and would silently revert it to non-LOO). **Reference / validation:** the authors' Stata `did_imputation` ships the same option; there is no CI-runnable anchor (R `didimputation` omits LOO, Stata is not in CI), so validation is the exact psi-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`). Source: arXiv:2108.12419v5 App. A.9 (the REStud Supplementary Material is canonical); the main-article review's A.9 GAP is now filled (see the `borusyak-jaravel-spiess-2024-review.md` provenance note). The stronger Prop. A8 variant that *also* leaves out for the `delta_hat` covariate estimation (exact-unbiased upper bound) is noted but not implemented — the Stata option is the `tau_tilde` residual rescale only.
 - **Note:** The Step-1 iterative FE solver (`_iterative_fe`) routes through the shared bincount Gauss-Seidel helper `diff_diff.utils._iterative_fe_solve`, and the covariate/pre-trend within-transformation routes through the shared MAP engine `diff_diff.utils.demean_by_groups` (factorize-once + `np.bincount`, optional Rust kernel; group order `[time, unit]` preserving the historical time-then-unit sweep) — the same convergence contract, accumulation-order numerics (~1e-10 vs the pre-3.7 pandas loops, not bit-for-bit), and `max_iter=10_000` budget documented under "Absorbed Fixed Effects with Survey Weights". Both surfaces emit `UserWarning` via `diff_diff.utils.warn_if_not_converged` when `max_iter` exhausts without reaching `tol` (the demean warning now carries the shared-engine label naming the affected variables rather than the estimator name). Silent return of the current iterate was classified as a silent failure under the Phase 2 audit and replaced with an explicit signal to match the logistic/Poisson IRLS pattern in `linalg.py`.
 - **Note:** Zero-total-weight groups (e.g. whole PSUs zeroed by JK1/BRR replicate weights, which reach Step 1 unmasked — `keep_mask` only drops always-treated units): a unit/period whose observations ALL carry zero weight has no identifying contribution and surfaces as `NaN` FE (key retained for the rank-condition membership check; matches the SpilloverDiD `_iterative_fe_subset` REGISTRY contract — never a silent finite `0.0`), and the shared demean engine's inert-row guard leaves those rows un-demeaned instead of NaN-poisoning the column. Before 3.7 the pandas loops divided 0/0 there: the covariate replicate path NaN-poisoned `y_dm`/`X_dm`, failed EVERY replicate refit inside `solve_ols(check_finite=True)`, and returned NaN SEs after a non-convergence warning storm; a main fit with zero-weight rows + covariates raised the same opaque `ValueError`. Both now produce finite results. `TwoStageDiD._mask_nan_ytilde`'s "non-finite imputed outcomes" `UserWarning` is suppressed (via `warn_nan=False`) ONLY inside the replicate-refit closures, where NaN FE for zeroed PSUs is expected mechanics — the main-fit warning is unchanged.
 - **Note:** `vcov_type` is permanently narrow to `{"hc1"}` per the Theorem 3 IF-based variance decomposition. Analytical-sandwich families `{classical, hc2, hc2_bm}` are rejected at `__init__` — the per-unit influence function aggregation has no equivalent single design matrix on which hat-matrix leverage or Bell-McCaffrey Satterthwaite DOF can be defined. `cluster=<col>` invokes per-cluster IF summation (Theorem 3 equation 7 conservative variance, `sigma_sq = (cluster_psi_sums**2).sum()` — plain CR1 with no Stata-style `(n-1)/(n-p)` finite-sample factor because the IF has no design-matrix `p` in the OLS sense); `cluster=None` (the default) routes the SAME Theorem 3 cluster-summed IF variance with `cluster_var = unit` (the unit column passed to `fit()`), so the summary renders `"CR1 cluster-robust at <unit>, G=<n_units>"` rather than the generic `"HC1"` label; `survey_design=` invokes TSL on the combined IF. Under bootstrap (`n_bootstrap > 0`) the analytical variance-family label is suppressed in `summary()` because `fit()` overwrites the reported SE/CI/p-value with bootstrap_results (mirrors the canonical `DiDResults` gate at `results.py:213-226`). `vcov_type='conley'` is deferred to the ImputationDiD Conley follow-up row in TODO.md.
@@ -1594,6 +1595,7 @@ where `W_it(h) = 1[K_it = h]` are lead indicators, estimated on `Omega_0` only.
 - [x] Step 3: Aggregate with researcher-chosen weights `w_it`
 - [x] Conservative clustered variance estimator (Theorem 3, Equation 7)
 - [x] Auxiliary model for treated residuals (unit-clustered Equation 8) with configurable partition (`aux_partition`)
+- [x] Leave-one-out finite-sample variance refinement (Supp. App. A.9, opt-in `leave_one_out`, default off)
 - [x] Supports unit FE, period FE, and time-varying covariates
 - [x] Refuses to estimate unidentified estimands (Proposition 5) — sets NaN with warning
 - [x] Pre-trend test uses only untreated observations (Test 1, Equation 9)
diff --git a/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md b/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md
index 5598b6ae..19ad8b7a 100644
--- a/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md
+++ b/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md
@@ -95,7 +95,7 @@ Partition `Ω₁ = ∪_g G_g`, impose `τ_it ≡ τ_g` within group, estimate:
                      Σ_i ( Σ_{t:it∈G_g} v_it )²
 ```
 - **Default partition (Stata `did_imputation`): groups by cohort × period** (when cohorts large); else by horizon relative to onset; single group = most conservative (minimal excess variance, Supp. App. A.8).
-- Finite-sample refinement: **leave-one-out** modification of `τ̃_it` (Supp. App. A.9 — GAP).
+- Finite-sample refinement: **leave-one-out** modification of `τ̃_it` (Supp. App. A.9). **✓ Implemented** as the opt-in `leave_one_out` parameter (default False); the efficient-rescale formula and derivation live in REGISTRY (see the A.9 provenance note under "Other notes").
 - **✓ Implemented (PR-B):** `_compute_auxiliary_residuals_treated` (`imputation.py`) computes the *unit-clustered* Eq. (8) above (per unit, the within-unit sums `a_{i,g}=Σ_t v_it` and `b_{i,g}=Σ_t v_it·τ̂*_it`, then `Σ_i a·b / Σ_i a²`). Non-target / unimputable (NaN `τ̂`, `v_it=0`) observations are excluded from the aggregation — exact for finite `τ̂` and NaN-safe. *Prior to the ImputationDiD methodology validation the code used the observation-level mean `Σ v·τ̂ / Σ v`, which coincides with Eq. (8) only under uniform within-group weights and diverges for survey-weighted / heterogeneity estimands (any partition) or the coarser `cohort` partition where a unit contributes multiple observations to a group.* The code docstring and the REGISTRY label now state the exact unit-clustered Eq. (8). Validated by white-box hand-calc + R `didimputation` parity in `tests/test_methodology_imputation.py`.
 - **No bootstrap is proposed** anywhere in the paper — SEs are analytical/cluster-robust. No explicit `G/(G−1)` or `(N−k)` DOF multiplier in the main text (the only finite-sample device is LOO).
 
@@ -149,7 +149,7 @@ Test  γ = 0  via heteroskedasticity- and cluster-robust Wald test.
 | Auxiliary partition `{G_g}` (eq. 8) | categorical | cohort × period | Within-group heterogeneity < across-group; coarser ⇒ more conservative; single group = max conservative |
 | Pre-trend lead count `k` (Test 1) | int | researcher choice | No automatic rule (fn 27); too-large `k` lowers power |
 | Step-1 weights (efficiency under known variances) | weights | unweighted OLS | `∝ σ_it^{−2}` (e.g. `∝ n_it` for aggregated cells) (fn 21) |
-| Leave-one-out variance refinement | bool | (library/supplement) | Supp. App. A.9 finite-sample improvement |
+| Leave-one-out variance refinement | bool | `leave_one_out=False` | Supp. App. A.9 finite-sample improvement — **✓ implemented** |
 | Anticipation horizon | int | 0 | Redefine event dates earlier by `k` periods |
 
 No regularization, bandwidth, factor-count, or cross-validation tuning (the estimator is unregularized OLS on untreated data).
@@ -177,7 +177,7 @@ No regularization, bandwidth, factor-count, or cross-validation tuning (the esti
 | **Supp. App. A.5** — heteroskedasticity/serial-correlation generalization | p. 3257, 3268 | Maps to any GLS-style variance extension |
 | **Supp. App. A.7 / Supp. Prop. A7** — low-level short-panel asymptotics | p. 3270–3272 | Conditions on researcher weights `w₁` |
 | **Supp. App. A.8** — minimal excess variance, single-group partition | p. 3272 | Justifies the conservative single-group default |
-| **Supp. App. A.9** — leave-one-out variance modification | p. 3256, 3272 | Finite-sample SE refinement |
+| **Supp. App. A.9** — leave-one-out variance modification | p. 3256, 3272 | Finite-sample SE refinement — **✓ implemented** (`leave_one_out`) |
 | **Supp. App. A.10** — efficient `v*_it` algorithm (iterative LS) | p. 3273 | The practical variance computation recipe |
 | **Supp. App. A.11** — full Monte Carlo study | p. 3268, 3282 | SD ratios 1.3–3.6×; coverage claims |
 | **Supp. App. B** — **ALL proofs** | p. 3262 | — |
@@ -188,4 +188,5 @@ No regularization, bandwidth, factor-count, or cross-validation tuning (the esti
 - The application (§5, Broda–Parker 2008 tax-rebate MPX) is context only, not implementation-relevant; preferred imputation estimate = first-month MPX $30.5 (3.4% of rebate); notional MPC 7.8–11.4%.
 - No contradictions were found between the three extraction passes or against the direct page reads; all numbered equations rendered cleanly (no `[UNREADABLE EQUATION]`).
 - **Resolved in PR-B:** the auxiliary `τ̃_g` aggregator now implements the exact unit-clustered Eq. (8) (`_compute_auxiliary_residuals_treated` + docstring + REGISTRY label), validated by white-box hand-calc + R `didimputation` parity (`tests/test_methodology_imputation.py`). The earlier obs-level mean was a valid conservative simplification but not the variance-minimizing form.
-- **Done in PR-B:** the **R `didimputation` parity fixture** and **`tests/test_methodology_imputation.py`** are on file. The `v*_it` weights (Supp. Prop. A3, absent from this PDF) were validated *empirically* against the reference: the exact two-way-FE projection `-A₀(A₀'A₀)⁻¹A₁'w` now matches `didimputation` (observed ~1e-10; tests assert SE `abs=1e-7`) — the prior FE-only closed form `-(w_i/n0_i + w_t/n0_t − w/N₀)` was a balanced-panel approximation that biased the SE ~27% on staggered (unbalanced-Ω₀) designs and was corrected in PR-B. The leave-one-out variance refinement (Supp. App. A.9) remains unimplemented (future work).
+- **Done in PR-B:** the **R `didimputation` parity fixture** and **`tests/test_methodology_imputation.py`** are on file. The `v*_it` weights (Supp. Prop. A3, absent from this PDF) were validated *empirically* against the reference: the exact two-way-FE projection `-A₀(A₀'A₀)⁻¹A₁'w` now matches `didimputation` (observed ~1e-10; tests assert SE `abs=1e-7`) — the prior FE-only closed form `-(w_i/n0_i + w_t/n0_t − w/N₀)` was a balanced-panel approximation that biased the SE ~27% on staggered (unbalanced-Ω₀) designs and was corrected in PR-B. The leave-one-out variance refinement (Supp. App. A.9) is now **implemented** as the opt-in `leave_one_out` parameter (default False, preserving R `didimputation` parity) — see the A.9 provenance note below.
+- **Supp. App. A.9 provenance (leave-one-out, added post-review):** A.9 was a GAP in this main-article review (the REStud Supplementary Material is not in the reviewed PDF). It was subsequently sourced from the arXiv preprint supplement (**arXiv:2108.12419v5, Appendix A.9 "Leave-Out Conservative Variance Estimation"**) to implement the opt-in `leave_one_out` parameter. The **REStud Supplementary Material is the canonical version**; per this doc's version-pin scope note, the efficient-rescale formula `ε̃^LO = ε̃/(1 − v_ig²/Σ_j v_jg²)`, its exact ψ-level equivalence to the direct leave-one-out `τ̃_it^LO`, Prop. A8 (unbiased for an upper bound), and the footnote-51 single-positive-weight-unit edge are recorded in `docs/methodology/REGISTRY.md` `## ImputationDiD` (NOT transcribed here). Validated by the exact ψ-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`); the authors' Stata `did_imputation` ships the same option (no CI-runnable anchor).
diff --git a/tests/test_methodology_imputation.py b/tests/test_methodology_imputation.py
index f70081b7..49f308ec 100644
--- a/tests/test_methodology_imputation.py
+++ b/tests/test_methodology_imputation.py
@@ -586,6 +586,350 @@ def test_public_se_regression_pin_cohort_partition(self) -> None:
 _EQ8_COHORT_SE_PIN = 0.042000264835
 
 
+# =============================================================================
+# Supplementary Appendix A.9 — leave-one-out conservative variance
+# =============================================================================
+
+
+def _synthetic_group(a: np.ndarray, T: np.ndarray):
+    """Build (per_unit, per_group) for a single group with per-unit weight sums
+    ``a_i`` and unit effects ``T_i``, matching the layout that
+    ``_compute_auxiliary_residuals_treated`` passes to ``_leave_one_out_factor``:
+    ``per_unit`` MultiIndex (g, u) with column ``a``; ``per_group`` index g with
+    column ``den = sum_i a_i^2``.
+    """
+    idx = pd.MultiIndex.from_arrays([[0] * len(a), list(range(len(a)))], names=["g", "u"])
+    per_unit = pd.DataFrame({"a": a, "b": a * T}, index=idx)
+    per_group = pd.DataFrame(
+        {"num": [float((a**2 * T).sum())], "den": [float((a**2).sum())]},
+        index=pd.Index([0], name="g"),
+    )
+    return per_unit, per_group
+
+
+class TestB2024AppendixA9LeaveOneOut:
+    """BJS 2024 Supplementary Appendix A.9 leave-one-out conservative variance.
+
+    The efficient rescale ``eps_tilde^LO = eps_tilde / (1 - v_ig^2/sum_j v_jg^2)``
+    reproduces the direct leave-one-out aggregate ``tau_tilde_it^LO`` exactly at
+    the per-unit cluster sum ``psi_i`` (App. A.9). Source: arXiv 2108.12419v5
+    App. A.9 (the REStud Supplementary Material is canonical); see REGISTRY.
+    """
+
+    _COMMON = dict(outcome="outcome", unit="unit", time="time", first_treat="first_treat")
+
+    def test_loo_factor_reproduces_direct_leave_one_out(self):
+        """PAPER-FIDELITY GATE: the rescaled per-unit cluster sum psi_i equals the
+        direct-LOO psi_i (recomputed with the paper's tau_tilde_it^LO), exactly."""
+        a = np.array([1.3, 0.7, 2.1, 0.9, 1.6])
+        T = np.array([2.0, 3.0, 1.5, 2.5, 1.9])
+        per_unit, per_group = _synthetic_group(a, T)
+        factor, n_single = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+        assert n_single == 0
+
+        N = float((a**2 * T).sum())
+        D = float((a**2).sum())
+        tau_g = N / D
+        # one obs per unit: eps_i = T_i - tau_tilde_g, psi_i = a_i * eps_i * factor_i
+        psi_rescale = a * (T - tau_g) * factor.to_numpy()
+        # direct LOO: tau_tilde_i^LO = (N - a_i^2 T_i)/(D - a_i^2)
+        psi_direct = np.array(
+            [a[i] * (T[i] - (N - a[i] ** 2 * T[i]) / (D - a[i] ** 2)) for i in range(len(a))]
+        )
+        np.testing.assert_allclose(psi_rescale, psi_direct, rtol=0, atol=1e-12)
+        # and the closed-form factor 1/(1 - a^2/D)
+        np.testing.assert_allclose(factor.to_numpy(), 1.0 / (1.0 - a**2 / D), rtol=0, atol=1e-14)
+
+    def test_loo_factor_equal_weight_is_k_over_k_minus_1(self):
+        """Equal-weight K-unit group -> per-unit factor = K/(K-1)."""
+        for K in (2, 3, 5, 10, 50):
+            a = np.ones(K)
+            T = 2.0 + 0.1 * np.arange(K)
+            per_unit, per_group = _synthetic_group(a, T)
+            factor, n_single = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+            assert n_single == 0
+            np.testing.assert_allclose(
+                factor.to_numpy(), np.full(K, K / (K - 1)), rtol=0, atol=1e-12
+            )
+
+    def test_loo_factor_unit_dominated_group_keeps_large_finite_factor(self):
+        """A >=2-unit group so dominated that `D - a_i^2` cancels to 0 in float64
+        (a=[1.0, 1e-8] -> D==1.0 exactly) must keep its large FINITE factor via
+        the exact sum-of-others denominator, NOT silently fall back to non-LOO."""
+        a = np.array([1.0, 1e-8])  # D = 1 + 1e-16 rounds to 1.0; D - a_0^2 == 0.0
+        T = np.array([2.0, 3.0])
+        per_unit, per_group = _synthetic_group(a, T)
+        factor, n_single = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+        assert n_single == 0  # 2 positive-weight units -> not a singleton
+        f0 = factor.to_numpy()[0]
+        assert np.isfinite(f0) and f0 > 1e10  # ~1/1e-16, NOT the fallback 1.0
+
+    def test_loo_factor_positive_near_cancellation_is_exact(self):
+        """POSITIVE near-cancellation (a=[1.0, 1e-6] -> D - a_0^2 stays > 0 but
+        loses ~4 digits) must use the exact sum-of-others denominator, not the
+        finite-but-wrong `D - a_i^2`. Exact factor = D / a_1^2."""
+        a = np.array([1.0, 1e-6])
+        T = np.array([2.0, 3.0])
+        per_unit, per_group = _synthetic_group(a, T)
+        factor, _ = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+        D = float((a**2).sum())
+        exact = D / (a[1] ** 2)  # = D / sum_{j!=0} a_j^2
+        # rel=1e-9 fails against the lossy `D - a_0^2` (off ~2e-4) but passes exact.
+        assert factor.to_numpy()[0] == pytest.approx(exact, rel=1e-9)
+
+    def test_loo_factor_effective_singleton_falls_back(self):
+        """A >=2-ROW group with only one positive-weight unit (a=[1.5, 0.0]) is an
+        effective singleton (fn. 51) -> factor 1.0 and counted as single."""
+        a = np.array([1.5, 0.0])
+        T = np.array([2.0, 3.0])
+        per_unit, per_group = _synthetic_group(a, T)
+        factor, n_single = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+        assert n_single == 1
+        np.testing.assert_allclose(factor.to_numpy(), [1.0, 1.0], rtol=0, atol=1e-14)
+
+    def test_loo_factor_single_unit_group_falls_back(self):
+        """Single positive-weight unit (App. A.9 fn. 51) -> factor 1.0, counted."""
+        a = np.array([1.5])
+        T = np.array([2.0])
+        per_unit, per_group = _synthetic_group(a, T)
+        factor, n_single = ImputationDiD._leave_one_out_factor(per_unit, per_group)
+        assert n_single == 1
+        np.testing.assert_allclose(factor.to_numpy(), [1.0], rtol=0, atol=1e-14)
+
+    def test_loo_se_geq_nonloo_at_unit_clustering(self):
+        """Prop. A8 direction: at the default unit clustering, LOO SE >= non-LOO;
+        strict on a multi-unit panel. ATT is unchanged."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 11)
+        panel = _make_staggered_panel(
+            rng, cohorts=[3, 4], n_per_cohort=40, n_periods=6, tau_constant=1.0
+        )
+        r0 = ImputationDiD(leave_one_out=False).fit(panel, **self._COMMON)
+        r1 = ImputationDiD(leave_one_out=True).fit(panel, **self._COMMON)
+        np.testing.assert_allclose(r1.overall_att, r0.overall_att, rtol=0, atol=1e-12)
+        assert r1.overall_se > r0.overall_se
+
+    def test_loo_false_is_byte_identical_to_default(self):
+        """leave_one_out=False changes nothing on the default path."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 12)
+        panel = _make_staggered_panel(rng, cohorts=[3], n_per_cohort=50, n_periods=5)
+        se_default = ImputationDiD().fit(panel, **self._COMMON).overall_se
+        se_loo_false = ImputationDiD(leave_one_out=False).fit(panel, **self._COMMON).overall_se
+        assert se_default == se_loo_false
+
+    def test_loo_single_unit_group_warns_and_returns_finite(self):
+        """A singleton treated cohort makes each cohort x horizon group a single
+        unit -> UserWarning + non-LOO fallback, finite SE."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 13)
+        rows: List[Dict[str, Any]] = []
+        uid = 0
+        for _ in range(8):  # never-treated controls
+            c_i = rng.standard_normal()
+            for t in range(1, 6):
+                rows.append(
+                    dict(
+                        unit=uid,
+                        time=t,
+                        first_treat=0,
+                        outcome=c_i + 0.5 * t + 0.1 * rng.standard_normal(),
+                    )
+                )
+            uid += 1
+        c_i = rng.standard_normal()  # ONE treated unit, cohort 3
+        for t in range(1, 6):
+            rows.append(
+                dict(
+                    unit=uid,
+                    time=t,
+                    first_treat=3,
+                    outcome=c_i + 0.5 * t + (1.0 if t >= 3 else 0.0) + 0.1 * rng.standard_normal(),
+                )
+            )
+        panel = pd.DataFrame(rows)
+        with pytest.warns(UserWarning, match="single positive-weight unit"):
+            res = ImputationDiD(leave_one_out=True).fit(panel, **self._COMMON)
+        assert np.isfinite(res.overall_se)
+
+    def test_loo_param_validation_and_roundtrip(self):
+        """leave_one_out is in get/set_params; a non-bool is rejected (TypeError)
+        in __init__ AND at fit-time (closing the naive-setattr set_params bypass)."""
+        assert ImputationDiD().get_params()["leave_one_out"] is False
+        assert ImputationDiD(leave_one_out=True).get_params()["leave_one_out"] is True
+        with pytest.raises(TypeError, match="leave_one_out must be a bool"):
+            ImputationDiD(leave_one_out="yes")  # type: ignore[arg-type]
+        # set_params is a naive setattr; the fit-time re-check must catch it
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 14)
+        panel = _make_staggered_panel(rng, cohorts=[3], n_per_cohort=20, n_periods=5)
+        est = ImputationDiD()
+        est.set_params(leave_one_out="yes")  # type: ignore[arg-type]
+        with pytest.raises(TypeError, match="leave_one_out must be a bool"):
+            est.fit(panel, **self._COMMON)
+
+    def test_loo_fit_is_idempotent_on_config(self):
+        """Repeated fits with leave_one_out=True give identical SE (no config mutation)."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 15)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=25, n_periods=6)
+        est = ImputationDiD(leave_one_out=True)
+        se_a = est.fit(panel, **self._COMMON).overall_se
+        se_b = est.fit(panel, **self._COMMON).overall_se
+        assert se_a == se_b
+        assert est.leave_one_out is True
+
+    def test_loo_composes_with_cluster_and_bootstrap(self):
+        """LOO applies the same residual rescale under a coarser cluster= and
+        under the multiplier bootstrap: finite SE (NOT asserting the >= direction
+        away from unit clustering)."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 16)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=30, n_periods=6)
+        panel["state"] = panel["unit"] % 5  # coarser cluster
+        se_cluster = (
+            ImputationDiD(leave_one_out=True, cluster="state").fit(panel, **self._COMMON).overall_se
+        )
+        assert np.isfinite(se_cluster)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            se_boot = (
+                ImputationDiD(leave_one_out=True, n_bootstrap=49, seed=7)
+                .fit(panel, **self._COMMON)
+                .overall_se
+            )
+        assert np.isfinite(se_boot)
+
+    def test_loo_replicate_weight_survey_raises(self):
+        """Replicate-weight variance bypasses the conservative IF path where LOO
+        lives, so leave_one_out=True must fail closed (not silently no-op)."""
+        from diff_diff.survey import SurveyDesign
+
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 17)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=20, n_periods=5)
+        panel["weight"] = 1.0
+        panel["rw1"] = 1.0
+        panel["rw2"] = 1.0
+        design = SurveyDesign(
+            weights="weight",
+            replicate_weights=["rw1", "rw2"],
+            replicate_method="JK1",
+            weight_type="pweight",
+        )
+        with pytest.raises(NotImplementedError, match="leave_one_out=True"):
+            ImputationDiD(leave_one_out=True).fit(panel, survey_design=design, **self._COMMON)
+
+    def test_loo_recorded_in_results_metadata(self):
+        """The result object is self-describing: leave_one_out is persisted on the
+        result, in to_dict(), and surfaced in summary() (it changes reported SEs)."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 18)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=20, n_periods=5)
+        r1 = ImputationDiD(leave_one_out=True).fit(panel, **self._COMMON)
+        r0 = ImputationDiD(leave_one_out=False).fit(panel, **self._COMMON)
+        assert r1.leave_one_out is True and r0.leave_one_out is False
+        assert r1.to_dict()["leave_one_out"] is True
+        assert r0.to_dict()["leave_one_out"] is False
+        assert "Leave-one-out" in r1.summary()
+        assert "Leave-one-out" not in r0.summary()
+
+    @staticmethod
+    def _assert_inference_consistent(effects):
+        """Every cell is either a genuine estimate (se > 0 -> effect/t/p/CI all
+        finite) or a degenerate/reference cell (se == 0 or NaN -> NaN t/p, per the
+        safe_inference contract), and at least one cell is a genuine estimate."""
+        assert effects is not None and len(effects) > 0
+        n_finite = 0
+        for eff in effects.values():
+            se = eff["se"]
+            if np.isfinite(se) and se > 0:
+                assert np.isfinite(eff["effect"])
+                assert np.isfinite(eff["t_stat"])
+                assert np.isfinite(eff["p_value"])
+                lo, hi = eff["conf_int"]
+                assert np.isfinite(lo) and np.isfinite(hi)
+                n_finite += 1
+            else:
+                # zero or non-finite SE -> undefined t-stat / p-value (NaN)
+                assert np.isnan(eff["t_stat"]) and np.isnan(eff["p_value"])
+        assert n_finite > 0
+
+    def test_loo_aggregate_all_analytical(self):
+        """aggregate="all": LOO routes through the event-study AND group
+        aggregators (not just overall). ATT unchanged vs non-LOO, overall LOO SE
+        >= non-LOO, and both aggregation surfaces have consistent finite inference."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 19)
+        panel = _make_staggered_panel(
+            rng, cohorts=[3, 4], n_per_cohort=40, n_periods=6, tau_constant=1.0
+        )
+        r0 = ImputationDiD(leave_one_out=False).fit(panel, aggregate="all", **self._COMMON)
+        r1 = ImputationDiD(leave_one_out=True).fit(panel, aggregate="all", **self._COMMON)
+        np.testing.assert_allclose(r1.overall_att, r0.overall_att, rtol=0, atol=1e-12)
+        assert r1.overall_se > r0.overall_se
+        self._assert_inference_consistent(r1.event_study_effects)
+        self._assert_inference_consistent(r1.group_effects)
+
+    def test_loo_aggregate_all_bootstrap(self):
+        """aggregate="all" under the multiplier bootstrap: event-study and group
+        inference fields are populated and consistent."""
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 20)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=30, n_periods=6)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            r = ImputationDiD(leave_one_out=True, n_bootstrap=99, seed=3).fit(
+                panel, aggregate="all", **self._COMMON
+            )
+        assert np.isfinite(r.overall_se)
+        self._assert_inference_consistent(r.event_study_effects)
+        self._assert_inference_consistent(r.group_effects)
+
+    def test_loo_aggregate_all_analytical_survey(self):
+        """aggregate="all" with an analytical survey design (weights + PSU): the
+        LOO rescale composes through the survey TSL variance with consistent
+        inference on every surface."""
+        from diff_diff.survey import SurveyDesign
+
+        rng = np.random.default_rng(_BASE_SEED_EQ8 + 21)
+        panel = _make_staggered_panel(rng, cohorts=[3, 4], n_per_cohort=40, n_periods=6)
+        panel["weight"] = 1.0 + (panel["unit"] % 3) * 0.1
+        panel["psu"] = panel["unit"] % 8
+        design = SurveyDesign(weights="weight", psu="psu")
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            r = ImputationDiD(leave_one_out=True).fit(
+                panel, aggregate="all", survey_design=design, **self._COMMON
+            )
+        assert np.isfinite(r.overall_se)
+        self._assert_inference_consistent(r.event_study_effects)
+        self._assert_inference_consistent(r.group_effects)
+
+    @pytest.mark.slow
+    def test_loo_coverage_geq_nominal(self):
+        """MC coverage: on an overfit-prone fine partition, LOO coverage is >=
+        non-LOO coverage and near/above nominal (LOO removes the downward bias)."""
+        true_att = 1.0
+        n_rep = 200
+        cov0 = cov1 = 0
+        for rep in range(n_rep):
+            rng = np.random.default_rng(20000 + rep)
+            panel = _make_staggered_panel(
+                rng,
+                cohorts=[3, 4, 5],
+                n_per_cohort=8,
+                n_periods=6,
+                tau_constant=true_att,
+                sigma=1.0,
+            )
+            r0 = ImputationDiD(leave_one_out=False).fit(panel, **self._COMMON)
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                r1 = ImputationDiD(leave_one_out=True).fit(panel, **self._COMMON)
+            for r, hit in ((r0, "cov0"), (r1, "cov1")):
+                lo = r.overall_att - 1.96 * r.overall_se
+                hi = r.overall_att + 1.96 * r.overall_se
+                if lo <= true_att <= hi:
+                    if hit == "cov0":
+                        cov0 += 1
+                    else:
+                        cov1 += 1
+        cov0f, cov1f = cov0 / n_rep, cov1 / n_rep
+        assert cov1f >= cov0f - 0.02  # LOO no worse than non-LOO
+        assert cov1f >= 0.90  # near/above nominal 95%
+
+
 # =============================================================================
 # Proposition 5 — non-identification without never-treated units
 # =============================================================================