Nonparametric Bayes dynamic modeling of relational data

Introduction

Relational data often take the form of a symmetric binary matrix, with entries indicating the presence or absence of links between pairs of individuals or entities. In dynamic settings, the links and the set of entities under consideration can change over time, and interest focuses on inferences on the time varying relational structure and in prediction. Examples include social network analysis, in which links encode friendship networks among individuals, and broader relational settings in which closeness between a pair of units (products, stimuli, countries, companies, etc) is expressed on a binary scale. Figure 1 shows an example of time-varying binary similarity matrices encoding dynamic co-movements in National Stock Market Indices from 2004 to 2013. Co-movements among a set of assets or market indices are typically analyzed via time-varying covariance or correlation matrices of their corresponding log-returns $Z_{t}=[z_{1,t},z_{2,t},...,z_{V,t}]^{\prime}$, $t=1,...,T$, (see e.g. Tsay, 2005, Wilson & Ghahramani, 2010, Nakajima & West, 2012, Durante et al., 2013); we instead provide a different and not yet fully explored direction of research by treating co-movements as dynamic relational data, shifting our attention from $Z_{t}$ to the $V\times V$ time-varying symmetric matrices $\{Y_{t},\ t\in\mathcal{T}\subset\Re^{+}\}$. The matrix $Y_{t}$ has entries $y_{ij,t}=y_{ji,t}=1$ if index $i$ and index $j$ move in the same direction at time $t$ (i.e. $z_{i,t}>0$ and $z_{j,t}>0$, or $z_{i,t}<0$ and $z_{j,t}<0$) and $y_{ij,t}=y_{ji,t}=0$ if they move in opposite directions (i.e. $z_{i,t}>0$ and $z_{j,t}<0$, or $z_{i,t}<0$ and $z_{j,t}>0$). Co-movements indicate similarity in the indices.

A rich literature is available on modeling similarity or dissimilarity matrices, with Multidimensional Scaling (MDS) providing a widely used technique for graphically representing units in a Euclidean space conditionally on their pairwise dissimilarity measures. General theory and applications are available for Euclidean distances and rank dissimilarities (see Cox & Cox, 2001), with subsequent developments in a Bayesian framework (Oh & Raftery, 2001, Oh & Raftery, 2007) improving the overall performance, but subject to possible issues due to non-identifiable latent coordinates, lack of full conditional conjugacy and absence of an automatic procedure for learning the dimension of the latent space. Moreover, generalizations in the dynamic case are lacking, with only few recent proposals restricted to specific applications for discrete time evolution (Jamali-Rad & Leus, 2012).

When binary similarity or dissimilarity matrices are analyzed, the previous procedures prove to be inappropriate or impractical (Holbrook et al., 1982), with predicted values outside the probability range and a large number of tied ranks for each unit in non-metric MDS applications. Spatial analysis of choice data (DeSarbo & Hoffman, 1987, DeSarbo et al., 1999) provides a possible generalization of MDS for binary variables, with recently developed algorithms available also in the dynamic case (Sarkar et al., 2007). However, questionable independence assumptions are required to ease maximum likelihood estimation, and Bayesian extensions (DeSarbo et al., 1999) to overcome this problem lack scalability in selecting the dimensionality of the latent space via cross-validation methods. Moreover, dynamic extensions via the Kalman filter rely on first and second order Taylor expansions for the observation model, providing difficulties in the derivation of theoretical properties for the exact formulation and requiring a sufficient number of observations to meet the Gaussian assumption. These models are specifically tailored for embedding problems in 2-mode co-occurrence data recording links between two different types of entities (i.e. consumer-products, author-words). Our focus is instead on dynamic modeling for one-mode binary matrices.

There is a growing body of literature in social networks on model-based statistical analysis of one-mode binary matrices, traditionally focusing on overly-restrictive models, such as Erdös & Rény (1959), the $p_{1}$ model (Holland & Leinhardt, 1981) and the Exponential Random Graph Model (ERGM) (Frank & Strauss, 1986), with generalizations for dynamic inference available via discrete temporal ERGM (Robins & Pattinson, 2001) and hidden temporal ERGM (Guo et al., 2007). ERGMs have had growing popularity, but have a number of drawbacks. Estimation relies on pseudo-likelihood (Strauss & Ikeda, 1990) and approximate MCMC methods (Snijders, 2002), due to the computational intractability in a fully likelihood approach. Solutions can be degenerate or nearly-degenerate (Handcock et al., 2003), and questions remain about coherence, inflexibility and other key issues.

