refactor(v4): PR 4c-i of the 4.0 cut - decouple the TWFE event-study mode from the MultiPeriodDiD containers and migrate the MPD-fixture suites

TODO.md

@@ -86,6 +86,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
+| TWFE event-study one-way HC2-BM absorbed-FE R parity against the `mpd_absorbed_fe_did` golden (`clubSandwich::vcovCR(cluster=1:n, type="CR2")`): the within spec with hc2_bm auto-clusters at unit, so the one-way lane needs an explicit singleton-cluster column (`cluster=<row id>` + hc2_bm = CR2 with G=n, PT2018 §3.3) once the G=n CR2 path is verified; the hc2 lane is already pinned by `TestTWFEEventStudyAbsorbedFERParity`, and the MultiPeriodDiD lanes die at the 4.0 cut (PR 4c-ii) | `tests/test_estimators_vcov_type.py` | #v4-4c-i | Mid | Medium |
 | Re-validate the plan-review dual engine on `gpt-6.1-sol` with `tools/plan-review-eval`: re-pin every arm's composite `cli_version` to the live Claude Code + codex CLI pair, add a gpt-6.1-sol arm, run the harness against the gpt-5.6-sol baseline (the 2026-10 swap carried the model over from the code-review campaign 2 verdict without a plan-review campaign). | `tools/plan-review-eval/config/configs.json` | reviewer-eval gpt61sol PR-2 | Heavy | Low |
 | TWFE-weights tutorial section: `attgt_weights` (all three `type=` values side by side) + `decompose_twfe_weights` + `plot_twfe_weights` on `mpdta`, either as a section of `02_staggered_did.ipynb` or a new notebook registered in `docs/tutorials/index.rst`; follow-on PR by decision, the feature ships without it | `docs/tutorials/02_staggered_did.ipynb` | #812 | Mid | Medium |
 | Geo-experiments tutorial: the SDiD fits on the `generate_factor_data(...)` panel (conversions ~1,500, sigma-hat ~310) hit the Frank-Wolfe iteration cap - the placebo fit warns once and the bootstrap cross-check warns on 100/100 draws - yet the narrative presents the bootstrap SE as a paper-faithful robustness check. Pre-existing on `main` under BOTH backends (surfaced when the notebook outputs were refreshed for the shape-only fit diagnostic); the committed April outputs predate the aggregated non-convergence warning. Either rescale / re-tune the DGP (or raise `max_iter` / loosen `min_decrease`) so the refits converge and re-execute, or label the bootstrap result unreliable and drop it from the "Validity evidence" claim. | `docs/tutorials/18_geo_experiments.ipynb` | local codex review of the SDiD shape-only fit-RMSE PR | Mid | Low |

changelog.d/20261024-v4-pr4c-i-twfe-es-decoupling.md

@@ -0,0 +1,26 @@
+### Internal
+- **The TWFE event-study surface is built directly from the shared core's payload**
+  (4.0 cut, PR 4c-i; ledger rows `M-010` / `M-011`): `DifferenceInDifferences._fit_event_study_core`
+  returns a private `_EventStudyCorePayload`, `TwoWayFixedEffects.fit(..., event_study=True)`
+  builds its `EventStudyResults` from it via `results_base._from_event_study_payload`, and only
+  `MultiPeriodDiD.fit` still wraps the payload into `MultiPeriodDiDResults`. `_from_mpd` and the
+  payload route share one calendar-surface builder, so they cannot drift; a committed oracle
+  captured on the pre-refactor tree (`tests/fixtures/v4_twfe_event_study_oracles.json`, capture
+  script `tests/_capture_v4_merge_mpd_oracles.py`) pins every field of the TWFE surface and of
+  its MultiPeriodDiD twin across eleven variance lanes. No user-visible behavior changes.
+- Spec item 4c is split into 4c-i (this PR) and 4c-ii (the removals); the 4c-ii bullet records
+  the `spec="within"` benchmark-arm relabel and the `M-016` view's record type (user decisions
+  2026-10-24). Ledger `M-010` / `M-011` notes updated; no status flips.
+
+### Testing
+- The MultiPeriodDiD-fixture suites (honest_did, pretrends, methodology pretrends/honest,
+  visualization, Conley, survey, replicate weights, vcov_type, df_convention, the estimators
+  event-study behaviour pins, the event-study surface roster, the DiagnosticReport /
+  BusinessReport generic fixtures) now run on `TwoWayFixedEffects().fit(..., event_study=True)`
+  (or a static TWFE fit for the report classes); MultiPeriodDiD's own tests, the MPD-only
+  designs (multi-`absorb`, `fixed_effects=`, wild bootstrap, container-field assertions, the
+  MPD-branch consumer tests) and the MPD-vs-TWFE parity oracles are unchanged until PR 4c-ii.
+  New `tests/test_twfe_event_study_decoupling.py`, shared helpers
+  `tests/helpers/twfe_es_oracle_lanes.py`, `event_study_average.py` and
+  `event_study_surface_builders.py`; the variance-conventions inventory gains the two TWFE
+  event-study rows; a TODO row tracks the one-way HC2-BM absorbed-FE R-parity lane.

