Time-Series Decomposition

"I am the trend, and I have been quietly taking credit for work the season did. That winter bump everyone admires as my steady climb? Half of it was a repeating swing I absorbed when nobody was watching the windows. Decompose me carefully or I will keep pocketing the season's applause."

A Trend Quietly Stealing Credit From the Season

The previous section, Section 3.4, taught you to recognize the building blocks of a series: a level, a trend, one or more seasonal cycles, and irregular fluctuation. Recognition is diagnosis; this section is treatment. We take the additive and multiplicative component models written there and turn them into estimators that actually return three numbered sequences from one observed sequence. The thread running through the whole section is a single equation and a single question. The equation is that an observed series $y_t$ decomposes into trend $T_t$, seasonal $S_t$, and remainder $R_t$. The question, asked of every method, is the same: how flexible is the seasonal pattern allowed to be, and how much does a single outlier get to distort the answer? Classical decomposition gives a rigid, fragile answer; STL gives a flexible, robust one; and the modern variants push robustness and multiple calendars further still.

1. The Decomposition Idea Beginner

Write the observed series as a combination of unobserved components. In the additive model the parts add,

where $T_t$ is a smoothly varying trend-cycle, $S_t$ is a seasonal component that repeats with period $m$ (so $m = 12$ for monthly data with a yearly cycle, $m = 7$ for daily data with a weekly cycle), and $R_t$ is the remainder, everything the first two components could not explain. The additive form is right when the size of the seasonal swing does not depend on the level of the series. When the seasonal swing grows in proportion to the level, as retail sales do every December, the multiplicative model fits instead,

and a logarithm converts one into the other, since $\log y_t = \log T_t + \log S_t + \log R_t$ turns a multiplicative decomposition into an additive one on the log scale. This is why "take logs first" is the standard reflex for series whose variance grows with level; it is the log-to-additive identity of Section 3.4, itself the Box-Cox/variance-stabilizing idea introduced in Section 3.2, applied so that an additive estimator can do multiplicative work. Throughout this section we use a monthly retail-style series with a strong yearly cycle as the running example, consistent with the finance-and-retail flavor of Part II. Figure 3.5.2 pictures the whole task in one image: patiently unbraiding one tangled series into separable strands.

What is decomposition for? Three uses recur. The first is understanding: a deseasonalized unemployment or sales series strips away the predictable calendar swing so that a genuine turning point becomes visible, which is exactly why statistical agencies publish "seasonally adjusted" numbers. The second is adjustment: many downstream models, and many human decisions, want $y_t - S_t$ (or $y_t / S_t$) rather than the raw series. The third is forecasting: if you can forecast each component separately and recombine them, you have a baseline forecaster that is transparent, fast, and surprisingly hard to beat, which is the construction of subsection five. Table 3.5.1 fixes the vocabulary before we estimate anything.

2. Classical Decomposition Beginner

The oldest recipe is mechanical and instructive. It proceeds in three steps. First, estimate the trend-cycle with a centered moving average of length equal to the seasonal period $m$, because averaging over exactly one full cycle cancels the seasonal pattern. When $m$ is even (as for monthly data, $m = 12$) a plain moving average is not centered on an integer index, so the standard fix is the $2 \times m$ moving average: average $m$ consecutive points, then average two adjacent such averages. Concretely, for even $m$,

which weights the two end points by one half so the window spans exactly $m$ intervals and stays centered. Second, detrend by subtracting the trend (additive) or dividing by it (multiplicative), leaving a detrended series that holds the seasonal pattern plus remainder. Third, estimate the seasonal component by averaging the detrended values within each season: all the Januaries together, all the Februaries together, and so on. These $m$ seasonal averages are then centered to sum to zero (additive) or to average one (multiplicative), and tiled across the whole series to form $\hat S_t$. The remainder is whatever is left, $\hat R_t = y_t - \hat T_t - \hat S_t$.

