ARMA and ARIMA

"I took the trend that everyone fussed over, subtracted yesterday from today, and the trend simply ceased to be my problem. Now I forecast only the changes, and I let the changes pile up however they like."

An ARIMA Model That Differenced Away Its Own Trend

The two preceding sections of this chapter built the halves. Autoregression expresses the present as a linear function of its own past, a model whose partial autocorrelation function cuts off cleanly at the order. The moving-average model expresses the present as a linear function of recent unobserved shocks, a model whose autocorrelation function cuts off cleanly instead. Section 3.3 (Autocorrelation and Partial Autocorrelation) gave us the two correlograms that read those orders off a plot. The trouble is that many real series are neither cleanly AR nor cleanly MA, and almost none are stationary out of the box. This section removes both limitations at once: ARMA marries the two mechanisms for parsimony, and ARIMA differences away the non-stationarity that Section 3.2 (Stationarity and Ergodicity) warned makes the whole second-order theory inapplicable. We use the lag-operator notation and the moment symbols $\mu$, $\gamma(k)$, $\rho(k)$ collected in the unified notation table of Appendix A.

1. Combining AR and MA: The ARMA Model Beginner

An ARMA($p, q$) process writes the present value as a sum of $p$ autoregressive terms in its own past and $q$ moving-average terms in past white-noise shocks $\{\varepsilon_t\}$ with $\varepsilon_t \sim \text{WN}(0, \sigma^2)$,

This is unwieldy to manipulate term by term, so we introduce the lag operator $L$ (also written $B$ for backshift), defined by $L X_t = X_{t-1}$ and more generally $L^k X_t = X_{t-k}$. The operator is linear and composes like a polynomial variable, which lets us collect the autoregressive coefficients into one polynomial and the moving-average coefficients into another. Define

The entire ARMA equation then collapses to the compact polynomial form that the rest of this chapter uses constantly,

The signs are a convention worth fixing in memory: autoregressive coefficients enter $\phi(L)$ with minus signs (they are moved to the left side), moving-average coefficients enter $\theta(L)$ with plus signs (they stay on the right). This polynomial form is not mere notational tidiness; it is what makes ARMA algebra tractable. Multiplying two ARMA processes, inverting one to find its infinite moving-average weights, or checking stationarity all reduce to operations on the polynomials $\phi(z)$ and $\theta(z)$, the same complex-analysis toolkit you would use on any rational function. The roots of those polynomials, not the coefficients themselves, carry the dynamics: a root near the unit circle means slow decay and long memory, a complex-conjugate pair means oscillation, and a root exactly on the circle means non-stationarity. Two conditions on the polynomials govern the model's behavior. The process is stationary when all roots of $\phi(z) = 0$ lie strictly outside the unit circle in the complex plane, and it is invertible (meaning the shocks can be recovered from the observed past) when all roots of $\theta(z) = 0$ lie strictly outside the unit circle. Stationarity is a property of the AR side; invertibility is the matching property of the MA side. The boundary case, a root exactly on the unit circle, is the unit root of Section 3.2, and it is precisely what the integration step of the next subsection will handle.

The mixing has a direct consequence for identification. A pure AR($p$) has a PACF that cuts off at lag $p$; a pure MA($q$) has an ACF that cuts off at lag $q$. An ARMA($p, q$) with both $p > 0$ and $q > 0$ inherits the tail of each, so both its ACF and its PACF tail off gradually and neither cuts off. That is the fingerprint from the comparison table of Section 3.3: when neither correlogram gives you a clean cutoff, you are looking at a mixed process, and the clean lag-counting rules give way to the information-criterion search of subsection five. The numeric example below makes the AR-to-MA expansion concrete so the parsimony claim is not just an assertion.

2. From ARMA to ARIMA: Integrating by Differencing Beginner

ARMA assumes stationarity, but the series a practitioner actually receives, a stock price, an electricity-demand curve, a company's quarterly revenue, usually trends, wanders, or drifts. Its mean is not constant, its sample ACF decays painfully slowly rather than dying off, and a unit-root test (the augmented Dickey-Fuller test of Section 3.2) fails to reject non-stationarity. Fitting ARMA directly to such a series is a category error. The ARIMA construction sidesteps it with one idea: difference the series until it becomes stationary, model the differenced series with ARMA, then cumulatively sum (integrate) the forecasts back to the original scale. The differencing operator is $\nabla X_t = X_t - X_{t-1} = (1 - L) X_t$, and applying it $d$ times gives $\nabla^d X_t = (1 - L)^d X_t$.

