Theme V – Models and Techniques forAnalyzing Seismicity
Basic models of seismicity: Temporal models

Jiancang Zhuang1 • David Harte2 • Maximilian J. Werner3 •Sebastian Hainzl4 • Shiyong Zhou5

1.
2.
3.
4.
5.

Institute of Statistical Mathematics
Statistics Research Associates
Department of Geosciences, Princeton University
Institute of Geosciences, University of Potsdam
Department of Geophysics, Peking University

How to cite this article:
Zhuang, J., D. Harte, M.J. Werner, S. Hainzl, and S. Zhou (2012), Basic models ofseismicity: temporal models, Community Online Resource for Statistical Seismicity
Analysis, doi:10.5078/corssa-79905851. Available at http://www.corssa.org.
Document Information:
Issue date: 14 August 2012 Version: 1.0

2

www.corssa.org

Contents
1
2
3
4
5
6
7

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Basic mathematical concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Long-term temporal models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Temporal clustering models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Transformed time sequence and residual analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Related software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
4
8
19
33
35
35

Basic models of seismicity: Temporal models

3

Abstract In this and subsequent articles, we present an overview of some models of seismicity that have been developed to describe, analyze and forecast theprobabilities of earthquake occurrences. The models that we focus on are not only
instrumental in the understanding of seismicity patterns, but also important tools
for time-independent and time-dependent seismic hazard analysis. We intend toprovide a general and probabilistic framework for the occurrence of earthquakes.
In this article, we begin with a survey of simple, one-dimensional temporal models
such as the Poisson and renewal models. Despite their simplicity, they remain highly
relevant to studies of the recurrence of large earthquakes on individual faults, to thedebate about the existence of seismic gaps, and also to probabilistic seismic hazard
analysis. We then continue with more general temporal occurrence models such as
the stress-release model, the Omori-Utsu formula, and the ETAS (Epidemic TypeAftershock Sequence) model.

1 Introduction
In this article, we present an overview of several temporal point-process models thatwere well studied during the 1990’s and 2000’s. We cannot go through all the existing interesting models, and we only focus on some generic models, since most of
the complicated models are constructed on the basis of these elegant generic models. These models include the Poisson model, the renewal models, the Omori-Utsu
formula, the ETAS (Epidemic Type Aftershock Sequence) model, and the stressrelease model. To address them systematically, we classify these models into thefollowing four classes: (1) 1D/temporal models, which include the stationary andnon-stationary Poisson model, the renewal model, the Omori-Utsu (ReseanbergJones) model, temporal ETAS models and single-region stress release models; (2)discrete multi-region models, including linked stress release models and self- andmutually exciting models; (3) spatiotemporal models, including homogeneous andnon-homogeneous Poisson models, space-time ETAS models; and, (4) models with
auxiliary information, which incorporate information from observation or/and calculation results of other physical variables. In this article, we only discuss models
from Class 1, i.e., temporal models that are formulated only based on the information of the catalogs. Please see other papers in this series for basic models fromthe 2nd and 3rd classes, and the CORSSA articles by Hainzl et al and Iwata (thistheme) for the 4th class of models.
To illustrate how to utilize these models for seismicity analysis, we give examples
that are extracted from various published papers. We encourage the reader to repeat
these examples and to apply similar analysis to their own datasets.

4

www.corssa.org

2 Basic mathematical concepts
Earthquakes are caused by the rapid energy release accompanying the movement
or fracture of faults in the crust or upper mantle of the earth, causing groundshaking of a duration from seconds to several minutes. However, generally only pointinformation is available for each earthquake, namely the origin time, hypocenter andmagnitude. If we compress the time axis of the seismograph to view all the recordsof earthquakes in a large time scale, the waveforms of each earthquake become apulse. We treat earthquakes as points in time and space. Such idealization enables
us to study the earthquake process by using point process models.
Mathematically, a point process is a random object that takes value of a countable
subset in a given spatiotemporal space. Here we can simply regard a point process asa type of stochastic model that defines probabilistic rules for the occurrence of points(i.e., earthquakes) in time and/or space. As well as purely temporal point processmodels, we consider in this article marked point processes that include earthquake
magnitudes. In a marked point process, a mark or size (e.g., magnitude) is assigned
to each point.
Denote a point process in time by N and a certain temporal location by t. We
assume in the following discussions that we have known observations up to timet. The most important characteristic is the waiting time, say u, to the next eventfrom time t. We consider the following cumulative probability distribution functionof waiting time
Ft (u) = Pr{from t onwards, waiting time to next event ≤ u },(1)

with corresponding probability density function (p.d.f.)ft (u) du = Pr{from t onward, waiting time is between u and u + du },(2)

the survival function
St (u) = Pr{from t on, waiting time > u }= Pr{no event occurs between t and t + u},(3)

and the hazard function
ht (u) du = Pr{next event occurs between t + u and t + u + du
| it does not occur between t and t + u.}.

(4)

Basic models of seismicity: Temporal models

5

One can easily prove that these functions are related in the following ways:

 Z u
Z ∞
ht (s) ds ,ft (s) ds = exp −St (u) = 1 − Ft (u) =
0
u

 Z u
dFt
dSt
ft (u) =ht (s) ds ,=−= ht (u) exp −dudu0ft (u)ddht (u) == − [log St (u)] = − [log(1 − Ft (u))].
St (u)
du
du

(5)
(6)
(7)

In the above equations, the first and second equality signs in each are obvious.
Here we only outline how to derive the third equality signs in each. Equation (7)obtained from (4), i.e.,Pr{next event occurs between t + u and t + u + du }Pr{the waiting time is greater than u.}
f (u) du== −d log S(u).
S(u)

ht (u) du =To obtain the exponential term in (5), replace u by x in the above equality andintegrate both sides with respect to x from 0 to u to get
Z u
Z u
ht (x) dx = −d log S(x) = log S(0) − log S(u) = − log S(u),
(8)
0

0

i.e.,

 Z u
ht (x) dx .
S(u) = exp −(9)0

The last equality sign in (6) can be proved by taking the derivative of the above
equation with respect to u.
In the above concepts, the hazard function ht is the key to understand our theory
because it has the nice property of additivity. Suppose that a point process N
consists of two independent sub-process, N1 and N2 . Using the notation in (1) to(7) for N and denoting the waiting time distribution function, the survival function
(i)
(i)
(i)
and the hazard function of the sub-processes by Ft , St and ht , respectively, withi = 1 or 2, respectively, for N1 or N2 . Then, we can prove
(1)

(2)

(1)

(2)

(1)

(2)

Ft (u) = Ft (u) + Ft (u) − Ft (u) Ft (u),
(1)

(2)

(10)
(1)

(2)

ft (u) = ft (u) + ft (u) − ft (u) Ft (u) − Ft (u) ft (u),

(11)

(1)
(2)
ft (u)ft (u)ft (u)(1)(2)=
ht (u) =+= ht (u) + ht (u).
(1)(2)
1 − Ft (u)1 − Ft (u) 1 − Ft (u)

(12)

6

www.corssa.org

Another property of the hazard function is the following: if there is no event
occurring between t and t + u and v ≥ u ≥ 0, then ht (v) = ht+u (v − u). This isfollowed by conditioning on no event occuring in (t, t + u], then, when x ≥ 0.
St+u (x) = Pr{no event occurs in (t + u, t + u + x) | no event occurs in (t, t + u]}
Pr{no event occurs in either (t + u, t + u + x) or (t, t + u]}=
Pr{no event occurs in (t, t + u]}
Pr{no event occurs in (t, t + u + x)}=
Pr{no event occurs in (t, t + u]}
St (u + x)=
(13)
St (u)
with x ≥ 0 and so
ht+u (x) =ddlog St+u (x) =log St (u + x) = ht (u + x).
dx
dx

(14)

If we set x = 0 in the above equation, we have ht (u) = ht+u (0) when no events occurin (t, t + u], i.e., h collapses to a function of one variable. Thus, as a function of t,
ht (0) is identical to ht0 (t − t0 ) from a time t0 , until the occurrence of another event,say t1 , where ht (0) is identical to ht1 (t − t1 ).
If an event occurs at t, then ht (0) has two possible explanations: one is the hazard
function at t − tprev from the occurrence time tprev of the last previous event, theother is the hazard function at 0 from the event occurring at t. To distinguish them,
we use ht− (0) and ht+ (0) for the first and the second explanations, respectively.
Conventionally, the first explanation is used for ht (0), i.e., ht (0) = ht− (0).
Example 1 If the distribution time from t is an exponential distribution, i.e.,
Ft (u) = 1 − exp[−λu], and ft = λ exp[−λu],(15)where λ is a positive constant, thenht (u) = λ.

(16)This is the case of the Poisson process, which will be discussed in detail in next
section.
Since the function ht (0) plays a key role in our understanding of the theory ofpoint processes, it is worthwhile to give it a specific name. Nowadays, it is usually
named as the conditional intensity and denoted by λ(t). Researchers also frequently
write it as λ(t|Ht ), where Ht represents the observation history up to time t but notincluding t, since the value of λ(t) also can depend on what is in the observation
history before t, but not on what happens at t.

Basic models of seismicity: Temporal models

7

From the definition of the hazard function, we haveλ(t) dt = Pr{one or more events occur in [t, t + dt) | Ht }.

(17)

This is usually used as the definition of the conditional intensity. When the point
process N is simple, i.e., there is at most 1 event occurring at the same time location,then Pr{N [t, t + dt) > 1} = o(dt). Using this relationship, λ(t) can be also defined
as follows:∞X
E [N [t, t + dt) | Ht ] =n Pr{N [t, t + dt) = n | Ht }n=0

≈ Pr{N [t, t + dt) = 1 | Ht } + o(dt)≈ λ(t) dt + o(dt).

(18)

The conditional intensity has a natural form for the purpose of forecasting.
Given an observed dataset of a point process N , say {t1 , t2 , · · · , tn }, in a given
time interval [S, T ]. The likelihood function L is the joint probability density of
waiting times to each of these events, i.e.,the joint probability that the waiting timefrom S is in the range (t1 − S, t1 − S + dt1 ) and the waiting time from t1 is in therange (t2 −t1 , t2 −t1 +dt2 ) and so on, with the final term being that the waiting timefrom tn is greater than (T − tn ). The probability each of these terms is conditionalon what has happened previously, so we can write the likelihood function as
L(N ; S, T ) dt1 dt2 · · · dtn= Pr{The waiting time from S is in (t1 − S, t1 − S + dt1 )}× Pr{The waiting time from t1 is in (t2 − t1 , t2 − t1 + dt2 ) | what happens before t1 }×···× Pr{The waiting time from tn is greater than T − tn | what happens before tn }= fS (t1 ) dt1 × ft1 (t2 − t1 ) dt2 × · · · × ftn−1 (tn − tn−1 ) dtn × Stn (T − tn )
 Z t1hS (u − S) du dt1= hS (t1 − S) exp −
