Diagnostic analysis of space-time
branching processes for earthquakes
Jiancang Zhuang1 , Yosihiko Ogata1 and David Vere-Jones2
1 Institute of Statistical Mathematics

4-6-7 Minami Azabu, Minato-Ku, Tokyo 106-8659, Japan
2 Victoria University of Wellington

P.O.Box 600, Wellington, New ZealandFebruary 27, 2004
Abstract
This paper discusses several methods to test the goodness-of-fit of a spacetime continuous-type branching process with immigration, namely the ETAS
(epidemic type aftershock sequence) model, to the occurrences of earthquakes,
including: 1) using thinning method to decompose the whole process into abackground process and its complementary process, to test the stationarityof the background process, 2) using the residual analysis to test the overall
goodness-of-fit, and 3) using the stochastic reconstruction method to compare
the characteristics of clustering features between the earthquake data and themodel.

1

Introduction

It is natural to use the branching process to describe the occurrence pattern ofearthquakes, which is clustered in both space and time. These clustering features
are important for seismological studies. For the purpose of long term earthquake
prediction, such as zoning, earthquake potential estimation, people try to remove
the temporary clustering to estimate the background seismicity; on the other hand,for the short-term or real-time prediction, we need a good understanding of theearthquake clusters. Thus separating the background seismicity from the earthquake
clusters is considered to be of central importance.
A number of methods have been proposed for declustering a catalogue or separating the background and the clustering seismicity. Most of these methods remove
the earthquakes in a space-time window around a large event called the mainshock;the differences among these methods relate to the choice of the window sizes andsome other detailed procedures: for example, see Gardner and Knopoff, 1974; KeilisBorok and Kossobokov, 1986; Utsu, 1969. In general, the bigger is the magnitudeof the mainshock, the larger is the window. An alternative to the window-based
declustering methods is the link method based on a space-time distance between
1the events (e.g., Reseanberg, 1985; Frohlich and Davis, 1990; Davis and Frohlich,1991).
All of those conventional methods contain arbitrary parameters, whether for
defining the sizes of the aftershock window or the link distances, in both space andtime. Each choice of parameters gives a different declustered catalogue, and so a
different estimate of the background seismicity. A trouble with the conventional
declustering algorithms is thus to make an optimal choice of the parameters. Fora data set from a specific (rather narrow) seismic region, and a given threshold
magnitude, these choices are usually tuned by the user’s subjective impression of
the declustered output. Also, given a data set, the concept of what constitutes anaftershock has not been uniquely defined. These are among the reasons why somany declustering algorithms have been proposed.
Alternative to these deterministic declustering methods, ideas of probability
treatment to the background component and clustering component first appeared in
Kagan and Knoppoff (1978). The assumed a constant background occurrence rate
throughout the whole studied region, which is in fact location dependent as we cansee in this article. Zhuang et al. (2002, 2004) proposed the stochastic declustering
method and brought the probability based declustering method into practice. The
core of the stochastic declustering method is an iterative approach to simultaneously estimate the background intensity, assumed to be a function of space but notof time, and the parameters associated with clustering structures. Making use of
these estimates and the thinning operation, one can obtain the probabilities for each
event being a background event or a triggered event. These probabilities are the key
to realizing stochastic versions of the clustering family trees in the catalogue, andof course to separating the background events from the earthquake clusters.
Because these probabilities are estimated through a particular model, the closeness between the formulation of the model and the reality is the essentially importantfactor influencing the output. The closer is the model to the real data, the morereliable is the output. Of course, some model selection procedures can be used tochoose the best model among many models fitted to the same set of data. But theseprocedures usually give us a number indicating the overall goodness-of-fit for themodel, and rarely tell us whether there are some good points in a model if its overall
fit is worse than another. It is also difficult to find clues about how to improve theformulation of the clustering models through model selection procedures. In thispaper, we are going to show some graphical diagnostic methods useful for improving
model formulation.

2

The space-time ETAS model

In the empirical studies, the Omori law (Omori, 1898; and see, Utsu et al., 1995, fora review), describes the decay of aftershock frequencies with time, i.e.,n(t) =

