{'fraction_non_alphanumeric': 0.06715237302248127, 'fraction_numerical': 0.019567027477102414, 'mean_word_length': 3.8783509341998377, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Random Networks with Heterogeneous Reciprocity

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.


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 [21,22,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 [6,26] 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 ρ m,r ∈ (0, 1), m, r ∈ {1, . . . , 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 inand out-degrees in Section 4, and give concluding remarks in Section 5. Technical proofs of results in Section 2 are collected in Section 6.

PA Model with Heterogeneous Reciprocity
The proposed model extends the directed preferential attachment (PA) model studied in [5,14,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 ≥ 1} are iid random variables, independent from the growth mechanism of the graph, with P(W v = r) = π r , for K r=1 π r = 1.
We imagine that when a node is born, it flips a multi-sided coin to determine its communication type. Think of π 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 ∈ 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 ≥ 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) = {(1, 1)}. 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 ρ m,r ∈ [0, 1], for m, r ∈ {1, 2, . . . , K}. Here ρ 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 ρ := (ρ m,r ) m,r is not necessarily a stochastic matrix, but can be an arbitrary matrix in M K×K ([0, 1]), the set of all K × K matrices with entries belonging to [0, 1]. Later in Section 2, we will give particular regularity conditions on ρ to facilitate theoretical analysis.
Let D in v (n), D out v (n) be the in-and out-degrees of node v ∈ V (n) in G(n), and we use the convention that . Let δ > 0 be an offset parameter. The evolution of the network G(n + 1) from G(n) is described as follows. 1. With probability α ∈ (0, 1), add a new node |V (n)| + 1 with a directed edge (|V (n)| + 1, v), where v ∈ V (n) is chosen with probability and update the node set V (n + 1) = V (n) ∪ {|V (n)| + 1} and W (n + 1) = W (n) ∪ {W |V (n)|+1 }. The new node |V (n)| + 1 belongs to group r with probability π r . If node v belongs to group m, then a reciprocal edge (v, |V (n)| + 1) is added with probability ρ m,r . Upon reciprocation, update the edge set as E(n and update the node set V (n The new node |V (n)| + 1 belongs to group r with probability π r . If node v belongs to group m, then a reciprocal edge (|V (n)| + 1, v) is added with probability ρ r,m . Upon reciprocation, update the edge set as E(n Note that |V (n)| = n + 1, n ≥ 0, since a new node is added at each step, and the offset parameter δ is assumed to be the same for both in-and out-degrees.

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. Let G n be denote the σ-field generated by observing the network evolution up to n steps, i.e.
For v ∈ V (n), we have Let V m (n) := {v ∈ V (n) : W v = m}, and note that and applying (3) and (5) gives For the third term in (6), we have that (v, |V (n)| + 1) is not due to reciprocation Therefore, combining (6) and (7) shows that Following a similar argument, we use the relationship in (4) to obtain that In addition, summing over m in (8) gives As remarked at the beginning of Section 2, the proof of Theorem 2.2 requires showing the convergence of E[|E(n)|]/n and for this it suffices to check E[|E in m (n)|]/n → x m and E[|E out m (n)|]/n → y m , for some x m , y m ∈ [0, 1], which is implied by the approximation E Gn |E in m (n + 1)| − |E in m (n)| ≈ x m , and E Gn (|E out m (n + 1)| − |E out m (n)|) ≈ y m . Therefore, we expect x m , y m , m = 1, . . . , K, to be solutions to the following system of equations: for m = 1, . . . , K. Summing over m in (10) and (11) shows K m=1 x m = K m=1 y m and from (10) we get and therefore K m=1 y m ≥ 1. Also, since x m + δπ m r x r + δ ≤ 1, and y m + δπ m r y r + δ ≤ 1, the right hand side of (10) is upper bounded by 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 , . . . , x K , y 1 , . . . , y K ) in The proof relies on the contraction mapping theorem [13, Theorem 1.2.2]. We proceed by defining a function f : Z → Z, using the right-hand-side expressions in (10) and (11). Applying the mean value theorem to f , we deduce a sufficient condition for 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 ⊙ and ⊗ denote the Hadamard and Kronecker products of matrices, respectively. Denote the identity matrix in R K as I K , and use 1 K to denote the vector in the system of equations in (10) and (11) has a unique solution in 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 in (n)|/n and |E out (n)|/n. Define a constant and a matrix The proof of Theorem 2.  (14),   m |E(∆ in m (n + 1))| m |E(∆ out m (n + 1))| |E(∆(n))| Hence, provided that all elements in H are strictly positive, applying the Perron-Frobenius theorem suggests restricting the largest eigenvalue of H will be sufficient to control the growth of the expected edge counts. In addition, the assumption of positive elements in H leads to the requirement of α, γ > 0, m ρ m• > 0, and m ρ •m > 0.
Recall we may interpret ρ m• and ρ •m , them as the likelihood of sending and attracting reciprocal edges for a user in group m, respectively. The assumption m ρ m• > 0 and m ρ •m > 0 requires that there exists at least one group whose probability of generating or attracting reciprocal edges is strictly positive.
is the solution to the system of equations given in (10) and (11). We have: Note that since (x m , y m ) satisfies (10) and (11), the a.s. limit of |E(n)|/n depends on the value of δ 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:

Growth of Degree Counts
We next focus on the asymptotic behavior of the joint inand out-degree counts: We extend the homogeneous reciprocity techniques in [6,26] 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).

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].
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 ξ δ (·, m) as a MBI process with group label m. The process ξ δ (·, 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 ξ (1) δ (·, m) and ξ (2) δ (·, m) to be α and γ ≡ 1 − α, respectively. For s = (s 1 , s 2 ) ∈ [0, 1] 2 , the branching structure of ξ δ (·, m) is specified through offspring generating functions: According to (16), when a group-m type I particle's lifetime ends, with probability 1 − ρ m• , it splits into two group-m type I particles, increasing the total number of group-m type I particles by 1, and with probability ρ m• , 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 δ > 0. When an immigration event happens, ξ δ (·, m) is incremented by (1, 0), (0, 1), or (1, 1) following the distribution 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 ξ δ (·, m): Applying the Perron-Frobenius theorem, A m has a largest positive eigenvalue with multiplicity 1, i.e.
Also, the smaller eigenvalue of A m is which plays an important role in the discussion of hidden regular variation in Theorem 4.2. Let be the left and right eigenvectors associated with λ m respectively, with all coordinates strictly positive, and u(m) Applying Theorem 1 in [26] gives that there exists some finite positive random variable Z(m) such that 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 and for t ≥ 0, the generating function of ξ δ (t, L) is Imagine L assigns a label to a MBI which remains unchanged throughout the evolution of ξ δ (·, ·). The initial value, ξ δ (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 (ξ δ (0, L), L), ξ δ (·, L) behaves like a MBI process with fixed group label L. In the sequel, we refer to {ξ δ (t, L) : t ≥ 0} as a two-type MBI process with random group label.

3.2.
Embedding Let {L v : v ≥ 1} be random variables with common pmf's as in (22). We now explain how to embed the degree sequence We assume {ξ v,δ (·, ·)} v≥1 have the same parameters but possibly different initializations.
k,δ (T n − T k−1 ) = n + 1 + n k=1 R k , we then have Hence, L n+2 and T n+1 − T n are independent under P FT n , and since In addition, from the four cases listed above, we also see that Similarly, we have Therefore, under P FT n , ξ n+2,δ (0, L n+2 ) is a random vector following the distribution Meanwhile, since {L k : 1 ≤ k ≤ n + 1} are F Tn -measurable, we then have Therefore, R n+1 , J n+1 is independent from T n+1 − T n under P FT n .

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.
where the second term goes to 0 a.s. as n → ∞. Then we only need to consider the limit of the first term in (36).

Then we have
From the proof of Theorem 2.2, we see that there exist constants C 1 (δ), C 2 (δ) > 0 such that Then the left hand side of (38) is bounded by which goes to 0 as n → ∞, thus proving the claim in (38).
Also, by the definition of ξ w,δ , w ≥ 1, we see that Since the function P[ ξ δ (t, L * ) = (k, l)] is bounded and continuous in t, then we conclude that A 4 (n) p −→ 0 by applying the Riemann integrability of P ξ δ (− log t/c * , L * ) = (k, l) , which completes the proof of (34).

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 P[(I m , O m ) ∈ ·], for each m = 1, . . . , K.

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 C 0 ⊂ C ⊂ R 2 + are two closed cones, and we provide the definition of Mconvergence in Definition 4.1 (cf. [4,7,10,15,16]) on C \ C 0 , which lays the theoretical foundation of regularly varying measures (cf. Definition 4.2). DEFINITION 4.1. Let M(C \ C 0 ) be the set of Borel measures on C \ C 0 which are finite on sets bounded away from C 0 , and C(C \ C 0 ) be the set of continuous, bounded, nonnegative functions on C \ C 0 whose supports are bounded away from C 0 . Then for µ n , µ ∈ Without loss of generality [16], we can and do take functions in C(C \ C 0 ) to be uniformly continuous as well. Denote the modulus of continuity of a uniformly continuous function f : where d(·, ·) is an appropriate metric on the domain of f . Uniform continuity means lim δ→0 ∆ f (δ) = 0. We now present the definition of multivariate regular variation with C = R 2 + and C 0 = {0}. Following Definition 2.1 in [20], we denote a regularly varying function f : R + → R + with index a ∈ R, as f ∈ RV a . Definition 4.2 gives the formal description of the MRV of distributions.
DEFINITION 4.2. The distribution P(Z ∈ ·) of a random vector Z on R 2 + , is (standard) regularly varying on R 2 + \ {0} with index c > 0 if there exists some regularly varying scaling function b(t) ∈ RV 1/c and a limit measure ν(·) ∈ M(R 2 + \ {0}) such that as t → ∞, When analyzing the asymptotic dependence between components of a bivariate random vector 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 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 Z is regularly varying on R 2 + \ {0} and has hidden regular variation on R 2 + \ C 0 if there exist 0 < c ≤ c 0 , scaling functions b(t) ∈ RV 1/c and b 0 (t) ∈ RV 1/c0 with b(t)/b 0 (t) → ∞ and limit measures ν, ν 0 , such that A convenient way to characterize HRV is through the generalized polar coordinate transformation for R 2 + \ C 0 and an associated metric d(·, ·) satisfying d(cx, cy) = cd(x, y) for scalars c > 0. The metric d(·, ·) that we use in practice is the L 1 -metric. When using generalized polar coordinates with respect to the forbidden zone C 0 , we define ℵ C0 := {x ∈ C \ C 0 : d(x, C 0 ) = 1}, the locus of points at distance 1 from C 0 . Then the generalized polar coordinates are specified through the transformation, GPOLAR : Let ν c (·) be a measure in M(R + \ {0}) satisfying ν c (x, ∞) = x −c x, c > 0, and S 0 (·) be a probability measure on ℵ C0 . Then generalized polar coordinates allow re-writing (42) as . See [7] and [16] for details. 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 ξ δ (·, m) is specified through the matrix A m given in (19), and the largest eigenvalue of A m is λ m given in (20). We assume parameters α, γ and the matrix ρ, are chosen such that λ m , m = 1, . . . , K, are all distinct values and without loss of generality, assume the behavioral group labels are chosen so that (43) λ 1 > λ 2 > · · · > λ K , Theorem 4.1 gives the multivariate regular variation of P[(I, O) ∈ · ] on R 2 + \ {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 λ 1 ≥ log 2 to guarantee that the moment condition in (69) is satisfied.
Theorem 4.1 shows that the limiting measure µ 1 concentrates Pareto mass on the ray L (1) := {(x, y) ∈ R 2 + : y = a(1)x} and concentration on L (1) suggests that at scale b(t) = t λ1/c * , large in-and out-degree pairs in the PA model with heterogeneous reciprocity levels satisfy O ≈ a(1)I.
PROOF. Similar to the proof of Theorem 6 in [26], we first see that P( Z 1 > 0) = 1 by the property of MBI processes. Also, since λ 1 > λ r for r = 2, . . . , K, Therefore, Then the proof of Theorem 4.1 is an application of Theorem A.1 after making the identifications The remaining piece is to show the moment condition (69) in this context and we will show any δ ≥ 0 and any q = 1, 2, . . ., there exists some constant C(δ, q) > 0 such that which is true by Proposition 2 in [26].
Next, Theorem 4.2 gives a second hidden regular variation (HRV) regime after removing L (1) [7,8,9,16,20]. 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 L (1) . Thus, we may seek a regular variation property on R 2 + \ 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 where d 1 (x, y) is a metric on R 2 + \ {0} chosen for convenience to be the L 1 -metric. The Define ℵ L(1) := {x ∈ R 2 + \ L (1) : d ′ 1 (x, L (1) ) = 1}, which are 2 line segments in R 2 + parallel to L (1) . Hidden regular variation will be present for (I, O) converges to a limit measure in M((R + \ {0}) × ℵ L(1) ). The next theorem explains the convergence on M(R 2 + \ L (1) ).

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 outdegree 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.
Let J * i,j be the (i, j)-th entry of the J * matrix. Equations (58)-(61) imply that Therefore, as long as δ > J * 1 − 1, Equation (57) gives indicating that f (z) = z has a unique solution in Z.