An alternative class of models focus on clustering the nodes, based on the pattern of inter-connections in the network. Stochastic Block Models (SBM) (Nowicki & Snijders, 2001) provide a common framework, with the Infinite Relational Model (IRM) (Kemp et al., 2006) allowing an unknown number of clusters via a Dirichlet process. Dynamic SBMs have been recently considered (Ishiguro et al., 2010,Yang et al., 2011, Xu & Hero, 2013). Ishiguro et al. (2010) focus on discrete dynamic evolution via a hidden Markov model. Xu & Hero (2013) accommodate continuous time analysis via a state space formulation, but require sufficient numbers of observations in each block to meet Gaussian assumptions for the sample mean. They use the extended Kalman filter to linearize the observation equation, leading to questions of accuracy.

We dynamically model binary relational matrices by embedding the nodes in a low-dimensional latent Euclidean space, with coordinates evolving in continuous time via Gaussian processes and edge probabilities constructed utilizing a logistic mapping function from the probability matrix space to the dot product of the latent coordinates. Hence, we are most closely related to the literature on latent space models (Hoff et al., 2002) and Mixed Membership Stochastic Block models (MMSB) (Airoldi et al., 2008), which allow each node to belong to multiple blocks with fractional membership. Dynamic latent space models (Sarkar & Moore, 2005) and MMSB models (Xing et al., 2010) incorporate Gaussian perturbations in discrete time and state space models, respectively. Posterior computation relies on several layers of approximation without theory available to justify accuracy. In contrast, we provide a simple Gibbs sampling algorithm for our model, which converges to the exact posterior and infers the dimension of the latent space automatically.

The paper is organized as follows. In Section 2, we describe the general model structure with particular attention to prior specification and theoretical properties. Section 3 provides the Gibbs sampling steps. A simulation study is examined in Section 4, and an application to quarterly co-movements in world financial markets is presented in Section 5.

Dynamic Latent Space Model

Let $Y_{t}$ be the symmetric binary similarity matrix at time $t\in\mathcal{T}$ and $\pi(t)$ be the corresponding symmetric probability matrix having entries $\pi_{ij}(t)=\pi_{ji}(t)=\mbox{pr}(y_{ij,t}=1)$ for every $i=1,\dots,V$ and $j=1,\dots,V$. Letting

independently for each $i=2,\ldots,V$ and $j=1,\ldots,i-1$, our aim is to define a prior $\Pi_{\pi}$ for the collection of time-varying probability matrices $\pi_{\mathcal{T}}=\{\pi(t),t\in\mathcal{T}\}$ with the goals being to (i) obtain a provably flexible specification, (ii) maintain simple computations, (iii) perform dimensionality reduction in order to scale to moderately large $V$ settings, (iv) allow unequal spacing and missing observations and (v) allow predictions including a measure of predictive uncertainty. Since the matrices are symmetric and the similarities or dissimilarities of a unit with itself are meaningless, we will focus on modeling the lower triangular part without taking into account the diagonal elements.

2 Latent space dynamic model formulation

We construct $\pi_{ij}(t)$ via a monotonic increasing link function $g(\cdot):\Re\rightarrow$ mapping a latent similarity measure among units $i$ and $j$ at time $t$, $s_{ij}(t)\in\Re$, into the probability space. Specifically, we choose $g(\cdot)$ to be the distribution function of the logistic random variable, obtaining

for $i=2,\ldots,V$, $j=1,\ldots,i-1,$ and $t\in\mathcal{T}$. Without further assumptions on $s_{ij}(t)$, one needs to model separately $\frac{1}{2}V(V-1)$ stochastic processes, one for each time-varying similarity measure $s_{ij}(t)$, with $i=2,...,V$, $j=1,...,i-1$ and $t\in\mathcal{T}$, leading to burdensome computations as $V$ increases and failing to borrow information exploiting the underlying process inducing similarities among the units. In order to reduce the dimensionality of the problem and to learn also the network structure among the units for every $t$, we express the similarity measures $s_{ij}(t)$ as a quadratic combination of a set of latent coordinates for unit $i$ and unit $j$. Specifically

where $x_{i}(t)=[x_{i1}(t),\ldots,x_{iH}(t)]^{\prime}$ for $i=2,\ldots,V$ and $x_{j}(t)=[x_{j1}(t),\ldots,x_{jH}(t)]^{\prime}$ for $j=1,\ldots,i-1$, are the vectors of latent coordinates of unit $i$ and $j$ respectively, giving rise, together with the baseline $\mu(t)$, to the similarity measure between the two units via a projection approach. According to this specification, units with latent coordinates in the same directions will have a higher probability of being similar (i.e. $y_{ij,t}=1$), while units with opposite coordinates are more likely to be dissimilar (i.e. $y_{ij,t}=0$).

This formulation is also intuitive in practical applications. Recall our motivating example of finance, and assume for simplicity $\mu(t)=0$ and only two latent coordinates representing for example unexpected inflation and industrial production, respectively. Then indices of countries with features in the same directions will have a higher probability of co-moving, while countries with opposite unexpected inflation and industrial production will more likely move on different directions.

