The ARMA Model: Derivation, Identification, and Estimation

The ARMA model is the workhorse of stationary time series analysis. Once you understand it, you understand much of what people do in applied econometrics with time-stamped data. The math is mostly polynomial algebra. Once you’ve seen it once, you can re-derive it on a napkin. The hard parts are identification (how do you know it’s an AR(2) versus an MA(3)?) and diagnostic checking (how do you know your fitted model is right?).

This article derives AR, MA, and ARMA from scratch. Hamilton (1994, Ch. 3) and Tsay (2010, Ch. 3) are the references. Brooks (2019, Ch. 5) has a more applied finance flavor if Hamilton feels too dry.

Starting Point: White Noise

We start with white noise $\varepsilon_t \sim \text{WN}(0, \sigma^2)$: zero mean, constant variance $\sigma^2$, no autocorrelation. We want to construct stationary processes from it. Two natural transformations give the two fundamental building blocks.

The Autoregressive (AR) Model

An AR($p$) process says: today’s value is a linear combination of the last $p$ values, plus a shock.

or, more compactly,

where $\Phi(L) = 1 - \phi_1 L - \phi_2 L^2 - \cdots - \phi_p L^p$ is the autoregressive polynomial in the lag operator $L$, defined by $L y_t = y_{t-1}$.

Stationarity condition: All roots of $\Phi(z) = 0$ must lie outside the unit circle in the complex plane. Equivalently, the polynomial $1 - \phi_1 z - \cdots - \phi_p z^p \neq 0$ for all $|z| \leq 1$.

If $\Phi(L)$ is invertible (which it is when the roots are outside the unit circle), we can write:

where $\mu = c / \Phi(1) = c / (1 - \phi_1 - \cdots - \phi_p)$ is the unconditional mean and the MA coefficients $\psi_j$ decay geometrically.

The Yule-Walker equations give the autocovariances in closed form. Multiply both sides of the AR equation by $y_{t-h}$ and take expectations:

Dividing through by $\gamma_0$ gives the autocorrelation recursion:

with $\rho_0 = 1$.

Special case. AR(1): $y_t = c + \phi_1 y_{t-1} + \varepsilon_t$ with $|\phi_1| < 1$.

Mean: $\mu = c / (1 - \phi_1)$. Variance: $\gamma_0 = \sigma^2 / (1 - \phi_1^2)$. Autocorrelation: $\rho_h = \phi_1^h$, which decays geometrically.

If $\phi_1 > 0$, autocorrelations are positive and decay smoothly. If $\phi_1 < 0$, they alternate in sign. Either way, the ACF tails off.

Special case. AR(2): $y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \varepsilon_t$.

The characteristic equation is $1 - \phi_1 z - \phi_2 z^2 = 0$. Stationarity requires both roots outside the unit circle. The roots determine whether the ACF is a mix of exponentials or a damped oscillation:

• Two real roots, whether positive, negative, or of opposite signs, produce a mix of two exponentials.
• Complex conjugate roots produce a damped sinusoidal oscillation.

This is the source of “pseudo-cycles” in macro data: a stationary AR(2) with complex roots looks cyclical even though it has no business cycle mechanism at all.

The Moving Average (MA) Model

An MA($q$) process says: today’s value is a linear combination of the current and last $q$ shocks.

or,

where $\Theta(L) = 1 + \theta_1 L + \cdots + \theta_q L^q$.

MA processes are always stationary for any values of $\theta_1, \ldots, \theta_q$ (because the MA representation has finite sum of $|\psi_j|$). The only restriction is invertibility: to recover the shocks from observed $y_t$, all roots of $\Theta(z) = 0$ must lie outside the unit circle.

For MA(1) with $|\theta_1| < 1$:

The ACF cuts off after lag 1: $\rho_0 = 1$, $\rho_1 = \theta_1 / (1 + \theta_1^2)$, $\rho_h = 0$ for $h \geq 2$.

This cutoff is the diagnostic signature of an MA process. It is exactly the opposite of the AR signature, where the ACF tails off.

The ARMA($p, q$) Model

Combining both:

or,

Stationarity requires all roots of $\Phi(z)$ outside the unit circle. Invertibility requires all roots of $\Theta(z)$ outside the unit circle.

