spb/wp5_uqo Public
UQO Working Paper No. 5 — Airbnb, residential rents and housing market pressure.
TeX 53.4%
Python 46.5%
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 # noqa: E4023435PROCESSED_HINT = "Run scripts/04_merge_data.py first."3637# ── attempt to import spatial libraries ──────────────────────────────────────38try:39 from libpysal.weights import KNN as Wknn40 import spreg41 HAS_SPATIAL = True42 print(" libpysal + spreg available — will run SAR / SEM")43except ImportError:44 HAS_SPATIAL = False45 print(" libpysal/spreg not available — using fallback spatial approach")464748def run_ols_robust(df, y_col, x_cols, label):49 """OLS with HC1 robust SE."""50 sub = df[[y_col] + x_cols].dropna()51 Y = sub[y_col]52 X = sm.add_constant(sub[x_cols])53 res = sm.OLS(Y, X).fit(cov_type="HC1")54 print(f" [{label}] N={int(res.nobs):,} R2={res.rsquared:.4f}")55 return res565758def add_manual_spatial_lags(rent_sp: pd.DataFrame, coords: np.ndarray) -> None:59 """Add k=5 nearest-neighbour spatial lags of rent and Airbnb count in place."""60 tree = cKDTree(coords)61 _, idx = tree.query(coords, k=6) # k+1 because first neighbour is self62 idx_neighbours = idx[:, 1:] # drop self63 rent_sp["spatial_lag_rent"] = rent_sp["log_rent"].values[idx_neighbours].mean(axis=1)64 rent_sp["spatial_lag_airbnb"] = rent_sp["airbnb_count_500m"].values[idx_neighbours].mean(axis=1)656667def spatial_models_table(rent_sp, coords, res_ols, all_x_baseline,68 controls, city_fe_cols, bt_cols, display_vars):69 """70 Estimate the spatial models (SAR/SEM when spreg is available, manual71 spatial-lag OLS otherwise) and return the LaTeX table lines.72 """73 if HAS_SPATIAL:74 # ── Build KNN weights ────────────────────────────────────────────────75 print("\n Building KNN(k=5) spatial weights ...")76 w = Wknn.from_array(coords, k=5)77 w.transform = "r" # row-standardise7879 Y_arr = rent_sp["log_rent"].values.reshape(-1, 1)80 X_arr = sm.add_constant(rent_sp[all_x_baseline].values)81 var_names = ["const"] + all_x_baseline8283 # ── SAR (Spatial Lag Model) ──────────────────────────────────────────84 sar = None85 sem = None86 try:87 print(" Estimating Spatial Lag Model (SAR) ...")88 sar = spreg.GM_Lag(89 Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names90 )91 print(f" [SAR] N={sar.n} pseudo-R2={sar.pr2:.4f} rho={sar.betas[-1][0]:.4f}")92 except Exception as e:93 print(f" [SAR] Failed: {e}")9495 # ── SEM (Spatial Error Model) ────────────────────────────────────────96 try:97 print(" Estimating Spatial Error Model (SEM) ...")98 sem = spreg.GM_Error(99 Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names100 )101 sem_lambda = sem.betas[-1][0]102 print(f" [SEM] N={sem.n} pseudo-R2={sem.pr2:.4f} lambda={sem_lambda:.4f}")103 except Exception as e:104 print(f" [SEM] Failed: {e}")105106 # ── If SAR/SEM failed, fall back to manual spatial-lag approach ─────107 if sar is None or sem is None:108 print(" SAR/SEM failed — falling back to manual spatial lag approach ...")109 add_manual_spatial_lags(rent_sp, coords)110 fb_x = all_x_baseline + ["spatial_lag_rent", "spatial_lag_airbnb"]111 res_fb = run_ols_robust(rent_sp, "log_rent", fb_x, "OLS + spatial lags")112 lines: list[str] = []113 lines.append(r"\begin{tabular}{lcc}")114 lines.append(r"\toprule")115 lines.append(r" & \textbf{OLS Baseline} & \textbf{OLS + Spatial Lags} \\")116 lines.append(r"\midrule")117 ols_p = dict(zip(res_ols.model.exog_names, res_ols.params))118 ols_s = dict(zip(res_ols.model.exog_names, res_ols.bse))119 ols_pv = dict(zip(res_ols.model.exog_names, res_ols.pvalues))120 fb_p = dict(zip(res_fb.model.exog_names, res_fb.params))121 fb_s = dict(zip(res_fb.model.exog_names, res_fb.bse))122 fb_pv = dict(zip(res_fb.model.exog_names, res_fb.pvalues))123 show = (["const", "airbnb_count_500m", "bedrooms", "bathrooms"]124 + bt_cols + ["spatial_lag_rent", "spatial_lag_airbnb"])125 for var in show:126 c1 = f"{ols_p[var]:.4f}{_star(ols_pv[var])}" if var in ols_p else ""127 s1 = f"({ols_s[var]:.4f})" if var in ols_s else ""128 c2 = f"{fb_p[var]:.4f}{_star(fb_pv[var])}" if var in fb_p else ""129 s2 = f"({fb_s[var]:.4f})" if var in fb_s else ""130 vn = var.replace("_", r"\_")131 lines.append(f"{vn} & {c1} & {c2}" + r" \\")132 lines.append(f" & {s1} & {s2}" + r" \\[4pt]")133 lines.append(r"\midrule")134 lines.append(f"Observations & {int(res_ols.nobs)} & {int(res_fb.nobs)}" + r" \\")135 lines.append(f"R$^2$ & {res_ols.rsquared:.4f} & {res_fb.rsquared:.4f}" + r" \\")136 lines.append(r"\bottomrule")137 lines.append(r"\end{tabular}")138 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.}")139 return lines140141 # ── Build LaTeX table: OLS vs SAR vs SEM ─────────────────────────────142 lines = []143 lines.append(r"\begin{tabular}{lccc}")144 lines.append(r"\toprule")145 lines.append(r" & \textbf{OLS} & \textbf{SAR (GM\_Lag)} & \textbf{SEM (GM\_Error)} \\")146 lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log\_rent}} \\")147 lines.append(r"\midrule")148 ols_params = dict(zip(["const"] + all_x_baseline, res_ols.params))149 ols_se = dict(zip(["const"] + all_x_baseline, res_ols.bse))150 ols_pval = dict(zip(["const"] + all_x_baseline, res_ols.pvalues))151 sar_var_names = var_names + ["W_log_rent"]152 sar_params = {v: sar.betas[i][0] for i, v in enumerate(sar_var_names)}153 sar_se = {v: sar.std_err[i] for i, v in enumerate(sar_var_names)}154 sar_z = {v: sar.z_stat[i] for i, v in enumerate(sar_var_names)}155 sem_params = {v: sem.betas[i][0] for i, v in enumerate(var_names)}156 sem_se = {v: sem.std_err[i] for i, v in enumerate(var_names)}157 sem_z = {v: sem.z_stat[i] for i, v in enumerate(var_names)}158 show_vars = display_vars + ["W_log_rent"]159 for var in show_vars:160 cells_c, cells_s = [], []161 if var in ols_params:162 cells_c.append(f"{ols_params[var]:.4f}{_star(ols_pval[var])}")163 cells_s.append(f"({ols_se[var]:.4f})")164 else:165 cells_c.append("")166 cells_s.append("")167 if var in sar_params:168 p_sar = 2 * (1 - scipy.stats.norm.cdf(abs(sar_z[var][0])))169 cells_c.append(f"{sar_params[var]:.4f}{_star(p_sar)}")170 cells_s.append(f"({sar_se[var]:.4f})")171 else:172 cells_c.append("")173 cells_s.append("")174 if var in sem_params:175 p_sem = 2 * (1 - scipy.stats.norm.cdf(abs(sem_z[var][0])))176 cells_c.append(f"{sem_params[var]:.4f}{_star(p_sem)}")177 cells_s.append(f"({sem_se[var]:.4f})")178 else:179 cells_c.append("")180 cells_s.append("")181 vn = var.replace("_", r"\_")182 lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\")183 lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]")184 lines.append(r"$\lambda$ (spatial error) & & & " + f"{sem.betas[-1][0]:.4f}" + r" \\")185 lines.append(r"\midrule")186 lines.append(f"Observations & {int(res_ols.nobs)} & {sar.n} & {sem.n} " + r"\\")187 lines.append(f"R$^2$ / pseudo-R$^2$ & {res_ols.rsquared:.4f} & {sar.pr2:.4f} & {sem.pr2:.4f} " + r"\\")188 lines.append(r"\bottomrule")189 lines.append(r"\end{tabular}")190 lines.append(r"\parbox{\textwidth}{\footnotesize Standard errors in parentheses. SAR estimated via GM\_Lag; SEM via GM\_Error. $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$. City FE included but not shown.}")191 return lines192193 # ── Fallback: manual spatial-lag OLS ─────────────────────────────────────194 print("\n Building KNN(k=5) manually with scipy ...")195 add_manual_spatial_lags(rent_sp, coords)196197 # Model with spatial lag of rent198 res_slag_rent = run_ols_robust(199 rent_sp, "log_rent",200 ["airbnb_count_500m", "spatial_lag_rent"] + controls + city_fe_cols,201 "OLS + spatial lag(rent)"202 )203204 # Model with both spatial lags205 res_slag_both = run_ols_robust(206 rent_sp, "log_rent",207 ["airbnb_count_500m", "spatial_lag_rent", "spatial_lag_airbnb"]208 + controls + city_fe_cols,209 "OLS + spatial lag(rent, airbnb)"210 )211212 show_vars_fb = [213 "const", "airbnb_count_500m", "spatial_lag_rent",214 "spatial_lag_airbnb", "bedrooms", "bathrooms",215 ] + bt_cols216217 lines = []218 lines.append(r"\begin{tabular}{lccc}")219 lines.append(r"\toprule")220 lines.append(r" & \textbf{OLS Baseline} & \textbf{+Lag(rent)} & \textbf{+Lag(rent, airbnb)} \\")221 lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log\_rent}} \\")222 lines.append(r"\midrule")223224 for var in show_vars_fb:225 cells_c, cells_s = [], []226 for res in [res_ols, res_slag_rent, res_slag_both]:227 if var in res.params.index:228 b = res.params[var]229 se = res.bse[var]230 p = res.pvalues[var]231 cells_c.append(f"{b:.4f}{_star(p)}")232 cells_s.append(f"({se:.4f})")233 else:234 cells_c.append("")235 cells_s.append("")236 vn = var.replace("_", r"\_")237 lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\")238 lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]")239240 lines.append(r"\midrule")241 for label, accessor in [242 ("Observations", lambda r: f"{int(r.nobs)}"),243 ("R$^2$", lambda r: f"{r.rsquared:.4f}"),244 ("Adj.\\ R$^2$", lambda r: f"{r.rsquared_adj:.4f}"),245 ]:246 row = [label]247 for res in [res_ols, res_slag_rent, res_slag_both]:248 row.append(accessor(res))249 lines.append(" & ".join(row) + r" \\")250251 lines.append(r"\bottomrule")252 lines.append(r"\end{tabular}")253 lines.append(254 r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in "255 r"parentheses. Spatial lag = mean of k=5 nearest neighbours. "256 r"City FE included but not shown. "257 r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"258 )259 return lines260261262def buffer_robustness(rent_m, controls, city_fe_cols) -> None:263 """Distance-buffer robustness: table + coefficient plot."""264 print("\n" + "-" * 72)265 print("Distance-buffer robustness (250m, 500m, 1km, 2km)")266 print("-" * 72)267268 buffers = {269 "250m": "airbnb_count_250m",270 "500m": "airbnb_count_500m",271 "1km": "airbnb_count_1000m",272 "2km": "airbnb_count_2000m",273 }274275 buffer_results = {}276 for buf_label, ab_var in buffers.items():277 if ab_var not in rent_m.columns:278 print(f" {ab_var} not found — skipping")279 continue280 x_cols = [ab_var] + controls + city_fe_cols281 res = run_ols_robust(rent_m, "log_rent", x_cols, f"buffer {buf_label}")282 buffer_results[buf_label] = (ab_var, res)283284 if not buffer_results:285 print(" No buffer variables found — skipping robustness table & plot.")286 return287288 # ── LaTeX table ──────────────────────────────────────────────────────────289 buf_labels = list(buffer_results.keys())290291 lines = []292 col_spec = "l" + "c" * len(buf_labels)293 lines.append(r"\begin{tabular}{" + col_spec + "}")294 lines.append(r"\toprule")295 header = " & ".join([""] + [f"\\textbf{{{b}}}" for b in buf_labels]) + r" \\"296 lines.append(header)297 lines.append(298 " & ".join(["Dep.\\ var:"] + [r"\textit{log\_rent}"] * len(buf_labels))299 + r" \\"300 )301 lines.append(r"\midrule")302303 # Show airbnb_count coefficient (the key variable differs per model)304 cells_c, cells_s = [], []305 for b in buf_labels:306 ab_var, res = buffer_results[b]307 bval = res.params[ab_var]308 se = res.bse[ab_var]309 p = res.pvalues[ab_var]310 cells_c.append(f"{bval:.6f}{_star(p)}")311 cells_s.append(f"({se:.6f})")312313 lines.append("Airbnb count & " + " & ".join(cells_c) + r" \\")314 lines.append(" & " + " & ".join(cells_s) + r" \\[4pt]")315316 # Controls row317 lines.append(318 "Controls & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\"319 )320 lines.append(321 "City FE & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\"322 )323324 lines.append(r"\midrule")325 # N, R2326 for lbl, acc in [327 ("Observations", lambda r: f"{int(r.nobs):,}"),328 ("R$^2$", lambda r: f"{r.rsquared:.4f}"),329 ]:330 cells = [lbl]331 for b in buf_labels:332 cells.append(acc(buffer_results[b][1]))333 lines.append(" & ".join(cells) + r" \\")334335 lines.append(r"\bottomrule")336 lines.append(r"\end{tabular}")337 lines.append(338 r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in "339 r"parentheses. Controls: bedrooms, bathrooms, building-type dummies. "340 r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"341 )342343 tex_buf = "\n".join(lines) + "\n"344 out_buf = TABLE_DIR / "buffer_robustness.tex"345 out_buf.write_text(tex_buf, encoding="utf-8")346 print(f" -> saved {out_buf}")347348 # ── Coefficient plot ─────────────────────────────────────────────────────349 fig, ax = plt.subplots(figsize=(6, 4))350 x_pos = np.arange(len(buf_labels))351 coefs = []352 ci_lo = []353 ci_hi = []354 for b in buf_labels:355 ab_var, res = buffer_results[b]356 beta = res.params[ab_var]357 se = res.bse[ab_var]358 coefs.append(beta)359 ci_lo.append(beta - 1.96 * se)360 ci_hi.append(beta + 1.96 * se)361362 coefs = np.array(coefs)363 ci_lo = np.array(ci_lo)364 ci_hi = np.array(ci_hi)365 err_lo = coefs - ci_lo366 err_hi = ci_hi - coefs367368 ax.errorbar(369 x_pos, coefs,370 yerr=[err_lo, err_hi],371 fmt="o", capsize=5, capthick=1.5, color="steelblue", markersize=8,372 )373 ax.axhline(0, color="grey", linestyle="--", linewidth=0.7)374 ax.set_xticks(x_pos)375 ax.set_xticklabels(buf_labels)376 ax.set_xlabel("Buffer distance")377 ax.set_ylabel(r"$\beta$ (Airbnb count)")378 ax.set_title("Airbnb Count Coefficient by Buffer Distance (95% CI)")379 fig.tight_layout()380381 fig_path = FIG_DIR / "coefficient_buffer_comparison.pdf"382 fig.savefig(fig_path, dpi=300)383 plt.close(fig)384 print(f" -> saved {fig_path}")385386387def main() -> None:388 print("=" * 72)389 print("07 SPATIAL MODELS")390 print("=" * 72)391392 rent = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))393 print(f"\nRent data: {rent.shape[0]:,} rows")394395 # ── Prepare common variables ────────────────────────────────────────────396 rent_m = rent.copy()397398 # Building-type dummies399 if "building_type" in rent_m.columns:400 bt_dummies = pd.get_dummies(401 rent_m["building_type"], prefix="bt", drop_first=True, dtype=float402 )403 rent_m = pd.concat([rent_m, bt_dummies], axis=1)404 bt_cols = list(bt_dummies.columns)405 else:406 bt_cols = []407408 # City FE — limit to top 15 cities to avoid singular matrices in spatial models409 if "city" in rent_m.columns:410 top_cities = rent_m["city"].value_counts().nlargest(15).index411 city_col_reduced = rent_m["city"].where(rent_m["city"].isin(top_cities), other="Other")412 city_dummies = pd.get_dummies(413 city_col_reduced, prefix="city", drop_first=True, dtype=float414 )415 rent_m = pd.concat([rent_m, city_dummies], axis=1)416 city_fe_cols = list(city_dummies.columns)417 else:418 city_fe_cols = []419420 controls = ["bedrooms", "bathrooms"] + bt_cols421 all_x_baseline = ["airbnb_count_500m"] + controls + city_fe_cols422423 # Key display variables (short list for tables)424 display_vars = ["const", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols425426 # ═════════════════════════════════════════════════════════════════════════427 # MODEL 4: Spatial Analysis428 # ═════════════════════════════════════════════════════════════════════════429 print("\n" + "-" * 72)430 print("MODEL 4: Spatial Models")431 print("-" * 72)432433 # Determine coordinate columns434 lat_col = "lat" if "lat" in rent_m.columns else "latitude"435 lon_col = "lon" if "lon" in rent_m.columns else ("long" if "long" in rent_m.columns else "longitude")436437 # Subset to complete cases for spatial models438 spatial_cols = ["log_rent", lat_col, lon_col] + all_x_baseline439 rent_sp = rent_m.dropna(subset=spatial_cols).reset_index(drop=True)440 print(f" Complete cases for spatial analysis: {len(rent_sp):,}")441442 coords = rent_sp[[lat_col, lon_col]].values443444 # OLS baseline on the spatial subsample445 res_ols = run_ols_robust(rent_sp, "log_rent", all_x_baseline,446 "OLS baseline (spatial sample)")447448 lines = spatial_models_table(rent_sp, coords, res_ols, all_x_baseline,449 controls, city_fe_cols, bt_cols, display_vars)450451 tex_spatial = "\n".join(lines) + "\n"452 out_spatial = TABLE_DIR / "spatial_models.tex"453 out_spatial.write_text(tex_spatial, encoding="utf-8")454 print(f" -> saved {out_spatial}")455456 # ═════════════════════════════════════════════════════════════════════════457 # DISTANCE-BUFFER ROBUSTNESS458 # ═════════════════════════════════════════════════════════════════════════459 buffer_robustness(rent_m, controls, city_fe_cols)460461 print("\n" + "=" * 72)462 print("07 DONE")463 print("=" * 72)464465466if __name__ == "__main__":467 main()468