Some Applications of Point
Processes in Seismicity Modellingand Prediction

Zhuang Jiancang

Doctor of Philosophy

Department of Statistical Science
School of Mathematical and Physical Science
The Graduate University for Advanced Studies

2003

ii

Some Applications of Point
Processes in Seismicity Modellingand Prediction

Zhuang Jiancang

A Thesissubmitted to Department of Statistical Science
the Graduate University for Advanced Studies
in ful llment of the requirements
for the degree of
the Philosophical Doctor in Statistics

March 2003

ii

Abstract

This thesis gives some results on the applications of point process in seismicity modellingand prediction, which can be divided into 3 parts.
The rst part discusses the general strategies of using point processes to model andpredict seismicity. The conditional intensity function is used for the model speci cation,
which has a natural form for prediction. Some models are discussed in this thesis including
the Poisson model, the stress release model, the self- and mutual-exciting process and theepidemic type aftershock sequence (ETAS) model. The parameter estimation, simulation,
prediction and scoring procedures are also outlined in this part.
The second part deals with the space-time ETAS models. I will rst analyze its mathematical properties, including criticality, moments and its simulation algorithm. In thesections on criticality, I will discuss not only the critical parameter but also the meaningof the eigenfunctions of the integral equation. In the simulation section, I will give analgorithm which can simulate a catalogue of events by using the space-time ETAS models
much faster than by the original thinning procedure given by Ogata (1998). Anotherfocus is placed on a technique called stochastic declustering. Stochastic declustering isconcerned with an objective estimation of the spatial intensity function of the backgroundearthquake occurrences from an earthquake catalogue. This algorithm can also be used
to produce declusterred catalogues, or to separate the earthquakes into explicit clusters
in an objective fashion. The background intensity function can be evaluated if the total
spatial intensity of seismicity and the branching structure are estimated. Or, the wholeiii
space-time process can be stochastically split into two sub-processes, the background andthe clustered ones. This algorithm combines a parametric estimation of the clustering
structures in the space-time ETAS model and a non-parametric intensity estimation of
the background seismicity that we call weighted variable kernel estimate. As a demonstration of the present methods, we estimate the background seismic rates in the central
New Zealand region, in the central and western Japan region and in the North China
region to produce the declustered catalogues.
The third part deals with the prediction of earthquakes based on precursory observations. I will give several models that may be useful in modelling the possible relations
between precursors and earthquake occurrences. An example is given for modelling therelationship between electric signal precursors and earthquakes around the Beijing region.
The precursory data are the underground electric signals from 4 stations monitoring theultra-low frequency components of electric eld in the vicinity of Beijing, and they areused as explanatory variables in predicting the occurrence of events with M  4 within a
300 km circle centered at Beijing. A modi ed version of Ogata's LIN-LIN model is used
to examine the in uence of a sequence of explanatory signals to the occurrence of eventsin a stochastic point process. The explanatory e ect is shown to be highly signi cant, andsubstantially superior to the explanatory e ect of the same signals applied to a randomized version of the earthquake data. All four stations show signi cant explanatory power,
although the combination of the two most e ective tends to dominate the forecasts. The
results are stable against perturbations in the time period or region of observation. The
prediction performance appear to be most e ective for the events with M  5, and forevents closer to the observing stations, although some of the smaller events appear to
produce detectable signals at distances of over 100 km from the source. Probability gains
over the simple Poisson process are in the range of 3 { 4 for the events of magnitudes 5or larger. A special study is made of predicted and unpredicted events in the subregion

ABSTRACT

iii

around the Ms 7.8 Tangshan earthquake of 1976, to reveal the common particular spatial
pattern of the classi ed events corresponding to all individual stations.

iv

Contents

Abstracti

List of Figuresx

List of Tables

xiv

1 Introduction and Motivation

1

I GENERAL IDEAS

5

2 Basic Concepts in Point Processes

7

3 Modelling strategies

13

3.1 Model formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1.1 Poisson models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Self-correcting models: stress release model and linked stress releasemodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.3 Self and mutually exciting models. . . . . . . . . . . . . . . . . . . 16v

vi

CONTENTS

3.2 Parameter estimation and model selection . . . . . . . . . . . . . . . . . . 17
3.3 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Prediction and scoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4.1 Prediction procedures . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4.2 Scoring procedures and information gain . . . . . . . . . . . . . . . 22

II SPACE-TIME ETAS MODELS AND THEIR APPLICATIONS
25
4 The Epidemic Type Aftershock Sequence Model27

4.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 De nition of space-time ETAS models . . . . . . . . . . . . . . . . . . . . 30
4.3 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.1 Criticality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.2 Moment properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.3 Spatial densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4 Simulation algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4.1 Thinning method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4.2 Composition method . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.5 Prediction performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

CONTENTS

vii

5 Stochastic Declustering

55

5.1 Why stochastic declustering? . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2 Thinning procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.3 Variable kernel estimates of seismicity . . . . . . . . . . . . . . . . . . . . . 58
5.4 Iteration algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.5.1 Central New Zealand region . . . . . . . . . . . . . . . . . . . . . . 62
5.5.2 Central and western Japan . . . . . . . . . . . . . . . . . . . . . . . 65
5.6 Probability-based declustering algorithm . . . . . . . . . . . . . . . . . . . 66
5.7 Results for the North China region . . . . . . . . . . . . . . . . . . . . . . 72
5.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.9 Conclusion of Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

III MODELLING EARTHQUAKES AND PRECURSORS 87
6 Formulation of Models on Earthquake and Precursors

89

6.1 On precursory observation data . . . . . . . . . . . . . . . . . . . . . . . . 89
6.2 The LINLIN model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3 Extensions of the LINLIN model . . . . . . . . . . . . . . . . . . . . . . . 92
6.4 On precursory observations at discrete times . . . . . . . . . . . . . . . . . 92

viii

CONTENTS

6.5 On the models for more multi-type precursors and multi-station observations 95

7 Statistical Analysis on the ULF Electric Signals

97

7.1 Data and model speci cation . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.2 Results for regional seismicity within the circles around the individual stations104
7.3 Results for the seismicity within 300 km circle around Beijing . . . . . . . 109
7.4 Results for the Tangshan region . . . . . . . . . . . . . . . . . . . . . . . . 116
7.5 Probability gains of forecasts based on electrical signals . . . . . . . . . . . 118
7.6 E ect of the change in the background seismicity . . . . . . . . . . . . . . 122
7.7 Testing whether the electric signals are postseismic phenomena . . . . . . . 124
7.8 Summaries of Chapter 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

IV APPENDIX131

A Observations on the ULF Ground Electric Field in China

133

A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
A.2 Monitoring network and observation system . . . . . . . . . . . . . . . . . 134
A.2.1 Network distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 134
A.2.2 Observation system . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
A.2.3 Methods for noise reduction . . . . . . . . . . . . . . . . . . . . . . 137

CONTENTS

ix

A.3 Possible noise sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
A.3.1 Recognition of interference signals . . . . . . . . . . . . . . . . . . . 141
A.3.2 Analysis of the e ects of air temperature, rainfall and lightning . . 143
A.3.3 Magnetic storms and electric signals . . . . . . . . . . . . . . . . . 144
A.4 Anomaly recognition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
A.4.1 About the data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
A.4.2 Classi cation of signals . . . . . . . . . . . . . . . . . . . . . . . . . 149
A.5 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
A.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

B On the Optimization Algorithm

153

Acknowledgement

159

Bibliography

161

The Author's Publications and Submissions During 2000 to 2003

171

x

CONTENTS

List of Figures

4.1 Illustration of the simulated catalogue. . . . . . . . . . . . . . . . . . . . . 51
4.2 Prediction probabilities for di erent magnitude classes . . . . . . . . . . . 52
5.1 Seismicity in the central New Zealand region. . . . . . . . . . . . . . . . . 62
5.2 Results from the central New Zealand region. . . . . . . . . . . . . . . . . 64
5.3 Seismicity the western and central Honshu area, Japan.

. . . . . . . . . . 66

5.4 Results from the central and western Honshu area of Japan. . . . . . . . . 67
5.5 Histograms of aftershock probabilities (c.f (5.2)) for (a) the central New
Zealand region; and (b) the western and central Honshu region. . . . . . . 68
5.6 Results from applying Algorithm 2 to the central New Zealand region. . . 70
5.7 Results from applying Algorithm 2 to the western and central Honshu area. 71
5.8 Seismicity and subregions in the north China region (1970{1999, M  3:5). 75
5.9 Estimated spatial intensity of the total process, the background process
and the clustering process.

. . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.10 (continued of Figure 5.9) Estimated spatial intensity of the total process,the background process and the clustering process. . . . . . . . . . . . . . 77
xi

xii

LIST OF FIGURES

5.11 Space-time distributions of the background events that are stochasticallydeclustered from the North China data set. . . . . . . . . . . . . . . . . . 79
5.12 Space-time distributions of the triggered events that are complements ofstochastically declustered events from the North China data set.

. . . . . 80

7.1 Epicenters of earthquakes occurred in the Beijing region and its vicinity
during the period 1982.1.1 { 1998.1.31. . . . . . . . . . . . . . . . . . . . . 99
7.2 Signals observed at each station and the occurrence times of the earthquakes within 300 km from the stations. . . . . . . . . . . . . . . . . . . . 100
7.3 Plots of cumulative numbers (top panel) and magnitudes (second panel) of
the earthquakes that occurred within 300 km from Beijing against occurrence times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.4 Estimated response functions Am;` by tting the model to the 300 km circles
around stations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
7.5 Plot of the conditional intensity function tted to the data from Sanhe (S)station. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
7.6 Plot of the tted conditional intensity function of the multi-predictor modelwith inputs from Qingxian (Q) and Sanhe (S) stations. . . . . . . . . . . . 113
7.7 Epicentral locations of predicted and unpredicted events. . . . . . . . . . . 114
7.8 Plots of magnitudes of the predicted earthquakes and unpredicted earthquakes against their distances to the stations. . . . . . . . . . . . . . . . . 115
7.9 Correspondence between the same conditional intensity as in Figure 7.6
and events with random time points. . . . . . . . . . . . . . . . . . . . . . 117

LIST OF FIGURES

xiii

7.10 Epicenter distribution of the predicted and unpredicted earthquakes in theTangshan region (the magni ed rectangular area in Figure 7.1). . . . . . . 119
7.11 Histogram of mean probability gains. The mean probability gain from theactual model is mark by a crossed circle. . . . . . . . . . . . . . . . . . . . 122
7.12 Results for Huailai Station using the electric signal events as the target
point process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
7.13 Results for Langfang using the electric signal events as the target point
process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
7.14 Results for Qingxian using the electric signal events as the target point
process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7.15 Results for Sanhe using the electric signal events as the target point process.129
A.1 Distribution of observation stations. The lines denote the locations of major faults . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
A.2 Illustration of the observation system (a) and the response function (b). . . 138
A.3 A photocopy of a record with noises (the upper panel) from ground strong
lightnings and (the bottom panel) from electricity leaked by power lines. . 142
A.4 3-day average of the electric signals (solid curve), the lightnings (bar plot)
and the earthquakes (circles) with MS > 5:2 in the monitoring region. . . . 145
A.5 3-day average of the electric signals (solid curve), the magnetic storms (barplot) and the earthquakes (circles) with MS  5:2 in the monitoring region. 147
A.6 Normal record (the top panel) and classi cation of anomalies: Class A (thesecond panel), Class B ( the third panel) and Class C (the bottom panel). 151

xiv

LIST OF FIGURESList of Tables

4.1 Prediction scores for a simulated data by using ETAS model . . . . . . . . 54
5.1 Convergence steps of the evaluation procedure in applying the space-time
ETAS model by Algorithm 1 to the earthquake data from the central New
Zealand region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.1 Summary of results for data within 200 km and 300 km circles aroundstations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
7.2 Performance comparisons for single and multiple predictor models in the
300 km region around Beijing . . . . . . . . . . . . . . . . . . . . . . . . . 110
7.3 Numbers and proportions of earthquakes, and days above and below themean occurrence rate of earthquake . . . . . . . . . . . . . . . . . . . . . . 112
7.4 Probability gains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.5 Comparisons of AIC values for the 1-stage and 2-stage background models
and the 1-stage background models . . . . . . . . . . . . . . . . . . . . . . 123
7.6 Fitted results by using the electric signals as the main process and earthquakes within 300km to the station as the external process . . . . . . . . . 124
A.1 Number of lightnings occurring in the monitoring region during 1987{1992 144
xv

xvi

LIST OF TABLES

A.2 Number of magnetic storms occurring during 1987{1982 recorded at Baijiatan Station, Beijing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
A.3 Correlation between magnetic storms and electric signals . . . . . . . . . . 146

Chapter 1

Introduction and Motivation

Point processes is the core subject in statistical seismology. The term \statistical seismology", rst appeared in Aki (1956), and then was adopted by Vere-Jones (1995), etc. The
main task in earthquake prediction, as pointed out by Aki (1984), is to develop statisticalmodels for analyzing the observations, based on which we can evaluate the probability ofearthquakes occurring in a certain space-time-magnitude window.
The idea of probability prediction was not easy to be accepted by seismologists andgeophysicists, who had been engaged in research on the deterministic earthquake prediction for many decades. However, they have gradually realized that the occurrences ofprecursors and earthquakes are not likely a simple one-to-one relationship. Some kinds ofprecursory phenomena may be observed before a large earthquake, but other earthquakes
occurred without these precursors, or even no earthquake occurred after we observed theseanomalies.
The central topic of this thesis is on the development of point process models to predict
the earthquake risk based on some observations in previous time, including earthquake
occurrences and precursors, i.e., to nd a suitable function

f (t1 ; t2 ) = risk of earthquakes in (t1 ; t2 ) given observation history
1

(1.1)

2

CHAPTER 1.

INTRODUCTION AND MOTIVATION

or in the derivative formf (t)dt = risk of equakes in (t; t + dt) given observation history(1.2)

to predict the earthquake occurrence in a certain time interval. Naturally, the history ofobservation is treated into two parts, records of earthquakes and records of precursors.
Modelling strategies based on stochastic point processes usually start from the completely random model or the Poisson process, and we then make some e orts to modify
this model to a new model with less randomness. Once a new model with less randomness
is proposed, some statistical principles are used to make a judgement whether the newmodel ts the observation better or not. If so, the new model becomes a new starting
point for developing further models.
Chapter 2 of this thesis gives some concepts in the theory of point processes. Since
this thesis will mainly deals with application proposes, only de nitions and properties ofnecessary concepts are summarized.
Chapters 3 discusses the strategies in modelling and predicting seismicity by using
point processes, including procedures of model formulation, parameterization, model selection, simulation, probability forecasting and scoring. In the model formulation part,several models, such as Poisson models, stress release models and self- and externalexciting models are also brie y discussed.
Of particular interests, the epidemic type aftershock sequence (ETAS) model is discussed in Chapter 4 in detail. The materials on the ETAS include its history, a generalde nition, its criticality, moment measure evaluations.
Chapter 5 focuses on the application of the ETAS model to the important problem ofcatalogue declustering. Clusters in earthquakes always complicate the analysis. Due tolack of clear knowledge of the underground seismicity, various deterministic declustering

3
methods have been developed by many researchers in subjective ways. The proposed
stochastic declustering method here allows the existence of uncertainties in the classi cation. Each earthquake cannot de nitely be said as a mainshock or an aftershock, but canbe said as a mainshock or an aftershock with a certain probability. The declustering canthen be realized by a re-sampling from the catalogue based on these probabilities.
Chapter 6 and Chapter 7 discuss the models for prediction based on some precursor
data, which is much useful but nevertheless less developed than models only based on thehistory of earthquakes themselves. In Chapter 7, the records of the ultra-low frequency
components of the underground electric signals are validated as the precursors of theearthquakes occurrences.

4

CHAPTER 1.

INTRODUCTION AND MOTIVATIONPart I

GENERAL IDEAS

5

Chapter 2

Basic Concepts in Point Processes

There are two traditional approaches in describing a point processes on the time axis or
the real line. One is based on counting of points on intervals. In such a way, the number ofpoints in an interval on the real line is a non-negative integer valued random variable. To
specify a point process is to give consistently de ned nite-dimensional joint distributionsof such random variables for any mutually disjoint intervals. The other approach is bydescribing a sequence of interval lengths between consecutive points, and then treating
the sequence as a positive valued stochastic process. To extend these approaches to higher
dimensional Euclidean spaces, the second approach have diÆculties. But when the process
has a time-like evolutionary character, the second approach leads us to the importantconcept, conditional intensity function, which is the key tool in modelling seismicity.
In this chapter, I will give the de nition of some useful concepts associated with point
processes, including the de nition of a point process on an Euclidean space, the momentmeasures and the conditional intensity.

Point measures and point processes Generally speaking, a point process is a randommeasure de ned on the considered space such as the real line. The strict de nition of thepoint process can be found, e.g., in Daley and Vere-Jones (1976), Karr (1991) or Reiss(1993).
7

8

CHAPTER 2.

BASIC CONCEPTS IN POINT PROCESSESLet E be an Euclidean space, with Borel  -algebra E . Let B denote the ring of bounded
Borel sets. A point measure  on E is a measure that (A) 2 N and (A) < 1 for every

A 2 B.
The fundamental point measures are the point masses Æx , de ned byÆx (A) = 1(x 2 A) =

1; if x 2 A; ;0; if x 2= A.

(2.1)

where the notation Æx is commonly called the Dirac measure.
A point measure  can be decomposed as

=

K
X
i=1

ai Æxi ;(2.2)

where 0  K  1, ai are positive integers and the xi are distinct points in E . If ai = 1
for each i, the point measure is called a simple point measure.
Let ( ; F ; P ) be a probability space. A point process on E is a measurable mapping

N of ( ; F ) into (Mp ; Mp) where Mp is the set of all the point measures on E and Mp
is the  -algebra. For a random measure M and A 2 E , the random variable M (A) is
termed as the number of points falling into A, and, for a function in CK , the integralR

N (f ) = E f (x)N (dx).

Intensity measures and moment measures The intensity measure m, also called
mean measure or rst moment measure, of a point process N is de ned by the expectationsm(B ) = E[N (B )]; for B 2 B:

(2.3)

Hence, m(B ) is the expected number of points in the set B .
More generally, for any positive integer k, let N k be the point process given by

N k (dx1 ;    ; dxk ) = N (dx1 )    N (dx2 );

(2.4)

9
then the moment measure of order k of the point process N is the measure mkN = E[N k ]or mk if no other processes than N are considered.
The covariance measure of N is the signed measurecN (A  B ) = Cov(N (A); N (B )) = m2N (A  B ) mN (A)mN (B )(2.5)

For many purposes, ordinary moment measure of point processes are unsatisfactory
because the products N k contain points on various diagonals. Without contributing
information not already present, such points engender spurious singularity of momentmeasures. Their removal is e ected by replacing N k by the point process!
k 1
X
N (k) (dx1 ;    ; dxk ) = N (dx1 )(N Æx )(dx2 )    N
Æxi (dxk );
i=1
1

(2.6)

which provided that N is simple, consists only k-tuples of points whose coordinates are
distinct. We refer to the measure m(k) = E[N (k) ] as factorial moment measures of N , orNm(k) if no other processes than N are considered..

Conditional intensities When the process has a time-like evolutionary character. The
most natural way to describe such a process is by successive conditioning: at each stage,
we condition on the times and other characteristics of the points which have gone before.
The intensity speci cation of point processes is developed closely relating to the information theory. Consider a prediction of an event occurring on a small time interval. That
is to say, suppose a point process on the half-line (0; 1), and divide it into small intervalsof length . Then, we get a stochastic process fk g, where k = N [(k1); k) is k-th

random variable on the subinterval [(k 1); k). For a small enough  we may assume
that fk g is a binary process (i.e., 0 or 1 valued). If the considered point process is a stationary Poisson, then fk g is identically and independently distributed, that is, a Bernoulli
series. Otherwise, it is easily seen that, in general, the joint probability of the sequence is

10

CHAPTER 2.

BASIC CONCEPTS IN POINT PROCESSES

essentially determined by a sequence of conditional probabilities P fk = 1j1 ; : : : ; k 1 g,

k = 1; 2; : : : ; on the history of events. Roughly speaking, as a derivative of the conditional
probability with respect to the time, the conditional intensity function (tjFt ) is de ned
by

P fN (t; t + ) = 1jFt g = (tjFt ) + o();in other words,

(tjFt ) = lim
Pfan event in [t; t + ) j Ft g=;
!0

(2.7)

where Ft is a set of observations up to time t, equivalent to the history of the occurrence
times of the events Ht = fti ; ti < tg (see Rubin, 1972, for example). With this reason,
we also usually write the conditional intensity as (tjHt ). The conditional intensity completely characterizes the corresponding point process (see Liptzer and Shiryaev, 1978).
From the de nition of the conditional intensity function, we can get
PrfAt least 1 event occurs in [t; t + s) j Ht g = 1 PrfNo event occurring in [t; t + s) j Ht g
Z t+s= 1 exp
(ujHt )dut= 1 S (t + s j Ht )(2.8)

The above conditional intensity is related with the hazard function. Suppose thatt1 < t2 <    < tn <    and an events occurs at each i; i = 1;    ; n;    , we express thehazard functions in terms ofhn (tjt1 ;    ; tn 1 ) = pn (tjt1 ;    ; tn 1 )=Sn(tjt1 ;    ; tn 1 );(2.9)

where pn (tjt1 ;    ; tn 1 )=Sn(tjt1 ;    ; tn 1 ) is the density function of an event occurring at

t given the history of previous occurrences, and the survival function Sn (tjt1 ;    ; tn 1 ) =1Rt
tn pn (s t1 ;j    ; tn 1 )ds. Thus we havepn (tjt1 ;    ; tn 1 ) = hn (tjt1 ;    ; tn 1 ) exp



Z t
tn 1



hn (ujt1 ;    ; tn 1 ) du :(2.10)

11
Given a sequence fti g with 0 < t1 < t2 <    < tn <    , we de ne an amalgam of thehazard function by

(tjHt ) =

h1 (t);
0 < t  t1 ;hn (tjt1 ;    ; tn 1 ); tn 1 < t  tn ; n  2:

(2.11)

Thus conditional intensity function for a regular point process on [0; 1) can be alsode ned piecewise by (2.11).
The above de nition can also be extended to the whole time axis. Suppose that   < t n <    < t 2 < t 1 < t0 < t1 < t2 <    < tn <    , we write the hazard functionaspn(tjt k ;    ; tn 1 );hn (tjHt ) = klim
!1 Sn (tjt k ;    ; tn 1 )(2.12)

provided that the limit exists. The conditional intensity function on a in nite history (seeDaley and Vere-Jones, 1984, Chapter 13) can be sequentially de ned by

(tjHt ) = hn (tjHt );tn 1 < t  tn :(2.13)

12

CHAPTER 2.

BASIC CONCEPTS IN POINT PROCESSES

Chapter 3

Modelling strategies

3.1 Model formulation
Formulation is the rst stage of modelling, for which there is no already-done or cure-all
methods. It depends on the experience, insight and imagination of the researcher. A goodmodel should be sensible enough to embody physical ideas about the processes supposed
to be represented, or to be able to t the processes suÆcient well. To satisfy such a
condition, it is helpful to use computer-aided methods to display the features of the data,in order to get an overview of the studied data.
In practice, there are quite a few models developed for modelling seismicity, which canbe roughly classi ed into two types of models: the inhibitory type and the exciting type.
The inhibitory type models include renewal processes and stress release models, while theexciting type models include epidemic type aftershock sequence models. In this chapter,
the Poisson model and the stress release model is explained in more detail. The ETAS
model will be discussed in next chapter.
13

14

CHAPTER 3. MODELLING STRATEGIES

3.1.1 Poisson models
This Poisson model is an example of the simplest point process, which is often used as
the null model in hypothesis tests.
Let  is a non-negative valued measure. A point process N on E is a Poisson processwith the mean measure  if(a) Whenever the sets A1 ; : : : ; Ak 2 E are disjoint, the random variables N (A1 ); : : : ; N (Ak )are independent;(b) For each A 2 E and K  0,
PrfN (A) = kg =(A)kexp[ (A)]:
k!

(3.1)

The distribution of the time interval  from the present time, say t, to the occurrenceof the next event is determined by the equation
Prf > xg = PrfN [(t; t +  )] = 0g
 Z t+x

= exp
(dt)
t

= exp[ ((t + x) (t))];

(3.2)

where
(t) =Z t

0

(ds):(3.3)

In the case of a stationary Poisson process, (t) = t, indicating that probability in(3.2) is just the exponential distribution. The above Poisson model in time can be easily
extended to space-time Poisson models.

3.1.

15

MODEL FORMULATION

3.1.2 Self-correcting models: stress release model and linked
stress release models
The stress release models and, more generally, self correcting models were studied by VereJones (1978), Ogata and Vere-Jones (1984), Vere-Jones and Ogata (1984), and Isham andWstcott (1971). The idea of this model is due to the elastic rebound theory. Typical selfcorrecting models have the conditional intensity

(tjHt ) = a + bt cN [0; t);(3.4)

where a, b and c are constants

The stress release model. The stress release model origins from the idea of the elasticrebound theory, and it has the conditional intensity function

(tjHt ) = exp[a + bt cS (t)];where

S (t) =

X
i:ti <t101:5Mi ;

(3.5)

(3.6)

representing the cumulative stress released by the events up to time t; a, b and c areconstants. Zheng and Vere-Jones (1992, 1994) applied this model the several earthquake
data sets including the historical earthquake data from North China, Japan and Iran.

The linked stress release model. The linked stress release models consider the interaction of stresses among di erent regions. Suppose that we consider K subregions. In eachregion, the earthquake risk is a nonnegative monotonous increasing function of its stresslevel, which increases with time. Once an earthquake occurs in one of the subregions,say, subregion A, it not only reduces the stress level in A but also may increase or
decrease the stress levels in the other subregions. If we denote the stress levels by a

16

CHAPTER 3. MODELLING STRATEGIES

vector X= (X1 ;    ; XK )T , where the ith component represents the stress level in the ithsubregion, then=a + t C s (t)X  bor in full form
0

1

0

1

0

1

X1 (t)a1b1
B X (t) C B a C B b C
B 2
C B 2 C B 2 C
B
C = B . C+B . Ct.
.
@
. A @ .. A @ .. A
XK (t)
aK
bK

(3.7)

0

1

c11 c12    c1K 0 s1 (t) 1
B c
CC
B 21 c22    c2K C B
B s2 (t) C ;
B ..
C
.
.
.
..
. . .. A @    A@ .
s4 (t)
cK 1 cK 2    cKK

where ai , bi and cij are parameters, and si (t) =

(3.8)

P
j : tj <t Si;j is the stress released by theearthquakes in the ith subregion up to time t. Thus the conditional intensity function is
0

1

0

1

1 (tjHt )exp[X1 (t)]B  (tjH ) C B exp[X (t)] C
t C B
2
B 2
C
=
B
C
B
C
.
.
..
..
@
A @
A
K (tjHt )exp[XK (t)](3.9)

The physical meaning of this model is that each subregion accumulates its own stressindependently, but once an earthquake occurs in a subregion, it not only releases thestress in this subregion, but also changes the stress level of other subregions, increasingor decreasing, depending on the sign of the non-diagonal elements cij .
The application of this model can be found in Liu, Chen et al. (1999), Bebbingtonand Harte (2001) and Lu, Harte and Bebbington (1999).

3.1.3 Self and mutually exciting models.
The original self- and mutual-exciting models were proposed by Hawkes (1971). This typeof point processes has three components as follows. The background process is controlled
by an intensity function of time representing the trend. Once an event occurs, it producesanother Poisson process which will decay o in a certain period. Beside these, the process
is also in uenced by an external point process, i.e., any events occurring in the external

3.2.

PARAMETER ESTIMATION AND MODEL SELECTION

17

point process also excites events with similar decay o in a certain period. The conditional
intensity of this process fti g can be written as

(tjHt ) = f (t) +X
i:ti <tg (t ti ) +X
k:uk <th(t uk );(3.10)

where f (t) is the trend function, g (t) is the impulse response function from an event inthis process, h(t) is the response function from an event in the external process fuig, and

Ht is the observation history of both the processes ftig and fuig.

3.2 Parameter estimation and model selection
The parameter estimation can be carry out be the maximum likelihood method. To
maximize the likelihood function, some standard nonlinear optimization subroutine canbe used, such as the Daviden-Fletcher-Porwell method. Incidentally, see Appendix B foran example of gradient free local optimization methods.
When plural models are tted to the same observation data, it is important to examine
which model best ts to the data. This can be implemented by using the AIC (Akaike
Information Criterion; Akaike, 1974)

AIC = 2 log L^ + 2K;(3.11)

where L^ represents the maximized likelihood and K is the number of parameters used in
the model. The model that has the smallest AIC value is taken as the best tted.

3.3 Simulation
There are many reasons why simulation is the important tool. First of all, it is a visible
realization of a particular model. For examples, simulation of a given model with di erent
parameters can help us to nd and understand di erent characteristics of the model.

18

CHAPTER 3. MODELLING STRATEGIES

Secondly, the simulation helps us with improving the model. It is more important toknow what cannot be described by the model than what have been described. After wet the model to the observational data, we can simulate synthetic data with the estimatedparameters. Then comparison of the simulation output with the original data helps us
with nding out the reasons of the discrepancies, which reveals some shortcomings ofthe model. These ndings will be taken into consideration of further step in the model
construction. Finally and most importantly, it links between simulation and prediction.
A complicated model has often no explicit prediction formulae, but we can still predict
the future of the process by simulating the same process for many times.

Simulation for Space-Time-Magnitude Poisson models
For a proper model for seismicity, it should include the components of the occurrencetime, the epicenter location and the magnitude. The complexities in the simulation of
space-time-magnitude Poisson models can be caused by several aspects.
1. Dependencies between components. If the above three components of the model are
independent of each other, we simulate each of them independently. Otherwise, it
may be more complicated.
2. Stationarity and homogeneity. A stationary Poisson model is always easier thana non-stationary one to be simulated. In the former case, the occurrence times of
events are uniformly distributed on the time axis. Similarly, when the process hasa homogeneous intensity, the locations of the events are uniformly distributed inspace.
3. Shape of occurrence regions.
Let us start from the most diÆcult case of a Poisson model with a intensity function

(t; x ; M ) on a region S . It can be simulated by sequential conditioning by generating

3.3.

19

SIMULATION

each component one by one, but we can use the rejection method to simulate for the mostgeneral cases.

Algorithm 3.1 (Rejection method)
1. Choose a density function  (t; x ; M ) that is easier to be simulated, such that C (t; x

; M ) > (t; x ; M ) on a regular region S 0  S , where C is a constant.

2. LetR

R

R= C M dM S 0 dx 0T  (t; x ; M )dt, where M represents the range of magni-

tude.
3. Generate a poisson random variable K with the mean .
4. Generate K random vectors (ti ; x i ; Mi ) with the density function  (t; x ; M ).
5. Generate K random numbers, U1 ; U2 ;    ; UK , uniformly distribution in [0; 1].
6. For each i, accept (ti ; x i ; Mi ) if CUi < (ti ; x i ; Mi )= (ti ; x i ; Mi ) and x i 2 S hold;otherwise, reject (ti ; x i ; Mi ).
7. Sort up all the accepted (ti ; x i ; Mi ) by ti with the ascending order.
This method is also sometimes called the thinning method, to distinguish from theoriginal rejection method for simulating random variables (for the references of the rejection method, see, e.g., von Neumann, 1951; Liu, 2001; and see Lewis and Shedler, 1979;
Ogata, 1991, for the thinning method), because this algorithm is for generating point
processes and must take into account the number of points to be generated.
The choice of the envelop function  (t; x ; M ) is the essentially technique in using therejection method. The simplest choice is to set  (t; x ; M ) a constant function. A more
complex choice is a nonconstant, separable function of the time, the location and themagnitude to be independent of each other, i.e., to set  (t; x ; M ) = t (t)x (x )M (M ), sothat we can generate each coordinate independently.

20

CHAPTER 3. MODELLING STRATEGIESThe choice of the covering region S 0 usually is among the ones such that the integral
in Step 2 can be easily obtained.
If (t; x ; M ) has a simple form, the above algorithm can be also simpli ed. For example, if (t; x ; M ) = uj (M ), where j (M ) is normalized as a probability function representing the magnitude distribution, and the simulation is on a regular region S , Algorithm 1
can be simpli ed as Algorithm 2 by choose the envelop function  (t; x ; M ) = (t; x ; M ) =uj (M ), and the covering region S 0 = S .

Algorithm 3.2
1. Generate a Poisson random variable K with the mean uT jS j.
2. Generate K random vectors (ti ; x i ; Mi ), where ti is uniformly on [0; T ] and x i uniformly on S , and Mi is generated from the density function j (M ).
3. Sort up all the accepted (ti ; x i ; Mi ) by ti with ascending order.
The omittance of Steps 5 and 6 in Algorithm 1 is self-evident.

3.4 Prediction and scoring
In this section, we consider the problems of assessing earthquake predictions and scoring
their performance based on point process models. The most common form of such a
prediction is to predict whether an event of the process will fall in a predetermined timeinterval, or more complicated, also in a certain space window and magnitude window.
Unlike \yes/no" prediction, the prediction based on point process give the probability
forecasts. The probabilities are associated with the occurrences of events in the given
time-space-magnitude windows. The former type of prediction is still important and can

3.4.

PREDICTION AND SCORING

21

be regarded as being close to a decision problem, i.e., the decision whether or not to declarethe occurrence of a damaging earthquake. This problem becomes more important if theconsequences of making such a yes/no prediction as serious as the economic consequencesof the occurrence of a damaging earthquake.
The most rational basis to make such a decision is the probabilities of the possible
outcomes based on the point process prediction, if available. If so, such decision can bemade within the knowledge of cost-bene t analysis, with the forecast probabilities used
to assess the expected cost.

3.4.1 Prediction procedures
Suppose that we are given observations of earthquakes and precursors up to a time t,
and try to predict the occurrence of at least one event in a given time interval [t; t + )and a space-time window S using a formulated model. The prediction procedures

were suggested by Vere-Jones (1999) based on simulation of point processes, which canbe outlined as following:
1. Fit a model to the observation data up to the given time t to estimate the parameters
by using, e.g., maximum likelihood method.
2. Simulate a time  to the next event by using the (tjHt ), or simulate the occurrencetime to the next event and its location and magnitude by using the (t; x ; M jHt )when the conditional intensity has a space-time-magnitude form. The simulation
method generally depends on the form of the conditional intensity function.
3. Add the new events to the history, and return to step (2), replacing t by t +  .
4. Continue until the end of the time interval for forecasting, (t; t + ), has exceed.

22

CHAPTER 3. MODELLING STRATEGIES

5. Repeat the above simulation in Steps (2) to (4) for a large number of times, the proportion containing at least one event in the target prediction space-time-magnitude
window can be determined and used as estimates of the required forecast probabilities.
The above prediction procedure has several characters. Firstly, it requires a modelwith a fully speci ed conditional intensity function. Secondly, when precursor data arealso used for the prediction, this procedure has some diÆculties to predict the risk at
a time quite long after the current time. In this case, we have to predict not only theearthquake process, but also the precursor processes. Thirdly, our present knowledge ofthe earthquake preparation precess makes it very diÆcult to develop a more practicalmodel at the present time.

3.4.2 Scoring procedures and information gain
Assume that pk are the forecasts for a sequence of intervals (k; (k + 1)), where k =0; 1; 2;    , produced by the preceding prediction method. These forecasts can be tested
when the actual results Xk are recorded, in such a way that

Xk =



1; at least one event occurring in (k; (k + 1));0;
otherwise.

(3.12)The natural score for the prediction performance is the likelihood for the sequence of thetrial, namely the binomial score

B=

X
k

[Xk log pk + (1Xk ) log(1 pk )]:(3.13)For comparison, the above procedures are also applied to the observation data fromthe Poisson model with the mean rate  = E[(tjHt )]. Let p = 1

exp(  ) be theprobability estimate for each interval from this Poisson process. The ratio pk =p is called

3.4.

23

PREDICTION AND SCORING

the probability gain for the kth interval. The average di erence between the binomial
scores and the Poisson score,

JT =

B

T

B





1X
p
1 pk
=
Xk log k + (1 Xk )
;
T k
p
1 p

(3.14)

is called the mean information gain per unit time, and

JN =

1 X
p
Xk log k ;N (T ) kp

(3.15)

is called the mean information gain per event. Vere-Jones (1999) applied the above scoringprocedure to the stress release model and the simple ETAS model.
Theoretical discussions about the above scoring procedures can be found in Vere-Jones(1999), Cziszar (1969) and Fritz (1973).

24

CHAPTER 3. MODELLING STRATEGIESPart IISPACE-TIME ETAS MODELSAND THEIR APPLICATIONS

25

Chapter 4

The Epidemic Type Aftershock
Sequence Model

4.1 Derivation
The ETAS model originates from the modi ed Omori law, which describes the decay of
the aftershock frequency as an inverse power law. For the aftershocks of the 1981 Nobi
earthquake of MS 8:0, Omori (1894) rst tried to t an exponential decay law to thedata, but no satisfactory results were obtained. Instead, he found that the number ofaftershocks occurring in the unit time interval was well tted with the equationn(t) = K (t + c) 1 ;(4.1)

where t is the elapsed time from the occurrence of the mainshock, K and c are constants.
Then Utsu (1957) found that the decay of the aftershock numbers was a bit faster than
(4.1) in a way such thatn(t) = K (t + c) p ;(4.2)

which is called the modi ed Omori law.
The the modi ed Omori law has frequently been used a lot to analyze the aftershock
activity of aftershocks. Utsu et al. (1995) reviewed the values of p for more than 200
27

28

CHAPTER 4.

THE ETAS MODEL

aftershock sequences. It ranges from 0.6 to 2.5 with a median of 1.1, and no relationshipwas found between the p-values and the mainshock magnitudes.
Some aftershock activities do not appear to be well tted to the modi ed Omori law
in the cases where strong aftershocks produce secondary aftershocks within the aftershock
sequence. Utsu (1970) and Ogata (1983) used the superposition of the modi ed Omori
law to describe those secondary aftershocks, i.e.,n(t) =nX
i=0

H (t Ti )Ki (t Ti + ci ) pk ;(4.3)

where H (t) is the unit step function such that H (t) = 0 for t  0 and H (t) = 1 for t > 0;T0 is the occurrence time of the mainshock; T1 ; T2 ; : : : ; Tn are the occurrence times ofsuch strong aftershocks, respectively; and Ki , ci and pi are constants, i = 0; 1; 2;    ; n.
Equation (4.3) indicates that the general aftershock sequence consists of the aftershocks
produced by the mainshock and the secondary aftershocks produced by strong aftershocks.
Ogata (1988, 1989) introduced a self-similarity in seismicity using the modi ed Omori
formula and suggested that not only the strong aftershocks but also each event in the
catalogue can produce its own aftershocks. He used the conditional intensity function todescribe this model, and suggested that the occurrence intensity conditional intensity ofsuch a model is given by

(tjHt ) = PrfN [t; t + dt) > 0jHt g =  +Ki;(t ti + c)p
ti <t

X

(4.4)

where the constant  represents the background occurrence rate, each term in the summation represents the contribution to aftershocks from each event to the occurrence rate,
and ti is the occurrence time of the ith event. The constant Ki depends on the magnitudeof the ith event in such a way that

Ki = K0 e (M M ) ;
0

(4.5)

4.1.

29

DERIVATION

where M0 is the reference magnitude and K0 is a constant representing productivity ofaftershocks. Equation (4.5) was hinted by the empirical formula of Utsu and Seki (1959)
log10 S = 1:02M

4:0;

(4.6)

where S is the area of the aftershock region and M is the magnitude of the mainshock.
Assuming the number of aftershocks being proportional to the area of the aftershock
region, the number of aftershock is proportional to the exponential of the the magnitude,i.e.,

N / e M;
where

(4.7)

is a constant. Thus, (4.4) becomes(tjHt ) =  + K0

f (Mi M0 )g :(t t + c)p

X exp
ti <t

i

(4.8)

Ogata called the seismicity satisfying (4.8) as the standard seismicity, and named thismodel "Epidemic Aftershock Sequence model" (ETAS model).
This model is similar to the ideas of Vere-Jones and Davies' triggering model (VereJones and Davies, 1976). In the trigger model, each event is classi ed to either triggering
event or triggered event. Only the triggering events can produce a number of o spring
and the numbers of o spring are non-negative-integer valued random variables. They
con rmed that the inverse power law (the modi ed Omori law) ts better than the exponential decay for the New Zealand earthquakes.
Ogata (1998) then adopted the following conditional intensity for the space-time ETAS
model

(t; x; y ) = (x; y ) +X
i:ti <tg(t ti ; x xi ; yyi ; Mi ):(4.9)

Here he assumed the separability of the response functiong (t; x; y; M ) = (M )g (t)f (x; y jM );(4.10)

30

CHAPTER 4.

THE ETAS MODEL

where(M ) = Ae (M M )0(4.11)

is the expected number of the aftershocks generated from a mainshock of magnitude M ,


t
p 1
1+
g (t) =cc

 p

(4.12)

is probability density function of the lagged time distribution of aftershocks, and 2

x + y2
1
f (x; y jM ) =f(4.13)
 (M )
 (M )
is the density function of the aftershock locations centred at a mainshock with a magnitudeM.
Before Ogata had extended the ETAS model to its space-time version, Musmeci andVere-Jones (1992) considered space-time di usion clustering models to analyze the historical seismicity in Italy. The conditional intensity function for their models are
Ae M e ct
1 x2 y 2g (t; x; y; M ) =exp+;(4.14)
2x y t
2t x2 y2
and

Ae M e ct t2 C Cg (t; x; y; M ) = 2 2 2 2 2 x y2 2 ;(4.15)
 (x + t Cx )(y + t Cy )where A, , x , y , Cx , and Cy are constants. For a xed point (x; y ), when t !
1, the functions in (4.14) and (4.15) decay with time according to t 1e ct and t 2 e ct,
respectively. Kagan (1991) and Rathbun (1993) also gave di erent forms for space-timecluster modelling.

4.2 De nition of space-time ETAS models
The space-time-magnitude ETAS model can be de ned by the following conditional intensity function of the form

(t; x ; M jHt ) = (t; x ; M ) +X
i:ti <th(t; x ; M ; ti ; x i ; Mi );(4.16)

4.2.

31

DEFINITION OF SPACE-TIME ETAS MODELS

where (t; x ; M ) represents the background intensity, h(t; x ; M ; ti ; x i ; Mi ) is the response
function from the i-th event de ned in the case t > ti . This model means that thebackground events is a Poisson process with a intensity function (t; x ; M ); and that once

an event of fti ; x i ; Mi g occurs, it generates a further Poisson process with the intensity ofh(t; x ; M ; ti ; x i ; Mi ).
If (t; x ; M ) is invariant in time, i.e., (t; x ; M ) = (x ; M ), the background intensity
is called stationary. If (t; x ; M ) = (t; x )ja (M ), the background intensity is called magnitude separable. Moreover, a magnitude separable stationary background rate is de nedby the relation (t; x ; M ) = (x ; M ) = (x )j (M ).
If h(t; x ; M ; ;  ; M 0 ) = h(t; x ; M ; M 0 ), we call the clustering response func-tion to be homogenous. Of course, we can de ne separately the homogenous clustering
response function in space or time. Ifh(t; x ; M ; ;  ; M 0 ) = g (tj )h0 (x ; M j  ; M 0 );we call that the clustering response function is time separable. If00h(t; x ; M ; ;  ; M 0 ) = h (t; x j;  ; M 0 )j (M jM 0 );we call that the clustering response function is magnitude separable. A space-time separable clustering response function is the case whereh(t; x ; M ; ;  ; M 0 ) = (M 0 )g (tj )f (x j  ; M 0 )j (M jM 0 ):In this chapter, the space-time separable homogenous clustering,h(t; x ; M ; ;  ; M 0 ) = (M 0 )g (t  )f (x jM 0 )j (M jM 0 );

will be specially studied.
Consider the model with a independent magnitude distributions

(t; x ; M ) = (x ; M ) = (x )ja (M );32

CHAPTER 4.

THE ETAS MODELh(t; x ; M ; ;  ; M 0 ) = (M 0 )g (tj )f (x j  ; M 0 )j (M jM 0 );andja (M ) = j (M jM 0 ) = j (M ):Then we have(t; x ; M jHt ) = (t; x jHt )j (M )where

(t; x jHt ) = (t; x ) +X 00
i:ti <th (t; x jti ; x i ; Mi ):Thus, in such a model, the magnitude distribution can be neglected, which leads us to
the space-time ETAS model.

4.3 Properties
4.3.1 Criticality
In the space-time ETAS model, the process of the immigrants or the background events is a
Poisson process with a intensity function (t; x ; M ). Once an event occurs, it produces ano spring process with the aforementioned intensity h. If we regard the background eventsas the immigrant generation G0 with the intensity function Q0 (; x ; M ) = (t; x ; M ), thenthe intensity function for the rst generation G1 is

Q1 (t; x ; M ) =ZZZSh(t; x ; M ; ;  ; M 0 )Q0 (;  ; M 0 )dd dM 0 ;(4.17)and the intensity function for the nth generation is

Qn (t; x ; M ) =ZZZSh(t; x ; M ; ;  ; M 0 )Qn 1 (;  ; M 0 )dd dM 0 :(4.18)

h[k](t; x ; M ; ;  ; M 0 )Q0 (;  ; M 0 )dd dM 0 :(4.19)

Thus we haveQn (t; x ; M ) =

ZZZ
S

4.3.

33

PROPERTIES

where

h[1](t; x ; M ) = h(t; x ; M );(4.20)andZZZ

h[k](t; x ; M ; ;  ; M 0 ) =Sh[k 1] (t; x ; M ; 1 ;  1 ; M1 ) h(1 ;  1 ; M1; ;  ; M 0 )d1d 1dM1 :(4.21)

Suppose that a(;  ; M 0 ) and b(t; x ; M ) are the left and right eigenfunctions of the response

function h(t; x ; M ; ;  ; M 0 ), respectively, corresponding to the maximum eigenvalue ,i.e.1 ,

a(;  ; M 0 ) =b(t; x ; M ) =ZZZZ Z ZSSa(t; x ; M )h(t; x ; M ; 1 ;  1 ; M1 )dtdx dM;(4.22)

h(t; x ; M ; ;  ; M 0 )b(;  ; M 0 )dd dM 0

(4.23)

with the restriction such that
ZZZ

a(t; x ; M )b(t; x ; M )dtdx dM = 1:

(4.24)

(t; x ; M ; ;  ; M 0 ) = a(;  ; M 0 )b(t; x ; M )(4.25)

S

Note thatis the projection operator of h corresponding to  satisfying the relationZZZ=

Z Z ZSS

(t; x ; M ; 1 ;  1 ; M1 )h(1 ;  1 ; M1 ; ;  ; M 0 )d1 d 1 dM1h(t; x ; M ; 1 ;  1 ; M1 )(1 ;  1 ; M1 ; ;  ; M 0 )d1 d 1 dM1= (t; x ; M ; ;  ; M 0 ):(4.26)

Considering (4.19), for the case where n ! 1, we have

h[n]
! ;
n

(4.27)

We suppose this maximum eigenvalue is isolated from other eigenvalues, i.e., there exists a positive
real number C , such that  > C >  , where  is any of other eigenvalues.
1

0

0

34

CHAPTER 4.

THE ETAS MODEL

and therefore

Qn (t; x ; M ) !

n

ZZZ
S

(t; x ; M ; ;  ; M 0 )Q0 (;  ; M 0 )dtd dM 0 :(4.28)

From (4.28), when n ! 1, we see that Qn ! 0 for  < 1, and that Qn ! 1 for  < 1.
Here  is called the critical parameter because if  < 1, the process is stable and otherwise
the population explosive to in nity.
Equation (4.28) can be rewritten as

Qn (t; x ; M ) !
=
=

n

ZZZ

b(t; x ; M )a(;  ; M 0 )Q0 (;  ; M 0 )dd dM 0S
ZZZn b(t; x ; M )a(;  ; M 0 )Q0 (;  ; M 0 )dd dM 0Sn b(t; x; M ) constant;(4.29)

indicating that b(t; x ; M ) is proportional to the asymptotic intensity of the population in
the nth generation, when n ! 1.

The eigenfunction a(;  ; M 0 ) can be interpreted as the ability of producing o spring,

including direct o spring and indirect o spring, from an ancestor of f;  ; M 0 g, because
lim

n!1

= nlim!1= nlim!1= nlim
!1

1 ZZZ
X
i=n
1
X
i=n
1
X

i

S
ZZZ

S

(t; x ; M ; ;  ; M 0 )dtdx dM

i a(;  ; M 0 )

i=n
n

1

h[k](t; x ; M ; ;  ; M 0 )dtdx dM





ZZZ
S

b(t; x ; M )dtdx dMa(;  ; M 0 )  constant:(4.30)For a space-time ETAS model with a separable magnitude distribution, a temporally
constant background rate and a space-time separable stationary and homogenous clustering response, i.e.,(t; x ; M jHt ) = (x )j (M ) +X
i:ti <t(Mi )g (t ti )f (xx i jMi )j (M jMi );(4.31)

4.3.

35

PROPERTIES

the eigenvalue equations area(;  ; M 0 ) =ZZZb(t; x ; M ) =

S

ZZZ
S

a(t; x ; M )(M 0 )g (t  )f (x  jM 0 ) j (M jM 0 )dtdx dM;(4.32)

(M 0 )g (t  )f (x  jM 0 )j (M jM 0 ) b(;  ; M 0 )dd dM 0:(4.33)

Integrate these with respect to  and  , we can prove that a(; x ; M 0 ) and b(t; x ; M ) areonly functions of M 0 and M , respectively. Namely, the equations (4.32) and (4.33) canbe simpli ed as

a(M 0 ) =andb(M ) =Za(M )(M 0 )j (M jM 0 )dM;

(4.34)

(M 0 )j (M jM 0 )b(M 0 )dM 0 :(4.35)M

Z

M

If j (M jM 0 ) = j (M ), then we havea(M 0 ) = (M 0 )andb(M ) = j (M )Za(M )j (M )dM(4.36)

(M )b(M )dM:

(4.37)M

Z

M

Thus we can see that a(M 0 ) = C1 (M 0 ) and b(M ) = C2 j (M ), where C1 and C2 are someconstants. Substituting a(M 0 ) = C1 (M 0 ) and b(M ) = C2 j (M ) into (4.36) and (4.37),
respectively. we can get

=

Z

M

(M )j (M )dM;

(4.38)a(M 0 ) = (M 0 );(4.39)b(M ) = j (M ):(4.40)and36

CHAPTER 4.

THE ETAS MODEL

Example 1: Consider the ETAS model with the conditional intensity function"#e (Mi M )(t; m) = e (M M )  + K
p ;
(
t
t
+
c
)
i
i:ti <t
X

0

where



p 1
t
g (t) =
1+
c
c
and

 p

0

(4.41)

;





c1 p
Ke (Mi M ) ;(M ) =
1 p
0

and assume that the locations are uniformly distributed on the studied area. Then,
according to our discussion above, we have the criticality parameter

=

Z 1
M0

j (M )(M ) dM =A

Example 2: Let h(M ) be a probability density function and H (M ) be the corresponding
cumulative distribution function. Set

j (M jMi ) =h(M ); 0  M < aMi ; 0 < a  1:
H (aMi )That is, the magnitudes of o springs are limited not to be greater than aMi , where Mi is
the magnitude of the ancestor of interest. Then we have, by Equation (4.35),
Z(Mi ) aMi
b(Mi ) =h(M )b(M ) dm:
H (aMi ) 0
In order to obtain the eigenvalue, we apply the limit operator to both sides of the above
equality. Notice that H (0) = 0 and h(M ) = H 0 (M ). Then we have
R aMi0b(0) = (0) Mlim
!0ih(M )b(M ) dm a(0)b(0)=:H (aMi )aThus, the criticality parameter  is

 = (0):

4.3.

37

PROPERTIES

Example 3: The coupled ETAS model treats the studied region as a combination ofseveral subregions. Each individual ancestor produces o spring in all the subregions andthe in uence to each subregion depends on the size and location of the ancestor. Denote
the rate in subregion R by (t; R; M ), where R = 1; 2; :::; N with N being the total
number of subregions. Then(t; R; M ) = (R; M ) +X X
R0 i:tR0 ;i <t(R; R0 ; MR0 ;i )g (t tR0 ;i )j (M jMR0 ;i );(4.42)

where (R; M ) is the background rate in the region R; (R; R0 ; M ) gives the in uence ofan event of magnitude M in region R0 to the region R; and MR0 ;i denotes the magnitudeof the i-th event in subregion R0 . If the interactions between subregions are independentof the magnitudes, that is, (R; R0 ; M ) = D(R; R0 )(M ) where D(R0 ; R) can be regarded
as the N  N in uence matrix, then (4.42) can be rewritten as(t; R; M ) = (R; M ) +

X
R0

2

4D (R; R0 )X
i:tR0 ;i <t3

(MR0 ;i)g (t tR0 ;i )j (M jMR0 ;i )5 : (4.43)

In this case, the eigenvalue equation is given by

V (Ri ; Mi ) =Z X
R

D(R; Ri )(Mi )V (R; M )j (M jMi )dM:

Assume thatV (M ) = fV (1; M ); V (2; M ); :::; V (N; M )gT
and de ne an N  N matrix D whose elements are

Dij = D(i; j ); i; j = 1; 2; :::N;
where N is the total number of the subregions. Then we have

 V (Mi ) = (Mi )DTZ

V (M )j (M jMi )dM:

Assuming j (M jMi ) = j (M ), substitute the form of V (M ) =B
(M ) into the aboveequality,

B


= DT

B


Z

(M )j (M )dM:

38

CHAPTER 4.

THE ETAS MODEL

Then we have

 = the maximum eigenvalues of D Z(M )j (M )dM:

4.3.2 Moment properties
First order moment. The rst moment, or the unconditional intensity, is de ned by
m1 (t; x ; M )dtdx dM = E[N (dt  dx  dM )]:(4.44)

From the de nition of the condition intensity function,

(t; x ; M jHt )dtdx dM = E[N (dt  dx  dM )jHt ];we havem1 (t; x ; M )dtd x dM = E[(t; x ; M jHt )]:(4.45)

From (4.16), we obtain

m1 (t; x ; M ) = E[(t; x ; M jHt )]"= (t; x ; M ) + E= (t; x ; M ) += (t; x ; M ) +X

i:t i<tZZZ

Z Z ZSS#h(t; x ; M ; ti ; x i ; Mi )h(t; x ; M ; ;  ; M 0 )N (dd dM 0 )h(t; x ; M ; ;  ; M 0 )m1 (;  ; M 0 )dd dM 0 ; (4.46)

where S = R2  ( 1; t)  M represents the whole region, the whole magnitude range
and the whole history up to time t.
If the model has a stationary background in time and a stationary clustering response
function in time, the unconditional intensity or the rst moment is invariant in time, i.e.,a function of only locations and magnitudes. This can be derived in the following way.
The solution for m1 (t; x ; M ) of the integral equation

m1 (t; x ; M ) = (x ; M ) +ZZZSh(t ; x ; M ;  ; M 0 )m1 (;  ; M 0 )dd dM 0 ; (4.47)

4.3.

39

PROPERTIES

if it exists, can be expanded by the so-called Neumann series (see, e.g., Kanwal R. P.,
1997),

m1 (t; x ; M ) = (x ; M ) +

1 ZZZ
X
S

k=1

h[k](t ; x ; M ;  ; M 0 )( ; M 0 )dd dM 0 ;(4.48)

whereh[k](tZZZ

; x ; M ;  ; M 0 ) =h[k 1] (t 1 ; x ; M ;  1 ; M1 )S

 h(1 ;  1 ; M1;  ; M 0 )d1d 1dM1 ;(4.49)

andh[1] (t ; x ; M ;  ; M 0 ) = h(t ; x ; M ;  ; M 0 ):(4.50)

Notice that S = R2  ( 1; t)  M, and set  0 = t  . Then we havem1 (t; x ; M ) = (x ; M ) +

1 ZZZ
X
S0

k=1

h[k]( 0 ; x ; M ;  ; M 0 )( ; M 0 )d 0 d dM 0 ;(4.51)

where S 0 = R2  (0; 1) M. The right hand side of (4.51) is independent of t, indicatingthat m1 (t; x ; M ) is independent of time t. Thus, the solution of the equation in (4.51)
can be rewritten as

m1 (x ; M ) = (x ; M ) +

1 ZZ
X
D

k=1

H [k](x ; M ;  ; M 0 )( ; M 0 )d dM 0 ;(4.52)

where D = R2  M andH [k](x; M ;  ; M 0 ) =




Z 1

0

h[k](t; x ; M ;  ; M 0 )dt:

(4.53)For a space-time ETAS model assuming separable magnitude, a temporally constantbackground rate and a space-time separable homogeneous clustering response function,such that,(t; x ; M jHt ) = (x)j (M ) +X
i:ti <t(Mi )g (t ti )f (x xi jMi )j (M jMi )(4.54)

40

CHAPTER 4.

THE ETAS MODEL

from (4.46), the unconditional intensity is given by

m1 (x ; M ) = (x )j (M ) += (x )j (M ) +

ZZZ
ZZ S
S

(M 0 )g (t  )f (x  jM 0 )j (M jM 0 )m1 ( ; M 0 )dd dM 0(Mi )f (x  jM 0 )j (M jM 0 )m1 ( ; M 0 )d dM 0 :(4.55)

Integrating both sides of (4.55) with respect to space, we obtain

(Mi )j (M jM 0 )mj
 t (M 0 )dM 0 ;(4.56)

j (M )+ (Mi )j (M jM 0 )jt (M 0 )dM 0 ;jt (M ) = m

M(4.57)mj t (M ) = j (M ) +that is to say,

Z

M
Z

where

m
=

Z

M

Z

dM

and

=

m1 (x ; M )dx ;(4.58)

(x )dx :(4.59)

R2

Z

R2

Equation (4.57) gives the relationship between the total magnitude distribution, the background magnitude distribution and the aftershock magnitude distributions.
Integrate again (4.57) with respect to magnitude. Then we havem
 = +m


Z

or

m
=

M

(M 0 )jt (M 0 )dM 0 ;;
1 E[(M )](4.60)(4.61)

where E[(M )] is interpreted as the average number of direct o spring generated by anancestor. Moreover, if the magnitude distribution is independent from the history ofmagnitudes, i.e., j (M ) = jt (M ) = j (M jM 0 ) = j (M ), thenE[(M )] =ZM(M )j (M )dM;

which is the same form as the critical parameter .

(4.62)

4.3.

41

PROPERTIES

Second order moment. The second order moment density formally is de ned by
m2 (t1 ; x 1 ; M1 ; t2 ; x 2 ; M2 ) =E[N (dt1 d  1 dM1 )N (dt2 d  2 dM2 )]dt1 d  1 dM1 dt2 d  2 dM2:(4.63)The de nition is derived from the following relation such that

E N (dt1 d  1 dM1 )N (dt2 d  2 dM2 ) =8< m 1 ( t1 ; x
 1 ; M1 )Æ (t1t2 ; x 1 x 2 ; M1 M2 )dt1 d  1 dM1 dt2 d  2 dM2 ; for t1 = t2 ;for t1 6= t2 ;(4.64): m 2 ( t1 ; x
 1 ; M1 ; t2 ; x
 2 ; M2 )dt1 d  1 dM1 dt2 d  2 dM2 ;where Æ (t1t2 ; x 1

x 2 ; M1

M2 ) is the Dirac function.

Similarly, the expectation in the above de nition also holds with respect to the history.
When t1 6= t2 , we have the followings:

m2 (t1 ; x 1 ; M1 ; t2 ; x 2 ; M2 )= E

(t2 ; x 2 ; M2 ) +ZZZh(t ; x2 ; M ; ; x; M 0 )N (ddxdM 0 )
2 

S





2






N (dt1 d 1 dM1) =dt1d 1 dM1= E (t2 ; x 2 ; M2 )N (dt1 d 1 dM1 ) =dt1 d 1 dM1+E

Z Z ZSh(t ; x2 ; M ; ; x; M 0 )N (dd x dM 0 )N (dt d dM ) =dt dx1 dM
2 

2





1 1

1

1 

 ; M1

M 0)

1= m1 (t1 ; x 1 ; M1 )(t2 ; x 2 ; M2 )
+

+

ZZZ

S

ZZZ
S

h(t2 ; x 2 ; M2 ; ; x ; M 0 )m1 (t1 ; x 1 ; M1 )Æ (t1; x 1 ddx dM 0h(t2 ; x 2 ; M2 ; ; x ; M 0 )m2 (;  ; M 0 ; t1 ;  1 ; M1 )ddx dM 0

= m1 (t1 ; x 1 ; M1 )(t2 ; x 2 ; M2 ) + h(t2 ; x 2 ; M2 ; t1 ; x 1 ; M1 )m1 (t1 ; x 1 ; M1 )
+

ZZZ

S

h(t2 ; x 2 ; M2 ; ; x ; M 0 )m2 (;  ; M 0 ; t1 ;  1 ; M1 )ddx dM 0 :(4.65)

42

CHAPTER 4.

THE ETAS MODEL

4.3.3 Spatial densities
When the process has a space-time separable homogeneous clustering response function
and an independent magnitude distribution such thath(t; x ; M ; ;  ; M 0 ) = (M 0 )g (t  )f (x jM 0 )j (M );(4.66)

the mean rate in (4.46) is of a convolution. Such a form is of special interest, because
all the forms suggested by Ogata (1998) are of this form and they can be investigated byusing the Fourier transformation.
R
If f (x ) denotes 0+1 (M 0 )f (x jM 0 )j (M 0 )dM 0 , which is called the average transfer
function or the average response function, then we can rewrite (4.46) in the form

m1 (x ) = (x ) +or(x ) = m1 (x )

ZZ
R2

ZZ
R2

f (x



f (x



 )m1 ( )d ;(4.67) )m1 ( )d ;(4.68)The background rate can be estimated by using (4.68) when the total rate m1 (x ) andthe average transfer function f (x ) are known, as will be discussed in next section. The
Fourier transform of (4.67) isI (! ) = U (! ) + F (! )I (! );(4.69)

where I , U and F are Fourier transform of m1 ,  and f respectively. Equation (4.69) is
important because it converts the convolution part in equation, (4.67) into a product.
From (4.69), we have

U (! ) = I (! )[1 F (! )];(4.70)

which represents the relationship between the background rate, the total rate and theaverage transfer function.

4.3.

43

PROPERTIESThe average transfer functions in the space-time ETAS model and their Fouriertransforms
I will discuss the spatial intensity of the three forms of space-time models suggested
by Ogata (1998) in this subsection. Consider that these three forms have space-time
separable clustering response functions and a background intensity constant in time (seeEquations (4.9) to (4.13)). Then we have the Fourier transform for each of the above
three branching structures below.

Form 1. In Form 1, Ogata considered rst that
(M ) = Ae M ;1expf (x; y jM ) =2de M
x2 + y 2
;
2de M

and

j (M ) = e

M:

The average transfer function for this case is given by
Z +1
f (x; y ) =(M )f (x; y jM )j (M )dM
0


Z +1
x2 + y 2Aexp
e M dM:
=
M
2
d
2
d
e
0
Its Fourier Transformation has the form:Z Z Z +1
F (!1 ; !2 ) =f (x; y jM )j (M )dM e i(! x+! y) dxdy
R
0



Z Z Z +1
A
x2 + y 2M=
exp
e dM e i(! x+! y) dxdy
(
M)
2
d
2
d
e
R
0


Z +1
1 M 2
(
)
M
2
de (!1 + !2 ) dM:
=A
e
exp
2
0
1

2

2

1

2

2

In the above integral, the Fourier transform of the standard normal density function,
expf  2 ! 2 =2g is used. Setting(!12 + !22 )e M d
= ;
2

44

CHAPTER 4.

we haveA

F (!1 ; !2 ) =THE ETAS MODEL

for !12 + !22 = 0;and
 = +1 Z +1
2
F ( !1 ; ! 2 ) = e  d22(!1 + !2 )d(! +! )d=2
 2



=
1
A (!1 + !22 )d(!12 + !22 )d
=
+ 1;
2
2


A

2
1

2
2

for !12 + !22 > 0, where (a; x) is the incomplete gamma function.
Before we discuss Form 2 and 3, we give a useful integral here,
Z +10 a q 1 K (ax)tJ0 (at)q 1
dt=(t2 + x2 )q
2x
(q )

1
for q > ;2

(4.71)

where J0 (x) is the Bessel function of the 0-th order and Kq (x) is the modi ed Bessel
function of the q -th order. This integral can be proved in the following way. Set

I = (q )Z +10

tJ0 (at)dt;(t2 + x2 )q

and consider the relation(q ) =

Z +1

qe

0

d = (t2 + x2 )q

Z +1

0

q 1 e (t2 +x2 ) d :Then we have

I =
=
Using

Z +1

0

Z +1

tJ0 (at)dt

0

q 1 exp[

Z +1

0

Z +1

0

x2 ]d

q 1 exp[

Z +1

0(t2 + x2 ) ]dtJ0 (at) exp[ t2 ]dt:

1 a
exp[ p2 t2 ]J0 (at)tdt = 2 e p ;
2p
2

4 2

(see Tranter 1968, page 75) or
Z +1

0

exp[

t2 ]J0 (at)tdt =

1 a
e ;
2
2

4

4.3.

45

PROPERTIES

we obtain


Z



a
1 +1 q 2
exp x2
d
I =
2 0
4



 

Z 
1
2x
a
2x
1  a q 1 1 2x q 2
=
exp
ax
+
d
2 2x
a
2
a
2x
a
0



Z 1


q
1
1 a
1
1
=
q 2 expax  +
d
2 2x
2

0
 a q 1
1=
Kq 1 (ax) for q > :2x
2The proof of the integral
Z 10

q 2 exp1
1
ax  +2d = 2Kq 1 (ax); for q >1
2

can be found in Page 75 of the book by Tranter (1968).

Form 2. In Form 2, Ogata then considered the second case where(M ) = A exp[ M ];f (x; y jM ) =B;(x2 + y 2 + d)q

andj (M ) = exp[

M ];where B is the normalising factor,B=(q1)dq 1:For this case we have
Z +1

ZZ



f (x; y jM )j (M )dM e (! x+! y) dxdy
R
0
ZZ
Z +1
e i(! x+! y)
B 2 2
A e ( )M dM=
q dxdy
(
x
+
y
+
d
)
R
0
ZZ
i
(
!
x
+
!
y
)
AB
e
=
B 2 2
dxdy:
(x + y + d)q
R

F (!1 ; !2 ) =

1

2

1

2

2

1

2

2

2

46

CHAPTER 4.

THE ETAS MODEL

Setting x = r cos x, y = r sin x, !1 = p cos , !2 = p sin , we have
Z 2 Z +1
re ipr cos( )ABdrdF (!1 ; !2 ) =(r2 + d)q
0
0
Z 2
Z +1
r
AB
ipr
cos(


)
e
d
dr
=
2
0Z
0 (r + d)q+1 rJ0 (pr)
AB
= 2
dr
(r2 + d)q0
pp q 1 Kq 1 (p d)
AB
p
= 2
;
(q )
2 d
p
where p = !12 + !22 and Kq 1 (x) denotes the modi ed Bessel function.

Form 3. In Form 3, Ogata nally considered the case where(M ) = Ae M ;f (x; y jM ) = Be (q 1)M x2 + y 2 + e M d q ;andM;j (M ) = e
where B is the normalising factor,B=(q1)dq 1:


The we have
Z +1

ZZ



f (x; y jM )j (M )dM e (! x+! y) dxdy
R
0

Z Z
Z +1
e i(! x+! y)(q)
M
dxdy dM:=Ae
B 2 2(x + y + e M d)q0
R

F (!1 ; !2 ) =

1

2

2

1

2

2

Setting x = r cos x, y = r sin x, !1 = p cos , !2 = p sin , we obtainZ 2 Z +1
Z +1
re ipr cos( )(q)
MF (!1 ; !2 ) = ABedrd dM(r2 + e M d)q
0
0
0
Z 2

Z +1
Z +1
r
(
q
)
M
jpr
cos(


)
= AB
e
e
d
dr dM
0Z
0 Z
0 (r2 + e M d)q+1+1 rJ0 (pr)dr dM
= 2ABe ( q )M00 (r2 + e M d)q p

q 1
Z +1
p
Kq 1 (p e M d)
(
q
)
M
p M
= 2AB
e
dM:
(q )
2 e d
0

4.4.

47

SIMULATION ALGORITHMSp

Setting e M d =  , we have  dM = 2d and
Z +1 
2AB2
 2q 2 =  p q 1
p
F (!1; !2 ) =Kq 1 (p ) d
p
(q )
2

d Z
d


+
1
q
1
4AB p=d = q p  q 2 = Kq 1 (p )d:(q ) 2d

4.4 Simulation algorithms
Here we consider the simulations for a space-time ETAS model with the intensity function

(t; x ; M jHt ) = (t; x ; M ) +X
i:ti <th(t; x ; M ; ti ; x i ; Mi );(4.72)

in a time interval [0; T ] and a region S . Except for the two basic techniques available forsimulation of a usual point process, i.e., the inverse and thinning methods, the composition
method seems to be the most eÆcient one for simulating an ETAS point process.
Since the simulation of the ETAS model by using the inverse method appears the mostineÆcient method, this section gives only the simulation algorithms based on the thinning
method and the composition method.

4.4.1 Thinning method
The thinning method for the simulation of point processes can be found in Lewis andShedler (1979), Ogata (1981, 1998), Musmeci and Vere-Jones (1992).
The following algorithm is a modi ed version from Ogata (1998) that uses the superposed intensity with respect to locations and magnitudes

0 (t) = 0 (t) +where

0 (t) =

Z

M

dM

Z
S

X

i (t);(4.73)

(t; x ; M )dx ;(4.74)

i: ti <t48
andi (t) =

Z

M

Z

dM

S

CHAPTER 4.

THE ETAS MODELh(t; x ; M ; ti ; x i ; Mi )dx ;(4.75)

and also the superposed intensity with respect to magnitudes

1 (t; x ) =ZM

(t; x ; M jHt )dM:

(4.76)

Algorithm 4.1 (Thinning method)
1. Set i = 0 and t0 = 0..
2. Set  = ti and nd a constant C such that 0 (tjHti ) < C holds for the interval[ti ; 1).
3. Starting from  , generate an event according to the stationary Poisson process witha mean rate C on the region S . Record its occurrence time as t0 .
4. If t0  T , stop.
5. Generate a uniform random variable U on [0; 1].
6. If U < 0 (t0 )=C , then set ti+1 = t0 and i = i + 1, and go to next step; otherwise, set

 = t0 and go to Step 3.
7. Generate the location x i of the ith event using the superposed intensity 1 (ti ; x ) in
(4.76).
8. Generate the magnitude Mi of the ith event according to the conditional intensity
function (ti ; x i ; M jHti ).
9. Record the coordinates of the ith event (ti ; x i ; Mi ); set i = i + 1; and then go to
Step 2.

4.4.

49

SIMULATION ALGORITHMS

Simulation by using this algorithm spends much of the CPU time in computing thesummation in (4.72) and the integrals in (4.73) and (4.76), especially when we calculatethem numerically. Another shortcoming of this algorithm is that the in uence of theevents outside of S to the occurrence rate (the boundary e ect) is not taken into account.

4.4.2 Composition method
The composition method is based on the branching-process interpretation of the ETAS
model, by which each events in the process generates its own descents independently.

Algorithm 4.2
1. Generate a Poisson process of intensity (t; x ; M ) in the time interval [0; T ] andregion S by using Algorithm 3.1, recored as G(0) = f(t(0) ; x(0) ; M (0) ) : i = 1; 2;   g,
i

i

i

i.e., population in Generation 0;
2. Set l = 0;
3. For each event i in G(l) , with (t(il) ; x (il) ; Mi(l) ), generate a non-stationary Poissonprocess of intensity h(t; x; M ; t(l) ; x(l) ; M (l) ) in the time interval [0; T ] and region S


i

i

i

by using Algorithm 3.1, and record the set of all the newly generated events as G(il) ,namely, the set of the events generated by the ith event in Generation G(l) .
[4. Set G(l+1) = G(il) .
i

5. If G(l+1) is the empty set, then set l = l + 1 and go to Step 3; otherwise, go to nextstep.
6. Set G =[lG(l) and sort all the events G in the ascending order of time. Then G isa simulated catalogue of the events.

50

CHAPTER 4.

THE ETAS MODEL

As we discussed before, this method may produce a slight less number of events owing
to the ignorance of the boundary e ects in both time and space. To avoid such a problem,
we can use an enlarged region and time interval for the simulation to select the events inthe target region and time interval to form the simulated catalogue.
When the model has space-time separable homogeneous clusters and an independent
magnitude distribution such that"(t; x ; M jHt ) = j (M ) (x ) +X
i:ti <t#(Mi )g (t  )f (x jMi ) ;(4.77)

Step 3 can be simpli ed as follows:
30 .

For each event in G(l) , for example, (t(il) ; x (il) ; Mi(l) ), generate a Poisson randomvariable N (l) of mean (M (l) ); then generate N (l) random variables, t ,    , t l ,
i

i

i

i

Ni( )

from the probability density function g (), and Ni(l) random variables, x 1 ,    , x Ni lfrom the density f (j). Then generate Ni(l) random variables, M1 ,    , MNi l fromthe density j (). Record the set of all the newly generated events as G(l) ; namely,
( )

( )

i

the set of the events generated by the ith event in Generation G(l) .

4.5 Prediction performance
Here is an examples for the prediction performance of the ETAS model. We simulated acatalogue according to the ETAS model of a conditional intensity function,

(t; M jHt ) = exp(
with parameters  = 1,

"



t ti
A(p 1) X
exp( M ) 1 +M)  +
cc
i:ti <t p#(4.78)= 0:8 log 10  1:8421, A = 0:16, c = 0:001 = log 10 = 2:3026

and p = 1:2. The same parameters were also used by Helmstetter and Sornette (2002) for
simulating a synthetic catalogue and testing the prediction performance. The simulated
catalogue is shown in Figure 4.1. The total number of simulated events is 257,618 for theinterval [0; 60000].

51

PREDICTION PERFORMANCE

(a)

6
3

4

5

Magnitude

7

8

9

4.5.

0

10000

20000

30000

40000

50000

40000

50000

40000

50000

60000

9.0

Time in days

7.5
6.0

6.5

7.0

Magnitude

8.0

8.5

(b)

0

10000

20000

30000

60000

100000 150000 200000 250000

(c)

0

50000

Cumulative numbers

Time in days

0

10000

20000

30000

60000

Time in days

Figure 4.1: Illustration of the simulated catalogue.
(a) Plot of magnitudes against times for all the events; (b) Plot of magnitudes against
times for the M  6:0 events; (c) Cumulative frequencies against times.

52

CHAPTER 4.

THE ETAS MODEL

20

(a)

10

12

0.6

0

0.0

2

4

0.2

6

8

0.4

Probability

14

16

0.8

18

1.0

M>=5.0

30000

35000

40000

45000

50000

55000

60000

50000

55000

60000

50000

55000

60000

Time in days

8

(b)

6
4

0.2

0

0.0

2

0.1

Probability

0.3

0.4

M>=6.0

30000

35000

40000

45000
Time in days

(c)

0.03

0

0.00

0.01

0.02

Probability

0.04

1

0.05

0.06

M>=7.0

30000

35000

40000

45000
Time in days

Figure 4.2: Prediction probabilities for di erent magnitude classes(a) M  5:0; (b) M  6:0; and (c) M  7:0. Two dished lines are the predictedprobabilities by using the Poisson model and the background intensity, respectively. The
bar plot indicating the number of target events, marked by the right hand axis.

4.5.

53

PREDICTION PERFORMANCE

The prediction is carried out for 50-day future intervals starting from the 30000th day
to the 60000th day. Following the prediction and scoring procedures in Section 3.4, forthe kth interval, I carried out 1000 times of simulation and took the proportion of thesimulations with at least one target event occurring to all the simulation as the predictedprobability, pk , by the ETAS model. The corresponding predicted probability by the
Poisson model is (M Mc ) ];p = exp[ ewhere  is the mean rate, M is the minimum magnitude to predict and Mc is the magnitude threshold. The predicted probabilities for the simulation data are shown in Figure
4.2. Then I calculated the total scores for the ETAS model over the intervals that haveat least one of the target event (event intervals), and for the non-event intervals, de ned
as

B1 =
and

X

Xk log pk(4.79)(1 Xi ) log(1 pK );(4.80)

k

X

B2 =

k

respectively, where Xk = 1 indicating a success and Xk = 0 indicating a failure. The
same scores are applied to the Poisson model, namely

B10 =
and

B20 =

X

Xk log pk(4.81)(1 Xk ) log(1 pk );(4.82)

kXkover the event intervals and the non-event intervals, respectively. Thus, the information
gain over the event intervals is given by

GE = B1

B10(4.83)and the information gain over the non-event intervals is given by

GN = B2

B20 :(4.84)

54

CHAPTER 4.

THE ETAS MODEL

Table 4.1: Prediction scores for a simulated data by using ETAS model

M
NI NE
GE
GN
GT
M  7:099-3.570301
5.854142
2.283841
M  6:5 22
28 -12.276903 19.406773
7.129871
M  6:0 73 100 -42.688716 60.792308 18.103591
M  5:5 193 321 -101.838535 164.126138 62.287604
M  5:0 405 1029 -126.807137 272.969109 146.161971

GP E-0.396700-0.558041-0.584777-0.527661-0.313104

The meanings of the symbol are de ned as following:M : thePtarget minimum magnitude;
NI = k Xk : number of intervals with events occurring in it (event intervals);
NE : number of events;
GE : information gain over the event intervals;
GN : information gain over the non-event intervals;
GT = GE + GN : information gain over all the intervals;
GP E : information gain per event interval.
The information gain per event interval is de ned as

GE:
GP E = Pk Xk(4.85)The results are shown in Table 4.1 and Figure 4.2. From the table we can that theoverall performance of the ETAS model is much better than the Poisson model. But for
the intervals containing at least one event, the prediction performance is worse than the
Poisson model. As pointed out by Vere-Jones (1999), the reason is as simple as that formost of the event-intervals, the conditional intensity starts at the background intensity.
After an event occurs, the intensity that it excites decay o quickly. But the Possoin
model su ers from the over-estimation of the probability in the non-event intervals, andget a lower score in the overall performance. To make an e ective prediction with thismodel, we need use a much shorter interval or update the forecast immediately after theoccurrence of each event.

Chapter 5

Stochastic Declustering

5.1 Why stochastic declustering?
The declustering problems arise from the prediction of large earthquakes using catalogues
including clusters of aftershocks. The earthquake clusters complicate the statistical analysis of the background seismic activity that might be related to stress changes in thetectonic eld. The cluster features di er from place to place, and lie in between two
extreme types of spatial clusters of earthquakes: those which eventually decrease in timesuch as aftershock sequences and swarms, and those which persist in time at the same location. Ultimately, the persistent background activity prevails over an aftershock activity.
To predict the location of the large earthquakes, it is necessary to analyze the backgroundseismicity, for which seismologists consider that elimination of temporal cluster members
is of central importance.
A number of methods have been proposed for declustering a catalogue to obtain the
background seismicity. Most of these methods remove the earthquakes in a space-time
window around a large event called the mainshock; the di erences among these methods
relate to the choice of the window sizes and some other detailed procedures: for example,
see Gardner and Knopo , 1974; Kellis-Borok and Kossobokov, 1986; Utsu, 1969. In
general, the bigger is the magnitude of the mainshock, the larger is the window. An
55

56

CHAPTER 5.

STOCHASTIC DECLUSTERING

alternative to the window-based declustering methods is the link method based on a
space-time distance between the events (e.g., Reseanberg, 1985; Frohlich and Davis, 1990;
Davis and Frohlich, 1991).
All of those conventional methods contain arbitrary parameters, either for de ning
the sizes of the aftershock window or the link distances, in both space and time. Eachchoice of parameters gives a di erent declustered catalogue, and so a di erent estimateof the background seismicity. A trouble with the conventional declustering algorithms isthus to make an optimal choice of the parameters. For a data set from a speci c (rather
narrow) seismic region, and for a given threshold magnitude, these choices are usually
tuned by the user's subjective impression of the declustered output. Also, given a dataset, the concept of what constitutes an aftershock has not been uniquely de ned. These
are the reasons why so many declustering algorithms have been proposed.
To avoid such diÆculties, it seems to me inevitable that the declustering method should
be based on a stochastic model to objectively quantify the observations, so that any given
event has a probability of being either a background event or an o spring (cluster) event
generated by others. Therefore the main task is to estimate this probability for each event.
There are now some papers which propose space-time models describing how earthquakes
generate clustered events. For Italian historical data, Musmeci and Vere-Jones (1992)
assumed that the background seismicity was proportional to a kernel estimate obtained
from the spatial locations of all earthquakes that include clusters. Ogata (1998) suggested
a space-time version of the epidemic type aftershock sequence model. Before tting thespace-time ETAS model to all the events, he used a conventional declustering algorithm forthe preliminary estimation of the background rate by tting bi-cubic B-spline functions.
However, as we will see later, the use of such a preliminary declustering method can be
avoided.
In the following I suggest an approach to objective methods for estimating the back-

5.2.

57

THINNING PROCEDURE

ground intensity. Several technical methods are discussed for carrying out such estimation,such as the thinning (random deletion) method, and a variable bandwidth kernel function method associated with the maximum likelihood estimation of the space-time ETAS
model. Then the background intensity is obtained by using an iterative algorithm, andalso the declustered catalogues can be generated.

5.2 Thinning procedure
The problem of estimating the background intensity can be solved if we can nd out
a method to realize or simulate a process which has the same intensity as the unknown
function (t; x ; M ) in (4.16). We make use of the following thinning operation (or random

deletion). Suppose that space-time events f(tk ; xk ; yk ) : k = 1; 2;    ; N g are associatedwith a series of probabilities fk : k = 1; 2;    ; N g. Delete each event k of this processwith the probability k . Then the remaining events compose a new point process called
the thinned process. The simplest form of the thinning operation is to delete each pointwith a xed probability (e.g., Daley and Vere-Jones, 1972). The thinning method has
also been used for the simulation of point processes (for examples, see Lewis and Shedler,
1979; Ogata, 1981 and 1998; and Musmeci and Vere-Jones, 1992).
Here we interpret k as the probability of the kth earthquake being an o spring (aftershock) in the process, de ned as follows. In the second term in (4.16), we set

i;k = Prfthe kth event is an o spring of ith eventjHtk gh(tk ; x k ; Mk ; ti ; x i ; Mi ):=(tk ; x k ; Mk jHtk )andk = Prfthe k-th event to be an o springjHtk g =

k 1
X
i=1

i;k ;

(5.1)

(5.2)

where k  1 since it represents the ratio of the sum term in (4.16) to the whole of (4.16).

58

CHAPTER 5.

STOCHASTIC DECLUSTERINGThe probability that the k-th event belongs to the background is therefore(tk ; x k ; Mk jHtk );'k = Prfthe k-th event is an immigrantjHtk g = 1 k =(tk ; x k ; Mk jHtk )(5.3)

If we delete the k-th event in the process with probability k for all k = 1; 2;    ; N ,then the thinned process should realize a non-homogeneous Poisson process of the spatialintensity (t; x ; M ) (see Ogata, 1981, for the justi cation). We call this process the
background sub-process, and the complementary sub-process the cluster sub-process or
the o spring process.

5.3 Variable kernel estimates of seismicity
For the declustering purpose, people usually assume that the background earthquake
intensity is constant in time, i.e., (t; x )  (x ). If this process is also stationary andergodic, then the total spatial intensity function ( rst order moment density) is equal toZ

1 T
m1 (x ) = Tlim
(t; x jHt ) dt;!1 T 0(5.4)

where T is the length of observation period. Replacing the limit in (5.4) by a nite
approximation, and substituting (4.16) in (5.4), we obtain
Z

(

)

X
1 T
(x ) +h(t; x ; M ; ti ; x i ; Mi ) dt
m1 (x ) 
T 0
k:tk <T
1 X= (x ) +h(t; x ; M ; ti ; x i ; Mi ):T k:t <T

(5.5)

kThus we have(x )  m1 (x )1 X
h(t; x ; M ; ti ; x i ; Mi ):T k:t <T

(5.6)

k

In practice, the above estimate for the mean spatial intensity can be carried out by
several methods, for example, using kernel estimates (e.g., Vere-Jones, 1992) or spline
functions (e.g., Ogata and Katsura, 1988). Here I would like to adopt the kernel estimate.

5.3.

59

VARIABLE KERNEL ESTIMATES OF SEISMICITY

However, the simple kernel estimate with a xed bandwidth has a serious disadvantage:
for a spatially clustered point data set, a small bandwidth gives a noisy or variable estimate
for the sparsely populated area, while a large bandwidth mixes up the boundaries between
the densely populated areas and the sparsely populated areas. Therefore, instead of thekernel estimate

N
1X
Z (x
m
^ 1 (x; y ) =T k=1 d x k )(5.7)

where Zd (x ) denotes the kernel function, for example, the 2 dimensional gaussian kernel

Zd (x; y ) =1 x dy
e
2d

2+ 2
2 2

with a xed bandwidth d, we adoptN1Xm^ 1 (x ) =Z (xT k=1 dk x k );(5.8)

where dk represents the varying bandwidth calculated for each event k in the following
way. Given a suitable integer np in between 10 and 100, nd the smallest disk centered at
the location of the k-th event which includes at least np other earthquakes and is larger
than some minimum value (e.g., a distance within 0.02 degree, which is of the order ofthe location error) and let this minimum radius be dk (e.g., Silverman, 1986, Chapter 5).
Similar ideas can be found in Choi and Hall (1998) or Musmeci and Vere-Jones (1986).
Using the same kernel function Zdk in (5.8), the spatial occurrence rates of the cluster
sub-process can then be estimated by the statistics^(x ) =1X
 Z (xT k k dk x k )(5.9)

where k is derived from (5.2) and spatial mean rate for the background subprocess by^(x ) = m^ 1 (x )^ (x ) =

1X
(1
T k

k )Zdk (xx k );(5.10)

respectively. We call the estimates in (5.9) and (5.10) the variable weighted kernel estimates.

60

CHAPTER 5.

STOCHASTIC DECLUSTERING

5.4 Iteration algorithm
In Sections 5.2{5.3, we have outlined methodological solutions to the following issues: (i)how to obtain the MLE of the space-time ETAS model for the branching structure when
the background rate is given; (ii) how to estimate the total spatial occurrence rate; and(iii) how to estimate the background seismicity in the case where the branching structure
is provided. However, actually, we know only the earthquake data, namely, the records ofthe occurrence times, the locations and the magnitudes of the observed events. In orderto estimate the background rate from this data only, we carry out the following iterativealgorithm that simultaneously estimates the background rate and the branching structure.

Algorithm 1
1. Given a preliminary parameter np , say 20, and calculate the bandwidth dk for each
event (tk ; x k ; Mk ), k = 1; 2; ::; N ;2. Set l = 0, and u(0) (x ) = 1.
3. By using the maximum likelihood procedure described e.g. in Ogata (1998), t theconditional intensity functionX(t; x jHt ) = u(l) (x ) +h(t; x ; M ; tk ; x k ; Mk )k: tk <tto the earthquake data.

(5.11)

4. Calculate k from Equations (5.1) and (5.2), for each k = 1; 2;    ; N .
5. Calculate ^(x ) from Equation (5.10), and record as u(l+1) (x ).
6. If max(x;y) ju(l+1) (x ) u(l) (x )j > ", where " is a small positive number, then set
l = l + 1, and go to step 3. Otherwise take ^u(l+1) (x) as the background rate, and


stop.

5.5.

61

EXAMPLES

Table 5.1: Convergence steps of the evaluation procedure in applying the space-time
ETAS model by Algorithm 1 to the earthquake data from the central New Zealand region
1
2
3
4
5
6
7
8

log L
-7987.49-7716.51-7710.75-7696.38-7695.60-7694.69-7694.88-7694.84


0.43594
1.25850
0.95684
1.009450.99755
1.000110.99955
0.99966

A
0.433780.31482
0.339440.333700.33492
0.33465
0.33471
0.33470c0.024737
0.016264
0.017698
0.017370
0.017441
0.017426
0.017429
0.017428

0.735540.93397
0.88535
0.89776
0.895180.89574
0.89562
0.89565p
1.1573
1.16711.16381.1644
1.1643
1.1643
1.1643
1.1643d0.39395e-02
0.19980e-02
0.23097e-02
0.22326e-02
0.22487e-02
0.22452e-02
0.22459e-02
0.22458e-02

The tuning of np can be accomplished by putting an outer loop of Algorithm 1, which
enumerates the values of np to nd the one able to minimize the di erence between the
low frequency part of IT (x; y ) and ^ (x; y ). In this study, we nd that the adjustment ofR

np is not so sensitive for the nal estimates of m
^ 1 (x; y ), comparing to 0T ^ (t; x; y )dt. The

estimates only change slightly when np changes in the range 15  100. In this paper,

np = 20 is adopted as a standard value for use with New Zealand and Japanese data.
This algorithm converges quickly. Table 1 shows an example of the convergence stepsof this algorithm when the space-time ETAS model of Form 1 in Section 4.4.3 was applied
to the earthquake data from the central New Zealand region, which will be discussed in
the next section.

5.5 Examples
In these applications, we make use of the speci c function forms. These constitute Form 1
in Section 4.4.1 of the space-time epidemic-type aftershock sequence (ETAS) model, i.e.,(M ) = Ae (M M ) ;pp 1
for t > 0;
g (t) = (p 1)c 0 (t + c) otherwise;
0(5.12)(5.13)

62

CHAPTER 5.

and
1
exp
f (x; y jM ) =2de (M M )STOCHASTIC DECLUSTERING
1 x2 + y 2:
2 de (M M )0(5.14)0

Given an estimated intensity function u(x; y ), we set(x; y ) = u(x; y )for the background rate, where  is a positive valued parameter. The the maximum^ ^ ; c^; p^; d^).
likelihood method is used to estimate the parameter vector ^ = (^ ; A;

5.5.1 Central New Zealand region
The data used in this calculation is from the New Zealand local catalogue recorded by
the Institute of Geology and Nuclear Sciences. For this study, we select the earthquake
data for the period January 1970  August 1999 from the rectangular area 38Æ  43Æ S
and 171Æ 179Æ E, the central New Zealand region, which extends from the boundary of
Bay of Plenty to Authur's Pass (see Figure 5.1). We take magnitude threshold ML =4.0,

38oS

172o

174o

176o

178oE

43o

42o

41o

40o

40o
41o
42o
43o

170o

(b)

39o

(a)

39o

38oS

and consider shallow events down to the depth of 40 km.

’70

’75

’80

’85

’90

’95

Figure 5.1: Seismicity in the central New Zealand region.
(a) Epicentral locations and (b) space-time plots for the latitudes against the times of all
the events.

5.5.

63

EXAMPLES

The nal MLE values of the space-time ETAS model for this sub-region are A^ =0:335 event=(deg 2  day ), ^ = 0:896, c^ = 0:0174 day , p^ = 1:164, and d^ = 0:00225=deg 2.
From here on, we use deg 2 for the unit of the geographical area, where the unit degree isthat of latitude (111.11 kilometers). The converging steps of the calculation are shown in
Table 1.
Figure 5.2(a) shows the temporal change of the functionRR^ (t; x; y jHt )dxdy
p(t) = SRR;(5.15)^(x; y )dxdy
SRRi.e., the integrated conditional intensity S ^ (t; x; y jHt)dxdy over the whole study region

S relative to the overall background rateRR^(x; y )dxdy . The pattern is similar to thatSfrom the simple ETAS model (Ogata, 1988). Figure 5.2(b) shows the total intensitym^ 1 (x; y ) of the central New Zealand Region. The spatial clustering intensity ^(x; y ),calculated by Equation (5.9), is shown in Figure 5.2(c). To enhance the contrast between
Figures 5.2(b) and (c), we introduce the concept of relative clustering e ect, which isshown by the contours in Figure 5.2(d). This is to measure the cluster e ect relative tothe overall rate, and de ned by the ratio of the spatial clustering intensity to the total
intensity, i.e.,Pj kdj (x xj ; y yj )^ (x; y )= Pj:
C (x; y ) =m^ (x; y )j kdj (x xj ; y yj )C (x; y ) can be regarded as a smoothed value of j .

(5.16)

From Figures 5.2(b) to (d), we see that most places in the study area have a total
intensity of less than 50 events=deg 2 throughout the whole period of 28.75 years (T =10500 days). The clustering intensity is less than 20 events=deg 2 throughout the period,
and the relative clustering e ect is less than 0.3. Some spots, however, get high values
for these three quantities. They are mainly located in the plate boundary between the
Paci c and Australian Plates or in the inland volcanic zones.
By subtracting the clustering intensity from the overall intensity as indicated in Equation (5.10), we can get the estimate of the background intensity function shown in Figure

64

(a)

10

100

STOCHASTIC DECLUSTERING

1

Relative intensity

CHAPTER 5.

’70

’72

’74

’76

’78

’80

’82

’84

’86

’88

’90

’92

’94

’96

’98

50

10

100

5

100

50

50

172o

(d)

170o

174o

176o

0.2

0.4

0.4

0.2
0.2
0.4

0.6
0.2

178o E

0.2
0.6

0.4

0.2

0.4
0.4

0.2

0.4

0.6
0.8

172o

174o

176o

178o E

5

5 20

50

50

20

10

20 10

10

5

100
20

50

170o

10

50

5

1

20

172o

174o

176o

178o E

20

1

(e)

170o

20

10

50

50

50

100

20

1

5
50

20

50

20

(c)

1

20

43 o 42 o 41 o 40 o 39o 38oS

50

50

170o

43 o 42 o 41 o 40 o 39o 38oS

50

5

10

10

43 o 42 o 41 o 40 o 39o 38oS

1

(b)

200

43 o 42 o 41 o 40 o 39o 38oS

Time (Year)

50

10

0

10

100

5
50

50

10
172o

20
o
174

176o

178o E

Figure 5.2: Results from the central New Zealand region.
(a) Plot of the conditional intensity function of time integrated over the studied region.
The horizontal axis represents time in years, and the vertical axis represents the ratio of theconditional intensity to the background intensity over the whole region; (b) Contours of
the estimated total rate m^ 1 (x; y ) (see Equation (5.8)); (c) the clustering intensity ^(x; y )(see Equation (5.9)); (d) the estimated relative clustering e ect C^ (x; y ) (see Equation(5.16)); and (e) the estimated background seismicity rate ^(x; y ) (see Equation (5.10)).
In (b), (c) and (e), the values are in events=deg 2 throughout the whole period (T =10500 days), and contours are drawn for approximately equal intervals in a logarithmic
scale. The thick solid lines represent the coast lines.

5.5.

EXAMPLES65

5.2(e). It can be seen that the background rate is much smoother than the total rate in
Figure 5.2(b). Most places have values in the range 2050 events=deg 2 except for theplate boundary, where the high values form up a strip with values 50100 events=deg 2
and some peaks up to 200 events=deg 2 . Also there are some high values in the volcanic
zones and o shore from Wanganui.

5.5.2 Central and western JapanThe second example used the data from the central and western Honshu area of Japan
from the hypocenter catalogue compiled by the Japan Meteorological Agency (JMA)(see Figure 5.3(a)). We selected data for the period 19261995 in the rectangular area34Æ 39Æ N and 131Æ  140Æ E with magnitudes MJ  4.0 and depths 100 km. Outside
the southern boundary there were two great earthquakes, namely the 1944 Tonankai
earthquake of MS 8.0 and 1946 Nankai earthquake of MS 8.2 with epicenters (136.62ÆE,33.80Æ N) and (135.62ÆE, 33.03Æ N), respectively. These events and their large aftershocks
apparently a ected the seismic activity in the above selected inland area. However, thedetection rate of earthquakes o shore is quite low compare to that of inland, and therefore
it is hard to widen the region of the analysis without raising the threshold magnitude
substantially. Instead, we keep the same area and same magnitude threshold as above forthe analysis, but allow for boundary e ects arising from the occurrence of the earthquakes
not more than 1 degree outside the study area, in order to include the in uence of thetwo great events with magnitude  7:0. Thus the likelihood calculation is implemented
for the enlarged history Ht containing data from the extended region.
The MLEs for this region are A^=0:204 event=(deg 2  day ), ^ =1:14, c^ = 0:0404 day ,p^=1:250, d^=0:00103=deg 2. Figures 5.4 (a) to (e) show the estimates of the time variationof the function p(t) in (5.15), the total intensity, the relative clustering e ect and thebackground intensity, respectively. We can see that the clustering e ects arising from

66

CHAPTER 5.

STOCHASTIC DECLUSTERING

typical aftershocks are eliminated in the estimates of the background intensity Figure
5.4(d). Those include the 1927 Kita-Tango, 1943 Tottori, 1945 Mikawa, 1948 Fukui,1964 Niigata, 1993 Noto-Hanto-Oki, and 1995 Kobe earthquakes. On the other hand,the background rates are almost the same as the total intensities in the Wakayama and

(a)

132o

134o

136o

138o

140oE

134o

136o

138o

140oE

(b)

’30

’45

’60

’75

’90

34o

35o

36o

37o

38o

39oN

Ibaragi regions which show almost pure spatial occurrence of non-clustered events.

132o

Figure 5.3: Seismicity the western and central Honshu area, Japan.
(a) Epicentral locations and (b) space-time plots for the longitudes against the times of
all the events.

5.6 Probability-based declustering algorithm
With the estimated background intensity and the branching structure (the space-time
ETAS model), we can carry out the following declustering algorithm, which is called the

5.6.

67

PROBABILITY-BASED DECLUSTERING ALGORITHM

Relative intensity
1
10 100

(a)

’25

’30

’35

’40

’45

’50

’55

’60

’65

’70

’75

’80

’85

’90 ’95’00

Time (Year)

0

(c)

38

50

50
20

10

0
10

34
133

134

135

136

137

138

(d)

8

0.

0.4

37
0.2

0.2

0.4

0.6

(e)

20

38

50

1
0.2

36
5

0.2
0.2

0.4

0.4

0.6

0.6

20

37

0.4

0.8

0.4 0.
2

0.4

50

34
131 132 133 134 135 136 137 138 139 140

10

20

100

50

35

50
20

0

20

35

140

39

38

36

139

6

39

132

5

10 5
20
34
131 132 133 134 135 136 137 138 139 140

0.

131

50

50

100

50

100

35

5

20

200

0

35

50

50

50

100

10
5

0
50

100

100

36

50

10

10

5

1
10

0

20

20

200

36

200

10
1

20

37

10

37

5

10

10

(b)
38

10

39

100

39

100

50
10

34
131 132 133 134 135 136 137 138 139 140

Figure 5.4: Results from the central and western Honshu area of Japan.
(a) Plot of the conditional intensity function of time integrated over the studied region.
The horizontal axis represents time in years, and the vertical axis represents the ratio of theconditional intensity to the background intensity over the whole region; (b) Contours of
the estimated total rate m^ 1 (x; y ) (see Equation (5.8)); (c) the clustering intensity ^(x; y )(see Equation (5.9)); (d) the estimated relative clustering e ect C^ (x; y ) (see Equation(5.16)); and (e) the estimated background seismicity rate ^(x; y ) (see Equation (5.10)).
In (b), (c) and (e), the values are in events=deg 2 throughout the whole period (T =26750 days), and contours are drawn for approximately equal intervals in a logarithmic
scale. The thick solid lines represent the coast lines.

68

CHAPTER 5.

STOCHASTIC DECLUSTERING

stochastic declustering algorithm.

Algorithm 2
1. For each earthquake j = 1; 2;    ; N , calculate probability j in (5.2) from the nal
solution in Algorithm 1;
2. Generate N uniform random numbers U1 ; U2 ;    ; UN in [0; 1].
3. If Uj < 1

j , keep the j -th event; otherwise, delete it from the catalogue as ano spring. The remaining events can be regarded as the background events.
The essence of the algorithm is the probability j associated with each event. Figure
5.5 shows the histogram of those probabilities of the two data sets. It can be seen thatwith very high probabilities most of events are either independent events or aftershocks.
Nevertheless, quite a few events are in between the two sides: about 31.7% of the New
Zealand events and 18.3% of the Japanese events have probabilities from 0.1 to 0.9 of

(b)
1000

(a)

0

100

500

200

300

400

500

1500

being an aftershock.

0.2

0.4

ρ

0.6

0.8

0.0

0.2

0.4

ρ

0.6

0.8

1.0

Figure 5.5: Histograms of aftershock probabilities (c.f (5.2)) for (a) the central New
Zealand region; and (b) the western and central Honshu region.
Using Algorithm 2, we have obtained examples of catalogues of background eventsfrom the two data sets. Typical examples are shown in Figures 5.6 and 5.7. Figures

5.6.

PROBABILITY-BASED DECLUSTERING ALGORITHM

69

5.6(a){(b) and 5.7(a){(b) show the epicentral locations of the independent events and theaftershock events, respectively; Figures 5.6(c){(d) and 5.7(c){(d) include the space-timeplots of the independent events and the aftershock events, respectively. Here the latitudeis used for the New Zealand data while the longitude is used for the Japanese data.
It should be noted that, in Figure 5.7(c) for the declustered events, two clusters of
linearly aligned events appear in 1945 and 1976 and in longitude ranges from 136ÆE to138Æ E and from 137ÆE to 140Æ E, respectively (see the arrows in the gure). The rst
cluster followed the December 1944 Tonankai earthquake of MJ 7.9 and preceded the January 1945 Mikawa earthquake of MJ 6.8 that occurred at 137.07Æ E and 34.68Æ N, We havelearned that these events include quite a few events remotely triggered by the Tonankai
earthquake whose epicenter (136.18Æ E, 33.57ÆN) is outside the studied area. In otherwords, the remotely triggered clusters in time have not been removed although the log
likelihood takes into account the e ect of the boundary condition of the Tonankai event.
The magnitudes of all the events in the second cluster are relatively small, and a possible triggering event is the 1975 deep great event of MJ 7.7 which occurred at (130.09ÆE,
38.79Æ N) at the depth of 549 km. The main reason for the detection of the remote triggering here is probably owing to the short-range range decay of the spatial response function(5.14) that conforms to the conventional understanding of aftershock area (Utsu, 1969).
Figures 5.6(d) and 5.7(d) include quite a few single events from which one may feelthat the background patterns have not been well separated from the clusters. However,the single or double events which remain in the cluster catalogue are the remains of thenumerous double or triple events in the original catalogue, most of which represent small
clusters from mainshock with magnitudes close to the lower magnitude threshold.
The next algorithm provides a stochastic realization of family trees based on theestimated space-time ETAS model.

Algorithm 3

70

38oS

(a)

41o
42o
178oE

43o

176o

’70

38oS

174o

’75

’80

’85

’90

’95

’80

’85

’90

’95

(d)

40o
41o
172o

174o

176o

178oE

43o

42o

41o

40o

39o

(b)

42o
43o

(c)

40o

40o
41o
42o
43o

172o

39o

38oS

170o

170o

STOCHASTIC DECLUSTERING

39o

39o

38oS

CHAPTER 5.

’70

’75

Figure 5.6: Results from applying Algorithm 2 to the central New Zealand region.
(a) epicenter locations of the background events; (b) epicenter locations of the o spring
events; (c) space-time plots of latitudes against times of the background events; and (d)space-time plots of latitudes against times of the o spring events.

39oN
37o

132o

134o

136o

138o

140oE

132o

134o

136o

138o

140oE

134o

136o

138o

140oE

(d)

’60
’45
’30

’30

’45

’60

’75

’90

(c)

34o

35o

36o

37o
36o
35o
34o

(b)

38o

(a)

’75

’90

71

PROBABILITY-BASED DECLUSTERING ALGORITHM

38o

39oN

5.6.

132o

134o

136o

138o

140oE

132o

Figure 5.7: Results from applying Algorithm 2 to the western and central Honshu area.
(a) epicenter locations of the background events; (b) epicenter locations of the o spring
events; (c) and (d) are space-time plots for the longitudes against the times of the background events and the o spring events, respectively. The two arrows mark the occurrencetime of the remote triggered events after 1944 and 1975.

72

CHAPTER 5.

STOCHASTIC DECLUSTERING

1. For each pair of earthquakes i; j = 1; 2;    ; N (i < j ), calculate probability i;j in(5.1) and 'j in (5.3) from the nal solution in Algorithm 1.
2. set j = 1.
3. Generate a uniform random number Uj in [0; 1].
4. If Uj < 'j , then the j -th event is regarded to be a immigrant (background event).
5. Otherwise, select the smallest I such that Uj < 'j +PI
i=1 i;j . Then the j -th eventis regarded to be a descendant of I -th event.
6. If j = N , terminate the algorithm; else set j = j + 1 and go to Step 3.
This algorithm enables us to identify a set of cluster families each starting from animmigrant. The earthquake dataset is stochastically separated into clusters and setsconsisting of single events; the separation is stochastic because it depends on the series of
uniform random numbers used in the above algorithm. It should be noted that the events
selected for the background here (and consequently those generated by Algorithm 2) arealways the rst occurring events in the cluster family, but not always the largest events inthe family. However, if preferred, one can take the set of the largest event in each cluster
as an alternative representation for the background seismicity.

5.7 Results for the North China region
The studied data are taken from the catalogue compiled by the Center for Analysis andPrediction of China Seismological Bureau. The considered region is 100Æ  125ÆE , 30Æ 
43Æ N . The magnitude threshold is ML 3:5. The data start from January 1st 1970 andends at December 31st 1999. Because all the earthquakes in this study are continental
ones, events of all the depths are used.

5.7.

RESULTS FOR THE NORTH CHINA REGION73

The seismicity in North China mainly distributed along several great seismic belts, as
shown in Figure 5.8. The main seismic belts are, from east to west, the Tanlu belt, theNorth China plain belt, the Shanxi belt and the North-South belt, of north-north-east
direction. According to these seismic belts, the selected seismic region can be divided into
four subregions: (a) the eastern subregion, (b) the Erdos subregion in the middle, (c) the
western subregion and (d) the southern subregion. The eastern subregion is separated
from the Erdos subregion by the Taihang Mountain Ranges and the west subregion by theboundary between the Erdos plateau and the North-South seismic belts. The southernsubregion is on the south of the Erdos subregion and the seashore subregion, separated
by approximately by the 34Æ N latitude, and is on the east of the west region. Figure
5.8(a) shows epicenters of the selected earthquakes, (b) and (c) shows their space-time
distribution.
The estimated parameters are A = 0:546 events/day, c = 0:00271 day, = 1:15,
p = 1:03, d = 0:000288 degree2 . Parameters for the variable kernel estimates are np = 5
and  = 0:05 degree. The estimates of the total space intensity, the background intensity,
the clustering intensity and the clustering coeÆcients are plotted in Figures 5.9 and 5.10.
From these gures, we nd the following features.

1. The western subregion has the most active seismicity in both total intensity andbackground intensity.
2. The eastern subregion is the second most active region, but seismicity is most clustered.
3. The Erdos subregion and the southern subregion are relatively aseismic. Only small
clusters occur in the southern subregion and at the boundaries of the Erdos subregion.

74

CHAPTER 5.

STOCHASTIC DECLUSTERING

4. In Figure 5.10(d), the areas of highest clustering intensities in the eastern subregion include the aftershock regions of the 1975/2/4 (MS 7:4) Haicheng earthquakeat (40Æ420 N, 122Æ410 E), the 1976/7/28 (MS 7:9) Tangshan earthquake at (39Æ 380 N,
118Æ 110E), and the 1998/1/10 (MS 6:2) Zhangbei earthquake at (41Æ 60 N, 114Æ 180 E);
and those in the western subregion include aftershock regions of the 1973/2/6(MS 7:7) Songpan earthquake at (31Æ 300 N, 100Æ240 E), the 1973/8/11 (MS 7:1) Luhuo
earthquake at (32Æ 550 N, 103Æ 550 E), the 1990/4/26 (MS 7:1) Gonghe earthquake at(36Æ 70 N, 100Æ 080 E) and the 1990/4/26 (MS 7:1) Gonghe earthquake at (36Æ 70 N,
100Æ 080E).
5. It is can be seen in Figure 5.9(b) that some seismic belts are not connected byhigh background seismicity regions. Such examples are the Fenwei belt in the east
subregion which is located on the east boundary of the Erdos region and the North
China plain belt in the middle of the east region. This is owing to the estimatesin the short time span of data. If we include the historical data before 1970 as thepast history, we may obtain connected high intensities in these belts. Meanwhile,
this may indicate the seismic gaps, where future large earthquakes will possibly llthem up.
Precisely saying, Figures 5.9 and 5.10 are the overview of the background and clustering seismicity in the only short time span during the years from 1970 to 1999. If the timespan of the data is long enough, then Figure 5.9(b) will provide the potential of future
large earthquake occurrences, and 5.10(c) shows how clustered events will be excited whenan event occurs in this region.
Applying Algorithm 2 to this data set, we get the a catalogue of the background events,
shown in Figure 5.11, and a catalogue of the clustering events, shown in Figure 5.12. From
the rst gure, the stochastic declustering procedure produces a quite temporally uniform

75

RESULTS FOR THE NORTH CHINA REGION

42

5.7.

(I)

40

(III)

36
34

Latitude

38

(II)

30

32

(IV)

100

105

110

115

120

125

36
30

32

34

Latitude

38

40

42

Longitude

0

2000

4000

6000

8000

10000

6000
0

2000

4000

Time

8000

10000

Time

100

105

110

115

120

125

Longitude

Figure 5.8: Seismicity and subregions in the north China region (1970{1999, M  3:5).
(a) Distribution of epicenters; (b) Latitudes against times; and (c) Times against longitudes. Subregions are (I) the eastern subregion, (II) the Erdos subregion, (III) the western
subregion and (IV) the southern subregion.

76

CHAPTER 5.

STOCHASTIC DECLUSTERING

1

104

1

10

10

42
1

10
100

10

100

40

1000
100

10

103

1

0.1
100 10

10

36

1
10

100
10

10

102

10

10
0

38

100

0
10

10
10

1

101

10

10

10
10

1

1

10

100

100

34

10

1
10
0
10

10

1

30
100

10

105

10
1

10

10

10
100

32

1

110

115

120

125

1

1

104

1

10

42

10

1
100

10

103

1

100

10

0.1

38

10

10

100

40

102

10

10

10

100 10

10

10

10

1

10

10

36

10−1

101

10

1
10

34

1

10

1

10

10

10

1

10

10

10

10
1

30
100

105

10

32

10

10−1

1

1

110

115

120

125

Figure 5.9: Estimated spatial intensity of the total process, the background process andthe clustering process.
(a) Total spatial intensity; (b) background intensity.

5.7.

77

RESULTS FOR THE NORTH CHINA REGION

104

0.1

42

1
10

1

1
10

1

10

10

10

0.1

1

36

10

1
0.1

1

0.1

10

101

1

1

1

1

1

10

10
0.1

10
10
0

0.1

10

0.1

1

1

0.1

10
0.1

1

0.1

1

10
0.1

1

10 1

0.1 1 10 0.1 10

1

10−1

1

1

30
100

102

10 10
1

1

10

10 10

10

1
0.1 1

0.1
10
10
0
10

10

10
0.1
10
1

10

34

0.1

1

1

10

1

1

0.

32

10

10

0.1

10 10
1

103

0.1

10
1

38

10
00
10 0
10

10

40

10
1001

1

10

10

0.1

10

10

105

110

115

120

125

1.0
42

0.2

0.2

0.4

0.2

0.4

0.2

0.4

0.2

8

0.4

30
100

0.2

0.4
0.2

0.2

0.4

0.2
0.2

0.4
0.2

0.6
0.2
0.2

0.4
0.6

0.2

0.4

0.2
0.2

0.2

0.2

0.2

0.6

0.8
0.4
0.2

2

0.2

0.2
0.2 0.2

0.2

0.6

0.4

0.

0.2
0.2

0.8

34

0.6

0.2

0.4

0.4 0.2
0.2

32

0.8

0.4

0.2

0.8

0.4

0.4

0.2

0.2

0.4

0.

36

0.2
0.2 0.4

0.2
0.4

0.20.4
0.6

0.6

0.2

2

0.6

0.2

0.4

0.

0.6
0.60.4

0.2

0.2

38

0.2 0.2

0.4

0.2

0.8

0.2
0.4

0.4
0.2

0.2

40

0.2

0.8
0.6

0.2 0.4

0.0
105

110

115

120

125

Figure 5.10: (continued of Figure 5.9) Estimated spatial intensity of the total process, thebackground process and the clustering process.
(a) Total spatial intensity; (b) background intensity; (c) clustering intensity; and (d)clustering coeÆcients.

78

CHAPTER 5.

STOCHASTIC DECLUSTERING

background catalogue, although, during the period from 1980 to 1984, the backgroundseismcity seems lower than other periods. The contrast between features of the east
subregion and the west subregion can be also views from these two gures: the east
subregion has a low background intensity than the west subregion, but with bigger cluster
sizes. The Erdos subregion is aseismic and most unclustered. In summary, the stochastic
declustering technique is shown to be useful in the study of seismicity in North China,
too, which helps us to nd di erent features in di erent subregions.

5.8 Discussion
In this section I would like to discuss various additional issues relating to the methods
proposed in the previous sections and also to the preliminary results obtained there.

(i) Justi cation of the underlying model
In the rst instance, the success of the proposed algorithms here depends on a successfulmodelling of the underlying clustering process. The ETAS model used has a long historyof use for this purpose, initially in the time domain (see. e.g. Ogata, 1988, 1999; Utsu
et al., 1995, for reviews of earlier work) and more recently in the space-time domain.
In particular, Ogata (1998) compared di erent possible forms of the spatial response
functions arising within the general framework of the space-time branching cluster modelusing Akaike's Information Criterion (AIC; Akaike, 1974).
More generally, in justifying a particular model on the basis of a declustering algorithm, it is worth bearing in mind that, in addition to such likelihood based tests, several
conventional graphical statistics are available for testing the absolute goodness-of- t of
the model to speci c space and time patterns of the data. These include the variancetime curve, the log-survivor function, etc. for the time sequences, (e.g., Ogata, 1988, and79

DISCUSSION

36
30

32

34

Latitude

38

40

42

5.8.

100

105

110

115

120

125

36
30

32

34

Latitude

38

40

42

Longitude

0

2000

4000

6000

8000

10000

6000
0

2000

4000

Time

8000

10000

Time

100

105

110

115

120

125

Longitude

Figure 5.11: Space-time distributions of the background events that are stochasticallydeclustered from the North China data set.
(a) Epicenter distribution; (b) plot of location latitudes against occurrence times; and(c) plot of occurrence times against location longitudes. Sizes of circles indicate themagnitudes of earthquakes from M 3:5 to M 7:9.

80

STOCHASTIC DECLUSTERING

36
30

32

34

Latitude

38

40

42

CHAPTER 5.

100

105

110

115

120

125

36
30

32

34

Latitude

38

40

42

Longitude

0

2000

4000

6000

8000

10000

6000
0

2000

4000

Time

8000

10000

Time

100

105

110

115

120

125

Longitude

Figure 5.12: Space-time distributions of the triggered events that are complements ofstochastically declustered events from the North China data set.
(a) Epicenter distribution; (b) a plot of location latitudes against occurrence times; and(c) a plot of occurrence times against location longitudes. Size of circle indicate themagnitudes of earthquakes from M 3:5 to M 7:9.

5.8.

DISCUSSION

81

Andersen et al., 1993) and K -statistic, nearest-neighbor statistics, etc. for spatial point
patterns (e.g., Ripley, 1981, and Diggle, 1983). Graphical comparisons of results from theoriginal data with those from repeated simulations from the model are particularly useful;see Ogata (1998) for details of simulation algorithms for space-time data. More sensitivegraphical comparisons can be carried out by using similar statistics for the transformed
times (e.g., Ogata, 1988) and possibly for the transformed higher dimensional data including locations (e.g., Schoenberg, 1999) in order to test against the stationary Poissonprocess.

(ii) Dealing with boundary e ects
As observed in Section 6.2, boundary e ects cannot be neglected. Taking these into
account, even if only in a rough way, can substantially improve the estimate of the background activity. In the present analysis, only the large events outside the region had
a signi cant e ect on estimates within the study region, because of the relatively shortrange of the spatial clustering terms. Thus it was suÆcient to include only the coordinatesof the larger events from outside the region in the augmented history. This represented
quite a useful compromise, insofar as it reduced both the the computational work and thedemands on the underlying catalogue, without causing signi cant loss of accuracy.

(iii) Interpretation of temporal uctuations in the estimated background activity
The term (x; y ) for the background seismicity implies that the space-time characteristicsof the declustered process are close to those of a Poisson process that is stationary in timeeven if it is non-homogeneous in space. In particular, the proposed algorithm assumes aconstant mean rate in time for the background seismicity at each location. If, then, thethinned outcomes for the background events do not appear uniform in time, this should

82

CHAPTER 5.

STOCHASTIC DECLUSTERING

indicate some departure from the base structure assumed for the clustering e ects. For
example, it may be caused by temporal activation or quiescence which has some physicalcause outside those which create the standard ETAS clusters, or it may be due to somehitherto undetected fault in the catalogue data. Revealing such phenomena and evaluating
their signi cance, or uncertainty, are the eventual aims of the present declustering method.
In e ect, the model used for declustering plays the role of a null model, and departuresfrom it represent e ects worthy of further study and interpretation.
For example, in our study, the aligned declustered events in space-time plots in Figure
5.7(c) that were discussed previously clearly indicated the presence of remotely triggered
events. This in turn suggests that the short-range decay function (5.14) adopted as thespatial response function of the model, may require modi cation, perhaps replacement by
a power law decay function, as was already indicated in Ogata (1998).
It is also possible that in some smaller regions, apparent di erences from stationarity
in time may be caused by artifacts of the data set, in particular, the substantial numbers of smaller aftershocks which may be missed (particularly in the early part of thecatalogue) in the time immediately after a mainshock. The rate of missed events aroundthe threshold magnitude decreases as time elapses after the mainshock; and, further, thistype of incompleteness of the data has been improving, as small earthquakes can now be
detected more rapidly. Consequently, the estimates of c and p in the modi ed Omori function in (5.14) for the aftershock sequences may be biased by the missing events in earlyyears, giving rise to some under-declustering of the later aftershocks. Some indications of
this e ect can be seen in the space-time plot of the declustered events in Figure 5.7(e) forthe earthquakes during the period 19431946 and with the longitude range 135Æ 138Æ E.

5.8.

DISCUSSION

83(iv) Non-uniqueness of the declustered catalogue
A stochastic declustered catalog is not unique, as it is dependent on the random numbers used in selecting the events to form the background seismicity. Unlike conventional
declustering methods, the stochastic declustering method does not make a xed judgment
on whether an event is an aftershock or not. Instead, it gives the probability of how likely
each event is to be an aftershock. Each catalogue produced by the stochastic declustering
method should be understood as a simulation, or re-sampling from the original catalogue.
This may be considered an advantage of the method rather than the reverse, as it allows uncertainty about the declustering to be quanti ed. Thus, by repeating the random
declustering, we can easily produce many stochastic copies of the declustered catalogue.
From these, graphical or other statistics can be calculated to discuss the uncertainty or
signi cance of a particular feature of the declustered catalogue.

(v) Preservation of information
One may argue that throwing aftershocks out by the thinning (random deletion) method
loses potentially useful information for carrying out seismicity analysis. While this may
be true of conventional declustering methods, or for one particular realization of thestochastically declustered catalogue, the essential output of the stochastic declustering
algorithm is not the particular realization of the declustered catalogue, but the probabilities associated with the decomposition of the intensity into background and clusteringterms. No data is lost in this process; rather, it adds a new item, the probability j , tothe parameters of the events in the original catalogue. It should not be forgotten thatthe basic purpose of declustering is not the production of declustered catalogues for theirown sake, but the investigation of features of the data which otherwise may be obscuredby the clustering e ects. For this purpose, the best representation of the combined information in both data and model is the ensemble of these probabilities for the original

84

CHAPTER 5.

STOCHASTIC DECLUSTERING

events, whether studied directly, or through the stochastically declustered catalogues towhich they give rise.

(vi) Avoiding the selection of mainshocks as aftershocks
As we can see in the declustering algorithm, the magnitude of an earthquake does notplay a role in the random deletion. If a mainshock is preceded shortly by one or moreforeshocks, it has some probability of being interpreted as an o spring of one of theforeshocks, and in this case may not be selected as a background event. If this is regarded
as a disadvantage, we can separate the clusters by Algorithm 3, and then adopt the largest
event in each cluster, instead of the initial event, as the representative for inclusion in the
catalogue of background events.

(vii) Prospects for future study
Possibilities for further study of the stochastic declustering method can be discussed under
two headings: theoretical problems and applications. On the theoretical side, there arestill important aspects of the model to be considered. For example, we did not consider
the important problem of deciding when the use of a non-homogemeous background isnecessary, either in space or in time, nor of how to estimate such time varying e ects
in the background if indeed these are present. On the applications side, the stochastic
declustering algorithm gives us a useful tool for the study of seismicity patterns. For
examples, in this paper, we have found some interesting di erences between the New
Zealand catalogue and and the Japanese catalogue: (i) the Japanese catalogue has denser
clusters in space than the New Zealand catalogue as can be seen from the parameters
estimated by MLE; (ii) the New Zealand catalogue has more events which cannot beclearly classi ed (31.7% events with probabilities 0.10.9 of being aftershocks) than theJapanese catalogue (only 18.3% events with probabilities 0.1 90 of being aftershocks).

5.9.

CONCLUSION OF CHAPTER 5

85

The reasons for such di erences need further study. Moreover, stochastic declustering
can also help to answer questions such as: does relative quiescence appear before a large
earthquake? does the occurrence of a large earthquake change the seismicity patterns in
the neighboring region? what is the probability for an event to produce a larger o spring,i.e., the probability of the occurrence of a foreshock? All of these important problems
require a careful evaluation of uncertainties, and while they all need further study, thestochastic declustering algorithm suggests one way that such problems might be tackled
in a more objective fashion than has been possible in the past.

5.9 Conclusion of Chapter 5
In this chapter we have outlined a procedure for objectively estimating the spatial intensityof the declustered or background seismicity from an earthquake catalogue. The method
is based on an assumed stochastic model (space-time ETAS model) for the clustering
patterns. Its core is a 2-stage iterative algorithm based on the thinning method in Section
5.2 and the estimation of the background seismicity in (5.8). The procedure stochasticallysplits the whole process into background events and o spring events grouped into clusters
or family trees. Each such splitting can be regarded as on possible realization of how theoriginal catalogue could have been produced, taking into account that individual events
can be classi ed as background events or cluster events only with a certain probability.
We can estimate the signi cance and uncertainties of the features of the catalogue notexplainable in terms of the assumed model for clustering by running the procedure manytimes.
Using this method, preliminary studies were made of earthquake data sets from thecentral New Zealand region, and the from the central and western region of Japan. Although preliminary, the results revealed several features of interests which could be thesubject of future study and interpretation.

86

CHAPTER 5.

STOCHASTIC DECLUSTERING

Part III

MODELLING EARTHQUAKESAND PRECURSORS

87

Chapter 6

Formulation of Models on
Earthquake and Precursors

6.1 On precursory observation data
Compared with point process models existing for seismicity only, there are few models
for examining the correlation between suspected precursors and earthquakes. The reason
is not only because we know very little about the mechanism of the earthquake processitself, but also because even less knowledge about the mechanisms how the precursors
are produced. Besides these, there are many other technical diÆculties in modelling therelationship between earthquake occurrences and the precursory observations. One ofthem is that the data formats of precursor observations may be quite di erent from thatof point processes. These formats include:
1. Continuous processes. At most stations, some precursory data are recorded continuously in time, such as tilt, geomagnetism, etc..
2. Point processes. Irregularly observed records describing times and various types ofanomalous events, such as animal anomaly behaviors.
3. Time series. Regularly observed records during a xed time interval, such as groundelectric resistances.
89

90CHAPTER 6.

FORMULATION OF MODELS ON EARTHQUAKE AND PRECURSORS

More unfortunately, even for these complicated data, there are also some arti cial problems, such as missing data, outliers, non-homogenous records due to changing monitoring
instruments and changing quanti cation methods etc. Despite these problems, the mostimportant strategy for statistical analysis is to nd a data set which has got a long enoughobservation.
To model a process of earthquake occurrences incorporated with the precursor data,
we introduce the conditional intensity function being dependent on the multi-histories ofobservations from several processes,

(tjHt(E ) ; Ht(P ) ) = E[N [t; t + dt)jHt(E ) ; Ht(P ) ];(6.1)

where Ht(E ) and Ht(P ) are the observed histories from an earthquake process and a precursor
process, respectively. To distinguish such a process from the previous processes whose
conditional intensity depends only on the history earthquake process itself, we named
it a partial process (Vere-Jones, 2002, personal communication), and the second case acomplete process.
As I described before, it is diÆcult to reveal the physical mechanism how the precursors
are produced in the earthquake preparation process. Nevertheless, we can still get somehints from the experiences in the practical earthquake predictions. After some precursory
phenomena are observed, people tend to increase their risk of earthquakes for the near
future. The length of the period with such a high risk depends on the experience of theexperts in earthquake prediction, as the same way of practical earthquake prediction in
China, and the more the precursors are observed, the higher the earthquake risk is.
This prediction method leads us to make use of the mutal-excitation process (Hawkes,
1971; Ogata et al., 1982). In such a process, there are two kinds of stochastic processes,
the target point process and external processes that may be further point processes,
continuous time stochastic processes or time series, and the occurrence rate of the target

6.2.

91

THE LINLIN MODEL

point process is at least partly controlled by the external processes. I would like to discussthese model and their extensions in the following sections.

6.2 The LINLIN model
The linear intensity model was rst introduced by Hawkes (1971), named as the self- andmutally-exciting model. Then, Ogata et al. (1982) used the latter to investigate the correlation between earthquake occurrences in di erent regions. In this model, the risk of theevents consists of the following three components. The rst component is called the background or the trend term. The second component is called the self-excitation component,
i.e, each event in the process excites further events, and also the newly occurred events
excites their own o spring. The third component, named the external-excitation term,
describes the contribution from an external process. In the Hawkes' model, the third
term is considered as a point process only (Hawkes, 1971). Ogata and Akaike (1982) andOgata et al. (1982) generalized it to be a stochastic process.
The model is de ned by the conditional intensity function such that

(tjH(E ) ; H(P ) ) =  +
t

Z t

t

1

g (t s)N (ds) +Z t

1h(t s)dXs ;(6.2)

where Xt is either a point process or a cumulative process

Xt =

Z t

1

x(s)ds(6.3)

of a stochastic process x(t). The response functions g (t) and h(t) are parameterized by
the Laguerre type polynomialsg (t) =

KX
k=1

ak tk 1 etandh(t) =

KX
k=1

bk tk 1 et:

Ogata et al. (1982) used this to examine the correlation between the intermediate depth
earthquakes (M  5:0) in the Hida region and the shallow earthquakes (M  5:5) in92CHAPTER 6.

FORMULATION OF MODELS ON EARTHQUAKE AND PRECURSORS

the central Kwanto region, Japan. Ma and Vere-Jones (2000) also used this model toinvestigate the correlation between shallow and deep earthquakes in the New Zealandregion.

6.3 Extensions of the LINLIN model
In the LINLIN model, the earthquake risk is divided into two part, one part being due toearthquake process itself, the other part can be interpreted by the external process. That
is to say, the conditional intensity function of the LINLIN model is the one among thegeneral form of linear models

(tjHt(E ) ; Ht(P ) ) = 0 (tjHt(E ) ) + 1 (tjHt(P ) );(6.4)

(tjHt(E ) ; Ht(P ) ) = 0 (tjHt(E ) )1 (tjHt(P ) );(6.5)

or log-linear models
where 0 (tjHt(E ) ) is the part that is conditional on the history of the earthquake process,which can be interpreted by the earthquake process itself, and 1 (tjHt(P ) ) is the part thatcan be interpreted by the precursor.
Similar to the LINLIN model, we can specify

1 (tjHt(P ) ) =Z t

1

g (t u)(dXu)(6.6)

for the model formulation, where Xu is the precursor process, in the form of either a point
process or a stochastic process, g (t) is the response function, and  makes a transformationof the marks.

6.4 On precursory observations at discrete times
If precursor records are sampled at discrete times, we can still deal with the precursor
process as a step function of times with values jumping at the observation times. However,6.4.

ON PRECURSORY OBSERVATIONS AT DISCRETE TIMES

93

this method is a sort of convenience. Instead, a more natural way to deal with the discrete
time observations is to convert the earthquake process into a discrete one.
Let ni be the number of earthquake events occurred in an equally divided interval1[i; (i + 1)); i = 0; 1; 2;    ; N , then fni ; i = 0; 1;    ; N g is a non-negative integer
series. Similar to the condition intensity, we can treat ni as a Poisson random variablewith mean i that depends on the previous observations such that

i = E[ni j nk ; Xk ; k = 1;    ; i 1]:(6.7)

As the discrete-time version of the models (6.4) and (6.5), we can respectively consider
a linear modelor a log-linear model

i 1
(
E) X
i = i + Ak (Xi k );
k=0

(6.8)

i 1
X
(
E)
i = i expf Ak (Xi k )g;
k=0

(6.9)

where (iE ) are the contribution from history of earthquake occurrences and Ak are themselves free parameters or a function of k represented by some free parameters. The above
two models indicates that the contribution of the precursor observations are through a
convolution of the sequence fAk : k = 1; 2;   g and the transformed values of the
precursory records fXk g by a non-linear function .
Parameterization of both functions fAk g and  can be carried out as described in
the forthcoming section to estimate by using the maximum likelihood estimates. The
likelihood of the model in (6.7) is given in the formN  ni
Y
i i
L=
e
:
n
!
i
i=1

(6.10)

Here we suppose that the precursor observations are given at times ; 2; ; N , in such a waythat the observation time of the precursor at the same time interval is not earlier than the observation
time of n .
1



i

94CHAPTER 6.

FORMULATION OF MODELS ON EARTHQUAKE AND PRECURSORS

Another method for prediction and modelling is to convert the earthquake occurrences
into a binomial series such that,

vi =



1; ni > 0;0; ni = 0;

(6.11)

where ni represents the number of earthquakes in the ith interval, with probabilities ofsuccess pi conditioned on the previous observations,pi = Prfvi = 1 j vk ; Xk ; k = 1; 2;    i 1g:(6.12)

Thus, our models are converted into a binomial regression problem. We can then t thefollowing models to these binary observations,
i 1
X
(
E)
 (pi ) =  (pi ) + Ak (Xi k );k=0

(6.13)

where  (p) can be one of the followings:

 (p) = p;

for a linear probability model;

 (p) = log p; (p) = log
 (p) =



p



1 p
log(1 p);for a log-linear;;for a logistic model;
for a complementary log-linear model; ;and p(iE ) is the probability conditional on the history of the earthquakes process.
The likelihood function for the binomial models is given by

L=

N
Y
i=1

pvi i (1 pi )1 vi ;(6.14)

where N is the total number of observations.
Vere-Jones (2002) used logistic regression to model the earthquake probabilities basedon the output from applying the M8 algorithm (Kellis-Borok and Kossobokov, 1986) tothe earthquake catalogue from the New Zealand region.

6.5. ON THE MODELS FOR MORE MULTI-TYPE PRECURSORS AND MULTI-STATION OBSE

6.5 On the models for more multi-type precursors
and multi-station observations
Besides seismicity, we may observe several physical variables at one station, or observed
the same variable at several stations. The more observations are obtained, the morereliable could be the prediction. Of course, the model will be more complicated than theone only for only a single station.
Before modelling for the data from multi-type precursors, exploratory analysis of eachtype of precursors from a single station is of importance. In revealing the underlying patterns in addition to the relationship between the earthquakes and each type of precursors
at a single station, we can combine di erent models used in the exploratory analysis. In
this stage, the AIC procedure can play an important role in nding the best models.
Another important aspect of multi-observations is to develop space-time or synoptic
models embodying physical or statistical insights concerning the relation between thelocal seismicity and neighbouring precursor observations. As Vere-Jones (2002) dealt
with outputs from applying M8 algorithm to the New Zealand earthquake data, we candivide the whole target region and time interval into small windows of equal sizes, possibly
overlapped, centered at grid points. Then the regression can be carried out over thesewindows.

96CHAPTER 6.

FORMULATION OF MODELS ON EARTHQUAKE AND PRECURSORSChapter 7

Statistical Analysis on the ULF
Electric Signals

7.1 Data and model speci cation
The low-frequency electrical signals referred to in the chapter title have been used as apart of the reports to regular earthquake-monitoring meetings of the China Seismological
Bureau for these several decades. The instruments for measuring the signals were installedin many provinces in China. The data for the present study is taken from the importantregion around Beijing. Although the methods are well-accepted by Chinese seismologists,
there have been few careful statistical analyses of the signals and their explanatory power.
The main aim of the present chapter is to supply such a study. A preliminary version ofsome of the results was presented in a report by Ogata and Zhuang (2001).
The signals used in the study are based on the uctuations in the low-frequency
component of the di erence in potential between two electrodes buried some 6-8 metres
below the ground surface, and approximately 40 meters apart. The electrodes are placedsuÆciently deeply below the surface and in suÆciently quiet locations to avoid mostambient e ects associated with local rainfall, electrical storms, and domestic or industrial
sources. On the other hand they appear to be excited by stress or strain changes preceding
even moderate-sized earthquakes at quite considerable distances from the station. A97

98

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

detailed description of the apparatus, the monitoring and measurement procedures used,
and the resulting data, can be found in Appendix A. A signal strength is an average of
two terms, the rst based on the mean deviation of the amplitudes recorded in durationof the signal, and the second based on the maximum of the amplitudes recorded in the
duration of the signal. The daily strength we used below is the sum of the strengths ofall signals during the day. Again we refer to Appendix A for details.
The data forming the basis for the present study are taken from four stations in thegeneral vicinity of Beijing, namely Huailai (H), Langfang (L), Qingxian (Q) and Sanhe(S). These basic data are held in the database of the Centre for Analysis and Predictionof the China Seismological Bureau in Beijing. The locations of these stations, as well as
the epicentres of events with M  4 over the period 1982/1/1 - 1998/1/30 are shown foran extended area around Beijing in Figure 7.1. The dashed circles shown on this map are
300 km radius circles centred at each of the recording stations. The electrical signal datafrom each station are shown in Figures 7.2, where the daily signal from a given station iseither represented as a zero value (if the signal strength falls below a minimum threshold)or as a vertical bar with height proportional to the fourth root of the signal strength,
which is close to the powerof the transformation (X ) = X in the model discussed

later. The observation time intervals are:Huailai: [1826, 5874] days, 1987/May/1 { 1998/Jan/31;
Langfang, Qingxian: [0, 5477] days, 1982/Jan/1 { 1996/Dec/31;
Sanhe: [0, 5874] days, 1982/Jan/1 { 1998/Jan/31.
These di erent time intervals are due to data availibility.
The earthquake locations and magnitudes are taken from the catalogue of eventsin the Beijing region, also prepared by the Centre for Analysis and Prediction of the

99

DATA AND MODEL SPECIFICATION

40

S

H
Beijing
L

Q

36

38

Latitude

42

7.1.

4

4.5

5

5.5

6

112

114

116

118

120

6.5

122

Longitude
Figure 7.1: Epicenters of earthquakes occurred in the Beijing region and its vicinity during
the period 1982.1.1 { 1998.1.31.
The solid large circle represents the 300 km distance from Beijing, the dashed circles
represent 300 km from individual stations. The shaded rectangle encloses the region
where seismicity has remained high since the great Tangshan earthquake. \H", \L", \Q"
and \S" represent the names of the stations Huailai, Langfang, Qingxian and Sanhe,respectively. The crosses mark the locations of stations. The sizes of the small circles
correspond to the magnitudes of the events, as shown in the inset.

100

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

China Seismological Bureau. Target events are those with magnitudes M  4 within a
speci ed circle or other region. Before commencing the analysis, the Beijing catalogue was
declustered according to the following simple rule, which was suÆcient to remove most of
the relatively small numbers of aftershocks/foreshocks from the events under study: an
event in the catalogue was removed if it occurred within a distance of 50 km and a timeinterval of 5 days of some larger event in the catalogue.
0

1000

2000

Time (Days)
3000

4000

5000

6000

6
4

0
1000
6000

Sanhe
6
4

0
1000
6000

Qingxian
6
4

0
1000
6000

Langfang
6
4

0
1000
6000

Huailai
1982

1984

1986

1988

1990
1992
Time (Years)

1994

1996

1998

Figure 7.2: Signals observed at each station and the occurrence times of the earthquakes
within 300 km from the stations.
Vertical lines above the horizontal line in each panel give the magnitudes of earthquakes
occurring in a 300 km region about the station (see scale at the left hand side); heights
below the horizontal line correpond to the signal strengths on a 4th-root scale (see scaleat the right hand side).
The clusters of epicentres which appear on the maps are caused mainly by events
nearby in space accumulating over the period of the study, and are not generally due toclustering in time. The most prominent feature is the high level of activity in the region
around the great Tangshan earthquake to the South-East of Beijing, which persisted

7.1.

101

DATA AND MODEL SPECIFICATION

throughout most if not all of the observation period. Another important cluster includes
the pair of Datong events (1989/Oct/19 and 1991/March/26), both with magnitude 6.1,
to the South-West of Beijing. A magnitude-time plot of the events in the study region isshown in Fig 7.3, alongside a cumulative frequency plot and a b-value plot. It is very clear
from the cumulative frequency plot that the seismicity in the region decreased markedly
after the Datong earthquakes, a feature which complicates model tting and testing. The
approach adopted in this study is to rst t the data with models which assume timeinvariant processes, and then examine the e ects of allowing for separate rates in the
periods before and after the Datong events.
The basic statistical model used in this study is a discrete-time version of the Lin-Lin
model described in Utsu and Ogata (1994), and itself based on the multivariate Hawkestype point process model (Hawkes, 1971) rst discussed in a seismological context by
Ogata and Akaike (1982) and Ogata et al. (1982). Other studies based on this or similar
models include Musmeci and Vere-Jones (1986) and Nishizawa, Lei and Nagao (1994).
In the discrete-time version proposed here, we treat the number of events occurring ina 1-day period, say nt for day t, as a Poisson random variable with mean (t) dependent
on the signals detected over the previous days. More speci cally, we suppose that (t) isrelated to signals Xt i at preceding times t

(t) =  +

K
X
m=1

cm

i (i = 1; 2; :::) by the equation

"t 1
X
`=1

#

Am;` m (Xt ` ) :(7.1)

Here the rst sum is taken over the number of participating stations (which may varyfrom 1 to 4 in our discussions), and the second sum over the number of time lags. and cm are constants to be tted, the rst representing the level of background activity(unpredictable events), and the second the weight to be given to the signal from the mthstation. For each m, the Am;` ; ` = 1; 2; : : : ; constitute the response function for stationm, which will commonly be taken as a geometric decay, normalised to form a probability

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

4.0

4.5

Magnitude

5.0

5.5

6.0

1982 1984 1986 1988 1990 1992 1994 1996 1998
Time

CHAPTER 7.

80
60
40
20
0

6.0
4.5
4.0
50
5 10 20
2
1Figure 7.3: Plots of cumulative numbers (top panel) and magnitudes (second panel) of
the earthquakes that occurred within 300 km from Beijing against occurrence times.
The third panel shows a frequency-magnitude plot for the same events.

102

Cumulative number

Magntiude
5.0 5.5

N(>=M)

7.1.

103

DATA AND MODEL SPECIFICATION

distribution Am;` = (1m )`m . The quantities

m (X ) appearing in the right side ofthe equation are transformed versions of the electrical signals. The original signals varyenormously in strength, and the transformation was introduced to dampen the size of thisvariation. Various possible transformations were tried; the one found to be most useful
was the power transformationm (X ) = X

where

m

m is a positive quantity, which in general may depend on the station : typicalvalues are around= 1=4: The quantities appearing in these equations are constrained

so that the means (t) are positive for all t .
Equation (7.1) is taken as de ning the successive means for a sequence of Poisson
random variables nt , which are conditionally independent given the values of the previous
signals. In practice, since the 1-day observation interval is very small compared to theaverage time interval between events, the nt nearly always take either the value 1, if an
event occurs within the target region and over the target magnitude threshold in day t,
or 0. Remember the declustering of events which we have made, the likelihood can be
written in the form

L=

T 
Y
(t)ntt=1

nt !

e (t)



;

(7.2)

where T is the length of the observation period in days. If no multiple events (morethan one event within a single day) are observed during the observation period, then thefactorial terms can be omitted, since they are equal to unity.
Estimates of the parameters can then be made nding the values of the parameters(maximum likelihood estimates) which maximize the expression in (7.2) or its logarithm.
Moreover, di erent models can be compared by comparing the values of the maximised
likelihoods, or more usefully by comparing the AIC (Akaike Information Criterion; Akaike,1974) values

AIC = 2 log L^ + 2K ;

(7.3)

104

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

where L^ is the maximum of the likelihood and K is the number of tted parameters.
In comparing models with di erent numbers of parameters, addition of the quantity 2K
roughly compensates for the additional exibility which the extra parameters confer upon
the models which possess them. The model with the lower AIC value is taken as giving
a better t to the data.
The following is a brief summary of the remaining sections of this Chapter. In Section
7.2, a preliminary examination of the e ectiveness of each station is given, based on theability of the station in question to forecast events in its associated circle. This is followed
by a discussion of the e ectiveness of the four stations jointly in providing forecasts forevents in a 300 km circle centred on Beijing itself (Section 7.3). In Section 7.4, a special
study is then made of the events falling within the extended aftershock region of theTangshan earthquake, namely, the rectangular region in Fig 7.1. Section 7.5 examines
the performance of one-day probability forecasts based on the output from the Lin-Lin
model. Section 7.6 includes a discussion of some further issues, including in particular the
problems caused by the rate change following the Datong events, and some suggestions for
additional studies. Section 7.7 examines whether or not the electric signals are postseismic
phenomena. The nal section draws the conclusion.

7.2 Results for regional seismicity within the circles
around the individual stations
We rst examine the ability of the signal from each station to predict earthquakes in the
300-km radius circle centred on the station itself. The choice of radius here is a compromise
between extending the region too far, and including large numbers of distant events which
are not well-predicted, and making it large enough for the number of included events togive statistically signi cant results. The selected value of 300 km was suggested to us in
Appendix A. (see also the discussion in Chen et al, 1998), and is also in accordance with

7.2. RESULTS FOR REGIONAL SEISMICITY WITHIN THE CIRCLES AROUND THE INDIVID

Table 7.1: Summary of results for data within 200 km and 300 km circles around stations
Station
H
H1
L
Q
S
H
H1
L
Q
S

Time span[1826,5874][1826,5874][0,5477][0,5477][0,5874][1826,5874][1826,5874][0,5477][0,5477][0,5874]

N
53
53
92
67
82
22
22
56
32
57

N+
26
26
64
38
59
13
13
40
24
33


c

AIC
6.681e-3 3.6159e-14 0.9818
3.4797 -22.03
6.143e-3 1.007e-3
0.9042 0.3( xed) -16.93
7.030e-3 1.473e-2
0.49550.1529 -32.57
4.732e-3 2.109e-3
0.88840.1722 -24.46
6.507e-3 3.026e-3
0.78810.3072 -48.84
3.503e-3 3.164e-19
0.87385.0781 -17.70
3.012e-3 5.553e-3
0.8699 0.3( xed) -11.84
5.299e-3 4.175e-3 2.2202e-3 0.0940 -15.18
1.738e-3 9.064e-4
0.92320.1316 -12.38
4.813e-3 1.677e-3
0.7432
0.3922 -32.97

The rst 5 rows are the result from tting model to the 300 km circles around stations,and the last 5 rows are results for the 200 km circle around stations. From this tableon, H , L, Q and S represent the station names, Huailai, Langfang, Qingxian and Sanhe,respectively. H1 is for the model tting the signal data from the Huailai station but witha xed parameter = 0:3. Time is in days since January 1st 1982. N is the numberof events in the given space-time window. N+ is the number of events occurring when
(t) >  (predicted events). The parameters, (event=day ), c,  and are interpreted
in Section 7.1. AIC s are the di erence between the AIC values of Lin-Lin and theconstant-rate Poisson models for the dataset from each station, and the negative values
show the better t of the Lin-Lin model than the Poisson model.
our own preliminary experiments to nd the radius which gives the best discrimination
against the Poisson model. The results are summarized in Tables 7.1 below, rst for thesuggested radius of 300 km, and then for a smaller radius of 200 km.
In each case the basic model tted was the version of Equation (7.1) with just one termin the rst summation, namely the term corresponding to the selected station, and witha geometric form assumed for the response function. The likelihood, and also the AIC
value, was then compared with that for the corresponding stationary, discrete, Poissonmodel, obtained by retaining only the background term in Equation (7.1). Thus, from
Table 7.1 it can be seen that for the 300 km circle around each station, the likelihoodfor the Lin-Lin model is substantially higher than that for the Poisson model, and that
although the Lin-Lin model with geometric response function involves four parameters,
as against the Poisson model's one, the AIC value for the Lin-Lin model is also better.

106

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

A di erence in 2 for the AIC values is generally considered suÆcient to demonstrate one
model's superiority over another. Alternatively, the likelihood ratio test on 4

1=3

degrees of freedom suggests that a di erence of 7:815 2  3 = 1:815 should be suÆcientto establish the superiority of the Lin-Lin model at the 5% level. Although both of theseresults are approximate, they are repeated for all stations and for both radii. This makes
the combined results very compelling.
Despite the general signi cance of the results, some stations appear to perform rather
better than others. Huailai appears to give the weakest result, Sanhe the strongest. The
parameters for Huailai are rather unstable, and di erent from those for the other stations,a feature which may be associated with the change in background rate in its vicinity.
Huailai is also rather distant from the large group of events associated with the Tangshan
aftershock region. Langfang station has a very noisy record, possibly because it is situatedin an area that has rapidly become built-up, but again geological or topographical factors
may be relevant. Despite the di erences, all stations appear to respond well to some formof signal which frequently precedes intermediate size earthquakes, and is generally absent
on other occasions.
The response functions for the di erent stations are shown in Figure 7.4 for the 300
km circles. With the exception of Langfang, which has a much more rapid fall-o , thegeometric decay factors  are around 0.80.95, suggesting that the signals may precede
the earthquakes by periods of up to a week or more. Unfortunately there is not suÆcient
data to investigate the shape of the response function in detail; informal experiments withother forms of response, however, suggest that the geometric form is as e ective as any;some further discussion is given in Appendix A. The rapid decay in the case of Langfang
may be a further consequence of its noisy record, as the more rapid decay will tend todiminish the e ect of spurious signals.
The form of model used in Lin-Lin implies a rapid build-up in the risk if the signals

0.2

0.3

HuailaiLangfang
Qingxian
Sanhe

0.1

event/day

0.4

0.5

7.2. RESULTS FOR REGIONAL SEISMICITY WITHIN THE CIRCLES AROUND THE INDIVID

5

10

15

20

Days

Figure 7.4: Estimated response functions Am;` by tting the model to the 300 km circles
around stations.
are repeated for several days in a row; thus the overall response to a burst of electric
signals is generally not a simple geometric form, but rather an initial pulse which builds
up to a maximum and then decays geometrically. A visual impression of the response as
quanti ed by Equation (7.1), and its relations to the observed signals and the earthquakes
occurring within the observation region, can be gained, for example, from Figure 7.5. Hereand in the following Figure, the horizontal dotted line corresponds to the background rate

 . Events which occur when (t) >  are marked by circles (named predicted events);
those where (t) <  by crosses (named unpredicted events). The numbers of events inthe former category are denoted by N + in Tables 7.1.
To summarize, the statistical analysis shows that, for each station individually, thesingle-signal Lin-Lin model gives a substantially better t to the earthquake data in its
area than the Poisson model. In other words, each station picks up important predictive
information about the events in its area, while most if not all respond to some events
which are not picked up by the other stations.
A similar performance is observed whether the analysis is restricted to events within

0.08
0.02

1200

1400

1600

1800

2200

2400

2600

2800

3200

3400

3600

3800

4200

4400

4600

4800

0.02

0.08

0.02

0.08

800

0.08

event/day

600

0.02

event/day

400

0.08

event/day

200

0.02

event/dayevent/day

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

0.08

event/day

CHAPTER 7.

0.02

108

sh.s0−5874.ps

5200

5400

5600

5800

Figure 7.5: Plot of the conditional intensity function tted to the data from Sanhe (S)station.
The signal strengthes are shown on the top of each panel in logrithmic scale from 1 to
10000. The dashed horizontal line represents the mean rate. The earthquake are marked
by circles and crosses, with di erent sizes representing the di erent magnitudes.

7.3. RESULTS FOR THE SEISMICITY WITHIN 300 KM CIRCLE AROUND BEIJING

109

200 km of the station in addition to within 300 km of the station. Even with the moredistant events (300-400 km from the station) a signi cant proportion are preceded byobserved anomalies in the electric signals (see the graphs in Figure 7.8).

7.3 Results for the seismicity within 300 km circle
around Beijing
We now consider the combined explanatory power of the four sets of signals for events
within a 300 km circle around Beijing. For this purpose, a hierarchical sequence of models
was considered, rst with no signals (Poisson model), then with the signal from one station
only (four single Lin-Lin models) then with the signals from selected pairs of stations, andso on. Table 7.2 summarizes the results from all the combination of stations, i.e., the single
stations, pairs of stations, triplets of stations, and the set of all four stations together.
In general, as successive stations are added into the model, the performance improves,
i.e. the AIC values decrease. However, it can be seen that there is very little increasein the log-likelihood (and hence a large penalty e ect in the AIC values) in moving fromthe best pair of stations - namely Sanhe and Qingxian - to the full set of 4 stations. The
implication is that while each station individually carries important predictive information
about the target events, most of the useable information is incorporated into the response
from the two best-placed stations; once they are in place, the other stations have little
further to add. The implication is that Sanhe and Qingxian are particularly well-placed
to pick up the bulk of the information reaching the other stations also. Even so, there
are occasional events picked up by one of the other stations which are not picked up by
these two; but presumably the number of such events is not suÆcient to o set the extra
noise incurred by incorporating these further stations into the model.
The results from sets of three stations perform more poorly (in the sense of AIC) thanthe best pair of stations. Thus the best overall performance is found by combining the

110

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

Table 7.2: Performance comparisons for single and multiple predictor models in the 300
km region around Beijing
Station
H
H1
L
Q
S
HL
HQ
HS
LQ
LS
QS
HLQ
HLS
HQS
LQS
HLQS

AIC N+ =0
-15.73 24 0.53
-10.72 24 0.53
-14.45 26 0.51
-30.11 32 0.30
-45.41 34 0.33
-19.85 31 0.30
-33.40 34 0.27
-44.33 35 0.185-29.99 34 0.12
-42.72 33 0.28
-50.36 39 0.15
-30.66 34 0.12
-41.07 37 0.15
-48.32 40 0.13
-44.92 39 0.14
-41.80 40 0.13

The common time period of the analysis is [1826; 5477] days with a number of earthquakes,
N = 50. The AIC values are relative to that of the constant-rate Poisson models, whose
log-likelihood is 264.55. N+ is the number of predicted events. The parameter  is the
background seismicity level and 0 is the average occurrence rate. The bold face number
shows the minimum AIC value. The maximum likelihood estimates of the best modelare ^ = 0:002182(events=day ), c^1 = 0:0007786, ^1 = 0:8613, ^1 = 0:5467, c^2 = 0:002548,
^2 = 0:7902 and ^2 = 0:3543, where subscripts 1 and 2 stand for Qingxian and Sanhestations, respectively.

7.3. RESULTS FOR THE SEISMICITY WITHIN 300 KM CIRCLE AROUND BEIJING

111

signals from the stations which also happen to performing best in the single signal studies,
namely Sanhe and Qingxian.
It should be noted here that there are many other possible ways in which the signals
from the four stations could be combined within a joint model. In particular, one might
consider their combination through the multi-predictor model described by Utsu (1977,1982) and Aki et al (1981), based on the assumption of independent precursors. A preliminary report on the performance of the multi-predictor model was given by Ogata andZhuang (2001). They showed that, within the area common to the four dashed circles
in Fig 7.1, it gave a substantially improved performance from that based on any single
predictor. A careful comparison of the performance of a range of multi-signal models, including both the Lin-Lin and independent precursor models, with the electric signals data,
is currently in preparation. SuÆce it to say here that the Lin-Lin model with signals fromthe best 2 stations gives the best performance out of all the models so far considered. Its
superiority to the independent precursor model may be because in fact the stations do notrespond independently to pending earthquakes; some earthquakes are detected well by all
stations; others are not detected by any of the four. Thus, in addition to the sensitivitiesof the individual stations, there seem to be other, common, factors which determine how
well the stations as a group are able to pick up a particular event. This feature shows up
particularly clearly in the more detailed examination of earthquakes from the Tangshan
region given in Section 7.4, although it is not restricted to events from that region.
As in the preceding section, a visual impression of the performance of the best- tting
model can be obtained from Figure 7.6. As in the previous gures, events occurring attimes where (t)   of equation are marked by circles, the others by crosses. Table 7.3
shows the number of circle events compared to the number of cross events for the best
single stations (Sanhe), for the best pair of stations (Qingxian and Sanhe), and for thesame pair of stations with earthquake occurrence times randomly located on the time axis

112

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

Table 7.3: Numbers and proportions of earthquakes, and days above and below the mean
occurrence rate of earthquake
Station period(days)
t   (percent) t <  (percent) total
S[0,5874] NoEQ
1688 (28.73%)
4100 (69.78%) 5788
EQ
54 ( 0.92%)
33 ( 0.56%)
87
QS[0,5477] NoEQ
1738 (31.73%)
3656 (66.74%) 5394
EQ
62 ( 1.13%)
22 ( 0.40%)
84
QSR[0,5477] NoEQ
EQ

1771 (32.33%)
29 ( 0.53%)

3620 (66.08%) 5394
58 ( 1.06%)
84

In this table, \NoEQ" means the days without earthquakes, \EQ" means the days withearthquakes. t   and t <  respectively represents the days of conditional intensity
value above and below the mean occurrence rate of earthquakes. The results \S" arefrom the best singal station, Sanhe, and \QS" are from the best combination of stations,
Qingxian and Sanhe. For comparison, the last two rows, \QSR" show the results of
estimating the same parameters for the same combination of stations and the same numberof earthquakes, but with earthquake occurrence times simulated randomly on the time
axis.
(noted by QSR in the table). It shows also the proportions of the total observation period
taken up by values of t where (t) lies above and below the background rate. Again the
superior performance of the two-signal model is clear.
It is of interest to examine more carefully the spatial patterns for the predicted andunpredicted events. Several examples of these are shown in the plots making up Figures
7.3 and 7.8, where the circles and crosses have the same interpretation as in the previous
gures. The rst group of ve plots (Figure 7.3) represent epicentral maps of events within
the 300 km Beijing circle. They show which events were picked up by which station ofobservation, and which events were picked up by the combined model.
The next group (Figure 7.8) shows magnitude/distance plots for the predicted andunpredicted events for each of the four stations separately. It is evident from these thatthe larger events are picked up more successfully, even at larger distances. In fact most of
the events with magnitude greater than 5 are picked up by all the stations, and all suchevents/day

events/day

events/day

events/day

events/day
0.02

0.02

0.02

0.02

0.02

0.02

0.08
0.08
0.08
0.08
0.08
0.08

3200

2200

1200

200

3400

2400

1400

400

3600

2600

1600

600

4800

3800

2800

1800

800

5800

4600

5600

4400

5400

4200

5200

113

Figure 7.6: Plot of the tted conditional intensity function of the multi-predictor modelwith inputs from Qingxian (Q) and Sanhe (S) stations.
The dashed horizontal line represents the mean rate. The earthquakes are marked bycircles and crosses, with sizes representing their di erent magnitudes.

qxsh.0−5477

7.3. RESULTS FOR THE SEISMICITY WITHIN 300 KM CIRCLE AROUND BEIJING

events/day

41

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

40

(B)

37

37

38

39

40

Latitude(A)

38

Latitude41

CHAPTER 7.

39

114

112

114

116

118

120

112

114

120

118

120

41
39

40

(D)

37

37

114

116

118

120

112

114

116

Longitude(F)

HuailaiSanhe
BeijingLangfang

39

39

40

Latitude

(E)

40

41

Longitude
41

118

38

39

40

Latitude(C)

112

Latitude116

Longitude

38

Latitude41

Longitude

38
37

37

38

Qingxian

112

114

116

Longitude

118

120

112

114

116

118

120

Longitude

Figure 7.7: Epicentral locations of predicted and unpredicted events.
Predicted events are taken to be those which occur when the conditional intensity lies
above its mean, and are marked by circle. Other events are taken as unpredicted andmarked by crosses. (A) { (E) show patterns for the estimates by the data from Huailai;
Langfang; Qingxian; Sanhe; and Qingxian and Sanhe, respectively. (F) shows the locations of the stations.

7.3. RESULTS FOR THE SEISMICITY WITHIN 300 KM CIRCLE AROUND BEIJING

115

events are picked up by the combined signal from the best Lin-Lin model. These panels
are for Stations Huailai, Langfang, Qingxian and Sanhe, respectively. The distance e ects
are less pronounced, but there appears to be some tendency for the smaller events to be5

Magnitude

4

4.5

5
4

4.5

Magnitude

(B)

5.5

6

(A)

5.5

6

successfully picked up only at closer distances to the station.

0

100

200

300

400

0

200

300

400

(D)

5
4

4

4.5

5

Magnitude

5.5

6

(C)

4.5

Magnitude

100

Distance from the station

5.5

6

Distance from the station

0

100

200

300

Distance from the station

400

0

100

200

300

400

Distance from the station

Figure 7.8: Plots of magnitudes of the predicted earthquakes and unpredicted earthquakes
against their distances to the stations.
The results are respectively shown for the signal data from (A)Huailai, (B) Langfang, (C)Qingxian, and (D) Sanhe stations.
Finally, Figures 7.9 show that the successfully and unsuccessfully predicted eventsoccur more or less at random throughout the observation period, without any marked

116

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

time uctuations in the proportions of successfully predicted events.
In summary, we conclude from the results in this section that the e ectiveness of theforecasts for the Beijing region can be substantially enhanced by combining the signals
from several stations. Within the family of Lin-Lin models, the most e ective model, in thesense of AIC, results from combining the signals from Sanhe and Qingxian stations only.
The results are con rmed by examining the proportions of events successfully predicted(in the sense of having expected intensities above the background) in relation to the
proportions of time occupied by signals above the background.

7.4 Results for the Tangshan region
As in many earthquake studies, the results from the analysis overall are strongly a ectedby the behaviour of one special group of events. In the present analysis, such a role isplayed by the events from the Tangshan aftershock region (see Figures 7.1), which constitute approximately half of the events recorded in the larger region over the observationinterval. These are not direct aftershocks of the 1976 Tangshan event of M7.8 (thesehave been removed by the declustering procedures referred to in the introduction) but
represent continuing high activity within the general vicinity of the Tangshan event.
The locations of the epicentres within this region are shown in detail in Figure 7.1.
The panel of plots in Figure 7.10 indicates which of these events were picked up by which
station, in the sense described above. The most interesting feature is the distinctive
cluster of events, roughly in the middle of the epicentral region, which were not picked up
by any of the stations (nor, as a consequence, by the combined model). It is this groupof events which shows, more clearly than any other feature that we have examined, that
there are di erences among the events in the catalogue, not to be explained by magnitude
or proximity to the observation stations, which a ect their capacity to be detected byevents/day

events/day

events/day

events/day
0.02

0.02

0.02

0.02

0.02

0.02

3200

2200

1200

200

3400

2400

1400

400

3600

2600

1600

600

3800

2800

1800

800

5800

4800

5600

4600

5400

4400

5200

4200

RESULTS FOR THE TANGSHAN REGION

events/day

0.08
0.08
0.08
0.08
0.08
0.08

qxshr.0−5477

117

Figure 7.9: Correspondence between the same conditional intensity as in Figure 7.6 andevents with random time points.
The dashed horizontal line represents the mean rate. The simulated earthquake aremarked by circles and crosses, dependign on being either predicted or unpredicted, respectively, with di erent sizes representing the di erent magnitudes.

7.4.

events/day

118

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

observations on the electrical signals. The most obvious explanation is that the events
di er in some way { in their mechanisms, in their orientations, in their geological or
topographical situations, or for other reasons { which a ects their capacity to generate
electric signals. There is little further light that the present analysis can throw on thise ect.
In most other respects, the events within this region appear to be little di erent fromthose in the rest of the catalogue in regard to their predictability by the present methods.
A Lin-Lin model was tted just to the data from this region, but the parameters are very
similar to those of the overall model.
From those results, one could work out the probability that 4 out of the 23 Tangshan earthquakes went undetected by all 4 stations as follows. The probability that any
one event goes undetected by all stations, assuming independent chances of detection, is
10/237/239/2313/231/80. We now have 23 trials, with probability of failure 1/80
and four \failure". If we calculate the signi cance as probability (four or more failure),
which must been pretty small, we can draw the conclusion that the stations can hardly
be behaving independently.

7.5 Probability gains of forecasts based on electricalsignals
So far we have considered the models from the viewpoint of establishing their statisticalsigni cance. It is equally important, and raises somewhat di erent issues, to consider theirperformance directly from a prediction viewpoint. The probability modelling framework
adapted in the present paper is very easily extended to the development of probability
forecasts, which can then be assessed directly, and in their turn can form the basis fora variety of decision procedures (see, for example, Vere-Jones (1995, 1999) for discussion

40.0

(B)

39.6

Latitude

39.0

39.2

39.4

39.6

118.2

118.4

118.6

118.8

117.8

118.2

118.4

118.6

118.8

118.6

118.8

(D)

39.6
39.2
39.0

39.0

39.2

39.4

Latitude

39.6

39.8

(C)

118.0

Longitude

40.0

Longitude

39.4

Latitude

39.4
39.2
39.0

118.0

39.8

40.0

117.8

Latitude

119

39.8

(A)

39.8

40.0

7.5. PROBABILITY GAINS OF FORECASTS BASED ON ELECTRICAL SIGNALS

117.8

118.0

118.2

118.4

118.6

118.8

117.8

118.0

118.2

118.4

Longitude

(E)

39.6
39.4
39.0

39.2

Latitude

39.8

40.0

Longitude

117.8

118.0

118.2

118.4

118.6

118.8

Longitude

Figure 7.10: Epicenter distribution of the predicted and unpredicted earthquakes in theTangshan region (the magni ed rectangular area in Figure 7.1).
The results are from tting the earthquakes in the 300 km circle centred at Beijing and thesignal data from (A)Huailai, (B) Langfang, (C) Qingxian, (D) Sanhe, and (E) Qingxian
and Sanhe, respectively.

120

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

and further references). In this section we present a brief discussion of the performanceof probability forecasts based on the best Lin-Lin model from Section 7.3.
The quantity (t) in Equation (7.1) represents the conditional expectation for thenumber of events in day t, given the signals observed up until time period t

1. Theconditional probability that an event over the target magnitude will occur in day t, as it
might be issued at the end of day t

1, is therefore given byp(t) = 1 e (t) ;(7.4)

which is suÆciently approximated by (t) itself if this is small (say less than .05). The
probability gain

for an event occurring on day t is then the ratio p(t)=p wherep = 1 e and  is the mean rate derived from tting a simple Poisson model to the same data.
The performance of the model on a particular data subset of interest, such as theearthquakes in the Tangshan region, or those lying above a certain magnitude threshold,can then be assessed by computing the average probability gains for the events in question.
In such a context, however, it may be more appropriate to use the geometric mean rather
than the arithmetic mean to compute the average probability gain. Thus we look at thearithmetic means of the logarithms of the probability gains, say

H = (1=M )MX

1[log p(ti )=p]

where the summation is taken over the time intervals (days) with an event occurring, andcompute the average gain as

p~ = eH :Quantities such as H lie behind the use of the entropy score as a means of scoring aprobability forecasting scheme, as suggested in Vere-Jones (1999).

7.5. PROBABILITY GAINS OF FORECASTS BASED ON ELECTRICAL SIGNALS

121

Table 7.4: Probability gains
Class

H M

H

p~

Sanhe Station in period [0, 5874]
 4:0 28.493 87 0.3275 1.3875
 4:5 20.039 28 0.5003 1.650
 5:0 10.078 12 0.8398 2.3160
 5:5 2.984 4 0.7460 2.1085
Qingxian and Sanhe in period [0, 5477]
 4:0 39.790 84 0.4737 1.6059
 4:5 26.337 37 0.7118 2.0377
 5:0 12.045 11 1.0950 2.9892
 5:5 4.034 3 1.3447 3.8369

 4:0
 4:5
 5:0
 5:5

randomized occurrences-28.412 84 -0.3382 0.7131-9.050 37 -0.2446 0.7830-7.322 11 -0.6656 0.5140-3.033 3 -1.0110 0.3638

The
PM de nitions of the symbols can be found in the text, i.e., H is de ned by H =
i=1 [log p(ti )=p], H is the log probability gain and p is the average gain. Non's are forthe non-event days. For comparison, the last four rows shows the results of estimating
the same parameters for the same stations and the same number of earthquakes, but withearthquake occurrence times distributed randomly on the time axis.
Table 7.4 shows the mean probability gains p~ computed in this way for various classes
of interest from the present data. The improvement in the performance of the forecastingmodel for larger events is very apparent, although the numbers of events at higher magnitudes may be too few to make this a truly reliable comparison. Similarly, we may seethat the performance of the model for the Tangshan events does not di er greatly fromits performance on the rest of the events.
Another way in which the same ideas can be used is to compare the average probability
gains recorded for the actual events, with the gains recorded for a similar number of events

122

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

randomly located on the time axis. In this case the average should be taken for predicting
both successes and failures, where for the days on which no events occur, the probability
gain is taken to be the ratio of the probability for no event occurring, as given by themodel, to the probability of no event occurring in the underlying Poisson model. By
repeating the simulations a large number of time, a histogram can be built up for thescores to be expected from a random model, and the score for the actual model can becompared to this. Such a histogram is shown in Figure 7.11. It provides further strong
con rmation of the reality of the explanatory power of the electric signals.

150
100
0

50

Relative Frequency

200

250

Histogram from 1000 simulations−0.020−0.015−0.010−0.005

0.000

0.005

0.010

Mean probability gain

Figure 7.11: Histogram of mean probability gains. The mean probability gain from theactual model is mark by a crossed circle.

7.6 E ect of the change in the background seismicity
In this section we pick up some additional questions relating to the analysis of the electricalsignals data. The rst among these is question noted in the introduction of the e ect on
the results of the change in seismicity following approximately the Datong events in 1991.
In the models so far considered we have assumed a constant rate for the null model(Poisson model) and a constant background rate  for the Lin-Lin model. The former
assumption, in particular, appears inappropriate in the light of Fig 7.3 and could beconsidered to give an unfair bias against the null model. To examine this e ect more

7.6.

EFFECT OF THE CHANGE IN THE BACKGROUND SEISMICITY

123

Table 7.5: Comparisons of AIC values for the 1-stage and 2-stage background models
and the 1-stage background models
StationL
Q
S
Q+S
Two-stage Poisson -9.52
-9.52
-9.52
-9.52
One-stage Linlin -31.41 -45.51 -53.78 -63.25
Two-stage Linlin -42.91 -40.54 -53.54 -60.89
These values are obtained by substracting the AIC value for the 1-stage Poisson modelfrom the AIC values of the three other models, all tted the earthquakes occurring in
the time inteval [0,5477]. The patterns of the conditional intensity for this 2-stage Linlin
model are similar to those for the 1-stage Linlin model.
carefully, we divided the observation region into two parts at around the time of the
Datong event, and allowed the Poisson rate to be determined separately for the two subperiods. Table 7.5 shows the AIC values for the two-stage models and one-stage models,
which are di erence from the AIC values of the constant rate Poisson process. Specially,
this table compare AIC values of the Lin-Lin models with a single background rate , andwith separate rates 1 and 2 for the two periods. From these results we see that even
if the background rate is held constant, the Lin-Lin model based on the electrical signals
still outperforms the two-stage Poisson model. Indeed, the performance of the two-stage
Lin-Lin model is no better than the one-stage Lin-Lin if we consider further penalty due
to the existence of a change-point in the background rate (cf. Ogata, 1999, Section 4.2).
More generally, we found that allowing for separate rates improved the performance also
in other respects - for example, it improved the stability in the t to the Huailai station,and brought its parameter values closer into line with those from the other stations. Thereis not enough data in the second period to allow a full re tting of the model for this period
only, but there is no indication at this stage that the performance of the electrical signals,
or the form of the response functions, have been a ected by the change of seismic regime
following the Datong earthquakes in 1991.

124

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

Table 7.6: Fitted results by using the electric signals as the main process and earthquakes
within 300km to the station as the external process
StationH
L
Q
S

Ns Orders AICs Orderd AICd Ordere AICe247
0-343.3
1/1 -359.18
0-34.50366
1-529.65 1/0 -568.07
0
-6.92
246
1
-525.53 1/0 -525.51
0
1.28
267
1
-488.98 1/0 -486.98
0-1.2

Ns is the number of signals. Orders is the order of the Laguerre polynomial used forthe best tted model of self-excitation processes. Orderd are the orders of the Laguerre
polynomials used for the self-excitation term and the external-excitation term in thebest tted model with both excitation terms. Orderd are the orders of the Laguerre
polynomials used in the best model with external-excitation term only. AICs, AICd and
AICe represent the AIC-values for the the models with only self-excitation terms, withboth excitation terms and with only external-excitation terms, respectively. Bold faces
give the best models for each station.

7.7 Testing whether the electric signals are postseismic phenomena
One may question whether the electric signals are of postseismic e ects or not. To investigate this, I selected the large signals with strengths  200mV/ms as the events attachedat the ends of the days of observation, and also selected the earthquakes within 300 km
as the external process. I tted the original LINLIN model in (6.2) to these data. The
tted result are shown in Table 7.6 and Figures 7.127.15. It is clear to see that for all
the stations except for the Huailai case, AIC selected the model with self-excitation but
without external-excitation. This means that, for the three stations, the signal anomalies
are clustered but not postseismic. For the Huailia case, in spite of the di erent choice,when we plot the condition intensity function as shown in Figure 7.12(b), we can nd thatexternal excitation has very long term response, expressing the trend of the frequencies ofthe signals detected at the Huailai station, is shown in Figure 7.12(d). This can be bettermodelled if the trend term is added to the model (6.2) as was implemented in Ogata(1983) and Ogata and Katsura (1986).

7.8.

SUMMARIES OF CHAPTER 7

125

In summary, these results imply that the electric signals are clustered, and earthquakes
seem to have no e ect to the occurrence rate of the large electric signals. Nishiziwa, Leiet al. (1994) concluded the similar relationship of earthquakes with the seismo-electric
signals observed in Greece (see, e.g., Varotsos and Kulhanek, (1993); Varotsos and et al,1996). These consistent conclusions show the changes of underground electric eld may
be caused by the stress changes before big earthquakes.
I also treated the earthquakes as the target sequence and the electric signals as theexternal sequence, and found that the AIC selected the models in which earthquakes are
not only clustered but also in uenced by the electric signals (Zhuang 1999, unpublishedmanuscript).

7.8 Summaries of Chapter 7
The main conclusions of the present chapter are summarized as follows.
1. For all four stations considered, the incorporation of information from the electric
signals, at least as formalized in the Lin-Lin model, results in a substantial improvement over the null hypothesis of a Poisson model, at both 200 km and 300 km
radii.
2. The performance of the model over the common region centred on Beijing is considerably enhanced by incorporating the information from several stations jointly. In
this particular case, although each station contained relevant information, the mosteconomical and e ective model, from the AIC point of view, was obtained by incorporating the signals from the two stations Sanhe and Qingxian. Once the signals
from these stations had been taken into account, only marginal improvements camefrom incorporating information from the other stations.

126

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

(a)
3
4
5
6
0

1000

2000

3000

4000

1000

2000

3000

4000

(b)
0.40

0.20

0.00
0

0.20

0.003

(c)

0.15

(d)

0.002

0.10
0.001

0.05
0.00

0.000
0

10

20

30

0

10

20

30

Figure 7.12: Results for Huailai Station using the electric signal events as the target point
process.
(a) Occurrence times of the electric signals events (upper) and earthquake events (lower);(b) conditional intensity; (c) the impulse response function of self-excitation; (d) the
impulse response function of external-excitation. The time axis in (a) and (b) representsdays from 1987/5/1.

7.8.

127

SUMMARIES OF CHAPTER 7

(a)

3
4
5
6
0

2000

4000

6000

2000

4000

6000

(b)
0.40

0.20

0.00
0

(c)

0.30
0.20
0.10
0.00
0

10

20

30

Figure 7.13: Results for Langfang using the electric signal events as the target point
process.
(a) Occurrence times of the electric signals events (upper) and earthquake events (lower);(b) conditional intensity; (c) the impulse response function of self-excitation. The time
axis in (a) and (b) represents days from 1982/1/1.

128

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

(a)

3
4
5
6
0

2000

4000

6000

2000

4000

6000

(b)
0.40

0.20

0.00
0

(c)

0.20

0.10

0.00
0

10

20

30

Figure 7.14: Results for Qingxian using the electric signal events as the target point
process.
(a) Occurrence times of the electric signals events (upper) and earthquake events (lower);(b) conditional intensity; (c) the impulse response function of self-excitation. The time
axis in (a) and (b) represents days from 1982/1/1.

7.8.

129

SUMMARIES OF CHAPTER 7

(a)
3
4
5
6
0

2000

4000

6000

2000

4000

6000

(b)
0.40

0.20

0.00
0

(c)
0.20

0.10

0.00
0

10

20

30

Figure 7.15: Results for Sanhe using the electric signal events as the target point process.
(a) Occurrence times of the electric signals events (upper) and earthquake events (lower);(b) conditional intensity; (c) the impulse response function of self-excitation. The time
axis in (a) and (b) represents days from 1982/1/1.

130

CHAPTER 7.

STATISTICAL ANALYSIS ON THE ULF ELECTRIC SIGNALS

3. Probability gains from the 1-day predictions based on the electrical signals and the
Lin-Lin varied from an average of around 3.8 for events with magnitudes M  5:5
to around 2.7 for events with magnitudes in the range 4  4:5. All of the eventswith M  5:5 occurred at times with enhanced probability estimates; approximately
70% of the events with magnitudes in the range 4  4:7 occurred in periods withenhanced probabilities.
4. While most events were picked up by most of the stations, some events were picked
up by none of the stations, more than could be accounted for by an independent
precursor model. The reasons for these changes in detectability, which stood out
particularly clearly among the Tangshan events, remain to be clari ed.
5. While there are several indications that the stations pick up signals more e ectively
at shorter rather than at longer distances, a systematic study of the spatial e ectsof the signals remains to be undertaken.
6. The ULF electric signals are clustered, but not excited by the occurrence of earthquakes.

Part IV

APPENDIX

131

Appendix A

Observations on the ULF Ground
Electric Field in China

A.1 Introduction
Seismologists began to study the phenomenon of electro-magnetic radiation before earthquakes in the 1970's, using both eld observations and laboratory experiments. The aims
were to try to develop an improved understanding of the physical mechanisms at work,and to nd better observable quantities for use in short-term earthquake prediction. Since
then, precursory electro-magnetic anomalies have become important topics in earthquake
prediction research throughout the world (Yoshino and Tomigawa, 1990; Oike and Ogawa,1994). Observations and research on electro-magnetic radiation before earthquakes werestarted in China in 1976, after the great Tangshan earthquake (MS 7.8). Up to the present,
more than 130 observation stations have been built in China. Most of those stations observe the electro-magnetic radioation in a ultra-low frequency (ULF) band of 0{20Hz,
although a few of them observe at point frequencies between several kHz and 100 kHZ.
At these stations, buried electrodes are used to observe the electric eld components,
while electric induction instruments are used to observe the magnetic eld components.
At most stations, the records are produced by analogue pen recorders, although in a few
stations records are made using a single-board computer or micro-computer.
133134APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA
The observation network described in the present paper is the network in the Beijing region and its vicinity. These stations observe the electro-magnetic radiation at
frequencies in the 0.1{10 Hz (ULF) range. During 20 years' observation, there have been
many instances where anomalous variations in the electro-magnetic eld were recorded
before the occurrence of a local earthquake with M  4. For example, during the period
1/1/1982{1/10/1998, evident anomalies were observed at the Sanhe station on 270 daysout of a total of 5874 days. During the same period, 82 earthquakes occurred with M  4,
and with epicenters within 300 km of the station. Of these, 54 were preceded by evident
anomalies. The accompanying paper (Zhuang, Ogata and et. al.,2002) gives a detailed
statistical analysis of the observation data. The present paper is intended to provide abackground to this analysis; it contains a brief description of the observation system anddata processing procedures used within the Beijing network.

A.2 Monitoring network and observation system
A.2.1 Network distribution
The network in the Beijing region and its vicinity was started in 1981. The 8 stationsin this network are located at Langfang, Sanhe, Wen'an, Qingxian, Xincheng, Ningjing,
Hualai and Changli (the locations shown in Figure A.1). The Langfang, Sanhe, Wen'an,
Qingxian, Xincheng, Ningjing stations were set up in 1981{1982, Huailian in 1987, andChangli in 1990. The stations were set up with the aim of monitoring uctuations in the
ULF electro-magnetic radiations around Beijing.

A.2.2 Observation system
The observation systems are similar at all the above stations. To illustrate the systems
in use, we brie y describe the E-EM system, designed by the Provincial Seismological

A.2.

135

MONITORING NETWORK AND OBSERVATION SYSTEM

116

114

118

120

41

Huailai
40

Beijing
Sanhe
Changli39

Langfang
Xincheng
Wen’anQingxian

38

Ningjing
Legends

37

main faults
stationFigure A.1: Distribution of observation stations. The lines denote the locations of major
faults136APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA
Bureau of Hebei. More detailed illustration and technical parameters can be found in
Chen, Du and et. al. (1998). The basic structure of the E-EM system is shown in Figure
A.2(a). The electric signals detected by the sensor are transmitted to the \pre-processor".
After removal of some of the noise, the signals are transferred to the recorder.

Sensor The electrodes of the HBEMD-3 type sensor used to detect the electric signals
are made of stainless steel (Gr18Ni9C), and are cylinder shaped with a height of 300mm
and a radius of 40 mm. They are usually buried from 6 meters to 12 meters deep and
40 meters apart. Usually two pairs of electrodes are installed at a given station, along
perpendicular axes. The transmission cables are buried 60{80 cm deep, and no power
supply wires, communication cables or metal cubes are present nearby.
The sensor for observing the magnetic eld consists of three electric inductives witha central frequency 0.8 Hz, installed in three perpendicular directions, north-south (NS),east-west (EW) and up-down (UD), usually buried 3 meters underground.

Terminal system The recording method is by automatic continuous pen record. The
recorder is modi ed from DJ-1 recorder of the DK-1 seismometer. The DK-1 seismometerand its DJ-1 record is designed to record the medium-long period components (0{10HZ) of ground movement, widely used in Chinese seismological stations. More detailed
information of the instrument can be found in Chen, Du and et al (1994).
After the ultra-low frequency electric signals are received, they are transmitted throughthe transmission cables and pre-processors ( including impedance matching unit, lter)to the ampli er. The main functions of the \pre-processor" are to prevent 50 Hz noiseproduced by industrial power supplies from coming into the ampli er, and to prevent thepolarization potential from changing the working point of the ampli er. Another of its
functions is to assist in impedance matching. Then the ampli er ampli es the signals

A.2.

MONITORING NETWORK AND OBSERVATION SYSTEM

137

in order to enable the record-pen to move. This is the way in which the waveforms ofelectric signals are recorded. The dominant frequency band for the overall system is 1{10Hz, while the input resistance is greater than 10M . The noise level is less than 0.15 mV.
The record sensitivity is 0.5 mV/mm. The record speed is 60mm/min.
The frequency response function of the system is shown in Figure A.2(b), from which
we can see the frequency band is 1{10 Hz.

A.2.3 Methods for noise reduction
Techniques used for reducing noise include carefully choosing the observation frequency
band, burying the electrodes and transmission cables deep underground, ltering thereceived signals, using a high input resistance, choosing appropriate instrument sensitivity,and improving the quality of the electric power supplies.

Choosing the observation frequency band
The main frequency bands used for electro-magnetic radiation observations in China are,according to the standard classi cation of wireless electric waves,: (1) ULF, <30 Hz; (2)ELF, <3 kHz; (3) VLF, 3{30 kHz; (4) LF, 30{300 kHz.
To obtain reliable precursory electro-magnetic signals, it is important to take into
account the frequencies of commonly occurring noise sources, so that the noise can be
avoided or reduced.
Observations and experiments on rock mechanics show (Cress and Braoly, 1987; Gershengon and Gokhberg, 1989) that electro-magnetic radiation before an earthquake occurs over an extremely wide range of frequencies, covering virtually the entire spectrumof electro-magnetic waves. In principle, it seems that almost any frequency band could be
used for the observation of such signals. But the frequency bands higher than 30 kHz are

138APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA

(b)

0.01

Response
0.05
0.2

0.5

1

(a)

0.05

0.2 0.5 1 2Periods in seconds

5

Figure A.2: Illustration of the observation system (a) and the response function (b).

A.2.

MONITORING NETWORK AND OBSERVATION SYSTEM

139

widely used for broadcast radio, television, communication, etc., which produce signi cant noise for seismic electro-magnetic radiation observations. Since signals of very low
frequency (< 30Hz) decay slowly during propagation and are not much used for the above
purposes, there is little noise in this frequency band and it is ideal for the observation ofseismo-electric signals. Within this band, signals of ultra-low frequencies, from 0.01{20Hz, have least noise. Experiments con rm that it is better to make observations at lower
frequencies than at higher ones (Li and et al, 1989).
Signals of ULF and VLF also have the feature of minimum decay. This is because thehigher the frequency of the electro-magnetic wave, the more it is absorbed by the mediaof the earth.
In summary, to avoid noise and so obtain reliable seismic electric signals, it is betterto observe in the ULF or ELF bands, or at points in the VLF and LF bands which haveno interference.

Deeply burying the sensor and screening the transmission lines
As the electro-magnetic signals come from deep under the ground, buried electrodes arechosen as the sensor to detect or extract the alternating component from the quasi-static
ground electric eld between the two electrodes. The main part of the sensor is buried toa depth > 6 meters, in order to use the ground to screen the sensor from noise sources inthe sky and on the earth's surface. In addition, the wires used for signal transmission arescreened by high quality metal nets, are buried 0.6{0.8 meters below the surface, and go
directly to the laboratory from underground. These policies help to improve the anti-noiseproperties of the transmission from the sensor to the recording instrument.

140APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA

Cutting o direct current, ltering, and impedance matching
As the body of the sensor of the electrodes is made of high quality Cr18Ni9C stainlesssteel, only very small polarization electric potential is produced between electrodes. To
prevent such polarized electric current from coming into the ampli er and then changing its working point, 1000 F electric capacities are connected serially into the systemfor cutting o the direct currents. Also dual T bandpass lters together with network
impedance matching methods are adopted to lter out the noises caused by the 50 Hz
signals from power supplies and by the unmatched impedance. These three steps areaccomplished by careful design of the pre-processing parts. They improve the anti-noiseproperties and stability of the recording system.

Increasing input resistance
Increasing the input resistance helps to decrease the in uence of ground resistance on the
instrument sensitivity.

Adjusting the sensitivity of the recording system
This is a diÆcult task and has to be done step by step over quite a long period. Becauseof di erences in geological or tectonic structures, the amplitudes of the received signals
may di er substantially in both size and occurrence frequency. Too high a sensitivityof the instrument may cause too much information to be processed, including a lot ofnoise. Thus, it takes time to adjust the sensitivity of the instruments. The aim, in ideal
circumstances, is to calibrate the system so that no signals are recorded if no earthquakes
with M 3.0 occur within 100 km of the station, while signals are recorded whenever an
event with M 3.5 occurs within this distance.

A.3.

POSSIBLE NOISE SOURCES

141

Improving the quality of the power supply
The commercial power supply usually has the problems of unstable voltage and electri ed
ground wire, which may distort the recorded signals. Thus, it is important to improve thequality of the power supply. Here we use a stabilization plant for the electricity supply,
while power from automatically recharged batteries is used when the power from thestabilization plant is cut o .

A.3 Possible noise sources
A.3.1 Recognition of interference signals
Many of the noise sources which may obscure or distort observations on the seismic
electro-magnetic radiation can be avoided or decreased by adopting suitable anti-noise
procedures. Experience has shown that moving iron objects, low-sky thunder and lightning, and strong magnetic storms may in uence the magnetic observations; while low-sky
or sky-to-ground lightning, and power lines that leak electricity to the ground near thestation, may in uence electric eld observations made using the buried electrode method.
With the exception of the e ects of magnetic storms on the magnetic component, which
are diÆcult to recognize, the other major sources of noise can be distinguished easily
from the seismic signals by the di erences in character of the recorded waveforms. As an
example, Figure A.3 shows the waveforms of two noise signals, one from ground strong
lightning and one from electricity leaked by power lines, recorded by the Langfang station.

142APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA

Figure A.3: A photocopy of a record with noises (the upper panel) from ground strong
lightnings and (the bottom panel) from electricity leaked by power lines.
In the upper panel, the short vertical bars, which form up a line intersecting the horizontal
line, are the timing marks at every minute, and irregular pulses caused by strong groundlightnings are marked by arrows. In bottom, the part in black represent the period of
power leakage.

A.3.

POSSIBLE NOISE SOURCES

143

A.3.2 Analysis of the e ects of air temperature, rainfall andlightning
In uence of air temperature and rainfall on the electro-magnetic signals
In the vicinity of Beijing, the climate is hot in summer and cold in winter, thus shows aclear 1-year cycle with a temperature di erence of 50Æ C between the maximum and theminimum. About 80% of the rainfall is concentrated in the summer months of June, Julyand August. However, from Figure A.4, the electro-magnetic signals do not show any
yearly cycle, nor is there any di erence between the appearance of the signals in summer
and in the other seasons.

Lightning
To study the in uence of relatively distant lightning on the signals observed at a particularstation, four stations - Sanhe, Langfang, Xincheng and Ningjin - were selected for study,
covering between them a region of about 320 km from north to south and about 200
km from east to west. Lightning data from meteorological stations within this region
were then collected. Before the relationship between lightning and electric signals wasanalysed, the lightning data was rst classi ed into three groups:(a) if lightning occurred near one station only, then it was called \casual lightning";(b) if, during the storm, lightning occurred near 2 or 3 stations, then it was called \local
lightning";(c) if, during the storm, lighting occurred near all four stations, it was called \regional
lightning".
Table A.1 shows the strengths of lightning storms in the 6 years from 1987 to 1992 in thestudied region.

144APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA
Table A.1: Number of lightnings occurring in the monitoring region during 1987{1992Year Casual Local Regional total months
1987
25
31
13
69
4{9
1988
31
31
6
68
4{10
1989
22
18
5
45
4{10
1990
27
34
10
71
4{9
1991
26
31
4
61
4{10
1992
39
32
3
64
4{10
The table shows that lightning mainly occurs in the months April-October each year.
It can be shown that the average number of days with lightning is about 63 days, witha standard deviation of 9 days. Thus the lightning days take up about 1/7 of the wholeyear.
A 3-day moving average of the electric-signals from 1987{1992, was then calculated,
and is shown with the lightning data in Figure A.4. The gure shows that the lightning
usually occurs in summer and autumn, but the variation of electric signals shows noseasonality, implying that there is no direct relationship between electric signals andlightning.

A.3.3 Magnetic storms and electric signals
Magnetic storms are a source of electro-magnetic noise from outer-space, and can produce strong e ects on the magnetic component of the electro-magnetic radiation. Here,however, we study the in uence of such storms only on the ULF ground electric eld.
Table A.2 shows magnetic storms recorded by Baijiatan Station, Beijing, over a 6-year
period. The storms have again been roughly classi ed according to their severity. In 1987{
1992, about 550 magnetic storms occurred, on average about 92 days every year.The e ectof such storms must be examined carefully, because they occur in about one quarter ofthe days in each year, and produce signals with a very wide range of frequencies, including

A.3.

POSSIBLE NOISE SOURCES

145

Figure A.4: 3-day average of the electric signals (solid curve), the lightnings (bar plot)
and the earthquakes (circles) with MS > 5:2 in the monitoring region.
The di erent lengthes of bars represents "local lightning", "evident lightning" and "extreme lightning" from shortest ones to longest ones.

146APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA
Table A.2: Number of magnetic storms occurring during 1987{1982 recorded at Baijiatan
Station, Beijing
Year Clear Evident Extreme total months
1987 58
12
0
70
1{12
1988 73
10
4
87
1{12
1989 68
14
13
95
1{12
1990 52
15
7
74
1{12
1991 89
9
24
122
1{11
1992 74
16
12
102
1{12
Table A.3: Correlation between magnetic storms and electric signals
Clear
Evident
Extreme
Total Related Unrelated Total Related Unrelated Total Related Unrelated
414
26
388
76
6
70
60
5
55
6%
94%
8%
92%
8%
92%
some very close to the observation frequency band. In the 6 year period of the study,
magnetic storms occurred in every month except December 1990. But the anomalies did
not occur regularly throughout the whole year.
To examine this e ect further, the following analysis was carried out: if an electromagnetic radiation anomaly occurred in the two days after the occurrence of the magnetic
storm, the two were counted as related events; otherwise, not. The results from 1987 to
1992 are shown in Table A.3. The table suggests that the occurrence times of the magnetic
storms and the anomalous electric signals are unrelated.
Figure A.5 shows in more detail the correlation between the uctuations of the electric
signals and the magnetic storms from 1987 to 1992. The moderate or strong earthquakes
occurring in this area are also plotted in this Figure, in order to show the links between
the di erent features. From the Figure it can be seen that the magnetic storms occur
frequently along the time axis, with their own occurrence pattern, no matter whether
earthquakes occur or not.

A.3.

POSSIBLE NOISE SOURCES

147

Figure A.5: 3-day average of the electric signals (solid curve), the magnetic storms (barplot) and the earthquakes (circles) with MS  5:2 in the monitoring region.
The di erent lengthes of bars represents "local magnetic storm", "evident magnetic
storm" and "extreme magnetic storm" from shortest ones to longest ones.

148APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA
From 1987 to 1992, the electromagnetic signals show ve groups of evident anomalies(peaks of 3-day averages of electric signals in Figure A.5), which were followed by amoderate or strong earthquake. On the other hand, \extremely strong" magnetic storms
occurred 60 times, and on only 5 of these occasions were clear electro-magnetic signals
observed in the two days after the storm. Further, in March 1989, April 1990, June, July,
August and November 1991, February, September and October 1992, many \extremely
strong" magnetic storms occurred, but not a single clear anomaly was observed in theelectric signals recorded during the intervals from 2 to 10 days after the magnetic storms.

A.4 Anomaly recognition
A.4.1 About the data
Over the 16 years from 1982 to the end of 1997, this observation network has worked
continuously, and has accumulated a large amount of data.
There is almost no noise in the background of the ULF ground electric eld unlessa source of strong low frequency electric signals is located within 100 meters from thestation.
On most days, the record is a straight line, with almost no noise component mixed in.
Usually clear or distinguishable electric signals start to appear several days or sometimes
even up to several weeks before the occurrence of a suitably located earthquake - either
a MS 3.5{4.4 event within 100 km of the station, or a MS 4.5{4.9 event within 200km of
the station, or a MS 5{6 event within 300km of the station.
In general, the durations of these anomalies tend to get longer, the amplitudes larger
and the oscillation frequencies higher, when the magnitude of the subsequent earthquake islarger, or the station is closer to the epicenter. At least half the stations showed anomalies

A.4.

ANOMALY RECOGNITION

149

before at least 10 of the earthquakes with MS > 5, and at some stations the anomalous
signals became more pronounced just before the earthquake occurred. Precursors before

MS 3.5{4.5 earthquakes tend to be seen only near the epicenter, and the duration times
and signal amplitudes tend to be smaller than for the larger events (Guan and Liu, 1995).
The observations show that the anomalies take up less than 2% of the whole recording
time and that most of them occur within a short time { usually only a few days { before
earthquakes. Again we refer to the companion paper (Zhuang et al, 2002) for a detailed
statistical analysis of these results.

A.4.2 Classi cation of signals
The waveforms of the electro-magnetic signals from the pen-record data are complicated.
We distinguish 4 classes according to the shape of the wave form, the dominant frequency,
and the relationship to subsequent earthquakes, namely, pulses of classes A, B, and C,
and other miscellaneous waveforms (Figure A.6). The rst three classes tend to have aclose relationship with moderate or strong earthquakes occurring within 300 km from thestations.
Their features are as follows:
Class A: Waveforms in Class A are oscillations within an exponential envelope, having
amplitudes of 1-4 mV and periods of 0.05-0.1 second; most are quasi-equal-amplitude
oscillations;Class B: Waveforms in Class B are positive or negative pulses with varying amplitudes
from 0.5mV up to 14 mV, and periods in the range 0.4{3 seconds;
Class C: Waveforms in Class C are quasi-pulses like the sine function, with smooth peaks
and valleys, and usually with amplitudes in the range 1{5 mV and periods in therange 3{6 seconds.

150APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA

A.5 Data processing
To date, the strength of an anomalous electric signal has been measured by two quantities,
the duration time of the anomaly, and its average amplitude. The details are as follows.
From the records, the following quantities are calculated:
1. V = average amplitude (mm)  instrument scaling factor (mV/mm), where average amplitude = (minimum amplitude + mean of the amplitude + maximum
amplitude)/3 corresponds to the voltage di erence between the two electrodes.
2. E = V=L, where L is the distance between electrodes (m) and E is the electric eld(mV/m).
3. A = T E (s mv/m), where T is the duration time of the signal in seconds.
4. The daily strength is Atotal = A1 + A2 +    + AN , where N is the total number ofsignals recorded in a day.

A.6 Summary
This paper describes the techniques used in the Beijing network to observe anomalousuctuations in electric signals in the ULF frequency range, and outlines the main methods used to reduce noise and to improve instrument stability. Only big lightning eventsoccurring close to the station, and electricity leaking to the ground from other sources
close to the station, appear able to in uence the signal record, and these noise sources
can be identi ed easily by their waveforms. The signals are independent of variationsin temperature, rainfall, distant lightning and magnetic storms. Although the magnetic
components of the signals are in uenced by strong magnetic storms and by low-sky lightning close to the station, the latter can also be recognized from the waveforms. To get the

A.6.

151

SUMMARY

most reliable electro-magnetic radiation signals, the observation network needs to be set
up with the capacity to record uctuations in both the electric and the magnetic elds.

Background

Class A waveforms

Class B waveforms

Class C waveforms

Figure A.6: Normal record (the top panel) and classi cation of anomalies: Class A (thesecond panel), Class B ( the third panel) and Class C (the bottom panel).

152APPENDIX A. OBSERVATIONS ON THE ULF GROUND ELECTRIC FIELD IN CHINA

Appendix B

On the Optimization Algorithm

This appendix describes an gradient free local optimization algorithm. The quasi-Newton
type methods are often used to solve a local optimization problem. In practice, computing program implementation of a quasi-Newton type optimizer always need to include asubroutine for calculating the gradients or even the Hessian matrix. For example, theoriginal Newton method need to calculating both the gradient and the Hessian matrix
explicitly, and the Davidon-Fletch-Powell method uses some technique to construct anquasi-Hessian matrix after the gradient of the target function is calculated explicitly. Of
course, we can use some numerical approximations instead of the explicit calculation tosave up time in programming. But a shortcoming for the conventional di erential algorithms is that, while the step get very small at the a local minimum or a local maximum,
the approximation diverges, and on the other side, a large step puts some limitation in
the precision of the solution.
Taken these into consideration, the following discusses an gradient free algorithm for
local optimization. Instead of calculating the Hessian matrix or quasi-Hessian matrix, thisalgorithm uses a quadric form to t the function near the minimizing solution, and use theminimizing solution for this quadratic form as an approximation around the minimizing
solution of the target function.
153

154

APPENDIX B. ON THE OPTIMIZATION ALGORITHM

The problem of local optimization is to search the minimum of the function f (x ) andits solution x0 . Here x denotes a n-dimensional vector, i.e., x= (x(1) ; x(2) ;    ; x(n) )T in

a n dimensional Euclidean space, x(i) being the i-th component of x , i = 1;    ; n, andupper-script T representing transposition of a vector.
If x is close enough to the minimizing vector x 0 , by Taylor series expansion,@ff (x )  f (x 0 ) +
@ x

!T

x =x 0(xx 0 )@2f@ x @ x T

+

1
(x
2 

x 0 )T= f (x 0 ) +

1
(x
2 

x 0 )T H0 (xx =x 0(xx 0 )x 0 );

(B.1)
(B.2)

where the rst order derivatives disappear because

@f
@ x

!

x =x 0= 0;

and H0 is the Hessian matrix such that

H0 =Equation (B.2) can be rewritten asf (x ) 
=

1
2

x

1
2

x

1
1

x 0
x 0

where

He =

@2f
@ x @ x T!T!T

x =x 0:

H0 0
T 2f (x)

0
He

x

1

H0 0
T 2f (x)

0

x 0

!

x

!

1

x 0

;

!

(B.3)

!

:

(B.4)

If S is a linear transformation on a n + 1 dimensional Euclidean space, which maps

Q:

U =

x

1

x 0

!

! V

=

x

1

x p

!

;

(B.5)

155
where x p is a given point, usually chosen as the point close to the minimum in the known
points.
Assuming that Q is invertible, substitute U = Q 1 V into (B.3), and
1 T e
H U
2 U
1
1V
e
= V T (Q 1 )T HQ

2
1 T e
= V G V ;2f (x ) 

(B.6)
(B.7)
(B.8)

where
1
e
Ge = (Q 1 )T HQ(B.9)

If we have calculated the function values of f at a set of m points, x 1 , x 2 ,    , x m , thefollowing equation can be obtain from (B.8),
8
>
>
>
>
>
>
>
>
>
>
>
>
<

f (x 1 ) =f (x 2 ) =

>
>
>
>
>
>
>
>
>
>
>
>
: f (x
m) =

x 1
x 2
x m

1
1

1

x p
x p
x p

!T

x 1

Ge

!T

x 2

Ge
..
.
!
T

x m

Ge

1
1

1

x p
x p
x p

!

+ deviation!

+ deviation:(B.10)!

+ deviation

The above equation is a linear equation with (n + 1)(1 + n=2) unknown parameters, allbeing the elements of the symmetric matrix Ge.
If m  (n +1)(1+ n=2) and all the deviations are set 0, we can get the general solutionof (B.10), which can be regarded as an approximation of Ge, denoted by G.
The next task is to nd the local minimum of the approximation form,1f (x ) 
2

x

1

x p

which can rewritten in the form of1
g (y ) 
2

y

1

!T

!T

G

G

x

y

1

1

x p

!

;

(B.11)!

;(B.12)

156

APPENDIX B. ON THE OPTIMIZATION ALGORITHM

where

y x

=

x p :(B.13)

Di erentiate (B.12) with respect to y , and at the minimum y  ,

@h
@ y

=G
y = y



y

1

!

=0 :

(B.14)

If the solution y of (B.14) is obtained, then let x  =x p + y  . This x  can be regarded

as a better approximation of the local minima x 0 . Put x  into the point set of the known

observation taking instead of the point having the largest function values. By repeating
the present algorithm, a better and better approximation of x 0 can be obtained.
There are two ideas to avoid (B.10) being ill-conditioned. One idea is that x p can bechosen as the best point (the point with the smallest function value) in the known points.
The other idea is that, after x m+1 is obtained, we also generate a number, namely k, ofn-dimensional random vectors from the probability density of the multi-variate normal
density with a mean of xm+1 and a standard variance matrix of s2 I , where s is themaximum distance between the m + 1 known points and I is the n  n unit matrix. We
put these k generated points and x m+1 into the point set of the known points, and selectthe m points with the smaller function values as the new point set for the next iteration.
In summary, the present algorithm can outline below.

Problem: Find x 0 minimizing f (x ).
Algorithm:
1. Given a number m greater than (1 + n)(1 + n=2), generate m n-dimensional multivariate random variables, x 1 ; x 2 ;    ; x m , from a normal density with a mean of
0and a reasonable variance. Let P = fx 1 ; x 2 ;    ; x m g. Calculate the function

values of f (x 1 ); f (x 2 );    ; f (x m ).

157
2. Find the point with the smallest function value in P , say x p .
3. Generate k n-dimensional multivariate random variables, x0 1 ; x0 2 ;    ; x0 k , from anormal density with a mean of xp and a variance matrix s2 I , wheres = x max k x i i ;x j 2PnXx j k = max

x i ;x j 2P

l=1

(x(il)

x(jl) )2

!1
2

(B.15)and I is the n  n unit matrix. Calculate f (x0 1 ); f (x0 2 );    ; f (x0 k ).
4. Find the m points with the smaller function values in the point set P;    ; x0 k g, and let P be the point set of these selected points.

S

fx0 1; x0 2

5. Find the point with the smallest function value in P , say x p .
6. Let y i =x ix p. Solve the general linear equations,
8
>
>
>
>
>
>
>
>
>
>
>
>
<

f (x 1 ) =f (x 2 ) =

>
>
>
>
>
>
>
>
>
>
>
>
: f (x
m) =

x 1
x 2
x m

1
1

1

!T

x p

!T

x p

G

x 1

G

x 2

G

x m

..
.!

T

x p

1
1

1

x p
x p
x p

!
!

;

(B.16)!

to get G.
7. Solve the linear equation,

G

y

1

!

=0 :

(B.17)

Denote the solution by y .
8. Let x  =y  + x p, calculate f (x  ), and let P = P
9. If s = x max k x i i ;x j 2P

go to Step 2.

Sfx g.

x j k smaller than a given small real number, stop; otherwise

158

APPENDIX B. ON THE OPTIMIZATION ALGORITHM

In the above algorithm, if k = 0 or Step 2 to Step 5 are omitted, the algorithm becomes a pure interplanative search method. This may speed up the convergence, but may
also cause problem on robustness. On the other hand, if Step 5 to Step 8 are omitted,
the algorithm because a pure random search method.

Acknowledgement

I would like to express my gratitude to Prof. Yosihiko Ogata for his supervision, supportand encouragement throughout the last three years during my researches and studies for
the PhD degree. I also thank Prof. David Vere-Jones, from School of Mathematical andComputer Sciences, Victoria University of Wellington, for his guidance and encouragementin the research eld of statistical seismology.
I thank Prof Li Ma, from Center for Analysis and Prediction, China Seismological
Bureau, for continuous encouragement and support. Prof. Masaharu Tanemura, Prof
Yoshiyasu Tamura, both from Institute of Statistical Mathematics, Prof. Shigeru Mase,from Tokyo Institute of Technology, and Prof. Kunihiko Shimazaki, from Tokyo University, carefully read the manuscript of this thesis and gave many helpful comments. The
assistance of Dr. Ray Brownrigg, Dr. David Harte, from Victoria University of Wellington, Prof. Makoto Taiji and Dr. Koichi Katsura, from Institute of Statistical Mathematics
was helpful. All of these are gratefully acknowledged here.
I thank the Ministry of Education, Science, Sport and Culture of Japan, for providingme the Mobusho scholarship in supporting my study program throughout last 3 years.
The researches associated with this thesis was mainly supported by grant-in-aid 11680334
for the Scienti c Research from the Ministry of Education, Science, Sport and Culture ofJapan, and partially supported by some funds from China Seismological Bureau, Asian
2000 from New Zealand, and a fund from Prof. Yoshiyasu Tamura. These nancial
supports are gratefully acknowledge here.
159

160

Acknowledgement

Finally, I express my thanks to my family for their support and encouragement.

Bibliography

Akaike H., (1974). A new look at the statistical model identi cation. IEEE Transaction
on Automatic Control, AC-19, 716-723.

Aki K., (1981). A probabilistic synthesis of precursory phenomena. Earthquake Prediction: An International Review (eds. D.W. Simpson and P. G. Richards), 566-574,

A.G.U., Washington D.C..
Andersen P., Borgan O., Gill R. and Keiding N. (1993).Statistical Models Based on
Counting Processes, Springer, New York.

Bebbington M. and Harte D., (2001).

On the statistics of the linked stress releasemodel. Journal of Applied Probabilities, 38A, 176 { 187.
Chen Z., Du X. and et al, (1998). Electromagnetic radiation and earthquakes (in Chinese)(Seismological Press, Beijing).
Choi E. and Hall P., (1999). Nonparametric approach to analysis of space-time data on
earthquake occurrences. Journal of Computational and Graphical Statistics, 8, 733{ 748.
Cress G. O., Braoly B. T. and Rowell G. A., (1987). Sources of electromagnetic radiationfrom fracture of rock samples in the laboratory. Geophysical Research Letters, 14(4),331 { 334.
Czisar I., (1969). On generalized entropy. Stud. Sci. Math. Hungar, 4, 401 { 419.
161

162

BIBLIOGRAPHY

Daley D. J. and Vere-Jones D., (1972). A summary of the theory of point process. In:
Stochastic Point Processes: Statistical Analysis, Theory and Applications, edited by

P. A. W. Lewis. Wiley Inter-science, a division of John Wiley & Sons, Inc., NewYork.
Davis S. D. and Frohlich C., (1991). Single-link cluster analysis, synthetic earthquake
catalogs, and aftershock identi cation. Geophysical Journal International, 104, 289{ 306.
Diggle P., (1983). Statistical Analysis of Spatial Point Patterns, Academic Press, London.
Fritz J., (1973). An approach to the entropy of point processes. Perio. Math. Hungar.,

3, 73 { 83.
Frohlich C. and Davis S. D., (1990) Single-link cluster analysis as a method to evaluate spatial and temporal properties of earthquake catalogs. Geophysical JournalInternational, 100, 19 { 32.

Gardner J. and Knopo L. (1974). Is the sequence of earthquakes in southern California, with aftershock removed, Poissonian? Bulletin of the Seismological Society ofAmerica, 64, 1363 { 1367.

Gershengon N. I., Gokhberg M. B., Karakin A. V., Petviashvili N. V. and Rykunov A.
L., (1989). Modelling the connection between earthquake preparation processes andcrustal electromagnetic emission. Physics of the Earth and Planetary Interiors, 57,
129 { 138.
Guan H. and Liu G., (1995). Electromagnetic radiation anomalies before moderate andstrong earthquakes.Acta Seismologica Sinica, 9(2), 289 { 299.
Guan H. and Chen Z., (1996). Relationship between electric-magnetic radiations andearthquakes. Acta Seismologica Sinica, 17(2), 237 { 246.

BIBLIOGRAPHY

163

Guan H., Chen Z. and Xiao W., (2002). Observation on the ULF ground electric eld.
Submitted to Pure and Applied Geophysics.
Harte D., Li. D., Vreede M. and Vere-Jones D., (2002a). Quantifying the M8 prediction
algorithm I: reduction to a single critical variable and stability results. Submitted.
Harte D., Li. D., Vreede M. and Vere-Jones D., (2002b). Quantifying the M8 prediction
algorithm II: model forecast and evaluation. Submitted.
Hawkes A. G., (1971). Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58, 83 { 90.
Helmstetter A. and Sornette D., (2002). Predictability in the ETAS model of interacting
triggered seismicity. Submitted to Journal of Geophysical Research.
Isham V. and Westcott M., (1979). A self-correct point process. Stochastic processes
and Their Applications, 8, 335 { 347.

Kagan Y., (1991). Likelihood analysis of earthquake catalogues. Journal of GeophysicalResearch, 106, Ser. B7, 135 { 148.

Kanwal Ram P., (1997). Linear Integral Equations: Theory and Technique. Birkhauser,
Boston.
Karr A. F., (1991). Point Processes and Their Statistical Inference. Dekker, New York.
Kellis-Borok V. I. and Kossobokov V. I., (1986). Time of increased probability for thegreat earthquakes of the world. Computational Seismology, 19, 45 { 58.
Kellis-Borok V. I. and Kossobokov, V. G., (1990). Premonitory activation of seismic
ow: Algorithm M8. Physics of the Earth and Plant. Interior, 61, 73 { 83.
Kossobokov, V. G. and Kellis-Borok, V. I. (1990). Localization of intermediate-term
earthquake prediction. Journal of Geophysical Research, 61, 73 { 83.

164

BIBLIOGRAPHY

Lewis P. A. W. and Shedler E., (1979). Simulation of non-homogeneous Poisson processesby thinning. Naval Res. Logistics Quart., 26, 403 { 413.
Li Y. and Xia Y., (1989) Experiments on electromagnetic wave radiation from rockfracture and earthquake prediction. In Collection of research papers on the observation and experiments on the electro-magnetic waves before earthquakes (ed China

Seismological Bureau, in Chinese), Seismological Press, Beijing.
Liptzer R. and Shiryev M., (1978). Statistics of Random Processes, Springer-Verlag,
Berlin.
Liu J., Chen Y., Shi Y. and Vere-Jones D, (1999). Couple stress release model for timedependent seismicity. Pure and Applied Geophysics, 155, 649 { 667.
Liu J. (2001). Monte Carlo Strategies in Scienti c Computing. Springer, New York.
Chapter 2.
Lu C., Harte D. and Bebbington M., (1999). A linked stress release model for historical Japanese earthquakes: Coupling among major seismic regions. Earth Planets
Space., 51, 907 { 916.

Ma L. and Vere-Jones D., (1997). Applications of M8 and Lin-Lin algorithms of New
Zealand earthquake data. New Zealand Journal of Geology and Geophysics, 40, 77{ 89.
Musmeci F. and Vere-Jones D., (1986). Orthogonal polynomial estimates for point
process. Report RT/PAS/86/2, Itialian Commission for Nuclear and Alternative
Energy Sources (ENEA).

Musmeci F. and Vere-Jones D., (1986). A variable-grid algorithm for smoothing clustereddata. Biometrics, 42, 483 { 494.

BIBLIOGRAPHY

165

Musmeci F. and Vere-Jones D., (1992). A space-time clustering model for historicalearthquakes. Annals of the Institute of Statistical Mathematics, 44, 1 { 11.
von Neumann J., (1951). Various techniques used in connection with random digits/
National Bureau of Standards Applied Mathematics Series, 12, 634 { -637.

Nishizawa O., Lei. X. and Nagao T., (1994). Hazard function analysis of seismo-electric
signals in Greece. In Electromagnetic Phenomena Related to Earthquake Prediction(edited by M. Hayakaw and Y. Fujinawa), Terra Scienti c Publishing Company,
Tokyo.
Ogata Y., (1981). On Lewis' simulation method for point processes. IEEE Translations
on Information Theory, IT-27, 23 { 31.

Ogata Y. and Akaike H., (1982). On linear intensity models for mixed doubly stochastic
Poisson and self-exciting point processes. Journal of the Royal Statistical Society,
Series B, 44, 102 { 107.

Ogata Y., Akaike H. and Katsura K., (1982). The application of linear intensity models
to the investigation of causal relations between a point process and another stochastic process. Annals of the Institute of Statistical Mathematics, 34, 373 { 387.
Ogat, Y., (1983). Likelihood analysis of point processes and its application to seismological data. Bulletin of the International Statistical Institute, 50(2), 943 { 961.
Ogata Y. and Vere-Jones D., (1984). Inference for earthquake models: a self-correcting
model. Stochastic Processes and Their Applications, 17, 337 { 347.
Ogata Y., (1988). Statistical models for earthquake occurrences and residual analysis
for point processes. Journal of American Statistical Association, 83, 9 { 27.
Ogata Y. and Katsura K., (1988). Likelihood analysis of spatial inhomogeneity formarked point patterns. Annals of the Institute of Statistical Mathematics, 40: 20 {
39.

166

BIBLIOGRAPHY

Ogata Y., (1989). Statistical model for standard seismicity and detection of anomalies
by residual analysis for point process. Tectonophysics, 169, 1 { 16.
Ogata Y., (1992). Detection of precursory relative quiescence before great earthquakes
through a statistical model. Journal of Geophysical Research, 97, 19845 { 19871.
Ogata Y., Utsu T. and Katsura K., (1995). Statistical discrimination of foreshocks fromother earthquake clusters. International Journal of Geophysics, 127, 17 { 30.
Ogata Y., (1998). Space-time point-process models for earthquake occurrences. Annalsof the Institute of Statistical Mathematics, 50(2), 379 { 402.

Ogata Y., (1999), Seismicity Analysis through point-process modelling: A Review. Pureand Applied Geophysics, 155, 471 { 507.

Ogata Y., (2001). Modelling of heterogeneous space-time seismic activity and its residual
analysis. Research Memorandum, No. 808, The Institute of Statistical Mathematics,
Tokyo.
Ogata Y. and Zhuang J., (2001). Statistical examination of anomalies for the precursor to earthquakes and the multi-element prediction formula: Hazard rate changesof strong earthquakes (M  4:0) around Beijing area based on the ultra-low frequency electric observation (1982 { 1997). Report of the Coordinating Committee
for Earthquake Prediction, 66, 562 { 570. Geophysical Survey Institute, Tsukuba,

Japan.
Oike K. and Ogawa T., (1994). Relationship between shallow earthquakes and electromagnetic noise. In Electromagnetic Phenomena Related to Earthquake Prediction,
Terra Scienti c Publishing Co., Tokyo.
Omori F., (1894). On after-shocks of earthquakes. J. Coll. Sci. Imp. Univ. Tokyo, 7,111-200.

BIBLIOGRAPHY

167

Rathbun S. L., (1993). Modeling marked spatio-temporal point patterns. Bulletin of theInternational Statistical Institute, 55, Book 2, 379 { 396.

Reiss R.-D., (1993). A Course on Point Processes. Springer-Verlag, New York.
Resenberg P., (1985). Second-order moment of Central California Seismicity, 1969 {1982. Journal of Geophysical Research. 90, Ser. B7, 5479 { 5495.
Rubin I., (1972). Regular point processes and their detection. IEEE Translation ofInformation Theory, IT-18, 547 { 557.

Schoenberg F., (1999). Transforming spatial point processes into Poisson processes,
Stochastic Processes and Their Applications, 81, 155 { 164.

Silverman B. W., (1986). Density Estimation for Statistics and Data Analysis (1st ed.),
Chapman & Hall, London.
Tranter C., (1968). Bessel functions with some physical applications. The English Universities Press Ltd.
Utsu T. and Seki A., (1955). Relation between the area of aftershock region and theenergy of the main shock (in Japanese with English summary). Zisin, Series 2, 7,
233 { 240.
Utsu T., (1957). Magnitudes of earthquakes and occurrence of their aftershocks (inJapanese with English summary). Zisin, Series 2, 10, 35 { 45.
Utsu T., (1962). On the nature of three Alaskan aftershock sequences, Bulletin of Seismological Society of America. 52, 279 { 297.

Utsu T., (1969). Aftershock and earthquake statistics (I): some parameters which characterize an aftershock sequence and their interrelations. Journal of the Faculty of
Science, Hokkaido University, 3, Ser. VII (Geophysics), 129 { 195.

168

BIBLIOGRAPHY

Utsu T., (1977). Probabilities in earthquake prediction. Zisin (in Japanese with English
abstract). 30(2), 179 { 185.
Utsu T., (1982). Probabilities in earthquake prediction (the second paper). Bulletin of
the Earthquake Research Institute, 57, 499 { -524.

Utsu T., Ogata Y. and Matsu'ura R. S., (1995). The centenary of the Omori formulafor a decay law of aftershock activity, J. Phys. Earth, 43, 1 { 33.
Utsu T. and Ogata Y., (1997). Statistical analysis of seismicity. In Algorithm for Earthquake Statistics and Prediction, IAPSPEI Software Library, International Associa-

tion of Seismologigy and Physics of the Earth's Interior, 6, 13 { 94.
Varotsos P. and Kulhanek O. (eds.), (1993). Measurements and Theoretical Models of
the Earth's Electric Field Variations related to Earthquakes, Tectonophysics, 224,

No. 1/3.
Varotsos P., Eftaxias K., and et al., (1996). Recent Earthquake Prediction Resultsin Greece Based on the Observation of Seismic Electric Signals. Acta Geophysica
Polonica, XLIV, 301-327.

Vere-Jones D., (1978). Earthquake prediction { a statistician's view. Journal of Physics
of the Earth., 26, 129 { 146.

Vere-Jones D. and Ogata Y., (1984). On the moments of a self-correcting process.
Journal of Applied Probabilities, 21, 335 { 342.

Vere-Jones D., (1988). A point process analysis of historical earthquakes from North
China. Earthquake research in China, 2(2), 165 { 181.
Vere-Jones D., (1992). Statistical methods for the description and display of earthquake
catalogues. In: Statistics in the environmental and earth sciences, edited by A.

169

BIBLIOGRAPHY

Walden, P. Guttorp. Edward Arnold, a division of Hodder & Stoughton, London,
220 { 236.
Vere-Jones D., (1995).

Forecasting earthquakes and earthquake risk. International

Journal of Forecasting, 11, 503 { 538.

Vere-Jones D., (1999).

Probability and information gain for earthquake forecasting.

Computational Seismology, 30, 248-263.

Yoshino T. and Tomizawa I., (1990). Observation Results of Low Frequency Electromagnetic Emissions as Precursors of Volcanic Eruption at Mt. Mihara in November,1986. Journal of Geomagnetism and Geoelectricity, 42, 225 { 235.
Zheng X. and Vere-Jones D., (1991). Application of stress release models to historicalearthquakes from North China. Pure and Applied Geophysics, 135, 559 { 576.
Zheng X. and Vere-Jones D., (1994). Further application of stress release models tohistorical earthquake data. Tectonophysics, 229, 101 { 121.
Zhuang J., Ogata Y. and Vere-Jones D., (2002). Stochastic declustering of space-timeearthquake occurrences. Journal of the American Statistical Association, 97, 369 {
380.

170

BIBLIOGRAPHYThe Author's Publications andSubmissions During 2000 to 2003

1. Zhuang J. and Yin X., (2000). Random distribution of the loading and unloading
response ratio under the assumptions of Poisson models. Earthquake research in
China, 14(2), 38 { 48.

2. Zhuang J., (2000). Statistical modelling of seismicity patterns before and after the1990 Oct 5 Cape Palliser earthquake, New Zealand. New Zealand Journal of Geology
and Geophysics, 43, 447 { 460.

3. Ma L. and Zhuang J., (2001). Relative quiescence within the Jiashi swarm in Xinjiang, China: an application of the ETAS point process model. Journal of Applied
Probability, 38(A): 213 { 221.

4. Zhuang J., Ogata Y. and Vere-Jones D., (2002). Stochastic declustering of spacetime earthquake occurrences. Journal of the American Statistical Association, 97,458, 369{380.
5. Ogata Y. and Zhuang J., (2001). Statistical examination of anomalies for the
precursor to earthquakes, and the multi-element prediction formula: Hazard ratechanges of strong earthquakes M34.0 around Beijing area based on the ultra-low
frequency electric observation (1982-1997). In Report of the Coordinating Committee for Earthquake Prediction (In Japanese), 66, 562 { 570. Geophysical Survey

Institute, Tsukuba, Japan.
171

172

BIBLIOGRAPHY

6. Zhuang J. and Ma L., (2000). Main shock and after shock { from the Omori law tothe epidemic aftershock model (In Chinese with English abstract). Recent Development in World Seismology. 2000(5), 12 { 18.

7. Zhuang J., Ogata Y., Vere-Jones D., Ma L. and Guan H., (2002). Statistical con rmation of a relationship between excitation of the ultra-low frequency ground electric signals and magnitude M  4 earthquakes in a 300 km radius region around
Beijing. Submitted to Pure and Applied Geophysics.