Marketing mix modeling in Python: a from-scratch tutorial with numpy and scipy

Marketing mix modeling in Python comes down to four steps: transform each channel’s spend for carryover (adstock) and diminishing returns (saturation), add controls for trend, seasonality and promotions, regress sales on all of it, and then read ROI and budget recommendations off the fitted model. You can do the whole thing with numpy, pandas and scipy in about 150 lines, which is what this tutorial does on a 104-week practice dataset.

When I built the earlier Vexpower MMM path, Python came after Excel for a reason: it’s much easier to trust a library like Meridian or PyMC-Marketing once you’ve built the same model by hand and seen where the judgment calls are. Every model result quoted on this page is printed by the downloadable notebook, and the code blocks below are the notebook’s cells, so you can run along and check you get the same answer. The practice data also has a known answer key, which almost never happens in real life, so at the end you’ll find out whether the method actually worked. If you need the background first, read what marketing mix modeling is.

What you need for marketing mix modeling in Python

You need Python 3.10 or later with numpy, pandas, scipy and matplotlib (pip install numpy pandas scipy matplotlib), and the practice CSV in the same folder as the notebook. If the file isn’t there, the first cell downloads it from this site. The dataset is invented for teaching: 104 weeks from January 2024 to December 2025, with revenue, spend on search, social and video, and a flag for promotion weeks. It’s deliberately small and clean, so it’s good for learning the mechanics and bad for ranking methods, and the data requirements page covers what a real dataset needs. It’s worth knowing up front that the data was generated from a straight linear equation with no noise, no carryover and no saturation, so the adstock and Hill steps below are there to teach the method: a correct model should find decay near zero and push saturation out to its bound, and that’s what happens.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize

DATA_URL = "https://vexpower.com/downloads/mmm/mmm-practice-synthetic.csv"
CHANNELS = ["search_spend", "social_spend", "video_spend"]

try:
    df = pd.read_csv("mmm-practice-synthetic.csv", parse_dates=["week"])
except FileNotFoundError:
    df = pd.read_csv(DATA_URL, parse_dates=["week"])

np.set_printoptions(suppress=True)
pd.set_option("display.width", 120)

Step 1: load and audit the data

The first job in any MMM is to check the data, because the most common way a model goes wrong is a gap, a duplicate week or a column in the wrong units, not a bad algorithm. Here that means the shape, missing values, the spacing between weeks, the range of each column and a first look at correlations.

print(df.shape)
print("Missing values:", int(df.isna().sum().sum()))
print("Week gaps:", df["week"].diff().dropna().value_counts().to_dict())
print(df["week"].min().date(), "to", df["week"].max().date())
print(df.drop(columns="week").describe().round(1).T[["mean", "min", "max"]])
print(df[["revenue"] + CHANNELS + ["promotion"]].corr().round(2)["revenue"])

(104, 6)
Missing values: 0
Week gaps: {Timedelta('7 days 00:00:00'): 103}
2024-01-01 to 2025-12-22
                 mean      min      max
revenue       25639.3  21720.0  30957.0
search_spend   2425.4   1565.0   3308.0
social_spend   1716.9   1000.0   2400.0
video_spend    1219.4    701.0   1700.0
promotion         0.1      0.0      1.0
revenue         1.00
search_spend    0.05
social_spend    0.19
video_spend    -0.10
promotion       0.22
Name: revenue, dtype: float64

The audit is clean: 104 consecutive weekly rows with no missing values. The correlations are the interesting part. Search has the biggest budget, yet its correlation with revenue is only 0.05, and video’s is negative. That doesn’t mean search does nothing; it means something else (here, a strong yearly seasonal cycle) is moving revenue by more than the media is, and a raw correlation can’t separate the two. This is the whole reason we build a regression with controls rather than eyeballing charts, and it’s a mild version of the confounding problem that makes MMM hard.

Step 2: geometric adstock and the Hill saturation curve

Two transformations do most of the work in a marketing mix model. Adstock captures carryover, the idea that some of this week’s advertising effect shows up next week and the week after. The geometric version keeps a fixed share of the previous week’s adstocked value, so a decay of 0.5 means half of last week’s effect survives into this week. The Hill saturation curve captures diminishing returns: it rises roughly in a straight line at low spend and flattens as spend grows, reaching half of its maximum at the half-saturation point.

