Data with an order

In ordinary tabular learning the rows are exchangeable — shuffle them and nothing is lost. A time series is the opposite: shuffle it and you destroy the very structure you care about. Each observation is tied to a timestamp, consecutive observations are correlated, and the questions are temporal — what's the underlying direction, does it repeat, is today's spike unusual, what comes next? Answering them starts with a simple decomposition.

Trend, season, noise

The classical model says any series is a sum of a few interpretable pieces:

Add the pieces up

Run the model forwards: dial a trend slope, a seasonal amplitude and period, and a dose of noise, and watch the three components stack into a single series. This is exactly how the synthetic data underneath every forecasting demo is made — and seeing it assembled makes the inverse problem (pulling it back apart) far less mysterious:

Turn the seasonal amplitude to zero and the repeating ripples vanish, leaving trend plus noise; kill the noise and the series becomes a clean deterministic curve. Real series are these same ingredients in unknown proportions — and decomposition is the act of estimating them.

Pull a real series apart

Here is the inverse, on a classic dataset: monthly international airline passengers, 1949–1960 — 144 points with an unmistakable upward trend and a yearly summer peak. Classical decomposition estimates the trend with a centered 12-month moving average, finds the seasonal shape by averaging each calendar month after detrending, and calls whatever remains the residual. The four stacked panels are the whole story of the series:

The same computation, laid out like a spreadsheet — each column derived from the ones before it, so you can trace a single month from raw count to residual (rows 6–13, mid-1949 to early-1950):

The trend panel shows steady growth; the seasonal panel shows the same yearly shape repeated twelve times; and the residual — crucially — is not flat. Its swings grow as the series rises, a tell that this series is multiplicative, not additive.

Do the swings grow?

The additive model y = T + S + R assumes the seasonal swing is a fixed number of passengers, the same in 1949 as in 1960. But airline seasonality is proportional — a 20% summer bump on a bigger base is a bigger absolute jump. That's the multiplicative model:

The fix is a trick you've seen in the feature-engineering world: take logs. Since log(T × S × R) = log T + log S + log R, a multiplicative series becomes additive on the log scale, with a flat-variance residual — which is why airline-passenger forecasts are almost always done in logs. The signature to watch for: if the seasonal swings fan out as the level rises, go multiplicative (or log-transform); if they stay a constant width, additive is fine.

See it directly. Toggle between the raw passenger counts and their logarithm, with each year's seasonal high and low traced as an envelope. On the raw scale the envelope fans open as traffic grows; on the log scale the two edges run parallel — the swing has become a constant proportion:

The slow direction

The trend is the low-frequency backbone. The simplest estimate is a moving average — slide a window and average — which smooths away the season and noise to leave the drift; a window equal to the seasonal period cancels the season exactly. Trends needn't be straight: they bend, saturate, and break. Fitting a line assumes constant growth; a moving average or local regression lets the trend curve with the data, which is why decomposition prefers the smoother. The trend is what you extrapolate when you forecast the long run.

The part that repeats

Seasonality is structure on a known, fixed period: 12 months in airline data, 7 days in web traffic, 24 hours in electricity demand — and real series often carry several at once (hour-of-day and day-of-week). You detect a candidate period from domain knowledge, from peaks in the autocorrelation, or from the calendar; you estimate its shape by averaging every point that shares a season-position after removing the trend. Subtracting the seasonal component gives the seasonally adjusted series — how economists report unemployment, so a January dip isn't mistaken for a real decline.

What's left over

After trend and season are removed, the residual is the part the decomposition can't explain — and ideally it's pure noise: no trend, no repeating pattern, constant spread, centred on zero (or one, for multiplicative). If structure remains — drift, leftover seasonality, fanning variance — the decomposition has missed something, or the model is the wrong type. The residual is your diagnostic: a flat, patternless residual means trend and season captured everything systematic, and what's left is genuine randomness. Checking it for leftover autocorrelation is the bridge to the next dive.

What it buys you

From scratch, then statsmodels

Classical decomposition is nothing but a moving average plus averaging-by-month — worth writing once in plain NumPy before reaching for a library:

# --- classical additive decomposition, from scratch ---
def decompose(y, P=12):
    n, h = len(y), P // 2
    trend = np.full(n, np.nan)
    for i in range(h, n - h):                         # centered 2xP moving average
        window = np.r_[0.5*y[i-h], y[i-h+1:i+h], 0.5*y[i+h]]   # half-weight the two ends
        trend[i] = window.sum() / P
    detrended = y - trend
    seasonal_idx = np.array([np.nanmean(detrended[m::P]) for m in range(P)])
    seasonal_idx -= seasonal_idx.mean()               # center so the season sums to zero
    seasonal = seasonal_idx[np.arange(n) % P]
    resid = y - trend - seasonal
    return trend, seasonal, resid

The first computable trend point averages months 1–13 with the two ends half-weighted: (½·112 + 118 + … + 104 + ½·115) / 12 = 126.79 passengers. statsmodels packages this exact recipe — plus the multiplicative and STL variants — and agrees with the hand-rolled version to machine precision:

from statsmodels.tsa.seasonal import seasonal_decompose, STL

res   = seasonal_decompose(y, model="additive", period=12)        # classical, additive
res_m = seasonal_decompose(y, model="multiplicative", period=12)  # or log-then-additive
stl   = STL(y, period=12, robust=True).fit()                      # lets the season drift over time

from-scratch vs statsmodels — identical
max|trend    - res.trend|    = 1.1e-13
max|seasonal - res.seasonal| = 2.8e-14
trend     climbs 127 -> 475 ;  seasonal peaks +64 (Jul), troughs -54 (Nov)
resid     swings GROW with the level -> prefer multiplicative / log

The shape of every series

Every time series is, to first approximation, a trend plus a repeating season plus noise — and the first thing to do with one is pull those apart. Decomposition gives you a baseline forecast, a seasonally adjusted view, an anomaly detector, and a residual to diagnose. Watch the residual's variance to choose additive vs multiplicative, smooth to estimate the trend, average-by-season to estimate the pattern, and treat any leftover structure as a clue you're not done. Everything later — smoothing, ARIMA, forecasting, anomaly detection — builds on this anatomy.

Correlation with your own past

Ordinary correlation measures how two different variables move together. Autocorrelation asks the same question of one series against a time-shifted version of itself: do high values tend to follow high values? The answer, computed at every lag, is a fingerprint. A slowly fading fingerprint means a trend; a fingerprint that spikes every twelve steps means a yearly season; a fingerprint that vanishes after one step means each value depends only on the one before. Reading that fingerprint is how you choose what to do next.

The autocorrelation function

For lag k, the sample autocorrelation is the lag-k covariance divided by the variance — correlation, restricted to pairs k steps apart:

Note r₀ = 1 always (a series is perfectly correlated with itself), and the denominator is fixed — every lag is normalised by the same total variance, so the r_k shrink naturally as k grows and fewer pairs overlap.

By hand, on a five-point series x = [2, 4, 6, 5, 3] with mean x̄ = 4, the deviations are [−2, 0, 2, 1, −1]. The lag-1 numerator pairs each deviation with the next: (−2)(0) + (0)(2) + (2)(1) + (1)(−1) = 1. The denominator is the sum of squared deviations: 4 + 0 + 4 + 1 + 1 = 10. So r₁ = 1 / 10 = 0.10 — almost no lag-1 memory.

Here is that arithmetic as a spreadsheet — the row number is t, and each column to the right is computed from the one before:

It's just a lagged scatter

Strip away the formula and autocorrelation is the most ordinary thing imaginable: shift the series by k, plot each value against its k-step-ago self, and measure the correlation of that scatter. Below is a series with real memory; drag the lag and watch the cloud tighten into a diagonal at small lags and dissolve into a shapeless blob as the memory fades — the slope of the cloud is r_k:

At lag 1 the points line up tightly — this series strongly resembles its previous value. Push the lag out and the correlation decays toward zero: the influence of the past leaks away. That decay curve, lag by lag, is the autocorrelation function.

Every lag at once

Collect r_k for k = 1, 2, 3, … and draw them as a stem plot: the correlogram. It's the workhorse diagnostic. Around zero sit dashed confidence bands at ±1.96/√n — under pure noise, 95% of the bars fall inside, so a stem poking past the band is a statistically real correlation. Three shapes recur and each means something definite: a slow linear decay signals a trend (every nearby point is similar); a spike at the seasonal period (lag 12 for monthly data) signals seasonality; and bars that stay inside the bands everywhere mean there's no linear structure left to exploit.

Two processes, two fingerprints

