Efficient Bayesian PARCOR Approaches for Dynamic Modeling of Multivariate Time Series

TL;DR

This paper introduces a Bayesian lattice filtering approach for efficient, interpretable modeling of multivariate non-stationary time series, enabling accurate inference of spectral densities and time-frequency relationships.

Contribution

It proposes a novel, lower-dimensional Bayesian PARCOR framework that improves computational efficiency and accuracy in multivariate time series analysis compared to traditional TV-VAR models.

Findings

Effective in neuroscience and environmental data analysis

Achieves good fit in time-frequency domain

Provides accurate short-term forecasts

Abstract

A Bayesian lattice filtering and smoothing approach is proposed for fast and accurate modeling and inference in multivariate non-stationary time series. This approach offers computational feasibility and interpretable time-frequency analysis in the multivariate context. The proposed framework allows us to obtain posterior estimates of the time-varying spectral densities of individual time series components, as well as posterior measurements of the time-frequency relationships across multiple components, such as time-varying coherence and partial coherence. The proposed formulation considers multivariate dynamic linear models (MDLMs) on the forward and backward time-varying partial autocorrelation coefficients (TV-VPARCOR). Computationally expensive schemes for posterior inference on the multivariate dynamic PARCOR model are avoided using approximations in the MDLM context. Approximate inference on the corresponding time-varying vector autoregressive (TV-VAR) coefficients is obtained via Whittle's algorithm. A key aspect of the proposed TV-VPARCOR representations is that they are of lower dimension, and therefore more efficient, than TV-VAR representations. The performance of the TV-VPARCOR models is illustrated in simulation studies and in the analysis of multivariate non-stationary temporal data arising in neuroscience and environmental applications. Model performance is evaluated using goodness-of-fit measurements in the time-frequency domain and also by assessing the quality of short-term forecasting.

Figures

Tables

Table 1. Table 1 : Mean ASE values and corresponding standard deviations (in parentheses) obtained from TV-VPARCOR and TV-VAR ( 2 ) 2 (2) models for the TV-VAR ( 2 ) 2 (2) simulated data.

	Case (i): ${\varphi }_{1,1,2}=0$
$Model$	$g_{11}$	$g_{22}$	${\rho }_{12}^2$
PARCOR	$0.0246\left(0.0183\right)$	$0.0255\left(0.0147\right)$	$0.0008\left(0.0006\right)$
TV-VAR(2)	$0.0171\left(0.0068\right)$	$0.0186\left(0.0080\right)$	$0.0009\left(0.0005\right)$
	Case (ii): ${\varphi }_{1,1,2}=-0.8$
$Model$	$g_{11}$	$g_{22}$	${\rho }_{12}^2$
PARCOR	$0.0284\left(0.0118\right)$	$0.0238\left(0.0086\right)$	$0.0027\left(0.0023\right)$
TV-VAR(2)	$0.0254\left(0.0073\right)$	$0.0253\left(0.0081\right)$	$0.0023\left(0.0011\right)$

Table 2. Table 2: Computation times (in seconds) for TV-VPARCOR and TV-VAR models.

Model	(i): ${\varphi }_{1,1,2}=0$	(ii): ${\varphi }_{1,1,2}=-0.8$
PARCOR	$558s$	$521s$
TV-VAR	$925s$	$828s$

Equations

$$\displaystyle\mbox{\boldmath$x$}_{t}=\mbox{\boldmath$A$}_{t,1}^{(P)}\mbox{\boldmath$x$}_{t-1}+\dots+\mbox{\boldmath$A$}_{t,P}^{(P)}\mbox{\boldmath$x$}_{t-P}+\mbox{\boldmath$\epsilon$}_{t},\quad\mbox{\boldmath$\epsilon$}_{t}\sim\mathcal{N}(\mbox{\boldmath$0$},\bm{\Sigma}_{t}),$$

$$\displaystyle\mbox{\boldmath$f$}^{(P)}_{t}=\mbox{\boldmath$x$}_{t}-\sum^{P}_{j=1}\bm{A}_{t,j}^{(P)}\mbox{\boldmath$x$}_{t-j},\;\;\;\mbox{and}\;\;\;\mbox{\boldmath$b$}^{(P)}_{t}=\mbox{\boldmath$x$}_{t}-\sum^{P}_{j=1}\mbox{\boldmath$D$}_{t,j}^{(P)}\mbox{\boldmath$x$}_{t+j}.$$

$$\displaystyle\mbox{\boldmath$f$}_{t}^{(m-1)}$$

$$\displaystyle\mbox{\boldmath$b$}_{t}^{(m-1)}$$

$$\displaystyle\mbox{\boldmath$A$}_{t,j}^{(m)}$$

$$\displaystyle\mbox{\boldmath$D$}_{t,j}^{(m)}$$

$$\displaystyle\mbox{\boldmath$f$}_{t}^{(m-1)}$$

$$\displaystyle\mbox{\boldmath$b$}_{t}^{(m-1)}$$

$$\text{mbox}vec({\mathbf{\Lambda }}_{t,m}^{(m)})$$

$$\text{mbox}vec({\mathbf{\Theta }}_{t,m}^{(m)})$$

$$\text{mbox}vec({\mathbf{\Lambda }}_{0,m}^{(m)})|{\mathcal{D}}_{0,f,m}$$

$$\text{mbox}vec({\mathbf{\Theta }}_{0,m}^{(m)})|{\mathcal{D}}_{0,b,m}$$

$${\mathbf{\Delta }}_{f,m}=diag({\mathbf{\delta }}_{f,m,1}^{-1/2},\dots ,{\mathbf{\delta }}_{f,m,K}^{-1/2}),\text{mbox}and{\mathbf{\Delta }}_{b,m}=diag({\mathbf{\delta }}_{b,m,1}^{-1/2},\dots ,{\mathbf{\delta }}_{b,m,K}^{-1/2}),$$

$$\displaystyle\mbox{\boldmath$S$}_{t,f,m}=\frac{1}{(n_{0,f,m}+t)}\left(n_{0,f,m}\mbox{\boldmath$S$}_{0,f,m}+\sum^{t}_{i=1}\mbox{\boldmath$S$}_{i-1,f,m}^{1/2}\bm{Q}_{i,f,m}^{-1/2}\mbox{\boldmath$e$}_{i,f,m}\mbox{\boldmath$e$}_{i,f,m}^{\prime}\bm{Q}_{i,f,m}^{-1/2}\mbox{\boldmath$S$}_{i-1,f,m}^{1/2}\right),$$

$$E(\mbox{\boldmath$f$}^{(m-1)}_{t}|\mathcal{D}_{t-1,f,m})\approx\mbox{\boldmath$B$}_{t-m}^{(m-1)}\mbox{\boldmath$m$}_{t-1,f,m}.$$

$$V(\mbox{\boldmath$f$}^{(m-1)}_{t}|\mathcal{D}_{t-1,f,m})\approx\bm{Q}_{t,f,m}=\mbox{\boldmath$B$}_{t-m}^{(m-1)}\mbox{\boldmath$R$}_{t,f,m}(\mbox{\boldmath$B$}_{t-m}^{(m-1)})^{\prime}+\mbox{\boldmath$S$}_{t-1,f,m},$$

$$\displaystyle\mbox{\boldmath$m$}_{t,f,m}$$

$$\displaystyle\mbox{\boldmath$C$}_{t,f,m}$$

$${\mathbf{U}}_{t,f,m}$$

$$\text{mbox}vec({\mathbf{\Lambda }}_{t,m}^{(m)})|{\mathcal{D}}_T\approx \mathcal{N}({\mathbf{a}}_{T,f,m}(t-T),{\mathbf{R}}_{T,f,m}(t-T)),$$

$${\mathbf{a}}_{T,f,m}(t-T)$$

$${\mathbf{R}}_{T,f,m}(t-T)$$

$$DIC$$

