Moving Average Models

"Give me a shock and I will repeat it for exactly $q$ steps, no more. After that I forget completely and return to my quiet average, pretending nothing ever happened. People mistake me for the smoother who lives next door; we are not the same."

A Moving Average Still Processing Last Week's Shocks

Chapter 3 gave us the language of stationarity and the correlogram; Section 3.3 in particular showed that a moving-average process leaves a sharp cut in its autocorrelation function and a gradual decay in its partial autocorrelation, the exact mirror of the autoregressive fingerprint. This section turns that diagnostic observation into a working model. We assume the weak-stationarity machinery of Section 3.2, the autocovariance notation $\gamma(k)$ and autocorrelation $\rho(k)$ of Section 3.3, and the unified symbol table of Appendix A, where the innovation sequence $\varepsilon_t$ and the lag operator $B$ are defined. The MA model is also the first half of the ARMA synthesis that the rest of this chapter assembles, and through the lag-operator algebra it is the classical ancestor of the linear-attention forecasters of Chapter 13, a thread we pick up at the end.

1. Two Senses of "Moving Average" Beginner

Before any mathematics, we must separate two ideas that share a name and almost nothing else. The first is the moving-average smoother, a deterministic filter you already met among the convolutional smoothers of Chapter 4. Given an observed series $x_1, \dots, x_n$, the simple moving average of window $w$ replaces each point by the mean of its neighbors,

and a weighted moving average attaches weights $a_j$ to those neighbors, $s_t = \sum_j a_j x_{t-j}$ with $\sum_j a_j = 1$. This is a transformation applied to data you can see: it averages observations, its purpose is to suppress high-frequency noise and expose trend, and it has no probabilistic content whatsoever. It is a low-pass filter, nothing more.

The second idea, the subject of this section, is the moving-average model MA($q$), a stochastic process that constructs the present value as a weighted sum of the last $q$ unobserved random shocks plus the current shock. The two objects differ on every axis that matters, and the comparison table below is the one to internalize before reading further, because nearly every confusion about MA models dissolves once the smoother and the model are kept apart.

The collision of names is historical: both average a sliding window, so both earned the adjective "moving." But the smoother averages what you observe to clean it, while the model averages what you cannot observe to explain what you do. As the comparison table records, the smoother needs no fitting (you pick a window), whereas the model demands the nonlinear estimation of subsection four precisely because its inputs are hidden. Keep the table's middle column in Chapter 4 and its right column here, and the rest of this section follows cleanly.

2. The MA($q$) Model and Its Always-Stationary Property Beginner

Let $\{\varepsilon_t\}$ be a white-noise innovation sequence: mean zero, constant variance $\sigma^2$, and uncorrelated across time, $\operatorname{Cov}(\varepsilon_t, \varepsilon_s) = 0$ for $t \ne s$. The moving-average model of order $q$ writes the observed series as the current shock plus a fixed linear combination of the previous $q$ shocks, around a level $\mu$,

Using the lag operator $B$ (defined by $B^k \varepsilon_t = \varepsilon_{t-k}$ in Appendix A) this compresses to $X_t = \mu + \theta(B)\varepsilon_t$, where $\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$ is the moving-average polynomial. The coefficient on $\varepsilon_t$ is fixed at one; any common scale is absorbed into $\sigma^2$, and $\mu$ is the process mean because every innovation has mean zero. The defining feature is finite memory: $X_t$ contains $\varepsilon_{t-q}$ but not $\varepsilon_{t-q-1}$, so a shock that arrives today influences the series for exactly $q$ steps and then vanishes from it completely. This is the literal meaning of the epigraph, pictured in Figure 5.1.2.

The first structural reward of finite memory is that an MA($q$) process is stationary for any choice of the parameters $\theta_1, \dots, \theta_q$, with no constraints whatsoever. This is the sharpest contrast with the AR model, whose stationarity requires its characteristic roots to lie outside the unit circle. The reason is immediate: a finite sum of stationary, uncorrelated terms is itself stationary, term by term. To see it, compute the moments directly. The mean is constant,

