# Graphical Assistant Grouped Network Autoregression Model: a Bayesian
Nonparametric Recourse

Yimeng Ren1, Xuening Zhu2, Guanyu Hu3
1School of Statistics, Renmin University of China
2School of Data Science, Fudan University
3Department of Statistics, University of Missouri Columbia  Yimeng Ren and
Xuening Zhu are supported by the National Natural Science Foundation of China
(nos. 11901105, 71991472, U1811461), the Shanghai Sailing Program for Youth
Science and Technology Excellence (19YF1402700), and funds from Fudan
University. Xuening Zhu is the corresponding author.

###### Abstract

Vector autoregression model is ubiquitous in classical time series data
analysis. With the rapid advance of social network sites, time series data
over latent graph is becoming increasingly popular. In this paper, we develop
a novel Bayesian grouped network autoregression model to simultaneously
estimate group information (number of groups and group configurations) and
group-wise parameters. Specifically, a graphically assisted Chinese restaurant
process is incorporated under framework of the network autoregression model to
improve the statistical inference performance. An efficient Markov chain Monte
Carlo sampling algorithm is used to sample from the posterior distribution.
Extensive studies are conducted to evaluate the finite sample performance of
our proposed methodology. Additionally, we analyze two real datasets as
illustrations of the effectiveness of our approach.

Keywords: Graphical Assistant Chinese Restaurant Process, MCMC, Mixture
Models, Network Data, Vector Autoregression.

##  1 Introduction

The network data is becoming increasingly popular, which has many important
applications in various disciplines, including sociology (Simmel, 1950;
Wasserman et al., 1994; Hanneman and Riddle, 2005), genomics (Wu et al., 2010;
Horvath, 2011), psychology (Borgatti et al., 2009; Cramer et al., 2010;
Borsboom and Cramer, 2013), finance and economics (Zou et al., 2017; Zhu et
al., 2019; Leung et al., 2017). In the field of sociology, the network data is
used to study the structure of relationships by linking individuals or other
social units and by modeling interdependencies in behaviors or attitudes
related to configurations of social relations (O’malley and Marsden, 2008). In
biological and genetic areas, the presence of an interaction between two genes
or proteins indicates a biologically functional relationship (Von Mering et
al., 2002), which makes the network data widely used in cancer and other
disease data analysis with promising results (Guda et al., 2009; Chuang et
al., 2007). A continuous response observed on each node at equally spaced time
points is very common in various applications such as biological studies,
economics, and finance. Building a network based model to recover the dynamics
of the responses is necessary.

We consider a network with $N$ nodes, which are indexed as $i=1,\cdots,N$. To
represent the network relationship among the network nodes, we employ an
adjacency matrix $\mathbf{A}=(a_{ij})\in\mathbb{R}^{N\times N}$, where
$a_{ij}=1$ indicates the $i$th node follows the $j$th node, otherwise
$a_{ij}=0$. Following the convention we let $a_{ii}=0$ for $1\leq i\leq N$.
Denote $Y_{it}$ as the continuous response collected from the node $i$ during
$1\leq t\leq T$. Correspondingly, we collect a number of nodal covariates
(e.g., user’s age, gender), which is denoted by $V_{i}\in\mathbb{R}^{p}$. To
model the dynamics of the response $Y_{it}$, Zhu et al. (2017) proposed a
network vector autoregression (NAR) model, which can be expressed as follows,

$\displaystyle
Y_{it}=\beta_{0}+\beta_{1}n_{i}^{-1}\sum_{j}a_{ij}Y_{j(t-1)}+\beta_{2}Y_{i(t-1)}+V_{i}^{\top}\gamma+\varepsilon_{it},$
(1.1)

where $n_{i}=\sum_{j}a_{ij}$ is the out-degree of node $i$ and
$\varepsilon_{it}$ is the random noise term. The NAR model as well as its
extensions have been implemented to a broad range of fields, including spatial
data modeling (Lee and Yu, 2009; Shi and Lee, 2017), social studies
(Sojourner, 2013; Liu et al., 2017; Zhu and Pan, 2018), financial risk
management (Härdle et al., 2016; Zou et al., 2017), and many others.

Despite the usefulness of the NAR model, it lacks flexibility in its model
formulation and suffers from model misspecification. Specifically, it
specifies a homogeneous network autoregression coefficient (i.e., $\beta_{1}$)
for all network nodes, which is highly restrictive in practice. To enhance the
model’s flexibility, a popular way is to consider a group structure of model
coefficients (Ke et al., 2015; Su et al., 2016; Ando and Bai, 2016;
Guðmundsson and Brownlees, 2021). With respect to the network data, Zhu and
Pan (2018) considered a grouped network autoregression (GNAR) model. To be
more specific, the nodes are clustered into several groups, and the nodes
within the same group share a set of common parameters. Both the group
memberships and group-specific regression coefficients are to be estimated.

Although the model’s flexibility is enhanced and the misspecification risk is
reduced, there are several unsolved challenges for the GNAR model. First, it
is necessary to incorporate the graphical information (i.e., $\mathbf{A}$)
into grouping process. Nevertheless, most existing methods ignore graphical
information in their grouping process. On the other hand, another line of
research takes advantages of the graphical information and uses nonparametric
Bayesian methods for clustering nodes (Sewell and Chen, 2017; Geng et al.,
2019; van der Pas and van der Vaart, 2018), but they totally ignore the
dynamics of $\\{Y_{it}\\}$. The graphical information, a connectivity pattern
among different nodes, will partially determines the group configurations. In
other words, the closely connected nodes should have higher probability
belonging to the same group. How to incorporate the graphical information into
the modelling procedure of the GNAR model is a challenging task. Second,
proposing a simultaneously inference procedure on both group configurations
and group-wised parameters is another important consideration for GNAR model.
The uncertainty of parameter estimation is often neglected in existing
literature, and it is difficult to get the probabilistic interpretations from
existing frequentist approaches. Lastly, another important problem in grouping
methods is how to determine the number of groups. Most existing methods
require specification of the number of groups first and then estimate group
configurations. The number of groups is determined by information criteria or
some ad hoc procedures. These procedures may ignore uncertainty of the
estimation of the number of groups in the first stage and lack uncertainty
quantification of the number of groups, which is also an important issue for
GNAR model.

Our methodology development is directly motivated by solving the
aforementioned challenges to fill the gap between nonparametric Bayesian
methods and network autoregression model. In this article, we propose a novel
graphical assistant grouped network autoregression (GAGNAR) model, considering
both the heterogeneous network effects across nodes and the graph information
in the grouping procedure under a novel nonparametric Bayesian mixture model
framework. Specifically, a graphically assistant Chinese Restaurant Process
(gaCRP) borrowed from Geng and Hu (2021) is introduced to the GNAR model
framework, which incorporates the graphical information and utilizes the
sampling procedure from the Chinese Restaurant Process (CRP) (Pitman, 1995).
The proposed prior will provide simultaneous inference on the number of groups
and group configurations. Moreover, the gaCRP has a Pólya urn scheme similar
to the traditional CRP. Due to this merit, we can develop a Gibbs sampler that
enables efficient full Bayesian inference on the number of groups, mixture
probabilities as well as group-specific parameters. We demonstrate its
excellent numerical performance through simulations and analysis of two real
datasets, compared with several benchmark methods.

Our proposed method is unique in the following aspects. First, the proposed
method provides a useful model-based solution for GNAR model that is able to
leverage graphical information without pre-specifying the number of groups.
Secondly, by adopting a full Bayesian framework, the grouping results yield
useful probabilistic interpretation. Thirdly, the developed posterior sampling
scheme also renders efficient computation and convenient inference since our
proposed collapsed Gibbs sampling algorithm by marginalizing over the number
of groups avoid complicated reversible jump Markov Chain Monte Carlo (MCMC)
algorithm (Green, 1995) or allocation samplers.

The rest of this paper is organized as follows. In Section 2, we briefly
review the NAR and GNAR models, followed by a review of Bayesian mixture model
and Dirichlet process. In addition, we describe the graphical assistant prior
model and introduce our proposed GAGNAR model and its theoretical properties.
The Bayesian inferences including the MCMC sampling algorithm, the post-MCMC
estimation and the model selection criterion for tuning parameter are
presented in Section 3. In Section 4, we study the finite sample performance
of our proposed method via extensive simulation studies. Two real datasets are
illustrated by our proposed method in Section 5. Section 6 concludes and
raises some in-depth discussions. All technique proofs and additional
numerical results are given in supplementary materials.

##  2 Methodology

###  2.1 Grouped Network Autoregression Model

Consider a network with $N$ nodes with the corresponding adjacency matrix
denoted by $\mathbf{A}$. Assume the nodes are clustered into $K$ latent
groups. Denote $\mathcal{F}_{t}$ as the $\sigma$-field generated by
$\\{Y_{it}:1\leq i\leq N,1\leq t\leq T\\}$. Given $\mathcal{F}_{t-1}$, the
observations at time point $t$ are assumed to be independent and follow a
Gaussian mixture distribution with $K$ components. Given the number of groups
$K$ and the adjacency matrix $\mathbf{A}$, the GNAR model can be expressed as
a mixture model, i.e.,

$\begin{split}&p(Y_{it}|\boldsymbol{\theta},\sigma^{2})=\sum_{k}\pi_{k}\mathcal{N}(\mu_{kt},\sigma_{k}^{2}),\\\
&\mathcal{N}(\mu_{kt},\sigma_{k}^{2})=\frac{1}{\sqrt{2\pi}\sigma_{k}}\exp\Big{\\{}-\frac{(Y_{it}-\mu_{kt})^{2}}{2\sigma_{k}^{2}}\Big{\\}},\\\
&\mu_{kt}=\beta_{0k}+\beta_{1k}n_{i}^{-1}\sum_{j}a_{ij}Y_{j(t-1)}+\beta_{2k}Y_{i(t-1)}+V_{i}^{\top}\gamma_{k},\end{split}$
(2.1)

where $\mathcal{N}(\mu_{kt},\sigma_{k}^{2})$ is the density function of normal
distribution with mean $\mu_{kt}$ and variance $\sigma_{k}^{2}$ in the $k$th
group, and $\pi_{k}$ is the mixing weight satisfying $\sum_{k}\pi_{k}=1$. Here
$\boldsymbol{\theta}_{k}=(\beta_{0k},\beta_{1k},\beta_{2k},\gamma_{k}^{\top})^{\top}$
collects the unknown parameters.

