ARIMA and SARIMA: differencing, forecasting and prediction intervals

Before this lesson: lesson 02 (stationarity, differencing, unit roots), lesson 04 (simple exponential smoothing) and lesson 05 (AR, MA, ARMA, $\psi$-weights, AIC). Lesson 03 gives the benchmarks and backtests we compare against. Time: about 90 minutes, plus the notebook.

---

Why this matters #

• The standard statistical forecaster. ARIMA, with automatic order selection (auto.arima in R, AutoARIMA in statsforecast, pmdarima), is a benchmark in almost every forecasting paper and in competitions such as M4 and M5. A new ML model that cannot beat it, or a seasonal naive forecast, is not useful.
• Honest uncertainty. The ARIMA fan shows how uncertainty grows with the horizon, and it is exact under the model's assumptions. It is the reference for the probabilistic forecasts of lesson 10.
• Links to the rest of the track. Simple exponential smoothing (lesson 04) is an ARIMA(0,1,1); every ARIMA is a state-space model whose likelihood the Kalman filter computes (lesson 07); residuals of a seasonal ARIMA are a common anomaly score (lesson 11).
• Science. The CO₂ record of Mauna Loa, river flows, electricity load and epidemic counts are all modelled with seasonal ARIMA or regression with ARIMA errors.

1. Predict first #

2. The picture #

Differencing removes a unit root. A random walk $y_t = y_{t-1} + \varepsilon_t$ is not stationary, but its difference $\nabla y_t = y_t - y_{t-1} = \varepsilon_t$ is white noise (lesson 02). A linear trend disappears after one difference; a stable seasonal pattern with period $m$ disappears after a seasonal difference $\nabla_m y_t = y_t - y_{t-m}$. ARIMA models the differences as ARMA and then "integrates" (sums) them back.

The forecast fan. Forecasts start at the forecast origin $n$ (yellow line) and go into the future. Around the point forecast (orange line) we draw the 80 % and 95 % prediction intervals (orange bands). The figure lets you compare model families.

Interactive figure: loading…

Why do the bands grow? Write the future value as a sum of shocks: $y_{n+h} = \text{(known part)} + \psi_0\varepsilon_{n+h} + \psi_1\varepsilon_{n+h-1} + \dots + \psi_{h-1}\varepsilon_{n+1}$. The $h$ future shocks are unknown, so they make up the forecast error. Each one adds $\psi_j^2\sigma^2$ to the variance. If the $\psi_j$ shrink quickly (a stationary model), the sum converges and the bands reach a fixed width. If the $\psi_j$ do not shrink (a differenced model, where every shock acts forever), the sum grows without limit.

Interactive figure: loading…

3. Definitions and main results #

ARIMA($p, d, q$) #

Let $\nabla = 1 - B$ be the difference operator, so $\nabla y_t = y_t - y_{t-1}$ and $\nabla^2 y_t = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2})$. An ARIMA($p, d, q$) model says that the $d$-times differenced series $w_t = \nabla^d y_t$ is a stationary, invertible ARMA($p, q$):

with $\phi(z) = 1 - \phi_1 z - \dots - \phi_p z^p$ and $\theta(z) = 1 + \theta_1 z + \dots + \theta_q z^q$ (plus-sign convention, as in lesson 05). The AR part of the whole model, $\phi(z)(1 - z)^d$, has $d$ roots exactly on the unit circle at $z = 1$: the unit roots. In practice $d \in \{0, 1, 2\}$.

Special cases.

Model	Equation	Name
ARIMA(0,0,0)	$y_t = c + \varepsilon_t$	white noise with mean $c$
ARIMA(0,1,0)	$y_t = y_{t-1} + \varepsilon_t$	random walk (the naive forecast of lesson 03)
ARIMA(0,1,0) with $c$	$y_t = c + y_{t-1} + \varepsilon_t$	random walk with drift (the drift method)
ARIMA($p$,0,0)	$\phi(B)y_t = c + \varepsilon_t$	AR($p$)
ARIMA(0,1,1)	$y_t = y_{t-1} + \varepsilon_t + \theta\varepsilon_{t-1}$	simple exponential smoothing

