# Author: Simon-Pierre Boucher — contact@spboucher.ai """ 07_spatial_models.py -------------------- Spatial regression analysis for the Airbnb-rent study. Model 4: Spatial models - If libpysal + spreg available: SAR, SEM, and OLS comparison - Fallback: manual spatial-lag OLS ("poor man's spatial model") - Distance-buffer robustness (250m, 500m, 1km, 2km) Outputs: results/tables/spatial_models.tex results/tables/buffer_robustness.tex figures/coefficient_buffer_comparison.pdf """ import sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 import numpy as np # noqa: E402 import pandas as pd # noqa: E402 import scipy.stats # noqa: E402 import statsmodels.api as sm # noqa: E402 from scipy.spatial import cKDTree # noqa: E402 from src.config import MERGED_ANALYSIS, TABLE_DIR, FIG_DIR, require # noqa: E402 from src.latex_tables import significance_star as _star # noqa: E402 PROCESSED_HINT = "Run scripts/04_merge_data.py first." # ── attempt to import spatial libraries ────────────────────────────────────── try: from libpysal.weights import KNN as Wknn import spreg HAS_SPATIAL = True print(" libpysal + spreg available — will run SAR / SEM") except ImportError: HAS_SPATIAL = False print(" libpysal/spreg not available — using fallback spatial approach") def run_ols_robust(df, y_col, x_cols, label): """OLS with HC1 robust SE.""" sub = df[[y_col] + x_cols].dropna() Y = sub[y_col] X = sm.add_constant(sub[x_cols]) res = sm.OLS(Y, X).fit(cov_type="HC1") print(f" [{label}] N={int(res.nobs):,} R2={res.rsquared:.4f}") return res def add_manual_spatial_lags(rent_sp: pd.DataFrame, coords: np.ndarray) -> None: """Add k=5 nearest-neighbour spatial lags of rent and Airbnb count in place.""" tree = cKDTree(coords) _, idx = tree.query(coords, k=6) # k+1 because first neighbour is self idx_neighbours = idx[:, 1:] # drop self rent_sp["spatial_lag_rent"] = rent_sp["log_rent"].values[idx_neighbours].mean(axis=1) rent_sp["spatial_lag_airbnb"] = rent_sp["airbnb_count_500m"].values[idx_neighbours].mean(axis=1) def spatial_models_table(rent_sp, coords, res_ols, all_x_baseline, controls, city_fe_cols, bt_cols, display_vars): """ Estimate the spatial models (SAR/SEM when spreg is available, manual spatial-lag OLS otherwise) and return the LaTeX table lines. """ if HAS_SPATIAL: # ── Build KNN weights ──────────────────────────────────────────────── print("\n Building KNN(k=5) spatial weights ...") w = Wknn.from_array(coords, k=5) w.transform = "r" # row-standardise Y_arr = rent_sp["log_rent"].values.reshape(-1, 1) X_arr = sm.add_constant(rent_sp[all_x_baseline].values) var_names = ["const"] + all_x_baseline # ── SAR (Spatial Lag Model) ────────────────────────────────────────── sar = None sem = None try: print(" Estimating Spatial Lag Model (SAR) ...") sar = spreg.GM_Lag( Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names ) print(f" [SAR] N={sar.n} pseudo-R2={sar.pr2:.4f} rho={sar.betas[-1][0]:.4f}") except Exception as e: print(f" [SAR] Failed: {e}") # ── SEM (Spatial Error Model) ──────────────────────────────────────── try: print(" Estimating Spatial Error Model (SEM) ...") sem = spreg.GM_Error( Y_arr, X_arr, w=w, name_y="log_rent", name_x=var_names ) sem_lambda = sem.betas[-1][0] print(f" [SEM] N={sem.n} pseudo-R2={sem.pr2:.4f} lambda={sem_lambda:.4f}") except Exception as e: print(f" [SEM] Failed: {e}") # ── If SAR/SEM failed, fall back to manual spatial-lag approach ───── if sar is None or sem is None: print(" SAR/SEM failed — falling back to manual spatial lag approach ...") add_manual_spatial_lags(rent_sp, coords) fb_x = all_x_baseline + ["spatial_lag_rent", "spatial_lag_airbnb"] res_fb = run_ols_robust(rent_sp, "log_rent", fb_x, "OLS + spatial lags") lines: list[str] = [] lines.append(r"\begin{tabular}{lcc}") lines.append(r"\toprule") lines.append(r" & \textbf{OLS Baseline} & \textbf{OLS + Spatial Lags} \\") lines.append(r"\midrule") ols_p = dict(zip(res_ols.model.exog_names, res_ols.params)) ols_s = dict(zip(res_ols.model.exog_names, res_ols.bse)) ols_pv = dict(zip(res_ols.model.exog_names, res_ols.pvalues)) fb_p = dict(zip(res_fb.model.exog_names, res_fb.params)) fb_s = dict(zip(res_fb.model.exog_names, res_fb.bse)) fb_pv = dict(zip(res_fb.model.exog_names, res_fb.pvalues)) show = (["const", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols + ["spatial_lag_rent", "spatial_lag_airbnb"]) for var in show: c1 = f"{ols_p[var]:.4f}{_star(ols_pv[var])}" if var in ols_p else "" s1 = f"({ols_s[var]:.4f})" if var in ols_s else "" c2 = f"{fb_p[var]:.4f}{_star(fb_pv[var])}" if var in fb_p else "" s2 = f"({fb_s[var]:.4f})" if var in fb_s else "" vn = var.replace("_", r"\_") lines.append(f"{vn} & {c1} & {c2}" + r" \\") lines.append(f" & {s1} & {s2}" + r" \\[4pt]") lines.append(r"\midrule") lines.append(f"Observations & {int(res_ols.nobs)} & {int(res_fb.nobs)}" + r" \\") lines.append(f"R$^2$ & {res_ols.rsquared:.4f} & {res_fb.rsquared:.4f}" + r" \\") lines.append(r"\bottomrule") lines.append(r"\end{tabular}") 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.}") return lines # ── Build LaTeX table: OLS vs SAR vs SEM ───────────────────────────── lines = [] lines.append(r"\begin{tabular}{lccc}") lines.append(r"\toprule") lines.append(r" & \textbf{OLS} & \textbf{SAR (GM\_Lag)} & \textbf{SEM (GM\_Error)} \\") lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log\_rent}} \\") lines.append(r"\midrule") ols_params = dict(zip(["const"] + all_x_baseline, res_ols.params)) ols_se = dict(zip(["const"] + all_x_baseline, res_ols.bse)) ols_pval = dict(zip(["const"] + all_x_baseline, res_ols.pvalues)) sar_var_names = var_names + ["W_log_rent"] sar_params = {v: sar.betas[i][0] for i, v in enumerate(sar_var_names)} sar_se = {v: sar.std_err[i] for i, v in enumerate(sar_var_names)} sar_z = {v: sar.z_stat[i] for i, v in enumerate(sar_var_names)} sem_params = {v: sem.betas[i][0] for i, v in enumerate(var_names)} sem_se = {v: sem.std_err[i] for i, v in enumerate(var_names)} sem_z = {v: sem.z_stat[i] for i, v in enumerate(var_names)} show_vars = display_vars + ["W_log_rent"] for var in show_vars: cells_c, cells_s = [], [] if var in ols_params: cells_c.append(f"{ols_params[var]:.4f}{_star(ols_pval[var])}") cells_s.append(f"({ols_se[var]:.4f})") else: cells_c.append("") cells_s.append("") if var in sar_params: p_sar = 2 * (1 - scipy.stats.norm.cdf(abs(sar_z[var][0]))) cells_c.append(f"{sar_params[var]:.4f}{_star(p_sar)}") cells_s.append(f"({sar_se[var]:.4f})") else: cells_c.append("") cells_s.append("") if var in sem_params: p_sem = 2 * (1 - scipy.stats.norm.cdf(abs(sem_z[var][0]))) cells_c.append(f"{sem_params[var]:.4f}{_star(p_sem)}") cells_s.append(f"({sem_se[var]:.4f})") else: cells_c.append("") cells_s.append("") vn = var.replace("_", r"\_") lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\") lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]") lines.append(r"$\lambda$ (spatial error) & & & " + f"{sem.betas[-1][0]:.4f}" + r" \\") lines.append(r"\midrule") lines.append(f"Observations & {int(res_ols.nobs)} & {sar.n} & {sem.n} " + r"\\") lines.append(f"R$^2$ / pseudo-R$^2$ & {res_ols.rsquared:.4f} & {sar.pr2:.4f} & {sem.pr2:.4f} " + r"\\") lines.append(r"\bottomrule") lines.append(r"\end{tabular}") 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.}") return lines # ── Fallback: manual spatial-lag OLS ───────────────────────────────────── print("\n Building KNN(k=5) manually with scipy ...") add_manual_spatial_lags(rent_sp, coords) # Model with spatial lag of rent res_slag_rent = run_ols_robust( rent_sp, "log_rent", ["airbnb_count_500m", "spatial_lag_rent"] + controls + city_fe_cols, "OLS + spatial lag(rent)" ) # Model with both spatial lags res_slag_both = run_ols_robust( rent_sp, "log_rent", ["airbnb_count_500m", "spatial_lag_rent", "spatial_lag_airbnb"] + controls + city_fe_cols, "OLS + spatial lag(rent, airbnb)" ) show_vars_fb = [ "const", "airbnb_count_500m", "spatial_lag_rent", "spatial_lag_airbnb", "bedrooms", "bathrooms", ] + bt_cols lines = [] lines.append(r"\begin{tabular}{lccc}") lines.append(r"\toprule") lines.append(r" & \textbf{OLS Baseline} & \textbf{+Lag(rent)} & \textbf{+Lag(rent, airbnb)} \\") lines.append(r"Dep.\ var: & \multicolumn{3}{c}{\textit{log\_rent}} \\") lines.append(r"\midrule") for var in show_vars_fb: cells_c, cells_s = [], [] for res in [res_ols, res_slag_rent, res_slag_both]: if var in res.params.index: b = res.params[var] se = res.bse[var] p = res.pvalues[var] cells_c.append(f"{b:.4f}{_star(p)}") cells_s.append(f"({se:.4f})") else: cells_c.append("") cells_s.append("") vn = var.replace("_", r"\_") lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\") lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]") lines.append(r"\midrule") for label, accessor in [ ("Observations", lambda r: f"{int(r.nobs)}"), ("R$^2$", lambda r: f"{r.rsquared:.4f}"), ("Adj.\\ R$^2$", lambda r: f"{r.rsquared_adj:.4f}"), ]: row = [label] for res in [res_ols, res_slag_rent, res_slag_both]: row.append(accessor(res)) lines.append(" & ".join(row) + r" \\") lines.append(r"\bottomrule") lines.append(r"\end{tabular}") lines.append( r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in " r"parentheses. Spatial lag = mean of k=5 nearest neighbours. " r"City FE included but not shown. " r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}" ) return lines def buffer_robustness(rent_m, controls, city_fe_cols) -> None: """Distance-buffer robustness: table + coefficient plot.""" print("\n" + "-" * 72) print("Distance-buffer robustness (250m, 500m, 1km, 2km)") print("-" * 72) buffers = { "250m": "airbnb_count_250m", "500m": "airbnb_count_500m", "1km": "airbnb_count_1000m", "2km": "airbnb_count_2000m", } buffer_results = {} for buf_label, ab_var in buffers.items(): if ab_var not in rent_m.columns: print(f" {ab_var} not found — skipping") continue x_cols = [ab_var] + controls + city_fe_cols res = run_ols_robust(rent_m, "log_rent", x_cols, f"buffer {buf_label}") buffer_results[buf_label] = (ab_var, res) if not buffer_results: print(" No buffer variables found — skipping robustness table & plot.") return # ── LaTeX table ────────────────────────────────────────────────────────── buf_labels = list(buffer_results.keys()) lines = [] col_spec = "l" + "c" * len(buf_labels) lines.append(r"\begin{tabular}{" + col_spec + "}") lines.append(r"\toprule") header = " & ".join([""] + [f"\\textbf{{{b}}}" for b in buf_labels]) + r" \\" lines.append(header) lines.append( " & ".join(["Dep.\\ var:"] + [r"\textit{log\_rent}"] * len(buf_labels)) + r" \\" ) lines.append(r"\midrule") # Show airbnb_count coefficient (the key variable differs per model) cells_c, cells_s = [], [] for b in buf_labels: ab_var, res = buffer_results[b] bval = res.params[ab_var] se = res.bse[ab_var] p = res.pvalues[ab_var] cells_c.append(f"{bval:.6f}{_star(p)}") cells_s.append(f"({se:.6f})") lines.append("Airbnb count & " + " & ".join(cells_c) + r" \\") lines.append(" & " + " & ".join(cells_s) + r" \\[4pt]") # Controls row lines.append( "Controls & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\" ) lines.append( "City FE & " + " & ".join(["Yes"] * len(buf_labels)) + r" \\" ) lines.append(r"\midrule") # N, R2 for lbl, acc in [ ("Observations", lambda r: f"{int(r.nobs):,}"), ("R$^2$", lambda r: f"{r.rsquared:.4f}"), ]: cells = [lbl] for b in buf_labels: cells.append(acc(buffer_results[b][1])) lines.append(" & ".join(cells) + r" \\") lines.append(r"\bottomrule") lines.append(r"\end{tabular}") lines.append( r"\parbox{\textwidth}{\footnotesize Robust (HC1) standard errors in " r"parentheses. Controls: bedrooms, bathrooms, building-type dummies. " r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}" ) tex_buf = "\n".join(lines) + "\n" out_buf = TABLE_DIR / "buffer_robustness.tex" out_buf.write_text(tex_buf, encoding="utf-8") print(f" -> saved {out_buf}") # ── Coefficient plot ───────────────────────────────────────────────────── fig, ax = plt.subplots(figsize=(6, 4)) x_pos = np.arange(len(buf_labels)) coefs = [] ci_lo = [] ci_hi = [] for b in buf_labels: ab_var, res = buffer_results[b] beta = res.params[ab_var] se = res.bse[ab_var] coefs.append(beta) ci_lo.append(beta - 1.96 * se) ci_hi.append(beta + 1.96 * se) coefs = np.array(coefs) ci_lo = np.array(ci_lo) ci_hi = np.array(ci_hi) err_lo = coefs - ci_lo err_hi = ci_hi - coefs ax.errorbar( x_pos, coefs, yerr=[err_lo, err_hi], fmt="o", capsize=5, capthick=1.5, color="steelblue", markersize=8, ) ax.axhline(0, color="grey", linestyle="--", linewidth=0.7) ax.set_xticks(x_pos) ax.set_xticklabels(buf_labels) ax.set_xlabel("Buffer distance") ax.set_ylabel(r"$\beta$ (Airbnb count)") ax.set_title("Airbnb Count Coefficient by Buffer Distance (95% CI)") fig.tight_layout() fig_path = FIG_DIR / "coefficient_buffer_comparison.pdf" fig.savefig(fig_path, dpi=300) plt.close(fig) print(f" -> saved {fig_path}") def main() -> None: print("=" * 72) print("07 SPATIAL MODELS") print("=" * 72) rent = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT)) print(f"\nRent data: {rent.shape[0]:,} rows") # ── Prepare common variables ──────────────────────────────────────────── rent_m = rent.copy() # Building-type dummies if "building_type" in rent_m.columns: bt_dummies = pd.get_dummies( rent_m["building_type"], prefix="bt", drop_first=True, dtype=float ) rent_m = pd.concat([rent_m, bt_dummies], axis=1) bt_cols = list(bt_dummies.columns) else: bt_cols = [] # City FE — limit to top 15 cities to avoid singular matrices in spatial models if "city" in rent_m.columns: top_cities = rent_m["city"].value_counts().nlargest(15).index city_col_reduced = rent_m["city"].where(rent_m["city"].isin(top_cities), other="Other") city_dummies = pd.get_dummies( city_col_reduced, prefix="city", drop_first=True, dtype=float ) rent_m = pd.concat([rent_m, city_dummies], axis=1) city_fe_cols = list(city_dummies.columns) else: city_fe_cols = [] controls = ["bedrooms", "bathrooms"] + bt_cols all_x_baseline = ["airbnb_count_500m"] + controls + city_fe_cols # Key display variables (short list for tables) display_vars = ["const", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols # ═════════════════════════════════════════════════════════════════════════ # MODEL 4: Spatial Analysis # ═════════════════════════════════════════════════════════════════════════ print("\n" + "-" * 72) print("MODEL 4: Spatial Models") print("-" * 72) # Determine coordinate columns lat_col = "lat" if "lat" in rent_m.columns else "latitude" lon_col = "lon" if "lon" in rent_m.columns else ("long" if "long" in rent_m.columns else "longitude") # Subset to complete cases for spatial models spatial_cols = ["log_rent", lat_col, lon_col] + all_x_baseline rent_sp = rent_m.dropna(subset=spatial_cols).reset_index(drop=True) print(f" Complete cases for spatial analysis: {len(rent_sp):,}") coords = rent_sp[[lat_col, lon_col]].values # OLS baseline on the spatial subsample res_ols = run_ols_robust(rent_sp, "log_rent", all_x_baseline, "OLS baseline (spatial sample)") lines = spatial_models_table(rent_sp, coords, res_ols, all_x_baseline, controls, city_fe_cols, bt_cols, display_vars) tex_spatial = "\n".join(lines) + "\n" out_spatial = TABLE_DIR / "spatial_models.tex" out_spatial.write_text(tex_spatial, encoding="utf-8") print(f" -> saved {out_spatial}") # ═════════════════════════════════════════════════════════════════════════ # DISTANCE-BUFFER ROBUSTNESS # ═════════════════════════════════════════════════════════════════════════ buffer_robustness(rent_m, controls, city_fe_cols) print("\n" + "=" * 72) print("07 DONE") print("=" * 72) if __name__ == "__main__": main()