The GNAR model (2.1) characterizes the conditional mean of the response
$Y_{it}$ by a combination of four components. The first one is the baseline
effect $\beta_{0k}$, which is a constant but varying among different groups.
The second one is the network component
($\beta_{1k}n_{i}^{-1}\sum_{j}a_{ij}Y_{j(t-1)}$), which reflects the average
impact of following nodes at the previous time point. Therefore, we refer to
$\beta_{1k}$ as the group network effect. The third one is the momentum
component ($\beta_{2k}Y_{i(t-1)}$). It characterizes how the focal node is
influenced by its historical behaviors. The corresponding parameter
$\beta_{2k}$ is then referred to as group momentum effect, quantifying the
node’s self-dependence. The last one is the nodal covariate component
($V_{i}^{\top}\gamma_{k}$), which includes the node specific characteristics
and it is invariant over time. As implied by the GNAR model (2.1), the nodes
within the same group share the same set of coefficients, therefore they
exhibit similar dynamic patterns. Compared to the NAR model (Zhu et al.,
2017), the GNAR model is able to capture the nodes’ heterogeneous pattern in
the group level. To estimate the model parameters simultaneously with the
group memberships, an EM algorithm could be designed (Zhu and Pan, 2018).

Despite the usefulness, the model has three main limitations. First, the
estimation can be unreliable when the observed time length is short. That
decreases the model’s stability and robustness in practice. Second, the
network structure information is not fully exploited for estimating the group
memberships. Third, the number of groups $K$ cannot be automatically estimated
but needs to be specified in advance. That increases the possibility that the
model is misspecified. As a consequence, we are motivated to consider the
problem in a Bayesian framework and introduce a graph assisted prior, where
the network structure information is employed when forming the groups. The
details are presented in the following section.

###  2.2 Bayesian Graph Assisted Prior

For each node $i$, suppose it carries a latent group membership
$z_{i}\in\\{1,\cdots,K\\}$, where $z_{i}=k$ implies that the node $i$ is from
the $k$th group. Given $\\{z_{i}\\}$ and
$\\{\boldsymbol{\theta}_{k},\sigma_{k}^{2}:1\leq k\leq K\\}$, we have a
Gaussian mixture model (Marin et al., 2005) as follows,

$\displaystyle(Y_{it}|z_{i},\boldsymbol{\theta}_{z_{i}},\sigma_{z_{i}}^{2})\sim\mathcal{N}(\mu_{z_{i}t},\sigma_{z_{i}}^{2}),\quad
i=1,\ldots,N,\quad t=1,\ldots,T,$

where
$\mu_{z_{i}t}=\beta_{0z_{i}}+\beta_{1z_{i}}n_{i}^{-1}\sum_{j}a_{ij}Y_{j(t-1)}+\beta_{2z_{i}}Y_{i(t-1)}+V_{i}^{\top}\gamma_{z_{i}}$.
Under the Bayesian hierarchical model, for the unknown group memberships, we
specify the model priors as

$\begin{split}&z_{i}\sim\text{Cat}(\pi_{1},\ldots,\pi_{K}),\quad
i=1,\ldots,N,\\\
&\pi_{1},\ldots,\pi_{K}\sim\text{Dirichlet}_{K}(\alpha,\ldots,\alpha),\\\
&\boldsymbol{\theta}_{j}\sim\pi(\boldsymbol{\theta}),\quad
j=1,\ldots,K,\end{split}$ (2.2)

where $z_{i}$ follows a categorical distribution with parameters
$\pi=\\{\pi_{1},\dots\pi_{K}\\}$, $\pi$ follows a Dirichlet distribution with
parameter $\alpha$, and $\pi(\boldsymbol{\theta})$ is joint prior for all
group-wise parameters. When $K$ is unknown, we can let $K$ go to infinity, and
then the model in (2.2) becomes the Dirichlet process mixture model (DPMM;
Ishwaran and Zarepour, 2002). The DPMM is very flexible and efficient for
density estimation. However, it suffers from several problems in estimating
the group structure. First, the partitions sampled from the DPMM posterior
tend to have multiple small transient groups, which decreases the
interpretability of the model. Second, the DPMM is inconsistent on the
estimation of the number of groups (Miller and Harrison, 2018). Third, the
DPMM does not take advantage of the network structure information: the
connected nodes do not have higher probability being in the same group.

In order to address the above problems, we employ a Bayesian graph assisted
prior. Specifically we assume that the group membership
$z_{i}\in\\{1,\cdots,K\\}$ is generated from the weighted Dirichlet process.
The process assumes that the nodes’ memberships are sequentially generated.
Suppose the nodes $1\leq j\leq i-1$ form $K^{\prime}$ groups, then the
conditional distribution for $z_{i}$ is expressed as

$\displaystyle P(z_{i}=k\mid
z_{1},\ldots,z_{i-1})\propto\Bigg{\\{}\begin{array}[]{ll}p_{k}(\mathbf{W}),&1\leq
k\leq K^{\prime}\\\ \alpha,&k=K^{\prime}+1\end{array}$ (2.5)

where $p_{k}(\mathbf{W})=\sum_{j=1}^{i-1}w_{ij}I(z_{j}=k)$ and
$\mathbf{W}=(w_{ij})$ is a weighting matrix calculated based on $\mathbf{A}$.
The $w_{ij}$ in the weight matrix denotes the connection strength between node
$i$ and $j$. The parameter $\alpha$ controls the probability for introducing a
new group. As a result, (2.5) implies that a node is more likely to fit in the
same group with its connected friends.

###### Remark 1

With the probabilistic Bayesian framework, simultaneous estimation of the
number of groups and the grouping configuration is typically achieved by
complicated searching algorithms (e.g., reversible jump MCMC (Green, 1995)).
However, it suffers from high computational complexity. In our framework, the
specification of (2.5) allows the nodes’ memberships to be generated
sequentially and therefore avoids setting the number of groups in advance.
According to (2.5), a new group will formulate once $z_{i}$ is assigned with a
new group label $K^{\prime}+1$. This enables a sequential generation of
groups. Once we generate all the group memberships $\\{z_{i}\\}$, we know
immediately the number of groups $K$.

Let $G_{0}$ be a continuous probability measure on
$\mathbb{R}^{p}\times\mathbb{R}^{+}$. We define the full conditional
distribution of $\boldsymbol{\theta}_{z_{i}}$ given
$\boldsymbol{\theta}_{z_{1}},\cdots,\boldsymbol{\theta}_{z_{i-1}}$ as

$\displaystyle
f(\boldsymbol{\theta}_{z_{i}}\mid\boldsymbol{\theta}_{z_{1}},\ldots,\boldsymbol{\theta}_{z_{i-1}})\propto\sum_{r=1}^{K^{*}}\sum_{j=1}^{i-1}w_{ij}I(\boldsymbol{\theta}_{z_{i}}=\boldsymbol{\theta}^{*}_{r})f_{\boldsymbol{\theta}_{r}^{*}}(\boldsymbol{\theta}_{z_{i}})+\alpha
G_{0}(\boldsymbol{\theta}_{z_{i}}),$ (2.6)

where $f(\cdot)$ is the density function, $K^{*}$ denotes the number of groups
excluding the $i$-th observation,
$\boldsymbol{\theta}^{*}_{1},\ldots,\boldsymbol{\theta}^{*}_{K^{*}}$ are
$K^{*}$ distinguished values of
$\boldsymbol{\theta}_{z_{1}},\ldots,\boldsymbol{\theta}_{z_{i-1}}$ and
$I(\cdot)$ is the identity function. In practice we set the base distribution
$G_{0}$ to be the normal inverse gamma distribution.

Inspired by the graph assisted Chinese restaurant process (gaCRP) (Geng and
Hu, 2021) for survival model, we define the graph distance $d_{ij}$ between
node $i$ and $j$ as the shortest path length between node $i$ and $j$. If node
$i$ and $j$ cannot be connected with a finite length of path, then
$d_{ij}=\infty$. The weight is then defined as

$\displaystyle w_{ij}=\Bigg{\\{}\begin{array}[]{ll}1,&\text{ if }d_{ij}\leq
1\\\ \exp(-d_{ij}\times h),&\text{ if }d_{ij}>1\end{array}\Bigg{.}.$ (2.9)

Therefore the weight $w_{ij}$ characterizes the closeness of the node $i$ and
$j$ in the network. The parameter $h$ denotes the smoothing scale of the
weighting matrix. A larger value of $h$ indicates the weight is more
concentrating on the local scale. For simplicity, we refer to gaCRP introduced
above as gaCRP($\alpha,h$). Although the idea of using a weighted version of
CRP has been explored in Geng and Hu (2021), their model is based on Cox
regression for survival data analysis, which significantly differs from our
model, since it lacks discussion on dynamic process of $Y_{it}$.

###### Remark 2

Implied by (2.9), if $h=0$, then the gaCRP is the same as the traditional CRP
and may tend to over cluster the nodes. If $h\rightarrow\infty$, the
probability for a new node to choose the existing groups tends to be small,
which may also lead to the over-clustering problem. Hence, appropriate
selection of $h$ is critical for the model performance. We discuss in details
about the selection of $h$ in Section 3.3.

In summary, we express the proposed model hierarchically as

$\begin{split}&\text{Data
model}:Y_{it}|z_{i},\boldsymbol{\theta}_{z_{i}},\sigma^{2}_{z_{i}}\sim\mathcal{N}(\mu_{z_{i}t},\sigma_{z_{i}}^{2}),\\\
&\text{Process
model}:p(z_{i}|z_{1},\cdots,z_{i-1})\sim\text{gaCRP}(\alpha,h),\\\
&\text{Parameter
model}:(\boldsymbol{\theta}_{z_{i}},\sigma_{z_{i}}^{2})\sim\text{NIG}(\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0}),\end{split}$
(2.10)

where $\boldsymbol{\tau}_{0}$, $\boldsymbol{\Sigma}_{0}$, $a_{0}$ and $b_{0}$
are hyper-parameters for base distribution of
$(\boldsymbol{\theta}_{k},\sigma_{k}^{2})$. Here the prior distribution of
$(\boldsymbol{\theta}_{z_{i}},\sigma_{z_{i}}^{2})$ is set to be normal inverse
gamma (NIG) distribution with parameters
$(\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0})$. Without loss of
generality, we choose $\boldsymbol{\tau}_{0}=\mathbf{0}$,
$\boldsymbol{\Sigma}_{0}=100\mathbf{I}$, $a_{0}=0.01$, $b_{0}=0.01$ in our
simulation study and real data applications. We refer to the model (2.10) as
graph assisted group network autoregression (GAGNAR) model.