Theorem (SES is ARIMA(0,1,1)). Simple exponential smoothing with parameter $\alpha$ gives the same one-step forecasts as ARIMA(0,1,1) with $\theta = \alpha - 1$, when both start from the same first forecast.

Proof.

1. SES in error-correction form (lesson 04): $\hat y_{t+1\mid t} = \ell_t = \ell_{t-1} + \alpha e_t$, where $e_t = y_t - \hat y_{t\mid t-1}$ and $\hat y_{t\mid t-1} = \ell_{t-1}$.
2. Replace $\ell_{t-1}$ by $y_t - e_t$: $\hat y_{t+1\mid t} = y_t - e_t + \alpha e_t = y_t + (\alpha - 1)e_t$.
3. For ARIMA(0,1,1), write the model at time $t + 1$: $y_{t+1} = y_t + \varepsilon_{t+1} + \theta\varepsilon_t$. The future shock $\varepsilon_{t+1}$ is replaced by 0 and the past shock $\varepsilon_t$ by the residual $e_t$: $\hat y_{t+1\mid t} = y_t + \theta e_t$.
4. Steps 2 and 3 are the same recursion when $\theta = \alpha - 1$. With the same start, all forecasts and residuals agree. $\blacksquare$

So $\alpha \in (0, 1)$ corresponds to $\theta \in (-1, 0)$, an invertible MA part. SES assumes that the data are a random walk plus noise; this explains its flat forecast with growing intervals.

Forecasting step by step #

Let $\hat y_{n+h\mid n}$ be the forecast made at origin $n$ for horizon $h$. The rules (FPP3 §9.8):

1. Write the model with $y_t$ alone on the left, everything else on the right.
2. Replace $t$ by $n + h$.
3. On the right, replace future observations ($y_{n+j}$ with $j \ge 1$) by their forecasts, future shocks ($\varepsilon_{n+j}$, $j \ge 1$) by 0, and past shocks ($\varepsilon_t$, $t \le n$) by the residuals $e_t$.

Why 0 for future shocks? The best forecast in mean squared error is the conditional mean $\mathrm E[y_{n+h} \mid y_1, \dots, y_n]$, and a future shock has mean 0 and is independent of the data. Past shocks are not random any more once we have seen the data: for an invertible model they are functions of the observations, and the residuals estimate them (lesson 05).

Example: AR(1) with mean $\mu$. $y_t - \mu = \phi(y_{t-1} - \mu) + \varepsilon_t$. Then $\hat y_{n+1\mid n} - \mu = \phi(y_n - \mu)$, and each further step multiplies by $\phi$ again:

The forecast decays geometrically from the last value to the mean. For any stationary ARMA the forecast tends to the mean $\mu$ as $h \to \infty$.

Interactive figure: loading…

Forecast variance and prediction intervals #

Theorem (forecast error of an ARIMA model). Let the model have the $\psi$-weights $\psi_0 = 1, \psi_1, \psi_2, \dots$, defined by $\psi(z) = \theta(z)/\bigl(\phi(z)(1 - z)^d\bigr)$. Assume the parameters are known and the shocks are white noise with variance $\sigma^2$. Then the $h$-step forecast error is

Proof idea.

1. For a stationary ARMA, the MA($\infty$) form gives $y_{n+h} = \mu + \sum_{j\ge0}\psi_j\varepsilon_{n+h-j}$.
2. The terms with $j \ge h$ contain shocks up to time $n$: they are known from the data, and they form the forecast. The terms with $j < h$ contain the $h$ future shocks: they have mean 0 and are independent of the data, so they form the error.
3. The shocks are uncorrelated, so the variances add: $\sigma^2\sum_{j<h}\psi_j^2$.
4. For $d \ge 1$ there is no stationary MA($\infty$) form, but the argument still works for a finite horizon. Unroll the model recursion from the origin $n$: $y_{n+h}$ equals a function of the data up to $n$ plus $\sum_{j<h}\psi_j\varepsilon_{n+h-j}$, where the $\psi_j$ are the first $h$ coefficients of the power series $\theta(z)/\bigl(\phi(z)(1 - z)^d\bigr)$. $\blacksquare$

The $\psi$-weights come from the recursion of lesson 05, applied to the full AR polynomial $\phi^*(z) = \phi(z)(1 - z)^d = 1 - \phi^*_1 z - \dots$: $\psi_0 = 1$, $\psi_j = \theta_j + \sum_i\phi^*_i\psi_{j-i}$.

Model	$\psi_j$ ($j \ge 1$)	$\operatorname{Var}(e_{n+h})/\sigma^2$	as $h \to \infty$
AR(1)	$\phi^j$	$(1 - \phi^{2h})/(1 - \phi^2)$	$\to 1/(1 - \phi^2)$, the variance of $y$
MA(1)	$\theta$ for $j = 1$, then 0	$1$ for $h = 1$, then $1 + \theta^2$	constant from $h = 2$
random walk	$1$	$h$	grows like $h$
ARIMA(0,1,1)	$1 + \theta$ ($= \alpha$)	$1 + (h - 1)(1 + \theta)^2$	grows like $h$

A $100(1 - a)\,\%$ prediction interval under Gaussian shocks is

At $h = 1$ every model gives $\hat y \pm 1.96\,\hat\sigma$. For a stationary model the half-width tends to $1.96$ times the standard deviation of $y$; for $d = 1$ it grows like $\sqrt h$; for $d = 2$ like $h^{3/2}$.

What $d$ and the constant do to the long run #

The constant $c$ and the number of differences $d$ fix the shape of the long-run forecast (FPP3 §9.5):

	$d = 0$	$d = 1$	$d = 2$
$c = 0$	goes to 0	goes to a constant	follows a straight line
$c \ne 0$	goes to the mean $c/\phi(1)$	follows a straight line	follows a quadratic

With $d = 1$, $c$ is the average change per step: the drift. The random walk with drift forecasts $\hat y_{n+h\mid n} = y_n + hc$, a straight line from the last value. With $d = 2$ a nonzero constant gives a quadratic trend, which is almost never sensible; software such as fable does not allow it. Higher $d$ also makes the intervals grow faster. So $d$ is a modelling decision about the long run, not only a way to pass a test.

Seasonal ARIMA #

For data with a seasonal period $m$ (12 for monthly data), the seasonal difference is $\nabla_m y_t = (1 - B^m)y_t = y_t - y_{t-m}$. A seasonal ARIMA $(p, d, q)(P, D, Q)_m$ multiplies non-seasonal and seasonal polynomials:

with $\Phi(z) = 1 - \Phi_1 z - \dots - \Phi_P z^P$ and $\Theta(z) = 1 + \Theta_1 z + \dots + \Theta_Q z^Q$. The capital letters act at multiples of $m$.

Seasonal patterns in the ACF/PACF. The seasonal part shows up at lags $m, 2m, 3m, \dots$:

• a seasonal MA(1), $(0,0,0)(0,0,1)_{12}$: one ACF spike at lag 12; the PACF decays over lags 12, 24, …;
• a seasonal AR(1), $(0,0,0)(1,0,0)_{12}$: the ACF decays over lags 12, 24, …; one PACF spike at lag 12.

The airline model $(0,1,1)(0,1,1)_{12}$ (Box and Jenkins, first used for monthly airline passenger numbers):

Expand both sides. The left side is $y_t - y_{t-1} - y_{t-12} + y_{t-13}$. The right side is $\varepsilon_t + \theta\varepsilon_{t-1} + \Theta\varepsilon_{t-12} + \theta\Theta\varepsilon_{t-13}$. So

It has MA terms at lags 1, 12 and 13, but only two free coefficients $\theta, \Theta$ (plus $\sigma^2$). The lag-13 coefficient is the product $\theta\Theta$: the multiplicative form ties it to the other two. The model says: "this month's change from last year is last month's change from last year, corrected by recent errors at both time scales".

