Bayesian Model Search for Nonstationary Periodic Time Series

Abstract

We propose a novel Bayesian methodology for analyzing nonstationary time series that exhibit oscillatory behavior. We approximate the time series using a piecewise oscillatory model with unknown periodicities, where our goal is to estimate the change-points while simultaneously identifying the potentially changing periodicities in the data. Our proposed methodology is based on a trans-dimensional Markov chain Monte Carlo algorithm that simultaneously updates the change-points and the periodicities relevant to any segment between them. We show that the proposed methodology successfully identifies time changing oscillatory behavior in two applications which are relevant to e-Health and sleep research, namely the occurrence of ultradian oscillations in human skin temperature during the time of night rest, and the detection of instances of sleep apnea in plethysmographic respiratory traces. Supplementary materials for this article are available online.

KEYWORDS:

• Bayesian spectral analysis
• Change-points
• Reversible-jump MCMC
• Sleep apnea
• Ultradian sleep cycles

1 Introduction

Identifying the periodicities present in a cyclical phenomenon allows us to gain insight into the sources of variability that drive the phenomenon. For example, respiratory traces obtained from a plethysmograph used on rodents in experimental sleep apnea research exhibit many abrupt changes in their periodic components as the rat naturally changes their breathing pattern in the course of its sleep-wake activities (Han et al. 2002; Nakamura, Fukuda, and Kuwaki 2003). Similarly, human temperature, as measured by a wearable sensing device over several days at relatively high temporal resolution, may be subject to a different periodic behavior during the night when the individual transitions between ultradian sleep stages (Carskadon and Dement 2005; Komarzynski et al. 2018). While the theory and methods for analyzing the periodicities of time series data whenever the time series is stationary are relatively well-developed, the task of modelling time series that show regime shifts in periodicity, amplitude, and phase remain challenging because the timing of changes and the relevant periodicities are usually unknown.

There have been many developments for modeling stationary oscillatory time series data. Rife and Boorstyn (1976) and Stoica et al. (1989) addressed the problem of estimating the frequencies, phases, and amplitudes of sinusoidal signals under the assumption of a known number of sinusoids, where inference is based on maximum likelihood frequency estimators. These models, however, require very long time series and a large separation in the frequencies that drive the process, which will not always be the case in practice (Djuric 1996; Andrieu and Doucet 1999). Quinn (1989), Yau and Bresler (1993), and Zhang and Wong (1993) tackled the problem of model selection on the number of sinusoidal signals by employing the Akaike information criterion (Akaike 1974) and the minimum description length principle (Rissanen 1978). Djuric (1996) showed that these procedures tend to estimate a wrong number of components when the sample size is small and the signal-to-noise ratio is low.

Bayesian approaches to modeling stationary oscillatory signals were explored for the first time by Bretthorst (1988, 1990) with applications to nuclear magnetic resonance spectroscopy. Dou and Hodgson (1995, 1996) presented a Bayesian approach that uses a Gibbs sampler to identify multiple frequencies that drive the signal. Their method required the number of frequencies to be fixed in advance, and model selection was achieved by choosing the most probable model based on the estimation of the parameters for all possible models. Bayesian model selection for stationary oscillatory signals based on posterior model probabilities were also investigated by Djuric (1996). Andrieu and Doucet (1999) introduced a more efficient reversible-jump Markov chain Monte Carlo (MCMC) method (Green 1995) that jointly tackles model selection and parameter estimation for an unknown number of stationary sinusoidal signals and avoids the computationally expensive numerical optimization of Dou and Hodgson (1995, 1996) by sampling the frequencies one-at-time via Metropolis–Hastings (M-H) steps. To the best of our knowledge, currently there is no extension of this methodology to analyze nonstationary oscillatory signals.

A formal statistical modeling framework for a specific class of nonstationary time series data, called locally stationary time series, was developed by Dahlhaus (1997). Extending this framework to a Bayesian setting, Rosen, Stoffer, and Wood (2009) proposed an approach to model the log of the time-varying spectral density using a mixture of smoothing splines. Rosen, Wood, and Stoffer (2012) improved on this by splitting the time series into an unknown but finite number of segments of variable lengths, thereby avoiding the need to preselect partitions, and to estimate the time-varying spectral density using a fixed number of smoothing splines. For a given partition of the time series, the likelihood function is approximated via a product of local Whittle likelihoods (Whittle 1957). The methodology was developed using a Bayesian framework and is based on the assumption that, conditional on the position and number of partitions, the time series are piecewise stationary, and the underlying spectral density for each partition is smooth over frequencies. However, exploratory analyses of the time series in both of our case studies revealed spectral densities with very sharp peaks, often at several nearby frequencies, thus invalidating the assumption that the spectral density is smooth over frequencies. In addition, the frequency location of these sharp peaks changed over time.

In this article, we propose a novel Bayesian methodology for modeling oscillatory data that show regime shifts in periodicity, amplitude and phase. In contrast to previous work our approach does not require prespecifying the number of regimes or the order of the model within a regime. We assume that, conditional on the position and number of change-points, the time series can be approximated by a piecewise changing sinusoidal regression model. The timing and number of changes are unknown, along with the number and values of relevant periodicities in each segment. We develop a reversible jump MCMC technique that jointly explores the parameter space of the change-points and sub-models for all segments.

The article is organized as follows. Sections 2 and 3 present the model, the prior specifications, and the general structure of our Bayesian approach. Sections 4 and 5 provide a detailed explanation of our sampling scheme and simulation studies to demonstrate the performance of our approach. In Section 6, we illustrate the use of our methodology in two data-rich scenarios related to sleep, circadian rhythm, and e-Health research, namely the identification of the spectral properties of experimental breathing traces arising in sleep apnea research, and the analysis of human temperature data measured over several days by a wearable sensor. We conclude and discuss our current findings in Section 7.

2 The Model

Consider a time series realization $y_1,\dots ,y_n$ whose periodic behavior may change at k unknown time-points ${\mathit{s}}_{ (k)}=(s_1,\dots ,s_k)\prime$ where k is also unknown. Assume that in each sub-interval ${\mathrm{I}}_j=[s_{j-1},s_j)$ there are m_j relevant frequencies ${\mathit{\omega }}_{ j}=({\omega }_{j, 1},\dots ,{\omega }_{j, m_j})\prime$, for $j=1,2,\dots ,k+1$. Setting $s_0=1$ and $s_{k+1}=n$, we can write the following sinusoidal model (Andrieu and Doucet 1999)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) where(2) $$\begin{array}{c}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )={\alpha }_j+{\mu }_j t\\ \hspace{1em}+\sum _{l=1}^{m_j}({\beta }_{j, l}^{ \left(1\right)}\hspace{0.17em}\text{cos}\hspace{0.17em}(2\pi {\omega }_{j, l} t)+{\beta }_{j, l}^{ \left(2\right)}\hspace{0.17em}\text{sin}\hspace{0.17em}(2\pi {\omega }_{j, l} t\left)\right),\end{array}$$(2)

${\mathit{\beta }}_{ j}=({\alpha }_j,{\mu }_{ j}, \mathit{\beta }{\prime }_{j, 1},\dots , \mathit{\beta }{\prime }_{j, m_j})\prime$, ${\mathit{\beta }}_{j, l}=( {\beta }_{j,l}^{ (1)}, {\beta }_{j,l}^{ (2)} )\prime$, $1_{[·]}$ denotes the indicator function, and μ_j and α_j may, if needed, account for a linear trend within each segment. For simplicity we assume independent zero-mean Gaussian errors with regime-specific variances(3) $${\epsilon }_t\sim \mathcal{N} (0,{\sigma }_j^2),\hspace{1em}\text{for} t\in {\mathrm{I}}_j \text{and} j=1,\dots ,k+1,$$(3) noting that the methodology can in principle be extended to the non-Gaussian case.

The dimension of the model is given by the number of change-points k and the number of frequency components in each regime denoted by ${\mathit{m}}_{ (k)}=(m_1,\dots ,m_{k+1})\prime$. Furthermore, let ${\mathit{\beta }}_{ (k)}=( \mathit{\beta }{\prime }_{ 1},\dots ,\mathit{\beta }{\prime }_{ k+1})\prime ,{\mathit{\omega }}_{ (k)}=(\mathit{\omega }{\prime }_{ 1},\dots ,\mathit{\omega }{\prime }_{ k+1})\prime$, ${\mathit{\sigma }}_{ (k)}^{ 2}=({\sigma }_{ 1}^2,\dots ,{\sigma }_{ k+1}^{ 2})\prime$, and ${\mathit{\theta }}_{ (k)}=( \mathit{\beta }{\prime }_{ (k)}, \mathit{\omega }{\prime }_{ (k)}, {\mathit{\sigma }}_{ (k)}^{ 2\prime } )\prime$. Using Equation (1)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) , the likelihood of $( k, {\mathit{m}}_{ (k)}, {\mathit{s}}_{ (k)},{\mathit{\theta }}_{ (k)} )$ given the data $\mathit{y}=( y_1,\dots ,y_n )\prime$ is(4) $$\begin{array}{c}\mathit{L}( k, {\mathit{m}}_{ \left(k\right)}, {\mathit{s}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)}, | \mathit{y} )=\prod _{j=1}^{k+1}\mathit{L}( m_j, {\mathit{\theta }}_{ j} | {\mathit{y}}_j ),\\ {\mathit{y}}_j=( y_t: t \in {\mathrm{I}}_j ),\end{array}$$(4) where(5) $$\begin{array}{c}\mathit{L}( m_j, {\mathit{\theta }}_{ j} | {\mathit{y}}_j )={( 2\pi {\sigma }_j^2 )}^{ -n_j/2}\\ \hspace{1em}\times \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }_j^2} \sum _{ t \in {\mathrm{I}}_j}^{}{\left\{ y_t-{\mathit{x}}_t \left({\mathit{\omega }}_j\right)\prime {\mathit{\beta }}_{ j} \right\}}^{ 2}\right],\end{array}$$(5)

${\mathit{\theta }}_{ j}=( \mathit{\beta }{\prime }_j, \mathit{\omega }{\prime }_j, {\sigma }_{ j}^{ 2\prime } )\prime$ is the vector of parameters, n_j the number of observations of the jth segment, and the vector of basis functions ${\mathit{x}}_t ( {\mathit{\omega }}_j )$ is defined as$$\begin{array}{c}{\mathit{x}}_t \left({\mathit{\omega }}_j\right)=(1, t, \hspace{0.17em}\text{cos}\hspace{0.17em}(2\pi {\omega }_{j,1}t), \hspace{0.17em}\text{sin}\hspace{0.17em}(2\pi {\omega }_{j,1}t),\dots ,\\ \hspace{1em}\hspace{0.17em}\text{cos}\hspace{0.17em}\left(2\pi {\omega }_{j,m_j}t\right), \hspace{0.17em}\text{sin}\hspace{0.17em}\left(2\pi {\omega }_{j,m_j}t\right))\prime .\end{array}$$

3 Bayesian Inference

Given some prefixed maximal numbers of change-points, $k_{\text{max}}$, and frequencies per regime, $m_{\text{max}}$, inference is achieved by assuming that the true model is unknown but comes from a finite class of models where each model ${\mathcal{M}}_k$, with k change-points, is parameterized by the vector$$( {\mathit{m}}_{ \left(k\right)}, {\mathit{s}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)} )\in {\Pi }_k,\hspace{1em}{\Pi }_k\in \Pi .$$