Two simple models generate most of the textbook intuition. An autoregressive AR(1) process feeds its own past forward, y_t = φ y_t−1 + ε_t; the influence compounds, so its ACF decays geometrically, r_k = φ^k — never quite reaching zero. A moving-average MA(1) process carries only one shock forward, y_t = ε_t + θ ε_t−1; values two or more steps apart share no shock, so its ACF cuts off sharply after lag 1, r₁ = θ/(1+θ²), r_k>1 = 0. Decay versus cut-off is the first thing to read off a correlogram.

Direct memory, echoes removed

The ACF has a blind spot. If y_t correlates with y_t−1, and y_t−1 with y_t−2, then y_t and y_t−2 will correlate too — purely by transitivity, even with no direct link. The partial autocorrelation (PACF) removes that echo: it is the correlation between y_t and y_t−k after regressing out everything in between. It answers "how much new information does lag k add?" The PACF is computed straight from the ACF by the Durbin–Levinson recursion, which solves the lag-k linear-prediction problem one order at a time and reads off the last coefficient.

Read the pair, name the model

Together the ACF and PACF identify a process — this is the heart of Box–Jenkins model selection. Switch between pure noise, an AR(1), and an MA(1), and drag the coefficient; watch the two plots swap roles. The rule is exact and worth memorising:

So an ACF that decays with a PACF that snaps to zero after one spike says "AR(1)"; the mirror image says "MA(1)". The lag where the cut-off happens is the order. That single reading turns a mystery series into a model.

The airline series, decoded

Here is the ACF of the airline-passenger series — and it tells the whole story at a glance. Raw, the bars decay agonisingly slowly (the trend: every month resembles its neighbours) with a clear bump at lag 12 (the yearly season). After first differencing — subtracting each month from the next to kill the trend — the slow decay collapses, but a tall spike survives at lag 12: the seasonality is still there, now unmistakable. Below the ACF sits the PACF — the direct lags with the echoes stripped out. Toggle and watch the trend's smear give way to the season's spike:

This is exactly how an analyst diagnoses what a series needs: the raw ACF screams "non-stationary — difference me," and the differenced ACF's surviving lag-12 spike says "and take a seasonal difference too." The correlogram doesn't just describe the series; it prescribes the next step.

The goal is no memory left

Decomposition and modelling aim at one target: residuals that are white noise — zero autocorrelation at every lag, every bar inside the bands. If structure remains, the model has left signal on the table. Eyeballing the correlogram is the quick check; the formal one is the Ljung–Box portmanteau test, which bundles the first h lags into a single statistic

A large Q (small p-value) rejects "the residuals are white" — there's leftover autocorrelation to model. A small Q means you're done: what remains is genuine noise, and a forecast can do no better than its mean.

See the test decide. Toggle between pure white noise and the AR(1) series from earlier: white noise keeps almost every bar inside the bands and Q below the χ² threshold; the AR series breaks both, again and again:

From scratch, then statsmodels

The ACF is two dot products; the PACF is the Durbin–Levinson recursion on top. Both are short enough to write plainly before trusting a library:

# --- ACF and PACF from scratch ---
def acf(x, K):
    x = x - x.mean()
    c0 = (x * x).sum()                                # total variance (the denominator)
    return np.array([ (x[k:] * x[:len(x)-k]).sum() / c0 for k in range(K+1) ])

def pacf(acf):                                        # Durbin-Levinson recursion from the ACF
    K = len(acf) - 1
    phi = np.zeros((K+1, K+1)); pac = np.zeros(K+1); pac[0] = 1
    phi[1, 1] = acf[1]; pac[1] = acf[1]
    for k in range(2, K+1):
        num = acf[k] - sum(phi[k-1, j] * acf[k-j] for j in range(1, k))
        den = 1      - sum(phi[k-1, j] * acf[j]   for j in range(1, k))
        phi[k, k] = num / den; pac[k] = phi[k, k]
        for j in range(1, k):
            phi[k, j] = phi[k-1, j] - phi[k, k] * phi[k-1, k-j]
    return pac

statsmodels does both in one call each — and agrees with the hand-rolled versions to machine precision on an AR(1) series (φ = 0.75):

from statsmodels.tsa.stattools import acf as sm_acf, pacf as sm_pacf
sm_acf(y, nlags=12, fft=False)              # sample ACF
sm_pacf(y, nlags=12, method="ywm")          # sample PACF

from-scratch vs statsmodels — identical
ACF  lag1 = 0.760   max|diff vs statsmodels| = 2.2e-16
PACF lag1 = 0.760 , lag2 = -0.090   (AR(1): PACF cuts off after lag 1)

One plot, the whole diagnosis

Autocorrelation is the memory of a series — correlation with its own lagged copy. The ACF across lags reveals trend (slow decay), seasonality (periodic spikes), and model order; the PACF strips out the transitive echoes to expose direct dependence. Read as a pair, they name the process: ACF decays + PACF cuts → AR; ACF cuts + PACF decays → MA. The bands tell you what's real, the Ljung–Box test confirms when nothing's left, and the airline ACF shows the workflow in miniature — diagnose, difference, repeat. Every model in the dives ahead is chosen by reading these two plots.

Statistics that don't drift

Take two windows of a series — one early, one late — and compute their mean, spread, and lag-1 correlation. If those summaries match, the series is stationary; the process generating it is the same in both eras, so what you learn from the past transfers to the future. The airline series fails immediately: its later years have a higher mean (the trend) and bigger swings (the growing season). A model fit to it learns numbers that won't hold next year. Stationarity is the property that makes "the past predicts the future" actually true.

Three reasons to insist on it

Weak (covariance) stationarity

The usable definition asks for three things to hold for all time t:

The mean is constant, the variance is constant, and the autocovariance between two points depends only on the gap k between them, never on when they occur. (This is weak stationarity; strict stationarity demands the entire distribution be time-invariant, which is stronger than anyone needs in practice.) Note the third condition is exactly what made the ACF well-defined in the previous dive — the autocorrelation depends on lag alone, which is only true for a stationary series.

Roll a window through

The quickest diagnostic is visual: slide a window along the series and plot the rolling mean and a rolling ±1σ band. On a stationary series both sit flat; on a non-stationary one the mean wanders and the band breathes. Toggle the raw airline series against its twice-differenced log — the difference is stark:

Raw, the rolling mean climbs from ~130 to ~470 and the band fans open — both the mean and the variance are functions of time, the textbook fingerprint of non-stationarity. After differencing, the rolling mean pins to zero and the band holds a near-constant width: the statistics have stopped drifting. The other tell you already know — a correlogram whose bars decay agonisingly slowly (last dive's raw airline ACF) is the autocorrelation signature of the same problem.

Subtract the recent past

The fix for a drifting mean is differencing — replace each value with its change since the step before:

A constant linear trend has a constant slope, so its first difference is a constant — the trend flattens to a level. The number of differences needed to reach stationarity is the order of integration d; most economic series need d = 1, occasionally 2, rarely more. (Over-differencing is a real mistake: difference a stationary series and you inject artificial negative autocorrelation, so stop as soon as the ADF test passes.) By hand, the first six first-differences of the airline series are 118−112, 132−118, 129−132, … = 6, 14, −3, −8, 14, 13 — the trend is gone, only the wiggles remain.

Here is differencing as a spreadsheet — each new column is a subtraction of a shifted copy of the input (rows 13–18):

Subtract last year

A first difference flattens a trend but leaves a fixed seasonal pattern untouched — every December still towers over its November. The cure is a seasonal difference: subtract the value from one full season ago,

with s = 12 for monthly data. Comparing each month to the same month last year cancels a stable yearly shape. Series with both trend and seasonality — like the airline data — typically need both a regular and a seasonal difference, ΔΔ₁₂y, applied in either order. That double-differenced log is exactly the stationary series the rolling-statistics panel above settled flat.

Logs for the variance

Differencing fixes a drifting mean; it does nothing for a growing variance. That's the multiplicative fanning from the anatomy dive, and the fix is the same: take logs (or a Box–Cox transform) first, to stabilise the spread, then difference to flatten the level. The standard recipe for a series like this is therefore log → seasonal difference → first difference. Each tool targets one defect: the log tames the variance, the differences tame the mean and the season.

A formal stationarity check

Eyeballing rolling statistics is good; a hypothesis test is better. The augmented Dickey–Fuller test fits the regression

and tests the null H₀: γ = 0, a unit root — meaning the series does not pull back toward any level, the hallmark of non-stationarity. If γ is reliably negative, shocks decay and the series is mean-reverting (stationary). The test statistic is the t-ratio γ̂/se, but it is not compared to the normal distribution — under a unit root it follows a special Dickey–Fuller distribution, so you compare against MacKinnon critical values (around −2.9 at 5%). A statistic more negative than the critical value, or equivalently p < 0.05, rejects the unit root: the series is stationary. The lagged Δy terms are the "augmented" part — they soak up autocorrelation so the test on γ stays valid.

What does γ actually control? Exactly how hard the series is yanked back toward its level. Write the same idea as an autoregression, y_t = φ·y_t−1 + ε_t with φ = 1 + γ, and drag φ through the three regimes the test distinguishes — same random shocks throughout, only the pull-back strength changes:

Below 1, every shock decays and the series oscillates tightly around zero — stationary, the case the ADF test wants to confirm. At exactly φ = 1 (γ = 0) the pull-back vanishes: shocks accumulate forever and the series wanders off as a random walk — the unit root, the null hypothesis. Above 1 the feedback compounds and the series explodes. The whole ADF test is just a careful way of asking which side of φ = 1 your series sits on.

Climb the ladder to stationarity

Put it together on the airline series and watch the ADF verdict flip. Step through the transform ladder — raw, then log, then a first difference, a seasonal difference, and finally both — and read the statistic and p-value at each rung:

The ladder tells the whole story: raw, p = 0.99 (hopelessly non-stationary); logging alone barely moves it; a single difference or a single seasonal difference each gets close but stalls around p ≈ 0.07; only with both differences applied does the statistic plunge to −4.4 and p to 0.000 — stationary at last. That ΔΔ₁₂ log series is what an ARIMA model would actually be fit on.

From scratch, then statsmodels

Differencing is one NumPy call; the Dickey–Fuller statistic is one small regression. Both are short enough to write before trusting the library:

import numpy as np

# --- transforms from scratch ---
ly   = np.log(y)
d1   = np.diff(ly)                       # first difference  (kills trend)
d12  = ly[12:] - ly[:-12]                # seasonal difference (kills season)
dd   = np.diff(d12)                       # both -> stationary

def df_tstat(s):                          # Dickey-Fuller: regress dY on Y_{t-1}
    Y, dY = s[:-1], np.diff(s)
    X = np.column_stack([np.ones(len(Y)), Y])
    b = np.linalg.lstsq(X, dY, rcond=None)[0]
    resid = dY - X @ b
    s2  = resid @ resid / (len(dY) - 2)
    se  = np.sqrt(s2 * np.linalg.inv(X.T @ X)[1, 1])
    return b[1] / se                      # t-ratio on gamma

statsmodels' adfuller wraps the same regression, adds the augmenting lags automatically, and supplies the MacKinnon p-value — agreeing with the hand-rolled core to machine precision:

from statsmodels.tsa.stattools import adfuller
stat, p, *_ = adfuller(dd)               # augmented, MacKinnon p-value

from-scratch vs statsmodels — final ΔΔ₁₂ log series
from-scratch DF t-stat   = -16.1911
statsmodels (maxlag=0)   = -16.1911      max|diff| = 1e-13
adfuller(autolag="AIC")  ADF = -4.44 , p = 0.000   ->  STATIONARY

Tame it before you model it

Stationarity — constant mean, variance, and lag-only autocovariance — is the assumption that makes the past a guide to the future, and most real series lack it. Spot non-stationarity with a wandering rolling mean, a breathing variance band, or a slowly-decaying ACF; remove a drifting mean with a first difference, a fixed season with a seasonal difference, and a growing variance with a log; and confirm the result with the augmented Dickey–Fuller test, reading p < 0.05 as "stationary." For the airline series the full recipe is log, then both differences — and the once-wild curve becomes the flat, modellable series that ARIMA is built to forecast.

Three letters, three jobs

ARIMA is a recipe with three dials. The I comes first in practice: difference the series d times until it's stationary (the previous dive's job). On that stationary series, AR(p) says today is a linear function of the last p values, and MA(q) says today is a linear function of the last q forecast errors. Add a constant, fit the coefficients to history, and iterate the equation forward to forecast. Almost every classical single-series forecast — and the baseline every fancier method must beat — is some ARIMA.

Regress on your own past

An autoregression predicts each value from a weighted sum of the previous ones:

AR(1) with 0 < φ < 1 decays smoothly back to its mean (the mean-reverting case from the last dive). AR(2) is richer: with the right coefficients its roots turn complex and the series oscillates, tracing pseudo-cycles no AR(1) can. Drag the two coefficients and watch the behaviour change — and watch it blow up when you leave the stationary region:

The key fact for fitting: because each y_t is a linear function of earlier values plus noise, fitting AR is nothing more than ordinary least squares — regress the series on lagged copies of itself. That makes the AR part trivially solvable, which is half of why ARIMA endures.

A memory of shocks

A moving-average term models today as a weighted sum of recent random shocks — not of past values:

This is a different beast from the moving-average smoother of the anatomy dive — here the ε are unobserved white-noise innovations, and an MA(q) process has memory exactly q steps long: a shock today influences the next q values, then is gone forever. That's why an MA(q)'s ACF cuts off sharply after lag q (the signature from the autocorrelation dive). Because the shocks are unobserved, MA terms can't be fit by plain least squares — they need iterative maximum likelihood, which is the other half of what a library like statsmodels does for you.

Integration: difference, then undo

AR and MA both assume stationarity, so ARIMA sandwiches them between differencing and its inverse. The I(d) step differences the series d times to remove trend; the model is fit on that stationary version; then forecasts are integrated back — cumulatively summed — to return to the original scale. "Integrated" is just the undo of "differenced." For the airline series we already found d = 1 (plus a seasonal difference) does the job, so it's an ARIMA with d = 1 at its core.

Read the orders off the plots

The three orders aren't guessed — they're read from the diagnostics you already built. The d comes from the stationarity work: difference until the ADF test passes. Then p and q come from the ACF and PACF of the differenced series, using the exact signatures from the autocorrelation dive:

An ACF that decays with a PACF that cuts at lag 1 says ARIMA(1, d, 0); the mirror says ARIMA(0, d, 1). In practice you shortlist a few candidates this way and let an information criterion (AIC) pick between them — but the plots get you to the right neighbourhood.

It's just regression on lags

Fitting the AR part means stacking the series against shifted copies of itself and solving least squares. Laid out as a spreadsheet, column B is simply column A shifted down one row — yesterday's value lined up next to today's (a stationary AR(1) series, mean ≈ 0):

Iterate the equation forward

To forecast, you just run the fitted equation forward, feeding each prediction back in as the next step's input. For a stationary AR, this means forecasts decay smoothly toward the mean — an AR(1) forecast h steps out is φ^h times today's deviation, fading to the long-run average as the past loses its grip. MA terms contribute only for the first q steps (after that all the relevant shocks are forecast as zero). Finally, every forecast made on the differenced scale is integrated back by cumulative summing, returning it to the units you started in. Uncertainty grows with the horizon, which is why forecast bands fan out.

Both behaviours are visible at once on a plain stationary AR(1). Drag φ and watch the forecast (orange) decay from the last value back to the mean, while the band fans out — not forever, but toward a fixed width set by the series' own variance. Higher φ means a longer memory: slower decay, wider final band:

This is the signature of every stationary forecast: the point prediction reverts to the mean at rate φ^h, and the uncertainty grows but saturates — once you're far enough out, the best guess is simply the average and the band is the full spread of the series. (A non-stationary forecast, by contrast, never settles; that's exactly why we difference first.)

SARIMA for repeating patterns

Plain ARIMA has no notion of "the same month last year." SARIMA(p, d, q)(P, D, Q)_s bolts on a second set of terms operating at the seasonal lag s: a seasonal difference (D), seasonal AR terms (P) linking y_t to y_t−s, and seasonal MA terms (Q). The famous "airline model" — SARIMA(0,1,1)(0,1,1)₁₂ on the log series — uses one ordinary and one seasonal difference, plus one ordinary and one seasonal MA term, and fits the airline data beautifully with just a handful of parameters. That's the model behind the forecast below.

Two years into the future

Here is the airline model fit to all 144 months and forecast 24 months ahead, on the original scale. The trend continues, the seasonal peaks repeat, and the shaded band shows the 80% interval widening with the horizon — exactly the behaviour the pieces above predict:

The forecast isn't magic — it's the trend extrapolated, the seasonal shape repeated, and the recent shocks faded out, all read from the same decomposition and autocorrelation structure this lab has been pulling apart. And the widening band is honest: two years out, the model knows much less than it does next month.

From scratch, then statsmodels