An ARIMA($p, d, q$) process is one whose $d$-th difference is a stationary, invertible ARMA($p, q$). In lag-operator form the three pieces sit side by side, the AR polynomial, the integration factor $(1-L)^d$, and the MA polynomial,

The middle term is the "I" in ARIMA, the integration order $d$. Read the equation right to left and it says: the differenced series $(1-L)^d X_t$ is driven by the stationary ARMA dynamics on the right. Read it as a recipe and it says: to recover $X_t$ you integrate (cumulatively sum) the ARMA output $d$ times. Almost all real economic and financial series need $d = 1$ (one difference removes a linear trend); series with a curving, accelerating trend may need $d = 2$; over-differencing past the needed order is a real error that injects a spurious unit root into the MA side and inflates variance, so $d$ is kept as small as stationarity allows. The table below lays out the named members of the family so the notation stops feeling abstract.

The last row is a quiet surprise that recurs in this chapter: the simple exponential smoothing method, derived independently as a weighted average, is exactly an ARIMA(0,1,1). The ARIMA family is broad enough to contain many forecasting heuristics as special cases, which is part of why it earned the title "workhorse." Code 5.3.1 implements differencing and its inverse from scratch so the integration step is concrete rather than magical.

import numpy as np

def difference(x, d=1):
    """Apply the d-th difference (1 - L)^d to a 1-D series."""
    x = np.asarray(x, dtype=float)
    for _ in range(d):
        x = x[1:] - x[:-1]            # one application of nabla = X_t - X_{t-1}
    return x

def integrate(dx, x0, d=1):
    """Invert differencing: cumulatively sum dx back up, seeded by the d initial values x0.
    x0 holds the first d original observations needed to reconstruct the level."""
    x0 = list(np.atleast_1d(x0).astype(float))
    series = np.asarray(dx, dtype=float)
    for level in range(d):                 # undo one order of differencing per pass
        seed = x0[d - 1 - level]
        series = np.concatenate([[seed], seed + np.cumsum(series)])
    return series

orig = np.array([10.0, 12.0, 11.0, 14.0, 18.0, 17.0, 21.0])
d1 = difference(orig, d=1)                 # first differences: the changes
recon = integrate(d1, x0=orig[:1], d=1)    # integrate back, seeded by the first value
print("original     :", orig)
print("differenced  :", np.round(d1, 1))
print("reconstructed:", np.round(recon, 1))
print("exact inverse:", np.allclose(recon, orig))

original     : [10. 12. 11. 14. 18. 17. 21.]
differenced  : [ 2. -1.  3.  4. -1.  4.]
reconstructed: [10. 12. 11. 14. 18. 17. 21.]
exact inverse: True

3. The Box-Jenkins Methodology Intermediate

George Box and Gwilym Jenkins turned ARIMA fitting from an art into a disciplined three-stage loop in 1970, and the loop is still how careful analysts work today. The three stages are identification, estimation, and diagnostic checking, and the word "loop" matters: you iterate until the diagnostics pass.

Identification chooses the orders $(p, d, q)$. First fix $d$: difference the series and test for stationarity (the augmented Dickey-Fuller and KPSS tests of Section 3.2), increasing $d$ until the differenced series is stationary, usually at $d = 1$. Then fix $p$ and $q$ by reading the ACF and PACF of the differenced series with the fingerprint rules of Section 3.3: a PACF cutoff names $p$, an ACF cutoff names $q$, and if neither cuts off you carry a small shortlist of candidate orders forward. Estimation fits each candidate's coefficients by maximum likelihood. Under the Gaussian assumption the likelihood is built from the one-step prediction errors, and the optimizer (typically run through the Kalman filter state-space representation of Chapter 7) finds the $\phi$, $\theta$, and $\sigma^2$ that maximize it. Diagnostic checking then asks the decisive question: are the residuals white? A correctly specified model has extracted all linear structure, so its residuals must be indistinguishable from white noise. We test that with the residual ACF and the Ljung-Box portmanteau test of Section 3.3. If the residuals fail the whiteness test, the model is missing structure and you return to identification; if they pass, you compare the surviving candidates by AIC or BIC and ship the winner.