def geometric_adstock(x, decay):
    """adstock[t] = x[t] + decay * adstock[t-1]"""
    out = np.zeros(len(x))
    carry = 0.0
    for t, value in enumerate(x):
        carry = value + decay * carry
        out[t] = carry
    return out


def hill(x, half_sat):
    """Share of the maximum effect reached at x. Equals 0.5 when x == half_sat."""
    return x / (x + half_sat)

I’ve fixed the Hill shape parameter at 1 to keep the search small, which makes it the same curve as Michaelis-Menten. The adstock here isn’t normalized, so a channel with high decay gets a bigger adstocked value; that only changes the scale of the coefficient, not the fit. If you want a vectorised version, scipy.signal.lfilter([1], [1, -decay], x) computes the same recursion without the loop (give it a float array; SciPy 1.18 raises an error on integer input).

Step 3: build the design matrix

The design matrix is every column the regression can use to explain revenue. The controls are an intercept (the baseline level of sales), a linear trend, a sine and cosine pair with a 52-week period to capture seasonality, and the promotion flag as a dummy variable. Each channel then gets one column: its spend, adstocked, then pushed through the Hill curve.

n = len(df)
t = np.arange(n)
controls = np.column_stack([
    np.ones(n),                   # intercept (baseline)
    t,                            # trend
    np.sin(2 * np.pi * t / 52),   # yearly seasonality
    np.cos(2 * np.pi * t / 52),
    df["promotion"].values,
])
control_names = ["intercept", "trend", "season_sin", "season_cos", "promotion"]
y = df["revenue"].values.astype(float)
spend = df[CHANNELS].values.astype(float)
mean_spend = spend.mean(axis=0)


def media_columns(params, spend=spend):
    decays, half_mults = params[:3], params[3:]
    cols = []
    for j in range(3):
        a = geometric_adstock(spend[:, j], decays[j])
        cols.append(hill(a, half_mults[j] * mean_spend[j]))
    return np.column_stack(cols)


def fit_ols(params, rows=slice(None)):
    X = np.column_stack([controls, media_columns(params)])
    beta, *_ = np.linalg.lstsq(X[rows], y[rows], rcond=None)
    return beta, X


def sse(params, rows=slice(None)):
    beta, X = fit_ols(params, rows)
    resid = y[rows] - X[rows] @ beta
    return float(resid @ resid)

Expressing half-saturation as a multiple of each channel’s average spend is a small trick that keeps the search well behaved. A multiple of 1 means the channel is already at half its maximum effect at average spend (heavily saturated); a multiple of 10 means average spend is far down the straight part of the curve. Once the six nonlinear parameters are fixed, the coefficients are just an ordinary least squares regression, which numpy solves in one line.

Step 4: search for decay and half-saturation

This is the part libraries hide from you. We have six nonlinear parameters (a decay and a half-saturation per channel) and, for any choice of them, OLS gives the best coefficients and a sum of squared errors. So we let scipy search the six parameters to minimize that error, starting from the best point on a coarse grid so the optimizer doesn’t get stuck somewhere silly.

bounds = [(0.0, 0.9)] * 3 + [(0.2, 10.0)] * 3

grid = [(d, k) for d in np.arange(0, 0.91, 0.1) for k in [0.5, 1, 2, 5, 10]]
start = min(([d] * 3 + [k] * 3 for d, k in grid), key=sse)

result = minimize(sse, x0=start, bounds=bounds, method="L-BFGS-B")
best = result.x
beta, X = fit_ols(best)
fitted = X @ beta

params_table = pd.DataFrame({
    "decay": best[:3].round(3),
    "half_sat_x_mean": best[3:].round(2),
    "half_sat_dollars": (best[3:] * mean_spend).round(0),
}, index=CHANNELS)
print(params_table)

r2 = 1 - ((y - fitted) ** 2).sum() / ((y - y.mean()) ** 2).sum()
mape = np.mean(np.abs((y - fitted) / y))
print(f"R-squared: {r2:.4f}   MAPE: {mape:.2%}")
print(pd.Series(beta, index=control_names + CHANNELS).round(1))

              decay  half_sat_x_mean  half_sat_dollars
