Fix CI review Round 1: NaN-consistency, overall surface, rank_deficient_action

P0: Cumulated DID^{fd} SE now requires ALL component SEs to be finite;
    non-finite SE at any horizon propagates NaN (was silently dropped).

P1: trends_linear + L_max>=2 overall surface now reports cumulated level
    effects from linear_trends_effects (was second-difference delta).
    cost_benefit_delta suppressed under trends_linear (meaningless on
    second-differences).

P2: rank_deficient_action threaded through _compute_covariate_residualization
    and _compute_heterogeneity_test (was hardcoded "warn").

P3: fit() docstrings updated for controls, trends_linear, trends_nonparam
    (were stale "Reserved for Phase 3" text).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: igerber

diff_diff/chaisemartin_dhaultfoeuille.py

@@ -525,14 +525,18 @@ def fit(
             Must be a positive integer not exceeding the number of
             post-baseline periods in the panel.
         controls : list of str, optional
-            **Reserved for Phase 3** (covariate adjustment via the
-            residualization-style ``DID^X`` from Web Appendix Section 1.2
-            of the dynamic paper).
+            Column names for covariate adjustment via residualization-style
+            ``DID^X`` (Web Appendix Section 1.2). Requires ``L_max >= 1``.
+            One ``theta_hat`` per baseline treatment value, estimated by
+            OLS on not-yet-treated observations. NOT doubly-robust.
         trends_linear : bool, optional
-            **Reserved for Phase 3** (group-specific linear trends via
-            ``DID^{fd}``).
-        trends_nonparam : Any, optional
-            **Reserved for Phase 3** (state-set-specific trends).
+            If ``True``, estimate group-specific linear trends via
+            ``DID^{fd}`` (Web Appendix Section 1.3, Lemma 6). Requires
+            ``L_max >= 1`` and at least 3 time periods.
+        trends_nonparam : str, optional
+            Column name for state-set membership. Restricts the control
+            pool to groups in the same set (Web Appendix Section 1.4).
+            Requires ``L_max >= 1`` and time-invariant values per group.
         honest_did : bool, default=False
             **Reserved for Phase 3** (HonestDiD integration on placebos).
         survey_design : Any, optional
@@ -959,6 +963,7 @@ def fit(
                 N_mat=N_mat,
                 baselines=baselines,
                 first_switch_idx=first_switch_idx_arr,
+                rank_deficient_action=self.rank_deficient_action,
             )
             # Keep raw Y_mat for the per-period DID path (which does not
             # support covariate residualization - it uses binary joiner/leaver
@@ -2073,12 +2078,20 @@ def fit(
                 cum_effect = float(
                     np.sum(S_arr[eligible] * running_per_group[eligible]) / N_l
                 )
-                # SE: conservative upper bound (sum of per-horizon SEs)
-                running_se_ub = sum(
-                    event_study_effects.get(ll, {}).get("se", 0.0)
-                    for ll in range(1, l_h + 1)
-                    if np.isfinite(event_study_effects.get(ll, {}).get("se", np.nan))
-                ) if event_study_effects is not None else float("nan")
+                # SE: conservative upper bound (sum of per-horizon SEs).
+                # NaN-consistency: if ANY component SE up to horizon l is
+                # non-finite, the cumulated SE is NaN (not 0.0).
+                if event_study_effects is not None:
+                    component_ses = [
+                        event_study_effects.get(ll, {}).get("se", np.nan)
+                        for ll in range(1, l_h + 1)
+                    ]
+                    if all(np.isfinite(s) for s in component_ses):
+                        running_se_ub = sum(component_ses)
+                    else:
+                        running_se_ub = float("nan")
+                else:
+                    running_se_ub = float("nan")
                 cum_t, cum_p, cum_ci = safe_inference(
                     cum_effect, running_se_ub, alpha=self.alpha, df=None
                 )
@@ -2091,6 +2104,22 @@ def fit(
                 }
             linear_trends_effects = cumulated if cumulated else None
 
+        # When trends_linear=True and L_max>=2, suppress cost_benefit_delta
+        # (which is computed on second-differences) and set overall_* from
+        # the cumulated level effects instead. This prevents the results
+        # surface from labeling a second-difference aggregate as delta^{fd}
+        # (a level-effect estimand).
+        if _is_trends_linear and L_max is not None and L_max >= 2:
+            cost_benefit_result = None
+            if linear_trends_effects:
+                max_h = max(linear_trends_effects.keys())
+                lt = linear_trends_effects[max_h]
+                effective_overall_att = lt["effect"]
+                effective_overall_se = lt["se"]
+                effective_overall_t = lt["t_stat"]
+                effective_overall_p = lt["p_value"]
+                effective_overall_ci = lt["conf_int"]
+
         # ------------------------------------------------------------------
         # Heterogeneity testing (Web Appendix Section 1.5, Lemma 7)
         # ------------------------------------------------------------------
@@ -2130,6 +2159,7 @@ def fit(
                 X_het=X_het,
                 L_max=L_max,
                 alpha=self.alpha,
+                rank_deficient_action=self.rank_deficient_action,
             )
 
         twfe_weights_df = None
@@ -2634,6 +2664,7 @@ def _compute_covariate_residualization(
     N_mat: np.ndarray,
     baselines: np.ndarray,
     first_switch_idx: np.ndarray,
+    rank_deficient_action: str = "warn",
 ) -> Tuple[np.ndarray, Dict[str, Any]]:
     """Residualize outcomes by partialling out covariates per baseline treatment.
 
@@ -2750,7 +2781,7 @@ def _compute_covariate_residualization(
             design,
             dY,
             return_vcov=True,
-            rank_deficient_action="warn",
+            rank_deficient_action=rank_deficient_action,
         )
 
         # Extract covariate coefficients (first n_covariates entries)
@@ -2837,6 +2868,7 @@ def _compute_heterogeneity_test(
     X_het: np.ndarray,
     L_max: int,
     alpha: float = 0.05,
+    rank_deficient_action: str = "warn",
 ) -> Dict[int, Dict[str, Any]]:
     """Test for heterogeneous treatment effects (Web Appendix Section 1.5).
 
@@ -2938,7 +2970,7 @@ def _compute_heterogeneity_test(
         coefs, _residuals, vcov = solve_ols(
             design, dep_arr,
             return_vcov=True,
-            rank_deficient_action="warn",
+            rank_deficient_action=rank_deficient_action,
         )
 
         beta_het = float(coefs[0])

tests/test_chaisemartin_dhaultfoeuille.py

@@ -2550,6 +2550,58 @@ def test_trends_with_covariates(self):
         assert r.covariate_residuals is not None
         assert r.linear_trends_effects is not None
 
+    def test_trends_linear_lmax2_overall_surface(self):
+        """Overall surface under trends_linear + L_max>=2 uses cumulated level effects."""
+        df = self._make_panel_with_trends()
+        r = ChaisemartinDHaultfoeuille(seed=1).fit(
+            df, "outcome", "group", "period", "treatment",
+            L_max=3, trends_linear=True,
+        )
+        # overall_att should equal the cumulated level effect at max horizon
+        assert r.linear_trends_effects is not None
+        max_h = max(r.linear_trends_effects.keys())
+        cum_effect = r.linear_trends_effects[max_h]["effect"]
+        assert r.overall_att == pytest.approx(cum_effect, abs=1e-10)
+        # cost_benefit_delta should be suppressed (not computed on second-diffs)
+        assert r.cost_benefit_delta is None
+
+    def test_cumulated_se_nan_propagation(self):
+        """Cumulated SE is NaN when a component horizon has NaN SE."""
+        # Create a panel where horizon 2 has no eligible switchers (NaN SE)
+        # but horizon 1 does. The cumulated effect at h=2 should have NaN SE.
+        rng = np.random.RandomState(77)
+        rows = []
+        for g in range(30):
+            group_fe = rng.normal(0, 1)
+            # Groups 0-9: switch at period 3 (enough pre-switch for trends)
+            # Groups 10-19: never switch (controls)
+            # Groups 20-29: switch at period 4 (only 1 post-switch period)
+            if g < 10:
+                switch_t = 3
+            elif g < 20:
+                switch_t = 99
+            else:
+                switch_t = 4
+            for t in range(5):
+                d = 1 if t >= switch_t else 0
+                y = group_fe + t + 3 * d + rng.normal(0, 0.3)
+                rows.append({"group": g, "period": t, "treatment": d, "outcome": y})
+        df = pd.DataFrame(rows)
+        r = ChaisemartinDHaultfoeuille(seed=1).fit(
+            df, "outcome", "group", "period", "treatment",
+            L_max=2, trends_linear=True,
+        )
+        # If SE at horizon 1 is finite but horizon 2 is NaN,
+        # cumulated h=2 SE must be NaN (not 0.0)
+        if r.linear_trends_effects is not None and 2 in r.linear_trends_effects:
+            cum_se = r.linear_trends_effects[2]["se"]
+            es = r.event_study_effects
+            if es and 2 in es and not np.isfinite(es[2]["se"]):
+                assert not np.isfinite(cum_se), (
+                    f"Cumulated SE should be NaN when component h=2 SE is NaN, "
+                    f"got {cum_se}"
+                )
+
 
 class TestStateSetTrends:
     """State-set-specific trends (ROADMAP item 3c)."""