Statistical forecasting: from baselines to ETS, SARIMA, and Theta

Statistical forecasting models are small. A simple exponential smoothing model has one parameter you tune and one number it remembers. A seasonal ARIMA model that fits a famous airline dataset well has two coefficients. These models are not glamorous, yet they remain the first thing to try on a new series. They train in milliseconds, they produce honest prediction intervals, and when a larger model loses to them in a backtest, you have learned something important about your data.

This guide builds the classical toolbox from first principles: baselines, exponential smoothing and its state-space form (ETS), stationarity and the ARIMA family, the Theta method, and Croston's method for intermittent demand. Each idea starts with what the symbols mean and a tiny numeric example before the general equation. Then we fit all of them with Nixtla's statsforecast on the AirPassengers series, check them with rolling-origin backtests, and hand-write the smoothing recursions in NumPy so nothing stays a black box. If forecast origins, rolling backtests, or MASE are new to you, read the mathematical intuition guide first. This guide assumes those ideas and goes deeper into the models.

The running example and four honest baselines

AirPassengers records monthly totals of international airline passengers, in thousands, from January 1949 to December 1960: 144 observations. It has a clear upward trend, a strong 12-month seasonal pattern with peaks in July and August, and one more feature that will matter a great deal later: the seasonal swings get larger as the level rises. In 1949 the gap between the busiest and quietest month is about 50 thousand passengers; by 1960 it is over 200 thousand.

Throughout the series of guides we hold out the last 12 months (1960) as a test set, forecast horizon $h = 12$, and season length $m = 12$. We report MAE, RMSE and MASE, where MASE divides the test MAE by the in-sample MAE of a seasonal naive forecast. A MASE below 1 means the model beat "repeat last year" on average over the test year, and above 1 means it lost.

Write $y_1, \dots, y_T$ for the training observations and $\hat{y}_{T+h \mid T}$ for the forecast of $y_{T+h}$ made with data up to time $T$. Four baselines need no fitting beyond simple arithmetic.

Naive. Tomorrow looks like today:

Seasonal naive. This month looks like the same month last year. With $k = \lfloor (h-1)/m \rfloor$ counting how many full seasons ahead we are:

For $h \le m$ this is simply $y_{T+h-m}$: the January 1960 forecast is the January 1959 value.

Mean. Everything is noise around one constant:

Drift. Draw a straight line from the first observation to the last and extend it:

A tiny example: if $y_1 = 100$, $y_T = 160$ and $T = 61$, the average step is $60/60 = 1$, so the 12-step drift forecast is $172$. The drift method is a random walk whose average step is estimated from the whole history.

Each baseline encodes one assumption: no change (naive), repeat the season (seasonal naive), no structure (mean), constant growth (drift). Any model you fit should be compared against the one whose assumption best fits your series. For monthly airline traffic, that is seasonal naive.

Simple exponential smoothing and Holt's trend

The naive forecast trusts only the last observation; the mean trusts all observations equally. Simple exponential smoothing (SES) sits between them. It keeps a running level $\ell_t$, an estimate of where the series currently sits, and updates it with each new observation:

Rearranged, the same update reads $\ell_t = \ell_{t-1} + \alpha(y_t - \ell_{t-1})$: move the level a fraction $\alpha$ of the way toward the surprise. The quantity $y_t - \ell_{t-1}$ is the one-step forecast error, so SES is an error-correction rule.

Deriving the geometric weights

Substitute the update into itself. One step back:

Keep going until you reach the initial level $\ell_0$:

The weight on an observation $j$ steps old is $\alpha(1-\alpha)^j$, a geometric sequence. Do these weights form a proper average? The geometric sum gives $\sum_{j=0}^{t-1}\alpha(1-\alpha)^j = 1-(1-\alpha)^t$, and adding the weight $(1-\alpha)^t$ on $\ell_0$ gives exactly 1. So the level is a weighted average of all past data and the starting value, with weights that decay exponentially with age. That is where the name comes from.

With $\alpha = 0.3$ the weights on the most recent six observations are $0.3, 0.21, 0.147, 0.103, 0.072, 0.050$, which together account for about $0.88$ of the total. With $\alpha = 0.9$ the newest point alone gets $0.9$ and the model behaves almost like naive. With $\alpha = 0.05$ the memory is long and the model behaves more like a slowly moving mean. Fitting SES means choosing $\alpha$ (and $\ell_0$) to minimise the sum of squared one-step errors, or equivalently maximising a likelihood, which we will meet in the ETS section.

The forecast is flat: every horizon gets $\ell_T$. SES has no notion of direction, which is fine for a series that wanders without trend and wrong for one that grows.

Holt's linear trend

Holt added a second state, the trend (slope) $b_t$, an estimate of how much the level changes per period:

