#!/usr/bin/env python3 """ 10_revision_reanalysis.py -- Editorial-board revision of UQO WP 4 (October 2026). Re-estimates every result reported in the revised paper from the only locally available copy of the Quebec Airbnb scrape, ``data/processed/airbnb_clean_wp5.parquet`` (3,456 priced listings out of 5,000 scraped; prices winsorised at p1/p99 and ratings median-imputed by the companion paper's cleaning step -- the imputation is reversed here by parsing the star rating from the structured title). Outputs (all regenerated, never typed): results/tables/rev_*.tex / .csv / .json paper/figures/rev_*.pdf (and results/figures/) Run time: about 3-4 minutes on an Apple-silicon laptop. """ from __future__ import annotations import json import re import sys import time import warnings from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm from scipy import stats warnings.filterwarnings("ignore") ROOT = Path(__file__).resolve().parent.parent DATA = ROOT / "data" / "processed" / "airbnb_clean_wp5.parquet" TAB = ROOT / "results" / "tables" FIG_R = ROOT / "results" / "figures" FIG_P = ROOT / "paper" / "figures" for d in (TAB, FIG_R, FIG_P): d.mkdir(parents=True, exist_ok=True) import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.family": "serif", "font.size": 10, "axes.spines.top": False, "axes.spines.right": False, "figure.dpi": 150}) SEED = 42 NIGHTLY_CAP = 1000.0 # main analysis sample: price <= $1,000 KEY = {} # numbers quoted in the text T0 = time.perf_counter() def log(msg): print(f"[{time.perf_counter() - T0:6.1f}s] {msg}", flush=True) def savefig(fig, stem): for d in (FIG_R, FIG_P): fig.savefig(d / f"{stem}.pdf", bbox_inches="tight") plt.close(fig) def stars(p): return "***" if p < 0.01 else "**" if p < 0.05 else "*" if p < 0.10 else "" def fmt(x, d=3): return f"{x:.{d}f}" def tex_escape(s): return s.replace("&", r"\&").replace("%", r"\%").replace("_", r"\_").replace("$", r"\$") def write_tex(name, body): (TAB / f"{name}.tex").write_text(body, encoding="utf-8") log(f"wrote {name}.tex") # -------------------------------------------------------------------------- # 1. Load and parse the structured title # -------------------------------------------------------------------------- df = pd.read_parquet(DATA) df["name"] = df["name"].fillna("").astype(str) N_PRICED = len(df) nm = df["name"] # star rating from title (reverses the median imputation of the parquet) star = nm.str.extract(r"★\s*(\d+(?:[.,]\d+)?)", expand=False).str.replace(",", ".", regex=False).astype(float) df["rating_raw"] = star df["rated"] = star.notna().astype(int) df["is_new"] = nm.str.contains(r"★\s*New\b|·\s*New\s*·", regex=True).astype(int) # 37 German-locale titles ("★4,97") were truncated to 4.0 in the scraped rating field; the title value is retained KEY["n_rating_truncated_fixed"] = int(((star.notna()) & ((star - df["rating"]).abs() > 0.011)).sum()) KEY["n_unrated"] = int((df["rated"] == 0).sum()) KEY["n_new"] = int(df["is_new"].sum()) # title category (English template "Category in City · ..."; 50 German-locale titles fall back to property_type) cat = nm.str.extract(r"^(.+?)\s+in\s+", expand=False) CAT_MAP = { "Rental unit": "Rental unit", "Chalet": "Chalet/cottage/cabin", "Cottage": "Chalet/cottage/cabin", "Cabin": "Chalet/cottage/cabin", "Vacation home": "Chalet/cottage/cabin", "Holiday home": "Chalet/cottage/cabin", "Nature lodge": "Chalet/cottage/cabin", "Home": "House", "Townhouse": "House", "Bungalow": "House", "Villa": "House", "Farm stay": "House", "Earthen home": "House", "Earth homes": "House", "Castle": "House", "Condo": "Condo/loft/apartment", "Loft": "Condo/loft/apartment", "Serviced apartment": "Condo/loft/apartment", "Apartment": "Condo/loft/apartment", "Guest suite": "Condo/loft/apartment", "Hotel": "Hotel/B&B/hostel", "Bed and breakfast": "Hotel/B&B/hostel", "Aparthotel": "Hotel/B&B/hostel", "Boutique hotel": "Hotel/B&B/hostel", "Hostel": "Hotel/B&B/hostel", "Guesthouse": "Hotel/B&B/hostel", "Guest house": "Hotel/B&B/hostel", "Resort": "Hotel/B&B/hostel", "casa particular": "Hotel/B&B/hostel", } PT_FALLBACK = {"Rental unit": "Rental unit", "Cabin/Chalet": "Chalet/cottage/cabin", "House": "House", "Condo": "Condo/loft/apartment", "Apartment": "Condo/loft/apartment", "Other": "Other/unconventional"} df["category"] = cat.map(CAT_MAP) df.loc[df["category"].isna() & cat.notna(), "category"] = "Other/unconventional" df.loc[cat.isna(), "category"] = df.loc[cat.isna(), "property_type"].map(PT_FALLBACK).fillna("Other/unconventional") CATS = ["Rental unit", "Chalet/cottage/cabin", "House", "Condo/loft/apartment", "Hotel/B&B/hostel", "Other/unconventional"] df["category"] = pd.Categorical(df["category"], categories=CATS) df["shared_bath"] = nm.str.contains(r"shared\s+(?:half-?)?bath|Gemeinschaftsbad", case=False, regex=True).astype(int) df["studio"] = nm.str.contains(r"\bStudio\b", case=False, regex=True).astype(int) df["beds"] = nm.str.extract(r"(\d+)\s*(?:beds?|Bett(?:en)?)\b", expand=False).astype(float) # the scraper stored bathrooms = 0 whenever the title qualified the bathroom ("1 shared bath", "1 private bath", # "Half-bath"); recover the count from the title; studios have bedrooms = 0 and are recoded to one sleeping room baths_title = nm.str.extract(r"(\d+(?:[.,]\d+)?)\s*(?:shared |private )?(?:bath|Badezimmer)", expand=False).str.replace(",", ".", regex=False).astype(float) half_bath = nm.str.contains(r"half-?bath", case=False, regex=True) df["bathrooms_raw"] = df["bathrooms"] fix = (df["bathrooms"] == 0) & baths_title.notna() df.loc[fix, "bathrooms"] = baths_title[fix] df.loc[(df["bathrooms"] == 0) & half_bath, "bathrooms"] = 0.5 KEY["n_bathrooms_recovered"] = int(fix.sum() + ((df["bathrooms_raw"] == 0) & half_bath & ~fix).sum()) bedr_title = nm.str.extract(r"(\d+)\s*(?:bedrooms?|Schlafzimmer)", expand=False).astype(float) fixb = (df["bedrooms"] == 0) & bedr_title.notna() & (df["studio"] == 0) df.loc[fixb, "bedrooms"] = bedr_title[fixb] df["bedrooms_raw"] = df["bedrooms"] df.loc[df["studio"] == 1, "bedrooms"] = 1 KEY["n_bedrooms_recovered"] = int(fixb.sum()); KEY["n_studio_recoded"] = int(df["studio"].sum()) # numeric prep for c in ("is_superhost", "is_guest_favorite", "pets_allowed"): df[c] = df[c].astype(int) df["price"] = df["price_numeric"].astype(float) df["lnp"] = np.log(df["price"]) df["ln_bedrooms"] = np.log(df["bedrooms"].clip(lower=1)) df["ln_bathrooms"] = np.log1p(df["bathrooms"]) df["ln_guests"] = np.log(df["guests_count"].clip(lower=1)) df["ln_amenities"] = np.log1p(df["amenities_count"]) df["ln_images"] = np.log1p(df["num_images"]) df["ln_reviews"] = np.log1p(df["num_reviews"]) df["rating_c"] = np.where(df["rated"] == 1, df["rating_raw"] - 4.8, 0.0) # centred at 4.8, zero when unrated df["is_montreal"] = (df["city"] == "Montreal").astype(int) # geography: projected km, distance to nearest centre, local density R = 6371.0 lat0, lon0 = df["lat"].mean(), df["long"].mean() df["x_km"] = R * np.radians(df["long"] - lon0) * np.cos(np.radians(lat0)) df["y_km"] = R * np.radians(df["lat"] - lat0) sys.path.insert(0, str(ROOT)) from src.config import CITY_CENTERS # noqa: E402 centres = {v.name: (v.lat, v.lon) for v in CITY_CENTERS.values()} cx = np.array([R * np.radians(lo - lon0) * np.cos(np.radians(lat0)) for la, lo in centres.values()]) cy = np.array([R * np.radians(la - lat0) for la, lo in centres.values()]) dd = np.sqrt((df["x_km"].values[:, None] - cx[None, :]) ** 2 + (df["y_km"].values[:, None] - cy[None, :]) ** 2) df["dist_centre_km"] = dd.min(axis=1) df["ln_dist"] = np.log1p(df["dist_centre_km"]) from sklearn.neighbors import KDTree # noqa: E402 tree = KDTree(df[["x_km", "y_km"]].values) df["density_1km"] = tree.query_radius(df[["x_km", "y_km"]].values, r=1.0, count_only=True) - 1 df["density_500m"] = tree.query_radius(df[["x_km", "y_km"]].values, r=0.5, count_only=True) - 1 df["ln_density_1km"] = np.log1p(df["density_1km"]) KEY["n_centres"] = len(centres) # -------------------------------------------------------------------------- # 2. The price field: bimodality, mixture, regime model # -------------------------------------------------------------------------- from sklearn.mixture import GaussianMixture # noqa: E402 gm = GaussianMixture(2, random_state=SEED, n_init=5).fit(df[["lnp"]].values) order = np.argsort(gm.means_.ravel()) mu = gm.means_.ravel()[order]; sd = np.sqrt(gm.covariances_.ravel()[order]); wt = gm.weights_[order] grid = np.linspace(df["lnp"].min(), df["lnp"].max(), 4000) post_hi = gm.predict_proba(grid[:, None])[:, order[1]] thr_ln = grid[np.argmax(post_hi >= 0.5)] KEY["mix"] = {"mu_lo": float(np.exp(mu[0])), "mu_hi": float(np.exp(mu[1])), "sd_lo": float(sd[0]), "sd_hi": float(sd[1]), "w_lo": float(wt[0]), "w_hi": float(wt[1]), "threshold": float(np.exp(thr_ln))} KEY["share_gt1000"] = float((df["price"] > 1000).mean()) KEY["n_gt1000"] = int((df["price"] > 1000).sum()) KEY["share_gt_mix"] = float((df["price"] > np.exp(thr_ln)).mean()) KEY["share_le600"] = float((df["price"] <= 600).mean()) KEY["price_median_all"] = float(df["price"].median()); KEY["price_mean_all"] = float(df["price"].mean()) KEY["bimodality_coefficient"] = float((stats.skew(df["lnp"]) ** 2 + 1) / (stats.kurtosis(df["lnp"]) + 3 * (N_PRICED - 1) ** 2 / ((N_PRICED - 2) * (N_PRICED - 3)))) KEY["hartigan_note"] = "dip test not run (package not installed)" # ratio of medians by bedrooms lo, hi = df[df.price <= 1000], df[df.price > 1000] ratio = (hi.groupby(hi.bedrooms.clip(upper=5))["price"].median() / lo.groupby(lo.bedrooms.clip(upper=5))["price"].median()) KEY["median_ratio_by_bedrooms"] = {str(k): float(v) for k, v in ratio.items()} rows = [] for b in range(6): l = lo[lo.bedrooms.clip(upper=5) == b]; h = hi[hi.bedrooms.clip(upper=5) == b] rows.append((f"{b}" if b < 5 else "5+", len(l), l.price.median(), len(h), h.price.median(), h.price.median() / l.price.median())) body = ["\\begin{tabular}{lrrrrr}", "\\toprule", "Bedrooms & $N$ & Median (\\$) & $N$ & Median (\\$) & Ratio \\\\", " & \\multicolumn{2}{c}{$\\leq$ \\$1{,}000} & \\multicolumn{2}{c}{$>$ \\$1{,}000} & \\\\", "\\midrule"] for r in rows: body.append(f"{r[0]} & {r[1]:,} & {r[2]:,.0f} & {r[3]:,} & {r[4]:,.0f} & {r[5]:.1f} \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_price_regimes", "\n".join(body)) # ---- two-regime mixture of hedonic regressions (EM) --------------------------------- # ln p_i = x_i'beta + delta * z_i + e_i, z_i in {0 (nightly), 1 (stay)}, e_i | z_i ~ N(0, sigma_z^2), # P(z_i = 1 | w_i) = Lambda(w_i'gamma). Common slopes; the stay regime shifts the intercept by delta # (a log stay length) and has its own variance (stay lengths vary). _CATD = [c for c in CATS if c != "Rental unit"] def _X_mix(d): X = pd.DataFrame({"ln_bedrooms": d["ln_bedrooms"], "ln_bathrooms": d["ln_bathrooms"], "ln_guests": d["ln_guests"], "shared_bath": d["shared_bath"], "studio": d["studio"], "ln_amenities": d["ln_amenities"], "ln_images": d["ln_images"], "is_superhost": d["is_superhost"], "is_guest_favorite": d["is_guest_favorite"], "ln_reviews": d["ln_reviews"], "rated": d["rated"], "rating_c": d["rating_c"], "ln_dist": d["ln_dist"], "ln_density_1km": d["ln_density_1km"], "is_montreal": d["is_montreal"]}, index=d.index) for c in _CATD: X[c] = (d["category"] == c).astype(float) return sm.add_constant(X.astype(float)) Xm = _X_mix(df); ym = df["lnp"].values Wm = sm.add_constant(pd.DataFrame({"is_montreal": df["is_montreal"], "ln_reviews": df["ln_reviews"], "ln_guests": df["ln_guests"], "is_superhost": df["is_superhost"], "hotel": (df["category"] == "Hotel/B&B/hostel").astype(float), "chalet": (df["category"] == "Chalet/cottage/cabin").astype(float)}, index=df.index).astype(float)) r = (df["price"] > 1000).astype(float).values # initialisation: the $1,000 rule delta, s0, s1 = 3.0, 0.5, 0.8; gamma = None for it in range(300): # M-step: WLS on stacked data for (beta, delta) Xs = np.vstack([np.column_stack([Xm.values, np.zeros(len(ym))]), np.column_stack([Xm.values, np.ones(len(ym))])]) ys = np.concatenate([ym, ym]); ws = np.concatenate([(1 - r) / s0 ** 2, r / s1 ** 2]) coef = np.linalg.lstsq(Xs * np.sqrt(ws)[:, None], ys * np.sqrt(ws), rcond=None)[0] beta, delta = coef[:-1], coef[-1] f0 = Xm.values @ beta; f1 = f0 + delta s0 = np.sqrt(np.sum((1 - r) * (ym - f0) ** 2) / np.sum(1 - r)); s1 = np.sqrt(np.sum(r * (ym - f1) ** 2) / np.sum(r)) glm = sm.GLM(r, Wm, family=sm.families.Binomial()).fit(); gamma = glm.params.values; pi = glm.fittedvalues.values # E-step l0 = (1 - pi) * stats.norm.pdf(ym, f0, s0); l1 = pi * stats.norm.pdf(ym, f1, s1) r_new = l1 / (l0 + l1) ll = np.sum(np.log(l0 + l1)) if np.max(np.abs(r_new - r)) < 1e-7: r = r_new; break r = r_new df["p_stay"] = r df["stay_regime"] = (r >= 0.5).astype(int) KEY["mixreg"] = {"iterations": it + 1, "loglik": float(ll), "delta": float(delta), "implied_nights": float(np.exp(delta)), "sigma_nightly": float(s0), "sigma_stay": float(s1), "share_stay": float(r.mean()), "n_stay": int(df["stay_regime"].sum()), "n_nightly": int((df["stay_regime"] == 0).sum()), "gamma": {k: [float(v), float(glm.bse[k]), float(glm.pvalues[k])] for k, v in glm.params.items()}, "beta_common": {k: float(v) for k, v in zip(Xm.columns, beta)}, "n_stay_le1000": int(((df["stay_regime"] == 1) & (df["price"] <= 1000)).sum()), "n_nightly_gt1000": int(((df["stay_regime"] == 0) & (df["price"] > 1000)).sum()), "n_ambiguous_25_75": int(((r > 0.25) & (r < 0.75)).sum()), "entropy_R2": float(1 + np.sum(r * np.log(np.clip(r, 1e-12, 1)) + (1 - r) * np.log(np.clip(1 - r, 1e-12, 1))) / (len(r) * np.log(2)))} # table: mixture parameters body = ["\\begin{tabular}{lr}", "\\toprule", "Parameter & Estimate \\\\", "\\midrule", f"Intercept shift of the stay regime, $\\delta$ (log points) & {delta:.3f} \\\\", f"Implied stay length, $\\exp(\\delta)$ (nights) & {np.exp(delta):.1f} \\\\", f"Residual s.d., nightly regime, $\\sigma_0$ & {s0:.3f} \\\\", f"Residual s.d., stay regime, $\\sigma_1$ & {s1:.3f} \\\\", f"Mean posterior probability of the stay regime & {r.mean():.3f} \\\\", f"Listings classified in the stay regime ($\\hat{{p}}_i \\geq 0.5$) & {int(df['stay_regime'].sum()):,} \\\\", f"\\quad of which displayed price $\\leq$ \\$1{{,}}000 & {KEY['mixreg']['n_stay_le1000']:,} \\\\", f"Listings classified nightly with displayed price $>$ \\$1{{,}}000 & {KEY['mixreg']['n_nightly_gt1000']:,} \\\\", f"Listings with $0.25 < \\hat{{p}}_i < 0.75$ & {KEY['mixreg']['n_ambiguous_25_75']:,} \\\\", f"Log-likelihood & {ll:,.1f} \\\\", f"EM iterations & {it + 1} \\\\", "\\midrule", "\\multicolumn{2}{l}{\\textit{Regime equation (logit coefficients, SE)}} \\\\"] GL = {"const": "Constant", "is_montreal": "Montreal", "ln_reviews": "ln(1 + reviews)", "ln_guests": "ln(guests)", "is_superhost": "Superhost", "hotel": "Hotel/B\\&B/hostel", "chalet": "Chalet/cottage/cabin"} for k, v in glm.params.items(): body.append(f"{GL[k]} & {v:.3f}{stars(glm.pvalues[k])} ({glm.bse[k]:.3f}) \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_mixture", "\n".join(body)) # regime logit: which listings are classified in the stay regime? Xl = pd.get_dummies(df[["ln_bedrooms", "ln_bathrooms", "ln_guests", "ln_amenities", "is_superhost", "is_guest_favorite", "rated", "ln_reviews", "is_montreal", "ln_dist", "category"]], columns=["category"], drop_first=True, dtype=float) Xl.columns = [c.replace("category_", "") for c in Xl.columns] logit = sm.Logit(df["stay_regime"], sm.add_constant(Xl)).fit(disp=0, cov_type="HC1") me = logit.get_margeff(at="overall") LAB = {"ln_bedrooms": "ln(bedrooms)", "ln_bathrooms": "ln(1 + bathrooms)", "ln_guests": "ln(guests)", "ln_amenities": "ln(1 + amenities)", "is_superhost": "Superhost", "is_guest_favorite": "Guest Favourite", "rated": "Rated (star shown)", "ln_reviews": "ln(1 + reviews)", "is_montreal": "Montreal", "ln_dist": "ln(1 + km to centre)", "rating_c": "Rating $-$ 4.8 (rated only)", "ln_images": "ln(1 + photos)", "pets_allowed": "Pets allowed", "shared_bath": "Shared bathroom", "studio": "Studio", "ln_density_1km": "ln(1 + listings within 1 km)", "Chalet/cottage/cabin": "Chalet/cottage/cabin", "House": "House", "Condo/loft/apartment": "Condo/loft/apartment", "Hotel/B&B/hostel": "Hotel/B\\&B/hostel", "Other/unconventional": "Other/unconventional", "const": "Constant"} body = ["\\begin{tabular}{lrrrr}", "\\toprule", "Variable & Coefficient & (SE) & Avg.\\ marginal effect & (SE) \\\\", "\\midrule"] me_tab = me.summary_frame() for v in Xl.columns: body.append(f"{LAB.get(v, v)} & {logit.params[v]:.3f}{stars(logit.pvalues[v])} & ({logit.bse[v]:.3f}) & " f"{me_tab.loc[v, 'dy/dx']:.3f} & ({me_tab.loc[v, 'Std. Err.']:.3f}) \\\\") body += ["\\midrule", f"Observations & \\multicolumn{{4}}{{c}}{{{int(logit.nobs):,}}} \\\\", f"Pseudo-$R^2$ & \\multicolumn{{4}}{{c}}{{{logit.prsquared:.3f}}} \\\\", f"Share in regime & \\multicolumn{{4}}{{c}}{{{df['stay_regime'].mean():.3f}}} \\\\", "\\bottomrule", "\\end{tabular}"] write_tex("rev_stay_logit", "\n".join(body)) KEY["logit"] = {v: [float(logit.params[v]), float(logit.bse[v]), float(logit.pvalues[v]), float(me_tab.loc[v, "dy/dx"])] for v in Xl.columns} # figure: log-price histogram with mixture fig, ax = plt.subplots(figsize=(6.5, 3.8)) ax.hist(df["lnp"], bins=60, density=True, color="#9ecae1", edgecolor="white") xs = grid for k, c in zip(range(2), ("#3182bd", "#de2d26")): ax.plot(xs, wt[k] * stats.norm.pdf(xs, mu[k], sd[k]), color=c, lw=1.5, label=f"component {k + 1}: mode \\${np.exp(mu[k]):,.0f}, weight {wt[k]:.2f}") ax.plot(xs, gm.score_samples(xs[:, None]).pipe(np.exp) if False else np.exp(gm.score_samples(xs[:, None])), "k--", lw=1, label="mixture") ax.axvline(np.log(1000), color="grey", ls=":", lw=1); ax.text(np.log(1000) + 0.05, ax.get_ylim()[1] * 0.9, "\\$1,000", color="grey") ax.axvline(thr_ln, color="#de2d26", ls=":", lw=1) ticks = [25, 50, 100, 200, 500, 1000, 2000, 5000, 10000] ax.set_xticks(np.log(ticks)); ax.set_xticklabels([f"{t:,}" for t in ticks]) ax.set_xlabel("Displayed price (CAD, log scale)"); ax.set_ylabel("Density"); ax.legend(frameon=False, fontsize=8) savefig(fig, "rev_fig_logprice_mixture") # -------------------------------------------------------------------------- # 3. Analysis sample and descriptives # -------------------------------------------------------------------------- A = df[df["stay_regime"] == 0].copy().reset_index(drop=True) # nightly regime (posterior < 0.5) KEY["N_A"] = len(A); KEY["N_priced"] = N_PRICED; KEY["N_scraped"] = 5000 KEY["A_price_median"] = float(A.price.median()); KEY["A_price_mean"] = float(A.price.mean()) KEY["A_share_montreal"] = float(A.is_montreal.mean()); KEY["A_n_cities"] = int(A.city.nunique()) KEY["A_unrated_share"] = float((A.rated == 0).mean()); KEY["A_superhost"] = float(A.is_superhost.mean()) KEY["A_guestfav"] = float(A.is_guest_favorite.mean()); KEY["A_cat_shares"] = A.category.value_counts(normalize=True).round(3).to_dict() KEY["A_cat_shares"] = {str(k): float(v) for k, v in KEY["A_cat_shares"].items()} KEY["A_p1_clip_n"] = int((A.price <= df.price.min() + 1e-6).sum()); KEY["A_price_max"] = float(A.price.max()); KEY["A_share_gt1000"] = float((A.price > 1000).mean()) DESC = [("price", "Displayed price (CAD)"), ("lnp", "ln(price)"), ("bedrooms", "Bedrooms"), ("bathrooms", "Bathrooms"), ("guests_count", "Guest capacity"), ("beds", "Beds (from title)"), ("amenities_count", "Amenities count"), ("num_images", "Photos"), ("num_reviews", "Reviews"), ("rating_raw", "Rating (rated listings)"), ("rated", "Rated (star shown)"), ("is_superhost", "Superhost"), ("is_guest_favorite", "Guest Favourite"), ("pets_allowed", "Pets allowed"), ("shared_bath", "Shared bathroom"), ("studio", "Studio"), ("dist_centre_km", "Distance to nearest centre (km)"), ("density_1km", "Listings within 1 km"), ("is_montreal", "Montreal"), ("p_stay", "Posterior prob.\\ of stay regime")] body = ["\\begin{tabular}{lrrrrrrr}", "\\toprule", "Variable & Mean & Std.\\ dev. & P10 & Median & P90 & $N$ & Mean, full \\\\", "\\midrule"] for c, lab in DESC: s = A[c].dropna().astype(float); f = df[c].dropna().astype(float) d = 0 if c in ("price", "bedrooms", "guests_count", "beds", "amenities_count", "num_images", "num_reviews", "density_1km") else 2 body.append(f"{lab} & {s.mean():,.{d}f} & {s.std():,.{d}f} & {s.quantile(.1):,.{d}f} & {s.median():,.{d}f} & {s.quantile(.9):,.{d}f} & {len(s):,} & {f.mean():,.{d}f} \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_summary", "\n".join(body)) body = ["\\begin{tabular}{lrrrrrr}", "\\toprule", "Title category & $N$ & Share & Mean (\\$) & Median (\\$) & P10 & P90 \\\\", "\\midrule"] for c in CATS: s = A.loc[A.category == c, "price"] body.append(f"{tex_escape(c)} & {len(s):,} & {len(s) / len(A):.3f} & {s.mean():,.0f} & {s.median():,.0f} & {s.quantile(.1):,.0f} & {s.quantile(.9):,.0f} \\\\") s = A.price body += ["\\midrule", f"All & {len(s):,} & 1.000 & {s.mean():,.0f} & {s.median():,.0f} & {s.quantile(.1):,.0f} & {s.quantile(.9):,.0f} \\\\", "\\bottomrule", "\\end{tabular}"] write_tex("rev_price_by_category", "\n".join(body)) # -------------------------------------------------------------------------- # 4. Nested OLS hedonic models # -------------------------------------------------------------------------- STRUCT = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "studio"] CATD = [c for c in CATS if c != "Rental unit"] AMEN = ["ln_amenities", "pets_allowed", "ln_images"] HOST = ["is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c"] LOC = ["ln_dist", "ln_density_1km"] # city fixed effects: cities with >= 15 listings in the analysis sample; the rest pooled by K-means cell MIN_CITY = 15 from sklearn.cluster import KMeans # noqa: E402 vc = A.city.value_counts() big = set(vc[vc >= MIN_CITY].index) small = A[~A.city.isin(big)] km_small = KMeans(n_clusters=25, random_state=SEED, n_init=10).fit(small[["x_km", "y_km"]].values) A["fe_city"] = np.where(A.city.isin(big), A.city, "cell_" + pd.Series(km_small.predict(A[["x_km", "y_km"]].values), index=A.index).astype(str)) KEY["n_city_fe"] = int(len(big)); KEY["n_cell_fe"] = 25; KEY["share_in_city_fe"] = float(A.city.isin(big).mean()) KEY["n_fe_groups"] = int(A.fe_city.nunique()) # finer grid-cell FE for robustness A["fe_cell60"] = KMeans(n_clusters=60, random_state=SEED, n_init=10).fit_predict(A[["x_km", "y_km"]].values) def design(d, cols, fe=None, drop_first_fe="Montreal"): X = pd.DataFrame(index=d.index) for c in cols: X[c] = (d.category == c).astype(float) if c in CATD else d[c].astype(float) if fe is not None: D = pd.get_dummies(d[fe].astype(str), prefix="fe", dtype=float) ref = f"fe_{drop_first_fe}" if f"fe_{drop_first_fe}" in D.columns else D.columns[0] D = D.drop(columns=[ref]) X = pd.concat([X, D], axis=1) return sm.add_constant(X) SPECS = { "M1": STRUCT + CATD, "M2": STRUCT + CATD + AMEN, "M3": STRUCT + CATD + AMEN + HOST, "M4": STRUCT + CATD + AMEN + HOST + LOC, "M5": STRUCT + CATD + AMEN + HOST + LOC, # + city FE } fits = {} for k, cols in SPECS.items(): X = design(A, cols, fe="fe_city" if k == "M5" else None) fits[k] = sm.OLS(A["lnp"], X).fit(cov_type="HC3") fit5_cl = sm.OLS(A["lnp"], design(A, SPECS["M5"], fe="fe_city")).fit(cov_type="cluster", cov_kwds={"groups": A["fe_city"].astype("category").cat.codes}) fit5_cell = sm.OLS(A["lnp"], design(A, SPECS["M5"], fe="fe_city")).fit(cov_type="cluster", cov_kwds={"groups": A["fe_cell60"]}) # VIF (M4, no FE) from statsmodels.stats.outliers_influence import variance_inflation_factor # noqa: E402 X4 = design(A, SPECS["M4"]) vif = pd.Series({c: variance_inflation_factor(X4.values, i) for i, c in enumerate(X4.columns) if c != "const"}) KEY["vif_max"] = float(vif.max()); KEY["vif_argmax"] = str(vif.idxmax()); KEY["vif"] = {k: float(v) for k, v in vif.items()} ROWS = STRUCT + CATD + AMEN + HOST + LOC body = ["\\begin{tabular}{l" + "r" * 6 + "}", "\\toprule", " & (1) & (2) & (3) & (4) & (5) & (5$'$) \\\\", "Dependent variable: ln(price) & M1 & M2 & M3 & M4 & M5 & M5 \\\\", "\\midrule"] for v in ROWS: line = [LAB.get(v, tex_escape(v))]; se_line = [""] for k in ("M1", "M2", "M3", "M4", "M5"): f = fits[k] if v in f.params.index: line.append(f"{f.params[v]:.3f}{stars(f.pvalues[v])}"); se_line.append(f"({f.bse[v]:.3f})") else: line.append(""); se_line.append("") if v in fit5_cl.params.index: line.append(f"{fit5_cl.params[v]:.3f}{stars(fit5_cl.pvalues[v])}"); se_line.append(f"[{fit5_cl.bse[v]:.3f}]") else: line.append(""); se_line.append("") body.append(" & ".join(line) + " \\\\"); body.append(" & ".join(se_line) + " \\\\") body += ["\\midrule", "City / cell fixed effects & No & No & No & No & Yes & Yes \\\\", "Standard errors & HC3 & HC3 & HC3 & HC3 & HC3 & Cluster \\\\"] for lab, fn in (("Observations", lambda f: f"{int(f.nobs):,}"), ("$R^2$", lambda f: f"{f.rsquared:.3f}"), ("Adjusted $R^2$", lambda f: f"{f.rsquared_adj:.3f}"), ("RMSE", lambda f: f"{np.sqrt(f.mse_resid * f.df_resid / f.nobs):.3f}"), ("AIC", lambda f: f"{f.aic:,.0f}")): body.append(lab + " & " + " & ".join(fn(fits[k]) for k in ("M1", "M2", "M3", "M4", "M5")) + f" & {fn(fit5_cl)} \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_ols_nested", "\n".join(body)) KEY["ols"] = {k: {"n": int(f.nobs), "r2": float(f.rsquared), "adj": float(f.rsquared_adj), "rmse": float(np.sqrt(f.mse_resid * f.df_resid / f.nobs)), "aic": float(f.aic), "params": {v: [float(f.params[v]), float(f.bse[v]), float(f.pvalues[v])] for v in ROWS if v in f.params.index}} for k, f in fits.items()} KEY["ols"]["M5_cluster_city"] = {v: [float(fit5_cl.params[v]), float(fit5_cl.bse[v]), float(fit5_cl.pvalues[v])] for v in ROWS} KEY["ols"]["M5_cluster_cell60"] = {v: [float(fit5_cell.params[v]), float(fit5_cell.bse[v]), float(fit5_cell.pvalues[v])] for v in ROWS} # Kennedy (1981) corrected premia for dummies in M5 f = fits["M5"] KEY["premia_M5"] = {v: float(100 * (np.exp(f.params[v] - 0.5 * f.bse[v] ** 2) - 1)) for v in CATD + ["shared_bath", "studio", "pets_allowed", "is_superhost", "is_guest_favorite", "rated"]} # F-tests of blocks (M5) def ftest(fit, vars_): R = np.zeros((len(vars_), len(fit.params))); for i, v in enumerate(vars_): R[i, list(fit.params.index).index(v)] = 1 t = fit.f_test(R); return float(np.squeeze(t.fvalue)), float(t.pvalue) KEY["ftests_M5"] = {"category": ftest(f, CATD), "host": ftest(f, HOST), "amenities": ftest(f, AMEN), "location": ftest(f, LOC)} fe_cols = [c for c in f.params.index if c.startswith("fe_")] KEY["ftests_M5"]["city_fe"] = ftest(f, fe_cols) KEY["r2_within_fe_only"] = float(sm.OLS(A["lnp"], design(A, [], fe="fe_city")).fit().rsquared) log(f"OLS done: R2 M1..M5 = {[round(fits[k].rsquared, 3) for k in fits]}") # -------------------------------------------------------------------------- # 5. Quantile regression (M5 specification) # -------------------------------------------------------------------------- TAUS = [0.10, 0.25, 0.50, 0.75, 0.90] X5 = design(A, SPECS["M5"], fe="fe_city") qfits = {t: sm.QuantReg(A["lnp"], X5).fit(q=t, max_iter=5000) for t in TAUS} QROWS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "Chalet/cottage/cabin", "House", "Condo/loft/apartment", "Hotel/B&B/hostel", "ln_amenities", "ln_images", "is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c", "ln_dist", "ln_density_1km"] body = ["\\begin{tabular}{l" + "r" * 6 + "}", "\\toprule", "Variable & OLS & $\\tau=0.10$ & $\\tau=0.25$ & $\\tau=0.50$ & $\\tau=0.75$ & $\\tau=0.90$ \\\\", "\\midrule"] for v in QROWS: line = [LAB.get(v, tex_escape(v)), f"{fits['M5'].params[v]:.3f}{stars(fits['M5'].pvalues[v])}"]; se = ["", f"({fits['M5'].bse[v]:.3f})"] for t in TAUS: q = qfits[t]; line.append(f"{q.params[v]:.3f}{stars(q.pvalues[v])}"); se.append(f"({q.bse[v]:.3f})") body.append(" & ".join(line) + " \\\\"); body.append(" & ".join(se) + " \\\\") def pseudo_r2(q, t): u = A["lnp"] - q.fittedvalues; rho = lambda u: np.sum(u * (t - (u < 0))) u0 = A["lnp"] - A["lnp"].quantile(t); return 1 - rho(u) / rho(u0) body += ["\\midrule", "Pseudo-$R^2$ / $R^2$ & " + f"{fits['M5'].rsquared:.3f} & " + " & ".join(f"{pseudo_r2(qfits[t], t):.3f}" for t in TAUS) + " \\\\", f"Observations & \\multicolumn{{6}}{{c}}{{{len(A):,}}} \\\\", "City / cell fixed effects & \\multicolumn{6}{c}{Yes} \\\\", "\\bottomrule", "\\end{tabular}"] write_tex("rev_quantile", "\n".join(body)) KEY["quantile"] = {str(t): {v: [float(qfits[t].params[v]), float(qfits[t].bse[v]), float(qfits[t].pvalues[v])] for v in QROWS} for t in TAUS} KEY["quantile_pseudoR2"] = {str(t): float(pseudo_r2(qfits[t], t)) for t in TAUS} # interquantile difference tests (0.90 - 0.10) by pairwise bootstrap (B = 200) rng = np.random.default_rng(SEED) B = 100; diffs = {v: [] for v in ("ln_bedrooms", "ln_guests", "is_superhost", "rating_c", "ln_reviews", "ln_dist")} cache = TAB / "rev_iq_bootstrap.json" if cache.exists(): diffs = json.loads(cache.read_text()); log("inter-quantile bootstrap loaded from cache") else: for b in range(B): idx = rng.integers(0, len(A), len(A)) q1 = sm.QuantReg(A["lnp"].values[idx], X5.values[idx]).fit(q=0.10, max_iter=3000) q9 = sm.QuantReg(A["lnp"].values[idx], X5.values[idx]).fit(q=0.90, max_iter=3000) for v in diffs: j = list(X5.columns).index(v); diffs[v].append(float(q9.params[j] - q1.params[j])) cache.write_text(json.dumps(diffs)) KEY["iq_bootstrap_B"] = B KEY["iq_diff_90_10"] = {} for v, arr in diffs.items(): arr = np.array(arr); est = qfits[0.90].params[v] - qfits[0.10].params[v] KEY["iq_diff_90_10"][v] = [float(est), float(arr.std(ddof=1)), float(2 * (1 - stats.norm.cdf(abs(est / arr.std(ddof=1)))))] log("quantile done") fig, axes = plt.subplots(2, 3, figsize=(9.5, 5.2)); axes = axes.ravel() for ax, v in zip(axes, ["ln_bedrooms", "ln_guests", "ln_bathrooms", "is_superhost", "rating_c", "ln_dist"]): est = [qfits[t].params[v] for t in TAUS]; se = [qfits[t].bse[v] for t in TAUS] ax.fill_between(TAUS, np.array(est) - 1.96 * np.array(se), np.array(est) + 1.96 * np.array(se), color="#c6dbef") ax.plot(TAUS, est, "o-", color="#08519c", ms=3) ax.axhline(fits["M5"].params[v], color="#de2d26", ls="--", lw=1); ax.axhline(0, color="grey", lw=0.5) ax.set_title(LAB[v], fontsize=9); ax.set_xticks(TAUS) fig.supxlabel("Quantile $\\tau$ of ln(price)"); fig.supylabel("Coefficient") fig.tight_layout(); savefig(fig, "rev_fig_quantile") # -------------------------------------------------------------------------- # 6. Spatial econometrics (M4 specification; W = k-nearest neighbours on projected km) # -------------------------------------------------------------------------- import libpysal, spreg # noqa: E402 from esda.moran import Moran # noqa: E402 coords = A[["x_km", "y_km"]].values def knn_w(k): w = libpysal.weights.KNN.from_array(coords, k=k); w.transform = "r"; return w W5, W10 = knn_w(5), knn_w(10) KEY["moran"] = {} for lab, f, w in (("M4_k5", fits["M4"], W5), ("M5_k5", fits["M5"], W5), ("lnp_k5", None, W5), ("M4_k10", fits["M4"], W10)): z = A["lnp"].values if f is None else f.resid.values mi = Moran(z, w, permutations=999) KEY["moran"][lab] = [float(mi.I), float(mi.z_norm), float(mi.p_norm)] Xs = design(A, SPECS["M4"]).drop(columns="const") y = A[["lnp"]].values; Xv = Xs.values; names = list(Xs.columns) sar = spreg.ML_Lag(y, Xv, W5, name_x=names, name_y="lnp", method="full") sem = spreg.ML_Error(y, Xv, W5, name_x=names, name_y="lnp", method="full") sdm = spreg.ML_Lag(y, Xv, W5, name_x=names, name_y="lnp", method="full", slx_lags=1) sar10 = spreg.ML_Lag(y, Xv, W10, name_x=names, name_y="lnp", method="full") ols4 = fits["M4"] def sp_row(mod, v, lagged=False): nm_ = [n for n in mod.name_x] try: j = nm_.index(("W_" + v) if lagged else v) except ValueError: return "", "" b = float(mod.betas[j, 0]); se_ = float(mod.std_err[j]); p = float(mod.z_stat[j][1]) return f"{b:.3f}{stars(p)}", f"({se_:.3f})" SROWS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "Chalet/cottage/cabin", "House", "Condo/loft/apartment", "Hotel/B&B/hostel", "ln_amenities", "ln_images", "is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c", "ln_dist", "ln_density_1km"] body = ["\\begin{tabular}{lrrrrr}", "\\toprule", "Variable & OLS (M4) & SAR & SEM & SDM & SDM, $\\mathbf{W}x$ \\\\", "\\midrule"] for v in SROWS: l = [LAB.get(v, tex_escape(v)), f"{ols4.params[v]:.3f}{stars(ols4.pvalues[v])}"]; s = ["", f"({ols4.bse[v]:.3f})"] for m in (sar, sem, sdm): a, b_ = sp_row(m, v); l.append(a); s.append(b_) a, b_ = sp_row(sdm, v, lagged=True); l.append(a); s.append(b_) body.append(" & ".join(l) + " \\\\"); body.append(" & ".join(s) + " \\\\") rho_p = lambda m: float(m.z_stat[-1][1]) body += ["\\midrule", f"$\\rho$ (spatial lag) & & {float(np.squeeze(sar.rho)):.3f}{stars(rho_p(sar))} & & {float(np.squeeze(sdm.rho)):.3f}{stars(rho_p(sdm))} & \\\\", f" & & ({float(sar.std_err[-1]):.3f}) & & ({float(sdm.std_err[-1]):.3f}) & \\\\", f"$\\lambda$ (spatial error) & & & {float(np.squeeze(sem.lam)):.3f}{stars(rho_p(sem))} & & \\\\", f" & & & ({float(sem.std_err[-1]):.3f}) & & \\\\", f"Log-likelihood & {ols4.llf:,.1f} & {float(np.squeeze(sar.logll)):,.1f} & {float(np.squeeze(sem.logll)):,.1f} & {float(np.squeeze(sdm.logll)):,.1f} & \\\\", f"AIC & {ols4.aic:,.0f} & {float(np.squeeze(sar.aic)):,.0f} & {float(np.squeeze(sem.aic)):,.0f} & {float(np.squeeze(sdm.aic)):,.0f} & \\\\", f"Pseudo-$R^2$ & {ols4.rsquared:.3f} & {float(np.squeeze(sar.pr2)):.3f} & {float(np.squeeze(sem.pr2)):.3f} & {float(np.squeeze(sdm.pr2)):.3f} & \\\\", f"Observations & \\multicolumn{{5}}{{c}}{{{len(A):,}}} \\\\", "\\bottomrule", "\\end{tabular}"] write_tex("rev_spatial", "\n".join(body)) KEY["spatial"] = {"rho_sar": [float(np.squeeze(sar.rho)), float(sar.std_err[-1]), rho_p(sar)], "lam_sem": [float(np.squeeze(sem.lam)), float(sem.std_err[-1]), rho_p(sem)], "rho_sdm": [float(np.squeeze(sdm.rho)), float(sdm.std_err[-1]), rho_p(sdm)], "rho_sar_k10": [float(np.squeeze(sar10.rho)), float(sar10.std_err[-1]), rho_p(sar10)], "ll": {"ols": float(ols4.llf), "sar": float(np.squeeze(sar.logll)), "sem": float(np.squeeze(sem.logll)), "sdm": float(np.squeeze(sdm.logll))}, "aic": {"ols": float(ols4.aic), "sar": float(np.squeeze(sar.aic)), "sem": float(np.squeeze(sem.aic)), "sdm": float(np.squeeze(sdm.aic))}} # LR test SDM vs SAR (restriction theta = 0) lr = 2 * (float(np.squeeze(sdm.logll)) - float(np.squeeze(sar.logll))); KEY["spatial"]["lr_sdm_vs_sar"] = [lr, int(len(names)), float(1 - stats.chi2.cdf(lr, len(names)))] # LeSage-Pace impacts (SAR and SDM), dense inverse Wd = W5.full()[0] I_n = np.eye(len(A)); S_sar = np.linalg.inv(I_n - float(np.squeeze(sar.rho)) * Wd); S_sdm = np.linalg.inv(I_n - float(np.squeeze(sdm.rho)) * Wd) imp = {} for v in SROWS: j = names.index(v); b_sar = float(sar.betas[j + 1, 0]) d = np.trace(S_sar) / len(A) * b_sar; t = S_sar.sum() / len(A) * b_sar b_sdm = float(sdm.betas[j + 1, 0]); th = float(sdm.betas[sdm.name_x.index("W_" + v), 0]) M = S_sdm @ (I_n * b_sdm + Wd * th) d2 = np.trace(M) / len(A); t2 = M.sum() / len(A) imp[v] = {"sar": [d, t - d, t], "sdm": [d2, t2 - d2, t2]} # simulated SEs for SAR and SDM impacts (300 draws from the ML asymptotic distribution; exact traces via the # complex eigenvalues of the asymmetric KNN matrix: tr[(I - rW)^-1] = sum_i 1/(1 - r*lambda_i)) rng = np.random.default_rng(SEED) evals = np.linalg.eigvals(Wd) def _sim(mod, has_theta): V = np.asarray(mod.vm); th_hat = mod.betas.ravel(); dr_all = rng.multivariate_normal(th_hat, V, size=300) out = {v: [] for v in SROWS} for dr in dr_all: r_ = dr[-1] if abs(r_) >= 1: continue tr = np.mean(1 / (1 - r_ * evals)).real; trw = np.mean(evals / (1 - r_ * evals)).real; tot = 1 / (1 - r_) for v in SROWS: b = dr[names.index(v) + 1]; th = dr[mod.name_x.index("W_" + v)] if has_theta else 0.0 direct = tr * b + trw * th; total = (b + th) * tot out[v].append([direct, total - direct, total]) return {v: np.array(a).std(axis=0, ddof=1) for v, a in out.items()} sd_sar = _sim(sar, False); sd_sdm = _sim(sdm, True) pz = lambda est, sd_: stars(2 * (1 - stats.norm.cdf(abs(est / sd_)))) if sd_ > 0 else "" body = ["\\begin{tabular}{lrrrrrr}", "\\toprule", " & \\multicolumn{3}{c}{SAR} & \\multicolumn{3}{c}{SDM} \\\\", "\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}", "Variable & Direct & Indirect & Total & Direct & Indirect & Total \\\\", "\\midrule"] KEY["impacts"] = {} for v in SROWS: a = imp[v]; s1_ = sd_sar[v]; s2_ = sd_sdm[v] body.append(f"{LAB.get(v, tex_escape(v))} & " + " & ".join(f"{a['sar'][i]:.3f}{pz(a['sar'][i], s1_[i])}" for i in range(3)) + " & " + " & ".join(f"{a['sdm'][i]:.3f}{pz(a['sdm'][i], s2_[i])}" for i in range(3)) + " \\\\") body.append(" & " + " & ".join(f"({s1_[i]:.3f})" for i in range(3)) + " & " + " & ".join(f"({s2_[i]:.3f})" for i in range(3)) + " \\\\") KEY["impacts"][v] = {"sar": a["sar"], "sar_se": s1_.tolist(), "sdm": a["sdm"], "sdm_se": s2_.tolist()} body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_impacts", "\n".join(body)) KEY["spatial"]["sar_multiplier"] = float(S_sar.sum() / len(A)); KEY["spatial"]["sar_direct_ratio"] = float(np.trace(S_sar) / len(A)) log("spatial done") # -------------------------------------------------------------------------- # 7. Machine learning: random vs spatially blocked cross-validation, SHAP # -------------------------------------------------------------------------- from sklearn.model_selection import KFold, GroupKFold # noqa: E402 from sklearn.linear_model import ElasticNetCV, LinearRegression # noqa: E402 from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor # noqa: E402 from sklearn.preprocessing import StandardScaler # noqa: E402 from sklearn.pipeline import make_pipeline # noqa: E402 import xgboost as xgb, lightgbm as lgb # noqa: E402 XML = design(A, SPECS["M5"], fe="fe_city").drop(columns="const") XML.columns = [re.sub(r"[^0-9a-zA-Z_]", "_", c) for c in XML.columns] yv = A["lnp"].values; pv = A["price"].values blocks = KMeans(n_clusters=20, random_state=SEED, n_init=10).fit_predict(coords) KEY["ml_n_blocks"] = 20 MODELS = { "OLS (M5)": lambda: LinearRegression(), "Elastic net": lambda: make_pipeline(StandardScaler(), ElasticNetCV(l1_ratio=[.1, .5, .9], cv=5, random_state=SEED, n_jobs=-1)), "Random forest": lambda: RandomForestRegressor(n_estimators=500, min_samples_leaf=3, max_features=0.5, random_state=SEED, n_jobs=-1), "Gradient boosting": lambda: GradientBoostingRegressor(n_estimators=500, learning_rate=0.05, max_depth=3, subsample=0.8, random_state=SEED), "XGBoost": lambda: xgb.XGBRegressor(n_estimators=600, learning_rate=0.03, max_depth=5, subsample=0.8, colsample_bytree=0.8, min_child_weight=3, random_state=SEED, n_jobs=-1, verbosity=0), "LightGBM": lambda: lgb.LGBMRegressor(n_estimators=600, learning_rate=0.03, num_leaves=15, min_child_samples=10, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, random_state=SEED, n_jobs=-1, verbose=-1), } def cv_eval(splitter, groups=None): out = {} for name, mk in MODELS.items(): yhat = np.zeros_like(yv) for tr, te in splitter.split(XML, yv, groups): m = mk(); m.fit(XML.iloc[tr], yv[tr]); yhat[te] = m.predict(XML.iloc[te]) r2 = 1 - np.sum((yv - yhat) ** 2) / np.sum((yv - yv.mean()) ** 2) rmse = np.sqrt(np.mean((yv - yhat) ** 2)); mae = np.mean(np.abs(yv - yhat)) mdape = np.median(np.abs(np.exp(yhat) - pv) / pv) * 100 out[name] = [r2, rmse, mae, mdape] return out res_rand = cv_eval(KFold(5, shuffle=True, random_state=SEED)) res_block = cv_eval(GroupKFold(5), groups=blocks) body = ["\\begin{tabular}{lrrrrrrrr}", "\\toprule", " & \\multicolumn{4}{c}{Random 5-fold CV} & \\multicolumn{4}{c}{Spatially blocked 5-fold CV} \\\\", "\\cmidrule(lr){2-5}\\cmidrule(lr){6-9}", "Model & $R^2$ & RMSE & MAE & MdAPE (\\%) & $R^2$ & RMSE & MAE & MdAPE (\\%) \\\\", "\\midrule"] for name in MODELS: a, b_ = res_rand[name], res_block[name] body.append(f"{name} & {a[0]:.3f} & {a[1]:.3f} & {a[2]:.3f} & {a[3]:.1f} & {b_[0]:.3f} & {b_[1]:.3f} & {b_[2]:.3f} & {b_[3]:.1f} \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_ml_cv", "\n".join(body)) KEY["ml"] = {"random": res_rand, "blocked": res_block} log("ML CV done") # SHAP on XGBoost fitted to the full analysis sample import shap # noqa: E402 xg = MODELS["XGBoost"](); xg.fit(XML, yv) # exact TreeSHAP values computed natively by XGBoost (pred_contribs); last column is the bias term contrib = xg.get_booster().predict(xgb.DMatrix(XML), pred_contribs=True) sv = contrib[:, :-1]; KEY['shap_bias'] = float(contrib[0, -1]) assert np.allclose(contrib.sum(axis=1), xg.predict(XML), atol=1e-3) imp_shap = pd.Series(np.abs(sv).mean(axis=0), index=XML.columns) fe_mask = imp_shap.index.str.startswith("fe_") imp_agg = imp_shap[~fe_mask].copy(); imp_agg["City / cell fixed effects (all)"] = imp_shap[fe_mask].sum() imp_agg = imp_agg.sort_values(ascending=False) LAB2 = {re.sub(r"[^0-9a-zA-Z_]", "_", k): v for k, v in LAB.items()} body = ["\\begin{tabular}{lrr}", "\\toprule", "Feature & Mean $|\\phi_k|$ & Share of total (\\%) \\\\", "\\midrule"] tot = imp_shap.sum() for k, v in imp_agg.head(16).items(): body.append(f"{LAB2.get(k, tex_escape(k))} & {v:.4f} & {100 * v / tot:.1f} \\\\") body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_shap_importance", "\n".join(body)) KEY["shap_top"] = {k: float(v) for k, v in imp_agg.head(10).items()}; KEY["shap_total"] = float(tot) KEY["shap_fe_share"] = float(imp_shap[fe_mask].sum() / tot) # beeswarm of the 12 most important non-FE features top = imp_shap[~fe_mask].sort_values(ascending=False).head(12).index Xd = XML[top].copy(); Xd.columns = [LAB2.get(c, c) for c in top] fig = plt.figure(figsize=(7, 4.6)) shap.summary_plot(sv[:, [list(XML.columns).index(c) for c in top]], Xd, show=False, plot_size=None, max_display=12) plt.gcf().set_size_inches(7, 4.6); plt.xlabel("SHAP value (contribution to predicted ln price)") savefig(plt.gcf(), "rev_fig_shap") log("SHAP done") # -------------------------------------------------------------------------- # 8. Robustness of the M5 specification # -------------------------------------------------------------------------- RKEYS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "is_superhost", "is_guest_favorite", "ln_reviews", "rating_c", "ln_dist"] def refit(d, fe="fe_city", label=""): d = d.reset_index(drop=True) if fe == "fe_cell60": X = design(d, SPECS["M5"], fe="fe_cell60", drop_first_fe="0") else: X = design(d, SPECS["M5"], fe=fe) X = X.loc[:, (X != 0).any(axis=0)] if label.startswith("Weighted"): return sm.WLS(d["lnp"], X, weights=1 - d["p_stay"]).fit(cov_type="HC3") return sm.OLS(d["lnp"], X).fit(cov_type="HC3") def _fe(d): return d.assign(fe_city=lambda d: np.where(d.city.isin(big), d.city, "cell_" + pd.Series(km_small.predict(d[["x_km", "y_km"]].values), index=d.index).astype(str))) variants = [ ("Baseline: nightly regime ($\\hat{p}_i < 0.5$)", A), ("Nightly regime, strict ($\\hat{p}_i < 0.25$)", _fe(df[df.p_stay < 0.25])), ("Nightly regime, lenient ($\\hat{p}_i < 0.75$)", _fe(df[df.p_stay < 0.75])), ("Price $\\leq$ \\$1{,}000 rule", _fe(df[df.price <= 1000])), ("Price $\\leq$ \\$600 rule", _fe(df[df.price <= 600])), ("All priced listings (bimodal)", _fe(df)), ("Excluding hotels, B\\&Bs, hostels", A[A.category != "Hotel/B&B/hostel"]), ("Rated listings only", A[A.rated == 1]), ("Montreal only", A[A.is_montreal == 1]), ("Outside Montreal", A[A.is_montreal == 0]), ("Trimmed 1\\%/99\\% within sample", A[(A.price > A.price.quantile(.01)) & (A.price < A.price.quantile(.99))]), ("Weighted by $1 - \\hat{p}_i$", A), ] body = ["\\begin{tabular}{l" + "r" * (len(RKEYS) + 2) + "}", "\\toprule", "Sample / specification & " + " & ".join(["Bedr.", "Bathr.", "Guests", "Superhost", "Guest Fav.", "ln reviews", "Rating", "ln dist."]) + " & $N$ & $R^2$ \\\\", "\\midrule"] KEY["robust"] = {} for lab, d in variants: f = refit(d, label=lab) cells = [f"{f.params[v]:.3f}{stars(f.pvalues[v])}" if v in f.params.index else "" for v in RKEYS] ses = [f"({f.bse[v]:.3f})" if v in f.params.index else "" for v in RKEYS] body.append(f"{lab} & " + " & ".join(cells) + f" & {int(f.nobs):,} & {f.rsquared:.3f} \\\\") body.append(" & " + " & ".join(ses) + " & & \\\\") KEY["robust"][lab] = {"n": int(f.nobs), "r2": float(f.rsquared), **{v: [float(f.params[v]), float(f.bse[v]), float(f.pvalues[v])] for v in RKEYS if v in f.params.index}} # grid-cell FE (60 cells) f = refit(A, fe="fe_cell60") body.append("60 spatial-cell fixed effects & " + " & ".join(f"{f.params[v]:.3f}{stars(f.pvalues[v])}" for v in RKEYS) + f" & {int(f.nobs):,} & {f.rsquared:.3f} \\\\") body.append(" & " + " & ".join(f"({f.bse[v]:.3f})" for v in RKEYS) + " & & \\\\") KEY["robust"]["cell60"] = {"n": int(f.nobs), "r2": float(f.rsquared), **{v: [float(f.params[v]), float(f.bse[v]), float(f.pvalues[v])] for v in RKEYS}} # Montreal x bedrooms heterogeneity body += ["\\bottomrule", "\\end{tabular}"] write_tex("rev_robustness", "\n".join(body)) log("robustness done") # -------------------------------------------------------------------------- # 9. Figures: map, price by category, binned partial relationships # -------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(7, 5.2)) sc = ax.scatter(A["long"], A["lat"], c=A["lnp"], s=5, cmap="viridis", alpha=0.7, linewidths=0) cb = fig.colorbar(sc, ax=ax, shrink=0.8); cb.set_label("Displayed price (CAD)") tk = [50, 100, 200, 400, 800]; cb.set_ticks(np.log(tk)); cb.set_ticklabels([str(t) for t in tk]) ax.set_xlabel("Longitude"); ax.set_ylabel("Latitude"); ax.set_xlim(-80, -61); ax.set_ylim(44.8, 50.5) for nmc, (la, lo) in centres.items(): if nmc in ("Montréal", "Québec", "Gatineau", "Mont-Tremblant", "Sherbrooke", "Saguenay", "Rimouski"): ax.annotate(nmc, (lo, la), fontsize=7, xytext=(3, 3), textcoords="offset points") savefig(fig, "rev_fig_map") fig, ax = plt.subplots(figsize=(7, 3.8)) data = [A.loc[A.category == c, "price"] for c in CATS] ax.boxplot(data, vert=False, showfliers=False, widths=0.6, medianprops={"color": "#de2d26"}) ax.set_yticks(range(1, len(CATS) + 1)); ax.set_yticklabels([f"{c} (n={len(d)})" for c, d in zip(CATS, data)], fontsize=8) ax.set_xscale("log"); ax.set_xlabel("Displayed price (CAD, log scale)"); savefig(fig, "rev_fig_price_by_category") fig, axes = plt.subplots(1, 2, figsize=(9, 3.4)) r4 = fits["M4"] # partial residual plots: residual + beta*x against x, binned for ax, v, xl in zip(axes, ("ln_dist", "rating_c"), ("Distance to nearest centre (km)", "Rating (rated listings)")): d = A.copy(); d["pr"] = r4.resid + r4.params[v] * d[v] if v == "rating_c": d = d[d.rated == 1]; xvar = d["rating_raw"] else: xvar = d["dist_centre_km"] bins = pd.qcut(xvar, 20, duplicates="drop"); g = d.groupby(bins, observed=True) xm = g.apply(lambda s: xvar.loc[s.index].mean()); ym = g["pr"].mean(); se = g["pr"].sem() ax.errorbar(xm, ym, yerr=1.96 * se, fmt="o", ms=3, color="#08519c", ecolor="#9ecae1") ax.set_xlabel(xl); ax.set_ylabel("Partial residual of ln(price)") if v == "ln_dist": ax.set_xscale("symlog", linthresh=1) fig.tight_layout(); savefig(fig, "rev_fig_partials") # -------------------------------------------------------------------------- # 10. Save key numbers # -------------------------------------------------------------------------- KEY["runtime_s"] = time.perf_counter() - T0 (TAB / "rev_key_numbers.json").write_text(json.dumps(KEY, indent=1, default=float), encoding="utf-8") A.drop(columns=["name"]).to_csv(ROOT / "data" / "processed" / "analysis_sample.csv", index=False) log("done")