##  3 Bayesian Estimation

In this section, we discuss the model estimation for the GAGNAR model. Define
group membership vector
$\mathbb{Z}=(z_{1},\cdots,z_{N})^{\top}\in\mathbb{R}^{N}$. If the group
memberships are known, then the model parameters in the GNAR model (2.1) can
be easily estimated with least square estimation method. However, in practice
$\mathbb{Z}$ is unobserved thus the least square estimation is infeasible. To
estimate the model parameters simultaneously with the latent group
memberships, we employ the Markov chain Monte Carlo (MCMC) algorithm with
details given as follows.

###  3.1 MCMC Algorithm

To conduct Bayesian estimation, we first derive the posterior distribution of
the unknown parameters and group memberships. Define
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{1},\cdots,\boldsymbol{\theta}_{K})$
and $\boldsymbol{\sigma}^{2}=(\sigma_{1}^{2},\cdots,\sigma_{K}^{2})$. In
addition, denote $\mathbb{V}=(V_{1},\cdots,V_{N})^{\top}\in\mathbb{R}^{N\times
p}$, $\mathbb{Y}_{t}=(Y_{1t},\cdots,Y_{Nt})^{\top}\in\mathbb{R}^{N}$, and
$\mathbb{Y}=(\mathbb{Y}_{1},\cdots,\mathbb{Y}_{T})\in\mathbb{R}^{N\times T}$.
Given the data $\\{\mathbb{Y},\mathbb{V}\\}$, the joint posterior distribution
of latent membership $\mathbb{Z}$ and unknown parameters
$\\{\boldsymbol{\theta},\boldsymbol{\sigma}^{2}\\}$ could be written as

$\displaystyle\pi(\mathbb{Z},\boldsymbol{\theta},\boldsymbol{\sigma}^{2}|\mathbb{Y},\mathbb{V})\propto\pi(\mathbb{Z})\pi(\boldsymbol{\theta})\pi(\boldsymbol{\sigma}^{2})f(\mathbb{Y}|\mathbb{Z},\mathbb{V},\boldsymbol{\theta},\boldsymbol{\sigma}^{2}),$
(3.1)

where $\propto$ denotes “proportional to", and $\pi(\mathbb{Z})$,
$\pi(\boldsymbol{\theta})$, $\pi(\boldsymbol{\sigma}^{2})$ stand for the prior
distributions for $\mathbb{Z}$, $\boldsymbol{\theta}$, and
$\boldsymbol{\sigma}^{2}$ respectively. For convenience, let
$\mathbb{X}_{t}={(X_{1t},\cdots,X_{Nt})^{\top}}=({\bf
1}_{N},\mathbf{W}\mathbb{Y}_{t-1},\mathbb{Y}_{t-1},{\bf
1}_{N}V_{i}^{\top})\in\mathbb{R}^{N\times(p+3)}$, then given
$\boldsymbol{\theta}_{k},\sigma_{k}^{2}$ we have

$\displaystyle(\mathbb{Y}_{t}|\boldsymbol{\theta}_{k},\sigma_{k}^{2})\sim\mathcal{N}(\boldsymbol{\mu}_{k},\sigma_{k}^{2}\mathbf{I}),\text{
where }\boldsymbol{\mu}_{k}=\mathbb{X}_{t}\boldsymbol{\theta}_{k}.$

The prior distribution for $\mathbb{Z}$ is specified as gaCRP$(\alpha,h)$ as
stated in Section 2.2. In addition, we specify the prior distribution for
$\boldsymbol{\theta}_{k}$ as multivariate normal distribution
$\mathcal{N}(\boldsymbol{\tau}_{0},\sigma_{k}^{2}\boldsymbol{\Sigma}_{0})$ and
we set the prior for $\sigma_{k}^{2}$ as the commonly used scaled inverse
gamma distribution $\text{IG}(a_{0},b_{0})$, where $\boldsymbol{\tau}_{0}$,
$\boldsymbol{\Sigma}_{0}$, $a_{0}$ and $b_{0}$ are the hyper parameters.

Given the joint posterior distribution of
$\\{\mathbb{Z},\boldsymbol{\theta},\boldsymbol{\sigma}^{2}\\}$ in (3.1), we
can obtain the conditional posterior distribution for each latent group
membership and model parameters respectively. This enables us to conduct the
MCMC algorithm for model estimation. We first give the conditional posterior
distribution for $z_{i}$ $(i=1,\cdots,N)$ as follows. For convenience we let
$\mathbb{Z}_{(-i)}=\\{z_{j}:j\neq i\\}$. Denote full conditional of $Z_{i}$ as
$\widetilde{f}(z_{i})=f(z_{i}|\mathbb{Z}_{(-i)},\mathbb{Y},\mathbb{V},\boldsymbol{\theta},\boldsymbol{\sigma}^{2})$.
Further define
$\mathbb{Y}_{(i)}=(Y_{i1},\cdots,Y_{iT})^{\top}\in\mathbb{R}^{T}$ and
$\mathbb{X}_{(i)}=(X_{i1},\cdots,X_{iT})^{\top}\in\mathbb{R}^{T\times(p+3)}$.
Then we have

$\displaystyle\widetilde{f}(z_{i})$ $\displaystyle\propto
f(\mathbb{Z},\boldsymbol{\theta},\boldsymbol{\sigma}^{2}|\mathbb{Y},\mathbb{V})\propto
f(\mathbb{Y}_{(i)}|z_{i},\mathbb{V},\boldsymbol{\theta},\boldsymbol{\sigma}^{2})f(z_{i}|\mathbb{Z}_{(-i)}).$

The analytical form of $\widetilde{f}(z_{i})$ is given in Proposition 1.

###### Proposition 1

Suppose currently the nodes are formed into $K^{\prime}$ groups. Then
according to (2.10) we have

$\displaystyle\widetilde{f}(z_{i})\propto\Bigg{\\{}\begin{array}[]{ll}\kappa_{k}f(\mathbb{Y}_{(i)};\mathbb{X}_{(i)},\boldsymbol{\theta}_{k},\sigma_{k}^{2}),&\text{
for }1\leq k\leq K^{\prime}\\\ \alpha
g(\mathbb{Y}_{(i)};\boldsymbol{\tau}_{0},\boldsymbol{\tau},\boldsymbol{\Sigma}_{0},\boldsymbol{\Sigma},a_{0},b_{0}),&\text{
for }k=K^{\prime}+1\end{array},$

where $\kappa_{k}=\sum_{j\neq i}w_{ij}I(z_{j}=k)$ and

$\displaystyle
f(\mathbb{Y}_{(i)};\mathbb{X}_{(i)},\boldsymbol{\theta}_{k},\sigma_{k}^{2})=(2\pi\sigma_{k}^{2})^{-(T-1)/2}\exp\Big{\\{}-\frac{1}{2\sigma_{k}^{2}}\|\mathbb{Y}_{(i)}-\mathbb{X}_{(i)}\boldsymbol{\theta}_{k}\|^{2}\Big{\\}},$
$\displaystyle
g(\mathbb{Y}_{(i)};\boldsymbol{\tau}_{0},\boldsymbol{\tau},\boldsymbol{\Sigma}_{0},\boldsymbol{\Sigma},a_{0},b_{0})=\varphi(\boldsymbol{\Sigma}_{0},\boldsymbol{\Sigma},a_{0},b_{0})\big{\\{}b_{0}+\frac{1}{2}(\boldsymbol{\tau}_{0}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\tau}_{0}+\mathbb{Y}_{(i)}^{\top}\mathbb{Y}_{(i)}-\boldsymbol{\tau}^{\top}{\boldsymbol{\Sigma}}^{-1}\boldsymbol{\tau})\big{\\}}^{-(a_{0}+\frac{T-1}{2})},$
$\displaystyle\varphi(\boldsymbol{\Sigma}_{0},\boldsymbol{\Sigma},a_{0},b_{0})=\frac{b_{0}^{a_{0}}\Gamma(a_{0}+\frac{T-1}{2})|{\boldsymbol{\Sigma}}|^{1/2}}{(2\pi)^{(T-1)/2}\Gamma(a_{0})|\boldsymbol{\Sigma}_{0}|^{1/2}},$
$\displaystyle{\boldsymbol{\Sigma}}=(\boldsymbol{\Sigma}_{0}^{-1}+\mathbb{X}_{(i)}^{\top}\mathbb{X}_{(i)})^{-1},~{}~{}~{}\boldsymbol{\tau}=\boldsymbol{\Sigma}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\tau}_{0}+\mathbb{X}_{(i)}^{\top}\mathbb{Y}_{(i)}).$

By the formulation of $\widetilde{f}(z_{i})$, we can observe that the
probability $z_{i}$ belonging to an existing group $k$ ($1\leq k\leq
K^{\prime}$) is closely related to the network weighting matrix $\mathbf{W}$.
Particularly $\kappa_{k}$ denotes a weighted average ratio of its connected
friends belonging to the group $k$, which characterizes a stickiness to the
group $k$ of $i$th neighbourhood. If the stickiness level is higher, then the
conditional probability of $z_{i}=k$ will be higher. That is how the network
topology information is involved in the sampling procedure. Next, the node is
allowed to fit into a new group with probability proportional to
$g(\mathbb{Y}_{(i)};\boldsymbol{\tau}_{0},\boldsymbol{\tau},\boldsymbol{\Sigma}_{0},\boldsymbol{\Sigma},a_{0},b_{0})$,
which relates to both prior information and data information.

Next, we investigate the full conditional distribution for model parameters
$\\{\boldsymbol{\theta},\boldsymbol{\sigma}^{2}\\}$. Denote the posterior
distribution of $(\boldsymbol{\theta}_{k},\sigma_{k}^{2})$ as
$\widetilde{f}(\boldsymbol{\theta}_{k},\sigma_{k}^{2})=f(\boldsymbol{\theta}_{k},\sigma_{k}^{2}|\mathbb{Z},\mathbb{Y},\mathbb{V})$,
then we derive $\widetilde{f}(\boldsymbol{\theta}_{k},\sigma_{k}^{2})$ in the
following Proposition 2.

###### Proposition 2

The full conditional distribution
$\widetilde{f}(\boldsymbol{\theta}_{k},\sigma_{k}^{2})(k=1,2,\cdots,K)$ is
given as

