1"""H5 — Cross-merchant price dispersion of matched products (cross-section at the snapshot).23Sample S_M: match groups (tiers A+B primary; A+B+C sensitivity) with comparable variant structure, both stores in S_R4(region known), stores not siblings, groups spanning >= 2 stores after filters.5Model: log P_im = alpha_i + gamma_r(m) + delta_platform(m) + b1 log n_m + b2 Reseller_m + b3 Certified_m + e_im,6group fixed effects absorbed by within-group demeaning; SE clustered by store; joint test of region effects and the7t-test on Reseller use a wild cluster bootstrap (Rademacher, 9 999 draws, restricted residuals).8Outputs experiments/h5_hedonic.json, data/processed/h5_groups.parquet, data/processed/h5_region_effects.csv9"""10from __future__ import annotations11import json12import numpy as np, pandas as pd13from src.common import ROOT, P, SEED, get_logger, ledger, register_output, write_json, PROC_DIR, EXP_DIR1415log = get_logger("models.h5")16PROC = PROC_DIR17B = 9999181920def build_sample(tiers):21 pr = pd.read_parquet(P["interim"] / "match_pairs.parquet")22 li = pd.read_parquet(P["interim"] / "listings_clean.parquet").set_index("uid")23 sp = pd.read_parquet(PROC / "store_panel.parquet").set_index("id")24 sib = pd.read_parquet(P["interim"] / "sibling_stores.parquet").set_index("id").sibling_rep25 pr = pr[pr.tier.isin(tiers) & pr.comparable]26 ok = lambda s: s.map(sp.region_final).fillna("Inconnue") != "Inconnue"27 pr = pr[ok(pr.store_a) & ok(pr.store_b)]28 pr = pr[pr.store_a.map(sib).fillna(pr.store_a) != pr.store_b.map(sib).fillna(pr.store_b)]29 import networkx as nx30 G = nx.Graph(); G.add_edges_from(pr[["uid_a", "uid_b"]].itertuples(index=False, name=None))31 rows = [(g, u) for g, comp in enumerate(nx.connected_components(G)) for u in comp]32 d = pd.DataFrame(rows, columns=["group", "uid"])33 d["store"] = d.uid.map(li.store_id); d["price"] = d.uid.map(li.price); d["category"] = d.uid.map(li.category)34 d["region"] = d.store.map(sp.region_final); d["platform"] = d.store.map(sp.platform).fillna("other")35 d["log_n"] = np.log(d.store.map(sp.n_listings)); d["reseller"] = d.store.map(sp.store_kind).isin(["revendeur", "collectif"]).astype(float)36 d["certified"] = d.store.map(sp.certified).fillna(False).astype(float); d["lp"] = np.log(d.price)37 gs = d.groupby("group").store.nunique(); d = d[d.group.map(gs) >= 2]38 return d394041def within(df, cols):42 return df[cols] - df.groupby("group")[cols].transform("mean")434445def fit(d, rng, label):46 # design: region dummies (drop the largest region), platform dummies (drop shopify), controls47 regions = sorted(d.region.unique()); base_r = d.region.value_counts().idxmax()48 plats = sorted(d.platform.unique()); base_p = "shopify" if "shopify" in plats else plats[0]49 X = pd.DataFrame(index=d.index)50 for r in regions:51 if r != base_r: X[f"r::{r}"] = (d.region == r).astype(float)52 for p_ in plats:53 if p_ != base_p: X[f"p::{p_}"] = (d.platform == p_).astype(float)54 X["log_n"] = d.log_n; X["reseller"] = d.reseller; X["certified"] = d.certified55 cols = list(X.columns); dd = pd.concat([d[["group", "store", "lp"]], X], axis=1)56 W = within(dd, cols + ["lp"]); Xw = W[cols].values; yw = W["lp"].values57 keep = Xw.std(0) > 1e-10; Xw = Xw[:, keep]; cols = [c for c, k in zip(cols, keep) if k]58 beta, *_ = np.linalg.lstsq(Xw, yw, rcond=None); e = yw - Xw @ beta59 cl = d.store.astype("category").cat.codes.values; G = cl.max() + 160 XtX_inv = np.linalg.pinv(Xw.T @ Xw)61 def cluster_cov(Xw, e):62 meat = np.zeros((Xw.shape[1], Xw.shape[1]))63 for g in range(G):64 m = cl == g; s = Xw[m].T @ e[m]; meat += np.outer(s, s)65 return XtX_inv @ meat @ XtX_inv * G / (G - 1)66 V = cluster_cov(Xw, e); se = np.sqrt(np.diag(V))67 ridx = [i for i, c in enumerate(cols) if c.startswith("r::")]68 def wald(beta, V, idx):69 b = beta[idx]; Vi = V[np.ix_(idx, idx)]; return float(b @ np.linalg.pinv(Vi) @ b)70 W_obs = wald(beta, V, ridx)71 # wild cluster bootstrap: restricted model (no region dummies) residuals, Rademacher by store72 nr = [i for i in range(len(cols)) if i not in ridx]73 b_r, *_ = np.linalg.lstsq(Xw[:, nr], yw, rcond=None); e_r = yw - Xw[:, nr] @ b_r; yhat_r = Xw[:, nr] @ b_r74 Wb = np.empty(B)75 i_res = cols.index("reseller") if "reseller" in cols else None76 nr2 = [i for i in range(len(cols)) if i != i_res] if i_res is not None else None77 if i_res is not None:78 b_r2, *_ = np.linalg.lstsq(Xw[:, nr2], yw, rcond=None); e_r2 = yw - Xw[:, nr2] @ b_r2; yhat_r2 = Xw[:, nr2] @ b_r279 tb = np.empty(B); t_obs = beta[i_res] / se[i_res] if i_res is not None else np.nan80 gcode = d.group.astype("category").cat.codes.values; ng = gcode.max() + 1; gsize = np.bincount(gcode)81 def gdemean(v): # the match-group fixed effects are absorbed by within-transformation: the bootstrap disturbance82 return v - (np.bincount(gcode, weights=v) / gsize)[gcode] # must be re-projected off the group means (referee 1, must-fix)83 for b in range(B):84 w = rng.choice([-1.0, 1.0], G)[cl]85 yb = yhat_r + gdemean(e_r * w); bb = XtX_inv @ (Xw.T @ yb); eb = yb - Xw @ bb86 Wb[b] = wald(bb, cluster_cov(Xw, eb), ridx)87 if i_res is not None:88 yb2 = yhat_r2 + gdemean(e_r2 * w); bb2 = XtX_inv @ (Xw.T @ yb2); eb2 = yb2 - Xw @ bb2; V2 = cluster_cov(Xw, eb2)89 tb[b] = bb2[i_res] / np.sqrt(V2[i_res, i_res])90 p_region = float((np.sum(Wb >= W_obs) + 1) / (B + 1))91 p_res = float((np.sum(np.abs(tb) >= abs(t_obs)) + 1) / (B + 1)) if i_res is not None else None92 # leave-one-region-out joint Wald (which region drives the rejection?) and single-region Wald contributions93 loo = {}94 for r in regions:95 if r == base_r: continue96 keep_idx = [i for i in ridx if cols[i] != f"r::{r}"]97 loo[r] = round(wald(beta, V, keep_idx), 2)98 contrib = {cols[i][3:]: round(float(beta[i] ** 2 / V[i, i]), 2) for i in ridx}99 n_groups = int(d.group.nunique())100 gaps = d.groupby("group").lp.agg(lambda s: s.max() - s.min())101 reg_stores = d.groupby("region").store.nunique()102 res = {"label": label, "obs": int(len(d)), "groups": n_groups, "stores": int(G), "regions": len(regions), "base_region": base_r,103 "regions_with_ge5_stores": int((reg_stores >= 5).sum()), "stores_by_region": reg_stores.to_dict(),104 "coef": {c: {"b": round(float(b), 4), "se_cluster": round(float(s), 4)} for c, b, s in zip(cols, beta, se)},105 "wald_region": round(W_obs, 3), "df_region": len(ridx), "p_region_wildboot": p_region, "t_reseller": round(float(t_obs), 3) if i_res is not None else None, "p_reseller_wildboot": p_res,106 "within_R2": round(float(1 - e.var() / yw.var()), 4), "wald_loo_region": loo, "wald_single_region_contrib": contrib, "groups_spanning_regions": int((d.groupby("group").region.nunique() >= 2).sum()),107 "gap_median": round(float(gaps.median()), 4), "gap_p90": round(float(gaps.quantile(0.9)), 4), "share_gap_gt5pct": round(float((gaps > 0.05).mean()), 4),108 "share_gap_zero": round(float((gaps < 1e-9).mean()), 4), "B": B}109 eff = pd.DataFrame({"region": [c[3:] for c in cols if c.startswith("r::")], "effect": [beta[i] for i in ridx], "se": [se[i] for i in ridx]})110 eff = pd.concat([eff, pd.DataFrame({"region": [base_r], "effect": [0.0], "se": [0.0]})]); eff["label"] = label; eff["n_stores"] = eff.region.map(reg_stores)111 return res, eff, d112113114def main():115 rng = np.random.default_rng(SEED); out = []; effs = []116 for label, tiers in (("primary_AB", ["A", "B"]), ("sens_ABC", ["A", "B", "C"])):117 d = build_sample(tiers)118 if d.group.nunique() < 20:119 out.append({"label": label, "status": "insufficient", "groups": int(d.group.nunique())}); continue120 res, eff, dd = fit(d, rng, label); out.append(res); effs.append(eff); log.info(json.dumps({k: v for k, v in res.items() if k not in ("coef", "stores_by_region")}))121 if label == "primary_AB":122 dd.to_parquet(PROC / "h5_groups.parquet", index=False); register_output(PROC / "h5_groups.parquet", "h5")123 pd.concat(effs).to_csv(PROC / "h5_region_effects.csv", index=False)124 # descriptive certification contrast (NOT H6): within-category, store-clustered125 li = pd.read_parquet(P["interim"] / "listings_clean.parquet"); sp = pd.read_parquet(PROC / "store_panel.parquet").set_index("id")126 s = li[li.f_analysis_price].copy(); s["cert"] = s.store_id.map(sp.cert_adq | sp.cert_lpq).fillna(False).astype(float)127 s["region"] = s.store_id.map(sp.region_final); s = s[s.region.notna() & (s.region != "Inconnue")]128 s["y"] = s.log_price - s.groupby("category").log_price.transform("mean"); s["x"] = s.cert - s.groupby("category").cert.transform("mean")129 b = float((s.x * s.y).sum() / (s.x ** 2).sum()); r = s.y - b * s.x130 g = s.assign(r=r).groupby("store_id").apply(lambda h: (h.x * h.r).sum()); se = float(np.sqrt((g ** 2).sum()) / (s.x ** 2).sum())131 cert = {"note": "descriptive store-certification contrast (Aliments du Québec / Les Produits du Québec), category FE, store-clustered SE — NOT a test of H6",132 "obs": int(len(s)), "stores": int(s.store_id.nunique()), "certified_stores": int(s.groupby('store_id').cert.first().sum()), "b_log_price": round(b, 4), "se": round(se, 4)}133 write_json({"models": out, "certification_contrast": cert}, EXP_DIR / "h5_hedonic.json", "h5"); ledger("models.h5_hedonic", models=[{k: v for k, v in m.items() if k not in ("coef", "stores_by_region")} for m in out], cert=cert)134 log.info(json.dumps(cert))135136137if __name__ == "__main__":138 main()139