1.Introduction SophieP´enisson andChristineJacob StochasticMethodologyfortheStudyofanEpidemicDecayPhase,BasedonaBranchingModel ResearchArticle

Volume 2012, Article ID 598701,32pages doi:10.1155/2012/598701

Research Article

Stochastic Methodology for

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

Sophie P ´enisson

and Christine Jacob

1Universit´e Paris-Est, LAMA (UMR 8050), UPEMLV, UPEC, CNRS, 94010 Cr´eteil, France

2INRA, MIA (UR 341), 78352 Jouy-en-Josas, France

Correspondence should be addressed to Sophie P´enisson,sophie.penisson@u-pec.fr Received 16 July 2012; Revised 10 October 2012; Accepted 10 October 2012

Academic Editor: Charles J. Mode

Copyrightq2012 S. P´enisson and C. Jacob. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 periodduring 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 in1 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}_nof 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”ktime 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}_nis then the following:

X_n d k 1

Xn−k

i 1

Y_n−k,n,i, 1.1

where the variable Y_n−k,n,i is the incidence of secondary cases produced at time n with a delayklatent periodby individualiinfectious at timen−k. The{Y_n−k,n,i}_i,kare assumed independent givenF_n−1: σ{X_n−k}_k1, and the{Y_n−k,n,i}_iare assumed i.i.d.identically and independently distributedgivenFn−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´e-Galton-Watsonbranching 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 modelsee, e.g.,3,4and citations therein.

The core of the paper lies inSection 2, where the whole methodology is presented. We first formulate the epidemic model{Xn}_nas 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 sizeSections2.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 theQ- process{X_n^∗}_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 {X_n^∗}_n, as well as the derived distributions above, we need to know the values of the parameters {Ψk}_k. We may write Ψk am−k

a 1 θ_aP_inc,akPageak, where θ_a is the mean number of individuals infected at ageaby an infective by direct or indirect transmission,a_m is the largest survival age,Pinc,akis the probability for the individual agedaat infection to have a latent period equal tokgiven his survival, andPageakis the probability to be aged akat the end of the latent period. Parameters{θa}_aare 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 θ0pa, wherepais constant over time, whileθ0may change with control measures. A typical example is whenθ₀ is the mean number of individuals infected by an infective by horizontal route at agea, assumed independent ofa2, andp1represents the maternal transmission probabilitypmat. In this case

Ψk θ₀

am

a k1

P_inc,a−kkPagea p_matP_inc,1kPagek1. 1.2

So we assume here that, except forθ₀, the other parameters of the{Ψk}_kare constant over time and are knowngenerally estimatedfrom 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θ₀ seeSection 3.1. In Sections2.5and2.6, we provide optimal WCLSE Weighted Conditional Least Squares Estimatorsofθ₀ in the decay phase, in the frame of {Xn}_n as well as in the frame of the associated conditioned process {X^∗_n}_n, and prove the strong consistency and the asymptotic normality of these estimators.

The finalSection 3is devoted to the application of this method to real epidemics. We first present inSection 3.1some general conditions under which the spread of a SEIR disease susceptible, exposed latent, infectious, removedcan be approximated by our epidemic process defined by1.1and give an explicit derivation of the parameters{Ψk}_k. We then illustrate the methodology inSection 3.2with the decay phase of the BSEBovine Spongiform Encephalopathyepidemic in Great Britain. According to the available data5, the epidemic is obviously fading out. We assume that the{Ψk}_ksatisfy1.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 process1.1defined inSection 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 orderd. Consequently, thed-dimensional process{Xn}_ndefined by

Xn: Xn, X_n−1, . . . , X_n−d1 2.1

is Markovian of order 1, and it stems directly from1.1that{Xn}_nis a multitype Bienaym´e- Galton-Watson BGW process withd typessee, 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 biologysee, e.g.,7, but simply correspond to the memory of the process{Xn}_n. The information provided by thed-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 theory6. The first basic tool is the generating function of the offspring distribution of{Xn}_n, f : f1, . . . , fd, defined on0,1^d byf_ir: Er^X1 |X₀ e_i, where ei : 0, . . . ,1, . . . ,0denotes theith basis vector ofN^dand u^v: _d

i 1u^v_iⁱfor u,v∈N^d. For all r∈0,1^d, we have here fir: e^−1−r1Ψir_i1, i 1· · ·d−1,

fdr: e^−1−r1Ψd.

2.2

