1"""Model B — promotion dynamics (H2, H3). Preregistration §2 H2–H3.23Storefront promotion state is inferred: week promoted if p_mode <= 0.90 x reference, reference = modal p_mode over the4trailing 8 weeks (incl. current). Circular items carry observed regular prices (original_price) and are analysed separately.5"""6import json, numpy as np, pandas as pd, pyfixest as pf7from scipy import stats8from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED9log = get_logger("modelB"); P = path("processed"); E = ROOT / "experiments"1011def rolling_mode(vals, window=8):12 out = []13 for i in range(len(vals)):14 w = vals[max(0, i - window + 1): i + 1]15 vc = pd.Series(w).value_counts()16 top = vc[vc == vc.max()].index17 out.append(max(top)) # ties -> highest price (reference = regular)18 return out1920def main():21 seed_everything(); out = {}22 pw = pd.read_parquet(P / "panel_weekly.parquet")23 sf = pw[pw.channel == "online_storefront"].sort_values(["uid", "week"]).copy()24 nw = sf.groupby("uid").week.transform("nunique"); sf = sf[nw >= 8].copy()25 log.info("H2/H3 storefront sample: %d listing-weeks, %d listings, %d banners", len(sf), sf.uid.nunique(), sf.banner.nunique())26 sf["ref"] = sf.groupby("uid").p_mode.transform(lambda s: pd.Series(rolling_mode(s.tolist()), index=s.index))27 sf["promo"] = sf.p_mode <= 0.90 * sf.ref28 sf["depth"] = np.where(sf.promo, np.log(sf.ref / sf.p_mode), 0.0)29 sf["log_ref"] = np.log(sf.ref)30 sf["dlog_p"] = sf.groupby("uid").log_p.diff(); sf["dlog_ref"] = sf.groupby("uid").log_ref.diff()31 sf["promo_lag"] = sf.groupby("uid").promo.shift(1)32 # validation vs current on_sale flag (listings active now; last complete week)33 last = sf[sf.week == sf.week.max()]; last = last[last.active == 1]34 agree = (last.promo == (last.on_sale_now == 1)).mean()35 out["promo_vs_on_sale_now_agreement"] = float(agree); out["promo_rate_last_week"] = float(last.promo.mean()); out["on_sale_now_rate"] = float((last.on_sale_now == 1).mean())36 out["promo_share_listing_weeks"] = float(sf.promo.mean()); out["promo_depth_median"] = float(sf.loc[sf.promo, "depth"].median()); out["promo_depth_mean"] = float(sf.loc[sf.promo, "depth"].mean())37 # H2: volatility raw vs reference38 chg = sf.dropna(subset=["dlog_p"])39 vol = chg.groupby("banner").agg(n=("dlog_p", "size"), sd_raw=("dlog_p", "std"), sd_ref=("dlog_ref", "std"), freq_change=("dlog_p", lambda s: (s.abs() > 1e-9).mean()),40 freq_ref_change=("dlog_ref", lambda s: (s.abs() > 1e-9).mean()),41 share_changes_promo=("promo", lambda s: np.nan)).reset_index()42 chg["is_change"] = chg.dlog_p.abs() > 1e-9; chg["promo_transition"] = chg.promo != chg.promo_lag43 vol["share_changes_promo"] = vol.banner.map(chg[chg.is_change].groupby("banner").promo_transition.mean())44 vol["ratio_raw_ref"] = vol.sd_raw / vol.sd_ref; vol.sort_values("n", ascending=False).to_csv(P / "modelB_volatility_by_banner.csv", index=False)45 volc = chg.groupby("category").agg(n=("dlog_p", "size"), sd_raw=("dlog_p", "std"), sd_ref=("dlog_ref", "std")).reset_index(); volc["ratio"] = volc.sd_raw / volc.sd_ref46 volc.to_csv(P / "modelB_volatility_by_category.csv", index=False)47 out["sd_raw_all"] = float(chg.dlog_p.std()); out["sd_ref_all"] = float(chg.dlog_ref.std()); out["ratio_all"] = out["sd_raw_all"] / out["sd_ref_all"]48 out["share_weekly_changes_promo_transition"] = float(chg[chg.is_change].promo_transition.mean()); out["freq_weekly_change"] = float(chg.is_change.mean())49 # bootstrap over listings for the ratio50 rng = np.random.default_rng(SEED); uids = chg.uid.unique(); byuid = {u: g for u, g in chg.groupby("uid")}51 rat = []52 for b in range(300):53 samp = rng.choice(uids, size=len(uids), replace=True); d = pd.concat([byuid[u] for u in samp])54 rat.append(d.dlog_p.std() / d.dlog_ref.std())55 out["ratio_ci95"] = [float(np.percentile(rat, 2.5)), float(np.percentile(rat, 97.5))]56 # Falsification: single-store online shops (Shopify/WooCommerce) vs anti-bot storefronts — share of multi-price weeks57 single = ["giant_tiger", "pa", "epipresto", "nuvo", "boite_a_grains", "bocoboco", "aliments_merci", "akhavan", "mayrand", "aubut", "maturin", "tau", "loco", "avril"]58 out["multi_price_week_share_single_store"] = float(sf[sf.source.isin(single)].multi_price_week.mean()); out["multi_price_week_share_antibot"] = float(sf[~sf.source.isin(single)].multi_price_week.mean())59 # Transition matrix60 tm = pd.crosstab(sf.promo_lag.dropna().astype(bool), sf.loc[sf.promo_lag.notna(), "promo"], normalize="index"); tm.to_csv(P / "modelB_transition_matrix.csv")61 out["P_stay_promo"] = float(tm.loc[True, True]) if True in tm.index else None; out["P_enter_promo"] = float(tm.loc[False, True]) if False in tm.index else None62 # Episodes63 eps = []64 for uid, g in sf.groupby("uid"):65 g = g.reset_index(drop=True); pr = g.promo.values; i = 066 while i < len(pr):67 if pr[i]:68 j = i69 while j + 1 < len(pr) and pr[j + 1]: j += 170 censored_start = i == 0; censored_end = j == len(pr) - 171 pre_ref = g.ref.iloc[i - 1] if i > 0 else np.nan; post_ref = g.ref.iloc[j + 1] if j + 1 < len(pr) else np.nan72 post_p = g.p_mode.iloc[j + 1] if j + 1 < len(pr) else np.nan73 eps.append(dict(uid=uid, banner=g.banner.iloc[0], category=g.category.iloc[0], start=g.week.iloc[i], end=g.week.iloc[j], duration=j - i + 1,74 depth=g.depth.iloc[i:j + 1].mean(), censored=censored_start or censored_end, pre_ref=pre_ref, post_ref=post_ref, post_p=post_p,75 dlog_ref=np.log(post_ref / pre_ref) if pre_ref and post_ref and not (censored_start or censored_end) else np.nan,76 dlog_post_price_vs_pre_ref=np.log(post_p / pre_ref) if pre_ref and post_p and not (censored_start or censored_end) else np.nan))77 i = j + 178 else: i += 179 ep = pd.DataFrame(eps); ep.to_parquet(P / "modelB_episodes.parquet", index=False)80 unc = ep[~ep.censored]81 out["n_episodes"] = int(len(ep)); out["n_uncensored"] = int(len(unc)); out["duration_median_unc"] = float(unc.duration.median()); out["duration_mean_unc"] = float(unc.duration.mean())82 out["share_duration_1w"] = float((unc.duration == 1).mean()); out["depth_median_unc"] = float(unc.depth.median())83 rec = ep.sort_values(["uid", "start"]); rec["prev_end"] = rec.groupby("uid").end.shift(1); rec["recur_weeks"] = (rec.start - rec.prev_end).dt.days / 784 out["recurrence_median_weeks"] = float(rec.recur_weeks.dropna().median()) if rec.recur_weeks.notna().any() else None; out["share_listings_with_2plus_episodes"] = float((ep.groupby("uid").size() >= 2).mean())85 # H3: reference change after episode86 d = unc.dlog_ref.dropna(); out["H3_dlog_ref_mean"] = float(d.mean()); out["H3_dlog_ref_se"] = float(d.std() / np.sqrt(len(d))); out["H3_n"] = int(len(d))87 out["H3_share_ref_up"] = float((d > 1e-9).mean()); out["H3_share_ref_down"] = float((d < -1e-9).mean()); out["H3_share_ref_same"] = float((d.abs() <= 1e-9).mean())88 d2 = unc.dlog_post_price_vs_pre_ref.dropna(); out["H3_post_price_vs_pre_ref_mean"] = float(d2.mean()); out["H3_post_price_vs_pre_ref_share_up"] = float((d2 > 1e-9).mean())89 # Event study around episode end (weeks -4..+4), listing FE + week FE; pretrend joint test90 ev = []91 for _, e in unc.iterrows():92 g = sf[sf.uid == e.uid]93 for k in range(-4, 5):94 wk = e.end + pd.Timedelta(weeks=k)95 r = g[g.week == wk]96 if len(r): ev.append(dict(uid=e.uid, week=wk, k=k, log_p=float(r.log_p.iloc[0]), log_ref=float(r.log_ref.iloc[0])))97 ev = pd.DataFrame(ev)98 if len(ev):99 ev["k_f"] = ev.k.astype(int)100 fit = pf.feols("log_p ~ i(k_f, ref=-1) | uid + week", data=ev, vcov={"CRV1": "uid"})101 es = pd.DataFrame({"term": fit.coef().index, "coef": fit.coef().values, "se": fit.se().values}); es.to_csv(P / "modelB_event_study.csv", index=False)102 pre = [t for t in fit.coef().index if any(t.endswith(f"[T.{k}]") for k in (-4, -3, -2))]103 try:104 w = fit.wald_test(R=None, q=None) if False else None105 except Exception: w = None106 b = fit.coef()[pre].values; V = fit._vcov; idx = [list(fit.coef().index).index(t) for t in pre]; Vs = V[np.ix_(idx, idx)]107 W = float(b @ np.linalg.solve(Vs, b)); pval = float(1 - stats.chi2.cdf(W, len(pre)))108 out["pretrend_wald"] = W; out["pretrend_p"] = pval; out["event_study_n"] = int(fit._N)109 fitr = pf.feols("log_ref ~ i(k_f, ref=-1) | uid + week", data=ev, vcov={"CRV1": "uid"})110 pd.DataFrame({"term": fitr.coef().index, "coef": fitr.coef().values, "se": fitr.se().values}).to_csv(P / "modelB_event_study_ref.csv", index=False)111 # Placebo: random pseudo-episode ends for non-promoted listings112 rng = np.random.default_rng(SEED); pl = []113 nonp = sf[~sf.promo].groupby("uid").filter(lambda g: len(g) >= 9)114 for uid, g in list(nonp.groupby("uid"))[:3000]:115 g = g.sort_values("week"); i = rng.integers(2, len(g) - 2); pre = g.ref.iloc[i - 1]; post = g.ref.iloc[i + 1]116 pl.append(np.log(post / pre))117 pl = np.array(pl); out["placebo_dlog_ref_mean"] = float(pl.mean()); out["placebo_dlog_ref_se"] = float(pl.std() / np.sqrt(len(pl))); out["placebo_n"] = int(len(pl))118 # Circular channel: observed regular prices119 circ = pd.read_parquet(P / "products.parquet"); ret = pd.read_parquet(P / "retailers.parquet")120 circ = circ.merge(ret[["feed_id", "channel", "banner", "is_grocery"]], left_on="source", right_on="feed_id")121 circ = circ[circ.channel.str.startswith("circular") & circ.is_grocery & circ.price.notna()]122 circ["depth_obs"] = np.log(circ.regular_price / circ.price)123 out["circular_items"] = int(len(circ)); out["circular_share_with_regular"] = float(circ.regular_price.notna().mean())124 out["circular_depth_median"] = float(circ.depth_obs.dropna().median()); out["circular_depth_mean"] = float(circ.depth_obs.dropna().mean())125 out["circular_share_recurring_items"] = float((circ.span_days > 8).mean())126 circ.groupby("banner").agg(items=("uid", "size"), with_regular=("regular_price", lambda s: s.notna().mean()), depth_median=("depth_obs", "median"), recurring=("span_days", lambda s: (s > 8).mean())).reset_index().to_csv(P / "modelB_circulars_by_banner.csv", index=False)127 (E / "modelB_summary.json").write_text(json.dumps(out, indent=2, default=float))128 log.info("H2: sd raw %.4f vs ref %.4f (ratio %.2f, CI %s); promo share %.3f depth med %.3f; agreement w/ on_sale %.2f", out["sd_raw_all"], out["sd_ref_all"], out["ratio_all"], np.round(out["ratio_ci95"], 2), out["promo_share_listing_weeks"], out["promo_depth_median"], agree)129 log.info("H3: episodes %d (uncensored %d), duration median %.1f, dlog_ref mean %.4f (se %.4f) up %.2f down %.2f same %.2f; placebo %.4f; pretrend p %.3f",130 out["n_episodes"], out["n_uncensored"], out["duration_median_unc"], out["H3_dlog_ref_mean"], out["H3_dlog_ref_se"], out["H3_share_ref_up"], out["H3_share_ref_down"], out["H3_share_ref_same"], out["placebo_dlog_ref_mean"], out.get("pretrend_p", float("nan")))131 ledger("4", "model_b", f"ratio {out['ratio_all']:.2f}; episodes {out['n_episodes']}; H3 mean {out['H3_dlog_ref_mean']:.4f}")132133if __name__ == "__main__":134 main()135