Weighted stochastic block model

We propose a weighted stochastic block model (WSBM) which extends the stochastic block model to the important case in which edges are weighted. We address the parameter estimation of the WSBM by use of maximum likelihood and variational approaches, and establish the consistency of these estimators. The problem of choosing the number of classes in a WSBM is addressed. The proposed model is applied to simulated data and an illustrative data set.


Introduction
Networks are used in many scientific disciplines to represent interactions among objects of interest. For example, in the social sciences, a network typically represents social ties between actors. In biological sciences, a network can represent interactions between proteins. The stochastic block model (SBM) (Holland et al. 1983;Snijders and Nowicki 1997) is a popular generative model which partitions vertices into latent classes. Conditional on the latent class allocations, the connection probability between two vertices depends only on the latent classes to which the two vertices belong. Many extensions of the class SBM have been proposed which include the degree correlated SBM (Karrer and Newman 2011;Peng and Carvalho 2016), mixed membership SBM (Airoldi et al. 2008) and overlapping SBM (Latouche et al. 2011).
The SBM and many of its variants are usually restricted to Bernoulli networks. However, many binary networks are produced after applying a threshold to a weighted relationship (Ghasemian et al. 2016) which results in the loss of potentially valuable information. Although most of the literature has focused on binary networks, there is a growing interest in weighted graphs (Barrat et al. 2004;Newman 2004;Peixoto 2018).
In particular, a number of clustering methods have been proposed for weighted graphs including algorithm based and model based methods. Algorithm based methods for clustering of weighted graphs can be further divided into two classes: algorithms which do not explicitly optimize any criteria (Pons and Latapy 2005;von Luxburg 2007) and those directly optimize a criterion (Clauset et al. 2004;Stouffer and Bascompte 2011). Model based methods (Mariadassou et al. 2010;Aicher et al. 2013Aicher et al. , 2015Ludkin 2020) attempt to take into account the random variability in the data. A recent review of graph clustering methods is given by (Leger et al. 2014). Mariadassou et al. (2010) presents a Poisson mixture random graph model for integer valued networks and proposes a variational inference approach for parameter estimation. The model can account for covariates via a regression model. In Zanghi et al. (2010), a mixture modelling framework is considered for random graphs with discrete or continuous edges. In particular, the edge distribution is assumed to follow an exponential family distribution. Aicher et al. (2013) proposed a general class of weighted stochastic block model for dense graphs where edge weights are assumed to be generated according to an exponential family distribution. In particular, their construction produces complete graphs, in which every pair of vertices is connected by some real-valued weight. Since most realworld networks are sparse, the constructed model cannot be applied directly. To address this shortcoming, Aicher et al. (2015) extends the work of Aicher et al. (2013) and models the edge existence using a Bernoulli distribution and the edge weights using an exponential family distribution. The contributions of edgeexistence distribution and edge-weight distribution in the likelihood function are then combined via a simple tuning parameter. However, their construction does not result in a generative model and it is not obvious how to simulate network observations from the proposed model. More recently, Ludkin (2020) presents a generalization of the SBM which allows artbitrary edge weight distributions and proposes a reversible jump Markov chain Monte Carlo sampler for estimating the parameters and the number of blocks. However, the use of continuous probability distribution to model the edge weights implies that the resulting graph is complete whereby every edge is present. This assumption is unrealistic for many applications whereby a certain proportion of the real-valued edges is 0. Haj et al. (2020) presents a binomial SBM for weighted graphs and proposes a variational expectation maximization algorithm for parameter estimation.
In this paper, we propose a weighted Stochastic Block model (WSBM) with gamma weights which aims to capture the information of weights directly using a generative model. Both maximum likelihood estimation and variational methods are considered for parameter estimation where consistency results are derived. We also address the problem of choosing the number of classes using the Integrated Completed Likelihood (ICL) criteria (Biernacki et al. 2000). The proposed models and inference methodology are applied to an illustrative data set.

Model specification
In this section, we present the weighted stochastic block model in detail and introduce the main notations and assumptions.
We let X ¼ ðV; X ; YÞ denote the set of directed weighted random graphs where V ¼ N is the set of countable vertices, X ¼ f0; 1g NÂN is the set of edge-existence adjacency matrix, and Y ¼ R NÂN þ is the set of weighted adjacency matrix. Given a random adjacency matrix X ¼ fX ij g i;j2N , X ij ¼ 1 if an edge exists from vertex i to vertex j and X ij ¼ 0 otherwise. The associated weighted random adjacency matrix is given by: Let P be a probability measure on X.