The second basic tool is the mean matrix M defined byEXn| F_n−1 X_n−1M, which is here

M

⎛

⎜⎜

⎜⎜

⎜⎜

⎝

Ψ1 1 0 · · · 0 Ψ2 0 1 · · · 0 ... ... . .. ...

Ψ_d−1 0 · · · 1 Ψd 0 · · · 0

⎞

⎟⎟

⎟⎟

⎟⎟

⎠

. 2.3

Let us notice that, sinceΨk>0 for eachk 1· · ·d, then{Xn}_nis nonsingular, positive regular see6and satisfies theXlogXcondition,

EX1lnX1 |X0 ei<∞, i 1, . . . , d, 2.4

where · denotes the sup norm inR^d. 2.1. Extinction of the Epidemic 2.1.1. Almost Sure Extinction

Since the single-type process{Xn}_nhas a memory of sized, it becomes extinct when it is null at dsuccessive times, or equivalently as soon as the d-dimensional process{Xn}_n reaches thed-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 eigenvaluealso called the Perron’s rootof the mean matrix M. Thusρis solution of _d

k 1Ψkρ^−k 1. In general ford > 1,ρ has no explicit expression. However,_d

k 1Ψkρ^−k 1 leads directly to the following explicit threshold criteria.

Proposition 2.1. The epidemic becomes extinct almost surely if and only ifR₀ 1, where R₀ : _d

k 1Ψk is the total mean number of secondary cases generated by one infective in aSEIRdisease.

We callR₀the basic reproduction number.

Moreover, whenρ 1, thenR0 ρwith equality if and only if eitherρ 1 ord 1, and whenρ >1, thenR₀ρwith equality if and only ifd 1.

Note that whend >1,R₀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 timen, fornlarge. Denoting by u and v the right and left eigenvectors of M associated to the Perron’s rootρ; that is, Mu^T ρu^T and vM ρv, with the normalization convention u·1 u·v 1, where u·v stands for the usual scalar product in R^d and where the superscript T denotes the transposition, then EXn | X₀ X₀Mⁿ ∼

n→ ∞ ρⁿX₀u^Tv. The first coordinate in the latter formula becomes for the epidemic processρⁿ_d

i 1X_−i1uiv1. Computing explicitly u and v, we obtain that for alli 1· · ·d,

ui

_d

k iΨkρ^−ki _d

j 1_d

k jΨkρ^−kj, vi ρ⁻ⁱ _d

j 1_d

k jΨkρ^−kj _d

j 1_d

k jΨkρ^−k , 2.5

which leads to the following asymptotic result:

EXn|X₀ ∼

n→ ∞ρⁿ d

i 1

X_{−i1 d}

k iΨkρ^−ki−1 _d

j 1_d

k jΨkρ^−k. 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 X₀ X0, X₋₁, . . . , X_−d1, 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. LetT : inf{n 1,X_n 0}denote the extinction time of the process {Xn}_n, and let fn : f◦f_n−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 PT n PXn 0 fn0^X0, that is to say,

PT n

fn,10_X0

fn,20_X−1

· · ·

fn,d0_X−d1

. 2.7

It can be immediately deduced from convergence results for fn0asn → ∞8, that ifρ 1, PT n ∼ 1−nη⁻¹X₀ ·u, while ifρ < 1, PT n ∼ 1−ρⁿγX₀ ·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.7enables the exact computationresp., estimationof PT nfor anynby the iterative computation of fn, X0being given, when the parameters Ψk of1.1are knownresp., estimated. Moreover, since forρ 1 the epidemic becomes extinct in an a.s. finite time andPT n _{n→ ∞}1, then for any given probabilityp ∈ 0,1 there existsn∈Nsuch thatPT np. So in practice, for anyp∈0,1,2.7enables us to compute thep-quantilen^T_p of the extinction time,

n^T_p : min

n1 : PT np

. 2.8

2.2. Total Size of the Epidemic

Under the assumption ρ 1 and the independence of the {Y_n−l,n,i}_i,l,n we previously assumed the independence of the{Y_n−l,n,i}_i,l, for eachn, we derive the distribution of the future total sizeN : _T

n 1X_n of the epidemic until its extinction, that is, the future total number of infectives until the extinction of the disease. It can be shown9that, given the initial value X0, the time origin being the same as inSection 2.1, the probability distribution ofNis

N ^D d k 1

X−k1

i 1

⎛

⎝^Yk,i

j 1

N_k,i,j

⎞

⎠, Y_k,i: ^d

