Maud Delattre

Valentine Genon-Catalot

Adeline Samson

ESTIMATION OF POPULATION PARAMETERS IN STOCHASTIC DIFFERENTIAL EQUATIONS WITH RANDOM EFFECTS IN THE DIFFUSION COEFFICIENT

Keywords: 1991 Mathematics Subject Classification. 62F10; 62F12. Approximate maximum likelihood estimator, asymptotic normality, consistency, estimating equations, random effects models, stochastic differential equations

come

Estimation of population parameters in stochastic differential equations with random effects in the diffusion coefficient

Maud Delattre, Valentine Genon-Catalot,

Introduction

Stochastic differential equations (SDEs) with random effects have been the subject of several recent contributions, with various applications such as pharmacokinetic/pharmacodynamic, neuronal modeling [START_REF] Delattre | Coupling the SAEM algorithm and the extended kalman filter for maximum likelihood estimation in mixed-effects diffusion models[END_REF][START_REF] Donnet | A review on estimation of stochastic differential equations for pharmacokinetic -pharmacodynamic models[END_REF]Picchini et al., 2010). Several estimation methods have been proposed to provide estimators in these complex models. For general mixed effects models (not only SDE), the main reference for a theoretical study of the exact maximum likelihood is (Nie, 2006). But the likelihood has no explicit expression except in some special cases and assumptions are then very difficult to verify. In Delattre et al. (2013), the case of a linear random effect in the drift together with a specific distribution for the random effects is investigated. In this case, the exact maximum likelihood estimator is explicit and studied. In the general case, Picchini et al. (2010); Picchini and Ditlevsen (2011) propose approximations of the likelihood based on Hermite expansion and Gaussian quadrature. All these references work with random effects in the drift, and not in the diffusion coefficient, except [START_REF] Delattre | Coupling the SAEM algorithm and the extended kalman filter for maximum likelihood estimation in mixed-effects diffusion models[END_REF] who incorporate measurement error and propose an approximation of the likelihood with the extended Kalman filter.

Here, we focus on discretely observed SDEs with a random effect in the diffusion coefficient. The distribution of the random effect depends on unknown parameters to be estimated. Applications of these models have been considered in population pharmacokinetics or pharmacodynamic where the diffusion coefficient represents the intra-individual variability, which may vary from one individual to another, thus being modeled as a random effect [START_REF] Delattre | Coupling the SAEM algorithm and the extended kalman filter for maximum likelihood estimation in mixed-effects diffusion models[END_REF]. We assume that there is a linear random effect in the diffusion coefficient. For simplicity, we assume that the drift is zero in Sections 2 and 3, and give the extension to the case of non null drift in Section 4.

Statistical inference for discretely observed SDEs with no random effects has been widely studied (see Kessler et al., 2012, and references therein). In Genon-Catalot and Jacod (1993) the estimation of unknown fixed parameters in the diffusion coefficient is studied with discrete observations of a single trajectory when the sampling interval tends to zero. The likelihood of these observations is not explicit, therefore estimating equations are built based on the Euler approximation of the SDE with drift set to zero. One of the strategies described below follows the same idea, but here the parameters are random. This complicates the definition and the theoretical study of the estimator.

More precisely, we consider N real valued stochastic processes (X i (t), t ≥ 0), i = 1, . . . , N , with dynamics ruled by the following SDEs:

d X i (t) = φ i σ(X i (t)) dW i (t), X i (0) = x 0 i , i = 1, . . . , N, (1)
where (W i ) 1≤i≤N are N independent Wiener processes, (φ i ) 1≤i≤N are N i.i.d. random variables taking values in (0, +∞), (φ i ) 1≤i≤N and (W i ) 1≤i≤N are independent. The function σ(x) is known and real-valued. Each process (X i (t)) represents an individual, the variable φ i represents the random effect of individual i. The variables (φ i ) 1≤i≤N have a common distribution g(ϕ, θ)dν(ϕ) on (0, +∞) where ν is a dominating measure and θ is a vector of unknown parameters called population parameters, belonging to a set Θ ⊂ R p . Our aim is to estimate θ from discrete observations {X i (t j ), j = 1, . . . , n, i = 1, . . . , N } where t j = jT /n and T is fixed. In the case of a linear random effect in the diffusion coefficient (1), choosing an inverse Gamma distribution leads to explicit estimators. Therefore, we consider the specific case

φ i = 1 Γ 1/2 i
with Γ i ∼ G(a, λ), a > 0, λ > 0, θ = (a, λ).

(2)

We study the exact maximum likelihood estimator in the case σ(.) ≡ 1. When σ(.) ≡ 1, we build estimating equations based on the Euler approximation of the fixed effect diffusion model. The difficulty of these estimating equations is that the Euler approximation has to be integrated out with respect to the distribution of the random effects. Moreover, we build another type of estimating equations, corresponding to the ideal likelihood of directly observed random effects where estimators of the random effects are plugged in. This second approach has the advantage to be easily generalized to any distribution for the random effects.

The paper is organized as follows. Section 2 introduces some assumptions and gives the exact likelihood and its approximation obtained by Euler scheme. Our asymptotic framework is when the number N of subjects tends to infinity. In Section 3, we study the asymptotic properties of the estimators. When σ(.) ≡ 1, the exact maximum likelihood estimator of θ is asymptotically Gaussian with rate √ N both for fixed number of measurements per subject n and for n tending to infinity. When σ(.) ≡ 1, we must assume that n depends on N and satisfies the constraint N/n → 0 for the first method, √ N /n → 0 for the second. Our estimators are asymptotically Gaussian with rate √ N . Section 4 deals with a non nul drift. Simulations illustrate the behavior of the estimators and results are presented in Section 5. Section 6 concludes the paper with some extensions. Proofs are gathered in Appendix.

Exact and approximate likelihoods

Consider N real valued stochastic processes (X i (t), t ≥ 0), i = 1, . . . , N , with dynamics ruled by (1). The processes (W i ) 1≤i≤N and the r.v.'s (φ i ) 1≤i≤N are defined on a common probability space (Ω, F, P). Consider the filtration (F t , t ≥ 0) defined by

F t = σ(φ i , W i (s), s ≤ t, i = 1, . . . , N ). As F t = σ(W i (s), s ≤ t) ∨ F i t , with F i t = σ(φ i , φ j , W j (s), s ≤ t, j = i) independent of W i , each process W i is a (F t , t ≥
0)-Brownian motion. Moreover, the random variables φ i are F 0 -measurable. In what follows, we assume that (H1) The function σ belongs to C 2 (R) and for all

x ∈ R, 0 < σ 2 0 ≤ σ 2 (x) ≤ σ 2 1 , |σ ′ (x)| + |σ ′′ (x)| ≤ K.
Under (H1), the process (X i (t)) is well-defined and (φ i , X i (t)) is strong Markov adapted to the filtration (F t , t ≥ 0). The N processes (φ i , X i (.)) 1≤i≤N are independent. For all ϕ, and all

x 0 i ∈ R, the fixed effect SDE d X ϕ,x 0 i i (t) = ϕ σ(X ϕ,x 0 i i (t)) dW i (t), X ϕ,x 0 i i (0) = x 0 i (3)
admits a unique strong solution process (X ϕ,x 0 i i (t), t ≥ 0) adapted to the filtration (F t , t ≥ 0). From the Markov property of (φ i , X i (t)), we deduce that the conditional distribution of X i given φ i = ϕ is identical to the distribution of X ϕ,x 0 i i (for more details, see Delattre et al., 2013). For i = 1, . . . , N , the process (X i (t), t ∈ [0, T ]) is discretely observed at times t j = jT /n, j = 0, . . . , n and we set for fixed

T . ∆ = T n , X i = (X i (t j ), j = 1, . . . , n), i = 1, . . . , N. (4)
We start with the exact likelihood of (4). The distribution of the observations

(X i ) 1≤i≤N on N i=1 R n has the form P θ = ⊗ N i=1 P i θ where P i θ is the distribution of X i on R n . If Q i ϕ,x 0 i denotes the distribution of X ϕ i = (X ϕ,x 0 i i (t j ), j = 1, . . . , n
) and p t (x, y, ϕ) the transition density of (3), then Q i ϕ,x 0 i admits the density n j=1 p ∆ (x i,j-1 , x i,j , ϕ) w.r.t. the Lebesgue measure of R n (with x i,0 = x 0 i ). Therefore, the density of P i θ w.r.t. the Lebesgue measure of R n is given by:

λ i (θ, x i ) = +∞ 0 g(ϕ, θ) n j=1 p ∆ (x i,j-1 , x i,j , ϕ)dν(ϕ),
where

x i = (x i,j , j = 1, . . . , n). The exact likelihood is Λ N (θ) = N i=1 λ i (θ, x i ).
Here, we are faced with two problems. First, the transition densities of (3) are generally not explicit. Second, and this is specific to random effect models (SDE or more generally non Gaussian models), even if these transition densities were explicit, it is generally not possible to get a closed-form expression for the marginal density of X i , which corresponds to the integral λ i (θ, x i ). Therefore the exact likelihood is not explicit and difficult to study theoretically and numerically.

Instead of using the exact transition densities of (3), it is now standard to use the approximation given by the transition densities of the corresponding Euler scheme, i.e. the one-step discretisation of (3) (see e.g. [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF]Genon-Catalot and Jacod, 1993;Kessler et al., 2012), although for small n, it may give a bias. Therefore, we introduce

Li (X i , ϕ) = 1 ϕ n n j=1 σ(X i (t j-1 )) exp (- S i 2ϕ 2 ) ∝ 1 ϕ n exp (- S i 2ϕ 2 ),
with, for i = 1, . . . , N ,

S i = 1 ∆ n j=1 (X i (t j ) -X i (t j-1 )) 2 σ 2 (X i (t j-1 )) (5)
Note that when σ(.) ≡ 1, Li (X i , ϕ) is the exact density of (X ϕ,x 0 i i (t j ), j = 1, . . . , n) (see ( 3)). To estimate θ, instead of the exact likelihood, we introduce the approximate likelihood, corresponding to the Euler scheme integrated with respect to the random effects distribution:

ΛN (θ) = N i=1 +∞ 0 ϕ -n exp - S i 2ϕ 2 g(ϕ, θ)dν(ϕ). (6)
When σ(.) ≡ 1, the problem reduces to

S i /n = φ 2 i 1 n n j=1 (Wi(tj )-Wi(tj-1)) 2 ∆ and 1 n∆ n j=1 (W i (t j )-W i (t j-1
)) 2 has a distribution χ 2 (n)/n, independent of φ i . Thus the model can be seen as an observation of φ i with a multiplicative error. In both cases of σ(.) ≡ 1 and σ(.) ≡ 1, a theoretical study of the estimators based on ΛN (θ) could be possible using the approach developed by Nie ( 2006) but his assumptions are generally difficult to verify. Below, as in Delattre et al. (2013), we rather introduce a specific distribution for the random effects allowing to obtain an explicit formula for (6). In Section 3, we are able to directly study the corresponding estimators.