Let ${\mathit{S}}_k=\{ {\mathit{s}}_{ (k)}\in {[1, n]}^{ k} : 1_{[1<s_1<\cdots <s_k<n]} \}$ and ${\Omega }_{m_j}={(0,0.5)}^{ m_j}$ denote, respectively, the sample space for the locations of change-points and the frequencies of the jth segment. The overall parameter space can be written as a finite union of subspaces$$\begin{array}{c}\Pi =\underset{k=0}{\overset{k_{\text{max}}}{\cup }}\{ k \}\times {\Pi }_k,\hspace{1em}\text{and}\\ {\Pi }_k={\mathit{S}}_k \times \prod _{j=1}^{k+1}\left\{m_j\right\} \times \underset{m_j=1}{\overset{m_{\text{max}}}{\cup }}\{ R^{2m_j+2}\times {\Omega }_{m_j}\times R^{+}\}.\end{array}$$

Bayesian inference on $k, {\mathit{m}}_{ (k)},{\mathit{s}}_{ (k)}$ and ${\mathit{\theta }}_{ (k)}$ is based on the following factorization of the joint posterior distribution$$\begin{array}{c}\pi ( k, {\mathit{m}}_{ \left(k\right)}, {\mathit{s}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)} | \mathit{y} )=\pi (k | \mathit{y}) \pi ({\mathit{m}}_{ \left(k\right)} | k, \mathit{y})\\ \hspace{1em}\pi ({\mathit{s}}_{ \left(k\right)} | {\mathit{m}}_{ \left(k\right)}, k, \mathit{y}) \pi ({\mathit{\theta }}_{ \left(k\right)} | {\mathit{s}}_{ \left(k\right)}, k, {\mathit{m}}_{ \left(k\right)}, \mathit{y}),\end{array}$$where we use $\pi ( · )$ as generic notation for probability density or mass function, whichever is appropriate. Sampling from it poses a multiple model selection problem, namely of the number of change-points and number of frequencies in each regime, which can be addressed by constructing a reversible-jump MCMC algorithm Green (1995). The algorithm in its basic structure iterates between the following two moves:

1. Segment model move: Given a partition of the data at k locations ${\mathit{s}}_{ (k)}$, inference on the parameters ${\mathit{m}}_{ (k)}$ and ${\mathit{\theta }}_{ (k)}$ is based on the conditional posterior$$\pi ({\mathit{m}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)} | k, {\mathit{s}}_{\left(k\right)}, \mathit{y})=\prod _{j=1}^{k+1}\pi (m_j, {\mathit{\theta }}_{ j} | k, {\mathit{s}}_{ \left(k\right)}, {\mathit{y}}_j).$$ A reversible-jump MCMC algorithm is performed in parallel on each of the k + 1 segments, where at each iteration the number of sinusoids m_j, the linear coefficients ${\mathit{\beta }}_{ j}$, the frequencies ${\mathit{\omega }}_{ j}$, and the residual variances ${\sigma }_j^2$ are sampled independently in each segment, for $j=1,\dots ,k+1$. Notice that at this stage the algorithm will explore subspaces of variable dimensionality regarding the number of frequencies per segment, while the change-point model remains fixed.

2. Change-point model move: This step performs a reversible-jump MCMC algorithm for change-point model search where the number k and locations of change-points ${\mathit{s}}_{ (k)}$ are sampled, along with the linear coefficients, number of frequencies and their values as well as the residual variances for any segments affected by the move.

Our prior specifications assume independent Poisson distributions for the number of break-points k and frequencies in each segment m_j, conditioned on $k\le k_{\text{max}}$ and $1\le m_j\le m_{\text{max}}$, respectively. Given k, a prior distribution for the positions of the change-points ${\mathit{s}}_{ (k)}$ can be chosen as in Green (1995)(6) $$\begin{array}{c}\pi ({\mathit{s}}_{ \left(k\right)} | k)=\frac{(2k+1)!}{{(n-1)}^{2k+1}}\prod _{j=0}^k(s_{j+1}-s_j) 1_{[s_0<s_1<\cdots <s_k<n]},\\ \hspace{1em}{s}_0=1,s_{k+1}=n.\end{array}$$(6)

Conditional on k and ${\mathit{m}}_{ (k)}$, we choose a uniform prior for the frequencies ${\omega }_{j, l}\sim \text{Uniform}(0,0.5),\hspace{1em}l=1,\dots ,m_j, \text{and} j=1,\dots ,k+1.$ Analogous to a Bayesian regression (Bishop 2006), a zero-mean isotropic Gaussian prior is assumed for the coefficients of the jth segment, ${\mathit{\beta }}_{ j}\sim {\mathit{N}}_{2m_j+2}( \mathit{0}, {\sigma }_{\beta }^2 \mathit{I} ),\hspace{1em}j=1,\dots ,k+1,$ where ${\sigma }_{\beta }^2$ is a prespecified large value, and the prior on the residual variance ${\sigma }_j^2$ of the jth partition is $\text{Inverse}-\text{Gamma} (\frac{{\nu }_0}{2},\frac{{\gamma }_0}{2})$, where η₀ and ν₀ are fixed at small values.

4 Sampling Scheme for Nonstationary Periodic Processes

Here we provide the sampling scheme associated with the nonstationary periodic processes that we wish to model. An outline of the overall procedure is as follows. Start with an initial configuration of number of change-points k, along with their locations ${\mathit{s}}_{ (k)}$; this yields a partition of the data $\mathit{y}=({\mathit{y}}_1,\dots ,{\mathit{y}}_{k+1})$. Initialize the number of frequencies in each regime ${\mathit{m}}_{ (k)}$ and their values ${\mathit{\omega }}_{ (k)}$, along with the coefficients ${\mathit{\beta }}_{ (k)}$ and residual variances ${\mathit{\sigma }}_{ (k)}^2.$ At each iteration of the algorithm a segment model and a change-point model move are estimated. A random choice with probabilities (7) based on the current number of parameters will determine whether to attempt a birth, death or a within-model move. In particular, let z denote the current number of parameters, that is, change-points k in the change-point model or frequencies m_j in the jth segment model; then, the dimension may increase by one (birth step) with probability b_z, decrease by one (death step) with probability d_z or remain unchanged (within step) with probability ${\mu }_z=1-b_z-d_z$, where(7) $$b_z=c \text{min}\{1,\frac{\pi (z+1)}{\pi \left(z\right)}\},\hspace{1em}{d}_{z+1}=c \text{min}\{1,\frac{\pi \left(z\right)}{\pi (z+1)}\},$$(7) for some constant $c\in [0,\frac{1}{2}]$, and $\pi (z)$ is the prior probability of the model including z. Reversibility of the Markov chain is guaranteed for move types that involve a change in dimensionality as $b_z \pi (z)=d_{z+1} \pi (z+1).$ Here we chose c = 0.4 but other values are legitimate as long as c is not larger than 0.5, to assure that the sum of the probabilities does not exceed 1 for some values of z. Naturally, $b_{k=k_{\text{max}}}=b_{m=m_{\text{max}}}=0$ and $d_{k=0}=d_{m=1}=0$. The pseudocode of the overall algorithm that describes an iteration of the sampler is given in Algorithm 1. We next describe in more detail the specific procedures needed to update the moves.

Algorithm 1:

1. For each segment $j=1,\dots ,k+1$, perform a segment model move (Section 4.1)

Draw $U\sim \text{Uniform} (0,1)$

if $U\le b_{m_j} \to$ birth-step

else if $b_{m_j}\le U\le d_{m_j} \to$ death-step

else $\to$ within-step

2. Perform a change-point model move (Section 4.2):

Draw $U\sim \text{Uniform} (0,1)$

if $U\le b_k \to$ birth-step

else if $b_k\le U\le d_k \to$ death-step

else $\to$ within-step

4.1 Updating a Segment Model

Given the number of change-points k and their locations ${\mathit{s}}_{ (k)}$, a segment model move is performed independently and in parallel on each of the k + 1 partitions. Hence, throughout this subsection the subscript relating to the jth segment may be dropped and a segment of interest is denoted by $\mathit{y}=(y_a,\dots ,y_b)\prime$, which contains n observations. Assume that the current number of frequencies is set at m; then, an independent random choice is made between attempting a birth, death or within-model step, with probabilities given in (7). An outline of these moves is as follows (further details are provided in the Appendix).

4.1.1 Within-Model Move

Conditioned on the number of frequencies m, we sample the vector of frequencies $\mathit{\omega }$ following Andrieu and Doucet (1999), that is, by sampling the frequencies one-at-time using a mixture of M-H steps, with target distribution(8) $$\begin{array}{c}\pi (\mathit{\omega } | \mathit{\beta }, {\sigma }^2, m, \mathit{y})\\ \hspace{1em}\propto \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }^2}\sum _{t = a}^b\{y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2\right]1_{[ \mathit{\omega } \in {\Omega }_m] }.\end{array}$$(8)

