A Dynamic Stochastic Block Model for Multidimensional Networks

A Dynamic Stochastic Block Model for Multidimensional Networks

Introduction

Network models are a tool for studying real-world complex systems and have become a convenient framework for describing social and economic interactions de2017econometrics. In terms of higher--order network properties, such as clustering, several stochastic models have been proposed apart from exponential random graphs, e.g., latent space models and Stochastic Block Models (SBMs) kim2018review,hoff2018. However, in systems with multiple types of interactions between nodes, a single graph does not fully describe the connectivity structure. Thus, graphs with multiple types of edges (layers) have been introduced. See Kiv14 for an introduction to multidimensional networks. To the best of our knowledge, relatively few works deal with dynamic clustering models for multidimensional graphs lee2019review,paul2020random,lei2022bias. The objective of the present work is to extend the Dynamic SBM (DSBM) to accommodate multiple edge types. \textcolor{black}{Our approach can be easily extended to a general multi--layer framework where edges between layers are allowed.}

In SBMs, edge clustering is driven by a probabilistic classification of the nodes into different communities. \textcolor{black}{In this paper, the terms group, cluster, block, and community are used interchangeably to describe a set of nodes sharing similar connectivity characteristics. These characteristics may include edge probability, strength, and variance of the strength.} An alternative dynamic specification involves using Hidden Markov Chains (HMCs) to capture temporal changes in the node's membership. See fruhwirth2006finite for an introduction to Hidden Markov models. yang2011detecting work with (un)directed and unweighted dynamic networks and matias2017statistical generalize their model to include weighted networks, discussing the identification issues that arise when block--dependent connectivity parameters are time--varying. Other extensions of SBMs address mixed--membership and further edge formation heterogeneity airoldi2008mixed,zhao2012consistency.

In these previous works, the DSBM has accounted for one type of edge. However, in social and economic relationships, there is usually more than one type of tie, such as multiple goods traded between regions or several assets exchanged by financial institutions. Therefore, we propose a DSBM for multidimensional networks (DSBMM) where each type of relationship is represented by a different layer, with no inter--layer edges and a time--invariant and layer--invariant node set \textcolor{black}{and with a fixed number of communities, potentially different between layers}.

To the best of our knowledge, extensions of the SBM in a multidimensional setting are static and correlational\textcolor{black}{---undirected dependence relationships}. jovanovski2019bayesian identify communities per layer (local clustering), and a consensus of all layers (global clustering) with unweighted edges. stanley2016clustering in a static framework, assume layer clustering, in such a way that the layers in the same group can share a common SBM. Other studies with no interest in community detection use multivariate distributions to make inference on the dependence between layers, but this restricts all layers to be of the same type, that is, either weighted or unweighted and directed or undirected salter2017latent. \textcolor{black}{Alternative models to the SBM family have been proposed to capture clustering in a dynamic setting, such as the latent space model proposed by durante2017bayesian, where layer--specific and shared latent positions evolve smoothly following Gaussian processes. Although this model is flexible enough to adapt to different network structures, connectivity dynamics, and to exploit information across layers, there are significant differences with the DSBMM. First, the community structure provides a clear interpretation in terms of node classification. In contrast, the latent position coordinates may require additional knowledge of the field of application to address identifiability issues, as in factor models. Second, although the shared latent positions allow for layer dependence, it does not provide information on the direction of the relationship between each pair of layers, only an undirected global dependence measure between all layers. Moreover, this dynamic latent position model assumes that all layers are unweighted and undirected, whereas the DSBMM encompasses layers with different relationship types simultaneously.} Although these alternative multidimensional models can provide a measure of association or clustering overlap between layers, it is not possible to infer if the nature of the relationship is unidirectional or bidirectional.

In this paper, we introduce a concept of nonlinear Granger block--causality (\textcolor{black}{NGB causality}) that identifies the directed dependence between layers' community structure. For instance, Layer $\ell$ may be \textcolor{black}{NGB} causing Layer $\mathfrak{m}$, but not vice--versa. Specifically, in the DSBMM, the nodes' membership dynamics in each layer is influenced by their respective membership in the other layers. By using a saturated multinomial specification to model the transition matrix of each layer, it is possible to test for Granger causality between layers. The approach also allows for layers of different types, i.e. (un)directed or (un)weighted, and a set of layer--specific covariates to control for the observed node or dyad heterogeneity.

We propose a Bayesian inference, based on a Multi--Laplacian prior, to induce a group LASSO penalty and address the overparameterization issue raman2009bayesian. A Pólya Gamma representation is introduced to make the multinomial model tractable polson2013bayesian. Our simulation results show that the Multi--Laplacian prior performs better than a Normal prior in retrieving the underlying causal structure of the DGP (\textcolor{black}{NGB causality}). This difference in performance is due to the twofold shrinkage of the Multi--Laplacian prior. In a Normal prior setting, information on the parameter grouping is excluded, thus the correlation between \textcolor{black}{lagged HMC variables} is not exploited, whereas in a Multi--Laplacian setting, different shrinkage effects, within groups of parameters and between groups, are allowed. A comparison between DSBMM and a benchmark model, that is a Bayesian Vector Autoregression (BVAR), shows that the latter is not able to detect the \textcolor{black}{NGB causality} under the great majority of scenarios.

We apply our DSBMM model to trade flows and free trade agreements (FTA) and contribute to the debate on the effectiveness of FTA as a policy instrument for global trade integration baier2007free,baier2019widely. The wide range of empirical results regarding the effects of FTAs on trade, combined with the complexity of overlapping FTAs between countries, suggests that considering the network structure can provide a complementary view of the FTA's influence on the topological properties of trade flows. Specifically, we test if FTA clustering helps predict the international trade community membership after controlling for country and dyad observed heterogeneity. Moreover, our model--based approach for community detection accounts for uncertainty in the parameters and community membership. Community detection in international trade flows has been analyzed mainly using modularity algorithms, which are descriptive and do not provide probabilistic statements. For example, bartesaghi2020communicability focuses on aggregate trade flows to infer groups of countries with dense connectivity, and barigozzi2011identifying uses commodity-specific layers to identify communities of countries per product. Furthermore, by allowing covariates to affect trade flows, our DSBMM integrates community detection into Gravity equation models, which are extensively used in international trade piermartini2016estimating.

The paper is organized as follows. Section (ref) introduces the DSBMM and presents some numerical illustrations of the model properties. Section (ref) addresses the inference procedure, including the choice of the prior distribution, the details of the Gibbs sampler, and a comparison with panel BVAR models {\color{black}and other dynamic network models} using synthetic data. Section (ref) provides an application to the Trade Flows--FTAs multidimensional network. Section (ref) summarizes the main conclusions.

A Dynamic SBM for multidimensional networks

The community detection task within a dynamic framework is defined in this section, and the multidimensional dependence is presented in detail, accompanied by examples.

Multidimensional DSBM

A graph can be defined as the ordered triplet $\mathcal{G}=\{\mathcal{V},\mathcal{E},Y\}$, where the node set $\mathcal{V}=\{1,\cdots,N\}$ is the same across the layer set $\mathcal{L}=\{1,\ldots,L\}$ and the edge set collection $\mathcal{E}=\{\mathcal{E}^{(1)},\ldots,\mathcal{E}^{(L)}\}$, $\mathcal{E}^{(\ell)}\subseteq \mathcal{V}\times \mathcal{V}$ allows only for intra--layer edges. The set $Y=\{Y^{(1)},\cdots, Y^{(L)}\}$ includes the layer--specific adjacency matrices with $Y^{(\ell)}\in \mathbb{R}^{N}\times\mathbb{R}^{N}$ where the $(i,j)$--th element of $Y^{(\ell)}$ is $Y^{(\ell)}_{ij}=0$ if $(i,j)\not\in \mathcal{E}^{(\ell)}$ and $Y^{(\ell)}_{ij}=a\in \mathbb{R}\backslash\{0\}$ if $(i,j)\in \mathcal{E}^{(\ell)}$. The multidimensional may jointly consider (un)directed and (un)weighted layers, resulting in (symmetric) asymmetric and (binary) real-valued adjacency matrices.

Considering that higher-order network structures, such as clubs or closed communities, are common features in real-world graphs csermely2013structure, an SBM can be used in each layer to model network topology through node grouping. We define the set $\mathfrak{V}=\{\mathfrak{V}^{(1)},\dots,\mathfrak{V}^{(L)}\}$, where each element is a layer--specific partition of the node set into $Q^{(\ell)}$ subsets, that is, $\mathfrak{V}^{(\ell)}=\{\mathcal{V}^{(\ell)}_{1},\ldots,\mathcal{V}^{(\ell)}_{Q^{(\ell)}}\}$ and each $\mathcal{V}^{(\ell)}_{q}\subseteq \mathcal{V}$ is called block or community in layer $\ell$, satisfying two properties: $\mathcal{V}^{(\ell)}_{q}\cap \mathcal{V}^{(\ell)}_{r}=\emptyset$, $q\neq r$, and $\cup_{q\in \mathcal{Q}^{(\ell)}}\mathcal{V}^{(l)}_{q}=\mathcal{V},\ \mathcal{Q}^{(\ell)}=\{1,\dots,Q^{(\ell)}\}$. \textcolor{black}{The number of blocks $Q^{(\ell)}$ is fixed.}