Remark 2.1. Except in the case σ(.) ≡ 1 where (6) is the exact likelihood, our approach based on an approximate likelihood imposes a double asymptotic framework where both N and n tend to infinity. As n → ∞, note that the statistic S i based on the i-th trajectory provides an estimator of the random effect φ i . Indeed, let

M i (t) = t 0 σ(X i (s))dW i (s), and
R i = n j=1 (M i (t j ) -M i (t j-1 )) 2 /σ 2 (X i (t j-1 )). (7)
By standard properties of quadratic variations, R i /T → 1 in probability as n → ∞. Thus, using X i (t j ) = X i (t j-1 ) + tj tj-1 φ i σ(X i (s))dW i (s), S i /n = φ 2 i R i /T tends to φ 2 i .

A specific distribution for the random effect

For a general distribution g(ϕ, θ)dν(ϕ) of the random effect φ i , the integral in (6) has no explicit expression. However, for the conjugate distribution, namely the inverse Gamma (2), an explicit expression is obtained. The unknown parameter is then θ = (a, λ) ∈ Θ = R + × R + . The true value is denoted by θ 0 .

Let us start with the ideal case of directly observed random effects φ i (or Γ i ). Then, the exact log-likelihood of (Γ 1 , . . . , Γ N ) is given by:

ℓ N (θ) = N a log λ -N log Γ(a) + (a -1) N i=1 log Γ i -λ N i=1 Γ i (8) with associated score function S N (θ) = ∂ ∂λ ℓ N (θ) ∂ ∂a ℓ N (θ) ′ where ∂ ∂λ ℓ N (θ) = N i=1 a λ -Γ i , ∂ ∂a ℓ N (θ) = N i=1 (-ψ(a) + log λ + log Γ i ) ,
where ψ(z) = Γ ′ (z) Γ(z) is the di-gamma function. By standard properties of Gamma distributions, we have, under the true value θ 0 , (1/ √ N )S N (θ 0 ) → D N 2 (0, I(θ 0 )), where I(θ) is

I(θ) = a λ 2 -1 λ -1 λ ψ ′ (a) . (9)
Note that using properties of the di-gamma function (see Section 8), I(θ) is invertible for all θ ∈ (0, +∞) 2 . The maximum likelihood estimator based on the observation of Γ 1 , . . . , Γ N , denoted θ N = θ N (Γ 1 , . . . , Γ N ) is consistent and satisfies √ N (θ Nθ 0 ) → D N 2 (0, I -1 (θ 0 )) as N tends to infinity. But the Γ i 's are not observed. Two different strategies are studied. Following Remark 2.1, a natural idea consists in plugging in ℓ N (θ) the estimator n/S i of Γ i . This reveals to be more complex than expected (Section 3.2) and we will need to truncate the estimator n/S i . The other strategy (Section 3.1) is based on (6). We provide asymptotic results when n is fixed and N → ∞ in the case σ(.) ≡ 1, and when both n, N → ∞ for a general σ(.).

Estimation based on the Euler approximation of the likelihood.

Let LN (θ) = log ΛN (θ) be the log Euler contrast (see ( 5) and ( 6)).

Proposition 3.1. Under (H1) and ( 2), we have:

LN (θ) = N i=1 log Γ(a + n/2) Γ(a) + aN log λ - N i=1 (a + n/2) log (λ + 1 2 S i ). ( 10
)
The associated gradient vector (pseudo-score function)

GN (θ) = ∂ ∂λ LN (θ) ∂ ∂a LN (θ) ′ (11) is given by ∂ ∂λ LN (θ) = N i=1 a λ -a+n/2 λ+Si/2 and ∂ ∂a LN (θ) = N i=1 (ψ(a + n/2) -ψ(a)) + N i=1 log λ λ+Si/2
. For the Hessian matrix (pseudo Fisher information matrix)

I N (θ) = - ∂ 2 ∂λ 2 LN (θ) ∂ 2 ∂λ∂a LN (θ) ∂ 2 ∂λ∂a LN (θ) ∂ 2 ∂a 2 LN (θ) , ( 12
) we get ∂ 2 ∂λ 2 LN (θ) = - N i=1 a λ 2 -a+n/2 (λ+Si/2) 2 , ∂ 2 ∂λ∂a LN (θ) = N i=1 1 λ -1 λ+Si/2 and ∂ 2 ∂a 2 LN (θ) = -N (ψ ′ (a) -ψ ′ (a + n/2)).
We study the estimators defined by the estimating equation:

GN ( θN ) = 0. ( 13
)
We consider two asymptotics: n fixed (section 3.1.1) and n → ∞ (section 3.1.2).

Fixed number of observations per subject

We assume that the number n of observations per subject is fixed and that the number of subjects N tends to infinity. The only model that can be studied with this asymptotic is the special case σ(.) ≡ 1. We denote by the upper index 1 all the quantities associated to this model:

dX 1 i (t) = φ i dW i (t)
, and the statistic is

S 1 i = n j=1 (X 1 i (t j ) -X 1 i (t j-1 )) 2 /∆. ( 14
)
The distribution of S 1 i can be explicitly computed. Proposition 3.2. Under P θ , the random variables

β 1 i (λ) = λ λ+S 1
i /2 , i = 1, . . . , N, are independent and β 1 i (λ) has distribution beta of the first kind on (0, 1) with parameters (a, n/2). The random variables S 1 i /(2λ) are independent with distribution on (0, +∞) beta of the second kind with parameters (n/2, a).

Then, LN (θ) = L1

N (θ) where S i is replaced by S 1 i is the exact log-likelihood. Define the associated exact maximum likelihood estimator as any solution of:

θ1 N = Argsup θ L1 N (θ). ( 15
)
Proposition 3.3. Assume that n is fixed. Then, the maximum likelihood estimator θ1 N (15) is consistent. Let

I n (θ) = a(n/2) λ 2 (a+1+n/2) -n/2 λ(a+n/2) -n/2 λ(a+n/2) ψ ′ (a) -ψ ′ (a + n/2) (16)
Then, the matrix I n (θ 0 ) is invertible and under

P θ0 √ N ( θ1 N -θ 0 ) → D N 2 (0, I -1 n (θ 0 )). Remark that I n (θ) = I(θ) + O( 1 n ).
The proof of Proposition 3.3 is standard and omitted. It is simply the asymptotic study of the maximum likelihood estimator based on the i.i.d. sample (S 1 i , i = 1, . . . , N ) whose distribution is specified by Proposition 3.2.

Number of observations per subject goes to infinity

Now, we assume that both n and N tend to infinity for some fixed T . The strategy consists in studying the case σ(.) ≡ 1 where computations can all be done explicitely and then studying the difference between the general case and the case σ(.) ≡ 1. Some preliminary results on the moments of R i /T are needed (Proposition 3.4). For these results specifically, we do not assume that the φ 2 i 's have inverse Gamma distribution. As already said before, S i /n = φ 2 i R i /T tends to φ 2 i as n tends to infinity (see ( 7) for the definition of R i ). Let us define the equivalent of R i for the model σ(.) ≡ 1:

R 1 i = n j=1 (W i (t j ) -W i (t j-1
)) 2 . We denote by P θ = ⊗ i≥1 P i θ the distribution of the sequence of processes (φ i , (X i (t), t ∈ [0, T ])), i ≥ 1 on (0, +∞) × i≥1 C([0, T ]), by E θ the corresponding expectation. Note that P θ is the marginal distribution of (X i , i = 1, . . . , N ) under P θ . Both R i /T and R 1 i /T tend to 1 in probability as n → ∞. Furthermore:

Proposition 3.4. Under (H1), for all θ, we have

E θ R 1 i T -1|φ i = 0, |E θ Ri T -1|φ i | ≤ C T n φ 2 i , and |E θ Ri T - R 1 i T |φ i | ≤ C T n φ 2 i , and for all p ≥ 1, E θ ( Ri T -1) 2p |φ i ≤ C( T p n p + T 2p n 2p (φ 2p i + φ 4p i )) and E θ ( Ri T - R 1 i T ) 2p |φ i ≤ C T 2p n 2p (φ 2p i + φ 4p i ).
We now study the score function ( 11) and the Fisher information matrix (12).

Proposition 3.5. Recall GN (θ 0 ) defined by ( 11) and I(θ 0 ) given in (9). For σ(.) ≡ 1, as N, n tend to infinity, under P θ0 , GN (θ 0 )/ √ N converges in distribution to N 2 (0, I(θ 0 )). In the general case, if E θ0 φ 8 i < +∞, i.e. if a 0 > 4 and N, n tend to infinity in such a way that N/n tends to 0, the same result holds.

The convergence of the Fisher information matrix is as follows: Proposition 3.6. In the case σ(.) ≡ 1 and the general case, the pseudo Fisher information matrix given in (12), I N (θ 0 )/N , converges in probability to I(θ 0 ) (see ( 9)) as N, n tend to infinity, under P θ0 . Now we study the estimator θN defined by (13).

Proposition 3.7. Assume that n, N → +∞ in such a way that N/n tends to 0. Then, an estimator θN which solves (13) exists with probability tending to one as N tends to infinity under P θ0 and is weakly consistent. The matrix I(θ 0 ) is invertible and under P θ0 , √ N ( θNθ 0 ) → D N 2 (0, I -1 (θ 0 )). Moreover, the estimator θN is asymptotically equivalent to the MLE θ N = θ N (Γ 1 , . . . , Γ N ) based on the direct observation of (Γ 1 , . . . , Γ N ).

If we discretely observe N trajectories of a fixed effect SDE, i.e. dX i (t) = ϕσ(X i (t))dW i (t), the rate of convergence for ϕ is √ N n with the same constraint N/n = o(1). This is because all the N discrete trajectories are used to estimate a unique parameter ϕ. In the case of a random effect SDE, the rate of convergence of θ is √ N , which is the same as the one obtained when the random effects Γ i are directly observed.

3.2. Approach based on estimators of the random effects.

In this section, we exploit directly the fact that the random effect φ 2 i = Γ -1 i can be estimated using the vector X i = (X i (t j ), j ≤ n) by S i /n. The idea is simply to replace the random variables Γ i by their estimator n/S i in the likelihood (8) of (Γ 1 , . . . , Γ N ). But this works only when σ(.) ≡ 1. More precisely, in this case, let us set:

U N (θ) = N a log λ -N log Γ(a) + (a -1) N i=1 log (n/S 1 i ) -λ N i=1 (n/S 1 i ).
and consider the estimators θ * N given by

∇U N (θ * N ) = 0. (17)
Otherwise, we must define estimators of log Γ i and Γ i which have appropriate moment properties. These properties are obtained by associating a truncation to the estimators log (n/S i ) and n/S i . Let us define:

logΓ i = log (n/S i )1 (Si/n≥k/ √ n) , Γ i = (n/S i )1 (Si/n≥k/ √ n) .
where k is a constant. Note that, by the above definitions, logΓ i and Γ i are set to 0 outside the set

(S i /n ≥ k/ √ n)
where S i /n is not bounded from below. The following holds.

Lemma 3.8. Assume that E θ0 φ 8 i < +∞, i.e. a 0 > 4 and that n > 4. Then,

E θ0 n S i 1 (Si/n≥k/ √ n) - n S 1 i 2 + E θ0 log S i n 1 (Si/n≥k/ √ n) -log S 1 i n 2 ≤ C n 2
Then we consider the function

V N (θ) = N a log λ -N log Γ(a) + (a -1) N i=1 logΓ i -λ N i=1 Γ i , (18)
and the associated estimator θ * * N defined by the estimating equation:

∇V N (θ * * N ) = 0. ( 19
)
Proposition 3.9. Assume that σ(.) ≡ 1. If N, n tend to infinity in such a way that √ N /n tends to 0, then an estimator θ * N which solves (17) exists with probability tending to 1 under P θ0 and is weakly consistent. Moreover, √ N (θ * Nθ 0 ) converges in distribution to N 2 (0, I -1 (θ 0 )) and θ * N is asymptotically equivalent to the exact MLE θ N associated to (Γ 1 , . . . , Γ N ), i.e.

√ N (θ * Nθ N ) = o P θ 0 (1). When σ(.) is not equal to 1, the same result holds for θ * * N under the condition E θ0 φ 8 i < +∞, i.e. a 0 > 4 .

Note that in this approach, even when σ(.) ≡ 1, the constraint √ N /n → 0 is required.

Extension to nonnul drift and unbounded coefficients.

The estimation results obtained above can be extended to more general SDEs with a drift term and unbounded function σ(.). Indeed, assume that the observed processes are given by:

d X i (t) = b(X i (t))dt + φ i σ(X i (t)) dW i (t), X i (0) = x 0 i , i = 1, . . . , N,
where the drift function b : R → R may be known or unknown. We define, for i = 1, . . . , N , S i with the same formula ( 5)

S i = 1 ∆ n j=1 (X i (t j ) -X i (t j-1 )) 2 σ 2 (X i (t j-1 )
) .

We consider the same estimation procedures using S i . As all proofs rely on Proposition 3.4, we only need to extend this proposition to φ -2 i S i /n instead of R i /T. This amounts to studying the properties of the difference

φ -2 i S i /n -R i /T or equivalently S i /n -φ 2 i R i /T
. This is given by the next proposition.

Proposition 4.1. Assume that the drift function b is C 1 (R), that σ is C 2 (R) with |b ′ | + |σ ′ | + |σ ′′ | ≤ K and σ 2 (x) ≥ σ 2 0 > 0.
Then, for all θ, we have

|E θ Si n -φ 2 i Ri T |φ i | ≤ Cφ 2 i T
n , and for all p ≥ 1,

E θ S i n -φ 2 i R i T 2p |φ i ≤ C T p n p (φ 2p i + 1) + T 2p n 2p φ 4p i
This proposition can be extended to more general drifts, especially to drifts including random effect.

Numerical simulation results.

We compare the performances of both estimation methods on simulated data for several models. Two sets of population parameters θ 0 are used: (a 0 = 6, λ 0 = 1) and (a 0 = 5, λ 0 = 3). In each case, 100 datasets are generated with an Euler scheme with sampling interval δ = 10 -4 T on time interval [0, T ], with T = 5, and N = 50, 100 subjects, n = 100, 500. The parameter θ 0 is estimated via θN (method 1) and via either θ * N or θ * * N (method 2). The empirical mean and standard deviation are computed from the 100 datasets. We consider:

Example 1. dX i (t) = φ i dW i (t) , X i (0) = 0. Example 2. dX i (t) = φ i 1 + X 2 i (t)dW i (t) , X i (0) = 0. Example 3. dX i (t) = -ρX i (t)dt + φ i dW i (t) , X i (0) = 0. Example 4. dX i (t) = -ρX i (t)dt + φ i 1 + X 2 i (t)dW i (t) , X i (0) = 0.
For Example 1, estimation method 1 leads to the exact MLE of θ 0 . Examples 2, 3, 4 rely on Section 4. Both methods require to optimize a criterion. The two criteria use the function gamma, which takes large values. Thus, the optimization is implemented vanishing the gradients of the two criteria, which are more stable.

The results for Examples 1 to 4 are displayed in Tables 1 to 4 respectively. The results are satisfactory overall and similar for the 4 models, even when the model includes a drift. Method 1 estimators are biased for n = 100. When σ(.) ≡ 1, this is expected due to the Euler approximation of the likelihood. Nevertheless, for fixed N , we observe the convergence of the estimators to the true value when n increases. We also observe the convergence of the estimators to the true values when both N and n increase. This clearly illustrates consistency of the estimators when both n and N tend to infinity. Finally, the implementation of the method 2 in Examples 2, 3, 4 requires to choose a value for the threshold k. The results are displayed for k = 0.5. Simulations with various values of k have not shown any significant impact of k on the estimators performances.

N = 50

N = 100 n = 100 n = 500 n = 100 n = 500 (a 0 = 5, λ 0 = 3) ã 5. 16 (1.33) Table 2. Example 2: σ 2 (x) = 1 + x 2 . Empirical mean and standard deviation (in brackets) of θ (method 1) and θ * * (method 2) computed from 100 datasets.

Extensions and concluding remarks.

In this paper, we study the estimation of population parameters in a SDE with a linear random effect in the diffusion coefficient from discrete observations of N i.i.d. trajectories on a fixed length time interval. We especially study the case of a null drift and of φ i = 1/Γ 1/2 i with Γ i ∼ G(a, λ). This leads to estimators using two different approaches. The first method is based on an approximation of the exact likelihood relying on the Euler scheme of the SDE. The second method uses a plug-in of estimators of the random effects in the likelihood of (φ 1 , . . . , φ N ).

Several extensions are possible. The theoretical results extend easily to N individual processes observed on [0, T i ], i = 1, . . . , N at n discrete time points t i,j = jT i /n with all T i 's fixed and T i ≤ T for some fixed T . Another direction for extensions is to look at other distributions for the random effects. In particular, the plug-in method applies for any distribution provided that we introduce appropriate truncations as is done here. The Euler approximation could also be applied with other distributions but numerical integration would then N = 50 N = 100 n = 100 n = 500 n = 100 n = 500 (a 0 = 5, λ 0 = 3) ã 5. 16 2.54 (0.39) 2.90 (0.44) (a 0 = 6, λ 0 = 1) ã 6.21 (1.28) 6.17 (1.15) 6.01 (0.81) 6.01 (0.76) λ 1.01 (0.22) 1.02 (0.20) 0.98 (0.15) 1.00 (0.14) a * * 5.48 (1.17) 6.02 (1.10) 5.38 (0.66) 5.88 (0.73) λ * * 0.88 (0.20) 1.00 (0.19) 0.86 (0.12) 0.98 (0.13) Table 4. Example 4. (ρ = 1) Empirical mean and standard deviation (in brackets) of θ (method 1) and θ * * (method 2) computed from 100 datasets. be necessary. For a more general model for the diffusion coefficient including a non linear random effect, the two approaches studied here could be extended, in particular the second method.

Genon-Catalot, V. and Jacod, J. (1993). On the estimation of the diffusion coefficient for multi-dimensional diffusion processes. Ann. Inst. H. Poincaré Probab. Statist., 29(1):119-151. Genon-Catalot, V., Jeantheau, T., and Laredo, C. (1999). Parameter estimation for discretely observed stochastic volatility models. Bernoulli, 5:855-872. Kessler, M., Lindner, A., and Sorensen, M., editors (2012)

Appendix: proofs

Proof of Proposition 3.1 Using the fact that φ -2 i has Gamma distribution G(a, λ), we get the result as:

λi (θ, X i ) = (0,+∞) λ a γ a-1+n/2 Γ(a) exp [-γ(λ + 1 2 S i )] dγ = λ a Γ(a + n/2) Γ(a) λ + 1 2 S i a+n/2 . Proof of Proposition 3.2 Let χ i = n j=1 (W i (t j ) -W i (t j-1 )) 2 /∆ = R 1 i /∆. As S 1 i = Γ -1 i χ i , (see 14) , β 1 i (λ) = λΓi λΓi+ 1 2 χi , S 1 i 2λ = χi/2
λΓi Under P θ , Γ i and χ i are independent, Γ i is G(a, λ) and χ i is χ 2 (n) = G(n/2, 1/2). Hence the results using Proposition 8.1.

Proof of Proposition 3.4 We need the following Lemma and Proposition:

Lemma 7.1. For all θ, E θ ((X i (t) -X i (s)) 2p |φ i ) ≤ C(2p)σ 2p 1 φ 2p i |t -s| p where C(2p) is a numerical constant. Proposition 7.2. (R 1 i /T ) -1 = T -1 T 0 H n i,1 (s)dW i (s) and (R i /T ) -1 = T -1 ( T 0 H n i (s)dW i (s) + T 0 K n i (s)dW i (s) + T 0 L n i (s)ds),
where, for j = 1, . . . , n and s ∈]t j-1 , t j ],

H n i,1 (s) = 2(W i (s) -W i (t j-1 )), H n i (s) = 2 (M i (s) -M i (t j-1 ))σ(X i (s)) σ 2 (X i (t j-1 )) , K n i (s) = 2φ i (t j -s) σ 2 (X i (s))σ ′ (X i (s)) σ 2 (X i (t j-1 )) , L n i (s) = φ 2 i κ(X i (s))(t j -s) σ 2 (X i (s)) σ 2 (X i (t j-1 ))
, with κ = σσ ′′ + (σ ′ ) 2 .

Lemma 7.1 and Proposition 7.2 yield E θ (

R 1 i T -1|F 0 ) = 0 and E θ ( R i T -1|F 0 ) = E θ ( R i T - R 1 i T |F 0 ) = 1 T T 0 E θ (L n i (s)|F 0 )ds.
Using (H1), we get

|L n i (s)| ≤ Cφ 2 i n j=1 1 ]tj-1,tj ] (s)(t j -s) ≤ Cφ 2 i ∆1 ]0,T ] (s)
, for C depending on σ 0 , σ 1 , K. Thus, the first inequality of Proposition 3.4.

As (

Ri T -1) 2p = (A 1 + A 2 + A 3 ) 2p ≤ 3 2p-1 3 i=1 A 2p
i , we study separately the three terms A 2p i . We have

A 2p 3 = 1 T T 0 L n i (s)ds 2p ≤ (Cφ 2 i ∆) 2p
. Next, we use the Burkholder-Davies-Gundy (BDG), the Hölder inequalities and (H1):

T 2p E θ (A 2p 2 |F 0 ) ≤ C(2p)E θ ( T 0 (K n i (s)) 2 ds p |F 0 ) ≤ C(2p)T p-1 E θ ( T 0 E θ ((K n i (s)) 2p |F 0 )ds)
where

(K n i (s)) 2p ≤ Cφ 2p i ∆ 2p 1 ]0,T ] (s). Finally, for T 2p E θ (A 2p 1 |F 0 ) we study T 0 (H n i (s)) 2 ds p
. By the Holder inequality, we have,

T 0 (H n i (s)) 2 ds p ≤ C p n p-1 n j=1 ∆ p-1 tj tj-1 (M i (s) -M i (t j-1 )) 2p ds.
Consequently, for constants C depending on σ 0 , σ 1 , K,

E θ ( T 0 (H n i (s)) 2 ds p |F 0 ) ≤ C n j=1 tj tj-1 E θ ( s tj-1 σ 2 (X i (u))du p |F 0 )ds ≤ C∆ p .
Finally, to study the difference R i -R 1 i , we only need to study the term:

T 0 (H n i (s) -H n i,1 (s))dW i (s) = T 0 2 3 k=1 n j=1 1 ]tj-1,tj ] (s)Z k i,j (s)dW i (s) with Z 1 i,j (s) = σ(X i (s)) -σ(X i (t j-1 )) σ 2 (X i (t j-1 )) s tj-1 (σ(X i (u)) -σ(X i (t j-1 )))dW i (u), Z 2 i,j (s) = σ(X i (s)) -σ(X i (t j-1 )) σ(X i (t j-1 )) s tj-1 dW i (u), Z 3 i,j (s) = 1 σ(X i (t j-1 )) s tj-1 (σ(X i (u)) -σ(X i (t j-1 )))dW i (u),
These terms are studied analogously using the BDG and Cauchy-Schwarz inequalities.

Proof of Lemma 7.1 Recall that φ i is F 0 -measurable and when dealing with the process (X i (t)), conditioning on φ i is equal to conditioning on F 0 . We have

E θ ((X i (t) -X i (s)) 2p |F 0 ) = φ 2p i E θ ((M i (t) -M i (s)) 2p |F 0 )
. By the BDG inequality and (H1), for s ≤ t,

E θ ((M i (t) -M i (s)) 2p |F 0 ) ≤ C(2p)E ( t s σ 2 (X i (u))du) p |F 0 ≤ C(2p)σ 2p 1 (t -s) 2p .
Proof of Proposition 7.2 By the Ito formula, we have:

(M i (t j ) -M i (t j-1 )) 2 = 2 tj tj-1 (M i (s) -M i (t j-1 ))σ(X i (s))dW i (s) + tj tj-1 σ 2 (X i (s))ds.
We split:

σ 2 (X i (s)) = σ 2 (X i (t j-1 ))+σ 2 (X i (s))-σ 2 (X i (t j-1
)) and use the Ito formula:

σ 2 (X i (s))-σ 2 (X i (t j-1 )) = φ i s tj-1 (σ 2 ) ′ (X i (u))σ(X i (u))dW i (u) + 1 2 φ 2 i s tj-1 (σ 2 ) ′′ (X i (u))σ 2 (X i (u))du.
Integrating, Fubini formula yields:

tj tj-1 (σ 2 (X i (s)) -σ 2 (X i (t j-1 )))ds = φ i tj tj-1 (t j -u)(σ 2 ) ′ (X i (u))σ(X i (u))dW i (u) + 1 2 φ 2 i tj tj-1 (t j -u)(σ 2 ) ′′ (X i (u))σ 2 (X i (u))du.
The result for R 1 i corresponds to σ(.) ≡ 1. Hence the results.

Proof of Proposition 3.5 Recall that S i = nΓ -1 i R i /T . We have by ( 11) GN (θ

0 ) = S N (θ 0 )+ N i=1 Y i (θ 0 ) N i=1 Z i (θ 0 ) ′ where Y i (θ 0 ) = Γ i - (a 0 + n 2 )Γ i λ 0 Γ i + C i and Z i (θ 0 ) = ψ(a 0 + n 2 ) -log (λ 0 Γ i + C i ).
Therefore, we have to prove that 1

√ N N i=1 Y i (θ 0 ) N i=1 Z i (θ 0 )
′ tends to 0 in P θ0 -probability as n, N tend to infinity. To distinguish the two cases σ(.) ≡ 1 and σ(.) = 1, we introduce the random variables Y 1 i (θ 0 ), Z 1 i (θ 0 ) where we replaced S i by S 1 i . We proceed on two steps:

(1) 1 √ N N i=1 Y 1 i (θ 0 ) N i=1 Z 1 i (θ 0 ) ′ = o P θ 0 (1) as N, n → ∞. (2) 1 √ N N i=1 (Y i (θ 0 ) -Y 1 i (θ 0 )) N i=1 (Z i (θ 0 ) -Z 1 i (θ 0 )) ′ = o P θ 0 (1) as N, n → ∞ under the constraints N/n → 0 and E θ0 φ 8 i < +∞. Proof of (1): Let C 1 i = nR 1 i /(2T ) = Γ i S 1 i /2 and G 0 i = λ 0 Γ i . We have Y 1 i (θ 0 ) = Γ i -a0+n/2 λ0 G 0 i G 0 i +C 1 i .
In what follows, we use repeatedly the fact that G 0 i and

C 1 i are independent, that G 0 i ∼ G(a 0 , 1) and C 1 i ∼ G(n/2, 1). Hence, G 0 i + C 1 i and G 0 i /(G 0 i + C 1 i )
are independent, the latter with distribution β (1) (a 0 , (n/2)), the former with distribution G(a 0 + (n/2), 1).

The r.v.

Y 1 i (θ 0 ), i = 1, . . . , N are i.i.d. with E θ0 (Y 1 i (θ 0 )) = 0, E θ0 (Y 1 i (θ 0 ) 2 ) = a0(a0+1) λ 2 0 (a0+1+n/2) . Therefore, N -1/2 N i=1 Y 1 i (θ 0 ) = o P θ 0 (1). Analogously, Z 1 i (θ 0 ) = ψ(a 0 + n/2) -log (G 0 i + C 1 i ) satisfies E θ0 (Z 1 i (θ 0 )) = 0, E θ0 (Z 1 i (θ 0 )) 2 = 1/(a 0 + n/2) + o(1/n). Thus, N -1/2 N i=1 Z 1 i (θ 0 ) = o P θ 0 (1).
Proof of (2): We introduce

C i = nR i /(2T ). We have Y i (θ 0 ) -Y 1 i (θ 0 ) = (a 0 + n 2 )Γ i G 0 i + C 1 i - (a 0 + n 2 )Γ i G 0 i + C i , and
Z i (θ 0 ) - Z 1 i (θ 0 ) = log (G 0 i + C 1 i ) -log (G 0 i + C i ). Thus, Y i (θ 0 ) -Y 1 i (θ 0 ) = a 0 + n 2 λ 0 G 0 i (G 0 i + C i )(G 0 i + C 1 i ) n 2 ( R i T - R 1 i T ) We introduce the set Ω i = {|(R i /T ) -1| ≤ 1/2}. On Ω i , we use G 0 i + C i ≥ (n/4). So E θ0 |Y i (θ 0 ) -Y 1 i (θ 0 )|1 Ωi ≤ 2 a 0 + n 2 λ 0 E θ0 G 0 i G 0 i + C 1 i | R i T - R 1 i T | .
Then, if E θ0 φ 4 i < +∞, we have by proposition 3.4

E θ0 G 0 i G 0 i + C 1 i 2 = O( 1 n 2 ) and E θ0 R i T - R 1 i T 2 ≤ (T /n) 2 E θ0 φ 4 i .
Therefore, using the Cauchy-Schwarz inequality,

E θ0 |Y i (θ 0 ) -Y 1 i (θ 0 )|1 Ωi ≤ Cn C n 4 1/2 ≤ C n . On Ω c i , we use G 0 i /(G 0 i + C i ) ≤ 1. Therefore, |Y i (θ 0 ) -Y 1 i (θ 0 )|1 Ω c i ≤ n 2 (a 0 + n 2 ) λ 0 1 G 0 i + C 1 i | R i T - R 1 i T |1 Ω c i .
We have if E θ0 φ 8 i < +∞, by Proposition 3.4, E θ0

1 G 0 i +C 1 i 4 = O( 1 n 4 ), E θ0 Ri T - R 1 i T 4 ≤ (T /n) 4 E θ0 (φ 4 i + φ 8 i ), and P θ0 (Ω c i ) ≤ 2 2p E θ0 | Ri T -1| 2p .
Using the Cauchy-Schwarz inequality twice, the above inequality with p = 2 and Proposition 3.4 with the condition E θ0 φ 8 i < +∞, we get:

E θ0 |Y i (θ 0 ) -Y 1 i (θ 0 )|1 Ω c i ≤ Cn 2 (C/n 4 ) 1/4 (C/n 4 ) 1/4 (P θ0 (Ω c i )) 1/2 ≤ C/n.
We can conclude that under the condition E θ0 φ 8 i < +∞,

E θ0 | 1 √ N N i=1 (Y i (θ 0 ) -Y 1 i (θ 0 ))| ≤ C √ N /n. ( 20
)
We now turn to the other difference. We have by the Taylor formula

Z i (θ 0 ) -Z 1 i (θ 0 ) = n 2 1 G 0 i + n/2 ( R 1 i T - R i T ) + n 2 ( R 1 i T - R i T ) 1 0 f i (s)ds,
where

f i (s) = - s n 2 ( R 1 i T -1)+(1-s) n 2 ( R i T -1) (G 0 i +n/2)(G 0 i +s n 2 R 1 i T +(1-s) n 2 R i T )
. We use that:

E θ0 n 2 1 G 0 i + n/2 ( R 1 i T - R 1 i T )|F 0 ≤ n 2 C n 1 G 0 i + n/2 φ 2 i = C 1 Γ 2 i (G 0 i + n/2) and if E θ0 φ 4 i < +∞, E θ0 1 Γ 2 i (G 0 i +n/2) ≤ C( 1 (a0+n/2-1)(a0+n/2-2) ) 1/2 , to obtain that , 1 √ N N i=1 E θ0 E θ0 n 2 1 G 0 i + n/2 ( R 1 i T - R i T )|F 0 ≤ C √ N n
On the other hand, noticing that for s ∈ [0, 1]:

|f i (s)| ≤ n/2 (G 0 i +n/2)G 0 i | R 1 i T -1| + | Ri T -1| , if E θ0 φ 8 i < +∞, E θ0 1 0 f i (s)ds 2 ≤ E θ0 1 (G 0 i ) 4 E θ0 | R 1 i T -1| 4 + | R i T -1| 4 1/2 ≤ C( T n ).
Finally, we have

E θ0 n 2 |( R 1 i T -Ri T ) 1 0 f i (s)ds| ≤ ( C n ) 1/2 . Therefore, E θ0 | 1 √ N N i=1 (Z i (θ 0 )-Z 1 i (θ 0 ))| ≤ C √ N (1/n+ 1/n 1/2 ).
Proof of Proposition 3.6: To obtain I N (θ 0 )/N = I(θ 0 ) + o P θ 0 (1), we have to prove that 1 N N i=1 A i (θ 0 ) → 0, and 1 N N i=1 B i (θ 0 ) → 0, in P θ0probability where

