# Stochastic Methodology for the Study of an Epidemic Decay Phase, Based on a Branching Model

Authors: Sophie Pénisson; Christine Jacob
Journal: International Journal of Stochastic Analysis
 (2012)

---

## Abstract

We present a stochastic methodology to study the decay phase of an epidemic. It is based on a general stochastic epidemic process with memory, suitable to model the spread in a large open population with births of any rare transmissible disease with a random incubation period and a Reed-Frost type infection. This model, which belongs to the class of multitype branching processes in discrete time, enables us to predict the incidences of cases and to derive the probability distributions of the extinction time and of the future epidemic size. We also study the epidemic evolution in the worst-case scenario of a very late extinction time, making use of theQ-process. We provide in addition an estimator of the key parameter of the epidemic model quantifying the infection and finally illustrate this methodology with the study of the Bovine Spongiform Encephalopathy epidemic in Great Britain after the 1988 feed ban law.


---

## Body

## 1. Introduction

Outbreaks of infectious diseases of animals or humans are subject, when possible, to control measures aiming at curbing their spread. Effective measures should force the epidemic to enter its decay phase and to reach extinction. The decay phase can then be simply detected by a decrease of the number of cases, when this decrease is obvious. However this is not always the case, and this rough qualitative information might not be sufficient to evaluate accurately the effectiveness of the proposed measures to reduce the final size and duration of the outbreak. The goal of this paper is to present a stochastic methodology in discrete time to study more accurately the decay phase of an epidemic. Our framework is the spread, in a large open population, of a rare transmissible disease such that the infection process may be assumed to follow a Reed-Frost type model, with a probability for a susceptible to become infected by a given dose of pathogens inversely proportional to the total population size. Moreover the latent period (during which an individual is infected but not yet infectious) may be random and long compared to the generation time. Questions about the decay phase include the following: which quantitative criteria can ensure that the disease has entered an extinction phase? What is the probability distribution of the epidemic extinction time, of the epidemic final size, and of the incidence of infected individuals? Finally, what would be the evolution of the epidemic in the event of a very late extinction of the disease?From a practical point of view, it is generally impossible to observe all infections. Susceptible and infected but not yet infectious individuals are most often not distinguishable, being both apparently healthy. This leads to the fact that the only available observations correspond to the incidence of individuals with clinical symptoms. One way to deal with this lack of information was proposed in [1] by Panaretos, who used a model taking into account two types of infected individuals, the observed and the unobserved. In order to answer the previous questions, we choose here a different approach, considering a stochastic model depending on the sole incidences {Xn}n of infectives at each time. We assume that an infective can transmit the disease during one given time unit at most. Therefore the incidence of infectives corresponds to the incidence of cases. The process then describes in a recursive way how one single infective can indirectly generate new infectives (so-called “secondary cases”) k time units later, where 1⩽k⩽d. We assume that this number of secondary cases follows a Poisson probability distribution with parameter Ψk>0. The recursive formula defining {Xn}n is then the following:
(1.1)Xn=∑k=1d‍∑i=1Xn-k‍Yn-k,n,i,
where the variable Yn-k,n,i is the incidence of secondary cases produced at time n with a delay k (latent period) by individual i infectious at time n-k. The {Yn-k,n,i}i,k are assumed independent given ℱn-1:=σ({Xn-k}k⩾1), and the {Yn-k,n,i}i are assumed i.i.d. (identically and independently distributed) given ℱn-1, with a common Poisson distribution with parameter Ψk. This model is therefore time homogeneous and is in this sense less general than the one introduced in [2], which describes the spread of infectious animal diseases in a varying environment. However since we focus on the extinction phase only, the assumption of a constant environment with no new control measure is well founded and enables us to describe more accurately the decay phase. This process is the generalization of the well-known single-type BGW (Bienaymé-Galton-Watson) branching process, which is the limit, as the total population size tends to infinity, of the process describing the spread in a closed population of an infectious disease with a negligible latent period and a probability to become infected following a Reed-Frost model (see, e.g., [3, 4] and citations therein).The core of the paper lies in Section2, where the whole methodology is presented. We first formulate the epidemic model {Xn}n as a multitype branching process with Poissonian transitions, the types representing the memory of the process. This formulation provides useful analytical results such as the extinction criteria, and the distributions of the extinction time and of the epidemic size (Sections 2.1–2.3). Then, in order to investigate the worst-case scenario of an extreme late extinction of the epidemic, we introduce in Section 2.4 the Q-process {Xn*}n, obtained by conditioning {Xn}n on a very late extinction. Using this process, we focus the study on the early behavior of the decay phase in the worst-case scenario, rather than on its long range behavior, which would have little meaning in our setting. Motivated by practical applications to real epidemics, for which we want to predict the processes {Xn}n and {Xn*}n, as well as the derived distributions above, we need to know the values of the parameters {Ψk}k. We may write Ψk=∑a=1am-k‍θaPinc,a(k)Page(a+k), where θa is the mean number of individuals infected at age a by an infective by direct or indirect transmission, am is the largest survival age, Pinc,a(k) is the probability for the individual aged a at infection to have a latent period equal to k (given his survival), and Page(a+k) is the probability to be aged a+k at the end of the latent period. Parameters {θa}a are the key quantities for the spread of the disease and can be subject to changes due to control measures during the epidemic. We assume here that θa=θ0+pa, where pa is constant over time, while θ0 may change with control measures. A typical example is when θ0 is the mean number of individuals infected by an infective by horizontal route at age a, assumed independent of a⩾2, and p1 represents the maternal transmission probability pmat. In this case
(1.2)Ψk=θ0∑a=k+1am‍Pinc,a-k(k)Page(a)+pmatPinc,1(k)Page(k+1).
So we assume here that, except for θ0, the other parameters of the {Ψk}k are constant over time and are known (generally estimated) from previous experiences or from the study of the whole epidemic evolution, in particular its growth phase. We moreover assume that each Ψk depends affinely on θ0 (see Section 3.1). In Sections 2.5 and 2.6, we provide optimal WCLSE (Weighted Conditional Least Squares Estimators) of θ0 in the decay phase, in the frame of {Xn}n as well as in the frame of the associated conditioned process {Xn*}n, and prove the strong consistency and the asymptotic normality of these estimators.The final Section3 is devoted to the application of this method to real epidemics. We first present in Section 3.1 some general conditions under which the spread of a SEIR disease (susceptible, exposed (latent), infectious, removed) can be approximated by our epidemic process defined by (1.1) and give an explicit derivation of the parameters {Ψk}k. We then illustrate the methodology in Section 3.2 with the decay phase of the BSE (Bovine Spongiform Encephalopathy) epidemic in Great Britain. According to the available data [5], the epidemic is obviously fading out. We assume that the {Ψk}k satisfy (1.2). Then thanks to the stochastic tools developed here, we provide in addition to this rough information, short- and long-term predictions about the future spread of the disease as well as an estimation of a potentially remaining horizontal infection route after the 1988 feed ban law.


## 2. Methodology for the Study of an Epidemic Decay Phase

In this section we present a general methodology to study the decay phase of a SEIR disease in a large population, modeled by the process (1.1) defined in Section 1. Our main goal is to provide analytical tools to evaluate the efficiency of the last control measures taken prior to the considered time period. Most of our results are derived from the fact that this epidemic model can be seen as a multitype branching process. Indeed, {Xn}n defined by (1.1) is a Markovian process of order d. Consequently, the d-dimensional process {Xn}n defined by
(2.1)Xn:=(Xn,Xn-1,…,Xn-d+1)
is Markovian of order 1, and it stems directly from (1.1) that {Xn}n is a multitype Bienaymé-Galton-Watson (BGW) process with d types (see, e.g., [6]). Note that the d types in this branching process do not correspond to any attribute of the individuals in the population, which is usually the case in mathematical biology (see, e.g., [7]), but simply correspond to the memory of the process {Xn}n. The information provided by the d-dimensional Markovian process {Xn}n is therefore the same as the one given by the 1-dimensional d-Markovian process {Xn}n, but the multitype branching process setting gives us powerful mathematical tools and results stemming from the branching processes theory [6]. The first basic tool is the generating function of the offspring distribution of {Xn}n, f:=(f1,…,fd), defined on [0,1]d by fi(r):=𝔼[rX1∣X0=ei], where ei:=(0,…,1,…,0) denotes the ith basis vector of ℕd and uv:=∏i=1d‍uivi for u,v∈ℕd. For all r∈[0,1]d, we have here
(2.2)fi(r):=e-(1-r1)Ψiri+1,i=1⋯d-1,fd(r):=e-(1-r1)Ψd.
The second basic tool is the mean matrix M defined by 𝔼(Xn∣ℱn-1)=Xn-1M, which is here
(2.3)M=(Ψ110⋯0Ψ201⋯0⋮⋮⋱⋮Ψd-10⋯⋯1Ψd0⋯⋯0).
Let us notice that, since Ψk>0 for each k=1⋯d, then {Xn}n is nonsingular, positive regular (see [6]) and satisfies the XlogX condition,
(2.4)𝔼[∥X1∥ln∥X1∥∣X0=ei]<∞,i=1,…,d,
where ∥·∥ denotes the sup norm in ℝd.


### 2.1. Extinction of the Epidemic

#### 2.1.1. Almost Sure Extinction