It is assumed that the multidimensional graph evolves over time, $\mathcal{G}_{1:T}=\{\mathcal{G}_{t}\}_{t\in \mathcal{T}}$, with $\mathcal{G}_{t}=\{\mathcal{V},\mathcal{E}_t,Y_t\}$, $\mathcal{T}=\{1,\ldots,T\}$, and a latent sequence of partitions $\mathfrak{V}_{1:T}=\{\mathfrak{V}_t\}_{t\in \mathcal{T}}$ drives its topology. For a layer $\ell$, the block membership of a node is indicated by the latent variable $Z^{(\ell)}_{it}=q\in \mathcal{Q}^{(\ell)}$, for $t\in \mathcal{T}$ \textcolor{black}{and $i\in \mathcal{V}$}. Notice that the partition sequence and the memberships are intrinsically related $\mathcal{V}_{qt}^{(\ell)}=\{i\in\mathcal{V}|Z^{(\ell)}_{it}=q\}$ \textcolor{black}{for $q\in \mathcal{Q}^{(\ell)}$, $t\in \mathcal{T}$ and $\ell \in \mathcal{L}$}. Then, the membership dynamics $Z^{(\ell)}_{{\color{black}i\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}}}=\{Z^{(\ell)}_{it}\}_{t\in \mathcal{T}}$, follows a hidden Markov--chain process \textcolor{black}{for $i \in \mathcal{V}$}. Specifically, for the initial membership $Z^{(\ell)}_{i1}$, the probabilities of belonging to a block are collected into the vector $\alpha^{(\ell)}=(\alpha^{(\ell)}_1,\dots,\alpha^{(\ell)}_{Q^{(\ell)}})'$, where $\sum_{q\in \mathcal{Q}^{(\ell)}} \alpha^{(\ell)}_q =1$. In other words, $Z_{i1}^{(\ell)}\sim\operatorname{Ca}(\alpha^{(\ell)})$, where $\operatorname{Ca}(p)$ denotes a categorical (or multinoulli) distribution with probability vector $p$. The membership changes for $t\in \mathcal{T}\backslash\{1\}$ are determined by the transition matrix $P_{it}^{(\ell)}$ \textcolor{black}{for $i \in \mathcal{V} $}, where \textcolor{black}{$P_{i t}^{(\ell)}=\left(P_{i t, r q}^{(\ell)}\right)_{r, q \in \mathcal{Q}^{(\ell)}}$,} $P^{(\ell)}_{it,qr}\in (0,1)$ for $q,r \in \mathcal{Q}^{(\ell)}$ and $\sum_{r\in\mathcal{Q}^{(\ell)}} P^{(\ell)}_{it,qr}=1$. The transition probabilities vary across nodes and time. Thus, $Z_{it}^{(\ell)}|Z_{it-1}^{(1:L)},\sim\operatorname{Ca}(P_{it,Z^{(\ell)}_{it-1}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}}^{(\ell)})$, where $P_{it,Z^{(\ell)}_{it-1}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}}^{(\ell)}=(P_{it,Z^{(\ell)}_{it-1}1}^{(\ell)},\ldots,P_{it,Z^{(\ell)}_{it-1}Q^{(\ell)}}^{(\ell)})'$.

The contemporaneous network $Y^{(\ell)}_t=\{Y^{(\ell)}_{ijt}\}_{i,j\in \mathcal{V}}$ only depends on $Z^{(\ell)}_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}}=\{Z^{(\ell)}_{it}\}_{i\in\mathcal{V}}$ and layer--specific covariates $X^{(\ell)}_{1:T}=\{X_{ijt}^{(\ell)}\}_{i,j\in \mathcal{V},\ t\in \mathcal{T}}$, i.e. given $Z^{(\ell)}_{1:T}\textcolor{black}{=\{Z^{(\ell)}_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}}\}_{t\in \mathcal{T}}}$ and $X^{(\ell)}_{1:T}$, $Y^{(\ell)}_{1:T}=\textcolor{black}{\{Y^{(\ell)}_t\}_{t\in \mathcal{T}}}$ components are independent \textcolor{black}{across time and edges}. Each entry of the $\ell$--th adjacency matrix follows the mixed distribution,

{\color{black} where $\vartheta^{(\ell)}=\{\nu^{(\ell)}, \theta^{(\ell)}\}$ are the set of connectivity parameters with the $Q^{(\ell)}$ square matrix $\nu^{(\ell)}=(\nu_{qr}^{(\ell)})_{q,r\in \mathcal{Q}^{(\ell)}}$ and the three--dimensional tensor collecting $(S+1)$ square matrices of dimension $Q^{(\ell)}$ $\theta^{(\ell)}=(\theta_{sqr}^{(\ell)})_{s \in \{0,1,\ldots,S\}\ , q,r\in \mathcal{Q}^{(\ell)}}$. The function $\delta(\cdot)$ denotes the Dirac function at zero, $\nu_{qr}^{(\ell)}$ the probability of having an active edge between two nodes, and $f^{(\ell)}(y|\theta^{(\ell)}_{qr}, X^{(\ell)}_{ijt})$ is a layer--specific probability (mass) density function with parameters $\theta^{(\ell)}_{qr}$. Notice that the DSBM can be interpreted as a dimensionality reduction of the tensor $Y_{1:T}^{(\ell)}$ with dimensions $N^2T$ that is summarized by the set of features in the set $\vartheta^{(\ell)}$ with dimensions $(S+2)Q^{(\ell)2}$. The dynamic network and the set of features are linked through the block memberships by the following relationships $$\nu_{Z^{(\ell)}_{it}Z^{(\ell)}_{jt}}^{(\ell)}=\tilde{Z}^{(\ell)\prime}_{it}\nu^{(\ell)}\tilde{Z}^{(\ell)}_{jt}, \quad \quad \theta^{(\ell)}_{Z^{(\ell)}_{it}Z^{(\ell)}_{jt}}=(\tilde{Z}^{(\ell)\prime}_{it}\theta_{0}^{(\ell)}\tilde{Z}^{(\ell)}_{jt},\ldots,\tilde{Z}^{(\ell)\prime}_{it}\theta_{S}^{(\ell)}\tilde{Z}^{(\ell)}_{jt})',$$ where $\theta_{s}^{(\ell)}=(\theta_{sqr}^{(\ell)})_{q,r\in \mathcal{Q}}$ is a slice of the tensor $\theta^{(\ell)}$ for $s\in \{0,\ldots,S\}$, and $\tilde{Z}^{(\ell)}_{it}=(0\ 0\ \ldots 1 \ 0 \ldots 0)'$ is a $Q^{(\ell)}$--vector of the standard basis with one in the $Z^ {(\ell)}_{it}$--th position, $i,j\in \mathcal{V}$, $t\in \mathcal{T}$, $Z^{(\ell)}_{it}\in \mathcal{Q}^{(\ell)}$. For example, if the weights are log--normally distributed, the $\theta^{(\ell)}_{Z^{(\ell)}_{it}Z^{(\ell)}_{jt}}=(\beta^{(\ell)\prime}_{qr},\sigma_{qr}^{(\ell)2})'$ includes: i) a vector $\beta^{(\ell)}_{Z^{(\ell)}_{it}Z^{(\ell)}_{jt}}=(\tilde{Z}^{(\ell)\prime}_{it}\beta_{0}^{(\ell)}\tilde{Z}^{(\ell)}_{jt},\ldots,\tilde{Z}^{(\ell)\prime}_{it}\beta_{S-1}^{(\ell)}\tilde{Z}^{(\ell)}_{jt})'$ of intercept and the coefficients of the covariates $(X^{(\ell)}_{ijt})$; and ii) the variance $\sigma^{(\ell)2}_{Z^{(\ell)}_{it}Z^{(\ell)}_{jt}}=\tilde{Z}^{(\ell)\prime}_{it}\sigma^{(\ell)}\tilde{Z}^{(\ell)}_{jt}$.} In this sense, the DSBMM model can also be interpreted as a regression model with parameters partially pooled across dyads and time. Consequently, ((ref)) is more flexible than a standard fixed parameters regression, and more feasible and parsimonious than a time--dyad varying case $\vartheta_{ijt}^{(\ell)}$. Thus, the DSBMM extends the dynamic single--layer SBM olivella2022dynamic and the static SBM with observed heterogeneity mariadassou2010uncovering to a dynamic multidimensional setup with time--varying node or dyad characteristics, $X^{(\ell)}_{ijt}$. In other words, apart from introducing inter--layer dependence, in the one--layer case $L=1$, the DSBMM in ((ref)) differs from yang2011detecting,matias2017statistical in accounting for node and dyad observed heterogeneity through $X^{(\ell)}_{ijt}$.

We introduce a different representation of the model, which offers many advantages. Let $D_{ijt}^{(\ell)}$ be an indicator variable following a Bernoulli distribution, such that $D_{ijt}^{(\ell)}=1$ if $(i,j)\in \mathcal{E}_t^{(\ell)}$ and $D_{ijt}^{(\ell)}=0$ if $(i,j)\not\in \mathcal{E}_t^{(\ell)}$, then ((ref)) can be rewritten as

