Research in Geophysics 2012; volume 2:e8

Long-term earthquake forecasts based on the epidemic-type aftershock sequence (ETAS) model for
short-term clustering
Jiancang Zhuang
Institute of Statistical Mathematics and Department of Statistical Science, the Graduate University
for Advanced Studies, Tokyo, Japan

ly

on

e

us

er

ci

Based on the ETAS (epidemic-type aftershock sequence) model,
which is used for describing the features of short-term clustering ofearthquake occurrence, this paper presents some theories and techniques related to evaluating the probability distribution of the maximum magnitude in a given space-time window, where theGutenberg-Richter law for earthquake magnitude distribution cannot
be directly applied. It is seen that the distribution of the maximum
magnitude in a given space-time volume is determined in the longterm by the background seismicity rate and the magnitude distribution of the largest events in each earthquake cluster. The techniques
introduced were applied to the seismicity in the Japan region in the
period from 1926 to 2009. It was found that the regions most likely tohave big earthquakes are along the Tohoku (northeastern Japan) Arc
and the Kuril Arc, both with much higher probabilities than the offshore Nankai and Tokai regions.

ity, for which removal of temporal cluster members is considered to
be of central importance. Both the background rate and the clustering structures are important to earthquake hazard estimation,
because the background rate tells us the potential risk of the occurrence of an earthquake cluster, and the clustering structure determines the size of the earthquake cluster and the magnitude of thebiggest event in a cluster.
In the literature of seismology, people traditionally use some window based methods1-3 or link-based methods4-6 to decluster the catalog or to identify earthquake clusters. The parameters in thesedeclustering rules depend to a large degree on the intuitions andexperience of researchers, sometimes even colored by the user’s subjective impression of the declustered output. The space-time epidemic-type aftershock sequence (ETAS) model describes the features ofearthquake clusters. It has become a standard model for testing
hypotheses and a starting point for short-term earthquake forecasts.7-12 Based on this, Zhuang et al.13 proposed the stochastic
declustering method to separate seismicity into a background component and a clustering component in a probability manner, following the initial ideas by Kagan and Knopoff.14 This method can also
produce stochastic versions of declustered catalogs similar to the
output of the conventional earthquake declustering methods.
The aim of this article is to show how to use the background rate
and the clustering structure to estimate the probability of future
large earthquakes in a particular region. A brief description of the
ETAS model is given and its use to separate the background rate fromoverall seismicity is discussed. Problems in estimating the background rate and the parameters of the clustering structure are examined. Formulas are then elaborated to evaluate the probabilities
related to the magnitude of the largest descendant from an initial
event with a given magnitude. Based on these probabilities, the probability density of the magnitude distribution of the largest earthquake in an arbitrary earthquake cluster is calculated. In order toillustrate this, we use the approach set out in the Japan
Meteorological Agency (JMA) catalog to estimate earthquake risk in
the Japan region.

al

AbstractommIntroduction

N

on

-c

The Poisson model, or other renewal models, are often used toevaluate the long-term risk of large earthquakes in a particularregion. However, clusters of aftershocks complicate the statisticalanalysis that evaluates the background seismic activity. The clustering features differ from place to place. To forecast the location of thelarge earthquakes, it is necessary to analyze the background seismic-

Correspondence: Jiancang Zhuang, Institute of Statistical Mathematics andDepartment of Statistical Science, the Graduate University for Advanced
Studies, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan.
E-mail: zhuangjc@ism.ac.jp
Key words: seismicity, statistical methods, background seismicity, ETAS
model, maximum magnitude.
Received for publication: 15 February 2012.
Revision received: 11 April 2012.
Accepted for publication: 18 April 2012.
This work is licensed under a Creative Commons Attribution
NonCommercial 3.0 License (CC BY-NC 3.0).
©Copyright J. Zhuang, 2012
Licensee PAGEPress, Italy
Research in Geophysics 2012; 2:e8
doi:10.4081/rg.2012.e8[page 52]

Definition of the epidemic-type aftershock
sequence model
In the ETAS model, there is no difference in triggering seismicity
among foreshocks, mainshocks and aftershocks. Each event triggers
its own offspring independently according to some probability rules.
In brief, the time-varying seismicity rate (mathematically termed asconditional intensity) of this model takes the form of:[Research in Geophysics 2012; 2:e8]

Articleus

represents the probability density functions for the occurrence times
and locations of direct offspring, respectively. In equations (2), (3) and(4), A, a, c, p, D, q and g are constant parameters.
Equation (3) is the probability density form of the Omori-Utsu formula.15-17 The form of equation (4) is based on the analysis in Zhuang
et al.18 and Ogata and Zhuang,19 which are improved versions of thosein Ogata.20 Many other forms can also be found in papers publishedover last 20 years.8,21

where i runs over all of the events in the whole process, T is the lengthof the time period of the process, and Zd is the Gaussian Kernel function
with a bandwidth d. The variable dj represents the varying bandwidth
calculated for each event j in the following way. Given a suitable integer
np, find the smallest disk centered at the location of the jth event which
includes at least np other earthquakes and whose radius is larger thansome minimum value (e.g. a distance within 0.02 degrees, which is theorder of the location error), d, and let this minimum radius be dj.
To find optimal values of np and d, the above variable kernel estimates
were applied to estimate the rate of a simulated inhomogeneous
Poisson point process. Through cross-validation,30 it is found that the
optimal np is in 3 ∼ 6. The parameter d does not influence the prediction
too much, as long as the locations of points are not rounded to a certain
precision. This parameter only becomes important when some points
happen to overlap at one location caused by rounding the number.

ly

and

on

is the mean number of direct offspring from an event sized m, mc being
the magnitude threshold, andprovides an estimate of the probability that the ith event is a background
event, since, if we keep the ith event with probability ϕi for all the eventsin the process, we can obtain a process with the occurrence rate of m(x,y).28,29 Thus, the background intensity can be estimated by applying
some smoothing techniques to the background catalog. Or, more efficiently, we can directly estimate the average by weighting all the eventswith their corresponding background probabilities. In this study, thevariable kernel estimates13 were chosen, i.e.

e

where
a) m(x, y) represents the background seismicity rate; b) the term κ(mi)g(t − ti) f(x − xi, y − yi, mi) is the contribution to seismicity rate by theith event occurring previously, andalModel estimationci

Iterative algorithm

To estimate the model parameters and the background seismicity rate
simultaneously,10,13 the following iterative algorithm can be introduced.
A1. Given a fixed np and d, (5 and 0.05 degree, equivalent to 5.56 km on
the earth’s surface, which is close to the location error of earthquakes),
calculate the bandwidth hj for each event (tj, xj, mj), j = 1, 2, ··· , N.
A2. Set ℓ←0, and u(0)(x, y)←1.
A3. Using the maximum likelihood procedure, fit the model with conditional intensity functionomMaximum likelihood estimatesmer

To forecast seismicity by using the model specified by equation (1),
we need to solve the following technical problems: i) estimating background seismicity rate; and ii) estimating the model parameters, (A, a,c, p, D, q, g).

N

on

-c

For an observation of the process, recorded as a list in the form of {(ti,xi, yi, mi): i=1, ···, N}, from a spatial region S and a time interval [0,T],
the likelihood has the standard form.22

If the background seismicity rate m(x, y) is known, the model parameters, (A, a, c, p, D, q, g), can be estimated through maximizing the likelihood function. In computation, the Davidon-Fletch-Powell method23 isused to minimize – ln L.

Estimating background seismicity
Many different approaches to assess background seismicity rate havealso been used: i) using a rate proportional to total seismicity rate of all
events or only of the big events in the catalog;20,24 ii) using a declustering method to decluster the catalog and use the total rate in the declustered catalog as background rate;20,25,26 iii) weighting each event by an
estimated probability that it is a background event;10,13 and iv) using themethod introduced by Ogata,27 which is a Bayesian smoothness prior on
a tessellation grids to estimate the spatial variation of the backgroundand the model parameter at the same time.
In this study, the third method is used because it is relatively simple
and gives an unbiased estimate of the intensity function. First of all,once the conditional intensity function is estimated, the contribution of
the spontaneous seismicity rate at the occurrence of the ith event,to the earthquake data, where κ, g and f are defined in equations (4),(5) and (6), and ν is the relaxing coefficient, which is introduced inorder to fasten the convergency speed of this algorithm. The model
parameters are (ν, A, a, c, p, D, q, g).
A4. Calculate ϕj for j = 1, 2, ··· , N by using equation (6).
A5. Calculate m(x, y) by using equation (7) and record it as u(ℓ+1)(x, y).
A6. If max|u(ℓ+1)(x, y)− u(ℓ)(x, y) | >ε, where ε is a given small positive
number, then set ℓ←ℓ+1 and go to step A3; otherwise, take νu(ℓ+1)(x, y)as the background rate and also output ϕi.

Theoretical distributions associated withthe largest magnitude of all the descendants
from a given event
The method for obtaining the background rate is based on the spacetime ETAS model, i.e. the occurrence rate of the clusters. However, thehazard caused by earthquake clusters not only depends on their occurrence rate, but also on the magnitudes of large events in each cluster.

[Research in Geophysics 2012; 2:e8][page 53]

Article
The largest magnitude in an arbitrary cluster and in a given space-time
volume must be evaluated.

Basic equationsFirst, assume that the probability density function (p.d.f), s(m), for all
the magnitudes obeys the Gutenberg-Richter law (exponential distribution), i.e.

equation (13). Considering that the background ground rate m(x, y)
gives the occurrence rate of earthquake clusters, a combination of F(m)and the background rate provides information on the maximum magnitude of future earthquakes.

Asymptotical properties of F(m)Because F(m) is a function determined by κ(m) and s(m), it isinfluenced only by the parameters A, a and β in the model. Substituting
equations (9) and (2) into equation (13) we obtain

where β is the constant linked with the Gutenberg-Richter’s b-value by
β=b ln 10 and mc is the magnitude threshold.
Note the number of children from an event of magnitude m is a
Poisson random variable with a mean of κ(m), i.e.
Pr {an event sized m has n direct offspring} =only

The probability that an event of magnitude m has no offspring greaterthan m¢, which is the same as the probability that the largest descendant from an event of magnitude m is smaller than m¢, can be derived
from the model analytically31 in the following way:

whereis the complementary incomplete gamma function. An analysis with moment generating functions32 also gave similar
forms of the above equations. But they did not consider the case inwhich the process is supercritical. The properties of F(m) are closely
related to the criticality of the process.31 According to the Appendix (seealso Zhuang et al.33), the critical parameter ρ is determined bye

Substituting equations (2) and (9) into the above equation,us

it is evident that β>a is required. According to the discussion in the
appendix, F has the following properties:31,34

ciali) subcritical case ( <1) where each family tree ultimately dies off andthe whole process is stable and stationary. When m is large enough,mer

It is evident that ζ(m, m’) has the form

om

where

and when a≤β/2

N

on

-c

represents the probability that the largest earthquake in an arbitrarycluster, including the initial event and all its descendants, is greaterthan m¢. Thus, F(m) gives the largest magnitude in an arbitrary cluster.
In other words, the magnitude distribution of the largest event in a cluster is no longer Gutenberg-Richter law, but determined by functionalii) critical case ( =1), where each family tree dies off with a long tailand the population of the whole process in unit time experiences anunlimited increase. In this case, the asymptotic property of F branches
according to the relative values of a and β.34 When β/2<a<β,Figure 1. Influence of A, a and β to function F(m). (A) Parameter A changes, but a=1.2 and β=2.4 are fixed. (B) Parameter a changes,
but A=0.3 and β=2.4 are fixed. (C) Parameter β changes, but A=0.3 and a=1.2 are fixed. The thin solid, thick solid and dashed curves
represent the subcritical, critical and supercritical regimes, respectively.
[page 54]

[Research in Geophysics 2012; 2:e8]

Articleiii) supercritical case ( ≥1), where some of the family trees may neverdie off and the population of the whole process explodes, F(m) ∼ C,
where C is linked with the extinct probability Pc(m) of the family starting from an event of magnitude m by Pc(m)=e–Cκ(m).
Using the following iteration algorithm, F(m) can be solved for eachvalue of m

In Figure 1, first two of A, a and β are fixed and then the other ischanged to see how F(m) changes. As discussed above, the influence of
all the three parameters, a, β and A, to F(m) is through their influence
to the criticality, i.e. the properties of F(m) can be divided into three
regimes: subcritical, critical and supercritical, where F(m) ∼ s(m), F(m)∼ 1/κ(m) and F(m) ∼ C, respectively.

Probabilities of a future large earthquake

m

er

ci

al

us

e

on

ly

It should be remembered that F(m) is the survivor probability function of the magnitude of the biggest event in a cluster, including the initial event and all its descendants. Given a space-time window V and amagnitude m, the probability that the largest event in V is greater thanm can be evaluated in the following way:

N

on

-c

om

where
is the mean number of clusters occurring in V.
Here it is assumed that all events of a cluster occur at the same location and at the same time, which may not be true. However, for longterm risk evaluation, such a simple approximation is reasonable andalso reduces the calculation tasks.
It should be noted that the above equation only holds true when theprocess is stationary, i.e. subcritical. In the case of supercritical or critical processes, the population of the process explodes to infinity. Except
for the trivial case of Λ(V)=0, the largest event in V can theoreticallyexceed any magnitude.

Results
To illustrate the above procedures, the JMA (Japan Meteorological
Agency) earthquake catalog was used to evaluate the example of thelong-term probabilities of large earthquakes in and around the Japanregion. First, data were selected from the JMA catalog in the ranges oflongitude 128°-148°E, latitude 30°-45°N, depth 0-100 km, over the period from 1 January 1926 to 31 December 2009, and magnitude MJ≥4.0.
There are 36,144 events in this dataset. Figure 2 shows the locations of
the epicenters and occurrence times of events in the dataset.
Before fitting the space-time ETAS model to the catalog, some preprocessing of data is necessary. For the earthquake catalog as in this study,
which covers records of a time period of 38 years, completeness andhomogeneity are usually problems for statistical analysis. If seismicity
in some regions or the whole region has an increasing trend, the fitting

Figure 2. (A) Epicenter map of the shallow earthquakes (MJ≥4.0
and depth ≤100, 1 January 1926 to 22 September 2003) from the
JMA catalog. (B) A space-time plot of latitudes against occurrence
times, and (C) a space-time plot of occurrence times against longitudes. The polygon in (A) represents the target region used in
the likelihood computation, and the black dots and red dots represent the target events and complementary events, respectively.

[Research in Geophysics 2012; 2:e8][page 55]

Article

Figure 4. Spatial variations of b-values estimated by variable kernel estimates.

m

er

ci

al

us

e

on

ly

Figure 3. Background seismicity rate (occurrence rate of earthquake clusters) [events/(day·deg2)] of earthquakes with magnitude
MJ≥4.0, obtained based on the space-time ETAS model.

-comFigure 5. Spatial variations of the occurrence probabilities that at least one earthquake in a 1°¥1°¥ l-year space-time volume is greaterthan MJ 5.0 (left panel), MJ 6.0 (middle panel) and MJ 7.0 (right panel).

Nonresults often do not converge or they converge to some unreasonablevalues. Thus, a space-time range was chosen, in which the seismicity
seems to be relatively and visually complete and homogeneous above
4.0, as the target space-time range of data to fit the model (Figure 2A
and Table 1 in Zhuang35). The events outside this target space-time
range are used as complementary events for calculating the conditional
occurrence rate in order to illuminate the boundary effect (see also
Zhuang35).
The same estimates as previously described35 were used in this study,i.e. =0.383 (events/day), =0.232 (events), =0.00578 day, =1.41(m−1), =1.08, =1.01×10−5 (deg2), =1.59 and =1.38 (m−1),
which is also similar to the results of Zhuang et al.9,10 The estimatedbackground seismicity rates are shown in Figure 3.
The spatial variations of the β-value (b-value) in the magnitude distribution are evaluated by using the MLE (maximum likelihood estimate) method: at the locations of each 0.1°×0.1° grid,where (x, y; k) represents the mean magnitude of the k closest earthquakes to (x, y), and 0.05 is the correction of the estimation error
caused by rounding the magnitude to 1 digit in the catalog. The results
are shown in Figure 4.
[page 56]

into equation
Substituting m(x, y), β(x, y), and the values of and(21), the probabilities that at least one event occurs in a 1°×1°×1-year
space-time window with a magnitude no less than MJ 5, MJ 6 and MJ 7 areevaluated and plotted in Figure 5. From Figure 5 a series of conclusions
can be drawn. i) The regions with highest probabilities are the Tohoku(northeastern Japan) Arc and the Kuril Arc, and the northern end of the
Ryukyu Arc, which is east to the Kyushu Island. ii) Comparing these locations to the Tokai and Nankai regions, which are both offshore south of
the Honshu and the Shikoku islands, such probabilities are relatively
lower than the Tohoku (northeastern Japan) Arc and the Kuril Arc
regions, and the northern end of the Ryukyu Arc. iii) Along the middle
ridge of Honshu island, which makes up the major part of the volcanic
front in the Japan region, occurrence probabilities in unit space-time volume are much lower than in the coastal and offshore regions. iv) The
coastal and offshore regions of the sea around Japan have similar earthquake probabilities to the Tokai and Nankai regions.

Conclusions
The ETAS model gives a natural description of earthquake clusters
that can be used as a starting point for studies of seismicity patterns.

[Research in Geophysics 2012; 2:e8]

Articleonlyground seismicity rate during the 1997–1998 Umbria-Marche, central Italy, sequence: apparent variation or fluid-driven triggering?
Bull Seismol Soc Am 2010;100:1138-52.
13. Zhuang J, Ogata Y, Vere-Jones D. Stochastic declustering of spacetime earthquake occurrences. J Am Stat Assoc 2002;97:369-80.
14. Kagan Y, Knopoff L. Statistical search for non-random features ofthe seismicity of strong earthquakes. Phys Earth Planet Inter
1976;12:291-318.
15. Omori F. On the aftershocks of earthquakes. J Coll Sci Imperial Univ
Tokyo 1894;7:111-200.
16. Utsu T. A method for determining the value of b in a formula logn =a − bM showing the magnitude-frequency relation for earthquakes.
Geophys Bull Hokkaido Univ 1965;13:99-103.
17. Utsu T, Ogata Y, Matsu’ura RS. The centenary of the Omori formulafor a decay law of aftershock activity. J Phys Earth 1995;43:1-33.
18. Zhuang J, Chang CP, Ogata Y, Chen YI. A study on the backgroundand clustering seismicity in the Taiwan region by using a point
process model. J Geophys Res 2005;110:B05S13.
19. Ogata Y, Zhuang J. Space–time ETAS models and an improved
extension. Tectonophysics 2006;413:13-23.
20. Ogata Y. Space-time point-process models for earthquake occurrences. Ann Inst Stat Math1998;50:379-402.
21. Console R, Murru M, Lombardi AM. Refining earthquake clusteringmodels. J Geophys Res 2003;108:2468.
22. Daley DD, Vere-Jones D. An introduction to theory of point processes – Vol 1: Elementary theory and methods (2nd ed). New York:
Springer; 2003.
23. Fletcher R, Powell MJD. A rapidly convergent descent method forminimization. Computer J 1963;6:163-8.
24. Musmeci F, Vere-Jones D. A variable-grid algorithm for smoothingclustered data. Biometrics 1986;42:483-94.
25. Helmstetter A, Kagan YY, Jackson DD. Comparison of short-term
and time-independent earthquake forecast models for southern
California. Bull Seismol Soc Am 2006;96:90-106.
26. Werner MJ, Helmstetter A, Jackson DD, Kagan YY. High-resolution
long-term and short-term earthquake forecasts for California. Bull
Seism Soc Am 2011;101:1630-48.
27. Ogata Y. Space-time model for regional seismicity and detection ofcrustal stress changes. J Geophys Res 2004;109:B03308.
28. Ogata Y. On Lewis’ simulation method for point processes. IEEE
Transact Informat Theory 1981;IT-27:23-31.
29. Zhuang J. Second-order residual analysis of spatiotemporal point
processes and applications in model evaluation. J Royal Stat Soc
Series B (Stat Methodol) 2006;68:635-53.
30. Picard R, Cook D. Cross-validation of regression models. J Am Stat
Assoc 1984;79:575-83.
31. Zhuang J, Ogata Y. Properties of the probability distribution associated with the largest event in an earthquake cluster and their implications to foreshocks. Phys Rev E 2006;73:046134.
32. Saichev A, Sornette D. Distribution of the largest aftershocks inbranching models of triggered seismicity: theory of the universalBåth law. Phys Rev E 2005;71:056127.
33. Zhuang J, Harte D, Hainzl S, et al. Basic models for seismicity.
Commun Online Resource Stat Seism Analysis 2012; ThemeV.
Available from: url:http://www.corssa.org/
34. Vere-Jones D, Zhuang J. Distribution of the largest event in the critical epidemic-type aftershock-sequence model. Phys Rev E 2008;
78:047102.
35. Zhuang J. Next-day earthquake forecasts by using the ETAS model.
Earth Planet Space 2011;63:207-16.