search_spend  0.000             10.0           24254.0
social_spend  0.007             10.0           17169.0
video_spend   0.000             10.0           12194.0
R-squared: 1.0000   MAPE: 0.04%
intercept       19477.3
trend               0.0
season_sin       3000.8
season_cos          1.7
promotion        1803.8
search_spend    41131.8
social_spend    16488.7
video_spend      8799.9
dtype: float64

Two things in that output deserve attention. First, the search found almost no carryover: decay is 0 for search and video and 0.007 for social, which is nothing. Second, all three half-saturation multiples finished at 10.0, the top of the range I allowed. A parameter sitting on its bound is the model telling you it wanted to go further, and here it means the data shows no sign of diminishing returns at the spend levels we’ve seen. The fit itself is essentially perfect (R-squared 1.0000, MAPE 0.04 percent) because this synthetic data has no noise beyond rounding. On real data, an R-squared above 0.9 for weekly revenue is normal and anything near 1.0 should make you suspicious of overfitting.

The widget below runs the same model in your browser, on the same CSV, with the same transforms and the same OLS. Start from the notebook’s best fit and check the ROI column matches the notebook (1.54, 0.87 and 0.65), then try adding carryover or saturation and watch the fit get worse.

Step 5: decompose revenue into baseline and channels

A decomposition multiplies each column of the design matrix by its coefficient, so every week’s revenue is split into the pieces that add up to the model’s prediction. Adding the intercept, trend and seasonality together gives the baseline, the revenue you’d expect with no marketing at all.

contrib = pd.DataFrame(X * beta, columns=control_names + CHANNELS, index=df["week"])
decomp = pd.DataFrame({
    "baseline": contrib[["intercept", "trend", "season_sin", "season_cos"]].sum(axis=1),
    "promotion": contrib["promotion"],
    "search": contrib["search_spend"],
    "social": contrib["social_spend"],
    "video": contrib["video_spend"],
})
share = decomp.sum() / decomp.sum().sum()
print(share.round(3))

fig, ax = plt.subplots(figsize=(10, 4))
ax.stackplot(decomp.index, decomp.T.values, labels=decomp.columns, alpha=0.85)
ax.plot(df["week"], y, color="black", lw=1, label="actual revenue")
ax.set_ylim(15000, None)
ax.legend(loc="upper left", ncol=6, fontsize=8)
ax.set_title("Weekly revenue decomposition (baseline axis starts at 15,000)")
plt.tight_layout()
plt.show()

baseline     0.760
promotion    0.005
search       0.145
social       0.058
video        0.031
dtype: float64

According to the model, 76.0 percent of two years’ revenue is baseline, 14.5 percent comes from search, 5.8 percent from social, 3.1 percent from video and 0.5 percent from promotions. A large baseline share is typical; most of a business’s sales would happen without this week’s ads. The chart in the notebook stacks these pieces week by week, and it’s the first thing I’d show a stakeholder, because they can sanity-check it against what they remember (did the summer peak really happen, are the promotion spikes in the right weeks).

Step 6: calculate ROI and marginal ROI per channel

Average ROI is total modeled contribution divided by total spend. Marginal ROI is the return on the next dollar at today’s spend level, which is what you actually need for a budget decision, and the gap between the two is a direct read of how saturated a channel is. The marginal ROAS page explains why those two numbers are so often confused.

decays, half = best[:3], best[3:] * mean_spend
media_beta = beta[-3:]


def weekly_response(s, j):
    """Steady-state weekly contribution if channel j spends s every week."""
    a = s / (1 - decays[j])
    return media_beta[j] * hill(a, half[j])


def marginal_roi(s, j):
    a = s / (1 - decays[j])
    return media_beta[j] * half[j] / (a + half[j]) ** 2 / (1 - decays[j])


roi = pd.DataFrame({
    "spend": spend.sum(axis=0).round(0),
    "contribution": contrib[CHANNELS].sum().values.round(0),
}, index=CHANNELS)
roi["roi"] = (roi["contribution"] / roi["spend"]).round(2)
roi["marginal_roi"] = [round(marginal_roi(mean_spend[j], j), 2) for j in range(3)]
print(roi)

                 spend  contribution   roi  marginal_roi
search_spend  252245.0      387889.0  1.54           1.4
social_spend  178558.0      155805.0  0.87           0.8
video_spend   126822.0       82614.0  0.65           0.6