and the variance, using $\theta_0 = 1$ and the fact that distinct innovations are uncorrelated so all cross terms drop, is finite and time-invariant,

The autocovariance at lag $k$ keeps only the products of innovations that the two windows share. Writing $X_t - \mu = \sum_{j=0}^q \theta_j \varepsilon_{t-j}$ and $X_{t+k} - \mu = \sum_{i=0}^q \theta_i \varepsilon_{t+k-i}$, a shared innovation requires $t - j = t + k - i$, that is $i = j + k$, which is possible only when $0 \le k \le q$. Hence

Dividing by $\gamma(0)$ gives the autocorrelation $\rho(k) = \gamma(k)/\gamma(0)$, which is exactly zero for every lag beyond $q$. This is the celebrated ACF cutoff: the autocorrelation function of an MA($q$) process is nonzero up to lag $q$ and then drops to zero and stays there, the precise mirror of the gradual PACF decay catalogued in Section 3.3. The cutoff is what lets you read $q$ straight off a correlogram: the order of an MA model is the last lag at which the sample ACF pierces the $\pm 1.96/\sqrt{n}$ band.

3. Invertibility and the Infinite-AR Representation Intermediate

The numeric example exposed a problem: two different MA(1) models, $\theta_1$ and $1/\theta_1$, produce identical autocorrelation structure and therefore identical second-order behavior. Since a Gaussian stationary process is fully described by its mean and autocovariance, these two parameterizations are observationally indistinguishable from data. We need a convention that picks one of each such pair, and the right convention is invertibility. An MA($q$) process is invertible when the roots of its characteristic polynomial

all lie strictly outside the unit circle, $|z| > 1$. For the MA(1) the single root is $z = -1/\theta_1$, which lies outside the unit circle exactly when $|\theta_1| < 1$. Of the pair $\{0.8, 1.25\}$ this selects $\theta_1 = 0.8$, and in general it selects the unique parameter set whose roots are all outside the circle. Invertibility is to MA models what stationarity-of-the-roots is to AR models: a convention that resolves identifiability and guarantees a well-behaved representation.

Why "invertible"? Because the condition is precisely what lets us solve the model for the unobserved shock $\varepsilon_t$ in terms of the observed past, turning the MA model inside out into an infinite autoregression. Formally, when all roots lie outside the unit circle the polynomial $\theta(B)$ has a convergent inverse $\theta(B)^{-1} = \pi(B) = 1 - \pi_1 B - \pi_2 B^2 - \cdots$, and applying it to $X_t - \mu = \theta(B)\varepsilon_t$ yields

which rearranges into the infinite-AR representation $X_t = \mu + \sum_{j=1}^{\infty} \pi_j (X_{t-j} - \mu) + \varepsilon_t$. An MA($q$) process is thus an AR($\infty$) process in disguise; this is the precise sense in which the ACF cutoff of the MA mirrors the PACF tail, since the AR($\infty$) form gives every lag a small direct weight. For the MA(1) the inverse expands geometrically, $\pi(B) = (1 + \theta_1 B)^{-1} = 1 - \theta_1 B + \theta_1^2 B^2 - \cdots$, so $\pi_j = -(-\theta_1)^j$. This series converges if and only if $|\theta_1| < 1$, which is exactly the invertibility condition: the non-invertible twin $\theta_1 = 1.25$ would make the $\pi_j$ explode, and the model could not be written as a sensible function of its own past.

4. Estimation and Forecasting Advanced

Because the innovations are latent, an MA model cannot be fit by ordinary least squares. The standard route is conditional sum of squares (CSS), which turns the invertible representation of subsection three into a recursive reconstruction of the shocks. Fix candidate parameters $(\mu, \theta_1, \dots, \theta_q)$, set the unobservable pre-sample innovations to their expectation zero ($\varepsilon_0 = \varepsilon_{-1} = \cdots = 0$, the "conditional" assumption), and then run the model forward solved for $\varepsilon_t$,