The ACF of an ARMA($p, q$) tails off after lag $q$ (like an MA), and the PACF tails off after lag $p$ (like an AR). This dual behavior is what makes ARMA flexible.

Identification: the $(p, q)$ orders must be identified from the data. There is no unique ARMA representation for any given process. Adding a factor to $\Phi(L)$ and removing the same factor from $\Theta(L)$ gives the same process. The convention is to enforce causality (stationarity + invertibility) by requiring roots of both polynomials outside the unit circle. Under this restriction, the representation is unique.

The Yule-Walker Equations for ARMA

For an ARMA($p, q$), the autocovariance satisfies the extended Yule-Walker recursion for $h > q$:

For $1 \leq h \leq q$, the relation involves both $\gamma$ and the cross-covariance between $y$ and $\varepsilon$. A closed-form solution requires inverting a block-Toeplitz matrix. In practice, software handles this calculation.

Identification: ACF and PACF

The Box-Jenkins (1976) identification rules are:

The PACF (partial autocorrelation function) at lag $h$ is the correlation between $y_t$ and $y_{t-h}$ after removing the linear effect of the intermediate lags. It is computed by fitting AR($k$) models for $k = 1, 2, \ldots$ and reading off the coefficient on $y_{t-k}$. Equivalently, solve the Yule-Walker system:

The last coefficient $\phi_{h,h}$ is the PACF at lag $h$.

The intuition: PACF tells you the additional predictive content of $y_{t-h}$ beyond $y_{t-1}, \ldots, y_{t-h+1}$. If the true model is AR($p$), lags beyond $p$ add nothing.

Estimation: OLS, Yule-Walker, and MLE

Three main estimation strategies:

1. OLS (least squares). For an AR($p$), regress $y_t$ on $y_{t-1}, \ldots, y_{t-p}$ plus a constant. Standard OLS theory applies because the regressors are predetermined (lagged dependent variables) and the error $\varepsilon_t$ is conditionally mean-independent. The estimates are consistent and asymptotically normal, but the small-sample bias is well-known and can be severe (Marriott & Pope, 1954; Shaman & Stine, 1988). For an MA($q$), OLS does not work because the unobserved shocks appear as regressors.

2. Yule-Walker (method of moments). Replace theoretical autocovariances with sample autocovariances in the Yule-Walker equations and solve. This is computationally cheap but inefficient because it uses only second-moment information, not the full distribution.

3. Maximum likelihood. For Gaussian ARMA, write the exact likelihood using the innovation recursion:

with backcasted values $\hat{\varepsilon}_0, \hat{\varepsilon}_{-1}, \ldots$ initialized at zero or via unconditional moments. The log-likelihood is:

Maximize over $\theta = (c, \phi_1, \ldots, \phi_p, \theta_1, \ldots, \theta_q, \sigma^2)$. Concentrating out $\sigma^2$ gives:

and reduces the optimization to $(p + q + 1)$ parameters. The conditional MLE uses $\hat{\varepsilon}_t = 0$ for $t \leq 0$. The exact MLE uses the Kalman filter or diffuse initialization to compute the exact likelihood. It is better for near-unit-root or short series.

The likelihood surface for MA and ARMA is famously rough, with multiple local maxima, ridges, and near-flat regions near the unit circle. Box-Jenkins recommended Yule-Walker for initial values, then refinement with MLE. Most modern software uses CSS (conditional sum of squares) starting values followed by MLE.

Diagnostic Checking

After fitting, check whether the residuals look like white noise. The standard battery:

1. Residual ACF / PACF. Compute the sample autocorrelations of $\hat{\varepsilon}_t$. Under the null that the model is correct, $\hat{\rho}_h \sim N(0, 1/T)$ for $h \geq 1$. Plot with $\pm 1.96/\sqrt{T}$ confidence bands. Any lag where the autocorrelation punches through is a sign of model misspecification.

2. Ljung-Box test. The Q statistic:

Tests the joint null that the first $m$ residual autocorrelations are zero. Reject if $Q$ exceeds the critical value. Use $m \approx \log(T)$ or $m = 10$ or $m = 20$.

3. Normality. Jarque-Bera, Shapiro-Wilk on residuals. ARMA assumes Gaussianity for inference; if the residuals are heavy-tailed, you need either fat-tailed innovations or a different model class.

