Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1#!/usr/bin/env python32"""310_revision_reanalysis.py -- Editorial-board revision of UQO WP 4 (October 2026).45Re-estimates every result reported in the revised paper from the only locally6available copy of the Quebec Airbnb scrape, ``data/processed/airbnb_clean_wp5.parquet``7(3,456 priced listings out of 5,000 scraped; prices winsorised at p1/p99 and8ratings median-imputed by the companion paper's cleaning step -- the imputation is9reversed here by parsing the star rating from the structured title).1011Outputs (all regenerated, never typed):12 results/tables/rev_*.tex / .csv / .json13 paper/figures/rev_*.pdf (and results/figures/)1415Run time: about 3-4 minutes on an Apple-silicon laptop.16"""17from __future__ import annotations1819import json20import re21import sys22import time23import warnings24from pathlib import Path2526import numpy as np27import pandas as pd28import statsmodels.api as sm29from scipy import stats3031warnings.filterwarnings("ignore")32ROOT = Path(__file__).resolve().parent.parent33DATA = ROOT / "data" / "processed" / "airbnb_clean_wp5.parquet"34TAB = ROOT / "results" / "tables"35FIG_R = ROOT / "results" / "figures"36FIG_P = ROOT / "paper" / "figures"37for d in (TAB, FIG_R, FIG_P):38 d.mkdir(parents=True, exist_ok=True)3940import matplotlib41matplotlib.use("Agg")42import matplotlib.pyplot as plt43plt.rcParams.update({"font.family": "serif", "font.size": 10, "axes.spines.top": False,44 "axes.spines.right": False, "figure.dpi": 150})4546SEED = 4247NIGHTLY_CAP = 1000.0 # main analysis sample: price <= $1,00048KEY = {} # numbers quoted in the text49T0 = time.perf_counter()505152def log(msg):53 print(f"[{time.perf_counter() - T0:6.1f}s] {msg}", flush=True)545556def savefig(fig, stem):57 for d in (FIG_R, FIG_P):58 fig.savefig(d / f"{stem}.pdf", bbox_inches="tight")59 plt.close(fig)606162def stars(p):63 return "***" if p < 0.01 else "**" if p < 0.05 else "*" if p < 0.10 else ""646566def fmt(x, d=3):67 return f"{x:.{d}f}"686970def tex_escape(s):71 return s.replace("&", r"\&").replace("%", r"\%").replace("_", r"\_").replace("$", r"\$")727374def write_tex(name, body):75 (TAB / f"{name}.tex").write_text(body, encoding="utf-8")76 log(f"wrote {name}.tex")777879# --------------------------------------------------------------------------80# 1. Load and parse the structured title81# --------------------------------------------------------------------------82df = pd.read_parquet(DATA)83df["name"] = df["name"].fillna("").astype(str)84N_PRICED = len(df)85nm = df["name"]8687# star rating from title (reverses the median imputation of the parquet)88star = nm.str.extract(r"★\s*(\d+(?:[.,]\d+)?)", expand=False).str.replace(",", ".", regex=False).astype(float)89df["rating_raw"] = star90df["rated"] = star.notna().astype(int)91df["is_new"] = nm.str.contains(r"★\s*New\b|·\s*New\s*·", regex=True).astype(int)92# 37 German-locale titles ("★4,97") were truncated to 4.0 in the scraped rating field; the title value is retained93KEY["n_rating_truncated_fixed"] = int(((star.notna()) & ((star - df["rating"]).abs() > 0.011)).sum())94KEY["n_unrated"] = int((df["rated"] == 0).sum())95KEY["n_new"] = int(df["is_new"].sum())9697# title category (English template "Category in City · ..."; 50 German-locale titles fall back to property_type)98cat = nm.str.extract(r"^(.+?)\s+in\s+", expand=False)99CAT_MAP = {100 "Rental unit": "Rental unit",101 "Chalet": "Chalet/cottage/cabin", "Cottage": "Chalet/cottage/cabin", "Cabin": "Chalet/cottage/cabin",102 "Vacation home": "Chalet/cottage/cabin", "Holiday home": "Chalet/cottage/cabin", "Nature lodge": "Chalet/cottage/cabin",103 "Home": "House", "Townhouse": "House", "Bungalow": "House", "Villa": "House", "Farm stay": "House",104 "Earthen home": "House", "Earth homes": "House", "Castle": "House",105 "Condo": "Condo/loft/apartment", "Loft": "Condo/loft/apartment", "Serviced apartment": "Condo/loft/apartment",106 "Apartment": "Condo/loft/apartment", "Guest suite": "Condo/loft/apartment",107 "Hotel": "Hotel/B&B/hostel", "Bed and breakfast": "Hotel/B&B/hostel", "Aparthotel": "Hotel/B&B/hostel",108 "Boutique hotel": "Hotel/B&B/hostel", "Hostel": "Hotel/B&B/hostel", "Guesthouse": "Hotel/B&B/hostel",109 "Guest house": "Hotel/B&B/hostel", "Resort": "Hotel/B&B/hostel", "casa particular": "Hotel/B&B/hostel",110}111PT_FALLBACK = {"Rental unit": "Rental unit", "Cabin/Chalet": "Chalet/cottage/cabin", "House": "House",112 "Condo": "Condo/loft/apartment", "Apartment": "Condo/loft/apartment", "Other": "Other/unconventional"}113df["category"] = cat.map(CAT_MAP)114df.loc[df["category"].isna() & cat.notna(), "category"] = "Other/unconventional"115df.loc[cat.isna(), "category"] = df.loc[cat.isna(), "property_type"].map(PT_FALLBACK).fillna("Other/unconventional")116CATS = ["Rental unit", "Chalet/cottage/cabin", "House", "Condo/loft/apartment", "Hotel/B&B/hostel", "Other/unconventional"]117df["category"] = pd.Categorical(df["category"], categories=CATS)118119df["shared_bath"] = nm.str.contains(r"shared\s+(?:half-?)?bath|Gemeinschaftsbad", case=False, regex=True).astype(int)120df["studio"] = nm.str.contains(r"\bStudio\b", case=False, regex=True).astype(int)121df["beds"] = nm.str.extract(r"(\d+)\s*(?:beds?|Bett(?:en)?)\b", expand=False).astype(float)122# the scraper stored bathrooms = 0 whenever the title qualified the bathroom ("1 shared bath", "1 private bath",123# "Half-bath"); recover the count from the title; studios have bedrooms = 0 and are recoded to one sleeping room124baths_title = nm.str.extract(r"(\d+(?:[.,]\d+)?)\s*(?:shared |private )?(?:bath|Badezimmer)", expand=False).str.replace(",", ".", regex=False).astype(float)125half_bath = nm.str.contains(r"half-?bath", case=False, regex=True)126df["bathrooms_raw"] = df["bathrooms"]127fix = (df["bathrooms"] == 0) & baths_title.notna()128df.loc[fix, "bathrooms"] = baths_title[fix]129df.loc[(df["bathrooms"] == 0) & half_bath, "bathrooms"] = 0.5130KEY["n_bathrooms_recovered"] = int(fix.sum() + ((df["bathrooms_raw"] == 0) & half_bath & ~fix).sum())131bedr_title = nm.str.extract(r"(\d+)\s*(?:bedrooms?|Schlafzimmer)", expand=False).astype(float)132fixb = (df["bedrooms"] == 0) & bedr_title.notna() & (df["studio"] == 0)133df.loc[fixb, "bedrooms"] = bedr_title[fixb]134df["bedrooms_raw"] = df["bedrooms"]135df.loc[df["studio"] == 1, "bedrooms"] = 1136KEY["n_bedrooms_recovered"] = int(fixb.sum()); KEY["n_studio_recoded"] = int(df["studio"].sum())137138# numeric prep139for c in ("is_superhost", "is_guest_favorite", "pets_allowed"):140 df[c] = df[c].astype(int)141df["price"] = df["price_numeric"].astype(float)142df["lnp"] = np.log(df["price"])143df["ln_bedrooms"] = np.log(df["bedrooms"].clip(lower=1))144df["ln_bathrooms"] = np.log1p(df["bathrooms"])145df["ln_guests"] = np.log(df["guests_count"].clip(lower=1))146df["ln_amenities"] = np.log1p(df["amenities_count"])147df["ln_images"] = np.log1p(df["num_images"])148df["ln_reviews"] = np.log1p(df["num_reviews"])149df["rating_c"] = np.where(df["rated"] == 1, df["rating_raw"] - 4.8, 0.0) # centred at 4.8, zero when unrated150df["is_montreal"] = (df["city"] == "Montreal").astype(int)151152# geography: projected km, distance to nearest centre, local density153R = 6371.0154lat0, lon0 = df["lat"].mean(), df["long"].mean()155df["x_km"] = R * np.radians(df["long"] - lon0) * np.cos(np.radians(lat0))156df["y_km"] = R * np.radians(df["lat"] - lat0)157sys.path.insert(0, str(ROOT))158from src.config import CITY_CENTERS # noqa: E402159centres = {v.name: (v.lat, v.lon) for v in CITY_CENTERS.values()}160cx = np.array([R * np.radians(lo - lon0) * np.cos(np.radians(lat0)) for la, lo in centres.values()])161cy = np.array([R * np.radians(la - lat0) for la, lo in centres.values()])162dd = np.sqrt((df["x_km"].values[:, None] - cx[None, :]) ** 2 + (df["y_km"].values[:, None] - cy[None, :]) ** 2)163df["dist_centre_km"] = dd.min(axis=1)164df["ln_dist"] = np.log1p(df["dist_centre_km"])165from sklearn.neighbors import KDTree # noqa: E402166tree = KDTree(df[["x_km", "y_km"]].values)167df["density_1km"] = tree.query_radius(df[["x_km", "y_km"]].values, r=1.0, count_only=True) - 1168df["density_500m"] = tree.query_radius(df[["x_km", "y_km"]].values, r=0.5, count_only=True) - 1169df["ln_density_1km"] = np.log1p(df["density_1km"])170KEY["n_centres"] = len(centres)171172# --------------------------------------------------------------------------173# 2. The price field: bimodality, mixture, regime model174# --------------------------------------------------------------------------175from sklearn.mixture import GaussianMixture # noqa: E402176gm = GaussianMixture(2, random_state=SEED, n_init=5).fit(df[["lnp"]].values)177order = np.argsort(gm.means_.ravel())178mu = gm.means_.ravel()[order]; sd = np.sqrt(gm.covariances_.ravel()[order]); wt = gm.weights_[order]179grid = np.linspace(df["lnp"].min(), df["lnp"].max(), 4000)180post_hi = gm.predict_proba(grid[:, None])[:, order[1]]181thr_ln = grid[np.argmax(post_hi >= 0.5)]182KEY["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]),183 "w_lo": float(wt[0]), "w_hi": float(wt[1]), "threshold": float(np.exp(thr_ln))}184KEY["share_gt1000"] = float((df["price"] > 1000).mean())185KEY["n_gt1000"] = int((df["price"] > 1000).sum())186KEY["share_gt_mix"] = float((df["price"] > np.exp(thr_ln)).mean())187KEY["share_le600"] = float((df["price"] <= 600).mean())188KEY["price_median_all"] = float(df["price"].median()); KEY["price_mean_all"] = float(df["price"].mean())189KEY["bimodality_coefficient"] = float((stats.skew(df["lnp"]) ** 2 + 1) / (stats.kurtosis(df["lnp"]) + 3 * (N_PRICED - 1) ** 2 / ((N_PRICED - 2) * (N_PRICED - 3))))190KEY["hartigan_note"] = "dip test not run (package not installed)"191192# ratio of medians by bedrooms193lo, hi = df[df.price <= 1000], df[df.price > 1000]194ratio = (hi.groupby(hi.bedrooms.clip(upper=5))["price"].median() / lo.groupby(lo.bedrooms.clip(upper=5))["price"].median())195KEY["median_ratio_by_bedrooms"] = {str(k): float(v) for k, v in ratio.items()}196rows = []197for b in range(6):198 l = lo[lo.bedrooms.clip(upper=5) == b]; h = hi[hi.bedrooms.clip(upper=5) == b]199 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()))200body = ["\\begin{tabular}{lrrrrr}", "\\toprule",201 "Bedrooms & $N$ & Median (\\$) & $N$ & Median (\\$) & Ratio \\\\",202 " & \\multicolumn{2}{c}{$\\leq$ \\$1{,}000} & \\multicolumn{2}{c}{$>$ \\$1{,}000} & \\\\", "\\midrule"]203for r in rows:204 body.append(f"{r[0]} & {r[1]:,} & {r[2]:,.0f} & {r[3]:,} & {r[4]:,.0f} & {r[5]:.1f} \\\\")205body += ["\\bottomrule", "\\end{tabular}"]206write_tex("rev_price_regimes", "\n".join(body))207208# ---- two-regime mixture of hedonic regressions (EM) ---------------------------------209# 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),210# P(z_i = 1 | w_i) = Lambda(w_i'gamma). Common slopes; the stay regime shifts the intercept by delta211# (a log stay length) and has its own variance (stay lengths vary).212_CATD = [c for c in CATS if c != "Rental unit"]213def _X_mix(d):214 X = pd.DataFrame({"ln_bedrooms": d["ln_bedrooms"], "ln_bathrooms": d["ln_bathrooms"], "ln_guests": d["ln_guests"],215 "shared_bath": d["shared_bath"], "studio": d["studio"], "ln_amenities": d["ln_amenities"],216 "ln_images": d["ln_images"], "is_superhost": d["is_superhost"], "is_guest_favorite": d["is_guest_favorite"],217 "ln_reviews": d["ln_reviews"], "rated": d["rated"], "rating_c": d["rating_c"], "ln_dist": d["ln_dist"],218 "ln_density_1km": d["ln_density_1km"], "is_montreal": d["is_montreal"]}, index=d.index)219 for c in _CATD:220 X[c] = (d["category"] == c).astype(float)221 return sm.add_constant(X.astype(float))222Xm = _X_mix(df); ym = df["lnp"].values223Wm = sm.add_constant(pd.DataFrame({"is_montreal": df["is_montreal"], "ln_reviews": df["ln_reviews"], "ln_guests": df["ln_guests"],224 "is_superhost": df["is_superhost"], "hotel": (df["category"] == "Hotel/B&B/hostel").astype(float),225 "chalet": (df["category"] == "Chalet/cottage/cabin").astype(float)}, index=df.index).astype(float))226r = (df["price"] > 1000).astype(float).values # initialisation: the $1,000 rule227delta, s0, s1 = 3.0, 0.5, 0.8; gamma = None228for it in range(300):229 # M-step: WLS on stacked data for (beta, delta)230 Xs = np.vstack([np.column_stack([Xm.values, np.zeros(len(ym))]), np.column_stack([Xm.values, np.ones(len(ym))])])231 ys = np.concatenate([ym, ym]); ws = np.concatenate([(1 - r) / s0 ** 2, r / s1 ** 2])232 coef = np.linalg.lstsq(Xs * np.sqrt(ws)[:, None], ys * np.sqrt(ws), rcond=None)[0]233 beta, delta = coef[:-1], coef[-1]234 f0 = Xm.values @ beta; f1 = f0 + delta235 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))236 glm = sm.GLM(r, Wm, family=sm.families.Binomial()).fit(); gamma = glm.params.values; pi = glm.fittedvalues.values237 # E-step238 l0 = (1 - pi) * stats.norm.pdf(ym, f0, s0); l1 = pi * stats.norm.pdf(ym, f1, s1)239 r_new = l1 / (l0 + l1)240 ll = np.sum(np.log(l0 + l1))241 if np.max(np.abs(r_new - r)) < 1e-7:242 r = r_new; break243 r = r_new244df["p_stay"] = r245df["stay_regime"] = (r >= 0.5).astype(int)246KEY["mixreg"] = {"iterations": it + 1, "loglik": float(ll), "delta": float(delta), "implied_nights": float(np.exp(delta)),247 "sigma_nightly": float(s0), "sigma_stay": float(s1), "share_stay": float(r.mean()), "n_stay": int(df["stay_regime"].sum()),248 "n_nightly": int((df["stay_regime"] == 0).sum()),249 "gamma": {k: [float(v), float(glm.bse[k]), float(glm.pvalues[k])] for k, v in glm.params.items()},250 "beta_common": {k: float(v) for k, v in zip(Xm.columns, beta)},251 "n_stay_le1000": int(((df["stay_regime"] == 1) & (df["price"] <= 1000)).sum()),252 "n_nightly_gt1000": int(((df["stay_regime"] == 0) & (df["price"] > 1000)).sum()),253 "n_ambiguous_25_75": int(((r > 0.25) & (r < 0.75)).sum()),254 "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)))}255# table: mixture parameters256body = ["\\begin{tabular}{lr}", "\\toprule", "Parameter & Estimate \\\\", "\\midrule",257 f"Intercept shift of the stay regime, $\\delta$ (log points) & {delta:.3f} \\\\",258 f"Implied stay length, $\\exp(\\delta)$ (nights) & {np.exp(delta):.1f} \\\\",259 f"Residual s.d., nightly regime, $\\sigma_0$ & {s0:.3f} \\\\",260 f"Residual s.d., stay regime, $\\sigma_1$ & {s1:.3f} \\\\",261 f"Mean posterior probability of the stay regime & {r.mean():.3f} \\\\",262 f"Listings classified in the stay regime ($\\hat{{p}}_i \\geq 0.5$) & {int(df['stay_regime'].sum()):,} \\\\",263 f"\\quad of which displayed price $\\leq$ \\$1{{,}}000 & {KEY['mixreg']['n_stay_le1000']:,} \\\\",264 f"Listings classified nightly with displayed price $>$ \\$1{{,}}000 & {KEY['mixreg']['n_nightly_gt1000']:,} \\\\",265 f"Listings with $0.25 < \\hat{{p}}_i < 0.75$ & {KEY['mixreg']['n_ambiguous_25_75']:,} \\\\",266 f"Log-likelihood & {ll:,.1f} \\\\", f"EM iterations & {it + 1} \\\\", "\\midrule",267 "\\multicolumn{2}{l}{\\textit{Regime equation (logit coefficients, SE)}} \\\\"]268GL = {"const": "Constant", "is_montreal": "Montreal", "ln_reviews": "ln(1 + reviews)", "ln_guests": "ln(guests)", "is_superhost": "Superhost",269 "hotel": "Hotel/B\\&B/hostel", "chalet": "Chalet/cottage/cabin"}270for k, v in glm.params.items():271 body.append(f"{GL[k]} & {v:.3f}{stars(glm.pvalues[k])} ({glm.bse[k]:.3f}) \\\\")272body += ["\\bottomrule", "\\end{tabular}"]273write_tex("rev_mixture", "\n".join(body))274275# regime logit: which listings are classified in the stay regime?276Xl = pd.get_dummies(df[["ln_bedrooms", "ln_bathrooms", "ln_guests", "ln_amenities", "is_superhost", "is_guest_favorite",277 "rated", "ln_reviews", "is_montreal", "ln_dist", "category"]], columns=["category"], drop_first=True, dtype=float)278Xl.columns = [c.replace("category_", "") for c in Xl.columns]279logit = sm.Logit(df["stay_regime"], sm.add_constant(Xl)).fit(disp=0, cov_type="HC1")280me = logit.get_margeff(at="overall")281LAB = {"ln_bedrooms": "ln(bedrooms)", "ln_bathrooms": "ln(1 + bathrooms)", "ln_guests": "ln(guests)",282 "ln_amenities": "ln(1 + amenities)", "is_superhost": "Superhost", "is_guest_favorite": "Guest Favourite",283 "rated": "Rated (star shown)", "ln_reviews": "ln(1 + reviews)", "is_montreal": "Montreal", "ln_dist": "ln(1 + km to centre)",284 "rating_c": "Rating $-$ 4.8 (rated only)", "ln_images": "ln(1 + photos)", "pets_allowed": "Pets allowed",285 "shared_bath": "Shared bathroom", "studio": "Studio", "ln_density_1km": "ln(1 + listings within 1 km)",286 "Chalet/cottage/cabin": "Chalet/cottage/cabin", "House": "House", "Condo/loft/apartment": "Condo/loft/apartment",287 "Hotel/B&B/hostel": "Hotel/B\\&B/hostel", "Other/unconventional": "Other/unconventional", "const": "Constant"}288body = ["\\begin{tabular}{lrrrr}", "\\toprule", "Variable & Coefficient & (SE) & Avg.\\ marginal effect & (SE) \\\\", "\\midrule"]289me_tab = me.summary_frame()290for v in Xl.columns:291 body.append(f"{LAB.get(v, v)} & {logit.params[v]:.3f}{stars(logit.pvalues[v])} & ({logit.bse[v]:.3f}) & "292 f"{me_tab.loc[v, 'dy/dx']:.3f} & ({me_tab.loc[v, 'Std. Err.']:.3f}) \\\\")293body += ["\\midrule", f"Observations & \\multicolumn{{4}}{{c}}{{{int(logit.nobs):,}}} \\\\",294 f"Pseudo-$R^2$ & \\multicolumn{{4}}{{c}}{{{logit.prsquared:.3f}}} \\\\",295 f"Share in regime & \\multicolumn{{4}}{{c}}{{{df['stay_regime'].mean():.3f}}} \\\\", "\\bottomrule", "\\end{tabular}"]296write_tex("rev_stay_logit", "\n".join(body))297KEY["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}298299# figure: log-price histogram with mixture300fig, ax = plt.subplots(figsize=(6.5, 3.8))301ax.hist(df["lnp"], bins=60, density=True, color="#9ecae1", edgecolor="white")302xs = grid303for k, c in zip(range(2), ("#3182bd", "#de2d26")):304 ax.plot(xs, wt[k] * stats.norm.pdf(xs, mu[k], sd[k]), color=c, lw=1.5,305 label=f"component {k + 1}: mode \\${np.exp(mu[k]):,.0f}, weight {wt[k]:.2f}")306ax.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")307ax.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")308ax.axvline(thr_ln, color="#de2d26", ls=":", lw=1)309ticks = [25, 50, 100, 200, 500, 1000, 2000, 5000, 10000]310ax.set_xticks(np.log(ticks)); ax.set_xticklabels([f"{t:,}" for t in ticks])311ax.set_xlabel("Displayed price (CAD, log scale)"); ax.set_ylabel("Density"); ax.legend(frameon=False, fontsize=8)312savefig(fig, "rev_fig_logprice_mixture")313314# --------------------------------------------------------------------------315# 3. Analysis sample and descriptives316# --------------------------------------------------------------------------317A = df[df["stay_regime"] == 0].copy().reset_index(drop=True) # nightly regime (posterior < 0.5)318KEY["N_A"] = len(A); KEY["N_priced"] = N_PRICED; KEY["N_scraped"] = 5000319KEY["A_price_median"] = float(A.price.median()); KEY["A_price_mean"] = float(A.price.mean())320KEY["A_share_montreal"] = float(A.is_montreal.mean()); KEY["A_n_cities"] = int(A.city.nunique())321KEY["A_unrated_share"] = float((A.rated == 0).mean()); KEY["A_superhost"] = float(A.is_superhost.mean())322KEY["A_guestfav"] = float(A.is_guest_favorite.mean()); KEY["A_cat_shares"] = A.category.value_counts(normalize=True).round(3).to_dict()323KEY["A_cat_shares"] = {str(k): float(v) for k, v in KEY["A_cat_shares"].items()}324KEY["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())325326DESC = [("price", "Displayed price (CAD)"), ("lnp", "ln(price)"), ("bedrooms", "Bedrooms"), ("bathrooms", "Bathrooms"),327 ("guests_count", "Guest capacity"), ("beds", "Beds (from title)"), ("amenities_count", "Amenities count"),328 ("num_images", "Photos"), ("num_reviews", "Reviews"), ("rating_raw", "Rating (rated listings)"),329 ("rated", "Rated (star shown)"), ("is_superhost", "Superhost"), ("is_guest_favorite", "Guest Favourite"),330 ("pets_allowed", "Pets allowed"), ("shared_bath", "Shared bathroom"), ("studio", "Studio"),331 ("dist_centre_km", "Distance to nearest centre (km)"), ("density_1km", "Listings within 1 km"), ("is_montreal", "Montreal"), ("p_stay", "Posterior prob.\\ of stay regime")]332body = ["\\begin{tabular}{lrrrrrrr}", "\\toprule",333 "Variable & Mean & Std.\\ dev. & P10 & Median & P90 & $N$ & Mean, full \\\\", "\\midrule"]334for c, lab in DESC:335 s = A[c].dropna().astype(float); f = df[c].dropna().astype(float)336 d = 0 if c in ("price", "bedrooms", "guests_count", "beds", "amenities_count", "num_images", "num_reviews", "density_1km") else 2337 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} \\\\")338body += ["\\bottomrule", "\\end{tabular}"]339write_tex("rev_summary", "\n".join(body))340341body = ["\\begin{tabular}{lrrrrrr}", "\\toprule", "Title category & $N$ & Share & Mean (\\$) & Median (\\$) & P10 & P90 \\\\", "\\midrule"]342for c in CATS:343 s = A.loc[A.category == c, "price"]344 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} \\\\")345s = A.price346body += ["\\midrule", f"All & {len(s):,} & 1.000 & {s.mean():,.0f} & {s.median():,.0f} & {s.quantile(.1):,.0f} & {s.quantile(.9):,.0f} \\\\",347 "\\bottomrule", "\\end{tabular}"]348write_tex("rev_price_by_category", "\n".join(body))349350# --------------------------------------------------------------------------351# 4. Nested OLS hedonic models352# --------------------------------------------------------------------------353STRUCT = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "studio"]354CATD = [c for c in CATS if c != "Rental unit"]355AMEN = ["ln_amenities", "pets_allowed", "ln_images"]356HOST = ["is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c"]357LOC = ["ln_dist", "ln_density_1km"]358# city fixed effects: cities with >= 15 listings in the analysis sample; the rest pooled by K-means cell359MIN_CITY = 15360from sklearn.cluster import KMeans # noqa: E402361vc = A.city.value_counts()362big = set(vc[vc >= MIN_CITY].index)363small = A[~A.city.isin(big)]364km_small = KMeans(n_clusters=25, random_state=SEED, n_init=10).fit(small[["x_km", "y_km"]].values)365A["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))366KEY["n_city_fe"] = int(len(big)); KEY["n_cell_fe"] = 25; KEY["share_in_city_fe"] = float(A.city.isin(big).mean())367KEY["n_fe_groups"] = int(A.fe_city.nunique())368# finer grid-cell FE for robustness369A["fe_cell60"] = KMeans(n_clusters=60, random_state=SEED, n_init=10).fit_predict(A[["x_km", "y_km"]].values)370371372def design(d, cols, fe=None, drop_first_fe="Montreal"):373 X = pd.DataFrame(index=d.index)374 for c in cols:375 X[c] = (d.category == c).astype(float) if c in CATD else d[c].astype(float)376 if fe is not None:377 D = pd.get_dummies(d[fe].astype(str), prefix="fe", dtype=float)378 ref = f"fe_{drop_first_fe}" if f"fe_{drop_first_fe}" in D.columns else D.columns[0]379 D = D.drop(columns=[ref])380 X = pd.concat([X, D], axis=1)381 return sm.add_constant(X)382383384SPECS = {385 "M1": STRUCT + CATD,386 "M2": STRUCT + CATD + AMEN,387 "M3": STRUCT + CATD + AMEN + HOST,388 "M4": STRUCT + CATD + AMEN + HOST + LOC,389 "M5": STRUCT + CATD + AMEN + HOST + LOC, # + city FE390}391fits = {}392for k, cols in SPECS.items():393 X = design(A, cols, fe="fe_city" if k == "M5" else None)394 fits[k] = sm.OLS(A["lnp"], X).fit(cov_type="HC3")395fit5_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})396fit5_cell = sm.OLS(A["lnp"], design(A, SPECS["M5"], fe="fe_city")).fit(cov_type="cluster", cov_kwds={"groups": A["fe_cell60"]})397# VIF (M4, no FE)398from statsmodels.stats.outliers_influence import variance_inflation_factor # noqa: E402399X4 = design(A, SPECS["M4"])400vif = pd.Series({c: variance_inflation_factor(X4.values, i) for i, c in enumerate(X4.columns) if c != "const"})401KEY["vif_max"] = float(vif.max()); KEY["vif_argmax"] = str(vif.idxmax()); KEY["vif"] = {k: float(v) for k, v in vif.items()}402403ROWS = STRUCT + CATD + AMEN + HOST + LOC404body = ["\\begin{tabular}{l" + "r" * 6 + "}", "\\toprule",405 " & (1) & (2) & (3) & (4) & (5) & (5$'$) \\\\",406 "Dependent variable: ln(price) & M1 & M2 & M3 & M4 & M5 & M5 \\\\", "\\midrule"]407for v in ROWS:408 line = [LAB.get(v, tex_escape(v))]; se_line = [""]409 for k in ("M1", "M2", "M3", "M4", "M5"):410 f = fits[k]411 if v in f.params.index:412 line.append(f"{f.params[v]:.3f}{stars(f.pvalues[v])}"); se_line.append(f"({f.bse[v]:.3f})")413 else:414 line.append(""); se_line.append("")415 if v in fit5_cl.params.index:416 line.append(f"{fit5_cl.params[v]:.3f}{stars(fit5_cl.pvalues[v])}"); se_line.append(f"[{fit5_cl.bse[v]:.3f}]")417 else:418 line.append(""); se_line.append("")419 body.append(" & ".join(line) + " \\\\"); body.append(" & ".join(se_line) + " \\\\")420body += ["\\midrule", "City / cell fixed effects & No & No & No & No & Yes & Yes \\\\",421 "Standard errors & HC3 & HC3 & HC3 & HC3 & HC3 & Cluster \\\\"]422for lab, fn in (("Observations", lambda f: f"{int(f.nobs):,}"), ("$R^2$", lambda f: f"{f.rsquared:.3f}"),423 ("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}"),424 ("AIC", lambda f: f"{f.aic:,.0f}")):425 body.append(lab + " & " + " & ".join(fn(fits[k]) for k in ("M1", "M2", "M3", "M4", "M5")) + f" & {fn(fit5_cl)} \\\\")426body += ["\\bottomrule", "\\end{tabular}"]427write_tex("rev_ols_nested", "\n".join(body))428KEY["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)),429 "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}}430 for k, f in fits.items()}431KEY["ols"]["M5_cluster_city"] = {v: [float(fit5_cl.params[v]), float(fit5_cl.bse[v]), float(fit5_cl.pvalues[v])] for v in ROWS}432KEY["ols"]["M5_cluster_cell60"] = {v: [float(fit5_cell.params[v]), float(fit5_cell.bse[v]), float(fit5_cell.pvalues[v])] for v in ROWS}433# Kennedy (1981) corrected premia for dummies in M5434f = fits["M5"]435KEY["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"]}436# F-tests of blocks (M5)437def ftest(fit, vars_):438 R = np.zeros((len(vars_), len(fit.params)));439 for i, v in enumerate(vars_): R[i, list(fit.params.index).index(v)] = 1440 t = fit.f_test(R); return float(np.squeeze(t.fvalue)), float(t.pvalue)441KEY["ftests_M5"] = {"category": ftest(f, CATD), "host": ftest(f, HOST), "amenities": ftest(f, AMEN), "location": ftest(f, LOC)}442fe_cols = [c for c in f.params.index if c.startswith("fe_")]443KEY["ftests_M5"]["city_fe"] = ftest(f, fe_cols)444KEY["r2_within_fe_only"] = float(sm.OLS(A["lnp"], design(A, [], fe="fe_city")).fit().rsquared)445log(f"OLS done: R2 M1..M5 = {[round(fits[k].rsquared, 3) for k in fits]}")446447# --------------------------------------------------------------------------448# 5. Quantile regression (M5 specification)449# --------------------------------------------------------------------------450TAUS = [0.10, 0.25, 0.50, 0.75, 0.90]451X5 = design(A, SPECS["M5"], fe="fe_city")452qfits = {t: sm.QuantReg(A["lnp"], X5).fit(q=t, max_iter=5000) for t in TAUS}453QROWS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "Chalet/cottage/cabin", "House", "Condo/loft/apartment",454 "Hotel/B&B/hostel", "ln_amenities", "ln_images", "is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c", "ln_dist", "ln_density_1km"]455body = ["\\begin{tabular}{l" + "r" * 6 + "}", "\\toprule",456 "Variable & OLS & $\\tau=0.10$ & $\\tau=0.25$ & $\\tau=0.50$ & $\\tau=0.75$ & $\\tau=0.90$ \\\\", "\\midrule"]457for v in QROWS:458 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})"]459 for t in TAUS:460 q = qfits[t]; line.append(f"{q.params[v]:.3f}{stars(q.pvalues[v])}"); se.append(f"({q.bse[v]:.3f})")461 body.append(" & ".join(line) + " \\\\"); body.append(" & ".join(se) + " \\\\")462463464def pseudo_r2(q, t):465 u = A["lnp"] - q.fittedvalues; rho = lambda u: np.sum(u * (t - (u < 0)))466 u0 = A["lnp"] - A["lnp"].quantile(t); return 1 - rho(u) / rho(u0)467468469body += ["\\midrule", "Pseudo-$R^2$ / $R^2$ & " + f"{fits['M5'].rsquared:.3f} & " + " & ".join(f"{pseudo_r2(qfits[t], t):.3f}" for t in TAUS) + " \\\\",470 f"Observations & \\multicolumn{{6}}{{c}}{{{len(A):,}}} \\\\", "City / cell fixed effects & \\multicolumn{6}{c}{Yes} \\\\",471 "\\bottomrule", "\\end{tabular}"]472write_tex("rev_quantile", "\n".join(body))473KEY["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}474KEY["quantile_pseudoR2"] = {str(t): float(pseudo_r2(qfits[t], t)) for t in TAUS}475# interquantile difference tests (0.90 - 0.10) by pairwise bootstrap (B = 200)476rng = np.random.default_rng(SEED)477B = 100; diffs = {v: [] for v in ("ln_bedrooms", "ln_guests", "is_superhost", "rating_c", "ln_reviews", "ln_dist")}478cache = TAB / "rev_iq_bootstrap.json"479if cache.exists():480 diffs = json.loads(cache.read_text()); log("inter-quantile bootstrap loaded from cache")481else:482 for b in range(B):483 idx = rng.integers(0, len(A), len(A))484 q1 = sm.QuantReg(A["lnp"].values[idx], X5.values[idx]).fit(q=0.10, max_iter=3000)485 q9 = sm.QuantReg(A["lnp"].values[idx], X5.values[idx]).fit(q=0.90, max_iter=3000)486 for v in diffs:487 j = list(X5.columns).index(v); diffs[v].append(float(q9.params[j] - q1.params[j]))488 cache.write_text(json.dumps(diffs))489KEY["iq_bootstrap_B"] = B490KEY["iq_diff_90_10"] = {}491for v, arr in diffs.items():492 arr = np.array(arr); est = qfits[0.90].params[v] - qfits[0.10].params[v]493 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)))))]494log("quantile done")495496fig, axes = plt.subplots(2, 3, figsize=(9.5, 5.2)); axes = axes.ravel()497for ax, v in zip(axes, ["ln_bedrooms", "ln_guests", "ln_bathrooms", "is_superhost", "rating_c", "ln_dist"]):498 est = [qfits[t].params[v] for t in TAUS]; se = [qfits[t].bse[v] for t in TAUS]499 ax.fill_between(TAUS, np.array(est) - 1.96 * np.array(se), np.array(est) + 1.96 * np.array(se), color="#c6dbef")500 ax.plot(TAUS, est, "o-", color="#08519c", ms=3)501 ax.axhline(fits["M5"].params[v], color="#de2d26", ls="--", lw=1); ax.axhline(0, color="grey", lw=0.5)502 ax.set_title(LAB[v], fontsize=9); ax.set_xticks(TAUS)503fig.supxlabel("Quantile $\\tau$ of ln(price)"); fig.supylabel("Coefficient")504fig.tight_layout(); savefig(fig, "rev_fig_quantile")505506# --------------------------------------------------------------------------507# 6. Spatial econometrics (M4 specification; W = k-nearest neighbours on projected km)508# --------------------------------------------------------------------------509import libpysal, spreg # noqa: E402510from esda.moran import Moran # noqa: E402511coords = A[["x_km", "y_km"]].values512513514def knn_w(k):515 w = libpysal.weights.KNN.from_array(coords, k=k); w.transform = "r"; return w516517518W5, W10 = knn_w(5), knn_w(10)519KEY["moran"] = {}520for lab, f, w in (("M4_k5", fits["M4"], W5), ("M5_k5", fits["M5"], W5), ("lnp_k5", None, W5), ("M4_k10", fits["M4"], W10)):521 z = A["lnp"].values if f is None else f.resid.values522 mi = Moran(z, w, permutations=999)523 KEY["moran"][lab] = [float(mi.I), float(mi.z_norm), float(mi.p_norm)]524Xs = design(A, SPECS["M4"]).drop(columns="const")525y = A[["lnp"]].values; Xv = Xs.values; names = list(Xs.columns)526sar = spreg.ML_Lag(y, Xv, W5, name_x=names, name_y="lnp", method="full")527sem = spreg.ML_Error(y, Xv, W5, name_x=names, name_y="lnp", method="full")528sdm = spreg.ML_Lag(y, Xv, W5, name_x=names, name_y="lnp", method="full", slx_lags=1)529sar10 = spreg.ML_Lag(y, Xv, W10, name_x=names, name_y="lnp", method="full")530ols4 = fits["M4"]531532533def sp_row(mod, v, lagged=False):534 nm_ = [n for n in mod.name_x]535 try:536 j = nm_.index(("W_" + v) if lagged else v)537 except ValueError:538 return "", ""539 b = float(mod.betas[j, 0]); se_ = float(mod.std_err[j]); p = float(mod.z_stat[j][1])540 return f"{b:.3f}{stars(p)}", f"({se_:.3f})"541542543SROWS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "shared_bath", "Chalet/cottage/cabin", "House", "Condo/loft/apartment", "Hotel/B&B/hostel",544 "ln_amenities", "ln_images", "is_superhost", "is_guest_favorite", "ln_reviews", "rated", "rating_c", "ln_dist", "ln_density_1km"]545body = ["\\begin{tabular}{lrrrrr}", "\\toprule", "Variable & OLS (M4) & SAR & SEM & SDM & SDM, $\\mathbf{W}x$ \\\\", "\\midrule"]546for v in SROWS:547 l = [LAB.get(v, tex_escape(v)), f"{ols4.params[v]:.3f}{stars(ols4.pvalues[v])}"]; s = ["", f"({ols4.bse[v]:.3f})"]548 for m in (sar, sem, sdm):549 a, b_ = sp_row(m, v); l.append(a); s.append(b_)550 a, b_ = sp_row(sdm, v, lagged=True); l.append(a); s.append(b_)551 body.append(" & ".join(l) + " \\\\"); body.append(" & ".join(s) + " \\\\")552rho_p = lambda m: float(m.z_stat[-1][1])553body += ["\\midrule",554 f"$\\rho$ (spatial lag) & & {float(np.squeeze(sar.rho)):.3f}{stars(rho_p(sar))} & & {float(np.squeeze(sdm.rho)):.3f}{stars(rho_p(sdm))} & \\\\",555 f" & & ({float(sar.std_err[-1]):.3f}) & & ({float(sdm.std_err[-1]):.3f}) & \\\\",556 f"$\\lambda$ (spatial error) & & & {float(np.squeeze(sem.lam)):.3f}{stars(rho_p(sem))} & & \\\\",557 f" & & & ({float(sem.std_err[-1]):.3f}) & & \\\\",558 f"Log-likelihood & {ols4.llf:,.1f} & {float(np.squeeze(sar.logll)):,.1f} & {float(np.squeeze(sem.logll)):,.1f} & {float(np.squeeze(sdm.logll)):,.1f} & \\\\",559 f"AIC & {ols4.aic:,.0f} & {float(np.squeeze(sar.aic)):,.0f} & {float(np.squeeze(sem.aic)):,.0f} & {float(np.squeeze(sdm.aic)):,.0f} & \\\\",560 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} & \\\\",561 f"Observations & \\multicolumn{{5}}{{c}}{{{len(A):,}}} \\\\", "\\bottomrule", "\\end{tabular}"]562write_tex("rev_spatial", "\n".join(body))563KEY["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)],564 "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)],565 "ll": {"ols": float(ols4.llf), "sar": float(np.squeeze(sar.logll)), "sem": float(np.squeeze(sem.logll)), "sdm": float(np.squeeze(sdm.logll))},566 "aic": {"ols": float(ols4.aic), "sar": float(np.squeeze(sar.aic)), "sem": float(np.squeeze(sem.aic)), "sdm": float(np.squeeze(sdm.aic))}}567# LR test SDM vs SAR (restriction theta = 0)568lr = 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)))]569# LeSage-Pace impacts (SAR and SDM), dense inverse570Wd = W5.full()[0]571I_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)572imp = {}573for v in SROWS:574 j = names.index(v); b_sar = float(sar.betas[j + 1, 0])575 d = np.trace(S_sar) / len(A) * b_sar; t = S_sar.sum() / len(A) * b_sar576 b_sdm = float(sdm.betas[j + 1, 0]); th = float(sdm.betas[sdm.name_x.index("W_" + v), 0])577 M = S_sdm @ (I_n * b_sdm + Wd * th)578 d2 = np.trace(M) / len(A); t2 = M.sum() / len(A)579 imp[v] = {"sar": [d, t - d, t], "sdm": [d2, t2 - d2, t2]}580# simulated SEs for SAR and SDM impacts (300 draws from the ML asymptotic distribution; exact traces via the581# complex eigenvalues of the asymmetric KNN matrix: tr[(I - rW)^-1] = sum_i 1/(1 - r*lambda_i))582rng = np.random.default_rng(SEED)583evals = np.linalg.eigvals(Wd)584def _sim(mod, has_theta):585 V = np.asarray(mod.vm); th_hat = mod.betas.ravel(); dr_all = rng.multivariate_normal(th_hat, V, size=300)586 out = {v: [] for v in SROWS}587 for dr in dr_all:588 r_ = dr[-1]589 if abs(r_) >= 1:590 continue591 tr = np.mean(1 / (1 - r_ * evals)).real; trw = np.mean(evals / (1 - r_ * evals)).real; tot = 1 / (1 - r_)592 for v in SROWS:593 b = dr[names.index(v) + 1]; th = dr[mod.name_x.index("W_" + v)] if has_theta else 0.0594 direct = tr * b + trw * th; total = (b + th) * tot595 out[v].append([direct, total - direct, total])596 return {v: np.array(a).std(axis=0, ddof=1) for v, a in out.items()}597sd_sar = _sim(sar, False); sd_sdm = _sim(sdm, True)598pz = lambda est, sd_: stars(2 * (1 - stats.norm.cdf(abs(est / sd_)))) if sd_ > 0 else ""599body = ["\\begin{tabular}{lrrrrrr}", "\\toprule", " & \\multicolumn{3}{c}{SAR} & \\multicolumn{3}{c}{SDM} \\\\",600 "\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}", "Variable & Direct & Indirect & Total & Direct & Indirect & Total \\\\", "\\midrule"]601KEY["impacts"] = {}602for v in SROWS:603 a = imp[v]; s1_ = sd_sar[v]; s2_ = sd_sdm[v]604 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)) + " & "605 + " & ".join(f"{a['sdm'][i]:.3f}{pz(a['sdm'][i], s2_[i])}" for i in range(3)) + " \\\\")606 body.append(" & " + " & ".join(f"({s1_[i]:.3f})" for i in range(3)) + " & " + " & ".join(f"({s2_[i]:.3f})" for i in range(3)) + " \\\\")607 KEY["impacts"][v] = {"sar": a["sar"], "sar_se": s1_.tolist(), "sdm": a["sdm"], "sdm_se": s2_.tolist()}608body += ["\\bottomrule", "\\end{tabular}"]609write_tex("rev_impacts", "\n".join(body))610KEY["spatial"]["sar_multiplier"] = float(S_sar.sum() / len(A)); KEY["spatial"]["sar_direct_ratio"] = float(np.trace(S_sar) / len(A))611log("spatial done")612613# --------------------------------------------------------------------------614# 7. Machine learning: random vs spatially blocked cross-validation, SHAP615# --------------------------------------------------------------------------616from sklearn.model_selection import KFold, GroupKFold # noqa: E402617from sklearn.linear_model import ElasticNetCV, LinearRegression # noqa: E402618from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor # noqa: E402619from sklearn.preprocessing import StandardScaler # noqa: E402620from sklearn.pipeline import make_pipeline # noqa: E402621import xgboost as xgb, lightgbm as lgb # noqa: E402622623XML = design(A, SPECS["M5"], fe="fe_city").drop(columns="const")624XML.columns = [re.sub(r"[^0-9a-zA-Z_]", "_", c) for c in XML.columns]625yv = A["lnp"].values; pv = A["price"].values626blocks = KMeans(n_clusters=20, random_state=SEED, n_init=10).fit_predict(coords)627KEY["ml_n_blocks"] = 20628MODELS = {629 "OLS (M5)": lambda: LinearRegression(),630 "Elastic net": lambda: make_pipeline(StandardScaler(), ElasticNetCV(l1_ratio=[.1, .5, .9], cv=5, random_state=SEED, n_jobs=-1)),631 "Random forest": lambda: RandomForestRegressor(n_estimators=500, min_samples_leaf=3, max_features=0.5, random_state=SEED, n_jobs=-1),632 "Gradient boosting": lambda: GradientBoostingRegressor(n_estimators=500, learning_rate=0.05, max_depth=3, subsample=0.8, random_state=SEED),633 "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),634 "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),635}636637638def cv_eval(splitter, groups=None):639 out = {}640 for name, mk in MODELS.items():641 yhat = np.zeros_like(yv)642 for tr, te in splitter.split(XML, yv, groups):643 m = mk(); m.fit(XML.iloc[tr], yv[tr]); yhat[te] = m.predict(XML.iloc[te])644 r2 = 1 - np.sum((yv - yhat) ** 2) / np.sum((yv - yv.mean()) ** 2)645 rmse = np.sqrt(np.mean((yv - yhat) ** 2)); mae = np.mean(np.abs(yv - yhat))646 mdape = np.median(np.abs(np.exp(yhat) - pv) / pv) * 100647 out[name] = [r2, rmse, mae, mdape]648 return out649650651res_rand = cv_eval(KFold(5, shuffle=True, random_state=SEED))652res_block = cv_eval(GroupKFold(5), groups=blocks)653body = ["\\begin{tabular}{lrrrrrrrr}", "\\toprule", " & \\multicolumn{4}{c}{Random 5-fold CV} & \\multicolumn{4}{c}{Spatially blocked 5-fold CV} \\\\",654 "\\cmidrule(lr){2-5}\\cmidrule(lr){6-9}", "Model & $R^2$ & RMSE & MAE & MdAPE (\\%) & $R^2$ & RMSE & MAE & MdAPE (\\%) \\\\", "\\midrule"]655for name in MODELS:656 a, b_ = res_rand[name], res_block[name]657 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} \\\\")658body += ["\\bottomrule", "\\end{tabular}"]659write_tex("rev_ml_cv", "\n".join(body))660KEY["ml"] = {"random": res_rand, "blocked": res_block}661log("ML CV done")662663# SHAP on XGBoost fitted to the full analysis sample664import shap # noqa: E402665xg = MODELS["XGBoost"](); xg.fit(XML, yv)666# exact TreeSHAP values computed natively by XGBoost (pred_contribs); last column is the bias term667contrib = xg.get_booster().predict(xgb.DMatrix(XML), pred_contribs=True)668sv = contrib[:, :-1]; KEY['shap_bias'] = float(contrib[0, -1])669assert np.allclose(contrib.sum(axis=1), xg.predict(XML), atol=1e-3)670imp_shap = pd.Series(np.abs(sv).mean(axis=0), index=XML.columns)671fe_mask = imp_shap.index.str.startswith("fe_")672imp_agg = imp_shap[~fe_mask].copy(); imp_agg["City / cell fixed effects (all)"] = imp_shap[fe_mask].sum()673imp_agg = imp_agg.sort_values(ascending=False)674LAB2 = {re.sub(r"[^0-9a-zA-Z_]", "_", k): v for k, v in LAB.items()}675body = ["\\begin{tabular}{lrr}", "\\toprule", "Feature & Mean $|\\phi_k|$ & Share of total (\\%) \\\\", "\\midrule"]676tot = imp_shap.sum()677for k, v in imp_agg.head(16).items():678 body.append(f"{LAB2.get(k, tex_escape(k))} & {v:.4f} & {100 * v / tot:.1f} \\\\")679body += ["\\bottomrule", "\\end{tabular}"]680write_tex("rev_shap_importance", "\n".join(body))681KEY["shap_top"] = {k: float(v) for k, v in imp_agg.head(10).items()}; KEY["shap_total"] = float(tot)682KEY["shap_fe_share"] = float(imp_shap[fe_mask].sum() / tot)683# beeswarm of the 12 most important non-FE features684top = imp_shap[~fe_mask].sort_values(ascending=False).head(12).index685Xd = XML[top].copy(); Xd.columns = [LAB2.get(c, c) for c in top]686fig = plt.figure(figsize=(7, 4.6))687shap.summary_plot(sv[:, [list(XML.columns).index(c) for c in top]], Xd, show=False, plot_size=None, max_display=12)688plt.gcf().set_size_inches(7, 4.6); plt.xlabel("SHAP value (contribution to predicted ln price)")689savefig(plt.gcf(), "rev_fig_shap")690log("SHAP done")691692# --------------------------------------------------------------------------693# 8. Robustness of the M5 specification694# --------------------------------------------------------------------------695RKEYS = ["ln_bedrooms", "ln_bathrooms", "ln_guests", "is_superhost", "is_guest_favorite", "ln_reviews", "rating_c", "ln_dist"]696697698def refit(d, fe="fe_city", label=""):699 d = d.reset_index(drop=True)700 if fe == "fe_cell60":701 X = design(d, SPECS["M5"], fe="fe_cell60", drop_first_fe="0")702 else:703 X = design(d, SPECS["M5"], fe=fe)704 X = X.loc[:, (X != 0).any(axis=0)]705 if label.startswith("Weighted"):706 return sm.WLS(d["lnp"], X, weights=1 - d["p_stay"]).fit(cov_type="HC3")707 return sm.OLS(d["lnp"], X).fit(cov_type="HC3")708709710def _fe(d):711 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)))712variants = [713 ("Baseline: nightly regime ($\\hat{p}_i < 0.5$)", A),714 ("Nightly regime, strict ($\\hat{p}_i < 0.25$)", _fe(df[df.p_stay < 0.25])),715 ("Nightly regime, lenient ($\\hat{p}_i < 0.75$)", _fe(df[df.p_stay < 0.75])),716 ("Price $\\leq$ \\$1{,}000 rule", _fe(df[df.price <= 1000])),717 ("Price $\\leq$ \\$600 rule", _fe(df[df.price <= 600])),718 ("All priced listings (bimodal)", _fe(df)),719 ("Excluding hotels, B\\&Bs, hostels", A[A.category != "Hotel/B&B/hostel"]),720 ("Rated listings only", A[A.rated == 1]),721 ("Montreal only", A[A.is_montreal == 1]),722 ("Outside Montreal", A[A.is_montreal == 0]),723 ("Trimmed 1\\%/99\\% within sample", A[(A.price > A.price.quantile(.01)) & (A.price < A.price.quantile(.99))]),724 ("Weighted by $1 - \\hat{p}_i$", A),725]726body = ["\\begin{tabular}{l" + "r" * (len(RKEYS) + 2) + "}", "\\toprule",727 "Sample / specification & " + " & ".join(["Bedr.", "Bathr.", "Guests", "Superhost", "Guest Fav.", "ln reviews", "Rating", "ln dist."]) + " & $N$ & $R^2$ \\\\", "\\midrule"]728KEY["robust"] = {}729for lab, d in variants:730 f = refit(d, label=lab)731 cells = [f"{f.params[v]:.3f}{stars(f.pvalues[v])}" if v in f.params.index else "" for v in RKEYS]732 ses = [f"({f.bse[v]:.3f})" if v in f.params.index else "" for v in RKEYS]733 body.append(f"{lab} & " + " & ".join(cells) + f" & {int(f.nobs):,} & {f.rsquared:.3f} \\\\")734 body.append(" & " + " & ".join(ses) + " & & \\\\")735 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}}736# grid-cell FE (60 cells)737f = refit(A, fe="fe_cell60")738body.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} \\\\")739body.append(" & " + " & ".join(f"({f.bse[v]:.3f})" for v in RKEYS) + " & & \\\\")740KEY["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}}741# Montreal x bedrooms heterogeneity742body += ["\\bottomrule", "\\end{tabular}"]743write_tex("rev_robustness", "\n".join(body))744log("robustness done")745746# --------------------------------------------------------------------------747# 9. Figures: map, price by category, binned partial relationships748# --------------------------------------------------------------------------749fig, ax = plt.subplots(figsize=(7, 5.2))750sc = ax.scatter(A["long"], A["lat"], c=A["lnp"], s=5, cmap="viridis", alpha=0.7, linewidths=0)751cb = fig.colorbar(sc, ax=ax, shrink=0.8); cb.set_label("Displayed price (CAD)")752tk = [50, 100, 200, 400, 800]; cb.set_ticks(np.log(tk)); cb.set_ticklabels([str(t) for t in tk])753ax.set_xlabel("Longitude"); ax.set_ylabel("Latitude"); ax.set_xlim(-80, -61); ax.set_ylim(44.8, 50.5)754for nmc, (la, lo) in centres.items():755 if nmc in ("Montréal", "Québec", "Gatineau", "Mont-Tremblant", "Sherbrooke", "Saguenay", "Rimouski"):756 ax.annotate(nmc, (lo, la), fontsize=7, xytext=(3, 3), textcoords="offset points")757savefig(fig, "rev_fig_map")758759fig, ax = plt.subplots(figsize=(7, 3.8))760data = [A.loc[A.category == c, "price"] for c in CATS]761ax.boxplot(data, vert=False, showfliers=False, widths=0.6, medianprops={"color": "#de2d26"})762ax.set_yticks(range(1, len(CATS) + 1)); ax.set_yticklabels([f"{c} (n={len(d)})" for c, d in zip(CATS, data)], fontsize=8)763ax.set_xscale("log"); ax.set_xlabel("Displayed price (CAD, log scale)"); savefig(fig, "rev_fig_price_by_category")764765fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))766r4 = fits["M4"]767# partial residual plots: residual + beta*x against x, binned768for ax, v, xl in zip(axes, ("ln_dist", "rating_c"), ("Distance to nearest centre (km)", "Rating (rated listings)")):769 d = A.copy(); d["pr"] = r4.resid + r4.params[v] * d[v]770 if v == "rating_c":771 d = d[d.rated == 1]; xvar = d["rating_raw"]772 else:773 xvar = d["dist_centre_km"]774 bins = pd.qcut(xvar, 20, duplicates="drop"); g = d.groupby(bins, observed=True)775 xm = g.apply(lambda s: xvar.loc[s.index].mean()); ym = g["pr"].mean(); se = g["pr"].sem()776 ax.errorbar(xm, ym, yerr=1.96 * se, fmt="o", ms=3, color="#08519c", ecolor="#9ecae1")777 ax.set_xlabel(xl); ax.set_ylabel("Partial residual of ln(price)")778 if v == "ln_dist": ax.set_xscale("symlog", linthresh=1)779fig.tight_layout(); savefig(fig, "rev_fig_partials")780781# --------------------------------------------------------------------------782# 10. Save key numbers783# --------------------------------------------------------------------------784KEY["runtime_s"] = time.perf_counter() - T0785(TAB / "rev_key_numbers.json").write_text(json.dumps(KEY, indent=1, default=float), encoding="utf-8")786A.drop(columns=["name"]).to_csv(ROOT / "data" / "processed" / "analysis_sample.csv", index=False)787log("done")788