$$p_{DIC}$$

$${\hat{p}}_{DIC}$$

$${\hat{p}}_{DIC}^m=\sum _{l=1}^m{\hat{p}}_{DIC,l}^l.$$

$$\hat{\mathbf{g}}(t,\omega )={\hat{\mathbf{\Phi }}}^{-1}(t,\omega )\times \hat{\mathbf{\Sigma }}\times {\hat{\mathbf{\Phi }}}^{\ast }(t,\omega {)}^{-1},$$

$${\hat{\rho }}_{i,j}^2(t,\omega )=|{\hat{g}}_{i,j}(t,\omega ){|}^2/\{{\hat{g}}_{i,i}(t,\omega ){\hat{g}}_{j,j}(t,\omega )\},$$

$${\hat{\gamma }}_{i,j}^2(t,\omega )=|{\hat{c}}_{i,j}(t,\omega ){|}^2/\{{\hat{c}}_{i,i}(t,\omega ){\hat{c}}_{j,j}(t,\omega )\}.$$

$${\mathbf{m}}_{T,f,m}(h)={\mathbf{m}}_{T,f,m};{\mathbf{C}}_{T,f,m}(h)={\mathbf{C}}_{T,f,m}+h\cdot {\mathbf{W}}_{T+1,f,m},$$

$${\hat{\mathbf{x}}}_{T+h}=\sum _{i=1}^P{\hat{\mathbf{A}}}_{T+h,i}^{(P)}{\hat{\mathbf{x}}}_{T+h-i}+{\widehat{\mathbf{ϵ}}}_{T+h}^{(P)},{\widehat{\mathbf{ϵ}}}_{T+h}^{(P)}\sim \mathcal{N}(0,{\mathbf{S}}_{T,f,P}).$$

$${\mathbf{x}}_t={\mathbf{\Phi }}_{1,t}{\mathbf{x}}_{t-1}+{\mathbf{\Phi }}_{2,t}{\mathbf{x}}_{t-2}+{\mathbf{ϵ}}_t,{\mathbf{ϵ}}_t\sim \mathcal{N}(0,{\mathbf{I}}_2),$$

$${\mathbf{\Phi }}_{1,t}=\left(\begin{array}{cc}{r}_{1,t}\cos (\frac{2\pi }{{\lambda }_{1,t}})& {\varphi }_{1,1,2}\\ 0& r_{2,t}\cos (\frac{2\pi }{{\lambda }_{2,t}})\end{array}\right)\text{mbox}and{\mathbf{\Phi }}_{2,t}=\left(\begin{array}{cc}-r_{1,t}^2& 0\\ 0& -r_{2,t}^2\end{array}\right),$$

$$\mathbf{g}(t,\omega )={\mathbf{\Phi }}^{-1}(t,\omega )\times \mathbf{\Sigma }\times {\mathbf{\Phi }}^{\ast }(t,\omega {)}^{-1},$$

Full text

The proposed formulation considers multivariate dynamic linear models (MDLMs) on the forward and backward time-varying partial autocorrelation coefficients (TV-VPARCOR). Computationally expensive schemes for posterior inference on the multivariate dynamic PARCOR model are avoided using approximations in the MDLM context. Approximate inference on the corresponding time-varying vector autoregressive (TV-VAR) coefficients is obtained via Whittle’s algorithm. A key aspect of the proposed TV-VPARCOR representations is that they are of lower dimension, and therefore more efficient, than TV-VAR representations. The performance of the TV-VPARCOR models is illustrated in simulation studies and in the analysis of multivariate non-stationary temporal data arising in neuroscience and environmental applications. Model performance is evaluated using goodness-of-fit measurements in the time-frequency domain and also by assessing the quality of short-term forecasting.

Keywords: Multivariate time series, Bayesian dynamic linear models, time-varying partial autocorrelations, time-varying vector autoregressions.

1. Introduction

Recent technological advances in scientific areas have led to multi-dimensional datasets with a complex temporal structure, often consisting of several time series components that are related over time. Inferring changes in the spectral content of each component, as well as time-varying relationships across components, is often relevant in applied areas. For example, understanding the interplay across temporal components derived from multi-channel/multi-location brain signals and brain imaging data is a key feature in brain connectivity studies (e.g., Yu et al., 2016; Chiang et al., 2017; Ting et al., 2017; Schmidt et al., 2016, among others). Multivariate time series analysis is also important for filtering, smoothing and prediction in environmental studies and finance where many variables are simultaneously measured over time (e.g., Zhang, 2017; Tsay, 2013).

Several time-domain, frequency-domain and time-frequency approaches are available for modeling and inferring spectral characteristics of univariate non-stationary time series. However, a much more limited number of approaches are available for computationally efficient and scientifically interpretable analysis of multivariate non-stationary time series. Furthermore, currently available statistical tools have important practical limitations. For instance, vector autoregressions (VARs) are often used in the analysis of multi-channel EEG (electroencephalogram) data, however, these models cannot capture the time-varying characteristics of these data. Alternative approaches incorporate more flexible and realistic dynamic structures (e.g., West et al., 1999; Prado et al., 2001; Nakajima and West, 2017), but are either not available for multivariate time series, or they are highly computationally intensive, requiring Markov chain Monte Carlo (MCMC) sampling for posterior inference. Frequency-domain and time-frequency approaches have also been developed, but typically these methods are only able to handle multiple (not multivariate) stationary time series (e.g., Cadonna et al., 2019), or are methodologically adequate and flexible (e.g., Li and Krafty, 2018; Bruce et al., 2018), but computationally unfeasible to jointly analyze more than a relatively small number of multivariate time series components.

In the univariate context Yang et al. (2016) considers a Bayesian lattice filter approach for analyzing a single time series which uses univariate dynamic linear models (DLMs) to describe the evolution of the forward and backward partial autocorrelation coefficients of such series. A key feature of this approach is its computational appeal. A DLM representation of a univariate time-varying autoregression requires a model with a state-parameter vector of dimension $P,$ where $P$ is the order of the autoregression (see e.g., West et al., 1999; Prado and West, 2010). Therefore, filtering and smoothing in this setting requires the inversion of $P\times P$ matrices at each time $t.$ Alternatively, the DLM formulation in the PARCOR domain requires fitting $2\times P$ DLMs, $P$ for the forward coefficients and $P$ for the backward coefficients, where each DLM has a univariate state-space parameter, fully avoiding matrix inversions and resulting in computational savings for cases in which $P\geq 2$.

In this paper we extend the Bayesian lattice filter approach of Yang et al. (2016) to the multivariate case. Our proposed models offer several advantages over currently available multivariate approaches for non-stationary time series including computational feasibility for joint analysis of relatively large-dimensional multivariate time series, and interpretable time-frequency analysis in the multivariate context. In particular, the proposed framework leads to posterior estimates of the time-varying spectral densities of each individual time series, as well as posterior measurements of the time-frequency relationships across multiple time series over time, such as time-varying coherence and partial coherence. We note that extending the approach Yang et al. (2016) to the multivariate case is non-trivial, as the closed-form inference used in the univariate DLM formulation of the lattice filter is not available for the multivariate case considered here. Multivariate DLM theory (West and Harrison, 1997; Prado and West, 2010) allows for full posterior inference in closed-form only when the covariance matrices of the innovations at the observation level and those at the system level are known, which is rarely the case in practice. Full posterior inference via MCMC can be obtained for more general multivariate DLM settings, but such posterior sampling schemes are very computationally expensive, making them only feasible when dealing with a small number of time series of small/moderate time lengths, and low-order TV-VAR models. We address these challenges by approximating the covariance matrices of the innovations at the observational level for the multivariate dynamic forward and backward PARCOR models using the approach of Triantafyllopoulos (2007). In addition, we use discount factors to specify the structure of the covariance matrices at the system levels. Our framework casts the time-varying multivariate representation of the input-output relations between the vectorial forward and backward predictions of a multivariate time series process –and their corresponding forward and backward matrices of PARCOR (partial autocorrelation) coefficients– as a Bayesian multivariate state-space model. Once approximate posterior inference is obtained for the multivariate time-varying PARCOR coefficients, posterior estimates for the implied time-varying vector autoregressive (TV-VAR) coefficient matrices and innovations covariance matrices can be obtained via Whittle’s algorithm (Zhou, 1992). Similarly, posterior estimates for any function of such matrices, such as the multivariate spectra and functions of the spectra, can also be obtained. A key feature of the proposed TV-VPARCOR representation is that it is more parsimonious and flexible than directly working with the TV-VAR state-space representation. We illustrate this in the analyses of simulated and real data presented in Sections 3 and 4. We also propose a method for selecting the number of stages in the TV-VPARCOR setting based on an approximate calculation of the Deviance Information Criterion (DIC).