Read the level equation as: blend the new observation with where last period's level plus slope said we would be. Read the trend equation as SES applied to the observed level changes. Example: $\ell_{t-1} = 200$, $b_{t-1} = 2$, $y_t = 206$, $\alpha = 0.5$, $\beta = 0.2$. The prediction was $202$, so $\ell_t = 0.5 \cdot 206 + 0.5 \cdot 202 = 204$. The level rose by $4$, so $b_t = 0.2 \cdot 4 + 0.8 \cdot 2 = 2.4$. The model nudged its growth estimate up because growth surprised it.

Damped trend

Extending a straight line 24 or 36 steps ahead usually overshoots, because real growth rarely continues at a constant rate forever. A damping parameter $0 \lt \phi \lt 1$ shrinks each future increment:

As $h \to \infty$ the sum approaches $\phi/(1-\phi)$, so forecasts flatten to $\ell_T + b_T\,\phi/(1-\phi)$. With $\phi = 0.9$ the long-run extra growth is at most $9 b_T$. The state recursions also gain a $\phi$: $\ell_t = \alpha y_t + (1-\alpha)(\ell_{t-1} + \phi b_{t-1})$ and $b_t = \beta(\ell_t - \ell_{t-1}) + (1-\beta)\phi\, b_{t-1}$. Values of $\phi$ between $0.8$ and $0.98$ are typical; $\phi = 1$ recovers Holt.

Holt-Winters: additive and multiplicative seasonality

Holt-Winters adds a third state: a seasonal index $s_t$ for each position in the cycle. With 12 months there are 12 seasonal values, and each one is updated only when its month comes around.

Additive seasonality

The additive version says the July effect is a fixed number of passengers added on top of the level:

The level update first removes last year's seasonal effect for this month ($y_t - s_{t-m}$, a seasonally adjusted observation), then smooths as before. The seasonal update compares the observation with the non-seasonal prediction and smooths that deviation into the month's index. The usual constraint is $0 \le \gamma \le 1 - \alpha$.

Multiplicative seasonality

The multiplicative version says July is a fixed percentage above the level, so the seasonal swing scales with the level:

Here the seasonal indices hover around 1: a July index of $1.22$ means July runs 22% above the deseasonalised level. The trend equation is unchanged.

Which one fits AirPassengers? The seasonal swing grows from about 50 to over 200 thousand while the level roughly quadruples, so the effect is proportional and multiplicative seasonality is the natural choice. There is a second route: take logarithms. Since $\log(\ell \cdot s) = \log \ell + \log s$, a multiplicative pattern becomes additive on the log scale, and any additive model can be used there. We will try both routes.

ETS: the state-space view and model selection

The smoothing recursions above are methods: rules that produce point forecasts. They become models once we say how randomness enters, and that step unlocks likelihoods, automatic model selection, and prediction intervals.

The simplest case

Let $\varepsilon_t$ be independent noise with mean zero and variance $\sigma^2$. The model ETS(A,N,N) reads:

Since $\varepsilon_t = y_t - \ell_{t-1}$ is the one-step error, the state equation is exactly SES in error-correction form. The same single noise term drives both the observation and the state, which is why these are called innovations state-space models.

The multiplicative-error version ETS(M,N,N) uses relative errors: $y_t = \ell_{t-1}(1 + \varepsilon_t)$ and $\ell_t = \ell_{t-1}(1 + \alpha\varepsilon_t)$. Its point forecasts are identical to ETS(A,N,N) but its noise grows with the level, so its intervals differ.

The taxonomy

ETS stands for Error, Trend, Seasonal. Each letter takes a value:

So SES is ETS(A,N,N) or ETS(M,N,N), Holt is ETS(A,A,N), damped Holt is ETS(A,Ad,N), and additive Holt-Winters is ETS(A,A,A). In general, with a state vector $x_t$ holding level, trend and seasonal values:

where $w$ maps states to a one-step prediction, $f$ says how states evolve without noise, and $r$, $g$ say how noise enters. For additive-error models $r = 1$; for multiplicative errors $r(x_{t-1}) = w(x_{t-1})$.

flowchart LR
  S0["State at t-1<br/>level, trend, seasons"] --> W["Predict y_t from states"]
  W --> Y["Observe y_t"]
  Y --> E["Innovation<br/>error = actual minus prediction"]
  E --> U["Update states<br/>using alpha, beta, gamma"]
  S0 --> U
  U --> S1["State at t"]
  S1 -.->|"next period"| S0

Likelihood and AICc

Because the model specifies a distribution for each one-step error, the parameters (smoothing weights, $\phi$, and the initial states) can be estimated by maximum likelihood. For additive Gaussian errors, maximising the likelihood is equivalent to minimising the sum of squared one-step errors.

To choose among the candidate ETS models, fit each one and compare the corrected Akaike information criterion:

where $L$ is the maximised likelihood, $k$ the number of estimated parameters including initial states, and $n$ the number of observations. AIC rewards fit and charges a price per parameter; the correction term matters when $n$ is small relative to $k$. A seasonal model with $m = 12$ has around 15 or more parameters, so the correction is not negligible on 132 points. AutoETS performs exactly this search. Two warnings: AICc is an in-sample criterion, so it can prefer a model that forecasts poorly, and AICc values are only comparable between models fitted to the same data on the same scale.

Stationarity, differencing and autocorrelation

ARIMA models take a different view. Instead of describing level, trend and season as states, they describe how each value correlates with its own past. That only makes sense if those correlations are stable over time.

A series is (weakly) stationary if its mean, its variance, and its covariance between values $k$ steps apart do not depend on $t$. AirPassengers is not: its mean rises, its variance rises, and the seasonal pattern repeats. White noise, independent draws with constant mean and variance, is the simplest stationary series.

Making a series stationary

Three standard moves, applied in order:

1. Stabilise the variance. When fluctuations scale with the level, take $\log y_t$. For AirPassengers this turns growing seasonal swings into roughly constant ones.
2. Remove the trend with a first difference: $y'_t = y_t - y_{t-1}$. The differenced series records monthly changes rather than levels.
3. Remove the season with a seasonal difference: $y_t - y_{t-m}$, the change from the same month last year.

ACF and PACF

The sample autocorrelation at lag $k$ measures how strongly the series correlates with itself shifted by $k$:

The partial autocorrelation at lag $k$ is the correlation between $y_t$ and $y_{t-k}$ after removing the part explained by the lags in between. A useful mental picture: ACF counts all paths of influence from $k$ steps back; PACF counts only the direct one.

Here are actual values on the log training data (computed with the formula above):

Slowly decaying correlations near 1 are the signature of a trend: the series is not stationary. After one difference the trend is gone, but lag 12 at $0.83$ shows the season is still there. After both differences the remaining structure is a negative spike at lag 1 and another at lag 12, with little else. Keep that pattern in mind; it points directly at a model.

Rules of thumb for identifying orders on a stationary series: an AR(p) process has a PACF that cuts off after lag $p$ and an ACF that decays gradually; an MA(q) process has an ACF that cuts off after lag $q$ and a PACF that decays. Real data rarely looks this clean, which is why automatic search exists.

AR, MA, ARIMA and SARIMA

Backshift notation

The backshift operator $B$ moves a series one step into the past: $B y_t = y_{t-1}$, and $B^k y_t = y_{t-k}$. It turns differencing into algebra: the first difference is $(1-B)y_t$, a second difference is $(1-B)^2 y_t$, and a seasonal difference is $(1-B^m)y_t$. Operators multiply like polynomials, which is what makes the seasonal model compact.

Autoregressive and moving-average pieces

An AR(p) model regresses the value on its own past $p$ values:

with $\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$. The intuition guide works through an AR(1) where a shock decays by a factor $\phi$ each step.

An MA(q) model uses past forecast errors, not past values:

Example with MA(1), $\mu = 50$, $\theta_1 = -0.4$: if last month's surprise was $\varepsilon_{t-1} = +10$, this month's expected value is $50 - 4 = 46$. Part of a surprise is reversed one step later and then forgotten entirely. An MA(1) has an ACF with a single spike at lag 1 and nothing after it, exactly as the rule of thumb said.

An ARMA(p,q) model combines both: $\phi(B)y_t = c + \theta(B)\varepsilon_t$.

ARIMA

ARIMA(p,d,q) applies ARMA to the series after $d$ differences:

The I stands for integrated: forecasts are built on the differenced scale and summed back up. A useful fact links this to the previous sections: ARIMA(0,1,1) produces the same point forecasts as SES, and ARIMA(0,2,2) matches Holt's linear method. The two families overlap; they are not rivals from different worlds.

Seasonal ARIMA

SARIMA adds seasonal AR, differencing and MA terms that operate at lag $m$, written $\text{SARIMA}(p,d,q)(P,D,Q)_m$:

Capital letters are seasonal polynomials in $B^m$. Return to the ACF table: after $d = 1$ and $D = 1$ on the log scale, we saw negative spikes at lag 1 and lag 12. A non-seasonal MA(1) explains the lag-1 spike and a seasonal MA(1) explains the lag-12 spike. That is $\text{ARIMA}(0,1,1)(0,1,1)_{12}$ on $\log y_t$, so famous it is called the airline model:

When we run AutoARIMA on the log training data below, it selects exactly this model, with $\theta_1 \approx -0.33$ and $\Theta_1 \approx -0.58$. On the raw scale it picks a different model, and the log version forecasts noticeably better. The log transform is not cosmetic: it makes the constant-variance assumption roughly true, and every interval and likelihood calculation leans on that assumption.

How AutoARIMA searches