4. Heteroskedasticity. Apply the ARCH-LM test (Engle, 1982) to residuals. If the test rejects, your ARMA has not absorbed all the conditional heteroskedasticity. Add a GARCH layer on top.

5. Parameter stability. Recursive least squares, CUSUM tests (Brown, Durbin, Evans, 1975). If parameters drift over the sample, you have structural breaks and should consider regime-switching models.

Code: Fitting an ARMA in Python

import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA


rng = np.random.default_rng(0)
n = 1000
phi = [0.5, -0.3]
theta = [0.4]
eps = rng.standard_normal(n)
y = np.zeros(n)
for t in range(2, n):
    y[t] = 0.5*y[t-1] - 0.3*y[t-2] + eps[t] + 0.4*eps[t-1]

# Fit ARMA(2,1) by MLE
model = ARIMA(y, order=(2, 0, 1), trend='c').fit()
print(model.summary())

# Diagnostic: residual ACF and Ljung-Box
from statsmodels.stats.diagnostic import acorr_ljungbox
lb = acorr_ljungbox(model.resid, lags=[10, 20], return_df=True)
print(lb)

The fitted $\hat{\phi}_1, \hat{\phi}_2, \hat{\theta}_1$ should be near the true $(0.5, -0.3, 0.4)$. Ljung-Box should not reject at conventional levels.

What ARMA Cannot Do

ARMA is for stationary, linear, short-memory processes. It cannot handle:

• Unit roots: $y_t = y_{t-1} + \varepsilon_t$ has $\phi_1 = 1$ and is non-stationary. ARMA on first differences works (ARIMA, next article).
• Long memory: autocorrelations that decay hyperbolically rather than geometrically (ARFIMA, fractional integration).
• Nonlinearity: GARCH for variance, Markov-switching for regimes, threshold models, neural networks.
• Structural breaks: parameters changing at known or unknown dates.
• Multivariate dynamics: cross-correlations between series (VAR, VECM).

For each of these, there is a specialized model. Most of the next 11 articles are about one of them.

The Small-Sample Bias Problem

The OLS estimator for AR coefficients is biased in finite samples, downward for positive AR coefficients. The bias is approximately $-2\phi/T$ for an AR(1) with coefficient $\phi$ (Kendall, 1954; Marriott & Pope, 1954). For $T = 100$ and $\phi = 0.9$, the bias is $-0.018$, small. For $T = 30$ and $\phi = 0.9$, the bias is $-0.06$, non-trivial.

For near-unit-root processes ($\phi \to 1$), the bias is severe. Shaman & Stine (1988) show that the bias can be more than half the true coefficient for highly persistent series in small samples. The implication: if you fit an AR(1) to 30 observations of a series with true $\phi = 0.95$, the OLS estimate might be 0.4 or 0.5, very far from the truth.

Bias correction is possible (see Marriott & Pope, 1954, or the jackknife bias correction of Shaman & Stine, 1988), but the standard practice is just to acknowledge the bias and use longer samples when possible.

MLE for ARMA has the same bias asymptotically. For MA and ARMA, the bias can be even worse than for AR because the MA polynomial can be near the unit circle, and the likelihood surface has a sharp peak near the boundary.

Forecasting with ARMA

The forecast at horizon $h$ from the end of the sample is:

This is computed recursively. For an AR(2):

and so on. For $h > p$, all lagged values are themselves forecasts, so the recursion uses $\hat{y}_{T+h-j}$ for $j > p$.

The forecast error variance:

where $\hat{\psi}_j$ are the MA coefficients from the infinite MA representation. For a stationary ARMA, $\sum \hat{\psi}_j^2$ converges to $\sigma^{-2} \gamma_0$ (the unconditional variance), so the forecast variance asymptotes.

For a near-unit-root AR, the forecast variance grows like $h \sigma^2$, much faster than in a stationary process.

The 95% prediction interval:

This assumes Gaussian innovations. With heavy-tailed innovations, the interval should be wider. Bootstrap prediction intervals handle this:

1. Fit the model, get residuals $\hat{\varepsilon}_t$.
2. Resample residuals with replacement to get $\varepsilon_t^*$.
3. Simulate $\hat{y}_t^* = \hat{c} + \hat{\phi}_1 \hat{y}_{t-1}^* + \cdots + \hat{\varepsilon}_t^*$.
4. Repeat $B = 1000$ times.
5. The 2.5% and 97.5% quantiles of $\hat{y}_{T+h}^*$ give the bootstrap interval.

