Autoregressive Networks with Dependent Edges
Jinyuan Chang Qin Fang Eric D. Kolaczyk Peter W. MacDonald Qiwei Yao





1,2]Jinyuan Chang
3]Qin Fang
4]Eric D. Kolaczyk
4]Peter W. MacDonald
5]Qiwei Yao
[1]Joint Laboratory of Data Science and Business Intelligence, Southwestern University of Finance and Economics, Chengdu, Sichuan 611130, China

[2]Academy of
Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
[3]Business School, University of Sydney, Sydney, 2008,  Australia
[4]Department of Mathematics and Statistics,
McGill University,
Montreal, Quebec H3A 0B9, Canada
[5]Department of Statistics, London School of Economics and Political Science, London, WC2A 2AE, U.K.




Autoregressive Networks with Dependent Edges
    [
    
============================================




We propose an autoregressive framework for modelling dynamic networks with dependent edges. It encompasses the models which accommodate, for example, transitivity, density-dependent and other stylized features often observed in real network data. By assuming the edges of network at each time are independent conditionally on their lagged values, the models, which exhibit a close connection with temporal ERGMs, facilitate both simulation and the maximum likelihood estimation in the straightforward manner. Due to the possible large number of parameters in the models, the initial MLEs may suffer from slow convergence rates. An improved estimator for each component parameter is proposed based on an iteration based on the projection which mitigates the impact of the other parameters <cit.>. Based on a martingale difference structure, the asymptotic distribution of the improved estimator is derived without the stationarity 
assumption. The limiting distribution is not normal in general, and it reduces to normal when the underlying process satisfies some mixing conditions. Illustration with a transitivity model was carried out in both simulation and a real network data set.


Key words: conditional independence, dynamic networks, maximum likelihood estimation, stylized features of network data, transitivity.





§ INTRODUCTION



Dynamic network modelling with dependent edges is  practically important and relevant but technicall challenging. Without dependence across different edges, it is impossible to incorporate into the models some stylized features often observed in real network data such as transitivity, density dependence. On the other hand, dependent edges make the dynamic structures of network processes complex and statistical inference challenging. Existing literature on modeling dynamic networks with dependent edges can be divided into two categories: latent process based models <cit.>, and temporal
exponential family random-graph models (ERGMs)  <cit.>. Inference and simulation for the latent
process based models rely on compute-intensive methods such as MCMC and EM algorithms, as their likelihood functions are not explicitly available. Furthermore the dependence depicted by latent process based models are implicit, and it cannot be configured easily to accommodate  stylized features of real network data. On the other hand, it is well documented that ERGMs without a proper control on the level of dependence may suffer from lack of computational scalability, instable estimation algorithms, and cancentration on extreme subspaces of graph space. See <cit.> and the references within. More recently <cit.> proposed a block logistic autoregressive network model with dependent edges, which was fitted using an EM algorithm. 

In this paper, we propose a new autoregressive based Markov chain model with dependent edge. Following <cit.>, we specify the transition probabilities of forming a new edge or dissolving an existing edge  between each pair of nodes explicitly depending on its history. Furthermore we allow those probabilities depending on the histories of other edge processes. This enlarged form of the transition probabilities make the model flexible enough to accommodate the stylized
features such as transitivity and density dependence. The resulting network processes have dependent edges, which is radically different from those considered in <cit.>.
Similar to <cit.>, <cit.> and <cit.>, we assume that the edges are conditionally independent given their joint histories. This makes both statistical inference and theoretical analysis more transparent. This conditional independence avoids the difficulties caused by the normalized constants in ERGMs; see <cit.> and <cit.>. Based on the conditional independence, we can build up a martingale difference structure which facilitates the asymptotic analysis for the maximum likelihood estimation (see Section <ref> below).
This, to our best knowledge, has never been done before in the context of dynamic networks with dependence edges.

The rest of the paper is organized as follows. Section <ref> presents the general AR network framework with dependent edges. 

We also discuss the relationship between the proposed AR models and temporal ERGMs. Section <ref> contains three concrete AR models which are designed to model, respectively, density-dependence, persistence, and transitivity — those are among the stylized features often observed in real network data. The  
two versions of the maximum likelihood estimation (MLE) for the parameters in the AR models  and the associated asymptotic theory are presented in Section <ref>. We introduce the concepts of local parameters and global parameters, which need to be identified and estimated separatly. They may also entertain different convergence rates. The initial MLE suffers from slower convergence rates due to the diverging number of local parameters. An improved MLE for each component parameter is obtained by projecting the score function onto the corresponding direction, which mitigates the impact of the other parameters <cit.>. Based on a martingale difference structure, the asymptotic distribution of the improved estimator is derived without the stationarity assumption. The limiting distribution is not normal in general. But it reduces to normal when the underlying process satisfies some mixing conditions which holds for many stationary processes. Section <ref> presents an illustration via simulation for the proposed transitivity model.
Further illustration with a real dynamic network data set is reported in Section <ref>.
  An online supplementary contains further more numerical results and all the technical proofs. 

Notation.  For any positive integer r, write [r] = {1,…, r} and ℝ_+^r={(x_1,…,x_r)^:x_i>0 for any i∈[r]}. For any x, y ∈ℝ, we write x ∨ y = max(x,y).  For two positive sequences {a_n} and {b_n}, we write a_n≪ b_n if  limsup_n →∞ a_n/b_n = 0, and 
 a_n ≲ b_n if limsup_n →∞ a_n/b_n <∞. 
For any vector = (b_1, …, b_r)^∈ℝ^r, we let _-l denote the sub-vector of  by removing the l-th component b_l. Given an index set ℳ⊂ [r], we let _ℳ denote the sub-vector of  that consists of the components of  indexed by ℳ.
For any r_1 × r_2 real matrix , denote by ^ its transpose. When r_1 = r_2, we use λ_min() to denote  the smallest eigenvalue of the square matrix . For any set
𝒰, |𝒰| denotes its cardinality.




§ AR(M) NETWORK FRAMEWORK
 


 §.§ Model


Consider a dynamic network process defined on  p nodes denoted by 1, …, p. Let _t ≡ (X_i,j^t)_p× p be the adjacency matrix at time t, where X_i,j^t =1 denotes the existence of an edge between nodes i and j at time t, and 0 otherwise. For simplicity, we only consider undirected networks without self-loops, i.e. X_i,i^t ≡ 0 and X_i,j^t= X_j,i^t. The main idea can be applied to directed networks.

 
 Conditionally on {_s}_s≤ t-1,
the edges {X_i,j^t}_1≤ i <j≤ p are mutually independent with 

    α_i,j^t-1   ≡ℙ(X_i,j^t=1 | X_i,j^t-1 =0,  _t-1∖ X_i,j^t-1, _t-k  for k≥ 2)
       = ℙ(X_i,j^t=1 | X_i,j^t-1 =0,  _t-1∖ X_i,j^t-1, _t-2, …, _t-m) ,
    β_i,j^t-1   ≡ℙ(X_i,j^t=0 | X_i,j^t-1 =1,  _t-1∖ X_i,j^t-1, _t-k  for k≥ 2)
       
    = ℙ(X_i,j^t=0 | X_i,j^t-1 =1,  _t-1∖ X_i,j^t-1,_t-2, …, _t-m) ,

where m≥ 1 is an integer. 


An AR(m) network process defined above is a Markov chain with
order m. Based on (<ref>) and (<ref>), we have 

    ℙ(X_i,j^t=1 | _t-1, …, _t-m) =
    α_i,j^t-1 +X_i,j^t-1 (1 - α_i,j^t-1 -β_i,j^t-1) ≡γ_i,j^t-1 ,

which implies that 

    X_i,j^t | _t-1,…,_t-m∼ Bernoulli(γ_i,j^t-1) ,    1≤ i<j≤ p .

Clearly 
edges
X_i,j^t, for different (i,j), are not independent with each other. We may impose various forms for the conditional probabilities α_i,j^t-1 and β_i,j^t-1 to reflect different stylized features of network data. Put 

    α_i,j^t-1 = f_i,j( _t-1∖ X_i,j^t-1, _t-2,…, _t-m; _0) ,
    β_i,j^t-1 = g_i,j (_t-1∖ X_i,j^t-1, _t-2, …, _t-m; _0) ,

where f_i,j's and g_i,j's are known functions, and _0∈⊂ℝ^q is a q-dimensional unknown true parameter vector. For any ∈, write 

    α_i,j^t-1() = f_i,j(_t-1∖ X_i,j^t-1, _t-2, …, _t-m; ) ,
    β_i,j^t-1() = g_i,j (_t-1∖ X_i,j^t-1, _t-2, …, _t-m; ) .

Then α_i,j^t-1=α_i,j^t-1(_0) and β_i,j^t-1=β_i,j^t-1(_0).

 Modelling dynamic networks by Markov or/and AR models is not new. See, for example, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. However, most available Markov models are designed for  Erdös-Renyi networks with independent edges. Our setting provides a general framework to accommodate various dependence structures across different edges. Some practical network models satisfy this general framework are introduced in Section <ref>.




For the special AR(1) processes (i.e. m=1), if both f_i,j and g_i,j in (<ref>) are always positive and smaller than 1 for all 1≤ i < j≤ p, {_t}_t≥1 is an irreducible homogeneous Markov chain with 2^p(p-1)/2 states. Therefore when p is  fixed,  (i) there exists a unique stationary distribution, and (ii)  if _0 is activated according to this stationary distribution, the process {_t}_t≥ 1 is strictly stationary and ergodic. See Theorems 3.1, 3.3 and 4.1 in Chapter 3 of <cit.>.  
Hence the density-dependent model introduced in Section <ref> 
and the transitivity model introduced in Section <ref> are strictly stationary for any fixed constant p if all the transition probability functions α_i,j^t-1 and β_i,j^t-1 are strictly between 0 and 1. 
It is worth pointing out that the ergodicity only holds for any fixed constant p. Hence we cannot take for granted that the sample means of _t and/or its summary statistics converge when p diverges together with the sample size, even when _t is stationary. Note that  stationarity is  not an asymptotic
property while ergodicity is.  





 §.§ Relationship to temporal ERGMs



Similar to temporal ERGMs explored in <cit.> and <cit.>, we assume that the edges are conditionally independent given their lagged values. But instead of specifying some exponential family distributions as the transition probabilities, we define, separately,  the probability for forming a new edge in (<ref>), and that for dissolving an existing edge in (<ref>). Those two probability functions can be in any desirable forms as presented in (<ref>).
This allows us to impose some closed-forms of parametric functions for α_i,j^t-1 and β_i,j^t-1, and those 
functions need only to be between 0 and 1. Hence, the likelihood functions are explicitly available, 

which allows us to depict more explicitly in our models some stylized features often observed in real network data. See Section <ref> for details. The numerical analysis with both simulated and real data 
in Sections <ref> and <ref> indicates that the proposed AR models are capable to simulate and to reflect some observed interesting dynamic network phenomena.

The temporal ERGMs with conditional independent edges <cit.> can be expressed in our AR(m) framework as

    α_i,j^t-1() =    exp{()^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}/
    1+ exp{()^_i,j( X_t-1∖ X_i,j^t-1, 
    X_t-2, …,  X_t-m )} , 
    β_i,j^t-1() =    exp{()^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}/1+ exp{()^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )} ,

where _i,j(·) and _i,j(·) have closed-form expressions based on the sufficient statistics of the original temporal ERGM. Furthermore, if () and () in the above expressions are replaced by, respectively, 
(_α) and (_β), where _α
and _β are two sets of different parameters, _t follows
the separable temporal ERGM with conditional independent
edges given the past networks <cit.>. See
 Section <ref> of the supplementary material for the detailed discussion on the relationship of the proposed AR models and temporal ERGMs. 




§ SOME INTERESTING AR NETWORK MODELS
 
To illustrate the usefulness of the AR(m) framework proposed above, we state three AR(m) network models which
reflect various stylized features in real network data.  In all three models, the parameters {ξ_i}_i=1^p and
{η_i}_i=1^p reflect node heterogeneity in, respectively, forming a new edge and dissolving an existing edge. Specifically, the larger ξ_i is, the more likely node i will form new edges with other nodes, and the larger η_i is, the more likely the existing edges between node i and the others will be dissolved.
Instances of these three models can be simulated using our development R package , available at <https://github.com/peterwmacd/arnetworks>.



 §.§ Density-dependent model
  
Let ϑ_i,j^t-1 = exp{ a_0 D_-i,-j^t-1
    + a_1 (D_i^t-1 +  D_j^t-1)} and ϖ_i,j^t-1 =  exp{ b_0 (1-D_-i,-j^t-1)
    + b_1 (2-D_i^t-1- D_j^t-1)}
with

    D_-i,-j^t-1= 1/(p-2)(p-3)∑_k,ℓ:  k,ℓ≠ i,j,  k ≠ℓ X_k, ℓ^t-1  and  D_i^t-1= 1/p-1∑_ℓ: ℓ i X_i, ℓ^t-1 ,

where D_i^t-1 and
D_j^t-1 are, respectively, the densities of node i and node j at time t-1, and D_-i,-j^t-1 is the network density excluding nodes i and j at time t-1.
We specify the transition probabilities 
as follows:

    α_i,j^t-1() =  ξ_i ξ_j ϑ_i,j^t-1 1 +ϑ_i,j^t-1 + ϖ_i,j^t-1 ,
            β_i,j^t-1() = η_i η_j ϖ_i,j^t-1 1 +ϑ_i,j^t-1 + ϖ_i,j^t-1 .
 
This is an AR(1) model with parameter vector = (a_0, a_1,  b_0, b_1,  ξ_1, …, ξ_p, η_1, …, η_p)^∈⊂ℝ_+^2p+4. 
 Then the propensity to form a new edge between nodes i and j at time t is positively impacted by 
D_-i,-j^t-1, D_i^t-1 and D_j^t-1, and the propensity to dissolve an existing edge between nodes i and j at time t is negatively impacted by these three densities. 

<cit.> proposed an ERGM with network density in its index function. In (<ref>) we explicitly specify the impact from the density functions on forming a new edge and dissolving an existing edge, while the model in Section  2.1 of <cit.> does not differentiate the representations for these two types of impact. Within a separable ERGM framework, the edge counts model of <cit.> assumes that the collection of all newly formed edges is conditionally independent of the collection of all newly dissolved edges given their history, and the two conditional distributions are controlled by different parameters.



 §.§ Persistence model
 
We define the transition probabilities 

 
 



    α_i,j^t-1()=    ξ_i ξ_j  exp [-  1 - a{ (1-X_i,j^t-2) +(1-X_i,j^t-2)(1-X_i,j^t-3) } ] ,
    β_i,j^t-1() =     η_i η_j exp{-  1 - b( X_i,j^t-2+X_i,j^t-2  X_i,j^t-3)} .


This is an AR(3) model with parameter vector = (a, b, ξ_1, …, ξ_p,
η_1, …, η_p)^∈⊂ℝ_+^2p+2. 
The probability to form a new edge between nodes i and j at time t is reduced if X_i,j^t-2 =0, and it is reduced further if, in addition, X_i,j^t-3=0. 
The probability to dissolve an existing edge is reduced if
X_i,j^t-2 =1, and it is reduced further if, in addition, X_i,j^t-3=1. Hence if the edge status between two nodes is unchanged for 2 or 3 time periods, the probability for it remaining unchanged next time is larger than that otherwise. 

