Bayesian nonparametric graphical models for time-varying parameters VAR

TL;DR

This paper introduces a Bayesian nonparametric framework for high-dimensional, time-varying VAR models that can infer dynamic Granger causality networks and handle non-linearities in macroeconomic data.

Contribution

It develops a novel BNP time-varying graphical model with a dependent Dirichlet Process prior, addressing overfitting and enabling flexible clustering of coefficients.

Findings

Effective inference of time-varying Granger causality networks.

Robustness demonstrated through comparison of prior specifications.

Model performs well on macroeconomic dataset.

Abstract

Over the last decade, big data have poured into econometrics, demanding new statistical methods for analysing high-dimensional data and complex non-linear relationships. A common approach for addressing dimensionality issues relies on the use of static graphical structures for extracting the most significant dependence interrelationships between the variables of interest. Recently, Bayesian nonparametric techniques have become popular for modelling complex phenomena in a flexible and efficient manner, but only few attempts have been made in econometrics. In this paper, we provide an innovative Bayesian nonparametric (BNP) time-varying graphical framework for making inference in high-dimensional time series. We include a Bayesian nonparametric dependent prior specification on the matrix of coefficients and the covariance matrix by mean of a Time-Series DPP as in Nieto-Barajas et al. (2012). Following Billio et al. (2019), our hierarchical prior overcomes over-parametrization and over-fitting issues by clustering the vector autoregressive (VAR) coefficients into groups and by shrinking the coefficients of each group toward a common location. Our BNP timevarying VAR model is based on a spike-and-slab construction coupled with dependent Dirichlet Process prior (DPP) and allows to: (i) infer time-varying Granger causality networks from time series; (ii) flexibly model and cluster non-zero time-varying coefficients; (iii) accommodate for potential non-linearities. In order to assess the performance of the model, we study the merits of our approach by considering a well-known macroeconomic dataset. Moreover, we check the robustness of the method by comparing two alternative specifications, with Dirac and diffuse spike prior distributions.

Equations

$${\mathbf{y}}_t=\sum _{i=1}^p{\mathbf{B}}_{t-i+1}{\mathbf{y}}_{t-i}+{\mathbf{ϵ}}_t,{\mathbf{ϵ}}_t\sim \mathcal{N}(0,\mathbf{\Sigma }).$$

$${\mathbf{y}}_t=\sum _{i=1}^p{\mathbf{X}}_{t-i+1}{\mathbf{\beta }}_{t-i+1}+{\mathbf{ϵ}}_t,{\mathbf{ϵ}}_t\sim \mathcal{N}(0,\mathbf{\Sigma }).$$

$${\mathbf{y}}_t={\mathbf{X}}_t{\mathbf{\beta }}_t+{\mathbf{ϵ}}_t,{\mathbf{ϵ}}_t\sim \mathcal{N}(0,\mathbf{\Sigma }).$$

$$P({\mathbf{\beta }}_t)=\prod _{j=1}^{n^2}P({\beta }_{j,t}),t=2,\dots ,T.$$

$${\beta }_{j,t}\sim {\pi }_{t}R({\beta }_{j,t})+(1-{\pi }_t)Q({\beta }_{j,t}),$$

$${\beta }_{j,t}|{\mu }_j,{\lambda }_j,{\pi }_t$$

$${\beta }_{j,t}|{\mu }_j,{\lambda }_j,{\pi }_t,{\tau }_0$$

$$P_t(\cdot )=\sum _{i=1}^{\infty }{w}_{i,t}{\delta }_{({\mathbf{\theta }}_{i,t})}(\cdot ),$$

$$v_{i,1}$$

$$z_{i,t}|v_{i,t}$$

$$v_{i,t+1}|z_{i,t}$$

$$P_t\sim DP(\alpha ,P_0).$$

$${\rho }_t(A)$$

$$a_{ti}=\frac{2(1+m_{ti})+\alpha }{(1+\alpha +m_{ti})(1+\alpha )(2+\alpha )},b_{ti}=\frac{\alpha -1}{1+\alpha }+a_{ti}$$

$$x\sim \mathcal{G}a(x|a,b)$$

$$x\sim \mathcal{I}\mathcal{G}(x|a,b)$$

$${\beta }_{j,t}|{\mu }_j,{\lambda }_j,{\pi }_t$$

$${\lambda }_j|{\tau }_j$$

$$({\mu }_j,{\tau }_j)|P_t$$

$$(P_2,\dots ,P_T)$$

$$P_0({\mu }_j,{\tau }_j)$$

$${\pi }_t$$

$${\beta }_{j,t}|{\mu }_j,{\tau }_j,{\pi }_t\sim {\pi }_{t}R({\beta }_{j,t})+(1-{\pi }_t)\mathcal{D}E({\beta }_{j,t}|{\mu }_j,{\tau }_j).$$

$$\mathbf{\Sigma }\sim \mathcal{I}\mathcal{W}(\nu ,\mathbf{\Psi })$$

$$\gamma :=\underset{t\in \mathbb{N}}{\inf }\mathbb{E}\left[\frac{1}{t+1}\log (\Vert {\mathbf{B}}_0{\mathbf{B}}_{-1}\cdots {\mathbf{B}}_{-t}\Vert )\right]$$

$${\mathbf{y}}_t=\sum _{i=0}^{\infty }{\mathbf{B}}_0{\mathbf{B}}_{t-1}\cdots {\mathbf{B}}_{t-i+1}{\mathbf{ϵ}}_{t-i}$$

$${\mathbf{y}}_{t+1}={\mathbf{B}}_{t+1}{\mathbf{y}}_t+{\mathbf{ϵ}}_{t+1},t\in \mathbb{Z}.$$

$$\begin{array}[]{cccccc}c=0&d=4.0&a_{1}=20.0&b_{1}=0.1&a_{0}=0.64&b_{0}=1.25,\\ \alpha=1.0&\eta=1.0&\nu=n+12&\bm{\Psi}=\mathbf{I}_{n}/n.\end{array}$$

$$\begin{array}{rl}P(\mathbf{B},\mathbf{U},\mathbf{V},\mathbf{Z},& \mathbf{D},\Gamma ,\mathbf{\pi },\mathbf{\mu },\mathbf{\tau },\mathbf{\lambda },\Sigma |\mathbf{Y})\propto L(\mathbf{Y}|\mathbf{B},\Sigma )\cdot P(\Sigma )\cdot \prod _{k=1}^{k^{\ast }}P({\mu }_k)P({\tau }_k)\\ & \cdot \prod _{t=2}^{T}P({\pi }_t)\prod _{j=1}^{n^2}P({\beta }_{j,t}|{\pi }_t,{\lambda }_j,\mathbf{\mu })P({\lambda }_j|\mathbf{\tau })P(d_{j,t}|{\mathbf{v}}_t,{\mathbf{u}}_t)P({\gamma }_{j,t}|{\pi }_t)\\ & \cdot \prod _{i=1}P(u_{i,t}|v_{i,t})P(v_{i,t}|z_{i,t-1})P(z_{i,t-1}|v_{i,t-1}),\end{array}$$