diff_diff/estimators.py

@@ -28,6 +28,7 @@
     solve_ols,
 )
 from diff_diff.results import DiDResults, MultiPeriodDiDResults, PeriodEffect
+from diff_diff.results_base import _EventStudyCorePayload
 from diff_diff.utils import (
     WildBootstrapResults,
     absorbed_fe_cr1_k_increment,
@@ -1336,9 +1337,16 @@ def _fit_event_study_core(
         cluster_override: Any = _USE_SELF_CLUSTER,
         estimator_name: str = "MultiPeriodDiD",
         _frame_offset: int = 0,
-    ) -> Tuple[MultiPeriodDiDResults, np.ndarray, np.ndarray]:
+    ) -> _EventStudyCorePayload:
         """Shared event-study estimation core (row M-010).
 
+        Returns a private :class:`~diff_diff.results_base._EventStudyCorePayload`
+        (PR 4c-i decoupling): ``MultiPeriodDiD.fit`` wraps it into the legacy
+        ``MultiPeriodDiDResults`` container via ``_mpd_results_from_payload``;
+        the TWFE event-study mode builds the unified ``EventStudyResults``
+        surface from it directly, so no MultiPeriodDiD container is
+        constructed on that path.
+
         The relocated body of ``MultiPeriodDiD.fit`` (verbatim through the
         3.9 merge; docs/v4-design.md section 4.1): validation, the
         pooled event-study design build, OLS via ``solve_ols``, survey /
@@ -2251,7 +2259,15 @@ def _refit_mp_absorb(w_r):
                     _bm_dof_avg = float(_dof_all[-1])
 
         # Extract period-specific treatment effects for ALL non-reference periods
-        period_effects = {}
+        # (index-aligned lists over non_ref_periods; the MultiPeriodDiD
+        # wrapper zips them into PeriodEffect records, the TWFE surface
+        # builder reads them as arrays).
+        pe_att: List[float] = []
+        pe_se: List[float] = []
+        pe_t: List[float] = []
+        pe_p: List[float] = []
+        pe_lo: List[float] = []
+        pe_hi: List[float] = []
         post_effect_values = []
         post_effect_indices = []
         # Per-period df PROVENANCE for the unified event-study surface
@@ -2284,14 +2300,12 @@ def _refit_mp_absorb(w_r):
                 else float("nan")
             )
 
-            period_effects[period] = PeriodEffect(
-                period=period,
-                effect=effect,
-                se=se,
-                t_stat=t_stat,
-                p_value=p_value,
-                conf_int=conf_int,
-            )
+            pe_att.append(effect)
+            pe_se.append(se)
+            pe_t.append(t_stat)
+            pe_p.append(p_value)
+            pe_lo.append(conf_int[0])
+            pe_hi.append(conf_int[1])
 
             if period in post_periods:
                 post_effect_values.append(effect)
@@ -2338,32 +2352,41 @@ def _refit_mp_absorb(w_r):
         n_treated = n_treated_raw
         n_control = n_control_raw
 
-        # Create coefficient dictionary (var_names uniqueness already enforced
-        # before the fit above).
-        coef_dict = {name: coef for name, coef in zip(var_names, coefficients)}
-
-        # Store results
-        _core_results = MultiPeriodDiDResults(
-            period_effects=period_effects,
+        # Hand off to the builders (PR 4c-i): MultiPeriodDiD wraps this into
+        # its container; the TWFE event-study mode builds the surface from it.
+        return _EventStudyCorePayload(
+            periods=list(non_ref_periods),
+            att=np.asarray(pe_att, dtype=float),
+            se=np.asarray(pe_se, dtype=float),
+            t_stat=np.asarray(pe_t, dtype=float),
+            p_value=np.asarray(pe_p, dtype=float),
+            ci_lower=np.asarray(pe_lo, dtype=float),
+            ci_upper=np.asarray(pe_hi, dtype=float),
+            df_by_period=es_df_used,
+            reference_period=reference_period,
+            pre_periods=pre_periods,
+            post_periods=post_periods,
+            alpha=self.alpha,
+            coefficients=coefficients,
+            coef_names=list(var_names),
+            vcov=vcov,
+            interaction_indices=interaction_indices,
+            residuals=residuals,
+            fitted_values=fitted,
+            r_squared=r_squared,
             avg_att=avg_att,
             avg_se=avg_se,
             avg_t_stat=avg_t_stat,
             avg_p_value=avg_p_value,
             avg_conf_int=avg_conf_int,
+            inference_df=(
+                float(_avg_df)
+                if _avg_df is not None and np.isfinite(_avg_df) and _avg_df > 0
+                else None
+            ),
             n_obs=len(y),
             n_treated=n_treated,
             n_control=n_control,
-            pre_periods=pre_periods,
-            post_periods=post_periods,
-            alpha=self.alpha,
-            coefficients=coef_dict,
-            vcov=vcov,
-            residuals=residuals,
-            fitted_values=fitted,
-            r_squared=r_squared,
-            reference_period=reference_period,
-            interaction_indices=interaction_indices,
-            survey_metadata=survey_metadata,
             # Report the family that actually produced the SE; may be the
             # remapped hc1 under a replicate-weight survey design.
             vcov_type=_fit_vcov_type,
@@ -2373,15 +2396,59 @@ def _refit_mp_absorb(w_r):
             ),
             conley_lag_cutoff=(self.conley_lag_cutoff if _fit_vcov_type == "conley" else None),
             df_convention=self.df_convention,
-            inference_df=(
-                float(_avg_df)
-                if _avg_df is not None and np.isfinite(_avg_df) and _avg_df > 0
-                else None
-            ),
-            event_study_df=es_df_used,
+            survey_metadata=survey_metadata,
         )
 
-        return _core_results, coefficients, vcov
+
+def _mpd_results_from_payload(p: _EventStudyCorePayload) -> MultiPeriodDiDResults:
+    """Wrap the event-study core's payload into the legacy MultiPeriodDiD
+    container (deleted with the class at the 4.0 cut, PR 4c-ii).
+
+    ``period_effects`` is built in ``p.periods`` order (the core's sorted
+    non-reference order, the same insertion order the core used before the
+    PR 4c-i decoupling); per-period scalars round-trip through float64
+    arrays (value-identical; ``np.float64`` is a ``float`` subclass).
+    """
+    period_effects = {
+        per: PeriodEffect(
+            period=per,
+            effect=p.att[i],
+            se=p.se[i],
+            t_stat=p.t_stat[i],
+            p_value=p.p_value[i],
+            conf_int=(p.ci_lower[i], p.ci_upper[i]),
+        )
+        for i, per in enumerate(p.periods)
+    }
+    return MultiPeriodDiDResults(
+        period_effects=period_effects,
+        avg_att=p.avg_att,
+        avg_se=p.avg_se,
+        avg_t_stat=p.avg_t_stat,
+        avg_p_value=p.avg_p_value,
+        avg_conf_int=p.avg_conf_int,
+        n_obs=p.n_obs,
+        n_treated=p.n_treated,
+        n_control=p.n_control,
+        pre_periods=p.pre_periods,
+        post_periods=p.post_periods,
+        alpha=p.alpha,
+        coefficients=dict(zip(p.coef_names, p.coefficients)),
+        vcov=p.vcov,
+        residuals=p.residuals,
+        fitted_values=p.fitted_values,
+        r_squared=p.r_squared,
+        reference_period=p.reference_period,
+        interaction_indices=p.interaction_indices,
+        survey_metadata=p.survey_metadata,
+        vcov_type=p.vcov_type,
+        cluster_name=p.cluster_name,
+        n_clusters=p.n_clusters,
+        conley_lag_cutoff=p.conley_lag_cutoff,
+        df_convention=p.df_convention,
+        inference_df=p.inference_df,
+        event_study_df=p.df_by_period,
+    )
 
 
 class MultiPeriodDiD(DifferenceInDifferences):
@@ -2624,7 +2691,7 @@ def fit(  # type: ignore[override]
             )
             effective_inference = "analytical"
 
-        results, coefficients, vcov = self._fit_event_study_core(
+        payload = self._fit_event_study_core(
             data,
             outcome,
             treatment,
@@ -2639,9 +2706,9 @@ def fit(  # type: ignore[override]
             effective_inference,
             _frame_offset=1,
         )
-        self.results_ = results
-        self._coefficients = coefficients
-        self._vcov = vcov
+        self.results_ = _mpd_results_from_payload(payload)
+        self._coefficients = payload.coefficients
+        self._vcov = payload.vcov
         self.is_fitted_ = True
         return self.results_