# Random Networks with Heterogeneous Reciprocity

Tiandong<EMAIL_ADDRESS>[ Sidney<EMAIL_ADDRESS>[ Department of Statistics, Texas A&M University, College Station, TX 77843,
USpresep=, ]e1  School of Operations Research and Information Engineering,
Cornell University, Ithaca, NY 14853, USpresep=, ]e2

###### Abstract

Users of social networks display diversified behavior and online habits. For
instance, a user’s tendency to reply to a post can depend on the user and the
person posting. For convenience, we group users into aggregated behavioral
patterns, focusing here on the tendency to reply to or _reciprocate_ messages.
The reciprocity feature in social networks reflects the information exchange
among users. We study the properties of a preferential attachment model with
heterogeneous reciprocity levels, give the growth rate of model edge counts,
and prove convergence of empirical degree frequencies to a limiting
distribution. This limiting distribution is not only multivariate regularly
varying, but also has the property of hidden regular variation.

60G70,

60J85,

Multivariate regular variation,

Hidden regular variation,

Preferential attachment,

Reciprocity,

###### keywords:

[class=MSC]

###### keywords:

and

##  1 Introduction

Social networks have grown rapidly and users are exhibiting different
interaction and behavioral patterns on platforms like Facebook and Twitter.
Reciprocity is one such pattern and helps characterize information exchange
between social network users [12, 18]. For instance, consider the network of
Facebook wall posts: a directed edge from user A to user B is formed when user
A leaves a message on the Facebook wall of user B, and a reciprocal edge from
user B to user A is created if and when user B replies to the message.
Depending on both users’ behavioral features and networking habits (e.g. the
closeness of friendship, whether user A is broadcasting, the obsessiveness of
user B’s replying habits), the probability of generating a reciprocal edge may
vary across different pairs of users, and network modeling should include such
heterogeneity. A study in [11] shows that online social networks tend to have
a high proportion of reciprocal edges, compared to other types of networks
such as biological networks, communication networks, software call graphs and
peer-to-peer networks. Due to the large and diversified user groups on social
networks, it is essential to incorporate heterogeneous reciprocity levels into
the modeling.

For the modeling of dynamic networks, the preferential attachment (PA) model
[5, 14] is an appealing starting point. This model captures the scale-free
property of complex networks, where both in- and out-degree distributions have
Pareto-like tails [22, 21, 23, 24]. Recently, [25] show that for a wide and
realistic range of model parameters, the standard PA model in [5] generates
networks with very small proportion of reciprocal edges, deviating from the
empirical observations in [11]. This led to two consecutive studies [26, 6] on
the theoretical properties and the estimation of a PA model with homogeneous
reciprocity, where reciprocal edges are generated by simply flipping a two-
sided coin. To capture the heterogeneous reciprocal patterns, we extend here
the model in [26] by dividing users into $K$ different behavioral groups, and
the probability of generating a reciprocal edge from a user in group $m$ to a
user in group $r$ is $\rho_{m,r}\in(0,1)$, $m,r\in\\{1,\ldots,K\\}$.

We study three properties of the proposed PA model with heterogeneous
reciprocity. Under modest assumptions, we first analyze the growth of the edge
counts by identifying the almost sure limit for the scaled number of edges
emanating from and pointing to nodes of group type $m$. Secondly, by embedding
in- and out-degree sequences into a family of multitype branching processes
whose particles are given group labels, we prove that the empirical
frequencies of nodes with in-degree $i$ and out-degree $j$ converge to a
limiting distribution $p_{i,j}$. Third, we show that the asymptotic limiting
distribution, $p_{i,j}$, is multivariate regularly varying with limit measure
concentrating on a ray and after removing large in- and out-degree pairs close
to the concentrating ray, we also detect hidden regular variation [16].

The rest of the paper is organized as follows. We start with a detailed
description of the PA model with heterogeneous reciprocity levels in Section
1.1. Section 2 studies the growth of edge counts, and in order to ensure the
convergence of scaled edge counts, sufficient assumptions are imposed. Also,
we extend the embedding technique in [26] to a family of multitype branching
processes with different group labels, and derive the limit of empirical
degree frequencies in Section 3. We then characterize the asymptotic
dependence structure of large in- and out-degrees in Section 4, and give
concluding remarks in Section 5. Technical proofs of results in Section 2 are
collected in Section 6.

###  1.1 PA Model with Heterogeneous Reciprocity

The proposed model extends the directed preferential attachment (PA) model
studied in [14, 5, 23, 24] by amending a mechanism that generates a reciprocal
edge with probabilities depending on characteristics of the node pair being
connected by the edge. Because these characteristics can be specific to the
node pair, the model incorporates heterogeneous responding patterns of users
in social networks.

We now specify a growing sequence of graphs. Let $G(n)$ be the graph after $n$
steps with $V(n)$ being the set of nodes and $E(n)$ being the set of edges in
$G(n)$. Attach to each node $v$ a communication type $W_{v}$, where
$\\{W_{v},v\geq 1\\}$ are iid random variables, independent from the growth
mechanism of the graph, with

$\mathbb{P}(W_{v}=r)=\pi_{r},\qquad\text{for}\quad\sum_{r=1}^{K}\pi_{r}=1.$

We imagine that when a node is born, it flips a multi-sided coin to determine
its communication type. Think of $\pi_{r}$ as the percent of the users in a
social network that have communication habits labeled as type $r$. Let
$W(n):=\\{W_{v}:v\in V(n)\\}$ denote the set of group types for all nodes in
$G(n)$. Throughout we assume that the communication group of node $v$ is
always observed upon its creation, and remains unchanged afterwards. In the
rest of this paper, we assume $G(n)=(V(n),E(n),W(n))$, for $n\geq 0$.

Denote the cardinality of a discrete set $S$ by $|S|$ and initialize the model
with graph $G(0)$, which consists of one node (labeled as node 1) and a self-
loop, with $V(0)=\\{1\\}$, $|V(0)|={1}$, $W(0)=\\{W_{1}\\}$, and
$E(0)=\left\\{(1,1)\right\\}$. For each new edge $(u,v)$ with
$W_{u}=r,W_{v}=m$, the reciprocity mechanism adds its reciprocal counterpart
$(v,u)$ instantaneously with probability $\rho_{m,r}\in[0,1]$, for
$m,r\in\\{1,2,\ldots,K\\}$. Here $\rho_{m,r}$ measures the probability of
adding a reciprocal edge from a node in group $m$ to a node in group $r$. Note
that the matrix $\bm{\rho}:=(\rho_{m,r})_{m,r}$ is not necessarily a
stochastic matrix, but can be an arbitrary matrix in $M_{K\times K}([0,1])$,
the set of all $K\times K$ matrices with entries belonging to $[0,1]$. Later
in Section 2, we will give particular regularity conditions on $\bm{\rho}$ to
facilitate theoretical analysis.

Let $\bigl{(}D^{\text{in}}_{v}(n),D^{\text{out}}_{v}(n)\bigr{)}$ be the in-
and out-degrees of node $v\in V(n)$ in $G(n)$, and we use the convention that
$D^{\text{in}}_{v}(n)=D^{\text{out}}_{v}(n)=0$ if $v\notin V(n)$. Let
$\delta>0$ be an offset parameter. The evolution of the network $G(n+1)$ from
$G(n)$ is described as follows.

  1. 1.

With probability $\alpha\in(0,1)$, add a new node $|V(n)|+1$ with a directed
edge $(|V(n)|+1,v)$, where $v\in V(n)$ is chosen with probability

$\frac{D^{\text{in}}_{v}(n)+\delta}{\sum_{v\in
V(n)}(D^{\text{in}}_{w}(n)+\delta)}=\frac{D^{\text{in}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|},$
(1)

and update the node set $V(n+1)=V(n)\cup\\{|V(n)|+1\\}$ and
$W(n+1)=W(n)\cup\\{W_{|V(n)|+1}\\}$. The new node $|V(n)|+1$ belongs to group
$r$ with probability $\pi_{r}$. If node $v$ belongs to group $m$, then a
reciprocal edge $(v,|V(n)|+1)$ is added with probability $\rho_{m,r}$. Upon
reciprocation, update the edge set as
$E(n+1)=E(n)\cup\\{(|V(n)|+1,v),(v,|V(n)|+1)\\}$. If the reciprocal edge is
not created, set $E(n+1)=E(n)\cup\\{(|V(n)|+1,v)\\}$.

  2. 2.

With probability $\gamma\equiv 1-\alpha\in(0,1)$, add a new node $|V(n)|+1$
with a directed edge $(v,|V(n)|+1)$, where $v\in V(n)$ is chosen with
probability

$\frac{D^{\text{out}}_{v}(n)+\delta}{\sum_{v\in
V(n)}(D^{\text{out}}_{v}(n)+\delta)}=\frac{D^{\text{out}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|},$
(2)

and update the node set $V(n+1)=V(n)\cup\\{|V(n)|+1\\}$,
$W(n+1)=W(n)\cup\\{W_{|V(n)|+1}\\}$. The new node $|V(n)|+1$ belongs to group
$r$ with probability $\pi_{r}$. If node $v$ belongs to group $m$, then a
reciprocal edge $(|V(n)|+1,v)$ is added with probability $\rho_{r,m}$. Upon
reciprocation, update the edge set as
$E(n+1)=E(n)\cup\\{(v,|V(n)|+1,v),(|V(n)|+1,v)\\}$. If the reciprocal edge is
not created, set $E(n+1)=E(n)\cup\\{(v,|V(n)|+1)\\}$.

Note that $|V(n)|=n+1$, $n\geq 0$, since a new node is added at each step, and
the offset parameter $\delta$ is assumed to be the same for both in- and out-
degrees.

##  2 Growth of Edge Counts

Heterogeneous reciprocity levels make analysis of the convergence of
$|E(n)|/n$ more complicated compared to the homogeneous case. In Theorem 2.2,
we show the concentration of $|E(n)|$ around $\mathbb{E}[|E(n)|]$, and then
prove the convergence of $\mathbb{E}[|E(n)|]/n$, thereby showing linear growth
of the number of edges. We now prepare for this result.

Write

$\displaystyle\rho_{m\bullet}$
$\displaystyle:=\sum_{r=1}^{K}\rho_{m,r}\pi_{r}=\mathbb{P}[\text{Type $m$ node
sends reciprocal edge}],$ $\displaystyle\rho_{\bullet m}$
$\displaystyle:=\sum_{r=1}^{K}\rho_{r,m}\pi_{r}=\mathbb{P}[\text{Type $m$ node
receives reciprocal edge}].$

Let $\mathcal{G}_{n}$ be denote the $\sigma$-field generated by observing the
network evolution up to $n$ steps, i.e.

$\mathcal{G}_{n}=\sigma\left\\{(V(k),E(k),W(k)):k=0,\ldots,n\right\\}.$

For $v\in V(n)$, we have

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}\left(D^{\text{in}}_{v}(n+1)-D^{\text{in}}_{v}(n)\right)$
$\displaystyle=\alpha\frac{D^{\text{in}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|}+\gamma\sum_{m=1}^{K}\rho_{\bullet
m}\frac{D^{\text{out}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|}\bm{1}_{\\{W_{v}=m\\}},$
(3)
$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}\left(D^{\text{out}}_{v}(n+1)-D^{\text{out}}_{v}(n)\right)$
$\displaystyle=\gamma\frac{D^{\text{out}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|}+\alpha\sum_{m=1}^{K}\rho_{m\bullet}\frac{D^{\text{in}}_{v}(n)+\delta}{|E(n)|+\delta|V(n)|}\bm{1}_{\\{W_{v}=m\\}}.$
(4)

Let $V_{m}(n):=\\{v\in V(n):W_{v}=m\\}$, and note that

$\displaystyle\sum_{v\in V(n)}$
$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}\left(D^{\text{in}}_{v}(n+1)-D^{\text{in}}_{v}(n)\right)\bm{1}_{\\{W_{v}=m\\}}$
$\displaystyle=$
$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}E(n+1)=E(n)\cup\\{(|V(n)|+1,v)\\},v\in
V_{m}(n)\bigr{)}$
$\displaystyle+\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}E(n+1)=E(n)\cup\\{(v,|V(n)|+1),(|V(n)|+1,v)\\},$
$\displaystyle\qquad\qquad v\in V_{m}(n),\text{$(|V(n)|+1,v)$ is due to
reciprocation.}\bigr{)}$ $\displaystyle=$
$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(|V(n)|+1,v)\in
E(n+1)\setminus E(n),v\in V_{m}(n)\bigr{)}.$ (5)

Write $E^{\text{in}}_{m}(n):=\\{(u,v)\in E(n):W_{v}=m\\}$, and
$E^{\text{out}}_{m}(n):=\\{(v,u)\in E(n):W_{v}=m\\}$. We see that

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}$
$\displaystyle\left(|E^{\text{in}}_{m}(n+1)|-|E^{\text{in}}_{m}(n)|\right)$
$\displaystyle=\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(|V(n)|+1,v)\in
E(n+1)\setminus E(n),v\in V_{m}(n)\bigr{)}$
$\displaystyle+\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n)\bigr{)},$ and applying (3) and
(5) gives $\displaystyle=\sum_{v\in
V(n)}\mathbb{E}^{\mathcal{G}_{n}}\left(D^{\text{in}}_{v}(n+1)-D^{\text{in}}_{v}(n)\right)\bm{1}_{\\{W_{v}=m\\}}$
$\displaystyle+\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n)\bigr{)}$
$\displaystyle=\alpha\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}+\gamma\rho_{\bullet
m}\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}$
$\displaystyle+\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n)\bigr{)}.$ (6)

For the third term in (6), we have that

$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n)\bigr{)}$
$\displaystyle=\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n),$
$\displaystyle\qquad\qquad\text{$(v,|V(n)|+1)$ is not due to
reciprocation}\bigr{)}$
$\displaystyle\quad+\mathbb{P}^{\mathcal{G}_{n}}\bigl{(}(v,|V(n)|+1)\in
E(n+1)\setminus E(n),W_{|V(n)|+1}=m,v\in V(n),$
$\displaystyle\quad\qquad\qquad\text{$(v,|V(n)|+1)$ is due to
reciprocation},W_{v}=r,\text{ for some $r$}\bigr{)}$
$\displaystyle=\gamma\pi_{m}+\alpha\pi_{m}\sum_{r=1}^{K}\rho_{r,m}\frac{|E^{\text{in}}_{r}(n)|+\delta|V_{r}(n)|}{|E(n)|+\delta|V(n)|}.$
(7)

Therefore, combining (6) and (7) shows that

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}$
$\displaystyle\left(|E^{\text{in}}_{m}(n+1)|-|E^{\text{in}}_{m}(n)|\right)$
$\displaystyle=$
$\displaystyle\alpha\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}+\gamma\rho_{\bullet
m}\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}$
$\displaystyle+\gamma\pi_{m}+\alpha\pi_{m}\sum_{r=1}^{K}\rho_{r,m}\frac{|E^{\text{in}}_{r}(n)|+\delta|V_{r}(n)|}{|E(n)|+\delta|V(n)|}.$
(8)

Following a similar argument, we use the relationship in (4) to obtain that

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}$
$\displaystyle\left(|E^{\text{out}}_{m}(n+1)|-|E^{\text{out}}_{m}(n)|\right)$
$\displaystyle=\gamma\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}+\alpha\rho_{m\bullet}\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}$
$\displaystyle+\alpha\pi_{m}+\gamma\pi_{m}\sum_{r=1}^{K}\rho_{m,r}\frac{|E^{\text{out}}_{r}(n)|+\delta|V_{r}(n)|}{|E(n)|+\delta|V(n)|}.$

In addition, summing over $m$ in (8) gives

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}$
$\displaystyle\left(|E(n+1)|-|E(n)|\right)$
$\displaystyle=1+\alpha\sum_{m=1}^{K}\rho_{m\bullet}\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}+\gamma\sum_{m=1}^{K}\rho_{\bullet
m}\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}.$ (9)

As remarked at the beginning of Section 2, the proof of Theorem 2.2 requires
showing the convergence of $\mathbb{E}[|E(n)|]/n$ and for this it suffices to
check $\mathbb{E}[|E^{\text{in}}_{m}(n)|]/n\to x_{m}$ and
$\mathbb{E}[|E^{\text{out}}_{m}(n)|]/n\to y_{m}$, for some
$x_{m},y_{m}\in[0,1]$, which is implied by the approximation
$\mathbb{E}^{\mathcal{G}_{n}}\left(|E^{\text{in}}_{m}(n+1)|-|E^{\text{in}}_{m}(n)|\right)\approx
x_{m}$, and
$\mathbb{E}^{\mathcal{G}_{n}}\left(|E^{\text{out}}_{m}(n+1)|-|E^{\text{out}}_{m}(n)|\right)\approx
y_{m}$. Therefore, we expect $x_{m},y_{m}$, $m=1,\ldots,K$, to be solutions to
the following system of equations:

$\displaystyle x_{m}$
$\displaystyle=\alpha\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\gamma\rho_{\bullet
m}\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}+\gamma\pi_{m}+\alpha\pi_{m}\sum_{r}\rho_{r,m}\frac{x_{r}+\delta\pi_{r}}{\sum_{r}x_{r}+\delta},$
(10) $\displaystyle y_{m}$
$\displaystyle=\gamma\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}+\alpha\rho_{m\bullet}\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\alpha\pi_{m}+\gamma\pi_{m}\sum_{r}\rho_{m,r}\frac{y_{r}+\delta\pi_{r}}{\sum_{r}y_{r}+\delta},$
(11)

