Bayesian Structural Time Series with CausalPy

47.1 Forecasting the counterfactual (40 min)

Concept. Synthetic control is one way to predict what R07 would have sold without the campaign. A time-series model is another: learn the relationship between R07's sales and things the campaign couldn't have affected (other regions' sales, trends, seasons) from the pre-period, then forecast the post-period. Google's CausalImpact (Brodersen et al., 2015) made this the standard tool in marketing analytics: a Bayesian structural time-series (BSTS) model with control series as predictors.

Being Bayesian gives two things synthetic control doesn't:

• a full posterior for the lift in every week, the cumulative lift, and anything computed from it, like profit;
• direct probability statements: "there's an 86% chance the campaign paid for itself".

We use CausalPy (by PyMC Labs), which runs these models on PyMC. Its InterruptedTimeSeries fits a Bayesian regression on the pre-period and forecasts the post-period. CausalPy also offers a full state-space version (StateSpaceTimeSeries, with a local level and trend that evolve over time); the regression version is faster and, with good control series, gives the same answers here.

PyMC samples with MCMC. The first run compiles the model, which takes a minute; later runs are faster. You'll see sampling messages: they're normal.

47.2 A time series without controls (40 min)

Create day47/causal_impact.py first: day47/its_no_controls.py imports its sampling settings.

# day47/causal_impact.py   (python -m day47.causal_impact)
"""
A Bayesian time-series model in the spirit of Google's CausalImpact (Brodersen et al., 2015), with CausalPy:
learn how R07's weekly sales move with other regions' sales BEFORE the campaign, forecast the post-period
from the other regions, and read the lift as actual minus forecast, with a posterior for every week.
"""
import causalpy as cp
import matplotlib.pyplot as plt
import polars as pl

from causalkit.regions import public, regions, truth, wide

SAMPLE = {"random_seed": 1, "progressbar": False, "draws": 1000, "tune": 1000, "chains": 2}


def causal_impact(df: pl.DataFrame, treated: str, start: int, controls: list[str]):
    """Fit CausalPy's InterruptedTimeSeries with control regions as regressors.
    Sales are centred on each region's pre-period mean and scaled to tens of $k, which makes sampling
    much faster; results are converted back to $k. Returns (result, scale)."""
    W = wide(df).to_pandas().set_index("week")
    Z = (W - W[W.index < start].mean()) / 10
    result = cp.InterruptedTimeSeries(Z, start, formula=f"{treated} ~ 1 + " + " + ".join(controls),
                                      model=cp.pymc_models.LinearRegression(sample_kwargs=SAMPLE))
    return result, 10.0


def summary(result, scale: float):
    tab = result.effect_summary().table
    return {k: {"mean": tab.loc[k, "mean"] * scale, "low": tab.loc[k, "hdi_lower"] * scale,
                "high": tab.loc[k, "hdi_upper"] * scale} for k in ["average", "cumulative"]}


if __name__ == "__main__":
    full = regions()
    df = public(full)
    t = truth(full)
    controls = ["R13", "R14", "R18", "R05", "R22"]          # Day 45's synthetic-control donors
    result, scale = causal_impact(df, "R07", 80, controls)
    s = summary(result, scale)
    print(f"Average lift per post week: {s['average']['mean']:.1f} $k  "
          f"[94% HDI {s['average']['low']:.1f} to {s['average']['high']:.1f}]   truth {t['avg_post_lift']:.1f}")
    print(f"Cumulative lift:            {s['cumulative']['mean']:.0f} $k  "
          f"[94% HDI {s['cumulative']['low']:.0f} to {s['cumulative']['high']:.0f}]   truth {t['cumulative_lift']:.0f}")
    fig, ax = result.plot()
    fig.savefig("figures/day47_causal_impact.png", dpi=110)
    plt.close(fig)
    print("Saved figures/day47_causal_impact.png")

Then create day47/its_no_controls.py:

# day47/its_no_controls.py   (python -m day47.its_no_controls)
"""
An interrupted time series with NO control regions: fit trend + seasonality to R07 before week 80 and
extend it. The model can't know that urban growth accelerated, so it credits the acceleration to the campaign.
"""
import causalpy as cp
import numpy as np

from causalkit.regions import public, regions, truth, wide
from day47.causal_impact import SAMPLE

full = regions()
W = wide(public(full)).to_pandas().set_index("week")
Z = (W - W[W.index < 80].mean()) / 10
Z["t"] = np.arange(len(Z)) / 52
Z["s1"], Z["c1"] = np.sin(2 * np.pi * Z.index / 52), np.cos(2 * np.pi * Z.index / 52)
result = cp.InterruptedTimeSeries(Z, 80, formula="R07 ~ 1 + t + s1 + c1",
                                  model=cp.pymc_models.LinearRegression(sample_kwargs=SAMPLE))
tab = result.effect_summary().table
print(f"Average lift per post week: {tab.loc['average', 'mean'] * 10:.1f} $k  "
      f"[94% HDI {tab.loc['average', 'hdi_lower'] * 10:.1f} to {tab.loc['average', 'hdi_upper'] * 10:.1f}]"
      f"   truth {truth(full)['avg_post_lift']:.1f}")

python -m day47.its_no_controls

Average lift per post week: 37.4 $k  [94% HDI 34.2 to 40.8]   truth 14.3

A confident and wrong answer: $37k a week, with a tight credible interval that excludes the truth. A trend fitted to R07's own past can't anticipate that urban growth accelerated after week 80, so the model credits the acceleration to the campaign. The credible interval only reflects uncertainty within the model; it can't warn you that the model is wrong. Without control series, an interrupted time series is a before-after comparison with a forecast attached.

(CausalPy reports a 94% highest-density interval, HDI, by default: the narrowest range containing 94% of the posterior.)

47.3 CausalImpact-style: control regions as predictors (65 min)

Now the control regions do the work: R07's sales are predicted from the five regions synthetic control chose on Day 45. Unlike synthetic control, the regression coefficients are unconstrained (any sign, any sum), and an intercept absorbs level differences.

One practical detail in causal_impact: each region's sales are centred on its pre-period mean and divided by 10 before fitting. On the original scale ($500k a week) the sampler needs many small steps and the run takes five minutes; centred and scaled, it takes twenty seconds. Results are converted back to $k.

python -m day47.causal_impact

Average lift per post week: 15.0 $k  [94% HDI 12.1 to 18.2]   truth 14.3
Cumulative lift:            375 $k  [94% HDI 302 to 454]   truth 357
Saved figures/day47_causal_impact.png

Within $1k a week of the truth, and the interval contains it. Open the figure: CausalPy's standard three panels (actual vs counterfactual, weekly impact, cumulative impact), in the scaled units (tens of $k, centred).

47.4 How many controls? (30 min)

Try it. Use all 29 other regions as predictors. Create day47/all_controls.py. Better or worse?

# day47/all_controls.py   (python -m day47.all_controls)
from causalkit.regions import public, regions, truth
from day47.causal_impact import causal_impact, summary

full = regions()
df = public(full)
controls = [r for r in df["region"].unique().sort().to_list() if r != "R07"]
s = summary(*causal_impact(df, "R07", 80, controls))
print(f"29 controls: {s['average']['mean']:.1f} $k/week [94% HDI {s['average']['low']:.1f}, {s['average']['high']:.1f}]"
      f"   truth {truth(full)['avg_post_lift']:.1f}")

29 controls: 13.1 $k/week [94% HDI 8.4, 18.3]   truth 14.3

47.5 Day 47 mini-project: did the campaign pay off? (65 min)

The campaign cost $100k. NorthMart keeps 30% of extra sales as margin. Every posterior draw of the cumulative lift gives a draw of profit = 0.30 × lift − 100. Create day47/posterior_roi.py:

# day47/posterior_roi.py   (python -m day47.posterior_roi)
"""
Mini-project: a posterior for the decision, not just for the lift.
The campaign cost $100k; NorthMart keeps 30% of extra sales as margin. Profit = 0.30 x cumulative lift - 100.
Every posterior draw of the lift gives a draw of profit, so we can state P(the campaign paid off).
"""
import matplotlib.pyplot as plt
import numpy as np

from causalkit.regions import public, regions, truth
from day47.causal_impact import causal_impact

MARGIN, COST = 0.30, 100.0
full = regions()
result, scale = causal_impact(public(full), "R07", 80, ["R13", "R14", "R18", "R05", "R22"])
lift = (result.post_impact.sum("obs_ind") * scale).values.ravel()        # cumulative lift, one per draw
profit = MARGIN * lift - COST
lo, hi = np.percentile(profit, [3, 97])
print(f"Profit: mean {profit.mean():+.1f} $k, 94% interval {lo:+.1f} to {hi:+.1f}")
print(f"P(profit > 0) = {(profit > 0).mean():.2f}")
print(f"P(return on spend above 20%) = {(profit > 0.2 * COST).mean():.2f}")
print(f"True profit: {MARGIN * truth(full)['cumulative_lift'] - COST:+.1f} $k")

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.hist(profit, bins=50, color="#2563eb", alpha=0.8)
ax.axvline(0, color="#dc2626", lw=1)
ax.set_xlabel("campaign profit ($k)")
ax.set_ylabel("posterior draws")
fig.tight_layout()
fig.savefig("figures/day47_profit_posterior.png", dpi=120)
print("Saved figures/day47_profit_posterior.png")

python -m day47.posterior_roi

Profit: mean +12.6 $k, 94% interval -9.5 to +34.5
P(profit > 0) = 0.86
P(return on spend above 20%) = 0.26
True profit: +7.2 $k
Saved figures/day47_profit_posterior.png

This is the output a decision-maker can use: "the campaign most likely made a small profit (86% probability), but it's unlikely to have returned 20% or more (26%)". Compare that with "the lift was significant", which says nothing about whether the money was well spent.

Write in day47/NOTES.md: if the marketing director asks whether to repeat the campaign at the same cost, what would you say, using these three numbers?

47.6 Save your work with Git (20 min)

git add .
git commit -m "Day 47: CausalPy interrupted time series, control regions, profit posterior"
git push

Day 47 checklist

• I can explain the forecasting approach to a counterfactual, and what CausalImpact adds
• I know that a time series without controls is a before-after comparison, whatever its interval says
• I can fit a CausalPy model with control series and read its HDIs
• I choose controls the treatment couldn't have affected
• I turn a posterior for the lift into a posterior for profit and probability statements