Autoregressive and moving-average models (AR, MA, ARMA)

Before this lesson: lesson 02 (stationarity, white noise, ACF, PACF, the AR(1) process) and lesson 03 (forecasts and residuals). From other tracks: eigenvalues (linear-algebra lesson 05), least squares (linear-algebra lesson 12), Toeplitz matrices (linear-algebra lesson 19) and maximum likelihood (the probability-and-statistics track). Time: about 90 minutes, plus the notebook.

---

Why this matters #

• ARIMA is built from these parts. Lesson 06 adds differencing and seasonality to ARMA and gives forecast intervals. Almost every forecasting paper uses ARIMA as a benchmark.
• Lag features. A regression on $y_{t-1}, \dots, y_{t-p}$ is an AR($p$) model. Gradient-boosted forecasters (lesson 08) use the same lags as input features.
• Linear recurrences in deep learning. A state-space layer (S4, Mamba) is a linear recurrence $x_t = A x_{t-1} + B u_t$. It is stable exactly when the eigenvalues of $A$ lie inside the unit circle, the same condition as for an AR model in companion form (deep-learning lesson 16).
• Science and engineering. Speech codecs model each short frame of speech as an AR process (linear predictive coding). Climate scientists use an AR(1) process ("red noise") as the null model for climate variability. Yule (1927) invented the AR(2) model to explain the sunspot cycle.

1. Predict first #

2. The picture #

Start with AR(1) from lesson 02: $y_t = \phi\, y_{t-1} + \varepsilon_t$ with white noise $\varepsilon_t \sim \mathrm{WN}(0, \sigma^2)$. Unroll it: $y_t = \varepsilon_t + \phi\varepsilon_{t-1} + \phi^2\varepsilon_{t-2} + \dots$. A shock that happened $j$ steps ago still has weight $\phi^j$. If $|\phi| < 1$ the weights shrink and the series is stationary; its ACF is $\rho(k) = \phi^k$. If $|\phi| \ge 1$ old shocks never fade.

AR(2): two lags of feedback. $y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \varepsilon_t$. Now the pair $(\phi_1, \phi_2)$ decides everything. In the figure, each point of the plane is one AR(2) model.

Interactive figure: loading…

Three regions matter:

• Inside the triangle, above the parabola: two real roots. The ACF decays like a mix of two geometric sequences. The path wanders smoothly ($\phi_1 > 0$) or flips sign ($\phi_1 < 0$).
• Inside the triangle, below the parabola: complex roots. The ACF is a damped cosine. The path shows a quasi-cycle: it oscillates with a typical period, but the amplitude and phase change at random. This is the sunspot picture.
• Outside the triangle: no stationary (causal) solution. Shocks are amplified and the path explodes, even though each coefficient may be smaller than 1.

Roots and eigenvalues. The stationarity rule has two equivalent forms. Write the AR polynomial $\phi(z) = 1 - \phi_1 z - \phi_2 z^2$. The rule "all roots of $\phi(z)$ lie outside the unit circle" is the same as "all eigenvalues of the companion matrix $\begin{pmatrix} \phi_1 & \phi_2 \\ 1 & 0 \end{pmatrix}$ lie inside the unit circle", because the eigenvalues are exactly the reciprocals $1/z$ of the roots. Section 3 proves it. The next figure lets you place the roots yourself.

Interactive figure: loading…

MA(1): a shock that lives for two steps. $y_t = \varepsilon_t + \theta\varepsilon_{t-1}$. The shock $\varepsilon_t$ enters $y_t$ and $y_{t+1}$ and then disappears. So $y_t$ and $y_{t+2}$ share no shock, and $\rho(2) = 0$. The ACF of an MA(1) cuts off after lag 1. Its PACF, however, decays slowly (section 3 explains why).

Interactive figure: loading…

3. Definitions and main results #

The backshift operator #

The backshift operator $B$ moves a series one step back: $B y_t = y_{t-1}$, and $B^k y_t = y_{t-k}$. Polynomials in $B$ act term by term: $(1 - 0.5B)\,y_t = y_t - 0.5\,y_{t-1}$. They multiply like ordinary polynomials: $(1 - aB)(1 - bB) = 1 - (a + b)B + ab\,B^2$. This turns a model into an algebra problem in one variable.

AR($p$) #

An autoregressive model of order $p$, AR($p$), is

or in short form $\phi(B)\,y_t = c + \varepsilon_t$ with the AR polynomial $\phi(z) = 1 - \phi_1 z - \dots - \phi_p z^p$. The shock $\varepsilon_t$ is also called the innovation: it is the part of $y_t$ that the past cannot predict.

Mean. If the process is stationary, $\mathrm E\,y_t = \mu$ for every $t$. Take expectations of both sides: $\mu = c + \phi_1\mu + \dots + \phi_p\mu$. So

The constant $c$ is not the mean. With $c = 2$ and $\phi_1 = 0.8$, the mean is $2/0.2 = 10$. Equivalently, $y_t - \mu = \phi_1(y_{t-1} - \mu) + \dots + \phi_p(y_{t-p} - \mu) + \varepsilon_t$.