In matrix notation, equation (3) can be rewritten as

where $S(t)$ is a $V\times V$ real symmetric matrix with latent similarity entries $s_{ij}(t)$ and $X(t)=[x_{1}(t),x_{2}(t),\dots,x_{V}(t)]^{\prime}$. Note that, assuming without loss of generality $\mu(t)=0$, the above decomposition is not unique. For example if we define $X(t)^{*}=X(t)Q$ with $Q$ a $H\times H$ orthogonal matrix, then $X(t)^{*}X(t)^{*^{\prime}}=X(t)QQ^{\prime}X(t)^{\prime}=X(t)X(t)^{\prime}$. If one is interested also in making inference on the latent coordinates matrix $X(t)$, different proposals are available in latent factor modeling to ensure identifiability via restrictions (see e.g. Bollen, 1989) or Procrustean transformations (Hoff et al., 2002). However since our focus is on making inference and prediction on the probability matrices, we follow Ghosh & Dunson (2009) in avoiding identifiability constraints, as such constraints are not necessary to ensure identifiability of the induced similarity matrix $S(t)$.

It is important to characterize the class of $\pi(t)$ matrices whose lower triangular elements can be represented as in (2) with latent similarities decomposed as in (4). Theorem 1 and the corresponding Corollary 2 state that for $H$ sufficiently large, the lower triangular matrix elements of any symmetric probability matrix have such a representation. For $H\geq V$, $\mathcal{X}_{X}$ denotes the space of all $V\times H$ dimensional matrices of arbitrary coordinate functions mapping from $\mathcal{X}\rightarrow\Re$ and $\mathcal{X}_{\mu}$ the space of all baseline mean functions.

Given a symmetric real matrix $S(t)$, $\forall\ t\in\mathcal{T}$, there exist $\{X(t),\mu(t)\}\in\mathcal{X}_{X}\otimes\mathcal{X}_{\mu}$ such that

Proof. Assume without loss of generality that $\mu(t)=0$ and take $H\geq V$. Consider

where $P(t)$ is the matrix of the eigenvectors of $S(t)$ and $\Lambda(t)$ the diagonal matrix with the corresponding eigenvalues. Then $S(t)=P(t)\Lambda(t)P(t)^{\prime}=X(t)X(t)^{\prime}$, for every $t\in\mathcal{T}$.

Given a symmetric probability matrix $\pi(t)$, $\forall\ t\in\mathcal{T}$, there exist $\{X(t),\mu(t)\}\in\mathcal{X}_{X}\otimes\mathcal{X}_{\mu}$ such that

Proof. The proof follows immediately from Theorem 1 and from the fact that the mapping from $s_{ij}(t)$ to $\pi_{ij}(t)$ is a one-to-one continuous increasing function.

This ensures that our specification is sufficiently flexible to characterize any true generating process, and hence can be viewed as nonparametric given sufficiently flexible priors for the components.

3 Prior Specification

We aim to specify independent prior distributions $\Pi_{X}$ and $\Pi_{\mu}$ for $X_{\mathcal{T}}=\{X(t),t\in\mathcal{T}\}$ and $\mu_{\mathcal{T}}=\{\mu(t),t\in\mathcal{T}\}$ in order to induce a prior $\Pi_{\pi}$ for $\pi_{\mathcal{T}}=\{\pi(t),t\in\mathcal{T}\}$ through (2) and (3). This prior is carefully defined to have large support, favor simple and efficient computation, allow missing values, induce a continuous time specification, and allow learning of the latent space dimension. Bhattacharya & Dunson (2011) proposed a useful approach for Bayesian learning of the number of latent factors in a model for a single large covariance matrix, and we extend their approach from independent Gaussian latent factors to Gaussian process latent factors. In particular, we let

independently for all $i=1,...,V$ and $h=1,\dots,H$, with $c_{X}$ a squared exponential correlation function $c_{X}(t,t^{\prime})=\exp(-\kappa_{X}||t-t^{\prime}||^{2}_{2})$, which allows for continuous time analysis and unequal spacing, and $\tau_{h}^{-1}$ a shrinkage parameter defined as

Note that if $a_{2}>1$ the expected value for $\vartheta_{k}$ is greater than $1$. As a result, as $h$ goes to infinity, $\tau_{h}$ tends to infinity, shrinking $x_{ih}(\cdot)$, for every $i=1,\ldots,V$ towards zero. This leads to a flexible prior for $x_{ih}(\cdot)$ with a local shrinkage parameter $\tau_{h}^{-1}$ that favors many stochastic processes of latent coordinates being close to as $h$ increases. To conclude prior specification we choose

with $c_{\mu}(t,t^{\prime})=\exp(-\kappa_{\mu}||t-t^{\prime}||^{2}_{2})$.