In particular, the proposal distribution is a combination of a normal random walk centered around the current frequency and a sample from values of the discrete Fourier transform of y. The corresponding vector of linear parameters $\mathit{\beta }$ is then updated in a M-H step, from the target posterior conditional distribution(9) $$\begin{array}{c}\pi ( \mathit{\beta } | \mathit{\omega }, {\sigma }^2, m, \mathit{y})\\ \hspace{1em}\propto \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }^2}\sum _{t = a}^b\{y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2-\frac{1}{2{\sigma }_{\beta }^2} \mathit{\beta }\prime \mathit{\beta }\right],\end{array}$$(9) where the proposed values are drawn from normal approximations to their posterior conditional distribution. Finally, the residual variance σ² is then updated in a Gibbs step from(10) $$\begin{array}{c}{\sigma }_{| \mathit{\omega }, \mathit{\beta }}^2\sim \text{Inverse}-\text{Gamma}\\ \hspace{1em}\left( \frac{n+{\nu }_0}{2}, \frac{{\gamma }_0+\sum _{t = a}^b\{ y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2}{2}\right).\end{array}$$(10)

4.1.2 Between-Model Moves

For this type of move, the number of frequencies is either proposed to increase by one (birth) or decrease by one (death). If a birth move is proposed, we have that $m^{ p}=m^{ c}+1$, where current and proposed values are denoted by the superscripts c and p, respectively. The proposed vector of frequencies is constructed by proposing an additional frequency to include in the current vector. Conditional on the frequencies, the corresponding vector of linear coefficients and the residual variance are sampled as in the within-model move. If a death move is proposed, we have that $m^p=m^c-1$. Hence, one of the current frequencies is randomly chosen to be removed. The proposed corresponding vector of linear coefficients is drawn, along with the residual variance. For both moves, the updates are jointly accepted or rejected in a M-H step.

4.2 Updating the Change-Point Model

This part of the algorithm identifies the number and locations of change-points. Suppose the number of change-points is currently set to some value k, then according to the probabilities given in (7) a random decision is made between adding, removing, or moving a change-point. The rules for updating these types of moves are described below and more details are given in the Appendix.

4.2.1 Within-Model Move

An existing change-point is proposed to be relocated with probability $\frac{1}{k}$, obtaining say $s_j^{ c}$. The update for the selected change-point is proposed from a mixture of a normal random walk centered on the current change-point $s_j^{ c}$ and a sample from a uniform distribution on the interval $[s_{j-1}^{ c}+{\psi }_s,s_{j+1}^{ c}-{\psi }_s]$. Here, we introduced ψ_s as a fixed minimum time between change-points avoiding change-points being too close to each other. Rosen, Wood, and Stoffer (2012) used a similar scheme, but on a discrete-scale. The number of frequencies and their values are kept fixed, and, conditional on the relocation, the linear coefficients for the segments affected by the relocation are sampled. These updates are jointly accepted or rejected in a M-H step and the residual variances are updated in a Gibbs step.

4.2.2 Between-Model Moves

For this type of move, the number of change-points may either increase (birth) or decrease (death) by one. If a birth move is proposed, we have that $k^{ p}=k^{ c}+1$. The new proposed change-point is drawn uniformly on $f ({\mathit{s}}_{ (k^{ c})}^{ c}, {\psi }_s)$, the support of ${\mathit{s}}_{ (k^{ c})}^{ c}$ given the constraints imposed by ψ_s, that is, $f ({\mathit{s}}_{ (k^{ c})}^{ c}, {\psi }_s)=[1+{\psi }_s, s_1^{ c}-{\psi }_s] \cup [s_1^{ c}+{\psi }_s, s_2^{ c}-{\psi }_s]\cup \cdots \cup [s_{k^{ c}}^{ c}+{\psi }_s, n-{\psi }_s]$. The latter involves splitting an existing segment. The number of frequencies and their values in the proposed segments are selected from the current states. Two residual variances for the new proposed segments are then constructed from the current single residual variance. Finally, two new vectors of linear parameters are sampled. If a death move is proposed, we have that $k^{ p}=k^{ c}-1$. Hence, a candidate change-point to be removed is selected from the vector of existing change-points, with probability $\frac{1}{k^{ c}}$. The latter involves merging two existing partitions. The number of frequencies and their values in the proposed segments are selected from the current states. A single residual variance is constructed from the current variances relative to the segments affected by the relocation. Finally, a new vector of linear coefficient is drawn. For both type of moves, these updates are jointly accepted or rejected.

5 Simulation Studies

We carry out simulation studies to explore the performance of our method, which will be referred to as Automatic Nonstationary Oscillatory Modeling (AutoNOM). In Section 5.1, we illustrate the performance of our methodology when the simulated data are generated from the proposed model. Section 5.2 deals with scenarios when the model is misspecified relative to the generating process. Our results are compared with two state-of-the-art existing methods.

5.1 Illustrative Example

In this simulation example, we generate a time series consisting of n = 900 data points from model (1) with k = 2 change-points located at positions ${\mathit{s}}_{ (2)}=(300,650)$, and fixed number of frequencies per regime ${\mathit{m}}_{ (2)}=(3,1,2)$. (Further details of the parameterization are available in supplementary materials, Section 1.1.) Figure 1 (top panel) shows a realization from this model. The prior means ${\lambda }_{\omega }$ and λ_s, say, on the number of frequencies and change-points, respectively, were set to 2, to reflect a fair degree of prior information on their numbers. We discuss in Section 1.2 of the supplementary materials that AutoNOM was relatively insensitive to these prior specifications, for this example. The maximum number of change-points $k_{\text{max}}$ was set to 15, and the maximum number of frequencies per regime $m_{\text{max}}$ was set to 10. Furthermore, we fixed ${\psi }_s=20$ and ${\varphi }_{\omega }$ = 0.25 (Appendix A.1.2) for the uniform distribution for sampling the frequencies. The full estimation algorithm was ran for 20,000 updates, 5000 of which are discarded as burn-in period. The estimation took 390 sec with a (serial) program written in Julia 0.62 on a Intel® Core™ i7-4790S Processor 16 GB RAM. The results, summarized in Table 1 clearly show that a model with two change-points has the highest estimated posterior probability (left panel) and that AutoNOM correctly identifies the right number of significant frequencies in each regime (right panel).

Fig. 1 Illustrative example. (Top) Simulated time series. (Middle) Estimated posterior distribution for the location of the change-points, conditioned on k = 2. The dotted vertical lines represent true location of change-points. (Bottom) Estimated posterior distribution of the frequencies for each different segment, conditioned on k = 2, $m_1=3,m_2=1$, and $m_3=2$. The dotted vertical lines represent true values of the frequencies.

Table 1 Illustrative example.

k	$\widehat{\pi } ( k | \mathit{y})$	m	$\widehat{\pi } ( m_1 | k=2, \mathit{y})$	$\widehat{\pi } ( m_2 | k=2, \mathit{y})$	$\widehat{\pi } ( m_3 | k=2, \mathit{y})$
0	0.00	1	0.00	0.99	0.00
1	0.02	2	0.00	0.01	0.98
2	0.97	3	0.98	0.00	0.02
3	0.01	4	0.02	0.00	0.00
4	0.00	5	0.00	0.00	0.00

NOTES: (Left panel) posterior probabilities for number of change-points; (right panel) posterior probabilities for number of frequencies in each regime, conditioned on k = 2.

Figure 1 (middle panel) shows the estimated posterior distribution for the location of the change-points, conditioned on three segments. The posterior means of the change-point locations are $\widehat{\mathbb{E}} (s_1 | k=2, \mathit{y})=298.7$ and $\widehat{\mathbb{E}} (s_2 | k=2, \mathit{y})=650.1$. Figure 1 (bottom panel) shows that the estimated posterior distributions are an excellent match to the true frequencies. In addition, we provide details about acceptance rates in supplementary materials, Section 2.

5.1.1 Detecting Spectral Peaks

We simulate a time series from the same simulation model as above with the only difference that the residual variances were set equal to one for all segments and thus are smaller than above. The performance of AutoNOM is compared with two existing methods, namely the Bayesian adaptive spectral estimation for nonstationary time series proposed by Rosen, Wood, and Stoffer (2012), referred to as AdaptSPEC, and the frequentist piecewise vector autoregressive method of Davis, Lee, and Rodriguez-Yam (2006), referred to as AutoPARM. Specifically, we explore the performances of these methodologies in identifying the number and location of change-points, and the number and location of frequency peaks in each estimated segment. AdaptSPEC requires the user to specify in advance the number of basis function J used for smoothing the periodogram in the segments. We run AdaptSPEC for two different specifications, namely J = 7 and J = 15 basis functions. The model is fitted with a total of 15,000 iterations, 5,000 of which are discarded as burn-in, by using the R package provided by the authors. Posterior samples of peak frequencies are obtained by considering the modes of the spectrum per MCMC iteration. AutoPARM is performed with default tuning parameters. We note that Davis, Lee, and Rodriguez-Yam (2006) do not discuss computation of confidence intervals for frequencies.

The modal number of change-points for AdaptSPEC is 2 for both J = 7 and J = 15, with posterior probability $\widehat{\pi } (k=2 | \mathit{y})$ of 76% and 88%, respectively; the modal number of change-points for AutoNOM is 2 and AutoPARM identifies 2 change-points as well. Conditioned on the modal number of change-points, Table 2 displays the estimated location of changes (left panel) and frequency peaks (right panel) for the different compared methods, where we report the standard deviation for the estimate obtained from the empirical distribution of the posterior samples. Similarly, we show in Figure 2 the estimated location of the frequency peaks and their 95% credible intervals, for each of the three identified segments; dotted vertical lines represents the true location of the frequency peaks. Results for AutoNOM are conditioned on the modal number of frequencies per regime.

Fig. 2 Illustrative example with unitary residual variances. Estimated frequency peaks for AutoNOM (AN), AdaptSPEC (AS, $J=7,15$), and AutoPARM (AP); 95% credible intervals (horizontal lines) are also reported for Bayesian methods. Dotted vertical lines are true locations of the frequency peaks.

Table 2 Illustrative example with unitary residual variances.

	s1	s2	${\omega }_{1,1}$	${\omega }_{1,2}$	${\omega }_{1,3}$	${\omega }_{2,1}$	${\omega }_{3,1}$	${\omega }_{3,2}$
True	300	650	0.042	0.067	0.143	0.083	0.046	0.067
AN	300.87	650.46	0.042	0.067	0.142	0.083	0.045	0.067
AN	(0.97)	(0.31)	(0.0004)	(0.0002)	(0.0005)	(0.0002)	(0.0003)	(0.0003)
AS J7	319.65	628.51	0.083	–	–	0.082	0.057	–
AS J7	(15.36)	(37.35)	(0.0030)	–	–	(0.0070)	(0.0010)	–
AS J15	298.5	647.01	0.057	0.141	–	0.088	0.056	–
AS J15	(1.21)	(0.14)	(0.0030)	(0.0020)	–	(0.0050)	(0.0030)	–
AP	299	648	0.060	0.140	–	0.080	0.055	–

NOTES: Estimated change-points locations (left panel) and frequency peaks (right panel) for AutoNOM (AN), AdaptSPEC (AS J = 7, 15), and AutoPARM (AP); posterior standard deviations are also reported for Bayesian methods.

It becomes clear that the detection of periodicities by AdaptSPEC is affected by the specification of the number of spline basis functions used for the smoothing, where increasing the number of basis function yields a better performance for AdaptSPEC. The example also shows that smoothing by splines may lead to peaks in the periodogram to be over-smoothed and neighboring close peaks to be merged. AutoPARM seems to also suffer from the latter problem.

When we increased the residual variance to the high levels set originally, AdaptSPEC failed to detect any change-points for both J = 7 and J = 15, with posterior probability $\widehat{\pi } (k=0 | \mathit{y})$ of 69% and 93%, respectively, while AutoPARM found 7 change-points and thus severely overestimates their number. Our conclusion from this comparison is that although AdaptSPEC and AutoPARM may be well suited for time series processes with smooth time-varying spectra with few or no peaks, both methods are severely challenged in detecting changes in spectra that exhibit pronounced peakedness, possibly at nearby frequencies, as can be expected to occur in reality for the type of time series that we wish to analyze.

5.2 Misspecified Model

We investigate the performance of our proposed method for identifying spectral peaks when the model is misspecified relative to the generating process. In particular, we explored simulation studies under three different settings. In the first two scenarios we generated data from two types of autoregressive (AR) processes, namely a piecewise AR process and a slowly varying AR process. We compare the performance of our procedure with AutoPARM and AdaptSPEC. In the third setting, we assumed that the innovations are t-distributed, and therefore violate the Gaussianity assumption of ${\epsilon }_t$ in Equation (3)(3) $${\epsilon }_t\sim \mathcal{N} (0,{\sigma }_j^2),\hspace{1em}\text{for} t\in {\mathrm{I}}_j \text{and} j=1,\dots ,k+1,$$(3) . For all models, our estimation algorithm was run for 20,000 iterations, 5000 of which were used as burn in, and the hyperparameters were chosen as ${\varphi }_s=40,{\lambda }_{\omega }=0.05$ and ${\lambda }_s=0.01$.

5.2.1 Piecewise Autoregressive Process

As pointed out by a referee, although modeling a time series as a linear combination of a finite number of sinusoids plus noise is common in the signal processing literature, such line-spectrum based models are rare in the statistics literature. In fact, it is commonly assumed that the power spectrum is continuous across frequencies. We investigate the performance of the proposed procedure when analyzing data generated from a piecewise AR process whose local spectral density functions show sharp peaks. Specifically, a realization is simulated from(11) $$y_t=\{\begin{array}{cc}1.9 y_{t-1}–0.975 y_{t-2}+{\epsilon }_t^{\left(1\right)}\hfill & \text{for} 1\le t\le 250\hfill \\ 1.9 y_{t-1}–0.991 y_{t-2}+{\epsilon }_t^{\left(2\right)}\hfill & \text{for} 251\le t\le 400\hfill \\ -1.35 y_{t-1}–0.37 y_{t-2}+0.36 y_{t-3}+{\epsilon }_t^{\left(3\right)}\hfill & \text{for} 401\le t\le 550,\hfill \end{array}$$(11) where ${\epsilon }_t^{(1)}\stackrel{\text{iid}}{\sim }\mathcal{N}(0,0.25)$ and ${\epsilon }_t^{(i)}\stackrel{\text{iid}}{\sim }\mathcal{N}(0,1)$ for i = 2, 3. Figure 3 (top panel) shows a realization from model (11). After applying our methodology AutoNOM, the posterior probability of two change-points is 97.93% and the posterior means of the change-point locations are $\widehat{\mathbb{E}} (s_1 | k=2, \mathit{y})=251.19$ and $\widehat{\mathbb{E}} (s_2 | k=2, \mathit{y})=401.56$. The estimated location of the frequency peaks for our proposed procedure in comparison to AdaptSPEC and AutoPARM and the true values are shown in Figure 3 (bottom panels). It is evident that the proposed and existing methodologies successfully identify the true location of the frequency peaks in each segment, with AdaptSPEC showing less precision.

Fig. 3 Piecewise AR process. (Top) A realization from model (11). Vertical dotted lines are the estimated locations of the change-points. (Bottom) Estimated frequency peaks for AutoNOM (AN), AdaptSPEC (AS J = 10), and AutoPARM (AP); 95% credible intervals (horizontal lines) are also reported for Bayesian methods. Dotted vertical lines are true locations of the frequency peaks.

5.2.2 Slowly Varying Autoregressive Process

In this section, we analyze an AR process whose continuous spectral density is changing slowly over time. We note though that this scenario is a large departure from the assumptions of our model. In particular, we consider the same slowly varying AR(2) process investigated by Ombao et al. (2001) and Davis, Lee, and Rodriguez-Yam (2006), namely(12) $$y_t=a_t y_{t-1}-0.81 y_{t-2}+{\epsilon }_t,\hspace{1em}t=1,\dots ,1031,$$(12) where $a_t=0.8 [1-0.5\hspace{0.17em}\text{cos}\hspace{0.17em} (\pi t/1031)]$ and ${\epsilon }_t\stackrel{\text{iid}}{\sim }\mathcal{N}(0,1)$. Notice that the parameter a_t is changing gradually over time whereas the coefficient associated with the second lag remains constant. A realization from model (12) is shown in Figure 4(a) and the corresponding time varying frequency peak is displayed in Figure 4(b) as a solid line. Figure 4(b) also shows the estimated time varying frequency peak for AutoNOM, AdaptSPEC, and AutoPARM. For AutoNOM and AdaptSpec, the time changing frequency peak has been averaged across the MCMC samples, giving a smoother estimate (especially for AdaptSPEC) than the one obtained by AutoPARM. For each method, we compute the residual sum of squares $\text{RSS}={\sum }_{t=1}^{1031}{({\omega }_t-{\widehat{\omega }}_t)}^2$ between the true time changing frequency peak ω_t and its estimate ${\widehat{\omega }}_t$. The RSS in this example was 0.111, 0.174, and 0.085 for AutoNOM, AdaptSPEC, and AutoPARM, respectively. It is clear that even in this scenario where the data generating model was very different from the underlying assumptions of our model, our approach seems to outperform AdaptSPEC and remains competitive with AutoPARM in estimating the time varying frequency peak.

Fig. 4 Slowly varying AR(2) process. (a) A realization from model (12). (b) True time varying frequency peak (solid line) and estimated time varying frequency peak for AutoNOM (AN), AdaptSPEC (AS J = 10), and AutoPARM (AP).

5.2.3 Non-Gaussian Time Series

We investigate the performance of our approach in the scenario when the innovations are t-distributed. We simulate a time series from the same simulation model presented in Section 5.1, where errors were generated from a t-distribution with 2, 3, and 2 degrees of freedom for the sequence of three segments, respectively. The degrees of freedom were chosen low such that the corresponding distributions show heavy tails. A realization of this time series is shown in Figure 5. Our proposed methodology correctly identifies the 2 change-points, as the estimated posterior probability $\widehat{\pi } (k=2 | \mathit{y})$ is 0.99. The posterior means of the change-point locations are $\widehat{\mathbb{E}} (s_1 | k=2, \mathit{y})=303.6$ and $\widehat{\mathbb{E}} (s_2 | k=2, \mathit{y})=650.5$, showing an excellent match to the true values ${\mathit{s}}_{ (2)}=(300,650)$. Furthermore, the posterior mode of the number of frequencies in each segment is ${\widehat{\mathit{m}}}_{ (2)}=(3,1,2)$, which is a correct estimate of ${\mathit{m}}_{ (2)}=(3,1,2)$. We also display in Figure 5 the estimated signal (using Equation (1)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) , supplementary materials) as a dotted line. We can conclude that, although our model assumes Gaussianity, AutoNOM seems to perform well even in the case where the oscillatory underlying process is t-distributed with heavy tails.

Fig. 5 Illustrative example with t-distributed residual variances. Simulated time series (solid line) and estimated signal (dotted line). The dotted vertical lines represent the estimated location of the change-points.

6 Case Studies

The development of our methodology was motivated by the following two case studies where dense physiological signals were observed which exhibit unknown periodicities whose role may change over time in a more or less abrupt manner and where their detection is of relevance to the health and well-being of the subject.

6.1 Analysis of Human Skin Temperature

The development of information and communication technologies, in particular widespread internet access and availability of mobile phones and tablets, allows considering new developments in the health care system. To address the issue of personalized medical treatment according to the circadian timing system of the patient, referred to as chronotherapy (Levi and Schibler 2007), a novel and validated noninvasive mobile e-Health platform pioneered by the French project PiCADo (Komarzynski et al. 2018) is used to record and teletransmit skin surface temperature as well as physical activity data (Huang et al. 2018) from an upper chest e-sensor. Figure 6(a) shows an example of 4 days of 5-min temperature recording for a healthy individual. The circadian rhythms in core and skin surface temperature are usually 8–12 hr out of phase, with respective maxima occurring near 16:00 at day time, and near 2:00 at night (Krauchi and Wirz-Justice 1994). The early night drop in core body temperature, which is critical for triggering the onset of sleep (Van Someren 2006), results from the vasodilatation of the skin vessels and associated rise in skin surface temperature (Kräuchi 2002). Under the assumption of stationarity, Komarzynski et al. (2018) analyzed the skin temperature time series identifying both strong 12 hr (circahemidian) and 24 hr (circadian) rhythms.

Fig. 6 Analysis of skin temperature of a healthy subject. Panel (a) are the time series of skin temperature and corresponding physical activity. Panel (b) is the estimated signal (solid line) along with its 95% credible interval; vertical lines are the estimated locations of the change-points. Panel (c) is the estimated posterior density histogram of the locations of the changes, conditioned on $k=7$ change-points. Rectangles on the time axis of each plot correspond to periods from 20.00 to 8.00. The variation in skin temperature finds analogies with the rest-activity pattern that alternates between day activity and night rest.

Here, we applied our methodology to the skin-temperature time series shown in Figure 6(a) for 300,000 iterations, discarding the first 100,000 updates as burn-in. The maximum number of change-points $k_{\text{max}}$ was set to 10, whereas the maximum number of frequencies per regime $m_{\text{max}}$ was set to 5. The estimated number of change-points had a mode at 7, with $\widehat{\pi } ( k=7 | \mathit{y})=0.97$ and their estimated posterior distributions are shown in Figure 6(c). Inspecting them alongside the physical activity data we can see that the change points mainly correspond to the start and endpoints of the prolonged rest periods at nights showing that skin temperature alternates between day activity and night rest including sleep. Figure 7 shows the estimated posterior distribution of the frequencies for the sleep segments (2, 4, 6, 8) along with the square root of the estimated power of the corresponding frequencies, where the power of each is frequency ${\omega }_{j,l}$ is summarized by the sum of squares of the corresponding linear coefficients, that is, $I ({\omega }_{j,l}) ={\beta }_{j,l}^{ {(1)}^{ 2}}+{\beta }_{j,l}^{ {(2)}^{ 2}}$ (Shumway and Stoffer 2005). Figure 6(b) shows the piecewise fitted signal, along with a 95% credible interval obtained from the 2.5 and 97.5 empirical percentiles of the posterior sample using Equation (1)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) , supplementary materials. Cycles of approximately 3 hr appear in segments 2, 4, and 6; cycles that range approximately 1–1.5 hr appear in segments 2, 4, 8 and cycles of around 2 hr appear in segments 4 and 6 while some longer periods identified in segments 4, 6, 8 indicate the presence of a trend.