Companion form. Stack the last $p$ values into a state vector $x_t = (y_t, y_{t-1}, \dots, y_{t-p+1})^\top$. Then the AR($p$) becomes a first-order vector recursion

with $e_1 = (1, 0, \dots, 0)^\top$. The first row is the AR equation; the other rows only shift the old values down by one place. $A$ is the companion matrix of the AR model.

Theorem (stationarity of AR($p$)). Assume $\phi_p \ne 0$. The following are equivalent:

1. all roots of $\phi(z)$ lie outside the unit circle, $|z| > 1$;
2. all eigenvalues of the companion matrix $A$ lie inside the unit circle, $|\lambda| < 1$.

Under these conditions the AR equation has exactly one stationary solution that depends only on present and past shocks (a causal solution): $y_t = \mu + \sum_{j \ge 0} \psi_j \varepsilon_{t-j}$ with $\psi_j = (A^j)_{11}$, and $\sum_j |\psi_j| < \infty$. If some eigenvalue has $|\lambda| \ge 1$, no causal stationary solution exists.

Proof of 1 ⇔ 2.

1. The characteristic polynomial of the companion matrix is $\det(\lambda I - A) = \lambda^p - \phi_1\lambda^{p-1} - \dots - \phi_{p-1}\lambda - \phi_p$. (Expand the determinant along the first row; for $p = 2$: $\det\begin{pmatrix} \lambda - \phi_1 & -\phi_2 \\ -1 & \lambda \end{pmatrix} = \lambda^2 - \phi_1\lambda - \phi_2$.)
2. Multiply $\phi(1/\lambda)$ by $\lambda^p$: $\lambda^p\,\phi(1/\lambda) = \lambda^p - \phi_1\lambda^{p-1} - \dots - \phi_p$. This is the same polynomial.
3. $\phi_p \ne 0$, so $\lambda = 0$ is not an eigenvalue. Therefore $\lambda$ is an eigenvalue of $A$ exactly when $z = 1/\lambda$ is a root of $\phi(z)$.
4. $|z| = 1/|\lambda|$, so $|z| > 1 \iff |\lambda| < 1$. $\blacksquare$

Why the eigenvalues decide stationarity

Eigenvalues inside the circle give a stationary solution.

1. Unroll the companion recursion $k$ times: $x_t = A^k x_{t-k} + \sum_{j=0}^{k-1} A^j e_1 (c + \varepsilon_{t-j})$.
2. If every eigenvalue has $|\lambda| < 1$, then $\|A^j\| \le C r^j$ for some $C$ and some $r < 1$ (by the Jordan form of linear-algebra lesson 05, the entries of $A^j$ are sums of terms $j^m\lambda^j$, and $j^m|\lambda|^j \to 0$ when $|\lambda| < 1$). So $A^k x_{t-k} \to 0$, and the sum converges absolutely when $k \to \infty$.
3. The limit $x_t = \sum_{j \ge 0} A^j e_1 (c + \varepsilon_{t-j})$ uses only present and past shocks and has the same distribution for every $t$: it is a causal stationary solution. Its first entry is $y_t = \mu + \sum_j \psi_j\varepsilon_{t-j}$ with $\psi_j = (A^j)_{11}$.

An eigenvalue on or outside the circle rules it out.

1. Let $w^\top A = \lambda w^\top$ with $|\lambda| \ge 1$ (a left eigenvector). For a companion matrix the first entry $w_1$ of a left eigenvector is never 0 (if $w_1 = 0$, the equations $w^\top A = \lambda w^\top$ force all other entries to 0, column by column).
2. Put $u_t = w^\top x_t$. Then $u_t = \lambda u_{t-1} + w_1(c + \varepsilon_t)$.
3. In a causal stationary solution, $u_{t-1}$ depends only on shocks up to $t - 1$, so it is uncorrelated with $\varepsilon_t$, and $\operatorname{Var}(u_t) = \operatorname{Var}(u_{t-1}) = v$. Then $v = |\lambda|^2 v + |w_1|^2\sigma^2 > |\lambda|^2 v \ge v$. This is impossible. $\blacksquare$

AR(2): the triangle and the quasi-cycle #

For AR(2) the companion eigenvalues solve $\lambda^2 - \phi_1\lambda - \phi_2 = 0$, so $\lambda_1 + \lambda_2 = \phi_1$ and $\lambda_1\lambda_2 = -\phi_2$.

Theorem (AR(2) stationarity triangle). An AR(2) is stationary if and only if

Its roots are complex when $\phi_1^2 + 4\phi_2 < 0$ (below the parabola).

Proof of the triangle

Let $g(\lambda) = \lambda^2 - \phi_1\lambda - \phi_2$. Note $g(1) = 1 - \phi_1 - \phi_2$ and $g(-1) = 1 + \phi_1 - \phi_2$.

Real roots ($\phi_1^2 + 4\phi_2 \ge 0$).