S

×
=

n−1
Y

 Z ti+1


 Z T

hti (ti+1 − ti ) exp −hti (u − ti ) du dti exp −htn (u − tn ) du

" ni=1Yhti (0) dti

i=1

ti

#

(∵ (6))

tn

#
 Z T
 "Y
 Z T

n
exp −
hu (0) du =λ(ti ) dti exp −λ(u) du .
S

We usually write the above formula as its logarithm, i.e.,
Z TnXlog L(N ; S, T ) =log λ(ti ) −λ(u) du.
i=1

(19)

S

i=1

S

(20)

8

www.corssa.org

A direct use of the likelihood function is for parameter estimation: when the model
involves some unknown regular parameters, say θ, we can estimate θ through maximizing the likelihood, i.e., the MLE (maximum likelihood estimate) isθ̂ = argθ max L(N ; S, T ; θ).

(21)

Given several models fit to the same dataset, the optimal one can be selected byusing the Akaike Information Criterion (AIC, see Akaike 1974). The statistic
AIC = −2 max log L(θ) + 2kp
θ(22)

is computed for each of the models fit to the data, where kp is the total number of
estimated parameters for a given model. Comparing models with different numbersof parameters, the addition of the quantity 2kp roughly compensates for the additional flexibility that the extra parameters provide. The model with the lowest AIC
value is selected as the optimal choice for forward prediction purposes.
3 Long-term temporal models
In this section, we present the one-dimensional temporal models for long-term earthquake forecasting, including Poisson models, renewal models and stress release models. Roughly speaking, these models are intended to be used at time scales of 10 to
100 years or even longer.
3.1 Stationary Poisson models
We suppose that the readers of this section already have some related backgroundknowledge of the Poisson process; if not, the readers are referred to Introduction toProbability Models (Ross 2003).
The Poisson process is named after the French mathematician Siméon-Denis
Poisson (1871-1840). The process is also known as the model of complete randomness. Thus, it naturally serves as the null model in many hypothesis tests to clarify
whether systematic structure is contained in the observations. For example, Gardner and Knopoff (1974) discussed whether the seismicity in California is Poissonian
or not after the aftershock clusters are removed. Furthermore, the stationary Poisson process forms the basis of most of today’s probabilistic seismic hazard analysis(Cornell 1968).
Let N (a, b) be the number of events in a point process N that fall in the timeinterval (a, b). A point process is called a stationary Poisson process if the system
has the following characteristics:

Basic models of seismicity: Temporal models

9a) Independent increments: The numbers of events occurring in two disjoint timeintervals are independent of each other. That is,
Pr{N (a, b) = m, N (c, d) = n} = Pr{N (a, b) = m} Pr{N (c, d) = n}(23)

for all non-negative integers m, n, and any interval pairs (a, b) ∩ (c, d) = ∅.
b) Stationarity: The probability distribution of the number of events falling in atime interval only depends on the length of the time interval.
c) Simplicity: There is never an occasion that two or more events occur simultaneously.
It can be shown that a stationary Poisson process has the following properties:a) The number of events occurring in an interval of length S has a Poisson distribution, i.e.,
λn S n −λS
Pr{N (a, a + S) = n} =
e ,
(24)
n!
where λ is the long term average of the number of events in an interval of unit
length.
b) Given a fixed time, the forward occurrence time tf (time to next event) and the
backward occurrence time tb (time back to the previous event) have the same
exponential distribution as the waiting time τ between two consecutive events,i.e.,Pr{tf < x} = Pr{tb < x} = Pr{τ < x} = 1 − e−λx .
(25)c) The mean and the variance of the number of observations are the same:
E [N (a, b)] = λ(b − a),Var [N (a, b)] = λ(b − a).

d) Binomial distribution. Given that n is the total number of events occurringin the interval [0, T ], then the number of events in [0, S], where S < T , is a
non-negative integer-valued random variable from a binomial distribution withparameters (n, S/T ), i.e.,   k 
n−k
S
S
n
Pr{N [0, S] = k | N [0, T ] = n} =
1−
(26)
k
T
T
e) Uniform distribution: Given that at least one event occurs in [0, T ], then thelocation of an arbitrary event is uniformly distributed in [0, T ].
f) Beta distribution: Given that there are a total of n events occurring in [0, 1],n ≥ 1, the location t(k) of the kth event (k < n), from the left to the right, has abeta distribution, i.e.,dx
Pr{t(k) ∈ (x, x + dx) | N [0, 1] = n} =xk−1 (1 − x)n−k−1 ,(27)
B(k, n − k)R1
where B(x, y) = 0 tx−1 (1 − t)y−1 dt is the beta function.

10

www.corssa.org

Conditional intensity Recall the definition of the conditional intensity function of apoint process, which is defined by
λ(t)dt = Pr {N [t, t + dt) = 1 | Ht } .

(28)

For the stationary Poisson processes with rate λ, the right-hand side is λ dt exp[−dt]+o(dt), i.e., λ(t) = λ. This is why we use λ to denote both the rate of the Poissonprocess and the conditional intensity in the case of more general point processes.
Likelihood functions and maximum likelihood estimates of the Poisson process Suppose that all observations of a realization of a Poisson process with a rate of λ ona time interval [S, T ] is {t1 , t2 , · · · , tn }, where S ≤ t1 < t2 < · · · < tn < T . The
likelihood of the observation is
L[S, T ; λ] = λn e−λ(T −S) ,(29)

and the logarithm of the likelihood function islog L(S, T ; λ) = n log λ − λ(T − S).

(30)

Usually, the rate λ is unknown. We can estimate it through maximizing the likelihood function, i.e., the maximum likelihood estimate (MLE) of λ isλ̂ = argλ max L(S, T ; λ)(31)

Taking the derivative of (30) and setting it to zero to obtain a maximum, we have∂log L(S, T ; λ)
= 0,
∂λ
λ=λ̂

(32)

λ̂ = n/(T − S),(33)

which yields
which is simply the average rate of events in the interval [S, T ].
3.2 Non-stationary Poisson models
A point process is a non-stationary Poisson process if it satisfies the independent
increments and simplicity conditions, but violates the stationarity condition. Consequently, the Poisson process has a rate that is a function of time, denoted by
λ(t).
Rt
R t+s
With the definition Λ(t) = 0 λ(x) dx and Λ(t+s) = 0 λ(x) dx, the probabilityof observing n events in the interval [t, t + s) is
Pr{N [t, t + s) = n} =e−(Λ(t+s)−Λ(t)) [Λ(t + s) − Λ(t)]n.
n!

Basic models of seismicity: Temporal models

11

R t+s
That is, N [t, t + s) is Poisson distributed with expectation value t λ(u) du where
λ(t) is the time dependent intensity. Also, given an event that falls in [S, T ], the probRT
ability density of its location is proportional to λ(·) on [S, T ], i.e., λ(t)/ S λ(u) du.
Likelihoods and maximum likelihood estimates of non-stationary Poisson processes
Suppose that the observation of a realization of a non-stationary Poisson processwith a rate of λ(t) on a time interval [S, T ] is {t1 , t2 , · · · , tn }, where S ≤ t1 < t2 <· · · < tn < T . In a similar way as for the stationary Poisson model, the likelihoodof the observations isYRTλ(ti ).
(34)
L(S, T ; λ) = e− S λ(u) dui

The corresponding logarithm of the likelihood function isZ TXλ(u) du.
log L(S, T ; λ) =log λ(ti ) −
i

(35)

S

There are relative few applications of non-stationary Poisson models in seismology. This type of model is usually used as a natural alternative to the null stationary
Poisson model when an apparent trend or seasonality is visible in the data. Some examples of non-stationary Poisson models are the exponential trend model in Zheng
and Vere-Jones (1994) and the polynomial trend models and the seasonality models
in Ogata and Katsura (2004) and Ma and Vere-Jones (1997).
3.3 Renewal/recurrence models
Renewal models are simple extensions of the Poisson model. One of their uses is
to model the characteristic recurrence of earthquakes on a particular fault or ina particular region. This class of models is widely used in seismicity and seismic
hazard analysis. For example, Field (2007) summarized how the Working Group on
California Earthquake Probabilities (WGCEP) estimates the recurrence probabilities of large earthquakes on major fault segments using various recurrence models
to produce the official California seismic hazard map. These models are sometimes
justified by the elastic rebound theory proposed by Reid (1910). According to his
theory, large earthquakes release the elastic strain that has built up since the lastlarge earthquake. Some seismologists deduce that the longer it has been since thelast earthquake, the more probable is an imminent event (e.g. Nishenko and Buland
1987; Nishenko 1991; Sykes and Menke 2006), while others contend that the data
contradict this view (e.g. Davis et al. 1989; Kagan and Jackson 1995). Renewal
models are often used to quantitatively demonstrate whether earthquakes occur
temporally in clusters or quasi-periodically.

12www.corssa.org

A renewal process is a generalization of the Poisson process. In essence, thePoisson process has independent identically distributed waiting times that are exponentially distributed before the occurrence of the next event. A renewal process
is defined as a point process with the waiting times having a more general distribution, which is not necessarily exponential distribution. In this section, we denotethe density function of the waiting times by f (·), which is also usually called therenewal density.
The conditional intensity of the renewal process can be derived as follows. Recall
that the time interval between any two adjacent events is independent from otherintervals. This means that the occurrence of the next event depends only on the timeof the last event but not on the full history. Thus, from Section 1, the conditional
intensity function is the same as the hazard function, i.e.,λ(t) =f (t − tN (t− ) ),1 − F (t − tN (t− ) )(36)

where tN (t− ) is the occurrence time of the last event before t and F is the cumulative
probability function of f .
The following probability functions are often chosen as the the renewal densities.
1. The gamma densitye−u/θ
f (u; k, θ) = uk−1 k,θ Γ (k)for u ≥ 0 and k, θ > 0 ,(37)with a hazard functionuuk−1 e− θ,h(u; k, θ) = kθ Γ uθ (k, uθ )(38)

where θ is called the scale parameter, k the shape parameter, and Γ and Γα
are the Rgamma and the incomplete
functions, respectively, defined by
R ∞gamma∞ x−1 −tx−1 −tΓ (x) = 0 t e dt and Γα (x) = α t e dt.
2. The log-normal density function isf (u; µ, σ) =

1
√

uσ 2π

e−

(ln(u)−µ)22σ 2, for u ≥ 0,(39)

with a hazard functionh(u; µ, σ) =

2
√

e−

(ln(u)−µ)22σ 2

uσ 2π 1 − erfhln(u)−µ
√
σ 2

i,

(40)

Basic models of seismicity: Temporal models

13