l k

Y_{−k1,−k1l,i}, 2.9

where^Ddenotes the equality in distribution, an empty sum is by convention 0, the{Yk,i}_k,i are independent, and the{Yk,i}_iand the{Nk,i,j}_k,i,jare i.i.d. with

Yk,iD

Poiss _d

l k

Ψl

, Nk,i,j D

Borel−Tanner _d

l 1

Ψl,1

, 2.10

that is, for eachn1,PNk,i,j n e^{−ndl 1Ψl}n_d

l 1Ψlⁿ⁻¹n!⁻¹. Consequently, under the convention that an empty product is 1, the probability distribution ofNis, for anyn∈N,

PN n

{0yk,in,{1nk,i,jn}j}i,k:

_d

k 1X−k1

i 1 yk,i

j 1nk,i,j n

d k 1

X−k1

i 1

e^{−dl kΨl} _d

l kΨl

_yk,i

yk,i!

×

yk,i

j 1

e^{−nk,i,jdl 1Ψl}

n_k,i,jd

l 1Ψl

_nk,i,j−1

n_k,i,j! ,

2.11

which may be calculated resp., estimated, replacing the Ψk by their values resp., estimations. In practice, for anyp ∈ 0,1,2.11enables to compute thep-quantilen^N_p of the total epidemic size,

n^N_p : min

n1 : PNnp

. 2.12

We obtain moreover an explicit formula for the mean value and variance of the size of the epidemic,

EN _d

k 1X_−k1d

l kΨl

1−_d

l 1Ψl

, VarN _d

k 1X_−k1d

l kΨl

1−_d

l 1Ψl

₃ . 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 {Zn}_n defined by the conditional distributionZ_n | X_n ^D PoissΨ0X_n, 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 epidemici.e., the incidence of infectives as well as of exposed individualsthanks 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 thep-quantilen^T_pof the extinction time with the probabilitypas large as wantedsee 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 theunlikelycase 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 diseasewould 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, wherek is very large. We therefore consider for anyn₁, n₂, . . . ∈ Nand any i₀,i₁,i₂, . . . ∈ N^d the conditioned probabilityPi0Xn1 i1, . . . ,Xnr ir | Xk/0, where the subscript i0

denotes the initial value. Ifkis finite this distribution cannot be easily handled due to its time

inhomogeneity. However, whenρ 1, it is known10that this conditioned distribution converges, ask → ∞, to the distribution of ad-dimensional Markov process{X^∗_n}_n:

k→ ∞limPi0Xn1 i₁, . . . ,X_nr i_r | X_k/0 Pi0

X^∗_n1 i₁, . . . ,X^∗_nr i_r

. 2.14

We will further discuss inProposition 2.5 the relevancy of approximating the conditioned probability forkfixed by the limiting object2.14. The conditioned process{X^∗_n}_ndefined by2.14is known in the literature as theQ-process associated with{Xn}_n, also described as the process conditioned on “not being extinct in the distant future.” It has the following transition probability10: for everyn1, i,j∈N^d, i/0,

P

X^∗_n j | X^∗_n−1 i 1 ρ

j·u

i·uPXn j | X_n−1 i, 2.15 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 processX_n^∗: X^∗_n,1. By construction we then have X_n,i^∗ X_n−i1^∗ , for eachnand eachi 1· · ·d.

Proposition 2.2. The stochastic process{X^∗_n}_nis, conditionally on its past, distributed as the sum of two independent Poisson and Bernoulli random variables:

X_n^∗ |X^∗_n−1^DPoiss

X^∗_n−1·Ψ

∗ B p

X^∗_n−1

, 2.16

whereΨ: Ψ1, . . . ,Ψd,∗is the convolution product symbol, and

p X^∗_n−1

: u₁X^∗_n−1·Ψ u₁X^∗_n−1·Ψ _d

k 2X^∗_n−k1u_k. 2.17

Proof. Applying1.1and2.15, we obtain that for allj ∈N, P

X_n^∗ j|X^∗_n−1 P

X^∗_n

j, X^∗_n−1, . . . , X_n−d−1^∗

|X^∗_n−1 ju_1d

k 2X^∗_n−k1u_k ρX^∗_n−1·u P

X_n

j, X_n−1^∗ , . . . , X_n−d−1^∗

|X_n−1 X^∗_n−1

ju1_d

k 2X^∗_n−k1uk

ρX^∗_n−1·u