$\displaystyle\widetilde{f}(\boldsymbol{\theta}_{k},\sigma_{k}^{2})$
$\displaystyle\propto\pi(\boldsymbol{\theta}_{k},\sigma_{k}^{2};\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0})\prod_{z_{i}=k}f(\mathbb{Y}_{(i)};\mathbb{X}_{(i)},\boldsymbol{\theta}_{k},\sigma_{k}^{2})$
$\displaystyle\propto
f(\boldsymbol{\theta}_{k},\sigma_{k}^{2};\boldsymbol{\tau}^{*},\boldsymbol{\Sigma}^{*},a^{*},b^{*}),$

where

$\displaystyle
f(\boldsymbol{\theta}_{k},\sigma_{k}^{2};\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0})$
$\displaystyle=\frac{b_{0}^{a_{0}}}{(2\pi)^{d/2}|\boldsymbol{\Sigma}_{0}|^{1/2}\Gamma(a_{0})}\Big{(}\frac{1}{\sigma_{k}^{2}}\Big{)}^{a_{0}+\frac{d}{2}+1}\exp\Big{[}-\frac{1}{\sigma_{k}^{2}}\Big{\\{}b_{0}+\frac{1}{2}(\boldsymbol{\theta}_{k}-\boldsymbol{\tau}_{0})^{\top}\boldsymbol{\Sigma}_{0}^{-1}(\boldsymbol{\theta}_{k}-\boldsymbol{\tau}_{0})\Big{\\}}\Big{]},$
$\displaystyle\boldsymbol{\tau}^{*}=(\boldsymbol{\Sigma}_{0}^{-1}+\sum_{z_{i}=k}\mathbb{X}_{(i)}^{\top}\mathbb{X}_{(i)})^{-1}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\tau}_{0}+\sum_{z_{i}=k}\mathbb{X}_{(i)}^{\top}\mathbb{Y}_{(i)}),~{}~{}~{}\boldsymbol{\Sigma}^{*}=(\boldsymbol{\Sigma}_{0}^{-1}+\sum_{z_{i}=k}\mathbb{X}_{(i)}^{\top}\mathbb{X}_{(i)})^{-1},$
$\displaystyle
a^{*}=a_{0}+(T-1)N_{k}/2,~{}~{}~{}b^{*}=b_{0}+\frac{1}{2}\Big{[}\boldsymbol{\tau}_{0}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\tau}_{0}+\sum_{z_{i}=k}\mathbb{Y}_{(i)}^{\top}\mathbb{Y}_{(i)}-\boldsymbol{\tau}^{*\top}\boldsymbol{\Sigma}^{*{-1}}\boldsymbol{\tau}^{*}\Big{]},$
$\displaystyle\text{and}\ N_{k}=\sum_{i}I(z_{i}=k).$

The proofs of Propositions are given in supplementary materials. As shown in
Proposition 2, the conditional distribution of
$\\{\boldsymbol{\theta}_{k},\sigma_{k}^{2}\\}$ follows normal inverse gamma
distribution (NIG) with density function
$f(\boldsymbol{\tau}^{*},\boldsymbol{\Sigma}^{*},a^{*},b^{*})$. With the
expressions of $(\boldsymbol{\tau}^{*},\boldsymbol{\Sigma}^{*},a^{*},b^{*})$,
we see that it unifies the prior information
$\\{\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0}\\}$ and the data
information from the $k$th group, i.e.,
$\\{\mathbb{X}_{(i)},\mathbb{Y}_{(i)}:z_{i}=k\\}$. Based on Proposition 1 and
Proposition 2, we could efficiently cycle through the full conditional
distribution $\widetilde{f}(z_{i})$ for $i=1,2,\cdots,N$ and
$\widetilde{f}(\boldsymbol{\theta}_{k},\sigma_{k}^{2})$ for $k=1,2,\cdots,K$
to conduct Gibbs sampling procedure. One should note that in the process of
node memberships generation, $K$ is marginalized over, which could avoid
complicated reversible jump MCMC algorithms or allocation samplers. We
summarize the collapsed Gibbs sampling procedure for the GAGNAR model in
Algorithm 1.

Algorithm 1 Collapsed Gibbs sampler for GAGNAR

1:Observations $\mathbb{Y}_{(1)},\mathbb{Y}_{(2)},\dots,\mathbb{Y}_{(N)}$.

2:Initialize the number of groups $k$, the group assignment $z_{i}$ and the
parameters $(\boldsymbol{\theta}_{z_{i}},\sigma_{z_{i}}^{2})$ for each
$\mathbb{Y}_{(i)}$

3:for $iter$ from $1$ to $M$ do

4: for $i$ from $1$ to $N$ do

5: Remove $\mathbb{Y}_{(i)}$ from group $z_{i}$. If $\mathbb{Y}_{(i)}$ is from
a new cluster, set $k=k-1$.

6: Update $(z_{1},\dots,z_{N})$ conditional on
$\boldsymbol{\theta},\sigma^{2}$ for each node $i\in(1,\dots,N)$ by the
conditional distribution $\widetilde{f}(z_{i})$ in Proposition 1.

7: end for

8: for $c$ from $1$ to $K$ do

9: Update $(\boldsymbol{\theta}_{c},\sigma_{c}^{2})$ for each
$c\in(1,\dots,K)$ conditional on $(\mathbb{Z},\mathbb{Y},\mathbb{V})$ by the
density function in Proposition 2 as
$\displaystyle\widetilde{f}(\boldsymbol{\theta}_{c},\sigma_{c}^{2})\sim\text{NIG}(\boldsymbol{\tau}^{*},\boldsymbol{\Sigma}^{*},a^{*},b^{*})$

10: end for

11:end for

###  3.2 Post MCMC Estimation

Let $\boldsymbol{\theta}_{k}^{(m)}$ and $\sigma_{k}^{2(m)}$ with $1\leq k\leq
K^{(m)}$ be the $m$th estimate after the burn-in iterations of the MCMC
algorithm. Correspondingly denote $\mathbb{Z}^{(m)}=(z_{i}^{(m)}:1\leq i\leq
N)^{\top}$ as the estimated memberships of the network nodes. Particularly we
note that the estimated number of groups $K^{(m)}$ can be different for $1\leq
m\leq M$. As a result, direct posterior estimation of the parameters would be
difficult since the number of estimated parameters are not the same across
different iterations. To address this issue, we adopt Dahl’s method (Dahl,
2006) to select the best post burn-in iteration with the least squares
criterion. The estimate output by the best post burn-in iteration is then
taken as the final post MCMC estimator.

Specifically, we take advantage of the co-membership matrix to help us select
the best post burn-in iteration. Define the co-membership matrix as
$\mathbf{B}=(b_{ij})\in\mathbb{R}^{N\times N}$, where $b_{ij}=I(z_{i}=z_{j})$.
Therefore the $(i,j)$th element of $\mathbf{B}$ denotes whether the $i$th node
is in the same group with the $j$th node. We can see that $\mathbf{B}$ is well
defined for different $K^{(m)}$s and it does not suffer from the label
switching issue. We then choose the best post burn-in iteration as

$m_{b}=\mbox{argmin}_{1\leq m\leq
M}\big{\|}\mathbf{B}^{(m)}-\overline{\mathbf{B}}\big{\|}_{F}^{2},$

where $\mathbf{B}^{(m)}$ is the estimated co-membership matrix in the $m$th
post burn-in iteration and
$\overline{}\mathbf{B}=M^{-1}\sum_{m=1}^{M}\mathbf{B}^{(m)}$. As a
consequence, the best post burn-in iteration is selected by the closest
$\mathbf{B}^{(m)}$ to the mean grouping co-membership $\overline{}\mathbf{B}$.
The number of groups is then determined as $K^{(m_{b})}$. Accordingly, the
post MCMC estimations of the parameters and memberships are given by
$(\boldsymbol{\theta}_{k}^{(m_{b})},\sigma_{k}^{2(m_{b})})$ with $1\leq k\leq
K^{(m_{b})}$ and $\mathbb{Z}^{(m_{b})}$ respectively.

###  3.3 Selection of Smoothing Parameter $h$

The selection of smoothing parameter $h$ is redesigned as a model selection
problem. Specifically, we use the logarithm of the pseudo marginal likelihood
(LPML) (Ibrahim et al., 2014) based on conditional predictive ordinate (CPO)
(Gelfand et al., 1992; Geisser, 1993; Gelfand and Dey, 1994) to select $h$.

The LPML given a specified $h$ is defined as

$\displaystyle\mathrm{LPML}(h)=\sum_{i=1}^{N}\log\\{\mathrm{CPO}_{i}(h)\\},$
(3.2)

where CPO${}_{i}(h)$ is the CPO for the node $i$ defined as
CPO${}_{i}(h)=p(\mathbb{Y}_{(i)}|\mathbb{Y}_{(-i)})$ (Pettit, 1990), and
$\mathbb{Y}_{(-i)}=\\{\mathbb{Y}_{(j)}:j\neq i\\}$. As one can see, LPML is a
pseudo log-likelihood function and we prefer a smoothing parameter $h$ to
maximize LPML. Borrowing the idea of Chen et al. (2012), we obtain the Monte
Carlo estimate of CPO within the Bayesian framework as

$\displaystyle\widehat{\mathrm{CPO}}_{i}(h)=\Big{\\{}\frac{1}{M}\sum_{m=1}^{M}\frac{1}{L(\boldsymbol{\theta}^{(m)}_{z_{i}},\sigma^{2(m)}_{z_{i}};h)}\Big{\\}}^{-1},$

where $M$ is the total number of Monte Carlo iterations, and

$\displaystyle L(\boldsymbol{\theta}_{z_{i}}^{(m)},\sigma_{z_{i}}^{2(m)};h)$
$\displaystyle=\prod_{t=2}^{T}p(Y_{it}|\boldsymbol{\theta}_{z_{i}}^{(m)},\sigma_{z_{i}}^{2(m)};h)$

is the likelihood function for the node $i$ with the specified smoothing
parameter $h$. Correspondingly we have
$\widehat{\mathrm{LPML}}(h)=\sum_{i=1}^{N}\log\\{\widehat{\mathrm{CPO}}_{i}(h)\\}$
and we choose the optimal $h$ by
$h_{b}=\arg\max_{h}\widehat{\mathrm{LPML}}(h)$.

##  4 Simulation

###  4.1 Simulation Models and Data