$$\begin{array}[]{cccc}\mathbf{B}_{t-1}&=\begin{bmatrix}0&0&0&0.8\\ 0&0&0.8&0.2\\ 0.8&0.2&0&0\\ 0&-0.4&0&0\end{bmatrix}&\quad\mathbf{D}_{t-1}&=\begin{bmatrix}0&0&0&2\\ 0&0&2&1\\ 2&1&0&0\\ 0&3&0&0\end{bmatrix}\\ \\ \mathbf{B}_{t}&=\begin{bmatrix}0&0&0&0.8\\ 0.2&0&0.8&-0.4\\ 0.2&0&0&0.8\\ 0&-0.4&0&0\end{bmatrix}&\quad\mathbf{D}_{t}&=\begin{bmatrix}0&0&0&2\\ 1&0&2&3\\ 1&0&0&2\\ 0&3&0&0\end{bmatrix}\\ \\ \mathbf{B}_{t+1}&=\begin{bmatrix}0&0&0&0\\ 0.2&0&0.8&-0.4\\ 0.8&0&0&0.8\\ 0.2&0.2&0&0\end{bmatrix}&\quad\mathbf{D}_{t+1}&=\begin{bmatrix}0&0&0&0\\ 1&0&2&3\\ 2&0&0&2\\ 1&1&0&0\end{bmatrix}.\end{array}$$

$$P(v_{k,t}|z_{k,t},z_{k,t-1},\dots )$$

$$\propto v_{k,t}^{z_{k,t-1}}(1-v_{k,t}{)}^{\alpha +m_{k,t-1}-z_{k,t-1}-1}\frac{m_{k,t}!}{z_{k,t}!(m_{k,t}-z_{k,t})!}{v}_{k,t}^{z_{k,t}}(1-v_{k,t}{)}^{m_{k,t}-z_{k,t}}$$

$$\cdot \left[\prod _{j:d_{j,t}=k,{\gamma }_{j,t}=1}{v}_{j,t}\right]\left[\prod _{j:d_{j,t}>k,{\gamma }_{j,t}=1}(1-v_{j,t})\right]$$

$$\propto \mathcal{B}e({\overline{\xi }}_t,{\overline{\alpha }}_t),$$

$${\overline{\xi }}_t$$

$${\overline{\alpha }}_t$$

Full text

In this paper, we provide an innovative Bayesian nonparametric (BNP) time-varying graphical framework for making inference in high-dimensional time series. We include a Bayesian nonparametric dependent prior specification on the matrix of coefficients and the covariance matrix by mean of a Time-Series DPP as in Nieto-Barajas et al., (2012). Following Billio et al., (2019), our hierarchical prior overcomes over-parametrization and over-fitting issues by clustering the vector autoregressive (VAR) coefficients into groups and by shrinking the coefficients of each group toward a common location. Our BNP time-varying VAR model is based on a spike-and-slab construction coupled with dependent Dirichlet Process prior (DPP) and allows to: (i) infer time-varying Granger causality networks from time series; (ii) flexibly model and cluster non-zero time-varying coefficients; (iii) accommodate for potential non-linearities.

In order to assess the performance of the model, we study the merits of our approach by considering a well-known macroeconomic dataset. Moreover, we check the robustness of the method by comparing two alternative specifications, with Dirac and diffuse spike prior distributions.

**Keywords: Bayesian Nonparametrics; Dependent Dirichlet process; Large vector autoregression; Sparsity; Time-Varying networks.

AMS 2000 subject classifications: Primary 62; secondary 91B84.

JEL Classification: C11, C32, C51, C53**

1 Introduction

Over the last decade, the availability of large datasets in economics and finance has allowed the introduction of high dimensional models. In particular, large datasets in macroeconomics help to improve the forecasts, while in finance some authors have investigated the use of large datasets to analyse financial crises, contagion effects and their impact on the real economy. In order to deal with high dimensional models, the introduction of Bayesian nonparametric techniques have become popular in different fields (such as statistics and machine learning), but only few attempts have been made in econometrics. In particular, Bayesian nonparametric approach allows to improve the estimation efficiency and the prediction accuracy in time series analysis.

Recently, time-varying parameter (TVP) models provide an interesting alternative to process multiple change points; for example, time-varying structural vector autoregressive (VAR) models have been used in Primiceri, (2005) for study monetary policy application; Dangl and Halling, (2012) forecast equity returns by mean of TVP models; and in Belmonte et al., (2014) the European inflation has been studied via a time-varying parameters model. As shown in Primiceri, (2005), Del Negro and Primiceri, (2015) and Bitto and Frühwirth-Schnatter, (2019), the advantage in capturing gradual changes is due to the flexibility of TVP models. We combine the ideas behind time-varying parameters models and Bayesian nonparametric techniques, thus allowing to model complex phenomena in a flexible and efficient manner. Moreover, we provide an innovative Bayesian nonparametric time-varying graphical framework for making inference in high-dimensional time series.

In this paper, we allow coefficients to be sparse, meaning that only a fraction of the time varying parameters have significant effects, but we retain flexibility in modelling non-zero coefficients, by including temporal dependence in the prior structure. In order to achieve these goals, we define a shrinkage prior on the VAR coefficients by means of a Bayesian nonparametric prior (BNP). This distribution is a spike-and-slab prior, where on the spike (parametric) distribution, we impose two different specifications: a Dirac spike and a “diffuse” spike. On the other hand, on the (non-parametric) slab distribution, we use a well know Bayesian nonparametric Lasso prior as in Billio et al., (2019).

The prior previously described groups the time-varying parameter vector autoregressive (TVP-VAR) coefficients into clusters and shrinks the coefficients within a cluster toward common notation. Differently from Markov-switching approach (Krolzig,, 1997) and random walk processes (Primiceri,, 2005; Del Negro and Primiceri,, 2015), we impose time variation on the distribution of the VAR coefficients. In the literature of time-varying coefficients, the VAR coefficients are represented as a direct dependence, in practice they can be represented as state-space models, where they are functions of the previous time. On the other hand, we introduce a different structure, the indirect dependence on the VAR coefficients. In this case, we have a dependence construct through the atoms of the Dirichlet process and not on a direct way. Thus, we include a Bayesian nonparametric dependent prior specification on the VAR coefficients and the covariance matrix by means of a time-series dependent Dirichlet process (tsDDP) as in Nieto-Barajas et al., (2012).

Following Billio et al., (2019), our hierarchical prior overcomes overparametrization and overfitting issues by clustering the VAR coefficients into groups and by shrinking the coefficients of each group toward a common location. This hierarchical prior allows to contemporaneously estimate the (potentially) sparse time-varying causal network structure and to cluster the corresponding coefficients. In our BNP-TVP-VAR model, time-varying coefficients allow to (i) estimate the temporal networks of contemporaneous and causal structures, (ii) identify different sources of time variation, from the size of shocks and/or the propagation mechanism, and (iii) accommodate for potential non-linearities.

We also contribute to the literature on financial and macroeconomic contagion (see Billio et al., (2012); Bianchi et al., (2019) and Barigozzi and Brownlees, (2019)) through the lens of Granger causality and graphs/networks representation. Our BNP prior is particularly suited for studying Granger causality from time series and in particular it allows to estimate the most significant time-varying dependence interrelationships between the variables of interest. As explained above, we can extract time-varying graphs by using the posterior random partition induced by the non-parametric (slab) distribution, which allows to cluster the edges into groups.

1.1 Literature

Since their introduction in macroeconomics (see Sims, (1980)), vector autoregessive (VAR) models have been extensively used in econometrics and time series statistics. Large VAR models have been used to analyse and forecast high-dimensional macroeconomic data (e.g., McCracken and Ng, (2016)) and financial panels (e.g., Barigozzi and Brownlees, (2019)). Moreover, in recent years VAR models have been used for studying financial and macroeconomic contagion (e.g., Cogley and Sargent, (2005), Stock and Watson, (2007), Diebold and Yilmaz, (2012) and Bianchi et al., (2019)). Although, VAR models have been extensively used for assessing the impact and spread of external shocks (i.e., to perform impulse-response analysis), forecasting, estimating networks from Granger-causal relationships and to study systemic risk and financial contagion (e.g., Diebold and Yilmaz, (2009), Billio et al., (2012) and Barigozzi and Brownlees, (2019)).

Despite being a potentially very flexible statistical tools, the high number of parameters and the typical limited length of standard macroeconomic datasets make unrestricted inference daunting as the cross-sectional size increases. This has favoured the use of penalised regression and Bayesian methods for dealing with the problem of over-parametrisation. The general idea is to use informative priors to shrink the unrestricted model towards a more parsimonious setting, thereby reducing parameter uncertainty and improving forecast accuracy (see Karlsson, (2013), Koop and Korobilis, (2010) for a survey).

In the Bayesian VAR (a.k.a. BVAR) literature, a plethora of different prior distributions have been proposed to perform sparse estimation (e.g., see Giannone et al., (2014)). Starting from the well-known Minnesota prior (see Doan et al., (1984), Litterman, (1986)), which specifies an objective prior on the coefficient and covariance matrices of a VAR, several parametric approaches have been developed exploiting hierarchical structures and finite mixtures (e.g., Kalli and Griffin, (2014), Gefang, (2014), Huber and Feldkircher, (2019), Kastner and Huber, (2018)).

Among the recent contributions for dealing with large dimensional models, we distinguish two approaches: the first attempts to reduce the size of the data to handle or to process during each step of the inferential algorithm, while the second is concerned with the reduction of the size of the parameter space. Within the first class, we mention the Bayesian compressed VAR of Koop et al., (2018), who tackled the dimensionality issue by using random projections to compress the data, and the Bayesian composite likelihood approach of Chan et al., (2018). On the other hand, Gefang et al., (2019) and Koop and Korobilis, (2018) adopted a variational Bayes approach for performing efficient approximate posterior inference in large parameter spaces. Also, Kastner and Huber, (2018) exploited factor models and hierarchical shrinkage priors for providing a parsimonious parametrisation of the covariance matrix which allows for equation-by-equation estimation. Additional contributions for estimating large VAR and VARMA models include Koop and Korobilis, (2013), Korobilis, (2016) and Chan et al., (2016).

In addition to the large cross-sectional dimensionality, also the temporal length of many economic and financial datasets is steadily increasing. Thus, the possible relations between different variables of interest can be described by static matrix of coefficients. This assumption can be elapsed by introducing a time-variation of the matrix of coefficients of the time series. In particular, the most common approach consists in specifying a process governing the evolution of the parameters of interest. According to the force driving this dynamics, we distinguish observation-driven and parameter-driven time-varying parameter (TVP) models. The first class is mainly represented by generalised autoregressive score models (GAS, see Creal et al., (2013)), while the second one includes Markov switching (e.g., Hamilton, (1989), Krolzig, (1997)), change point (e.g. Pesaran et al., (2006)) and random walk models (e.g. Del Negro and Primiceri, (2015), Primiceri, (2005)). These processes are able to describe parameters whose evolution is subject to switching regimes, structural breaks or smooth changes, respectively.

In the Bayesian and frequentist literature, the use of parametric models has been widely studied by applying different shrinkage methods (such as the Least Absolute Shrinkage and Selection Operator, known as LASSO). In particular, important papers focus on sparse and efficient estimation in high-dimensional datasets. However, more recently increasing attention is being devoted to the issue of over-shrinkage and to the modelling of non-zero coefficients (e.g., Giannone et al., (2018)). Consequently, there is an increasing need for adequate statistical tools capable of flexibly model the dynamics described by a VAR process, allowing for sparsity without incurring into over-shrinking.

In this paper, we aim to contribute to the growing literature on the use of Bayesian nonparametrics in time series analysis. In particular, Bayesian nonparametric techniques are widely in statistics, machine learning and data analysis as powerful tools for flexible modelling of complex data structure. Only recently, Bayesian nonparametrics has increased popularity in econometrics and in economic time series modelling to capture observation clustering effects (see e.g. Bassetti et al., (2014), Kalli and Griffin, (2018) and Billio et al., (2019)).

Up to our knowledge, our paper is the first to provide sparse Bayesian nonparametric VAR model when the coefficients are time-varying and the proposed two-stage prior specification can be easily extended to other classes, such as the seemingly unrelated regression (SUR) models. We propose a novel Bayesian nonparametric prior structure, which provides a sparse estimation of the coefficient matrix of a VAR model. This representation allows to manage the flexibility of non-zero entries and most importantly, to manage the time-variation in the matrix of coefficients through the atoms of the Dirichlet process and not through a state-space representation.

Our approach substantially differs from the existing literature in two aspects as described below. First, we consider a spike-and-slab prior distribution for each entry of the coefficient matrix, where on the spike we have a parametric prior specification by mean of Dirac or diffuse prior. On the hand, the slab component has random nonparametric prior. Second, we impose prior dependence on the coefficients by specifying a Markov process for their random distribution. As a by-product of the estimation procedure, we are able to extract a time series of dependent Granger-causality graphs. This shows how the BNP-TVP-VAR contributes to the literature on the estimation of time-varying networks from economic and/or financial series.

The paper is organized as follows. Section 2 introduces the modelling framework and presents the BNP-TVP prior structure, then Section 3 presents posterior approximation and describes how to extract Granger-causal time varying graphs from time series. Finally, Section 4 draws the conclusions.

2 A Bayesian Time-Varying VAR Model

2.1 TVP-VAR models

Let $n$ be the number of units in a dataset and $\mathbf{y}_{t}=(y_{1,t},\dots,y_{n,t})$ a vector of $n$ variables available at time $t$. A time-varying parameters vector autoregressive model of order $p$ (TVP-VAR($p$)) is defined as

[TABLE]

where $B_{t}$ is the $(n\times n)$ matrix of time-varying coefficients an $t=p,\ldots,T$ is the time period. We assume that the error terms $\bm{\epsilon}_{t}=(\epsilon_{1,t},\ldots,\epsilon_{n,t})^{\prime}$ are i.i.d. for $t$ with Gaussian distribution $\mathcal{N}(\mathbf{0},\bm{\Sigma})$. Eq. (1) can be written in the more compact form as

[TABLE]

where we define $\bm{\beta}_{t}=\operatorname{vec}(\mathbf{B}_{t})$; $\mathbf{X}_{t}=(\mathbf{y}_{t-1}^{\prime}\otimes\mathbf{I}_{n})$; $\otimes$ is the Kronecker product and $\operatorname{vec}(\cdot)$ the column-wise vectorization operator that stacks the columns of a matrix into a column vector.

