[

0.05,

0.1,

0.15000000000000002,

0.2,

0.25,

0.3,

0.35000000000000003,

0.4,

0.45,

0.5,

0.55,

0.6000000000000001,

0.6500000000000001,

0.7000000000000001,

0.7500000000000001,

0.8,

0.8500000000000001,

0.9000000000000001,

0.95,

]

"QDiD is deprecated and will be removed in 4.0; use "

"ChangesInChanges(method='qdid') instead - the same engine, so the numbers "

"are unchanged. Only the class spelling is deprecated: method='qdid' is a "

"fully supported comparison mode and emits no warning of its own. The "

"QDiDResults alias is deprecated with it."

[

0.01,

0.02,

0.03,

0.04,

0.05,

0.060000000000000005,

0.06999999999999999,

0.08,

0.09,

0.09999999999999999,

0.11,

0.12,

0.13,

0.14,

0.15000000000000002,

0.16,

0.17,

0.18000000000000002,

0.19,

0.2,

0.21000000000000002,

0.22,

0.23,

0.24000000000000002,

0.25,

0.26,

0.27,

0.28,

0.29000000000000004,

0.3,

0.31,

0.32,

0.33,

0.34,

0.35000000000000003,

0.36000000000000004,

0.37,

0.38,

0.39,

0.4,

0.41000000000000003,

0.42000000000000004,

0.43,

0.44,

0.45,

0.46,

0.47000000000000003,

0.48000000000000004,

0.49,

0.5,

0.51,

0.52,

0.53,

0.54,

0.55,

0.56,

0.5700000000000001,

0.5800000000000001,

0.59,

0.6,

0.61,

0.62,

0.63,

0.64,

0.65,

0.66,

0.67,

0.68,

0.6900000000000001,

0.7000000000000001,

0.7100000000000001,

0.72,

0.73,

0.74,

0.75,

0.76,

0.77,

0.78,

0.79,

0.8,

0.81,

0.8200000000000001,

0.8300000000000001,

0.8400000000000001,

0.85,

0.86,

0.87,

0.88,

0.89,

0.9,

0.91,

0.92,

0.93,

0.9400000000000001,

0.9500000000000001,

0.9600000000000001,

0.97,

0.98,

0.99,

]

"y00": "control pre-period (treatment=0, time=0)",

"y01": "control post-period (treatment=0, time=1)",

"y10": "treated pre-period (treatment=1, time=0)",

"y11": "treated post-period (treatment=1, time=1)",

y: np.ndarray, g: np.ndarray, t: np.ndarray, X: Optional[np.ndarray] = None

) -> Optional[Dict[str, np.ndarray]]:

"""Split outcomes into the four sorted (group, period) cells.

cells = {}

for key, (gv, tv) in {"y00": (0, 0), "y01": (0, 1), "y10": (1, 0), "y11": (1, 1)}.items():

mask = (g == gv) & (t == tv)

cell = y[mask]

if cell.size == 0:

return None

if X is None:

cells[key] = np.sort(cell)

else:

order = np.argsort(cell, kind="stable")

cells[key] = cell[order]

cells["x" + key[1:]] = X[mask][order]

y: np.ndarray, g: np.ndarray, t: np.ndarray, X: Optional[np.ndarray] = None

) -> Dict[str, np.ndarray]:

"""Like :func:`_build_cells` but raises on empty cells (Assumption 5.1(ii))."""

for key, (gv, tv) in {"y00": (0, 0), "y01": (0, 1), "y10": (1, 0), "y11": (1, 1)}.items():

if not np.any((g == gv) & (t == tv)):

raise ValueError(

f"Empty (group, period) cell: no observations in the {_CELL_LABELS[key]} cell. "

"All four 2x2 cells must be non-empty (Athey-Imbens Assumption 5.1(ii))."

)

cells = _build_cells(y, g, t, X)

assert cells is not None

"""Empirical CDF of ``sorted_sample`` evaluated at ``x`` (eq. 34, ``<=`` semantics).

n = sorted_sample.shape[0]

"""R ``quantile(x, probs, type=1)`` - exact port including R's fuzz arithmetic.

n = sorted_sample.shape[0]

p = np.clip(np.asarray(probs, dtype=float), 0.0, 1.0)

fuzz = 4.0 * np.finfo(float).eps

nppm = n * p

j = np.floor(nppm * (1.0 + fuzz)).astype(np.int64)

h = (nppm > j).astype(np.int64)

idx = np.clip(j + h, 1, n)

"""Linear quantile regression via the Koenker-Bassett LP, one solve per tau.