The paper is organized as follows. Section 2 presents the models and discusses approximate posterior inference. Section 3 illustrates the performance of the proposed TV-VPARCOR models in simulation studies. Comparisons with results obtained from DLM representations of TV-VAR models are also provided. Section 4 presents analyses of two real multivariate temporal datasets. The first application considers joint analysis of multi-channel electroencephalogram data and the second one illustrates the analysis and forecasting of multi-location bivariate wind components. Finally, Section 5 presents a discussion and briefly describes potential future developments.

2. Models and Methods for Posterior Inference

2.1 Time-Varying Vector Autoregressive Models (TV-VAR) and Lattice Filters

Let \mbox{\boldmathx}_{t} be a $K\times 1$ vector time series for $t=1,\dots,T.$ A time-varying vector autoregressive model of order $P,$ referred to as TV-VAR$(P),$ is given by

[TABLE]

where \mbox{\boldmathA}_{t,j}^{(P)} is the $K\times K$ matrix of time-varying coefficients at lag $j,$ $j=1,\dots,P,$ and $\bm{\Sigma}_{t}$ is the $K\times K$ innovations variance-covariance matrix at time $t.$ The \mbox{\boldmath\epsilon}_{t}s are assumed to be independent over time.

Yang et al. (2016) considers a Bayesian lattice filter dynamic linear modeling approach for the case of univariate time-varying autoregressions (TVAR), i.e., when $K=1$ above. Such approach is based on a lattice structure formulation of the univariate Durbin-Levinson algorithm (see, e.g., Brockwell and Davis, 1991; Shumway and Stoffer, 2017) used in Kitagawa (2010). The idea is to obtain posterior estimation on the forward and backward time-varying PARCOR coefficients using a computationally efficient lattice filter representation. Once dynamic PARCOR estimation is obtained, estimates of the TVAR coefficients can be derived using the Durbin-Levinson recursion. The main advantage of using the dynamic PARCOR lattice filter representation instead of a dynamic linear model TVAR representation such as that used in West et al. (1999), is that the former avoids the inversion of $P\times P$ matrices required in the TVAR DLM filtering and smoothing equations. Instead, the PARCOR approach considers $2P$ dynamic linear models with univariate state parameters (e.g., $P$ DLMs with univariate state parameters for the forward coefficients and $P$ DLMs with univariate state parameters for the backward coefficients), completely avoiding matrix inversions. This is important for computational efficiency when considering models with $P>2$ and large $T$. The PARCOR approach also offers additional modeling advantages due to the fact that considering $2P$ DLMs with univariate state parameters generally provides more flexibility, than using a single DLM TVAR with $P$-dimensional state parameters.

We extend the approach of Yang et al. (2016) to consider multivariate non-stationary time series. More specifically, we consider Bayesian multivariate DLMs that use the multivariate Whittle algorithm (Zhou, 1992), also known as the multivariate Durbin-Levinson algorithm (Brockwell and Davis, 1991), to obtain a representation of the TV-VAR coefficient matrices in terms of time-varying PARCOR matrices as follows. Let \mbox{\boldmathf}_{t}^{(P)} and \mbox{\boldmathb}_{t}^{(P)} be the $K$-dimensional prediction error vectors at time $t$ for the forward and backward TV-VAR$(P)$ model, respectively, where,

[TABLE]

\mbox{\boldmathA}_{t,j}^{(P)} and \mbox{\boldmathD}_{t,j}^{(P)} denote, respectively, the $K\times K$ time-varying matrices of forward and backward TV-VAR$(P)$ coefficients for $j=1,\ldots,P.$ Similarly, \mbox{\boldmathA}_{t,j}^{(m)} and \mbox{\boldmathD}_{t,j}^{(m)} denote the time-varying matrices of forward and backward TV-VAR$(m)$ coefficients for $j=1,\ldots,m.$ Then, we write the $m$-stage of the lattice filter in terms of the pair of input-output relations between the forward and backward $K$-dimensional vector predictions, as follows,

[TABLE]

where $\bm{\Lambda}_{t,m}^{(m)}$ and $\bm{\Theta}_{t,m}^{(m)}$ are, respectively, the $K\times K$ matrices of time-varying forward and backward PARCOR coefficients for $m=1,\ldots,P.$ Note that for stationary AR$(P)$, i.e., models with $K=1$ and static AR coefficients in the stationary region, the forward and backward PARCOR coefficients are equal, i.e., $\lambda_{m}^{(m)}=\theta_{m}^{(m)}$ for all $m.$ For general $K$ and non-stationary processes the forward and backward PARCOR coefficients are not the same.

For each stage $m$ of the lattice structure above, we obtain the forward and backward TV-VAR coefficient matrices, \mbox{\boldmathA}_{t,m}^{(P)} and \mbox{\boldmathD}_{t,m}^{(P)}, from the time-varying forward and backward PARCOR coefficient matrices, $\bm{\Lambda}_{t,m}^{(m)}$ and $\bm{\Theta}_{t,m}^{(m)}$, using Whittle’s algorithm (see, e.g., Zhou, 1992), i.e.,

[TABLE]

with \mbox{\boldmathA}_{t,m}^{(m)}=\bm{\Lambda}_{t,m}^{(m)} and \mbox{\boldmathD}_{t,m}^{(m)}=\bm{\Theta}_{t,m}^{(m)}, for $m=1,\ldots,P.$

2.2 Model specification and inference

Our proposed model specification uses equations (1) and (2) as observational level equations of multivariate DLMs (West and Harrison, 1997; Prado and West, 2010) on the forward and backward PARCOR time-varying coefficients. These multivariate DLMs are specified as follows. For each $t$, let $\mbox{vec}(\bm{\Lambda}^{(m)}_{t,m})$ and $\mbox{vec}(\bm{\Theta}^{(m)}_{t,m})$ be the vectorized forward and backward PARCOR coefficients, i.e., these are $K^{2}$ vectors obtained by stacking the forward and backward PARCOR coefficient matrices at time $t$, $\bm{\Lambda}^{(m)}_{t,m}$ and $\bm{\Theta}^{(m)}_{t,m},$ by columns, respectively. In addition, define the forward and backward $K\times K^{2}$ matrices \mbox{\boldmathF}^{(m-1)}_{t+m}=(\mbox{\boldmathf}_{t+m}^{(m-1)})\otimes\mbox{\boldmathI}_{K\times K} and \mbox{\boldmathB}^{(m-1)}_{t-m}=(\mbox{\boldmathb}_{t-m}^{(m-1)})\otimes\mbox{\boldmathI}_{K\times K}, where \mbox{\boldmathI}_{K\times K} denotes the $K\times K$ identity matrix and $\otimes$ denotes the Kronecker product. Then, equations (1) and (2) can be rewritten as

[TABLE]

which correspond to the observational equations of two multivariate dynamic linear regressions on \mbox{\boldmathf}_{t}^{(m-1)} and \mbox{\boldmathb}_{t}^{(m-1)}, with dynamic coefficients $\mbox{vec}(\bm{\Lambda}_{t,m}^{(m)})$ and $\mbox{vec}(\bm{\Theta}_{t,m}^{(m)}),$ respectively. In order to complete the MDLM structure we specify random walk evolution equations for $\mbox{vec}(\bm{\Lambda}_{t,m}^{(m)})$ and $\mbox{vec}(\bm{\Theta}_{t,m}^{(m)})$ as follows,

[TABLE]

where \mbox{\boldmathW}_{t,f,m} and $\bm{W}_{t,b,m}$ are time dependent system covariance matrices. Finally, we specify prior distributions for $\mbox{vec}(\bm{\Lambda}_{0,m}^{(m)})$ and $\mbox{vec}({\bm{\Theta}}_{0,m}^{(m)})$ and all $m.$ We use conjugate normal priors for these parameters, i.e., we assume

[TABLE]

where $\mathcal{D}_{0,f,m}$ and $\mathcal{D}_{0,b,m}$ denote the information available at time $t=0$ for the forward and backward state parameter vectors, respectively.

Given $\bm{\Sigma}_{t,f,m},$ and \mbox{\boldmathW}_{t,f,m} for all $t=1,\ldots,T,$ and all $m=1,\ldots,P,$ equations (5), (7) and (9) define a normal MDLM (see, e.g., Prado and West, 2010, Chapter 10) for the forward time-varying PARCOR. Similarly, given $\bm{\Sigma}_{t,b,m}$ and \mbox{\boldmathW}_{t,b,m} for all $t$ and all $m,$ equations (6), (8) and (10) define a normal MDLM for the backward time-varying PARCOR.

Note that posterior inference in the case of univariate models with $K=1$ is available in closed form via the DLM filtering and smoothing equations. This is used in Yang et al. (2016) to obtain posterior inference in this univariate case. However, posterior inference in the general multivariate setting proposed here is not available in closed form when the observational and system covariance matrices are unknown, which is typically the case in practical settings. Therefore, as explained below, we use discount factors to specify \mbox{\boldmathW}_{t,f,m}, and \mbox{\boldmathW}_{t,b,m}. We also assume $\bm{\Sigma}_{t,f,m}=\bm{\Sigma}_{f,m}$ and $\bm{\Sigma}_{t,b,m}=\bm{\Sigma}_{b,m}$ for all $t,$ and use the approach of Triantafyllopoulos (2007) to obtain estimates of $\bm{\Sigma}_{f,m}$ and $\bm{\Sigma}_{b,m},$ which allows us to get approximate posterior inference in the multivariate case.

We first define the $K^{2}\times K^{2}$ system covariance matrices using discount factors by setting

[TABLE]

where each component, $\bm{\delta}_{.,m,i},$ is a $K$-dimensional vector that contains the discount factors for each of the $K$ components at stage $m$. Although we can assume different discount factors for different elements of $\bm{\delta}_{\cdot,m,k}$ and also across different $k$s, in practice we usually set all the elements of $\bm{\delta}_{f,m,k}$ equal to $\delta_{f,m}$ and all the elements of $\bm{\delta}_{b,m,k}$ equal to $\delta_{b,m}$ for all $k=1,\ldots,K,$ and then choose $\delta_{f,m}$ and $\delta_{b,m}$ optimally according to some criterion for each stage $m$ (this is discussed in Section 2.3). This structure for $\bm{\Delta}_{f,m}$ and $\bm{\Delta}_{b,m}$ allows us to obtain closed form expressions for \mbox{\boldmathW}_{t,f,m} and \mbox{\boldmathW}_{t,b,m} sequentially over time.

We now describe the full algorithm for approximate posterior inference in the forward TV-VPARCOR model. The algorithm for the backward model is similar. Let $\mathcal{D}_{t,f,m}$ denote all the information available up to time $t$ at stage $m$ for the forward model, with \mathcal{D}_{t,f,m}=\{\mathcal{D}_{t-1,f,m},\mbox{\boldmathf}_{t}^{(m-1)}\}. Consider the posterior expectation of $\bm{\Sigma}_{f,m}$ up to time $t,$ i.e., $E(\bm{\Sigma}_{f,m}|\mathcal{D}_{t,f,m}),$ and assume that $\lim_{t\rightarrow\infty}E(\bm{\Sigma}_{f,m}|\mathcal{D}_{t,f,m})=\bm{\Sigma}_{f,m}.$ Let $n_{0,f,m}$ be a positive scalar and \mbox{\boldmathS}_{0,f,m} be the prior expectation of $\bm{\Sigma}_{f,m}.$ Assume that at time $t-1,$ we have that $\mbox{vec}(\bm{\Lambda}_{t-1,m}^{(m)})|\mathcal{D}_{t-1,f,m}$ is approximately distributed as N(\mbox{\boldmathm}_{t-1,f,m},\mbox{\boldmathC}_{t-1,f,m}), and so, E(\mbox{\boldmathf}_{t}^{(m-1)}|\mathcal{D}_{t-1,f,m}) is approximated by \mbox{\boldmathB}_{t-m}^{(m-1)}\mbox{\boldmathm}_{t-1,f,m} and V(\mbox{\boldmathf}_{t}^{(m-1)}|\mathcal{D}_{t-1,f,m}) is approximated by \bm{Q}_{t-1,f,m}=\mbox{\boldmathB}_{t-m}^{(m-1)}\mbox{\boldmathR}_{t,f,m}(\mbox{\boldmathB}_{t-1}^{(m-1)})^{\prime}+\mbox{\boldmathS}_{t-1,f,m}, with \mbox{\boldmathR}_{t,f,m}=\mbox{\boldmathC}_{t-1,f,m}+\mbox{\boldmathW}_{t,f,m}, for some \mbox{\boldmathS}_{t-1,f,m}. Then, following Theorem 1 of Triantafyllopoulos (2007) we have that, if $\bm{\Sigma}_{f,m}$ is bounded, \mbox{\boldmathS}_{t,f,m} will approximate $\bm{\Sigma}_{f,m}$ for $t$ large, with

[TABLE]

where in our case \mbox{\boldmathe}_{t,f,m}=\mbox{\boldmathf}_{t}^{(m-1)}-\mbox{\boldmathB}_{t-m}^{(m-1)}\mbox{\boldmathm}_{t-1,f,m}, and \mbox{\boldmathS}_{i-1,f,m}^{1/2},\bm{Q}_{i,f,m}^{-1/2} are symmetric square roots of the matrices \mbox{\boldmathS}_{i-1,f,m} and $\bm{Q}_{i,f,m}^{-1},$ respectively, based on the spectral decomposition factorization of symmetric positive definite matrices for all $i=1,\ldots,t.$

Using the approximation above we obtain the filtering equations below for approximate inference in the forward TV-VPARCOR model.

The one-step ahead forecast mean and covariance at time $t$ are given by:

[TABLE]

and

[TABLE]

where \mbox{\boldmathR}_{t,f,m}=\mbox{\boldmathC}_{t-1,f,m}+\mbox{\boldmathW}_{t,f,m} and \mbox{\boldmathW}_{t,f,m}=\Delta_{f,m}\mbox{\boldmathC}_{t-1,f,m}\Delta_{f,m}-\mbox{\boldmathC}_{t-1,f,m}.

• -

The one-step forecast error vector is given by \mbox{\boldmathe}_{t,f,m}=\mbox{\boldmathf}^{(m-1)}_{t}-\mbox{\boldmathB}_{t-m}^{(m-1)}\mbox{\boldmathm}_{t-1,f,m}.

• -

Using Bayes’ theorem and the equations above we can obtain the approximate posterior distribution at time $t$ as \mbox{vec}(\bm{\Lambda}^{(m)}_{t,m})|\mathcal{D}_{t,f,m}\approx\mathcal{N}(\mbox{\boldmathm}_{t,f,m},\mbox{\boldmathC}_{t,f,m}), where