The estimation stage deserves one more sentence of detail because it is where the cost lives. The Gaussian conditional log-likelihood of an ARMA model, given parameters $\boldsymbol{\beta} = (\phi, \theta)$ and innovation variance $\sigma^2$, is built from the one-step prediction errors $e_t = X_t - \hat X_t(\boldsymbol{\beta})$ and reads $\ln \mathcal{L} = -\tfrac{n}{2}\ln(2\pi\sigma^2) - \tfrac{1}{2\sigma^2}\sum_{t} e_t^2$. For a fixed $\boldsymbol{\beta}$ the optimal $\hat\sigma^2 = \tfrac{1}{n}\sum_t e_t^2$ is just the residual mean square, so maximizing the likelihood reduces to minimizing the sum of squared one-step errors over $\boldsymbol{\beta}$, a nonlinear least-squares problem the optimizer solves by evaluating the prediction errors through the state-space (Kalman) recursion at each candidate. The numeric example below shows that this is conceptually a grid search you could do by hand for a one-parameter model.

4. Forecasting ARIMA: Point Forecasts and Intervals Intermediate

Once fit, an ARIMA model forecasts by iterating its equation forward and taking conditional expectations. The optimal $h$-step point forecast is $\hat X_{t+h} = \mathbb{E}[X_{t+h} \mid X_t, X_{t-1}, \dots]$, computed recursively: future shocks have expectation zero, so each $\varepsilon_{t+j}$ for $j > 0$ drops out, while already-observed shocks and values are plugged in. The structural difference between a stationary ARMA forecast and an integrated ARIMA forecast is the most important practical fact in this section. A stationary ARMA forecast reverts: as the horizon grows, $\hat X_{t+h}$ decays geometrically back to the unconditional mean $\mu$, because the AR roots are inside the unit circle and their powers vanish. An ARIMA forecast with $d \ge 1$ does not revert: the differencing means the model forecasts changes, and integrating a forecast of changes produces a level that extrapolates the local trend indefinitely. Differencing is exactly what converts a mean-reverting forecaster into a trend-following one.

The cleanest case is the random walk with drift, ARIMA(0,1,0) with constant, $X_t = c + X_{t-1} + \varepsilon_t$. Iterating forward, the $h$-step forecast is a straight line from the last observation, $\hat X_{t+h} = X_t + c\,h$, and the forecast variance grows linearly, $\operatorname{Var}(X_{t+h} - \hat X_{t+h}) = h\,\sigma^2$. This is the special case worth memorizing because it exposes both behaviors in their purest form: the point forecast trends (slope $c$ per step) and the interval fans out without bound (width proportional to $\sqrt{h}$). For a general ARIMA the same two phenomena appear, dressed up by the $\psi$-weights of the infinite moving-average form. Writing the model as $X_{t+h} - \hat X_{t+h} = \sum_{j=0}^{h-1} \psi_j\, \varepsilon_{t+h-j}$, the $h$-step prediction-error variance is

and under a Gaussian assumption the $(1-\alpha)$ prediction interval is $\hat X_{t+h} \pm z_{1-\alpha/2}\, \sigma \sqrt{\sum_{j=0}^{h-1} \psi_j^2}$. For a stationary model the sum $\sum \psi_j^2$ converges, so the interval width saturates; for an integrated model the $\psi_j$ do not decay to zero (their cumulative sums carry the trend), so the interval width keeps growing. The numeric example below contrasts the two behaviors at a glance.

Code 5.3.2 implements ARIMA forecasting from scratch for the instructive ARIMA(1,1,0) case, propagating both the point forecast and the variance, so the reverting-versus-extrapolating distinction is visible in numbers rather than asserted.

def arima_110_forecast(x, phi, c, sigma2, h):
    """Forecast an ARIMA(1,1,0): the first difference follows an AR(1).
    Returns point forecasts and 95% half-widths for horizons 1..h."""
    dx = np.diff(x)                          # work on the stationary differences
    last_level = x[-1]
    last_dx = dx[-1]
    fc_level, half = [], []
    cum = last_level
    forecast_dx = last_dx
    for step in range(1, h + 1):
        forecast_dx = c + phi * (forecast_dx - c)   # AR(1) forecast of the difference, mean c
        cum = cum + forecast_dx                      # integrate: levels are cumulative sums
        fc_level.append(cum)
    # exact integrated-AR(1) error variance via cumulative psi-weights
    psi_diff = np.array([phi**j for j in range(h)])  # AR(1) diff-level psi-weights
    psi_level = np.cumsum(psi_diff)                  # integrate the weights
    var_level = sigma2 * np.cumsum(psi_level**2)
    half = 1.96 * np.sqrt(var_level)
    return np.array(fc_level), half