Fig. 7 Spectral properties for segments corresponding to night rest. Estimated posterior distribution of the frequencies along with square root of the estimated power of the corresponding frequencies. The results are conditional on the modal number of frequencies per segment.

Stages of sleep are characterized by ultradian oscillations between rapid eye movements (REM) and non-REM. The biological functionality that regulates the alternations between these two types of sleep is not yet much understood (Altevogt and Colten 2006). However, several physiological changes that occur over night differ between REM and non-REM phases, such as heart rate, brain activity, muscle tone and body temperature (Berlad et al. 1993; Pace-Schott and Hobson 2002). The body cycles between REM and non-REM sleep stages with an average length that ranges approximately between 70 and 120 min, and there are usually four to six of these sleep cycles each night (Carskadon and Dement 2005; Shneerson 2009). Our analysis was able to use skin temperature data alone to detect periods of sleep throughout the day and identify oscillatory behavior during the night, whose frequencies are compatible with ultradian oscillations between REM and different non-REM sleep stages.

A comparison with the current state-of-the-art methods, AutoPARM and AdaptSPEC, is provided in the supplementary materials, Section 4.1. Circadian and ultradian rhythmicity are expected because body temperature is known to be a circadian biomarker (Krauchi and Wirz-Justice 1994), but these existing methods fail. Furthermore, we notice that in the framework of analyzing circadian biomarker data, such as body temperature, a change in acrophase may be of interest to the clinician as this may be indicative of a disruption of the bodyclock. The methodology can indeed be used to investigate phase which can be computed from the sinusoidal function that characterizes the jth segment (see supplementary materials, Section 3).

6.2 Characterizing Instances of Sleep Apnea in Rodents

Sleep apnea is the temporary ($\ge$ 2 breaths) interruption of breathing during sleep. Moderate or severe ($\ge$ 15 events per hour) sleep apnea, occurs in about 50% of men and 25% of women over the age of 40 (Heinzer et al. 2015), with 91% of people with sleep apnea being undiagnosed (Tan et al. 2016). Sleep apnea is linked to many diseases. Patients with sleep apnea are at increased risk of: cardiovascular events (Lanfranchi et al. 1999), cancer (Nieto et al. 2012), liver disease (Sundaram et al. 2016), diabetes (Harsch et al. 2004), metabolic syndrome (Parish, Adam, and Facchiano 2007), cognitive decline (Osorio et al. 2016), and increased risk of dementia in the elderly (Lal, Strange, and Bachman 2012). The motivation of this research is to provide a statistical methodology that can be applied to analyze large breathing datasets resulting from in vivo plethysmograph studies in rats to characterize the occurrence of sleep apnea under different experimental conditions. If this could be attained, a concrete aid to the understanding of the pathological implications of this status could be provided to clinicians and experimental biologists.

An unrestrained whole-body plethysmograph is used to produce a breathing trace from freely behaving rats for periods of up to 3 hr. Plethysmographs were made using an 2 L air-tight box connected to a pressure transducer, with an air pump and outlet valve producing a flow rate of 2 L/min. Airflow pressure signals were amplified using Neurolog system (Digitimer) connected to a 1401 interface and acquired on a computer using Spike2 software (CED).

Apneas are subclassified as post-sigh apneas, if the preceding breath was at least 25% above the average amplitude of prior breaths, or spontaneous apneas, if there was no manifestation of a previous sigh (Davis and O’Donnell 2013). Airflow traces from the plethysmograph are shown in Figure 8 (left panels) and consist of three time series, which will be referred to as (a), (b), and (c). They correspond to different actions for this rat: (a) an alternation of sniffing and normal breathing; (b) spontaneous apnea followed by normal breathing; (c) normal breathing followed by a sigh, and a post-sigh apnea. We note that these actions were classified by eye by an experienced experimental researcher. Each time series contains 20,000 observations where the signal was sampled at 2000 Hz so that we have 2000 observations per second.