\textcolor{black}{A first advantage is that ((ref)) and ((ref)) account for both weighted and unweighted layers. In case of unweighted layers, equations ((ref)) and ((ref)) simplify to ((ref)). The second advantage concerns the model's flexibility regarding the edge's existence. If $Y^{(\ell)}_{i j t}=0$ indicates the absence of an edge and not a missing value, then this representation provides a two--part network formation process and the common stochastic membership accounts for possible dependence between the two parts. If all zeros are missing data, $D^{(\ell)}_{ijt}=1$ refers to an observed $Y^{(\ell)}_{i j t}$ and $D^{(\ell)}_{ijt}=0$ denotes missing values. This representation can be interpreted as a random process for the missing values, driven by a common stochastic membership. Assuming dependence between the missing value mechanism and the observable data is common in many applied network areas, from sample selection or common factor approaches, and extending pattern mixture helpman2008estimating, little1993pattern,song2004imputation. The third advantage is that this representation can be used to include covariates in (ref), i.e. $\tau\left(\nu^{(\ell)}_{qr}\right)= X^{(\ell)'}_{ijt}\theta^{(\ell)}_{qr}$, where $\tau:(0,1)\to\mathbb{R}$ is a link function. Indeed, other latent components, apart from the stochastic membership, can be used as in Heckman's procedure or common factor approaches van2011bayesian. }

{\color{black} The specification ((ref)) implies not only a time dependence through the HMCs, but a cross--sectional edge dependence. In Proposition (ref), this property is exemplified under specific distributional assumptions. The edges within the same community pair interaction have a covariance different from zero, while the edges with different community pair interactions have zero covariance, i.e. a block covariance structure. Moreover, the discrete mixture also provides more flexibility in the relationship between the first and second moments of the weights. In particular, while the log--normal distribution of the weights implies a positive quadratic relationship between mean and variance within community pair interaction, the mixture allows for discontinuities that accommodate monotonic and non--monotonic jumps in the relationship mean--variance.

The conditional expected weights follow a Pareto--log--normal distribution, a heavy--tail distribution used in network analysis to describe the degree distribution and also in international trade theory to capture the productivity heterogeneity of the firms at a granular level nigai2017tale,fang2011double. As in the Pareto distribution, the Pareto--log--normal distribution's tail properties would be preserved under convolutions, meaning that the behavior of the extreme values in the network strength would inherit similar characteristics. The variance--to--mean ratio (VM) measures the overdispersion features and is given in the following.

As $\underline{e}^{(\ell)}$ goes to infinity, the moments and VM ratio gets closer to the log--normal case, while as $\underline{e}^{(\ell)}\to2^+$ the VM ratio explodes due to the heavy tails, with tail index of $1/\underline{e}^{(\ell)}$. These properties can be used to calibrate the hyperparameters of the prior in applications, where prior information on the tails of the strength is available.

\textcolor{black}{Notice that by specifying inter--layer edge conditional distribution given the layer blocks, the model framework and properties can be easily extended to a general multi--layer networks, where layers are node aligned and edges between layers are allowed. Since this specification is not relevant to our application, we leave it for further research.} }

Layer Dependence

A statistical model for multidimensional networks should account for edge redundancy, that is, the persistence of edges between the same pair of nodes across network layers, and more structurally, clustering redundancy han2015consistent,stanley2016clustering,jovanovski2019bayesian. For instance, if a set of nodes is densely connected in one layer, a similar structure is highly likely to be present in the second layer. Nevertheless, existing approaches are static and correlational and only account for clustering overlap. To capture dependence between the node partitions in the different layers, including non--overlap dependence, and identify a \textcolor{black}{NGB} causal structure, we can exploit the time dimension and assume dependent Markov Chain processes otranto2005multi,agudze2021markov.

We are interested in a Granger non--causality relationship across layers. In this respect, the transition probabilities for each layer only depend on the collection of memberships $\mathfrak{V}_{t-1}$, that is $Z_{it-1}=\textcolor{black}{(Z^{(\ell)}_{it-1})_{\ell \in \mathcal{L}}}$. This assumption simplifies as follows the joint probability $\mathbb{P}(Z_{it}|Z_{it-1})$, $i\in \mathcal{V}$,

In contrast to previous studies, which typically cover dependent Markov chains with only two states, the DSBMM requires a more general approach. Therefore, we use a multinomial logit form to express the Markov chain dependence through the transition matrix $\mathbb{P}(Z^{(\ell)}_{it}=r|Z_{it-1})=P^{(\ell)}_{it,Z^{(\ell)}_{it-1}r}, \ r=1,\dots, Q^{(\ell)}$ with $Q^{(\ell)}$ states. Similarly to the log--linear models used for contingency tables, $Z_{it-1}$ is represented by a full set of dummy variables, \textcolor{black}{from now on lagged HMC variables}, $W^{(\ell)}_{it,q}=\mathbb{I}_{\{q\}}\left(Z^{(\ell)}_{it}\right)$, and $W^{(\ell)}_{it}=(W^{(\ell)}_{it,1},\ldots, W^{(\ell)}_{it, Q^{(\ell)}-1})'$, including the interactions and the main effects of the different layers, that is

where $\kappa^{(\ell)}_{\mathcal{U},r}=\left(\kappa^{(\ell)}_{\mathcal{U},r,1},\ldots,\kappa^{(\ell)}_{\mathcal{U},r,s(\mathcal{U})}\right)'$ \textcolor{black}{is the collection of parameters measuring the effect of past membership values on the transition to community $r$ in layer $\ell$, and $s(\mathcal{U})=\prod_{\mathfrak{m}\in\mathcal{U}}\left(Q^{(\mathfrak{m})}-1\right)$ is the size of the vector, and $ \mathcal{U}\in \mathcal{P}$ identifies the set of interacting layers with $\mathcal{P}$ the power set of the layer set $\mathcal{L}$}. \textcolor{black}{For instance, the main effects of the layer 1 on the log odds of changing to block $r$ in the layer $\ell$, i.e. $\mathcal{U}=\{1\}$, is captured by the vector $\kappa^{(\ell)}_{\{1\},r}=\left(\kappa^{(\ell)}_{\{1\},r,1},\ldots,\kappa^{(\ell)}_{\{1\},r,s(\{1\})}\right)'$, where the number of elements is related to the number of communities in layer $s(\{1\})=Q^{(1)}-1$. Following the example, the first--order interaction effects of layer 1 and layer $\ell$ are collected in the vector $\kappa^{(\ell)}_{\{\ell,1\},r}=\left(\kappa^{(\ell)}_{\{\ell,1\},r,1},\ldots,\kappa^{(\ell)}_{\{\ell,1\},r,s(\{\ell,1\})}\right)'$, where $s(\{\ell,1\})=(Q^{(1)}-1)(Q^{(\ell)}-1)$.} On the left side, $C^{(\ell)}_{it}=\sum_{k\in \mathcal{Q}^{(\ell)} }\exp\left(\widetilde{W}_{it-1}\kappa_k^{(\ell)}\right)$ is the normalizing constant of the multinomial model, where $\widetilde{W}_{it-1}=$ $(1,W^{(\ell)'}_{it-1}$, $\ldots$, $\bigotimes_{\mathfrak{m}=1}^{L}(W^{(\mathfrak{m})'}_{it-1}))$ and $\kappa_{r}^{(\ell)}=\left(\kappa^{(\ell)}_{0,r},\ldots,\kappa^{(\ell)'}_{1\ldots L,r}\right)'$ are a row and a column vector, respectively, with $p+1$ elements, with $p=\sum_{\mathcal{U}\in \mathcal{P}} s(\mathcal{U})$. Equation ((ref)) can be rewritten in an equivalent, more compact form:

For identification, $Z^{(\ell)}_{it}=Q^{(\ell)}$ can be used as reference state/community, i.e. $\kappa^{(\ell)}_{Q^{(\ell)}}=\mathbf{0}$. This multivariate logistic specification is more flexible than the linear specification proposed by agudze2021markov for a panel Markov switching model because it does not have to impose any restriction on the parameter space of $\kappa_{r}^{(\ell)}$. Consequently, ((ref)) captures the special case of clustering overlap, as well as more complex dependence structures, such as the decoupling one. In this modeling framework, it is straightforward to add other exogenous variables that affect memberships and the transition matrix kaufmann2015k,holsclaw2017bayesian. This flexibility in ((ref)) comes at a cost of lower tractability and a larger number of parameters. In this paper, we provide a solution to both issues based on a suitable inference approach.

We extend to the block membership processes the definition of non--causality introduced by mosconi2006non. Let $\{Z_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}}=(Z^{(1)}_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}},\dots,Z^{(L)}_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}})', \ t\in \mathcal{T}\}$, or simply $\{Z_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}}\}$, be a random sequence in the probability space with the triplet $(\Xi, \mathcal{A},\mathbb{P})$\textcolor{black}{, where $\Xi$ refers to the sample space, that is the potential memberships in all layers, $\mathcal{A}$ the sigma algebra of events, in other words the potential trajectories of the memberships and $\mathbb{P}$ its probability measure}. An information set and a restricted information set are required to define a non--causality property. The available information up to time $t$, in terms of $\{Z_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}t}}\}$, is referred to as canonical filtration $\left\{\mathcal{F}_t, \ t\in \mathcal{T}\right\}$, which is a sub--$\sigma$--algebra of $\mathcal{A}$. The reduced information sets are $\mathcal{R}^{(-\ell)}_t=\sigma((Z_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}s}}^{(-\ell)}),\ 1\leq s\leq t)$ and $\mathcal{Z}^{(\ell)}_t=\sigma(Z_{{\color{black}\mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}s}}^{(\ell)},\ 1\leq s\leq t)$, for a layer $\ell\in \mathcal{L}$, so that $\mathcal{Z}^{(\ell)}_t\subseteq\mathcal{R}^{(-\mathfrak{m})}_t\subseteq\mathcal{F}_t$ for all $\mathfrak{m}\neq \ell$.

Following the specification ((ref)), Layer $\mathfrak{m}$ does not \textcolor{black}{NGB} cause Layer $\ell$ if all the main effects and interactions parameters $\kappa^{(\ell)}_{\mathcal{U},r}$ involving $\mathfrak{m}$ in the equation ((ref)) are null, otherwise $\mathfrak{m}$ does \textcolor{black}{NGB} cause Layer $\ell$. \textcolor{black}{Granger causality hypothesis testing based on this definition and the parameterization used will be introduced later on in Section (ref).}