[TABLE]

Approximate filtering and predictive distributions for $\mbox{vec}(\bm{\Lambda}_{t,m}^{m})|\mathcal{D}_{t,f,m}$, $\bm{f}^{(m-1)}_{t}|\mathcal{D}_{t,f,m}$ and $\bm{f}^{(m-1)}_{t+h}|\mathcal{D}_{t,f,m}$ for a positive integer $h>0$ can also be obtained by taking $\bm{\Sigma}_{f,m}=\bm{S}_{t,f,m}.$

After applying the filtering equations up to time $T$, it is possible to compute approximate smoothing distributions for the forward PARCOR model by setting $\bm{\Sigma}_{f,m}=\bm{S}_{T,f,m}.$ This leads to approximate smoothing distributions

[TABLE]

where the mean and covariance are computed recursively via

[TABLE]

for $t=(T-1),\ldots,1,$ with $\bm{J}_{t,f,m}=\bm{C}_{t,f,m}\bm{R}_{t+1,f,m}^{-1},$ and starting values $\bm{a}_{T}(0)=\bm{m}_{T,f,m}$ and \mbox{\boldmathR}_{T,f,m}(0)=\bm{C}_{T,f,m}. Filtering and smoothing equations can be obtained for the backward PARCOR model in a similar manner. Finally, the algorithm for approximate posterior estimation is as follows.

Algorithm

Given hyperparameters $\{P,\bm{\Delta}_{f,m},\bm{\Delta}_{b,m};m=1,\dots,P\},$ set $\bm{f}_{t}^{(0)}=\bm{b}_{t}^{(0)}=\bm{x}_{t},\;\textrm{for }t=1,\dots,T.$

• 2.

Use $\{\bm{f}^{(0)}_{t}\}$ and $\{\bm{b}^{(0)}_{t}\}$ as vectors of responses in the observational level equations (1) and (2), respectively, which, combined with the random walk evolution equations (7) and (8), and the priors (9) and (10), define the multivariate PARCOR forward and backward models. Then, use the sequential filtering equations (12) to (14) to obtain the estimated $\{\bm{S}_{T,f,1}\}$ and $\{\bm{S}_{T,b,1}\}$. Use the sequential filtering equations (12) to (14) along with the smoothing equations (15) and (16) to obtain a series of estimated parameters $\{\mbox{vec}(\hat{\bm{\Lambda}}^{(1)}_{t,1})\}$, $\{\mbox{vec}(\hat{\bm{\Theta}}^{(1)}_{t,1})\}$ for $t=1:T$. These estimated parameters are set at the posterior means of the smoothing distributions, i.e., the values in (16) for the forward case and a similar equation in the backward case.

• 3.

Use the observational equations (1) and (2) to obtain the new series of forward and backward prediction errors, $\{\bm{f}_{t}^{(1)}\}$ and $\{\bm{b}_{t}^{(1)}\},$ for $t=1,\dots,T.$

• 4.

Repeat steps 2-3 above until $\{\mbox{vec}(\hat{\bm{\Lambda}}^{(m)}_{t,m})\}$, $\{\mbox{vec}(\hat{\bm{\Theta}}^{(m)}_{t,m})\}$, $\{\bm{S}_{T,f,m}\}$ and $\{\bm{S}_{T,b,m}\}$ have been obtained for all $m=1,\ldots,P.$

• 5.

Finally, use $\{\mbox{vec}(\hat{\bm{\Lambda}}^{(m)}_{t,m})\}$ and $\{\mbox{vec}(\hat{\bm{\Theta}}^{(m)}_{t,m})\}$, for $m=1,\dots,P,$ and equations (3) and (4) to obtain the forward and backward TV-VAR coefficient matrices via Whittle’s algorithm.

2.3 Model selection and time-frequency representation

In order to select the optimal model order and discount factors, we begin by specifying a potential maximum value of $P,$ say $P_{\max},$ for the model order. At level $m$ we search for the optimal values of $\bm{\Delta}_{f,m}$ and $\bm{\Delta}_{b,m}.$ In other words, at level $m=1$ we search for the combination of values of $\bm{\Delta}_{f,1}$ and $\bm{\Delta}_{b,1}$ maximizing the log-likelihood resulting from $(\ref{eqn1})$ with $m=1$. Using the selected optimal $\bm{\Delta}_{f,1}$ and $\bm{\Delta}_{b,1}$, we can obtain the corresponding series $\{\bm{f}_{t}^{(2)}\}$ and $\{\bm{b}_{t}^{(2)}\}$, for $t=1,\dots,T$, as well as the maximum log-likelihood value $\mathcal{L}_{f,1}.$ Then, we repeat the above search procedure for stage two, i.e., $m=2$, using the output $\{\bm{f}_{t}^{(2)}\}$ and $\{\bm{b}_{t}^{(2)}\}$ obtained from implementing the filtering and smoothing equations with the previously selected hyperparameters $\bm{\Delta}_{f,1}$ and $\bm{\Delta}_{b,1}$. We obtain optimal $\bm{\Delta}_{f,2}$, $\bm{\Delta}_{b,2}$ as well as $\{\bm{f}_{t}^{(3)}\}$ and $\{\bm{b}_{t}^{(3)}\}$, for $t=1,\dots,T.$ We also obtain the value of the corresponding maximum log-likelihood $\mathcal{L}_{f,2}$. We repeat the procedure until the set $\{\bm{\Delta}_{f,m},\bm{\Delta}_{b,m},\mathcal{L}_{f,m}\},m=1,\dots P_{\max},$ has been selected. We then consider two different methods for selecting the optimal model order as described below. Note that one can also obtain the optimal likelihood values from the backward model, $\mathcal{L}_{b,m},$ for $m=1,\ldots P_{\max}.$ For all the examples and real data analyses presented below we choose the optimal model orders based on the optimal likelihood values for the forward model. Similar results were obtained based on the optimal likelihood values for the backward models.

Method 1: Scree plots. This method was used by Yang et al. (2016) to select the model order visually by plotting $\mathcal{L}_{f,m}$ against the order $m.$ The idea is that, when the observed vector of time series truly follows a TV-VAR model, the values of $\mathcal{L}_{f,m}$ will stop increasing after a specific lag and this lag is then chosen to be the model order. A numerical version of this method can also be implemented by computing the percent of change in the likelihood going from $\mathcal{L}_{f,m-1}$ to $\mathcal{L}_{f,m},$ however, here we use scree plots as a visualization tool and use the model selection criterion below to numerically find an optimal model order.

Method 2: DIC model selection criterion. We consider an approach based on the deviance information criterion (DIC) to choose the model order (see Gelman et al., 2014, and references therein). In general, for a model with parameters denoted as $\bm{\theta},$ the DIC is defined as

[TABLE]

where $\bm{y}$ denotes the data, $\hat{\bm{\theta}}_{Bayes}$ is the Bayes estimator of $\bm{\theta}$ and $p_{DIC}$ is the effective number of parameters. The effective number of parameters is given by

[TABLE]

where the expectation in the second term is an average of $\bm{\theta}$ over its posterior distribution. The expression above is typically estimated using samples $\bm{\theta}^{s},s=1,\dots,S,$ from the posterior distribution as

[TABLE]

Note, however, that in our case we do not have samples from the exact posterior distribution of the parameters since we are using approximate inference to avoid computationally costly exact inference via MCMC. Therefore, for a given model order $m$ we compute the likelihood term in the DIC calculation approximately using the forward filtering distributions as explained below. Also, note that, fitting a PARCOR model at stage $m$ requires fitting all the models of the previous $m-1$ stages. Therefore, the effective number of parameters at stage $m$ is computed by adding the estimated effective number of parameters of stage $m$ plus the estimated effective number of parameters for the previous $m-1$ stages. In other words, for each stage $m:$

