igerber · igerber · Jul 6, 2026 · Jul 6, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -22,6 +22,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - **StackedDiD** event-study coefficient point estimates; **PlaceboTests** leave-one-out
     `t_stat`/`p_value`; **DIDHAD** QUG order-statistic test (`t_stat`/`p_value`, previously
     unasserted against the golden). All at ~1e-10.
+  - **estimatr 2SLS intercept SE** (HAD mass-point path): the `HC1`/`CR1` intercept SE now
+    matches `estimatr::iv_robust` `se_intercept` to ~1e-15 (the 2×2 sandwich already computed
+    `V[0,0]`; an opt-in `return_intercept_se=` hook on the private `_fit_mass_point_2sls`
+    surfaces it — the default 3-tuple return is byte-unchanged, no public-API change). The
+    `classical` intercept is excluded (same documented `O(1/n)` projection/DOF deviation as the slope).
+  - **CR2 Satterthwaite DOF** via CI-inversion (no new result fields): `TwoWayFixedEffects(hc2_bm)`
+    `dof_hc2_bm` and `MultiPeriodDiD(hc2_bm)` per-period `dof_per_coef` are pinned by reconstructing
+    `conf_int` from the golden DOF — with ATT+SE already locked, the CI matches iff the DOF matches.
+  - **PlaceboTests** leave-one-out `df` (recovered from the public `leave_one_out_effects` count);
+    **DIDHAD** Yatchew linearity test `p_value` and both `sigma2` components (through the documented
+    `N/(N-1)` sample-vs-population convention shift). All at ~1e-10.
+  - **`fixest` cluster-robust SE band**: the DiD/TWFE cluster-at-unit SE is pinned within the
+    documented ~0.25% fixest-CR1 small-sample DOF-convention band (guards an unintended SE-formula
+    change; the machine-precision hetero/cluster lock is deferred — needs an unbalanced-DGP golden).
 
 ### Added
 - **`ContinuousDiD` lowest-dose-as-control** (`control_group="lowest_dose"`, CGBS 2024 Remark 3.1) for

diff --git a/TODO.md b/TODO.md
@@ -66,7 +66,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
-| SE-audit CI-lock follow-up ("fiddly bits" deferred from the G2+C1-C3/C5 coverage batch). The `fixest_did_twfe_golden.json` stores `cluster_unit` ATT+SE for DiD and TWFE but only the ATT is asserted (the cluster SE carries the documented CR1 DOF-convention difference vs fixest; the hetero SE collapses to iid on these balanced 2-group designs). Also deferred, each needing plumbing or convention work: surface the 2SLS intercept SE for estimatr `se_intercept` parity (C4); surface StackedDiD intercept SEs (`se_cr1/cr2_intercept`, C1); add a `df` field + `boundary_gap` computation for PlaceboTests (C3); assert Yatchew `p`/`sigma2_lin/diff` through the N/(N-1) convention shift + event-map (C5); assert `dof_hc2_bm`/`dof_per_coef` via CI-inversion (C2). Plus the tolerance-tightenings C6-C8. | `tests/test_fixest_did_twfe_parity.py`, `tests/test_estimatr_iv_robust_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py`, `tests/test_did_had_parity.py` | SE-audit | Mid | Low |
+| SE-audit CI-lock — remaining "fiddly bits" after the second coverage batch landed (that batch pinned C2 `dof_hc2_bm`/`dof_per_coef` via CI-inversion, C3 LOO `df`, C4 estimatr HC1/CR1 intercept SE, C5 Yatchew `p`/`sigma2_lin`/`sigma2_diff`, and the G2 fixest cluster-SE band). Still deferred, each needing a golden regeneration, new computation, or a documented-deviation call: **(a) G2 machine-precision hetero/cluster SE** — needs an unbalanced/heteroskedastic-DGP regen of `fixest_did_twfe_golden.json` (hetero collapses to iid on the current balanced design; the cluster SE is the documented ~0.25% fixest-CR1 DOF-convention deviation, currently band-pinned only); **(b) PlaceboTests `boundary_gap`** — a permutation randomization-inference margin NOT computed anywhere in code (a new feature + result field, not a coverage lock); **(c) StackedDiD intercept SEs** (`se_cr1/cr2_intercept`, C1) — MEASURED to diverge ~0.3% from R: a nuisance-parameter reference-cell/parameterization gap, NOT machine-precision lockable (the event-study interaction SEs already match ~2e-13; surfacing it would add an unasserted, R-divergent public field); **(d) estimatr `classical` intercept SE** — same documented `O(1/n)` projection/DOF deviation as the slope (reference-only, excluded from parity). Plus the tolerance-tightenings C6-C8. | `benchmarks/R/generate_fixest_did_twfe_golden.R`, `tests/test_fixest_did_twfe_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py` | SE-audit | Mid | Low |
 | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low |
 | `ImputationDiD` covariate-path variance lacks a dedicated parity anchor — only the no-covariate staggered panel is R-parity'd, though the covariate path shares the same validated projection code. Add a small dense-design **hand-calc** for the covariate projection (no external tooling), or a covariate (time-varying X) R `didimputation` golden asserting overall/ES SE parity (the golden variant needs local R). | `tests/test_methodology_imputation.py`, `benchmarks/R/generate_didimputation_golden.R` | imputation-validation | Mid | Low |
 | Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low |