K,(t + c)p(1)which n(t) is the occurrence rate of events at the time t after the occurrence of themainshock, K, c and p are constants. Another commonly accepted empirical law is2the Gutenberg-Ritcher law, describing the relationship between the magnitudes andoccurrence frequencies of earthquakes and taking the form oflog10 N (≥ M ) = a − bM,(2)i.e., the number of earthquake events decreases exponentially when we increase themagnitude threshold.
These above empirical laws has been considered as the foundation for statisticalmodelling in seismicity. Several point process models, in different forms (Ogata,
1988, 1998; Kagan, 1991; Musmeci and Vere-Jones, 1992; Rathbun, 1993; Console
and Marru, 2001; Console et al., 2003). Among them, Ogata (1998) suggested amodel with a component of Omori type (inverse power law) decay in time and acomponent for the locations of triggered events, independent from the time component. In general, all of those models classify the seismicity into two components, thebackground and the cluster, where each earthquake event, no matter it is a background event or generated by another event, produces (triggers) its own offspring(aftershocks) according to some branching rules.
All these models can be formulated in the form of the conditional intensity
function, i.e., the process is controlled by an intensity conditional on the observation
history (see, e.g., Daley and Vere-Jones, 2002, Chapter 7, for details of mathematicaltreatment on the conditional intensity), namely,
λ(t, x, y, M |Ht )dt dx dy dM = E[N (dt dx dy dM )|Ht ],(3)

where Ht is the observation history up to time t. In this study, we base our analysis
on the formulation of the space-time epidemic type aftershock sequence (ETAS)
models by Ogata (1998),
λ(t, x, y, M ) = λ(t, x, y|Ht )J(M )X
λ(t, x, y|Ht ) = µ(x, y) +κ(Mi )g(t − ti )f (x − xi , y − yi |Mi ),

(4)
(5)

i: ti <twhere

1. µ(x, y) is the background intensity, which is a function of space, but not oftime;
2. κ(M ) is the expected number of events triggered from an event of a magnitude
M , in the formκ(M ) = A exp[α(M − MC )],(6)
where A and α are constant, and MC is the magnitude threshold (see Utsu,
1969; Yamanaka and Shimazaki, 1990);
3. g(t) is the probability density function of the occurrence times of the triggered
events, taking the form

−p
t
p−1
,
(7)
1+
g(t) =
cci.e., the p.d.f form of the modified Omori law (Utsu, 1969);
3

4. f (x, y|M ) is the location distribution of the triggered events, which is formulated in either of two ways, a short range Gaussian decay and a long range
inverse power decay, or explicitly,
1x2 + y 2f (x, y|M ) =
exp −
(8)
2πDeα(M −MC )2Deα(M −MC )for the short-range decay, and(q − 1)D (q−1) eα(q−1)(M −MC )f (x, y|M ) =q ;π [x2 + y 2 + Deα(M −MC ) ]q > 1,(9)for the long-range decay; and5. J(M ) is the probability density of magnitudes of all the events, independent
from other components and taking the form of the Gentenburg-Richter law,J(M ) = βe−β(M −Mc ) , for M ≥ Mc

(10)

where β is linked with the Gutenberg-Richter’s b-value by β = b log 10 and Mc
is the magnitude threshold considered.
In the text below, we call an ETAS model equipped with (8) as Model I, andModel II if equipped with (9). The expected number of offspring that an event can
trigger, namely κ(M ), is also called its triggering ability.
In (6), (8) and (9), the spatial scaling factor for the direct aftershock region
is proportional to the triggering ability of the ancestor. This judgment is from
Kanamori and Anderson (1975) when they tried to explain the exponential law
for the number of aftershocks. In the coming section, we will use the stochastic
reconstruction method to verify whether it is a good choice.
Thinning procedure One of the key points of the stochastic declustering method
is the thinning operation on a point process determined by (5), which can split thewhole process into several sub-processes (Lewis and Shedler 1979; Ogata, 1981).
Suppose that events are numbered in time order from 1 to N . The probability ofan event j being triggered by the ith event is ζi (tj , xj , yj ) ,when j > i,
ρij =
(11)
λ(tj , xj , yj |Htj )0,otherwise,whereζi (t, x, y) = κ(Mi )g(t − ti )f (x − xi , y − yi |Mi ),(12)