Before proceeding with posterior computation, we focus on the support of the induced prior $\Pi_{\pi}$ based on priors $\Pi_{X}$ and $\Pi_{\mu}$. Specifically we are interested in proving whether the prior can generate a time-varying symmetric probability matrix that is arbitrarily close to any function $\{\pi(t),t\in\mathcal{T}\}$. Intuitively, large support on continuous symmetric similarity matrix functions $\{S(t),t\in\mathcal{T}\}$ relies on the continuity of the Gaussian process coordinate functions. Since for each fixed $t=t_{0}$, $x_{ih}(t_{0})$ are independently Gaussian distributed, $X(t_{0})X(t_{0})^{\prime}$ is distributed according to a sum of independent Wishart random variables. Combining the large support of the Wishart distribution with the one of the Gaussian for the baseline $\mu(t_{0})$, provides large support for the induced prior $\Pi_{S}$. Since $\pi(t)$ is obtained via a one to one continuous increasing function of $S(t)$, we will map non-null probability subsets of the space of $S(t)$ into non-null probability subsets of the space of $\pi(t)$, providing the desired large support for the induced prior $\Pi_{\pi}$. Theorem 3 states the large support property for $\Pi_{S}$, while Corollary 4 provides the same property for $\Pi_{\pi}$ by combining results in the previous Theorem with the fact that $\pi(t_{0})$ is defined as a monotonic increasing continuous function of $S(t_{0})$. Proof of Theorem 3 is provided in Appendix.

Let $\Pi_{S}$ denote the induced prior on $\{S(t),t\in\mathcal{T}\}$ based on the specified prior $\Pi_{X}\otimes\Pi_{\mu}$ on $\mathcal{X}_{X}\otimes\mathcal{X}_{\mu}$. Assuming $\mathcal{T}$ compact, for all continuous $S^{*}(t)$ and for all $\epsilon>0$

Let $\Pi_{\pi}$ denote the induced prior on $\{\pi(t),t\in\mathcal{T}\}$ based on the specified prior $\Pi_{X}\otimes\Pi_{\mu}$ on $\mathcal{X}_{X}\otimes\mathcal{X}_{\mu}$. Assuming $\mathcal{T}$ compact, for all continuous $\pi^{*}(t)$ and for all $\delta>0$

Proof. Since the elements of $\pi(t)$ are defined as a one to one continuous mapping of the elements of $S(t)$ through the function $g(\cdot)$, by definition of continuity we have that for every $\delta>0$ there exists an $\epsilon>0$ such that

for all $S(t)$ such that $\sup_{t\in\mathcal{T}}||S(t)-S^{*}(t)||_{2}<\epsilon$, where $g\left(S(t)\right)$ means that the function $g(\cdot)$ is applied to every element of $S(t)$. Finally, since by Theorem 3 the event $\sup_{t\in\mathcal{T}}||S(t)-S^{*}(t)||_{2}<\epsilon$ has non-null probability, it follows that the same holds for the event $\sup_{t\in\mathcal{T}}||\pi(t)-\pi^{*}(t)||_{2}<\delta$, completing the proof.

Posterior computation

Posterior computation is performed adapting a recently proposed data-augmentation scheme based on a new class of Pólya-Gamma distributions; for a detailed description see Polson et al. (2013). The approach provides a strategy for fully Bayesian inference in models with binomial likelihoods, which bypasses the need for analytic approximations, while allowing us to exploit conjugacy for block updating.

The main result is that binomial likelihoods parameterized by log-odds can be represented as a mixture of Gaussians with respect to Pólya-Gamma distributions. Specifically

where $z=a-b/2$ and $\omega\sim\mbox{PG}(b,0)$, with $\mbox{PG}(b,c)$ denoting the Pólya-Gamma random variable with parameters $c\in\Re$ and $b>0$. When $\psi=x^{\prime}\beta$ is a linear predictor, and a Gaussian prior is considered for $\beta$, full conditional conjugacy is ensured for Bayesian inference on the coefficients. Moreover the implied conditional distribution for $\omega$, given $\psi$, is again Pólya-Gamma, providing a simple Gibbs sampler alternating between two main steps. Specifically, letting $y_{i}$ be the number of successes and $x_{i}=[x_{i1},...,x_{ip}]^{\prime}$ the vector of regressors for every observation $i=1,...,N$, and assuming a Bayesian logistic regression setting where $y_{i}\sim\mbox{Bern}(1/[1+e^{-\psi_{i}}])$, $\psi_{i}=x_{i}^{\prime}\beta$ and $\beta$ having Gaussian prior $\beta\sim\mbox{N}_{p}(b,B)$, the Gibbs alternates between