diff --git a/diff_diff/had.py b/diff_diff/had.py
@@ -66,7 +66,7 @@
 
 import warnings
 from dataclasses import dataclass
-from typing import Any, Dict, List, Optional, Tuple, Union
+from typing import Any, Dict, List, Literal, Optional, Tuple, Union, overload
 
 import numpy as np
 import pandas as pd
@@ -2326,6 +2326,34 @@ def _sup_t_multiplier_bootstrap(
 # =============================================================================
 
 
+@overload
+def _fit_mass_point_2sls(
+    d: np.ndarray,
+    dy: np.ndarray,
+    d_lower: float,
+    cluster: Optional[np.ndarray],
+    vcov_type: str,
+    *,
+    weights: Optional[np.ndarray] = ...,
+    return_influence: bool = ...,
+    return_intercept_se: Literal[False] = ...,
+) -> Tuple[float, float, Optional[np.ndarray]]: ...
+
+
+@overload
+def _fit_mass_point_2sls(
+    d: np.ndarray,
+    dy: np.ndarray,
+    d_lower: float,
+    cluster: Optional[np.ndarray],
+    vcov_type: str,
+    *,
+    weights: Optional[np.ndarray] = ...,
+    return_influence: bool = ...,
+    return_intercept_se: Literal[True],
+) -> Tuple[float, float, Optional[np.ndarray], float]: ...
+
+
 def _fit_mass_point_2sls(
     d: np.ndarray,
     dy: np.ndarray,
@@ -2335,7 +2363,11 @@ def _fit_mass_point_2sls(
     *,
     weights: Optional[np.ndarray] = None,
     return_influence: bool = False,
-) -> Tuple[float, float, Optional[np.ndarray]]:
+    return_intercept_se: bool = False,
+) -> Union[
+    Tuple[float, float, Optional[np.ndarray]],
+    Tuple[float, float, Optional[np.ndarray], float],
+]:
     """Wald-IV point estimate and structural-residual 2SLS sandwich SE.
 
     The just-identified binary instrument ``Z_g = 1{D_{g,2} > d_lower}``
@@ -2390,6 +2422,13 @@ def _fit_mass_point_2sls(
         ``compute_survey_if_variance(IF, trivial)`` ≈ ``V_HC1[1, 1]``
         at ``atol=1e-10`` (PR #359 convention; see the "IF scale
         convention" section of the Phase 4.5 B plan for derivation).
+    return_intercept_se : bool
+        When True, the returned tuple gains a 4th element: the intercept
+        standard error ``sqrt(V[0, 0])`` (``V[1, 1]`` is the slope
+        variance that drives ``se_beta``). ``NaN`` whenever ``se_beta``
+        is ``NaN`` (early-exit / singular paths). Test-only surfacing hook
+        for estimatr-parity coverage; production callers leave it False and
+        the return is the unchanged 3-tuple.
 
     Returns
     -------
@@ -2398,7 +2437,9 @@ def _fit_mass_point_2sls(
         ``return_influence=True``, else ``None``. NaN for SE when the
         dose-gap vanishes (``Dbar_{Z=1} == Dbar_{Z=0}``) or the
         sandwich is singular; in those cases ``psi`` is returned as a
-        length-``n`` zero array when ``return_influence=True``.
+        length-``n`` zero array when ``return_influence=True``. When
+        ``return_intercept_se=True`` a 4th element ``se_intercept`` is
+        appended.
     """
     d = np.asarray(d, dtype=np.float64)
     dy = np.asarray(dy, dtype=np.float64)
@@ -2443,8 +2484,25 @@ def _fit_mass_point_2sls(
 
     _null_psi = np.zeros(n, dtype=np.float64) if return_influence else None
 
+    def _pack(
+        beta: float,
+        se_b: float,
+        psi_out: Optional[np.ndarray],
+        se_int: float = float("nan"),
+    ) -> Union[
+        Tuple[float, float, Optional[np.ndarray]],
+        Tuple[float, float, Optional[np.ndarray], float],
+    ]:
+        # Centralize the optional 4th (intercept-SE) element so every one of
+        # the function's early-exit and success returns keeps a consistent
+        # arity. The @overload stubs give callers the precise per-flag type;
+        # the True branch appends se_int for the test-only surfacing hook.
+        if return_intercept_se:
+            return (beta, se_b, psi_out, se_int)
+        return (beta, se_b, psi_out)
+
     if n_above == 0 or n_at_or_below == 0:
-        return float("nan"), float("nan"), _null_psi
+        return _pack(float("nan"), float("nan"), _null_psi)
 
     # Point estimate: weighted Wald-IV ratio (reduces to unweighted at w=1).
     if weighted:
@@ -2453,7 +2511,7 @@ def _fit_mass_point_2sls(
         w_Z1 = float(w_arr[Z1_idx].sum())
         w_Z0 = float(w_arr[Z0_idx].sum())
         if w_Z1 <= 0.0 or w_Z0 <= 0.0:
-            return float("nan"), float("nan"), _null_psi
+            return _pack(float("nan"), float("nan"), _null_psi)
         dose_gap = float(
             (w_arr[Z1_idx] * d[Z1_idx]).sum() / w_Z1 - (w_arr[Z0_idx] * d[Z0_idx]).sum() / w_Z0
         )
@@ -2470,7 +2528,7 @@ def _fit_mass_point_2sls(
 
     if abs(dose_gap) < 1e-12 * max(1.0, abs(d_bar)):
         # No dose variation around d_lower -> beta undefined.
-        return float("nan"), float("nan"), _null_psi
+        return _pack(float("nan"), float("nan"), _null_psi)
 
     # dy_gap / dy_bar were computed inside the weighted/unweighted blocks above.
     beta_hat = float(dy_gap / dose_gap)
@@ -2495,7 +2553,7 @@ def _fit_mass_point_2sls(
     try:
         ZtWX_inv = np.linalg.inv(ZtWX)
     except np.linalg.LinAlgError:
-        return beta_hat, float("nan"), _null_psi
+        return _pack(beta_hat, float("nan"), _null_psi)
 
     vcov_type = vcov_type.lower()
     if vcov_type in _MASS_POINT_VCOV_UNSUPPORTED:
@@ -2525,25 +2583,25 @@ def _fit_mass_point_2sls(
         n_clusters = len(clusters_unique)
         if n_clusters < 2:
             # Cluster-robust SE undefined with a single cluster.
-            return beta_hat, float("nan"), _null_psi
+            return _pack(beta_hat, float("nan"), _null_psi)
         Omega *= (n_clusters / (n_clusters - 1)) * ((n - 1) / (n - k))
     elif vcov_type == "classical":
         if weighted:
             dof = w_sum - k
             if dof <= 0:
-                return beta_hat, float("nan"), _null_psi
+                return _pack(beta_hat, float("nan"), _null_psi)
             sigma2 = float((w_arr * w_arr * u * u).sum()) / dof
             Omega = sigma2 * (Zd.T @ ((w_arr * w_arr)[:, None] * Zd))
         else:
             dof = n - k
             if dof <= 0:
-                return beta_hat, float("nan"), _null_psi
+                return _pack(beta_hat, float("nan"), _null_psi)
             sigma2 = float((u * u).sum()) / dof
             Omega = sigma2 * (Zd.T @ Zd)
     elif vcov_type == "hc1":
         dof = n - k
         if dof <= 0:
-            return beta_hat, float("nan"), _null_psi
+            return _pack(beta_hat, float("nan"), _null_psi)
         if weighted:
             # Pweight HC1: meat = Z' diag(w² u²) Z (Wooldridge 2010 Eq 12.37,
             # matches linalg.py:1141 convention and estimatr::iv_robust
@@ -2561,11 +2619,19 @@ def _fit_mass_point_2sls(
     V = ZtWX_inv @ Omega @ ZtWX_inv.T
     var_beta = float(V[1, 1])
     if not np.isfinite(var_beta) or var_beta < 0:
-        return beta_hat, float("nan"), _null_psi
+        return _pack(beta_hat, float("nan"), _null_psi)
     se_beta = float(np.sqrt(var_beta))
+    # Intercept SE = sqrt(V[0, 0]); surfaced only when return_intercept_se
+    # (see _pack). NaN if the intercept variance is non-finite/negative.
+    var_intercept = float(V[0, 0])
+    se_intercept = (
+        float(np.sqrt(var_intercept))
+        if np.isfinite(var_intercept) and var_intercept >= 0
+        else float("nan")
+    )
 
     if not return_influence:
-        return beta_hat, se_beta, None
+        return _pack(beta_hat, se_beta, None, se_intercept)
 
     # Per-unit influence function on β̂-scale, scaled so that
     # compute_survey_if_variance(psi, trivial_resolved) ≈ V_HC1[1, 1]
@@ -2583,7 +2649,7 @@ def _fit_mass_point_2sls(
         psi0 = (Zd @ bread_row) * u
     dof_psi = max(n - k, 1)
     psi = psi0 * np.sqrt((n - 1) / dof_psi)
-    return beta_hat, se_beta, psi
+    return _pack(beta_hat, se_beta, psi, se_intercept)
 
 
 # =============================================================================

diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
@@ -3167,6 +3167,7 @@ Under `survey=SurveyDesign(weights, strata, psu, fpc)`, the variance composes vi
 - **CR1 pweight-cluster meat**: for each cluster c, `s_c = Z'_c·(w·u)_c`; `Ω_CR1 = (G/(G-1))·((n-1)/(n-k))·Σ_c s_c s_c'` (`w¹` inside cluster score). Bit-exact with `estimatr::iv_robust(..., weights=, clusters=, se_type="stata")` at `atol=1e-10`.
 - **Classical**: sandwich form `Ω_cl = σ²·Z'·diag(w²)·Z` with `σ² = Σw²u²/(Σw-k)`. Deviates from `estimatr` classical (projection-form + `n-k` DOF) by `O(1/n)` at non-uniform weights; unweighted path is bit-exact by equivalence. Skipped in cross-language parity tests.
 - **Per-unit IF on β̂-scale** (for Binder-TSL survey composition): `psi_g = [(Z'WX)^{-1} · z_g · w_g · u_g][1] · sqrt((n-1)/(n-k))`. The scaling factor absorbs DOF / small-sample differences so `compute_survey_if_variance(psi, trivial_resolved) ≈ V_HC1[1,1]` at `atol=1e-10` (mirrors PR #359 convention; asserted by `TestIFScaleInvariant` and bit-exact against estimatr HC1 on 4 DGPs). Fixture: `benchmarks/R/generate_estimatr_iv_robust_golden.R` → `benchmarks/data/estimatr_iv_robust_golden.json`.
+- **Intercept SE (`return_intercept_se=True`)**: the 2×2 sandwich `V` already carries the intercept variance `V[0,0]`; the opt-in hook surfaces `sqrt(V[0,0])` (default off — the production 3-tuple return is byte-unchanged, so no public-API change). HC1 and CR1 intercept SEs are bit-exact with `estimatr` `se_intercept` at `atol=1e-10` (`tests/test_estimatr_iv_robust_parity.py`); the **classical** intercept carries the same `O(1/n)` projection/DOF deviation as the slope (bullet above) and is likewise excluded from the parity lock.
 
 *Event-study survey composition (Phase 4.5 B):* The per-horizon loop in `_fit_event_study` threads `weights_unit_full` + `resolved_survey_unit_full` through to both `_fit_continuous` and `_fit_mass_point_2sls` (the latter with `return_influence=True` under weighted fits). The returned IF matrix `Psi ∈ R^{G × H}` has a shared construction contract across paths — each column on the β̂-scale, such that `compute_survey_if_variance(Psi[:, e], resolved) ≈ V_β[e]`. Per-horizon analytical variance uses Binder-TSL via `compute_survey_if_variance` (survey= path) or the weighted-robust HC1 sandwich (weights= shortcut). `survey_metadata`, `variance_formula` (`"survey_binder_tsl"` / `"survey_binder_tsl_2sls"` / `"pweight"` / `"pweight_2sls"`), and `effective_dose_mean` populate identically to the static path. Pre-PR numerical output is preserved bit-exactly on the unweighted path when `cband=False` (stability invariant; Phase 2b convention unchanged for unweighted fits).
 

diff --git a/tests/test_did_had_parity.py b/tests/test_did_had_parity.py
@@ -431,6 +431,17 @@ def test_yatchew_t_stat_parity(
         r_yatchew_t = _as_list(r_result["yatchew_t"])
         r_yatchew_n = _as_list(r_result["yatchew_n"])
         r_event_ids = _as_list(r_result["event_id"])
+        # SE-audit C5: also lock the yatchew p-value and both sigma^2 columns
+        # (previously only T_hr was asserted). All three carry the SAME
+        # documented N/(N-1) convention shift as T_hr: R's YatchewTest uses
+        # base-R `var()` (1/(N-1)) where our estimator uses the population
+        # (1/G, 1/2G) convention, so R = ours x G/(G-1). Measured to lock at
+        # ~1e-15 across every row (incl. the literal-2G sigma2_diff divisor).
+        r_yatchew_p = _as_list(r_result["yatchew_p"])
+        r_yatchew_s2_lin = _as_list(r_result["yatchew_sigma2_lin"])
+        r_yatchew_s2_diff = _as_list(r_result["yatchew_sigma2_diff"])
+        from scipy.stats import norm as _norm
+
         for r_idx, r_id in enumerate(r_event_ids):
             e = _r_id_to_event_time(int(r_id), trends_lin)
             if e not in dy_dict:
@@ -459,6 +470,35 @@ def test_yatchew_t_stat_parity(
                     f"after N/(N-1) convention shift"
                 ),
             )
+            # p-value: shift T first, THEN evaluate the normal CDF (a p-value
+            # cannot be rescaled directly). 1 - Phi(T_hr x G/(G-1)) vs golden.
+            _fac = G_horizon / (G_horizon - 1)
+            np.testing.assert_allclose(
+                1.0 - _norm.cdf(py_t_in_r_convention),
+                float(r_yatchew_p[r_idx]),
+                atol=YATCHEW_ATOL,
+                rtol=0,
+                err_msg=(
+                    f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: "
+                    f"p-value mismatch after N/(N-1) shift"
+                ),
+            )
+            # sigma^2 columns carry the same G/(G-1) convention shift as T_hr
+            # (sigma2_diff holds despite our literal-2G divisor — measured).
+            np.testing.assert_allclose(
+                float(r.sigma2_lin) * _fac,
+                float(r_yatchew_s2_lin[r_idx]),
+                atol=YATCHEW_ATOL,
+                rtol=0,
+                err_msg=(f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: sigma2_lin mismatch"),
+            )
+            np.testing.assert_allclose(
+                float(r.sigma2_diff) * _fac,
+                float(r_yatchew_s2_diff[r_idx]),
+                atol=YATCHEW_ATOL,
+                rtol=0,
+                err_msg=(f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: sigma2_diff mismatch"),
+            )
 
 
 class TestFixtureMetadata: