ARIMA and Exponential Smoothing Explained With Real Numbers

Every forecast is a bet that the past contains a pattern that will keep going. But which pattern? Do you trust the last value, the recent average, a trend, or a rhythm? Today we study the two classic answers: ARIMA and exponential smoothing. They look very different, and by the end you will see how closely related they are.

• ARIMA
• Exponential smoothing

In lesson 45 we tuned hyperparameters with Bayesian optimization. Here we return to time series from lesson 25, where we met trend, seasonality, stationarity, and the ACF and PACF plots. We will use those ideas only as a quick reference, then build the actual models, with numbers you can check by hand.

• Stationarity and differencing
• ARIMA family
• Exponential smoothing
• Python evaluation

Quick recall from lesson 25. A series is weakly stationary when its mean and variance stay constant and its correlation between two points depends only on the gap between them, called the lag. ARMA models need this property, because they learn one fixed set of coefficients. If the mean drifts upward, one set of coefficients cannot describe both the early and late data.

• Constant mean
• Constant variance
• Correlation depends on lag only

The fix is differencing: replace each value with its change from the previous value. Take the series one hundred, one hundred three, one hundred seven, one hundred eight, one hundred twelve. The differences are three, four, one, four. The upward climb has become a flat series wiggling around about three point places. One round of differencing removes a linear-looking trend.

• Replace value with its change
• 100, 103, 107, 108, 112
• Differences: 3, 4, 1, 4

A random walk is the classic case. Each value equals the previous value plus a fresh random shock. It has no fixed level to return to, so shocks never fade. We say it has a unit root. Differencing a random walk gives back pure shocks, which are stationary. The number of differences needed is called d in ARIMA.

• Shocks never fade
• Unit root
• Difference once: d = 1

How do we test for a unit root? The augmented Dickey-Fuller test, or ADF. It regresses the change in the series on its previous level, plus some lagged changes. The key coefficient, gamma, is zero when there is a unit root. So the null hypothesis is non-stationary. A small p-value, say below zero point zero five, lets us reject it and treat the series as stationary.

• Null: unit root exists
• Small p-value: reject, stationary
• Gamma equals zero means unit root

Two common mistakes here. First, a large p-value does not prove a unit root; it only says you lack evidence against one. Second, do not difference blindly. Over-differencing injects artificial negative correlation and inflates variance. A practical rule: difference until the ADF test passes and the series looks flat, and almost never go beyond d equals two.

• Large p is not proof
• Over-differencing adds noise
• Rarely exceed d = 2

Now the first building block: the autoregressive model. An AR of order p predicts today as a weighted sum of the last p values plus noise. For AR one: y at t equals phi times y at t minus one, plus epsilon at t. Phi is the memory strength. If phi is zero point eight, each value keeps eighty percent of yesterday's deviation from the mean.

• Weighted sum of past values
• Phi is memory strength
• Plus white noise

The ACF of an AR one with phi equal to zero point eight decays geometrically: zero point eight, zero point six four, zero point five one two. Correlation at lag k is phi to the k. The influence echoes forward, shrinking each step. For stationarity we need the absolute value of phi below one. At phi equal to one, the echo never fades, and that is the random walk.

• ACF at lag k is phi^k
• Needs |phi| less than 1
• phi = 1 is a random walk

How do we estimate AR coefficients? Yule-Walker. Multiply the model by a past value, take expectations, and the equations link autocorrelations to coefficients. For AR two: rho one equals phi one plus phi two times rho one, and rho two equals phi one times rho one plus phi two. Autocorrelations go in, coefficients come out, by solving a small linear system.

• Link autocorrelations to coefficients
• A small linear system

A hand example. Suppose the sample autocorrelations are rho one equals zero point six and rho two equals zero point two. The equations are: phi one plus zero point six phi two equals zero point six, and zero point six phi one plus phi two equals zero point two. Substitute phi one equals zero point six minus zero point six phi two into the second. You get zero point six four phi two equals negative zero point one six, so phi two is negative zero point two five, and phi one is zero point seven five.

• rho1 = 0.6, rho2 = 0.2
• phi2 = -0.25
• phi1 = 0.75

Always verify. Plug back in: zero point seven five plus zero point six times negative zero point two five gives zero point seven five minus zero point one five, which is zero point six. Correct. And zero point six times zero point seven five minus zero point two five is zero point four five minus zero point two five, which is zero point two. Correct again. This is exactly what the statsmodels function yule_walker computes, just with more lags.

• Plug back into both equations
• Both check out

The moving average model is different. An MA of order q predicts today from the last q shocks, not past values. For MA one: y at t equals epsilon at t plus theta times epsilon at t minus one. A shock hits now, echoes once at the next step, and then vanishes completely. The memory is finite, unlike AR.

• Weighted sum of past shocks
• Memory lasts exactly q steps

This gives MA a signature. Its ACF cuts off after lag q. For MA one with theta equal to zero point five, rho one is theta over one plus theta squared, which is zero point five over one point two five, or zero point four. Every lag beyond one is exactly zero. Note this has nothing to do with the smoothing moving average from lesson 25. Here the average is of noise.

• rho1 = theta / (1 + theta^2)
• theta = 0.5 gives 0.4
• Zero after lag q

Combine them and you get ARMA p, q: past values and past shocks together. This is more economical than either alone, since a few parameters can mimic a long pattern. To write all of this compactly, we use the lag operator B. B applied to y at t gives y at t minus one. Then the AR part is a polynomial in B on the left, and the MA part is a polynomial in B on the right.

• B shifts back one step
• phi(B) y = theta(B) epsilon

The lag operator makes differencing easy to write. One minus B applied to y gives y at t minus y at t minus one, the first difference. So ARIMA p, d, q is simply ARMA applied to the d times differenced series. Phi of B, times one minus B to the d, times y, equals theta of B times epsilon. The I stands for integrated, the reverse of differencing.

• p: autoregressive order
• d: differences
• q: moving average order

How do we choose p and q? Read the plots. For a pure AR p, the PACF cuts off after lag p while the ACF decays. For a pure MA q, the ACF cuts off after lag q while the PACF decays. If both decay gradually, you likely need a mixed ARMA. The PACF at lag k is the correlation at that lag after removing the influence of the shorter lags in between.

• AR: PACF cuts off at p
• MA: ACF cuts off at q
• ARMA: both tail off

Real plots are messy, so treat the reading as a starting guess. A sensible workflow is to fit a few candidates near your guess and compare them with AIC, which rewards fit but penalizes extra parameters. Lower AIC is better. Remember lesson 24: only compare AIC across models fit to the same differenced data. Changing d changes the data itself.

• Fit several candidates
• Lower AIC wins
• Same d only

Estimation is by maximum likelihood, as in lesson 10. Assume the shocks are independent normal noise. Given coefficients, we can compute each one step prediction error, then the likelihood of the whole series is a product of normal densities of those errors. An optimizer searches the coefficients that maximize it. For pure AR models this is close to least squares, but MA terms need numerical optimization, since the past shocks are not observed.

• Normal shock assumption
• Maximize likelihood of errors
• MA needs numerical search

After fitting, check the residuals. A good model leaves behind white noise: no remaining correlation. The Ljung-Box test checks this. Its statistic Q is n times n plus two, times the sum over lags of squared residual autocorrelation divided by n minus k. The null says residuals are uncorrelated. We want a large p-value here, the opposite of what we wanted for ADF.

• Residuals should be white noise
• Ljung-Box: want large p-value
• Check ACF and histogram

If Ljung-Box rejects, the model missed structure. Look at where the residual ACF spikes. A spike at lag one suggests adding an MA or AR term, and a spike at lag twelve in monthly data suggests seasonality. Also remember to adjust the degrees of freedom by the number of fitted ARMA parameters. Chasing a perfect p-value by adding terms blindly, however, leads to overfitting.

• Spike at lag 1: add term
• Spike at lag 12: seasonality
• Do not overfit

Now the second family. Simple exponential smoothing forecasts with a level that is updated each step. The new level equals alpha times the newest observation, plus one minus alpha times the old level. Alpha between zero and one controls how fast the model reacts. Large alpha trusts recent data. Small alpha averages over a long history.

• New level = blend
• Alpha controls reaction speed

A numeric example with alpha equal to zero point five. Data: ten, twelve, eleven, thirteen. Start the level at ten. After twelve: half of twelve plus half of ten is eleven. After eleven: half of eleven plus half of eleven is eleven. After thirteen: half of thirteen plus half of eleven is twelve. The forecast for every future step is twelve.

• Data: 10, 12, 11, 13
• Levels: 10, 11, 11, 12
• Forecast: 12 forever

Why is it called exponential? Unroll the recursion. The level is a weighted average of all past observations with weights alpha, alpha times one minus alpha, alpha times one minus alpha squared, and so on. For alpha of one half: zero point five, zero point two five, zero point one two five. These weights shrink geometrically. Check: half times thirteen, plus a quarter times eleven, plus an eighth times twelve, plus an eighth of the starting ten, gives six point five plus two point seven five plus one point five plus one point two five, which is twelve.

• Weights: alpha(1 - alpha)^k
• 0.5, 0.25, 0.125, ...
• Sum of weights is 1

Simple smoothing has no trend, so its forecast is flat. Holt's method adds a trend component b. The level update uses the previous level plus the previous trend as its baseline, and the trend is itself smoothed from changes in the level, with its own parameter beta. Suppose level is one hundred, trend is two, new observation is one hundred five, alpha is zero point five, beta is zero point three.

