From 1a355265ea871158eb04ddb6f3423000681124b7 Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Sat, 4 Jul 2026 13:54:15 -0400
Subject: [PATCH] feat(continuous-did): covariate support (reg + dr) under
 conditional parallel trends

Add `covariates=` and `estimation_method in {"reg","dr"}` to ContinuousDiD for
dose-response estimation under conditional parallel trends
(E[dY(0)|D=d,X] = E[dY(0)|D=0,X]). Each (g,t) cell's control counterfactual
becomes a covariate-adjusted prediction:
- reg: outcome regression on controls (dY_i - X_i'gamma_hat)
- dr (default): doubly-robust (DRDID drdid_panel), plus a scalar augmentation

Scalar overall_att + SE match DRDID reg_did_panel / drdid_panel to ~1e-8;
analytical, multiplier-bootstrap, and event-study inference all compose with
covariates. reg and dr share the ATT(d) shape and ACRT(d), differing only in
the overall_att/ATT(d) level and the doubly-robust SE.

estimation_method="ipw" with covariates raises NotImplementedError (pure IPW's
covariate adjustment is a scalar level shift and cannot adjust the curve shape);
covariates + survey_design deferred. The no-covariate path is unchanged.

Validated: DRDID parity (reg/dr, skip-guarded), R-free NumPy influence-function
cross-check at p>=2 (in CI), MC coverage (reg 96%, dr 95%), DGP recovery, and
reg-vs-dr ACRT identity.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LHDijzf8zHXk5T8ahS2mKi
---
 CHANGELOG.md                             |  10 +
 TODO.md                                  |   2 +-
 diff_diff/continuous_did.py              | 435 +++++++++++++++++-
 diff_diff/continuous_did_results.py      |  13 +-
 diff_diff/guides/llms-full.txt           |  12 +
 docs/methodology/REGISTRY.md             |  20 +-
 docs/methodology/continuous-did.md       |   9 +-
 tests/test_continuous_did.py             | 163 ++++++-
 tests/test_methodology_continuous_did.py | 560 ++++++++++++++++++++---
 9 files changed, 1129 insertions(+), 95 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a452c55ce..b047c41ac 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **`ContinuousDiD` covariate support** (`covariates=`, `estimation_method ∈ {"reg", "dr"}`) for
+  dose-response estimation under **conditional** parallel trends
+  (`E[ΔY(0) | D=d, X] = E[ΔY(0) | D=0, X]`). `reg` uses an outcome-regression control counterfactual;
+  `dr` (default) is doubly-robust (DRDID `drdid_panel`). The scalar `overall_att` + standard error
+  match `DRDID::reg_did_panel` / `drdid_panel` to ~1e-8; analytical, multiplier-bootstrap, and
+  event-study inference all compose with covariates. `reg` and `dr` share the dose-response *shape*
+  and `ACRT(d)`, differing only in the `overall_att` / ATT(d) level and the doubly-robust SE.
+  `estimation_method="ipw"` with covariates raises `NotImplementedError` (pure IPW's covariate
+  adjustment is a scalar level shift and cannot adjust the curve shape); `covariates=` +
+  `survey_design=` is deferred (`NotImplementedError`).
 - **`HeterogeneousAdoptionDiD` event-study `cluster=` support** (both designs). On
   `aggregate="event_study"`, `cluster=` now produces cluster-robust per-horizon pointwise
   confidence intervals AND a cluster-robust simultaneous sup-t confidence band (`cband=True`),
diff --git a/TODO.md b/TODO.md
index 69c80cfa7..e29cf4419 100644
--- a/TODO.md
+++ b/TODO.md
@@ -29,7 +29,7 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
 | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
-| `ContinuousDiD` CGBS-2024 extensions (matches R `contdid` v0.1.0 deferral set): (a) `covariates=` kwarg; (b) discrete-treatment saturated regression (integer dose currently warned, not routed to per-level coefficients); (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. REGISTRY `## ContinuousDiD` → Implementation Checklist marks these `[ ]`. | `continuous_did.py` | CGBS-2024 | Heavy | Low |
+| `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; remaining: (b) discrete-treatment saturated regression (integer dose currently warned, not routed to per-level coefficients); (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. Also deferred from the covariate work: `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate — documented `NotImplementedError`), and `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF). | `continuous_did.py` | CGBS-2024 | Heavy | Low |
 | `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low |
 | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low |
 | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low |
diff --git a/diff_diff/continuous_did.py b/diff_diff/continuous_did.py
index cb4f8e4e1..8c86968df 100644
--- a/diff_diff/continuous_did.py
+++ b/diff_diff/continuous_did.py
@@ -29,7 +29,7 @@
     ContinuousDiDResults,
     DoseResponseCurve,
 )