Model (<ref>) defines an AR(3) network process _t = (X_i,j^t)_p× p with p(p-1)/2 independent edge processes. Although the conclusion on the AR(1) stationarity in the last paragraph of Section <ref> does not apply directly, this AR(3) network process is also strictly stationary, which is implied by the fact that {X_i,j^t}_t≥ 1 is strictly stationary for each 1≤ i<j≤ p. Formally, for given (i,j) such that 1≤ i<j≤ p, let
_t = (X_i,j^t, X_i,j^t-1, X_i,j^t-2)^.
Then {_t}_t≥1 is a homogeneous Markov chain with 2^3 =8 states. Let  denote
the transition probability matrix of {_t}_t≥1. Then 
 is a 8× 8 matrix with only 2 positive elements in each row and each column, provided that ξ_i ξ_j,  η_i η_j ∈ (0, e). It is straightforward to check that each row or column of ^2 has only 4 positive elements, and, more importantly, all the elements of ^3 is positive. Hence, the Markov chain {_t}_t≥1 is irreducible. By Theorems 3.1 and 3.3 in Chapter 3 of <cit.>,  the process {_t}_t≥1 is strictly stationary, and so is {X_i,j^t}_t≥1.


The persistent connectivity or non-connectivity is widely observed in, for example, brain networks, gene connections and social networks. The stability ERGM of <cit.> does not differentiate between the propensity for retaining an existing edge and that for retaining a no-edge status.



 §.§ Transitivity model
 

We propose an AR(1) model to reflect the feature of transitivity which
refers to the phenomenon that nodes are more likely to link if they share
links in common (i.e. ‘the friend of my friend is also my friend’).
 To this end, we specify the transition probabilities as follows:
 

    α_i,j^t-1()    = ξ_i ξ_jexp(a U_i,j^t-1) /1+ exp(a U_i,j^t-1)+ exp(b V_i,j^t-1) ,
    β_i,j^t-1()    =η_i η_j exp(b V_i,j^t-1)/1+ exp(a U_i,j^t-1)+ exp(b V_i,j^t-1) ,

where = (a, b, ξ_1, …, ξ_p, η_1, …, η_p)^⊤∈ℝ_+^2p+2, and

    U_i,j^t-1 =    1 p-2∑_k: k≠ i,j X_i,k^t-1 X_j,k^t-1 ,
     V_i,j^t-1 =    1  p-2∑_k: k≠ i, j{X_i,k^t-1(1-X_j,k^t-1)+ (1-X_i,k^t-1)X_j,k^t-1} .


 The pair (U_i,j^t-1, V_i,j^t-1) characterizes the  number of
nodes with which both nodes i and j are connected,  and the  number of nodes
with which only one of i and j is connected at  time t-1.  The larger U_i,j^t-1 is (i.e. the more common friends i and j share at time t-1), the more likely X_i,j^t = 1. The larger V_i,j^t-1 is, the more likely X_i,j^t = 0. This reflects the transitivity of the networks. High levels of transitivity are found in various networks including friendship networks, industrial supply-chains, international trade flows, and alliances across firms and nations.   Note that the quantity U_i,j^t-1, used in <cit.> to define the edge status of X_i,j^t, reflects
the information based on which companies such as Facebook and LinkedIn have recommended new links to their customers. See also the transitivity ERGM of <cit.>.

We may use different parameters a and b in defining α_i,j^t-1() and β_i,j^t-1() in (<ref>). We do not pursue this more general form as (i) using different ξ_i and η_i  reflects already the differences in the propensity between forming a new edge and dissolving an existing edge, and, perhaps more importantly, (ii) since most practical networks are sparse, the effective sample size for estimating the transition probability from the state of an existing edge is small. Therefore estimating the parameters only occurring in β_i,j^t-1() will be harder than those in α_i,j^t-1(). Using the same a and b in both α_i,j^t-1() and β_i,j^t-1() improves the estimation by pulling the information together. See also the relevant simulation results in Section <ref> in the online supplementary.





 

 
 


  
 



























 






















 
 
 
   
  
 








§ ESTIMATION
 


 §.§ General approach


The natural units of observation in our model are the X^t_i,j, indicating presence or absence of an edge between nodes i and j at time t.  Intuitively, the extent to which these observations can contribute useful information to the estimation of a given element θ_l of ∈Θ depends in turn on the extent to which that element plays a consistent role over time t in the corresponding probabilities

    γ_i,j^t-1() =α_i,j^t-1() +X_i,j^t-1{1 - α_i,j^t-1() -β_i,j^t-1()} .
 
By (<ref>), we have γ_i,j^t-1=γ_i,j^t-1(_0). 

We formalize the above intuition as follows.
 
Write =(θ_1,…,θ_q)^, where q≥ 1 is the total number of parameters.  Let 

    𝒢={l∈[q]:γ_i,j^t-1()  involves θ_l  for  all 1≤ i< j≤ p  and t∈[n]∖[m]} .

We call _𝒢 and _𝒢^,
respectively, the global parameter vector and the local parameter vector. 

In all three models presented in Section <ref>, {ξ_i}_i=1^p and {η_i }_i=1^p are local parameters, while all the other parameters in the models are global parameters. 
As we discuss below in Section <ref>, the global parameter vector _0,𝒢 and the local parameter vector _0,𝒢^ need to be treated differently, which we accomplish via partial likelihoods. The resulting estimators may also entertain different convergence rates.

The asymptotic theory on the convergence rates and the limiting distributions of the proposed estimators will be developed here under the scenario where the sample size n→∞ while the number of nodes p can be either fixed or diverge together with n. When p diverges with n, both the ergodicity and the central limit theorem for stationary Markov chains no longer apply even when _t is stationary  (see the last paragraph in Section <ref>). Based on the conditional independence in our models, regardless of whether p is fixed or diverges with n, we can construct some martingale difference sequences, appropriate partial sums of which are amendable to the required asymptotic analysis without the stationarity assumption.

We develop the estimation theory for our models in three stages below.  Sufficient conditions for identification of _0 are established with respect to an expected partial log-likelihood 
ℓ^(l)_n,p(), defined in (<ref>).  An initial estimator  results from maximizing the corresponding partial log-likelihoods ℓ̂^(l)_n,p() defined in (<ref>), for each l∈[q]. Finally, because of the potential high-dimensionality of our models (number of local parameters increasing with number of nodes), these estimators can suffer from slow rates of convergence.  We offer estimators with improved rate of convergence, derived as a refinement of the initial estimator via the notion of projected score functions.




 §.§ Identification of _0


Let ℱ_t be the σ-field generated by {_1,…,_t}. For any l∈[q], define

    𝒮_l={(i,j):1≤ i< j≤ p and γ_i,j^t-1() involves θ_l for any t∈[n]∖[m]} .

If θ_l is a global parameter, 𝒮_l={(i,j): 1≤ i<j ≤ p}.  
For estimating θ_l for l∈ [q],
put

    ℓ_n,p^(l)()=1/(n-m)|𝒮_l|∑_t=m+1^n∑_  (i,j) ∈𝒮_l𝔼_ℱ_t-1{log[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t]} ,
 




where 𝔼_ℱ_t-1(·) denotes the conditional expectation given ℱ_t-1 with the unknown true  parameter vector _0. 



For any t∈[n]∖[m] and 1≤ i<j≤ p, due to log x≤ x-1 for any x>0, we have

    𝔼_ℱ_t-1{log[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t]}
                   -𝔼_ℱ_t-1{log[{γ_i,j^t-1(_0)}^X_i,j^t{1-γ_i,j^t-1(_0)}^1-X_i,j^t]}
    
            ≤𝔼_ℱ_t-1[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t/{γ_i,j^t-1(_0)}^X_i,j^t{1-γ_i,j^t-1(_0)}^1-X_i,j^t]-1
            =γ_i,j^t-1()/γ_i,j^t-1(_0)·γ_i,j^t-1(_0)+1-γ_i,j^t-1()/1-γ_i,j^t-1(_0)·{1-γ_i,j^t-1(_0)}-1=0 ,

which implies
ℓ_n,p^(l)()≤ℓ_n,p^(l)(_0) for any ∈. Notice that 

    𝔼_ℱ_t-1{log[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t]}
                   -𝔼_ℱ_t-1{log[{γ_i,j^t-1(_0)}^X_i,j^t{1-γ_i,j^t-1(_0)}^1-X_i,j^t]}=0