Fig. 8 Plots of the respiratory traces of a rat (left panels) and corresponding estimated posterior power (right panels). Panel (a) is characterized by an alternation of sniffing and normal breathing. Panel (b) is a plot of the trace of a spontaneous apnea, followed by normal breathing. Panel (c) shows normal breathing followed by a sigh, and a post-sigh apnea. Dotted vertical lines correspond to the estimated locations of the change-points.

Our procedure allows us to set an upper bound, ${\varphi }_{\omega }$, (Appendix A.1.2) for the uniform interval where the new frequencies are sampled. As the periodogram ordinates for these data were approximately zero for all frequencies larger than 0.01, we decided accordingly to set ${\varphi }_{\omega }=0.01$. The locations of the changes (vertical lines) are displayed in Figure 8 (left panels). The posterior power of the frequencies, for each time series, is shown in Figure 8 (right panels). These results are conditional on the modal number of change-points and the modal number of frequencies per segment. For each dataset, we summarize in Table 3 the spectral properties of each partition by displaying the periodicities corresponding to the first two largest values of the estimated power. When the rat is sniffing, (a), the air flow trace oscillates with a dominant period of approximately 0.2 sec, namely 5 cycles per second. Normal breathing, (a) and (b), is characterized by lower frequencies and lower magnitude than sniffing, by oscillating with a dominant period of around 0.5 sec, namely around 2 cycles every second. Apneas, (b) and (c), appear to be characterized by higher frequencies than normal breathing but with a lower power, with dominant periods of around 0.25 and 0.35 sec. Notice that in the first partition of (c), the highest value of the power corresponds to the frequency responsible for a sigh before apnea. Moreover, our methodology identifies different frequencies that explain the variation between the third and fourth partition of (c), leading to the hypothesis that there might be a time changing spectrum during the occurrence of an apnea instance. A comparison of our results with the results from AutoPARM and AdaptSPEC is provided in the supplementary materials, Section 4.2.

Table 3 Spectral properties of respiratory traces of a rat.

	(a)		(b)		(c)
	Period	Power	Period	Power	Period	Power
Segment 1	0.19	0.2272	0.20	0.0015	0.25	5.2540
	0.15	0.0883	0.25	0.0010	0.21	0.0915
Segment 2	0.20	0.2613	0.49	0.0117	0.27	0.0107
	0.28	0.0686	0.34	0.0099	0.42	0.0043
Segment 3	0.48	0.0692	0.48	0.0145	0.24	0.0365
	0.32	0.0491	0.28	0.0122	0.37	0.0095
Segment 4	0.58	0.0400	–	–	0.36	0.0253
	0.47	0.0251			0.24	0.0155
Segment 5	0.19	0.3231	–	–	–	–
	0.16	0.1044

NOTES: Periodicities (in seconds) corresponding to the first two largest values of the estimated power, for each time series (a), (b), and (c).

7 Summary and Discussion

We developed a novel Bayesian methodology for analyzing nonstationary time series that exhibit oscillatory behavior. Our approach is based on the assumption that, conditional on the position and number of change-points, the time series can be approximated by a piecewise changing sinusoidal regression model. The timing and number of changes are unknown, along with the number and values of relevant periodicities in each regime. Bayesian inference is performed via a reversible jump MCMC algorithm that can simultaneously estimate both the number and location of change-points, as well as the number, frequency and magnitude of sinusoids within each segment. Our methodology can be seen as a novel and relevant extension of the work in Andrieu and Doucet (1999) to the nonstationary setting.

We illustrated the utility of our methodology in two case studies. First, we analyzed human skin temperature time series data obtained from a wearable device, which exhibited unknown periodicities that changed over time in an abrupt manner. Our proposed methodology identified interesting oscillations whose frequencies are consistent with ultradian oscillations between REM and non-REM sleep stages. Second, we characterized the occurrence of sleep apnea in large breathing datasets resulting from in vivo plethysmograph studies on rodents. Our spectral investigation was able to distinguish very sharp peaks, corresponding to different nearby frequencies, that are responsible for the different actions of the rodent.

Although we have not discussed this in detail here, several diagnostics for monitoring convergence were carried out in both simulation and case studies. In particular, we verified that the target posterior distribution reached a stable regime by analyzing the trace plot of the log-likelihood across MCMC iterations (Marin and Robert 2007). We are aware that assessing convergence only based on this simple tool may sometimes be misleading since stable values of the log-likelihood could simply mean that the Markov chain is stuck in some local mode of the posterior distribution. Additionally, conditioned on the modal number of change-points and modal number of frequencies per regime, we have also monitored (within-model) convergence by analyzing the traces and running averages plots for all parameters across MCMC iterations, with satisfactory results. Comparable results were also obtained when running several chains starting at over-dispersed initial values. We notice that the diagnostic tool used by Bruce et al. (2018) and Li and Krafty (2019) to assess convergence for reversible jump MCMC samplers appears relevant. In the context of adaptive spectral analysis of nonstationary time series, they point out that although the number of partitions change across models, a power spectrum is defined at each time point. The power spectrum is modeled with a fixed number of splines, yielding a vector of summary measures of parameters that maintain the same interpretation across models in their samplers. However, our proposed sampler has a further layer of variable dimensionality, as not only the number and locations of the change-points may change from one iteration to the next, but also the number of frequencies in each segment are not fixed throughout the simulation.

We conclude this article by noticing that, although a Gaussian distribution is assumed, it is conceivable that our model can be extended to allow for other error distributions in Equation (1)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) . For example, a generalized linear model (McCullagh and Nelder 1989) may be used to model periodic count data by assuming that the observed data follows a Poisson distribution, that is, $y_t\sim \text{Poisson} ({\mu }_t)$. The logarithmic link function of the expected value μ_t of the response variable y_t may be expressed as $\text{log}\hspace{0.17em} ( {\mu }_t )={\sum }_{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}$, where the definitions of the variables are the same as for Equation (1)(1) $$y_t=\sum _{j=1}^{k+1}f (t, {\mathit{\beta }}_{ j}, {\mathit{\omega }}_{ j} )1_{[ t \in I_j ]}+{\epsilon }_t,$$(1) . Bayesian inference can, in principle, be achieved in a similar way as described in the article, namely by iterating between segment and change-point model moves, where the formulation of the acceptance probabilities and some proposal distributions need to be modified accordingly. We believe that such an extension would find use in several ranges of applications, for example, in studying population cycles in ecology and epidemiology, where the abundance of species are measured as count variables (White and Bennetts 1996; Bhaskaran et al. 2013; Bramness et al. 2015).

Supplementary Materials

Supplementary materials are available and include further details about simulation studies, acceptance rates, phase investigation, and performances of existing methods in the case studies. Code that implements the methodology and the data used in the case studies are available as online supplemental material and can be also found at https://github.com/Beniamino92/AutoNOM.

Acknowledgments

We wish to thank the referees, the associate editor, Dr. Paul Jenkins, Zeda Li, and Jack Jewson for their insightful and valuable comments. The work presented in this article was developed as part of the first author’s PhD thesis at the University of Warwick.

Additional information

Funding

B. Hadj-Amar was supported by the Oxford-Warwick Statistics Programme (OxWaSP) and the Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/L016710/1. B. Finkenstädt and F. Lévi were supported by the Medical Research Council (MRC), grant reference: MR/M013170/1. F. Lévi was partly supported by the Conseil Régional d’Île de France, the Conseil Régional de Champagne-Ardenne, Mairie de Paris and the Banque Publique d’Investissement (BPI France) through the Fonds Unique Interministériel 12 (PiCADo, contract 11017951), and the Institut de Recherches en Santé Publique from France (CLOCK-DOM1, grant 2014-BDCR-EC). R. Huckstepp was supported by the MRC, grant reference: MC/PC/15070.

Appendix A:

Details of the Sampling Scheme

A.1 Updating the Segment Model

A.1.1 Within-Model Move

Sampling $\mathit{\omega }$: To obtain samples from the conditional posterior distribution $\pi (\mathit{\omega } | \mathit{\beta }, {\sigma }^2, m, \mathit{y})$ (see Equation (8)(8) $$\begin{array}{c}\pi (\mathit{\omega } | \mathit{\beta }, {\sigma }^2, m, \mathit{y})\\ \hspace{1em}\propto \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }^2}\sum _{t = a}^b\{y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2\right]1_{[ \mathit{\omega } \in {\Omega }_m] }.\end{array}$$(8) ), we draw the frequencies one-at-time using a mixture of M-H steps. To explore the parameter space efficiently, we design a mixture distribution $q ( {\omega }_l^{ p} | {\omega }_l^{ c} )$, so that(A.1) $$\begin{array}{c}q ( {\omega }_l^{ p} | {\omega }_l^{ c} )={\delta }_{\omega } q_1 ( {\omega }_l^{ p} | {\omega }_l^{ c} )\\ \hspace{1em}+(1-{\delta }_{\omega }) q_2 ( {\omega }_l^{ p} | {\omega }_l^{ c} ),\hspace{1em}l=1,\dots ,m,\end{array}$$(A.1) where q₁ is defined in Equation (A.2), q₂ is the density of a univariate normal $\mathcal{N} ({\omega }_l^{ c},{\sigma }_{\omega }^2),{\delta }_{\omega }$ is a positive real number such that $0\le {\delta }_{\omega }\le 1$, and c and p refer to current and proposed values, respectively. According to Equation (A.1) we carry out with probability ${\delta }_{\omega }$ a M-H step with proposal distribution $q_1 ( {\omega }_l^{ p} | {\omega }_l^{ c} )$,(A.2) $$q_1 ( {\omega }_l^{ p} | {\omega }_l^{ c} )\propto \sum _{h = 0}^{\tilde{n}-1}{I}_h 1_{[ h/n \le {\omega }_l^{ p} < (h+1)/n ] },$$(A.2) where $\tilde{n}=\lfloor n/2\rfloor$ and I_h is the value of the squared modulus of the discrete Fourier transform of the observations y evaluated at frequency h/n $$I_h={\left| \sum _{j = a}^b{y}_j \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-i 2\pi \frac{h}{n} \right) \right|}^{ 2}.$$

In this way frequencies are proposed from regions in parameter space with high posterior density, yielding a Markov chain which converges quickly to its invariant distribution. The proposal distribution $q_1 ( {\omega }_l^{ p} | {\omega }_l^{ c} )$ is independent of the current state ${\omega }_l^{ c}$. The acceptance probability for this move is$$\alpha =\min \left\{1,\frac{\pi ({\mathit{\omega }}^{ p} | \mathit{\beta }, {\sigma }^2, m, \mathit{y})}{\pi ({\mathit{\omega }}^{ c} | \mathit{\beta }, {\sigma }^2, m, \mathit{y})}\times \frac{q_1 ( {\omega }_l^{ c} )}{q_1 ( {\omega }_l^{ p} )}\right\},$$where ${\mathit{\omega }}^{ p}=({\omega }_1^{ c},\dots ,{\omega }_{l-1}^{ c},{\omega }_l^{ p},{\omega }_{l+1}^{ c},\dots ,{\omega }_m^{ c})\prime$. On the other hand, with probability $1-{\delta }_{\omega }$, we perform a random walk M-H step with proposal distribution $q_2 ( {\omega }_l^{ p} | {\omega }_l^{ c} )$, whose density is a univariate normal distribution with mean ${\omega }_l^{ c}$ and variance ${\sigma }_{\omega }^2$, that is, ${\omega }_l^{ p} | {\omega }_l^{ c} \sim \mathcal{N}( {\omega }_l^{ c}, {\sigma }_{\omega }^2 )$. This perturbation around the current value ${\omega }_l^{ c}$ allows a local exploration of the conditional posterior distribution. The acceptance probability for this move is$$\alpha =\min \left\{1,\frac{\pi ({\mathit{\omega }}^{ p} | \mathit{\beta }, {\sigma }^2, m, \mathit{y})}{\pi ({\mathit{\omega }}^{ c} | \mathit{\beta }, {\sigma }^2, m, \mathit{y})}\right\}.$$

