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

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 the Q-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.

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 {𝑋_𝑛}_𝑛 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”) 𝑘 time units later, where 1⩽𝑘⩽𝑑. We assume that this number of secondary cases follows a Poisson probability distribution with parameter Ψ_𝑘>0. The recursive formula defining {𝑋_𝑛}_𝑛 is then the following: 𝑋_𝑛=^𝑑𝑋_𝑘=1^𝑛−𝑘_𝑖=1𝑌_{𝑛−𝑘,𝑛,𝑖}, (1.1) where the variable 𝑌_{𝑛−𝑘,𝑛,𝑖} is the incidence of secondary cases produced at time 𝑛 with a delay 𝑘 (latent period) by individual 𝑖 infectious at time 𝑛−𝑘. The {𝑌_{𝑛−𝑘,𝑛,𝑖}}_𝑖,𝑘 are assumed independent given ℱ_𝑛−1∶=𝜎({𝑋_𝑛−𝑘}_𝑘⩾1), and the {𝑌_{𝑛−𝑘,𝑛,𝑖}}_𝑖 are assumed i.i.d. (identically and independently distributed) given ℱ_𝑛−1, with a common Poisson distribution with parameter Ψ_𝑘. 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 Section 2, where the whole methodology is presented. We first formulate the epidemic model {𝑋_𝑛}_𝑛 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 𝑄-process {𝑋^∗_𝑛}_𝑛, obtained by conditioning {𝑋_𝑛}_𝑛 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 {𝑋_𝑛}_𝑛 and {𝑋^∗_𝑛}_𝑛, as well as the derived distributions above, we need to know the values of the parameters {Ψ_𝑘}_𝑘. We may write Ψ_𝑘=∑^{𝑎𝑚−𝑘}_𝑎=1𝜃_𝑎𝑃_inc,𝑎(𝑘)𝑃_age(𝑎+𝑘), where 𝜃_𝑎 is the mean number of individuals infected at age 𝑎 by an infective by direct or indirect transmission, 𝑎_𝑚 is the largest survival age, 𝑃_inc,𝑎(𝑘) is the probability for the individual aged 𝑎 at infection to have a latent period equal to 𝑘 (given his survival), and 𝑃_age(𝑎+𝑘) is the probability to be aged 𝑎+𝑘 at the end of the latent period. Parameters {𝜃_𝑎}_𝑎 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 𝜃_𝑎=𝜃₀+𝑝_𝑎, where 𝑝_𝑎 is constant over time, while 𝜃₀ may change with control measures. A typical example is when 𝜃₀ is the mean number of individuals infected by an infective by horizontal route at age 𝑎, assumed independent of 𝑎⩾2, and 𝑝₁ represents the maternal transmission probability 𝑝_mat. In this case Ψ_𝑘=𝜃₀^𝑎𝑚_𝑎=𝑘+1𝑃_{inc,𝑎−𝑘}(𝑘)𝑃_age(𝑎)+𝑝_mat𝑃_inc,1(𝑘)𝑃_age(𝑘+1). (1.2) So we assume here that, except for 𝜃₀, the other parameters of the {Ψ_𝑘}_𝑘 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 Ψ_𝑘 depends affinely on 𝜃₀ (see Section 3.1). In Sections 2.5 and 2.6, we provide optimal WCLSE (Weighted Conditional Least Squares Estimators) of 𝜃₀ in the decay phase, in the frame of {𝑋_𝑛}_𝑛 as well as in the frame of the associated conditioned process {𝑋^∗_𝑛}_𝑛, and prove the strong consistency and the asymptotic normality of these estimators.

The final Section 3 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 {Ψ_𝑘}_𝑘. 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 {Ψ_𝑘}_𝑘 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, {𝑋_𝑛}_𝑛 defined by (1.1) is a Markovian process of order 𝑑. Consequently, the 𝑑-dimensional process {𝐗_𝑛}_𝑛 defined by 𝐗_𝑛∶=𝑋_𝑛,𝑋_𝑛−1,…,𝑋_{𝑛−𝑑+1} (2.1) is Markovian of order 1, and it stems directly from (1.1) that {𝐗_𝑛}_𝑛 is a multitype Bienaymé-Galton-Watson (BGW) process with 𝑑 types (see, e.g., [6]). Note that the 𝑑 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 {𝑋_𝑛}_𝑛. The information provided by the 𝑑-dimensional Markovian process {𝐗_𝑛}_𝑛 is therefore the same as the one given by the 1-dimensional 𝑑-Markovian process {𝑋_𝑛}_𝑛, 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 {𝐗_𝑛}_𝑛, 𝐟∶=(𝑓₁,…,𝑓_𝑑), defined on [0,1]^𝑑 by 𝑓_𝑖(𝐫)∶=𝔼[𝐫^𝐗1∣𝐗₀=𝐞_𝑖], where 𝐞_𝑖∶=(0,…,1,…,0) denotes the 𝑖th basis vector of ℕ^𝑑 and 𝐮^𝐯∶=∏^𝑑_𝑖=1𝑢^𝑣𝑖_𝑖 for 𝐮,𝐯∈ℕ^𝑑. For all 𝐫∈[0,1]^𝑑, we have here 𝑓_𝑖(𝐫)∶=𝑒^{−(1−𝑟1)Ψ𝑖}𝑟_𝑖+1𝑓,𝑖=1⋯𝑑−1,_𝑑(𝐫)∶=𝑒^{−(1−𝑟1)Ψ𝑑}. (2.2) The second basic tool is the mean matrix 𝐌 defined by 𝔼(𝐗_𝑛∣ℱ_𝑛−1)=𝐗_𝑛−1𝐌, which is here Ψ𝐌=¹Ψ10⋯0₂Ψ01⋯0⋮⋮⋱⋮_𝑑−1Ψ0⋯⋯1_𝑑0⋯⋯0. (2.3) Let us notice that, since Ψ_𝑘>0 for each 𝑘=1⋯𝑑, then {𝐗_𝑛}_𝑛 is nonsingular, positive regular (see [6]) and satisfies the 𝑋log𝑋 condition, 𝔼𝐗₁ln𝐗₁∣𝐗₀=𝐞_𝑖<∞,𝑖=1,…,𝑑, (2.4) where ‖⋅‖ denotes the sup norm in ℝ^𝑑.

2.1. Extinction of the Epidemic

2.1.1. Almost Sure Extinction

Since the single-type process {𝑋_𝑛}_𝑛 has a memory of size 𝑑, it becomes extinct when it is null at 𝑑 successive times, or equivalently as soon as the 𝑑-dimensional process {𝐗_𝑛}_𝑛 reaches the 𝑑-dimensional null vector 𝟎. According to the theory of multitype positive regular and nonsingular BGW processes [6], the extinction of the process {𝐗_𝑛}_𝑛 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 𝐌. Thus 𝜌 is solution of ∑^𝑑_𝑘=1Ψ_𝑘𝜌^−𝑘=1. In general for 𝑑>1, 𝜌 has no explicit expression. However, ∑^𝑑_𝑘=1Ψ_𝑘𝜌^−𝑘=1 leads directly to the following explicit threshold criteria.

Proposition 2.1. The epidemic becomes extinct almost surely if and only if 𝑅₀⩽1, where 𝑅₀∶=∑^𝑑_𝑘=1Ψ_𝑘 is the total mean number of secondary cases generated by one infective in a 𝑆𝐸𝐼𝑅 disease. We call 𝑅₀ the basic reproduction number.
Moreover, when 𝜌⩽1, then 𝑅₀⩽𝜌 with equality if and only if either 𝜌=1 or 𝑑=1, and when 𝜌>1, then 𝑅₀⩾𝜌 with equality if and only if 𝑑=1.

Note that when 𝑑>1, 𝑅₀ 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 𝑛, for 𝑛 large. Denoting by 𝐮 and 𝐯 the right and left eigenvectors of 𝐌 associated to the Perron's root 𝜌; that is, 𝐌𝐮^𝑇=𝜌𝐮^𝑇 and 𝐯𝐌=𝜌𝐯, with the normalization convention 𝐮⋅𝟏=𝐮⋅𝐯=1, where 𝐮⋅𝐯 stands for the usual scalar product in ℝ^𝑑 and where the superscript 𝑇 denotes the transposition, then 𝔼(𝐗_𝑛∣𝐗₀)=𝐗₀𝐌^𝑛∼_𝑛→∞𝜌^𝑛𝐗₀𝐮^𝑇𝐯. The first coordinate in the latter formula becomes for the epidemic process 𝜌^𝑛∑^𝑑_𝑖=1𝑋_−𝑖+1𝑢_𝑖𝑣₁. Computing explicitly 𝐮 and 𝐯, we obtain that for all 𝑖=1⋯𝑑, 𝑢_𝑖=∑^𝑑_𝑘=𝑖Ψ_𝑘𝜌^−𝑘+𝑖∑^𝑑_𝑗=1∑^𝑑_𝑘=𝑗Ψ_𝑘𝜌^−𝑘+𝑗,𝑣_𝑖=𝜌^−𝑖∑^𝑑_𝑗=1∑^𝑑_𝑘=𝑗Ψ_𝑘𝜌^−𝑘+𝑗∑^𝑑_𝑗=1∑^𝑑_𝑘=𝑗Ψ_𝑘𝜌^−𝑘, (2.5) which leads to the following asymptotic result: 𝔼𝑋_𝑛∣𝐗₀∼_𝑛→∞𝜌^𝑛𝑑_𝑖=1𝑋_−𝑖+1∑^𝑑_𝑘=𝑖Ψ_𝑘𝜌^{−𝑘+𝑖−1}∑^𝑑_𝑗=1∑^𝑑_𝑘=𝑗Ψ_𝑘𝜌^−𝑘. (2.6) 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 value 𝐗₀=(𝑋₀,𝑋₋₁,…,𝑋_−𝑑+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 𝑇∶=inf{𝑛⩾1,𝐗_𝑛=𝟎} denote the extinction time of the process {𝐗_𝑛}_𝑛, and let 𝐟_𝑛∶=𝐟∘𝐟_𝑛−1 be the 𝑛th iterate of the generating function 𝐟 given by (2.2). We denote 𝐟_𝑛∶=(𝑓_𝑛,1,…,𝑓_𝑛,𝑑). Then, by the branching property of the process [6], the probability of extinction of the epidemic before time 𝑛 is immediately given by ℙ(𝑇⩽𝑛)=ℙ(𝐗_𝑛=𝟎)=𝐟_𝑛(𝟎)^𝐗0, that is to say, ℙ(𝑇⩽𝑛)=𝑓_𝑛,1(𝟎)^𝑋0𝑓_𝑛,2(𝟎)^𝑋−1⋯𝑓_𝑛,𝑑(𝟎)^{𝑋−𝑑+1}. (2.7) It can be immediately deduced from convergence results for 𝐟_𝑛(𝟎) as 𝑛→∞ [8], that if 𝜌=1, ℙ(𝑇⩽𝑛)∼1−(𝑛𝜂)⁻¹𝐗₀⋅𝐮, while if 𝜌<1, ℙ(𝑇⩽𝑛)∼1−𝜌^𝑛𝛾𝐗₀⋅𝐮, 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 ℙ(𝑇⩽𝑛) for any 𝑛 by the iterative computation of 𝐟_𝑛, 𝐗₀ being given, when the parameters Ψ_𝑘 of (1.1) are known (resp., estimated). Moreover, since for 𝜌⩽1 the epidemic becomes extinct in an a.s. finite time and ℙ(𝑇⩽𝑛)↗_𝑛→∞1, then for any given probability 𝑝∈[0,1[ there exists 𝑛∈ℕ such that ℙ(𝑇⩽𝑛)⩾𝑝. So in practice, for any 𝑝∈[0,1[, (2.7) enables us to compute the 𝑝-quantile 𝑛^𝑇_𝑝 of the extinction time, 𝑛^𝑇_𝑝∶=min𝑛⩾1∶ℙ(𝑇⩽𝑛)⩾𝑝. (2.8)