if and only if 

    {γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t/{γ_i,j^t-1(_0)}^X_i,j^t{1-γ_i,j^t-1(_0)}^1-X_i,j^t≡1 ,

where (<ref>) is equivalent to γ_i,j^t-1()=γ_i,j^t-1(_0). Hence, for any ∈∖{_0}, ℓ_n,p^(l)()=ℓ_n,p^(l)(_0) if and only if γ_i,j^t-1()=γ_i,j^t-1(_0) for any t∈[n]∖[m] and (i,j)∈𝒮_l. To guarantee the identification of _0, we impose the following regularity conditions.


(i) There exists some universal constant C_1>0 such that

    min_t∈[n]∖[m]min_i,j: 1≤ i<j≤ pinf_∈γ_i,j^t-1(){1-γ_i,j^t-1()}≥ C_1 .

(ii) For any 1≤ i<j≤ p and t∈[n]∖[m], γ_i,j^t-1() is thrice continuously
differentiable with respect to ∈. Furthermore, there exists some universal constant
C_2>0 such that

    max_t∈[n]∖[m]max_i,j: 1≤ i<j≤ psup_∈|∂^k γ_i,j^t-1()/∂^k|_∞≤ C_2

for any k∈[3].



Condition <ref> specifies conditions for  the parameter space
. Recall that γ_i,j^t-1() =α_i,j^t-1() +X_i,j^t-1{1 - α_i,j^t-1() -β_i,j^t-1()}. Due to X_i,j^t-1∈{0,1}, Condition <ref>(i) holds if there exist two universal constants c_1,c_2∈(0,1) with c_1<c_2 such that

    c_1≤α_i,j^t-1()≤ c_2   and   1-c_2≤β_i,j^t-1()≤ 1-c_1

for any ∈, t∈[n]∖[m] and 1≤ i<j≤ p. Also, Condition <ref>(ii) holds provided that

    |∂^k α_i,j^t-1()/∂^k|_∞≤ C_2   and   |∂^k β_i,j^t-1()/∂^k|_∞≤ C_2

for any ∈, t∈[n]∖[m] and 1≤ i<j≤ p. Based on the explicit forms of α_i,j^t-1() and β_i,j^t-1() in the specific models, we can identify the associated restrictions for the parameter space .





For any 1≤ i<j≤ p, we define 

    ℐ_i,j={l∈ [q]:γ_i,j^t-1() involves θ_l for any t∈[n]∖[m]} .










There exists a universal constant s≥1 such that
max_1≤ i<j≤ p|ℐ_i,j|≤ s .



Condition <ref> requires that the dynamics of each edge process {X_i,j^t}_t≥ 1 be driven by a finite number of  parameters. Hence the number of global parameters is finite while the total number of local parameters may diverge together with p. For the density-dependent model introduced in Section <ref>, we have _ℐ_i,j  = (a_0,a_1,b_0,b_1,ξ_i,ξ_j, η_i,η_j)^ with s = 8. For both the persistence model and transitivity model introduced in Sections <ref> and <ref>, we have _ℐ_i,j  = (a,b,ξ_i,ξ_j,η_i,η_j)^ with s=6.



There exists a universal constant C_3>0 such that

    min_i,j: 1≤ i<j≤ pλ_min{1/n-m∑_t=m+1^n∂γ_i,j^t-1(_0)/∂_ℐ_i,j∂γ_i,j^t-1(_0)/∂_ℐ_i,j^}≥ C_3





with probability approaching one when n→∞. 




























Let Conditions <ref>–<ref> hold, and C_*=2(2C_1^-2+C_1^-3)C_2^3+3(C_1^-1+C_1^-2)C_2^2+C_1^-1C_2 with (C_1,C_2) specified in Condition <ref>. Assume sup_∈|-_0|_∞<2C_3/(C_*s^3).
As n→∞, it holds 
with probability approaching one that 

    ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()≥C̅/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_2^2

for any ∈ and l∈[q], where C̅>0 is a universal constant.  





































The proof of Proposition <ref> is given in Section <ref> of the supplementary material.
Notice that |𝒮_ l∩𝒮_l'|=|𝒮_ l| for any l' ∈𝒢∪{l}.
By Proposition <ref>, it holds with probability approaching one that for any ∈ and l ∈ [q],

    ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()≥    C̅/|𝒮_l|∑_(i,j) ∈𝒮_l∑_l'∈ℐ_i,j|θ_l'-θ_0,l'|_2^2
    = C̅/|𝒮_l|∑_l'=1^q∑_(i,j)∈𝒮_l∩𝒮_l'|θ_l'-θ_0,l'|^2
    
    =    C̅∑_l'∈𝒢∪{l}|θ_l'-θ_0,l'|^2+C̅∑_l'∈𝒢^∖{l}|𝒮_l∩𝒮_l'||θ_l'-θ_0,l'|^2/|𝒮_l| .



 Hence, for any l∈[q], the function  ℓ_n,p^(l)(·) defined as (<ref>) is a good candidate for identifying θ_0,l and the global parameter vector _0,𝒢 but is powerless in identifying θ_0,l' with l'∈𝒢^∖{l} if |𝒮_l∩𝒮_l'|≪|𝒮_l|.





 §.§ Initial estimation for _0
 

With available observations _1, …, _n, since {_t}_t≥1 is a Markov chain with order m, the likelihood function for , conditionally on
_1, …, _m, admits the form

    ℒ_n,p(_n, …, _m+1 | _m, …, _1; ) =
    ∏_t=m+1^n L_t,p(_t  | _t-1, …, _t-m;) ,

where L_t,p(_t  | _t-1, …, _t-m;) is the transition probability of _t given _t-1,…,_t-m. 
By (<ref>), the (normalized) log-likelihood admits the form

    2/(n-m)p(p-1)logℒ_n,p(_n, …, _m+1 | _m, …, _1; )
                =2/(n-m)p(p-1)∑_t=m+1^n∑_i,j: 1≤ i<j≤ plog[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t] ,

which is the sample version of ℓ_n,p^(l)() defined as (<ref>) with l∈𝒢.
As pointed out below (<ref>), we should not estimate the local parameters based on this full log-likelihood. Therefore, for each l∈[q], we define 



















 










    ℓ̂_n,p^(l)() =1/(n-m)|𝒮_l|∑_t=m+1^n∑_ (i,j) ∈𝒮_llog[{γ_i,j^t-1()}^X_i,j^t{1-γ_i,j^t-1()}^1-X_i,j^t] ,

which contains only the terms depending on θ_l on the right-hand side of (<ref>) (with a rescaled normalized constant).

For any l∈[q], Lemma <ref> in the supplementary material shows that ℓ̂_n,p^(l)() converges in probability to ℓ_n,p^(l)() defined as (<ref>) uniformly over ∈. 
Together with Proposition <ref>, we can estimate the global parameter vector _0,𝒢 by maximizing the full log-likelihood ℓ̂_n,p^(l')() with some l'∈𝒢, and estimate the local parameter θ_0,l with l∈𝒢^ by maximizing the corresponding (partial) log-likelihood ℓ̂_n,p^(l)(). More specifically, letting (θ̂_*,1^(l), …,θ̂_*,q^(l) )^ =max_∈ℓ̂_n,p^(l)() for each l∈[q], we define the initial estimator =(_𝒢^,_𝒢^^)^ for _0 as

    _𝒢=(θ̂_*,l^(l'))_l∈𝒢   and   _𝒢^=(θ̂_*,l^(l))_l∈𝒢^

for some l'∈𝒢. Due to 𝒮_l={(i,j):1≤ i < j ≤ p}
for any ℓ∈𝒢, we know ℓ̂_n,p^(l_1)()=ℓ̂_n,p^(l_2)() for any l_1,l_2∈𝒢, which implies that the estimator _𝒢 given in (<ref>) does not depend on the selection of l'∈𝒢. 





To investigate the theoretical properties of the  estimator =(_𝒢^,_𝒢^^)^, we define 
    {
       c_n,𝒢^2=    qlog(np)/√(n)p+q^3/2log^3/2(np)/√(n)p^2 , 
    
       c_n,𝒢^^2=    qlog(nS_𝒢^,min)/√(nS_𝒢^,min)+q^3/2log^3/2(nS_𝒢^,min)/√(n)S_𝒢^,min ,
    .

where S_𝒢^,min=min_l∈𝒢^|𝒮_l|. 
Theorem <ref> shows that the convergence rate of the
initial estimator for
the local parameters is slower than that of the global parameters if
S_𝒢^,min≪ p^2. The proof of Theorem <ref> is given in Section <ref> of the supplementary material.









 
















Let the conditions of Proposition <ref> hold. Then |_𝒢-_0,𝒢|_2=O_ p(c_n,𝒢) and |_𝒢^-_0,𝒢^|_∞=O_ p(c_n,𝒢^). 


 
By Theorem <ref>,  the initial estimator _𝒢 for the global parameters is consistent provided that 

    q≪min{√(n)p/log(np), n^1/3p^4/3/log(np)} ,

and  the initial estimator _𝒢^ for the local parameters is consistent provided that

    q≪min{√(nS_𝒢^,min)/log(nS_𝒢^,min), n^1/3S_𝒢^,min^2/3/log(nS_𝒢^,min)} .

For the density-dependent model introduced in Section <ref>, we have q=2p+4 and S_𝒢^,min=p-1. For both the persistence model and transitivity model introduced in Sections <ref> and <ref>, we have q=2p+2 and S_𝒢^,min=p-1. Hence, for these three models, Theorem <ref> gives the convergence rates of _𝒢 and _𝒢^ as follows: 

    |_𝒢-_0,𝒢|_2=     O_ p{log^1/2(np)/n^1/4∨log^3/4(np)/(np)^1/4} ,  
    
    |_𝒢^-_0,𝒢^|_∞=     O_ p{p^1/4log^1/2(np)/n^1/4∨p^1/4log^3/4(np)/n^1/4} ,

which implies the consistency of _𝒢 provided that log p ≪ n^1/2, and the consistency of _𝒢^ provided that p≪ n(log n)^-3. 









 §.§ Improved estimation for _0
 

Recall =(θ_1,…,θ_q)^. The initial estimator  specified in (<ref>) suffers from slow convergence rates due to the high dimensionality of .
In this section, we improve the estimation for each component θ_0,l by projecting the score function onto certain direction. See (<ref>) below for details. An improved estimator for θ_0,l is then obtained by solving the projected score function while  letting _-l =
_-l. The projection mitigates the impact 
of _-l in the improved estimation for θ_0,l. This strategy was initially proposed by <cit.> and <cit.> for constructing the valid confidence regions of some low-dimensional subvector of the whole parameters in high-dimensional models with removing the impact of the high-dimensional nuisance parameter.



For (c_n,𝒢,c_n,𝒢^) defined as (<ref>), put

    Δ_n=max{|𝒢|c_n,𝒢^2,|𝒢^|^2c_n,𝒢^^2} .

For any t∈[n]∖[m], l∈[q] and ∈, we define

    _t^(l)()=1/|𝒮_l|∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂γ_i,j^t-1()/∂ .

Then the score function can be written as

    ∂ℓ̂_n,p^(l)()/∂=1/n-m∑_t=m+1^n^(l)_t() .

To estimate θ_0,l, _-l 
can be treated as a nuisance parameter vector. Following  <cit.> and <cit.>, 
we  project _t^(l)() to form a new estimating function:

    f̂_t^(l)()=_l^_t^(l)(),

where _l is defined as

    _l=min_∈ℝ^q||_1    s.t.  |{1/n-m∑_t=m+1^n∂_t^(l)()/∂}^-_l|_∞≤τ .

In the above expression, τ>0 is a tuning parameter satisfying τ≲Δ_n^1/2 with Δ_n defined as (<ref>), =(θ̃_1,…,θ̃_q)^ is the initial estimator defined as (<ref>), and _l is a q-dimensional vector with the l-th component being 1 and other components being 0. Then we can re-estimate _0 by =(θ̌_1,…,θ̌_q)^, where 

    θ̌_l=min_θ_l∈ B(θ̃_l, r̃)|1/n-m∑_t=m+1^nf̂_t^(l)(θ_l,_-l)|^2

for some r̃>0 satisfying max{c_n,𝒢,c_n,𝒢^}≪r̃≪ 1. 

To construct the convergence rate of |-_0|_∞, we need the following regularity condition, which is analogous to Condition 1 of <cit.> and Condition 7 of <cit.>. See the discussion there for the validity of such condition.

  For each l∈[q], there is a nonrandom vector _l∈ℝ^q 
  such that |_l|_1≤ C_4 for some universal constant C_4>0, and max_l∈[q]|_l-_l|_1=O_ p(ω_n) for some ω_n→0 satisfying ω_n(log q)^1/2log(qn)=o(1).   


 Proposition <ref> shows that  has faster convergence rate than the initial estimator  given in (<ref>).  The proof of Proposition <ref>
is given in Section <ref> of the supplementary material.

Let the conditions of Proposition <ref> and Condition <ref> holds. Then
|-_0|_∞=O_ p(Δ_n), where Δ_n is defined as (<ref>).



Based on the obtained , we consider the final estimate =(θ̂_1,…,θ̂_q)^ for _0 defined as follows:

    θ̂_l=min_θ_l∈ B(θ̌_l,ř)|1/n-m∑_t=m+1^nf̂_t^(l)(θ_l,_-l)|^2

for some ř>0 satisfying qΔ_n≪ř≪ 1 with Δ_n defined as (<ref>).

 
Given the initial estimate , there are three tuning parameters (τ,r̃,ř) for deriving our final estimate . 
For the density-dependent model introduced in Section <ref>, we have |𝒢| = 4 and |𝒢^| = 2p. For both the persistence model and transitivity model introduced in Sections <ref> and <ref>, we have |𝒢| = 2 and |𝒢^| = 2p.  Together with Remark <ref>, we have Δ_n = n^-1/2p^5/2log^3/2(np) for these three models. The improved estimation procedure thus requires τ≲ n^-1/4p^5/4log^3/4(np), n^-1/4p^1/4log^3/4(np)≪r̃≪ 1 and n^-1/2p^7/2log^3/2(np)≪ř≪ 1, which suggests p≪ n^1/7(log n)^-3/7.
In practice, for the three models introduced in Section <ref>, we compute the final estimate  with τ proportional to n^-1/4p^5/4log^3/4(np) and adopting reasonably large r̃ and ř. Numerical experiments in Sections <ref> and <ref>  validate the robustness of our proposed estimation procedure regarding the selections of r̃ and ř as long as _0 falls within the defined search range.






For any ∈ and l∈[q], define

    ζ_n,l()=1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l1/γ_i,j^t-1(){1-γ_i,j^t-1()}{_l^∂γ_i,j^t-1()/∂}^2 ,

where _l is given in Condition <ref>. Under Conditions <ref> and <ref>, we have |ζ_n,l()|≤ C_1^-1C_2^2C_4^2, which implies that, for any ∈, ζ_n,l() is a bounded random variable. To construct the asymptotic distribution of each θ̂_l, we require the following condition. 


For each l∈[q], there exists some random variable κ_l≥0 such that
ζ_n,l(_0)→κ_l in probability as n→∞.



For each l∈ [q] and t≥ m+1,
let 
    v_l^t-1=1/|𝒮_l|∑_(i,j)∈𝒮_l1/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}{_l^∂γ_i,j^t-1(_0)/∂}^2 .

As {ζ_n,l(_0)}_n≥ m+1 is a bounded sequence of random variables for each l∈[q], Condition <ref> is mild and κ_l is a random variable in general. Generally speaking, the asymptotic distribution
of θ̂_l is a mixture of normal distributions. See Theorem <ref> below for details. However, if the long-run variance of {v_l^t-1}_l=m+1^n satisfies the condition

    Var(1/√(n-m)∑_t=m+1^nv_l^t-1)=o(√(n)) ,

κ_l is reduced to a constant 

    κ_l=lim_n→∞𝔼(1/n-m∑_t=m+1^nv_l^t-1) .

Then Theorem <ref> implies that θ̂_l is asymptotically normal distributed. When the sequence
{v_l^t}_t≥ m is α-mixing with the mixing coefficients 
attaining certain convergence rates, (<ref>) holds automatically. 



















Let the conditions of Proposition <ref> and Conditions <ref> and <ref> hold. For each l∈[q], if √(n|𝒮_l|)max{qΔ_n^3/2,q^2Δ_n^2}=o(1) with Δ_n defined as (<ref>), it then holds that

    √(n|𝒮_l|)(θ̂_l-θ_0,l)→√(κ_l)· Z

in distribution as n→∞, where Z is a standard normally distributed random variable independent of κ_l specified in Condition <ref>.


The proof of Theorem <ref>
is given in Section <ref> of the supplementary material.


    
    
(i) Theorem <ref> shows that, for the global parameter θ_l with l ∈𝒢,
    
    |θ̂_l-θ_0,l| = O_ p(1/√(n)p) ,
   
provided that  
 
    q≪   min{n^1/10p^1/5/|𝒢|^3/5log^3/5(np), n^1/13p^8/13/|𝒢|^6/13log^9/13(np),
    
                n^1/10S_𝒢^,min^3/10/|𝒢^|^6/5p^2/5log^3/5(nS_𝒢^,min), n^1/13S_𝒢^,min^6/13/|𝒢^|^12/13p^4/13log^9/13(nS_𝒢^,min)} ,

 and for the local parameter θ_l with l ∈𝒢^,
        
    |θ̂_l-θ_0,l| = O_ p(1/√(n |𝒮_l|)) ,

  provided that 





  
    q≪   min{n^1/10p^3/5/|𝒢|^3/5|𝒮_l|^1/5log^3/5(np), n^1/13p^12/13/|𝒢|^6/13|𝒮_l|^2/13log^9/13(np), n^1/8p^1/2/|𝒢|^1/2|𝒮_l|^1/8log^1/2(np),
    
                n^1/10S_𝒢^,min^3/10/|𝒢^|^6/5|𝒮_l|^1/5log^3/5(nS_𝒢^,min), n^1/13S_𝒢^,min^6/13/|𝒢^|^12/13|𝒮_l|^2/13log^9/13(nS_𝒢^,min)} .

In particular, for the three models introduced in Section <ref>,  the estimators satisfy |θ̂_l-θ_0,l| = O_ p(n^-1/2p^-1) for l ∈𝒢 if p ≪ n^1/23 (log n)^-9/23, and |θ̂_l-θ_0,l| = O_ p{(np)^-1/2} for l ∈𝒢^ if p ≪ n^1/21 (log n)^-3/7. Compared with the results in Theorem <ref>, the improved estimator  achieves a faster convergence rate than the initial estimator . 

(ii) For each l∈[q], write 

    ζ̂_n,l()=1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l1/γ_i,j^t-1(){1-γ_i,j^t-1()}{_l^∂γ_i,j^t-1()/∂}^2

with _l defined as (<ref>). 
Since ζ̂_n,l()-ζ_n,l(_0)→0 in probability as n→∞, by Corollary 3.2 of <cit.>, it holds that 

    √(n|𝒮_l|/ζ̂_n,l())(θ̂_l-θ_0,l)→𝒩(0,1)

in distribution as n→∞, provided that ℙ(κ_l>0)=1. We can use (<ref>) to construct the confidence interval for each θ_l. 
    













§ SIMULATION WITH TRANSITIVITY MODELS
 

In this section, we use the transitivity model introduced in Section <ref> as an example to illustrate numerical behaviour of both the initial estimation proposed in Section <ref>
and the improved estimation suggested in Section <ref>. 



 §.§ Implementation details



Network data {_1, …, _n} used in the experiments described below are generated according to (<ref>), (<ref>) and (<ref>). 
For each sample, we generate a sequence of  length n+200, and discard the first 200 observations.

Regarding implementation of our estimation procedures, recall  that 𝒢 and 𝒢^ are, respectively, the index sets of the global parameters and the local parameters. For the transitivity model (<ref>), we have =(a,b,ξ_1,…,ξ_p,η_1,…,η_p)^, where a and b are the global parameters and {ξ_i}_i=1^p and {η_i}_i=1^p are the local parameters.  Hence, for this model we have |𝒮_l| = p(p-1)/2  and |𝒮_l∩𝒮_l'| = p-1 when  l∈𝒢 and l'∈𝒢^. By (<ref>), for each given l∈𝒢 and ∈, it holds with probability approaching one that ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()≥C̅|_𝒢-_0,𝒢|_2^2+2C̅p^-1|_𝒢^-_0,𝒢^|_2^2 for some universal constant C̅>0 independent of , which means that the function  ℓ_n,p^(l)(·) defined as (<ref>) exhibits robustness against fluctuations in the values of local parameters when l∈𝒢 and p is large. 

Motivated by this fact, when we compute the initial estimator _𝒢  for the global parameter vector _0,𝒢, we can just approximate _𝒢 by _𝒢^( app)=max__𝒢ℓ̂_n,p^(l)(_𝒢,_𝒢^), for some given _𝒢^ and ℓ̂_n,p^(l)() defined as (<ref>) for some l∈𝒢. This simple idea can significantly improve the computational efficiency. Specifically, note that computing the original _𝒢 requires solving an optimization problem with 2p+2 variables while this alternative approach only requires solving an optimization problem with two variables. Our above discussion guarantees _𝒢^( app) can approximate _𝒢 well. 
Similarly, when we compute the initial estimator θ̃_l for the local parameter θ_0,l with l∈𝒢^, 
we can  approximate it by θ̃_l^( app)=max_θ_lℓ̂_n,p^(l)(_𝒢^( app),  θ_l,  _𝒢^∖{l}) with some given _𝒢^∖{l}.

In practice, we first estimate the global parameters a and b via the Quasi-Newton method, 
 given certain initial values for the local parameters {ξ_i}_i=1^p and {η_i}_i=1^p. To be specific, we consider 9 different sets of the initial values between 0.5 and 0.9 for {ξ_i}_i=1^p and {η_i}_i=1^p, and compute ã^(ν) and b̃^(ν)  for the ν-th initial setting.






With a=ã^(ν) and b=b̃^(ν), we then compute ξ̃_1^(ν), …, ξ̃_p^(ν), η̃_1^(ν), …, η̃_p^(ν) through maximizing each of the associated ℓ̂_n,p^(l)() with l∈𝒢^. 

Subsequently, the improved estimates â^(ν), b̂^(ν), ξ̂_1^(ν), …, ξ̂_p^(ν), η̂_1^(ν), …, η̂_p^(ν) are obtained according to (<ref>)
with (τ, r̃, ř) = (0.5 Δ_n^1/2,0.5,0.1) for the local parameters and (τ, r̃, ř) =  (0.01 Δ_n^1/2, 10, 2) for the global parameters. 
The simulations in this section utilise our development R package , which provides a user-friendly implementation of the practical estimation algorithm described above.




 §.§ Estimation errors


Here we report results on experiments exploring the behavior of the
initial estimator  given in (<ref>) and the improved estimator  given in (<ref>).  For simplicity, we again set all the true values for  {ξ_i}_i=1^p to be the same, and those for {η_i}_i=1^p also to be the same. The same four sets of parameter values were used as in Section <ref>.
We set n∈{100, 200} and p∈{ 50, 100, 150}.
For each setting, we replicate the estimation 400 times.

Table <ref> presents the means  and the standard errors, over the 400 replications, of the relative mean absolute errors (rMAE):


    rMAE(ξ̂_i)= 1/9∑_ν = 1^91/p∑_i =1^p|ξ̂_i^(ν) - ξ_i/ξ_i|   and  rMAE(â)= 1/9∑_ν = 1^9|â^(ν) - a/a|,

where the sum over ν corresponds to taking the average over the 9 initial values discussed in Section <ref>.


As indicated in the table, the improved estimator (<ref>) is significantly more accurate than the initial estimator (<ref>). For example, for setting (0.6,0.7,15,10), we observe an approximate 70% improvement in the estimation accuracy of a when p = 150. Furthermore, the setting (0.7,0.8,30,15) attains the lowest overall estimation errors. 
This is well-expected, as this is the most dynamic setting among the four settings considered. 










§ REAL DATA ANALYSIS: EMAIL INTERACTIONS
 


In this section, we apply the transitivity model (<ref>) to a dynamic network dataset of email interactions in a medium-sized Polish manufacturing company, from January to September 2010 <cit.>.
We analyze a subset of the data among p=106 of the most active participants out of an original 167 employees.
The organizational tree of direct reports in the company is also available for these employees.
Each of the n=39 network snapshots corresponds to a non-overlapping time window, with X_i,j^t=1 if participants i and j exchanged at least one email in the previous seven days. This accounts for periodic weekly effects.

We first present some preliminary summaries of the data to inspect the stationarity of the network and the effective sample size. 
The behavior shown in Figure <ref> suggests a change point in the network behavior, in terms of both edge density D_t and two dynamics density measures D_1,t and D_0,t (see (<ref>)). 
Hence in the following analysis, we fit the model separately to the first 13 and last 26 snapshots, referred to as “period 1” and “period 2”. 
In the right panel, about 4% of node pairs see a grown edge or a dissolved edge between consecutive snapshots. 
After accounting for the low edge density, the relative frequency of growing a new edge is about 5%, while relative frequency of an existing edge to dissolve is only about 45%, clear evidence of temporal edge dependence.



We also identify empirical evidence in the data for transitivity effects. This is demonstrated in Figure <ref>.
 To construct these plots, we partition the edge variables as follows: for each integer ℓ≥0, define

    𝒰_ℓ = {(i,j,t) : 1 ≤ i < j ≤ p ,  t ∈ [n]∖{1} ,  X_i,j^t-1=0 ,  U_i,j^t-1 = ℓ/(p-2) } , 
       𝒱_ℓ = {(i,j,t) : 1 ≤ i < j ≤ p ,  t ∈[n]∖{1} ,  X_i,j^t-1=1 ,  V_i,j^t-1 = ℓ/(p-2) } ,
       𝒰_ℓ^1 = {(i,j,t) ∈𝒰_ℓ,   X_i,j^t=1 },   𝒱_ℓ^0 = {(i,j,t) ∈𝒱_ℓ,  X_i,j^t=0
        } ,

where U_i,j^t-1 and V_i,j^t-1 are given in (<ref>). 

The left panel plots the relative frequency |𝒰_ℓ^1|/|𝒰_ℓ| against ℓ for ℓ =0,1, ⋯, showing that this frequency of grown edges tends to be higher for node pairs with more common neighbours in the previous snapshot. 
The right panel analogously plots the relative frequency  |𝒱_ℓ^0|//|𝒱_ℓ| against ℓ, and shows a similar increasing relationship between disjoint neighbours and frequency of dissolved edges.



This is confirmed by the fit of our model parameters, using the estimation algorithm described in Section <ref>, and implemented in our development R package . For period 1, we estimate the global parameters â=13.12 and b̂=9.34, suggesting a tendency towards edge growth given more common neighbours, and edge dissolution given more distinct neighbours, which agrees with the empirical evidence in Figure <ref>. 
We interpret the estimates of the local parameters {ξ_i}_i=1^106 and {η_i}_i=1^106 in the left panel of Figure <ref>.
The estimates {ξ̂_i}_i=1^106 have mean 0.61 and skew towards the right, implying degree heterogeneity in the edge growth. Conversely, the estimates {η̂_i}_i=1^106 have mean 0.89 and skew towards the left. 
There is a decreasing relationship between the paired parameters: employees who tend to grow new edges also tend to maintain existing edges.
Finally, there is an observed relationship between email behavior and company hierarchy: managers (non-leaf nodes in the organizational tree) tend to have larger estimates {ξ̂_i}_i=1^106 compared to non-managers (means 0.74 and 0.57 respectively), implying that managers are more likely to grow edges. 
However, this increasing pattern does not continue at higher levels of the organizational tree.



The model fit to period 2 shows many of the same patterns. We estimate â=21.69 and b̂=9.84 and summarize the estimates {ξ̂_i}_i=1^106 and {η̂_i}_i=1^106 in the right  panel of Figure <ref>.
Relative to period 1, the larger estimate of a implies a stronger transitivity effect in this time period.
The estimates {ξ̂_i}_i=1^106 now have mean 0.49 and the estimates {η̂_i}_i=1^106 have mean 0.92, to model overall lower edge density.
The decreasing relationship between the paired parameters is stronger, and the means of ξ̂_i for managers and
non-managers are, respectively, 0.68 and 0.43.
Along with the stronger transitivity effect, we interpret that the decreased network density in period 2 has led to a concentration of email activity among a smaller group of employees, many of them managers.

We compare our model to some competing models from the literature in terms of Akaike and Bayesian information criteria (AIC, BIC). 
To briefly describe these competitors: the “global AR model” and “edgewise AR model” fit the model of <cit.>, with either two global switching parameters or two parameters for each edge.  The “edgewise mean model” assumes 

    X_i,j^t iid∼Bernoulli(P_i,j)

with no temporal dependence, and estimates the edge probability {P_i,j}_i,j: i<j for each node pair by its relative frequency in the training set; and the “degree parameter mean model” 
assumes 

    X_i,j^t iid∼Bernoulli(ν_iν_j)

and estimates the degree parameters {ν_i}_i=1^106 by fitting 1-dimensional adjacency spectral embedding <cit.> to the mean adjacency matrix over the training set. 
Note that the edgewise mean model has O(p^2) parameters, while the degree parameter model has O(p) parameters like our AR network model with transitivity.
All of these models can be directly compared using the AR network model likelihood, although only our AR network model with transitivity incorporates edge dependence, and the final two models do not incorporate any temporal dependence. 
Results for both regimes are reported in Table <ref>.



In period 1, our AR network model with transitivity achieves the lowest AIC and BIC, while in period 2 it is outperformed slightly by the edgewise mean model in terms of AIC, but achieves a lower BIC as it uses fewer parameters. 
This reduction of the parameter space is important for modeling sparse dynamic network data: although there is clear temporal edge dependence in this data, the edgewise mean model outperforms the edgewise AR model, as there is low effective sample size to estimate the edge dissolution parameters. 

Finally, we compare the performances of those  models in an edge forecasting task on the final 26 network snapshots (period 2). 
For n_train=10, …, 23, we train these models on the first n_train snapshots of period 2, then forecast the state of each edge n_step steps forward, for n_step=1, 2, 3. 
The combined results are presented in Figure <ref> as receiver operating characteristic (ROC) curves.
We also include a single point summarizing the naive forecasting performance using the most recent observation of that edge in the training set.



Our AR network model with transitivity dominates or is competitive with all the models besides these highly parameterized edgewise models for all choices of n_step. 
The good performance of the edgewise models suggests the presence of higher order structure in this network that cannot be modeled with only two parameters per node.
Note that the edgewise mean and edgewise AR models give very similar, but not identical edge predictions; as mentioned above, due to network sparsity the edgewise AR model has a low effective sample size to estimate the dissolution parameters, leading to slightly worse link prediction performance.




















 
 
 
 
 
 
 
 
 

























































































































































1.1
dcu


1.4

			




	Jinyuan Chang, Qin Fang, Eric D. Kolaczyk, Peter W. MacDonald, and Qiwei Yao



	
	
	

 








	





This supplementary material contains a detailed analysis on the relationship
between the proposed AR models and temporal ERGMs (Section <ref>), all the technical proofs (Section <ref>), additional numerical simulation results (Section <ref>), and the analysis of an additional dynamic network dataset (Section <ref>).




§ RELATIONSHIP TO TEMPORAL ERGMS



A dynamic network sequence follows a temporal ERGM of order m if it satisfies

    ℙ( X_t  |  X_t-1, …,  X_t-m ; ) ∝exp{()^( X_t,  X_t-1, …,  X_t-m)} ,

where : ℝ^q →ℝ^p maps the parameter vector  to the vector of natural parameters, and  maps the data, including the past network snapshots, to the corresponding sufficient statistics.

As in Equation (2) of <cit.>, suppose  factors over the edges of the present snapshot,

    ( X_t,  X_t-1, …,  X_t-m) = ∑_i,j: i < j_i,j(X^t_i,j ;  X_t-1, …,  X_t-m) .

Then 

    ℙ( X_t  |  X_t-1, …,  X_t-m ; ) 
                       ∝∏_i,j: i < jexp{()^_i,j(X^t_i,j ; X_i,j^t-1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )} ,

which implies X_t will have mutually independent edges conditional on the past snapshots. We refer to this as the edge conditional independence assumption, which is a property of AR network models defined in Definition 1 of the main document.

We will show that any edge conditionally independent temporal ERGM can be rewritten as an AR network model.






Denote the logit function by σ(x) = log{x/(1-x)}, and specify an AR network model defined in Definition 1 by setting

    α_i,j^t-1
      =    σ^-1[ ()^{_i,j(1 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                                - _i,j(0 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}] 
    
      :=    σ^-1[ ()^(_i,j,10^t-1 - _i,j,00^t-1)] ,
    β_i,j^t-1
      =    σ^-1[ ()^{_i,j(0 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m)
                                - _i,j(1 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}]
    
      :=    σ^-1[ ()^(_i,j,01^t-1 - _i,j,11^t-1)] .

With renormalizing, we have

    α_i,j^t-1   = exp{()^_i,j,10^t-1}/exp{()^_i,j,10^t-1} + exp{()^_i,j,00^t-1} ,
    
      1 - α_i,j^t-1   = exp{()^_i,j,00^t-1}/exp{()^_i,j,10^t-1} + exp{()^_i,j,00^t-1} , 
    β_i,j^t-1   = exp{()^_i,j,01^t-1}/exp{()^_i,j,01^t-1} + exp{()^_i,j,11^t-1} , 
    
      1 - β_i,j^t-1   = exp{()^_i,j,11^t-1}/exp{()^_i,j,01^t-1} + exp{()^_i,j,11^t-1} .

For the AR network model with (α_i,j^t-1,β_i,j^t-1) specified in (<ref>) and (<ref>), it holds that

    ℙ( X_t |  X_t-1, …,  X_t-m ; )
      = ∏_i,j: i < jℙ(X_i,j^t  |  X_t-1, …,  X_t-m ; ) 
    
              = ∏_i,j: i < j,  X_i,j^t=0,   X_i,j^t-1=0 (1 - α_i,j^t-1)·∏_i,j: i < j,  X_i,j^t=1,  X_i,j^t-1=0α_i,j^t-1
                      ·∏_i,j: i < j, X_i,j^t=0,   X_i,j^t-1=1β_i,j^t-1·∏_i,j: i < j,  X_i,j^t=1,  X_i,j^t-1=1 (1- β_i,j^t-1) 
             ∝∏_i,j: i < j,  X_i,j^t=0,   X_i,j^t-1=0exp{()^_i,j,00^t-1}·∏_i,j: i < j,  X_i,j^t=1,  X_i,j^t-1=0exp{()^_i,j,10^t-1}
                      ·∏_i,j: i < j,  X_i,j^t=0,  X_i,j^t-1=1exp{()^_i,j,01^t-1}·∏_i,j: i < j,  X_i,j^t=1,  X_i,j^t-1=1exp{()^_i,j,11^t-1}
            = ∏_i,j: i < jexp{()^_i,j(X^t_i,j ; X_i,j^t-1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )} ,

which shares the same form as (<ref>).
Thus for any edge conditionally independent temporal ERGM, we can specify an AR network model with the same distribution.

Conversely, suppose we have specified an AR network model such that

    σ(α_i,j^t-1)    = ()^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )  ,
    σ(β_i,j^t-1)    = ()^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m ) ,

for some functions of the  parameter vector , and the past network behavior. We claim that this AR network model can be written as an edge conditionally independent temporal ERGM (<ref>) with () = (()^, ()^)^ and ( X_t ,  X_t-1, …,  X_t-m) = (_α^ , _β^)^, where
_α = ∑_i,j: i < j_α,i,j and _β = ∑_i,j: i < j_β,i,j with

    _α,i,j(X_i,j^t ; X_i,j^t-1,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m)
                = _i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t=1, X_i,j^t-1=0)  ,
       _β,i,j(X_i,j^t ; X_i,j^t-1,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                =_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t=0, X_i,j^t-1=1) .





<cit.> define the concept of separability of a dynamic model for binary networks.


Define two subnetworks X^+_t and X^-_t, where

    X_i,j^t,+   = 1 - I(X_i,j^t=0,X_i,j^t-1=0), 
    
        X_i,j^t,-   = I(X_i,j^t=1,X_i,j^t-1=1),

for all i,j.
A dynamic network model is said to be separable if X^+_t and X^-_t are independent conditional on the past, and do not share any parameters.
In particular, a separable temporal ERGM (STERGM) can be specified by a product of a formation model and a dissolution model:

    ℙ( X_t  |  X_t-1, …,  X_t-m ; ) 
             ∝exp{^+(^+)^^+( X^+_t,  X_t-1, …,  X_t-m) + ^-(^-)^^-( X_t^-,  X_t-1, …,  X_t-m)} ,

where ^+ and ^- map parameter vectors ^+ and ^- to the vectors of natural parameters, and ^+ and ^- map the data, including the past network snapshots, to the corresponding sufficient statistics.













Under the edge conditional independence assumption, we can write

    ^+( X^+_t,  X_t-1, …,  X_t-m) = ∑_i < j_i,j^+(X_i,j^t,+; X_i,j^t-1,  X_t-1∖ X_i,j^t-1, …,  X_t-m),

or equivalently

    ^+( X^+_t,  X_t-1, …,  X_t-m) = ∑_i < j_i,j^+(X_i,j^t; X_i,j^t-1,  X_t-1∖ X_i,j^t-1, …,  X_t-m),

since for all i < j, X_i,j^t,+ can be recovered from X_i,j^t and X_i,j^t-1.
Define ^-_i,j analogously for the dissolution model.

In this way X_t is an edge conditionally independent TERGM with parameter vector (^+,^-), natural parameter

    (ς^+(^+),ς^-(^-)),

and sufficient statistic

    ( ∑_i < j^+_i,j, ∑_i < j^-_i,j).





Following the above construction to rewrite this as an AR network model, we can write

    α_i,j^t-1
      =    σ^-1[ ^+(^+)^{^+_i,j(1 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                       - ^+_i,j(0 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}
          +^-(^-)^{^-_i,j(1 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                       - ^-_i,j(0 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}],
 
and note that 

    ^-_i,j(1 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m)    = ^-_i,j(0 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m ) 
       = ^-_i,j(0 ; 0 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )

for all i < j and t, since X_i,j^t,-=0 whenever X_i,j^t-1=0.
Thus α_i,j^t-1 is free of ^-. Similarly,
we can write 

    β_i,j^t-1
      =    σ^-1[ ^+(^+)^{^+_i,j(0 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                       - ^+_i,j(1 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}
          +^-(^-)^{^-_i,j(0 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                       - ^-_i,j(1 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )}],
 
and

    ^+_i,j(0 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m)    = ^+_i,j(1 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m ) 
       = ^+_i,j(1 ; 1 ,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )

for all i < j and t, since X_i,j^t,+=1 whenever X_i,j^t-1=1.
Then β_i,j^t-1 is free of ^+.
Hence any edge conditionally independent STERGM can be written as an AR network model with separable parameters.

Conversely, suppose we have specified an AR network model such that

    σ(α_i,j^t-1)    = (_α)^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m ) 
    σ(β_i,j^t-1)    = (_β)^_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )

with separable parameters _α and _β.

We follow the same construction as above to rewrite this model as an edge conditionally independent TERGM, with parameter vector () = ((_α)^, (_β)^)^ and sufficient statistics ( X_t ,  X_t-1, …,  X_t-m) = (_α^ , _β^)^, where
_α = ∑_i,j: i < j_α,i,j and _β = ∑_i,j: i < j_β,i,j with

    _α,i,j(X_i,j^t ; X_i,j^t-1,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                = _i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t=1, X_i,j^t-1=0)  , 
                = _i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t,+=1, X_i,j^t-1=0)  , 
       _β,i,j(X_i,j^t ; X_i,j^t-1,  X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m) 
                =_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t=0, X_i,j^t-1=1) . 
                =_i,j( X_t-1∖ X_i,j^t-1,  X_t-2, …,  X_t-m )I(X_i,j^t,-=0, X_i,j^t-1=1) .


Note that in (<ref>) and (<ref>), the sufficient statistics depend on X_i,j^t only through X_i,j^t,+ and X_i,j^t,- respectively. 
It follows that when written as a conditionally independent TERGM, the distribution factors into a product of a formation model and dissolution model, as in (<ref>). 
Thus this AR network model is an edge conditionally independent STERGM.








§ TECHNICAL PROOFS
 



 §.§ Proof of Proposition <ref>
 

Recall 

    ℓ_n,p^(l)()=1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_lh_i,j^t-1()

with

    h_i,j^t-1()
    =log{1- γ_i,j^t-1()}+γ_i,j^t-1(_0)log{γ_i,j^t-1()/1- γ_i,j^t-1()} .

Then

    ∂ h_i,j^t-1()/∂θ_l_1=    γ_i,j^t-1(_0)-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂γ_i,j^t-1()/∂θ_l_1 ,
    ∂^2 h_i,j^t-1()/∂θ_l_1∂θ_l_2=   -[1/γ_i,j^t-1(){1-γ_i,j^t-1()}+{γ_i,j^t-1(_0)-γ_i,j^t-1()}{1-2γ_i,j^t-1()}/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2]
            ×∂γ_i,j^t-1()/∂θ_l_1∂γ_i,j^t-1()/∂θ_l_2
       +γ_i,j^t-1(_0)-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂^2γ_i,j^t-1()/∂θ_l_1∂θ_l_2 ,
    ∂^3h_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3=    2[1-2γ_i,j^t-1()/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2+γ_i,j^t-1(_0)-γ_i,j^t-1()/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2
            +{γ_i,j^t-1(_0)-γ_i,j^t-1()}{1-2γ_i,j^t-1()}^2/{γ_i,j^t-1()}^3{1-γ_i,j^t-1()}^3]∂γ_i,j^t-1()/∂θ_l_1∂γ_i,j^t-1()/∂θ_l_2∂γ_i,j^t-1()/∂θ_l_3
       -[1/γ_i,j^t-1(){1-γ_i,j^t-1()}+{γ_i,j^t-1(_0)-γ_i,j^t-1()}{1-2γ_i,j^t-1()}/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2]
            ×{∂^2γ_i,j^t-1()/∂θ_l_1∂θ_l_2∂γ_i,j^t-1()/∂θ_l_3+∂^2γ_i,j^t-1()/∂θ_l_2∂θ_l_3∂γ_i,j^t-1()/∂θ_l_1+∂^2γ_i,j^t-1()/∂θ_l_3∂θ_l_1∂γ_i,j^t-1()/∂θ_l_2}
       +γ_i,j^t-1(_0)-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂^3γ_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3 .

By the triangle inequality and Condition <ref>,

    |∂^3h_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3|≤    2(2C_1^-2+C_1^-3)|∂γ_i,j^t-1()/∂|_∞^3+C_1^-1|∂^3γ_i,j^t-1()/∂^3|_∞
       +3(C_1^-1+C_1^-2)|∂^2γ_i,j^t-1()/∂^2|_∞|∂γ_i,j^t-1()/∂|_∞
    ≤    2(2C_1^-2+C_1^-3)C_2^3+3(C_1^-1+C_1^-2)C_2^2+C_1^-1C_2=:C_* .

Write =(θ_1,…,θ_q)^ and _0=(θ_0,1,…,θ_0,q)^. By Taylor's theorem, (<ref>) and (<ref>), 

    h_i,j^t-1()-h_i,j^t-1(_0)=    ∂ h_i,j^t-1(_0)/∂^(-_0)+1/2(-_0)^∂^2h_i,j^t-1(_0)/∂∂^(-_0)+R_i,j^t-1()
    
    =   -1/2γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}|∂γ_i,j^t-1(_0)/∂^(_0-)|^2+R_i,j^t-1() ,

where 

    R_i,j^t-1()=    1/2∑_l_1=1^q(θ_l_1-θ_0,l_1)^3∫_0^1(1-v)^2∂^3h_i,j^t-1(_v)/∂θ_l_1^3  dv
       +3/2∑_l_1≠ l_2(θ_l_1-θ_0,l_1)^2(θ_l_2-θ_0,l_2)∫_0^1(1-v)^2∂^3h_i,j^t-1(_v)/∂θ_l_1^2∂θ_l_2  dv
       +3∑_l_1≠ l_2≠ l_3(θ_l_1-θ_0,l_1)(θ_l_2-θ_0,l_2)(θ_l_3-θ_0,l_3)∫_0^1(1-v)^2∂^3h_i,j^t-1(_v)/∂θ_l_1∂θ_l_2∂θ_l_3  dv

with _v=_0+v(-_0). Recall ℐ_i,j={l⊂[q]:γ_i,j^t-1() involves θ_l for any t∈[n]∖[m]}. We have 
∂γ_i,j^t-1()/∂θ_l_1≡0 if l_1∉ℐ_i,j, 
∂^2γ_i,j^t-1()/∂θ_l_1∂θ_l_2≡0 if l_1∉ℐ_i,j or l_2∉ℐ_i,j, and ∂^3γ_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3≡0 if l_1∉ℐ_i,j or l_2∉ℐ_i,j or l_3∉ℐ_i,j. By (<ref>), it holds that

    ∂^3h_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3≡0 if l_1∉ℐ_i,j or  l_2∉ℐ_i,j or l_3∉ℐ_i,j .

By Condition <ref>, it holds that 
|R_i,j^t-1()|≤ C_*s^3|_ℐ_i,j-_0,ℐ_i,j|_∞^3, 
which implies

    ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()=1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l{h_i,j^t-1(_0)-h_i,j^t-1()}
            ≥1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l1/2γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}|∂γ_i,j^t-1(_0)/∂^(_0-)|^2
                -C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_∞^3
            ≥2/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l|∂γ_i,j^t-1(_0)/∂^_ℐ_i,j(_0,ℐ_i,j-_ℐ_i,j)|^2-C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_∞^3

for any ∈. By  Condition <ref>, it holds with probability approaching one that 

    1/n-m∑_t=m+1^n∑_(i,j)∈𝒮_l|∂γ_i,j^t-1(_0)/∂^_ℐ_i,j(_0,ℐ_i,j-_ℐ_i,j)|^2≥ C_3∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2

for any l∈[q], which implies

    ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()≥    2C_3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2-C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_∞^3

for any ∈ and l∈[q] with probability approaching one. Since d:=sup_∈|-_0|_∞<2C_3/(C_*s^3), there exists a universal constant C̃>0 such that d<C̃<2C_3/(C_*s^3). Hence, it holds with probability approaching one that

    ℓ_n,p^(l)(_0)-ℓ_n,p^(l)()≥    2C_3-C̃C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2
       +C̃C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2-C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_∞^3
    ≥    2C_3-C̃C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2
       +C̃C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2-C_*s^3d/|𝒮_l|∑_(i,j)∈𝒮_l|_ℐ_i,j-_0,ℐ_i,j|_∞^2
    ≥    2C_3-C̃C_*s^3/|𝒮_l|∑_(i,j)∈𝒮_l|_0,ℐ_i,j-_ℐ_i,j|_2^2

for any ∈ and l∈[q]. 
We then have Proposition <ref> by selecting C̅=2C_3-C̃C_*s^3. 









 






































 §.§ Proof of Theorem <ref>
 
Recall S_ℋ,min=min_l∈ℋ|𝒮_l| for any ℋ⊂[q], and 
    {
       c_n,𝒢^2=    qlog(n|𝒮_l'|)/√(n|𝒮_l'|)+q^3/2log^3/2(n|𝒮_l'|)/√(n)|𝒮_l'| , 
    
       c_n,𝒢^^2=    qlog(nS_𝒢^,min)/√(nS_𝒢^,min)+q^3/2log^3/2(nS_𝒢^,min)/√(n)S_𝒢^,min ,
    .

where l'∈𝒢. 
To show Theorem <ref>, we first present the following lemma whose proof is given in Section <ref>. 


 Assume Conditions <ref> and <ref> hold. Then  

    max_l∈ℋsup_∈|ℓ̂_n,p^(l)() - ℓ_n,p^(l)()|=O_ p{qlog(nS_ℋ,min)/√(nS_ℋ,min)}+O_ p{q^3/2log^3/2(nS_ℋ,min)/√(n)S_ℋ,min}

for any ℋ⊂[q].







Notice that _*^(l)=(θ̂_*,1^(l), …,θ̂_*,q^(l) )^=max_∈ℓ̂_n,p^(l)() for any l∈[q]. Then

    ℓ_n,p^(l')(_0)-sup_∈|ℓ̂_n,p^(l')()-ℓ_n,p^(l')()|≤    ℓ̂_n,p^(l')(_0)
    ≤    ℓ̂_n,p^(l')(_*^(l'))≤ℓ_n,p^(l')(_*^(l'))+sup_∈|ℓ̂_n,p^(l')()-ℓ_n,p^(l')()| ,

which implies

    0≤ℓ_n,p^(l')(_0)-ℓ_n,p^(l')(_*^(l'))≤ 2sup_∈|ℓ̂_n,p^(l')()-ℓ_n,p^(l')()| .

Recall _𝒢=_*,𝒢^(l'). Selecting ℋ={l'} in Lemma <ref>, we have sup_∈|ℓ̂_n,p^(l')()-ℓ_n,p^(l')()|=O_ p(c_n,𝒢^2), which implies 

    ℓ_n,p^(l')(_0)-ℓ_n,p^(l')(_*^(l'))=O_ p(c_n,𝒢^2) .

For any diverging ϵ_n,p>0, if |_𝒢-_0,𝒢|_2≥ϵ_n,pc_n,𝒢, Proposition <ref> yields that

    ℓ_n,p^(l')(_0)-ℓ_n,p^(l')(_*^(l))≥C̅ϵ_n,p^2c_n,𝒢^2

with probability approaching one, which contradicts with (<ref>) and then implies |_𝒢-_0,𝒢|_2=O_ p( ϵ_n,pc_n,𝒢). Notice that we can select arbitrary slowly diverging ϵ_n,p. Following a standard result from
probability theory, we have |_𝒢-_0,𝒢|_2=O_ p(c_n,𝒢). 

For any l∈𝒢^, 






we also have 

    0≤ℓ_n,p^(l)(_0)-ℓ_n,p^(l)(_*^(l))≤ 2sup_∈|ℓ̂_n,p^(l)()-ℓ_n,p^(l)()| .

Recall _𝒢^=(θ̂_*,l^(l))_l∈𝒢^. Selecting ℋ=𝒢^ in Lemma <ref>, we have max_l∈𝒢^sup_∈|ℓ̂_n,p^(l)()-ℓ_n,p^(l)()|=O_ p(c_n,𝒢^^2), which implies 

    max_l∈𝒢^{ℓ_n,p^(l)(_0)-ℓ_n,p^(l)(_*^(l))}=O_ p(c_n,𝒢^^2) .

For any diverging ϵ_n,p>0, if |_𝒢^-_0,𝒢^|_∞≥ϵ_n,pc_n,𝒢^, Proposition <ref> yields that

    max_l∈𝒢^{ℓ_n,p^(l)(_0)-ℓ_n,p^(l)(_*^(l))}≥C̅ϵ_n,p^2c_n,𝒢^^2

with probability approaching one, which contradicts with (<ref>) and then implies |_𝒢^-_0,𝒢^|_∞=O_ p( ϵ_n,pc_n,𝒢^). Notice that we can select arbitrary slowly diverging ϵ_n,p. Following a standard result from
probability theory, we have |_𝒢^-_0,𝒢^|_∞=O_ p(c_n,𝒢^). We complete the proof of Theorem <ref>.  





























 §.§ Proof of Proposition <ref>
 

Given _l specified in Condition <ref>, define

    f_t^(l)()=_l^_t^(l)()

for any t∈[n]∖[m] and ∈. 
To show Proposition <ref>, we need the following lemmas whose proofs are given in Sections <ref>–<ref>.


Assume Condition <ref> holds. Then

    |1/n-m∑_t=m+1^n{f̂_t^(l)(θ_0,l,_-l)-f̂_t^(l)(_0)}|≤τ |-_0|_1+C_*|_l|_1 |-_0|_1^2

for any l∈[q], where C_* is a universal constant specified in (<ref>). 



Assume Conditions <ref> and <ref> hold. Then

    max_l∈ℋ|1/n-m∑_t=m+1^nf̂_t^(l)(_0)|
              =O_ p{√(log(1+ |ℋ|)log(n|ℋ|))/√(nS_ℋ,min)}+O_ p{√(log(1+ |ℋ|))log(n|ℋ|)/√(n)S_ℋ,min}

for any ℋ⊂[q]. 



Assume Condition <ref> holds. Then

    sup_θ_l∈ B(θ̃_l,r̃)|1/n-m∑_t=m+1^nf̂_t^(l)(θ_l,_-l)/∂θ_l-1|≤τ+C_*r̃|_l|_1

for any l∈[q], where C_* is a universal constant specified in (<ref>). 


Recall |_𝒢-_0,𝒢|_2=O_ p(c_n,𝒢) and |_𝒢^-_0,𝒢^|_∞=O_ p(c_n,𝒢^) with 
    {
       c_n,𝒢^2=    qlog(n|𝒮_l'|)/√(n|𝒮_l'|)+q^3/2log^3/2(n|𝒮_l'|)/√(n)|𝒮_l'| , 
    
       c_n,𝒢^^2=    qlog(nS_𝒢^,min)/√(nS_𝒢^,min)+q^3/2log^3/2(nS_𝒢^,min)/√(n)S_𝒢^,min ,
    .

where l'∈𝒢. Then |-_0|_1=|_𝒢-_0,𝒢|_1+|_𝒢^-_0,𝒢^|_1=O_ p(|𝒢|^1/2c_n,𝒢)+O_ p(|𝒢^|c_n,𝒢^)=O_ p(Δ_n^1/2). Based on Lemmas <ref> and <ref>, due to τ≲Δ_n^1/2 and |𝒢|+|𝒢^|=q, we have

    max_l∈𝒢|1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)|
                =O_ p(Δ_n)+O_ p{√(log(1+|𝒢|)log(n|𝒢|))/√(n|𝒮_l'|)}+O_ p{√(log(1+|𝒢|))log(n|𝒢|)/√(n)|𝒮_l'|}
                =O_ p(Δ_n)

for some l'∈𝒢, and 

    max_l∈𝒢^|1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)|
                =O_ p(Δ_n)+O_ p{√(log(1+|𝒢^|)log(n|𝒢^|))/√(nS_𝒢^,min)}+O_ p{√(log(1+|𝒢^|))log(n|𝒢^|)/√(n)S_𝒢^,min}
                =O_ p(Δ_n) .

Due to θ̌_l=min_θ_l∈ B(θ̃_l,r̃)|(n-m)^-1∑_t=m+1^nf̂_t^(l)(θ_l,_-l)|^2 and θ_0,l∈ B(θ̃_l,r̃) with probability approaching one, then

    |1/n-m∑_t=m+1^nf̂_t^(l)(θ̌_l,_-l)|≤|1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)|

with probability approaching one, 
which implies

    max_l∈[q]|1/n-m∑_t=m+1^n{f̂_t^(l)(θ̌_l,_-l)-f̂_t^(l)(θ_0,l,_-l)}|=O_ p(Δ_n) .

Notice that

    1/n-m∑_t=m+1^n{f̂_t^(l)(θ̌_l,_-l)-f̂_t^(l)(θ_0,l,_-l)}={1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}(θ̌_l-θ_0,l)

for some θ̅_l on the joint line between θ̌_l and θ_0,l. 
By Lemma <ref> and Condition <ref>, due to τ=o(1) and r̃=o(1), we know

    min_l∈[q]{1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}>1/2

with probability approaching one. Hence, by (<ref>), we have |-_0|_∞=O_ p(Δ_n). We complete the proof of Proposition <ref>.  




 §.§ Proof of Theorem <ref>
 
Recall f̂_t^(l)()=_l^_t^(l)(). By the definition of _l given in (<ref>), we have 

    max_l∈[q]|1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l|_∞≤τ .

It follows from the Taylor expansion that

    |1/n-m∑_t=m+1^n{∂f̂_t^(l)()/∂_-l-∂f̂_t^(l)()/∂_-l}|_∞≤|1/n-m∑_t=m+1^n∂^2f̂_t^(l)()/∂_-l∂^|_∞|-|_1 ,

where  is on the joint line between  and . 





Notice that ∈ B(θ̃_1,r̃)×⋯× B(θ̃_q,r̃). 
Following the same arguments for deriving (<ref>) in Section <ref> for the proof of Lemma <ref>, we know 

    sup_∈ B(θ̃_1,r̃)×⋯× B(θ̃_q,r̃)max_l∈[q]|1/n-m∑_t=m+1^n∂^2f̂_t^(l)()/∂_-l∂^|_∞≤ C_*max_l∈[q]|_l|_1

for some universal constant C_* specified in (<ref>), which implies

    max_l∈[q]|1/n-m∑_t=m+1^n{∂f̂_t^(l)()/∂_-l-∂f̂_t^(l)()/∂_-l}|_∞≤ C_*|-|_1max_l∈[q]|_l|_1 .

By Condition <ref>, we have max_l∈[q]|_l|_1=O_ p(1), which implies 

    max_l∈[q]|1/n-m∑_t=m+1^n{∂f̂_t^(l)()/∂_-l-∂f̂_t^(l)()/∂_-l}|_∞≤|-|_1· O_ p(1) .

Hence, it holds that 

    max_l∈[q]|1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l|_∞≤τ+ |-|_1· O_ p(1) .

Repeating the arguments for deriving Lemma <ref> in Section <ref>, we can also show

    max_l∈[q]|1/n-m∑_t=m+1^n{f̂_t^(l)(θ_0,l,_-l)-f̂_t^(l)(_0)}|
                        ≤{τ+|-|_1· O_ p(1)} |-_0|_1+ |-_0|_1^2· O_ p(1)
                        ≤{τ+|-_0|_1· O_ p(1)} |-_0|_1+ |-_0|_1^2· O_ p(1) .

Together with Lemma <ref>, it holds that

    max_l∈[q]|1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)|≤ O_ p(qΔ_n) .

 Since θ̂_l=min_θ_l∈ B(θ̌_l,ř)|(n-m)^-1∑_t=m+1^nf̂_t^(l)(θ_l,_-l)|^2 and θ_0,l∈ B(θ̌_l,ř) with probability approaching one, then
 
    |1/n-m∑_t=m+1^nf̂_t^(l)(θ̂_l,_-l)|≤|1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)|

with probability approaching one, which implies 

    max_l∈[q]|1/n-m∑_t=m+1^n{f̂_t^(l)(θ̂_l,_-l)-f̂_t^(l)(θ_0,l,_-l)}|=O_ p(qΔ_n) .

Following the same arguments for deriving Lemma <ref> in Section <ref> and noting (<ref>), we can also have

    sup_θ_l∈ B(θ̌_l,ř)|1/n-m∑_t=m+1^nf̂_t^(l)(θ_l,_-l)/∂θ_l-1|≤τ+o_ p(1)+C_*ř|_l|_1

for any l∈[q]. Due to τ=o(1), ř=o(1) and max_l∈[q]|_l|_1=O_ p(1), we have

    max_l∈[q]sup_θ_l∈ B(θ̌_l,ř)|1/n-m∑_t=m+1^nf̂_t^(l)(θ_l,_-l)/∂θ_l-1|=o_ p(1) .

Due to

    1/n-m∑_t=m+1^n{f̂_t^(l)(θ̂_l,_-l)-f̂_t^(l)(θ_0,l,_-l)}={1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}(θ̂_l-θ_0,l)

for some θ̅_l on the joint line between θ̂_l and θ_0,l, by (<ref>) and (<ref>), we have |-_0|_∞=O_ p(qΔ_n), which implies |-|_∞≤ |-_0|_∞+|-_0|_∞=O_ p(qΔ_n). Hence, |-|_∞≪ř with probability approaching one. Therefore, 

    0={1/n-m∑_t=m+1^nf̂_t^(l)(θ̂_l,_-l)}{1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̂_l,_-l)/∂θ_l}

with probability approaching one. Together with (<ref>), we have

    1/n-m∑_t=m+1^nf̂_t^(l)(θ̂_l,_-l)=0

with probability approaching one.


By (<ref>) and (<ref>), due to |-_0|_1=O_ p(qΔ_n), |-_0|_1=O_ p(Δ_n^1/2) and τ≲Δ_n^1/2, then

    θ̂_l-θ_0,l=   -{1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}^-1{1/n-m∑_t=m+1^nf̂_t^(l)(θ_0,l,_-l)}
    
    =   -{1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}^-1{1/n-m∑_t=m+1^nf̂_t^(l)(_0)}
       +O_ p(qΔ_n^3/2)+O_ p(q^2Δ_n^2)

with probability approaching one. As shown in (<ref>) in Section <ref> for the proof of Lemma <ref>, 

    |1/n-m∑_t=m+1^n{f̂_t^(l)(_0)-f_t^(l)(_0)}|
                =O_ p{ω_n√((log q)log(qn))/√(n|𝒮_l|)}+O_ p{ω_n(log q)^1/2log(qn)/√(n)|𝒮_l|}=o_ p(1/√(n|𝒮_l|)) ,

where the last step is based on Condition <ref>. 
Hence, 
    θ̂_l-θ_0,l
    =   -{1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l}^-1{1/n-m∑_t=m+1^nf_t^(l)(_0)}
       +O_ p(qΔ_n^3/2)+O_ p(q^2Δ_n^2)+o_ p(1/√(n|𝒮_l|))

with probability approaching one. 

Write 

    Q_l,t=1/|𝒮_l|^1/2∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_0)/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}{_l^∂γ_i,j^t-1(_0)/∂} .

Then

    1/n-m∑_t=m+1^nf_t^(l)(_0)
    =1/(n-m)|𝒮_l|^1/2∑_t=m+1^nQ_l,t .

In the sequel, we will use the martingale central limit theorem to establish the asymptotic distribution of (n-m)^-1/2∑_t=m+1^nQ_l,t. 
Denote by ℙ_ℱ_t-1(·) and 𝔼_ℱ_t-1(·), respectively, the conditional probability measure and the conditional expectation given ℱ_t-1 with the unknown true parameter vector _0. By Conditions <ref> and <ref>, it holds that

    max_(i,j)∈𝒮_l|X_i,j^t-γ_i,j^t-1(_0)/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}{_l^∂γ_i,j^t-1(_0)/∂}|≤ C_1^-1C_2C_4=:C_** .
 




It follows from the Bernstein inequality that

    ℙ_ℱ_t-1(|Q_l,t|≥ x)≤ 2exp(-3|𝒮_l|^1/2x^2/6|𝒮_l|^1/2C_**^2+2xC_**)

for any x>0, 




which implies, for any δ>0,

    ∑_t=m+1^n𝔼_ℱ_t-1{(Q_l,t/√(n-m))^2I(|Q_l,t|/√(n-m)≥δ)}
              =1/n-m∑_t=m+1^n𝔼_ℱ_t-1{Q_l,t^2I(|Q_l,t|≥δ√(n-m))}→0

in probability as n→∞. Meanwhile, by Condition <ref>, we also have

    1/n-m∑_t=m+1^n𝔼_ℱ_t-1(Q_l,t^2)
                  =1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_l1/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}{_l^∂γ_i,j^t-1(_0)/∂}^2→κ_l





in probability as n→∞. By Conditions <ref> and <ref>, we know κ_l is a almost surely bounded random variable. Corollary 3.1 of <cit.> implies

    1/√(n-m)∑_t=m+1^nQ_l,t→√(κ_l)· Z

in distribution as n→∞, where Z is a standard normally distributed random variable independent of κ_l. By (<ref>) and (<ref>), due to  

    1/n-m∑_t=m+1^n∂f̂_t^(l)(θ̅_l,_-l)/∂θ_l→1

in probability which is obtained in (<ref>), it holds that
 
    √(n|𝒮_l|)(θ̂_l-θ_0,l)→√(κ_l)· Z

in distribution as n→∞, provided that √(n|𝒮_l|)max{qΔ_n^3/2,q^2Δ_n^2}=o(1). We complete the proof of Theorem <ref>. 







 §.§ Proofs of auxiliary lemmas




  §.§.§ Proof of Lemma <ref>

Without loss of generality, we assume =[-C,C]^q for some constant C>0. For given ϵ>0 which will be specified later, we partition [-C,C] into K=⌈ 2C/ϵ⌉ sub-intervals B_1,…,B_K with equal length, where the length of each B_k does not exceed ϵ. Based on such defined B_1,…,B_K, we can partition  as follows:

    =⋃_k_1=1^K⋯⋃_k_q=1^KB_k_1×⋯× B_k_q ,

which includes K^q hyper-rectangles ℬ_1,…,ℬ_K^q. For each given u∈[K^q], there exists (k_1,u,…,k_q,u)∈[K]^q such that ℬ_u=B_k_1,u×⋯× B_k_q,u. Let _u be the center of ℬ_u. 

For each ∈ℬ_u, since γ_i,j^t-1() only depends on θ_l with l∈ℐ_i,j, it follows from the Taylor expansion that  

    ℓ̂_n,p^(l)()-ℓ̂_n,p^(l)(_u)
              =[1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_u)/γ_i,j^t-1(_u){1-γ_i,j^t-1(_u)}∂γ_i,j^t-1(_u)/∂^](-_u)
              =1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_u)/γ_i,j^t-1(_u){1-γ_i,j^t-1(_u)}∂γ_i,j^t-1(_u)/∂_ℐ_i,j^(_ℐ_i,j-_u,ℐ_i,j) ,

