Fast Machine-Precision Spectral Likelihoods for Stationary Time Series

We provide in this work an algorithm for approximating a very broad class of symmetric Toeplitz matrices to machine precision in O(n log n) time. In particular, for a Toeplitz matrix Σ with values Σj,k = h|j-k| = ∫-1/21/2 e2 π i |j-k| ω S(ω) d ω where S(ω) is piecewise smooth, we give an approximation F Σ FH ≈ D + U VH, where F is the DFT matrix, D is diagonal, and the matrices U and V are in Cn × r with r 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 from $3$ to $14$ digits for a matrix Σ ∈ R10⁵ × 10⁵. The method and analysis of this work applies well beyond time series analysis, providing an algorithm for extremely accurate direct solves with a wide variety of symmetric Toeplitz matrices. 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.