X^∗_n−1·Ψ_j j! e^{−X∗n−1·Ψ} u₁X^∗_n−1·Ψ

ρX^∗_n−1·u

X^∗_n−1·Ψ_j−1 j−1

! e^{−X∗n−1·Ψ} _d

k 2X_n−k1^∗ u_k ρX^∗_n−1·u

X^∗_n−1·Ψj

j! e^{−X∗n−1·Ψ}.

2.18

The equality Mu^T ρu^T implies that for allk 1· · ·d−1, ρuk Ψku1u_k1, and that ρu_d Ψdu₁. Consequently,ρX^∗_n−1·u u₁X^∗_n−1·Ψ _d

k 2X^∗_n−k1u_k, and thus P

X^∗_n j|X^∗_n−1

u₁X^∗_n−1·Ψ u1X^∗_n−1·Ψ _d

k 2X_n−k1^∗ uk

X^∗_n−1·Ψ_j−1 j−1

! e^{−X∗n−1·Ψ}

1− u₁X^∗_n−1·Ψ u1X^∗_n−1·Ψ _d

k 2X_n−k1^∗ uk

X^∗_n−1·Ψ_j j! e^{−X∗n−1·Ψ}

Poiss

X^∗_n−1·Ψ

∗ B

u1X^∗_n−1·Ψ u₁X^∗_n−1·Ψ _d

k 2X^∗_n−k1u_k

j .

2.19

Remark 2.3. Note that if one compares 2.16 with the transition probability of the unconditioned processX_n | X_n−1 ^D PoissXn−1·Ψ, it appears that{X^∗_n}_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 X_n−1^∗ · · · X^∗_n−d−1 0, then according to2.17,pX^∗_n−1 1, which implies that at timen, the probability to add one unit is equal to one, thus preventing the extinction of the process.

Proposition 2.4. The process{X^∗_n}_nadmits a stationary probability measureπwith finite first- and second-order moments.

Proof. Since the multitype branching process{Xn}_n satisfies property2.4, it is known10 that theQ-process{X^∗_n}_nis positive recurrent with a stationary probability measureπ given by

πi: i·uνi

k∈N^dk·uνk, i∈N^d, 2.20

whereνis the Yaglom distribution of the process{Xn}_n, uniquely defined by the following property: for all i,j∈N^d\ {0}, lim_{n→ ∞}PXn i|X₀ j,X_n/0 νi. In the literature, this stationary measure for the conditioned process{X^∗_n}_nis also referred to as the doubly limiting conditional probability. Moreover, byProposition 2.2,

EX^∗_n E E

X^∗_n|X^∗_n−1 E

X^∗_n−1·Ψ p X^∗_n−1

^d

k 1

E X^∗_n−k

Ψk1, 2.21

which implies that lim_{n→ ∞}EX_n^∗1−_d

k 1Ψk⁻¹ <∞. We consequently obtain by means of Fatou’s lemma that, for everyi 1· · ·d,

j∈N^d

j_iπj E

nlim→ ∞X_n,i^∗

E

nlim→ ∞X_n−i1^∗

E

nlim→ ∞X_n^∗

lim

n→ ∞EX_n^∗<∞. 2.22

We similarly prove thatπhas finite second-order moments by writing

Var

X_n^∗|X^∗_n−1

X^∗_n−1·Ψ u₁X^∗_n−1·Ψ_d

k 2X^∗_n−k1u_k

u₁X^∗_n−1·Ψ _d

k 2X_n−k1^∗ u_k2 X^∗_n−1·Ψ 1

4. 2.23

Let us discuss the relevancy of approximating the epidemic process{Xn}_nconditioned on nonextinction at some finite timek, forklarge, by theQ-process{X^∗_n}_nobtained 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 valuekfor the survival of the disease in the population. We therefore might consider thatk is chosen large enough such that the approximation of the process{Xn}_nconditioned on the event{Xk/0}by the process{X^∗_n}_nis valid. Of course, the order of magnitude of suchkwill depend on the rate of convergence of the conditioned process to{X^∗_n}_n.

Proposition 2.5. Let n1 · · · nr k and i0, . . . ,ir ∈ N^d \ {0}. Then the difference

|Pi0Xn1 i₁, . . . ,X_nr i_r | X_k/0−Pi0X^∗_n1 i₁, . . . ,X^∗_nr i_r|decreases, ask → ∞, with max{k^s−1|λ|ρ^−1k/2, ρ^k/2}, whereλis an eigenvalue of M such thatρ > |λ| |λ3| |λ4| · · ·, with theλibeing the other eigenvalues of M. In case|λ| |λ3|we stipulate that the multiplicitysofλ is at least as great as the multiplicity ofλ₃.