where µ and σ are the mean and standard deviation of the variable’s natural
logarithm1 . and erf is the error function.
3. The Weibull distribution is one of the most widely used lifetime distributions
in reliability engineering. It is named after Waloddi Weibull who described it
in detail in 1951 (Weibull 1951). The probability density function of a Weibull
random variable X is  
 k u k−1 −(u/µ)keu ≥ 0,
f (u; µ, k) = µ µ(41)0u < 0,
where k > 0 is the shape parameter and µ > 0 is the scale parameter of the distribution. Its cumulative distribution function is a stretched exponential function
F (u; µ, k) = 1 − e−(u/µ)

k

and its hazard function iskh(u; k, µ, k) =µ

 k−1u.
µ4. Kagan and Knopoff (1987) used the inverse Gaussian distribution to model theevolution of stress as a random walk with a tectonic drift. Later on, a particularly interesting renewal model, called the Brownian passage-time model, wasintroduced by Matthews et al. (2002) and Ellsworth et al. (1999) based on theproperties of the Brownian relaxation oscillator (BRO). In this conceptual modelof Matthews et al. (2002), the loading of the system has two components: (1) aconstant-rate loading component, λt, and (2) a random component, (t) = σW (t),that is defined as a Brownian motion (where W is a standard Brownian motionand σ is a nonnegative scale parameter). The Brownian perturbation process forthe state variable X(t) (see Figure 2 in Matthews et al. (2002)) is defined as:
X(t) = λt + σW (t).
An event will occur when X(t) ≥ Xf ; event times are seen as “first passage”
or “hitting” times of Brownian motion with drift. The BRO are a family of
stochastic renewal processes defined by four parameters: the drift or mean loading(λ), the perturbation rate (σ 2 ), the ground state (X0 ), and the failure state(Xf ). On the other hand, the recurrence properties of the BRO (response times)
are described by a Brownian passage-time distribution which is characterizedby two parameters: (1) the mean time or period between events, (µ), and (2)1 If X is a normal random variable with mean µ and variance σ 2 , then Y = exp X has a log-normal distribution
with a density of (39).

14

www.corssa.org

the aperiodicity of the mean time, α, which is equivalent to the coefficient ofvariation (defined in as the ratio of the variance to the mean occurrence time).
The probability density for the Brownian passage-time (BPT) model is given by: µ  12(t − µ)2
exp −
, t≥0
(42)
f (t; µ, α) =2πα2 t3
2α2 µtwith a cumulative probability function


  

t−µ
2
t+µ
F (t; µ, α) = Φ
+ exp
Φ − √
,
√
α µx
α2
α µx

(43)

where Φ is the cumulative probability distribution function of a standard normal
random variable, i.e., F is the inverse Gaussian or Wald distribution.
Forward recurrence time In practice, we usually do not start our observation at theoccurrence of an event. Then we are interested in the distribution of waiting times
until the next event starting from this starting time, i.e., the forward recurrencetime from an arbitrary time. This problem is equivalent to finding the distributionof the distance between an arbitrary point to the next event in the time line, andcan be solved in three steps.
1. Assume that the
R ∞length of each inter-event interval ∆ has a density f with finite
mean µ0 , i.e., 0 u f (u) du = µ0 .
2. The arbitrary time point falls in each inter-event interval ∆ with a probability
proportional to its length. That is to say, this arbitrary time is uniformly distributed on the time axis and has more chances to fall in larger intervals; or,Pr{Arbitrary point falls in ∆ | ∆ ∈ [z, z + dz]} ∝ z.
3. From 2, given the length ∆ = z of the inter-event interval that the arbitrary time
falls in, the probability that the forward recurrence time X is less than t is 1 ifz < t, or t/z if z ≥ t. In other words, the probability that the forward recurrencetime X is greater than t is 0 if z < t, or (z − t)/z if z ≥ t. In summary,
Z ∞h
Pr{X ≥ t} =Pr{Arbitrary point falls in ∆ | ∆ ∈ [z, z + dz]}ti
× Pr{∆ ∈ [z, z + dz]} × Pr{X > t | ∆ ∈ [z, z + dz)}Z ∞(z − t)z f (z)
∝
dz
z
t
Z ∞
=
(z − t)f (z) dzt= µc (t) − t[1 − F (t)](44)

Basic models of seismicity: Temporal models

15

where
Z ∞

c

z f (z) dz.

µ (t) =tSimilarly, it can be shown that
Pr{X < t} ∝ t[1 − F (t)] + µ(t)

(45)

where
Z t
µ(t) =z f (z) dz.
0

We find the normalizing factor in (44) by using
1 = Pr{X < ∞} = Pr{X ≥ 0},(46)

t[1 − F (t)] µ(t)
+
µ0
µ0

(47)

µc (t) t [1 − F (t)]
−
,
µ0
µ0

(48)

which gives
Pr{X < t} =and
Pr{X ≥ t} =where µ0 = E ∆ =R∞0z f (z) dz.

It is easy to see that the phrase “the forward recurrence time” can be replaced
by “the backward recurrence time” in the above discussion. That is to say, thebackward recurrence time Y has the same distribution as the forward distributionfunction. The density function wf of the forward recurrence times is obtained bytaking the derivative of (47), i.e.,
1 − F (t)
,
µ0

(49)

wf (t)
1 − F (t)= c.
Pr{X > t}µ (t) − t[1 − F (t)](50)

wf (t) =with corresponding hazard function
hf (t) =16

www.corssa.org

Table 1 Historical Earthquakes at Nankai Trough (DATA FILE: example2.data)
Occurrence times
t1
t2
t3
t4
t5
t6
t7
t8
t9
t10

Years
684
887
1099
1233
1361
1498
1605
1707
1854
1946

Table 2 Results from fitting renewal models to the Nankai Trough data
Model
PoissonGamma
log-normal
WeibullBPT

MLE
λ = 7.092 × 10−3
k = 7.768, θ = 16.865µ = 4.813, σ = 0.365
k = 2.886, µ = 144.849µ = 123.881, α = 0.419

log L-59.488-54.156-54.580-55.176-55.978

AIC
120.975
112.312113.160
114.351113.954

Likelihood function Suppose that the realization of a renewal process is {t1 , t2 , · · · , tn }in [0, T ], and we do not know what happens before time 0. The likelihood function
can be written as
L(0, T ; ·) = wf (t1 ) ×nYf (ti − ti−1 ) × F (T − tn ).

(51)

i=2

Another form of the likelihood is through the conditional intensity function,
Z T
X
λ(u) du.
log L(0, T ; ·) =log λ(ti ) −i: ti ∈[0,T ]

(52)

0

It can be shown that both forms are equivalent.
Example 2 Table 1 [Data file: example2.data] lists 10 great earthquakes of M 8+
that have occurred in the Nankai Trough south of Japan since 600 A.D., taken from
Table 10.5 in Utsu (1999). Ogata (2002) analyzed this data set by using recurrence
models and Bayesian methods to evaluate the uncertainties in the estimated occurrence times. Here, as a simple illustration, we fit the Poisson and the four renewal
models discussed in this section. The estimated parameters and the corresponding likelihood and AIC scores are listed in Table 2, indicating that, in this simple
example, the Gamma model fits the data the best.

Basic models of seismicity: Temporal models

17

Among the models mentioned in this section, the BPT model has become themost popular one and has been used in obtaining some interesting results (see, e.g.,
Ogata 2002 and Nomura et al. 2011). However, there is no clear evidence showingthat one model is superior to the others. This might be caused by the relatively
large standard errors of the parameter estimates that are due to very few data
points in the historical catalogs. Another reason may due to the nature of renewal
processes (e.g., MacFadden and Weissblum 1965; Daley and Vere-Jones 2003, Pages
77–79): The superposition of independent nontrivial renewal processes does not
result in a renewal process unless all of the superposed processes are simply Poisson.
This property requires that the renewal process is applied to a fault system thatis relatively disjoint in behavior from the surrounding environment, or is loaded bya constant tectonic stressing rate. However, the identification of such fault systems
seems technically difficult in practice.
3.4 Stress release models
Based on the elastic rebound theory (Reid 1910), the stress release model was introduced in a series of papers by Vere-Jones and others (see, e.g., Zheng and Vere-Jones
1991, 1994). It was extended to the linked stress release model by Shi et al. (1998),
Liu et al. (1998), Lu et al. (1999), Lu and Vere-Jones (2000), and Bebbington andHarte (2001, 2003). This model assumes that the stress level in a certain region
is gradually built up, linearly with time by tectonic movements, and drops down
suddenly coinciding with earthquakes, i.e., the stress level at t, X(t), can be writtenasX(t) = X(0) + ρt − S(t),(53)
where X(0) is the initial stress level, ρ is the constant loading rate from the tectonic
movements, and S(t) is the accumulated stress release from earthquakes in the
period [0, t), namely
X
S(t) =si .
(54)
i: ti <t

It is assumed that the stress release related to individual earthquakes scales withthe magnitude mi according tosi = 100.75mi ,(55)

which is hinted from the Benioff strain. The conditional intensity is a monotonicallyincreasing function of the stress level, usually chosen as
λ(t) = Ψ (X) = eγX = exp[γ(X(0) + ρt − S(t))] = exp[a + bt − c S(t)],

(56)

18

www.corssa.org

where γ is a single parameter representing the strength and heterogeneity of theregional crust, and a = γ X(0), b = ργ, and c = γ are the final model parameters. The assumption of an exponential distribution is in agreement with laboratory
experiments of stress corrosion where the mean waiting time until fracture can beapproximated by an exponential function of negative applied stress (Scholz 2002).
The amount of stress si released by an earthquake has a distribution with adensity function of ξ(s) and cumulative probability function Ξ(s), which is, forsimplification, assumed to be independent from the stress level. That is to say, thefull model has a conditional intensity function in the form of
λ(t, s) = Ψ [X(t)] ξ(s).

(57)

Given the observations {(ti , si ) : i = 1, 2, · · · , n} in a time period [T1 , T2 ], the loglikelihood function for the full model is
Z T2 Z ∞n
X
log L(T1 , T2 ; ·) =λ(u, v) dv dulog λi (ti , si ) −
=

i=1
n
X
i=1

T1

0

Z T2
log Ψ (X(ti )) −

Ψ (X(u)) du +T1nXlog ξ(si ).

(58)

i=1

In the above equation, the first two terms are from the usual likelihood for pure
temporal processes, and the third term represents the joint distribution of stress
reductions caused by the earthquakes.
Usually, we consider the following three types of distributions for S: the Pareto,
the truncated Pareto and the tapered Pareto distributions:
1. The Pareto distribution of S corresponds to the case where the magnitude distribution follows the Gutenberg-Richter magnitude-frequency relation, i.e., theexponential distribution. It has a probability density function
 −k
k−1 s
ξ(s) =; s ≥ S0 , k > 1,
(59)
S0
S0
where k is linked with the Gutenberg-Richter b-value by k = 43 b + 1.
2. The truncated Pareto distribution of s corresponds to the case where the magnitude distribution follows the truncated Gutenberg-Richter magnitude-frequency
relation. It has a probability density function
 −k