n = y.shape[0]

X_design = np.column_stack([np.ones(n), X])

p = X_design.shape[1]

A_eq = np.hstack([X_design, np.eye(n), -np.eye(n)])

bounds = [(None, None)] * p + [(0.0, None)] * (2 * n)

coefs = np.empty((taus.shape[0], p))

for j, tau in enumerate(taus):

c = np.concatenate([np.zeros(p), np.full(n, tau), np.full(n, 1.0 - tau)])

res = linprog(c, A_eq=A_eq, b_eq=y, bounds=bounds, method="highs")

if res.status != 0 or res.x is None or not np.all(np.isfinite(res.x[:p])):

return None

coefs[j] = res.x[:p]

"""quantreg ``predict.rqs(type="Fhat", stepfun=TRUE)`` evaluated at ``y``.

n_obs = y.shape[0]

taus_ext = np.concatenate([taus[:1], taus])

ranks = np.empty(n_obs)

for i in range(n_obs):

o = np.argsort(preds[i], kind="stable")

yvals = np.concatenate([taus_ext[:1], taus_ext[o]])

ranks[i] = yvals[np.searchsorted(preds[i][o], y[i], side="right")]

"""quantreg ``predict.rqs(type="Qhat", stepfun=TRUE)`` evaluated at ``ranks``.

n_obs = ranks.shape[0]

out = np.empty(n_obs)

for i in range(n_obs):

yvals = np.concatenate([preds[i, :1], preds[i]])

out[i] = yvals[np.searchsorted(taus, ranks[i], side="right")]

cells: Dict[str, np.ndarray], quantiles: np.ndarray

) -> Tuple[float, np.ndarray, np.ndarray]:

"""CiC ATT (eq. 36) and quantile effects (eq. 18), qte::CiC() arithmetic.

ranks = _ecdf(cells["y00"], cells["y10"])

cf = _quantile_type1(cells["y01"], ranks)

att = float(np.mean(cells["y11"]) - np.mean(cf))

cf_sorted = np.sort(cf)

qte = _quantile_type1(cells["y11"], quantiles) - _quantile_type1(cf_sorted, quantiles)

"""QDiD ATT and quantile effects, matching qte::QDiD() exactly.

q1 = _quantile_type7(cells["y11"], quantiles)

q0 = (

_quantile_type7(cells["y10"], quantiles)

+ _quantile_type7(cells["y01"], quantiles)

- _quantile_type7(cells["y00"], quantiles)

)

ranks = _ecdf(cells["y10"], cells["y10"])

att = float(

np.mean(cells["y11"])

- (

np.mean(cells["y10"])

+ np.mean(_quantile_type7(cells["y01"], ranks))

- np.mean(_quantile_type7(cells["y00"], ranks))

)

)

cells: Dict[str, np.ndarray], quantiles: np.ndarray

) -> Optional[Tuple[float, np.ndarray, np.ndarray, np.ndarray]]:

"""CiC with covariates - exact port of qte::CiC()'s ``xformla`` branch.

coefs00 = _rq_fit(cells["y00"], cells["x00"], _QR_TAU_GRID)

coefs01 = _rq_fit(cells["y01"], cells["x01"], _QR_TAU_GRID)

if coefs00 is None or coefs01 is None:

return None

x10_design = _design_matrix(cells["x10"])

preds00 = x10_design @ coefs00.T

preds01 = x10_design @ coefs01.T

ranks = _fhat_eval(preds00, _QR_TAU_GRID, cells["y10"])

y0t = _qhat_eval(preds01, _QR_TAU_GRID, ranks)

att = float(np.mean(cells["y11"]) - np.mean(y0t))

y0t_sorted = np.sort(y0t)

qte = _quantile_type1(cells["y11"], quantiles) - _quantile_type1(y0t_sorted, quantiles)

envelope_flags = (cells["y10"] < preds00.min(axis=1)) | (cells["y10"] >= preds00.max(axis=1))

cells: Dict[str, np.ndarray], quantiles: np.ndarray

) -> Optional[Tuple[float, np.ndarray]]:

"""QDiD with covariates - exact port of qte::QDiD()'s ``xformla`` branch.

coefs00 = _rq_fit(cells["y00"], cells["x00"], _QR_TAU_GRID)

