Solver for the true auto-covariance function in AR(p)

Question

Suppose I have the following $AR(p)$ model.

$$X_t = \sum_{i=1}^{p} \phi_i X_{t-i} + \epsilon_t\,, $$

where $\epsilon_t$ has mean 0 variance $\sigma^2$. I am in the situation where the $\phi$s are known and my goal is to obtain the true auto-covariance $$ \gamma(k) = \text{Cov}(X_1, X_{1+k})\,. $$ (I am not interested in estimating $\gamma(k)$). The spectral density at $0$ for AR$(p)$ model is $$ f(0) = \sum_{k=-\infty}^{\infty} \gamma(k) = \dfrac{\sigma^2}{(1 - \sum_{i=1}^{p} \phi_i)^2}\, $$ which is available in closed-form to me, since all of $\sigma^2$ and $\phi_i$ are known to me. Now, by the Yule-Walker equations, for $k = 1, \dots, p$ $$ \gamma(k) = \sum_{i=1}^{p} \phi_i \gamma(k-i) $$ and $\gamma(0) = \sum_{i=1}^{p} \phi_i \gamma(k-i) + \sigma^2$.

Thus obtaining the true $\gamma(k)$ for $k = 0, \dots, p$ will require solving the above system of equations. I have two questions:

1. Is there an off-the-shelf R/Python/Matlab function available that outputs $\gamma(k)$ if I give it $\phi_i$ and $\sigma^2$?
2. What are the higher lag covariances: $\gamma(k)$ for $k > p$?

Answer 1

You can use ARMA.autocov in the ts.extend package

The ts.extend package contains a number of functions that compute theoretical aspects of stationary ARMA models, including the auto-correlation and auto-covariance functions, the auto-covariance matrix, and the standard probability functions for the stationary Gaussian ARMA model (i.e., density, distribution, and random-generation function). The ARMA.autocov function takes an auto-regression vector ar and a moving-average vector ma (so long as the AR polynomial is for a stationary model) and gives the auto-correlation or auto-covariance function up to any specified length n. (This function extends a similar function in the stats package, which lets you get the auto-correlation, but not the auto-covariance.$^\dagger$) This can be implemented using the following syntax.

#Set the parameters
AR    <- c(0.8, -0.4, 0.2, 0.1)
sigma <- 5

#Compute auto-covariance function up to length n
AUTOCOV <- (sigma^2)*ts.extend::ARMA.autocov(n = 10, ar = AR)

#Plot the auto-correlation function
barplot(AUTOCOV, 
        main = 'Autocovariance Function\n(for an AR(4) model)',
        ylab = 'Autocovariance')

Note here that the parameter $\sigma$ is not part of the ARMA.autocov function, but it affects the auto-covariance by using multiplication by $\sigma^2$.

---

$^\dagger$ The function stats::ARMAacf computes the auto-correlation function, but there does not appear to be any function in the stats package to compute the auto-covariance function, which is quite annoying. If you can separately compute the variance of the time-series values (i.e., the first element of the auto-covariance function) then you can multiply the auto-correlation function to get the auto-covariance function, but the stats package does not appear to give you a way to do this.