A i (θ 0 ) = a 0 + n/2 (λ 0 + S i /2) 2 = a 0 + n/2 λ 2 0 (G 0 i ) 2 (G 0 i + C i ) 2 , B i (θ 0 ) = 1 λ 0 + S i /2 = Γ i G 0 i + C i . ( 21
)
As in the previous proposition, we separate the cases σ(.) ≡ 1 and σ(.) = 1 and define the random variables A 1 i (θ 0 ), B 1 i (θ 0 ) where S i is replaced by S 1 i :

A 1 i (θ 0 ) = a 0 + n/2 λ 2 0 (G 0 i ) 2 (G 0 i + C 1 i ) 2 , B 1 i (θ 0 ) = Γ i G 0 i + C 1 i .
Recall that C 1 i ∼ G(n/2, 1) and is independent of Γ i . Thus,

E θ0 A 1 i (θ 0 ) = a 0 (a 0 + 1) λ 2 0 (a 0 + n/2 + 1) = O( 1 n ), E θ0 B 1 i (θ 0 ) = a 0 λ 0 (a 0 + n/2) = O( 1 n ).
This implies 1

N N i=1 E θ0 A 1 i (θ 0 ) = O( 1 n ), 1 N N i=1 E θ0 B 1 i (θ 0 ) = O( 1 n ).
Next, we study the differences A i (θ 0 ) -

A 1 i (θ 0 ), B i (θ 0 ) -B 1 i (θ 0 ). A i (θ 0 ) -A 1 i (θ 0 ) = a 0 + n/2 λ 2 0 (C i -C 1 i ) (G 0 i ) 2 (G 0 i + C i )(G 0 i + C 1 i ) 2 + (G 0 i ) 2 (G 0 i + C i ) 2 (G 0 i + C 1 i ) Thus: |A i (θ 0 )-A 1 i (θ 0 )| ≤ a0+n/2 λ 2 0 |C 1 i -C i | 2G 0 i (G 0 i +Ci)(G 0 i +C 1 i ) .
We introduce again the set Ω i = {|(R i /T )-1| ≤ 1/2}. On Ω i , using G 0 i + C i ≥ (n/4), we have

E θ0 |A i (θ 0 ) -A 1 i (θ 0 )|1 Ωi ≤ 2 a 0 + n/2 λ 2 0 (n/4) E θ0 G 0 i G 0 i + C 1 i 2 E θ0 C 1 i -C i 2 1/2 ≤ C/n. Next, using G 0 i /(G 0 i + C i ) ≤ 1 on Ω c i , we have |A i (θ 0 ) -A 1 i (θ 0 )|1 Ω c i ≤ a0+n/2 λ 2 0 |C 1 i -C i | 2 (G 0 i +C 1 i )
. Thus, using the same arguments than in proof of proposition 3.5

E θ0 (|A i (θ 0 ) -A 1 i (θ 0 )|1 Ω c i ) ≤ C (P θ0 (Ω c i )) 1/2 ≤ C n . Analogously, B i (θ 0 ) -B 1 i (θ 0 ) = 1 λ0 (C 1 i -C i ) G 0 i (G 0 i +Ci)(G 0 i +C 1 i ) .
Introducing the set Ω i again, we obtain |B i (θ 0 ) -

B 1 i (θ 0 )| ≤ C G 0 i G 0 i +C 1 i + 1 Ω c i n G 0 i +C 1 i |R i /T -R 1 i /T |. We conclude E θ0 |B i (θ 0 ) -B 1 i (θ 0 )| ≤ C n + (P θ0 (Ω c i )) 1/2 . So
the proof is complete.

Proof of Proposition 3.7: To obtain the weak consistency of θN and its asymptotic normality, we follow the scheme described in Barndorff-Nielsen and Sorensen (1991) (Theorem 3.4 and Lemma 3.5) and Genon-Catalot et al. (1999) (Theorem 4.1), see also Sweeting (1980). We must prove that:

(1) Under P θ0 , GN (θ

0 )/ √ N → D N 2 (0, I(θ 0 )) as N → ∞. (2) I N (θ 0 )/N → I(θ 0 ) in P θ0 -probability. (3) sup θ∈M c,N | I N (θ)/N -I(θ 0 )| → 0 in P θ0 -probability, where M c,N = {θ ∈ (0, +∞) 2 , θ -θ 0 ≤ c/ √ N }. (Uniformity condition)
Points ( 1) and ( 2) are directly implied by Propositions 3.5 and 3.6. It remains to prove (3). We will prove

(a) E θ0 (sup θ∈M c,N | I 1 N (θ)/N -I 1 N (θ 0 )/N |) → 0. (b) E θ0 (sup θ∈M c,N | I 1 N (θ)/N -I N (θ)/N |) → 0. Point (a) Let ε > 0 be such that a 0 -ε > 0, λ 0 -ε > 0. Choose N large enough to ensure that M c,N = {(a, λ) ∈ (0, +∞) 2 , |a -a 0 | ≤ c/ √ N , |λ -λ 0 | ≤ c/ √ N } ⊂ [a 0 -ε, a 0 + c] × [λ 0 -ε, λ 0 + c]
and n > 8. We have

I 1 N (θ)/N -I 1 N (θ 0 )/N = I(θ) -I(θ 0 ) + D 11 N (θ, θ 0 ) D 12 N (θ, θ 0 ) D 12 N (θ, θ 0 ) D 22 N (θ, θ 0 ) where D 11 N (θ, θ 0 ) = 1 N N i=1 (A 1 i (θ) -A 1 i (θ 0 )), D 12 N (θ, θ 0 ) = 1 N N i=1 (B 1 i (θ) -B 1 i (θ 0 )), D 22 N (θ, θ 0 ) = -(ψ ′ (a + n/2)-ψ ′ (a 0 +n/2))
. We only study D 11 and D 12 which are the most difficult. We can write D

11 N (θ, θ 0 ) = c N +d N , with G i = λΓ i and c N = (λ -λ 0 )(a 0 + n/2) 1 N N i=1 Γ i (G i + G 0 i + 2C 1 i ) (G i + C 1 i ) 2 (G 0 i + C 1 i ) 2 , d N = (a 0 -a) 1 N N i=1 Γ i (G i + C 1 i ) 2 .
For θ ∈ M c,N , we have the bounds

|d N | ≤ c (λ0-ε) 2 √ N , and
|c N | ≤ c(a 0 + n/2) √ N (2λ 0 + c) Γ i C 1 i 4 + 2 Γ i C 1 i 3 We have, for n > 2k, E θ0 Γi 2C 1 i k = (2λ 0 ) -k (a0+k-1)(a0+k-2)...a0 (n/2-1)(n/2-2)...(n/2-k) . Thus, E θ0 sup θ∈M c,N |c N | ≤ C √ N n 2 . Thus E θ0 sup θ∈M c,N |D 11 N (θ, θ 0 )| = O(1/ √ N ).
For the other term, (a) is proved as we have

sup θ∈M c,N |D 12 N (θ, θ 0 )| ≤ |λ -λ 0 | 1 λ 0 (λ 0 -ε) = O(1/ √ N ).
We now prove (b). For this, we prove the convergence to 0 of (see ( 21):

E θ0 sup θ∈M c,N |A i (θ) -A 1 i (θ)| and E θ0 sup θ∈M c,N |B i (θ) -B 1 i (θ)|. We have |A i (θ) -A 1 i (θ)| ≤ a0+n/2+c (λ0-ε) 2 |C 1 i -C i | 2Gi (Gi+Ci)(Gi+C 1 i ) . On the set Ω i = {|(R i /T ) -1| ≤ 1/2}, R i /T ≥ 1/2 and G i + C i ≥ n/4. Thus, as G i + C 1 i > C 1 i and G i and C 1 i are independent, for n > 4, E θ0 sup θ∈M c,N |A i (θ) -A 1 i (θ)|1 Ωi ≤ 2 a 0 + n/2 + c (λ 0 -ε) 2 (n/4) (n/2) E θ0 (R 1 i /T -R i /T ) 2 E θ0 ((λ 0 + c)Γ i ) 2 E θ0 (1/C 1 i ) 2 ) 1/2 ≤ C T n E θ0 (φ 2 i + φ 4 i ) 1/2 = O( 1 n ). Next (n > 8, E θ0 φ 8 i < +∞) E θ0 sup θ∈M c,N |A i (θ) -A 1 i (θ)|1 Ω c i ≤ 2(a 0 + n/2 + c)n/2 (λ 0 -ε) 2 (P θ0 (Ω c i )) 1/2 E θ0 (R 1 i /T -R i /T ) 4 E θ0 1 (nR 1 i /T ) 4 1/4 = o(1) We have B i (θ) -B 1 i (θ) = Γi (λΓi+ nR i 2T )(λΓi+ nR 1 i 2T ) ( nR 1 i 2T -nRi 2T ). Using that λ ≥ λ 0 -ε, E θ0 sup θ∈M c,N |B i (θ) -B 1 i (θ)| ≤ C n 2(λ 0 -ε) T n 1 (n/2 -1)(n/2 -2) E θ0 (φ 2 i + φ 4 i ) 1/2 = O( 1 n ).
Therefore, the proof of the first part of Proposition 3.7 is complete. The fact that √ N ( θnθ N ) = o P θ 0 (1) can be deduced from the above proof.

Proof of Proposition 3.9: We first consider the case σ(.) ≡ 1 and the estimating function ∇U N (θ). To get the result, it is enough to prove that:

1 √ N N i=1 ( n S 1 i -Γ i ) = o P θ 0 (1), 1 √ N N i=1 (log n S 1 i -log Γ i ) = o P θ 0 (1), (22)
where we recall that S 1 i /n = Γ -1 i R 1 i /T and nR 1 i /T is independent of Γ i and has distribution χ 2 (n). Using results recalled in Section 8, for n > 2, we have

E θ0 ( n S 1 i -Γ i ) = E θ0 Γ i E θ0 ( n/2 C 1 i -1) = a0 λ0 O(n -1 ). Analogously, for n > 4, E θ0 ( n S 1 i -Γ i ) 2 = E θ0 Γ 2 i E θ0 ( n/2 C 1 i -1) 2 = O(n -1
). This implies: E θ0 ( 1

√ N N i=1 ( n S 1 i -Γ i )) 2 = O(n -1 ) + N -1 2 (O(n -1 )) 2 .
Hence, the first part of ( 22) holds provided that √ N /n = o(1). For the second assertion, we compute E θ0 (log n S 1 i log Γ i ) = -ψ(n/2) + log (n/2) = O(n -1 ), and Var θ0 (log n

S 1 i - log Γ i ) = E θ0 (-log C 1 i + ψ(n/2)) 2 = O(n -1
). Therefore, the second part of ( 22) also holds for √ N /n = o(1).

