1"""All paper figures (vector PDF) from processed outputs. Every figure states N, period and source in its caption file2(paper/figures/<name>.caption.txt) which the LaTeX sections quote. Figures whose inputs are missing are skipped with a log line."""3from __future__ import annotations4import json5import numpy as np, pandas as pd, matplotlib6matplotlib.use("Agg")7import matplotlib.pyplot as plt8from matplotlib.colors import TwoSlopeNorm, LinearSegmentedColormap9from src.common import ROOT, P, REGIONS, get_logger, ledger, PROC_DIR, EXP_DIR, BASE_PROC, BASE_EXP1011log = get_logger("reports.figures")12PROC = PROC_DIR; FIG = P["figures"]; FIG.mkdir(exist_ok=True)13C = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4"]; GRAY = "#8a8a85"14SEQ = LinearSegmentedColormap.from_list("seqblue", ["#cde2fb", "#2a78d6", "#0d366b"])15DIV = LinearSegmentedColormap.from_list("divob", ["#eb6834", "#f2f2ef", "#2a78d6"])16SRC = "Source: Fabri-Ka research snapshot (2026-11-01), author's computations."17plt.rcParams.update({"font.family": "serif", "font.size": 9, "axes.spines.top": False, "axes.spines.right": False, "axes.edgecolor": GRAY, "xtick.color": "#444", "ytick.color": "#444",18 "axes.grid": True, "grid.color": "#e6e6e2", "grid.linewidth": 0.5, "axes.axisbelow": True, "legend.frameon": False, "pdf.fonttype": 42})19SHORT = {"Saguenay–Lac-Saint-Jean": "Saguenay–Lac-St-Jean", "Gaspésie–Îles-de-la-Madeleine": "Gaspésie–Îles", "Abitibi-Témiscamingue": "Abitibi-Témis.", "Chaudière-Appalaches": "Chaudière-App.", "Capitale-Nationale": "Capitale-Nat.", "Centre-du-Québec": "Centre-du-Qc", "Bas-Saint-Laurent": "Bas-St-Laurent", "Nord-du-Québec": "Nord-du-Qc"}20CAT = {"erable": "Maple & honey", "epicerie": "Fine grocery", "cafe_the": "Coffee & tea", "chocolat": "Chocolate", "boulangerie": "Bakery", "viande_poisson": "Meat & fish", "fromage_laitier": "Cheese & dairy", "boissons": "Beverages", "alcool": "Alcohol", "sante_beaute": "Beauty & care", "maison": "Home & decor", "mode": "Fashion", "bijoux": "Jewellery", "art": "Art & craft", "enfants": "Children", "animaux": "Pets", "plein_air": "Outdoor & sport", "jardin": "Garden", "industriel": "Industrial", "autre": "Other"}212223def caption(name, text):24 (FIG / f"{name}.caption.txt").write_text(text)252627def save(fig, name):28 fig.tight_layout(); fig.savefig(FIG / f"{name}.pdf"); plt.close(fig); log.info("figure %s", name)293031def fig_raa_heatmap():32 m = pd.read_csv(PROC / "raa_matrix.csv"); m = m[(m["sample"] == "S_R") & (m.scheme == "M")]33 piv = m.pivot(index="region", columns="category", values="raa"); sig = m.pivot(index="region", columns="category", values="sig_gt1") | m.pivot(index="region", columns="category", values="sig_lt1")34 order_r = piv.index.tolist(); order_c = [c for c in CAT if c in piv.columns]35 piv = piv.loc[order_r, order_c]; sig = sig.loc[order_r, order_c]36 fig, ax = plt.subplots(figsize=(8.2, 5.6))37 im = ax.imshow(np.log2(piv.values), cmap=DIV, norm=TwoSlopeNorm(vcenter=0, vmin=-2.5, vmax=2.5), aspect="auto")38 ax.set_xticks(range(len(order_c))); ax.set_xticklabels([CAT[c] for c in order_c], rotation=55, ha="right"); ax.set_yticks(range(len(order_r))); ax.set_yticklabels([SHORT.get(r, r) for r in order_r])39 for i in range(len(order_r)):40 for j in range(len(order_c)):41 if sig.values[i, j]: ax.text(j, i, "•", ha="center", va="center", color="#222", fontsize=7)42 ax.grid(False); cb = fig.colorbar(im, ax=ax, fraction=0.03, pad=0.02); cb.set_label("log2 RAA (merchant-equal weights)")43 save(fig, "raa_heatmap"); caption("raa_heatmap", f"Revealed assortment advantage by region and category, merchant-equal weights, regional sample (1 214 stores, 216 861 listings). Dots mark cells whose store-bootstrap 95% interval excludes 1. {SRC}")444546def fig_rank_size():47 rs = pd.read_csv(PROC / "rank_size.csv"); h3 = json.load(open(EXP_DIR / "h3_longtail.json"))48 fig, ax = plt.subplots(figsize=(5.2, 3.8)); ax.loglog(rs["rank"], rs.n_listings, color=C[0], lw=1.6)49 pl = h3["powerlaw"]; xm = pl["xmin"]; ax.axhline(xm, color=GRAY, lw=0.8, ls="--"); ax.text(1.2, xm * 1.15, f"power-law tail from {xm} listings (α = {pl['alpha']:.2f})", color="#444", fontsize=8)50 ax.set_xlabel("Store rank (by catalogue size)"); ax.set_ylabel("Listings in the assortment sample")51 save(fig, "rank_size"); caption("rank_size", f"Rank–size distribution of merchant catalogue sizes, {h3['stores']} stores, {h3['listings']:,} listings. Gini {h3['gini']:.3f} (95% CI {h3['gini_ci'][0]:.3f}–{h3['gini_ci'][1]:.3f}); top 10% of stores hold {h3['top10pct_share']:.1%} of listings. {SRC}".replace(",", " "))525354def fig_diversity():55 d = pd.read_csv(PROC / "diversity_region.csv").dropna(subset=["H_rarefied15"]).sort_values("H_rarefied15")56 fig, ax = plt.subplots(figsize=(5.6, 4.6)); y = np.arange(len(d))57 ax.errorbar(d.H_rarefied15, y, xerr=[(d.H_rarefied15 - d.H_rarefied15_lo).clip(lower=0), (d.H_rarefied15_hi - d.H_rarefied15).clip(lower=0)], fmt="o", color=C[0], ms=5, capsize=2, lw=1, label="rarefied (15 stores), 95% range")58 ax.scatter(d.H_norm_M, y, marker="s", s=18, color=C[1], label="merchant-equal, all stores", zorder=3)59 ax.set_yticks(y); ax.set_yticklabels([SHORT.get(r, r) for r in d.region]); ax.set_xlabel("Normalised Shannon entropy of category shares"); ax.legend(loc="lower right", fontsize=8)60 save(fig, "diversity_region"); caption("diversity_region", f"Regional assortment diversity: rarefied entropy (mean over 2 000 draws of 15 stores; bars = 2.5–97.5 percentiles of the draws) and merchant-equal entropy on all stores; 16 regions with ≥ 15 stores (Nord-du-Québec, 3 stores, omitted). {SRC}")616263def fig_price_gaps():64 g = pd.read_parquet(PROC / "h5_groups.parquet"); gaps = g.groupby("group").lp.agg(lambda s: s.max() - s.min())65 pr = pd.read_parquet(P["interim"] / "match_pairs.parquet")66 fig, axes = plt.subplots(1, 2, figsize=(8, 3.4))67 axes[0].hist(gaps.clip(upper=1.0), bins=40, color=C[0], edgecolor="white", lw=0.4); axes[0].set_xlabel("Within-group max–min log price gap (matched groups, tiers A+B)"); axes[0].set_ylabel("Groups")68 axes[0].axvline(gaps.median(), color=C[1], lw=1.2); axes[0].text(gaps.median() + 0.02, axes[0].get_ylim()[1] * 0.9, f"median {gaps.median():.3f}", color=C[1], fontsize=8)69 for i, t in enumerate("ABCD"):70 s = pr[pr.tier == t].log_gap.clip(upper=1.0); axes[1].plot(np.sort(s), np.linspace(0, 1, len(s)), color=C[i] if i < 4 else GRAY, lw=1.4, label=f"tier {t} (n = {len(s)})")71 axes[1].set_xlabel("|log price gap| of candidate pairs"); axes[1].set_ylabel("Cumulative share"); axes[1].legend(fontsize=8, loc="lower right")72 save(fig, "price_gaps"); caption("price_gaps", f"Left: distribution of within-group price gaps for {len(gaps)} matched product groups (tiers A+B, comparable variants, both stores with a known region, {g.store.nunique()} stores). Right: cumulative distribution of absolute log price gaps for all candidate pairs by matching tier (tier D is never used for inference). Prices are posted prices at the snapshot. {SRC}")737475def fig_region_effects():76 e = pd.read_csv(PROC / "h5_region_effects.csv"); e = e[e.label == "primary_AB"].sort_values("effect")77 fig, ax = plt.subplots(figsize=(5.4, 4.2)); y = np.arange(len(e))78 ax.errorbar(e.effect, y, xerr=1.96 * e.se, fmt="o", color=C[0], ms=5, capsize=2, lw=1); ax.axvline(0, color=GRAY, lw=0.8)79 ax.set_yticks(y); ax.set_yticklabels([f"{SHORT.get(r, r)} ({int(n)})" for r, n in zip(e.region, e.n_stores)]); ax.set_xlabel("Region effect on log posted price (vs Montréal), ±1.96 cluster SE")80 save(fig, "h5_region_effects"); caption("h5_region_effects", f"Region fixed effects from the matched-product hedonic regression (group FE, platform FE, log catalogue size, reseller and certification controls; SE clustered by store; number of stores per region in parentheses). Point estimates are identified only through groups spanning regions and should be read with the wild-cluster-bootstrap joint test in the text. {SRC}")818283def fig_network():84 import networkx as nx85 h7 = json.load(open(EXP_DIR / "h7_complexity.json")); phi = pd.read_csv(PROC / "category_proximity.csv", index_col=0)86 G = nx.Graph()87 for a in phi.index:88 for b in phi.columns:89 if a < b and phi.loc[a, b] > 0: G.add_edge(a, b, weight=phi.loc[a, b])90 mst = nx.maximum_spanning_tree(G, weight="weight"); comm = h7["V4_community_stability"].get("communities_run0", {})91 pos = nx.kamada_kawai_layout(mst, weight="weight"); fig, ax = plt.subplots(figsize=(6.4, 5))92 strong = [(a, b) for a, b, w in G.edges(data="weight") if w >= np.percentile([w for *_, w in G.edges(data="weight")], 80)]93 nx.draw_networkx_edges(G, pos, edgelist=strong, ax=ax, edge_color="#d9d9d4", width=0.8); nx.draw_networkx_edges(mst, pos, ax=ax, edge_color=GRAY, width=1.6)94 cols = [C[comm.get(n, 0) % 5] for n in mst.nodes]; sizes = [120 + 10 * phi.loc[n].sum() for n in mst.nodes]95 nx.draw_networkx_nodes(mst, pos, ax=ax, node_color=cols, node_size=sizes, edgecolors="white", linewidths=1); nx.draw_networkx_labels(mst, pos, {n: CAT[n] for n in mst.nodes}, ax=ax, font_size=7.5)96 ax.axis("off"); ax.grid(False)97 save(fig, "category_network"); caption("category_network", f"Category relatedness: maximum spanning tree (dark edges) and top-quintile proximity edges (light) of the Hidalgo-type proximity computed from co-occurrence of categories across {int(pd.read_parquet(PROC / 'store_category.parquet').store_id.nunique())} stores; node colour = Louvain community (two communities, NMI = {h7['V4_community_stability'].get('nmi_median_vs_run0')} across 200 seeds). {SRC}")9899100def fig_attrition():101 flow = json.load(open(BASE_EXP / "sample_flow.json"))["flow"]102 fig, ax = plt.subplots(figsize=(7.4, 3.2)); labels = [k.replace("… ", "").replace(" [ASSORTMENT SAMPLE]", "").replace(" [PRICE SAMPLE]", "") for k, _ in flow]; vals = [v for _, v in flow]103 y = np.arange(len(vals))[::-1]; ax.barh(y, vals, color=[C[0] if "SAMPLE" in k else "#9ec5f4" for k, _ in flow], height=0.6)104 for yi, v, (k, _) in zip(y, vals, flow): ax.text(v + 3000, yi, f"{v:,}".replace(",", " "), va="center", fontsize=8)105 ax.set_yticks(y); ax.set_yticklabels([l[:70] for l in labels], fontsize=7.5); ax.set_xlabel("Listings"); ax.set_xlim(0, max(vals) * 1.18)106 save(fig, "attrition"); caption("attrition", f"Sample flow from database rows to the assortment and price samples (see Table on sample flow). {SRC}")107108109def fig_h7():110 h7 = json.load(open(EXP_DIR / "h7_complexity.json")); ind = pd.read_csv(PROC / "h7_region_indicators.csv")111 fig, axes = plt.subplots(1, 2, figsize=(8, 3.4))112 axes[0].scatter(ind.diversity, ind.eci, color=C[0], s=28)113 for _, r in ind.iterrows(): axes[0].annotate(SHORT.get(r.region, r.region), (r.diversity, r.eci), fontsize=6.5, xytext=(3, 2), textcoords="offset points")114 axes[0].set_xlabel("Diversity (categories with RAA ≥ 1)"); axes[0].set_ylabel("ECI (standardised)")115 v1 = h7["V1_bootstrap_rank_stability"]; axes[1].bar(["ECI", "Fitness", "Diversity"], [v1["eci_rho_median"], v1["fitness_rho_median"], v1["diversity_rho_median"]], color=[C[0], C[1], C[2]], width=0.55)116 axes[1].axhline(0.8, color=GRAY, ls="--", lw=0.8); axes[1].text(2.3, 0.82, "criterion 0.8", fontsize=7.5, color="#444", ha="right"); axes[1].set_ylabel("Median Spearman ρ, 500 store-bootstrap draws"); axes[1].set_ylim(0, 1)117 save(fig, "h7_validation"); caption("h7_validation", f"Left: ECI against diversity across the 17 regions (binary matrix M = 1[RAA ≥ 1], merchant-equal weights). Right: rank stability of the three indicators under store bootstrap; the preregistered criterion (median ρ ≥ 0.8) is not met by the ECI. {SRC}")118119120def fig_h4():121 h4 = json.load(open(EXP_DIR / "h4_dynamics.json")); d = pd.read_parquet(PROC / "h4_panel.parquet")122 d["bin"] = pd.qcut(d.density, 10, labels=False, duplicates="drop"); g = d.groupby("bin").agg(x=("density", "mean"), y=("entry", "mean"), n=("entry", "size"))123 fig, ax = plt.subplots(figsize=(5, 3.4)); ax.scatter(g.x, g.y, color=C[0], s=30); ax.set_xlabel("Relatedness density of the absent category (baseline weeks)"); ax.set_ylabel("Share of (store, category) cells with entry")124 ax.set_title(f"LPM slope {h4['lpm_density_coef']:.3f} (cluster SE {h4['se_cluster_store']:.3f}); placebo 95th pct {h4['placebo_p95']:.3f}", fontsize=8)125 save(fig, "h4_density"); caption("h4_density", f"Exploratory: entry of stores into categories absent in their baseline portfolio (ISO weeks 36–39) during weeks 40–44, binned by relatedness density; {h4['entries']} entries among {h4['store_category_cells_absent_at_baseline']:,} absent cells, {h4['stores']} stores. {SRC}".replace(",", " "))126127128def fig_robustness():129 rb = json.load(open(EXP_DIR / "robustness.json")); keys = [k for k in rb if isinstance(rb[k], dict) and "h1_M" in rb[k]]130 fig, ax = plt.subplots(figsize=(6.2, 3.8)); y = np.arange(len(keys))[::-1]131 for j, (sch, col, lab) in enumerate((("h1_M", C[0], "merchant-equal (M)"), ("h1_T", C[1], "trimmed listings (T)"))):132 v = [rb[k][sch]["V"] for k in keys]; nm = [rb[k][sch]["null_mean"] for k in keys]133 ax.scatter(v, y + (0.15 if j == 0 else -0.15), color=col, s=26, label=f"observed V, {lab}", zorder=3); ax.scatter(nm, y + (0.15 if j == 0 else -0.15), marker="|", color=col, s=60, label=f"permutation-null mean, {lab}")134 ax.set_yticks(y); ax.set_yticklabels([k.replace("_", " ") for k in keys], fontsize=8); ax.set_xlabel("Cramér's V (region × category)"); ax.legend(fontsize=7, loc="lower right")135 save(fig, "robustness_h1"); caption("robustness_h1", f"H1 across preregistered sample and category variants: observed Cramér's V against the mean of the store-reshuffle permutation null, merchant-equal and trimmed weighting. {SRC}")136137138def fig_map():139 import glob140 try:141 import geopandas as gpd142 except ImportError:143 return144 cands = glob.glob(str(ROOT / "data/external/**/*.gpkg"), recursive=True) + glob.glob(str(ROOT / "data/external/**/*.shp"), recursive=True) + glob.glob(str(ROOT / "data/external/**/*.geojson"), recursive=True)145 layer = None146 for c in cands:147 try: g = gpd.read_file(c)148 except Exception: continue149 for col in g.columns:150 if g[col].dtype == object and sum(r in set(map(str, g[col].dropna())) or r.replace("–", "-") in set(map(str, g[col].dropna())) for r in REGIONS) >= 15:151 g["geometry"] = g.geometry.make_valid(); layer = g.to_crs(32198).dissolve(by=col).reset_index().rename(columns={col: "region"}); break152 if layer is not None: break153 if layer is None:154 log.warning("no region layer: map skipped"); return155 layer["region"] = layer.region.map(lambda v: next((r for r in REGIONS if r == v or r.replace("–", "-") == v), v)); layer["geometry"] = layer.geometry.simplify(1500)156 sp = pd.read_parquet(PROC / "store_panel.parquet"); cnt = sp[sp.region_final != "Inconnue"].groupby("region_final").agg(stores=("id", "size"), listings=("n_listings", "sum"))157 layer = layer.merge(cnt, left_on="region", right_index=True, how="left")158 pts = sp[sp.lat.notna()]; gp = gpd.GeoDataFrame(pts, geometry=gpd.points_from_xy(pts.lng, pts.lat), crs=4326).to_crs(32198)159 fig, axes = plt.subplots(1, 2, figsize=(9, 3.6))160 layer.plot(column="stores", cmap=SEQ, ax=axes[0], edgecolor="white", linewidth=0.5, legend=True, legend_kwds={"shrink": 0.5, "label": "stores with known region"})161 layer.plot(color="#f2f2ef", edgecolor="white", linewidth=0.5, ax=axes[1]); gp.plot(ax=axes[1], markersize=np.sqrt(gp.n_listings.clip(1, 5000)) * 0.6, color=C[1], alpha=0.5, linewidth=0)162 xs, ys = gp.geometry.x, gp.geometry.y; x0, x1 = np.percentile(xs, [0.5, 99.5]); y0, y1 = np.percentile(ys, [0.5, 99.5])163 for ax in axes: ax.set_axis_off(); ax.set_xlim(x0 - 6e4, x1 + 6e4); ax.set_ylim(y0 - 6e4, y1 + 6e4); ax.set_aspect("equal")164 axes[0].set_title("Stores per administrative region", fontsize=9); axes[1].set_title("Store locations (FSA centroids), size ∝ √listings", fontsize=9)165 save(fig, "coverage_map"); caption("coverage_map", f"Coverage of the regional sample: {int(cnt.stores.sum())} stores with a known region (left) and FSA-centroid locations of {len(gp)} geocoded stores (right); the southern part of the province is shown. Official administrative-region boundaries. {SRC}")166167168def main():169 for f in (fig_attrition, fig_raa_heatmap, fig_rank_size, fig_diversity, fig_price_gaps, fig_region_effects, fig_network, fig_h7, fig_h4, fig_robustness, fig_map):170 try: f()171 except Exception as e: # noqa: BLE001172 log.error("%s failed: %s", f.__name__, e)173 ledger("reports.make_figures", figures=sorted(p.name for p in FIG.glob("*.pdf")))174175176if __name__ == "__main__":177 main()178