2.2 Prior specification

Let us consider the problem of defining a flexible prior for a time-varying parameter model and we define the following TVP-VAR($p$) with $p$ equal to $1$ as

[TABLE]

In Eq. (3), the number of parameters of the $n$-dimensional TVP-VAR($1$) model is $(T-1)n^{2}+n(n+1)/2=O(n^{2})$, thus scales quadratically in $n$. In macroeconomic and financial applications, the number of variables of interest ranges from $n=3$ (small size model) to $n=20$ (large model) and even $n=100$ (huge model). This highlights the twofold need for shrinkage estimation methods and in particular, for the introduction of sparse estimation of the coefficient matrix. In fact, it is very hard both to provide a meaningful interpretation for a large VAR with full time-varying matrix $\mathbf{B}_{t}$ and to have an efficient and computationally feasible algorithm for an unrestricted estimation.

Motivated by this fact, we provide a prior distribution, which allows for sparse estimation in a time-varying parameter setting. For each coefficient of the matrix of parameters $B_{t}$, we introduce a mixture prior with independent location and scale parameters:

[TABLE]

where $P(\beta_{j,t})$ is the probability distribution of the vector matrix of coefficients (for example, we can choose it as a Double Exponential or Laplace distribution). One of the most successful and widespread approach in the Bayesian literature consists in the use of (independent) spike-and-slab prior distributions (e.g., Mitchell and Beauchamp, (1988), George and McCulloch, (1993), Smith and Kohn, (1996), George and McCulloch, (1997)) for each coefficient $\beta_{j,t}$. Based on it, we specify a spike-and-slab prior distribution for each $\beta_{j,t}$, with $j=1,\ldots,n^{2}$ and $t=2,\ldots,T$, of the form

[TABLE]

where $R,Q$ correspond to the spike and slab distributions, respectively, and $\pi_{t}$ is the time-varying mixing probability (i.e., the prior probability of the spike component). In the literature, we have two commonly choices for $D$: a Dirac mass at [math] such that $R(\beta_{j,t})=\delta_{\{0\}}(\beta_{j,t})$; and a centered (in zero) Normal distribution $R(\beta_{j,t})=\mathcal{N}(\beta_{j,t}|0,\tau_{0})$. The Dirac spike is a degenerate distribution that allows for variable selection as a by-product of the estimation. Instead, the choice of a continuous, diffuse prior (like a Gaussian) allows for shrinkage of the coefficients and is computationally faster, but requires the post-processing specification of threshold for the sake of variable selection.

The standard choice for the slab component $Q$ is a heavy-tailed distribution belonging to the family of Generalised Hyperbolic distribution (e.g., Double Exponential, Cauchy, $t$-Student), since the aim of this component is to capture potentially large non-zero coefficients. In the case of Dirac spike, the prior for each coefficient, for $j=1,\ldots,n^{2}$ and $t=2,\ldots,T$, is given by

[TABLE]

while for a diffuse (Gaussian) spike we have

[TABLE]

Example 1**.**

In Fig. 1 we report an example of spike-and-slab prior, with centred Gaussian spike distribution (in blue) and centred double exponential slab distribution (in red). From the left panel, which shows the two distributions, we can see that the Gaussian accounts for most of the prior mass on [math] while the double exponential governs the tails. This is reflected in the plot on the right, which shows that the mixture distribution (with equal weights) has fatter tails than the Gaussian and more mass in a neighborhood of [math] than the double exponential.

As previously described, we study the performance of the two different specification of the spike component by comparing the performances of the two constructions in extracting time-varying Granger causality networks from time series data.

The literature on time-varying parameter (TVP) models is vast. Some common parametric specifications include the threshold AR (TAR, e.g., Tong and Lim, (1980)), smooth transition AR (STAR, e.g.,Teräsvirta, (1994)), along with their multivariate generalisations, Markov switching process (e.g., Hamilton, (1989), Krolzig, (1997)), change point process (e.g. Pesaran et al., (2006)) and random walk process (e.g. Del Negro and Primiceri, (2015), Primiceri, (2005)). The choice of the particular specifications has been motivated by the intent to capture a particular feature of the dynamic evolution of the coefficients, such as changing regimes, structural breaks or smooth variations.

Differently from the existing literature on TVP-VAR models, we model the temporal dependence of the autoregressive parameters $\beta_{j,t}$ via assuming that the underlying prior (random) distributions evolve according to a discrete-time Markov process. A standard approach in Bayesian nonparametrics involves the specification of a Dirichlet Process (a.k.a. DP, see Ferguson, (1973)) or a Dirichlet Process mixture (a.k.a. DPM, see e.g., Lo, (1984)) prior for the distribution of the parameters of interest. The use of DP and related priors for a random probability measure $P$ allows for clustering of the variables $x_{i}\mathbin{\overset{iid}{\kern 0.0pt\sim}}P$.

We proceed by introducing the prior temporal dependence between the random measures $P_{2},\ldots,P_{T}$ via the time series Dirichlet Process (tsDDP) of Nieto-Barajas et al., (2012). As stated in the paper, for time series models, it is convenient to use dependence on the weights and common location. In practice, we apply a common discretization over the sequence of random measures, while the assumption of common weights and dependent location will lead to a discretization over the probability scale. In opposite to Taddy, (2010), who was working with equally spaced time points, we accommodate for unequal time points. In our analysis a latent binomial process to induce the desired correlation has been used, differently from the stick-breaking random probability measures as in Taddy, (2010), which use a beta autoregression on the fractions of the stick-breaking constructions by mean of two sets of latent variables.

By exploiting the stick-breaking construction of Sethuraman, (1994), the time series Dirichlet Process imposes a dependence for the random probability measures

[TABLE]

where the locations are fixed $\bm{\theta}_{i,t}=\bm{\theta}_{i}$ and the weights $w_{i,t}$ vary over time. The dependence is described by a Markov process for each un-normalised stick-breaking weight $v_{i,t}$, with $i=1,\ldots,\infty$, via auxiliary variables $z_{i,t}$ (in the spirit of Pitt et al., (2002), Pitt and Walker, (2005)), as follows

[TABLE]

The hyper-parameter $m_{i,t}$ tunes the strength of the dependence between $P_{t}$ and $P_{t+1}$, such that $m_{i,t}=0$ implies $P_{t}\perp P_{t+1}$, while $m_{i,t}\to\infty$ implies $P_{t}=P_{t+1}$ with probability 1 (see Nieto-Barajas et al., (2012)). Note that this construction implies that at each time $t=2,\ldots,T$, the marginal distribution of each random measure is a Dirichlet Process, that is

[TABLE]

with total mass parameter $\alpha$ and base measure $P_{0}$, such that the base measure defines the expectation and the mass parameter is interpreted as the precision parameter.

Eq. (8) explains the joint distribution of $z_{i,t}$ and $v_{i,t}$ and it allows us to define the joint model for $(P_{1},\ldots,P_{T})$ as a $tsDDP(\alpha,P_{0},\mathbf{m})$, where $\mathbf{m}$ is the sequence of the strength of dependence, $m_{i,t}$ for $i=1,2,\ldots$ and $t=1,\ldots,T$. In Figure 2, we show a single draw of $(P_{1},\ldots,P_{T})\sim tsDDP(\alpha,P_{0},\mathbf{m})$, with total mass parameter $\alpha$ equal to $10$; base measure $P_{0}$ as a normal distribution with zero mean and variance $\sqrt{2}$, i.e. $P_{0}\sim\mathcal{N}(0,\sqrt{2})$ and total timing $T=6$.