The method has three limitations that motivate everything after it. The seasonal pattern is forced to be identical in every cycle, because a single average per season is reused for all years; a December effect that strengthens over a decade cannot be represented. The trend is undefined at the two ends, because a centered window of width $m$ loses $m/2$ points at each end, so the first and last half-year of a monthly series have no trend estimate at all. And it has no robustness: one freak observation contaminates the moving average around it and the seasonal average for its season. Figure 3.5.1 shows the three steps on the running series.

Code 3.5.1 implements the full classical additive decomposition from scratch, so that the moving average, the per-season averaging, and the centering are all explicit rather than hidden behind a call.

import numpy as np

def classical_decompose(y, m):
    """Additive classical decomposition: trend by 2xm centered MA,
    seasonal by per-season averaging, remainder by subtraction."""
    n = len(y)
    # --- Step 1: trend via centered moving average of length m ---
    trend = np.full(n, np.nan)
    half = m // 2
    if m % 2 == 0:                                   # even m: 2xm MA (half weights at ends)
        w = np.concatenate(([0.5], np.ones(m - 1), [0.5])) / m
        for t in range(half, n - half):
            trend[t] = np.dot(w, y[t - half: t + half + 1])
    else:                                            # odd m: plain centered MA
        w = np.ones(m) / m
        for t in range(half, n - half):
            trend[t] = np.dot(w, y[t - half: t + half + 1])
    # --- Step 2: detrend ---
    detrended = y - trend
    # --- Step 3: seasonal = average of detrended values within each season ---
    seasonal_means = np.array([np.nanmean(detrended[s::m]) for s in range(m)])
    seasonal_means -= seasonal_means.mean()          # center so the season sums to zero
    seasonal = np.tile(seasonal_means, n // m + 1)[:n]
    remainder = y - trend - seasonal
    return trend, seasonal, remainder

# Build a synthetic monthly series: linear trend + yearly season + noise.
rng = np.random.default_rng(0)
t = np.arange(72)
true_trend = 100 + 0.5 * t
true_season = 8 * np.sin(2 * np.pi * t / 12)
y = true_trend + true_season + rng.normal(0, 1.5, size=72)

T, S, R = classical_decompose(y, m=12)
print("first valid trend index:", np.argmax(~np.isnan(T)))
print("seasonal pattern (one cycle):", np.round(S[:12], 2))
print("remainder std (where defined):", round(np.nanstd(R), 3))

first valid trend index: 6
seasonal pattern (one cycle): [ 0.13  4.06  6.97  8.02  6.84  4.   -0.07 -4.08 -6.93 -8.   -6.96 -4.05]
remainder std (where defined): 1.43

3. STL: Seasonal-Trend Decomposition Using Loess Intermediate

STL, introduced by Cleveland and colleagues in 1990, replaces every rigid step of the classical method with a loess smoother and so removes all three limitations at once. Loess (locally estimated scatterplot smoothing) fits a low-degree polynomial in a moving neighborhood, weighting nearby points more than distant ones; sliding that local fit along the series produces a smooth curve whose flexibility you control through the neighborhood width. STL uses loess in two nested loops. The inner loop alternates between estimating the seasonal component and the trend: it detrends, then smooths the values for each season position across cycles with loess (so the December effect may drift slowly from year to year rather than being a single fixed number), removes any low-frequency leakage from that seasonal estimate, subtracts the seasonal from the data, and smooths the result with loess to get the trend. Because the seasonal sub-series is itself smoothed over time, STL gives evolving seasonality, the single feature classical decomposition cannot provide.

The outer loop adds robustness. After an inner-loop pass, STL computes the remainder and assigns each observation a robustness weight that shrinks toward zero for points with large remainders, using the bisquare function

so an observation whose remainder exceeds six times the median absolute remainder gets weight zero and is effectively ignored on the next inner-loop pass. Repeating the inner loop with these weights makes the trend and seasonal estimates resistant to outliers, a spike no longer drags the smoother toward it. Two parameters dominate in practice. The seasonal window (called seasonal or $n_{(s)}$) sets how many cycles the seasonal loess looks across: a large window means a nearly constant seasonal pattern (close to classical), a small window lets it change quickly. The trend window ($n_{(t)}$) sets the trend smoother's span; larger means a smoother, stiffer trend. A useful default is to make the seasonal window large unless you have evidence the pattern truly evolves, and to set the trend window to roughly $1.5 m / (1 - 1.5/n_{(s)})$ rounded to the next odd number, the rule the original authors recommend, which most libraries apply automatically.

import numpy as np
from statsmodels.tsa.seasonal import STL

# Reuse the synthetic monthly series y from Code 3.5.1, but make the
# seasonal AMPLITUDE grow over time so classical decomposition would fail.
rng = np.random.default_rng(1)
t = np.arange(120)
trend = 100 + 0.4 * t
growing_season = (1 + 0.012 * t) * 8 * np.sin(2 * np.pi * t / 12)  # amplitude rises
y = trend + growing_season + rng.normal(0, 1.5, size=120)

# robust=True turns on the outer-loop bisquare robustness weights.
stl = STL(y, period=12, seasonal=13, robust=True).fit()
T, S, R = stl.trend, stl.seasonal, stl.resid

# The seasonal pattern is NOT identical each year: compare year 1 vs year 9.
amp_year1 = S[0:12].max() - S[0:12].min()
amp_year9 = S[96:108].max() - S[96:108].min()
print(f"seasonal swing year 1: {amp_year1:5.2f}")
print(f"seasonal swing year 9: {amp_year9:5.2f}   (evolving seasonality)")
print(f"remainder std: {R.std():.3f}")

seasonal swing year 1:  16.41
seasonal swing year 9:  21.83   (evolving seasonality)
remainder std: 1.52

4. Robust and Multiple-Seasonality Variants Advanced

Two extensions matter for real data. The first is robustness, which we have already met as STL's outer loop. Setting robust=True is the right default whenever the series contains outliers, level shifts, or promotional spikes: the bisquare weights down-weight contaminated points so a single Black-Friday spike does not bleed into the trend or distort the December seasonal for every other year. The cost is a few extra loop passes and a slight loss of efficiency when the data are in fact clean, which is almost always worth paying. The second extension is multiple seasonality. High-frequency series rarely have a single cycle: hourly electricity demand has a daily cycle ($m_1 = 24$), a weekly cycle ($m_2 = 168$), and a yearly cycle ($m_3 = 8766$) all at once, and a single seasonal component cannot hold three periods. MSTL (Multiple STL) solves this by applying STL iteratively, extracting one seasonal component per period and passing the deseasonalized remainder to the next period, so the additive model becomes

with one seasonal term $S_t^{(k)}$ per period $m_k$. MSTL, proposed by Bandara, Hyndman, and Bergmeir in 2021 and now the standard tool for several calendars, cycles through the periods from shortest to longest, refining each seasonal estimate over a few iterations. Table 3.5.2 is the decision guide for choosing among the three decomposers.

import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import MSTL

# Hourly series with TWO seasonalities: daily (24) and weekly (168).
rng = np.random.default_rng(2)
h = np.arange(24 * 7 * 8)                          # eight weeks of hourly data
daily  = 10 * np.sin(2 * np.pi * h / 24)           # within-day cycle
weekly = 6 * np.sin(2 * np.pi * h / 168)           # within-week cycle
level  = 50 + 0.002 * h
# Wrap as a pandas Series so MSTL returns labeled seasonal columns.
y = pd.Series(level + daily + weekly + rng.normal(0, 1.0, size=h.size))

# One call, a list of periods, two separate seasonal components returned.
res = MSTL(y, periods=(24, 168)).fit()
S = res.seasonal                                   # DataFrame: one column per period
print("seasonal components:", list(S.columns))
# np.ptp is the peak-to-peak (max - min) function; the ndarray .ptp() method
# was removed in NumPy 2.0, so call the function form on the column values.
print("recovered daily swing :", round(np.ptp(S['seasonal_24'].values), 2))
print("recovered weekly swing:", round(np.ptp(S['seasonal_168'].values), 2))

seasonal components: ['seasonal_24', 'seasonal_168']
recovered daily swing : 20.07
recovered weekly swing: 11.83

5. Using Decomposition in Practice Intermediate

Decomposition earns its keep in three concrete jobs. The first is seasonal adjustment: the deseasonalized series is $y_t - \hat S_t$ in the additive case (or $y_t / \hat S_t$ multiplicatively), and it is what you plot when you want to see the trend without the calendar shouting over it. The second is residual diagnostics. The remainder $\hat R_t$ is a direct verdict on the decomposition: if trend and seasonal captured all the structure, $\hat R_t$ should look like white noise, with no autocorrelation at any lag. This is precisely the Ljung-Box test of Section 3.3, applied not to a fitted forecasting model but to the leftover of a decomposition. A significant Ljung-Box statistic on $\hat R_t$ (especially a spike at lag $m$) is the tell-tale sign that the seasonal estimate is too rigid and leftover seasonality is hiding in the remainder, which is exactly the failure classical decomposition suffers on the evolving-amplitude series of subsection three.

The third job is forecasting, and it gives the simplest honest baseline in all of time-series analysis: decompose, forecast each part separately, and recombine. Seasonal-naive handles the seasonal component (next year's December effect equals this year's), a simple drift or linear extrapolation handles the trend, and the remainder, being white noise, is forecast as its mean (zero in the additive case). The recombined forecast $\hat y_{T+h} = \hat T_{T+h} + \hat S_{T+h}$ is transparent and fast, and the statsmodels STLForecast wrapper formalizes it by deseasonalizing, fitting any model you like to the seasonally adjusted series, forecasting it, and adding the seasonal pattern back. Code 3.5.4 builds this baseline explicitly so the recombination is visible, then states the library shortcut.

import numpy as np
from statsmodels.tsa.seasonal import STL
from statsmodels.stats.diagnostic import acorr_ljungbox

rng = np.random.default_rng(3)
t = np.arange(120)
y = 100 + 0.4 * t + 8 * np.sin(2 * np.pi * t / 12) + rng.normal(0, 1.5, size=120)

# 1) Decompose with STL.
stl = STL(y, period=12, robust=True).fit()
T, S, R = stl.trend, stl.seasonal, stl.resid

# 2) Residual diagnostic: the remainder should pass Ljung-Box (Section 3.3).
lb = acorr_ljungbox(R, lags=[12], return_df=True)
print(f"Ljung-Box(12) p-value on remainder: {lb['lb_pvalue'].iloc[0]:.3f}")

# 3) Forecast h steps: trend by linear drift, seasonal by seasonal-naive,
#    remainder by its mean (~0). Then recombine.
h, m = 12, 12
slope = (T[-1] - T[-13]) / 12.0                      # local trend drift per step
trend_fc = T[-1] + slope * np.arange(1, h + 1)       # extrapolate the trend
season_fc = S[-m:][ (np.arange(h)) % m ]             # repeat last full seasonal cycle
y_fc = trend_fc + season_fc                          # recombine; remainder mean is ~0
print("12-step forecast:", np.round(y_fc, 1))

Ljung-Box(12) p-value on remainder: 0.612
12-step forecast: [148.5 152.1 154.7 155.6 154.2 150.8 146.5 142.5 140.0 139.8 141.9 145.6]

6. Looking Back: How Chapter 3 Sets Up the Models Ahead Beginner

This section closes Chapter 3, so it is worth seeing the chapter as one argument. Section 3.1 made stationarity the central concept; Section 3.2 gave the transformations (differencing, logs, Box-Cox) that coax a series toward it; Section 3.3 built the autocorrelation tools and the Ljung-Box test that we just used to grade a remainder; Section 3.4 taught you to read a series as level, trend, season, and noise; and this section turned that reading into estimators that return the components and into the simplest honest forecasting baseline. Together these five sections are the statistical vocabulary that the rest of Part II speaks. Every model in the chapters ahead is a more powerful way to do something you have now done by hand.