1"""Model-based figures: 2 (depreciation surface), 5 (hazards), 6 (event study), 7 (regional premia map),28 (AMLI series + components), 9 (AMLI internal validation). Descriptive figures 1, 3, 4 live in descriptives.py;3figure 10 in prediction.py. Every plotted number is read from results/ files produced by the pipeline."""4from __future__ import annotations56import json78import matplotlib.pyplot as plt9import numpy as np10import pandas as pd11import polars as pl1213from . import config as C14from . import plotstyle as ps1516R = C.RESULTS171819def fig02_surface():20 g = pd.read_csv(R / "hedonic" / "surface_A2t.csv")21 # mask grid cells without data support (fewer than 25 S1 listings within +-0.75 y and +-1.5 (10k km))22 L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s1_hedonic")).select(["age_years_entry", "mileage_km"]).to_pandas()23 a = L["age_years_entry"].values; k = L["mileage_km"].values / 1e424 H, ae, ke = np.histogram2d(a, k, bins=[np.arange(0, 16.5, 0.5), np.arange(0, 31.5, 1.0)])25 def support(age, km):26 i = min(int(age / 0.5), H.shape[0] - 1); j = min(int(km / 1.0), H.shape[1] - 1)27 return H[max(0, i - 1):i + 2, max(0, j - 1):j + 2].sum()28 g["support"] = [support(x, y) for x, y in zip(g["age"], g["km"])]29 g.loc[g["support"] < 25, "logp_rel"] = np.nan30 A = g.pivot(index="km", columns="age", values="logp_rel")31 fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.9), gridspec_kw={"width_ratios": [1.25, 1]})32 ax = axes[0]33 X, Y = np.meshgrid(A.columns.values, A.index.values * 10)34 cs = ax.contourf(X, Y, np.ma.masked_invalid(A.values), levels=np.linspace(np.nanmin(A.values), 0, 16), cmap=ps.cmap_seq.reversed())35 cl = ax.contour(X, Y, np.ma.masked_invalid(A.values), levels=[-1.5, -1.25, -1.0, -0.75, -0.5, -0.25], colors=ps.TEXT2, linewidths=0.5)36 ax.clabel(cl, fmt=lambda v: f"{100*(np.exp(v)-1):.0f} %", fontsize=7)37 ax.set_xlabel("Vehicle age (years)"); ax.set_ylabel("Mileage (thousand km)"); ax.set_title("(a) Quality-adjusted asking price relative to 1 year / 10 000 km")38 cb = fig.colorbar(cs, ax=ax, shrink=0.9, pad=0.02, ticks=[-1.5, -1.0, -0.5, 0]); cb.set_label("log price difference"); cb.ax.set_yticklabels([f"{100*(np.exp(v)-1):.0f} %" for v in [-1.5, -1.0, -0.5, 0]])39 ax.grid(False)40 ax = axes[1]41 sub = g[np.isclose(g["km"], g["km"].iloc[(g["km"] - 8).abs().argmin()])]42 for km_t, col in [(2, ps.CAT[0]), (8, ps.CAT[1]), (16, ps.CAT[2])]:43 s = g[np.isclose(g["km"], g["km"].iloc[(g["km"] - km_t).abs().argmin()])].sort_values("age")44 ax.plot(s["age"], 100 * (np.exp(s["logp_rel"]) - 1), color=col, label=f"{km_t*10000:,} km".replace(",", "\u2009"))45 ax.set_xlabel("Vehicle age (years)"); ax.set_ylabel("Price relative to 1 year / 10 000 km (%)"); ax.set_title("(b) Age profiles at fixed mileage")46 ax.legend(title="Mileage", loc="upper right")47 ps.save(fig, "fig02_depreciation_surface")484950def fig05_hazards():51 hm = pd.read_csv(R / "hazards" / "empirical_hazard_first_markdown.csv"); he = pd.read_csv(R / "hazards" / "empirical_hazard_exit.csv")52 fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.5))53 for ax, d, title, col in [(axes[0], hm, "(a) First sustained markdown", ps.CAT[0]), (axes[1], he, "(b) Observed listing exit", ps.CAT[1])]:54 d = d[d["listing_age"] >= 1]55 ax.fill_between(d["listing_age"], 100 * d["lo"], 100 * d["hi"], color=col, alpha=0.2, linewidth=0)56 ax.plot(d["listing_age"], 100 * d["hazard"], color=col)57 ax.set_xlabel("Listing age (days)"); ax.set_ylabel("Daily hazard (%)"); ax.set_title(title); ax.set_xlim(0, 70)58 ax = axes[2]59 order = ["direct_dealer", "full_portal", "banner"]60 for i, sc in enumerate(order):61 p = R / "hazards" / f"empirical_hazard_exit_{sc}.csv"62 if p.exists():63 d = pd.read_csv(p); d = d[(d["listing_age"] >= 1) & (d["at_risk"] >= 200)]64 S = np.cumprod(1 - d["hazard"].values)65 ax.plot(d["listing_age"], 100 * S, color=ps.color_for(ps.CLASS_ORDER, sc), label=ps.CLASS_LABEL[sc])66 ax.set_xlabel("Listing age (days)"); ax.set_ylabel("Still listed (%)"); ax.set_title("(c) Survival by source class"); ax.legend(); ax.set_xlim(0, 70)67 ps.save(fig, "fig05_hazard_curves")686970def fig06_eventstudy():71 main = pd.read_csv(R / "competition" / "eventstudy_main.csv"); pla = pd.read_csv(R / "competition" / "eventstudy_placebo_permuted.csv")72 fc = pd.read_csv(R / "competition" / "eventstudy_firstcut.csv")73 dd = R / "competition" / "eventstudy_dealer_dow_fe.csv"74 fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.7), sharey=False)75 ax = axes[0]76 for d, col, lab, off in [(main, ps.CAT[0], "Main (market×week, model-year×week FE)", -0.15), (pla, ps.NEUTRAL, "Placebo (competitors permuted within make-model)", 0.15)]:77 ax.errorbar(d["k"] + off, 100 * d["Estimate"], yerr=100 * (d["97.5%"] - d["Estimate"]), fmt="o", ms=3, color=col, ecolor=col, elinewidth=0.8, capsize=0, label=lab)78 if dd.exists():79 d = pd.read_csv(dd); ax.plot(d["k"], 100 * d["Estimate"], color=ps.CAT[1], lw=1.2, label="+ dealer×weekday FE")80 ax.axhline(0, color=ps.TEXT2, lw=0.6); ax.axvline(-0.5, color=ps.GRID, lw=0.8)81 ax.set_xlabel("Days relative to competitor's price cut (k)"); ax.set_ylabel("Own daily cut probability (pp)"); ax.set_title("(a) Any sustained own cut")82 ax.legend(fontsize=6.5, loc="upper left")83 ax = axes[1]84 ax.errorbar(fc["k"], 100 * fc["Estimate"], yerr=100 * (fc["97.5%"] - fc["Estimate"]), fmt="o", ms=3, color=ps.CAT[2], ecolor=ps.CAT[2], elinewidth=0.8, label="First markdown risk set")85 ax.axhline(0, color=ps.TEXT2, lw=0.6); ax.set_xlabel("Days relative to competitor's price cut (k)"); ax.set_ylabel("First-cut daily probability (pp)"); ax.set_title("(b) First sustained markdown"); ax.legend(fontsize=6.5)86 ps.save(fig, "fig06_competitor_event_study")878889def fig07_region_map():90 prem = pd.read_csv(R / "hedonic" / "region_premia_A4.csv"); prem = prem[prem["region"] != "(missing)"]91 prem = pd.concat([prem, pd.DataFrame([{"region": "Montréal", "Estimate": 0.0, "Std. Error": 0.0, "2.5%": 0.0, "97.5%": 0.0}])])92 L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s1_hedonic") & pl.col("lat_r").is_between(44, 63) & pl.col("lng_r").is_between(-80, -56))93 cen = L.group_by("region").agg(pl.col("lat_r").median().alias("lat"), pl.col("lng_r").median().alias("lng"), pl.len().alias("n")).to_pandas()94 d = prem.merge(cen, on="region")95 fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.2), gridspec_kw={"width_ratios": [1.1, 1]})96 ax = axes[0]97 pts = L.select(["lat_r", "lng_r"]).to_pandas()98 ax.hexbin(pts["lng_r"], pts["lat_r"], gridsize=45, cmap=ps.cmap_seq, mincnt=1, linewidths=0.1, bins="log", alpha=0.55)99 vmax = max(abs(d["Estimate"].min()), abs(d["Estimate"].max()), 0.03)100 sc = ax.scatter(d["lng"], d["lat"], s=18 + 260 * d["n"] / d["n"].max(), c=100 * d["Estimate"], cmap=ps.cmap_div, vmin=-100 * vmax, vmax=100 * vmax, edgecolor=ps.TEXT, linewidth=0.5, zorder=3)101 for _, r in d.iterrows():102 ax.annotate(r["region"].split("–")[0][:12], (r["lng"], r["lat"]), fontsize=5, color=ps.TEXT2, xytext=(4, 2), textcoords="offset points")103 cb = fig.colorbar(sc, ax=ax, shrink=0.8, pad=0.02); cb.set_label("Quality-adjusted premium vs Montréal (%)")104 ax.set_xlabel("Longitude"); ax.set_ylabel("Latitude"); ax.set_title("(a) Listing locations and regional premia"); ax.grid(False)105 ax = axes[1]106 d2 = d.sort_values("Estimate")107 ax.errorbar(100 * d2["Estimate"], range(len(d2)), xerr=100 * (d2["97.5%"] - d2["Estimate"]), fmt="o", ms=3, color=ps.CAT[0], ecolor=ps.CAT[0], elinewidth=0.8)108 ax.set_yticks(range(len(d2))); ax.set_yticklabels(d2["region"], fontsize=6.5); ax.axvline(0, color=ps.TEXT2, lw=0.6)109 ax.set_xlabel("Premium vs Montréal (%), 95 % CI"); ax.set_title("(b) Model A4 region effects")110 ps.save(fig, "fig07_regional_premia_map")111112113def fig08_amli():114 Z = pd.read_csv(R / "amli" / "index_region_week.csv"); P = pd.read_csv(R / "amli" / "index_province_week.csv")115 B = R / "amli" / "bootstrap_region_week.csv"116 comps = ["exit", "mdfreq", "mddepth", "disp", "arrival"]117 labels = {"exit": "Exit intensity", "mdfreq": "Markdown frequency", "mddepth": "Markdown depth", "disp": "Residual dispersion (−)", "arrival": "Arrival rate"}118 fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.6))119 ax = axes[0]120 for i, c in enumerate(comps):121 ax.plot(P["week"], P[f"z_{c}"], color=ps.CAT[i], label=labels[c], lw=1.2)122 ax.set_xlabel("Panel week"); ax.set_ylabel("Standardised component (signed)"); ax.set_title("(a) Province-level components"); ax.legend(fontsize=6)123 ax = axes[1]124 ax.plot(P["week"], P["amli_equal"], color=ps.CAT[0], label="Equal weights"); ax.plot(P["week"], P["amli_pca"], color=ps.CAT[1], label="PCA weights")125 if "amli_supervised" in P: ax.plot(P["week"], P["amli_supervised"], color=ps.CAT[2], label="Supervised (ex-exit)")126 ax.set_xlabel("Panel week"); ax.set_ylabel("AMLI (higher = more liquid)"); ax.set_title("(b) Province AMLI"); ax.legend(fontsize=6)127 ax = axes[2]128 top = Z.groupby("region")["ld"].sum().sort_values(ascending=False).index[:6]129 for i, r in enumerate(top):130 z = Z[Z["region"] == r].sort_values("week"); ax.plot(z["week"], z["amli_equal"], color=ps.color_for(ps.REGIONS, r), lw=1.1, label=r)131 if B.exists():132 b = pd.read_csv(B); b = b[b["region"] == r].sort_values("week"); ax.fill_between(b["week"], b["amli_equal_lo"], b["amli_equal_hi"], color=ps.color_for(ps.REGIONS, r), alpha=0.12, linewidth=0)133 ax.set_xlabel("Panel week"); ax.set_ylabel("AMLI (equal weights)"); ax.set_title("(c) Six largest regions, 95 % bootstrap bands"); ax.legend(fontsize=5.5, ncol=2, loc="lower right")134 ps.save(fig, "fig08_amli_series")135136137def fig09_amli_validation():138 Z = pd.read_csv(R / "amli" / "index_region_week.csv")139 others = [f"z_{c}" for c in ["mdfreq", "mddepth", "disp", "arrival"]]140 Z["amli_noexit"] = Z[others].mean(axis=1); Z = Z.sort_values(["region", "week"]); Z["exit_next"] = Z.groupby("region")["exit"].shift(-1)141 d = Z.dropna(subset=["amli_noexit", "exit_next"]); d = d[d["week"] >= 6]142 S = json.load(open(R / "amli" / "summary.json"))["region_validation"]143 fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.7))144 ax = axes[0]145 dm = d.copy(); dm["x"] = dm["amli_noexit"] - dm.groupby("region")["amli_noexit"].transform("mean"); dm["y"] = dm["exit_next"] - dm.groupby("region")["exit_next"].transform("mean")146 ax.scatter(dm["x"], 100 * dm["y"], s=10, color=ps.CAT[0], alpha=0.7)147 b = S["noexit_fe"]["b"]; xs = np.linspace(dm["x"].min(), dm["x"].max(), 10); ax.plot(xs, 100 * b * xs, color=ps.CAT[1], label=f"slope {100*b:.2f} pp per unit [{100*S['noexit_fe']['lo']:.2f}, {100*S['noexit_fe']['hi']:.2f}]")148 ax.set_xlabel("Leave-exit-out AMLI (within-region)"); ax.set_ylabel("Next-week exit intensity (pp, within-region)"); ax.set_title("(a) Held-out component prediction, weeks 6–10"); ax.legend(fontsize=6)149 ax = axes[1]150 corr = Z[["z_exit", "z_mdfreq", "z_mddepth", "z_disp", "z_arrival"]].corr().values151 im = ax.imshow(corr, cmap=ps.cmap_div, vmin=-1, vmax=1); names = ["Exit", "Md. freq.", "Md. depth", "Disp. (−)", "Arrival"]152 ax.set_xticks(range(5)); ax.set_xticklabels(names, rotation=35, ha="right", fontsize=7); ax.set_yticks(range(5)); ax.set_yticklabels(names, fontsize=7)153 for i in range(5):154 for j in range(5):155 ax.text(j, i, f"{corr[i, j]:.2f}", ha="center", va="center", fontsize=6, color=ps.TEXT if abs(corr[i, j]) < 0.6 else "white")156 ax.set_title("(b) Component correlations (region-weeks)"); ax.grid(False); fig.colorbar(im, ax=ax, shrink=0.8)157 ps.save(fig, "fig09_amli_validation")158159160def run(which=None):161 ps.setup()162 fns = {"02": fig02_surface, "05": fig05_hazards, "06": fig06_eventstudy, "07": fig07_region_map, "08": fig08_amli, "09": fig09_amli_validation}163 done, failed = [], {}164 for k, f in fns.items():165 if which and k not in which:166 continue167 try:168 f(); done.append(k)169 except Exception as e: # keep going; missing inputs are reported170 failed[k] = repr(e)171 return {"done": done, "failed": failed}172173174if __name__ == "__main__":175 import sys176 print(run(sys.argv[1:] or None))177