Compute the estimated implied log-likelihood from equation (5) for $t=1,\ldots,T,$ using ${\mbox{vec}(\hat{\bm{\Lambda}}_{t,m}^{(m)}})$ and $\bm{S}_{T,m,f}$. In this way we obtain the first term in the calculation of the DIC for model order $m.$

• -

Obtain samples, $\mbox{vec}(\bm{\Lambda}^{(m)}_{t,m,s})$, for $s=1,\ldots,S,$ from the approximate sequential filtering equations with distributions $\mathcal{N}(\bm{m}_{t,f,m},\bm{C}_{t,f,m}),$ and use these samples to compute the estimated number of parameters related only to stage $m$ which we denote as $\hat{p}_{DIC,m}^{m}$. Note that, as mentioned above, stage $m$ requires fitting all the PARCOR models for the previous $(m-1)$ stages and so, in the final DIC calculation at stage $m$ the total estimated effective number of parameters is computed as

[TABLE]

We denote the final estimated DIC for model order $m$ as $\widehat{DIC}_{m}.$

2.4 Posterior summaries

Once an optimal TV-VPARCOR model is chosen we can obtain posterior summaries of any quantities associated to such model. For instance, we can obtain posterior summaries of the TV-VPARCOR coefficients over time at each stage, and consequently summaries of the corresponding TV-VAR$(P)$ coefficients over time.

Time-frequency representations are generally more useful in practice, and these can be obtained by computing the spectral density matrix, $\bm{g}(t,\omega),$ for any time $t$ and frequency $\omega\in(0,1/2),$ as well as measurements derived from this matrix such as coherence or partial coherence. The spectral density matrix is estimated as

[TABLE]

where $\bm{\hat{\Phi}}(t,\omega)=\bm{I}-\sum_{m=1}^{P}\bm{\hat{A}}^{(P)}_{t,m}\exp\{-2\pi im\omega\},$ with $i=\sqrt{-1}$ (see e.g., Shumway and Stoffer, 2017, Chapter 4). $\bm{\hat{\Sigma}}$ can be set at $\mathbf{S}_{T,f,P}.$ Note that the spectral density matrix ${\bm{g}}(t,\omega)$ consists of individual spectra ${g}_{j,j}(t,\omega)$ for each component $j=1,\dots,K$ of \mbox{\boldmathx}_{t}, and the cross-spectra ${g}_{i,j}(t,\omega)$ between components $i$ and $j.$ From these we can compute the estimated squared coherence between components $i$ and $j$ as

[TABLE]

for all $i\neq j.$ This measure is used to estimate the power transfer between two components of the time series. Similarly, the partial squared coherence between components $i$ and $j$ can be estimated as follows. Let $\bm{c}(t,\omega)={\bm{g}}^{-1}(t,\omega)$ be the inverse of the spectral density matrix with elements $c_{i,j}(t,\omega)$ for $i,j=1,\dots,K.$ Then, the estimated squared partial coherence between components $i$ and $j$ is given by

[TABLE]

The squared partial coherence is essentially the frequency domain squared correlation coefficient between components $i$ and $j$ after the removal of the linear effects of all the remaining components of \mbox{\boldmathx}_{t}. Finally, uncertainty measures for the spectral density matrix, and any functions of this matrix, can be obtained from the approximate filtering and smoothing posterior distributions of the forward and backward TV-VPARCOR models.

2.5 Forecasting

In this section, we show how to obtain $h$-steps ahead forecasts. In order to have a non-explosive behavior in the forecasts, we assume the series is locally stationary in the future, i.e, $\bm{\Lambda}^{(m)}_{t,m}=\bm{\Theta}^{(m)}_{t,m}$ at time $t=T+1,\dots,T+h$. Then, the approximate $h$-steps ahead forecast posterior distribution of the PARCOR coefficients, with $h>0,$ is approximated as $(\bm{\Lambda}_{T+h,m}^{(m)}|\mathcal{D}_{T,f,m})\approx\mathcal{N}(\bm{m}_{T,f,m}(h),\bm{C}_{T,f,m}(h))$, where

[TABLE]

with \bm{W}_{T+1,f,m}=\Delta_{f,m}\mbox{\boldmathC}_{T,f,m}\Delta_{f,m}-\mbox{\boldmathC}_{T,f,m}, for $m=1,\ldots,P.$ Then, we apply Whittle’s algorithm to transform the PARCOR coefficients, $\bm{\Lambda}_{T+h,P}^{(P)}$, into TV-VAR coefficients $\bm{A}_{T+h,j}^{(P)}$ and $\bm{D}_{T+h,j}^{(P)}$, for $j=1,\dots,P.$ Finally, we obtain the $h$-steps ahead forecasts using

[TABLE]

3. Simulation Studies

In this section we illustrate our proposed approach in the analysis of simulated data. The relative performances of the models considered here, including that of the proposed TV-VPARCOR, were assessed by computing the average squared error (ASE) between the estimated spectral density matrix and the true spectral density matrix.

3.1 Bivariate TV-VAR$(2)$ processes

We simulated $50$ bivariate time series of length $T=1024$ from the following TV-VAR$(2)$ model:

[TABLE]

with

[TABLE]

where $r_{1,t}=\frac{0.1}{1024}t+0.85,$ $r_{2,t}=-\frac{0.1}{1024}t+0.95,$ $\lambda_{1,t}=\frac{15}{1024}t+5,$ and $\lambda_{2,t}=-\frac{10}{1024}t+15.$ We also considered two values for $\phi_{1,1,2},$ namely (i) $\phi_{1,1,2}=0$; (ii) $\phi_{1,1,2}=-0.8.$ The true $2\times 2$ spectral matrix of this process is given by

[TABLE]

where $\bm{\Phi}(t,\omega)=\bm{I}_{2}-\bm{\Phi}_{1,t}\exp\{-2\pi i\omega\}-\bm{\Phi}_{2,t}\exp\{-4\pi i\omega\}$, and $\bm{\Sigma}=\bm{I}_{2}.$ The spectral matrix ${\bm{g}}(t,\omega)$ is symmetric, with corresponding components ${g}_{11}(t,\omega),$ ${g}_{12}(t,\omega)$ and ${g}_{22}(t,\omega)$, representing, respectively, the spectrum of the first component, the co-spectrum between the first and the second components, and the spectrum of the second component. The squared coherence between the first and second components is given by

[TABLE]

Note that when $\phi_{1,1,2}=0$, the two processes are uncorrelated and ${g}_{1,2}(t,\omega)=\rho_{12}^{2}(t,\omega)=0$ for all $t$ and $\omega.$ Figure 1 shows the true log spectral densities ${g}_{11}(t,\omega)$ and ${g}_{22}(t,\omega).$ The true log spectral densities and square coherences for the case of $\phi_{1,1,2}=-0.8$ are shown in the top row plots of Figure 4.

We fit bivariate TV-VPARCOR models to each of the 50 simulated bivariate time series under cases (i) and (ii). We set a maximum of order $P_{\max}=5$. The elements of the diagonal component of discount factor matrices $\bm{\Delta}_{f,m}$ and $\bm{\Delta}_{b,m},$ $\delta_{f,m}$ and $\delta_{b,m}$ respectively, were chosen from a grid of values in $(0.995,1)$. We set the hyperparameters $n_{f,m,0}=n_{b,m,0}=1$, $\bm{S}_{0,f,m}=\bm{S}_{0,b,m}=\bm{I}_{2}$, $\bm{m}_{0,f,m}=\bm{m}_{0,b,m}=(0,0,0,0)^{\prime}$ and $\bm{C}_{0,f,m}=\bm{C}_{0,b,m}=\bm{I}_{4}$. For comparison, we also fit TV-VAR models to the simulated bivariate data with model orders ranging from $1$ to $5$. Multivariate DLM representations of bivariate TV-VAR$(m)$ process were considered for each $m=1,\ldots,5$. Each TV-VAR representation has an $(m\times 4)$-dimensional state parameter vector. For each model order a single optimal discount factor, $\delta_{m}$ was chosen from a grid of values in $(0.995,1).$ Furthermore, in order to provide a similar model setting to the one we used in our TV-VPARCOR approach, the covariance matrix at the observational level in the DLM formulation for each TV-VAR$(m)$ was also specified following the approach of Triantafyllopoulos (2007).

