Chapter 46: Time Series Fundamentals

The Four Classical Components

Every time series can be conceptually separated into four additive (or multiplicative) components:

• Trend (<!--MATHBLOCK2-->): the long-run direction — upward, downward, or flat. Monthly global CO₂ concentrations have an unmistakable upward trend.
• Seasonality (<!--MATHBLOCK3-->): periodic fluctuations with a known, fixed period — daily, weekly, annual. Retail sales spike every December.
• Cycle (<!--MATHBLOCK4-->): quasi-periodic swings with variable duration, typically longer than a year. Business cycles and sunspot cycles are canonical examples. In practice, cycles are often lumped with the trend.
• Irregular / Noise (<!--MATHBLOCK5-->): the residual after the other components are removed.

Additive vs. Multiplicative Models

In the additive model the components simply sum:

This is appropriate when the amplitude of the seasonal swings stays roughly constant regardless of the level of the series.

In the multiplicative model the components multiply:

This is appropriate when the seasonal swings scale with the level — a common pattern in economic and demand data. Multiplicative decomposition can be converted to additive by taking logarithms: $\log(Y_t) = \log(T_t) + \log(S_t) + \ldots$

Classical Decomposition with statsmodels

Classical decomposition estimates the trend by a centered moving average and recovers the seasonal component by averaging over corresponding periods:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# Monthly airline passengers (classic Box-Jenkins dataset)
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.columns = ["passengers"]

result = seasonal_decompose(df["passengers"], model="multiplicative", period=12)

fig, axes = plt.subplots(4, 1, figsize=(10, 8), sharex=True)
for ax, component, label in zip(
    axes,
    [df["passengers"], result.trend, result.seasonal, result.resid],
    ["Observed", "Trend", "Seasonal", "Residual"],
):
    ax.plot(component)
    ax.set_ylabel(label)
plt.tight_layout()
plt.show()

The multiplicative model is chosen here because the height of the seasonal spikes grows proportionally with the passenger counts — a tell-tale sign visible in the raw plot.

STL Decomposition

Classical decomposition has two weaknesses: it cannot handle every period length cleanly, and it assumes the seasonal pattern is constant over time. STL (Seasonal-Trend decomposition using Loess) removes both restrictions:

from statsmodels.tsa.seasonal import STL

stl = STL(df["passengers"], period=12, robust=True)
stl_result = stl.fit()
stl_result.plot()
plt.tight_layout()
plt.show()

print(stl_result.seasonal.head(12))  # first year of seasonal factors

The robust=True flag uses iteratively reweighted least squares to down-weight outliers. STL is the decomposition method of choice for modern practitioners.

Interpreting Residuals

After decomposition, inspect the residual component carefully:

• Residuals should look like white noise — no obvious pattern, no heteroskedasticity.
• Remaining autocorrelation in residuals signals that the decomposition has not captured all structure (see the ACF/PACF section).
• Large isolated spikes often correspond to known events (strikes, pandemics, data entry errors).

A common pitfall is treating decomposition output as a forecast. It is not. Decomposition is a diagnostic tool that helps you understand what you are dealing with before choosing a model.

What Is Stationarity?

A time series $\{Y_t\}$ is strictly stationary if the joint distribution of $(Y_{t_1}, Y_{t_2}, \ldots, Y_{t_k})$ is the same as $(Y_{t_1+h}, \ldots, Y_{t_k+h})$ for all $h$ and all choices of time points. This is too strong a requirement in practice.

Weak (covariance) stationarity requires only:

1. Constant mean: <!--MATHBLOCK7--> for all <!--MATHBLOCK8-->.
2. Constant, finite variance: <!--MATHBLOCK9--> for all <!--MATHBLOCK10-->.
3. Autocovariance depends only on lag: <!--MATHBLOCK11--> for all <!--MATHBLOCK12-->.

Most classical models (AR, MA, ARMA) assume weak stationarity. A trending series violates condition 1; a series whose variance explodes over time violates condition 2.

Why Stationarity Matters