where $\Sigma_{\beta}=(X^{\prime}\Omega X+B^{-1})^{-1}$ and $\mu_{\beta}=\Sigma_{\beta}(X^{\prime}z+B^{-1}b)$; with $z=[y_{1}-1/2,....,y_{N}-1/2]^{\prime}$ and $\Omega$ is the diagonal matrix with $\omega_{i}$’s entries.

Recalling model (1), with probabilities defined as in (2) and latent similarities from (3), for $i=2,...,V$, $j=1,...,i-1$ and $t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}$, and taking a fixed truncation level $H^{*}$ for the number of latent coordinates, the Gibbs sampler for our model, is:

Update each augmented data $\omega_{ij,t}$ from the full conditional Pólya-Gamma posterior:

for every $i=2,...,V$, $j=1,...,i-1$ and $t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}$.

Given $\{y_{ij,t}\}$, $X(t)$ and $\{\omega_{ij,t}\}$, the Pólya-Gamma data augmentation scheme ensures full conditional Gaussian posterior for $\mu(t)$ with $t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}$, of the form

With $\Sigma_{\mu}=\left[\mbox{diag}\left(\sum_{i=2}^{V}\sum_{j=1}^{i-1}\omega_{ij,t_{1}},\dots,\sum_{i=2}^{V}\sum_{j=1}^{i-1}\omega_{ij,t_{T}}\right)+K_{\mu}^{-1}\right]^{-1}$, and $K_{\mu}$ the Gaussian process covariance matrix with $[K_{\mu}]_{ij}=\exp(-\kappa_{\mu}||t_{i}-t_{j}||^{2}_{2})$.

Update the time-varying latent coordinate vector $\{x_{v}(t)=[x_{v1}(t),...,x_{vH^{*}}(t)]^{\prime}\}_{t=t_{1}}^{t_{T}}$ for every unit $v=1,...,V$ from its conditional posterior. Specifically, conditionally on $X^{(-v)}=\{x_{j}(t):j\neq v,t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}\}$, $\mu=[\mu(t_{1}),\dots,\mu(t_{T})]^{\prime}$, $\{y_{ij,t}\}$, $\{\omega_{ij,t}\}$, $\{\tau_{h}\}$ and defining $y_{ij}=[y_{ij,t_{1}},\dots,y_{ij,t_{T}}]^{\prime}$ and $\pi_{ij}=[\pi_{ij,t_{1}},\dots,\pi_{ij,t_{T}}]^{\prime}$, let $Y^{(v)}$ be the vector obtained by stacking sub-vectors $y_{ij}$ for all the couples $(i,j)$ such that $i=v$ or $j=v$, with $i>j$; and $\pi^{(v)}$ the corresponding vector of probabilities, then

for all the probabilities $\pi_{ij}(t)$ such that $i=v$ or $j=v$, with $i>j$ and $t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}$. Model (6) is a proper logistic regression with linear predictor, therefore, according to our Pólya-Gamma sampling scheme, we update the vector of time-varying coordinates $\{x_{v}(t)=[x_{v1}(t),...,x_{vH^{*}}(t)]^{\prime}\}_{t=t_{1}}^{t_{T}}$, represented by $\beta_{x_{v}(t)}$ by sampling from:

and $\Omega_{x_{v}(t)}$ is the diagonal matrix with the corresponding Pólya-Gamma augmented data.

Conditioned on $X(t)$ and $\{\tau_{h}\}$, sample the global shrinkage hyperparameters from

Where $\tau_{l}^{(-h)}=\prod_{t=1,t\neq h}^{l}\vartheta_{t}$ for $h=1,...,H^{*}$ and $x_{il}=[x_{il}(t_{1}),\dots,x_{il}(t_{T})]^{\prime}$.

We can easily handle missing values by adding a further step imputing the unobserved binary similarities from their conditional distribution given the current state of the chain. Specifically:

5. Given $X(t)$ and $\mu(t)$ sample each missing value from its conditional distribution

Step 5 provides also a strategy for predicting new outcomes. Specifically, if we are interested in making inference on future $\pi(t_{T+1})$ with $t_{T+1}>t_{T}$ given the observed similarity matrices $Y_{t}$, $t\in\mathcal{T}_{0}=\{t_{1},...,t_{T}\}$, then we can simply perform the previous posterior computations adding to the observed dataset $\{Y_{t}\}_{t\in\mathcal{T}_{0}}$ a new matrix $Y_{t_{T+1}}$ of missing values and make inference on the predictive posterior distribution using the samples of the Markov chain for $\pi(t_{T+1})$.

Simulation Study

We provide a simulation study with the aim to evaluate the performance of the proposed model in analyzing a dataset constructed to mimic also a possible generating process in the finance application. The focus is on the ability to correctly reconstruct the true underlying processes, and also on the performance with respect to out of sample predictions. We also provide a comparison between our proposed approach and the estimated probability process for each time-varying binary outcome when using only temporal information without exploring matrix structure, showing graphically the sub-optimality of the latter in terms of efficiency and bias.