Proof. Thanks to2.15and to the Markov property together with the fact thatPiXn 0 f_n0ⁱ, we have

Pi0Xn1 i1, . . . ,Xnr ir | Xk/0−Pi0

X^∗_n1 i1, . . . ,X^∗_nr ir

1−f_k−nr0^ir 1−f_k0ⁱ⁰ − 1

ρ^nr ir·u i₀·u

Pi0Xn1 i1, . . . ,Xnr ir. 2.24

The right term of2.24is known to converge to 0, ask → ∞, thanks to the property that lim_ka_k γu, for someγ > 0, where a_k : ρ^−k1−f_k0 see 8. This stems from two convergences, namely, lim_kb_k γ, whereb_k : ρ^−kv·1−f_k0, and limka_kb⁻¹_k u. Let us write ak γuεk, where limkεk 0. Sincebk v·a_kand u·v 1, it comes

akb⁻¹_k γuεk

γv·εk∼k→ ∞uγ⁻¹εk−v·εku. 2.25

It thus appears that the rate of convergence of a_k toγu is of the same order of magnitude as the one of a_kb⁻¹_k to u. Let us determine this rate in an accurate way. We use the following inequality produced by Joffe and Spitzer in8: for eachkn1,

a_kb⁻¹_k −u

1−f_k0

v·1−fk0−u

2δn_k

j k−n1αj

1−δ_n−_k

j k−n1α_j, 2.26

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 M^k, ask → ∞, presented for instance in11; we have M^k ρ^kROk^s−1|λ|^k, where R u^Tv. For the sake of clarity the symbolO·will denote either a scalar or a matrix with all the entries satisfying the associated property. This implies the existence of some constanta >0 such that, for allk∈N,1−δkR ρ^−kM^k 1δ_kR, whereδ_k : ak^s−1ρ⁻¹|λ|^k. Moreover, following8, let us write, for all r ∈ 0,1^d, 1−fr M−Er1−r, where 0 Er M, and Er O1 −r as r → 1. Thenρ⁻¹Ef_k−10 ρ⁻¹O1−f_k−10 Oρ^k−2, which implies the existence of some constantb > 0 such that, for allk ∈N, 0ρ⁻¹Efk−10 αkR, withαk : bρ^k−2. We thus have provided an explicit formula for the sequencesδk_kandαk_kintroduced by Joffe and Spitzer in8. Finally let us apply2.26ton k/2and replace in this inequalityδ_n andαnby their explicit expressions that we got. We obtain that, for allk∈N,

1−fk0

v·1−f_k0−u

2^2−sak^s−1

|λ|/ρ_k/2

b/ρ 1−ρ

ρ^k/2 1−2^1−sak^s−1

|λ|/ρ_k/2

− b/ρ

1−ρ

ρ^k/2. 2.27

Consequently, the right member of 2.24 will decrease with max{k^s−1|λ|ρ^−1k/2, bρ1 − ρ⁻¹ρ^k/2}, ask → ∞.

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 longseeSection 2.1.

2.5. Estimation of the Infection Parameter

We assume for this subsection that the parametersΨkof the epidemic model1.1–2.1are not entirely known. More precisely, we assume that theΨkare of the form, for allk 1· · ·d,

Ψkθ0 a_kθ₀b_k, 2.28

wherea_k>0 andb_k0 are constants, andθ₀is an unknown real parameter. We will write in what followsΨθ0 aθ₀b, whereΨθ0: Ψ1θ0, . . . ,Ψdθ0and so forth. This general assumption corresponds in particular to the case where theΨkare of the form1.2.

We estimate θ₀ by the following WCLSE in model 1.1–2.1. This estimator generalizes the well-known Harris estimator12in a BGW process. LetΘ : θ1, θ₂,θ₂ >

θ1 >0, such thatθ0 ∈Θ. The WCLSE is based on the normalized processYn: Xn/ a·X_n−1 and is defined by

θ_|X0|: arg min

θ∈Θ

n k 1

Yk−EθYk|X_k−1² arg min

θ∈Θ

n k 1

Xk−Ψθ·X_k−1²

a·X_k−1 . 2.29

We easily derive the following explicit form:

θ_{|X0| n}

k 1Xk−b·X_{k−1 n}

k 1a·X_k−1 . 2.30