represents the intensity triggered by the ith event. Moreover, the probability of theevent j being a triggered event is
X
ρj =
ρij ,
(13)
i<j

4

and the probability that the jth event belongs to the background is
ϕj = 1 − ρ j =µ(xj , yj ).
λ(tj , xj , yj |Htj )(14)

If we select each event j with probabilities ρij , ρj or ϕj , we can form up a newprocess being the triggered process by the ith event, the clustering process or the
background process, respectively.
Variable kernel estimates Once a background process is obtained, we can estimate the background intensity by some smoothing techniques. Rather than repeat
the thinning procedure and the kernel estimation procedure many times to get an
average estimate of the background intensity, we directly estimate the average by
weighting, i.e.,
1Xµ̂(x, y) =ϕi Zhi (x − xi , y − yi ),(15)
T i
where i runs over all of the events in the whole process, T is the length of the time
period of the process, and Zh is the Gaussian Kernel function with a bandwidth h.
The variable bandwidth hj is determined by
hj = max{, inf(r : N [B(xi , yi ; r)] > np )},(16)

where  is a small real number, B(x, y; r) is the disk centered at (x, y) with a radiusof r, and np is positive integer, i.e., hj is the distance to np th closest event.
Maximum likelihood estimates Given an estimated intensity function u(x, y),we setµ(x, y) = νu(x, y)(17)
for the background rate, where ν is a positive-value parameter, and then obtain the
maximum likelihood estimates, θ̂ = (ν̂, Â, α̂, ĉ, p̂, D̂) for Mode I and θ̂ = (ν̂, Â, α̂, ĉ, p̂, D̂, q̂)for Mode II, by maximizing the log likelihood
Z T ZZX
log L(θ) =log λθ (tk , xk , yk | Htk ) −λθ (t, x, y|Ht )dxdydt
(18)
0

k

S

where k runs over all the events in the study regions S and the study time interval[0, T ] (NOT necessary all the events observed).
Thus, the background seismicity and the parameters in the clustering structures
simultaneously can be estimated by an iterative approach (Zhuang et al., 2002).
We first assume an initial rough background intensity, calculate the MLE of the
branching parameters, calculate ϕj , then use these ϕj to estimate a better evaluationof the background intensity. Replacing the initial background intensity by the newly
estimated one, we repeat these procedures until the results convergences.

5

(b)

45
40
35

latitude

25

30

35
25

30

latitude

40

45

(a)

125

130

135

140

145

150

155

0

5000

longitude

10000

15000

20000

25000

time in days

Figure 1: Seismicity in the Japan region and nearby during 1926–1999 (MJ ≥ 4.2).
(a) Epicenter locations; (b) Latitudes of epicenter locations against occurrence times.
The shaded region represents the study space-time range.

3

Data and preliminary results