To demonstrate the finite sample performance of our proposed method, we
conduct a number of numerical studies in this section. Specifically, we
utilize three types of graph and assign two parameter settings to each of
them.

For each simulation setting, we generate the random noise $\varepsilon_{it}$
from a standard normal distribution. In addition, node covariates
$V_{i}\in\mathbb{R}^{p}$ are independently sampled from a multivariate normal
distribution $\mathcal{N}(\mathbf{0},\mathbf{I}_{p})$. For each graph, two
different scenarios of parameters
$(\beta_{0k},\beta_{1k},\beta_{2k},\gamma_{k},\sigma_{k}^{2})$, which are
listed in Table 1. Given the initial value $\mathbb{Y}_{0}=\mathbf{0}$, the
time series $\mathbb{Y}_{t}$ is generated according to the grouped network
vector auto-regression model in (2.1). In each scenario, 100 replicated
datasets are generated, and in each replicate we set the time length $T=20$.
In addition, for the prior distribution gaCRP$(\alpha,h)$, we set $\alpha=1$
and select the smoothing parameters $h\in\\{0,0.2,0.4,\cdots,5.0\\}$ by LPML
(3.2). A total of 1500 MCMC iterations are run for each replicate, with the
first 500 iterations treated as burn-in. Besides, the hyper-parameters set for
the base distribution in (2.10) are
$(a_{0},b_{0},\boldsymbol{\tau}_{0},\boldsymbol{\Sigma}_{0})=(0.01,0.01,\mathbf{0},100\mathbf{I})$,
which are noninformative priors.

Table 1: Simulation parameters setting for three examples

| $\sigma^{2}$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\gamma$ | $\sigma^{2}$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\gamma$
---|---|---|---|---|---|---|---|---|---|---
Example 1
| Scenario 1 | Scenario 2
Group 1 | 2.0 | $5.0$ | $0.2$ | $0.1$ | $(0.5,0.7,1.0)^{\top}$ | 2.0 | $0.0$ | $0.1$ | $0.3$ | $(0.5,0.7,1.0)^{\top}$
Group 2 | 1.0 | $-5.0$ | $-0.4$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$ | 4.0 | $0.2$ | $-0.3$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$
Group 3 | 3.0 | $0.0$ | $0.2$ | $0.4$ | $(0.2,-1.0,2.0)^{\top}$ | 3.0 | $0.5$ | $0.2$ | $0.7$ | $(0.2,-0.2,1.4)^{\top}$
Example 2
| Scenario 1 | Scenario 2
Group 1 | 2.0 | $5.0$ | $0.2$ | $0.1$ | $(0.5,0.7,1.0)^{\top}$ | 2.0 | $0.0$ | $0.1$ | $0.3$ | $(0.5,0.7,1.0)^{\top}$
Group 2 | 1.0 | $-5.0$ | $-0.4$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$ | 1.0 | $0.2$ | $-0.3$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$
Group 3 | 3.0 | $0.0$ | $0.2$ | $0.4$ | $(0.2,-1.0,2.0)^{\top}$ | 3.0 | $0.5$ | $0.2$ | $0.7$ | $(0.2,-0.2,1.4)^{\top}$
Group 4 | 4.0 | $-0.1$ | $0.1$ | $0.2$ | $(1.0,-1.0,1.5)^{\top}$ | 4.0 | $-0.1$ | $0.1$ | $0.2$ | $(1.0,-1.0,1.5)^{\top}$
Group 5 | 2.0 | $3.0$ | $0.5$ | $0.2$ | $(0.8,0.5,-2.0)^{\top}$ | 2.0 | $0.8$ | $0.5$ | $0.2$ | $(0.8,0.5,-1.0)^{\top}$
Example 3
| Scenario 1 | Scenario 2
Group 1 | 2.0 | $5.0$ | $0.2$ | $0.1$ | $(0.5,0.7,1.0)^{\top}$ | 2.0 | $0.0$ | $0.1$ | $0.3$ | $(0.5,0.7,1.0)^{\top}$
Group 2 | 1.0 | $-5.0$ | $-0.4$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$ | 1.0 | $3.0$ | $-0.3$ | $0.2$ | $(0.1,0.9,0.4)^{\top}$
Group 3 | 3.0 | $0.0$ | $0.2$ | $0.4$ | $(0.2,-1.0,2.0)^{\top}$ | 3.0 | $-3.0$ | $0.2$ | $0.7$ | $(0.2,-0.2,1.4)^{\top}$
Group 4 | 4.0 | $3.0$ | $0.1$ | $0.2$ | $(1.0,-1.0,1.5)^{\top}$ | 1.5 | $4.5$ | $0.1$ | $0.2$ | $(1.0,-1.0,1.5)^{\top}$
Group 5 | 2.0 | $-3.0$ | $0.5$ | $0.2$ | $(0.8,0.5,-2.0)^{\top}$ | 2.5 | $-2.0$ | $0.5$ | $0.2$ | $(0.8,0.5,-1.0)^{\top}$
Group 6 | 3.0 | $2.0$ | $-0.6$ | $-0.2$ | $(-0.8,0.5,2.0)^{\top}$ | 1.0 | $2.0$ | $-0.6$ | $-0.2$ | $(-0.8,0.5,2.0)^{\top}$

Example 1. (Stochastic Block Model) We first consider the stochastic block
model (SBM) (Wang and Wong, 1987; Nowicki and Snijders, 2001; Zhao et al.,
2012). The SBM assumes that nodes in the same group (block) are more likely to
be connected, when compared with nodes from different groups. The model is
widely used to discover community structures for network data. We follow
Nowicki and Snijders (2001) to randomly assign each node a group label
$k=1,\cdots,K$ with equal probability $1/K$ and we set $K=3$. Next, let
$P(a_{ij}=1)=20N^{-1}$ if node $i$ and node $j$ are in the same group, and
$P(a_{ij}=1)=2N^{-1}$ otherwise. For this example we set the network size
$N=100$.

Example 2. (Chinese Cities Graph) In this example we construct the graph of
Chinese cities by using the geographical information. We collect 151 cities of
mainland China and treat each city as a network node in the graph. The edge
between two cities is defined as whether they share the common boarder (Cao et
al., 2017; Zhang et al., 2018). Specifically, $a_{ij}=1$ illustrates that two
cities are connected, otherwise $a_{ij}=0$. Lastly, we let $K=5$ be the total
number of groups and assign the group labels using $k$-means clustering on the
adjacency matrix.

Example 3. (Common Shareholder Network) We lastly construct the stock graph
from Chinese A stock market, which contains $N=180$ actively traded stocks in
the Shanghai and the Shenzhen Stock Exchange. The stocks used in this
simulation example is a subset of our empirical study for computational
convenience. Specifically, $a_{ij}=1$ indicates that two stocks share at least
one of the top two shareholders, otherwise $a_{ij}=0$. The graph structure is
shown in the bottom right panel of Figure 1. We set the number of groups $K=6$
and assign the group memberships by implementing the SBM to fit the adjacency
matrix.

The three network structures with the true group labels are visualized in
Figure 1.

Figure 1: Three graphs with true group labels, each color representing a
group. Left: graph generated from SBM with $K=3$. Middle: geographical graph
of $N=151$ Chinese cities with $K=5$. Right: common shareholder graph of
$N=180$ stocks with $K=6$.

###  4.2 Performance Measurements and Simulation Results

Let
$(\hat{\beta}_{0k}^{(r)},\hat{\beta}_{1k}^{(r)},\hat{\beta}_{2k}^{(r)},\hat{\gamma}_{k}^{(r)},\hat{\sigma}_{k}^{2(r)})$
be the estimated parameters of the $k$th group by the $r$th replicate ($1\leq
r\leq R$). For each node $i$, we obtain its group label as
$\hat{z}_{i}^{(r)}(i=1,\cdots,N)$ by Dalh’s method. Subsequently, the
estimated parameters for each node $i$ is given as
$\big{(}\hat{\beta}_{0\widehat{z}_{i}^{(r)}}^{(r)},\hat{\beta}_{1\widehat{z}_{i}^{(r)}}^{(r)},\hat{\beta}_{2\widehat{z}_{i}^{(r)}}^{(r)},\hat{\gamma}_{\widehat{z}_{i}^{(r)}}^{(r)},\hat{\sigma}_{\widehat{z}_{i}}^{2(r)}\big{)}$.
We consider the following measurements to evaluate the finite sample
performance. First, we employ the root mean square error (RMSE) to evaluate
the estimation accuracy. For $\beta_{sk}$ ($0\leq s\leq 2,1\leq k\leq K$), the
RMSE is defined as

$\displaystyle\operatorname{RMSE}_{\beta_{s}}=\Big{\\{}(RN)^{-1}\sum_{r=1}^{R}\sum_{i=1}^{N}\big{(}\hat{\beta}_{s\widehat{z}_{i}^{(r)}}^{(r)}-\beta_{sz_{i}}\big{)}^{2}\Big{\\}}^{1/2},$

and the RMSE for $\gamma$ is calculated as

$\displaystyle\operatorname{RMSE}_{\gamma}=\Big{\\{}(RN)^{-1}\sum_{r=1}^{R}\sum_{i=1}^{N}\big{\|}\hat{\gamma}_{\widehat{z}_{i}^{(r)}}^{(r)}-\gamma_{z_{i}}\big{\|}^{2}\Big{\\}}^{1/2}.$

To measure the similarity between the estimated group memberships and the true
group memberships, we use Adjusted Rand Index (ARI) (Rand, 1971; Hubert and
Arabie, 1985). Define $\\{z_{i}:1\leq i\leq N\\}$ as the set of the true group
memberships. Then ARI is defined as

$\displaystyle\text{ARI}^{(r)}=\frac{\text{RI}^{(r)}-E\big{(}\text{RI}^{(r)}\big{)}}{\max\big{(}\text{RI}^{(r)}\big{)}-E\big{(}\text{RI}^{(r)}\big{)}},~{}~{}~{}\text{RI}^{(r)}=\frac{a^{(r)}+b^{(r)}}{C_{N}^{2}},$

where
$a^{(r)}=\sum_{i,j}I(z_{i}=z_{j},\widehat{z}_{i}^{(r)}=\widehat{z}_{j}^{(r)})$,
$b^{(r)}=\sum_{i,j}I(z_{i}\neq
z_{j},\widehat{z}_{i}^{(r)}\neq\widehat{z}_{j}^{(r)})$, and
$C_{N}^{2}=N(N-1)/2$ is the total number of possible node pairs formed by the
$N$ nodes. Hence ARI measures the alignment level of two grouping results and
a higher ARI value implies the estimated group memberships are more consistent
with the true memberships. The ARI can be calculated using the R-package
mclust. We compare the parameter estimation of proposed model with that of EM
algorithm and two-step estimation for the GNAR model (Zhu and Pan, 2018). One
should note that, since the number of groups could not be inferred in the two
baseline methods, we apply the true value of $K$ as the input of EM and two-
step method.