The AR part is pure least squares — write it once to see there's no magic — and it matches statsmodels' estimator to machine precision:

# --- fit AR(p) from scratch: regress y_t on its p lags ---
def ar_fit(x, p):
    n = len(x)
    Y = x[p:]                                          # targets
    X = np.column_stack([np.ones(n-p)] +
                        [x[p-i-1 : n-i-1] for i in range(p)])   # const + lagged columns
    beta = np.linalg.lstsq(X, Y, rcond=None)[0]
    return beta                                        # [c, phi_1, ..., phi_p]

The full model — including the MA terms that need iterative maximum likelihood, and the seasonal structure — is one call to statsmodels:

from statsmodels.tsa.statespace.sarimax import SARIMAX
res = SARIMAX(np.log(y), order=(0,1,1),
              seasonal_order=(0,1,1,12)).fit()
forecast = np.exp(res.get_forecast(24).predicted_mean)   # back to original scale

from-scratch AR vs statsmodels AutoReg
AR(1):  phi_1 = 0.7604   (== lag-1 ACF)   max|diff vs statsmodels| = 1e-16
AR(2):  phi_1 = 0.8221 , phi_2 = -0.0842                max|diff| = 2e-16
SARIMA(0,1,1)(0,1,1)12 on log airline:  AIC = -483.4 ,  24-mo forecast peaks ~738

The model the whole lab was for

ARIMA is the assembly of everything prior: difference to stationarity (the I, from the stationarity dive), regress on lagged values (the AR, fit by least squares), and model the shock memory (the MA, whose order you read from the ACF). Choose d with the ADF test, p and q from the PACF and ACF, let AIC settle ties, add seasonal terms for repeating patterns, and forecast by iterating forward and integrating back. It's not the newest tool, but it's the one every forecaster should reach for first and the baseline every neural method is measured against.

In-sample fit lies

Any flexible model can be made to hug its training data — add enough AR and MA terms and the residuals vanish. That tells you nothing about forecasting, which is prediction of the unseen. The honest question is always: held out the last stretch of history, how close were the forecasts to what actually happened? Every number in this dive is computed on data the model never touched during fitting.

You cannot shuffle

In ordinary supervised learning you split data randomly. Do that to a time series and you commit a subtle but fatal sin: the model trains on June while testing on May, learning from the future to predict the past. That leakage inflates every score and collapses the moment you deploy. The only valid split respects the arrow of time — train on an initial segment, test on what comes strictly after. Any preprocessing (scaling, differencing parameters, feature statistics) must likewise be fit on the training portion alone.

Hold out the future

The simplest honest evaluation is a single time-ordered split: train on the first stretch, forecast the held-out tail, and compare to what really happened. For the airline data we train on the first 120 months (ten years) and forecast the final 24 — a genuine two-year-ahead test. The held-out portion is sacred: it's looked at exactly once, to score the final model, never to choose it.

Three ways to measure error

Given actuals y_t and forecasts f_t over a test set of size n, three metrics dominate:

MAE is the average error in the original units — interpretable and robust. RMSE squares first, so it punishes a few large misses far more than many small ones; it's the right choice when big errors are disproportionately costly. MAPE reports error as a percentage, making it scale-free and comparable across series — but it explodes when actuals near zero and penalises over- and under-prediction asymmetrically. Here is the per-month computation laid out as a spreadsheet (first 6 of the 24 test months, SARIMA forecasts):

Earn your complexity

A metric in isolation is meaningless — RMSE of 43 is only "good" relative to something. That something is a naive baseline. Two are standard: the naive forecast predicts the last observed value for everything ahead (a flat line), and the seasonal naive predicts the value from one full season ago (last January for this January). On a trending, seasonal series the seasonal naive is a strong opponent — it gets the shape right for free. Any model that can't beat it isn't earning its complexity. This is the time-series cousin of the dummy-classifier baseline from the ML lab.

Refit as you go

A single split uses the test set once; a walk-forward backtest wrings far more signal from it. Start at the split, forecast one step, then reveal the true value, slide the window forward to include it, refit, and forecast the next step — marching through the entire test period the way the model would actually be used in production. With an expanding window the training set grows each step; with a sliding window it stays a fixed length (better when the process drifts). Walk-forward gives many test forecasts instead of one batch, a more stable error estimate, and an honest picture of one-step-ahead performance — at the cost of refitting repeatedly.

Step through it. Drag the backtest forward: the training window (shaded) grows month by month, and at each step the model forecasts just one step ahead (orange) and is graded against the truth (cyan) before sliding on. Because every forecast is only a month out — the easiest horizon — the errors are tiny:

Note the punchline: this walk-forward one-step RMSE settles near 15 — a third of the 43 from the two-year batch forecast above. They're not contradictory; they answer different questions. "How well can I predict next month, refitting as I go?" is a far easier task than "how well can I predict the next two years from today?" Always report which horizon a number describes.

Grade them on the holdout

Now the payoff: forecast the 24-month holdout three ways and look. Toggle the model and watch the forecast track (or miss) the actuals, with each month's error drawn as a vertical drop:

And the same three models as a leaderboard — switch the metric and the ranking holds: the naive flat line is hopeless, the seasonal naive is a genuinely strong baseline, and SARIMA edges past it by learning the trend on top of the season:

SARIMA wins on every metric (MAE 39.5, RMSE 43.2, MAPE 8.5%), but notice how much it beats the seasonal naive (MAE 47.6) — modestly. That margin is the real value the model adds over a free baseline, and it's the honest story to report: not "8.5% error," but "8.5% versus the 10.5% you'd get for nothing."

Where, not just how much

A single number hides structure. Plot the forecast errors over time and interrogate them: do they drift (the model is biased, systematically high or low)? Do they swell toward the end (accuracy decays with horizon, as it must)? Are they autocorrelated — a big miss followed by another? Leftover autocorrelation in the errors means the model missed exploitable structure, which is exactly the Ljung–Box check from the autocorrelation dive applied to forecast residuals. Good forecast errors look like white noise centred on zero; anything else is a clue for the next iteration.

Ways to fool yourself

From scratch, then the loop

The metrics are three one-liners; the walk-forward backtest is a short loop that refits as it marches forward:

# --- metrics from scratch ---
mae  = np.mean(np.abs(y - f))
rmse = np.sqrt(np.mean((y - f)**2))
mape = np.mean(np.abs((y - f) / y)) * 100

# --- walk-forward backtest (expanding window, one-step-ahead) ---
errors = []
for t in range(split, len(series)):
    model = SARIMAX(series[:t], order=(0,1,1),
                    seasonal_order=(0,1,1,12)).fit(disp=False)   # train on the past only
    f = model.forecast(1)[0]                                     # predict the next step
    errors.append(series[t] - f)                                 # then reveal the truth and move on
rmse_wf = np.sqrt(np.mean(np.square(errors)))

airline holdout — train 120, test 24
naive (last value)   MAE 115.2   RMSE 137.3   MAPE 23.6%
seasonal naive       MAE  47.6   RMSE  50.0   MAPE 10.5%   <- strong baseline
SARIMA               MAE  39.5   RMSE  43.2   MAPE  8.5%   <- wins on every metric

Honest grading, start to finish

Evaluating a forecast is mostly about not cheating. Split on time so the test set is the future; fit every transform on the training portion; measure with MAE, RMSE, and MAPE while knowing what each rewards; always compare to a naive and seasonal-naive baseline; backtest walk-forward for a stable, deployment-realistic estimate; and inspect the errors for bias, horizon decay, and leftover autocorrelation. Report the margin over the baseline, not the raw metric. Do all that and the number you publish is one you can actually stand behind.

Recent data matters more

If you had to forecast tomorrow from a noisy history, you'd weight yesterday more than a year ago. Exponential smoothing makes that instinct precise: every forecast is a weighted average of all past values where the weights decay geometrically into the past. One smoothing parameter controls how fast. Layer one such average for the level, another for the trend, another for the season, and you've reconstructed the whole anatomy of a series — but as a live, updating forecaster rather than a static decomposition.

One level, updated

The simplest case — no trend, no season — tracks a single level ℓ_t, nudging it toward each new observation:

Each level is a blend of the newest value and the previous level; the forecast for every future step is just the last level (a flat line). By hand with α = 0.3 on the first airline values, starting ℓ₁ = 112: ℓ₂ = 0.3·118 + 0.7·112 = 113.8, then ℓ₃ = 0.3·132 + 0.7·113.8 = 119.26, and so on. Laid out as a spreadsheet, each level cell reads from the value beside it and the level above it:

How fast to forget

