1"""Automobile Market Liquidity Index (AMLI) — Phase 5 (analysis_plan §4, DECISIONS D-013).23Higher AMLI = more liquid. Market-week panel (region x ISO-like week of the panel,4weeks 0..10 complete; week 11 partial is dropped). Components computed with5information available in week w, from the reliable-source listing-day universe (S2L):6 exit = confirmed exits / listing-days (+)7 mdfreq = sustained price cuts / listing-days (+)8 mddepth = median sustained cut size (%) (+)9 disp = IQR of hedonic residuals among active listings (-)10 arrival = new entries (S2) / active stock (+)11 stale = share of active stock with listing age > 30 d (-) [5-component variant only: mechanically censored before week 5]12Weights: equal, PCA (first component, signed), supervised ridge on next-week exit (fit weeks 0-5).13"""14from __future__ import annotations1516import json17import time1819import numpy as np20import pandas as pd21import polars as pl2223from . import config as C2425OUT = C.RESULTS / "amli"26WEEKS = 11 # complete weeks 0..10 (days 0..76)27WEEK_MIN = 2 # weeks 0-1 excluded: portal start-up (otogo day 6, autousagee day 11), DECISIONS D-01528MIN_LISTING_DAYS = 3029PRIMARY = ["exit", "mdfreq", "mddepth", "disp", "arrival"]30SIGN = {"exit": 1, "mdfreq": 1, "mddepth": 1, "disp": -1, "arrival": 1, "stale": -1}313233def listing_days() -> pd.DataFrame:34 RS = pl.read_parquet(C.PROCESSED / "risk_set.parquet").select(["listing_id", "t", "listing_age", "cut", "exit", "s2_dynamic", "left_censored", "logp"])35 L = pl.read_parquet(C.PROCESSED / "listings.parquet").select(["listing_id", "dealer_id", "region", "entry_price", "body_filled"])36 R = pl.read_parquet(C.RESULTS / "hedonic" / "residuals.parquet").select(["listing_id", "resid_A2"])37 E = pl.read_parquet(C.PROCESSED / "price_events.parquet").filter(pl.col("kind") == "cut").with_columns(38 ((pl.col("ts") / 86400).floor() - (C.WINDOW_START_TS // 86400)).cast(pl.Int32).alias("t")).select(["listing_id", "t", "change_pct"])39 d = RS.join(L, on="listing_id", how="left").join(R, on="listing_id", how="left").join(E, on=["listing_id", "t"], how="left")40 d = d.with_columns([(pl.col("t") // 7).alias("week"), (pl.col("listing_age") == 0).cast(pl.Int8).alias("entry"),41 (pl.col("listing_age") > 30).cast(pl.Int8).alias("stale_flag"), pl.col("region").fill_null("(missing)")])42 # price tier (terciles of entry price over S2L)43 q = d.select(pl.col("entry_price").quantile(0.333).alias("q1"), pl.col("entry_price").quantile(0.667).alias("q2")).row(0)44 d = d.with_columns(pl.when(pl.col("entry_price") < q[0]).then(pl.lit("low")).when(pl.col("entry_price") < q[1]).then(pl.lit("mid")).otherwise(pl.lit("high")).alias("tier"))45 return d.filter((pl.col("week") < WEEKS) & (pl.col("week") >= WEEK_MIN)).to_pandas()464748def components(d: pd.DataFrame, by: list[str], w: np.ndarray | None = None) -> pd.DataFrame:49 """Aggregate listing-days to market-week components; w = optional row weights (bootstrap)."""50 d = d.copy()51 d["w"] = 1.0 if w is None else w52 d["w_exit"] = d["w"] * d["exit"]; d["w_cut"] = d["w"] * d["cut"]; d["w_entry"] = d["w"] * d["entry"] * d["s2_dynamic"]53 d["w_stale"] = d["w"] * d["stale_flag"]54 g = d.groupby(by + ["week"])55 out = g.agg(ld=("w", "sum"), exits=("w_exit", "sum"), cuts=("w_cut", "sum"), entries=("w_entry", "sum"), stale=("w_stale", "sum")).reset_index()56 out["exit"] = out["exits"] / out["ld"]; out["mdfreq"] = out["cuts"] / out["ld"]; out["arrival"] = out["entries"] / out["ld"]; out["stale"] = out["stale"] / out["ld"]57 # depth: weighted median of cut sizes; dispersion: weighted IQR of residuals (approximate by unweighted when w is None)58 cutd = d[d["cut"] == 1].dropna(subset=["change_pct"])59 if w is None:60 dep = cutd.groupby(by + ["week"])["change_pct"].median().rename("mddepth").reset_index()61 rr = d.dropna(subset=["resid_A2"]).groupby(by + ["week"])["resid_A2"].agg(lambda s: s.quantile(0.75) - s.quantile(0.25)).rename("disp").reset_index()62 else:63 def wq(s, ww, q):64 o = np.argsort(s.values); cs = np.cumsum(ww.values[o]); return s.values[o][np.searchsorted(cs, q * cs[-1])]65 dep = cutd.groupby(by + ["week"]).apply(lambda x: wq(x["change_pct"], x["w"], 0.5)).rename("mddepth").reset_index()66 rr = d.dropna(subset=["resid_A2"]).groupby(by + ["week"]).apply(lambda x: wq(x["resid_A2"], x["w"], 0.75) - wq(x["resid_A2"], x["w"], 0.25)).rename("disp").reset_index()67 out = out.merge(dep, on=by + ["week"], how="left").merge(rr, on=by + ["week"], how="left")68 out["mddepth"] = -out["mddepth"] * 100 # positive % cut69 out.loc[out["ld"] < MIN_LISTING_DAYS, PRIMARY + ["stale"]] = np.nan70 return out717273def build_index(comp: pd.DataFrame, cols: list[str], fit_weeks=range(WEEK_MIN, 6)) -> tuple[pd.DataFrame, dict]:74 Z = comp.copy()75 mu, sd = {}, {}76 for c in cols:77 mu[c] = float(Z[c].mean()); sd[c] = float(Z[c].std()); Z[f"z_{c}"] = SIGN[c] * (Z[c] - mu[c]) / sd[c]78 zc = [f"z_{c}" for c in cols]79 ok = Z[zc].notna().sum(axis=1) >= min(3, len(cols))80 Z["amli_equal"] = Z.loc[ok, zc].mean(axis=1)81 # PCA on complete rows82 comp_rows = Z[zc].dropna()83 X = comp_rows.values - comp_rows.values.mean(0)84 U, S, Vt = np.linalg.svd(X, full_matrices=False)85 v = Vt[0]; v = v * np.sign(v[cols.index("exit")] if "exit" in cols else v[0])86 expl = float(S[0] ** 2 / (S ** 2).sum())87 Z["amli_pca"] = (Z[zc].fillna(Z[zc].mean()) @ v) / np.abs(v).sum()88 Z.loc[~ok, "amli_pca"] = np.nan89 # supervised: ridge of next-week exit on the other components (fit weeks)90 weights = {"equal": {c: 1 / len(cols) for c in cols}, "pca": dict(zip(cols, v.tolist())), "pca_explained": expl}91 if "exit" in cols:92 others = [f"z_{c}" for c in cols if c != "exit"]93 Z = Z.sort_values([*[c for c in Z.columns if c in ("region", "tier", "body_filled", "market")], "week"])94 key = [c for c in ("region", "tier", "body_filled", "market") if c in Z.columns]95 Z["exit_next"] = Z.groupby(key)["exit"].shift(-1) if key else Z["exit"].shift(-1)96 fitrows = Z[Z["week"].isin(list(fit_weeks))].dropna(subset=others + ["exit_next"])97 A = fitrows[others].values; y = fitrows["exit_next"].values; lam = 1.098 beta = np.linalg.solve(A.T @ A + lam * np.eye(A.shape[1]), A.T @ (y - y.mean()))99 Z["amli_supervised"] = Z[others].fillna(0) @ beta / np.abs(beta).sum()100 weights["supervised"] = dict(zip([c for c in cols if c != "exit"], beta.tolist()))101 weights["supervised_fit_weeks"] = list(fit_weeks)102 return Z, {"mu": mu, "sd": sd, "weights": weights}103104105def predictive_validation(Z: pd.DataFrame, key: list[str], eval_weeks=range(6, WEEKS)) -> dict:106 """Leave-exit-out index (equal weights on non-exit components) vs next-week exit intensity, with market FE."""107 import pyfixest as pf108 others = [f"z_{c}" for c in PRIMARY if c != "exit"]109 Z = Z.copy(); Z["amli_noexit"] = Z[others].mean(axis=1)110 Z["mkt"] = Z[key].astype(str).agg("|".join, axis=1) if key else "all"111 Z["exit_next"] = Z.groupby("mkt")["exit"].shift(-1)112 d = Z[Z["week"].isin(list(eval_weeks))].dropna(subset=["amli_noexit", "exit_next", "exit"])113 res = {}114 for name, f in {"noexit_fe": "exit_next ~ amli_noexit | mkt", "noexit_fe_lagexit": "exit_next ~ amli_noexit + exit | mkt", "noexit_pooled": "exit_next ~ amli_noexit"}.items():115 m = pf.feols(f, data=d, vcov={"CRV1": "mkt"} if "|" in f or len(d["mkt"].unique()) > 1 else "hetero")116 t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); r = t[t["term"] == "amli_noexit"].iloc[0]117 res[name] = {"b": float(r["Estimate"]), "se": float(r["Std. Error"]), "lo": float(r["2.5%"]), "hi": float(r["97.5%"]), "n": int(m._N), "r2": float(m._r2)}118 return res119120121def bootstrap(d: pd.DataFrame, key: list[str], reps: int = 200, seed: int = 36) -> pd.DataFrame:122 """Block bootstrap over dealers: index values per market-week (equal weights, 5 components)."""123 rng = np.random.default_rng(seed)124 dealers = d["dealer_id"].unique(); idx = pd.factorize(d["dealer_id"])[0]125 rows = []126 for r in range(reps):127 draw = rng.integers(0, len(dealers), len(dealers)); mult = np.bincount(draw, minlength=len(dealers)).astype(float)128 w = mult[idx]129 comp = components(d[w > 0], key, w[w > 0])130 Z, _ = build_index(comp, PRIMARY)131 Z["rep"] = r; rows.append(Z[key + ["week", "amli_equal", "amli_pca", "exit"]])132 if r % 20 == 0:133 print("boot", r, flush=True)134 B = pd.concat(rows)135 q = B.groupby(key + ["week"]).agg(amli_equal_lo=("amli_equal", lambda s: s.quantile(0.025)), amli_equal_hi=("amli_equal", lambda s: s.quantile(0.975)),136 amli_pca_lo=("amli_pca", lambda s: s.quantile(0.025)), amli_pca_hi=("amli_pca", lambda s: s.quantile(0.975))).reset_index()137 return q138139140def run(reps: int = 200) -> dict:141 OUT.mkdir(parents=True, exist_ok=True); t0 = time.time()142 d = listing_days()143 out = {"listing_days": int(len(d)), "weeks": WEEKS}144 # --- region x week (primary) --------------------------------------------145 comp = components(d, ["region"])146 Z, meta = build_index(comp, PRIMARY)147 Z.to_csv(OUT / "index_region_week.csv", index=False); out["region"] = meta148 out["region_validation"] = predictive_validation(Z, ["region"])149 # 5-component variant with staleness (weeks >= 6)150 Z5, meta5 = build_index(comp[comp["week"] >= 6], PRIMARY + ["stale"]); Z5.to_csv(OUT / "index_region_week_5comp.csv", index=False); out["region_5comp"] = meta5151 # leave-one-component-out152 loo = {}153 for c in PRIMARY:154 cols = [x for x in PRIMARY if x != c]; Zc, mc = build_index(comp, cols)155 loo[c] = {"corr_with_full_equal": float(Zc["amli_equal"].corr(Z["amli_equal"])), "pca_explained": mc["weights"]["pca_explained"]}156 out["leave_one_out"] = loo157 out["corr_equal_pca"] = float(Z["amli_equal"].corr(Z["amli_pca"])); out["corr_equal_supervised"] = float(Z["amli_equal"].corr(Z["amli_supervised"]))158 # --- province x week, daily province -------------------------------------159 d["market"] = "Quebec"; compP = components(d, ["market"]); ZP, metaP = build_index(compP, PRIMARY); ZP.to_csv(OUT / "index_province_week.csv", index=False); out["province"] = metaP160 # --- segments -----------------------------------------------------------161 for seg in ["tier", "body_filled"]:162 dd = d.dropna(subset=[seg]); cs = components(dd, [seg]); Zs, ms = build_index(cs, PRIMARY); Zs.to_csv(OUT / f"index_{seg}_week.csv", index=False)163 out[f"segment_{seg}"] = {"pca_explained": ms["weights"]["pca_explained"], "validation": predictive_validation(Zs, [seg])}164 # --- component correlations ------------------------------------------------165 out["component_corr"] = comp[PRIMARY].corr().round(3).to_dict()166 # --- bootstrap --------------------------------------------------------------167 if reps:168 q = bootstrap(d, ["region"], reps=reps); q.to_csv(OUT / "bootstrap_region_week.csv", index=False)169 out["seconds"] = round(time.time() - t0, 1)170 json.dump(out, open(OUT / "summary.json", "w"), indent=2, default=float)171 return out172173174if __name__ == "__main__":175 print(json.dumps(run(), indent=2, default=float))176