for $m=1,\ldots,K$. Summing over $m$ in (10) and (11) shows
$\sum_{m=1}^{K}x_{m}=\sum_{m=1}^{K}y_{m}$ and from (10) we get

$\sum_{m=1}^{K}x_{m}\geq\sum_{m=1}^{K}\left(\alpha\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\gamma\pi_{m}\right)=1,$

and therefore $\sum_{m=1}^{K}y_{m}\geq 1$. Also, since

$\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\leq
1,\qquad\text{and}\qquad\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\leq
1,$

the right hand side of (10) is upper bounded by

$\alpha+\gamma\rho_{\bullet
m}+\gamma\pi_{m}+\alpha\pi_{m}\text{max}_{r}\rho_{r,m}\leq 2.$

Similarly, the right hand side of (11) is bounded by 2.

Lemma 2.1 shows that under some regularity conditions, the system of equations
(10) and (11) has a unique solution, $(x_{1},\ldots,x_{K},y_{1},\ldots,y_{K})$
in

$\mathcal{Z}:=\\{\bm{z}\in[0,2]^{2K}:\sum_{i=1}^{K}z_{i}\geq
1,\sum_{i=K+1}^{2K}z_{i}\geq 1\\}.$

The proof relies on the contraction mapping theorem [13, Theorem 1.2.2]. We
proceed by defining a function $\bm{f}:\mathcal{Z}\mapsto\mathcal{Z}$, using
the right-hand-side expressions in (10) and (11). Applying the mean value
theorem to $\bm{f}$, we deduce a sufficient condition for $\bm{f}$ to be a
contraction by deriving an upper bound for the 1-norm of the Jacobian matrix.
Due to the linked definition of $(x_{m})$ and $(y_{m})$ in (10) and (11), we
write the Jacobian matrix as a block matrix, and find the upper bounds block
by block.

######  Lemma 2.1.

Let $\odot$ and $\otimes$ denote the Hadamard and Kronecker products of
matrices, respectively. Denote the identity matrix in $\mathbb{R}^{K}$ as
$\bm{I}_{K}$, and use $\bm{1}_{K}$ to denote the vector in $\mathbb{R}^{K}$
with all entries being 1. Set $\bm{\pi}:=(\pi_{1},\ldots,\pi_{K})$, and the
reciprocity matrix, $\bm{\rho}$, whose $(i,j)$-th entry is $\rho_{i,j}$.
Consider the four matrices:

$\displaystyle
J^{*}({1,1})=\alpha\left(\bm{I}_{K}+(\bm{\pi}\bm{1}_{K}^{T})\odot\bm{\rho}^{T}\right),$
$\displaystyle\quad
J^{*}({1,2})=\gamma(\bm{\rho}^{T}\bm{\pi})\otimes\bm{1}_{K}^{T},$
$\displaystyle J^{*}({2,1})=\alpha(\bm{\rho}\bm{\pi})\otimes\bm{1}_{K}^{T},$
$\displaystyle\quad
J^{*}({2,2})=\gamma\left(\bm{I}_{K}+(\bm{\pi}\bm{1}_{K}^{T})\odot\bm{\rho}\right),$

and the block matrix

$J^{*}:=\begin{bmatrix}J^{*}({1,1})&J^{*}({1,2})\\\
J^{*}({2,1})&J^{*}({2,2})\end{bmatrix}.$

Let $\|J^{*}\|_{1}$ denote the 1-norm of the matrix $J^{*}$. Then as long as

$\displaystyle\delta>\text{max}\\{\|J^{*}\|_{1}-1,0\\},$ (12)

the system of equations in (10) and (11) has a unique solution in
$\mathcal{Z}$.

The proof of Lemma 2.1 is deferred to Section 6.1. Continuing to assume (12)
holds, the next Theorem 2.2 gives the a.s. convergence of $|E(n)|/n$,
$|E^{\text{in}}(n)|/n$ and $|E^{\text{out}}(n)|/n$. Define a constant

$\displaystyle
C_{\delta}=\alpha\sum_{m=1}^{K}\rho_{m\bullet}\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\gamma\sum_{m=1}^{K}\rho_{\bullet
m}\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta},$ (13)

and a matrix

$\displaystyle\bm{H}:=\begin{bmatrix}\frac{\alpha}{1+\delta}\left(1+\bigvee_{m}\rho_{m\bullet}\right)&\frac{\gamma}{1+\delta}\bigvee_{m}\rho_{\bullet
m}&\frac{1}{1+\delta}\left(\alpha+C_{\delta}\right)\\\
\frac{\alpha}{1+\delta}\bigvee_{m}\rho_{m\bullet}&\frac{\gamma}{1+\delta}\left(1+\bigvee_{m}\rho_{\bullet
m}\right)&\frac{1}{1+\delta}\left(\gamma+C_{\delta}\right)\\\
\frac{\alpha}{1+\delta}\bigvee_{m}\rho_{m\bullet}&\frac{\gamma}{1+\delta}\bigvee_{m}\rho_{\bullet
m}&\frac{C_{\delta}}{1+\delta}\end{bmatrix}.$ (14)

The proof of Theorem 2.2 requires supposing $\lambda_{H}$, the largest
eigenvalue of $\bm{H}$, satisfies $\lambda_{H}<1$ to ensure the convergence of
$\mathbb{E}[|E(n)|]/n$, $\mathbb{E}[|E^{\text{in}}(n)|]/n$ and
$\mathbb{E}[|E^{\text{out}}(n)|]/n$.

Define
$\Delta^{\text{in}}_{m}(n):=|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|-n(x_{m}+\delta\pi_{m})$,
$\Delta^{\text{out}}_{m}(n):=|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|-n(y_{m}+\delta\pi_{m})$,
and $\Delta(n):=|E(n)|-n\sum_{m}x_{m}$. We will prove Theorem 2.2 by deriving
iterative upper bounds such that for the matrix $\bm{H}$ given in (14),

$\begin{bmatrix}\sum_{m}|\mathbb{E}(\Delta^{\text{in}}_{m}(n+1))|\\\
\sum_{m}|\mathbb{E}(\Delta^{\text{out}}_{m}(n+1))|\\\
|\mathbb{E}(\Delta(n))|\end{bmatrix}\leq\left(\bm{I}+\frac{1}{n}\bm{H}\right)\begin{bmatrix}\sum_{m}|\mathbb{E}(\Delta^{\text{in}}_{m}(n))|\\\
\sum_{m}|\mathbb{E}(\Delta^{\text{out}}_{m}(n))|\\\
|\mathbb{E}(\Delta(n))|\end{bmatrix}.$

Hence, provided that all elements in $\bm{H}$ are strictly positive, applying
the Perron-Frobenius theorem suggests restricting the largest eigenvalue of
$\bm{H}$ will be sufficient to control the growth of the expected edge counts.
In addition, the assumption of positive elements in $\bm{H}$ leads to the
requirement of $\alpha,\gamma>0$, $\bigvee_{m}\rho_{m\bullet}>0$, and
$\bigvee_{m}\rho_{\bullet m}>0$.

Recall we may interpret $\rho_{m\bullet}$ and $\rho_{\bullet m}$, them as the
likelihood of sending and attracting reciprocal edges for a user in group $m$,
respectively. The assumption $\bigvee_{m}\rho_{m\bullet}>0$ and
$\bigvee_{m}\rho_{\bullet m}>0$ requires that there exists at least one group
whose probability of generating or attracting reciprocal edges is strictly
positive.

######  Theorem 2.2.

Assume $\alpha,\gamma>0$, $\bigvee_{m}\rho_{m\bullet}>0$,
$\bigvee_{m}\rho_{\bullet m}>0$, $\delta>\|J^{*}\|_{F}-1$ and $\lambda_{H}<1$.
Suppose that $(x_{1},\ldots,x_{K},y_{1},\ldots,y_{K})$ is the solution to the
system of equations given in (10) and (11). We have:

  1. (i)

For $m=1,\ldots,K$,
$\frac{|E^{\text{in}}_{m}(n)|}{n}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}x_{m}$,
and
$\frac{|E^{\text{out}}_{m}(n)|}{n}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}y_{m}$.

  2. (ii)

In addition,
$\frac{|E(n)|}{n}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\sum_{m=1}^{K}x_{m}=\sum_{m=1}^{K}y_{m}$.

Note that since $(x_{m},y_{m})$ satisfies (10) and (11), the a.s. limit of
$|E(n)|/n$ depends on the value of $\delta$ as well. This is different from
the model with a homogeneous reciprocity level, where only the reciprocity
parameter determines the limit of $|E(n)|/n$.

The proof of Theorem 2.2 is in Section 6.2. Throughout the rest of the paper,
assume the following regularity conditions:

$\alpha,\gamma>0,\quad\delta>\|J^{*}\|_{1}-1,\quad\bigvee_{m}\rho_{m\bullet}>0,\,\bigvee_{m}\rho_{\bullet
m}>0,\quad\lambda_{H}<1.$ (15)

##  3 Growth of Degree Counts

We next focus on the asymptotic behavior of the joint in- and out-degree
counts:

$N_{k,l}(n):=\sum_{v\in
V(n)}\bm{1}_{\left\\{(D^{\text{in}}_{v}(n),D^{\text{out}}_{v}(n))=(k,l)\right\\}}.$

We extend the homogeneous reciprocity techniques in [26, 6] to obtain the
convergence of $N_{k,l}(n)/n$ under heterogeneous reciprocity. The way forward
is via embedding the in- and out-degree sequences in a family of multi-type
Markov branching processes with immigration (MBI processes).

###  3.1 Markov Branching with Immigration

To pursue count asymptotics using embedding, we need a linked family of two-
type Markov branching processes ([3]) where each process is a Markov Branching
with Immigration (MBI) process. See Section 2.1.1 of [26]. For each
$m\in\\{1,\ldots,K\\}$,

$\\{\bm{\xi}_{\delta}(t,m)\equiv\bigl{(}\xi^{(1)}_{\delta}(t,m),\xi^{(2)}_{\delta}(t,m)\bigr{)}:t\geq
0\\},$

is an MBI process whose branching structure depends on $m$; this will be
specified more precisely in (16) and (17). In the rest of the paper, we refer
to $\bm{\xi}_{\delta}(\cdot,m)$ as a _MBI process with group label_ $m$. The
process $\bm{\xi}_{\delta}(\cdot,m)$ is designed to mimic evolution of in- and
out-degrees of a fixed node with communication type $m$. The general setup of
continuous-time multitype branching processes without immigration is reviewed
in [2, Chapter V] and discussions on the MBI process are included, for
instance, in [19, 26].

We assume life time parameters of $\xi^{(1)}_{\delta}(\cdot,m)$ and
$\xi^{(2)}_{\delta}(\cdot,m)$ to be $\alpha$ and $\gamma\equiv 1-\alpha$,
respectively. For $\bm{s}=(s_{1},s_{2})\in[0,1]^{2}$, the branching structure
of $\bm{\xi}_{\delta}(\cdot,m)$ is specified through offspring generating
functions:

$\displaystyle f^{(1)}_{v}(\bm{s},m)$
$\displaystyle=\left(1-\rho_{m\bullet}\right)s_{1}^{2}+\rho_{m\bullet}s_{1}^{2}s_{2},$
(16) $\displaystyle f^{(2)}_{v}(\bm{s},m)$
$\displaystyle=\left(1-\rho_{\bullet m}\right)s_{2}^{2}+\rho_{\bullet
m}s_{1}s_{2}^{2},$ (17)

According to (16), when a group-$m$ type I particle’s lifetime ends, with
probability $1-\rho_{m\bullet}$, it splits into two group-$m$ type I
particles, increasing the total number of group-$m$ type I particles by 1, and
with probability $\rho_{m\bullet}$, it splits into two group-$m$ type I
particles and one group-$m$ type II particle. This last eventuality increases
the total numbers of group-$m$ type I and II particles both by 1. Immigration
events arrive following a homogeneous Poisson process with rate $\delta>0$.
When an immigration event happens, $\bm{\xi}_{\delta}(\cdot,m)$ is incremented
by $(1,0)$, $(0,1)$, or $(1,1)$ following the distribution

$\displaystyle
p_{0}(\bm{x},m)=\left(\alpha(1-\rho_{m\bullet})\right)^{\bm{1}_{\\{\bm{x}=(1,0)\\}}}\left(\gamma(1-\rho_{\bullet
m})\right)^{\bm{1}_{\\{\bm{x}=(0,1)\\}}}\left(\alpha\rho_{m\bullet}+\gamma\rho_{\bullet
m}\right)^{\bm{1}_{\\{\bm{x}=(1,1)\\}}}.$ (18)

All immigrants are of group $m$, and have the same branching structure as
given in (16) and (17). Following the discussion in Chapter V.7.2 of [2], we
use (16)–(17) to obtain a matrix that helps specify the branching structure of
$\bm{\xi}_{\delta}(\cdot,m)$:

$A_{m}=\begin{bmatrix}\alpha&\alpha\rho_{m\bullet}\\\ \gamma\rho_{\bullet
m}&\gamma\end{bmatrix},\qquad m\in\\{1,\ldots,K\\}.$ (19)

Applying the Perron-Frobenius theorem, $A_{m}$ has a largest positive
eigenvalue with multiplicity 1, i.e.

$\displaystyle\lambda_{m}=\frac{1}{2}\left(1+\sqrt{(\alpha-\gamma)^{2}+4\alpha\gamma\rho_{m\bullet}\rho_{\bullet
m}}\right)=:\frac{1}{2}\left(1+\sqrt{D_{0}(m)}\right).$ (20)

Also, the smaller eigenvalue of $A_{m}$ is

$\lambda^{\prime}_{m}=\frac{1}{2}\left(1-\sqrt{D_{0}(m)}\right),$

which plays an important role in the discussion of hidden regular variation in
Theorem 4.2. Let $\bm{v}{(m)}\equiv(v^{(1)}{(m)},v^{(2)}{(m)})$,
$\bm{u}{(m)}\equiv(u^{(1)}{(m)},u^{(2)}{(m)})$ be the left and right
eigenvectors associated with $\lambda_{m}$ respectively, with all coordinates
strictly positive, and $\bm{u}(m)^{T}\bm{1}=1$, $\bm{u}(m)^{T}\bm{v}(m)=1$.
Applying Theorem 1 in [26] gives that there exists some finite positive random
variable $Z(m)$ such that

$\displaystyle
e^{-\lambda_{m}t}\bm{\xi}_{\delta}(t,m)\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}Z(m)\bm{v}(m),\qquad
m\in\\{1,\dots,K\\}.$ (21)

The PA model with heterogeneous reciprocity assumes that the communication
group of a fixed node is determined at the node’s creation by flipping a
$K$-sided coin and a corresponding feature must also be present in the
embedding process. Let $L$ be a random variable with

$\displaystyle\mathbb{P}(L=m)=\pi_{m},\qquad m=1,\ldots,K,$ (22)

and for $t\geq 0$, the generating function of $\bm{\xi}_{\delta}(t,L)$ is

$\mathbb{E}\left(s_{1}^{\xi^{(1)}_{\delta}(t,L)}s_{2}^{\xi^{(2)}_{v,\delta}(t,L)}\right)=\sum_{m=1}^{K}\pi_{m}\mathbb{E}\left(s_{1}^{\xi^{(1)}_{\delta}(t,m)}s_{2}^{\xi^{(2)}_{\delta}(t,m)}\right),\qquad
s_{1},s_{2}\in[0,1].$

Imagine $L$ assigns a label to a MBI which remains unchanged throughout the
evolution of $\bm{\xi}_{\delta}(\cdot,\cdot)$. The initial value,
$\bm{\xi}_{\delta}(0,L)$, is a random vector on $\\{(0,1),(1,0),(1,1)\\}$,
whose distribution depends on the label $L$ and will be specified during the
embedding process. Conditioning on $(\bm{\xi}_{\delta}(0,L),L)$,
$\bm{\xi}_{\delta}(\cdot,L)$ behaves like a MBI process with fixed group label
$L$. In the sequel, we refer to $\\{\bm{\xi}_{\delta}(t,L):t\geq 0\\}$ as a
two-type MBI process with random group label.

###  3.2 Embedding

Let $\\{L_{v}:v\geq 1\\}$ be random variables with common pmf’s as in (22). We
now explain how to embed the degree sequence
$\\{(D^{\text{in}}_{v}(n),D^{\text{out}}_{v}(n)):v\in V(n)\\}_{n\geq 1}$ into
a family of linked two-type MBI processes with random group labels,
$\\{\bm{\xi}_{v,\delta}(t,L_{v}):t\geq 0\\}_{v\geq 1}$. We assume
$\\{\bm{\xi}_{v,\delta}(\cdot,\cdot)\\}_{v\geq 1}$ have the same parameters
but possibly different initializations.

At $T_{0}=0$, flip a $K$-sided coin with outcome $L_{1}$. If $L_{1}=m_{1}$,
which occurs with probability $\pi_{m_{1}}$, initiate the MBI process
$\bm{\xi}_{1,\delta}(\cdot,m_{1})$ with $\bm{\xi}_{1,\delta}(0,m_{1})=(1,1)$.
Let $T_{1}$ be the first jump time of $\bm{\xi}_{1,\delta}(\cdot,m_{1})$.Then
for $t\geq 0$,