Hyndman and Khandakar's algorithm, which statsforecast reimplements, works in three steps. First choose $D$ with a seasonal strength test and $d$ with repeated KPSS unit-root tests, because AICc cannot compare models fitted to differently differenced data. Then fit a few starting models and take the best by AICc. Finally, do a stepwise search: vary $p$, $q$, $P$, $Q$ by plus or minus one, toggle the constant, keep any improvement, and stop when no neighbour improves. This visits a few dozen models instead of hundreds.

That automates the classical Box-Jenkins loop:

flowchart TD
  A["Plot the series"] --> B{"Variance grows with level?"}
  B -->|"yes"| C["Log or Box-Cox transform"]
  B -->|"no"| D["Check stationarity"]
  C --> D
  D --> E["Difference: d and seasonal D"]
  E --> F["Identify p, q, P, Q from ACF and PACF"]
  F --> G["Estimate by maximum likelihood"]
  G --> H{"Residuals look like white noise?"}
  H -->|"no"| F
  H -->|"yes"| I["Forecast and backtest"]

The Theta method and Croston for intermittent demand

Theta

The Theta method was introduced by Assimakopoulos and Nikolopoulos and became well known after performing strongly in the M3 forecasting competition. The original description decomposes a series into "theta lines" that exaggerate or flatten its local curvature: the $\theta = 0$ line is a straight regression line through the data (pure long-run trend), and the $\theta = 2$ line doubles the short-term wiggles. The forecast extrapolates the $\theta = 0$ line, forecasts the $\theta = 2$ line with SES, and averages the two.

Hyndman and Billah later showed that, for the standard version, this is equivalent to SES with drift, where the drift equals half the slope of a linear trend fitted to the data:

where $\ell_T$ is the SES level and $b_0$ the regression slope. In plain words: a smoothed level plus a conservative, half-strength trend. Seasonality is handled by deseasonalising first (classically with multiplicative decomposition) and reseasonalising the forecast. AutoTheta chooses among standard, optimised and dynamic variants by in-sample fit; on our training data it chose the optimised Theta model with multiplicative decomposition.

Croston and intermittent demand

Spare parts, slow-moving retail items and many inventory series are mostly zeros with occasional non-zero values. SES on such a series is pulled down after every zero and spikes after every sale, so its forecast immediately after a sale is biased high. Croston's method splits the series into two parts:

• $z$: the size of each non-zero demand,
• $p$: the number of periods between non-zero demands.

Both are smoothed with SES, updated only when a demand occurs, and the forecast per period is their ratio:

Example: typical order size $5$, typical gap $4$ weeks, forecast $1.25$ units per week. That is a rate for planning stock, not a prediction that next week will have $1.25$ sales. The Syntetos-Boylan approximation (SBA) multiplies by $(1 - \alpha/2)$ to correct a known upward bias, and TSB smooths the probability of a demand occurring each period instead of the interval, so it can decay toward zero for an item that has become obsolete. All three are available in statsforecast as CrostonClassic, CrostonSBA and TSB.

Prediction intervals from statistical models

A prediction interval states a range that should contain the future value with a given probability, if the model's assumptions hold. With Gaussian errors and $h$-step forecast standard deviation $\sigma_h$:

The models differ in how $\sigma_h$ grows with the horizon. With $\sigma$ the standard deviation of one-step residuals:

Read the SES row slowly. Each future step the level absorbs a fraction $\alpha$ of a new shock, so uncertainty accumulates at rate $\alpha^2$. With $\alpha$ near 1 it grows almost like the naive random walk; with small $\alpha$ the intervals stay narrow because the level barely moves. For ARIMA models $\sigma_h$ comes from writing the model as an infinite MA and summing the squared weights. For multiplicative ETS models there is often no closed form, and software uses approximations or simulation.

Two practical consequences. First, if you fit on the log scale, transform interval endpoints back with $\exp$. Quantiles survive monotone transforms, so this gives valid intervals, and $\exp$ of the point forecast is the median on the original scale, not the mean. The mean is larger, approximately $\exp(\mu + \sigma^2/2)$ for log-normal errors. Second, these intervals ignore parameter and model-selection uncertainty, so they tend to be too narrow in practice. Always check empirical coverage in a backtest, as we do below.

Implement it: every model in statsforecast

Install into a fresh virtual environment. The results below were produced with Python 3.14, statsforecast 2.1.1, utilsforecast 0.2.17, NumPy 2.5.4 and pandas 2.3.3 on CPU.

python3 -m venv .venv && source .venv/bin/activate
pip install statsforecast utilsforecast

Load and split

The data comes in Nixtla's long format: one row per series and timestamp, with columns unique_id, ds and y. Many series can live in one frame; here there is just one. The timestamps are month ends, so the frequency is "ME".

from functools import partial

