Stationary Time Series & ARMA

Simulate an autoregression with unit marginal variance, $X_t = \phi X_{t-1} + \sqrt{1 - \phi^2}\,Z_t$ with $Z_t$ i.i.d. $\mathcal{N}(0, 1)$, the parameterization of Topic 26’s Example 8. Every $X_t$ is $\mathcal{N}(0, 1)$, whatever $\phi$ is; only the dependence between neighbors changes. Take $n = 200$ observations and form the textbook 95% interval for the mean, $\bar x \pm 1.96\,s/\sqrt n$, where $\bar x$ is the sample mean and $s$ the sample standard deviation. Over 4,000 replications the interval covers the true mean $\mu = 0$:

• 0.946 of the time at $\phi = 0$ (independent data: the interval works);
• 0.331 of the time at $\phi = 0.9$, against the theoretical value $2\Phi\bigl(1.96\sqrt{(1 - \phi)/(1 + \phi)}\bigr) - 1 = 0.347$, where $\Phi$ is the standard Normal CDF;
• 0.9995 of the time at $\phi = -0.5$: negative dependence makes neighbors cancel, so the naive interval over-covers.

An interval built with the true long-run variance (Theorem 2) covers 0.958 at $\phi = 0.9$. The standard econometric fix, a Newey–West interval with bandwidth $M = 4$ (Remark 6), recovers only to 0.621. At $\phi = 0.9$ the variance of $\bar X$ is about $19/n$, not $1/n$, so 200 observations carry the information of about 10.5 independent ones.

Two forecasting models are compared on two years of daily losses. The loss difference $d_t = \ell^{(A)}_t - \ell^{(B)}_t$ is positively autocorrelated, because both models miss the same regime shifts. If the differences are as persistent as Example 1’s $\phi = 0.9$ series, the naive standard error $s/\sqrt n$ of their mean is too small by a factor $\sqrt{\sigma^2_{\text{LR}}/\gamma(0)} = \sqrt{19} \approx 4.36$, where $\sigma^2_{\text{LR}}$ is Theorem 2’s long-run variance. A “significant” improvement can be pure dependence; Remark 21 says how to estimate the right standard error.

A time series is a sequence of random variables $\{X_t : t \in \mathbb{Z}\}$ indexed by the integers; we call it the process. The data are one realization $x_1, \dots, x_n$ of $n$ consecutive values. We write $\mathbb{E}[X]$ for the expectation (the center of mass of $X$‘s distribution), $\mathrm{Var}(X) = \mathbb{E}[(X - \mathbb{E}X)^2]$ for the variance, and $\mathrm{Cov}(X, Y) = \mathbb{E}[(X - \mathbb{E}X)(Y - \mathbb{E}Y)]$ for the covariance. Time is discrete: $t$ counts steps.

Everything second-order lives in $\gamma(h)$. Sections 35.2–35.6 read it in time: averages, prediction, Wold, ARMA, estimation. Sections 35.7–35.8 read it in frequency: the spectral density, the periodogram, smoothing. Section 35.9 shows what $\gamma$ cannot see: a GARCH series has the autocovariance of white noise, yet predictable volatility.

A process $\{X_t\}$ is strictly stationary if shifting time leaves every finite-dimensional distribution unchanged: for every $k \ge 1$, all times $t_1, \dots, t_k$ and every shift $h$, the vector $(X_{t_1 + h}, \dots, X_{t_k + h})$ has the same joint distribution as $(X_{t_1}, \dots, X_{t_k})$.

A process $\{X_t\}$ is weakly stationary if $\mathbb{E}X_t^2 < \infty$, the mean $\mathbb{E}X_t = \mu$ does not depend on $t$, and the covariance between two observations depends only on how far apart they are:

The function $\gamma$ is the autocovariance function (ACVF) and $\rho(h) = \gamma(h)/\gamma(0)$ is the autocorrelation function (ACF). Directly from the definition, $\gamma(0) = \mathrm{Var}(X_t) \ge 0$, $\gamma(-h) = \gamma(h)$, and $|\gamma(h)| \le \gamma(0)$ by the Cauchy–Schwarz inequality.

A process $\{\varepsilon_t\}$ is white noise with variance $\sigma^2$, written $\varepsilon_t \sim \mathrm{WN}(0, \sigma^2)$, if $\mathbb{E}\varepsilon_t = 0$, $\mathrm{Var}(\varepsilon_t) = \sigma^2$, and $\mathrm{Cov}(\varepsilon_s, \varepsilon_t) = 0$ for $s \ne t$. The stronger assumption that the $\varepsilon_t$ are independent and identically distributed is written $\mathrm{IID}(0, \sigma^2)$. White noise is uncorrelated, not necessarily independent; GARCH (§35.9) will be white noise whose squares are correlated.