Interactive figure: loading…

The Box–Jenkins workflow #

Box and Jenkins proposed an iterative cycle: identify, estimate, check, forecast. In the modern form of FPP3 §9.7:

1. Plot the series. Look for outliers, level shifts and changes of variance.
2. Transform if the variance grows with the level: a log or a Box–Cox transform (lesson 01).
3. Difference until the series is stationary. Choose $D$ first (seasonal strength, or a clear seasonal pattern), then $d$ (a unit-root test such as KPSS, lesson 02). Rarely use $d + D > 2$.
4. Identify candidate orders from the ACF/PACF of the differenced series (lesson 05 table, plus the seasonal patterns above).
5. Estimate candidates by maximum likelihood and compare them by AICc. Remember: only models with the same $d$ and $D$ are comparable by AICc.
6. Check the residuals: ACF inside the band, Ljung–Box test with $h - p - q - P - Q$ degrees of freedom. If they fail, go back to step 4.
7. Forecast, and back-transform if you transformed in step 2.

Back-transformation. If you model $w_t = \log y_t$ and forecast $\hat w$ with forecast variance $\sigma_h^2$, then $\exp(\hat w)$ is the forecast median of $y$, not its mean. For a Gaussian $w$ the mean is $\exp(\hat w + \sigma_h^2/2) \approx \exp(\hat w)(1 + \sigma_h^2/2)$ (FPP3 §5.6 calls this the bias adjustment). Intervals back-transform directly: $[\exp(\text{lower}), \exp(\text{upper})]$.

Automatic ARIMA. The Hyndman–Khandakar algorithm (auto.arima) automates steps 3–5: it picks $d$ by repeated KPSS tests and $D$ by a seasonal-strength measure, then searches $p, q$ (and $P, Q$) stepwise. It starts from four models, for example ARIMA(2,$d$,2) and ARIMA(0,$d$,0), and moves $p$ and $q$ by $\pm 1$, with and without the constant, as long as AICc improves. It is a strong default, but it still needs the residual check and a backtest.

Interactive figure: loading…

Regression with ARIMA errors #

Often other variables help: temperature for electricity load, prices for sales, holidays. A regression with ARIMA errors (also called dynamic regression or ARIMAX) is

The regression explains the effect of the covariates, and the ARIMA error $\eta_t$ models the autocorrelation that is left. Three points matter:

• Fit both parts together by maximum likelihood. Ordinary least squares with autocorrelated errors gives wrong standard errors and can find spurious relations between unrelated trending series (lesson 02). If $d = 1$, the software differences $y$ and all $x$ together.
• To forecast, you need the future covariates. Calendar effects are known in advance; weather must itself be forecast, and its forecast error adds to yours.
• Fourier terms for long or non-integer seasons. For daily data with a yearly season ($m = 365.25$) or several seasons, a seasonal ARIMA with $m = 365$ is impractical. Use covariates $\sin(2\pi kt/m)$ and $\cos(2\pi kt/m)$ for $k = 1, \dots, K$ (a dynamic harmonic regression), choose $K$ by AICc, and let a short ARIMA model the rest.

Changing variance: ARCH and GARCH

ARIMA assumes a constant shock variance $\sigma^2$. Financial returns violate this: calm periods and turbulent periods alternate (volatility clustering). The returns themselves are nearly uncorrelated, but their squares are strongly autocorrelated. Engle's ARCH model (1982) and Bollerslev's GARCH(1,1) (1986) model the conditional variance:

$$\varepsilon_t = \sigma_t z_t, \quad z_t \overset{iid}{\sim} \mathcal N(0, 1), \qquad \sigma_t^2 = \omega + a_1\varepsilon_{t-1}^2 + b_1\sigma_{t-1}^2 ,$$

with $\omega > 0$, $a_1, b_1 \ge 0$ and $a_1 + b_1 < 1$, so the unconditional variance $\omega/(1 - a_1 - b_1)$ is finite. A large shock raises the next variance; $b_1$ makes the effect persist. ARIMA can model the mean and GARCH the variance of the same series. The point forecast does not change, but the intervals widen after turbulent days and narrow after calm ones.