Setting ${\delta }_{\omega }$ to a relative low value integrates a fairly high acceptance rate with a quick exploration of the parameter space. For our experiments, we set ${\sigma }_{\omega }^2={(1/50n)}^2$ and ${\delta }_{\omega }=0.2$.

Sampling $\mathit{\beta }$: Given values of $\mathit{\omega }$ and ${\sigma }^2$, the vector of linear parameters $\mathit{\beta }$ can be sampled via a M-H step from the target posterior conditional distribution $\pi ( \mathit{\beta } | \mathit{\omega }, {\sigma }^2, m, \mathit{y})$ (see Equation (9)(9) $$\begin{array}{c}\pi ( \mathit{\beta } | \mathit{\omega }, {\sigma }^2, m, \mathit{y})\\ \hspace{1em}\propto \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }^2}\sum _{t = a}^b\{y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2-\frac{1}{2{\sigma }_{\beta }^2} \mathit{\beta }\prime \mathit{\beta }\right],\end{array}$$(9) ). The proposed vector of coefficients ${\mathit{\beta }}^{ p}$ can be drawn from normal approximations to their posterior conditional distribution (e.g., Gilks, Richardson, and Spiegelhalter 1995; Rosen, Wood, and Stoffer 2012),(A.3) $${\mathit{\beta }}^{ p}\sim {\mathcal{N}}_{2m+2} ( {\widehat{\mathit{\beta }}}^{ p}, {\widehat{\Sigma }}^{ p}),$$(A.3) where$${\widehat{\mathit{\beta }}}^{ p}=\underset{{\mathit{\beta }}^{ p}}{\text{arg} \text{max}} \pi ( {\mathit{\beta }}^{ p} | \mathit{\omega }, {\sigma }^2, m, \mathit{y}),$$ and$${\widehat{\Sigma }}^{ p}={\left\{-\frac{{\partial }^2\hspace{0.17em}\text{log}\hspace{0.17em} \pi ( {\mathit{\beta }}^{ p} | \mathit{\omega }, {\sigma }^2, m, \mathit{y})}{\partial {\mathit{\beta }}^{ p} \partial {\mathit{\beta }}^{ p \prime }}{|}_{{\mathit{\beta }}^{ p} = {\widehat{\mathit{\beta }}}^{ p}}\right\}}^{-1}.$$

The proposal for ${\mathit{\beta }}^{ p}$ is independent of the current values ${\mathit{\beta }}^{ c}$, and the acceptance probability for this move is$$\alpha =\min \left\{1,\frac{\pi ( {\mathit{\beta }}^{ p} | \mathit{\omega }, {\sigma }^2, m, \mathit{y})}{\pi ( {\mathit{\beta }}^{ c} | \mathit{\omega }, {\sigma }^2, m, \mathit{y})}\times \frac{q ( {\mathit{\beta }}^{ c} )}{q ( {\mathit{\beta }}^{ p} )}\right\},$$where $q ( {\mathit{\beta }}^{ c} )$ and $q ( {\mathit{\beta }}^{ p} )$ denote the normal proposal densities ${\mathcal{N}}_{2m+2} ( {\widehat{\mathit{\beta }}}^{ c}, {\widehat{\Sigma }}^{ c})$ and ${\mathcal{N}}_{2m+2} ( {\widehat{\mathit{\beta }}}^{ p}, {\widehat{\Sigma }}^{ p})$, respectively.

A.1.2 Between-Model Moves

The number of frequencies on a segment is proposed to either increase or decrease by one. Let ${\mathit{\theta }}^{ c}=( {\mathit{\beta }}^{ c \prime }, {\mathit{\omega }}^{ c \prime }, {\sigma }^{ 2c} )\prime$ and assume the Markov chain is currently at $( m^{ c}, {\mathit{\theta }}^{ c} )$. We propose a move to $( m^{ p}, {\mathit{\theta }}^{ p} )$ by drawing from a proposal density $q ( m^{ p}, {\mathit{\theta }}^{ p} | m^{ c}, {\mathit{\theta }}^{ c} )$ and accepting this update with probability(A.4) $$\begin{array}{c}\alpha =\min \{1,\frac{\mathit{L}( m^{ p}, {\mathit{\theta }}^{ p} | \mathit{y})}{\mathit{L}( m^{ c}, {\mathit{\theta }}^{ c} | \mathit{y})}\times \frac{\pi \left(m^{ p}\right) \pi ({\mathit{\theta }}^{ p} | m^{ p})}{\pi \left(m^{ c}\right) \pi ({\mathit{\theta }}^{ c} | m^{ c})}\\ \hspace{1em}\times \frac{q (m^{ c}, {\mathit{\theta }}^{ c} | m^{ p}, {\mathit{\theta }}^{ p} )}{q (m^{ p}, {\mathit{\theta }}^{ p} | m^{ c}, {\mathit{\theta }}^{ c} )}\}.\end{array}$$(A.4)

The proposed state $( m^{ p}, {\mathit{\theta }}^{ p} )$ is drawn by first drawing $m^{ p}$, followed by ${\mathit{\omega }}^{ p}$, ${\mathit{\beta }}^{ p}$ and ${\sigma }^{ 2p}$. In fact, we can rewrite the proposal density as$$\begin{array}{c}q (m^{ p}, {\mathit{\theta }}^{ p} | m^{ c}, {\mathit{\theta }}^{ c} )=q (m^{ p} | m^{ c})\times q ({\mathit{\theta }}^{ p} | m^{ p}, m^{ c}, {\mathit{\theta }}^{ c})\hfill \\ =q (m^{ p} | m^{ c})\times q ({\mathit{\omega }}^{ p} | m^{ p}, m^{ c}, {\mathit{\theta }}^{ c})\hfill \\ \hspace{1em}\times q ( {\mathit{\beta }}^{ p} | {\mathit{\omega }}^{ p}, m^{ p}, m^{ c}, {\mathit{\theta }}^{ c})\hfill \\ \hspace{1em}\times q ( {\sigma }^{ 2p} | {\mathit{\beta }}^{ p}, {\mathit{\omega }}^{ p}, m^{ p}, m^{ c}, {\mathit{\theta }}^{ c}).\hfill \end{array}$$

Birth move: If a birth move is proposed, we have that $m^{ p}=m^{ c}+1$. The proposed frequency vector ${\mathit{\omega }}^{ p}$ is constructed as$${\mathit{\omega }}^{ p}=({\omega }_1^c,\dots , {\omega }_{m^c}^c, {\omega }_{m^p}^{ p})\prime ,$$namely by keeping the current vector of frequencies and proposing an additional frequency ${\omega }_{m^p}^{ p}$. Alternatively to Andrieu and Doucet (1999), we choose to sample ${\omega }_{m^p}^{ p}$ uniformly on the interval $(0, {\varphi }_{\omega })$, where $0<{\varphi }_{\omega }<0.5$. The value of ${\varphi }_{\omega }$ can be chosen to reflect prior information about the significant frequencies that drive the variation in the data, for example, by choosing ${\varphi }_{\omega }$ in the low frequencies range ($0<{\varphi }_{\omega }<0.1$). Additionally, for computational and/or modelling reasons, we would like not to sample frequencies that are too close to each other. Hence, we choose to draw a candidate value ${\omega }_{m^p}^p$ uniformly from the union of intervals of the form $[{\omega }_l^c+{\psi }_{\omega }, {\omega }_{l+1}^c-{\psi }_{\omega }]$, for $l=0,\dots ,m_c$ and denoting ${\omega }_0^{ c}=0$ and ${\omega }_{ m^c+1}^c={\varphi }_{\omega }$. Here, ${\psi }_{\omega }$ is a fixed minimum distance between frequencies, which is chosen larger than $\frac{1}{n}$; in fact, when the separation of two frequencies is less than the Nyquist step (Priestley 1981), that is, $| {\omega }_{ l}-{\omega }_{ l+1} |<\frac{1}{n}$, the two frequencies are indistinguishable (Dou and Hodgson 1995). Moreover, we sort the proposed vector of frequencies ${\mathit{\omega }}^{ p}$ to ensure identifiability and perform practical estimation, as suggested in Andrieu and Doucet (1999). For proposed ${\mathit{\omega }}^{ p}$ and given ${\sigma }^{ 2c}$, the proposed vector of linear coefficients ${\mathit{\beta }}^{ p}$ is sampled following the same procedure of Section A.1.1, namely we draw ${\mathit{\beta }}^{ p}$ from a normal approximation to their posterior conditional distribution $\pi ( {\mathit{\beta }}^{ p} | {\mathit{\omega }}^{ p}, {\sigma }^{ 2c}, m^{ p}, \mathit{y} )$. Finally, the residual variance ${\sigma }^{ 2p}$ is sampled directly from its posterior conditional distribution $\pi ({\sigma }^{ 2p} | {\mathit{\omega }}^{ p}, {\mathit{\beta }}^{ p}, m^{ p}, \mathit{y})$ (see Equation (10)(10) $$\begin{array}{c}{\sigma }_{| \mathit{\omega }, \mathit{\beta }}^2\sim \text{Inverse}-\text{Gamma}\\ \hspace{1em}\left( \frac{n+{\nu }_0}{2}, \frac{{\gamma }_0+\sum _{t = a}^b\{ y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2}{2}\right).\end{array}$$(10) ). The proposed state $(m^p, {\mathit{\theta }}^{ p} )$ is accepted or reject in a M-H step with probability$$\begin{array}{c}\alpha =\min \{1, \frac{\mathit{L}( {\mathit{\theta }}^{ p}, m^{ p} | \mathit{y})}{\mathit{L}( {\mathit{\theta }}^{ c}, m^{ c} | \mathit{y})}\times \frac{\pi \left(m^{ p}\right) \pi ({\mathit{\theta }}^{ p} | m^{ p})}{\pi \left(m^{ c}\right) \pi ({\mathit{\theta }}^{ c} | m^{ c})}\\ \hspace{1em}\times \frac{d_{m^{ p}}·\left(\frac{1}{m^{ p}}\right)·q ( {\mathit{\beta }}^{ c} )·q ( {\sigma }^{ 2c} )}{b_{m^{ c}}·q \left({\omega }_{m^{ p}}^{ p}\right)·q ( {\mathit{\beta }}^{ p} )·q ( {\sigma }^{ 2p} )}\},\end{array}$$ where the likelihood function is given in Equation (5)(5) $$\begin{array}{c}\mathit{L}( m_j, {\mathit{\theta }}_{ j} | {\mathit{y}}_j )={( 2\pi {\sigma }_j^2 )}^{ -n_j/2}\\ \hspace{1em}\times \hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\frac{1}{2{\sigma }_j^2} \sum _{ t \in {\mathrm{I}}_j}^{}{\left\{ y_t-{\mathit{x}}_t \left({\mathit{\omega }}_j\right)\prime {\mathit{\beta }}_{ j} \right\}}^{ 2}\right],\end{array}$$(5) , $\pi ( m )$ is the density of the Poisson distribution truncated at $m_{\text{max}},b_{m^{ c}}$ and $d_{m^{ p}}$ are defined in Equation (7)(7) $$b_z=c \text{min}\{1,\frac{\pi (z+1)}{\pi \left(z\right)}\},\hspace{1em}{d}_{z+1}=c \text{min}\{1,\frac{\pi \left(z\right)}{\pi (z+1)}\},$$(7) , $q ({\omega }_{m^p}^p)$ is the density of the uniform proposal for sampling the additional frequency, $q ({\mathit{\beta }}^{ c})$ and $q ({\mathit{\beta }}^{ p})$ are the multivariate normal proposal densities ${\mathcal{N}}_{2m^{ c}+2} ( {\widehat{\mathit{\beta }}}^{ c}, {\widehat{\Sigma }}^{ c})$ and ${\mathcal{N}}_{2m^{ p}+2} ( {\widehat{\mathit{\beta }}}^{ p}, {\widehat{\Sigma }}^{ p})$, respectively; $q ( {\sigma }_p^2 )$ and $q ( {\sigma }_c^2 )$ are the Inverse-Gamma proposal densities defined in Equation (10)(10) $$\begin{array}{c}{\sigma }_{| \mathit{\omega }, \mathit{\beta }}^2\sim \text{Inverse}-\text{Gamma}\\ \hspace{1em}\left( \frac{n+{\nu }_0}{2}, \frac{{\gamma }_0+\sum _{t = a}^b\{ y_t-{\mathit{x}}_t ( \mathit{\omega } )\prime \mathit{\beta } {\}}^2}{2}\right).\end{array}$$(10) .

