Bayesian multiscale feature detection of log-spectral densities

Abstract

A fully-automatic Bayesian visualization tool to identify periodic components of evenly sampled stationary time series, is presented. The given method applies the multiscale ideas of the SiZer-methodology to the log-spectral density of a given series. The idea is to detect significant peaks in the true underlying curve viewed at different resolutions or scales. The results are displayed in significance maps, illustrating for which scales and for which frequencies, peaks in the log-spectral density are detected as significant. The inference involved in producing the significance maps is performed using the recently developed simplified Laplace approximation. This is a Bayesian deterministic approach used to get accurate estimates of posterior marginals for latent Gaussian Markov random fields at a low computational cost, avoiding the use of Markov chain Monte Carlo techniques. Application of the given exploratory tool is illustrated analyzing both synthetic and real time series.

Introduction

Many natural phenomena show a repetitive or regular behavior over time, best described as functions of periodic variations. Periodicities of stationary time series are most easily studied applying a frequency-domain approach, analyzing the series by means of their spectral representations. The spectral representation of stationary time series essentially decomposes the series into a sum of sines and cosines with uncorrelated random coefficients. The spectral density function, also named the power spectrum, reflects the periodic nature of a time series and describes the variance of the series in terms of frequency. Peaks in the power spectrum indicate which periodic components the time series is composed of, giving valuable information about the underlying physical process. Spectral density analysis is an important subject in various fields, for example statistical signal processing, economic analyses or climatology. For a thorough mathematical introduction to spectral density analysis, see for example Brockwell and Davis (1991) and Brillinger (2001).

The power spectrum of a zero-mean stationary real-valued stochastic process $x_t,t\in Z$, is defined by

$$f\left(\omega \right)=\sum _{h=-\infty }^{\infty }\gamma \left(h\right)e^{-2\pi i\omega h}\text{,}-1/2\le \omega \le 1/2\text{,}$$

where the autocovariance function (ACF) of the series, $\gamma \left(h\right)=E\left(x_t{x}_{t+h}\right),h\in Z$ is assumed to be absolutely summable, i.e. $\sum \left|\gamma \left(h\right)\right|<\infty$. We will assume stationarity in a wide-sense, having a stationary ACF. The power spectrum, being the Fourier transform of the ACF, provides exactly the same information on the second-order properties of the series as the ACF. We also note that the total variance equals the spectrum integrated over all frequencies.

An aim in spectral density analysis is to estimate the power spectrum of a random signal based on an observed time series. The field of suggested estimators is extensive, including both parametric, nonparametric and Bayesian methods. The classical nonparametric periodogram estimator is given by

$$I\left(\omega \right)=\frac{1}{2n}{\left|\sum _{t=1}^{2n}{x}_t{e}^{-2\pi i\omega t}\right|}^2\text{,}-1/2\le \omega \le 1/2\text{,}$$

where $x_1,\dots ,x_{2n}$ denote the observed time series. The raw periodogram is known to be asymptotically unbiased but not consistent, giving a highly fluctuating estimate. In practice, the periodogram is usually smoothed, both to obtain a consistent estimator and to make the resulting estimate more easy to interpret.

In most of the early literature on finding a consistent estimator of the power spectrum, the periodogram is smoothed directly, applying weighted local average approaches, see for example Blackman and Tukey (1958), Parzen (1961) and Priestley (1981). Other well-known consistent nonparametric estimators include procedures based on the Whittle (1962) likelihood, and also least-squares based methods for the log-periodogram, see for example Wahba (1980). The presented method is based on smoothing the log-periodogram applying a Bayesian approach, including an integrated Wiener process as a prior. Other methods for smoothing the periodogram within a Bayesian framework include, among others, Carter and Khon (1997), Liseo et al. (2001), Choudhuri et al. (2004) and Tonellato (2007), all eliciting different choices of prior models.

The aim of the presented method is to identify significant peaks in the log-spectral density of stationary time series, using a multiscale approach. Multiscale methods are often referred to as scale-space methods which are well-known in computer vision, see Lindeberg (1994) for an introduction. Multiscale ideas applied to nonparametric function estimation were introduced in Chaudhuri and Marron (1999), presenting the SiZer-methodology (SIgnificant ZERo-crossings of derivatives). The idea of SiZer is to detect significant trends in the true underlying curve viewed at different levels of smoothing, also referred to as scales. The results are visualized in SiZer maps, illustrating significant trends of the true underlying curve for all scales simultaneously, see Erästö (2005) for a discussion of SiZer maps. A theoretical justification of the SiZer-methodology is given in Chaudhuri and Marron (2000).

Originally, quantiles needed to evaluate the SiZer map were computed applying Gaussian distributional assumptions or bootstrapping, see Chaudhuri and Marron (1999). Some computational improvements, based on asymptotic theory, were developed in Hannig and Marron (2006). Bayesian versions of the SiZer-methodology were given in Erästö and Holmström (2005), presenting both an analytically solvable regression model and an approach using Gibbs sampling. In Øigård et al. (2006), a Bayesian approach was presented for Gaussian observational models, calculating quantiles exactly.

We extend the application of the SiZer-methodology to spectral density estimation and to a non-Gaussian observational model, using Bayesian approximate inference. Applying the integrated Wiener process as a prior model, the vector of discretely observed log-spectral densities augmented with first derivatives is a Gaussian Markov random field (GMRF). The simplified Laplace approximation, given in Rue et al. (2009), can then be applied to obtain accurate estimates of the posterior marginals needed to estimate the log-spectral density and its first derivative for different scales. Inference is summarized in color maps similar to SiZer maps, here named significance maps.

The computations involved are performed efficiently, applying numerical algorithms that take advantage of the special conditional independence structure of GMRFs. This class of models is widely used in Bayesian analysis, in both temporal and spatial contexts. Traditionally, Bayesian inference for such models has been performed by Markov chain Monte Carlo (MCMC) techniques. However, MCMC algorithms might have a high computational cost, especially in finding accurate estimates for tail probabilities, having an additive error. In comparison, the simplified Laplace approximation has a relative error. By applying this method, also the tails of a distribution are accurately estimated.

Details on prior choices, the simplified Laplace approximation, implementation and evaluation of probabilities to find significant peaks, are given in Section 2. Examples, illustrating the performance of the given method, are presented in Section 3. Some concluding remarks are given in Section 4.