Approximations of the allelic frequency spectrum in general supercritical branching populations

TL;DR

This paper develops approximation methods for the allelic frequency spectrum in general supercritical branching populations, analyzing their asymptotic errors and applying the results to genetic data interpretation.

Contribution

It introduces new approximation techniques for the allelic frequency spectrum in complex branching models and evaluates their accuracy through asymptotic analysis.

Findings

Derived central limit theorems for approximation errors

Provided numerical analysis of approximation performance

Applied methods to genetic data analysis, such as EHH

Abstract

We consider a general branching population where the lifetimes of individuals are i.i.d.\ with arbitrary distribution and where each individual gives birth to new individuals at Poisson times independently from each other. In addition, we suppose that individuals experience mutations at Poissonian rate $\theta$ under the infinitely many alleles assumption assuming that types are transmitted from parents to offspring. This mechanism leads to a partition of the population by type, called the allelic partition. The main object of this work is the frequency spectrum $A(k,t)$ which counts the number of families of size $k$ in the population at time $t$. The process $(A(k,t),\ t\in\mathbb{R}_+)$ is an example of non-Markovian branching process belonging to the class of general branching processes counted by random characteristics. In this work, we propose methods of approximation to replace the frequency spectrum by simpler quantities. Our main goal is study the asymptotic error made during these approximations through central limit theorems. In a last section, we perform several numerical analysis using this model, in particular to analyze the behavior of one of these approximations with respect to Sabeti's Extended Haplotype Homozygosity [18].

Equations

$$(A(k,t){)}_{k\ge 1}\approx (c_k{)}_{\ge 1}{e}^{\alpha t}\mathcal{E}.$$

$$A(k,t)\approx (c_k{)}_{\ge 1}{N}_t.$$

$$W(t):=\mathbb{E}\left[N_t\mid N_t>0\right].$$

$${\int }_{[0,\infty )}W(s)e^{-\lambda s}ds=\frac{1}{\psi (\lambda )},\forall \lambda >\alpha ,$$

$$\psi (x)=x-{\int }_{(0,\infty ]}\left(1-e^{-rx}\right)b{\mathbb{P}}_V(dr),x\in {\mathbb{R}}_{+},$$

$${\int }_{{\mathbb{R}}_{+}}{e}^{-\lambda v}{\mathbb{P}}_V(dv)=1+\frac{\psi (\lambda )-\lambda }{b}.$$

$${\int }_{{\mathbb{R}}_{+}}{e}^{-\alpha v}{\mathbb{P}}_V(dv)=1-\frac{\alpha }{b}.$$

$$\mathbb{P}\left(H_1>s\right)=\frac{1}{W(s)},s\in {\mathbb{R}}_{+}.$$

$$\left(X_1,\dots ,X_N\right)\stackrel{d}{=}\left(\widetilde{X_1},\dots ,{\tilde{X}}_N\right),$$

$$N_t=\inf \{i\ge 1\mid H_i>t\},$$

$$\mathbb{P}\left(N_t=k\mid N_t>0\right)=\frac{1}{W(t)}{\left(1-\frac{1}{W(t)}\right)}^{k-1}.$$

$$\mathbb{E}\left[N_t\mid N_t>0\right]=W(t).$$

$$\mathbb{E}{N}_t=W(t)-W\star {\mathbb{P}}_V(t),$$

$$\mathbb{P}\left(N_t>0\right)=1-\frac{W\star {\mathbb{P}}_V(t)}{W(t)},$$

$$W\star {\mathbb{P}}_V(t):={\int }_{[0,t]}W(t-s){\mathbb{P}}_V(ds).$$

$$\mathcal{N}\left(da,di\right)=\text{mathds}{1}_{H_i>t-a}\text{mathds}{1}_{i<{\mathcal{N}}_t}\mathcal{P}\left(di,da\right),$$

$${\psi }_{\theta }(x)=x-{\int }_{(0,\infty ]}\left(1-e^{-rx}\right){\Lambda }_{\theta }(dr)=\frac{x\psi (x+\theta )}{x+\theta }.$$

$$\mathbb{P}\left(Z_0(t)=k\mid Z_0(t)>0\right)=\frac{1}{W_{\theta }(t)}{\left(1-\frac{1}{W_{\theta }(t)}\right)}^{k-1}.$$

$$\mathbb{P}\left(Z_0(t)=k\mid N_t>0\right)=\frac{e^{-\theta t}W(t)}{W_{\theta }(t{)}^2}{\left(1-\frac{1}{W_{\theta }(t)}\right)}^{k-1},\forall k\ge 1.$$

$$e^{-\alpha t}{\psi }^{\prime }(\alpha )W(t)-1=\mathcal{O}\left(e^{-\gamma t}\right).$$

$$W_{\theta }(t)\underset{t\to \infty }{\sim }\frac{e^{\left(\alpha -\theta \right)t}}{{\psi }_{\theta }(\alpha -\theta )}.$$

$$W_{\theta }(t)=\frac{\theta }{\psi (\theta )}+\mathcal{O}\left(e^{-\left(\theta -\alpha \right)t}\right).$$

$$W_{\theta }(t)\underset{t\to \infty }{\sim }\frac{\theta t}{{\psi }^{\prime }(\alpha )}.$$

$$W(t)=\frac{e^{\alpha t}}{{\psi }^{\prime }(\alpha )}-e^{\alpha t}F(t),t\ge 0,$$

