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 # all PUBLISHED listings of live stores, including listings inactive at the snapshot (listings_clean keeps active=1 only)48 raw = pd.read_parquet(P["interim"] / "listings_raw.parquet"); raw = raw[raw.listing_status == "published"]49 lc = pd.read_parquet(P["interim"] / "listings_clean.parquet")50 live = set(lc.loc[lc.f_live_store, "store_id"]); bad = set(lc.loc[lc.f_hors_mission | lc.f_sibling_dup, "store_id"])51 li = raw[raw.store_id.isin(live) & ~raw.store_id.isin(bad)].copy()52 if LABEL == "model":53 cm = pd.read_parquet(BASE_PROC / "category_model.parquet").set_index("uid").category_model54 li = li.assign(category=li.uid.map(cm).fillna(li.category))55 sp = pd.read_parquet(PROC / "store_panel.parquet").set_index("id")56 base = present(li, *BASE); out = present(li, *OUT)57 phi, pv = proximity(base)58 cats = list(pv.columns); stores = pv.index59 dens = (pv.values @ phi.values) / np.maximum(phi.values.sum(0), 1e-9) # density of c around store's baseline portfolio60 D = pd.DataFrame(dens, index=stores, columns=cats).stack().rename("density").reset_index().rename(columns={"level_0": "store_id", "level_1": "category"})61 D = D.merge(pv.stack().rename("base_present").reset_index().rename(columns={"level_0": "store_id", "level_1": "category"}), on=["store_id", "category"])62 outp = out.assign(out_present=1)[["store_id", "category", "out_present"]]63 D = D.merge(outp, on=["store_id", "category"], how="left"); D["out_present"] = D.out_present.fillna(0)64 D = D[D.base_present == 0].copy(); D["entry"] = D.out_present.astype(int)65 b, se = lpm(D, rng)66 # placebo: relabel categories in the proximity matrix (random permutation of rows+columns of phi) so that the67 # store's portfolio is scored against a wrong relatedness structure; 200 draws68 plac = []69 for k in range(200):70 pi = rng.permutation(len(cats)); php = phi.values[np.ix_(pi, pi)]71 dp = (pv.values @ php) / np.maximum(php.sum(0), 1e-9)72 Dp = D.copy(); Dp["density"] = pd.DataFrame(dp, index=stores, columns=cats).stack().reindex(pd.MultiIndex.from_frame(D[["store_id", "category"]])).values73 plac.append(lpm(Dp, rng)[0])74 # region-level entries (descriptive)75 sr = sp.region_final76 rb = base.assign(region=base.store_id.map(sr)).query("region!='Inconnue'").groupby(["region", "category"]).n.sum()77 ro = out.assign(region=out.store_id.map(sr)).query("region!='Inconnue'").groupby(["region", "category"]).n.sum()78 all_rc = pd.MultiIndex.from_product([sorted(sr[sr != "Inconnue"].unique()), cats])79 rb = rb.reindex(all_rc).fillna(0); ro = ro.reindex(all_rc).fillna(0)80 region_entries = [(r, c) for (r, c) in all_rc if rb[(r, c)] == 0 and ro[(r, c)] > 0]81 region_exits = [(r, c) for (r, c) in all_rc if rb[(r, c)] > 0 and ro[(r, c)] == 0]82 # connector audit: entering stores with an 'added' burst (> 50 % of baseline catalogue added in one sync) in the outcome window83 sl = pd.read_parquet(P["interim"] / "sync_log.parquet"); slo = sl[(sl.ts_utc >= OUT[0]) & (sl.ts_utc < OUT[1])]84 basesize = base.groupby("store_id").n.sum()85 burst = slo.groupby("store_id").added.max() / basesize.reindex(slo.store_id.unique()).fillna(1)86 entering = D[D.entry == 1].store_id.unique()87 share_burst = float((burst.reindex(entering).fillna(0) > 0.5).mean()) if len(entering) else None88 inactive_n = int((li.active == 0).sum())89 res = {"status": "exploratory", "listings_used": int(len(li)), "inactive_listings_included": inactive_n, "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())],90 "stores": int(pv.shape[0]), "store_category_cells_absent_at_baseline": int(len(D)), "entries": int(D.entry.sum()), "entering_stores": int(len(entering)),91 "lpm_density_coef": round(b, 4), "se_cluster_store": round(se, 4), "t": round(b / se, 3) if se > 0 else None,92 "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)),93 "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),94 "region_level_entries": region_entries, "region_level_exits": region_exits, "n_region_entries": len(region_entries),95 "share_entering_stores_with_added_burst": share_burst}96 D.to_parquet(PROC / "h4_panel.parquet", index=False); register_output(PROC / "h4_panel.parquet", "h4")97 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})98 log.info(json.dumps({k: v for k, v in res.items() if "region_level" not in k}))99100101if __name__ == "__main__":102 main()103