Answers to section 1. (1) The random walk. The AR(1) interval tends to $\pm 1.96/\sqrt{1 - \phi^2}$ and stops growing; the random-walk interval is $\pm 1.96\sqrt h$ and grows without limit ($\pm 19.6$ at $h = 100$). (2) A straight line: $\hat y_{n+h\mid n} = 10 + 0.5h$. (3) At lags 1, 12 and 13 (the last one with coefficient $\theta\Theta$).

4. Worked example #

(a) ARMA(1,1) forecasts and intervals by hand. Model: $y_t = 0.6\,y_{t-1} + \varepsilon_t + 0.3\,\varepsilon_{t-1}$ with $\sigma = 1$ and mean 0. At the origin, $y_n = 2.0$ and the last residual is $e_n = 0.5$.

1. $h = 1$: $y_{n+1} = 0.6y_n + \varepsilon_{n+1} + 0.3\varepsilon_n$. Replace $\varepsilon_{n+1}$ by 0 and $\varepsilon_n$ by $e_n$: $\hat y_{n+1\mid n} = 0.6 \cdot 2.0 + 0 + 0.3 \cdot 0.5 = 1.2 + 0.15 = 1.35$.
2. $h = 2$: $y_{n+2} = 0.6y_{n+1} + \varepsilon_{n+2} + 0.3\varepsilon_{n+1}$. Both shocks are in the future: $\hat y_{n+2\mid n} = 0.6 \cdot 1.35 = 0.81$.
3. $h = 3$: $\hat y_{n+3\mid n} = 0.6 \cdot 0.81 = 0.486$. From now on each forecast is $0.6$ times the previous one.
4. $\psi$-weights: $\psi_0 = 1$, $\psi_1 = \theta + \phi = 0.9$, $\psi_2 = \phi\psi_1 = 0.54$, $\psi_3 = 0.324$.
5. Variances: $h = 1$: $1$; $h = 2$: $1 + 0.81 = 1.81$; $h = 3$: $1.81 + 0.2916 = 2.1016$.
6. 95 % intervals: $1.35 \pm 1.96$; $0.81 \pm 1.96\sqrt{1.81} = 0.81 \pm 2.64$; $0.486 \pm 1.96\sqrt{2.1016} = 0.486 \pm 2.84$. As $h \to \infty$ the half-width tends to $1.96\sqrt{\sum_j\psi_j^2} = 1.96\sqrt{1 + 0.81/(1 - 0.36)} = 1.96 \cdot 1.505 = 2.95$.

(b) SES as ARIMA(0,1,1). SES with $\alpha = 0.3$ is ARIMA(0,1,1) with $\theta = -0.7$. All $\psi_j = 1 + \theta = 0.3$ for $j \ge 1$, so $\operatorname{Var}(e_{n+h}) = \sigma^2(1 + 0.09(h - 1))$. With $\sigma = 2$: the 95 % half-width is $1.96 \cdot 2 = 3.92$ at $h = 1$ and $1.96 \cdot 2\sqrt{1 + 0.09 \cdot 11} = 3.92 \cdot 1.411 = 5.53$ at $h = 12$. The point forecast is flat.

(c) The airline model on CO₂. Monthly mean CO₂ at Mauna Loa, March 1958 to December 2001 ($n = 526$; the 5 months without data are filled by linear interpolation). The variance does not grow with the level, so no log is needed. The ACF of $\nabla\nabla_{12}y$ has spikes at lag 1 ($-0.30$) and lag 12 ($-0.41$). Maximum likelihood (statsmodels SARIMAX) gives

Model	AICc	BIC
$(0,1,1)(0,1,1)_{12}$ airline	259.2	271.9
$(1,1,1)(0,1,1)_{12}$	256.9	273.8
$(0,1,2)(0,1,1)_{12}$	258.3	275.2
$(0,1,1)(1,1,1)_{12}$	260.9	277.8
$(2,1,0)(0,1,1)_{12}$	265.6	282.5
$(1,1,0)(1,1,0)_{12}$	404.0	416.7