Unrolling the recursion shows what α really does: the weight on the value k steps back is α(1−α)^k — a geometric decay. Large α forgets quickly and hugs the latest point (responsive but jumpy); small α forgets slowly and approaches the long-run average (smooth but sluggish). Drag α and watch the smoothed level chase a noisy, stepping series, with the weight-decay shown alongside:

There's no universally right α — it's fit from the data by minimising forecast error. Notice the tension: crank α up and the level reacts instantly to the step change but jitters on the noise; turn it down and the noise vanishes but the level takes many steps to catch a real shift. That trade-off between responsiveness and smoothness is the whole game.

Add a slope

Simple smoothing forecasts a flat line, so it lags any trending series badly (try it on the airline data and it trails the climb forever). Holt's method fixes this by smoothing a second component — the trend b_t — with its own parameter β:

The level now updates around the previous level plus trend, and the trend is itself an exponentially-smoothed estimate of how much the level changes each step. The forecast is no longer flat but a straight line: ŷ_T+h = ℓ_T + h·b_T. That single addition turns a lagging tracker into something that extrapolates growth.

See the difference on a trending series. Simple smoothing (amber) forecasts a flat line from its last level — hopeless against a climb. Holt (orchid) carries the slope forward. Drag β to control how quickly the trend itself adapts:

The flat SES forecast falls further behind every step, while Holt's sloped forecast stays on the trend line. A higher β lets the slope chase recent changes (useful if growth is accelerating) at the cost of more sensitivity to noise — the same responsiveness-versus-smoothness trade-off as α, now for the trend.

Add the season

One more component and one more parameter captures the yearly pattern. Holt–Winters smooths a seasonal term s_t with rate γ, recycled every m = 12 months, and combines all three for the forecast. Seasonality comes in two flavours:

This is the same additive-vs-multiplicative choice from the anatomy dive, now baked into the forecaster. For the airline series the multiplicative form is the right call, because the summer peaks grow in proportion to the rising baseline.

Three smoothers, working together

The full multiplicative Holt–Winters is just three coupled exponential averages plus a forecast rule:

Each line is the same pattern — a fraction of the newest evidence plus the remainder of the old estimate — applied to the level, the per-step change, and the seasonal factor respectively. The level equation deseasonalises y_t by dividing out the season before updating; the season equation re-estimates the factor from how far y_t sits above or below the level.

Pick the smoothing rates

The parameters α, β, γ (each in [0,1]) and the initial level, trend, and seasonal states are all chosen at once to minimise the sum of squared one-step errors over the history — a numerical optimisation, just as ARIMA's coefficients are fit by maximum likelihood. On the airline series the fit lands near α = 0.32, β ≈ 0 (the trend is steady, so it barely updates), and γ = 0.60 (the seasonal shape drifts moderately). A near-zero β is informative: it says the growth rate is essentially constant across the decade.

Two families, similar power

ETS and ARIMA are the two pillars of classical forecasting, and they overlap more than they differ: many exponential-smoothing models have an exact ARIMA equivalent (simple smoothing is an ARIMA(0,1,1), for instance). The mental split is in framing — ETS thinks in explicit, interpretable components you can read off; ARIMA thinks in autocorrelation and differencing. In practice you try both and let out-of-sample error decide. On the airline holdout from the evaluation dive, Holt–Winters actually edges ahead:

Neither family is universally best — the lesson from the evaluation dive stands: don't assume, measure.

Holt–Winters on the airline

Here is the multiplicative Holt–Winters model fit to all 144 months and forecast 24 ahead, with the SARIMA forecast from the ARIMA dive overlaid for comparison. The two classical families trace nearly the same future — trend continued, seasonal peaks growing — from completely different machinery:

That the two methods agree so closely is reassuring, not redundant: when models built on different principles converge on the same forecast, you can trust it more. Where they diverge is where the forecast is genuinely uncertain.

From scratch, then statsmodels

Simple smoothing is a three-line loop; the level it produces matches statsmodels exactly:

# --- simple exponential smoothing from scratch ---
def ses(y, alpha):
    level = y[0]
    out = [level]
    for t in range(1, len(y)):
        level = alpha * y[t] + (1 - alpha) * level    # blend new obs with old level
        out.append(level)
    return out                                         # forecast = out[-1], held flat

The full Holt–Winters model — all three components, additive or multiplicative, fit by optimisation — is one statsmodels call:

from statsmodels.tsa.holtwinters import ExponentialSmoothing
hw = ExponentialSmoothing(y, trend='add', seasonal='mul',
                          seasonal_periods=12).fit()
forecast = hw.forecast(24)