Each reconstructed $\hat\varepsilon_t$ depends on the parameters through every earlier $\hat\varepsilon$, so the residuals are a nonlinear function of $\theta$. The CSS estimator minimizes their sum of squares, $S(\theta) = \sum_{t=1}^n \hat\varepsilon_t^2$, over the invertible region, using a numerical optimizer (Gauss-Newton, L-BFGS, or Nelder-Mead) rather than a closed-form normal equation. Under Gaussian innovations CSS is a close approximation to full maximum likelihood: the exact Gaussian log-likelihood is

where $\mathbf{R}(\theta)$ is the correlation matrix built from the MA autocovariances of subsection two. Exact MLE evaluates this via the Kalman filter (the state-space recursion of Chapter 7) for an $O(n)$ likelihood; CSS drops the $\ln|\mathbf{R}|$ determinant term and the pre-sample correction, which is negligible for large $n$. The practical message is that MA fitting is iterative and can land in local minima, so good starting values (often from the method-of-moments inversion of the $\rho(k)$ formulas) and the invertibility constraint both matter.

Forecasting an invertible MA($q$) is where finite memory shows its hand most vividly. The minimum-mean-squared-error forecast is the conditional expectation $\hat X_{t}(h) = \mathbb{E}[X_{t+h} \mid X_t, X_{t-1}, \dots]$. Future innovations have expectation zero, and innovations up to time $t$ are recoverable by the invertible recursion, so

The sum runs only over innovations dated at or before $t$ (those with index $t+h-j \le t$, i.e. $j \ge h$); every future shock drops out. The consequence is stark: for horizons within the memory window ($h \le q$) the forecast still carries information from the last observed shocks, but the instant the horizon exceeds $q$ every term vanishes and the forecast collapses to the unconditional mean $\mu$. An MA($q$) model has nothing to say beyond $q$ steps; it reverts to the mean and stays there. The forecast-error variance behaves symmetrically, growing from $\sigma^2$ at $h=1$ up to the full process variance $\gamma(0) = \sigma^2 \sum_{j=0}^q \theta_j^2$ for all $h > q$, since beyond the memory horizon the prediction error is just the entire process fluctuating around its mean.

5. Worked Example: Simulate, Recover, Fit, Forecast Advanced

We now put the whole section to work on the finance running dataset of Chapter 1: daily log returns of an equity index, the canonical case where MA terms earn their keep because microstructure effects (bid-ask bounce, non-synchronous trading) induce a short, sharp lag-1 autocorrelation in returns that an MA(1) captures exactly. We simulate a realistic MA(1)/MA(2) return series, confirm the ACF cutoff predicted in subsection two, fit the model from scratch with a CSS estimator, and produce the mean-reverting forecast of subsection four. Code 5.1.1 simulates the process and recovers its ACF cutoff.

import numpy as np

def simulate_ma(n, thetas, mu=0.0, sigma=1.0, rng=None, burn=200):
    """Simulate X_t = mu + e_t + sum_j theta_j e_{t-j}, an MA(q) process."""
    rng = rng or np.random.default_rng(0)
    q = len(thetas)
    e = rng.normal(0.0, sigma, size=n + burn + q)   # innovations, with burn-in
    theta = np.concatenate(([1.0], np.asarray(thetas, float)))   # [1, theta_1..theta_q]
    x = np.empty(n + burn)
    for t in range(q, n + burn):
        x[t] = mu + np.dot(theta, e[t - np.arange(q + 1)])  # current + q past shocks
    return x[burn:]

