1"""Model A — matched-product cross-banner price dispersion (H1). Preregistration §2 H1.23Outputs: data/processed/modelA_*.csv and experiments/modelA_summary.json4"""5import json, re, numpy as np, pandas as pd, pyfixest as pf6from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED7log = get_logger("modelA"); P = path("processed"); E = ROOT / "experiments"; E.mkdir(exist_ok=True)89def group_week_stats(pm):10 g = pm.groupby(["group_id", "week"])11 st = g.agg(n_banners=("banner", "nunique"), n_listings=("uid", "size"), gap=("log_p", lambda s: s.max() - s.min()),12 cv=("p_mode", lambda s: s.std(ddof=0) / s.mean()), iqr=("log_p", lambda s: s.quantile(.75) - s.quantile(.25)),13 same_company=("company", lambda s: s.nunique() == 1), category=("category", "first"), ambiguous=("ambiguous", "max")).reset_index()14 return st1516def pair_gaps(pm):17 """Banner-pair log gaps per group-week (one listing per banner-channel; ambiguous group-weeks dropped upstream)."""18 rows = []19 for (gid, w), g in pm.groupby(["group_id", "week"]):20 g = g.drop_duplicates("banner")21 b = g.set_index("banner")22 bs = sorted(b.index)23 for i in range(len(bs)):24 for j in range(i + 1, len(bs)):25 rows.append((gid, w, bs[i], bs[j], b.at[bs[i], "log_p"] - b.at[bs[j], "log_p"], b.at[bs[i], "company"] == b.at[bs[j], "company"], g.category.iloc[0]))26 return pd.DataFrame(rows, columns=["group_id", "week", "banner_r", "banner_s", "gap", "same_company", "category"])2728def persistence(pg):29 pg = pg.sort_values(["group_id", "banner_r", "banner_s", "week"])30 pg["gap_lag"] = pg.groupby(["group_id", "banner_r", "banner_s"]).gap.shift(1)31 pg["week_lag"] = pg.groupby(["group_id", "banner_r", "banner_s"]).week.shift(1)32 ok = pg[(pg.week - pg.week_lag).dt.days == 7].dropna(subset=["gap_lag"])33 fit = pf.feols("gap ~ gap_lag", data=ok, vcov={"CRV1": "group_id"})34 rho = float(fit.coef()["gap_lag"]); se = float(fit.se()["gap_lag"])35 sign = pg.dropna(subset=["gap"]).groupby(["group_id", "banner_r", "banner_s"]).gap.agg(lambda s: (np.sign(s).nunique() == 1) and (s.abs() > 1e-9).all() and len(s) >= 4)36 nflip = pg.groupby(["group_id", "banner_r", "banner_s"]).gap.agg(lambda s: ((np.sign(s) != np.sign(s.shift(1))) & s.shift(1).notna()).sum())37 nweeks = pg.groupby(["group_id", "banner_r", "banner_s"]).size()38 stable = sign[nweeks >= 4].mean()39 return rho, se, float(stable), int((nweeks >= 4).sum()), ok4041def main():42 seed_everything()43 out = {}44 pm_all = pd.read_parquet(P / "panel_matched_AB.parquet")45 pm = pm_all[(pm_all.channel == "online_storefront") & (~pm_all.ambiguous)].copy()46 nb = pm.groupby(["group_id", "week"]).banner.transform("nunique"); pm = pm[nb >= 2]47 log.info("H1 primary sample: %d listing-weeks, %d groups, %d group-weeks", len(pm), pm.group_id.nunique(), pm.groupby(["group_id", "week"]).ngroups)48 st = group_week_stats(pm); st.to_csv(P / "modelA_groupweek_stats.csv", index=False)49 def q(s): return dict(n=int(len(s)), median=float(s.median()), p25=float(s.quantile(.25)), p75=float(s.quantile(.75)), mean=float(s.mean()), share_ge_5pct=float((s >= 0.05).mean()), share_ge_10pct=float((s >= 0.10).mean()))50 out["gap_all"] = q(st.gap); out["gap_same_company"] = q(st[st.same_company].gap); out["gap_cross_company"] = q(st[~st.same_company].gap)51 out["cv_all"] = q(st.cv)52 bycat = st.groupby("category").gap.agg(["size", "median", "mean"]).reset_index().sort_values("median", ascending=False); bycat.to_csv(P / "modelA_gap_by_category.csv", index=False)53 # Banner fixed effects: log_p ~ 1 | group_week + banner (reference banner = most frequent)54 pm["gw"] = pm.group_id + "_" + pm.week.dt.strftime("%Y%m%d")55 ref = pm.banner.value_counts().index[0]56 pm["banner_f"] = pd.Categorical(pm.banner, categories=[ref] + sorted(set(pm.banner) - {ref}))57 fit = pf.feols("log_p ~ i(banner_f, ref=" + repr(ref) + ") | gw", data=pm, vcov={"CRV1": "group_id"})58 co = fit.coef(); se = fit.se()59 fe = pd.DataFrame({"term": co.index, "coef": co.values, "se": se.values}); fe["banner"] = fe.term.map(lambda t: re.sub(r"^.*?::|^.*\[T\.|\]$", "", t))60 fe["n_listing_weeks"] = fe.banner.map(pm.banner.value_counts()); fe = fe.sort_values("coef"); fe.to_csv(P / "modelA_banner_effects.csv", index=False)61 out["banner_fe_reference"] = ref; out["resid_sd"] = float(np.sqrt(np.mean(fit.resid() ** 2))); out["n_obs_fe"] = int(fit._N)62 # two-way cluster sensitivity (group + banner)63 try:64 fit2 = pf.feols("log_p ~ i(banner_f, ref=" + repr(ref) + ") | gw", data=pm, vcov={"CRV1": "group_id + banner"})65 fe["se_twoway"] = fe.term.map(fit2.se()); fe.to_csv(P / "modelA_banner_effects.csv", index=False)66 except Exception as e:67 log.warning("two-way clustering failed: %s", e)68 # Persistence69 pg = pair_gaps(pm); pg.to_parquet(P / "modelA_pair_gaps.parquet", index=False)70 rho, se_rho, stable, npairs, ok = persistence(pg)71 out["ar1_rho"] = rho; out["ar1_se"] = se_rho; out["share_pairs_sign_stable_ge4w"] = stable; out["n_pairs_ge4w"] = npairs72 rho_same, se_same, st_same, n_same, _ = persistence(pg[pg.same_company]); rho_x, se_x, st_x, n_x, _ = persistence(pg[~pg.same_company])73 out["ar1_same_company"] = [rho_same, se_same, st_same, n_same]; out["ar1_cross_company"] = [rho_x, se_x, st_x, n_x]74 # Placebo: permute banner labels within group-week, recompute AR(1)75 rng = np.random.default_rng(SEED); rhos = []76 for b in range(20):77 pp = pm.copy(); pp["banner"] = pp.groupby(["group_id", "week"]).banner.transform(lambda s: rng.permutation(s.values))78 pgp = pair_gaps(pp); r, _, _, _, _ = persistence(pgp); rhos.append(r)79 out["placebo_ar1_mean"] = float(np.mean(rhos)); out["placebo_ar1_sd"] = float(np.std(rhos))80 # Banner-pair table (median gap and persistence) for the heatmap81 pairtab = pg.groupby(["banner_r", "banner_s"]).agg(n=("gap", "size"), median_gap=("gap", "median"), mean_gap=("gap", "mean"), share_r_cheaper=("gap", lambda s: (s < 0).mean())).reset_index()82 pairtab.to_csv(P / "modelA_banner_pairs.csv", index=False)83 # Within-banner store-mixture bound (storefront, grocery universe)84 pw = pd.read_parquet(P / "panel_weekly.parquet", columns=["banner", "channel", "n_distinct", "rel_range", "p_min", "p_max", "p_mode"])85 sf = pw[pw.channel == "online_storefront"]86 sm = sf.groupby("banner").agg(listing_weeks=("rel_range", "size"), share_multi=("n_distinct", lambda s: (s > 1).mean()),87 median_range_when_multi=("rel_range", lambda s: s[s > 0].median()), median_log_range_multi=("p_max", lambda s: np.nan)).reset_index()88 sf2 = sf[sf.n_distinct > 1].assign(lr=lambda d: np.log(d.p_max / d.p_min)); sm["median_log_range_multi"] = sm.banner.map(sf2.groupby("banner").lr.median())89 sm.sort_values("share_multi", ascending=False).to_csv(P / "modelA_within_banner_mixture.csv", index=False)90 # Sensitivity: tiers and price definitions91 sens = []92 for tiers in ["A", "AB", "ABC", "ABCD"]:93 for pdef in ["p_mode", "p_min", "p_max", "p_twmean"]:94 x = pd.read_parquet(P / f"panel_matched_{tiers}.parquet"); x = x[(x.channel == "online_storefront") & (~x.ambiguous)]95 x = x[x.groupby(["group_id", "week"]).banner.transform("nunique") >= 2]96 lg = np.log(x[pdef]); gap = lg.groupby([x.group_id, x.week]).agg(lambda s: s.max() - s.min())97 sens.append(dict(tiers=tiers, price_def=pdef, group_weeks=len(gap), median_gap=gap.median(), share_ge_5pct=(gap >= .05).mean()))98 for inc in ["with_circulars"]:99 x = pm_all[~pm_all.ambiguous]; x = x[x.groupby(["group_id", "week"]).banner.transform("nunique") >= 2]100 gap = x.groupby(["group_id", "week"]).log_p.agg(lambda s: s.max() - s.min()); sens.append(dict(tiers="AB", price_def="p_mode+circulars", group_weeks=len(gap), median_gap=gap.median(), share_ge_5pct=(gap >= .05).mean()))101 pd.DataFrame(sens).to_csv(P / "modelA_sensitivity.csv", index=False)102 (E / "modelA_summary.json").write_text(json.dumps(out, indent=2, default=float))103 log.info("H1: median gap %.3f (same-co %.3f, cross-co %.3f); share>=5%% %.2f; AR1 %.3f (se %.3f), placebo %.3f; sign-stable pairs %.2f (n=%d)",104 out["gap_all"]["median"], out["gap_same_company"]["median"], out["gap_cross_company"]["median"], out["gap_all"]["share_ge_5pct"], rho, se_rho, out["placebo_ar1_mean"], stable, npairs)105 ledger("4", "model_a", f"median gap {out['gap_all']['median']:.3f}; AR1 {rho:.3f}")106107if __name__ == "__main__":108 main()109