Trend and Seasonality

"They keep calling me noise. I am not noise. I arrive every December, on schedule, with a peak you could set your calendar by. Subtract me and you will finally see what is actually going on."

A Seasonal Component That Peaks Punctually Every December

The previous section, Section 3.3, studied autocorrelation and the way a series remembers its own past. That memory is what makes the layers of a series visible: a slow trend shows up as autocorrelation that decays almost not at all, and a seasonal cycle shows up as autocorrelation that spikes at the season length. Here we give those layers names and a generative model. The organizing idea is simple and runs through the whole chapter: before you model the unpredictable part of a series you must first account for the parts that are entirely predictable, the steady climb and the calendar-driven swing, because a model asked to learn structure it could have been told will spend its capacity badly and forecast worse.

1. The Components of a Series Beginner

Classical time-series analysis treats an observed series $X_t$ as the visible result of several unobserved components that evolve on different timescales. There are five worth naming. The level is the baseline value the series fluctuates around at a given moment. The trend $T_t$ is the long-run direction, the slow climb or decline over many periods. The seasonality $S_t$ is the pattern that repeats over a fixed, known period: a day, a week, a year. The cycle is a longer wave-like swing whose period is not fixed and not known in advance, the business cycle being the classic example, which is precisely what distinguishes it from seasonality. The remainder $R_t$ (also called the irregular or residual component) is what is left once the structured parts are removed, the part this chapter's stationary-process machinery from Section 3.2 will model. In most decompositions the cycle is folded into the trend, leaving the trio of trend, seasonality, and remainder that organizes the rest of this section.

The components must be combined into the observation somehow, and there are two dominant rules. The additive model adds them,

and the multiplicative model multiplies them,

The distinction is not cosmetic. In the additive model the seasonal swing has a constant absolute size: if December always adds two hundred units, it adds two hundred whether the level is one thousand or ten thousand. In the multiplicative model the seasonal swing is a constant fraction of the level: if December always adds twenty percent, the absolute size of the December bump grows as the series grows. The deciding question is empirical and visual, illustrated in Figure 3.4.1: does the height of the seasonal wiggle stay flat as the series rises, or does it fan out? Flat means additive; fanning means multiplicative.

The link between the two models is the logarithm, and this is the single most useful fact in the section. Taking logs of the multiplicative model turns the product into a sum,

so a multiplicative series becomes additive in log space. This is exactly why analysts log-transform series whose variance grows with their level (sales, populations, prices, energy demand): the transform stabilizes the variance and lets every additive tool apply. We will use it directly in the worked example of subsection five.

2. Trend: Deterministic or Stochastic Intermediate

The word "trend" hides a deep and dangerous ambiguity, and resolving it is the most important idea in this subsection. A deterministic trend is a fixed function of the clock: the mean of the series is literally $\mu_t = f(t)$ for some known form, most simply a line $\mu_t = \beta_0 + \beta_1 t$, and the random part fluctuates around it,

Here the future is anchored: at any horizon the expected level is $\beta_0 + \beta_1 t$ and the uncertainty band has constant width. A stochastic trend is the opposite. It is a trend with no fixed equation, produced by accumulating random shocks, the random walk of Section 3.2,

which carries a unit root. Both can look like a rising line on a chart, yet they behave in opposite ways and demand opposite treatment. The deterministic trend returns to its line after a shock; the stochastic trend never returns, because each shock $\varepsilon_s$ is permanently baked into the running sum and the variance $\operatorname{Var}(X_t) = \sigma^2 t$ grows without bound. The right cure differs accordingly: you remove a deterministic trend by subtracting the fitted function (detrending by regression), and you remove a stochastic trend by differencing, $\nabla X_t = X_t - X_{t-1}$, which turns the random walk back into stationary noise. Apply the wrong cure and you make things worse, as the warning below records.

Estimating a trend, once you know which kind you have, is straightforward. A deterministic trend is fit by regression on time (a line, a polynomial, or a spline). Either kind can be smoothed out of view by a moving average, which estimates the local level by averaging a window of neighbors. A centered moving average of odd window $2k+1$ is

and its window length is the key knob: too short and seasonality bleeds through, too long and genuine bends in the trend are flattened. The standard trick for an even seasonal period $m$ (so the window cannot be centered on an integer) is the $m$-term centered moving average that half-weights the two endpoints, which we will use in the next section's classical decomposition. Code 3.4.1 fits a linear deterministic trend by least squares and estimates the same trend with a moving average, then contrasts how each behaves on a noisy upward series.

import numpy as np

