feat: add LWDiD estimator (Lee & Wooldridge 2025, 2026)

CHANGELOG.md

@@ -5,6 +5,76 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+- **`LWDiD` (Lee & Wooldridge 2025, 2026 rolling-transformation DiD).** Unit-specific
+  demean/detrend (plus quarterly `demeanq`/`detrendq`) converts panel data to
+  cross-sectional transformed outcomes; supports common timing and staggered
+  adoption with never-treated / not-yet-treated controls,
+  `estimation_method` in `'reg'`/`'ipw'`/`'dr'`/`'psm'`, analytical
+  (`vcov_type` in `'classical'`/`'hc1'`/`'hc2'`/`'hc3'`) and cluster-robust
+  (constructor `cluster=`) inference, multiplier bootstrap, wild cluster
+  bootstrap, and randomization inference. Common-timing fits expose the same
+  post-fit event-study surface as staggered ones —
+  `results.aggregate('event_study')` returns per-period effects on the
+  calendar-time axis, so no separate per-period fit option exists.
+
+### Changed
+- **`LWDiD` API canonicalized to the v4 vocabulary agreed in PR #588's review**
+  (renames relative to the PR's earlier review rounds; nothing here was ever
+  released): `estimator=` -> `estimation_method=` with values `'ra'` -> `'reg'`
+  and `'ipwra'` -> `'dr'`; `vce=` -> `vcov_type=` with no `'cluster'` value —
+  cluster-robust (CR1) inference activates via the constructor `cluster=`
+  column instead; `bootstrap_seed=` -> `seed=` (default `None`);
+  `trim_threshold=` -> `pscore_trim=`. `vcov_type='hc3'` is computed through
+  the shared `diff_diff.linalg` HC machinery used by the other estimators,
+  and the `'hc0'`/`'hc4'` values are removed from the surface.
+  Unit-constancy validation is centralized and applies uniformly to
+  covariates and the cluster column across all estimation paths.
+
+### Removed
+- **`LWDiD` pre-v4 review-round surface** (never released): the `LW` alias,
+  the functional `lwdid()` wrapper, the `lwdid_trend_diagnostics` module
+  (including `recommend_transformation`), and the `overall_att` /
+  `period_effects` result fields together with the `period_specific` fit
+  option — per-period effects are served by the post-fit
+  `results.aggregate('event_study')` surface instead.
+
+### Fixed
+- **`LWDiD` review-round fixes** (staggered contract and inference tightenings):
+  - Staggered classical/HC SEs now come from the joint influence function
+    across cohort-time cells (the LW 2026 eq. 7.19 pooled-regression basis),
+    accounting for correlation among cohort effects that share controls
+    instead of assuming independence.
+  - On unbalanced panels the overall ATT point estimate is unified to the
+    eq. (7.18) composite-regression estimand `tau_omega`; the joint
+    influence function contributes the standard error only, so switching
+    variance options no longer moves the point estimate (gated to
+    `rolling` in `'demean'`/`'detrend'` with `control_group='never_treated'`,
+    `estimation_method='reg'`, and no covariates; the quarterly variants
+    keep their previous behavior).
+  - t-test degrees of freedom are computed from one design-based rule across
+    common-timing and staggered paths instead of two inconsistent ones.
+  - All-eventually-treated panels under `control_group='not_yet_treated'`
+    raise `ValueError` instead of silently truncating the sample; staggered
+    `covariates` must be unit-constant, time-varying columns raise
+    `ValueError`.
+  - Randomization inference uses the inclusive Phipson-Smyth rule
+    p = (c+1)/(B+1) and counts ties as extreme (`>=`), so p is never 0 and
+    an all-tie permutation distribution yields p = 1.0.
+  - `estimation_method='dr'` without covariates warns (`UserWarning`) that it
+    reduces to regression adjustment instead of silently doing so.
+  - `sensitivity_analysis` gains a `not_estimable` robustness level (with a
+    warning) when the ratio cannot be computed — including the zero-baseline
+    case — instead of mislabeling it.
+  - `to_dict()` output is fully JSON-native, including datetime/Period
+    cohort and time labels (ISO-8601 / period strings, NaT -> None).
+  - Staggered fits accept datetime64 and Period time scales, and panels
+    mixing the two time families are rejected in both directions with a
+    clear `ValueError`; cluster variable equal to the unit column no longer
+    raises a spurious column-lookup error.
+
 ## [3.9.1] - 2026-08-17
 
 ### Added

README.md

