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 RS = RS.with_columns(pl.col("first_cut").fill_null(0), pl.col("cut").fill_null(0), pl.col("exit").fill_null(0))40 df = RS.to_pandas()41 assert df["first_cut"].isna().sum() == 042 df["age_bin"] = pd.cut(df["listing_age"], bins=AGE_BINS, labels=AGE_LABELS, right=False)43 df["ldens"] = np.log1p(df["dens_25"]); df["ldens_entry"] = np.log1p(df["dens25_entry"])44 df["lkm"] = np.log(df["mileage_km"]); df["lprice"] = np.log(df["entry_price"])45 df["body_filled"] = df["body_filled"].fillna("unknown"); df["region"] = df["region"].fillna("(missing)")46 # price position is undefined when no competitor is active: set to 0 and flag (keeps the full risk set)47 for c in ["pricepos_25", "pricepos_entry"]:48 df[c] = df[c].replace([np.inf, -np.inf], np.nan)49 df["nocomp"] = df["pricepos_25"].isna().astype(float); df["nocomp_entry"] = df["pricepos_entry"].isna().astype(float)50 df["pricepos_25"] = df["pricepos_25"].fillna(0.0); df["pricepos_entry"] = df["pricepos_entry"].fillna(0.0)51 return df525354def design(df: pd.DataFrame, numeric: list[str], cats: list[str]) -> pd.DataFrame:55 X = df[numeric].astype(float).copy()56 for c in cats:57 d = pd.get_dummies(df[c].astype(str), prefix=c, drop_first=True, dtype=float)58 X = pd.concat([X, d], axis=1)59 X = sm.add_constant(X, has_constant="add")60 return X616263def cloglog(df: pd.DataFrame, y: str, numeric: list[str], cats: list[str], name: str) -> dict:64 t0 = time.time()65 d = df.dropna(subset=numeric + [y]).copy()66 X = design(d, numeric, cats)67 # drop all-zero / collinear dummies68 X = X.loc[:, (X != 0).any(axis=0)]69 model = sm.GLM(d[y].values.astype(float), X, family=sm.families.Binomial(link=sm.families.links.CLogLog()))70 res = model.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["dealer_id"])[0]}, maxiter=100)71 tab = pd.DataFrame({"term": X.columns, "coef": res.params.values, "se": res.bse.values, "p": res.pvalues.values,72 "ci_lo": res.conf_int().iloc[:, 0].values, "ci_hi": res.conf_int().iloc[:, 1].values})73 tab.to_csv(OUT / f"cloglog_{name}.csv", index=False)74 # joint tests75 def joint(prefix):76 idx = [i for i, c in enumerate(X.columns) if c.startswith(prefix)]77 if not idx:78 return None79 b = res.params.values[idx]; V = res.cov_params().values[np.ix_(idx, idx)]80 W = float(b @ np.linalg.solve(V, b)); from scipy.stats import chi281 return {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}82 out = {"n_rows": int(len(d)), "n_events": int(d[y].sum()), "n_listings": int(d["listing_id"].nunique()),83 "llf": float(res.llf), "joint_age_bins": joint("age_bin_"), "seconds": round(time.time() - t0, 1),84 "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_"))}}85 return out868788def empirical_hazard(df: pd.DataFrame, y: str, risk_col: str | None, name: str) -> pd.DataFrame:89 d = df if risk_col is None else df[df[risk_col]]90 g = d.groupby("listing_age").agg(events=(y, "sum"), at_risk=(y, "size")).reset_index()91 g = g[g["listing_age"] <= 70]92 p = g["events"] / g["at_risk"]; n = g["at_risk"]; z = norm.ppf(0.975)93 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)) / den94 g["hazard"] = p; g["lo"] = ctr - half; g["hi"] = ctr + half95 g.to_csv(OUT / f"empirical_hazard_{name}.csv", index=False)96 return g979899def listing_level(df: pd.DataFrame) -> pd.DataFrame:100 first = df.sort_values(["listing_id", "t"]).groupby("listing_id").first().reset_index()101 first["T_exit"] = first["duration_days"].clip(lower=0.01); first["E_exit"] = first["exit_confirmed"].astype(int)102 # time to first cut: from first_cut_day (day index) else censored at last day103 last = df.groupby("listing_id").agg(last_t=("t", "max"), fc=("first_cut_day", "first")).reset_index()104 first = first.merge(last, on="listing_id")105 first["E_cut"] = first["fc"].notna().astype(int)106 first["T_cut"] = np.where(first["E_cut"] == 1, first["fc"] - first["e_day"], first["last_t"] - first["e_day"]).astype(float) + 0.5107 return first108109110def cox_and_aft(first: pd.DataFrame, T: str, E: str, name: str) -> dict:111 from lifelines import CoxPHFitter, LogLogisticAFTFitter112 from lifelines.statistics import proportional_hazard_test113 cov = ["age_years_entry", "lkm", "lprice", "resid_A2", "ldens_entry", "pricepos_entry", "nocomp_entry"]114 d = first.dropna(subset=cov + [T, E]).copy()115 X = pd.concat([d[cov], pd.get_dummies(d["powertrain"], prefix="pt", drop_first=True, dtype=float),116 pd.get_dummies(d["body_filled"], prefix="body", drop_first=True, dtype=float),117 pd.get_dummies(d["source_class"], prefix="src", drop_first=True, dtype=float)], axis=1)118 X[T] = d[T].values; X[E] = d[E].values119 X = X.loc[:, (X != 0).any(axis=0)]120 X["dealer_c"] = pd.factorize(d["dealer_id"])[0]121 cph = CoxPHFitter(penalizer=0.001)122 cph.fit(X, duration_col=T, event_col=E, cluster_col="dealer_c", robust=True) # dealer-clustered sandwich SEs123 s = cph.summary.reset_index().rename(columns={"covariate": "term"})124 s.to_csv(OUT / f"cox_{name}.csv", index=False)125 Xs = X.drop(columns=["dealer_c"])126 sub = Xs.sample(min(30000, len(Xs)), random_state=36)127 cph2 = CoxPHFitter(penalizer=0.001).fit(sub, duration_col=T, event_col=E)128 ph = proportional_hazard_test(cph2, sub, time_transform="rank").summary.reset_index()129 ph.to_csv(OUT / f"cox_ph_test_{name}.csv", index=False)130 aft = LogLogisticAFTFitter(penalizer=0.001).fit(Xs, duration_col=T, event_col=E)131 a = aft.summary.reset_index(); a.to_csv(OUT / f"aft_loglogistic_{name}.csv", index=False)132 return {"n": int(len(X)), "events": int(X[E].sum()), "concordance": float(cph.concordance_index_),133 "ph_test_min_p": float(ph["p"].min()), "ph_violations_p<0.01": int((ph["p"] < 0.01).sum()),134 "cox_coef": {t: {"b": float(b), "p": float(p)} for t, b, p in zip(s["term"], s["coef"], s["p"])},135 "aft_ll": float(aft.log_likelihood_)}136137138def cut_size(df: pd.DataFrame) -> dict:139 import pyfixest as pf140 ev = df[df["first_cut"] == 1].copy()141 ev["lcut"] = np.log(-ev["first_cut_pct"].astype(float))142 ev["age_at_cut"] = ev["listing_age"].astype(float)143 ev = ev.dropna(subset=["lcut", "ldens", "resid_A2", "pricepos_25"])144 ev["mm"] = pd.factorize(ev["make"] + "|" + ev["model"])[0]145 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",146 data=ev, vcov={"CRV1": "dealer_id"}, fixef_rm="none")147 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); t.to_csv(OUT / "cutsize_ols.csv", index=False)148 return {"n": int(m._N), "r2_within": float(m._r2_within), "median_first_cut_pct": float(-ev["first_cut_pct"].median() * 100),149 "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(")}}150151152def run() -> dict:153 OUT.mkdir(parents=True, exist_ok=True)154 df = load()155 out = {}156 num_B = ["ldens", "resid_A2", "pricepos_entry", "nocomp_entry", "lprice", "lkm", "age_years_entry", "own_stock"]157 cats = ["age_bin", "powertrain", "body_filled", "source_class", "region", "cal_week"]158 print("B cloglog", flush=True)159 out["B_first_markdown_cloglog"] = cloglog(df[df["in_first_cut_risk"]], "first_cut", num_B, cats, "B_first_markdown")160 print("C cloglog", flush=True)161 num_C = ["ldens", "resid_A2", "pricepos_25", "nocomp", "prior_cuts", "lprice", "lkm", "age_years_entry", "own_stock"]162 out["C_exit_cloglog"] = cloglog(df, "exit", num_C, cats, "C_exit")163 empirical_hazard(df, "first_cut", "in_first_cut_risk", "first_markdown")164 empirical_hazard(df, "exit", None, "exit")165 # by source class, for heterogeneity166 for sc in df["source_class"].unique():167 empirical_hazard(df[df["source_class"] == sc], "exit", None, f"exit_{sc}")168 first = listing_level(df)169 print("Cox/AFT exit", flush=True)170 out["C_exit_cox"] = cox_and_aft(first, "T_exit", "E_exit", "exit")171 print("Cox/AFT cut", flush=True)172 out["B_cut_cox"] = cox_and_aft(first, "T_cut", "E_cut", "first_markdown")173 out["B_cut_size"] = cut_size(df)174 json.dump(out, open(OUT / "summary.json", "w"), indent=2, default=float)175 return out176177178if __name__ == "__main__":179 print(json.dumps(run(), indent=2, default=float))180