$\displaystyle\mathbb{P}\left(T_{1}>t,L_{1}=m_{1}\right)$
$\displaystyle=\pi_{m_{1}}\exp\left\\{-t\left(\alpha\left(\bm{\xi}^{(1)}_{1,\delta}(0,m_{1})+\delta\right)+\gamma\left(\bm{\xi}^{(2)}_{1,\delta}(0,m_{1})+\delta\right)\right)\right\\}$
$\displaystyle=\pi_{m_{1}}e^{-t(1+\delta)},$

and $\mathbb{P}(T_{1}>t)=e^{-t(1+\delta)}$.

At $T_{1}$, flip another $K$-sided coin with outcome $L_{2}$ and start the
process $\\{\bm{\xi}_{2,\delta}(t-T_{1},L_{2}):t\geq T_{1}\\}$. The initial
value of $\bm{\xi}_{2,\delta}(0,L_{2})$ is set depending on which of the
following four cases happens:

  1. (i)

If $L_{1}=m_{1}$ and at time $T_{1}$ the process
$\bm{\xi}_{1,\delta}(\cdot,{m_{1}})$ increases by $(1,0)$, then we have
$L_{2}=m_{2}$ with (conditional) probability
$\pi_{m_{2}}(1-\rho_{m_{1},m_{2}})/(1-\rho_{m_{1}\bullet})$, and set
$\bm{\xi}_{2,\delta}(0,m_{2})=(0,1)$.

  2. (ii)

If the $\bm{\xi}_{1,\delta}(\cdot,L_{1})$ with $L_{1}=m_{1}$ is increased by
$(0,1)$, then we have $L_{2}=m_{2}$ with (conditional) probability
$\pi_{m_{2}}(1-\rho_{m_{2},m_{1}})/(1-\rho_{\bullet m_{1}})$, and set
$\bm{\xi}_{2,\delta}(0,m_{2})=(1,0)$.

  3. (iii)

If one type I particle in $\bm{\xi}_{1,\delta}(\cdot,L_{1})$ with
$L_{1}=m_{1}$ splits into 2 type I and 1 type II particles at $T_{1}$, then we
have $L_{2}=m_{2}$ with (conditional) probability
$\pi_{m_{2}}\rho_{m_{1},m_{2}}/\rho_{m_{1}\bullet}$, and set
$\bm{\xi}_{2,\delta}(0,m_{2})=(1,1)$.

  4. (iv)

If one type II particle in $\bm{\xi}_{1,\delta}(\cdot,L_{1})$ with
$L_{1}=m_{1}$ splits into 1 type I and 2 type II particles at $T_{1}$, then we
have $L_{2}=m_{2}$ with (conditional) probability
$\pi_{m_{2}}\rho_{m_{2},m_{1}}/\rho_{\bullet m_{1}}$, and set
$\bm{\xi}_{2,\delta}(0,m_{2})=(1,1)$.

Therefore, we see that

$\displaystyle\mathbb{P}(L_{1}=m_{1},L_{2}=m_{2},T_{2}-T_{1}>t)$
$\displaystyle=\pi_{m_{1}}e^{-t(1+\delta)}\left(\alpha(1-\rho_{m_{1}\bullet})\frac{\pi_{m_{2}}(1-\rho_{m_{1},m_{2}})}{1-\rho_{m_{1}\bullet}}+\alpha\rho_{m_{1}\bullet}\frac{\pi_{m_{2}}\rho_{m_{1},m_{2}}}{\rho_{m_{1}\bullet}}\right.$
$\displaystyle\left.\quad+\gamma(1-\rho_{\bullet
m_{1}})\frac{\pi_{m_{2}}(1-\rho_{m_{2},m_{1}})}{1-\rho_{\bullet
m_{1}}}+\gamma\rho_{\bullet
m_{1}}\frac{\pi_{m_{2}}\rho_{m_{2},m_{1}}}{\rho_{\bullet m_{1}}}\right)$
$\displaystyle=\pi_{m_{1}}\pi_{m_{2}}e^{-t(1+\delta)},$

which shows that $(L_{1},L_{2})$ are independent. Also,
$\bm{\xi}_{2,\delta}(0,L_{2})$ is a 2-dimensional random vector whose
generating function satisfies

$\displaystyle\mathbb{E}\left(s_{1}^{\xi^{(1)}_{2,\delta}(0,L_{2})}s_{2}^{\xi^{(2)}_{2,\delta}(0,L_{2})}\middle|L_{1}\right)$
$\displaystyle=\alpha(1-\rho_{L_{1}\bullet})s_{2}+\gamma(1-\rho_{\bullet
L_{1}})s_{1}+(\alpha\rho_{L_{1}\bullet}+\gamma\rho_{\bullet
L_{1}})s_{1}s_{2},$

so that for
$\rho_{0}:=\sum_{m}\pi_{m}\rho_{m\bullet}=\sum_{m}\pi_{m}\rho_{\bullet m}$,

$\displaystyle\mathbb{E}\left(s_{1}^{\xi^{(1)}_{2,\delta}(0,L_{2})}s_{2}^{\xi^{(2)}_{2,\delta}(0,L_{2})}\right)$
$\displaystyle=\sum_{m=1}^{K}\pi_{m}\left(\alpha(1-\rho_{m\bullet})s_{2}+\gamma(1-\rho_{\bullet
m})s_{1}+(\alpha\rho_{m\bullet}+\gamma\rho_{\bullet m})s_{1}s_{2}\right)$
$\displaystyle=\alpha(1-\rho_{0})s_{2}+\gamma(1-\rho_{0})s_{1}+\rho_{0}s_{1}s_{2},\qquad
s_{1},s_{2}\in[0,1].$ (23)

Define
$R_{1}:=\bm{1}_{\\{\bm{\xi}_{1,\delta}(T_{1},L_{1})=\bm{\xi}_{1,\delta}(0,L_{1})+(1,1)\\}}$,
so that

$\mathbb{P}(R_{1}=1|L_{1}=m)=\alpha\rho_{m\bullet}+\gamma\rho_{\bullet
m}=1-\mathbb{P}(R_{1}=0|L_{1}=m),$

and $\mathbb{P}(R_{1}=1)=\rho_{0}=1-\mathbb{P}(R_{1}=0)$. Set

$\displaystyle\mathcal{F}_{T_{1}}:=$
$\displaystyle\sigma\left(\\{L_{k}\\}_{k=1,2};\left\\{\bm{\xi}_{k,\delta}(t-T_{k-1},L_{k}):t\in[T_{k-1},T_{1}]\right\\}_{k=1,2}\right)$
$\displaystyle=$
$\displaystyle{\sigma\left(\\{L_{k}\\}_{k=1,2};\bm{\xi}_{1,\delta}(t,L_{1}),t\in[0,T_{1}];\bm{\xi}_{2,\delta}(0,L_{2})\right)},$

we see that $R_{1}$ is $\mathcal{F}_{T_{1}}$-measurable.

In general, for $n\geq 1$, suppose that we have initiated $n+1$ MBI processes
with group labels $\\{L_{k}:1\leq k\leq n+1\\}$ at time $T_{n}$,

$\displaystyle\\{\bm{\xi}_{k,\delta}(t-T_{k-1},L_{k}):t\geq T_{k-1}\\}_{1\leq
k\leq n+1}.$ (24)

Define $T_{n+1}$ as the first time when one of the processes in (24) jumps,
and let $J_{n+1}$ be the index of the process that jumps at $T_{n+1}$. Define
the $\sigma$-field

$\mathcal{F}_{T_{n}}:=\sigma\left(\\{L_{k}\\}_{k=1}^{n+1};\left\\{\bm{\xi}_{k,\delta}(t-T_{k-1},L_{k}):t\in[T_{k-1},T_{n}]\right\\}_{1\leq
k\leq n+1}\right),$

and

$R_{n+1}:=\bm{1}_{\left\\{\bm{\xi}_{J_{n+1},\delta}(T_{n+1}-T_{J_{n+1}-1},L_{J_{n+1}})=\bm{\xi}_{J_{n+1},\delta}(T_{n}-T_{J_{n+1}-1},L_{J_{n+1}})+(1,1)\right\\}}.$

Then we see that $\\{R_{k}:1\leq k\leq n\\}$ are
$\mathcal{F}_{T_{n}}$-measurable.

At $T_{n+1}$, we initiate the MBI process with group label $L_{n+2}$,
$\\{\bm{\xi}_{n+2,\delta}(t-T_{n+1},L_{n+2}):t\geq T_{n+1}\\}$, and one of the
following four cases happens:

  1. (i)

If the $\bm{\xi}_{J_{n+1},\delta}(\cdot,L_{J_{n+1}})$ with
$L_{J_{n+1}}=m_{n+1}$ is increased by $(1,0)$, then we have $L_{n+2}=m_{n+2}$
with (conditional) probability
$\pi_{m_{n+2}}(1-\rho_{m_{n+1},m_{n+2}})/(1-\rho_{m_{n+1}\bullet})$, and set
$\bm{\xi}_{n+2,\delta}(0,m_{n+2})=(0,1)$.

  2. (ii)

If the $\bm{\xi}_{J_{n+1},\delta}(\cdot,L_{J_{n+1}})$ with
$L_{J_{n+1}}=m_{n+1}$ is increased by $(0,1)$, then we have $L_{n+2}=m_{n+2}$
with (conditional) probability
$\pi_{m_{n+2}}(1-\rho_{m_{n+2},m_{n+1}})/(1-\rho_{\bullet m_{n+1}})$, and set
$\bm{\xi}_{n+2,\delta}(0,m_{n+2})=(1,0)$.

  3. (iii)

If one type I particle in $\bm{\xi}_{J_{n+1},\delta}(\cdot,L_{J_{n+1}})$ with
$L_{J_{n+1}}=m_{n+1}$ splits into 2 type I and 1 type II particles at
$T_{n+1}$, then we have $L_{n+2}=m_{n+2}$ with (conditional) probability
$\pi_{m_{n+2}}\rho_{m_{n+1},m_{n+2}}/\rho_{m_{n+1}\bullet}$, and set
$\bm{\xi}_{n+2,\delta}(0,m_{n+2})=(1,1)$.

  4. (iv)

If one type II particle in $\bm{\xi}_{J_{n+1},\delta}(\cdot,L_{J_{n+1}})$ with
$L_{J_{n+1}}=m_{n+1}$ splits into 1 type I and 2 type II particles at
$T_{n+1}$, then we have $L_{n+2}=m_{n+2}$ with (conditional) probability
$\pi_{m_{n+2}}\rho_{m_{n+2},m_{n+1}}/\rho_{\bullet m_{n+1}}$, and set
$\bm{\xi}_{n+2,\delta}(0,m_{n+2})=(1,1)$.

