1"""Models B (first markdown hazard, cut size) and C (listing-exit hazard) on the S2 listing-day risk set.23Discrete-time complementary log-log hazards (statsmodels GLM, cluster-robust by dealer),4Cox PH with Schoenfeld test and a log-logistic AFT at the listing level (lifelines),5empirical hazard curves by listing age with Wilson intervals (for Figure 5).6"""7from __future__ import annotations89import json10import time1112import numpy as np13import pandas as pd14import polars as pl15import statsmodels.api as sm16from scipy.stats import norm1718from . import config as C1920OUT = C.RESULTS / "hazards"21AGE_BINS = [0, 1, 4, 8, 15, 22, 31, 46, 61, 10_000]22AGE_LABELS = ["0", "1-3", "4-7", "8-14", "15-21", "22-30", "31-45", "46-60", "61+"]232425def load(sample: str = "s2_dynamic") -> pd.DataFrame:26 RS = pl.read_parquet(C.PROCESSED / "risk_set.parquet")27 RS = RS.filter(pl.col(sample)) if sample == "s2_dynamic" else RS28 L = pl.read_parquet(C.PROCESSED / "listings.parquet").select(29 ["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"])30 R = pl.read_parquet(C.RESULTS / "hedonic" / "residuals.parquet").select(["listing_id", "resid_A2"])31 RS = RS.join(L, on="listing_id", how="left").join(R, on="listing_id", how="left")32 # entry-day competition (static) and cumulative prior cuts (time-varying)33 RS = RS.sort(["listing_id", "t"]).with_columns([34 pl.col("dens_25").first().over("listing_id").alias("dens25_entry"),35 pl.col("pricepos_25").first().over("listing_id").alias("pricepos_entry"),36 (pl.col("cut").cum_sum().over("listing_id") - pl.col("cut")).alias("prior_cuts"),37 (pl.col("t") // 7).alias("cal_week"),38 ])39 df = RS.to_pandas()40 df["age_bin"] = pd.cut(df["listing_age"], bins=AGE_BINS, labels=AGE_LABELS, right=False)41 df["ldens"] = np.log1p(df["dens_25"]); df["ldens_entry"] = np.log1p(df["dens25_entry"])42 df["lkm"] = np.log(df["mileage_km"]); df["lprice"] = np.log(df["entry_price"])43 df["body_filled"] = df["body_filled"].fillna("unknown"); df["region"] = df["region"].fillna("(missing)")44 # price position is undefined when no competitor is active: set to 0 and flag (keeps the full risk set)45 for c in ["pricepos_25", "pricepos_entry"]:46 df[c] = df[c].replace([np.inf, -np.inf], np.nan)47 df["nocomp"] = df["pricepos_25"].isna().astype(float); df["nocomp_entry"] = df["pricepos_entry"].isna().astype(float)48 df["pricepos_25"] = df["pricepos_25"].fillna(0.0); df["pricepos_entry"] = df["pricepos_entry"].fillna(0.0)49 return df505152def design(df: pd.DataFrame, numeric: list[str], cats: list[str]) -> pd.DataFrame:53 X = df[numeric].astype(float).copy()54 for c in cats:55 d = pd.get_dummies(df[c].astype(str), prefix=c, drop_first=True, dtype=float)56 X = pd.concat([X, d], axis=1)57 X = sm.add_constant(X, has_constant="add")58 return X596061def cloglog(df: pd.DataFrame, y: str, numeric: list[str], cats: list[str], name: str) -> dict:62 t0 = time.time()63 d = df.dropna(subset=numeric + [y]).copy()64 X = design(d, numeric, cats)65 # drop all-zero / collinear dummies66 X = X.loc[:, (X != 0).any(axis=0)]67 model = sm.GLM(d[y].values.astype(float), X, family=sm.families.Binomial(link=sm.families.links.CLogLog()))68 res = model.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["dealer_id"])[0]}, maxiter=100)69 tab = pd.DataFrame({"term": X.columns, "coef": res.params.values, "se": res.bse.values, "p": res.pvalues.values,70 "ci_lo": res.conf_int().iloc[:, 0].values, "ci_hi": res.conf_int().iloc[:, 1].values})71 tab.to_csv(OUT / f"cloglog_{name}.csv", index=False)72 # joint tests73 def joint(prefix):74 idx = [i for i, c in enumerate(X.columns) if c.startswith(prefix)]75 if not idx:76 return None77 b = res.params.values[idx]; V = res.cov_params().values[np.ix_(idx, idx)]78 W = float(b @ np.linalg.solve(V, b)); from scipy.stats import chi279 return {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}80 out = {"n_rows": int(len(d)), "n_events": int(d[y].sum()), "n_listings": int(d["listing_id"].nunique()),81 "llf": float(res.llf), "joint_age_bins": joint("age_bin_"), "seconds": round(time.time() - t0, 1),82 "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_"))}}83 return out848586def empirical_hazard(df: pd.DataFrame, y: str, risk_col: str | None, name: str) -> pd.DataFrame:87 d = df if risk_col is None else df[df[risk_col]]88 g = d.groupby("listing_age").agg(events=(y, "sum"), at_risk=(y, "size")).reset_index()89 g = g[g["listing_age"] <= 70]90 p = g["events"] / g["at_risk"]; n = g["at_risk"]; z = norm.ppf(0.975)91 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)) / den92 g["hazard"] = p; g["lo"] = ctr - half; g["hi"] = ctr + half93 g.to_csv(OUT / f"empirical_hazard_{name}.csv", index=False)94 return g959697def listing_level(df: pd.DataFrame) -> pd.DataFrame:98 first = df.sort_values(["listing_id", "t"]).groupby("listing_id").first().reset_index()99 first["T_exit"] = first["duration_days"].clip(lower=0.01); first["E_exit"] = first["exit_confirmed"].astype(int)100 # time to first cut: from first_cut_day (day index) else censored at last day101 last = df.groupby("listing_id").agg(last_t=("t", "max"), fc=("first_cut_day", "first")).reset_index()102 first = first.merge(last, on="listing_id")103 first["E_cut"] = first["fc"].notna().astype(int)104 first["T_cut"] = np.where(first["E_cut"] == 1, first["fc"] - first["e_day"], first["last_t"] - first["e_day"]).astype(float) + 0.5105 return first106107108def cox_and_aft(first: pd.DataFrame, T: str, E: str, name: str) -> dict:109 from lifelines import CoxPHFitter, LogLogisticAFTFitter110 from lifelines.statistics import proportional_hazard_test111 cov = ["age_years_entry", "lkm", "lprice", "resid_A2", "ldens_entry", "pricepos_entry", "nocomp_entry"]112 d = first.dropna(subset=cov + [T, E]).copy()113 X = pd.concat([d[cov], pd.get_dummies(d["powertrain"], prefix="pt", drop_first=True, dtype=float),114 pd.get_dummies(d["body_filled"], prefix="body", drop_first=True, dtype=float),115 pd.get_dummies(d["source_class"], prefix="src", drop_first=True, dtype=float)], axis=1)116 X[T] = d[T].values; X[E] = d[E].values117 X = X.loc[:, (X != 0).any(axis=0)]118 X["dealer_c"] = pd.factorize(d["dealer_id"])[0]119 cph = CoxPHFitter(penalizer=0.001)120 cph.fit(X, duration_col=T, event_col=E, cluster_col="dealer_c", robust=True) # dealer-clustered sandwich SEs121 s = cph.summary.reset_index().rename(columns={"covariate": "term"})122 s.to_csv(OUT / f"cox_{name}.csv", index=False)123 Xs = X.drop(columns=["dealer_c"])124 sub = Xs.sample(min(30000, len(Xs)), random_state=36)125 cph2 = CoxPHFitter(penalizer=0.001).fit(sub, duration_col=T, event_col=E)126 ph = proportional_hazard_test(cph2, sub, time_transform="rank").summary.reset_index()127 ph.to_csv(OUT / f"cox_ph_test_{name}.csv", index=False)128 aft = LogLogisticAFTFitter(penalizer=0.001).fit(Xs, duration_col=T, event_col=E)129 a = aft.summary.reset_index(); a.to_csv(OUT / f"aft_loglogistic_{name}.csv", index=False)130 return {"n": int(len(X)), "events": int(X[E].sum()), "concordance": float(cph.concordance_index_),131 "ph_test_min_p": float(ph["p"].min()), "ph_violations_p<0.01": int((ph["p"] < 0.01).sum()),132 "cox_coef": {t: {"b": float(b), "p": float(p)} for t, b, p in zip(s["term"], s["coef"], s["p"])},133 "aft_ll": float(aft.log_likelihood_)}134135136def cut_size(df: pd.DataFrame) -> dict:137 import pyfixest as pf138 ev = df[df["first_cut"] == 1].copy()139 ev["lcut"] = np.log(-ev["first_cut_pct"].astype(float))140 ev["age_at_cut"] = ev["listing_age"].astype(float)141 ev = ev.dropna(subset=["lcut", "ldens", "resid_A2", "pricepos_25"])142 ev["mm"] = pd.factorize(ev["make"] + "|" + ev["model"])[0]143 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",144 data=ev, vcov={"CRV1": "dealer_id"}, fixef_rm="none")145 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); t.to_csv(OUT / "cutsize_ols.csv", index=False)146 return {"n": int(m._N), "r2_within": float(m._r2_within), "median_first_cut_pct": float(-ev["first_cut_pct"].median() * 100),147 "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(")}}148149150def run() -> dict:151 OUT.mkdir(parents=True, exist_ok=True)152 df = load()153 out = {}154 num_B = ["ldens", "resid_A2", "pricepos_entry", "nocomp_entry", "lprice", "lkm", "age_years_entry", "own_stock"]155 cats = ["age_bin", "powertrain", "body_filled", "source_class", "region", "cal_week"]156 print("B cloglog", flush=True)157 out["B_first_markdown_cloglog"] = cloglog(df[df["in_first_cut_risk"]], "first_cut", num_B, cats, "B_first_markdown")158 print("C cloglog", flush=True)159 num_C = ["ldens", "resid_A2", "pricepos_25", "nocomp", "prior_cuts", "lprice", "lkm", "age_years_entry", "own_stock"]160 out["C_exit_cloglog"] = cloglog(df, "exit", num_C, cats, "C_exit")161 empirical_hazard(df, "first_cut", "in_first_cut_risk", "first_markdown")162 empirical_hazard(df, "exit", None, "exit")163 # by source class, for heterogeneity164 for sc in df["source_class"].unique():165 empirical_hazard(df[df["source_class"] == sc], "exit", None, f"exit_{sc}")166 first = listing_level(df)167 print("Cox/AFT exit", flush=True)168 out["C_exit_cox"] = cox_and_aft(first, "T_exit", "E_exit", "exit")169 print("Cox/AFT cut", flush=True)170 out["B_cut_cox"] = cox_and_aft(first, "T_cut", "E_cut", "first_markdown")171 out["B_cut_size"] = cut_size(df)172 json.dump(out, open(OUT / "summary.json", "w"), indent=2, default=float)173 return out174175176if __name__ == "__main__":177 print(json.dumps(run(), indent=2, default=float))178