C15 · Should Parkay run a promotion? — Bayesian Modelling in Practice

The brief

You work for the brand manager of Parkay margarine. Parkay sells sticks (its big seller) and tubs. The retailer has offered a slot for a temporary price reduction of 20% on Parkay sticks, funded entirely by Parkay. The brand manager asks:

"If we cut the stick price by 20% for a promotion period, how much extra share do Parkay sticks get? How much of that comes out of our own tubs, and how much do we steal from Blue Bonnet, the house brand and the rest? And after all that - does it make money? My team says promotions also win us loyal customers who keep buying Parkay afterwards, so count that too. I want to be at least 80% sure a promotion pays before I sign."

The economics, as the brand manager gives them (treat them as fixed):

Regular prices are each product's median shelf price in the panel (computed below). The promotion cuts the Parkay stick price by 20% for one shopping trip of every household (the promotion period); all other products stay at their regular prices.
Parkay's contribution margin at regular prices is $0.25 per lb on sticks and $0.45 per lb on tubs. Parkay funds the cut, so every promoted pound of sticks earns 20% of the regular shelf price less.
The region sells about 200,000 lb of margarine (all brands) in a promotion period. Treat every purchase as one pound (the panel has no quantities).
Carry-over counts: purchases on each household's next 6 shopping trips, at regular prices, that the promotion changed.

This is the pricing sequel to example E84 (market structure from the same panel). Do E84 first: it introduces the data, the pooled logit and its IIA property, and the idea that apparent loyalty can be heterogeneity in disguise. Here the question is a price decision, the model is a full random-coefficient logit, and the output is a profit distribution.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", price_sd=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import logging
import time

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from scipy import special

from pymc_challenges import Hints, data

RANDOM_SEED = 1991
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)  # several fits: no sampler banner per fit

h = Hints("C15")
h.tasks()

C15 - Should Parkay run a promotion?

task0 - Look at Parkay's prices (3 hints, 2 check(s))
task1 - The pooled logit (3 hints, 2 check(s))
task2 - Households have their own tastes (3 hints, 3 check(s))
task3 - Loyalty (3 hints, 2 check(s))
task4 - Simulate the promotion (3 hints, 5 check(s))
task5 - Does it pay? (3 hints, 3 check(s))

The data

One row per purchase, in time order within each household (the panel has no dates). choice is the product bought; price_<product> is the shelf price in US$ per lb of every product on that trip, bought or not. Products: Pk = Parkay, BB = Blue Bonnet, Fl = Fleischmann's, Hse = the store's house brand, Gen = generic, Imp = Imperial, SS = Shedd's Spread; _Stk = sticks, _Tub = tub.

The cell below adds integer codes (y, hh), each purchase's previous purchase by the same household (prev), and the evaluation split used throughout: a household's first purchase is only used as the "previous purchase" of its second (so that every model sees the same rows), and the last purchase of every household with at least three is held out for scoring. The regular prices of the brief are in REG.

data.describe("margarine")
df = data.load("margarine")
PRODUCTS = ["Pk_Stk", "BB_Stk", "Fl_Stk", "Hse_Stk", "Gen_Stk", "Imp_Stk",
            "SS_Tub", "Pk_Tub", "Fl_Tub", "Hse_Tub"]
J = len(PRODUCTS)
STK, TUB = PRODUCTS.index("Pk_Stk"), PRODUCTS.index("Pk_Tub")
PRICE_COLS = [f"price_{p}" for p in PRODUCTS]

df["y"] = df["choice"].map({p: j for j, p in enumerate(PRODUCTS)})
df["hh"] = pd.factorize(df["hhid"])[0]
H = df["hh"].max() + 1
df["prev"] = df.groupby("hh")["y"].shift(1)
df["T"] = df.groupby("hh")["trip"].transform("max")
rows = df[df["trip"] >= 2].copy()
rows["test"] = (rows["trip"] == rows["T"]) & (rows["T"] >= 3)
train, test = rows[~rows["test"]], rows[rows["test"]]

REG = df[PRICE_COLS].median().to_numpy()          # regular shelf prices (US$ per lb)
DEPTH = 0.20                                      # the promotion: 20% off Parkay sticks
MARGIN_STK, MARGIN_TUB = 0.25, 0.45               # Parkay's margin per lb at regular prices
VOLUME = 200_000                                  # category purchases (lb) in the promotion period
CARRY_TRIPS = 6                                   # later trips that count as carry-over

print(f"{len(df):,} purchases by {H} households; {len(train):,} training rows, "
      f"{len(test)} held-out last purchases")
