The Mathematics of Seasonality & Trend

A marketing-mix model decomposes a KPI into media-driven contributions plus a baseline — the sales the business would have made anyway. That baseline is rarely flat: it breathes with the calendar and drifts over the years. This page derives the two basis expansions the framework uses for those motions — Fourier harmonics for seasonality, and linear / piecewise / B-spline bases for trend — with the same math, the same assert-encoded identities, and the same charts as nbs/math/math_03_seasonality_trend.ipynb. If the model has no language for these motions, they don't disappear — they leak into the media coefficients and corrupt every ROI estimate.

This is notebook 03 of the math series. The full additive mean \( \mu_t \) that these terms slot into is assembled in Math 00 — The Generative Model; the priors named below are sampled for real in Math 04 — The Bayesian Model. The recurring idea on this page is one recipe used twice:

Definition: basis expansion

Build a fixed design matrix \( \Phi \in \mathbb{R}^{T \times M} \) of basis functions, then let the Bayesian model learn a coefficient vector \( \boldsymbol\theta \) under a regularizing prior: \( \text{component}_t = \Phi(t)\,\boldsymbol\theta \). The shape vocabulary lives in \( \Phi \) (this page); the amplitude is controlled by the priors. Get the basis right and the prior does the rest.

Motion	Basis	Framework function	Shape it draws
Seasonal (periodic)	Fourier harmonics	`create_fourier_features`	smooth, repeating
Trend (drift)	B-spline or piecewise-linear	`create_bspline_basis`, `create_piecewise_trend_matrix`	smooth / kinked

🧭 Before you start

Prerequisites: comfort with calculus & probability; the Workshop series recommended first. · Time: ~30–45 min.

You will learn:

How one basis-expansion recipe builds both seasonality (truncated Fourier series, zero-mean interleaved sin/cos columns) and trend (linear, piecewise, B-spline families).
Why Fourier order \( K \) is a bias–variance dial, why the piecewise trend changes slope rather than level, and why B-spline rows form a partition of unity.
Why baseline structure is an attribution control, not a fit knob — motions the model has no language for leak into the media coefficients.

Run it: nbs/math/math_03_seasonality_trend.ipynb on GitHub.

Part I — Seasonality via Fourier features

A seasonal pattern with period \( P \) (say \( P = 52 \) weeks for a yearly cycle) is a periodic function: \( g(t + P) = g(t) \). The Fourier theorem says any reasonable periodic function is a sum of harmonics — sines and cosines whose frequencies are integer multiples of the fundamental \( 1/P \):

\[ g(t) \;=\; \sum_{k=1}^{\infty} \Big[\, a_k \sin\!\Big(\tfrac{2\pi k\, t}{P}\Big) + b_k \cos\!\Big(\tfrac{2\pi k\, t}{P}\Big) \,\Big]. \]

Truncating at a finite order \( K \) gives a smooth approximation that sharpens as \( K \) grows. That truncation is exactly what create_fourier_features(t, period, order) builds: a design matrix \( \Phi \in \mathbb{R}^{T \times 2K} \) whose columns are the harmonics, interleaved as

\[ \Phi(t) \;=\; \big[\, \underbrace{\sin\tfrac{2\pi t}{P}}_{k=1},\; \underbrace{\cos\tfrac{2\pi t}{P}}_{k=1},\; \underbrace{\sin\tfrac{4\pi t}{P}}_{k=2},\; \underbrace{\cos\tfrac{4\pi t}{P}}_{k=2},\; \dots,\; \sin\tfrac{2\pi K t}{P},\; \cos\tfrac{2\pi K t}{P} \,\big]. \]

The column order is \( [\sin_1, \cos_1, \sin_2, \cos_2, \dots] \) — sine then cosine within each harmonic, not all sines followed by all cosines. In the model the seasonal coefficients carry a tight prior, \( [a_k, b_k] \sim \mathcal{N}(0,\, 0.3^2) \), and \( \text{seasonality}_t = \Phi(t)\,\boldsymbol\theta \) (seasonality = features @ coef). That shrinking prior is the regularization that tames the overfitting demonstrated below.