import numpy as np
from statsforecast import StatsForecast
from statsforecast.models import (
    Naive, SeasonalNaive, HistoricAverage, RandomWalkWithDrift,
    AutoETS, HoltWinters, AutoARIMA, AutoTheta,
)
from statsforecast.utils import AirPassengersDF
from utilsforecast.evaluation import evaluate
from utilsforecast.losses import mae, rmse, mase

df = AirPassengersDF.copy()
h, m = 12, 12
train = df.iloc[:-h]          # Jan 1949 .. Dec 1959, 132 rows
test = df.iloc[-h:]           # Jan 1960 .. Dec 1960, 12 rows
print(df.ds.min(), df.ds.max(), len(train), len(test))

Fit and forecast with intervals

StatsForecast takes a list of model objects and a frequency. forecast fits every model on every series and returns point forecasts plus interval columns named like AutoETS-lo-95 and AutoETS-hi-95. Use alias to run two configurations of the same class side by side.

models = [
    Naive(),
    SeasonalNaive(season_length=m),
    HistoricAverage(),
    RandomWalkWithDrift(),
    HoltWinters(season_length=m, error_type="A", alias="HW_add"),   # ETS(A,A,A)
    HoltWinters(season_length=m, error_type="M", alias="HW_mult"),  # ETS(M,A,M)
    AutoETS(season_length=m),    # searches the ETS taxonomy by AICc
    AutoARIMA(season_length=m),  # stepwise SARIMA search by AICc
    AutoTheta(season_length=m),
]
sf = StatsForecast(models=models, freq="ME", n_jobs=1)
fcst = sf.forecast(df=train, h=h, level=[80, 95])

The log-scale SARIMA

This statsforecast version raises NotImplementedError if you pass a Box-Cox blambda to AutoARIMA, so do the transform by hand: fit on $\log y$, forecast, and exponentiate every column. The exponentiated point forecast is a median.

log_train = train.assign(y=np.log(train["y"]))
sf_log = StatsForecast(
    models=[AutoARIMA(season_length=m, alias="AutoARIMA_log")], freq="ME"
)
fl = sf_log.forecast(df=log_train, h=h, level=[80, 95])
for col in fl.columns:
    if col.startswith("AutoARIMA_log"):
        fl[col] = np.exp(fl[col])
fcst = fcst.merge(fl, on=["unique_id", "ds"])

Score on the test year

utilsforecast.evaluate computes metrics per series and model. MASE needs the training frame for its seasonal-naive scale.

res = test.merge(fcst, on=["unique_id", "ds"])
names = [c for c in fcst.columns
         if c not in ("unique_id", "ds") and "-lo-" not in c and "-hi-" not in c]
scores = evaluate(
    res,
    metrics=[mae, rmse, partial(mase, seasonality=m)],
    models=names,
    train_df=train,
)
table = (scores.drop(columns="unique_id").set_index("metric").T
         .round(3).sort_values("mase"))
print(table)

See what the automatic models chose

forecast does not keep the fitted models, so call fit to inspect them. In this version each fitted model stores a dictionary in model_. For ARIMA, arma is ordered (p, q, P, Q, m, d, D).

sf_fit = StatsForecast(
    models=[AutoETS(season_length=m), AutoARIMA(season_length=m),
            AutoTheta(season_length=m)],
    freq="ME",
).fit(train)
for fitted in sf_fit.fitted_[0]:
    info = fitted.model_
    print(fitted.alias, info.get("method"), info.get("arma"), info.get("modeltype"))

sf_log_fit = StatsForecast(models=[AutoARIMA(season_length=m)], freq="ME").fit(log_train)
print("log scale:", sf_log_fit.fitted_[0][0].model_["arma"],
      sf_log_fit.fitted_[0][0].model_["coef"])

On this run: AutoETS chose ETS(M,N,M), multiplicative error and season with no trend; raw-scale AutoARIMA chose $\text{ARIMA}(1,1,0)(0,1,0)_{12}$; log-scale AutoARIMA chose the airline model $\text{ARIMA}(0,1,1)(0,1,1)_{12}$ with MA coefficients $-0.327$ and seasonal MA $-0.578$; AutoTheta chose the optimised Theta model.

Check interval coverage

for name in ["SeasonalNaive", "AutoETS", "AutoARIMA", "AutoARIMA_log", "HW_mult"]:
    for level in (80, 95):
        lo, hi = res[f"{name}-lo-{level}"], res[f"{name}-hi-{level}"]
        covered = ((res.y >= lo) & (res.y <= hi)).mean()
        print(f"{name:14s} {level}%  coverage={covered:.2f}  width={(hi - lo).mean():.0f}")

Rolling-origin cross-validation

One test year is one draw. cross_validation refits at several forecast origins inside the training data. With n_windows=3 and step_size=12, the cutoffs are December 1956, 1957 and 1958, each followed by a 12-month forecast. Only training data is used, so the 1960 test year stays untouched.

cv_models = [SeasonalNaive(season_length=m), AutoETS(season_length=m),
             AutoARIMA(season_length=m), AutoTheta(season_length=m)]