Death move: If a death move is proposed, then $m^{ p}=m^{ c}-1$. A vector of frequencies ${\mathit{\omega }}^{ p}$ is constructed by randomly selecting with probability $\frac{1}{m^{ c}}$ one of the current frequencies as the candidate frequency for removal. Given ${\mathit{\omega }}^{ p}$ and ${\sigma }^{ 2c}$, a vector of linear coefficients ${\mathit{\beta }}^{ p}$ is sampled from a normal approximation to its posterior conditional distribution. Finally, conditioned on ${\mathit{\omega }}^{ p}$ and ${\mathit{\beta }}^{ p}$, the residual variance ${\sigma }^{ 2p}$ is drawn from its posterior Inverse-Gamma distribution. It is straightforward to see that the acceptance probability for the death step has the same form as the birth step, with the proper change of labelling of the variables, and the ratio terms inverted.

A.2 Updating the Change-Point Model

A.2.1 Within-Model Move

Let ${\mathit{s}}_{ (k)}^{ c}=(s_1^{ c},\dots ,s_k^{ c})\prime$ be the current vector of change-points locations, ${\mathit{m}}_{ (k)}^{ c}=(m_1^{ c},\dots ,m_{k+1}^{ c})\prime$ be the current vector of number of frequencies, ${\mathit{\omega }}_{ (k)}^{ c}=({\mathit{\omega }}_1^{ c \prime },\dots ,{\mathit{\omega }}_{k+1}^{ c \prime })\prime$ be the current vector of frequencies. Let ${\mathit{\beta }}_{ (k)}^{ c}=( {\mathit{\beta }}_1^{ c \prime },\dots ,{\mathit{\beta }}_{k+1}^{ c \prime })\prime$ and ${\mathit{\sigma }}_{(k)}^{ 2}=( {\sigma }_1^{ 2c},\dots ,{\sigma }_{ k+1}^{ 2c})\prime$ be the current vectors of linear coefficients and residual variances, respectively.

Let us also define ${\mathit{\theta }}_{ (k)}^{ c}=({\mathit{\beta }}_{ (k)}^{ c \prime }, {\mathit{\omega }}_{ (k)}^{ c \prime }, {\mathit{\sigma }}_{(k)}^{ 2c \prime })\prime .$ Following Green (1995), a change-point, $s_j^{ c}$ say, is randomly selected with probability $\frac{1}{k}$ from the existing set of change-points. To explore the parameter space in an efficient way and similar to above we construct a mixture distribution $q ( s_j^{ p} | s_j^{ c} )$, as(A.5) $$q ( s_j^{ p} | s_j^{ c} )={\delta }_s q_1 ( s_j^{ p} | s_j^{ c} )+(1-{\delta }_s) q_2 ( s_j^{ p} | s_j^{ c} ),$$(A.5) where q₁ is the density of a Uniform $[s_{j-1}^{ c}+\psi ,s_{j+1}^{ c}-\psi ]$, q₂ is the density of a univariate normal $\mathcal{N} (s_j^{ c},{\sigma }_s^2)$ and δ_s is a positive real number such that $0\le {\delta }_s\le 1$. We propose with probability δ_s a candidate value $s_j^{ p}$ from the above uniform distribution where ψ is a fixed minimum time between change-points avoiding change-points being too close to each other. On the other hand, with probability $(1-{\delta }_s),s_j^{ p}$ arises as a normal random walk proposal centered at the current change-point $s_j^{ c}$. The proposed vector of change-points locations is denoted by$${\mathit{s}}_{ \left(k\right)}^{ p}=(s_1^{ c},\dots ,s_{j-1}^{ c},s_j^{ p},s_{j+1}^{ c},\dots ,s_k^{ c})\prime ,$$ and hence the proposed value $s_j^{ p}$ induces a new proposed data partition on $[s_{j-1}^{ c},s_{j+1}^{ c}]$ corresponding to $[s_{j-1}^{ c},s_j^{ p})$ and $[s_j^{ p},s_{j+1}^{ c})$. We denote the vectors of observations belonging to these two proposed segments as ${\mathit{y}}_j^{ p}$ and ${\mathit{y}}_{j+1}^{ p}$, which include $n_j^{ p}$ and $n_{j+1}^{ p}$ observations, respectively. Given ${\mathit{s}}_{ (k)}^{ p}$, the proposed number of frequencies $m_j^{ p},m_{j+1}^{ p}$ are set equal to the current ones $m_j^{ c},m_{j+1}^{ c}$, so that ${\mathit{m}}_{ (k)}^{ p}={\mathit{m}}_{ (k)}^{ c}$. Similarly, the proposed pair of frequency vectors ${\mathit{\omega }}_j^{ p},{\mathit{\omega }}_{j+1}^{ p}$ is chosen equal to the current pair ${\mathit{\omega }}_j^{ c},{\mathit{\omega }}_{j+1}^{ c}$ in the corresponding segments, that is, ${\mathit{\omega }}_{ (k)}^{ p}={\mathit{\omega }}_{ (k)}^{ c}$. The proposed vectors ${\mathit{\beta }}_j^{ p},{\mathit{\beta }}_{j+1}^{ p}$ are sampled from normal approximations to their posterior conditional distributions $\pi ( {\mathit{\beta }}_j^{ p} | {\mathit{\omega }}_j^{ p}, {\sigma }_{ j}^{ 2c}, m_j^{ p}, {\mathit{y}}_j^{ p} )$ and $\pi ( {\mathit{\beta }}_{j+1}^{ p} | {\mathit{\omega }}_{j+1}^{ p}, {\sigma }_{ j+1}^{ 2c}, m_{j+1}^{ p}, {\mathit{y}}_{j+1}^{ p} )$ and are accepted in a M-H step with probability$$\begin{array}{c}\alpha =\min \{ 1,\frac{\mathit{L}( k, {\mathit{m}}_{ \left(k\right)}^{ p}, {\mathit{s}}_{ \left(k\right)}^{ p}, {\mathit{\theta }}_{ \left(k\right)}^{ p} | \mathit{y} )}{\mathit{L}( k, {\mathit{m}}_{ \left(k\right)}^{ c}, {\mathit{s}}_{ \left(k\right)}^{ c}, {\mathit{\theta }}_{ \left(k\right)}^{ c} | \mathit{y} )}\\ \hspace{1em}\times \frac{\pi ( {\mathit{s}}_{ \left(k\right)}^{ p} | k) \pi ( {\mathit{\theta }}_{ \left(k\right)}^{ p} | {\mathit{m}}_{ \left(k\right)}^{ p}, k)}{\pi ( {\mathit{s}}_{ \left(k\right)}^{ c} | k) \pi ( {\mathit{\theta }}_{ \left(k\right)}^{ c} | {\mathit{m}}_{ \left(k\right)}^{ c}, k )}\times \frac{\prod _{h=j}^{j+1}q ( {\mathit{\beta }}_h^{ c} ) }{\prod _{h=j}^{j+1}q ( {\mathit{\beta }}_h^{ p} ) } \},\end{array}$$ where the likelihood is specified in Equation (4)(4) $$\begin{array}{c}\mathit{L}( k, {\mathit{m}}_{ \left(k\right)}, {\mathit{s}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)}, | \mathit{y} )=\prod _{j=1}^{k+1}\mathit{L}( m_j, {\mathit{\theta }}_{ j} | {\mathit{y}}_j ),\\ {\mathit{y}}_j=( y_t: t \in {\mathrm{I}}_j ),\end{array}$$(4) and $q ( {\mathit{\beta }}_h^{ c} )$ and $q ( {\mathit{\beta }}_h^{ p} )$ are multivariate Gaussian as in Equation (A.3). Note that the likelihood ratio and the prior ratio differ from one only for the two segments affected by the move of the change-points. Finally, the residual variances ${\sigma }_j^{ 2p}$, ${\sigma }_{j+1}^{ 2p}$ are drawn from their posterior conditional distributions $\pi ({\sigma }_j^{ 2p} | {\mathit{\omega }}_j^{ p}, {\mathit{\beta }}_j^{ p}, m_j^{ p}, {\mathit{y}}_j^{ p}),\pi ({\sigma }_{j+1}^{ 2p} | {\mathit{\omega }}_{j+1}^{ p}, {\mathit{\beta }}_{j+1}^{ p}, m_{j+1}^{ p}, {\mathit{y}}_{j+1}^{ p})$ in a Gibbs step.

A.2.2 Between-Model Moves

Let ${\mathit{\xi }}_{ (k^{ c})}^{ c}={( {\mathit{s}}_{ (k^{ c})}^{ c \prime }, {\mathit{m}}_{ (k^{ c})}^{ c \prime }, {\mathit{\theta }}_{ (k^{ c})}^{ c \prime } )}^{ \prime }$ and assume the Markov chain is at $( k^{ c}, {\mathit{\xi }}_{ (k^{ c})}^{ c} )$. We propose a move to $( k^{ p}, {\mathit{\xi }}_{ (k^{ p})}^{ p} )$ by first drawing $k^{ p}$, followed by sampling the change-point locations ${\mathit{s}}_{ (k^{ p})}^{ p}$. The latter involves either merging two segments (death) or splitting a segment (birth). The number of frequencies and their values in the proposed segments are selected from the current state. We draw ${\mathit{\beta }}_{ (k^{ p})}^{ p}$ and jointly update the entire state $( k^{ p}, {\mathit{\xi }}_{ (k^{ p})}^{ p} )$. Hence, we propose a move to $( k^{ p}, {\mathit{\xi }}_{ (k^{ p})}^{ p} )$ by drawing from a proposal density of the form$$\begin{array}{c}q ( k^{ p}, {\mathit{\xi }}_{ \left(k^{ p}\right)}^{ p} | k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} )\hfill \\ \hspace{1em}=q ( k^{ p} | k^{ c} )\times q ( {\mathit{\xi }}_{ \left(k^{ p}\right)}^{ p} | k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} )\hfill \\ \hspace{1em}=q ( k^{ p} | k^{ c} ) \times q ( {\mathit{s}}_{ \left(k^{ p}\right)}^{ p} | k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} )\hfill \\ \hspace{1em}\times q ( {\mathit{m}}_{ \left(k^{ p}\right)}^{ p}, {\mathit{\omega }}_{ \left(k^{ p}\right)}^{ p} | {\mathit{s}}_{ \left(k^{ p}\right)}^{ p}, k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} )\hfill \\ \hspace{1em}\times q ({\mathit{\sigma }}_{ \left(k^{ p}\right)}^{ 2p} | {\mathit{m}}_{ \left(k^{ p}\right)}^{ p}, {\mathit{\omega }}_{ \left(k^{ p}\right)}^{ p}, {\mathit{s}}_{ \left(k^{ p}\right)}^{ p}, k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} )\hfill \\ \hspace{1em}\times q ( {\mathit{\beta }}_{ \left(k^{ p}\right)}^{ p} | {\mathit{\sigma }}_{ \left(k^{ p}\right)}^{ 2p}, {\mathit{m}}_{ \left(k^{ p}\right)}^{ p}, {\mathit{\omega }}_{ \left(k^{ p}\right)}^{ p}, {\mathit{s}}_{ \left(k^{ p}\right)}^{ p}, k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} ).\hfill \end{array}$$

