Round 18: fix sigma_fe to use paper w_{g,t} (not contribution weights)

P1: Round 17's sigma_fe fix still used the contribution weights
w_gt = s * w_paper as if they were the paper's w_{g,t} object. The
paper's Corollary 1 defines w_{g,t} = eps / E_treated[eps], which
does NOT include the observation-share factor in the numerator. The
contribution weights (which include the share factor) are the right
object for the TWFE decomposition itself, but sigma_fe needs the
paper's w_{g,t} centered at 1.

The fix computes w_paper separately:
  w_paper = eps_treated / sum(s * eps_treated)
  sigma(w) = sqrt(sum(s * (w_paper - 1)^2))
  sigma_fe = |beta_fe| / sigma(w)

The exported weights column still contains the contribution weights
(the TWFE decomposition object), not the paper weights.

P2: added test_twfe_diagnostic_hand_checkable_sigma_fe in
TestTWFEDiagnostic with a 4-group 3-period staggered panel.
Asserts beta_fe = 3.5, fraction_negative = 0.0, and
sigma_fe = 6.8641 numerically. This test would have caught both
the Round 17 and pre-Round-17 sigma_fe formula errors.

Test counts: 119 -> 120. Black, ruff clean.

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

Author: igerber

diff_diff/chaisemartin_dhaultfoeuille.py

@@ -2091,18 +2091,30 @@ def _compute_twfe_diagnostic(
 
     # Step 6: sigma_fe per Corollary 1 of AER 2020
     #
-    # sigma_fe = |beta_fe| / sigma(w), where sigma(w) is the
-    # observation-share-weighted standard deviation of w_{g,t} across
-    # treated cells, using shares s_{g,t} = N_{g,t} / N_1 where
-    # N_1 = sum of N_{g,t} over treated cells. This matches the
-    # paper's population-expectation formula for sigma(w).
-    w_treated = w_gt[treated_mask]
+    # The paper defines w_{g,t} = eps_{g,t} / E_treated[eps], which
+    # is DIFFERENT from the contribution weights w_gt exported in the
+    # weights DataFrame (contribution_weight = s * w_paper). The paper
+    # weight has the property that sum(s * w_paper) = 1 (centered at
+    # 1 under observation-share weighting). sigma_fe uses the paper
+    # weight:
+    #
+    #   w_paper = eps / sum_treated(s * eps)
+    #   sigma(w) = sqrt(sum_treated(s * (w_paper - 1)^2))
+    #   sigma_fe = |beta_fe| / sigma(w)
+    #
+    # where s_{g,t} = N_{g,t} / N_1 are observation shares.
+    eps_treated = eps[treated_mask]
     n_treated_arr = n_arr[treated_mask]
     n1 = float(n_treated_arr.sum())  # total treated observations
     if n1 > 0:
-        shares = n_treated_arr / n1  # observation shares
-        w_bar = float((shares * w_treated).sum())
-        var_w = float((shares * (w_treated - w_bar) ** 2).sum())
+        shares = n_treated_arr / n1  # s_{g,t} = N_{g,t} / N_1
+        denom_paper = float((shares * eps_treated).sum())
+        if abs(denom_paper) > 0:
+            w_paper = eps_treated / denom_paper  # paper's w_{g,t}
+            # Weighted variance around 1 (the weighted mean of w_paper is 1 by construction)
+            var_w = float((shares * (w_paper - 1.0) ** 2).sum())
+        else:
+            var_w = 0.0
     else:
         var_w = 0.0
     if var_w > 0 and np.isfinite(beta_fe):

tests/test_methodology_chaisemartin_dhaultfoeuille.py

@@ -445,6 +445,58 @@ def test_twfe_diagnostic_runs_on_real_data(self):
         # For simplicity, just verify the weights array is not all zero
         assert (weights_df["weight"] != 0).any()
 
+    def test_twfe_diagnostic_hand_checkable_sigma_fe(self):
+        """
+        Hand-checkable TWFE diagnostic on a 4-group 3-period panel with
+        staggered treatment (g1 at t=1, g2 at t=2, g3-g4 never).
+
+        Expected values computed analytically (equal cell sizes):
+        - beta_fe = 3.5 (TWFE coefficient from OLS of y on FE + d)
+        - Treated cells: (g1,t1), (g1,t2), (g2,t2) with contribution
+          weights [0.4, 0.1, 0.5]
+        - Paper weights w_{g,t} (Corollary 1): [1.2, 0.3, 1.5]
+          (contribution_weight / share, centered at 1.0)
+        - sigma(w) = sqrt(sum(s * (w_paper - 1)^2)) = 0.5099
+        - sigma_fe = |3.5| / 0.5099 = 6.8641
+        - fraction_negative = 0.0 (all treated weights positive)
+        """
+        df = pd.DataFrame(
+            {
+                "group": [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4],
+                "period": [0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2],
+                "treatment": [0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+                "outcome": [
+                    10,
+                    14,
+                    15,
+                    10,
+                    11,
+                    16,
+                    10,
+                    11,
+                    12,
+                    10,
+                    11,
+                    12,
+                ],
+            }
+        )
+        from diff_diff import twowayfeweights
+
+        result = twowayfeweights(
+            df,
+            outcome="outcome",
+            group="group",
+            time="period",
+            treatment="treatment",
+        )
+        # beta_fe: the plain TWFE coefficient
+        assert result.beta_fe == pytest.approx(3.5, abs=0.01)
+        # fraction_negative: all treated weights positive
+        assert result.fraction_negative == pytest.approx(0.0)
+        # sigma_fe: the Corollary 1 sign-flip threshold
+        assert result.sigma_fe == pytest.approx(6.8641, abs=0.01)
+
     def test_twfe_disabled_means_none(self):
         data = generate_reversible_did_data(n_groups=30, n_periods=4, seed=1)
         est = ChaisemartinDHaultfoeuille(twfe_diagnostic=False)