Generative model
We now describe the procedure of generating a sample of random graph (V, X, Y) with n vertices from X.
• Let p ¼ ðp ql Þ Q q;l¼1 be a Q Â Q matrix with entries in [0, 1]. Conditional on the block allocations Z ½n , the entries X ij for i 6 ¼ j of the edge-existence adjacency matrix X ½n is generated from independently a Bernoulli distribution • Let a ¼ ða ql Þ Q q;l¼1 and b ¼ ðb ql Þ Q q;l¼1 be Q Â Q matrices with entries taking values in the positive reals. Conditional on the latent block allocations Z ½n and edge-existence adjacency matrix X, the weighted adjacency matrix Y ½n is generated independently from where Ga ðÁ; ÁÞ denotes the gamma distribution and d fÁg is the Dirac delta function.
The generative framework described above is a straightforward extension of the binary stochastic block model whereby a positive weight is generated according to a gamma distribution for each edge. In particular, (X, Z) is a realization of the binary directed SBM. The gamma distribution is chosen due to its flexibility in the sense that, depending on the value of its shape parameter, it can represent distributions of different shapes.
The log-likelihood of the observations X ½n and Y ½n is given by L 2 ðY ½n ; X ½n ; h; p; a; bÞ ¼ log X z ½n e L 1 ðY ½n ;X ½n ;z ½n ;p;a;bÞ PfZ ½n ¼ z ½n g ; where the sum is over all possible latent block allocations, PfZ ½n ¼ z ½n g ¼ Q n i¼1 h z i is the probability of latent block allocation z ½n , and L 1 ðY ½n ; X ½n ; z ½n ; p; a; bÞ is the complete data log-likelihood, where f ðÁ; a; bÞ is the gamma probability density function with shape parameter a and rate parameter b,

Assumptions
We present several assumptions needed for identifiability and consistency of maximum likelihood estimates. The following four assumptions were presented in Celisse et al. (2012) and are needed in this paper.
(A1) requires that no two classes have the same connectivity probabilities. If this assumption is violated, the resulting model has too many classes and is nonidentifiable. (A2) requires that the connectivity probability between any two classes strictly lies within a closed subset of the unit interval. Note that this assumption is slightly more restrictive compared to assumption 2 of Celisse et al. (2012) in that we do not consider the boundary cases where p q;l 2 f0; 1g. The boundary cases require special treatment and are not pursued in this paper. (A3) ensures that no class is empty with high probability while (A4) is the empirical version of (A3). We note that (A4) is satisfied asymptotically under the generative framework in Sect. 2.1 since the block allocations are generated according to a multinomial distribution.
In addition to the four assumptions above, we also have the following constraints on the gamma parameters.
The log-likelihood function (1) contains degeneracies that prevent the direct estimation of parameters h; p; a; b. To see this, we note that the probability density function of a gamma distribution Ga ða; bÞ is given by f ðy; a; bÞ ¼ y aÀ1 expðÀbyÞb a CðaÞ .
By Stirling's formula, we have Therefore, letting a ! 1 while keeping ab ¼ y, we have f ðy; a; bÞ ! 1. One can therefore show that the log-likelihood function is unbounded above. To avoid likelihood degeneracy, we compactify the parameter space. That is, we restrict the parameter space to a compact subset which contains the true paraemters. Therefore, we have the following assumption.
Assumption 6 There exists 0\a c \a C \1 and 0\b c \b C \1 such that for all ðq; lÞ 2 f1; . . .; Qg, With this assumption, it is easy to see that the log-likelihood function is bounded for any sample size.

Identifiability
Sufficient conditions for identifiability of binary SBM with two classes have been first obtained by Allman et al. (2009). Celisse et al. (2012) show that the SBM parameters are identifiable up to a permutation of class labels under the conditions that ph has distinct coordinates and n ! 2Q. The condition on ph is mild since the set of vectors violating this assumption has Lebesgue measure 0. The identifiability of weighted SBM is more challenging where the only known result (Section 4 of Allman et al. (2011)) requires all entries of ðp; a; bÞ to be distinct. We note that the assumptions in the previous section are not necessarily sufficient but are necessary to ensure that the identifiability of the parameters.