Search returns 1.54 in revenue per dollar on average and 1.40 at the margin, social 0.87 and 0.80, and video 0.65 and 0.60. Remember these are revenue returns, not profit, so with a 50 percent gross margin even search would be losing money on its own; the breakeven ROAS page covers that conversion. The fact that average ROI is about 10 percent higher than marginal ROI for every channel is a consequence of the Hill curve bending slightly at a half-saturation of 10, and we’ll come back to whether that bend is real.

Step 7: draw the response curves

A response curve shows the weekly revenue contribution a channel would produce at a steady weekly spend level. In steady state, geometric adstock settles at spend divided by one minus decay, so the curve is just the coefficient times the Hill function of that.

fig, ax = plt.subplots(figsize=(8, 4))
for j, name in enumerate(CHANNELS):
    s = np.linspace(0, 3 * mean_spend[j], 100)
    ax.plot(s, weekly_response(s, j), label=name)
    ax.scatter(mean_spend[j], weekly_response(mean_spend[j], j), zorder=3)
ax.set_xlabel("Weekly spend")
ax.set_ylabel("Weekly revenue contribution")
ax.set_title("Response curves (dots = average spend; beyond the highest observed spend is extrapolation)")
ax.legend()
plt.tight_layout()
plt.show()

On this data the curves are nearly straight lines across the range of spend we’ve observed, which is consistent with what the parameter search told us. The further right you read, the more you’re extrapolating: the maximum observed search spend was 3,308 a week, so anything the curve says about 6,000 a week is the Hill function’s opinion, not the data’s. That’s the most important caveat on any response curve, and it’s why the next step limits how far the budget can move.

Step 8: reallocate the same budget across channels

Budget optimization asks how to split a fixed total to get the most revenue. With the response curves in hand, it’s a constrained optimization: maximize the sum of the three channels’ weekly responses, keep the total equal to today’s average weekly spend, and limit each channel to plus or minus 30 percent of its current level so we don’t extrapolate wildly. The budget allocation page goes deeper on choosing those constraints.

total = mean_spend.sum()


def neg_total_response(s):
    return -sum(weekly_response(s[j], j) for j in range(3))


plan = minimize(
    neg_total_response,
    x0=mean_spend,
    bounds=[(0.7 * m, 1.3 * m) for m in mean_spend],
    constraints=[{"type": "eq", "fun": lambda s: s.sum() - total}],
    method="SLSQP",
).x

allocation = pd.DataFrame({
    "current": mean_spend.round(0),
    "proposed": plan.round(0),
    "change": (plan - mean_spend).round(0),
}, index=CHANNELS)
print(allocation)
uplift = -neg_total_response(plan) + neg_total_response(mean_spend)
print(f"Extra weekly revenue for the same budget: {uplift:,.0f}")

              current  proposed  change
search_spend   2425.0    3153.0   728.0
social_spend   1717.0    1355.0  -362.0
video_spend    1219.0     854.0  -366.0
Extra weekly revenue for the same budget: 474

The optimizer moves 728 a week into search, taking it to its 30 percent ceiling of 3,153, and funds that by cutting social by 362 and video by 366 (video hits its 30 percent floor). The model says that’s worth an extra 474 in weekly revenue for the same total spend of about 5,362 a week. Because the response curves are almost linear, the answer is close to “move as much as you’re allowed into the channel with the highest marginal ROI”, which is exactly what you’d expect. On real data with genuine saturation, the optimizer stops moving money once marginal returns are equal across channels, and that’s where it earns its keep.

Step 9: check the model on weeks it hasn’t seen

A model that explains the past perfectly can still be useless for planning. The simplest test is a holdout: refit everything on the first 78 weeks (18 months) and predict the last 26 weeks, which the model never saw. The model validation page covers the other checks I’d run.

train = slice(0, 78)
test = slice(78, n)
holdout_fit = minimize(sse, x0=start, args=(train,), bounds=bounds, method="L-BFGS-B")
beta_train, X_all = fit_ols(holdout_fit.x, train)
pred = X_all @ beta_train
test_mape = np.mean(np.abs((y[test] - pred[test]) / y[test]))
train_mape = np.mean(np.abs((y[train] - pred[train]) / y[train]))
print(f"Train MAPE: {train_mape:.2%}   Holdout MAPE: {test_mape:.2%}")