1. If both roots lie in $(-1, 1)$, then $g(1) > 0$ and $g(-1) > 0$ (an upward parabola is positive outside its roots), and $|\phi_2| = |\lambda_1\lambda_2| < 1$.
2. Conversely, let $g(1) > 0$, $g(-1) > 0$ and $|\lambda_1\lambda_2| < 1$. A root in $(-\infty, -1]$ and another in $[1, \infty)$ would make $g(1) \le 0$. So both roots are on the same side of $1$ and on the same side of $-1$: both in $(-1, 1)$, both $> 1$, or both $< -1$. The last two cases give $|\lambda_1\lambda_2| > 1$. So both roots lie in $(-1, 1)$.

Complex roots ($\phi_1^2 + 4\phi_2 < 0$). Then $\lambda_2 = \bar\lambda_1$ and $|\lambda_1|^2 = \lambda_1\bar\lambda_1 = -\phi_2$. So $|\lambda| < 1 \iff \phi_2 > -1$. The other two conditions hold automatically: $\phi_2 < -\phi_1^2/4$ gives $\phi_1 + \phi_2 < \phi_1 - \phi_1^2/4 = 1 - (1 - \phi_1/2)^2 \le 1$, and the same for $-\phi_1$. $\blacksquare$

Quasi-cycle. With complex roots, write the eigenvalues as $\lambda = r e^{\pm i\omega}$. From $|\lambda|^2 = -\phi_2$ and $\lambda_1 + \lambda_2 = 2r\cos\omega = \phi_1$:

The ACF solves the recursion $\rho(k) = \phi_1\rho(k-1) + \phi_2\rho(k-2)$, so it is a combination of $\lambda_1^k$ and $\lambda_2^k$. With complex eigenvalues this combination is real and equals $a\,r^k\cos(k\omega + b)$ for some constants $a, b$: a cosine with period $2\pi/\omega$ that decays by the factor $r$ per step. Random shocks keep restarting the oscillation, so the path never settles.

Yule–Walker equations and the ACF of AR($p$) #

Theorem (Yule–Walker). For a stationary AR($p$) with autocovariance $\gamma(k)$ and ACF $\rho(k) = \gamma(k)/\gamma(0)$:

Proof.

1. Subtract the mean: $\tilde y_t = y_t - \mu$ satisfies $\tilde y_t = \phi_1\tilde y_{t-1} + \dots + \phi_p\tilde y_{t-p} + \varepsilon_t$.
2. Multiply both sides by $\tilde y_{t-k}$ and take expectations. On the left: $\mathrm E[\tilde y_t\tilde y_{t-k}] = \gamma(k)$. On the right: $\phi_i\,\gamma(k - i)$ for each lag.
3. The last term is $\mathrm E[\varepsilon_t\tilde y_{t-k}]$. For $k \ge 1$ it is 0, because in the causal solution $\tilde y_{t-k}$ depends only on shocks up to $t - k$. For $k = 0$ it is $\mathrm E[\varepsilon_t\tilde y_t] = \sigma^2$ (only the $\psi_0\varepsilon_t$ term of $\tilde y_t$ contributes).
4. So $\gamma(k) = \sum_i\phi_i\gamma(k - i)$ for $k \ge 1$ and $\gamma(0) = \sum_i\phi_i\gamma(i) + \sigma^2$. Divide by $\gamma(0)$ and use $\gamma(-i) = \gamma(i)$. $\blacksquare$

The equations for $k = 1, \dots, p$ form a linear system with a symmetric Toeplitz matrix (constant along each diagonal, see linear-algebra lesson 19):

Two uses. Forward: given $\phi$, solve for $\rho(1), \dots, \rho(p)$ and then run the recursion to get the whole ACF. Backward (estimation): replace $\rho(k)$ by the sample ACF $\hat\rho(k)$ and solve for $\hat\phi$. This is the Yule–Walker estimator. With the sample ACF that divides by $n$ (lesson 02), $R_p$ is positive definite and the fitted model is always stationary.

The ACF of an AR($p$) decays. For $k > p$ the recursion has no forcing term, so $\rho(k)$ is a combination of $\lambda_i^k$ over the companion eigenvalues. All $|\lambda_i| < 1$, so $\rho(k)$ decays geometrically: like $\phi^k$, like a sum of such terms, or like a damped cosine. It never cuts off.

The PACF of an AR($p$) cuts off after lag $p$. Recall from lesson 02 that the PACF value $\phi_{kk}$ is the last coefficient of the best linear predictor of $y_t$ from $y_{t-1}, \dots, y_{t-k}$. For $k \ge p$, the best predictor is the AR equation itself, because the remaining error $\varepsilon_t$ is uncorrelated with all past values. Its coefficient vector is $(\phi_1, \dots, \phi_p, 0, \dots, 0)$. So $\phi_{pp} = \phi_p$ and $\phi_{kk} = 0$ for every $k > p$.

MA($q$) #

A moving-average model of order $q$, MA($q$), is

with the MA polynomial $\theta(z) = 1 + \theta_1 z + \dots + \theta_q z^q$. It is a weighted sum of the last $q + 1$ shocks; it is not a moving average of the data (that is a smoother, lesson 01).