Table 2 summarises the average RMSE of estimated parameters under the optimal
$h$ selected by LPML. For all three simulation examples, we see that the
overall RMSEs of estimated parameters by the GAGNAR model are much lower than
those of the EM and two-step methods. This is consistent with the weakness we
mentioned that the GNAR model does not utilize the network structure
information. In addition, we see that graphs in example 2 and 3 are more
complicated than that generated from SBM, due to the larger number of groups
and the group patterns which are not easily captured (see Figure 1). In both
examples, the estimation accuracy of our model is shown to be much higher than
that EM and two-step estimators. This further confirms the usefulness of the
proposed GAGNAR model in the intricate reality.

Figure 2 shows the estimated number of groups in the left panels, together
with ARIs in the right panels. When $h=0$, the gaCRP method reduces to the
traditional CRP method, which always tends to over-cluster and shows smaller
ARI than results by LPML method. We see that, when $h$ increases, the
estimated number of groups firstly decreases and then increases from
histograms in Figure 2 . The reason is that as $h\rightarrow\infty$, the
graphical weight $w_{ij}$ for disconnected nodes becomes particularly small,
which means only connected nodes can be classified into the same group,
therefore leading to the over-cluster problem, as we discussed in Remark 2.
The group concordance under optimal $h$ selected by LPML generally performs
better than those of EM method, with the proportion of selecting true number
of groups always greater than 80%.

Table 2: RMSE of parameter estimates. For each parameter, the left column
comes from scenario 1, and the right column is from scenario 2.

| $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\gamma$ | $\sigma^{2}$
---|---|---|---|---|---
Example 1
LPML | 0.626 | 0.471 | 0.179 | 0.377 | 0.057 | 0.108 | 0.264 | 0.532 | 0.330 | 0.624
EM | 1.706 | 0.233 | 0.164 | 0.167 | 0.132 | 0.123 | 0.514 | 0.279 | 1.450 | 1.466
2-step | 3.774 | 0.226 | 0.318 | 0.242 | 0.496 | 0.258 | 0.717 | 0.380 | 1.128 | 1.474
Example 2
LPML | 0.722 | 0.265 | 0.088 | 0.176 | 0.085 | 0.135 | 0.404 | 0.429 | 0.450 | 0.596
EM | 2.452 | 0.317 | 0.329 | 0.186 | 0.157 | 0.149 | 0.963 | 0.482 | 1.245 | 1.241
2-step | 3.191 | 0.370 | 0.329 | 0.282 | 0.508 | 0.371 | 0.956 | 0.726 | 1.188 | 1.309
Example 3
LPML | 0.920 | 0.644 | 0.158 | 0.180 | 0.093 | 0.089 | 0.460 | 0.318 | 0.443 | 0.246
EM | 1.942 | 1.416 | 0.300 | 0.288 | 0.110 | 0.120 | 0.748 | 0.573 | 1.183 | 0.576
2-step | 2.880 | 1.890 | 0.356 | 0.315 | 0.568 | 0.480 | 0.878 | 0.723 | 0.948 | 0.576

Figure 2: Grouping results of three examples, with results of scenario 1 in
the first line and scenario 2 in the second line. Left panel: histogram of
estimates of $K$ under different $h$. Right panel: boxplots of ARI under EM
method, $h=0$ and optimal $h$. The red texts show the average accuracy of
group memberships. The optimal $h$s selected by LPML are listed at the top of
each example.

##  5 Empirical Case Studies

In this section, we conduct two empirical case studies based on the graphs as
Examples 2–3 in Section 4. For the first real data example, we study the local
government economic competition by using city-level fiscal revenue (FR) data.
For the second real data example, we investigate the dynamic patterns of the
stocks traded in the Shanghai and Shenzhen Stock Exchange. We first conduct a
descriptive analysis of the real datasets. We calculate the average ratio of
fiscal revenue to GDP from 2002 to 2019 and the average weekly stock return
rate during 2019, as shown in Figure 3. The average FR/GDP ratio is 0.068 and
the weekly return of the stocks is 0.43%.

Figure 3: Left panel: histogram of $N=151$ cities’ yearly average fiscal
FR/GDP. Right panel: histogram of $N=384$ stocks’ weekly average returns.

###  5.1 China Fiscal Revenue

####  5.1.1 Background and Data Description

We first investigate how the city fiscal revenue relates to its spatial
adjacent cities and its historical information. The model with space–time lags
has been widely used on economic panel data research, including analysis about
equation systems with time dynamics and spatial spillover effects in regional
science (De Graaff et al., 2012; Gebremariam et al., 2011), fiscal policy with
government competition (Hauptmeier et al., 2012; Allers and Elhorst, 2011),
and many related fields (Korniotis, 2010; Brown and Laschever, 2012). In
addition, the latent group structure has been considered in recent macro-
economic researches. For instance, Bonhomme and Manresa (2015) proposed the
time-varying grouped patterns of heterogeneity in linear panel data models,
and model the relationship between the degree of democracy and GDP. They found
the group patterns across countries. A simple and fast approach was proposed
by Liu et al. (2020) to identify and estimate the unknown group structure in
panel data models, which is used to analyze the aggregate production function.
They discovered both individual-level and group-level heterogeneity for
women’s labor force participation.

We collect $N=151$ city public financial statements from Fiscal Statistics of
Cities and Counties in China during the period 2002-2019. The yearbook is
published by China Financial and Economic Publishing House, a state-owned
press under the supervision of the Ministry of Finance of the People’s
Republic of China. It collects detailed information on city-level fiscal
statistics, such as fiscal revenues and expenditures, fiscal accounting
balances, transfer payments, and the fiscally-supported population (Yu et al.,
2016). The response $Y_{it}$ is the fiscal revenue divided by local GDP.
Following Zhang and Zou (1998), Devereux et al. (2007) and Lv et al. (2020),
we consider four covariates, which are POP (population at the end of year),
GDP1st (the proportion of primary industry to GDP), SAV (year-end savings of
urban and rural residents), and FOR (actual foreign investment). Histogram in
the left panel of Figure 3 shows the average FR/GDP across all cities. A right
skewed distribution pattern can be observed.

####  5.1.2 Model Estimation

We next implement the GAGNAR model (2.10) to the China fiscal revenue data. We
use the data from 2005 to 2016 as the training set, and make prediction from
2006 to 2017. The smoothing parameter $h$ is chosen from
$\\{0,0.2,0.4,\cdots,2.0,3.0,4.0,5.0\\}$ and we set $\alpha=1$. For each $h$,
we run 1500 MCMC iterations and drop the first 500 as burn-in. The best $h$
selected by LPML is 2, and the estimated number of groups is $K=4$. The city
group pattern under $h=2$ is displayed in the top left panel of Figure 4.
Cities in the first group are mainly from Guangdong, Jiangxi, Zhejiang and
Anhui province, which are located in the south China. The second group
contains mainly the cities in Henan and Shandong province (east China), with
the average FR/GDP at the lowest level as the top right panel of Figure 4
shows. Most cities in Liaoning, Jilin, Heilongjiang provinces and others
constitute the third group, which are mostly in the Northeast China. Some
cities in Jiangsu are in the fourth group, with the highest average FR/GDP
values. According to the top right panel of Figure 4, cities in four groups
exhibit different economic features, as the average FR/GDPs among four groups
show different levels.

Table 3 reports the estimated parameter for each group under the optimal LPML
criterion. As one could see, the network effect and momentum effect are
different among four groups. Specifically, the cities in the second and the
third group have positive network effects, which means there might exist
imitation economic strategies for those cities with their spatial neighbors.
On contrary, cities in the first and the last group are negatively related to
their neighbors. This implies that they tend to take opposite macroeconomic
decisions to their adjacent cities. All of the four groups have positive
momentum effects, which shows that the responses are positively related to
their historical performances. Furthermore, the covariates tend to have
different effects on the response. For example, the population has negative
effect on the fiscal revenue for cities in the first and the second group
while its influence is positive for other cities. The foreign investment,
which represents the openness of a city, negatively affects the fiscal revenue
for cities in the fourth group, while it is positive for other cities. To
summarize the posterior distribution of each parameter, the Highest Posterior
Density (HPD) (Chen and Shao, 1999; Kruschke, 2014) interval that spans the
majority of the posterior distribution is specified. We calculate the 95% HPD
intervals of parameters of four typical cities (Beijing, Linyi, Fushun and
Taizhou) from four groups using the R package HDInterval. The results are
reported in the bottom left panel of Figure 4. Here the HPD interval is
calculated for each city specific parameters. According to the result, one
could see that though the coefficients of POP are similar among four groups,
the other parameters are quite different across the four cities.

Lastly, we compare the prediction accuracy of our model with CRP as well as
the GNAR model. Denote $T_{train}$ and $T_{test}$ as the time length of
training set and testing set respectively. To evaluate the prediction
accuracy, we calculate the mean square prediction error
$\text{MSPE}=(NT_{test})^{-1}\sum_{t=t_{test,0}}^{t_{test,0}+T_{test}}\sum_{i}(\widehat{Y}_{it}-Y_{it})^{2}$
where $\widehat{Y}_{it}$ is the predicted response for node $i$ at time point
$t$ and $t_{test,0}$ is the starting time point of the testing set. Define
$\text{MSPE}_{0}=(NT_{test})^{-1}\sum_{t=t_{test,0}}^{t_{test,0}+T_{test}}\sum_{i}(Y_{it}-\widehat{\mu}_{i,train})^{2}$
as the baseline MSPE, where
$\widehat{\mu}_{i,train}=T_{train}^{-1}\sum_{t=t_{test,0}-T_{train}}^{t_{test,0}-1}Y_{it}$
as the mean response of the $i$th node in the training set. The relative
prediction error is defined as

$\displaystyle\text{ReMSPE}=\text{MSPE}_{pred}/\text{MSPE}_{0},$ (5.1)

which is used to evaluate the model performance. As shown in the bottom right
panel of Figure 4, when the time length of training data is short, the
prediction accuracy of GNAR model is lower and not robust. The prediction
performances of our model and the CRP model are comparable, while our model
achieves slightly higher accuracy level.