from-scratch SES vs statsmodels (alpha=0.3)
final level = 138.388   (matches SimpleExpSmoothing exactly)
Holt-Winters fit:  alpha=0.319 , beta=0.000 , gamma=0.601
24-month forecast peaks ~705   ·   holdout MAE 29.0 (beats SARIMA's 39.5)

The components, made live

Exponential smoothing forecasts a series by maintaining a handful of exponentially-weighted components and projecting them forward. Start with a single smoothed level (simple exponential smoothing); add a smoothed trend for a sloped forecast (Holt); add a smoothed seasonal factor, additive or multiplicative, for repeating patterns (Holt–Winters). The smoothing rates α, β, γ set how fast each component forgets, and they're fit by minimising error. It's the interpretable cousin of ARIMA — different machinery, comparable accuracy — and on the airline data it wins. Two families, one lesson: build the components you can see, then let the held-out error pick the winner.

A series is a regression in disguise

Supervised learning needs a table: rows of features, each with a target. A time series doesn't look like that — until you build the table yourself. Take each month, use the few months before it as the input features and the month itself as the target, and you've converted one long sequence into hundreds of ordinary (x, y) rows. Fit any regressor to that table and you have a forecaster. This is how the same gradient boosting and random forests from the ML lab end up winning forecasting competitions.

The sliding window

The core trick is the lag feature: to predict y_t, use y_t−1, y_t−2, …, y_t−k as inputs. Slide a window of width k along the series and read off one row at each stop. Drag the window below and watch the design matrix assemble itself — each position contributes the highlighted lags as features and the next point as the target:

Stacking every window gives the feature matrix X and target vector y — a perfectly ordinary regression dataset (here with three lags on the airline series):

Beyond just lags

This is where the ML framing pulls decisively ahead of univariate ARIMA: once the series is a feature table, you can bolt on any column. Feature engineering becomes the main lever, and four families do most of the work:

Which features actually carry the forecast? Fit a model and read its feature importances. With lags [1, 2, 3, 12, 13] on the airline series, the answer is emphatic — and it depends on the model. Toggle between a random forest and a linear model:

The random forest pours 97% of its importance into the seasonal lag y_t−12 — for monthly data, last year's same month is overwhelmingly the best single predictor, and a tree happily keys on it alone. The linear model spreads its weight across the recent lag, the seasonal lag, and their interaction y_t−13, because it must combine them additively. Same data, same features, two very different stories about how the forecast is made — a reminder that importance is a property of the model, not just the series.

Forecasting more than one step

A lag model naturally predicts one step. For a longer horizon there are two strategies. Recursive: forecast one step, then feed that prediction back in as a lag to forecast the next, and so on — simple, but errors compound. Direct: train a separate model for each horizon (one for "+1 month", another for "+2 months", …) — no compounding, but more models and no coherence between them. Recursive is the common default and the one used for the forecast later in this dive.

From linear to boosted trees

Any regressor fits the table, and the choice echoes the ML lab. Linear regression on lag features is special: it's exactly an AR model — the same least-squares-on-lags fit from the ARIMA dive, just assembled by hand. Tree ensembles (random forests, gradient boosting) earn their keep when the dynamics are nonlinear or when many heterogeneous features interact — they capture thresholds and regime changes a linear model can't, which is why gradient boosting dominates tabular forecasting competitions with many related series.

That last point is the modern twist: the global model. Instead of one model per series, train a single regressor on the stacked feature tables of thousands of related series at once — every product in a catalogue, every sensor in a fleet. The model learns patterns that recur across series and lets data-rich series lend strength to sparse ones, something no per-series ARIMA can do. This is how gradient-boosted trees (LightGBM in particular) won the M5 forecasting competition, and it's ML forecasting's genuine structural edge: classical methods fit each series alone, while one global model exploits all of them together.

Why trees flatline on a trend

There's a catch that bites beginners. A tree can only predict averages of values it saw in training, so it cannot extrapolate beyond the range of the data. Point a raw gradient-boosting model at a rising series and its forecast flattens at the highest value it ever saw — useless for a trend. On the airline holdout this is stark: linear-on-lags scores RMSE 37.7, but the same-featured gradient booster scores a worse 54.8, precisely because it can't follow the climb. The fix is the lesson from the stationarity dive: difference or detrend first so the model predicts changes (which stay in range), or include an explicit trend feature and a model that can use it. Classical methods extrapolate trends for free; ML models must be set up to.

Watch it happen. The same lag features feed a linear model and a gradient booster, each forecasting 24 months past the end of the airline series. The linear forecast sails upward with the trend; the booster's forecast bends over and flatlines, unable to print a value above the dashed ceiling it saw in training:

The booster isn't broken — it's doing exactly what trees do, averaging training targets, none of which exceeded 622. Difference the series first and the problem vanishes: the trees then predict changes, which stay comfortably inside the range they trained on, and the cumulative sum rebuilds the trend. It's the stationarity dive's lesson with real teeth.

The same rules apply

Everything from the evaluation dive carries over, with teeth. You still must not shuffle: use a time-ordered split and, for tuning, rolling-origin cross-validation (expanding-window backtests) rather than random k-fold. Feature engineering is a new leakage hazard — a rolling mean that peeks at the target month, or scaling fit on the whole series, quietly leaks the future. Every feature for row t must be computable using only data available strictly before t. The flexibility of ML buys power and a dozen new ways to fool yourself.

Lag-features vs SARIMA

Here is a linear model on lag features [1, 2, 3, 12, 13] of the log series, forecast recursively 24 months and laid over the SARIMA forecast. Built from generic regression machinery — no time-series model at all — it lands almost on top of SARIMA:

On the holdout this lag-feature regression actually scores MAE 33.9 — a touch better than SARIMA's 39.5 — because with the seasonal lag y_t−12 in the feature set, the linear model is essentially a seasonal autoregression. It's a reminder that the boundary between "classical" and "ML" forecasting is blurry: a regression on the right lags is both.

When to reach for which

Build the table, fit any model

The entire ML approach is two steps: window the series into a matrix, then hand it to a regressor:

# --- turn a series into a supervised table ---
def make_Xy(s, lags):
    X, y = [], []
    for t in range(max(lags), len(s)):
        X.append([s[t - L] for L in lags])    # lag features for row t
        y.append(s[t])                         # the value to predict
    return np.array(X), np.array(y)

X, y = make_Xy(np.log(series), lags=[1, 2, 3, 12, 13])
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(X, y)           # linear-on-lags == an AR model

# --- recursive multi-step forecast ---
hist = list(np.log(series))
for _ in range(24):
    x = [[hist[-L] for L in lags]]
    hist.append(model.predict(x)[0])           # feed prediction back in
forecast = np.exp(hist[-24:])

airline holdout (train 120, test 24), lags [1,2,3,12,13] on log
LinearRegression       MAE 33.9   RMSE 37.7   <- beats SARIMA (39.5 / 43.2)
GradientBoosting       MAE 41.5   RMSE 54.8   <- can't extrapolate the trend
swap in any sklearn regressor: same two-step recipe

Forecasting is regression with the right features

The machine-learning route to forecasting is a single idea — slide a window to turn the series into a feature table — followed by the entire supervised toolbox. Use lag features (a seasonal lag is gold), add calendar and exogenous columns ARIMA can't easily take, pick a regressor (linear is an AR model; trees handle nonlinearity), and forecast recursively or directly. Mind the two traps: tree models can't extrapolate a trend, so difference first; and feature engineering must never peek across the time split. On one short series the classical methods hold their own — but give ML many series and rich features, and it pulls ahead. Either way, the discipline from the evaluation dive is what keeps the number honest.

A regression, with time as data

Ordinary regression explains y from predictors x. Time series regression does the same, but the predictors are things you can build for any date: a trend term (the row index t), seasonal indicators (which month it is), and any external variable you have a value for. The model is plain linear regression — y_t = β₀ + β₁t + (seasonal) + (drivers) + ε_t — and its coefficients are directly interpretable: "+1.0% per month of trend," "+21% in July." That interpretability, and the ability to fold in drivers, is why it's the workhorse of econometric forecasting.

Trend plus season

The simplest useful version regresses (log) airline traffic on a linear trend and a set of monthly indicators. Toggle the features and watch the fit: trend alone is a sloped line that captures the growth but misses every seasonal swing; add the seasonal terms and the fitted curve snaps onto the data:

Trend-only explains 90% of the variance (R² = 0.90) but leaves the whole seasonal pattern in the residuals; adding eleven monthly indicators lifts it to R² = 0.98 and the fitted line tracks each peak and trough. The fitted trend coefficient is +1.0% per month — about 12.8% annual growth — and the July indicator sits highest, the November lowest, recovering the seasonal shape as plain regression coefficients you can read off and report.

The deterministic features

Two feature types do the work here. The trend is just the time index t = 0, 1, 2, … (or a polynomial, or a piecewise spline if growth bends). The season is encoded as dummy variables: for monthly data, eleven 0/1 columns (one month is the baseline, absorbed into the intercept), each coefficient measuring that month's average offset from baseline after trend. Unlike the differencing in ARIMA, these are deterministic — they're known for any future date, so forecasting just means evaluating the formula forward. The cost: eleven parameters for the season alone, which the next idea trims dramatically.

Seasonality in a few waves

Eleven dummies is a lot of parameters. Fourier terms express the same seasonality as a sum of sine and cosine waves — each "harmonic" is one sin/cos pair at a multiple of the seasonal frequency. A handful of harmonics reconstruct a smooth seasonal shape with far fewer parameters. Drag the number of harmonics and watch them build up the monthly pattern:

With K harmonics you spend just 2K parameters instead of eleven, and for the airline season two or three harmonics already capture R² ≈ 0.97–0.98 of what the full dummy set achieves. Here is the design matrix — the deterministic features computed straight from the time index t:

The real power: outside drivers

Trend and season you could get from ARIMA. What time series regression adds for free is exogenous predictors — any external series you believe drives the target. Retail sales regressed on advertising spend and a price index; electricity demand on temperature; bookings on a holiday calendar. Each is just another column in the feature matrix, and its coefficient is the estimated effect of that driver, holding the others fixed. This is the bridge to causal-flavoured questions ("how much does a 1° rise in temperature add to demand?") that pure autoregression can't answer. The price: to forecast, you now need future values of the drivers — known for a calendar, but themselves a forecasting problem for things like temperature.

Why plain OLS isn't enough

Here's the catch that makes time series regression its own topic. Ordinary least squares assumes the errors are independent — but in a time series, a residual above the line is almost always followed by another above the line. Plot the autocorrelation of the trend-plus-season residuals and the violation is blatant:

The OLS residuals have a lag-1 autocorrelation of 0.78 and a Durbin–Watson statistic of 0.43 (2.0 means no autocorrelation) — the errors carry enormous leftover structure. The coefficients themselves stay unbiased, but their standard errors are badly wrong, so every p-value and confidence interval the regression reports is untrustworthy, and the forecasts ignore exploitable signal. You cannot stop here.

Model the leftover structure

The fix is elegant: keep the regression, but stop pretending the error is white noise. Regression with ARIMA errors writes y_t = (regression on features) + η_t, where the error η_t follows an ARIMA process, and fits both parts together:

Toggle the residual plot above to "+ ARIMA errors": the autocorrelation collapses inside the bands — the ARIMA part has absorbed exactly the structure plain OLS left behind. Now the standard errors are valid, inference is honest, and the forecasts use the error dynamics too. This combination — interpretable regression coefficients plus a properly-modelled stochastic error — is dynamic regression, and it's the right way to put external predictors into a forecast.

The practical recipe

Put the pieces together for a seasonal series and you get dynamic harmonic regression: Fourier terms for a compact deterministic season, plus ARIMA errors (with differencing) to handle trend and leftover autocorrelation. For the airline series, Fourier(K=3) regressors with ARIMA(1,1,1)(1,0,1)₁₂ errors leaves residuals that pass the Ljung–Box test (p = 0.79 — white) and forecasts a two-year peak near 715 passengers — squarely among SARIMA (738), Holt–Winters (705), and the ML lag model (702). Four methods, four philosophies, one consensus future. It's especially handy when the season is long (weekly data, period 52) where a SARIMA would need an unwieldy number of seasonal terms but a few Fourier harmonics suffice.

Seeing is believing. Here are all four forecasting families from this lab — dynamic regression, SARIMA, Holt–Winters, and the ML lag model — each fit to the airline series and projected two years ahead. Built on four different philosophies, they trace nearly the same future:

The four lines bunch together through the forecast — that agreement is the signal, and it should raise your confidence in the consensus shape. Where they fan apart, at the far summer peaks two years out, is exactly where the future is genuinely uncertain and no method can rescue you. Four philosophies converging is worth more than any single model's claim.

Four ways to forecast, compared

Dynamic regression is the one to reach for when you have drivers — variables outside the series that explain it — and when you need to report the effect of each, not just a forecast. For a pure self-driven series, ARIMA or ETS is simpler.

OLS, then ARIMA errors

The linear model is one matrix solve; the upgrade to ARIMA errors is one statsmodels call with an exog argument:

import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX

t = np.arange(len(y))
def fourier(t, K, P=12):
    cols = [f(2*np.pi*k*t/P) for k in range(1, K+1) for f in (np.sin, np.cos)]
    return np.column_stack(cols)

X = np.column_stack([np.ones_like(t), t, fourier(t, 3)])   # const + trend + season
beta = np.linalg.lstsq(X, np.log(y), rcond=None)[0]         # plain OLS (biased SEs!)

# --- dynamic regression: same features, ARIMA error ---
m = SARIMAX(np.log(y), exog=fourier(t, 3),
            order=(1,1,1), seasonal_order=(1,0,1,12)).fit()
future = fourier(np.arange(len(y), len(y)+24), 3)
forecast = np.exp(m.get_forecast(24, exog=future).predicted_mean)

airline · trend + seasonal dummies (OLS)
R2 = 0.983 , trend = +1.0%/mo (12.8%/yr)
residual autocorrelation: lag-1 = 0.78 , Durbin-Watson = 0.43   <- errors NOT independent
+ ARIMA(1,1,1)(1,0,1)12 errors (dynamic harmonic regression)
residual Ljung-Box p = 0.79  -> white   ·   24-month forecast peak ~715

Explain the series, then model the rest

Time series regression flips the forecasting question from "what does the past imply?" to "what explains this?" Build deterministic features — a trend term, seasonal dummies or compact Fourier harmonics — and add any external drivers you have, then fit a linear model whose coefficients you can interpret and report. But never trust plain OLS on a time series: its residuals are autocorrelated, its standard errors are wrong. Upgrade to regression with ARIMA errors, which keeps the interpretable regression while modelling the leftover structure, giving honest inference and forecasts that exploit both the drivers and the error dynamics. When you have outside variables that move your series — and you want to know how much — this is the tool.

Learnable nonlinear features

Linear regression on lags can only add up its inputs; a tree can only split them. A neural network instead passes the inputs through layers of weighted sums and nonlinear "activations," learning its own intermediate features as it trains. With enough neurons it can represent essentially any smooth relationship between the recent past and the next value. That flexibility is the appeal — and the danger: it needs a lot of data to pin down all those weights, and on a short series it will happily memorise noise.

One unit, from scratch

The building block is a single neuron: multiply each input by a weight, add a bias, and pass the total through a nonlinear activation. That's the whole computation — here it is by hand for inputs [1.0, 0.5, −0.2] with weights [0.6, −0.3, 0.2] and bias 0.1:

From inputs to prediction, by hand

A network is just neurons in layers. Push three inputs x = [1.0, 0.5, −0.2] through a two-neuron hidden layer (tanh), then one output neuron. Each hidden neuron does its own weighted sum of all three inputs, adds its bias, and applies tanh — every number by hand:

In matrix form it's two lines, and the from-scratch code prints exactly the by-hand number — this is what every neural-net library computes internally:

import numpy as np
x  = np.array([1.0, 0.5, -0.2])                 # input (scaled lags)
W1 = np.array([[0.5,-0.4,0.3],[-0.2,0.6,0.5]])  # hidden weights (2x3)
b1 = np.array([0.1, -0.2])
W2 = np.array([0.7, -0.5]);  b2 = 0.2           # output weights

h    = np.tanh(W1 @ x + b1)    # hidden layer:  [0.3275, -0.1974]
yhat = W2 @ h + b2             # prediction out

from-scratch forward pass
h    = [0.3275, -0.1974]
yhat = 0.52792          <- matches the by-hand value exactly

Why stacking neurons works

One neuron draws one soft bend. Put many in a layer and sum their outputs, and you can approximate any continuous function — the more neurons, the finer the detail. Drag the hidden-unit count and watch a network fit a wiggly target: one neuron manages only a near-flat line, while sixteen trace every curve (and start chasing noise):

This is the universal approximation property: a single hidden layer, wide enough, can represent any continuous function to arbitrary accuracy. It's why neural networks are so expressive — and why capacity control (regularisation, early stopping, enough data) matters so much. Too few neurons underfit; too many, with too little data, overfit. The whole art is finding the middle.

Gradient descent and backprop

So far the weights were given. Training finds them by shrinking a loss. Say the true value was 1.0 but the forward pass predicted 0.528 — the squared-error loss is L = ½(0.528 − 1)² = 0.111. Backpropagation sends that error backward through the layers by the chain rule, getting each weight's gradient — how much nudging it moves the loss:

Everything hinges on the step size η. Drag it on this loss bowl (the loss as one weight varies): too small and training crawls; near the critical value it overshoots and bounces; too large and it diverges entirely:

This is why learning rate is the most fiddled-with knob in deep learning. The same logic scales from this one weight to the millions in a real network — backprop computes every gradient in one backward sweep, and an optimiser (SGD, Adam) takes the step.

NNAR: the simplest forecaster

The most direct way to forecast with a network needs nothing new: take the very same lag-feature table from the ML dive — y_t−1, y_t−2, … as inputs, y_t as target — and fit a feed-forward network instead of a linear model. This is neural network autoregression (NNAR). It inherits everything from the ML dive: the sliding window, recursive multi-step forecasting, and — sharply — the extrapolation problem. Fed a trending series and its own predictions, an unconstrained network doesn't just flatline like a tree; it can spiral to absurd values. The cure is the same: difference or seasonally anchor the series first so the network predicts something that stays in range.

Windows in, windows out

The architectures built for sequences usually skip hand-made lag columns and consume the raw window. You slide a fixed input window of the last L steps and ask the network to emit an output window of the next H steps — a "sequence-to-sequence" setup. An encoder compresses the input window into a summary vector; a decoder expands that into the forecast horizon. This framing forecasts all H steps at once (no recursive error compounding) and generalises cleanly to multivariate input (many series in) and probabilistic output (a distribution, not a point).

Memory that flows through time

A recurrent neural network (RNN) reads the sequence one step at a time, maintaining a hidden state — a memory vector — that it updates at each step from the previous state and the new input. Step through the unrolled network below: the same cell fires at every timestep, passing its memory rightward:

In principle the hidden state can carry information across the whole sequence, but plain RNNs struggle: gradients shrink as they propagate back through many steps (the vanishing gradient problem), so long-range memory fades. The fix is gated cells — the LSTM and GRU — which add learnable "gates" that decide what to keep, what to forget, and what to output. Those gates let a path for the gradient stay open across hundreds of steps, and they were the workhorse of sequence modelling for years.

Scanning with a wider and wider eye

A different idea borrows from image models: slide a small convolutional filter along the time axis, detecting local patterns (a spike, a ramp). To see far back without huge filters, temporal convolutional networks (TCNs) stack dilated causal convolutions — each layer skips exponentially more steps (1, 2, 4, 8…), so the receptive field doubles per layer and reaches deep into the past with few parameters. "Causal" means each output depends only on the present and past, never the future. TCNs train fast (every position computed in parallel, unlike an RNN's sequential pass) and often match or beat recurrent nets on long sequences.

Look anywhere at once

The dominant modern idea drops recurrence entirely. Attention lets the model, when predicting a step, weight every past position by how relevant it is and read a blend of them directly — no information has to squeeze through a single hidden state. For seasonal data the weights light up exactly where you'd hope: the recent points and the same month a year ago. The bars below sketch a plausible attention pattern over the recent history:

And the weights are real arithmetic, not magic. For a query vector q and three past keys k, attention scores each by the dot product q·k, scales by √d, exponentiates, and normalises (softmax) — then blends the values by those weights:

Stacks of attention layers make a Transformer, the architecture behind large language models and now the state of the art for long-horizon forecasting (Informer, Autoformer, PatchTST, and relatives). Attention's strength is reaching arbitrarily far back in one hop; its cost is computation that grows with the square of the sequence length, which the forecasting variants work hard to tame.

Forecasting's pretrained era

The newest turn mirrors language models: train one enormous network on millions of time series from every domain, then forecast a series it has never seen — zero-shot, with no fitting at all. These time-series foundation models (Chronos, TimesFM, Lag-Llama, Moirai, TimeGPT) tokenise a series much like text and predict the continuation. They're remarkable when you have little history for the target but the pattern resembles something in their vast training set, and they hand you calibrated probabilistic forecasts out of the box. Whether they consistently beat a well-tuned classical model on ordinary business series is still an open, actively-debated question.

When neural nets help — and when they don't

Here is the result that surprises newcomers. Neural networks are data-hungry: they have thousands of weights to estimate, and a single series of 144 points starves them. Across the famous M3 and M4 forecasting competitions, simple statistical methods and their ensembles repeatedly beat standalone deep networks on individual series. Deep learning pulls decisively ahead in a different regime — many related series (thousands of products, sensors, or households), long histories, rich exogenous inputs, or complex nonlinear and multivariate structure. The M5 competition, run on such data, was won by gradient boosting with neural nets close behind. Match the tool to the regime, not to the hype.

NNAR on the airline — a cautionary tale

Proof on our running series. A feed-forward NNAR (lags 1, 12, 13; one small hidden layer), forecast two years ahead, against the strong classical model:

The network's forecast is plausible in shape but visibly undershoots the trend, and its holdout error is more than double SARIMA's. With only 144 points it cannot learn the growth and seasonality reliably — exactly the data-starvation lesson, made concrete. This isn't a failure of neural networks; it's a failure to match method to data size. Hand the same architecture millions of points across thousands of series and the verdict flips.

A network is a few lines

The feed-forward NNAR is one estimator swap from the ML dive; the sequence models live in deep-learning frameworks:

# --- from scratch: forward pass + one backprop step ---
h    = np.tanh(W1 @ x + b1);  yhat = W2 @ h + b2   # forward (data in -> out)
t, lr = 1.0, 0.1
dy  = yhat - t                  # dL/dyhat
W2 -= lr * dy * h               # output-layer gradient + step
dz  = (dy * W2) * (1 - h**2)    # backprop through tanh
W1 -= lr * np.outer(dz, x)      # hidden-layer gradient + step
# loss 0.111 -> 0.061 after this single step

In practice you reach for a library, which adds optimisers, batching, and GPU support — the feed-forward NNAR is one estimator swap from the ML dive:

from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

X, y_t = make_Xy(np.log(series), lags=[1, 12, 13])   # same table as the ML dive
net = make_pipeline(StandardScaler(),                 # networks need scaled inputs!
                    MLPRegressor(hidden_layer_sizes=(8,), activation='tanh',
                                 alpha=0.1, max_iter=6000))
net.fit(X, y_t)                                       # then recursive forecast as before

# RNN / LSTM / Transformer live in torch or keras, e.g.:
# model = nn.LSTM(input_size=1, hidden_size=32); ... train on (window_in, window_out) pairs

airline holdout (train 120, test 24)
NNAR (MLP, lags 1/12/13)   MAE 93    <- data-starved on 144 points
SARIMA                     MAE 39    ·  Holt-Winters  29  ·  lag regression  34
verdict: on one short series, classical wins; NNs need scale

Powerful, hungry, situational

Neural networks bring genuine new capability to forecasting: learnable nonlinear features (the neuron and universal approximation), and architectures purpose-built for sequences — recurrent nets and LSTMs that carry memory, dilated convolutions that widen their view cheaply, attention and Transformers that reach anywhere in the past, and foundation models pretrained across millions of series for zero-shot use. But all of it rides on data. The same lag-feature framing from the ML dive gets you an NNAR in one line; the sequence models need a deep-learning framework and far more data. On a single short series the classical methods usually win, cheaply and interpretably — so the skill is knowing which regime you're in.

One series, seven steps

A forecasting project isn't one model — it's a pipeline. Each stage is one of the dives in this lab, and the output of each feeds the next:

The discipline is to go in order. Skip the decomposition and you'll pick the wrong seasonal period; skip stationarity and the model misbehaves; skip evaluation and you'll ship a number you can't defend. Let's walk it.

Look before you model

Always plot first. The classical decomposition splits CO₂ into trend, season, and residual — and immediately tells us the structure to model:

The trend climbs almost linearly from ~354 to ~371 ppm; the seasonal component repeats a fixed ±3.5 ppm swing (CO₂ dips each northern-hemisphere summer as plants grow, peaks in spring); and the residual is small and structureless. Crucially, the seasonal band doesn't grow with the level the way the airline's did — so here we use additive seasonality and skip the log transform. The anatomy dive's first question — additive or multiplicative? — already changes our recipe.

Read the orders off the plots

The ACF of the raw series decays with a huge spike at lag 12 — the textbook signature of trend plus annual seasonality, and a loud demand to difference. After differencing, what survives at the seasonal lags tells us how many AR and MA terms, ordinary and seasonal, to shortlist. As in the autocorrelation dive: a PACF cutting off suggests the AR order, an ACF cutting off the MA order. For CO₂ the diagnostics point to one ordinary and one seasonal difference, with a small number of MA terms — candidates like SARIMA(1,1,1)(0,1,1)₁₂.

Difference until the test passes

The ADF test makes "difference how much?" precise. Climbing the same ladder from the stationarity dive:

Exactly as with the airline: neither difference alone suffices, but a regular and a seasonal difference together drive p below 0.05. That fixes d = 1 and D = 1 in the seasonal model — and tells the ML approach to predict the doubly-differenced series so its trees don't have to extrapolate.

Three rivals and a baseline

With the structure known, we fit the three families this lab built, plus the baseline they must beat:

Grade them on a holdout

Train on the first 120 months, forecast the last 24, and score against the truth — a time-ordered split, never shuffled. Switch the metric; the story holds:

Two lessons leap out. First, every real model demolishes the seasonal-naive baseline (MAE ~1.3 ppm) — sub-third-of-a-ppm errors, because CO₂ is extraordinarily regular. Second, the three sophisticated methods finish in a tight cluster, with the lag-feature regression nosing ahead (MAE 0.25) over Holt–Winters (0.28) and SARIMA (0.37). On a clean, strongly-structured series like this, the choice between good models barely matters — which is itself worth knowing before you agonise over it.

Two years of CO₂

Refit the chosen model on all 144 months and forecast 24 ahead. The trend continues, the annual cycle repeats at constant amplitude, and the band is tight because the series is so predictable:

The forecast says what we'd expect physically — CO₂ keeps rising, cresting near 377 ppm in the spring of the second year — and the narrow band reflects genuine confidence, not optimism: this series simply has very little noise to be uncertain about.

Did we miss anything?

The final discipline from the evaluation dive: interrogate the residuals. If structure remains in the model's errors, there's signal left to capture. Running the Ljung–Box test on the fitted model's residuals gives Q(12) = 5.9, p = 0.92 — comfortably failing to reject whiteness. The errors are indistinguishable from noise: the model has wrung out the trend, the season, and the short-term autocorrelation, leaving nothing exploitable behind. That's the signal to stop — the forecast is as good as this series allows.

What each step decided

The whole project, as a paper trail — each dive contributing one decision:

The pipeline, condensed

import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox

# 1 decompose — additive? then no log
seasonal_decompose(y, model='additive', period=12)
# 2-3 how many differences? climb until ADF passes
adfuller(np.diff(np.diff(y, 12)))            # p = 0.000 -> d=1, D=1
# 4-5 fit candidates, score on a time-ordered holdout
m = SARIMAX(y[:120], order=(1,1,1), seasonal_order=(0,1,1,12)).fit()
mae = np.mean(np.abs(y[120:] - m.forecast(24)))
# 6 refit on everything and forecast
final = SARIMAX(y, order=(1,1,1), seasonal_order=(0,1,1,12)).fit()
forecast = final.get_forecast(24)
# 7 are the residuals white?
acorr_ljungbox(final.resid, lags=[12])       # p = 0.92 -> done

CO2 holdout (train 120, test 24) — MAE in ppm
seasonal naive  1.31    SARIMA  0.37    Holt-Winters  0.28    lag-features  0.25
final forecast peaks ~377 ppm   ·   residual Ljung-Box p = 0.92 (white)

What you can now do

You've taken an unfamiliar series from raw numbers to a defensible, validated forecast using nothing but the tools in this lab. Decompose to see the structure and choose additive or multiplicative. Read the ACF and PACF to identify orders. Use the ADF test to decide how much to difference. Fit the classical families — ARIMA and exponential smoothing — and the ML reframing, all on the same footing. Split on time, score against a baseline with the right metric, and report the margin. Forecast the winner, then prove it's done by checking the residuals are white. That sequence is the entire discipline of classical forecasting, and it transfers to any series you'll meet.