s
k−1

ξ(s) =, S 0 ≤ s ≤ Su ,(60)
 1−k  S
0Su
S0 1 − S0
where Su and S0 are the upper and lower thresholds, respectively.

Basic models of seismicity: Temporal models

19

3. The tapered Pareto (or Kagan) distribution (see, Vere-Jones et al. 2001; Kaganand Schoenberg 2001) has a cumulative distribution function s −k+1
S0 − sΞ(s) = 1 −, S0 ≤ s < ∞,
(61)
exp
S0
Sc
and a densityξ(s) =

k−1
1
+
s
Sc



s
S0

−k+1


exp

S0 − sSc

, S0 ≤ s < ∞,(62)

where Sc is a parameter governing the strength of the exponential taper affecting
frequency of large events.
All three distributions are illustrated in Fig. 1.
The magnitude-frequency distribution of earthquake is well described by the
Gutenberg-Richter magnitude-frequency relation (Gutenberg and Richter 1944),
which takes the form of an exponential distribution, i.e.,Pr{magnitude > m} = 10−b(m−mc ) , m ≥ mc ,(63)

where b is the so-called Gutenberg-Richter b-value. To reduce the possibility of
extremely large earthquakes in simulation, the following truncated form is also
often adopted:
 −b(m−m )c10= e−β(m−mc ) , if mc ≤ m ≤ Mmax ,Pr{magnitude > m} =0otherwise,
where Mmax is the upper bound of magnitudes.

Example 3 The following example is taken from Zheng and Vere-Jones (1994) where
the stress release model was fit to the historical catalog from North China during
the period from 1480 to 1992. Table 3 (Data file example3.data) lists this historicalcatalog from North China, Figure 2 shows the division of subregions, and Table
4 gives the corresponding results. Figure 3 shows the temporal variations of theconditional intensity functions in each subregion.
4 Temporal clustering models
On short time scales, seismicity has been shown to be clustered. The terms foreshock,main shock, aftershock, and earthquake swarm are concepts related the clustering of20

ξ(s) S0

10−5

1www.corssa.org

10

−15

10

−10

ParetoTruncated Pareto
Tapered Pareto

1

10

103

100

104

105

s S0

Fig. 1 Probability density functions ξ(s) of the Pareto, the truncated Pareto, and the tapered Pareto distribution.
Here we take k = 7/3(b = 1), Sc = 104 S0 , and Su = 103.5 S0 .

Fig. 2 Subregions and earthquake locations for the North China dataset [from Zheng and Vere-Jones (1994)].

the earthquake events. They are often used loosely in the many studies on seismicity.
However, a precise definition of each is generally not agreed, and indeed, it is often
difficult to correctly label observed events prospectively as they occur.
Short-term statistical seismicity models developed in the last century, including
the Omori-Utsu formula, the multiple Omori-Utsu formula, and the epidemic-typeBasic models of seismicity: Temporal models

21

Table 3 DATA FILE: example3.data Historical large earthquake from north China, 1480-1996, reprinted from
Zheng and Vere-Jones (1994).
Date1484.01.29
1501.01.19
1536.10.22
1556.01.23
1568.04.25
1573.01.10
1597.10.06
1614.10.23
1618.11.16
1622.10.25
1624.04.17
1626.06.28
1634.01.1654.07.21
1665.04.16
1679.09.02
1695.05.18
1709.10.14
1720.07.12
1739.01.03
1820.08.03
1830.06.12
1852.05.26
1879.07.01
1885.01.14
1888.11.02
1922.09.29
1945.09.23
1967.03.27
1975.02.04
1976.07.28
1979.08.25

Lat.
40.40
34.80
39.60
34.50
34.40
34.40
38.50
37.20
39.80
36.50
39.80
39.40
34.10
34.30
39.90
40.00
36.00
37.40
40.40
38.80
34.10
36.40
37.50
33.20
34.50
37.10
39.20
39.50
38.50
40.70
39.40
41.20

Long.
116.10110.10116.80109.70109.00104.10120.00112.50114.50106.30118.80114.20105.30105.50116.60117.00111.50105.30115.50106.50113.90114.20105.20104.70105.70104.20120.50119.00116.50122.80118.00108.10

Mag.
6.70
7.00
6.00
8.00
6.70
6.70
7.00
6.50
6.50
7.00
6.20
7.00
6.00
8.00
6.50
8.00
8.00
7.50
6.70
8.00
6.00
7.50
6.00
8.00
6.00
6.20
6.50
6.20
6.30
7.30
7.80
6.00

Reg.
W1
E2
W1
E2
E2
E1
W2
E2
W1
E1
W1
W1
E1
E1
W1
W1
E2
E1
W1
E1
W1
W1
E1
E1
E1
E1
W2
W1
W1
W2
W1
E1

Date
1487.08.10
1502.10.17
1548.09.13
1561.07.25
1568.05.15
1587.04.10
1604.10.25
1618.05.20
1622.03.18
1624.02.10
1624.07.04
1627.02.15
1642.06.30
1658.02.03
1668.07.25
1683.11.22
1704.09.28
1718.06.29
1730.09.30
1815.10.23
1829.11.19
1831.09.28
1861.07.19
1882.12.02
1888.06.13
1920.12.06
1937.08.01
1966.03.22
1969.07.18
1976.04.06
1976.09.23
1989.10.18

Lat.
34.30
35.70
38.00
37.50
39.00
35.20
34.20
37.00
35.50
32.40
35.40
37.50
35.10
39.40
35.30
38.70
34.90
35.00
40.00
34.80
36.60
32.80
39.70
38.10
38.50
36.70
35.40
37.50
38.20
40.20
39.90
40.00

Long.
108.90115.30121.00106.20119.00113.80105.00111.90116.00119.50105.90105.50111.10115.70118.60112.70106.80105.20116.20111.20118.50116.80121.70115.50119.00104.90115.10115.10119.40111.10106.40113.70

Mag.
6.20
6.50
7.00
7.20
6.00
6.00
6.00
6.50
6.00
6.00
6.00
6.00
6.00
6.00
8.60
7.00
6.00
7.50
6.50
6.70
6.00
6.20
6.00
6.00
7.50
8.50
7.00
7.20
7.40
6.20
6.20
6.00

Reg.
E2
W1
W2
E1
W2
W1
E1
E2
W1
W2
E1
E1
E2
W1
W2
E2
E1
E1
W1
E2
W2
W2
W2
W1
W2
E1
W1
W1
W2
E2
E1
E2

Table 4 Results from fitting the stress release model to the historical earthquakes from each of the four subregionsin North China. N represents the number of events in each subregion.

E1
E2
W1
W2

a
-5.397-3.900-2.764-5.152b0.02090.01760.00350.0238

c
0.571
0.880
1.731
0.410

log L (SRM)-52.04-83.34-50.36-91.97log L (Poisson)-56.76-87.58-52.99-77.84

∆AIC
-5.44
-4.48
-1.24
-7.75

N
11
22
11
18

aftershock sequence (ETAS) model, all emphasize the evolutionary nature of a developing earthquake cluster. In this section, we give an introduction to these models.
The models described in this section do not have a spatial component. Hence, there
is an implicit assumption that the spatial area is sufficiently small so that any given

22www.corssa.org

Fig. 3 Conditional intensity curves for the four subregions of the North China dataset [from Zheng and Vere-Jones(1994)].

event can conceivably interact with all following events, regardless of their spatial
locations.
4.1 The Omori-Utsu formula
The Omori-Utsu formula describes the decay of the aftershock frequency with timeafter a mainshock as an inverse power law. When Omori (1894) was studying theaftershocks of the 1891 Ms 8.0 Nobi earthquake, he first tried to use an exponential
decay function to fit the data, but unsatisfactory results were obtained. Then hefound that the number of aftershocks occurring each day can be well described by
the equationn(t) = K(t + c)−1 ,(64)

Basic models of seismicity: Temporal models

23

where t is the time from the occurrence of the mainshock, K and c are constants.
Utsu (1957) postulated that the decay of the aftershock numbers could vary, andshowed thatn(t) = K(t + c)−p(65)yields better fitting results. Equation (65) is now called the modified Omori formulaor Omori-Utsu formula.
The Omori-Utsu formula has been used extensively to analyze, model and forecast
aftershock activity. Utsu et al. (1995) reviewed the values of p for more than 200
aftershock sequences and found that it ranges between 0.6 and 2.5 with a medianof 1.1. He found no clear relationship between estimates of p-values and mainshockmagnitudes.
In the context of statistical seismology, we usually make use of (65) in the formof a conditional intensity function, given by
λ(t) =

K.
(t + c)p(66)

When the magnitude distribution of the aftershocks is also considered, the conditional intensity for the full model is written as
λ(t, m) =K s(m),(t + c)p(67)

where s(m) is the magnitude probability density function and usually takes the formof the Gutenberg-Richter magnitude-frequency relation. This model has been used
by Reasenberg and Jones for forecasting aftershock activity (see, e.g., Reasenberg
and Jones 1989, 1994). Therefore, some researchers also call (67) the ReasenbergJones model.
Likelihood of the single and multiple Omori-Utsu formulas Suppose that the occurrence time of the mainshock is 0. As discussed in Section 2, given observations
N = {(ti , mi ) : i = 1, 2, · · · , n} in the the time interval [S, T ], T > S ≥ 0, thelog-likelihood function for the model of the Omori-Utsu formula can be written in
the formZ Z T
Xlog L(N ; S, T ) =log λ(ti , mi ) −λ(u, m) du dm
M

i: ti ∈[S,T ]=X
i: ti ∈[S,T ]S

Z Tlog λ(ti ) −λ(u) du +S

X

log s(mi ),(68)

i: ti ∈[S,T ]where the first two terms on the righthand side represent the contribution to thelikelihood from the occurrence times and the third term represents the contribution

24

www.corssa.org

Fig. 4 A plot of magnitudes against indices of the aftershocks of the 2008 Wenchuan earthquake in China. The
red line indicates the estimation of the completeness magnitude.

from the magnitudes. The parameters in the Omori-Utsu formula can be estimatedfrom maximizing the temporal part of the likelihood (first two terms in the righthand
side)
X

Z T
log λ(ti ) −= N log K − pλ(u) duS

i: ti ∈[S,T ]X
i: ti ∈[S,T ]log(ti + c) −K (S + c)1−p − (T + c)1−p .
p−1

(69)