Two datasets are considered in this study. The first dataset is the JMA catalogue.
We select the events in the ranges of longitude 121◦ – 155◦ E, latitude: 21◦ – 48◦ N ,depth 0 – 100 km, time 1926/Jan/1 – 1999/Dec/31 and magnitude ≥ MJ 4.2. There
are 19,139 events in this data set.
As other earthquake catalogues covering a long history, completeness and homogeneity are always problems causing troubles for statistical analysis. The incompleteness of the early period and inhomogeneity in this data set can be easily seenfrom Figures 1. To tackle these problems, we choose a study space–time range, inwhich the seismicity seems to be relatively complete and homogeneous, with rangesof longitude 130◦ – 146◦ E, latitude 33◦ – 42.5◦ N , a time period of 10,000 – 26,814
days after 1926/Jan/1, and same depth and magnitude ranges as the whole data
set. There are 8283 events occurring in the study space–time range.
The earthquake events not in the target region and time period but in the dataset, namely the boundary catalogue, are used as the boundary effects in the estimation of the clustering parameters and the background rate in the model, i.e., these
events are not counted in computing the likelihood function (18), but are put in theobservation history when evaluating the conditional intensity (5).
The second dataset is the Taiwan CWB (Central Weather Bureau) catalogue,
we take the events of M ≥ 5.3 in the period from 1941-1-1 (14975 days from 19001-1) to 2001-12-31 and in the region (120∼123◦ E, 21.5∼25.2◦ N) for estimating theparameters, non-target events with a magnitude no less than 5.3 in the catalogue
are used as boundary events, i.e., they do not directly contribute to the likelihood
function, but influence the conditional intensity function. According to the local
structures, we then divide each of these three main subregions into smaller ones, as
shown in Figure 2.
The background intensity of the study region of JMA catalogue is shown in
6

25

NC
NW

NE1

24

NE2

NE3

CW
SW
23

Latitude

CC
SE1

SC

21

22

SE2

120.0

120.5

121.0

121.5

122.0

122.5

123.0

Longitude

Figure 2: Epicenter locations of earthquakes in the Taiwan region during 1900–2000(M ≥ 5.3) and subdivision. Sizes of the circle indicate magnitudes from 5.3 to 8.2.

7Figure 3. As a feature complementary to the background, we defineω(x, y) = 1 −µ(x, y)Λ(x, y)(19)

as the clustering coefficient to measure the clustering effect relative to the total
intensity at a location (x, y). This function is as important as µ(x, y) in earthquake
hazard estimation, for damages caused by aftershocks cannot be ignored.

4

Verifying stationarity of the backgroundIn the ETAS model, we assume that the background process is a stationary one.
It is not so easy to test the stationarity of this background seismicity by using
conventional declustering methods because of the difficulties in choosing declusteringparameters. Using the thinning method, we can easily test whether the backgroundseismicity rate changes with the time in a specified region, by constructing thefunction
X
S(t) =
ϕi ,
(20)
ti <t

where the subscript i runs over all of the events in this region. The above function
S(t) represents the cumulative number of events in the background process in theregion considered.
The results of S(t) for some subregions of the Taiwan region are in shown Figure 4, we can see that the earthquake clusters in each subregion are removed, whichmakes it easier for us to study the background seismicity. A clear feature of the background seismicity in these subregions is the quiescence before the Chichi earthquake
occurring subregion CW in 1999. Including Subregion CW itself, all the neighboring
subregions show quiescence in seismicity from early 1960 to the occurrence time of
the Chichi earthquake. On the other hand, the eastern subregions do not show sucha clear quiescence in seismicity.
In summary, the background seismicity seems not to be stationary, and its occurrence rate usually changes before or after great earthquakes.

5

Residual Analysis

The goodness-of-fit of the space-time ETAS model to the earthquake data can beconveniently investigated by using residual analysis. There have been several methods for carrying out the residual analysis. For example, Ogata (2003) re-fitted thedata with a time-variant version of the ETAS model and then took the ratio between
the conditional intensities of the two models as the residual process. Schoenberg(2003) proposed the thinning residuals, i.e., to keep each event i with probability min λ(t, x, y|Ht )/λ(ti , xi , yi |Hti ) to obtain a homogeneous Poisson process. Thismethod is the essentially the same as Stoyan and Grabarnik (1991). More general
discussions on transforming a point process into a Poisson process can be found in
Schoenberg (1991), Vere-Jones and Schoenberg (2004).

8

(a)

42

103
40

100

10
38
1
36

0.1

0.01

34

130

135

140

145

(b)

42

1.0

0.8
40

0.6
38

0.4
36

0.2
34

130

135

140

145

Figure 3: (a) Estimated background rate µ̂ in Equation (15) in the study region(unit: events/(degree2 ×74 years)) and (b) the clustering coefficient in Equation
(19).

9

1920

(a) NW

1900

1920(d) NC

1900

