A generating function approach to Markov chains undergoing binomial catastrophes
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n A GENERATING FUNCTION APPROACH TO MARKOVCHAINS UNDERGOING BINOMIAL CATASTROPHES
BRANDA GONCALVES AND THIERRY HUILLET
Abstract.
In a Markov chain population model subject to catastrophes, ran-dom immigration events (birth), promoting growth, are in balance with theeffect of binomial catastrophes that cause recurrent mass removal (death). Us-ing a generating function approach, we study two versions of such populationmodels when the binomial catastrophic events are of a slightly different ran-dom nature. In both cases, we describe the subtle balance between the twobirth and death conflicting effects.
Keywords:
Markov chain, random walk, random population growth, bi-nomial catastrophe events, recurrence/transience transition, height and lengthof excursions, extinction events, time to extinction, first passage, harmonicsequence, Green kernel.
PACS
MSC primary 60J10, secondary 92A15 Introduction, motivations and background
In some simple growth process, some organism changes the amount of its constitu-tive cells at random as follows. In the course of its lifetime, the organism alternatesat random busy and idle periods. In a busy growth period, it produces a random orfixed amount of new cells which are being added to its current stock of cells (say, theincremental or batch cells). In an idle catastrophic period, the organism stands atrisk being subject say to virus/bacteria attacks or radiation, resulting in each of itsconstitutive cells being either lysed or lost after mutations with some fixed mortal-ity probability, independently of each other. Under such a catastrophic event, thepopulation size is thus reduced according to a binomial distribution with survivalprobability say c ; hence the name binomial catastrophes. It will happen that thecurrent size of the organism is reduced to zero at some random time. In a worstdisaster scenario for instance, all cells can be lysed in a single idle period leadinginstantaneously to a first disastrous extinction event. From a first extinction event,the organism can then either recover taking advantage of a subsequent busy epochand starting afresh from zero, or not, being stuck to zero for ever. In the first case,we shall speak of a local extinction, state zero being reflecting, while in the lattercase eventual or global extinction is at stake, state zero being absorbing.As just described in words, the process under concern turns out to be a Markovchain on the non-negative integers, displaying a subtle balance between generalizedbirth and death events. Whenever it is ergodic, there are infinitely many localextinctions and the pieces of sample paths separating consecutive passages to zero are called excursions, the height and length of which are of some relevance in theunderstanding of the population size evolution. Excursions indeed form indepen-dent and identically distributed (iid) blocks of this Markov chain.Stochastic models subject to catastrophes have a wide application in different fieldsviz. bioscience, marketing, ecology, computer and natural sciences, etc. For in-stance,- In recent years, the study of catastrophe models has attracted much attentiondue to their wide application in computer-communication networks and digitaltelecommunication systems where interruptions due to various types of virus attacksare referred to as catastrophes. Here, unfriendly events, like virus attacks, resultin abrupt changes in the state of the system caused by the removal of some of itselements (packets), posing a major threat to these queueing systems. Continuous-time queueing models subject to total disasters have also been analyzed by a fewresearchers such as [5], [6] . A detailed survey of catastrophic events occurring incommunication networks has been carried out in [11].- In oversimplistic earthquakes’ physics models, the state of the process may beviewed as the total energy embodied in the earth’s crust system, obtained as thesum of the energies of its constitutive blocks, each carrying say a fixed unit quan-tum of energy. A busy period corresponds to an accumulation of energy epochin the system, while an idle period yields a stress release event. More realisticgrowth/collapse or decay/surge models in the same spirit, although in the contin-uum, were considered in [15], [16] and [17], with physical applications in mind,including stress release mechanisms in the physics of earthquakes.- In forestry, a good management of the biomass depends on the prediction of howsteady growth periods alternate with tougher periods, due to the occurrence ofcataclysms such as hurricanes or droughts or floods, hitting each tree in a similarway.- In Economy, pomp periods often alternate with periods of scarcity when a crisisequally strikes all the economic agents (the Joseph and Noah effects).Some aspects of related catastrophe models were recently addressed in [4], [1], [13],[14] and [8], chiefly in continuous times. For instance, [13] analyzed a continuous-time binomial catastrophe model wherein individuals arrive according to a com-pound Poisson process while catastrophes occur according to a renewal process. Heobtained the distribution of the population size at post-catastrophe, arbitrary andpre-catastrophe epochs. In [1], a population model with the immigration processsubject to binomial and geometric catastrophes has been investigated; the authorsobtained the removed population size and maximum population size between twoextinctions (height of excursions). They also discussed the first extinction time andthe survival time of a tagged individual. In [14] also, several approaches for thetransient analysis of a total disaster model were introduced with an extension ofthe methodologies to the binomial catastrophe model. The authors of [23] carriedout the transient analysis of a binomial catastrophe model with birth/immigration-death processes and obtained explicit expressions for factorial moments, which are
INOMIAL CATASTROPHES 3 further used to develop the population size distribution. Results for the equilib-rium moments and population size distributions are also given. Finally, recently,[38] obtained the stationary queue length of the M X /M/ ∞ queue with binomialcatastrophes in the heavy traffic system, using a central limit theorem.Let us finally mention that not all aspects of binomial catastrophes are included inthe above formulation. In [35], a density-dependent catastrophe model of threat-ened population is investigated. Density-dependence is shown here to have a signif-icant effect on population persistence (mean time to extinction), with a decreasingmean persistence time at large initial population sizes and causing a relative in-crease at intermediate sizes. A density-dependent model involving total disasterswas designed in [21].Discrete-time random population dynamics with catastrophes balanced by randomgrowth has a long history in the literature, starting with [25]. Mathematically, bino-mial catastrophe models are Markov chains (MCs) which are random walks on thenon-negative integers (as a semigroup), so differing from standard random walkson the integers (as a group) in that a one-step move down from some positive inte-ger cannot take one to a negative state, resulting in transition probabilities beingstate-dependent, [12]. Such MCs may thus be viewed as generalized birth and deathchains. The transient and equilibrium behaviour of such stochastic population pro-cesses with either disastrous or mild binomial catastrophes is one of the purposesof this work. We aim at studying the equilibrium distribution of this process andderiving procedures for its approximate computation. Another issue of importanceconcerns the measures of the risk of extinction, first extinction time, time elapsedbetween two consecutive extinction times and maximum population size reached inbetween.The detailed structure of the manuscript, attempting to realize this program, canbe summarized as follows: • Section 2 is designed to introduce the model in probabilistic terms. We developthree particular important cases:- Survival probability c = 1. In this case, on a catastrophic event,the chain remains in its current state with no depletion of cells atall.- Survival probability c = 0. This is a case of total disastersfor which, on a catastrophic event, the chain is instantaneouslypropelled to state 0.- the semi-stochastic growth/collapse scenario when the adjunc-tion of incremental cells, on a growth event, is deterministic beingreduced to a single element, and c = { , } so that binomial mor-tality is not degenerate. • In Section 3, we discuss the conditions under which the chain is recurrent (positiveor null) or transient. We emphasize that positive recurrence is generic, unless thebatch random variables have unrealistic very heavy tails. We use a generatingfunction approach which is well-suited to this analysis. When, as in the positive
BRANDA GONCALVES AND THIERRY HUILLET recurrent case, it is non trivial, we discuss the shape of the invariant probabilitymass function of the chain. Specifically,- When c ∈ [0 , . We refer to [31] for the notions of discrete infi-nite divisibility (compounding Poisson), self-decomposability andunimodality.- When c = 0, (total disasters), the invariant probability massfunction is the one of a geometric sum of the incremental (batch)random variables, so in the compound Poisson class, but not neces-sarily self-decomposable. We give some sufficient conditions underwhich it is self-decomposable, so unimodal and then with modenecessarily at the origin. • In Section 4, we analyze some important characteristics of the chain in the positiverecurrent regime, making use of its Green (potential) kernel, including:- from the Green (potential) kernel at (0 , x , x > . We show that this randomvariable has geometric tails. Some representation results on themean (persistence) time to extinction are supplied.- from the Green (potential) kernel at ( x , y ): first hitting timeof state y when starting from x . - Average position of the chain started at x > • In Section 5, we sketch some brief insight on the way the parameters of the chaincould be estimated from an N − sample of observations. • In Section 6, we study an absorbed version of the binomial catastrophe chain inview of further analyzing the quasi-distribution of the process conditioned on beingcurrently not extinct. • Finally, a variant of the binomial catastrophe model is introduced and studied inSection 7. It was motivated by the following observation: suppose at each step of itsevolution, the size of some organism either grows by a random number of batch cells(not reduced to a fixed number) or shrinks deterministically by only one unit (asemi-stochastic decay/surge scenario ‘dual’ to the semi-stochastic growth/collapseone). The above standard binomial catastrophe model is not able to represent a sce-nario where at least one individual is removed from the population at catastrophicevents. To remedy this, we therefore define and analyze a variant of the standardbinomial model where, on a catastrophic event, the binomial shrinking mechanismapplies to all currently alive cells but one which is systematically removed. Wealso take control on the future of the population once it hits state zero. Despite
INOMIAL CATASTROPHES 5 the apparent small changes in the definition of this modified binomial catastropheMarkov chain, the impact on its asymptotic behavior is shown to be drastic.2.
The binomial catastrophe model
We first describe the model as a toy model one.2.1.
The model.
Consider a discrete-time MC X n taking values in N := { , , , ... } . With b n ( c ), n = 1 , , ... an independent identically distributed (iid) sequence ofBernoulli random variables (rvs) with success parameter c ∈ (0 , c ◦ X n = P X n m =1 b m ( c ) denote the Bernoulli thinning of X n . Let β n , n = 1 , , ... be an iidbirth sequence of rvs with values in N = { , , ... } . The dynamics of the MC underconcern here is a balance between birth and death events according as ( p + q = 1): X n +1 = 1 ◦ X n + β n +1 = X n + β n +1 with probability pX n +1 = c ◦ X n with probability q = 1 − p. This model was considered by [25], [2]. We have c ◦ X n ∼ bin( X n , c ) hence the namebinomial catastrophe. The binomial effect is appropriate when, on a catastrophicevent, the individuals of the current population each die (with probability 1 − c ) orsurvive (with probability c ) in an independent and even way, resulting in a drasticdepletion of individuals at each step. Owing to: c ◦ X n = X n − (1 − c ) ◦ X n , thenumber of stepwise removed individuals is (1 − c ) ◦ X n with probability (wp) q. This way of depleting the population size (at shrinkage times) by a fixed randomfraction c of its current size is very drastic, especially if X n happens to be large.Unless c is very close to 1 in which case depletion is modest (the case c = 1 isdiscussed below), it is very unlikely that the size of the upward moves will be largeenough to compensate depletion while producing a transient chain drifting at ∞ .We will make this very precise below.Note also X n = 0 ⇒ X n +1 = β n +1 wp p , X n +1 = 0 wp q (translating reflection of X n at 0 if q > . Let B n ( p ), n = 1 , , ... be an iid sequence of Bernoulli rvs with P ( B ( p ) = 1) = p , P ( B ( p ) = 0) = q. Let C n ( p, c ) = B n ( p ) + cB n ( p ) = c B n ( p ) where B n ( p ) = 1 − B n ( p ). The above process’ dynamics (driven by B n ( p ) β n ) is compactly equivalentto X n +1 = c B n +1 ( p ) ◦ X n + B n +1 ( p ) β n +1 , X = 0= B n +1 ( p ) · c ◦ X n + B n +1 ( p ) (cid:0) X n + β n +1 (cid:1) , (cid:0) X n , β n +1 (cid:1) being mutually independent. The thinning coefficients are now c B n +1 ( p ) , so random.With b x := P ( β = x ), x ≥ d x,y := (cid:0) xy (cid:1) c y (1 − c ) x − y the binomial probabilitymass function (pmf), the one-step-transition matrix P of the MC X n is given by: P (0 ,
0) = q, P (0 , y ) = p P ( β = y ) = pb y , y ≥ In words, the thinning operation acting on a discrete random variable is the natural discreteanalog of scaling a continuous variable, i.e., multiplying it by a constant in [0 , BRANDA GONCALVES AND THIERRY HUILLET (1) P ( x, y ) = q (cid:18) xy (cid:19) c y (1 − c ) x − y = qd x,y , x ≥ ≤ y ≤ xP ( x, y ) = p P ( β = y − x ) = pb y − x , x ≥ y > x. If β has first and second moment finite, with c = 1 − c , c ◦ x ∼ bin( x, c ), as x getslarge m ( x ) = E (( X n +1 − X n ) | X n = x ) = p E ( β ) − qcx ∼ − qcxm ( x ) = E (cid:16) ( X n +1 − X n ) | X n = x (cid:17) = p E (cid:0) β (cid:1) + q (cid:0) cx + c x ( x − (cid:1) ∼ qc x with m ( x ) m ( x ) ∼ − cx , x large.Note the variance of the increment is σ ( X n +1 − X n | X n = x ) = pσ ( β ) + qccx ∼ qccx. Special cases. ( i ) When c = 1, the lower triangular part of P vanishes leadingto P (0 ,
0) = q , P (0 , y ) = p P ( β = y ) = pb y , y ≥ P ( x, y ) = 0, x ≥ ≤ y < x , P ( x, x ) = q , x ≥ P ( x, y ) = p P ( β = y − x ) = pb y − x , x ≥ y > x. The transition matrix P is upper-triangular with diagonal terms. The process X n is non-decreasing, so it drifts to ∞ .( ii ) When c = 0 (total disasters), P (0 ,
0) = q , P (0 , y ) = p P ( β = y ) = pb y , y ≥ P ( x, y ) = 0, x ≥ < y ≤ x , P ( x,
0) = q , x ≥ P ( x, y ) = p P ( β = y − x ) = pb y − x , x ≥ y > x. When a downward move occurs, it takes instantaneously X n to zero (a case of totaldisasters), independently of the value of X n . This means that, defining τ x , =inf ( n ≥ X n = 0 | X = x ), the first extinction time of X n , P ( τ x , = x ) = qp x − , x ≥
1, a geometric distribution with success parameter q , with mean E ( τ x , ) = 1 /q , independently of x ≥ . Note that τ , , as the length of anyexcursion between consecutive visits to 0 , also has a geometric distribution withsuccess parameter q and finite mean 1 /q . In addition, the height H of an excursionis clearly distributed like P τ , − x =1 β x (with the convention P x =1 β x = 0) . Finally,this particular Markov chain clearly is always positive recurrent, whatever the dis-tribution of β . Consecutive excursions are the iid pieces of this random walk onthe non-negative integers.Some Markov catastrophe models involving total disasters are described in [33], [21].( iii ) If β ∼ δ a move up results in the addition of only one individual, which is thesimplest deterministic drift upwards. In this case, the transition matrix P is lower-Hessenberg. This model constitutes a simple discrete version of a semi-stochasticgrowth/collapse model in the continuum, [15]. INOMIAL CATASTROPHES 7 Recurrence versus transience
Using a generating function approach, we start with the transient analysis beforeswitching to the question of equilibrium.3.1.
The transient analysis.
Let B ( z ) = E (cid:0) z β (cid:1) be the pgf of β = β , asan absolutely monotone function on [0 , . Let π ′ n = ( π n (0) , π n (1) , ... ) where π n ( x ) = P ( X n = x ) and ′ denotes the transposition. With z = (cid:0) , z, z , ... (cid:1) ′ , acolumn vector obtained as the transpose ′ of the row vector (cid:0) , z, z , ... (cid:1) , defineΦ n ( z ) = E (cid:0) z X n (cid:1) = π ′ n z , the pgf of X n . With D z =diag (cid:0) , z, z , ... (cid:1) , z = D z , the time evolution π ′ n +1 = π ′ n P yields Φ n +1 ( z ) = π ′ n +1 z = π ′ n P z = π ′ n P D z , leading to the transient dynamics(2) Φ n +1 ( z ) = pB ( z ) Φ n ( z ) + q Φ n (1 − c (1 − z )) , Φ ( z ) = 1 . The fixed point pgf of X ∞ , if it exists, solves(3) Φ ∞ ( z ) = pB ( z ) Φ ∞ ( z ) + q Φ ∞ (1 − c (1 − z )) .( i ) When c = 1, there is no move down possible. The only solution to Φ ∞ ( z ) = pB ( z ) Φ ∞ ( z ) + q Φ ∞ ( z ) is Φ ∞ ( z ) = 0 , corresponding to X ∞ ∼ δ ∞ . Indeed, when c = 1, combined to Φ ∞ (1) = 1 , Φ n +1 ( z ) = ( q + pB ( z )) Φ n ( z ) , Φ ( z ) = 1Φ n ( z ) /n = q + pB ( z ) , showing that, if 1 ≤ ρ := B ′ (1) = E ( β ) < ∞ , n − X n → q + pρ ≥ n → ∞ . The process X n is transient in that, after a finite number of passagesin state 0, it drifts to ∞ .( ii ) When c = 0 (total disasters), combined to Φ ∞ (1) = 1, (3) yields(4) Φ ∞ ( z ) = q − pB ( z ) =: φ ∆ ( z ) , as an admissible pgf solution. In (4), we introduced the integral-valued rv ∆ whosepgf is φ ∆ ( z ) = E (cid:0) z ∆ (cid:1) = q/ (1 − pB ( z )) obtained while compounding a shifted-geometric pgf q/ (1 − pz ) with the pgf B ( z ) of the β s . We conclude that in thetotal disaster setup when c = 0, the law of X ∞ is obtained as a compound shifted-geometric sum ∆ of the β s, whatever the distribution of β . Note that the rv ∆clearly is the height of any total disaster excursion, as the sample path between any A function B is said to be absolutely monotone on (0 ,
1) if it has all its derivatives B ( n ) ( z ) ≥ z ∈ (0 , A geometric( q ) rv with success probability q takes values in N = { , , ... } . A shiftedgeometric( q ) rv with success probability q takes values in N = { , , , ... } . It is obtained whileshifting the former one by one unit. BRANDA GONCALVES AND THIERRY HUILLET two consecutive visits of X n to 0. The length of such excursions clearly is geometricwith success probability q , in that case.3.2. Existence and shape of an invariant pmf ( c ∈ [0 , ). We shall distinguishtwo cases. • The case c ∈ (0 , . From (3) and (4), the limit law pgf Φ ∞ ( z ), if it exists, solves the functional equation(5) Φ ∞ ( z ) = φ ∆ ( z ) Φ ∞ (1 − c (1 − z )) , so that, formally(6) Φ ∞ ( z ) = Y n ≥ φ ∆ (1 − c n (1 − z )) , as an infinite product pgf. Proposition . The invariant measure exists for all c ∈ (0 , if and only if E (cid:0) log + ∆ (cid:1) < ∞ . Proof (Theorem 2 in [25]): By a comparison argument, we need to check theconditions under which π (0) = Φ ∞ (0) converges to a positive number. We getΦ ∞ (0) = Y n ≥ φ ∆ (1 − c n ) > ⇔ X n ≥ (1 − φ ∆ (1 − c n )) < ∞⇔ Z − φ ∆ ( z )1 − z dz < ∞ ⇔ X x ≥ log x P (∆ = x ) = E (cid:0) log + ∆ (cid:1) < ∞ , meaning that ∆ has a finite logarithmic first moment.For most β s therefore, the process X n is positive recurrent, in particular if β hasfinite mean.When β has finite first and second order moments, so do ∆ and X ∞ which exist.Indeed:If B ′ (1) = E β = ρ < ∞ , (with E ∆ = ( pρ ) /q )Φ ′∞ (1) = q (cid:18) pρq + cq Φ ′∞ (1) (cid:19) ⇒ Φ ′∞ (1) = E ( X ∞ ) =: µ = pρq (1 − c ) < ∞ If B ′′ (1) < ∞ , X ∞ has finite variance:Φ ′′ ∞ (1) − Φ ′∞ (1) = pq − c (cid:18) B ′′ (1) + 2 pq c − − c ρ (cid:19) Counter-example:
With β, C >
0, suppose that P ( β > x ) ∼ x ↑∞ C (log x ) − β translating that β has very heavy logarithmic tails (any other than logarithmicslowly varying function would do the job as well). Then E β q = ∞ for all q > β has no moments of arbitrary positive order. Equivalently, B ( z ) ∼ z ↓ − INOMIAL CATASTROPHES 9 C ( − log(1 − z )) β . Therefore, with C ′ = pC/q,φ ∆ ( z ) = q − pB ( z ) ∼ z ↓ − C ′ ( − log (1 − z )) β translating that P (∆ > y ) ∼ y ↑∞ C ′ (log y ) − β shares the same tail behavior as β .From this, P (log ∆ > x ) ∼ j ↑∞ C ′ x − β so that log ∆ has a first moment if and onlyif β > . For such a (logarithmic tail) model of β , we conclude that X remainspositive recurrent if β > β < . The case β = 1is a critical null-recurrent situation.Being strongly attracted to 0, the binomial catastrophe model exhibits a recur-rence/transience transition but only for such very heavy-tailed choices of β . Recallthat:When positive recurrent, the chain visits state 0 infinitely often and the expectedreturn time to 0 has finite mean.When null recurrent, the chain visits state 0 infinitely often but the expected returntime to 0 has infinite mean.When transient, the chain visits state 0 a finite number of times before drifting to ∞ for ever after an infinite number of steps (no finite time explosion is possible fordiscrete-time Markov chains). Corollary.
If the process is null recurrent or transient, no non-trivial ( = ′ )invariant measure exists.Proof: This is because , Φ ∞ ( z ) being an absolutely monotone function on [0 ,
1] ifit exists, π (0) = Φ ∞ (0) = 0 ⇒ Φ ∞ (1 − c ) = 0 ⇒ π ( x ) = 0 for all x ≥ . Clustering (sampling at times when thinning occurs, time change): Let G =inf ( n ≥ B n ( p ) = 0), with P ( G = k ) = p k − q , E (cid:0) z G − (cid:1) = q − pz . G is the timeelapsed between two consecutive catastrophic events. So long as there is no thin-ning of X (a catastrophic event), the process grows of ∆ = P G − k =1 β k individuals.Consider a time-changed process of X whereby one time unit is the time elapsed be-tween consecutive catastrophic events. During this laps of time, the original process X n grew of ∆ individuals, before shrinking to a random amount of its current sizeat catastrophe times. We are thus led to consider the time-changed integral-valuedOrnstein-Uhlenbeck process [ also known as an Integer-Valued Autoregressive ofOrder 1 (in short INAR(1)) process , see [24] ]: X k +1 = c ◦ X k + ∆ k +1 , X = 0 , with ∆ k , k = 1 , , ... an iid sequence of compound shifted geometric rvs . In this form, X k is also a pure-death subcritical branching process with immigration,∆ k +1 being the number of immigrants at generation k +1, independent of X k . With Φ k ( z ) = E (cid:0) z X k (cid:1) , we haveΦ k +1 ( z ) = q − pB ( z ) Φ k (1 − c (1 − z )) , Φ ( z ) = 1 . The limit law (if it exists) Φ ∞ ( z ) also solves (5). ThusΦ ∞ ( z ) = Y n ≥ φ ∆ (1 − c n (1 − z )) , corresponding to X ∞ d = X n ≥ c n ◦ ∆ n +1 . As conventional wisdom suggests, the time-changed process has the same limitlaw as the original binomial catastrophe model, so if and only if the condition E log + ∆ < ∞ holds. Proposition.
When the law of X ∞ exists, it is discrete self-decomposable (SD).Proof: This follows, for example, from Theorem 3 . X k ). See [31] for an account on discrete SD distributions, as aremarkable subclass of compound Poisson ones.The rv X ∞ being SD, it is unimodal, with mode at the origin if π (1) < π (0), orwith two modes at { , } if π (1) π (0) = 1 (see [32], Theorem 4 . P (∆ = 1) = φ ′ ∆ (0) = pq P ( β = 1), we have π (0) = Φ ∞ (0) = q Φ ∞ (1 − c ) π (1) = Φ ′∞ (0) = P (∆ = 1) Φ ∞ (1 − c ) + qc Φ ′∞ (1 − c )= P (∆ = 1) q π (0) + qc Φ ′∞ (1 − c ) > p P ( β = 1) π (0)A condition for unimodality at 0 is thus(7) (log Φ ∞ ) ′ (1 − c ) < − p P ( β = 1) c . Note also π (1) = Φ ′∞ (0) = X m ≥ c m φ ′ ∆ (1 − c m ) Y n = m φ ∆ (1 − c n )= π (0) X m ≥ c m (log φ ∆ ) ′ (1 − c m )= π (0) X m ≥ c m pB ′ (1 − c m )1 − pB (1 − c m )giving a closed-form condition for unimodality at 0. For instance, if B ( z ) = z,π (1) < π (0) if and only if X m ≥ pc m q + pc m < . INOMIAL CATASTROPHES 11
Tails of X ∞ . The probabilities π ( x ) = [ z x ] Φ ∞ ( z ), x ≥ z x − coefficientin the power series expansion of Φ ∞ ( z )) are hard to evaluate. However, someinformation on the large x tails P y>x π ( y ) can be estimated in some cases.- Consider a case where B ( z ) ∼ C · (1 − z/z c ) − as z → z c > P ( β > x ) ∼ C · z − xc has geometric tails. As detailed below, φ ∆ ( z ) ∼ C ′ · (1 − z/z ′ c ) − as z → z ′ c > C ′ = pqC/ (1 − pC ) < P (∆ > x ) ∼ C ′ · z ′− xc also has geometric(heavier) tails but with modified rate z ′ c = z c (1 − pC ) < z c . Then Φ ∞ ( z ) ∼ C ′′ · (1 − z/z ′ c ) − as z → z ′ c > P ( X ∞ > x ) ∼ C ′ · z ′− xc also has geometric tails with rate z ′ c . Indeed, with z ′′ c = ( z ′ c − (1 − c )) /c > z ′ c ,C ′′ · (1 − z/z ′ c ) − ∼ C ′ C ′′ · (1 − z/z ′ c ) − (1 − (1 − c (1 − z )) /z ′ c ) − = C ′ C ′′ (cid:18) z ′ c cz ′′ c (cid:19) · (1 − z/z ′ c ) − (1 − z/z ′′ c ) − showing thatΦ ∞ ( z ) ∼ C ′ C ′′ (cid:18) z ′ c cz ′′ c (cid:19) (1 − z ′ c /z ′′ c ) − · (1 − z/z ′ c ) − , as z → z ′ c > . - Consider a positive recurrent case with B ′ (1) = E β = ∞ . This is the case if,with α ∈ (0 , B ( z ) ∼ − (1 − z ) α as z → , or if (unscaled Sibuya, see [30]) B ( z ) = 1 − (1 − z ) α . In this case φ ∆ ( z ) ∼ − pq (1 − z ) α as a scaled Sibuya rv(with scale factor pq = E ( G − q − p (1 − (1 − z ) α ) = 11 + pq (1 − z ) α ∼ z → − pq (1 − z ) α Then, X n being recurrent, in view of Φ ∞ ( z ) = q − pB ( z ) Φ ∞ (1 − c (1 − z )) , Φ ∞ ( z ) ∼ z → − γ (1 − z ) α where γ = p/ [ q (1 − c ) α ]. Indeed, as z → , − γ (1 − z ) α ∼ (cid:20) − pq (1 − z ) α (cid:21) [1 − γc α (1 − z ) α ]allowing to identify the scale parameter γ. The three rvs β, ∆ and X ∞ have powerlaw tails with index α. • The case c = 0 (total disasters) . In that case, combined to Φ ∞ (1) = 1 , Φ ∞ ( z ) = q − pB ( z ) = φ ∆ ( z )is an admissible pgf solution. The law of X ∞ is a compound shifted-geometric ofthe β s, whatever the distribution of β .The probabilities π ( x ) = [ z x ] Φ ∞ ( z ), x ≥ X ∞ . - If in particular, with C ∈ (0 , B ( z ) ∼ C · (1 − z/z c ) − as z → z c > P ( β > x ) ∼ C · z − xc has geometric tails, thenΦ ∞ ( z ) ∼ C ′ · (1 − z/z ′ c ) − as z → z ′ c > C ′ = pqC/ (1 − pC ) < P ( X ∞ > x ) ∼ C ′ · z ′− xc also has geometric(heavier) tails but with modified rate z ′ c = z c (1 − pC ) < z c . - If B ( z ) = (cid:0) e θz − (cid:1) / (cid:0) e θ − (cid:1) ( β is Poisson conditioned to be positive) is entire.Φ ∞ ( z ) has a simple pole at z c > (cid:0) e θz c − (cid:1) / (cid:0) e θ − (cid:1) = 1 /p and X ∞ has geometric tails with rate z c .- If, with α ∈ (0 , B ( z ) ∼ − (1 − z ) α as z → , or if (unscaled Sibuya) B ( z ) = 1 − (1 − z ) α , then Φ ∞ ( z ) ∼ − pq (1 − z ) α scaled Sibuya (with scale factor pq = E ( G −
1) : q − p (1 − (1 − z ) α ) = 11 + pq (1 − z ) α ∼ z → − pq (1 − z ) α - Suppose, with z > B ( z ) = 1 − (1 − z/z ) α − (1 − /z ) α If there exists z c > B ( z c ) = 1 /p else if B ( z ) = 1 / (1 − (1 − /z ) α ) > /p ((1 − /z ) α > q or z > / (cid:0) − q /α (cid:1) ), then Φ ∞ ( z ) has a simple algebraic poleat z c . If no such z c exists, Φ ∞ ( z ) is entire. This is reminiscent of a condensationphenomenon. We finally observe that1 − φ ∆ ( z )1 − z = p (1 − B ( z ))(1 − z ) (1 − pB ( z )) = pq − B ( z )1 − z φ ∆ ( z ) , P (∆ > n ) = pq n X m =0 P ( β > n − m ) P (∆ = m ) , P (∆ = n ) = pq n − X m =0 P ( β = n − m ) P (∆ = m ) + P (∆ = n ) ! . When is ∆ with φ ∆ ( z ) = q − pB ( z ) itself SD? In any case, ∆ is at least infinitelydivisible (ID) else compound Poisson because φ ∆ ( z ) = exp − r (1 − ψ ( z )) where r > ψ ( z ) is a pgf with ψ (0) = 0. Indeed, with q = e − r , ψ ( z ) = − log (1 − pB ( z )) − log q is a pgf (the one of a Fisher-log-series rv). Proposition.
With b x = [ z x ] B ( z ) , x ≥ , the condition (8) b x +1 b x ≤ x − pb x + 1 for any x ≥ , entails that ∆ is SD. INOMIAL CATASTROPHES 13
Proof:
If ∆ is SD then (see [29], Lemma 2 . φ ∆ ( z ) = e − r R z − h ( z ′ ) − z ′ dz ′ , for some r > h ( z ) obeying h (0) = 0 . We are led to check if pB ′ ( z )1 − pB ( z ) = r − h ( z )1 − z , for some pgf h and r = pb where h ( z ) = 1 − b (1 − z ) B ′ ( z )1 − pB ( z ) = 1 b b (1 − pB ( z )) − (1 − z ) B ′ ( z )1 − pB ( z ) . Denoting the numerator N ( z ), a sufficient condition is that[ z x ] N ( z ) ≥ x ≥ N ( z ) = X x ≥ z x [( x − pb ) b x − ( x + 1) b x +1 ] . (cid:3) Let us show on four examples that these conditions can be met.1. Suppose B ( z ) = z . Then 0 = b x +1 b x ≤ x − px +1 for all x ≥
1. The simple shifted-geometric rv ∆ is SD.2. Suppose B ( z ) = b z + b z with b = 1 − b . We need to check conditionsunder which b b ≤ − pb . This condition is met if and only if the polynomial pb − b + 2 ≤ b ≥ b ∗ where b ∗ ∈ (0 ,
1) is the zero ofthis polynomial in (0 , B ( z ) = αz/ (1 − αz ), α ∈ (0 , α ) rv, with b x = αα x − . The condition reads: α ≤ x − pαx +1 . It is fulfilled if α ≤ x − px + q for all x ≥ α ≤ q/ (1 + q ) < α = b ≥ b ∗ = 1 / (1 + q )).4. Sibuya. Suppose B ( z ) = 1 − (1 − z ) α , α ∈ (0 , b x = α [ α ] x − /x !, x ≥ α ] x = α ( α + 1) ... ( α + x − x ≥ α and[ α ] := 1). The condition reads: α + x − x +1 ≤ x − pαx +1 which is always fulfilled. Theshifted-geometric rv with Sibuya distributed compounding rv is always SD. Proposition . Under the condition that X ∞ is SD and so unimodal, X ∞ has alwaysmode at the origin.Proof: The condition is that, with π (0) = Φ ∞ (0) = q and π (1) = Φ ′∞ (0) = pqB ′ (0) = pqb , π (1) /π (0) = pb < Green (potential) kernel analysis
In this Section, we analyze some important characteristics of the chain making useof its Green (potential) kernel.
Green (potential) kernel at (0 ,
0) :
Contact probability at and firstreturn time to . Suppose X = 0. With then Φ ( z ) = 1, define the double gen-erating function Φ ( u, z ) = P n ≥ u n Φ n ( z ), obeying Φ ( u,
1) = 1 / (1 − u ), see [36].[Note that log z was Laplace-conjugate to X n and now log u is Laplace-conjugateto n ]. Then, 1 u (Φ ( u, z ) −
1) = pB ( z ) Φ ( u, z ) + q Φ ( u, − c (1 − z ))(9) Φ ( u, z ) = 1 + qu Φ ( u, − c (1 − z ))1 − puB ( z )Φ ( u,
0) = 1 + qu Φ ( u, − c ) . With H ( u, z ) := 1 / (1 − puB ( z )) , upon iterating, with Φ (0 , z ) = Φ ( z ) = 1(10) Φ ( u, z ) = X n ≥ ( qu ) n n Y m =0 H ( u, − c m (1 − z )) . In particular(11) Φ ( u,
0) = X n ≥ ( qu ) n n Y m =0 H ( u, − c m ) = 1 + X n ≥ n Y m =1 qu − puB (1 − c m ) . Note that G , ( u ) := Φ ( u,
0) = X n ≥ u n P ( X n = 0) = ( I − uP ) − (0 , ,
0) (the matrix element (0 ,
0) of the resolventof P ). Consequently, Proposition.
The Green kernel G , ( u ) is given by (11).With h m ′ := h u m ′ i Q n − m =0 (1 − puB (1 − c m )) − , we have the following expressionfor the contact probability at u n ] Φ ( u,
0) = Φ n (0) = P ( X n = 0) = n − X m ′ =0 q n − m ′ h m ′ . Remarks. ( i ) Let us have a quick check of this formula. When n = 1 , this leads to P ( X = 0) = q, and, if n = 2 to P ( X = 0) = q + q (cid:2) u (cid:3) (1 − puB (1 − c )) − = q + pqB (1 − c ) . The second part is not quite trivial because it accounts for any movement up in thefirst step (wp p ) immediately followed by a movement down to . This is consis-tent however with the binomial formula P ( X = 0 | X = x ) = q P ( c ◦ x = 0) = q (1 − c (1 − z )) x | z =0 = q (1 − c ) x so that the second part is p X x ≥ P ( β = x ) P ( X = 0 | X = x ) = pq X x ≥ b x (1 − c ) x = pqB (1 − c ) . INOMIAL CATASTROPHES 15 ( ii ) With B m := pB (1 − c m ) , m ≥ , an increasing [0 , − valued sequence con-verging to p , decomposing the product Q n − m =1 (1 − uB m ) − into simple fraction el-ements, we get if n > h m ′ = n − X m =1 A m,n B m ′ m , where A m,n = B n − m / Q m ′ ∈{ ,..,n − }\{ m } ( B m − B m ′ ) . The term A n − B m ′ n − con-tributes most to the sum. Hence, P ( X = 0) = q and if n > P ( X n = 0) = q n n − X m ′ =0 q − m ′ n − X m =1 A m,n B m ′ m = q n n − X m =1 A m,n n − X m ′ =0 (cid:18) B m q (cid:19) m ′ = q n − X m =1 A m,n q n − B nm q − B m is an alternative representation of (12) . In the positive recurrent case, as n → ∞ , P ( X n = 0) → π (0) = Y n ≥ φ ∆ (1 − c n ) > . The Green kernel at (0 ,
0) is thus G , ( u ) = Φ ( u, n ≥ , from the recurrence P ( X n = 0) = P n (0 ,
0) = P nm =0 P ( τ , = m ) P n − m (0 , φ , ( u ) = E ( u τ , ) of the first return time to 0, τ , and G , ( u )are related by the Feller relation (see [3] pp 3 − G , ( u ) = 11 − φ , ( u ) and φ , ( u ) = G , ( u ) − G , ( u ) . Hence,
Proposition.
The pgf φ , ( u ) of the first return time τ , is (13) φ , ( u ) = 1 − G , ( u ) , where G , ( u ) is given by (11). Note G , (1) = X n ≥ P ( X n = 0) = 1 + X n ≥ q n n Y m =0 H (1 , − c m ) = ∞ if and only if X is recurrent, [26], [28]. And in that case, φ , (1) = P ( τ , < ∞ ) =1 − G , (1) = 1 . Positive-(- null) recurrence is when φ ′ , (1) = E ( τ , ) = 1 /π < ∞ (= ∞ ), [22]. Note finally G , (0) = 1 so that φ , (0) = P (cid:0) τ , = 0 (cid:1) = 0.4.2. Starting from x > Green kernel at ( x , and first extinctiontime τ x , . Suppose now X = x > . After shifting X n of x , with Φ ( u, z ) = P n ≥ u n E (cid:0) z x + X n (cid:1) , we now get1 u (Φ ( u, z ) − z x ) = pB ( z ) Φ ( u, z ) + q Φ ( u, − c (1 − z )) . Then(14) Φ ( u, z ) = z x + qu Φ ( u, − c (1 − z ))1 − puB ( z ) , entailingΦ ( u, z ) = X n ≥ ( qu ) n (1 − c n (1 − z )) x n Y m =0 H ( u, − c m (1 − z )) , Φ ( u,
0) = X n ≥ ( qu ) n (1 − c n ) x n Y m =0 H ( u, − c m ) . We obtained:
Proposition.
The contact probability at for the chain started at x > is givenby (15) [ u n ] Φ ( u,
0) = Φ n (0) = P x ( X n = 0) = n − X m ′ =0 (cid:16) − c n − m ′ (cid:17) x q n − m ′ h m ′ . Let us give the first two terms as compared to when x = 0. As required, when n = 1, P x ( X = 0) = q (1 − c ) x and when n = 2, P x ( X = 0) = q (cid:0) − c (cid:1) x + pq (1 − c ) x B (1 − c )a weighted sum of the two terms appearing in the above expression of P ( X = 0) . Corollary. ( i ) When x is large and n fixed, the small but dominant term is when m ′ = 0 whichis q n (1 − c n ) x . So P x ( X n = 0) decays geometrically with x . This expressionquantifies the probability that the population is in an early state of extinction giventhe initial population size was large. Early is when c n ≫ /x (so that (1 − c n ) x ≪ (1 − /x ) x ≪ e − ), so when n ≪ − log c x . ( ii ) In the transient case, when n is large and x is fixed, the dominant term is when m ′ = n − which is q (1 − c ) x ( pB (1 − c )) n . So P x ( X n = 0) decays geometricallywith n at rate pB (1 − c ) . In the positive recurrent case, P x ( X n = 0) → π (0) > as n → ∞ , independently of x . • The Green kernel at ( x ,
0) is thus G x , ( u ) = (cid:2) z (cid:3) Φ ( u, z ) (the matrix element( x ,
0) of the resolvent of P ) . It is related to the pgf of the first extinction time τ x , by the Feller relation(16) φ x , ( u ) = E ( u τ x , ) = G x , ( u ) G , ( u ) . Therefore,
Proposition.
With x > , the pgf of the first extinction time τ x , is (17) φ x , ( u ) = E ( u τ x , ) = P n ≥ ( qu ) n (1 − c n ) x Q nm =1 H ( u, − c m ) P n ≥ ( qu ) n Q nm =0 H ( u, − c m ) . INOMIAL CATASTROPHES 17
In the recurrent case, state 0 is visited infinitely often and so both G , (1) and G x , (1) = ∞ , and P ( τ x , < ∞ ) = φ x , (1) = G x , (1) G , (1) = 1 . We finally note that, because state 0 is reflecting, τ x , is only a local extinctiontime, followed by subsequent extinction times after τ , . In the sequel, we shall let P stand for the substochastic transition matrix obtained from P while deleting itsfirst row and column.There are two alternative representations of φ x , ( u ) . • One alternative representation follows from the following classical first-step anal-ysis:Let X ( x ) be the position of X n started at x .Let X + ( x ) be a positive rv with P ( X + ( x ) = y ) = P ( x , y ) / P y ≥ P ( x , y ), y ≥
1. With τ ′ X + ( x ) , a statistical copy of τ X + ( x ) , , first-step analysis yields, (see[27], [36]): τ x , d = 1 · { X ( x )=0 } + { X ( x ) > } · (cid:16) τ ′ X + ( x ) , (cid:17) . Clearly, P ( X ( x ) = 0) = P ( x ,
0) = q (1 − c ) x =: q x , P ( X ( x ) >
0) =: p x = P y ≥ P ( x , y ) = 1 − q x . Therefore φ x , ( u ) := E ( u τ x , ) , obeys the recurrence φ x , ( u ) = q x u + up x E φ X + ( x ) , ( u ).With φ ( u ) = (cid:0) φ , ( u ) , φ , ( u ) , ... (cid:1) ′ the column-vector of the φ x , ( u ) = E u τ x , ,and q = ( q , q , ... ) ′ the first column-vector of the matrix P, φ ( u ) then solves:(18) φ ( u ) = u q + uP φ ( u ) , whose formal solution is (compare with the explicit expression (17))(19) φ ( u ) = u (cid:0) I − uP (cid:1) − q =: uG ( u ) q , involving the resolvent matrix G ( u ) of P .
Note φ (1) = (cid:0) I − P (cid:1) − q gives thecolumn-vector of the probabilities of eventual extinction φ (1) := (cid:0) φ , (1) , φ , (1) , ... (cid:1) ′ ,so with φ x , (1) = P ( τ x , < ∞ ) if G (1) q < ∞ . Clearly φ (1) = (the all-one col-umn vector) in the recurrent case. In that case, from (18), introducing the columnvector E ( τ ., ) := ( E ( τ , ) , E ( τ , ) , ... ) ′ where E ( τ x , ) = φ ′ x , (1) and observing q + P φ (1) = , we get φ ′ (1) := E ( τ ., ) = + P E ( τ ., ) , equivalently E ( τ ., ) = (cid:0) I − P (cid:1) − = G (1) , or E ( τ x , ) = X y ≥ G x ,y (1) , where G x ,y is the matrix element ( x , y ) of the resolvent of P . • Yet another alternative representation φ ( u ) := (cid:0) φ , ( u ) , φ , ( u ) , ... (cid:1) ′ is as fol-lows. From the identity P n ( x , y ) = P x ( X n = y, τ x , > n ) , we get P ( τ ., > n ) = P n , as the column vector of ( P ( τ , > n ) , P ( τ , > n ) , ... ) ′ ,and so, X n ≥ u n P ( τ ., > n ) = G ( u ) . This leads in particular, as expected, to E ( τ ., ) = G (1) and (compare with (19)and the explicit expression (17)) to(20) φ ( u ) = X n ≥ u n P ( τ ., = n ) = − (1 − u ) G ( u ) . We obtained:
Proposition:
We have φ (1) = P ( τ ., < ∞ ) and so G (1) < ∞ ⇒ P ( τ ., < ∞ ) = , meaning recurrence of X n . In fact, positive recurrence is precisely when E ( τ ., ) = G (1) < ∞ . If G (1) = ∞ , the chain is null-recurrent if (1 − u ) G ( u ) → as u → transient if (1 − u ) G ( u ) → P ( τ ., = ∞ ) as u → , a non-zero limit. The matrix P is substochatic with spectral radius ρ ∈ (0 , . With r and l ′ the cor-responding right and left positive eigenvectors of P , so with P r = ρ r and l ′ P = ρ l ′ , P n ∼ ρ n · rl ′ (as n is large) where rl ′ is the projector onto the first eigenspace. ByPerron-Frobenius theorem, [34], [20], we can normalize l to be of l − norm one to get Proposition:
In the positive recurrent case for X n ( E log + ∆ < ∞ ):( i ) With r ( x ) the x − entry of r , (21) ρ − n P ( τ x , > n ) → r ( x ) , as n → ∞ , showing that P ( τ x , > n ) has geometric tails with rate ρ (extinction is fast) .( ii ) With l ( y ) the y − entry of l , for all x > , (22) P x ( X n = y | τ x , > n ) → l ( y ) , as n → ∞ , showing that the left eigenvector l is the quasi-stationary distribution of X n (orYaglom limit, [37] ), [9] .Proof: In this case, with R = ρ − > , the convergence radius of G , G ( R ) = ∞ and P is R − positive recurrent. ( i ) follows from P ( τ ., > n ) = P n and ( ii ) from P x ( X n = y | τ x , > n ) = P n ( x , y ) /P n . Remark.
The full Green kernel at ( x , y ) is G x ,y ( u ) = [ z y ] Φ ( u, z ) . Hence (23) G x ,y ( u ) = X n ≥ ( qu ) n [ z y ] (1 − c n (1 − z )) x n Y m =0 H ( u, − c m (1 − z ))= X n ≥ ( qu ) n y X y =0 h n,y ( u ) g n,y − y where g n,y = [ z y ] (1 − c n (1 − z )) x = (cid:18) x y (cid:19) c ny (1 − c n ) x − y INOMIAL CATASTROPHES 19 and h n,y ( u ) = [ z y ] n Y m =0 H ( u, − c m (1 − z )) , which can be obtained from a decomposition into simple elements of the inner prod-uct.Using P n ( x , y ) = P nm =1 P ( τ x ,y = m ) P n − m ( y , y ) , n ≥ , we easily get theexpression of the pgf of the first hitting times τ x ,y = inf ( n ≥ X n = y | X = x ) ,as (24) φ x ,y ( z ) = G x ,y ( z ) G y ,y ( z ) . Average position of X n started at x . The double-generating functionΦ ( u, z ) can be used to compute the evolution of E x ( X n ). With Φ ′ ( u, z ) = ∂ z Φ ( u, z ) , Φ ′ ( u,
1) = x X n ≥ ( qcu ) n H ( u, n +1 + pρ − ( p + qc ) (cid:18) − u − − u ( p + qc ) (cid:19) = x H ( u, − qcuH ( u,
1) + pρ − ( p + qc ) (cid:18) − u − − u ( p + qc ) (cid:19) = x − ( p + qc ) u + pρ − ( p + qc ) (cid:18) − u − − u ( p + qc ) (cid:19) This shows that, in the positive recurrent case, with µ n = E x ( X n ) , (25) µ n = [ u n ] Φ ′ ( u,
1) = x ( p + qc ) n + pρq (1 − c ) (1 − (1 − q (1 − c )) n )= x ( p + qc ) n + E ( X n ) → pρq (1 − c ) whatever x , solving(26) µ n +1 = ( p + qc ) µ n + pρ, µ = x . A similar analysis making use of the second derivative Φ ′′ ( u,
1) would yield thetransient evolution of the variance of X n , started at x . We skip the details.4.4.
The height H of an excursion. Assume X = x and consider a version of X n which is absorbed at 0 . Let X n ∧ τ x , stopping X n when it first hits 0 . Let usdefine the scale (or harmonic) function or sequence ϕ of X n as the function whichmakes Y n ≡ ϕ (cid:16) X n ∧ τ x , (cid:17) a martingale. The function ϕ is important because, forall 0 ≤ x < h, with τ x = τ x , ∧ τ x ,h the first hitting time of { , h } starting from x (assuming ϕ (0) ≡ P (cid:0) X τ x = h (cid:1) = P ( τ x ,h < τ x , ) = ϕ ( x ) ϕ ( h ) , resulting from E ϕ (cid:16) X n ∧ τ x (cid:17) = ϕ ( x ) = ϕ ( h ) P ( τ x ,h < τ x , ) + ϕ (0) P ( τ x ,h > τ x , ) . • The case β ∼ δ . Let us consider the height H of an excursion of the originalMC X n first assuming β ∼ δ (a birth event adds only one individual) . Withprobability q , H = 0 and with probability p , starting from X = 1, it is the heightof a path from state 1 to 0 of the absorbed process X n . Using this remark, theevent H = h is realized when τ ,h < τ , and τ h,h +1 > τ h, , the latter two eventsbeing independent. Thus (with P ( H = 0) = q ):(27) P ( H = h ) = p ϕ (1) ϕ ( h ) (cid:18) − ϕ ( h ) ϕ ( h + 1) (cid:19) , h ≥ . We clearly have P h ≥ P ( H = h ) = p because partial sums form a telescopingseries. But (27) is also(28) P ( H ≥ h ) = 1 /ϕ ( h ) , h ≥ , with ϕ (1) = 1 /p .It remains to compute ϕ with ϕ (0) = 0 and P ( H ≥
1) = 1 /ϕ (1) = p . We wish tohave: E x ( Y n +1 | Y n = x ) = x , leading to ϕ ( x ) = pϕ ( x + 1) + q x X y =1 (cid:18) xy (cid:19) c y (1 − c ) x − y ϕ ( y ) , x ≥ . The vector ϕ is the right eigenvector associated to the eigenvalue 1 of the modifiedversion P ∗ of the stochastic matrix P having 0 as an absorbing state: ( P ∗ (0 ,
0) =1), so with: ϕ = P ∗ ϕ , ϕ (0) = 0, [27]. The searched ‘harmonic’ function is increas-ing and given by recurrence, ϕ (1) = 1 /p and(29) ϕ ( x + 1) = 1 p ϕ ( x ) [1 − qc x ] − q x − X y =1 (cid:18) xy (cid:19) c y (1 − c ) x − y ϕ ( y ) ! , x ≥ ϕ (2) = 1 p (1 − qc ) ϕ (1) = 1 p (1 − qc ) ,ϕ (3) = 1 p (cid:0) ϕ (2) (cid:2) − qc (cid:3) − qc (1 − c ) ϕ (1) (cid:1) , = 1 p (1 − qc ) (cid:0) − qc (cid:1) − qp c (1 − c ) . The sequence ϕ ( x ) is diverging when the chain X is recurrent. Proposition.
When β ∼ δ , equations (28) and (29) characterize the law of theexcursion height H of the random walker in the recurrent case. In the transientcase, ϕ ( x ) converges to a value ϕ ∗ and P ( H = ∞ ) = 1 /ϕ ∗ = P ( τ , = ∞ ) . • General β . Whenever the law of β is general, the matrix P ∗ is no longer lowerHessenberg and the harmonic vector ϕ = P ∗ ϕ , with ϕ (0) = 0 , cannot be obtainedby a recurrence. However, the event H ≥ h ≥ β ≥ h or, if β < h , whenever for all states h ′ ≥ h being hitwhen the amplitude β of a last upper jump is larger than h ′ − h , then τ β ,h ′ < τ β , . Hence,
INOMIAL CATASTROPHES 21
Proposition.
For a recurrent walker with general β , P ( H = ∞ ) = 0 where, when h ≥ , P ( H ≥ h ) = p P ( β ≥ h ) + p h − X x =1 P ( β = x ) X h ′ ≥ h ϕ ( x ) ϕ ( h ′ ) P ( β > h ′ − h )= p P ( β ≥ h ) + p · h − X x =1 P ( β = x ) ϕ ( x ) · X h ′′ ≥ ϕ ( h + h ′′ ) P ( β > h ′′ ) generalizing (28). Estimation from an N -sample of X , say ( x , x , ..., x N )We briefly sketch here how (in presence of real data which are suspected to be inthe binomial catastrophe framework), to estimate its constitutive parameters.From the transition matrix P ( x, y ), the log-likelihood function of the N -sample is L ( x , x , ..., x N ) = N X n =1 (cid:2) log (cid:0) qd x n − ,x n (cid:1) ( x n ≤ x n − ) + log (cid:0) pb x n − x n − (cid:1) log ( x n > x n − ) (cid:3) . If one knows that some population grows and decays according to the binomialcatastrophe model with E β = ρ < ∞ , we propose the following estimators: theMaximum-Likelihood-Estimator of p while setting ∂ p L = 0 is(30) b p = 1 N N X n =1 ( x n > x n − ) . With ρ = E β < ∞ , if the law of β is a known one-parameter ρ − family of pmfs,(31) b ρ = 1 P Nn =1 ( x n > x n − ) N X n =1 ( x n − x n − ) ( x n > x n − ) . Note that if c ρp = 1 N N X n =1 ( x n − x n − ) ( x n > x n − ) , then c ρp = b ρ b p. Also, in view of P ( X n +1 = x | X n = x ) = qc x , x ≥ , (32) b c = 1 P Nn =1 ( x n = 1) N X n =1 ( x n = 1 , x n − = 1) . Analysis of the absorbed version of X n and quasi-stationarity We consider again a version of X n started at x > , but which is now absorbedwhen first hitting 0. We aim at having additional insight into the quasi-stationarydistribution, which is known to be a difficult issue. The absorbed process wasconsidered in [19]. We work under the condition that the original (non-absorbed)process is not transient, because if it were, so would its absorbed version with apositive probability to drift to ∞ . The only thing which changes in the transition matrix P as from (1) is its first row which becomes P (0 , y ) = δ ,y , y ≥
0, leading to P ∗ . Letting Φ n ( z ) = E (cid:0) z X n (cid:1) with Φ ( z ) = z x , taking into account the behaviorof X n at 0, with Φ ( z ) = z x ,the recurrence (2) now becomesΦ n +1 ( z ) = Φ n (0) + pB ( z ) (Φ n ( z ) − Φ n (0)) + q [Φ n (1 − c (1 − z )) − Φ n (0)] , Φ n +1 ( z ) = p (1 − B ( z )) Φ n (0) + pB ( z ) Φ n ( z ) + q Φ n (1 − c (1 − z )) , Φ n +1 ( z ) − Φ n ( z ) = p ( B ( z ) −
1) (Φ n ( z ) − Φ n (0)) + q (Φ n (1 − c (1 − z )) − Φ n ( z )) . One can check that for all z ∈ [0 , , Φ ∞ ( z ) = 1 , consistently with X ∞ = 0 (eventual extinction) . Note also Φ n +1 (1) = Φ n (1) = Φ (1) = 1 (no mass loss),Φ n +1 (0) = p Φ n (0) + q Φ n (1 − c ) else,Φ n +1 (0) − Φ n (0) = q (Φ n (1 − c ) − Φ n (0)) > n ≥ , showing that Φ n (0) = P ( X n = 0 | X = x ) is increasing tending to 1 under therecurrence condition. Finally,- Φ n (0) = 1 ⇔ X n = 0 ⇔ Φ n ( z ) = 1 ⇒ Φ n (1 − c ) = 1 ⇒ Φ n +1 (0) = 1 (0 isabsorbing)- X n = 0 ⇔ τ x , ≤ n ⇒ Φ n (0) = P ( τ x , ≤ n ).- P ( τ x , = n ) = P ( τ x , ≤ n ) − P ( τ x , ≤ n −
1) = Φ n (0) − Φ n − (0) > τ x , be the first (and last) hitting time of 0 for X n given X = x >
0. Ithas the same distribution as the one obtained for the original non-absorbed chain.Considering the substochastic transition matrix P of X n where the first row andcolumn have been removed, the dynamics ofΨ n ( z ) = E (cid:0) z X n τ x , >n (cid:1) is(33)Ψ n +1 ( z ) = pB ( z ) Ψ n ( z )+ q [Ψ n (1 − c (1 − z )) − Ψ n (1 − c )] , Ψ ( z ) = z x , Ψ (0) = 0 . Note Ψ n (0) = 0 entails Ψ n +1 (0) = 0 : as required, there is no probability mass at0, because Ψ n ( z ) is the pgf of X n , on the event τ x , > n .Furthermore, Ψ n +1 (1) = p Ψ n (1) + q [Ψ n (1) − Ψ n (1 − c )] = Ψ n (1) − q Ψ n (1 − c ) , so that Ψ n +1 (1) − Ψ n (1) < n (1) = P ( τ x , > n ) (note Ψ n (1) = 1 − Φ n (0)). Conditioning, withΦ cn ( z ) = E (cid:0) z X n | τ x , > n (cid:1) , upon normalizing, we get Proposition.
Under the positive recurrence condition for X n , as n → ∞ (34) Φ cn ( z ) := Ψ n ( z )Ψ n (1) → Ψ ∞ ( z )Ψ ∞ (1) , the pgf of the quasi-stationary distribution l (as a left eigenvector of P ). An explicitexpression of Ψ n ( z ) follows from (35), (36) below. INOMIAL CATASTROPHES 23
The double generating function of Ψ n ( z ) is1 u (Ψ ( u, z ) − z x ) = pB ( z ) Ψ ( u, z ) + q Ψ ( u, − c (1 − z )) − q Ψ ( u, − c )Its iterated version is(35) Ψ ( u, z ) = X n ≥ ( qu ) n (1 − c n (1 − z )) x n Y m =0 H ( u, − c m (1 − z )) − Ψ ( u, − c ) X n ≥ ( qu ) n n − Y m =0 H ( u, − c m (1 − z )) . If z = 1 − c ,Ψ ( u, − c ) = X n ≥ ( qu ) n (cid:0) − c n +1 (cid:1) x n Y m =0 H (cid:0) u, − c m +1 (cid:1) − Ψ ( u, − c ) X n ≥ ( qu ) n n − Y m =0 H (cid:0) u, − c m +1 (cid:1) , so that(36) Ψ ( u, − c ) = P n ≥ ( qu ) n (cid:0) − c n +1 (cid:1) x Q nm =0 H (cid:0) u, − c m +1 (cid:1) P n ≥ ( qu ) n Q nm =1 H ( u, − c m ) . Plugging (36) into (35) yields a closed form expression of Ψ ( u, z ) and then ofΨ n ( z )Ψ n (1) = [ u n ] Ψ ( u, z )[ u n ] Ψ ( u, . The value of Ψ ( u, z ) at z = 1 is Ψ ( u,
1) = P n ≥ u n P ( τ x , > n ) so that Ψ ( u,
1) =(1 − E u τ x , ) / (1 − u ). With H ( u,
1) = 1 / (1 − pu ) , from (35)Ψ ( u,
1) = X n ≥ ( qu ) n n Y m =0 H ( u, − Ψ ( u, − c ) X n ≥ ( qu ) n n − Y m =0 H ( u, − pu − qu/ (1 − pu ) − Ψ ( u, − c ) qu − pu − qu/ (1 − pu )= 11 − u − Ψ ( u, − c ) qu − u = 11 − u (1 − qu Ψ ( u, − c )) . The pgf of τ x , is thus also (a third representation of φ x , ( u ) := E u τ x , , see (19)and (20))(37) φ x , ( u ) = qu Ψ ( u, − c ) , with Ψ ( u, − c ) given by (36). This expression of φ x , ( u ) is consistent with (17). A variant of the binomial catastrophe model
Suppose we are interested in the following simple semi-stochastic decay/surge model:at each step of its evolution, the size of some population either grows by a randomnumber of individuals or shrinks by only one unit. The above binomial catastrophemodel is not able to represent this scenario where at least one individual is removedfrom the population at catastrophic events. To remedy this, we therefore define andstudy a variant of the above binomial model whereby the transition probabilities inthe bulk and at 0 are slightly modified in order to account for the latter decay/surgesituation, [17].- If X n ≥
1, define X n +1 = 1 ◦ X n + β n +1 = X n + β n +1 wp pX n +1 = c ◦ ( X n −
1) wp q Given a move down wp q : One individual out of X n is systematically removedfrom the population ( X n → X n − X n − c / − c ) in thenext generation.- If X n = 0 ( p + q = 1) then, X n +1 = β n +1 wp p , = 0 wp q . Unless p = p , our model yields some additional control on the future of the popu-lation once it hits 0 (extinction event).With b x := P ( β = x ), x ≥ , the one-step-transition matrix P of the modified MC X n now is given by: P (0 ,
0) = q , P (0 , y ) = p δ , if y ≥ ,P ( x, x ) = 0 if x ≥ ,P ( x, y ) = q (cid:18) x − y (cid:19) c y (1 − c ) x − − y , if x ≥ ≤ y < x,P ( x, y ) = p P ( β = y − x ) = pb y − x , if x ≥ y > x. Note that because at least one individual dies out in a shrinking event, the diagonalterms P ( x, x ) of P are now 0 for all x ≥ . Remark . One can introduce a holding probability r x to stay in state x given X n = x , filling up now the diagonal of P . This corresponds to a time changewhile considering a modified transition matrix e P where: p → p x = p (1 − r x ) and q → q x = q (1 − r x ) ( p x + q x + r x = 1 ). So, with ρ : = − r , a column vector withentries ρ x = 1 − r x , P → e P = I + D ρ ( P − I ) . If the invariant measure π obeying π ′ = π ′ P, (the fixed point of π ′ n +1 = π ′ n P )exists, then the one of e P also exists and obeys: e π ′ D ρ = π ′ . Suppose first X = 0. Let E z X n := Φ n ( z ) = Φ n ( z ) + Φ n (0) with Φ ( z ) = 1translating X = 0 . Then, with p + q = 1, INOMIAL CATASTROPHES 25 Φ n +1 ( z ) = ( q + p z ) Φ n (0) + pB ( z ) Φ n ( z ) + q − c (1 − z ) Φ n (1 − c (1 − z )) , Φ ( z ) = 0 , Φ n +1 (0) = q Φ n (0) + q − c Φ n (1 − c ) , Φ (0) = 1 . Note Φ n (0) = P ( X n = 0) and Φ n (0) = 0 for each n ≥ . Thus, if these fixed point quantities existΦ ∞ ( z ) = p ( z −
1) Φ ∞ (0) + pB ( z ) Φ ∞ ( z ) + q − c (1 − z ) Φ ∞ (1 − c (1 − z )) , Φ ∞ (0) = qp (1 − c ) Φ ∞ (1 − c ) . We shall iterate the first fixed point equation which makes sense only when c = 0.With C ( z ) = p ( z − − pB ( z ) and D ( z ) = q − c (1 − z ) 11 − pB ( z ) , we getΦ ∞ ( z ) = Φ ∞ (0) X n ≥ C (1 − c n (1 − z )) n − Y m =0 D (1 − c m (1 − z ))= Φ ∞ (0) C ( z ) + D ( z ) X n ≥ C (1 − c n (1 − z )) n − Y m =1 D (1 − c m (1 − z )) = p Φ ∞ (0) ( z − − pB ( z ) X n ≥ ( cq ) n n Y m =1 − c m (1 − z ) 11 − pB (1 − c m (1 − z )) . Except when c = 1, the term inside the bracket has no pole at z = 1. ThenΦ ∞ (1) = 0 and so, assuming Φ ∞ (0) = 0, Φ ∞ ( z ) = 1 for all z ∈ [0 ,
1] is the onlypossible solution to the first fixed point equation. Recalling from the second onethat Φ ∞ (0) = qp (1 − c ) Φ ∞ (1 − c ), we conclude that Φ ∞ (0) = 0 and so Φ ∞ ( z ) = 0for all z ∈ [0 , . The only solution Φ ∞ ( z ) is the trivial null one.It remains to study the cases c = 1 and c = 0 . • If c = 1. In this case, only a single individual can stepwise be removed fromthe population; the transition matrix P is upper- Hessenberg. This constitutes asimple discrete version of a decay/surge model (some kind of time-reversed versionof the simple growth/collapse model).Φ ∞ ( z ) = p z ( z − z (1 − pB ( z )) − q Φ ∞ (0)= p z ( z − z − p (1 − zB ( z )) Φ ∞ (0) . with Φ ∞ (0) = 0 . Letting B ( z ) = (1 − B ( z )) / (1 − z ) the tail generating functionof β , this is also Φ ∞ ( z ) = p z − p (cid:0) zB ( z ) (cid:1) Φ ∞ (0)Φ ∞ ( z ) = p z − p (cid:0) zB ( z ) (cid:1) ! Φ ∞ (0) . With B (1) = B ′ (1) = E β =: ρ , Φ ∞ (1) = p q − pρ Φ ∞ (0) = − c Φ ∞ (1 − c ) . We conclude :
Proposition (phase transition) . - Subcritical case: Φ ∞ (1) + Φ ∞ (0) = 1 ⇒ Φ ∞ (0) = ( q − pρ ) / ( q − pρ + p ) , well-defined as a probability only if pρ < q. In this case, the chain is positive-recurrent.The term pρ is the average size of a move-up which has to be smaller than theaverage size q of a move-down.- Critical case: If q = pρ , then Φ ∞ (0) = 0 ⇒ Φ ∞ ( z ) = 0 for each z . The chain isnull-recurrent and it has no non-trivial ( = ) invariant measure π . - Supercritical case: If ∞ ≥ pρ > q , the chain is transient at ∞ . Examples: ( i ) In addition to c = 1, assume B ( z ) = αz/ (1 − αz ) (a geometric model for β withsuccess probability α ). Then, if p < qα , Φ ∞ (0) = π (0) = ( qα − p ) / ( α ( q + p ) − p )is a probability andΦ ∞ ( z ) = p z − p (1 + z/ (1 − αz )) Φ ∞ (0) = p z (1 − αz ) q − z ( α + pα ) Φ ∞ (0) . Thus, with ( α < a := ( α + pα ) /q < π ( x ) = [ z x ] Φ ∞ ( z ) = π (0) p q ( a − α ) a x − , x ≥ a .( ii ) If in addition to c = 1, we assume B ( z ) = z , X is reduced to a simplebirth and death chain (random walk) on the non-negative integers, reflected atthe origin. In this case, we get Φ ∞ ( z ) = p z − p (1+ z ) Φ ∞ (0) with ( ρ = 1): Φ ∞ (0) =( q − p ) / ( q − p + p ) . The corresponding chain is positive recurrent if p < /
2, null-recurrent if p = 1 / ∞ when p > / . In the positive-recurrentcase, with π (0) = (1 − p ) / (1 − p + p ), we have π ( x ) = [ z x ] Φ ∞ ( z ) = π (0) p q (cid:18) pq (cid:19) x − , x ≥ , a well-known result, [18].In both examples, whenever the process is positive recurrent, the invariant measurehaving geometric decay at different rates a and p/q < INOMIAL CATASTROPHES 27 • If c = 0 (total disasters)Φ ∞ ( z ) = p ( z −
1) Φ ∞ (0) + pB ( z ) Φ ∞ ( z ) + q Φ ∞ (1)Φ ∞ (0) = qp Φ ∞ (1) , leading to Φ ∞ (0) = q/ ( p + q ) = π (0) andΦ ∞ ( z ) = p z − pB ( z ) Φ ∞ (0) , or(38) Φ ∞ ( z ) = (cid:18) p z − pB ( z ) (cid:19) Φ ∞ (0) . We conclude:
Proposition.
In the total disaster case, the modified chain is always positive re-current with invariant measure having pgf (38), as a shifted compound geometricdistribution.
Examples: ( i ) In addition to c = 0, assume B ( z ) = αz/ (1 − αz ) (a geometric model for β with success probability α ). ThenΦ ∞ ( z ) = (cid:18) p z − pB ( z ) (cid:19) Φ ∞ (0) . ( ii ) In addition to c = 0, assume B ( z ) = z. Then, with π (0) = q/ ( p + q ) , Φ ∞ ( z ) = (cid:18) p z − pz (cid:19) Φ ∞ (0) ,π ( x ) = [ z x ] Φ ∞ ( z ) = π (0) p p x − , x ≥ , a geometric distribution with decay rate p . Acknowledgments:
T. Huillet acknowledges partial support from the “Chaire
Mod´elisation math´ematique et biodiversit´e”.
B. Goncalves and T. Huillet acknowl-edge support from the labex MME-DII Center of Excellence (
Mod`eles math´ematiqueset ´economiques de la dynamique, de l’incertitude et des interactions , ANR-11-LABX-0023-01 project). Finally, this work was also funded by CY Initiative of Ex-cellence (grant “Investissements d’Avenir” ANR- 16-IDEX-0008), Project EcoDepPSI-AAP 2020-0000000013. We warmly thank our Referees for suggesting signifi-cant improvements in the presentation of this manuscript.
References [1] Artalejo, J. R., Economou, A., & Lopez-Herrero M. J. (2007). Evaluating growth measures inan immigration process subject to binomial and geometric catastrophes. Math. Biosci. Eng.,4(4), 573-594.[2] Ben-Ari, I., Roitershtein, A., & Schinazi, R. B. (2019). A random walk with catastrophes.Electron. J. Probab. Volume 24, paper no. 28, 21 pp.[3] Bingham, N.H. (2001). Random walk and fluctuation theory. Stochastic Processes: Theoryand Methods. Handbook of Statistics 19 (ed. C. R. Rao & D. N. Shanbhag) 171-213, Elsevier).[4] Brockwell, P. J., Gani, J. & Resnick, S. I. (1982). Birth, immigration and catastrophe pro-cesses. Advances in Applied Probability, Vol. 14, No. 4, 709-731. [5] Baumann, H., Sandmann, W. (2012). Steady state analysis of level dependent quasi-birth-and-death processes with catastrophes. Computers & Operations Research, Volume 39, Issue2, 413-423.[6] Boudali, O., Economou, A. (2013). The effect of catastrophes on the strategic customerbehavior in queueing systems. Naval Research Logistics, Volume 60, Issue 7.[7] Bouzar, N., Satheesh, S. (2008). Comments on α -decomposability. METRON - InternationalJournal of Statistics, vol. LXVI, n. 2, 243-252.[8] Cairns, B., Pollett, P. K. (2004). Extinction times for a general birth, death and catastropheprocess. J. Appl. Probab., Vol. 41, No. 4, 1211-1218.[9] Collet P., Mart´ınez S. & Mart´ın J. S. Quasi-Stationary Distributions.
Springer, New York,2013.[10] Comtet, L.
Analyse combinatoire.
Tome 1. Presses Universitaires de France, Paris, 1970.[11] Dabrowski, C. (2015). Catastrophic event phenomena in communication networks: A survey.Computer Science Review. Volume 18, 10-45.[12] Dette, H. (1996). On the generating functions of a random walk on the non-negative integers.J. of Appl. Probab., Vol. 33, No. 4, pp. 1033-1052.[13] Economou, A. (2004). The compound Poisson immigration process subject to binomial catas-trophes. Journal of Applied Probability, Volume 41, Issue 2, 508-523.[14] Economou, A., Fakinos, D. (2008). Alternative approaches for the transient analysis ofMarkov chains with catastrophes, Journal of Statistical Theory and Practice, 2:2, 183-197.[15] Eliazar, I., Klafter, J. (2006). Growth-collapse and decay-surge evolutions, and geometricLangevin equations. Physica A, 387, 106-128.[16] Eliazar, I., Klafter, J. (2007). Nonlinear Shot Noise: From aggregate dynamics to maximaldynamics. Euro-Phys. Lett., 78, 40001.[17] Eliazar, I., Klafter, J. (2009). The maximal process of nonlinear shot noise. Physica A, 388,1755-1779.[18] Feller, W.
An introduction to probability theory and its applications,
2, Wiley, New York,1971.[19] Fontes, L. R., Schinazi, R. B. (2019). Metastability of a random walk with catastrophes.Electron. Commun. Probab. Volume 24, paper no. 70, 8 pp.[20] Glynn, P. W., Desai, P. Y. (2018). A probabilistic proof of the Perron-Frobenius Theorem.arXiv:1808.04964.[21] Goncalves, B., Huillet, T. (2020). Scaling features of two special Markov chains involvingtotal disasters.
J. Stat. Phys.
Markov chains.
Cambridge University Press, 1998.[28] Sato, K. I. (1972). Potential operators for Markov processes. Proc. Sixth Berkeley Symp. onMath. Statist. and Prob., Vol. 3 (Univ. of Calif. Press), 193-211.[29] Schreiber, K. (1999). Discrete self-decomposable distributions. Dr. rer. nat. Thesis disserta-tion. Otto-von-Guericke-Universit¨at Magdeburg.[30] Sibuya, M. (1979). Generalized hypergeometric, digamma and trigamma distributions. Annalsof the Institute of Statistical Mathematics, 31, 373-390.[31] Steutel, F.W., van Harn, K. (1979). Discrete analogues of self-decomposability and stability.Annals of Probability, 7, 893-899.[32] Steutel, F. W., van Harn, K.
Infinite divisibility of probability distributions on the real line.
Monographs and Textbooks in Pure and Applied Mathematics, 259. Marcel Dekker, Inc.,New York, 2004.
INOMIAL CATASTROPHES 29 [33] Swift, R. J. (2001). Transient probabilities for a simple birth-death-immigration process underthe influence of total catastrophes, Internat. J. Math. Math. Sci., 25, 689-692.[34] Vere-Jones, D. (1967). Ergodic properties of nonnegative matrices. I. Pacific J. Math. 22,361-386.[35] Wilcox, C., Elderd, B. (2003). The effect of density-dependent catastrophes on populationpersistence time. Journal of Applied Ecology, Vol. 40, No. 5, 859-871.[36] Woess, W.
Denumerable Markov chains. Generating functions, boundary theory, randomwalks on trees.