Next, the result on θ * N will follow analogously from the fact that T N = ( √ N ) -1 N i=1 ( logΓ ilog Γ i ) and τ N = ( √ N ) -1 N i=1 ( Γ i -Γ i ) both tend to 0 in probability. The result follows from the first part of the proof and from Lemma 3.8 which we prove now. Proof of Lemma 3.8: We write: Γ

i -n S 1 i = n Si 1 (Si/n≥k/ √ n) -n S 1 i = ν i + ν ′ i with ν i = ( n S i - n S 1 i ) 1 (Si/n≥k/ √ n) , ν ′ i = - n S 1 i 1 (Si/n<k/ √ n) . (23)
And analogously log Si n 1 (Si/n≥k/ √ n)log

S 1 i n = τ i + τ ′ i with τ i = (log S i n -log S 1 i n ) 1 (Si/n≥k/ √ n) , τ ′ i = -log S 1 i n 1 (Si/n<k/ √ n) . (24)
For n > 4 and E θ0 φ 8 i < +∞, using explicit computations, E θ0 ( n S 1 i ) 2 = O(1), E θ0 log 2 S 1 i n = O(1). To obtain that E θ0 (ν ′ i + τ ′ i ) 2 = O(n -2 ), we now prove that:

P θ0 (S i /n < k/ √ n) ≤ C/n 2 . ( 25
)
Proof of (25): We remark:

( S i n < k √ n ) ⊂ (|φ 2 i - S i n | > φ 2 i - k √ n , φ 2 i ≥ 2 k √ n ) ∪ (φ 2 i < 2 k √ n ) ⊂ (|φ 2 i - S i n | > φ 2 i 2 ) ∪ (φ -2 i > √ n 2k ) = (|1 - R i T | > 1 2 ) ∪ (φ -2 i > √ n 2k
).

Consequently, using the Markov inequality and Proposition 3.4 yields:

P θ0 (S i /n < k/ √ n) ≤ C 2 4 ( T n ) 2 (1 + E θ0 (φ 4 i + φ 8 i ) + (2k) 4 n 2 E θ0 φ -4 i ≤ c ′ /n 2 .
So the proof of (25) is complete.

It remains to study the terms ν i , τ i . We have on (S i /n ≥ k/ √ n):

|ν i | = | 1 S 1 i /n (S 1 i /n -S i /n)( 1 S i /n -Γ i + Γ i )| ≤ 1 R 1 i /T |( R 1 i T - R i T )(1 - R i T )|Γ -1 i √ n k + Γ i R 1 i /T | R 1 i T - R i T |.
We use Proposition 3.4, the Cauchy-Schwarz inequality and the exact distribution of Γ i and 1 R 1 i /T to obtain, for n > 4 and E θ0 φ 8 i < +∞:

E θ0 ν 2 i ≤ C n k 2 E θ0 ( R 1 i T - R i T ) 4 E θ0 (1 - R i T ) 4 1/2 + E θ0 ( R 1 i T - R i T ) 2 ≤ C ′ n 2 .
For the term τ i , we use the Taylor formula and get:

τ i = (S i /n -S 1 i /n) 1 0 ds s(S i /n) + (1 -s)(S 1 i /n) .
Then, we split the integral:

1 0 ds s(S i /n) + (1 -s)(S 1 i /n) = 1 S 1 i /n + 1 0 s(S 1 i /n -S i /n) (s(S i /n) + (1 -s)(S 1 i /n))S 1 i /n ds.
Thus, on (S i /n ≥ k/ √ n), we obtain, after simplifications:

|τ i | ≤ 1 R 1 i /T |R i /T -R 1 i /T | + √ n k 1 R 1 i /T Γ -1 i (R i /T -R 1 i /T ) 2 .
This yields E θ0 τ 2 i ≤ C/n 2 . The proof of Lemma 3.8 is now complete. Applying Lemma 3.8, we obtain the result of Proposition 3.9.

Proof of Proposition 4.1 We only sketch the proof. As in Genon-Catalot and Jacod (1993) (p. 128-133), we proceed in two steps. First, we prove the result assuming that b, σ are bounded by K. We develop S i /n as follows:

S i n = φ 2 i R i T + T 1 + T 2 + T 3 ,
where, for k = 1, 2, 3, T k = T -1 n j=1 T k,j σ -2 (X i (t j-1 )) and (b(X i (s))b(X i (t j-1 )) + b(X i (t j-1 ))) (M i (s) -M i (t j-1 ))ds.

Moreover, we easily check that, for all θ, E θ ((X i (t) -X i (s)) 2p |φ i ) ≤ C(1 + φ 2p i )|t -s| p .

We have E θ (T 2 |φ) = 0. We proceed as in the proof of Proposition 3.4 to bound the conditional expectations of all terms and get the result of Proposition 4.1 when b, σ are bounded.

  i (u))du b(X i (s))ds, T 2,j = 2φ

Table 1 .

 1 Example 1: σ(.) ≡ 1. Empirical mean and standard deviation (in brackets) of θ (method 1) and θ * (method 2) computed from 100 datasets.

		5.21 (1.17)	5.05 (0.75) 5.04 (0.64)
	λ 3.12 (0.87) 3.14 (0.76)	3.06 (0.47) 3.05 (0.41)
	a * 4.70 (0.99) 5.10 (1.12)	4.60 (0.62) 4.94 (0.62)
	λ * 2.79 (0.64) 3.07 (0.73)	2.73 (0.39) 2.98 (0.39)
		(a 0 = 6, λ 0 = 1)
	ã 6.11 (1.61) 6.09 (1.31)	5.94 (0.86) 5.93 (0.94)
	λ 1.02 (0.29) 1.01 (0.23)	0.99 (0.15) 0.99 (0.14)
	a * 5.43 (1.31) 5.95 (1.25)	5.32 (0.68) 5.80 (0.79)
	λ * 0.89 (0.23) 0.98 (0.22)	0.87 (0.12) 0.96 (0.14)
	N = 50	N = 100
	n = 100	n = 500	n = 100	n = 500
		(a 0 = 5, λ 0 = 3)
	ã 4.80 (1.16) 4.95 (0.99)	4.89 (0.81) 5.08 (0.74)
	λ 2.94 (0.75) 2.99 (0.64)	2.96 (0.49) 3.04 (0.43)
	a * * 4.35 (1.08) 4.85 (0.95)	4.47 (0.68) 4.99 (0.71)
	λ * * 2.61 (0.68) 2.93 (0.61)	2.65 (0.40) 2.97 (0.42)
		(a 0 = 6, λ 0 = 1)
	ã 6.11 (1.51) 6.20 (1.23)	5.80 (0.93) 5.83 (0.77)
	λ 1.01 (0.26) 1.02 (0.22)	0.98 (0.17) 0.98 (0.14)
	a * * 5.43 (1.18) 6.05 (1.17)	5.22 (0.75) 5.71 (0.96)
	λ * * 0.88 (0.21) 1.00 (0.20)	0.86 (0.13) 0.96 (0.13)

  . Statistical Methods for Stochastic Differential Equations. CRC Press. Nie, L. (2006). Strong consistency of the maximum likelihood estimator in generalized linear and nonlinear mixed-effects models. Metrika, 63(2):123-143. Picchini, U., De Gaetano, A., and Ditlevsen, S. (2010). Stochastic differential mixed-effects models. Scandinavian Journal of Statistics, 37:67-90. Picchini, U. and Ditlevsen, S. (2011). Practicle estimation of high dimensional stochastic differential mixedeffects models. Computational Statistics & Data Analysis, 55:1426-1444. Sweeting, T. (1980). Uniform asymptotic normality of the maximum likelihood estimator. Annals of Statistics, 8:1179-1403.

AgroParisTech, UMR518, 75005 Paris, France; maud.delattre@agroparistech.fr

UMR CNRS 8145, Laboratoire MAP5, Université Paris Descartes, Sorbonne Paris Cité, France; valentine.genon-catalot@parisdescartes.fr

Laboratoire Jean Kuntzmann, UMR CNRS 5224, Université Grenoble-Alpes, France; adeline.leclercq-samson@imag.fr c EDP Sciences, SMAI 1999

Then, we use a localization procedure by introducing appropriate stopping times (see Genon-Catalot and Jacod, 1993, p.123). Let

and Proposition 4.1 holds for this process. Afterwards, we conclude using that the sequence of stopping times (R i,ℓ , ℓ > 0) is increasing and tends a.s. to infinity as ℓ tends to infinity.

Auxiliary results

We recall some properties of Gamma and related distributions. The Gamma distribution with parameters (a, λ) (a > 0, λ > 0) , G(a, λ), has density γ a,λ (x) = (λ a /Γ(a))x a-1 e -x 1 (0,+∞) (x), where Γ(a) is the Gamma fonction. The digamma function ψ(a) = Γ ′ (a)/Γ(a) admits the following integral representation:

Consequently, using an integration by part, -aψ ′ (a) = -1 -1 0 t a g(t)dt, where g(t) = (log t/(1t)) ′ . A simple study yields that t a g(t) integrable on (0, 1) and positive except at t = 1. Thus, 1aψ ′ (a) = 0. The following asymptotic expansions as a tends to infinity hold:

), ( 26)

The following results are classical.

