# Author: Simon-Pierre Boucher — contact@spboucher.ai """ 12_neighbourhood_fe_inference.py -------------------------------- Added 2026-10-17 (editorial-board revision). Two questions left open by the original design: 1. Are the results robust to neighbourhood fixed effects finer than the city? The rental listings carry a postal code; its first three characters define the forward sortation area (FSA, Canada Post), a neighbourhood-scale unit (333 FSAs in the sample, 180 with at least ten listings). Within-FSA variation in the Airbnb count is a much stricter comparison than within-city variation and absorbs most time-invariant neighbourhood amenities that could drive both Airbnb location and rents. 2. Which standard errors? The exposure variable is a spatial moving average, so neighbouring rental listings have mechanically correlated regressors and plausibly correlated errors. We compare HC1, city-clustered, FSA-clustered and Conley (1999) spatial HAC standard errors with a Bartlett kernel at several distance cutoffs. Specifications (table neighbourhood_fe.tex): (F1) city FE (reference, identical to Model 1c) (F2) FSA FE (F3) FSA FE, 250 m count (the most local exposure measure) (F4) FSA FE, extensive + intensive margin (any Airbnb within 500 m, count) (F5) FSA FE, log(1 + count) (F6) FSA FE, count + interior size (subsample reporting size) Inference table (inference_comparison.tex): coefficient on the 500 m count under city FE and FSA FE with HC1, cluster-city, cluster-FSA and Conley standard errors (1, 2, 5, 10 km cutoffs). Outputs: results/tables/neighbourhood_fe.tex results/tables/inference_comparison.tex results/tables/neighbourhood_fe.json """ import json import sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) import numpy as np # noqa: E402 import pandas as pd # noqa: E402 import statsmodels.api as sm # noqa: E402 from scipy import stats # noqa: E402 from src.config import MERGED_ANALYSIS, TABLE_DIR, CHUNK_SIZE, require # noqa: E402 from src.geo import haversine_matrix # noqa: E402 from src.latex_tables import significance_star as _star, label # noqa: E402 PROCESSED_HINT = "Run scripts/04_merge_data.py first." CONLEY_CUTOFFS_KM = [1.0, 2.0, 5.0, 10.0] # ───────────────────────────────────────────────────────────────────────────── # Data preparation # ───────────────────────────────────────────────────────────────────────────── def prepare(df: pd.DataFrame) -> tuple[pd.DataFrame, list[str]]: d = df.copy() d["fsa"] = d["postal_code"].astype(str).str[:3].str.upper() d["log1p_airbnb_500m"] = np.log1p(d["airbnb_count_500m"]) d["any_airbnb_500m"] = (d["airbnb_count_500m"] > 0).astype(float) 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 complete_cases(d: pd.DataFrame, x_main: list[str], bt_cols: list[str], fe: str, extra: list[str] | None = None) -> pd.DataFrame: controls = ["bedrooms", "bathrooms"] + (extra or []) needed = ["log_rent", "lat", "lon", "city", "fsa"] + x_main + controls + bt_cols sub = d[needed].dropna().reset_index(drop=True) # drop singleton groups (they contribute nothing to within-group variation) counts = sub[fe].map(sub[fe].value_counts()) sub = sub[counts > 1].reset_index(drop=True) return sub # ───────────────────────────────────────────────────────────────────────────── # Estimation with absorbed fixed effects (FWL) + alternative covariances # ───────────────────────────────────────────────────────────────────────────── def within_transform(sub: pd.DataFrame, cols: list[str], fe: str) -> np.ndarray: g = sub.groupby(fe) return (sub[cols] - g[cols].transform("mean")).values.astype(float) def fit_fe(sub: pd.DataFrame, x_cols: list[str], fe: str) -> dict: """ OLS of log_rent on x_cols with fixed effects for `fe` absorbed by the within transformation. Returns coefficients, residuals, bread, HC1 SEs, cluster SEs (city, FSA) and the R-squared of the full dummy model. """ Xw = within_transform(sub, x_cols, fe) yw = within_transform(sub, ["log_rent"], fe).ravel() n, k = Xw.shape n_groups = sub[fe].nunique() bread = np.linalg.inv(Xw.T @ Xw) beta = bread @ Xw.T @ yw e = yw - Xw @ beta # full-model R2: fitted = group means + Xw beta fitted = sub["log_rent"].values - e ss_tot = ((sub["log_rent"] - sub["log_rent"].mean()) ** 2).sum() r2 = 1 - (e ** 2).sum() / ss_tot within_r2 = 1 - (e ** 2).sum() / (yw ** 2).sum() dof_adj = n / (n - k - n_groups) # HC1-type correction with absorbed dummies S = Xw * e[:, None] hc1 = np.sqrt(np.diag(bread @ (S.T @ S) @ bread) * dof_adj) out = {"beta": dict(zip(x_cols, beta)), "se_hc1": dict(zip(x_cols, hc1)), "n": n, "k": k, "n_groups": n_groups, "r2": r2, "within_r2": within_r2, "_Xw": Xw, "_e": e, "_bread": bread, "_S": S} return out def cluster_se(fit: dict, groups: pd.Series) -> np.ndarray: S, bread = fit["_S"], fit["_bread"] G = groups.nunique() n, k = S.shape meat = np.zeros((k, k)) for _, idx in groups.groupby(groups).indices.items(): sg = S[idx].sum(axis=0) meat += np.outer(sg, sg) adj = (G / (G - 1)) * ((n - 1) / (n - k)) return np.sqrt(np.diag(bread @ meat @ bread) * adj) def conley_se(fit: dict, lat: np.ndarray, lon: np.ndarray, cutoff_km: float) -> np.ndarray: """ Conley (1999) spatial HAC covariance with a Bartlett kernel K(d) = max(0, 1 - d / cutoff) on great-circle distances; pairs with d = 0 (units in the same building) receive full weight. """ S, bread = fit["_S"], fit["_bread"] n, k = S.shape meat = np.zeros((k, k)) for start in range(0, n, CHUNK_SIZE): end = min(start + CHUNK_SIZE, n) d = haversine_matrix(lat[start:end], lon[start:end], lat, lon) K = np.clip(1.0 - d / cutoff_km, 0.0, None) meat += S[start:end].T @ (K @ S) return np.sqrt(np.diag(bread @ meat @ bread)) # ───────────────────────────────────────────────────────────────────────────── # Tables # ───────────────────────────────────────────────────────────────────────────── def star(b, se): p = 2 * (1 - stats.norm.cdf(abs(b / se))) return _star(p) def table_specs(results: list[dict], headers: list[str], fe_rows: list[str], display: list[str]) -> None: n = len(results) lines = [r"\begin{tabular}{l" + "c" * n + "}", r"\toprule", " & ".join([""] + [f"\\textbf{{{h}}}" for h in headers]) + r" \\", "Dep.\\ var: & \\multicolumn{" + str(n) + r"}{c}{\textit{log rent}} \\", r"\midrule"] for v in display: c, s_ = [], [] for r in results: if v in r["beta"]: b, se = r["beta"][v], r["se_hc1"][v] c.append(f"{b:.4f}{star(b, se)}") s_.append(f"({se:.4f})") else: c.append(""); s_.append("") lines.append(f"{label(v)} & " + " & ".join(c) + r" \\") lines.append(" & " + " & ".join(s_) + r" \\[4pt]") lines += [r"\midrule", "Fixed effects & " + " & ".join(fe_rows) + r" \\", "FE groups & " + " & ".join(f"{r['n_groups']:,}" for r in results) + r" \\", "Observations & " + " & ".join(f"{r['n']:,}" for r in results) + r" \\", "R$^2$ & " + " & ".join(f"{r['r2']:.4f}" for r in results) + r" \\", "Within R$^2$ & " + " & ".join(f"{r['within_r2']:.4f}" for r in results) + r" \\", r"\bottomrule", r"\end{tabular}", r"\parbox{\textwidth}{\footnotesize Heteroskedasticity-robust (HC1) standard " r"errors in parentheses. All columns control for bedrooms, bathrooms and " r"building type; fixed effects are absorbed (within transformation) and " r"singleton groups are dropped. FSA = forward sortation area (first three " r"characters of the postal code). Column (F3) uses the 250\,m count; column " r"(F6) adds interior size on the subsample reporting it. $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.10$.}"] (TABLE_DIR / "neighbourhood_fe.tex").write_text("\n".join(lines) + "\n", encoding="utf-8") print(f" -> saved {TABLE_DIR / 'neighbourhood_fe.tex'}") def table_inference(rows: list[dict]) -> None: """rows: dicts with keys spec, beta, se_hc1, se_city, se_fsa, se_conley{cut}. Layout (2026-10-17): one column block per specification, one row per covariance estimator, reporting the standard error and the t-statistic.""" cuts = CONLEY_CUTOFFS_KM n = len(rows) lines = [r"\begin{tabular}{l" + "cc" * n + "}", r"\toprule", " & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{\\textbf{{{r['spec']}}}}}" for r in rows) + r" \\", " & " + " & ".join(r"\multicolumn{2}{c}{$\hat\beta = " + f"{r['beta']:.4f}" + "$}" for r in rows) + r" \\", r"\cmidrule(lr){2-3}" + (r"\cmidrule(lr){4-5}" if n > 1 else ""), "Standard errors & " + " & ".join(["SE & $t$"] * n) + r" \\", r"\midrule"] def row(label, key_fn): cells = [] for r in rows: se = key_fn(r) cells += [f"({se:.4f})", f"{r['beta'] / se:.1f}"] lines.append(f"{label} & " + " & ".join(cells) + r" \\") row("Heteroskedasticity-robust (HC1)", lambda r: r["se_hc1"]) row("Clustered by city", lambda r: r["se_city"]) row("Clustered by FSA", lambda r: r["se_fsa"]) for c in cuts: row(f"Conley spatial HAC, {c:g}\\,km cutoff", lambda r, c=c: r["se_conley"][c]) lines += [r"\midrule", "City clusters & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{{r['n_city']}}}" for r in rows) + r" \\", "FSA clusters & " + " & ".join(f"\\multicolumn{{2}}{{c}}{{{r['n_fsa']}}}" for r in rows) + r" \\", r"\bottomrule", r"\end{tabular}", r"\parbox{\textwidth}{\footnotesize Coefficient on the Airbnb count within " r"500\,m from the specification in the column (controls: bedrooms, bathrooms, " r"building type; fixed effects absorbed, singleton groups dropped). Each row " r"reports the standard error (in parentheses) and the implied $t$-statistic " r"under a different covariance estimator: heteroskedasticity-robust (HC1); " r"clustered by city and by FSA \citep{cameron2015practitioner}; and Conley " r"spatial HAC with a Bartlett kernel declining linearly to zero at the stated " r"great-circle cutoff \citep{conley1999gmm, kelejian2007hac}. Listings in the " r"same building (identical coordinates) receive full weight in the Conley " r"estimator.}"] (TABLE_DIR / "inference_comparison.tex").write_text("\n".join(lines) + "\n", encoding="utf-8") print(f" -> saved {TABLE_DIR / 'inference_comparison.tex'}") # ───────────────────────────────────────────────────────────────────────────── # Main # ───────────────────────────────────────────────────────────────────────────── def main() -> None: print("=" * 72) print("12 NEIGHBOURHOOD FIXED EFFECTS AND SPATIAL INFERENCE") print("=" * 72) df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT)) d, bt_cols = prepare(df) print(f"\nSample: {len(d):,} rows | FSAs: {d['fsa'].nunique()} | cities: {d['city'].nunique()}") rec: dict = {} # ── Specifications F1–F5 ───────────────────────────────────────────────── print("\n--- Specifications ---") sub_city = complete_cases(d, ["airbnb_count_500m"], bt_cols, "city") f1 = fit_fe(sub_city, ["airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "city") sub_fsa = complete_cases(d, ["airbnb_count_500m"], bt_cols, "fsa") f2 = fit_fe(sub_fsa, ["airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa") # any_airbnb / log1p are derived from the count: compute on the FSA sample sub_fsa["any_airbnb_500m"] = (sub_fsa["airbnb_count_500m"] > 0).astype(float) sub_fsa["log1p_airbnb_500m"] = np.log1p(sub_fsa["airbnb_count_500m"]) sub_250 = complete_cases(d, ["airbnb_count_250m"], bt_cols, "fsa") f3 = fit_fe(sub_250, ["airbnb_count_250m", "bedrooms", "bathrooms"] + bt_cols, "fsa") f4 = fit_fe(sub_fsa, ["any_airbnb_500m", "airbnb_count_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa") f5 = fit_fe(sub_fsa, ["log1p_airbnb_500m", "bedrooms", "bathrooms"] + bt_cols, "fsa") sub_size = complete_cases(d, ["airbnb_count_500m"], bt_cols, "fsa", extra=["size_100sqft"]) f6 = fit_fe(sub_size, ["airbnb_count_500m", "bedrooms", "bathrooms", "size_100sqft"] + bt_cols, "fsa") for name, f in [("F1 city FE", f1), ("F2 FSA FE", f2), ("F3 FSA FE 250m", f3), ("F4 FSA FE any+count", f4), ("F5 FSA FE log1p", f5), ("F6 FSA FE + size", f6)]: key = list(f["beta"].keys())[0] print(f" [{name}] N={f['n']:,} groups={f['n_groups']} R2={f['r2']:.4f} " f"within R2={f['within_r2']:.4f} beta({key})={f['beta'][key]:.6f} " f"SE={f['se_hc1'][key]:.6f}") rec[name] = {"n": f["n"], "groups": f["n_groups"], "r2": f["r2"], "within_r2": f["within_r2"], "beta": {k: float(v) for k, v in f["beta"].items()}, "se_hc1": {k: float(v) for k, v in f["se_hc1"].items()}} # bivariate reference (no controls, no FE) for the coefficient-movement discussion biv = sm.OLS(d["log_rent"], sm.add_constant(d["airbnb_count_500m"].astype(float))).fit(cov_type="HC1") rec["bivariate"] = {"beta": float(biv.params["airbnb_count_500m"]), "se": float(biv.bse["airbnb_count_500m"]), "r2": float(biv.rsquared), "n": int(biv.nobs)} print(f" [bivariate] beta={rec['bivariate']['beta']:.6f} R2={rec['bivariate']['r2']:.4f}") table_specs([f1, f2, f3, f4, f5, f6], ["(F1)", "(F2)", "(F3)", "(F4)", "(F5)", "(F6)"], ["City", "FSA", "FSA", "FSA", "FSA", "FSA"], ["any_airbnb_500m", "airbnb_count_500m", "airbnb_count_250m", "log1p_airbnb_500m", "bedrooms", "bathrooms", "size_100sqft"]) # ── Inference comparison ───────────────────────────────────────────────── print("\n--- Inference comparison ---") rows = [] for spec, f, sub in [("City FE (Model 1c)", f1, sub_city), ("FSA FE (F2)", f2, sub_fsa)]: lat = sub["lat"].values.astype(float) lon = sub["lon"].values.astype(float) merged = sub[["city", "fsa"]].reset_index(drop=True) se_city = cluster_se(f, merged["city"]) se_fsa = cluster_se(f, merged["fsa"]) se_con = {c: conley_se(f, lat, lon, c) for c in CONLEY_CUTOFFS_KM} j = 0 # airbnb_count_500m is the first regressor row = {"spec": spec, "beta": float(f["beta"]["airbnb_count_500m"]), "se_hc1": float(f["se_hc1"]["airbnb_count_500m"]), "se_city": float(se_city[j]), "se_fsa": float(se_fsa[j]), "se_conley": {c: float(se_con[c][j]) for c in CONLEY_CUTOFFS_KM}, "n_city": int(merged["city"].nunique()), "n_fsa": int(merged["fsa"].nunique())} rows.append(row) print(f" [{spec}] beta={row['beta']:.6f} HC1={row['se_hc1']:.6f} city={row['se_city']:.6f} " f"({row['n_city']} cl.) fsa={row['se_fsa']:.6f} ({row['n_fsa']} cl.) " + " ".join(f"Conley{c:g}km={row['se_conley'][c]:.6f}" for c in CONLEY_CUTOFFS_KM)) rec["inference"] = rows table_inference(rows) (TABLE_DIR / "neighbourhood_fe.json").write_text( json.dumps(rec, indent=2, default=float), encoding="utf-8") print("\n" + "=" * 72) print("12 DONE") print("=" * 72) if __name__ == "__main__": main()