Asymptotic recovery of class labels
We study the posterior probability distribution of the class labels Z ½n given the random adjacency matrix X ½n and weight matrix Y ½n , which is denoted by PðZ ½n jX ½n ; Y ½n Þ. Since X ½n and Y ½n are random, PðZ ½n jX ½n ; Y ½n Þ is also random. Let P Ã ðX ½n ; Y ½n Þ :¼ PðX ½n ; Y ½n jZ ½n ¼ z Ã ½n Þ be the true conditional distribution of ðX ½n ; Y ½n Þ which depends on the true parameters ðh Ã ; p Ã ; a Ã ; b Ã Þ. We study the convergence rate of PðZ ½n jX ½n ; Y ½n Þ towards 1 with respect to P Ã .
The matrices p; a; b are permutation-invariant if one permutes both its rows and columns according to some permutation r : f1; . . .; Qg ! f1; . . .; Qg. Let p r ; a r ; b r be the matrices defined by p r ql ¼ p rðqÞ;rðlÞ , a r ql ¼ a rðqÞ;rðlÞ , b r ql ¼ b rðqÞ;rðlÞ and define the set Two vectors of class labels z and z 0 are equivalent if there exists r 2 R such that denote the equivalence class of z and will omit the square-brackets in the equivalence class notation as long as no confusion arises.
The following result extends Theorem 3.1 of Celisse et al. (2012) to the case of WSBM.

Maximum likelihood estimation of WSBM parameters
For the binary SBM, consistency of parameter estimation have been shown for profile likelihood maximization (Bickel and Chen 2009), spectral clustering method (Rohe et al. 2011), method of moments approach (Bickel et al. 2011), method based on empirical degrees (Channarond et al. 2012), and others (Choi et al. 2012). Consistency of both maximum likelihood estimation and variational approximation method are established in Celisse et al. (2012) and Bickel et al. (2013) where asymptotic normality is also established in Bickel et al. (2013). Abbe (2018) reviews recent development in the stochastic block model and community detections.
Ambroise and Matias (2012) proposes a general class of sparse and weighted SBM where the edge distribution may exhibit any parametric form and studies the consistency and convergence rates of various estimators considered in their paper. However, their model requires the edge existence parameter to be constant across the graph, that is, p ql ¼ p; for all q, l, or p ql can be modelled as p ql ¼ a 1 I q¼l þ a 2 I q6 ¼l where I Á is the indicator function and a 1 6 ¼ a 2 . Furthermore, they also assume that conditional on the block assignments, the edge weight Y ij jX i ¼ q; X j ¼ l is modelled using a parametric distribution with a single parameter h ql . They further impose the restriction that These assumptions are more restrictive than those imposed in this paper. Jog and Loh (2015) studies the problem of characterizing the boundary between success and failure of MLE when edge weights are drawn from discrete distributions. More recently, Brault et al. (2020) studies the consistency and asymptotic normality of the MLE and variational estimators for the latent block model which is a generalization of the SBM. However, the model considered in Brault et al. (2020) is restricted to the dense setting and requires the observations in the data matrix to be modelled by univariate exponential family distributions. This section addresses the consistency of the MLE of WSBM. In particular, we extend the results obtained in the pioneering paper of Celisse et al. (2012) to the case of weighted graphs. Our proof closely follows the proof of consistency of the MLE in Celisse et al. (2012). The MLE consistency proof of ðp; a; bÞ and h require different treatments since there are nðn À 1Þ edges but only n vertices. The following result established the MLE consistency of ðp; a; bÞ.
Theorem 3 Let ðĥ;p;â;bÞ denote the MLE of ðh Ã ; p Ã ; a Ã ; b Ã Þ and assume that for any metric d in R Q .