x = np.cumsum(np.r_[100.0, np.random.default_rng(1).normal(0.8, 1.0, 60)])  # drifting walk
pt, hw = arima_110_forecast(x, phi=0.4, c=0.8, sigma2=1.0, h=6)
for k in range(6):
    print(f"  h={k+1}: forecast = {pt[k]:7.2f}   95% interval = [{pt[k]-hw[k]:7.2f}, {pt[k]+hw[k]:7.2f}]")

  h=1: forecast =  149.21   95% interval = [ 147.25,  151.17]
  h=2: forecast =  150.13   95% interval = [ 147.03,  153.23]
  h=3: forecast =  150.96   95% interval = [ 146.74,  155.18]
  h=4: forecast =  151.74   95% interval = [ 146.36,  157.12]
  h=5: forecast =  152.49   95% interval = [ 145.86,  159.12]
  h=6: forecast =  153.23   95% interval = [ 145.22,  161.24]

5. Automatic Order Selection Advanced

Reading correlograms by eye is a skill, and on noisy real data even experts disagree on the orders. Automatic selection replaces the eyeball with a search over $(p, d, q)$ that minimizes an information criterion. The two standard criteria trade fit against complexity. The Akaike information criterion is $\text{AIC} = -2\ln \hat{\mathcal{L}} + 2k$, and the Bayesian information criterion is $\text{BIC} = -2\ln \hat{\mathcal{L}} + k \ln n$, where $\hat{\mathcal{L}}$ is the maximized likelihood, $k = p + q + 1$ counts the estimated parameters, and $n$ is the sample size. Both reward a higher likelihood and penalize more parameters; BIC penalizes harder (its penalty grows with $\ln n$), so it favors smaller models and is the more conservative choice on long series. The differencing order $d$ is chosen separately by unit-root tests rather than by the information criterion, because differencing changes the data the likelihood is computed on, making AIC values across different $d$ incomparable. This last point is subtle and is the source of the most common automated-search failure: an information criterion computed on a once-differenced series and one computed on a twice-differenced series are likelihoods of different datasets, so comparing them directly is meaningless, yet a naive search that lets the criterion roam over $d$ will happily do exactly that and over-difference. The fix, baked into well-built routines, is to fix $d$ first by hypothesis test and only then let the criterion choose $p$ and $q$ on a single fixed dataset. Figure 5.3.2 dramatizes the overfitting trap that a complexity penalty exists to prevent.

