spb/wp5_uqo

1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""311_extended_robustness.py4-------------------------5Extended robustness checks for the preferred hedonic specification6(log_rent ~ airbnb_count_500m + bedrooms + bathrooms + building type + city FE):78  (R1) Baseline with HC1 standard errors (reference column)9  (R2) Same coefficients with standard errors clustered by city10  (R3) Adding interior size (in 100 sq ft) on the subsample reporting it11  (R4) log(1 + count) functional form for the exposure variable12  (R5) Ring exposure: count in the 500m-1km annulus alongside the 500m count1314plus a leave-one-city-out sensitivity figure (dropping each of the largest15cities in turn).1617Outputs:18    results/tables/extended_robustness.tex19    figures/leave_one_city_out.pdf20"""2122import sys23from pathlib import Path2425sys.path.insert(0, str(Path(__file__).resolve().parent.parent))2627import matplotlib28matplotlib.use("Agg")29import 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 _star  # noqa: E4023637PROCESSED_HINT = "Run scripts/04_merge_data.py first."3839CONTROL_LABELS = {40    "airbnb_count_500m": r"Airbnb count (500m)",41    "log1p_airbnb_500m": r"log(1 + Airbnb count 500m)",42    "airbnb_ring_500m_1km": r"Airbnb count (500m--1km ring)",43    "bedrooms": "Bedrooms",44    "bathrooms": "Bathrooms",45    "size_100sqft": r"Interior size (100 sq ft)",46}474849def prepare(df: pd.DataFrame) -> tuple[pd.DataFrame, list[str]]:50    """Add derived regressors and building-type dummies; return (frame, bt_cols)."""51    d = df.copy()52    d["log1p_airbnb_500m"] = np.log1p(d["airbnb_count_500m"])53    d["airbnb_ring_500m_1km"] = d["airbnb_count_1000m"] - d["airbnb_count_500m"]54    d["size_100sqft"] = d["size_interior_sqft"] / 100.05556    bt = pd.get_dummies(d["building_type"], prefix="bt", drop_first=True, dtype=float)57    d = pd.concat([d, bt], axis=1)58    return d, list(bt.columns)596061def fit(d: pd.DataFrame, x_main: list[str], bt_cols: list[str],62        label: str, cluster: bool = False, extra_controls: list[str] | None = None):63    """64    OLS of log_rent on x_main + bedrooms/bathrooms (+extras) + building type65    + city FE, with HC1 or city-clustered standard errors.66    """67    controls = ["bedrooms", "bathrooms"] + (extra_controls or [])68    needed = ["log_rent", "city"] + x_main + controls69    sub = d[needed + bt_cols].dropna().reset_index(drop=True)7071    city_fe = pd.get_dummies(sub["city"], prefix="city", drop_first=True, dtype=float)72    X = sm.add_constant(pd.concat(73        [sub[x_main + controls + bt_cols].astype(float), city_fe], axis=174    ).astype(float))75    y = sub["log_rent"].astype(float)7677    if cluster:78        res = sm.OLS(y, X).fit(cov_type="cluster",79                               cov_kwds={"groups": sub["city"]})80        res.n_clusters = int(sub["city"].nunique())81        print(f"  [{label}]  city clusters = {res.n_clusters}")82    else:83        res = sm.OLS(y, X).fit(cov_type="HC1")8485    print(f"  [{label}]  N={int(res.nobs):,}  R2={res.rsquared:.4f}  "86          f"beta({x_main[0]})={res.params[x_main[0]]:.6f}  "87          f"SE={res.bse[x_main[0]]:.6f}")88    return res899091def build_table(results: list, col_headers: list[str], se_row: list[str],92                display_vars: list[str]) -> None:93    """Write the extended-robustness table fragment."""94    n = len(results)95    lines: list[str] = []96    lines.append(r"\begin{tabular}{l" + "c" * n + "}")97    lines.append(r"\toprule")98    lines.append(" & ".join([""] + [f"\\textbf{{{h}}}" for h in col_headers]) + r" \\")99    lines.append(100        " & ".join(["Dep.\\ var:"] + [r"\textit{log rent}"] * n) + r" \\"101    )102    lines.append(r"\midrule")103104    for var in display_vars:105        cells_c, cells_s = [], []106        for res in results:107            if var in res.params.index:108                b, se, p = res.params[var], res.bse[var], res.pvalues[var]109                cells_c.append(f"{b:.4f}{_star(p)}")110                cells_s.append(f"({se:.4f})")111            else:112                cells_c.append("")113                cells_s.append("")114        label = CONTROL_LABELS.get(var, var.replace("_", r"\_"))115        lines.append(f"{label} & " + " & ".join(cells_c) + r" \\")116        lines.append(" & " + " & ".join(cells_s) + r" \\[4pt]")117118    lines.append(r"\midrule")119    lines.append("Standard errors & " + " & ".join(se_row) + r" \\")120    lines.append("Clusters & " + " & ".join(121        f"{r.n_clusters}" if hasattr(r, "n_clusters") else "" for r in results) + r" \\")122    lines.append("Building type \\& city FE & " + " & ".join(["Yes"] * n) + r" \\")123    lines.append("Observations & "124                 + " & ".join(f"{int(r.nobs):,}" for r in results) + r" \\")125    lines.append("R$^2$ & "126                 + " & ".join(f"{r.rsquared:.4f}" for r in results) + r" \\")127    lines.append(r"\bottomrule")128    lines.append(r"\end{tabular}")129    lines.append(130        r"\parbox{\textwidth}{\footnotesize Standard errors in parentheses: "131        r"HC1 except column (R2), which clusters by city. All columns control "132        r"for bedrooms, bathrooms, building type, and city fixed effects; "133        r"column (R3) adds interior size (per 100 sq ft) on the subsample "134        r"reporting it. Column (R5) includes the 500m count and the count in "135        r"the 500m--1km annulus jointly. "136        r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"137    )138139    out = TABLE_DIR / "extended_robustness.tex"140    out.write_text("\n".join(lines) + "\n", encoding="utf-8")141    print(f"  -> saved {out}")142143144def leave_one_city_out(d: pd.DataFrame, bt_cols: list[str], n_cities: int = 8) -> None:145    """Coefficient on the 500m count when each large city is excluded in turn."""146    print("\n--- Leave-one-city-out ---")147    top = d["city"].value_counts().head(n_cities).index.tolist()148149    res_full = fit(d, ["airbnb_count_500m"], bt_cols, "Full sample")150    beta_full = res_full.params["airbnb_count_500m"]151    se_full = res_full.bse["airbnb_count_500m"]152153    labels, betas, ses = [], [], []154    for city in top:155        res = fit(d[d["city"] != city], ["airbnb_count_500m"], bt_cols,156                  f"excl. {city}")157        labels.append(f"excl. {city}")158        betas.append(res.params["airbnb_count_500m"])159        ses.append(res.bse["airbnb_count_500m"])160161    betas = np.array(betas)162    ses = np.array(ses)163164    fig, ax = plt.subplots(figsize=(7, 0.45 * len(labels) + 2))165    y_pos = np.arange(len(labels))166    ax.errorbar(betas, y_pos, xerr=1.96 * ses,167                fmt="o", color="#2166ac", ecolor="#92c5de", capsize=4,168                markersize=6, elinewidth=1.5)169    ax.axvline(beta_full, color="firebrick", linestyle="--", linewidth=1.2,170               label=f"Full sample ({beta_full:.4f})")171    ax.axvspan(beta_full - 1.96 * se_full, beta_full + 1.96 * se_full,172               color="firebrick", alpha=0.08)173    ax.axvline(0, color="grey", linestyle=":", linewidth=0.8)174    ax.set_yticks(y_pos)175    ax.set_yticklabels(labels)176    ax.set_xlabel(r"Coefficient on Airbnb count within 500m ($\beta$)")177    ax.set_title("Leave-One-City-Out Sensitivity")178    ax.legend(fontsize=8, loc="best")179    ax.invert_yaxis()180    fig.tight_layout()181182    out = FIG_DIR / "leave_one_city_out.pdf"183    fig.savefig(out, dpi=300)184    plt.close(fig)185    print(f"  -> saved {out}")186187188def main() -> None:189    print("=" * 72)190    print("11  EXTENDED ROBUSTNESS")191    print("=" * 72)192193    df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))194    d, bt_cols = prepare(df)195    print(f"\nSample: {len(d):,} rows | "196          f"interior size non-missing: {d['size_100sqft'].notna().sum():,}")197198    print("\n--- Specifications R1-R5 ---")199    res_r1 = fit(d, ["airbnb_count_500m"], bt_cols, "R1 baseline HC1")200    res_r2 = fit(d, ["airbnb_count_500m"], bt_cols, "R2 clustered", cluster=True)201    res_r3 = fit(d, ["airbnb_count_500m"], bt_cols, "R3 + size",202                 extra_controls=["size_100sqft"])203    res_r4 = fit(d, ["log1p_airbnb_500m"], bt_cols, "R4 log(1+count)")204    res_r5 = fit(d, ["airbnb_count_500m", "airbnb_ring_500m_1km"], bt_cols,205                 "R5 ring")206207    build_table(208        [res_r1, res_r2, res_r3, res_r4, res_r5],209        ["(R1)", "(R2)", "(R3)", "(R4)", "(R5)"],210        ["HC1", "City cluster", "HC1", "HC1", "HC1"],211        display_vars=[212            "airbnb_count_500m", "log1p_airbnb_500m", "airbnb_ring_500m_1km",213            "bedrooms", "bathrooms", "size_100sqft",214        ],215    )216217    leave_one_city_out(d, bt_cols)218219    print("\n" + "=" * 72)220    print("11  DONE")221    print("=" * 72)222223224if __name__ == "__main__":225    main()226