In order to assess the dependence structure induced by the time series Dirichlet Process, we consider the correlation between two random probability measures $P_{t}$ and $P_{t+1}$. The following proposition is explaining this correlation:

Proposition 2.1** (Nieto-Barajas et al., (2012)).**

Let $A\subset\mathbb{R}$ be measurable. For $(P_{1},\ldots,P_{T})\sim tsDDP(\alpha,P_{0},\mathbf{m})$ and any $t=1,\ldots,T-1$ let $\rho_{t}(A)\coloneqq\textnormal{Corr}(P_{t}(A),P_{t+1}(A))$. Then

[TABLE]

where

[TABLE]

Remark 1*.*

The correlation between $(P_{t},P_{t+1})$ is larger in regions where the prior mean $P_{0}$ assigns more probability, meaning that the $tsDDP$ prior places strongest dependence in $P_{0}$-most probable regions. Note that strong dependence between $(P_{t},P_{t+1})$ does not imply strong dependence between their outcomes $(\beta_{i,t},\beta_{j,t+1})$.

We can summarize what we have described above in the following prior structure111We use the shape-scale parametrisation of the Gamma distribution (thus $\mathbb{E}[x]=ab$ and $\mathbb{V}[x]=ab^{2}$) and Inverse Gamma distribution, whose probability density functions are, respectively

$\displaystyle x\sim\mathcal{G}a(x|a,b)$ $\displaystyle\iff p(x|a,b)=\frac{1}{b^{a}\Gamma(a)}x^{a-1}e^{-x/b},\qquad x\in(0,\infty),$

$\displaystyle x\sim\mathcal{IG}(x|a,b)$ $\displaystyle\iff p(x|a,b)=\frac{b^{a}}{\Gamma(a)}x^{-a-1}e^{-b/x},\qquad x\in(0,\infty).$

The exponential distribution is obtained as a particular case when $a=1$, that is $\mathcal{E}xp(x|b)=\mathcal{G}a(x|1,b)$., for $j=1,\ldots,n^{2}$ and $t=2,\ldots,T$ as follows

[TABLE]

where $R(\beta_{j,t})$ is either a Dirac mass at [math] or a Normal distribution centered in zero and with variance $\tau_{0}$, which marginally has Inverse Gamma prior distribution $\tau_{0}\sim\mathcal{IG}(\tau_{0}|a_{0},b_{0})$. If we marginalize over $\lambda_{j}$, we have a Double Exponential slab distribution for each entry of the coefficient matrix. Following the notation of Eq. (4) we have $Q(\beta_{j,t})=\mathcal{D}E(\beta_{j,t}|\mu_{j},\tau_{j})$, resulting in

[TABLE]

For the covariance matrix $\Sigma$, we assume an Inverse Wishart prior distribution that is:

[TABLE]

where $\nu$ and $\bm{\Psi}$ are the degrees of freedom and scale hyperparameters, respectively.

In summary, the observational model in Eq. (2) together with the prior structure in Eqs. (9) to (14) lead to the BNP-TVP-VAR(1) model. Eq. (9)–(14) represent our hierarchical prior and in Figure 3 we represent them through a Directed Acyclic Graph (DAG) for the Normal spike specification. The observable and non-observable random variables are indicated through shadow and empty circles, respectively. On the left side we have the prior for $\Sigma$, while on the right side, we have the hierarchial prior for $\beta_{t}$, with a description of the first and second stage of the hierarchy by means of the shrinking parameters $\mu$, $\tau$ and $\lambda$.

Sufficient conditions for stationarity of TVP autoregressive models are given for the univariate case (despite the proof is valid also in the multivariate setting) in Brandt, (1986), while Bourgerol and Picard, (1992) provides conditions for multivariate regressions where the coefficients are independent and identically distributed. The sufficient conditions given by Brandt, (1986) is reported below.

Theorem 2.1** (Brandt, (1986)).**

Let $\{(\mathbf{B}_{t},\bm{\epsilon}_{t}),\;t\in\mathbb{Z}\}$ be a strictly stationary ergodic process such that both $\mathbb{E}[\log^{+}(\left\lVert\mathbf{B}_{0}\right\rVert)]$ and $\mathbb{E}[\log^{+}(\left\lVert\bm{\epsilon}_{0}\right\rVert)]$ are finite. Suppose that the top Lyapunov exponent $\gamma$ defined by

[TABLE]

is strictly negative. Then, for all $t\in\mathbb{Z}$, the series

[TABLE]

converges a.s., and the process $\{\mathbf{y}_{t},\;t\in\mathbb{Z}\}$ is the unique strictly stationary solution of

[TABLE]

2.3 Hyper-parameter elicitation

Following Nieto-Barajas et al., (2012), we assume $m_{i,t}=m$, for each $i=1,2,\ldots$ and $t=1,\ldots,T$. Higher values of $m$ strengthen the dependence between the un-normalised weights $v_{i,t}$, however when big $m$ may induce the prior to overcome the likelihood, especially when the sample size is small. For this reason they specify a Poisson prior distribution for $m$, truncated on $\{1,\ldots,5\}$. Given the complexity of our prior specification, we prefer to fix the value of $m=5$, which is sufficiently small to avoid overweighting of the prior222In our empirical application, the sample size is $T=248$, while Nieto-Barajas et al., (2012) have $T=8$. and then check the robustness of the results to alternative values of $m$. We choose the following values for the hyper-parameters:

[TABLE]

This choice amounts to assuming a uniform prior on each $\pi_{t}$ and a rather uninformative prior on the covariance matrix, $\Sigma$. The value of the concentration parameter $\alpha$ is set according to standard practice in Dirichlet Process literature. The hyper-parameters $(c,d,a_{1},b_{1})$ imply that for each new component of the Dirichlet Process the prior distribution of $\mu_{j}$ is centered at zero mean with medium-high variance, whereas the prior for $\tau_{j}$ has mean $2$. Instead, the values of $(a_{0},b_{0})$ imply that the prior variance of the (diffuse) spike distribution is $0.8$, reflecting that this component should account for coefficients $\beta_{j,t}$ not significantly different from zero.

3 Posterior computation

3.1 Sampling method

Since the joint posterior distribution is not tractable and it is complex to be sample from, Bayesian estimator cannot be obtained analytically. In this paper, we rely on simulation based inference methods, and develop a Gibbs sampler algorithm for approximating the posterior distribution.

In order to deal with the finite mixture provided by the spike-and-slab prior and the infinite mixture given by the DPM, we exploited a data augmentation approach. For each $j=1,\ldots,n^{2}$ and $t=2,\ldots,T$, we introduce two sets of allocation variables $\gamma_{j,t},d_{j,t}$; a set of stick-breaking variables, $\mathbf{v}_{t}=\{v_{i,t}:i=1,2,\ldots\}$; a set of auxiliary variables $z_{i,t}$ (for $i=1,2,\ldots$) and a set of slice variables, $\mathbf{u}_{t}=\{u_{j,t}:j=1,\ldots,n^{2}\}$. The allocation variables, $\gamma_{j,t}$, selects the spike component $R(\cdot)$, when $\gamma_{j,t}$ is equal to zero and the slab component, when it is equal to one. The second allocation variable, $d_{j,t}$, selects the component of the Dirichlet mixture to which each single coefficient $\beta_{j,t}$ is allocated to. The sequence of stick-breaking variables defines the mixture weights, whereas the slice variable, $u_{j,t}$, allows us to deal with the infinite mixture components by identifying a finite number of stick-breaking variables to be sampled and an upper bound for the allocation variables $d_{j,t}$.