{\color{black} Based on Definition (ref), the potential scenarios of dependence, as in a VAR, are: a unidirectional causality ($Z^{(\mathfrak{m})}\stackrel{G}{\Rightarrow} Z^{(\ell)}$ and $Z^{(\ell)}\not \stackrel{G}{\Rightarrow} Z^{(\mathfrak{m})}$ or vice--versa), bidirectional causality ($Z^{(\mathfrak{m})}\stackrel{G}{\Rightarrow} Z^{(\ell)}$ and $Z^{(\ell)}\stackrel{G}{\Rightarrow} Z^{(\mathfrak{m})}$) or absence of causality ($Z^{(\mathfrak{m})}\not\stackrel{G}{\Rightarrow} Z^{(\ell)}$ and $Z^{(\ell)}\not\stackrel{G}{\Rightarrow} Z^{(\mathfrak{m})}$)}. An illustration of the unidirectional \textcolor{black}{NGB causality} is presented in Figure (ref). It shows a two--layer undirected and unweighted network with nine nodes clustered into three communities, indicated by different gray shades. At time $t_1$ (panel a), the partition of the nodes in Layer 1 differs significantly from that in Layer 2. However, over time the partition in Layer 2 aligns with the one in Layer 1, producing a clustering/community overlap. Based on the definition, Layer 1 is \textcolor{black}{NGB}, causing Layer 2. This membership coupling effect is only an example, as other coupling or decoupling scenarios are possible within our DSBMM (see Section (ref)).

\textcolor{black}{Our definition of causality shares some similarities with and differs in many aspects from other notions of blockwise Granger causality dufour2010short,hu2015shortcomings. These latter are linear since they have been proposed in a VAR setup. Moreover, they refer to a relationship between groups of variables (blocks) at different lags. In our DSBMM there are collections of membership variables within each layer and a non--linear relationship between them is assumed. Thus, our NGB for DSBMM is naturally a non--linear causality notion, and regards the relationship between groups of membership variables in different layers (blocks) at different lags. The blocks of membership processes should not be confused with the blocks of nodes, nevertheless one can expect that the NGB in the membership influences the dynamics in the node block structure. As in the conditional Granger causality for linear models, our NGB controls for other blocks (i.e., layers) mediated effects, including in the conditional relationships and in the testing procedure the other--layer node memberships with lags. Finally, in our NGB, unobserved factors influencing node clustering in layer $\mathfrak{m}$ are related to their counterparts in Layer $\ell$. Thus, a substantial difference between NGB and standard Granger causality notions, is that NGB is defined for latent processes, whereas standard notions refer to observable processes.}

We illustrate our DSBMM through simulations in two settings: unidirectional and bidirectional \textcolor{black}{NGB causality} (see Settings (ref) and (ref) in Section (ref) of the Supplementary Material for further details). In these two settings, the network has three layers. Layers 1 and 2 are weighted and directed and have two and three communities, respectively. Layer 3 is unweighted and undirected and has three communities. The number of nodes is $N=50$, and the time horizon is $T=15$. For simplicity no regressors are included and the weighted layers are log--normal distributed with mean $\beta^{(\ell)}_{qr0}$ and variance $\sigma^{(\ell)2}_{qr}$, that is $f^{(\ell)}(y| \theta^{(\ell)}_{qr},X^{(\ell)}_{ijt})$ is equivalent to $\operatorname{LN}(\beta^{(\ell)}_{qr0},\sigma^{(\ell)2}_{qr})$.

The unidirectional configuration returns the membership dynamics given in Panel (a) of Figure (ref). The alluvial plots show that memberships remain stable in Layer 3, and after 15 periods, the blocks in Layer 2 are almost aligned with those in Layer 3; that is, Layer 3 is \textcolor{black}{NGB}, which causes Layer 2. In the bidirectional case (Panel b) the memberships change significantly in both layers and after 15 periods the blocks in the two layers are aligned.

Bayesian Inference

The prior choices are discussed in this section, and the full conditional posterior distributions are derived together with the Markov chain Monte Carlo (MCMC) approximation algorithm. The proofs of the results are given in Appendix (ref)

Prior specification

Given the presence of latent variables in the DSBMM, the observed likelihood is not tractable. Moreover, the specification ((ref)) introduces all possible interaction effects between layers, causing over--parameterization. A Bayesian approach addresses these issues and provides uncertainty measures of the latent variables and the parameters such as credible intervals yang2011detecting. The prior distributions used for the connectivity parameters are

which belong to the Normal, Inverse--Gamma, Beta, and Dirichlet families, respectively. These distributions are assumed to be independent and are commonly used in the SBM literature, given their convenience as conditional conjugate yang2011detecting,lee2019review.

Regarding the transition parameters of a Markov chain, the standard choice is a Dirichlet prior. Nevertheless, in our DSBMM, the multinomial representation of the HMC dependence calls for the use of an informative prior over the real line, such as the Normal distribution held2006bayesian, fruhwirth2010data,billio2016interconnections. The Normal prior induces a ridge shrinkage, pushing the estimates of the transition probabilities away from zero and one, which is particularly useful for the entries of the transition matrix with few to no observations. However, this prior choice may still be unsatisfactory because it does not cope with the over--parameterization resulting from including all layers and order effects in ((ref)). Since the main interest of DSBMM is to select the layers that have a significant influence on transition entries and not to select individual variables, applying a global shrinkage or a global variable selection method would omit important information in considering the correlation between \textcolor{black}{lagged HMC variables}---i.e., variable grouping. We follow a Bayesian group--LASSO, through a Multi--Laplacian prior raman2009bayesian.

The Multi--Laplacian distribution selects among groups of covariates and shrinks the parameters of highly correlated variables within each subset $\mathcal{U}$. The precision $c(\mathcal{U})$ governs the sparsity level in the parameters $\kappa_{\mathcal{U}}$. If all lagged HMC variable groups are singletons, i.e., $s(\mathcal{U})=1$ for all $\mathcal{U} \in \mathcal{K}$, the group LASSO is equivalent to a standard Bayesian LASSO. As in the classical LASSO, the objective is to identify the most relevant variables, but at a group level.

We assume a Multi--Laplacian prior for the parameters of the higher--order effects and of the main effect of other layers. The Normal prior is assumed for the remaining parameters. As stated in the following, our prior distribution for $\kappa^ {(\ell)}$ can be represented as a continuous scale mixture of Normal distributions with a Gamma mixing distribution park2008bayesian.

The shrinking in ((ref)) is twofold: shrinking order and layer effects following a $l_1$--penalty to choose the most relevant groups and within--group shrinking following a $l_2$--penalty to deal with highly correlated variables and entries with few to no transitions. The diagonal elements of the prior variance $\underline{K}_r^{(\ell)}$ are groupwise sampled from a Gamma distribution conditionally on $\rho^{(l)}$. Comparing ((ref)) with the original proposal by raman2009bayesian, two changes are introduced. First, contrary to a Poisson model for contingency tables, a multinomial logit implies $Q^{(\ell)}-1$ regressions, all of them sharing the same parameter $\rho^{(\ell)}$. In this way, variable selection applies even across equations. This cross--equation information improves the estimation of $\rho^{(\ell)}$ by jointly contrasting more groups of variables. Second, not only is the intercept of each of the $Q^{(\ell)}-1$ regression excluded from the Multi--Laplacian prior, but also the main effects of the own layer in order to capture the HMC persistence. As a consequence, the variable selection does not apply to the intercept and the own lagged membership parameters, $Z^{(\ell)}_{it-1}$.

The variable selection can be sensitive to the choice of $\rho^{(\ell)}$ because this latter has a direct relationship with the Lagrange parameter of the LASSO regression, $c^{(\ell)}$. Thus, a Gamma prior is used for $\rho^{(\ell)}$ with hyperparameters $\iota^{(\ell)}_{1}$ and $\iota^{(\ell)}_{2}$

The Bayesian interpretation of penalizations as shrinkage priors allows for the estimation of the shrinkage parameter $c^{(\ell)}$ (through $\rho^{(\ell)}$) as a natural hierarchical extension and it provides valid standard errors and other measures of uncertainty casella2010penalized. {\color{black} Since the focus of the paper is on lagged causality, we assume independent priors for the layer parameters. Nevertheless, cross--layer contemporaneous dependence can be incorporated by assuming partial parameter pooling or hierarchical priors.}

Figure (ref) summarizes the DAG of the DSBMM (only for a layer $\ell$ for the sake of simplicity). This representation underlines the influence of changes in community membership on the topological properties observed in the network, which in turn is affected by the community membership in the rest of the layers $Z^{(-\ell)}_{it}$. This illustration excludes the hyperparameter of the prior distribution of the connectivity parameters and the initial partition of the nodes $\alpha^{(\ell)}$ to focus on the transition parameters. The priors for these latter, that are the Normal and Multi--Laplacian distributions, are presented as a heteroschedastic Normal prior.

Posterior approximation

Consistently with our application to trade flows, we assume $f^{(\ell)}(y| \theta^{(\ell)}_{qr},X^{(\ell)}_{ijt})$ is a log--normal distribution $\operatorname{LN}(\beta^{(\ell)}_{qr},\sigma^{(\ell)2}_{qr})$. Nevertheless, our modeling framework is general and can be easily modified to account for other distribution assumptions, such as a zero--truncated Poisson for count data matias2017statistical. Since the likelihood function is not tractable, a Gibbs sampling procedure is derived to approximate the posterior distribution.

