SPB Git
spboucher.ai

spb/wp37_uqo

Public

phase4: models A (dispersion), C (spillovers), D (indices), E (basket MIP), F (forecasts); pulp pinned <3

Simon-Pierre Boucher committed (Nov 1, 2026) parent 8931742 Browse files

17 changed files +930 −2

17 changed files
Unified Split
added experiments/modelA_summary.json +59 −0
@@ -0,0 +1,59 @@
1 +{
2 + "gap_all": {
3 + "n": 30042,
4 + "median": 0.15415067982725827,
5 + "p25": 0.05901417143705534,
6 + "p75": 0.27866231037551104,
7 + "mean": 0.18936608290604118,
8 + "share_ge_5pct": 0.7851341455295919,
9 + "share_ge_10pct": 0.6325144797283803
10 + },
11 + "gap_same_company": {
12 + "n": 27939,
13 + "median": 0.15128944754622609,
14 + "p25": 0.057746822212323856,
15 + "p75": 0.26844510500133234,
16 + "mean": 0.18576393643510483,
17 + "share_ge_5pct": 0.7821682952145746,
18 + "share_ge_10pct": 0.6273667633057733
19 + },
20 + "gap_cross_company": {
21 + "n": 2103,
22 + "median": 0.18254405095440074,
23 + "p25": 0.08017111858864046,
24 + "p75": 0.3370446464822916,
25 + "mean": 0.23722170356771002,
26 + "share_ge_5pct": 0.8245363766048502,
27 + "share_ge_10pct": 0.7009034712315739
28 + },
29 + "cv_all": {
30 + "n": 30042,
31 + "median": 0.07341658510817389,
32 + "p25": 0.029429075927015894,
33 + "p75": 0.13034226381318856,
34 + "mean": 0.08951038693130896,
35 + "share_ge_5pct": 0.6284867851674323,
36 + "share_ge_10pct": 0.3644231409360229
37 + },
38 + "banner_fe_reference": "Metro",
39 + "resid_sd": 0.11456460507159695,
40 + "n_obs_fe": 67732,
41 + "ar1_rho": 0.663411776499718,
42 + "ar1_se": 0.013989188479587735,
43 + "share_pairs_sign_stable_ge4w": 0.457002457002457,
44 + "n_pairs_ge4w": 4884,
45 + "ar1_same_company": [
46 + 0.6675381113692006,
47 + 0.013865058825916333,
48 + 0.45115007970849463,
49 + 4391
50 + ],
51 + "ar1_cross_company": [
52 + 0.6302714331760938,
53 + 0.03769326769881746,
54 + 0.5091277890466531,
55 + 493
56 + ],
57 + "placebo_ar1_mean": 0.0035750286222568906,
58 + "placebo_ar1_sd": 0.011048954753771625
59 +}
\ No newline at end of file
added experiments/modelC_summary.json +71 −0
@@ -0,0 +1,71 @@
1 +{
2 + "rival_all": {
3 + "n": 22213,
4 + "coefs": {
5 + "rival_all": 0.036983611990011225,
6 + "rival_all_l1": 0.03674370006018662,
7 + "rival_all_l2": 0.014867128258090253,
8 + "rival_all_l3": 0.04341238029758149,
9 + "rival_all_f1": 0.008899934370312498,
10 + "rival_all_f2": 0.013343162704560697
11 + },
12 + "ses": {
13 + "rival_all": 0.018207550537640733,
14 + "rival_all_l1": 0.012276621362321993,
15 + "rival_all_l2": 0.013107075427690282,
16 + "rival_all_l3": 0.014656351032885128,
17 + "rival_all_f1": 0.012048599868225932,
18 + "rival_all_f2": 0.014684867882272816
19 + },
20 + "P_own_change_given_rival_change_l1": 0.2527131782945736,
21 + "P_own_change_given_no_rival_change_l1": 0.1711969770843491
22 + },
23 + "rival_same": {
24 + "n": 21192,
25 + "coefs": {
26 + "rival_same": 0.027632240932606675,
27 + "rival_same_l1": 0.04166870039441454,
28 + "rival_same_l2": 0.01145626003174879,
29 + "rival_same_l3": 0.034260009805555276,
30 + "rival_same_f1": 0.0132625155508785,
31 + "rival_same_f2": 0.014714728442878371
32 + },
33 + "ses": {
34 + "rival_same": 0.01831796970903035,
35 + "rival_same_l1": 0.01234955835151791,
36 + "rival_same_l2": 0.013246905947066992,
37 + "rival_same_l3": 0.015028588991760062,
38 + "rival_same_f1": 0.012251592712678097,
39 + "rival_same_f2": 0.014337631789336497
40 + },
41 + "P_own_change_given_rival_change_l1": 0.25367852486496556,
42 + "P_own_change_given_no_rival_change_l1": 0.17588320798837137
43 + },
44 + "rival_cross": {
45 + "n": 1784,
46 + "coefs": {
47 + "rival_cross": 0.08084154607798574,
48 + "rival_cross_l1": -0.0012395703633851167,
49 + "rival_cross_l2": -0.056863578932101456,
50 + "rival_cross_l3": 0.046743672971606504,
51 + "rival_cross_f1": -0.10799246348174829,
52 + "rival_cross_f2": -0.064020894931646
53 + },
54 + "ses": {
55 + "rival_cross": 0.06613440280922223,
56 + "rival_cross_l1": 0.04779603312485358,
57 + "rival_cross_l2": 0.028928467120645663,
58 + "rival_cross_l3": 0.06013080127915635,
59 + "rival_cross_f1": 0.048390891586986506,
60 + "rival_cross_f2": 0.04882941084495431
61 + },
62 + "P_own_change_given_rival_change_l1": 0.2744360902255639,
63 + "P_own_change_given_no_rival_change_l1": 0.18370607028753994
64 + },
65 + "promo_comovement": {
66 + "group_weeks": 30042,
67 + "share_any_low": 0.16027561414020372,
68 + "share_two_plus_low": 0.017209240396777844,
69 + "expected_two_plus_if_independent": 0.006207116594008007
70 + }
71 +}
\ No newline at end of file
added experiments/modelD_summary.json +59 −0
@@ -0,0 +1,59 @@
1 +{
2 + "storefront_posted": {
3 + "naive_mean_change": 0.0036769234266957884,
4 + "naive_median_change": 0.0,
5 + "jevons_chained_change": 0.0028498796180160824,
6 + "fixed_jevons_change": 0.000961438796386993,
7 + "fixed_dutot_change": -0.0005377093740115413,
8 + "fixed_carli_change": 0.026867158956678905,
9 + "n_balanced": 44423,
10 + "n_total_first_week": 51434
11 + },
12 + "storefront_unit": {
13 + "naive_mean_change": -0.012446313570385947,
14 + "naive_median_change": -0.006544530339934296,
15 + "jevons_chained_change": 0.0027664755512584183,
16 + "fixed_jevons_change": 0.0007395269633583279,
17 + "fixed_dutot_change": -0.004284566992962091,
18 + "fixed_carli_change": 0.03566130323976746,
19 + "n_balanced": 29671,
20 + "n_total_first_week": 36003
21 + },
22 + "all_channels_posted": {
23 + "naive_mean_change": -0.02927666201945356,
24 + "naive_median_change": 0.0,
25 + "jevons_chained_change": 0.0023785412285990404,
26 + "fixed_jevons_change": 0.0005360460992416272,
27 + "fixed_dutot_change": -0.0008775599602732775,
28 + "fixed_carli_change": 0.027756326994421432,
29 + "n_balanced": 45214,
30 + "n_total_first_week": 57942
31 + },
32 + "bootstrap": {
33 + "gap_mean": 0.0006271123459269828,
34 + "gap_ci95": [
35 + -0.007830383632211136,
36 + 0.008431239690139423
37 + ],
38 + "jevons_ci95": [
39 + 0.001509196870153784,
40 + 0.004221513313892511
41 + ]
42 + },
43 + "assortment": {
44 + "listings": 71694,
45 + "balanced": 44515,
46 + "entrants": 20260,
47 + "exits": 17295,
48 + "weekly_entry_rate": 0.03595173117340451
49 + },
50 + "size_label_history_available": false,
51 + "statcan": {
52 + "cpi_qc_food_stores_aug": 194.4,
53 + "sep": 193.1,
54 + "cpi_change_aug_sep": -0.006687242798353976,
55 + "ours_jevons_aug_sep": 0.000991433026574784,
56 + "cpi_latest": "2026-09-01",
57 + "cpi_yoy_sep": 0.016315789473684283
58 + }
59 +}
\ No newline at end of file
added experiments/modelE_summary.json +38 −0
@@ -0,0 +1,38 @@
1 +{
2 + "items_in_common_set_mean": 7.181818181818182,
3 + "banners_available_mean": 7.7272727272727275,
4 + "frontier_mean_cost_by_K": {
5 + "1": 23.12194937122785,
6 + "2": 20.888101262563136,
7 + "3": 20.204123130937315,
8 + "4": 20.054139457868775,
9 + "5": 20.034210164939484,
10 + "99": 20.028755619484937
11 + },
12 + "savings_unconstrained_vs_best_single": 3.0931937517429127,
13 + "savings_share_unconstrained": 0.1337773775939487,
14 + "share_of_unconstrained_savings_at_K2": 0.7221817603265276,
15 + "single_banner_cost_spread_mean": 39.69805844961497,
16 + "single_banner_mean_cost": 32.88295457993973,
17 + "most_expensive_vs_cheapest_single_share": 1.856238474704022,
18 + "travel_cost_sensitivity": [
19 + {
20 + "fixed_cost": 5.0,
21 + "n_used": 1.0,
22 + "cost_goods": 18.12194937122785,
23 + "cost_incl": 23.12194937122785
24 + },
25 + {
26 + "fixed_cost": 10.0,
27 + "n_used": 1.0,
28 + "cost_goods": 13.12194937122785,
29 + "cost_incl": 23.12194937122785
30 + },
31 + {
32 + "fixed_cost": 20.0,
33 + "n_used": 1.0,
34 + "cost_goods": 3.121949371227848,
35 + "cost_incl": 23.12194937122785
36 + }
37 + ]
38 +}
\ No newline at end of file
added experiments/modelF_summary.json +104 −0
@@ -0,0 +1,104 @@
1 +{
2 + "summary": [
3 + {
4 + "target": "log_price",
5 + "model": "banner_cat_drift",
6 + "origins": 4,
7 + "n": 215264,
8 + "mae": 0.024807352788469867,
9 + "rmse": 0.10842948605770154,
10 + "brier": NaN,
11 + "auc": NaN
12 + },
13 + {
14 + "target": "log_price",
15 + "model": "lightgbm",
16 + "origins": 4,
17 + "n": 215264,
18 + "mae": 0.030072224891589998,
19 + "rmse": 0.10001695748926878,
20 + "brier": NaN,
21 + "auc": NaN
22 + },
23 + {
24 + "target": "log_price",
25 + "model": "mean_reversion_ref",
26 + "origins": 4,
27 + "n": 215264,
28 + "mae": 0.02837850853015459,
29 + "rmse": 0.10249483142934793,
30 + "brier": NaN,
31 + "auc": NaN
32 + },
33 + {
34 + "target": "log_price",
35 + "model": "naive_last",
36 + "origins": 4,
37 + "n": 215264,
38 + "mae": 0.022605082941485766,
39 + "rmse": 0.10825508687878052,
40 + "brier": NaN,
41 + "auc": NaN
42 + },
43 + {
44 + "target": "log_price",
45 + "model": "ridge",
46 + "origins": 4,
47 + "n": 215264,
48 + "mae": 0.029885922259458478,
49 + "rmse": 0.09917790589117646,
50 + "brier": NaN,
51 + "auc": NaN
52 + },
53 + {
54 + "target": "promotion",
55 + "model": "base_rate",
56 + "origins": 4,
57 + "n": 215264,
58 + "mae": NaN,
59 + "rmse": NaN,
60 + "brier": 0.0568232667303276,
61 + "auc": 0.5
62 + },
63 + {
64 + "target": "promotion",
65 + "model": "hist_rate",
66 + "origins": 4,
67 + "n": 215264,
68 + "mae": NaN,
69 + "rmse": NaN,
70 + "brier": 0.057192760775280654,
71 + "auc": 0.6778171962496753
72 + },
73 + {
74 + "target": "promotion",
75 + "model": "lightgbm",
76 + "origins": 4,
77 + "n": 215264,
78 + "mae": NaN,
79 + "rmse": NaN,
80 + "brier": 0.03127230289344953,
81 + "auc": 0.9359091370790923
82 + },
83 + {
84 + "target": "promotion",
85 + "model": "logit",
86 + "origins": 4,
87 + "n": 215264,
88 + "mae": NaN,
89 + "rmse": NaN,
90 + "brier": 0.04392957451269847,
91 + "auc": 0.8939986946425332
92 + },
93 + {
94 + "target": "promotion",
95 + "model": "persist",
96 + "origins": 4,
97 + "n": 215264,
98 + "mae": NaN,
99 + "rmse": NaN,
100 + "brier": 0.050992947231181564,
101 + "auc": 0.7518564786841642
102 + }
103 + ]
104 +}
\ No newline at end of file
modified pyproject.toml +1 −1
@@ -14,7 +14,7 @@ dependencies = [
14 14 "lightgbm>=4.3",
15 15 "matplotlib>=3.9",
16 16 "duckdb>=1.0",
17 − "pulp>=2.8",
17 + "pulp>=2.8,<3",
18 18 "pyyaml>=6.0",
19 19 "pytest>=8.0",
20 20 "requests>=2.32",
modified requirements-lock.txt +1 −1
@@ -57,7 +57,7 @@ pluggy==1.7.0
57 57 prompt-toolkit==3.0.53
58 58 psutil==7.2.2
59 59 ptyprocess==0.7.0
60 −pulp==4.2.0
60 +pulp==2.9.0
61 61 pure-eval==0.2.4
62 62 pyarrow==26.0.0
63 63 pycparser==3.12
added src/econometrics/__init__.py +0 −0
added src/econometrics/build_matched_panel.py +58 −0
@@ -0,0 +1,58 @@
1 +"""Phase 4.0 — analysis-ready weekly panels (preregistration §1).
2 +
3 +Outputs (data/processed):
4 + panel_weekly.parquet — grocery-universe listing × complete circular week, posted price = weekly mode, unit prices,
5 + banner/company/channel/category, size, flags
6 + panel_matched_<TIERS>.parquet — panel_weekly restricted to listings in cross-banner match groups (A, AB, ABC, ABCD),
7 + with group_id and `ambiguous` flag (group has >1 simultaneously posted listing from the
8 + same banner-channel in that week)
9 +"""
10 +import numpy as np, pandas as pd
11 +from src.common import path, get_logger, ledger, SETTINGS, write_manifest
12 +log = get_logger("panel4")
13 +P = path("processed")
14 +
15 +def main():
16 + wk = pd.read_parquet(P / "offers_weekly.parquet")
17 + prod = pd.read_parquet(P / "products.parquet")
18 + ret = pd.read_parquet(P / "retailers.parquet")
19 + prod = prod.merge(ret[["feed_id", "banner", "company", "channel", "store_format", "is_grocery"]], left_on="source", right_on="feed_id", how="left")
20 + sc = SETTINGS["scope"]
21 + uni = prod[(~prod.source.isin(sc["exclude_sources"])) & prod.is_grocery & (~prod.category.isin(sc["non_grocery_categories"])) & (~prod.flag_unrealistic)]
22 + log.info("grocery universe listings: %d (of %d)", len(uni), len(prod))
23 + weeks = sorted(wk.week.unique())
24 + complete = [w for w in weeks if w >= pd.Timestamp("2026-08-13").date() and w <= pd.Timestamp("2026-10-22").date()]
25 + pw = wk[wk.week.isin(complete) & wk.uid.isin(uni.uid)].copy()
26 + cols = ["uid", "source", "banner", "company", "channel", "store_format", "category", "name", "brand", "size_label", "size_qty", "size_base", "size_pack",
27 + "per_weight_label", "flag_weighted", "flag_multibuy", "flag_loyalty", "flag_online_only", "gtin_valid", "active", "on_sale", "regular_price", "price"]
28 + pw = pw.merge(uni[cols].rename(columns={"source": "source_p", "price": "price_now", "on_sale": "on_sale_now", "regular_price": "regular_now"}), on="uid", how="left").drop(columns=["source_p"])
29 + pw["week"] = pd.to_datetime(pw.week); pw["week_idx"] = pw.week.rank(method="dense").astype(int)
30 + pw["log_p"] = np.log(pw.p_mode)
31 + qty = pw.size_qty.where(pw.size_base != "un", pw.size_qty)
32 + pw["unit_p"] = np.where(pw.size_base == "un", pw.p_mode / pw.size_qty, pw.p_mode / pw.size_qty * 100.0)
33 + pw.loc[pw.per_weight_label.eq("kg"), "unit_p"] = pw.p_mode / 10.0
34 + pw.loc[pw.per_weight_label.eq("lb"), "unit_p"] = pw.p_mode / 4.53592
35 + pw["unit_base"] = np.where(pw.size_base == "un", "un", np.where(pw.size_base == "ml", "100 ml", "100 g"))
36 + pw.loc[pw.per_weight_label.isin(["kg", "lb"]), "unit_base"] = "100 g"
37 + pw["log_unit_p"] = np.log(pw.unit_p)
38 + pw["multi_price_week"] = pw.n_distinct > 1
39 + pw.to_parquet(P / "panel_weekly.parquet", index=False)
40 + log.info("panel_weekly: %d listing-weeks, %d listings, %d weeks (%s → %s), %d banners", len(pw), pw.uid.nunique(), pw.week.nunique(), pw.week.min().date(), pw.week.max().date(), pw.banner.nunique())
41 + files = [P / "panel_weekly.parquet"]
42 + for tiers in ["A", "AB", "ABC", "ABCD"]:
43 + g = pd.read_parquet(P / f"match_groups_{tiers}.parquet")[["group_id", "uid"]]
44 + pm = pw.merge(g, on="uid", how="inner")
45 + # keep group-weeks with >= 2 banners posting
46 + nb = pm.groupby(["group_id", "week"]).banner.transform("nunique"); pm = pm[nb >= 2].copy()
47 + dup = pm.groupby(["group_id", "week", "banner", "channel"]).uid.transform("size") > 1
48 + amb = pm.assign(d=dup).groupby(["group_id", "week"]).d.transform("max"); pm["ambiguous"] = amb.astype(bool)
49 + pm["same_company_group"] = pm.groupby(["group_id", "week"]).company.transform("nunique") == 1
50 + pm.to_parquet(P / f"panel_matched_{tiers}.parquet", index=False); files.append(P / f"panel_matched_{tiers}.parquet")
51 + log.info("panel_matched_%-4s: %d listing-weeks, %d groups, %d group-weeks, ambiguous group-weeks %.1f%%, same-company-only group-weeks %.1f%%",
52 + tiers, len(pm), pm.group_id.nunique(), pm.groupby(["group_id", "week"]).ngroups,
53 + 100 * pm.drop_duplicates(["group_id", "week"]).ambiguous.mean(), 100 * pm.drop_duplicates(["group_id", "week"]).same_company_group.mean())
54 + write_manifest("panel4", files)
55 + ledger("4", "build_matched_panel", f"{len(pw)} listing-weeks; weeks {pw.week.min().date()}..{pw.week.max().date()}")
56 +
57 +if __name__ == "__main__":
58 + main()
added src/econometrics/model_a_dispersion.py +108 −0
@@ -0,0 +1,108 @@
1 +"""Model A — matched-product cross-banner price dispersion (H1). Preregistration §2 H1.
2 +
3 +Outputs: data/processed/modelA_*.csv and experiments/modelA_summary.json
4 +"""
5 +import json, re, numpy as np, pandas as pd, pyfixest as pf
6 +from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED
7 +log = get_logger("modelA"); P = path("processed"); E = ROOT / "experiments"; E.mkdir(exist_ok=True)
8 +
9 +def 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 st
15 +
16 +def 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"])
27 +
28 +def 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()), ok
40 +
41 +def 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 + # Persistence
69 + 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"] = npairs
72 + 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 heatmap
81 + 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 definitions
91 + 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}")
106 +
107 +if __name__ == "__main__":
108 + main()
added src/econometrics/model_b_promotions.py +134 −0
@@ -0,0 +1,134 @@
1 +"""Model B — promotion dynamics (H2, H3). Preregistration §2 H2–H3.
2 +
3 +Storefront promotion state is inferred: week promoted if p_mode <= 0.90 x reference, reference = modal p_mode over the
4 +trailing 8 weeks (incl. current). Circular items carry observed regular prices (original_price) and are analysed separately.
5 +"""
6 +import json, numpy as np, pandas as pd, pyfixest as pf
7 +from scipy import stats
8 +from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED
9 +log = get_logger("modelB"); P = path("processed"); E = ROOT / "experiments"
10 +
11 +def 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()].index
17 + out.append(max(top)) # ties -> highest price (reference = regular)
18 + return out
19 +
20 +def 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.ref
28 + 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 reference
38 + 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_lag
43 + 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_ref
46 + 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 ratio
50 + 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 weeks
57 + 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 matrix
60 + 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 None
62 + # Episodes
63 + eps = []
64 + for uid, g in sf.groupby("uid"):
65 + g = g.reset_index(drop=True); pr = g.promo.values; i = 0
66 + while i < len(pr):
67 + if pr[i]:
68 + j = i
69 + while j + 1 < len(pr) and pr[j + 1]: j += 1
70 + censored_start = i == 0; censored_end = j == len(pr) - 1
71 + 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.nan
72 + post_p = g.p_mode.iloc[j + 1] if j + 1 < len(pr) else np.nan
73 + 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 + 1
78 + else: i += 1
79 + 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 / 7
84 + 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 episode
86 + 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 test
90 + 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 None
105 + except Exception: w = None
106 + 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 listings
112 + 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 prices
119 + 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}")
132 +
133 +if __name__ == "__main__":
134 + main()
added src/econometrics/model_c_spillovers.py +63 −0
@@ -0,0 +1,63 @@
1 +"""Model C — competitive spillovers (H4), descriptive lead–lag only. Preregistration §2 H4.
2 +
3 +Within matched groups (A+B, storefront), regress a banner's weekly log-price change on rivals' current and lagged
4 +changes (k = 0..3), with group-week and banner fixed effects; leads (k = -1, -2) as placebo. Rival change = mean change of
5 +the other banners in the group that week. Same-company vs cross-company rivals estimated separately. No causal claim:
6 +contemporaneous terms absorb common shocks; lagged terms are lead–lag associations.
7 +"""
8 +import json, numpy as np, pandas as pd, pyfixest as pf
9 +from src.common import path, ROOT, get_logger, ledger, seed_everything
10 +log = get_logger("modelC"); P = path("processed"); E = ROOT / "experiments"
11 +
12 +def main():
13 + seed_everything(); out = {}
14 + pm = pd.read_parquet(P / "panel_matched_AB.parquet"); pm = pm[(pm.channel == "online_storefront") & (~pm.ambiguous)].copy()
15 + pm = pm.drop_duplicates(["group_id", "week", "banner"]).sort_values(["group_id", "banner", "week"])
16 + pm["dlp"] = pm.groupby(["group_id", "banner"]).log_p.diff()
17 + pm["gap_w"] = (pm.week - pm.groupby(["group_id", "banner"]).week.shift(1)).dt.days
18 + pm.loc[pm.gap_w != 7, "dlp"] = np.nan
19 + # rival mean change, same-company vs cross-company
20 + tot = pm.groupby(["group_id", "week"]).dlp.transform(lambda s: s.sum(min_count=1)); cnt = pm.groupby(["group_id", "week"]).dlp.transform("count")
21 + pm["rival_all"] = (tot - pm.dlp.fillna(0)) / (cnt - pm.dlp.notna().astype(int)).replace(0, np.nan)
22 + def rival_by(same):
23 + rows = []
24 + for (g, w), d in pm.groupby(["group_id", "week"]):
25 + for i, r in d.iterrows():
26 + oth = d[(d.index != i) & ((d.company == r.company) == same)].dlp.dropna()
27 + rows.append((i, oth.mean() if len(oth) else np.nan))
28 + return pd.Series(dict(rows))
29 + pm["rival_same"] = rival_by(True); pm["rival_cross"] = rival_by(False)
30 + for col in ["rival_all", "rival_same", "rival_cross"]:
31 + for k in (1, 2, 3):
32 + pm[f"{col}_l{k}"] = pm.groupby(["group_id", "banner"])[col].shift(k)
33 + for k in (1, 2):
34 + pm[f"{col}_f{k}"] = pm.groupby(["group_id", "banner"])[col].shift(-k)
35 + pm["gw"] = pm.group_id + "_" + pm.week.dt.strftime("%Y%m%d")
36 + res = []
37 + for col in ["rival_all", "rival_same", "rival_cross"]:
38 + terms = [col] + [f"{col}_l{k}" for k in (1, 2, 3)] + [f"{col}_f{k}" for k in (1, 2)]
39 + d = pm.dropna(subset=["dlp"] + terms)
40 + if len(d) < 500: log.warning("%s: only %d obs", col, len(d)); continue
41 + fit = pf.feols("dlp ~ " + " + ".join(terms) + " | banner + week", data=d, vcov={"CRV1": "group_id"})
42 + for t in terms:
43 + res.append(dict(rival=col, term=t, coef=float(fit.coef()[t]), se=float(fit.se()[t]), n=int(fit._N)))
44 + out[col] = dict(n=int(fit._N), coefs={t: float(fit.coef()[t]) for t in terms}, ses={t: float(fit.se()[t]) for t in terms})
45 + # descriptive: probability a banner changes price in w given rival changed in w-1
46 + d["rival_changed_l1"] = d[f"{col}_l1"].abs() > 1e-9; d["own_changed"] = d.dlp.abs() > 1e-9
47 + ct = d.groupby("rival_changed_l1").own_changed.mean()
48 + out[col]["P_own_change_given_rival_change_l1"] = float(ct.get(True, np.nan)); out[col]["P_own_change_given_no_rival_change_l1"] = float(ct.get(False, np.nan))
49 + pd.DataFrame(res).to_csv(P / "modelC_distributed_lags.csv", index=False)
50 + # Promotion co-movement: share of group-weeks in which >=2 banners are simultaneously at a low (p <= 0.9 x own 8-week mode)
51 + pm["ref"] = pm.groupby(["group_id", "banner"]).p_mode.transform(lambda s: s.rolling(8, min_periods=1).apply(lambda x: pd.Series(x).mode().max()))
52 + pm["low"] = pm.p_mode <= 0.9 * pm.ref
53 + gl = pm.groupby(["group_id", "week"]).agg(n=("banner", "size"), nlow=("low", "sum"))
54 + gl = gl[gl.n >= 2]; out["promo_comovement"] = dict(group_weeks=int(len(gl)), share_any_low=float((gl.nlow >= 1).mean()), share_two_plus_low=float((gl.nlow >= 2).mean()),
55 + expected_two_plus_if_independent=float(((gl.nlow / gl.n).mean()) ** 2))
56 + (E / "modelC_summary.json").write_text(json.dumps(out, indent=2, default=float))
57 + for col in out:
58 + if col.startswith("rival"): log.info("%s: n=%d, k0 %.3f (%.3f), l1 %.3f (%.3f), l2 %.3f (%.3f), f1 %.3f (%.3f)", col, out[col]["n"], out[col]["coefs"][col], out[col]["ses"][col], out[col]["coefs"][col + "_l1"], out[col]["ses"][col + "_l1"], out[col]["coefs"][col + "_l2"], out[col]["ses"][col + "_l2"], out[col]["coefs"][col + "_f1"], out[col]["ses"][col + "_f1"])
59 + log.info("promo co-movement: %s", out["promo_comovement"])
60 + ledger("4", "model_c", "distributed lags estimated (descriptive)")
61 +
62 +if __name__ == "__main__":
63 + main()
added src/forecasting/__init__.py +0 −0
added src/forecasting/rolling_forecast.py +65 −0
@@ -0,0 +1,65 @@
1 +"""Model F — rolling-origin forecasts of next-week posted log price and promotion indicator (preregistration §2 F).
2 +
3 +Features use only information up to week t (strict lags). Origins t = 6..10 (train weeks <= t, test week t+1).
4 +Models: seasonal naive (last value), banner×category mean reversion (shrink toward 8-week listing mode), ridge, LightGBM.
5 +Promotion target: next week p <= 0.9 x trailing 8-week mode. Metrics: MAE/RMSE (log price), Brier/AUC (promotion).
6 +"""
7 +import json, numpy as np, pandas as pd, lightgbm as lgb
8 +from sklearn.linear_model import Ridge, LogisticRegression
9 +from sklearn.metrics import roc_auc_score, brier_score_loss
10 +from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED
11 +log = get_logger("modelF"); P = path("processed"); E = ROOT / "experiments"
12 +
13 +def rolling_mode_series(s, window=8):
14 + out = []
15 + v = s.tolist()
16 + for i in range(len(v)):
17 + w = pd.Series(v[max(0, i - window + 1): i + 1]).value_counts(); out.append(max(w[w == w.max()].index))
18 + return pd.Series(out, index=s.index)
19 +
20 +def main():
21 + seed_everything(); out = {}
22 + pw = pd.read_parquet(P / "panel_weekly.parquet"); sf = pw[pw.channel == "online_storefront"].sort_values(["uid", "week"]).copy()
23 + sf = sf[sf.groupby("uid").week.transform("nunique") >= 6]
24 + sf["t"] = sf.week.rank(method="dense").astype(int)
25 + g = sf.groupby("uid")
26 + sf["ref"] = g.p_mode.transform(rolling_mode_series); sf["promo"] = (sf.p_mode <= 0.9 * sf.ref).astype(int)
27 + for k in (1, 2, 3): sf[f"lp_l{k}"] = g.log_p.shift(k)
28 + sf["lp_ref"] = np.log(sf.ref); sf["promo_l1"] = g.promo.shift(1); sf["promo_l2"] = g.promo.shift(2)
29 + sf["n_changes_l0"] = sf.n_changes; sf["rel_range_l0"] = sf.rel_range; sf["multi_l0"] = sf.multi_price_week.astype(int)
30 + sf["promo_rate_hist"] = g.promo.transform(lambda s: s.shift(1).expanding().mean())
31 + sf["y_lp"] = g.log_p.shift(-1); sf["y_promo"] = g.promo.shift(-1); sf["t_next"] = g.t.shift(-1)
32 + sf = sf[(sf.t_next == sf.t + 1)]
33 + sf["banner_c"] = sf.banner.astype("category").cat.codes; sf["cat_c"] = sf.category.astype("category").cat.codes
34 + feats = ["log_p", "lp_l1", "lp_l2", "lp_l3", "lp_ref", "promo", "promo_l1", "promo_l2", "n_changes_l0", "rel_range_l0", "multi_l0", "promo_rate_hist", "banner_c", "cat_c", "t"]
35 + res = []
36 + for t in range(6, int(sf.t.max())):
37 + tr = sf[(sf.t <= t - 1)].dropna(subset=feats + ["y_lp", "y_promo"]); te = sf[sf.t == t].dropna(subset=feats + ["y_lp", "y_promo"])
38 + if len(te) == 0 or len(tr) == 0: continue
39 + # price models
40 + preds = {"naive_last": te.log_p.values, "mean_reversion_ref": 0.5 * te.log_p.values + 0.5 * te.lp_ref.values}
41 + bc = tr.groupby(["banner", "category"]).apply(lambda d: (d.y_lp - d.log_p).mean()).rename("drift").reset_index()
42 + preds["banner_cat_drift"] = te.log_p.values + te.merge(bc, on=["banner", "category"], how="left").drift.fillna(0).values
43 + X = tr[feats].values; Xt = te[feats].values
44 + rd = Ridge(alpha=1.0).fit(X, tr.y_lp.values); preds["ridge"] = rd.predict(Xt)
45 + gb = lgb.LGBMRegressor(n_estimators=400, learning_rate=0.05, num_leaves=63, min_child_samples=50, subsample=0.8, colsample_bytree=0.8, random_state=SEED, verbose=-1).fit(X, tr.y_lp.values)
46 + preds["lightgbm"] = gb.predict(Xt)
47 + for m, p in preds.items():
48 + e = p - te.y_lp.values; res.append(dict(origin_t=t, target="log_price", model=m, n=len(te), mae=float(np.mean(np.abs(e))), rmse=float(np.sqrt(np.mean(e ** 2)))))
49 + # promotion models
50 + pp = {"persist": te.promo.values.astype(float), "base_rate": np.full(len(te), tr.y_promo.mean()), "hist_rate": te.promo_rate_hist.fillna(tr.y_promo.mean()).values}
51 + lr = LogisticRegression(max_iter=2000).fit(X, tr.y_promo.values); pp["logit"] = lr.predict_proba(Xt)[:, 1]
52 + gc = lgb.LGBMClassifier(n_estimators=400, learning_rate=0.05, num_leaves=63, min_child_samples=50, subsample=0.8, colsample_bytree=0.8, random_state=SEED, verbose=-1).fit(X, tr.y_promo.values)
53 + pp["lightgbm"] = gc.predict_proba(Xt)[:, 1]
54 + for m, p in pp.items():
55 + y = te.y_promo.values
56 + res.append(dict(origin_t=t, target="promotion", model=m, n=len(te), brier=float(brier_score_loss(y, np.clip(p, 0, 1))), auc=float(roc_auc_score(y, p)) if len(set(y)) > 1 else np.nan, base_rate=float(y.mean())))
57 + r = pd.DataFrame(res); r.to_csv(P / "modelF_rolling_results.csv", index=False)
58 + summ = r.groupby(["target", "model"]).agg(origins=("origin_t", "nunique"), n=("n", "sum"), mae=("mae", "mean"), rmse=("rmse", "mean"), brier=("brier", "mean"), auc=("auc", "mean")).reset_index()
59 + summ.to_csv(P / "modelF_summary.csv", index=False); out["summary"] = summ.to_dict(orient="records")
60 + (E / "modelF_summary.json").write_text(json.dumps(out, indent=2, default=float))
61 + log.info("\n%s", summ.round(4).to_string(index=False))
62 + ledger("4", "model_f", "rolling-origin forecasts done")
63 +
64 +if __name__ == "__main__":
65 + main()
added src/indices/__init__.py +0 −0
added src/indices/basket_optimization.py +83 −0
@@ -0,0 +1,83 @@
1 +"""Model E — fixed-basket shopping optimisation (H6). Preregistration §2 H6 and config/basket.yaml (frozen).
2 +
3 +For each complete week: map each basket item to candidate listings per banner (category + name regex, parsable size,
4 +unit price), take the cheapest unit price per banner, and solve min-cost assignment with at most K banners (K=1..5) as a
5 +MIP (PuLP/CBC). Travel costs are unobserved: an illustrative fixed per-banner cost grid is added as a labelled assumption.
6 +"""
7 +import json, re, numpy as np, pandas as pd, pulp, yaml
8 +from src.common import path, ROOT, get_logger, ledger
9 +log = get_logger("modelE"); P = path("processed"); E = ROOT / "experiments"
10 +
11 +def candidates(pw, items):
12 + rows = []
13 + for it in items:
14 + d = pw[(pw.category == it["category"]) & pw.name.fillna("").str.contains(it["name_regex"], case=False, regex=True)]
15 + if it.get("exclude_regex"): d = d[~d.name.fillna("").str.contains(it["exclude_regex"], case=False, regex=True)]
16 + d = d[(d.unit_base == ("un" if it["base"] == "un" else f"100 {it['base']}")) & d.unit_p.notna() & (d.unit_p > 0)]
17 + if len(d) == 0: continue
18 + d = d.assign(item=it["id"], item_cost=lambda x: x.unit_p * (it["qty"] if it["base"] == "un" else it["qty"] / 100.0))
19 + rows.append(d[["item", "week", "banner", "company", "uid", "name", "size_label", "p_mode", "unit_p", "item_cost"]])
20 + return pd.concat(rows, ignore_index=True)
21 +
22 +def solve(cost_tbl, K, fixed_cost=0.0):
23 + """cost_tbl: item x banner matrix of min item cost (NaN = unavailable). Returns (total cost, banners used, items covered)."""
24 + items = list(cost_tbl.index); banners = list(cost_tbl.columns)
25 + prob = pulp.LpProblem("basket", pulp.LpMinimize)
26 + x = {(i, b): pulp.LpVariable(f"x_{n}_{m}", cat="Binary") for n, i in enumerate(items) for m, b in enumerate(banners) if not np.isnan(cost_tbl.at[i, b])}
27 + y = {b: pulp.LpVariable(f"y_{m}", cat="Binary") for m, b in enumerate(banners)}
28 + prob += pulp.lpSum(cost_tbl.at[i, b] * v for (i, b), v in x.items()) + fixed_cost * pulp.lpSum(y.values())
29 + for i in items:
30 + prob += pulp.lpSum(v for (ii, b), v in x.items() if ii == i) == 1
31 + for (i, b), v in x.items():
32 + prob += v <= y[b]
33 + prob += pulp.lpSum(y.values()) <= K
34 + prob.solve(pulp.PULP_CBC_CMD(msg=False))
35 + if pulp.LpStatus[prob.status] != "Optimal": return np.nan, [], 0
36 + used = [b for b, v in y.items() if v.value() > 0.5]
37 + tot = sum(cost_tbl.at[i, b] * v.value() for (i, b), v in x.items())
38 + return float(tot), used, len(items)
39 +
40 +def main():
41 + out = {}
42 + items = yaml.safe_load((ROOT / "config" / "basket.yaml").read_text())["items"]
43 + pw = pd.read_parquet(P / "panel_weekly.parquet"); sf = pw[pw.channel == "online_storefront"]
44 + cand = candidates(sf, items); cand.to_parquet(P / "modelE_candidates.parquet", index=False)
45 + cov = cand.groupby("item").agg(listings=("uid", "nunique"), banners=("banner", "nunique")).reset_index(); cov.to_csv(P / "modelE_item_coverage.csv", index=False)
46 + log.info("basket items with candidates: %d / %d ; median banners per item %.0f", cov.item.nunique(), len(items), cov.banners.median())
47 + results = []
48 + for w, d in cand.groupby("week"):
49 + tbl = d.groupby(["item", "banner"]).item_cost.min().unstack()
50 + # keep banners carrying >= 60 % of covered items, and items available at >= 2 banners (so that K=1 is feasible on a common set)
51 + tbl = tbl.loc[:, tbl.notna().mean() >= 0.6]
52 + tbl = tbl[tbl.notna().sum(axis=1) >= 2]
53 + # for K=1 feasibility, restrict to items available at every retained banner
54 + full = tbl.dropna()
55 + for K in (1, 2, 3, 4, 5, len(tbl.columns)):
56 + tot, used, n = solve(full, K)
57 + results.append(dict(week=w, K=K if K < len(tbl.columns) else 99, cost=tot, banners=";".join(used), n_items=n, n_banners_avail=len(tbl.columns)))
58 + # single-banner cost for each banner (full item set)
59 + for b in full.columns:
60 + results.append(dict(week=w, K=1, cost=float(full[b].sum()), banners=b, n_items=len(full), n_banners_avail=len(tbl.columns), single=True))
61 + for fc in (5.0, 10.0, 20.0):
62 + tot, used, n = solve(full, len(full.columns), fixed_cost=fc)
63 + results.append(dict(week=w, K=99, cost=tot - fc * len(used), banners=";".join(used), n_items=n, n_banners_avail=len(tbl.columns), fixed_cost=fc, n_used=len(used), cost_incl_travel=tot))
64 + r = pd.DataFrame(results); r.to_csv(P / "modelE_results.csv", index=False)
65 + base = r[(r.get("single").isna()) & (r.fixed_cost.isna())] if "fixed_cost" in r else r[r.get("single").isna()]
66 + front = base.groupby("K").cost.mean(); singles = r[r.get("single") == True].groupby("week").cost.agg(["min", "max", "mean"])
67 + out["items_in_common_set_mean"] = float(base.n_items.mean()); out["banners_available_mean"] = float(base.n_banners_avail.mean())
68 + out["frontier_mean_cost_by_K"] = {int(k): float(v) for k, v in front.items()}
69 + c1, c2, cinf = front.get(1), front.get(2), front.get(99)
70 + out["savings_unconstrained_vs_best_single"] = float(c1 - cinf); out["savings_share_unconstrained"] = float((c1 - cinf) / c1)
71 + out["share_of_unconstrained_savings_at_K2"] = float((c1 - c2) / (c1 - cinf)) if c1 - cinf > 0 else None
72 + out["single_banner_cost_spread_mean"] = float((singles["max"] - singles["min"]).mean()); out["single_banner_mean_cost"] = float(singles["mean"].mean())
73 + out["most_expensive_vs_cheapest_single_share"] = float(((singles["max"] - singles["min"]) / singles["min"]).mean())
74 + if "fixed_cost" in r:
75 + tc = r[r.fixed_cost.notna()].groupby("fixed_cost").agg(n_used=("n_used", "mean"), cost_goods=("cost", "mean"), cost_incl=("cost_incl_travel", "mean")).reset_index()
76 + tc.to_csv(P / "modelE_travel_sensitivity.csv", index=False); out["travel_cost_sensitivity"] = tc.to_dict(orient="records")
77 + (E / "modelE_summary.json").write_text(json.dumps(out, indent=2, default=float))
78 + log.info("E: common items %.1f, banners %.1f ; cost K=1 %.2f, K=2 %.2f, unconstrained %.2f ; unconstrained savings %.1f%% ; share captured at K=2 %.2f ; single-banner spread %.2f (%.1f%%)",
79 + out["items_in_common_set_mean"], out["banners_available_mean"], c1, c2, cinf, 100 * out["savings_share_unconstrained"], out["share_of_unconstrained_savings_at_K2"] or float("nan"), out["single_banner_cost_spread_mean"], 100 * out["most_expensive_vs_cheapest_single_share"])
80 + ledger("4", "model_e", f"K1 {c1:.2f} K2 {c2:.2f} Kinf {cinf:.2f}")
81 +
82 +if __name__ == "__main__":
83 + main()
added src/indices/price_indices.py +86 −0
@@ -0,0 +1,86 @@
1 +"""Model D — price indices (H5). Preregistration §2 H5. No expenditure weights exist: all indices are unweighted and are
2 +NOT consumer price indices. Size-label history is not observed (only current state) → within-listing package changes are
3 +unobservable; this is reported as a limitation, not estimated.
4 +"""
5 +import json, numpy as np, pandas as pd
6 +from src.common import path, ROOT, get_logger, ledger, seed_everything, SEED
7 +log = get_logger("modelD"); P = path("processed"); E = ROOT / "experiments"
8 +
9 +def jevons_chained(df, col):
10 + """Chained Jevons (matched-model): geometric mean of week-on-week ratios over listings present in both weeks."""
11 + weeks = sorted(df.week.unique()); idx = [100.0]; n = [np.nan]
12 + piv = df.pivot_table(index="uid", columns="week", values=col)
13 + for a, b in zip(weeks, weeks[1:]):
14 + both = piv[[a, b]].dropna(); r = np.exp((np.log(both[b]) - np.log(both[a])).mean()) if len(both) else np.nan
15 + idx.append(idx[-1] * r); n.append(len(both))
16 + return pd.DataFrame({"week": weeks, "index": idx, "n_matched": n})
17 +
18 +def fixed_basket(df, col):
19 + """Balanced listings (present in every week), equal weights: Dutot (ratio of means), Carli (mean of ratios), Jevons (geo mean of ratios)."""
20 + weeks = sorted(df.week.unique()); piv = df.pivot_table(index="uid", columns="week", values=col).dropna()
21 + base = piv[weeks[0]]
22 + rows = []
23 + for w in weeks:
24 + rows.append(dict(week=w, n=len(piv), dutot=100 * piv[w].mean() / base.mean(), carli=100 * (piv[w] / base).mean(), jevons=100 * np.exp(np.log(piv[w] / base).mean())))
25 + return pd.DataFrame(rows)
26 +
27 +def naive_mean(df, col):
28 + g = df.groupby("week")[col].agg(["mean", "median", "size"]).reset_index(); g["index_mean"] = 100 * g["mean"] / g["mean"].iloc[0]; g["index_median"] = 100 * g["median"] / g["median"].iloc[0]
29 + return g
30 +
31 +def main():
32 + seed_everything(); out = {}
33 + pw = pd.read_parquet(P / "panel_weekly.parquet"); sf = pw[pw.channel == "online_storefront"].copy()
34 + res = {}
35 + for name, d, col in [("storefront_posted", sf, "p_mode"), ("storefront_unit", sf[sf.unit_p.notna()], "unit_p"), ("all_channels_posted", pw, "p_mode")]:
36 + nm = naive_mean(d, col); jc = jevons_chained(d, col); fb = fixed_basket(d, col)
37 + m = nm.merge(jc, on="week").merge(fb, on="week", suffixes=("", "_fb")); m["sample"] = name; res[name] = m
38 + out[name] = dict(naive_mean_change=float(nm.index_mean.iloc[-1] / 100 - 1), naive_median_change=float(nm.index_median.iloc[-1] / 100 - 1), jevons_chained_change=float(jc["index"].iloc[-1] / 100 - 1),
39 + fixed_jevons_change=float(fb.jevons.iloc[-1] / 100 - 1), fixed_dutot_change=float(fb.dutot.iloc[-1] / 100 - 1), fixed_carli_change=float(fb.carli.iloc[-1] / 100 - 1), n_balanced=int(fb.n.iloc[0]), n_total_first_week=int(nm["size"].iloc[0]))
40 + pd.concat(res.values()).to_csv(P / "modelD_indices.csv", index=False)
41 + # by category (storefront posted)
42 + bc = []
43 + for c, d in sf.groupby("category"):
44 + if d.uid.nunique() < 100: continue
45 + nm = naive_mean(d, "p_mode"); jc = jevons_chained(d, "p_mode"); fb = fixed_basket(d, "p_mode")
46 + bc.append(dict(category=c, listings=d.uid.nunique(), naive_mean_change=nm.index_mean.iloc[-1] / 100 - 1, jevons_chained_change=jc["index"].iloc[-1] / 100 - 1, fixed_jevons_change=fb.jevons.iloc[-1] / 100 - 1, n_balanced=fb.n.iloc[0]))
47 + pd.DataFrame(bc).to_csv(P / "modelD_indices_by_category.csv", index=False)
48 + # by banner (storefront posted): Jevons chained
49 + bb = []
50 + for b, d in sf.groupby("banner"):
51 + if d.uid.nunique() < 100: continue
52 + jc = jevons_chained(d, "p_mode"); fb = fixed_basket(d, "p_mode"); nm = naive_mean(d, "p_mode")
53 + bb.append(dict(banner=b, listings=d.uid.nunique(), jevons_chained_change=jc["index"].iloc[-1] / 100 - 1, fixed_jevons_change=fb.jevons.iloc[-1] / 100 - 1, naive_mean_change=nm.index_mean.iloc[-1] / 100 - 1, n_balanced=fb.n.iloc[0]))
54 + pd.DataFrame(bb).to_csv(P / "modelD_indices_by_banner.csv", index=False)
55 + # Bootstrap (listing-level, B=999) of the naive-vs-Jevons gap, storefront posted
56 + rng = np.random.default_rng(SEED); uids = sf.uid.unique(); piv = sf.pivot_table(index="uid", columns="week", values="p_mode"); weeks = list(piv.columns)
57 + gaps = []; jev = []
58 + for b in range(999):
59 + s = rng.choice(uids, size=len(uids), replace=True); pb = piv.loc[s]
60 + naive = pb[weeks[-1]].mean() / pb[weeks[0]].mean() - 1
61 + lj = 0.0
62 + for a, c in zip(weeks, weeks[1:]):
63 + both = pb[[a, c]].dropna(); lj += (np.log(both[c]) - np.log(both[a])).mean()
64 + jev.append(np.exp(lj) - 1); gaps.append(naive - (np.exp(lj) - 1))
65 + out["bootstrap"] = dict(gap_mean=float(np.mean(gaps)), gap_ci95=[float(np.percentile(gaps, 2.5)), float(np.percentile(gaps, 97.5))], jevons_ci95=[float(np.percentile(jev, 2.5)), float(np.percentile(jev, 97.5))])
66 + # Assortment turnover (entry/exit) — the observable part of H5
67 + first = sf.groupby("uid").week.min(); last = sf.groupby("uid").week.max(); weeks_s = sorted(sf.week.unique())
68 + out["assortment"] = dict(listings=int(sf.uid.nunique()), balanced=int(((first == weeks_s[0]) & (last == weeks_s[-1])).sum()), entrants=int((first > weeks_s[0]).sum()), exits=int((last < weeks_s[-1]).sum()),
69 + weekly_entry_rate=float(sf.groupby("week").uid.nunique().pipe(lambda s: (first.value_counts().reindex(weeks_s).fillna(0)[1:] / s.values[:-1]).mean())))
70 + # Unit-price vs posted: cross-banner unit price gap for matched groups is covered in Model A; here, size-missing share
71 + out["size_label_history_available"] = False
72 + # StatCan benchmark: Quebec CPI 'Food purchased from stores' Aug->Sep 2026 ; our storefront Jevons Aug (weeks 13,20,27 Aug) -> Sep (3,10,17,24 Sep)
73 + cpi = pd.read_csv(P / "statcan_cpi_food.csv"); q = cpi[(cpi.geo == "Quebec") & (cpi["product"] == "Food purchased from stores")].sort_values("ref_date")
74 + aug = q[q.ref_date.str.startswith("2026-08")].value.iloc[0]; sep = q[q.ref_date.str.startswith("2026-09")].value.iloc[0]
75 + jc = jevons_chained(sf, "p_mode"); jc["month"] = pd.to_datetime(jc.week).dt.month
76 + ours = jc.groupby("month")["index"].mean(); out["statcan"] = dict(cpi_qc_food_stores_aug=float(aug), sep=float(sep), cpi_change_aug_sep=float(sep / aug - 1), ours_jevons_aug_sep=float(ours.loc[9] / ours.loc[8] - 1),
77 + cpi_latest=str(q.ref_date.iloc[-1]), cpi_yoy_sep=float(sep / q[q.ref_date.str.startswith("2025-09")].value.iloc[0] - 1))
78 + (E / "modelD_summary.json").write_text(json.dumps(out, indent=2, default=float))
79 + s = out["storefront_posted"]
80 + log.info("D storefront: naive mean %.3f, naive median %.3f, Jevons chained %.3f (CI %s), fixed Jevons %.3f (n=%d), Dutot %.3f, Carli %.3f; gap naive−Jevons %.3f CI %s",
81 + s["naive_mean_change"], s["naive_median_change"], s["jevons_chained_change"], np.round(out["bootstrap"]["jevons_ci95"], 3), s["fixed_jevons_change"], s["n_balanced"], s["fixed_dutot_change"], s["fixed_carli_change"], out["bootstrap"]["gap_mean"], np.round(out["bootstrap"]["gap_ci95"], 3))
82 + log.info("StatCan QC food-from-stores Aug→Sep %.3f vs ours %.3f ; CPI yoy Sep %.3f", out["statcan"]["cpi_change_aug_sep"], out["statcan"]["ours_jevons_aug_sep"], out["statcan"]["cpi_yoy_sep"])
83 + ledger("4", "model_d", f"Jevons {s['jevons_chained_change']:.3f} vs naive {s['naive_mean_change']:.3f}")
84 +
85 +if __name__ == "__main__":
86 + main()
87
Menu

Theme
↑↓ navigate ↵ open esc close / or ⌘K anywhere