Maximum Likelihood Estimation for ARMA Models (Gaussian MLE)

1) The goal: “best-fitting” ARMA coefficients for a chosen order $(p,q)$

Suppose you have observed a time series dataset $x_1,\dots,x_n$, and you have decided to try an ARMA(p,q) model as a useful approximation. The model has:

• AR coefficients: $\phi_1,\dots,\phi_p$
• MA coefficients: $\theta_1,\dots,\theta_q$
• Noise variance: $\sigma^2$ (variance of the white-noise shocks)

For a fixed order $(p,q)$, you want to pick coefficient values that make the observed data “most plausible” under that model. The approach here defines “most plausible” using maximum likelihood.

In practice, software (e.g., R’s arima) does two connected things:

1. Produces preliminary coefficient estimates (for example, Yule–Walker for AR parts).
2. Uses those as a starting point for a numerical optimization that maximizes a likelihood function.

So preliminary estimation and maximum likelihood are not competing ideas—they are complementary parts of the same workflow.

---

2) Maximum likelihood in a simple, familiar setting (IID Gaussian)

This section builds intuition for what “likelihood” means before moving to time series.

Setup

Assume $X_1,\dots,X_n$​ are IID Normal with mean 0 and unknown variance $\sigma^2$:

$$X_i \sim \mathcal{N}(0,\sigma^2)$$

The joint density for observing values $x_1,\dots,x_n$​ is:

$$f(x_1,\dots,x_n \mid \sigma^2)= \frac{1}{(2\pi\sigma^2)^{n/2}} \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^{n} x_i^2\right)$$

Likelihood idea

• In probability, $f(\cdot\mid\sigma^2)$ is seen as a function of the data given a parameter.
• In likelihood, we treat the observed data as fixed and see the same expression as a function of $\sigma^2$:

$$L(\sigma^2 \mid x_1,\dots,x_n) = f(x_1,\dots,x_n \mid \sigma^2)$$

The maximum likelihood estimate $\hat{\sigma}^2$ is the value of $\sigma^2$ that maximizes $L(\sigma^2)$.

Why use the log-likelihood

Because $\log(\cdot)$ is increasing, maximizing $L$ is equivalent to maximizing $\log L$, and logs simplify products into sums.

Log-likelihood:$$\log L(\sigma^2)= -\frac{n}{2}\log(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum_{i=1}^{n}x_i^2$$

Differentiate with respect to $\sigma^2$, set to 0, solve, and you get:

$$\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}x_i^2$$

This gives a clear interpretation of MLE:

• you choose the parameter values that make the observed data most likely under the assumed probabilistic model.

---

3) Extending likelihood to a Gaussian time series (multivariate normal)

A time series is not independent across time. That dependence is captured by a covariance structure.

Assume the observed vector$$X_n = (X_1,\dots,X_n)^\top$$

is multivariate normal with mean 0 and covariance matrix $\Gamma_n$:

$$X_n \sim \mathcal{N}(0,\Gamma_n)$$

Then the Gaussian likelihood for $\Gamma_n$​ is:

$$L(\Gamma_n)=\frac{1}{(2\pi)^{n/2}(\det \Gamma_n)^{1/2}} \exp\left(-\frac{1}{2}X_n^\top \Gamma_n^{-1} X_n\right)$$

This looks intimidating because it involves:

• $\Gamma_n^{-1}$​ (matrix inverse)
• $\det(\Gamma_n)$ (determinant)

Both are expensive or messy to compute directly when $n$ is large.

---

4) The key simplification: rewrite the likelihood using one-step prediction errors

A Gaussian time series has a very important property:

If you predict each $X_{j}$ using the best linear predictor based on past observations $X_1,\dots,X_{j-1}$​, then the prediction errors behave in a particularly convenient way.

Define:

• $\hat{X}_j = P(X_j \mid X_1,\dots,X_{j-1})$ (best one-step predictor)
• prediction error: $e_j = X_j – \hat{X}_j$
• prediction error variance: $v_{j-1} = E(e_j^2)$

A nontrivial linear algebra result shows:

$$X_n^\top \Gamma_n^{-1} X_n = \sum_{j=1}^{n}\frac{(X_j-\hat{X}_j)^2}{v_{j-1}}$$

and$$\det(\Gamma_n)= v_0 v_1 \cdots v_{n-1}$$

Substituting these into the likelihood gives a simpler expression:

$$L(\Gamma_n)=\frac{1}{(2\pi)^{n/2}\sqrt{v_0v_1\cdots v_{n-1}}} \exp\left(-\frac{1}{2}\sum_{j=1}^{n}\frac{(X_j-\hat{X}_j)^2}{v_{j-1}}\right)$$

Conceptually, this says:

• the likelihood can be expressed using standardized one-step-ahead errors and their variances
• instead of directly using $\Gamma_n^{-1}$ and $\det(\Gamma_n)$

---

5) Specializing to ARMA(p,q): parameters are $\phi,\theta,\sigma^2$

For a causal ARMA(p,q) model, the parameters of interest are:

$$\phi=(\phi_1,\dots,\phi_p),\quad \theta=(\theta_1,\dots,\theta_q),\quad \sigma^2$$

Here $\sigma^2$ is the variance of the white noise shocks $W_t$.

A crucial normalization is introduced:

$$r_j := \frac{v_j}{\sigma^2} \quad\text{for } j=0,\dots,n-1$$

This is useful because:

• $r_j$ does not depend on $\sigma^2$ (it depends on $\phi,\theta$)
• the predictors $\hat{X}_j$​ also do not depend on $\sigma^2$