We generate a set of $15\times 15$ time varying $Y_{t}$ matrices with $t$ in the discrete set $\mathcal{T}_{0}=\{1,2,\dots,40\}$. Each $y_{ij,t}$ is simulated according to (1) with probabilities obtained from (2) and (3), generating $\{\mu(t)\}_{t=1}^{40}$ from a $\mbox{GP}(0,c_{\mu})$ with length scale $\kappa_{\mu}=0.01$ and choosing $2$ time-varying latent coordinates $\{x_{i1}(t)\}_{t=1}^{40}$, $\{x_{i2}(t)\}_{t=1}^{40}$ from Gaussian processes with length scale $\kappa_{x}=0.01$, independently for each unit $i=1,...,15$. To evaluate the out of sample predictive performance we take $Y_{40}$ to be a matrix of missing values, and assume similarities between units $10$ and $11$ and all the others, missing at times $t=20,...,25$ to assess the behavior with respect to missing data. For inference we choose a truncation level $H^{*}=10$, length scales $\kappa_{\mu}=\kappa_{x}=0.05$ and set $a_{1}=a_{2}=2$ for the shrinkage parameters. We ran $5{,}000$ Gibbs iterations which proved to be enough for reaching convergence and discarded the first $1{,}000$. Mixing was assessed by analyzing the effective sample sizes of the MCMC chains for the quantities of interest (i.e. $\pi_{ij}(t)$, for $i=2,...,V$, $j=1,...,i-1$ and $t\in\mathcal{T}_{0}$) after burn-in. We found most of these values concentrating around $\approx 1{,}700$ effective samples on a total of $4{,}000$, providing a good mixing result.

The comparison in Figure 2 between true probability matrices and their corresponding posterior mean for some selected time $t$, highlights the good performance of our approach in correctly estimating the true latent process and making predictions. The latter can be noticed by comparing true and estimated probability matrices at $t=40$, recalling that in our simulation we assumed $Y_{40}$ having missing entries and we were interested in analyzing the predictive performance of our model with respect to $\pi(40)$. Similar results are provided by the plot of true $\pi_{ij}(t)$ against the corresponding estimates $\hat{\pi}_{ij}(t)$ and by the ROC curve in Figure 3 having an area underneath of $0.87$.

Figure 4 shows a graphical comparison between the performance of our model with respect to $\mu(t)$ and some selected probability trajectories $\pi_{ij}(t)$ (top), and the inferential results when the mean process and probability process $\pi_{ij}(t)$ are estimated with the same setting of our model but using only the time series of the corresponding $y_{ij,t}$ without borrowing information across the network (bottom). The sub-optimality of the independent approach is apparent in terms of both bias (over-smoothed trajectories) and variance (larger hpd intervals). When network structure is taken into account, the model provides accurate estimates, with posterior distributions rapidly concentrating around the true corresponding processes, while accurately selecting the dimension of the latent space. In particular, we find that the estimated $\hat{\tau}_{h}^{-1}$ values start at 0.8 and 0.7 for $h=1$ and $2$, respectively, but then drop to small values for the later factors. This implies that these later factor trajectories are quite flat and have limited influence. Borrowing information across the network over time has the additional advantage of reducing hyperparameter sensitivity, in particular with respect to the length scale in GP prior. We obtain, in fact, similar results when instead letting $\kappa_{\mu}=\kappa_{x}=0.03$, $\kappa_{\mu}=\kappa_{x}=0.1$ and $\kappa_{\mu}=\kappa_{x}=0.5$ in sensitivity analyses.

Application to co-movements among National Stock Market Indices

National Stock Indices represent technical tools constructed by a synthesis of numerous data on the evolution of the various stocks, and represent important indicators of the financial condition in a given country. Modeling co-variations among these quantities, and in general among assets, represents a fundamental issue in many financial applications, such as the Arbitrage Pricing Theory (APT) of Ross (1976) and the Capital Asset Pricing Model (CAPM) developed by Sharpe (1964), and the correlations or covariances among asset’s returns are the typical measures of co-movements employed in this framework.

A rich literature is available in modeling dynamic covariance or correlation matrices, covering multivariate generalizations of ARCH and GARCH models (see e.g. Tsay, 2005, Engle, 2002, Alexander, 2001, Bollerslev et al., 1988), Stochastic volatility models (Harvey et al., 1994) and recent Bayesian extensions (see e.g. Wilson & Ghahramani, 2010, Nakajima & West, 2012, Durante et al., 2013). In this application, we instead provide a different and not fully explored measure of co-movement exploiting the network structure among financial indices and giving exactly the probability that such event happens at a given time. This is accomplished by applying our model to the time-varying $Y_{t}$ matrices having entries $y_{ij,t}=y_{ji,t}=1$ if index $i$ and index $j$ co-move at time $t$ (indices are similar), and $y_{ij,t}=y_{ji,t}=0$ if opposite increments are recorded (indices are dissimilar).

