Fast Machine-Precision Spectral Likelihoods for Stationary Time Series

TL;DR

The method and analysis 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)$.

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Σ$ with values $\mathbfΣ_{j,k} = h_{|j-k|} = \int_{-1/2}^{1/2} e^{2 πi |j-k| ω} S(ω) \mathrm{d} ω$ where $S(ω)$ is piecewise smooth, we give an approximation $\mathbf{\mathcal{F}} \mathbfΣ \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Σ \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(ω)$. 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.

Figures (5)

• Figure 1: The runtime cost of assembling the approximation $\mathbb{\mathcal{F}} \bm{\Sigma} \mathbb{\mathcal{F}}^H \approx \bm{D} + \bm{U} \bm{V}^H$ and evaluating (\ref{['eq:appx_nll']}) for various different ranks $r$. We also provide the runtime cost of assembling just the diagonal matrix of $S(\omega)$ at Fourier frequencies ($r=0$) and the debiased Whittle estimator.
• Figure 2: The relative errors $|\ell(\bm{\theta} \, | \, \bm{y}) - \tilde{\ell}(\bm{\theta} \, | \, \tilde{\bm{y}})|/|\ell(\bm{\theta} \, | \, \bm{y})|$ for various models, ranks, and data sizes. The middle column shows the error in naively applying the asymptotic expansion method of Proposition \ref{['thm:asexp']} to obtain $\left\{ h_k \right\}_{k=0}^{n-1}$ despite $S_2(\omega)$ having zero derivatives at the origin, and the right column applies the split expansion method from Corollary \ref{['cor:asexp']}.
• Figure 3: Runtime costs for evaluating the approximated gradient (\ref{['eq:grad']}), the exact Fisher information matrix (\ref{['eq:fish']}), and the symmetrized stochastic expected Fisher information matrix for various data sizes and ranks using $S_2(\omega)$ with asymptotic expansions split at $0$.
• Figure 4: Relative error-type metrics for the approximated gradient and both exact and stochastically approximated Fisher information matrices for various ranks and sizes.
• Figure 5: A comparison of the $50$ trials of fitting a process simulated with a true SDF of $S_2(\omega)$ with $\theta_1 = 5$ and $\theta_2 = 35$. The first row shows a scatter plot of points $(\hat{\theta}_1^{\text{MLE}}, \hat{\theta}_1)$ with $\hat{\theta}_1$ denoting the alternative estimator, and the lower row shows the analog for the estimates of $\theta_2$.

Theorems & Definitions (7)

• Proposition 1
• Proof 1
• Proposition 1
• Proof 2
• Corollary 2
• Proposition 1
• Proof 3