The top plots in Figure 2 show the BLF-scree plots obtained from the PARCOR approach for each of the 50 datasets under the two scenarios (i) $\phi_{1,1,2}=0$ and (ii) $\phi_{1,1,2}=-0.8$ for model orders $m=1,\ldots,5.$ We see that in both scenarios the BLF-scree plots indicate that the optimal model order is $P=2.$ We also computed the DIC as explained in the previous section for each model order $m=1,\ldots,5$ and each dataset under scenarios (i) and (ii). The bottom left plot in Figure 2 shows the distributions of the optimal model orders chosen by the TV-VPARCOR and TV-VAR approaches for scenario (i), while the right bottom plot shows the distribution for scenario (ii). We see that both, the TV-VPARCOR and TV-VAR approaches lead to very similar results and model order 2 is adequately chosen as the optimal model order under the two scenarios for most of the 50 datasets.

Figures 3 and 4 summarize posterior inference obtained from the TV-VPARCOR approach using a model order of 2 for the two scenarios (i) and (ii), respectively. Estimated spectral densities were obtained from the posterior means of the approximate smoothing distributions of the forward and backward PARCOR coefficient matrices over time. The estimated log spectral densities displayed in the figures were obtained by averaging the estimated log spectral densities over the 50 simulated datasets. The bivariate TV-VPARCOR model is able to adequately capture the structure of the individual spectral densities and also that of the squared coherences. From these figures we also see that when $\phi_{1,1,2}=-0.8$, the second series has stronger impact on the first one and therefore their coherence is stronger. The TV-VPARCOR model is able to adapt and adequately capture this feature.

In order to compare the performance of the TV-VPARCOR and TV-VAR models in estimating the various time-frequency representations, we computed the mean and standard deviations of the average squared error (ASE) for each of the models in each of the two simulation scenarios. The ASE is defined as follows (Ombao et al., 2001)

[TABLE]

where $\omega_{l}=0,0.001,0.011,\dots,0.5$. Note that we have $n=50$ simulated datasets for each of the two scenarios. Table 1 summarizes the mean and standard deviations of the ASE based on $\mbox{ASE}_{n}$ for scenarios (i) and (ii). Note that the simulated data are actually generated from TV-VAR$(2)$ models, not from TV-VPARCOR models, so we expect TV-VAR models to do better in terms of ASE for this specific simulation study. Nevertheless the proposed TV-VPARCOR approach has comparable performance in terms of estimating the time-frequency characteristics of the original process while being computationally more efficient.

In fact, Table 2 presents the computation times for both models again averaging over the 50 realizations in each case. We see that even for this example with only two time series components and a model order of 2, the TV-VPARCOR models require almost half of the computation time required by the TV-VAR$(2)$ models. As the model order and the number of time series components increase, differences in computational will be more pronounced, making the TV-VPARCOR approach more efficient for modeling large temporal datasets.

3.2 20-dimensional TV-VAR(1)

We analyze data simulated from a 20-dimensional non-stationary TV-VAR$(1)$ process with $T=300$ in which the $(i,j)$ elements of the matrix of VAR coefficients at time $t,$ $\bm{\Phi}_{t},$ are given as follows:

[TABLE]

for $t=1,\ldots,300.$ In addition, we assume $\bm{\Sigma}=0.1\bm{I}_{20}.$

We fit TV-VPARCOR models considering $P_{max}=3.$ Note that the PARCOR approach with $P_{max}=3$ requires fitting 6 multivariate DLMs with state-space parameter vectors of dimension $400.$ Alternatively, working directly with TV-VAR representations with $P_{max}=3$ requires fitting 3 multivariate DLMs with state-space parameter vectors of dimension $400$ for model order 1, $800$ for model order 2, and $1200$ for model order 3. The TV-VAR model representation leads to a rapid increase of the dimension of the state-space vector with the model order, which significantly reduces the computational efficiency, particularly for large and even moderate $T.$ The TV-VPARCOR approach requires fitting more multivariate DLMs, but the dimensionality of the state-space vectors remains constant with the model order. This is an important advantage of the TV-VPARCOR approach. In fact, the TV-VPARCOR model required 585s of computation time for $P_{\max}=3,$ while the TV-VAR model required 3379s with the same $P_{\max}=3$ value. Posterior computations were completed in both cases using a MacBookPro13 with Intel Core i5, with a 2 GHz (1 Processor). Note also that, for a given model order the PARCOR approach can be further optimized in terms of computational efficiency, as the forward and backward DLMs can run in parallel.

We assumed prior hyperparameters $\bm{m}_{0,\cdot,m}=\bm{0}$ and $\bm{C}_{0,\cdot,m}=\bm{I}_{400}$ for the forward and backward PARCOR models. The elements of the diagonal component of discount factor matrices, $\delta_{f,m}$ and $\delta_{b,m},$ were chosen from a grid of values in $(0.99,1)$. As mentioned above we also fit TV-VAR models with model orders going from $1$ to $3$ using similar prior hyperparameters and discount factors. For both types of models the DIC picked model order $1$ as the optimal model order, which is the corresponding true model order in this case. Both types of models led to similar posterior inference of the time-frequency spectra.

Here we only show the results from the TV-VPARCOR approach. Figure 5 shows the true and estimated log spectral densities from the TV-VPARCOR model for 4 components of the 20-dimensional time series, namely, components $1$, $2$, $8$ and $15$. Figure 6 shows the true and estimated coherences between components $1$ and $5$, components $2$ and $15$, components $5$ and $12$, and components $15$ and $20$. Overall we see that the TV-VPARCOR approach adequately captures the space-time characteristics of the original multivariate non-stationary time series process. Furthermore, the TV-VPARCOR approach led to similar posterior estimates of the VAR coefficients over time to those obtained from using a DLM representation of a TV-VAR (see Figure 1 in the Supplementary Material).

4. Case studies

4.1 Analysis of multi-channel EEG data

We consider the analysis of multi-channel EEG data recorded on a patient that received electroconvulsive therapy (ECT) as a treatment for major depression. These data are part of a larger dataset, code named Ictal19, that corresponds to recordings of $19$ EEG channels from one subject during ECT. As an illustration, we use our multivariate TV-VPARCOR model to analyze $9$ of the channels, specifically channels $F_{3},F_{z},F_{4},C_{3},C_{z}$, $C_{4},P_{3},P_{z},P_{4}$ shown in Figure 7. We chose these channels because the are closely located and because based on previous analyses we expect strong similarities in their temporal structure over time. The full multi-channel dataset was analyzed in West et al. (1999) and Prado et al. (2001) using univariate TVARs separately for each channel, and also using dynamic regression models. The original recordings of about $26,000$ observations per channel were subsampled every sixth observation from the highest amplitude portion of the seizure, leading to a set of series of $3,600$ observations (corresponding to $83.72s$) per channel (Prado et al., 2001).

We analyzed the $K=9$ series listed above jointly using a multivariate TV-VPARCOR model. We considered a maximum model order of $P_{\max}=20$ and discount factor values on a grid in the $(0.99,1]$ range (with equal spacing of $0.001$). We further assumed that the discount factor values were the same across channels. This assumption was based on previous analyses of the individual channels using univariate TVAR models that showed similar optimal discount values for the different channels. We set $n_{0,f,m}=n_{0,b,m}=1,$ and $S_{0,f,m}=S_{0,b,m}=2000\bm{I}_{9}$ for all $m.$ In addition, we set the same initial prior parameters $\bm{m}_{0,f,m}=\bm{m}_{0,b,m}=\bm{0}$ and $\bm{C}_{0,f,m}=\bm{C}_{0,b,m}=1000\bm{I}_{81}.$ The optimal model order was found to be 5 (see Figure 2 in the Supplementary Material) and so, the results presented here correspond to a TV-VPARCOR model with this order. Higher order models were also fitted leading to similar but slightly smoother results in terms of the estimated spectral density, coherence and partial coherence.