An MA($q$) is always stationary. It is a finite sum of white noise terms, so its mean $c$ and its autocovariances do not depend on $t$. With $\theta_0 = 1$:

Proof. $\operatorname{Cov}(y_t, y_{t+k}) = \sum_{i}\sum_{j}\theta_i\theta_j\operatorname{Cov}(\varepsilon_{t-i}, \varepsilon_{t+k-j})$. The covariance is $\sigma^2$ when $t - i = t + k - j$, that is $j = i + k$, and 0 otherwise. If $k > q$, no pair $(i, i + k)$ has both indices in $\{0, \dots, q\}$. $\blacksquare$

So the ACF of an MA($q$) cuts off after lag $q$. For MA(1):

Since $(1 - |\theta|)^2 \ge 0$ gives $1 + \theta^2 \ge 2|\theta|$, we get $|\rho(1)| \le 1/2$, with equality only at $\theta = \pm 1$. An MA(1) can never have a lag-1 autocorrelation above 0.5.

Invertibility. Replace $\theta$ by $1/\theta$ and $\sigma^2$ by $\theta^2\sigma^2$: $\gamma(0) = \theta^2\sigma^2(1 + 1/\theta^2) = \sigma^2(1 + \theta^2)$ and $\gamma(1) = \theta^2\sigma^2 \cdot (1/\theta) = \sigma^2\theta$. Same autocovariances. With Gaussian shocks the two models give exactly the same distribution of the data, so no data set can tell them apart. We pick one by a rule:

An MA (or ARMA) model is invertible if all roots of $\theta(z)$ lie outside the unit circle. For MA(1) this means $|\theta| < 1$. Then the shocks can be recovered from past data. From $\varepsilon_t = y_t - \theta\varepsilon_{t-1}$ (take $c = 0$) and repeated substitution:

The sum converges because $|\theta| < 1$. So an invertible MA(1) is an AR($\infty$): $y_t = \theta y_{t-1} - \theta^2 y_{t-2} + \theta^3 y_{t-3} - \dots + \varepsilon_t$. This has two consequences. First, residuals computed from data converge to the true shocks, which forecasting needs (lesson 06). Second, the best linear predictor uses all past lags, so the PACF of an MA(1) never cuts off: it decays like $(-\theta)^k$.

ARMA($p, q$) #

An ARMA($p, q$) model combines both parts:

• It is causal (stationary) if the roots of $\phi(z)$ lie outside the unit circle. Then $y_t = \mu + \sum_{j\ge0}\psi_j\varepsilon_{t-j}$, the MA($\infty$) form, with $\psi(z) = \theta(z)/\phi(z)$. Comparing coefficients in $\phi(z)\psi(z) = \theta(z)$ gives the $\psi$-weights by a recursion: $\psi_0 = 1$ and $\psi_j = \theta_j + \sum_{i=1}^{\min(j, p)}\phi_i\psi_{j-i}$, with $\theta_j = 0$ for $j > q$. For ARMA(1,1): $\psi_j = \phi^{j-1}(\phi + \theta)$ for $j \ge 1$.
• It is invertible if the roots of $\theta(z)$ lie outside the unit circle. Then $\varepsilon_t = \pi(B)(y_t - \mu)$ with $\pi(z) = \phi(z)/\theta(z)$, the AR($\infty$) form.
• Parameter redundancy. If $\phi(z)$ and $\theta(z)$ share a root, it cancels. Example: $y_t = 0.5y_{t-1} + \varepsilon_t - 0.5\varepsilon_{t-1}$ is $(1 - 0.5B)y_t = (1 - 0.5B)\varepsilon_t$, which is white noise. Fitting software can return such a pair with large, opposite coefficients; it is a sign that the model is too big.

Both the ACF and the PACF of an ARMA($p, q$) with $p, q \ge 1$ decay and do not cut off. For ARMA(1,1) the ACF is $\rho(1) = \dfrac{(1 + \phi\theta)(\phi + \theta)}{1 + 2\phi\theta + \theta^2}$ and $\rho(k) = \phi\,\rho(k-1)$ for $k \ge 2$.

Wold's theorem: why ARMA is general

Wold's decomposition theorem. Every zero-mean, covariance-stationary process can be written as $y_t = \sum_{j\ge0}\psi_j\varepsilon_{t-j} + d_t$, with $\psi_0 = 1$, $\sum_j\psi_j^2 < \infty$, $\varepsilon_t$ the white-noise one-step prediction errors (the innovations), and $d_t$ a deterministic part that the past predicts perfectly (for example a pure sinusoid). So every stationary process without a deterministic part is an MA($\infty$). An ARMA model approximates the infinite sequence $\psi_j$ by a ratio of two short polynomials $\theta(z)/\phi(z)$. Few parameters can then describe a long memory pattern, which is why ARMA models work well in practice.

Estimation #

