"""Model F — prediction and economic value (analysis_plan §3 "Model F"). Targets on the S2 dynamic sample, all defined at a 30-day horizon from entry: B markdown30 first sustained markdown at listing age <= 30 days (listings that exit before any cut count as 0); B' markdown30cr competing-risk-aware variant: same event, but listings that exit within 30 days *without* a prior cut are removed from the sample (markdown among listings still exposed to the markdown risk); C exit30 observed disappearance confirmed by >= 7 days of absence at listing age <= 30 days (an exit is NOT a sale). Target A (log entry price) is covered by the hedonic grouped CV (hedonic.grouped_cv) and is not redone here. Observability. An exit at age a is confirmable only if entry_ts + a + 7 d <= snapshot; the 30-day outcome of a listing is therefore determined only when entry_ts <= SNAPSHOT_TS - 37 d (entry day <= 43, 2026-09-25 UTC). The analysis plan's chronological cutoff (2026-10-01) leaves no such listing in the test fold, so the primary chronological split uses 2026-09-15 UTC; the Oct-1 cutoff is reported as a flagged sensitivity on the looser rule "entry day <= 81 - 30" (exit labels partly right-censored there, see metrics.json notes). Features are measured AT ENTRY only. dealer_id is never a feature; dealer-level outcome rates are target-encoded on the training fold only (inner 5-fold out-of-fold encoding for training rows, smoothed towards the training prior). """ from __future__ import annotations import datetime as dt import json import math import time from pathlib import Path import numpy as np import pandas as pd import polars as pl from sklearn.linear_model import LogisticRegression from sklearn.metrics import brier_score_loss, log_loss, roc_auc_score from . import config as C OUT = C.RESULTS / "prediction" HORIZON = 30.0 SEED = 36 N_BOOT = 200 N_BINS = 10 DAY0 = math.floor(C.WINDOW_START_TS / 86400) CHRONO_CUTOFF_TS = dt.datetime(2026, 9, 15, tzinfo=dt.timezone.utc).timestamp() OCT1_TS = dt.datetime(2026, 10, 1, tzinfo=dt.timezone.utc).timestamp() T_LAST = int(math.floor(C.SNAPSHOT_TS / 86400) - DAY0) # 81 STRICT_ENTRY_MAX_TS = C.SNAPSHOT_TS - (HORIZON + C.EXIT_CONFIRM_DAYS) * 86400 N_REGIONS_LOO = 6 N_THREADS = 8 NUM_BASE = ["age_years_entry", "km10k", "year", "logp", "resid_A2", "week", "log_dens_10", "log_dens_25", "log_dens_50", "log_own_stock", "pricepos_25"] CAT_COARSE = ["powertrain", "body", "region", "source_class"] CAT_FINE = ["make", "model", "trim", "transmission", "drivetrain"] CAT_ALL = CAT_COARSE + CAT_FINE DEALER_FEATS = ["dealer_rate_markdown30", "dealer_rate_exit30", "dealer_log_n_train"] TARGETS = {"markdown30": "y_markdown30", "markdown30cr": "y_markdown30cr", "exit30": "y_exit30"} LGB_PARAMS = dict(n_estimators=2000, learning_rate=0.05, num_leaves=31, min_child_samples=50, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=1.0, cat_l2=10.0, cat_smooth=50.0, max_cat_to_onehot=4, n_jobs=N_THREADS, verbose=-1, random_state=SEED) def log(*a) -> None: print(time.strftime("%H:%M:%S"), *a, flush=True) # -------------------------------------------------------------------------------------- # pure helpers (unit-tested on synthetic arrays) # -------------------------------------------------------------------------------------- def make_targets(entry_ts: np.ndarray, first_cut_ts: np.ndarray, duration_days: np.ndarray, exit_confirmed: np.ndarray, horizon: float = HORIZON) -> dict[str, np.ndarray]: """30-day outcome labels from listing-level timestamps. first_cut_ts may be NaN (no sustained cut). Returns float arrays; markdown30cr is NaN for listings excluded by the competing-risk rule (exit before any cut within the horizon).""" entry_ts = np.asarray(entry_ts, float); fc = np.asarray(first_cut_ts, float) dur = np.asarray(duration_days, float); ex = np.asarray(exit_confirmed, bool) cut_age = (fc - entry_ts) / 86400.0 has_cut_h = np.isfinite(cut_age) & (cut_age <= horizon) exit_h = ex & (dur <= horizon) y_b = has_cut_h.astype(float) y_c = exit_h.astype(float) exit_before_cut = exit_h & ~has_cut_h y_bcr = np.where(exit_before_cut, np.nan, y_b) return {"y_markdown30": y_b, "y_markdown30cr": y_bcr, "y_exit30": y_c} def calibration_table(y: np.ndarray, p: np.ndarray, n_bins: int = N_BINS) -> pd.DataFrame: """Reliability table on equal-width probability bins [0,1].""" y = np.asarray(y, float); p = np.asarray(p, float) edges = np.linspace(0, 1, n_bins + 1) b = np.clip(np.digitize(p, edges[1:-1], right=False), 0, n_bins - 1) rows = [] for k in range(n_bins): m = b == k rows.append({"bin": k, "lo": round(float(edges[k]), 6), "hi": round(float(edges[k + 1]), 6), "n": int(m.sum()), "mean_pred": float(p[m].mean()) if m.any() else np.nan, "frac_pos": float(y[m].mean()) if m.any() else np.nan}) return pd.DataFrame(rows) def expected_calibration_error(y: np.ndarray, p: np.ndarray, n_bins: int = N_BINS) -> float: t = calibration_table(y, p, n_bins) t = t[t["n"] > 0] return float((t["n"] / t["n"].sum() * (t["frac_pos"] - t["mean_pred"]).abs()).sum()) def _metrics(y, p) -> dict: p = np.clip(p, 1e-6, 1 - 1e-6) auc = float(roc_auc_score(y, p)) if 0 < y.mean() < 1 else float("nan") return {"log_loss": float(log_loss(y, p, labels=[0, 1])), "brier": float(brier_score_loss(y, p)), "auc": auc} def bootstrap_metrics(y: np.ndarray, p: np.ndarray, groups: np.ndarray | None = None, n_boot: int = N_BOOT, seed: int = SEED) -> dict: """Point estimates + percentile bootstrap CIs (resampling test rows). When `groups` is given, rows are resampled within each group and the metric is the unweighted mean over groups (used for leave-one-region-out).""" y = np.asarray(y, float); p = np.asarray(p, float) rng = np.random.default_rng(seed) if groups is None: point = _metrics(y, p) idx_sets = [None] else: groups = np.asarray(groups) idx_sets = [np.flatnonzero(groups == g) for g in np.unique(groups)] point = {k: float(np.mean([_metrics(y[ix], p[ix])[k] for ix in idx_sets])) for k in ("log_loss", "brier", "auc")} draws = {k: [] for k in point} for _ in range(n_boot): if groups is None: ix = rng.integers(0, len(y), len(y)); m = _metrics(y[ix], p[ix]) else: ms = [] for ix in idx_sets: jx = ix[rng.integers(0, len(ix), len(ix))]; ms.append(_metrics(y[jx], p[jx])) m = {k: float(np.mean([x[k] for x in ms])) for k in point} for k in point: draws[k].append(m[k]) out = {} for k in point: d = np.asarray(draws[k], float) out[k] = point[k]; out[f"{k}_ci"] = [float(np.nanpercentile(d, 2.5)), float(np.nanpercentile(d, 97.5))] out["ece"] = expected_calibration_error(y, p) out["n_test"] = int(len(y)); out["base_rate"] = float(y.mean()); out["n_boot"] = int(n_boot) return out def target_encode(groups_tr: np.ndarray, y_tr: np.ndarray, groups_te: np.ndarray, m: float = 20.0, k: int = 5, seed: int = SEED) -> tuple[np.ndarray, np.ndarray]: """Smoothed group mean of y computed on the training fold only. Training rows get an inner k-fold out-of-fold encoding (their own label never enters their feature); test rows get the full-training encoding; unseen groups get the training prior.""" g_tr = pd.Series(np.asarray(groups_tr)); y_tr = np.asarray(y_tr, float); g_te = pd.Series(np.asarray(groups_te)) prior = float(np.nanmean(y_tr)) ok = ~np.isnan(y_tr) def enc(gs: pd.Series, ys: np.ndarray, target: pd.Series) -> np.ndarray: d = pd.DataFrame({"g": gs.values, "y": ys}).dropna() agg = d.groupby("g")["y"].agg(["sum", "count"]) s = target.map(agg["sum"]).fillna(0.0).values; n = target.map(agg["count"]).fillna(0.0).values return (s + m * prior) / (n + m) enc_tr = np.full(len(g_tr), prior) rng = np.random.default_rng(seed); fold = rng.integers(0, k, len(g_tr)) for f in range(k): tr_in = (fold != f) & ok; te_in = fold == f enc_tr[te_in] = enc(g_tr[tr_in], y_tr[tr_in], g_tr[te_in]) enc_te = enc(g_tr[ok], y_tr[ok], g_te) return enc_tr, enc_te # -------------------------------------------------------------------------------------- # data # -------------------------------------------------------------------------------------- def build_frame() -> pd.DataFrame: """S2 listings with entry-time features, risk-set competition on the entry day and the 30-day labels (all S2 entries; observability masks are applied per split).""" L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s2_dynamic")) L = L.select(["listing_id", "dealer_id", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km", "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "duration_days", "exit_confirmed", "first_cut_ts", "age_years_entry"]) R = pl.read_parquet(C.PROCESSED / "risk_set.parquet", columns=["listing_id", "t", "listing_age", "dens_10", "dens_25", "dens_50", "own_stock", "pricepos_25"]) first = R.filter(pl.col("listing_age") == 0).select(["listing_id", "t", "dens_10", "dens_25", "dens_50", "own_stock", "pricepos_25"]) res = pl.read_parquet(C.RESULTS / "hedonic" / "residuals.parquet").select(["listing_id", "resid_A2"]) df = L.join(first, on="listing_id", how="left").join(res, on="listing_id", how="left").to_pandas() df = df.rename(columns={"body_filled": "body"}) for c, fill in [("body", "unknown"), ("trim", "(none)"), ("region", "(missing)"), ("transmission", "unknown"), ("drivetrain", "unknown"), ("powertrain", "unknown")]: df[c] = df[c].fillna(fill).astype(str) df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60) df["logp"] = np.log(df["entry_price"]); df["km10k"] = df["mileage_km"] / 1e4 df["entry_day"] = (np.floor(df["entry_ts"] / 86400) - DAY0).astype(int) df["week"] = df["entry_day"] // 7 for r in ("10", "25", "50"): df[f"log_dens_{r}"] = np.log1p(df[f"dens_{r}"].astype(float)) df["log_own_stock"] = np.log1p(df["own_stock"].astype(float)) df["pricepos_25"] = df["pricepos_25"].astype(float).replace([np.inf, -np.inf], np.nan) # NaN = no competitor df["has_coords"] = df["t"].notna() y = make_targets(df["entry_ts"].values, df["first_cut_ts"].values, df["duration_days"].values, df["exit_confirmed"].values) for k, v in y.items(): df[k] = v df["strict_observable"] = df["entry_ts"] <= STRICT_ENTRY_MAX_TS df["naive_observable"] = df["entry_day"] <= T_LAST - int(HORIZON) return df # -------------------------------------------------------------------------------------- # models # -------------------------------------------------------------------------------------- def _design_logistic(tr: pd.DataFrame, te: pd.DataFrame, num: list[str]) -> tuple[np.ndarray, np.ndarray, list[str]]: """Standardised numerics (train statistics, NaN -> train median + missing indicator) + one-hot of coarse categoricals with training levels.""" cols, Xtr, Xte = [], [], [] for c in num: a, b = tr[c].astype(float).values.copy(), te[c].astype(float).values.copy() miss_a, miss_b = np.isnan(a), np.isnan(b) if miss_a.any() or miss_b.any(): med = np.nanmedian(a); a[miss_a] = med; b[miss_b] = med Xtr.append(miss_a.astype(float)); Xte.append(miss_b.astype(float)); cols.append(f"{c}_missing") mu, sd = a.mean(), a.std() sd = sd if sd > 0 else 1.0 Xtr.append((a - mu) / sd); Xte.append((b - mu) / sd); cols.append(c) for c in CAT_COARSE: levels = sorted(tr[c].unique())[1:] # first level = reference for lv in levels: Xtr.append((tr[c].values == lv).astype(float)); Xte.append((te[c].values == lv).astype(float)); cols.append(f"{c}={lv}") return np.column_stack(Xtr), np.column_stack(Xte), cols def fit_logistic(tr: pd.DataFrame, te: pd.DataFrame, y: str, num: list[str]) -> np.ndarray: Xtr, Xte, _ = _design_logistic(tr, te, num) m = LogisticRegression(C=1.0, max_iter=2000, solver="lbfgs") # L2 (sklearn default) m.fit(Xtr, tr[y].values.astype(int)) return m.predict_proba(Xte)[:, 1] def _lgb_frame(tr: pd.DataFrame, te: pd.DataFrame, num: list[str]) -> tuple[pd.DataFrame, pd.DataFrame, list[str]]: feats = num + CAT_ALL Ztr, Zte = tr[feats].copy(), te[feats].copy() for c in CAT_ALL: cats = pd.Index(sorted(Ztr[c].unique())) Ztr[c] = pd.Categorical(Ztr[c], categories=cats); Zte[c] = pd.Categorical(Zte[c], categories=cats) return Ztr, Zte, feats def fit_lgbm(tr: pd.DataFrame, te: pd.DataFrame, y: str, num: list[str], inner: str) -> tuple[np.ndarray, pd.Series, int]: """LightGBM with early stopping on an inner validation split carved from the training fold only: `inner="chrono"` = latest 15 % of training entries; otherwise 15 % of training dealers (seed 36).""" import lightgbm as lgb if inner == "chrono": cut = np.quantile(tr["entry_ts"], 0.85); va = tr["entry_ts"].values >= cut else: rng = np.random.default_rng(SEED); dl = np.array(tr["dealer_id"].unique(), dtype=object); rng.shuffle(dl) hold = set(dl[: max(1, int(0.15 * len(dl)))]); va = tr["dealer_id"].isin(hold).values Ztr, Zte, feats = _lgb_frame(tr, te, num) ytr = tr[y].values.astype(int) m = lgb.LGBMClassifier(**LGB_PARAMS) m.fit(Ztr[~va], ytr[~va], eval_X=Ztr[va], eval_y=ytr[va], eval_metric="binary_logloss", categorical_feature=CAT_ALL, callbacks=[lgb.early_stopping(100, verbose=False)]) p = m.predict_proba(Zte, num_iteration=m.best_iteration_)[:, 1] imp = pd.Series(m.booster_.feature_importance(importance_type="gain"), index=feats) imp = imp / imp.sum() if imp.sum() > 0 else imp return p, imp, int(m.best_iteration_ or 0) def run_split(df: pd.DataFrame, tr_mask: np.ndarray, te_mask: np.ndarray, target: str, inner: str, use_dealer_feats: bool = True) -> dict: """Fit both models on one train/test partition; returns predictions + metadata.""" y = TARGETS[target] tr = df[tr_mask & df[y].notna().values].copy(); te = df[te_mask & df[y].notna().values].copy() num = list(NUM_BASE) if use_dealer_feats: for tgt in ("markdown30", "exit30"): a, b = target_encode(tr["dealer_id"].values, tr[TARGETS[tgt]].values, te["dealer_id"].values) tr[f"dealer_rate_{tgt}"] = a; te[f"dealer_rate_{tgt}"] = b cnt = tr["dealer_id"].value_counts() tr["dealer_log_n_train"] = np.log1p(tr["dealer_id"].map(cnt).fillna(0).values) te["dealer_log_n_train"] = np.log1p(te["dealer_id"].map(cnt).fillna(0).values) num += DEALER_FEATS p_lr = fit_logistic(tr, te, y, num) p_gb, imp, best_it = fit_lgbm(tr, te, y, num, inner) return {"y": te[y].values.astype(float), "p_logistic": p_lr, "p_lightgbm": p_gb, "importance": imp, "best_iteration": best_it, "n_train": int(len(tr)), "n_test": int(len(te)), "base_rate_train": float(tr[y].mean()), "test_index": te.index.values, "region": te["region"].values} # -------------------------------------------------------------------------------------- # survival (Target D, exploratory) # -------------------------------------------------------------------------------------- def survival_cox(df: pd.DataFrame) -> dict: """lifelines CoxPHFitter for time-to-exit on all S2 listings; chronological split at 2026-10-01 (censoring handled by the likelihood, so no observability rule).""" from lifelines import CoxPHFitter from lifelines.utils import concordance_index d = df.copy() d["pricepos_missing"] = d["pricepos_25"].isna().astype(float); d["pricepos_25f"] = d["pricepos_25"].fillna(0.0) d["log_dens_25"] = d["log_dens_25"].fillna(d["log_dens_25"].median()); d["coords_missing"] = (~d["has_coords"]).astype(float) cols = ["age_years_entry", "logp", "resid_A2", "km10k", "log_dens_25", "pricepos_25f", "pricepos_missing", "coords_missing"] X = d[cols].copy() for c, ref in [("powertrain", "ICE"), ("body", "SUV"), ("region", "Montréal")]: for lv in sorted(d[c].unique()): if lv != ref: X[f"{c}={lv}"] = (d[c] == lv).astype(float) X["T"] = d["duration_days"].values; X["E"] = d["exit_confirmed"].astype(int).values tr = (d["entry_ts"] < OCT1_TS).values; te = ~tr Xtr, Xte = X[tr], X[te] # drop test-only zero-variance dummies from both frames (fitting requires variation) keep = [c for c in X.columns if c in ("T", "E") or Xtr[c].std() > 0] Xtr, Xte = Xtr[keep], Xte[keep] t0 = time.time() cph = CoxPHFitter(penalizer=0.01) cph.fit(Xtr, duration_col="T", event_col="E", robust=False) ph_te = cph.predict_partial_hazard(Xte).values; ph_tr = cph.predict_partial_hazard(Xtr).values out = {"model": "CoxPHFitter(penalizer=0.01)", "split": "chronological, entry before vs on/after 2026-10-01 UTC", "n_train": int(len(Xtr)), "events_train": int(Xtr["E"].sum()), "n_test": int(len(Xte)), "events_test": int(Xte["E"].sum()), "c_index_train": float(concordance_index(Xtr["T"], -ph_tr, Xtr["E"])), "c_index_test": float(concordance_index(Xte["T"], -ph_te, Xte["E"])), "dynamic_auc_30d": None, "dynamic_auc_note": "scikit-survival (sksurv) is not installed in the project environment; the 30-day " "cumulative/dynamic AUC was not computed.", "covariates": [c for c in keep if c not in ("T", "E")], "seconds": round(time.time() - t0, 1)} summ = cph.summary.reset_index().rename(columns={"covariate": "term"}) summ.to_csv(OUT / "cox_exit_summary.csv", index=False) return out # -------------------------------------------------------------------------------------- # orchestration # -------------------------------------------------------------------------------------- def _store(res: dict, target: str, split: str, store: dict, cal_tag: str | None = None, groups=None) -> None: rec = {"n_train": res["n_train"], "n_test": res["n_test"], "base_rate_train": res["base_rate_train"], "models": {}} for mname in ("logistic", "lightgbm"): p = res[f"p_{mname}"] rec["models"][mname] = bootstrap_metrics(res["y"], p, groups=groups) if cal_tag: calibration_table(res["y"], p).to_csv(OUT / f"calibration_{target}_{cal_tag}_{mname}.csv", index=False) rec["models"]["lightgbm"]["best_iteration"] = res.get("best_iteration") store.setdefault(target, {})[split] = rec def _pool(parts: list[dict]) -> dict: return {"y": np.concatenate([p["y"] for p in parts]), "p_logistic": np.concatenate([p["p_logistic"] for p in parts]), "p_lightgbm": np.concatenate([p["p_lightgbm"] for p in parts]), "region": np.concatenate([p["region"] for p in parts]), "n_train": int(np.mean([p["n_train"] for p in parts])), "n_test": int(sum(p["n_test"] for p in parts)), "base_rate_train": float(np.mean([p["base_rate_train"] for p in parts])), "best_iteration": int(np.mean([p["best_iteration"] for p in parts]))} def run_all(n_boot: int = N_BOOT) -> dict: global N_BOOT N_BOOT = n_boot OUT.mkdir(parents=True, exist_ok=True) t_start = time.time() df = build_frame() S = df[df["strict_observable"]].copy() log(f"S2 listings {len(df):,}; strictly observable 30-day sample {len(S):,} (entry day <= {int(S['entry_day'].max())})") metrics: dict = {"targets": {}, "notes": {}, "samples": {}} metrics["samples"] = { "n_S2": int(len(df)), "n_strict_observable": int(len(S)), "strict_rule": f"entry_ts <= SNAPSHOT_TS - {int(HORIZON)} - {C.EXIT_CONFIRM_DAYS} days (entry day <= {int(S['entry_day'].max())}); " "every listing has a determined 30-day exit and markdown outcome", "n_naive_observable": int(df["naive_observable"].sum()), "naive_rule": f"entry day <= {T_LAST - int(HORIZON)} (analysis-plan rule; exit at ages > snapshot-7d-entry cannot be confirmed)", "n_without_coordinates": int((~S["has_coords"]).sum()), "base_rates_strict": {k: float(S[v].mean()) for k, v in TARGETS.items()}, "n_markdown30cr": int(S["y_markdown30cr"].notna().sum()), "n_exit_before_cut_removed_in_cr": int(S["y_markdown30cr"].isna().sum()), } cal_rows = [] # ---- (i) chronological ------------------------------------------------------------ tr = (S["entry_ts"] < CHRONO_CUTOFF_TS).values; te = ~tr for target in TARGETS: r = run_split(S, tr, te, target, inner="chrono") log(f"[chrono] {target}: train {r['n_train']:,} / test {r['n_test']:,}") _store(r, target, "chronological", metrics["targets"], cal_tag="chronological") r["importance"].sort_values(ascending=False).rename("gain_share").to_csv(OUT / f"feature_importance_{target}_chronological.csv", header=True) cal_rows.append((target, r)) log(f" logistic {metrics['targets'][target]['chronological']['models']['logistic']['log_loss']:.4f} " f"lightgbm {metrics['targets'][target]['chronological']['models']['lightgbm']['log_loss']:.4f} (best_it {r['best_iteration']})") # ---- (ii) dealer holdout: 5-fold GroupKFold by dealer (dealers shuffled, seed 36) ----------------- rng = np.random.default_rng(SEED); dealers = np.array(S["dealer_id"].unique(), dtype=object); rng.shuffle(dealers) fold_of = {d: i % 5 for i, d in enumerate(dealers)}; fold = S["dealer_id"].map(fold_of).values for target in TARGETS: parts, fold_metrics, imps = [], [], [] for f in range(5): r = run_split(S, fold != f, fold == f, target, inner="group") parts.append(r); imps.append(r["importance"]) fold_metrics.append({m: _metrics(r["y"], r[f"p_{m}"]) | {"n_test": r["n_test"]} for m in ("logistic", "lightgbm")}) log(f"[dealer] {target} fold {f}: logistic {fold_metrics[-1]['logistic']['log_loss']:.4f} lightgbm {fold_metrics[-1]['lightgbm']['log_loss']:.4f}") pooled = _pool(parts) _store(pooled, target, "dealer_holdout", metrics["targets"], cal_tag="dealer_holdout") metrics["targets"][target]["dealer_holdout"]["folds"] = fold_metrics 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")} 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")} 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) # ---- (iii) leave-one-region-out over the 6 largest regions ---------------------------- big = [r for r in S["region"].value_counts().index if r != "(missing)"][:N_REGIONS_LOO] metrics["notes"]["region_holdout_regions"] = big for target in TARGETS: parts, per_region = [], {} for reg in big: te = (S["region"] == reg).values r = run_split(S, ~te, te, target, inner="group") parts.append(r) 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")} 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}") pooled = _pool(parts) _store(pooled, target, "region_holdout", metrics["targets"], cal_tag="region_holdout", groups=pooled["region"]) metrics["targets"][target]["region_holdout"]["per_region"] = per_region metrics["targets"][target]["region_holdout"]["aggregation"] = "unweighted mean over the 6 held-out regions; CI = bootstrap within region, then averaged" # ---- (iv) drift: train weeks 0-3, evaluate weeks 4-5 vs week 6 (observable part) ---------------- tr = (S["week"] <= 3).values for target in TARGETS: drift = {} for name, lo, hi in [("weeks_4_5", 4, 5), ("weeks_6_7", 6, 7)]: te = ((S["week"] >= lo) & (S["week"] <= hi)).values if te.sum() < 200: drift[name] = {"n_test": int(te.sum()), "skipped": "too few observable test rows"}; continue r = run_split(S, tr, te, target, inner="chrono") 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()), "models": {m: bootstrap_metrics(r["y"], r[f"p_{m}"]) for m in ("logistic", "lightgbm")}} 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}") metrics["targets"][target]["drift"] = drift metrics["notes"]["drift"] = ("train = entries in weeks 0-3 (days 0-27); the week 6-7 window only contains observable entries " f"up to day {int(S['entry_day'].max())} (week 6), week 7 has no determinable 30-day outcome") # ---- (v) sensitivity: plan's Oct-1 cutoff on the looser observability rule ----------------------- N = df[df["naive_observable"]].copy() tr = (N["entry_ts"] < OCT1_TS).values; te = ~tr for target in ("markdown30", "exit30"): r = run_split(N, tr, te, target, inner="chrono") _store(r, target, "chronological_oct1_naive", metrics["targets"]) log(f"[oct1-naive] {target}: train {r['n_train']:,} test {r['n_test']:,} base {r['y'].mean():.3f} " f"logistic {metrics['targets'][target]['chronological_oct1_naive']['models']['logistic']['log_loss']:.4f} " f"lightgbm {metrics['targets'][target]['chronological_oct1_naive']['models']['lightgbm']['log_loss']:.4f}") metrics["notes"]["chronological_oct1_naive"] = ( "Analysis-plan cutoff (entry before vs on/after 2026-10-01 UTC) applied to the looser rule entry day <= 51. " "Test entries (days 50-51) have only ~22-24 days of confirmable follow-up: exits at ages between ~23 and 30 days are " "coded 0 (exit30 base rate is a lower bound) and the last-day markdowns may be unclassified. Reported as a sensitivity only.") metrics["notes"]["chronological"] = (f"primary chronological split: train entries before 2026-09-15 UTC, test 2026-09-15 to day " f"{int(S['entry_day'].max())} (2026-09-25); the plan's 2026-10-01 cutoff leaves no listing with a " "determinable 30-day outcome in the test fold") metrics["notes"]["targets"] = { "markdown30": "first sustained markdown (>=1 % and >=100 $) at listing age <= 30 d; exits before any cut count as 0", "markdown30cr": "competing-risk-aware: listings that exit within 30 d without a prior cut are removed; event = markdown within 30 d among the rest", "exit30": "disappearance confirmed by >= 7 d of absence at listing age <= 30 d (NOT a sale)", } metrics["notes"]["features"] = { "numeric": NUM_BASE, "dealer_target_encoding": DEALER_FEATS + ["(training fold only; inner 5-fold OOF for training rows; smoothing m=20)"], "categorical_lightgbm": CAT_ALL, "one_hot_logistic": CAT_COARSE, "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)", "excluded": "dealer_id (identity), anything measured after entry", } metrics["notes"]["models"] = {"logistic": "sklearn LogisticRegression L2, C=1, standardised numerics + missing indicators + one-hot coarse categoricals", "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)"}} 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|"} # ---- (vi) survival ---------------------------------------------------------------------------- log("[survival] Cox PH on S2 time-to-exit") metrics["survival"] = survival_cox(df) log(f" c-index train {metrics['survival']['c_index_train']:.4f} test {metrics['survival']['c_index_test']:.4f}") metrics["seconds_total"] = round(time.time() - t_start, 1) (OUT / "metrics.json").write_text(json.dumps(metrics, indent=2, default=float)) make_figure(cal_rows, metrics) write_table(metrics) log(f"done in {metrics['seconds_total']} s") return metrics # -------------------------------------------------------------------------------------- # figure + table # -------------------------------------------------------------------------------------- def make_figure(cal_rows: list[tuple[str, dict]], metrics: dict) -> None: import matplotlib.pyplot as plt from . import plotstyle as PS PS.setup() fig, (a, b) = plt.subplots(1, 2, figsize=(7.0, 3.1), gridspec_kw={"width_ratios": [1, 1.35]}) col = {"logistic": PS.CAT[0], "lightgbm": PS.CAT[1]} style = {"markdown30": ("-", "o", "Markdown within 30 d"), "exit30": ("--", "s", "Exit within 30 d")} a.plot([0, 1], [0, 1], color=PS.NEUTRAL, lw=0.8, zorder=1) for target, r in cal_rows: if target not in style: continue ls, mk, lab = style[target] for m in ("logistic", "lightgbm"): t = calibration_table(r["y"], r[f"p_{m}"]); t = t[t["n"] >= 20] a.plot(t["mean_pred"], t["frac_pos"], ls, marker=mk, color=col[m], lw=1.0, ms=3, label=f"{lab} — {'logistic' if m == 'logistic' else 'LightGBM'}") a.set_xlim(0, 1); a.set_ylim(0, 1); a.set_xlabel("Predicted probability"); a.set_ylabel("Observed frequency") a.set_title("(a) Reliability, chronological split"); a.legend(loc="upper left", fontsize=6.2) splits = [("chronological", "Chrono."), ("dealer_holdout", "Dealer"), ("region_holdout", "Region")] labels, xs = [], [] x = 0.0 for target, tlab in [("markdown30", "Markdown"), ("exit30", "Exit")]: for s, slab in splits: for j, m in enumerate(("logistic", "lightgbm")): mm = metrics["targets"][target][s]["models"][m] lo, hi = mm["log_loss_ci"] 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) 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) labels.append(f"{slab}\n{tlab}"); xs.append(x + 0.19); x += 1.2 x += 0.35 b.set_xticks(xs); b.set_xticklabels(labels, fontsize=6.5); b.set_ylabel("Log loss (test)") b.set_title("(b) Out-of-sample log loss by split (30-day targets)"); b.legend(loc="upper left", fontsize=7) b.grid(axis="x", visible=False) PS.save(fig, "fig10_prediction") def _fmt_int(x) -> str: return f"{int(round(float(x))):,}".replace(",", "\\,") def write_table(metrics: dict) -> None: rows = ["\\begin{tabular}{llllrrrr}", "\\toprule", "Target & Split & Model & Log loss [95\\% CI] & Brier & AUC & $n_{\\text{test}}$ & Base rate \\\\", "\\midrule"] tl = {"markdown30": "Markdown $\\le$ 30 d", "exit30": "Exit $\\le$ 30 d"} sl = {"chronological": "Chronological", "dealer_holdout": "Dealer holdout", "region_holdout": "Region LORO"} for target in ("markdown30", "exit30"): for s in ("chronological", "dealer_holdout", "region_holdout"): rec = metrics["targets"][target][s] for k, m in enumerate(("logistic", "lightgbm")): mm = rec["models"][m] lead = f"{tl[target]} & {sl[s]}" if k == 0 else " & " 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}] " f"& {mm['brier']:.3f} & {mm['auc']:.3f} & {_fmt_int(mm['n_test'])} & {mm['base_rate']:.3f} \\\\") if target == "markdown30": rows.append("\\addlinespace") rows += ["\\bottomrule", "\\end{tabular}"] C.TABLES.mkdir(parents=True, exist_ok=True) (C.TABLES / "tab08_prediction.tex").write_text("\n".join(rows) + "\n", encoding="utf-8")