If a series is non-stationary, sample statistics (mean, variance, autocorrelations) computed on one part of the series may not generalise to another part — the very foundation of statistical inference collapses. Regressing one non-stationary series on another produces spurious regression: you can get a statistically significant relationship between two completely unrelated trending series (e.g., the number of films Nicolas Cage appeared in vs. swimming pool drownings).

The Unit Root and Random Walk

The most common form of non-stationarity in economic and financial data is the unit root. A simple AR(1) process:

is stationary when $|\phi| < 1$ and has a unit root when $\phi = 1$, yielding a random walk:

The variance of a random walk grows without bound — $\text{Var}(Y_t) = t \sigma^2$ — so it is non-stationary.

Augmented Dickey-Fuller Test

The Augmented Dickey-Fuller (ADF) test formalises the unit root hypothesis. It fits:

The null hypothesis is $H_0: \gamma = 0$ (unit root present, non-stationary). A small p-value (typically $< 0.05$) provides evidence to reject $H_0$ and conclude stationarity.

The lag order $p$ controls for serial correlation in the residuals. statsmodels can select it automatically using AIC or BIC.

import pandas as pd
from statsmodels.tsa.stattools import adfuller

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
y = df["passengers"]

def adf_report(series, title=""):
    result = adfuller(series, autolag="AIC")
    print(f"--- ADF Test: {title} ---")
    print(f"  Test statistic : {result[0]:.4f}")
    print(f"  p-value        : {result[1]:.4f}")
    print(f"  Lags used      : {result[2]}")
    for key, val in result[4].items():
        print(f"  Critical value ({key}): {val:.4f}")
    conclusion = "Non-stationary (fail to reject H0)" if result[1] > 0.05 else "Stationary (reject H0)"
    print(f"  Conclusion     : {conclusion}\n")

adf_report(y, "Raw passengers")
adf_report(y.diff().dropna(), "First-differenced passengers")
adf_report(y.diff(12).dropna(), "Seasonally differenced passengers (lag=12)")

KPSS Test: A Complementary Perspective

The ADF test has low power against near-unit-root processes. The KPSS (Kwiatkowski-Phillips-Schmidt-Shin) test reverses the hypothesis: $H_0$ is stationarity. Using both tests together reduces the chance of error:

• ADF rejects + KPSS does not reject → evidence of stationarity.
• ADF does not reject + KPSS rejects → evidence of unit root.
• Both reject → possibly trend-stationary (stationary around a deterministic trend).
• Neither rejects → insufficient data.

from statsmodels.tsa.stattools import kpss

stat, p_value, lags, critical = kpss(y, regression="c", nlags="auto")
print(f"KPSS stat={stat:.4f}, p={p_value:.4f}")
print("Stationary" if p_value > 0.05 else "Non-stationary")

Visual Checks Before Formal Tests

Never rely solely on formal tests. Plot the rolling mean and standard deviation over a window of, say, 12 periods. A drifting rolling mean or rising rolling standard deviation are immediate visual evidence of non-stationarity. Formal tests should confirm — not replace — visual inspection.

Understanding Autocorrelation

Ordinary correlation measures the linear relationship between two different variables. Autocorrelation measures the linear relationship of a time series with itself at different lags. The autocorrelation at lag $h$ is:

By definition $\rho(0) = 1$. The full set of autocorrelations $\{\rho(h)\}_{h \geq 0}$ is called the autocorrelation function (ACF).

Partial Autocorrelation (PACF)

The ACF at lag $h$ includes indirect correlations mediated through intermediate lags. For example, $Y_t$ is correlated with $Y_{t-2}$ partly because both are correlated with $Y_{t-1}$. The partial autocorrelation at lag $h$, denoted $\phi_{hh}$, removes those intermediate effects — it is the correlation between $Y_t$ and $Y_{t-h}$ after partialling out $Y_{t-1}, \ldots, Y_{t-h+1}$.

Practically, $\phi_{hh}$ is the coefficient on $Y_{t-h}$ in a regression of $Y_t$ on its $h$ most recent lags.

Reading ACF/PACF Plots