AICc prefers $(1,1,1)(0,1,1)_{12}$ by 2.3 units; BIC prefers the airline model. A difference of about 2 is small, so we keep the simpler airline model. Its residuals pass the Ljung–Box test with $h = 24$ ($Q = 20.0$, 22 degrees of freedom, p-value 0.58). Forecasts from December 2001, with 95 % intervals:

Horizon	Month	Forecast (ppm)	$\hat\sigma_h$	95 % interval
1	Jan 2002	371.98	0.30	371.38 – 372.57
12	Dec 2002	372.54	0.72	371.12 – 373.95
24	Dec 2003	374.13	1.07	372.03 – 376.22
36	Dec 2004	375.72	1.39	372.99 – 378.44

The forecast continues the trend and the season. The interval widens like that of a model with two differences.

5. In code #

NumPy checks the worked example (a) and the $\psi$-weight variance; the airline fit uses statsmodels.

import numpy as np

def psi_weights(phi_star, theta, m):
    """psi-weights for the full AR polynomial phi*(z) = 1 - sum phi*_i z^i (plus-sign MA)."""
    psi = np.zeros(m); psi[0] = 1.0
    for j in range(1, m):
        psi[j] = (theta[j - 1] if j <= len(theta) else 0.0) + sum(
            phi_star[i - 1] * psi[j - i] for i in range(1, min(j, len(phi_star)) + 1))
    return psi

# (a) ARMA(1,1): phi = 0.6, theta = 0.3, y_n = 2.0, e_n = 0.5
f = [0.6 * 2.0 + 0.3 * 0.5]
for h in range(2, 4):
    f.append(0.6 * f[-1])
psi = psi_weights([0.6], [0.3], 3)
half = 1.96 * np.sqrt(np.cumsum(psi**2))
print(np.round(f, 3), psi, np.round(half, 2))   # [1.35 0.81 0.486] [1. 0.9 0.54] [1.96 2.64 2.84]

# (b) ARIMA(0,1,1) with theta = -0.7: phi*(z) = 1 - z, so phi* = [1]
psi = psi_weights([1.0], [-0.7], 12)
print(psi[:4], round(1.96 * 2 * np.sqrt(np.sum(psi**2)), 2))   # [1. 0.3 0.3 0.3] 5.53

# (c) airline model on monthly CO2 (the data ship with statsmodels)
import statsmodels.api as sm
from statsmodels.tsa.statespace.sarimax import SARIMAX
co2 = sm.datasets.co2.load_pandas().data["co2"].resample("MS").mean().interpolate().to_numpy()
res = SARIMAX(co2, order=(0, 1, 1), seasonal_order=(0, 1, 1, 12)).fit(disp=False)
print(np.round(res.params, 3))                  # [-0.353 -0.855  0.093]
fc = res.get_forecast(36)
print(np.round(fc.predicted_mean[[0, 11, 23, 35]], 2))   # [371.98 372.54 374.13 375.72]

6. Where it shows up #

• Forecasting practice. Retail demand, call-centre volumes, energy load, tourism and macroeconomic series are forecast with (seasonal) ARIMA or with its automatic versions; statistical agencies use X-13ARIMA-SEATS, a seasonal ARIMA, to remove the season from official statistics.
• Benchmarks for ML. The M4 and M5 competitions and most deep-learning forecasting papers report AutoARIMA and ETS as baselines. Lesson 08 and 09 compare against them with the backtests of lesson 03.
• Hybrid models. Gradient boosting on the residuals of an ARIMA, or ARIMA errors in a regression on weather covariates, combine the strengths of both worlds.
• Science. Seasonal ARIMA describes CO₂ and temperature records, river discharge and epidemic surveillance counts; residuals above the upper interval flag outbreaks or sensor faults (lesson 11).

7. Common mistakes #

