1"""Model A — hedonic asking prices (H1, H4) on the S1 sample.23log P_i = f(age_i, km_i) + X_i'b + mu_{make x model (x trim)} + tau_{entry week} + theta_{source class} + region + e_i45Identification note (DECISIONS D-010): with an 81-day window, model-year fixed6effects would absorb vehicle age (age varies only with the entry date within a7make-model-year cell); quality is therefore held fixed at the make x model (A2) or8make x model x trim (A4) level and age is identified across model years.9"""10from __future__ import annotations1112import json13import time14from pathlib import Path1516import numpy as np17import pandas as pd18import polars as pl19import pyfixest as pf20from patsy import dmatrix2122from . import config as C2324OUT = C.RESULTS / "hedonic"25AGE_DF, KM_DF = 5, 526CONTROLS = ["transmission", "drivetrain", "powertrain", "body", "source_class"]272829def load() -> pd.DataFrame:30 L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s1_hedonic"))31 df = L.select(["listing_id", "dealer_id", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km",32 "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "age_years_entry",33 "s2_dynamic", "s2l_dynamic_with_stock"]).to_pandas()34 df = df.rename(columns={"body_filled": "body"})35 df["body"] = df["body"].fillna("unknown"); df["trim"] = df["trim"].fillna("(none)")36 df["region"] = df["region"].fillna("(missing)")37 df["logp"] = np.log(df["entry_price"]); df["age"] = df["age_years_entry"]; df["km"] = df["mileage_km"] / 1e438 df["week"] = ((df["entry_ts"] - C.WINDOW_START_TS) // (7 * 86400)).astype(int)39 df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60)40 mm_s = df["make"] + "|" + df["model"]; mmt_s = mm_s + "|" + df["trim"]41 df["mm"] = pd.factorize(mm_s)[0]; df["mmt"] = pd.factorize(mmt_s)[0] # integer FE codes (robust parsing)42 df["mm_label"] = mm_s; df["mmt_label"] = mmt_s43 df["region"] = pd.Categorical(df["region"], categories=["Montréal"] + sorted(r for r in df["region"].unique() if r != "Montréal"))44 # reference categories: ICE, Automatic, AWD/4x4, direct_dealer45 def cat(col, ref):46 levels = [ref] + sorted(x for x in df[col].unique() if x != ref)47 df[col] = pd.Categorical(df[col], categories=levels)48 cat("powertrain", "ICE"); cat("transmission", "Automatic"); cat("drivetrain", "AWD/4x4"); cat("source_class", "direct_dealer")49 return df505152def spline_basis(df: pd.DataFrame, fit_df: pd.DataFrame | None = None) -> tuple[pd.DataFrame, list[str], list[str], list[str]]:53 """Age/km B-spline bases orthogonalised against the linear term (so that the54 'nonlinear' columns test curvature) + tensor product."""55 ref = fit_df if fit_df is not None else df56 qa = np.quantile(ref["age"], [0.2, 0.4, 0.6, 0.8]); qk = np.quantile(ref["km"], [0.2, 0.4, 0.6, 0.8])57 lo_a, hi_a = ref["age"].min(), ref["age"].max(); lo_k, hi_k = ref["km"].min(), ref["km"].max()58 Ba = np.asarray(dmatrix(f"bs(x, knots={list(qa)}, degree=3, include_intercept=False, lower_bound={lo_a}, upper_bound={hi_a}) - 1",59 {"x": np.clip(df["age"], lo_a, hi_a)}))60 Bk = np.asarray(dmatrix(f"bs(x, knots={list(qk)}, degree=3, include_intercept=False, lower_bound={lo_k}, upper_bound={hi_k}) - 1",61 {"x": np.clip(df["km"], lo_k, hi_k)}))62 out = df.copy()63 # orthogonalise each basis column on [1, x] using the fitting sample's projection64 def nonlin(B, x, xref, Bref):65 X = np.column_stack([np.ones(len(xref)), xref]); beta = np.linalg.lstsq(X, Bref, rcond=None)[0]66 return B - np.column_stack([np.ones(len(x)), x]) @ beta67 Ba_ref = Ba if fit_df is None else np.asarray(dmatrix(f"bs(x, knots={list(qa)}, degree=3, include_intercept=False, lower_bound={lo_a}, upper_bound={hi_a}) - 1", {"x": ref["age"]}))68 Bk_ref = Bk if fit_df is None else np.asarray(dmatrix(f"bs(x, knots={list(qk)}, degree=3, include_intercept=False, lower_bound={lo_k}, upper_bound={hi_k}) - 1", {"x": ref["km"]}))69 Na = nonlin(Ba, df["age"].values, ref["age"].values, Ba_ref); Nk = nonlin(Bk, df["km"].values, ref["km"].values, Bk_ref)70 acols = [f"na{j}" for j in range(Na.shape[1])]; kcols = [f"nk{j}" for j in range(Nk.shape[1])]71 out[acols] = Na; out[kcols] = Nk72 tcols = []73 for i in range(Ba.shape[1]):74 for j in range(Bk.shape[1]):75 c = f"t{i}_{j}"; out[c] = Ba[:, i] * Bk[:, j]; tcols.append(c)76 return out, acols, kcols, tcols777879def rhs(acols, kcols, tcols, linear_only=False, tensor=False, powertrain_int=False):80 terms = ["age", "km"]81 if not linear_only:82 terms += acols + kcols83 if tensor:84 terms += tcols85 if powertrain_int:86 terms += [f"i(powertrain, {c}, ref='ICE')" for c in ["age"] + acols]87 terms += [f"C({c})" for c in CONTROLS]88 return " + ".join(terms)899091def fit_all(df: pd.DataFrame) -> dict:92 OUT.mkdir(parents=True, exist_ok=True)93 X, acols, kcols, tcols = spline_basis(df)94 specs = {95 "A1_linear_mm": (rhs(acols, kcols, tcols, linear_only=True) + " + C(region)", "mm + week"),96 "A2_spline_mm": (rhs(acols, kcols, tcols) + " + C(region)", "mm + week"),97 "A2t_tensor_mm": (rhs(acols, kcols, tcols, tensor=True) + " + C(region)", "mm + week"),98 "A3_spline_x_powertrain": (rhs(acols, kcols, tcols, powertrain_int=True) + " + C(region)", "mm + week"),99 "A4_spline_mmt": (rhs(acols, kcols, tcols) + " + C(region)", "mmt + week"),100 "A5_spline_dealerFE": (rhs(acols, kcols, tcols), "mmt + week + dealer_id"),101 }102 res, models = {}, {}103 for name, (r, fe) in specs.items():104 t0 = time.time()105 m = pf.feols(f"logp ~ {r} | {fe}", data=X, vcov={"CRV1": "dealer_id"}, fixef_rm="none")106 models[name] = m107 tab = m.tidy().reset_index(); tab = tab.rename(columns={tab.columns[0]: "term"})108 tab.to_csv(OUT / f"coefs_{name}.csv", index=False)109 res[name] = {"n": int(m._N), "r2": float(m._r2), "r2_within": float(m._r2_within) if m._r2_within is not None else None,110 "rmse": float(np.sqrt(np.mean(m.resid() ** 2))), "seconds": round(time.time() - t0, 1)}111 print(name, res[name], flush=True)112 # --- H1 Wald tests (cluster-robust) -------------------------------------113 tests = {}114 m2 = models["A2_spline_mm"]115 names = list(m2._coefnames)116 def wald(m, cols):117 idx = [names.index(c) for c in cols if c in names]118 b = np.asarray(m.coef())[idx]; V = np.asarray(m._vcov)[np.ix_(idx, idx)]119 W = float(b @ np.linalg.solve(V, b)); k = len(idx)120 from scipy.stats import chi2121 return {"stat": W, "df": k, "p": float(chi2.sf(W, k))}122 tests["H1_nonlinear_age"] = wald(m2, acols)123 tests["H1_nonlinear_km"] = wald(m2, kcols)124 tests["H1_nonlinear_joint"] = wald(m2, acols + kcols)125 m2t = models["A2t_tensor_mm"]; names = list(m2t._coefnames)126 tests["H1_tensor_interaction"] = wald(m2t, tcols)127 m3 = models["A3_spline_x_powertrain"]; names = list(m3._coefnames)128 tests["H1_age_x_powertrain"] = wald(m3, [n for n in names if "powertrain::" in n])129 # --- H4 region premia from A4 (ref Montréal) ----------------------------130 m4 = models["A4_spline_mmt"]; names = list(m4._coefnames)131 reg = [n for n in names if n.startswith("C(region)")]132 tests["H4_region_joint"] = wald(m4, reg)133 t4 = m4.tidy().reset_index(); t4 = t4.rename(columns={t4.columns[0]: "term"})134 rp = t4[t4["term"].str.startswith("C(region)")].copy()135 rp["region"] = rp["term"].str.extract(r"\[T\.(.*)\]")136 rp[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].to_csv(OUT / "region_premia_A4.csv", index=False)137 t2 = models["A2_spline_mm"].tidy().reset_index(); t2 = t2.rename(columns={t2.columns[0]: "term"})138 rp2 = t2[t2["term"].str.startswith("C(region)")].copy(); rp2["region"] = rp2["term"].str.extract(r"\[T\.(.*)\]")139 rp2[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].to_csv(OUT / "region_premia_A2.csv", index=False)140 # residuals (A4) for H5 and dispersion; A2 residuals as alternative141 pd.DataFrame({"listing_id": X["listing_id"], "resid_A4": m4.resid(), "resid_A2": m2.resid(), "region": X["region"].astype(str),142 "mmt": X["mmt_label"], "s2_dynamic": X["s2_dynamic"], "s2l": X["s2l_dynamic_with_stock"]}).to_parquet(OUT / "residuals.parquet")143 # within-specification dispersion by region (IQR of A4 residuals; cells with >= 5 listings)144 # dispersion uses A2 residuals: A4 trim cells with a single listing have residual 0 by construction145 R = pd.DataFrame({"region": X["region"].astype(str), "mmt": X["mmt_label"], "r": m2.resid()})146 disp = (R.groupby("region")["r"].agg(iqr=lambda s: s.quantile(0.75) - s.quantile(0.25), sd="std", n="size")).reset_index()147 disp.to_csv(OUT / "dispersion_by_region.csv", index=False)148 # --- depreciation surface (A2t): prediction grid relative to age=1, km=1 (10k km) ----149 grid = pd.DataFrame([(a, k) for a in np.linspace(0.25, 15, 60) for k in np.linspace(0.5, 30, 60)], columns=["age", "km"])150 base = X.iloc[0:1].copy()151 G = pd.concat([base] * len(grid), ignore_index=True); G["age"] = grid["age"].values; G["km"] = grid["km"].values152 for c in CONTROLS:153 G[c] = pd.Categorical([X[c].mode()[0]] * len(G), categories=X[c].cat.categories) if hasattr(X[c], "cat") else X[c].mode()[0]154 G["region"] = "Montréal"155 Gx, _, _, _ = spline_basis(G, fit_df=df)156 # prediction without FE: use coefficients directly157 coef = models["A2t_tensor_mm"].coef()158 def lin_pred(m, D):159 s = np.zeros(len(D))160 for term, b in m.coef().items():161 if term in D.columns:162 s += b * D[term].values163 return s164 surf = lin_pred(models["A2t_tensor_mm"], Gx)165 grid["logp_rel"] = surf - surf[np.argmin((grid["age"] - 1) ** 2 + (grid["km"] - 1) ** 2)]166 grid.to_csv(OUT / "surface_A2t.csv", index=False)167 # age profile by powertrain (A3): marginal log-price vs age at km = median168 json.dump({"fits": res, "tests": tests, "n_S1": int(len(X)), "age_knots": list(np.quantile(df["age"], [0.2, 0.4, 0.6, 0.8])),169 "km_knots": list(np.quantile(df["km"], [0.2, 0.4, 0.6, 0.8]))}, open(OUT / "summary.json", "w"), indent=2, default=float)170 return {"fits": res, "tests": tests}171172173def grouped_cv(df: pd.DataFrame, k: int = 5, seed: int = 36) -> dict:174 """5-fold CV grouped by dealer for A1, A2, A2t, A4 and LightGBM (A6)."""175 import lightgbm as lgb176 rng = np.random.default_rng(seed)177 dealers = df["dealer_id"].unique(); rng.shuffle(dealers)178 fold_of = {d: i % k for i, d in enumerate(dealers)}179 df = df.assign(fold=df["dealer_id"].map(fold_of))180 out = {m: [] for m in ["A1_linear_mm", "A2_spline_mm", "A2t_tensor_mm", "A4_spline_mmt", "A6_lightgbm"]}181 for f in range(k):182 tr, te = df[df.fold != f], df[df.fold == f]183 Xtr, acols, kcols, tcols = spline_basis(tr); Xte, _, _, _ = spline_basis(te, fit_df=tr)184 for name, (r, fe) in {"A1_linear_mm": (rhs(acols, kcols, tcols, linear_only=True) + " + C(region)", "mm + week"),185 "A2_spline_mm": (rhs(acols, kcols, tcols) + " + C(region)", "mm + week"),186 "A2t_tensor_mm": (rhs(acols, kcols, tcols, tensor=True) + " + C(region)", "mm + week"),187 "A4_spline_mmt": (rhs(acols, kcols, tcols) + " + C(region)", "mmt + week")}.items():188 m = pf.feols(f"logp ~ {r} | {fe}", data=Xtr, vcov="iid", fixef_rm="none")189 pred = np.asarray(m.predict(newdata=Xte), dtype=float)190 ok = ~np.isnan(pred) # unseen FE levels -> NaN; report coverage191 err = Xte["logp"].values[ok] - pred[ok]192 out[name].append({"fold": f, "rmse": float(np.sqrt(np.mean(err ** 2))), "mae": float(np.mean(np.abs(err))), "coverage": float(ok.mean()), "n_test": int(ok.sum())})193 cat = ["make", "model", "trim", "powertrain", "body", "transmission", "drivetrain", "region", "source_class"]194 feats = ["age", "km", "year", "week"] + cat195 Ztr = tr[feats].copy(); Zte = te[feats].copy()196 for c in cat:197 cats = pd.Categorical(Ztr[c]).categories198 Ztr[c] = pd.Categorical(Ztr[c], categories=cats); Zte[c] = pd.Categorical(Zte[c], categories=cats)199 g = lgb.LGBMRegressor(n_estimators=1500, learning_rate=0.03, num_leaves=63, min_child_samples=30, subsample=0.8, subsample_freq=1,200 colsample_bytree=0.8, verbose=-1, random_state=seed)201 g.fit(Ztr, tr["logp"], categorical_feature=cat)202 err = te["logp"].values - g.predict(Zte)203 out["A6_lightgbm"].append({"fold": f, "rmse": float(np.sqrt(np.mean(err ** 2))), "mae": float(np.mean(np.abs(err))), "coverage": 1.0, "n_test": int(len(te))})204 print("fold", f, {m: round(v[-1]["rmse"], 4) for m, v in out.items()}, flush=True)205 summ = {m: {"rmse_mean": float(np.mean([x["rmse"] for x in v])), "rmse_sd": float(np.std([x["rmse"] for x in v])),206 "mae_mean": float(np.mean([x["mae"] for x in v])), "coverage": float(np.mean([x["coverage"] for x in v]))} for m, v in out.items()}207 json.dump({"folds": out, "summary": summ}, open(OUT / "cv_grouped_dealer.json", "w"), indent=2)208 return summ209210211if __name__ == "__main__":212 df = load()213 print(json.dumps(fit_all(df), indent=2, default=float))214 print(json.dumps(grouped_cv(df), indent=2))215