• AR($p$) by conditional least squares. For $t = p + 1, \dots, n$, regress $y_t$ on $1, y_{t-1}, \dots, y_{t-p}$. This is ordinary least squares (linear-algebra lesson 12); "conditional" means we condition on the first $p$ values. Estimate $\sigma^2$ by the mean squared residual. The Yule–Walker estimator gives almost the same answer on long series.
• MA and ARMA by maximum likelihood. The residuals $e_t = y_t - c - \sum_i\phi_i y_{t-i} - \sum_j\theta_j e_{t-j}$ depend on $\theta$ through the recursion, so the problem is not linear. Software maximises the exact Gaussian likelihood numerically. The likelihood is computed with the Kalman filter (lesson 07). A cheaper version minimises the conditional sum of squares $\sum_t e_t^2$ with $e_t = 0$ before the start.

Choosing the order #

Model	ACF	PACF
AR($p$)	decays (geometric or damped cosine)	cuts off after lag $p$
MA($q$)	cuts off after lag $q$	decays
ARMA($p, q$)	decays (after lag $q$)	decays (after lag $p$)

A sample ACF or PACF value counts as "zero" if it lies inside the band $\pm 1.96/\sqrt n$ (lesson 02). The table is a good start on long, clean series. On short or noisy series the patterns are hard to see. Practice the reading in the next figure.

Interactive figure: loading…

Therefore, in practice, fit a small grid of candidate models and compare them with an information criterion. With $\hat L$ the maximised likelihood and $k$ the number of estimated parameters (including $c$ and $\sigma^2$):

The first term rewards fit; the second penalises the number of parameters, so the criterion trades fit against complexity. Smaller is better. AICc corrects AIC for small samples; BIC penalises more and picks smaller models. Compare criteria only between models fitted to the same data: the same time points and the same transformation and differencing.

Then check the residuals. If the model is adequate, its residuals look like white noise: no spikes in their ACF. The Ljung–Box test summarises the first $h$ residual autocorrelations:

If the residuals are white noise, $Q$ approximately follows a $\chi^2$ distribution with $h - p - q$ degrees of freedom. A small p-value says that correlation is left in the residuals.

Answers to section 1. (1) It explodes. $\phi_1 + \phi_2 = 1.1 > 1$ breaks the first triangle condition. Equivalently, $g(\lambda) = \lambda^2 - 0.5\lambda - 0.6$ has roots $\lambda \approx 1.06$ and $-0.56$; one eigenvalue is outside the unit circle. Small coefficients do not guarantee stability. (2) At most $0.5$, reached at $\theta = 1$. (3) Yes. With complex roots an AR(2) produces a quasi-cycle with period $2\pi/\omega$. Section 4 fits it to the sunspots and finds a period of about 10.9 years.

4. Worked example #

(a) An AR(2) by hand. Take $y_t = y_{t-1} - 0.5\,y_{t-2} + \varepsilon_t$, so $\phi_1 = 1$, $\phi_2 = -0.5$.

1. Triangle: $\phi_1 + \phi_2 = 0.5 < 1$, $\phi_2 - \phi_1 = -1.5 < 1$, $|\phi_2| = 0.5 < 1$. Stationary.
2. Roots of $\phi(z) = 1 - z + 0.5z^2$: $z = \dfrac{1 \pm \sqrt{1 - 2}}{1} = 1 \pm i$, with modulus $\sqrt 2 \approx 1.414 > 1$.
3. Companion eigenvalues: $\lambda^2 - \lambda + 0.5 = 0$ gives $\lambda = (1 \pm i)/2$, with modulus $\sqrt{0.5} \approx 0.707 = 1/\sqrt 2$. They are the reciprocals of the roots: $1/(1 + i) = (1 - i)/2$.
4. Period: $\cos\omega = \phi_1/(2\sqrt{-\phi_2}) = 1/(2 \cdot 0.7071) = 0.7071$, so $\omega = \pi/4$ and the period is $2\pi/(\pi/4) = 8$ steps.
5. ACF by Yule–Walker. $\rho(1) = \phi_1/(1 - \phi_2) = 1/1.5 = 0.667$. Then $\rho(k) = \rho(k-1) - 0.5\rho(k-2)$. For example $\rho(2) = 0.667 - 0.5 \cdot 1 = 0.167$ and $\rho(4) = -0.167 - 0.5 \cdot 0.167 = -0.250$. The table lists the first values.

$k$	0	1	2	3	4	5	6	7	8
$\rho(k)$	1	0.667	0.167	−0.167	−0.250	−0.167	−0.042	0.042	0.063

The sign pattern repeats every 8 lags and the size shrinks by the factor $r = 0.707$ per lag: a damped cosine.

6. PACF: $\phi_{11} = \rho(1) = 0.667$, $\phi_{22} = \phi_2 = -0.5$, and $\phi_{kk} = 0$ for $k \ge 3$.
7. Variance: $\gamma(0) = \sigma^2/(1 - \phi_1\rho(1) - \phi_2\rho(2)) = \sigma^2/(1 - 0.667 + 0.083) = 2.4\,\sigma^2$.