Since
$\sum_{k=1}^{n+1}\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})=\sum_{k=1}^{n+1}\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1})=n+1+\sum_{k=1}^{n}R_{k}$,
we then have

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(L_{n+2}=m,T_{n+1}-T_{n}>t\right)$
$\displaystyle=\sum_{r=1}^{K}\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(L_{n+2}=m,L_{J_{n+1}}=r,T_{n+1}-T_{n}>t\right)$
$\displaystyle=\exp\left\\{-t\left((1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}\right)\right\\}$
$\displaystyle\times\sum_{r=1}^{K}\left[\frac{\sum_{k=1}^{n+1}(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta)\bm{1}_{\\{L_{k}=r\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\left(\alpha(1-\rho_{r\bullet})\frac{\pi_{m}(1-\rho_{r,m})}{1-\rho_{r\bullet}}+\alpha\rho_{r\bullet}\frac{\pi_{m}\rho_{r,m}}{\rho_{r\bullet}}\right)\right.$
$\displaystyle\left.\quad+\frac{\sum_{k=1}^{n+1}(\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta)\bm{1}_{\\{L_{k}=r\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\left(\gamma(1-\rho_{\bullet
r})\frac{\pi_{m}(1-\rho_{m,r})}{1-\rho_{\bullet r}}+\gamma\rho_{\bullet
r}\frac{\pi_{m}\rho_{m,r}}{\rho_{\bullet r}}\right)\right]$
$\displaystyle=\pi_{m}\exp\left\\{-t\left((1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}\right)\right\\}$
$\displaystyle=\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(L_{n+2}=m\right)\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(T_{n+1}-T_{n}>t\right).$

Hence, $L_{n+2}$ and $T_{n+1}-T_{n}$ are independent under
$\mathbb{P}^{\mathcal{F}_{T_{n}}}$, and since

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(L_{n+2}=m\right)=\pi_{m}=\mathbb{P}\left(L_{n+2}=m\right),$
(25)

$L_{n+2}$ is independent from $\\{L_{k}:1\leq k\leq n+1\\}$.

In addition, from the four cases listed above, we also see that

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}(\bm{\xi}_{n+2,\delta}(0,L_{n+2})=(0,1),L_{J_{n+1}}=m,L_{n+2}=r)$
$\displaystyle=\frac{\alpha\pi_{r}(1-\rho_{m,r})\sum_{k=1}^{n+1}\left(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{\alpha\sum_{k=1}^{n+1}\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\gamma\sum_{k=1}^{n-1}\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+(n+1)\delta}$
$\displaystyle=\alpha\pi_{r}(1-\rho_{m,r})\frac{\sum_{k=1}^{n+1}\left(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}.$

Similarly, we have

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}$
$\displaystyle(\bm{\xi}_{n+2,\delta}(0,L_{n+2})=(1,0),L_{J_{n+1}}=m,L_{n+2}=r)$
$\displaystyle=\gamma\pi_{r}(1-\rho_{r,m})\frac{\sum_{k=1}^{n+1}\left(\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}},$

and

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}$
$\displaystyle(\bm{\xi}_{n+2,\delta}(0,L_{n+2})=(1,1),L_{J_{n+1}}=m,L_{n+2}=r)$
$\displaystyle=$
$\displaystyle\frac{\alpha\pi_{r}\rho_{m,r}\sum_{k=1}^{n+1}\left(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}$
$\displaystyle+\frac{\gamma\pi_{r}\rho_{r,m}\sum_{k=1}^{n+1}\left(\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}.$

Therefore, under $\mathbb{P}^{\mathcal{F}_{T_{n}}}$,
$\bm{\xi}_{n+2,\delta}(0,L_{n+2})$ is a random vector following the
distribution

$\displaystyle p_{n+2,0}(\bm{x})$
$\displaystyle=\left(\sum_{m=1}^{K}\frac{\alpha(1-\rho_{m\bullet})\sum_{k=1}^{n+1}\left(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\right)^{\bm{1}_{\\{\bm{x}=(0,1)\\}}}$
$\displaystyle\times\left(\sum_{m=1}^{K}\frac{\gamma(1-\rho_{\bullet
m})\sum_{k=1}^{n+1}\left(\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\right)^{\bm{1}_{\\{\bm{x}=(1,0)\\}}}$
$\displaystyle\times\left(\sum_{m=1}^{K}\left(\frac{\alpha\rho_{m\bullet}\sum_{k=1}^{n+1}\left(\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\right.\right.$
$\displaystyle\left.\left.\quad+\frac{\gamma\rho_{\bullet
m}\sum_{k=1}^{n+1}\left(\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\delta\right)\bm{1}_{\\{L_{k}=m\\}}}{(1+\delta)(n+1)+\sum_{k=1}^{n}R_{k}}\right)\right)^{\bm{1}_{\\{\bm{x}=(1,1)\\}}}.$

Meanwhile, since $\\{L_{k}:1\leq k\leq n+1\\}$ are
$\mathcal{F}_{T_{n}}$-measurable, we then have

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(R_{n+1}=1,J_{n+1}=w,T_{n+1}-T_{n}>t\right)$
$\displaystyle=\frac{\alpha\rho_{L_{w}\bullet}\left(\xi^{(1)}_{w,\delta}(T_{n}-T_{w-1},L_{w})+\delta\right)+\gamma\rho_{\bullet
L_{w}}\left(\xi^{(2)}_{w,\delta}(T_{n}-T_{w-1},L_{w})+\delta\right)}{\alpha\sum_{k=1}^{n+1}\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\gamma\sum_{k=1}^{n+1}\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+(n+1)\delta}$
$\displaystyle\quad\times\exp\left\\{-t\left(\alpha\sum_{k=1}^{n+1}\xi^{(1)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+\gamma\sum_{k=1}^{n+1}\xi^{(2)}_{k,\delta}(T_{n}-T_{k-1},L_{k})+(n+1)\delta\right)\right\\}$
$\displaystyle=\frac{\alpha\rho_{L_{w}\bullet}\left(\xi^{(1)}_{w,\delta}(T_{n}-T_{w-1},L_{w})+\delta\right)+\gamma\rho_{\bullet
L_{w}}\left(\xi^{(2)}_{w,\delta}(T_{n}-T_{w-1},L_{w})+\delta\right)}{(n+1)(1+\delta)+\sum_{k=1}^{n}R_{k}}e^{-t((n+1)(1+\delta)+\sum_{k=1}^{n}R_{k})}$
$\displaystyle=\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(R_{n+1}=1,J_{n+1}=w\right)\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(T_{n+1}-T_{n}>t\right).$
(26)

Therefore, $\bigl{(}R_{n+1},J_{n+1}\bigr{)}$ is independent from
$T_{n+1}-T_{n}$ under $\mathbb{P}^{\mathcal{F}_{T_{n}}}$.

Define for $n\geq 0$,

${\bm{\xi}}^{*}_{\delta}\left(T_{n},L_{[n+1]}\right):=\left(\bm{\xi}_{1,\delta}(T_{n},L_{1}),\bm{\xi}_{2,\delta}(T_{n}-T_{1},L_{2}),\ldots,\bm{\xi}_{n+1,\delta}(0,L_{n+1}),(0,0),\ldots\right),$

then the embedding framework just described shows that
$\\{{\bm{\xi}}^{*}_{\delta}\left(T_{n},L_{[n+1]}\right):n\geq 0\\}$ is Markov
with state space $\left(\mathbb{N}^{2}\right)^{\infty}$. The next theorem
gives the embedding of in- and out-degree sequences in the PA model with
heterogeneous reciprocity into a linked system of delayed MBI processes.

######  Theorem 3.1.

In $\left(\mathbb{N}^{2}\right)^{\infty}$, define the in- and out-degree
sequences as

$\bm{D}(n):=\left(\bigl{(}D^{\text{in}}_{1}(n),D^{\text{out}}_{1}(n)\bigr{)},\ldots,\bigl{(}D^{\text{in}}_{n+1}(n),D^{\text{out}}_{n+1}(n)\bigr{)},(0,0),\ldots\right)\qquad
n\geq 0.$

Then for $\\{T_{k}:k\geq 0\\}$ and
$\\{\bm{\xi}_{k,\delta}(t-T_{k-1},L_{k}):t\geq T_{k-1}\\}_{k\geq 1}$
constructed above, we have that in
$\left(\left(\mathbb{N}^{2}\right)^{\infty}\right)^{\infty}$,

$\displaystyle\bigl{\\{}\bm{D}(n):n\geq 0\bigr{\\}}\stackrel{{\scriptstyle
d}}{{=}}\left\\{{\bm{\xi}}^{*}_{\delta}\left(T_{n},L_{[n+1]}\right):n\geq
0\right\\}.$

###### Proof.

By the model description in Section 1.1, $\\{\bm{D}(n):n\geq 0\\}$ is
Markovian on $\left(\mathbb{N}^{2}\right)^{\infty}$, so it suffices to check
the transition probability from $\bm{D}(n)$ to $\bm{D}(n+1)$ agrees with that
from ${\bm{\xi}}^{*}_{\delta}\left(T_{n},L_{[n+1]}\right)$ to
${\bm{\xi}}^{*}_{\delta}\left(T_{n+1},L_{[n+2]}\right)$. Write

$\displaystyle\bm{e}_{v}^{(1)}$
$\displaystyle:=\left(\bigl{(}0,0\bigr{)},\ldots,\bigl{(}0,0\bigr{)},\underbrace{\bigl{(}1,0\bigr{)}}_{\text{$v$-th
entry}},\bigl{(}0,0\bigr{)},\ldots\right),$ $\displaystyle\bm{e}_{v}^{(2)}$
$\displaystyle:=\left(\bigl{(}0,0\bigr{)},\ldots,\bigl{(}0,0\bigr{)},\underbrace{\bigl{(}0,1\bigr{)}}_{\text{$v$-th
entry}},\bigl{(}0,0\bigr{)},\ldots\right),$ $\displaystyle\bm{e}_{v}^{(3)}$
$\displaystyle:=\left(\bigl{(}0,0\bigr{)},\ldots,\bigl{(}0,0\bigr{)},\underbrace{\bigl{(}1,1\bigr{)}}_{\text{$v$-th
entry}},\bigl{(}0,0\bigr{)},\ldots\right),$

and we have

$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\left(\bm{D}(n+1)=\bm{D}(n)+\bm{e}_{v}^{(1)}+\bm{e}_{|V(n)|+1}^{(2)}\right)$
$\displaystyle=\frac{\alpha(D^{\text{in}}_{v}(n)+\delta)}{|E(n)|+\delta|V(n)|}\sum_{m=1}^{K}(1-\rho_{W_{v},m}\pi_{m})$
$\displaystyle=\frac{\alpha(D^{\text{in}}_{v}(n)+\delta)}{|E(n)|+\delta|V(n)|}(1-\rho_{W_{v}\bullet});$
(27) similarly,
$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\left(\bm{D}(n+1)=\bm{D}(n)+\bm{e}_{v}^{(2)}+\bm{e}_{|V(n)|+1}^{(1)}\right)$
$\displaystyle=\frac{\gamma(D^{\text{out}}_{v}(n)+\delta)}{|E(n)|+\delta|V(n)|}(1-\rho_{\bullet
W_{v}}),$ (28)
$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}\left(\bm{D}(n+1)=\bm{D}(n)+\bm{e}_{v}^{(3)}+\bm{e}_{|V(n)|+1}^{(3)}\right)$
$\displaystyle=\frac{\alpha\rho_{W_{v}\bullet}(D^{\text{in}}_{v}(n)+\delta)+\gamma\rho_{\bullet
W_{v}}(D^{\text{out}}_{v}(n)+\delta)}{|E(n)|+\delta|V(n)|}.$ (29)

Note that $|V(n)|=n+1$ for all $n\geq 0$, and from (26), we have

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}$
$\displaystyle\left({\bm{\xi}}^{*}_{\delta}\left(T_{n+1},L_{[n+2]}\right)={\bm{\xi}}^{*}_{\delta}\left(T_{n},L_{[n+1]}\right)+\bm{e}_{v}^{(3)}+\bm{e}_{|V(n)|+1}^{(3)}\right)$
$\displaystyle=\mathbb{P}^{\mathcal{F}_{T_{n}}}\left(R_{n+1}=1,J_{n+1}=v\right)$
$\displaystyle=\frac{\alpha\rho_{L_{v}\bullet}\left(\xi^{(1)}_{v,\delta}(T_{n}-T_{v-1},L_{v})+\delta\right)+\gamma\rho_{\bullet
L_{v}}\left(\xi^{(2)}_{v,\delta}(T_{n}-T_{v-1},L_{v})+\delta\right)}{(n+1)(1+\delta)+\sum_{k=1}^{n}R_{k}}.$
(30)

So it remains to check whether $\sum_{k=1}^{n}R_{k}$ has the same distribution
as $|E(n)|-(n+1)$.

Again, applying (26) gives that

$\displaystyle\mathbb{P}^{\mathcal{F}_{T_{n}}}$
$\displaystyle\left(R_{n+1}=1\right)$
$\displaystyle=\sum_{v=1}^{n}\frac{\alpha\rho_{L_{v}\bullet}\left(\xi^{(1)}_{v,\delta}(T_{n}-T_{v-1},L_{v})+\delta\right)+\gamma\rho_{\bullet
L_{v}}\left(\xi^{(2)}_{v,\delta}(T_{n}-T_{v-1},L_{v})+\delta\right)}{(n+1)(1+\delta)+\sum_{k=1}^{n}R_{k}}.$

Also, we obtain from (8) that $|E(n)|-(n+1)$ satisfies

$\displaystyle\mathbb{P}^{\mathcal{G}_{n}}$
$\displaystyle\left(|E(n+1)|-|E(n)|-1=1\right)$
$\displaystyle=\alpha\sum_{m=1}^{K}\rho_{m\bullet}\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta(n+1)}+\gamma\sum_{m=1}^{K}\rho_{\bullet
m}\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta(n+1)}$
$\displaystyle=\sum_{v=1}^{n}\frac{\alpha\rho_{W_{v}\bullet}\left(D^{\text{in}}_{v}(n)+\delta\right)+\gamma\rho_{\bullet
W_{v}}\left(D^{\text{out}}_{v}(n)+\delta\right)}{|E(n)|+\delta(n+1)}.$ (31)

By (25), we see that $\\{L_{v}:v\geq 1\\}$ are iid random variables with
$\mathbb{P}(L_{v}=m)=\pi_{m}$, agreeing with the distributional property of
$\\{W_{v}:v\geq 1\\}$. Therefore, we conclude from (3.2) that
$\sum_{k=1}^{n}R_{k}$ has the same distribution as $|E(n)|-(n+1)$, which
implies the agreement between the transition probabilities in (29) and (30). ∎

###  3.3 Asymptotic Growth of Empirical Degree Frequencies

We now use the embedding results in Theorem 3.1 to prove the convergence of
$N_{k,l}(n)/n$.

######  Theorem 3.2.

Let $\\{\widetilde{\bm{\xi}}_{\delta}(t,L^{*}):t\geq 0\\}$ be a MBI process
with random group label $L^{*}$, where $L^{*}$ satisfies
$\mathbb{P}(L^{*}=m)=\pi_{m}$, $m=1,\ldots,K$. Suppose that the regularity
condition (15) holds, and the initialization,
$\widetilde{\bm{\xi}}_{\delta}(0,L^{*})$, satisfies that for $s_{i}\in[0,1]$,
$i=1,2$,

$\displaystyle\mathbb{E}$
$\displaystyle\left(s_{1}^{\widetilde{\xi}^{(1)}_{\delta}(0,L^{*})}s_{2}^{\widetilde{\xi}^{(2)}_{\delta}(0,L^{*})}\right)$
$\displaystyle=\sum_{r=1}^{K}\pi_{r}\left[\alpha\left(1-\sum_{m=1}^{K}\rho_{m,r}\frac{x_{m}+\delta\pi_{m}}{\sum_{s}x_{s}+\delta}\right)s_{2}+\gamma\left(1-\sum_{m=1}^{K}\rho_{r,m}\frac{y_{m}+\delta\pi_{m}}{\sum_{s}y_{s}+\delta}\right)s_{1}\right.$
$\displaystyle\left.+\left(\alpha\sum_{m=1}^{K}\rho_{m,r}\frac{x_{m}+\delta\pi_{m}}{\sum_{s}x_{s}+\delta}+\gamma\sum_{m=1}^{K}\rho_{r,m}\frac{y_{m}+\delta\pi_{m}}{\sum_{s}y_{s}+\delta}\right)s_{1}s_{2}\right].$
(32)

For $L^{*}=m$, the branching structure of
$\\{\widetilde{\bm{\xi}}_{\delta}(t,m):t\geq 0\\}$ is given by $A_{m}$ (cf.
(19)). Write

$\rho^{*}:=\sum_{m=1}^{K}x_{m}-1=\sum_{m=1}^{K}y_{m}-1>0\text{ and
}c^{*}:=1+\rho^{*}+\delta,$ (33)

then as $n\to\infty$, we have for $k,l\geq 0$,

$\displaystyle\frac{N_{k,l}(n)}{n}$ $\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\sum_{m=1}^{K}\pi_{m}\int_{0}^{\infty}c^{*}e^{-c^{*}t}\mathbb{P}\left(\widetilde{\bm{\xi}}_{\delta}(t,m)=(k,l)\right)\differential
t.$ (34)

Let $T^{*}$ be an exponential random variable with rate $c^{*}$, independent
from $\\{\widetilde{\bm{\xi}}_{\delta}(\cdot,m):m=1,\ldots,K\\}$, then the
integral on the right hand side of (34) is

$\mathbb{P}\left(\widetilde{\bm{\xi}}_{\delta}(T^{*},m)=(k,l)\right),\qquad
m=1,\ldots,K,$

representing the limiting in- and out-degree frequencies for nodes of
communication group $m$. In the sequel, set
$(\mathcal{I}_{m},\mathcal{O}_{m})=\widetilde{\bm{\xi}}_{\delta}(T^{*},m)$ to
denote the limit random variables so that (34) becomes

$\displaystyle\frac{N_{k,l}(n)}{n}$ $\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\sum_{m=1}^{K}\pi_{m}\mathbb{P}\Bigl{[}\big{(}\mathcal{I}_{m},\mathcal{O}_{m}\bigr{)}=(k,l)\Bigr{]}.$
(35)

###### Proof.

By the embedding results in Theorem 3.1, we have

$\displaystyle\frac{N_{k,l}(n)}{n}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n}\sum_{w=2}^{n+1}\bm{1}_{\left\\{\bm{\xi}_{w,\delta}(T_{n}-T_{w-1},L_{w})=(k,l)\right\\}}+\frac{1}{n}\bm{1}_{\left\\{\bm{\xi}_{1,\delta}(T_{n},L_{1})=(k,l)\right\\}},$
(36)

where the second term goes to 0 a.s. as $n\to\infty$. Then we only need to
consider the limit of the first term in (36).

Let $\\{\widetilde{\bm{\xi}}_{w,\delta}(t,\widetilde{L}_{w}):t\geq 0\\}_{w\geq
1}$ be a sequence of iid MBI processes with random labels, which have the same
distributional properties as $\widetilde{\bm{\xi}}_{\delta}(\cdot,L^{*})$ and
satisfy the initialization condition in (32). Then we divide the first term in
(36) into different parts:

$\displaystyle\frac{1}{n}\sum_{w=2}^{n+1}\bm{1}_{\left\\{\bm{\xi}_{w,\delta}(T_{n}-T_{w-1},L_{w})=(k,l)\right\\}}$
$\displaystyle=\frac{1}{n}\sum_{w=2}^{n+1}\left(\bm{1}_{\left\\{\bm{\xi}_{w,\delta}\left(T_{n}-T_{w-1},L_{w}\right)=(k,l)\right\\}}-\bm{1}_{\left\\{\bm{\xi}_{w,\delta}\bigl{(}\log(n/w)/c^{*},L_{w}\bigr{)}=(k,l)\right\\}}\right)$
$\displaystyle\,+\frac{1}{n}\sum_{w=2}^{n+1}\left(\bm{1}_{\left\\{\bm{\xi}_{w,\delta}\bigl{(}\log(n/w)/c^{*}\bigr{)}=(k,l)\right\\}}-\mathbb{P}^{\mathcal{F}_{T_{w-1}}}\left[\bm{\xi}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),L_{w}\right)=(k,l)\right]\right)$
$\displaystyle\,+\frac{1}{n}\sum_{w=2}^{n+1}\left(\mathbb{P}^{\mathcal{F}_{T_{w-1}}}\left[\bm{\xi}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),L_{w}\right)=(k,l)\right]\right.$
$\displaystyle\left.\qquad\qquad\qquad-\mathbb{P}\left[\widetilde{\bm{\xi}}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),\widetilde{L}_{w}\right)=(k,l)\middle|\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right]\right)$
$\displaystyle\,+\left(\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{P}\left[\widetilde{\bm{\xi}}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),\widetilde{L}_{w}\right)=(k,l)\middle|\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right]\right.$
$\displaystyle\left.\qquad-\int_{0}^{1}\mathbb{P}\left[\widetilde{\bm{\xi}}_{\delta}\left(-\frac{1}{c^{*}}\log
t,L^{*}\right)=(k,l)\right]\differential t\right)$
$\displaystyle\,+\int_{0}^{1}\mathbb{P}\left[\widetilde{\bm{\xi}}_{\delta}\left(-\frac{1}{c^{*}}\log
t,L^{*}\right)=(k,l)\right]\differential t$
$\displaystyle=:A_{1}(n)+A_{2}(n)+A_{3}(n)+A_{4}(n)+A_{5}.$

Note that by a change of variable argument, $A_{5}$ is identical to the right
hand side of (34), and we now show that $A_{j}(n)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0$, for $j=1,2,3,4$.

For $A_{1}(n)$, we have

$\displaystyle\mathbb{E}|A_{1}(n)|$
$\displaystyle\leq\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{E}\left|\bm{1}_{\left\\{\xi^{(1)}_{w,\delta}(T_{n}-T_{w-1},L_{w})=k\right\\}}-\bm{1}_{\\{\xi^{(1)}_{w,\delta}\bigl{(}\log(n/w)/c^{*},L_{w}\bigr{)}=k\\}}\right|$
$\displaystyle\quad+\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{E}\left|\bm{1}_{\\{\xi^{(2)}_{w,\delta}(T_{n}-T_{w-1},L_{w})=l\\}}-\bm{1}_{\\{\xi^{(2)}_{w,\delta}\bigl{(}\log(n/w)/c^{*},L_{w}\bigr{)}=l\\}}\right|.$
(37)

Since both $\xi^{(1)}_{w,\delta}(\cdot,L_{w})$ and
$\xi^{(2)}_{w,\delta}(\cdot,L_{w})$ have finite number of jumps in any finite
time interval $[0,K]$ a.s., then applying Lemma 3.1 in [1] gives that for all
$K>0$ and $w\geq 2$,

$\lim_{\epsilon\downarrow
0}\sup_{t\in[0,K]}\mathbb{P}\left(\xi^{(i)}_{w,\delta}(t+\epsilon,L_{w})-\xi^{(i)}_{w,\delta}\bigl{(}(t-\epsilon)\wedge
0,L_{w})\bigr{)}\geq 1\right)=0,\qquad i=1,2.$

Also, we see from [1, Corollary 2.1(iii)] that for $\eta>0$,

$\sup_{n\eta\leq w\leq
n}\left|T_{n}-T_{w-1}-\frac{1}{c^{*}}\log(n/w)\right|\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}0.$

Then using techniques from [1, Theorem 1.2, pp 489–490] further gives

$\displaystyle\left|\bm{1}_{\\{\xi^{(1)}_{w,\delta}(T_{n}-T_{w-1},L_{w})=k\\}}-\bm{1}_{\\{\xi^{(1)}_{w,\delta}\bigl{(}\log(n/w)/c^{*},m\bigr{)}=k\\}}\right|$
$\displaystyle\leq\sup_{w\geq 1}\sup_{1\leq m\leq
K}\sup_{t\in[0,-\log\eta/c^{*}]}\mathbb{P}\left(\xi^{(1)}_{w,\delta}(t+\epsilon,m)-\xi^{(1)}_{w,\delta}\bigl{(}(t-\epsilon)\wedge
0,m)\bigr{)}\geq 1\right)$ $\displaystyle\quad+\mathbb{P}\left(\sup_{n\eta\leq
w\leq
n}\left|T_{n}-T_{w-1}-\frac{1}{c^{*}}\log(n/w)\right|\geq\epsilon\right)=:p_{1}(\epsilon,\eta).$