Variational estimators
Direct maximization of the log-likelihood function is intractable except for very small graphs since it involves a sum over Q n terms. In practice, approximate algorithms such as Markov Chain Monte Carlo (MCMC) and variational inference algorithms are often used for parameter inference. For the SBM, both MCMC and variational inference approaches have been proposed (Snijders and Nowicki 1997;Daudin et al. 2008). Variational inference algorithms have also been developed for mixed membership SBM (Airoldi et al. 2008), overlapping SBM (Latouche et al. 2011), and the weighted SBM proposed in Aicher et al. (2013). This section develops a variational inference algorithm for the WSBM which can be considered a natural extension of the algorithm proposed in Daudin et al. (2008) for the SBM. The variational method consists in approximating PðZ ½n ¼ ÁjX ½n ; Y ½n Þ by a product of n multinomial distributions. Let D n denote a set of product multinomial distributions For any D s ½n 2 D n , the variational log-likelihood is defined by J ðY ½n ; X ½n ; s ½n ; h; p; a; bÞ ¼ L 2 ðY ½n ; X ½n ; h; p; a; bÞ À KLðD s ½n ; PðÁjX ½n ; Y ½n ÞÞ.
Here KLðÁ; ÁÞ denotes the Kullback-Leibler divergence between two probability distributions, which is nonnegative. Therefore, J provides a lower bound on the log-likelihood function. We have that where b ðÁ; pÞ denotes the probability mass function of a Bernoulli distribution with parameter p, and recall that f ðÁ; a ql ; b ql Þ denotes the density function of a gamma distribution Ga ða ql ; b ql Þ. The variational algorithm works by iteratively maximizing the lower bound J with respect to the approximating distribution D s ½n , and estimating the model parameters. Maximization of J with respect to D s ½n consists of solvinĝ s ½n :¼ argmax s ½n J ðY ½n ; X ½n ; s ½n ; h; p; a; bÞ, where h; p; a; b can be replaced by plug-in estimates. This has a closed-form solution given byŝ Conditional onŝ ½n , the variational estimators of ðh; p; a; bÞ are found by solving, ðh;p;ã;bÞ ¼ argmax h;p;a;b J ðY ½n ; X ½n ; s ½n ; h; p; a; bÞ: Closed-form updates forh andp exist and are given bỹ On the other hand, updates forã andb do not have a closed form since the maximum likelihood estimators of the two parameters of a gamma distribution do not have closed forms. However, using the fact that a gamma distribution is a special case of a generalized gamma distribution, Ye and Chen (2017) derived simple closed-form estimators for the two parameters of gamma distribution. The estimators were shown to be strongly consistent and asymptotically normal. For q; l ¼ 1; . . .; Q, let us define the quantities then the updates for a ql ; b ql are given bỹ We obtain the variational estimators ðh;p;ã;bÞ by computing (3), (4), (5), (6), (7) until convergence. We now address the consistency of the variational estimators derived above. The following two propositions are the counterpart of Theorem 2 and 3 for variational estimators. We omit the proof since they follow similar arguments as the proof of Proposition 1 Assume that assumptions (A1), (A2), (A3), (A5) and (A6), and let ðh;p;ã;bÞ be the variational estimators defined above. Then for any distance dðÁ; ÁÞ on ðp; a; bÞ, dððp;ã;bÞ; ðp Ã ; a Ã ; b Ã ÞÞ À! ½ n ! 1P0: Proposition 2 Assume that the variational estimators ðp;ã;bÞ converge at rate 1/n to ðp Ã ; a Ã ; b Ã Þ, respectively, and assumptions (A1), (A2), (A3), (A5), (A6) hold. We have dðh; h Ã Þ À! ½ n ! 1P0; where d denotes any distance between vectors in R Q .
Note a stronger assumption on the convergence rate 1/n ofp;ã;b is assumed for Proposition 2 compared to ffiffiffiffiffiffiffiffiffi ffi log n p =n in Theorem 3. The same assumption is also used in Theorem 4.4. of Celisse et al. (2012).

