"""Model-based figures: 2 (depreciation surface), 5 (hazards), 6 (event study), 7 (regional premia map), 8 (AMLI series + components), 9 (AMLI internal validation). Descriptive figures 1, 3, 4 live in descriptives.py; figure 10 in prediction.py. Every plotted number is read from results/ files produced by the pipeline.""" from __future__ import annotations import json import matplotlib.pyplot as plt import numpy as np import pandas as pd import polars as pl from . import config as C from . import plotstyle as ps R = C.RESULTS def fig02_surface(): g = pd.read_csv(R / "hedonic" / "surface_A2t.csv") # mask grid cells without data support (fewer than 25 S1 listings within +-0.75 y and +-1.5 (10k km)) L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s1_hedonic")).select(["age_years_entry", "mileage_km"]).to_pandas() a = L["age_years_entry"].values; k = L["mileage_km"].values / 1e4 H, ae, ke = np.histogram2d(a, k, bins=[np.arange(0, 16.5, 0.5), np.arange(0, 31.5, 1.0)]) def support(age, km): i = min(int(age / 0.5), H.shape[0] - 1); j = min(int(km / 1.0), H.shape[1] - 1) return H[max(0, i - 1):i + 2, max(0, j - 1):j + 2].sum() g["support"] = [support(x, y) for x, y in zip(g["age"], g["km"])] g.loc[g["support"] < 25, "logp_rel"] = np.nan A = g.pivot(index="km", columns="age", values="logp_rel") fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.9), gridspec_kw={"width_ratios": [1.25, 1]}) ax = axes[0] X, Y = np.meshgrid(A.columns.values, A.index.values * 10) 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()) 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) ax.clabel(cl, fmt=lambda v: f"{100*(np.exp(v)-1):.0f} %", fontsize=7) 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") 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]]) ax.grid(False) ax = axes[1] sub = g[np.isclose(g["km"], g["km"].iloc[(g["km"] - 8).abs().argmin()])] for km_t, col in [(2, ps.CAT[0]), (8, ps.CAT[1]), (16, ps.CAT[2])]: s = g[np.isclose(g["km"], g["km"].iloc[(g["km"] - km_t).abs().argmin()])].sort_values("age") ax.plot(s["age"], 100 * (np.exp(s["logp_rel"]) - 1), color=col, label=f"{km_t*10000:,} km".replace(",", "\u2009")) 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") ax.legend(title="Mileage", loc="upper right") ps.save(fig, "fig02_depreciation_surface") def fig05_hazards(): hm = pd.read_csv(R / "hazards" / "empirical_hazard_first_markdown.csv"); he = pd.read_csv(R / "hazards" / "empirical_hazard_exit.csv") fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.5)) 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])]: d = d[d["listing_age"] >= 1] ax.fill_between(d["listing_age"], 100 * d["lo"], 100 * d["hi"], color=col, alpha=0.2, linewidth=0) ax.plot(d["listing_age"], 100 * d["hazard"], color=col) ax.set_xlabel("Listing age (days)"); ax.set_ylabel("Daily hazard (%)"); ax.set_title(title); ax.set_xlim(0, 70) ax = axes[2] order = ["direct_dealer", "full_portal", "banner"] for i, sc in enumerate(order): p = R / "hazards" / f"empirical_hazard_exit_{sc}.csv" if p.exists(): d = pd.read_csv(p); d = d[(d["listing_age"] >= 1) & (d["at_risk"] >= 200)] S = np.cumprod(1 - d["hazard"].values) ax.plot(d["listing_age"], 100 * S, color=ps.color_for(ps.CLASS_ORDER, sc), label=ps.CLASS_LABEL[sc]) 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) ps.save(fig, "fig05_hazard_curves") def fig06_eventstudy(): main = pd.read_csv(R / "competition" / "eventstudy_main.csv"); pla = pd.read_csv(R / "competition" / "eventstudy_placebo_permuted.csv") fc = pd.read_csv(R / "competition" / "eventstudy_firstcut.csv") dd = R / "competition" / "eventstudy_dealer_dow_fe.csv" fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.7), sharey=False) ax = axes[0] 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)]: 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) if dd.exists(): d = pd.read_csv(dd); ax.plot(d["k"], 100 * d["Estimate"], color=ps.CAT[1], lw=1.2, label="+ dealer×weekday FE") ax.axhline(0, color=ps.TEXT2, lw=0.6); ax.axvline(-0.5, color=ps.GRID, lw=0.8) 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") ax.legend(fontsize=6.5, loc="upper left") ax = axes[1] 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") 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) ps.save(fig, "fig06_competitor_event_study") def fig07_region_map(): prem = pd.read_csv(R / "hedonic" / "region_premia_A4.csv"); prem = prem[prem["region"] != "(missing)"] prem = pd.concat([prem, pd.DataFrame([{"region": "Montréal", "Estimate": 0.0, "Std. Error": 0.0, "2.5%": 0.0, "97.5%": 0.0}])]) 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)) 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() d = prem.merge(cen, on="region") fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.2), gridspec_kw={"width_ratios": [1.1, 1]}) ax = axes[0] pts = L.select(["lat_r", "lng_r"]).to_pandas() ax.hexbin(pts["lng_r"], pts["lat_r"], gridsize=45, cmap=ps.cmap_seq, mincnt=1, linewidths=0.1, bins="log", alpha=0.55) vmax = max(abs(d["Estimate"].min()), abs(d["Estimate"].max()), 0.03) 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) for _, r in d.iterrows(): ax.annotate(r["region"].split("–")[0][:12], (r["lng"], r["lat"]), fontsize=5, color=ps.TEXT2, xytext=(4, 2), textcoords="offset points") cb = fig.colorbar(sc, ax=ax, shrink=0.8, pad=0.02); cb.set_label("Quality-adjusted premium vs Montréal (%)") ax.set_xlabel("Longitude"); ax.set_ylabel("Latitude"); ax.set_title("(a) Listing locations and regional premia"); ax.grid(False) ax = axes[1] d2 = d.sort_values("Estimate") 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) ax.set_yticks(range(len(d2))); ax.set_yticklabels(d2["region"], fontsize=6.5); ax.axvline(0, color=ps.TEXT2, lw=0.6) ax.set_xlabel("Premium vs Montréal (%), 95 % CI"); ax.set_title("(b) Model A4 region effects") ps.save(fig, "fig07_regional_premia_map") def fig08_amli(): Z = pd.read_csv(R / "amli" / "index_region_week.csv"); P = pd.read_csv(R / "amli" / "index_province_week.csv") B = R / "amli" / "bootstrap_region_week.csv" comps = ["exit", "mdfreq", "mddepth", "disp", "arrival"] labels = {"exit": "Exit intensity", "mdfreq": "Markdown frequency", "mddepth": "Markdown depth", "disp": "Residual dispersion (−)", "arrival": "Arrival rate"} fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.6)) ax = axes[0] for i, c in enumerate(comps): ax.plot(P["week"], P[f"z_{c}"], color=ps.CAT[i], label=labels[c], lw=1.2) ax.set_xlabel("Panel week"); ax.set_ylabel("Standardised component (signed)"); ax.set_title("(a) Province-level components"); ax.legend(fontsize=6) ax = axes[1] 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") if "amli_supervised" in P: ax.plot(P["week"], P["amli_supervised"], color=ps.CAT[2], label="Supervised (ex-exit)") ax.set_xlabel("Panel week"); ax.set_ylabel("AMLI (higher = more liquid)"); ax.set_title("(b) Province AMLI"); ax.legend(fontsize=6) ax = axes[2] top = Z.groupby("region")["ld"].sum().sort_values(ascending=False).index[:6] for i, r in enumerate(top): 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) if B.exists(): 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) 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") ps.save(fig, "fig08_amli_series") def fig09_amli_validation(): Z = pd.read_csv(R / "amli" / "index_region_week.csv") others = [f"z_{c}" for c in ["mdfreq", "mddepth", "disp", "arrival"]] Z["amli_noexit"] = Z[others].mean(axis=1); Z = Z.sort_values(["region", "week"]); Z["exit_next"] = Z.groupby("region")["exit"].shift(-1) d = Z.dropna(subset=["amli_noexit", "exit_next"]); d = d[d["week"] >= 6] S = json.load(open(R / "amli" / "summary.json"))["region_validation"] fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.7)) ax = axes[0] 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") ax.scatter(dm["x"], 100 * dm["y"], s=10, color=ps.CAT[0], alpha=0.7) 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}]") 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) ax = axes[1] corr = Z[["z_exit", "z_mdfreq", "z_mddepth", "z_disp", "z_arrival"]].corr().values im = ax.imshow(corr, cmap=ps.cmap_div, vmin=-1, vmax=1); names = ["Exit", "Md. freq.", "Md. depth", "Disp. (−)", "Arrival"] 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) for i in range(5): for j in range(5): 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") ax.set_title("(b) Component correlations (region-weeks)"); ax.grid(False); fig.colorbar(im, ax=ax, shrink=0.8) ps.save(fig, "fig09_amli_validation") def run(which=None): ps.setup() fns = {"02": fig02_surface, "05": fig05_hazards, "06": fig06_eventstudy, "07": fig07_region_map, "08": fig08_amli, "09": fig09_amli_validation} done, failed = [], {} for k, f in fns.items(): if which and k not in which: continue try: f(); done.append(k) except Exception as e: # keep going; missing inputs are reported failed[k] = repr(e) return {"done": done, "failed": failed} if __name__ == "__main__": import sys print(run(sys.argv[1:] or None))