Similarly,

$\displaystyle\left|\bm{1}_{\\{\xi^{(2)}_{w,\delta}(T_{n}-T_{w-1},L_{w})=l\\}}-\bm{1}_{\\{\xi^{(2)}_{w,\delta}\bigl{(}\log(n/w)/c^{*},L_{w}\bigr{)}=l\\}}\right|$
$\displaystyle\leq\sup_{w\geq 1}\sup_{1\leq m\leq
K}\sup_{t\in[0,-\log\eta/c^{*}]}\mathbb{P}\left(\xi^{(2)}_{w,\delta}(t+\epsilon,m)-\xi^{(2)}_{w,\delta}\bigl{(}(t-\epsilon)\wedge
0),m\bigr{)}\geq 1\right)$ $\displaystyle\quad+\mathbb{P}\left(\sup_{n\eta\leq
w\leq
n}\left|T_{n}-T_{w-1}-\frac{1}{c^{*}}\log(n/w)\right|\geq\epsilon\right)=:p_{2}(\epsilon,\eta).$

Then by (37), we see that

$\displaystyle\mathbb{E}|A_{1}(n)|$ $\displaystyle\leq 2\cdot\frac{1}{n}\cdot
n\eta+\frac{1}{n}(1-\eta)n\bigl{(}p_{1}(\epsilon,\eta)+p_{2}(\epsilon,\eta)\bigr{)},$

which implies $\lim_{n\to\infty}\mathbb{E}|A_{1}(n)|=0$. Therefore,
$A_{1}(n)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$.

For $A_{2}(n)$, define

$X_{w}:=\bm{1}_{\left\\{\bm{\xi}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),L_{w}\right)=(k,l)\right\\}}-\mathbb{P}^{\mathcal{F}_{T_{w-1}}}\left(\bm{\xi}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),L_{w}\right)=(k,l)\right),$

then we see that
$\mathbb{E}(X_{w})=\mathbb{E}\left(\mathbb{E}^{\mathcal{F}_{T_{w-1}}}(X_{w})\right)=0$.
Also, for $u<w$, since

$\mathbb{E}(X_{w}X_{u})=\mathbb{E}\left(\mathbb{E}^{\mathcal{F}_{T_{w-1}}}(X_{w}X_{u})\right)=\mathbb{E}\left(\mathbb{E}^{\mathcal{F}_{T_{w-1}}}(X_{w})\mathbb{E}^{\mathcal{F}_{T_{w-1}}}(X_{u})\right)=0,$

then by the weak law of large numbers, we have
$A_{2}(n)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$.

For $A_{3}(n)$, we first note that for $w\geq 2$ and $t\geq 0$,

$\mathbb{P}^{\mathcal{F}_{T_{w-1}}}\left(\bm{\xi}_{w,\delta}(t,L_{w})=(k,l)\right)=\mathbb{P}\left(\bm{\xi}_{w,\delta}(t,L_{w})=(k,l)\middle|\bm{\xi}_{w,\delta}(0,L_{w}),L_{w}\right).$

Then we have

$\displaystyle|A_{3}(n)|$
$\displaystyle\leq\frac{1}{n}\sum_{w=2}^{n+1}\bm{1}_{\left\\{\left(\bm{\xi}_{w,\delta}(0,L_{w}),L_{w}\right)\neq\left(\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right)\right\\}}.$

Hence, to prove $A_{3}(n)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$, it
suffices to show

$\displaystyle\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{P}\left[\left(\bm{\xi}_{w,\delta}(0,L_{w}),L_{w}\right)\neq\left(\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right)\right]\to
0.$ (38)

Let $\mathcal{X}:=\\{(0,1),(1,0),(1,1)\\}\times\\{1,\ldots,K\\}$, and we note
that

$\displaystyle\mathbb{P}$
$\displaystyle\left[\left(\bm{\xi}_{w,\delta}(0,L_{w}),L_{w}\right)\neq\left(\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right)\right]$
$\displaystyle\leq\sum_{\bm{x}\in\mathcal{X}}\left|\mathbb{P}\left[\left(\bm{\xi}_{w,\delta}(0,L_{w}),L_{w}\right)=\bm{x}\right]-\mathbb{P}\left[\left(\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),L_{w}\right)=\bm{x}\right]\right|$
$\displaystyle\leq
2\alpha\sum_{m=1}^{K}\rho_{m\bullet}\left|\mathbb{E}\left(\frac{|E^{\text{in}}_{m}(w-1)|+\delta|V_{m}(w-1)|}{|E(w-1)|+\delta|V(w-1)|}\right)-\frac{x_{m}+\delta\pi_{m}}{\sum_{s}x_{s}+\delta}\right|$
$\displaystyle\quad+2\gamma\sum_{m=1}^{K}\rho_{\bullet
m}\left|\mathbb{E}\left(\frac{|E^{\text{out}}_{m}(w-1)|+\delta|V_{m}(w-1)|}{|E(w-1)|+\delta|V(w-1)|}\right)-\frac{y_{m}+\delta\pi_{m}}{\sum_{s}y_{s}+\delta}\right|.$

From the proof of Theorem 2.2, we see that there exist constants
$C_{1}(\delta),C_{2}(\delta)>0$ such that

$\left|\mathbb{E}\left(\frac{|E^{\text{in}}_{m}(w-1)|+\delta|V_{m}(w-1)|}{|E(w-1)|+\delta|V(w-1)|}\right)-\frac{x_{m}+\delta\pi_{m}}{\sum_{s}x_{s}+\delta}\right|\leq
C_{1}(\delta)w^{\lambda_{H}-1},$

and

$\left|\mathbb{E}\left(\frac{|E^{\text{out}}_{m}(w-1)|+\delta|V_{m}(w-1)|}{|E(w-1)|+\delta|V(w-1)|}\right)-\frac{y_{m}+\delta\pi_{m}}{\sum_{s}y_{s}+\delta}\right|\leq
C_{2}(\delta)w^{\lambda_{H}-1}.$

Then the left hand side of (38) is bounded by

$\frac{2}{n}\sum_{w=2}^{n+1}\sum_{m=1}^{K}\left(\rho_{m\bullet}C_{1}(\delta)+\rho_{\bullet
m}C_{2}(\delta)\right)w^{\lambda_{H}-1},\qquad\lambda_{H}<1,$

which goes to 0 as $n\to\infty$, thus proving the claim in (38).

For $A_{4}(n)$, since
$\\{(\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}):w\geq
1\\}$ are iid random vectors in $\mathcal{X}$, then

$\displaystyle\frac{1}{n}\sum_{w=2}^{n+1}$
$\displaystyle\left(\mathbb{P}\left[\widetilde{\bm{\xi}}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),\widetilde{L}_{w}\right)=(k,l)\middle|\widetilde{\bm{\xi}}_{w,\delta}(0,\widetilde{L}_{w}),\widetilde{L}_{w}\right]\right.$
$\displaystyle\left.\qquad-\mathbb{P}\left[\widetilde{\bm{\xi}}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),\widetilde{L}_{w}\right)=(k,l)\right]\right)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0.$

Also, by the definition of $\widetilde{\bm{\xi}}_{w,\delta}$, $w\geq 1$, we
see that

$\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{P}\left[\widetilde{\bm{\xi}}_{w,\delta}\left(\frac{1}{c^{*}}\log(n/w),\widetilde{L}_{w}\right)=(k,l)\right]=\frac{1}{n}\sum_{w=2}^{n+1}\mathbb{P}\left[\widetilde{\bm{\xi}}_{\delta}\left(\frac{1}{c^{*}}\log(n/w),L^{*}\right)=(k,l)\right].$

Since the function $\mathbb{P}[\widetilde{\bm{\xi}}_{\delta}(t,L^{*})=(k,l)]$
is bounded and continuous in $t$, then we conclude that
$A_{4}(n)\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$ by applying the
Riemann integrability of $\mathbb{P}\left[\widetilde{\bm{\xi}}_{\delta}(-\log
t/c^{*},L^{*})=(k,l)\right]$, which completes the proof of (34). ∎

##  4 Power Laws and Asymptotic Dependence of Degree Frequencies

In this section, we study the dependence between large in- and out-degrees by
examining the asymptotic behavior of the distribution
$\mathbb{P}[(\mathcal{I}_{m},\mathcal{O}_{m})\in\cdot]$, for each
$m=1,\ldots,K$.

###  4.1 Multivariate and Hidden Regular Variation

To formalize our analysis, we provide some useful definitions related to
_multivariate regular variation_ (MRV) and _hidden regular variation_ (HRV) of
distributions.

Suppose that $\mathbb{C}_{0}\subset\mathbb{C}\subset\mathbb{R}_{+}^{2}$ are
two closed cones, and we provide the definition of $\mathbb{M}$-convergence in
Definition 4.1 (cf. [16, 10, 7, 15, 4]) on
$\mathbb{C}\setminus\mathbb{C}_{0}$, which lays the theoretical foundation of
regularly varying measures (cf. Definition 4.2).

######  Definition 4.1.

Let $\mathbb{M}(\mathbb{C}\setminus\mathbb{C}_{0})$ be the set of Borel
measures on $\mathbb{C}\setminus\mathbb{C}_{0}$ which are finite on sets
bounded away from $\mathbb{C}_{0}$, and
$\mathcal{C}(\mathbb{C}\setminus\mathbb{C}_{0})$ be the set of continuous,
bounded, non-negative functions on $\mathbb{C}\setminus\mathbb{C}_{0}$ whose
supports are bounded away from $\mathbb{C}_{0}$. Then for
$\mu_{n},\mu\in\mathbb{M}(\mathbb{C}\setminus\mathbb{C}_{0})$, we say
$\mu_{n}\to\mu$ in $\mathbb{M}(\mathbb{C}\setminus\mathbb{C}_{0})$, if $\int
f\differential\mu_{n}\to\int f\differential\mu$ for all
$f\in\mathcal{C}(\mathbb{C}\setminus\mathbb{C}_{0})$.

Without loss of generality [16], we can and do take functions in
$\mathcal{C}(\mathbb{C}\setminus\mathbb{C}_{0})$ to be uniformly continuous as
well. Denote the modulus of continuity of a uniformly continuous function
$f:\mathbb{R}_{+}^{p}\mapsto\mathbb{R}_{+}$ by

$\Delta_{f}(\delta)=\sup\\{|f(\bm{x})-f(\bm{y})|:d(\bm{x},\bm{y})<\delta\\}$
(39)

where $d(\cdot,\cdot)$ is an appropriate metric on the domain of $f$. Uniform
continuity means $\lim_{\delta\to 0}\Delta_{f}(\delta)=0.$ We now present the
definition of multivariate regular variation with
$\mathbb{C}=\mathbb{R}_{+}^{2}$ and $\mathbb{C}_{0}=\\{\bm{0}\\}$.

Following Definition 2.1 in [20], we denote a regularly varying function
$f:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$ with index $a\in\mathbb{R}$, as $f\in
RV_{a}$. Definition 4.2 gives the formal description of the MRV of
distributions.

######  Definition 4.2.

The distribution $\mathbb{P}(\bm{Z}\in\cdot)$ of a random vector $\bm{Z}$ on
$\mathbb{R}_{+}^{2}$, is (standard) regularly varying on
$\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\}$ with index $c>0$ if there exists
some regularly varying scaling function $b(t)\in\text{RV}_{1/c}$ and a limit
measure $\nu(\cdot)\in\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\})$
such that as $t\to\infty$,

$t\mathbb{P}(\bm{Z}/b(t)\in\cdot)\rightarrow\nu(\cdot),\qquad\text{in
}\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\}).$ (40)

It is convenient to write
$\mathbb{P}(\bm{Z}\in\cdot)\in\text{MRV}(c,b(t),\nu,\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\})$.

When analyzing the asymptotic dependence between components of a bivariate
random vector $\bm{Z}$ satisfying (40), it is often informative to make a
polar coordinate transform and consider the transformed points located on the
$L_{1}$ unit sphere

$\displaystyle(x,y)\mapsto\left(\frac{x}{x+y},\frac{y}{x+y}\right),$ (41)

after thresholding the data according to the $L_{1}$ norm.

When a limit measure concentrates on a subcone of the full state space, to
improve estimates of probabilities in the complement of the subcone, we can
seek a second hidden regular variation regime after removing the subcone.

######  Definition 4.3.

The vector $\bm{Z}$ is regularly varying on
$\mathbb{R}^{2}_{+}\setminus\\{\bm{0}\\}$ and has hidden regular variation on
$\mathbb{R}^{2}_{+}\setminus\mathbb{C}_{0}$ if there exist $0<c\leq c_{0}$,
scaling functions $b(t)\in RV_{1/c}$ and $b_{0}(t)\in RV_{1/c_{0}}$ with
$b(t)/b_{0}(t)\to\infty$ and limit measures $\nu$, $\nu_{0}$, such that

$\displaystyle\mathbb{P}(\bm{Z}\in\cdot)\in\text{MRV}(c,b(t),\nu,\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\})\cap\text{MRV}(c_{0},b_{0}(t),\nu_{0},\mathbb{R}_{+}^{2}\setminus\\{\mathbb{C}_{0}\\}).$
(42)

A convenient way to characterize HRV is through the _generalized polar
coordinate transformation_ for $\mathbb{R}_{+}^{2}\setminus\mathbb{C}_{0}$ and
an associated metric $d(\cdot,\cdot)$ satisfying $d(cx,cy)=cd(x,y)$ for
scalars $c>0$. The metric $d(\cdot,\cdot)$ that we use in practice is the
$L_{1}$-metric. When using generalized polar coordinates with respect to the
forbidden zone $\mathbb{C}_{0}$, we define
$\aleph_{\mathbb{C}_{0}}:=\\{\bm{x}\in\mathbb{C}\setminus\mathbb{C}_{0}:d(\bm{x},\mathbb{C}_{0})=1\\}$,
the locus of points at distance 1 from $\mathbb{C}_{0}$. Then the generalized
polar coordinates are specified through the transformation,
$\text{GPOLAR}:\mathbb{R}_{+}^{2}\setminus\mathbb{C}_{0}\mapsto(0,\infty)\times\aleph_{\mathbb{C}_{0}}$
with

$\text{GPOLAR}(\bm{x})=\left(d(\bm{x},\mathbb{C}_{0}),\frac{\bm{x}}{d(\bm{x},\mathbb{C}_{0})}\right).$

Let $\nu_{c}(\cdot)$ be a measure in
$\mathbb{M}(\mathbb{R}_{+}\setminus\\{0\\})$ satisfying
$\nu_{c}(x,\infty)=x^{-c}$ $x,c>0$, and $S_{0}(\cdot)$ be a probability
measure on $\aleph_{\mathbb{C}_{0}}$. Then generalized polar coordinates allow
re-writing (42) as

$t\mathbb{P}\left[\left(\frac{d(\bm{Z},\mathbb{C}_{0})}{b_{0}(t)},\frac{\bm{Z}}{d(\bm{Z},\mathbb{C}_{0})}\right)\in\cdot\right]\to(\nu_{c_{0}}\times
S_{0})(\cdot)$

in
$\mathbb{M}((\mathbb{R}_{+}\setminus\\{0\\})\times\aleph_{\mathbb{C}_{0}})$.
See [7] and [16] for details.

###  4.2 Degree Frequencies and HRV

Let $L^{*}$ be a random variable with pmf $\mathbb{P}(L^{*}=m)=\pi_{m}$,
$m=1,\ldots,K$, independent from $T^{*}$ and
$\\{\widetilde{\bm{\xi}}_{\delta}(\cdot,m):m=1,\ldots,K\\}$. Using (34), the
limit in (35) becomes

$\sum_{m=1}^{K}\pi_{m}\mathbb{P}\Bigl{[}\big{(}\mathcal{I}_{m},\mathcal{O}_{m}\bigr{)}=(k,l)\Bigr{]}=\mathbb{P}\left(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*})=(k,l)\right)=:{\mathbb{P}\Big{(}\bigl{(}\mathcal{I},\mathcal{O}\bigr{)}=(k,l)\Bigr{)}},$

the limiting empirical proportion of nodes with in-degree $k$ and out-degree
$l$.

We now discuss the regular variation properties of this limit distribution.
The branching structure of $\widetilde{\bm{\xi}}_{\delta}(\cdot,m)$ is
specified through the matrix $A_{m}$ given in (19), and the largest eigenvalue
of $A_{m}$ is $\lambda_{m}$ given in (20). We assume parameters $\alpha$,
$\gamma$ and the matrix $\bm{\rho}$, are chosen such that $\lambda_{m}$,
$m=1,\ldots,K$, are all distinct values and without loss of generality, assume
the behavioral group labels are chosen so that

$\lambda_{1}>\lambda_{2}>\cdots>\lambda_{K},$ (43)

Theorem 4.1 gives the multivariate regular variation of
$\mathbb{P}[(\mathcal{I},\mathcal{O})\in\cdot\,]$ on
$\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\}$. The proof is based on an extended
Breiman’s theorem [26, Theorem 3], reviewed in Theorem A.1. Theorem 4.1
requires an additional assumption that $\lambda_{1}\geq\log 2$ to guarantee
that the moment condition in (69) is satisfied.

