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 n) time. In particular, for a Toeplitz matrix with values ₉, ₊ = h|₉-₊| = -₁/₂^1/2 e^2 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 C^n 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 R^10⁵ 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.