Finally, we obtain the following joint posterior distribution

[TABLE]

where $\mathbf{U}=\{u_{j,t}:j=1,\ldots,n^{2};\mbox{ and }t=2,\ldots,T\}$ and $\mathbf{V}=\{v_{i,t}:i=1,2,\ldots\mbox{ and }t=2,\ldots,T\}$ are the collections of slice variables and stick-breaking components, respectively; $\mathbf{Z}=\{z_{i,t}:i=1,2,\ldots\mbox{ and }t=2,\ldots,T\}$ and $\bm{\lambda}=\{\lambda_{j}:j=1,\ldots,n^{2}\}$ are the auxiliary and latent variables, respectively; $\mathbf{D}=\{d_{j,t}:j=1,\ldots,n^{2};\mbox{ and }t=2,\ldots,T\}$ and $\Gamma=\{\gamma_{j,t}:j=1,\ldots,n^{2};\mbox{ and }t=2,\ldots,T\}$ are the allocation variables; $(\bm{\mu},\bm{\tau})=\{(\mu_{k},\tau_{k}):k=1,\ldots,k^{\ast}\}$ are the atoms, where $k$ ranges from $1$ to the number $k^{*}$ of allocated DP components; $\mathbf{B}=\{\bm{\beta}_{t}:t=2,\ldots,T\}$ is the vector of VAR coefficients and $\bm{\pi}=\{\pi_{t}:t=2,\ldots,T\}$ are the specific probabilities of shrinking coefficients to zero.

We obtain random samples from the posterior distributions by Gibbs sampling. The Gibbs sampler is based on the algorithm of Hatjispyros et al., (2011) and on the slice sampler approach of Walker, (2007) and Kalli et al., (2011) for estimating the weights and locations of each random measure $P_{t}$. For improving the mixing of the MCMC, we introduced some Hamiltonian Monte Carlo (see Neal, (2011)) steps in spite of drawing from the full conditional posterior distribution. Hereafter, we show the iterative steps by using the conditional independence between variables, for $k=1,\ldots,k^{\ast}$, $j=1,\ldots,n^{2}$, $i=1,2,\ldots$ and $t=2,\ldots,T$:

1. (1)

the slice and stick-breaking variables $u_{j,t}$ and $v_{i,t}$ are updated along with the auxiliary variable $z_{i,t}$ given $\left[d_{j,t},\gamma_{j,t}\right]$; 2. (2)

the latent scale variables $\lambda_{j}$ are updated given $\left[\bm{\mu},\bm{\tau},(\beta_{j,t},d_{j,t},\gamma_{j,t})_{t}\right]$; 3. (3)

the parameters of the stick-breaking locations $(\mu_{k},\tau_{k})$ are updated given $\left[\bm{\lambda},\mathbf{B},\mathbf{D},\Gamma\right]$; 4. (4)

the allocation variables $d_{j,t},\gamma_{j,t}$ are jointly updated given $\left[\mu_{k},\tau_{k},\beta_{j,t},u_{j,t},v_{i,t},\pi_{t}\right]$; 5. (5)

the VAR coefficients $\bm{\beta}_{t}$ are jointly updated given $\left[\bm{\mu},\bm{\tau},\bm{\lambda},\Sigma,(d_{j,t},\gamma_{j,t})_{j},\mathbf{y}_{t}\right]$; 6. (6)

The covariance matrix $\Sigma$ is updated given $\left[\mathbf{B},\mathbf{Y}\right]$; 7. (7)

the mixing probability $\pi_{t}$ of having sparse coefficients is updated given $\left[(\gamma_{j,t})_{j}\right]$.

The detailed Gibbs sampler is described in Appendix A and Appendix C.

3.2 Graph extraction

Based on the Gibbs sampler previously described, we are able to extract time-varying Granger-causal graphs. In the literature, linkages and networks describing the relationships between variables of interest, such as macroeconomics and financial linkages (e.g. Billio et al., (2012) and Barigozzi and Brownlees, (2019)) can be used to extract pairwise Granger causality. This approach is generating spurious causality effects and does not consider conditioning on variables of interest. The main problem relies on the high number of variables available relative to the number of data, thus it could lead to overparametrization and inefficiency in gauging the causal relationships. Our proposed prior can be used to extract the networks and pairwise Granger causality while reducing the overfitting and curse of dimensionality problems. Moreover, the introduction of our prior could lead to the extraction of edge-colored graphs, that allows us to identify stylized facts in financial or macroeconomics networks and to show the presence of communities, hubs and linkage heterogeneity.

From the MCMC output of the time-varying coefficient matrix $\mathbf{B}_{t}$, we are able to extract time-varying Granger-causal graphs. At each time $t=2,\ldots,T$, we use the posterior random partition induced by the nonparametric (slab) distribution to cluster the edges of the graph (i.e., the entries $\beta_{j,t}$, $j=1,\ldots,n^{2}$ of the vectorised coefficient matrix $\bm{\beta}_{t}$) into groups.

Formally, a graph $G$ is a pair $(V,E)$, where $V$ is a set of nodes and $E$ is a set of nodes pairs, named links or edges. The nodes are labeled and a link/edge is identified by the pair of nodes it connects, $(i,j)$. In particular, we have the existence of an edge if and only if the time-varying VAR coefficients of the variable $y_{i,t-1}$ in the equation of $y_{j,t}$ is not null. In our network analysis, we focus on the adjacency matrix constructed a posterior from the allocation variables and it allows to take both values between [math] and $1$ if we apply a threshold, while if the values are allowed to vary between [math] and $1$, we have a weighted graph. The purpose is to estimate the most significant time-varying dependence interrelationships (in terms of Granger-causality) between the $n$ variables of interest.

Example 2**.**

Consider the TVP-VAR(1) model in Eq. (2) and let $n=4$. Without loss of generality, focus on the coefficient matrices at three consecutive times $t-1,t$ and $t+1$. Suppose the posterior estimates of the coefficient matrices and allocation variables $d_{j,t}$, respectively, are as follows

[TABLE]

The corresponding Granger-causal graphs are given in Fig. 4, where colours have been used to denote the cluster assignment encoded in the matrices $\mathbf{D}_{t-1},\mathbf{D}_{t}$ and $\mathbf{D}_{t+1}$.

4 Conclusions

We proposed the BNP-TVP-VAR model for sparse, nonparametric inference in time-varying VAR models. The use of spike-and-slab priors with time-series dependent Dirichlet Process prior for the slab component allows to contemporaneously shrink the autoregressive coefficients and flexibly modelling time-varying non-zero entries. We applied the proposed methodology using two alternative spike distributions: a Dirac and a Normal distribution. The performance of the resulting models has been compared in terms of: (i) sparse estimation and variable selection, and (ii) clustering structure. Moreover, we showed how the BNP-TVP-VAR model can be used for extracting Granger-causal time-dependent graphs from multivariate time series.