Example 4 The aftershock sequence of the 2008 MS 8.0 Wenchuan earthquake, China:
Table 5 [Data file example4.data] lists the mainshock and the aftershocks (MS ≥ 4.0)of the Wenchuan earthquake (2009-5-12, Ms 8.0), occurring within 25 days after themainshock. Since aftershocks immediately after a large mainshock are usually missing due to detection and recording problems (see, e.g., Kagan and Jackson 1995),
the Omori-Utsu formula should only be fit to a period where the records in theconsidered magnitude range are complete. A simple and direct way to detect thisstarting point is by plotting each event magnitude against its catalog sequential
number of event. As shown in Figure 4, we arbitrarily take T0 = 0.3 days, corresponding to the occurrence time of the 40th event. The estimated parameters fromfitting the Omori-Utsu formula are K̂ = 45.228, ĉ = 0.129 days, and p̂ = 1.107; withlog L = 270.575.

Basic models of seismicity: Temporal models

25

Fig. 5 Comparison of the observed and modeled cumulative number of M ≥ 4-aftershocks of the Wenchuan
earthquake: observation (black) and Omori-Utsu formula (red).

Multiple Omori-Utsu formula It is often observed that not only mainshocks trigger
aftershocks, but also large aftershocks may trigger their own aftershocks. To modelsuch phenomena, Utsu (1970) used the following model−p

λ(t) = K/(t − t0 + c)NTXKi H(t − ti )+,(t − ti + ci )−pi
i=1

(70)

where t0 is the occurrence time of the mainshock, ti , i = 1, . . . NT , define the occurrence times of the triggering aftershocks and H is the Heaviside function. The
likelihood for the multiple Omori-Utsu formula is sightly more complicated than forthe simple Omori-Ustu formula, but can be written in a similar way.
Example 5 This example is on the 1965 Rat Island earthquake of MW 8.7 and its
aftershocks, and is taken from Ogata and Shimazaki (1984). Earthquakes within therange 170◦ − 180◦ E and 48◦ − 55◦ N and with magnitudes mb ≥ 4.7 are selected.
An ordinary log-log plot of the daily number of events against time is time shownin Figure 6. The results from fitting a simple Omori-Utsu formula areK̂ = 85.088,ĉ = 0.204,p̂ = 1.055,

AIC = −811.4,

and from the multiple formula, by assuming that the largest aftershock of MW 7.6
also has its own aftershocks, areK̂ = 82.284,

K̂1 = 6.117,ĉ = ĉ1 = 0.204,p̂ = p̂1 = 1.055,

AIC = −873.4.

26

www.corssa.org

Fig. 6 An ordinary log-log plot of the number of aftershocks of the 1965 MW 8.7 Rat Island event per day against
time (Ogata and Shimazaki 1984). The curve represents the intensity function fitted to observed data points marked
by the squares.

Here the AIC criteria selects the model where c = c1 and p = p1 as the best model
among the class of models with NT = 1.
One difficulty in applying the multiple Omori-Utsu formula is to determine which
earthquakes are triggering events. The largest aftershocks often have secondary aftershocks, but not always. We can use the techniques of residual analysis in Section
5, as a diagnostic tool, to find out which events have secondary aftershocks.

4.2 The Epidemic-Type Aftershock Sequence (ETAS) model
In most cases, aftershock activity consists of secondary aftershock clustering as described by the multiple Omori-Utsu formula (70). However, even though possible
triggering events could potentially be determined through a residual analysis, separating triggering earthquakes from the others is difficult. Ogata (1988) generalized
this model by not needing to make a distinction between triggering events and theother events. He proposed that each event, irrespective of whether it is a small ora big event, can in principle trigger its own offspring. More precisely, the temporal

Basic models of seismicity: Temporal models

27

conditional intensity of this model isλ(t) = µ +Xκ(mi ) g(t − ti ),(71)

i: ti <tor, the full form of the condition intensity for this marked point process is"#X
λ(t, m) = s(m) µ +κ(mi ) g(t − ti ) ,(72)

i: ti <twhere s(m) = β e−β(m−m0 ) , m ≥ m0 , is the p.d.f. form of the Gutenberg-Richter
relation, m0 being the magnitude threshold, κ(m) = A exp[α(m − m0 )] is the mean
number of events directly triggered by an event of magnitude m and g(u) =(p − 1)(1 + t/c)−p /c is the probability density function (p.d.f.) of the time difference between the parent event and its children, i.e., the p.d.f. form of the OmoriUtsu formula. Ogata (1988) named (71) the Epidemic-Type Aftershock Sequence(ETAS) model, based on an analogy with the spread of epidemics. The ETAS model
also belongs to the more general class of self-exciting Hawkes processes (Hawkes
1971a,1971b; Hawkes and Oakes 1974).
Example 6 [Data file: example6.data] The aftershock sequence of the Darfield earthquake, New Zealand (MW 7.1, 2010-9-4). The dataset is taken from the New Zealandcatalog compiled by GNS, in the region (43.2 − 44.0◦ S, 171.0 − 173.5◦ E) and in thetime interval from 1985-1-1 to 2011-7-13. The magnitude threshold is M 3.8. By using the ETAS fitting program in the SASeis software (Ogata 2006), the MLEs of themodel parameters are µ = 0.001414 (events/day), K = 10.32, c = 0.006006 (day),
α = 1.832, p = 1.0167. Figure 7 shows the conditional intensity for the fitted ETAS
model.
Criticality and branching ratio A condition for stability of the ETAS model is thatthe process must be subcritical. In data analysis this important condition justifies
whether estimated model parameters are reasonable or irrational. Also, when extending the ETAS model or specifying an explicit form of an ETAS-like branchingmodel, the stability conditions must be satisfied. Forecasting with model parametersthat specify a supercritical model definitely overestimate the earthquake risk in themedium- or long-term future.
The stability of the ETAS model or a more general branching process is closely
related to the concepts of criticality and branching ratio. These two concepts are
identical for the ETAS model, but not the same for general branching processes.
To illustrate the differences between them, we start our discussions on criticality
and the branching ratio with more general assumptions. Readers who would like toknow more details are recommended to read through the following box.

28

www.corssa.org

10
1
6
4

5

Magnitude

7

10

−3

10

−2

0.1

Events/day

10
1
0.1
6
5
4

Magnitude

7

10

−3

10

−2

Events/day

100 103(b)

100 103

(a)

6000

8000

10000 12000 14000

Days from 1970−1−1

14800

14900

15000

15100

Days from 1970−1−1

Fig. 7 Temporal variation of the conditional intensity from fitting the ETAS model to the Darfield earthquake
sequence for (a) the whole time period, and for (b) the time period of the aftershock sequence only. The magnitudes
against occurrence times are plotted in lower panels.

Basic models of seismicity: Temporal models

29

We start our discussions with more general branching models with the following assumptions:(1) The mean number of direct children produced by an event of magnitude m0
is a Poisson random variable with a mean of κ(m0 );(2) The magnitudes of the children from a parent of magnitude m0 are independent and identical distributed (i.i.d.) random variable samples (r.v.s) with adensity s(· | m0 );(3) The intensity function of background events is G0 (m) = C s0 (m), wheres0 (m) is the p.d.f. of the magnitude density for background events and C isa constant.
It is easy to see that the intensity function for the first generation (events ofdirect children from background events) isZ
G1 (m) =κ(m0 ) s(m | m0 ) G0 (m0 ) dm0 ,S

where S represents the range of magnitudes, and for the second generation isZ
κ(m0 ) s(m | m0 ) G1 (m0 ) dm0 .
G2 (m) =S

Similarly,
Z

0

0

0

0

Z

κ(m ) s(m | m ) Gn (m ) dm =Gn+1 (m) =S

K [n] (m, m0 ) G0 (m0 ) dm,

S

where K [n] (n ≥ 1) is defined by the following recursion:
K [1] (m, m0 ) = κ(m0 ) s(m | m0 ),Z[n]0
K (m, m ) = κ(m∗ ) s(m | m∗ ) K [n−1] (m∗ , m0 ) dm∗ ,for n = 2, 3, · · · .

Suppose that a(m0 ) and b(m) are the left and right eigenfunctions of K corresponding to the maximum eigenvalue %, i.e.,Z0% a(m ) =a(m) K(m; m0 ) dm,
(73)
S

and

Z
% b(m) =S

K(m; m0 ) b(m0 ) dm0 ,

(74)

30

www.corssa.org

respectively, satisfyingZa(m) b(m) dm = 1.

(75)Ω(m; m0 ) = a(m0 ) b(m),

(76)

S

Let
i.e., Ω is the projection operator of K corresponding to %, or
Z
Z
∗
∗
0
∗
Ω(m; m ) K(m ; m ) dm = K(m; m∗ ) Ω(m∗ ; m0 ) dm∗
SS= % Ω(m; m0 ).

(77)

When both κ and s(m | m0 ) are stepwise continuous functions, i.e., linear integral
equations (73) and (74) can be viewed as the continuum limit of eigenvalue
equations of the form
X
Mi,j vj = %vi ,jthe maximum eigenvalue is separated from the others, and when n → ∞,
K [n]
→ Ω,
%n

(78)

and thusGn (m) → %nZΩ(m; m0 ) G0 (m0 ) dm0 .

(79)

S

From (79), we can see that if % < 1, then Gn → 0 when n → ∞, and that if% > 1, then Gn → ∞ when n → ∞. Here % is called the criticality parameterbecause if % < 1, the process is stable; otherwise the population might explodeand become infinitely large.
Equation (79) can be rewritten asZnGn (m) → %b(m) a(m0 ) G0 (m0 ) dm0SZn= % b(m) a(m0 ) G0 (m0 ) dm0
S= %n b(m) × const,(80)

implying that b(m) is asymptotically proportional to the intensity of the population when n → ∞.

Basic models of seismicity: Temporal models

31

The eigenfunction a(m0 ) can be interpreted as the asymptotic ability in producing offspring, directly and indirectly, from an ancestor {m0 } because
Z
∞
∞ Z
X
X
[i]
0
i
K (m; m ) dm = lim%Ω(m; m0 ) dm
lim
n→∞

i=n

n→∞

S

= lim

n→∞

i=n
∞
X

S
i

0

Z

% a(m )

i=nnb(m) dmS%a(m0 ) × const.
n→∞ 1 − %

= lim

(81)

For the ETAS model, where the magnitude density is separable and the background rate is constant, the eigenvalue equations areZ00% a(m ) = κ(m ) a(m) s(m) dm, and(82)
SZ% b(m) = s(m) κ(m0 ) b(m0 ) dm0 ,(83)
Swhere M is the magnitude range. We can see thata(m0 ) = C1 κ(m0 ),(84)b(m) = C2 s(m).

(85)

and
0

Substituting a(m ) and b(m) back into (82) and (83), we getZ%=κ(m) s(m) dm.

(86)

S

The branching ratio is defined as the proportion of triggered events in all theevents. To obtain the branching ratio for general branching models, consider amodel with the following intensityXκ(mi ) g(t − ti ) s(m | mi ).
(87)
λ(t, m) = µ s0 (m) +ti : ti <tTaking expectations on both sides with respect to time, we get,"#X
λ̄ s1 (m) = E [λ(t, m)] = µ s0 (m) + E
κ(mi ) g(t − ti ) s(m | mi ) ,ti : ti <t