where _u ∈ℬ_u is on the joint line between  and _u. Write =(θ_1,…,θ_q)^ and _u=(θ_u,1,…,θ_u,q)^. 
By Conditions <ref> and <ref>, it holds that 

    |ℓ̂_n,p^(l)() - ℓ̂_n,p^(l)(_u)|
           
             ≤1/(n-m)|𝒮_l|∑_t = m+1^n ∑_(i,j)∈𝒮_l∑_l∈ℐ_i,j| X_i,j^t- γ_i,j^t-1(_u)/γ_i,j^t-1(_u){1- γ_i,j^t-1(_u)}||∂γ_i,j^t-1(_u)/∂θ_l| |θ_l - θ_u,l| 
           
             ≤C_1^-1C_2ϵ/(n-m)|𝒮_l|∑_t = m+1^n ∑_(i,j)∈𝒮_l∑_l∈ℐ_i,j  |X_i,j^t- γ_i,j^t-1(_u)|≤ sC_1^-1C_2ϵ

for any ∈ℬ_u and l∈[q]. Analogously, we also have 

    sup_∈ℬ_u|ℓ_n,p^(l)()-ℓ_n,p^(l)(_u)|≤ sC_1^-1C_2ϵ .

Therefore, by the triangle inequality, it holds that

    max_l∈ℋsup_∈|ℓ̂_n,p^(l)() - ℓ_n,p^(l)()| =max_l∈ℋmax_u∈[K^q]sup_∈ℬ_u |ℓ̂_n,p^(l)() - ℓ_n,p^(l)()| 
             ≤max_l∈ℋmax_u∈[K^q] |ℓ̂_n,p^(l)(_u) - ℓ_n,p^(l)(_u)| +max_l∈ℋmax_u∈[K^q]sup_∈ℬ_u|ℓ̂_n,p^(l)() - ℓ̂_n,p^(l)(_u)|
                 +max_l∈ℋmax_u∈[K^q]sup_∈ℬ_u|ℓ_n,p^(l)()-ℓ_n,p^(l)(_u)|
             ≤max_l∈ℋmax_u∈[K^q] |ℓ̂_n,p^(l)(_u) - ℓ_n,p^(l)(_u)| +  2sC_1^-1C_2ϵ .

 For each l∈[q], u∈[K^q] and t∈[n]∖[m], define

    Q_u,t^(l)=    1/|𝒮_l|^1/2∑_(i,j)∈𝒮_l[ {X_i,j^t - γ_i,j^t-1(_0)}log{γ_i,j^t-1(_u)/1- γ_i,j^t-1(_u)}] .

