RetailPulse

import pandas as pd, numpy as np
import matplotlib.pyplot as plt
from scipy import stats

fact   = pd.read_csv("data/fact_sales_weekly.csv", parse_dates=["week_start"])
prod   = pd.read_csv("data/dim_product.csv")
region = pd.read_csv("data/dim_region.csv")
promo  = pd.read_csv("data/dim_promo.csv", parse_dates=["week_start"])

df = (fact.merge(prod, on="product_id")
          .merge(region, on="region_id")
          .merge(promo[["promo_id", "discount_pct", "promo_type", "vendor_funding_pct"]], on="promo_id", how="left"))
df["week_idx"] = ((df.week_start - df.week_start.min()).dt.days // 7).astype(int)
df["woy"]      = df.week_idx % 52
df["on_promo"] = df.promo_id.notna().astype(int)
# net unit cost after vendor funding (national-brand promos are 50% supplier-funded)
df["net_unit_cost"] = np.where(df.units_sold > 0, df.cogs / df.units_sold.clip(lower=1), df.unit_cost)

print(f"{len(df):,} rows | {df.product_id.nunique()} products | "
      f"{df.region_id.nunique()} regions | {df.week_idx.nunique()} weeks")
print(f"Revenue ${df.revenue.sum()/1e6:,.1f}M | promo weeks stocked out: "
      f"{df.loc[df.on_promo==1,'stockout_flag'].mean():.0%} vs "
      f"{df.loc[df.on_promo==0,'stockout_flag'].mean():.1%} on regular weeks")

# Each item is sold as a FreshCart and a National-brand version ("siblings").
sib = prod.merge(prod[["product_id", "item"]], on="item", suffixes=("", "_sib"))
sib = sib[sib.product_id != sib.product_id_sib][["product_id", "product_id_sib"]]

promo_weeks = df[df.on_promo == 1][["product_id", "region_id", "week_idx"]].drop_duplicates()

# Was this product on promo LAST week? (captures pantry-loading / pull-forward)
post = promo_weeks.assign(week_idx=promo_weeks.week_idx + 1, post_promo=1)
df = df.merge(post, on=["product_id", "region_id", "week_idx"], how="left")

# Is the sibling product on promo THIS week? (captures cannibalization)
sp = promo_weeks.rename(columns={"product_id": "product_id_sib"}).merge(sib, on="product_id_sib")
sp = sp[["product_id", "region_id", "week_idx"]].assign(sib_promo=1)
df = df.merge(sp, on=["product_id", "region_id", "week_idx"], how="left")
df[["post_promo", "sib_promo"]] = df[["post_promo", "sib_promo"]].fillna(0).astype(int)
print(df[["on_promo", "post_promo", "sib_promo"]].mean().round(3).to_string())

PROMO_TYPES = ["Weekly ad", "In-store display", "Digital coupon"]

def design(d):
    X = pd.DataFrame({"const": 1.0, "log_price": np.log(d.price_paid)}, index=d.index)
    for t in PROMO_TYPES:                        # promo lift beyond the price cut
        X["promo_" + t] = ((d.promo_type == t) & (d.on_promo == 1)).astype(float)
    X["post_promo"] = d.post_promo
    X["sib_promo"]  = d.sib_promo
    X["trend"]      = d.week_idx / 52
    for r in (2, 3, 4):  X[f"region_{r}"] = (d.region_id == r).astype(float)
    for w in range(1, 52): X[f"woy_{w}"] = (d.woy == w).astype(float)   # seasonality
    return X

def fit(d):
    d = d[(d.stockout_flag == 0) & (d.units_sold > 0)]  # stockouts censor demand: drop them
    X, y = design(d), np.log(d.units_sold)
    beta, *_ = np.linalg.lstsq(X.values, y.values, rcond=None)
    resid = y.values - X.values @ beta
    cov = (resid @ resid / (len(y) - X.shape[1])) * np.linalg.pinv(X.T.values @ X.values)
    return pd.Series(beta, X.columns), pd.Series(np.sqrt(np.clip(np.diag(cov), 0, None)), X.columns)

models, rows = {}, []
for pid, d in df.groupby("product_id"):
    b, se = fit(d)
    models[pid] = b
    rows.append({"product_id": pid, "elasticity": b.log_price, "se": se.log_price,
                 **{t: np.exp(b["promo_" + t]) for t in PROMO_TYPES},
                 "post_promo_effect": np.exp(b.post_promo) - 1,
                 "sib_promo_effect": np.exp(b.sib_promo) - 1})
el = prod.merge(pd.DataFrame(rows), on="product_id")
print(el.groupby("category").elasticity.agg(["mean", "min", "max"]).round(2))

# Because the data is synthetic, we know the true elasticity used to generate it.
truth = pd.read_csv("truth_elasticity.csv")
chk = el.merge(truth, on="product_id")
mae  = (chk.elasticity - chk.true_elasticity).abs().mean()
r    = np.corrcoef(chk.elasticity, chk.true_elasticity)[0, 1]
cover = ((chk.true_elasticity > chk.elasticity - 1.96*chk.se) &
         (chk.true_elasticity < chk.elasticity + 1.96*chk.se)).mean()
print(f"Mean absolute error: {mae:.3f} | correlation: {r:.3f} | 95% CI coverage: {cover:.0%}")

fig, ax = plt.subplots(figsize=(6.2, 4.4))
for c, g in chk.groupby("category"):
    ax.scatter(g.true_elasticity, g.elasticity, s=26, label=c, alpha=.85)
lim = [chk[["elasticity", "true_elasticity"]].min().min()-.2, -0.3]
ax.plot(lim, lim, "--", color="#9CA3AF", lw=1)
ax.set(xlabel="True elasticity", ylabel="Estimated elasticity",
       title="Estimated vs. true price elasticity (120 products)")
ax.legend(fontsize=8, frameon=False); plt.tight_layout(); plt.savefig("charts/elasticity_validation.png", dpi=150)

cat = el.groupby("category").agg(e=("elasticity", "mean"), se=("elasticity", "sem")).sort_values("e")
fig, ax = plt.subplots(figsize=(6.2, 3.6))
ax.barh(cat.index, cat.e, xerr=1.96*cat.se, color="#14B8A6", ecolor="#0B1F3A", capsize=3)
ax.axvline(-1, color="#DC2626", ls="--", lw=1)
ax.text(-1.03, 5.45, "elastic  ←  |  →  inelastic", ha="center", fontsize=8, color="#DC2626")
ax.set(xlabel="Average price elasticity (95% CI)", title="How price-sensitive is each category?")
plt.tight_layout(); plt.savefig("charts/elasticity_by_category.png", dpi=150)
print(cat.round(2))

def predict(pid, d):
    return np.exp(design(d).values @ models[pid].values)

pw_ = {k: g for k, g in df.groupby(["product_id", "week_idx"])}   # fast lookup
promo_week = df.dropna(subset=["promo_id"]).groupby("promo_id").week_idx.first()
ev = []
for p in promo.itertuples():
    w  = promo_week[p.promo_id]
    me = pw_[(p.product_id, w)]
    # Baseline: what the model says would have sold at list price, no promo
    base = me.copy(); base["price_paid"] = base.list_price; base["on_promo"] = 0; base["promo_type"] = None
    base_units = predict(p.product_id, base)
    demand     = predict(p.product_id, me)              # uncensored promo demand
    sold       = me.units_sold.values
    promo_margin = (sold * (me.price_paid.values - me.net_unit_cost.values)).sum()   # after vendor funding
    base_margin  = (base_units * (me.list_price.values - me.unit_cost.values)).sum()
    # Cannibalization: sibling units lost this week
    sid = sib.product_id_sib[sib.product_id == p.product_id].iloc[0]
    s = pw_[(sid, w)]
    s0 = s.copy(); s0["sib_promo"] = 0
    cann = ((predict(sid, s0) - predict(sid, s)) * (s.price_paid.values - s.unit_cost.values)).sum()
    # Pull-forward: own units lost next week
    n = pw_.get((p.product_id, w + 1), df.iloc[:0])
    n0 = n.copy(); n0["post_promo"] = 0
    pull = ((predict(p.product_id, n0) - predict(p.product_id, n)) * (n.price_paid.values - n.unit_cost.values)).sum() if len(n) else 0
    ev.append({"promo_id": p.promo_id, "product_id": p.product_id, "week": str(p.week_start.date()),
               "discount_pct": p.discount_pct, "promo_type": p.promo_type,
               "units": sold.sum(), "baseline_units": base_units.sum(),
               "lift_pct": sold.sum() / base_units.sum() - 1,
               "lost_sales_units": np.clip(demand - sold, 0, None)[me.stockout_flag.values == 1].sum(),
               "gross_incr_margin": promo_margin - base_margin,
               "cannibalization": cann, "pull_forward": pull,
               "discount_spend": (sold * (me.list_price.values - me.price_paid.values)).sum()
                                 * (1 - p.vendor_funding_pct)})            # retailer-funded share only
sc = pd.DataFrame(ev).merge(prod[["product_id", "product_name", "category", "brand"]], on="product_id")
sc["net_incr_margin"] = sc.gross_incr_margin - sc.cannibalization - sc.pull_forward
sc["roi"] = sc.net_incr_margin / sc.discount_spend
sc["verdict"] = np.select([sc.net_incr_margin > 0, sc.gross_incr_margin > 0], ["Repeat", "Fix"], "Stop")
print(sc.verdict.value_counts().to_string())
print(f"Net incremental margin from all promos: ${sc.net_incr_margin.sum()/1e3:,.0f}K "
      f"on ${sc.discount_spend.sum()/1e6:,.2f}M of discounts")
print(sc.groupby("discount_pct").agg(promos=("roi", "size"), avg_lift=("lift_pct", "mean"),
      avg_roi=("roi", "mean"), pct_repeat=("verdict", lambda v: (v == "Repeat").mean())).round(2))

piv = sc.pivot_table(index="category", columns="discount_pct", values="roi", aggfunc="mean")
fig, ax = plt.subplots(figsize=(6.4, 3.8))
im = ax.imshow(piv.values, cmap="RdYlGn", vmin=-1, vmax=1, aspect="auto")
ax.set_xticks(range(len(piv.columns)), [f"{c:.0%}" for c in piv.columns])
ax.set_yticks(range(len(piv.index)), piv.index)
for i in range(piv.shape[0]):
    for j in range(piv.shape[1]):
        ax.text(j, i, f"{piv.values[i, j]:.0%}", ha="center", va="center", fontsize=8)
ax.set(xlabel="Discount depth", title="Promo ROI (net incremental margin ÷ discount spend)")
plt.colorbar(im, fraction=.03); plt.tight_layout(); plt.savefig("charts/promo_roi_heatmap.png", dpi=150)

train, test = df[df.week_idx < 92], df[df.week_idx >= 92].copy()
ly = df[["product_id", "region_id", "week_idx", "units_sold"]].rename(columns={"units_sold": "naive"})
test = test.merge(ly.assign(week_idx=ly.week_idx + 52), on=["product_id", "region_id", "week_idx"], how="left")
test["model"] = np.nan
for pid, d in train.groupby("product_id"):
    b, _ = fit(d)
    t = test[test.product_id == pid]
    test.loc[t.index, "model"] = np.exp(design(t).values @ b.values)
test["actual_demand"] = test.units_sold   # note: stockout weeks understate true demand

def mape(a, f): return (np.abs(a - f) / a).mean()
cw = test.groupby(["category", "week_start"])[["actual_demand", "naive", "model"]].sum()
pw = test[test.on_promo == 1]
print(f"Category-week MAPE  naive: {mape(cw.actual_demand, cw.naive):.1%} | model: {mape(cw.actual_demand, cw.model):.1%}")
print(f"Promo-week MAPE     naive: {mape(pw.units_sold, pw.naive):.1%} | model: {mape(pw.units_sold, pw.model):.1%}")

tw = test.groupby("week_start")[["actual_demand", "naive", "model"]].sum()
fig, ax = plt.subplots(figsize=(6.4, 3.4))
ax.plot(tw.index, tw.actual_demand/1e3, "o-", color="#0B1F3A", label="Actual")
ax.plot(tw.index, tw.model/1e3, "o-", color="#14B8A6", label="Model forecast")
ax.plot(tw.index, tw.naive/1e3, "o--", color="#9CA3AF", label="Naive (same week last year)")
ax.set(ylabel="Units (000s)", title="Holdout: Oct–Dec 2025, all products")
ax.legend(fontsize=8, frameon=False); fig.autofmt_xdate(); plt.tight_layout()
plt.savefig("charts/forecast_holdout.png", dpi=150)

ab = pd.read_csv("data/ab_test_digital_coupon.csv")
s = ab.pivot_table(index=["store_id", "group"], columns="period", values=["sales", "gross_margin"]).reset_index()
s.columns = ["store_id", "group", "gm_pre", "gm_test", "sales_pre", "sales_test"]
s["sales_change"] = s.sales_test / s.sales_pre - 1    # each store vs. its own pre-period
s["gm_change"]    = s.gm_test / s.gm_pre - 1

def readout(col):
    t, c = s[s.group == "treatment"][col], s[s.group == "control"][col]
    diff = t.mean() - c.mean()
    se = np.sqrt(t.var()/len(t) + c.var()/len(c))
    p = stats.ttest_ind(t, c, equal_var=False).pvalue
    return diff, diff - 1.96*se, diff + 1.96*se, p

for col, name in [("sales_change", "Sales"), ("gm_change", "Gross margin $")]:
    d, lo, hi, p = readout(col)
    print(f"{name:15s} lift {d:+.2%}  (95% CI {lo:+.2%} to {hi:+.2%})  p = {p:.3f}")

w = ab.groupby(["week_start", "group"]).sales.mean().unstack()
w = w / w.iloc[:8].mean()
fig, ax = plt.subplots(figsize=(6.4, 3.3))
ax.plot(range(16), w.treatment, "o-", color="#14B8A6", label="Treatment (coupons)")
ax.plot(range(16), w.control, "o-", color="#9CA3AF", label="Control")
ax.axvline(7.5, color="#0B1F3A", ls="--", lw=1); ax.text(7.7, w.values.max(), "test starts", fontsize=8)
ax.set(xlabel="Week", ylabel="Sales index (pre = 1.00)", title="Digital coupon test: 20 vs. 20 stores")
ax.legend(fontsize=8, frameon=False); plt.tight_layout(); plt.savefig("charts/ab_test.png", dpi=150)