1"""Phase 6 — robustness suite (analysis_plan §6; identification threats and sensitivity).23Every variant is re-derived in memory from data/interim (price_log, vehicles,4listings_full) and data/processed (listings, risk_set); no existing module or5processed file is modified. Blocks:67 1. Hedonic H1/H4 (Model A2: splines + make-model FE, dealer-clustered):8 sample variants (sampled portals, Kijiji, 1 000-500 000 $ price filter),9 region joint Wald under each, heterogeneity of the linear age/km slopes,10 composition reweighting to the otogo make-model-year mix (pyfixest weights).11 2. Hazards H2/H5 (discrete-time cloglog, GLM, dealer-clustered): markdown12 thresholds 2 %/200 $ and 0.5 %/50 $, transient rule 24 h / 72 h, exit13 confirmation 3 / 14 d, S2L vs S2. Listing-days are subsampled to <= MAX_ROWS14 with a fixed seed (reported).15 3. AMLI H6: rank standardisation, MIN_LISTING_DAYS = 100, 2-week aggregation,16 province vs region, leave-one-component-out predictive validation.17 4. Holm correction over the six primary preregistered tests.18 5. paper/tables/tab09_robustness.tex (booktabs fragment).19"""20from __future__ import annotations2122import json23import time24import warnings2526import numpy as np27import pandas as pd28import polars as pl29import pyfixest as pf30import statsmodels.api as sm31from scipy.stats import chi2, norm3233from . import amli, config as C, features, hazards, hedonic, pricepaths3435OUT = C.RESULTS / "robustness"36warnings.filterwarnings("ignore", message=".*multicollinearity.*") # na6/nk6 (redundant orthogonalised spline columns) as in hedonic.py37SEED = 3638MAX_ROWS = 800_00039HAZ_TERMS = ["ldens", "resid_A2", "lprice", "lkm"]40NUM_B = ["ldens", "resid_A2", "pricepos_entry", "nocomp_entry", "lprice", "lkm", "age_years_entry", "own_stock"]41NUM_C = ["ldens", "resid_A2", "pricepos_25", "nocomp", "prior_cuts", "lprice", "lkm", "age_years_entry", "own_stock"]42CATS = ["age_bin", "powertrain", "body_filled", "source_class", "region", "cal_week"]43_T0 = time.time()444546def _log(msg: str) -> None:47 print(f"[{time.time() - _T0:7.1f}s] {msg}", flush=True)484950# =============================================================================== 1. hedonic5152def _wald(m, cols: list[str]) -> dict:53 names = list(m._coefnames)54 idx = [names.index(c) for c in cols if c in names]55 if not idx:56 return {"stat": None, "df": 0, "p": None}57 b = np.asarray(m.coef())[idx]; V = np.asarray(m._vcov)[np.ix_(idx, idx)]58 W = float(b @ np.linalg.solve(V, b))59 return {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}606162def _term(m, term: str) -> dict:63 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"})64 r = t[t["term"] == term]65 if r.empty:66 return {"b": None, "se": None, "lo": None, "hi": None, "p": None}67 r = r.iloc[0]68 return {"b": float(r["Estimate"]), "se": float(r["Std. Error"]), "lo": float(r["2.5%"]), "hi": float(r["97.5%"]), "p": float(r["Pr(>|t|)"])}697071def _region_premia(m) -> pd.DataFrame:72 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"})73 rp = t[t["term"].str.startswith("C(region)")].copy()74 rp["region"] = rp["term"].str.extract(r"\[T\.(.*)\]")75 return rp[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].rename(columns={"Estimate": "est", "Std. Error": "se", "2.5%": "lo", "97.5%": "hi"})767778def hedonic_frame(flt: pl.Expr) -> pd.DataFrame:79 """Same preparation as hedonic.load() but with an arbitrary sample filter."""80 L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(flt)81 df = L.select(["listing_id", "dealer_id", "source", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km",82 "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "age_years_entry"]).to_pandas()83 df = df.rename(columns={"body_filled": "body"})84 df["body"] = df["body"].fillna("unknown"); df["trim"] = df["trim"].fillna("(none)"); df["region"] = df["region"].fillna("(missing)")85 df["logp"] = np.log(df["entry_price"]); df["age"] = df["age_years_entry"]; df["km"] = df["mileage_km"] / 1e486 df["week"] = ((df["entry_ts"] - C.WINDOW_START_TS) // (7 * 86400)).astype(int)87 df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60)88 mm_s = df["make"] + "|" + df["model"]; mmt_s = mm_s + "|" + df["trim"]89 df["mm"] = pd.factorize(mm_s)[0]; df["mm_label"] = mm_s; df["mmt"] = pd.factorize(mmt_s)[0]90 df["mmy_label"] = mm_s + "|" + df["year"].astype(str)91 df["region"] = pd.Categorical(df["region"], categories=["Montréal"] + sorted(r for r in df["region"].unique() if r != "Montréal"))92 for col, ref in (("powertrain", "ICE"), ("transmission", "Automatic"), ("drivetrain", "AWD/4x4"), ("source_class", "direct_dealer")):93 levels = ([ref] if ref in set(df[col]) else []) + sorted(x for x in df[col].unique() if x != ref)94 df[col] = pd.Categorical(df[col], categories=levels)95 return df.dropna(subset=["age", "km", "logp"])969798def _rhs(acols, kcols, controls) -> str:99 return " + ".join(["age", "km"] + acols + kcols + [f"C({c})" for c in controls])100101102def fit_A2(df: pd.DataFrame, controls=hedonic.CONTROLS, weights: str | None = None, region: bool = True, fe: str = "mm") -> tuple[dict, pd.DataFrame]:103 """A2 (fe='mm': make-model FE) or A4 (fe='mmt': make-model-trim FE) with the same spline RHS."""104 X, acols, kcols, _ = hedonic.spline_basis(df)105 fml = f"logp ~ {_rhs(acols, kcols, controls)}" + (" + C(region)" if region else "") + f" | {fe} + week"106 kw = {"weights": weights} if weights else {}107 m = pf.feols(fml, data=X, vcov={"CRV1": "dealer_id"}, fixef_rm="none", **kw)108 res = {"n": int(m._N), "n_dealers": int(X["dealer_id"].nunique()), "r2": float(m._r2), "rmse": float(np.sqrt(np.mean(m.resid() ** 2))),109 "age": _term(m, "age"), "km": _term(m, "km"),110 "H1_nonlinear_joint": _wald(m, acols + kcols), "H1_nonlinear_age": _wald(m, acols), "H1_nonlinear_km": _wald(m, kcols)}111 names = list(m._coefnames)112 res["H4_region_joint"] = _wald(m, [n for n in names if n.startswith("C(region)")]) if region else None113 return res, (_region_premia(m) if region else pd.DataFrame())114115116def otogo_weights(df: pd.DataFrame, cap: float = 10.0) -> tuple[pd.Series, dict]:117 """Post-stratification weights: direct-dealer + banner listings reweighted to the otogo118 make-model-year distribution (fallback to the make-model cell; cells absent from otogo119 get weight 0 and are dropped); other sources keep weight 1."""120 tgt = df[df["source"] == "otogo"]; smp = df[df["source_class"].isin(["direct_dealer", "banner"])]121 t_mmy = tgt["mmy_label"].value_counts(normalize=True); s_mmy = smp["mmy_label"].value_counts(normalize=True)122 t_mm = tgt["mm_label"].value_counts(normalize=True); s_mm = smp["mm_label"].value_counts(normalize=True)123 w = pd.Series(1.0, index=df.index)124 r_mmy = (t_mmy / s_mmy).reindex(smp["mmy_label"]).values125 r_mm = (t_mm / s_mm).reindex(smp["mm_label"]).values126 ws = np.where(np.isfinite(r_mmy), r_mmy, np.where(np.isfinite(r_mm), r_mm, 0.0))127 n_fallback_mm = int((~np.isfinite(r_mmy) & np.isfinite(r_mm)).sum()); n_zero = int((ws == 0).sum()); n_capped = int((ws > cap).sum())128 ws = np.minimum(ws, cap)129 pos = ws > 0130 ws[pos] = ws[pos] / ws[pos].mean() # mean 1 among retained dealer/banner listings131 w.loc[smp.index] = ws132 info = {"n_target_otogo": int(len(tgt)), "n_reweighted": int(len(smp)), "n_fallback_mm_cell": n_fallback_mm, "n_zero_weight_dropped": n_zero,133 "n_capped_at_10": n_capped, "weight_max": float(ws.max()), "weight_cv": float(ws[pos].std() / ws[pos].mean()), "cap": cap}134 return w, info135136137def run_hedonic() -> dict:138 OUT.mkdir(parents=True, exist_ok=True)139 out, rows, prem = {}, [], []140 base_flt = pl.col("s1_hedonic")141 s1core = (pl.col("f_auto") & pl.col("f_canonical") & pl.col("f_year") & pl.col("f_km") & pl.col("f_mm") & pl.col("f_nocollision"))142 variants = {143 "baseline_S1": base_flt,144 "plus_sampled_portals": pl.col("s1_all_sources") & (pl.col("source_class") != "private_marketplace"),145 "plus_kijiji": pl.col("s1_all_sources") & (pl.col("source_class") != "sampled_portal"),146 "plus_sampled_and_kijiji": pl.col("s1_all_sources"),147 "price_1000_500000": s1core & pl.col("entry_price").is_between(1000.0, 500000.0) & pl.col("f_not_sampled_or_private"),148 }149 base_prem = None150 for name, flt in variants.items():151 df = hedonic_frame(flt)152 res, rp = fit_A2(df)153 r4, rp4 = fit_A2(df, fe="mmt"); res["A4_mmt"] = {k: r4[k] for k in ("n", "r2", "age", "km", "H4_region_joint")}; res["H4_region_joint_A4"] = r4["H4_region_joint"]154 res["source_classes"] = sorted(df["source_class"].unique().tolist())155 if base_prem is None:156 base_prem = rp.set_index("region")["est"]157 res["region_premia_corr_with_baseline"] = float(rp.set_index("region")["est"].reindex(base_prem.index).corr(base_prem))158 res["region_premia_max_abs"] = float(rp["est"].abs().max())159 out[name] = res160 rows.append({"variant": name, **_flat_hed(res)}); rp.insert(0, "variant", name); prem.append(rp)161 _log(f"hedonic {name}: n={res['n']} age={res['age']['b']:.4f} km={res['km']['b']:.4f} H1p={res['H1_nonlinear_joint']['p']:.2e} H4p(A2)={res['H4_region_joint']['p']:.3f} H4p(A4)={res['H4_region_joint_A4']['p']:.3f}")162 # --- (f) composition reweighting --------------------------------------------163 df = hedonic_frame(base_flt)164 w, winfo = otogo_weights(df)165 df["w_otogo"] = w.values166 dfw = df[df["w_otogo"] > 0].copy()167 res_w, rp_w = fit_A2(dfw, weights="w_otogo")168 r4w, _ = fit_A2(dfw, weights="w_otogo", fe="mmt"); res_w["A4_mmt"] = {k: r4w[k] for k in ("n", "r2", "age", "km", "H4_region_joint")}; res_w["H4_region_joint_A4"] = r4w["H4_region_joint"]169 res_w["weights"] = winfo170 res_w["region_premia_corr_with_baseline"] = float(rp_w.set_index("region")["est"].reindex(base_prem.index).corr(base_prem)); res_w["region_premia_max_abs"] = float(rp_w["est"].abs().max())171 out["reweighted_to_otogo_mmy"] = res_w; rows.append({"variant": "reweighted_to_otogo_mmy", **_flat_hed(res_w)}); rp_w.insert(0, "variant", "reweighted_to_otogo_mmy"); prem.append(rp_w)172 _log(f"hedonic reweighted: n={res_w['n']} age={res_w['age']['b']:.4f} H4p={res_w['H4_region_joint']['p']:.3f} {winfo}")173 # direct+banner only: unweighted vs reweighted174 dd = df[df["source_class"].isin(["direct_dealer", "banner"])].copy()175 res_dd, rp_dd = fit_A2(dd); out["direct_banner_only_unweighted"] = res_dd; rows.append({"variant": "direct_banner_only_unweighted", **_flat_hed(res_dd)})176 ddw = dd[dd["w_otogo"] > 0].copy(); res_ddw, rp_ddw = fit_A2(ddw, weights="w_otogo"); out["direct_banner_only_reweighted"] = res_ddw; rows.append({"variant": "direct_banner_only_reweighted", **_flat_hed(res_ddw)})177 for nm, rp_ in (("direct_banner_only_unweighted", rp_dd), ("direct_banner_only_reweighted", rp_ddw)):178 rp_.insert(0, "variant", nm); prem.append(rp_)179 out[nm]["region_premia_corr_with_baseline"] = float(rp_.set_index("region")["est"].reindex(base_prem.index).corr(base_prem))180 _log(f"hedonic direct+banner: unweighted age={res_dd['age']['b']:.4f} H4p={res_dd['H4_region_joint']['p']:.3f}; reweighted age={res_ddw['age']['b']:.4f} H4p={res_ddw['H4_region_joint']['p']:.3f}")181 pd.DataFrame(rows).to_csv(OUT / "hedonic_variants.csv", index=False)182 pd.concat(prem).to_csv(OUT / "hedonic_region_premia_variants.csv", index=False)183 # --- (e) heterogeneity -------------------------------------------------------184 het_rows, het = [], {}185 df["yband"] = pd.cut(df["year"], bins=[0, 2015, 2020, 9999], labels=["<=2015", "2016-2020", "2021+"]).astype(str)186 groups = {"powertrain": ["ICE", "BEV", "HEV"], "yband": ["<=2015", "2016-2020", "2021+"], "body": ["SUV", "Sedan", "Pickup"],187 "source_class": ["direct_dealer", "full_portal", "banner"]}188 for gvar, levels in groups.items():189 controls = [c for c in hedonic.CONTROLS if c != gvar]190 for lv in levels:191 sub = df[df[gvar] == lv]192 if len(sub) < 1000:193 continue194 r, _ = fit_A2(sub, controls=controls)195 het[f"{gvar}={lv}"] = {k: r[k] for k in ("n", "n_dealers", "r2", "age", "km", "H1_nonlinear_joint", "H4_region_joint")}196 het_rows.append({"group_var": gvar, "group": lv, "n": r["n"], "age_b": r["age"]["b"], "age_se": r["age"]["se"], "age_lo": r["age"]["lo"], "age_hi": r["age"]["hi"],197 "km_b": r["km"]["b"], "km_se": r["km"]["se"], "km_lo": r["km"]["lo"], "km_hi": r["km"]["hi"], "H1_joint_p": r["H1_nonlinear_joint"]["p"], "H4_region_p": r["H4_region_joint"]["p"]})198 _log(f"hedonic het {gvar}={lv}: n={r['n']} age={r['age']['b']:.4f} ({r['age']['se']:.4f}) km={r['km']['b']:.4f} ({r['km']['se']:.4f})")199 out["heterogeneity"] = het200 pd.DataFrame(het_rows).to_csv(OUT / "hedonic_heterogeneity.csv", index=False)201 pub = json.load(open(C.RESULTS / "hedonic" / "summary.json"))202 out["published"] = {"H1_nonlinear_joint_A2": pub["tests"]["H1_nonlinear_joint"], "H4_region_joint_A4": pub["tests"]["H4_region_joint"], "n_S1": pub["n_S1"]}203 return out204205206def _flat_hed(r: dict) -> dict:207 return {"n": r["n"], "n_dealers": r["n_dealers"], "r2": r["r2"], "rmse": r["rmse"], "age_b": r["age"]["b"], "age_se": r["age"]["se"], "age_lo": r["age"]["lo"], "age_hi": r["age"]["hi"],208 "km_b": r["km"]["b"], "km_se": r["km"]["se"], "km_lo": r["km"]["lo"], "km_hi": r["km"]["hi"],209 "H1_joint_stat": r["H1_nonlinear_joint"]["stat"], "H1_joint_p": r["H1_nonlinear_joint"]["p"],210 "H4_region_stat": r["H4_region_joint"]["stat"], "H4_region_df": r["H4_region_joint"]["df"], "H4_region_p": r["H4_region_joint"]["p"],211 "H4_region_p_A4_mmt": (r.get("H4_region_joint_A4") or {}).get("p"), "age_b_A4": ((r.get("A4_mmt") or {}).get("age") or {}).get("b"),212 "region_premia_max_abs": r.get("region_premia_max_abs"), "region_premia_corr_with_baseline": r.get("region_premia_corr_with_baseline")}213214215# =============================================================================== 2. hazards216217def rebuild_events(transient_seconds: float, min_pct: float, min_abs: float) -> tuple[pd.DataFrame, pd.DataFrame]:218 """Re-run the sustained-path algorithm under alternative rules for the S2L universe.219 Returns (per-listing: first_cut_day, n_cuts, oscillating, entry_price) and (listing_id, t) cut days."""220 LF = pl.read_parquet(C.INTERIM / "listings_full.parquet").filter(pl.col("s2l_dynamic_with_stock")).select(["uid", "listing_id"])221 V = pl.read_parquet(C.INTERIM / "vehicles.parquet").select(["uid", "last_seen"])222 P = pl.read_parquet(C.INTERIM / "price_log.parquet").filter(pl.col("uid").is_in(LF["uid"].to_list()))223 ls = dict(zip(V["uid"].to_list(), V["last_seen"].to_list()))224 lid = dict(zip(LF["uid"].to_list(), LF["listing_id"].to_list()))225 grouped = P.sort(["uid", "ts"]).group_by("uid", maintain_order=True).agg([pl.col("ts"), pl.col("price")])226 rows, cut_rows = [], []227 for uid, ts, pr in grouped.iter_rows():228 ps = pricepaths.sustain(ts, pr, ls.get(uid, max(ts)), transient_seconds=transient_seconds)229 ev = pricepaths.classify(ps.levels, min_pct=min_pct, min_abs=min_abs)230 cuts = [e for e in ev if e["kind"] == "cut"]231 l = lid[uid]232 for e in cuts:233 cut_rows.append((l, int(features.day_of(np.array([e["ts"]]))[0])))234 rows.append({"listing_id": l, "first_cut_day_v": int(features.day_of(np.array([cuts[0]["ts"]]))[0]) if cuts else np.nan,235 "n_cuts_v": len(cuts), "oscillating_v": ps.n_transients >= C.OSCILLATION_MAX_TRANSIENTS,236 "entry_price_v": ev[0]["price"] if ev else np.nan, "first_cut_pct_v": cuts[0]["change_pct"] if cuts else np.nan})237 per = pd.DataFrame(rows)238 cd = pd.DataFrame(cut_rows, columns=["listing_id", "t"]).drop_duplicates(); cd["cut_v"] = 1.0239 return per, cd240241242def _cloglog(d: pd.DataFrame, y: str, numeric: list[str], cats: list[str]) -> dict:243 """Same estimator as hazards.cloglog (GLM cloglog, dealer-clustered) without writing to results/hazards."""244 t0 = time.time()245 d = d.dropna(subset=numeric + [y]).copy()246 X = hazards.design(d, numeric, cats); X = X.loc[:, (X != 0).any(axis=0)]247 model = sm.GLM(d[y].values.astype(float), X, family=sm.families.Binomial(link=sm.families.links.CLogLog()))248 res = model.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["dealer_id"])[0]}, maxiter=100)249 idx = [i for i, c in enumerate(X.columns) if c.startswith("age_bin_")]250 b = res.params.values[idx]; V = res.cov_params().values[np.ix_(idx, idx)]; W = float(b @ np.linalg.solve(V, b))251 ci = res.conf_int()252 return {"n_rows": int(len(d)), "n_events": int(d[y].sum()), "n_listings": int(d["listing_id"].nunique()), "n_dealers": int(d["dealer_id"].nunique()),253 "llf": float(res.llf), "joint_age_bins": {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}, "seconds": round(time.time() - t0, 1),254 "coef": {c: {"b": float(res.params[c]), "se": float(res.bse[c]), "p": float(res.pvalues[c]), "lo": float(ci.loc[c, 0]), "hi": float(ci.loc[c, 1])} for c in HAZ_TERMS + ["prior_cuts", "pricepos_25", "pricepos_entry"] if c in X.columns}}255256257def _subsample(d: pd.DataFrame, max_rows: int = MAX_ROWS) -> pd.DataFrame:258 if len(d) <= max_rows:259 return d260 return d.sample(n=max_rows, random_state=SEED)261262263def _apply_events(d: pd.DataFrame, per: pd.DataFrame, cd: pd.DataFrame, drop_new_oscillating: bool) -> pd.DataFrame:264 """Rebuild first-cut risk set, first_cut, cut and prior_cuts from re-derived events."""265 d = d.drop(columns=[c for c in ("first_cut_day_v", "n_cuts_v", "oscillating_v", "entry_price_v", "first_cut_pct_v", "cut_v") if c in d.columns])266 d = d.merge(per, on="listing_id", how="left").merge(cd, on=["listing_id", "t"], how="left")267 if drop_new_oscillating:268 d = d[~d["oscillating_v"].fillna(False).astype(bool)]269 d["first_cut_day"] = d["first_cut_day_v"]270 d["in_first_cut_risk"] = d["first_cut_day"].isna() | (d["t"] <= d["first_cut_day"])271 d["first_cut"] = (d["t"] == d["first_cut_day"]).astype(float)272 d["cut"] = d["cut_v"].fillna(0.0)273 d = d.sort_values(["listing_id", "t"])274 d["prior_cuts"] = d.groupby("listing_id")["cut"].cumsum() - d["cut"]275 return d276277278def run_hazards() -> dict:279 OUT.mkdir(parents=True, exist_ok=True)280 out = {"subsample": {"max_rows": MAX_ROWS, "seed": SEED, "unit": "listing-days, simple random sample without replacement"}}281 _log("hazards: loading S2L risk set")282 full = hazards.load("s2l") # all rows of the risk set (S2L universe)283 # NOTE: in risk_set.parquet, `first_cut` is NULL for listings that never get a cut (t == NULL first_cut_day).284 # hazards.cloglog drops NaN outcomes, so the published Model B was fitted on listings that eventually got a cut only285 # (n_listings == n_events in results/hazards/summary.json). Here NULL first_cut is coded 0 (listing-day without a cut).286 out["note_first_cut_null"] = {"risk_set_first_cut_null_rows_S2": int(full.loc[full["s2_dynamic"], "first_cut"].isna().sum()),287 "published_B_n_rows": None, "published_B_n_listings": None, "published_B_n_events": None}288 pub = json.load(open(C.RESULTS / "hazards" / "summary.json"))289 out["note_first_cut_null"].update({"published_B_n_rows": pub["B_first_markdown_cloglog"]["n_rows"], "published_B_n_listings": pub["B_first_markdown_cloglog"]["n_listings"],290 "published_B_n_events": pub["B_first_markdown_cloglog"]["n_events"]})291 full["first_cut"] = full["first_cut"].fillna(0).astype(float)292 full["exit"] = full["exit"].astype(float); full["cut"] = full["cut"].astype(float)293 L = pl.read_parquet(C.PROCESSED / "listings.parquet").select(["listing_id", "last_seen_ts", "active"]).to_pandas()294 full = full.merge(L, on="listing_id", how="left")295 full["last_t"] = full.groupby("listing_id")["t"].transform("max")296 base_events = pl.read_parquet(C.PROCESSED / "price_events.parquet")297 S2 = full[full["s2_dynamic"]].copy()298 # reproduction check of exit coding (7-day rule): must equal the stored `exit`299 ex7 = ((S2["t"] == S2["last_t"]) & (S2["active"] == 0) & (S2["last_seen_ts"] <= C.SNAPSHOT_TS - 7 * 86400)).astype(float)300 out["exit_recode_check"] = {"stored_exits_S2": int(S2["exit"].sum()), "recoded_7d_exits_S2": int(ex7.sum()), "disagreements": int((ex7 != S2["exit"]).sum())}301 _log(f"hazards: S2 rows={len(S2)} exits={int(S2['exit'].sum())} recode check={out['exit_recode_check']}")302303 results = {"B": {}, "C": {}}304 def fitBC(name: str, d: pd.DataFrame, note: str):305 dB = _subsample(d[d["in_first_cut_risk"]]); dC = _subsample(d)306 rB = _cloglog(dB, "first_cut", NUM_B, CATS); rB["note"] = note; rB["n_rows_before_subsample"] = int(d["in_first_cut_risk"].sum())307 rB["n_events_before_subsample"] = int(d.loc[d["in_first_cut_risk"], "first_cut"].sum())308 results["B"][name] = rB309 _log(f"B {name}: n={rB['n_rows']} ev={rB['n_events']} ageP={rB['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rB['coef'][k]['b']:+.3f}({rB['coef'][k]['se']:.3f})" for k in HAZ_TERMS))310 rC = _cloglog(dC, "exit", NUM_C, CATS); rC["note"] = note; rC["n_rows_before_subsample"] = int(len(d)); rC["n_events_before_subsample"] = int(d["exit"].sum())311 results["C"][name] = rC312 _log(f"C {name}: n={rC['n_rows']} ev={rC['n_events']} ageP={rC['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rC['coef'][k]['b']:+.3f}({rC['coef'][k]['se']:.3f})" for k in HAZ_TERMS))313314 # published baselines (full sample)315 results["B"]["published_full"] = {k: pub["B_first_markdown_cloglog"][k] for k in ("n_rows", "n_events", "n_listings", "joint_age_bins")}316 results["B"]["published_full"]["coef"] = {k: pub["B_first_markdown_cloglog"]["coef"][k] for k in HAZ_TERMS}317 results["B"]["published_full"]["note"] = "results/hazards/summary.json as published (fitted on listings with a NULL-free first_cut, i.e. eventually-cut listings only)"318 results["C"]["published_full"] = {k: pub["C_exit_cloglog"][k] for k in ("n_rows", "n_events", "n_listings", "joint_age_bins")}319 results["C"]["published_full"]["coef"] = {k: pub["C_exit_cloglog"]["coef"][k] for k in HAZ_TERMS}320 results["C"]["published_full"]["note"] = "results/hazards/summary.json as published (full S2 risk set)"321 # corrected B on the full first-cut risk set (no subsample) — reference for the NULL-coding issue322 rBfull = _cloglog(S2[S2["in_first_cut_risk"]], "first_cut", NUM_B, CATS); rBfull["note"] = "full S2 first-cut risk set, NULL first_cut coded 0, no subsample"323 results["B"]["baseline_full_nullfix"] = rBfull324 _log(f"B baseline_full_nullfix: n={rBfull['n_rows']} ev={rBfull['n_events']} ageP={rBfull['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rBfull['coef'][k]['b']:+.3f}({rBfull['coef'][k]['se']:.3f})" for k in HAZ_TERMS))325326 # --- baseline rules on the subsample ---------------------------------------------327 fitBC("baseline_sub", S2, "S2; 1 %/100 $; 48 h; exit 7 d; subsample")328 # reproduction check: re-deriving events under the baseline rules must give the stored first_cut_day329 per0, _ = rebuild_events(C.TRANSIENT_SECONDS, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS)330 stored = S2.drop_duplicates("listing_id").set_index("listing_id")["first_cut_day"]331 chk = per0.set_index("listing_id")["first_cut_day_v"].reindex(stored.index)332 out["events_baseline_reproduction"] = {"listings_S2": int(len(stored)), "first_cut_day_disagreements": int(((chk.fillna(-1) != stored.fillna(-1))).sum()),333 "listings_with_cut_S2L": int(per0["first_cut_day_v"].notna().sum())}334 _log(f"events baseline reproduction: {out['events_baseline_reproduction']}")335 # --- (a) markdown thresholds, (b) transient rule -------------------------------------336 for name, (tr, pct, ab), drop_osc in [("threshold_2pct_200", (C.TRANSIENT_SECONDS, 0.02, 200.0), False), ("threshold_0.5pct_50", (C.TRANSIENT_SECONDS, 0.005, 50.0), False),337 ("transient_24h", (24 * 3600, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS), True), ("transient_72h", (72 * 3600, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS), True)]:338 per, cd = rebuild_events(tr, pct, ab)339 d = _apply_events(S2, per, cd, drop_new_oscillating=drop_osc)340 info = {"transient_seconds": tr, "min_pct": pct, "min_abs": ab, "listings_with_cut": int(per["first_cut_day_v"].notna().sum()), "oscillating_under_rule": int(per["oscillating_v"].sum()),341 "listings_dropped_as_newly_oscillating": int(S2["listing_id"].nunique() - d["listing_id"].nunique()) if drop_osc else 0,342 "entry_price_changed_share": float((np.abs(per["entry_price_v"] - per["listing_id"].map(S2.drop_duplicates("listing_id").set_index("listing_id")["entry_price"])) > 0.5).mean())}343 out[f"events_{name}"] = info344 _log(f"events {name}: {info}")345 fitBC(name, d, f"S2; rebuilt events ({tr/3600:.0f} h, {pct*100:g} %/{ab:g} $); lprice/resid_A2 from baseline entry price; subsample")346 # --- (c) exit confirmation 3 / 14 days ----------------------------------------------347 for k in (3, 14):348 d = S2.copy()349 d["exit"] = ((d["t"] == d["last_t"]) & (d["active"] == 0) & (d["last_seen_ts"] <= C.SNAPSHOT_TS - k * 86400)).astype(float)350 out[f"exit_confirm_{k}d"] = {"confirmed_exits_S2": int(d["exit"].sum()), "baseline_7d": int(S2["exit"].sum())}351 fitBC(f"exit_confirm_{k}d", d, f"S2; exit = inactive and last_seen <= snapshot - {k} d (exit on last observed listing-day); subsample")352 # --- (d) S2L ---------------------------------------------------------------------353 out["S2L"] = {"rows": int(len(full)), "listings": int(full["listing_id"].nunique()), "left_censored_listings": int(full.loc[full["left_censored"], "listing_id"].nunique()),354 "note": "left-censored listings enter with listing_age measured from first_seen (lower bound)"}355 fitBC("S2L_with_stock", full, "S2L (S2 + left-censored stock, age from first_seen); subsample")356 out["results"] = results357 # flat CSV358 rows = []359 for mdl in ("B", "C"):360 for nm, r in results[mdl].items():361 row = {"model": mdl, "variant": nm, "n_rows": r["n_rows"], "n_events": r["n_events"], "n_listings": r["n_listings"], "joint_age_stat": r["joint_age_bins"]["stat"], "joint_age_p": r["joint_age_bins"]["p"]}362 for k in HAZ_TERMS:363 c = r["coef"].get(k, {})364 row.update({f"{k}_b": c.get("b"), f"{k}_se": c.get("se"), f"{k}_p": c.get("p")})365 row["note"] = r.get("note"); rows.append(row)366 pd.DataFrame(rows).to_csv(OUT / "hazards_variants.csv", index=False)367 return out368369370# =============================================================================== 3. AMLI371372def _validate(Z: pd.DataFrame, key: list[str], zcols: list[str], eval_weeks) -> dict:373 """Leave-exit-out (equal-weight mean of zcols) vs next-period exit, market FE, cluster by market."""374 Z = Z.copy(); Z["idx"] = Z[zcols].mean(axis=1)375 Z["mkt"] = Z[key].astype(str).agg("|".join, axis=1)376 Z = Z.sort_values(["mkt", "week"]); Z["exit_next"] = Z.groupby("mkt")["exit"].shift(-1)377 d = Z[Z["week"].isin(list(eval_weeks))].dropna(subset=["idx", "exit_next", "exit"])378 nm = d["mkt"].nunique()379 if nm > 1:380 m = pf.feols("exit_next ~ idx | mkt", data=d, vcov={"CRV1": "mkt"})381 else:382 m = pf.feols("exit_next ~ idx", data=d, vcov="hetero")383 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); r = t[t["term"] == "idx"].iloc[0]384 b, se = float(r["Estimate"]), float(r["Std. Error"])385 return {"b": b, "se": se, "lo": float(r["2.5%"]), "hi": float(r["97.5%"]), "p_normal": float(2 * norm.sf(abs(b / se))) if se > 0 else None, "n": int(m._N), "n_markets": int(nm), "r2": float(m._r2)}386387388def _rank_index(comp: pd.DataFrame, cols: list[str]) -> pd.DataFrame:389 """Rank standardisation: signed normal scores of the within-panel rank (replaces z-scores)."""390 Z = comp.copy()391 for c in cols:392 r = Z[c].rank(); n = Z[c].notna().sum()393 Z[f"z_{c}"] = amli.SIGN[c] * norm.ppf((r - 0.5) / n)394 zc = [f"z_{c}" for c in cols]395 ok = Z[zc].notna().sum(axis=1) >= min(3, len(cols))396 Z["amli_equal"] = Z.loc[ok, zc].mean(axis=1)397 return Z398399400def run_amli() -> dict:401 OUT.mkdir(parents=True, exist_ok=True)402 out, rows = {}, []403 d = amli.listing_days()404 others = [f"z_{c}" for c in amli.PRIMARY if c != "exit"]405 ev_w = range(6, amli.WEEKS)406 comp = amli.components(d, ["region"])407 Z, meta = amli.build_index(comp, amli.PRIMARY)408 base = _validate(Z, ["region"], others, ev_w)409 out["baseline"] = {"validation_noexit_fe": base, "pca_explained": meta["weights"]["pca_explained"], "n_market_weeks": int(Z["amli_equal"].notna().sum()),410 "min_listing_days": amli.MIN_LISTING_DAYS, "weeks": [int(Z["week"].min()), int(Z["week"].max())]}411 rows.append({"variant": "baseline_z_region_weekly", **base, "corr_equal_with_baseline": 1.0})412 _log(f"AMLI baseline: slope={base['b']:.5f} ({base['se']:.5f}) n={base['n']}")413 # (a) rank standardisation414 Zr = _rank_index(comp, amli.PRIMARY)415 vr = _validate(Zr, ["region"], others, ev_w)416 out["rank_standardisation"] = {"validation_noexit_fe": vr, "corr_equal_with_baseline": float(Zr["amli_equal"].corr(Z["amli_equal"])),417 "spearman_equal_with_baseline": float(Zr["amli_equal"].corr(Z["amli_equal"], method="spearman"))}418 rows.append({"variant": "rank_standardisation", **vr, "corr_equal_with_baseline": out["rank_standardisation"]["corr_equal_with_baseline"]})419 _log(f"AMLI rank: slope={vr['b']:.5f} ({vr['se']:.5f}) corr={out['rank_standardisation']['corr_equal_with_baseline']:.3f}")420 # (b) MIN_LISTING_DAYS = 100421 comp100 = comp.copy(); comp100.loc[comp100["ld"] < 100, amli.PRIMARY + ["stale"]] = np.nan422 Z100, m100 = amli.build_index(comp100, amli.PRIMARY)423 v100 = _validate(Z100, ["region"], others, ev_w)424 out["min_listing_days_100"] = {"validation_noexit_fe": v100, "n_market_weeks": int(Z100["amli_equal"].notna().sum()), "n_market_weeks_baseline": int(Z["amli_equal"].notna().sum()),425 "cells_between_30_and_100_listing_days": int(((comp["ld"] >= 30) & (comp["ld"] < 100)).sum()), "min_listing_days_in_panel": float(comp["ld"].min()),426 "ld_quantiles": {q: float(comp["ld"].quantile(q)) for q in (0.0, 0.05, 0.25, 0.5)},427 "corr_equal_with_baseline": float(Z100["amli_equal"].corr(Z["amli_equal"])), "pca_explained": m100["weights"]["pca_explained"]}428 rows.append({"variant": "min_listing_days_100", **v100, "corr_equal_with_baseline": out["min_listing_days_100"]["corr_equal_with_baseline"]})429 _log(f"AMLI minld100: slope={v100['b']:.5f} ({v100['se']:.5f}) mw={out['min_listing_days_100']['n_market_weeks']}")430 # (c) 2-week aggregation (periods: weeks 2-3, 4-5, 6-7, 8-9, 10)431 d2 = d.copy(); d2["week"] = (d2["week"] - amli.WEEK_MIN) // 2432 comp2 = amli.components(d2, ["region"])433 Z2, m2 = amli.build_index(comp2, amli.PRIMARY, fit_weeks=range(0, 2))434 v2 = _validate(Z2, ["region"], others, range(2, 5))435 # comparable weekly index averaged over the same 2-week periods436 Zw = Z.copy(); Zw["period"] = (Zw["week"] - amli.WEEK_MIN) // 2437 avg = Zw.groupby(["region", "period"])["amli_equal"].mean().rename("amli_weekly_avg").reset_index().rename(columns={"period": "week"})438 cmp2 = Z2.merge(avg, on=["region", "week"], how="left")439 out["two_week_aggregation"] = {"validation_noexit_fe": v2, "periods": "weeks 2-3, 4-5, 6-7, 8-9, 10 (last period = 1 week)", "fit_periods": [0, 1], "eval_periods": [2, 3, 4],440 "corr_with_weekly_index_averaged": float(cmp2["amli_equal"].corr(cmp2["amli_weekly_avg"])), "pca_explained": m2["weights"]["pca_explained"],441 "n_market_periods": int(Z2["amli_equal"].notna().sum())}442 rows.append({"variant": "two_week_aggregation", **v2, "corr_equal_with_baseline": out["two_week_aggregation"]["corr_with_weekly_index_averaged"]})443 _log(f"AMLI 2-week: slope={v2['b']:.5f} ({v2['se']:.5f}) n={v2['n']}")444 # (d) province vs region445 dP = d.copy(); dP["market"] = "Quebec"446 compP = amli.components(dP, ["market"]); ZP, mP = amli.build_index(compP, amli.PRIMARY)447 wavg = Z.dropna(subset=["amli_equal"]).assign(wx=lambda x: x["amli_equal"] * x["ld"]).groupby("week").agg(wx=("wx", "sum"), ld=("ld", "sum"))448 wavg["region_weighted_equal"] = wavg["wx"] / wavg["ld"]449 cmpP = ZP.set_index("week")[["amli_equal", "amli_pca"]].join(wavg["region_weighted_equal"])450 vP = _validate(ZP, ["market"], others, ev_w)451 vP_all = _validate(ZP, ["market"], others, range(amli.WEEK_MIN, amli.WEEKS))452 out["province_vs_region"] = {"corr_province_equal_vs_region_weighted_equal": float(cmpP["amli_equal"].corr(cmpP["region_weighted_equal"])),453 "spearman": float(cmpP["amli_equal"].corr(cmpP["region_weighted_equal"], method="spearman")), "n_weeks": int(cmpP["amli_equal"].notna().sum()),454 "province_validation_pooled_eval_weeks": vP, "province_validation_pooled_all_weeks": vP_all, "region_validation": base,455 "note": "province index has one market: pooled OLS with HC SEs on 4 (eval) / 8 (all) weekly points, descriptive only",456 "series": cmpP.reset_index().to_dict(orient="list")}457 rows.append({"variant": "province_level_pooled", **vP, "corr_equal_with_baseline": out["province_vs_region"]["corr_province_equal_vs_region_weighted_equal"]})458 _log(f"AMLI province: corr={out['province_vs_region']['corr_province_equal_vs_region_weighted_equal']:.3f} slope(pooled)={vP['b']:.5f} n={vP['n']}")459 # (e) leave-one-component-out predictive validation (on top of leave-exit-out)460 loo = {}461 for c in [x for x in amli.PRIMARY if x != "exit"]:462 cols = [z for z in others if z != f"z_{c}"]463 loo[c] = {"baseline_z": _validate(Z, ["region"], cols, ev_w), "rank": _validate(Zr, ["region"], cols, ev_w), "min_ld_100": _validate(Z100, ["region"], cols, ev_w),464 "two_week": _validate(Z2, ["region"], cols, range(2, 5))}465 rows.append({"variant": f"loo_drop_{c}", **loo[c]["baseline_z"], "corr_equal_with_baseline": None})466 _log(f"AMLI LOO drop {c}: slope={loo[c]['baseline_z']['b']:.5f} ({loo[c]['baseline_z']['se']:.5f})")467 out["leave_one_component_out"] = loo468 pubA = json.load(open(C.RESULTS / "amli" / "summary.json"))469 out["published"] = {"region_validation_noexit_fe": pubA["region_validation"]["noexit_fe"], "listing_days": pubA["listing_days"]}470 pd.DataFrame(rows).to_csv(OUT / "amli_variants.csv", index=False)471 return out472473474# =============================================================================== 4. Holm475476def holm_primary_tests() -> pd.DataFrame:477 hed = json.load(open(C.RESULTS / "hedonic" / "summary.json")); haz = json.load(open(C.RESULTS / "hazards" / "summary.json"))478 comp = json.load(open(C.RESULTS / "competition" / "summary.json")); am = json.load(open(C.RESULTS / "amli" / "summary.json"))479 v = am["region_validation"]["noexit_fe"]; p6 = float(2 * norm.sf(abs(v["b"] / v["se"])))480 tests = [481 ("H1", "hedonic A2: joint Wald on nonlinear age/km spline terms", "results/hedonic/summary.json tests.H1_nonlinear_joint", hed["tests"]["H1_nonlinear_joint"]["p"]),482 ("H2", "Model B cloglog: joint test of listing-age bins", "results/hazards/summary.json B_first_markdown_cloglog.joint_age_bins", haz["B_first_markdown_cloglog"]["joint_age_bins"]["p"]),483 ("H3", "Model D LPM 25 km: joint test of competitor-cut lags", "results/competition/summary.json main_25km.lags_joint", comp["main_25km"]["lags_joint"]["p"]),484 ("H4", "hedonic A4: joint Wald on region effects", "results/hedonic/summary.json tests.H4_region_joint", hed["tests"]["H4_region_joint"]["p"]),485 ("H5", "Model C cloglog: coefficient on entry hedonic residual (resid_A2)", "results/hazards/summary.json C_exit_cloglog.coef.resid_A2", haz["C_exit_cloglog"]["coef"]["resid_A2"]["p"]),486 ("H6", "AMLI leave-exit-out index vs next-week exit (region FE), normal p from b/se", "results/amli/summary.json region_validation.noexit_fe", p6),487 ]488 df = pd.DataFrame(tests, columns=["hypothesis", "test", "source", "p_raw"])489 m = len(df); df = df.sort_values("p_raw").reset_index(drop=True); df["rank"] = np.arange(1, m + 1)490 # Holm step-down: adj_(i) = max_{j<=i} min(1, (m-j+1) p_(j))491 adj = np.maximum.accumulate(np.minimum(1.0, (m - df["rank"].values + 1) * df["p_raw"].values))492 df["p_holm"] = adj; df["reject_0.05_raw"] = df["p_raw"] < 0.05; df["reject_0.05_holm"] = df["p_holm"] < 0.05493 df = df.sort_values("hypothesis").reset_index(drop=True)494 df.to_csv(OUT / "holm_primary_tests.csv", index=False)495 return df496497498# =============================================================================== 5. table499500def _fp(p) -> str:501 if p is None or not np.isfinite(p):502 return "--"503 return "$<$0.001" if p < 0.001 else f"{p:.3f}"504505506def _fe(b, se=None, d=3) -> str:507 if b is None or not np.isfinite(b):508 return "--"509 s = f"{b:.{d}f}"510 return s + (f" ({se:.{d}f})" if se is not None and np.isfinite(se) else "")511512513def _delta(a, b, d=3) -> str:514 if a is None or b is None or not (np.isfinite(a) and np.isfinite(b)):515 return "--"516 return f"{b - a:+.{d}f}"517518519def _verdict(pa, pb) -> str:520 if pa is None or pb is None:521 return "--"522 return "same" if (pa < 0.05) == (pb < 0.05) else "flips"523524525def _esc(s: str) -> str:526 return s.replace("%", r"\%").replace("$", r"\$").replace("_", r"\_")527528529def write_table(S: dict | None = None) -> None:530 if S is None:531 S = json.load(open(OUT / "summary.json"))532 H, Z, A, holm = S["hedonic"], S["hazards"]["results"], S["amli"], S["holm"]533 holm = pd.DataFrame(holm) if not isinstance(holm, pd.DataFrame) else holm534 b = H["baseline_S1"]535 rows = []536 def hed_row(label, key, what):537 v = H[key]538 if what == "age":539 rows.append((label, "age slope (log \\$/yr)", _fe(b["age"]["b"], b["age"]["se"], 4), _fe(v["age"]["b"], v["age"]["se"], 4), _delta(b["age"]["b"], v["age"]["b"], 4)))540 elif what == "km":541 rows.append((label, "km slope (log \\$/10k km)", _fe(b["km"]["b"], b["km"]["se"], 4), _fe(v["km"]["b"], v["km"]["se"], 4), _delta(b["km"]["b"], v["km"]["b"], 4)))542 elif what == "H4":543 rows.append((label, "H4 region joint $p$ (A4 trim FE / A2)", f"{_fp(b['H4_region_joint_A4']['p'])} / {_fp(b['H4_region_joint']['p'])}", f"{_fp(v['H4_region_joint_A4']['p'])} / {_fp(v['H4_region_joint']['p'])}", _verdict(b["H4_region_joint_A4"]["p"], v["H4_region_joint_A4"]["p"])))544 elif what == "H1":545 rows.append((label, "H1 nonlinearity Wald $p$", _fp(b["H1_nonlinear_joint"]["p"]), _fp(v["H1_nonlinear_joint"]["p"]), _verdict(b["H1_nonlinear_joint"]["p"], v["H1_nonlinear_joint"]["p"])))546 rows.append(("\\textit{Hedonic (Model A2 splines, make--model FE, dealer-clustered)}",))547 hed_row("+ sampled portals (autotrader, cargurus)", "plus_sampled_portals", "age"); hed_row("+ sampled portals", "plus_sampled_portals", "H4")548 hed_row("+ Kijiji private sellers", "plus_kijiji", "age"); hed_row("+ Kijiji private sellers", "plus_kijiji", "H4")549 hed_row("Price filter 1\\,000--500\\,000\\,\\$", "price_1000_500000", "age"); hed_row("Price filter 1\\,000--500\\,000\\,\\$", "price_1000_500000", "H4")550 hed_row("Reweighted to otogo make--model--year mix", "reweighted_to_otogo_mmy", "age"); hed_row("Reweighted to otogo make--model--year mix", "reweighted_to_otogo_mmy", "H4")551 het = H["heterogeneity"]552 def het_row(label, k1, k2, what):553 a, c = het.get(k1), het.get(k2)554 if a and c:555 l1, l2 = _esc(k1.split("=", 1)[1]).replace("<=", "$\\leq$"), _esc(k2.split("=", 1)[1]).replace("<=", "$\\leq$")556 rows.append((label, f"{what} slope: {l1} vs {l2}", _fe(a[what]["b"], a[what]["se"], 4), _fe(c[what]["b"], c[what]["se"], 4), _delta(a[what]["b"], c[what]["b"], 4)))557 het_row("Heterogeneity: powertrain", "powertrain=ICE", "powertrain=BEV", "age")558 het_row("Heterogeneity: model-year band", "yband=<=2015", "yband=2021+", "age")559 het_row("Heterogeneity: body", "body=SUV", "body=Pickup", "km")560 het_row("Heterogeneity: source class", "source_class=direct_dealer", "source_class=full_portal", "age")561 n_hed = len(rows)562 # hazards563 rows.append(("\\textit{Hazards (cloglog, dealer-clustered; $\\leq$800k listing-day subsample)}",))564 bB, bC = Z["B"]["baseline_sub"], Z["C"]["baseline_sub"]565 def haz_row(label, model, key, term):566 base = bB if model == "B" else bC; v = Z[model][key]567 if term == "age":568 rows.append((label, f"{model}: listing-age bins joint $p$", _fp(base["joint_age_bins"]["p"]), _fp(v["joint_age_bins"]["p"]), _verdict(base["joint_age_bins"]["p"], v["joint_age_bins"]["p"])))569 else:570 rows.append((label, f"{model}: coef.\\ {_esc(term)}", _fe(base["coef"][term]["b"], base["coef"][term]["se"]), _fe(v["coef"][term]["b"], v["coef"][term]["se"]), _delta(base["coef"][term]["b"], v["coef"][term]["b"])))571 # (rows comparing the pre-E-016 survivor-only Model B with the full risk set were removed once the572 # pipeline bug was fixed: the two estimates are now identical by construction)573 haz_row("Markdown $\\geq$2\\,\\%/200\\,\\$", "B", "threshold_2pct_200", "ldens")574 haz_row("Markdown $\\geq$0.5\\,\\%/50\\,\\$", "B", "threshold_0.5pct_50", "ldens")575 haz_row("Transient rule 24\\,h", "B", "transient_24h", "ldens"); haz_row("Transient rule 72\\,h", "B", "transient_72h", "ldens")576 haz_row("Exit confirmation 3\\,d", "C", "exit_confirm_3d", "resid_A2")577 haz_row("Exit confirmation 14\\,d", "C", "exit_confirm_14d", "resid_A2")578 haz_row("S2L (left-censored stock included)", "C", "S2L_with_stock", "resid_A2"); haz_row("S2L (left-censored stock included)", "C", "S2L_with_stock", "ldens")579 # AMLI580 rows.append(("\\textit{AMLI (leave-exit-out index $\\to$ next-week exit intensity, region FE)}",))581 bA = A["baseline"]["validation_noexit_fe"]582 def am_row(label, v):583 rows.append((label, "slope $\\hat\\beta$ (s.e.)", _fe(bA["b"], bA["se"], 4), _fe(v["b"], v["se"], 4), _delta(bA["b"], v["b"], 4)))584 am_row("Rank standardisation", A["rank_standardisation"]["validation_noexit_fe"])585 am_row("Min.\\ 100 listing-days per cell", A["min_listing_days_100"]["validation_noexit_fe"])586 am_row("Two-week aggregation", A["two_week_aggregation"]["validation_noexit_fe"])587 loo = A["leave_one_component_out"]; worst = min(loo, key=lambda c: loo[c]["baseline_z"]["b"]); best = max(loo, key=lambda c: loo[c]["baseline_z"]["b"])588 am_row(f"Leave-one-component-out (min: drop {_esc(worst)})", loo[worst]["baseline_z"]); am_row(f"Leave-one-component-out (max: drop {_esc(best)})", loo[best]["baseline_z"])589 pv = A["province_vs_region"]590 rows.append(("Province vs region aggregation", "corr.\\ equal-weight index (weekly)", "1.000", f"{pv['corr_province_equal_vs_region_weighted_equal']:.3f}", "--"))591 # Holm592 rows.append(("\\textit{Multiple testing (six primary tests)}",))593 hm = holm594 rows.append(("Holm correction", "tests rejected at 5\\,\\% (raw $\\to$ Holm)", str(int(hm["reject_0.05_raw"].sum())), str(int(hm["reject_0.05_holm"].sum())), f"{int(hm['reject_0.05_raw'].sum() - hm['reject_0.05_holm'].sum())} lost"))595 assert sum(1 for r in rows if len(r) == 5) <= 30, "table exceeds 30 data rows"596 lines = ["\\begin{tabular}{llrrr}", "\\toprule", "Check & Quantity & Baseline & Variant & Change \\\\", "\\midrule"]597 for i, r in enumerate(rows):598 if len(r) == 1:599 if i > 0:600 lines.append("\\addlinespace")601 lines.append(f"\\multicolumn{{5}}{{l}}{{{r[0]}}} \\\\")602 else:603 lines.append(" & ".join(r) + " \\\\")604 lines += ["\\bottomrule", "\\end{tabular}"]605 C.TABLES.mkdir(parents=True, exist_ok=True)606 (C.TABLES / "tab09_robustness.tex").write_text("\n".join(lines) + "\n")607608609# =============================================================================== driver610611def run() -> dict:612 OUT.mkdir(parents=True, exist_ok=True)613 t0 = time.time(); S = {"generated": time.strftime("%Y-%m-%d %H:%M:%S"), "seed": SEED}614 _log("=== 1. hedonic"); t1 = time.time(); S["hedonic"] = run_hedonic(); S["hedonic"]["seconds"] = round(time.time() - t1, 1)615 _log("=== 2. hazards"); t1 = time.time(); S["hazards"] = run_hazards(); S["hazards"]["seconds"] = round(time.time() - t1, 1)616 _log("=== 3. AMLI"); t1 = time.time(); S["amli"] = run_amli(); S["amli"]["seconds"] = round(time.time() - t1, 1)617 _log("=== 4. Holm"); S["holm"] = holm_primary_tests().to_dict(orient="records")618 S["seconds_total"] = round(time.time() - t0, 1)619 json.dump(S, open(OUT / "summary.json", "w"), indent=2, default=float)620 _log("=== 5. table"); write_table(S)621 _log(f"done in {S['seconds_total']} s")622 return S623624625if __name__ == "__main__":626 run()627