Then 

    ℓ̂_n,p^(l)(_u) - ℓ_n,p^(l)(_u) =1/(n-m)|𝒮_l|^1/2∑_t = m+1^n Q_u,t^(l) .




Notice that   

    |X_i,j^t-γ_i,j^t-1(_0)||log{γ_i,j^t-1(_u)/1-γ_i,j^t-1(_u)}|≤max_i≠ j|log{γ_i,j^t-1(_u)/1-γ_i,j^t-1(_u)}| .

Due to 

    log[γ_i,j^t-1(_u){1-γ_i,j^t-1(_u)}]<log{γ_i,j^t-1(_u)/1-γ_i,j^t-1(_u)}<-log[γ_i,j^t-1(_u){1-γ_i,j^t-1(_u)}] ,

by Condition <ref>, we have

    max_u∈[K^q]max_t∈[n]∖[m]|X_i,j^t-γ_i,j^t-1(_0)||log{γ_i,j^t-1(_u)/1-γ_i,j^t-1(_u)}| 
                      ≤max_u∈[K^q]max_t∈[n]∖[m]max_i≠ j|log[γ_i,j^t-1(_u){1-γ_i,j^t-1(_u)}]|≤log C_1^-1 .

Denote by ℙ_ℱ_t-1(·) the conditional probability measure given ℱ_t-1 with the unknown true parameter vector _0. It follows from the Bernstein inequality that

    max_u∈[K^q]max_t∈[n]∖[m]ℙ_ℱ_t-1{|Q_u,t^(l)|≥ x}
                     ≤ 2exp{-3|𝒮_l|^1/2x^2/6|𝒮_l|^1/2(log C_1^-1)^2+2xlog C_1^-1}