2.2. Total Size of the Epidemic

Under the assumption 𝜌⩽1 and the independence of the {𝑌_{𝑛−𝑙,𝑛,𝑖}}_{𝑖,𝑙,𝑛} (we previously assumed the independence of the {𝑌_{𝑛−𝑙,𝑛,𝑖}}_𝑖,𝑙, for each 𝑛), we derive the distribution of the future total size 𝑁∶=∑^𝑇_𝑛=1𝑋_𝑛 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 𝐗₀, the time origin being the same as in Section 2.1, the probability distribution of 𝑁 is 𝑁^𝒟=^𝑑𝑋_𝑘=1^−𝑘+1_𝑖=1^{𝑌𝑘,𝑖}_𝑗=1𝑁_{𝑘,𝑖,𝑗},𝑌_𝑘,𝑖∶=^𝑑_𝑙=𝑘𝑌_{−k+1,−𝑘+1+𝑙,𝑖}, (2.9) where ^𝒟= denotes the equality in distribution, an empty sum is by convention 0, the {𝑌_𝑘,𝑖}_𝑘,𝑖 are independent, and the {𝑌_𝑘,𝑖}_𝑖 and the {𝑁_{𝑘,𝑖,𝑗}}_{𝑘,𝑖,𝑗} are i.i.d. with 𝑌^𝒟_𝑘,𝑖=Poiss^𝑑_𝑙=𝑘Ψ_𝑙,𝑁^𝒟_{𝑘,𝑖,𝑗}=Borel−Tanner^𝑑_𝑙=1Ψ_𝑙,1, (2.10) that is, for each 𝑛⩾1, ℙ(𝑁_{𝑘,𝑖,𝑗}=𝑛)=𝑒^{−𝑛∑𝑑𝑙=1Ψ𝑙}(𝑛∑^𝑑_𝑙=1Ψ_𝑙)^𝑛−1(𝑛!)⁻¹. Consequently, under the convention that an empty product is 1, the probability distribution of 𝑁 is, for any 𝑛∈ℕ, ℙ(𝑁=𝑛)=_{𝑗}∶{0⩽𝑦𝑘,𝑖⩽𝑛,{1⩽𝑛𝑘,𝑖,𝑗⩽𝑛}𝑖,𝑘∑𝑑𝑘=1∑𝑋−𝑘+1𝑖=1∑𝑦𝑘,𝑖𝑗=1𝑛𝑘,𝑖,𝑗}^𝑑_=𝑛^𝑋_𝑘=1^−𝑘+1_𝑖=1𝑒^{−∑𝑑𝑙=𝑘Ψ𝑙}∑^𝑑_𝑙=𝑘Ψ_𝑙^{𝑦𝑘,𝑖}𝑦_𝑘,𝑖!×^{𝑦𝑘,𝑖}_𝑗=1𝑒^{−𝑛𝑘,𝑖,𝑗∑𝑑𝑙=1Ψ𝑙}𝑛_{𝑘,𝑖,𝑗}∑^𝑑_𝑙=1Ψ_𝑙^{𝑛𝑘,𝑖,𝑗−1}𝑛_{𝑘,𝑖,𝑗}!, (2.11) which may be calculated (resp., estimated), replacing the Ψ_𝑘 by their values (resp., estimations). In practice, for any 𝑝∈[0,1[, (2.11) enables to compute the 𝑝-quantile 𝑛^𝑁_𝑝 of the total epidemic size, 𝑛^𝑁_𝑝∶=min𝑛⩾1∶ℙ(𝑁⩽𝑛)⩾𝑝. (2.12) We obtain moreover an explicit formula for the mean value and variance of the size of the epidemic, 𝔼(𝑁)=∑^𝑑_𝑘=1𝑋_−𝑘+1∑^𝑑_𝑙=𝑘Ψ_𝑙1−∑^𝑑_𝑙=1Ψ_𝑙(𝑁)=∑,𝕍ar^𝑑𝑘=1𝑋_−𝑘+1∑^𝑑_𝑙=𝑘Ψ_𝑙1−∑^𝑑_𝑙=1Ψ_𝑙³. (2.13)

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 {𝑍_𝑛}_𝑛 defined by the conditional distribution 𝑍_𝑛∣𝑋_𝑛^𝒟=Poiss(Ψ₀𝑋_𝑛), where Ψ₀ is the mean number of individuals infected at time 𝑛 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 {𝑋_𝑛}_𝑛.

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 𝑝-quantile 𝑛^𝑇_𝑝 of the extinction time with the probability 𝑝 as large as wanted (see (2.8)), the epidemic might become extinct after this time with a small but nonnull probability of order 1−𝑝. 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 {𝐗_𝑛}_𝑛 conditionally on the event that the epidemic has not become extinct at time 𝑘, where 𝑘 is very large. We therefore consider for any 𝑛₁,𝑛₂,…∈ℕ and any 𝐢₀,𝐢₁,𝐢₂,…∈ℕ^𝑑 the conditioned probability ℙ_𝐢0(𝐗_𝑛1=𝐢₁,…,𝐗_𝑛𝑟=𝐢_𝑟∣𝐗_𝑘≠𝟎), where the subscript 𝐢₀ denotes the initial value. If 𝑘 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 𝑘→∞, to the distribution of a 𝑑-dimensional Markov process {𝐗^∗_𝑛}_𝑛: lim_𝑘→∞ℙ_𝐢0𝐗_𝑛1=𝐢₁,…,𝐗_𝑛𝑟=𝐢_𝑟∣𝐗_𝑘≠𝟎=ℙ_𝐢0𝐗^∗_𝑛1=𝐢₁,…,𝐗^∗_𝑛𝑟=𝐢_𝑟. (2.14) We will further discuss in Proposition 2.5 the relevancy of approximating the conditioned probability for 𝑘 fixed by the limiting object (2.14). The conditioned process {𝐗^∗_𝑛}_𝑛 defined by (2.14) is known in the literature as the 𝑄-process associated with {𝐗_𝑛}_𝑛, also described as the process conditioned on “not being extinct in the distant future.” It has the following transition probability [10]: for every 𝑛⩾1, 𝐢,𝐣∈ℕ^𝑑, 𝐢≠𝟎, ℙ𝐗^∗_𝑛=𝐣∣𝐗^∗_𝑛−1=1=𝐢𝜌𝐣⋅𝐮ℙ𝐗𝐢⋅𝐮^𝑛=𝐣∣𝐗_𝑛−1=𝐢, (2.15) where 𝐮 is the normalized right eigenvector of 𝐌 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 {𝐗_𝑛}_𝑛, we define the 1-dimensional process 𝑋^∗_𝑛∶=𝑋^∗_𝑛,1. By construction we then have 𝑋^∗_𝑛,𝑖=𝑋^∗_{𝑛−𝑖+1}, for each 𝑛 and each 𝑖=1⋯𝑑.

Proposition 2.2. The stochastic process {𝑋^∗_𝑛}_𝑛 is, conditionally on its past, distributed as the sum of two independent Poisson and Bernoulli random variables: 𝑋^∗_𝑛∣𝐗^∗𝒟_𝑛−1=Poiss𝐗^∗_𝑛−1⋅Ψ∗ℬ𝑝𝐗^∗_𝑛−1, (2.16) where Ψ∶=(Ψ₁,…,Ψ_𝑑), ∗ is the convolution product symbol, and 𝑝𝐗^∗_𝑛−1𝑢∶=¹𝐗^∗_𝑛−1⋅Ψ𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘. (2.17)

Proof. Applying (1.1) and (2.15), we obtain that for all 𝑗∈ℕ, ℙ𝑋^∗_𝑛=𝑗∣𝐗^∗_𝑛−1=ℙ𝐗^∗_𝑛=𝑗,𝑋^∗_𝑛−1,…,𝑋^∗_{𝑛−(𝑑−1)}∣𝐗^∗_𝑛−1=𝑗𝑢₁+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝜌𝐗^∗_𝑛−1ℙ𝐗⋅𝐮^𝑛=𝑗,𝑋^∗_𝑛−1,…,𝑋^∗_{𝑛−(𝑑−1)}∣𝐗_𝑛−1=𝐗^∗_𝑛−1=𝑗𝑢₁+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝜌𝐗^∗_𝑛−1𝐗⋅𝐮^∗𝑛−1⋅Ψ^𝑗𝑒𝑗!^{−𝐗∗𝑛−1⋅Ψ}=𝑢₁𝐗^∗_𝑛−1⋅Ψ𝜌𝐗^∗_𝑛−1𝐗⋅𝐮^∗𝑛−1⋅Ψ^𝑗−1!𝑒𝑗−1^{−𝐗∗𝑛−1⋅Ψ}+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝜌𝐗^∗_𝑛−1𝐗⋅𝐮^∗𝑛−1⋅Ψ^𝑗𝑒𝑗!^{−𝐗∗𝑛−1⋅Ψ}. (2.18) The equality 𝐌𝐮^𝑇=𝜌𝐮^𝑇 implies that for all 𝑘=1⋯𝑑−1, 𝜌𝑢_𝑘=Ψ_𝑘𝑢₁+𝑢_𝑘+1, and that 𝜌𝑢_𝑑=Ψ_𝑑𝑢₁. Consequently, 𝜌𝐗^∗_𝑛−1⋅𝐮=𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘, and thus ℙ𝑋^∗_𝑛=𝑗∣𝐗^∗_𝑛−1=𝑢₁𝐗^∗_𝑛−1⋅Ψ𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝐗^∗_𝑛−1⋅Ψ^𝑗−1!𝑒𝑗−1^{−𝐗∗𝑛−1⋅Ψ}+𝑢1−¹𝐗^∗_𝑛−1⋅Ψ𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝐗^∗_𝑛−1⋅Ψ^𝑗𝑒𝑗!^{−𝐗∗𝑛−1⋅Ψ}=Poiss𝐗^∗_𝑛−1𝑢⋅Ψ∗ℬ¹𝐗^∗_𝑛−1⋅Ψ𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝑗. (2.19)

Remark 2.3. Note that if one compares (2.16) with the transition probability of the unconditioned process 𝑋_𝑛∣𝐗^𝒟_𝑛−1=Poiss(𝐗_𝑛−1⋅Ψ), it appears that {𝑋^∗_𝑛}_𝑛 behaves at each time step like {𝑋_𝑛}_𝑛, 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 𝑋^∗_𝑛−1=⋯=𝑋^∗_{𝑛−(𝑑−1)}=0, then according to (2.17), 𝑝(𝐗^∗_𝑛−1)=1, which implies that at time 𝑛, the probability to add one unit is equal to one, thus preventing the extinction of the process.

Proposition 2.4. The process {𝐗^∗_𝑛}_𝑛 admits a stationary probability measure 𝜋 with finite first- and second-order moments.

Proof. Since the multitype branching process {𝐗_𝑛}_𝑛 satisfies property (2.4), it is known [10] that the 𝑄-process {𝐗^∗_𝑛}_𝑛 is positive recurrent with a stationary probability measure 𝜋 given by 𝜋(𝐢)(𝐢)∶=𝐢⋅𝐮𝜈∑_𝐤∈ℕ^𝑑𝐤⋅𝐮𝜈(𝐤),𝐢∈ℕ^𝑑, (2.20) where 𝜈 is the Yaglom distribution of the process {𝐗_𝑛}_𝑛, uniquely defined by the following property: for all 𝐢,𝐣∈ℕ^𝑑⧵{𝟎}, lim_𝑛→∞ℙ(𝐗_𝑛=𝐢∣𝐗₀=𝐣,𝐗_𝑛≠𝟎)=𝜈(𝐢). In the literature, this stationary measure for the conditioned process {𝐗^∗_𝑛}_𝑛 is also referred to as the doubly limiting conditional probability. Moreover, by Proposition 2.2, 𝔼𝑋^∗_𝑛=𝔼𝔼𝑋^∗_𝑛∣𝐗^∗_𝑛−1=𝔼𝐗^∗_𝑛−1⋅Ψ+𝑝𝐗^∗_𝑛−1⩽^𝑑_𝑘=1𝔼𝑋^∗_𝑛−𝑘Ψ_𝑘+1, (2.21) which implies that lim_𝑛→∞𝔼(𝑋^∗_𝑛)⩽(1−∑^𝑑_𝑘=1Ψ_𝑘)⁻¹<∞. We consequently obtain by means of Fatou's lemma that, for every 𝑖=1⋯𝑑, _{𝐣∈ℕ𝑑}𝑗_𝑖𝜋(𝐣)=𝔼lim_𝑛→∞𝑋^∗_𝑛,𝑖=𝔼lim_𝑛→∞𝑋^∗_{𝑛−𝑖+1}=𝔼lim_𝑛→∞𝑋^∗_𝑛⩽lim_𝑛→∞𝔼𝑋^∗_𝑛<∞. (2.22) We similarly prove that 𝜋 has finite second-order moments by writing 𝕍ar𝑋^∗_𝑛∣𝐗^∗_𝑛−1=𝐗^∗_𝑛−1𝑢⋅Ψ+¹𝐗^∗_𝑛−1⋅Ψ∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘𝑢₁𝐗^∗_𝑛−1⋅Ψ+∑^𝑑_𝑘=2𝑋^∗_{𝑛−𝑘+1}𝑢_𝑘²⩽𝐗^∗_𝑛−11⋅Ψ+4. (2.23)

Let us discuss the relevancy of approximating the epidemic process {𝐗_𝑛}_𝑛 conditioned on nonextinction at some finite time 𝑘, for 𝑘 large, by the 𝑄-process {𝐗^∗_𝑛}_𝑛 obtained by letting 𝑘→∞. 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 𝑘 for the survival of the disease in the population. We therefore might consider that 𝑘 is chosen large enough such that the approximation of the process {𝐗_𝑛}_𝑛 conditioned on the event {𝐗_𝑘≠𝟎} by the process {𝐗^∗_𝑛}_𝑛 is valid. Of course, the order of magnitude of such 𝑘 will depend on the rate of convergence of the conditioned process to {𝐗^∗_𝑛}_𝑛.

Proposition 2.5. Let 𝑛₁⩽⋯⩽𝑛_𝑟⩽𝑘 and 𝐢₀,…,𝐢_𝑟∈ℕ^𝑑⧵{𝟎}. Then the difference |ℙ_𝐢0(𝐗_𝑛1=𝐢₁,…,𝐗_𝑛𝑟=𝐢_𝑟∣𝐗_𝑘≠𝟎)−ℙ_𝐢0(𝐗^∗_𝑛1=𝐢₁,…,𝐗^∗_𝑛𝑟=𝐢_𝑟)| decreases, as 𝑘→∞, with max{𝑘^𝑠−1(|𝜆|𝜌⁻¹)^𝑘/2,𝜌^𝑘/2}, where 𝜆 is an eigenvalue of 𝐌 such that 𝜌>|𝜆|⩾|𝜆₃|⩾|𝜆₄|⩾⋯, with the 𝜆_𝑖 being the other eigenvalues of 𝐌. In case |𝜆|=|𝜆₃| we stipulate that the multiplicity 𝑠 of 𝜆 is at least as great as the multiplicity of 𝜆₃.

Proof. Thanks to (2.15) and to the Markov property together with the fact that ℙ_𝐢(𝐗_𝑛=𝟎)=𝐟_𝑛(𝟎)^𝐢, we have ℙ_𝐢0𝐗_𝑛1=𝐢₁,…,𝐗_𝑛𝑟=𝐢_𝑟∣𝐗_𝑘≠𝟎−ℙ_𝐢0𝐗^∗_𝑛1=𝐢₁,…,𝐗^∗_𝑛𝑟=𝐢_𝑟=1−𝐟_{𝑘−𝑛𝑟}(𝟎)^𝐢𝑟1−𝐟_𝑘(𝟎)^𝐢0−1𝜌^𝑛𝑟𝐢_𝑟⋅𝐮𝐢₀ℙ⋅𝐮^𝐢0𝐗_𝑛1=𝐢₁,…,𝐗_𝑛𝑟=𝐢_𝑟. (2.24) The right term of (2.24) is known to converge to 0, as 𝑘→∞, thanks to the property that lim_𝑘𝐚_𝑘=𝛾𝐮, for some 𝛾>0, where 𝐚_𝑘∶=𝜌^−𝑘(𝟏−𝐟_𝑘(𝟎)) (see [8]). This stems from two convergences, namely, lim_𝑘𝑏_𝑘=𝛾, where 𝑏_𝑘∶=𝜌^−𝑘(𝐯⋅(𝟏−𝐟_𝑘(𝟎))), and lim_𝑘𝐚_𝑘𝑏_𝑘⁻¹=𝐮. Let us write 𝐚_𝑘=𝛾𝐮+𝜀_𝑘, where lim_𝑘𝜀_𝑘=𝟎. Since 𝑏_𝑘=𝐯⋅𝐚_𝑘 and 𝐮⋅𝐯=1, it comes 𝐚_𝑘𝑏_𝑘⁻¹=𝛾𝐮+𝜀_𝑘𝛾+𝐯⋅𝜀_𝑘∼_𝑘→∞𝐮+𝛾⁻¹𝜀_𝑘−𝐯⋅𝜀_𝑘𝐮. (2.25) It thus appears that the rate of convergence of 𝐚_𝑘 to 𝛾𝐮 is of the same order of magnitude as the one of 𝐚_𝑘𝑏_𝑘⁻¹ to 𝐮. Let us determine this rate in an accurate way. We use the following inequality produced by Joffe and Spitzer in [8]: for each 𝑘⩾𝑛⩾1, 𝐚_𝑘𝑏_𝑘⁻¹−𝐮=𝟏−𝐟_𝑘(𝟎)𝐯⋅𝟏−𝐟_𝑘(𝟎⩽)−𝐮2𝛿_𝑛+∑^𝑘_{𝑗=𝑘−𝑛+1}𝛼_𝑗1−𝛿_𝑛−∑^𝑘_{𝑗=𝑘−𝑛+1}𝛼_𝑗, (2.26) where we are going to replace 𝛿_𝑛 and 𝛼_𝑛 by some explicit formulae function of 𝜌 and 𝑛. For this purpose, we use a detailed asymptotic behavior of 𝐌^𝑘, as 𝑘→∞, presented for instance in [11]; we have 𝐌^𝑘=𝜌^𝑘𝐑+𝑂(𝑘^𝑠−1|𝜆|^𝑘), where 𝐑=𝐮^𝑇𝐯. For the sake of clarity the symbol 𝑂(⋅) will denote either a scalar or a matrix with all the entries satisfying the associated property. This implies the existence of some constant 𝑎>0 such that, for all 𝑘∈ℕ, (1−𝛿_𝑘)𝐑⩽𝜌^−𝑘𝐌^𝑘⩽(1+𝛿_𝑘)𝐑, where 𝛿_𝑘∶=𝑎𝑘^𝑠−1(𝜌⁻¹|𝜆|)^𝑘. Moreover, following [8], let us write, for all 𝐫∈[0,1]^𝑑, 𝟏−𝐟(𝐫)=(𝐌−𝐄(𝐫))(𝟏−𝐫), where 𝟎⩽𝐄(𝐫)⩽𝐌, and 𝐄(𝐫)=𝑂(‖𝟏−𝐫‖) as 𝐫→𝟏. Then 𝜌⁻¹𝐄(𝐟_𝑘−1(𝟎))=𝜌⁻¹𝑂(‖𝟏−𝐟_𝑘−1(𝟎)‖)=𝑂(𝜌^𝑘−2), which implies the existence of some constant 𝑏>0 such that, for all 𝑘∈ℕ, 𝟎⩽𝜌⁻¹𝐄(𝐟_𝑘−1(𝟎))⩽𝛼_𝑘𝐑, with 𝛼_𝑘∶=𝑏𝜌^𝑘−2. We thus have provided an explicit formula for the sequences (𝛿_𝑘)_𝑘 and (𝛼_𝑘)_𝑘 introduced by Joffe and Spitzer in [8]. Finally let us apply (2.26) to 𝑛=⌊𝑘/2⌋ and replace in this inequality 𝛿_𝑛 and 𝛼_𝑛 by their explicit expressions that we got. We obtain that, for all 𝑘∈ℕ, 𝟏−𝐟_𝑘(𝟎)𝐯⋅𝟏−𝐟_𝑘(𝟎)⩽2−𝐮^2−𝑠𝑎𝑘^𝑠−1|𝜆|/𝜌^𝑘/2+𝑏/𝜌1−𝜌𝜌^𝑘/21−2^1−𝑠𝑎𝑘^𝑠−1|𝜆|/𝜌^𝑘/2−𝑏/𝜌1−𝜌𝜌^𝑘/2. (2.27) Consequently, the right member of (2.24) will decrease with max{𝑘^𝑠−1(|𝜆|𝜌⁻¹)^𝑘/2,𝑏(𝜌(1−𝜌))⁻¹𝜌^𝑘/2}, as 𝑘→∞.

Hence the concept of the 𝑄-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 Ψ_𝑘 of the epidemic model (1.1)–(2.1) are not entirely known. More precisely, we assume that the Ψ_𝑘 are of the form, for all 𝑘=1⋯𝑑, Ψ_𝑘𝜃₀=𝑎_𝑘𝜃₀+𝑏_𝑘, (2.28) where 𝑎_𝑘>0 and 𝑏_𝑘⩾0 are constants, and 𝜃₀ is an unknown real parameter. We will write in what follows Ψ(𝜃₀)=𝐚𝜃₀+𝐛, where Ψ(𝜃₀)∶=(Ψ₁(𝜃₀),…,Ψ_𝑑(𝜃₀)) and so forth. This general assumption corresponds in particular to the case where the Ψ_𝑘 are of the form (1.2).

We estimate 𝜃₀ by the following WCLSE in model (1.1)–(2.1). This estimator generalizes the well-known Harris estimator [12] in a BGW process. Let Θ∶=]𝜃₁,𝜃₂[, 𝜃₂>𝜃₁>0, such that 𝜃₀∈Θ. The WCLSE is based on the normalized process 𝑌_𝑛∶=𝑋_𝑛/𝐚⋅𝐗_𝑛−1 and is defined by 𝜃_∣𝐗0∣∶=argmin^𝑛_{𝜃∈Θ𝑘=1}𝑌_𝑘−𝔼_𝜃𝑌_𝑘∣𝐗_𝑘−1²=argmin^𝑛_{𝜃∈Θ𝑘=1}𝑋_𝑘−Ψ(𝜃)⋅𝐗_𝑘−1²𝐚⋅𝐗_𝑘−1. (2.29) We easily derive the following explicit form: 𝜃_∣𝐗0∣=∑^𝑛_𝑘=1𝑋_𝑘−𝐛⋅𝐗_𝑘−1∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1. (2.30) The normalization of the process 𝑋_𝑛 by 𝐚⋅𝐗_𝑛−1 appears to be the most natural and suitable for the following reasons. First, this normalization generalizes the normalization 𝑋_𝑛/𝑎𝑋_𝑛−1 in the monotype case, which is the one leading to the Harris estimator 𝑚 of 𝑚₀=𝑎𝜃₀+𝑏 since we have, for 𝑑=1, 𝑎𝜃_𝑋0+𝑏=𝑚. It also corresponds, in the linear case 𝐛=𝟎, to the maximum likelihood estimator of 𝜃₀. In addition, defining for any vector 𝐱, 𝐱∶=min_𝑖𝑥_𝑖 and 𝐱∶=max_𝑖𝑥_𝑖, we have 𝜃₀+𝐛𝐚⩽𝔼_𝜃0𝑌_𝑘−𝔼_𝜃0𝑌_𝑘∣𝐗_𝑘−1²∣𝐗_𝑘−1=𝜃₀+𝐛⋅𝐗_𝑘−1𝐚⋅𝐗_𝑘−1⩽𝜃₀+𝐛𝐚, (2.31) hence the conditional variance of the error term 𝑌_𝑘−𝔼_𝜃0(𝑌_𝑘∣𝐗_𝑘−1) in the stochastic regression equation 𝑌_𝑘=𝔼_𝜃0(𝑌_𝑘∣𝐗_𝑘−1)+𝑌_𝑘−𝔼_𝜃0(𝑌_𝑘∣𝐗_𝑘−1) is invariant under multiplication of the whole process, and bounded respectively to {𝐗_𝑛}_𝑛, leading to the quasi-optimality of 𝜃_∣𝐗0∣ at finite |𝐗₀| and 𝑛, in the sense of [13].

Let us provide asymptotic results for the estimator 𝜃_∣𝐗0∣ defined by (2.30), as the initial population size |𝐗₀|=𝑋₀+𝑋₋₁+⋯+𝑋_−𝑑+1 tends to infinity. We denote by 𝑚^(𝑘)_𝑖𝑗(𝜃) the (𝑖,𝑗)th entry in the 𝑘th power of the matrix 𝐌(𝜃) given by (2.3).

Theorem 2.6. Let us assume that, for each 𝑖=1⋯𝑑, there exists some 𝛼_𝑖∈[0,1] such that lim_{∣𝐗0∣→∞}𝑋_0,𝑖|𝐗₀|⁻¹=𝛼_𝑖. Then 𝜃_∣𝐗0∣ is strongly consistent, that is, lim_{∣𝐗0∣→∞}𝜃=_∣𝐗0∣^𝑎.𝑠.𝜃₀, and is asymptotically normally distributed: lim_{∣𝐗0∣→∞}∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1𝜎²𝜃_∣𝐗0∣𝜃_∣𝐗0∣−𝜃₀^𝒟=𝒩(0,1), (2.32) where 𝜎²(𝜃)∑∶=𝜃+^𝑛𝑘=1∑^𝑑_𝑗=1∑^𝑑_𝑖=1𝛼_𝑗𝑏_𝑖𝑚^(𝑘−1)_𝑗𝑖(𝜃)∑^𝑛_𝑘=1∑^𝑑_𝑗=1∑^𝑑_𝑖=1𝛼_𝑗𝑏_𝑖𝑚^(𝑘−1)_𝑗𝑖(𝜃). (2.33)

Proof. Let us first prove that, for each 𝑘=1⋯𝑛 and each 𝑖=1⋯𝑑, lim_{∣𝐗0∣→∞}𝑋_𝑘,𝑖𝐗₀^a.s.=^𝑑_𝑗=1𝛼_𝑗𝑚^(𝑘)_𝑗𝑖𝜃₀. (2.34) Using the branching property of the process {𝐗_𝑛}_𝑛 ([6]) we write 𝑋_𝑘,𝑖=^𝑋0,1_𝑗=1𝑋⁽¹⁾_{𝑘,𝑖,𝑗}+⋯+^𝑋0,𝑑_𝑗=1𝑋^(𝑑)_{𝑘,𝑖,𝑗}, (2.35) where, for all 𝑙=1⋯𝑑 and 𝑗=1⋯𝑋_0,𝑙, 𝑋^(𝑙)_{𝑘,𝑖,𝑗} is the 𝑖th coordinate of a 𝑑-type branching process at time 𝑘 initialized by a single particle of type 𝑙. For 𝑘, 𝑖, and 𝑙 fixed the random variables {𝑋^(𝑙)_{𝑘,𝑖,𝑗}}_𝑗 are i.i.d. with mean value 𝑚^(𝑘)_𝑙𝑖(𝜃₀). According to the strong law of large numbers and under the theorem assumption, we have, for every 𝑙=1⋯𝑑 such that 𝑋_0,𝑙≠0, lim_{0||𝐗→∞}∑^𝑋0,𝑙_𝑗=1𝑋^(𝑙)_{𝑘,𝑖,𝑗}𝑋_0,𝑙^a.s.=𝑚^(𝑘)_𝑙𝑖𝜃₀, (2.36) which together with the theorem assumption leads to (2.34).
To prove the consistency of 𝜃_∣𝐗0∣ we apply (2.34) to (2.30), using the fact that 𝑋_𝑘=𝑋_𝑘,1 and 𝑋_𝑘−𝑖=𝑋_{𝑘−1,𝑖}, and obtain lim_{∣𝐗0∣→∞}𝜃=_∣𝐗0∣^a.s.∑^𝑛_𝑘=1∑^𝑑_𝑗=1𝛼_𝑗𝑚^(𝑘)_𝑗1𝜃₀−∑^𝑑_𝑖=1𝑏_𝑖𝑚^(𝑘−1)_𝑗𝑖𝜃₀∑^𝑛_𝑘=1∑^𝑑_𝑖=1∑^𝑑_𝑗=1𝑎_𝑖𝛼_𝑗𝑚^(𝑘−1)_𝑗𝑖𝜃₀. (2.37) By definition, 𝑚^(𝑘)_𝑗1𝜃₀=^𝑑_𝑖=1𝑚^(𝑘−1)_𝑗𝑖𝜃₀𝑚_𝑖1𝜃₀=^𝑑_𝑖=1𝑚^(𝑘−1)_𝑗𝑖𝜃₀𝑎_𝑖𝜃₀+𝑏_𝑖, (2.38) hence (2.37) immediately leads to the strong consistency.
We are now interested in the asymptotic distribution of 𝜃_∣𝐗0∣−𝜃₀. We derive from (2.30) that ^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1𝜃_∣𝐗0∣−𝜃₀=∑^𝑛_𝑘=1𝑋_𝑘−Ψ𝜃₀⋅𝐗_𝑘−1∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1. (2.39) By (1.1), 𝑋_𝑘−Ψ𝜃₀⋅𝐗_𝑘−1=^𝑑𝑋_𝑖=1^𝑘−𝑖_𝑗=1𝑌_{𝑘−𝑖,𝑘,𝑗}−Ψ_𝑖𝜃₀=∶^𝑑𝑋_𝑖=1^{𝑘−𝑖∘}_𝑗=1𝑌_{𝑘−𝑖,𝑘,𝑗}, (2.40) where the {𝑌_{𝑘−𝑖,𝑘,𝑗}}_𝑗 are i.i.d. given ℱ_𝑘−1, following a Poisson distribution with parameter Ψ_𝑖(𝜃₀), and the {𝑌_{𝑘−𝑖,𝑘,𝑗}}_𝑖,𝑗 are independent given ℱ_𝑘−1. Renumbering the 𝑌^∘_{𝑘−𝑖,𝑘,𝑗} we then obtain ^𝑛_𝑘=1𝑋_𝑘−Ψ𝜃₀⋅𝐗_𝑘−1=^𝑑∑_𝑖=1^{𝑛𝑘=1𝑋𝑘−𝑖∘}_𝑗=1𝑌_{𝑘−𝑖,𝑘,𝑗}. (2.41) 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 𝑖=1⋯𝑑, lim_{∣𝐗0∣→∞}∑^{∑𝑛𝑘=1𝑋𝑘−𝑖∘}_𝑗=1𝑌_{𝑘−𝑖,𝑘,𝑗}∑^𝑛_𝑘=1𝑋^𝒟_𝑘−𝑖=𝒩0,𝑎_𝑖𝜃₀+𝑏_𝑖. (2.42) We have used the fact that |𝐗₀| is a real positive sequence growing to infinity, and ∑^𝑛_𝑘=1𝑋_𝑘−𝑖 is a sequence of integer-valued random variables such that ∑^𝑛_𝑘=1𝑋_𝑘−𝑖/|𝐗₀| converges in probability to a finite random variable. In our case the limit is actually deterministic, since we have shown in (2.34) that lim_𝐗0→∞∑^𝑛_𝑘=1𝑋_𝑘−𝑖𝐗₀^a.s.=^𝑛𝑑_𝑘=1𝑗=1𝛼_𝑗𝑚^(𝑘−1)_𝑗𝑖𝜃₀. (2.43) Using (2.41) in (2.39), we write ^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1𝜃_∣𝐗0∣−𝜃₀=^𝑑_𝑖=1∑^{∑𝑛𝑘=1𝑋𝑘−𝑖∘}_𝑗=1𝜁_{𝑘−𝑖,𝑘,𝑗}∑^𝑛_𝑘=1𝑋_𝑘−𝑖∑^𝑛_𝑘=1𝑋_𝑘−𝑖∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1. (2.44) Using again (2.34), lim_{∣𝐗0∣→∞}∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−𝑖∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1^a.s.=∑^𝑛_𝑘=1∑^𝑑_𝑗=1𝛼_𝑗𝑚^(𝑘−1)_𝑗𝑖𝜃₀∑^𝑛_𝑘=1∑^𝑑_𝑗=1∑^𝑑_𝑙=1𝛼_𝑗𝑎_𝑙𝑚^(𝑘−1)_𝑗𝑙𝜃₀, (2.45) which, combined to (2.42) and (2.44), implies by Slutsky's theorem that lim_{∣𝐗0∣→∞}^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1𝜃_∣𝐗0∣−𝜃₀^𝒟=𝒩0,𝜎²𝜃₀. (2.46) By (2.33) and the strong consistency, lim_∣𝐗0∣𝜎²(𝜃₀)/𝜎²(𝜃_∣𝐗0∣)=^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 {𝑋^∗_𝑛}_𝑛 introduced in Section 2.4, we need to estimate the parameter 𝜃₀ in the setting of this conditioned process. We point out that 𝜃₀ does not play the same role in the conditioned process {𝑋^∗_𝑛}_𝑛 and in the unconditioned process {𝑋_𝑛}_𝑛, 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 𝜃₀ with an estimator aimed for the unconditioned process, such as 𝜃_∣𝐗0∣. Let us notice that, according to (2.16), the process {𝐗^∗_𝑛}_𝑛 could be written as a multitype branching process with state dependent immigration. Because of this state-dependence, and since the parameter 𝜃₀ 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 𝑌^∗_𝑛∶=𝑋^∗_𝑛/𝐚⋅𝐗^∗_𝑛−1, namely, 𝜃^∗_𝑛∶=argmin_𝜃∈Θ𝑆_𝑛(𝜃),𝑆_𝑛(𝜃)∶=^𝑛_𝑘=1𝑌^∗_𝑘−𝑓𝜃,𝐗^∗_𝑘−1², (2.47) where Θ is defined in Section 2.5, and where 𝑓𝜃₀,𝐗^∗_𝑘−1∶=𝔼_𝜃0𝑌^∗_𝑘∣𝐗^∗_𝑘−1=𝐗^∗_𝑘−1⋅Ψ𝜃₀+𝑝𝜃₀,𝐗^∗_𝑘−1𝐚⋅𝐗^∗_𝑘−1. (2.48) Let 𝜀^∗_𝑘∶=𝑌^∗_𝑘−𝑓(𝜃₀,𝐗^∗_𝑘−1) be the error term between the normalized process and its conditional expectation. We obtain that 𝑔𝜃₀,𝐗^∗_𝑘−1∶=𝔼_𝜃0𝜀^∗_𝑘²∣𝐗^∗_𝑘−1=𝐗^∗_𝑘−1⋅Ψ𝜃₀+𝑝𝜃₀,𝐗^∗_𝑘−11−𝑝𝜃₀,𝐗^∗_𝑘−1𝐚⋅𝐗^∗_𝑘−1, (2.49) which implies 𝜃₀+𝐛𝐚⩽𝔼_𝜃0𝜀^∗_𝑘²∣𝐗^∗_𝑘−1⩽𝜃₀+𝐛+1𝐚. (2.50) In what follows, we denote by 𝑓^′ the derivative of 𝑓 with respect to 𝜃, and similarly for the other quantities depending on 𝜃.

Theorem 2.7. The estimator 𝜃^∗_𝑛 is strongly consistent, that is, lim_𝑛→∞𝜃=^∗_𝑛^𝑎.𝑠.𝜃₀, and has the following asymptotic distribution: lim_𝑛→∞∑^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²∑^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗^∗_𝑘𝜃^∗_𝑛−𝜃₀^𝒟=𝒩(0,1). (2.51)

Remark 2.8. Note that (2.51) involves the function 𝑓^′ and thus requires the knowledge of the derivative of the function 𝑢_𝑗 given in (2.5), which is not an explicit function of 𝜃 since 𝜌 is not either. However, 𝜌^′ satisfies 𝜌^′=∑^𝑑_𝑘=1𝑎_𝑘𝜌^−𝑘[∑^𝑑_𝑘=1𝑘(𝑎_𝑘𝜃+𝑏_𝑘)𝜌^−𝑘−1]⁻¹; hence 𝑢^′_𝑗 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 𝜌(𝜃^∗obs_𝑛) is known; for this purpose one can for instance numerically approximate the largest solution 𝜌 of ∑^𝑑_𝑘=1Ψ_𝑘(𝜃^∗obs_𝑛)𝜌^−𝑘=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 {𝐗^∗_𝑛}_𝑛 and its stationary distribution 𝜋_𝜃0 (given by Proposition 2.4), which states [16]: for every 𝜋_𝜃0-integrable function ℎ∶ℕ^𝑑⧵{𝟎}→ℝ, lim_𝑛→∞1𝑛^𝑛−1_𝑘=0ℎ𝐗=^∗_𝑘^a.s._{𝐣∈ℕ𝑑}ℎ(𝐣)𝜋_𝜃0(𝐣). (2.52) Note that the Perron's root 𝜌(𝜃) of 𝐌(𝜃) and the associated right normalized eigenvector 𝐮(𝜃) are 𝒞^∞-functions of 𝜃.
In order to prove the strong consistency of 𝜃^∗_𝑛 as 𝑛→∞, since 𝑓^′′(⋅,𝐗^∗_𝑘−1) is not linear, we cannot use the general standard method based on the first-order expansion of 𝑆^′_𝑛(𝜃^∗_𝑛) at 𝜃₀, and on the strong law of large numbers for martingales applied to 𝑆^′_𝑛(𝜃₀) correctly normalized. We consequently use the following conditions given in [17]:(i)𝑓(.,𝐗^∗_𝑘−1) is Lipschitz on Θ; that is, there exists a nonnegative ℱ^∗_𝑘−1-measurable function 𝐴_𝑘 (where ℱ^∗_𝑘−1∶=𝜎(𝐗^∗₀,…,𝐗^∗_𝑘−1)), satisfying, for all 𝛿₁, 𝛿₂∈Θ, |𝑓(𝛿₁,𝐗^∗_𝑘−1)−𝑓(𝛿₂,𝐗^∗_𝑘−1)|⩽𝐴_𝑘|𝛿₁−𝛿₂| a.s., (ii)lim_𝑘→∞𝔼_𝜃0((𝜀^∗_𝑘)²∣𝐗^∗_𝑘−1)<^a.s.∞,(iii)lim_𝑛→∞inf_{𝜃∈Θ∣𝜃−𝜃0∣⩾𝛿}∑^𝑛_𝑘=1(𝑓(𝜃₀,𝐗^∗_𝑘−1)−𝑓(𝜃,𝐗^∗_𝑘−1))=^2a.s.∞. Let us note that, in the frame of a general model 𝑓, since Θ is compact, (i) is satisfied as soon as 𝑓(𝜃,𝐗^∗_𝑘−1) has a first derivative in 𝜃 with sup_𝜃∈Θ𝑓^′(𝜃,𝐗^∗_𝑘−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 𝐣∈ℕ^𝑑, 𝐣≠𝟎, 𝑓^′(𝜃,𝐣)=(𝐚⋅𝐣+𝑝^′(𝜃,𝐣))(𝐚⋅𝐣)^−1/2, where 𝑝^′ denotes the derivative of 𝑝 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 ∣𝜃−𝜃₀∣⩾𝛿, applying the mean value theorem to the 𝒞¹-function 𝑝(⋅,𝐗^∗_𝑘−1), ^𝑛_𝑘=1𝑓(𝜃₀,𝐗^∗_𝑘−1)−𝑓(𝜃,𝐗^∗_𝑘−1)²=𝜃₀−𝜃^2𝑛_𝑘=1𝐚⋅𝐗^∗_𝑘−1𝑝𝜃1+⁰,𝐗^∗_𝑘−1−𝑝𝜃,𝐗^∗_𝑘−1𝜃₀−𝜃𝐚⋅𝐗^∗_𝑘−1²⩾𝛿^2𝑛_𝑘=1𝐚⋅𝐗^∗_𝑘−1inf_𝜃∈Θ𝑝1+^′(𝜃,𝐗^∗_𝑘−1)𝐚⋅𝐗^∗_𝑘−1². (2.53) Let us show that the function 𝐣↦𝐚⋅𝐣inf_𝜃∈Θ(1+𝑝^′(𝜃,𝐣)/𝐚⋅𝐣)² is 𝜋_𝜃0-integrable. For every 𝜃∈Θ, 𝐣∈ℕ^𝑑 and 𝐣≠𝟎, denoting Ξ∶=sup_{𝑖,𝜃∈Θ}{𝑢_𝑖(𝜃),∣𝑢^′_𝑖(𝜃)∣}<∞ (by continuity of 𝑢^′_𝑖(𝜃) on Θ_𝛾⊃Θ), and 𝐮∶=min_{𝑖,𝜃∈Θ}𝑢_𝑖(𝜃)>0, 𝑝^{′(𝜃,𝐣)}=𝑢^′_𝑖(𝜃)𝐣⋅Ψ(𝜃)+𝑢₁(𝜃)𝐣⋅𝐚∑^𝑑_𝑖=2𝑗_𝑖𝑢_𝑖(𝜃)−𝑢₁(𝜃)𝐣⋅Ψ(𝜃)∑^𝑑_𝑖=2𝑗_𝑖𝑢^′_𝑖(𝜃)𝑢₁(𝜃)𝐣⋅Ψ(𝜃)+∑^𝑑_𝑖=2𝑗_𝑖𝑢_𝑖(𝜃)²⩽Ξ²𝐮²3𝐣⋅Ψ(𝜃)∑^𝑑_𝑖=2𝑗_𝑖𝐣⋅Ψ(𝜃)+∑^𝑑_𝑖=2𝑗_𝑖²⩽3Ξ²4𝐮²=∶𝐶₁. (2.54) Hence, for all 𝐣≠𝟎, ∣𝐚⋅𝐣inf_𝜃∈Θ(1+𝑝^′(𝜃,𝐣)/𝐚⋅𝐣)²∣⩽(1+𝐶₁/𝐚)²𝐚⋅𝐣, and applying (2.52) together with (2.22) we obtain that lim_𝑛→∞1𝑛+1^𝑛_𝑘=1𝐚⋅𝐗^∗_𝑘−1inf_𝜃∈Θ𝑝1+^′𝜃,𝐗^∗_𝑘−1𝐚⋅𝐗^∗_𝑘−1^2a.s.=_{𝐣∈ℕ𝑑}𝐚⋅𝐣inf_𝜃∈Θ𝑝1+^′(𝜃,𝐣)𝐚⋅𝐣²𝜋_𝜃0(𝐣). (2.55) Let 𝐣≠𝟎 fixed. Since, for all 𝜃∈Θ, 𝑝^′(𝜃,𝐣)≠−𝐚⋅𝐣, the extreme value theorem implies that inf_𝜃∈Θ(1+𝑝^′(𝜃,𝐣)/𝐚⋅𝐣)²>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 𝜃^∗_𝑛−𝜃₀. For this purpose, we follow the steps of the proof of Proposition 6.1 in [17]. Writing the Taylor expansion of 𝑆^′_𝑛(𝜃^∗_𝑛) at 𝜃₀ we obtain that 𝜃^∗_𝑛−𝜃₀=−𝑆^′_𝑛(𝜃₀)/𝑆_𝑛^′′(𝜃_𝑛), for some 𝜃_𝑛=𝜃₀+𝑡_𝑛(𝜃^∗_𝑛−𝜃₀), with 𝑡_𝑛∈]0,1[. Since 𝑆^′_𝑛(𝜃₀)=−2∑^𝑛_𝑘=1𝜀^∗_𝑘𝑓^′(𝜃₀,𝐗^∗_𝑘−1), we can write √𝑛𝜃^∗_𝑛−𝜃₀=∑^𝑛_𝑘=1𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛𝐹_𝑛𝑛⁻¹12𝑆_𝑛^′′𝜃_𝑛𝐹_𝑛⁻¹, (2.56) where 𝐹_𝑛∶=∑^𝑛_𝑘=1(𝑓^′(𝜃₀,𝐗^∗_𝑘−1))². Let us first show that lim_𝑛→∞𝐹_𝑛𝑛^a.s.=_{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝜋_𝜃0(𝐣). (2.57) This is an application of (2.52) and (2.22), since, for all 𝐣∈ℕ^𝑑, 𝐣≠𝟎, 𝑓^′𝜃₀,𝐣²⩽𝐚⋅𝐣+2𝐶₁+𝐶²₁𝐚. (2.58) In view of (2.56), we now prove that lim_𝑛→∞𝑆_𝑛^′′𝜃_𝑛𝐹_𝑛^a.s.=2. (2.59) Computing 𝑆_𝑛^′′ thanks to the formula 𝑆_𝑛(𝜃)=∑^𝑛_𝑘=1(𝜀^∗_𝑘+𝑓(𝜃₀,𝐗^∗_𝑘−1)−𝑓(𝜃,𝐗^∗_𝑘−1))², it appears that (2.59) is true, as soon as the following holds: lim_𝑛→∞sup_𝜃∈Θ∑^𝑛_𝑘=1𝜀^∗_𝑘𝑓^′′𝜃,𝐗^∗_𝑘−1𝐹_𝑛^a.s.=0,lim_𝑛→∞∑^𝑛_𝑘=1𝑓^′𝜃_𝑛,𝐗^∗_𝑘−1²𝐹_𝑛^a.s.=1,lim_𝑛→∞∑^𝑛_𝑘=1𝑓𝜃₀,𝐗^∗_𝑘−1−𝑓𝜃_𝑛,𝐗^∗_𝑘−1𝑓^′′𝜃_𝑛,𝐗^∗_𝑘−1𝐹_𝑛^a.s.=0. (2.60)(2.61)(2.62) 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 lim_𝑛𝐹=_𝑛^a.s.∞ (as an immediate consequence of the stronger result (2.57)), and since 𝑓^′′(⋅,𝐗^∗_𝑘−1) fulfills the required Lipschitz condition. Indeed we proved earlier that 𝑢_𝑖^′′′(𝜃) is continuous on the compact set Θ and is thus bounded on Θ, which implies that 𝑓^′′′(⋅,𝐗^∗_𝑘−1)=𝑝^′′′(⋅,𝐗^∗_𝑘−1)(𝐚⋅𝐗^∗_𝑘−1)^−1/2 is bounded by a ℱ^∗_𝑘−1-measurable function. In view of (2.61), we consider the function 𝑓(𝜃,𝐣)² and its derivative 2𝑓^′(𝜃,𝐣)𝑓^′′(𝜃,𝐣). Similarly as for (2.54), one can show that there exists a constant 𝐶₂>0 such that for all 𝜃∈Θ, and all 𝐣≠𝟎, ∣𝑝^′′(𝜃,𝐣)∣⩽𝐶₂. This implies 2𝑓^′(𝜃,𝐣)𝑓^′′(𝜃,𝐣)⩽2𝐚⋅𝐣+𝑝^′(𝜃,𝐣)𝑝^′′(𝜃,𝐣)𝐚⋅𝐣⩽2𝐶₂𝐶1+¹𝐚. (2.63) Consequently, ∑^𝑛_𝑘=1𝑓^′𝜃_𝑛,𝐗^∗_𝑘−1²−𝑓^′𝜃₀,𝐗^∗_𝑘−1²𝐹_𝑛⩽2𝐶₂𝐶1+¹𝐚𝜃^∗_𝑛−𝜃₀𝐹_𝑛𝑛⁻¹, (2.64) which by (2.57) and the strong consistency of 𝜃^∗_𝑛 almost surely tends to 0. Writing ∑^𝑛_𝑘=1𝑓𝜃_𝑛,𝐗^∗_𝑘−1²𝐹_𝑛∑=1+^𝑛𝑘=1𝑓^′𝜃_𝑛,𝐗^∗_𝑘−1²−𝑓^′𝜃₀,𝐗^∗_𝑘−1²𝐹_𝑛, (2.65) this implies (2.61). It now remains to prove (2.62). We write ∑^𝑛_𝑘=1𝑓𝜃₀,𝐗^∗_𝑘−1−𝑓𝜃_𝑛,𝐗^∗_𝑘−1𝑓^′′𝜃_𝑛,𝐗^∗_𝑘−1𝐹_𝑛⩽1𝐹_𝑛^𝑛_𝑘=1𝐚⋅𝐗^∗_𝑘−1𝜃₀−𝜃_𝑛+𝑝𝜃₀,𝐗^∗_𝑘−1−𝑝𝜃_𝑛,𝐗^∗_𝑘−1𝐚⋅𝐗^∗_𝑘−1𝑝^′′𝜃_𝑛,𝐗^∗_𝑘−1⩽1𝐹_𝑛^𝑛_𝑘=1𝜃₀−𝜃_𝑛+𝐶₁𝜃₀−𝜃_𝑛𝑎𝐶₂⩽𝜃₀−𝜃^∗_𝑛𝐶₂𝐶1+¹𝑎𝐹_𝑛𝑛⁻¹, (2.66) which thanks to (2.57) and the strong consistency of 𝜃^∗_𝑛 implies (2.62). In view of (2.56), we finally want to prove that ∑^𝑛_𝑘=1𝜀^∗_𝑘𝑓^′(𝜃₀,𝐗^∗_𝑘−1)/√𝑛 converges in distribution. For this purpose we make use of a central limit theorem for martingale difference arrays [18, 19]; if {𝑀_𝑘^(𝑛),ℱ_𝑘^(𝑛),1⩽𝑘⩽𝑛}, 𝑛⩾1 is a sequence of square integrable martingales with associated Meyer process ⟨𝑀⟩^(𝑛)=(⟨𝑀⟩_𝑘^(𝑛))_{1⩽𝑘⩽𝑛} satisfying lim_𝑛→∞⟨𝑀_𝑛⟩^𝑃(𝑛)=𝑐² for some constant 𝑐, and such that, for all 𝜀>0, lim^𝑛_{𝑛→∞𝑘=1}𝔼𝑀_𝑘^(𝑛)−𝑀^(𝑛)_𝑘−1²𝟏_{|𝑀𝑘^(𝑛)_−𝑀^(𝑛)_{𝑘−1|⩾𝜀}}∣ℱ=^(𝑛)_𝑘−1^𝑃0, (2.67) then lim_𝑛→∞𝑀_𝑛^𝒟(𝑛)=𝒩(0,𝑐²). Let us define 𝑀_𝑘^(𝑛)∶=∑^𝑘_𝑙=1𝜀^∗_𝑙𝑓^′(𝜃₀,𝐗^∗_𝑙−1)/√𝑛, for every 𝑘⩽𝑛. First, for any 𝑘⩽𝑛, 𝔼_𝜃0(𝜀^∗_𝑘𝑓^′(𝜃₀,𝐗^∗_𝑘−1)/√𝑛∣𝐗^∗_𝑘−1)=0. Second, 𝔼_𝜃0𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛²∣𝐗^∗_𝑘−1=𝑓^′𝜃₀,𝐗^∗_𝑘−1²𝑔𝜃₀,𝐗^∗_𝑘−1𝑛, (2.68) hence {𝑀_𝑘^(𝑛)}_𝑘⩽𝑛 is a sequence of square integrable martingales. Moreover, using inequalities (2.50) and (2.58), we obtain, by (2.22), _{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣)⩽𝜃₀+𝐛+1𝐚_{𝐣∈ℕ𝑑}𝐚⋅𝐣𝜋_𝜃0(𝐣)+2𝐶₁+𝐶²₁𝐚<∞. (2.69) So, by means of (2.52), lim_𝑛→∞𝑀_𝑛^(𝑛)=lim^𝑛_{𝑛→∞𝑘=1}𝔼_𝜃0𝜀_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛²∣𝐗^∗_𝑘−1a.s.=_{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣). (2.70) Third, using Cauchy-Schwarz and Bienaymé-Chebyshev inequalities, ^𝑛_𝑘=1𝔼_𝜃0𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛²𝟏_{{|𝜀∗𝑘𝑓′(𝜃0,𝐗∗𝑘−1√)/𝑛|⩾𝜀}}∣𝐗^∗_𝑘−1⩽^𝑛_𝑘=1𝔼_𝜃0𝜀^∗_k𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛⁴∣𝐗^∗_𝑘−1^1/2ℙ_𝜃0𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛⩾𝜀∣𝐗^∗_𝑘−1^1/2⩽1𝑛^3/2𝜀^𝑛_𝑘=1𝑓^′𝜃₀,𝐗^∗_𝑘−1³𝔼_𝜃0𝜀^∗_𝑘⁴∣𝐗^∗_𝑘−1^1/2𝔼_𝜃0𝜀^∗_𝑘²∣𝐗^∗_𝑘−1^1/2. (2.71) Let us compute 𝔼_𝜃0((𝜀^∗_𝑘)⁴∣𝐗^∗_𝑘−1). We can show that the 4th central moment of the independent sum of a Poisson and a Bernoulli random variables equals 𝜇₄+6𝜇₂𝛾₂+𝛾₄, where 𝜇_𝑖 and 𝛾_𝑖 denote the 𝑖th central moment of the Poisson and of the Bernoulli variable. If these variables have parameter 𝜆 and 𝑝, respectively, then 𝜇₄=𝜆(1+3𝜆), 𝜇₂=𝜆, 𝛾₄=𝑝(1−𝑝)(3𝑝²−3𝑝+1)∈[0,1], and 𝛾₂=𝑝(1−𝑝)∈[0,1]. We thus obtain 𝔼_𝜃0𝜀^∗_𝑘⁴∣𝐗^∗_𝑘−1⩽Ψ𝜃₀⋅𝐗^∗_𝑘−17+3Ψ𝜃₀⋅𝐗^∗_𝑘−1+1𝐚⋅𝐗^∗_𝑘−1². (2.72) Hence 𝑓^∗′𝜃₀,𝐗^∗_𝑘−1³𝔼_𝜃0𝜀⁴_𝑘∣𝐗^∗_𝑘−1^1/2𝔼_𝜃0𝜀²_𝑘∣𝐗^∗_𝑘−1^1/2⩽𝐚⋅𝐗^∗_𝑘−1+𝐶₁√𝐚³Ψ𝜃₀⋅𝐗^∗_𝑘−17+3Ψ𝜃₀⋅𝐗^∗_𝑘−1+1𝐚⋅𝐗^∗_𝑘−1𝜃₀+𝐛+1𝐚^1/2. (2.73) Since the highest power of 𝑋^∗_𝑛 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 lim^𝑛_{𝑛→∞𝑘=1}𝔼_𝜃0𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛²𝟏_{{|𝜀∗𝑘𝑓′(𝜃0,𝐗∗𝑘−1√)/𝑛|⩾𝜀}}∣𝐗=^∗_𝑘−1^a.s.0. (2.74) It then ensues from the central limit theorem mentioned above that lim_𝑛→∞∑^𝑛_𝑘=1𝜀^∗_𝑘𝑓^′𝜃₀,𝐗^∗_𝑘−1√𝑛^𝒟=𝒩0,_{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣). (2.75) Finally, (2.56), together with (2.57), (2.59), (2.75), and Slutsky's theorem, implies that lim_𝑛→∞√𝑛𝜃^∗_𝑛−𝜃₀^𝒟∑=𝒩0,^{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣)∑_𝐣∈ℕ^𝑑𝑓^′𝜃₀,𝐣²𝜋_𝜃0(𝐣)², (2.76) which leads to lim_𝑛→∞𝑛∑_𝐣∈ℕ^𝑑𝑓^′𝜃₀,𝐣²𝜋_𝜃0(𝐣)²∑_𝐣∈ℕ^𝑑𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣)𝜃^∗_𝑛−𝜃₀^𝒟=𝒩(0,1). (2.77)
It now remains to prove the asymptotic distribution (2.51). Thanks to what precedes, this result is immediate as soon as we prove that lim_𝑛→∞1𝑛+1^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗=^∗_𝑘^a.s._{𝐣∈ℕ𝑑}𝑓^′𝜃₀,𝐣²𝑔𝜃₀,𝐣𝜋_𝜃0(𝐣), (2.78) as well as the equivalent result for the numerator. For this purpose, we write ^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗^∗_𝑘=^𝑛_𝑘=0𝑓^′𝜃₀,𝐗^∗_𝑘²𝑔𝜃₀,𝐗^∗_𝑘+^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗^∗_𝑘−𝑓^′𝜃₀,𝐗^∗_𝑘²𝑔𝜃₀,𝐗^∗_𝑘, (2.79) and show that (𝑓^′(⋅,𝐣))²𝑔(⋅,𝐣) has a bounded derivative and is thus Lipschitz. We have indeed ∣𝑔^′(𝜃,𝐣)∣=∣𝐚⋅𝐣+𝑝^′(𝜃,𝐣)∣(𝐚⋅𝐣)⁻¹⩽1+𝐶₁𝐚⁻¹, hence 2𝑓^′′(𝜃,𝐣)𝑓^′(𝜃,𝐣)𝑔(𝜃,𝐣)+𝑓^′(𝜃,𝐣)²𝑔^′(𝜃,𝐣)⩽2𝐶₂1+𝐶₁𝐚⁻¹𝐚⁻¹(Ψ(𝜃)⋅𝐣+1)+1+𝐶₁𝐚⁻¹²𝐚⋅𝐣+𝐶₁⩽𝐶₃Ψ𝜃₂⋅𝐣, (2.80) for some constant 𝐶₃>0. This enables us to write 1𝑛+1^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗^∗_𝑘−𝑓^′𝜃₀,𝐗^∗_𝑘²𝑔𝜃₀,𝐗^∗_𝑘⩽𝜃^∗_𝑛−𝜃₀𝐶₃1𝑛+1^𝑛_𝑘=0Ψ𝜃₂⋅𝐗^∗_𝑘. (2.81) By the strong consistency of 𝜃^∗_𝑛 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 𝐼 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 {𝑋_𝑛}_𝑛 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 𝐸 and 𝑆 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 {𝐼_𝑛}_𝑛 as the initial total population size |𝐍₀|→∞, where 𝐼_𝑛 denotes the number of new infectives at time 𝑛. More precisely, denoting by 𝐸_𝑛 the incidence of exposed individuals at time 𝑛 and by 𝑎_𝑚 the largest individual survival age, we obtain that {𝑋_𝑛,𝑍_𝑛}_𝑛^𝒟=lim_𝑁0→∞{𝐼_𝑛,𝐸_𝑛}_𝑛, with 𝑋_𝑛∣𝑋_𝑛−1,…,𝑋_{𝑛−𝑎𝑚+1}^𝒟=Poiss^𝑎𝑚−1_𝑘=1Ψ_𝑘𝑋_𝑛−𝑘,𝑍_𝑛∣𝑋_𝑛^𝒟=PoissΨ₀𝑋_𝑛, (3.1)(3.2) and where for all 𝑘=1⋯𝑎_𝑚−1, Ψ_𝑘 satisfies (1.2) in Section 1, with the latent period distribution given survival, 𝑃_inc,𝑎, assumed independent of 𝑎, being denoted by 𝑃_inc. Thus Ψ_𝑘=𝜃₀𝑃_inc(𝑘)^𝑎𝑚_𝑎=𝑘+1𝑃_age(𝑎)+𝑝_mat𝑃_inc(𝑘)𝑃_age(𝑘+1), (3.3) and similarly, Ψ₀=𝜃₀+𝑝_mat𝑃_age(1). (3.4) Moreover, when assuming a relatively stable population size and denoting by 𝑆(𝑎) the individual survival probability until age 𝑎 at least, then 𝑃_age(𝑎) is equal to 𝑆(𝑎)(∑^𝑎𝑚_𝑎′=1𝑆(𝑎^′))⁻¹.

3.2. Illustration: The BSE Epidemic in Great Britain

In this section we provide an illustration of the methodology developed in Section 2 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 {𝑋_𝑛}_𝑛 defined by (3.1)–(3.4) appears relevant to model the spread of this fatal disease after the 1988 ban, where the infectious state 𝐼 corresponds to the end of the incubation period and the clinical state, and 𝑅 corresponds to the death (assumed to be either by routine slaughtering for 𝐸 and 𝑆 cattle, and control slaughtering for 𝐼 cattle). Choosing a time step of one year, X_𝑛 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 𝜃₀ 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 𝜃₀. We provide an estimation of 𝜃₀ 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 𝑎_𝑚 is set to 𝑎_𝑚=10, which corresponds to the largest reported survival age with a nonnegligible probability. Next, we set 𝑝_mat=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 𝑘⩾1, 𝑃_inc(𝑘)∶=𝑒^{−((𝛼−1)/𝛼𝛽𝛼)(𝑘−1)𝛼}−𝑒^{−((𝛼−1)/𝛼𝛽𝛼)𝑘𝛼}, 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 𝑆(𝑘) (probability for an apparently healthy animal to survive at least until 𝑘 years) corresponds to the probability of being slaughtered after the age of 𝑘 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 (Figure 1) provided by the World Organisation for Animal Health [5], denoted by 𝑋^obs_𝑛 for each year 𝑛. We set 𝑑=𝑎_𝑚−1=9 and denote by 𝐗^obs_𝑛 the 𝑑-dimensional vector (𝑋^obs_𝑛,…,𝑋^obs_{𝑛−𝑑+1}). We choose the first time 𝑛=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 𝑋_−𝑑+1=𝑋^obs₁₉₈₉, and 𝑛=1 corresponds to the year 1989+𝑑=1998. We estimate 𝜃₀ by the WCLSE (2.30) in model (1.1)–(3.3), where {𝑝_mat,𝛼,𝛽,{𝑆(𝑘)}_𝑘} are given by the previous values, and where 𝐗₀=𝐗^obs₁₉₉₇. According to [5], |𝐗₀|=|𝐗^obs₁₉₉₇|=167977. We are thus close to the asymptotic |𝐗₀|→∞. The number of observations is 𝑛=14. The estimator (2.30) provides the estimation 𝜃^obs_∣𝐗0∣=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∶=𝜃_∣𝐗0∣−1.96𝑐̂₁⁻¹, 𝜃_max∶=𝜃_∣𝐗0∣+1.96𝑐̂₁⁻¹, and 𝑐₁̂∶=(∑^𝑛_𝑘=1𝐚⋅𝐗_𝑘−1/𝜎²(𝜃_∣𝐗0∣))^1/2. Assuming 𝛼_𝑖=𝑋^obs_{1997−𝑖+1}|𝐗^obs₁₉₉₇|⁻¹, we get 𝑐^obs₁̂=40.7343 (observed value of 𝑐₁̂), and ℙ𝜃₀∈𝜃_min,𝜃_max≃95%,𝜃^obs_min=2.3842,𝜃^obs_max=2.4805. (3.5) Although this confidence interval is an asymptotic one, as |𝐗₀|→∞, it is a very good approximation of the true confidence interval for a finite |𝐗₀|, since |𝐗₀| is here very large. Since the estimation of 𝜃₀ relies on the values given to the other model parameters {𝑝_mat,𝛼,𝛽}, we evaluate in addition its sensitivity to the values of these parameters. For this purpose, we compute the estimation of 𝜃₀ and the associated confidence interval with asymptotic probability 95%, for different values of (𝑝_mat,𝛼,𝛽) (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 𝜃₀ 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 𝜃₀ 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, 𝜃₀ 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 𝜃₀, and its confidence interval [𝜃^obs_min,𝜃^obs_max] with asymptotic probability 95%, for different values of the maternal infection parameter 𝑝_mat 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 𝑝_mat. The estimations of 𝜃₀ are based on the observed data over the years 1989–2011.

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 Proposition 2.1 that {𝑋_𝑛}_𝑛 becomes extinct almost surely if and only if 𝑅₀=∑^𝑑_𝑘=1Ψ_𝑘(𝜃₀)⩽1. The estimated basic reproduction is here 𝑅₀(𝜃^obs_∣𝐗0∣)=0.1072. Moreover, solving with a computing program the equation ∑^𝑑_𝑘=1Ψ_𝑘(𝜃^obs_∣𝐗0∣)𝜌^−𝑘=1, we obtain the following value for the Perron's root 𝜌(𝜃^obs_∣𝐗0∣)=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 {𝑋_𝑛}_𝑛, where 𝜃₀ is replaced by its previous estimation 𝜃^obs_∣𝐗0∣=2.4324, and where the initial time of the model is 2011, that is, 𝐗₀=𝐗^obs₂₀₁₁. The simulations are done recursively using the transition law (3.1). We point out that the model initialized by 𝐗₀=𝐗^obs₁₉₉₇ 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 {𝑋_𝑛}_𝑛 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 {𝑋_𝑛}_𝑛 initialized by the observed values 𝐗^obs₂₀₁₁, with the estimated infection parameter 𝜃^obs_∣𝐗0∣=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 𝑍_𝑛 of infected cattle at time 𝑛, conditionally on the number 𝑋_𝑛 of cases at that time, is given by the Poisson distribution (3.2). For every 𝑛⩾2012 and for each of the 1000 previously simulated values 𝑋_𝑛, we generate one realization of 𝑍_𝑛. We then illustrate in Figure 2(d), the yearly maximum, minimum, median, 2.5%, and 97.5% quantiles associated with the 1000 realizations.

(a)

(b)

(c)

(d)

Figure 2 

Figure 2(a) 10 simulations of {𝑋𝑛}𝑛 initialized by 𝐗obs1997, and comparison with the observations on the period 1998–2011. Figure 2(b) 5 simulations of {𝑋𝑛}𝑛 initialized by 𝐗obs2011. 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 𝜃obs∣𝐗0∣=2.4324.

3.2.6. Prediction of the Year of Extinction

Let 𝑇∶=2011+inf{𝑛⩾1,𝐗_𝑛=𝟎} denote the extinction year of the epidemic process {𝑋_𝑛}_𝑛. According to (2.7) and denoting by 𝐟_𝜃0 the offspring generating function of {𝐗_𝑛}_𝑛 defined in (2.2) (from now on we let the dependence in 𝜃₀ appear in the notation), we have ℙ_𝜃0(𝑇⩽2011+𝑛)=𝐟_𝜃0,𝑛(𝟎)^𝐗obs2011, for every 𝑛⩾1, which by iterating 𝐟_𝜃0 can be computed explicitly. We obtain in particular, for the estimated value 𝜃^obs_∣𝐗0∣=2.4324, the following 𝑝-quantiles for the extinction time (see (2.8)): 𝑛^𝑇_0.5=2028, 𝑛^𝑇_0.95=2035 and 𝑛^𝑇_0.99=2039. Keeping in mind that 𝑇 corresponds to the complete extinction of the epidemic (i.e., 𝑑=9 consecutive years without any case), these results actually mean that for the infection parameter 𝜃₀=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 𝜃^obs_∣𝐗0∣ of the infection parameter 𝜃₀, we make use of the asymptotic confidence interval (3.5) of 𝜃₀ and of the fact that 𝜃↦ℙ_𝜃(𝑇⩽𝑛) is a decreasing function of 𝜃, which implies that for every 𝑛⩾2011, ℙℙ_𝜃0()∈ℙ𝜃𝑇⩽𝑛_max()𝜃𝑇⩽𝑛,ℙ_min(𝑇⩽n)≃95%. (3.6) We collect in Table 2 the observed interval 𝜃[ℙ^obs_max𝜃(𝑇⩽𝑛),ℙ^obs_min(𝑇⩽𝑛)], for each 𝑛⩾2020 (if 𝑛<2019 we have, conditionally on the initial value 𝐗^obs₂₀₁₁, ℙ(𝐗_𝑛=𝟎)=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(𝑇⩽𝑛).

Table 2 

Cumulative distribution function of the year of extinction computed with the infection parameters 𝜃^obs_min=2.3842 and 𝜃^obs_max=2.4805 defined by (3.5). The values in bold character correspond to the 𝑝-quantiles 𝑛^𝑇_𝑝 for 𝑝=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙ_𝜃0(𝑇⩽𝑛^𝑇_𝑝) based on (3.6).

3.2.7. Prediction of the Epidemic Size

Let 𝑁∶=∑^𝑇−2011_𝑛=1𝑋_𝑛 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 𝑁 using (2.11), conditionally on the event {𝐗₀=𝐗^obs₂₀₁₁}. We obtain in particular, for the estimated value 𝜃^obs_∣𝐗0∣=2.4324, the following 𝑝-quantiles for the total epidemic size (see (2.12)): 𝑛^𝑁_0.5=22, 𝑛^𝑁_0.95=31 and 𝑛^𝑁_0.99=35. From (2.13) we deduce that the mean value and variance of 𝑁 for the parameter 𝜃^obs_min (resp., 𝜃^obs_max) 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𝜃(𝑁⩽𝑛),ℙ_min(𝑁⩽𝑛)] of ℙ_𝜃0(𝑁⩽𝑛), for every 𝑛∈ℕ, with asymptotic probability 95%, collected in Table 3.

Table 3 

Cumulative distribution function of the total size of the epidemic for the infection parameters 𝜃^obs_min=2.3842 and 𝜃^obs_max=2.4805 defined by (3.5). The values in bold character correspond to the 𝑝-quantiles 𝑛^𝑁_𝑝 for 𝑝=0.5, 0.95, and 0.99, and to the asymptotic confidence intervals of ℙ_𝜃0(𝑁⩽𝑛^𝑁_𝑝).

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 𝜃^∗obs_𝑛 (2.47) based on the data in Great Britain. The number of available observations is only 𝑛=14, so we are far from the asymptotic setting 𝑛→∞ of Theorem 2.7. However, the large value of |𝐗₀| can make us hope for a good accuracy. We point out that, by making use of the estimator 𝜃^∗_𝑛 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 𝜃^∗obs_𝑛 should a priori be a bit smaller than the value 2.4324 provided by 𝜃^obs_∣𝐗0∣. Indeed we obtain 𝜃^∗obs_𝑛=2.4305. Using this value, we deduce from (2.51) the confidence interval [𝜃^∗_min,𝜃^∗_max] of 𝜃₀ with asymptotic probability 95%, where 𝜃^∗_min∶=𝜃^∗_𝑛−1.96𝑐̂₂⁻¹, 𝜃^∗_max∶=𝜃^∗_𝑛+1.96𝑐̂₂⁻¹, and 𝑐₂̂∶=^𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘^2𝑛_𝑘=0𝑓^′𝜃^∗_𝑛,𝐗^∗_𝑘²𝑔𝜃^∗_𝑛,𝐗^∗_𝑘^−1/2. (3.7) Therefore, 𝑐^obs₂̂=40.6988, and ℙ𝜃₀∈𝜃^∗_min,𝜃^∗_max≃95%,𝜃^∗obs_min=2.3823,𝜃^∗obs_max=2.4787, (3.8) 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 𝑘 to the 𝑄-process as 𝑘→∞ is not very fast. As a consequence the study of the 𝑄-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 𝐗₀=𝐗^obs₁₉₉₇, and where 𝜃₀ is estimated by 𝜃^∗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 𝐗^∗₀=𝐗^obs₂₀₁₁. It appears that the values of this simulated trajectory are rapidly very small and of course are never equal to 0 for 𝑑=9 consecutive times. For a finer prediction, we simulate 1000 realizations of this process from 2012 until 2050, with 𝐗^∗₀=𝐗^obs₂₀₁₁. Moreover, for every 𝑛⩾2012 and for each of the 1000 simulated values 𝑋^∗_𝑛, we make one realization of the incidence 𝑍_𝑛 of infected cattle at time 𝑛, 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.

(a)

(b)

(c)

(d)

Figure 3 

Figure 3(a) 10 simulations of {𝑋∗𝑛}𝑛 initialized by 𝐗obs1997, and comparison with the observations on the period 1998–2011. Figure 3(b) one simulation of {𝑋∗𝑛}𝑛 initialized by 𝐗obs2011. 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 𝜃∗obs∣𝐗0∣=2.4305.

4. Conclusion

The methodology developed in Section 2 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.