1920(g) NE1

1900

1940

1960

1960

Time

1940

1960

Time

1940

1980

1980

1980

2000

2000

2000

1920

(b) CW

1900

1920(e) CC

1900

1920(h) NE2

1900

1960

Time

1940

1960

1960

Time

1940

1940

1980

1980

1980

2000

2000

2000

20
15
10
5
0
7

Cummulative #Magnitude
Cummulative #

6
5
12
10
8
6
4
2
0
7
6
5
60
50
40

8
6
4
0

2
7
6
5
100
80

0

10

20

Magnitude

5

10

1920

(c) SW

1900

1920(f) SC

1900

1920(i) SE1

1900

1960

Time

1940

1960

Time

1940

1960

Time

1940

1980

1980

1980

2000

2000

2000

Figure 4: Cumulative background seismicity S(t) in (20) (dark step lines) and clustering seismicity #(t) − S(t) (gray step lines) in some Taiwan subregions. The
magnitudes against the occurrence times of the events are plotted in the lower partof each panel.

Time

7

Time

40

6

20
5

6

0

30

25
20
15
10
0
5

6

7

5

14
12
10

7

60

Magnitude
Cummulative #Cummulative #Magnitude
Cummulative #Magnitude
Cummulative #Magnitude

Cummulative #Magnitude

8
6
4
0

2
7
5

Cummulative #

6
5
25
20
0

Magnitude

15
10
0

5
7
6
5
20
15
10
7
6
5

Cummulative #Magnitude

Following Stoyan and Grabarnik (1991), to each event, we attach a weight
wi =

1
,
λ(ti , xi , yi |Hti )then the total weights wi in a space-time volume B has expectation
XE[wi ] = |B|,

(21)

(22)

(ti ,xi ,yi )∈B

where |B| denotes the volume of B.
From the residual analysis (Figure 5), we can that there are some deviation in theresidual of the data to the ETAS model. In Figure 5(a), the slope of the cumulative
wights against the time axis changes at the years around 1940 and 1975, or evenaround 1994. These changes mainly correspond to the changes of the observation
system. On the latitude axis, the cumulative weights change the accumulation rate,showing that this seismic region can be basically divided into three subregion fromsouth to north; while along the longitude axis, the accumulation rate of the weights
changes around 120.5, 121.5 and 122.5, showing the difference in seismicity from thewest to the east. All these results indicating the difference in clustering seismicity
caused by local tectonic environments.

6

Verifying model formulation by stochastic reconstruction

As mentioned in Section 2, if we select each event j with probabilities ρij , ρj or
ϕj , we can get a new process being the triggered process by the ith event, theclustering process or the background processes, respectively. That is to say, we
can separate the whole catalogue into different family trees. This method tackles
the difficulties in testing hypotheses associated with earthquake clustering features,
which are caused by the complicated overlapping of the background seismicity anddifferent earthquake clusters in both space and time. We can repeat such thinning operations for many times to get different stochastic versions of separationsof the earthquake clusters. The non-uniqueness of such realizations illustrates theuncertainty in determining earthquake clusters, and thus repetition can help us toevaluate the significance of some properties of seismicity clustering patterns. However, we can also implement these tests by working with the probabilities ϕj andρij directly. In the coming sections, we will show how to use these probabilities toreconstruct the characteristics associated with earthquake clustering features, using
the Japanese JMA catalogue as an example.

6.1

Location distributions

Define the standardized distance between a triggered event j and its direct ancestor,
assumed i, bys(xj − xi )2 + (yj − yi )2
rij =
.
(23)
D exp[α(Mi − Mc )]11(b)2500001500000

50000

Cumumlative residuals

15000050000
0Cumumlative residuals

250000

(a)

1925

1940

1955

1970

1985

2000

Time (year)

21.5

22.0

22.5

23.0

23.5

24.0

24.5

Latitude15000050000
0Cumumlative residuals

250000

(c)

120.0

120.5

121.0

121.5

122.0

122.5

123.0

Longitude