So $\sigma^2$ becomes easier to “separate out” analytically.

Define the function:

$$S(\phi,\theta) := \sum_{j=1}^{n}\frac{(X_j-\hat{X}_j)^2}{r_{j-1}}$$

Then the Gaussian likelihood for an ARMA model can be written as:

$$L(\phi,\theta,\sigma^2)= \frac{1}{(2\pi\sigma^2)^{n/2}\sqrt{r_0r_1\cdots r_{n-1}}} \exp\left(-\frac{1}{2\sigma^2}S(\phi,\theta)\right)$$

This is the central likelihood expression for estimation.

---

6) Why numerical search is needed for $\phi$ and $\theta$

Even though the likelihood formula looks explicit, the difficulty is that:

• $\hat{X}_j$​ depends on $\phi,\theta$ in a complicated way
• each $r_j$ depends on $\phi,\theta$ in a complicated way

Therefore the likelihood surface over $(\phi,\theta)$ generally has no closed-form solution. You must use numerical optimization.

However, $\sigma^2$ can be handled analytically.

---

7) Eliminating $\sigma^2$: closed-form MLE for the noise variance once $\phi,\theta$ are known

Take logs:

$$\log L(\phi,\theta,\sigma^2)= -\frac{n}{2}\ln(2\pi\sigma^2) -\frac{1}{2}\sum_{j=0}^{n-1}\ln r_j -\frac{1}{2\sigma^2}S(\phi,\theta)$$

Differentiate with respect to $\sigma^2$, set to zero at the optimum, solve:

$$\hat{\sigma}^2 = \frac{1}{n}S(\hat{\phi},\hat{\theta})$$

This is directly analogous to the IID Gaussian case:

• “variance estimate = average standardized squared errors”

Then you plug this back into the log-likelihood, which yields an objective function depending only on $(\phi,\theta)$. In other words:

• First, find $(\hat{\phi},\hat{\theta})$ by minimizing a function $\ell(\phi,\theta)$ derived from the concentrated likelihood.
• Then compute $\hat{\sigma}^2$ from the formula above.

This reduces the optimization dimension and makes computation more stable.

---

8) What R’s arima() is doing in practice

When you run arima(x, order=c(p,d,q), ...), it typically:

1. Uses a preliminary procedure to get starting values (for AR components, Yule–Walker is a common choice).
2. Runs a numerical optimizer to maximize the Gaussian likelihood (equivalently minimize the negative log-likelihood).
3. Reports:
   • coefficient estimates (AR and MA terms, plus an intercept if included)
   • estimated noise variance $\hat{\sigma}^2$
   • standard errors of coefficient estimates
   • maximized log-likelihood value
   • AIC (used later for comparing models across different orders)

---

9) Example: AR(2) simulation — comparison of preliminary estimates and MLE

A simulated AR(2):$$X_t = -0.3X_{t-1} + 0.4X_{t-2} + W_t,\quad \text{Var}(W_t)=1$$

• Preliminary (Yule–Walker) produced approximately:
  • $\hat{\phi}_1 \approx -0.287$
  • $\hat{\phi}_2 \approx 0.407$
  • $\hat{\sigma}^2 \approx 1.024$

Running:

• arima(X, order=c(2,0,0), include.mean=FALSE)

produces nearly the same coefficient values:

• ar1 $\approx -0.2866$, ar2 $\approx 0.4074$
• $\hat{\sigma}^2 \approx 1.024$

Interpretation: with a large sample size, Yule–Walker and Gaussian MLE can be extremely close for AR models, because both are using essentially the same correlation information, and the likelihood is sharply peaked near the true parameter.

Standard errors

The output includes s.e. for each coefficient. These are estimates of how much the coefficient estimate would vary across repeated samples from the same process.

A rough 95% confidence interval is:

$$\hat{\phi}_j \pm 1.96\cdot \text{s.e.}(\hat{\phi}_j)$$

In the example, the interval width is small, reflecting the large $n$.

---

10) Example: LakeHuron — why an intercept appears automatically

Running:

• arima(LakeHuron, order=c(2,0,0))

returns AR coefficients and also an intercept (a constant level term).

This happens because:

• by default, the model includes a mean level unless you explicitly exclude it
• real-world series often have a nonzero mean

The model reported can be read as an AR(2) model around a mean level:

$$Y_t – \mu = \phi_1(Y_{t-1}-\mu) + \phi_2(Y_{t-2}-\mu) + W_t$$

and the output gives estimates of $\phi_1,\phi_2,\mu$ and $\sigma^2$.

Numerically, the AR coefficients and variance estimates are close to the earlier method-of-moments fit:

• both methods suggest similar dynamics
• MLE is typically preferred as the final fit because it directly optimizes the probability of the observed data under the model

---

11) What to remember conceptually

1. Likelihood is a way to score parameter values by how plausible they make the observed dataset under the model.
2. In Gaussian time series, the likelihood for $X_1,\dots,X_n$​ depends on the covariance structure; for ARMA, that covariance structure is controlled by $\phi,\theta,\sigma^2$.
3. A major simplification rewrites the likelihood in terms of one-step prediction errors and their variances.
4. $\sigma^2$ can be estimated explicitly once $\phi,\theta$ are chosen: $$\hat{\sigma}^2 = \frac{1}{n}S(\hat{\phi},\hat{\theta})$$
5. $\phi,\theta$ usually require numerical optimization, which is why software is essential.
6. Preliminary estimates (like Yule–Walker for AR terms) are often used to initialize the numerical search; MLE then refines.