Choosing the number of classes
In real world applications, the number of classes is typically unknown and needs to be estimated from the data. For the SBM, a number of methods have been developed to determine the number of classes, including log-likelihood ratio statistic (Wang and Bickel 2017), composite likelihood (Saldaña et al. 2017), exact integrated complete data likelihood (Côme and Latouche 2015) and Bayesian framework (Yan 2016) based methods. Model selection for variants of SBM have also been investigated (Latouche et al. 2014).
We apply the integrated classification likelihood (ICL) criterion developed by Biernacki et al. (2000) to choose the number of classes for the WSBM. The ICL is an approximation of the complete-date integrated likelihood. The ICL criterion for SBM have been derived by Daudin et al. (2008) under the assumptions that the prior distribution of ðh; pÞ factorizes and a non-informative Dirichlet prior on h. Here we follow the approach of Daudin et al. (2008) to derive an approximate ICL for the WSBM.
Let m Q denote the model with Q blocks, the ICL criterion is an approximation of the complete-data integrated likelihood: LðY ½n ; X ½n ; Z ½n ; h; p; a; b; m Q Þgðh; p; a; bÞdhdpdadb where gðh; p; a; bÞ is the prior distribution of the parameters. Assuming a noninformative Jeffreys prior on h, a Stirling approximation to the gamma function, and finally a BIC approximation to the conditional log-likelihood function, the approximate ICL can be derived. For a model m Q with Q blocks, and assuming a non-informative Jeffreys prior Dirð0:5; . . .; 0:5Þ on h, the approximate ICL criterion is: log LðY ½n ; X ½n ;Z ½n ; h; p; a; b; m Q Þ À 3 2 ðQðQ þ 1ÞÞ log nðn À 1Þ À Q À 1 2 log n whereZ ½n is the estimate of Z ½n . The derivation of above follows exactly the same lines as the proof of Proposition 8 of Daudin et al. (2008).

Simulation
We validate the theoretical results developed in previous sections by conducting simulation studies. In particular, we investigate how fast parameter estimates of WSBM converge to their true values and the accuracy of posterior block allocations. Additionally, we investigate the performance of ICL in choosing the number of blocks.

Experiment 1 (two-class model)
For each fixed number of nodes n, 50 realizations of WSBM are generated based on the fixed parameter setting given in (8 and 9). The variational inference algorithm derived in Sect. 5 is then applied to estimate the model parameters and class allocations. We can see from Table 1 that the estimated model parameters converge to their true values as the number of nodes increases while the posterior class assignment is accurate across any number of nodes. Table 2 shows the ICL criterion tends to select the correct number of classes, especially when the number of nodes is large.
We can see from Table 3 that the estimated parameters converge to their true values quickly as the number of nodes increases. The ICL criterion tends to overestimate the number of classes when the number of nodes is small, but consistently selects the correct model when the number of nodes is large (Table 4).

Computational complexity
The computational complexity of the variational estimators derived in Sect. 5 scales as Oðn 2 Þ, which has the same complexity as the variational algorithm developed by Daudin et al. (2008). Therefore, the algorithm may be prohibitively expensive for networks with more than 1000 nodes. The estimated computational time for the two-class model in Sect. 7.1 and for the three-class model in Sect. 7.2 for various values of n are shown in Fig. 1. The estimated computing time at each n is the average running time of the variational algorithm over 20 replications.

Application: Washington bike data set
We apply the WSBM to analyse the Washington bike sharing scheme data set. 1 Information with respect to start stations and end stations of trips as well as length (travel time) of trips are available in the data set. We select a time window of one week staring from January 10th, 2016 and construct the adjacency matrix X and weight matrix Y as follows: • X ij ¼ 1 if there is trip starting from station i and finishing at station j.
• Y ij is the total length of trips (in minutes) from station i to station j.
The resulting network consists of 370 nodes with an average out-degree of 36.14. The average total length of trips between any pair of stations is 42.15 minutes. We apply the ICL criterion to select the number of classes for the WSBM. For each number of classes Q, the variational inference algorithm is fitted to the network 20 times with 20 random initializations and the highest value of ICL is recorded, and the six-class model is chosen (Table 5). Each bike station is plotted on the map in   Fig. 2 where its colour represents the estimated class assignment. We observe that bike stations in class 6 (colored in brown) tend to be concentrated in the central area of Washington wheraas stations in class 3 (colored in red) tend to be located further from the center. Figure 2 shows some spatial effect in the class assignment of bike stations whereby stations that are close in distance tend to be in the same cluster, with the exception of class 3 (colored in red). One potential extension of the model is to take into account the spatial locations of the bike stations as covariates.
The estimated class proportionsĥ shown in 12 indicate that class 3 has the largest number of stations whereas class 4 has the smallest number of stations. The estimatedp shows that within class connectivity is generally higher compared to between class connectivity. We further observe that the connection probabilities  between bike stations in class 1, 4, and 6 are substantially higher. Interestingly, we observe a near symmetry in the matrixp indicating that the probability of having a trip from a station in class k to another station in class l is similar with the probability of having a trip from class l to class k. The estimated densities of travel time between each pair of classes are shown in Fig. 3. We observe that the majority of the estimated densities have mode and mean close to 0, particularly for the estimated densities in the diagonal of Fig. 3. This implies that the total travel times between stations in the same class are quite short. In comparison, the total travel time between stations in different classes tend to be longer. This is reasonable as the distance between bike stations in different classes tend to be longer which in turn requires longer travel time.