A function $\gamma : \mathbb{Z} \to \mathbb{R}$ is the autocovariance function of a weakly stationary series if and only if it is even and nonnegative definite: $\sum_{i,j=1}^n a_i\gamma(t_i - t_j)a_j \ge 0$ for all $n$, all real $a_1, \dots, a_n$, and all integers $t_1, \dots, t_n$.

Simulate $n = 400$ observations from four processes and compare the sample ACF (Remark 5) with the theory.

• White noise: $\rho(h) = 0$ for every $h \ne 0$. The sample autocorrelations scatter inside the band $\pm 1.96/\sqrt n = \pm 0.098$, the approximate 95% range when the true value is zero (Theorem 9).
• AR(1), $\phi = \pm 0.8$: the autoregression $X_t = \phi X_{t-1} + \varepsilon_t$ has $\rho(h) = \phi^{|h|}$, alternating in sign when $\phi < 0$ (Example 5 derives it). At $\phi = 0.8$, $\rho(5) = 0.3277$, but the sample value is 0.4053: sample autocorrelations wander from the truth in runs (Remark 5).
• MA(1), $\theta = 0.8$: the moving average $X_t = \varepsilon_t + \theta\varepsilon_{t-1}$ has $\rho(1) = \theta/(1 + \theta^2) = 0.4878$ (sample 0.4808) and $\rho(h) = 0$ for $h \ge 2$: two observations more than one step apart share no shock. Under this MA(1) the correct band beyond lag 1 is wider than the white-noise band, $\pm 0.1185$, because Bartlett’s formula (Theorem 9) adds the variance contributed by $\rho(1)$.

Let $\{X_t\}$ be weakly stationary with mean $\mu$ and autocovariance $\gamma$ satisfying $\sum_{h \in \mathbb{Z}}|\gamma(h)| < \infty$, and let $\bar X_n = n^{-1}\sum_{t=1}^n X_t$. Then

Consequently $\mathbb{E}(\bar X_n - \mu)^2 \to 0$ and $\bar X_n \xrightarrow{P} \mu$.

Topic 26’s AR(1) chain $X_{t+1} = 0.9X_t + \sqrt{0.19}\,Z_t$ is stationary with $\gamma(0) = 1$ and $\gamma(h) = 0.9^{|h|}$. Theorem 2’s long-run variance is

which is exactly the integrated autocorrelation time $\tau$ of Topic 26 Def 7. A chain of 5,000 draws is worth $N_{\text{eff}} = 5000/19 \approx 263$ independent draws, matching Topic 26 Example 8’s Geyer estimate.

Topic 10 §10.3 proved the weak law for uncorrelated sequences, which is the special case $\gamma(h) = 0$ for $h \ne 0$, where $\sigma^2_{\text{LR}} = \gamma(0)$. Theorem 2 needs only summable correlations, and consistency alone needs only $\gamma(h) \to 0$, since the Cesàro average $n^{-1}\sum_{|h|<n}(1 - |h|/n)\gamma(h)$ then vanishes. Summability buys the $1/n$ rate.

Strict stationarity plus a finite second moment implies weak. The converse fails in general, except for a Gaussian process, whose finite-dimensional laws are fixed by means and covariances (Topic 8). Most theorems here need only weak stationarity; Theorem 15 (GARCH) works with strict.

The sample autocovariance is

and $\hat\rho(h) = \hat\gamma(h)/\hat\gamma(0)$ is the sample ACF. Dividing by $n$ rather than $n - h$ biases $\hat\gamma(h)$ toward zero at large lags, but makes the sample covariance matrix $\hat\Gamma_n = [\hat\gamma(i - j)]$ nonnegative definite (it equals $n^{-1}TT^\top$ for a matrix $T$ of shifted copies of the demeaned data), which Yule–Walker needs (Remark 15). Sample autocorrelations at nearby lags are strongly correlated (Theorem 9), which is why Example 3’s AR(1) sample ACF stays above the theory for a run of lags.

Example 1’s oracle interval used the true $\sigma^2_{\text{LR}}$. In practice we estimate it, and the plug-in $\sum_{|h| < n}\hat\gamma(h)$ is inconsistent — it equals $2\pi$ times the periodogram at frequency zero, which §35.8 shows never settles down. The Newey–West estimator (NEW1987) truncates and down-weights the sample autocovariances:

With the common rule $M = \lfloor 4(n/100)^{2/9}\rfloor = 4$ at $n = 200$, it lifts Example 1’s coverage from 0.331 to only 0.621 at $\phi = 0.9$: the truncation lag $M$ is a bandwidth, too small for a series this persistent (Remark 21). This is the HAC covariance used in GMM with dependent data: formalML: generalized-method-of-moments .