ACF and PACF plots are the primary diagnostic tools for selecting ARIMA model orders. The 95 % confidence bands (shown as blue shading or dashed lines) are approximately $\pm 1.96 / \sqrt{n}$ under the white-noise null. Spikes outside these bands are statistically significant.

The classical signatures are:

• AR(p) process: ACF decays geometrically (may oscillate), PACF cuts off sharply after lag <!--MATHBLOCK19-->.
• MA(q) process: ACF cuts off sharply after lag <!--MATHBLOCK20-->, PACF decays geometrically.
• ARMA(p, q): both ACF and PACF decay geometrically — no sharp cutoff.
• Random walk: ACF decays very slowly, often appearing nearly flat — a hallmark of non-stationarity.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

np.random.seed(42)
n = 300

# Pure AR(2): Y_t = 0.6*Y_{t-1} - 0.3*Y_{t-2} + noise
ar2 = np.zeros(n)
for t in range(2, n):
    ar2[t] = 0.6 * ar2[t-1] - 0.3 * ar2[t-2] + np.random.randn()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ar2, lags=30, ax=axes[0], title="ACF of AR(2)")
plot_pacf(ar2, lags=30, ax=axes[1], title="PACF of AR(2)", method="ywm")
plt.tight_layout()
plt.show()
# Expected: PACF spikes at lags 1 and 2, zero thereafter
# ACF decays geometrically (possibly oscillating)

Always difference a non-stationary series before interpreting ACF/PACF. Running ACF on a random walk produces a near-flat, slowly decaying ACF that provides no useful model-order information.

Seasonal ACF Patterns

With seasonal data, you will see significant spikes at multiples of the seasonal period. Monthly data with annual seasonality will show ACF spikes at lags 12, 24, 36, …. After seasonal differencing, these spikes should disappear or greatly reduce.

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
y = df["passengers"]

# Log-transform, then first + seasonal difference
y_transformed = np.log(y).diff().diff(12).dropna()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(y_transformed, lags=36, ax=axes[0],
         title="ACF after log + d=1 + D=1 (m=12)")
plot_pacf(y_transformed, lags=36, ax=axes[1],
          title="PACF after log + d=1 + D=1 (m=12)", method="ywm")
plt.tight_layout()
plt.show()

After differencing, the residual ACF/PACF pattern guides selection of the non-seasonal $(p, q)$ and seasonal $(P, Q)$ orders for a SARIMA model.

Ljung-Box Test for Residual Autocorrelation

After fitting a model, the Ljung-Box test checks whether any autocorrelation remains in the residuals. The test statistic is:

Under $H_0$ (no autocorrelation up to lag $m$), $Q \sim \chi^2_m$. A significant result means the model has not captured all temporal structure.

from statsmodels.stats.diagnostic import acorr_ljungbox

residuals = y_transformed  # using differenced series as a placeholder
lb = acorr_ljungbox(residuals, lags=range(1, 13), return_df=True)
print(lb[["lb_stat", "lb_pvalue"]].round(4))

In production model validation, always run the Ljung-Box test on model residuals at several lag lengths.

Why Differencing?

Differencing is the primary technique for removing non-stationarity induced by trends and unit roots. Rather than modelling the level of a series, we model its changes. The first difference is:

For a random walk $Y_t = Y_{t-1} + \varepsilon_t$, the first difference $\Delta Y_t = \varepsilon_t$ is white noise — perfectly stationary. In general, if a series requires $d$ rounds of differencing to become stationary, it is said to be integrated of order $d$, written $I(d)$.

The backshift operator $B$ provides clean notation: $B Y_t = Y_{t-1}$, so $\Delta Y_t = (1 - B) Y_t$ and $d$-order differencing is $(1 - B)^d Y_t$.

Seasonal Differencing

For series with period $m$, seasonal differencing removes the seasonal pattern:

For monthly data with annual seasonality, $m = 12$. A common strategy for seasonal time series is to apply both non-seasonal and seasonal differencing:

This is the transformation used in SARIMA$(p,1,q)(P,1,Q)_{12}$ models.