Figure 5: Results from residual analysis for the Taiwan catalogue. The cumulative
weights against times, latitudes and longitudes are plotted in (a), (b), and (c),respectively.

12

25.0

1.0

MODEL 2

0.6
0.0

0.2

0.4

Prob. Density

0.6
0.4
0.0

0.2

Prob. Density

(b)

0.8

(a)

0.8

1.0

MODEL 1

0

2

4

6

8

10

0

Standardized distance

2

4

6

8

10

Standardized distance

Figure 6: Reconstruction results for the distribution of the standardized triggering
distances fˆR (r) in (26) (gray circles) by using Model I (a) and Model II (b) for the
JMA catalogue. The theoretical curves of fR (r) in (24) and (25) are plotted in solidlines in (a) and (b), respectively.
From (8) and (9), rij has a density function of2fR (r) = 2re−r , r ≥ 0;
for Model I andfR (r) =2r(q − 1), r ≥ 0;(1 + r 2 )q

(24)

(25)

for Model II, respectively. The distribution with a density of (24) is called a Rayleighdistribution. On the other hand, fR (r) can be reconstructed throughPi,j ρij I(|rij − r| < ∆r/2)
P,(26)
fˆR (r) =∆r i,j ρij

where ∆r is a small positive number, and I is the index function such that

1, if the logical statement x holds,
I(x) =
0, else.

(27)

The comparison between fˆR and fR for the two models are shown in Figure 6.
It can be seen that, if Model I is used, the reconstructed probability density of thetransformed distances between the ancestors and the direct offspring is quite different
from the theoretical one. When Model II is used, the reconstructed probability
density is very close to the theoretical one. These results confirm that the aftershocks
decay in a long range in space rather than a short range (Ogata, 1998; Console andMurru, 2001; Console et al. 2003). These results also imply the robustness of thereconstruction method, for we can get a reconstructed probability density function
very close to the corresponding function in Model II, even if an improper model like
Model I is employed.
Since Model II fits the seismicity much better than Model I, we only consider
Model II for reconstruction in the following sections.
13

50
10 20
2

5

(b)

1

Triggering abilities

All events
Triggered
Background
Theoretical

0.2

0.5

50
0.2

0.5

1

2

5

10 20

(a)

0.05

0.05

Triggering abilities

All events
Triggered
Background
Theoretical

5

6
Magnitude

7

8

5

6
Magnitude

7

8

Figure 7: Reconstruction of the triggering abilities, κ̂b (M ) in (29) of the background
events and κ̂t (M ) in (30) of the triggered events for (a) the JMA catalogue and (b)the simulated catalogue. For comparison, the empirical functions of the triggering
abilities, κ̂(M ) in (28) for all the events are plotted in gray circles, and the corresponding theoretical functions, κ(M ) = Aeα(M −MC ) , are represented by the straightlines.

6.2

Difference in triggering ability between the background
events and the triggered events

The triggering abilities of an event sized M from all the events, the background
events and the triggered events for the JMA catalogue and the simulated catalogue
can be reconstructed by using
P P
ρij I(|Mi − M | < ∆M/2)iPjκ̂(M ) =(28)
i I(|Mi − M | < ∆M/2)
P P
ϕi ρij I(|Mi − M | < ∆M/2)iPjκ̂b (M ) =(29)i ϕi I(|Mi − M | < ∆M/2)

and

P P
(1 − ϕi )ρij I(|Mi − M | < ∆M/2)iPj,κ̂t (M ) =i (1 − ϕi )I(|Mi − M | < ∆M/2)(30)

respectively, as shown in Figure 7. Both the background events and the triggered
events generate offspring approximately according to different exponential laws. For
the same ancestor magnitude, a triggered event generates more offspring than abackground event. The higher is the magnitude, the smaller is the difference. The
results for the simulated catalogue show that there is no differences in triggering
abilities between these two types of events, indicating that these differences are notcaused by numerical procedures.
In the ETAS model, the background events and triggered events generate offspring in the same way. The reconstruction results show that we should have different exponential laws for the triggering abilities of these two types of event, which is14