We constructed $Y_{t}$ using the quarterly log-returns of the $23$ main National Stock Market Indices ($V=23$) from 2004 to 2013 ($T=39$, with the last empty matrix $Y_{39}$ used for prediction), available at http://finance.yahoo.com/ and applied model (1), with probabilities specified as in (2) and latent similarity measures obtained via the projection approach defined in (3). For posterior computation we run $5{,}000$ Gibbs iterations with a burn-in of $1{,}000$, setting a truncation level $H^{*}=15$, length scales $\kappa_{\mu}=0.03$, $\kappa_{x}=0.01$ and $a_{1}=a_{2}=2$. Similarly to the simulation study, most of the chains have effective sample sizes around $1{,}600$ on a total of $4{,}000$ after burn-in, showing good mixing. We find that the first two latent factors are the most informative, with the remaining $13$ latent processes being concentrated near zero. A similar result was obtained in the seminal work of Fama & French (1993), providing three main common risk factors in the returns of stocks.

The estimated trajectory of the baseline process $\mu(t)$ together with the point-wise 0.95 hpd intervals in Figure 5, provide important insights on the overall financial market behavior, in agreement with other theories on financial crises (see, e.g., Baig and Goldfaijn, 1999, and Claessens & Forbes, 2009) and recent applications (Durante et al., 2013, Kastner et al., 2013). Increasing and persistent level of the baseline process, inducing higher probability of co-movements, are recorded during the growth and burst of USA housing bubble and the initial turmoils before the 2008 global financial crisis (A). This result provides an empirical proof in favor of the increasing inter-connection among financial markets due to the proliferation of risky loans between 2004 and 2007, and the growing demand by foreign countries for financial assets built from the real estate market, such as residential mortgage-backed securities (RMBS) and collateralized debt obligations (CDO). As expected the global financial crisis between late-2008 and end-2009 (B), and the following, Greek debt crisis together with the worsening of European sovereign-debt crisis (C), are manifested through a further increase of the co-movement probabilities, highlighting a clear financial contagion effect.

Figure 6 shows the estimated (blue lines) and predicted (red lines) co-movement probability trajectories among USA and some selected European countries, pointing out the good performance of the proposed model in adaptively learning the data structure, confirmed also by a ROC curve having an area underneath of $0.79$. It is worth noticing that the local adaptivity of the estimated trajectories is not due to an over-parameterization of the model since the shrinkage prior on $\tau_{h}$ and the choice of small length scales in the GP covariance functions, imply smooth trajectories and a parsimonious model formulation. Thus adaptivity is provided by the information borrowed in the financial network for each time $t$. Co-movement probabilities among USA and Greece register a sharp drop in correspondence of the Greek debt crisis, differently from what happens with other European countries such as Germany and France, which instead evolve on similar patterns. We found this result reasonable in providing an empirical proof on the attempt to reduce the inter-connection with a country in crisis.

Finally, Figure 7 provides interesting insights on the financial network structure among the countries under investigation. Specifically we represent three different weighted networks, with weights given by the average estimated co-movement probability over all the time window considered (a), the estimated probability averaged over the period of the global financial crisis (b), and the Greek debt crisis (c). A reasonable global network structure with countries having similar financial economies most closely related among each other is provided in plot (a). As expected Japan appears to be closer to Western economies than Asian financial markets, while China has lower inter-connections with other countries. Stronger networks are estimated for European markets and Asian Tigers. International financial contagion effect is highlighted through strong inter-connections among all financial markets during the 2008 global financial crisis (b), with a still evident clustering effect, and Greece already showing a slightly different behavior. Finally, when the network during the Greek debt crisis is analyzed, we register evident low connections among Greece and almost all the other financial markets considered, and interestingly learn a strong network between Greece, Spain and Italy, representing the countries most affected by the European sovereign-debt crisis.

Discussion

We proposed a Bayesian nonparametric dynamic model for binary similarity matrices, borrowing information across time and the network structure of the data under investigation and allowing for dimensionality reduction. The model has been constructed using latent similarity measures defined by the dot product of latent coordinate vectors, with entries evolving in continuous time via Gaussian process priors. The shrinkage hyperprior allows us to automatically learn the dimension of the latent space and ensures a parsimonious definition of the model, with the risk of over-parameterization due to a higher number of latent features avoided. The Pólya-Gamma data augmentation strategy allows us to define a simple and efficient Gibbs sampler for posterior computations based on full conditional conjugate posterior distributions, which is promising in terms of scaling to moderately large $V$, and easily handling missing values as well as forecasting problems. Scalability to large $T$ could be, instead, improved via stochastic differential equations models approximating the GP prior on the latent coordinate processes (Zhu and Dunson, 2012). We provided also theoretical results on the flexibility of the model, illustrated its performance via a simulation study and obtained interesting insights on the network among financial markets during the recent crisis, by applying the model to time-varying co-movement data.