cv = StatsForecast(models=cv_models, freq="ME").cross_validation(
    df=train, h=h, n_windows=3, step_size=h
)
cv_names = [c for c in cv.columns if c not in ("unique_id", "ds", "cutoff", "y")]
cv_scores = evaluate(cv, metrics=[mae, rmse], models=cv_names)
print(cv_scores.round(2))                                   # one row per cutoff
print(cv_scores.drop(columns=["unique_id", "cutoff"])
      .groupby("metric").mean().round(2))                   # averaged over windows

flowchart LR
  D["Training history"] --> C1["Cutoff 1: fit on data up to cutoff"]
  C1 --> F1["Forecast next h steps"]
  F1 --> S1["Score against actuals"]
  S1 --> N{"More cutoffs?"}
  N -->|"yes, move cutoff forward by step_size"| C1
  N -->|"no"| A["Average errors by model and horizon"]
  A --> P["Pick model, then refit on all training data"]
  P --> T["Score once on the untouched test year"]

Implement it: the recursions in NumPy

Libraries hide the most instructive part: the loop. Here are SES and damped Holt in a few lines each. Every value in fitted is computed before the observation it predicts is read, which is the information boundary from the intuition guide made concrete.

import numpy as np
from statsforecast.utils import AirPassengersDF


def ses(y, alpha, level0):
    """Simple exponential smoothing: one-step fitted values and final level."""
    level = level0
    fitted = np.empty(len(y))
    for t, obs in enumerate(y):
        fitted[t] = level                        # forecast of y_t made at t-1
        level = alpha * obs + (1 - alpha) * level
    return fitted, level


def holt(y, alpha, beta, level0, trend0, phi=1.0):
    """Holt's linear trend method, damped when phi < 1."""
    level, trend = level0, trend0
    fitted = np.empty(len(y))
    for t, obs in enumerate(y):
        fitted[t] = level + phi * trend
        prev_level = level
        level = alpha * obs + (1 - alpha) * (prev_level + phi * trend)
        trend = beta * (level - prev_level) + (1 - beta) * phi * trend
    return fitted, level, trend


def holt_forecast(level, trend, h, phi=1.0):
    steps = np.cumsum(phi ** np.arange(1, h + 1))   # phi + phi^2 + ... + phi^k
    return level + steps * trend


y = AirPassengersDF["y"].to_numpy()[:-12]

# SES: choose alpha by the sum of squared one-step errors
alphas = np.linspace(0.01, 0.99, 99)
sse = lambda a: np.sum((y - ses(y, a, y[0])[0]) ** 2)
a_best = min(alphas, key=sse)
_, last_level = ses(y, a_best, y[0])
print(f"SES alpha={a_best:.2f}  flat forecast={last_level:.1f}")

# Damped Holt with phi fixed at 0.98
grid = [(a, b) for a in np.linspace(0.05, 0.95, 19) for b in np.linspace(0.01, 0.5, 50)]
def holt_sse(p, phi=0.98):
    fitted, *_ = holt(y, p[0], p[1], y[0], y[1] - y[0], phi)
    return np.sum((y - fitted) ** 2)
a_h, b_h = min(grid, key=holt_sse)
_, lv, tr = holt(y, a_h, b_h, y[0], y[1] - y[0], phi=0.98)
print(f"Holt alpha={a_h:.2f} beta={b_h:.2f} level={lv:.1f} trend={tr:.2f}")
print(np.round(holt_forecast(lv, tr, 12, phi=0.98), 1))

# Geometric weights for alpha = 0.3
print(np.round(0.3 * 0.7 ** np.arange(6), 4))

Output from this run: SES chose $\alpha = 0.99$ with a flat forecast of $404.6$; damped Holt chose $\alpha = 0.95$, $\beta = 0.01$ and forecasts that creep from $403.6$ to $409.7$. A cross-check with statsforecast's SimpleExponentialSmoothing(alpha=0.99) gives the same $404.6$.

That $\alpha = 0.99$ is the lesson. SES has no seasonal state, so the only way it can reduce one-step errors on a strongly seasonal series is to chase the last observation: it collapses into the naive method. A smoothing parameter near 1 is often a symptom of missing structure rather than a sign that the data is genuinely that volatile.

So add the seasonal state, and use the log transform to make the seasonality additive:

import itertools

import numpy as np
from statsforecast.utils import AirPassengersDF


def holt_winters_additive(y, alpha, beta, gamma, m, phi=1.0):
    """Additive Holt-Winters: one-step fitted values and final states."""
    level = y[:m].mean()
    trend = (y[m:2 * m].mean() - y[:m].mean()) / m
    season = list(y[:m] - level)               # season[t % m] holds s_(t-m)
    fitted = np.empty(len(y))
    for t, obs in enumerate(y):
        s = season[t % m]
        fitted[t] = level + phi * trend + s
        prev_level, prev_trend = level, trend
        level = alpha * (obs - s) + (1 - alpha) * (prev_level + phi * trend)
        trend = beta * (level - prev_level) + (1 - beta) * phi * trend
        season[t % m] = gamma * (obs - prev_level - phi * prev_trend) + (1 - gamma) * s
    return fitted, level, trend, season


