AdaptSPEC-X: Covariate Dependent Spectral Modeling of Multiple Nonstationary Time Series

TL;DR

AdaptSPEC-X is a Bayesian covariate-dependent spectral model for analyzing multiple nonstationary time series, capable of handling missing data and estimating time-varying means and spectra with predictive capabilities.

Contribution

It extends AdaptSPEC to incorporate covariates, handle missing data, and model time-varying means, enabling flexible joint analysis and prediction of nonstationary time series.

Findings

Effective in simulated data

Applied successfully to Australian rainfall data

Demonstrated on US measles incidence

Abstract

We present a method for the joint analysis of a panel of possibly nonstationary time series. The approach is Bayesian and uses a covariate-dependent infinite mixture model to incorporate multiple time series, with mixture components parameterized by a time varying mean and log spectrum. The mixture components are based on AdaptSPEC, a nonparametric model which adaptively divides the time series into an unknown number of segments and estimates the local log spectra by smoothing splines. We extend AdaptSPEC to handle missing values, a common feature of time series which can cause difficulties for nonparametric spectral methods. A second extension is to allow for a time varying mean. Covariates, assumed to be time-independent, are incorporated via the mixture weights using the logistic stick breaking process. The model can estimate time varying means and spectra at observed and unobserved covariate values, allowing for predictive inference. Estimation is performed by Markov chain Monte Carlo (MCMC) methods, combining data augmentation, reversible jump, and Riemann manifold Hamiltonian Monte Carlo techniques. We evaluate the methodology using simulated data, and describe applications to Australian rainfall data and measles incidence in the US. Software implementing the method proposed in this paper is available in the R package BayesSpec.

Tables

Table 1. Table 1 : Estimated posterior probabilities p ^ ( ⋅ ∣ 𝒙 ) \hat{p}(\cdot\mid\bm{x}) of various events for monthly rainfall at the four example sites. The events correspond to the WW2 drought (first four rows), the Millenium drought (second two rows), and long term change (last two rows). For each event and site ( 𝒖 𝒖 \bm{u} ), the table presents probabilities that the mean μ ​ ( t , 𝒖 ) 𝜇 𝑡 𝒖 \mu(t,\bm{u}) and the variance σ 2 ​ ( t , 𝒖 ) superscript 𝜎 2 𝑡 𝒖 \sigma^{2}(t,\bm{u}) changed before or after the event.

		Site ($𝒖$)
Event	$\hat{p}(\cdot \mid 𝒙)$	10525	14042	47053	69018
WW2 drought		0.539	0.442	0.700	0.802
		0.572	0.390	0.508	0.833
		0.915	0.980	0.997	0.994
		0.997	0.996	1.000	0.997
Millenium drought		0.878	0.906	0.975	0.942
		0.747	0.932	0.999	0.978
Long term		0.818	0.907	0.940	0.972
		0.961	0.954	0.999	0.989

Equations

$$p(\mathbf{x}\mid \mu ,f)=\frac{1}{(2\pi {)}^{n/2}}\frac{1}{{\prod }_{k=1}^{n}f({\omega }_k{)}^{1/2}}\exp \left\{-\frac{1}{2}\sum _{k=1}^n\frac{I_k}{f({\omega }_k)}\right\},$$

$$d_k=\frac{1}{\sqrt{n}}\sum _{t=1}^n(x_t-\mu )e^{-2\pi i{\omega }_k(t-1)}$$

$$\log f(\omega )={\alpha }_0+h(\omega ),$$

$$X_t=\sum _{s=1}^m{X}_t^s{\delta }_{s,m}(t),$$

$$g(\mathbf{x}\mid \Theta )=\prod _{s=1}^{m}p({\mathbf{x}}_{s,m}\mid {\mu }_{s,m},f_{s,m}),$$

$$\{X_t\}\sim \sum _{h=1}^H{\pi }_h(\mathbf{u})g_h(\{X_t\}\mid {\Theta }_h),$$