Since the single-type process{Xn}n has a memory of size d, it becomes extinct when it is null at d successive times, or equivalently as soon as the d-dimensional process {Xn}n reaches the d-dimensional null vector 0. According to the theory of multitype positive regular and nonsingular BGW processes [6], the extinction of the process {Xn}n occurs almost surely (a.s.), if and only if ρ⩽1, where ρ is the dominant eigenvalue (also called the Perron's root) of the mean matrix M. Thus ρ is solution of ∑k=1d‍Ψkρ-k=1. In general for d>1, ρ has no explicit expression. However, ∑k=1d‍Ψkρ-k=1 leads directly to the following explicit threshold criteria.Proposition 2.1.
The epidemic becomes extinct almost surely if and only ifR0⩽1, where R0:=∑k=1d‍Ψk is the total mean number of secondary cases generated by one infective in a SEIR disease. We call R0 the basic reproduction number.
Moreover, whenρ⩽1, then R0⩽ρ with equality if and only if either ρ=1 or d=1, and when ρ>1, then R0⩾ρ with equality if and only if d=1.Note that whend>1, R0 only provides information about the threshold level, whereas ρ provides an additional information about the speed of extinction of the process, as shown in the next two paragraphs.


#### 2.1.2. Speed of Extinction

Thanks to well-known results in the literature about multitype branching processes and more particularly to the Perron-Frobenius theorem (see, e.g., [6]), we can deduce the expected incidence of infectives in the population at time n, for n large. Denoting by u and v the right and left eigenvectors of M associated to the Perron's root ρ; that is, MuT=ρuT and vM=ρv, with the normalization convention u·1=u·v=1, where u·v stands for the usual scalar product in ℝd and where the superscript T denotes the transposition, then 𝔼(Xn∣X0)=X0Mn~n→∞ρnX0uTv. The first coordinate in the latter formula becomes for the epidemic process ρn∑i=1d‍X-i+1uiv1. Computing explicitly u and v, we obtain that for all i=1⋯d,
(2.5)ui=∑k=idΨkρ-k+i∑j=1d∑k=jdΨkρ-k+j,vi=ρ-i∑j=1d∑k=jdΨkρ-k+j∑j=1d∑k=jdΨkρ-k,
which leads to the following asymptotic result:
(2.6)𝔼(Xn∣X0)~n→∞ρn∑i=1d‍X-i+1∑k=idΨkρ-k+i-1∑j=1d∑k=jdΨkρ-k.
Hence if ρ<1, the mean number of infectives decreases exponentially at the rate ρ. In the following section, we provide a much finer result on the estimation of the disease extinction time in the population.


#### 2.1.3. Extinction Time of the Epidemic

The extinction time distribution can be derived as a function of the offspring generating function. As usual in stochastic processes, this quantity is calculated conditionally on the initial valueX0=(X0,X-1,…,X-d+1), but for the sake of simplicity we do not let it appear in the notations. Note that since we are building tools for the prediction of the spread of the disease, the time origin 0 corresponds here to the time of the last available data (generally the current date). Let T:=inf{n⩾1,Xn=0} denote the extinction time of the process {Xn}n, and let fn:=f∘fn-1 be the nth iterate of the generating function f given by (2.2). We denote fn:=(fn,1,…,fn,d). Then, by the branching property of the process [6], the probability of extinction of the epidemic before time n is immediately given by ℙ(T⩽n)=ℙ(Xn=0)=fn(0)X0, that is to say,
(2.7)ℙ(T⩽n)=(fn,1(0))X0(fn,2(0))X-1⋯(fn,d(0))X-d+1.
It can be immediately deduced from convergence results for fn(0) as n→∞ [8], that if ρ=1, ℙ(T⩽n)~1-(nη)-1X0·u, while if ρ<1, ℙ(T⩽n)~1-ρnγX0·u, for some constants η,γ>0. As a consequence, the closer ρ is to unity, the longer the time to extinction will be in most realizations. More specifically, (2.7) enables the exact computation (resp., estimation) of ℙ(T⩽n) for any n by the iterative computation of fn, X0 being given, when the parameters Ψk of (1.1) are known (resp., estimated). Moreover, since for ρ⩽1 the epidemic becomes extinct in an a.s. finite time and ℙ(T⩽n)↗n→∞1, then for any given probability p∈[0,1[ there exists n∈ℕ such that ℙ(T⩽n)⩾p. So in practice, for any p∈[0,1[, (2.7) enables us to compute the p-quantile npT of the extinction time,
(2.8)npT:=min{n⩾1:ℙ(T⩽n)⩾p}.


### 2.2. Total Size of the Epidemic

Under the assumptionρ⩽1 and the independence of the {Yn-l,n,i}i,l,n (we previously assumed the independence of the {Yn-l,n,i}i,l, for each n), we derive the distribution of the future total size N:=∑n=1T‍Xn of the epidemic until its extinction, that is, the future total number of infectives until the extinction of the disease. It can be shown [9] that, given the initial value X0, the time origin being the same as in Section 2.1, the probability distribution of N is
(2.9)N=𝒟∑k=1d‍∑i=1X-k+1‍(∑j=1Yk,i‍Nk,i,j),Yk,i:=∑l=kd‍Y-k+1,-k+1+l,i,
where =𝒟 denotes the equality in distribution, an empty sum is by convention 0, the {Yk,i}k,i are independent, and the {Yk,i}i and the {Nk,i,j}k,i,j are i.i.d. with
(2.10)Yk,i=𝒟Poiss(∑l=kd‍Ψl),Nk,i,j=𝒟Borel-Tanner(∑l=1d‍Ψl,1),
that is, for each n⩾1, ℙ(Nk,i,j=n)=e-n∑l=1d‍Ψl(n∑l=1d‍Ψl)n-1(n!)-1. Consequently, under the convention that an empty product is 1, the probability distribution of N is, for any n∈ℕ,
(2.11)ℙ(N=n)=∑∑k=1d∑i=1X-k+1∑j=1yk,ink,i,j=n{0⩽yk,i⩽n,{1⩽nk,i,j⩽n}j}i,k:‍∏k=1d‍∏i=1X-k+1‍e-∑l=kdΨl(∑l=kdΨl)yk,iyk,i!×∏j=1yk,i‍e-nk,i,j∑l=1dΨl(nk,i,j∑l=1dΨl)nk,i,j-1nk,i,j!,
which may be calculated (resp., estimated), replacing the Ψk by their values (resp., estimations). In practice, for any p∈[0,1[, (2.11) enables to compute the p-quantile npN of the total epidemic size,
(2.12)npN:=min{n⩾1:ℙ(N⩽n)⩾p}.
We obtain moreover an explicit formula for the mean value and variance of the size of the epidemic,
(2.13)𝔼(N)=∑k=1dX-k+1∑l=kdΨl1-∑l=1dΨl,𝕍ar(N)=∑k=1dX-k+1∑l=kdΨl(1-∑l=1dΨl)3.


### 2.3. Exposed Population

Depending on the disease, it might also be crucial to study and predict the evolution of the incidence of exposed individuals in the population, which is generally unobservable. We assume that this information is given by the process{Zn}n defined by the conditional distribution Zn∣Xn=𝒟Poiss(Ψ0Xn), where Ψ0 is the mean number of individuals infected at time n by an infective of this time (see Section 3.1). This property enables on the one hand to reconstruct the whole past epidemic (i.e., the incidence of infectives as well as of exposed individuals) thanks to the observable data. On the other hand, it allows to simulate the evolution of the incidence of exposed individuals in the future, based on predictions of the evolution of the epidemic process {Xn}n.


### 2.4. Worst-Case Scenario: Very Late Extinction of the Epidemic

Even in the case when the epidemic dies out almost surely (ρ⩽1), and although one can provide the p-quantile npT of the extinction time with the probability p as large as wanted (see (2.8)), the epidemic might become extinct after this time with a small but nonnull probability of order 1-p. This raises the following question: how would the incidences of infectious and exposed individuals evolve in the (unlikely) case of a very late extinction? In terms of risk analysis, this issue appears to be crucial to evaluate the risks associated with this worst-case scenario. The tools developed in the previous subsections allow to evaluate the probability of all possible outcomes. But since the worst ones, typically a very late extinction, have a negligible probability, these tools do not bring any information in these worst cases and in particular do not inform on the evolution at each time step of the spread of the disease (would it decrease extremely slowly, stay at a constant rate for a very long time, present several peaks in its evolution, etc.?). In order to study the propagation of the epidemic in the decay phase, assuming that extinction occurs very late, we are interested in the distribution of the process {Xn}n conditionally on the event that the epidemic has not become extinct at time k, where k is very large. We therefore consider for any n1,n2,…∈ℕ and any i0,i1,i2,…∈ℕd the conditioned probability ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0), where the subscript i0 denotes the initial value. If k is finite this distribution cannot be easily handled due to its time inhomogeneity. However, when ρ⩽1, it is known ([10]) that this conditioned distribution converges, as k→∞, to the distribution of a d-dimensional Markov process {Xn*}n:
(2.14)limk→∞ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)=ℙi0(Xn1*=i1,…,Xnr*=ir).
We will further discuss in Proposition 2.5 the relevancy of approximating the conditioned probability for k fixed by the limiting object (2.14). The conditioned process {Xn*}n defined by (2.14) is known in the literature as the Q-process associated with {Xn}n, also described as the process conditioned on “not being extinct in the distant future.” It has the following transition probability [10]: for every n⩾1, i,j∈ℕd, i≠0,
(2.15)ℙ(Xn*=j∣Xn-1*=i)=1ρj·ui·uℙ(Xn=j∣Xn-1=i),
where u is the normalized right eigenvector of M associated to the Perron's root ρ as introduced in Section 2.1, and computed explicitly in (2.5). In the same way as for the process {Xn}n, we define the 1-dimensional process Xn*:=Xn,1*. By construction we then have Xn,i*=Xn-i+1*, for each n and each i=1⋯d.Proposition 2.2.
The stochastic process{Xn*}n is, conditionally on its past, distributed as the sum of two independent Poisson and Bernoulli random variables:
(2.16)Xn*∣Xn-1*=𝒟Poiss(Xn-1*·Ψ)*ℬ(p(Xn-1*)),
where Ψ:=(Ψ1,…,Ψd), * is the convolution product symbol, and
(2.17)p(Xn-1*):=u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk.Proof.
Applying (1.1) and (2.15), we obtain that for all j∈ℕ,
(2.18)ℙ(Xn*=j∣Xn-1*)=ℙ(Xn*=(j,Xn-1*,…,Xn-(d-1)*)∣Xn-1*)=ju1+∑k=2dXn-k+1*ukρXn-1*·uℙ(Xn=(j,Xn-1*,…,Xn-(d-1)*)∣Xn-1=Xn-1*)=ju1+∑k=2dXn-k+1*ukρXn-1*·u(Xn-1*·Ψ)jj!e-Xn-1*·Ψ=u1Xn-1*·ΨρXn-1*·u(Xn-1*·Ψ)j-1(j-1)!e-Xn-1*·Ψ+∑k=2dXn-k+1*ukρXn-1*·u(Xn-1*·Ψ)jj!e-Xn-1*·Ψ.
The equality MuT=ρuT implies that for all k=1⋯d-1, ρuk=Ψku1+uk+1, and that ρud=Ψdu1. Consequently, ρXn-1*·u=u1Xn-1*·Ψ+∑k=2d‍Xn-k+1*uk, and thus
(2.19)ℙ(Xn*=j∣Xn-1*)=u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk(Xn-1*·Ψ)j-1(j-1)!e-Xn-1*·Ψ+(1-u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk)(Xn-1*·Ψ)jj!e-Xn-1*·Ψ=[Poiss(Xn-1*·Ψ)*ℬ(u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk)](j).Remark 2.3.
Note that if one compares (2.16) with the transition probability of the unconditioned process Xn∣Xn-1=𝒟Poiss(Xn-1·Ψ), it appears that {Xn*}n behaves at each time step like {Xn}n, according to a Poisson distribution, except that it has the possibility at each time step to add one unit or not, according to a Bernoulli random variable. Moreover, if Xn-1*=⋯=Xn-(d-1)*=0, then according to (2.17), p(Xn-1*)=1, which implies that at time n, the probability to add one unit is equal to one, thus preventing the extinction of the process.Proposition 2.4.
The process{Xn*}n admits a stationary probability measure π with finite first- and second-order moments.Proof.
Since the multitype branching process{Xn}n satisfies property (2.4), it is known [10] that the Q-process {Xn*}n is positive recurrent with a stationary probability measure π given by
(2.20)π(i)≔i·uν(i)∑k∈ℕdk·uν(k),i∈ℕd,
where ν is the Yaglom distribution of the process {Xn}n, uniquely defined by the following property: for all i,j∈ℕd∖{0}, limn→∞ℙ(Xn=i∣X0=j,Xn≠0)=ν(i). In the literature, this stationary measure for the conditioned process {Xn*}n is also referred to as the doubly limiting conditional probability. Moreover, by Proposition 2.2,
(2.21)𝔼(Xn*)=𝔼[𝔼(Xn*∣Xn-1*)]=𝔼[Xn-1*·Ψ+p(Xn-1*)]⩽∑k=1d‍𝔼(Xn-k*)Ψk+1,
which implies that lim¯n→∞𝔼(Xn*)⩽(1-∑k=1d‍Ψk)-1<∞. We consequently obtain by means of Fatou's lemma that, for every i=1⋯d,
(2.22)∑j∈ℕd‍jiπ(j)=𝔼(limn→∞Xn,i*)=𝔼(limn→∞Xn-i+1*)=𝔼(limn→∞Xn*)⩽lim_n→∞𝔼(Xn*)<∞.
We similarly prove that π has finite second-order moments by writing
(2.23)𝕍ar(Xn*∣Xn-1*)=Xn-1*·Ψ+u1Xn-1*·Ψ∑k=2dXn-k+1*uk(u1Xn-1*·Ψ+∑k=2dXn-k+1*uk)2⩽Xn-1*·Ψ+14.Let us discuss the relevancy of approximating the epidemic process{Xn}n conditioned on nonextinction at some finite time k, for k large, by the Q-process {Xn*}n obtained by letting k→∞. When considering the case of late extinction, one works under an hypothetical assumption based on the unknown future, hence in practice one does not focus on a specific value k for the survival of the disease in the population. We therefore might consider that k is chosen large enough such that the approximation of the process {Xn}n conditioned on the event {Xk≠0} by the process {Xn*}n is valid. Of course, the order of magnitude of such k will depend on the rate of convergence of the conditioned process to {Xn*}n.Proposition 2.5.
Letn1⩽⋯⩽nr⩽k and i0,…,ir∈ℕd∖{0}. Then the difference |ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)-ℙi0(Xn1*=i1,…,Xnr*=ir)| decreases, as k→∞, with max{ks-1(|λ|ρ-1)k/2,ρk/2}, where λ is an eigenvalue of M such that ρ>|λ|⩾|λ3|⩾|λ4|⩾⋯, with the λi being the other eigenvalues of M. In case |λ|=|λ3| we stipulate that the multiplicity s of λ is at least as great as the multiplicity of λ3.Proof.
Thanks to (2.15) and to the Markov property together with the fact that ℙi(Xn=0)=fn(0)i, we have
(2.24)|ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)-ℙi0(Xn1*=i1,…,Xnr*=ir)|=|1-fk-nr(0)ir1-fk(0)i0-1ρnrir·ui0·u|ℙi0(Xn1=i1,…,Xnr=ir).
The right term of (2.24) is known to converge to 0, as k→∞, thanks to the property that limkak=γu, for some γ>0, where ak:=ρ-k(1-fk(0)) (see [8]). This stems from two convergences, namely, limkbk=γ, where bk:=ρ-k(v·(1-fk(0))), and limkakbk-1=u. Let us write ak=γu+ɛk, where limkɛk=0. Since bk=v·ak and u·v=1, it comes
(2.25)akbk-1=γu+ɛkγ+v·ɛk~k→∞u+γ-1[ɛk-(v·ɛk)u].
It thus appears that the rate of convergence of ak to γu is of the same order of magnitude as the one of akbk-1 to u. Let us determine this rate in an accurate way. We use the following inequality produced by Joffe and Spitzer in [8]: for each k⩾n⩾1,
(2.26)∥akbk-1-u∥=∥1-fk(0)v·(1-fk(0))-u∥⩽2δn+∑j=k-n+1kαj1-δn-∑j=k-n+1kαj,
where we are going to replace δn and αn by some explicit formulae function of ρ and n. For this purpose, we use a detailed asymptotic behavior of Mk, as k→∞, presented for instance in [11]; we have Mk=ρkR+O(ks-1|λ|k), where R=uTv. For the sake of clarity the symbol O(·) will denote either a scalar or a matrix with all the entries satisfying the associated property. This implies the existence of some constant a>0 such that, for all k∈ℕ, (1-δk)R⩽ρ-kMk⩽(1+δk)R, where δk:=aks-1(ρ-1|λ|)k. Moreover, following [8], let us write, for all r∈[0,1]d, 1-f(r)=(M-E(r))(1-r), where 0⩽E(r)⩽M, and E(r)=O(∥1-r∥) as r→1. Then ρ-1E(fk-1(0))=ρ-1O(∥1-fk-1(0)∥)=O(ρk-2), which implies the existence of some constant b>0 such that, for all k∈ℕ, 0⩽ρ-1E(fk-1(0))⩽αkR, with αk:=bρk-2. We thus have provided an explicit formula for the sequences (δk)k and (αk)k introduced by Joffe and Spitzer in [8]. Finally let us apply (2.26) to n=⌊k/2⌋ and replace in this inequality δn and αn by their explicit expressions that we got. We obtain that, for all k∈ℕ,
(2.27)∥1-fk(0)v·(1-fk(0))-u∥⩽22-saks-1(|λ|/ρ)k/2+(b/ρ(1-ρ))ρk/21-21-saks-1(|λ|/ρ)k/2-(b/ρ(1-ρ))ρk/2.
Consequently, the right member of (2.24) will decrease with max{ks-1(|λ|ρ-1)k/2,b(ρ(1-ρ))-1ρk/2}, as k→∞.Hence the concept of theQ-process will have most practical relevance to approximate the very late extinction case if ρ is near to zero and if |λ| is small compared with ρ. Note however that the very late extinction scenario is more likely to happen if ρ is near to unity because the time to extinction in most realizations will then be long (see Section 2.1).


### 2.5. Estimation of the Infection Parameter

We assume for this subsection that the parametersΨk of the epidemic model (1.1)–(2.1) are not entirely known. More precisely, we assume that the Ψk are of the form, for all k=1⋯d,
(2.28)Ψk(θ0)=akθ0+bk,
where ak>0 and bk⩾0 are constants, and θ0 is an unknown real parameter. We will write in what follows Ψ(θ0)=aθ0+b, where Ψ(θ0):=(Ψ1(θ0),…,Ψd(θ0)) and so forth. This general assumption corresponds in particular to the case where the Ψk are of the form (1.2).We estimateθ0 by the following WCLSE in model (1.1)–(2.1). This estimator generalizes the well-known Harris estimator [12] in a BGW process. Let Θ:=]θ1,θ2[, θ2>θ1>0, such that θ0∈Θ. The WCLSE is based on the normalized process Yn:=Xn/a·Xn-1 and is defined by
(2.29)θ^∣X0∣:=argminθ∈Θ∑k=1n‍(Yk-𝔼θ(Yk∣Xk-1))2=argminθ∈Θ∑k=1n‍(Xk-Ψ(θ)·Xk-1)2a·Xk-1.
We easily derive the following explicit form:
(2.30)θ^∣X0∣=∑k=1n(Xk-b·Xk-1)∑k=1na·Xk-1.
The normalization of the process Xn by a·Xn-1 appears to be the most natural and suitable for the following reasons. First, this normalization generalizes the normalization Xn/aXn-1 in the monotype case, which is the one leading to the Harris estimator m^ of m0=aθ0+b since we have, for d=1, aθ^X0+b=m^. It also corresponds, in the linear case b=0, to the maximum likelihood estimator of θ0. In addition, defining for any vector x, x_:=minixi and x¯:=maxixi, we have
(2.31)θ0+b_a¯⩽𝔼θ0((Yk-𝔼θ0(Yk∣Xk-1))2∣Xk-1)=θ0+b·Xk-1a·Xk-1⩽θ0+b¯a_,
hence the conditional variance of the error term Yk-𝔼θ0(Yk∣Xk-1) in the stochastic regression equation Yk=𝔼θ0(Yk∣Xk-1)+Yk-𝔼θ0(Yk∣Xk-1) is invariant under multiplication of the whole process, and bounded respectively to {Xn}n, leading to the quasi-optimality of θ^∣X0∣ at finite |X0| and n, in the sense of [13].Let us provide asymptotic results for the estimatorθ^∣X0∣ defined by (2.30), as the initial population size |X0|=X0+X-1+⋯+X-d+1 tends to infinity. We denote by mij(k)(θ) the (i,j)th entry in the kth power of the matrix M(θ) given by (2.3).Theorem 2.6.
Let us assume that, for eachi=1⋯d, there exists some αi∈[0,1] such that lim∣X0∣→∞X0,i|X0|-1=αi. Then θ^∣X0∣ is strongly consistent, that is, lim∣X0∣→∞θ^∣X0∣=a.s.θ0, and is asymptotically normally distributed:
(2.32)lim∣X0∣→∞∑k=1na·Xk-1σ2(θ^∣X0∣)(θ^∣X0∣-θ0)=𝒟𝒩(0,1),
where
(2.33)σ2(θ):=θ+∑k=1n∑j=1d∑i=1dαjbimji(k-1)(θ)∑k=1n∑j=1d∑i=1dαjbimji(k-1)(θ).Proof.
Let us first prove that, for eachk=1⋯n and each i=1⋯d,
(2.34)lim∣X0∣→∞Xk,i|X0|=a.s.∑j=1d‍αjmji(k)(θ0).
Using the branching property of the process {Xn}n ([6]) we write
(2.35)Xk,i=∑j=1X0,1‍Xk,i,j(1)+⋯+∑j=1X0,d‍Xk,i,j(d),
where, for all l=1⋯d and j=1⋯X0,l, Xk,i,j(l) is the ith coordinate of a d-type branching process at time k initialized by a single particle of type l. For k, i, and l fixed the random variables {Xk,i,j(l)}j are i.i.d. with mean value mli(k)(θ0). According to the strong law of large numbers and under the theorem assumption, we have, for every l=1⋯d such that X0,l≠0,
(2.36)lim|X0|→∞∑j=1X0,lXk,i,j(l)X0,l=a.s.mli(k)(θ0),
which together with the theorem assumption leads to (2.34).
To prove the consistency ofθ^∣X0∣ we apply (2.34) to (2.30), using the fact that Xk=Xk,1 and Xk-i=Xk-1,i, and obtain
(2.37)lim∣X0∣→∞θ^∣X0∣=a.s.∑k=1n∑j=1dαj(mj1(k)(θ0)-∑i=1dbimji(k-1)(θ0))∑k=1n∑i=1d∑j=1daiαjmji(k-1)(θ0).
By definition,
(2.38)mj1(k)(θ0)=∑i=1d‍mji(k-1)(θ0)mi1(θ0)=∑i=1d‍mji(k-1)(θ0)(aiθ0+bi),
hence (2.37) immediately leads to the strong consistency.
We are now interested in the asymptotic distribution ofθ^∣X0∣-θ0. We derive from (2.30) that
(2.39)∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=∑k=1n(Xk-Ψ(θ0)·Xk-1)∑k=1na·Xk-1.
By (1.1),
(2.40)Xk-Ψ(θ0)·Xk-1=∑i=1d‍∑j=1Xk-i‍(Yk-i,k,j-Ψi(θ0))=:∑i=1d‍∑j=1Xk-i‍Y°k-i,k,j,
where the {Yk-i,k,j}j are i.i.d. given ℱk-1, following a Poisson distribution with parameter Ψi(θ0), and the {Yk-i,k,j}i,j are independent given ℱk-1. Renumbering the Y°k-i,k,j we then obtain
(2.41)∑k=1n‍(Xk-Ψ(θ0)·Xk-1)=∑i=1d‍∑j=1∑k=1nXk-i‍Y°k-i,k,j.
Applying a central limit theorem for the sum of a random number of independent random variables (see, e.g., [14]), we obtain that, for all i=1⋯d,
(2.42)lim∣X0∣→∞∑j=1∑k=1nXk-iY°k-i,k,j∑k=1nXk-i=𝒟𝒩(0,aiθ0+bi).
We have used the fact that |X0| is a real positive sequence growing to infinity, and ∑k=1n‍Xk-i is a sequence of integer-valued random variables such that ∑k=1n‍Xk-i/|X0| converges in probability to a finite random variable.  In our case the limit is actually deterministic, since we have shown in (2.34) that
(2.43)lim|X0|→∞∑k=1nXk-i|X0|=a.s.∑k=1n‍∑j=1d‍αjmji(k-1)(θ0).
Using (2.41) in (2.39), we write
(2.44)∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=∑i=1d‍∑j=1∑k=1nXk-iζ°k-i,k,j∑k=1nXk-i∑k=1nXk-i∑k=1na·Xk-1.
Using again (2.34),
(2.45)lim∣X0∣→∞∑k=1na·Xk-i∑k=1na·Xk-1=a.s.∑k=1n∑j=1dαjmji(k-1)(θ0)∑k=1n∑j=1d∑l=1dαjalmjl(k-1)(θ0),
which, combined to (2.42) and (2.44), implies by Slutsky's theorem that
(2.46)lim∣X0∣→∞∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=𝒟𝒩(0,σ2(θ0)).
By (2.33) and the strong consistency, lim∣X0∣σ2(θ0)/σ2(θ^∣X0∣)=a.s.1, from which we finally deduce (2.32).


### 2.6. Estimation of the Infection Parameter in the Worst-Case Scenario

In order to make predictions of the evolution of the epidemic in case of a very late extinction, that is, in order to make predictions of the behavior of the conditioned process{Xn*}n introduced in Section 2.4, we need to estimate the parameter θ0 in the setting of this conditioned process. We point out that θ0 does not play the same role in the conditioned process {Xn*}n and in the unconditioned process {Xn}n, since, as shown in Proposition 2.2, this parameter interferes not only in the Poisson random variable but also in the Bernoulli one. It would thus be irrelevant to estimate θ0 with an estimator aimed for the unconditioned process, such as θ^∣X0∣. Let us notice that, according to (2.16), the process {Xn*}n could be written as a multitype branching process with state dependent immigration. Because of this state-dependence, and since the parameter θ0 acts in a nonlinear way in the immigration, the methods developed in estimation theory for branching processes with immigration (see, e.g., [15]) cannot be directly applied here. Similarly as in Section 2.5 we consider the WCLSE based on the process Yn*:=Xn*/a·Xn-1*, namely,
(2.47)θ^n*:=argminθ∈ΘSn(θ),Sn(θ):=∑k=1n‍(Yk*-f(θ,Xk-1*))2,
where Θ is defined in Section 2.5, and where
(2.48)f(θ0,Xk-1*):=𝔼θ0(Yk*∣Xk-1*)=Xk-1*·Ψ(θ0)+p(θ0,Xk-1*)a·Xk-1*.
Let ɛk*:=Yk*-f(θ0,Xk-1*) be the error term between the normalized process and its conditional expectation. We obtain that
(2.49)g(θ0,Xk-1*):=𝔼θ0((ɛk*)2∣Xk-1*)=Xk-1*·Ψ(θ0)+p(θ0,Xk-1*)(1-p(θ0,Xk-1*))a·Xk-1*,
which implies
(2.50)θ0+b_a¯⩽𝔼θ0((ɛk*)2∣Xk-1*)⩽θ0+b¯+1a_.
In what follows, we denote by f' the derivative of f with respect to θ, and similarly for the other quantities depending on θ.Theorem 2.7.
The estimatorθ^n* is strongly consistent, that is, limn→∞θ^n*=a.s.θ0, and has the following asymptotic distribution:
(2.51)limn→∞∑k=0nf′(θ^n*,Xk*)2∑k=0nf′(θ^n*,Xk*)2g(θ^n*,Xk*)(θ^n*-θ0)=𝒟𝒩(0,1).Remark 2.8.
Note that (2.51) involves the function f' and thus requires the knowledge of the derivative of the function uj given in (2.5), which is not an explicit function of θ since ρ is not either. However, ρ' satisfies ρ'=∑k=1d‍akρ-k[∑k=1d‍k(akθ+bk)ρ-k-1]-1; hence uj′ is known as soon as ρ can be computed. Consequently, denoting ρ by ρ(θ) when θ is the parameter of the model, (2.51) can be used as soon as ρ(θ^n*obs) is known; for this purpose one can for instance numerically approximate the largest solution ρ of ∑k=1d‍Ψk(θ^n*obs)ρ-k=1.Proof.
The proof heavily relies on a strong law of large numbers for homogeneous irreducible positive recurrent Markov chains applied to the conditioned process{Xn*}n and its stationary distribution πθ0 (given by Proposition 2.4), which states [16]: for every πθ0-integrable function h:ℕd∖{0}→ℝ,
(2.52)limn→∞1n∑k=0n-1‍h(Xk*)=a.s.∑j∈ℕd‍h(j)πθ0(j).
Note that the Perron's root ρ(θ) of M(θ) and the associated right normalized eigenvector u(θ) are 𝒞∞-functions of θ.
In order to prove the strong consistency ofθ^n* as n→∞, since f′′(·,Xk-1*) is not linear, we cannot use the general standard method based on the first-order expansion of Sn′(θ^n*) at θ0, and on the strong law of large numbers for martingales applied to Sn′(θ0) correctly normalized. We consequently use the following conditions given in [17]:(i)
f(.,Xk-1*) is Lipschitz on Θ; that is, there exists a nonnegative ℱk-1*-measurable function Ak (where ℱk-1*:=σ(X0*,…,Xk-1*)), satisfying, for all δ1, δ2∈Θ, |f(δ1,Xk-1*)-f(δ2,Xk-1*)|⩽Ak|δ1-δ2| a.s.,(ii)
lim¯k→∞𝔼θ0((ɛk*)2∣Xk-1*)<a.s.∞,(iii)
limn→∞infθ∈Θ∣θ-θ0∣⩾δ∑k=1n‍(f(θ0,Xk-1*)-f(θ,Xk-1*))2=a.s.∞. Let us note that, in the frame of a general model f, since Θ¯ is compact, (i) is satisfied as soon as f(θ,Xk-1*) has a first derivative in θ with supθ∈Θf'(θ,Xk-1*)<∞, (ii) is satisfied for any optimal estimator in the sense of [13] and moreover could be weakened (see [17]), and (iii) is a necessary condition. First, for all θ∈Θ and j∈ℕd, j≠0, f'(θ,j)=(a·j+p'(θ,j))(a·j)-1/2, where p' denotes the derivative of p with respect to θ, which thanks to (2.17) is bounded on Θ. Condition (i) is thus satisfied. Condition (ii) follows from (2.50). Condition (iii) comes from the fact that, for every δ>0 and every θ∈Θ such that ∣θ-θ0∣⩾δ, applying the mean value theorem to the 𝒞1-function p(·,Xk-1*),
(2.53)∑k=1n‍(f(θ0,Xk-1*)-f(θ,Xk-1*))2=(θ0-θ)2∑k=1n‍a·Xk-1*(1+p(θ0,Xk-1*)-p(θ,Xk-1*)(θ0-θ)a·Xk-1*)2⩾δ2∑k=1n‍a·Xk-1*infθ∈Θ¯(1+p'(θ,Xk-1*)a·Xk-1*)2.
Let us show that the function j↦a·jinfθ∈Θ¯(1+p'(θ,j)/a·j)2 is πθ0-integrable. For every θ∈Θ¯, j∈ℕd and j≠0, denoting Ξ:=supi,θ∈Θ¯{ui(θ),∣ui′(θ)∣}<∞ (by continuity of ui′(θ) on Θγ⊃Θ¯), and u_:=mini,θ∈Θ¯ui(θ)>0,
(2.54)|p'(θ,j)|=|(ui′(θ)j·Ψ(θ)+u1(θ)j·a)∑i=2djiui(θ)-u1(θ)j·Ψ(θ)∑i=2djiui′(θ)(u1(θ)j·Ψ(θ)+∑i=2djiui(θ))2|⩽Ξ2u_23j·Ψ(θ)∑i=2dji(j·Ψ(θ)+∑i=2dji)2⩽3Ξ24u_2=:C1.
Hence, for all j≠0, ∣a·jinfθ∈Θ¯(1+p'(θ,j)/a·j)2∣⩽(1+C1/a_)2a·j, and applying (2.52) together with (2.22) we obtain that
(2.55)limn→∞1n+1∑k=1n‍a·Xk-1*infθ∈Θ¯(1+p′(θ,Xk-1*)a·Xk-1*)2=a.s.∑j∈ℕd‍a·jinfθ∈Θ¯(1+p′(θ,j)a·j)2πθ0(j).
Let j≠0 fixed. Since, for all θ∈Θ¯, p'(θ,j)≠-a·j, the extreme value theorem implies that infθ∈Θ¯(1+p'(θ,j)/a·j)2>0. Hence the right term in (2.55) is strictly positive, which together with (2.55) leads to (iii).
Let us now consider the asymptotic distribution ofθ^n*-θ0. For this purpose, we follow the steps of the proof of Proposition 6.1 in [17]. Writing the Taylor expansion of Sn′(θ^n*) at θ0 we obtain that θ^n*-θ0=-Sn′(θ0)/Sn′′(θ~n), for some θ~n=θ0+tn(θ^n*-θ0), with tn∈]0,1[. Since Sn′(θ0)=-2∑k=1n‍ɛk*f'(θ0,Xk-1*), we can write
(2.56)n(θ^n*-θ0)=∑k=1nɛk*f′(θ0,Xk-1*)n(Fnn)-1(12Sn′′(θ~n)Fn)-1,
where Fn:=∑k=1n‍(f'(θ0,Xk-1*))2. Let us first show that
(2.57)limn→∞Fnn=a.s.∑j∈ℕd‍(f′(θ0,j))2πθ0(j).
This is an application of (2.52) and (2.22), since, for all j∈ℕd, j≠0,
(2.58)(f′(θ0,j))2⩽a·j+2C1+C12a_.
In view of (2.56), we now prove that
(2.59)limn→∞Sn′′(θ~n)Fn=a.s.2.
Computing Sn′′ thanks to the formula Sn(θ)=∑k=1n‍(ɛk*+f(θ0,Xk-1*)-f(θ,Xk-1*))2, it appears that (2.59) is true, as soon as the following holds:
(2.60)limn→∞supθ∈Θ|∑k=1nɛk*f′′(θ,Xk-1*)|Fn=a.s.0,(2.61)limn→∞∑k=1n(f′(θ~n,Xk-1*))2Fn=a.s.1,(2.62)limn→∞∑k=1n(f(θ0,Xk-1*)-f(θ~n,Xk-1*))f′′(θ~n,Xk-1*)Fn=a.s.0.
Let us prove (2.60)–(2.62). First, (2.60) is given by a strong law of large numbers proved in [17, Proposition 5.1]. The latter can be indeed applied since limnFn=a.s.∞ (as an immediate consequence of the stronger result (2.57)), and since f′′(·,Xk-1*) fulfills the required Lipschitz condition. Indeed we proved earlier that ui′′′(θ) is continuous on the compact set Θ¯ and is thus bounded on Θ¯, which implies that f′′′(·,Xk-1*)=p′′′(·,Xk-1*)(a·Xk-1*)-1/2 is bounded by a ℱk-1*-measurable function. In view of (2.61), we consider the function f(θ,j)2 and its derivative 2f'(θ,j)f′′(θ,j). Similarly as for (2.54), one can show that there exists a constant C2>0 such that for all θ∈Θ¯, and all j≠0, ∣p′′(θ,j)∣⩽C2. This implies
(2.63)|2f′(θ,j)f′′(θ,j)|⩽2|(a·j+p′(θ,j))p′′(θ,j)a·j|⩽2C2(1+C1a_).
Consequently,
(2.64)|∑k=1n(f′(θ~n,Xk-1*))2-(f′(θ0,Xk-1*))2|Fn⩽2C2(1+C1a_)|θ^n*-θ0|(Fnn)-1,
which by (2.57) and the strong consistency of θ^n* almost surely tends to 0. Writing
(2.65)∑k=1nf(θ~n,Xk-1*)2Fn=1+∑k=1n((f′(θ~n,Xk-1*))2-(f′(θ0,Xk-1*))2)Fn,
this implies (2.61). It now remains to prove (2.62). We write
(2.66)|∑k=1n(f(θ0,Xk-1*)-f(θ~n,Xk-1*))f′′(θ~n,Xk-1*)|Fn⩽1Fn∑k=1n‍a·Xk-1*|θ0-θ~n|+|p(θ0,Xk-1*)-p(θ~n,Xk-1*)|a·Xk-1*|p′′(θ~n,Xk-1*)|⩽1Fn∑k=1n‍(|θ0-θ~n|+C1|θ0-θ~n|a_)C2⩽|θ0-θ^n*|C2(1+C1a_)(Fnn)-1,
which thanks to (2.57) and the strong consistency of θ^n* implies (2.62). In view of (2.56), we finally want to prove that ∑k=1n‍ɛk*f'(θ0,Xk-1*)/n converges in distribution. For this purpose we make use of a central limit theorem for martingale difference arrays [18, 19]; if {Mk(n),ℱk(n),1⩽k⩽n}, n⩾1 is a sequence of square integrable martingales with associated Meyer process 〈M〉(n)=(〈M〉k(n))1⩽k⩽n satisfying limn→∞〈Mn〉(n)=Pc2 for some constant c, and such that, for all ɛ>0,
(2.67)limn→∞∑k=1n‍𝔼[|Mk(n)-Mk-1(n)|21{|Mk(n)-Mk-1(n)|⩾ɛ}∣ℱk-1(n)]=P0,
then limn→∞Mn(n)=𝒟𝒩(0,c2). Let us define Mk(n):=∑l=1k‍ɛl*f'(θ0,Xl-1*)/n, for every k⩽n. First, for any k⩽n, 𝔼θ0(ɛk*f'(θ0,Xk-1*)/n∣Xk-1*)=0. Second,
(2.68)𝔼θ0((ɛk*f′(θ0,Xk-1*)n)2∣Xk-1*)=(f′(θ0,Xk-1*))2g(θ0,Xk-1*)n,
hence {Mk(n)}k⩽n is a sequence of square integrable martingales. Moreover, using inequalities (2.50) and (2.58), we obtain, by (2.22),
(2.69)∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j)⩽(θ0+b¯+1a_)(∑j∈ℕd‍a·jπθ0(j)+2C1+C12a_)<∞.
So, by means of (2.52),
(2.70)limn→∞〈Mn〉(n)=limn→∞∑k=1n‍𝔼θ0((ɛkf′(θ0,Xk-1*)n)2∣Xk-1*)=a.s.∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j).
Third, using Cauchy-Schwarz and Bienaymé-Chebyshev inequalities,
(2.71)∑k=1n‍𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|21{|ɛk*f′(θ0,Xk-1*)/n|⩾ɛ}∣Xk-1*]⩽∑k=1n‍(𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|4∣Xk-1*])1/2(ℙθ0[|ɛk*f′(θ0,Xk-1*)n|⩾ɛ∣Xk-1*])1/2⩽1n3/2ɛ∑k=1n‍|f′(θ0,Xk-1*)|3(𝔼θ0[(ɛk*)4∣Xk-1*])1/2(𝔼θ0[(ɛk*)2∣Xk-1*])1/2.
Let us compute 𝔼θ0((ɛk*)4∣Xk-1*). We can show that the 4th central moment of the independent sum of a Poisson and a Bernoulli random variables equals μ4+6μ2γ2+γ4, where μi and γi denote the ith central moment of the Poisson and of the Bernoulli variable. If these variables have parameter λ and p, respectively, then μ4=λ(1+3λ), μ2=λ, γ4=p(1-p)(3p2-3p+1)∈[0,1], and γ2=p(1-p)∈[0,1]. We thus obtain
(2.72)|𝔼θ0((ɛk*)4∣Xk-1*)|⩽Ψ(θ0)·Xk-1*(7+3Ψ(θ0)·Xk-1*)+1(a·Xk-1*)2.
Hence
(2.73)|f*′(θ0,Xk-1*)|3(𝔼θ0[ɛk4∣Xk-1*])1/2(𝔼θ0[ɛk2∣Xk-1*])1/2⩽(a·Xk-1*+C1a_)3Ψ(θ0)·Xk-1*(7+3Ψ(θ0)·Xk-1*)+1a·Xk-1*(θ0+b¯+1a_)1/2.
Since the highest power of Xn* involved in (2.73) is 3/2, and since by Proposition 2.4 the stationary distribution πθ0 has finite second moments, we can apply (2.52) to (2.71) and obtain that
(2.74)limn→∞∑k=1n‍𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|21{|ɛk*f′(θ0,Xk-1*)/n|⩾ɛ}∣Xk-1*]=a.s.0.
It then ensues from the central limit theorem mentioned above that
(2.75)limn→∞∑k=1nɛk*f′(θ0,Xk-1*)n=𝒟𝒩(0,∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j)).
Finally, (2.56), together with (2.57), (2.59), (2.75), and Slutsky's theorem, implies that
(2.76)limn→∞n(θ^n*-θ0)=𝒟𝒩(0,∑j∈ℕd(f′(θ0,j))2g(θ0,j)πθ0(j)(∑j∈ℕd(f′(θ0,j))2πθ0(j))2),
which leads to
(2.77)limn→∞n(∑j∈ℕdf′(θ0,j)2πθ0(j))2∑j∈ℕdf′(θ0,j)2g(θ0,j)πθ0(j)(θ^n*-θ0)=𝒟𝒩(0,1).
It now remains to prove the asymptotic distribution (2.51). Thanks to what precedes, this result is immediate as soon as we prove that
(2.78)limn→∞1n+1∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)=a.s.∑j∈ℕd‍(f'(θ0,j))2g(θ0,j)πθ0(j),
as well as the equivalent result for the numerator. For this purpose, we write
(2.79)∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)=∑k=0n‍(f′(θ0,Xk*))2g(θ0,Xk*)+∑k=0n‍[(f′(θ^n*,Xk*))2g(θ^n*,Xk*)-(f′(θ0,Xk*))2g(θ0,Xk*)],
and show that (f'(·,j))2g(·,j) has a bounded derivative and is thus Lipschitz. We have indeed ∣g'(θ,j)∣=∣a·j+p'(θ,j)∣(a·j)-1⩽1+C1a_-1, hence
(2.80)|2f′′(θ,j)f′(θ,j)g(θ,j)+(f′(θ,j))2g′(θ,j)|⩽2C2(1+C1a_-1)a_-1(Ψ(θ)·j+1)+(1+C1a_-1)2(a·j+C1)⩽C3Ψ(θ2)·j,
for some constant C3>0. This enables us to write
(2.81)1n+1∑k=0n‍|(f′(θ^n*,Xk*))2g(θ^n*,Xk*)-(f′(θ0,Xk*))2g(θ0,Xk*)|⩽|θ^n*-θ0|C31n+1∑k=0n‍Ψ(θ2)·Xk*.
By the strong consistency of θ^n* together with (2.52) and (2.22), (2.81) almost surely tends to zero. Combined with (2.52) and (2.69) in (2.79), this implies (2.78).


## 2.1. Extinction of the Epidemic

### 2.1.1. Almost Sure Extinction

Since the single-type process{Xn}n has a memory of size d, it becomes extinct when it is null at d successive times, or equivalently as soon as the d-dimensional process {Xn}n reaches the d-dimensional null vector 0. According to the theory of multitype positive regular and nonsingular BGW processes [6], the extinction of the process {Xn}n occurs almost surely (a.s.), if and only if ρ⩽1, where ρ is the dominant eigenvalue (also called the Perron's root) of the mean matrix M. Thus ρ is solution of ∑k=1d‍Ψkρ-k=1. In general for d>1, ρ has no explicit expression. However, ∑k=1d‍Ψkρ-k=1 leads directly to the following explicit threshold criteria.Proposition 2.1.
The epidemic becomes extinct almost surely if and only ifR0⩽1, where R0:=∑k=1d‍Ψk is the total mean number of secondary cases generated by one infective in a SEIR disease. We call R0 the basic reproduction number.
Moreover, whenρ⩽1, then R0⩽ρ with equality if and only if either ρ=1 or d=1, and when ρ>1, then R0⩾ρ with equality if and only if d=1.Note that whend>1, R0 only provides information about the threshold level, whereas ρ provides an additional information about the speed of extinction of the process, as shown in the next two paragraphs.


### 2.1.2. Speed of Extinction

Thanks to well-known results in the literature about multitype branching processes and more particularly to the Perron-Frobenius theorem (see, e.g., [6]), we can deduce the expected incidence of infectives in the population at time n, for n large. Denoting by u and v the right and left eigenvectors of M associated to the Perron's root ρ; that is, MuT=ρuT and vM=ρv, with the normalization convention u·1=u·v=1, where u·v stands for the usual scalar product in ℝd and where the superscript T denotes the transposition, then 𝔼(Xn∣X0)=X0Mn~n→∞ρnX0uTv. The first coordinate in the latter formula becomes for the epidemic process ρn∑i=1d‍X-i+1uiv1. Computing explicitly u and v, we obtain that for all i=1⋯d,
(2.5)ui=∑k=idΨkρ-k+i∑j=1d∑k=jdΨkρ-k+j,vi=ρ-i∑j=1d∑k=jdΨkρ-k+j∑j=1d∑k=jdΨkρ-k,
which leads to the following asymptotic result:
(2.6)𝔼(Xn∣X0)~n→∞ρn∑i=1d‍X-i+1∑k=idΨkρ-k+i-1∑j=1d∑k=jdΨkρ-k.
Hence if ρ<1, the mean number of infectives decreases exponentially at the rate ρ. In the following section, we provide a much finer result on the estimation of the disease extinction time in the population.


### 2.1.3. Extinction Time of the Epidemic

The extinction time distribution can be derived as a function of the offspring generating function. As usual in stochastic processes, this quantity is calculated conditionally on the initial valueX0=(X0,X-1,…,X-d+1), but for the sake of simplicity we do not let it appear in the notations. Note that since we are building tools for the prediction of the spread of the disease, the time origin 0 corresponds here to the time of the last available data (generally the current date). Let T:=inf{n⩾1,Xn=0} denote the extinction time of the process {Xn}n, and let fn:=f∘fn-1 be the nth iterate of the generating function f given by (2.2). We denote fn:=(fn,1,…,fn,d). Then, by the branching property of the process [6], the probability of extinction of the epidemic before time n is immediately given by ℙ(T⩽n)=ℙ(Xn=0)=fn(0)X0, that is to say,
(2.7)ℙ(T⩽n)=(fn,1(0))X0(fn,2(0))X-1⋯(fn,d(0))X-d+1.
It can be immediately deduced from convergence results for fn(0) as n→∞ [8], that if ρ=1, ℙ(T⩽n)~1-(nη)-1X0·u, while if ρ<1, ℙ(T⩽n)~1-ρnγX0·u, for some constants η,γ>0. As a consequence, the closer ρ is to unity, the longer the time to extinction will be in most realizations. More specifically, (2.7) enables the exact computation (resp., estimation) of ℙ(T⩽n) for any n by the iterative computation of fn, X0 being given, when the parameters Ψk of (1.1) are known (resp., estimated). Moreover, since for ρ⩽1 the epidemic becomes extinct in an a.s. finite time and ℙ(T⩽n)↗n→∞1, then for any given probability p∈[0,1[ there exists n∈ℕ such that ℙ(T⩽n)⩾p. So in practice, for any p∈[0,1[, (2.7) enables us to compute the p-quantile npT of the extinction time,
(2.8)npT:=min{n⩾1:ℙ(T⩽n)⩾p}.


## 2.1.1. Almost Sure Extinction

Since the single-type process{Xn}n has a memory of size d, it becomes extinct when it is null at d successive times, or equivalently as soon as the d-dimensional process {Xn}n reaches the d-dimensional null vector 0. According to the theory of multitype positive regular and nonsingular BGW processes [6], the extinction of the process {Xn}n occurs almost surely (a.s.), if and only if ρ⩽1, where ρ is the dominant eigenvalue (also called the Perron's root) of the mean matrix M. Thus ρ is solution of ∑k=1d‍Ψkρ-k=1. In general for d>1, ρ has no explicit expression. However, ∑k=1d‍Ψkρ-k=1 leads directly to the following explicit threshold criteria.Proposition 2.1.
The epidemic becomes extinct almost surely if and only ifR0⩽1, where R0:=∑k=1d‍Ψk is the total mean number of secondary cases generated by one infective in a SEIR disease. We call R0 the basic reproduction number.
Moreover, whenρ⩽1, then R0⩽ρ with equality if and only if either ρ=1 or d=1, and when ρ>1, then R0⩾ρ with equality if and only if d=1.Note that whend>1, R0 only provides information about the threshold level, whereas ρ provides an additional information about the speed of extinction of the process, as shown in the next two paragraphs.


## 2.1.2. Speed of Extinction

Thanks to well-known results in the literature about multitype branching processes and more particularly to the Perron-Frobenius theorem (see, e.g., [6]), we can deduce the expected incidence of infectives in the population at time n, for n large. Denoting by u and v the right and left eigenvectors of M associated to the Perron's root ρ; that is, MuT=ρuT and vM=ρv, with the normalization convention u·1=u·v=1, where u·v stands for the usual scalar product in ℝd and where the superscript T denotes the transposition, then 𝔼(Xn∣X0)=X0Mn~n→∞ρnX0uTv. The first coordinate in the latter formula becomes for the epidemic process ρn∑i=1d‍X-i+1uiv1. Computing explicitly u and v, we obtain that for all i=1⋯d,
(2.5)ui=∑k=idΨkρ-k+i∑j=1d∑k=jdΨkρ-k+j,vi=ρ-i∑j=1d∑k=jdΨkρ-k+j∑j=1d∑k=jdΨkρ-k,
which leads to the following asymptotic result:
(2.6)𝔼(Xn∣X0)~n→∞ρn∑i=1d‍X-i+1∑k=idΨkρ-k+i-1∑j=1d∑k=jdΨkρ-k.
Hence if ρ<1, the mean number of infectives decreases exponentially at the rate ρ. In the following section, we provide a much finer result on the estimation of the disease extinction time in the population.


## 2.1.3. Extinction Time of the Epidemic

The extinction time distribution can be derived as a function of the offspring generating function. As usual in stochastic processes, this quantity is calculated conditionally on the initial valueX0=(X0,X-1,…,X-d+1), but for the sake of simplicity we do not let it appear in the notations. Note that since we are building tools for the prediction of the spread of the disease, the time origin 0 corresponds here to the time of the last available data (generally the current date). Let T:=inf{n⩾1,Xn=0} denote the extinction time of the process {Xn}n, and let fn:=f∘fn-1 be the nth iterate of the generating function f given by (2.2). We denote fn:=(fn,1,…,fn,d). Then, by the branching property of the process [6], the probability of extinction of the epidemic before time n is immediately given by ℙ(T⩽n)=ℙ(Xn=0)=fn(0)X0, that is to say,
(2.7)ℙ(T⩽n)=(fn,1(0))X0(fn,2(0))X-1⋯(fn,d(0))X-d+1.
It can be immediately deduced from convergence results for fn(0) as n→∞ [8], that if ρ=1, ℙ(T⩽n)~1-(nη)-1X0·u, while if ρ<1, ℙ(T⩽n)~1-ρnγX0·u, for some constants η,γ>0. As a consequence, the closer ρ is to unity, the longer the time to extinction will be in most realizations. More specifically, (2.7) enables the exact computation (resp., estimation) of ℙ(T⩽n) for any n by the iterative computation of fn, X0 being given, when the parameters Ψk of (1.1) are known (resp., estimated). Moreover, since for ρ⩽1 the epidemic becomes extinct in an a.s. finite time and ℙ(T⩽n)↗n→∞1, then for any given probability p∈[0,1[ there exists n∈ℕ such that ℙ(T⩽n)⩾p. So in practice, for any p∈[0,1[, (2.7) enables us to compute the p-quantile npT of the extinction time,
(2.8)npT:=min{n⩾1:ℙ(T⩽n)⩾p}.


## 2.2. Total Size of the Epidemic

Under the assumptionρ⩽1 and the independence of the {Yn-l,n,i}i,l,n (we previously assumed the independence of the {Yn-l,n,i}i,l, for each n), we derive the distribution of the future total size N:=∑n=1T‍Xn of the epidemic until its extinction, that is, the future total number of infectives until the extinction of the disease. It can be shown [9] that, given the initial value X0, the time origin being the same as in Section 2.1, the probability distribution of N is
(2.9)N=𝒟∑k=1d‍∑i=1X-k+1‍(∑j=1Yk,i‍Nk,i,j),Yk,i:=∑l=kd‍Y-k+1,-k+1+l,i,
where =𝒟 denotes the equality in distribution, an empty sum is by convention 0, the {Yk,i}k,i are independent, and the {Yk,i}i and the {Nk,i,j}k,i,j are i.i.d. with
(2.10)Yk,i=𝒟Poiss(∑l=kd‍Ψl),Nk,i,j=𝒟Borel-Tanner(∑l=1d‍Ψl,1),
that is, for each n⩾1, ℙ(Nk,i,j=n)=e-n∑l=1d‍Ψl(n∑l=1d‍Ψl)n-1(n!)-1. Consequently, under the convention that an empty product is 1, the probability distribution of N is, for any n∈ℕ,
(2.11)ℙ(N=n)=∑∑k=1d∑i=1X-k+1∑j=1yk,ink,i,j=n{0⩽yk,i⩽n,{1⩽nk,i,j⩽n}j}i,k:‍∏k=1d‍∏i=1X-k+1‍e-∑l=kdΨl(∑l=kdΨl)yk,iyk,i!×∏j=1yk,i‍e-nk,i,j∑l=1dΨl(nk,i,j∑l=1dΨl)nk,i,j-1nk,i,j!,
which may be calculated (resp., estimated), replacing the Ψk by their values (resp., estimations). In practice, for any p∈[0,1[, (2.11) enables to compute the p-quantile npN of the total epidemic size,
(2.12)npN:=min{n⩾1:ℙ(N⩽n)⩾p}.
We obtain moreover an explicit formula for the mean value and variance of the size of the epidemic,
(2.13)𝔼(N)=∑k=1dX-k+1∑l=kdΨl1-∑l=1dΨl,𝕍ar(N)=∑k=1dX-k+1∑l=kdΨl(1-∑l=1dΨl)3.


## 2.3. Exposed Population

Depending on the disease, it might also be crucial to study and predict the evolution of the incidence of exposed individuals in the population, which is generally unobservable. We assume that this information is given by the process{Zn}n defined by the conditional distribution Zn∣Xn=𝒟Poiss(Ψ0Xn), where Ψ0 is the mean number of individuals infected at time n by an infective of this time (see Section 3.1). This property enables on the one hand to reconstruct the whole past epidemic (i.e., the incidence of infectives as well as of exposed individuals) thanks to the observable data. On the other hand, it allows to simulate the evolution of the incidence of exposed individuals in the future, based on predictions of the evolution of the epidemic process {Xn}n.


## 2.4. Worst-Case Scenario: Very Late Extinction of the Epidemic

Even in the case when the epidemic dies out almost surely (ρ⩽1), and although one can provide the p-quantile npT of the extinction time with the probability p as large as wanted (see (2.8)), the epidemic might become extinct after this time with a small but nonnull probability of order 1-p. This raises the following question: how would the incidences of infectious and exposed individuals evolve in the (unlikely) case of a very late extinction? In terms of risk analysis, this issue appears to be crucial to evaluate the risks associated with this worst-case scenario. The tools developed in the previous subsections allow to evaluate the probability of all possible outcomes. But since the worst ones, typically a very late extinction, have a negligible probability, these tools do not bring any information in these worst cases and in particular do not inform on the evolution at each time step of the spread of the disease (would it decrease extremely slowly, stay at a constant rate for a very long time, present several peaks in its evolution, etc.?). In order to study the propagation of the epidemic in the decay phase, assuming that extinction occurs very late, we are interested in the distribution of the process {Xn}n conditionally on the event that the epidemic has not become extinct at time k, where k is very large. We therefore consider for any n1,n2,…∈ℕ and any i0,i1,i2,…∈ℕd the conditioned probability ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0), where the subscript i0 denotes the initial value. If k is finite this distribution cannot be easily handled due to its time inhomogeneity. However, when ρ⩽1, it is known ([10]) that this conditioned distribution converges, as k→∞, to the distribution of a d-dimensional Markov process {Xn*}n:
(2.14)limk→∞ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)=ℙi0(Xn1*=i1,…,Xnr*=ir).
We will further discuss in Proposition 2.5 the relevancy of approximating the conditioned probability for k fixed by the limiting object (2.14). The conditioned process {Xn*}n defined by (2.14) is known in the literature as the Q-process associated with {Xn}n, also described as the process conditioned on “not being extinct in the distant future.” It has the following transition probability [10]: for every n⩾1, i,j∈ℕd, i≠0,
(2.15)ℙ(Xn*=j∣Xn-1*=i)=1ρj·ui·uℙ(Xn=j∣Xn-1=i),
where u is the normalized right eigenvector of M associated to the Perron's root ρ as introduced in Section 2.1, and computed explicitly in (2.5). In the same way as for the process {Xn}n, we define the 1-dimensional process Xn*:=Xn,1*. By construction we then have Xn,i*=Xn-i+1*, for each n and each i=1⋯d.Proposition 2.2.
The stochastic process{Xn*}n is, conditionally on its past, distributed as the sum of two independent Poisson and Bernoulli random variables:
(2.16)Xn*∣Xn-1*=𝒟Poiss(Xn-1*·Ψ)*ℬ(p(Xn-1*)),
where Ψ:=(Ψ1,…,Ψd), * is the convolution product symbol, and
(2.17)p(Xn-1*):=u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk.Proof.
Applying (1.1) and (2.15), we obtain that for all j∈ℕ,
(2.18)ℙ(Xn*=j∣Xn-1*)=ℙ(Xn*=(j,Xn-1*,…,Xn-(d-1)*)∣Xn-1*)=ju1+∑k=2dXn-k+1*ukρXn-1*·uℙ(Xn=(j,Xn-1*,…,Xn-(d-1)*)∣Xn-1=Xn-1*)=ju1+∑k=2dXn-k+1*ukρXn-1*·u(Xn-1*·Ψ)jj!e-Xn-1*·Ψ=u1Xn-1*·ΨρXn-1*·u(Xn-1*·Ψ)j-1(j-1)!e-Xn-1*·Ψ+∑k=2dXn-k+1*ukρXn-1*·u(Xn-1*·Ψ)jj!e-Xn-1*·Ψ.
The equality MuT=ρuT implies that for all k=1⋯d-1, ρuk=Ψku1+uk+1, and that ρud=Ψdu1. Consequently, ρXn-1*·u=u1Xn-1*·Ψ+∑k=2d‍Xn-k+1*uk, and thus
(2.19)ℙ(Xn*=j∣Xn-1*)=u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk(Xn-1*·Ψ)j-1(j-1)!e-Xn-1*·Ψ+(1-u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk)(Xn-1*·Ψ)jj!e-Xn-1*·Ψ=[Poiss(Xn-1*·Ψ)*ℬ(u1Xn-1*·Ψu1Xn-1*·Ψ+∑k=2dXn-k+1*uk)](j).Remark 2.3.
Note that if one compares (2.16) with the transition probability of the unconditioned process Xn∣Xn-1=𝒟Poiss(Xn-1·Ψ), it appears that {Xn*}n behaves at each time step like {Xn}n, according to a Poisson distribution, except that it has the possibility at each time step to add one unit or not, according to a Bernoulli random variable. Moreover, if Xn-1*=⋯=Xn-(d-1)*=0, then according to (2.17), p(Xn-1*)=1, which implies that at time n, the probability to add one unit is equal to one, thus preventing the extinction of the process.Proposition 2.4.
The process{Xn*}n admits a stationary probability measure π with finite first- and second-order moments.Proof.
Since the multitype branching process{Xn}n satisfies property (2.4), it is known [10] that the Q-process {Xn*}n is positive recurrent with a stationary probability measure π given by
(2.20)π(i)≔i·uν(i)∑k∈ℕdk·uν(k),i∈ℕd,
where ν is the Yaglom distribution of the process {Xn}n, uniquely defined by the following property: for all i,j∈ℕd∖{0}, limn→∞ℙ(Xn=i∣X0=j,Xn≠0)=ν(i). In the literature, this stationary measure for the conditioned process {Xn*}n is also referred to as the doubly limiting conditional probability. Moreover, by Proposition 2.2,
(2.21)𝔼(Xn*)=𝔼[𝔼(Xn*∣Xn-1*)]=𝔼[Xn-1*·Ψ+p(Xn-1*)]⩽∑k=1d‍𝔼(Xn-k*)Ψk+1,
which implies that lim¯n→∞𝔼(Xn*)⩽(1-∑k=1d‍Ψk)-1<∞. We consequently obtain by means of Fatou's lemma that, for every i=1⋯d,
(2.22)∑j∈ℕd‍jiπ(j)=𝔼(limn→∞Xn,i*)=𝔼(limn→∞Xn-i+1*)=𝔼(limn→∞Xn*)⩽lim_n→∞𝔼(Xn*)<∞.
We similarly prove that π has finite second-order moments by writing
(2.23)𝕍ar(Xn*∣Xn-1*)=Xn-1*·Ψ+u1Xn-1*·Ψ∑k=2dXn-k+1*uk(u1Xn-1*·Ψ+∑k=2dXn-k+1*uk)2⩽Xn-1*·Ψ+14.Let us discuss the relevancy of approximating the epidemic process{Xn}n conditioned on nonextinction at some finite time k, for k large, by the Q-process {Xn*}n obtained by letting k→∞. When considering the case of late extinction, one works under an hypothetical assumption based on the unknown future, hence in practice one does not focus on a specific value k for the survival of the disease in the population. We therefore might consider that k is chosen large enough such that the approximation of the process {Xn}n conditioned on the event {Xk≠0} by the process {Xn*}n is valid. Of course, the order of magnitude of such k will depend on the rate of convergence of the conditioned process to {Xn*}n.Proposition 2.5.
Letn1⩽⋯⩽nr⩽k and i0,…,ir∈ℕd∖{0}. Then the difference |ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)-ℙi0(Xn1*=i1,…,Xnr*=ir)| decreases, as k→∞, with max{ks-1(|λ|ρ-1)k/2,ρk/2}, where λ is an eigenvalue of M such that ρ>|λ|⩾|λ3|⩾|λ4|⩾⋯, with the λi being the other eigenvalues of M. In case |λ|=|λ3| we stipulate that the multiplicity s of λ is at least as great as the multiplicity of λ3.Proof.
Thanks to (2.15) and to the Markov property together with the fact that ℙi(Xn=0)=fn(0)i, we have
(2.24)|ℙi0(Xn1=i1,…,Xnr=ir∣Xk≠0)-ℙi0(Xn1*=i1,…,Xnr*=ir)|=|1-fk-nr(0)ir1-fk(0)i0-1ρnrir·ui0·u|ℙi0(Xn1=i1,…,Xnr=ir).
The right term of (2.24) is known to converge to 0, as k→∞, thanks to the property that limkak=γu, for some γ>0, where ak:=ρ-k(1-fk(0)) (see [8]). This stems from two convergences, namely, limkbk=γ, where bk:=ρ-k(v·(1-fk(0))), and limkakbk-1=u. Let us write ak=γu+ɛk, where limkɛk=0. Since bk=v·ak and u·v=1, it comes
(2.25)akbk-1=γu+ɛkγ+v·ɛk~k→∞u+γ-1[ɛk-(v·ɛk)u].
It thus appears that the rate of convergence of ak to γu is of the same order of magnitude as the one of akbk-1 to u. Let us determine this rate in an accurate way. We use the following inequality produced by Joffe and Spitzer in [8]: for each k⩾n⩾1,
(2.26)∥akbk-1-u∥=∥1-fk(0)v·(1-fk(0))-u∥⩽2δn+∑j=k-n+1kαj1-δn-∑j=k-n+1kαj,
where we are going to replace δn and αn by some explicit formulae function of ρ and n. For this purpose, we use a detailed asymptotic behavior of Mk, as k→∞, presented for instance in [11]; we have Mk=ρkR+O(ks-1|λ|k), where R=uTv. For the sake of clarity the symbol O(·) will denote either a scalar or a matrix with all the entries satisfying the associated property. This implies the existence of some constant a>0 such that, for all k∈ℕ, (1-δk)R⩽ρ-kMk⩽(1+δk)R, where δk:=aks-1(ρ-1|λ|)k. Moreover, following [8], let us write, for all r∈[0,1]d, 1-f(r)=(M-E(r))(1-r), where 0⩽E(r)⩽M, and E(r)=O(∥1-r∥) as r→1. Then ρ-1E(fk-1(0))=ρ-1O(∥1-fk-1(0)∥)=O(ρk-2), which implies the existence of some constant b>0 such that, for all k∈ℕ, 0⩽ρ-1E(fk-1(0))⩽αkR, with αk:=bρk-2. We thus have provided an explicit formula for the sequences (δk)k and (αk)k introduced by Joffe and Spitzer in [8]. Finally let us apply (2.26) to n=⌊k/2⌋ and replace in this inequality δn and αn by their explicit expressions that we got. We obtain that, for all k∈ℕ,
(2.27)∥1-fk(0)v·(1-fk(0))-u∥⩽22-saks-1(|λ|/ρ)k/2+(b/ρ(1-ρ))ρk/21-21-saks-1(|λ|/ρ)k/2-(b/ρ(1-ρ))ρk/2.
Consequently, the right member of (2.24) will decrease with max{ks-1(|λ|ρ-1)k/2,b(ρ(1-ρ))-1ρk/2}, as k→∞.Hence the concept of theQ-process will have most practical relevance to approximate the very late extinction case if ρ is near to zero and if |λ| is small compared with ρ. Note however that the very late extinction scenario is more likely to happen if ρ is near to unity because the time to extinction in most realizations will then be long (see Section 2.1).


## 2.5. Estimation of the Infection Parameter

We assume for this subsection that the parametersΨk of the epidemic model (1.1)–(2.1) are not entirely known. More precisely, we assume that the Ψk are of the form, for all k=1⋯d,
(2.28)Ψk(θ0)=akθ0+bk,
where ak>0 and bk⩾0 are constants, and θ0 is an unknown real parameter. We will write in what follows Ψ(θ0)=aθ0+b, where Ψ(θ0):=(Ψ1(θ0),…,Ψd(θ0)) and so forth. This general assumption corresponds in particular to the case where the Ψk are of the form (1.2).We estimateθ0 by the following WCLSE in model (1.1)–(2.1). This estimator generalizes the well-known Harris estimator [12] in a BGW process. Let Θ:=]θ1,θ2[, θ2>θ1>0, such that θ0∈Θ. The WCLSE is based on the normalized process Yn:=Xn/a·Xn-1 and is defined by
(2.29)θ^∣X0∣:=argminθ∈Θ∑k=1n‍(Yk-𝔼θ(Yk∣Xk-1))2=argminθ∈Θ∑k=1n‍(Xk-Ψ(θ)·Xk-1)2a·Xk-1.
We easily derive the following explicit form:
(2.30)θ^∣X0∣=∑k=1n(Xk-b·Xk-1)∑k=1na·Xk-1.
The normalization of the process Xn by a·Xn-1 appears to be the most natural and suitable for the following reasons. First, this normalization generalizes the normalization Xn/aXn-1 in the monotype case, which is the one leading to the Harris estimator m^ of m0=aθ0+b since we have, for d=1, aθ^X0+b=m^. It also corresponds, in the linear case b=0, to the maximum likelihood estimator of θ0. In addition, defining for any vector x, x_:=minixi and x¯:=maxixi, we have
(2.31)θ0+b_a¯⩽𝔼θ0((Yk-𝔼θ0(Yk∣Xk-1))2∣Xk-1)=θ0+b·Xk-1a·Xk-1⩽θ0+b¯a_,
hence the conditional variance of the error term Yk-𝔼θ0(Yk∣Xk-1) in the stochastic regression equation Yk=𝔼θ0(Yk∣Xk-1)+Yk-𝔼θ0(Yk∣Xk-1) is invariant under multiplication of the whole process, and bounded respectively to {Xn}n, leading to the quasi-optimality of θ^∣X0∣ at finite |X0| and n, in the sense of [13].Let us provide asymptotic results for the estimatorθ^∣X0∣ defined by (2.30), as the initial population size |X0|=X0+X-1+⋯+X-d+1 tends to infinity. We denote by mij(k)(θ) the (i,j)th entry in the kth power of the matrix M(θ) given by (2.3).Theorem 2.6.
Let us assume that, for eachi=1⋯d, there exists some αi∈[0,1] such that lim∣X0∣→∞X0,i|X0|-1=αi. Then θ^∣X0∣ is strongly consistent, that is, lim∣X0∣→∞θ^∣X0∣=a.s.θ0, and is asymptotically normally distributed:
(2.32)lim∣X0∣→∞∑k=1na·Xk-1σ2(θ^∣X0∣)(θ^∣X0∣-θ0)=𝒟𝒩(0,1),
where
(2.33)σ2(θ):=θ+∑k=1n∑j=1d∑i=1dαjbimji(k-1)(θ)∑k=1n∑j=1d∑i=1dαjbimji(k-1)(θ).Proof.
Let us first prove that, for eachk=1⋯n and each i=1⋯d,
(2.34)lim∣X0∣→∞Xk,i|X0|=a.s.∑j=1d‍αjmji(k)(θ0).
Using the branching property of the process {Xn}n ([6]) we write
(2.35)Xk,i=∑j=1X0,1‍Xk,i,j(1)+⋯+∑j=1X0,d‍Xk,i,j(d),
where, for all l=1⋯d and j=1⋯X0,l, Xk,i,j(l) is the ith coordinate of a d-type branching process at time k initialized by a single particle of type l. For k, i, and l fixed the random variables {Xk,i,j(l)}j are i.i.d. with mean value mli(k)(θ0). According to the strong law of large numbers and under the theorem assumption, we have, for every l=1⋯d such that X0,l≠0,
(2.36)lim|X0|→∞∑j=1X0,lXk,i,j(l)X0,l=a.s.mli(k)(θ0),
which together with the theorem assumption leads to (2.34).
To prove the consistency ofθ^∣X0∣ we apply (2.34) to (2.30), using the fact that Xk=Xk,1 and Xk-i=Xk-1,i, and obtain
(2.37)lim∣X0∣→∞θ^∣X0∣=a.s.∑k=1n∑j=1dαj(mj1(k)(θ0)-∑i=1dbimji(k-1)(θ0))∑k=1n∑i=1d∑j=1daiαjmji(k-1)(θ0).
By definition,
(2.38)mj1(k)(θ0)=∑i=1d‍mji(k-1)(θ0)mi1(θ0)=∑i=1d‍mji(k-1)(θ0)(aiθ0+bi),
hence (2.37) immediately leads to the strong consistency.
We are now interested in the asymptotic distribution ofθ^∣X0∣-θ0. We derive from (2.30) that
(2.39)∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=∑k=1n(Xk-Ψ(θ0)·Xk-1)∑k=1na·Xk-1.
By (1.1),
(2.40)Xk-Ψ(θ0)·Xk-1=∑i=1d‍∑j=1Xk-i‍(Yk-i,k,j-Ψi(θ0))=:∑i=1d‍∑j=1Xk-i‍Y°k-i,k,j,
where the {Yk-i,k,j}j are i.i.d. given ℱk-1, following a Poisson distribution with parameter Ψi(θ0), and the {Yk-i,k,j}i,j are independent given ℱk-1. Renumbering the Y°k-i,k,j we then obtain
(2.41)∑k=1n‍(Xk-Ψ(θ0)·Xk-1)=∑i=1d‍∑j=1∑k=1nXk-i‍Y°k-i,k,j.
Applying a central limit theorem for the sum of a random number of independent random variables (see, e.g., [14]), we obtain that, for all i=1⋯d,
(2.42)lim∣X0∣→∞∑j=1∑k=1nXk-iY°k-i,k,j∑k=1nXk-i=𝒟𝒩(0,aiθ0+bi).
We have used the fact that |X0| is a real positive sequence growing to infinity, and ∑k=1n‍Xk-i is a sequence of integer-valued random variables such that ∑k=1n‍Xk-i/|X0| converges in probability to a finite random variable.  In our case the limit is actually deterministic, since we have shown in (2.34) that
(2.43)lim|X0|→∞∑k=1nXk-i|X0|=a.s.∑k=1n‍∑j=1d‍αjmji(k-1)(θ0).
Using (2.41) in (2.39), we write
(2.44)∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=∑i=1d‍∑j=1∑k=1nXk-iζ°k-i,k,j∑k=1nXk-i∑k=1nXk-i∑k=1na·Xk-1.
Using again (2.34),
(2.45)lim∣X0∣→∞∑k=1na·Xk-i∑k=1na·Xk-1=a.s.∑k=1n∑j=1dαjmji(k-1)(θ0)∑k=1n∑j=1d∑l=1dαjalmjl(k-1)(θ0),
which, combined to (2.42) and (2.44), implies by Slutsky's theorem that
(2.46)lim∣X0∣→∞∑k=1n‍a·Xk-1(θ^∣X0∣-θ0)=𝒟𝒩(0,σ2(θ0)).
By (2.33) and the strong consistency, lim∣X0∣σ2(θ0)/σ2(θ^∣X0∣)=a.s.1, from which we finally deduce (2.32).


## 2.6. Estimation of the Infection Parameter in the Worst-Case Scenario

In order to make predictions of the evolution of the epidemic in case of a very late extinction, that is, in order to make predictions of the behavior of the conditioned process{Xn*}n introduced in Section 2.4, we need to estimate the parameter θ0 in the setting of this conditioned process. We point out that θ0 does not play the same role in the conditioned process {Xn*}n and in the unconditioned process {Xn}n, since, as shown in Proposition 2.2, this parameter interferes not only in the Poisson random variable but also in the Bernoulli one. It would thus be irrelevant to estimate θ0 with an estimator aimed for the unconditioned process, such as θ^∣X0∣. Let us notice that, according to (2.16), the process {Xn*}n could be written as a multitype branching process with state dependent immigration. Because of this state-dependence, and since the parameter θ0 acts in a nonlinear way in the immigration, the methods developed in estimation theory for branching processes with immigration (see, e.g., [15]) cannot be directly applied here. Similarly as in Section 2.5 we consider the WCLSE based on the process Yn*:=Xn*/a·Xn-1*, namely,
(2.47)θ^n*:=argminθ∈ΘSn(θ),Sn(θ):=∑k=1n‍(Yk*-f(θ,Xk-1*))2,
where Θ is defined in Section 2.5, and where
(2.48)f(θ0,Xk-1*):=𝔼θ0(Yk*∣Xk-1*)=Xk-1*·Ψ(θ0)+p(θ0,Xk-1*)a·Xk-1*.
Let ɛk*:=Yk*-f(θ0,Xk-1*) be the error term between the normalized process and its conditional expectation. We obtain that
(2.49)g(θ0,Xk-1*):=𝔼θ0((ɛk*)2∣Xk-1*)=Xk-1*·Ψ(θ0)+p(θ0,Xk-1*)(1-p(θ0,Xk-1*))a·Xk-1*,
which implies
(2.50)θ0+b_a¯⩽𝔼θ0((ɛk*)2∣Xk-1*)⩽θ0+b¯+1a_.
In what follows, we denote by f' the derivative of f with respect to θ, and similarly for the other quantities depending on θ.Theorem 2.7.
The estimatorθ^n* is strongly consistent, that is, limn→∞θ^n*=a.s.θ0, and has the following asymptotic distribution:
(2.51)limn→∞∑k=0nf′(θ^n*,Xk*)2∑k=0nf′(θ^n*,Xk*)2g(θ^n*,Xk*)(θ^n*-θ0)=𝒟𝒩(0,1).Remark 2.8.
Note that (2.51) involves the function f' and thus requires the knowledge of the derivative of the function uj given in (2.5), which is not an explicit function of θ since ρ is not either. However, ρ' satisfies ρ'=∑k=1d‍akρ-k[∑k=1d‍k(akθ+bk)ρ-k-1]-1; hence uj′ is known as soon as ρ can be computed. Consequently, denoting ρ by ρ(θ) when θ is the parameter of the model, (2.51) can be used as soon as ρ(θ^n*obs) is known; for this purpose one can for instance numerically approximate the largest solution ρ of ∑k=1d‍Ψk(θ^n*obs)ρ-k=1.Proof.
The proof heavily relies on a strong law of large numbers for homogeneous irreducible positive recurrent Markov chains applied to the conditioned process{Xn*}n and its stationary distribution πθ0 (given by Proposition 2.4), which states [16]: for every πθ0-integrable function h:ℕd∖{0}→ℝ,
(2.52)limn→∞1n∑k=0n-1‍h(Xk*)=a.s.∑j∈ℕd‍h(j)πθ0(j).
Note that the Perron's root ρ(θ) of M(θ) and the associated right normalized eigenvector u(θ) are 𝒞∞-functions of θ.
In order to prove the strong consistency ofθ^n* as n→∞, since f′′(·,Xk-1*) is not linear, we cannot use the general standard method based on the first-order expansion of Sn′(θ^n*) at θ0, and on the strong law of large numbers for martingales applied to Sn′(θ0) correctly normalized. We consequently use the following conditions given in [17]:(i)
f(.,Xk-1*) is Lipschitz on Θ; that is, there exists a nonnegative ℱk-1*-measurable function Ak (where ℱk-1*:=σ(X0*,…,Xk-1*)), satisfying, for all δ1, δ2∈Θ, |f(δ1,Xk-1*)-f(δ2,Xk-1*)|⩽Ak|δ1-δ2| a.s.,(ii)
lim¯k→∞𝔼θ0((ɛk*)2∣Xk-1*)<a.s.∞,(iii)
limn→∞infθ∈Θ∣θ-θ0∣⩾δ∑k=1n‍(f(θ0,Xk-1*)-f(θ,Xk-1*))2=a.s.∞. Let us note that, in the frame of a general model f, since Θ¯ is compact, (i) is satisfied as soon as f(θ,Xk-1*) has a first derivative in θ with supθ∈Θf'(θ,Xk-1*)<∞, (ii) is satisfied for any optimal estimator in the sense of [13] and moreover could be weakened (see [17]), and (iii) is a necessary condition. First, for all θ∈Θ and j∈ℕd, j≠0, f'(θ,j)=(a·j+p'(θ,j))(a·j)-1/2, where p' denotes the derivative of p with respect to θ, which thanks to (2.17) is bounded on Θ. Condition (i) is thus satisfied. Condition (ii) follows from (2.50). Condition (iii) comes from the fact that, for every δ>0 and every θ∈Θ such that ∣θ-θ0∣⩾δ, applying the mean value theorem to the 𝒞1-function p(·,Xk-1*),
(2.53)∑k=1n‍(f(θ0,Xk-1*)-f(θ,Xk-1*))2=(θ0-θ)2∑k=1n‍a·Xk-1*(1+p(θ0,Xk-1*)-p(θ,Xk-1*)(θ0-θ)a·Xk-1*)2⩾δ2∑k=1n‍a·Xk-1*infθ∈Θ¯(1+p'(θ,Xk-1*)a·Xk-1*)2.
Let us show that the function j↦a·jinfθ∈Θ¯(1+p'(θ,j)/a·j)2 is πθ0-integrable. For every θ∈Θ¯, j∈ℕd and j≠0, denoting Ξ:=supi,θ∈Θ¯{ui(θ),∣ui′(θ)∣}<∞ (by continuity of ui′(θ) on Θγ⊃Θ¯), and u_:=mini,θ∈Θ¯ui(θ)>0,
(2.54)|p'(θ,j)|=|(ui′(θ)j·Ψ(θ)+u1(θ)j·a)∑i=2djiui(θ)-u1(θ)j·Ψ(θ)∑i=2djiui′(θ)(u1(θ)j·Ψ(θ)+∑i=2djiui(θ))2|⩽Ξ2u_23j·Ψ(θ)∑i=2dji(j·Ψ(θ)+∑i=2dji)2⩽3Ξ24u_2=:C1.
Hence, for all j≠0, ∣a·jinfθ∈Θ¯(1+p'(θ,j)/a·j)2∣⩽(1+C1/a_)2a·j, and applying (2.52) together with (2.22) we obtain that
(2.55)limn→∞1n+1∑k=1n‍a·Xk-1*infθ∈Θ¯(1+p′(θ,Xk-1*)a·Xk-1*)2=a.s.∑j∈ℕd‍a·jinfθ∈Θ¯(1+p′(θ,j)a·j)2πθ0(j).
Let j≠0 fixed. Since, for all θ∈Θ¯, p'(θ,j)≠-a·j, the extreme value theorem implies that infθ∈Θ¯(1+p'(θ,j)/a·j)2>0. Hence the right term in (2.55) is strictly positive, which together with (2.55) leads to (iii).
Let us now consider the asymptotic distribution ofθ^n*-θ0. For this purpose, we follow the steps of the proof of Proposition 6.1 in [17]. Writing the Taylor expansion of Sn′(θ^n*) at θ0 we obtain that θ^n*-θ0=-Sn′(θ0)/Sn′′(θ~n), for some θ~n=θ0+tn(θ^n*-θ0), with tn∈]0,1[. Since Sn′(θ0)=-2∑k=1n‍ɛk*f'(θ0,Xk-1*), we can write
(2.56)n(θ^n*-θ0)=∑k=1nɛk*f′(θ0,Xk-1*)n(Fnn)-1(12Sn′′(θ~n)Fn)-1,
where Fn:=∑k=1n‍(f'(θ0,Xk-1*))2. Let us first show that
(2.57)limn→∞Fnn=a.s.∑j∈ℕd‍(f′(θ0,j))2πθ0(j).
This is an application of (2.52) and (2.22), since, for all j∈ℕd, j≠0,
(2.58)(f′(θ0,j))2⩽a·j+2C1+C12a_.
In view of (2.56), we now prove that
(2.59)limn→∞Sn′′(θ~n)Fn=a.s.2.
Computing Sn′′ thanks to the formula Sn(θ)=∑k=1n‍(ɛk*+f(θ0,Xk-1*)-f(θ,Xk-1*))2, it appears that (2.59) is true, as soon as the following holds:
(2.60)limn→∞supθ∈Θ|∑k=1nɛk*f′′(θ,Xk-1*)|Fn=a.s.0,(2.61)limn→∞∑k=1n(f′(θ~n,Xk-1*))2Fn=a.s.1,(2.62)limn→∞∑k=1n(f(θ0,Xk-1*)-f(θ~n,Xk-1*))f′′(θ~n,Xk-1*)Fn=a.s.0.
Let us prove (2.60)–(2.62). First, (2.60) is given by a strong law of large numbers proved in [17, Proposition 5.1]. The latter can be indeed applied since limnFn=a.s.∞ (as an immediate consequence of the stronger result (2.57)), and since f′′(·,Xk-1*) fulfills the required Lipschitz condition. Indeed we proved earlier that ui′′′(θ) is continuous on the compact set Θ¯ and is thus bounded on Θ¯, which implies that f′′′(·,Xk-1*)=p′′′(·,Xk-1*)(a·Xk-1*)-1/2 is bounded by a ℱk-1*-measurable function. In view of (2.61), we consider the function f(θ,j)2 and its derivative 2f'(θ,j)f′′(θ,j). Similarly as for (2.54), one can show that there exists a constant C2>0 such that for all θ∈Θ¯, and all j≠0, ∣p′′(θ,j)∣⩽C2. This implies
(2.63)|2f′(θ,j)f′′(θ,j)|⩽2|(a·j+p′(θ,j))p′′(θ,j)a·j|⩽2C2(1+C1a_).
Consequently,
(2.64)|∑k=1n(f′(θ~n,Xk-1*))2-(f′(θ0,Xk-1*))2|Fn⩽2C2(1+C1a_)|θ^n*-θ0|(Fnn)-1,
which by (2.57) and the strong consistency of θ^n* almost surely tends to 0. Writing
(2.65)∑k=1nf(θ~n,Xk-1*)2Fn=1+∑k=1n((f′(θ~n,Xk-1*))2-(f′(θ0,Xk-1*))2)Fn,
this implies (2.61). It now remains to prove (2.62). We write
(2.66)|∑k=1n(f(θ0,Xk-1*)-f(θ~n,Xk-1*))f′′(θ~n,Xk-1*)|Fn⩽1Fn∑k=1n‍a·Xk-1*|θ0-θ~n|+|p(θ0,Xk-1*)-p(θ~n,Xk-1*)|a·Xk-1*|p′′(θ~n,Xk-1*)|⩽1Fn∑k=1n‍(|θ0-θ~n|+C1|θ0-θ~n|a_)C2⩽|θ0-θ^n*|C2(1+C1a_)(Fnn)-1,
which thanks to (2.57) and the strong consistency of θ^n* implies (2.62). In view of (2.56), we finally want to prove that ∑k=1n‍ɛk*f'(θ0,Xk-1*)/n converges in distribution. For this purpose we make use of a central limit theorem for martingale difference arrays [18, 19]; if {Mk(n),ℱk(n),1⩽k⩽n}, n⩾1 is a sequence of square integrable martingales with associated Meyer process 〈M〉(n)=(〈M〉k(n))1⩽k⩽n satisfying limn→∞〈Mn〉(n)=Pc2 for some constant c, and such that, for all ɛ>0,
(2.67)limn→∞∑k=1n‍𝔼[|Mk(n)-Mk-1(n)|21{|Mk(n)-Mk-1(n)|⩾ɛ}∣ℱk-1(n)]=P0,
then limn→∞Mn(n)=𝒟𝒩(0,c2). Let us define Mk(n):=∑l=1k‍ɛl*f'(θ0,Xl-1*)/n, for every k⩽n. First, for any k⩽n, 𝔼θ0(ɛk*f'(θ0,Xk-1*)/n∣Xk-1*)=0. Second,
(2.68)𝔼θ0((ɛk*f′(θ0,Xk-1*)n)2∣Xk-1*)=(f′(θ0,Xk-1*))2g(θ0,Xk-1*)n,
hence {Mk(n)}k⩽n is a sequence of square integrable martingales. Moreover, using inequalities (2.50) and (2.58), we obtain, by (2.22),
(2.69)∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j)⩽(θ0+b¯+1a_)(∑j∈ℕd‍a·jπθ0(j)+2C1+C12a_)<∞.
So, by means of (2.52),
(2.70)limn→∞〈Mn〉(n)=limn→∞∑k=1n‍𝔼θ0((ɛkf′(θ0,Xk-1*)n)2∣Xk-1*)=a.s.∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j).
Third, using Cauchy-Schwarz and Bienaymé-Chebyshev inequalities,
(2.71)∑k=1n‍𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|21{|ɛk*f′(θ0,Xk-1*)/n|⩾ɛ}∣Xk-1*]⩽∑k=1n‍(𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|4∣Xk-1*])1/2(ℙθ0[|ɛk*f′(θ0,Xk-1*)n|⩾ɛ∣Xk-1*])1/2⩽1n3/2ɛ∑k=1n‍|f′(θ0,Xk-1*)|3(𝔼θ0[(ɛk*)4∣Xk-1*])1/2(𝔼θ0[(ɛk*)2∣Xk-1*])1/2.
Let us compute 𝔼θ0((ɛk*)4∣Xk-1*). We can show that the 4th central moment of the independent sum of a Poisson and a Bernoulli random variables equals μ4+6μ2γ2+γ4, where μi and γi denote the ith central moment of the Poisson and of the Bernoulli variable. If these variables have parameter λ and p, respectively, then μ4=λ(1+3λ), μ2=λ, γ4=p(1-p)(3p2-3p+1)∈[0,1], and γ2=p(1-p)∈[0,1]. We thus obtain
(2.72)|𝔼θ0((ɛk*)4∣Xk-1*)|⩽Ψ(θ0)·Xk-1*(7+3Ψ(θ0)·Xk-1*)+1(a·Xk-1*)2.
Hence
(2.73)|f*′(θ0,Xk-1*)|3(𝔼θ0[ɛk4∣Xk-1*])1/2(𝔼θ0[ɛk2∣Xk-1*])1/2⩽(a·Xk-1*+C1a_)3Ψ(θ0)·Xk-1*(7+3Ψ(θ0)·Xk-1*)+1a·Xk-1*(θ0+b¯+1a_)1/2.
Since the highest power of Xn* involved in (2.73) is 3/2, and since by Proposition 2.4 the stationary distribution πθ0 has finite second moments, we can apply (2.52) to (2.71) and obtain that
(2.74)limn→∞∑k=1n‍𝔼θ0[|ɛk*f′(θ0,Xk-1*)n|21{|ɛk*f′(θ0,Xk-1*)/n|⩾ɛ}∣Xk-1*]=a.s.0.
It then ensues from the central limit theorem mentioned above that
(2.75)limn→∞∑k=1nɛk*f′(θ0,Xk-1*)n=𝒟𝒩(0,∑j∈ℕd‍(f′(θ0,j))2g(θ0,j)πθ0(j)).
Finally, (2.56), together with (2.57), (2.59), (2.75), and Slutsky's theorem, implies that
(2.76)limn→∞n(θ^n*-θ0)=𝒟𝒩(0,∑j∈ℕd(f′(θ0,j))2g(θ0,j)πθ0(j)(∑j∈ℕd(f′(θ0,j))2πθ0(j))2),
which leads to
(2.77)limn→∞n(∑j∈ℕdf′(θ0,j)2πθ0(j))2∑j∈ℕdf′(θ0,j)2g(θ0,j)πθ0(j)(θ^n*-θ0)=𝒟𝒩(0,1).
It now remains to prove the asymptotic distribution (2.51). Thanks to what precedes, this result is immediate as soon as we prove that
(2.78)limn→∞1n+1∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)=a.s.∑j∈ℕd‍(f'(θ0,j))2g(θ0,j)πθ0(j),
as well as the equivalent result for the numerator. For this purpose, we write
(2.79)∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)=∑k=0n‍(f′(θ0,Xk*))2g(θ0,Xk*)+∑k=0n‍[(f′(θ^n*,Xk*))2g(θ^n*,Xk*)-(f′(θ0,Xk*))2g(θ0,Xk*)],
and show that (f'(·,j))2g(·,j) has a bounded derivative and is thus Lipschitz. We have indeed ∣g'(θ,j)∣=∣a·j+p'(θ,j)∣(a·j)-1⩽1+C1a_-1, hence
(2.80)|2f′′(θ,j)f′(θ,j)g(θ,j)+(f′(θ,j))2g′(θ,j)|⩽2C2(1+C1a_-1)a_-1(Ψ(θ)·j+1)+(1+C1a_-1)2(a·j+C1)⩽C3Ψ(θ2)·j,
for some constant C3>0. This enables us to write
(2.81)1n+1∑k=0n‍|(f′(θ^n*,Xk*))2g(θ^n*,Xk*)-(f′(θ0,Xk*))2g(θ0,Xk*)|⩽|θ^n*-θ0|C31n+1∑k=0n‍Ψ(θ2)·Xk*.
By the strong consistency of θ^n* together with (2.52) and (2.22), (2.81) almost surely tends to zero. Combined with (2.52) and (2.69) in (2.79), this implies (2.78).


## 3. Application to Real Epidemics

### 3.1. Explicit Epidemic Model

Although the epidemic model (1.1) has a clear interpretation in terms of the disease propagation, it can be also proved [9] that it is the limit process, as the initial population size tends to infinity, of the incidence of infectives described by a thorough process in a large branching population, taking into account all the health states of the disease and describing in detail the infection and latent processes for each individual in the population. This detailed process is a multitype branching process with age and population size dependence, taking into account the variability of many individual factors such as the reproduction, the survival, the transmission of the disease, and the latent time. We assume a disease of the general SEIR type. We also assume a random latent period and, for the sake of simplicity, a duration of the state I of one time unit only. Hence the incidence of infectives exactly corresponds to the incidence of cases.The main assumptions for the construction of the limit model{Xn}n are, concerning the disease as follows: (i) at the initial time the disease is rare and the total population size is large; (ii) the infection via horizontal route is of Reed-Frost type, with the probability for a susceptible individual to become infected by a given dose of pathogens being inversely proportional to the total population size; (iii) the individual survival law is the same for E and S individuals and is independent of the population and of the time; (iv) the latent time law given the individual survival is independent of the time, the individual age, and the infectives population during the latent period. This last property is possible if we assume that overinfection during the incubation has a negligible effect on the latent time. We moreover assume that the whole healthy population size is relatively stable over time.The limit process (1.1) is then obtained in an inductive way as the limit in distribution of {In}n as the initial total population size |N0|→∞, where In denotes the number of new infectives at time n. More precisely, denoting by En the incidence of exposed individuals at time n and by am the largest individual survival age, we obtain that {Xn,Zn}n=𝒟limN0→∞{In,En}n, with
(3.1)Xn∣(Xn-1,…,Xn-am+1)=𝒟Poiss(∑k=1am-1‍ΨkXn-k),(3.2)Zn∣Xn=𝒟Poiss(Ψ0Xn),
and where for all k=1⋯am-1, Ψk satisfies (1.2) in Section 1, with the latent period distribution given survival, Pinc,a, assumed independent of a, being denoted by Pinc. Thus
(3.3)Ψk=θ0Pinc(k)∑a=k+1am‍Page(a)+pmatPinc(k)Page(k+1),
and similarly,
(3.4)Ψ0=θ0+pmatPage(1).
Moreover, when assuming a relatively stable population size and denoting by S(a) the individual survival probability until age a at least, then Page(a) is equal to S(a)(∑a'=1am‍S(a'))-1.


### 3.2. Illustration: The BSE Epidemic in Great Britain

In this section we provide an illustration of the methodology developed in Section2 by studying the decay phase of the BSE epidemic in Great Britain, based on the observations of the yearly diagnosed cases from 1989 until 2011, with 1988 being the date of the main control measure. The disease that was first officially identified in 1986 [20] reached its peak in 1992 (36682 cases) and is obviously now in its decay phase. Since only a very few cases were recently reported [5], namely, 11 in 2010 and 5 in 2011, the spread of the disease should a priori “soon” come to an end. It has been accepted that the epidemic is fading out, with a very low level of risk for cattle and humans. However it might be interesting to have a more precise idea of the extinction phase of the disease, for instance, of its length and of its intensity.


#### 3.2.1. Choice of the Epidemic Model

BSE is a transmissible disease through ingestion of prions (horizontal route) and through maternal route. It is commonly accepted (see, e.g., [21]) that the infectious phase (including the clinical state) lasts at most one year, while the latent state can be of several years. If one compares the epidemic data with the total population of around 9 million cattle in Great Britain, it is reasonable to assume that, even at the peak time, BSE may be considered as a rare disease in a large population. We moreover assume that the probability for an animal to become infected by a horizontal route follows a Reed-Frost type model as in (ii) of Section 3.1, and that the probability for a calf to become infected by its dam is nonnull and constant over time. Moreover, for the sake of simplicity and parsimony, we do not take into account potential heterogeneity factors such as the different regulations from 1989 (the main regulation was the feed ban of July 1988 and was shown to be quite efficient [22]), the different types of breeding and of races, the age for animals older than one year, and the evolution of the surveillance system and of diagnosis tests. So the process {Xn}n defined by (3.1)–(3.4) appears relevant to model the spread of this fatal disease after the 1988 ban, where the infectious state I corresponds to the end of the incubation period and the clinical state, and R corresponds to the death (assumed to be either by routine slaughtering for E and S cattle, and control slaughtering for I cattle). Choosing a time step of one year, Xn then represents the yearly incidence of notified cases, which corresponds to the available data provided by the World Organisation for Animal Health [5].


#### 3.2.2. Choice of the Model Parameters

In order to predict the future epidemic evolution, we need to evaluate all the parameters involved in (3.3). These parameters were already studied in previous works [22–25] using different models, amounts of observations, and kinds of estimators. Since it might have a crucial role in the decay phase of the disease, we focus here on θ0 the infection parameter from 1989 via a horizontal route of transmission (i.e., the mean number per infective and per year of newly infected). Until the 1988 feed ban regulation, the main routes of transmission of BSE were horizontal via protein supplements (Meat and Bone Meal, milk replacers), and maternal from a dam to its calf. Since a previous statistical study [22] concluded to the full efficiency of the 1988 ban and since most of cattle are slaughtered before the age of 10 years, the fact that cases of BSE are still observed more than 20 years later could suggest the existence of a remaining source of infection, either via a maternal transmission route, or via a horizontal one, for example, via the ingestion of excreted prions from alive infected animals, or from prions left in the environment. The prediction of the future disease spread thus strongly relies on the intensity of this infection, quantified by the parameter θ0. We provide an estimation of θ0 together with a confidence interval, which will be used for prediction estimations. We provide in addition a sensitivity analysis to the other parameters of the model, which are chosen as follows. First, the maximal age of the individuals in the population am is set to am=10, which corresponds to the largest reported survival age with a nonnegligible probability. Next, we set pmat=0.1, which is the largest order of magnitude commonly accepted for the maternal infection probability [23, 25]. Following a previous work, [22], we assume here a discretized Weibull distribution with parameters α,β for the distribution of the latent period; that is to say, for each k⩾1, Pinc(k):=e-((α-1)/αβα)(k-1)α-e-((α-1)/αβα)kα, where α is a shape parameter, and β is the mode of the probability density of the corresponding continuous Weibull distribution. The parameters α and β have a very bad identifiability with the infection parameters on a sole monotonous phase and are consequently estimated in [22] on the whole epidemic series (growth and decay), by Bayesian maximum a posteriori (MAP) estimations. We set α=α^MAP=3.84 and β=β^MAP=7.46. Since the death of susceptible and exposed animals is mostly due to routine slaughtering, the survival probability S(k) (probability for an apparently healthy animal to survive at least until k years) corresponds to the probability of being slaughtered after the age of k years. It is derived from existing literature [24].


#### 3.2.3. Estimation of the Infection Parameter

Our study is based on the yearly number of cases of BSE reported in Great Britain from 1989 until 2011 (Figure1) provided by the World Organisation for Animal Health [5], denoted by Xnobs for each year n. We set d=am-1=9 and denote by Xnobs the d-dimensional vector (Xnobs,…,Xn-d+1obs). We choose the first time n=1 of the epidemic model such that the model with its initial values covers the period starting from 1989, which we consider here as a time-homogeneous period. Hence X-d+1=X1989obs, and n=1 corresponds to the year 1989+d=1998. We estimate θ0 by the WCLSE (2.30) in model (1.1)–(3.3), where {pmat,α,β,{S(k)}k} are given by the previous values, and where X0=X1997obs. According to [5], |X0|=|X1997obs|=167977. We are thus close to the asymptotic |X0|→∞. The number of observations is n=14. The estimator (2.30) provides the estimation θ^∣X0∣obs=2.4324. We point out that this estimation is of the same order of magnitude as the maximum a posteriori Bayesian estimation θ^MAP=2.43 based on the whole epidemic until 2007, assuming a uniform prior probability [22]. Using (2.32) we obtain the following confidence interval [θ^min,θ^max] with asymptotic probability 95%, where θ^min:=θ^∣X0∣-1.96c^1-1, θ^max:=θ^∣X0∣+1.96c^1-1, and c^1:=(∑k=1n‍a·Xk-1/σ2(θ^∣X0∣))1/2. Assuming αi=X1997-i+1obs|X1997obs|-1, we get c^1obs=40.7343 (observed value of c^1), and
(3.5)ℙ(θ0∈[θ^min,θ^max])≃95%,θ^minobs=2.3842,θ^maxobs=2.4805.
Although this confidence interval is an asymptotic one, as |X0|→∞, it is a very good approximation of the true confidence interval for a finite |X0|, since |X0| is here very large. Since the estimation of θ0 relies on the values given to the other model parameters {pmat,α,β}, we evaluate in addition its sensitivity to the values of these parameters. For this purpose, we compute the estimation of θ0 and the associated confidence interval with asymptotic probability 95%, for different values of (pmat,α,β) (Table 1). The first line of the table corresponds to the parameters chosen for the model. In each of the four following lines, we fix two coordinates and choose an extremal (unrealistic) value for the third one. It appears that the estimation of θ0 is almost independent of the value of the maternal infection parameter. However, the estimation seems more strongly dependent on the parameters of the latent period distribution. Nevertheless, even for very unrealistic values (α,β), all the estimations of θ0 remain in the same order of magnitude of several units. This is really small compared to estimations obtained for the infection via Meat and Bone Meal or lactoreplacers (before 1989) which are of the order of 1000 [22]. However, although these estimations are all very small, θ0 seems nonnull. This could suggest the existence of a minor but nonnull infection source which is not of maternal type.Table 1
Sensitivity analysis. Estimation of the infection parameter θ0, and its confidence interval [θ^minobs,θ^maxobs] with asymptotic probability 95%, for different values of the maternal infection parameter pmat and of the latency parameters (α,β). The values (3.84,7.46) correspond to the Bayesian MAP estimations (α^MAP,β^MAP), and 0.1 to the largest commonly accepted order of magnitude for pmat. The estimations of θ0 are based on the observed data over the years 1989–2011.
p
mat
α
β
θ
^
|
X
0
|
obs
[
θ
^
min
obs
,
θ
^
max
obs
]
0.1
3.84
7.46
2.4324
[
2.3842,2.4805
]
0
3.84
7.46
2.4860
[
2.4379,2.5342
]
1
3.84
7.46
1.9492
[
1.9014,1.9970
]
0.1
2
7.46
2.7835
[
2.7287,2.8382
]
0.1
20
7.46
4.0186
[
3.9395,4.0977
]
0.1
3.84
1
1.0127
[
0.9925,1.0329
]
0.1
3.84
10
6.2128
[
6.0914,6.3341
]
0.1
3
6
1.5402
[
1.5095,1.5710
]
0.1
4
5
1.0227
[
1.0020,1.0434
]Figure 1
Yearly number of cases of BSE reported in Great Britain from 1987 to 2011 [5]. Our study only takes into account the data from 1989, after the 1988 feed ban law.


#### 3.2.4. Extinction of the Epidemic

We know thanks to Proposition2.1 that {Xn}n becomes extinct almost surely if and only if R0=∑k=1d‍Ψk(θ0)⩽1. The estimated basic reproduction is here R0(θ^∣X0∣obs)=0.1072. Moreover, solving with a computing program the equation ∑k=1d‍Ψk(θ^∣X0∣obs)ρ-k=1, we obtain the following value for the Perron's root ρ(θ^∣X0∣obs)=0.6665, which provides the speed of decay of the expected yearly incidence of cases (see (2.6)); from a certain time, the expected number of new cases will decrease from around 33% every year.


#### 3.2.5. Prediction of the Incidences of Cases and Incidences of Infected Cattle

Let us predict the spread of the disease from 2012 by means of simulations of{Xn}n, where θ0 is replaced by its previous estimation θ^∣X0∣obs=2.4324, and where the initial time of the model is 2011, that is, X0=X2011obs. The simulations are done recursively using the transition law (3.1). We point out that the model initialized by X0=X1997obs provides quite realistic simulations on the period 1998–2011 compared to the real observations on the same period, as illustrated in Figure 2(a). The epidemic process {Xn}n thus seems to provide a satisfying prediction of the overall evolution of the real epidemic. In order to predict the incidences of futures cases, we simulate 1000 trajectories of {Xn}n initialized by the observed values X2011obs, with the estimated infection parameter θ^∣X0∣obs=2.4324. We illustrate in Figure 2(c), for each year from 2012, the maximum, minimum, median, 2.5%, and 97.5% quantiles associated with these 1000 realizations. It is also relevant to study and predict the evolution of the incidence of infected cattle in the population, which represents the hidden face of the epidemic. The incidence Zn of infected cattle at time n, conditionally on the number Xn of cases at that time, is given by the Poisson distribution (3.2). For every n⩾2012 and for each of the 1000 previously simulated values Xn, we generate one realization of Zn. We then illustrate in Figure 2(d), the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations.Figure2(a) 10 simulations of {Xn}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 2(b) 5 simulations of {Xn}n initialized by X2011obs. Figures 2(c) and 2(d) prediction, based on 1000 simulations of the process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣obs=2.4324.
(a)
(b)
(c)
(d)


#### 3.2.6. Prediction of the Year of Extinction

LetT:=2011+inf{n⩾1,Xn=0} denote the extinction year of the epidemic process {Xn}n. According to (2.7) and denoting by fθ0 the offspring generating function of {Xn}n defined in (2.2) (from now on we let the dependence in θ0 appear in the notation), we have ℙθ0(T⩽2011+n)=fθ0,n(0)X2011obs, for every n⩾1, which by iterating fθ0 can be computed explicitly. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the extinction time (see (2.8)): n0.5T=2028, n0.95T=2035 and n0.99T=2039. Keeping in mind that T corresponds to the complete extinction of the epidemic (i.e., d=9 consecutive years without any case), these results actually mean that for the infection parameter θ0=2.4324, with probability larger than 50% (resp., 95% and 99%), no case will arise in the population from year 2020 (resp., 2027 and 2031). Moreover, in order to take into account the uncertainty around the estimation θ^∣X0∣obs of the infection parameter θ0, we make use of the asymptotic confidence interval (3.5) of θ0 and of the fact that θ↦ℙθ(T⩽n) is a decreasing function of θ, which implies that for every n⩾2011,
(3.6)ℙ(ℙθ0(T⩽n)∈[ℙθ^max(T⩽n),ℙθ^min(T⩽n)])≃95%.
We collect in Table 2 the observed interval [ℙθ^maxobs(T⩽n),ℙθ^minobs(T⩽n)], for each n⩾2020 (if n<2019 we have, conditionally on the initial value X2011obs, ℙ(Xn=0)=0 because of the memory which is not equal to 0). Note that these intervals are very narrow, leading to an accurate estimation of ℙθ0(T⩽n).Table 2
Cumulative distribution function of the year of extinction computed with the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npT for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(T⩽npT) based on (3.6).
n
ℙ
θ
^
max
obs
(
T
⩽
n
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
2020
0.0000
0.0000
2030
0.7610
0.7807
2040
0.9950
0.9957
2021
0.0000
0.0000
2031
0.8310
0.8465
2041
0.9967
0.9972
2022
0.0010
0.0014
2032
0.8816
0.8934
2042
0.9978
0.9981
2023
0.0121
0.0152
2033
0.9186
0.9273
2043
0.9985
0.9988
2024
0.0496
0.0592
2034
0.9451
0.9513
2044
0.9990
0.9992
2025
0.1211
0.1390
2035
0.9633
0.9677
2045
0.9995
0.9994
2026
0.2325
0.2579
2036
0.9755
0.9786
2046
0.9996
0.9996
2027
0.3756
0.4047
2037
0.9835
0.9857
2047
0.9997
0.9998
2028
0.5303
0.5584
2038
0.9889
0.9904
2048
0.9998
0.9998
2029
0.6619
0.6862
2039
0.9925
0.9936
2049
0.9999
0.9999


#### 3.2.7. Prediction of the Epidemic Size

LetN:=∑n=1T-2011‍Xn be the total size of the future epidemic from 2012 (total number of cases from 2012 until the extinction of the epidemic). We compute the distribution of N using (2.11), conditionally on the event {X0=X2011obs}. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the total epidemic size (see (2.12)): n0.5N=22, n0.95N=31 and n0.99N=35. From (2.13) we deduce that the mean value and variance of N for the parameter θ^minobs (resp., θ^maxobs) are 21 and 27 (resp., 22 and 28). Moreover, we obtain similarly as for the extinction time and thanks to (3.5) a confidence interval [ℙθ^max(N⩽n),ℙθ^min(N⩽n)] of ℙθ0(N⩽n), for every n∈ℕ, with asymptotic probability 95%, collected in Table 3.Table 3
Cumulative distribution function of the total size of the epidemic for the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npN for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(N⩽npN).
n
ℙ
θ
^
max
obs
(
N
⩽
n
)
ℙ
θ
^
min
obs
(
N
⩽
n
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
⋯
⋯
⋯
16
0.1293
0.1669
28
0.8702
0.9048
5
0.0001
0.0000
17
0.1795
0.2263
29
0.8942
0.9302
6
0.0002
0.0003
18
0.2373
0.2915
30
0.9283
0.9495
7
0.0005
0.0008
19
0.3025
0.3637
31
0.9480
0.9646
8
0.0013
0.0020
20
0.3741
0.4409
32
0.9635
0.9759
9
0.0027
0.0046
21
0.4486
0.5189
33
0.9742
0.9838
10
0.0060
0.0094
22
0.5231
0.5934
34
0.9824
0.9892
11
0.0119
0.0178
23
0.5975
0.6639
35
0.9881
0.9928
12
0.0214
0.0310
24
0.6655
0.7287
36
0.9917
0.9952
13
0.0365
0.0518
25
0.7282
0.7840
37
0.9945
0.9968
14
0.0585
0.0801
26
0.7829
0.8323
38
0.9964
0.9980
15
0.0891
0.1181
27
0.8304
0.8726
39
0.9977
0.9989


#### 3.2.8. Study of the Very Late Extinction Case

In order to predict the behavior of the “most dangerous” evolution of the epidemic, we first need to compute the estimationθ^n*obs (2.47) based on the data in Great Britain. The number of available observations is only n=14, so we are far from the asymptotic setting n→∞ of Theorem 2.7. However, the large value of |X0| can make us hope for a good accuracy. We point out that, by making use of the estimator θ^n* on the real data, we make an unverifiable assumption on the future of the epidemic; we consider the data as if they were the beginning of a trajectory with very late extinction. This should have the following consequence: the estimation provided by θ^n*obs should a priori be a bit smaller than the value 2.4324 provided by θ^∣X0∣obs. Indeed we obtain θ^n*obs=2.4305. Using this value, we deduce from (2.51) the confidence interval [θ^min*,θ^max*] of θ0 with asymptotic probability 95%, where θ^min*:=θ^n*-1.96c^2-1, θ^max*:=θ^n*+1.96c^2-1, and
(3.7)c^2:=∑k=0n‍(f′(θ^n*,Xk*))2[∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)]-1/2.
Therefore, c^2obs=40.6988, and
(3.8)ℙ(θ0∈[θ^min*,θ^max*])≃95%,θ^min*obs=2.3823,θ^max*obs=2.4787,
which is of the same magnitude order as the confidence interval [2.3842,2.4805] obtained with the unconditioned process (see (3.5)). Let us predict the most dangerous evolution thanks to the transition law of the conditioned process given by Proposition 2.2. Note that if one computes the eigenvalues ρ and λ introduced in Proposition 2.5, one obtains ρ(θ^*obs)=0.6664 and |λ(θ^*obs)|=0.5570. It thus appears that the convergence of the epidemic process conditioned on nonextinction at time k to the Q-process as k→∞ is not very fast. As a consequence the study of the Q-process only provides information about the behavior of the disease spread in the case of an extremely late extinction. First, we see thanks to Figure 3(a) that the simulations provided by the conditioned process initialized by X0=X1997obs, and where θ0 is estimated by θ^n*obs=2.4305, are true to the real observations on the period 1999–2011. Figure 3(b) is an example of one simulation in the period 2012–2040 of the conditioned process, for X0*=X2011obs. It appears that the values of this simulated trajectory are rapidly very small and of course are never equal to 0 for d=9 consecutive times. For a finer prediction, we simulate 1000 realizations of this process from 2012 until 2050, with X0*=X2011obs. Moreover, for every n⩾2012 and for each of the 1000 simulated values Xn*, we make one realization of the incidence Zn of infected cattle at time n, according to the law given by (3.2). Figures 3(c) and 3(d) represent the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations of, respectively, the incidence of cases and infected cattle, in case of an extreme late extinction. It appears thanks to Figures 3(c) and 3(d) that the supposedly “most dangerous” trajectories nevertheless do not reach high values and do not present a new peak of epidemic.Figure3(a) 10 simulations of {Xn*}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 3(b) one simulation of {Xn*}n initialized by X2011obs. Figures 3(c) and 3(d) prediction, based on 1000 simulations of the conditioned process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣*obs=2.4305.
(a)
(b)
(c)
(d)


## 3.1. Explicit Epidemic Model

Although the epidemic model (1.1) has a clear interpretation in terms of the disease propagation, it can be also proved [9] that it is the limit process, as the initial population size tends to infinity, of the incidence of infectives described by a thorough process in a large branching population, taking into account all the health states of the disease and describing in detail the infection and latent processes for each individual in the population. This detailed process is a multitype branching process with age and population size dependence, taking into account the variability of many individual factors such as the reproduction, the survival, the transmission of the disease, and the latent time. We assume a disease of the general SEIR type. We also assume a random latent period and, for the sake of simplicity, a duration of the state I of one time unit only. Hence the incidence of infectives exactly corresponds to the incidence of cases.The main assumptions for the construction of the limit model{Xn}n are, concerning the disease as follows: (i) at the initial time the disease is rare and the total population size is large; (ii) the infection via horizontal route is of Reed-Frost type, with the probability for a susceptible individual to become infected by a given dose of pathogens being inversely proportional to the total population size; (iii) the individual survival law is the same for E and S individuals and is independent of the population and of the time; (iv) the latent time law given the individual survival is independent of the time, the individual age, and the infectives population during the latent period. This last property is possible if we assume that overinfection during the incubation has a negligible effect on the latent time. We moreover assume that the whole healthy population size is relatively stable over time.The limit process (1.1) is then obtained in an inductive way as the limit in distribution of {In}n as the initial total population size |N0|→∞, where In denotes the number of new infectives at time n. More precisely, denoting by En the incidence of exposed individuals at time n and by am the largest individual survival age, we obtain that {Xn,Zn}n=𝒟limN0→∞{In,En}n, with
(3.1)Xn∣(Xn-1,…,Xn-am+1)=𝒟Poiss(∑k=1am-1‍ΨkXn-k),(3.2)Zn∣Xn=𝒟Poiss(Ψ0Xn),
and where for all k=1⋯am-1, Ψk satisfies (1.2) in Section 1, with the latent period distribution given survival, Pinc,a, assumed independent of a, being denoted by Pinc. Thus
(3.3)Ψk=θ0Pinc(k)∑a=k+1am‍Page(a)+pmatPinc(k)Page(k+1),
and similarly,
(3.4)Ψ0=θ0+pmatPage(1).
Moreover, when assuming a relatively stable population size and denoting by S(a) the individual survival probability until age a at least, then Page(a) is equal to S(a)(∑a'=1am‍S(a'))-1.


## 3.2. Illustration: The BSE Epidemic in Great Britain

In this section we provide an illustration of the methodology developed in Section2 by studying the decay phase of the BSE epidemic in Great Britain, based on the observations of the yearly diagnosed cases from 1989 until 2011, with 1988 being the date of the main control measure. The disease that was first officially identified in 1986 [20] reached its peak in 1992 (36682 cases) and is obviously now in its decay phase. Since only a very few cases were recently reported [5], namely, 11 in 2010 and 5 in 2011, the spread of the disease should a priori “soon” come to an end. It has been accepted that the epidemic is fading out, with a very low level of risk for cattle and humans. However it might be interesting to have a more precise idea of the extinction phase of the disease, for instance, of its length and of its intensity.


### 3.2.1. Choice of the Epidemic Model

BSE is a transmissible disease through ingestion of prions (horizontal route) and through maternal route. It is commonly accepted (see, e.g., [21]) that the infectious phase (including the clinical state) lasts at most one year, while the latent state can be of several years. If one compares the epidemic data with the total population of around 9 million cattle in Great Britain, it is reasonable to assume that, even at the peak time, BSE may be considered as a rare disease in a large population. We moreover assume that the probability for an animal to become infected by a horizontal route follows a Reed-Frost type model as in (ii) of Section 3.1, and that the probability for a calf to become infected by its dam is nonnull and constant over time. Moreover, for the sake of simplicity and parsimony, we do not take into account potential heterogeneity factors such as the different regulations from 1989 (the main regulation was the feed ban of July 1988 and was shown to be quite efficient [22]), the different types of breeding and of races, the age for animals older than one year, and the evolution of the surveillance system and of diagnosis tests. So the process {Xn}n defined by (3.1)–(3.4) appears relevant to model the spread of this fatal disease after the 1988 ban, where the infectious state I corresponds to the end of the incubation period and the clinical state, and R corresponds to the death (assumed to be either by routine slaughtering for E and S cattle, and control slaughtering for I cattle). Choosing a time step of one year, Xn then represents the yearly incidence of notified cases, which corresponds to the available data provided by the World Organisation for Animal Health [5].


### 3.2.2. Choice of the Model Parameters

In order to predict the future epidemic evolution, we need to evaluate all the parameters involved in (3.3). These parameters were already studied in previous works [22–25] using different models, amounts of observations, and kinds of estimators. Since it might have a crucial role in the decay phase of the disease, we focus here on θ0 the infection parameter from 1989 via a horizontal route of transmission (i.e., the mean number per infective and per year of newly infected). Until the 1988 feed ban regulation, the main routes of transmission of BSE were horizontal via protein supplements (Meat and Bone Meal, milk replacers), and maternal from a dam to its calf. Since a previous statistical study [22] concluded to the full efficiency of the 1988 ban and since most of cattle are slaughtered before the age of 10 years, the fact that cases of BSE are still observed more than 20 years later could suggest the existence of a remaining source of infection, either via a maternal transmission route, or via a horizontal one, for example, via the ingestion of excreted prions from alive infected animals, or from prions left in the environment. The prediction of the future disease spread thus strongly relies on the intensity of this infection, quantified by the parameter θ0. We provide an estimation of θ0 together with a confidence interval, which will be used for prediction estimations. We provide in addition a sensitivity analysis to the other parameters of the model, which are chosen as follows. First, the maximal age of the individuals in the population am is set to am=10, which corresponds to the largest reported survival age with a nonnegligible probability. Next, we set pmat=0.1, which is the largest order of magnitude commonly accepted for the maternal infection probability [23, 25]. Following a previous work, [22], we assume here a discretized Weibull distribution with parameters α,β for the distribution of the latent period; that is to say, for each k⩾1, Pinc(k):=e-((α-1)/αβα)(k-1)α-e-((α-1)/αβα)kα, where α is a shape parameter, and β is the mode of the probability density of the corresponding continuous Weibull distribution. The parameters α and β have a very bad identifiability with the infection parameters on a sole monotonous phase and are consequently estimated in [22] on the whole epidemic series (growth and decay), by Bayesian maximum a posteriori (MAP) estimations. We set α=α^MAP=3.84 and β=β^MAP=7.46. Since the death of susceptible and exposed animals is mostly due to routine slaughtering, the survival probability S(k) (probability for an apparently healthy animal to survive at least until k years) corresponds to the probability of being slaughtered after the age of k years. It is derived from existing literature [24].


### 3.2.3. Estimation of the Infection Parameter

Our study is based on the yearly number of cases of BSE reported in Great Britain from 1989 until 2011 (Figure1) provided by the World Organisation for Animal Health [5], denoted by Xnobs for each year n. We set d=am-1=9 and denote by Xnobs the d-dimensional vector (Xnobs,…,Xn-d+1obs). We choose the first time n=1 of the epidemic model such that the model with its initial values covers the period starting from 1989, which we consider here as a time-homogeneous period. Hence X-d+1=X1989obs, and n=1 corresponds to the year 1989+d=1998. We estimate θ0 by the WCLSE (2.30) in model (1.1)–(3.3), where {pmat,α,β,{S(k)}k} are given by the previous values, and where X0=X1997obs. According to [5], |X0|=|X1997obs|=167977. We are thus close to the asymptotic |X0|→∞. The number of observations is n=14. The estimator (2.30) provides the estimation θ^∣X0∣obs=2.4324. We point out that this estimation is of the same order of magnitude as the maximum a posteriori Bayesian estimation θ^MAP=2.43 based on the whole epidemic until 2007, assuming a uniform prior probability [22]. Using (2.32) we obtain the following confidence interval [θ^min,θ^max] with asymptotic probability 95%, where θ^min:=θ^∣X0∣-1.96c^1-1, θ^max:=θ^∣X0∣+1.96c^1-1, and c^1:=(∑k=1n‍a·Xk-1/σ2(θ^∣X0∣))1/2. Assuming αi=X1997-i+1obs|X1997obs|-1, we get c^1obs=40.7343 (observed value of c^1), and
(3.5)ℙ(θ0∈[θ^min,θ^max])≃95%,θ^minobs=2.3842,θ^maxobs=2.4805.
Although this confidence interval is an asymptotic one, as |X0|→∞, it is a very good approximation of the true confidence interval for a finite |X0|, since |X0| is here very large. Since the estimation of θ0 relies on the values given to the other model parameters {pmat,α,β}, we evaluate in addition its sensitivity to the values of these parameters. For this purpose, we compute the estimation of θ0 and the associated confidence interval with asymptotic probability 95%, for different values of (pmat,α,β) (Table 1). The first line of the table corresponds to the parameters chosen for the model. In each of the four following lines, we fix two coordinates and choose an extremal (unrealistic) value for the third one. It appears that the estimation of θ0 is almost independent of the value of the maternal infection parameter. However, the estimation seems more strongly dependent on the parameters of the latent period distribution. Nevertheless, even for very unrealistic values (α,β), all the estimations of θ0 remain in the same order of magnitude of several units. This is really small compared to estimations obtained for the infection via Meat and Bone Meal or lactoreplacers (before 1989) which are of the order of 1000 [22]. However, although these estimations are all very small, θ0 seems nonnull. This could suggest the existence of a minor but nonnull infection source which is not of maternal type.Table 1
Sensitivity analysis. Estimation of the infection parameter θ0, and its confidence interval [θ^minobs,θ^maxobs] with asymptotic probability 95%, for different values of the maternal infection parameter pmat and of the latency parameters (α,β). The values (3.84,7.46) correspond to the Bayesian MAP estimations (α^MAP,β^MAP), and 0.1 to the largest commonly accepted order of magnitude for pmat. The estimations of θ0 are based on the observed data over the years 1989–2011.
p
mat
α
β
θ
^
|
X
0
|
obs
[
θ
^
min
obs
,
θ
^
max
obs
]
0.1
3.84
7.46
2.4324
[
2.3842,2.4805
]
0
3.84
7.46
2.4860
[
2.4379,2.5342
]
1
3.84
7.46
1.9492
[
1.9014,1.9970
]
0.1
2
7.46
2.7835
[
2.7287,2.8382
]
0.1
20
7.46
4.0186
[
3.9395,4.0977
]
0.1
3.84
1
1.0127
[
0.9925,1.0329
]
0.1
3.84
10
6.2128
[
6.0914,6.3341
]
0.1
3
6
1.5402
[
1.5095,1.5710
]
0.1
4
5
1.0227
[
1.0020,1.0434
]Figure 1
Yearly number of cases of BSE reported in Great Britain from 1987 to 2011 [5]. Our study only takes into account the data from 1989, after the 1988 feed ban law.


### 3.2.4. Extinction of the Epidemic

We know thanks to Proposition2.1 that {Xn}n becomes extinct almost surely if and only if R0=∑k=1d‍Ψk(θ0)⩽1. The estimated basic reproduction is here R0(θ^∣X0∣obs)=0.1072. Moreover, solving with a computing program the equation ∑k=1d‍Ψk(θ^∣X0∣obs)ρ-k=1, we obtain the following value for the Perron's root ρ(θ^∣X0∣obs)=0.6665, which provides the speed of decay of the expected yearly incidence of cases (see (2.6)); from a certain time, the expected number of new cases will decrease from around 33% every year.


### 3.2.5. Prediction of the Incidences of Cases and Incidences of Infected Cattle

Let us predict the spread of the disease from 2012 by means of simulations of{Xn}n, where θ0 is replaced by its previous estimation θ^∣X0∣obs=2.4324, and where the initial time of the model is 2011, that is, X0=X2011obs. The simulations are done recursively using the transition law (3.1). We point out that the model initialized by X0=X1997obs provides quite realistic simulations on the period 1998–2011 compared to the real observations on the same period, as illustrated in Figure 2(a). The epidemic process {Xn}n thus seems to provide a satisfying prediction of the overall evolution of the real epidemic. In order to predict the incidences of futures cases, we simulate 1000 trajectories of {Xn}n initialized by the observed values X2011obs, with the estimated infection parameter θ^∣X0∣obs=2.4324. We illustrate in Figure 2(c), for each year from 2012, the maximum, minimum, median, 2.5%, and 97.5% quantiles associated with these 1000 realizations. It is also relevant to study and predict the evolution of the incidence of infected cattle in the population, which represents the hidden face of the epidemic. The incidence Zn of infected cattle at time n, conditionally on the number Xn of cases at that time, is given by the Poisson distribution (3.2). For every n⩾2012 and for each of the 1000 previously simulated values Xn, we generate one realization of Zn. We then illustrate in Figure 2(d), the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations.Figure2(a) 10 simulations of {Xn}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 2(b) 5 simulations of {Xn}n initialized by X2011obs. Figures 2(c) and 2(d) prediction, based on 1000 simulations of the process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣obs=2.4324.
(a)
(b)
(c)
(d)


### 3.2.6. Prediction of the Year of Extinction

LetT:=2011+inf{n⩾1,Xn=0} denote the extinction year of the epidemic process {Xn}n. According to (2.7) and denoting by fθ0 the offspring generating function of {Xn}n defined in (2.2) (from now on we let the dependence in θ0 appear in the notation), we have ℙθ0(T⩽2011+n)=fθ0,n(0)X2011obs, for every n⩾1, which by iterating fθ0 can be computed explicitly. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the extinction time (see (2.8)): n0.5T=2028, n0.95T=2035 and n0.99T=2039. Keeping in mind that T corresponds to the complete extinction of the epidemic (i.e., d=9 consecutive years without any case), these results actually mean that for the infection parameter θ0=2.4324, with probability larger than 50% (resp., 95% and 99%), no case will arise in the population from year 2020 (resp., 2027 and 2031). Moreover, in order to take into account the uncertainty around the estimation θ^∣X0∣obs of the infection parameter θ0, we make use of the asymptotic confidence interval (3.5) of θ0 and of the fact that θ↦ℙθ(T⩽n) is a decreasing function of θ, which implies that for every n⩾2011,
(3.6)ℙ(ℙθ0(T⩽n)∈[ℙθ^max(T⩽n),ℙθ^min(T⩽n)])≃95%.
We collect in Table 2 the observed interval [ℙθ^maxobs(T⩽n),ℙθ^minobs(T⩽n)], for each n⩾2020 (if n<2019 we have, conditionally on the initial value X2011obs, ℙ(Xn=0)=0 because of the memory which is not equal to 0). Note that these intervals are very narrow, leading to an accurate estimation of ℙθ0(T⩽n).Table 2
Cumulative distribution function of the year of extinction computed with the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npT for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(T⩽npT) based on (3.6).
n
ℙ
θ
^
max
obs
(
T
⩽
n
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
2020
0.0000
0.0000
2030
0.7610
0.7807
2040
0.9950
0.9957
2021
0.0000
0.0000
2031
0.8310
0.8465
2041
0.9967
0.9972
2022
0.0010
0.0014
2032
0.8816
0.8934
2042
0.9978
0.9981
2023
0.0121
0.0152
2033
0.9186
0.9273
2043
0.9985
0.9988
2024
0.0496
0.0592
2034
0.9451
0.9513
2044
0.9990
0.9992
2025
0.1211
0.1390
2035
0.9633
0.9677
2045
0.9995
0.9994
2026
0.2325
0.2579
2036
0.9755
0.9786
2046
0.9996
0.9996
2027
0.3756
0.4047
2037
0.9835
0.9857
2047
0.9997
0.9998
2028
0.5303
0.5584
2038
0.9889
0.9904
2048
0.9998
0.9998
2029
0.6619
0.6862
2039
0.9925
0.9936
2049
0.9999
0.9999


### 3.2.7. Prediction of the Epidemic Size

LetN:=∑n=1T-2011‍Xn be the total size of the future epidemic from 2012 (total number of cases from 2012 until the extinction of the epidemic). We compute the distribution of N using (2.11), conditionally on the event {X0=X2011obs}. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the total epidemic size (see (2.12)): n0.5N=22, n0.95N=31 and n0.99N=35. From (2.13) we deduce that the mean value and variance of N for the parameter θ^minobs (resp., θ^maxobs) are 21 and 27 (resp., 22 and 28). Moreover, we obtain similarly as for the extinction time and thanks to (3.5) a confidence interval [ℙθ^max(N⩽n),ℙθ^min(N⩽n)] of ℙθ0(N⩽n), for every n∈ℕ, with asymptotic probability 95%, collected in Table 3.Table 3
Cumulative distribution function of the total size of the epidemic for the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npN for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(N⩽npN).
n
ℙ
θ
^
max
obs
(
N
⩽
n
)
ℙ
θ
^
min
obs
(
N
⩽
n
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
⋯
⋯
⋯
16
0.1293
0.1669
28
0.8702
0.9048
5
0.0001
0.0000
17
0.1795
0.2263
29
0.8942
0.9302
6
0.0002
0.0003
18
0.2373
0.2915
30
0.9283
0.9495
7
0.0005
0.0008
19
0.3025
0.3637
31
0.9480
0.9646
8
0.0013
0.0020
20
0.3741
0.4409
32
0.9635
0.9759
9
0.0027
0.0046
21
0.4486
0.5189
33
0.9742
0.9838
10
0.0060
0.0094
22
0.5231
0.5934
34
0.9824
0.9892
11
0.0119
0.0178
23
0.5975
0.6639
35
0.9881
0.9928
12
0.0214
0.0310
24
0.6655
0.7287
36
0.9917
0.9952
13
0.0365
0.0518
25
0.7282
0.7840
37
0.9945
0.9968
14
0.0585
0.0801
26
0.7829
0.8323
38
0.9964
0.9980
15
0.0891
0.1181
27
0.8304
0.8726
39
0.9977
0.9989


### 3.2.8. Study of the Very Late Extinction Case

In order to predict the behavior of the “most dangerous” evolution of the epidemic, we first need to compute the estimationθ^n*obs (2.47) based on the data in Great Britain. The number of available observations is only n=14, so we are far from the asymptotic setting n→∞ of Theorem 2.7. However, the large value of |X0| can make us hope for a good accuracy. We point out that, by making use of the estimator θ^n* on the real data, we make an unverifiable assumption on the future of the epidemic; we consider the data as if they were the beginning of a trajectory with very late extinction. This should have the following consequence: the estimation provided by θ^n*obs should a priori be a bit smaller than the value 2.4324 provided by θ^∣X0∣obs. Indeed we obtain θ^n*obs=2.4305. Using this value, we deduce from (2.51) the confidence interval [θ^min*,θ^max*] of θ0 with asymptotic probability 95%, where θ^min*:=θ^n*-1.96c^2-1, θ^max*:=θ^n*+1.96c^2-1, and
(3.7)c^2:=∑k=0n‍(f′(θ^n*,Xk*))2[∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)]-1/2.
Therefore, c^2obs=40.6988, and
(3.8)ℙ(θ0∈[θ^min*,θ^max*])≃95%,θ^min*obs=2.3823,θ^max*obs=2.4787,
which is of the same magnitude order as the confidence interval [2.3842,2.4805] obtained with the unconditioned process (see (3.5)). Let us predict the most dangerous evolution thanks to the transition law of the conditioned process given by Proposition 2.2. Note that if one computes the eigenvalues ρ and λ introduced in Proposition 2.5, one obtains ρ(θ^*obs)=0.6664 and |λ(θ^*obs)|=0.5570. It thus appears that the convergence of the epidemic process conditioned on nonextinction at time k to the Q-process as k→∞ is not very fast. As a consequence the study of the Q-process only provides information about the behavior of the disease spread in the case of an extremely late extinction. First, we see thanks to Figure 3(a) that the simulations provided by the conditioned process initialized by X0=X1997obs, and where θ0 is estimated by θ^n*obs=2.4305, are true to the real observations on the period 1999–2011. Figure 3(b) is an example of one simulation in the period 2012–2040 of the conditioned process, for X0*=X2011obs. It appears that the values of this simulated trajectory are rapidly very small and of course are never equal to 0 for d=9 consecutive times. For a finer prediction, we simulate 1000 realizations of this process from 2012 until 2050, with X0*=X2011obs. Moreover, for every n⩾2012 and for each of the 1000 simulated values Xn*, we make one realization of the incidence Zn of infected cattle at time n, according to the law given by (3.2). Figures 3(c) and 3(d) represent the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations of, respectively, the incidence of cases and infected cattle, in case of an extreme late extinction. It appears thanks to Figures 3(c) and 3(d) that the supposedly “most dangerous” trajectories nevertheless do not reach high values and do not present a new peak of epidemic.Figure3(a) 10 simulations of {Xn*}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 3(b) one simulation of {Xn*}n initialized by X2011obs. Figures 3(c) and 3(d) prediction, based on 1000 simulations of the conditioned process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣*obs=2.4305.
(a)
(b)
(c)
(d)


## 3.2.1. Choice of the Epidemic Model

BSE is a transmissible disease through ingestion of prions (horizontal route) and through maternal route. It is commonly accepted (see, e.g., [21]) that the infectious phase (including the clinical state) lasts at most one year, while the latent state can be of several years. If one compares the epidemic data with the total population of around 9 million cattle in Great Britain, it is reasonable to assume that, even at the peak time, BSE may be considered as a rare disease in a large population. We moreover assume that the probability for an animal to become infected by a horizontal route follows a Reed-Frost type model as in (ii) of Section 3.1, and that the probability for a calf to become infected by its dam is nonnull and constant over time. Moreover, for the sake of simplicity and parsimony, we do not take into account potential heterogeneity factors such as the different regulations from 1989 (the main regulation was the feed ban of July 1988 and was shown to be quite efficient [22]), the different types of breeding and of races, the age for animals older than one year, and the evolution of the surveillance system and of diagnosis tests. So the process {Xn}n defined by (3.1)–(3.4) appears relevant to model the spread of this fatal disease after the 1988 ban, where the infectious state I corresponds to the end of the incubation period and the clinical state, and R corresponds to the death (assumed to be either by routine slaughtering for E and S cattle, and control slaughtering for I cattle). Choosing a time step of one year, Xn then represents the yearly incidence of notified cases, which corresponds to the available data provided by the World Organisation for Animal Health [5].


## 3.2.2. Choice of the Model Parameters

In order to predict the future epidemic evolution, we need to evaluate all the parameters involved in (3.3). These parameters were already studied in previous works [22–25] using different models, amounts of observations, and kinds of estimators. Since it might have a crucial role in the decay phase of the disease, we focus here on θ0 the infection parameter from 1989 via a horizontal route of transmission (i.e., the mean number per infective and per year of newly infected). Until the 1988 feed ban regulation, the main routes of transmission of BSE were horizontal via protein supplements (Meat and Bone Meal, milk replacers), and maternal from a dam to its calf. Since a previous statistical study [22] concluded to the full efficiency of the 1988 ban and since most of cattle are slaughtered before the age of 10 years, the fact that cases of BSE are still observed more than 20 years later could suggest the existence of a remaining source of infection, either via a maternal transmission route, or via a horizontal one, for example, via the ingestion of excreted prions from alive infected animals, or from prions left in the environment. The prediction of the future disease spread thus strongly relies on the intensity of this infection, quantified by the parameter θ0. We provide an estimation of θ0 together with a confidence interval, which will be used for prediction estimations. We provide in addition a sensitivity analysis to the other parameters of the model, which are chosen as follows. First, the maximal age of the individuals in the population am is set to am=10, which corresponds to the largest reported survival age with a nonnegligible probability. Next, we set pmat=0.1, which is the largest order of magnitude commonly accepted for the maternal infection probability [23, 25]. Following a previous work, [22], we assume here a discretized Weibull distribution with parameters α,β for the distribution of the latent period; that is to say, for each k⩾1, Pinc(k):=e-((α-1)/αβα)(k-1)α-e-((α-1)/αβα)kα, where α is a shape parameter, and β is the mode of the probability density of the corresponding continuous Weibull distribution. The parameters α and β have a very bad identifiability with the infection parameters on a sole monotonous phase and are consequently estimated in [22] on the whole epidemic series (growth and decay), by Bayesian maximum a posteriori (MAP) estimations. We set α=α^MAP=3.84 and β=β^MAP=7.46. Since the death of susceptible and exposed animals is mostly due to routine slaughtering, the survival probability S(k) (probability for an apparently healthy animal to survive at least until k years) corresponds to the probability of being slaughtered after the age of k years. It is derived from existing literature [24].


## 3.2.3. Estimation of the Infection Parameter

Our study is based on the yearly number of cases of BSE reported in Great Britain from 1989 until 2011 (Figure1) provided by the World Organisation for Animal Health [5], denoted by Xnobs for each year n. We set d=am-1=9 and denote by Xnobs the d-dimensional vector (Xnobs,…,Xn-d+1obs). We choose the first time n=1 of the epidemic model such that the model with its initial values covers the period starting from 1989, which we consider here as a time-homogeneous period. Hence X-d+1=X1989obs, and n=1 corresponds to the year 1989+d=1998. We estimate θ0 by the WCLSE (2.30) in model (1.1)–(3.3), where {pmat,α,β,{S(k)}k} are given by the previous values, and where X0=X1997obs. According to [5], |X0|=|X1997obs|=167977. We are thus close to the asymptotic |X0|→∞. The number of observations is n=14. The estimator (2.30) provides the estimation θ^∣X0∣obs=2.4324. We point out that this estimation is of the same order of magnitude as the maximum a posteriori Bayesian estimation θ^MAP=2.43 based on the whole epidemic until 2007, assuming a uniform prior probability [22]. Using (2.32) we obtain the following confidence interval [θ^min,θ^max] with asymptotic probability 95%, where θ^min:=θ^∣X0∣-1.96c^1-1, θ^max:=θ^∣X0∣+1.96c^1-1, and c^1:=(∑k=1n‍a·Xk-1/σ2(θ^∣X0∣))1/2. Assuming αi=X1997-i+1obs|X1997obs|-1, we get c^1obs=40.7343 (observed value of c^1), and
(3.5)ℙ(θ0∈[θ^min,θ^max])≃95%,θ^minobs=2.3842,θ^maxobs=2.4805.
Although this confidence interval is an asymptotic one, as |X0|→∞, it is a very good approximation of the true confidence interval for a finite |X0|, since |X0| is here very large. Since the estimation of θ0 relies on the values given to the other model parameters {pmat,α,β}, we evaluate in addition its sensitivity to the values of these parameters. For this purpose, we compute the estimation of θ0 and the associated confidence interval with asymptotic probability 95%, for different values of (pmat,α,β) (Table 1). The first line of the table corresponds to the parameters chosen for the model. In each of the four following lines, we fix two coordinates and choose an extremal (unrealistic) value for the third one. It appears that the estimation of θ0 is almost independent of the value of the maternal infection parameter. However, the estimation seems more strongly dependent on the parameters of the latent period distribution. Nevertheless, even for very unrealistic values (α,β), all the estimations of θ0 remain in the same order of magnitude of several units. This is really small compared to estimations obtained for the infection via Meat and Bone Meal or lactoreplacers (before 1989) which are of the order of 1000 [22]. However, although these estimations are all very small, θ0 seems nonnull. This could suggest the existence of a minor but nonnull infection source which is not of maternal type.Table 1
Sensitivity analysis. Estimation of the infection parameter θ0, and its confidence interval [θ^minobs,θ^maxobs] with asymptotic probability 95%, for different values of the maternal infection parameter pmat and of the latency parameters (α,β). The values (3.84,7.46) correspond to the Bayesian MAP estimations (α^MAP,β^MAP), and 0.1 to the largest commonly accepted order of magnitude for pmat. The estimations of θ0 are based on the observed data over the years 1989–2011.
p
mat
α
β
θ
^
|
X
0
|
obs
[
θ
^
min
obs
,
θ
^
max
obs
]
0.1
3.84
7.46
2.4324
[
2.3842,2.4805
]
0
3.84
7.46
2.4860
[
2.4379,2.5342
]
1
3.84
7.46
1.9492
[
1.9014,1.9970
]
0.1
2
7.46
2.7835
[
2.7287,2.8382
]
0.1
20
7.46
4.0186
[
3.9395,4.0977
]
0.1
3.84
1
1.0127
[
0.9925,1.0329
]
0.1
3.84
10
6.2128
[
6.0914,6.3341
]
0.1
3
6
1.5402
[
1.5095,1.5710
]
0.1
4
5
1.0227
[
1.0020,1.0434
]Figure 1
Yearly number of cases of BSE reported in Great Britain from 1987 to 2011 [5]. Our study only takes into account the data from 1989, after the 1988 feed ban law.


## 3.2.4. Extinction of the Epidemic

We know thanks to Proposition2.1 that {Xn}n becomes extinct almost surely if and only if R0=∑k=1d‍Ψk(θ0)⩽1. The estimated basic reproduction is here R0(θ^∣X0∣obs)=0.1072. Moreover, solving with a computing program the equation ∑k=1d‍Ψk(θ^∣X0∣obs)ρ-k=1, we obtain the following value for the Perron's root ρ(θ^∣X0∣obs)=0.6665, which provides the speed of decay of the expected yearly incidence of cases (see (2.6)); from a certain time, the expected number of new cases will decrease from around 33% every year.


## 3.2.5. Prediction of the Incidences of Cases and Incidences of Infected Cattle

Let us predict the spread of the disease from 2012 by means of simulations of{Xn}n, where θ0 is replaced by its previous estimation θ^∣X0∣obs=2.4324, and where the initial time of the model is 2011, that is, X0=X2011obs. The simulations are done recursively using the transition law (3.1). We point out that the model initialized by X0=X1997obs provides quite realistic simulations on the period 1998–2011 compared to the real observations on the same period, as illustrated in Figure 2(a). The epidemic process {Xn}n thus seems to provide a satisfying prediction of the overall evolution of the real epidemic. In order to predict the incidences of futures cases, we simulate 1000 trajectories of {Xn}n initialized by the observed values X2011obs, with the estimated infection parameter θ^∣X0∣obs=2.4324. We illustrate in Figure 2(c), for each year from 2012, the maximum, minimum, median, 2.5%, and 97.5% quantiles associated with these 1000 realizations. It is also relevant to study and predict the evolution of the incidence of infected cattle in the population, which represents the hidden face of the epidemic. The incidence Zn of infected cattle at time n, conditionally on the number Xn of cases at that time, is given by the Poisson distribution (3.2). For every n⩾2012 and for each of the 1000 previously simulated values Xn, we generate one realization of Zn. We then illustrate in Figure 2(d), the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations.Figure2(a) 10 simulations of {Xn}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 2(b) 5 simulations of {Xn}n initialized by X2011obs. Figures 2(c) and 2(d) prediction, based on 1000 simulations of the process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣obs=2.4324.
(a)
(b)
(c)
(d)


## 3.2.6. Prediction of the Year of Extinction

LetT:=2011+inf{n⩾1,Xn=0} denote the extinction year of the epidemic process {Xn}n. According to (2.7) and denoting by fθ0 the offspring generating function of {Xn}n defined in (2.2) (from now on we let the dependence in θ0 appear in the notation), we have ℙθ0(T⩽2011+n)=fθ0,n(0)X2011obs, for every n⩾1, which by iterating fθ0 can be computed explicitly. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the extinction time (see (2.8)): n0.5T=2028, n0.95T=2035 and n0.99T=2039. Keeping in mind that T corresponds to the complete extinction of the epidemic (i.e., d=9 consecutive years without any case), these results actually mean that for the infection parameter θ0=2.4324, with probability larger than 50% (resp., 95% and 99%), no case will arise in the population from year 2020 (resp., 2027 and 2031). Moreover, in order to take into account the uncertainty around the estimation θ^∣X0∣obs of the infection parameter θ0, we make use of the asymptotic confidence interval (3.5) of θ0 and of the fact that θ↦ℙθ(T⩽n) is a decreasing function of θ, which implies that for every n⩾2011,
(3.6)ℙ(ℙθ0(T⩽n)∈[ℙθ^max(T⩽n),ℙθ^min(T⩽n)])≃95%.
We collect in Table 2 the observed interval [ℙθ^maxobs(T⩽n),ℙθ^minobs(T⩽n)], for each n⩾2020 (if n<2019 we have, conditionally on the initial value X2011obs, ℙ(Xn=0)=0 because of the memory which is not equal to 0). Note that these intervals are very narrow, leading to an accurate estimation of ℙθ0(T⩽n).Table 2
Cumulative distribution function of the year of extinction computed with the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npT for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(T⩽npT) based on (3.6).
n
ℙ
θ
^
max
obs
(
T
⩽
n
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
2020
0.0000
0.0000
2030
0.7610
0.7807
2040
0.9950
0.9957
2021
0.0000
0.0000
2031
0.8310
0.8465
2041
0.9967
0.9972
2022
0.0010
0.0014
2032
0.8816
0.8934
2042
0.9978
0.9981
2023
0.0121
0.0152
2033
0.9186
0.9273
2043
0.9985
0.9988
2024
0.0496
0.0592
2034
0.9451
0.9513
2044
0.9990
0.9992
2025
0.1211
0.1390
2035
0.9633
0.9677
2045
0.9995
0.9994
2026
0.2325
0.2579
2036
0.9755
0.9786
2046
0.9996
0.9996
2027
0.3756
0.4047
2037
0.9835
0.9857
2047
0.9997
0.9998
2028
0.5303
0.5584
2038
0.9889
0.9904
2048
0.9998
0.9998
2029
0.6619
0.6862
2039
0.9925
0.9936
2049
0.9999
0.9999


## 3.2.7. Prediction of the Epidemic Size

LetN:=∑n=1T-2011‍Xn be the total size of the future epidemic from 2012 (total number of cases from 2012 until the extinction of the epidemic). We compute the distribution of N using (2.11), conditionally on the event {X0=X2011obs}. We obtain in particular, for the estimated value θ^∣X0∣obs=2.4324, the following p-quantiles for the total epidemic size (see (2.12)): n0.5N=22, n0.95N=31 and n0.99N=35. From (2.13) we deduce that the mean value and variance of N for the parameter θ^minobs (resp., θ^maxobs) are 21 and 27 (resp., 22 and 28). Moreover, we obtain similarly as for the extinction time and thanks to (3.5) a confidence interval [ℙθ^max(N⩽n),ℙθ^min(N⩽n)] of ℙθ0(N⩽n), for every n∈ℕ, with asymptotic probability 95%, collected in Table 3.Table 3
Cumulative distribution function of the total size of the epidemic for the infection parametersθ^minobs=2.3842 and θ^maxobs=2.4805 defined by (3.5). The values in bold character correspond to the p-quantiles npN for p=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙθ0(N⩽npN).
n
ℙ
θ
^
max
obs
(
N
⩽
n
)
ℙ
θ
^
min
obs
(
N
⩽
n
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
n
ℙ
θ
^
max
obs
(
⋯
)
ℙ
θ
^
min
obs
(
⋯
)
⋯
⋯
⋯
16
0.1293
0.1669
28
0.8702
0.9048
5
0.0001
0.0000
17
0.1795
0.2263
29
0.8942
0.9302
6
0.0002
0.0003
18
0.2373
0.2915
30
0.9283
0.9495
7
0.0005
0.0008
19
0.3025
0.3637
31
0.9480
0.9646
8
0.0013
0.0020
20
0.3741
0.4409
32
0.9635
0.9759
9
0.0027
0.0046
21
0.4486
0.5189
33
0.9742
0.9838
10
0.0060
0.0094
22
0.5231
0.5934
34
0.9824
0.9892
11
0.0119
0.0178
23
0.5975
0.6639
35
0.9881
0.9928
12
0.0214
0.0310
24
0.6655
0.7287
36
0.9917
0.9952
13
0.0365
0.0518
25
0.7282
0.7840
37
0.9945
0.9968
14
0.0585
0.0801
26
0.7829
0.8323
38
0.9964
0.9980
15
0.0891
0.1181
27
0.8304
0.8726
39
0.9977
0.9989


## 3.2.8. Study of the Very Late Extinction Case

In order to predict the behavior of the “most dangerous” evolution of the epidemic, we first need to compute the estimationθ^n*obs (2.47) based on the data in Great Britain. The number of available observations is only n=14, so we are far from the asymptotic setting n→∞ of Theorem 2.7. However, the large value of |X0| can make us hope for a good accuracy. We point out that, by making use of the estimator θ^n* on the real data, we make an unverifiable assumption on the future of the epidemic; we consider the data as if they were the beginning of a trajectory with very late extinction. This should have the following consequence: the estimation provided by θ^n*obs should a priori be a bit smaller than the value 2.4324 provided by θ^∣X0∣obs. Indeed we obtain θ^n*obs=2.4305. Using this value, we deduce from (2.51) the confidence interval [θ^min*,θ^max*] of θ0 with asymptotic probability 95%, where θ^min*:=θ^n*-1.96c^2-1, θ^max*:=θ^n*+1.96c^2-1, and
(3.7)c^2:=∑k=0n‍(f′(θ^n*,Xk*))2[∑k=0n‍(f′(θ^n*,Xk*))2g(θ^n*,Xk*)]-1/2.
Therefore, c^2obs=40.6988, and
(3.8)ℙ(θ0∈[θ^min*,θ^max*])≃95%,θ^min*obs=2.3823,θ^max*obs=2.4787,
which is of the same magnitude order as the confidence interval [2.3842,2.4805] obtained with the unconditioned process (see (3.5)). Let us predict the most dangerous evolution thanks to the transition law of the conditioned process given by Proposition 2.2. Note that if one computes the eigenvalues ρ and λ introduced in Proposition 2.5, one obtains ρ(θ^*obs)=0.6664 and |λ(θ^*obs)|=0.5570. It thus appears that the convergence of the epidemic process conditioned on nonextinction at time k to the Q-process as k→∞ is not very fast. As a consequence the study of the Q-process only provides information about the behavior of the disease spread in the case of an extremely late extinction. First, we see thanks to Figure 3(a) that the simulations provided by the conditioned process initialized by X0=X1997obs, and where θ0 is estimated by θ^n*obs=2.4305, are true to the real observations on the period 1999–2011. Figure 3(b) is an example of one simulation in the period 2012–2040 of the conditioned process, for X0*=X2011obs. It appears that the values of this simulated trajectory are rapidly very small and of course are never equal to 0 for d=9 consecutive times. For a finer prediction, we simulate 1000 realizations of this process from 2012 until 2050, with X0*=X2011obs. Moreover, for every n⩾2012 and for each of the 1000 simulated values Xn*, we make one realization of the incidence Zn of infected cattle at time n, according to the law given by (3.2). Figures 3(c) and 3(d) represent the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations of, respectively, the incidence of cases and infected cattle, in case of an extreme late extinction. It appears thanks to Figures 3(c) and 3(d) that the supposedly “most dangerous” trajectories nevertheless do not reach high values and do not present a new peak of epidemic.Figure3(a) 10 simulations of {Xn*}n initialized by X1997obs, and comparison with the observations on the period 1998–2011. Figure 3(b) one simulation of {Xn*}n initialized by X2011obs. Figures 3(c) and 3(d) prediction, based on 1000 simulations of the conditioned process of the yearly incidences of cases (resp., infected cattle) from 2012. 95% of the trajectories remain in the band delimited by the blue dotted lines. All the simulations are done with the infection parameter θ^∣X0∣*obs=2.4305.
(a)
(b)
(c)
(d)


## 4. Conclusion

The methodology developed in Section2 thus applies well here and enables not only to confirm mathematically what is commonly accepted, namely, that BSE is fading out, but also to predict, with a very large probability, that the last BSE case will occur before 2027, and that until the complete extinction of the epidemic in the population, there will be less than 31 cases to come. We obtain moreover from Figure 2(d) the order of magnitude of the number of infected cattle in the population. In addition, the estimation of the infection parameter concludes to the possible existence of a minor but nonnull infection source which is not of maternal type, and which is very small (only around 3 newly infected animal per year and per infective) compared to the main source of horizontal infection until 1988 due to protein supplements. Finally, the study of the worst-case scenario shows that even in the case of an extreme late extinction of the disease in the population, the incidence of cases will decrease quite rapidly to 0, with afterwards only 1 or 2 yearly cases occurring regularly but sparsely, with no appearance of a new peak of epidemic. We have shown with this example that the methodology developed in Section 2 provides accurate tools to study the decay phase of an epidemic under the current sanitary measures, which would help to make new policy decisions. This evaluation is all the more relevant since it is obtained not by simply computing what should most probably happen, but also by taking into account the variability of many factors (infection, incubation, and survival), and by studying the potentially most dangerous evolution.



---

*Source: 598701-2012-12-04.xml*