Table 3: Parameter estimated by China Fiscal Revenue dataset with optimal
$h=2$.

| Group 1 | Group 2 | Group 3 | Group 4
---|---|---|---|---
$N_{k}$ | 83 | 38 | 25 | 5
$\widehat{\beta}_{0k}$ | 0.003 | 0.037 | 0.038 | 0.047
$\widehat{\beta}_{1k}$ | network effect | -0.011 | 0.101 | 0.212 | -0.179
$\widehat{\beta}_{2k}$ | momentum effect | 0.996 | 0.229 | 0.345 | 0.177
$\widehat{\gamma}_{k}$ | POP | -0.001 | -0.025 | 0.022 | 0.113
GDP1st | -0.001 | 0.175 | -0.128 | 0.186
SAV | 0.008 | -0.066 | 0.001 | -0.074
FOR | 0.008 | 0.071 | 0.078 | -0.383
$\widehat{\sigma}_{k}^{2}$ | 0.0001 | 0.0003 | 0.0006 | 0.0078

Figure 4: The top left panel shows the group pattern of cities estimated by
GAGNAR under $h=2$. The nodes in group 1–4 are marked in yellow, dark blue,
green and pink. The top right panel shows the average FR/GDP of cities
corresponding to four estimated groups. Bottom left panel shows the 95% HPD
intervals of four cities from different groups, and the ReMSPE of response in
2006–2017 is displayed in the bottom right panel.

###  5.2 Stock Return Analysis with Shareholder Network

####  5.2.1 Background and Data Description

The network information has been taken into consideration as an important
factor recently in analyzing financial time series. For instance, Zhu et al.
(2019) investigate the financial contagion phenomenon in the stock market with
a common shareholder network. To improve predictability of stocks return,
Leung et al. (2017) define a search-based cluster using co-search network of
stocks on Yahoo.com. In this study, we aim to analyze the stock returns using
a common shareholder network. Specifically, we collect $N=384$ stocks traded
in the Shanghai and the Shenzhen Stock Exchange. The response $Y_{it}$ is the
weekly stock return for 52 weeks in the year of 2019. Motivated by Fama and
French (2015), we consider six covariates related to corporations’
fundamentals. They are, SIZE (the logarithm of market value), BM (book to
market value), PR (yearly incremental profit ratio), AR (yearly incremental
asset ratio), LEV (leverage ratio), and CASH (cash flow). The histogram in the
right panel of Figure 3 shows the average return during the whole year,
calculated by $T^{-1}\sum_{t}Y_{it}$. The mean weekly return rate of all
stocks is given by 0.43%.

####  5.2.2 Model Estimation

To characterize the dynamics of the stock return, the GAGNAR model is
implemented to investigate the group patterns among the stocks. We use data in
the 52 traded weeks, with $T_{train}=40$ weeks used as the training set and
the remaining $T_{test}=12$ as the test set. The smoothing parameters are
selected from $h\in\\{0,0.2,0.4,0.6,\dots,3.0,4.0,5.0\\}$ For each $h$, 1500
MCMC iterations are run, and we drop the first 500 as burn-in. The best $h$
selected by LPML is 0.6, under which the number of groups is estimated as
$K=5$. The group pattern is shown in the top left panel of Figure 5. We then
look into the industry information of the five groups. Stocks in the first
group are mainly from the real estate industry. The second group has the
majority of stocks in information technology industry. Stocks from electrical
machinery and equipment manufacturing industry constitute the major part of
the third group, while stocks from the production and supply of electricity,
steam and hot water form the fourth group. Lastly, the chemical related
industries as well as civil engineering construction are included in the fifth
group. In addition, the top right panel of Figure 5 shows that the stock
return levels are diversified in the five groups. Particularly, the third
group achieves the highest average weekly return rate, while stock returns of
the first and fourth groups are relatively low.

The parameter estimation is shown in Table 4. As presented in the result, the
network effect and momentum effect are different among five groups.
Specifically, the stocks in the second group are positively related to their
connected stocks, while the stocks in the rest four group behave oppositely.
The opposite network effects are also observed in empirical studies in
financial stock market (Borgatti et al., 2009; Peng and Dey, 2013; Chen et
al., 2019). The momentum effects among the first four groups are negative,
which shows that they are conversely related to their historical return rate.
The last group, on contrary, is positively related to its historical
behaviors. Furthermore, we observe that the covariates have different effects
on the response. For example, the market value (SIZE) of stocks in the first,
second and the fourth group positively contribute to the weekly traded return.
The cash flow (CASH) of the stocks in the second group is positively related
to the return rate, while it is not the case for other four groups. Next, the
95% HPD intervals are depicted in the bottom left panel of Figure 5, from
which one could see that there are obvious dissimilarity in the network effect
and momentum effect among stocks in different groups. Though there is slight
difference in the covariates coefficients among five groups, our proposed
model could capture the imperceptible distinction for different nodes.

Lastly, to evaluate the out-sample prediction performance of the GAGNAR model,
we also use ReMSPE (5.1) to compare the prediction accuracy of GAGNAR with CRP
and the GNAR model. The results are shown in bottom right panel of Figure 5.
As one could observe, the GAGNAR model achieves slightly higher and more
robust prediction accuracy than other competing methods. This illustrates the
usefulness of the proposed methodology.

Table 4: Parameter estimated ($\times$10) by Stock Return dataset with optimal
$h=0.6$.

| Group 1 | Group 2 | Group 3 | Group 4 | Group 5
---|---|---|---|---|---
$N_{k}$ | 92 | 94 | 124 | 71 | 3
$\widehat{\beta}_{0k}$ | 0.008 | 0.154 | 0.029 | -0.007 | 6.164
$\widehat{\beta}_{1k}$ | network effect | -0.313 | 0.343 | -0.678 | -0.172 | -7.776
$\widehat{\beta}_{2k}$ | momentum effect | -1.066 | -0.541 | -0.692 | -1.032 | 1.795
$\widehat{\gamma}_{k}$ | SIZE | 0.006 | 0.009 | -0.006 | 0.001 | -4.968
BM | -0.003 | -0.01 | 0.03 | -0.012 | 9.136
PR | -0.006 | 0.007 | -0.014 | 0.011 | -2.988
AR | -0.016 | 0.077 | 0.008 | 0.001 | -4.028
| LEV | 0.011 | 0.033 | -0.018 | 0.003 | 11.44
| CASH | -0.001 | 0.023 | -0.002 | -0.009 | -13.783
$\widehat{\sigma}_{k}^{2}$ | 0.025 | 0.074 | 0.037 | 0.012 | 0.228

Figure 5: The top left panel shows the group pattern of stocks estimated by
GAGNAR under $h=0.6$. The nodes in group 1–5 are marked in yellow, dark blue,
green, pink, and skyblue. The top right panel shows the average return rate of
stocks corresponding to five estimated groups. Bottom left panel shows the 95%
HPD intervals of five stocks, which are the first stock of the five groups
sorted by the stock code respectively. The ReMSPE in week 41–52 is displayed
in the bottom right panel.

##  6 Discussion

In this paper, we propose a novel Bayesian nonparametric grouping approach for
learning the heterogeneity of dynamic response for network data. The proposed
method is able to conduct model estimation, estimate the number of groups, and
infer group configurations simultaneously. Building upon the gaCRP framework,
we develop a collapsed Gibbs sampler for efficient Bayesian inference.
Specifically we adopt Dahl’s method for post MCMC inference and introduced
LPML for smoothing parameter selection. Numerical results have confirmed that
the proposed method is able to simultaneously infer the number of groups and
the group-wise parameters with high accuracy. Comparing to the traditional
techniques such as EM algorithm and two-step approach in Zhu and Pan (2018),
the proposed method is able to improve the grouping performance, especially
when grouping configurations contain certain graphical information. Lastly, we
illustrate the usefulness of the proposed method by studying two real data
examples, including China fiscal revenue analysis and stock return prediction.
The results indicate that incorporating graphical information on grouping
process will improve the interpretability and prediction power of the
methodology.

A few topics beyond the scope of this paper are worth further investigation. A
natural extension for this paper is for large scale networks. Designing an
efficient algorithm for massive datasets with lower computational complexity
devotes an interesting future work. Two promising solutions for tackling this
challenge are discovering low rank structure and imposing sparsity. Second,
the proposed method is restricted to continuous responses and thus can be
extended to binary or counted data. Finally, incorporating the GNAR model with
the stochastic block model (SBM) will be another interesting topic for future
study.

