SPB Git

spb/wp5_uqo Public

UQO Working Paper No. 5 — Airbnb, residential rents and housing market pressure.

TeX 53.4% Python 46.5%
8.0 KB · 255 lines python
Raw Blame History
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""308_quantile_models.py4---------------------5Quantile regression analysis for the Airbnb-rent study.67Model 5: Quantile regression across the conditional rent distribution.89  Q_tau(log_rent | X) = alpha_tau + beta_tau * airbnb_count_500m + X'gamma_tau1011for tau in {0.10, 0.25, 0.50, 0.75, 0.90}1213Outputs:14    results/tables/quantile_regression.tex15    figures/quantile_coefficients.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 statsmodels.api as sm  # noqa: E40229from statsmodels.regression.quantile_regression import QuantReg  # noqa: E4023031from src.config import MERGED_ANALYSIS, TABLE_DIR, FIG_DIR, require  # noqa: E40232from src.latex_tables import significance_star as _star  # noqa: E4023334PROCESSED_HINT = "Run scripts/04_merge_data.py first."3536TAUS = [0.10, 0.25, 0.50, 0.75, 0.90]373839def prepare_data(rent: pd.DataFrame):40    """Build dummies and return the complete-case design (Y, X, x_cols, bt_cols)."""41    rent_m = rent.copy()4243    # Building-type dummies44    if "building_type" in rent_m.columns:45        bt_dummies = pd.get_dummies(46            rent_m["building_type"], prefix="bt", drop_first=True, dtype=float47        )48        rent_m = pd.concat([rent_m, bt_dummies], axis=1)49        bt_cols = list(bt_dummies.columns)50    else:51        bt_cols = []5253    # City dummies — use only the top 5 cities to keep manageable54    if "city" in rent_m.columns:55        top5 = rent_m["city"].value_counts().nlargest(5).index.tolist()56        print(f"  Top-5 cities for quantile FE: {top5}")57        rent_m["city_top5"] = rent_m["city"].where(rent_m["city"].isin(top5), other="Other")58        city_dum = pd.get_dummies(59            rent_m["city_top5"], prefix="city", drop_first=True, dtype=float60        )61        rent_m = pd.concat([rent_m, city_dum], axis=1)62        city_fe_cols = list(city_dum.columns)63    else:64        city_fe_cols = []6566    controls = ["bedrooms", "bathrooms"] + bt_cols67    x_cols = ["airbnb_count_500m"] + controls + city_fe_cols6869    # Drop missing70    use_cols = ["log_rent"] + x_cols71    rent_q = rent_m[use_cols].dropna().reset_index(drop=True)72    print(f"  Complete cases: {len(rent_q):,}")7374    Y = rent_q["log_rent"]75    X = sm.add_constant(rent_q[x_cols])76    return Y, X, bt_cols777879def build_table(qr_results: dict, res_ols, bt_cols) -> None:80    """Write the quantile-regression LaTeX table."""81    print("\n  Building LaTeX table ...")8283    # Variables to display (omit city dummies for readability)84    display_vars = ["const", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols8586    n_models = len(TAUS) + 1  # quantiles + OLS87    col_spec = "l" + "c" * n_models88    all_res = [qr_results[t] for t in TAUS] + [res_ols]89    col_headers = [f"$\\tau={t:.2f}$" for t in TAUS] + ["OLS"]9091    lines: list[str] = []92    lines.append(r"\begin{tabular}{" + col_spec + "}")93    lines.append(r"\toprule")94    lines.append(95        " & ".join([""] + [f"\\textbf{{{h}}}" for h in col_headers]) + r" \\"96    )97    lines.append(98        " & ".join(["Dep.\\ var:"] + [r"\textit{log\_rent}"] * n_models) + r" \\"99    )100    lines.append(r"\midrule")101102    for var in display_vars:103        cells_c = []104        cells_s = []105        for res in all_res:106            if var in res.params.index:107                b = res.params[var]108                se = res.bse[var]109                p = res.pvalues[var]110                cells_c.append(f"{b:.4f}{_star(p)}")111                cells_s.append(f"({se:.4f})")112            else:113                cells_c.append("")114                cells_s.append("")115        vn = var.replace("_", r"\_")116        lines.append(f"{vn} & " + " & ".join(cells_c) + r" \\")117        lines.append(f" & " + " & ".join(cells_s) + r" \\[4pt]")118119    lines.append(r"\midrule")120121    # City FE indicator122    lines.append(123        "City FE (top 5) & " + " & ".join(["Yes"] * n_models) + r" \\"124    )125126    # N127    lines.append(128        "Observations & "129        + " & ".join([f"{int(res.nobs)}" for res in all_res])130        + r" \\"131    )132133    # Pseudo-R2 / R2134    r2_cells = []135    for i, res in enumerate(all_res):136        if i < len(TAUS):137            r2_cells.append(f"{res.prsquared:.4f}")138        else:139            r2_cells.append(f"{res.rsquared:.4f}")140    lines.append("(Pseudo-)R$^2$ & " + " & ".join(r2_cells) + r" \\")141142    lines.append(r"\bottomrule")143    lines.append(r"\end{tabular}")144    lines.append(145        r"\parbox{\textwidth}{\footnotesize Standard errors in parentheses. "146        r"OLS uses HC1 robust SE. "147        r"$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$. "148        r"City FE limited to the 5 largest cities (others grouped).}"149    )150151    tex = "\n".join(lines) + "\n"152    out_table = TABLE_DIR / "quantile_regression.tex"153    out_table.write_text(tex, encoding="utf-8")154    print(f"  -> saved {out_table}")155156157def build_plot(qr_model: QuantReg, res_ols) -> None:158    """Coefficient plot: beta_tau(airbnb_count_500m) across a fine tau grid."""159    print("  Building coefficient plot ...")160161    fine_taus = np.arange(0.05, 0.96, 0.05)162    fine_betas = []163    fine_ci_lo = []164    fine_ci_hi = []165166    for tau in fine_taus:167        res = qr_model.fit(q=tau)168        beta = res.params["airbnb_count_500m"]169        se = res.bse["airbnb_count_500m"]170        fine_betas.append(beta)171        fine_ci_lo.append(beta - 1.96 * se)172        fine_ci_hi.append(beta + 1.96 * se)173174    fine_betas = np.array(fine_betas)175    fine_ci_lo = np.array(fine_ci_lo)176    fine_ci_hi = np.array(fine_ci_hi)177178    fig, ax = plt.subplots(figsize=(7, 4.5))179180    # 95% CI band181    ax.fill_between(fine_taus, fine_ci_lo, fine_ci_hi, alpha=0.2, color="steelblue",182                    label="95% CI")183184    # Quantile regression line185    ax.plot(fine_taus, fine_betas, "o-", color="steelblue", markersize=4,186            label=r"$\beta_\tau$ (Quantile Reg.)")187188    # OLS reference189    ols_beta = res_ols.params["airbnb_count_500m"]190    ols_se = res_ols.bse["airbnb_count_500m"]191    ax.axhline(ols_beta, color="firebrick", linestyle="--", linewidth=1.2,192               label=f"OLS estimate ({ols_beta:.5f})")193    ax.axhspan(ols_beta - 1.96 * ols_se, ols_beta + 1.96 * ols_se,194               color="firebrick", alpha=0.08)195196    ax.axhline(0, color="grey", linestyle=":", linewidth=0.7)197    ax.set_xlabel(r"Quantile ($\tau$)")198    ax.set_ylabel(r"$\beta_\tau$ (airbnb\_count\_500m)")199    ax.set_title("Effect of Airbnb Count (500m) on Log Rent Across Quantiles")200    ax.legend(fontsize=8, loc="best")201    fig.tight_layout()202203    fig_path = FIG_DIR / "quantile_coefficients.pdf"204    fig.savefig(fig_path, dpi=300)205    plt.close(fig)206    print(f"  -> saved {fig_path}")207208209def main() -> None:210    print("=" * 72)211    print("08  QUANTILE REGRESSION MODELS")212    print("=" * 72)213214    rent = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))215    print(f"\nRent data: {rent.shape[0]:,} rows")216217    Y, X, bt_cols = prepare_data(rent)218219    qr_model = QuantReg(Y, X)220    qr_results: dict[float, object] = {}221222    print("\n" + "-" * 72)223    print("MODEL 5: Quantile Regressions")224    print("-" * 72)225226    for tau in TAUS:227        res = qr_model.fit(q=tau)228        qr_results[tau] = res229        beta = res.params["airbnb_count_500m"]230        se = res.bse["airbnb_count_500m"]231        p = res.pvalues["airbnb_count_500m"]232        print(233            f"  tau={tau:.2f}  β(airbnb)={beta:.6f}  SE={se:.6f}  "234            f"p={p:.4f}{_star(p)}  pseudo-R2={res.prsquared:.4f}"235        )236237    # Also run OLS for comparison238    res_ols = sm.OLS(Y, X).fit(cov_type="HC1")239    print(240        f"  OLS      β(airbnb)={res_ols.params['airbnb_count_500m']:.6f}  "241        f"SE={res_ols.bse['airbnb_count_500m']:.6f}  "242        f"R2={res_ols.rsquared:.4f}"243    )244245    build_table(qr_results, res_ols, bt_cols)246    build_plot(qr_model, res_ols)247248    print("\n" + "=" * 72)249    print("08  DONE")250    print("=" * 72)251252253if __name__ == "__main__":254    main()255