Birth move: If a birth move is proposed, we have that $k^{ p}=k^{ c}+1$. We draw a new change-point uniformly on $f ({\mathit{s}}_{ (k^{ c})}^{ c}, {\psi }_s)$, the support of ${\mathit{s}}_{ (k^{ c})}^{ c}$ given the constraints imposed by ψ_s, that is, $f ({\mathit{s}}_{ (k^{ c})}^{ c}, {\psi }_s)=[1+{\psi }_s, s_1^{ c}-{\psi }_s] \cup [s_1^{ c}+{\psi }_s, s_2^{ c}-{\psi }_s]\cup \cdots \cup [s_{k^{ c}}^{ c}+{\psi }_s, n-{\psi }_s].$ Hence, the new proposed location ${\tilde{s}}_j$ is sampled from a $\text{Uniform} \{ f ({\mathit{s}}_{ (k^{ c})}, {\psi }_s)\}$, where the proposal density is given by(A.6) $$q ({\mathit{s}}_{ \left(k^{ p}\right)}^{ p} | k^{ p}, k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c})=\frac{1}{(n-2{\psi }_s (k^{ c}+1)-1)}.$$(A.6)

As the proposed location ${\tilde{s}}_j$ will lie within an existing interval $(s_j^{ c}, s_{j+1}^{ c})$ with probability one, we can define the proposed change-points location vector as$${\mathit{s}}_{ \left(k^{ p}\right)}^{ p}=(s_1^{ c},\dots ,s_j^{ c},{\tilde{s}}_j,s_{j+1}^{ c},\dots ,s_{k^{ c}}^{ c})\prime .$$

The number of frequencies $m_j^{ p},m_{j+1}^{ p}$ corresponding to the two newly proposed segments $[s_j^{ c},{\tilde{s}}_j)$ and $[{\tilde{s}}_j,s_{j+1}^{ c})$ are set equal to the current number of frequencies on the whole segment $(s_j^{ c}, s_{j+1}^{ c})$. Therefore, we can construct the proposed vector of the number of frequencies ${\mathit{m}}_{ (k^{ p})}^{ p}$ and the proposed vector of frequencies ${\mathit{\omega }}_{ (k^{ p})}^{ p}$ as$$\begin{array}{c}{\mathit{m}}_{ \left(k^{ p}\right)}^{ p}=( m_1^{ c},\dots ,m_{j-1}^{ c},m_j^{ c},m_j^{ c},m_{j+1}^{ c},\dots ,m_{k^{ c}+1}^{ c})\prime ,\hfill \\ {\mathit{\omega }}_{ \left(k^{ p}\right)}^{ p}=({\mathit{\omega }}_1^{ c \prime },\dots ,{\mathit{\omega }}_{j-1}^{ c \prime },{\mathit{\omega }}_j^{ c \prime },{\mathit{\omega }}_j^{ c \prime },{\mathit{\omega }}_{j+1}^{ c \prime },\dots ,{\mathit{\omega }}_{k^{ c}+1}^{ c \prime })\prime .\hfill \end{array}$$

The proposed vector of residual variances ${\mathit{\sigma }}_{ (k^{ p})}^{ 2p}$ is$${\mathit{\sigma }}_{ \left(k^{ p}\right)}^{ 2p}={( {\sigma }_1^{ 2c},\dots ,{\sigma }_{j-1}^{ 2c},{\sigma }_j^{ 2p},{\sigma }_{j+1}^{ 2p},{\sigma }_{j+1}^{ 2c},\dots ,{\sigma }_{k^{ c}+1}^{ 2c} )}^{ \prime },$$where the residual variances ${\sigma }_j^{ 2p},{\sigma }_{j+1}^{ 2p}$ for the split partition are constructed following Green (1995), namely as a perturbation of the current variance ${\sigma }_j^{ 2c}$. Specifically, we draw $u\sim \text{Uniform} (0,1)$ and let ${\sigma }_j^{ 2p},{\sigma }_{j+1}^{ 2p}$ be deterministic transformations of ${\sigma }_j^{ 2c}$, that is(A.7) $${\sigma }_j^{ 2p}=\frac{u}{1-u}{\sigma }_j^{ 2c},\hspace{1em}{\sigma }_{j+1}^{ 2p}=\frac{1-u}{u}{\sigma }_j^{ 2c}.$$(A.7)

Finally, the proposed vector of linear coefficients ${\mathit{\beta }}_{ (k^{ p})}^{ p}$ is$${\mathit{\beta }}_{ \left(k^{ p}\right)}^{ p}={( {\mathit{\beta }}_1^{ c \prime },\dots ,{\mathit{\beta }}_{j-1}^{ c \prime },{\mathit{\beta }}_j^{ p \prime },{\mathit{\beta }}_{j+1}^{ p \prime },{\mathit{\beta }}_{j+1}^{ c \prime },\dots ,{\mathit{\beta }}_{k^{ c}+1}^{ c \prime } )}^{ \prime },$$where the vectors ${\mathit{\beta }}_j^{ p},{\mathit{\beta }}_{j+1}^{ p}$ are drawn from normal approximations of their posterior conditional distribution, as in Section A.1.1. The proposed move to the state $(k^{ p}, {\mathit{\xi }}_{ (k^{ p})}^{ p})$ is accepted with probability$$\begin{array}{c}\alpha =\min \{ 1,\frac{\mathit{L}( k^{ p}, {\mathit{\xi }}_{ \left(k^{ p}\right)}^{ p} | \mathit{y} )}{\mathit{L}( k^{ c}, {\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} | \mathit{y} )}\times \frac{\pi \left(k^{ p}\right) \pi ({\mathit{\xi }}_{ \left(k^{ p}\right)}^{ p} | k^{ p})}{\pi \left(k^{ c}\right) \pi ({\mathit{\xi }}_{ \left(k^{ c}\right)}^{ c} | k^{ c})}\\ \hspace{1em}\times \frac{d_{k^{ p}}·\frac{1}{k^{ p}}·\frac{1}{2}·q \left({\mathit{\beta }}_j^{ c}\right)}{b_{k^{ c}}·q \left({\mathit{s}}_{ \left(k^{ p}\right)}^{ p}\right)·\prod _{h=j}^{j+1}q \left({\mathit{\beta }}_h^{ p}\right)}\times {\mathit{J}}_{{\sigma }^{ 2}} \},\end{array}$$ where the likelihood function is provided in Equation (4)(4) $$\begin{array}{c}\mathit{L}( k, {\mathit{m}}_{ \left(k\right)}, {\mathit{s}}_{ \left(k\right)}, {\mathit{\theta }}_{ \left(k\right)}, | \mathit{y} )=\prod _{j=1}^{k+1}\mathit{L}( m_j, {\mathit{\theta }}_{ j} | {\mathit{y}}_j ),\\ {\mathit{y}}_j=( y_t: t \in {\mathrm{I}}_j ),\end{array}$$(4) , $q ({\mathit{s}}_{ (k^{ p})}^{ p})$ is the uniform density defined in Equation (A.6); $q ({\mathit{\beta }}_h^{ c})$, $q ({\mathit{\beta }}_h^{ p})$ are the multivariate normal proposal densities, and the Jacobian ${\mathit{J}}_{{\sigma }^{ 2}}$ is$${\mathit{J}}_{{\sigma }^{ 2}}=| \frac{\partial ({\sigma }_j^{ 2p},{\sigma }_{j+1}^{ 2p})}{\partial ({\sigma }_j^{ 2c},u)}|=2 {(\sqrt{{\sigma }_j^{ 2p}}+\sqrt{{\sigma }_{j+1}^{ 2p}} )}^{ 2}.$$

The numerator of the proposal ratio is better understood by looking at the details of the death step, which are given below.

Death move: If a death step is proposed, then $k^{ p}=k^{ c}-1$. A candidate change-point $s_j^{ c}$ to be removed is sampled uniformly from the vector of existing change-points; that is, we propose to remove $s_j^{ c}$ with probability $\frac{1}{k^{ c}}$. Then, the proposed vector of change-points locations ${\mathit{s}}_{ (k)}^{ p}$ is defined as$${\mathit{s}}_{ \left(k\right)}^{ p}={(s_1^{ c},\dots ,s_{j-1}^{ c},s_{j+1}^{ c},\dots ,s_{k^{ c}}^{ c})}^{ \prime }.$$

The number $m_j^{ p}$ and the vector of relevant frequencies ${\mathit{\omega }}_j^{ p}$ of the newly merged segment $[s_{j-1}^{ c},s_{j+1}^{ c})$ are selected by drawing an index at random from $\{j,j+1\}$, obtaining say $j^{ *}$, and setting the proposed parameters equal to the current ones relative to the selected index. That is, we set $m_j^{ p}=m_{j^{ *}}^{ c}$ and ${\mathit{\omega }}_j^{ p}={\mathit{\omega }}_{j^{ *}}^{ c}$. Hence, the proposed vectors of number of frequencies ${\mathit{m}}_{ (k^{ p})}^{ p}$ and their values ${\mathit{\omega }}_{ (k^{ p})}^{ p}$ are constructed as follows$$\begin{array}{c}{\mathit{m}}_{ \left(k^{ p}\right)}^{ p}={( m_1^{ c},\dots ,m_{j-1}^{ c},m_j^{ p},m_{j+2}^{ c},\dots ,m_{k^{ c}+1}^{ c} )}^{ \prime },\hfill \\ {\mathit{\omega }}_{ (k^{ p)}}^{ p}={( {\mathit{\omega }}_1^{ c \prime },\dots ,{\mathit{\omega }}_{j-1}^{ c \prime },{\mathit{\omega }}_j^{ p \prime },{\mathit{\omega }}_{j+2}^{ c \prime },\dots ,{\mathit{\omega }}_{k^{ c}+1}^{ c \prime })}^{ \prime }.\hfill \end{array}$$

The residual variance ${\sigma }_j^{ 2p}$ of the newly merged segment is obtained by inverting the transformation of Equation (A.7). Specifically, we construct ${\sigma }_j^{ 2p}=\sqrt{{\sigma }_j^{ 2c} {\sigma }_{j+1}^{ 2c}}$ and set the proposed vector of residual variances ${\mathit{\sigma }}_{ (k^{ p})}^{ 2p}$ as$${\mathit{\sigma }}_{ \left(k^{ p}\right)}^{ 2p}={( {\sigma }_1^{ 2c},\dots ,{\sigma }_{j-1}^{ 2c},{\sigma }_j^{ 2p},{\sigma }_{j+2}^{ 2c},\dots ,{\sigma }_{k^{ c}+1}^{ 2c} )}^{ \prime }.$$

The proposed vector of linear coefficients ${\mathit{\beta }}_{ (k^{ p})}^{ p}$ is$${\mathit{\beta }}_{ \left(k^{ p}\right)}^{ p}={( {\mathit{\beta }}_1^{ c \prime },\dots ,{\mathit{\beta }}_{j-1}^{ c \prime },{\mathit{\beta }}_j^{ p \prime },{\mathit{\beta }}_{j+2}^{ c \prime },\dots ,{\mathit{\beta }}_{k^{ c}+1}^{ c \prime })}^{ \prime },$$where the vector of coefficients ${\mathit{\beta }}_j^{ p}$ is drawn from normal approximation to its posterior conditional distribution. The acceptance probability for the death step has the same form of the birth step, with the proper change of labelling of the variables, and the ratio terms inverted.