The normalization of the processXnby

a·X_n−1appears to be the most natural and suitable for the following reasons. First, this normalization generalizes the normalizationX_n/

aX_n−1 in the monotype case, which is the one leading to the Harris estimatormofm₀ aθ₀bsince we have, ford 1,aθX0b 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 : minixiand x : maxixi, we have

θ0b a Eθ0

Yk−Eθ0Yk|X_k−1² |X_k−1

θ0b·X_k−1

a·X_k−1 θ0b

a, 2.31

hence the conditional variance of the error termYk−Eθ0Yk|X_k−1in the stochastic regression equationY_k Eθ0Yk|X_k−1Yk−Eθ0Yk|X_k−1is invariant under multiplication of the whole process, and bounded respectively to{Xn}_n, leading to the quasi-optimality ofθ_|X0|at finite

|X0|andn, in the sense of13.

Let us provide asymptotic results for the estimatorθ_|X0|defined by2.30, as the initial population size|X0| X₀X₋₁· · ·X_−d1tends to infinity. We denote bym^k_ij θthei, jth entry in thekth power of the matrix Mθgiven by2.3.

Theorem 2.6. Let us assume that, for each i 1· · ·d, there exists some αi ∈ 0,1 such that lim_{|X0| → ∞}X0,i|X0|⁻¹ αi. Then θ_|X0| is strongly consistent, that is, lim_{|X0| → ∞}θ_|X0| ^a.s. θ0, and is asymptotically normally distributed:

|Xlim0| → ∞

_n

k 1a·X_k−1 σ²

θ_|X0|

θ_|X0|−θ0

_D

N0,1, 2.32

where

σ²θ: θ _n

k 1

_d

j 1_d

i 1α_jb_im^k−1_ji θ _n

k 1_d

j 1_d

i 1αjbim^k−1_ji θ. 2.33

Proof. Let us first prove that, for eachk 1· · ·nand eachi 1· · ·d,

|Xlim0| → ∞

X_k,i

|X0|

a.s.d j 1

α_jm^k_ji θ0. 2.34

Using the branching property of the process{Xn}_n6we write

X_k,i

X0,1

j 1

X_k,i,j¹ · · ·

X0,d

j 1

X_k,i,j^d, 2.35

where, for all l 1· · ·dandj 1· · ·X0,l,X^l_k,i,j is theith coordinate of ad-type branching process at timek initialized by a single particle of typel. For k,i, and l fixed the random

variables{X^l_k,i,j}_j are i.i.d. with mean valuem^k_li θ0. According to the strong law of large numbers and under the theorem assumption, we have, for everyl 1· · ·dsuch thatX_0,l/0,

lim

|X0|→ ∞

X0,l

j 1X_k,i,j^l X0,l

a.s.m^k_li θ0, 2.36

which together with the theorem assumption leads to2.34.

To prove the consistency ofθ_|X0|we apply2.34to2.30, using the fact thatX_k X_k,1 andX_k−i X_k−1,i, and obtain

|Xlim0| → ∞θ_|X0|^a.s.

_n

k 1_d

j 1α_j

m^k_j1 θ0−_d

i 1b_im^k−1_ji θ0 _n

k 1_d

i 1_d

j 1a_iα_jm^k−1_ji θ0 . 2.37

By definition,

m^k_j1 θ0 ^d

i 1

m^k−1_ji θ0mi1θ0 ^d

i 1

m^k−1_ji θ0aiθ0bi, 2.38

hence2.37immediately leads to the strong consistency.

We are now interested in the asymptotic distribution ofθ_|X0|−θ0. We derive from2.30 that

ⁿ

k 1

a·X_k−1

θ_|X0|−θ_{0 n}

k 1Xk−Ψθ0·X_k−1

!_n

k 1a·X_k−1 . 2.39

By1.1,

Xk−Ψθ0·X_k−1 ^d

i 1 Xk−i

j 1

Y_k−i,k,j−Ψiθ0 :^d

i 1 Xk−i

j 1

Y◦_k−i,k,j, 2.40

where the{Y_k−i,k,j}_j are i.i.d. givenF_k−1, following a Poisson distribution with parameter Ψiθ0, and the {Y_k−i,k,j}_i,j are independent given F_k−1. Renumbering theY^◦_k−i,k,j we then obtain

n k 1

Xk−Ψθ₀·X_k−1 ^d

i 1 n

k 1Xk−i

j 1

Y◦_k−i,k,j. 2.41