• Over-differencing. Differencing a stationary series introduces an MA unit root: the differenced series has $\rho(1)$ near $-0.5$ and a non-invertible MA part, and the intervals become too wide. If the ACF of the differenced series starts with a large negative spike, you may have differenced once too often.
• Comparing AICc across different $d$ or $D$. Differencing changes the data that the likelihood is computed on. Choose $d$ and $D$ first, then compare orders.
• Reading intervals as guaranteed. They ignore parameter uncertainty and model error. Check their coverage in a backtest.
• Forgetting the back-transform bias. $\exp(\text{forecast of }\log y)$ is a median. Sums of such forecasts (for example monthly forecasts added up to a yearly total) are then biased low.
• Adding a constant with $d = 2$, or by habit with $d = 1$. With $d = 1$ the constant is a drift that continues forever; include it only if the trend is real.
• Using future covariates you will not have. A regression with ARIMA errors evaluated with the observed future temperature looks better than it will be in production (leakage, lesson 03).

8. Check yourself #

1. Write ARIMA(1,1,0) with $\phi = 0.4$ as an equation for $y_t$ in terms of past $y$.
   Answer

   $(1 - 0.4B)(1 - B)y_t = \varepsilon_t$. The left side is $(1 - 1.4B + 0.4B^2)y_t$, so $y_t = 1.4y_{t-1} - 0.4y_{t-2} + \varepsilon_t$. It is an AR(2) with a unit root ($1 - 1.4z + 0.4z^2 = (1 - z)(1 - 0.4z)$).
2. For that model, compute the forecasts at $h = 1, 2$ from $y_n = 20$, $y_{n-1} = 18$.
   Answer

   $\hat y_{n+1} = 1.4 \cdot 20 - 0.4 \cdot 18 = 28 - 7.2 = 20.8$. $\hat y_{n+2} = 1.4 \cdot 20.8 - 0.4 \cdot 20 = 29.12 - 8 = 21.12$. The forecast changes are $0.8 = 0.4 \cdot 2$ and $0.32 = 0.4^2 \cdot 2$, where $2$ is the last observed change $y_n - y_{n-1}$. They decay geometrically, so the forecast levels off at $20 + 2 \cdot 0.4/(1 - 0.4) = 21.33$.
3. What are the $\psi$-weights of the random walk with drift, and the 80 % half-width at $h = 9$ for $\sigma = 1$?
   Answer

   All $\psi_j = 1$ (the drift does not change them). $\operatorname{Var}(e_{n+9}) = 9$, so the half-width is $1.28 \cdot 3 = 3.84$.
4. Why does ARIMA(0,1,1) with $\theta$ close to $-1$ behave almost like a constant mean model?
   Answer

   $(1 - B)y_t = (1 + \theta B)\varepsilon_t$ with $\theta \to -1$ makes both sides share the factor $(1 - B)$, which cancels: $y_t \approx c + \varepsilon_t$. In SES terms, $\alpha = 1 + \theta \to 0$: the level hardly moves, so the forecast is close to the long-run average.
5. A monthly series shows, after $\nabla\nabla_{12}$, a single significant ACF spike at lag 12 and PACF spikes at lags 12, 24 and 36 that shrink. Which seasonal term do you add?
   Answer

   A seasonal MA(1): $Q = 1$, the factor $(1 + \Theta B^{12})$. A cut-off in the ACF at the seasonal lag with a decaying seasonal PACF is the MA signature, one season level up.
6. You fit a model to $\log y$ and get $\hat w = 5$ with $\sigma_h^2 = 0.04$. Give the median and the mean forecast of $y$.
   Answer

   Median $e^5 = 148.4$. Mean $e^{5 + 0.02} = 151.4$ (the approximation $e^5(1 + 0.02)$ also gives 151.4).

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the ARIMA($p, d, q$) model with the backshift operator.

   Show answer

2. Which ARIMA model is simple exponential smoothing with parameter $\alpha$?

   Show answer

3. When you forecast with an ARIMA model, what do you put in for future shocks and for past shocks?

   Show answer