Why both sine and cosine?

A single sine \( \sin(2\pi k t / P) \) is pinned to a fixed phase (zero-crossing at \( t = 0 \)). Real seasonal peaks land wherever the calendar puts them. The identity

\[ A\sin(\omega t + \varphi) \;=\; \underbrace{A\cos\varphi}_{a}\,\sin(\omega t) + \underbrace{A\sin\varphi}_{b}\,\cos(\omega t) \]

shows that a sine/cosine pair spans every amplitude \( A \) and every phase \( \varphi \) at that frequency. Two columns per harmonic buy arbitrary phase for free — which is why the basis always comes in pairs and the matrix is \( (T, 2K) \), never \( (T, K) \).

Computed live in your browser — the same math as the notebook. The sine columns of the first three harmonics over one 52-week year: harmonic \( k \) completes exactly \( k \) cycles per period. Higher harmonics oscillate faster and let the basis carve out sharper features.

💡 Assert-encoded identities (verified by the notebook run)

The baked notebook verifies, column by column against np.sin/np.cos, that the order-3 matrix has shape (104, 6) with column \( 2(k{-}1) \) equal to the sine harmonic and \( 2(k{-}1){+}1 \) the cosine — and that the whole matrix is exactly periodic: \( \Phi(t{+}P) = \Phi(t) \), i.e. features[0] == features[52] holds to floating point (the docstring's promise, generalized to every row). Fourier columns are also zero-mean: they carry no constant, so the baseline level always needs a separate intercept.

Order is a bias–variance dial

What does order \( K \) buy? Take a deliberately non-sinusoidal target — a smooth wave plus a sharp recurring "holiday" spike — and reconstruct it by least squares at increasing order:

\[ \hat{\boldsymbol\theta}_K = \arg\min_{\boldsymbol\theta} \big\| \, g - \Phi_K\,\boldsymbol\theta \, \big\|_2^2 , \qquad \hat g_K = \Phi_K\,\hat{\boldsymbol\theta}_K . \]

On the noiseless target the notebook measures the residual \( \|g - \hat g_K\| \) falling monotonically: 3.408 at \( K{=}1 \), 1.473 at \( K{=}2 \), 0.858 at \( K{=}3 \), 0.616 at \( K{=}6 \) (order-6 correlation with the target 0.9927). More basis vectors can only help the in-sample fit. The question is whether the extra wiggles are signal or overfit — and noisy data answers it. Drive the dial yourself:

Fourier order \( K \) K = 4

about right

Computed live in your browser — the same math as the notebook (least squares on create_fourier_features-style harmonics + intercept, fit on two years of seeded noisy data, one year shown; the JS noise realization is seeded but differs from numpy's). Training residual – keeps falling as \( K \) rises; distance to the true signal – falls, bottoms out around \( K \approx 3\text{–}6 \), then climbs again as high orders chase noise between the data points. Measured in the notebook (seeded run): training residual \( K{=}1: 4.38 \ge K{=}4: 3.19 \ge K{=}24: 2.23 \), yet against the truth \( K{=}4 \) errs 1.11 vs \( K{=}24 \)'s 2.31 — the in-sample winner is the out-of-sample loser.

3.41 → 0.62

noiseless residual, order 1 → 6 (monotone, assert-checked)

1.11

error vs truth at K = 4 on noisy data (measured)

2.31

error vs truth at K = 24 — better in-sample, worse in fact

⚠️ The dial is a prior, not a search

The framework does not pick \( K \) by minimizing training error — that procedure is the red curve above. It pairs a modest default order with the shrinking \( \mathcal{N}(0, 0.3) \) prior on every \( [a_k, b_k] \): regularization, not raw flexibility, controls variance. What this cannot absorb is seasonality whose shape changes over the years — a static low-order basis fit to evolving seasonality scrambles the channel split while totals stay green (worst channel +104%); see Stress 02 — Time Structure.

Part II — The three trend families

Seasonality is periodic; trend is the non-periodic drift of the baseline — the slow climb or decline that does not repeat. The framework offers three families, ordered by flexibility (a fourth, a Gaussian-process / Hilbert-space trend, exists for fully nonparametric drift but is not exercised here). All trend bases are defined on a rescaled time axis \( t \in [0, 1] \) — a different convention from the integer weeks of Part I; feeding raw week numbers degenerates both bases.

Linear — the default. One slope: \( \text{trend}(t) = k + m\,t \) with \( m \sim \mathcal{N}(0, 0.5) \) (trend_slope). Maximal shrinkage, zero ability to bend. Right when the baseline is roughly stable; helpless at a structural break.
Piecewise-linear (Prophet-style) — a line whose slope changes at a handful of changepoints. Built by create_piecewise_trend_matrix. Derived next.
B-spline — a smooth, locally-supported polynomial basis. Built by create_bspline_basis. Derived after that.

Trend family 2 — Piecewise linear, done correctly

create_piecewise_trend_matrix(t, n_changepoints, changepoint_range=0.8) returns a pair \( (s, A) \):

\( s \in \mathbb{R}^{C} \) — the changepoint locations, evenly spaced inside the first changepoint_range of the axis (default 0.8, so no changepoints in the last 20% — a guard against overfitting the most recent, most influential observations). For \( C = 8 \) the notebook measures \[ s = [0.089,\ 0.178,\ 0.267,\ 0.356,\ 0.444,\ 0.533,\ 0.622,\ 0.711]. \]
\( A \in \{0,1\}^{T \times C} \) — a step indicator matrix, \( A_{t,j} = \mathbf{1}[\,t \ge s_j\,] \): column \( j \) switches on at \( s_j \) and stays on.

Computed live in your browser — the same math as the notebook. The design matrix \( A_{t,j} = \mathbf{1}[t \ge s_j] \) as a heatmap: each column turns on (dark) at its changepoint and stays on. The staircase of activation fronts is the piecewise basis. Assert-checked in the notebook: \( A \) is exactly binary, every column equals its indicator, and \( \max_j s_j \le 0.8 \).

The Prophet trend uses \( A \) not as a level shifter but as a switch on slope changes. Write a base growth rate \( k \) (the slope before any changepoint), an offset \( m \) (the level at \( t = 0 \)), and let each \( \delta_j \) be a change in slope that activates at \( s_j \). The cumulative slope at time \( t \) is \( k + (A\boldsymbol\delta)_t \), and the trend is

\[ \text{trend}(t) \;=\; \big(\,\underbrace{k}_{\text{base rate}} + (A\,\boldsymbol\delta)_t\,\big)\, t \;+\; \big(\,\underbrace{m}_{\text{offset}} + (A\,\boldsymbol\gamma)_t\,\big), \qquad \gamma_j \;=\; -\,s_j\,\delta_j . \]

Why the cumulative slope multiplies \( t \) — and what \( \boldsymbol\gamma \) is for. The whole point is that \( \delta_j \) changes the slope, not the level; that only happens if \( A\boldsymbol\delta \) enters as a coefficient on \( t \). But multiplying by \( t \) alone would make the trend jump at each changepoint — the new slope would suddenly apply to all the elapsed time before \( s_j \) too. The offset correction \( \gamma_j = -s_j \delta_j \) cancels that jump exactly. Expanding one active changepoint:

\[ \delta_j\, t + \gamma_j \;=\; \delta_j\,(t - s_j), \quad\text{which is zero at } t = s_j \;\;\Longrightarrow\;\; \text{trend}(t) \;=\; k\,t + m + \sum_{j=1}^{C} \delta_j\,(t - s_j)_{+}. \]

The pair (slope-change column \( A_{\cdot j}\, t \), level-correction column \( A_{\cdot j} \)) spans the cumulative ramp basis \( (t - s_j)_{+} \) — continuous, kinking only after its changepoint. The notebook asserts this identity to \( 10^{-12} \): the constructed trend equals the ramp form exactly, the largest single grid step is 0.0352 against a trend range of 1.18 (no level jumps), the slope after the last changepoint equals \( k + \sum_j \delta_j = 2.600 \), and the instantaneous slope jumps by exactly \( \delta_j \) across each active \( s_j \).

⚠️ Framework fix: slope changes, not level shifts

Using \( A\boldsymbol\delta \) as a bare additive term without the \( \cdot\,t \) — i.e. \( k + m\,t + A\boldsymbol\delta \) — makes it a sum of step functions: piecewise-constant level shifts, with the slope never changing. The framework's _build_piecewise_trend (in model/base.py) was found doing exactly that, and was fixed to the slope-change form above; the math notebook was rewritten (2026-06-02) to derive the corrected formula. The tell was the lone \( \boldsymbol\gamma \) continuity term sitting in the code: it is meaningless under level shifts (there is no slope jump to cancel) and only makes sense when slope changes were intended all along. In the model, sparse Laplace priors on \( \boldsymbol\delta \) (trend_delta) keep most changepoints at zero, so a slope bend is "spent" only where the data demand one.

Build a piecewise trend yourself

Fix the base rate \( k \) and offset \( m = 1 \), then drive three of the eight \( \delta_j \) sliders (the notebook's active set: \( s_1, s_3, s_6 \)). Watch the top panel stay continuous — no level jumps, only kinks — while the bottom panel shows the instantaneous slope \( k + (A\boldsymbol\delta)_t \) as a step function jumping by exactly \( \delta_j \) at each changepoint.

base rate \( k \) 0.6

\( \delta_1 \) at \( s_1{=}0.178 \) +5.0

\( \delta_3 \) at \( s_3{=}0.356 \) −9.0

\( \delta_6 \) at \( s_6{=}0.622 \) +6.0

Computed live in your browser — the same math as the notebook (\( \text{trend} = (k + A\boldsymbol\delta)\,t + (m + A\boldsymbol\gamma) \) with \( \gamma_j = -s_j\delta_j \), evaluated on the notebook's 160-point grid and its 8 measured changepoints). Final slope \( k + \sum\delta \) = 2.60. At the notebook's defaults (\( k{=}0.6, \delta = [+5, -9, +6] \)) this reproduces the baked figure, including the assert-checked final slope of 2.600. Set all \( \delta_j = 0 \) and the trend collapses to the dashed base line \( k\,t + m \).

Trend family 3 — B-splines as a partition of unity

create_bspline_basis(t, n_knots, degree=3) returns \( B \in \mathbb{R}^{T \times M} \) with \( M = \texttt{n\_knots} + \texttt{degree} + 1 \) columns (the notebook's \( n\_knots = 6 \), cubic, gives the measured shape (160, 10)). Each column is a B-spline: a piecewise polynomial that is non-zero only over a few adjacent knot intervals. That local support is what makes splines well-behaved — moving one coefficient nudges the curve locally instead of globally. Two structural properties are guaranteed by the construction and assert-encoded in the notebook:

Definition: partition of unity

At every \( t \), the basis functions sum to one: \[ \sum_{j=1}^{M} B_j(t) \;=\; 1 \qquad \text{for all } t \in [0,1]. \] The notebook asserts every row of \( B \) sums to 1 (measured max deviation 2.22 × 10⁻¹⁶ — machine epsilon), with every entry non-negative. A constant coefficient vector therefore reproduces a constant function: the basis never injects or removes "mass", it only redistributes it. This is the spline analogue of the normalized adstock kernel from Math 01 — Adstock (weights sum to 1).

Clamping. The boundary knots are repeated \( \text{degree}+1 \) times, pinning the first and last basis functions to 1 at the endpoints — assert-checked as \( B_{0,0} = 1 \) at \( t = 0 \) and \( B_{T-1,M-1} = 1 \) at \( t = 1 \). The fitted curve honors the endpoints cleanly rather than drifting off the edge of the data. The trend is \( B\boldsymbol\beta \) with \( \boldsymbol\beta \) given Normal priors whose scale is itself HalfNormal — hierarchical shrinkage that lets the data decide how wiggly the trend should be.

Computed live in your browser — the same construction as the notebook (clamped cubic B-splines via the Cox–de Boor recursion, 6 interior knots → 10 basis functions). Ten smooth local bumps, the first and last rising to 1 at the clamped boundaries — and their sum, the dark flat line pinned at 1 everywhere: the partition of unity made visible.

Head-to-head: three families, one kinked truth

Now the payoff. Synthesize one trending baseline with two structural breaks — a gentle climb (\( 2 + t \)), a sharp acceleration after \( t = 0.35 \) (\( +6\,(t-0.35)_+ \)), then a flattening-into-decline after \( t = 0.70 \) (\( -7\,(t-0.70)_+ \)) — plus noise. Fit each family by least squares:

Linear \( X = [\mathbf 1,\; t] \) — two parameters; must underfit the kinks.
Piecewise \( X = [\mathbf 1,\; t,\; A\!\cdot\! t,\; A] \) — the global line plus, per changepoint, a slope-change column \( A_{\cdot j}\,t \) and its level-correction column \( A_{\cdot j} \); together they span the ramp basis \( (t - s_j)_+ = A_{\cdot j}\,t - s_j A_{\cdot j} \), so least squares recovers the continuous slope-change trend exactly.
B-spline \( X = B \) — partition of unity already spans the constant, so no extra intercept column.

Computed live in your browser — the same math as the notebook (seeded JS noise realization; the notebook's seeded numpy run is the measured record). Measured correlations with the true kinked trend: linear 0.9613, piecewise 0.9983, B-spline 0.9990 — assert-checked that both flexible families exceed 0.97 and beat the linear fit. The linear trend slices straight through both breaks; the flexible bases bend with them. (The notebook does not assert the linear fit is good — by design it cannot be.)

0.961

linear trend × truth correlation on the kinked world (measured)

0.998

piecewise (Prophet-style) — catches both breaks

0.999

B-spline — same place, no event knowledge needed

Why this is about attribution, not fit

A 0.96 correlation sounds harmless — until you remember the trend is fit jointly with the media coefficients. Whatever baseline structure the trend family cannot represent does not vanish from the regression: the residual time-shape gets absorbed by whichever channel's spend happens to co-move with it. Misfit time structure is not a fit problem; it is an attribution problem wearing a good R².

⚠️ Measured consequence: TV at −63% with green diagnostics

Stress 02 — Time Structure builds exactly this trap: a structural trend break plus a reactive media ramp. Under the default linear trend, TV's contribution is misattributed by ~−63%, truth outside the 90% interval, r-hat clean. The same series also surfaced a second framework fix: the old piecewise changepoint prior, Laplace scale 0.05, was too stiff to pay for the break (fitted trend × truth correlation 0.33 — no better than the line). The default changepoint_prior_scale was changed 0.05 → 0.5 (2026-06-10), after which the trend tracks the shock-and-recovery at correlation 0.86 and median channel error drops below 20%. A flexible trend family is not enough — the prior on its flexibility is part of the spec.

The doctrine, stated once: the basis matrix decides what baseline shapes are representable; the prior decides how much of that flexibility gets used; and anything neither can express is silently reassigned to media. The full assembly — baseline plus adstocked, saturated media response, with every prior sampled and the decomposition graded — is Math 04 — The Bayesian Model.

Recap

What to take away

Seasonality is a truncated Fourier series: \( \Phi \in \mathbb{R}^{T\times 2K} \), columns interleaved \( [\sin_1, \cos_1, \sin_2, \cos_2, \dots] \) and zero-mean (the level needs its own intercept); each sin/cos pair buys arbitrary amplitude and phase.
Order \( K \) is a bias–variance dial: training residual falls monotonically with \( K \) (measured 4.38 → 3.19 → 2.23 for \( K = 1, 4, 24 \)), but error against the truth bottoms out near \( K \approx 4 \) (1.11) and worsens at \( K = 24 \) (2.31). The framework controls variance with the \( \mathcal{N}(0, 0.3) \) coefficient prior, not by chasing in-sample fit.
The piecewise trend changes slope, not level: \( \text{trend}(t) = (k + A\boldsymbol\delta)\,t + (m + A\boldsymbol\gamma) \) with \( \gamma_j = -s_j\delta_j \) is identically the continuous ramp sum \( k\,t + m + \sum_j \delta_j (t - s_j)_+ \) (assert-checked to \( 10^{-12} \)). Dropping the \( \cdot\,t \) gives level shifts — a real framework bug in _build_piecewise_trend, now fixed.
B-splines are a partition of unity: rows of \( B \) sum to 1 (measured max deviation 2.2 × 10⁻¹⁶), columns are non-negative with local support, boundaries clamped — the trend-side sibling of the normalized adstock kernel.
Baseline structure is an attribution control: on a kinked truth the flexible families hit correlation 0.998–0.999 vs linear's 0.961, and what the trend can't represent lands on whichever channel co-moves with it — TV −63% in Stress 02.

Check your understanding

Commit to an answer before opening each one.

Why does the Fourier seasonality basis always come in sine/cosine pairs, making the design matrix (T, 2K) rather than (T, K)?

A lone sine at a given frequency is pinned to a fixed phase, but real seasonal peaks land wherever the calendar puts them. The identity A sin(ωt + φ) = a sin(ωt) + b cos(ωt) shows that a sine/cosine pair spans every amplitude and every phase at that frequency, so two columns per harmonic buy arbitrary phase. See "Why both sine and cosine?" in Part I.

In the piecewise trend, why must the changepoint effects Aδ multiply t, and what does the correction γ_j = −s_jδ_j accomplish?

Each δ_j is meant to change the slope, which only happens if Aδ enters as a coefficient on t; used as a bare additive term it produces piecewise-constant level shifts instead — a real framework bug that was found and fixed. Multiplying by t alone would make the trend jump at each changepoint, because the new slope would retroactively apply to all elapsed time; γ_j cancels that jump exactly, collapsing each term to a continuous ramp that kinks only after its changepoint. The "Piecewise linear, done correctly" section derives this.

Why does this page insist that baseline structure is an attribution control rather than a fit knob?

The baseline is fit jointly with the media coefficients, so whatever time-shape the trend family cannot represent does not vanish — it is absorbed by whichever channel's spend happens to co-move with it. On the kinked truth a linear trend still correlates 0.961 with the truth, yet in Stress 02 the same misfit misattributed TV's contribution by about −63% with green diagnostics. The "Why this is about attribution, not fit" section states the doctrine: the basis decides what is representable, the prior decides how much flexibility is used, and anything neither can express lands on media.

💡 Run it yourself

This page mirrors nbs/math/math_03_seasonality_trend.ipynb, authored by nbs/builders/build_math_03_seasonality_trend.py (pure numpy/scipy — it bakes in seconds, no MCMC). Every plotted claim ends in an assert on the real framework functions (create_fourier_features, create_piecewise_trend_matrix, create_bspline_basis), so "the notebook executes clean" means "the math on this page is right". Series guide: nbs/README_math.md. Source: github.com/redam94/mmm-framework.