A linear process is a series of the form

It is causal if $\psi_j = 0$ for $j < 0$, so that $X_t$ depends only on present and past shocks. The infinite sum converges in $L^2$ and almost surely (Topic 9), and because the covariance is continuous in $L^2$ we may compute it term by term: by bilinearity and the uncorrelatedness of white noise, $\gamma(h) = \sigma^2\sum_j\psi_j\psi_{j+|h|}$.

A moving average of order $q$, $X_t = \varepsilon_t + \theta_1\varepsilon_{t-1} + \dots + \theta_q\varepsilon_{t-q}$, is a linear process with $\psi_0 = 1$, $\psi_j = \theta_j$ for $1 \le j \le q$, and $\psi_j = 0$ otherwise. Since $\psi_j\psi_{j+|h|} = 0$ for $|h| > q$, its autocovariance vanishes beyond lag $q$.

An AR(1) $X_t = \phi X_{t-1} + \varepsilon_t$ with $|\phi| < 1$ becomes a linear process by substituting the equation into itself $k$ times:

The remainder has $\mathbb{E}(\phi^kX_{t-k})^2 = \phi^{2k}\gamma(0) \to 0$, so in $L^2$, $X_t = \sum_{j \ge 0}\phi^j\varepsilon_{t-j}$: $\psi_j = \phi^j$. Definition 4’s formula then gives

which discharges Example 3’s promise that $\rho(h) = \phi^{|h|}$. For $\phi = 0.9$ and $\sigma^2 = 1$, $\gamma(0) = 1/0.19 = 5.263158$.

Let $\{X_t\}$ be zero-mean and weakly stationary. Write $\mathcal{H}_n = \overline{\mathrm{sp}}\{X_s : s \le n\}$ for the closed linear span of the infinite past (all $L^2$ limits of finite linear combinations). The best linear predictor of $X_{n+1}$ is the projection $\hat X_{n+1} = P_{\mathcal{H}_n}X_{n+1}$, the innovation is $\varepsilon_{n+1} = X_{n+1} - \hat X_{n+1}$, and the one-step mean squared error is $\sigma^2 = \mathbb{E}\varepsilon_{n+1}^2$, which does not depend on $n$ by stationarity. The series is deterministic if $\sigma^2 = 0$ (the past predicts the future perfectly) and purely nondeterministic if $\bigcap_n\mathcal{H}_n = \{0\}$ (the remote past carries no information).

With only finitely many observations, $P_nX_{n+1}$ denotes the projection onto $\mathrm{sp}\{X_1, \dots, X_n\}$, with mean squared error $v_n$.

Let $\{X_t\}$ be zero-mean weakly stationary with $\gamma(0) > 0$ and $\gamma(h) \to 0$ as $h \to \infty$. Write $\Gamma_n = [\gamma(i - j)]_{i,j=1}^n$ and $\boldsymbol\gamma_n = (\gamma(1), \dots, \gamma(n))^\top$.

(i) The best linear predictor of $X_{n+1}$ from $X_1, \dots, X_n$ is $P_nX_{n+1} = \sum_{j=1}^n\phi_{nj}X_{n+1-j}$, where $\boldsymbol\phi_n = (\phi_{n1}, \dots, \phi_{nn})^\top$ is the unique solution of $\Gamma_n\boldsymbol\phi_n = \boldsymbol\gamma_n$, and its mean squared error is $v_n = \gamma(0) - \boldsymbol\phi_n^\top\boldsymbol\gamma_n$.

(ii) With $v_0 = \gamma(0)$ and $\phi_{11} = \gamma(1)/\gamma(0)$, for $n \ge 2$:

For a causal linear process in its own innovations, the $h$-step forecast drops the shocks not yet seen, so its mean squared error is $\sigma^2\sum_{j<h}\psi_j^2$. For the AR(1), $P_nX_{n+h} = \phi^hX_n$ and the MSE is $\gamma(0)(1 - \phi^{2h})$. As a fraction of the variance $\gamma(0)$:

• AR(1), $\phi = 0.9$: 0.19 at $h = 1$, so the past explains 81% of the variance (0.81); 0.651322 at $h = 5$.
• AR(2) with roots $1 \pm i$: 0.416667 at $h = 1$, 0.963542 at $h = 5$.
• MA(2), $\theta = (0.6, 0.4)$: 0.657895 at $h = 1$, and exactly 1 from $h = 3$ on. An MA($q$) is unpredictable beyond $q$ steps.