(88)

32

www.corssa.org

where λ̄ is the total average rate, and s1 (m) is the magnitude density for overall
events. The expectation of the summation in the right-hand side can be writtenas"#X
κ(mi ) g(t − ti ) s(m | mi )Eti : ti <tZ Z t

∗

∗

∗

∗



=E
κ(m ) g(t − u) s(m | m ) × λ̄ s1 (m ) dut dm
S −∞Z= λ̄ κ(m∗ ) s(m | m∗ ) s1 (m∗ ) dm∗ ,Sand henceZκ(m∗ ) s(m | m∗ ) s1 (m∗ ) dm∗ .

λ̄ s1 (m) = µ s0 (m) + λ̄

(89)

S

Integrating on both sides with respect to m givesZλ̄ = µ + λ̄ κ(m∗ ) s1 (m∗ ) dm∗ .

(90)

S

The branching ratio is obtained byµω =1− =λ̄Zκ(m∗ ) s1 (m∗ ) dm∗

(91)

S

which is also the average number of events that are triggered by an arbitrary
event.
For the ETAS model, the magnitude is completely separable from the whole
intensity, and, thus, ω is identical to the criticality parameter %. But this does
not hold for general cases.
Example 7 : Here we assume that h(m) is a probability density function andH(m) is the corresponding c.d.f. Sets(m | m0 ) =h(m),H(am0 )

0 < m < am0 ,0 < a ≤ 1.,

and set mc = 0 to abbreviate the notation. That is, the magnitudes of offspring
are limited to not be greater than am0 , where m0 is the magnitude of the ancestorof interest. By Equation (82),Z a m0κ(m0 )
0h(m)V (m) dm.
ρV (m ) =H(a m0 ) 0Basic models of seismicity: Temporal models

33

To obtain the eigenvalue, we apply the limit operation to both side of the above
equality. When Mi → 0, notice that H(0) = 0 and h(M ) = H 0 (M ), thenR a m0h(m)V (m) dmaκ(0)V (0)0ρV (0) = κ(0) lim
=
.
0
m0 →0
H(am )aThus the criticality parameter ρ is
ρ = κ(0).
In this example, if κ(m) is a monotonically increasing function of m, then, by(91), the branching ratio is
Z
Z
∗
∗
∗
ω=
κ(m ) s1 (m ) dm >κ(0) s1 (m∗ ) dm∗ = κ(0) = %.
M

M

For the ETAS model, substituting κ(m) = A eα(m−m0 ) and s(m) = β e−β(m−m0 )into (86), the criticality parameter is
Z ∞
Aβ
.
%=
s(m) κ(m) dm =β−αm0When % < 1, the ETAS model is stable and stationary, which requires β ≥ α andA ≤ 1 − α/β. When % ≥ 1, there is a finite probability that the number of events ina unit time interval becomes infinite as t increases to infinity. More details on the
behavior of the ETAS model were discussed by Helmstetter and Sornette (2002),
Zhuang and Ogata (2006) and Saichev and Sornette (2007).
Applications There have been many applications of the ETAS model. For example:(i) detection of anomalous seismicity patterns (e.g. Ogata 2005) (ii) characterizationof clustering characteristics by (regional) variations of the ETAS parameters (e.g.
Enescu et al. 2009); (iii) detection of fluid-related forcing signals (e.g. Hainzl andOgata 2005; Lombardi et al. 2010); (iv) regional probabilistic earthquake forecasting(e.g., Helmstetter et al. 2006), and so forth. Here we refer to Zhuang et al. (2011)
for further discussions.
5 Transformed time sequence and residual analysis
Residual analysis is a powerful tool for assessing the fit of a particular model to aset of occurrence times (e.g., Papangelou 1972; Ogata 1988; Daley and Vere-Jones34

www.corssa.org

2003). Assume that we have a realization of a point process with event times denoted by t1 , t2 , · · · , tn . We then calculate transformed event times, denoted byτ1 , τ2 , · · · , τn , in such a way that the transformed times have the same distributional properties as the a homogeneous Poisson process with unit rate parameter.
Suppose that the conditional intensity is λ(t). The transformed time sequence, fori = 1, 2, · · · , n, is calculated as
Z ti
τi =

λ(u) du.

(92)

0Then the sequence {τi : i = 1, 2, · · · , n} forms a Poisson process with unit rate.
To prove
R t the above statement, we need only to prove that, from any transformed
time τ = 0 λ(u) du, the waiting transformed time Y to next event is an exponential distribution with the unit rate. If the true model for the point process has a
conditional λ(t), by (5) and the equivalence between the conditional intensity andthe hazard function, the survival function of the waiting time X to next event fromt is
 Z t+x
St (x) = exp −λ(u) du .
(93)t

To continue our proof, we require the following theorem: If X is a continuous
random variable with a strictly increasing distribution function F , then F has aninverse F −1 defined on the open interval (0,1), and F (X) and 1 − F (U ) are both
uniformly distributed on (0,1). Moreover, if U is a uniform r.v. on (0, 1), then F −1 (U )is a r.v. with the distributionh R function Fi.
t+X
Thus, St (X) = exp − tλ(u) du is a uniform r.v. on (0, 1). Since Y =Rt
λ(u) du, we have St (X) = exp[−Y ], i.e., Y is exponentially distributed with
0
rate 1.
This method can be used to test the goodness-of-fitof the model. If theo fittedn
R ti
model λ̂(t) is close to the true model, then τ̂i = 0 λ̂(u) du : i = 1, 2, · · · , is close
to the standard Poisson model. The departure of the rate in the transformed time
sequence from the unit rate indicates either increased activity or quiescence, relative
to the seismic rate expected by the original model.
Example 8 (Continuation of Example 5.) A comparison between the cumulative
numbers of earthquakes (black curve) and the corresponding fit values (red curve,
Figure 5) shows that quiescence may have started on the 17th day after the mainshock. This is also confirmed in the plot of the transformed time domain (Figure 8).

Basic models of seismicity: Temporal models

35

Fig. 8 Comparison of the observed and modeled cumulative number of M ≥ 4-aftershocks of the Wenchuan
earthquake: observation (black) and Omori-Utsu formula (red) in the transformed time domain.