fig, ax = plt.subplots(figsize=(10, 3.5))
ax.plot(df["week"], y, color="black", lw=1, label="actual")
ax.plot(df["week"], pred, color="tab:green", label="model trained on first 78 weeks")
ax.axvline(df["week"].iloc[78], ls="--", color="gray")
ax.legend()
plt.tight_layout()
plt.show()

Train MAPE: 0.03%   Holdout MAPE: 0.09%

Training MAPE is 0.03 percent and holdout MAPE is 0.09 percent, so the model predicts unseen weeks almost as well as the weeks it was trained on. On real data you’d expect both to be in the low single digits, and the useful signal is the gap: if holdout error is several times the training error, the model has memorized noise. Holdout accuracy is necessary but not sufficient, though; a model can forecast well and still split the credit between channels wrongly, which is why the next step matters.

Did the model get the right answer?

This practice dataset was generated with known effects. Each dollar of search, social and video added 1.4, 0.8 and 0.6 of revenue in the same week, with no carryover and no saturation, on top of a baseline of 20,000, a seasonal swing of 3,000 and a promotion lift of 1,800. That makes it a fair test of whether the pipeline works.

truth = pd.Series([1.4, 0.8, 0.6], index=CHANNELS)
wide = minimize(sse, x0=start, bounds=[(0.0, 0.9)] * 3 + [(0.2, 100.0)] * 3, method="L-BFGS-B")
beta_wide, X_wide = fit_ols(wide.x)
roi_wide = (X_wide[:, -3:] * beta_wide[-3:]).sum(axis=0) / spend.sum(axis=0)
print(pd.DataFrame({
    "roi_bound_10": roi["roi"],
    "roi_bound_100": roi_wide.round(2),
    "true_roi": truth,
}))
print("Half-saturation multiples with the wider bound:", wide.x[3:].round(1))

              roi_bound_10  roi_bound_100  true_roi
search_spend          1.54           1.41       1.4
social_spend          0.87           0.81       0.8
video_spend           0.65           0.61       0.6
Half-saturation multiples with the wider bound: [100. 100. 100.]

The marginal ROIs from Step 6 (1.40, 0.80 and 0.60) match the truth exactly, and the promotion coefficient of 1,803.8 is close to 1,800. The average ROIs, however, were about 10 percent too high, and the reason is the bound I set on half-saturation. Forcing the Hill curve to bend a little at a multiple of 10 makes the early dollars look more productive than they are, and the intercept drops to 19,477 to compensate. Widening the bound to 100 lets the curve straighten out, and the average ROIs fall to 1.41, 0.81 and 0.61, within a cent of the truth. The lesson is general: the bounds and priors you put on saturation are assumptions, they change the ROI you report, and a parameter sitting on its bound is a flag you should always mention.

• The audit found no gaps, duplicates or obviously wrong units
• Seasonality and promotions are in the model as controls
• No parameter is sitting on its bound, or you've said so if one is
• Holdout error is close to training error
• Response curves are only read within (or just beyond) the observed spend range
• ROI is converted to profit before anyone makes a budget decision

Where to go next: PyMC-Marketing and Meridian

The from-scratch model gives you one best guess for each number. What it can’t do well is tell you how uncertain those numbers are, include what you already know from experiments, or pool information across regions. Those are the reasons to move to a Bayesian library, and the two I’d look at in Python are PyMC-Marketing and Google’s Meridian. The open-source MMM comparison has a picker if you’re not sure which.

The concepts map across directly. Your geometric_adstock becomes PyMC-Marketing’s GeometricAdstock (with an l_max for the maximum lag) or Meridian’s geometric adstock decay. Your hill becomes HillSaturation in PyMC-Marketing or Meridian’s Hill curve. The bounds you gave scipy become priors, which are a softer and more honest version of the same assumption: instead of “half-saturation can’t exceed 10 times average spend”, you say “it’s probably in this range, and the data can argue”. The OLS coefficients become posterior distributions, so ROI comes with an interval. And the answer-key step becomes calibration, where a lift test result is fed in as a prior (Meridian) or a likelihood term (PyMC-Marketing) so the model has to agree with it.

If you’d rather stay in another language, the same model is built step by step in R with base lm and in Excel with LINEST and Solver, so you can compare their numbers with this notebook’s. The Vexpower MMM course takes the ideas further on a messier invented brand, where the answer isn’t handed to you at the end.