def hw_forecast(level, trend, season, n_obs, h, m, phi=1.0):
    steps = np.cumsum(phi ** np.arange(1, h + 1))
    seas = np.array([season[(n_obs + k) % m] for k in range(h)])
    return level + steps * trend + seas


y_all = AirPassengersDF["y"].to_numpy().astype(float)
train, test = y_all[:-12], y_all[-12:]
z = np.log(train)
m = 12

grid = np.linspace(0.05, 0.95, 19)
def sse(p):
    fitted, *_ = holt_winters_additive(z, *p, m)
    return np.sum((z[m:] - fitted[m:]) ** 2)    # skip the initialisation year

candidates = [p for p in itertools.product(grid, [0.01, 0.02, 0.05, 0.1, 0.2], grid)
              if p[2] <= 1 - p[0]]              # usual constraint gamma <= 1 - alpha
best = min(candidates, key=sse)
_, lv, tr, se = holt_winters_additive(z, *best, m)
pred = np.exp(hw_forecast(lv, tr, se, len(z), 12, m))

mae_val = np.mean(np.abs(test - pred))
scale = np.mean(np.abs(train[m:] - train[:-m]))   # in-sample seasonal-naive MAE
print("alpha, beta, gamma =", np.round(best, 2))
print("MAE", round(mae_val, 2),
      "RMSE", round(np.sqrt(np.mean((test - pred) ** 2)), 2),
      "MASE", round(mae_val / scale, 3))

On this run the grid search picks $\alpha = 0.35$, $\beta = 0.01$, $\gamma = 0.55$, and the test-year scores are MAE $17.44$, RMSE $22.35$, MASE $0.573$. Forty lines of NumPy with a crude grid beat every untransformed model in statsforecast, which says more about the value of the log transform than about the code.

For completeness, here is Croston on a synthetic spare-parts series where roughly one week in four has a sale:

import numpy as np
import pandas as pd
from statsforecast import StatsForecast
from statsforecast.models import CrostonClassic, CrostonSBA, TSB, Naive

rng = np.random.default_rng(0)
n = 104
demand = np.where(rng.random(n) < 0.25, rng.poisson(4, n) + 1, 0)
spare = pd.DataFrame({"unique_id": "part_17",
                      "ds": pd.date_range("2023-01-01", periods=n, freq="W"),
                      "y": demand.astype(float)})
sf = StatsForecast(models=[CrostonClassic(), CrostonSBA(),
                           TSB(alpha_d=0.2, alpha_p=0.2), Naive()], freq="W")
print(sf.forecast(df=spare, h=4).round(2))

With 78% zero weeks in this sample, the run gives flat rates of $1.07$ (Croston), $1.01$ (SBA) and $0.98$ (TSB) units per week, while naive forecasts $0$ because the last week happened to be empty: a forecast that would never reorder.

Reading the results

Test-year scores (January to December 1960) from the statsforecast run above and the NumPy Holt-Winters. Numbers come from one run with the library versions listed earlier; they can shift slightly across versions because optimisers and search defaults change.

What to take from this:

• Seasonal naive is the bar. Every model that captures both trend and season clears it comfortably. Naive, drift and mean each miss one of the two patterns and lose badly.
• The log transform helps the most. The same AutoARIMA search improves MAE from 18.5 to 13.2 when the variance is stabilised first. Multiplicative Holt-Winters beats additive Holt-Winters for the same reason.
• AICc chose a model without a trend. AutoETS picked ETS(M,N,M), which had the lowest AICc in-sample but cannot extrapolate growth, so it underforecasts 1960. It still beats the seasonal naive forecast for the test year, but its MASE above 1 and its large gap to the trended models show the cost. In-sample criteria select models; backtests validate them.
• Interval coverage is uneven. At the nominal 95% level, AutoARIMA, log AutoARIMA and seasonal naive each covered 11 of 12 test months; AutoETS and multiplicative Holt-Winters covered 10 of 12. At 80%, seasonal naive covered only 3 of 12, because its intervals are centred on last year's level and 1960 was higher. Twelve points is far too few to judge calibration; this is a sanity check, not a verdict.

The rolling-origin backtest on the training data tells a more cautious story than the single test year. Averaged over the three 12-month windows (cutoffs December 1956, 1957 and 1958), MAE was 33.36 for seasonal naive, 31.92 for AutoETS, 25.67 for AutoARIMA and 25.08 for AutoTheta. In the window after the 1957 cutoff, seasonal naive had the lowest MAE of the four (12.58), because 1958 grew unusually little. A ranking from one window can flip in the next, which is why you look at several origins before you trust a choice.

