UQO Working Paper No. 5 — Airbnb, residential rents and housing market pressure.
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""312_neighbourhood_fe_inference.py4--------------------------------5Added 2026-10-17 (editorial-board revision).67Two questions left open by the original design:891. Are the results robust to neighbourhood fixed effects finer than the city?10 The rental listings carry a postal code; its first three characters define11 the forward sortation area (FSA, Canada Post), a neighbourhood-scale unit12 (333 FSAs in the sample, 180 with at least ten listings). Within-FSA13 variation in the Airbnb count is a much stricter comparison than14 within-city variation and absorbs most time-invariant neighbourhood15 amenities that could drive both Airbnb location and rents.16172. Which standard errors? The exposure variable is a spatial moving average,18 so neighbouring rental listings have mechanically correlated regressors19 and plausibly correlated errors. We compare HC1, city-clustered,20 FSA-clustered and Conley (1999) spatial HAC standard errors with a21 Bartlett kernel at several distance cutoffs.2223Specifications (table neighbourhood_fe.tex):24 (F1) city FE (reference, identical to Model 1c)25 (F2) FSA FE26 (F3) FSA FE, 250 m count (the most local exposure measure)27 (F4) FSA FE, extensive + intensive margin (any Airbnb within 500 m, count)28 (F5) FSA FE, log(1 + count)29 (F6) FSA FE, count + interior size (subsample reporting size)3031Inference table (inference_comparison.tex): coefficient on the 500 m count32under city FE and FSA FE with HC1, cluster-city, cluster-FSA and Conley33standard errors (1, 2, 5, 10 km cutoffs).3435Outputs:36 results/tables/neighbourhood_fe.tex37 results/tables/inference_comparison.tex38 results/tables/neighbourhood_fe.json39"""4041import json42import sys43from pathlib import Path4445sys.path.insert(0, str(Path(__file__).resolve().parent.parent))4647import numpy as np # noqa: E40248import pandas as pd # noqa: E40249import statsmodels.api as sm # noqa: E40250from scipy import stats # noqa: E4025152from src.config import MERGED_ANALYSIS, TABLE_DIR, CHUNK_SIZE, require # noqa: E40253from src.geo import haversine_matrix # noqa: E40254from src.latex_tables import significance_star as _star, label # noqa: E4025556PROCESSED_HINT = "Run scripts/04_merge_data.py first."57CONLEY_CUTOFFS_KM = [1.0, 2.0, 5.0, 10.0]585960# ─────────────────────────────────────────────────────────────────────────────61# Data preparation62# ─────────────────────────────────────────────────────────────────────────────6364def prepare(df: pd.DataFrame) -> tuple[pd.DataFrame, list[str]]:65 d = df.copy()66 d["fsa"] = d["postal_code"].astype(str).str[:3].str.upper()67 d["log1p_airbnb_500m"] = np.log1p(d["airbnb_count_500m"])68 d["any_airbnb_500m"] = (d["airbnb_count_500m"] > 0).astype(float)69 d["size_100sqft"] = d["size_interior_sqft"] / 100.070 bt = pd.get_dummies(d["building_type"], prefix="bt", drop_first=True, dtype=float)71 d = pd.concat([d, bt], axis=1)72 return d, list(bt.columns)737475def complete_cases(d: pd.DataFrame, x_main: list[str], bt_cols: list[str],76 fe: str, extra: list[str] | None = None) -> pd.DataFrame:77 controls = ["bedrooms", "bathrooms"] + (extra or [])78 needed = ["log_rent", "lat", "lon", "city", "fsa"] + x_main + controls + bt_cols79 sub = d[needed].dropna().reset_index(drop=True)80 # drop singleton groups (they contribute nothing to within-group variation)81 counts = sub[fe].map(sub[fe].value_counts())82 sub = sub[counts > 1].reset_index(drop=True)83 return sub848586# ─────────────────────────────────────────────────────────────────────────────87# Estimation with absorbed fixed effects (FWL) + alternative covariances88# ─────────────────────────────────────────────────────────────────────────────8990def within_transform(sub: pd.DataFrame, cols: list[str], fe: str) -> np.ndarray:91 g = sub.groupby(fe)92 return (sub[cols] - g[cols].transform("mean")).values.astype(float)939495def fit_fe(sub: pd.DataFrame, x_cols: list[str], fe: str) -> dict:96 """97 OLS of log_rent on x_cols with fixed effects for `fe` absorbed by the98 within transformation. Returns coefficients, residuals, bread, HC1 SEs,99 cluster SEs (city, FSA) and the R-squared of the full dummy model.100 """101 Xw = within_transform(sub, x_cols, fe)102 yw = within_transform(sub, ["log_rent"], fe).ravel()103 n, k = Xw.shape104 n_groups = sub[fe].nunique()105 bread = np.linalg.inv(Xw.T @ Xw)106 beta = bread @ Xw.T @ yw107 e = yw - Xw @ beta108 # full-model R2: fitted = group means + Xw beta109 fitted = sub["log_rent"].values - e110 ss_tot = ((sub["log_rent"] - sub["log_rent"].mean()) ** 2).sum()111 r2 = 1 - (e ** 2).sum() / ss_tot112 within_r2 = 1 - (e ** 2).sum() / (yw ** 2).sum()113 dof_adj = n / (n - k - n_groups) # HC1-type correction with absorbed dummies114 S = Xw * e[:, None]115 hc1 = np.sqrt(np.diag(bread @ (S.T @ S) @ bread) * dof_adj)116 out = {"beta": dict(zip(x_cols, beta)), "se_hc1": dict(zip(x_cols, hc1)),117 "n": n, "k": k, "n_groups": n_groups, "r2": r2, "within_r2": within_r2,118 "_Xw": Xw, "_e": e, "_bread": bread, "_S": S}119 return out120121122def cluster_se(fit: dict, groups: pd.Series) -> np.ndarray:123 S, bread = fit["_S"], fit["_bread"]124 G = groups.nunique()125 n, k = S.shape126 meat = np.zeros((k, k))127 for _, idx in groups.groupby(groups).indices.items():128 sg = S[idx].sum(axis=0)129 meat += np.outer(sg, sg)130 adj = (G / (G - 1)) * ((n - 1) / (n - k))131 return np.sqrt(np.diag(bread @ meat @ bread) * adj)132133134def conley_se(fit: dict, lat: np.ndarray, lon: np.ndarray, cutoff_km: float) -> np.ndarray:135 """136 Conley (1999) spatial HAC covariance with a Bartlett kernel137 K(d) = max(0, 1 - d / cutoff) on great-circle distances; pairs with d = 0138 (units in the same building) receive full weight.139 """140 S, bread = fit["_S"], fit["_bread"]141 n, k = S.shape142 meat = np.zeros((k, k))143 for start in range(0, n, CHUNK_SIZE):144 end = min(start + CHUNK_SIZE, n)145 d = haversine_matrix(lat[start:end], lon[start:end], lat, lon)146 K = np.clip(1.0 - d / cutoff_km, 0.0, None)147 meat += S[start:end].T @ (K @ S)148 return np.sqrt(np.diag(bread @ meat @ bread))149150151# ─────────────────────────────────────────────────────────────────────────────152# Tables153# ─────────────────────────────────────────────────────────────────────────────154155def star(b, se):156 p = 2 * (1 - stats.norm.cdf(abs(b / se)))157 return _star(p)158159160def table_specs(results: list[dict], headers: list[str], fe_rows: list[str],161 display: list[str]) -> None:162 n = len(results)163 lines = [r"\begin{tabular}{l" + "c" * n + "}", r"\toprule",164 " & ".join([""] + [f"\\textbf{{{h}}}" for h in headers]) + r" \\",165 "Dep.\\ var: & \\multicolumn{" + str(n) + r"}{c}{\textit{log rent}} \\",166 r"\midrule"]167 for v in display:168 c, s_ = [], []169 for r in results:170 if v in r["beta"]:171 b, se = r["beta"][v], r["se_hc1"][v]172 c.append(f"{b:.4f}{star(b, se)}")173 s_.append(f"({se:.4f})")174 else:175 c.append(""); s_.append("")176 lines.append(f"{label(v)} & " + " & ".join(c) + r" \\")177 lines.append(" & " + " & ".join(s_) + r" \\[4pt]")178 lines += [r"\midrule",179 "Fixed effects & " + " & ".join(fe_rows) + r" \\",180 "FE groups & " + " & ".join(f"{r['n_groups']:,}" for r in results) + r" \\",181 "Observations & " + " & ".join(f"{r['n']:,}" for r in results) + r" \\",182 "R$^2$ & " + " & ".join(f"{r['r2']:.4f}" for r in results) + r" \\",183 "Within R$^2$ & " + " & ".join(f"{r['within_r2']:.4f}" for r in results) + r" \\",184 r"\bottomrule", r"\end{tabular}",185 r"\parbox{\textwidth}{\footnotesize Heteroskedasticity-robust (HC1) standard "186 r"errors in parentheses. All columns control for bedrooms, bathrooms and "187 r"building type; fixed effects are absorbed (within transformation) and "188 r"singleton groups are dropped. FSA = forward sortation area (first three "189 r"characters of the postal code). Column (F3) uses the 250\,m count; column "190 r"(F6) adds interior size on the subsample reporting it. $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"]191 (TABLE_DIR / "neighbourhood_fe.tex").write_text("\n".join(lines) + "\n", encoding="utf-8")192 print(f" -> saved {TABLE_DIR / 'neighbourhood_fe.tex'}")193194195def table_inference(rows: list[dict]) -> None:196 """rows: dicts with keys spec, beta, se_hc1, se_city, se_fsa, se_conley{cut}.197 Layout (2026-10-17): one column block per specification, one row per198 covariance estimator, reporting the standard error and the t-statistic."""199 cuts = CONLEY_CUTOFFS_KM200 n = len(rows)201 lines = [r"\begin{tabular}{l" + "cc" * n + "}", r"\toprule",202 " & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{\\textbf{{{r['spec']}}}}}" for r in rows) + r" \\",203 " & " + " & ".join(r"\multicolumn{2}{c}{$\hat\beta = " + f"{r['beta']:.4f}" + "$}" for r in rows) + r" \\",204 r"\cmidrule(lr){2-3}" + (r"\cmidrule(lr){4-5}" if n > 1 else ""),205 "Standard errors & " + " & ".join(["SE & $t$"] * n) + r" \\", r"\midrule"]206207 def row(label, key_fn):208 cells = []209 for r in rows:210 se = key_fn(r)211 cells += [f"({se:.4f})", f"{r['beta'] / se:.1f}"]212 lines.append(f"{label} & " + " & ".join(cells) + r" \\")213214 row("Heteroskedasticity-robust (HC1)", lambda r: r["se_hc1"])215 row("Clustered by city", lambda r: r["se_city"])216 row("Clustered by FSA", lambda r: r["se_fsa"])217 for c in cuts:218 row(f"Conley spatial HAC, {c:g}\\,km cutoff", lambda r, c=c: r["se_conley"][c])219 lines += [r"\midrule",220 "City clusters & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{{r['n_city']}}}" for r in rows) + r" \\",221 "FSA clusters & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{{r['n_fsa']}}}" for r in rows) + r" \\",222 r"\bottomrule", r"\end{tabular}",223 r"\parbox{\textwidth}{\footnotesize Coefficient on the Airbnb count within "224 r"500\,m from the specification in the column (controls: bedrooms, bathrooms, "225 r"building type; fixed effects absorbed, singleton groups dropped). Each row "226 r"reports the standard error (in parentheses) and the implied $t$-statistic "227 r"under a different covariance estimator: heteroskedasticity-robust (HC1); "228 r"clustered by city and by FSA \citep{cameron2015practitioner}; and Conley "229 r"spatial HAC with a Bartlett kernel declining linearly to zero at the stated "230 r"great-circle cutoff \citep{conley1999gmm, kelejian2007hac}. Listings in the "231 r"same building (identical coordinates) receive full weight in the Conley "232 r"estimator.}"]233 (TABLE_DIR / "inference_comparison.tex").write_text("\n".join(lines) + "\n", encoding="utf-8")234 print(f" -> saved {TABLE_DIR / 'inference_comparison.tex'}")235236237# ─────────────────────────────────────────────────────────────────────────────238# Main239# ─────────────────────────────────────────────────────────────────────────────240241def main() -> None:242 print("=" * 72)243 print("12 NEIGHBOURHOOD FIXED EFFECTS AND SPATIAL INFERENCE")244 print("=" * 72)245 df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))246 d, bt_cols = prepare(df)247 print(f"\nSample: {len(d):,} rows | FSAs: {d['fsa'].nunique()} | cities: {d['city'].nunique()}")248249 rec: dict = {}250251 # ── Specifications F1–F5 ─────────────────────────────────────────────────252 print("\n--- Specifications ---")253 sub_city = complete_cases(d, ["airbnb_count_500m"], bt_cols, "city")254 f1 = fit_fe(sub_city, ["airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "city")255 sub_fsa = complete_cases(d, ["airbnb_count_500m"], bt_cols, "fsa")256 f2 = fit_fe(sub_fsa, ["airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa")257 # any_airbnb / log1p are derived from the count: compute on the FSA sample258 sub_fsa["any_airbnb_500m"] = (sub_fsa["airbnb_count_500m"] > 0).astype(float)259 sub_fsa["log1p_airbnb_500m"] = np.log1p(sub_fsa["airbnb_count_500m"])260 sub_250 = complete_cases(d, ["airbnb_count_250m"], bt_cols, "fsa")261 f3 = fit_fe(sub_250, ["airbnb_count_250m", "bedrooms", "bathrooms"] + bt_cols, "fsa")262 f4 = fit_fe(sub_fsa, ["any_airbnb_500m", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa")263 f5 = fit_fe(sub_fsa, ["log1p_airbnb_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa")264 sub_size = complete_cases(d, ["airbnb_count_500m"], bt_cols, "fsa", extra=["size_100sqft"])265 f6 = fit_fe(sub_size, ["airbnb_count_500m", "bedrooms", "bathrooms", "size_100sqft"] + bt_cols, "fsa")266267 for name, f in [("F1 city FE", f1), ("F2 FSA FE", f2), ("F3 FSA FE 250m", f3),268 ("F4 FSA FE any+count", f4), ("F5 FSA FE log1p", f5), ("F6 FSA FE + size", f6)]:269 key = list(f["beta"].keys())[0]270 print(f" [{name}] N={f['n']:,} groups={f['n_groups']} R2={f['r2']:.4f} "271 f"within R2={f['within_r2']:.4f} beta({key})={f['beta'][key]:.6f} "272 f"SE={f['se_hc1'][key]:.6f}")273 rec[name] = {"n": f["n"], "groups": f["n_groups"], "r2": f["r2"],274 "within_r2": f["within_r2"],275 "beta": {k: float(v) for k, v in f["beta"].items()},276 "se_hc1": {k: float(v) for k, v in f["se_hc1"].items()}}277278 # bivariate reference (no controls, no FE) for the coefficient-movement discussion279 biv = sm.OLS(d["log_rent"], sm.add_constant(d["airbnb_count_500m"].astype(float))).fit(cov_type="HC1")280 rec["bivariate"] = {"beta": float(biv.params["airbnb_count_500m"]),281 "se": float(biv.bse["airbnb_count_500m"]), "r2": float(biv.rsquared),282 "n": int(biv.nobs)}283 print(f" [bivariate] beta={rec['bivariate']['beta']:.6f} R2={rec['bivariate']['r2']:.4f}")284285 table_specs([f1, f2, f3, f4, f5, f6], ["(F1)", "(F2)", "(F3)", "(F4)", "(F5)", "(F6)"],286 ["City", "FSA", "FSA", "FSA", "FSA", "FSA"],287 ["any_airbnb_500m", "airbnb_count_500m", "airbnb_count_250m",288 "log1p_airbnb_500m", "bedrooms", "bathrooms", "size_100sqft"])289290 # ── Inference comparison ─────────────────────────────────────────────────291 print("\n--- Inference comparison ---")292 rows = []293 for spec, f, sub in [("City FE (Model 1c)", f1, sub_city), ("FSA FE (F2)", f2, sub_fsa)]:294 lat = sub["lat"].values.astype(float)295 lon = sub["lon"].values.astype(float)296 merged = sub[["city", "fsa"]].reset_index(drop=True)297 se_city = cluster_se(f, merged["city"])298 se_fsa = cluster_se(f, merged["fsa"])299 se_con = {c: conley_se(f, lat, lon, c) for c in CONLEY_CUTOFFS_KM}300 j = 0 # airbnb_count_500m is the first regressor301 row = {"spec": spec, "beta": float(f["beta"]["airbnb_count_500m"]),302 "se_hc1": float(f["se_hc1"]["airbnb_count_500m"]),303 "se_city": float(se_city[j]), "se_fsa": float(se_fsa[j]),304 "se_conley": {c: float(se_con[c][j]) for c in CONLEY_CUTOFFS_KM},305 "n_city": int(merged["city"].nunique()), "n_fsa": int(merged["fsa"].nunique())}306 rows.append(row)307 print(f" [{spec}] beta={row['beta']:.6f} HC1={row['se_hc1']:.6f} city={row['se_city']:.6f} "308 f"({row['n_city']} cl.) fsa={row['se_fsa']:.6f} ({row['n_fsa']} cl.) " +309 " ".join(f"Conley{c:g}km={row['se_conley'][c]:.6f}" for c in CONLEY_CUTOFFS_KM))310 rec["inference"] = rows311 table_inference(rows)312313 (TABLE_DIR / "neighbourhood_fe.json").write_text(314 json.dumps(rec, indent=2, default=float), encoding="utf-8")315 print("\n" + "=" * 72)316 print("12 DONE")317 print("=" * 72)318319320if __name__ == "__main__":321 main()322