print("regular prices:", dict(zip(PRODUCTS, REG.round(2))))
df.head()

margarine: 4,470 margarine purchases by 516 households, in time order within household (no dates): hhid, trip (1, 2, ... within household), choice (one of 10 products: Pk = Parkay, BB = Blue Bonnet, Fl = Fleischmann's, Hse = the store's house brand, Gen = generic, Imp = Imperial, SS = Shedd's Spread; _Stk = sticks, _Tub = tub), and price_<product>: the shelf price (US$ per pound) of every product on that shopping trip.
  source: Allenby & Rossi (1991), 'Quality perceptions and asymmetric switching between brands', Marketing Science 10:185-205: an A.C. Nielsen scanner panel, distributed as `margarine` in the R package bayesm (Rossi; GPL >= 2; the URL is its CRAN source tarball). CONVERTED from .rda by tools/build_e84_margarine.py - keep the cached file.
  url:    https://cran.r-project.org/src/contrib/Archive/bayesm/bayesm_3.1-6.tar.gz
4,470 purchases by 516 households; 3,498 training rows, 456 held-out last purchases
regular prices: {'Pk_Stk': np.float64(0.58), 'BB_Stk': np.float64(0.58), 'Fl_Stk': np.float64(0.99), 'Hse_Stk': np.float64(0.45), 'Gen_Stk': np.float64(0.33), 'Imp_Stk': np.float64(0.75), 'SS_Tub': np.float64(0.85), 'Pk_Tub': np.float64(1.09), 'Fl_Tub': np.float64(1.19), 'Hse_Tub': np.float64(0.59)}

	hhid	trip	choice	price_Pk_Stk	price_BB_Stk	price_Fl_Stk	price_Hse_Stk	price_Gen_Stk	price_Imp_Stk	price_SS_Tub	price_Pk_Tub	price_Fl_Tub	price_Hse_Tub	prev	T
0	2100016	1	Pk_Stk	0.66	0.67	1.09	0.57	0.36	0.93	0.85	1.09	1.19	0.33	NaN	7
1	2100016	2	Pk_Stk	0.63	0.67	0.99	0.57	0.36	1.03	0.85	1.09	1.19	0.37	0.0	7
2	2100016	3	Pk_Stk	0.29	0.50	0.99	0.57	0.36	0.69	0.79	1.09	1.19	0.59	0.0	7
3	2100016	4	Pk_Stk	0.62	0.61	0.99	0.57	0.36	0.75	0.85	1.09	1.19	0.59	0.0	7
4	2100016	5	Pk_Stk	0.50	0.58	0.99	0.45	0.33	0.72	0.85	1.07	1.19	0.59	0.0	7

Task 0 · Look at Parkay's prices

Deliver

The distribution of the Parkay stick price over trips. Is a 20% cut from the regular price something the data have seen, or an extrapolation? Where does most of the price variation come from?
The raw promotion response: Parkay stick share and Parkay tub share on trips grouped by the Parkay stick price, with uncertainty. What does it suggest about cannibalisation?
One reason this raw comparison is not the answer to the brief.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A model learns its price coefficient from price variation. Where is the Parkay stick price most of the time, and how far below that does it go? Then ask which way the raw share comparison could be biased: what else tends to happen on a deal trip?

2 Approach

Histogram df["price_Pk_Stk"] with the regular price and the promotion price (0.8 * REG[STK]) marked. Group the trips into price bands with pd.cut and compute, per band, the share of purchases that are Parkay stick and Parkay tub, with binomial standard errors sqrt(p (1 - p) / n).

3 Skeleton

pk = df["price_Pk_Stk"]
df["band"] = pd.cut(pk, [0, 0.32, 0.47, 0.55, 1.0])
band = df.groupby("band", observed=True).agg(
    n=("y", "size"),
    stk=("y", lambda s: (s == ___).mean()),
    tub=("y", lambda s: (s == ___).mean()))
share_promo = (df.loc[pk <= 0.47, "y"] == STK).mean()
share_regular = (df.loc[pk >= ___, "y"] == STK).mean()

Self-checks in the notebook (h.check("task0", …)):

share_promo — share of purchases that are Parkay stick on trips where it costs at most $0.47, about 0.71
share_regular — share of purchases that are Parkay stick on trips where it costs at least $0.55, about 0.26

Task 1 · The pooled logit, and what it says about the promotion

Fit a multinomial logit in which every household shares the same product intercepts and the same price coefficient (no previous-purchase term yet).

Deliver

The price coefficient, with a sanity check: its sign, and the own-price elasticity of Parkay sticks it implies at regular prices. Is that plausible for a grocery brand?
The promotion at regular prices: Parkay stick share before and after the 20% cut, the gain in percentage points, and where the gain comes from - for each other product, its share of the gain, compared with its share of the non-Parkay-stick market. Explain the pattern you find.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A purchase is one choice among 10 products, with the prices of all 10 on that trip. What happens to choice probabilities if you add the same number to every product's utility - so how many intercepts can you estimate? For the source of the gain, write down the logit's ratio of two products' probabilities and ask what a third product's price does to it.

2 Approach

Utility V[i, j] = alpha[j] + beta * price[i, j] with alpha fixed at 0 for Parkay sticks, and pm.Categorical("y", logit_p=V, observed=y) (the softmax is done for you). Elasticity of product k's share with respect to its own price in a logit: beta * p_k * (1 - s_k).

For the promotion, compute choice probabilities at the regular price vector REG and with REG[STK] multiplied by 0.8 (scipy.special.softmax over the product axis, for every posterior draw). The share of the gain coming from product j is -(s1[j] - s0[j]) / (s1[STK] - s0[STK]); compare it with s0[j] / (1 - s0[STK]).

3 Skeleton

with pm.Model(coords={"alt": PRODUCTS[1:]}) as pooled:
    a = pm.Normal("alpha", 0, 2, dims="alt")
    alpha = pt.concatenate([pt.zeros(1), a])
    beta = pm.Normal("beta", ___, ___)
    V = alpha + beta * X_tr                       # (purchases, 10)
    pm.Categorical("y", logit_p=___, observed=y_tr)
    idata_pool = pm.sample(random_seed=RANDOM_SEED)

post = az.extract(idata_pool, var_names=["alpha", "beta"])
A = np.concatenate([np.zeros((post.sizes["sample"], 1)),
                    post["alpha"].transpose("sample", ...).to_numpy()], axis=1)
b = post["beta"].to_numpy()
X1 = REG.copy(); X1[STK] *= 0.8
s0 = special.softmax(A + b[:, None] * REG, axis=-1)
s1 = special.softmax(A + b[:, None] * ___, axis=-1)

Self-checks in the notebook (h.check("task1", …)):

beta — posterior median of the pooled price coefficient, about -6.4 per $ per lb
gain_pp — Parkay stick share gain from the 20% cut at regular prices, about 18 percentage points

Task 2 · Households have their own tastes

Give every household its own product intercepts and its own price sensitivity, drawn from a population distribution that lets them be correlated (do households who like the premium brands care less about price?). Keep the price coefficient negative for everybody.

Deliver

The model written in its most direct form first. Diagnose what is wrong with the sampling, then fix it. Report what you changed and the diagnostics before and after.
The population median price coefficient and the household-to-household spread of price sensitivity. Which products' fans are least price-sensitive?
A comparison with the pooled logit on the held-out last purchases: the gain in log predictive density, with a standard error.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Each household has about eight purchases. How well can eight purchases pin down ten household coefficients - and what does that do to the joint posterior of those coefficients and their population spread? A price sensitivity that must be positive for everybody suggests a scale on which to put the population distribution.

2 Approach

Stack the 9 intercepts and the LOG price sensitivity into one 10-vector per household, theta_h ~ MvNormal(mu, Sigma), with pm.LKJCholeskyCov(..., compute_corr=True) for Sigma and b_h = -exp(theta_h[-1]). Index by household: V = A[hh] + b[hh, None] * X.

The direct form (pm.MvNormal("theta", mu, chol=chol, dims=("hh", "coef"))) is the centred parameterisation: look at r_hat, ESS and tree depth of the sds. The non-centred form is theta = mu + z @ chol.T with z ~ Normal(0, 1). For the held-out score, for every test purchase average its probability over posterior draws (log-mean-exp of the log probabilities) using that household's own theta, then sum; the standard error of a difference is sqrt(n) * std(pointwise difference).

3 Skeleton

with pm.Model(coords=coords) as hier:
    a = pm.Normal("alpha", 0, 2, dims="alt")
    mu_b = pm.Normal("mu_b", ___, 1.0)                 # mean LOG price sensitivity
    chol, _, _ = pm.LKJCholeskyCov("chol", n=10, eta=2.0,
                                   sd_dist=pm.HalfNormal.dist(___), compute_corr=True)
    mu = pt.concatenate([a, mu_b[None]])
    z = pm.Normal("z", 0, 1, dims=("hh", "coef"))
    theta = pm.Deterministic("theta", mu + z @ ___, dims=("hh", "coef"))
    A = pt.concatenate([pt.zeros((H, 1)), theta[:, :-1]], axis=1)
    b = -pt.exp(theta[:, -1])
    pm.Categorical("y", logit_p=A[hh_tr] + b[hh_tr, None] * X_tr, observed=y_tr)
    idata_hier = pm.sample(random_seed=RANDOM_SEED)

# held-out: th = draws of theta (S, H, 10)
V = A_draws[:, hh_te] + b_draws[:, hh_te, None] * X_te
lp = (V - special.logsumexp(V, -1, keepdims=True))[:, np.arange(len(y_te)), y_te]
pointwise = special.logsumexp(lp, 0) - np.log(___)

Self-checks in the notebook (h.check("task2", …)):

median_beta — price coefficient of the median household, -exp(posterior median of mu_b), about -9 per $ per lb
price_sd — posterior median of the sd of LOG price sensitivity across households, about 0.40
elpd_gain — held-out log predictive density of the hierarchical model minus the pooled one (no lag in either), about 260 nats

Task 3 · Loyalty: does buying Parkay make you buy Parkay?

The brand manager's team believes promotions "win loyal customers". In a choice model that is state dependence: buying a product last time raises the chance of buying it again, beyond a household's fixed tastes.

Deliver

Add the previous purchase to the pooled and to the hierarchical model. How big is the loyalty effect in each, as an odds ratio, and why do they disagree?
What adding it does to the estimated heterogeneity (the sds of Task 2), and what the held-out purchases say about the four models.
Which loyalty estimate should go into the promotion's carry-over, and why.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

If households simply have stable tastes, which products will follow which in their purchase histories - even with no habit at all? Which of the two models can tell "same product because same household" from "same product because bought last time"?

2 Approach

Add gamma * last to the utilities, where last is a (purchases, 10) one-hot matrix of the previous purchase (np.eye(J)[train["prev"].astype(int)]). Compare gamma between the pooled and the hierarchical model as exp(gamma), the odds ratio. Compare the posterior means of chol_stds with and without the lag, and add both lag models to the held-out comparison.

3 Skeleton

last_tr = np.eye(J)[train["prev"].astype(int).to_numpy()]
# inside each model:
gamma = pm.Normal("gamma", 0, 2)
V = ___ + gamma * last_tr
# afterwards
gamma_pooled = np.median(az.extract(idata_pool_lag, var_names="gamma"))
gamma_hier = np.median(az.extract(idata_hier_lag, var_names=___))

Self-checks in the notebook (h.check("task3", …)):

gamma_pooled — posterior median of the previous-purchase effect in the pooled logit, about 1.8 (odds x6)
gamma_hier — posterior median of the previous-purchase effect in the hierarchical logit, about 0.4 (odds x1.5)

Task 4 · Simulate the promotion

Use the households' current state: each household's most recent purchase in the panel is its "previous purchase" when the promotion trip comes.

Deliver, for the pooled-with-lag and the hierarchical-with-lag model:

The Parkay stick share gain on the promotion trip (percentage points) and the fraction of the gain that is cannibalised from Parkay tub, compared with what IIA would have said.
Which rivals lose the most, relative to their size, under the hierarchical model, and why. (Call a product's share of the gain divided by its IIA share its "lift".)
The carry-over: the extra Parkay (stick + tub) purchases per household on each of the next 6 trips at regular prices. Plot the decay for both models.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

With a previous-purchase effect, a promoted purchase changes the household's state for the next trip, which changes the trip after, and so on. For one household and one posterior draw that is a Markov chain over the 10 products. What is its transition matrix at regular prices?

2 Approach

Promotion trip: per draw and household, probabilities P0 (regular prices) and P1 (stick cut), using each household's latest purchase as its state. Transition matrix Q[s, h, k, j] = softmax_j(A[s, h, j] + b[s, h] * REG[j] + gamma[s] * (j == k)). Then propagate both: q <- einsum("shk,shkj->shj", q, Q) six times, and average q1 - q0 over households. Use a few hundred draws: Q is draws x households x 10 x 10.

Product j's share of the gain is -(s1[j] - s0[j]) / (s1[STK] - s0[STK]); IIA would give s0[j] / (1 - s0[STK]). Their ratio is the product's "lift"; cannibalisation is the Parkay tub entry.

3 Skeleton

state = np.eye(J)[df.groupby("hh")["y"].last().sort_index().to_numpy()]   # (H, J)
V0 = A + b[..., None] * REG + g[:, None, None] * state                    # (S, H, J)
V1 = A + b[..., None] * X1 + g[:, None, None] * ___
P0, P1 = special.softmax(V0, -1), special.softmax(V1, -1)
Q = special.softmax(A[:, :, None, :] + b[:, :, None, None] * REG
                    + g[:, None, None, None] * np.eye(J), axis=-1)       # (S, H, from, to)
q0, q1, carry = P0, P1, []
for _ in range(6):
    q0 = np.einsum("shk,shkj->shj", q0, Q)
    q1 = np.einsum("shk,shkj->shj", ___, Q)
    carry.append((q1 - q0).mean(1))

Self-checks in the notebook (h.check("task4", …)):

gain_pp — Parkay stick share gain on the promotion trip under the hierarchical + lag model, about 18.4 points
cannibalised — fraction of the stick gain taken from Parkay tub (hierarchical + lag), about 0.07
bb_lift — Blue Bonnet stick's share of the gain divided by its IIA share-proportional share (hierarchical + lag), about 1.4
carry_pooled — extra Parkay (stick + tub) purchases per household over the next 6 trips, pooled + lag model, about 0.14
carry_hier — the same under the hierarchical + lag model, about 0.02

Task 5 · Does it pay?

Deliver

The posterior distribution of Parkay's incremental profit from the 20% promotion over the whole period (promotion trip plus carry-over, sticks and tubs, at the brief's margins and volume), and the probability that it pays. Show what the pooled-with-lag model would have told the brand manager.
The break-even stick margin: how large would Parkay's stick margin have to be for the 20% promotion to break even?
Is there a discount depth that the hierarchical model says pays? Plot expected incremental profit against depth (5% to 30%) with uncertainty, for both models.
A recommendation in plain words, including what this analysis cannot see.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Write the profit change per category purchase as three pieces: what the extra stick pounds earn, what the cut gives away on EVERY promoted pound (including those that would have been sold anyway), and what the lost tub pounds cost. Then add the carry-over pounds at regular margins. Do it per posterior draw - the probability that it pays is a fraction of draws.

2 Approach

Per draw: d_stk and d_tub = change in share on the promotion trip plus the summed carry-over. Profit = VOLUME * (MARGIN_STK * d_stk + MARGIN_TUB * d_tub - 0.2 * REG[STK] * s1[STK]). Setting it to zero and solving for the stick margin gives the break-even margin. For the depth curve, rerun the simulation for depths 0.05 ... 0.30.

3 Skeleton

def profit(r, depth=0.2):
    d_stk = r["s1"][:, STK] - r["s0"][:, STK] + r["carry"][:, :, STK].sum(1)
    d_tub = r["s1"][:, TUB] - r["s0"][:, TUB] + r["carry"][:, :, TUB].sum(1)
    cost = depth * REG[STK] * ___
    return VOLUME * (MARGIN_STK * d_stk + MARGIN_TUB * d_tub - cost)

breakeven = (0.2 * REG[STK] * r["s1"][:, STK] - MARGIN_TUB * d_tub) / ___
p_pays = (profit(r) > 0).mean()

Self-checks in the notebook (h.check("task5", …)):

profit — posterior median incremental profit of the 20% promotion (hierarchical + lag), about -3300 US$
p_pays — posterior probability that the 20% promotion pays (hierarchical + lag), essentially 0
breakeven_margin — posterior median break-even stick margin (hierarchical + lag), about $0.33 per lb

Going further

Non-linear price response. Utility linear in price learns its slope from the half-price deals and applies it to a 20% cut. Let the response bend (log price, a separate coefficient for the discount from the regular price, or a deal-depth spline) and see whether the shallow-cut profit changes.
Demographics in the hierarchy. margarine_demos has income, family size and education for every household. Put them in the population mean of the household coefficients: are poorer or larger households more price-sensitive, and could Parkay target the promotion (a coupon) at them?
Reference prices. Households who have seen Parkay at $0.29 may treat $0.46 as "not a deal". Add a reference-price (loss-aversion) term built from each household's past prices.
Who should be promoted? Run the same calculation for a 20% cut on Parkay tubs or Blue Bonnet sticks. Which product's promotion steals most from rivals per dollar of margin given away?

h.progress()

C15: 0/18 hints used, 17/17 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Look at Parkay's prices

#Task 1 · The pooled logit, and what it says about the promotion

#Task 2 · Households have their own tastes

#Task 3 · Loyalty: does buying Parkay make you buy Parkay?

#Task 4 · Simulate the promotion

#Task 5 · Does it pay?

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Look at Parkay's prices

Task 1 · The pooled logit, and what it says about the promotion

Task 2 · Households have their own tastes

Task 3 · Loyalty: does buying Parkay make you buy Parkay?

Task 4 · Simulate the promotion

Task 5 · Does it pay?

Going further