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"""310_robustness_tables_figures.py4-------------------------------5Generate additional publication-quality robustness tables and figures.67Robustness tables8-----------------91. Winsorization / outlier robustness -> results/tables/robustness_outliers.tex102. Subsample analysis -> results/tables/robustness_subsamples.tex1112Additional figures13------------------141. figures/rent_airbnb_heatmap.pdf — hexbin of Montreal listings152. figures/coefficient_robustness.pdf — coefficient comparison plot163. figures/rent_by_airbnb_bins.pdf — mean rent by airbnb quintiles17"""1819import sys20import warnings21from pathlib import Path2223sys.path.insert(0, str(Path(__file__).resolve().parent.parent))2425from src.plotting import use_publication_style2627use_publication_style()2829import matplotlib.pyplot as plt # noqa: E40230import numpy as np # noqa: E40231import pandas as pd # noqa: E40232import statsmodels.api as sm # noqa: E4023334from src.config import MERGED_ANALYSIS, TABLE_DIR, FIG_DIR, require # noqa: E40235from src.latex_tables import significance_star as _stars_raw # noqa: E4023637warnings.filterwarnings("ignore")3839PROCESSED_HINT = "Run scripts/04_merge_data.py first."404142def _stars(pval: float) -> str:43 """Return significance stars (empty for NaN p-values)."""44 if pd.isna(pval):45 return ""46 return _stars_raw(pval)474849# ─────────────────────────────────────────────────────────────────────────────50# Helpers: OLS with city fixed effects51# ─────────────────────────────────────────────────────────────────────────────5253def run_ols_city_fe(data: pd.DataFrame, label: str,54 airbnb_var: str = "airbnb_count_500m",55 with_fit_stats: bool = True) -> dict:56 """57 OLS: log_rent ~ airbnb_var + bedrooms + bathrooms + city_FE.58 Returns dict with coefficient, std error, p-value (and N, R² when59 with_fit_stats is True).60 """61 needed = ["log_rent", airbnb_var, "bedrooms", "bathrooms", "city"]62 sub = data[needed].dropna()63 if len(sub) < 30:64 out = {"label": label, "coef": np.nan, "se": np.nan, "pval": np.nan}65 if with_fit_stats:66 out.update({"N": len(sub), "R2": np.nan})67 return out68 city_dummies = pd.get_dummies(sub["city"], prefix="city", drop_first=True, dtype=float)69 X = pd.concat([70 sub[[airbnb_var, "bedrooms", "bathrooms"]].astype(float).reset_index(drop=True),71 city_dummies.reset_index(drop=True),72 ], axis=1)73 X = sm.add_constant(X.astype(float))74 y = sub["log_rent"].astype(float).reset_index(drop=True)75 model = sm.OLS(y, X).fit(cov_type="HC1")76 out = {77 "label": label,78 "coef": model.params[airbnb_var],79 "se": model.bse[airbnb_var],80 "pval": model.pvalues[airbnb_var],81 }82 if with_fit_stats:83 out.update({"N": int(model.nobs), "R2": model.rsquared})84 return out858687# ═════════════════════════════════════════════════════════════════════════════88# TABLE 1: Winsorization / outlier robustness89# ═════════════════════════════════════════════════════════════════════════════9091def table_outlier_robustness(df: pd.DataFrame) -> None:92 print("\n--- Winsorization robustness ---")9394 # (a) Full sample95 res_full = run_ols_city_fe(df, "Full sample")9697 # (b) Winsorized at 5th / 95th percentile on monthly_rent98 p5 = df["monthly_rent"].quantile(0.05)99 p95 = df["monthly_rent"].quantile(0.95)100 df_win = df[(df["monthly_rent"] >= p5) & (df["monthly_rent"] <= p95)].copy()101 res_win = run_ols_city_fe(df_win, "Rent winsorized (5/95)")102103 # (c) Excluding top / bottom 1% of airbnb_count_500m104 p1 = df["airbnb_count_500m"].quantile(0.01)105 p99 = df["airbnb_count_500m"].quantile(0.99)106 df_trim = df[(df["airbnb_count_500m"] >= p1) & (df["airbnb_count_500m"] <= p99)].copy()107 res_trim = run_ols_city_fe(df_trim, "Airbnb trimmed (1/99)")108109 outlier_results = [res_full, res_win, res_trim]110 for r in outlier_results:111 print(f" {r['label']:30s} coef={r['coef']:.6f} N={r['N']}")112113 # Build LaTeX table114 lines = []115 lines.append(r"\begin{tabular}{lccc}")116 lines.append(r"\toprule")117 lines.append(r" & (1) Full Sample & (2) Rent 5/95 & (3) Airbnb 1/99 \\")118 lines.append(r"\midrule")119 # Coefficient row120 coef_cells = " & ".join(121 f"{r['coef']:.6f}{_stars(r['pval'])}" for r in outlier_results122 )123 lines.append(f"Airbnb count (500m) & {coef_cells} \\\\")124 # SE row125 se_cells = " & ".join(f"({r['se']:.6f})" for r in outlier_results)126 lines.append(f" & {se_cells} \\\\")127 lines.append(r"\midrule")128 # N129 n_cells = " & ".join(f"{r['N']:,}" for r in outlier_results)130 lines.append(f"N & {n_cells} \\\\")131 # R²132 r2_cells = " & ".join(f"{r['R2']:.4f}" for r in outlier_results)133 lines.append(f"$R^2$ & {r2_cells} \\\\")134 lines.append(r"City FE & Yes & Yes & Yes \\")135 lines.append(r"\bottomrule")136 lines.append(r"\end{tabular}")137 lines.append(r"\begin{tablenotes}\small")138 lines.append(r"\item \textit{Notes:} Robust standard errors (HC1) in parentheses. "139 r"* $p<0.10$, ** $p<0.05$, *** $p<0.01$. "140 r"Column (2) winsorizes monthly rent at the 5th and 95th percentiles. "141 r"Column (3) trims the top and bottom 1\% of Airbnb listing counts.")142 lines.append(r"\end{tablenotes}")143144 (TABLE_DIR / "robustness_outliers.tex").write_text("\n".join(lines) + "\n",145 encoding="utf-8")146 print(f"Saved {TABLE_DIR / 'robustness_outliers.tex'}")147148149# ═════════════════════════════════════════════════════════════════════════════150# TABLE 2: Subsample analysis151# ═════════════════════════════════════════════════════════════════════════════152153def table_subsamples(df: pd.DataFrame) -> list[dict]:154 print("\n--- Subsample analysis ---")155156 # Identify Montreal157 montreal_mask = df["city"].str.contains("Montr", case=False, na=False)158159 sub_results = []160161 # (a) Montreal only162 res_mtl = run_ols_city_fe(df[montreal_mask], "Montreal only")163 sub_results.append(res_mtl)164165 # (b) Outside Montreal166 res_non = run_ols_city_fe(df[~montreal_mask], "Outside Montreal")167 sub_results.append(res_non)168169 # (c) Apartments only170 apt_mask = df["building_type"].str.contains("Apartment", case=False, na=False)171 res_apt = run_ols_city_fe(df[apt_mask], "Apartments only")172 sub_results.append(res_apt)173174 # (d) Houses only175 house_mask = df["building_type"].str.contains("House", case=False, na=False)176 res_house = run_ols_city_fe(df[house_mask], "Houses only")177 sub_results.append(res_house)178179 for r in sub_results:180 print(f" {r['label']:30s} coef={r['coef']:.6f} N={r['N']}")181182 # Build LaTeX table183 lines2 = []184 lines2.append(r"\begin{tabular}{lcccc}")185 lines2.append(r"\toprule")186 lines2.append(r" & (1) Montreal & (2) Outside Mtl & (3) Apartments & (4) Houses \\")187 lines2.append(r"\midrule")188 coef_cells2 = " & ".join(189 f"{r['coef']:.6f}{_stars(r['pval'])}" if not pd.isna(r['coef']) else "---"190 for r in sub_results191 )192 lines2.append(f"Airbnb count (500m) & {coef_cells2} \\\\")193 se_cells2 = " & ".join(194 f"({r['se']:.6f})" if not pd.isna(r['se']) else ""195 for r in sub_results196 )197 lines2.append(f" & {se_cells2} \\\\")198 lines2.append(r"\midrule")199 n_cells2 = " & ".join(200 f"{r['N']:,}" for r in sub_results201 )202 lines2.append(f"N & {n_cells2} \\\\")203 r2_cells2 = " & ".join(204 f"{r['R2']:.4f}" if not pd.isna(r['R2']) else "---"205 for r in sub_results206 )207 lines2.append(f"$R^2$ & {r2_cells2} \\\\")208 lines2.append(r"City FE & Yes & Yes & Yes & Yes \\")209 lines2.append(r"\bottomrule")210 lines2.append(r"\end{tabular}")211 lines2.append(r"\begin{tablenotes}\small")212 lines2.append(r"\item \textit{Notes:} Robust standard errors (HC1) in parentheses. "213 r"* $p<0.10$, ** $p<0.05$, *** $p<0.01$. "214 r"All specifications include city fixed effects, bedrooms, and bathrooms "215 r"as controls.")216 lines2.append(r"\end{tablenotes}")217218 (TABLE_DIR / "robustness_subsamples.tex").write_text("\n".join(lines2) + "\n",219 encoding="utf-8")220 print(f"Saved {TABLE_DIR / 'robustness_subsamples.tex'}")221222 return sub_results223224225# ═════════════════════════════════════════════════════════════════════════════226# FIGURE 1: Hexbin heatmap — Montreal area227# ═════════════════════════════════════════════════════════════════════════════228229def figure_heatmap(df: pd.DataFrame) -> None:230 print("\n--- Hexbin heatmap (Montreal) ---")231232 mtl_area = df[233 (df["lat"] >= 45.4) & (df["lat"] <= 45.6)234 & (df["lon"] >= -73.8) & (df["lon"] <= -73.5)235 ].dropna(subset=["lat", "lon", "airbnb_count_500m"])236237 fig, ax = plt.subplots(figsize=(8, 6))238 hb = ax.hexbin(239 mtl_area["lon"], mtl_area["lat"],240 C=mtl_area["airbnb_count_500m"],241 reduce_C_function=np.mean,242 gridsize=40, cmap="YlOrRd", mincnt=1,243 )244 fig.colorbar(hb, ax=ax, label="Mean Airbnb count (500m)")245 ax.set_xlabel("Longitude")246 ax.set_ylabel("Latitude")247 ax.set_title("Airbnb Density — Montreal Area")248 fig.savefig(FIG_DIR / "rent_airbnb_heatmap.pdf")249 plt.close(fig)250 print(f"Saved {FIG_DIR / 'rent_airbnb_heatmap.pdf'}")251252253# ═════════════════════════════════════════════════════════════════════════════254# FIGURE 2: Coefficient robustness plot255# ═════════════════════════════════════════════════════════════════════════════256257def figure_coefficient_robustness(df: pd.DataFrame, sub_results: list[dict]) -> None:258 print("\n--- Coefficient robustness plot ---")259260 res_mtl, res_non, res_apt = sub_results[0], sub_results[1], sub_results[2]261262 specs = []263 # Baseline (500m)264 specs.append(run_ols_city_fe(df, "Baseline (500m)"))265266 # Buffer variants267 for buf, label in [("250m", "250m buffer"), ("1km", "1km buffer"), ("2km", "2km buffer")]:268 var = f"airbnb_count_{buf}"269 if var in df.columns:270 specs.append(run_ols_city_fe(df, label, airbnb_var=var,271 with_fit_stats=False))272273 # Subsamples274 specs.append(res_mtl)275 specs.append(res_non)276 specs.append(res_apt)277278 # Filter out NaN results279 specs = [s for s in specs if not pd.isna(s.get("coef", np.nan))]280281 labels = [s["label"] for s in specs]282 coefs = np.array([s["coef"] for s in specs])283 ses = np.array([s["se"] for s in specs])284285 fig, ax = plt.subplots(figsize=(7, 0.5 * len(specs) + 2))286 y_pos = np.arange(len(specs))287 ax.errorbar(288 coefs, y_pos, xerr=1.96 * ses,289 fmt="o", color="#2166ac", ecolor="#92c5de", capsize=4,290 markersize=6, elinewidth=1.5,291 )292 ax.axvline(0, color="grey", linestyle="--", linewidth=0.8)293 ax.set_yticks(y_pos)294 ax.set_yticklabels(labels)295 ax.set_xlabel(r"Coefficient on Airbnb count ($\beta$)")296 ax.set_title("Coefficient Robustness Across Specifications")297 ax.invert_yaxis()298 fig.savefig(FIG_DIR / "coefficient_robustness.pdf")299 plt.close(fig)300 print(f"Saved {FIG_DIR / 'coefficient_robustness.pdf'}")301302303# ═════════════════════════════════════════════════════════════════════════════304# FIGURE 3: Mean rent by Airbnb count quintiles305# ═════════════════════════════════════════════════════════════════════════════306307def figure_rent_by_bins(df: pd.DataFrame) -> None:308 print("\n--- Rent by Airbnb quintiles ---")309310 sub_q = df[["monthly_rent", "airbnb_count_500m"]].dropna()311 sub_q["quintile"] = pd.qcut(312 sub_q["airbnb_count_500m"], q=5, labels=False, duplicates="drop"313 )314 quintile_means = sub_q.groupby("quintile")["monthly_rent"].mean()315 quintile_labels = []316 for q in sorted(sub_q["quintile"].unique()):317 lo = sub_q.loc[sub_q["quintile"] == q, "airbnb_count_500m"].min()318 hi = sub_q.loc[sub_q["quintile"] == q, "airbnb_count_500m"].max()319 quintile_labels.append(f"Q{int(q)+1}\n[{lo:.0f}-{hi:.0f}]")320321 fig, ax = plt.subplots(figsize=(7, 5))322 bars = ax.bar(323 range(len(quintile_means)),324 quintile_means.values,325 color="#2166ac",326 edgecolor="white",327 width=0.65,328 )329 # Add value labels on bars330 for bar, val in zip(bars, quintile_means.values):331 ax.text(332 bar.get_x() + bar.get_width() / 2, bar.get_height() + 10,333 f"${val:,.0f}", ha="center", va="bottom", fontsize=9,334 )335 ax.set_xticks(range(len(quintile_means)))336 ax.set_xticklabels(quintile_labels)337 ax.set_xlabel("Airbnb Count (500m) Quintile")338 ax.set_ylabel("Mean Monthly Rent ($)")339 ax.set_title("Mean Rent by Airbnb Exposure Quintile")340 fig.savefig(FIG_DIR / "rent_by_airbnb_bins.pdf")341 plt.close(fig)342 print(f"Saved {FIG_DIR / 'rent_by_airbnb_bins.pdf'}")343344345def main() -> None:346 print("Loading data ...")347 df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))348 print(f"Full sample: {len(df):,} observations")349350 table_outlier_robustness(df)351 sub_results = table_subsamples(df)352 figure_heatmap(df)353 figure_coefficient_robustness(df, sub_results)354 figure_rent_by_bins(df)355356 print("\n=== 10_robustness_tables_figures.py complete ===")357358359if __name__ == "__main__":360 main()361