1"""Phase 6 — robustness and falsification checks that are not already inside the model scripts.23Runs: balanced common-product panel (H1); unit-price normalisation (H1 on log unit prices); outlier trimming (1/99 %);4category heterogeneity (H1); exclusion of same-company pairs; promotion-week exclusion (H1 gaps on non-promoted weeks only);5Holm correction over the preregistered primary tests; attrition table across stages. Results → data/processed/robustness_*.csv6and experiments/robustness_summary.json.7"""8import json, numpy as np, pandas as pd9from scipy import stats10from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED11log = get_logger("robust"); P = path("processed"); E = ROOT / "experiments"1213def gaps(pm, col="log_p"):14 return pm.groupby(["group_id", "week"])[col].agg(lambda s: s.max() - s.min())1516def main():17 seed_everything(); out = {}18 pm = pd.read_parquet(P / "panel_matched_AB.parquet"); pm = pm[(pm.channel == "online_storefront") & (~pm.ambiguous)]19 pm = pm[pm.groupby(["group_id", "week"]).banner.transform("nunique") >= 2].copy()20 base = gaps(pm); out["baseline_median_gap"] = float(base.median())21 # 1 balanced panel: groups with >= 2 banners in every one of the 11 weeks22 nwk = pm.groupby("group_id").week.nunique(); bal_groups = nwk[nwk == pm.week.nunique()].index23 bal = pm[pm.group_id.isin(bal_groups)]; gb = gaps(bal)24 out["balanced"] = dict(groups=int(len(bal_groups)), median_gap=float(gb.median()), share_ge_5pct=float((gb >= .05).mean()))25 # 2 unit-price normalisation26 pu = pm[pm.log_unit_p.notna()]; gu = gaps(pu, "log_unit_p"); out["unit_price"] = dict(group_weeks=int(len(gu)), median_gap=float(gu.median()), share_ge_5pct=float((gu >= .05).mean()))27 # 3 trimming28 lo, hi = base.quantile(.01), base.quantile(.99); tr = base[(base >= lo) & (base <= hi)]; out["trimmed_1_99"] = dict(median_gap=float(tr.median()), mean_gap=float(tr.mean()), mean_untrimmed=float(base.mean()))29 # 4 category heterogeneity (with Holm-adjusted test that median gap >= 5 % via sign test on gap-0.05)30 rows = []31 for c, g in pm.groupby("category"):32 gg = gaps(g)33 if len(gg) < 50: continue34 k = int((gg > .05).sum()); n = int(len(gg)); p = stats.binomtest(k, n, 0.5, alternative="greater").pvalue35 rows.append(dict(category=c, group_weeks=n, median_gap=float(gg.median()), share_ge_5pct=float((gg >= .05).mean()), p_sign_test=p))36 cat = pd.DataFrame(rows).sort_values("median_gap", ascending=False)37 m = len(cat); order = np.argsort(cat.p_sign_test.values); adj = np.empty(m)38 for rank, i in enumerate(order): adj[i] = min(1.0, cat.p_sign_test.values[i] * (m - rank))39 cat["p_holm"] = np.maximum.accumulate(adj[order])[np.argsort(order)]; cat.to_csv(P / "robustness_category_gaps.csv", index=False)40 # 5 same-company exclusion (cross-company only) and same-company only41 xc = pm[pm.groupby(["group_id", "week"]).company.transform("nunique") >= 2]; gx = gaps(xc)42 out["cross_company_only"] = dict(group_weeks=int(len(gx)), median_gap=float(gx.median()))43 # 6 promotion exclusion: drop listing-weeks at <= 0.9 x trailing mode44 def rolling_mode(s):45 v = s.tolist(); o = []46 for i in range(len(v)):47 vc = pd.Series(v[max(0, i - 7): i + 1]).value_counts(); o.append(max(vc[vc == vc.max()].index))48 return pd.Series(o, index=s.index)49 pm = pm.sort_values(["uid", "week"]); pm["ref"] = pm.groupby("uid").p_mode.transform(rolling_mode); pm["promo"] = pm.p_mode <= 0.9 * pm.ref50 np_ = pm[~pm.promo]; np_ = np_[np_.groupby(["group_id", "week"]).banner.transform("nunique") >= 2]; gn = gaps(np_)51 out["non_promo_weeks_only"] = dict(group_weeks=int(len(gn)), median_gap=float(gn.median()), share_ge_5pct=float((gn >= .05).mean()))52 po = pm[pm.groupby(["group_id", "week"]).promo.transform("max")]; gp = gaps(po); out["group_weeks_with_a_promo"] = dict(group_weeks=int(len(gp)), median_gap=float(gp.median()))53 # 7 Holm over preregistered primary tests (one p per hypothesis family where a p-value exists)54 A = json.load(open(E / "modelA_summary.json")); prim = {}55 # H1: sign test that median gap > 5 %56 k = int((base > .05).sum()); prim["H1_median_gap_gt_5pct"] = stats.binomtest(k, len(base), 0.5, alternative="greater").pvalue57 prim["H1_AR1_gt_0.5"] = float(1 - stats.norm.cdf((A["ar1_rho"] - 0.5) / A["ar1_se"]))58 if (E / "modelB_summary.json").exists():59 B = json.load(open(E / "modelB_summary.json"))60 lo_, hi_ = B["ratio_ci95"]; prim["H2_ratio_ne_1"] = float(2 * min(np.mean(np.array([lo_, hi_]) <= 1), 0.5)) if not (lo_ <= 1 <= hi_) else 1.061 prim["H2_ratio_ne_1"] = 0.0 if not (lo_ <= 1 <= hi_) else 1.0 # CI-based (bootstrap): reject at 5 % iff CI excludes 162 prim["H3_return_episodes_dlog_ref_ne_0"] = float(2 * (1 - stats.norm.cdf(abs(B["H3_return_dlog_ref_mean"] / B["H3_return_dlog_ref_se"]))))63 prim["H3_all_low_price_episodes_dlog_ref_ne_0"] = float(2 * (1 - stats.norm.cdf(abs(B["H3_dlog_ref_mean"] / B["H3_dlog_ref_se"]))))64 if (E / "modelC_summary.json").exists():65 C = json.load(open(E / "modelC_summary.json"))66 c = C["rival_all"]; z = c["coefs"]["rival_all_l1"] / c["ses"]["rival_all_l1"]; prim["H4_lag1_ne_0"] = float(2 * (1 - stats.norm.cdf(abs(z))))67 if (E / "modelD_summary.json").exists():68 D = json.load(open(E / "modelD_summary.json")); lo_, hi_ = D["bootstrap"]["gap_ci95"]; prim["H5_naive_vs_jevons_gap_ne_0"] = 0.0 if not (lo_ <= 0 <= hi_) else 1.069 if (E / "modelE_summary.json").exists():70 Ee = json.load(open(E / "modelE_summary.json")); prim["H6_K2_captures_all"] = np.nan # descriptive share; no test71 ps = pd.Series({k: v for k, v in prim.items() if not (isinstance(v, float) and np.isnan(v))}).sort_values()72 holm = pd.DataFrame({"test": ps.index, "p": ps.values}); mm = len(holm); holm["p_holm"] = np.minimum(1, np.maximum.accumulate(holm.p.values * (mm - np.arange(mm))))73 holm.to_csv(P / "robustness_holm.csv", index=False); out["holm"] = holm.to_dict(orient="records")74 (E / "robustness_summary.json").write_text(json.dumps(out, indent=2, default=float))75 log.info("%s", json.dumps({k: v for k, v in out.items() if k != "holm"}, indent=None, default=float)[:900]); log.info("\n%s", holm.to_string(index=False))76 ledger("6", "robustness", "robustness tables written")7778if __name__ == "__main__":79 main()80