a possible direction on improving the current ETAS models, i.e., to treat the background events and triggered events separately in the model. The higher triggering
abilities of the triggered events may be because they occur in an environment where
the stress fields is adjusting to the stress changes caused by their ancestors.

6.3

Distributions of offspring locations from different magnitude classes

Because of the historical reasons mentioned in Section 2, we model the distributionof locations of the direct offspring from an earthquake as an inverse power one with ascaling factor associated with the ancestor’s magnitude, i.e., Deα(M −MC ) . Immediate
questions for such a choice are:(a) Is this scaling factor necessary? i.e., Can we use a constant D0 to replace thescaling factor Deα(M −MC ) in the model?
(b) Is it necessary to link the scaling factor to the triggering ability, or should we
introduce a new parameter γ instead of α for the scaling factor?
To answer the above questions, for a small interval M of magnitudes, we selectthe pairs {(i, j)} such that Mi ∈ M and then estimate the scaling factor DM for
M, in a way similarly to the estimation of the parameters in mixture models byusing the expectation-maximization algorithm (see, e.g., Eggermont and LaRiccia,
2001, Section 2.4). Given ρij , consider the following log-likelihood function,X X(q − 1)D (q−1) Rij.
(31)
log L =
ρij log
2
(Rij
+ D)q
M ∈M ji

To maximize it, let∂ log L= 0,
∂D D=DM
i.e.,

X X
ρij
q−1 X X= 0.
ρij − q2
DM M ∈M j
R
+
D
M
ij
M ∈M ji(32)

i

Thus, we can construct the following iteration to solve the above equationPP(q − 1) Mi ∈M j ρij(n+1)
DM = X X
.
ρij
q
(n)
2
Mi ∈M j Rij + DM

(33)

Figure 8 shows the values of DM against the magnitude classes. We can see thatvalues of DM have a slope different from κ(M ), which is 0.5008 by the least squarefit. Thus, it is not suitable use κ(M ) as the scaling factor, a better choice is tointroduce another parameter γ as the coefficient in the exponential part .

15

0.1
0.001

0.01

Scaling Factor

0.1
0.01

Scaling Factor

0.001

5

6

7

8

5

Magnitude

6

7

8

Magnitude

Figure 8: Re-estimated DM for the JMA catalogue (a) and the simulated catalogue.
Theoretical fitting curves, Deα(M −MC ) , are represented by the straight lines.

7

Conclusion

In this paper, we discussed using the ETAS model to descibe the clustering featuresof earthquake processes and some methods on how to evaluate the goodness-of-fit.
The results show that the ETAS model is a fairly good starting point for modellingthe earthquake clusters. Then using different diagnostic methods, we tested thestationary of the background process, the overall goodness-of-fit and each clustering
component in the model formulation. These methods are not only useful in justifying
the models fitted to the data, but also helpful for us to better understand the data.

References[1] Console, R., and M. Murru, A simple and testable model for earthquake clustering, Journal of Geophysical Research, 106, 8699-8711, 2001.
[2] Console, R., M. Murru, and A. M. Lombardi, Refining earthquake
clustering models. Journal of Geophysical Research, 108(B10), 2468,
doi:10.1029/2002JB002130, 2003.
[3] Daley, D., and D. Vere-Jones, An Introduction to the Theory of Point Process,Vol. 1, Springer, New York, 2002.
[4] Davis, S. D., and C. Frohlich, Single-link cluster analysis, synthetic earthquake
catalogs, and aftershock identification, Geophysical Journal International, 104,
289–306, 1991.
[5] Eggermont, P.P.B., and V.N. LaRiccia, Maximum Penalized Likelihood Estimation, Vol I: Density Estimation, Springer, New York, 2001.

16