Let $\{X_t\}$ be zero-mean weakly stationary, $\mathcal{H}_t = \overline{\mathrm{sp}}\{X_s : s \le t\}$, and suppose the one-step prediction error $\sigma^2 = \mathbb{E}(X_t - P_{\mathcal{H}_{t-1}}X_t)^2$ is positive. Then

where (i) $\varepsilon_t = X_t - P_{\mathcal{H}_{t-1}}X_t$ is $\mathrm{WN}(0, \sigma^2)$ and $\varepsilon_t \in \mathcal{H}_t$; (ii) $\psi_0 = 1$, $\psi_j = \mathbb{E}[X_t\varepsilon_{t-j}]/\sigma^2$ does not depend on $t$, and $\sum_j\psi_j^2 < \infty$; (iii) $\mathbb{E}[\varepsilon_sV_t] = 0$ for all $s, t$; (iv) $V_t \in \mathcal{H}_{-\infty} = \bigcap_s\mathcal{H}_s$; (v) $\{V_t\}$ is deterministic: $V_t$ lies in the closed span of $\{V_s : s < t\}$.

Let $X_t = A\cos(\lambda t) + B\sin(\lambda t)$, where $A$ and $B$ are uncorrelated with mean 0 and common variance $\sigma_A^2$, and $\lambda \in (0, \pi)$ is a fixed frequency. The angle-addition formula gives $\gamma(h) = \sigma_A^2\cos(\lambda h)$, so the series is stationary. It is also exactly predictable: the trigonometric identity $\cos(\lambda t) = 2\cos\lambda\cos(\lambda(t-1)) - \cos(\lambda(t-2))$, and the same for sines, gives

with no error term, so $\sigma^2 = 0$ and Wold’s $V_t$ is the whole series. Its spectrum (§35.7) is two point masses at $\pm\lambda$, a spectral distribution without a density.

The innovations are uncorrelated, not independent. A GARCH series (§35.9) is white noise, so its Wold representation is trivial ($\psi_0 = 1$, all other $\psi_j = 0$), yet its squared returns are predictable. Linear prediction is blind to volatility; that blindness is §35.9’s subject.

For a Gaussian series the joint density of $x_1, \dots, x_n$ factors into one-step predictions:

where $\hat x_t$ is the best linear predictor of $x_t$ from $x_1, \dots, x_{t-1}$ and $v_{t-1}$ its mean squared error. The exact likelihood is a product of one-step prediction errors. Theorem 3 computes the $\hat x_t$ and $v_t$ in $O(n^2)$ operations; Topic 36 computes them in $O(n)$ with the Kalman filter, after writing an ARMA model in state-space form.

If $\{X_t\}$ is a Gaussian process, $P_nX_{n+1} = \mathbb{E}[X_{n+1} \mid X_1, \dots, X_n]$: the best linear predictor is the best predictor of any kind. This is the conditional multivariate Normal formula of Topic 8 §8.4, with $\Gamma_n$ as the conditioning covariance block. Outside the Gaussian case the conditional mean can beat every linear predictor.

Write $B$ for the backshift operator, $BX_t = X_{t-1}$, so $B^kX_t = X_{t-k}$. Define the AR polynomial and MA polynomial

A zero-mean series $\{X_t\}$ is an ARMA($p$, $q$) process if it is stationary and

that is, $X_t - \phi_1X_{t-1} - \dots - \phi_pX_{t-p} = \varepsilon_t + \theta_1\varepsilon_{t-1} + \dots + \theta_q\varepsilon_{t-q}$. We require that $\phi$ and $\theta$ have no common zeros: a shared factor cancels from both sides and the model is over-parameterized. A series with mean $\mu$ is ARMA if $X_t - \mu$ is.

An ARMA process is causal if $X_t = \sum_{j \ge 0}\psi_j\varepsilon_{t-j}$ with $\sum_j|\psi_j| < \infty$: the present is built from present and past shocks only. It is invertible if $\varepsilon_t = \sum_{j \ge 0}\pi_jX_{t-j}$ with $\sum_j|\pi_j| < \infty$: the shocks can be recovered from present and past observations. The $\psi_j$ are the MA(∞) weights and the $\pi_j$ the AR(∞) weights.

Let $\phi(B)X_t = \theta(B)\varepsilon_t$ with $\varepsilon \sim \mathrm{WN}(0, \sigma^2)$, where $\phi$ and $\theta$ have no common zeros. A causal stationary solution $X_t = \sum_{j \ge 0}\psi_j\varepsilon_{t-j}$ with $\sum_j|\psi_j| < \infty$ exists if and only if $\phi(z) \ne 0$ for all complex $|z| \le 1$. In that case the solution is unique and $\sum_j\psi_jz^j = \theta(z)/\phi(z)$ for $|z| \le 1$.