How Much Differencing is Enough?

Over-differencing is a real hazard. If a series is already stationary and you difference it, you introduce autocorrelation — specifically, you create an MA(1) component. The signature of over-differencing is a first-lag ACF near $-0.5$.

Guidelines:

• If the first-lag ACF is positive (common with trending series), the series likely needs differencing.
• If the first-lag ACF is more negative than <!--MATHBLOCK18-->, you may have over-differenced.
• Formal ADF/KPSS tests (Section 2) guide the decision systematically.
• Usually <!--MATHBLOCK19--> and <!--MATHBLOCK20--> suffice.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
y = df["passengers"]

# --- Step 1: log transform to stabilise variance (multiplicative -> additive)
log_y = np.log(y)

# --- Step 2: seasonal difference to remove seasonality
d12 = log_y.diff(12)  # Y_t - Y_{t-12}

# --- Step 3: first difference to remove trend
d1d12 = d12.diff(1)  # (1-B)(1-B^12) log Y_t

fig, axes = plt.subplots(4, 1, figsize=(10, 10), sharex=True)
for ax, series, title in zip(
    axes,
    [y, log_y, d12, d1d12],
    ["Original", "Log transform", "Log + seasonal diff (D=1)", "Log + seasonal + first diff (d=1, D=1)"]
):
    ax.plot(series, linewidth=0.8)
    ax.set_title(title)
plt.tight_layout()
plt.show()

Common Variance-Stabilising Transformations

Differencing addresses the mean; transformations address non-constant variance.

• Log transform <!--MATHBLOCK21-->: suitable when variance grows proportionally with level (multiplicative model). Requires <!--MATHBLOCK22-->.
• Square root <!--MATHBLOCK23-->: appropriate for count data following a Poisson-like distribution.
• Box-Cox transformation:

The Box-Cox parameter $\lambda$ can be estimated from data by maximum likelihood.

from scipy.stats import boxcox
import numpy as np

# Estimate optimal lambda from airline data
bc_transformed, lambda_opt = boxcox(y.values)
print(f"Optimal Box-Cox lambda: {lambda_opt:.4f}")
# A lambda near 0 confirms the log transform is appropriate

Pitfalls

Pitfall 1: Differencing removes the absolute level of the series. To recover forecasts on the original scale, you must integrate (cumulatively sum) the differenced forecasts, then undo any log transform — in the correct order.

Pitfall 2: Each order of differencing removes one observation from the beginning of the series. First-differencing loses one; seasonal differencing loses $m$; both lose $m + 1$ observations. For short series this matters.

Pitfall 3: Do not difference purely to match the data to a model. If you remove the trend via differencing but the trend is actually deterministic (linear in $t$), you have over-specified your model. A simple linear regression detrend may be preferable.

Pandas as a Time Series Engine

Pandas was designed with time series in mind. A DatetimeIndex unlocks a rich API for slicing, resampling, and rolling operations. Mastering this API is a prerequisite for any real-world forecasting project.

Parsing and Indexing

Always parse dates at load time and set the datetime column as the index:

import pandas as pd

df = pd.read_csv("data.csv", parse_dates=["timestamp"], index_col="timestamp")
df = df.sort_index()          # ensure chronological order
df = df.asfreq("D")           # fill gaps to make the frequency explicit
print(df.index.freq)          # <Day>

Common frequency strings: "T" (minute), "H" (hour), "D" (day), "B" (business day), "W" (week-end), "MS" (month start), "ME" (month end), "QS" (quarter start), "YS" (year start).

Date Slicing

Pandas supports partial string indexing — you can slice by year, month, or specific ranges:

df["2020"]               # all of 2020
df["2020-01":"2020-06"]  # Jan-Jun 2020
df.loc["2020-03-15"]     # single date

Resampling

Resampling changes the frequency of a series. Downsampling aggregates (e.g., hourly to daily); upsampling interpolates.

import pandas as pd
import numpy as np

np.random.seed(1)
idx = pd.date_range("2020-01-01", periods=365, freq="D")
daily = pd.Series(np.random.randn(365) + 10, index=idx, name="value")