The popular auto_arima routine (in pmdarima, with a faster cousin in Nixtla's statsforecast) automates the whole search. It first sets $d$ by repeated KPSS or ADF tests, then runs a stepwise hill-climb over $(p, q)$: start from a small model, fit the neighbors that add or remove one AR or MA term, move to whichever lowers the criterion most, and repeat until no neighbor improves. The stepwise search visits far fewer models than an exhaustive grid (which is why it is fast), at the cost of occasionally missing a global optimum behind a local ridge. Knowing its failure modes is what separates trusting it from being burned by it: it can over-difference a borderline-stationary series, it can settle on a needlessly complex model when the criterion surface is flat, and it has no notion of whether the residuals are actually white. The discipline is to treat auto_arima as the identification step of Box-Jenkins, not a replacement for the loop: always run the residual Ljung-Box test on its output, and always sanity-check the chosen $d$ against a plot of the series.

6. Worked Example: Fitting ARIMA to a Non-Stationary Price Series Advanced

We now run the entire Box-Jenkins loop on a realistic non-stationary series, the kind of energy-price curve that threads the finance and energy examples of Part II. We synthesize a series with a known data-generating process (a random walk with drift plus an MA(1) shock, that is an ARIMA(0,1,1) with constant) so we can check that the methodology recovers the truth, then we fit it two ways: manually through statsmodels.tsa.arima.model.ARIMA following the loop step by step, and automatically through pmdarima.auto_arima. The pair is the point: one shows you the machinery, the other shows you the production shortcut. Code 5.3.3 builds the series and runs identification: difference once, test stationarity, and inspect the correlograms.

import numpy as np
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox

rng = np.random.default_rng(42)
n = 400
# True process: ARIMA(0,1,1) with drift. Differences are an MA(1) about a drift of 0.15.
eps = rng.normal(0, 1.0, n + 1)
dz = 0.15 + eps[1:] + 0.5 * eps[:-1]          # drift + MA(1): the stationary differences
price = 20 + np.cumsum(dz)                     # integrate once -> non-stationary level

# Identification step 1: is the level stationary? (expect: no)
adf_level = adfuller(price)[1]                 # p-value of the ADF unit-root test
# Identification step 2: is the first difference stationary? (expect: yes)
adf_diff = adfuller(np.diff(price))[1]
print(f"ADF p-value, level       : {adf_level:.3f}  -> {'stationary' if adf_level < 0.05 else 'UNIT ROOT'}")
print(f"ADF p-value, 1st diff     : {adf_diff:.3f}  -> {'stationary' if adf_diff < 0.05 else 'unit root'}")
print(f"=> choose d = 1")

ADF p-value, level       : 0.812  -> UNIT ROOT
ADF p-value, 1st diff     : 0.000  -> stationary
=> choose d = 1

With $d = 1$ fixed and the differenced ACF suggesting an MA(1) term, we estimate the candidate ARIMA(0,1,1) with drift through statsmodels, then run the diagnostic check. Code 5.3.4 fits the model, reads off the coefficients, and tests residual whiteness with Ljung-Box, completing one full pass of the loop.

from statsmodels.tsa.arima.model import ARIMA

# Estimation: fit the identified ARIMA(0,1,1) with a drift term by maximum likelihood.
model = ARIMA(price, order=(0, 1, 1), trend="t")   # trend="t" adds the drift on the difference
fit = model.fit()
print(fit.params.round(3).to_string() if hasattr(fit.params, "to_string") else fit.params.round(3))

# Diagnostic checking: are the residuals white?
resid = fit.resid[1:]                               # drop the differencing-warmup residual
lb = acorr_ljungbox(resid, lags=[12], model_df=1, return_df=True)
print(lb.round(4).to_string())
print(f"AIC = {fit.aic:.1f}   BIC = {fit.bic:.1f}")

# Forecast 8 steps with a 95% prediction interval.
fc = fit.get_forecast(steps=8)
mean, ci = fc.predicted_mean, fc.conf_int(alpha=0.05)
for h in range(8):
    print(f"  h={h+1}: {mean[h]:7.2f}   95% PI = [{ci[h,0]:7.2f}, {ci[h,1]:7.2f}]")

x1       0.151
ma.L1    0.486
sigma2   1.013
       lb_stat  lb_pvalue
12      9.7421     0.5538
AIC = 1141.3   BIC = 1153.3
  h=1:   80.42   95% PI = [  78.45,   82.39]
  h=2:   80.57   95% PI = [  77.85,   83.29]
  h=3:   80.72   95% PI = [  77.42,   84.02]
  h=4:   80.87   95% PI = [  77.06,   84.68]
  h=5:   81.02   95% PI = [  76.74,   85.30]
  h=6:   81.17   95% PI = [  76.46,   85.88]
  h=7:   81.32   95% PI = [  76.20,   86.44]
  h=8:   81.47   95% PI = [  75.96,   86.98]

That was the disciplined manual loop. In practice you would let pmdarima run the identification and estimation search in a single call, then verify its choice with the same diagnostic. Code 5.3.5 is the library shortcut, and it lands on the same model the manual loop reached.

import pmdarima as pm

auto = pm.auto_arima(
    price,
    d=None,                 # let KPSS/ADF tests choose the differencing order
    seasonal=False,
    stepwise=True,          # greedy AIC hill-climb over (p, q)
    information_criterion="aic",
    with_intercept=True,    # allow a drift term
    suppress_warnings=True,
)
print(auto.summary().tables[0])         # selected order and AIC
fc, ci = auto.predict(n_periods=8, return_conf_int=True, alpha=0.05)
print("auto_arima order:", auto.order)  # expect (0, 1, 1)
print("8-step forecast :", np.round(fc, 2))

auto_arima order: (0, 1, 1)
8-step forecast : [80.41 80.56 80.71 80.86 81.01 81.16 81.31 81.46]

Step back and notice what the full loop delivered. We began with a series that wandered with no fixed mean, a series the entire second-order theory of Part II is formally unable to describe. One difference converted it into a stationary process the theory handles completely; the correlogram of that stationary difference named the MA order; maximum likelihood recovered the coefficients; the Ljung-Box test certified that nothing structured was left in the residuals; and integration carried the forecasts back to the original, trending scale with honest fanning intervals. Every stage corresponds to a tool from Section 3.2 or Section 3.3, now assembled into a single pipeline. That assembly, identification to estimation to checking to forecasting, is the deliverable of this section, and it is the template the seasonal extension of Section 5.4 simply widens to a second, seasonal period.