Appendix A Posterior distributions: diffuse DE spike

A.1 Posterior for stick breaking un-normalised weights $v_{i,t}$

Posterior distribution for $v_{i,t}$, for all $t=1,\ldots,T$ and all $k=1,\ldots,k^{*}$, where $k^{*}$ is the number of ties. We use the convention $z_{i,0}=m_{i,0}=0$, for all $i=1,2,\ldots$.

[TABLE]

where

[TABLE]

The normalised weights are then computed by $w_{k,t}=v_{k,t}\prod_{j<k}(1-v_{j,t})$. The posterior distribution of the latent auxiliary variables $z_{k,t}$, for all $t=1,\ldots,T-1$ and all $k=1,\ldots,k^{*}$ is given by

[TABLE]

Finally, the auxiliary variable for the slice sampler has posterior distribution given by

[TABLE]

A.2 Posterior for $\lambda_{j}$

Posterior distribution for $\lambda_{j}$, for $j=1,\ldots,n^{2}$. Define $\mu_{k_{j}}^{*},\tau_{k_{j}}^{*}$ are the location and scale of $\beta_{j,t}$, respectively, when sparse or non-sparse component is chosen.

[TABLE]

where

[TABLE]

The probability density function of a generalised inverse Gaussian, for $p\in\mathbb{R},a>0,b>0$ and $K_{p}(\cdot)$ being a modified Bessel function of the second kind, is

[TABLE]

The vector $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{n^{2}})$ represents the diagonal of the diagonal matrix $\bm{\Lambda}$ of size $(n^{2}\times n^{2})$.

A.3 Posterior of the stick breaking locations $(\mu_{k},\tau_{k})$

The posterior distribution for the stick breaking locations in the sparse case are given by:

[TABLE]

and

[TABLE]

with

[TABLE]

where $m_{0}=\sum_{(j,t)}\mathbb{I}(\gamma_{j,t}=0)=\sum_{j=1}^{n^{2}}\sum_{t=1}^{T}(1-\gamma_{j,t})$ and we used the parametrisation of the Gamma distribution with shape $a>0$ and scale $b>0$, that is

[TABLE]

Regarding the sample in the non-sparse case, we generate $(\mu_{k},\tau_{k})$ for $k=1,\ldots,N^{\ast}$ and we have the following full conditional separately for $\mu_{k}$ and $\tau_{k}$. The posterior distribution for the stick breaking locations in the non-sparse case, is given by:

[TABLE]

where

[TABLE]

and for the scale parameter we have the following representation:

[TABLE]

with

[TABLE]

where $m_{1}=\sum_{(j,t)}\mathbb{I}(\gamma_{j,t}=1)\mathbb{I}(d_{j,t}=k)=\sum_{j=1}^{n^{2}}\sum_{t=1}^{T}\gamma_{j,t}\mathbb{I}(d_{j,t}=k)$.

A.4 Posterior of the allocation variables $(d_{j,t},\gamma_{j,t})$

Regarding the allocation variables $d_{j,t},\gamma_{j,t}$, we obtain the following full conditional for the sparse and non-sparse case. Let $u_{j,t}\sim\mathcal{U}([0,1])$ for all $j=1,\ldots,n^{2}$ and all $t=1,\ldots,T$, be the auxiliary slice sampling variable. In the non-sparse case, we have:

[TABLE]

where $A_{w}(u_{j,t})=\{k\in\{1,2,\ldots,k^{*}\}:w_{k}>u_{j,t}\}$, where $k^{*}$ is the number of ties. On the other hand, in the sparse case we have the non-normalised posterior probability

[TABLE]

and full conditional distribution has the following representation:

[TABLE]

A.5 Posterior for $\beta_{j,t}$

Posterior distribution for $\bm{\beta}_{t}$, for all $t=2,\ldots,T$. Let $\bm{\mu}_{t}^{*}=(\mu_{d_{1,t}},\dots,\mu_{d_{n^{2},t}})^{\prime}$ and $\bm{\gamma}_{t}=(\gamma_{1,t},\dots,\gamma_{n^{2},t})^{\prime}$, and denote the Hadamard product by $\odot$. We have:

[TABLE]

where

[TABLE]

A.6 Posterior for covariance matrix $\Sigma$

Posterior distribution for the covariance matrix $\bm{\Sigma}$.

[TABLE]

where

[TABLE]

A.7 Posterior for mixing probability $\pi_{t}$

Posterior distribution for the mixing probability $\pi_{t}$, for all $t=1,\ldots,T$.

[TABLE]

where

[TABLE]

Appendix B Posterior distributions: diffuse Normal spike

B.1 Posterior for stick breaking un-normalised weights $v_{i,t}$

Posterior distribution for $v_{i,t}$, for all $t=1,\ldots,T$ and all $k=1,\ldots,k^{*}$, where $k^{*}$ is the number of ties. We use the convention $z_{i,0}=m_{i,0}=0$, for all $i=1,2,\ldots$.

[TABLE]

where

[TABLE]

The normalised weights are then computed by $w_{k,t}=v_{k,t}\prod_{j<k}(1-v_{j,t})$. The posterior distribution of the latent auxiliary variables $z_{k,t}$, for all $t=1,\ldots,T-1$ and all $k=1,\ldots,k^{*}$ is given by

[TABLE]

Finally, the auxiliary variable for the slice sampler has posterior distribution given by

[TABLE]

B.2 Posterior for $\lambda_{j}$

Posterior distribution for $\lambda_{j}$, for $j=1,\ldots,n^{2}$. Define $\mu_{k_{j}}^{*},\tau_{k_{j}}^{*}$ are the location and scale of $\beta_{j,t}$, respectively, when sparse or non-sparse component is chosen.

[TABLE]

where

[TABLE]

The probability density function of a generalised inverse Gaussian, for $p\in\mathbb{R},a>0,b>0$ and $K_{p}(\cdot)$ being a modified Bessel function of the second kind, is

[TABLE]

The vector $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{n^{2}})$ represents the diagonal of the diagonal matrix $\bm{\Lambda}$ of size $(n^{2}\times n^{2})$.

B.3 Posterior of the stick breaking locations $(\mu_{k},\tau_{k})$

The posterior distribution for the stick breaking locations in the sparse case are given by:

[TABLE]

and

[TABLE]

with

[TABLE]

where $m_{0}=\sum_{(j,t)}\mathbb{I}(\gamma_{j,t}=0)=\sum_{j=1}^{n^{2}}\sum_{t=1}^{T}(1-\gamma_{j,t})$ and we used the parametrisation of the Inverse Gamma distribution with shape $a>0$ and scale $b>0$, that is

[TABLE]

Regarding the sample in the non-sparse case, we generate $(\mu_{k},\tau_{k})$ for $k=1,\ldots,N^{\ast}$ and we have the following full conditional separately for $\mu_{k}$ and $\tau_{k}$. The posterior distribution for the stick breaking locations in the non-sparse case, is given by:

[TABLE]

where

[TABLE]

and for the scale parameter we have the following representation:

[TABLE]

with

[TABLE]