6 Related software
There are a couple of versions of available software for fitting the models discussed
in this paper. The first is IASPEI software SASeis implemented by Utsu and Ogata(1997) which is implemented using FORTRAN code, and later is revised as SASeis2006 by Ogata (2006) (URL: http://www.ism.ac.jp/∼ogata/Ssg/softwares.html).
An alternative is provided by Harte (2010, 2012) in the R package “PtProcess”(URL: http://cran.r-project.org/web/packages/PtProcess/), which uses an object
orientated approach. R is a statistical computing language developed by R Development Core Team (2012).
7 SummaryIn this article, we have illustrated a developmental approach to modeling earthquake data using statistical models. By applying a sequence of progressively moresophisticated and possibly competing models, we can refine our understanding of
the earthquake process. Starting from the model of complete randomness, the Poisson model, new models are constructed by adding deterministic factors of physicalhypotheses or empirical observations relations. For example, the hypothesis of characteristic recurrence of earthquakes leads to the renewal models, the elastic-reboundtheory leads to the stress release model, and empirical Omori-Utsu formula leads
to the ETAS model. Many interesting studies, extensions and applications of the36

www.corssa.org

above-mentioned models exist. We cannot go through all of them in the scope ofour discussions. Instead, we provide the reader with a general overview of some of
the most important point-process models of seismicity and the relevant methods formodel fitting and inference.
Acknowledgements The authors thank the editor, Jeremy D. Zechar, and two
anonymous reviewers for their constructive comments.

References
Akaike, H. (1974), A new look at the statistical model identification, Automatic Control, IEEE Transactions on,
19 (6), 716 – 723, doi:10.1109/TAC.1974.1100705. 8
Bebbington, M., and D. Harte (2001), On the statistics of the linked stress release process, Journal of Applied
Probability, 38A, 176–187. 17
Bebbington, M., and D. Harte (2003), The linked stress release model for spatio-temporal seismicity: formulations,procedures and applications, Geophysical Journal International, 154, 925–946. 17
Cornell, C. A. (1968), Engineering seismic risk analysis, Bullitin of the Seismogological Society of America, 58 (5),1583–1606. 8
Daley, D. D., and D. Vere-Jones (2003), An Introduction to Theory of Point Processes – Volume 1: Elementrary
Theory and Methods (2nd Edition), Springer, New York, NY. 17, 33
Davis, P. M., D. D. Jackson, and Y. Y. Kagan (1989), The longer it has been since the last earthquake, the longer
the expected time till the next?, Bullitin of the Seismogological Society of America, 79 (5), 1439–1456. 11
Ellsworth, W. L., M. V. Matthews, R. M. Nadeau, S. P. Nishenko, P. A. Reasenberg, and R. W. Simpson (1999),
A physically-based earthquake recurrence model for estimation of long-term earthquake probabilities, U.S.
Geological Survey Open-File Report, pp. 99–522. 13
Enescu, B., S. Hainzl, and Y. Ben-Zion (2009), Correlations of seismicity patterns in southern california with surface
heat flow data, Bullitin of the Seismogological Society of America, 99, 3114 – 3123, doi:10.1785/0120080038. 33
Field, E. H. (2007), A summary of previous Working Groups on California Earthquake Probabilities, Bullitin of the
Seismogological Society of America, 97 (4), 1033–1053, doi:10.1785/0120060048. 11
Gardner, J. K., and L. Knopoff (1974), Is the sequence of earthquakes in southern California, with aftershocks
removed, Poissonian?, Bullitin of the Seismogological Society of America, 64 (5), 1363–1367. 8
Gutenberg, B., and C. F. Richter (1944), Frequency of earthquakes in California, Bull. Seis. Soc. Am., 34, 184–188.
19
Hainzl, S., and Y. Ogata (2005), Detecting fluid signals in seismicity data through statistical earthquake modeling,
Journal of Geophysical Research, 110 (B05), B05S07, doi:10.1029/2004JB003247. 33
Harte, D. (2010), PtProcess: An R package for modelling marked point processes indexed by time, Journal of
Statistical Software, 35 (8), 1–32. 35
Harte, D. (2012), PtProcess: Time Dependent Point Process Modelling, R package version 3.3-1. 35
Hawkes, A. G. (1971a), Spectra of some self-exciting and mutually exciting point processes, Biometrika, 58 (1),83–90, doi:10.1093/biomet/58.1.83. 27
Hawkes, A. G. (1971b), Point spectra of some mutually exciting point processes, J. Royal Stat. Soc. Series B(Meth.), 33 (3), 438–443. 27
Hawkes, A. G., and D. Oakes (1974), A cluster process representation of a self-exciting process, J. of Appl. Prob.,
11 (3), 493–503. 27
Helmstetter, A., and D. Sornette (2002), Subcritical and supercritical regimes in epidemic models of earthquake
aftershocks, Journal of Geophysical Research, 107 (B10), 2237. 33
Helmstetter, A., Y. Y. Kagan, and D. D. Jackson (2006), Comparison of short-term and time-independent earthquake
forecast models for Southern California, Bullitin of the Seismogological Society of America, 96 (1), 90–106, doi:
10.1785/0120050067. 33
Kagan, Y. Y., and D. D. Jackson (1995), New seismic gap hypothesis: Five years after, J. Geophys. Res., 100 (B3),3943–3959. 11, 24

Basic models of seismicity: Temporal models

37

Kagan, Y. Y., and L. Knopoff (1987), Random stress and earthquake statistics - time-dependence, GeophysicalJournal of the Royal Astronomical Socieity, 88 (3), 723–731. 13
Kagan, Y. Y., and F. Schoenberg (2001), Estimation of the upper cutoff parameter for the tapered pareto distribution, Journal of Applied Probability, 38A, 158–175. 19
Liu, J., D. Vere-Jones, L. Ma, Y. Shi, and J. Zhuang (1998), The principal of coupled stress release model and its
application, Acta Seismologica Sinica, 11, 273–281. 17
Lombardi, A. M., M. Cocco, and W. Marzocchi (2010), On the increase of background seismicity rate during the
1997–1998 Umbria-Marche, central Italy, sequence: apparent variation or fluid-driven triggering?, Bullitin of the
Seismogological Society of America, 100 (3), 1138–1152, doi:10.1785/0120090077. 33
Lu, C., and D. Vere-Jones (2000), Application of linked stress release model to historical earthquake data:comparison between two kinds of tectonic seismicity, Pure and Applied Geophysics, 157, 2351–2364, doi:
10.1007/PL00001087. 17
Lu, C., D. Harte, and M. Bebbington (1999), A linked stress release model for historical japanese earthquakes:
Coupling among major seismic regions,, Earth Planets Space, 51, 907–918. 17
Ma, L., and D. Vere-Jones (1997), Application of m8 and lin-lin algorithms to new zealand earthquake data, New
Zealand Journal of Geology and Geophysics, 40, 77–89. 11
MacFadden, J. A., and W. Weissblum (1965), Higher-oder properties of a stationary point process, Journal of theRoyal Statistical Society, Ser. B (Methodological), 25, 413–431. 17
Matthews, M. V., W. L. Ellsworth, and P. A. Reasenberg (2002), A Brownian model for recurrent earthquakes,
Bulletin of the Seismological Society of America, 92 (6), 2233–2250, doi:10.1785/0120010267. 13
Nishenko, S. P. (1991), Circum-pacific seismic potential 1989-1999, Pure Appl. Geophys., 135, 169–259. 11
Nishenko, S. P., and R. Buland (1987), A generic recurrence interval distribution for earthquake forecasting, Bulletinof the Seismological Society of America, 77 (4), 1382–1399. 11
Nomura, S., Y. Ogata, F. Komaki, and S. Toda (2011), Bayesian forecasting of recurrent earthquakes and predictive
performance for a small sample size, Journal of Geophysical Research, 116, B04,315, doi:10.1029/2010JB007917.
17
Ogata, Y. (1988), Statistical models for earthquake occurrences and residual analysis for point processes, Journalof the American Statistical Association, 83, 9 – 27. 26, 27, 33
Ogata, Y. (2002), Slip-size-dependent renewal processes and bayesian inferences for uncertainties, Journal of Geophysical Research, 107 (B11), 2268, doi:10.1029/2001JB000668. 16, 17
Ogata, Y. (2005), Detection of anomalous seismicity as a stress change sensor, Journal of Geophysical Research,
110, B05S06, doi:10.1029/2004JB003245. 33
Ogata, Y. (2006), Statistical Analysis of Seismicity - Updated Version (SASeis2006), Computer Science Monograph,
33, 1-28. The Institute of Statistical Mathematics, Tokyo, Japan. 27, 35
Ogata, Y., and K. Katsura (2004), Point-process models with linearly parameterized intensity for application toearthquake data, in Essays in Time Series and Allied Processes (Papers in honour of E.J. Hannan), vol. 23A,
edited by J. Gani and M. B. Priestley, pp. 291–310, Journal of Applied Probability, New York. 11
Ogata, Y., and K. Shimazaki (1984), Transition from aftershock to normal activity, Bulletin of the SeismologicalSociety of America, 74 (5), 1757–1765. 25, 26
Omori, F. (1894), On the aftershocks of earthquakes, Journal of the College of Science, Imperial University ofTokyo, 7, 111–200. 22
Papangelou, F. (1972), Integrability of expected increments of point processes and a related random change of scale,Trans. Amer. Math. Soc., 164, 438–506. 33
R Development Core Team (2012), R: A Language and Environment for Statistical Computing, R Foundation for
Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0. 35
Reasenberg, P. A., and L. M. Jones (1989), Earthquake hazard after a mainshock in California, Science, 243, 1173
– 1176. 23
Reasenberg, P. A., and L. M. Jones (1994), Earthquake aftershocks: Update, Science, 265, 1251 – 1252. 23
Reid, H. (1910), The Mechanics of the Earthquake, The California Earthquake of April 18, 1906, Report of theState Investigation Commission, Vol. 2, Carnegie Institution of Washington, Washington, DC, pp. 16–28. 11,
17
Ross, S. M. (2003), Introduction to Probability Models, Eighth Edition, 773 pp., Academic Press. 8
Saichev, A., and D. Sornette (2007), Theory of earthquake recurrence times, J. Geophys. Res., 112 (B04313), doi:
10.1029/2006JB004536. 33
Scholz, C. H. (2002), The mechanics of earthquakes and faulting (2nd ed.), Cambridge University Press. 18

38

www.corssa.org

Shi, Y., J. Liu, D. Vere-Jones, J. Zhuang, , and L. Ma (1998), Application of mechanical and statistical models
to study of seismicity of synthetic earthquakes and the prediction of natural ones, Acta Seismol. Sinica, 11,421–430. 17
Sykes, L. R., and W. Menke (2006), Repeat times of large earthquakes: Implications for earthquake mechanics andlong-term prediction, Bull. Seismol. Soc. Am., 96 (5), 1569–1596, doi:10.1785/0120050083. 11
Utsu, T. (1957), Magnitude of earthquakes and occurrence of their aftershocks, Zisin (J. Seismol. Soc. Jap.), 10,
35–45 (in Japanese). 23
Utsu, T. (1970), Aftershock and earthquake statistics (ii) – further investigation of aftershocks and other earthquake
sequences based on a new classification of earthquake sequences, Journal of the Faculty the Science, Hokkaido
Univesity, Ser. VII (Geophysics), 3, 197 – 266. 25
Utsu, T. (1999), Seismicity Studies: A Comprehensive Review, 876 pp., Univ. of Tokyo Press. 16
Utsu, T., and Y. Ogata (1997), Statistical analysis of seismicity., in Algorithms for Earthquake Statistics andPrediction. International Association of Seismology and Physics of the Earth’s Interior (IASPEI) Library
Volume 6., edited by J. Healy, V. Keilis-Borok, and W. Lee, pp. 13–94, IASPEI, Menlo Park CA. 35
Utsu, T., Y. Ogata, and R. S. Matsu’ura (1995), The centenary of the Omori formula for a decay law of aftershock
activity, Journal of Physics of the Earth, 43, 1–33. 23
Vere-Jones, D., R. Robinson, and W. Yang (2001), Remarks on the accelerated moment release model: problems ofmodel formulation, simulation and estimation, Geophysical Journal International, 144 (3), 517–531. 19
Weibull, W. (1951), A statistical distribution function of wide applicability, Journal of Applied MechanicsTransactions of the ASME, 18 (3), 293–297. 13
Zheng, X., and D. Vere-Jones (1991), Application of stress release models to historical earthquakes from North
China, Pure and Applied Geophysics, 135 (4), 559–576, doi:10.1007/BF01772406. 17
Zheng, X., and D. Vere-Jones (1994), Further applications of the stochastic stress release model to historical earthquake data, Tectonophysics, 229, 101–121. 11, 17, 19, 20, 22
Zhuang, J., and Y. Ogata (2006), Properties of the probability distribution associated with the largest event in an
earthquake cluster and their implications to foreshocks, Physical Review, E, 73, 046,134, doi:10.1103/PhysRevE.
73.046134. 33
Zhuang, J., M. Werner, S. Hainzl, D. Harte, and S. Zhou (2011), Basic models of seismicity: spatiotemporal models,Community Online Resource for Statistical Seismicity Analysis, doi:10.5078/corssa-07487583. 33

Appendix: List of aftershocks of the Wenchuan earthquake
Table 5: DATA FILE: example4.data
#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

year
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008
2008

month
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5

day
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12

hour
14
14
15
15
16
16
16
16
16
16
16
17
17
17
17
17
17

minute
28
43
34
40
10
21
26
35
36
47
50
1
7
23
30
42
54

lat.
31
31
31
31
31.2
31.3
31.5
31.4
31
32.2
32.6
32.2
31.3
32.3
32.4
31.4
31

long.
103.4
103.5
103.5
103.6
103.4
104.1
103.8
103.5
103.2
105.4
105.2
104.7
103.8
104.8
104.9
104
103.2

mag.
8
6
5
4.7
4.8
5.2
4.3
4.6
4.3
4.8
4.5
4.1
5
5.1
4.1
5.2
4.3

days since mainshock00.01040.04580.05
0.07080.07850.08190.08820.08890.09650.09860.10620.11040.12150.12640.13470.1431
Continued on next page ...

Basic models of seismicity: Temporal models(Continuation from previous page ...)
#
year
month
day
hour
minute
18
2008
5
12
18
2
19
2008
5
12
18
23
20
2008
5
12
19
10
21
2008
5
12
19
33
22
2008
5
12
19
41
23
2008
5
12
19
45
24
2008
5
12
19
52
25
2008
5
12
20
4
26
2008
5
12
20
6
27
2008
5
12
20
11
28
2008
5
12
20
15
29
2008
5
12
20
23
30
2008
5
12
20
29
31
2008
5
12
20
33
32
2008
5
12
20
54
33
2008
5
12
21
2
34
2008
5
12
21
7
35
2008
5
12
21
32
36
2008
5
12
21
36
37
2008
5
12
21
40
38
2008
5
12
21
55
39
2008
5
12
22
6
40
2008
5
12
22
9
41
2008
5
12
22
15
42
2008
5
12
22
26
43
2008
5
12
22
37
44
2008
5
12
22
46
45
2008
5
12
22
55
46
2008
5
12
23
5
47
2008
5
12
23
5
48
2008
5
12
23
16
49
2008
5
12
23
28
50
2008
5
13
0
28
51
2008
5
13
0
34
52
2008
5
13
1
1
53
2008
5
13
1
29
54
2008
5
13
1
54
55
2008
5
13
2
46
56
2008
5
13
2
55
57
2008
5
13
3
53
58
2008
5
13
4
8
59
2008
5
13
4
45
60
2008
5
13
4
51
61
2008
5
13
5
8
62
2008
5
13
5
51
63
2008
5
13
6
19
64
2008
5
13
6
24
65
2008
5
13
6
47
66
2008
5
13
7
38
67
2008
5
13
7
46
68
2008
5
13
7
54
69
2008
5
13
8
22
70
2008
5
13
8
54

lat.
32.2
31
31.4
32.6
32.4
32.4
32.6
32.6
32.2
31.4
32
32.7
31.4
31.4
31.3
31.1
31
31.2
32.9
31
32
32.5
31.9
32.2
31.3
32.2
32.7
32.4
31.3
31.3
30.9
31
31.2
32.5
30.9
31.3
31.3
32.4
31.9
31.3
31.4
31.7
32.4
31.3
32.5
31.9
32.2
31.3
31.9
31.2
31.3
31.3
32.6

long.
105.1
103.3
103.6
105.4
105.1
105
105.4
105.2
105.5
103.8
104.4
105.3
103.9
104.1
103.4
103.5
103.4
103.9
105.5
103.5
104.3
105.1
104.7
104.9
103.9
104.5
105.5
105
103.6
103.5
103.2
103.5
103.8
105
103.4
103.4
103.4
105
105.1
103.6
104
104.5
105.2
103.2
105.3
104.2
105
103.4
104.5
103.4
103.6
104
105.2

mag.
4.8
4.9
6
4.5
4.6
4
4.7
4.2
4.1
4.5
4.9
4.5
4.1
4.2
4.3
4.6
4.3
4.3
4
5.1
4.2
4
4.5
4.4
4
4.3
5.1
4.2
5
5.2
4.1
5
4.4
4.2
4
4.6
5
4.4
4.4
4.3
5.7
5.2
4.7
4.5
4.6
4.1
4
4.5
4
5.3
5.1
4.2
4.3

days since mainshock
0.1486
0.1632
0.19580.21180.21740.22010.225
0.23330.23470.23820.241
0.24650.25070.25350.26810.2736
0.27710.29440.29720.3
0.31040.31810.32010.32430.33190.3396
0.34580.35210.359
0.359
0.36670.375
0.41670.42080.43960.459
0.47640.51250.51880.559
0.56940.59510.59930.61110.641
0.6604
0.66390.67990.71530.72080.72640.74580.7681
Continued on next page ...

39

40

www.corssa.org(Continuation from previous page ...)
#
year
month
day
hour
minute
71
2008
5
13
9
7
72
2008
5
13
10
15
73
2008
5
13
10
33
74
2008
5
13
10
59
75
2008
5
13
11
0
76
2008
5
13
11
48
77
2008
5
13
12
45
78
2008
5
13
12
50
79
2008
5
13
13
25
80
2008
5
13
13
36
81
2008
5
13
13
37
82
2008
5
13
14
38
83
2008
5
13
15
7
84
2008
5
13
15
19
85
2008
5
13
15
51
86
2008
5
13
15
53
87
2008
5
13
16
11
88
2008
5
13
16
20
89
2008
5
13
17
41
90
2008
5
13
18
16
91
2008
5
13
18
36
92
2008
5
13
20
51
93
2008
5
13
21
13
94
2008
5
13
21
31
95
2008
5
13
23
10
96
2008
5
13
23
54
97
2008
5
14
0
23
98
2008
5
14
3
30
99
2008
5
14
3
51
100 2008514
6
3
101 200851488102 200851499103 2008514956
104 20085
14
10
54
105 20085141354
106 20085
141433
107 20085
14
1726
108 20085
14
1751
109 20085
14
1757
110 20085
14
18
0
111 20085
14
18
11
112 20085
14
18
18
113 20085
14
18
30
114 20085
14
21
29
115 20085151
17116 20085
151
33
117 20085
153
59
118 2008515
5
1
119 2008515
6
10120 200851589121 2008515
850
122 20085151227
123 2008
5
15
13
27

lat.
31.4
31.6
31.3
31
31.2
31.2
31
31.3
32.6
32.4
31
31.4
30.9
32.3
32.5
32.3
32.5
31.4
32.1
31.8
31.3
32.3
32.5
32.4
32.6
32.1
31.7
31.1
31
31.3
31.1
31.4
31.1
31.3
31.9
31.4
31.4
32.4
32.3
32.2
32.2
32.4
32.4
32.3
31.5
31.4
31.1
31.6
31.2
31.8
31.3
31.3
32

long.
103.7
103.9
103.6
103.3
103.5
103.7
103.3
103.4
105.2
105.2
103.5
103.8
103.4
105
105.3
105
105.2
103.9
104.4
104.3
103.6
105
105.5
105.1
105.5
104.9
104.3
103.3
103.3
103.6
103.4
103.8
103.5
103.4
104.2
103.9
104
104.2
104.8
104.7
104.5
105.1
105.2
105.1
103.8
103.5
103.5
104.2
103.6
104.4
103.4
103.7
104.3

mag.
4
4.5
4.2
4.3
4.7
4.5
4.1
4.1
4.2
4.4
4.6
4.3
6.1
4.8
4.4
4.8
4.1
4.8
4.1
4
4.1
4.6
4.4
4.5
4.1
4.5
4
4.1
4.2
4.2
4
4.8
4.1
5.6
4.9
4.1
5.1
4.8
4.2
4.7
4.2
4.8
4.9
4.5
4.8
4.6
4.2
5
4.6
4.3
4
4.1
4.6

days since mainshock
0.77710.82430.83680.85490.85560.88890.92850.93190.95620.96390.96461.00691.02711.0354
1.05761.059
1.0715
1.07781.134
1.1583
1.17221.266
1.2812
1.29371.3625
1.39311.4132
1.5431
1.55761.6493
1.73611.77851.8111
1.8514
1.97642.0035
2.12362.141
2.14512.14722.15492.15972.16812.29242.45072.46182.56322.60632.65422.73682.76532.916
2.9576
Continued on next page ...

Basic models of seismicity: Temporal models(Continuation from previous page ...)
#
year
month
day
hour
minute124 200851520
10
125 2008515
21
4
126 2008516
555
127 2008516
6
10128 2008516
634
129 200851611
10130 2008516
1134
131 2008516
13
25
132 2008516
1434
133 2008516
18
17
134 20085
16
18
20
135 20085
16
1851
136 20085
16
21
21
137 20085
170
14138 20085
171
22
139 20085
173
59
140 2008517
4
0
141 2008517
4
15142 2008517
4
16
143 2008517
4
29
144 2008517
633
145 20085177
23
146 2008517
8
28
147 2008517838
148 2008517
1538
149 20085
17
21
32
150 20085181
8151 2008518
426
152 2008518
845
153 2008518
94154 20085
18
1151
155 2008518
1725
156 2008518
2037
157 2008519128
158 2008519146
159 2008520
1
52
160 2008520
857
161 2008520
1142
162 2008520
1217163 2008520
1454
164 20085
21
038
165 20085
21
1640
166 20085
21
1733
167 20085
21
2159
168 20085
21
23
29
169 2008522436
170 200852215
18
171 20085221922
172 2008522
230
173 20085
231
37
174 20085238
5175 2008523
923
176 2008
5
23
11
12

lat.
31.4
32.6
32.3
31.4
31.9
32.5
31.4
31.4
32.4
31.3
32.5
31.4
31.9
31.2
31.2
31
32.6
32.2
31.3
31.4
32.2
31.3
31.6
32
32
32.2
32.1
31.2
31.8
31.1
31
31.2
31.3
32.1
32.5
32.3
31.7
32.6
30.8
31.8
30.9
31.4
32.3
31.4
32.4
32.2
31.2
32.6
31.9
31.3
31.2
31.2
31.2

long.
103.8
105.6
104.7
103.9
104.4
105.1
104.1
103.2
105.2
103.5
105.1
103.6
104.2
103.5
103.6
103.5
105.4
104.4
103.5
103.3
105.1
103.8
104
104
104.4
104.7
105
103.5
104
103.5
103.4
103.1
103.2
105
105.3
104.9
104
105.4
103.3
104.2
103.3
103.3
105.2
103.9
105.1
104.8
103.6
105.4
104.3
103.6
103.6
103.5
103.2

mag.
4.2
4.4
4.5
4.6
4.2
4
4.9
5.9
4.3
4.3
4
4.2
4
5.1
4.5
4.3
4.1
4.8
5
4.5
4.2
4.2
4.1
4
4.1
4.7
6
4.1
4.2
4.2
4.2
4
4.2
4.8
5.4
5
4.1
4
4.3
4
4
4.1
4.5
4.5
4.3
4.6
4.7
4
4.3
4.6
4.7
4.2
4.2

days since mainshock
3.23753.275
3.64373.65423.67083.86253.8792
3.95624.0042
4.159
4.1611
4.18264.28684.40694.45424.56324.56394.57434.575
4.584
4.67014.70494.75
4.75695.04865.29445.44445.58195.76185.775
5.891
6.12296.25626.90286.98477.475
7.7701
7.88477.909
8.01818.42369.09179.12859.31329.37579.588910.0347
10.2042
10.3556
10.4646
10.73410.7882
10.8639
Continued on next page ...

41

42

www.corssa.org(Continuation from previous page ...)
#
year
month
day
hour
minute177 20085240
10
178 20085
241
53
179 2008524
110
180 2008525
1227
181 20085
25
16
21
182 2008525
1734
183 2008526
839
184 2008527
16
3
185 20085
27
16
37
186 20085
27
21
59
187 2008528
046
188 2008528
1
35
189 2008529
1248
190 2008529
1510
191 2008531
1422
192 2008531
1534
193 200861
1123
194 200863
119
195 2008651
26
196 200865521
197 20086
51241
198 2008
6
5
14
2

lat.
32.2
32.5
31.1
32
32.6
33
30.8
32.7
32.8
32.5
32.2
32.7
32.6
31.4
32.4
32.6
31.6
32
32.3
31.2
32.3
32.7

long.
105
105.2
103.4
104.6
105.4
104.9
103.3
105.6
105.6
105.2
104.6
105.4
105.5
103.7
105
105.4
104
104.5
105.1
103.4
105
105.5

mag.
4
4.1
4
4.2
6.4
4.7
4.3
5.4
5.7
4.7
4.2
4.7
4.6
4.5
4
4
4.5
4.3
4
4.2
5
4.3

days since mainshock
11.4042
11.4757
11.8556
12.91613.0785
13.1292
13.7576
15.06615.089615.3132
15.4292
15.4632
16.9306
17.0292
18.9958
19.0458
19.8715
21.861823.4569
23.6201
23.9257
23.9819