1"""H4 — Related diversification (EXPLORATORY). Weekly (store, category) presence reconstructed from first_seen/last_seen.2Baseline weeks W36–W39 (2026-08-31 … 2026-09-27): portfolio and relatedness (proximity across stores) computed here only.3Outcome: entry into a category absent in W36–W39 and present (>= 1 listing) in W40–W44 (2026-09-28 … 2026-11-01).4LPM Entry_sc = b * Density_sc + a_s + g_c + e, SE clustered by store; placebo with category-relabelled proximity.5Region-level entries counted descriptively. Connector-change audit: share of entering stores whose sync_log shows a6status change or an 'added' burst in the outcome window.7"""8from __future__ import annotations9import json10import numpy as np, pandas as pd11from src.common import ROOT, P, SEED, get_logger, ledger, register_output, write_json, PROC_DIR, EXP_DIR, BASE_PROC, BASE_EXP, LABEL1213log = get_logger("models.h4")14PROC = PROC_DIR15BASE = (pd.Timestamp("2026-08-31", tz="UTC"), pd.Timestamp("2026-09-28", tz="UTC"))16OUT = (pd.Timestamp("2026-09-28", tz="UTC"), pd.Timestamp("2026-11-02", tz="UTC"))171819def present(li, start, end):20 """listings alive at some point in [start, end): first_seen < end and last_seen >= start"""21 m = (li.first_seen_utc < end) & (li.last_seen_utc >= start)22 return li[m].groupby(["store_id", "category"]).size().rename("n").reset_index()232425def proximity(pres):26 pv = pres.assign(v=1).pivot_table(index="store_id", columns="category", values="v", fill_value=0)27 co = pv.T.values @ pv.values; k = np.diag(co).astype(float)28 phi = co / np.maximum(np.minimum.outer(k, k), 1); np.fill_diagonal(phi, 0); return pd.DataFrame(phi, index=pv.columns, columns=pv.columns), pv293031def lpm(df, rng):32 # within transformation for two-way FE (store, category) via iterative demeaning33 y = df.entry.values.astype(float); x = df.density.values.astype(float)34 s = df.store_id.astype("category").cat.codes.values; c = df.category.astype("category").cat.codes.values35 def demean(v):36 v = v.copy()37 for _ in range(50):38 v -= pd.Series(v).groupby(s).transform("mean").values; v -= pd.Series(v).groupby(c).transform("mean").values39 return v40 yw, xw = demean(y), demean(x); b = float((xw * yw).sum() / (xw ** 2).sum()); r = yw - b * xw41 g = pd.Series(xw * r).groupby(s).sum(); se = float(np.sqrt((g ** 2).sum()) / (xw ** 2).sum())42 return b, se434445def main():46 rng = np.random.default_rng(SEED)47 li = pd.read_parquet(P["interim"] / "listings_clean.parquet"); li = li[li.f_published | (li.active == 0)] # all live-store listings incl. inactive48 li = li[li.f_live_store & ~li.f_hors_mission & ~li.f_sibling_dup]49 if LABEL == "model":50 cm = pd.read_parquet(BASE_PROC / "category_model.parquet").set_index("uid").category_model51 li = li.assign(category=li.uid.map(cm).fillna(li.category))52 sp = pd.read_parquet(PROC / "store_panel.parquet").set_index("id")53 base = present(li, *BASE); out = present(li, *OUT)54 phi, pv = proximity(base)55 cats = list(pv.columns); stores = pv.index56 dens = (pv.values @ phi.values) / np.maximum(phi.values.sum(0), 1e-9) # density of c around store's baseline portfolio57 D = pd.DataFrame(dens, index=stores, columns=cats).stack().rename("density").reset_index().rename(columns={"level_0": "store_id", "level_1": "category"})58 D = D.merge(pv.stack().rename("base_present").reset_index().rename(columns={"level_0": "store_id", "level_1": "category"}), on=["store_id", "category"])59 outp = out.assign(out_present=1)[["store_id", "category", "out_present"]]60 D = D.merge(outp, on=["store_id", "category"], how="left"); D["out_present"] = D.out_present.fillna(0)61 D = D[D.base_present == 0].copy(); D["entry"] = D.out_present.astype(int)62 b, se = lpm(D, rng)63 # placebo: relabel categories in the proximity matrix (random permutation of rows+columns of phi) so that the64 # store's portfolio is scored against a wrong relatedness structure; 200 draws65 plac = []66 for k in range(200):67 pi = rng.permutation(len(cats)); php = phi.values[np.ix_(pi, pi)]68 dp = (pv.values @ php) / np.maximum(php.sum(0), 1e-9)69 Dp = D.copy(); Dp["density"] = pd.DataFrame(dp, index=stores, columns=cats).stack().reindex(pd.MultiIndex.from_frame(D[["store_id", "category"]])).values70 plac.append(lpm(Dp, rng)[0])71 # region-level entries (descriptive)72 sr = sp.region_final73 rb = base.assign(region=base.store_id.map(sr)).query("region!='Inconnue'").groupby(["region", "category"]).n.sum()74 ro = out.assign(region=out.store_id.map(sr)).query("region!='Inconnue'").groupby(["region", "category"]).n.sum()75 all_rc = pd.MultiIndex.from_product([sorted(sr[sr != "Inconnue"].unique()), cats])76 rb = rb.reindex(all_rc).fillna(0); ro = ro.reindex(all_rc).fillna(0)77 region_entries = [(r, c) for (r, c) in all_rc if rb[(r, c)] == 0 and ro[(r, c)] > 0]78 region_exits = [(r, c) for (r, c) in all_rc if rb[(r, c)] > 0 and ro[(r, c)] == 0]79 # connector audit: entering stores with an 'added' burst (> 50 % of baseline catalogue added in one sync) in the outcome window80 sl = pd.read_parquet(P["interim"] / "sync_log.parquet"); slo = sl[(sl.ts_utc >= OUT[0]) & (sl.ts_utc < OUT[1])]81 basesize = base.groupby("store_id").n.sum()82 burst = slo.groupby("store_id").added.max() / basesize.reindex(slo.store_id.unique()).fillna(1)83 entering = D[D.entry == 1].store_id.unique()84 share_burst = float((burst.reindex(entering).fillna(0) > 0.5).mean()) if len(entering) else None85 res = {"status": "exploratory", "baseline_window": [str(BASE[0].date()), str((BASE[1] - pd.Timedelta(days=1)).date())], "outcome_window": [str(OUT[0].date()), str((OUT[1] - pd.Timedelta(days=1)).date())],86 "stores": int(pv.shape[0]), "store_category_cells_absent_at_baseline": int(len(D)), "entries": int(D.entry.sum()), "entering_stores": int(len(entering)),87 "lpm_density_coef": round(b, 4), "se_cluster_store": round(se, 4), "t": round(b / se, 3) if se > 0 else None,88 "placebo_coef_mean": round(float(np.mean(plac)), 4), "placebo_p95": round(float(np.percentile(plac, 95)), 4), "placebo_p": float((np.sum(np.array(plac) >= b) + 1) / (len(plac) + 1)),89 "density_mean_entry": round(float(D[D.entry == 1].density.mean()), 4), "density_mean_no_entry": round(float(D[D.entry == 0].density.mean()), 4),90 "region_level_entries": region_entries, "region_level_exits": region_exits, "n_region_entries": len(region_entries),91 "share_entering_stores_with_added_burst": share_burst}92 D.to_parquet(PROC / "h4_panel.parquet", index=False); register_output(PROC / "h4_panel.parquet", "h4")93 write_json(res, EXP_DIR / "h4_dynamics.json", "h4"); ledger("models.h4_dynamics", **{k: v for k, v in res.items() if "region_level" not in k})94 log.info(json.dumps({k: v for k, v in res.items() if "region_level" not in k}))959697if __name__ == "__main__":98 main()99