Discussion
This paper proposes a weighted stochastic block model (WSBM) for networks. The proposed model is an extension of the stochastic block model. A variational inference strategy is developed for parameter estimation. Asymptotic properties of maximum likelihood estimators and variational estimators are derived, and the problem of choosing the number of classes is addressed by using an ICL criteria. Simulation studies are conducted to evaluate the performance of variational estimators and the use of ICL to determine the number of classes. The proposed model and inference methods are an illustrative data set. It is straightforward to extend the WSBM to allow node covariates. Let w i 2 R d be the covariates for each node i ¼ 1; . . .; n, and let w ij be the covariates for each pair of nodes, i; j ¼ 1; . . .; n; i 6 ¼ j. The edge probability p ij between a pair of nodes i, j with node i in block q and node j in block l can be modelled as log p ij 1 À p ij ¼ n ql;0 þ n T ql;1 w ij þ n T ql;2 w i þ n T ql;3 w j ; where n ql;0 2 R and n ql;1 ; n ql;2 ; n ql;3 2 R d . Conditional on block assignments and existence of an edge between a pair of nodes i, j, Y ij can be modelled as a Gamma random variable with mean l ij and variance r ij where where / ql;0 2 R and / ql;1 ; / ql;2 ; / ql;3 2 R d . Many possible future extensions are possible. First, it is desirable to investigate further theoretical properties of maximum likelihood and variational estimators of WSBM parameters such as asymptotic normality of the estimators. Furthermore, some of the assumptions imposed in this work in order to ensure consistency of the estimators maybe relaxed. Moreover, the number of blocks is assumed to be fixed in the asymptotic analysis of the estimators. It would be interesting to allow the number of blocks to grows as the number of nodes grows.

Auxiliary results
Definition 1 A random variable X with mean l ¼ EðXÞ is sub-exponential if there are non-negative parameters ðm; bÞ such that Eðe tðXÀlÞ Þ e m 2 t 2 2 for all jtj\ 1 b : The following lemma is a straight forward consequence of the definition.
Lemma 1 If the independent random variables fX i g n i¼1 are sub-exponential with parameters ðm i ; b i Þ, i ¼ 1; . . .; n, then P n i¼1 X i is sub-exponential with parameters ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P n i¼1 m 2 i p ; max n i¼1 b i Þ. The following results show that for a Gamma random variable Y and a Bernoulli random variable X, both XY and X log Y are sub-exponential random variables. They are useful for the proof of the main theorems.
Proof The expectation of YX is given by and thus The moment generating function of YX is given by Àa which is defined for t\b. Taylor series expansion of ð1 À t b Þ Àa around t ¼ 0 gives that Similary, Taylor series expansion of exp Àt a b p À Á around 0 gives that exp Àt Thus, we have This leads to Hence, we can choose suitable m;b such that Eðe tðYXÀlÞ Þ exp m 2 t 2 2 for all jtj\ 1 b : h Proposition 4 If Y $ Gaða; bÞ and X $ BerðpÞ, X log Y is a sub-exponential random variable.
Proof The proof is analogous to the proof of Proposition 3. It is straightforward to show that Eðe tX log Y Þ ¼ ð1 À pÞ þ p Cða þ tÞ CðaÞb t and l :¼ EðX log YÞ ¼ pðWðaÞ À logðbÞÞ where WðÁÞ is the digamma function. The rest of the proof follows the same line as the proof of Proposition 3 by considering the Taylor series expansion of e Àlt and Eðe tX log Y Þ around t ¼ 0. h The following lemma provides an upper bound for the tail probability of a subexponential random variable.
Lemma 2 (Sub-exponential tail bound)Suppose that X with mean l is subexponential with parameters ðm; bÞ. Then i .
The following inequality for suprema of random processes is needed for the proof of Theorem 2.