######  Theorem 4.1.

Recall the definition of $\rho^{*},\,c^{*}$ in (33), suppose
$\lambda_{1}\geq\log 2$ and that the regularity conditions (15) hold. Then as
$t\to\infty$,

$t\mathbb{P}\left[\frac{(\mathcal{I},\mathcal{O})}{t^{\lambda_{1}/c^{*}}}\in\cdot\right]\longrightarrow\mu_{1},\qquad\text{in}\quad\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\}),$

where the limiting measure
$\mu_{1}\in\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\})$ satisfies that
for $f\in\mathcal{C}(\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\})$,

$\displaystyle\mu_{1}(f)=$
$\displaystyle\mathbb{P}(L^{*}=1)\int_{0}^{\infty}\mathbb{E}\left(f(y\widetilde{Z}(1)\bm{v}(1))\right)\nu_{c^{*}/\lambda_{1}}(\mathrm{d}y).$
$\displaystyle=$
$\displaystyle\mathbb{P}(L^{*}=1)\mathbb{E}\left(\widetilde{Z}(1)^{c^{*}/\lambda_{1}}\right)\int_{0}^{\infty}f\bigl{(}s\bm{v}(1)\bigr{)}\nu_{c^{*}/\lambda_{1}}(\mathrm{d}s).$

Here $\widetilde{Z}(1)$ satisfies
$e^{-\lambda_{1}t}\widetilde{\bm{\xi}}_{\delta}(t,1)\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\widetilde{Z}(1)\bm{v}(1)$,
and

$a(1):=\frac{v^{(2)}{(1)}}{v^{(1)}{(1)}}=\frac{\gamma-\alpha+\sqrt{D_{0}(1)}}{2\gamma\rho_{\bullet
1}}.$ (44)

Write

$a(m):={v^{(2)}{(m)}}/{v^{(1)}{(m)}},\quad m=1,\ldots,K.$ (45)

Theorem 4.1 shows that the limiting measure $\mu_{1}$ concentrates Pareto mass
on the ray $\mathcal{L}_{(1)}:=\\{(x,y)\in\mathbb{R}_{+}^{2}:y=a{(1)}x\\}$ and
concentration on $\mathcal{L}_{(1)}$ suggests that at scale
$b(t)=t^{\lambda_{1}/c^{*}}$, large in- and out-degree pairs in the PA model
with heterogeneous reciprocity levels satisfy $\mathcal{O}\approx
a{(1)}\mathcal{I}$.

###### Proof.

Similar to the proof of Theorem 6 in [26], we first see that
$\mathbb{P}(\widetilde{Z}_{1}>0)=1$ by the property of MBI processes. Also,
since $\lambda_{1}>\lambda_{r}$ for $r=2,\ldots,K$,

$e^{-\lambda_{1}t}\widetilde{\bm{\xi}}_{\delta}(t,{r})\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}0.$

Therefore,

$e^{-\lambda_{1}t}\widetilde{\bm{\xi}}_{\delta}(t,L^{*})\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\widetilde{Z}(1)\bm{v}(1)\bm{1}_{\\{L^{*}=1\\}}$

Then the proof of Theorem 4.1 is an application of Theorem A.1 after making
the identifications

$\displaystyle\bm{\xi}(t)=t^{-1}\widetilde{\bm{\xi}}_{\delta}\left(\frac{1}{\lambda_{1}}\log
t,L^{*}\right),$
$\displaystyle\bm{\xi}_{\infty}=\widetilde{Z}(1)\bm{v}(1)\bm{1}_{\\{L^{*}=1\\}},$
$\displaystyle X=e^{\lambda_{1}T^{*}},$ $\displaystyle
b(t)=t^{\lambda_{1}/c^{*}},$ $\displaystyle c=c^{*}/\lambda_{1}.$

The remaining piece is to show the moment condition (69) in this context and
we will show any $\delta\geq 0$ and any $q=1,2,\ldots$, there exists some
constant $C(\delta,q)>0$ such that

$\displaystyle\sup_{t\geq
0}e^{-\lambda_{1}qt}\mathbb{E}\left[\left(\widetilde{\xi}^{(1)}_{\delta}(t,1)\right)^{q}\right]\leq
C(\delta,q),$ (46)

which is true by Proposition 2 in [26]. ∎

Next, Theorem 4.2 gives a second hidden regular variation (HRV) regime after
removing $\mathcal{L}_{(1)}$ [8, 9, 20, 7, 16]. The existence of HRV has been
detected empirically in network data [9], and here we theoretically prove HRV
present in the PA model with heterogeneous reciprocity.

The limit measure given in Theorem 4.1 concentrates on $\mathcal{L}_{(1)}$.
Thus, we may seek a regular variation property on
$\mathbb{R}_{+}^{2}\setminus\mathcal{L}_{(1)}$ using a weaker scaling function
$b_{0}(t)$. A convenient way to seek the hidden regular variation is by using
generalized polar coordinates which in this case amount to the transformation

$\bm{x}\to\Bigl{(}d_{1}(\bm{x},\mathcal{L}_{(1)}),\frac{\bm{x}}{d_{1}(\bm{x},\mathcal{L}_{(1)})}\Bigr{)},$

where $d_{1}(\bm{x},\bm{y})$ is a metric on
$\mathbb{R}_{+}^{2}\setminus\\{\bm{0}\\}$ chosen for convenience to be the
$L_{1}$-metric. The $L_{1}$-distance of a point $(x,y)$ to the line
$\mathcal{L}_{(1)}$ is readily computed to be

$d_{1}\bigl{(}(x,y),\mathcal{L}_{(1)}\bigr{)}={|y-a(1)x|}/{\text{max}\\{1,a(1\\}},$

and we use a scaled version

$d^{\prime}_{1}\bigl{(}(x,y),\mathcal{L}_{(1)}\bigr{)}=|y-a(1)x|.$ (47)

Define
$\aleph_{\mathcal{L}_{(1)}}:=\\{\bm{x}\in\mathbb{R}_{+}^{2}\setminus\mathcal{L}_{(1)}:d^{\prime}_{1}(\bm{x},\mathcal{L}_{(1)})=1\\}$,
which are 2 line segments in $\mathbb{R}_{+}^{2}$ parallel to
$\mathcal{L}_{(1)}$. Hidden regular variation will be present for
$(\mathcal{I},\mathcal{O})\stackrel{{\scriptstyle
d}}{{=}}\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*})$ if

$tP\left[\left(\frac{d^{\prime}_{1}(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*}),\mathcal{L}_{(1)})}{b_{0}(t)},\frac{\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*})}{d^{\prime}_{1}\left(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*}),\mathcal{L}_{(1)}\right)}\right)\in\cdot\,\right]$

converges to a limit measure in
$\mathbb{M}((\mathbb{R}_{+}\setminus\\{0\\})\times\aleph_{\mathcal{L}_{(1)}})$.
The next theorem explains the convergence on
$\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\mathcal{L}_{(1)})$.

######  Theorem 4.2.

Assume that $\lambda_{2}>\lambda_{1}/2$, and $\lambda_{2}\geq\log 2$. Let
$\widetilde{Z}(2)$ be the limiting random variable satisfying
$e^{-t\lambda_{2}}\widetilde{\bm{\xi}}_{\delta}(t,2)\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\widetilde{Z}(2)\bm{v}(2)$
as $t\to\infty$. Then we have in
$\mathbb{M}(\mathbb{R}^{2}_{+}\setminus\mathcal{L}_{(1)})$ that

$\displaystyle
t\mathbb{P}\left[\frac{(\mathcal{I},\mathcal{O})}{t^{\lambda_{2}/c^{*}}}\in\cdot\right]\to\mu_{2},$
(48)

where the limit measure
$\mu_{2}\in\mathbb{M}(\mathbb{R}_{+}^{2}\setminus\mathcal{L}_{(1)})$
concentrates on the ray $y=a(2)x,$ $x>0$ (recall (45)) in the first quadrant
and satisfies for
$g\in\mathcal{C}(\mathbb{R}_{+}^{2}\setminus\mathcal{L}_{(1)})$,

$\displaystyle\mu_{2}(g)=$
$\displaystyle\mathbb{P}(L^{*}=2)\int_{0}^{\infty}\mathbb{E}\left(g(y\widetilde{Z}(2)\bm{v}(2))\right)\nu_{c^{*}/\lambda_{2}}(\mathrm{d}y).$
$\displaystyle=$
$\displaystyle\mathbb{P}(L^{*}=2)\mathbb{E}\left(\widetilde{Z}(2)^{c^{*}/\lambda_{2}}\right)\int_{0}^{\infty}g(y\bm{v}(2))\nu_{c^{*}/\lambda_{2}}(\mathrm{d}y).$

Theorem 4.2 suggests that after deleting large in- and out-degree pairs close
to $\mathcal{L}_{(1)}$, at the scale $t^{\lambda_{2}/c^{*}}$ the remaining
large observations of in- and out-degrees tend to concentrate around another
line

$\mathcal{L}_{(2)}:=\\{(x,y)\in(0,\infty)^{2}:y=a(2)x\\}.$

Define also that $\mathcal{L}_{(0)}:=\\{\bm{0}\\}$, and combining Theorems 4.1
and 4.2 gives

$\mathbb{P}((\mathcal{I},\mathcal{O})\in\cdot)\in\bigcap_{m=1}^{2}\text{MRV}\left(c^{*}/\lambda_{m},t^{\lambda_{m}/c^{*}},\mu_{m},\mathbb{R}_{+}^{2}\setminus\left(\bigcup_{i=0}^{m-1}\mathcal{L}_{(i)}\right)\right).$

######  Remark 4.1.

Results in Theorem 4.2 can be extended as follows. Let $\lambda_{m_{0}}$,
$m_{0}\geq 2$, denote the $m_{0}$-th largest eigenvalue such that for all
$m\in\\{2,\ldots,m_{0}\\}$,

$\lambda_{m}>\lambda_{m-1}/2,\qquad\text{and}\qquad\lambda_{m}\geq\log 2.$

For $m=1,\ldots,K$, write
$\mathcal{L}_{(m)}:=\\{(x,y)\in(0,\infty)^{2}:y=a(m)x\\}$, and set
$\widetilde{Z}(m)$ to be the limiting random variable satisfying
$e^{-t\lambda_{m}}\widetilde{\bm{\xi}}_{\delta}(t,m)\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\widetilde{Z}(m)\bm{v}(m)$
as $t\to\infty$. Define also the measure
$\mu_{m}\in\mathbb{M}\left(\mathbb{R}_{+}^{2}\setminus\bigl{(}\bigcup_{i=1}^{m}\mathcal{L}_{(i)}\bigr{)}\right)$
such that for
$g\in\mathcal{C}\left(\mathbb{R}_{+}^{2}\setminus\bigl{(}\bigcup_{i=1}^{m}\mathcal{L}_{(i)}\bigr{)}\right)$,

$\mu_{m}(g)=\mathbb{P}(L=m)\int_{0}^{\infty}\mathbb{E}\left(g(y\widetilde{Z}(m)\bm{v}(m\right)\nu_{c^{*}/\lambda_{m}}(\mathrm{d}y).$

Then applying the proof technique of Theorem 4.2 for $m_{0}$ times gives

$\mathbb{P}\bigl{(}(\mathcal{I},\mathcal{O})\in\cdot\bigr{)}\in\bigcap_{m=1}^{m_{0}}\text{MRV}\left(c^{*}/\lambda_{m},t^{\lambda_{m}/c^{*}},\mu_{m},\mathbb{R}_{+}^{2}\setminus\left(\bigcup_{i=0}^{m-1}\mathcal{L}_{(i)}\right)\right).$

###### Proof.

Define
$\bm{\theta}_{(2)}:=(1,a(2))/\left|a(1)-a(2)\right|)\in\aleph_{\mathcal{L}_{1}}$.
Applying the generalized polar transformation shows that verifying (48) is
equivalent to justifying

$\displaystyle t\mathbb{P}\Bigl{[}\Bigl{(}$
$\displaystyle\frac{d_{1}^{\prime}(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*}),\mathcal{L}_{1})}{t^{\lambda_{2}/c^{*}}},\frac{\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*})}{d_{1}^{\prime}(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*}),\mathcal{L}_{1})}\Bigr{)}\in\cdot\,\Bigr{]}$
$\displaystyle=$ $\displaystyle
t\mathbb{P}\left[\left(\frac{\left|a{(1)}\widetilde{\xi}^{(1)}_{\delta}(T^{*},L^{*})-\widetilde{\xi}^{(2)}_{\delta}(T^{*},L^{*})\right|}{t^{\lambda_{2}/c^{*}}},\frac{\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*})}{\left|a{(1)}\widetilde{\xi}^{(1)}_{\delta}(T^{*},L^{*})-\widetilde{\xi}^{(2)}_{\delta}(T^{*},L^{*})\right|}\right)\in\cdot\right]$
$\displaystyle\to$ $\displaystyle
C_{2}\nu_{c^{*}/\lambda_{2}}(\cdot)\times\epsilon_{\bm{\theta}_{(2)}}(\cdot),$
(49)

in $\mathbb{M}\left((0,\infty)\times\aleph_{[\mathcal{L}_{(1)}]}\right)$,
where

$\displaystyle C_{2}=\mathbb{P}(L^{*}=2)\frac{c^{*}}{\lambda_{2}}$
$\displaystyle\times\int_{0}^{\infty}z^{-1-c^{*}/\lambda_{2}}\mathbb{P}\left(\widetilde{Z}(2)>\frac{1/z}{v^{(1)}{(2)}\left|a(1)-a(2)\right|}\right)\mathrm{d}z.$

To prove (4.2), we first claim that for $\lambda_{2}>\lambda_{1}/2$, as
$t\to\infty$,

$\displaystyle
e^{-t\lambda_{2}}d_{1}^{\prime}(\widetilde{\bm{\xi}}_{\delta}(T^{*},L^{*}),\mathcal{L}_{1})=$
$\displaystyle
e^{-t\lambda_{2}}\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,L^{*})-\widetilde{\xi}^{(2)}_{\delta}(t,L^{*})\right|$
$\displaystyle\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}$
$\displaystyle\left|a(1)-a(2)\right|v^{(1)}{(2)}\widetilde{Z}(2)\bm{1}_{\\{L^{*}=2\\}}.$
(50)

In addition, since

$\displaystyle\frac{\widetilde{\bm{\xi}}_{\delta}(t,m)}{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,m)-\widetilde{\xi}^{(2)}_{\delta}(t,m)\right|}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\frac{(1,a(m))}{\bigl{|}a(1)-a(m)\bigr{|}},$

we have

$\displaystyle\left(\frac{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,L^{*})-\widetilde{\xi}^{(2)}_{\delta}(t,L^{*})\right|}{e^{t\lambda_{2}}},\,\frac{\widetilde{\bm{\xi}}_{\delta}(t,L^{*})}{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,L^{*})-\widetilde{\xi}^{(2)}_{\delta}(t,L^{*})\right|}\right)$
$\displaystyle=\sum_{m=1}^{K}\left(\frac{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,m)-\widetilde{\xi}^{(2)}_{\delta}(t,m)\right|}{e^{t\lambda_{2}}},\,\frac{\widetilde{\bm{\xi}}_{\delta}(t,m)}{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,m)-\widetilde{\xi}^{(2)}_{\delta}(t,m)\right|}\right)\bm{1}_{\\{L^{*}=m\\}}$
$\displaystyle\quad\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\left(\bigl{|}a(1)-a(2)\bigr{|}v^{(1)}{(2)}\widetilde{Z}(2)\bm{1}_{\\{L^{*}=2\\}},\,\sum_{m=1}^{K}\frac{(1,a(m))}{\bigl{|}a(1)-a(m)\bigr{|}}\bm{1}_{\\{L^{*}=m\\}}\right).$
(51)

Equation (51) further gives that in
$\mathbb{M}\left((0,\infty)\times\left((0,\infty)\times\aleph_{[\mathcal{L}_{(1)}]}\right)\right)$,

$\displaystyle
t\mathbb{P}\left[\left(\frac{e^{T^{*}}}{t^{1/c^{*}}},\right.\right.$
$\displaystyle\left.\left.\sum_{m=1}^{K}\left(\frac{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(T^{*},m)-\widetilde{\xi}^{(2)}_{\delta}(T^{*},m)\right|}{e^{\lambda_{2}T^{*}}},\frac{\widetilde{\bm{\xi}}_{\delta}(T^{*},m)}{\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(T^{*},m)-\widetilde{\xi}^{(2)}_{\delta}(T^{*},m)\right|}\right)\bm{1}_{\\{L^{*}=m\\}}\right)\in\cdot\right]$
$\displaystyle\longrightarrow\mathbb{P}(L^{*}=2)\,\nu_{c^{*}/\lambda_{2}}(\cdot)\times\mathbb{P}\left(\bigl{|}a(1)-a(2)\bigr{|}v^{(1)}{(2)}\widetilde{Z}(2)\in\cdot\right)\times\epsilon_{\bm{\theta}_{(2)}}(\cdot).$

Hence, as long as we check the moment condition that for $q=1,2,\ldots$ and
$\delta>0$,

$\displaystyle\sup_{t\geq
0}e^{-qt\lambda_{2}}\mathbb{E}\left[\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,2)-\widetilde{\xi}^{(2)}_{\delta}(t,2)\right|^{q}\right]<\infty,$
(52)

