Repository Public Research preview
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Checks and corrected counts computed from the STORED v1.2 estimates (no5refits). November 2026 review, batch L1 (items WP03-THEO-01 preliminary,6ADVE-05, EDIT-04, THEO-06, LITE-08, ECON-03 preliminary, ADVE-12, ECON-10,7ADVE-13, ECON-12, ADVE-22, THEO-14, ADVE-09, THEO-07, ADVE-15, ADVE-11,8ECON-09, ADVE-21, THEO-15, ADVE-23, THEO-10, THEO-08, DATA-05, ECON-11).910Writes results/stored_checks/*.csv and the first version of11results/inline_numbers.json (every entry: value + source artefact).1213Usage: python scripts/10_stored_checks.py14"""1516import json17import sys18from pathlib import Path1920import numpy as np21import pandas as pd22from scipy import stats23from sklearn.metrics import r2_score2425sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))26from wp3 import config, data # noqa: E4022728OUT = config.RESULTS_DIR / "stored_checks"29REV = config.RESULTS_DIR / "revision"30TAB = config.RESULTS_DIR / "tables"31INLINE = config.RESULTS_DIR / "inline_numbers.json"32numbers = {}333435def put(key, value, source, note=""):36 numbers[key] = {"value": value, "source": source, "note": note}373839def holm(p):40 p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 0.041 for i, j in enumerate(order):42 v = max(min(1.0, (m - i) * p[j]), prev); adj[j] = v; prev = v43 return adj444546def bh(p):47 p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 1.048 for i in range(m - 1, -1, -1):49 j = order[i]; v = min(min(1.0, m / (i + 1) * p[j]), prev); adj[j] = v; prev = v50 return adj515253def within_decomposition(ml, df):54 """(i) oracle within-state R2 and SSE shares on the stored random test set (THEO-01)."""55 y = ml["y_test"].values56 st = df.loc[ml["X_test"].index, "address_state"].values57 rows = []58 for name, pred in [("OLS", ml["y_pred_ols"]), ("LightGBM", ml["y_pred_lgb"]), ("XGBoost", ml["y_pred_xgb"])]:59 g = pd.Series(st)60 yd = y - pd.Series(y).groupby(g.values).transform("mean").values61 pdm = pred - pd.Series(pred).groupby(g.values).transform("mean").values62 e = pred - y; eb = pd.Series(e).groupby(g.values).transform("mean").values63 sst = ((y - y.mean()) ** 2).sum()64 rows.append({"model": name, "design": "random 80/20 (stored v1.2 predictions)", "n_test": len(y),65 "r2": r2_score(y, pred), "r2_within_state_oracle": r2_score(yd, pdm),66 "between_share_of_variance": 1 - yd.var() / y.var(),67 "sse_within_share": ((e - eb) ** 2).sum() / sst, "sse_between_share": (eb ** 2).sum() / sst,68 "sse_total_share": (e ** 2).sum() / sst})69 out = pd.DataFrame(rows); out.to_csv(OUT / "within_decomposition_random.csv", index=False)70 for _, r in out.iterrows():71 put(f"r2_within_state_random_{r.model.lower()}", round(r.r2_within_state_oracle, 4), "stored_checks/within_decomposition_random.csv")72 return out737475def significance_counts():76 """(ii) counts over the 62 slopes with |t| > 1.96, t(50), Holm, BH (ADVE-05, EDIT-04, THEO-06, LITE-08)."""77 c = pd.read_csv(REV / "cluster_se.csv", index_col=0)78 slopes = c.drop(index="const")79 rows = []80 for se_kind, G in [("hc3", None), ("cluster_zip3", 886), ("cluster_state", 51)]:81 t = slopes[f"t_{se_kind}"].values82 dfree = (G - 1) if G else np.inf83 p = 2 * (1 - (stats.t.cdf(np.abs(t), dfree) if G else stats.norm.cdf(np.abs(t))))84 crit = stats.t.ppf(0.975, dfree) if G else 1.9685 rows.append({"se": se_kind, "n_slopes": len(t), "n_abs_t_gt_196": int((np.abs(t) > 1.96).sum()),86 "critical_t_G_minus_1": crit, "n_abs_t_gt_crit": int((np.abs(t) > crit).sum()),87 "n_holm_5pct": int((holm(p) < 0.05).sum()), "n_bh_5pct": int((bh(p) < 0.05).sum()),88 "n_incl_const_abs_t_gt_196": int((c[f"t_{se_kind}"].abs() > 1.96).sum())})89 out = pd.DataFrame(rows); out.to_csv(OUT / "significance_counts.csv", index=False)90 r = out.set_index("se")91 put("n_slopes", 62, "cluster_se.csv")92 put("n_sig_state_slopes_196", int(r.loc["cluster_state", "n_abs_t_gt_196"]), "stored_checks/significance_counts.csv", "v1.2 said '33 of 63' incl. the constant")93 put("n_sig_state_slopes_t50", int(r.loc["cluster_state", "n_abs_t_gt_crit"]), "stored_checks/significance_counts.csv", "critical t(50) = 2.009")94 put("n_sig_state_slopes_holm", int(r.loc["cluster_state", "n_holm_5pct"]), "stored_checks/significance_counts.csv")95 put("n_sig_state_slopes_bh", int(r.loc["cluster_state", "n_bh_5pct"]), "stored_checks/significance_counts.csv")96 put("n_sig_zip3_slopes_196", int(r.loc["cluster_zip3", "n_abs_t_gt_196"]), "stored_checks/significance_counts.csv")97 put("n_sig_hc3_slopes_196", int(r.loc["hc3", "n_abs_t_gt_196"]), "stored_checks/significance_counts.csv")98 # SE ratio ranges (ADVE-11, DATA-05)99 rz, rs = slopes["ratio_zip3_hc3"], slopes["ratio_state_hc3"]100 put("se_ratio_state_hc3_min", round(float(rs.min()), 3), "cluster_se.csv"); put("se_ratio_state_hc3_max", round(float(rs.max()), 1), "cluster_se.csv")101 put("se_ratio_state_hc3_median", round(float(rs.median()), 1), "cluster_se.csv"); put("se_ratio_zip3_hc3_median", round(float(rz.median()), 1), "cluster_se.csv")102 put("se_ratio_state_hc3_p10", round(float(rs.quantile(0.1)), 1), "cluster_se.csv"); put("se_ratio_state_hc3_p90", round(float(rs.quantile(0.9)), 1), "cluster_se.csv")103 put("n_clusters_zip3", int(c["n_clusters_zip3"].iloc[0]), "cluster_se.csv", "885 ZIP3 prefixes + 1 missing-ZIP code '000'")104 # Kennedy effect with clustered variance for the dummies (THEO-08)105 dummies = ["on_waterfront", "has_pool", "has_garage", "has_central_air", "has_hardwood", "is_condo",106 "tag_foreclosure", "tag_new_construction", "has_hoa", "has_basement", "has_spa", "has_fireplace"]107 kr = []108 for v in dummies:109 b = slopes.loc[v, "coef"]110 for k in ["hc3", "cluster_zip3", "cluster_state"]:111 s = slopes.loc[v, f"se_{k}"]112 kr.append({"variable": v, "se": k, "coef": b, "se_value": s, "naive_pct": 100 * (np.exp(b) - 1),113 "kennedy_pct": 100 * (np.exp(b - 0.5 * s ** 2) - 1), "t": b / s})114 kdf = pd.DataFrame(kr); kdf.to_csv(OUT / "kennedy_clustered.csv", index=False)115 return out116117118def iqr_adjusted():119 """(ii-b) Holm / BH on the stored paired-bootstrap z (ECON-03 preliminary, EDIT-08)."""120 d = pd.read_csv(REV / "iqr_wald_bootstrap.csv")121 d["p_normal_zip3"] = 2 * (1 - stats.norm.cdf(d["z_boot_zip3"].abs()))122 d["p_holm_zip3"] = holm(d["p_normal_zip3"]); d["p_bh_zip3"] = bh(d["p_normal_zip3"])123 d["sig_percentile_ci_zip3"] = (d["ci_low_zip3"] > 0) | (d["ci_high_zip3"] < 0)124 d["sig_independent"] = d["z_independent"].abs() > 1.96125 d[["variable", "difference", "z_independent", "z_boot_zip3", "p_normal_zip3", "p_holm_zip3", "p_bh_zip3",126 "significant_5pct_zip3", "sig_percentile_ci_zip3", "sig_independent"]].to_csv(OUT / "iqr_adjusted_counts.csv", index=False)127 put("iqr_n_sig_independent", int(d["sig_independent"].sum()), "stored_checks/iqr_adjusted_counts.csv", "of 13")128 put("iqr_n_sig_zip3_normal", int(d["significant_5pct_zip3"].sum()), "stored_checks/iqr_adjusted_counts.csv")129 put("iqr_n_sig_zip3_holm", int((d["p_holm_zip3"] < 0.05).sum()), "stored_checks/iqr_adjusted_counts.csv")130 put("iqr_n_sig_zip3_bh", int((d["p_bh_zip3"] < 0.05).sum()), "stored_checks/iqr_adjusted_counts.csv")131 put("iqr_n_sig_zip3_percentile_ci", int(d["sig_percentile_ci_zip3"].sum()), "stored_checks/iqr_adjusted_counts.csv")132 put("iqr_n_draws_v12", int(d["n_draws"].iloc[0]), "iqr_wald_bootstrap.csv")133134135def zip_cluster(df):136 """(iii) the '000' pseudo-cluster (ADVE-12, ECON-10)."""137 z3 = df["address_zipcode"].fillna("000").astype(str).str[:3]138 n000 = int((z3 == "000").sum())139 pd.DataFrame([{"cluster": "000 (missing ZIP)", "n_listings": n000, "n_zip3_prefixes": int(z3.nunique() - 1),140 "n_clusters_total": int(z3.nunique())}]).to_csv(OUT / "zip3_missing_cluster.csv", index=False)141 put("n_listings_missing_zip", n000, "stored_checks/zip3_missing_cluster.csv")142 put("n_zip3_prefixes", int(z3.nunique() - 1), "stored_checks/zip3_missing_cluster.csv")143 put("n_zip5", int(df["address_zipcode"].nunique()), "analytical sample")144145146def moran_sd():147 """(iv) Moran's I SDs with ddof = 0 and 1 (ADVE-13, ECON-12)."""148 m = pd.read_csv(TAB / "morans_i_robustness.csv")149 rows = [{"table": "Table 11 (OLS residuals, 3 x 5,000)", "model": "OLS", "mean_I": m.morans_i.mean(),150 "sd_ddof0": m.morans_i.std(ddof=0), "sd_ddof1": m.morans_i.std(ddof=1), "values": " ".join(f"{v:.4f}" for v in m.morans_i)}]151 mm = pd.read_csv(REV / "moran_xgb_vs_ols.csv")152 for model, g in mm.groupby("model"):153 rows.append({"table": "Table 12 (test residuals, 3 x 5,000)", "model": model, "mean_I": g.I.mean(),154 "sd_ddof0": g.I.std(ddof=0), "sd_ddof1": g.I.std(ddof=1), "values": " ".join(f"{v:.4f}" for v in g.I),155 "z_min": g.z.min(), "z_max": g.z.max()})156 try:157 ext = data.load_result("extended_results.pkl")["morans_i"]158 rows.append({"table": "extended_results.pkl morans_i (1 x 5,000, k=10)", "model": "OLS", "mean_I": float(ext["I"]),159 "sd_ddof0": np.nan, "sd_ddof1": np.nan, "values": f"z={float(ext['z']):.2f} n={ext['n']} k={ext['k']}"})160 except Exception:161 pass162 out = pd.DataFrame(rows); out.to_csv(OUT / "moran_sd_convention.csv", index=False)163 put("moran_ols_mean", round(float(m.morans_i.mean()), 4), "tables/morans_i_robustness.csv")164 put("moran_ols_sd_ddof0", round(float(m.morans_i.std(ddof=0)), 4), "stored_checks/moran_sd_convention.csv")165 put("moran_ols_sd_ddof1", round(float(m.morans_i.std(ddof=1)), 4), "stored_checks/moran_sd_convention.csv")166167168def loso_summaries():169 """(v)-(vi) LOSO summaries and within-state ranges (ADVE-22, THEO-14, ADVE-09, THEO-07, ADVE-15)."""170 l = pd.read_csv(REV / "loso_by_state.csv")171 rows = []172 for variant, g in l.groupby("variant"):173 rows.append({"variant": variant, "median": g.r2.median(), "listing_weighted_mean": np.average(g.r2, weights=g.n_test),174 "unweighted_mean": g.r2.mean(), "p10": g.r2.quantile(0.1), "p90": g.r2.quantile(0.9),175 "n_states_gt_05": int((g.r2 > 0.5).sum()), "n_negative": int((g.r2 < 0).sum()),176 "median_100_trees": g.r2_100_trees.median(),177 "listing_weighted_mean_100_trees": np.average(g.r2_100_trees, weights=g.n_test)})178 out = pd.DataFrame(rows); out.to_csv(OUT / "loso_summaries_v12.csv", index=False)179 f = out[out.variant.str.startswith("Full")].iloc[0]180 put("loso_median_v12", round(float(f["median"]), 3), "stored_checks/loso_summaries_v12.csv")181 put("loso_listing_weighted_mean_v12", round(float(f["listing_weighted_mean"]), 3), "stored_checks/loso_summaries_v12.csv")182 put("loso_n_states_gt_05_v12", int(f["n_states_gt_05"]), "stored_checks/loso_summaries_v12.csv", "v1.2 text said 44 in one place, 43 in another")183 put("loso_p10_v12", round(float(f["p10"]), 3), "stored_checks/loso_summaries_v12.csv"); put("loso_p90_v12", round(float(f["p90"]), 3), "stored_checks/loso_summaries_v12.csv")184 v = pd.read_csv(REV / "variants_two_designs.csv")185 put("within_state_r2_min_variants_v12", round(float(v.r2_geo_within_state.min()), 3), "variants_two_designs.csv")186 put("within_state_r2_max_variants_v12", round(float(v.r2_geo_within_state.max()), 3), "variants_two_designs.csv")187 put("within_state_r2_full_v12", round(float(v.r2_geo_within_state.iloc[0]), 4), "variants_two_designs.csv", "0.6179 -> '0.618'")188 s = pd.read_csv(REV / "state_heterogeneity_summary.csv").set_index("variable")189 put("ln_sqft_state_p10_v12", round(float(s.loc["ln_sqft", "p10"]), 2), "state_heterogeneity_summary.csv")190 put("ln_sqft_state_p90_v12", round(float(s.loc["ln_sqft", "p90"]), 2), "state_heterogeneity_summary.csv")191 put("ln_sqft_state_min_v12", round(float(s.loc["ln_sqft", "min"]), 2), "state_heterogeneity_summary.csv")192 put("ln_sqft_state_max_v12", round(float(s.loc["ln_sqft", "max"]), 2), "state_heterogeneity_summary.csv")193194195def trace_untracked():196 """(vii) 0.547 and README RF 0.542 trace to extended_results.pkl (ECON-09, ADVE-21)."""197 ext = data.load_result("extended_results.pkl")198 g = ext["geo_holdout"]199 rows = [{"quantity": "ten-state holdout, v1.1 run, XGBoost", "value": g["xgb"]["r2"], "artefact": "extended_results.pkl['geo_holdout']['xgb']['r2']", "quoted_as": "0.547 (§5.6.2 'v1.1 run')"},200 {"quantity": "ten-state holdout, v1.1 run, Random Forest", "value": g["rf"]["r2"], "artefact": "extended_results.pkl['geo_holdout']['rf']['r2']", "quoted_as": "0.542 (README)"},201 {"quantity": "ten-state holdout, v1.1 run, OLS (untrimmed)", "value": g["ols"]["r2"], "artefact": "extended_results.pkl['geo_holdout']['ols']['r2']", "quoted_as": "-1.96e6"},202 {"quantity": "random split, RF", "value": ext["r2_rf"], "artefact": "extended_results.pkl['r2_rf']", "quoted_as": "0.784 (Table 17)"},203 {"quantity": "random split, state-FE OLS", "value": ext["r2_state_test"], "artefact": "extended_results.pkl['r2_state_test']", "quoted_as": "0.678 (Table 9)"}]204 pd.DataFrame(rows).to_csv(OUT / "untracked_numbers_traced.csv", index=False)205 put("geo_xgb_v11_run", round(float(g["xgb"]["r2"]), 4), "extended_results.pkl geo_holdout", "round count of that run not stored")206 put("geo_rf_v11_run", round(float(g["rf"]["r2"]), 4), "extended_results.pkl geo_holdout")207 # realised rounds per ablation stage (ECON-09)208 v3 = data.load_result("v3_spatial_results.pkl")209 ab = pd.DataFrame([{"stage": k, "best_iter_random": v["best_iter_random"], "best_iter_geo": v["best_iter_geo"],210 "r2_random": v["r2_random"], "r2_geo": v["r2_geo"]} for k, v in v3["ablation"].items()])211 ab.to_csv(OUT / "ablation_realised_rounds_v11.csv", index=False)212213214def duplicates(ml, df):215 """(viii) exact duplicates across the random split; re-score without them (THEO-15)."""216 keys = ["address_zipcode", "living_area_sqft", "bedrooms", "bathrooms", "year_built", "price"]217 h = pd.util.hash_pandas_object(df[keys].fillna(-1), index=False).values218 test_idx = ml["X_test"].index.values219 is_test = np.zeros(len(df), bool); is_test[test_idx] = True220 train_keys = set(h[~is_test])221 dup_in_train = np.array([k in train_keys for k in h[test_idx]])222 y = ml["y_test"].values223 rows = []224 for name, pred in [("OLS", ml["y_pred_ols"]), ("LightGBM", ml["y_pred_lgb"]), ("XGBoost", ml["y_pred_xgb"])]:225 rows.append({"model": name, "r2_all": r2_score(y, pred), "r2_without_test_rows_duplicated_in_train": r2_score(y[~dup_in_train], pred[~dup_in_train]),226 "r2_duplicated_rows_only": r2_score(y[dup_in_train], pred[dup_in_train]),227 "n_test": len(y), "n_test_with_train_duplicate": int(dup_in_train.sum())})228 out = pd.DataFrame(rows); out.to_csv(OUT / "duplicates_random_split.csv", index=False)229 put("n_test_rows_with_train_duplicate", int(dup_in_train.sum()), "stored_checks/duplicates_random_split.csv")230 put("share_test_rows_with_train_duplicate_pct", round(100 * float(dup_in_train.mean()), 2), "stored_checks/duplicates_random_split.csv")231 for _, r in out.iterrows():232 put(f"r2_random_{r.model.lower()}_without_duplicates", round(r.r2_without_test_rows_duplicated_in_train, 4), "stored_checks/duplicates_random_split.csv")233 return out234235236def zestimate_selection(df):237 """(ix) listings with vs without a Zestimate on observables (ADVE-23, THEO-10)."""238 has = df["zestimate"].notna()239 cols = ["price", "living_area_sqft", "bedrooms", "bathrooms", "age", "lot_size_sqft", "is_condo", "has_pool",240 "tag_new_construction", "tag_foreclosure", "walk_score", "avg_school_rating", "luxury_score", "days_on_zillow"]241 rows = []242 for c in cols:243 a, b = df.loc[has, c].astype(float), df.loc[~has, c].astype(float)244 rows.append({"variable": c, "mean_with_zestimate": a.mean(), "mean_without": b.mean(), "median_with": a.median(),245 "median_without": b.median(), "std_diff": (a.mean() - b.mean()) / np.sqrt((a.var() + b.var()) / 2)})246 reg = df.groupby("region")["zestimate"].apply(lambda s: s.notna().mean())247 ht = df.groupby("is_condo")["zestimate"].apply(lambda s: s.notna().mean())248 out = pd.DataFrame(rows); out.to_csv(OUT / "zestimate_selection.csv", index=False)249 pd.concat([reg.rename("share_with_zestimate"), ht.rename("share_with_zestimate")]).to_csv(OUT / "zestimate_coverage_by_group.csv")250 put("share_with_zestimate_pct", round(100 * float(has.mean()), 1), "stored_checks/zestimate_selection.csv")251 put("n_with_zestimate", int(has.sum()), "stored_checks/zestimate_selection.csv")252253254def inline_basics(df):255 """Inline numbers that are recomputable from the sample (DATA-03)."""256 ld = df["days_on_zillow"].dropna()257 put("days_on_zillow_median", float(ld.median()), "analytical sample"); put("days_on_zillow_p25", float(ld.quantile(0.25)), "analytical sample")258 put("days_on_zillow_p75", float(ld.quantile(0.75)), "analytical sample", "v1.2 said IQR 35-134")259 put("n_tax_rates_distinct", int(df["property_tax_rate"].nunique()), "analytical sample")260 put("share_school_rating_imputed_pct", round(100 * 733929 / 788842, 1), "build/missingness_rules.csv")261 ms = pd.read_csv(config.RESULTS_DIR / "build" / "missingness_rules.csv").set_index("raw_column")262 for c in ["avg_school_rating", "lot_size_sqft", "year_built", "stories", "parking_spaces", "walk_score", "has_garage", "transit_score"]:263 put(f"missing_share_{c}_pct", round(100 * float(ms.loc[c, "missing_share"]), 1), "build/missingness_rules.csv")264 lot = pd.read_csv(config.RESULTS_DIR / "build" / "lot_reconciliation.csv").set_index("class")["n"]265 put("lot_null", int(lot.iloc[0]), "build/lot_reconciliation.csv"); put("lot_zero", int(lot.iloc[1]), "build/lot_reconciliation.csv")266 put("lot_one", int(lot.iloc[2]), "build/lot_reconciliation.csv"); put("lot_negative", int(lot.iloc[4]), "build/lot_reconciliation.csv")267 put("ln_lot_zero_total", int(lot.iloc[5]), "build/lot_reconciliation.csv", "v1.2: '6,329 zero lots' = 6,115 + 211 + 3")268 put("n_implausible", 4970, "build/build_meta.json"); put("n_trimmed_sample", 783872, "build/build_meta.json")269 put("n_sample", 788842, "build/attrition.csv"); put("n_raw", 839313, "build/attrition.csv")270 put("n_duplicate_rows", 6787, "build/duplicates.csv"); put("n_duplicate_groups", 3152, "build/duplicates.csv")271 put("county_fips_null_pct", 97.0, "build/build_meta.json")272 lc = json.loads((config.RESULTS_DIR / "primary" / "legacy_ols_check.json").read_text())273 put("ols_r2_baseline", round(lc["fit_std"]["r2"], 4), "primary/legacy_ols_check.json"); put("ols_f_baseline", round(lc["fit_std"]["f"], 0), "primary/legacy_ols_check.json")274 put("ols_resid_skew", round(lc["fit_std"]["resid_skew"], 3), "primary/legacy_ols_check.json"); put("ols_resid_kurtosis", round(lc["fit_std"]["resid_kurtosis"], 3), "primary/legacy_ols_check.json")275 put("ols_jarque_bera", round(lc["fit_std"]["jarque_bera"], 0), "primary/legacy_ols_check.json")276 put("ols_rmse_baseline_mse_resid", round(lc["checks"]["rmse_model_data"], 4), "model_data.pkl", "sqrt(SSR/(n-k))")277 put("r2_ols_random_test_v12", round(lc["checks"]["r2_ols_test_rebuilt"], 4), "primary/legacy_ols_check.json")278279280def main():281 OUT.mkdir(parents=True, exist_ok=True)282 df = data.load_analytical_sample()283 ml = data.load_result("ml_results.pkl")284 print("(i) within decomposition"); print(within_decomposition(ml, df).round(4).to_string(index=False))285 print("(ii) significance counts"); print(significance_counts().to_string(index=False))286 iqr_adjusted(); zip_cluster(df); moran_sd(); loso_summaries(); trace_untracked()287 print("(viii) duplicates"); print(duplicates(ml, df).round(4).to_string(index=False))288 zestimate_selection(df); inline_basics(df)289 INLINE.write_text(json.dumps(numbers, indent=1, default=float))290 print(f"{len(numbers)} inline numbers -> {INLINE}")291292293if __name__ == "__main__":294 main()295