Pitfalls, and how to choose a model

Pitfalls

• Comparing AICc across transformations or differencing orders. A model on $\log y$ and a model on $y$ have likelihoods on different scales. Compare them with out-of-sample errors on the original scale.
• Treating the exponentiated forecast as a mean. It is a median. If you are forecasting a total that will be summed or budgeted, apply a bias adjustment.
• Forgetting the seasonal period. Every seasonal model in statsforecast defaults to season_length=1, which silently means no seasonality. Pass the period explicitly.
• Wrong frequency string. AirPassengers stamps are month ends. Using a month-start frequency produces forecast dates such as 1960-01-01 that do not match the test frame's 1960-01-31, so a merge on ds finds no rows to score.
• Fitting on the full series before evaluation. Any parameter, transformation choice or model order picked after looking at the test year is tuned on it. Make choices with cross_validation on training data.
• Over-reading a smoothing parameter near 0 or 1. As the SES run showed, $\alpha$ near 1 often means the model is missing a component.
• Long-horizon straight lines. An undamped trend extrapolated far ahead tends to overshoot. Damped trends usually behave better beyond one season.
• Croston on smooth demand, or SES on intermittent demand. Check the share of zeros before choosing.

Choosing

flowchart TD
  S["New series"] --> Z{"Mostly zeros?"}
  Z -->|"yes"| CR["Croston, SBA or TSB"]
  Z -->|"no"| V{"Seasonal swings grow with level?"}
  V -->|"yes"| LG["Log transform, or use multiplicative seasonality"]
  V -->|"no"| SE{"Clear seasonality?"}
  LG --> SE
  SE -->|"yes"| SM["Candidates: SeasonalNaive, AutoETS,<br/>AutoARIMA with season_length, AutoTheta"]
  SE -->|"no"| TR{"Trend?"}
  TR -->|"yes"| TM["Candidates: Drift, damped Holt,<br/>AutoARIMA, Theta"]
  TR -->|"no"| LV["Candidates: Naive, SES, Mean"]
  SM --> BT["Rolling-origin backtest, keep what beats the baseline"]
  TM --> BT
  LV --> BT
  CR --> BT

Statistical models are the right first choice when you have one or a few series, short histories, a need for intervals and explanations, or thousands of series where each needs a cheap, robust fit. ETS is natural when you think in terms of level, trend and season; ARIMA when autocorrelation structure is the main signal; Theta as a robust middle ground; and an ensemble that averages ETS, ARIMA and Theta is often hard to beat. Reach for neural models or transformer models when you have many related series to learn across, rich covariates, or nonlinear patterns that a few states cannot express. Consider pretrained foundation models for zero-shot forecasts on new series. In every case, keep a fitted statistical model in the comparison: if a larger model cannot beat it in a backtest, the larger model has not earned its place.

Practice exercises

1. Derive the damped limit. Show that $\sum_{j=1}^{h}\phi^j = \phi(1-\phi^h)/(1-\phi)$ and compute the long-run forecast increment for $\phi = 0.95$ and $b_T = 2$.
2. SES equals ARIMA(0,1,1). Start from $y_t = y_{t-1} + \varepsilon_t + \theta\varepsilon_{t-1}$ and show that its one-step forecast satisfies the SES recursion with $\alpha = 1 + \theta$.
3. Damped trend in practice. Rerun AutoETS with damped=True and compare it with the undamped search. Does damping fix the missing-trend choice?
4. Bias adjustment. Estimate the residual variance of the log airline model and multiply its forecasts by $\exp(\sigma^2/2)$. Does the MAE improve, worsen, or stay the same? Why might MAE prefer the median?
5. More origins. Increase n_windows to 6 with step_size=6 and report the error by horizon, from 1 to 12 months. Which model degrades fastest with horizon?
6. Ensemble. Average the AutoETS, log AutoARIMA and AutoTheta forecasts and score the average. Compare it with each member across the cross-validation windows, not just on 1960.
7. Intermittent demand. Increase the zero share in the Croston example to 90% and compare Croston, SBA, TSB and SES with a rolling backtest using MAE and the mean error (bias).

Further reading

• Hyndman and Athanasopoulos, Forecasting: Principles and Practice, 3rd edition, especially the chapters on exponential smoothing and ARIMA models.
• Hyndman, Koehler, Snyder and Grose, A state space framework for automatic forecasting using exponential smoothing methods, International Journal of Forecasting, 2002.
• Hyndman and Khandakar, Automatic time series forecasting: the forecast package for R, Journal of Statistical Software, 2008: the AutoARIMA search.
• Hyndman and Billah, Unmasking the Theta method, International Journal of Forecasting, 2003.
• Nixtla, statsforecast documentation.