then applying the generalized Breiman’s theorem (cf. Theorem A.1) completes
the proof of (4.2). To verify (52), we notice that since

$\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,2)-\widetilde{\xi}^{(2)}_{\delta}(t,2)\right|\leq
a(1)\widetilde{\xi}^{(1)}_{\delta}(t,2)+\widetilde{\xi}^{(2)}_{\delta}(t,2),$

it suffices to show

$\sup_{t\geq
0}e^{-qt\lambda_{2}}\mathbb{E}\left[\left(a(1)\widetilde{\xi}^{(1)}_{\delta}(t,2)+\widetilde{\xi}^{(2)}_{\delta}(t,2)\right)^{q}\right]<\infty,\qquad\lambda_{2}\geq\log
2,$

which is true by Proposition 2 in [26].

It remains to prove the claim in (4.2). Since
$\left|a(m)\widetilde{\xi}^{(1)}_{\delta}(t,m)-\widetilde{\xi}^{(2)}_{\delta}(t,m)\right|\leq
a(m)\widetilde{\xi}^{(1)}_{\delta}(t,m)+\widetilde{\xi}^{(2)}_{\delta}(t,m)$,
then for $m\notin\\{1,2\\}$, we have

$e^{-t\lambda_{2}}\left|a(m)\widetilde{\xi}^{(1)}_{\delta}(t,m)-\widetilde{\xi}^{(2)}_{\delta}(t,m)\right|\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}0,\qquad\text{as
}t\to\infty.$

Hence, we only need to consider $m=1,2$. Consider the case when $m=1$, and
define $\\{\bm{\psi}_{i}(t,1):t\geq 0\\}_{i\geq 1}$ as a sequence of iid two-
type branching processes (without immigration) with group label $1$, whose
branching structure is specified through $A_{1}$. Also, assume that
$\bm{\psi}_{i}(0,1)$ is a 2-dimensional random vector with distribution
$p_{0}(\bm{r},1)$ (cf. (18)). Let $0<\tau_{1}<\tau_{2}<\ldots$ be arrival
times of points in a homogeneous Poisson process with rate $\delta>0$, which
is independent from $\\{\bm{\psi}_{i}(t,1):t\geq 0\\}_{i\geq 1}$. Then by the
distributional construction of the MBI process in [19], we write

$\displaystyle\widetilde{\bm{\xi}}_{\delta}(t,1)=\sum_{i=1}^{\infty}\bm{\psi}_{i}(t-\tau_{i},1)\bm{1}_{\\{t\geq\tau_{i}\\}}.$
(53)

Next, we use a Borel Cantelli argument to show

$\displaystyle
e^{-t\lambda_{2}}\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(t,1)-\widetilde{\xi}^{(2)}_{\delta}(t,1)\right|\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}0,$
(54)

as $t\to\infty$, i.e. we will show that

$\displaystyle\sum_{n=1}^{\infty}\mathbb{P}\left(e^{-n\lambda_{2}}\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(n,1)-\widetilde{\xi}^{(2)}_{\delta}(n,1)\right|>\epsilon\right)<\infty.$
(55)

By Markov’s inequality, we have for $n\geq 1$,

$\displaystyle\mathbb{P}$
$\displaystyle\left(e^{-n\lambda_{2}}\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(n,1)-\widetilde{\xi}^{(2)}_{\delta}(n,1)\right|>\epsilon\right)$
$\displaystyle\leq\epsilon^{-2}\mathbb{E}\left[e^{-2n\lambda_{2}}\left|a(1)\widetilde{\xi}^{(1)}_{\delta}(n,1)-\widetilde{\xi}^{(2)}_{\delta}(n,1)\right|^{2}\right]$
$\displaystyle\leq\epsilon^{-2}e^{-n(2\lambda_{2}-\lambda_{1})}\sup_{t\geq
0}\mathbb{E}\left[\left(e^{-t\lambda_{1}/2}\left(a(1)\widetilde{\xi}^{(1)}_{\delta}(t,1)-\widetilde{\xi}^{(2)}_{\delta}(t,1\right)\right)^{2}\right].$
(56)

Note that the vector $[a(1),-1]^{T}$ is the right eigenvector associated with
the smaller eigenvalue of $A_{1}$. Recall the discussion in Section 3.1 that
$\lambda^{\prime}_{1}$ is the smaller eigenvalue of $A_{1}$, then the
condition $\lambda_{1}\geq\log 2$ guarantees that
$\lambda_{1}>2\lambda^{\prime}_{1}$. Hence, we apply Theorem V.7.1(i) in [2]
to conclude that for $i\geq 1$,

$\sup_{t\geq
0}\mathbb{E}\left[\left(e^{-t\lambda_{1}/2}\left(a(1)\psi^{(1)}_{i}(t,1)-\psi^{(2)}_{i}(t,1\right)\right)^{2}\right]<\infty.$

We then conclude from (53) that

$\displaystyle\left(e^{-t\lambda_{1}/2}\left(a(1)\widetilde{\xi}^{(1)}_{\delta}(t,1)-\widetilde{\xi}^{(2)}_{\delta}(t,1\right)\right)^{2}$
$\displaystyle=\left(e^{-t\lambda_{1}/2}\sum_{i=1}^{\infty}\left(a(1)\psi_{i}^{(1)}(t-\tau_{i},1)-\psi_{i}^{(2)}(t-\tau_{i},1\right)\bm{1}_{\\{t\geq\tau_{i}\\}}\right)^{2}$
$\displaystyle\leq\left(\sum_{i=1}^{\infty}e^{-\tau_{i}\lambda_{1}/2}\right)$
$\displaystyle\quad\times\left(\sum_{i=1}^{\infty}e^{-\tau_{i}\lambda_{1}/2}\left(e^{-(t-\tau_{i})\lambda_{1}/2}\left(a(1)\psi_{i}^{(1)}(t-\tau_{i},1)-\psi_{i}^{(2)}(t-\tau_{i},1)\right)\bm{1}_{\\{t\geq\tau_{i}\\}}\right)^{2}\right).$

Therefore,

$\displaystyle\mathbb{E}$
$\displaystyle\left(e^{-t\lambda_{1}/2}\left(a(1)\widetilde{\xi}^{(1)}_{\delta}(t,1)-\widetilde{\xi}^{(2)}_{\delta}(t,1\right)\right)^{2}$
$\displaystyle\leq\sup_{t\geq
0}\mathbb{E}\left[\left(e^{-t\lambda_{1}/2}\left(a(1)\psi^{(1)}_{i}(t,1)-\psi^{(2)}_{i}(t,1\right)\right)^{2}\right]$
$\displaystyle\quad\times\mathbb{E}\left[\sum_{i=1}^{\infty}e^{-\tau_{i}\lambda_{1}/2}\left(\sum_{i=1}^{\infty}e^{-\tau_{i}\lambda_{1}/2}\right)\right]<\infty.$

Also, since we assume $\lambda_{2}>\lambda_{1}/2$, then (55) follows from
(56), leading to (54). Then the claim in (4.2) follows by realizing that

$\displaystyle
e^{-t\lambda_{2}}\left(a(1)\widetilde{\xi}^{(1)}_{\delta}(t,2)-\widetilde{\xi}^{(2)}_{\delta}(t,2)\right)$
$\displaystyle\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\left(a(1)v^{(1)}(2)-v^{(2)}(2)\right)\widetilde{Z}(2)$
$\displaystyle=\left(a(1)-a(2)\right)v^{(1)}(2)\widetilde{Z}(2).$

∎

##  5 Concluding Remarks

In this paper, we propose a preferential attachment model with heterogeneous
reciprocity levels and study its theoretical properties. Using the MBI
embedding technique, we find that the distribution of large in- and out-
degrees is jointly regularly varying, concentrating along a specific ray.
Additionally, after deleting large in- and out-degree pairs close to the ray,
we further have hidden regular variation with limit measure concentrating
around another ray.

We now outline some open problems related to this model, which will be left as
future research.

#### Estimation.

The fitting of the proposed model remains open at this point. For a given
dataset, we need to first decide how many communication groups ($K$) should be
assumed. One may consult classical clustering methods such as $k$-means and
$k$-nearest neighbors to determine a proper $K$ beforehand, but their
applications to the network framework needs rigorous justification. Once $K$
is chosen, we can also derive proper tools (similar to those in [9]) to detect
the existence of hidden regular variation under the network setup.

#### Varying group labels.

So far we have assumed that the communication group of a user (node) is
determined upon its creation and remains unchanged afterwards. For real-world
applications, however, this may be a naive assumption, since users’
interaction patterns can change over time. A possible extension is to assume
that the group label for each node follows another Markov chain with finite
state space. One may consider applying variational Bayesian inference methods
(cf. [17]) to these extended models.

#### Covariate-dependent reciprocity probabilities.

The model studied in this paper assumes a deterministic matrix, $\bm{\rho}$,
to characterize the reciprocation between different communication groups. If
more node-specific information (covariates) is available, e.g. various
demographic, socio-economic, and behavioral factors, then we can generalize
the definition of the matrix $\bm{\rho}$ to be covariate-dependent. This may
possibly lead to heterogeneous extremal behaviors.

##  6 Proofs of Results in Section 2

###  6.1 Proof of Lemma 2.1

Define a function $\bm{f}:\mathcal{Z}\mapsto\mathcal{Z}$ such that for
$i=1,\ldots,K$,

$\displaystyle
f_{i}(\bm{z})=\alpha\frac{z_{i}+\delta\pi_{i}}{\sum_{j=1}^{K}z_{j}+\delta}+\gamma\rho_{\bullet
i}\frac{z_{K+i}+\delta\pi_{i}}{\sum_{j=1}^{K}z_{K+j}+\delta}+\gamma\pi_{i}+\alpha\pi_{i}\sum_{r=1}^{K}\rho_{r,i}\frac{z_{i}+\delta\pi_{i}}{\sum_{j=1}^{K}z_{j}+\delta},$

and for $i=K+1,\ldots,2K$,

$\displaystyle f_{i}(\bm{z})=$
$\displaystyle\gamma\frac{z_{i}+\delta\pi_{i-K}}{\sum_{j=K+1}^{2K}z_{j}+\delta}+\alpha\rho_{(i-K)\,\bullet}\frac{z_{i-K}+\delta\pi_{i-K}}{\sum_{j=1}^{K}z_{j}+\delta}+\alpha\pi_{i-K}$
$\displaystyle+\gamma\pi_{i-K}\sum_{r=1}^{K}\rho_{i,r}\frac{z_{i}+\delta\pi_{i-K}}{\sum_{j=K+1}^{2K}z_{j}+\delta}.$

Let $J(\bm{z}^{\prime})$ be the Jacobian matrix evaluated at some
$\bm{z}^{\prime}$ between $\bm{z}_{1}$ and $\bm{z}_{2}$. Use $\|\bm{z}\|_{1}$
to denote the $L_{1}$-norm for some $\bm{z}\in\mathcal{Z}$, and set
$\|\cdot\|_{(1,\infty)}$ to be the $L_{1,\infty}$ norm of a matrix. Then by
the mean value theorem, we see that for $\bm{z}_{1},\bm{z}_{2}\in\mathcal{Z}$,

$\displaystyle\|\bm{f}(\bm{z}_{1})-\bm{f}(\bm{z}_{2})\|_{1}$
$\displaystyle\leq\sup_{\bm{z}^{\prime}\in\mathcal{Z}}\|J(\bm{z}^{\prime})\|_{(1,\infty)}\|\bm{z}_{1}-\bm{z}_{2}\|_{\infty}$
$\displaystyle\leq\sup_{\bm{z}^{\prime}\in\mathcal{Z}}\|J(\bm{z}^{\prime})\|_{(1,\infty)}\|\bm{z}_{1}-\bm{z}_{2}\|_{1}.$
(57)

If we can show
$\sup_{\bm{z}^{\prime}\in\mathcal{Z}}\|J(\bm{z}^{\prime})\|_{(1,\infty)}<1$,
then by the contraction mapping theorem (cf. Theorem 1.2.2 in [13]) we are
able to show the existence of a unique solution to $\bm{f}(\bm{z})=\bm{z}$.

To find an upper bound for
$\sup_{\bm{z}\in\mathcal{Z}}\|J(\bm{z})\|_{(1,\infty)}$, we now give upper
bounds for the absolute value of each entry in $J(\bm{z})$. Let
$J_{i,j}(\bm{z})$ be the $(i,j)$-th entry in $J(\bm{z})$, and we have for
$1\leq i,j\leq K$,

$\displaystyle\left|J_{i,j}(\bm{z})\right|$
$\displaystyle=\left|\frac{\partial f_{i}}{\partial
z_{j}}(\bm{z})\right|\leq\alpha(\pi_{i}\rho_{j,i}+\bm{1}_{\\{i=j\\}})\frac{1}{\sum_{m=1}^{K}z_{m}+\delta}\leq\frac{\alpha(\pi_{i}\rho_{j,i}+\bm{1}_{\\{i=j\\}})}{1+\delta}.$
(58)

For $1\leq i\leq K$ and $K+1\leq j\leq 2K$,

$\displaystyle\left|J_{i,j}(\bm{z})\right|$
$\displaystyle=\left|\frac{\partial f_{i}}{\partial
z_{j}}(\bm{z})\right|=\gamma\rho_{\bullet
i}\left|-\frac{z_{K+i}+\delta\pi_{i}}{\left(\sum_{m=K+1}^{2K}z_{m}+\delta\right)^{2}}+\frac{\bm{1}_{\\{j=K+i\\}}}{\sum_{m=K+1}^{2K}z_{m}+\delta}\right|\leq\frac{\gamma\rho_{\bullet
i}}{1+\delta}.$ (59)

For $K+1\leq i\leq 2K$ and $1\leq j\leq K$, we see that

$\displaystyle\left|J_{i,j}(\bm{z})\right|$
$\displaystyle=\left|\frac{\partial f_{i}}{\partial
z_{j}}(\bm{z})\right|=\alpha\rho_{i\bullet}\left|-\frac{z_{i-K}+\delta\pi_{i-K}}{\left(\sum_{m=1}^{K}z_{m}+\delta\right)^{2}}+\frac{\bm{1}_{\\{j=i-K\\}}}{\sum_{m=1}^{K}z_{m}+\delta}\right|\leq\frac{\alpha\rho_{i\bullet}}{1+\delta}.$
(60)

Also, for $K+1\leq i,j\leq 2K$, we have

$\displaystyle\left|J_{i,j}(\bm{z})\right|$
$\displaystyle=\left|\frac{\partial f_{i}}{\partial
z_{j}}(\bm{z})\right|\leq\gamma(\pi_{i-K}\rho_{i,j}+\bm{1}_{\\{i=j\\}})\frac{1}{\sum_{m=K+1}^{2K}z_{m}+\delta}\leq\frac{\gamma(\pi_{i-K}\rho_{i,j}+\bm{1}_{\\{i=j\\}})}{1+\delta}.$
(61)

Let $J^{*}_{i,j}$ be the $(i,j)$-th entry of the $J^{*}$ matrix. Equations
(58)–(61) imply that

$\displaystyle\sup_{\bm{z}^{\prime}\in\mathcal{Z}}\|J(\bm{z}^{\prime})\|_{(1,\infty)}$
$\displaystyle\leq\frac{1}{1+\delta}\bigvee_{j=1}^{K}\left(\sum_{i=1}^{K}|J^{*}_{i,j}|\right)=\frac{1}{1+\delta}\|J^{*}\|_{1}.$

Therefore, as long as

$\delta>\|J^{*}\|_{1}-1,$

Equation (57) gives

$\|\bm{f}(\bm{z}_{1})-\bm{f}(\bm{z}_{2})\|_{1}\leq\frac{1}{1+\delta}\|J^{*}\|_{1}\|\bm{z}_{1}-\bm{z}_{2}\|_{1}<\|\bm{z}_{1}-\bm{z}_{2}\|_{1},$

indicating that $\bm{f}(\bm{z})=\bm{z}$ has a unique solution in
$\mathcal{Z}$.

###  6.2 Proof of Theorem 2.2

We first show the concentration of $|E(n)|/n$ around its expectation,
$\mathbb{E}[|E(n)|]/n$. Suppose we have a graph $G(n)=(V(n),E(n))$,
constructed from the PA model with heterogeneous reciprocity levels. We claim
that for $C_{0}>2$,

$\displaystyle\mathbb{P}\left(\bigl{|}|E(n)|-\mathbb{E}[|E(n)|]\bigr{|}\geq
C_{0}\sqrt{n\log n}\right)\leq 2n^{-2}.$ (62)

We prove (62) using the Azuma-Hoeffding inequality. For $k\leq n$, define
$M_{k}:=\mathbb{E}\bigl{(}|E(n)|\big{|}G(k)\bigr{)}$, then
$\mathbb{E}(M_{k+1}|G(k))=M_{k}$. So we need to consider the martingale
difference:

$M_{k+1}-M_{k}=\mathbb{E}\bigl{(}|E(n)|\big{|}G(k+1)\bigr{)}-\mathbb{E}\bigl{(}|E(n)|\big{|}G(k)\bigr{)}.$

Define another graph $G^{\prime}(n)\equiv(V^{\prime}(n),E^{\prime}(n))$ such
that $G^{\prime}(s)=G(s)$ for $s\leq k$, while $G^{\prime}(s)$ evolves
independently of $\\{G(s)\\}_{s\geq k+1}$ for $s\geq k+1$, according to the
same evolution rules as in Section 1.1. Then we have

$\displaystyle M_{k+1}-M_{k}$
$\displaystyle=\mathbb{E}\bigl{(}|E(n)|\big{|}G(k+1)\bigr{)}-\mathbb{E}\bigl{(}|E^{\prime}(n)|\big{|}G(k+1)\bigr{)}$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\bigl{(}|E(n)|-|E^{\prime}(n)|\big{|}G(k+1),G^{\prime}(k+1)\bigr{)}\middle|G(k+1)\right]$
$\displaystyle=\mathbb{E}\left[|E(k+1)|-|E^{\prime}(k+1)|\middle|G(k+1)\right],$

which gives

$|M_{k+1}-M_{k}|\leq 1.$

Since $|E(n)|-\mathbb{E}[|E(n)|]=\sum_{k=0}^{n-1}(M_{k+1}-M_{k})$, then
applying the Azuma-Hoeffding inequality gives that for $\epsilon>0$,

$\mathbb{P}\left(\bigl{|}|E(n)|-\mathbb{E}[|E(n)|]\bigr{|}\geq\epsilon\right)\leq
2e^{-\frac{\epsilon^{2}}{2n}},$

and (62) follows by setting $\epsilon=C_{0}\sqrt{n\log n}$ for $C_{0}>2$.

Next, we study the convergence of $\mathbb{E}[|E(n)|]/n$. Since $x_{m}$
satisfies (10), then we see from (8) that

$\displaystyle\mathbb{E}^{\mathcal{G}_{n}}(\Delta^{\text{in}}_{m}(n+1))$
$\displaystyle=\Delta^{\text{in}}_{m}(n)+\left(\mathbb{E}^{\mathcal{G}_{n}}\left(|E^{\text{in}}_{m}(n+1)|-|E^{\text{in}}_{m}(n)|\right)-x_{m}\right)$
$\displaystyle=\Delta^{\text{in}}_{m}(n)+\alpha\left(\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\right)$
$\displaystyle+\gamma\rho_{\bullet
m}\left(\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\right)$
$\displaystyle+\alpha\pi_{m}\sum_{r}\rho_{r,m}\left(\frac{|E^{\text{in}}_{r}(n)|+\delta|V_{r}(n)|}{|E(n)|+\delta|V(n)|}-\frac{x_{r}+\delta\pi_{r}}{\sum_{r}x_{r}+\delta}\right),$

which implies

$\displaystyle\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n+1))\right|$
$\displaystyle\leq\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|+\alpha\left|\mathbb{E}\left(\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\right)\right|$
$\displaystyle+\gamma\rho_{\bullet
m}\left|\mathbb{E}\left(\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\right)\right|$
$\displaystyle+\alpha\pi_{m}\sum_{r}\rho_{r,m}\left|\mathbb{E}\left(\frac{|E^{\text{in}}_{r}(n)|+\delta|V_{r}(n)|}{|E(n)|+\delta|V(n)|}-\frac{x_{r}+\delta\pi_{r}}{\sum_{r}x_{r}+\delta}\right)\right|.$
(63)

Also, we have that for $m=1,\ldots,K$,

$\displaystyle\left|\mathbb{E}\left(\frac{|E^{\text{in}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\right)\right|$
$\displaystyle=\left|\mathbb{E}\left(\frac{\Delta^{\text{in}}_{m}(n)}{|E(n)|+\delta|V(n)|}-\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\frac{\Delta(n)}{|E(n)|+\delta|V(n)|}\right)\right|$
$\displaystyle\leq\frac{|\mathbb{E}(\Delta^{\text{in}}_{m}(n))|}{(1+\delta)n}+\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}\frac{|\mathbb{E}(\Delta(n))|}{(1+\delta)n}.$

Similarly, we see that

$\displaystyle\left|\mathbb{E}\left(\frac{|E^{\text{out}}_{m}(n)|+\delta|V_{m}(n)|}{|E(n)|+\delta|V(n)|}-\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\right)\right|\leq\frac{|\mathbb{E}(\Delta^{\text{out}}_{m}(n))|}{(1+\delta)n}+\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\frac{|\mathbb{E}(\Delta(n))|}{(1+\delta)n}.$

Therefore, it follows from (63) that

$\displaystyle\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n+1))\right|$
$\displaystyle\leq\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|\left(1+\frac{\alpha}{(1+\delta)n}\right)+\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n))\right|\frac{\gamma\rho_{\bullet
m}}{(1+\delta)n}+$
$\displaystyle+\frac{\alpha\pi_{m}}{(1+\delta)n}\sum_{r}\rho_{r,m}\left|\mathbb{E}(\Delta^{\text{in}}_{r}(n))\right|$
$\displaystyle+\frac{\left|\mathbb{E}(\Delta(n))\right|}{(1+\delta)n}\left(\alpha\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\gamma\rho_{\bullet
m}\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}+\alpha\pi_{m}\sum_{r}\rho_{r,m}\frac{x_{r}+\delta\pi_{r}}{\sum_{r}x_{r}+\delta}\right).$