# Downsample to monthly sums, means, and last values
monthly_mean = daily.resample("ME").mean()
monthly_sum  = daily.resample("ME").sum()
monthly_last = daily.resample("ME").last()

# Upsample from monthly to daily, forward-fill
monthly_ff = monthly_mean.resample("D").ffill()

print(monthly_mean.head(3).round(3))

For OHLC financial data use .resample("W").ohlc() to get Open/High/Low/Close bars.

Rolling Windows

Rolling statistics compute a statistic over a sliding window of fixed size. They are the workhorse of signal processing and feature engineering for time-series ML models.

roll = daily.rolling(window=30, min_periods=1)

roll_mean = roll.mean()    # 30-day moving average
roll_std  = roll.std()     # 30-day rolling volatility
roll_max  = roll.max()

import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
axes[0].plot(daily, alpha=0.4, label="Daily")
axes[0].plot(roll_mean, linewidth=2, label="30-day MA")
axes[0].fill_between(daily.index, roll_mean - roll_std, roll_mean + roll_std,
                     alpha=0.3, label="±1 std band")
axes[0].legend()
axes[1].plot(roll_std, color="orange", label="30-day rolling std")
axes[1].legend()
plt.tight_layout()
plt.show()

Expanding Windows

An expanding window grows from the beginning of the series to the current point — it computes a statistic over all observations up to and including $t$. This is appropriate when you want a running total or cumulative mean:

expanding_mean = daily.expanding(min_periods=1).mean()
expanding_std  = daily.expanding(min_periods=1).std()

Expanding windows are commonly used in walk-forward validation (Section 6) to simulate how a model would have performed if retrained at each step.

Time-Based Rolling Windows

Pandas also supports time-offset windows (e.g., the last 30 days) rather than a fixed number of observations. This handles irregular timestamps correctly:

irregular_idx = idx[np.sort(np.random.choice(365, 200, replace=False))]
irregular = pd.Series(np.random.randn(200) + 10, index=irregular_idx)

# Rolling mean over the trailing 30 calendar days
time_roll = irregular.rolling("30D").mean()
print(time_roll.tail(5).round(3))

Lag Features for ML Models

Tree-based and neural network models for time series typically require explicit lag features. Pandas .shift() creates them:

df_features = pd.DataFrame({"y": daily})
for lag in [1, 7, 14, 30]:
    df_features[f"lag_{lag}"] = df_features["y"].shift(lag)

# Rolling summary features
df_features["rolling_mean_7"] = df_features["y"].shift(1).rolling(7).mean()
df_features["rolling_std_7"]  = df_features["y"].shift(1).rolling(7).std()
df_features = df_features.dropna()
print(df_features.head(3))

Note the .shift(1) before the rolling aggregation. This ensures that when computing the rolling mean for row $t$, you are using observations from $t-1, t-2, \ldots$ only — not $t$ itself. Forgetting this shift is one of the most common sources of data leakage in time series ML.

The Cardinal Rule: Never Shuffle Time Series

In standard supervised learning, randomly shuffling and splitting data is good practice — it ensures the train and test sets have similar distributions. In time series, random shuffling is data leakage. If a test observation from 2023 is predicted using training data that includes 2024, your model has access to the future. Performance estimates will be wildly optimistic.

The only valid split strategy is chronological: train on the past, evaluate on the future.

Simple Holdout Split

Reserve the last $k$ time steps as the test set:

import pandas as pd
import numpy as np

np.random.seed(0)
n = 200
idx = pd.date_range("2019-01-01", periods=n, freq="ME")
y = pd.Series(np.cumsum(np.random.randn(n)) + 100, index=idx)

# Holdout: last 24 months as test
train = y.iloc[:-24]
test  = y.iloc[-24:]

print(f"Train: {train.index[0].date()} to {train.index[-1].date()} ({len(train)} obs)")
print(f"Test : {test.index[0].date()} to {test.index[-1].date()} ({len(test)} obs)")