coefs01 = _rq_fit(cells["y01"], cells["x01"], _QR_TAU_GRID)

coefs10 = _rq_fit(cells["y10"], cells["x10"], _QR_TAU_GRID)

if coefs00 is None or coefs01 is None or coefs10 is None:

return None

x10_design = _design_matrix(cells["x10"])

preds10 = x10_design @ coefs10.T

preds01 = x10_design @ coefs01.T

preds00 = x10_design @ coefs00.T

ranks = _fhat_eval(preds10, _QR_TAU_GRID, cells["y10"])

y0t = (

cells["y10"]

+ _qhat_eval(preds01, _QR_TAU_GRID, ranks)

- _qhat_eval(preds00, _QR_TAU_GRID, ranks)

)

att = float(np.mean(cells["y11"]) - np.mean(y0t))

q1 = _quantile_type7(cells["y11"], quantiles)

q0 = _quantile_type1(np.sort(y0t), quantiles)

"""Eq. (17) plug-in interior range for CiC quantile effects.

y10 = cells["y10"]

q_lower = float(_ecdf(y10, np.array([cells["y00"][0]]))[0])

q_upper = float(_ecdf(y10, np.array([cells["y00"][-1]]))[0])

formula: str, data: pd.DataFrame

) -> Tuple[str, str, str, Optional[List[str]]]:

"""Parse ``"outcome ~ treatment * time [+ covariates]"`` style 2x2 formulas.

if "~" not in formula:

raise ValueError("Formula must contain '~' to separate outcome from predictors")

lhs, rhs = formula.split("~", 1)

outcome = lhs.strip()

rhs = rhs.strip()

covariates: Optional[List[str]] = None

if "*" in rhs:

parts = [p.strip() for p in rhs.split("*")]

if len(parts) != 2:

raise ValueError("Currently only supports single interaction (treatment * time)")

treatment, time = parts

if "+" in time:

time_parts = [p.strip() for p in time.split("+")]

time = time_parts[0]

covariates = time_parts[1:]

elif ":" in rhs:

terms = [t.strip() for t in rhs.split("+")]

interaction = None

mains: List[str] = []

for term in terms:

if ":" in term:

if interaction is not None:

raise ValueError("Formula must contain exactly one interaction term")

interaction = term

else:

mains.append(term)

if interaction is None:

raise ValueError(

"Formula must include an interaction term (treatment * time or treatment:time)"

)

pair = [p.strip() for p in interaction.split(":")]

if len(pair) != 2:

raise ValueError("Interaction term must involve exactly two variables")

if pair[0] == pair[1]:

raise ValueError("Interaction term must involve two distinct variables")

if pair[0] not in mains or pair[1] not in mains:

raise ValueError(

"Both variables in the interaction term must also appear as main effects "

"('outcome ~ treatment + time + treatment:time [+ covariates]'); roles are "

"taken from the main-effect order."

)

# Roles come from the MAIN-EFFECT order, not the interaction-term order;

# remaining main effects are covariates. First occurrences only, so a

# duplicated main effect cannot corrupt the role assignment.

role_mains: List[str] = []

for m in mains:

if m in pair and m not in role_mains:

role_mains.append(m)

treatment, time = role_mains[0], role_mains[1]

extras = [m for m in mains if m not in pair]

covariates = extras if extras else None

else:

raise ValueError(

"Formula must include an interaction term (treatment * time or treatment:time)"

)

for name in (outcome, treatment, time, *(covariates or [])):

if name not in data.columns:

raise ValueError(f"Column '{name}' from formula not found in data")

"""Warn on treated pre-period support outside the control pre-period range (CiC)."""

if cells["y10"][0] < cells["y00"][0] or cells["y10"][-1] > cells["y00"][-1]:

warnings.warn(

"Treated pre-period outcomes fall outside the control pre-period support "

"(Athey-Imbens Assumption 3.4 violated). The counterfactual distribution is "

"only partially identified (Corollary 3.1): quantile effects are reliable "

"only inside the reported (q_lower, q_upper) interior range, and the ATT "

"involves extrapolation at the support edges.",

UserWarning,

stacklevel=2,

"""Warn on conditional support failure under covariates (CiC).

share = float(np.mean(envelope_flags))

if share > _ENVELOPE_SHARE_WARN:

n_out = int(np.count_nonzero(envelope_flags))

warnings.warn(

f"{n_out} of {envelope_flags.size} treated pre-period outcomes ({share:.0%}) fall "

"outside their conditional quantile envelope (the span of the 99 predicted "

"control pre-period grid quantiles, taus 0.01-0.99, at their own covariates). "

"This suggests the conditional support/overlap condition "

"(Melly-Santangelo 2015, Assumption 4 - the covariate analogue of Athey-Imbens "

"Assumption 3.4) fails: conditional ranks for these observations are extrapolated "

"tail plateaus, and the counterfactual involves out-of-support extrapolation.",

UserWarning,

stacklevel=2,

"""Warn on heavy ties (discrete-looking outcomes) in any outcome cell."""

max_share = 0.0

for key in ("y00", "y01", "y10", "y11"):

cell = cells[key]

share = 1.0 - np.unique(cell).size / cell.size

max_share = max(max_share, share)

if max_share > _TIE_SHARE_WARN:

warnings.warn(

f"Outcome has heavy ties (up to {max_share:.0%} duplicate values within a "

"(group, period) cell), suggesting a discrete or mixed distribution. The "

"continuous CiC/QDiD machinery assumes continuously distributed outcomes "

"(Athey-Imbens Assumption 5.1(iii)); with discrete outcomes only bounds are "

"point-identified (Section 4, deferred) and the continuous formulas silently "

"deliver one endpoint of those bounds.",

UserWarning,

stacklevel=2,

"""Warn when QDiD's counterfactual quantile curve is non-monotone (footnote 21)."""

cq = (

_quantile_type7(cells["y10"], quantiles)

+ _quantile_type7(cells["y01"], quantiles)

- _quantile_type7(cells["y00"], quantiles)

)

if np.any(np.diff(cq) < -1e-12):

warnings.warn(

"QDiD's implied counterfactual quantile function is non-monotone on the "

"requested grid (Athey-Imbens footnote 21: the QDiD model places testable "

"restrictions on the data, and they appear violated here). Interpret the "

"quantile effects with caution; ChangesInChanges does not impose this "

"restriction and is the recommended estimator (p. 447).",

UserWarning,

stacklevel=2,

point_fn: Callable[[Dict[str, np.ndarray], np.ndarray], Optional[Tuple[Any, ...]]],

y: np.ndarray,

g: np.ndarray,

t: np.ndarray,

unit_ids: Optional[np.ndarray],

panel: bool,

n_bootstrap: int,

quantiles: np.ndarray,

rng: np.random.Generator,

X: Optional[np.ndarray] = None,

) -> np.ndarray:

"""Bootstrap replicate matrix, shape ``(n_bootstrap, 1 + K)`` (col 0 = ATT).

n_cols = 1 + quantiles.shape[0]

out = np.full((n_bootstrap, n_cols), np.nan)

if panel:

assert unit_ids is not None

# Pre-pivot to unit-level arrays: one (y_pre, y_post, group) triple per unit.

order = np.argsort(unit_ids, kind="stable")

uid, y_o, g_o, t_o = unit_ids[order], y[order], g[order], t[order]

X_o = X[order] if X is not None else None

pre_mask = t_o == 0

# Balanced panel (enforced in fit): each unit has exactly one pre and one post row.

y_pre = y_o[pre_mask]

y_post = y_o[~pre_mask]

g_unit = g_o[pre_mask]

X_pre = X_o[pre_mask] if X_o is not None else None

X_post = X_o[~pre_mask] if X_o is not None else None

# uid is sorted, so pre/post slices align unit-by-unit.

assert np.array_equal(uid[pre_mask], uid[~pre_mask])

n_units = y_pre.shape[0]

for b in range(n_bootstrap):

idx = rng.integers(0, n_units, n_units)

gb = g_unit[idx]

yb = np.concatenate([y_pre[idx], y_post[idx]])

tb = np.concatenate([np.zeros(n_units), np.ones(n_units)])

# Stacking order matches yb: pre-period rows first, then post.

Xb = (

np.vstack([X_pre[idx], X_post[idx]])

if X_pre is not None and X_post is not None

else None

)

cells = _build_cells(yb, np.concatenate([gb, gb]), tb, Xb)

if cells is None:

continue

res_b = point_fn(cells, quantiles)

if res_b is None:

continue

out[b, 0] = res_b[0]

out[b, 1:] = res_b[1]

else:

n_rows = y.shape[0]

for b in range(n_bootstrap):

idx = rng.integers(0, n_rows, n_rows)

Xb = X[idx] if X is not None else None

cells = _build_cells(y[idx], g[idx], t[idx], Xb)

if cells is None:

continue

res_b = point_fn(cells, quantiles)

if res_b is None:

continue

out[b, 0] = res_b[0]

out[b, 1:] = res_b[1]

replicates: np.ndarray,

qte_hat: np.ndarray,

n_bootstrap: int,

context: str,

) -> Tuple[float, np.ndarray, float, int]:

"""SEs and the sup-t critical value from the replicate matrix (qte conventions).

finite_rows = np.all(np.isfinite(replicates), axis=1)

n_valid = int(np.count_nonzero(finite_rows))

warn_bootstrap_failure_rate(n_valid, n_bootstrap, context)

if n_valid < max(2, _MIN_VALID_REPLICATE_SHARE * n_bootstrap):

k = replicates.shape[1] - 1

return np.nan, np.full(k, np.nan), np.nan, n_valid

good = replicates[finite_rows]

with np.errstate(divide="ignore", invalid="ignore"):

ses = np.std(good, axis=0, ddof=1)

att_se = float(ses[0])

qte_ses = ses[1:]

qte_cols = good[:, 1:]

z_iqr = stats.norm.ppf(0.75) - stats.norm.ppf(0.25)

q75 = np.array([_quantile_type1(np.sort(col), np.array([0.75]))[0] for col in qte_cols.T])

q25 = np.array([_quantile_type1(np.sort(col), np.array([0.25]))[0] for col in qte_cols.T])

sigmahalf = (q75 - q25) / z_iqr

if np.any(sigmahalf == 0):

sigmahalf = np.maximum(qte_ses, 1e-9)

sup_stats = np.max(np.abs(qte_cols - qte_hat[None, :]) / sigmahalf[None, :], axis=1)

sup_t_crit = float(_quantile_type1(np.sort(sup_stats), np.array([0.95]))[0])

est: Any,

data: pd.DataFrame,

outcome: Optional[str],

treatment: Optional[str],

time: Optional[str],

formula: Optional[str],

covariates: Optional[List[str]],

unit: Optional[str],

kind: str,

) -> ChangesInChangesResults:

"""Shared fit pipeline for ChangesInChanges and QDiD (``kind`` in {"cic", "qdid"})."""

# Re-validate hyperparameters (set_params may have mutated them since __init__).

_validate_all_params(est.get_params())

quantiles = np.sort(

np.asarray(_DEFAULT_QUANTILES if est.quantiles is None else est.quantiles, dtype=float)

)

# Derived from the INSTANCE, not from ``kind``: with method= on the merged class,

# a kind-based map would name "QDiD" in errors raised by a fit the user made

# through ChangesInChanges(method="qdid") - a class they never constructed. This

# is byte-identical to the old mapping on both 3.x surfaces (QDiD -> "QDiD",

# ChangesInChanges -> "ChangesInChanges").

estimator_name = type(est).__name__

# ---- column resolution -------------------------------------------------

if formula is not None:

# Uniform strictness (deliberately stricter than DifferenceInDifferences,

# which silently lets the formula win over explicit kwargs): mixing

# formula with explicit column arguments is ambiguous - reject it.

supplied = [

name

for name, value in (

("outcome", outcome),

("treatment", treatment),

("time", time),

("covariates", covariates),

)

if value is not None

]

if supplied:

raise ValueError(

"Provide either 'formula' or explicit column arguments, not both "

f"(got formula= together with {supplied})."

)

outcome, treatment, time, covariates = _parse_2x2_formula(formula, data)

elif outcome is None or treatment is None or time is None:

raise ValueError(

"Must provide either 'formula' or all of 'outcome', 'treatment', and 'time'"

)

# ---- covariate validation ------------------------------------------------

if isinstance(covariates, (str, bytes)):

# A bare string would iterate character-wise ("x1" -> ["x", "1"]) and

# could silently fit the wrong covariate set if such columns exist.

raise ValueError(

f"covariates must be a list of column names, got the bare string "

f"{covariates!r} - did you mean covariates=[{covariates!r}]?"

)

if covariates is not None and len(covariates) == 0:

covariates = None

if covariates is not None:

covariates = [str(c) for c in covariates]

reserved = {outcome, treatment, time}

if est.panel and unit is not None:

reserved.add(unit)