Proposition 8.1. If X has distribution G(a, λ), then λX has distribution G(a, 1). For all integer k, E(λX

. If X, Y are independent, X having distribution G(a, 1) and Y having distribution G(b, 1) (a, b > 0), then, X + Y and X/(X + Y ) are independent, X + Y has distribution G(a + b, 1), T = X X+Y has distribution beta of the first kind with parameters (a, b), denoted by β (1) (a, b), and density (1 + z) a+b 1 (0,+∞) (z).

We have E(T ) = a a+b , E(T 2 ) =


Young-Ho Eom

Hang-Hyun Jo

Tail-scope: Using friends to estimate heavy tails of degree distributions in large-scale complex networks

Many complex networks in natural and social phenomena have often been characterized by heavytailed degree distributions. However, due to rapidly growing size of network data and concerns on privacy issues about using these data, it becomes more difficult to analyze complete data sets. Thus, it is crucial to devise effective and efficient estimation methods for heavy tails of degree distributions in large-scale networks only using local information of a small fraction of sampled nodes. Here we propose a tail-scope method based on local observational bias of the friendship paradox. We show that the tail-scope method outperforms the uniform node sampling for estimating heavy tails of degree distributions, while the opposite tendency is observed in the range of small degrees. In order to take advantages of both sampling methods, we devise the hybrid method that successfully recovers the whole range of degree distributions. Our tail-scope method shows how structural heterogeneities of large-scale complex networks can be used to effectively reveal the network structure only with limited local information.

I. INTRODUCTION

Complex networks have served as a powerful mathematical framework to describe complex systems of nature, society, and technology [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Newman | Networks: An Introduction[END_REF][START_REF] Barabási | Network biology: understanding the cells functional organization[END_REF][START_REF] Lazer | Computational social science[END_REF][START_REF] Vespignani | Modelling dynamical processes in complex socio-technical systems[END_REF]. Most complex networks obtained from complex systems are known to be heterogeneous in various aspects [START_REF] Barabási | Emergence of scaling in random networks[END_REF][START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF][START_REF] Newman | Assortative mixing in networks[END_REF][START_REF] Fortunato | Community detection in graphs[END_REF]. One of distinctive heterogeneous features in complex networks is the heavytailed degree distribution: A small number of highly connected nodes coexist with the large number of lowly connected nodes. Highly connected nodes or hubs found in heavy tails have significant roles on the evolution of complex networks and dynamics on such networks. For examples, the existence of hubs leads networks to endemic states in epidemic spreading [START_REF] Pastor-Satorras | Epidemic spreading in scale-free networks[END_REF][START_REF] Castellano | Competing activation mechanisms in epidemics on networks[END_REF], makes networks vulnerable to intended attacks [START_REF] Albert | Error and attack tolerance of complex networks[END_REF], and contributes to the key functions of biological systems [START_REF] Jeong | Lethality and centrality in protein networks[END_REF][START_REF] Han | Evidence for dynamically organized modularity in the yeast protein-protein interaction network[END_REF][START_REF] Zotenko | Why do hubs in the Yeast protein interaction network tend to be essential: Reexamining the connec-tion between the network topology and essentiality[END_REF]. Therefore, identifying the degree distribution and particularly hubs in the heavy tail of degree distribution is the essential step for the network analysis [START_REF] Clauset | Powerlaw distributions in empirical data[END_REF].

Owing to the rapid development of digital technologies, a huge amount of network data is being generated and recorded. In particular, the network data from social media like Twitter and Wikipedia contain tens of millions to billion nodes (users or articles). The role of social media on social dynamics such as public opinion formation, information diffusion, and popularity [START_REF] Centola | The spread of behavior in an online social network experiment[END_REF][START_REF] Bakshy | The role of social networks in information diffusion[END_REF][START_REF] Christakis | The spread of obesity in a large social network over 32 years[END_REF] is getting more crucial, requiring us to timely monitor the large-scale dynamics and to identify the network structure underlying these dynamics [START_REF] Lazer | Computational social science[END_REF][START_REF] Castello | Statistical physics of social dynamics[END_REF]. However, since the social media are constantly growing and changing, the acquisition and analysis of complete network data is an extremely tricky task. Further, increasing public concerns on privacy issues about using these data can inhibit us from analyzing the complete network data [START_REF] Garcia-Herranz | Using friends as sensors to detect global-scale contagious outbreaks[END_REF].

Because of the above difficulties, degree distributions of complex networks need to be estimated based on partial information or by sampling nodes from networks. The simplest method could be to sample nodes randomly, which is called uniform node sampling (UNS). Since the number of nodes corresponding to the tail part of distribution is typically very small, those nodes are rarely sampled, limiting the sampling resolution. Accordingly, much larger statistical fluctuations are expected for the tail part of degree distribution estimated by UNS, when compared to its body part.

The friendship paradox (FP) [START_REF] Feld | Why Your Friends Have More Friends Than Yo Do[END_REF][START_REF] Hodas | Friendship paradox redux: Your friends are more interesting than you[END_REF][START_REF] Eom | Generalized friendship paradox in complex networks: The case of scientific collaboration[END_REF] can shed light on how to effectively estimate the heavy tails of degree distributions. The FP states that the degree of an individual is on average smaller than the average degree of its friends or neighbors. The underlying mechanism behind the FP is the observational bias such that highly connected nodes are more likely to be observed by their neighbors. One can take advantage of this observational bias for the effective sampling of highly connected nodes. Indeed, the group made of friends of randomly chosen nodes turns out to contain highly connected nodes more than the group made of uniformly sampled nodes [START_REF] Eom | Generalized friendship paradox in complex networks: The case of scientific collaboration[END_REF][START_REF] Avrachenkov | Quick Detection of High-degree Entities in Large Directed Networks[END_REF]. Further, the FP has also been used for early detection of contagious outbreaks [START_REF] Garcia-Herranz | Using friends as sensors to detect global-scale contagious outbreaks[END_REF][START_REF] Christakis | Social network sensors for early detection of contagious outbreaks[END_REF] and natural disaster [START_REF] Kryvasheyeu | Performance of Social Network Sensors During Hurricane Sandy[END_REF], and for designing efficient immunization strategy [START_REF] Cohen | Efficient immunization strategies for computer networks and populations[END_REF]. These are mainly based on the observation of activities of highly connected nodes via the FP rather than uniformly sampled nodes.

In this paper, we devise a novel sampling method, called tail-scope, to effectively estimate the heavy tails of degree distributions in large-scale complex networks. We exploit the observational bias of FP as a magnifying glass to observe heavy tails with better resolution and to overcome the resolution limit in the UNS. It is shown that the tail-scope method estimates heavy tails of empirical degree distributions in large-scale networks more accurately than the UNS. Finally, we suggest a hybrid sampling method taking advantages of both UNS and tail-scope methods to recover the whole range of degree distribution.

II. RESULTS

A. Tail-scope: Estimating the tail of degree distribution using the friendship paradox

We consider a directed network G = G(N, L) with N nodes and L directed links. In case of undirected networks, each undirected link is considered as two directed links in both directions. For a node i, the in-degree k i represents the number of incoming links to i from i's inneighbors, and the in-degree distribution is denoted by P (k). Similarly, one can define the out-degree as the number of out-neighbors.

Our goal is to effectively estimate the heavy tail of indegree distribution, i.e., the region of k ≫ 1, by using partial information such as by sampling n nodes with n ≪ N . The observational bias of friendship paradox (FP) indicates that observation via friends can lead to the larger number of high degree nodes than that by the uniform node sampling (UNS), because the chance of a node being observed by its neighbors is proportional to the degree of the node. For this, we randomly choose n directed links and construct a set of nodes reached by following those links. The probability of finding a node of in-degree k in the set is proportional to kP (k) not to P (k), which we denote by Q(k):

Q(k) ∝ kP (k). (1)
Then we obtain the estimated in-degree distribution as

P (k) ∝ Q(k) k . (2)
Thanks to the observational bias of FP, the estimated P (k) has the larger number of highly connected nodes and hence less statistical fluctuation for the tail part than when the UNS is used. Our method can be called tailscope. Precisely, the sampling resolution characterized by the cutoff k c (n) of the distribution is higher for the tail-scope method than for the UNS.

In order to demonstrate the effectiveness of tail-scope method for estimating the heavy tail of the distribution, we consider a network showing the power-law in-degree distribution with power-law exponent α > 2 and minimum in-degree k min :

P (k) = (α -1)k α-1 min k -α , (3)
where we have assumed for convenience that the indegree k is a continuous variable. At first, by randomly choosing n nodes (i.e., by UNS) we obtain the estimated in-degree distribution P UN (k) that is expected to be ∝ k -α . Due to the finiteness of n, we find the natural cutoff to the power-law tail as

P UN (k) ∝ k -α e -k/kc , (4)
where k c can be characterized by the condition

1 n = ∞ kc P (k)dk, (5)
leading to

k c = k min n 1/(α-1) . (6)
Next, for the tail-scope method, we expect from

Q(k) ∝ kP (k) that Q(k) ∝ k -(α-1) e -k/k ′ c , (7)
k ′ c = k min n 1/(α-2) . (8)
Then one gets the estimated in-degree distribution in Eq. ( 2):

P (k) ∝ k -α e -k/k ′ c . (9)
It is evident that the sampling resolution k ′ c for the tailscope case is higher than k c for the UNS, precisely,

k ′ c k c = n 1/[(α-2)(α-1)] > 1. (10)
Therefore, our tail-scope method indeed outperforms the UNS for estimating the tail of the distribution. Since the tail-scope method is based on the uniform link sampling, it can also be called link tail-scope, mainly in order to distinguish from node tail-scope to be discussed in the next Subsection. We numerically test our calculations by constructing the Barabási-Albert (BA) scale-free network [START_REF] Barabási | Emergence of scaling in random networks[END_REF] with N = 10 6 , k min = 2, and α = 3, and then by sampling n = 500 nodes. From the calculations, we expect that k c ≈ 45 and k ′ c ≈ 1000, which are numerically confirmed as shown in Fig 1(A). In the figures, we have used the complementary cumulative distribution function (CCDF), defined as

F (k) = ∞ k ′ =k P (k ′ ), for clearer vi- sualization.

B. Node-based tail-scope method

Our tail-scope method is based on the uniform link sampling. However, in many realistic situations, we can use only the node-based sampling not the link-based sampling. For instance, most application programming interfaces (APIs) of social media like Twitter allow us to retrieve only the user-specific information rather than the relationship-based ones. Thus it is necessary to develop a sampling method using node-based data but aimed to simulate the link tail-scope method.

As social media APIs allow to get only user-specific local information in most cases, we assume that whenever a node is sampled or retrieved, we get the set of in-and out-neighbors for the sampled node. These constraints inevitably introduce correlations between sampled links, implying that any node-based tail-scope methods cannot be exactly mapped to the link tail-scope method. In addition, we assume that the number of retrievals, i.e., sampling size, is strictly limited to n for the fair comparison to other sampling methods, e.g., the UNS. We propose the node tail-scope method as follows.

Node tail-scope method:

• Step 1. Randomly choose n/2 nodes (called primary nodes) from the network and retrieve their out-neighbors to construct a set A of those outneighbors.

• Step 2. Randomly choose n/2 nodes from the set A and retrieve their in-degrees to construct the distribution Q NT (k).

• Step 3. Obtain the estimated in-degree distribution

P NT (k) from Q NT (k)/k.
Here the subscript NT of distributions is the abbreviation of node tail-scope. Note that as the total number of retrievals is limited to n, we use n/2 retrievals for getting out-neighbors, and the rest n/2 retrievals for getting in-degrees. However, there are more high degree nodes sampled than when the UNS is used, leading to the higher resolution for the tail-scope method. For a node sampled several times in Step 2, we consider each sampling as a different case.

By using the same BA network in the previous Subsection, we compare the performance of node tail-scope, shown in Fig. 1(B) to that of link tail-scope in Fig. 1(A). It is observed that there is no significant difference between two results.

C. Performance of the node tail-scope method

In order to empirically compare the performance of node tail-scope method to the UNS, we consider several large-scale complex networks: three undirected networks and four directed networks. For details of these networks, see the Method Section and Table I. From now on, we use the sample size n = 1000 in all cases.

Figure 2 shows estimated in-degree distributions P NT (k) (node tail-scope) and P UN (k) (UNS), in comparison to the original in-degree distribution P (k) obtained from the complete set of nodes in the network. The agreements between original distributions and the distributions by node tail-scope method in the tail parts are remarkable, while some fluctuations are observed in the body parts. On the other hand, the distributions by the UNS show good agreements with the original distributions in the body parts, not in the tail parts. Note that the sample size n = 1000 is much smaller than the network size N ranging from hundreds of thousands to tens of millions nodes (see Table I). We find that the results using n = 2000 and n = 4000 are qualitatively the same as the case of n = 1000.

For the quantitative comparison of performance by different sampling methods, we use Kolmogorov-Smirnov (KS) static D, defined as the maximum difference between two CCDFs. The KS D-static is mainly used as a part of KS test to reject null hypothesis. For example, it has been used to test if a given distribution has a power-law tail [START_REF] Clauset | Powerlaw distributions in empirical data[END_REF]. In this paper, we simply use D-static to measure the agreement between the original in-degree distribution and the estimated in-degree distribution by each sampling method. The D-static for the node tailscope method is obtained as

D NT = max k |F NT (k) -F (k)|, (11)
where F (k) denotes the CCDF of the original in-degree distribution, and F NT (k) denotes the CCDF of P NT (k).

Similarly, D UN is defined for the UNS. The smaller Dstatic implies the better agreement to the original distribution.

Then, we define a p-value to compare the two considered sampling methods. The p-value represents the probability that the distribution by node tail-scope method has the smaller D-static with the original distribution than the distribution by the UNS, i.e., p = Pr(D NT < D UN ).

(

) 12
To focus on the tail part of the distribution, we compare the CCDFs only for the region of k ≥ k 0 , or equivalently for the fraction γ of high degree nodes, where γ = F (k 0 ). The case of γ = 1 corresponds to the comparison for the entire range of in-degree. Figure 3 shows the values of p(γ) for different ranges of in-degree and for each considered network. It is found for all networks that the node tail-scope method clearly outperforms the UNS for the tail parts. The opposite tendency is observed when the entire range of the distribution is compared, because the UNS outperforms the node tail-scope for estimating the body part of the distribution. Since the sample size n is limited, the larger number of high degree nodes for the node tail-scope method results in the smaller number of low degree nodes and hence the larger fluctuations than the case of UNS.

D. Hybrid method for recovering the whole distribution

It is evident that the UNS and the node tail-scope method are good at sampling low and high degree nodes, respectively. In order to take advantages of both methods, we suggest the hybrid method for recovering the whole range of the distribution. It is notable that at Step 1 in our node tail-scope method, n/2 primary nodes are randomly chosen and hence their in-degrees can be utilized for the low degree region. From the primary nodes, we get the in-degree distribution P NT0 (k). Then the hybrid distribution is obtained by

P H (k) = aP NT0 (k) + (1 -a)P NT (k). ( 13
)
The weight parameter a ∈ [0, 1] can be chosen according to which part of the distribution is focused. Here we set as a = 0.5. The hybrid method performs well for the BA network in Fig. 1(B) as well as for empirical networks, two of which are shown in Fig. 4. As expected, the distributions estimated by the hybrid method fit the original distributions better than the UNS for the tail parts, and better than the node tail-scope method for the body parts (see insets in Fig. 4). These findings are also consistent with the values of p(γ) shown in Fig. 5: The larger values of p(γ) for small values of γ in Fig. 5(A) imply the better performance of the hybrid method than the UNS for the tail parts. The larger values of p(γ) for large values of γ in Fig. 5(B) imply the better performance of the hybrid method than the node tail-scope for the body parts. Therefore, we conclude that the hybrid method successfully recovers the whole range of in-degree distributions, by taking advantages of both the UNS and the node tailscope methods. Other values of a = 0.25 and a = 0.75 have been also tested and all results are as expected.

III. DISCUSSION

Modern societies have been shaped by large-scale networked systems like World Wide Web, social media, and transportation systems. Monitoring global activities and identifying the network structure of these systems are of utmost importance in better understanding collective social dynamics. However, increasing size of data from these systems and growing concerns on privacy issues about using these data make the exhausted analysis of complete data sets infeasible. Thus, effective and efficient estimation of large-scale networks based on the small sample size or partial information is necessary. One of the simplest method could be uniform node sampling (UNS). The UNS has drawbacks in particular for estimating the heavy tails of degree distributions, due to the limited sampling resolution and large statistical fluctuations. Since high degree nodes found in the heavy tails are in many cases very important to characterize the structure and dynamics of complex networks, we propose the tail-scope method, which is the effective and efficient sampling method for estimation of heavy tails of degree distributions.

Provided that the sample size is limited, it is inevitable that the larger number of high degree nodes by the tail-scope method leads to the smaller number of low degree nodes than when the UNS is used. In order to take advantages of both the tail-scope and the UNS, we propose the hybrid method to recover the whole range of degree distributions. In this paper, we have considered a very simple form of hybrid method by superposing the estimated degree distributions of the UNS and the tail-scope. It turns out that the hybrid method performs better than the UNS for the tail parts, and better than the tail-scope for the body parts. Devising more general and better hybrid methods will be interesting as a future work, e.g., one can use the degree-dependent weight parameter a in Eq. [START_REF] Jeong | Lethality and centrality in protein networks[END_REF].

Our tail-scope method can be also used for estimating high attribute nodes found in the heavy tail of attribute distribution. The attribute of a node can be its activity, income, happiness, and so on. Recently, the generalized friendship paradox (GFP) has been observed and analyzed in complex networks [START_REF] Eom | Generalized friendship paradox in complex networks: The case of scientific collaboration[END_REF][START_REF] Jo | Generalized friendship paradox in networks with tunable degree-attribute correlation[END_REF]. The GFP states that the attribute of a node is on average lower than the average attribute of its neighbors. In the network showing the positive correlation between degrees and attributes, high degree nodes tend to have higher attributes. It implies that the high attribute nodes are more likely to be observed by their neighbors. Such generalized observational bias can be exploited to effectively estimate high attribute nodes who play important roles, e.g., in early detection of new trends or in designing efficient immunization strategies. Thus, it would be very interesting to generalize our tail-scope method to other attributes of nodes, especially for the large-scale complex networks.

Our tail-scope method shows how structural heterogeneities can help us reveal the network structure only with limited information. By exploiting such heterogeneities of complex networks we can properly evaluate priority and importance of each node in the networks. It is getting more important to better understand the heterogeneities since they are key features characterizing the complexity of large-scale networks.

IV. METHODS

A. Data description

In this paper, we consider seven empirical networks: three of them are undirected and the others are directed. The summary of the networks is presented in Table I. The detailed feature of each network is as following.

Coauthorship. We used a coauthorship network constructed in [START_REF] Eom | Generalized friendship paradox in complex networks: The case of scientific collaboration[END_REF]. The nodes are scientists and the links are formed whenever two scientists coauthored the paper. The network size is N = 242592.

Gowalla. We used a Gowalla friendship network constructed in [START_REF] Cho | Friendship and Mobility: User Movement in Location-Based Social Networks[END_REF]. Gowalla is a location-based social networking service. Each user defines a node. The network size is N = 196562.

LiveJournal. We used a LiveJournal friendship network constructed in [START_REF] Yang | Defining and Evaluating Network Communities based on Ground-truth[END_REF]. Livejournal.com is a social networking service for blog, journal, and diary. The nodes are users of LiveJournal and the users can declare friendship to another user, defining a link. The network size is N = 3997962.

Citation. We used a citation network constructed in [START_REF] Eom | Characterizing and modeling citation dynamics[END_REF]. The network is based on the bibliographic database from 1893 to 2009 provided by American Physical Society (APS). The nodes are articles published in APS journal such as Physical Review Letters or Physical Review E and the directed links represent the citation relation between articles. The network size is N = 463349.

Web graph.

We used a web graph constructed in [START_REF] Leskovec | Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters[END_REF]. The nodes represent webpages in the domains of berkely.edu and stanford.edu domains, and the links are hyperlink between webpages. The network size is N = 685230.

Wikipedia. We used an English Wikipedia network constructed in [START_REF] Eom | Interactions of cultures and top people of Wikipedia from ranking of 24 language editions[END_REF]. The Wikipedia data set was collected in February 2013. The nodes are English Wikipedia articles and the links are hyperlinks between those articles. The network size is N = 4212493.

Twitter. We used a Twitter users network constructed in [START_REF] Kwak | What is Twitter, a social network or a news media[END_REF]. The nodes are Twitter users and the links between users represent the following relations in Twitter. The network size is N = 41652230. FIG. 3: Performance of the node tail-scope method compared to the uniform node sampling for all the considered networks. p(γ) is calculated by Eq. ( 12) but with Kolmogorov-Smirnov D-statics defined only for the range of k ≥ k0, where γ = F (k0).

The smaller γ corresponds to the larger k0. The larger p(γ)-values imply the better performance of the node tail-scope method than the uniform node sampling. To get p-values, we used 1000 realizations of sampling, for each of which the sample size is n = 1000. FIG. 5: Performance of the hybrid method compared to the uniform node sampling (A) and to the node tail-scope method (B) for all the considered networks. p(γ) is calculated by Eq. ( 12) but with Kolmogorov-Smirnov D-statics defined only for the range of k ≥ k0, where γ = F (k0). The smaller γ corresponds to the larger k0. The larger p(γ)-values imply the better performance of the hybrid method than the uniform node sampling (A) or the node tail-scope method (B). To get p-values, we used 1000 realizations of sampling, for each of which the sample size is n = 1000.

FIG. 1 :FIG. 2 :

 12 FIG.1: Comparison of in-degree distributions estimated by uniform node sampling, link tail-scope, node tail-scope, and hybrid methods to the original distribution for the Barabási-Albert scale-free network with N = 10 6 and minimum in-degree kmin = 2. The sample size is n = 500. In all cases, complementary cumulative distribution functions (CCDFs) are presented.

FIG. 4 :

 4 FIG.4: Comparison of in-degree distributions estimated by uniform node sampling, node tail-scope, and hybrid methods to the original distributions for networks of LiveJournal (A) and Twitter (B). The insets show results for the range of k ≤ 20. In all cases, complementary cumulative distribution functions (CCDFs) are presented.

TABLE I :

 I Basic statistics of empirical undirected and directed networks. N denotes the total number of nodes and k denotes the average in-degree. The isolated nodes have been excluded for the analysis.

		1 0.9				BA Coauthor
		0.8				Gowalla LiveJournal
	p(γ)-value	0.4 0.5 0.6 0.7				Citation Web graph Wikipedia Twitter
		0.3
		0.2
		0.1
		0
		0.001	0.01	γ	0.1	1
	Undirected network	N	k Directed network	N	k
	Coauthorship	242592 59.6		Citation	463349 12.2
	Gowalla	196562 9.7		Web graph	685230 12.3
	LiveJournal	3997962 17.3		Wikipedia	4212493 26.4
						Twitter	41652230 36.6

Acknowledgements

The authors thank American Physical Society for providing Physical Review bibliographic data. Y.-H.E. acknowledges support from the EC FET Open project "New tools and algorithms for directed network analysis" (NADINE number 288956). H.-H.J. acknowledges financial support by the Aalto University postdoctoral programme.