Fertility Heterogeneity as a Mechanism for Power Law Distributions of Recurrence Times
We study the statistical properties of recurrence times in the self-excited Hawkes conditional Poisson process, the simplest extension of the Poisson process that takes into account how the past events influence the occurrence of future events. Specifically, we analyze the impact of the power law distribution of fertilities with exponent , where the fertility of an event is the number of aftershocks of first generation that it triggers, on the probability distribution function (pdf) of the recurrence times between successive events. The other input of the model is an exponential Omori law quantifying the pdf of waiting times between an event and its first generation aftershocks, whose characteristic time scale is taken as our time unit. At short time scales, we discover two intermediate power law asymptotics, for and for , where is associated with the self-excited cascades of aftershocks. For , we find a constant plateau , while at long times, , has an exponential tail controlled by the arrival rate of exogenous events. These results demonstrate a novel mechanism for the generation of power laws in the distribution of recurrence times, which results from a power law distribution of fertilities in the presence of self-excitation and cascades of triggering.
The statistics of ‘recurrence times’, defined as the random durations of time intervals between two consequent events, is widely used to characterize systems punctuated by short-duration occurrences interspersed between quiet phases. For instance, the statistics of recurrence times between earthquakes is the basis for hazard assessment in seismology. The statistics of recurrence times has recently been the focus of researchers interested in the properties of different natural Zaslavsky91 (); Corral2003 (); Corral2004 (); Davisenetal07 () and social systems Barabasi_Nature05 (); Vasquez_et_al_06 ().
The study of recurrence between earthquakes is perhaps the most advanced quantitatively due to the availability of data and the high involved stakes. The statistics of earthquake recurrence times in large geographic domains have been reported to be characterized by universal intermediate power law asymptotics, both for single homogeneous regions Bak2002 (); Corral2003 () and when averaged over multiple regions Bak2002 (); Corral2004 (). These intermediate power laws, as well as the scaling properties of the distribution of recurrence times, were theoretically explained by the present authors SaiSor2006 (); SaiSor2007 (), using the parsimonious ETAS model of earthquakes triggering Ogata88 (), which is presently the benchmark model in statistical seismology. We recall that the acronym ETAS stands for Epidemic-Type Aftershocks and the ETAS model is an incarnation of the Hawkes self-excited conditional Poisson process Hawkes1 (); Hawkes2 (); Hawkes3 (); Hawkes4 (). Our previous works SaiSor2006 (); SaiSor2007 () have shown that one does not need to invoke new universal laws or fancy scaling in order to explain quantitatively with high accuracy the previously reported scaling laws Bak2002 (); Corral2003 (); Bak2002 (); Corral2004 (). In other words, the finding reported in Bak2002 (); Corral2003 (); Corral2004 () do not contain evidence for any new physics/geophysics laws but constitute just a reformulation of the following laws:
Earthquakes tend to trigger other earthquakes according to the same triggering mechanism, independently of their magnitudes.
The Omori-Utsu law for aftershocks, generalized into the phenomenon of earthquake triggering where all earthquakes are treated on the same footing, states that the rate of events that are triggered by a preceding event that occurred at time decays as
The function can also be interpreted as the probability density function (pdf) of the durations of the waiting time intervals between the reference “mother” event and the triggered events of first generation corresponding to direct triggered by the mother event. The constant describes a characteristic microscopic time scale of the generalized Omori law that ensures regularization at small time and normalization.
Our previous analytical derivations SaiSor2006 (); SaiSor2007 (), found in excellent agreement with empirical data Bak2002 (); Corral2003 (); Corral2004 (), was essentially based on the long-memory of the Omori law (1), for large with . However, it did not take into account the impact of heterogenous fertilities, which come in wildly varying values. Indeed, the number of daughters triggered by an earthquake of a given magnitude grows exponentially with its magnitude. For instance, a magnitude 8-earthquake may have tens of thousands of aftershocks of magnitude larger than 2 while a magnitude 2-earthquake may generate no more than 0.1 earthquake on average of magnitude larger than 2 Helmstteterfertility (). Given the fact that the distribution of magnitudes is itself an exponentially decaying function of magnitudes (called the Gutenberg-Richter law), this translates into a heavy tail distribution of fertilities Saihelmsor05 (), i.e. the distribution of the number of first generation events triggered by a given event has the following power law asymptotic:
Precisely, is the probability that the random number of first generation aftershocks triggered independently by a given mother event is equal to a given integer .
In fact, the main approximation in our previous work SaiSor2006 (); SaiSor2007 () was to consider that, for the estimation of the distribution of recurrence times, it is sufficient to assume that each mother event triggers at most one event, so that the power law (2) is completely irrelevant. This surprising approximation was justified by the focus on the tail of the distribution of recurrence times, for which typically only one event, among the set of events triggered by a given earthquake, does contribute.
The goal of the present paper is to reexamine this approximation and present an exact analysis of the impact of the power law form (2) of the distribution of fertilities on the distribution of the recurrence times. To make the analysis feasible and exact, we consider the case where the Omori law is no more heavy-tailed but has a shorter memory in the form of an exponential distribution, expressed by a suitable choice of time units in the form
In addition to getting exact analytical expressions, our study of the case of an exponential memory kernel (3) is motivated by the fact that, for many applications, this is the default assumption Chavezetal05 (); BauwensHautsch09 (); Eymanetal10 (); Azizetal10 (); Aitsahaliaetal10 (); SalmonTham08 (); Filisor12 (). Therefore, this parameterization has an genuine interest and intrinsic value. This exponential memory function should not be confused with the Poisson model, which has no memory. In contrast, the Hawkes process takes into account the full set of interactions between all past events and the future events, mediated by the influence function given by .
The main result of the present paper is the exact expression for the full distribution of recurrence times that results from all the possible cascades of triggering of events over all generations. We make explicit the substantial dependence of the power law exponents characterizing the distribution on the exponent of the power law tail (2) for the distribution of fertilities and on the branching ratio , defined as the mean number of events of first generation triggered per event:
The paper is organized as follows. Section 2 present the self-excited conditional Hawkes Poisson process and the main exact equations obtained using generating probability functions. Section 3 studies the probability of quiescence within a given fixed time interval and derives different statistical properties of earthquake clusters, such as the mean number of earthquakes in clusters and their mean duration. Section 4 presents the main results concerning the probability density functions (pdf) of the waiting times between successive earthquakes. Section 5 summarizes our main results and concludes. The Appendix makes more specific the analytical model used in our derivations and describes useful statistical properties of first generation aftershocks. While we use the language of seismology and events are named ‘earthquakes’, our results obviously apply to the many natural and social-economic-financial systems in which self-excitation occurs.
Ii Branching model of earthquake triggering
Before discussing the statistical properties of recurrence times, let us develop the statistical description of the random number of earthquakes occurring within the time window . The Hawkes model that we consider assumes that there are exogenous events (called “immigrants” in the literature on branching processes, or “noise events”), occurring spontaneously according to a Poissonian stationary flow statistics. Thus, successive instants of the noise earthquakes belong to the stationary Poisson point process with mean rate . Then, the Hawkes model assumes that any given noise earthquake, occurring at , triggers a total of first generation earthquakes. We assume that the set of the total numbers of first generation aftershocks triggered by noise earthquakes are iid random integers, with the same generating probability function (GPF) . Furthermore, the distribution of waiting times between the time of a given noise earthquake and the occurrences of its aftershocks is assumed to be defined by expression (3). In turn, each of the aftershocks of first generation triggered by a given noise earthquake also triggers independently its own first generation aftershocks, with the same statistical properties (same GPF) as the first generation aftershocks.
Specifically, the Hawkes self-excited conditional Poisson process is defined by the following form of the Intensity function
where the history includes all events that occurred before the present time and the sum in expression (5) runs over all past triggered events. The set of parameters is denoted by the symbol . The term means that there are some external noise (or immigrant, exogenous) sources occurring according to a Poisson process with intensity , which may be a function of time, but all other events can be both triggered by previous events and can themselves trigger their offsprings. This gives rise to the existence of many generations of events. In the sequel, we will consider only the case where is constant.
Under the above definitions and assumptions, one can prove that the GPF of the total number of all earthquakes (including noise and triggered events of all generations) occurring within the time window can be expressed as a product of two terms:
The first term is the GPF of the number of all earthquakes in triggered by noise and triggered earthquakes that occurred up to time . The second term is the GPF of the number of all noise earthquakes that occurred within and of all earthquakes triggered in that window by events also in . The factorization of given by expression (6) simply expresses the independence between the branching processes starting outside and within the time window .
where the functions and satisfy the following nonlinear integral equations:
The symbol represents the convolution operator with respect to the repeating time arguments or . We have introduced the auxiliary function
and the cumulative distribution function (cdf) of the first generation aftershocks instants
The functions and have an intuitive geometric sense, illustrated by figure 1. is the GPF of the random number of aftershocks of all generations triggered till the current time by some noise earthquake (or some aftershock) that occurred at the origin of time . In turn, is the GPF of the random number of all generations aftershocks triggered within the window (for ) by some noise earthquake that occurred at the origin of time. The GPFs and are related by their definition as follows:
Iii Probability of quiescence, mean number of earthquakes and mean duration of clusters
iii.1 Probability of quiescence
For , relation (6) reduces to
is the probability that there are no earthquakes (including the noise earthquakes and their aftershocks of all generations) within the window . can be decomposed as the product of the probability that no noise earthquakes occur within and the probability that no aftershocks occur within that could have been triggered by noise earthquakes and their aftershocks that occurred before and until time .
iii.2 Solution of equation (22) and determination of for the exponential pdf
Using the exponential form of the Omori law (3), it is possible to obtain an exact analytical solution of equation (22), which gives us the possibility to explore in detail the probabilistic properties of recurrence times.
Using expression (15) for leads to
whose solution is given by
We have used the fact that
We can now rewrite probability (19) in the form
Taking into account expression (25) and the equality
the explicit calculation of integral (29) yields
where and is the incomplete beta function
In view of the key role played by the function in the following, it is useful to describe some of its properties. A first result of interest is its limit behavior as the branching ratio tends to . Recall that this limit corresponds to the critical regime of the Hawkes process, separating the subcritical phase and the supercritical phase . For , each noise earthquake has only a finite number of aftershocks. For , there is a non-zero probability that a single noise earthquake may generate an infinite number and infinitely long-lives sequence of aftershocks. The relevant physical regime is thus and the boundary value plays a special role, especially when one remembers that can also be interpreted as the ratio of the total number of triggered events to the total number of events HelmstteterSornette03 (). Hence, when , most of the observed activity is endogenous, i.e., triggered by past activity.
Thus, the limit of as reads
For , it is convenient to choose values of that take the form
so that it is possible to express the function
under the form
The auxiliary function is defined as the sum of the first terms of the Taylor series expansion of the logarithm function with respect to :
For (), we have
For (), we have
For (), we have
iii.3 Mean duration of seismic clusters
iii.3.1 Seismic clusters of all types
Let us consider the random duration of the aftershock cluster triggered by the th noise earthquake that occurred at time . By definition,
where is the occurrence time of the last of its triggered aftershocks over all generations. The mean value of is by definition , where is the pdf of the random cluster durations .
By hypothesis, the instants of the noise earthquakes form a Poissonian point process with mean rate . Moreover, within the Hawkes branching process model, the clusters durations are iid random variables. One can easily show that the probability given by (28), that no aftershocks occur within that could have been triggered by noise earthquakes and their aftershocks that occurred before and until time , take the following value in the limit :
Now, can be obtained as the value of the function given by (31) at :
In particular, in the critical case , we have
iii.3.2 Seismic clusters with at least aftershock
The mean duration of clusters given by expression (43) with (45) includes the contribution of the empty clusters for which the noise earthquake does not trigger any aftershock. It is thus interesting to evaluate another derived quantity , defined as the mean duration of clusters that contain at least one aftershock. To get , we divide by the probability that the number of aftershocks is strictly positive,
where is the probability that the number of aftershocks is positive. Accordingly, one introduce the mean rate of the non-empty clusters equal to
where is the mean rate of noise earthquakes.
Obviously, is given by
where is the probability that a noise earthquake does not trigger any first generation aftershock at all. Using the parameterization defined in the Appendix for the Hawkes model, is given by expression (85), leading to
In particular, in the critical case , the mean duration of the non-empty clusters is given by
As an example, taking and yields a mean duration of the non-empty clusters in the critical case equal to . Recall that the unit time is the characteristic decay time of the Omori law (3).
This allows us to define regimes of low seismicity as characterized by the following inequalities
which means that clusters are well individualized, being separated by comparatively long quiet time intervals.
It is useful to generalize the mean duration of clusters that contain at least one aftershock to the mean durations of the clusters that contain at least aftershocks. We now derive the general equation allowing one to calculate these ’s. Using the total probability formula, one can represent the pdf of clusters durations in the form
where is the probability that a given noise earthquake triggers aftershocks of all generations, and is the conditional pdf of cluster durations under the condition that the number of aftershocks is equal to . Accordingly, the pdf of the durations of the clusters that have or more aftershocks is equal to
The corresponding conditional expectation is equal to
In particular, taking into account that
The dependences of the cluster durations (47), (52) and (59) as functions of the exponent are shown in figure 3. Note the large jump in mean durations of clusters containing at least two aftershocks compared with clusters containing at least one aftershock.
Iv Pdf of recurrence intervals
iv.1 General relations
The knowledge of the exact probability (16) allows one to calculate exactly the pdf of the random waiting times between subsequent earthquakes. Indeed, a general result of the theory of point processes states that
where denotes the mean waiting time between subsequent earthquakes. Therefore, the complementary cumulative distribution function (ccdf) of the random waiting times is equal to
By normalization, , so that
Making explicit with expression (31) yields
Using (62), we have
and finally obtain
Differentiating this last expression (66) with respect to yields the pdf of waiting times between successive earthquakes:
Recall that represents , which is defined by expression (12).
In the following subsections, we analyze in detail expressions (68) and (69) in order to derive the properties of the pdf (67). Expression (68) suggests that it is natural to decompose the analysis of into two discussions, one centered on the term surviving in the limit and the other one. The two next subsections analyze these two terms in turn.
Taking the limit amounts to neglecting the occurrence of any noise earthquake within the window of analysis. As shown in the next subsection, this first case already reveals interesting properties of the pdf , which remain valid in the general case .
where the function is given by expression (69).
The main asymptotics of the function are respectively
This regime corresponds to and thus . Relation (71) thus leads to
The denominator of expression (72) determines a new characteristic time scale
from which one can define a critical value for the branching ratio (4)
such that for and for . We shall refer to the case as the subcritical regime, while is called the near-critical regime.
For (subcritical regime), , and one may replace the asymptotics (72) by the pure power law:
In contrast, in the near-critical case , expression (72) leads to two power asymptotics:
Figure 4 shows the pdf (70) of the waiting times between successive earthquakes for , , (), and in the near-critical case . The two power law asymptotics (78) and the exponential asymptotics (74) for large ’s are clearly visible. Figure 5 is the same as figure 4, except for the value . Although, formally, this value belongs also to the near-critical case (), the intermediate power law asymptotics is barely visible and, for any , the subcritical power asymptotics dominates at short times.
Figure 6 shows log-log plots of the pdf for , , and different values of the branching ratio belonging to the subcritical regime. The subcritical power law asymptotics is dominant and the different pdf’s are similar.
iv.3 General case
We now take into account the occurrence of noise earthquakes within the window of analysis. The main qualitative difference between this case and the previous one is the replacement of the large asymptotics (74) by
The interesting regime of well-defined clusters occurs for , for which the typical waiting time between noise earthquakes is much larger than the characteristic decay time of the Omori law, which is the typical waiting time between a noise earthquake and its aftershocks. In the interval , expression (79) reduces to (74) and this exponential decay can be replaced by the plateau
For , expression (79) simplifies into
These different regimes are illustrated in figure 7, which depicts the pdf of waiting times between successive earthquakes for , , and three values of the positive parameter . One can clearly observe the subcritical and near-critical power law asymptotics (78), as well as the plateau (80) joining the exponential asymptotics (74) and (81). The plateau is especially visible for the smallest values and .
In this paper, we have studied the impact of the power law form (2) of the distribution of fertilities on the distribution of the recurrence times. In contrast with previous studies, we have considered the simple case where the Omori law is not heavy-tailed but is described by an exponential function. Motivated by real applications, this choice allows us to develop an exact analytical treatment using the formalism of generating probability functions.
Our analysis emphasizes the importance of three time scales controlling the different regimes of the probability density function (pdf) of waiting times between successive earthquakes:
the characteristic decay time of the exponential Omori law (3) describing the occurrence of first generation aftershocks, taken as our time unit,
the average waiting time between two successive “noise earthquakes”, which constitute the exogenous sources of the self-excited processes followed by their aftershocks,
the characteristic time defined by expression (75) associated with the self-excited cascades of aftershocks.
In the interesting and relevant regime where earthquake clusters are well defined, namely when the typical waiting time till the first aftershocks (unit time) is smaller than the waiting times between noise earthquakes, we have found that the pdf of recurrence times exhibits several intermediate power law asymptotics:
For , .
For , .
For , .
For , .
In these formulas, is the exponent of the power law distribution of fertilities (2), which is the pdf of the number of first generation aftershocks triggered by a given event of any type. In turn, stands for the branching ratio defined by equation (4), i.e. the average number of daughters of first generation per mother event, and is the mean duration of a cluster that starts with a noise earthquake and ends with its last aftershock over all generation. It is given by expression (43).
Only the first two intermediate asymptotics and at short time scales reflect the influence of the power law distribution of fertilities (2), which is revealed by the remarkable effect of the cascade of triggering over the population of aftershocks of many different generations.
Finally, let us stress the differences between the present investigation and our previous work SaiSor2006 (); SaiSor2007 () on the same problem. In Refs. SaiSor2006 (); SaiSor2007 (), we determined the asymptotic behavior at long times of the distribution of recurrence times, under the approximation that it was sufficient to consider only one aftershock at most per mother event (‘noise earthquake’ in the present terminology). In addition, we considered the standard power law Omori law (1) and not the exponential law (3). By performing a detailed analysis made possible by the use of the exactly tractable exponential Omori law (3), the present paper has thus demonstrated the existence of additional short-time intermediate asymptotics that reveal the distribution of fertilities. This opens the possibility to estimate the exponent of the distribution of cluster sizes from purely dynamic measures of activity.
Appendix: Statistics of first generation aftershocks
Given the power law (2) for the right tail of the pdf of the number of the first generation aftershocks probabilities, we show that the leading relevant terms of the expansion of the GPF in powers of take the form
so that the corresponding auxiliary function (11) takes the form
They depend on the branching ratio defined in (4), the exponent of the power law distribution (2) of the number of first-generation aftershocks. The additional scale parameter satisfies the following inequalities
which ensure the necessary constraints
Rather than deriving the form (82) from (2), it is more convenient to show that the tail of the pdf whose GPF is given by (82) is the power law (2). Given (82) and the definition linking to , namely , we obtain
Using the properties of gamma functions and in particular the well-known equality
we obtain that
Accordingly, expression (86) takes the form
using the asymptotic relation
we finally recover the power law (2).
The case where the exponent in expression (82) of requires a special mention. Indeed, this form describes the special situation in which each noise earthquake (and any aftershock as well) can trigger not more than two first generation aftershocks. Accordingly, there are, in general, only three nonzero probabilities
This special situation arises due to the fact that the GPF has been truncated beyond the quadratic order . It will not be considered further in this paper.
- (1) Aït-Sahalia, Y., J. Cacho-Diaz and R.J.A. Laeven, Modeling financial contagion using mutually exciting jump processes, Working Paper 15850, http://www.nber.org/papers/w15850 (2010).
- (2) Azizpour, S., K. Giesecke and G. Schwenkler, Exploring the Sources of Default Clustering, working paper, Stanford University (2010).
- (3) Bak, P., K. Christensen, L. Danon, and T. Scanlon, Unified Scaling Law for Earthquakes, Phys. Rev. Lett. 88, 178501 (2002).
- (4) Barabási, A.-L., The origin of bursts and heavy tails in humans dynamics, Nature 435, 207-211 (2005).
- (5) Bauwens, L. and N. Hautsch, Modelling Financial High Frequency Data Using Point Processes, Handbook of Financial Time Series, Part 6, 953-979, DOI:10.1007/978-3-540-71297-8_41 (2009).
- (6) Chavez-Demoulin, V., A.C. Davison and A.J. McNeil, Estimating value-at-risk: a point process approach, Quantitative Finance, 5 (2), 227-234 (2005).
- (7) Corral, A. (2003), Local distributions and rate fluctuations in a unified scaling law for earthquakes, Phys. Rev. E 68, 035102(R) (2003).
- (8) Corral, A., Universal local versus unified global scaling laws in the statistics of seismicity, Physica A, 340, 590-597 (2004).
- (9) Davidsen, J., S. Stanchits and G. Dresen, Scaling and universality in rock fracture Phys. Rev. Lett. 98, 125502 (2007).
- (10) Errais, E., K. Giesecke and L.R. Goldberg, Affine Point Processes and Portfolio Credit Risk (June 7, 2010). Available at SSRN: http://ssrn.com/abstract=908045.
- (11) Filimonov, V. and D. Sornette, Quantifying reflexivity in financial markets: towards a prediction of flash crashes, Phys. Rev. E 85 (5): 056108 (2012).
- (12) Hawkes, A.G., Point spectra of some mututally exciting point processes. Journal of Royal Statistical Society, series B, 33, 438-443 (1971).
- (13) Hawkes, A.G., Spectra of some mutually exciting point processes with associated variables, In Stochastic Point Processes, ed. P.A.W. Lewis, Wiley, 261-271 (1972).
- (14) Hawkes, A.G. and Adamopoulos, L., Cluster models for earthquakes - regional comparisons, Bull Internat. Stat. Inst. 45, 454-461 (1973).
- (15) Hawkes, A.G. and Oakes D., A cluster representation of a self-exciting process, Journal Apl. Prob. 11, 493-503 (1974).
- (16) Helmstetter, A., Is earthquake triggering driven by small earthquakes? Phys. Res. Lett. 91, 058501 (2003).
- (17) Helmstetter, A. and D. Sornette, Importance of direct and indirect triggered seismicity in the ETAS model of seismicity, Geophys. Res. Lett. 30 (11) doi:10.1029/2003GL017670 (2003).
- (18) Ogata, Y., Statistical models for earthquake occurrence and residual analysis for point processes, J. Am. stat. Assoc. 83, 9-27 (1988).
- (19) Saichev, A., A. Helmstetter and D. Sornette, Anomalous Scaling of Offspring and Generation Numbers in Branching Processes, Pure and Applied Geophysics 162, 1113-1134 (2005).
- (20) Saichev A. and Sornette D., Universal distribution of Interearthquake Times Explained, Physical Review Letters, 97, 078501 (2006).
- (21) Saichev A. and Sornette D., Theory of Earthquake Recurrence Times, Journal of Geophysical Research, 112, B04313-1-26 (2007).
- (22) Saichev, A. and Sornette D., Power law distribution of seismic rates: theory and data analysis, Eur. Phys. J. B 49, 377-401 (2006).
- (23) Salmon, M. and W. W. Tham, Preferred Habitat, Time Deformation and the Yield Curve , working paper (Revised and resubmit - Journal of Financial Markets) (2008).
- (24) Vazquez, A., J. G. Oliveira, Z. Dezso, K. I. Goh, I. Kondor, and A. L. Barabasi, Modeling bursts and heavy tails in human dynamics, Physical Review E 73, 036127 (2006).
- (25) Zaslavsky, G.M. and M.K. Tippett, Connection between recurrence-time statistics and anomalous transport, Phys. Rev. Lett. 67, 3251-3254 (1991).