Maximum value of ACF at lag 1 for $\text{MA}(q)$ process

Question

Given an $\text{MA}(1)$ process with parameter $\theta_1$, we know that $$ \rho(1) = \frac{\theta_1}{1+\theta_1^2} $$ which has a maximum value of 0.5 when $\theta_1 = 1$. I saw somewhere (I don't remember where) that for an $\text{MA}(q)$ process $\rho(h)\leq \cos(\pi/(q+2))$ for any set of parameters, given the equation $$ Y_t = \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q} + e_t $$ Then, for some set of parameters, it is true that $\rho(h) = \cos(\pi/(q+2))$. It is also true that for an $\text{MA}(q)$ process one can calculate $\rho(1)$ as $$ \rho(1) = \frac{\sum_{k=0}^{q-1} \theta_k \theta_{k+1}}{\sum_{k=0}^{q}\theta_k^2} $$ with $\theta_0 = 1$.

Checking for the MA(1) we indeed have that $\rho(1) = \cos(\pi/3) = 0.5$ when $\theta_1 = 1$. Moreover, for $q = 2$ using $\theta_1 = \sqrt{2}$ and $\theta_2 = 1$ we get $$ \begin{aligned} \rho(1) &= \frac{\sqrt{2} + \sqrt{2}\times 1}{1+2+1} \\ &= \frac{2\sqrt{2}}{4} \\ &= \frac{\sqrt{2}}{2} = \cos(\pi /(2+2)) \end{aligned} $$ I went up until the MA(4) finding that for $\theta = (\sqrt{3},2,\sqrt{3},1)$ one has $$ \begin{aligned} \rho(1) &= \frac{\sqrt{3} + \sqrt{3}\times 2 + 2\times \sqrt{3} + \sqrt{3}\times 1}{1+3+4 + 3+1} \\ &= \frac{6\sqrt{3}}{12} \\ &= \frac{\sqrt{3}}{2} = \cos(\pi /(4+2)) \end{aligned} $$ I believe the solution for the MA(5) case is given by $\theta \approx (1.80194, 2.24698,2.24698,1.80194,1)$. I think also the solution always has $\theta_q = 1$, and the intermediate solutions can be obtained through polynomials. For instance, for this last case the first and (almost) last values $\theta_1$ and $\theta_4$ are close to the roots of the polynomial $x^3 - x^2 - 2x +1$ and the values of $\theta_2$ and $\theta_3$ are close to the roots of the polynomial $x^3 - 2x^2 - x +1$, which are pretty similar.

My question is, is there a general result for this? Basically, what are the set of parameters $\theta = (\theta_1, \theta_2,...,\theta_q)$ of an $\text{MA}(q)$ process such that $\rho(1) = \cos(\pi/(q+2))$.

EDIT: Had $2/\sqrt{2}$ instead of $\sqrt{2}/2$ when calculating $\rho(1)$ for the MA(2). Changed it.

EDIT: I also confirmed at least up to 6 decimal points that indeed for the MA(5) process using $\theta \approx (1.80194, 2.24698,2.24698,1.80194,1)$ gives out $$ \rho(1) = 0.9009689 \approx \cos(\pi / (5+2)). $$

EDIT: As pointed by @whuber, the general formula for $\rho(1)$ of a $\text{MA}(q)$ process was incorrect. Had $$ \rho(1) = \frac{\sum_{k=0}^{q-1} \theta_k \theta_{k+1}}{\sum_{k=0}^{q}\theta_k}. $$ The correct version is $$ \rho(1) = \frac{\sum_{k=0}^{q-1} \theta_k \theta_{k+1}}{\sum_{k=0}^{q}\theta_k^2}. $$

Answer 1

This question concerning what weights $\theta_i$ realize the most extreme possible values of the lag-1 correlation is connected to many parts of mathematics -- but we can answer it using familiar results in statistics concerning covariance, correlation and Principal Components Analysis (PCA).

---

The sequence $(e_t) = \ldots, e_{-2}, e_{-1}, e_0, e_1, \ldots, e_t, e_{t+1}, \ldots$ consists of uncorrelated unit-variance, zero-mean random variables. The MA$(q)$ process $(Y_t)$ (MA = "Moving Average") simply is the windowed linear combination (generalizing the weighted average) of these values with coefficients $\theta = (\theta_0=1,\theta_1,\theta_2, \ldots, \theta_q).$ That is,

$$Y_t = \theta_0 e_t + \theta_1 e_{t-1} + \cdots + \theta_q e_{t-q} = \theta \, (e_t, e_{t-1}, \ldots, e_{t-q})^\prime = \theta\, \mathbf e_t^\prime$$

in matrix notation where $\mathbf e_t = (e_t,e_{t-1},\ldots,e_{t-q})$ is the "window" of length $q+1$ ending at $e_t.$ Because $(e_t)$ is second-order stationary, so is $(Y_t).$

Because $(Y_t)$ is second-order stationary, the lag-1 covariance $R(1) = \operatorname{Cov}(Y_t, Y_{t+1})$ does not depend on $t.$ Since covariance is bilinear,

$$R(1) = \operatorname{Cov}(\theta\, \mathbf e_t^\prime, \theta\, \mathbf e_{t+1}\prime) = \theta \operatorname{Cov}(\mathbf e_t^\prime, \mathbf e_{t+1})\theta^\prime = \frac{1}{2}\theta \pmatrix{0&1&0&0&\cdots&0\\1&0&1&0&\cdots&0\\0&1&0&1&\cdots&0\\\vdots&\vdots&\vdots&\ddots & \cdots & 0\\ 0&0&\cdots & 1 & 0 & 1\\0&0&\cdots&0&1&0}\theta^\prime. $$ Let's call this matrix $\mathbb C_{q},$ so that $2R(1)=\theta\,\mathbb C_q\,\theta^\prime.$

Similarly, the variance of $(Y_t)$ does not depend on $t$ and works out to $$\operatorname{Var}(Y_t) = ||\theta||^2 = \theta_0^2 + \theta_1^2 + \cdots + \theta_q^2.$$

A standard formula for the correlation divides the covariance by the norms of its components:

$$2\rho(1) = \frac{2\operatorname{Cov}(Y_t,Y_{t+1})}{\sqrt{\operatorname{Var}(Y_t)\operatorname{Var}(Y_{t+1})}} = 2\frac{\theta_0\theta_1 + \theta_1\theta_2 + \cdots + \theta_{q-1}\theta_q}{\theta_0^2+\theta_1^2+\cdots \theta_q^2}= \frac{\theta\, \mathbb C_q\, \theta^\prime}{||\theta||^2}.$$

The question asks how to choose $\theta$ to extremize $\rho(1).$ This is the basis of PCA: it tells us the maximum is the largest eigenvalue of $\mathbb C_{q}$ and the minimum is half the smallest eigenvalue of $\mathbb C_q.$

(If you wish to object that $\mathbb C_q$ is not positive definite, you would be correct -- but it makes no difference because you can perform PCA on the symmetric positive-definite matrix $2 - \mathbb C_q$ (with 2's on the diagonal and -1's just off the diagonal) instead.)

There are many ways to find the eigenvalues of $\mathbb C_q.$ Because it is a tridiagonal Toeplitz matrix with $a=0$ on the diagonal and $b=c=1$ on the sub- and super-diagonals, its eigenvalues are

$$\lambda_k = a + 2\sqrt{bc}\cos\left(\frac{\pi k}{q+2}\right) = 2\cos\left(\frac{\pi k}{q+2}\right); \ k = 1, 2, \ldots, q+1.$$

Dividing these values by 2, we find

The extremes of $\rho(1)$ for an MA$(q)$ process are $\lambda_1$ and $\lambda_q,$ equal to $\pm \cos\left(\frac{\pi}{q+2}\right).$ The extremes are attained when $\theta$ is a corresponding eigenvector of $\mathbb C_q,$ normalized to $\theta_0 = 1.$

(These eigenvalues are the roots of Chebyshev polynomials of the second kind,, $S_{q+1}(x) = U_{q+1}(x/2).$)

From the eigenvalue equation $\mathbb C_q \theta = \lambda \theta$ for any eigenvalue $\lambda$ and the normalization $\theta_0=1$ we may recursively find the eigenvectors (when $\lambda \ne 0$) via

$$\theta_0 = 1;\ \theta_1 = \lambda;\ \theta_{i} = \lambda \theta_{i-1} - \theta_{i-2}, i=2, 3, \ldots, q.$$

When $\lambda = 0$ (which occurs when $q$ is even), the eigenvector is $(1,\ldots,1,0,-1,\ldots,-1)^\prime.$

This answers the question. (Although it the recursion doesn't look like a closed formula, it is both simple and efficient, taking only $O(q)$ operations to compute the weights.)

For example, with $q = 5$ and $\lambda = \pm 2\cos\left(\frac{\pi}{q+2}\right) \approx \pm 1.801938$ we obtain

$$\theta \approx (1.000000, \pm 1.801938, 2.246980, \pm 2.246980, 1.801938, \pm 1.000000)$$

for the weights in the MA$(5)$ processes with extremal values of $\rho(1).$

---

Here is the ACF of a long realization of such a process along with shorter examples of various eigenvalues, including the two extremal values in the right column. The realization was long enough to demonstrate that the correlations at lags $0,1,\ldots, q$ are all nonzero and is consistent with the correlations at larger lags being zero.

The examples of this process use a common realization of the random variables $e_t:$ that is, they differ only in the weights used for the moving average.

The following R implementation produced these examples and, with the function theta, provides general-purpose code to compute the weights $\theta.$ The implementation of the MA$(q)$ process in MA.max might be of separate interest because it's just a single short line.

#
# Return the matrix C_q, q >= 1.
#
C <- function(q) {
  zero <- c(1/2, rep(0, q+1))
  X <- matrix(c(0, rep(zero, q)), q+1)
  X + t(X)
}
#
# Return the weights for eigenvector `j` of C_q.
#
theta <- Vectorize(function(q, j = 1) {
  # j can be 1, 2, ..., q+1.
  lambda <- 2 * cos(pi * j / (q+2))
  if (zapsmall(c(lambda, 1), floor(-log10(.Machine$double.eps))-3)[1] == 0) {
    e <- c(rep(1, q/2), 0, rep(-1, q/2))
  } else {
    e <- c(1, lambda, rep(NA_real_, q-1))
    if (isTRUE(q >= 2)) for (i in 3:(q+1)) 
      e[i] <- lambda * e[i-1] - e[i-2]
  }
  e
}, "j")
#
# How to find the eigenvalues.
#
q <- 5
(C(q) %*% theta(q, 1:(q+1)) / theta(q, 1:(q+1)))[1,] # Eigenvalues
#
# Generate an MA(q) process of length `n` using eigenvector `j`.
#
MA.max <- function(n, q = 1, j = 1) {
  convolve(rnorm(n + q), theta(q, j), type = "filter")
}
#
# Graphical examples.
#
par(mfrow = c(2,2))
q <- 5
n <- min(1e5, ceiling(sum(theta(q, 1)^2)^2 * 40))
set.seed(17)
acf(MA.max(n, q), lag.max = q + 1, main = bquote(paste("MA(", .(q), ")")))
set.seed(17)
plot(MA.max(200, q), type = "l", 
     xlab = "t", ylab = "y", main=bquote(rho==.(cos(pi/ (q+2)))))
rho <- zapsmall(c(cos(piceiling((q+1)/2)/ (q+2)), 1))[1]
set.seed(17)
plot(MA.max(200, q, ceiling((q+1)/2)), type = "l", 
     xlab = "t", ylab = "y", main=bquote(rho==.(rho)))
set.seed(17)
plot(MA.max(200, q, q+1), type = "l", 
     xlab = "t", ylab = "y", main=bquote(rho==.(cos(pi(q+1)/ (q+2)))))
par(mfrow = c(1,1))