rng = np.random.default_rng(7)
T = 120
t = np.arange(T)
# A genuine deterministic line plus noise: the trend is a fixed function of t.
true_trend = 5.0 + 0.30 * t
x = true_trend + rng.normal(0, 4.0, size=T)

# (a) Deterministic trend by least-squares regression on time.
A = np.column_stack([np.ones(T), t])            # design matrix [1, t]
beta, *_ = np.linalg.lstsq(A, x, rcond=None)    # beta = [intercept, slope]
trend_reg = A @ beta
print(f"fitted intercept={beta[0]:.2f}, slope={beta[1]:.3f}")

# (b) Trend by a centered moving average of window 2k+1.
def centered_ma(series, k):
    out = np.full_like(series, np.nan, dtype=float)
    for i in range(k, len(series) - k):
        out[i] = series[i - k : i + k + 1].mean()
    return out

trend_ma = centered_ma(x, k=6)                  # window 13
# Compare each estimate to the known truth on the interior points.
interior = slice(6, T - 6)
print(f"regression RMSE vs truth: {np.sqrt(np.mean((trend_reg[interior]-true_trend[interior])**2)):.3f}")
print(f"moving-avg  RMSE vs truth: {np.sqrt(np.nanmean((trend_ma[interior]-true_trend[interior])**2)):.3f}")

fitted intercept=5.36, slope=0.297
regression RMSE vs truth: 0.74
moving-avg  RMSE vs truth: 1.36

3. Seasonality: Dummies and Fourier Harmonics Intermediate

Seasonality is a pattern that repeats over a fixed, known period $m$: $m=7$ for day-of-week structure in daily data, $m=24$ for hour-of-day structure in hourly data, $m=12$ for month-of-year structure in monthly data. Two representations dominate, and the choice between them governs how many parameters you spend. The first is seasonal dummies: one indicator variable per season, so the model carries a separate, freely estimated offset for each of the $m$ phases,

where $D_{j,t}$ is one when $t$ falls in season $j$ and zero otherwise (one dummy is dropped to avoid collinearity with the intercept). Dummies are maximally flexible, capturing any seasonal shape, but they cost $m-1$ parameters, which is fine for $m=12$ and ruinous for $m=365$ or $m=168$ (the weekly cycle in hourly data).

Whichever encoding you pick, the thing being encoded is the punctual, calendar-locked swing that Figure 3.4.2 personifies: a component that peaks at the same phase every cycle, ready to be subtracted before the remainder is modeled.

The second representation is Fourier terms, and it is the workhorse for long and multiple seasonalities. Any periodic shape can be approximated by a sum of sines and cosines at harmonic frequencies of the base period. Using $K$ harmonics for a season of length $m$,

which costs only $2K$ parameters, and $K$ can be far smaller than $m$ when the seasonal shape is smooth. A yearly cycle in daily data ($m=365$) is often captured well by $K=4$ to $10$ harmonics, eight to twenty parameters instead of three hundred sixty-four. The smoothness is the catch: Fourier terms impose a smooth periodic shape and cannot reproduce a spike that lasts a single season, where a dummy can. The two are complementary, not rivals.

The decisive advantage of Fourier terms appears with multiple seasonalities, which are everywhere in modern data. An hourly electricity or sensor series has a daily cycle ($m_1 = 24$), a weekly cycle ($m_2 = 168$), and a yearly cycle ($m_3 = 8766$) layered on top of one another. Dummies cannot handle this gracefully; Fourier terms simply add one harmonic block per period,

each with its own period $m_p$ and its own number of harmonics $K_p$. This additive-in-frequency construction is exactly what Prophet and the TBATS family use under the hood, and it is the engine of the worked example. The frequency-domain reading of these harmonics, why each $(\sin, \cos)$ pair is a single tone in a spectrum, is the subject of Chapter 4.

Beyond these regular cycles sit calendar and holiday effects: a public holiday, a sales event, the leap-year shift, the moving date of Easter or Ramadan. These do not fit a fixed period and are handled as explicit regressors (holiday indicator columns), which Fourier terms and dummies alike accommodate by simply appending the column.

4. Removing Components Before Modeling the Remainder Intermediate

Why pull the series apart at all? Because the remainder is the only part worth handing to a stochastic model, and it is hidden until the structured parts are gone. The strategy has a name, decompose-then-model, and it underlies the entire classical pipeline: estimate and remove the trend, estimate and remove the seasonality, and feed the stationary remainder to an ARMA-style model that captures the short-range dependence the structured terms cannot. There are two distinct ways to remove each component, and they are not interchangeable.

Detrending removes the trend, either by subtracting a fitted deterministic function (regression detrending) or by first differencing $\nabla X_t = X_t - X_{t-1}$ to kill a stochastic trend, the unit-root cure of subsection two. Seasonal differencing removes seasonality by subtracting the value one full season ago,