ARMA and the Efficient Market Hypothesis

A subtle but important point: ARMA models have nothing to say about whether asset returns are predictable. If returns are iid (white noise), the best ARMA fit is ARMA(0, 0): a constant mean plus white-noise residuals. The coefficient estimates are zero, the residual variance equals the sample variance, and forecasts are flat.

If the ARMA fit gives significant AR or MA coefficients, it suggests some predictability: past returns help forecast future returns. Whether this predictability is exploitable for profit depends on transaction costs, its magnitude, and whether it persists out of sample.

The efficient market hypothesis says returns are not predictable from public information. ARMA fitting is one way to test this. But the EMH is a statement about the joint distribution of returns and conditioning information, not just about linear predictability. Markets can be efficient (no exploitable patterns) even if ARMA coefficients are non-zero, if the predictability is too small to overcome trading costs.

Seasonal ARMA

For data with strong seasonality (monthly, quarterly), the ARMA model extends naturally. The simplest form is the multiplicative seasonal ARMA:

where $s$ is the seasonal period (12 for monthly, 4 for quarterly), and $\tilde{\Phi}, \tilde{\Theta}$ are the seasonal AR and MA polynomials.

A seasonal AR(1) at lag $s$ looks like:

This says: today’s value depends on the value $s$ periods ago. For monthly retail sales with annual cycles, this would relate January of one year to January of the next year.

Seasonal MA(1) at lag $s$:

This says: today’s value depends on a shock that happened $s$ periods ago, persisting for $s$ periods.

Identification: the ACF of a seasonal ARMA has spikes at the seasonal lags (12, 24, 36 for monthly data with annual cycles). If only lag 12 is significant, you have a seasonal MA(1). If lags 12 and 24 are significant, you might have a seasonal AR(2).

For seasonal ARIMA with differencing, the seasonal differencing $(1 - L^s)$ is added:

This is the SARIMA$(p, d, q) \times (P, D, Q)_s$ model. It handles trends (via $d$), seasonal trends (via $D$), and ARMA dynamics at both regular and seasonal lags.

For most monthly economic data, SARIMA$(0, 1, 1)(0, 1, 1)_{12}$ is a strong default. For quarterly data with strong seasonality, SARIMA$(0, 1, 1)(0, 1, 1)_4$ often fits well.

Fractional ARIMA: A Brief Note

The ARIMA class assumes $d$ is an integer. The ARFIMA (Autoregressive Fractionally Integrated Moving Average) allows non-integer $d$:

where $(1 - L)^d$ is the fractional differencing operator, defined by the binomial series:

For $0 < d < 0.5$, the process is stationary but long-memory: autocorrelations decay hyperbolically ($\rho_h \sim h^{2d-1}$) rather than geometrically. This matches empirical evidence on volatility, river flows, internet traffic.

The catch: ARFIMA is harder to estimate than ARIMA, and the asymptotic distribution of the $d$ estimator is non-standard (slow convergence to normal). For most applications, ARIMA with integer $d$ is sufficient.

References

• Box, G. E. P., & Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control. Holden-Day.
• Brooks, C. (2019). Introductory Econometrics for Finance (4th ed.). Cambridge University Press. Ch. 5-6.
• Brown, R. L., Durbin, J., & Evans, J. M. (1975). Techniques for testing the constancy of regression relationships over time. Journal of the Royal Statistical Society B, 37(2), 149-192.
• Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007.
• Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 3.
• Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303.
• Marriott, F. H. C., & Pope, J. A. (1954). Bias in the estimation of autocorrelations. Biometrika, 41(3/4), 390-402.
• Shaman, P., & Stine, R. A. (1988). The bias of autoregressive coefficient estimators. Journal of the American Statistical Association, 83(403), 842-848.
• Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 3.
• Kendall, M. G. (1954). Note on bias in the estimation of autocorrelation. Biometrika, 41(3/4), 403-404.
• Granger, C. W. J., & Joyeux, R. (1980). An introduction to long-memory time series models and fractional differencing. Journal of Time Series Analysis, 1(1), 15-29.