(b) Two MA(1) models with the same ACF. $\theta = 0.5$ gives $\rho(1) = 0.5/1.25 = 0.4$. $\theta = 2$ gives $\rho(1) = 2/5 = 0.4$. With $\sigma^2 = 1$ for the first and $\sigma^2 = 0.25$ for the second, also $\gamma(0) = 1.25$ and $\gamma(1) = 0.5$ agree. Only $\theta = 0.5$ is invertible: the root of $\theta(z) = 1 + 0.5z$ is $z = -2$, outside the unit circle; the root of $1 + 2z$ is $z = -0.5$, inside. The PACF of the invertible model is $0.4, -0.190, 0.094, -0.047, \dots$: it alternates in sign and decays, roughly by the factor $-\theta$ per lag.

(c) The sunspot cycle with AR(2). The data are the yearly sunspot numbers 1700–2008 ($n = 309$, from the statsmodels data sets). Regress $y_t$ on $(1, y_{t-1}, y_{t-2})$ for the 307 years 1702–2008:

1. Stationary? $1.392 - 0.690 = 0.702 < 1$, $-0.690 - 1.392 < 1$, $|-0.690| < 1$. Yes.
2. Complex? $\phi_1^2 + 4\phi_2 = 1.937 - 2.761 = -0.824 < 0$. Yes: a quasi-cycle.
3. Eigenvalue modulus $r = \sqrt{0.690} = 0.831$; root modulus $1/0.831 = 1.204$.
4. $\cos\omega = 1.392/(2 \cdot 0.831) = 0.838$, so $\omega = 0.578$ and the period is $2\pi/0.578 = 10.9$ years: close to the 11-year solar cycle.
5. Mean: $\hat\mu = 14.91/(1 - 1.392 + 0.690) = 14.91/0.298 = 49.9$; the sample mean is 49.8.

Is AR(2) the best order? On the common years 1710–2008, the AIC is 2732.5 for AR(1), 2536.3 for AR(2), 2533.1 for AR(3) and smallest, 2485.8, for AR(9). The Ljung–Box test on the AR(2) residuals with $h = 10$ gives $Q = 32.4$ with 8 degrees of freedom, p-value 0.0001: some correlation is left. AR(2) captures the cycle with two parameters; AR(9) captures more detail. (The real cycle is also asymmetric, rising faster than it falls, which no linear Gaussian model reproduces.)

Interactive figure: loading…

5. In code #

NumPy checks the hand computations. The notebook embeds the sunspot data; here we load it from statsmodels.

import numpy as np
import statsmodels.api as sm

# (a) AR(2) with phi = (1, -0.5): roots, companion eigenvalues, period, ACF
phi1, phi2 = 1.0, -0.5
roots = np.roots([-phi2, -phi1, 1])            # phi(z) = 1 - phi1 z - phi2 z^2
eig = np.linalg.eigvals([[phi1, phi2], [1, 0]])
print(roots, np.abs(roots))                    # [1.+1.j 1.-1.j] [1.414 1.414]
print(np.abs(eig), 2 * np.pi / np.abs(np.angle(eig[0])))   # [0.707 0.707] 8.0
rho = [1.0, phi1 / (1 - phi2)]
for k in range(2, 9):
    rho.append(phi1 * rho[-1] + phi2 * rho[-2])
print(np.round(rho, 3))   # [ 1. 0.667 0.167 -0.167 -0.25 -0.167 -0.042 0.042 0.063]

# (b) MA(1): theta and 1/theta give the same rho(1)
for th in (0.5, 2.0):
    print(th, th / (1 + th**2))                # 0.4 for both

# (c) AR(2) on the sunspots by least squares
y = sm.datasets.sunspots.load_pandas().data["SUNACTIVITY"].to_numpy()
X = np.column_stack([np.ones(len(y) - 2), y[1:-1], y[:-2]])
b, *_ = np.linalg.lstsq(X, y[2:], rcond=None)
print(np.round(b, 3))                          # [14.907  1.392 -0.69 ]
w = np.arccos(b[1] / (2 * np.sqrt(-b[2])))
print(round(2 * np.pi / w, 2))                 # 10.87

statsmodels.tsa.ar_model.AutoReg(y, lags=2).fit().params returns the same three numbers.

6. Where it shows up #

• Classical forecasting. ARIMA and SARIMA (lesson 06) are ARMA models of a differenced series. Exact likelihood estimation uses the state-space form and the Kalman filter (lesson 07).
• Machine learning. Lag features $y_{t-1}, \dots, y_{t-p}$ in gradient boosting and in linear baselines such as DLinear are nonlinear or multi-output cousins of AR($p$) (lesson 08, lesson 09). AR($p$) residuals are a common anomaly score (lesson 11).
• Deep learning. Linear recurrent layers and state-space models (S4, Mamba) keep the eigenvalues of their transition matrix inside the unit circle for the same reason as the AR stationarity theorem. Autoregressive generation in language models uses the same word: each output depends on the previous outputs.
• Signal processing. Linear predictive coding fits an AR model of order about 10–16 to each speech frame by Yule–Walker (the Levinson–Durbin recursion solves the Toeplitz system in $O(p^2)$ steps). AR spectra give high-resolution spectral estimates.
• Science. AR(1) "red noise" is the null model for climate variability; AR models describe sunspots, EEG rhythms and economic cycles; VAR models describe several linked series (macroeconomics, neuroscience).