Choose the holdout size to match your actual forecast horizon. If you will forecast 3 months ahead in production, a 3-month (or longer) test set is appropriate.

Walk-Forward (Rolling Origin) Validation

A single holdout split gives only one performance estimate. Walk-forward validation — also called expanding-window or rolling-origin cross-validation — provides many estimates by sliding the test window forward in time:

Fold 1: [===Train===|Test]
Fold 2: [====Train====|Test]
Fold 3: [=====Train=====|Test]

In the expanding window variant the training set grows; in the sliding window variant it stays fixed in size (simulating a model retrained on only the most recent data).

from sklearn.model_selection import TimeSeriesSplit
import numpy as np
import matplotlib.pyplot as plt

X = np.arange(n).reshape(-1, 1)
y_arr = y.values

tscv = TimeSeriesSplit(n_splits=5, test_size=24)

fig, ax = plt.subplots(figsize=(12, 5))
for fold, (train_idx, test_idx) in enumerate(tscv.split(X)):
    ax.scatter(train_idx, [fold] * len(train_idx), c="steelblue", s=6, label="Train" if fold == 0 else "")
    ax.scatter(test_idx,  [fold] * len(test_idx),  c="tomato",   s=6, label="Test"  if fold == 0 else "")
    print(f"Fold {fold+1}: train {len(train_idx)}, test {len(test_idx)}")
ax.set_xlabel("Time step")
ax.set_ylabel("Fold")
ax.legend()
ax.set_title("TimeSeriesSplit — Expanding Window")
plt.tight_layout()
plt.show()

Gap Between Train and Test

For multi-step forecasting with a horizon $h$, you should leave a gap of $h$ periods between the last training point and the first test point. This ensures that any features derived from recent observations (rolling windows, lag features) could have actually been computed in production at forecast time. TimeSeriesSplit does not do this automatically:

gap = 3  # 3-month forecast horizon
tscv_gap = TimeSeriesSplit(n_splits=5, test_size=12, gap=gap)
for fold, (tr, te) in enumerate(tscv_gap.split(X)):
    print(f"Fold {fold+1}: last train idx={tr[-1]}, first test idx={te[0]}, gap={te[0]-tr[-1]-1}")

Evaluation Metrics for Time Series

Standard regression metrics (MAE, RMSE, MAPE) apply, but beware:

• MAPE (Mean Absolute Percentage Error) is undefined when actual values are zero and heavily penalises underpredictions when values are small. Use sMAPE or MASE for intermittent demand series.
• MASE (Mean Absolute Scaled Error) scales errors by the in-sample naive forecast error, making it scale-independent and comparable across series.
• Report metrics per forecast horizon if you care about accuracy at specific lead times. A model with low average error may be terrible at horizon <!--MATHBLOCK3--> and great at <!--MATHBLOCK4-->.

import numpy as np

def mase(y_train, y_test, y_pred, m=1):
    """Mean Absolute Scaled Error. m=seasonality period (1=non-seasonal naive)."""
    naive_errors = np.abs(np.diff(y_train, n=m))
    scale = naive_errors.mean()
    return np.abs(y_test - y_pred).mean() / scale

# Example: naive forecast (last value)
y_train_arr = y.iloc[:-24].values
y_test_arr  = y.iloc[-24:].values
y_naive     = np.full(24, y_train_arr[-1])

mae_val  = np.abs(y_test_arr - y_naive).mean()
rmse_val = np.sqrt(((y_test_arr - y_naive) ** 2).mean())
mase_val = mase(y_train_arr, y_test_arr, y_naive)

print(f"Naive forecast — MAE: {mae_val:.2f}, RMSE: {rmse_val:.2f}, MASE: {mase_val:.2f}")

A Note on Stationarity and the Test Set

Even if you have correctly split by time, a non-stationary series can still cause problems: the test set may be in a different distributional regime than the training set (e.g., post-COVID demand patterns vs. pre-COVID). No amount of careful splitting fixes distributional shift — that requires domain knowledge, regime detection, or models that adapt online. The split only prevents look-ahead leakage; it does not guarantee that the future looks like the past.