Module 2 — Decomposition: trend, seasonality, noise

Before fitting any model, we look at the series. Not "look" in a poetic sense — a decomposition is a formal split of the series into three additive or multiplicative pieces, each of which has a name, a plot and consequences for the model choice we will make in later modules.

The pharmacy series has two obvious seasonalities: a weekly cycle (Saturday spikes, Sunday dips because most stores close), and a yearly cycle (winter respiratory-illness peaks, quiet Augusts). Any model that ignores either will lose to a baseline that merely copies the same day one year and one week ago.

The additive model

The simplest decomposition writes the observed value as a sum of three components:

$$y_t = T_t + S_t + R_t$$

where $T_t$ is the trend, $S_t$ the seasonal component (which we require to sum to zero over one period), and $R_t$ the residual. If a series has multiple seasonalities, $S_t$ splits into $S_t^{(7)}$ (weekly), $S_t^{(365)}$ (yearly) and so on.

An additive decomposition is appropriate when the amplitude of the seasonality is independent of the level of the trend. In our pharmacy data, a store selling 40 units per day on average has weekly swings of about $\pm 20$ units; a store selling 200 units per day also has weekly swings of $\pm 20$ units in absolute terms only if its behavior really is additive. In practice bigger stores swing more, which is where the multiplicative model comes in.

The multiplicative model

The multiplicative form writes:

$$y_t = T_t \cdot S_t \cdot R_t$$

with $S_t$ centered on 1. Now the seasonal effect scales with the level of the trend: a 20 % Saturday spike remains a 20 % spike as the store grows.

The practical trick to move between the two is a log transform. Taking $\log(y_t)$ turns the multiplicative decomposition into an additive one on the log scale:

$$\log(y_t) = \log(T_t) + \log(S_t) + \log(R_t)$$

For strictly positive series (unit counts, dollar revenue), we routinely fit models on $\log(1 + y_t)$ and back-transform predictions. Zeros stay in place, and multiplicative behavior is captured for free. This trick will resurface in modules 5 and 6.

STL: seasonal-trend decomposition using LOESS

Classical decomposition (moving averages) is fragile: it cannot handle changing seasonality, misses outliers badly, and assumes a single period. STL — Seasonal and Trend decomposition using LOESS — is the workhorse we will use, because it handles all three problems.

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import STL

df = pd.read_parquet("pharmacy_daily.parquet")
series = df.set_index("date")["units"].asfreq("D")

stl = STL(series, period=7, robust=True).fit()

fig, axes = plt.subplots(4, 1, figsize=(10, 8), sharex=True)
axes[0].plot(series); axes[0].set_title("observed")
axes[1].plot(stl.trend); axes[1].set_title("trend")
axes[2].plot(stl.seasonal); axes[2].set_title("weekly seasonal")
axes[3].plot(stl.resid); axes[3].set_title("residual")
plt.tight_layout()

STL fits a smooth trend, a periodic seasonal (here weekly), and returns the residual as what is left over. The robust=True flag downweights outliers when estimating the components — critical whenever a promotional day quadruples sales for one Saturday and would otherwise distort the seasonal curve for months.

Multiple seasonalities: weekly and yearly together

Our pharmacy series has both a 7-day and a 365-day period. Single-period STL cannot handle that; the successor MSTL (multiple seasonal-trend decomposition using LOESS) can.

from statsmodels.tsa.seasonal import MSTL

mstl = MSTL(series, periods=(7, 365)).fit()
weekly = mstl.seasonal["seasonal_7"]
yearly = mstl.seasonal["seasonal_365"]

The two seasonal components come out separately, and we can plot each one. On the pharmacy data, weekly amplitude sits around $\pm 15$ units, yearly amplitude around $\pm 25$ units with a peak in January and a trough in August. Removing both from the observed series leaves a residual whose behavior we can now think about cleanly.

Reading the residual

The residual is where most of the diagnostic value lives, and looking at it is worth ten minutes for every hour spent modeling.

Three things to check on the residual plot.

Zero mean. The residual should average to roughly zero. A drifting residual means the trend estimate is too smooth and missed a movement.

Constant variance. The residual should look like white noise of roughly constant amplitude. Widening variance is the classic sign that a multiplicative model was needed — try again with $\log y_t$.

No visible pattern. Any leftover pattern is a missed component: another seasonality, a change point, a promotional effect that should be modeled explicitly. In our data, decomposing without the yearly period leaves a residual that dips every August — an obvious signature the model must not miss.

resid = mstl.resid.dropna()
print("mean:", resid.mean().round(3))
print("std:", resid.std().round(3))

# Ljung-Box test: is the residual white noise?
from statsmodels.stats.diagnostic import acorr_ljungbox
print(acorr_ljungbox(resid, lags=[10, 20], return_df=True))

A Ljung-Box p-value below 0.05 says there is significant autocorrelation left in the residual: our decomposition has missed something.

Change points and level shifts

The pharmacy chain went through a rebrand in mid-2023 that lifted average sales by roughly 8 % overnight. STL will fit a smoother trend that partly absorbs the shift, but the seasonal and residual components will look bad for six months on each side of the event.

A quick check: fit STL, then look at the trend around the suspected date. A near-vertical step is a change point, and the honest answer is to model it explicitly — either an additive after_rebrand dummy variable (module 7), or splitting the training window to start after the rebrand if you have enough history.

Prophet, in module 6, includes automatic change-point detection; ARIMA and gradient boosting need you to declare them.

Decomposition is a diagnostic, not a model

STL and MSTL are excellent at showing you what is in your data, but by themselves they are not forecasting models. They produce the seasonal component only where the observed values exist; extending it into the future is a separate step (STLForecast in statsmodels does one version of this). Use decomposition to know what you are dealing with, then pick a model from the next modules.

Summary

• A time series can be decomposed as trend + seasonality + noise (additive) or trend × seasonality × noise (multiplicative); a log transform bridges the two.
• STL is the flexible workhorse for a single seasonality; MSTL handles multiple periods (weekly and yearly for our pharmacy).
• The residual is the diagnostic: zero mean, constant variance, no pattern. Anything else points to a missed component.
• Change points show up as vertical steps in the trend; ignore them and every downstream model will underperform for months around the event.

Next module: stationarity, differencing and autocorrelation — the language ARIMA speaks, and the tests that tell us if a series is ready for it.