Using the prior distributions presented in Section (ref), the full conditional posterior distributions of the connectivity parameters and initial proportions are given in the following propositions. A full Gibbs sampling algorithm using the conditional distributions is then used to approximate the posterior distribution.

The multinomial representation of DSBMM does not allow for a conditional conjugate prior, and the use of the standard Metropolis--Hastings (MH) algorithm can be highly inefficient. As suggested by fruhwirth2010data, data augmentation can be applied to avoid MHs steps. Therefore, in the following propositions we extend the Pólya gamma representation proposed by holsclaw2017bayesian and polson2013bayesian to the saturated multinomial logit and obtain tractable full conditional distributions for the transition parameters $\kappa_{r}^{(\ell)}$, the auxiliary variables of the two data augmentations, $\omega^{(\ell)}_{it,r}$ and $\zeta_{r\mathcal{U}}^{(l)2}$, and the hyperparameter related to the level of shrinkage $\rho^{(l)}$.

Regarding the HMC, we compute iteratively the conditional prediction probability $\mathbb{P}(Z^{(\ell)}_{it}=r|Z^{^{(-\ell)}}_{i,1:t}, Z^{^{(\ell)}}_{-i,1:t-1},\psi^{(\ell)}_{t-1})$, the conditional filtered probability $\mathbb{P}(Z^{(\ell)}_{it}=r|Z^{^{(-\ell)}}_{i,1:t+1},Z^{^{(\ell)}}_{-i,1:t},\psi^{(\ell)}_{t})$ and sample $Z^{(\ell)}_{i,t}$ from the smoothed probability $\mathbb{P}(Z^{(\ell)}_{i,1:T}|Z^{(-\ell)}_{i,1:T},$ $Z^{(\ell)}_{-i,1:T},\psi^{(\ell)}_{T})$, where $\psi^{(\ell)}_t=[Y^{(\ell)}_{1:t},D^{(\ell)}_{1:t},X^{(\ell)}_{1:t}]$ fruhwirth2006finite,touloupou2020scalable. See Section (ref) for further details. The last step in the Gibbs sampler generates samples for the allocation variables $Z^{(\ell)}_{{\color{black} i \mathbin{\vcenter{\hbox{\scalebox{.4}{$\,\bullet\,$}}}}}}$ by using a conditional forward filtering and backward sampling (FFBS). \textcolor{black}{For applications where the node set is time--varying, the membership dynamics can be reconstructed by assuming a node--centered sampling design and by following the missing handling strategy based on FFBS proposed in tabouy2020variational,hamaker2012regime.}

In mixture models and HMC models, state labels are not identifiable, and label changes may occur between Gibbs iterations, which is referred to as label switching. The DSBMs belong to the class HMC models, and the MCMC approximation can be affected by this issue matias2017statistical. Prior constraints that fix a specific cluster order can be used to achieve identification. We follow an alternative approach based on post--processing MCMC samples fruhwirth2006finite.

Model selection and causality testing

Regarding variable selection, a procedure based on the posterior estimates of the group LASSO is not straightforward because it does not lead to exact zeroes for the posterior modes or means. The alternative used in this paper is a credible interval criterion, that is, a parameter whose credible interval includes zero is considered to be null van2019shrinkage.

From a Bayesian perspective, a credible interval criterion has three main drawbacks: the results depend on the chosen credible level, it does not express the different sources of uncertainty, and does not account for the joint posterior distribution of the parameters. However, in a Bayesian LASSO framework, calculating Bayes factors or posterior variable inclusion probabilities is computationally intensive for a saturated multinomial logit hans2010model. Alternatively, spike-and-slab prior distributions can be used to explore the model set george1997approaches. Nevertheless, this space is very large, and the search procedure is either computationally intensive or inefficient. Additionally, a threshold is required to induce sparsity, and the results are typically sensitive to the choice of this threshold. This selective inference issue can be even more severe for highly correlated \textcolor{black}{lagged HMC variables} since it does not consider the joint distribution of the parameters bondell2012consistent. Therefore, considering its computational convenience and given that its drawbacks are also challenging for alternative methods, the credible interval criterion appears to be a suitable option for large models, such as the DSBMM.

{\color{black}Regarding NGB causality testing, an advantage of our approach is that NGB causality can be tested through a few parameter restrictions and a theoretically valid procedure based on the Bayes factor (BF). In the multiple-layer case, blockwise testing is still feasible and can be applied to pairs of layers. Let us denote with $\kappa^{(\ell)}_{\mathcal{U}}=(\kappa^{(\ell)\prime}_{\mathcal{U},1},\ldots, \kappa^{(\ell)\prime}_{\mathcal{U}, Q^{(\ell)}-1})'$ and $\kappa^{(\ell)}_{\mathcal{U},r}=(\kappa^{(\ell)}_{\mathcal{U},r,1},\ldots,\kappa^{(\ell)}_{\mathcal{U},r,s(\mathcal{U})})'$ the collections of transition parameters in the membership probabilities of the nodes in layer $\ell$. Following Definition 1 and the Markov--chain conditionally independence assumption in (ref), the layer $\mathfrak{m}$ does NGB cause layer $\ell$, i.e. $Z^{(\mathfrak{m})} \stackrel{G}{\Rightarrow} Z^{(\ell)}$, if the null hypothesis $H_0:\kappa^{(\ell)}_{\mathcal{U}}={\bf 0}, \ \forall\ \mathcal{U}\in \mathcal{P}^{(\mathfrak{m})}$ is not satisfied, where $\mathcal{P}^{(\mathfrak{m})}$ denotes the parameters subscript indexes associated to the main effects and interactions of layer $\mathfrak{m}$ \textcolor{black}{lagged HMC variables}, i.e. $\mathcal{P}^{(\mathfrak{m})}=\{ \mathcal{U}\in \mathcal{P}| \mathfrak{m}\in \mathcal{U}\}$. For example, in specifications with three blocks in layer $\ell$ and two blocks in layer $\mathfrak{m}$, layer $\mathfrak{m}$ does not NGB causes layer $\ell$ if the null $H_0:\kappa^{(\ell)}_{\{\mathfrak{m}\},1,1}=\kappa^{(\ell)}_{\{\mathfrak{m}\},2,1}=0\, (\text{main effects}),\ \kappa^{(\ell)}_{\{\ell,\mathfrak{m}\},1,1}=\kappa^{(\ell)}_{\{\ell,\mathfrak{m}\},1,2}=\kappa^{(\ell)}_{\{\ell,\mathfrak{m}\},2,1}=\kappa^{(\ell)}_{\{\ell,\mathfrak{m}\},2,2}=0\ (\text{interaction effects})$. The Bayes factor can be used for Bayesian testing, which involves the ratio of the posterior probability of the hypothesis to its prior probability. By assuming the same prior probabilities for the two hypotheses, the Bayes factor can be expressed as a ratio of the marginal likelihoods under the null and the alternative, i.e. $$ BF_{1,0}(Y)=\frac{m(Y|H_1)}{m(Y|H_0)}, $$ where $m(Y|H_i)=\int_{\Theta \times K} \sum_{\mathcal{Z}}h_{i}d\vartheta d\kappa$ for $i=0,1$ and $h_{i}=L(Y,Z|X,\vartheta,\kappa) \pi_{i}(\vartheta,\kappa)$. Notice that $\pi_{0}(\vartheta,\kappa)$ corresponds to the prior under the null hypothesis, where the prior of the restricted parameters $\kappa^{(\ell)}_{\mathcal{U}}={\bf 0}, \ \forall\ \mathcal{U}\in \mathcal{P}^{(\mathfrak{m})}$ reduces to a Dirac distribution at zero, while in the $\pi_{1}(\vartheta,\kappa)$ (alternative hypothesis) corresponds to the Multi--Laplacian prior.

The marginal posterior can be approximated using MCMC draws (see llorente2023marginal for a more extensive review). As in carallo2024generalized, the method proposed by geyer1991estimating can be used to approximate log--Bayes factor by expressing the posterior draws of the restricted and unrestricted models as a mixture, which can be written as a reverse logistic regression (RLR) in Proposition (ref). \setcounter{proposition}{5}

In other words, it is equivalent to maximizing the quasi--likelihood of the mixture with respect to $d$, which can be implemented using a standard logit regression with an intercept and a regressor that is the same for $i=\{0,1\}$ and whose slope parameter is restricted to be one. }

Simulation results

{\color{black} We study the effectiveness of the model under different DSBMM parameter prior assumptions (Normal and group LASSO priors), different non--linear Granger causality settings (no causality and unidirectional or bidirectional causality), and causality signs (coupling and decoupling layers). We also considered different network topologies (core--periphery and assortative structures) and sample size scenarios (a reference scenario, $N=50$ and $T=15$, a large $N$ scenario, $N=100$ and $T=15$, and a large $T$ scenario, $N=50$ and $T=30$). In all causality settings and sample scenarios, the DSBMM combined with a Bayesian group LASSO prior (Multi--Laplacian prior) outperforms the Normal prior in correctly retrieving the NGB causality structure (see Table (ref) in Section (ref) of the Supplementary Material). We also provide some guidelines to choose the hyper--parameters of the Multi--Laplacian prior.}