[6] Frohlich, C., and Davis, S. D. Single-link cluster analysis as a method to evaluate spatial and temporal properties of earthquake catalogs, Geophysical JournalInternational, 100, 19–32, 100, 1990.
[7] Gardner, J., and L. Knopoff, Is the sequence of earthquakes in southern California, with aftershock removed, Poissonian? Bulletin of the Seismological society ofAmerica, 64, 1363 – 1367, 1974.
[8] Kagan, Y. Likelihood analysis of earthquake catalogues,” Journal of GeophysicalResearch, 106(B7), 135–148, 1991.
[9] Kagan, Y. Y., and L. Knopoff Statistical study of the occurrence of shallow
earthquakes, Geophys. J.Roy. astr. Soc., 55, 67–86, 1978.
[10] Kanamori H., and Anderson D., Theoretical basis of some empirical relationsin Seismology, Bulletin of Seismological Society of America, 65, 1073–1095, 1975.
[11] Keilis-Borok, V. I., and V. I. Kossobokov, Time of increased probability for thegreat earthquakes of the world,” Computational Seismology, 19, 45–58, 1986.
[12] Lewis, P. A. W., and E. Shedler, Simulation of non-homogeneous Poisson processes by thinning, Naval Res. Logistics Quart., 26, 403 – 413, 1979.
[13] Musmeci, F., and D. Vere-Jones, A variable-grid algorithm for smoothing clustered data, Biometrics, 42, 483–494, 1986.
[14] Musmeci, F., and D. Vere-Jones, A space-time clustering model for historicalearthquakes, Annals of the Institute of Statistical Mathematics, 44, 1–11, 1992.
[15] Ogata, Y., On Lewis’ simulation method for point processes, IEEE translations
on information theory, IT-27, 23–31, 1981.
[16] Ogata, Y., Statistical models for earthquake occurrences and residual analysis
for point processes, Journal of American Statistical Association, 83, 9–27, 1988.
[17] Ogata, Y., Space-time point-process models for earthquake occurrences, Annalsof the Institute of Statistical Mathematics, 50, 379–402, 1998.
[18] Ogata, Y., K. Katsura, and M. Tanemura, Modelling heterogeneous space-time
occurrences of earthquake and its residual analysis, Journal of the Royal StatisticalSociety: Applied Statistics, 52(Part 4), 499–509, 2003.
[19] Omori, F., On after-shocks of earthquakes, J. Coll. Sci. Imp. Univ. Tokyo, 7,
111–200, 1894.
[20] Rathbun, S. L., Modeling marked spatio-temporal point patterns, Bulletin of
the International Statistical Institute, 55, Book 2, 379–396, 1993.
[21] Reasenberg, P., Second-order moment of Central California Seismicity, 1969 –
1982, Journal of Geophysical Research, 90(B7), 5479–5495, 1985.

17

[22] Schoenberg, F. P., Multi-dimensional residual analysis of point process models
for earthquake occurrences, To appear on Journal of the American Statistical
Association, 2003.
[23] Stoyan D. and Grabarnik P., Second-order characteristics for stochastic structures connected with Gibbs point processes. Mathematische Nachritchten, 151,
95-100, 1991.
[24] Utsu, T., Aftershock and earthquake statistics (I): Some parameters which characterize an aftershock sequence and their interrelations,” Journal of the Facultyof Science, Hokkaido University, 3, Ser. VII (Geophysics), 129–195, 1969.
[25] Utsu, T., Y. Ogata and R. S. Matsu’ura, The centenary of the Omori formulafor a decay law of aftershock activity, J. Phys. Earth, 43, 1–33, 1995.
[26] Vere-Jones, D. and Schoenberg, F.P., Rescaling marked point processes. Australian and New Zealand Journal of Statistics, to appear, 2004.
[27] Yamanaka, Y. and K. Shimazaki, Scaling relationship between the number ofaftershocks and the size of the main shock, J. Phys. Earth, 38, 305–324, 1990.
[28] Zhuang, J., Y. Ogata and D. Vere-Jones, Stochastic declustering of space-timeearthquake occurrences, Journal of the American Statistical Association, 97, 369–
380, 2002.
[29] Zhuang J., Ogata Y. and Vere-Jones D., Analyzing earthquake clustering features by using stochastic reconstruction. Journal of Geophysical Research, to appear, 2004.

18