Vector autoregression (VAR) and Granger causality

For $K$ series stacked in a vector $\mathbf y_t \in \R^K$, a VAR($p$) model is $\mathbf y_t = \mathbf c + \Phi_1\mathbf y_{t-1} + \dots + \Phi_p\mathbf y_{t-p} + \boldsymbol\varepsilon_t$ with $K \times K$ matrices $\Phi_i$ and a noise covariance $\Sigma$. Each equation is a regression on the lags of all series, so least squares works equation by equation. Stability: the $Kp \times Kp$ block companion matrix (first block row $\Phi_1, \dots, \Phi_p$, identity blocks below) must have all eigenvalues inside the unit circle. The number of parameters, $K^2p$, grows fast; large VARs need regularisation.

Series $x$ Granger-causes series $y$ if the past of $x$ improves the forecast of $y$ beyond the past of $y$ alone; in a VAR it means some lag coefficients from $x$ to $y$ are nonzero (an F-test). This is a statement about prediction, not about causation: a third series that drives both, or different reaction delays, can produce it without any causal link (see the causal-inference track).

7. Common mistakes #

• Checking each coefficient instead of the roots. $|\phi_i| < 1$ for every $i$ does not give stationarity: $(0.5, 0.6)$ explodes. Use the roots, the companion eigenvalues or, for AR(2), the triangle.
• Inside or outside? The roots of $\phi(z)$ must lie outside the unit circle; the companion eigenvalues must lie inside. They are reciprocals, so both statements say the same thing. Software such as statsmodels often reports arroots (outside) and some books plot inverse roots (inside); read the label.
• Reading $c$ as the mean. In $y_t = c + \phi y_{t-1} + \varepsilon_t$ the mean is $c/(1 - \phi)$. Many packages report the mean $\mu$ under the name "const" or "intercept"; check which one you get.
• Mixing MA sign conventions. A fitted $\theta = -0.4$ in statsmodels is $\theta = 0.4$ in Box–Jenkins notation.
• Trusting the ACF/PACF table blindly. On short series the sample ACF is noisy, and about 1 lag in 20 crosses the $\pm 1.96/\sqrt n$ band by chance. Use the table to pick candidates, then compare them by AICc and residual checks.
• Comparing AIC across different data. An AR(2) on years 1702–2008 and an AR(9) on 1709–2008 do not have comparable likelihoods. Fit all candidates on the same time points (and the same differencing).

8. Check yourself #

1. Is $y_t = 1.2\,y_{t-1} - 0.5\,y_{t-2} + \varepsilon_t$ stationary? Are its roots real or complex? If complex, give the period.
   Answer

   Triangle: $1.2 - 0.5 = 0.7 < 1$, $-0.5 - 1.2 = -1.7 < 1$, $|-0.5| < 1$: stationary. $\phi_1^2 + 4\phi_2 = 1.44 - 2 = -0.56 < 0$: complex. $\cos\omega = 1.2/(2\sqrt{0.5}) = 0.849$, $\omega = 0.556$, period $2\pi/0.556 \approx 11.3$ steps.
2. A stationary AR(1) has $c = 3$ and $\phi = 0.7$. What is its mean? What is its variance if $\sigma^2 = 1$?
   Answer

   $\mu = 3/(1 - 0.7) = 10$. $\gamma(0) = \sigma^2/(1 - \phi^2) = 1/0.51 \approx 1.96$.
3. Compute the ACF of the MA(2) $y_t = \varepsilon_t + 0.5\varepsilon_{t-1} + 0.2\varepsilon_{t-2}$.
   Answer

   $\gamma(0)/\sigma^2 = 1 + 0.25 + 0.04 = 1.29$; $\gamma(1)/\sigma^2 = 1 \cdot 0.5 + 0.5 \cdot 0.2 = 0.6$; $\gamma(2)/\sigma^2 = 1 \cdot 0.2 = 0.2$. So $\rho(1) = 0.6/1.29 = 0.465$, $\rho(2) = 0.2/1.29 = 0.155$, and $\rho(k) = 0$ for $k \ge 3$.
4. Is the MA(1) $y_t = \varepsilon_t + 1.25\,\varepsilon_{t-1}$ invertible? Give the invertible model with the same ACF.
   Answer

   No: the root of $1 + 1.25z$ is $z = -0.8$, inside the unit circle. Use $\theta = 1/1.25 = 0.8$ with the noise variance multiplied by $1.25^2 = 1.5625$. Both give $\rho(1) = 1.25/2.5625 = 0.8/1.64 \approx 0.488$.
