UQO Working Paper No. 5 — Airbnb, residential rents and housing market pressure.
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""307_spatial_models.py4--------------------5Spatial regression analysis for the Airbnb-rent study.67Model 4: Spatial models8 - If libpysal + spreg available: SAR, SEM, and OLS comparison9 - Fallback: manual spatial-lag OLS ("poor man's spatial model")10 - Distance-buffer robustness (250m, 500m, 1km, 2km)1112Outputs:13 results/tables/spatial_models.tex14 results/tables/buffer_robustness.tex15 figures/coefficient_buffer_comparison.pdf16"""1718import sys19from pathlib import Path2021sys.path.insert(0, str(Path(__file__).resolve().parent.parent))2223import matplotlib24matplotlib.use("Agg")25import matplotlib.pyplot as plt # noqa: E40226import numpy as np # noqa: E40227import pandas as pd # noqa: E40228import scipy.stats # noqa: E40229import statsmodels.api as sm # noqa: E40230from scipy.spatial import cKDTree # noqa: E4023132from src.config import MERGED_ANALYSIS, TABLE_DIR, FIG_DIR, require # noqa: E40233from src.latex_tables import significance_star as _star, label # noqa: E40234import json # noqa: E40235from scipy.sparse import identity as sp_identity # noqa: E40236from scipy.sparse.linalg import splu # noqa: E4023738PROCESSED_HINT = "Run scripts/04_merge_data.py first."3940# ── attempt to import spatial libraries ──────────────────────────────────────41try:42 from libpysal.weights import KNN as Wknn43 import spreg44 HAS_SPATIAL = True45 print(" libpysal + spreg available — will run SAR / SEM")46except ImportError:47 HAS_SPATIAL = False48 print(" libpysal/spreg not available — using fallback spatial approach")495051def project_km(lat: np.ndarray, lon: np.ndarray) -> np.ndarray:52 """53 Equirectangular projection to kilometres around the sample centroid.5455 2026-10-17 correction: the k-nearest-neighbour search was previously run56 on raw (latitude, longitude) degrees. At 45-47 degrees North one degree of57 longitude spans only ~0.69 of a degree of latitude, so the neighbourhoods58 were anisotropic (stretched east-west). Distances are now in kilometres.59 """60 lat0 = np.radians(np.nanmean(lat))61 x = (lon - np.nanmean(lon)) * 111.32 * np.cos(lat0)62 y = (lat - np.nanmean(lat)) * 111.3263 return np.column_stack([x, y])646566def lesage_pace_impacts(w_sparse, rho: float, beta: dict, vm: np.ndarray,67 idx: dict, n_draws: int = 2000, seed: int = 42) -> dict:68 """69 Average direct, indirect and total impacts of the SAR model70 (LeSage and Pace, 2009, ch. 2) for a row-standardised W.7172 total = beta / (1 - rho) (exact for row-standardised W)73 direct = beta * (1/n) tr[(I - rho W)^{-1}]74 indirect = total - direct7576 The trace is evaluated exactly at the point estimate with a sparse LU77 factorisation, and as a power series tr[(I - rho W)^{-1}] = sum_k rho^k tr(W^k)78 inside a Monte-Carlo loop over draws of (beta, rho) from their asymptotic79 normal distribution to obtain simulated standard errors. tr(W^k) for k >= 280 is estimated with Hutchinson's stochastic trace estimator (50 Rademacher81 vectors); tr(W^0) = n and tr(W) = 0 are exact.82 """83 n = w_sparse.shape[0]84 W = w_sparse.tocsc()85 # exact mean diagonal at the point estimate86 lu = splu((sp_identity(n, format="csc") - rho * W).tocsc())87 diag_sum = 0.088 chunk = 50089 for start in range(0, n, chunk):90 cols = np.arange(start, min(start + chunk, n))91 rhs = np.zeros((n, len(cols)))92 rhs[cols, np.arange(len(cols))] = 1.093 sol = lu.solve(rhs)94 diag_sum += sol[cols, np.arange(len(cols))].sum()95 mean_diag_exact = diag_sum / n9697 # Hutchinson traces for the power series98 rng = np.random.default_rng(seed)99 K = 60100 traces = np.zeros(K + 1)101 traces[0] = n102 m = 50103 V = rng.choice([-1.0, 1.0], size=(n, m))104 Wk = V.copy()105 for k in range(1, K + 1):106 Wk = W @ Wk107 traces[k] = np.einsum("ij,ij->", V, Wk) / m108 traces[1] = 0.0 # zero diagonal by construction109110 def mean_diag(r):111 return np.sum((r ** np.arange(K + 1)) * traces) / n112113 # Monte-Carlo over (beta_j, rho)114 names = list(beta.keys())115 j_rho = idx["rho"]116 sel = [idx[v] for v in names] + [j_rho]117 mu = np.array([beta[v] for v in names] + [rho])118 cov = vm[np.ix_(sel, sel)]119 draws = rng.multivariate_normal(mu, cov, size=n_draws)120 out = {"rho": rho, "mean_diag_exact": mean_diag_exact,121 "mean_diag_series": mean_diag(rho)}122 for i, v in enumerate(names):123 b = beta[v]124 direct = b * mean_diag_exact125 total = b / (1 - rho)126 d_sim = np.array([dr[i] * mean_diag(dr[-1]) for dr in draws])127 t_sim = draws[:, i] / (1 - draws[:, -1])128 out[v] = {129 "direct": direct, "indirect": total - direct, "total": total,130 "direct_se": d_sim.std(ddof=1), "indirect_se": (t_sim - d_sim).std(ddof=1),131 "total_se": t_sim.std(ddof=1),132 }133 return out134135136def run_ols_robust(df, y_col, x_cols, label):137 """OLS with HC1 robust SE."""138 sub = df[[y_col] + x_cols].dropna()139 Y = sub[y_col]140 X = sm.add_constant(sub[x_cols])141 res = sm.OLS(Y, X).fit(cov_type="HC1")142 print(f" [{label}] N={int(res.nobs):,} R2={res.rsquared:.4f}")143 return res144145146def add_manual_spatial_lags(rent_sp: pd.DataFrame, coords: np.ndarray) -> None:147 """Add k=5 nearest-neighbour spatial lags of rent and Airbnb count in place."""148 tree = cKDTree(coords)149 _, idx = tree.query(coords, k=6) # k+1 because first neighbour is self150 idx_neighbours = idx[:, 1:] # drop self151 rent_sp["spatial_lag_rent"] = rent_sp["log_rent"].values[idx_neighbours].mean(axis=1)152 rent_sp["spatial_lag_airbnb"] = rent_sp["airbnb_count_500m"].values[idx_neighbours].mean(axis=1)153154155def spatial_models_table(rent_sp, coords, res_ols, all_x_baseline,156 controls, city_fe_cols, bt_cols, display_vars):157 """158 Estimate the spatial models (SAR/SEM when spreg is available, manual159 spatial-lag OLS otherwise) and return the LaTeX table lines.160 """161 if HAS_SPATIAL:162 # ── Build KNN weights ────────────────────────────────────────────────163 print("\n Building KNN(k=5) spatial weights ...")164 w = Wknn.from_array(coords, k=5)165 w.transform = "r" # row-standardise166167 Y_arr = rent_sp["log_rent"].values.reshape(-1, 1)168 X_arr = sm.add_constant(rent_sp[all_x_baseline].values)169 var_names = ["const"] + all_x_baseline170171 # ── SAR (Spatial Lag Model) ──────────────────────────────────────────172 sar = None173 sem = None174 try:175 print(" Estimating Spatial Lag Model (SAR) ...")176 # 2026-10-17: heteroskedasticity-robust (White) standard errors,177 # consistent with the HC1 inference used elsewhere in the paper.178 sar = spreg.GM_Lag(179 Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names,180 robust="white",181 )182 print(f" [SAR] N={sar.n} pseudo-R2={sar.pr2:.4f} rho={sar.betas[-1][0]:.4f}")183 except Exception as e:184 print(f" [SAR] Failed: {e}")185186 # ── SEM (Spatial Error Model) ────────────────────────────────────────187 try:188 print(" Estimating Spatial Error Model (SEM) ...")189 sem = spreg.GM_Error(190 Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names191 )192 sem_lambda = sem.betas[-1][0]193 print(f" [SEM] N={sem.n} pseudo-R2={sem.pr2:.4f} lambda={sem_lambda:.4f}")194 except Exception as e:195 print(f" [SEM] Failed: {e}")196197 # ── If SAR/SEM failed, fall back to manual spatial-lag approach ─────198 if sar is None or sem is None:199 print(" SAR/SEM failed — falling back to manual spatial lag approach ...")200 add_manual_spatial_lags(rent_sp, coords)201 fb_x = all_x_baseline + ["spatial_lag_rent", "spatial_lag_airbnb"]202 res_fb = run_ols_robust(rent_sp, "log_rent", fb_x, "OLS + spatial lags")203 lines: list[str] = []204 lines.append(r"\begin{tabular}{lcc}")205 lines.append(r"\toprule")206 lines.append(r" & \textbf{OLS Baseline} & \textbf{OLS + Spatial Lags} \\")207 lines.append(r"\midrule")208 ols_p = dict(zip(res_ols.model.exog_names, res_ols.params))209 ols_s = dict(zip(res_ols.model.exog_names, res_ols.bse))210 ols_pv = dict(zip(res_ols.model.exog_names, res_ols.pvalues))211 fb_p = dict(zip(res_fb.model.exog_names, res_fb.params))212 fb_s = dict(zip(res_fb.model.exog_names, res_fb.bse))213 fb_pv = dict(zip(res_fb.model.exog_names, res_fb.pvalues))214 show = (["const", "airbnb_count_500m", "bedrooms", "bathrooms"]215 + bt_cols + ["spatial_lag_rent", "spatial_lag_airbnb"])216 for var in show:217 c1 = f"{ols_p[var]:.4f}{_star(ols_pv[var])}" if var in ols_p else ""218 s1 = f"({ols_s[var]:.4f})" if var in ols_s else ""219 c2 = f"{fb_p[var]:.4f}{_star(fb_pv[var])}" if var in fb_p else ""220 s2 = f"({fb_s[var]:.4f})" if var in fb_s else ""221 vn = label(var)222 lines.append(f"{vn} & {c1} & {c2}" + r" \\")223 lines.append(f" & {s1} & {s2}" + r" \\[4pt]")224 lines.append(r"\midrule")225 lines.append(f"Observations & {int(res_ols.nobs)} & {int(res_fb.nobs)}" + r" \\")226 lines.append(f"R$^2$ & {res_ols.rsquared:.4f} & {res_fb.rsquared:.4f}" + r" \\")227 lines.append(r"\bottomrule")228 lines.append(r"\end{tabular}")229 lines.append(r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in parentheses. $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$. City FE included but not shown.}")230 return lines231232 # ── Build LaTeX table: OLS vs SAR vs SEM ─────────────────────────────233 lines = []234 lines.append(r"\begin{tabular}{lccc}")235 lines.append(r"\toprule")236 lines.append(r" & \textbf{OLS} & \textbf{SAR (GM\_Lag)} & \textbf{SEM (GM\_Error)} \\")237 lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log rent}} \\")238 lines.append(r"\midrule")239 ols_params = dict(zip(["const"] + all_x_baseline, res_ols.params))240 ols_se = dict(zip(["const"] + all_x_baseline, res_ols.bse))241 ols_pval = dict(zip(["const"] + all_x_baseline, res_ols.pvalues))242 sar_var_names = var_names + ["W_log_rent"]243 sar_params = {v: sar.betas[i][0] for i, v in enumerate(sar_var_names)}244 sar_se = {v: sar.std_err[i] for i, v in enumerate(sar_var_names)}245 sar_z = {v: sar.z_stat[i] for i, v in enumerate(sar_var_names)}246 sem_params = {v: sem.betas[i][0] for i, v in enumerate(var_names)}247 sem_se = {v: sem.std_err[i] for i, v in enumerate(var_names)}248 sem_z = {v: sem.z_stat[i] for i, v in enumerate(var_names)}249 show_vars = display_vars + ["W_log_rent"]250 for var in show_vars:251 cells_c, cells_s = [], []252 if var in ols_params:253 cells_c.append(f"{ols_params[var]:.4f}{_star(ols_pval[var])}")254 cells_s.append(f"({ols_se[var]:.4f})")255 else:256 cells_c.append("")257 cells_s.append("")258 if var in sar_params:259 p_sar = 2 * (1 - scipy.stats.norm.cdf(abs(sar_z[var][0])))260 cells_c.append(f"{sar_params[var]:.4f}{_star(p_sar)}")261 cells_s.append(f"({sar_se[var]:.4f})")262 else:263 cells_c.append("")264 cells_s.append("")265 if var in sem_params:266 p_sem = 2 * (1 - scipy.stats.norm.cdf(abs(sem_z[var][0])))267 cells_c.append(f"{sem_params[var]:.4f}{_star(p_sem)}")268 cells_s.append(f"({sem_se[var]:.4f})")269 else:270 cells_c.append("")271 cells_s.append("")272 vn = label(var)273 lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\")274 lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]")275 lines.append(r"$\lambda$ (spatial error) & & & " + f"{sem.betas[-1][0]:.4f}" + r" \\")276277 # ── LeSage–Pace impacts for the SAR model (added 2026-10-17) ─────────278 print(" Computing LeSage-Pace impacts ...")279 idx = {v: i for i, v in enumerate(sar_var_names)}280 idx["rho"] = idx["W_log_rent"]281 imp_vars = ["airbnb_count_500m", "bedrooms", "bathrooms"]282 imp = lesage_pace_impacts(283 w.sparse, sar_params["W_log_rent"],284 {v: sar_params[v] for v in imp_vars}, np.asarray(sar.vm), idx,285 )286 print(f" rho={imp['rho']:.4f} mean diag exact={imp['mean_diag_exact']:.5f} "287 f"series={imp['mean_diag_series']:.5f}")288 ilines = [r"\begin{tabular}{lccc}", r"\toprule",289 r" & \textbf{Direct} & \textbf{Indirect} & \textbf{Total} \\", r"\midrule"]290 for v in imp_vars:291 r = imp[v]292 cells = []293 for key in ("direct", "indirect", "total"):294 est, se = r[key], r[f"{key}_se"]295 pz = 2 * (1 - scipy.stats.norm.cdf(abs(est / se))) if se > 0 else 1.0296 cells.append(f"{est:.4f}{_star(pz)}")297 ilines.append(f"{label(v)} & " + " & ".join(cells) + r" \\")298 ilines.append(" & " + " & ".join(f"({r[k + '_se']:.4f})" for k in ("direct", "indirect", "total")) + r" \\[4pt]")299 ilines += [r"\midrule",300 r"$\hat\rho$ & \multicolumn{3}{c}{" + f"{imp['rho']:.4f}" + r"} \\",301 r"$n^{-1}\,\mathrm{tr}[(I-\hat\rho W)^{-1}]$ & \multicolumn{3}{c}{" + f"{imp['mean_diag_exact']:.4f}" + r"} \\",302 r"Observations & \multicolumn{3}{c}{" + f"{sar.n:,}" + r"} \\",303 r"\bottomrule", r"\end{tabular}",304 r"\parbox{\textwidth}{\footnotesize Average direct, indirect and total "305 r"impacts of the SAR model in Table~\ref{tab:spatial_results} "306 r"\citep{lesage2009introduction}: total $=\hat\beta/(1-\hat\rho)$, direct "307 r"$=\hat\beta\, n^{-1}\mathrm{tr}[(I-\hat\rho W)^{-1}]$, indirect $=$ total $-$ "308 r"direct. Standard errors in parentheses from 2,000 draws of $(\hat\beta,\hat\rho)$ "309 r"from their estimated asymptotic normal distribution. "310 r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"]311 (TABLE_DIR / "spatial_impacts.tex").write_text("\n".join(ilines) + "\n", encoding="utf-8")312 print(f" -> saved {TABLE_DIR / 'spatial_impacts.tex'}")313 # machine-readable record of the spatial estimates314 rec = {"rho": sar_params["W_log_rent"], "rho_se": sar_se["W_log_rent"],315 "lambda": float(sem.betas[-1][0]),316 "sar": {v: {"beta": sar_params[v], "se": sar_se[v]} for v in display_vars if v in sar_params},317 "sem": {v: {"beta": sem_params[v], "se": sem_se[v]} for v in display_vars if v in sem_params},318 "ols": {v: {"beta": ols_params[v], "se": ols_se[v]} for v in display_vars if v in ols_params},319 "impacts": {v: imp[v] for v in imp_vars},320 "mean_diag_exact": imp["mean_diag_exact"], "n": int(sar.n),321 "sar_pr2": float(sar.pr2), "sem_pr2": float(sem.pr2), "ols_r2": float(res_ols.rsquared)}322 (TABLE_DIR / "spatial_models.json").write_text(json.dumps(rec, indent=2), encoding="utf-8")323 lines.append(r"\midrule")324 lines.append(f"Observations & {int(res_ols.nobs)} & {sar.n} & {sem.n} " + r"\\")325 lines.append(f"R$^2$ / pseudo-R$^2$ & {res_ols.rsquared:.4f} & {sar.pr2:.4f} & {sem.pr2:.4f} " + r"\\")326 lines.append(r"\bottomrule")327 lines.append(r"\end{tabular}")328 lines.append(r"\parbox{\textwidth}{\footnotesize Standard errors in parentheses: HC1 (OLS), White-robust (SAR, spatial two-stage least squares, \texttt{GM\_Lag}), and the classical GM standard errors of \texttt{GM\_Error} (SEM; no standard error is reported for $\lambda$). $W$ is a row-standardised $k=5$ nearest-neighbour matrix on projected coordinates. All columns include building-type dummies and fixed effects for the 15 largest cities (remaining cities grouped), so the OLS column differs slightly from column (1c) of Table~\ref{tab:hedonic_rent}. $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}")329 return lines330331 # ── Fallback: manual spatial-lag OLS ─────────────────────────────────────332 print("\n Building KNN(k=5) manually with scipy ...")333 add_manual_spatial_lags(rent_sp, coords)334335 # Model with spatial lag of rent336 res_slag_rent = run_ols_robust(337 rent_sp, "log_rent",338 ["airbnb_count_500m", "spatial_lag_rent"] + controls + city_fe_cols,339 "OLS + spatial lag(rent)"340 )341342 # Model with both spatial lags343 res_slag_both = run_ols_robust(344 rent_sp, "log_rent",345 ["airbnb_count_500m", "spatial_lag_rent", "spatial_lag_airbnb"]346 + controls + city_fe_cols,347 "OLS + spatial lag(rent, airbnb)"348 )349350 show_vars_fb = [351 "const", "airbnb_count_500m", "spatial_lag_rent",352 "spatial_lag_airbnb", "bedrooms", "bathrooms",353 ] + bt_cols354355 lines = []356 lines.append(r"\begin{tabular}{lccc}")357 lines.append(r"\toprule")358 lines.append(r" & \textbf{OLS Baseline} & \textbf{+Lag(rent)} & \textbf{+Lag(rent, airbnb)} \\")359 lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log\_rent}} \\")360 lines.append(r"\midrule")361362 for var in show_vars_fb:363 cells_c, cells_s = [], []364 for res in [res_ols, res_slag_rent, res_slag_both]:365 if var in res.params.index:366 b = res.params[var]367 se = res.bse[var]368 p = res.pvalues[var]369 cells_c.append(f"{b:.4f}{_star(p)}")370 cells_s.append(f"({se:.4f})")371 else:372 cells_c.append("")373 cells_s.append("")374 vn = label(var)375 lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\")376 lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]")377378 lines.append(r"\midrule")379 for lbl, accessor in [380 ("Observations", lambda r: f"{int(r.nobs)}"),381 ("R$^2$", lambda r: f"{r.rsquared:.4f}"),382 ("Adj.\\ R$^2$", lambda r: f"{r.rsquared_adj:.4f}"),383 ]:384 row = [lbl]385 for res in [res_ols, res_slag_rent, res_slag_both]:386 row.append(accessor(res))387 lines.append(" & ".join(row) + r" \\")388389 lines.append(r"\bottomrule")390 lines.append(r"\end{tabular}")391 lines.append(392 r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in "393 r"parentheses. Spatial lag = mean of k=5 nearest neighbours. "394 r"City FE included but not shown. "395 r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"396 )397 return lines398399400def buffer_robustness(rent_m, controls, city_fe_cols) -> None:401 """Distance-buffer robustness: table + coefficient plot."""402 print("\n" + "-" * 72)403 print("Distance-buffer robustness (250m, 500m, 1km, 2km)")404 print("-" * 72)405406 buffers = {407 "250m": "airbnb_count_250m",408 "500m": "airbnb_count_500m",409 "1km": "airbnb_count_1000m",410 "2km": "airbnb_count_2000m",411 }412413 buffer_results = {}414 for buf_label, ab_var in buffers.items():415 if ab_var not in rent_m.columns:416 print(f" {ab_var} not found — skipping")417 continue418 x_cols = [ab_var] + controls + city_fe_cols419 res = run_ols_robust(rent_m, "log_rent", x_cols, f"buffer {buf_label}")420 buffer_results[buf_label] = (ab_var, res)421422 if not buffer_results:423 print(" No buffer variables found — skipping robustness table & plot.")424 return425426 # ── LaTeX table ──────────────────────────────────────────────────────────427 buf_labels = list(buffer_results.keys())428429 lines = []430 col_spec = "l" + "c" * len(buf_labels)431 lines.append(r"\begin{tabular}{" + col_spec + "}")432 lines.append(r"\toprule")433 header = " & ".join([""] + [f"\\textbf{{{b}}}" for b in buf_labels]) + r" \\"434 lines.append(header)435 lines.append(436 " & ".join(["Dep.\\ var:"] + [r"\textit{log rent}"] * len(buf_labels))437 + r" \\"438 )439 lines.append(r"\midrule")440441 # Show airbnb_count coefficient (the key variable differs per model)442 cells_c, cells_s = [], []443 for b in buf_labels:444 ab_var, res = buffer_results[b]445 bval = res.params[ab_var]446 se = res.bse[ab_var]447 p = res.pvalues[ab_var]448 cells_c.append(f"{bval:.6f}{_star(p)}")449 cells_s.append(f"({se:.6f})")450451 lines.append("Airbnb count & " + " & ".join(cells_c) + r" \\")452 lines.append(" & " + " & ".join(cells_s) + r" \\[4pt]")453454 # Controls row455 lines.append(456 "Controls & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\"457 )458 lines.append(459 "City FE & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\"460 )461462 lines.append(r"\midrule")463 # N, R2464 for lbl, acc in [465 ("Observations", lambda r: f"{int(r.nobs):,}"),466 ("R$^2$", lambda r: f"{r.rsquared:.4f}"),467 ]:468 cells = [lbl]469 for b in buf_labels:470 cells.append(acc(buffer_results[b][1]))471 lines.append(" & ".join(cells) + r" \\")472473 lines.append(r"\bottomrule")474 lines.append(r"\end{tabular}")475 lines.append(476 r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in "477 r"parentheses. Controls: bedrooms, bathrooms, building-type dummies; "478 r"fixed effects for the 15 largest cities (remaining cities grouped), as in "479 r"Table~\ref{tab:spatial_results}. "480 r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"481 )482483 tex_buf = "\n".join(lines) + "\n"484 out_buf = TABLE_DIR / "buffer_robustness.tex"485 out_buf.write_text(tex_buf, encoding="utf-8")486 print(f" -> saved {out_buf}")487488 # ── Coefficient plot ─────────────────────────────────────────────────────489 fig, ax = plt.subplots(figsize=(6, 4))490 x_pos = np.arange(len(buf_labels))491 coefs = []492 ci_lo = []493 ci_hi = []494 for b in buf_labels:495 ab_var, res = buffer_results[b]496 beta = res.params[ab_var]497 se = res.bse[ab_var]498 coefs.append(beta)499 ci_lo.append(beta - 1.96 * se)500 ci_hi.append(beta + 1.96 * se)501502 coefs = np.array(coefs)503 ci_lo = np.array(ci_lo)504 ci_hi = np.array(ci_hi)505 err_lo = coefs - ci_lo506 err_hi = ci_hi - coefs507508 ax.errorbar(509 x_pos, coefs,510 yerr=[err_lo, err_hi],511 fmt="o", capsize=5, capthick=1.5, color="steelblue", markersize=8,512 )513 ax.axhline(0, color="grey", linestyle="--", linewidth=0.7)514 ax.set_xticks(x_pos)515 ax.set_xticklabels(buf_labels)516 ax.set_xlabel("Buffer distance")517 ax.set_ylabel(r"$\beta$ (Airbnb count)")518 ax.set_title("Airbnb Count Coefficient by Buffer Distance (95% CI)")519 fig.tight_layout()520521 fig_path = FIG_DIR / "coefficient_buffer_comparison.pdf"522 fig.savefig(fig_path, dpi=300)523 plt.close(fig)524 print(f" -> saved {fig_path}")525526527def main() -> None:528 print("=" * 72)529 print("07 SPATIAL MODELS")530 print("=" * 72)531532 rent = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))533 print(f"\nRent data: {rent.shape[0]:,} rows")534535 # ── Prepare common variables ────────────────────────────────────────────536 rent_m = rent.copy()537538 # Building-type dummies539 if "building_type" in rent_m.columns:540 bt_dummies = pd.get_dummies(541 rent_m["building_type"], prefix="bt", drop_first=True, dtype=float542 )543 rent_m = pd.concat([rent_m, bt_dummies], axis=1)544 bt_cols = list(bt_dummies.columns)545 else:546 bt_cols = []547548 # City FE — limit to top 15 cities to avoid singular matrices in spatial models549 if "city" in rent_m.columns:550 top_cities = rent_m["city"].value_counts().nlargest(15).index551 city_col_reduced = rent_m["city"].where(rent_m["city"].isin(top_cities), other="Other")552 city_dummies = pd.get_dummies(553 city_col_reduced, prefix="city", drop_first=True, dtype=float554 )555 rent_m = pd.concat([rent_m, city_dummies], axis=1)556 city_fe_cols = list(city_dummies.columns)557 else:558 city_fe_cols = []559560 controls = ["bedrooms", "bathrooms"] + bt_cols561 all_x_baseline = ["airbnb_count_500m"] + controls + city_fe_cols562563 # Key display variables (short list for tables)564 display_vars = ["const", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols565566 # ═════════════════════════════════════════════════════════════════════════567 # MODEL 4: Spatial Analysis568 # ═════════════════════════════════════════════════════════════════════════569 print("\n" + "-" * 72)570 print("MODEL 4: Spatial Models")571 print("-" * 72)572573 # Determine coordinate columns574 lat_col = "lat" if "lat" in rent_m.columns else "latitude"575 lon_col = "lon" if "lon" in rent_m.columns else ("long" if "long" in rent_m.columns else "longitude")576577 # Subset to complete cases for spatial models578 spatial_cols = ["log_rent", lat_col, lon_col] + all_x_baseline579 rent_sp = rent_m.dropna(subset=spatial_cols).reset_index(drop=True)580 print(f" Complete cases for spatial analysis: {len(rent_sp):,}")581582 coords = project_km(rent_sp[lat_col].values.astype(float),583 rent_sp[lon_col].values.astype(float))584585 # OLS baseline on the spatial subsample586 res_ols = run_ols_robust(rent_sp, "log_rent", all_x_baseline,587 "OLS baseline (spatial sample)")588589 lines = spatial_models_table(rent_sp, coords, res_ols, all_x_baseline,590 controls, city_fe_cols, bt_cols, display_vars)591592 tex_spatial = "\n".join(lines) + "\n"593 out_spatial = TABLE_DIR / "spatial_models.tex"594 out_spatial.write_text(tex_spatial, encoding="utf-8")595 print(f" -> saved {out_spatial}")596597 # ═════════════════════════════════════════════════════════════════════════598 # DISTANCE-BUFFER ROBUSTNESS599 # ═════════════════════════════════════════════════════════════════════════600 buffer_robustness(rent_m, controls, city_fe_cols)601602 print("\n" + "=" * 72)603 print("07 DONE")604 print("=" * 72)605606607if __name__ == "__main__":608 main()609