$$\underset{t\to \infty }{\lim }{e}^{\alpha t}F(t)=\{\begin{array}{ll}\frac{1}{b\mathbb{E}V-1}& \text{mbox}if\mathbb{E}V<\infty ,\\ 0& \text{mbox}otherwise.\end{array}$$

$$\mathbb{P}\left(\text{NonEx}\right)=\underset{t\to \infty }{\lim }\mathbb{P}\left(N_t>0\right)=\frac{\alpha }{b},$$

$$\prod _{i=1}^{l}A(k_i,t)=\sum _{i=1}^l{\int }_{[0,t]\times \mathbb{N}}\text{mathds}{1}_{Z_0^{(u)}(a)=k_i}\sum _{u_{1:l-1}=1}^{N_{t-a}^{(t)}}\prod _{\underset{i\ne j}{j=1}}^{l-1}{A}^{(u_j)}(k_j,a)\mathcal{N}\left(da,du\right),$$

$$\sum _{u_1=1}^{N_{t-a}^{(t)}}\cdots \sum _{u_{l-1}=1}^{N_{t-a}^{(t)}}.$$

$${\mathbb{E}}_{t}A(k,t)=W(t){\int }_0^t\frac{\theta e^{-\theta s}}{W_{\theta }(s{)}^2}{\left(1-\frac{1}{W_{\theta }(s)}\right)}^{k-1}ds,$$

$$\mathbb{E}A(k,t)A(l,t)$$

Full text

MSC 2000 subject classifications: Primary 60J80; secondary 92D10, 60J85, 60G51, 60K15, 60F05.

Key words and phrases. branching process – splitting tree – Crump–Mode–Jagers process – linear birth–death process – Central Limit Theorem.

1 Introduction

In this paper, we consider a general branching population where the lifetimes of the individuals and their reproductions processes are i.i.d. Moreover, we assume that their lifetimes are distributed according to an arbitrary probability distribution $\mathbb{P}_{V}$ and that the births occur, during their lifetime, according to a Poisson process with rate $b$. The tree underlying this dynamics is called a splitting tree. This class of random trees was introduced in [11] by Geiger and Kersting and has been widely studied in the last decade [14, 15, 16].

We suppose, in addition, that mutations occur on individuals and that each new mutation confers to its holder a brand new type (i.e. never seen in the population): this is the infinitely many alleles assumption. This allows modeling the occurrence of a new type in a population (such as a new species or a new phenotype in a given species). We also suppose that every individual inherits the type of its parent. This model leads to a partition of the population by types. The frequency spectrum of the population alive at time $t$ is defined as the sequence of number $\left(A(k,t)\right)_{k\geq 1}$ where, for each $k$, $A(k,t)$ is the number of families of size $k$ in the population. The famous example of Ewens sampling formula gives explicit expression for the law of the frequency spectrum [9] when the genealogy is given by the Kingman’s coalescent. Other works studied similar quantities in the case of Galton-Waston branching processes (see [4] or [12]). In our model, the frequency spectrum has also been widely studied in the past [6, 7, 8, 5].

Another object of interest is the process $(N_{t},\ t\in\mathbb{R}_{+})$ which counts the number of living individuals in the population at a given time $t$. This process is known as binary homogeneous Crump-Mode-Jagers process. One of the main result of the theory of such process is the law of large number which gives in our particular case that $e^{-\alpha t}N_{t}$ converges almost surely to a random variable $\mathcal{E}$ which is exponential conditionally on non-extinction (for some positive constant $\alpha$).

As for $e^{-\alpha t}N_{t}$, it is also known that the quantities $e^{-\alpha t}A(k,t)$ converge almost surely to $c_{k}\mathcal{E}$, where $c_{k}$ is an explicit constant. This result can be easily obtained by conjunction of the works of [6] and [17] using the theory of general branching processes counted by random characteristics (a complete statement can be found in [8]). An alternative proof avoiding the use of the general branching processes theory can be found in [5].

It appears that the frequency spectrum $(A(k,t))_{k\geq 1}$ is a quantity which is hard to manipulate from the probabilistic point of view (see [6, 7, 5]). This implies that such a model is inconvenient for practical applications. In this work we propose to use the laws of large numbers in order to replace $(A(k,t))_{k\geq 1}$ by more manipulable quantities and propose to investigate the error made during this approximation. The first possible approximation is the following.

Approximation 1:

[TABLE]

However, this is unsatisfactory for practical applications since the random variable $\mathcal{E}$ is not observable at finite times. Another idea is to exploit the fact that the random variable appearing in the law of large numbers for $A(k,t)$ and for $N_{t}$ is the same. This leads to the second approximation.

Approximation 2:

[TABLE]

In order to investigate the errors made during this approximation (at least asymptotically), one would like to have central limit theorems associated to the law of large numbers for the frequency spectrum. In a previous work [13], we have showed that the error in the convergence of $e^{-\alpha t}N_{t}$ is of order $e^{\alpha t/2}$ and obtained a central limit theorem for this error. An important aspect of the method introduced in [13] is that it can be used to derive CLTs for other branching processes counted by random characteristics. In particular, the main goal of this work is to obtain central limit theorems for the convergence of the frequency spectrum. We also study the Markovian cases (when $\mathbb{P}_{V}$ is exponential) where we can obtain more information on the limit distribution.

The original motivation of this study (and of other works on this model [6, 7, 5]) comes from the works of Sabeti and al. [18] where the frequency spectrum is used to detecte positive selection of an allele in an increasing population. More specifically, suppose that you want to detect the positive selection of an allele on a given gene. The main idea is that, under neutral evolution, the allele under consideration needs a long time to reach a high frequency in the population. Hence, if the frequency of the allele w.r.t. its age is significantly higher than the expected frequency (w.r.t. its age and under neutral growth), this anomaly would suggest a positive selection of this allele. The main problem is now to be able to estimate how old the allele is. Sabeti and al. remarked that the allelic partition can be used as a clock to estimate the age of an allele. More precisely, their study begins by selecting a small region of chromosome which characterized the presence of the allele under consideration. Now, the type of an individual, at a distance $x$ (measured in kb) from the core region, is the sequence of $x$ bases following the core region (excluded). As a consequence, the allelic partition of the subpopulation carrying the allele becomes thinner as $x$ increases (because the higher $x$ is, the higher is the probability that a mutation occurred on the sequence of $x$ bases). Finally, the speed of fragmentation of the allelic partition, when $x$ increases, gives clues on the age of the allele. One of the purposes of this model is to understand how the frequency spectrum evolves under neutral evolution. In this work, we discuss some aspects of this method and give some directions in order to construct rigorous tests for the positive selection (see Section 7).

The paper is organized as follows. Section 2 is devoted to the mathematical description of the model and to preliminary results which are used in the sequel. Section 3 gives the mains theoretical results of this work and, in particular, a central limit theorems which allow to study the error in our proposed approximations. Section 4, 5, 6 are devoted to the proofs of Theorem 3.1, 3.3 and 3.5 respectively. Finally, in Section 7 we perform some numerical studies on the model to stress the quality of our approximation. The discussions about the method of Sabeti and al. are given in this last section. An appendix contains some technical proofs and a section which is a reminder of renewal theory.

2 Model and preliminaries

In this work, we consider a branching population with the following dynamic: starting with a single individual (called the ancestor) whose lifetime is distributed according to an arbitrary probability distribution $\mathbb{P}_{V}$, this ancestor gives birth to new individuals at a Poissonian rate $b$. Each birth event giving a single new individual. From this point, each child of the ancestor lives and gives birth according to the same mechanism independently from the other individuals in the population. This formal description can be made rigorous through the definition of a probability distribution on the set of chronological trees. For the details of such construction, we refer the reader to [14]. The first quantity of interest when studying such population is the number $N_{t}$ of alive individuals in the population at a fixed time $t$ (assuming that the time $t=0$ is birth-date of the ancestor). The process $(N_{t},\ t\in\mathbb{R}_{+})$ is known as binary homogeneous Crump-Mode-Jagers process and is a simple example of non-Markovian branching process. In the sequel, we denote by $W(t)$ the expectation of $N_{t}$ conditionally on the non-extinction at time $t$. That is

[TABLE]

In [14], the author shows that the random variable $N_{t}$ is geometrically distributed under $\mathbb{P}_{t}$ with parameter $\frac{1}{W(t)}$. In addition, the author of [14] showed that the Laplace transform of $W$ can be linked to the Laplace transform of $\mathbb{P}_{V}$ through the relation

[TABLE]

where

[TABLE]

and $\alpha$ is the largest root of $\psi$. In particular, the Laplace transform of $\mathbb{P}_{V}$ can be expressed in terms of $\psi$,

[TABLE]

In this work, we assume that $\alpha$ is a strictly positive real number. This case is called the supercritical case and is equivalent to $b\mathbb{E}[V]>1$. In the supercritical case, the real number $\alpha$ is called the Malthusian parameter of the population because it corresponds to the mean exponential growth rate of the population. Before gong further, let us remark that equation (2.2) leads easily to the following identity:

[TABLE]

Many previous works demonstrate [6, 7, 8] that some properties of the splitting tree were easier to study on the tree describing only the genealogical relation between the lineages of the individuals alive at time $t$. For instance, in the model with mutations, the difference between two individuals in term of type lies only on the time past since their lineages has diverged. Hence, this particular genealogical tree, known as coalescent point processes (CPP), contains the essential information to study the allelic partition. In order to derive the law of that genealogical tree, one needs to characterize the joint law of the times of coalescence between pairs of individuals in the population, which are the times since their lineages have split.

In [14], the author defines an order on the set of individuals alive at a fixed time $t$ and consider the sequence of times of coalescences $(H_{i})_{0\leq i\leq N_{t}-1}$ between two consecutive individuals (that is $H_{i}$ is the time passed since the lineage of individuals $i$ and $i+1$ have diverged) with the convention that the older lineage is the first one (i.e. $H_{0}=t$). Moreover, in [14], the author shows that the random vector $(H_{i})_{0\leq i\leq N_{t}-1}$ can be produced from a sequence $(H_{i})_{i\geq 1}$ of i.i.d. random variable stopped at its first value greater than $t$ and such that

[TABLE]

To summarize, given the population is still alive at time $t$, one can forget about the details of the splitting tree and code the genealogy by a new object called the coalescent point process (CPP). Its law is the law of a sequence $\left(H_{i}\right)_{0\leq i\leq N_{t}-1}$, where the family $\left(H_{i}\right)_{i\geq 1}$ is i.i.d. with the same law as $H$, stopped before its first value $H_{N_{t}}$ greater than $t$, and $H_{0}$ is deterministic equal to $t$ (see Figure 1).

Although we do not use directly the CPP in this work, this object allowed us to obtain [5] formulas for the moments of the frequency spectrum which are widely used in the sequel.

Remark 2.1**.**

Let $N$ be a integer valued random variable. In the sequel we said that a random vector with random size $\left(X_{i}\right)_{1\leq i\leq N}$ form an i.i.d. family of random variables independent of $N$, if and only if

[TABLE]

where $\left(\tilde{X}_{i}\right)_{i\geq 1}$ is a sequence of i.i.d. random variables distributed as $X_{1}$ independent of $N$.

Before going further, let us point out that if we define $N_{t}$ as the first value of the sequence $\left(H_{i}\right)_{i\geq 1}$ greater than $t$, i.e.

[TABLE]

then $N_{t}$ is indeed geometric with the expected parameter. More precisely, for a positive integer $k$,

[TABLE]

In particular,

[TABLE]

Moreover, it can be showed (see [17]), that

[TABLE]

and

[TABLE]

where

[TABLE]

Now, let us introduce the mathematical formalism for the mutation process used in this work (this formalism comes from [5]). Since only the mutations occurring on the lineages of living individuals at time $t$ can be observed, it follows from standard properties on Poisson point processes, that the mutation process can be defined directly on the CPP. So, let $\mathcal{P}$ be a Poisson random measure on $[0,t]\times\mathbb{N}$ with intensity measure $\theta\lambda\otimes C$, where $C$ is the counting measure on $\mathbb{N}$, then the mutation random measure $\mathcal{N}$ on the CPP is defined by

[TABLE]

where an atom at $(a,i)$ means that the $i$th branch experiences a mutation at time $t-a$. We suppose that each individual inherits the type of its parent. This rule yields a partition of the population by types. The distribution of the sizes of the families in the population is called the frequency spectrum and is defined as the sequence $\left(A(k,t)\right)_{k\geq 1}$ where $A(k,t)$ is the number of types carried by exactly $k$ individuals in the alive population at time $t$, excluding the family holding the ancestral type of the population (i.e. individuals holding the same type as the root at time [math]). This last family is called clonal, as the ancestral type.

In the study of the frequency spectrum, an important role is played by the law of the clonal family. We denote by $Z_{0}(t)$ the size of this family at time $t$.

To study this family, it is easier to consider the clonal splitting tree constructed from the original splitting tree by cutting every branches beyond mutations. This clonal splitting tree is a standard splitting tree without mutations, where individuals are killed as soon as they die or experience a mutation. The new lifespan law is therefore the minimum between an exponential random variable of parameter $\theta$ and an independent copy of $V$. It is straightforward by simple manipulations of Laplace transforms that the Laplace exponent of the corresponding contour process is

[TABLE]

We denote by $W_{\theta}$ the corresponding scale function. This leads to,

[TABLE]

When $\alpha>\theta$ (resp. $\alpha=\theta$, $\alpha<\theta$), this new tree is supercritical (resp. critical, sub-critical) and we talk about clonal supercritical case (resp. critical, sub-critical case).

Moreover, the law of $Z_{0}$ conditionally on the event $\{N_{t}>0\}$ can be obtained, and is given by

[TABLE]

For the rest of this paper, unless otherwise stated, the notation $\mathbb{P}_{t}$ refers to $\mathbb{P}\left(.\mid N_{t}>0\right)$ whereas $\mathbb{P}_{\infty}$ refers to the probability measure conditioned on the non-extinction event (which has positive probability in the supercritical case).

Finally, we recall the asymptotic behavior of the scale functions $W(t)$ and $W_{\theta}(t)$ which is widely used in the sequel,

Lemma 2.2**.**

*([6, Thm. 3.21])

There exist a positive constant $\gamma$ such that,*

[TABLE]

In the case that $\theta<\alpha$ (clonal supercritical case),

[TABLE]

In the case that $\theta>\alpha$ (clonal sub-critical case),

[TABLE]

In the case where $\theta=\alpha$ (clonal critical case),

[TABLE]

For a purpose, a more precise description of the asymptotic behavior of $W$ is needed. It is given by the following result.

Lemma 2.3**.**

[13*, Prop. 5.1]**

There exists a positive non-increasing càdlàg function $F$ such that*

[TABLE]

and

[TABLE]

From this Lemma and (2.7), one can easily deduce that

[TABLE]

where NonEx refer to the non-extinction event.

In [5], we show that a CPP stopped at time $t$ with scale function $W$ can be constructed by grafting independent CPP stopped at a fixed time $a\leq t$ on a CPP stopped at time $t-a$ with an explicit scale function different of $W$ (see Figure 2).

Moreover, we showed that the frequency spectrum can be expressed as an integral with respect to the random measure $\mathcal{N}$ along the CPP, that is

[TABLE]

where $A^{(u)}(k,a)$ (resp. $Z^{(u)}_{0}$) refers to the frequency spectrum (resp. clonal family) of the $u$th grafted sub-CPP, and $\sum_{u_{1:l-1}=1}^{N^{(t)}_{t-a}}$ denotes for the multi-sum

[TABLE]

Moreover, in [5, Thm, 3.1] we show that the expectation of such integral can be computed easily when the integrand presents local independence properties with the random measure as in formula (2.10). Equation (2.10) is used later to obtain some moments estimates useful to prove our theorems. In particular, this allows to prove that (see [5]) for any positive integer $k$ and $l$,

[TABLE]

and

[TABLE]

These tools also allow, for instance, to prove next two results [14, 8, 5].

Theorem 2.4**.**

There exists a random variable $\mathcal{E}$, such that

[TABLE]

Moreover, under $\mathbb{P}_{\infty}$, $\mathcal{E}$ is exponentially distributed with parameter one.

Theorem 2.5**.**

For any positive integer $k$,

[TABLE]

where $\mathcal{E}$ is the random variable of the Theorem 2.4 and

[TABLE]

3 Main results

The a.s. convergence stated in Section 2 suggests studying the second order properties of the convergence to get central limit theorem. Our main result, Theorem 3.5, allows to study the asymptotic error in the approximation $2$ proposed in the introduction of this work. In addition, we prove more standard central limit theorems which are interesting from the theoretical point of view.

Before going further, we recall that the Laplace distribution with mean $\mu\in\mathbb{R}^{n}$ and covariance matrix $K$ is the probability distribution whose characteristic function is given, for all $\lambda\in\mathbb{R}^{n}$ by

[TABLE]

We denote this law by $\mathcal{L}\left(\mu,K\right)$. We also recall that, if $G$ is a Gaussian random vector with mean $\mu$ and covariance matrix $K$ and $\mathcal{E}$ is an exponential random variable with parameter $1$ independent of $G$, then $\sqrt{\mathcal{E}}G$ is Laplace $\mathcal{L}\left(\mu,K\right)$.

3.1 CLT for the convergence of Theorem 2.5

Theorem 3.1**.**

Suppose that $\theta>\alpha$ and $\int_{[0,\infty)}e^{\left(\theta-\alpha\right)v}\mathbb{P}_{V}(dv)>1$ . Then, we have, under $\mathbb{P}_{\infty}$,

[TABLE]

*where $K$ is some covariance matrix and the constants $c_{k}$ are defined in (2.13). *

The proof of this result can be found in Section 4.

Remark 3.2**.**

We are not able to compute explicitly the covariance matrix $K$ in the general case due to our method of demonstration. However, all our other results give explicit formulas. In particular, the case where $\mathbb{P}_{V}$ is exponential is given by the next theorem. The Yule case is also covered in the following theorem for $d=0$ although it does not satisfy the hypothesis of Theorem 3.1.

Theorem 3.3**.**

Suppose that $V$ is exponentially distributed with parameter $d\in[0,b)$. In this case, $\alpha=b-d$. We still suppose that $\alpha<\theta$, then

[TABLE]

where $K$ is given by

[TABLE]

and

[TABLE]

where $W$, $W_{\theta}$, $\psi^{\prime}(\alpha)$ are defined in the Section 2.

The proof of this result can be found in Section 6. Note that an explicit formula for $\mathbb{E}_{t}A(k,t)$ is given by (2.11). Explicit formulas for $\mathbb{E}_{t}\left[A(k,t)\mathds{1}_{Z_{0}(t)=l}\right]$ can also be found in Proposition 4.5 of [5], and a formula for $\mathbb{E}_{t}\left[N_{a}\mathds{1}_{Z_{0}(t)=k}\right]$ can be found in Proposition 4.1 of [6].

Remark 3.4**.**

The condition on $V$ in Theorem 3.1 is required only to ensure controls of the moments of the considered quantities. However, although the Yule case does not satisfy this condition ($V=\infty$ p.s.) it is included in this last theorem (d=0). This suggests that the condition on $V$ may not be needed.

3.2 CLT for the error between $A(k,t)$ and $c_{k}N_{t}$

The next theorem concerns the error between $A(k,t)$ and $c_{k}N_{t}$. Once again, we have an explicit expression of the covariance matrix of the limit.

Theorem 3.5**.**

Suppose that $\theta>\alpha$, then

[TABLE]

where $M$ is defined in relation (3.3).

The proof of this result can be found in Section 5.

Remark 3.6**.**

We do not known yet if the exponential random variable appearing the Gaussian mixing leading to a Laplace distribution is the same as the exponential limit of $e^{-\alpha t}A(k,t)$. However, the CLT for Markov branching processes in [3] suggest that it is, actually, the case. If, this is true in our case, it would be enough to know the correlations between the limits involved in Theorem 3.1 and 3.5 to obtain an explicit expression for the covariance matrix in Theorem 3.1.

4 Proof of Theorem 3.1

The proof of this theorem is based on the proof of the central limit theorem for the process $(N_{t},\ t\in\mathbb{R}_{+})$ given in [13]. The structure of the proof follows the same lines and is detailed in Section 4 of [13]. In a sake of conciseness, we only highlight the difficulties arising in our new context. The results which are straightforward rewording of the proofs given in [13] are left to the reader. However, we think it is necessary to recall some aspects of [13], in particular from [13, Section 4]. First, we recall that there exists a family $(N^{(i)}_{t},\ t\in\mathbb{R}_{+})_{i\geq 1}$ of i.i.d. population counting processes with the same law as $\left(N_{t},\ t\in\mathbb{R}_{+}\right)$, and a Poisson random measure $\xi$ on $\mathbb{R}_{+}$ with intensity $b\,da$ such that

[TABLE]

where $\xi_{u}=\xi\left([0,u]\right)$. In addition, we have that $t\to\mathbb{E}\left[N_{t}\mathcal{E}\right]$ is the unique solution bounded on finite intervals of the renewal equation,

[TABLE]

and it is given by

[TABLE]

We also recall that equation (4.2) is obtained by taking the product $N_{t}N_{s}$, for some real number $t$ and $s$. Now, equation (4.1) allows to obtain a renewal equation for $\mathbb{E}[N_{t}N_{s}]$ which leads to (4.2) when taking the limit in $s$ of the renormalized equation. We also recall that Lemma 2.3 and equation (2.7) gives

[TABLE]

This also leads, in conjunction with equation (4.3), to

[TABLE]

Finally, let us recall that for any fixed time $u$, there is a natural order (for instance given by the contour process [14]) of the individuals alive at this time. Moreover, we denote, for $1\leq i\leq N_{t}$, $O_{i}^{(u)}$ the residual lifetime of the $i$th individual alive at time $u$. The law of the vector $(O^{(u)}_{2},\dots,O^{(u)}_{N_{u}})$ is given by the following lemma which comes from [13].

Lemma 4.1**.**

Let $u$ in $\mathbb{R}_{+}$, we denote by $O_{i}$ for $i$ an integer between $1$ and $N_{u}$ the residual lifetime of the $i$th individuals alive at time $u$. Then under $\mathbb{P}_{u}$, the family $\left(O_{i},\ i\in\{1,\dots,N_{u}\}\right)$ form a family of independent random variables, independent of $N_{u}$, and, expect $O_{1}$, having the same distribution, given by, for $2\leq i\leq N_{t}$,

[TABLE]

Moreover, it follows that the family $\left(N_{s}(O_{i}),s\in\mathbb{R}_{+}\right)_{1\leq i\leq N_{u}}$ is an independent family of process, i.i.d. for $i\geq 2$, and independent of $N_{u}$.

To end this reminder, let us recall the decomposition of the limiting random variable $\mathcal{E}$ (given for instance in Theorem 2.4) at a fixed time $u$.

Lemma 4.2**.**

[13*, Lemma 6.8]**

We have the following decomposition of $\mathcal{E}$,*

[TABLE]

Moreover, under $\mathbb{P}_{u}$, the random variables $\left(\mathcal{E}_{i}\left(O_{i}\right)\right)_{i\geq 1}$ are independent, independent of $N_{u}$, and identically distributed for $i\geq 2$.

We can now start the proof of theorem 3.1. As in [13], the proof begins by some estimate on moments.

4.1 Preliminary moments estimates

We start by computing the moment in the case of a standard splitting tree. According to [13, Section 4], the next step is to obtain the same kind of estimates in the case of a splitting tree whose ancestor individual has a lifetime distribution which can be different from the rest of the population.

4.1.1 Case $V_{\emptyset}\overset{\mathcal{L}}{=}V$

One of the main difficulties to extend the preceding proof to the frequency spectrum is to get estimates on

[TABLE]

We first study the renewal equation satisfied by $\mathbb{E}A(k,t)\mathcal{E}$ similarly as in [13, Lemma 6.1].

Lemma 4.3** (Joint moment of $\mathcal{E}$ and $A(k,t)$).**

$\mathbb{E}\left[A(k,t)\mathcal{E}\right]$* is the unique solution bounded on finite intervals of the renewal equation,*

[TABLE]

with $X_{t}$ the number of families of size $k$ alive at time $t$ whose original mutation has taken place during the lifetime of the ancestor individual.

Proof.

We recall that $A(k,t)$ is the number of non-ancestral families of size $k$ at time $t$. Similarly, as for $N_{t}$, $A(k,t)$ can be obtained as the sum of the contributions of all the trees grafted on the lifetime of the ancestor individual in addition to the mutations which take place on the ancestral branch, that is,

[TABLE]

where $\left(A(k,t,i),t\in\mathbb{R}_{+}\right)_{i\geq 1}$ is a family of independent processes having the same law as $A(k,t)$. Now, taking the product $A(k,t)N_{s}$ and using the same arguments as in the proof of lemma [13, Lemma 6.1] to take the limit in $s$ leads to the result. In particular, the last term is obtained using that

[TABLE]

∎

The result of Lemma 4.3 is quite disappointing since the presence of the mysterious process $X_{t}$ prevents any explicit resolution of equation (4.3). However, one may note that equation (4.3) is quite similar to equation (4.2) driving $\mathbb{E}N_{t}\mathcal{E}$, so if the contribution of $X_{t}$ in the renewal structure of the process is small enough, one can expect the same asymptotic behavior for $\mathbb{E}A(k,t)\mathcal{E}$ as for $\mathbb{E}N_{t}\mathcal{E}$. Moreover, we clearly have on $X_{t}$ the following a.s. estimate,

[TABLE]

where $Z^{(i)}_{0}$ denote for the ancestral families on the $i$th trees grafted on the ancestral branch. Hence, if we take $\theta>\alpha$ and we suppose $V<\infty$ a.s., one can expect that $X_{t}$ decreases very fast. These are the ideas the following Lemma is based on. Moreover, as it is seen in the proof of the following lemma, the hypothesis $V<\infty$ a.s. can be weakened.

Lemma 4.4**.**

Under the hypothesis of Theorem 3.1, for all $k\geq 1$, there exists a constant $\gamma_{k}\in\mathbb{R}$ such that,

[TABLE]

Proof.

Combining equations $\eqref{eq:renewNt}$ and $\eqref{eq:renewAk}$, we get that,

[TABLE]

which is also a renewal equation. On one hand, using equations (2.5) and (2.11) imply that

[TABLE]

which leads using Lemma 2.2, to

[TABLE]

for some positive real constant $\mathcal{C}$.

The derivative of the r.h.s. of (4.1.1) is given by

[TABLE]

which is equal to

[TABLE]

using Stieljes integration by parts. Now, since,

[TABLE]

this shows that the right hand side of (4.1.1) is decreasing for $t$ large enough. Moreover, it is straightforward to shows that the r.h.s. of (4.1.1) is also integrable. This implies that $\xi_{1}^{(k)}$ is DRI from the same Lemma. On the other hand, it follows from (4.8) that

[TABLE]

Then, we obtain using Cauchy-Schwarz inequality, that

[TABLE]

It follows that we need to investigate the behavior of

[TABLE]

which is equal to

[TABLE]

using [13, Lemma 2.6]. Then, since, from (2.8) and Lemma 2.2,

[TABLE]

it follows, using that the right hand side of (4.1.1) is DRI and Lemma A.1, that $\xi^{(k)}_{2}$ is DRI. Finally, it comes from Theorem A.2, that

[TABLE]

∎

Using the preceding lemma, we can now get the quadratic error in the convergence of the frequency spectrum.

Lemma 4.5** (Quadratic error for the convergence of $A(k,t)$.).**

Let $k$ and $l$ two positive integers. Then under the hypothesis of Theorem 3.1, there exists a family of real numbers $\left(a_{k,l}\right)_{l,k\geq 1}$ such that,

[TABLE]

where the sequence $\left(c_{k}\right)_{k\geq 1}$ is defined by (2.13).

Proof.

Now, noting

[TABLE]

we have, from (2) and Lemma 2.3,

[TABLE]

with

[TABLE]

and $F$, $\mu$ are defined in Lemma 2.3. Now, using (4.4), we have

[TABLE]

which leads to

[TABLE]

Since, by Lemma 2.2

[TABLE]

it follows, combining (4.13), (4.5), and Lemma 4.4, that

[TABLE]

The result follows readily from the fact that $\mathbb{P}\left(N_{t}>0\right)\sim\frac{\alpha}{b}$.

∎

Lemma 4.6** (Boundedness of the third moment).**

Let $k_{1},k_{2},k_{3}$ three positive integers, then

[TABLE]

Proof.

We have,

[TABLE]

Hence, we only have to prove the Lemma for $k_{1}=k_{2}=k_{3}=k$. Hence,

[TABLE]

The last two terms have been treated in the proof of [13, Lemma 6.4], and the boundedness of

[TABLE]

follows from the following Lemma 4.7 and Hölder’s inequality.

∎

Lemma 4.7**.**

For all $k\geq 1$,

[TABLE]

is bounded.

Due to technicality, the proof of this lemma is postponed to the end in appendix.

4.1.2 Arbitrary initial distribution case

The following Lemmas are the counter part of Lemmas 6.5, 6.6, and 6.7 of [13]. They play the same role in the proof of Theorem 3.1 as in the proof of the central limit theorem given in [13]. In the sequel, we denote by $\left(A(k,t,\Xi)\right)_{k\geq 1}$, the frequency spectrum of the splitting tree where the lifetime of the ancestral individual is $\Xi$, in the same manner as for $N_{t}\left(\Xi\right)$ in [13].

[TABLE]

and, let $\mathcal{E}\left(\Xi\right)$ be the random variable defined by

[TABLE]

Lemma 4.8** ($L^{2}$ convergence in the general case).**

Consider the general frequency spectrum $\left(A(k,t,\Xi)\right)_{k\geq 1}$, then, for all $k$, $\psi^{\prime}(\alpha)e^{-\alpha t}A(k,t,\Xi)$ converge to $\mathcal{E}\left(\Xi\right)$ (see 4.17) in $L^{2}$ as $t$ goes to infinity and

[TABLE]

where the convergence is uniform w.r.t. the random variable $\Xi$. In the case where $\Xi$ is distributed as $O^{(\beta t)}_{2}$, for $0<\beta<\frac{1}{2}$, we get

[TABLE]

Lemma 4.9** (First moment).**

The first moments are asymptotically bounded, that is, for all $k\geq 1$,

[TABLE]

uniformly with respect to the random variable $\Xi$.

Lemma 4.10** (Boundedness in the general case.).**

Let $k_{1},k_{2},k_{3}$ three positive integers, then

[TABLE]

uniformly with respect to the random variable $\Xi$.

We do not detail the proofs of these results since they are direct adaptations of the proofs of Lemmas 6.5, 6.6, and 6.7 of [13].

4.2 Proof of the result

The following result is based on the fact that, in the clonal sub-critical case, the lifetime of a family is expected to be small. It follows that one can expect that all the family of size $k$ live in different subtrees as soon as $t>>u$. This is the point of the following lemma.

Lemma 4.11**.**

Suppose that $\alpha<\theta$. If we denote by $\Gamma_{u,t}$ the event,

[TABLE]

then, for all $\beta$ in $(0,1-\frac{\alpha}{\theta})$, we have

[TABLE]

Proof.

The proof of this Lemma, as the calculation of the moments of $A(k,t)$ relies on the representation of the genealogy of the living population at time $t$ as a coalescent point process [5]. Moreover, we denote by $\widetilde{N}^{(t)}_{u}$ the number of living individuals at time $u$ who have alive descent at time $t$. In [5], we showed that, under $\mathbb{P}_{t}$, $\widetilde{N}^{(t)}_{u}$ is geometrically distributed with parameter $\frac{W(t-u)}{W(t)}$.

Now, $\mathds{1}_{\Gamma_{u,t}}$ can be rewritten as

[TABLE]

where $Z_{0}^{i}(t-u)$ denotes the number of individuals alive at time $t$ descending from the $i$th individual alive at time $u$ and carrying its type (the clonal type of the sub-CPP). Moreover, from Proposition 4.3 of [5], we know that that under $\mathbb{P}_{t}$, the family $Z^{(i)}_{0}(t-u)$ is an i.i.d. family of random variables distributed as $Z_{0}(t-u)$ under $\mathbb{P}_{t-u}$, and $\widetilde{N}^{(t)}_{u}$ is independent of $Z^{(i)}_{0}(t-u)$ (still under $\mathbb{P}_{t}$).

Then,

[TABLE]

Using (2.8), some calculus leads to,

[TABLE]

Now, since,

[TABLE]

taking $u=\beta t$, we obtain, using Lemma 2.2 and

[TABLE]

the desired result. ∎

Proof of Theorem 3.1.

Fix $0<u<t$. Note that the event $\Gamma_{u,t}$ of Lemma 4.11 can be rewritten as

[TABLE]

where $Z_{0}^{i}(t-u,O_{i})$ denote the number of individuals alive at time $t$ carrying the same type as the $i$th alive individual at time $u$, that is the ancestral family of the splitting constructed from the residual lifetime of the $i$th individual (see Section 4 in [13]).

Let $K$ be a multi-integer, we denote by $\mathcal{L}^{(K)}$ (resp. $A(K,t)$) the random vector $\left(\mathcal{L}^{k_{1}},\dots,\mathcal{L}^{k_{N}}\right)$ (resp. $\left(A(k_{1},t),\dots,A(k_{N},t)\right)$) with

[TABLE]

On the event $\Gamma_{u,t}$, we have a.s.,

[TABLE]

where the family $\left(A^{(i)}\left(k_{l},t-u,O_{i}\right)\right)_{i\geq 1}$ stand for the frequency spectrum for each subtree, which are independent from Lemma 4.1. Hence, using Lemma 4.2,

[TABLE]

By Lemma 4.1, that the family $\left(A^{i}(k_{l},t-u,O_{i})\right)_{2\leq i\leq N_{u}}$ is i.i.d. under $\mathbb{P}_{u}$.

In the sequel, we denote, for all $l$ and $i\geq 1$,

[TABLE]

Now, let

[TABLE]

From this point, following closely the proof of Theorem 3.2 of [13]. Taking $u=\beta$ in $\left(0,\frac{1}{2}\wedge(1-\frac{\alpha}{\theta})\right)$, the only difficulty is to handle the indicator function $\mathds{1}_{Z_{0}(t-u,O_{i})>0}$ in the Taylor development of $\varphi_{K}$. We show how it can be done for one of the second order terms, and leave the rest of the details to the reader.

It follows from H lder’s inequality that

[TABLE]

from which it follows, using Lemma 4.10, that the r.h.s. of this last inequality is $\mathcal{O}\left(\mathbb{P}\left(Z^{2}_{0}(t-u,O_{2})>0\right)^{\frac{1}{3}}\right)$. Now, using (4.18) and Lemma 4.11, it is easily seen that

[TABLE]

Finally, using Lemma 4.5, we get

[TABLE]

These allow us to conclude that

[TABLE]

where $K_{i,j}$ is given by

[TABLE]

with $K$ is the multi-integer $(k_{1},\dots,k_{N})$, and the $a_{l,k}$s are defined in Lemma 4.5.

To end the proof, note that,

[TABLE]

thanks to Lemma 4.11. ∎

5 Proof of Theorem 3.5

Since all the ideas of the proof of this theorem have been developed the preceding sections, we do not detail all the proof. The only step which needs clarification is the computation of the covariance matrix of the Laplace limit law $\mathcal{M}$. According to the proof of Theorem 3.1, it is given by

[TABLE]

which is equal, thanks to (4.19) and an easy adaptation of Lemma 6.6 in [13], to

[TABLE]

So it remains to get the limit of

[TABLE]

as $t$ goes to infinity. We recall that using the calculus made in the proof of Theorem 6.3 of [5], we have

[TABLE]

Moreover, (4.15) entails

[TABLE]

with

[TABLE]

These identities allow us to obtain

[TABLE]

Taking the limit as $t$ goes to infinity leads to

[TABLE]

Finally, since $\mathbb{P}\left(N_{t}>0\right)\sim\frac{\alpha}{b}$,

[TABLE]

6 Markovian cases

We can get more information on the unknown covariance matrix $K$ in the case where the life duration distribution is exponential. Our study also cover the case $\mathbb{P}_{V}=\delta_{\infty}$ (Yule case), although it does not fit the conditions required by the Theorem 3.1. The reason comes from our method of calculation for $\mathbb{E}\left[A(k,t)\mathcal{E}\right]$. Let us consider the filtration $\left(\mathcal{F}_{t}\right)_{t\in\mathbb{R}_{+}}$, where $\mathcal{F}_{t}$ is the $\sigma$-field generated by the tree truncated above $t$ and the restriction of the mutation measure on $[0,t)$.

Then $N_{t}$ is Markovian with respect to $\mathcal{F}_{t}$ and for all positive real numbers $t\leq s$,

[TABLE]

So that,

[TABLE]

By making a renormalization by $e^{-\alpha s}$ and taking the limit as $s$ goes to infinity, we get,

[TABLE]

since, in the Markovian case, it is known from [8] that

[TABLE]

Suppose first that $d>0$. It follows that,

[TABLE]

By (4.4),

[TABLE]

so

[TABLE]

Finally, since, using Proposition 2.3,

[TABLE]

it follows from (5.2),

[TABLE]

using that $\mu=\frac{1}{b\mathbb{E}V-1}$. In the Yule case, an easy adaptation of the preceding proof leads to

[TABLE]

7 Numerical studies

The purpose of this section is to analyze our approximation method and the estimation of the error by virtue of numerical experiments. There are several practical difficulties appearing when one tries to perform such study.

The first problem, which involves no conceptual difficulties, lies only on the implementation of the formulas appearing in Theorems 3.1, 3.3 and 3.5. In particular, the computation of the moments of type $\mathbb{E}[A(k,t)\mathds{1}_{Z_{0}(t)=l}]$ are particularly complicated (see Proposition 5.4 in [5]).

Another difficulty is to obtain numerical approximations of the scale functions $W$ and $W_{\theta}$. For instance, these functions appear in the computation of the covariance matrix of Theorems 3.3 and 3.5 or when one wants to simulate the coalescent point process. To obtain such approximations, we need to apply numerically the Laplace inversion operator to the functions $\frac{1}{\psi}$ and $\frac{1}{\psi_{\theta}}$.

Unfortunately, the Laplace numerical inversion is a rather difficult problem (see for instance [1] or [2]) which is often computationally expensive. As a consequence, the computational cost of performing multiple numerical integration involving $W$ or $W_{\theta}$ can be important when done with a crude method. Moreover, these methods presents rough numerical instabilities when the original function is exponentially increasing (inverting $\lambda\to\frac{1}{1-\lambda}$, whose inverse is $x\to e^{x}$, is already a tough numerical problem).

For all these reasons, we provide with this work a Matlab toolbox which handle all these difficulties and allows users to perform numerical experiments without having to take care of these issues.

In this whole section, we are interested in the approximation of the frequency spectrum at a fixed time $t$ by the sequence $N_{t}(c_{k})_{k\geq 1}$ (we recall that $c_{k}$ was defined in equation (2.13)). As a consequence, the error in this approximation are computed thanks to Theorem 3.5. The parameters of the model are set as follows:

• •

$\mathbb{P}_{V}$ is a Rice distribution with shape parameter 1 and scale parameter 1.

• •

$b=1$.

• •

$\theta=1$.

For such parameters $\alpha$ approximately equals to $0.5$. Figure 3 shows the evolution of the frequency spectrum (for $k$ between $1$ and $10$) through time. The different quantities seem to growth exponentially with rate $\alpha$ with a time-shift which depend on $k$. An interesting open question would be to understand the behavior of these shifts.

In order to stress our methods of approximation, the first idea is to look to the renormalized frequency spectrum $\left(\frac{A(k,t)}{c_{k}}\right)_{k\geq 1}$ which is expected to look like $(N_{t},\ t\in\mathbb{R}_{+})$. As showed in Figure 4, the approximation seems to be quite accurate for $k=1,2$. However, a more quantitative analysis is required. Figure 5 shows the absolute error in the approximation of $A(1,t)$ by $c_{1}N_{t}$. This error is a little disappointing since it since to diverge when $t$ goes to infinity. However, even if, according to Figure 5, the absolute error at time $20$ if of order $10^{3}$, the relative error shows that this error is quite small with respect to the value of $A(1,20)$.

Another question is about the speed of convergence in the central limit theorem stated in Theorem 3.5. The red curve of Figure 6 shows the density of the Laplace distribution given in Theorem 3.5 in the case of $A(1,t)$ whereas the blue histogram shows the distribution of $\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ for $t=10$ ($\alpha t\sim 5$ and $\mathbb{E}_{t}[N_{t}]\sim 300$) from 10000 simulations. This Figure highlights the fact that even if the taken time $t$ is quite small the distribution $\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ seems already close to the limiting distribution. Figure 7 shows the same kind of behavior in the multidimensional case. To be more quantitative, Figure 8 shows the evolution in time of the distance between the density of limit distribution given in Theorem 3.5 and a kernel estimation of the distribution of $\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ (the estimation is made from $10000$ simulations at each time). This suggest an exponential rate of convergence in Theorem 3.5. In the view of Figure 8, one may think that Berry-Essen type results for Theorem 3.5 would be quite interesting, in particular to understand how the speed of convergence is related to choice of the parameters.

Another interesting question which could be partially probed by simulation is the study of the behavior of the error in the clonal supercritical case. Figure 9 shows a kernel estimation (from $10000$ simulation) of the density of $\psi^{\prime}(\alpha)(e^{\alpha\frac{t}{2}}(A(1,t)-c_{k}N_{t}))$ in the clonal supercritical case ($\theta=0.2<\alpha$). Figure 9 suggest a totally different behavior with a limit distribution which is asymmetric with respect to [math]. In particular, in the view of the shape of the distribution, one could conjecture that the limit is a skew stable distribution.

To end this section, let us goes back to one of the motivation of this work. The following discussion dot not claim to be rigorous and is essentially formal. We recall that the Extended Haplotype Homozygosity (EHH) can be used to detect positive selection in a population [18]. In particular, the behavior of the frequency spectrum in this model gives a standard for the behavior of a subpopulation carrying a common allele under neutral evolution. In order to have a rigorous model to describe this phenomenon, we need to introduce a new mutation measure which is different from the one given in Section 2. We define it directly on the CPP but this could be equivalently defined on the splitting tree. So let $\mathcal{P}$ be a Poisson random measure on $[0,t]\times\mathbb{N}\times\mathbb{R}_{+}$ with intensity measure $\lambda\otimes C\otimes\lambda$, where $C$ is the counting measure on $\mathbb{N}$, then, for any mutation rate $\theta$ in $\mathbb{R}_{+}$, we define the $\theta$-mutation random measure $\mathcal{N}_{\theta}$ by

[TABLE]

where, as before, an atom at $(a,i)$ means that the $i$th branch experiences a mutation at time $t-a$. This construction allows to increase the mutation in consistent manner. This allows to model the type of an individual at a distance $x$ (such that the mutation rate is a function of $x$) from the core haplotype (we refer the reader to [18] for more details). Now, following [5], we can define the frequency spectrum at mutation rate $\theta$ by

[TABLE]

where $Z_{0}(i,a)$ is the number of individual at time $t$ carrying the type of the $i$th individual at time $t-a$ (see [5] for more details). Let us also define $Z_{0}^{\theta}(t)$ the number of individuals carrying the type of the first individual at time [math] when the mutation measure is given by $\mathcal{N}_{\theta}$. As expected, the allelic partition of the population becomes thinner as $\theta$ growth.

Now, the definition of the EHH $G_{\theta}(t)$ is the probability that two uniformly sampled individuals in the population have the same type, that is

[TABLE]

Using that

[TABLE]

this rewrite

[TABLE]

Finally, using the approximation

[TABLE]

proposed in this work, one could expect that

[TABLE]

We stress the fact that the above expression make sens only in the clonal subcritical case (in the other cases the integral in note finite). Now, we can look at the accuracy of this approximation in view of numerical simulation. Figure 10 shows the value of the EHH (when $\theta$ increase) from a simulation of the model (blue curve) and the one obtained using our approximation (red curve). In view of Figure 10, the approximation seems pretty accurate. In order to be more quantitative, Figure 11 shows the relative error between the EHH and its approximation for one simulation. This shows that the error, as least for sufficiently large $\theta$, remains under $8$%. To end, let us highlight that Theorem 3.5 can be used to give confidence intervals for fixed $\theta$ but in order the construct tests of selection from curves like these of Figure 10 one would need to have functional CLT in long time for the process $((A^{\theta}(k,t)-c^{\theta}_{k}N_{t})_{k\geq 1},\ \theta\in\mathbb{R}_{+})$.

Appendix A A bit of renewal theory

The purpose of this part is to recall some facts on renewal equations borrowed from [10]. Let $h:\mathbb{R}\to\mathbb{R}$ be a function bounded on finite intervals with support in $\mathbb{R}_{+}$ and $\Gamma$ a probability measure on $\mathbb{R}_{+}$. The equation

[TABLE]

called a renewal equation, is known to admit a unique solution finite on bounded interval.

Here, our interest is focused on the asymptotic behavior of $F$. We said that the function $h$ is DRI (directly Riemann integrable) if for any $\delta>0$, the quantities

[TABLE]

and

[TABLE]

converge as $n$ goes to infinity respectively to some real number $I_{sup}^{\delta}$ and $I_{inf}^{\delta}$, and

[TABLE]

In the sequel, we use the two following criteria for the DRI property:

Lemma A.1**.**

Let $h$ a function as defined previously. If $h$ satisfies one of the next two conditions, then $h$ is DRI:

$h$* is non-negative decreasing and classically Riemann integrable on $\mathbb{R}_{+}$,*

• 2.

$h$* is càdlàg and bounded by a DRI function.*

We can now state the next result, which is constantly used in the sequel.

Theorem A.2**.**

Suppose that $\Gamma$ is non-lattice, and $h$ is DRI, then

[TABLE]

with

[TABLE]

if the above integral is finite, and zero otherwise.

Remark A.3**.**

In particular, if we suppose that $\Gamma$ is a measure with mass lower than $1$, and that there exists a constant $\alpha\geq 0$ such that

[TABLE]

then, one can perform the change a measure

[TABLE]

in order to apply Theorem A.2 to a new renewal equation to obtain the asymptotic behavior of $F$. (See [10] for details). This method is also used in the sequel.

Appendix B Formula for the fourth moment of the error

Lemma B.1**.**

[TABLE]

Proof.

The proof of this Lemma lies on the calculation of the expectation of each term in the development of

[TABLE]

To make this, we intensively use the relation (2.10) and the method developed in [5]. We begin by computing

[TABLE]

Formula (2.10) gives us,

[TABLE]

The decomposition of the sum in form

[TABLE]

has then been made to distinguish independence properties in our calculation. Actually, as soon as, $i\neq j$, $A^{i}(k,a)$ is independent from $A^{i}(k,a)$ (see [5] for details). It is essential to note that the expectation of these integrals with respect to the random measure $\mathcal{N}$ are all calculated thanks to Theorem 3.1 of [5]. So, taking the expectation now leads to,

[TABLE]

Using the same method for all the other terms and that, for any positive real number $a$ lower than $t$,

[TABLE]

we get Lemma B.1 by reassembling similar terms together. The last term is obtained using the geometric distribution of $N_{t}$ under $\mathbb{P}_{t}$. ∎

Appendix C Boundedness of the fourth moment

Lemma C.1**.**

We begin the proof of the boundedness of the fourth moment by some estimates.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proof.

Relation (i) is easily obtained using the expectation of $N_{t}$ and $A(k,t)$ using (2.11), (2.13) and the behaviour of $W$ provided by Proposition 2.3. The relation (iii) has been obtained in the proof of Theorem 6.1 in [5]. The two last relations are easily obtained from (2.4), (2.8) and Lemma 2.2. The relation (ii) is obtained using the following estimation,

[TABLE]

We begin the proof by computing the r.h.s. of the previous inequality using the same techniques as in Appendix A.

[TABLE]

[TABLE]

Finally,

[TABLE]

Now, an analysis similar to the one of Lemma 4.7 leads to the result. ∎

Proof of Lemma 4.7.

The ideas of the proof, is to analyses one to one every terms of the expression of

[TABLE]

given by Lemma B.1 using Lemma C.1 to show that they behave as $\mathcal{O}\left(W(t)^{2}\right)$. Since the ideas are the same for every terms, we just give a few examples.

First of all, we consider

[TABLE]

Using Lemma C.1 (ii), we have

[TABLE]

Now take the term

[TABLE]

we have from Lemma C.1 (i) and (iv),

[TABLE]

Every term in $W(t)$ or $W(t)^{2}$ are treated this way. Now, we consider the term in $W(t)^{4}$ which is

[TABLE]

since $N_{t}$ is geometrically distributed under $\mathbb{P}_{t}$, and that

[TABLE]

On the other hand, using the law of $Z_{0}(t)$ given by (2.8) and the expectation of $A(k,t)$ given by (2.11) (under $\mathbb{P}_{t}$), we have,

[TABLE]

Finally,

[TABLE]

The last example is the most technical and relies with the term in $W(t)^{3}$, which is, using (C.1) and Lemma B.1,

[TABLE]

On the other hand, using the calculus made in the proof of Theorem 6.3 of [5], we have

[TABLE]

Substituting this last expression in $J$ leads to

[TABLE]

Using many times that,

[TABLE]

thanks to (2.8), (2.11), and (2.6), we finally get

[TABLE]

This shows that $J$ is $\mathcal{O}\left(W(t)^{2}\right)$. ∎