def sample_acf(x, max_lag):
    """Sample ACF rho_hat(k), k = 0..max_lag, divisor n (as in Section 3.3)."""
    x = np.asarray(x, float) - np.mean(x)
    g0 = np.dot(x, x) / len(x)
    return np.array([np.dot(x[:len(x) - k], x[k:]) / len(x) / g0
                     for k in range(max_lag + 1)])

rng = np.random.default_rng(42)
# Finance framing: a daily log-return series with a single negative microstructure lag.
returns = simulate_ma(n=1000, thetas=[0.6], mu=0.0004, sigma=0.011, rng=rng)
acf = sample_acf(returns, max_lag=6)
band = 1.96 / np.sqrt(len(returns))
print(f"95% band: +/- {band:.3f}")
for k in range(1, 7):
    flag = "  <-- significant" if abs(acf[k]) > band else ""
    print(f"  lag {k}: rho_hat = {acf[k]:+.3f}{flag}")

95% band: +/- 0.062
  lag 1: rho_hat = +0.434  <-- significant
  lag 2: rho_hat = +0.013
  lag 3: rho_hat = -0.028
  lag 4: rho_hat = +0.021
  lag 5: rho_hat = -0.006
  lag 6: rho_hat = +0.019

Having recovered the order, we fit the MA(1) from scratch. Code 5.1.2 implements the conditional-sum-of-squares estimator of subsection four directly: it reconstructs the latent innovations by the invertible recursion, sums their squares, and hands that objective to a scalar optimizer constrained to the invertible region $|\theta_1| < 1$.

from scipy.optimize import minimize_scalar

def css_objective(theta1, x, mu):
    """Conditional sum of squares for an MA(1): reconstruct shocks, return sum e^2."""
    eps_prev = 0.0                       # pre-sample innovation set to its mean, zero
    ssq = 0.0
    for xt in x:
        eps = (xt - mu) - theta1 * eps_prev   # invertible recursion solves for e_t
        ssq += eps * eps
        eps_prev = eps
    return ssq

mu_hat = returns.mean()                  # MA(q) mean is the sample mean
# Optimize CSS over the invertible interval (-1, 1); MA fitting is genuinely nonlinear.
res = minimize_scalar(css_objective, bounds=(-0.99, 0.99),
                      args=(returns, mu_hat), method="bounded")
theta1_hat = res.x
sigma2_hat = res.fun / len(returns)      # residual variance estimate
print(f"from-scratch CSS:  theta1_hat = {theta1_hat:+.3f}   "
      f"sigma2_hat = {sigma2_hat:.2e}   (true theta1 = 0.600)")

from-scratch CSS:  theta1_hat = +0.598   sigma2_hat = 1.21e-04   (true theta1 = 0.600)

Forty-odd lines of recursion and optimization above exist so the machinery is transparent. In practice the same fit, with exact maximum likelihood, standard errors, and forecasting, is three lines of statsmodels. Code 5.1.3 fits the identical model as an ARIMA(0,0,1) and produces the mean-reverting forecast.

from statsmodels.tsa.arima.model import ARIMA

# MA(1) is ARIMA(p=0, d=0, q=1); exact MLE via the Kalman filter (Chapter 7).
model = ARIMA(returns, order=(0, 0, 1), trend="c").fit()
print(model.summary().tables[1])         # const (mu) and ma.L1 (theta1) with std errors

fc = model.get_forecast(steps=5)         # forecast 5 steps ahead
print("\nstep   forecast      se")
for h, (m, s) in enumerate(zip(fc.predicted_mean, fc.se_mean), start=1):
    print(f"  {h}   {m:+.5f}   {s:.5f}")

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0004    0.000      1.07     0.285      -0.000       0.001
ma.L1          0.5996    0.025     24.18     0.000       0.551       0.648
==============================================================================

step   forecast      se
  1   +0.00185   0.01101
  2   +0.00040   0.01283
  3   +0.00040   0.01283
  4   +0.00040   0.01283
  5   +0.00040   0.01283