validate_covariate_names(covariates, reserved, estimator=estimator_name)

for col in covariates:

if col not in data.columns:

raise ValueError(f"Covariate column '{col}' not found in data")

if not pd.api.types.is_numeric_dtype(data[col]):

raise ValueError(

f"Covariate column '{col}' is not numeric. {estimator_name} accepts "

"numeric covariates only - dummy-encode categorical variables first "

"(e.g. pd.get_dummies(data, columns=[...]))."

)

used_cols = [outcome, treatment, time] + (covariates or [])

if est.panel:

if unit is None:

raise ValueError("'unit' is required when panel=True (unit identifier column)")

used_cols.append(unit)

for col in used_cols:

if col not in data.columns:

raise ValueError(f"Column '{col}' not found in data")

# ---- NA handling -------------------------------------------------------

frame = data[used_cols].copy()

n_before = len(frame)

frame = frame.dropna()

n_dropped = n_before - len(frame)

if n_dropped > 0:

warnings.warn(

f"Dropped {n_dropped} row(s) with missing values in " f"{used_cols} before estimation.",

UserWarning,

stacklevel=2,

)

if len(frame) == 0:

raise ValueError("No observations remain after dropping missing values")

y_check = frame[outcome].to_numpy(dtype=float)

if not np.all(np.isfinite(y_check)):

n_nonfinite = int(np.count_nonzero(~np.isfinite(y_check)))

raise ValueError(

f"Outcome column '{outcome}' contains {n_nonfinite} non-finite value(s) "

"(inf/-inf). Clean or drop these observations before fitting - they would "

"silently corrupt the empirical CDFs, quantiles, and bootstrap."

)

if covariates is not None:

x_check = frame[covariates].to_numpy(dtype=float)

if not np.all(np.isfinite(x_check)):

n_nonfinite = int(np.count_nonzero(~np.isfinite(x_check)))

raise ValueError(

f"Covariate column(s) {covariates} contain {n_nonfinite} non-finite "

"value(s) (inf/-inf). Clean or drop these observations before fitting - "

"they would silently corrupt the quantile regressions and bootstrap."

)

validate_binary(frame[treatment].to_numpy(dtype=float), "treatment")

validate_binary(frame[time].to_numpy(dtype=float), "time")

# ---- panel hygiene -----------------------------------------------------

unit_ids: Optional[np.ndarray] = None

if est.panel:

if frame.duplicated(subset=[unit, time]).any():

raise ValueError(

"panel=True requires at most one row per (unit, period); found duplicate "

f"('{unit}', '{time}') combinations."

)

g_nunique = frame.groupby(unit)[treatment].nunique()

if (g_nunique > 1).any():

bad = g_nunique[g_nunique > 1].index.tolist()[:5]

raise ValueError(

"The treatment-group indicator must be constant within unit in the 2x2 "

f"design (it marks group membership, not treatment receipt); units with "

f"varying values include {bad}."

)

counts = frame.groupby(unit)[time].count()

incomplete = counts[counts < 2].index

if len(incomplete) > 0:

warnings.warn(

f"Dropped {len(incomplete)} unit(s) not observed in both periods "

"(balanced-panel requirement, matching qte's makeBalancedPanel).",

UserWarning,

stacklevel=2,

)

frame = frame[~frame[unit].isin(incomplete)]

if len(frame) == 0:

raise ValueError("No balanced units remain after panel balancing")

unit_ids = frame[unit].to_numpy()

y = frame[outcome].to_numpy(dtype=float)

g = frame[treatment].to_numpy(dtype=float).astype(np.int64)

t = frame[time].to_numpy(dtype=float).astype(np.int64)

# Extracted AFTER panel balancing so dropped-unit rows leave X consistently.

X = frame[covariates].to_numpy(dtype=float) if covariates is not None else None

# ---- cells + diagnostics -----------------------------------------------

cells = _split_cells(y, g, t, X)

if X is not None:

# Fail closed at fit: below k+2 rows the quantile-regression LP has

# exact-fit degeneracy (rq produces garbage or errors there too).

# Bootstrap replicates instead NaN-row on LP failure, matching qte,

# which has no size guard.

k = X.shape[1]

qr_cells = ("y00", "y01") if kind == "cic" else ("y00", "y01", "y10")

for key in qr_cells:

if cells[key].size < k + 2:

raise ValueError(

f"Too few observations for quantile regression in the "