# Author: Simon-Pierre Boucher — contact@spboucher.ai """ 11_extended_robustness.py ------------------------- Extended robustness checks for the preferred hedonic specification (log_rent ~ airbnb_count_500m + bedrooms + bathrooms + building type + city FE): (R1) Baseline with HC1 standard errors (reference column) (R2) Same coefficients with standard errors clustered by city (R3) Adding interior size (in 100 sq ft) on the subsample reporting it (R4) log(1 + count) functional form for the exposure variable (R5) Ring exposure: count in the 500m-1km annulus alongside the 500m count plus a leave-one-city-out sensitivity figure (dropping each of the largest cities in turn). Outputs: results/tables/extended_robustness.tex figures/leave_one_city_out.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 statsmodels.api as sm # 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." CONTROL_LABELS = { "airbnb_count_500m": r"Airbnb count (500m)", "log1p_airbnb_500m": r"log(1 + Airbnb count 500m)", "airbnb_ring_500m_1km": r"Airbnb count (500m--1km ring)", "bedrooms": "Bedrooms", "bathrooms": "Bathrooms", "size_100sqft": r"Interior size (100 sq ft)", } def prepare(df: pd.DataFrame) -> tuple[pd.DataFrame, list[str]]: """Add derived regressors and building-type dummies; return (frame, bt_cols).""" d = df.copy() d["log1p_airbnb_500m"] = np.log1p(d["airbnb_count_500m"]) d["airbnb_ring_500m_1km"] = d["airbnb_count_1000m"] - d["airbnb_count_500m"] d["size_100sqft"] = d["size_interior_sqft"] / 100.0 bt = pd.get_dummies(d["building_type"], prefix="bt", drop_first=True, dtype=float) d = pd.concat([d, bt], axis=1) return d, list(bt.columns) def fit(d: pd.DataFrame, x_main: list[str], bt_cols: list[str], label: str, cluster: bool = False, extra_controls: list[str] | None = None): """ OLS of log_rent on x_main + bedrooms/bathrooms (+extras) + building type + city FE, with HC1 or city-clustered standard errors. """ controls = ["bedrooms", "bathrooms"] + (extra_controls or []) needed = ["log_rent", "city"] + x_main + controls sub = d[needed + bt_cols].dropna().reset_index(drop=True) city_fe = pd.get_dummies(sub["city"], prefix="city", drop_first=True, dtype=float) X = sm.add_constant(pd.concat( [sub[x_main + controls + bt_cols].astype(float), city_fe], axis=1 ).astype(float)) y = sub["log_rent"].astype(float) if cluster: res = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": sub["city"]}) else: res = sm.OLS(y, X).fit(cov_type="HC1") print(f" [{label}] N={int(res.nobs):,} R2={res.rsquared:.4f} " f"beta({x_main[0]})={res.params[x_main[0]]:.6f} " f"SE={res.bse[x_main[0]]:.6f}") return res def build_table(results: list, col_headers: list[str], se_row: list[str], display_vars: list[str]) -> None: """Write the extended-robustness table fragment.""" n = len(results) lines: list[str] = [] lines.append(r"\begin{tabular}{l" + "c" * n + "}") lines.append(r"\toprule") lines.append(" & ".join([""] + [f"\\textbf{{{h}}}" for h in col_headers]) + r" \\") lines.append( " & ".join(["Dep.\\ var:"] + [r"\textit{log\_rent}"] * n) + r" \\" ) lines.append(r"\midrule") for var in display_vars: cells_c, cells_s = [], [] for res in results: if var in res.params.index: b, se, p = res.params[var], res.bse[var], res.pvalues[var] cells_c.append(f"{b:.4f}{_star(p)}") cells_s.append(f"({se:.4f})") else: cells_c.append("") cells_s.append("") label = CONTROL_LABELS.get(var, var.replace("_", r"\_")) lines.append(f"{label} & " + " & ".join(cells_c) + r" \\") lines.append(" & " + " & ".join(cells_s) + r" \\[4pt]") lines.append(r"\midrule") lines.append("Standard errors & " + " & ".join(se_row) + r" \\") lines.append("Building type \\& city FE & " + " & ".join(["Yes"] * n) + r" \\") lines.append("Observations & " + " & ".join(f"{int(r.nobs):,}" for r in results) + r" \\") lines.append("R$^2$ & " + " & ".join(f"{r.rsquared:.4f}" for r in results) + r" \\") lines.append(r"\bottomrule") lines.append(r"\end{tabular}") lines.append( r"\parbox{\textwidth}{\footnotesize Standard errors in parentheses: " r"HC1 except column (R2), which clusters by city. All columns control " r"for bedrooms, bathrooms, building type, and city fixed effects; " r"column (R3) adds interior size (per 100 sq ft) on the subsample " r"reporting it. Column (R5) includes the 500m count and the count in " r"the 500m--1km annulus jointly. " r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}" ) out = TABLE_DIR / "extended_robustness.tex" out.write_text("\n".join(lines) + "\n", encoding="utf-8") print(f" -> saved {out}") def leave_one_city_out(d: pd.DataFrame, bt_cols: list[str], n_cities: int = 8) -> None: """Coefficient on the 500m count when each large city is excluded in turn.""" print("\n--- Leave-one-city-out ---") top = d["city"].value_counts().head(n_cities).index.tolist() res_full = fit(d, ["airbnb_count_500m"], bt_cols, "Full sample") beta_full = res_full.params["airbnb_count_500m"] se_full = res_full.bse["airbnb_count_500m"] labels, betas, ses = [], [], [] for city in top: res = fit(d[d["city"] != city], ["airbnb_count_500m"], bt_cols, f"excl. {city}") labels.append(f"excl. {city}") betas.append(res.params["airbnb_count_500m"]) ses.append(res.bse["airbnb_count_500m"]) betas = np.array(betas) ses = np.array(ses) fig, ax = plt.subplots(figsize=(7, 0.45 * len(labels) + 2)) y_pos = np.arange(len(labels)) ax.errorbar(betas, y_pos, xerr=1.96 * ses, fmt="o", color="#2166ac", ecolor="#92c5de", capsize=4, markersize=6, elinewidth=1.5) ax.axvline(beta_full, color="firebrick", linestyle="--", linewidth=1.2, label=f"Full sample ({beta_full:.4f})") ax.axvspan(beta_full - 1.96 * se_full, beta_full + 1.96 * se_full, color="firebrick", alpha=0.08) ax.axvline(0, color="grey", linestyle=":", linewidth=0.8) ax.set_yticks(y_pos) ax.set_yticklabels(labels) ax.set_xlabel(r"Coefficient on Airbnb count within 500m ($\beta$)") ax.set_title("Leave-One-City-Out Sensitivity") ax.legend(fontsize=8, loc="best") ax.invert_yaxis() fig.tight_layout() out = FIG_DIR / "leave_one_city_out.pdf" fig.savefig(out, dpi=300) plt.close(fig) print(f" -> saved {out}") def main() -> None: print("=" * 72) print("11 EXTENDED ROBUSTNESS") print("=" * 72) df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT)) d, bt_cols = prepare(df) print(f"\nSample: {len(d):,} rows | " f"interior size non-missing: {d['size_100sqft'].notna().sum():,}") print("\n--- Specifications R1-R5 ---") res_r1 = fit(d, ["airbnb_count_500m"], bt_cols, "R1 baseline HC1") res_r2 = fit(d, ["airbnb_count_500m"], bt_cols, "R2 clustered", cluster=True) res_r3 = fit(d, ["airbnb_count_500m"], bt_cols, "R3 + size", extra_controls=["size_100sqft"]) res_r4 = fit(d, ["log1p_airbnb_500m"], bt_cols, "R4 log(1+count)") res_r5 = fit(d, ["airbnb_count_500m", "airbnb_ring_500m_1km"], bt_cols, "R5 ring") build_table( [res_r1, res_r2, res_r3, res_r4, res_r5], ["(R1)", "(R2)", "(R3)", "(R4)", "(R5)"], ["HC1", "City cluster", "HC1", "HC1", "HC1"], display_vars=[ "airbnb_count_500m", "log1p_airbnb_500m", "airbnb_ring_500m_1km", "bedrooms", "bathrooms", "size_100sqft", ], ) leave_one_city_out(d, bt_cols) print("\n" + "=" * 72) print("11 DONE") print("=" * 72) if __name__ == "__main__": main()