useBased on the stochastic declustering method, we can objectively decompose the seismicity into the background and clustering components.
The background component can be used to estimate long-term earthquake hazards as declustered catalogs (declustered by window-based or
link-based declustering methods). However, only using a simple combination of the background rate and the Gutenberg-Richter law for themagnitude distribution might underestimate such risks, since the maximum magnitude in an earthquake cluster no longer obeys the
Gutenberg-Richter law, but a different distribution as described aboveby equation (14). Instead, more viable forecasts can be made by using
the equation (21) to correct the bias. For data analysis, the seismicity inand around the Japan region was used to show how to apply the procedures and techniques introduced in this study. The results show that the
Tohoku Arc and Kuril Arc are the regions at the greatest potential riskof having large earthquakes in the future and that this risk is muchhigher than in the Tokai and Nankai regions where scientists and the
public at large expect one or more earthquakes of M7.5 or over.
The data analysis in this study is based on the hypothesis that future
seismicity follows similar patterns seen in the past. This may not be true.
For example, some areas currently in a quiescent period may be re-activated in the future. Certainly, any practicable models for long-term earthquake forecast should consider these possible large-scale changes in seismicity. However, before considering more advanced conditions, a solidreference model should be built to provide a null hypothesis by which totest whether these more complicated models offer any advantages.

alReferences

N

on

-c

om

m

er

ci

1. Utsu T. Aftershock and earthquake statistics (l): Some parameterswhich characterize an aftershock sequence and their interrelations.
J Fac Sci Hokkaido Univ Ser VII 1969;3:129-95.
2. Gardner J, Knopoff L. Is the sequence of earthquakes in southern
California, with aftershock removed, Poissonian? Bull Seismol Soc
Am 1974;64:1363-7.
3. Keilis-Borok VI, Kossobokov VI. Time of increased probability for thegreat earthquakes of the world. Comput Seismol 1986;19:45-58.
4. Reasenberg P. Second-order moment of central California seismicity, 1969-1982. J Geophys Res 1985;90:5479-95.
5. Frohlich C, Davis SD. Single-link cluster analysis as a method toevaluate spatial and temporal properties of earthquake catalogs.
Geophys J Int 1990;100:19-32.
6. Davis SD, Frohlich C. Single-link cluster analysis, synthetic earthquake catalogs, and aftershock identification. Geophys J Int
1991;104:289-306.
7. Console R, Murru M. A simple and testable model for earthquake
clustering. J Geophys Res 2001;106:8699-711.
8. Helmstetter A, Ouillon G, Sornette D. Are aftershocks of large
Californian earthquakes diffusing? J Geophys Res 2003;108:2483.
9. Zhuang J, Christophosen A, Savage MK, et al. Differences between
spontaneous and triggered earthquakes: their influences on foreshock probabilities. J Geophys Res 2008;113:B11302.
10. Zhuang J, Ogata Y, Vere-Jones D. Analyzing earthquake clusteringfeatures by using stochastic reconstruction. J Geophys Res 2004;
B5:B05301.
11. Hainzl S, Ogata Y. Detecting fluid signals in seismicity data throughstatistical earthquake modeling. J Geophys Res 2005;110:B05S07.
12. Lombardi AM, Cocco M, Marzocchi W. On the increase of back-[Research in Geophysics 2012; 2:e8][page 57]