for any x>0. Furthermore, due to ℙ{|Q_u,t^(l)|≥ x}=𝔼[ℙ_ℱ_t-1{|Q_u,t^(l)|≥ x}], we also have 

    max_u∈[K^q]max_t∈[n]∖[m]ℙ{|Q_u,t^(l)|≥ x}≤ 2exp{-3|𝒮_l|^1/2x^2/6|𝒮_l|^1/2(log C_1^-1)^2+2xlog C_1^-1}

for any x>0. Let

    Q̃_u,t^(l)=    Q_u,t^(l)I{|Q_u,t^(l)|≤ M}-𝔼_ℱ_t-1[Q_u,t^(l)I{|Q_u,t^(l)|≤ M}]

for some diverging M>0 specified later. Notice that {Q̃_u,t^(l)}_t∈[n]∖[m] is a martingale difference sequence with max_t∈[n]∖[m]|Q̃_u,t^(l)|≤ 2M. By the Azuma's inequality (, ; see also Theorem 3.1 of <cit.>), we have

    max_l∈[q]max_u∈[K^q]ℙ{|1/n-m∑_t=m+1^nQ̃_u,t^(l)|≥ x}≤ 2exp{-(n-m)x^2/8M^2}

for any x>0. By (<ref>), 

    max_u∈[K^q]ℙ[|1/n-m∑_t=m+1^nQ_u,t^(l)I{|Q_u,t^(l)|>M}|>0]≤max_u∈[K^q]ℙ{max_t∈[n]∖[m]|Q_u,t^(l)|>M}
                ≤ 2nexp{-3|𝒮_l|^1/2M^2/6|𝒮_l|^1/2(log C_1^-1)^2+2Mlog C_1^-1} .

