1"""Phase F — figures generated only from pipeline outputs (data/metadata, data/final, results/econometrics).2Writes results/figures/*.{pdf,png} and copies PDFs to paper/figures/. Run: python scripts/visualization/make_figures.py"""3import json, shutil4from pathlib import Path5import matplotlib; matplotlib.use("Agg")6import matplotlib.pyplot as plt7from matplotlib.patches import FancyBboxPatch8import numpy as np, pandas as pd, yaml910ROOT = Path(__file__).resolve().parents[2]; META = ROOT / "data" / "metadata"; FINAL = ROOT / "data" / "final"; ECO = ROOT / "results" / "econometrics"11FIG = ROOT / "results" / "figures"; PFIG = ROOT / "paper" / "figures"; FIG.mkdir(parents=True, exist_ok=True); PFIG.mkdir(parents=True, exist_ok=True)12PAL = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"]; INK = "#0b0b0b"; INK2 = "#52514e"; GRID = "#e6e5e1"13plt.rcParams.update({"font.family": "serif", "font.size": 9, "axes.titlesize": 9.5, "axes.labelsize": 9, "legend.fontsize": 8, "xtick.labelsize": 8, "ytick.labelsize": 8,14 "axes.edgecolor": INK2, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False,15 "lines.linewidth": 1.6, "legend.frameon": False, "figure.dpi": 150, "savefig.dpi": 300, "axes.prop_cycle": matplotlib.cycler(color=PAL)})16LABELS = yaml.safe_load((ROOT / "scripts" / "harmonization" / "taxonomy.yaml").read_text())17LAB = {leaf: v["label"] for grp in LABELS.values() for leaf, v in grp.items()}; GRP = {leaf: g for g, grp in LABELS.items() for leaf in grp}18FAMLAB = {"pct_premium": "% premium (binary attribute, log price)", "semi_elasticity": "% per unit (log price, level X)", "elasticity": "elasticity (log-log)"}192021def save(fig, name):22 fig.tight_layout(); fig.savefig(FIG / f"{name}.pdf", bbox_inches="tight"); fig.savefig(FIG / f"{name}.png", bbox_inches="tight"); shutil.copy(FIG / f"{name}.pdf", PFIG / f"{name}.pdf"); plt.close(fig); print("figure", name)232425def fig01_prisma():26 c = json.load(open(META / "prisma_counts.json"))27 fig, ax = plt.subplots(figsize=(7.4, 7.2)); ax.axis("off"); ax.set_xlim(0, 10); ax.set_ylim(0, 12)28 def box(x, y, w, h, text, col="#eef3fb"):29 ax.add_patch(FancyBboxPatch((x, y), w, h, boxstyle="round,pad=0.02,rounding_size=0.08", fc=col, ec=INK2, lw=0.8)); ax.text(x + w / 2, y + h / 2, text, ha="center", va="center", fontsize=7.6, color=INK)30 def arrow(x1, y1, x2, y2): ax.annotate("", (x2, y2), (x1, y1), arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))31 src = c["by_source"]32 box(0.3, 10.6, 5.4, 1.1, f"Records identified from databases (n = {c['records_identified_raw']:,})\nSemantic Scholar {src.get('semanticscholar_bulk',0):,} · citations {src.get('semanticscholar_citations',0):,} · Crossref {src.get('crossref',0):,}\narXiv {src.get('arxiv',0):,} · OpenAlex {src.get('openalex',0):,}")33 box(6.0, 10.6, 3.7, 1.1, f"Duplicates removed\n(n = {c['records_identified_raw'] - c['records_after_dedup']:,})", "#f6f6f4")34 arrow(3.0, 10.6, 3.0, 9.6)35 box(0.3, 8.5, 5.4, 1.1, f"Unique records screened by rules\n(title/abstract keyword scoring; n = {c['records_after_dedup']:,})")36 box(6.0, 8.5, 3.7, 1.1, f"Excluded by rules\n(no residential price outcome or no\nhedonic/regression signal; n = {c['excluded_by_rules']:,})", "#f6f6f4"); arrow(5.7, 9.05, 6.0, 9.05)37 arrow(3.0, 8.5, 3.0, 7.5)38 box(0.3, 6.4, 5.4, 1.1, f"Stage-1 candidates screened by independent agents\n(title/abstract; n = {c['screened_by_llm']:,})")39 ex = c["llm_exclusion_reasons"]; top = "; ".join(f"{k.split('_')[0]} {v}" for k, v in sorted(ex.items(), key=lambda kv: -kv[1])[:5])40 box(6.0, 6.0, 3.7, 1.5, f"Excluded with reason (n = {c['llm_excluded']:,})\n{top}\nunclear, flagged for review: {c['llm_unclear']:,}", "#f6f6f4"); arrow(5.7, 6.95, 6.0, 6.95)41 arrow(3.0, 6.4, 3.0, 5.4)42 box(0.3, 4.3, 5.4, 1.1, f"Included at stage 1: hedonic coefficient studies n = {c['included_hedonic_stage1']:,}\nML/AVM predictive studies n = {c['included_ml_predictive_stage1']:,} (AI synthesis)")43 arrow(3.0, 4.3, 3.0, 3.3)44 box(0.3, 2.2, 5.4, 1.1, f"Open-access full texts sought (hedonic + unclear): n = {c['fulltext_sought_hedonic_or_unclear']:,}\nretrieved n = {c['fulltext_retrieved']:,}")45 box(6.0, 2.2, 3.7, 1.1, f"Not retrievable legally\n(paywalled / no OA copy / download failed;\nn = {c['fulltext_not_retrieved']:,})", "#f6f6f4"); arrow(5.7, 2.75, 6.0, 2.75)46 arrow(3.0, 2.2, 3.0, 1.2)47 box(0.3, 0.1, 5.4, 1.1, f"Full texts assessed and extracted: n = {c['fulltext_assessed_extracted']:,}\nstudies eligible for coefficient synthesis: n = {c['fulltext_eligible']:,}", "#e3f1ea")48 box(6.0, 0.1, 3.7, 1.1, f"Excluded at full text\n(n = {c['fulltext_ineligible']:,})", "#f6f6f4"); arrow(5.7, 0.65, 6.0, 0.65)49 save(fig, "fig01_prisma_flow")505152def fig02_trends():53 rec = pd.read_csv(META / "records_dedup.csv", low_memory=False, usecols=["record_id", "year", "stage1_candidate"]); llm = pd.read_csv(META / "screening_stage1_llm.csv")54 m = rec.merge(llm[["record_id", "reason_code", "ai_ml_flag"]], on="record_id", how="left")55 yrs = np.arange(1970, 2027)56 inc = m[m.reason_code == "ELIGIBLE_HEDONIC"].groupby("year").size().reindex(yrs, fill_value=0); ml = m[m.reason_code == "ELIGIBLE_ML_PREDICTIVE_ONLY"].groupby("year").size().reindex(yrs, fill_value=0)57 fig, ax = plt.subplots(1, 2, figsize=(8.2, 3.0))58 ax[0].bar(yrs, inc.values, color=PAL[0], width=0.8, label="hedonic coefficient studies (stage-1 include)"); ax[0].bar(yrs, ml.values, bottom=inc.values, color=PAL[1], width=0.8, label="ML / AVM predictive studies")59 ax[0].set_xlabel("publication year"); ax[0].set_ylabel("records included at stage 1"); ax[0].set_title("(a) Annual publication trend of included records"); ax[0].legend(fontsize=7); ax[0].set_xlim(1969, 2027)60 share = (ml / (inc + ml)).rolling(3, center=True, min_periods=1).mean()61 ax[1].plot(yrs, 100 * share.values, color=PAL[1]); ax[1].set_xlabel("publication year"); ax[1].set_ylabel("share of ML/AVM studies (%)"); ax[1].set_title("(b) Rise of machine-learning valuation studies"); ax[1].set_xlim(1995, 2027)62 save(fig, "fig02_publication_trends")636465def fig03_map():66 S = pd.read_csv(FINAL / "studies.csv"); S = S[S.eligible == True]67 llm = pd.read_csv(META / "screening_stage1_llm.csv"); cc = yaml.safe_load((ROOT / "scripts" / "harmonization" / "countries.yaml").read_text())68 g1 = llm[llm.reason_code == "ELIGIBLE_HEDONIC"].country_guess.fillna("unknown").map(lambda c: str(c).split(",")[0].strip()).value_counts()69 fig, ax = plt.subplots(figsize=(8.4, 4.4)); ax.set_xlim(-170, 180); ax.set_ylim(-58, 80); ax.set_xlabel("longitude"); ax.set_ylabel("latitude"); ax.set_aspect("equal")70 for name, (cont, inc, iso, lat, lon) in cc.items():71 n1 = int(g1.get(name, 0)); n2 = int((S.country == name).sum())72 if n1 > 0: ax.scatter(lon, lat, s=6 + 2.0 * n1, color=PAL[0], alpha=0.35, lw=0)73 if n2 > 0: ax.scatter(lon, lat, s=6 + 14 * n2, color=PAL[1], alpha=0.85, lw=0.4, edgecolor=INK)74 if n1 >= 20 or n2 >= 3: ax.text(lon + 3, lat + 2, f"{iso} {n1}/{n2}", fontsize=6.5, color=INK2)75 ax.scatter([], [], color=PAL[0], alpha=0.35, s=60, label="stage-1 included hedonic records (abstract-based country)"); ax.scatter([], [], color=PAL[1], alpha=0.85, s=60, label="studies extracted for coefficient synthesis")76 ax.legend(loc="lower left", fontsize=7); ax.set_title("Geographic coverage (bubble area proportional to counts; label = ISO3 stage-1/extracted)")77 save(fig, "fig03_world_coverage")787980def fig04_heatmap():81 C = pd.read_csv(FINAL / "coefficients_raw.csv"); S = pd.read_csv(FINAL / "studies.csv")82 C = C[C.taxonomy_id != "other"].merge(S[["record_id", "continent"]], on="record_id", how="left"); C["continent"] = C.continent.fillna("unknown")83 st = C.groupby(["taxonomy_id", "continent"]).record_id.nunique().unstack(fill_value=0)84 st = st.loc[st.sum(1).sort_values(ascending=False).index[:30]]85 fig, ax = plt.subplots(figsize=(7.2, 8.0)); im = ax.imshow(st.values, cmap="Blues", aspect="auto")86 ax.set_xticks(range(len(st.columns))); ax.set_xticklabels(st.columns, rotation=30, ha="right"); ax.set_yticks(range(len(st))); ax.set_yticklabels([LAB.get(i, i) for i in st.index]); ax.grid(False)87 for i in range(st.shape[0]):88 for j in range(st.shape[1]):89 v = st.values[i, j]90 if v: ax.text(j, i, str(v), ha="center", va="center", fontsize=6.5, color="white" if v > st.values.max() * 0.6 else INK)91 fig.colorbar(im, ax=ax, label="number of extracted studies reporting the characteristic"); ax.set_title("Variable frequency by continent (extracted studies)")92 save(fig, "fig04_variable_frequency_heatmap")939495def forest(fam, name, title):96 p = pd.read_csv(ECO / "pooled_effects.csv"); p = p[p.family == fam].sort_values("est_re")97 if len(p) == 0: return98 fig, ax = plt.subplots(figsize=(7.0, 0.32 * len(p) + 1.6)); y = np.arange(len(p))99 ax.errorbar(p.est_re, y, xerr=[p.est_re - p.ci_lo, p.ci_hi - p.est_re], fmt="o", color=PAL[0], ms=4, capsize=2, lw=1, label="REML random effects, 95% CI")100 ax.hlines(y, p.pi_lo, p.pi_hi, color=PAL[0], alpha=0.25, lw=5, label="95% prediction interval")101 if p.est_ml.notna().any(): ax.scatter(p.est_ml, y + 0.18, marker="s", s=14, color=PAL[1], label="three-level (study/model) estimate")102 ax.axvline(0, color=INK2, lw=0.8); ax.set_yticks(y); ax.set_yticklabels([f"{LAB.get(t, t)} (k={k}, n={n})" for t, k, n in zip(p.taxonomy_id, p.k_studies, p.n_estimates)])103 ax.set_xlabel(FAMLAB.get(fam, fam)); ax.set_title(title); ax.legend(fontsize=7, loc="lower right")104 lim = np.nanpercentile(np.r_[p.pi_lo, p.pi_hi], [2, 98]); ax.set_xlim(lim[0] - 0.1 * abs(lim[0]) - 1, lim[1] + 0.1 * abs(lim[1]) + 1)105 save(fig, name)106107108def fig05_forests():109 forest("pct_premium", "fig05_forest_pct_premium", "Pooled percentage premiums of binary attributes"); forest("semi_elasticity", "fig06_forest_semi_elasticity", "Pooled semi-elasticities of continuous attributes"); forest("elasticity", "fig07_forest_elasticity", "Pooled elasticities (log-log models)")110111112def fig08_geo_heterogeneity():113 E = pd.read_csv(FINAL / "effects.csv"); top = E.groupby(["family", "taxonomy_id"]).study_id.nunique().sort_values(ascending=False)114 top = [k for k in top.index if top[k] >= 8][:6]115 if not top: return116 fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()117 for a_, (fam, tx) in zip(axes, top):118 s = E[(E.family == fam) & (E.taxonomy_id == tx)]; s = s[s.continent.notna()]119 order = s.groupby("continent").study_id.nunique().sort_values(ascending=False).index120 data = [s[s.continent == c].es.values for c in order]121 a_.boxplot(data, vert=False, labels=[f"{c} (k={s[s.continent == c].study_id.nunique()})" for c in order], widths=0.6, showfliers=False, medianprops=dict(color=PAL[1]))122 a_.axvline(0, color=INK2, lw=0.8); a_.set_title(f"{LAB.get(tx, tx)} — {FAMLAB.get(fam, fam).split(' (')[0]}", fontsize=8); a_.set_xlabel("harmonised effect size")123 for a_ in axes[len(top):]: a_.axis("off")124 save(fig, "fig08_geographic_heterogeneity")125126127def fig09_temporal():128 E = pd.read_csv(FINAL / "effects.csv"); E = E[E.data_mid_year.notna()]129 top = E.groupby(["family", "taxonomy_id"]).study_id.nunique().sort_values(ascending=False); top = [k for k in top.index if top[k] >= 8][:6]130 if not top: return131 fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()132 for a_, (fam, tx) in zip(axes, top):133 s = E[(E.family == fam) & (E.taxonomy_id == tx)]134 a_.scatter(s.data_mid_year, s.es, s=np.clip(8 / np.maximum(s.es_se, 1e-6), 4, 60), color=PAL[0], alpha=0.5, lw=0)135 if len(s) > 5:136 w = (1 / s.es_se ** 2).values; X = np.column_stack([np.ones(len(s)), s.data_mid_year.values - 2000]); b = np.linalg.lstsq(X * np.sqrt(w)[:, None], s.es.values * np.sqrt(w), rcond=None)[0]137 xs = np.linspace(s.data_mid_year.min(), s.data_mid_year.max(), 50); a_.plot(xs, b[0] + b[1] * (xs - 2000), color=PAL[1], label=f"WLS trend {b[1]:+.3f}/yr")138 a_.legend(fontsize=7)139 a_.axhline(0, color=INK2, lw=0.8); a_.set_title(f"{LAB.get(tx, tx)}", fontsize=8); a_.set_xlabel("midpoint of the transaction sample period"); a_.set_ylabel("harmonised effect")140 for a_ in axes[len(top):]: a_.axis("off")141 save(fig, "fig09_temporal_evolution")142143144def fig10_funnel():145 E = pd.read_csv(FINAL / "effects.csv"); top = E.groupby(["family", "taxonomy_id"]).size().sort_values(ascending=False); top = [k for k in top.index if top[k] >= 15][:6]146 if not top: return147 fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()148 for a_, (fam, tx) in zip(axes, top):149 s = E[(E.family == fam) & (E.taxonomy_id == tx)]; w = 1 / s.es_se ** 2; mu = (w * s.es).sum() / w.sum()150 a_.scatter(s.es, s.es_se, s=10, color=PAL[0], alpha=0.6, lw=0); se = np.linspace(0, s.es_se.max(), 50)151 a_.plot(mu - 1.96 * se, se, color=INK2, lw=0.7, ls="--"); a_.plot(mu + 1.96 * se, se, color=INK2, lw=0.7, ls="--"); a_.axvline(mu, color=PAL[1], lw=0.9)152 a_.invert_yaxis(); a_.set_title(f"{LAB.get(tx, tx)} (n={len(s)})", fontsize=8); a_.set_xlabel("harmonised effect"); a_.set_ylabel("standard error")153 for a_ in axes[len(top):]: a_.axis("off")154 save(fig, "fig10_funnel_plots")155156157def fig11_moderators():158 f = ECO / "meta_regression.csv"159 if not f.exists(): return160 m = pd.read_csv(f); m = m[m.term != "intrcpt"]161 if len(m) == 0: return162 keys = m.groupby(["family", "taxonomy_id"]).size().index[:6]163 fig, axes = plt.subplots(1, len(keys), figsize=(2.6 * len(keys) + 1, 3.6), sharey=True); axes = np.atleast_1d(axes)164 for a_, (fam, tx) in zip(axes, keys):165 s = m[(m.family == fam) & (m.taxonomy_id == tx)]; y = np.arange(len(s))166 a_.errorbar(s.est, y, xerr=1.96 * s.se_rve.fillna(s.se), fmt="o", color=PAL[0], ms=4, capsize=2, lw=1); a_.axvline(0, color=INK2, lw=0.8)167 a_.set_yticks(y); a_.set_yticklabels(s.term); a_.set_title(f"{LAB.get(tx, tx)}", fontsize=8); a_.set_xlabel("moderator effect (95% RVE CI)")168 save(fig, "fig11_moderator_effects")169170171def fig12_evidence_density():172 C = pd.read_csv(FINAL / "coefficients_raw.csv"); S = pd.read_csv(FINAL / "studies.csv")173 C = C[C.taxonomy_id != "other"].merge(S[["record_id", "sample_start", "sample_end", "year"]], on="record_id", how="left")174 for c in ("sample_start", "sample_end", "year"):175 C[c] = pd.to_numeric(C[c], errors="coerce"); C.loc[(C[c] < 1900) | (C[c] > 2027), c] = np.nan176 C["mid"] = C[["sample_start", "sample_end"]].mean(axis=1).fillna(C.year - 3); C = C[C.mid.notna()]; C["decade"] = (C.mid // 10 * 10).astype("Int64")177 st = C.groupby(["taxonomy_id", "decade"]).record_id.nunique().unstack(fill_value=0); st = st.loc[st.sum(1).sort_values(ascending=False).index[:30]]178 fig, ax = plt.subplots(figsize=(6.4, 8.0)); im = ax.imshow(st.values, cmap="Oranges", aspect="auto"); ax.grid(False)179 ax.set_xticks(range(len(st.columns))); ax.set_xticklabels([f"{int(d)}s" for d in st.columns]); ax.set_yticks(range(len(st))); ax.set_yticklabels([LAB.get(i, i) for i in st.index])180 for i in range(st.shape[0]):181 for j in range(st.shape[1]):182 if st.values[i, j]: ax.text(j, i, str(st.values[i, j]), ha="center", va="center", fontsize=6.5, color=INK)183 fig.colorbar(im, ax=ax, label="studies"); ax.set_title("Evidence density: characteristic × decade of the transaction data")184 save(fig, "fig12_evidence_density")185186187def fig13_ai_evolution():188 llm = pd.read_csv(META / "screening_stage1_llm.csv"); rec = pd.read_csv(META / "records_dedup.csv", low_memory=False, usecols=["record_id", "year", "title", "abstract"])189 m = llm.merge(rec, on="record_id"); m = m[m.decision == "include"]; t = (m.title.fillna("") + " " + m.abstract.fillna("")).str.lower()190 fams = {"random forest / boosting": r"random forest|gradient boost|xgboost|lightgbm|catboost|boosting", "neural networks / deep learning": r"neural network|deep learning|cnn|lstm|transformer", "computer vision / images": r"street view|satellite|image|computer vision|photo", "text / NLP": r"text|nlp|natural language|sentiment|listing description|word", "spatial econometrics": r"spatial lag|spatial error|geographically weighted|gwr|spatial autoregress|spatial durbin|spatial econometric"}191 yrs = np.arange(2000, 2027); fig, ax = plt.subplots(figsize=(7.2, 3.2))192 for i, (lab, pat) in enumerate(fams.items()):193 cnt = m[t.str.contains(pat)].groupby("year").size().reindex(yrs, fill_value=0).rolling(3, center=True, min_periods=1).mean(); ax.plot(yrs, cnt.values, color=PAL[i], label=lab)194 ax.set_xlabel("publication year"); ax.set_ylabel("included records mentioning the method (3-yr moving average)"); ax.set_title("Methodological evolution in the included literature (title/abstract mentions)"); ax.legend(fontsize=7)195 save(fig, "fig13_ai_methodology_evolution")196197198def fig14_summary():199 p = pd.read_csv(ECO / "pooled_effects.csv"); p = p[p.family.isin(["pct_premium", "semi_elasticity"])].copy(); p["label"] = p.taxonomy_id.map(lambda t: LAB.get(t, t)) + p.family.map({"pct_premium": " (premium %)", "semi_elasticity": " (% per unit)"})200 p = p.sort_values("est_re"); fig, ax = plt.subplots(figsize=(7.4, 0.26 * len(p) + 1.5)); y = np.arange(len(p))201 cols = [PAL[GRP.get(t, "A")[0] in "A" and 0 or (1 if GRP.get(t, "")[0] in "BC" else 2 if GRP.get(t, "")[0] == "D" else 3 if GRP.get(t, "")[0] == "E" else 4)] for t in p.taxonomy_id]202 ax.barh(y, p.est_re, color=cols, height=0.6); ax.errorbar(p.est_re, y, xerr=[p.est_re - p.ci_lo, p.ci_hi - p.est_re], fmt="none", ecolor=INK2, capsize=2, lw=0.8)203 ax.set_yticks(y); ax.set_yticklabels(p.label); ax.axvline(0, color=INK2, lw=0.8); ax.set_xlabel("pooled estimate (REML) with 95% CI"); ax.set_title("Summary of the pooled evidence (colour = taxonomy group)")204 save(fig, "fig14_summary_evidence")205206207ALL = [fig01_prisma, fig02_trends, fig03_map, fig04_heatmap, fig05_forests, fig08_geo_heterogeneity, fig09_temporal, fig10_funnel, fig11_moderators, fig12_evidence_density, fig13_ai_evolution, fig14_summary]208if __name__ == "__main__":209 import sys, traceback210 bad = []211 for f in ALL:212 try: f()213 except Exception as e: bad.append(f.__name__); print("FAILED", f.__name__, repr(e)); traceback.print_exc()214 sys.exit(1 if bad else 0)215