## References

  * Allers and Elhorst (2011) Allers, M. A. and J. P. Elhorst (2011).  A simultaneous equations model of fiscal policy interactions.  Journal of regional science 51(2), 271–291.
  * Ando and Bai (2016) Ando, T. and J. Bai (2016).  Panel data models with grouped factor structure under unknown group membership.  Journal of Applied Econometrics 31(1), 163–191.
  * Bonhomme and Manresa (2015) Bonhomme, S. and E. Manresa (2015).  Grouped patterns of heterogeneity in panel data.  Econometrica 83(3), 1147–1184.
  * Borgatti et al. (2009) Borgatti, S. P., A. Mehra, D. J. Brass, and G. Labianca (2009).  Network analysis in the social sciences.  science 323(5916), 892–895.
  * Borsboom and Cramer (2013) Borsboom, D. and A. O. Cramer (2013).  Network analysis: an integrative approach to the structure of psychopathology.  Annual review of clinical psychology 9, 91–121.
  * Brown and Laschever (2012) Brown, K. M. and R. A. Laschever (2012).  When they’re sixty-four: Peer effects and the timing of retirement.  American Economic Journal: Applied Economics 4(3), 90–115.
  * Cao et al. (2017) Cao, Q., Y. Liang, and X. Niu (2017).  China’s air quality and respiratory disease mortality based on the spatial panel model.  International journal of environmental research and public health 14(9), 1081.
  * Chen and Shao (1999) Chen, M.-H. and Q.-M. Shao (1999).  Monte carlo estimation of bayesian credible and hpd intervals.  Journal of Computational and Graphical Statistics 8(1), 69–92.
  * Chen et al. (2012) Chen, M.-H., Q.-M. Shao, and J. G. Ibrahim (2012).  Monte Carlo methods in Bayesian computation.  Springer Science & Business Media.
  * Chen et al. (2019) Chen, X., W. Shangguan, Y. Liu, and S. Wang (2019).  Can network structure predict cross-sectional stock returns? evidence from co-attention networks in china.  Finance Research Letters, 101422.
  * Chuang et al. (2007) Chuang, H.-Y., E. Lee, Y.-T. Liu, D. Lee, and T. Ideker (2007).  Network-based classification of breast cancer metastasis.  Molecular systems biology 3(1), 140.
  * Cramer et al. (2010) Cramer, A. O., L. J. Waldorp, H. L. Van Der Maas, and D. Borsboom (2010).  Comorbidity: A network perspective.  Behavioral and brain sciences 33(2-3), 137–150.
  * Dahl (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a dirichlet process mixture model.  Bayesian inference for gene expression and proteomics 4, 201–218.
  * De Graaff et al. (2012) De Graaff, T., F. G. Van Oort, and R. J. Florax (2012).  Regional population–employment dynamics across different sectors of the economy.  Journal of Regional Science 52(1), 60–84.
  * Devereux et al. (2007) Devereux, M. P., B. Lockwood, and M. Redoano (2007).  Horizontal and vertical indirect tax competition: Theory and some evidence from the usa.  Journal of Public Economics 91(3-4), 451–479.
  * Fama and French (2015) Fama, E. F. and K. R. French (2015).  A five-factor asset pricing model.  Journal of Financial Economics 116(1), 1–22.
  * Gebremariam et al. (2011) Gebremariam, G. H., T. G. Gebremedhin, and P. V. Schaeffer (2011).  Employment, income, and migration in appalachia: A spatial simultaneous equations approach.  Journal of Regional Science 51(1), 102–120.
  * Geisser (1993) Geisser, S. (1993).  Predictive inference, Volume 55.  CRC press.
  * Gelfand and Dey (1994) Gelfand, A. E. and D. K. Dey (1994).  Bayesian model choice: asymptotics and exact calculations.  Journal of the Royal Statistical Society: Series B (Methodological) 56(3), 501–514.
  * Gelfand et al. (1992) Gelfand, A. E., D. K. Dey, and H. Chang (1992).  Model determination using predictive distributions with implementation via sampling-based methods.  Technical report, STANFORD UNIV CA DEPT OF STATISTICS.
  * Geng et al. (2019) Geng, J., A. Bhattacharya, and D. Pati (2019).  Probabilistic community detection with unknown number of communities.  Journal of the American Statistical Association 114(526), 893–905.
  * Geng and Hu (2021) Geng, L. and G. Hu (2021).  Bayesian spatial homogeneity pursuit for survival data with an application to the seer respiratory cancer data.  Biometrics.
  * Green (1995) Green, P. J. (1995).  Reversible jump markov chain monte carlo computation and bayesian model determination.  Biometrika 82(4), 711–732.
  * Guda et al. (2009) Guda, P., S. V. Chittur, and C. Guda (2009).  Comparative analysis of protein-protein interactions in cancer-associated genes.  Genomics, proteomics & bioinformatics 7(1-2), 25–36.
  * Guðmundsson and Brownlees (2021) Guðmundsson, G. S. and C. Brownlees (2021).  Detecting groups in large vector autoregressions.  Journal of Econometrics.
  * Hanneman and Riddle (2005) Hanneman, R. A. and M. Riddle (2005).  Introduction to social network methods.
  * Härdle et al. (2016) Härdle, W. K., W. Wang, and L. Yu (2016).  Tenet: Tail-event driven network risk.  Journal of Econometrics 192(2), 499–513.
  * Hauptmeier et al. (2012) Hauptmeier, S., F. Mittermaier, and J. Rincke (2012).  Fiscal competition over taxes and public inputs.  Regional science and urban economics 42(3), 407–419.
  * Horvath (2011) Horvath, S. (2011).  Weighted network analysis: applications in genomics and systems biology.  Springer Science & Business Media.
  * Hubert and Arabie (1985) Hubert, L. and P. Arabie (1985).  Comparing partitions.  Journal of classification 2(1), 193–218.
  * Ibrahim et al. (2014) Ibrahim, J. G., M.-H. Chen, and D. Sinha (2014).  Bayesian survival analysis.  Wiley StatsRef: Statistics Reference Online.
  * Ishwaran and Zarepour (2002) Ishwaran, H. and M. Zarepour (2002).  Exact and approximate sum representations for the dirichlet process.  Canadian Journal of Statistics 30(2), 269–283.
  * Ke et al. (2015) Ke, Z. T., J. Fan, and Y. Wu (2015).  Homogeneity pursuit.  Journal of the American Statistical Association 110(509), 175–194.
  * Korniotis (2010) Korniotis, G. M. (2010).  Estimating panel models with internal and external habit formation.  Journal of Business & Economic Statistics 28(1), 145–158.
  * Kruschke (2014) Kruschke, J. (2014).  Doing bayesian data analysis: A tutorial with r, jags, and stan.
  * Lee and Yu (2009) Lee, L.-f. and J. Yu (2009).  Spatial nonstationarity and spurious regression: the case with a row-normalized spatial weights matrix.  Spatial Economic Analysis 4(3), 301–327.
  * Leung et al. (2017) Leung, A. C. M., A. Agarwal, P. Konana, and A. Kumar (2017).  Network analysis of search dynamics: The case of stock habitats.  Management Science 63(8), 2667–2687.
  * Liu et al. (2020) Liu, R., Z. Shang, Y. Zhang, and Q. Zhou (2020).  Identification and estimation in panel models with overspecified number of groups.  Journal of Econometrics 215(2), 574–590.
  * Liu et al. (2017) Liu, X., E. Patacchini, and E. Rainone (2017).  Peer effects in bedtime decisions among adolescents: a social network model with sampled data.  The econometrics journal 20(3), S103–S125.
  * Lv et al. (2020) Lv, B., Y. Liu, and Y. Li (2020).  Fiscal incentives, competition, and investment in china.  China Economic Review 59, 101371.
  * Marin et al. (2005) Marin, J.-M., K. Mengersen, and C. P. Robert (2005).  Bayesian modelling and inference on mixtures of distributions.  Handbook of statistics 25, 459–507.
  * Miller and Harrison (2018) Miller, J. W. and M. T. Harrison (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association 113(521), 340–356.
  * Nowicki and Snijders (2001) Nowicki, K. and T. A. B. Snijders (2001).  Estimation and prediction for stochastic blockstructures.  Journal of the American statistical association 96(455), 1077–1087.
  * O’malley and Marsden (2008) O’malley, A. J. and P. V. Marsden (2008).  The analysis of social networks.  Health services and outcomes research methodology 8(4), 222–269.
  * Peng and Dey (2013) Peng, G. and D. Dey (2013).  Research note—a dynamic view of the impact of network structure on technology adoption: the case of oss development.  Information Systems Research 24(4), 1087–1099.
  * Pettit (1990) Pettit, L. (1990).  The conditional predictive ordinate for the normal distribution.  Journal of the Royal Statistical Society: Series B (Methodological) 52(1), 175–184.
  * Pitman (1995) Pitman, J. (1995).  Exchangeable and partially exchangeable random partitions.  Probability theory and related fields 102(2), 145–158.
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical association 66(336), 846–850.
  * Sewell and Chen (2017) Sewell, D. K. and Y. Chen (2017).  Latent space approaches to community detection in dynamic networks.  Bayesian analysis 12(2), 351–377.
  * Shi and Lee (2017) Shi, W. and L.-f. Lee (2017).  Spatial dynamic panel data models with interactive fixed effects.  Journal of Econometrics 197(2), 323–347.
  * Simmel (1950) Simmel, G. (1950).  The sociology of georg simmel, Volume 92892.  Simon and Schuster.
  * Sojourner (2013) Sojourner, A. (2013).  Identification of peer effects with missing peer data: Evidence from project star.  The Economic Journal 123(569), 574–605.
  * Su et al. (2016) Su, L., Z. Shi, and P. C. Phillips (2016).  Identifying latent structures in panel data.  Econometrica 84(6), 2215–2264.
  * van der Pas and van der Vaart (2018) van der Pas, S. and A. van der Vaart (2018).  Bayesian community detection.  Bayesian Analysis 13(3), 767–796.
  * Von Mering et al. (2002) Von Mering, C., R. Krause, B. Snel, M. Cornell, S. G. Oliver, S. Fields, and P. Bork (2002).  Comparative assessment of large-scale data sets of protein–protein interactions.  Nature 417(6887), 399–403.
  * Wang and Wong (1987) Wang, Y. J. and G. Y. Wong (1987).  Stochastic blockmodels for directed graphs.  Journal of the American Statistical Association 82(397), 8–19.
  * Wasserman et al. (1994) Wasserman, S., K. Faust, et al. (1994).  Social network analysis: Methods and applications.
  * Wu et al. (2010) Wu, G., X. Feng, and L. Stein (2010).  A human functional protein interaction network and its application to cancer data analysis.  Genome biology 11(5), 1–23.
  * Yu et al. (2016) Yu, J., L.-A. Zhou, and G. Zhu (2016).  Strategic interaction in political competition: Evidence from spatial effects across chinese cities.  Regional Science and Urban Economics 57, 23–37.
  * Zhang et al. (2018) Zhang, J., Y. Chang, L. Zhang, and D. Li (2018).  Do technological innovations promote urban green development?—a spatial econometric analysis of 105 cities in china.  Journal of cleaner production 182, 395–403.
  * Zhang and Zou (1998) Zhang, T. and H.-f. Zou (1998).  Fiscal decentralization, public spending, and economic growth in china.  Journal of public economics 67(2), 221–240.
  * Zhao et al. (2012) Zhao, Y., E. Levina, J. Zhu, et al. (2012).  Consistency of community detection in networks under degree-corrected stochastic block models.  The Annals of Statistics 40(4), 2266–2292.
  * Zhu and Pan (2018) Zhu, X. and R. Pan (2018).  Grouped network vector autoregression.  Statistica Sinca.
  * Zhu et al. (2017) Zhu, X., R. Pan, G. Li, Y. Liu, H. Wang, et al. (2017).  Network vector autoregression.  The Annals of Statistics 45(3), 1096–1123.
  * Zhu et al. (2019) Zhu, X., W. Wang, H. Wang, and W. K. Härdle (2019).  Network quantile autoregression.  Journal of econometrics 212(1), 345–358.
  * Zou et al. (2017) Zou, T., W. Lan, H. Wang, and C.-L. Tsai (2017).  Covariance regression analysis.  Journal of the American Statistical Association 112(517), 266–281.