Together with (<ref>), by the Bonferroni inequality, we have

    ℙ(max_l∈ℋmax_u∈[K^q]|1/n-m∑_t=m+1^n[Q̃_u,t^(l)+Q_u,t^(l)I{|Q_u,t^(l)|>M}]|>x)
                  ≤ 2|ℋ|K^qexp{-(n-m)x^2/8M^2}
                     +2|ℋ|K^qnmax_l∈ℋexp{-3|𝒮_l|^1/2M^2/6|𝒮_l|^1/2(log C_1^-1)^2+2Mlog C_1^-1}

for any x>0. Recall S_ℋ,min=min_l∈ℋ|𝒮_l| and |ℋ|≤ q. Selecting 

    M= Cmax(√(qlog K), qlog K/√(S_ℋ,min), √(log n), log n/√(S_ℋ,min))

for some sufficiently large constant C>0, it holds that 

    max_l∈ℋmax_u∈[K^q]|1/n-m∑_t=m+1^n[Q̃_u,t^(l)+Q_u,t^(l)I{|Q_u,t^(l)|>M}]|=O_ p(M√(qlog K/n))
               =O_ p(qlog K/√(n))+O_ p{(qlog K)^3/2/√(nS_ℋ,min)}
                  +O_ p{√(q(log K)(log n))/√(n)}+O_ p{(qlog K)^1/2(log n)/√(nS_ℋ,min)} .

On the other hand, by (<ref>), it holds that

    𝔼_ℱ_t-1[|Q_u,t^(l)|I{|Q_u,t^(l)|>M}]=Mℙ_ℱ_t-1{|Q_u,t^(l)|>M}+∫_M^∞ℙ_ℱ_t-1(|Q_u,t|>x)  dx
             ≤ 2Mexp{-3|𝒮_l|^1/2M^2/6|𝒮_l|^1/2(log C_1^-1)^2+2Mlog C_1^-1}
                +2∫_M^∞exp{-x^2/4(log C_1^-1)^2}  dx+2∫_M^∞exp(-3|𝒮_l|^1/2x/4log C_1^-1)  dx
             ≤ 2Mexp{-3|𝒮_l|^1/2M^2/6|𝒮_l|^1/2(log C_1^-1)^2+2Mlog C_1^-1}
                +4(log C_1^-1)^2M^-1exp{-M^2/4(log C_1^-1)^2}+8log C_1^-1/3|𝒮_l|^1/2exp(-3|𝒮_l|^1/2M/4log C_1^-1) ,

which implies

    max_u∈[K^q]|1/n-m∑_t=m+1^n𝔼_ℱ_t-1[Q_u,t^(l)I{|Q_u,t^(l)|>M}]|
               ≤ 2Mexp{-3|𝒮_l|^1/2M^2/6|𝒮_l|^1/2(log C_1^-1)^2+2Mlog C_1^-1}
                  +4(log C_1^-1)^2exp{-M^2/4(log C_1^-1)^2}+8log C_1^-1/3|𝒮_l|^1/2exp(-3|𝒮_l|^1/2M/4log C_1^-1) .

Due to S_ℋ,min=min_l∈ℋ|𝒮_l|→∞ and M→∞ satisfying (<ref>), by (<ref>), we have

    max_l∈ℋmax_u∈[K^q]|1/n-m∑_t=m+1^n𝔼_ℱ_t-1[Q_u,t^(l)I{|Q_u,t^(l)|>M}]|
                  ≲exp(-C̃_1M^2)+exp(-C̃_2MS_ℋ,min^1/2)=o(n^-1/2) ,

where C̃_1>0 and C̃_2>0 are two universal constants. Due to Q_u,t^(l)=Q̃_u,t^(l)+Q_u,t^(l)I{|Q_u,t^(l)|>M}-𝔼_ℱ_t-1,_0[Q_u,t^(l)I{|Q_u,t^(l)|>M}], together with (<ref>), it holds that 

    max_l∈ℋmax_u∈[K^q]|1/n-m∑_t=m+1^nQ_u,t^(l)|=    O_ p(qlog K/√(n))+O_ p{(qlog K)^3/2/√(nS_ℋ,min)}
       +O_ p{√(q(log K)(log n))/√(n)}+O_ p{(qlog K)^1/2(log n)/√(nS_ℋ,min)} .

By (<ref>), we have 

    max_l∈ℋmax_u∈[K^q]|ℓ̂_n,p^(l)(_u) - ℓ_n,p^(l)(_u)|=    O_ p(qlog K/√(nS_ℋ,min))+O_ p{(qlog K)^3/2/√(n)S_ℋ,min}
       +O_ p{√(q(log K)(log n))/√(nS_ℋ,min)}+O_ p{(qlog K)^1/2(log n)/√(n)S_ℋ,min} .

Recall K≍ϵ^-1. Together with (<ref>), it holds that 

    max_l∈ℋsup_∈|ℓ̂_n,p^(l)() - ℓ_n,p^(l)()|≤    2sC_1^-1C_2ϵ+O_ p(qlog K/√(nS_ℋ,min))+O_ p{(qlog K)^3/2/√(n)S_ℋ,min}
       +O_ p{√(q(log K)(log n))/√(nS_ℋ,min)}+O_ p{(qlog K)^1/2(log n)/√(n)S_ℋ,min} .

Due to s≤ q, with selecting ϵ≍ n^-1/2S_ℋ,min^-1/2, we have

    max_l∈ℋsup_∈|ℓ̂_n,p^(l)() - ℓ_n,p^(l)()|=O_ p{qlog(nS_ℋ,min)/√(nS_ℋ,min)}+O_ p{q^3/2log^3/2(nS_ℋ,min)/√(n)S_ℋ,min} .

We complete the proof of Lemma <ref>. 






  §.§.§ Proof of Lemma <ref>


By the Taylor expansion, we have

    1/n-m∑_t=m+1^n{f̂_t^(l)(θ_0,l,_-l)-f̂_t^(l)(_0)}
                 =1/n-m∑_t=m+1^n∂f̂_t^(l)(θ_0,l,_-l)/∂_-l^(_-l-_0,-l)
                 =1/n-m∑_t=m+1^n∂f̂_t^(l)(θ_0,l,_-l)/∂_-l^(_-l-_0,-l)_R_1,l
                    +1/n-m∑_t=m+1^n{∂f̂_t^(l)(θ_0,l,_-l)/∂_-l^-∂f̂_t^(l)(θ_0,l,_-l)/∂_-l^}(_-l-_0,-l)_R_2,l ,

where _-l is on the joint line between _0,-l and _-l. 

For R_1,l, by the Taylor expansion, it holds that

    R_1,l=    1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l^(_-l-_0,-l)
       +(θ_0,l-θ̃_l){1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂_-l^}(_-l-_0,-l) ,

where θ̅_l is on the joint line between θ_0,l and θ̃_l. Recall f̂_t^(l)()=_l^_t^(l)(). By the definition of _l given in (<ref>), we have 

    |1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l|_∞≤τ ,

which implies

    |1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l^(_-l-_0,-l)|≤    |1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂_-l|_∞ |_-l-_0,-l|_1
    ≤    τ |_-l-_0,-l|_1 .

For any k∈[q], due to

    ∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂θ_k=_l^∂^2_t^(l)(θ̅_l,_-l)/∂θ_l∂θ_k ,

we then have

    |∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂θ_k|≤     |_l|_1|∂^2_t^(l)(θ̅_l,_-l)/∂θ_l∂θ_k|_∞
    ≤    |_l|_1/|𝒮_l|∑_(i,j)∈𝒮_l|∂^2/∂θ_l∂θ_k[X_i,j^t-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂γ_i,j^t-1()/∂]|_∞

Notice that

    ∂^2/∂θ_l_1∂θ_l_2[X_i,j^t-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂γ_i,j^t-1()/∂θ_l_3]
             =2[1-2γ_i,j^t-1()/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2+X_i,j^t-γ_i,j^t-1()/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2
                     +{X_i,j^t-γ_i,j^t-1()}{1-2γ_i,j^t-1()}^2/{γ_i,j^t-1()}^3{1-γ_i,j^t-1()}^3]∂γ_i,j^t-1()/∂θ_l_1∂γ_i,j^t-1()/∂θ_l_2∂γ_i,j^t-1()/∂θ_l_3
                -[1/γ_i,j^t-1(){1-γ_i,j^t-1()}+{X_i,j^t-γ_i,j^t-1()}{1-2γ_i,j^t-1()}/{γ_i,j^t-1()}^2{1-γ_i,j^t-1()}^2]
    
    
                      ×{∂^2γ_i,j^t-1()/∂θ_l_1∂θ_l_2∂γ_i,j^t-1()/∂θ_l_3+∂^2γ_i,j^t-1()/∂θ_l_2∂θ_l_3∂γ_i,j^t-1()/∂θ_l_1+∂^2γ_i,j^t-1()/∂θ_l_3∂θ_l_1∂γ_i,j^t-1()/∂θ_l_2}
                +X_i,j^t-γ_i,j^t-1()/γ_i,j^t-1(){1-γ_i,j^t-1()}∂^3γ_i,j^t-1()/∂θ_l_1∂θ_l_2∂θ_l_3 .

By the triangle inequality and Condition <ref>,  we know 

    max_k∈[q]|∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂θ_k|≤ C_*|_l|_1

for some universal constant C_* specified in (<ref>), which implies

    |(θ_0,l-θ̃_l){1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂_-l^}(_-l-_0,-l)|
               ≤ |θ_0,l-θ̃_l||1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l∂_-l^|_∞|_-l-_0,-l|_1
               ≤ C_*|_l|_1 |θ_0,l-θ̃_l||_-l-_0,-l|_1 .

Together with (<ref>), it holds that 
|R_1,l|≤τ |_-l-_0,-l|_1+C_*|_l|_1 |θ_0,l-θ̃_l||_-l-_0,-l|_1. 

For R_2,l, by the Taylor expansion, we have

    R_2,l=(_-l-_-l)^{1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ_0,l,_-l)/∂_-l∂_-l^}(_-l-_0,-l)

for some _-l on the joint line between _-l and _-l, which implies

    |R_2,l|≤ |_-l-_-l|_1|_-l-_0,-l|_1|1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ_0,l,_-l)/∂_-l∂_-l^|_∞ .

Since _-l is on the joint line between _0,-l and _-l, then |_-l-_-l|_1≤ |_-l-_0,-l|_1. Parallel to (<ref>), we can also show

    |1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ_0,l,_-l)/∂_-l∂_-l^|_∞≤ C_*|_l|_1 .