{\color{black} We also present a comparison between the DSBMM and alternative models from time--series analysis canova2013panel,koop2013forecasting and dynamic network literature: the Bayesian Learning for Dynamic network (BLDMN) introduced by durante2017bayesian and the layer--independent DSBM introduced by matias2017statistical. The results in Table (ref) show that BVAR models have difficulties in detecting NGB causality. Additionally, the DSBMM outperformed other network models, especially out--of--sample (see Table (ref)). This outcome is due to the DSBMM's ability to incorporate information across layers, while capturing the topological change between an assortative and a core--periphery structure. More importantly, the DSBMM provides three modeling advantages: a conditional directed relationship by layer--pair (NGB causality), a general framework that jointly considers (un)directed and (un)weighted layers, and an intuitive interpretation in terms of clustering and classification of the nodes.}

FTAs and trade networks

Our DSBMM model and inference provide a unified and coherent framework for international trade analysis, merging Gravity models and community detection algorithms. Since the General Agreement of Tariffs and Trade, FTAs have been an important tool conceived as a means to increase trade flows among countries. We use the DSBMM to infer the {\color{black}NGB} causal relationship between FTAs and trade flows.

Gravity model, communities and the DSBMM

Empirical studies have used the Gravity equation to explain bilateral trade and infer the effect of trade policies. This equation suggests that the economy's size (e.g., GDP, population) and the frictions (e.g., physical distance, language similarity, colonial links) of the countries involved determine trade flows. Although this gravity framework started as an intuitive empirical regularity, a theoretical foundation is possible using a general equilibrium model in a static framework à la Heckscher--Ohlin with monopolistic competition anderson2011gravity, kabir2017gravity. These models extend the interpretation of trade flows beyond dyad relationships and focus on the structure of the trade network. In its reduced form, this family of models can be written as

where the trade flow $Y^{(1)}_{ijt}$ increases with the node characteristics $GDP_{it}$ and $GDP_{jt}$, i.e. $\beta^{(1)}_{1},\beta^{(1)}_{2}>0$, and the frictions $t_{ijt}$ have a negative effect on trade. The two fully unobserved elements $\Upsilon_{jt}$ and $\Pi_{it}$ positively affecting the trade flow between $i$ and $j$ at time $t$, called inward and outward multi--resistance respectively, encompass an aggregated measure of the barriers between all the pairs of nodes in the trade network. {\color{black}The barriers are partially observed and are usually approximated by a set of observable variables ($T_{ijt}$), such as export subsidies and tariffs between countries or by unobserved effects ($\tilde{t}_{ijt}$), i.e., $t_{ijt}=\exp(T'_{ijt}\lambda+\tilde{t}_{ijt})$. To evaluate specific economic policies the main parameter of interest may include elements of the vector $(1-\varrho)\lambda$ and the reconstruction of the barrier network $(t_{ijt})^{1-\varrho}$, which allows to estimate general equilibrium impact of trade policies after solving for the multi--resistance terms in (\refeq{eq:grav}) yotov2016advanced. The parameter $\varrho$ can be interpreted as an elasticity of substitution and is usually estimated if the tariff data is available or calibrated based on previous studies anderson2020transitional. While the DSBMM can be used to estimate the structural parameter of interest, our main objective is to address parameter heterogeneity and its relationship with FTAs through the multidimensional dependence and the NGB causality.}

In the DSBMM, the parameter heterogeneity can be present across the network and time. Previous works have underlined the importance of both dimensions. For instance, yotov2012simple demonstrates the decreasing influence of geographical distance on trade over time due to globalization, but assumes the marginal effect of distance remains constant for all countries. On the other hand, based on firms' behavior, technology, market demand, and economic policies, the literature on trade has suggested specific network structures and potential dyad clustering. In particular, under some circumstances, the multi--resistance terms may not include information on all the network, but only on the barriers of the neighboring nodes of country $i$ and $j$, as in the Gravity model proposed by magerman2020pecking\textcolor{black}{, where the neighboring countries of country $j$ importing from country $i$ refers to the set of active exporters $\{k\in \mathcal{V}:Y_{kj\textcolor{black}{t}}>0\}$}. This latter model predicts the existence of two communities in the trade network: a core and a periphery, resulting from firms' behavior, which determines their export destinations based on a pecking order. Other authors suggest that firms choose new destinations that are similar to previous ones, given the high cost of adapting products to heterogeneous markets, which may implicitly produce blocks of countries. Based on these considerations, a block structure in the network of barriers $t_{ijt}$ can be directly creating communities on $Y^{(1)}_{ij\textcolor{black}{t}}$, and/or indirectly through $\Upsilon_{j\textcolor{black}{t}}$ and $\Pi_{i\textcolor{black}{t}}$, while heterogeneous policy effects across the clusters may be also a possibility. For instance, this motivates the empirical work of sopranzetti2018overlapping that studies the overlap between FTAs and trade networks from a dyad perspective.

Similarly, the empirical literature on community detection has provided evidence of clustering in trade and FTA networks, but it has mainly employed algorithmic methods, such as maximizing a modularity measure, to identify groups of countries with dense subgraphs, as in barigozzi2011identifying,bartesaghi2020communicability. In this sense, DSBMM provides a statistical model, and subsequently, its inference includes a measure of uncertainty of the parameters and the latent membership. More importantly, it allows for dynamic membership of nodes and covariates that control for node characteristics, which are essential to link the community structure with the Gravity equation and parameter instability.

Specifically, we suggest that detecting communities in the trade flows after controlling for the observed covariates in ((ref)) concentrates the attention on the block structure in the unobserved components, that is, the multi--resistance terms and the unobserved barriers to trade. In order to focus on the multidimensional dependence, we assume only geographical distance is available to approximate the observed barriers ($dist_{ij}$), i.e. $t_{ij\textcolor{black}{t}}=\exp(\lambda_{Z_{it}Z_{jt}} \textcolor{black}{\log} \ dist_{ij}+\tilde{t}_{Z_{it}Z_{jt}})$. The effect of distance ($\lambda_{Z_{it}Z_{jt}}$) and the unobserved barriers are block--dependent ($\tilde{t}_{Z_{it}Z_{jt}}$) and can vary across dyad and time, resulting in a partial pooling of the linear model, that is the first layer of a multidimensional network with $L=2$ is given by

where \textcolor{black}{the covariates are} $X_{ijt}=(1,\textcolor{black}{\log}\ GDP_{it},\textcolor{black}{\log} \ GDP_{jt},\textcolor{black}{\log} \ dist_{ij})'$, $\beta_{qr}^{(1)}=(\check{\beta}^{(1)}_{0qr},\beta^{(1)}_{1qr},\beta^{(1)}_{2qr},\beta^{(1)}_{3qr})'$, $\beta^{(1)}_{3qr}=(1-\varrho)\lambda_{qr}$, $\check{\beta}^{(1)}_{0qr}=\beta^{(1)}_{0qr}+(1-\varrho) \tilde{t}_{qr}-(1-\varrho)\textcolor{black}{\log}\ \Upsilon_r-(1-\varrho)\textcolor{black}{\log}\ \Pi_q$ and $q,r \in \mathcal{Q}^{(1)}$. Notice that zero--trade cases are allowed with probability $\nu^{(1)}_{qr}$ depending on the community structure.

The second layer corresponds to the unweighted and directed FTA network, and it has its own community structure

Since FTAs are part of the unobserved frictions $\tilde{t}_{qr}$, this specification ((ref)) and ((ref)) can be used to test if the community structure in the FTA network is {\color{black}NGB} causing the block membership in the trade network. For example, if $Q^{(1)}=Q^{(2)}=2$, ((ref)) becomes

and the \textcolor{black}{NGB causality} tests can be stated as: FTAs does not \textcolor{black}{NGB} cause unobserved barriers and multi--resistance terms network, $H_0:\kappa^{(1)}_{\{2\},1,1}=\kappa^{(1)}_{\{1,2\},1,1}=0$; or, trade network does not \textcolor{black}{NGB} cause FTAs network, $H_0:\kappa^{(2)}_{\{1\},1,1}=\kappa^{(2)}_{\{1,2\},1,1}=0$. A third scenario would imply a bidirectional causality, $H_0:\kappa^{(1)}_{\{2\},1,1}=\kappa^{(1)}_{\{1,2\},1,1}=\kappa^{(2)}_{\{1\},1,1}=\kappa^{(2)}_{\{1,2\},1,1}=0$.

The FTA network is a trade policy instrument and may help predict the trade flows layer. However, this is not the only possible scenario because FTAs are not fully exogenous. As posed by baier2004economic's theoretical model, the welfare expectations of consumers and firms may depend on the pre--policy levels of trade flows and affect the likelihood of an FTA. Moreover, integration processes, domestic regulations, and production fragmentation are inducing changes in both networks and may \textcolor{black}{NGB} cause clustering patterns baier2007free,baltagi2008estimating,ivlevs2010fdi,behrens2012dual. Thus, it is not clear if FTAs and trade flows are unidirectional or bidirectional related in the \textcolor{black}{NGB} sense, and the empirical results obtained from panel analysis and impact evaluation methods on FTAs effectiveness at a dyad level show mixed evidence baier2004economic.

{\color{black}The DSBMM accommodates the existing challenges in estimating structural Gravity models. The latent block structure allows for heteroschedasticity ($\sigma_{qr}^{2}$) based on the pair block interacting and potential dependence between the probability of exporting and the volume of trade. Censoring in economic data has resulted in different approaches, such as sample selection, hurdle, and two--part models jones2000health. In particular, in international trade, helpman2008estimating assumes the zeros represent missing data, suggesting a Heckman procedure in which the selection and trade flows follow different mechanisms with unobserved correlated components. In the case of silva2006log, the zeros are “genuinely” a no--trade situation between pairs of countries, and a single mechanism is assumed, that is, the covariates and the effect are the same for the level of trade and the cases of absence of trade. In this setting, the application of DSBMM to trade is equivalent to silva2006log in assuming that zeros imply no trade. Still, it differs by considering a two--part specification with separate mechanisms for the excess of zeros and positive trade flows similar to the one in metulini2018spatial,egger2011trade. More importantly, in the DSBMM both mechanisms are not independent, they share a set of common unobserved factors captured in the latent components that induce clustering in the network, i.e. membership variables $Z_{it}$ and $Z_{jt}$, that can capture positive or negative (nonlinear) relationship between the probability of exporting, the trade flow and the variance of the trade flows. Indeed, contrary to standard two--part models, the two mechanisms in the DSBMM cannot be independently estimated.

The simplified specification of the Gravity in (ref) can also be easily extended to relax some assumptions and focus on other objectives in the trade literature such as estimating the parameter of interest for trade policy, i.e. $(1-\varrho)\lambda$ and $(t_{ijt})^{1-\varrho}$. For instance, more covariates can be included to approximate trade barriers, explanatory variables can be incorporated into the probability of trade, or alternative distributions for trade flows can be used. The assumption of constant unobserved components within a pair of interacting blocks can be relaxed by introducing node--time and pair fixed effects or random effects. This latter option constitutes a mixed--effects alternative approach to deal with the multi--resistance terms, which can be interpreted as clustering the fixed effects, while allowing for the inclusion of time--invariant and node--time invariant covariates such as policy variables. To exemplify the flexibility of the DSBMM, some of these extensions are considered in Section (ref) in Supplementary Material. Therefore, the application of DSBMM to the trade and FTA network provides a unified framework for two streams of literature: Gravity models and community detection.}

Data description

The dataset used in this application includes 159 countries detailed in Table (ref), for the period 1995--2017. The trade network corresponds to the value of all goods imported (CIF) between pairs of countries collected by the International Monetary Fund (IMF) in U.S. dollars, which are transformed into constant terms by using the GDP PPP deflator from the same source with base year 2010, as in helpman2008estimating. The information reported to the IMF is not always complete due to missing reporting countries (nodes) or specific flows (edges) at a period $t$. We rely on the imputation techniques suggested by IMF dippelsman2018new.

Regarding the FTAs, the network is constructed using the database on Economic Integration Agreements maintained by the NSF-Kellogg Institute. The information provided on the type of FTA is dichotomized into the existence or non--existence of any trade agreement, from a non-reciprocal preferential trade arrangement (NR--PTA) to an Economic Union (EUN). The former type of agreement makes the adjacency matrix of the FTAs non--symmetric. This data is collected from various sources, including the CIA World Factbook and the WTO database, among others. If there is no information on any of the sources for a pair of countries, the database includes FTAs based on trade flows. For instance, if two countries have never traded and there is no international agreement between them in the databases reviewed, it is assumed that they have not signed an FTA. Although these imputations covered 40% of the pairs between 1950--2012, the assumptions are reasonable and commonly accepted in the literature on FTA effectiveness baier2019widely. Finally, the covariates used for the first layer are from the dynamic gravity dataset of the U.S. International Trade Commission (USITC), which collects GDPs in real terms from the World Development Indicators (WDI), and the distance ($dist_{ij}$) refers to the population--weighted distance between each pair of countries gurevich2018dynamic.

A general description of the resulting multi-layer FTAs--Trade Flows is presented in Table (ref), which confirms some of the empirical stylized facts in international trade. Since the beginning of the period, both networks have been dense, and without considering the direction of the edges, there is only one component (WCC), indicating that a path always exists to connect any pair of countries. Suppose the direction is accounted (SCC). In that case, the FTA has more than one component due to the presence of regional agreements and non--reciprocal agreements between developed and developing countries, which makes some nodes unilaterally unreachable; i.e., there exists a path from $i$ to $j$, but no path in the opposite direction. Over time, the FTAs network converges into a single component. In the trade network, following WCC and SCC, there is only one component during the period under investigation. It can be seen that the trade network is significantly more connected than the FTAs.

{ {8pt}

}

Both layers share common trends in most of the indicators in Table (ref). The average degrees and densities tend to increase over time, which can be the first indication of dependence. Consistent with this higher connectivity, the average path length for the unweighted layer has decreased; however, this is not the case for the weighted one. Despite its general decline, the trend experienced a reversal in the 2000s. The number of triangles in the FTA and trade networks has steadily increased, resulting in higher clustering coefficients. Moreover, the betweenness centrality is decreasing for both networks, which can be interpreted as a loss of influence of intermediary nodes given the presence of new one--step shortest paths and an average increase in the connectivity of the nodes.

The different topologies of the two layers and the temporal changes in connectivity call for the use of DSBM. Moreover, the simultaneous changes in the connectivity features may suggest a significant influence of trade policies (FTAs) on the trade network. To find evidence of this possibility a DSBMM framework can be applied. The DSBMM offers the possibility of testing for the existence of community structure in trade flows after controlling for GDP and distance, which, to the best of our knowledge, has not been directly examined in the literature. Moreover, the DSBMM can reveal whether changes in the community structure of the (unobserved) multi--resistance terms are related to changes in trade policy through their membership dynamics.

Results

To apply the DSBMM in ((ref)) and ((ref)) to the FTAs--Trade flows multidimensional network, the number of blocks $Q^{(\ell)}$, $\ell=1,2$, needs to be selected. \textcolor{black}{Standard criteria provide a large number of communities since either they are not designed for DSBMM or they are not consistent (see Section (ref) in Supplementary Material for a discussion and further results). In addition, assuming time--varying number of communities can make inference computationally demanding. Thus, in this paper, we assume a constant number of communities and choose it following the results from previous studies. Some robustness checks of the main results are then discussed. We should note that assuming a constant number of blocks is not too restrictive, as some blocks can be empty in certain periods, which compensates for the reduced flexibility of the model. Our work can be extended by assuming a random number of communities and building on infinite Markov models dufays2016infinite, or mixtures of finite mixture models combined with the telescopic sampling method recently introduced by fruhwirth2021generalized.}

Regarding the trade flows, based on magerman2020pecking, a core--periphery structure is a parsimonious two-community representation of the network: a core, that is, a subset of highly connected nodes, and a periphery, the rest of the nodes that are more densely linked to the core than to the other nodes in the periphery. This can be interpreted using a pecking order scheme, firms in each country export following a destination list ordered according to market size and lower trade costs. Hence, countries with less productive firms (periphery) are only able to trade with the top of the list, while countries with more productive firms (core) can trade with large and smaller markets and/or less and more costly destinations. In order to verify the core--periphery structure, we assume $Q^{(1)}=2$. Regarding the FTAs layer, the literature provides little information on the number of communities. In this paper, following previous findings in barigozzi2011identifying, we select four communities in the second layer, i.e., $Q^{(2)}=4$. Some robustness checks are then considered assuming $Q^{(1)}=Q^{(2)}=2$ and $Q^{(1)}=Q^{(2)}=4$ (see Figure (ref) and (ref) in the Supplementary Material).

The results under $Q^{(1)}=2$ and $Q^{(2)}=4$ confirm a core--periphery structure in the Trade layer (\textcolor{black}{top--left panel of} Figure (ref)). Community 1 corresponds to the core with almost a fully connected sub--graph in terms of intra--block edges (value of $\nu_{11}^{(1)}$ close to 0.98 and is highly dense in the inter-block edges to community 2 (value of $\nu_{12}^{(1)}$ close to 0.84). Community 2 is the periphery since it is better connected to the core (value of $\nu_{21}^{(1)}$ close to 0.79) than to other nodes within its block (value of $\nu_{22}^{(1)}$ close to 0.36). \textcolor{black}{At the beginning of the sample, the community sizes are unbalanced and the core has a smaller membership probability compared to the periphery ($\hat{\alpha}^{(1)}_{1}=0.29$ and $\hat{\alpha}^{(1)}_{2}=0.71$, respectively). For a list of countries by block at the beginning and end of the period, see Figure (ref) and Table (ref) in the Supplementary Material}. Additionally, \textcolor{black}{the flows with highest variance expressed as $ \mathbb{E}(\nu^{(1)2}_{qr}\exp(2(\beta_{qr}^{(1)})'X_{ijt}+\sigma_{qr}^{(1)2})[\exp(\sigma^2)/\nu_{qr}-1]|Y,Z_{it}=q,Z_{jt}=r)$, are those from the periphery to the core (bottom--left panel of Figure (ref)) with a value close to 53.5 in log--scale. The latent membership variables in the DSBMM capture a non--linear and non--monotonic relationship between variance and the first moment of the trade flows, i.e. $ \mathbb{E}(\nu^{(1)}_{qr}\exp((\beta_{qr}^{(1)})'X_{ijt}+\sigma_{qr}^{(1)2}/2)|Y,Z_{it}=q,Z_{jt}=r)$. Although, within each community pair the relationship between mean and variance is positive, the variance in the intra--community trade flows of the core are lower than the variance in the flows from periphery to core. This result contradicts the models used in the literature that assume larger variance in larger trade flows silva2006log. In the top--left panel of Figure (ref), a positive relationship, between the probability to trade and first moment of the trade volume, $\mathbb{E}(\nu^{(1)}_{qr}\exp(X^{(1)}_{ijt}\beta_{qr}^{(2)}+\sigma_{qr}^{(1)2}/2)|Y>0,Z_{it}=q,Z_{jt}=r)$, is observed suggesting a close dependence between both mechanisms, as underlined by previous studies helpman2008estimating,egger2011trade.}.

In the FTA layer, the structure seems more complex, except for community 3, where intra--block connectivity $\nu_{qq}^{(2)}$ is higher than its inter-block counterpart, the other blocks have a non--assortative feature, that is community 1, 2 and 4 are highly connected to block 3 and less densely linked to the other communities (see right panel of Figure (ref)). This configuration is compatible with a core--periphery graph but with some particularities. \textcolor{black}{ The core is Community 3, which includes the countries with the largest number of treaties (high degree nodes) such as Western Europe and countries with a large number of indirect connections due to non--reciprocal agreements (highest eigen--centrality), such as the U.S. under the Generalized System of Preferences (GSP). At the beginning of the sample, the block sizes are unbalanced and two blocks, including the core, have a smaller initial membership probability ($\hat{\alpha}^{(2)}_{3}=0.14$ and $\hat{\alpha}^{(2)}_{4}=0.11$) compared to the other two blocks ($\hat{\alpha}^{(2)}_{1}=0.41$ and $\hat{\alpha}^{(2)}_{2}=0.34$). See Table (ref) for further details on the community composition. See Table (ref) for further details.} For instance, the intra--edges of community 4 are almost as dense as their in--degree from the core, which is the case of countries from eastern Africa (see Figure (ref) and (ref)).

The block structure in the trade flows is consistent with the hypothesis of the picking order destination list. In particular, the most productive firms in the periphery (and in the core), which are those able to export, prioritize destinations with large market sizes and lower trade costs (the core) magerman2020pecking. However, it is unclear from this theory why the flows from the periphery to the core are more unstable (higher variance). It could be due to supply factors such as productivity heterogeneity across countries in the periphery or demand factors from the core, since densely connected countries may have more alternative suppliers to consider.

Regarding the connectivity parameters in the Gravity equation, the results evidence significant differences between blocks, reflecting contrasting unobserved barriers and multi--resistance terms. \textcolor{black}{The lowest intercept is observed between the core countries, while the highest corresponds to the peripheral nodes (see Figure (ref))}. At the same time, this is compensated by a less negative distance parameter between the core countries and an almost twice as negative marginal effect within community 2. A similar pattern applies to the exporter and importer GDP elasticities, both are higher in the core. In the inter--block interaction, when the core (periphery) is the exporter, trade flows are more (less) responsive to the GDP, in contrast to its role as an importer. Hence, the heterogeneity in the Gravity equation parameters is coherent with the differences in the connectivity parameters of the unweighted part of the trade network ($\nu^{(1)}$) and the international trade literature studying intensive and extensive margins bernard2007firms. The gaps in the distance parameter can be attributed to differences in technology and/or transportation costs. In the case of the GDP elasticities, the differences suggest the possibility of persistent trade deficits for the periphery, which is not fully accounted for by the standard trade theory that assumes the elasticities should be close to one for all countries. \textcolor{black}{The interpretation of the GDP coefficients is not straightforward because trade flows are measured in gross value, while GDP is measured in value added anderson2004trade. Nevertheless, these discrepancies in income elasticity have raised discussions on the influence of home market effects and oligopolistic market structures feenstra2001using,fujiwara2024firm and they are also relevant in macroeconomic theory, in particular the literature on balance of payment constrained growth and structural change thirlwall2012balance,matsuyama2019engel and consistent with the findings in silva2006log, based on Poisson Gravity models.}

This network structure is not static, and countries may change membership through time and inherit the corresponding set of connectivity parameters, which to the extent of our review, is not taken into account in the standard partial or general equilibrium models used for ex-post and ex-ante trade policy evaluation. In this sense, for the trade network, 22.0% of countries changed membership at least one time between 1995--2017 mostly from the periphery to core, while in the FTA layer, all countries deviate from their initial community (see Figure (ref) and (ref)). This new result complements existing community detection studies on trade flows, which deduce an increase in network density by identifying blocks at each point in time. However, given their static approach, they are not able to recover specific node changes.

The block dynamic allows us to test for \textcolor{black}{NGB causality} and check if the changes in the membership structure of the FTAs and in the memberships of the unobserved components of the Trade network are related. In Figure (ref), the membership dependence between the two layers of the multidimensional network is summarized into four graphs, each of them associated with a vector $\kappa_{q}^{(\ell)}$. For instance, the first graph shows the transition probability of moving/staying in community 1 (core) of the trade network, and it can be noticed that its diagonal elements are significant, i.e. $W^{(1)}_{it-1,1}$ and the intercept, which denotes persistence in membership. More importantly, this persistence is attenuated by the influence of the country's membership in the FTA network. Essentially, a change in membership to the core of the FTA network induces a change in the trade network, potentially producing a decrease in the multi-resistance terms and unobserved barriers. The three graphs of the second layer describe the FTA network transition parameters and indicate that membership persistence is lower compared to the Trade network. From a variable selection approach, it appears not to be associated with the block structure of the latter. \textcolor{black}{In other words, there is evidence that FTA core membership $W_{2,3}$ influences the probability of moving to the Trade core (see Panel $\kappa^{1}_{1}$). The evidence of a predictive impact of Trade core membership $W_{1,1}$ on the transition to FTA core is weaker in all FTA communities (see panels $\kappa^{2}_{1}$,$\kappa^{2}_{2}$ and $\kappa^{2}_{3}$). Indeed, there are some cases of transitions to the Trade core (see boxes in Figure (ref)), which include countries in the FTA core from the beginning of the sample (e.g. HUN) and countries that entered the FTA core later in the sample (e.g. BLR, EST, LTU, LVA). See Figure (ref) and Table (ref) for more details. However, a proper test of the NGB requires estimating the log--Bayes factor among the competing hypotheses. Following the reverse logistic approach, the log--factor in favor of no NGB from Trade to FTA, from FTA to Trade, and no NGB causality is -6136.21,-8853.52, and -5424.32, respectively. This suggests that all models incorporating some form of NGB causality perform better than models with no NGB causality. In particular, the posterior probability of the model allowing bilateral NGB causality is significantly higher than unilateral alternatives or no dependence---i.e., close to one.} The same result applies for $Q^{(1)}=Q^{(2)}=4$, but for $Q^{(1)}=Q^{(2)}=2$ there is evidence of bidirectional causality (see Figure (ref) and (ref) in the Supplementary Material).

In summary, the bidirectional block--causality suggests that the FTAs have effectively modified the trade flows network structure, favoring the transition of countries from the periphery to the core and increasing global trade integration and vice--versa. In other words, the level of trade for a specific pair of countries also helps predict their probability of signing or not a FTA, through the clustering structure. This is consistent with the endogeneity argument suggested by the trade literature, where expected gains from signing a treaty affect the trade policies. This result is robust to the number of blocks.

The core--periphery structure is also robust to the number of blocks. However, some countries are not part of the core of the FTA network and are densely connected in the trade network (e.g. China and Costa Rica). This could be due to the type of exported products, technological particularities, and the level of integration to the Global Value Chains (GVC), among other factors that can be easily included in the DSBMM as extra covariates or as extra layers (e.g. the global investment network). Moreover, the time--varying membership captured by our DSBMM implies evidence of parameter heterogeneity in the Gravity equation, which may be used by the trade literature to evaluate the impact of changes in tariffs. This heterogeneity is especially relevant in the evaluation of trade policies aiming at a structural change of the network beyond a specific pair of countries in the same community, for example, the theoretical scenarios of full autarky or in ambitious treaties such as the Transatlantic Trade and Investment Partnership (TTIP) or the Trans--Pacific Partnership (TPP), in these cases the counterfactual scenario may also imply different parameters.

Conclusion

This work has focused on extending the DSBM to a multidimensional setting with dependent Hidden Markov Chains (HMC). This task is achieved assuming node membership is driven by a different HMC in each layer. The HMC dependence across layers is specified as a saturated multinomial logit model. The resulting model, DSBMM, enables us to study the relationship between the block structure of layers and provides a framework for a \textcolor{black}{NGB causality} test. Compared to the literature on multidimensional SBM, our DSBMM is novel in many aspects. First, it can capture coupling trends between the membership structure of the layers, decoupling tendencies, or more complex relationships. Second, it provides information on the sign, magnitude, and direction of the dependence. Moreover, node- and dyad-specific covariates are included to capture heterogeneity that may affect the weights and the sparsity patterns of the networks. Finally, the dependence between the four types of layers --- weighted and unweighted, directed and undirected --- can be easily incorporated into our model.

To cope with over-paramerization, we propose a new Bayesian shrinkage approach for DSBMM based on Multi--Laplacian prior distributions. This prior class allows for different degrees of shrinkage between and within groups of \textcolor{black}{lagged HMC variables} and accounts for the correlation structure of the \textcolor{black}{hidden state variables}. A Pólya Gamma representation makes the posterior distribution more tractable, and an efficient Gibbs sampler is designed to approximate the posterior distribution. The efficiency of the sampler, the effectiveness of the Multi-Laplace prior in recovering the true {\color{black}NGB} causal structure, and the superior fitting ability of the DSBMM compared to benchmark models are studied through simulation experiments under different network configurations and parameter settings.

We apply our DSBMM model to international trade networks and contribute to the debate on the effectiveness of FTAs beyond the country pairs level. Moreover, to the best of our knowledge, this is the first econometric model to integrate two commonly used approaches in international trade: community detection and gravity models. A temporal, multidimensional network of Trade Flows and FTAs is considered among 159 countries over 23 years. We found evidence of core--periphery structure in FTAs and trade networks, and of dyad and time heterogeneity of the Gravity equation parameters. The standard theoretical models do not account for this result, which may impact existing methods for FTA evaluations. Regarding the debate on the predictive power of FTAs, we find evidence of a \textcolor{black}{NGB} causal relationship between FTAs and trade. Moreover, the direction is sensitive to the choice of the number of communities. There is bidirectional block--causality between the FTAs and the unobserved barriers.

Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.