Our model has a broad range of applicability, with dynamic social network analysis and time-varying binary evaluations among units providing two natural fields of application. Further directions of research could be devoted to the definition of similar models for discrete valued dynamic matrices, which could provide useful tools for analyzing edge valued dynamic social networks or datasets with comparison among units expressed on a Likert scale.

References

Appendix

Proof of Theorem 3: Since $\mathcal{T}$ is compact, for every $\epsilon_{0}>0$ there exists an open covering of $\epsilon_{0}$-balls $B_{\epsilon_{0}}(t_{0}):\{t:||t-t_{0}||_{2}<\epsilon_{0}\}$ with a finite subcover such that $\mathcal{T}\subset\cup_{t_{0}\in\mathcal{T}_{0}}B_{\epsilon_{0}}(t_{0})$, where $|\mathcal{T}_{0}|=T$. Then:

Define $Z(t_{0})=\sup_{t\in B_{\epsilon_{0}}(t_{0})}||S(t)-S^{*}(t)||_{2}$. Since

we only need to look at each $\epsilon_{0}$-ball independently as follow:

Where the first inequality comes from repeated uses of triangle inequality, and the second follows from the fact that each of these terms is an independent event. We evaluate each of these terms in turn.

Based on the continuity of $S^{*}(\cdot)$, for all $\epsilon/3>0$, there exists an $\epsilon_{0,1}>0$ such that:

Therefore, $\Pi_{S}\left(\sup_{t\in B_{\epsilon_{0,1}}(t_{0})}||S^{*}(t_{0})-S^{*}(t)||_{2}<\frac{\epsilon}{3}\right)=1$.

Given the GP prior on the elements of $X(\cdot)$ and letting $x_{ih}(t)=[X(t)]_{ih}$, the equation

represents a finite sum over pairwise products of almost surely continuous functions (recalling GP assumption on the elements $x_{ih}$) and thus result in a matrix $X(t)X(t)^{\prime}$ with elements almost surely continuous on $\mathcal{T}$. Therefore $S(t)=\mu(t)\times 1_{V}1_{V}^{\prime}+X(t)X(t)^{\prime}$ is almost surely continuous on $\mathcal{T}$ since the baseline $\mu(\cdot)$ is itself almost surely continuous given the GP prior assumption. Therefore, similarly as before, for all $\epsilon/3>0$, there exists and $\epsilon_{0,2}>0$ such that:

Where $\{X(t_{0})^{*},\mu^{*}(t_{0})\}$ is any element of $\mathcal{X}_{X}\otimes\mathcal{X}_{\mu}$ such that $S^{*}(t_{0})=\mu^{*}(t_{0})\times 1_{V}1_{V}^{\prime}+X(t_{0})^{*}X(t_{0})^{*^{\prime}}$. Such a factorization exists by Corollary 2. Thus, using triangle inequality, we can bound this probability by:

Based on the support of the Gaussian prior,

For studying the first term of the previous decomposition note that:

where $x_{h}(t_{0})=[x_{1h}(t_{0}),\dots x_{Vh}(t_{0})]^{\prime}$ is distributed, according to our prior specification, as $\mbox{N}_{V}(0,\tau_{h}^{-1}I_{V})$, implying that $x_{h}(t_{0})x_{h}(t_{0})^{\prime}|\tau_{h}\sim\mbox{W}_{V}(1,\tau_{h}^{-1}I_{V})$ independently for all $h=1,...,H$, where $\mbox{W}_{V}(\cdot,\cdot)$ denotes the Wishart random variable. Using the triangle inequality we obtain:

Since $x_{h}(t_{0})^{*}x_{h}(t_{0})^{*^{\prime}}$ is an arbitrary rank-1 symmetric matrix in $\Re^{V\times V}$, and based on the support of the Wishart distribution:

Thus $\Pi_{S}\left(||X(t_{0})X(t_{0})^{\prime}-X(t_{0})^{*}X(t_{0})^{*^{\prime}}||_{2}<\frac{\epsilon}{6}\right)>0$ and combining it with the large support property previously proved for the prior on the baseline $\mu(\cdot)$, we have:

For every $S^{*}(\cdot)$ and $\epsilon>0$, let $\epsilon_{0}=\min(\epsilon_{0,1},\epsilon_{0,2})$, with $\epsilon_{0,1}$ and $\epsilon_{0,2}$ defined as above. Then, combining the positivity results of each of the three terms in 7 completes the proof.