@@ -120,6 +120,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`.
 - [WooldridgeDiD](https://diff-diff.readthedocs.io/en/stable/api/wooldridge_etwfe.html) - Wooldridge (2023, 2025) ETWFE: saturated OLS, logit/Poisson QMLE (ASF-based ATT). Alias `ETWFE`.
 - [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html) - Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting), variance- or equally-weighted ATT, for absorbing or non-absorbing (reversible) treatment
 - [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html) - Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: full counterfactual distribution and quantile treatment effects via CDF transformation, plus the QDiD comparison estimator via `method="qdid"`; bootstrap inference; R qte parity. Alias `CiC`
+- [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html) - Lee & Wooldridge (2025, 2026) rolling-transformation DiD: unit-specific demean/detrend converts panel to cross-section, staggered adoption, `estimation_method` in `reg`/`ipw`/`dr`/`psm` (the papers' RA/IPW/IPWRA plus propensity-score matching), exact small-N inference
 - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings
 
 ## Diagnostics & Sensitivity

diff_diff/__init__.py

@@ -160,6 +160,8 @@
 )
 from diff_diff.lpdid import LPDiD
 from diff_diff.lpdid_results import LPDiDResults
+from diff_diff.lwdid import LWDiD
+from diff_diff.lwdid_results import LWDiDResults
 from diff_diff.mmm import (
     MeridianROIPrior,
     meridian_calibration_mask,
@@ -459,6 +461,9 @@ def __getattr__(name: str) -> _Any:
     # LPDiD (Local Projections DiD)
     "LPDiD",
     "LPDiDResults",
+    # LWDiD (Lee & Wooldridge rolling transformation DiD)
+    "LWDiD",
+    "LWDiDResults",
     # Visualization
     "plot_bacon",
     "plot_event_study",

diff_diff/guides/llms.txt

@@ -80,6 +80,7 @@ The site is organized into 5 sections, each with a landing page:
 - [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html): Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting); variance- or equally-weighted ATT, premean differencing, pooled pre/post, fast. Absorbing by default; non-absorbing (reversible) treatment via `non_absorbing="first_entry"` (Eq. 12) or `"effect_stabilization"` (Eq. 13, window `L`). Complex-survey designs (pweight + stratified-PSU TSL SEs) on the default path via `fit(survey_design=...)`.
 - [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html): Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: recovers the treated group's full counterfactual outcome distribution and quantile treatment effects (ATT + QTE grid) via the CDF transformation `F_10(F_00^{-1}(F_01(y)))`; invariant to monotone outcome transformations (unconditional fits; the covariate QR branch is not); bootstrap inference (panel or repeated cross-section resampling); point parity with R `qte::CiC()`, including its covariate branch (`covariates=` -> per-cell linear quantile regression, Melly-Santangelo-style conditional CiC). Continuous outcomes, numeric covariates. Alias `CiC`.
 - [QDiD](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html): **Deprecated 3.9, removed 4.0 - use `ChangesInChanges(method="qdid")`.** Athey & Imbens (2006) quantile DiD comparison estimator (additive quantile-by-quantile DiD, matching R `qte::QDiD()` including its covariate branch via `covariates=`); same bootstrap machinery as ChangesInChanges. The paper recommends CiC over QDiD (scale-dependent model with testable restrictions; a non-monotonicity warning fires when violated - unconditional fits only, the covariate-path counterfactual quantile curve is monotone by construction).
+- [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html): Lee & Wooldridge (2025, 2026) rolling-transformation DiD — unit-specific demean/detrend converts panel to cross-section, supports staggered adoption with flexible control groups. Signature: `LWDiD(rolling='demean', estimation_method='reg', vcov_type='hc1', cluster=None, control_group='not_yet_treated', alpha=0.05, n_bootstrap=0, seed=None, pscore_trim=0.01, n_neighbors=1, caliper=None, with_replacement=True, n_jobs=1).fit(data, outcome, unit, time, treatment, first_treat=None, covariates=None)`. `estimation_method` values: `reg` (papers' RA), `ipw`, `dr` (papers' IPWRA, doubly robust), `psm`; `vcov_type` values: `classical`/`hc1`/`hc2`/`hc3`; cluster-robust inference via the constructor's `cluster=` column (not a `vcov_type` value). Per-period effects: post-fit `results.aggregate('event_study')`.
 - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html): Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings
 
 ## Diagnostics and Sensitivity Analysis

diff_diff/linalg.py

@@ -1913,7 +1913,7 @@ def _solve_ols_numpy(
         return coefficients, residuals, vcov
 
 
-_VALID_VCOV_TYPES = frozenset({"classical", "hc1", "hc2", "hc2_bm", "conley"})
+_VALID_VCOV_TYPES = frozenset({"classical", "hc1", "hc2", "hc2_bm", "hc3", "conley"})
 
 
 def _validate_vcov_args(
@@ -1936,7 +1936,7 @@ def _validate_vcov_args(
     ValueError
         If ``vcov_type`` is not in the allowed set, or if ``cluster_ids`` is
         combined with a ``vcov_type`` that is one-way only (``classical``,
-        ``hc2``).
+        ``hc2``, ``hc3``).
     NotImplementedError
         If ``vcov_type == "conley"`` is combined with ``weights`` (regardless
         of ``weight_type``: weighted Conley is not implemented on the
@@ -1956,7 +1956,7 @@ def _validate_vcov_args(
     # Mirrored K_reference-adjustment contract for direct compute_robust_vcov
     # / kernel callers (solve_ols routes enforce it at its own front door).
     _validate_cluster_k_adjustment(cluster_k_adjustment, cluster_ids, vcov_type)
-    if vcov_type in ("classical", "hc2") and cluster_ids is not None:
+    if vcov_type in ("classical", "hc2", "hc3") and cluster_ids is not None:
         msg = {
             "classical": (
                 "classical SEs are one-way only; pass vcov_type='hc1' or "
@@ -1965,6 +1965,10 @@ def _validate_vcov_args(
             "hc2": (
                 "hc2 is one-way only. Use vcov_type='hc2_bm' for " "cluster-robust Bell-McCaffrey."
             ),
+            "hc3": (
+                "hc3 is one-way only. Use vcov_type='hc1' (CR1) or "
+                "'hc2_bm' (CR2 Bell-McCaffrey) for cluster-robust."
+            ),
         }[vcov_type]
         raise ValueError(msg)
     # Weighted Bell-McCaffrey (both one-way and cluster) is now supported via
@@ -2068,7 +2072,7 @@ def resolve_vcov_type(
       ``"hc1"`` and ``robust=False`` to ``"classical"``.
     - If ``vcov_type`` is supplied: it must be one of the values in the
       module-level ``_VALID_VCOV_TYPES`` set, namely
-      ``{"classical", "hc1", "hc2", "hc2_bm", "conley"}``.
+      ``{"classical", "hc1", "hc2", "hc2_bm", "hc3", "conley"}``.
     - If ``robust=False`` is supplied together with a non-``"classical"`` ``vcov_type``,
       raise ``ValueError`` - the combination is ambiguous.
 
@@ -2086,7 +2090,8 @@ def resolve_vcov_type(
     Returns
     -------
     str
-        One of ``"classical"``, ``"hc1"``, ``"hc2"``, ``"hc2_bm"``, ``"conley"``.
+        One of ``"classical"``, ``"hc1"``, ``"hc2"``, ``"hc2_bm"``,
+        ``"hc3"``, ``"conley"``.
 
     Raises
     ------
@@ -2134,7 +2139,7 @@ def compute_robust_vcov(
     conley_lag_cutoff: Optional[int] = None,
 ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
     """
-    Compute variance-covariance matrix under one of five `vcov_type` variants.
+    Compute variance-covariance matrix under one of six `vcov_type` variants.
 
     Uses the sandwich estimator: (X'X)^{-1} * meat * (X'X)^{-1}, with the meat
     matrix determined by the ``vcov_type`` dispatch:
@@ -2149,6 +2154,9 @@ def compute_robust_vcov(
       ``sum_i (u_i^2 / (1 - h_ii)) x_i x_i'`` where ``h_ii`` are hat-matrix
       diagonals. No DOF adjustment beyond ``n - k``. One-way only; errors with
       ``cluster_ids``.
+    - ``"hc3"``: jackknife-style leverage correction, meat
+      ``sum_i (u_i^2 / (1 - h_ii)^2) x_i x_i'`` (matches ``sandwich::vcovHC``
+      type="HC3": no DOF factor). One-way only; errors with ``cluster_ids``.
     - ``"hc2_bm"``: one-way HC2 meat plus Imbens-Kolesar (2016) Bell-McCaffrey
       Satterthwaite degrees of freedom per coefficient when ``cluster_ids`` is
       ``None``. When ``cluster_ids`` is supplied, dispatches to the
@@ -2195,8 +2203,8 @@ def compute_robust_vcov(
         Weight type: "pweight", "fweight", or "aweight".
     vcov_type : str, default "hc1"
         One of ``"classical"``, ``"hc1"``, ``"hc2"``, ``"hc2_bm"``,
-        ``"conley"`` (see top-level docstring above for the dispatch
-        contract).
+        ``"hc3"``, ``"conley"`` (see top-level docstring above for the
+        dispatch contract).
     conley_coords : ndarray of shape (n, 2), optional, keyword-only
         Required when ``vcov_type="conley"``. Two-column array of
         ``[lat, lon]`` (degrees, for ``conley_metric="haversine"``) or
@@ -2221,8 +2229,9 @@ def compute_robust_vcov(
     return_dof : bool, default False
         When True, returns ``(vcov, dof_vec)`` tuple. ``dof_vec`` is a length-k
         array of per-coefficient degrees of freedom. For ``classical``,
-        ``hc1``, ``hc2``: every element is ``n_eff - k``. For ``hc2_bm``
-        one-way: Imbens-Kolesar (2016) Satterthwaite DOF per contrast.
+        ``hc1``, ``hc2``, ``hc3``: every element is ``n_eff - k``. For
+        ``hc2_bm`` one-way: Imbens-Kolesar (2016) Satterthwaite DOF per
+        contrast.
     cluster_k_adjustment : int, default 0, keyword-only
         Signed K_reference adjustment added to the visible column count in
         the CLUSTERED CR1 finite-sample factor only (absorbed FE not nested
@@ -3550,9 +3559,9 @@ def _compute_robust_vcov_numpy(
         return vcov_cr2
 
     # ------------------------------------------------------------------
-    # HC2 / HC2+BM one-way (no cluster).
+    # HC2 / HC2+BM / HC3 one-way (no cluster).
     # ------------------------------------------------------------------
-    if vcov_type in ("hc2", "hc2_bm"):
+    if vcov_type in ("hc2", "hc2_bm", "hc3"):
         # cluster path handled above; here cluster_ids is None by construction.
         # **Weighted hc2_bm one-way**: clubSandwich's CR2 with singleton clusters
         # uses the bias-corrected adjustment `A_i = 1 / sqrt(G_i)` where
@@ -3598,26 +3607,29 @@ def _compute_robust_vcov_numpy(
                 return_dof=return_dof,
             )
         one_minus_h = np.maximum(1.0 - h_diag, 1e-10)
-        # HC2 meat: sum_i (u_i^2 / (1 - h_ii)) x_i x_i', with pweight scaling
-        # matching the HC1 convention (w_i * u_i / sqrt(1 - h_ii) as score).
+        # HC2 meat: sum_i (u_i^2 / (1 - h_ii)) x_i x_i'; HC3 squares the
+        # leverage denominator (jackknife-style, sandwich::vcovHC type="HC3").
+        # pweight scaling matches the HC1 convention (w_i * u_i / sqrt(denom)
+        # as score).
+        lev_denom = one_minus_h**2 if vcov_type == "hc3" else one_minus_h
         if weights is not None and weight_type == "fweight":
-            factor = weights * (residuals**2) / one_minus_h
+            factor = weights * (residuals**2) / lev_denom
             meat = X.T @ (X * factor[:, np.newaxis])
         elif weights is not None and weight_type == "pweight":
-            # pweight scores carry w in the score, so meat = sum (w u / sqrt(1-h))^2 x x'
-            scaled = weights * residuals / np.sqrt(one_minus_h)
+            # pweight scores carry w in the score, so meat = sum (w u / sqrt(denom))^2 x x'
+            scaled = weights * residuals / np.sqrt(lev_denom)
             scores_hc2 = X * scaled[:, np.newaxis]
             meat = scores_hc2.T @ scores_hc2
         else:
-            # aweight / unweighted: meat = sum_i (u_i^2 / (1 - h_ii)) x_i x_i'
-            factor = (residuals**2) / one_minus_h
+            # aweight / unweighted: meat = sum_i (u_i^2 / denom_i) x_i x_i'
+            factor = (residuals**2) / lev_denom
             # Zero out zero-weight rows under aweight (subpopulation invariance)
             if weights is not None and np.any(weights == 0):
                 factor = factor * (weights > 0)
             meat = X.T @ (X * factor[:, np.newaxis])
 
-        # Sandwich without DOF adjustment for HC2 (matches sandwich::vcovHC
-        # type="HC2" convention: no (n/(n-k)) factor).
+        # Sandwich without DOF adjustment for HC2/HC3 (matches sandwich::vcovHC
+        # type="HC2"/"HC3" convention: no (n/(n-k)) factor).
         try:
             temp = np.linalg.solve(bread_matrix, meat)
             vcov = np.linalg.solve(bread_matrix, temp.T).T
@@ -3632,7 +3644,7 @@ def _compute_robust_vcov_numpy(
 
         if not return_dof:
             return vcov
-        if vcov_type == "hc2":
+        if vcov_type in ("hc2", "hc3"):
             dof_vec = np.full(k, n_eff - k, dtype=np.float64)
         else:  # hc2_bm
             dof_vec = _compute_bm_dof_oneway(X, bread_matrix, h_diag, weights=weights)