which annihilates a constant seasonal pattern just as ordinary differencing annihilates a constant drift. This is distinct from seasonal adjustment, which estimates the seasonal component $\hat S_t$ and subtracts it, $X_t - \hat S_t$, leaving the trend intact rather than differencing it away. The distinction matters: seasonal differencing removes both the trend and the season in one stroke and changes the data's interpretation (you are now modeling year-over-year changes); seasonal adjustment isolates and removes only the season, preserving the level and trend you may still want to forecast. The seasonal-difference operator $\nabla_m$ is the seasonal half of the SARIMA models in Chapter 5; the estimate-and-subtract route is the path the decomposition methods of Section 3.5 take.

5. Worked Example: Detrending an Hourly Sensor Series With Multiple Seasonalities Advanced

We now run the full pipeline on the sensor/IoT running dataset introduced in Chapter 2: an hourly telemetry stream that carries a slow upward trend (a sensor drifting as it ages), a strong daily cycle ($m=24$), and a weekly cycle ($m=168$), all buried in noise. The plan is exactly the decompose-then-model recipe: build a design matrix of a linear trend plus Fourier harmonics for both seasonalities, fit it by least squares, subtract the fitted structure, and inspect the remainder to confirm the seasonality is gone. Code 3.4.2 implements the Fourier seasonal regression entirely from scratch so every column of the design matrix is visible.

import numpy as np

rng = np.random.default_rng(11)
N = 24 * 7 * 8                       # eight weeks of hourly data
t = np.arange(N)

# --- Synthesize a series with trend + daily(24) + weekly(168) + noise ---
trend  = 20.0 + 0.004 * t
daily  = 6.0 * np.sin(2*np.pi*t/24)  + 2.0 * np.cos(2*np.pi*2*t/24)   # 2 harmonics
weekly = 4.0 * np.cos(2*np.pi*t/168) + 1.5 * np.sin(2*np.pi*2*t/168)  # 2 harmonics
x = trend + daily + weekly + rng.normal(0, 1.0, size=N)

# --- Build Fourier design columns for one seasonal period ---
def fourier_terms(t, period, K):
    cols = []
    for k in range(1, K + 1):
        cols.append(np.sin(2*np.pi*k*t/period))
        cols.append(np.cos(2*np.pi*k*t/period))
    return np.column_stack(cols)

# Design matrix: [intercept, t, daily harmonics, weekly harmonics]
A = np.column_stack([
    np.ones(N), t,
    fourier_terms(t, 24, K=3),       # daily, 3 harmonics -> 6 columns
    fourier_terms(t, 168, K=3),      # weekly, 3 harmonics -> 6 columns
])
beta, *_ = np.linalg.lstsq(A, x, rcond=None)   # one least-squares solve
fitted = A @ beta
remainder = x - fitted                          # decompose-then-inspect

print(f"design matrix shape: {A.shape}")        # (N, 14)
print(f"recovered trend slope: {beta[1]:.4f}  (true 0.0040)")
print(f"remainder std: {remainder.std():.3f}  (noise std was 1.0)")
# A clean decomposition leaves almost no daily structure in the remainder:
lag24_corr = np.corrcoef(remainder[24:], remainder[:-24])[0, 1]
print(f"remainder autocorr at lag 24: {lag24_corr:.3f}  (near 0 = season removed)")

design matrix shape: (1344, 14)
recovered trend slope: 0.0040  (true 0.0040)
remainder std: 0.987  (noise std was 1.0)
remainder autocorr at lag 24: -0.012  (near 0 = season removed)

The from-scratch version makes the mechanics transparent, but in practice you reach for a library that builds the Fourier features, fits the regression, and forecasts in a handful of lines. Code 3.4.3 reproduces the decomposition with statsmodels for the ordinary-least-squares fit and a compact Fourier helper, the production path.

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.deterministic import DeterministicProcess, Fourier

N = 24 * 7 * 8
idx = pd.RangeIndex(N)
# (reuse x from Code 3.4.2 as the observed series)
series = pd.Series(x, index=idx)

# DeterministicProcess assembles trend + multiple Fourier seasonalities for you.
dp = DeterministicProcess(
    index=idx, constant=True, order=1,                    # intercept + linear trend
    additional_terms=[Fourier(24, order=3),               # daily, 3 harmonics
                      Fourier(168, order=3)],             # weekly, 3 harmonics
)
design = dp.in_sample()
fitted = sm.OLS(series, design).fit()
remainder = series - fitted.predict(design)
print(f"remainder std: {remainder.std():.3f}")