Hence, 
|R_2,l|≤ C_*|_l|_1|_-l-_0,-l|_1^2. Then

    |1/n-m∑_t=m+1^n{f̂_t^(l)(θ_0,l,_-l)-f̂_t^(l)(_0)}|≤ |R_1,l|+|R_2,l|
                ≤τ |_-l-_0,-l|_1+C_*|_l|_1 |θ_0,l-θ̃_l||_-l-_0,-l|_1+C_*|_l|_1|_-l-_0,-l|_1^2
                = τ |_-l-_0,-l|_1+C_*|_l|_1 |-_0|_1|_-l-_0,-l|_1
                ≤τ |-_0|_1+C_*|_l|_1 |-_0|_1^2 .

We complete the proof of Lemma <ref>. 




  §.§.§ Proof of Lemma <ref>


Due to f̂_t^(l)()=_l^_t^(l)() and f_t^(l)()=_l^_t^(l)(), then

    |1/n-m∑_t=m+1^n{f̂_t^(l)(_0)-f_t^(l)(_0)}|≤ |_l-_l|_1|1/n-m∑_t=m+1^n_t^(l)(_0)|_∞ .

Write (G_1^(l),…,G_q^(l))^=(n-m)^-1∑_t=m+1^n_t^(l)(_0). Due to 

    _t^(l)(_0)=1/|𝒮_l|∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_0)/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}∂γ_i,j^t-1(_0)/∂ ,

we then have

    G_k^(l)=    1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_0)/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}∂γ_i,j^t-1(_0)/∂θ_k
    
    :=    1/(n-m)|𝒮_l|^1/2∑_t=m+1^nQ̌_k,t^(l)

for any k∈[q]. Using the same arguments for deriving (<ref>), we can also show

    max_l∈ℋmax_k∈[q]|G_k^(l)|=O_ p{√((log q)log(qn))/√(nS_ℋ,min)}+O_ p{(log q)^1/2log(qn)/√(n)S_ℋ,min} ,

where S_ℋ,min=min_l∈ℋ|𝒮_l|. Together with Condition <ref>, it holds that

    max_l∈ℋ|1/n-m∑_t=m+1^n{f̂_t^(l)(_0)-f_t^(l)(_0)}|
               =O_ p{ω_n√((log q)log(qn))/√(nS_ℋ,min)}+O_ p{ω_n(log q)^1/2log(qn)/√(n)S_ℋ,min} .

Notice that

    1/n-m∑_t=m+1^nf_t^(l)(_0)=1/(n-m)|𝒮_l|∑_t=m+1^n∑_(i,j)∈𝒮_lX_i,j^t-γ_i,j^t-1(_0)/γ_i,j^t-1(_0){1-γ_i,j^t-1(_0)}{_l^∂γ_i,j^t-1(_0)/∂} .

By Conditions <ref> and <ref>, we know

    max_l∈[q]max_t∈[n]∖[m]max_(i,j)∈𝒮_l|_l^∂γ_i,j^t-1(_0)/∂|≤ C_2C_4 .

Using the same arguments for deriving (<ref>), we can also show

    max_l∈ℋ|1/n-m∑_t=m+1^nf_t^(l)(_0)|
                =O_ p{√(log(1+|ℋ|)log(n|ℋ|))/√(nS_ℋ,min)}+O_ p{√(log(1+|ℋ|))log(n|ℋ|)/√(n)S_ℋ,min} .

Together with (<ref>), due to ω_n(log q)^1/2log(qn)=o(1), we have Lemma <ref>. 



  §.§.§ Proof of Lemma <ref>


For any θ_l∈ B(θ̃_l,r̃), by the Taylor expansion, we have

    1/n-m∑_t=m+1^n{∂f̂_t^(l)(θ_l,_-l)/∂θ_l-∂f̂_t^(l)()/∂θ_l}={1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l^2}(θ_l-θ̃_l)

for some θ̅_l on the joint line between θ_l and θ̃_l, which implies

    |1/n-m∑_t=m+1^n{∂f̂_t^(l)(θ_l,_-l)/∂θ_l-∂f̂_t^(l)()/∂θ_l}|=|1/n-m∑_t=m+1^n∂^2f̂_t^(l)(θ̅_l,_-l)/∂θ_l^2||θ_l-θ̃_l| .

Due to θ_l∈ B(θ̃_l,r̃), by (<ref>), it holds that

    sup_θ_l∈ B(θ̃_l,r̃)|1/n-m∑_t=m+1^n{∂f̂_t^(l)(θ_l,_-l)/∂θ_l-∂f̂_t^(l)()/∂θ_l}|≤ C_*r̃|_l|_1 .

Since

    1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂θ_l=_l^/n-m∑_t=m+1^n∂_t^(l)()/∂θ_l ,

by (<ref>), we know

    |1/n-m∑_t=m+1^n∂f̂_t^(l)()/∂θ_l-1|≤τ ,

which implies

    sup_θ_l∈ B(θ̃_l,r̃)|1/n-m∑_t=m+1^n∂f̂_t^(l)(θ_l,_-l)/∂θ_l-1|≤τ+C_*r̃|_l|_1 .





We complete the proof of Lemma <ref>.  
















































































































































































































































































































































































































































































































































§ ADDITIONAL SIMULATION RESULTS
  

Further to Section <ref>, we report more simulation results on the transitivity model.




 §.§ More plots for the simulated transitivity models

 
We provide two additional figures for the simulated transitivity models
in Section <ref>.
Figure <ref> and <ref> display the time series plots of simulated  {D_t}_t=2^200, {D_1,t}_t=2^200, {D_0,t}_t=2^200, {D̅_t}_t=2^200, {D̅_1,t}_t=2^200 and {D̅_0,t}_t=2^200 when p=100 and 150, respectively. All settings are stationary, among which (ξ_i, η_i, a, b) = (0.7,0.8,30,15) shows the most frequent edge-changing behaviour.












 §.§ Stationarity and ergodicity


In this section, we investigate the stationarity and ergodicity of our proposed dynamic models with more synthetic data. We simulate {_1, …, _n} from two dependent AR(1) networks, i.e. density-dependent model as in (<ref>), and transitivity model in (<ref>).
For each sample, we generate a sequence of length n +200 and discard the first 200 observations. For the density-dependent model, we consider 625 settings with a_0=a_1∈{0.2,0.4,0.6,0.8,1}, b_0=b_1∈{0.2,0.4,0.6,0.8,1}, ξ_1 = …= ξ_p ∈{0.1,0.3,0.5,0.7,0.9} and η_1 = … =η_p ∈{0.1,0.3,0.5,0.7,0.9}. For transitivity model, we focus on another 625 settings with a ∈{5,15,25,35,45}, b ∈{5,15,25,35,45}, ξ_1 = … =  ξ_p ∈{0.1,0.3,0.5,0.7,0.9} and η_1 = … =η_p ∈{0.1,0.3,0.5,0.7,0.9}. 



For each setting, we present the  time series plots of simulated  {D_t}_t=2^500, {D_1,t}_t=2^500, {D_0,t}_t=2^500, {D̅_t}_t=2^500, {D̅_1,t}_t=2^500 and {D̅_0,t}_t=2^500, defined as in Section <ref>, when p = 50. Figures <ref>–<ref> and Figures <ref>–<ref> report the dynamic patterns of the density-dependent model and transitivity models, respectively. Several conclusions can be drawn. First,  all simulated series {D_t}_t≥2, {D_1,t}_t≥2 and {D_0,t}_t≥2 exhibit stationary behavior. Similar patterns are discussed in Section <ref>. Second,  in contrast to many exponential random graph models (ERGMs), the proposed autoregressive (AR) model design is free from degeneracy problem when α_i,j^t-1 and β_i,j^t-1 are strictly between 0 and 1.
Note that we do observe instances of nearly complete or empty graphs when ξ_i = 0.1 or η_i = 0.1. Such degenerate graphs fall outside the scope of our discussion as the node-specific parameters are set to values as low as 0.01. Examples of this can be seen in Figures <ref>–<ref> and <ref>–<ref>. 


































































































 §.§ Stationarity and ergodicity
 

As stated in the last  paragraph of Section <ref>, for each fixed constant p, {_t}_t≥1 defined by the transitivity model (<ref>) is stationary and ergodic as long as α_i,j^t-1 and β_i,j^t-1 are strictly between 0 and 1. Nevertheless it is a Markov chain with 2^p(p-1)/2 states. When p is a fixed constant, the ergodicity of _t 
(i.e. the average in time converges to the average over the state space) may take a long time to be observed; see Theorem 4.1 in Chapter 3 of <cit.>. However the ergodicity of some scalar summary statistics of _t can be observed in much short time spans, as indicated in the simulation reported below.

We  consider
the following three network density measures at each time t:

    D_t= ∑_i,j: i < j X_i,j^t  p(p-1)/2 ,   
    D_1,t = ∑_i,j: i < j (1-X_i,j^t-1) X_i,j^t/p(p-1)/2 ,
       
    D_0,t = ∑_i,j: i< j X_i,j^t-1 (1-X_i,j^t)/p(p-1)/2 ,

where D_t is the network density at time t, and D_1,t and D_0,t are, 
respectively, the densities of newly formed edges and newly dissolved edges at time t. If {_t}_t≥1 is stationary, all three density sequences {D_t}_t≥1, {D_1,t}_t≥2 and {D_0,t}_t≥2 are also  stationary. We also plot 

    D̅_t = 1  t∑_u=1^t D_u ,     D̅_1,t = 1  t∑_u=1^t D_1,u ,     D̅_0,t = 1  t∑_u=1^t D_0,u,

against t for t≥2, to see how quickly the ergodicity can be observed. These are sample means of one-dimensional network summaries. We expect that their convergences are much faster than that for the sample mean of p× p network _t itself.


Setting ξ_1=⋯=ξ_p and η_1=⋯=η_p, we let
(ξ_i, η_i, a, b) take four different sets of values: 
(0.7,0.8,30,15), (0.6,0.7,20,20), (0.6,0.7,15,10) and (0.6,0.7,10,10).

Figure <ref> displays the time series plots of simulated  {D_t}_t=2^200, {D_1,t}_t=2^200, {D_0,t}_t=2^200, {D̅_t}_t=2^200, {D̅_1,t}_t=2^200 and {D̅_0,t}_t=2^200 when p=50. 
As expected, all simulated series {D_t}_t≥2, {D_1,t}_t≥2 and {D_0,t}_t≥2 exhibit patterns in line with stationarity. The convergence of their sample means is observed with the sample sizes greater than 50. 


In particular, (ξ_i, η_i, a, b) = (0.7,0.8,30,15) displays the most dynamic edge changing behaviour, while (ξ_i, η_i, a, b) =(0.6,0.7,20,20) is the least dynamic among the four settings.  









 §.§ A more general model



As stated towards the end of Section <ref>,
 a more general transitivity model admits the form:

    α_i,j^t-1() = ξ_i ξ_je^a_1 U_i,j^t-1/1+ e^a_1 U_i,j^t-1+ e^b_1 V_i,j^t-1 ,
            β_i,j^t-1() =η_i η_j e^b_2 V_i,j^t-1/1+ e^a_2 U_i,j^t-1+ e^b_2 V_i,j^t-1 ,

with =(a_1, b_1, a_2, b_2,ξ_1,…,ξ_p, η_1,…,η_p)^∈ℝ_+^2p+4 >0. Different from the transitivity model (<ref>) introduced in Section <ref>, we allow a_1≠ a_2 and b_1≠ b_2 in (<ref>). We adopt the same simulation settings as above, i.e. 
(ξ_i,  a_1, b_1, η_i, a_2, b_2) ∈{(0.7,30,15, 0.8,30,15), (0.6,20,20,0.7,20,20), (0.6,15,10,0.7,15,10), (0.6,10,10,0.7,10,10)}. 




The two estimation methods are implemented in the same manner as in Section <ref>.
Table <ref> reports the resulting rMAEs over 400 replications with p = 50 and n ∈{100,200,400}. 
The estimation for the parameters in  α_i,j^t-1(), namely ξ_i, a_1 and b_1, exhibits the similar patterns as in Table <ref>, where (ξ_i,  a_1, b_1, η_i, a_2, b_2)=(0.7,30,15, 0.8,30,15) achieves the best estimation accuracy. In contrast, the estimation for parameters in β_i,j^t-1() deteriorates significantly,  and especially for a_2 and b_2.

Note that only some components of _t-1 with X_i,j^t-1=1, t∈[n]∖[m] and j≠ i,  were used in estimating parameters a_2 and b_2. For sparse networks, the total number of those data points is small. This is the intrinsic difficulty in estimating the parameters in β_i,j^t-1().
See also the relevant discussion at the end of Section <ref>.






























§ CONFERENCE INTERACTIONS



We apply our model to an additional dynamic network dataset, in this case of face-to-face interactions among attendees of an academic conference.
The conference in question was the 2009 congress of the Société Française d'Hygiène Hospitalière (SFHH) <cit.>.
The original data was collected automatically by RFID badges worn by the conference participants. 
We analyze a subset of the data corresponding to an active portion of the first day of the congress (June 4, 2009) from about 11:00AM to 6:00PM among the p=200  most active participants out of the total 403.
Each of the n=22 network snapshots corresponds to a non-overlapping time window, with X_i,j^t=1 if participants were in close proximity at any time during the prior 20 minutes.



Similar to Section <ref>, we summarize some key features of this dataset.
Figure <ref> (left panel) shows that while there are some spikes in edge density, there is no clear increasing or decreasing pattern, so we choose to model this dataset with a single AR network model.
Figures <ref> and <ref> show empirical evidence of temporal edge dependence, as well as transitivity effects: after accounting for edge density, edges persist at a higher rate than they grow, they more often grow for node pairs which had more common neighbours, and they more often dissolve for node pairs which had more disjoint neighbours.







Fitting our AR network model with transitivity, we estimate â=26.75 and b̂=15.14, confirming these empirical dynamic effects of common and disjoint neighbours. 
We summarize the estimates of the local parameters {ξ_i}_i=1^200 and {η_i}_i=1^200 in Figure <ref>.



The estimates {ξ̂_i}_i=1^200 have mean 0.18 and a longer right tail, while the estimates {η̂_i}_i=1^200 have mean 1.16 and a longer left tail. 
Moreover, their scatter plot shows that there is a negative relationship between these estimates for a given node. All of these observations are consistent with overall degree heterogeneity of the network.
The estimates {ξ̂_i}_i=1^200 have a wide range from 0.05 to 0.64, while the estimates {η̂_i}_i=1^200 are all between 0.97 and 1.24.
This implies that conference attendees are more heterogeneous in their propensity to form new connections, than in their propensity to extend the length of existing ones.








Finally, we compare our model to the same competing models described in Section <ref> in terms of AIC and BIC. These results are reported in Table <ref>. Our AR network model with transitivity achieves the smallest AIC, followed closely by the global AR model, which only has 2 parameters.
The global AR model achieves the smallest BIC, thus in both cases the best model incorporates temporal edge dependence.
Under either criterion, the two edgewise models require O(p^2) parameters and thus perform poorly, as there are relatively many nodes (p=200) compared to network samples (n=22).

Although the transitivity model has O(p) parameters, it achieves the smallest AIC and the 2nd smallest BIC,  suggesting that degree heterogeneity parameters and the imposed transitivity form are an effective parameterization to summarize the structure in this dataset.

