1"""Model F — prediction and economic value (analysis_plan §3 "Model F").23Targets on the S2 dynamic sample, all defined at a 30-day horizon from entry:4 B markdown30 first sustained markdown at listing age <= 30 days (listings that5 exit before any cut count as 0);6 B' markdown30cr competing-risk-aware variant: same event, but listings that exit7 within 30 days *without* a prior cut are removed from the sample8 (markdown among listings still exposed to the markdown risk);9 C exit30 observed disappearance confirmed by >= 7 days of absence at10 listing age <= 30 days (an exit is NOT a sale).11Target A (log entry price) is covered by the hedonic grouped CV (hedonic.grouped_cv)12and is not redone here.1314Observability. An exit at age a is confirmable only if entry_ts + a + 7 d <= snapshot;15the 30-day outcome of a listing is therefore determined only when16entry_ts <= SNAPSHOT_TS - 37 d (entry day <= 43, 2026-09-25 UTC). The analysis plan's17chronological cutoff (2026-10-01) leaves no such listing in the test fold, so the18primary chronological split uses 2026-09-15 UTC; the Oct-1 cutoff is reported as a19flagged sensitivity on the looser rule "entry day <= 81 - 30" (exit labels partly20right-censored there, see metrics.json notes).2122Features are measured AT ENTRY only. dealer_id is never a feature; dealer-level23outcome rates are target-encoded on the training fold only (inner 5-fold24out-of-fold encoding for training rows, smoothed towards the training prior).25"""26from __future__ import annotations2728import datetime as dt29import json30import math31import time32from pathlib import Path3334import numpy as np35import pandas as pd36import polars as pl37from sklearn.linear_model import LogisticRegression38from sklearn.metrics import brier_score_loss, log_loss, roc_auc_score3940from . import config as C4142OUT = C.RESULTS / "prediction"43HORIZON = 30.044SEED = 3645N_BOOT = 20046N_BINS = 1047DAY0 = math.floor(C.WINDOW_START_TS / 86400)48CHRONO_CUTOFF_TS = dt.datetime(2026, 9, 15, tzinfo=dt.timezone.utc).timestamp()49OCT1_TS = dt.datetime(2026, 10, 1, tzinfo=dt.timezone.utc).timestamp()50T_LAST = int(math.floor(C.SNAPSHOT_TS / 86400) - DAY0) # 8151STRICT_ENTRY_MAX_TS = C.SNAPSHOT_TS - (HORIZON + C.EXIT_CONFIRM_DAYS) * 8640052N_REGIONS_LOO = 653N_THREADS = 85455NUM_BASE = ["age_years_entry", "km10k", "year", "logp", "resid_A2", "week",56 "log_dens_10", "log_dens_25", "log_dens_50", "log_own_stock", "pricepos_25"]57CAT_COARSE = ["powertrain", "body", "region", "source_class"]58CAT_FINE = ["make", "model", "trim", "transmission", "drivetrain"]59CAT_ALL = CAT_COARSE + CAT_FINE60DEALER_FEATS = ["dealer_rate_markdown30", "dealer_rate_exit30", "dealer_log_n_train"]61TARGETS = {"markdown30": "y_markdown30", "markdown30cr": "y_markdown30cr", "exit30": "y_exit30"}62LGB_PARAMS = dict(n_estimators=2000, learning_rate=0.05, num_leaves=31, min_child_samples=50, subsample=0.8,63 subsample_freq=1, colsample_bytree=0.8, reg_lambda=1.0, cat_l2=10.0, cat_smooth=50.0,64 max_cat_to_onehot=4, n_jobs=N_THREADS, verbose=-1, random_state=SEED)656667def log(*a) -> None:68 print(time.strftime("%H:%M:%S"), *a, flush=True)697071# --------------------------------------------------------------------------------------72# pure helpers (unit-tested on synthetic arrays)73# --------------------------------------------------------------------------------------74def make_targets(entry_ts: np.ndarray, first_cut_ts: np.ndarray, duration_days: np.ndarray, exit_confirmed: np.ndarray,75 horizon: float = HORIZON) -> dict[str, np.ndarray]:76 """30-day outcome labels from listing-level timestamps.7778 first_cut_ts may be NaN (no sustained cut). Returns float arrays; markdown30cr is79 NaN for listings excluded by the competing-risk rule (exit before any cut within80 the horizon)."""81 entry_ts = np.asarray(entry_ts, float); fc = np.asarray(first_cut_ts, float)82 dur = np.asarray(duration_days, float); ex = np.asarray(exit_confirmed, bool)83 cut_age = (fc - entry_ts) / 86400.084 has_cut_h = np.isfinite(cut_age) & (cut_age <= horizon)85 exit_h = ex & (dur <= horizon)86 y_b = has_cut_h.astype(float)87 y_c = exit_h.astype(float)88 exit_before_cut = exit_h & ~has_cut_h89 y_bcr = np.where(exit_before_cut, np.nan, y_b)90 return {"y_markdown30": y_b, "y_markdown30cr": y_bcr, "y_exit30": y_c}919293def calibration_table(y: np.ndarray, p: np.ndarray, n_bins: int = N_BINS) -> pd.DataFrame:94 """Reliability table on equal-width probability bins [0,1]."""95 y = np.asarray(y, float); p = np.asarray(p, float)96 edges = np.linspace(0, 1, n_bins + 1)97 b = np.clip(np.digitize(p, edges[1:-1], right=False), 0, n_bins - 1)98 rows = []99 for k in range(n_bins):100 m = b == k101 rows.append({"bin": k, "lo": round(float(edges[k]), 6), "hi": round(float(edges[k + 1]), 6), "n": int(m.sum()),102 "mean_pred": float(p[m].mean()) if m.any() else np.nan,103 "frac_pos": float(y[m].mean()) if m.any() else np.nan})104 return pd.DataFrame(rows)105106107def expected_calibration_error(y: np.ndarray, p: np.ndarray, n_bins: int = N_BINS) -> float:108 t = calibration_table(y, p, n_bins)109 t = t[t["n"] > 0]110 return float((t["n"] / t["n"].sum() * (t["frac_pos"] - t["mean_pred"]).abs()).sum())111112113def _metrics(y, p) -> dict:114 p = np.clip(p, 1e-6, 1 - 1e-6)115 auc = float(roc_auc_score(y, p)) if 0 < y.mean() < 1 else float("nan")116 return {"log_loss": float(log_loss(y, p, labels=[0, 1])), "brier": float(brier_score_loss(y, p)), "auc": auc}117118119def bootstrap_metrics(y: np.ndarray, p: np.ndarray, groups: np.ndarray | None = None, n_boot: int = N_BOOT,120 seed: int = SEED) -> dict:121 """Point estimates + percentile bootstrap CIs (resampling test rows). When122 `groups` is given, rows are resampled within each group and the metric is the123 unweighted mean over groups (used for leave-one-region-out)."""124 y = np.asarray(y, float); p = np.asarray(p, float)125 rng = np.random.default_rng(seed)126 if groups is None:127 point = _metrics(y, p)128 idx_sets = [None]129 else:130 groups = np.asarray(groups)131 idx_sets = [np.flatnonzero(groups == g) for g in np.unique(groups)]132 point = {k: float(np.mean([_metrics(y[ix], p[ix])[k] for ix in idx_sets])) for k in ("log_loss", "brier", "auc")}133 draws = {k: [] for k in point}134 for _ in range(n_boot):135 if groups is None:136 ix = rng.integers(0, len(y), len(y)); m = _metrics(y[ix], p[ix])137 else:138 ms = []139 for ix in idx_sets:140 jx = ix[rng.integers(0, len(ix), len(ix))]; ms.append(_metrics(y[jx], p[jx]))141 m = {k: float(np.mean([x[k] for x in ms])) for k in point}142 for k in point:143 draws[k].append(m[k])144 out = {}145 for k in point:146 d = np.asarray(draws[k], float)147 out[k] = point[k]; out[f"{k}_ci"] = [float(np.nanpercentile(d, 2.5)), float(np.nanpercentile(d, 97.5))]148 out["ece"] = expected_calibration_error(y, p)149 out["n_test"] = int(len(y)); out["base_rate"] = float(y.mean()); out["n_boot"] = int(n_boot)150 return out151152153def target_encode(groups_tr: np.ndarray, y_tr: np.ndarray, groups_te: np.ndarray, m: float = 20.0, k: int = 5,154 seed: int = SEED) -> tuple[np.ndarray, np.ndarray]:155 """Smoothed group mean of y computed on the training fold only. Training rows get156 an inner k-fold out-of-fold encoding (their own label never enters their feature);157 test rows get the full-training encoding; unseen groups get the training prior."""158 g_tr = pd.Series(np.asarray(groups_tr)); y_tr = np.asarray(y_tr, float); g_te = pd.Series(np.asarray(groups_te))159 prior = float(np.nanmean(y_tr))160 ok = ~np.isnan(y_tr)161 def enc(gs: pd.Series, ys: np.ndarray, target: pd.Series) -> np.ndarray:162 d = pd.DataFrame({"g": gs.values, "y": ys}).dropna()163 agg = d.groupby("g")["y"].agg(["sum", "count"])164 s = target.map(agg["sum"]).fillna(0.0).values; n = target.map(agg["count"]).fillna(0.0).values165 return (s + m * prior) / (n + m)166 enc_tr = np.full(len(g_tr), prior)167 rng = np.random.default_rng(seed); fold = rng.integers(0, k, len(g_tr))168 for f in range(k):169 tr_in = (fold != f) & ok; te_in = fold == f170 enc_tr[te_in] = enc(g_tr[tr_in], y_tr[tr_in], g_tr[te_in])171 enc_te = enc(g_tr[ok], y_tr[ok], g_te)172 return enc_tr, enc_te173174175# --------------------------------------------------------------------------------------176# data177# --------------------------------------------------------------------------------------178def build_frame() -> pd.DataFrame:179 """S2 listings with entry-time features, risk-set competition on the entry day and180 the 30-day labels (all S2 entries; observability masks are applied per split)."""181 L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s2_dynamic"))182 L = L.select(["listing_id", "dealer_id", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km",183 "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "duration_days",184 "exit_confirmed", "first_cut_ts", "age_years_entry"])185 R = pl.read_parquet(C.PROCESSED / "risk_set.parquet",186 columns=["listing_id", "t", "listing_age", "dens_10", "dens_25", "dens_50", "own_stock", "pricepos_25"])187 first = R.filter(pl.col("listing_age") == 0).select(["listing_id", "t", "dens_10", "dens_25", "dens_50", "own_stock", "pricepos_25"])188 res = pl.read_parquet(C.RESULTS / "hedonic" / "residuals.parquet").select(["listing_id", "resid_A2"])189 df = L.join(first, on="listing_id", how="left").join(res, on="listing_id", how="left").to_pandas()190 df = df.rename(columns={"body_filled": "body"})191 for c, fill in [("body", "unknown"), ("trim", "(none)"), ("region", "(missing)"), ("transmission", "unknown"),192 ("drivetrain", "unknown"), ("powertrain", "unknown")]:193 df[c] = df[c].fillna(fill).astype(str)194 df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60)195 df["logp"] = np.log(df["entry_price"]); df["km10k"] = df["mileage_km"] / 1e4196 df["entry_day"] = (np.floor(df["entry_ts"] / 86400) - DAY0).astype(int)197 df["week"] = df["entry_day"] // 7198 for r in ("10", "25", "50"):199 df[f"log_dens_{r}"] = np.log1p(df[f"dens_{r}"].astype(float))200 df["log_own_stock"] = np.log1p(df["own_stock"].astype(float))201 df["pricepos_25"] = df["pricepos_25"].astype(float).replace([np.inf, -np.inf], np.nan) # NaN = no competitor202 df["has_coords"] = df["t"].notna()203 y = make_targets(df["entry_ts"].values, df["first_cut_ts"].values, df["duration_days"].values, df["exit_confirmed"].values)204 for k, v in y.items():205 df[k] = v206 df["strict_observable"] = df["entry_ts"] <= STRICT_ENTRY_MAX_TS207 df["naive_observable"] = df["entry_day"] <= T_LAST - int(HORIZON)208 return df209210211# --------------------------------------------------------------------------------------212# models213# --------------------------------------------------------------------------------------214def _design_logistic(tr: pd.DataFrame, te: pd.DataFrame, num: list[str]) -> tuple[np.ndarray, np.ndarray, list[str]]:215 """Standardised numerics (train statistics, NaN -> train median + missing216 indicator) + one-hot of coarse categoricals with training levels."""217 cols, Xtr, Xte = [], [], []218 for c in num:219 a, b = tr[c].astype(float).values.copy(), te[c].astype(float).values.copy()220 miss_a, miss_b = np.isnan(a), np.isnan(b)221 if miss_a.any() or miss_b.any():222 med = np.nanmedian(a); a[miss_a] = med; b[miss_b] = med223 Xtr.append(miss_a.astype(float)); Xte.append(miss_b.astype(float)); cols.append(f"{c}_missing")224 mu, sd = a.mean(), a.std()225 sd = sd if sd > 0 else 1.0226 Xtr.append((a - mu) / sd); Xte.append((b - mu) / sd); cols.append(c)227 for c in CAT_COARSE:228 levels = sorted(tr[c].unique())[1:] # first level = reference229 for lv in levels:230 Xtr.append((tr[c].values == lv).astype(float)); Xte.append((te[c].values == lv).astype(float)); cols.append(f"{c}={lv}")231 return np.column_stack(Xtr), np.column_stack(Xte), cols232233234def fit_logistic(tr: pd.DataFrame, te: pd.DataFrame, y: str, num: list[str]) -> np.ndarray:235 Xtr, Xte, _ = _design_logistic(tr, te, num)236 m = LogisticRegression(C=1.0, max_iter=2000, solver="lbfgs") # L2 (sklearn default)237 m.fit(Xtr, tr[y].values.astype(int))238 return m.predict_proba(Xte)[:, 1]239240241def _lgb_frame(tr: pd.DataFrame, te: pd.DataFrame, num: list[str]) -> tuple[pd.DataFrame, pd.DataFrame, list[str]]:242 feats = num + CAT_ALL243 Ztr, Zte = tr[feats].copy(), te[feats].copy()244 for c in CAT_ALL:245 cats = pd.Index(sorted(Ztr[c].unique()))246 Ztr[c] = pd.Categorical(Ztr[c], categories=cats); Zte[c] = pd.Categorical(Zte[c], categories=cats)247 return Ztr, Zte, feats248249250def fit_lgbm(tr: pd.DataFrame, te: pd.DataFrame, y: str, num: list[str], inner: str) -> tuple[np.ndarray, pd.Series, int]:251 """LightGBM with early stopping on an inner validation split carved from the252 training fold only: `inner="chrono"` = latest 15 % of training entries;253 otherwise 15 % of training dealers (seed 36)."""254 import lightgbm as lgb255 if inner == "chrono":256 cut = np.quantile(tr["entry_ts"], 0.85); va = tr["entry_ts"].values >= cut257 else:258 rng = np.random.default_rng(SEED); dl = np.array(tr["dealer_id"].unique(), dtype=object); rng.shuffle(dl)259 hold = set(dl[: max(1, int(0.15 * len(dl)))]); va = tr["dealer_id"].isin(hold).values260 Ztr, Zte, feats = _lgb_frame(tr, te, num)261 ytr = tr[y].values.astype(int)262 m = lgb.LGBMClassifier(**LGB_PARAMS)263 m.fit(Ztr[~va], ytr[~va], eval_X=Ztr[va], eval_y=ytr[va], eval_metric="binary_logloss",264 categorical_feature=CAT_ALL, callbacks=[lgb.early_stopping(100, verbose=False)])265 p = m.predict_proba(Zte, num_iteration=m.best_iteration_)[:, 1]266 imp = pd.Series(m.booster_.feature_importance(importance_type="gain"), index=feats)267 imp = imp / imp.sum() if imp.sum() > 0 else imp268 return p, imp, int(m.best_iteration_ or 0)269270271def run_split(df: pd.DataFrame, tr_mask: np.ndarray, te_mask: np.ndarray, target: str, inner: str,272 use_dealer_feats: bool = True) -> dict:273 """Fit both models on one train/test partition; returns predictions + metadata."""274 y = TARGETS[target]275 tr = df[tr_mask & df[y].notna().values].copy(); te = df[te_mask & df[y].notna().values].copy()276 num = list(NUM_BASE)277 if use_dealer_feats:278 for tgt in ("markdown30", "exit30"):279 a, b = target_encode(tr["dealer_id"].values, tr[TARGETS[tgt]].values, te["dealer_id"].values)280 tr[f"dealer_rate_{tgt}"] = a; te[f"dealer_rate_{tgt}"] = b281 cnt = tr["dealer_id"].value_counts()282 tr["dealer_log_n_train"] = np.log1p(tr["dealer_id"].map(cnt).fillna(0).values)283 te["dealer_log_n_train"] = np.log1p(te["dealer_id"].map(cnt).fillna(0).values)284 num += DEALER_FEATS285 p_lr = fit_logistic(tr, te, y, num)286 p_gb, imp, best_it = fit_lgbm(tr, te, y, num, inner)287 return {"y": te[y].values.astype(float), "p_logistic": p_lr, "p_lightgbm": p_gb, "importance": imp, "best_iteration": best_it,288 "n_train": int(len(tr)), "n_test": int(len(te)), "base_rate_train": float(tr[y].mean()),289 "test_index": te.index.values, "region": te["region"].values}290291292# --------------------------------------------------------------------------------------293# survival (Target D, exploratory)294# --------------------------------------------------------------------------------------295def survival_cox(df: pd.DataFrame) -> dict:296 """lifelines CoxPHFitter for time-to-exit on all S2 listings; chronological split297 at 2026-10-01 (censoring handled by the likelihood, so no observability rule)."""298 from lifelines import CoxPHFitter299 from lifelines.utils import concordance_index300 d = df.copy()301 d["pricepos_missing"] = d["pricepos_25"].isna().astype(float); d["pricepos_25f"] = d["pricepos_25"].fillna(0.0)302 d["log_dens_25"] = d["log_dens_25"].fillna(d["log_dens_25"].median()); d["coords_missing"] = (~d["has_coords"]).astype(float)303 cols = ["age_years_entry", "logp", "resid_A2", "km10k", "log_dens_25", "pricepos_25f", "pricepos_missing", "coords_missing"]304 X = d[cols].copy()305 for c, ref in [("powertrain", "ICE"), ("body", "SUV"), ("region", "Montréal")]:306 for lv in sorted(d[c].unique()):307 if lv != ref:308 X[f"{c}={lv}"] = (d[c] == lv).astype(float)309 X["T"] = d["duration_days"].values; X["E"] = d["exit_confirmed"].astype(int).values310 tr = (d["entry_ts"] < OCT1_TS).values; te = ~tr311 Xtr, Xte = X[tr], X[te]312 # drop test-only zero-variance dummies from both frames (fitting requires variation)313 keep = [c for c in X.columns if c in ("T", "E") or Xtr[c].std() > 0]314 Xtr, Xte = Xtr[keep], Xte[keep]315 t0 = time.time()316 cph = CoxPHFitter(penalizer=0.01)317 cph.fit(Xtr, duration_col="T", event_col="E", robust=False)318 ph_te = cph.predict_partial_hazard(Xte).values; ph_tr = cph.predict_partial_hazard(Xtr).values319 out = {"model": "CoxPHFitter(penalizer=0.01)", "split": "chronological, entry before vs on/after 2026-10-01 UTC",320 "n_train": int(len(Xtr)), "events_train": int(Xtr["E"].sum()), "n_test": int(len(Xte)), "events_test": int(Xte["E"].sum()),321 "c_index_train": float(concordance_index(Xtr["T"], -ph_tr, Xtr["E"])),322 "c_index_test": float(concordance_index(Xte["T"], -ph_te, Xte["E"])),323 "dynamic_auc_30d": None,324 "dynamic_auc_note": "scikit-survival (sksurv) is not installed in the project environment; the 30-day "325 "cumulative/dynamic AUC was not computed.",326 "covariates": [c for c in keep if c not in ("T", "E")], "seconds": round(time.time() - t0, 1)}327 summ = cph.summary.reset_index().rename(columns={"covariate": "term"})328 summ.to_csv(OUT / "cox_exit_summary.csv", index=False)329 return out330331332# --------------------------------------------------------------------------------------333# orchestration334# --------------------------------------------------------------------------------------335def _store(res: dict, target: str, split: str, store: dict, cal_tag: str | None = None, groups=None) -> None:336 rec = {"n_train": res["n_train"], "n_test": res["n_test"], "base_rate_train": res["base_rate_train"], "models": {}}337 for mname in ("logistic", "lightgbm"):338 p = res[f"p_{mname}"]339 rec["models"][mname] = bootstrap_metrics(res["y"], p, groups=groups)340 if cal_tag:341 calibration_table(res["y"], p).to_csv(OUT / f"calibration_{target}_{cal_tag}_{mname}.csv", index=False)342 rec["models"]["lightgbm"]["best_iteration"] = res.get("best_iteration")343 store.setdefault(target, {})[split] = rec344345346def _pool(parts: list[dict]) -> dict:347 return {"y": np.concatenate([p["y"] for p in parts]), "p_logistic": np.concatenate([p["p_logistic"] for p in parts]),348 "p_lightgbm": np.concatenate([p["p_lightgbm"] for p in parts]), "region": np.concatenate([p["region"] for p in parts]),349 "n_train": int(np.mean([p["n_train"] for p in parts])), "n_test": int(sum(p["n_test"] for p in parts)),350 "base_rate_train": float(np.mean([p["base_rate_train"] for p in parts])),351 "best_iteration": int(np.mean([p["best_iteration"] for p in parts]))}352353354def run_all(n_boot: int = N_BOOT) -> dict:355 global N_BOOT356 N_BOOT = n_boot357 OUT.mkdir(parents=True, exist_ok=True)358 t_start = time.time()359 df = build_frame()360 S = df[df["strict_observable"]].copy()361 log(f"S2 listings {len(df):,}; strictly observable 30-day sample {len(S):,} (entry day <= {int(S['entry_day'].max())})")362 metrics: dict = {"targets": {}, "notes": {}, "samples": {}}363 metrics["samples"] = {364 "n_S2": int(len(df)), "n_strict_observable": int(len(S)),365 "strict_rule": f"entry_ts <= SNAPSHOT_TS - {int(HORIZON)} - {C.EXIT_CONFIRM_DAYS} days (entry day <= {int(S['entry_day'].max())}); "366 "every listing has a determined 30-day exit and markdown outcome",367 "n_naive_observable": int(df["naive_observable"].sum()),368 "naive_rule": f"entry day <= {T_LAST - int(HORIZON)} (analysis-plan rule; exit at ages > snapshot-7d-entry cannot be confirmed)",369 "n_without_coordinates": int((~S["has_coords"]).sum()),370 "base_rates_strict": {k: float(S[v].mean()) for k, v in TARGETS.items()},371 "n_markdown30cr": int(S["y_markdown30cr"].notna().sum()),372 "n_exit_before_cut_removed_in_cr": int(S["y_markdown30cr"].isna().sum()),373 }374 cal_rows = []375 # ---- (i) chronological ------------------------------------------------------------376 tr = (S["entry_ts"] < CHRONO_CUTOFF_TS).values; te = ~tr377 for target in TARGETS:378 r = run_split(S, tr, te, target, inner="chrono")379 log(f"[chrono] {target}: train {r['n_train']:,} / test {r['n_test']:,}")380 _store(r, target, "chronological", metrics["targets"], cal_tag="chronological")381 r["importance"].sort_values(ascending=False).rename("gain_share").to_csv(OUT / f"feature_importance_{target}_chronological.csv", header=True)382 cal_rows.append((target, r))383 log(f" logistic {metrics['targets'][target]['chronological']['models']['logistic']['log_loss']:.4f} "384 f"lightgbm {metrics['targets'][target]['chronological']['models']['lightgbm']['log_loss']:.4f} (best_it {r['best_iteration']})")385 # ---- (ii) dealer holdout: 5-fold GroupKFold by dealer (dealers shuffled, seed 36) -----------------386 rng = np.random.default_rng(SEED); dealers = np.array(S["dealer_id"].unique(), dtype=object); rng.shuffle(dealers)387 fold_of = {d: i % 5 for i, d in enumerate(dealers)}; fold = S["dealer_id"].map(fold_of).values388 for target in TARGETS:389 parts, fold_metrics, imps = [], [], []390 for f in range(5):391 r = run_split(S, fold != f, fold == f, target, inner="group")392 parts.append(r); imps.append(r["importance"])393 fold_metrics.append({m: _metrics(r["y"], r[f"p_{m}"]) | {"n_test": r["n_test"]} for m in ("logistic", "lightgbm")})394 log(f"[dealer] {target} fold {f}: logistic {fold_metrics[-1]['logistic']['log_loss']:.4f} lightgbm {fold_metrics[-1]['lightgbm']['log_loss']:.4f}")395 pooled = _pool(parts)396 _store(pooled, target, "dealer_holdout", metrics["targets"], cal_tag="dealer_holdout")397 metrics["targets"][target]["dealer_holdout"]["folds"] = fold_metrics398 metrics["targets"][target]["dealer_holdout"]["fold_mean"] = {m: {k: float(np.mean([x[m][k] for x in fold_metrics])) for k in ("log_loss", "brier", "auc")} for m in ("logistic", "lightgbm")}399 metrics["targets"][target]["dealer_holdout"]["fold_sd"] = {m: {k: float(np.std([x[m][k] for x in fold_metrics])) for k in ("log_loss", "brier", "auc")} for m in ("logistic", "lightgbm")}400 pd.concat(imps, axis=1).mean(axis=1).sort_values(ascending=False).rename("gain_share").to_csv(OUT / f"feature_importance_{target}_dealer_holdout.csv", header=True)401 # ---- (iii) leave-one-region-out over the 6 largest regions ----------------------------402 big = [r for r in S["region"].value_counts().index if r != "(missing)"][:N_REGIONS_LOO]403 metrics["notes"]["region_holdout_regions"] = big404 for target in TARGETS:405 parts, per_region = [], {}406 for reg in big:407 te = (S["region"] == reg).values408 r = run_split(S, ~te, te, target, inner="group")409 parts.append(r)410 per_region[reg] = {m: _metrics(r["y"], r[f"p_{m}"]) | {"n_test": r["n_test"], "base_rate": float(r["y"].mean())} for m in ("logistic", "lightgbm")}411 log(f"[region] {target} {reg}: n {r['n_test']:,} logistic {per_region[reg]['logistic']['log_loss']:.4f} lightgbm {per_region[reg]['lightgbm']['log_loss']:.4f}")412 pooled = _pool(parts)413 _store(pooled, target, "region_holdout", metrics["targets"], cal_tag="region_holdout", groups=pooled["region"])414 metrics["targets"][target]["region_holdout"]["per_region"] = per_region415 metrics["targets"][target]["region_holdout"]["aggregation"] = "unweighted mean over the 6 held-out regions; CI = bootstrap within region, then averaged"416 # ---- (iv) drift: train weeks 0-3, evaluate weeks 4-5 vs week 6 (observable part) ----------------417 tr = (S["week"] <= 3).values418 for target in TARGETS:419 drift = {}420 for name, lo, hi in [("weeks_4_5", 4, 5), ("weeks_6_7", 6, 7)]:421 te = ((S["week"] >= lo) & (S["week"] <= hi)).values422 if te.sum() < 200:423 drift[name] = {"n_test": int(te.sum()), "skipped": "too few observable test rows"}; continue424 r = run_split(S, tr, te, target, inner="chrono")425 drift[name] = {"n_train": r["n_train"], "n_test": r["n_test"], "weeks_observable": sorted(int(w) for w in S.loc[te, "week"].unique()),426 "models": {m: bootstrap_metrics(r["y"], r[f"p_{m}"]) for m in ("logistic", "lightgbm")}}427 log(f"[drift] {target} {name}: n {r['n_test']:,} logistic {drift[name]['models']['logistic']['log_loss']:.4f} lightgbm {drift[name]['models']['lightgbm']['log_loss']:.4f}")428 metrics["targets"][target]["drift"] = drift429 metrics["notes"]["drift"] = ("train = entries in weeks 0-3 (days 0-27); the week 6-7 window only contains observable entries "430 f"up to day {int(S['entry_day'].max())} (week 6), week 7 has no determinable 30-day outcome")431 # ---- (v) sensitivity: plan's Oct-1 cutoff on the looser observability rule -----------------------432 N = df[df["naive_observable"]].copy()433 tr = (N["entry_ts"] < OCT1_TS).values; te = ~tr434 for target in ("markdown30", "exit30"):435 r = run_split(N, tr, te, target, inner="chrono")436 _store(r, target, "chronological_oct1_naive", metrics["targets"])437 log(f"[oct1-naive] {target}: train {r['n_train']:,} test {r['n_test']:,} base {r['y'].mean():.3f} "438 f"logistic {metrics['targets'][target]['chronological_oct1_naive']['models']['logistic']['log_loss']:.4f} "439 f"lightgbm {metrics['targets'][target]['chronological_oct1_naive']['models']['lightgbm']['log_loss']:.4f}")440 metrics["notes"]["chronological_oct1_naive"] = (441 "Analysis-plan cutoff (entry before vs on/after 2026-10-01 UTC) applied to the looser rule entry day <= 51. "442 "Test entries (days 50-51) have only ~22-24 days of confirmable follow-up: exits at ages between ~23 and 30 days are "443 "coded 0 (exit30 base rate is a lower bound) and the last-day markdowns may be unclassified. Reported as a sensitivity only.")444 metrics["notes"]["chronological"] = (f"primary chronological split: train entries before 2026-09-15 UTC, test 2026-09-15 to day "445 f"{int(S['entry_day'].max())} (2026-09-25); the plan's 2026-10-01 cutoff leaves no listing with a "446 "determinable 30-day outcome in the test fold")447 metrics["notes"]["targets"] = {448 "markdown30": "first sustained markdown (>=1 % and >=100 $) at listing age <= 30 d; exits before any cut count as 0",449 "markdown30cr": "competing-risk-aware: listings that exit within 30 d without a prior cut are removed; event = markdown within 30 d among the rest",450 "exit30": "disappearance confirmed by >= 7 d of absence at listing age <= 30 d (NOT a sale)",451 }452 metrics["notes"]["features"] = {453 "numeric": NUM_BASE, "dealer_target_encoding": DEALER_FEATS + ["(training fold only; inner 5-fold OOF for training rows; smoothing m=20)"],454 "categorical_lightgbm": CAT_ALL, "one_hot_logistic": CAT_COARSE,455 "competition_at_entry": "risk_set row with listing_age = 0 (dens_10/25/50, own_stock, pricepos_25; NaN when no competitor or no coordinates)",456 "excluded": "dealer_id (identity), anything measured after entry",457 }458 metrics["notes"]["models"] = {"logistic": "sklearn LogisticRegression L2, C=1, standardised numerics + missing indicators + one-hot coarse categoricals",459 "lightgbm": {k: v for k, v in LGB_PARAMS.items()} | {"early_stopping": "100 rounds on inner validation (latest 15 % of training entries for chronological splits, 15 % of training dealers otherwise)"}}460 metrics["notes"]["metrics"] = {"ci": f"{N_BOOT} bootstrap resamples of test rows, percentile 2.5/97.5", "calibration": f"{N_BINS} equal-width probability bins; ECE = n-weighted |frac_pos - mean_pred|"}461 # ---- (vi) survival ----------------------------------------------------------------------------462 log("[survival] Cox PH on S2 time-to-exit")463 metrics["survival"] = survival_cox(df)464 log(f" c-index train {metrics['survival']['c_index_train']:.4f} test {metrics['survival']['c_index_test']:.4f}")465 metrics["seconds_total"] = round(time.time() - t_start, 1)466 (OUT / "metrics.json").write_text(json.dumps(metrics, indent=2, default=float))467 make_figure(cal_rows, metrics)468 write_table(metrics)469 log(f"done in {metrics['seconds_total']} s")470 return metrics471472473# --------------------------------------------------------------------------------------474# figure + table475# --------------------------------------------------------------------------------------476def make_figure(cal_rows: list[tuple[str, dict]], metrics: dict) -> None:477 import matplotlib.pyplot as plt478 from . import plotstyle as PS479 PS.setup()480 fig, (a, b) = plt.subplots(1, 2, figsize=(7.0, 3.1), gridspec_kw={"width_ratios": [1, 1.35]})481 col = {"logistic": PS.CAT[0], "lightgbm": PS.CAT[1]}482 style = {"markdown30": ("-", "o", "Markdown within 30 d"), "exit30": ("--", "s", "Exit within 30 d")}483 a.plot([0, 1], [0, 1], color=PS.NEUTRAL, lw=0.8, zorder=1)484 for target, r in cal_rows:485 if target not in style:486 continue487 ls, mk, lab = style[target]488 for m in ("logistic", "lightgbm"):489 t = calibration_table(r["y"], r[f"p_{m}"]); t = t[t["n"] >= 20]490 a.plot(t["mean_pred"], t["frac_pos"], ls, marker=mk, color=col[m], lw=1.0, ms=3,491 label=f"{lab} — {'logistic' if m == 'logistic' else 'LightGBM'}")492 a.set_xlim(0, 1); a.set_ylim(0, 1); a.set_xlabel("Predicted probability"); a.set_ylabel("Observed frequency")493 a.set_title("(a) Reliability, chronological split"); a.legend(loc="upper left", fontsize=6.2)494 splits = [("chronological", "Chrono."), ("dealer_holdout", "Dealer"), ("region_holdout", "Region")]495 labels, xs = [], []496 x = 0.0497 for target, tlab in [("markdown30", "Markdown"), ("exit30", "Exit")]:498 for s, slab in splits:499 for j, m in enumerate(("logistic", "lightgbm")):500 mm = metrics["targets"][target][s]["models"][m]501 lo, hi = mm["log_loss_ci"]502 b.bar(x + j * 0.38, mm["log_loss"], width=0.36, color=col[m], label=("Logistic" if m == "logistic" else "LightGBM") if (target == "markdown30" and s == "chronological") else None)503 b.errorbar(x + j * 0.38, mm["log_loss"], yerr=[[mm["log_loss"] - lo], [hi - mm["log_loss"]]], color=PS.TEXT, lw=0.7, capsize=2)504 labels.append(f"{slab}\n{tlab}"); xs.append(x + 0.19); x += 1.2505 x += 0.35506 b.set_xticks(xs); b.set_xticklabels(labels, fontsize=6.5); b.set_ylabel("Log loss (test)")507 b.set_title("(b) Out-of-sample log loss by split (30-day targets)"); b.legend(loc="upper left", fontsize=7)508 b.grid(axis="x", visible=False)509 PS.save(fig, "fig10_prediction")510511512def _fmt_int(x) -> str:513 return f"{int(round(float(x))):,}".replace(",", "\\,")514515516def write_table(metrics: dict) -> None:517 rows = ["\\begin{tabular}{llllrrrr}", "\\toprule",518 "Target & Split & Model & Log loss [95\\% CI] & Brier & AUC & $n_{\\text{test}}$ & Base rate \\\\", "\\midrule"]519 tl = {"markdown30": "Markdown $\\le$ 30 d", "exit30": "Exit $\\le$ 30 d"}520 sl = {"chronological": "Chronological", "dealer_holdout": "Dealer holdout", "region_holdout": "Region LORO"}521 for target in ("markdown30", "exit30"):522 for s in ("chronological", "dealer_holdout", "region_holdout"):523 rec = metrics["targets"][target][s]524 for k, m in enumerate(("logistic", "lightgbm")):525 mm = rec["models"][m]526 lead = f"{tl[target]} & {sl[s]}" if k == 0 else " & "527 rows.append(f"{lead} & {'Logistic' if m == 'logistic' else 'LightGBM'} & {mm['log_loss']:.3f} [{mm['log_loss_ci'][0]:.3f}, {mm['log_loss_ci'][1]:.3f}] "528 f"& {mm['brier']:.3f} & {mm['auc']:.3f} & {_fmt_int(mm['n_test'])} & {mm['base_rate']:.3f} \\\\")529 if target == "markdown30":530 rows.append("\\addlinespace")531 rows += ["\\bottomrule", "\\end{tabular}"]532 C.TABLES.mkdir(parents=True, exist_ok=True)533 (C.TABLES / "tab08_prediction.tex").write_text("\n".join(rows) + "\n", encoding="utf-8")534