Proof
In this section, we prove Theorems 1, 2, 3 for the special case a q;l ¼ a 11 for all ðq; lÞ 2 f1; . . .; Qg. The more general case can be proved similary by using the fact that X ij log Y ij is a sub-exponential random variable.

Proof of Theorem 1
Proof Our proof is adapted from Celisse et al. (2012). Using the fact that for every z 0 ½n 2 ½z ½n , Pðz 0 ½n jX ½n ; Y ½n Þ ¼ Pðz ½n jX ½n ; Y ½n Þ, we have that where jjz ½n À z Ã ½n jj 0 is the number of difference between z ½n and z Ã ½n . To simplify notation, we write z :¼ z ½n , X :¼ X ½n , and Y ¼ Y ½n . We have that Consider the first term on the RHS of the last inequality, log PðYjX; zÞ The summand in the summation above vanishes when b Ã For two vectors z and z 0 , we define for all i, j. An application of Hoeffding's inequality yields that for some constant L 1 [ 0.
By Proposition 3, the random variable Y ij X ij is sub-exponential and the tail bound for sub-exponential random variables (Lemma 2) implies that for some constants L 2 ; L 3 [ 0. Proposition B.4 of Celisse et al. (2012) shows that N r ðzÞ is bounded below by N r ðzÞ ! c 2 2 njjz ½n À z Ã ½n jj 0 .
Combining inequalities (14), (15)  where c is some negative constant. Therefore, the first condition of Proposition 5 is satisfied. Fix some b ð0Þ , (A6) implies that there exist u 1 ¼ OðnÞ and u 2 ¼ Oð1Þ such that jjb À b ð0Þ jj 2 u 1 , and jjb À b ð0Þ jj 2 u 2 . Therefore the second condition of Proposition 5 is also satisfied with D ¼ Q 2 . Now, we have Since P i6 ¼j X ij Y ij b ð0Þ z i ;z j is a sub-exponential random variable with parameters ðnm; bÞ, the second term on the RHS of inequality 18 can be bounded by To bound the first term on the RHS of inequality 18, we introduce the set Pz ½n ¼ fðp; a; bÞ : ðz ½n ; p; a; bÞ 2 Pg for every z ½n , and define the event X n ðz ½n Þ ¼ n sup Pðz ½n Þ q n jS n ðbÞ À S n ðb ð0Þ Þj q n j ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 2 ðQ 2 þ x n Þ p þ q n bðQ 2 þ x n Þ o .

Proof of Thereom 3
Proof We denotePðZ ½n ¼ z ½n jX ½n ; Y ½n Þ as the conditional distribution of Z ½n under the parameters ðĥ;p;â;bÞ. The following result is an extension of Proposition 3.8 of Celisse et al. (2012) and is needed to establish the consistency ofĥ.
The proof of Proposition 3.8 of Celisse et al. (2012) shows that P Ã ½fT 0 [ À 10r log ng \ X n g C 1 & exp 8n log n À C 2 ðlog nÞ 2 nv 2 ÀKL G ðða Ã ql ; b Ã ql Þjjða Ã q 0 ;l 0 b Ã q 0 ;l 0 ÞÞ ¼ ÀjD Ã jK Ã Since X ij Y ij is a sub-exponential random variable and X ij is bounded for all i, j, one can show that Using similar techniques as in the proof of Proposition 3.8 of Celisse et al. (2012), one can show that P Ã ½fT 1;4 [ tg \ X n X k X D;jDj¼k Q 2 exp À C 4 t 2 v 2 n ðk þ jD Ã jÞ , ð28Þ i . Since ð1 þ ðQ À 1Þu 0 n Þ n ! 1 as n ! 1, and the proof is completed. h The proof of the consistency ofĥ is a consequence of the proposition above and follows the same lines as the proof of Theorem 3.9 of Celisse et al. (2012). h Funding Open Access funding provided by the IReL Consortium. The funding was provided by Science Foundation Ireland (Grant No. SFI/12/RC/2289-2).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/.