Fast Machine-Precision Spectral Likelihoods for Stationary Time Series

TL;DR

This paper introduces a fast, highly accurate algorithm for approximating symmetric Toeplitz matrices in spectral likelihood computations, significantly improving time series analysis accuracy and efficiency.

Contribution

The authors develop a novel $ ext{O}(n ext{log} n)$ algorithm for machine-precision approximation of symmetric Toeplitz matrices, enhancing spectral likelihood methods in time series modeling.

Findings

Achieves $ ext{O}(n ext{log} n)$ approximation in machine precision.

Improves Whittle likelihood accuracy from 3 to 14 digits with low-rank correction.

Demonstrates efficiency on large matrices up to size $10^5 imes 10^5$.

Abstract

We provide in this work an algorithm for approximating a very broad class of symmetric Toeplitz matrices to machine precision in $\mathcal{O}(n \log n)$ time with applications to fitting time series models. In particular, for a symmetric Toeplitz matrix $\mathbf{\Sigma}$ with values $\mathbf{\Sigma}_{j,k} = h_{|j-k|} = \int_{-1/2}^{1/2} e^{2 \pi i |j-k| \omega} S(\omega) \mathrm{d} \omega$ where $S(\omega)$ is piecewise smooth, we give an approximation $\mathbf{\mathcal{F}} \mathbf{\Sigma} \mathbf{\mathcal{F}}^H \approx \mathbf{D} + \mathbf{U} \mathbf{V}^H$, where $\mathbf{\mathcal{F}}$ is the DFT matrix, $\mathbf{D}$ is diagonal, and the matrices $\mathbf{U}$ and $\mathbf{V}$ are in $\mathbb{C}^{n \times r}$ with $r \ll n$. Studying these matrices in the context of time series, we offer a theoretical explanation of this structure and connect it to existing spectral-domain approximation frameworks. We then give a complete discussion of the numerical method for assembling the approximation and demonstrate its efficiency for improving Whittle-type likelihood approximations, including dramatic examples where a correction of rank $r = 2$ to the standard Whittle approximation increases the accuracy of the log-likelihood approximation from $3$ to $14$ digits for a matrix $\mathbf{\Sigma} \in \mathbb{R}^{10^5 \times 10^5}$. The method and analysis of this work applies well beyond time series analysis, providing an algorithm for extremely accurate solutions to linear systems with a wide variety of symmetric Toeplitz matrices whose entries are generated by a piecewise smooth $S(\omega)$. The analysis employed here largely depends on asymptotic expansions of oscillatory integrals, and also provides a new perspective on when existing spectral-domain approximation methods for Gaussian log-likelihoods can be particularly problematic.