where $m_{1}=\sum_{(j,t)}\mathbb{I}(\gamma_{j,t}=1)\mathbb{I}(d_{j,t}=k)=\sum_{j=1}^{n^{2}}\sum_{t=1}^{T}\gamma_{j,t}\mathbb{I}(d_{j,t}=k)$.

B.4 Posterior of the allocation variables $(d_{j,t},\gamma_{j,t})$

Regarding the allocation variables $d_{j,t},\gamma_{j,t}$, we obtain the following full conditional for the sparse and non-sparse case. Let $u_{j,t}\sim\mathcal{U}([0,1])$ for all $j=1,\ldots,n^{2}$ and all $t=1,\ldots,T$, be the auxiliary slice sampling variable. In the non-sparse case, we have:

[TABLE]

where $A_{w}(u_{j,t})=\{k\in\{1,2,\ldots,k^{*}\}:w_{k}>u_{j,t}\}$, where $k^{*}$ is the number of ties. On the other hand, in the sparse case we have the non-normalised posterior probability

[TABLE]

and full conditional distribution has the following representation:

[TABLE]

B.5 Posterior for $\beta_{j,t}$

Posterior distribution for $\bm{\beta}_{t}$, for all $t=2,\ldots,T$. Let $\bm{\mu}_{t}^{*}=(\mu_{d_{1,t}},\dots,\mu_{d_{n^{2},t}})^{\prime}$ and $\bm{\gamma}_{t}=(\gamma_{1,t},\dots,\gamma_{n^{2},t})^{\prime}$, and denote the Hadamard product by $\odot$. We have:

[TABLE]

where

[TABLE]

and $\bm{\Lambda}_{\tau_{0},\gamma}$ is a $n^{2}\times n^{2}$ diagonal matrix whose $jj$-th entry is

[TABLE]

B.6 Posterior for covariance matrix $\Sigma$

Posterior distribution for the covariance matrix $\bm{\Sigma}$.

[TABLE]

where

[TABLE]

B.7 Posterior for mixing probability $\pi_{t}$

Posterior distribution for the mixing probability $\pi_{t}$, for all $t=1,\ldots,T$.

[TABLE]

where

[TABLE]

Appendix C Posterior distributions: Dirac spike

C.1 Posterior for stick breaking un-normalised weights $v_{i,t}$

Posterior distribution for $v_{i,t}$, for all $t=2,\ldots,T$ and all $k=1,\ldots,k^{*}$, where $k^{*}$ is the number of ties. We use the convention $z_{i,0}=m_{i,0}=0$, for all $i=1,2,\ldots$.

[TABLE]

where

[TABLE]

The normalised weights are then computed by $w_{k,t}=v_{k,t}\prod_{j<k}(1-v_{j,t})$. The posterior distribution of the latent auxiliary variables $z_{k,t}$, for all $t=1,\ldots,T-1$ and all $k=1,\ldots,k^{*}$ is given by

[TABLE]

Finally, the auxiliary variable for the slice sampler has posterior distribution given by

[TABLE]

C.2 Posterior for $\lambda_{j}$

Posterior distribution for $\lambda_{j}$, for $j=1,\ldots,n^{2}$. Define $\mu_{k_{j}}^{*},\tau_{k_{j}}^{*}$ are the location and scale of $\beta_{j,t}$, respectively, when sparse or non-sparse component is chosen.

[TABLE]

where

[TABLE]

The probability density function of a generalised inverse Gaussian, for $p\in\mathbb{R},a>0,b>0$ and $K_{p}(\cdot)$ being a modified Bessel function of the second kind, is

[TABLE]

The vector $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{n^{2}})$ represents the diagonal of the diagonal matrix $\bm{\Lambda}$ of size $(n^{2}\times n^{2})$.

C.3 Posterior of the stick breaking locations $(\mu_{k},\tau_{k})$

In the non-sparse case, we generate $(\mu_{k},\tau_{k})$ for $k=1,\ldots,k^{*}$ and we have the following full conditional separately for $\mu_{k}$ and $\tau_{k}$. The posterior distribution for the stick breaking locations in the non-sparse case, is given by:

[TABLE]

where

[TABLE]

and for the scale parameter we have the following representation:

[TABLE]

with

[TABLE]

where $m_{1}=\sum_{(j,t)}\mathbb{I}(\gamma_{j,t}=1)\mathbb{I}(d_{j,t}=k)=\sum_{j=1}^{n^{2}}\sum_{t=2}^{T}\gamma_{j,t}\mathbb{I}(d_{j,t}=k)$.

C.4 Posterior of the allocation variable $(d_{j,t},\gamma_{j,t})$

Regarding the allocation variables $(d_{j,t},\gamma_{j,t})$, we obtain the following full conditional for the sparse and non-sparse case. The (conditional) prior distribution for the allocation variable $d_{j,t}$, for each $k=1,\dots,k^{*}$, $j=1,\ldots,n^{2}$, $t=2,\ldots,T$, is given by:

[TABLE]

while the prior for $\gamma_{j,t}$ is

[TABLE]

Define $\bm{\gamma}_{-j,t}=(\gamma_{1,t},\ldots,\gamma_{j-1,t},\gamma_{j+1,t},\ldots,\gamma_{n^{2},t})^{\prime}$, $\mathbf{d}_{-j,t}=(d_{1,t},\ldots,d_{j-1,t},d_{j+1,t},\ldots,d_{n^{2},t})^{\prime}$. The joint posterior distribution of the allocation variables $(d_{j,t},\gamma_{j,t})$, for each $k=1,\dots,k^{*}$, $j=1,\ldots,n^{2}$, $t=2,\ldots,T$, is obtained as:

[TABLE]

where the (conditional) marginal likelihood obtained by integrating out the $\bm{\beta}_{t}$, for each $t=2,\ldots,T$, is given by

[TABLE]

where $\mathbf{X}_{\bm{\gamma},\mathbf{d},t}$ contains the columns of $\mathbf{X}_{t}$ corresponding to the $j$ such that $\gamma_{j,t}=1$ and $d_{j,t}=\mathbf{d}_{j}$ (similarly for $\bm{\Lambda}_{\bm{\gamma},\mathbf{d},t}$), $\bm{\mu}_{\bm{\gamma},\mathbf{d},t}$ is the sub-vector of $(\mu_{d_{1,t}},\ldots,\mu_{d_{n^{2},t}})^{\prime}$ containing the elements $j$ such that $\gamma_{j,t}=1$, and

[TABLE]

C.5 Posterior for $\beta_{j,t}$

Posterior distribution for $\bm{\beta}_{t}$, for all $t=2,\ldots,T$. In the sparse case we have

[TABLE]

In the non-sparse case, denoting $\bm{\beta}_{\bm{\gamma},t}$ the sub-vector of $\bm{\beta}_{t}$ corresponding to the coefficients such that $\gamma_{j,t}=1$, we have

[TABLE]

where

[TABLE]

C.6 Posterior for covariance matrix $\Sigma$

Posterior distribution for the covariance matrix $\bm{\Sigma}$.

[TABLE]

where

[TABLE]

C.7 Posterior for mixing probability $\pi_{t}$

Posterior distribution for the mixing probability $\pi_{t}$, for all $t=2,\ldots,T$.

[TABLE]

where

[TABLE]