5. Write the companion matrix of the AR(3) $y_t = 0.5y_{t-1} + 0.2y_{t-2} - 0.1y_{t-3} + \varepsilon_t$. How do you decide stationarity with it?
   Answer

   $A = \begin{pmatrix} 0.5 & 0.2 & -0.1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}$. Compute its eigenvalues (for example with `np.linalg.eigvals`); the model is stationary if all have modulus below 1. Here the characteristic polynomial $\lambda^3 - 0.5\lambda^2 - 0.2\lambda + 0.1 = (\lambda - 0.5)(\lambda^2 - 0.2)$ gives $\lambda = 0.5$ and $\lambda = \pm 0.447$: all inside the unit circle, so the model is stationary.
6. Why does the PACF of an AR(3) vanish at lag 5?
   Answer

   The PACF at lag 5 is the last coefficient of the best linear predictor of $y_t$ from five lags. For an AR(3) that predictor is the AR equation itself, with coefficients $(\phi_1, \phi_2, \phi_3, 0, 0)$, because the error $\varepsilon_t$ is uncorrelated with every past value. So the fifth coefficient is 0.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. When is an AR($p$) model stationary (and causal)?

   Show answer

2. How are the roots of $\phi(z)$ and the companion eigenvalues related?

   Show answer

3. When does an AR(2) have complex roots, and what does the path look like then?

   Show answer

4. State the Yule–Walker equations of an AR($p$).

   Show answer

5. What is $\rho(1)$ of an MA(1), and what is its largest possible value?

   Show answer

6. Which model has an ACF that cuts off after lag $q$, and which has a PACF that cuts off after lag $p$?

   Show answer

7. What does invertibility mean for an MA(1), and why do we require it?

   Show answer

8. What is the mean of a stationary AR(1) with $c = 2$ and $\phi = 0.8$?

   Show answer

9. What does AIC trade off?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you check that an AR model is stationary?" — Roots of $\phi(z)$ outside the unit circle, or companion eigenvalues inside. For AR(2), the triangle. Coefficients below 1 are not enough.
• "How do you choose $p$ and $q$?" — Look at the ACF/PACF for a first guess (cut-offs), fit a small grid, compare AICc on the same data, and check that the residuals are white noise (Ljung–Box). Finally, compare forecasts with a backtest (lesson 03).
• "What is the difference between an MA model and a moving-average smoother?" — The MA model is a stochastic process built from past shocks; the smoother averages past observations.
• "Why is an MA model not identifiable without invertibility?" — $\theta$ and $1/\theta$ (with rescaled noise variance) give the same autocovariances and, for Gaussian data, the same likelihood.
• "How does an AR model relate to a linear RNN?" — An AR($p$) in companion form is a linear RNN with a fixed transition matrix; stability is the same eigenvalue condition.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. simulate_arma(phi, theta, n, burn, rng): simulate an ARMA with burn-in and compare its sample ACF with the theoretical ACF from the $\psi$-weights.
2. is_stationary(phi): decide stationarity with the companion eigenvalues; test it on the corners of the AR(2) triangle.
3. yule_walker(y, p): estimate an AR($p$) with scipy.linalg.solve_toeplitz and compare with least squares on a long AR(2) sample.
4. fit_ar_ols(y, p) and AIC for $p = 1, \dots, 10$ on the sunspots; the implied period of the AR(2).
5. MA(1) invertibility: show that $(\theta, \sigma^2)$ and $(1/\theta, \theta^2\sigma^2)$ give the same autocovariances, and recover the shocks with the invertible model.

An optional cell fits the same models with statsmodels ARIMA(order=(p, 0, q)) if it is installed.

11. Watch and read #

Time Series Talk : Autoregressive Modelritvikmath · 9 min

The AR(p) equation and its intuition, as in section 3.

Time Series Talk : Moving Average Modelritvikmath · 7 min

The MA(q) model as a sum of past shocks; compare with the cut-off ACF in the figure.

Time Series Talk : Moving Average and ACFritvikmath · 9 min

Derives why the ACF of an MA model cuts off.

Time Series Talk : ARMA Modelritvikmath · 7 min

How the AR and MA parts combine.

Vector Auto Regression : Time Series Talkritvikmath · 8 min

For the VAR block in section 6.

• Hyndman and Athanasopoulos, Forecasting: Principles and Practice, 3rd ed. (free online): §9.1 Stationarity, §9.2 Backshift notation, §9.3 Autoregressive models, §9.4 Moving average models, §9.5 Non-seasonal ARIMA models (ACF/PACF reading) and §9.6 Estimation and order selection.
• Shumway and Stoffer, Time Series Analysis and Its Applications, 4th ed., Springer 2017, ch. 3.1–3.5: causality, invertibility, the $\psi$-weights, Yule–Walker and the PACF, with proofs. R code for every example is on the authors' page.
• Penn State, STAT 510 Applied Time Series Analysis, lessons 1–3: free course notes on AR, MA and identification with the ACF and PACF.
• G. U. Yule (1927), "On a method of investigating periodicities in disturbed series, with special reference to Wolfer's sunspot numbers", Phil. Trans. R. Soc. A 226, 267–298 (doi:10.1098/rsta.1927.0007): the paper that introduced autoregression.