Summing over $m$ gives

$\displaystyle\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n+1))\right|\leq$
$\displaystyle\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|\left(1+\frac{\alpha}{(1+\delta)n}\right)+\sum_{m}\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n))\right|\frac{\gamma\bigvee_{m}\rho_{\bullet
m}}{(1+\delta)n}+$
$\displaystyle+\frac{\alpha\bigvee_{m}\rho_{m\bullet}}{(1+\delta)n}\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|+\frac{\left|\mathbb{E}(\Delta(n))\right|}{(1+\delta)n}\left(\alpha+C_{\delta}\right)$
$\displaystyle=$
$\displaystyle\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|\left(1+\frac{\alpha}{(1+\delta)n}\left(1+\bigvee_{m}\rho_{m\bullet}\right)\right)$
$\displaystyle+\sum_{m}\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n))\right|\frac{\gamma\bigvee_{m}\rho_{\bullet
m}}{(1+\delta)n}+\frac{\left|\mathbb{E}(\Delta(n))\right|}{(1+\delta)n}\left(\alpha+C_{\delta}\right).$
(64)

Following a similar reasoning, we have

$\displaystyle\sum_{m}\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n+1))\right|\leq$
$\displaystyle\sum_{m}\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n))\right|\left(1+\frac{\gamma}{(1+\delta)n}\left(1+\bigvee_{m}\rho_{\bullet
m}\right)\right)$
$\displaystyle+\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|\frac{\alpha\bigvee_{m}\rho_{m\bullet}}{(1+\delta)n}+\frac{\left|\mathbb{E}(\Delta(n))\right|}{(1+\delta)n}\left(\gamma+C_{\delta}\right),$
(65)

and

$\displaystyle\sum_{m}$ $\displaystyle\left|\mathbb{E}(\Delta(n+1))\right|$
(66) $\displaystyle\leq$
$\displaystyle\left|\mathbb{E}(\Delta(n))\right|\left(1+\frac{\alpha}{(1+\delta)n}\sum_{m}\rho_{m\bullet}\frac{x_{m}+\delta\pi_{m}}{\sum_{r}x_{r}+\delta}+\frac{\gamma}{(1+\delta)n}\sum_{m}\rho_{\bullet
m}\frac{y_{m}+\delta\pi_{m}}{\sum_{r}y_{r}+\delta}\right)$
$\displaystyle+\sum_{m}\left|\mathbb{E}(\Delta^{\text{in}}_{m}(n))\right|\frac{\alpha\bigvee_{m}\rho_{m\bullet}}{(1+\delta)n}+\sum_{m}\left|\mathbb{E}(\Delta^{\text{out}}_{m}(n))\right|\frac{\gamma\bigvee_{m}\rho_{\bullet
m}}{(1+\delta)n}.$ (67)

Combining Equations (64), (65) and (67) gives:

$\displaystyle\begin{bmatrix}\sum_{m}|\mathbb{E}(\Delta^{\text{in}}_{m}(n+1))|\\\
\sum_{m}|\mathbb{E}(\Delta^{\text{out}}_{m}(n+1))|\\\
|\mathbb{E}(\Delta(n+1))|\end{bmatrix}\leq\left(\bm{I}+\frac{1}{n}\bm{H}\right)\begin{bmatrix}\sum_{m}|\mathbb{E}(\Delta^{\text{in}}_{m}(n))|\\\
\sum_{m}|\mathbb{E}(\Delta^{\text{out}}_{m}(n))|\\\
|\mathbb{E}(\Delta(n))|\end{bmatrix}.$

Then by the Perron-Frobenius theorem, we see that

$\frac{1}{n}\sum_{m}|\mathbb{E}(\Delta^{\text{in}}_{m}(n))|\to
0,\quad\frac{1}{n}\sum_{m}|\mathbb{E}(\Delta^{\text{out}}_{m}(n))|\to
0,\quad\text{and}\quad\frac{1}{n}|\mathbb{E}(\Delta(n))|\to 0,$

if the largest eigenvalue of $\bm{H}$, $\lambda_{H}$, is less than 1, thus
completing the proof of the theorem.

##  Appendix A Generalized Breiman’s Theorem

######  Theorem A.1.

Suppose $\\{\bm{\xi}(t):t\geq 0\\}$ is an $\mathbb{R}_{+}^{p}$-valued
stochastic process for some $p\geq 1$. Let $X$ be a positive random variable
with regularly varying distribution satisfying for some scaling function
$b(t)$,

$\lim_{t\to\infty}t\mathbb{P}(X/b(t)>x)=x^{-c}=:\nu_{c}\bigl{(}(x,\infty)\bigr{)},\quad
x>0,c>0.$

Further suppose

  1. 1.

For some finite and positive random vector $\bm{\xi}_{\infty}$,

$\lim_{t\to\infty}{\bm{\xi}(t)}=\bm{\xi}_{\infty}\quad(\text{almost surely});$

  2. 2.

The random variable $X$ and the process $\bm{\xi}(\cdot)$ are independent.

Then:

(i) In $\mathbb{M}(\mathbb{R}_{+}^{p}\times(\mathbb{R}_{+}\setminus\\{0\\}))$,

$t\mathbb{P}\Bigl{[}\Bigl{(}{\bm{\xi}(X)},\frac{X}{b(t)}\Bigr{)}\in\cdot\,\Bigr{]}\longrightarrow\mathbb{P}(\bm{\xi}_{\infty}\in\cdot\,)\times\nu_{c}(\cdot)=:\eta(\cdot).$
(68)

If $\bm{\xi}_{\infty}$ is of the form $\bm{\xi}_{\infty}=:L\bm{v}$ where $L>0$
almost surely and $\bm{v}\in(0,\infty)^{p}$, then $\eta(\cdot)$ concentrates
on the subcone $\mathcal{L}\times(\mathbb{R}_{+}\setminus\\{0\\})$ where
$\mathcal{L}=\\{\theta\bm{v}:\theta>0\\}$.

(ii) If additionally, for some $c^{\prime}>c$ we have the condition

$\kappa:=\sup_{t\geq
0}\mathbb{E}\left[\Bigl{(}{\|\bm{\xi}(t)\|}\Bigr{)}^{c^{\prime}}\right]<\infty,$
(69)

for some $L_{p}$ norm $\|\cdot\|$, then the product of components in (68),
$\bm{\xi}(X)X$, has a regularly varying distribution with scaling function
$b(t)$ and in $\mathbb{M}(\mathbb{R}_{+}^{p}\setminus\\{\bm{0}\\})$,

$t\mathbb{P}\Bigl{[}\frac{X\bm{\xi}(X)}{b(t)}\in\cdot\,\Bigr{]}\longrightarrow\left(\mathbb{P}(\bm{\xi}_{\infty}\in\cdot\,)\times\nu_{c}\right)\circ
h^{-1},$ (70)

where $h(\bm{y},x)=x\bm{y}$.

For the classical Breiman Theorem where $p=1$ and
$\bm{\xi}(t)\equiv\bm{\xi}_{\infty}$, (69) is the expected moment condition.

## References

  * [1] [author] Athreya, K. B.K. B., Ghosh, A. P.A. P. and Sethuraman, S.S. (2008). Growth of preferential attachment random graphs via continuous-time branching processes. Proceedings Mathematical Sciences 118 473–494.
  * [2] [author] Athreya, K.K. and Ney, P.P. (1972). Branching Processes. Springer-Verlag, New York.
  * [3] [author] Athreya, K. B.K. B. and Ney, P.P. (2004). Branching processes. Reprint of the 1972 original. Springer, New York.
  * [4] [author] Basrak, B.B. and Planinić, H.H. (2019). A note on vague convergence of measures. Statist. Probab. Lett. 153 180–186. 10.1016/j.spl.2019.06.004 3979308
  * [5] [author] Bollobás, B.B., Borgs, C.C., Chayes, J.J. and Riordan, O.O. (2003). Directed scale-free graphs. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms (Baltimore, 2003) 132-139. ACM, New York.
  * [6] [author] Cirkovic, D.D., Wang, T.T. and Resnick, S. I.S. I. (2022). Preferential Attachment with Reciprocity: Properties and Estimation. arXiv preprint arXiv:2201.03769.
  * [7] [author] Das, B.B., Mitra, A.A. and Resnick, S.S. (2013). Living on the multi-dimensional edge: Seeking hidden risks using regular variation. Advances in Applied Probability 45 139–163.
  * [8] [author] Das, B.B. and Resnick, S. I.S. I. (2015). Models with hidden regular variation: Generation and detection. Stochastic Systems 5 195-238. 10.1214/14-SSY141
  * [9] [author] Das, B.B. and Resnick, S. I.S. I. (2017). Hidden regular variation under full and strong asymptotic dependence. Extremes 20 873–904.
  * [10] [author] Hult, H.H. and Lindskog, F.F. (2006). Regular variation for measures on metric spaces. Publ. Inst. Math. (Beograd) (N.S.) 80(94) 121–140. 10.2298/PIM0694121H 2281910 (2008g:28016)
  * [11] [author] Jiang, B.B., Zhang, Z. L.Z. L. and Towsley, D.D. (2015). Reciprocity in Social Networks with Capacity Constraints. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. KDD ’15 457–466. Association for Computing Machinery, New York, NY, USA. 10.1145/2783258.2783410
  * [12] [author] Kilduff, M.M., Tsai, W.W. and Hanke, R.R. (2006). A paradigm too far? A dynamic stability reconsideration of the social network research program. Academy of Management Review 31 1031–1048.
  * [13] [author] Kirk, W. A.W. A. (2001). Handbook of Metric Fixed Point Theory, 1 ed. Springer Netherlands, Dordrecht.
  * [14] [author] Krapivsky, P. L.P. L. and Redner, S.S. (2001). Organization of growing random networks. Physical Review E 63 066123:1–14.
  * [15] [author] Kulik, R.R. and Soulier, P.P. (2020). Heavy-Tailed Time Series. Springer Series in Operations Research and Financial Engineering. Springer, New York, NY.
  * [16] [author] Lindskog, F.F., Resnick, S. I.S. I. and Roy, J.J. (2014). Regularly varying measures on metric spaces: Hidden regular variation and hidden jumps. Probab. Surv. 11 270-314. 10.1214/14-PS231
  * [17] [author] Matias, C.C. and Miele, V.V. (2017). Statistical clustering of temporal networks through a dynamic stochastic block model. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 79 1119–1141.
  * [18] [author] Molm, L. D.L. D., Collett, J. L.J. L. and Schaefer, D. R.D. R. (2007). Building solidarity through generalized exchange: A theory of reciprocity. American journal of sociology 113 205–242.
  * [19] [author] Rabehasaina, L.L. and Woo, J. K.J. K. (2021). Multitype branching process with nonhomogeneous Poisson and contagious Poisson immigration. Journal of Applied Probability.
  * [20] [author] Resnick, S. I.S. I. (2007). Heavy Tail Phenomena: Probabilistic and Statistical Modeling. Springer Series in Operations Research and Financial Engineering. Springer-Verlag, New York. ISBN: 0-387-24272-4.
  * [21] [author] Resnick, S. I.S. I. and Samorodnitsky, G.G. (2015). Tauberian theory for multivariate regularly varying distributions with application to preferential attachment networks. Extremes 18(3) 349–367. 10.1007/s10687-015-0216-2
  * [22] [author] Samorodnitsky, G.G., Resnick, S.S., Towsley, D.D., Davis, R.R., Willis, A.A. and Wan, P.P. (2016). Nonstandard regular variation of in-degree and out-degree in the preferential attachment model. Journal of Applied Probability 53(1) 146–161. 10.1017/jpr.2015.15
  * [23] [author] Wan, P.P., Wang, T.T., Davis, R. A.R. A. and Resnick, S. I.S. I. (2017). Fitting the linear preferential attachment model. Electron. J. Statist. 11 3738-3780. 10.1214/17-EJS1327
  * [24] [author] Wang, T.T. and Resnick, S. I.S. I. (2020). Degree growth rates and index estimation in a directed preferential attachment model. Stochastic Processes and their Applications 130 878–906.
  * [25] [author] Wang, T.T. and Resnick, S. I.S. I. (2022). Measuring Reciprocity in a Directed Preferential Attachment Network. Advances in Applied Probability 1–25. 10.1017/apr.2021.52
  * [26] [author] Wang, T.T. and Resnick, S. I.S. I. (2022). Asymptotic Dependence of In-and Out-Degrees in a Preferential Attachment Model with Reciprocity. Extremes 25 417–450. 10.1007/s10687-022-00439-5