Figure 8 displays estimated log spectral densities of channels Cz, Pz and F4. We note that the multi-channel EEG data are dominated by frequency components in the lower frequency band (below 18Hz). Furthermore, each EEG channel shows a decrease in of the dominant frequency over time, starting around 5Hz and ending around approximately 3Hz. This decrease in the dominant frequency was also found in West et al. (1999). Channels Cz and Pz are more similar to each other than to channel F4 in terms of their log-spectral densities. The three channels show the largest power around the same frequencies, however, channel F4 displays smaller values in the power log-spectra than those for channels Cz and Pz. The remaining channels also show similarities in their spectral content (not shown).

Figure 9 shows estimated squared coherences (top) and estimated squared partial coherences (bottom) between channels Pz and Cz, F4 and Cz, and F4 and Pz. Channels Pz and Cz show a very strong coherence over time across almost all the frequency bands under 35Hz. On the other hand, channel F4 shows strong coherence with channels Pz and Cz across frequencies below 15-18 Hz at the beginning of the seizure. After the initial 10-15s, and approximately until about 50s, there is a strong coherence between F4 and Pz and Cz only at the dominant frequency of 3-5Hz that dissipates towards the end of the seizure. The partial coherence across pairs of channels is the frequency domain version of the squared correlation coefficient between relationship between pairs of components after the removal of the effects of all the other components. Figure 9 shows that the estimated squared partial coherences between Pz and Cz, F4 and Cz and F4 and Pz are essentially negligible for most frequency bands over the seizure course. This makes sense due to the fact that most of the 9 EEG channels are so strongly coherent across different frequency bands over the entire period of recording. The estimated squared partial coherence between channels Pz and Cz is large for frequencies below 5Hz only at the very beginning of the seizure. These findings are consistent with results from the analysis of these data in West et al. (1999) and Prado et al. (2001).

4.2 Analysis of multi-location wind data

We analyze wind component data derived from median wind speed and direction measurements taken every 4 hours from June 1st 2010 to August 15th in 3 stations in Northern California. These data were obtained from the Iowa Environmental Mesonet (IEM) Automated Surface Observing System (ASOS) Network, a publicly available database (see http://mesonet.agron.iastate.edu/ASOS/). ASOS stations are located at airports and take observations and basic reports from the National Weather Service (NWS), the Federal Aviation Administration (FAA), and the Department of Defense (DOD). For additional information about the ASOS measurements see NOA (1998). Here we analyze time series data from Monterey, Salinas and Watsonville, 3 stations located near the Monterey Bay.

We use the TV-VPARCOR approach for joint analysis of the six-dimensional time series corresponding to the wind time series components for the 3 stations. We set $P_{max}=10$ and consider discount factor values on a grid in the $(0.9,1]$ range. We assume that discount factor values were the same across components for the 3 stations. We set the prior hyperparameters as follows: $n_{0,f,m}=n_{0,b,m}=1,$ and $S_{0,f,m}=S_{0,b,m}=5\bm{I}_{6},$ $\bm{m}_{0,f,m}=\bm{m}_{0,b,m}=\bm{0}$ and $\bm{C}_{0,f,m}=\bm{C}_{0,b,m}=10\bm{I}_{36}$ for all $m.$ The optimal model order chosen by the approximate DIC calculation is $P=3$ (see Figure 3 in the Supplementary Material). For this model order we found that the optimal discount factors were $0.97,0.97$ and $0.99,$ respectively, for each of the 3 levels of the forward PARCOR model, and $0.98,$ $0.98$ and $0.99$ for each of the 3 levels of the backward PARCOR model.

Figure 10 shows the estimated log spectral densities of the East-West component (X component) and the North-South component (Y component) for each location. We can observe that there is a dominant quasi-periodic behavior around the 24 hour period for the East-West (X) components in Monterey and Salinas, as well as the North-South (Y) component in Watsonville. This quasi-periodic behavior is also present, although is less persistent over time, in the East-West component in Watsonville and the North-South components in Monterey and Salinas. The observed quasi-periodic pattern observed in the estimated log-spectral for these three locations is consistent with the fact that stronger winds are usually observed in the afternoons/evenings during the summer in these locations, while calmer winds are observed during the rest of the day. Note also that the quasi-periodic daily behavior is more persistent over the entire set of summer months for the North-South component than the East-West component in Watsonville, while the quasi-periodic behavior is more persistent in the East-West component than in the North-South component in Monterey and Salinas.

Figure 11 shows the estimated squared coherences between each pair of wind components across the three locations. There is a very strong coherence between Monterey and Salinas in the East-West (X) components for periods above 15 hours, with the strongest relationship observed around 24 hours. We also observe that in general, there is a strong coherence between all the components around the 24 hours period. This coherence relationship tends to be more marked across some locations during the month of June (e.g., between the North-South components of Monterey and Salinas).

Furthermore, the estimated squared partial coherence (see Figure 4 in the Supplementary Material) between the East-West components of Monterey and Salinas also shows that there is a relatively large linear relationship between these components for periods above 13 hours even after removing of the effect of all the other components for these locations and also after removing the effect of the wind components in Watsonville.

The TV-VPARCOR model can also be used for forecasting as described in Section 2.5. Figure 12 shows 72 hours forecasts obtained from the TV-VPARCOR model for the North-South wind component in Monterey. We see that the model adequately captures the general future behavior of this time series component.

5. Discussion

We consider a computationally efficient approach for analysis and forecasting of non-stationary multivariate time series. We propose a multivariate dynamic linear modeling framework to describe the evolution of the PARCOR coefficients of a multivariate time series process over time. We use approximations in this multivariate TV-VPARCOR setting to obtain computationally efficient and stable inference and forecasting in the time and time-frequency domains. The approximate posterior distributions derived from our approach are all of standard form. We also provide a method to choose the optimal number of stages in the TV-VPARCOR model based on an approximate DIC calculation. In addition, our model can provide reliable short term forecasting.

The proposed framework provides computational efficiency and excellent performance in terms of the average squared error between the true and estimated time-varying spectral densities as shown in simulation studies and in the analysis of two multivariate time series datasets. The TV-VPARCOR model representations also lead to very significant reduction in computational time when compared to TV-VAR model representations, particularly for cases in which we have model orders larger than 2-3 and more than a handful of time series components.

In addition to simulation studies we have shown that the TV-VPARCOR approach can be successfully used to analyze real multivariate non-stationary time series data. We presented the analysis of non-stationary multi-channel EEG data and also the analysis and forecasting of multi-location wind data. In the EEG case, our model was able to adequately detect the main time-frequency characteristics of individual EEG channels as well as the relationships across multiple channels over time. For the multi-location wind component data, our model detected a quasi-periodic pattern through the estimated spectral densities of each time series component which is consistent with the expected behavior of these components during the summer for locations near the Monterey Bay area. The model was also able to describe the time-varying relationships across multiple components and locations and led to reasonable short term forecasting.

The proposed dynamic multivariate PARCOR approach is computationally efficient when compared to state-space representations TV-VAR models. However, in many practical settings we may expect sparsity in the model parameters or situations in which some parameters change over time and others do not. Future work will explore inducing, possibly time-varying, sparsity and dimension reduction in these multivariate TV-VPARCOR models while maintaining computational efficiency and accuracy in inference and forecasting.

Supplementary Material

Supplementary figures for the 20-dimensional TV-VAR(1) example

Supplementary figures for the multi-channel EEG analysis

Supplementary figures for the multi-location wind component data