"""Models B (first markdown hazard, cut size) and C (listing-exit hazard) on the S2 listing-day risk set. Discrete-time complementary log-log hazards (statsmodels GLM, cluster-robust by dealer), Cox PH with Schoenfeld test and a log-logistic AFT at the listing level (lifelines), empirical hazard curves by listing age with Wilson intervals (for Figure 5). """ from __future__ import annotations import json import time import numpy as np import pandas as pd import polars as pl import statsmodels.api as sm from scipy.stats import norm from . import config as C OUT = C.RESULTS / "hazards" AGE_BINS = [0, 1, 4, 8, 15, 22, 31, 46, 61, 10_000] AGE_LABELS = ["0", "1-3", "4-7", "8-14", "15-21", "22-30", "31-45", "46-60", "61+"] def load(sample: str = "s2_dynamic") -> pd.DataFrame: RS = pl.read_parquet(C.PROCESSED / "risk_set.parquet") RS = RS.filter(pl.col(sample)) if sample == "s2_dynamic" else RS L = pl.read_parquet(C.PROCESSED / "listings.parquet").select( ["listing_id", "dealer_id", "source_class", "make", "model", "year", "entry_price", "mileage_km", "powertrain", "body_filled", "region", "age_years_entry", "exit_confirmed", "duration_days", "first_cut_pct", "n_cuts"]) R = pl.read_parquet(C.RESULTS / "hedonic" / "residuals.parquet").select(["listing_id", "resid_A2"]) RS = RS.join(L, on="listing_id", how="left").join(R, on="listing_id", how="left") # entry-day competition (static) and cumulative prior cuts (time-varying) RS = RS.sort(["listing_id", "t"]).with_columns([ pl.col("dens_25").first().over("listing_id").alias("dens25_entry"), pl.col("pricepos_25").first().over("listing_id").alias("pricepos_entry"), (pl.col("cut").cum_sum().over("listing_id") - pl.col("cut")).alias("prior_cuts"), (pl.col("t") // 7).alias("cal_week"), ]) RS = RS.with_columns(pl.col("first_cut").fill_null(0), pl.col("cut").fill_null(0), pl.col("exit").fill_null(0)) df = RS.to_pandas() assert df["first_cut"].isna().sum() == 0 df["age_bin"] = pd.cut(df["listing_age"], bins=AGE_BINS, labels=AGE_LABELS, right=False) df["ldens"] = np.log1p(df["dens_25"]); df["ldens_entry"] = np.log1p(df["dens25_entry"]) df["lkm"] = np.log(df["mileage_km"]); df["lprice"] = np.log(df["entry_price"]) df["body_filled"] = df["body_filled"].fillna("unknown"); df["region"] = df["region"].fillna("(missing)") # price position is undefined when no competitor is active: set to 0 and flag (keeps the full risk set) for c in ["pricepos_25", "pricepos_entry"]: df[c] = df[c].replace([np.inf, -np.inf], np.nan) df["nocomp"] = df["pricepos_25"].isna().astype(float); df["nocomp_entry"] = df["pricepos_entry"].isna().astype(float) df["pricepos_25"] = df["pricepos_25"].fillna(0.0); df["pricepos_entry"] = df["pricepos_entry"].fillna(0.0) return df def design(df: pd.DataFrame, numeric: list[str], cats: list[str]) -> pd.DataFrame: X = df[numeric].astype(float).copy() for c in cats: d = pd.get_dummies(df[c].astype(str), prefix=c, drop_first=True, dtype=float) X = pd.concat([X, d], axis=1) X = sm.add_constant(X, has_constant="add") return X def cloglog(df: pd.DataFrame, y: str, numeric: list[str], cats: list[str], name: str) -> dict: t0 = time.time() d = df.dropna(subset=numeric + [y]).copy() X = design(d, numeric, cats) # drop all-zero / collinear dummies X = X.loc[:, (X != 0).any(axis=0)] model = sm.GLM(d[y].values.astype(float), X, family=sm.families.Binomial(link=sm.families.links.CLogLog())) res = model.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["dealer_id"])[0]}, maxiter=100) tab = pd.DataFrame({"term": X.columns, "coef": res.params.values, "se": res.bse.values, "p": res.pvalues.values, "ci_lo": res.conf_int().iloc[:, 0].values, "ci_hi": res.conf_int().iloc[:, 1].values}) tab.to_csv(OUT / f"cloglog_{name}.csv", index=False) # joint tests def joint(prefix): idx = [i for i, c in enumerate(X.columns) if c.startswith(prefix)] if not idx: return None b = res.params.values[idx]; V = res.cov_params().values[np.ix_(idx, idx)] W = float(b @ np.linalg.solve(V, b)); from scipy.stats import chi2 return {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))} out = {"n_rows": int(len(d)), "n_events": int(d[y].sum()), "n_listings": int(d["listing_id"].nunique()), "llf": float(res.llf), "joint_age_bins": joint("age_bin_"), "seconds": round(time.time() - t0, 1), "coef": {c: {"b": float(res.params[c]), "se": float(res.bse[c]), "p": float(res.pvalues[c])} for c in X.columns if not c.startswith(("region_", "cal_week_", "age_bin_", "make_"))}} return out def empirical_hazard(df: pd.DataFrame, y: str, risk_col: str | None, name: str) -> pd.DataFrame: d = df if risk_col is None else df[df[risk_col]] g = d.groupby("listing_age").agg(events=(y, "sum"), at_risk=(y, "size")).reset_index() g = g[g["listing_age"] <= 70] p = g["events"] / g["at_risk"]; n = g["at_risk"]; z = norm.ppf(0.975) den = 1 + z ** 2 / n; ctr = (p + z ** 2 / (2 * n)) / den; half = z * np.sqrt(p * (1 - p) / n + z ** 2 / (4 * n ** 2)) / den g["hazard"] = p; g["lo"] = ctr - half; g["hi"] = ctr + half g.to_csv(OUT / f"empirical_hazard_{name}.csv", index=False) return g def listing_level(df: pd.DataFrame) -> pd.DataFrame: first = df.sort_values(["listing_id", "t"]).groupby("listing_id").first().reset_index() first["T_exit"] = first["duration_days"].clip(lower=0.01); first["E_exit"] = first["exit_confirmed"].astype(int) # time to first cut: from first_cut_day (day index) else censored at last day last = df.groupby("listing_id").agg(last_t=("t", "max"), fc=("first_cut_day", "first")).reset_index() first = first.merge(last, on="listing_id") first["E_cut"] = first["fc"].notna().astype(int) first["T_cut"] = np.where(first["E_cut"] == 1, first["fc"] - first["e_day"], first["last_t"] - first["e_day"]).astype(float) + 0.5 return first def cox_and_aft(first: pd.DataFrame, T: str, E: str, name: str) -> dict: from lifelines import CoxPHFitter, LogLogisticAFTFitter from lifelines.statistics import proportional_hazard_test cov = ["age_years_entry", "lkm", "lprice", "resid_A2", "ldens_entry", "pricepos_entry", "nocomp_entry"] d = first.dropna(subset=cov + [T, E]).copy() X = pd.concat([d[cov], pd.get_dummies(d["powertrain"], prefix="pt", drop_first=True, dtype=float), pd.get_dummies(d["body_filled"], prefix="body", drop_first=True, dtype=float), pd.get_dummies(d["source_class"], prefix="src", drop_first=True, dtype=float)], axis=1) X[T] = d[T].values; X[E] = d[E].values X = X.loc[:, (X != 0).any(axis=0)] X["dealer_c"] = pd.factorize(d["dealer_id"])[0] cph = CoxPHFitter(penalizer=0.001) cph.fit(X, duration_col=T, event_col=E, cluster_col="dealer_c", robust=True) # dealer-clustered sandwich SEs s = cph.summary.reset_index().rename(columns={"covariate": "term"}) s.to_csv(OUT / f"cox_{name}.csv", index=False) Xs = X.drop(columns=["dealer_c"]) sub = Xs.sample(min(30000, len(Xs)), random_state=36) cph2 = CoxPHFitter(penalizer=0.001).fit(sub, duration_col=T, event_col=E) ph = proportional_hazard_test(cph2, sub, time_transform="rank").summary.reset_index() ph.to_csv(OUT / f"cox_ph_test_{name}.csv", index=False) aft = LogLogisticAFTFitter(penalizer=0.001).fit(Xs, duration_col=T, event_col=E) a = aft.summary.reset_index(); a.to_csv(OUT / f"aft_loglogistic_{name}.csv", index=False) return {"n": int(len(X)), "events": int(X[E].sum()), "concordance": float(cph.concordance_index_), "ph_test_min_p": float(ph["p"].min()), "ph_violations_p<0.01": int((ph["p"] < 0.01).sum()), "cox_coef": {t: {"b": float(b), "p": float(p)} for t, b, p in zip(s["term"], s["coef"], s["p"])}, "aft_ll": float(aft.log_likelihood_)} def cut_size(df: pd.DataFrame) -> dict: import pyfixest as pf ev = df[df["first_cut"] == 1].copy() ev["lcut"] = np.log(-ev["first_cut_pct"].astype(float)) ev["age_at_cut"] = ev["listing_age"].astype(float) ev = ev.dropna(subset=["lcut", "ldens", "resid_A2", "pricepos_25"]) ev["mm"] = pd.factorize(ev["make"] + "|" + ev["model"])[0] m = pf.feols("lcut ~ age_at_cut + I(age_at_cut**2) + ldens + resid_A2 + lprice + lkm + age_years_entry + C(powertrain) + C(body_filled) + C(source_class) | mm + region + cal_week", data=ev, vcov={"CRV1": "dealer_id"}, fixef_rm="none") t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); t.to_csv(OUT / "cutsize_ols.csv", index=False) return {"n": int(m._N), "r2_within": float(m._r2_within), "median_first_cut_pct": float(-ev["first_cut_pct"].median() * 100), "coef": {r.term: {"b": float(r.Estimate), "se": float(r["Std. Error"]), "p": float(r["Pr(>|t|)"])} for _, r in t.iterrows() if not r.term.startswith("C(")}} def run() -> dict: OUT.mkdir(parents=True, exist_ok=True) df = load() out = {} num_B = ["ldens", "resid_A2", "pricepos_entry", "nocomp_entry", "lprice", "lkm", "age_years_entry", "own_stock"] cats = ["age_bin", "powertrain", "body_filled", "source_class", "region", "cal_week"] print("B cloglog", flush=True) out["B_first_markdown_cloglog"] = cloglog(df[df["in_first_cut_risk"]], "first_cut", num_B, cats, "B_first_markdown") print("C cloglog", flush=True) num_C = ["ldens", "resid_A2", "pricepos_25", "nocomp", "prior_cuts", "lprice", "lkm", "age_years_entry", "own_stock"] out["C_exit_cloglog"] = cloglog(df, "exit", num_C, cats, "C_exit") empirical_hazard(df, "first_cut", "in_first_cut_risk", "first_markdown") empirical_hazard(df, "exit", None, "exit") # by source class, for heterogeneity for sc in df["source_class"].unique(): empirical_hazard(df[df["source_class"] == sc], "exit", None, f"exit_{sc}") first = listing_level(df) print("Cox/AFT exit", flush=True) out["C_exit_cox"] = cox_and_aft(first, "T_exit", "E_exit", "exit") print("Cox/AFT cut", flush=True) out["B_cut_cox"] = cox_and_aft(first, "T_cut", "E_cut", "first_markdown") out["B_cut_size"] = cut_size(df) json.dump(out, open(OUT / "summary.json", "w"), indent=2, default=float) return out if __name__ == "__main__": print(json.dumps(run(), indent=2, default=float))