-from diff_diff.linalg import _rank_guarded_inv, solve_ols
+from diff_diff.linalg import _rank_guarded_inv, solve_logit, solve_ols
 from diff_diff.survey import (
     _resolve_survey_for_fit,
     _validate_unit_constant_survey,
@@ -72,6 +72,36 @@ class ContinuousDiD:
         Random seed for reproducibility.
     rank_deficient_action : str, default="warn"
         Action for rank-deficient B-spline OLS: ``"warn"``, ``"error"``, or ``"silent"``.
+    covariates : list of str, optional
+        Column names of covariates for **conditional** parallel trends
+        (``E[ΔY(0) | D=d, X] = E[ΔY(0) | D=0, X]``). When ``None`` (default) the
+        estimator uses unconditional parallel trends. Covariates enter through a
+        covariate-adjusted per-cell control counterfactual (see ``estimation_method``).
+        Covariates are read from the base period of each ``(g, t)`` comparison.
+        Not currently composable with ``survey_design=`` (raises ``NotImplementedError``).
+    estimation_method : str, default="dr"
+        Covariate-adjustment method (only used when ``covariates`` is set):
+        ``"reg"`` (outcome regression) or ``"dr"`` (doubly-robust, the default).
+        ``"ipw"`` is **not supported on the dose / event-study aggregation** — pure
+        IPW's covariate adjustment is a single scalar level shift, so it cannot
+        adjust the dose-response *shape* (ACRT(d) would be identical to the
+        unconditional fit); it raises ``NotImplementedError``. ``reg`` and ``dr``
+        share the dose-response shape and ACRT(d); ``dr`` differs only in the
+        ``overall_att`` / ATT(d) level and in its doubly-robust standard errors.
+    pscore_trim : float, default=0.01
+        Propensity-score trimming bound for the ``dr`` path (scores clipped to
+        ``[pscore_trim, 1 - pscore_trim]``).
+    epv_threshold : float, default=10.0
+        Events-per-variable threshold for the ``dr`` propensity logit diagnostics.
+    pscore_fallback : str, default="error"
+        Action when ``dr`` propensity estimation raises (the logit IRLS fails
+        with a ``LinAlgError`` / ``ValueError``, e.g. perfect separation or rank
+        deficiency): ``"error"`` (re-raise — the default, fail-closed so a `dr`
+        fit never silently degrades to a non-DR estimate) or ``"unconditional"``
+        (fall back to an unconditional propensity with a warning; the affected
+        cells are then reg-like — use only when you knowingly accept that). Note:
+        low events-per-variable emits a diagnostic warning but does not itself
+        trigger the fallback.
 
     Examples
     --------
@@ -86,6 +116,7 @@ class ContinuousDiD:
 
     _VALID_CONTROL_GROUPS = {"never_treated", "not_yet_treated"}
     _VALID_BASE_PERIODS = {"varying", "universal"}
+    _VALID_ESTIMATION_METHODS = {"reg", "dr", "ipw"}
 
     def __init__(
         self,
@@ -100,6 +131,11 @@ def __init__(
         bootstrap_weights: str = "rademacher",
         seed: Optional[int] = None,
         rank_deficient_action: str = "warn",
+        covariates: Optional[List[str]] = None,
+        estimation_method: str = "dr",
+        pscore_trim: float = 0.01,
+        epv_threshold: float = 10.0,
+        pscore_fallback: str = "error",
     ):
         self.degree = degree
         self.num_knots = num_knots
@@ -112,10 +148,15 @@ def __init__(
         self.bootstrap_weights = bootstrap_weights
         self.seed = seed
         self.rank_deficient_action = rank_deficient_action
+        self.covariates = covariates
+        self.estimation_method = estimation_method
+        self.pscore_trim = pscore_trim
+        self.epv_threshold = epv_threshold
+        self.pscore_fallback = pscore_fallback
         self._validate_constrained_params()
 
     def _validate_constrained_params(self) -> None:
-        """Validate control_group and base_period values."""
+        """Validate control_group, base_period, and estimation_method values."""
         if self.control_group not in self._VALID_CONTROL_GROUPS:
             raise ValueError(
                 f"Invalid control_group: '{self.control_group}'. "
@@ -126,6 +167,24 @@ def _validate_constrained_params(self) -> None:
                 f"Invalid base_period: '{self.base_period}'. "
                 f"Must be one of {self._VALID_BASE_PERIODS}."
             )
+        if self.estimation_method not in self._VALID_ESTIMATION_METHODS:
+            raise ValueError(
+                f"Invalid estimation_method: '{self.estimation_method}'. "
+                f"Must be one of {self._VALID_ESTIMATION_METHODS}."
+            )
+        if self.pscore_fallback not in {"unconditional", "error"}:
+            raise ValueError(
+                f"Invalid pscore_fallback: '{self.pscore_fallback}'. "
+                "Must be 'unconditional' or 'error'."
+            )
+        if not (np.isfinite(self.pscore_trim) and 0.0 <= self.pscore_trim < 0.5):
+            raise ValueError(
+                f"Invalid pscore_trim: {self.pscore_trim}. " "Must be finite and in [0, 0.5)."
+            )
+        if not (np.isfinite(self.epv_threshold) and self.epv_threshold > 0):
+            raise ValueError(
+                f"Invalid epv_threshold: {self.epv_threshold}. Must be finite and > 0."
+            )
 
     def get_params(self) -> Dict[str, Any]:
         """Return estimator parameters as a dictionary."""
@@ -141,15 +200,29 @@ def get_params(self) -> Dict[str, Any]:
             "bootstrap_weights": self.bootstrap_weights,
             "seed": self.seed,
             "rank_deficient_action": self.rank_deficient_action,
+            "covariates": self.covariates,
+            "estimation_method": self.estimation_method,
+            "pscore_trim": self.pscore_trim,
+            "epv_threshold": self.epv_threshold,
+            "pscore_fallback": self.pscore_fallback,
         }
 
     def set_params(self, **params) -> "ContinuousDiD":
-        """Set estimator parameters and return self."""
-        for key, value in params.items():
+        """Set estimator parameters and return self.
+
+        Transactional: the candidate configuration is validated against a clone
+        before any attribute on ``self`` is mutated, so an invalid update leaves
+        the estimator unchanged (no half-applied state).
+        """
+        for key in params:
             if not hasattr(self, key):
                 raise ValueError(f"Invalid parameter: {key}")
+        # Validate the candidate config on a throwaway before mutating self.
+        candidate = ContinuousDiD(**{**self.get_params(), **params})
+        # candidate.__init__ ran _validate_constrained_params() successfully.
+        del candidate
+        for key, value in params.items():
             setattr(self, key, value)
-        self._validate_constrained_params()
         return self
 
     # ------------------------------------------------------------------
@@ -215,10 +288,42 @@ def fit(
         # Bootstrap + survey supported via PSU-level multiplier bootstrap.
 
         df = data.copy()
-        for col in [outcome, unit, time, first_treat, dose]:
+        cov_cols = list(self.covariates) if self.covariates else []
+        for col in [outcome, unit, time, first_treat, dose, *cov_cols]:
             if col not in df.columns:
                 raise ValueError(f"Column '{col}' not found in data.")
 
+        # Covariate-path guards (conditional parallel trends).
+        if cov_cols:
+            if survey_design is not None:
+                raise NotImplementedError(
+                    "ContinuousDiD does not yet support covariates= together with "
+                    "survey_design= (weighted covariate outcome-regression / "
+                    "propensity influence functions are a follow-up). Use one or "
+                    "the other for now."
+                )
+            if self.estimation_method == "ipw":
+                raise NotImplementedError(
+                    "estimation_method='ipw' is not supported with covariates on the "
+                    "dose-response / event-study aggregation. Pure IPW's covariate "
+                    "adjustment is a single scalar (a propensity-reweighted control "
+                    "mean), which shifts only the ATT(d) level and leaves ACRT(d) "
+                    "identical to the unconditional fit — it cannot adjust the "
+                    "dose-response shape. Use estimation_method='reg' or 'dr'."
+                )
+            # Fail closed on missing/non-finite covariates: a per-cell fallback to
+            # unconditional estimation would silently mix conditional-PT and
+            # unconditional-PT cells in the aggregate (no-silent-failures).
+            cov_nonfinite = ~np.isfinite(df[cov_cols].to_numpy(dtype=float))
+            if cov_nonfinite.any():
+                n_bad = int(cov_nonfinite.any(axis=1).sum())
+                raise ValueError(
+                    f"{n_bad} row(s) have missing/non-finite covariate values. "
+                    "ContinuousDiD requires complete covariates (a per-cell fallback "
+                    "would mix conditional and unconditional estimands). Drop or "
+                    "impute the affected rows, or fit without covariates."
+                )
+
         # Verify dose is time-invariant
         dose_nunique = df.groupby(unit)[dose].nunique()
         if dose_nunique.max() > 1:
@@ -377,6 +482,7 @@ def fit(
             dose,
             time_periods,
             survey_weights=survey_weights,
+            covariates=self.covariates,
         )
 
         # Compute dvals (evaluation grid)
@@ -658,6 +764,16 @@ def fit(
                             if not b_info:
                                 continue
                             w = ws[idx_cell]
+                            # Covariate path: the binarized event-study effect is
+                            # att_glob, whose per-unit cell IF is precomputed.
+                            cov_if = b_info.get("cov_if")
+                            if cov_if is not None:
+                                np.add.at(
+                                    if_es,
+                                    cov_if["cell_indices"],
+                                    w * cov_if["if_att_glob"],
+                                )
+                                continue
                             treated_idx = b_info["treated_indices"]
                             control_idx = b_info["control_indices"]
                             n_t = b_info["n_treated"]
@@ -768,6 +884,11 @@ def fit(
             n_control_units=n_control,
             alpha=self.alpha,
             control_group=self.control_group,
+            covariates=self.covariates,
+            estimation_method=self.estimation_method,
+            pscore_trim=self.pscore_trim,
+            epv_threshold=self.epv_threshold,
+            pscore_fallback=self.pscore_fallback,
             degree=self.degree,
             num_knots=self.num_knots,
             base_period=self.base_period,
@@ -794,6 +915,7 @@ def _precompute_structures(
         dose: str,
         time_periods: List[Any],
         survey_weights: Optional[np.ndarray] = None,
+        covariates: Optional[List[str]] = None,
     ) -> Dict[str, Any]:
         """Pivot to wide format and build lookup structures."""
         all_units = sorted(df[unit].unique())
@@ -809,6 +931,17 @@ def _precompute_structures(
             j = period_to_col[row[time]]
             outcome_matrix[i, j] = row[outcome]
 
+        # Covariate cube: {period -> (n_units, n_cov)} read from the given
+        # period (the cell reads its base period). Covariates may be
+        # time-varying; the cell uses the base-period slice.
+        covariate_by_period = None
+        if covariates:
+            cov_cube = np.full((n_units, n_periods, len(covariates)), np.nan)
+            ui = df[unit].map(unit_to_idx).to_numpy()
+            ti = df[time].map(period_to_col).to_numpy()
+            cov_cube[ui, ti, :] = df[list(covariates)].to_numpy(dtype=float)
+            covariate_by_period = {t: cov_cube[:, period_to_col[t], :] for t in time_periods}
+
         # Per-unit cohort and dose
         unit_cohorts = np.zeros(n_units, dtype=float)
         dose_vector = np.zeros(n_units, dtype=float)
@@ -850,6 +983,211 @@ def _precompute_structures(
             "n_units": n_units,
             "unit_survey_weights": unit_survey_weights,
             "unit_first_panel_row": unit_first_panel_row,
+            "covariate_by_period": covariate_by_period,
+        }
+
+    # ------------------------------------------------------------------
+    # Covariate adjustment (conditional parallel trends)
+    # ------------------------------------------------------------------
+
+    def _fit_covariate_adjustment(
+        self,
+        delta_y_treated: np.ndarray,
+        delta_y_control: np.ndarray,
+        X_treated_raw: np.ndarray,
+        X_control_raw: np.ndarray,
+        g: Any,
+        t: Any,
+    ) -> Optional[Dict[str, Any]]:
+        """Covariate-adjusted control counterfactual for one (g,t) cell.
+
+        Returns the per-treated counterfactual ``mu_0_vec`` and the nuisance
+        pieces the influence function needs (reg: outcome-regression only; dr:
+        + propensity and the DRDID doubly-robust per-unit IF). Missing/non-finite
+        covariate values raise ``ValueError`` (fail-closed; ``fit()`` also rejects
+        them up front — a per-cell fallback would mix conditional and
+        unconditional estimands). ``reg`` and ``dr`` share the same OLS outcome
+        regression on controls; ``dr`` adds a propensity model and a scalar
+        augmentation ``eta_cont`` (a level term).
+        """
+        if np.any(np.isnan(X_treated_raw)) or np.any(np.isnan(X_control_raw)):
+            # Defensive: fit() rejects non-finite covariates up front, so this
+            # is an internal-invariant guard (no silent per-cell fallback).
+            raise ValueError(
+                f"Missing covariate values reached cell (g={g}, t={t}). "
+                "ContinuousDiD requires complete covariates."
+            )
+
+        n_t = len(delta_y_treated)
+        n_c = len(delta_y_control)
+        Xt = np.column_stack([np.ones(n_t), X_treated_raw])
+        Xc = np.column_stack([np.ones(n_c), X_control_raw])
+
+        # Outcome regression on controls (OLS) — shared by reg and dr.
+        gamma, _, _ = solve_ols(
+            Xc,
+            delta_y_control,
+            return_vcov=False,
+            rank_deficient_action=self.rank_deficient_action,
+        )
+        gamma = np.where(np.isnan(gamma), 0.0, gamma)
+        mu_treated = Xt @ gamma
+        or_resid_c = delta_y_control - Xc @ gamma  # per-control OR residual
+
+        # OR-nuisance influence function: IF_gamma[k] = M_c^{-1} X_c[k] resid_c[k]
+        Mc = (Xc.T @ Xc) / n_c
+        Mc_inv, _, _ = _rank_guarded_inv(Mc)
+        IF_gamma = (Xc @ Mc_inv) * or_resid_c[:, np.newaxis]  # (n_c, p), Mc_inv symmetric
+        x_bar_treated = Xt.mean(axis=0)
+        n_total = n_t + n_c
+
+        eta_cont = 0.0
+        dr_inf = None
+        if self.estimation_method == "dr":
+            D = np.concatenate([np.ones(n_t), np.zeros(n_c)])
+            X_raw_all = np.vstack([X_treated_raw, X_control_raw])
+            try:
+                _, ps = solve_logit(
+                    X_raw_all,
+                    D,
+                    rank_deficient_action=self.rank_deficient_action,
+                    epv_threshold=self.epv_threshold,
+                )
+            except (np.linalg.LinAlgError, ValueError):
+                # Fail closed by default; also honor an explicit
+                # rank_deficient_action="error" request (mirrors CS).
+                if self.pscore_fallback == "error" or self.rank_deficient_action == "error":
+                    raise
+                warnings.warn(
+                    f"Propensity score estimation failed for (g={g}, t={t}); "
+                    "falling back to an unconditional propensity for this cell "
+                    "(this cell is now reg-like, not doubly-robust). "
+                    "Consider estimation_method='reg' to avoid propensity scores.",
+                    UserWarning,
+                    stacklevel=3,
+                )
+                ps = np.full(n_total, n_t / n_total)
+            ps = np.clip(ps, self.pscore_trim, 1 - self.pscore_trim)
+            odds_c = ps[n_t:] / (1 - ps[n_t:])
+            eta_cont = float(np.sum(odds_c * or_resid_c) / np.sum(odds_c))
+            dY_all = np.concatenate([delta_y_treated, delta_y_control])
+            dr_inf = self._dr_cell_inf_func(dY_all, D, np.vstack([Xt, Xc]), gamma, ps)
+
+        return {
+            "mu_0_vec": mu_treated + eta_cont,  # per-treated counterfactual
+            "Xt": Xt,  # (n_t, p) treated covariates incl. intercept
+            "IF_gamma": IF_gamma,  # (n_c, p)
+            "x_bar_treated": x_bar_treated,  # (p,)
+            "eta_cont": eta_cont,
+            "dr_inf": dr_inf,  # (n_total,) DRDID att.inf.func, treated-then-control
+            "n_total": n_total,
+        }
+
+    def _dr_cell_inf_func(
+        self,
+        dY: np.ndarray,
+        D: np.ndarray,
+        X: np.ndarray,
+        gamma: np.ndarray,
+        ps: np.ndarray,
+    ) -> np.ndarray:
+        """DRDID ``drdid_panel`` doubly-robust per-unit influence function.
+
+        Direct port of ``DRDID::drdid_panel$att.inf.func`` (unit weights = 1,
+        propensity already clipped so no drop-trimming). Units are ordered
+        treated-then-control. Validated to ~1e-13 against DRDID in the spike.
+        """
+        n = len(D)
+        out_delta = X @ gamma
+        w_treat = D
+        w_cont = ps * (1 - D) / (1 - ps)
+        dr_treat = w_treat * (dY - out_delta)
+        dr_cont = w_cont * (dY - out_delta)
+        eta_treat = dr_treat.mean() / w_treat.mean()
+        eta_cont = dr_cont.mean() / w_cont.mean()
+        weights_ols = 1.0 - D
+        wols_eX = (weights_ols * (dY - out_delta))[:, np.newaxis] * X
+        XpX = ((weights_ols[:, np.newaxis] * X).T @ X) / n
+        XpX_inv, _, _ = _rank_guarded_inv(XpX)
+        asy_wols = wols_eX @ XpX_inv
+        W = ps * (1 - ps)
+        score_ps = (D - ps)[:, np.newaxis] * X
+        Hess, _, _ = _rank_guarded_inv((X.T @ (W[:, np.newaxis] * X)) / n)
+        asy_ps = score_ps @ Hess
+        inf_treat_1 = dr_treat - w_treat * eta_treat
+        M1 = (w_treat[:, np.newaxis] * X).mean(axis=0)
+        inf_treat = (inf_treat_1 - asy_wols @ M1) / w_treat.mean()
+        inf_cont_1 = dr_cont - w_cont * eta_cont
+        M2 = (w_cont[:, np.newaxis] * (dY - out_delta - eta_cont)[:, np.newaxis] * X).mean(axis=0)
+        M3 = (w_cont[:, np.newaxis] * X).mean(axis=0)
+        inf_control = (inf_cont_1 + asy_ps @ M2 - asy_wols @ M3) / w_cont.mean()
+        return inf_treat - inf_control
+
+    def _covariate_cell_influence(
+        self,
+        cov_adj: Dict[str, Any],
+        delta_tilde_y: np.ndarray,
+        att_glob: float,
+        residuals: np.ndarray,
+        Psi: np.ndarray,
+        bread: np.ndarray,
+        Psi_eval: np.ndarray,
+        dPsi_eval: np.ndarray,
+        dpsi_bar: np.ndarray,
+        treated_indices: np.ndarray,
+        control_indices: np.ndarray,
+        n_t: int,
+        n_c: int,
+    ) -> Dict[str, Any]:
+        """Per-cell-unit influence functions for the covariate-adjusted estimands.
+
+        Returns arrays aligned to ``cell_indices = [treated_indices,
+        control_indices]`` in the 1/n convention (the aggregator scatters them
+        with the cell weight and the SE is ``sqrt(sum(IF^2))``). At p=1 (intercept
+        only) this reduces exactly to the unconditional control-IF path.
+        """
+        Xt = cov_adj["Xt"]  # (n_t, p)
+        IF_gamma = cov_adj["IF_gamma"]  # (n_c, p)
+        x_bar_treated = cov_adj["x_bar_treated"]  # (p,)
+
+        # Treated: beta perturbation from the B-spline residual score.
+        beta_if_t = (Psi * residuals[:, np.newaxis]) @ bread / n_t  # (n_t, K), bread symmetric
+        att_d_if_t = beta_if_t @ Psi_eval.T  # (n_t, n_grid)
+        acrt_d_if_t = beta_if_t @ dPsi_eval.T
+        acrt_glob_if_t = beta_if_t @ dpsi_bar  # (n_t,)
+        att_glob_if_t = (delta_tilde_y - att_glob) / n_t  # (n_t,)
+
+        # Control: enters through the outcome-regression nuisance gamma_hat.
+        # E_T[Psi X'] (K x p) replaces the scalar psi_bar; E_T[X'] the scalar 1.
+        Psi_X_bar = (Psi.T @ Xt) / n_t  # (K, p)
+        beta_if_c = -((IF_gamma @ Psi_X_bar.T) @ bread) / n_c  # (n_c, K)
+        att_d_if_c = beta_if_c @ Psi_eval.T
+        acrt_d_if_c = beta_if_c @ dPsi_eval.T
+        acrt_glob_if_c = beta_if_c @ dpsi_bar
+        att_glob_if_c = -(IF_gamma @ x_bar_treated) / n_c  # (n_c,)
+
+        if_att_glob = np.concatenate([att_glob_if_t, att_glob_if_c])
+        if_acrt_glob = np.concatenate([acrt_glob_if_t, acrt_glob_if_c])
+        if_att_d = np.vstack([att_d_if_t, att_d_if_c])  # (n_total, n_grid)
+        if_acrt_d = np.vstack([acrt_d_if_t, acrt_d_if_c])
+
+        if self.estimation_method == "dr":
+            # dr shares the reg curve shape (ACRT identical); att_glob / ATT(d)
+            # differ by the augmentation eta_cont. Ground att_glob's IF in the
+            # validated DRDID doubly-robust IF and shift ATT(d) uniformly by the
+            # augmentation IF (= reg att_glob IF - dr att_glob IF). ACRT untouched.
+            n_total = cov_adj["n_total"]
+            dr_att_glob_if = cov_adj["dr_inf"] / n_total  # (n_total,)
+            if_eta = if_att_glob - dr_att_glob_if  # augmentation IF, per unit
+            if_att_d = if_att_d - if_eta[:, np.newaxis]
+            if_att_glob = dr_att_glob_if
+
+        return {
+            "cell_indices": np.concatenate([treated_indices, control_indices]),
+            "if_att_glob": if_att_glob,
+            "if_acrt_glob": if_acrt_glob,
+            "if_att_d": if_att_d,
+            "if_acrt_d": if_acrt_d,
         }
 
     def _compute_dose_response_gt(
@@ -940,14 +1278,34 @@ def _compute_dose_response_gt(
                     "acrt_d": np.full(len(dvals), np.nan),
                 }
 
-        # Control counterfactual (weighted mean when survey weights present)
-        if w_control is not None:
+        # Control counterfactual.
+        #   - No covariates: scalar control mean mu_0 (unconditional PT).
+        #   - Covariates (reg/dr): per-treated covariate-adjusted counterfactual
+        #     mu_0_vec = X_i'gamma_hat (+ scalar DR augmentation eta_cont for dr),
+        #     under conditional PT. `cov_adj` carries the nuisance pieces the
+        #     influence function needs. Survey weights are rejected upstream on
+        #     the covariate path, so w_control/w_treated are None here.
+        covariate_by_period = precomp.get("covariate_by_period")
+        cov_adj = None
+        if covariate_by_period is not None:
+            cov_adj = self._fit_covariate_adjustment(
+                delta_y_treated,
+                delta_y_control,
+                covariate_by_period[base_t][treated_mask],
+                covariate_by_period[base_t][control_mask],
+                g,
+                t,
+            )
+
+        if cov_adj is not None:
+            mu_0 = float(np.mean(delta_y_control))  # retained for metadata only
+            delta_tilde_y = delta_y_treated - cov_adj["mu_0_vec"]
+        elif w_control is not None:
             mu_0 = float(np.average(delta_y_control, weights=w_control))
+            delta_tilde_y = delta_y_treated - mu_0
         else:
             mu_0 = float(np.mean(delta_y_control))
-
-        # Demean
-        delta_tilde_y = delta_y_treated - mu_0
+            delta_tilde_y = delta_y_treated - mu_0
 
         # Treated doses
         treated_doses = dose_vector[treated_mask]
@@ -1015,8 +1373,12 @@ def _compute_dose_response_gt(
         att_d = Psi_eval @ beta_pred
         acrt_d = dPsi_eval @ beta_pred
 
-        # Summary parameters
-        if w_treated is not None:
+        # Summary parameters. With covariates, att_glob = mean_T(delta_tilde_y)
+        # (reg: mean_T(dY - X'gamma); dr: additionally minus the augmentation
+        # eta_cont, already folded into delta_tilde_y via mu_0_vec).
+        if cov_adj is not None:
+            att_glob = float(np.mean(delta_tilde_y))
+        elif w_treated is not None:
             att_glob = float(np.average(delta_y_treated, weights=w_treated) - mu_0)
         else:
             att_glob = float(np.mean(delta_y_treated) - mu_0)
@@ -1090,6 +1452,29 @@ def _compute_dose_response_gt(
         else:
             dpsi_bar = np.mean(dPsi_treated, axis=0)
 
+        # Covariate influence-function arrays (per cell unit, treated-then-control,
+        # in the 1/n "sum-of-squares SE" convention). Reg uses the p-vector OR
+        # generalization of the unconditional control IF (which is its p=1 case);
+        # dr reuses the reg curve shape and grounds att_glob + the augmentation IF
+        # in the validated DRDID doubly-robust per-unit IF. See REGISTRY.
+        cov_if = None
+        if cov_adj is not None:
+            cov_if = self._covariate_cell_influence(
+                cov_adj,
+                delta_tilde_y,
+                att_glob,
+                residuals,
+                Psi,
+                bread,
+                Psi_eval,
+                dPsi_eval,
+                dpsi_bar,
+                treated_indices,
+                control_indices,
+                n_treated,
+                n_control,
+            )
+
         bootstrap_info = {
             "bread": bread,
             "ee_treated": ee_treated,
@@ -1110,6 +1495,7 @@ def _compute_dose_response_gt(
             "mu_0": mu_0,
             "att_glob": att_glob,
             "acrt_glob": acrt_glob,
+            "cov_if": cov_if,
         }
 
         # Store survey-weighted masses and per-unit arrays for IF linearization
@@ -1205,6 +1591,16 @@ def _compute_analytical_se(
             info = gt_bootstrap_info[gt]
             if not info:
                 continue
+            # Covariate path: scatter the pre-computed per-unit cell IFs
+            # (unconditional / survey path is untouched below).
+            cov_if = info.get("cov_if")
+            if cov_if is not None:
+                idx = cov_if["cell_indices"]
+                np.add.at(if_att_glob, idx, w * cov_if["if_att_glob"])
+                np.add.at(if_acrt_glob, idx, w * cov_if["if_acrt_glob"])
+                np.add.at(if_att_d, idx, w * cov_if["if_att_d"])
+                np.add.at(if_acrt_d, idx, w * cov_if["if_acrt_d"])
+                continue
             treated_idx = info["treated_indices"]
             control_idx = info["control_indices"]
             n_t = info["n_treated"]
@@ -1485,6 +1881,19 @@ def _run_bootstrap(
         # Helper to bootstrap a single (g,t) cell
         def _bootstrap_gt_cell(gt, info):
             """Returns att_glob_b array (B,) for this cell."""
+            # Covariate path: the multiplier bootstrap perturbs the same per-unit
+            # cell influence functions the analytical SE uses.
+            cov_if = info.get("cov_if")
+            if cov_if is not None:
+                xi = all_weights[:, cov_if["cell_indices"]]  # (B, n_cell)
+                beta_pred_c = info["beta_pred"]
+                cell_att_d = info["Psi_eval"] @ beta_pred_c
+                cell_acrt_d = info["dPsi_eval"] @ beta_pred_c
+                att_d_b = cell_att_d[np.newaxis, :] + xi @ cov_if["if_att_d"]
+                acrt_d_b = cell_acrt_d[np.newaxis, :] + xi @ cov_if["if_acrt_d"]
+                att_glob_b = info["att_glob"] + xi @ cov_if["if_att_glob"]
+                acrt_glob_b = xi @ cov_if["if_acrt_glob"]
+                return att_d_b, acrt_d_b, att_glob_b, acrt_glob_b, info.get("acrt_glob", 0.0)
             treated_idx = info["treated_indices"]
             control_idx = info["control_indices"]
             n_t = info["n_treated"]
diff --git a/diff_diff/continuous_did_results.py b/diff_diff/continuous_did_results.py
index f60ea5ead..9a5bad1ef 100644
--- a/diff_diff/continuous_did_results.py
+++ b/diff_diff/continuous_did_results.py
@@ -139,6 +139,14 @@ class ContinuousDiDResults:
     bootstrap_weights: str = "rademacher"
     seed: Optional[int] = None
     rank_deficient_action: str = "warn"
+    # Covariate adjustment (conditional parallel trends). ``covariates`` is None
+    # for the unconditional path; ``estimation_method`` is only meaningful when
+    # covariates are used (``"reg"`` or ``"dr"``).
+    covariates: Optional[List[str]] = field(default=None)
+    estimation_method: str = "dr"
+    pscore_trim: float = 0.01
+    epv_threshold: float = 10.0
+    pscore_fallback: str = "error"
     event_study_effects: Optional[Dict[int, Dict[str, Any]]] = field(default=None)
     # Survey design metadata (SurveyMetadata instance from diff_diff.survey)
     survey_metadata: Optional[Any] = field(default=None)
@@ -224,8 +232,11 @@ def summary(self, alpha: Optional[float] = None) -> str:
             f"{'Interior knots:':<30} {self.num_knots:>10}",
             f"{'Base period:':<30} {self.base_period:>10}",
             f"{'Anticipation:':<30} {self.anticipation:>10}",
-            "",
         ]
+        if self.covariates:
+            lines.append(f"{'Covariates:':<30} {', '.join(self.covariates):>10}")
+            lines.append(f"{'Estimation method:':<30} {self.estimation_method:>10}")
+        lines.append("")
 
         # Add survey design info
         if self.survey_metadata is not None:
diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt
index 08193daba..7d3177b05 100644
--- a/diff_diff/guides/llms-full.txt
+++ b/diff_diff/guides/llms-full.txt
@@ -696,9 +696,21 @@ ContinuousDiD(
     bootstrap_weights: str = "rademacher",
     seed: int | None = None,
     rank_deficient_action: str = "warn",
+    covariates: list[str] | None = None,      # Conditional parallel trends (X). None = unconditional
+    estimation_method: str = "dr",            # "reg" or "dr" (used only with covariates); "ipw" not
+                                              #   supported on the dose curve (raises NotImplementedError)
+    pscore_trim: float = 0.01,                # dr propensity trimming
+    epv_threshold: float = 10.0,              # dr propensity events-per-variable
+    pscore_fallback: str = "error",           # dr propensity-failure action: "error" (fail-closed
+                                              #   default) or "unconditional" (opt-in reg-like fallback)
 )
 ```
 
+Covariates give conditional parallel trends: each (g,t) cell's control counterfactual becomes a
+covariate-adjusted prediction. `reg`/`dr` share the ATT(d) *shape* and ACRT(d); `dr` (doubly-robust,
+default) differs only in the overall_att/ATT(d) level and SE. `overall_att`+SE match DRDID
+reg_did_panel/drdid_panel. `covariates=` + `survey_design=` is not yet supported.
+
 **Alias:** `CDiD`
 
 **fit() parameters:**
diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
index a9d07b7a5..bd05cdaf7 100644
--- a/docs/methodology/REGISTRY.md
+++ b/docs/methodology/REGISTRY.md
@@ -889,6 +889,12 @@ No selection into dose groups on the basis of treatment effects.
 Implies `ATT(d|d) = ATT(d)` for all d.
 Additionally identifies: `ATT(d)`, `ACRT(d)`, `ACRT^{glob}`, and cross-dose comparisons.
 
+**Conditional Parallel Trends (with covariates).** When `covariates=` is passed the PT/SPT
+assumptions are conditional on covariates `X`:
+`E[Y_t(0) - Y_{t-1}(0) | D = d, X] = E[Y_t(0) - Y_{t-1}(0) | D = 0, X]`. The per-`(g,t)` cell's
+control counterfactual becomes a covariate-adjusted prediction instead of the unconditional control
+mean (see Key Equations and the covariate Note below).
+
 See `docs/methodology/continuous-did.md` Section 4 for full details.
 
 ### Key Equations
@@ -907,6 +913,17 @@ See `docs/methodology/continuous-did.md` Section 4 for full details.
 3. OLS: `beta = (Psi'Psi)^{-1} Psi' Delta_tilde_Y`
 4. `ATT(d) = Psi(d)' beta`, `ACRT(d) = dPsi(d)/dd' beta`
 
+**Covariate-adjusted `Delta_tilde_Y` (conditional PT).** With `covariates=`, step 1's scalar control
+mean is replaced by a per-treated-unit covariate-adjusted counterfactual (`X_i` includes an intercept):
+- `reg` (outcome regression): fit `gamma_hat = (X_C'X_C)^{-1} X_C' Delta_Y_C` on controls;
+  `Delta_tilde_Y_i = Delta_Y_i - X_i' gamma_hat`.
+- `dr` (doubly-robust, DRDID `drdid_panel`): same OLS `gamma_hat`, plus a propensity model and a scalar
+  augmentation `eta_cont = odds_weighted_mean_C(Delta_Y - X' gamma_hat)`;
+  `Delta_tilde_Y_i = Delta_Y_i - X_i' gamma_hat - eta_cont`.
+Steps 2-4 are unchanged. Because the augmentation is a constant, `reg` and `dr` share the same
+`ACRT(d)` (a constant only shifts the B-spline intercept, which `dPsi` annihilates); they differ only
+in the `ATT(d)` / `ATT^{glob}` level (by `-eta_cont`) and in the doubly-robust SE.
+
 ### Edge Cases
 
 - **No untreated group**: Remark 3.1 (lowest-dose-as-control) not implemented; requires P(D=0) > 0.
@@ -939,6 +956,7 @@ labels.*
 2. **Note:** `bspline_derivative_design_matrix` derivative-failure `UserWarning` — Phase 2 axis-C #12 silent-failures audit fix. No R correspondence; `contdid` v0.1.0 does not implement an equivalent warning. Cross-references the § Edge Cases `**Note:**` bullet above (`bspline_derivative_design_matrix` entry) and `METHODOLOGY_REVIEW.md` § ContinuousDiD Deviations #2. Locked in `tests/test_continuous_did.py::TestBSplineDerivativeDegenerateBasis` (3 tests); source-level aggregate-warning block at `diff_diff/continuous_did_bspline.py:150-187`.
 3. **Note:** `+inf` → `0` never-treated recoding emits `UserWarning` reporting the affected row count; negative `first_treat` (including `-inf`) raises `ValueError`. Axis-E silent-coercion fix per Phase 2 audit. No R correspondence; `contdid` v0.1.0 silently absorbs `+inf` without a signal. Cross-references the § Implementation Checklist `**Note:**` below and `METHODOLOGY_REVIEW.md` § ContinuousDiD Deviations #3.
 4. **Note:** Zero-`first_treat` rows with nonzero `dose` are force-zeroed with `UserWarning` reporting the affected row count (axis-E silent-coercion). No R correspondence; `contdid` v0.1.0 has the same `first_treat = 0` → `D = 0` invariant but silently coerces without a warning. Cross-references the § Implementation Checklist `**Note:**` below and `METHODOLOGY_REVIEW.md` § ContinuousDiD Deviations #4.
+5. **Note (covariate support — library extension beyond `contdid` v0.1.0):** `covariates=` with `estimation_method ∈ {reg, dr}` adds conditional-parallel-trends adjustment. This is a **library extension**: `contdid` v0.1.0 hard-stops on any covariate (`stop("covariates not currently supported…")`), so there is **no external R anchor for the covariate-adjusted dose *curve***. Validation instead: (a) the **scalar `overall_att` + SE** map *exactly* onto `DRDID::reg_did_panel` (reg) / `DRDID::drdid_panel` (dr) — a tight (~1e-8) component anchor, skip-guarded since DRDID is not in CI (`tests/test_methodology_continuous_did.py::TestCovariateReg`); (b) an **R-free NumPy reconstruction** of the reg/dr `att`+SE runs *in CI* at p≥2 (`test_dr_reg_numpy_crosscheck_p2`) — the guard the p=1 reduction cannot provide (at p=1 the intercept-only propensity is constant, so `eta_cont ≡ 0` and dr collapses to reg); (c) DGP recovery + MC coverage (reg 96%, dr 95%). **`ipw` restricted:** `estimation_method="ipw"` with covariates raises `NotImplementedError` — pure IPW's covariate adjustment is a single scalar (a propensity-reweighted control mean) that shifts only the ATT(d) level and leaves `ACRT(d)` identical to the unconditional fit, so it cannot adjust the dose-response *shape*. **Deviations from DRDID:** unit weights are 1 (unweighted; `covariates=` + `survey_design=` raises `NotImplementedError`, deferred); propensity trimming uses clip semantics (`pscore_trim`) rather than DRDID's drop-trimming — identical on moderate-overlap data (the anchor regime), diverging only at extreme propensities. **Fail-closed policies (no-silent-failures):** (i) missing/non-finite covariate values raise `ValueError` up front — a per-cell fallback to unconditional estimation would silently mix conditional-PT and unconditional-PT cells in the aggregate; (ii) `dr` propensity-estimation failure raises by default (`pscore_fallback="error"`) so a `dr` fit never silently degrades to a non-DR estimate — `pscore_fallback="unconditional"` opts into the graceful (warned, reg-like) fallback. Cross-references `docs/methodology/continuous-did.md` § Covariates.
 
 ### Implementation Checklist
 
@@ -948,7 +966,7 @@ labels.*
 - [x] Multiplier bootstrap for inference
 - [x] Analytical SEs via influence functions
 - [x] Equation verification tests (linear, quadratic, multi-period)
-- [ ] Covariate support (deferred, matching R v0.1.0)
+- [x] Covariate support (reg / dr) — conditional parallel trends. **ipw restricted** (see Note below); survey × covariate deferred; discrete-saturated & Remark 3.1 still deferred.
 - [ ] Discrete treatment saturated regression
 - [ ] Lowest-dose-as-control (Remark 3.1)
 - [x] Survey design support (Phase 3): weighted B-spline OLS, TSL on influence functions; bootstrap+survey supported (Phase 6)
diff --git a/docs/methodology/continuous-did.md b/docs/methodology/continuous-did.md
index c1db94342..76d8378a1 100644
--- a/docs/methodology/continuous-did.md
+++ b/docs/methodology/continuous-did.md
@@ -434,7 +434,14 @@ can be meaningfully aggregated.
 ### Phase 3 (Advanced)
 8. CCK nonparametric estimation
 9. Uniform confidence bands
-10. Covariates support (DR/IPW/OR)
+10. Covariates support — **implemented** for outcome-regression (`reg`) and doubly-robust (`dr`)
+   under conditional parallel trends (`covariates=`, `estimation_method=`). Each `(g,t)` cell replaces
+   the unconditional control mean with a covariate-adjusted counterfactual (`reg`:
+   `ΔY_i − X_i'γ̂`; `dr`: additionally minus the DRDID augmentation `η̄_cont`); the B-spline dose
+   layer is unchanged. Scalar `overall_att` + SE match `DRDID::reg_did_panel` / `drdid_panel`. `ipw`
+   is not offered on the dose curve (its covariate adjustment is a scalar level shift → `ACRT(d)`
+   unchanged), and `covariates=` + `survey_design=` is deferred. See REGISTRY § ContinuousDiD
+   Note #5.
 
 ### Defer
 - Discrete treatment (saturated regression — simpler, add later)
diff --git a/tests/test_continuous_did.py b/tests/test_continuous_did.py
index 4cffb060d..ba7c35db4 100644
--- a/tests/test_continuous_did.py
+++ b/tests/test_continuous_did.py
@@ -1559,8 +1559,165 @@ def force_drop(A, **kwargs):
                     data, "outcome", "unit", "period", "first_treat", "dose"
                 )
         msgs = [str(w.message) for w in caught]
-        assert any(
-            "ContinuousDiD ACRT variance" in m and "rank-deficient" in m for m in msgs
-        ), msgs
+        assert any("ContinuousDiD ACRT variance" in m and "rank-deficient" in m for m in msgs), msgs
         # rank-reduced bread still yields finite ACRT SEs.
         assert np.all(np.isfinite(results.dose_response_acrt.se))
+
+
+# =============================================================================
+# Covariate adjustment — API, params, and fail-closed behavior
+# =============================================================================
+
+
+def _cov_data(seed=5, n_units=120):
+    """Continuous-DiD panel with one covariate column added."""
+    data = generate_continuous_did_data(n_units=n_units, n_periods=3, seed=seed, noise_sd=0.5)
+    rng = np.random.default_rng(seed)
+    # time-invariant covariate (one value per unit)
+    uc = data.groupby("unit").ngroup().to_numpy()
+    per_unit = rng.normal(size=data["unit"].nunique())
+    data["x1"] = per_unit[uc]
+    return data
+
+
+class TestCovariateAPI:
+    def test_covariates_and_method_in_params(self):
+        est = ContinuousDiD(covariates=["x1"], estimation_method="reg")
+        p = est.get_params()
+        assert p["covariates"] == ["x1"]
+        assert p["estimation_method"] == "reg"
+        assert "pscore_trim" in p and "epv_threshold" in p and "pscore_fallback" in p
+
+    def test_default_estimation_method_is_dr(self):
+        assert ContinuousDiD().estimation_method == "dr"
+
+    def test_invalid_estimation_method_raises(self):
+        with pytest.raises(ValueError, match="estimation_method"):
+            ContinuousDiD(estimation_method="bogus")
+
+    def test_set_params_transactional(self):
+        est = ContinuousDiD(estimation_method="reg")
+        with pytest.raises(ValueError):
+            est.set_params(estimation_method="nope")
+        # invalid update left the config unmutated
+        assert est.estimation_method == "reg"
+        est.set_params(estimation_method="dr")
+        assert est.estimation_method == "dr"
+
+    def test_ipw_with_covariates_raises(self):
+        data = _cov_data()
+        est = ContinuousDiD(covariates=["x1"], estimation_method="ipw")
+        with pytest.raises(NotImplementedError, match="ipw"):
+            est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+
+    def test_ipw_without_covariates_ok(self):
+        # estimation_method only matters with covariates; ipw default must not
+        # break the unconditional path.
+        data = _cov_data()
+        est = ContinuousDiD(estimation_method="ipw")
+        res = est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+        assert np.isfinite(res.overall_att)
+
+    def test_survey_with_covariates_raises(self):
+        from diff_diff.survey import SurveyDesign
+
+        data = _cov_data()
+        data["w"] = 1.0
+        est = ContinuousDiD(covariates=["x1"], estimation_method="reg")
+        with pytest.raises(NotImplementedError, match="survey_design"):
+            est.fit(
+                data,
+                "outcome",
+                "unit",
+                "period",
+                "first_treat",
+                "dose",
+                aggregate="dose",
+                survey_design=SurveyDesign(weights="w"),
+            )
+
+    def test_missing_covariate_column_raises(self):
+        data = _cov_data()
+        est = ContinuousDiD(covariates=["not_a_col"], estimation_method="reg")
+        with pytest.raises(ValueError, match="not found"):
+            est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+
+    def test_missing_covariate_values_raise(self):
+        # Fail closed: a per-cell fallback would silently mix conditional and
+        # unconditional estimands in the aggregate.
+        data = _cov_data()
+        data.loc[data.index[:4], "x1"] = np.nan
+        est = ContinuousDiD(covariates=["x1"], estimation_method="reg")
+        with pytest.raises(ValueError, match="missing/non-finite covariate"):
+            est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+
+    def test_default_pscore_fallback_is_error(self):
+        assert ContinuousDiD().pscore_fallback == "error"
+
+    def test_invalid_nuisance_params_raise(self):
+        with pytest.raises(ValueError, match="pscore_fallback"):
+            ContinuousDiD(pscore_fallback="maybe")
+        with pytest.raises(ValueError, match="pscore_trim"):
+            ContinuousDiD(pscore_trim=0.6)
+        with pytest.raises(ValueError, match="pscore_trim"):
+            ContinuousDiD(pscore_trim=-0.1)
+        with pytest.raises(ValueError, match="epv_threshold"):
+            ContinuousDiD(epv_threshold=0.0)
+
+    def test_covariate_metadata_on_results(self):
+        data = _cov_data()
+        est = ContinuousDiD(
+            covariates=["x1"],
+            estimation_method="dr",
+            pscore_trim=0.02,
+            epv_threshold=8.0,
+            pscore_fallback="unconditional",
+        )
+        res = est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+        assert res.covariates == ["x1"]
+        assert res.estimation_method == "dr"
+        assert res.pscore_trim == 0.02
+        assert res.epv_threshold == 8.0
+        assert res.pscore_fallback == "unconditional"
+        assert "Covariates" in res.summary()
+
+    def test_covariate_eventstudy_and_bootstrap(self):
+        """Covariate paths compose with aggregate='eventstudy' and n_bootstrap>0
+        (reg + dr), with finite inference and bootstrap SE near analytical."""
+        data = _cov_data(n_units=200)
+        for method in ("reg", "dr"):
+            es = ContinuousDiD(covariates=["x1"], estimation_method=method).fit(
+                data, "outcome", "unit", "period", "first_treat", "dose", aggregate="eventstudy"
+            )
+            assert np.isfinite(es.overall_att) and np.isfinite(es.overall_att_se)
+            ana = ContinuousDiD(covariates=["x1"], estimation_method=method).fit(
+                data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose"
+            )
+            boot = ContinuousDiD(
+                covariates=["x1"], estimation_method=method, n_bootstrap=199, seed=3
+            ).fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+            assert np.isfinite(boot.overall_att_se)
+            # bootstrap SE within ~30% of analytical (same linearized IF)
+            assert abs(boot.overall_att_se - ana.overall_att_se) / ana.overall_att_se < 0.3
+
+    def test_clone_refit_idempotent(self):
+        data = _cov_data()
+        est = ContinuousDiD(covariates=["x1"], estimation_method="dr", seed=1)
+        r1 = est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+        clone = ContinuousDiD(**est.get_params())
+        r2 = clone.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+        assert abs(float(r1.overall_att) - float(r2.overall_att)) < 1e-12
+        assert abs(float(r1.overall_att_se) - float(r2.overall_att_se)) < 1e-12
+
+    def test_no_covariate_path_unchanged(self):
+        """Passing covariates=None runs the unchanged unconditional path."""
+        data = _cov_data()
+        base = ContinuousDiD(seed=1).fit(
+            data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose"
+        )
+        # estimation_method has no effect without covariates
+        alt = ContinuousDiD(estimation_method="reg", seed=1).fit(
+            data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose"
+        )
+        assert float(base.overall_att) == float(alt.overall_att)
+        assert float(base.overall_att_se) == float(alt.overall_att_se)
diff --git a/tests/test_methodology_continuous_did.py b/tests/test_methodology_continuous_did.py
index 89e60c745..b50200181 100644
--- a/tests/test_methodology_continuous_did.py
+++ b/tests/test_methodology_continuous_did.py
@@ -48,18 +48,14 @@ def linear_data(self):
     def test_linear_att_recovery(self, linear_data):
         """With degree=1 and linear truth, ATT(d) should be exactly 2d."""
         est = ContinuousDiD(degree=1, num_knots=0, dvals=np.array([1.0, 3.0, 5.0]))
-        results = est.fit(
-            linear_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(linear_data, "outcome", "unit", "period", "first_treat", "dose")
         expected = np.array([2.0, 6.0, 10.0])
         np.testing.assert_allclose(results.dose_response_att.effects, expected, atol=1e-10)
 
     def test_linear_acrt(self, linear_data):
         """ACRT(d) should be constant = 2 for linear truth."""
         est = ContinuousDiD(degree=1, num_knots=0, dvals=np.array([1.5, 3.0, 4.5]))
-        results = est.fit(
-            linear_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(linear_data, "outcome", "unit", "period", "first_treat", "dose")
         # Derivative of 2d is 2
         np.testing.assert_allclose(results.dose_response_acrt.effects, 2.0, atol=1e-6)
 
@@ -71,17 +67,13 @@ def test_att_glob_binarized(self, linear_data):
         expected_att_glob = mean_delta_treated - mean_delta_control
 
         est = ContinuousDiD(degree=1, num_knots=0)
-        results = est.fit(
-            linear_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(linear_data, "outcome", "unit", "period", "first_treat", "dose")
         np.testing.assert_allclose(results.overall_att, expected_att_glob, atol=1e-10)
 
     def test_acrt_glob_plugin(self, linear_data):
         """ACRT_glob = mean(ACRT(D_i)) over treated = 2."""
         est = ContinuousDiD(degree=1, num_knots=0)
-        results = est.fit(
-            linear_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(linear_data, "outcome", "unit", "period", "first_treat", "dose")
         np.testing.assert_allclose(results.overall_acrt, 2.0, atol=1e-6)
 
 
@@ -109,13 +101,9 @@ def test_quadratic_recovery(self, quadratic_data):
         """Cubic basis should recover d^2 exactly."""
         eval_grid = np.array([1.5, 2.5, 3.5, 4.5])
         est = ContinuousDiD(degree=3, num_knots=0, dvals=eval_grid)
-        results = est.fit(
-            quadratic_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(quadratic_data, "outcome", "unit", "period", "first_treat", "dose")
         expected = eval_grid**2
-        np.testing.assert_allclose(
-            results.dose_response_att.effects, expected, atol=1e-6
-        )
+        np.testing.assert_allclose(results.dose_response_att.effects, expected, atol=1e-6)
 
 
 class TestMultiPeriodAggregation:
@@ -136,18 +124,14 @@ def staggered_data(self):
 
     def test_multiple_groups(self, staggered_data):
         est = ContinuousDiD(degree=1, num_knots=0)
-        results = est.fit(
-            staggered_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(staggered_data, "outcome", "unit", "period", "first_treat", "dose")
         assert len(results.groups) == 2
         assert 2 in results.groups
         assert 3 in results.groups
 
     def test_gt_cell_count(self, staggered_data):
         est = ContinuousDiD(degree=1, num_knots=0)
-        results = est.fit(
-            staggered_data, "outcome", "unit", "period", "first_treat", "dose"
-        )
+        results = est.fit(staggered_data, "outcome", "unit", "period", "first_treat", "dose")
         # Group 2: periods 1(pre-via-varying),2,3,4; Group 3: periods 2(pre),3,4
         # Exact count depends on base period logic
         assert len(results.group_time_effects) >= 4
@@ -167,8 +151,12 @@ def test_all_same_dose(self):
             rows.append({"unit": i, "period": 2, "outcome": 0.0, "first_treat": 0, "dose": 0.0})
         for j in range(n_treated):
             uid = n_control + j
-            rows.append({"unit": uid, "period": 1, "outcome": 0.0, "first_treat": 2, "dose": dose_val})
-            rows.append({"unit": uid, "period": 2, "outcome": 5.0, "first_treat": 2, "dose": dose_val})
+            rows.append(
+                {"unit": uid, "period": 1, "outcome": 0.0, "first_treat": 2, "dose": dose_val}
+            )
+            rows.append(
+                {"unit": uid, "period": 2, "outcome": 5.0, "first_treat": 2, "dose": dose_val}
+            )
 
         data = pd.DataFrame(rows)
         est = ContinuousDiD(degree=1, num_knots=0, rank_deficient_action="silent")
@@ -202,17 +190,34 @@ def test_all_same_dose(self):
         for i in range(n_control):
             y_pre = rng.normal(0, 0.3)
             y_post = rng.normal(0, 0.3)
-            rows_hetero.append({"unit": i, "period": 1, "outcome": y_pre, "first_treat": 0, "dose": 0.0})
-            rows_hetero.append({"unit": i, "period": 2, "outcome": y_post, "first_treat": 0, "dose": 0.0})
+            rows_hetero.append(
+                {"unit": i, "period": 1, "outcome": y_pre, "first_treat": 0, "dose": 0.0}
+            )
+            rows_hetero.append(
+                {"unit": i, "period": 2, "outcome": y_post, "first_treat": 0, "dose": 0.0}
+            )
         for j in range(n_treated):
             uid = n_control + j
             y_pre = rng.normal(0, 0.3)
-            rows_hetero.append({"unit": uid, "period": 1, "outcome": y_pre, "first_treat": 2, "dose": dose_val})
-            rows_hetero.append({"unit": uid, "period": 2, "outcome": y_pre + 5.0, "first_treat": 2, "dose": dose_val})
+            rows_hetero.append(
+                {"unit": uid, "period": 1, "outcome": y_pre, "first_treat": 2, "dose": dose_val}
+            )
+            rows_hetero.append(
+                {
+                    "unit": uid,
+                    "period": 2,
+                    "outcome": y_pre + 5.0,
+                    "first_treat": 2,
+                    "dose": dose_val,
+                }
+            )
         data_hetero = pd.DataFrame(rows_hetero)
         est_boot = ContinuousDiD(
-            degree=1, num_knots=0, n_bootstrap=199,
-            rank_deficient_action="silent", seed=42,
+            degree=1,
+            num_knots=0,
+            n_bootstrap=199,
+            rank_deficient_action="silent",
+            seed=42,
         )
         with pytest.warns(UserWarning, match="[Ii]dentical"):
             results_boot = est_boot.fit(
@@ -245,7 +250,9 @@ def _check_r_contdid():
     try:
         result = subprocess.run(
             ["Rscript", "-e", "library(contdid); cat('OK')"],
-            capture_output=True, text=True, timeout=10,
+            capture_output=True,
+            text=True,
+            timeout=10,
         )
         return result.stdout.strip() == "OK"
     except (FileNotFoundError, subprocess.TimeoutExpired):
@@ -260,8 +267,14 @@ def _check_r_contdid():
 )
 
 
-def _run_r_contdid(csv_path, degree=3, num_knots=0, control_group="nevertreated",
-                    aggregation="dose", staggered=False):
+def _run_r_contdid(
+    csv_path,
+    degree=3,
+    num_knots=0,
+    control_group="nevertreated",
+    aggregation="dose",
+    staggered=False,
+):
     """Run R's cont_did() and return results for comparison.
 
     For 2-period data (staggered=False): recomputes ATT(d)/ACRT(d) with consistent
@@ -381,7 +394,9 @@ def _run_r_contdid(csv_path, degree=3, num_knots=0, control_group="nevertreated"
         """
     result = subprocess.run(
         ["Rscript", "-e", r_code],
-        capture_output=True, text=True, timeout=120,
+        capture_output=True,
+        text=True,
+        timeout=120,
     )
     if result.returncode != 0:
         pytest.skip(f"R contdid failed: {result.stderr[:500]}")
@@ -392,17 +407,28 @@ def _run_r_contdid(csv_path, degree=3, num_knots=0, control_group="nevertreated"
 class TestRBenchmark:
     """R `contdid` v0.1.0 benchmark comparisons."""
 
-    def _compare_with_r(self, data, degree=3, num_knots=0,
-                        control_group="never_treated", aggregation="dose",
-                        staggered=False, att_tol=0.01, acrt_tol=0.02):
+    def _compare_with_r(
+        self,
+        data,
+        degree=3,
+        num_knots=0,
+        control_group="never_treated",
+        aggregation="dose",
+        staggered=False,
+        att_tol=0.01,
+        acrt_tol=0.02,
+    ):
         """Helper: run both Python and R, compare."""
         with tempfile.NamedTemporaryFile(suffix=".csv", mode="w", delete=False) as f:
             data.to_csv(f, index=False)
             csv_path = f.name
 
         r_out = _run_r_contdid(
-            csv_path, degree=degree, num_knots=num_knots,
-            control_group=control_group, aggregation=aggregation,
+            csv_path,
+            degree=degree,
+            num_knots=num_knots,
+            control_group=control_group,
+            aggregation=aggregation,
             staggered=staggered,
         )
 
@@ -416,11 +442,18 @@ def _compare_with_r(self, data, degree=3, num_knots=0,
         # Python estimation using R's dvals for exact grid match
         dvals = np.array(r_out["dvals"])
         est = ContinuousDiD(
-            degree=degree, num_knots=num_knots, dvals=dvals,
+            degree=degree,
+            num_knots=num_knots,
+            dvals=dvals,
             control_group=control_group,
         )
         results = est.fit(
-            data, "outcome", "unit", "period", "first_treat", "dose",
+            data,
+            "outcome",
+            "unit",
+            "period",
+            "first_treat",
+            "dose",
             aggregate=py_aggregate,
         )
 
@@ -461,24 +494,33 @@ def _compare_with_r(self, data, degree=3, num_knots=0,
     def test_benchmark_1_basic_cubic(self):
         """2 periods, 1 cohort, degree=3, no knots, never_treated."""
         data = generate_continuous_did_data(
-            n_units=300, n_periods=2, cohort_periods=[2],
-            seed=100, noise_sd=0.5,
+            n_units=300,
+            n_periods=2,
+            cohort_periods=[2],
+            seed=100,
+            noise_sd=0.5,
         )
         self._compare_with_r(data, degree=3, num_knots=0)
 
     def test_benchmark_2_linear(self):
         """2 periods, 1 cohort, degree=1 (linear), never_treated."""
         data = generate_continuous_did_data(
-            n_units=300, n_periods=2, cohort_periods=[2],
-            seed=101, noise_sd=0.5,
+            n_units=300,
+            n_periods=2,
+            cohort_periods=[2],
+            seed=101,
+            noise_sd=0.5,
         )
         self._compare_with_r(data, degree=1, num_knots=0)
 
     def test_benchmark_3_interior_knots(self):
         """2 periods, 1 cohort, degree=3, 2 interior knots."""
         data = generate_continuous_did_data(
-            n_units=300, n_periods=2, cohort_periods=[2],
-            seed=102, noise_sd=0.5,
+            n_units=300,
+            n_periods=2,
+            cohort_periods=[2],
+            seed=102,
+            noise_sd=0.5,
         )
         self._compare_with_r(data, degree=3, num_knots=2)
 
@@ -532,36 +574,54 @@ def test_benchmark_4_staggered_dose(self):
             """
             result = subprocess.run(
                 ["Rscript", "-e", r_code],
-                capture_output=True, text=True, timeout=120,
+                capture_output=True,
+                text=True,
+                timeout=120,
             )
             if result.returncode != 0:
                 pytest.skip(f"R contdid failed: {result.stderr[:500]}")
             r_out = json.loads(result.stdout)
 
             data = pd.read_csv(tmp_path)
-            data = data.rename(columns={
-                "id": "unit", "time_period": "period",
-                "Y": "outcome", "G": "first_treat", "D": "dose",
-            })
+            data = data.rename(
+                columns={
+                    "id": "unit",
+                    "time_period": "period",
+                    "Y": "outcome",
+                    "G": "first_treat",
+                    "D": "dose",
+                }
+            )
             dvals = np.array(r_out["dvals"])
             est = ContinuousDiD(
-                degree=3, num_knots=0, dvals=dvals,
+                degree=3,
+                num_knots=0,
+                dvals=dvals,
                 control_group="never_treated",
             )
             results = est.fit(
-                data, "outcome", "unit", "period", "first_treat", "dose",
+                data,
+                "outcome",
+                "unit",
+                "period",
+                "first_treat",
+                "dose",
                 aggregate="dose",
             )
 
             # Overall ATT
-            att_diff = abs(results.overall_att - r_out["overall_att"]) / (abs(r_out["overall_att"]) + 1e-10)
+            att_diff = abs(results.overall_att - r_out["overall_att"]) / (
+                abs(r_out["overall_att"]) + 1e-10
+            )
             assert att_diff < 0.01, (
                 f"Overall ATT diff: {att_diff:.4f} "
                 f"(R={r_out['overall_att']:.6f}, Py={results.overall_att:.6f})"
             )
 
             # Overall ACRT
-            acrt_diff = abs(results.overall_acrt - r_out["overall_acrt"]) / (abs(r_out["overall_acrt"]) + 1e-10)
+            acrt_diff = abs(results.overall_acrt - r_out["overall_acrt"]) / (
+                abs(r_out["overall_acrt"]) + 1e-10
+            )
             assert acrt_diff < 0.01, (
                 f"Overall ACRT diff: {acrt_diff:.4f} "
                 f"(R={r_out['overall_acrt']:.6f}, Py={results.overall_acrt:.6f})"
@@ -612,34 +672,52 @@ def test_benchmark_5_not_yet_treated(self):
             """
             result = subprocess.run(
                 ["Rscript", "-e", r_code],
-                capture_output=True, text=True, timeout=120,
+                capture_output=True,
+                text=True,
+                timeout=120,
             )
             if result.returncode != 0:
                 pytest.skip(f"R contdid failed: {result.stderr[:500]}")
             r_out = json.loads(result.stdout)
 
             data = pd.read_csv(tmp_path)
-            data = data.rename(columns={
-                "id": "unit", "time_period": "period",
-                "Y": "outcome", "G": "first_treat", "D": "dose",
-            })
+            data = data.rename(
+                columns={
+                    "id": "unit",
+                    "time_period": "period",
+                    "Y": "outcome",
+                    "G": "first_treat",
+                    "D": "dose",
+                }
+            )
             dvals = np.array(r_out["dvals"])
             est = ContinuousDiD(
-                degree=3, num_knots=0, dvals=dvals,
+                degree=3,
+                num_knots=0,
+                dvals=dvals,
                 control_group="not_yet_treated",
             )
             results = est.fit(
-                data, "outcome", "unit", "period", "first_treat", "dose",
+                data,
+                "outcome",
+                "unit",
+                "period",
+                "first_treat",
+                "dose",
                 aggregate="dose",
             )
 
-            att_diff = abs(results.overall_att - r_out["overall_att"]) / (abs(r_out["overall_att"]) + 1e-10)
+            att_diff = abs(results.overall_att - r_out["overall_att"]) / (
+                abs(r_out["overall_att"]) + 1e-10
+            )
             assert att_diff < 0.01, (
                 f"Overall ATT diff: {att_diff:.4f} "
                 f"(R={r_out['overall_att']:.6f}, Py={results.overall_att:.6f})"
             )
 
-            acrt_diff = abs(results.overall_acrt - r_out["overall_acrt"]) / (abs(r_out["overall_acrt"]) + 1e-10)
+            acrt_diff = abs(results.overall_acrt - r_out["overall_acrt"]) / (
+                abs(r_out["overall_acrt"]) + 1e-10
+            )
             assert acrt_diff < 0.01, (
                 f"Overall ACRT diff: {acrt_diff:.4f} "
                 f"(R={r_out['overall_acrt']:.6f}, Py={results.overall_acrt:.6f})"
@@ -686,31 +764,363 @@ def test_benchmark_6_event_study(self):
             """
             result = subprocess.run(
                 ["Rscript", "-e", r_code],
-                capture_output=True, text=True, timeout=120,
+                capture_output=True,
+                text=True,
+                timeout=120,
             )
             if result.returncode != 0:
                 pytest.skip(f"R contdid failed: {result.stderr[:500]}")
             r_out = json.loads(result.stdout)
 
             data = pd.read_csv(tmp_path)
-            data = data.rename(columns={
-                "id": "unit", "time_period": "period",
-                "Y": "outcome", "G": "first_treat", "D": "dose",
-            })
+            data = data.rename(
+                columns={
+                    "id": "unit",
+                    "time_period": "period",
+                    "Y": "outcome",
+                    "G": "first_treat",
+                    "D": "dose",
+                }
+            )
             est = ContinuousDiD(
-                degree=3, num_knots=0,
+                degree=3,
+                num_knots=0,
                 control_group="never_treated",
             )
             results = est.fit(
-                data, "outcome", "unit", "period", "first_treat", "dose",
+                data,
+                "outcome",
+                "unit",
+                "period",
+                "first_treat",
+                "dose",
                 aggregate="eventstudy",
             )
 
             # Compare overall ATT (binarized)
-            att_diff = abs(results.overall_att - r_out["overall_att"]) / (abs(r_out["overall_att"]) + 1e-10)
+            att_diff = abs(results.overall_att - r_out["overall_att"]) / (
+                abs(r_out["overall_att"]) + 1e-10
+            )
             assert att_diff < 0.01, (
                 f"Overall ATT diff: {att_diff:.4f} "
                 f"(R={r_out['overall_att']:.6f}, Py={results.overall_att:.6f})"
             )
         finally:
             os.unlink(tmp_path)
+
+
+# =============================================================================
+# Phase 3: Covariate adjustment (conditional parallel trends) — reg + dr
+#
+# The covariate-adjusted dose curve has NO external anchor (`contdid` v0.1.0
+# hard-stops on covariates). Validation strategy:
+#   - scalar overall_att + SE map EXACTLY onto DRDID reg_did_panel / drdid_panel
+#     (skip-guarded, DRDID not in CI);
+#   - an R-free NumPy reconstruction of the reg/dr att + SE runs IN CI (the guard
+#     the p=1 reduction cannot provide for the dr propensity/augmentation terms);
+#   - reg/dr ACRT(d) identity, DGP recovery, and MC coverage (R-free).
+# =============================================================================
+
+
+def _check_r_drdid():
+    try:
+        result = subprocess.run(
+            ["Rscript", "-e", "library(DRDID); cat('OK')"],
+            capture_output=True,
+            text=True,
+            timeout=10,
+        )
+        return result.stdout.strip() == "OK"
+    except (FileNotFoundError, subprocess.TimeoutExpired):
+        return False
+
+
+_HAS_R_DRDID = _check_r_drdid()
+require_drdid = pytest.mark.skipif(not _HAS_R_DRDID, reason="R or DRDID package not installed")
+
+
+def _covariate_cell_data(seed=101, n=500):
+    """2-period single-cohort panel with 2 covariates and conditional PT.
+
+    Returns (df, dY, D, X_with_intercept) so the estimator's single (g=2,t=2)
+    cell overall_att maps onto one DRDID reg/dr att.
+    """
+    rng = np.random.default_rng(seed)
+    x1 = rng.normal(size=n)
+    x2 = rng.normal(size=n)
+    treated = (rng.uniform(size=n) < 1 / (1 + np.exp(-(0.4 * x1 - 0.3 * x2)))).astype(int)
+    dose = np.where(treated == 1, rng.uniform(0.2, 2.0, n), 0.0)
+    g = np.where(treated == 1, 2, 0)
+    y1 = 1.0 + 0.5 * x1 + rng.normal(0, 0.3, n)
+    y2 = (
+        y1 + (0.8 * x1 - 0.5 * x2) + np.where(treated == 1, 1.5 * dose, 0.0) + rng.normal(0, 0.3, n)
+    )
+    rows = []
+    for i in range(n):
+        for period, y in [(1, y1[i]), (2, y2[i])]:
+            rows.append(
+                {
+                    "unit": i,
+                    "period": period,
+                    "outcome": y,
+                    "first_treat": g[i],
+                    "dose": dose[i],
+                    "x1": x1[i],
+                    "x2": x2[i],
+                }
+            )
+    df = pd.DataFrame(rows)
+    dY = y2 - y1
+    D = treated.astype(float)
+    X = np.column_stack([np.ones(n), x1, x2])
+    return df, dY, D, X
+
+
+def _numpy_reg_dr(dY, D, X):
+    """Independent NumPy reconstruction of DRDID reg_did_panel / drdid_panel
+    att + SE (unit weights=1, moderate overlap => no trimming). Separate
+    implementation from the estimator internals — the in-CI transcription guard.
+    """
+    n = len(D)
+    cf = D == 0
+    gamma = np.linalg.lstsq(X[cf], dY[cf], rcond=None)[0]
+    out = X @ gamma
+
+    def se(inf):
+        return inf.std(ddof=1) * np.sqrt(n - 1) / n
+
+    # reg
+    eta_t = (D * dY).mean() / D.mean()
+    eta_c = (D * out).mean() / D.mean()
+    reg_att = eta_t - eta_c
+    XpX = ((1 - D)[:, None] * X).T @ X / n
+    XpXi = np.linalg.inv(XpX)
+    asy = ((1 - D) * (dY - out))[:, None] * X @ XpXi
+    inf_t = (D * dY - D * eta_t) / D.mean()
+    inf_c1 = D * out - D * eta_c
+    M1 = (D[:, None] * X).mean(0)
+    reg_inf = inf_t - (inf_c1 + asy @ M1) / D.mean()
+    # dr
+    ps = 1 / (1 + np.exp(-X @ np.linalg.lstsq(X, D, rcond=None)[0]))  # rough; refit below
+    # proper logit via IRLS
+    beta = np.zeros(X.shape[1])
+    for _ in range(100):
+        mu = 1 / (1 + np.exp(-X @ beta))
+        W = mu * (1 - mu)
+        z = X @ beta + (D - mu) / np.clip(W, 1e-10, None)
+        beta_new = np.linalg.solve((X * W[:, None]).T @ X, (X * W[:, None]).T @ z)
+        if np.max(np.abs(beta_new - beta)) < 1e-12:
+            beta = beta_new
+            break
+        beta = beta_new
+    ps = np.clip(1 / (1 + np.exp(-X @ beta)), 0.01, 0.99)
+    wt = D
+    wc = ps * (1 - D) / (1 - ps)
+    drt = wt * (dY - out)
+    drc = wc * (dY - out)
+    et = drt.mean() / wt.mean()
+    ec = drc.mean() / wc.mean()
+    dr_att = et - ec
+    asy_w = ((1 - D) * (dY - out))[:, None] * X @ XpXi
+    Wm = ps * (1 - ps)
+    Hess = np.linalg.inv(X.T @ (Wm[:, None] * X)) * n
+    asy_ps = ((D - ps)[:, None] * X) @ Hess
+    it = (drt - wt * et - asy_w @ (wt[:, None] * X).mean(0)) / wt.mean()
+    M2 = (wc[:, None] * (dY - out - ec)[:, None] * X).mean(0)
+    M3 = (wc[:, None] * X).mean(0)
+    ic = (drc - wc * ec + asy_ps @ M2 - asy_w @ M3) / wc.mean()
+    dr_inf = it - ic
+    return {"reg_att": reg_att, "reg_se": se(reg_inf), "dr_att": dr_att, "dr_se": se(dr_inf)}
+
+
+def _fit_cov(df, method, **kw):
+    est = ContinuousDiD(covariates=["x1", "x2"], estimation_method=method, **kw)
+    import warnings
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        return est.fit(df, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose")
+
+
+class TestCovariateReg:
+    """reg + dr covariate adjustment (conditional parallel trends)."""
+
+    @require_drdid
+    def test_reg_matches_drdid(self):
+        df, dY, D, X = _covariate_cell_data()
+        with tempfile.TemporaryDirectory() as tmp:
+            np.savetxt(f"{tmp}/dY.csv", dY, delimiter=",")
+            np.savetxt(f"{tmp}/D.csv", D, delimiter=",")
+            np.savetxt(f"{tmp}/X.csv", X, delimiter=",")
+            r = subprocess.run(
+                [
+                    "Rscript",
+                    "-e",
+                    f"""
+                suppressMessages(library(DRDID)); library(jsonlite)
+                dY<-as.numeric(read.csv("{tmp}/dY.csv",header=F)[,1])
+                D<-as.numeric(read.csv("{tmp}/D.csv",header=F)[,1])
+                X<-as.matrix(read.csv("{tmp}/X.csv",header=F)); n<-length(D)
+                o<-reg_did_panel(dY,rep(0,n),D,covariates=X)
+                cat(toJSON(list(att=o$ATT,se=o$se),auto_unbox=T,digits=12))
+            """,
+                ],
+                capture_output=True,
+                text=True,
+            )
+            if r.returncode != 0:
+                pytest.skip(f"R DRDID failed: {r.stderr[:300]}")
+            ref = json.loads(r.stdout)
+        res = _fit_cov(df, "reg")
+        assert abs(float(res.overall_att) - ref["att"]) < 1e-8
+        assert abs(float(res.overall_att_se) - ref["se"]) < 1e-8
+
+    @require_drdid
+    def test_dr_matches_drdid(self):
+        df, dY, D, X = _covariate_cell_data()
+        with tempfile.TemporaryDirectory() as tmp:
+            np.savetxt(f"{tmp}/dY.csv", dY, delimiter=",")
+            np.savetxt(f"{tmp}/D.csv", D, delimiter=",")
+            np.savetxt(f"{tmp}/X.csv", X, delimiter=",")
+            r = subprocess.run(
+                [
+                    "Rscript",
+                    "-e",
+                    f"""
+                suppressMessages(library(DRDID)); library(jsonlite)
+                dY<-as.numeric(read.csv("{tmp}/dY.csv",header=F)[,1])
+                D<-as.numeric(read.csv("{tmp}/D.csv",header=F)[,1])
+                X<-as.matrix(read.csv("{tmp}/X.csv",header=F)); n<-length(D)
+                o<-drdid_panel(dY,rep(0,n),D,covariates=X)
+                cat(toJSON(list(att=o$ATT,se=o$se),auto_unbox=T,digits=12))
+            """,
+                ],
+                capture_output=True,
+                text=True,
+            )
+            if r.returncode != 0:
+                pytest.skip(f"R DRDID failed: {r.stderr[:300]}")
+            ref = json.loads(r.stdout)
+        res = _fit_cov(df, "dr")
+        # dr att/se match DRDID (~1e-6; Python IRLS vs fastglm method=3)
+        assert abs(float(res.overall_att) - ref["att"]) < 1e-5
+        assert abs(float(res.overall_att_se) - ref["se"]) < 1e-5
+
+    def test_dr_reg_numpy_crosscheck_p2(self):
+        """R-free: estimator reg/dr att+se match an independent NumPy
+        reconstruction at p>=2. This is the CI guard the p=1 reduction cannot
+        provide for the dr propensity/augmentation terms (at p=1 dr collapses
+        to reg)."""
+        df, dY, D, X = _covariate_cell_data(seed=202)
+        ref = _numpy_reg_dr(dY, D, X)
+        reg = _fit_cov(df, "reg")
+        dr = _fit_cov(df, "dr")
+        assert abs(float(reg.overall_att) - ref["reg_att"]) < 1e-9
+        assert abs(float(reg.overall_att_se) - ref["reg_se"]) < 1e-9
+        assert abs(float(dr.overall_att) - ref["dr_att"]) < 1e-7
+        assert abs(float(dr.overall_att_se) - ref["dr_se"]) < 1e-7
+
+    def test_reg_vs_dr_acrt_identical(self):
+        """reg and dr share the dose-response shape: ACRT(d) point AND SE are
+        identical (the dr augmentation enters only the intercept direction,
+        which dPsi annihilates). ATT(d) differs by a single constant."""
+        df, _, _, _ = _covariate_cell_data(seed=303)
+        reg = _fit_cov(df, "reg")
+        dr = _fit_cov(df, "dr")
+        np.testing.assert_allclose(
+            reg.dose_response_acrt.effects, dr.dose_response_acrt.effects, atol=1e-10
+        )
+        np.testing.assert_allclose(reg.dose_response_acrt.se, dr.dose_response_acrt.se, atol=1e-10)
+        att_diff = reg.dose_response_att.effects - dr.dose_response_att.effects
+        assert np.ptp(att_diff) < 1e-9  # constant shift across the grid
+        assert abs(float(reg.overall_att) - float(dr.overall_att)) > 1e-4  # levels differ
+
+    def test_covariate_dgp_recovery(self):
+        """Conditional-PT DGP where unconditional DiD is biased but covariate
+        reg/dr recover the truth."""
+        tau = 2.0
+        rng = np.random.default_rng(7)
+        n = 600
+        x1 = rng.normal(size=n)
+        x2 = rng.normal(size=n)
+        tr = (rng.uniform(size=n) < 1 / (1 + np.exp(-(0.7 * x1 - 0.5 * x2)))).astype(int)
+        dose = np.where(tr == 1, rng.uniform(0.2, 2, n), 0.0)
+        g = np.where(tr == 1, 2, 0)
+        y1 = 1 + 0.5 * x1 + rng.normal(0, 0.3, n)
+        y2 = y1 + (1.2 * x1 - 0.8 * x2) + np.where(tr == 1, tau, 0.0) + rng.normal(0, 0.3, n)
+        rows = []
+        for i in range(n):
+            for period, y in [(1, y1[i]), (2, y2[i])]:
+                rows.append(
+                    {
+                        "unit": i,
+                        "period": period,
+                        "outcome": y,
+                        "first_treat": g[i],
+                        "dose": dose[i],
+                        "x1": x1[i],
+                        "x2": x2[i],
+                    }
+                )
+        df = pd.DataFrame(rows)
+        import warnings
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            uncond = ContinuousDiD().fit(
+                df, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose"
+            )
+        assert abs(float(uncond.overall_att) - tau) > 0.5  # unconditional biased
+        for method in ("reg", "dr"):
+            res = _fit_cov(df, method)
+            assert abs(float(res.overall_att) - tau) < 0.1  # recovered
+
+    @pytest.mark.slow
+    def test_covariate_coverage(self):
+        """Analytical SE achieves nominal coverage under conditional PT (reg &
+        dr). R-free validation of the curve/scalar SE (esp. the dr SE). Slow —
+        a fixed rep count is needed for a meaningful coverage estimate; the fast
+        CI guard is test_dr_reg_numpy_crosscheck_p2."""
+        tau = 2.0
+        n = 400
+        reps = 150
+        import warnings
+
+        for method in ("reg", "dr"):
+            cover = 0
+            total = 0
+            for s in range(reps):
+                rng = np.random.default_rng(10_000 + s)
+                x1 = rng.normal(size=n)
+                x2 = rng.normal(size=n)
+                tr = (rng.uniform(size=n) < 1 / (1 + np.exp(-(0.7 * x1 - 0.5 * x2)))).astype(int)
+                dose = np.where(tr == 1, rng.uniform(0.2, 2, n), 0.0)
+                g = np.where(tr == 1, 2, 0)
+                y1 = 1 + 0.5 * x1 + rng.normal(0, 0.3, n)
+                y2 = (
+                    y1 + (1.2 * x1 - 0.8 * x2) + np.where(tr == 1, tau, 0.0) + rng.normal(0, 0.3, n)
+                )
+                rows = []
+                for i in range(n):
+                    for period, y in [(1, y1[i]), (2, y2[i])]:
+                        rows.append(
+                            {
+                                "unit": i,
+                                "period": period,
+                                "outcome": y,
+                                "first_treat": g[i],
+                                "dose": dose[i],
+                                "x1": x1[i],
+                                "x2": x2[i],
+                            }
+                        )
+                df = pd.DataFrame(rows)
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore")
+                    res = _fit_cov(df, method)
+                lo, hi = res.overall_att_conf_int
+                cover += int(lo <= tau <= hi)
+                total += 1
+            rate = cover / total
+            # Wide band: MC noise at 150 reps; catches a broken (too small/large) SE.
+            assert 0.88 <= rate <= 0.99, f"{method} coverage {rate:.3f} off nominal"