"""Model A — hedonic asking prices (H1, H4) on the S1 sample. log P_i = f(age_i, km_i) + X_i'b + mu_{make x model (x trim)} + tau_{entry week} + theta_{source class} + region + e_i Identification note (DECISIONS D-010): with an 81-day window, model-year fixed effects would absorb vehicle age (age varies only with the entry date within a make-model-year cell); quality is therefore held fixed at the make x model (A2) or make x model x trim (A4) level and age is identified across model years. """ from __future__ import annotations import json import time from pathlib import Path import numpy as np import pandas as pd import polars as pl import pyfixest as pf from patsy import dmatrix from . import config as C OUT = C.RESULTS / "hedonic" AGE_DF, KM_DF = 5, 5 CONTROLS = ["transmission", "drivetrain", "powertrain", "body", "source_class"] def load() -> pd.DataFrame: L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(pl.col("s1_hedonic")) df = L.select(["listing_id", "dealer_id", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km", "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "age_years_entry", "s2_dynamic", "s2l_dynamic_with_stock"]).to_pandas() df = df.rename(columns={"body_filled": "body"}) df["body"] = df["body"].fillna("unknown"); df["trim"] = df["trim"].fillna("(none)") df["region"] = df["region"].fillna("(missing)") df["logp"] = np.log(df["entry_price"]); df["age"] = df["age_years_entry"]; df["km"] = df["mileage_km"] / 1e4 df["week"] = ((df["entry_ts"] - C.WINDOW_START_TS) // (7 * 86400)).astype(int) df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60) mm_s = df["make"] + "|" + df["model"]; mmt_s = mm_s + "|" + df["trim"] df["mm"] = pd.factorize(mm_s)[0]; df["mmt"] = pd.factorize(mmt_s)[0] # integer FE codes (robust parsing) df["mm_label"] = mm_s; df["mmt_label"] = mmt_s df["region"] = pd.Categorical(df["region"], categories=["Montréal"] + sorted(r for r in df["region"].unique() if r != "Montréal")) # reference categories: ICE, Automatic, AWD/4x4, direct_dealer def cat(col, ref): levels = [ref] + sorted(x for x in df[col].unique() if x != ref) df[col] = pd.Categorical(df[col], categories=levels) cat("powertrain", "ICE"); cat("transmission", "Automatic"); cat("drivetrain", "AWD/4x4"); cat("source_class", "direct_dealer") return df def spline_basis(df: pd.DataFrame, fit_df: pd.DataFrame | None = None) -> tuple[pd.DataFrame, list[str], list[str], list[str]]: """Age/km B-spline bases orthogonalised against the linear term (so that the 'nonlinear' columns test curvature) + tensor product.""" ref = fit_df if fit_df is not None else df 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]) lo_a, hi_a = ref["age"].min(), ref["age"].max(); lo_k, hi_k = ref["km"].min(), ref["km"].max() Ba = np.asarray(dmatrix(f"bs(x, knots={list(qa)}, degree=3, include_intercept=False, lower_bound={lo_a}, upper_bound={hi_a}) - 1", {"x": np.clip(df["age"], lo_a, hi_a)})) Bk = np.asarray(dmatrix(f"bs(x, knots={list(qk)}, degree=3, include_intercept=False, lower_bound={lo_k}, upper_bound={hi_k}) - 1", {"x": np.clip(df["km"], lo_k, hi_k)})) out = df.copy() # orthogonalise each basis column on [1, x] using the fitting sample's projection def nonlin(B, x, xref, Bref): X = np.column_stack([np.ones(len(xref)), xref]); beta = np.linalg.lstsq(X, Bref, rcond=None)[0] return B - np.column_stack([np.ones(len(x)), x]) @ beta 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"]})) 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"]})) Na = nonlin(Ba, df["age"].values, ref["age"].values, Ba_ref); Nk = nonlin(Bk, df["km"].values, ref["km"].values, Bk_ref) acols = [f"na{j}" for j in range(Na.shape[1])]; kcols = [f"nk{j}" for j in range(Nk.shape[1])] out[acols] = Na; out[kcols] = Nk tcols = [] for i in range(Ba.shape[1]): for j in range(Bk.shape[1]): c = f"t{i}_{j}"; out[c] = Ba[:, i] * Bk[:, j]; tcols.append(c) return out, acols, kcols, tcols def rhs(acols, kcols, tcols, linear_only=False, tensor=False, powertrain_int=False): terms = ["age", "km"] if not linear_only: terms += acols + kcols if tensor: terms += tcols if powertrain_int: terms += [f"i(powertrain, {c}, ref='ICE')" for c in ["age"] + acols] terms += [f"C({c})" for c in CONTROLS] return " + ".join(terms) def fit_all(df: pd.DataFrame) -> dict: OUT.mkdir(parents=True, exist_ok=True) X, acols, kcols, tcols = spline_basis(df) specs = { "A1_linear_mm": (rhs(acols, kcols, tcols, linear_only=True) + " + C(region)", "mm + week"), "A2_spline_mm": (rhs(acols, kcols, tcols) + " + C(region)", "mm + week"), "A2t_tensor_mm": (rhs(acols, kcols, tcols, tensor=True) + " + C(region)", "mm + week"), "A3_spline_x_powertrain": (rhs(acols, kcols, tcols, powertrain_int=True) + " + C(region)", "mm + week"), "A4_spline_mmt": (rhs(acols, kcols, tcols) + " + C(region)", "mmt + week"), "A5_spline_dealerFE": (rhs(acols, kcols, tcols), "mmt + week + dealer_id"), } res, models = {}, {} for name, (r, fe) in specs.items(): t0 = time.time() m = pf.feols(f"logp ~ {r} | {fe}", data=X, vcov={"CRV1": "dealer_id"}, fixef_rm="none") models[name] = m tab = m.tidy().reset_index(); tab = tab.rename(columns={tab.columns[0]: "term"}) tab.to_csv(OUT / f"coefs_{name}.csv", index=False) 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, "rmse": float(np.sqrt(np.mean(m.resid() ** 2))), "seconds": round(time.time() - t0, 1)} print(name, res[name], flush=True) # --- H1 Wald tests (cluster-robust) ------------------------------------- tests = {} m2 = models["A2_spline_mm"] names = list(m2._coefnames) def wald(m, cols): idx = [names.index(c) for c in cols if c in names] b = np.asarray(m.coef())[idx]; V = np.asarray(m._vcov)[np.ix_(idx, idx)] W = float(b @ np.linalg.solve(V, b)); k = len(idx) from scipy.stats import chi2 return {"stat": W, "df": k, "p": float(chi2.sf(W, k))} tests["H1_nonlinear_age"] = wald(m2, acols) tests["H1_nonlinear_km"] = wald(m2, kcols) tests["H1_nonlinear_joint"] = wald(m2, acols + kcols) m2t = models["A2t_tensor_mm"]; names = list(m2t._coefnames) tests["H1_tensor_interaction"] = wald(m2t, tcols) m3 = models["A3_spline_x_powertrain"]; names = list(m3._coefnames) tests["H1_age_x_powertrain"] = wald(m3, [n for n in names if "powertrain::" in n]) # --- H4 region premia from A4 (ref Montréal) ---------------------------- m4 = models["A4_spline_mmt"]; names = list(m4._coefnames) reg = [n for n in names if n.startswith("C(region)")] tests["H4_region_joint"] = wald(m4, reg) t4 = m4.tidy().reset_index(); t4 = t4.rename(columns={t4.columns[0]: "term"}) rp = t4[t4["term"].str.startswith("C(region)")].copy() rp["region"] = rp["term"].str.extract(r"\[T\.(.*)\]") rp[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].to_csv(OUT / "region_premia_A4.csv", index=False) t2 = models["A2_spline_mm"].tidy().reset_index(); t2 = t2.rename(columns={t2.columns[0]: "term"}) rp2 = t2[t2["term"].str.startswith("C(region)")].copy(); rp2["region"] = rp2["term"].str.extract(r"\[T\.(.*)\]") rp2[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].to_csv(OUT / "region_premia_A2.csv", index=False) # residuals (A4) for H5 and dispersion; A2 residuals as alternative pd.DataFrame({"listing_id": X["listing_id"], "resid_A4": m4.resid(), "resid_A2": m2.resid(), "region": X["region"].astype(str), "mmt": X["mmt_label"], "s2_dynamic": X["s2_dynamic"], "s2l": X["s2l_dynamic_with_stock"]}).to_parquet(OUT / "residuals.parquet") # within-specification dispersion by region (IQR of A4 residuals; cells with >= 5 listings) # dispersion uses A2 residuals: A4 trim cells with a single listing have residual 0 by construction R = pd.DataFrame({"region": X["region"].astype(str), "mmt": X["mmt_label"], "r": m2.resid()}) disp = (R.groupby("region")["r"].agg(iqr=lambda s: s.quantile(0.75) - s.quantile(0.25), sd="std", n="size")).reset_index() disp.to_csv(OUT / "dispersion_by_region.csv", index=False) # --- depreciation surface (A2t): prediction grid relative to age=1, km=1 (10k km) ---- 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"]) base = X.iloc[0:1].copy() G = pd.concat([base] * len(grid), ignore_index=True); G["age"] = grid["age"].values; G["km"] = grid["km"].values for c in CONTROLS: 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] G["region"] = "Montréal" Gx, _, _, _ = spline_basis(G, fit_df=df) # prediction without FE: use coefficients directly coef = models["A2t_tensor_mm"].coef() def lin_pred(m, D): s = np.zeros(len(D)) for term, b in m.coef().items(): if term in D.columns: s += b * D[term].values return s surf = lin_pred(models["A2t_tensor_mm"], Gx) grid["logp_rel"] = surf - surf[np.argmin((grid["age"] - 1) ** 2 + (grid["km"] - 1) ** 2)] grid.to_csv(OUT / "surface_A2t.csv", index=False) # age profile by powertrain (A3): marginal log-price vs age at km = median 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])), "km_knots": list(np.quantile(df["km"], [0.2, 0.4, 0.6, 0.8]))}, open(OUT / "summary.json", "w"), indent=2, default=float) return {"fits": res, "tests": tests} def grouped_cv(df: pd.DataFrame, k: int = 5, seed: int = 36) -> dict: """5-fold CV grouped by dealer for A1, A2, A2t, A4 and LightGBM (A6).""" import lightgbm as lgb rng = np.random.default_rng(seed) dealers = df["dealer_id"].unique(); rng.shuffle(dealers) fold_of = {d: i % k for i, d in enumerate(dealers)} df = df.assign(fold=df["dealer_id"].map(fold_of)) out = {m: [] for m in ["A1_linear_mm", "A2_spline_mm", "A2t_tensor_mm", "A4_spline_mmt", "A6_lightgbm"]} for f in range(k): tr, te = df[df.fold != f], df[df.fold == f] Xtr, acols, kcols, tcols = spline_basis(tr); Xte, _, _, _ = spline_basis(te, fit_df=tr) for name, (r, fe) in {"A1_linear_mm": (rhs(acols, kcols, tcols, linear_only=True) + " + C(region)", "mm + week"), "A2_spline_mm": (rhs(acols, kcols, tcols) + " + C(region)", "mm + week"), "A2t_tensor_mm": (rhs(acols, kcols, tcols, tensor=True) + " + C(region)", "mm + week"), "A4_spline_mmt": (rhs(acols, kcols, tcols) + " + C(region)", "mmt + week")}.items(): m = pf.feols(f"logp ~ {r} | {fe}", data=Xtr, vcov="iid", fixef_rm="none") pred = np.asarray(m.predict(newdata=Xte), dtype=float) ok = ~np.isnan(pred) # unseen FE levels -> NaN; report coverage err = Xte["logp"].values[ok] - pred[ok] 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())}) cat = ["make", "model", "trim", "powertrain", "body", "transmission", "drivetrain", "region", "source_class"] feats = ["age", "km", "year", "week"] + cat Ztr = tr[feats].copy(); Zte = te[feats].copy() for c in cat: cats = pd.Categorical(Ztr[c]).categories Ztr[c] = pd.Categorical(Ztr[c], categories=cats); Zte[c] = pd.Categorical(Zte[c], categories=cats) g = lgb.LGBMRegressor(n_estimators=1500, learning_rate=0.03, num_leaves=63, min_child_samples=30, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, verbose=-1, random_state=seed) g.fit(Ztr, tr["logp"], categorical_feature=cat) err = te["logp"].values - g.predict(Zte) 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))}) print("fold", f, {m: round(v[-1]["rmse"], 4) for m, v in out.items()}, flush=True) summ = {m: {"rmse_mean": float(np.mean([x["rmse"] for x in v])), "rmse_sd": float(np.std([x["rmse"] for x in v])), "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()} json.dump({"folds": out, "summary": summ}, open(OUT / "cv_grouped_dealer.json", "w"), indent=2) return summ if __name__ == "__main__": df = load() print(json.dumps(fit_all(df), indent=2, default=float)) print(json.dumps(grouped_cv(df), indent=2))