• Level plus trend
• Alpha smooths level
• Beta smooths trend

Walk through it. The one step forecast is level plus trend, one hundred two. The new level is half of one hundred five plus half of one hundred two, which is one hundred three point five. The level moved by three point five, so the new trend is zero point three times three point five plus zero point seven times two, equals one point zero five plus one point four, which is two point four five. A three step forecast is one hundred three point five plus three times two point four five, or one hundred ten point eight five.

• New level: 103.5
• New trend: 2.45
• 3-step forecast: 110.85

Holt-Winters adds a third component for seasonality, with a third parameter, gamma. Seasonality can be additive, where the seasonal swing has constant size, or multiplicative, where it grows with the level. Pick additive when the peaks keep the same height, and multiplicative when they grow with the series. A related caution: a trend extrapolated forever often overshoots, so damped trends are a popular safer choice.

• Level, trend, seasonal
• Additive vs multiplicative
• Damped trend is safer

Here is the connection promised at the start. Simple exponential smoothing is exactly the same as an ARIMA zero, one, one model, with theta equal to alpha minus one. Holt's linear method corresponds to ARIMA zero, two, two. So the families overlap. Exponential smoothing is a restricted, very robust corner of the ARIMA world, and later state space formulations put both under one framework.

• SES equals ARIMA(0,1,1)
• Holt equals ARIMA(0,2,2)

For seasonality in ARIMA, we add seasonal terms working at multiples of the season length m. SARIMA is written p, d, q times capital P, D, Q with subscript m. Seasonal differencing subtracts the value one season ago, so for monthly data with m of twelve, y at t minus y at t minus twelve. The famous airline model, zero, one, one times zero, one, one at twelve, is a strong baseline for monthly seasonal series.

• Seasonal AR, MA, differencing
• m is season length
• Airline model baseline

In Python, statsmodels does the heavy lifting. Import ARIMA from statsmodels dot tsa dot arima dot model, create it with your series and order, for example one, one, one, and call fit. The summary shows coefficients, standard errors, AIC, and a Ljung-Box test. For smoothing, use ExponentialSmoothing or ETSModel with trend and seasonal options. Calling get_forecast gives predictions and confidence intervals.

• ARIMA(y, order=(p,d,q))
• fit, summary, get_forecast
• ExponentialSmoothing for ETS

Evaluation must respect time. A random split leaks the future into training, a mistake we warned about in lesson 35. Use rolling origin evaluation. Fit on the first chunk, forecast the next horizon, record the error, then slide the origin forward, refit, and repeat. Average the errors, using MAE or a scaled error, across many origins. A single train and test split can look lucky by chance.

• Fit, forecast, record error
• Move origin forward
• Average over many origins

A forecast without uncertainty is only half a forecast. ARIMA gives a prediction interval from the model's variance. Intervals widen with the horizon, because uncertainty accumulates. With d equal to one, shocks never fade, so the band keeps growing. Plot the history, the central forecast, and a shaded ninety five percent band. If actual values land outside the band far more than five percent of the time, your intervals are too narrow.

• Intervals widen with horizon
• Shaded 95 percent band
• Check coverage

How do these compare to machine learning with lag features? You can turn forecasting into regression: build columns for the last k values, calendar features, and rolling averages, then train gradient boosting from lesson 42. This handles many series, extra predictors, and nonlinearity. But trees cannot extrapolate a trend beyond the training range, so they often need differenced targets.

• Lags become columns
• Handles extra predictors
• Trees cannot extrapolate trends

So which should you pick? For one short series with clear structure, ARIMA or exponential smoothing is quick, interpretable, and gives principled intervals. For thousands of related series or many external drivers, lag feature models often win. Whatever you choose, always include a naive baseline, such as the last value or last season's value, and use rolling origin evaluation. If a fancy model cannot beat the baseline, it is not worth deploying.

• Small and clean: ARIMA or ETS
• Large and rich: ML with lags
• Always beat the baseline

Recap. Difference to reach stationarity and confirm with the ADF test. AR uses past values, MA uses past shocks, and ARIMA combines them with differencing, written compactly using the lag operator. Read the PACF and ACF, fit by maximum likelihood, compare by AIC, and check residuals with Ljung-Box. Exponential smoothing is a weighted average with geometric weights, extended by Holt and Holt-Winters.

• Stationarity and ADF
• ARIMA and diagnostics
• Smoothing as geometric weights
• Rolling-origin evaluation

That is lesson 46. Next, in lesson 47, we move into deep learning training: weight initialization, vanishing gradients, and optimizers like momentum and Adam.

• Initialization and optimizers
• Subscribe for the full course