4. Give the $h$-step forecast error variance of an ARIMA model.

   Show answer

5. How does the 95 % interval of a random walk grow with the horizon $h$?

   Show answer

6. Compute the forecast $\hat y_{n+2\mid n}$ of an AR(1) with $\mu = 10$, $\phi = 0.5$, $y_n = 14$.

   Show answer

7. What is the airline model, and at which lags does it have MA terms?

   Show answer

8. What does a seasonal difference with $D = 1$ do?

   Show answer

9. You forecast $\log y$. Is $\exp(\hat w)$ the mean or the median forecast of $y$?

   Show answer

10. What do you need to forecast with a regression with ARIMA errors?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Explain ARIMA in one minute." — Difference the series $d$ times to make it stationary, model the differences as ARMA($p, q$), and sum back to forecast. Choose $d$ by unit-root tests, $p$ and $q$ by AICc, and check the residuals.
• "Why do ARIMA intervals widen?" — The forecast error is a sum of the future shocks weighted by the $\psi$-weights. With differencing, every shock has a permanent effect, so the variance grows linearly with the horizon.
• "ARIMA or exponential smoothing?" — They overlap (SES = ARIMA(0,1,1), and some Holt models are ARIMA too), but neither contains the other. Choose by a backtest, not by habit.
• "How do you handle multiple seasonalities or daily data with a yearly cycle?" — Fourier terms in a regression with ARIMA errors, or a state-space / ML model; not a SARIMA with $m = 365$.
• "Your 95 % intervals cover only 85 % in the backtest. Why?" — Parameter uncertainty, model misspecification, non-Gaussian or heteroscedastic errors. Fixes: bootstrap residuals, GARCH, or conformal calibration.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. arma_forecast(y, resid, phi, theta, c, h): the forecast recursion, checked on AR(1) against $\mu + \phi^h(y_n - \mu)$ and on the worked example (a).
2. psi_weights(phi_star, theta, m): the recursion, checked on AR(1) ($\phi^j$) and MA(1).
3. forecast_intervals for ARIMA(0,1,1) and their coverage by simulation.
4. Fit a seasonal AR model by least squares on $\nabla\nabla_{12}$ CO₂, undo the differences to forecast, and compare it with the seasonal naive forecast on the last two years.
5. SES = ARIMA(0,1,1): the same one-step forecasts with $\theta = \alpha - 1$.

An optional cell fits the airline model with statsmodels SARIMAX if it is installed.

11. Watch and read #

Time Series Talk : ARIMA Modelritvikmath · 9 min

Differencing plus ARMA, as in section 3.

Time Series Talk : Seasonal ARIMA Modelritvikmath · 12 min

The seasonal polynomials and the multiplicative form.

What are ARIMA ModelsAric LaBarr · 5 min

A short overview of the whole ARIMA idea.

GARCH Model : Time Series Talkritvikmath · 10 min

For the GARCH block: a model for changing variance.

• Hyndman and Athanasopoulos, Forecasting: Principles and Practice, 3rd ed. (free online): §9.5 Non-seasonal ARIMA models (the constant and $d$ table), §9.7 ARIMA modelling in fable (the workflow and the Hyndman–Khandakar algorithm), §9.8 Forecasting, §9.9 Seasonal ARIMA models, §9.10 ARIMA vs ETS, §5.6 Forecasting using transformations, and chapter 10 on regression with ARIMA errors and dynamic harmonic regression.
• R. J. Hyndman and Y. Khandakar (2008), "Automatic time series forecasting: the forecast package for R", Journal of Statistical Software 27(3): the auto.arima algorithm.
• Shumway and Stoffer, Time Series Analysis and Its Applications, 4th ed., Springer 2017, ch. 3.7–3.9: forecasting, integrated and seasonal ARIMA, with proofs.
• Robert Nau (Duke University), Introduction to ARIMA models: free notes with practical rules for identifying orders.
• T. Bollerslev (1986), "Generalized autoregressive conditional heteroskedasticity", Journal of Econometrics 31: the GARCH model.