Estimating ARMA Models

You have written down a model — an ARMA(p, q) — and you have a stretch of data. Between the two sits the question this whole page answers: what numbers do the parameters actually take? An ARMA model is a template with holes in it,

x_t = \phi_1 x_{t-1} + \dots + \phi_p x_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q},

and estimation is the craft of choosing \hat\phi, \hat\theta and the shock variance \hat\sigma^2 so that the template best reproduces the series in front of you. Everything downstream — forecasts, intervals, diagnostics — is only as trustworthy as these estimates. There are three classic routes to them, of increasing sophistication and payoff.

Route 1 — method of moments (Yule–Walker), AR only

The oldest trick: match the moments the model predicts to the moments you measure. For a pure autoregressive model this is beautifully clean. The theoretical autocorrelations satisfy the Yule–Walker equations — a linear system tying the \phi's to the autocorrelations \rho_1, \dots, \rho_p. Replace each theoretical \rho_k with its sample version \hat\rho_k, solve the linear system, and out drop the coefficients:

\hat{\boldsymbol\phi} = \hat R^{-1}\, \hat{\boldsymbol\rho}, \qquad \hat R = \big[\hat\rho_{|i-j|}\big]_{i,j=1}^{p}.

No iteration, no starting values — just invert a p \times p matrix of sample correlations. It is fast, always returns a stationary fit, and makes a superb starting point for the fancier methods. Its fatal limitation: it works only for AR. The moment equations for a moving-average part are nonlinear in \theta, and solving them by hand is a mess — so the moment method quietly retires the moment a single MA term appears.

Route 2 — least squares

Treat the model as a regression and minimise the sum of squared one-step errors. For an AR(p) this is genuinely ordinary least squares: regress x_t on its own lags and read off the coefficients — this is least squares with the design matrix built from shifted copies of the series. We call it conditional least squares because it conditions on the first few observations (it has to start somewhere).

For an AR(1) the objective collapses to a single-variable parabola,

S(\phi) = \sum_{t=2}^{T}\big(x_t - \phi\, x_{t-1}\big)^2 = \sum x_t^2 - 2\phi \sum x_t x_{t-1} + \phi^2 \sum x_{t-1}^2,

whose minimum you can read off by setting the derivative to zero. That gives the estimator we work by hand below. When the model has MA terms the past shocks \varepsilon_{t-1}, \dots are never observed, so the sum of squares becomes a nonlinear function of the parameters and must be minimised numerically — nonlinear least squares. Exact least squares goes further and refuses to throw away the information in the first observations, reconstructing them too.

The SSE surface for an AR(1)

The estimate is wherever the sum-of-squares curve bottoms out. Below is S(\phi) = 8\phi^2 - 12\phi + 10 for a series with \sum x_t x_{t-1} = 6 and \sum x_{t-1}^2 = 8. Setting S'(\phi) = 0 gives \hat\phi = 6/8 = 0.75 — the single low point of the bowl. For a one-parameter AR this surface is a perfect parabola with one guaranteed minimum; that convexity is exactly why AR estimation is so easy and MA estimation, as we will see, is not.

Worked example — the AR(1) least-squares estimate

The formula that falls out of the parabola is worth memorising:

Concretely, with (mean-centred) sums \sum x_t x_{t-1} = 6 and \sum x_{t-1}^2 = 8, the estimate is \hat\phi = 6/8 = 0.75: a strongly persistent, mean-reverting series. Always centre the data first (subtract the sample mean), or the intercept leaks into the slope.

Route 3 — maximum likelihood, and why it wins

The gold standard is maximum likelihood. Assume the shocks are Gaussian; then the whole vector (x_1, \dots, x_T) is multivariate normal, and we pick the parameters that make the observed series most probable. The trick that makes this computable is the prediction-error (innovations) decomposition: the joint density factorises into a product of one-step-ahead conditional densities,

\ell(\phi,\theta,\sigma^2) = -\tfrac{1}{2}\sum_{t=1}^{T}\left[\log(2\pi v_t) + \frac{e_t^2}{v_t}\right],

where e_t is the one-step prediction error and v_t its variance. Those e_t, v_t are produced, recursively and exactly, by running the Kalman filter over the series — which is why in practice "fit an ARMA by MLE" and "run a Kalman filter inside an optimiser" mean the same thing. Maximum likelihood is preferred because it is asymptotically efficient (smallest possible variance among reasonable estimators), it handles MA terms that defeat the moment method, and it comes with a ready-made theory of uncertainty.

Standard errors come for free

Maximum likelihood does not just hand you a point estimate; it hands you its precision. The curvature of the log-likelihood at its peak — the Fisher information — sets the asymptotic covariance:

\hat{\boldsymbol\beta} \;\approx\; \mathcal N\!\big(\boldsymbol\beta,\; \mathcal I(\boldsymbol\beta)^{-1}\big), \qquad \mathcal I = -\,\mathbb E\!\left[\frac{\partial^2 \ell}{\partial \boldsymbol\beta\,\partial \boldsymbol\beta^{\!\top}}\right].

A sharply-peaked likelihood (high curvature) means a tightly-determined parameter and a small standard error; a flat ridge means the data barely pin the parameter down. The reported standard errors, and hence the "is this coefficient significant?" verdicts, are square roots of the diagonal of that inverse information matrix. This is the same machinery that gives you AIC and BIC, which reuse the maximised log-likelihood to compare whole models.

For a pure AR the sum-of-squares surface is convex (that clean parabola), so any optimiser marches straight to the unique answer. The instant you add a moving-average term, that guarantee evaporates: the likelihood as a function of \theta can be multimodal, riddled with local optima and flat plateaus, and near the invertibility boundary (|\theta| \to 1) it flattens into a ridge where the optimiser stalls. Fire a naïve gradient search from a random start and it may converge to a decoy peak — a plausible-looking but wrong fit. The fix is to seed the optimiser well: use the fast Yule–Walker or Hannan–Rissanen estimates as starting values, and be suspicious of a fit whose coefficients sit right on the |\theta| = 1 edge. Never trust a single run of the optimiser without checking it landed somewhere sensible.

For a long, clean AR series the two barely differ — conditional least squares and MLE agree to the third decimal, and least squares is faster. The gap opens in three places. First, short samples: exact likelihood uses the information in the first few observations that conditional least squares throws away, which matters when T is small or \phi is near one. Second, moving-average terms: least squares has to reconstruct unobserved shocks and loses its closed form, while the likelihood handles them natively through the innovations recursion. Third, efficiency and inference: only the likelihood delivers the asymptotically minimum-variance estimate and the Fisher-information standard errors in one shot. Least squares is the sketch; maximum likelihood is the finished drawing.