$$p({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N)=\prod _{j=1}^N\sum _{h=1}^H{\pi }_h({\mathbf{u}}_j)g_h({\mathbf{x}}_j\mid {\Theta }_h).$$

$${\pi }_h(\mathbf{u})=v_h(\mathbf{u})\prod _{h^{\prime }=1}^{h-1}(1-v_{h^{\prime }}(\mathbf{u})),$$

$${\mathbf{w}}_h\sim \text{N}({\beta }_{0,h}{1}_N+U{\mathbf{\beta }}_h,{\tau }_h^2{\Sigma }_w),$$

$$\begin{array}[]{lll}\bm{w}_{h}&=&U^{\dagger}\bm{\beta}^{\dagger}_{h},\\ \bm{\beta}^{\mathrm{GP}}_{h}&\sim&\text{N}(\bm{0}_{B},\tau^{2}_{h}I_{B}),\end{array}$$

$${\mathbf{\beta }}_h^{†}\sim \text{N}\left(\left(\begin{array}{c}{\mathbf{\mu }}_{\beta }\\ 0_B\end{array}\right),\left(\begin{array}{cc}{\Sigma }_{\beta }& 0\\ 0& {\tau }_h^2{I}_B\end{array}\right)\right).$$

$$p({\tau }_h)\propto {\left(1+\frac{1}{{\nu }_{\tau }}{\left(\frac{{\tau }_h}{A_{\tau }}\right)}^2\right)}^{-({\nu }_{\tau }+1)/2},{\tau }_h>0,$$

$$p(\mathbf{x}\mid \mu ,f)=\frac{1}{(2\pi {)}^{n/2}}|R{|}^{1/2}\exp \left\{-\frac{1}{2}(\mathbf{x}-\mu {)}^{\prime }VR{V}^{\ast }(\mathbf{x}-\mu )\right\},$$

$$({\mathbf{x}}_{\text{mis}}\mid {\mathbf{x}}_{\text{obs}},\mu ,f)\sim \text{N}({\mu }_{\text{mis}\mid \text{obs}},{\Lambda }_{\text{mis}\mid \text{obs}}^{-1}),$$

$$\begin{array}{rl}{\mu }_{\text{mis}\mid \text{obs}}& =\mu -{\Lambda }_{\text{mis},\text{mis}}^{-1}{\Lambda }_{\text{mis},\text{obs}}({\mathbf{x}}_{\text{obs}}-\mu ),\\ {\Lambda }_{\text{mis}\mid \text{obs}}& ={\Lambda }_{\text{mis},\text{mis}},\end{array}$$

$$\Lambda =\left[\begin{array}{cc}{\Lambda }_{\text{mis},\text{mis}}& {\Lambda }_{\text{mis},\text{obs}}\\ {\Lambda }_{\text{mis},\text{obs}}^{\prime }& {\Lambda }_{\text{obs},\text{obs}}\end{array}\right]=\left[\begin{array}{cc}{V}_{\text{mis}}R{V}_{\text{mis}}^{\ast }& V_{\text{mis}}R{V}_{\text{obs}}^{\ast }\\ V_{\text{obs}}R{V}_{\text{mis}}^{\ast }& V_{\text{obs}}R{V}_{\text{obs}}^{\ast }\end{array}\right],$$

$$z_j^{\text{swap}}=\{\begin{array}{ll}{h}_2& \text{if }{z}_j^{[l+0.5]}=h_1,\\ h_1& \text{if }{z}_j^{[l+0.5]}=h_2,\\ z_j^{[l+0.5]}& \text{otherwise,}\end{array}$$

$${\tau }_h^{\text{swap}}=\{\begin{array}{ll}{\tau }_{h_2}^{[l+0.5]}& \text{if }h=h_1\\ {\tau }_{h_1}^{[l+0.5]}& \text{if }h=h_2\\ {\tau }_h^{[l+0.5]}& \text{otherwise,}\end{array}$$

$$q({\mathbf{\beta }}_h^{†\text{swap}}\mid {\mathbf{z}}^{[l+0.5]},{\mathbf{\tau }}^{[l+0.5]})\sim \text{N}({\mu }_h^{\text{mode}},{\Sigma }_h^{\text{mode}}),$$

$$\displaystyle\min\Bigg{\{}1,\frac{p(\bm{\beta}^{\dagger\text{swap}},\bm{z}^{\text{swap}},\Theta^{\text{swap}},\bm{\tau}^{\text{swap}}\mid\bm{x}_{\text{mis}}^{\text{all}[l+1]},\bm{x}^{\text{all}}_{\text{obs}})}{p(\bm{\beta}^{\dagger[l+0.5]},\bm{z}^{[l+0.5]},\Theta^{[l+0.5]},\bm{\tau}^{[l+0.5]}\mid\bm{x}_{\text{mis}}^{\text{all}[l+1]},\bm{x}^{\text{all}}_{\text{obs}})}$$

$$\displaystyle\qquad\qquad\frac{q(\bm{\beta}^{\dagger[l+0.5]}_{h_{1}}\mid\bm{z}^{\text{swap}},\bm{\tau}^{\text{swap}})}{q(\bm{\beta}^{\dagger\text{swap}}_{h_{1}}\mid\bm{z}^{[l+0.5]},\bm{\tau}^{[l+0.5]})}\frac{q(\bm{\beta}^{\dagger[l+0.5]}_{h_{2}}\mid\bm{z}^{\text{swap}},\bm{\tau}^{\text{swap}})}{q(\bm{\beta}^{\dagger\text{swap}}_{h_{2}}\mid\bm{z}^{[l+0.5]},\bm{\tau}^{[l+0.5]})}\Bigg{\}}.$$

$$(x_{j,t}-{\mu }_{z_j,t})={\varphi }_{z_j,1,t}(x_{j,t-1}-{\mu }_{z_j,t-1})+{\varphi }_{z_j,2,t}(x_{j,t-2}-{\mu }_{z_j,t-2})+{ϵ}_{j,t},$$

$${\text{MSE}}_{\text{mean}}(\mathbf{u})=\frac{1}{n}\sum _{t=1}^n{\left[\hat{\mu }(t,\mathbf{u})-\mu (t,\mathbf{u})\right]}^2,$$

$${\text{MSE}}_{\text{spec}}(\mathbf{u})=\frac{1}{n}\frac{1}{k_{\text{max}}}\sum _{t=1}^n\sum _{k=1}^{k_{\text{max}}}{\left[\log \hat{f}\left(t,\frac{k-1}{2k_{\text{max}}-2},\mathbf{u}\right)-\log f\left(t,\frac{k-1}{2k_{\text{max}}-2},\mathbf{u}\right)\right]}^2,$$

$$p({\Theta }_h\mid \mathbf{z},{\mathbf{x}}^{\text{all}})\propto p({\Theta }_h)\prod _{\{j:z_j=h\}}{g}_h({\mathbf{x}}_j\mid {\Theta }_h).$$

$$\log f({\omega }_k)={\alpha }_0+h({\omega }_k)={\mathbf{q}}_k^{\prime }\mathbf{b},$$

$${\mathbf{q}}_k={\left(1,\frac{\sqrt{2}\cos (2\pi {\omega }_k)}{\pi },\frac{\sqrt{2}\cos (4\pi {\omega }_k)}{2\pi },\dots ,\frac{\sqrt{2}\cos (2J\pi {\omega }_k)}{J\pi }\right)}^{\prime }$$

$${\Sigma }_b=\left(\begin{array}{cc}{\sigma }_{\alpha }^2& 0\\ 0& {\tau }_b^2{I}_J\end{array}\right),$$

$$\begin{array}{rl}p(\mu \mid \mathbf{x},f)& \propto p(\mathbf{x}\mid f,\mu )p(\mu )\\ & \propto \exp \left\{-\frac{1}{2}\sum _{k=1}^n\frac{I_k}{f({\omega }_k)}\right\}p(\mu )\\ & \propto \exp \left\{-\frac{1}{2}\frac{I_1}{f({\omega }_1)}\right\}\text{mathbbm}1({\mu }_{-}<\mu <{\mu }_{+})\\ & \propto \exp \left\{-\frac{1}{2}\frac{(\bar{x}-\mu {)}^2}{f({\omega }_1)/n}\right\}\text{mathbbm}1({\mu }_{-}<\mu <{\mu }_{+}),\end{array}$$

$$(\mu \mid \mathbf{x},f)\sim N_{{\mu }_{-},{\mu }_{+}}(\bar{x},f({\omega }_1)/n),$$

$$-E_{\mathbf{x}}\left[\frac{{\partial }^2}{\partial \mathbf{b}\partial {\mathbf{b}}^{\prime }}\log p(\mathbf{b}\mid \mathbf{x},\mu ,{\tau }_b^2)\right],$$

Full text

**Key words: Locally stationary time series; Measles; Multiple time series; Rainfall; Reversible jump Markov chain Monte Carlo; Whittle likelihood; **

1 Introduction

When the available data are multiple time series thought to be realizations of nonstationary random processes, estimation of their time varying mean and spectrum offers insight into the behavior of the processes, including whether and how they have changed over time. For example, an analysis of the spatial distribution of the frequency domain characteristics of time series of rainfall for sites spanning a wide spatial field can quantify the cyclical variability of the underlying process, while allowing for nonstationarity can suggest ways in which the climate has changed over the observation period. Extending the work of Priestley (1965), Dahlhaus (1997) established an asymptotic framework for locally stationary processes in which the spectral density of the process is allowed to evolve over time. Key to these asymptotics is the notion that better estimates of the local spectrum for a single stochastic process come from observing the stochastic process at finer time intervals, rather than observing it over a longer time period. This is problematic for historic time series, e.g., rainfall, for which no further observations can be made. Joint modeling of multiple time series with similar or identical local spectra is one way to ameliorate this problem, as estimates for different time series can borrow strength from each other. If additional information is available, such as covariates, it should be incorporated into the model to improve estimation. This is the motivation behind the approach taken in this article.

In particular, this article presents methodology for analyzing a panel of possibly nonstationary time series using a covariate-dependent infinite mixture model, with mixture components parameterized by their time varying mean and spectrum. The mixture components are based on AdaptSPEC (Rosen et al., 2012), which partitions a (centered) time series into an unknown but finite number of segments, estimating the spectral density within each segment by smoothing splines. As part of the proposed method, AdaptSPEC is extended to handle missing values, a common feature of time series which can cause difficulties for nonparametric spectral methods. A second extension is the incorporation of a time varying mean, which avoids having to de-mean (center) the time series as a preliminary step. The covariates, which are assumed to be time-independent, are incorporated via the mixture using the logistic stick breaking process (LSBP) of Rigon and Durante (2020), where the log odds for each ‘stick break’ are modeled using a thin plate spline Gaussian process over the covariates. The model is formulated in a Bayesian framework, where Markov chain Monte Carlo (MCMC) methods are used for parameter estimation and to deal with missing values. Specifically, as in AdaptSPEC, reversible jump MCMC (RJMCMC) is used to estimate the mixture component parameters, while the LSBP parameters are estimated via the Pólya-Gamma based latent variable expansion of Rigon and Durante (2020) (see also Polson et al., 2013, for the original latent variable expansion in the finite mixture case). The model and sampling scheme are capable of handling large panels, such as that of the measles application which has nearly 200,000 observations. In addition to estimating time varying spectra for each time series in the panel, the covariate-dependent mixture structure allows inference about the underlying process at unobserved covariate values, enabling predictive inference. For instance, in this work, we use longitude and latitude as covariates when modeling Australian rainfall data, and are able to infer the predictive time varying spectrum of the rainfall process at unobserved locations.

Many methods have been proposed for the spectral analysis of time series. As the focus of this paper is multiple time series, we provide background in the form of an overview of methods for nonstationary single time series. We then review methods for multiple time series, stationary or otherwise. This excludes methods for multivariate time series, and we refer readers to Li and Krafty (2018) for a review of past and recent work in this active research area.

Approaches to spectral estimation for a single nonstationary time series include fitting parametric time series models with time varying parameters (Kitagawa and Gersch, 1996; Dahlhaus, 1997; West et al., 1999; Yang et al., 2016), smoothing the log periodogram (Ombao et al., 2001; Guo et al., 2003), dividing the time series into locally stationary segments (Adak, 1998; Davis et al., 2006; Rosen et al., 2009, 2012), and using short time Fourier transforms (Yang and Zhou, 2020). For a recent and extensive overview of methods for single nonstationary time series, see Yang et al. (2016). Most directly relevant to this paper, Rosen et al. (2009) estimate the log of the spectral density using a Bayesian mixture of splines. The time series is partitioned into small sections, and it is assumed that the log spectral density within each partition is given by a mixture of smoothing splines. The mixture weights are assumed to be time varying. Rosen et al. (2012) present the AdaptSPEC method, which avoids the fixed partitions of Rosen et al. (2009). AdaptSPEC partitions the time series into one or more variable length segments in an adaptive manner, modeling the log spectral density within each segment via a smoothing spline. This results in better estimates than those obtained from the method of Rosen et al. (2009). Furthermore, the method can accommodate both slowly and abruptly varying processes, as well as identify stationary processes. AdaptSPEC forms the basis of our proposed model for multiple time series.

For multiple time series, Krafty et al. (2011) construct a covariate-dependent model for multiple stationary series in which the log spectrum has a mixed effects representation, where the effects are functions over the frequency domain. Macaro and Prado (2014) propose a Bayesian model for multiple stationary time series with a covariate-dependent spectral density composed as a sum of spectral densities corresponding to different levels of two or more factors. Following Choudhuri et al. (2004), they model the spectral density associated with the factors via Bernstein-Dirichlet priors. Krafty et al. (2017) present a Bayesian model for stationary multivariate time series based on the work of Rosen and Stoffer (2007), where multiple (multivariate) time series from different subjects are available, and subjects have an associated single covariate. Bruce et al. (2018) present a method for multiple nonstationary time series with a single covariate that is referred to as conditional adaptive Bayesian spectrum analysis, which adaptively partitions both time- and covariate-space, modeling the spectrum within each partition by smoothing splines (as in Rosen et al., 2012). Cadonna et al. (2018) model multiple stationary time series using a Bayesian hierarchical model. The log-periodogram of a single stationary time series is modeled as a mixture of Gaussian distributions where the mixture weights and mean functions are frequency-dependent. The hierarchical model for multiple time series is constructed by setting the mean functions to be common to all time series while letting the weights vary between time series.

Adding to this literature, we present methodology, referred to as AdaptSPEC-X, which combines four features: multiple time series, nonstationarity in both mean and spectrum, multiple covariates, and missing data. We demonstrate the method on simulated data, and show how it can be used to estimate the mean and spectra in two application areas: Australian rainfall data, and measles incidence in the United States. Software implementing AdaptSPEC-X is available in the R package BayesSpec111Available from the authors; latest CRAN version does not contain AdaptSPEC-X..

The paper proceeds as follows. Section 2 describes AdaptSPEC, the model for single nonstationary time series forming the basis for the analysis of multiple nonstationary time series. AdaptSPEC-X, a covariate-dependent infinite mixture model, is presented in Section 3. Section 4 outlines the MCMC scheme used to estimate the model parameters. Section 5 presents a study of the performance of AdaptSPEC-X on replicated simulated data from multiple time series. Section 6 describes the application areas. Australian rainfall is analyzed in Section 6.1, and measles incidence in the US is discussed in Section 6.2. Appendix A provides details of the conditional distributions necessary for the sampling scheme, and Appendix B expands on the covariance structure used to derive the conditional distribution of the missing values.

2 Model for single time series

Rosen et al. (2012) (henceforth RWS12) present the AdaptSPEC method for modeling single nonstationary time series which we summarize in this section. Let $\bm{x}=(x_{1},\ldots,x_{n})^{\prime}$ be a time series of length $n$. Assume for ease of notation that $n$ is even, and suppose initially that $\bm{x}$ is a realization from a stationary process $\{X_{t}\}$ with constant mean $\mu$ and a bounded positive spectral density $f(\omega)$ for $\omega\in(-\frac{1}{2},\frac{1}{2}]$. Whittle (1957) showed that, for large $n$, the likelihood of $\bm{x}$ can be approximated as

[TABLE]

where $\omega_{k}=\frac{k-1}{n}$ for $k=1,\ldots,n$ are the Fourier frequencies, $I_{k}=|d_{k}|^{2}$ is the periodogram at $\omega_{k}$, and

[TABLE]

is the discrete Fourier transform (DFT) at $\omega_{k}$, with $i=\sqrt{-1}$. RWS12 follow Wahba (1990) by expressing $\log f$ as

[TABLE]

and placing a smoothing spline prior on $h(\omega)$. Due to the evenness of $f(\omega)$ and the periodogram, $h(\omega)$ is modeled on the domain $\omega\in[0,0.5]$, corresponding to the first $\frac{n}{2}+1$ Fourier frequencies. This prior is expressed via a linear combination of $J$ basis functions, where $J<\frac{n}{2}+1$ is chosen to balance prior flexibility and computational resources. See Appendix A.1 and RWS12 for details.

Next we allow the underlying process $\{X_{t}\}$ to be nonstationary. Let a time series consist of a number of segments, $m$, and let $\xi_{s,m}$ be the end of the $s$th segment, $s=1,\ldots,m$, where $\xi_{0,m}=0$ and $\xi_{m,m}=n$. Then we assume that $\{X_{t}\}$ is piecewise stationary, with

[TABLE]

where the processes $\{X_{t}^{s}\}$ are independent and stationary with means $\mu_{s,m}$, spectral densities $f_{s,m}(\omega)$, and $\delta_{s,m}(t)=1$ iff $t\in(\xi_{s-1,m},\xi_{s,m}]$. Consider a realization $\bm{x}$ from (4). RWS12 approximate the likelihood of $\bm{x}$ by

[TABLE]

where $\Theta=\{m,\bm{\xi}_{m},\bm{\mu}_{m},f_{1,m},\ldots,f_{m,m}\}$, $\bm{\xi}_{m}=(\xi_{1,m},\ldots,\xi_{m,m})^{\prime}$, $\bm{\mu}_{m}=(\mu_{1,m},\ldots,\mu_{m,m})^{\prime}$, $\bm{x}_{s,m}=\{x_{t}:\delta_{s,m}(t)=1\}$ are the data for the $s$th segment, and $p(\bm{x}_{s,m}\mid\mu_{s,m},f_{s,m})$ is the Whittle likelihood (1). The values of $m$, $\bm{\xi}_{m}$ and $f_{s,m}$ for $s=1,\ldots,m$ are considered unknown and are assigned priors. For $f_{s,m}$, the prior in (3) is used, while $m$ is given the discrete uniform prior between $1$ and $M$. (See RWS12 for the details of the prior on $\bm{\xi}_{m}$). RWS12 consider $\bm{\mu}_{m}$ to be known and equal to zero, but in this work we consider it unknown and assign to $\mu_{s,m}$ a uniform prior with support $\mu_{-}<\mu_{s,m}<\mu_{+}$. A minimum segment length $t_{\text{min}}$ is set to ensure that there are sufficient time periods within each segment so that the Whittle likelihood approximation is appropriate.

RWS12 develop a reversible jump Markov chain Monte Carlo algorithm (Green, 1995) that samples from the posterior distribution of this model. They show that AdaptSPEC can handle both abruptly and slowly varying nonstationary time series, as well as identify whether a time series is stationary. AdaptSPEC forms the basis of our spectral estimation technique for multiple time series.

3 Model for multiple time series

In this section we extend the model of Section 2 to multiple time series. Suppose now that the stochastic process $\{X_{t}\}$ has associated covariates $\bm{u}=(u_{1},\ldots,u_{P})^{\prime}$. We model $\{X_{t}\}$ with a covariate-dependent mixture structure

[TABLE]

where the mixture component distributions $g_{h}$ are instances of AdaptSPEC (Equation (5)) with parameters $\Theta_{h}=\{m^{h},\bm{\xi}^{h}_{m},\bm{\mu}^{h}_{m},f^{h}_{1,m},\ldots,f^{h}_{m,m}\}$, $2\leq H\leq\infty$, and the mixture weights $\pi_{h}(\cdot)$ satisfy $0\leq\pi_{h}(\cdot)\leq 1$ and $\sum_{h=1}^{H}\pi_{h}(\bm{u})=1$. Equation (6) implies that $\{X_{t}\}$’s distribution is determined by its covariates $\bm{u}$, which, importantly, do not vary with time. The purpose of the mixture structure in Equation (6) is to induce covariate-dependence in a flexible, semi-parametric manner, and we do not use this structure to perform inference about clustering or the number of clusters among multiple time series.

Let $\{\bm{x}_{1},\ldots,\bm{x}_{N}\}$ be a finite collection of $N$ time series, of length $n$ each, where each time series $\bm{x}_{j}=(x_{1,j},\ldots,x_{n,j})^{\prime}$ has covariates $\bm{u}_{j}=(u_{1,j},\ldots,u_{P,j})^{\prime}$ for $j=1,\ldots,N$. Assuming independence conditional on $\pi_{h}(\cdot)$ and $\Theta_{h}$, it follows from Equation (6) that the joint distribution of the collection is

[TABLE]

3.1 Model for mixture weights

For the mixture weights $\pi_{h}(\bm{u})$ in Equation $\eqref{eqn:single_mixture}$, we use the logit stick-breaking prior (LSBP) developed by Rigon and Durante (2020), according to which $\pi_{h}(\bm{u})$ is given by

[TABLE]

where $\operatorname{logit}v_{h}(\bm{u})=w_{h}(\bm{u})$, so that $w_{h}(\bm{u})$ are the covariate-dependent log odds. This prior allows for $2\leq H\leq\infty$ mixture components, with $v_{H}(\bm{u})=1$ when $H<\infty$.

As described by Rigon and Durante (2020), this construction can be interpreted via sequential (continuation-ratio) logits (Agresti, 2018). Let $z_{j}\in\{1,2,\ldots,H\}$ be a latent indicator such that $(\bm{x}_{j}\mid z_{j}=h)\sim g_{h}(\bm{x}_{j}\mid\Theta_{h})$. Then the LSBP can be represented in a generative manner as a sequence of decisions, where $p(z_{j}=1)=v_{1}(\bm{u}_{j})=(1+\exp(-w_{1}(\bm{u}_{j})))^{-1}$, $p(z_{j}=2\mid z_{j}>1)=v_{2}(\bm{u}_{j})=(1+\exp(-w_{2}(\bm{u}_{j})))^{-1}$, and so on, such that in general, $p(z_{j}=h\mid z_{j}>h-1)=(1+\exp(-w_{h}(\bm{u}_{j})))^{-1}$.

The model in equations (7) and (8) has a similar structure to the classic mixture of experts model (Jacobs et al., 1991), in which the weights of a finite mixture depend on covariates through multinomial logits. Our motivation for choosing the LSBP (Equation (8)) over multinomial logits is to obviate the choice of the number of mixture components. The model in equations (7) and (8) is in principle an infinite mixture and so the question of the number of components becomes irrelevant. In practice, however, it is common to truncate the infinite representation at a suitably high but finite $K$ (Ishwaran and James, 2001). The LSBP is analogous to the probit stick-breaking process (Chung and Dunson, 2009), where a probit link function is used in place of the logit.

3.2 Model for log odds

We model the log odds by a Gaussian process (GP) prior $w_{h}(\bm{u})\sim\text{GP}(\beta_{0,h}+\bm{u}^{\prime}\bm{\beta}_{h},\tau_{h}^{2}\Omega(\cdot,\cdot))$, where $\beta_{0,h}$ is an intercept, $\bm{\beta}_{h}=(\beta_{1,h},\ldots,\beta_{P,h})^{\prime}$ is a vector of regression coefficients, $\tau_{h}^{2}$ is a smoothing parameter, and $\Omega(\bm{u},\bm{u}^{\prime})$ is the covariance kernel constructed via the reproducing kernel Hilbert space defined by a $P$-dimensional thin-plate Gaussian process prior (see Wood, 2013). For a finite collection of $N$ time series $\{\bm{x}_{1},\ldots,\bm{x}_{N}\}$, with associated covariates $\{\bm{u}_{1},\ldots,\bm{u}_{N}\}$, the log odds vector, $\bm{w}_{h}=(w_{h}(\bm{u}_{1}),\ldots,w_{h}(\bm{u}_{N}))^{\prime}$, has a multivariate normal distribution

[TABLE]

where $\bm{1}_{N}$ is an $n\times 1$ vectors of ones, $U=(\bm{u}_{1},\ldots,\bm{u}_{N})^{\prime}$ is an $N\times P$ matrix, and $\Sigma_{w}$ is an $N\times N$ matrix whose $j_{1}j_{2}$th entry is equal to $\Omega(\bm{u}_{j_{1}},\bm{u}_{j_{2}})$. To facilitate the posterior sampling scheme in Section 4, we transform the problem via a basis expansion (Wood, 2013). Let $\Sigma_{w}=QDQ^{\prime}$ be the eigenvalue decomposition of $\Sigma_{w}$, where $Q$ is an $N\times N$ orthogonal matrix whose columns are the eigenvectors of $\Sigma_{w}$, and $D$ is a diagonal matrix containing the eigenvalues of $\Sigma_{w}$. Define $U^{\dagger}=(\bm{1}_{N},U,QD^{1/2})$ by columnwise concatenation, let $\bm{\beta}^{\mathrm{GP}}_{h}$ be an $N\times 1$ vector, and $\bm{\beta}^{\dagger}_{h}=(\beta_{0,h},\bm{\beta}_{h}^{\prime},\bm{\beta}^{\mathrm{GP}\prime}_{h})^{\prime}$. The first column of $U^{\dagger}$ is a vector of ones, the next $P$ columns of $U^{\dagger}$ (equal to $U$) are the original covariates, while the last $N$ columns (equal to $QD^{1/2}$) are basis functions, with associated coefficients $\bm{\beta}^{\mathrm{GP}}_{h}$. For computational convenience and parsimony, we truncate the basis expansion to the first $B<N$ basis functions, so that $U^{\dagger}$ is $N\times(P+B+1)$ and $\bm{\beta}^{\dagger}_{h}$ is $(P+B+1)\times 1$. Equation (9) now takes the form

[TABLE]

where $\bm{0}_{B}$ is an $B\times 1$ vector of zeros, and $I_{B}$ is the $B\times B$ identity matrix. This basis expansion, combined with the interpretation via sequential logits given in the previous section, facilitate the development of the posterior sampling scheme presented in Section 4.

The expression for $v_{h}(\bm{u}_{j})$ becomes $\operatorname{logit}v_{h}(\bm{u}_{j})=\bm{u}^{\dagger\prime}_{j}\bm{\beta}^{\dagger}_{h}$, where $\bm{u}^{\dagger\prime}_{j}$ is the $j$th row of $U^{\dagger}$. The prior placed on $(\beta_{0,h},\bm{\beta}_{h}^{\prime})^{\prime}$ is $\text{N}(\bm{\mu}_{\beta},\Sigma_{\beta})$, where $\bm{\mu}_{\beta}$ is a $(P+1)\times 1$ vector and $\Sigma_{\beta}$ is a $(P+1)\times(P+1)$ covariance matrix, so that

[TABLE]

Finally, to complete the model specification, we assign $\tau_{h}$ a half-$t$ distribution (Gelman, 2006) with density

[TABLE]

where $A_{\tau}$ and $\nu_{\tau}$ are scale and degrees of freedom parameters, respectively. As described in Appendix A, the half-$t$ distribution can be expressed as a scale mixture of inverse Gamma distributions, which simplifies the sampling of $\tau_{h}$. For the results in this paper, we set $\bm{\mu}_{\beta}=\bm{0}$, $\Sigma_{\beta}=100I_{P+1}$, $\nu_{\tau}=3$ and $A_{\tau}=10$.

Figure 1 displays a graphical summary of the model in Equations (5) to (10), showing the dependence between the data, covariates, parameters, and hyperparameters.

3.3 Missing values

The model can accommodate missing values by exploiting the fact that the Whittle likelihood describes a multivariate normal distribution (Whittle, 1953). As is shown below, this implies that the conditional distribution of the missing values is also multivariate normal, and we use this fact to accommodate missing values by integrating them out as part of the MCMC scheme in Section 4. For ease of exposition, in this section we first return to the case where $\bm{x}$ is a single stationary time series with spectral density $f$, then later describe how this is extended to multiple nonstationary time series. Define the $n\times n$ matrix $V$ with entries $V_{tk}=\frac{1}{\sqrt{n}}\exp\left(-2\pi i(t-1)\omega_{k}\right)$ for $t=1,\ldots,n$ and $k=0,\ldots,n-1$. It then follows that $V^{\prime}(\bm{x}-\mu)=(d_{1},\ldots,d_{n})^{\prime}$, where $d_{k}$ is given in (2). Noting that $V$ is a unitary matrix, (1) may be rewritten as

[TABLE]

where $R=\operatorname{diag}(\bm{r})$, $\bm{r}=(1/f(\omega_{1}),\ldots,1/f(\omega_{n}))^{\prime}$, and $V^{*}$ is the conjugate transpose of $V$. The precision matrix $\Lambda=VRV^{*}$ is symmetric and circulant, and is thus defined by its first column, with entries $\Lambda_{t,1}=\frac{1}{n}\sum_{k=1}^{n}\frac{1}{f(\omega_{k})}e^{-2\pi i(t-1)\omega_{k}}$ (see Appendix B for derivation). Suppose some values of $\bm{x}$ are missing, and write $\bm{x}=(\bm{x}_{\text{mis}}^{\prime},\bm{x}_{\text{obs}}^{\prime})^{\prime}$, where $\bm{x}_{\text{mis}}$ and $\bm{x}_{\text{obs}}$ are the missing and observed values, respectively. From standard multivariate normal conditioning results,

[TABLE]

where

[TABLE]

and the quantities in Equation (12) can be obtained from expressing $\Lambda$ as

[TABLE]

in which $V_{\text{mis}}$ and $V_{\text{obs}}$ are matrices made up of the rows of $V$ corresponding to the missing and observed times, respectively. The quantities in Equation (13) can be computed efficiently by the fast Fourier transform. When simulating from Equation (11), the most computationally intensive step is the inversion, $\Lambda_{\text{mis},\text{mis}}^{-1}$, in Equation (12). In a general framework for spectral estimation with stationary time series, Guinness (2019) describes computationally efficient methods for missing data imputation that could be used to simulate from Equation (11). For this article we compute $\Lambda_{\text{mis},\text{mis}}^{-1}$ the usual way, i.e., via its Cholesky decomposition.

Missing values may be accommodated in AdaptSPEC by partitioning the data in each segment, $\bm{x}_{s,m}$, into missing and observed times as above, and sampling from (11) for each segment as part of the MCMC scheme. As shown in the next section, this is extended to the full AdaptSPEC-X model by conditioning on the latent indicators $\bm{z}$.

The Whittle likelihood is used in AdaptSPEC-X in two ways: first, as a nonparametric technique based on its asymptotic properties (see Section 2 and Rosen et al., 2012), and second, in the assumption that the missing values follow a multivariate normal distribution. We justify the latter assumption in several ways. First, Guinness (2019), who makes a similar assumption regarding missing values, found through both theory and numerical experiments that the assumption did not have a deleterious effect on spectral estimation. Second, it is hard to avoid making some assumption about the distribution of the missing values, and multivariate normality is arguably the most parsimonious choice, as it requires only the first two central moments to be specified (the minimum necessary to have mean and spectrum consistent with the observed data). Finally, as shown above and by Guinness (2019), the assumption is computationally convenient.

4 Sampling scheme

Define $\bm{z}=(z_{1},\ldots,z_{N})^{\prime}$, $\bm{x}^{\text{all}}=(\bm{x}_{1}^{\prime},\ldots,\bm{x}_{N}^{\prime})^{\prime}$, $\bm{\beta}^{\dagger}=\{\bm{\beta}^{\dagger}_{1},\ldots,\bm{\beta}^{\dagger}_{H-1}\}$, $\Theta=\{\Theta_{1},\ldots,\Theta_{H}\}$, and $\bm{\tau}=(\tau_{1},\ldots,\tau_{H-1})^{\prime}$. Let $\bm{x}^{\text{all}}=(\bm{x}^{\text{all}\prime}_{\text{mis}}\bm{x}^{\text{all}\prime}_{\text{obs}})^{\prime}$ be the decomposition of $\bm{x}^{\text{all}}$ into missing and observed times, respectively. Values produced in steps 2, 3, 4, and 5 of the MCMC sampling scheme below are indicated by the superscript ‘$[l+0.5]$’, and are then used in a label swapping move in Step 6 to produce the $(l+1)$th iteration. As described below, this move improves convergence of the following sampling scheme.

1. Step 1.

$(\bm{x}_{\text{mis}}^{\text{all}[l+1]}\mid\bm{\beta}^{\dagger[l]},\Theta^{[l]},\bm{z}^{[l]},\bm{x}^{\text{all}}_{\text{obs}})$ as per Section 3.3. 2. Step 2.

$(\Theta^{[l+0.5]}_{h}\mid\bm{x}_{\text{mis}}^{\text{all}[l+1]},\bm{z}^{[l]},\bm{x}^{\text{all}}_{\text{obs}})$ for each $h=1,\ldots,H$. This step is potentially transdimensional, and uses the reversible-jump MCMC scheme of RWS12, with two modifications. The first one samples the segment means, $\bm{\mu}^{h}_{m}$, as described in Appendix A.1.1. The second modification incorporates a Riemann manifold Hamiltonian Monte Carlo (RMHMC for short, see Girolami and Calderhead, 2011) step to accelerate convergence, as described in Appendix A.1.2. 3. Step 3.

$(\bm{z}^{[l+0.5]}\mid\bm{\beta}^{\dagger[l]},\Theta^{[l+0.5]},\bm{x}_{\text{mis}}^{\text{all}[l+1]},\bm{x}^{\text{all}}_{\text{obs}})$. 4. Step 4.

$(\bm{\beta}^{\dagger[l+0.5]}_{h}\mid\bm{\tau}^{[l]},\bm{z}^{[l+0.5]})$ for each $h=1,\ldots,H$. This uses the Polya-Gamma data augmentation scheme developed by Polson et al. (2013), as applied to the LSBP by Rigon and Durante (2020). 5. Step 5.

$(\tau^{[l+0.5]}_{h}\mid\bm{\beta}^{\dagger[l+0.5]})$ for each $h=1,\ldots,H$. 6. Step 6.

$(\Theta^{[l+1]},\bm{z}^{[l+1]},\bm{\beta}^{\dagger[l+1]},\bm{\tau}^{[l+1]}\mid\bm{x}_{\text{mis}}^{\text{all}[l+1]},\bm{x}^{\text{all}}_{\text{obs}})$ using a label swapping step, described below.

The details of steps 2, 3, 4, and 5 are presented in Appendix A, while Step 1, for $\bm{x}_{\text{mis}}$, is described in Section 3.3. For Step 6, we adapt a label swapping move from Hastie et al. (2015), who find that it improves convergence in the context of MCMC samplers for Dirichlet process mixture models. The label swapping step is composed of the following substeps:

1. Step 6a.

Pick uniformly at random components $h_{1},h_{2}\in\{1,\ldots,H\}$, $h_{1}<h_{2}$, to swap. 2. Step 6b.

Construct proposal component indicators $\bm{z}^{\text{swap}}$ and $\Theta^{\text{swap}}$ such that

[TABLE] 3. Step 6c.

Construct proposal $\bm{\tau}^{\text{swap}}$ by setting

[TABLE]

and sample proposal $\bm{\beta}^{\dagger\text{swap}}_{h_{1}},\bm{\beta}^{\dagger\text{swap}}_{h_{2}}$ from

[TABLE]

where $\mu^{\text{mode}}_{h}$ and $\Sigma^{\text{mode}}_{h}$ are the mode and the negative inverse of the Hessian, respectively, of $\log p(\bm{\beta}^{\dagger\text{swap}}_{h}\mid\bm{z}^{\text{swap}},\bm{\tau}^{\text{swap}})$. 4. Step 6d.

Accept the swap with probability equal to the Metropolis-Hastings ratio

[TABLE]

If accepted, set $\bm{\beta}^{\dagger[l+1]},\bm{z}^{[l+1]},\Theta^{[l+1]}$ and $\bm{\tau}^{[l+1]}$ equal to $\bm{\beta}^{\dagger\text{swap}},\bm{z}^{\text{swap}},\Theta^{\text{swap}}$ and $\bm{\tau}^{\text{swap}}$, respectively. Otherwise, $\bm{\beta}^{\dagger[l+1]},\bm{z}^{[l+1]},\Theta^{[l+1]}$ and $\bm{\tau}^{[l+1]}$ are set to $\bm{\beta}^{\dagger[l+0.5]},\bm{z}^{[l+0.5]},\Theta^{[l+0.5]}$ and $\bm{\tau}^{[l+0.5]}$, respectively.

In Step 6b, the labels of the component indicators for the chosen pair $h_{1},h_{2}$ are swapped, as are the corresponding mixture component parameters $\Theta_{h}$, leaving the likelihood unchanged. Step 6c swaps the smoothing spline parameters $\tau_{h_{1}},\tau_{h_{2}}$ and samples new values of $\bm{\beta}^{\dagger}_{h}$ from a normal approximation centered on its conditional mode. The latter new values are necessary because due to the sequential nature of the LSBP, merely swapping the values of $\bm{\beta}^{\dagger}_{h_{1}}$ and $\bm{\beta}^{\dagger}_{h_{2}}$ is unlikely to result in an acceptable proposal.

AdaptSPEC-X is implemented in the latest version of the R package BayesSpec.222Available from the authors; latest CRAN version does not contain AdaptSPEC-X. The implementation is in R and C++, and can take advantage of multiple processor cores to reduce the running time of the analysis.

5 Simulation study

We now demonstrate AdaptSPEC-X using replicated simulated data including covariates and multiple time series with known time varying mean and spectrum. We are interested in the model’s ability to recover the means and spectra at both observed and unobserved covariate values. Let $U=(\bm{u}_{1},\ldots,\bm{u}_{100})^{\prime}$ be a $100\times 2$ design matrix corresponding to $N=100$ subjects, each with two covariates, where the $\bm{u}_{j}$ are sampled uniformly from $[0,1]\times[0,1]$. Each $\bm{u}_{j}$ is mapped deterministically to $z_{j}\in\{1,2,3,4\}$, according to the plot shown in Figure 2(a), which also includes the locations of the $100$ sampled points, denoted by crosses. Four locations are chosen as example time series, marked in green circles, and labeled D1 through to D4 (corresponding to $z_{j}=1$ through $z_{j}=4$, respectively). Four more locations are marked in red diamonds labeled T1 through to T4 (again for $z_{j}=1$ to $4$). These values of $\bm{u}$ have no corresponding time series and are used as test points to evaluate the predictive inferences. The four different regions correspond to four different data generating processes. Each time series $\bm{x}_{j}$, within region $z_{j}$, $j=1,\ldots,100$, is a realization of length $n=256$ from

[TABLE]

where $\epsilon_{j,t}\sim\text{N}(0,1)$, and the values of $\mu_{z_{j},t}$ and $\phi_{z_{j},p,t}$ are given in the following table:

[TABLE]

Thus, time series with $z_{j}=1$ have two segments with different means and different spectra, those with $z_{j}=2$ have two segments with different means but with the same spectra, time series with $z_{j}=3$ have only one stationary segment, and those with $z_{j}=4$ have two segments with the same mean but with different spectra. In each time series, 10% of the times are set as missing. Figure 2(b) displays example realizations from Process (14) for $z_{j}=1,2,3$ and $4$, showing the time series values, underlying time varying mean, and the times at which values are missing.

We sample 100 replicates from Process $\eqref{eqn:multiple_simulation_study_model}$, and to each fit AdaptSPEC-X using the MCMC sampling scheme of Section 4. We run 50,000 iterations of the MCMC scheme, where the first 10,000 are discarded as burn-in. Each mixture component has $M=4$, as the maximum number of segments, $t_{\text{min}}=40$, as the minimum segment length, $J=25$, as the number of basis functions for the smoothing spline prior on the log spectra, and $(\mu_{-},\mu_{+})=(-10,10)$, as the support of the prior on $\mu^{h}_{s,m}$. The LSBP is truncated at $H=25$ components, and has $B=10$ basis functions. To assess the quality of the estimated time varying mean, we define the mean squared error (MSE) for the mean as

[TABLE]

where $\hat{\mu}(t,\bm{u})$ is the estimate of $\mu(t,\bm{u})$, the true time varying mean at covariates $\bm{u}$. Similarly, we define the MSE for the spectrum as

[TABLE]

where $k_{\text{max}}=128$, and $\log\hat{f}(t,\omega,\bm{u})$ is the estimate of $\log f(t,\omega,\bm{u})$, the true time varying log spectral density at location $\bm{u}$. Figure 3 presents boxplots of $\text{MSE}_{\text{mean}}$ (top) and $\text{MSE}_{\text{spec}}$ (bottom), at each observed location, D1–D4, and unobserved test location, T1–T4, from left to right, respectively. The median $\text{MSE}_{\text{mean}}$ is less than 0.02 at all covariate values except for D2 and T2, for which it is 0.09. Similarly, the median $\text{MSE}_{\text{spec}}$ is less than 0.08 at all locations except for D2 and T2, for which it is 0.34 and 0.32, respectively. Estimates of the time varying mean and spectrum corresponding to the median MSE values are shown in figures 4 and 5, respectively. These qualitatively match the $\text{MSE}_{\text{mean}}$ and $\text{MSE}_{\text{spec}}$ scores, in that the estimates for points other than D2 and T2 are visually very close to the truth, while for D2 and T2 some differences are visible.

There are several explanations for the worse performance for D2 and T2 relative to the other points. One reason is that the spatial model has relatively little information to identify the existence of a cluster: these points belong to the category $h=2$, which is represented by only 8 time series (in contrast to $h=1,3$ and $4$, which have 41, 18, and 33 members, respectively), and which has the smallest spatial area in the study. Another reason is that the process mean for $h=2$ (equal to $1$ in the first half and $-1$ in the second half) is similar to that of the surrounding category, $h=4$ (constant mean of $1$), making it harder to distinguish between the categories. For both reasons, it is not surprising to see worse performance for T2 and D2, which corresponds to good model behavior in the sense that it represents genuine model uncertainty. This is seen in Figure 4, where the $\hat{\mu}(t,\bm{u})$ for T2 and D2 in the second half of the time series is shrunk towards that of $h=4$, the enclosing category. It is also worth noting that the median $\text{MSE}_{\text{mean}}$ and $\text{MSE}_{\text{spec}}$ for T2 and D2 are small relative to the scales of their mean and log spectrum, respectively. This can also be seen qualitatively by the similarity between the true and estimated mean and spectra in figures 4 and 5.

6 Applications

In this section, we describe two applications. Section 6.1 considers Australian rainfall, while incidence counts of measles in the United States are analyzed in Section 6.2.

6.1 Australian rainfall data

Rainfall is governed in large part by cyclical processes, in particular the seasonal cycle driven by the Earth’s orbit around the sun. It has therefore historically been a natural application area for spectral methods. These have been used to study interannual variation (see, among many others, Alter, 1924; Rajagopalan and Lall, 1998; Ansell et al., 2000), intraannual or intraseasonal variation (Joshi and Pandey, 2011), and the connections between rainfall and other climatic processes (Rajagopalan and Lall, 1998; Ansell et al., 2000). Here we focus on identifying changes in both the mean and spectrum of Australian rainfall. As part of a report on climate change tendered by several Australian government agencies, Cai et al. (2007) found that, since 1950, the Australian north has seen increased annual rainfall, while the southeast and southwest have experienced the opposite. The causes of these trends have been the subject of study and debate (see, among others, Hope et al., 2006; Ummenhofer et al., 2009; Pook et al., 2012; Risbey et al., 2013). Apart from trends in overall rainfall, several authors have reported relative increases in heavy rainfall events, indicating changes in the variability of rainfall (Cai et al., 2007; Gallant et al., 2013). We contribute to this literature by using AdaptSPEC-X to analyze the time varying mean and spectrum of Australian rainfall from sites dispersed over a wide spatial field, addressing simultaneously the question of whether changes have occurred in rainfall levels and rainfall variability.

We use data from Bertolacci et al. (2019), who studied the climatology of Australian daily rainfall using measurements from 17,606 sites across the continent. In particular, we use 151 of these sites characterized by having long and nearly contiguous rainfall records, the locations of which are displayed in Figure 6(a). These sites are among those identified by Lavery et al. (1992) as having high quality records suitable for monitoring and assessing climate change. The raw time series are daily, and observations are typically made at 9 am local time, recording the total rainfall in millimeters (mm) for the previous 24 hours. Aggregation to monthly data is performed by calculating the average daily rainfall for the month. To avoid artifacts, we consider as missing any month with fewer than fifteen days of measurements available (that is, not missing).

The resulting time series span the 1,078 months from September, 1914 to June, 2004 (inclusive), for a total of 162,778 observations, of which 4,095 are missing. The smallest possible measurement is 0 mm, corresponding to no rainfall for the month; this is true for 9,933 months. Time series for four example sites are displayed in Figure 6(b), and their locations are marked on inset maps (also marked in green in Figure 6(a)). The four time series span a wide range of average rainfall levels from around 1mm at site 47053 to 4mm at site 14042. They also exhibit varying levels of seasonality, where sites 10525 and 14042 have highly seasonal rainfall, while rainfall at sites 47053 and 69018 is less seasonal.

We fit AdaptSPEC-X to these data, setting $\bm{u}_{j}=(\text{lon}_{j},\text{lat}_{j})$, the longitude and latitude of the sites. Each mixture component has $t_{\text{min}}=60$ (5 years), chosen to ensure several observations of the dominant annual cycle are available. This constrains the maximum number of segments to be $M_{\mathrm{max}}=17$. For the number of basis functions for the smoothing spline prior on the log spectra, we choose $J=60$, giving sufficient flexibility to represent the spike at the dominant annual frequency. The support of the prior on $\mu^{h}_{s,m}$ is set to $(\mu_{-},\mu_{+})=(0,30)$. The lower bound reflects the positivity of rainfall, while the upper bound is three times larger than the largest empirical mean attained in any 5-year period in the data. The LSBP is truncated at $H=25$ components, which was found to be large enough that higher values made no difference. Finally, the thin-plate GP prior for the LSBP has $B=20$ basis functions, which captures more than 95% of the variation implied by the prior. In addition to the locations with measurements, we estimate the predictive time varying mean and spectrum at four locations on the Australian landmass without measurements. These locations are indicated by red diamonds in Figure 6(a).

Two major droughts occurred during the study period: the World War II drought of 1937–1945, and the Millennium drought that started in 1996 and was still ongoing by the end of the study period in 2004 (Ummenhofer et al., 2009). Table 1 presents estimated posterior probabilities of changes in the mean $\mu(t,\bm{u})$ or variance $\sigma^{2}(t,\bm{u})=2\int_{0}^{1/2}f(t,\omega,\bm{u})d\omega$ around these times. Specifically, it shows that $\hat{P}(\mu_{1940}<\mu_{1950})>0.9$ and $\hat{P}(\sigma^{2}_{1940}<\sigma^{2}_{1950})>0.9$ at all four sites. However, $\hat{P}(\mu_{1940}<\mu_{1930})<0.7$ and $\hat{P}(\sigma^{2}_{1940}<\sigma^{2}_{1930})<0.57$ except for site 69018 for which these probabilities are greater than 0.8. The Millennium drought is associated with drops in $\mu(t,\bm{u})$, $\sigma^{2}(t,\bm{u})$ or both at sites 14042, 47053 and 69018 as can be seen from the estimated probabilities at these sites: $\hat{P}(\mu_{2004}<\mu_{1990})>0.9$ and $\hat{P}(\sigma^{2}_{2004}<\sigma^{2}_{1990})>0.93$. Site 10525 in the southwest of the continent does not exhibit a drop with probability greater than $0.9$, consistent with the fact that the drought principally affected southeastern Australia (Ummenhofer et al., 2009).

Cai et al. (2007) report large trends in rainfall since 1950. For southeast Australia, they report reductions in annual rainfall corresponding to 10–15mm/decade for the site 47053 and 50mm/decade for 69018. Consistent with this, Table 1 shows that, for these sites, AdaptSPEC-X estimates that both $\mu(t,\bm{u})$ and $\sigma^{2}(t,\bm{u})$ declined between January, 1950 and January, 2004 with probability greater than $0.94$. The estimated drop in $\mu(t,\bm{u})$ for site 47053 corresponds to a reduction of 1–16mm/decade333Calculated as $10\times 365.25\times\text{difference in daily average}\ /\ (2004-1950)$ (10th–90th percentile), consistent with Cai et al.’s estimate. However, for site 69018, the estimated reduction is 6–42mm/decade, substantially less than 50mm/decade. For the site 14042 in the tropical north, Cai et al. estimate increases of around 40mm/decade since 1950, while the estimates of Table 1 indicate a decline over the same period. Finally, Cai et al. report a decline of 20–30mm/decade in the region of southwest Australia containing the site 10525, but Table 1 does not indicate a significant change in $\mu(t,\bm{u})$ at this location (though the variance does appear to have declined). In contrast to Cai et al. and Gallant et al. (2013), the reduction in $\sigma^{2}(t,\bm{u})$ since 1950 at all locations suggests that rainfall variability has declined. This could have resulted from the use of different definitions of variability, e.g., counts of extreme events as in Cai et al. (2007), versus our use of change in variance.

6.2 Measles incidence in the United States

Measles is a highly contagious disease that causes fever, cough, runny nose, and a rash (Moreno, 2018). Complications of measles can include pneumonia, deafness, or death (Moreno, 2018). Prior to the licensing of a vaccine in 1963, the incidence rate of measles averaged 318 cases per 100,000 population per year, with outbreaks occurring annually or every other year (van Panhuis et al., 2013). After vaccine licensure, incidence declined dramatically, and endemic measles transmission was declared eliminated from the United States in 2000 (Katz and Hinman, 2004). As part of a study on the impact of vaccination for a variety of contagious diseases, van Panhuis et al. (2013) collated a unique data set by digitizing weekly surveillance reports from the United States of several nationally notifiable diseases, including measles, and have made these data available online at the Project Tycho website444https://www.tycho.pitt.edu/. In this section we analyze Project Tycho’s measles data using AdaptSPEC-X.

The data comprise weekly time series of measles incidence for each state, reporting the weekly (where the week starts on a Sunday) incidence rate per 100,000 population. In this work we use time series from the continental United States (that is, excluding Hawaii and Alaska), plus the District of Columbia. These span the 3,914 weeks from week one of 1928 (which we write as 1928-01) to week one of 2003 (2003-01). Across all 49 time series there are 191,786 observations. Of these, 50,067 (26%) are missing, and 30,439 (16%) have incidence equal to zero. The time series are shown in Figure 10, where each panel presents the series for a state, and the incidence axis is on a square-root scale (note however that the spectral analysis performed later is applied to the untransformed data). The layout of the panels roughly matches the geographic distribution of the states. The most striking aspect of these plots is the dramatic decline in both the level and volatility of incidence starting in 1963, the year of vaccine licensure.

We set $\bm{u}_{j}=(\text{lon}_{j},\text{lat}_{j})$, the longitude and latitude of the centroid of each state, and fit AdaptSPEC-X to the measles time series. Each mixture component has $t_{\text{min}}=208$ (4 years) as the minimum segment length. This was chosen to ensure four observations of the annual cycle in each segment. The maximum number of segments was set as $M_{\mathrm{max}}=18$, the maximum allowed by the combination of $t_{\text{min}}$ and the number of weeks in the data. The support of the prior on $\mu^{h}_{s,m}$ was set to $(\mu_{-},\mu_{+})=(0,20)$; the lower bound represents the known positivity of the incidence rates (and is particularly helpful to constrain $\mu^{h}_{s,m}$ in periods with very small counts), while the upper bound is twice as large as the posterior mean of this parameter. As in the rainfall application, the number of basis functions for the smoothing spline prior on the log spectra was set to $J=60$. This was found to be high enough to represent the spectrum accurately, particularly the spike corresponding to the annual cycle. The LSBP is truncated at $H=10$ components, as we found higher values did not change the results. Finally, the thin-plate GP prior for the LSBP has $B=20$ basis functions, which captures more than 95% of the variation implied by the prior. The resulting estimated time varying means and spectra are shown in figures 11 and 12, respectively. The top two rows of Figure 13 display zoomed-in plots of the time varying spectra for four states: Arizona, Florida, Maine, and Washington.

As in Figure 10, the post-vaccine drop in the mean and power of incidence is the most obvious feature. This drop occurs in steps, starting with a dramatic drop following the licensing of the Edmonston vaccine in 1963, stalling around 1970, then dropping again around 1980. This corresponds to the waxing and waning of government funding and effort targeted at measles elimination, which culminated in an intensified elimination drive (Hinman et al., 1979; Atkinson et al., 1992). An outbreak during the early 1990s is visible as an increase in power in Figure 12; this outbreak received much attention and resulted in changes to the immunization schedule for children (Atkinson et al., 1992).

In the pre-vaccine period, the spectra in Figure 12 have peaks around frequencies $1/52$ and 0, indicating annual seasonality and long-term dependence, respectively. After the introduction of the vaccine, the annual peak disappears. Grenfell et al. (2001) identified a biennial cycle in similar measles data for the UK, but this does not appear to be a feature of the US data. Figures 11 and 12 also indicate changes in the mean and spectrum during the pre-vaccine years, where all states exhibit periods of increased mean incidence and power centered around 1940 and 1955. This concords with van Panhuis et al. (2013), who noted that incidence rates had variable patterns in the pre-vaccine time period, speculating that these may have been due to sanitation, hygiene, or demographic factors.

Because of the extreme nonstationarity introduced by the vaccine, the time varying spectra in Figure 12 span such a wide range of powers that the spectra for all states look almost identical. This is not the case, as shown in Figure 14 and in the bottom two rows of Figure 13, which display time varying spectra for the pre-vaccine years only. These spectra highlight the existence of geographic heterogeneity between states, where higher power is more typical of the west and north, compared to the south and east.

Since the elimination of endemic measles in the US in 2000, there have been a number of outbreaks associated with individuals ‘importing’ measles by acquiring the disease while outside the US and spreading it upon their return (Parker et al., 2006; CDC, 2019). Phadke et al. (2016) associated several of these outbreaks with individuals unvaccinated for nonmedical reasons, which they term as vaccine refusal. Some authors have even declared that a resurgence of measles has occurred (Lynfield and Daum, 2014). Using the AdaptSPEC-X fits, we tested whether the mean $\mu(t,\bm{u})$ or variance $\sigma^{2}(t,\bm{u})$ of measles increased from 1995-01, a few years after the big outbreak in the early 1990s, to 2003-01, the last time period in the data. We find no evidence of increase in $\mu(t,\bm{u})$ in any state, which can be seen from the fact that the highest posterior probability of an increase equals 0.68 in Oklahoma. As for $\sigma^{2}(t,\bm{u})$, the posterior probabilities of an increase in all states range between 0.81–0.87, which we consider to be weak evidence of change. Unfortunately, because the data end in 2003, it is not possible to assess changes to the mean or spectrum of measles incidence in more recent years.

7 Discussion

This article has presented AdaptSPEC-X, a Bayesian method for analyzing a panel of possibly nonstationary time series using a covariate-dependent infinite mixture model, with mixture components parameterized by their time varying mean and spectrum. AdaptSPEC-X extends AdaptSPEC to accommodate multiple time series, each with its own covariate values. Specifically, the covariates, which are assumed to be time-independent, are incorporated via the mixture weights using the logistic stick breaking process. The mixture components are based on AdaptSPEC, which handles a single nonstationary time series. In particular, it partitions a time series into an unknown but finite number of segments, and estimates the spectral density within each segment by smoothing splines. New features which have been added to the AdaptSPEC components include estimation of time varying means and handling of missing observations. The model and sampling scheme can accommodate large panels, such as that of the measles application. In addition to estimating time varying spectra for each time series in the panel, AdaptSPEC-X allows inference about the underlying process at unobserved covariate values, enabling predictive inference. Efficient software implementing AdaptSPEC-X is available in the R package BayesSpec.

In Section 3.2, the log odds of the LSBP, which determine the mixture weights, are modeled using a thin-plate GP prior. While this prior is flexible, it is also smooth and stationary. This property may be inappropriate in settings where changes in the mean or spectrum of the individual time series occur abruptly over the covariate space. An extension to a nonstationary prior for the log odds, or a piecewise prior as in Bruce et al. (2018), may be of interest in these cases.

AdaptSPEC (and therefore AdaptSPEC-X) relies on Whittle’s approximation to the likelihood (Equation (1)). The Whittle likelihood is asymptotically correct for both Gaussian and non-Gaussian time series (Hannan, 1973), but is known to be inefficient for small sample sizes (Contreras-Cristán et al., 2006). Several methods exist in the literature to ameliorate this problem (see, for example, Sykulski et al., 2019), and these methods might produce useful extensions to AdaptSPEC for settings with short time series or small segment lengths (i.e., small $t_{\text{min}}$).

Neither AdaptSPEC nor AdaptSPEC-X account explicitly for measurement error. For i.i.d. Gaussian measurement error, this should not cause a problem, as the added variance would appear as an added constant (that is, white noise) in the spectrum. On the other hand, if the measurement error is changing over time, this may be estimated as false nonstationary, in the sense that the underlying process is not actually changing. Adding an explicit layer for measurement error to AdaptSPEC-X would be a useful extension for these cases. More specialized forms of measurement error may be useful in other settings. For example, in the measles application of Section 6.2, the data are incidence rates per 100,000 population. While rates are continuous quantities, they are constructed from counts, so the measurement process is actually discrete. A hierarchical extension accounting for this would be an interesting addition to AdaptSPEC-X.

AdaptSPEC-X allows for covariates that do not vary with time, and an extension to a more general framework accommodating time varying covariates would facilitate new types of inference. For example, external time varying climate indicators such as the Southern Oscillation Index are known to influence rainfall patterns (Bertolacci et al., 2019). A time varying covariate influencing the mean or the spectrum might introduce useful shrinkage, improving predictive performance. One challenge with such an approach would be to determine how a time varying covariate could interact with the segmentation approach used by AdaptSPEC to handle nonstationarity. It would also be interesting to allow the mean to be time varying and covariate-dependent within segments, not only between segments. Future research will focus on these extensions.

Reproducibility

Data and code reproducing the figures and tables in this manuscript are available online at https://github.com/mbertolacci/adaptspecx.

Funding

E. Cripps and M. Bertolacci were supported by the ARC Industrial Transformation Research Hub for Offshore Floating Facilities which is funded by the Australian Research Council, Woodside Energy, Shell, Bureau Veritas and Lloyds Register (Grant No. 140100012). S. Cripps is the recipient of an Australian Research Council Australian Future Fellowship (140101266) funded by the Australian Government. O. Rosen was supported in part by grants NSF DMS-1512188 and NIH 2R01GM113243-05.

Appendix A Conditional distributions used for the MCMC scheme

This appendix provides the details of the conditional distributions sampled in the MCMC scheme of Section 4, with the exception of the distribution of $(\bm{x}^{\text{all}}_{\text{mis}}\mid\bm{\beta}^{\dagger},\Theta,\bm{z},\bm{x}^{\text{all}}_{\text{obs}})$, which is described in Section 3.3.

A.1 Drawing $\Theta_{h}$

The conditional distribution is

[TABLE]

This step is potentially transdimensional as the number of segments $m$ may increase or decrease, so that $\Theta_{h}$ may change dimension. RWS12 describe a reversible-jump transition kernel for $\Theta_{h}$; the interested reader can find the details therein.

The next two sections describe modifications to the transition kernel of RWS12. The first modification samples the segment means, $\bm{\mu}^{h}_{m}$, and is described in Section A.1.1. The second modification is the incorporation of an RMHMC step to accelerate convergence, described in Section A.1.2. Some additional details from RWS12 are needed. For ease of exposition, in what follows we return to the case of $\bm{x}$ stationary, omitting the component indicator superscript $h$ and the segment subscript $s,m$. RWS12 express the smoothing spline prior on $\log f(\omega)$ (see Section 2) using the linear basis expansion,

[TABLE]

where

[TABLE]

and $\bm{b}=(\alpha_{0},b_{1},\ldots,b_{J})^{\prime}$. The prior on $\bm{b}$ is $(\bm{b}\mid\tau_{b}^{2})\sim N(\bm{0},\Sigma_{b})$, where

[TABLE]

$\sigma^{2}_{\alpha}$ is a hyperparameter, and $\tau_{b}^{2}$ is a smoothing parameter (having its own prior, see RWS12). Note that the symmetry of the spectral density is such that $f(\omega_{k})=f(\omega_{n-k+2})$. Thus, later in this section when we write $\bm{q}_{k}$ for $k>n/2+1$, we refer to $\bm{q}_{n-k+2}$. In general, each segment $(s,m)$ of each component $h$ has an associated $\bm{b}^{h}_{s,m}$ and $\tau^{2h}_{b,s,m}$.

A.1.1 Drawing $\mu$

The conditional distribution of $\mu$ is

[TABLE]

where $\bar{x}$ is the arithmetic mean of $\bm{x}$, and $\operatorname{\mathbbm{1}}(\cdot)$ is an indicator function. The third line of Equation (A.1) follows from the fact that only $I_{1}$ depends on $\mu$, and the final line follows from $I_{1}=n(\bar{x}-\mu)^{2}$. Therefore,

[TABLE]

where $N_{\mu_{-},\mu_{+}}$ denotes the truncated normal distribution on $(\mu_{-},\mu_{+})$.

This distribution is incorporated in two places in the transition kernel of RWS12. First, the proposal distribution for $\bm{b}$ based on the mode of its full conditional used both for the between- and within- model steps (see Appendix of RWS12) is adjusted to propose from the mode of $(\mu,\bm{b}^{\prime})^{\prime}$, where the modal value for $\mu$ is $\bar{x}$. Second, an additional step is added to each iteration of the MCMC scheme, in which a segment $s\in\{1,\ldots,m\}$ is chosen uniformly at random, and the corresponding $\mu_{s,m}$ is updated based on Equation (A.2). This latter step is simply to aid convergence of the resulting Markov chain, and is necessary to ensure ergodicity.

A.1.2 Drawing $\bm{b}$ using RMHMC

Girolami and Calderhead (2011) propose a variant of the Hamiltonian Monte Carlo (HMC) sampling method (Duane et al., 1987; Neal, 1993) for sampling from an unnormalized probability density: Riemann manifold HMC (RMHMC). RMHMC exploits the Riemann geometry of the parameter space of the density in order to efficiently traverse the space and thereby provide an efficient proposal mechanism without the need for tuning. These authors propose the use of the Fisher information plus the negative Hessian of the log prior to define a metric tensor over the parameter space, and show how to modify plain HMC to incorporate this geometry. They show that the performance of RMHMC is excellent, comparable or better than that of competing methods, while requiring less tuning. In this section, we show how an RMHMC step for the smoothing spline parameters $\bm{b}$ may be incorporated into the sampling scheme of RWS12, improving the convergence of the scheme.

The appropriate metric tensor for the conditional distribution of $\bm{b}$ is

[TABLE]

where $E_{\bm{x}}$ means that the expectation is with respect to $\bm{x}$. We now derive this expectation analytically. The conditional distribution for $\bm{b}$ is

[TABLE]

The first and second derivatives of $\log p(\bm{b}\mid\bm{x},\mu,\tau^{2}_{b})$ are

[TABLE]

respectively. The negative expectation of the last line of Equation (A.3) with respect to $\bm{x}$ is

[TABLE]

where the last step follows from the fact that the Whittle likelihood (Equation (1)) implies that $I_{k}\sim f(\omega_{k})\chi^{2}_{1}$ for $k=1$ and $k=\frac{n}{2}+1$, and $I_{k}\sim f(\omega_{k})\chi^{2}_{2}/2$ otherwise, so that $E_{\bm{x}}[I_{k}/e^{\bm{q}_{k}^{\prime}\bm{b}}]=E_{\bm{x}}[I_{k}/f(\omega_{k})]=1$.

Equation (A.4) shows that the metric tensor, $-E_{\bm{x}}\left[\frac{\partial^{2}}{\partial\bm{b}\partial\bm{b}^{\prime}}\log p(\bm{b}\mid\bm{x},\mu,\tau^{2}_{b})\right]$, is constant with respect to $\bm{b}$. A constant metric tensor corresponds to a special case of RMHMC, avoiding computational complexity required in general. This special case is equivalent to plain HMC with a mass matrix equal to the metric tensor (see Neal, 2011, for an explanation of the term mass matrix, and the terms leap frog step size and number of steps used below). We add an RMHMC step to each iteration of the sampling scheme of RWS12, picking a segment $s\in\{1,\ldots,m\}$ uniformly at random, and updating $\bm{b}_{s,m}$ using RMHMC. The HMC leap frog step size and number of steps used are chosen uniformly at random from $[0.1,1]$ and $[1,10]$, respectively.

A.2 Drawing $\bm{z}$

The mixture component indicators have the conditional probability mass function

[TABLE]

which can be sampled directly for each $j=1,2,\ldots,N$.

A.3 Drawing $\bm{\beta}^{\dagger}_{h}$

To sample $\bm{\beta}^{\dagger}_{h}$ we draw on the work of Rigon and Durante (2020), who extend the data augmentation scheme of Polson et al. (2013). Let $Z_{h}=\{j:j\in\{1,2,\ldots,N\},z_{j}\geq h\}$, and $U^{\dagger}_{h}$ be the matrix with rows given by $\bm{u}^{\dagger\prime}_{j}$ for $j\in Z_{h}$. The model is augmented with variables $\eta_{j,h}$ so that

[TABLE]

where PG denotes the Pólya-Gamma distribution (see Polson et al., 2013),

[TABLE]

$\bm{\eta}_{h}=(\eta_{j,h})_{j\in Z_{h}}^{\prime}$, $\kappa_{j,h}=\operatorname{\mathbbm{1}}(z_{j}=h)-0.5$, $\operatorname{\mathbbm{1}}(z_{j}=h)$ is 1 if $z_{j}=h$ and 0 otherwise, $\bm{\kappa}_{h}=(\kappa_{j,h})_{j\in Z_{h}}^{\prime}$, $\bm{\mu}^{\dagger}=(\bm{\mu}_{\beta}^{\prime},\bm{0}_{B}^{\prime})^{\prime}$, and $\Sigma_{\beta^{\dagger}}=\begin{pmatrix}\Sigma_{\beta}&0\\ 0&\tau^{2}_{h}I_{B}\end{pmatrix}$. Sampling then proceeds for each $k\in\{1,2,\ldots,H\}$ by sampling $\eta_{j,h}$ for each $j\in Z_{h}$, then sampling $\bm{\beta}^{\dagger}_{h}$.

A.4 Drawing $\tau_{h}$

We follow Wand et al. (2012), expanding the half-$t$ prior by augmenting the model with a latent variable $a_{h}$, using the hierarchical structure

[TABLE]

so that the full conditional distributions are

[TABLE]

The sampling scheme therefore proceeds by first sampling from $(a_{h}^{[l+1]}\mid\tau_{h}^{2[l]})$, then $(\tau_{h}^{2[l+1]}\mid a_{h}^{[l+1]},\bm{\beta}^{\mathrm{GP}[l+1]})$.

Appendix B The structure of $\Lambda$

This appendix presents the derivation for the structure of the matrix $\Lambda$ of Section 3.3. The matrix $\Lambda$ has entries

[TABLE]

where the notation in the final line expresses the fact that $\Lambda_{t_{1},t_{2}}$ depends only on $t_{1}-t_{2}$. The matrix $\Lambda$ is therefore Toeplitz. Letting $\bm{\lambda}=(\lambda_{0},\lambda_{1},\ldots,\lambda_{n-1})^{\prime}$, we have $\bm{\lambda}=\frac{1}{\sqrt{n}}V^{\prime}\bm{r}$, that is, $\bm{\lambda}$ is the DFT of $\bm{r}$ divided by $\sqrt{n}$. In fact, $\bm{\lambda}$ is real valued, as

[TABLE]

and the fact that $f(\omega_{k})=f(\omega_{n-k+2})$ and $-\sin(2\pi t\omega_{k})=\sin(2\pi t\omega_{n-k+2})$ ensure that the complex parts of the sum cancel. Finally, $\lambda_{t}=\lambda_{n-t}$, so that $\Lambda$ has a symmetric circulant structure.