Ayten Yiğiter

*Hacettepe University, Faculty of Science, Department of Statistics Beytepe-Ankara Turkey* 

#### **1. Introduction**

22 Earthquake Research and Analysis / Book 5

24 Earthquake Research and Analysis – Statistical Studies, Observations and Planning

[30] G. Yakovlev, D. L. Turcotte, J. B. Rundle, and P. B. Rundle, "Simulation-Based Distributions of Earthquake Recurrence Times on the San Andreas Fault System",

[31] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in C*

[32] T. Seno, S. Seth, and E. G. Alice, "A model for the motion of the Philippine Sea plate consistent with NUVEL-1 and geological data", *Journal of Geophysics Research*, vol 98, pp.

[33] C. M. R. Fowler, *The Solid Earth: An Introduction to Global Geophysics*, Cambridge

[34] L. Ruff and H. Kanamori, "Seismicity and the subduction process", *Physics of the Earth*

[35] T. Wong, R. H. C. Wong, K. T. Chau, and, C. A. Tang, "Microcrack statistics. Weibull distribution and micromechanical modeling of compressive failure in rock", *Mechanics*

[36] A. Ghosh, "A FORTRAN program for fitting Weibull distribution and generating

[37] D. L. Turcotte, W. I. Newman, R. Shcherbakov, "Micro and macroscopic models of rock

*Bulletin of the Seismological Society of America*, vol. 96, pp. 1995-2007, 2006.

*2nd edition*, Cambridge University Press, Cambridge, 1995.

samples", *Computational Geosciences*, vol. 25, pp. 729-738, 1999.

fracture", *Geophysical Journal International*, vol. 152, pp. 718-728, 2003.

17941-17948, 1993.

University Press, New York, 1990.

*of Materials*. vol. 38, pp. 664-681, 2006.

[38] Y. Aizawa and T. Hasumi, in preparation (2011).

*and Planetary Interiors*, vol. 23, pp. 240-252, 1980.

Earthquake forecasts are very important in human life in terms of estimating hazard and managing emergency systems. Defining of earthquake characteristics plays an important role in these forecasts. Of these characteristics, one is the frequency distribution of earthquakes and and the other is the magnitude distribution of the earthquakes. Each statistical distribution has many parameters describing the actual distribution.

There are various statistical distributions used to model the earthquake numbers. As is well known, these are binomial, Poisson, geometric and negative binomial distributions. It is generally assumed that earthquake occurrences are well described by the Poisson distribution because of its certain characteristics (for some details, see Kagan, 2010; Leonard & Papasouliotis, 2001). In their study, Rydelek & Sacks (1989) used the Poisson distribution of earthquakes at any magnitude level. The Poisson distribution is generally used for earthquakes of a large magnitude, and the earthquake occurrences with time/space can be modeled with the Poisson process in which, as is known, the Poisson distribution is one that counts the events that have occurred over a certain period of time.

There is a significant amount of research on the change point as applied to earthquake data. Amorese (2007) used a nonparametric method for the detection of change points in frequency magnitude distributions. Yigiter & İnal (2010) used earthquake data for their method developed for the estimator of the change point in Poisson process. Aktaş et al. (2009) investigated a change point in Turkish earthquake data. Rotondi & Garavaglia (2002) applied the hierarchical Bayesian method for the change point in data, taken from the Italian NT4.1.1 catalogue.

Recently, much research in the literature has focused on whether there is an increase in the frequency of earthquake occurrences. It is further suggested that any increase in the frequency of earthquakes, in some aspects, is due to climate change in the world. There is considerable debate on whether climate change really does increase the frequency of natural disasters such as earthquakes and volcano eruptions. In many studies, it is emphasized that there is serious concern about impact of climate change on the frequencies of hazardous events (Peduzzi, 2005; Lindsey, 2007; Mandewille, 2007 etc.). In Peduzzi's study (2005), there are some indicators about increasing number of the earthquakes especially affecting human settlements, and it is also reported that there is an increase in the percentage of earthquakes affecting human settlements from 1980 onwards. The change point analysis can be used to study the increase or decrease in the frequency of the earthquake occurrences.

Change Point Analysis in Earthquake Data 27

 Nt t i i 1

> 

2

in which {N ,t 0} <sup>t</sup> is a Poisson process, and {Y ,i 1,2, } <sup>i</sup> are independent identically distributed (i.i.d.) random variables. The compound Poisson distribution has stationary

> t t

**3.1 Change point estimation with the maximum likelihood method (frequentist** 

t t

The Maximum likelihood method and the Bayesian method are basic methods in statistical

Let X ,X , ,X 12 n be a sequence of the random variables such that Xi (i 1, , ) has probability density function <sup>1</sup> f(x, ) and Xi (i 1, ,n) has probability density function <sup>2</sup> f(x, ) where the change point is unknown discrete parameter and the parameters <sup>1</sup> , 2 can be assumed to be either known or unknown. The single change point model in the

X ,X , ,X ~ f(x, )

 1 2 1 12 n 2

Under model Eq.(5), for the observed values of the sequence of the random variables, the

12 1 n 12 i1 i2

 v n 1 2 i 1 i 2 i 1 iv1

1 2 i1 i2

When the parameters, <sup>1</sup> , 2 are unknown, for any fixed values of change point , the maximum likelihood estimator of the parameters <sup>1</sup> , 2 are found to be the derivative of

L( , , x , ,x ) L( , , ) f(x , ) f(x , ) (6)


X ,X , ,X ~ f(x, ). (5)

nL( , , ) nf(x , ) nf(x , ) (7)

(8)

nf(x , ) in Eq.(7) then Eq.(7) can be re-written as follows:

i1 iv1

v n

i 2

v

i 1 nL( , , ) nf(x , ) nf(x , )

E(X ) E N E Y V(X ) E N E(Y ).

independent increments, and the mean and variance of the Xt,

See Parzen (1962) for some details about Poisson processes.

change point analyses for point estimation and hypothesis testing.

**3. Methods for Estimating Change Point** 

sequence of the random variables is written:

and the logarithm of the likelihood function is:

 v

i 1

**method)** 

likelihood function is

After adding the expression

Eq.(8) 1 and 2 respectively:

X Y (3)

(4)

The change point analysis is a very useful statistical tool to detect an abrupt or a structural change in the characteristic of the data observed sequentially or chronologically. Many different statistical methods are available to detect the change point in the distribution of a sequence of random variables (Smith, 1975; Hinkley, 1970; Boudjellaba et al., 2001 in many others). Multiple change points have also been investigated for many years (Hendersen & Matthews, 1993; Yao, 1993; Chen & Gupta, 1997, among others.).

In this study, we aimed to find an evidence of increased or decreased in world seismic activities after 1901. Earthquake frequencies are modeled by the Poisson distribution as the most commonly used discrete distribution. We modeled the earthquakes in the world with the Poisson process since the number of earthquakes is counted with Poisson distribution. We investigated any abrupt change point(s) in parameters of the model using the frequentist (maximum likelihood) and Bayesian method. When the magnitude of the earthquakes is taken into account in the process, it is then called the compound Poisson process. We also investigated a change in the magnitudes of the world earthquakes. For this purpose, Poisson process and compound Poisson process are introduces in Section 2. The frequentist and Bayesian methods used for change point estimates are explained in Section 3. Worldwide earthquake data and change point analysis of this data are given Section 4.

#### **2. Poisson process**

The Poisson process has an important place in stochastic processes. It is a Markov chain with a continuous parameter.

A stochastic process {N ,t 0} <sup>t</sup> is said to be a homogeneous Poisson process if


In the homogeneous Poisson Process, events occur independently throughout time. The number of events occurred time interval (0, t] has Poisson distribution with parameter t,

$$\mathbf{P(N\_t = i)} = \mathbf{e^{\cdot \lambda t}} \frac{\left(\lambda t\right)^i}{i!}, \qquad i = 0, 1, \ldots \tag{1}$$
 
$$= 0, \qquad \text{otherwise}$$

where is called the occurrence rate of the event in unit time. The occurrence rate is assumed to be a constant throughout the process.

Let S ,S ,S , <sup>012</sup> be occurrence times of the events where S 0 <sup>0</sup> and TSS <sup>110</sup> , T S S, 2 21 be time intervals between the events. T ,T , 1 2 are independent identically exponential random variables with parameter . The distribution of the random variables T ,T , 1 2 is:

$$\begin{aligned} \mathbf{f}(\mathbf{t}) &= \lambda \, \mathbf{e}^{\cdot \lambda \mathbf{t}}, \quad \mathbf{t} \ge \mathbf{0} \\ &= \mathbf{0}, \qquad \mathbf{t} < \mathbf{0} \end{aligned} \tag{2}$$

#### **2.1 Compound poisson process**

A stochastic process {X ,t 0} <sup>t</sup> is said to be a compound Poisson process if it can be presented, for t 0, by

The change point analysis is a very useful statistical tool to detect an abrupt or a structural change in the characteristic of the data observed sequentially or chronologically. Many different statistical methods are available to detect the change point in the distribution of a sequence of random variables (Smith, 1975; Hinkley, 1970; Boudjellaba et al., 2001 in many others). Multiple change points have also been investigated for many years (Hendersen &

In this study, we aimed to find an evidence of increased or decreased in world seismic activities after 1901. Earthquake frequencies are modeled by the Poisson distribution as the most commonly used discrete distribution. We modeled the earthquakes in the world with the Poisson process since the number of earthquakes is counted with Poisson distribution. We investigated any abrupt change point(s) in parameters of the model using the frequentist (maximum likelihood) and Bayesian method. When the magnitude of the earthquakes is taken into account in the process, it is then called the compound Poisson process. We also investigated a change in the magnitudes of the world earthquakes. For this purpose, Poisson process and compound Poisson process are introduces in Section 2. The frequentist and Bayesian methods used for change point estimates are explained in Section 3. Worldwide

The Poisson process has an important place in stochastic processes. It is a Markov chain

ii. for any times s and t such that s<t the number of events occurred in time interval (s, t)

In the homogeneous Poisson Process, events occur independently throughout time. The number of events occurred time interval (0, t] has Poisson distribution with parameter t,


where is called the occurrence rate of the event in unit time. The occurrence rate is

Let S ,S ,S , <sup>012</sup> be occurrence times of the events where S 0 <sup>0</sup> and TSS <sup>110</sup> , T S S, 2 21 be time intervals between the events. T ,T , 1 2 are independent identically exponential random variables with parameter . The distribution of the random variables

> - <sup>t</sup> f(t) e , t 0 0, t 0

A stochastic process {X ,t 0} <sup>t</sup> is said to be a compound Poisson process if it can be

( t) P(N i) e , i 0,1, i!

0, otherwise

(1)

(2)

Matthews, 1993; Yao, 1993; Chen & Gupta, 1997, among others.).

earthquake data and change point analysis of this data are given Section 4.

A stochastic process {N ,t 0} <sup>t</sup> is said to be a homogeneous Poisson process if

i. {N ,t 0} <sup>t</sup> has stationary independent increments

has Poisson distribution with parameter (t-s).

assumed to be a constant throughout the process.

**2.1 Compound poisson process** 

presented, for t 0, by

t

**2. Poisson process** 

T ,T , 1 2 is:

with a continuous parameter.

$$\mathbf{X}\_{\mathbf{t}} = \sum\_{\mathbf{i}=1}^{N\_{\mathbf{t}}} \mathbf{Y}\_{\mathbf{i}} \tag{3}$$

in which {N ,t 0} <sup>t</sup> is a Poisson process, and {Y ,i 1,2, } <sup>i</sup> are independent identically distributed (i.i.d.) random variables. The compound Poisson distribution has stationary independent increments, and the mean and variance of the Xt,

$$\begin{aligned} \mathrm{E(X\_t)} &= \mathrm{E(N\_t)E(Y)}\\ \mathrm{V(X\_t)} &= \mathrm{E(N\_t)E(Y^2)}. \end{aligned} \tag{4}$$

See Parzen (1962) for some details about Poisson processes.

#### **3. Methods for Estimating Change Point**

The Maximum likelihood method and the Bayesian method are basic methods in statistical change point analyses for point estimation and hypothesis testing.

#### **3.1 Change point estimation with the maximum likelihood method (frequentist method)**

Let X ,X , ,X 12 n be a sequence of the random variables such that Xi (i 1, , ) has probability density function <sup>1</sup> f(x, ) and Xi (i 1, ,n) has probability density function <sup>2</sup> f(x, ) where the change point is unknown discrete parameter and the parameters <sup>1</sup> , 2 can be assumed to be either known or unknown. The single change point model in the sequence of the random variables is written:

$$\begin{aligned} \mathbf{X}\_{1}, \mathbf{X}\_{2}, \dots, \mathbf{X}\_{\mathbf{v}} &\sim \mathbf{f}(\mathbf{x}, \boldsymbol{\theta}\_{1}) \\ \mathbf{X}\_{\mathbf{v}+1}, \mathbf{X}\_{\mathbf{v}+2}, \dots, \mathbf{X}\_{\mathbf{n}} &\sim \mathbf{f}(\mathbf{x}, \boldsymbol{\theta}\_{2}). \end{aligned} \tag{5}$$

Under model Eq.(5), for the observed values of the sequence of the random variables, the likelihood function is

$$\mathbf{L}(\boldsymbol{\theta}\_{1}, \boldsymbol{\theta}\_{2}, \mathbf{v} \mid \mathbf{x}\_{1}, \dots, \mathbf{x}\_{n}) = \mathbf{L}(\boldsymbol{\theta}\_{1}, \boldsymbol{\theta}\_{2}, \mathbf{v}) = \prod\_{i=1}^{\mathbf{v}} \mathbf{f}(\mathbf{x}\_{i}, \boldsymbol{\theta}\_{1}) \prod\_{i=\mathbf{v}+1}^{n} \mathbf{f}(\mathbf{x}\_{i}, \boldsymbol{\theta}\_{2}) \tag{6}$$

and the logarithm of the likelihood function is:

$$\ell \ln \text{L} \{ \theta\_1, \theta\_2, \mathbf{v} \} = \sum\_{i=1}^{\text{v}} \ell \text{nf} \{ \mathbf{x}\_i, \theta\_1 \} + \sum\_{i=\text{v}+1}^{\text{n}} \ell \text{nf} \{ \mathbf{x}\_i, \theta\_2 \} \tag{7}$$

After adding the expression v i 2 i 1 nf(x , ) in Eq.(7) then Eq.(7) can be re-written as follows:

$$\ell \ln \text{L} (\theta\_1, \theta\_2, \mathbf{v}) \propto \sum\_{i=1}^{\text{v}} \left[ \ell \text{nf} (\mathbf{x}\_i, \theta\_1) \cdot \ell \text{nf} (\mathbf{x}\_i, \theta\_2) \right] \tag{8}$$

When the parameters, <sup>1</sup> , 2 are unknown, for any fixed values of change point , the maximum likelihood estimator of the parameters <sup>1</sup> , 2 are found to be the derivative of Eq.(8) 1 and 2 respectively:

$$\frac{\partial \ell \ln \text{L}(\Theta\_1, \Theta\_2; \mathbf{v})}{\partial \theta\_1} = 0 \,, \tag{9}$$

Change Point Analysis in Earthquake Data 29

To detect a change point in earthquakes, the relevant data is taken from the website of the U.S. Geological Survey (2011), consisting of 819 earthquakes worldwide of magnitude 4.0 or above covering the period from 3-March-1901 until 11-March-2011. It is assumed that the earthquake occurrences follow the homogeneous Poisson process. The earthquakes are observed at times S ,S , ,S 1 2 818 in the continuous time interval (0,T]. For the earthquakes data, the starting point S0 is in 3-March-1901 and the end point S818=T is in 11-March-2011. T ,T , ,T 1 2 818 are exponentially distributed random variables as given in Eq.(2). Under the assumption that is equal to the occurrence time of an event, the likelihood function can be written under the change point model in the sequence of exponentially distributed random


Where shows change point and is a continuous parameter defined time interval (0,T], <sup>2</sup> is the occurrence rate in the time interval ( , T] and N is the number of events that occurred in the time interval (0, ] and T is the sum of the time interval between the

After the derivation of the log likelihood function given by Eq.(15) for 1 and <sup>2</sup> , the

For the possible values of , s ,s , ,s 12 n , and the variable N takes the values 1,2, ,n , and so the estimations ( , ),( , ), ,( , ) are calculated from Eq.(16) and then ˆˆ ˆˆ ˆˆ

The scatter plot of the time intervals between the earthquakes is shown in Figure 1. From Figure 1, it is clearly seen that there is a change (or at least one change) in the occurrence

When we assume a single change point in the data and the method explained above is

As is shown in Table 1, the occurrence rate of the earthquake in unit time (in day or in year) increased considerably after February 3, 2002, approximately ten times of the occurrence

From the log likelihood function of the change point in Figure 2, it can be seen that there are at least two change points in the data. There many methods are used to detect multiple change points in the data. One of the basic methods used for this purpose is binary

1 2 - λ λ - 

11 21 12 22 1n 2n

these estimations are substituted into Eq.(15). One of 1i 2i i

 818 i i 1

maximum likelihood estimations of 1 and <sup>2</sup> are obtained respectively:

2

n N (T ) e . (14)

T t (Raftery & Akman, 1986; Akman & Raftery, 1986).

ˆ ˆ N nN , <sup>T</sup> (16)

ˆ ˆ , , . Hence ˆ 1 2 ˆ ˆ , and ˆ are the maximum likelihood

ˆ ˆ ( , ,s ) , i 1,2 , ,n , gives the

<sup>2</sup> -- - nL( , , ) N n( ) ( ) n n T 1 2 12 1 2 2 (15)

variables:

earthquake occurrences, that is

The logarithm of the likelihood function is

maximum of this function, say 1 2

applied this data and the results are given in Table 1.

rate of the earthquakes in unit time before February 3, 2002.

estimation for 1, 2 and .

rate of the earthquakes.

$$\frac{\partial \ell \ln \text{L}(\Theta\_1, \Theta\_2; \mathbf{v})}{\partial \Theta\_2} = \mathbf{0} \; . \tag{10}$$

After solving the equations system given in Eq.(9) and (10), the maximum likelihood estimator of the parameters <sup>1</sup> , <sup>2</sup> , <sup>1</sup> <sup>ˆ</sup> , <sup>2</sup> ˆ are obtained. The maximum likelihood estimate of the change point is:

$$\hat{\mathbf{v}} = \arg\max\_{\mathbf{k}=1,\ldots,n-1} \sum\_{i=1}^{k} \left[ \ell \text{nf}(\mathbf{x}\_{i}, \hat{\boldsymbol{\theta}}\_{1}) \cdot \ell \text{nf}(\mathbf{x}\_{i}, \hat{\boldsymbol{\theta}}\_{2}) \right] \tag{11}$$

#### **3.2 Change point estimation with the Bayesian method**

The Bayesian method differs from the frequentist method in that each parameter is assumed to be a random variable and each one has a probability function called prior distribution. The estimate of the unknown parameter is obtained by deriving a posterior distribution on the basis of the prior distributions and the likelihood function. The posterior distribution is obtained:

#### Posterior likelihood prior.

Under change point model given in Eq.(5), let p ( , ) be the joint prior distribution of 012 parameters <sup>1</sup> , 2 and let p( ) 0 be the prior distribution of change point . The likelihood function is given by Eq.(6) and so the joint posterior distribution is written as follows:

$$\mathbf{p}\_1(\theta\_1, \theta\_2, \mathbf{v}) \propto \mathcal{L}(\theta\_1, \theta\_2, \mathbf{v}) \mathbf{p}\_0(\theta\_1, \theta\_2 \mid \mathbf{v}) \mathbf{p}\_0(\mathbf{v}) \,. \tag{12}$$

Integrate Eq.(12) with respect to the parameters, the marginal posterior distribution of change point is proportional to

$$\mathbf{p}\_1(\mathbf{v}) \propto \int\_{\theta\_2 \mid \theta\_1} \mathbf{L}(\theta\_1, \theta\_2, \mathbf{v}) \mathbf{p}\_0(\theta\_1, \theta\_2 \mid \mathbf{v}) \mathbf{p}\_0(\mathbf{v}) \, \mathrm{d}\theta\_1 \mathrm{d}\theta\_2 \tag{13}$$

and the Bayesian estimate ˆ of the change point is found by maximizing the marginal posterior distribution given by Eq.(13). Assuming uniform priors, the joint posterior mode, which gives the maximum likelihood estimates, is at ˆ , <sup>1</sup> <sup>ˆ</sup> , <sup>2</sup> ˆ (Smith, 1975).

#### **4. Detecting change point in worldwide earthquake data and findings**

This study investigates whether there is a change point in worldwide earthquake activities such as the number of earthquake occurrences and their magnitude. At first, assuming that the occurrences of the earthquakes follow the homogeneous Poisson process, a change point is investigated in the occurrence rate of the earthquake in unit time; secondly, the magnitudes of the earthquake are assumed to be normally distributed random variables; a change point is investigated in the mean of the normally distributed random variables describing the earthquake magnitudes.

( , ;) <sup>0</sup> , (9)

( , ;) <sup>0</sup> . (10)

ˆ are obtained. The maximum likelihood estimate

 nL 1 2 1

 nL 1 2 2


The Bayesian method differs from the frequentist method in that each parameter is assumed to be a random variable and each one has a probability function called prior distribution. The estimate of the unknown parameter is obtained by deriving a posterior distribution on the basis of the prior distributions and the likelihood function. The posterior

Posterior likelihood prior. Under change point model given in Eq.(5), let p ( , ) be the joint prior distribution of 012 parameters <sup>1</sup> , 2 and let p( ) 0 be the prior distribution of change point . The likelihood function is given by Eq.(6) and so the joint posterior distribution is written as follows:

Integrate Eq.(12) with respect to the parameters, the marginal posterior distribution of

p ( ) L( , , )p ( , )p ( )d d <sup>1</sup> 12 012 0 1 2

and the Bayesian estimate ˆ of the change point is found by maximizing the marginal posterior distribution given by Eq.(13). Assuming uniform priors, the joint posterior mode,

This study investigates whether there is a change point in worldwide earthquake activities such as the number of earthquake occurrences and their magnitude. At first, assuming that the occurrences of the earthquakes follow the homogeneous Poisson process, a change point is investigated in the occurrence rate of the earthquake in unit time; secondly, the magnitudes of the earthquake are assumed to be normally distributed random variables; a change point is investigated in the mean of the normally distributed random variables

**4. Detecting change point in worldwide earthquake data and findings** 

ˆ ˆ ˆ arg max nf(x , ) nf(x , ) (11)

p ( , , ) L( , , )p ( , )p ( ) 112 12 012 0 . (12)

(13)

<sup>ˆ</sup> , <sup>2</sup>

ˆ (Smith, 1975).

i1 i2 k 1, n 1 i 1

<sup>ˆ</sup> , <sup>2</sup>

k


**3.2 Change point estimation with the Bayesian method** 

2 1

which gives the maximum likelihood estimates, is at ˆ , <sup>1</sup>

estimator of the parameters <sup>1</sup> , <sup>2</sup> , <sup>1</sup>

of the change point is:

distribution is obtained:

change point is proportional to

describing the earthquake magnitudes.

After solving the equations system given in Eq.(9) and (10), the maximum likelihood

To detect a change point in earthquakes, the relevant data is taken from the website of the U.S. Geological Survey (2011), consisting of 819 earthquakes worldwide of magnitude 4.0 or above covering the period from 3-March-1901 until 11-March-2011. It is assumed that the earthquake occurrences follow the homogeneous Poisson process. The earthquakes are observed at times S ,S , ,S 1 2 818 in the continuous time interval (0,T]. For the earthquakes data, the starting point S0 is in 3-March-1901 and the end point S818=T is in 11-March-2011. T ,T , ,T 1 2 818 are exponentially distributed random variables as given in Eq.(2). Under the assumption that is equal to the occurrence time of an event, the likelihood function can be written under the change point model in the sequence of exponentially distributed random variables:

$$\mathbf{L}(\boldsymbol{\lambda}\_1, \boldsymbol{\lambda}\_2, \boldsymbol{\tau}) = \boldsymbol{\lambda}\_1^{\text{N}\_{\boldsymbol{\tau}}} \mathbf{e}^{\text{-}\boldsymbol{\lambda}\_1 \cdot \boldsymbol{\tau}} \boldsymbol{\lambda}\_{\boldsymbol{\tau}^{\text{n-N}\_{\boldsymbol{\tau}}}}^{\text{n-N}\_{\boldsymbol{\tau}}} \mathbf{e}^{\text{-}\boldsymbol{\lambda}\_2(\boldsymbol{T} - \boldsymbol{\tau})}.\tag{14}$$

Where shows change point and is a continuous parameter defined time interval (0,T], <sup>2</sup> is the occurrence rate in the time interval ( , T] and N is the number of events that occurred in the time interval (0, ] and T is the sum of the time interval between the earthquake occurrences, that is 818 i i 1 T t (Raftery & Akman, 1986; Akman & Raftery, 1986). The logarithm of the likelihood function is

$$\text{tr}\,\text{l}\,\text{nL}(\lambda\_1, \lambda\_2, \tau) = \text{N}\_\tau \ell \text{n}(\lambda\_1/\lambda\_2) \text{ -}\,\text{\textdegree{c}}(\lambda\_1 \text{ -}\lambda\_2) + \text{n}\,\text{l}\,\text{n}\lambda\_2 \text{ -}\lambda\_2\text{T} \tag{15}$$

After the derivation of the log likelihood function given by Eq.(15) for 1 and <sup>2</sup> , the maximum likelihood estimations of 1 and <sup>2</sup> are obtained respectively:

$$
\hat{\lambda}\_1 = \frac{\mathbf{N}\_\tau}{\mathbf{\tau}}, \quad \hat{\lambda}\_2 = \frac{\mathbf{n} \cdot \mathbf{N}\_\tau}{\mathbf{T} \cdot \mathbf{\tau}} \tag{16}
$$

For the possible values of , s ,s , ,s 12 n , and the variable N takes the values 1,2, ,n , and so the estimations ( , ),( , ), ,( , ) are calculated from Eq.(16) and then ˆˆ ˆˆ ˆˆ 11 21 12 22 1n 2n these estimations are substituted into Eq.(15). One of 1i 2i i ˆ ˆ ( , ,s ) , i 1,2 , ,n , gives the maximum of this function, say 1 2 ˆ ˆ , , . Hence ˆ 1 2 ˆ ˆ , and ˆ are the maximum likelihood estimation for 1, 2 and .

The scatter plot of the time intervals between the earthquakes is shown in Figure 1. From Figure 1, it is clearly seen that there is a change (or at least one change) in the occurrence rate of the earthquakes.

When we assume a single change point in the data and the method explained above is applied this data and the results are given in Table 1.

As is shown in Table 1, the occurrence rate of the earthquake in unit time (in day or in year) increased considerably after February 3, 2002, approximately ten times of the occurrence rate of the earthquakes in unit time before February 3, 2002.

From the log likelihood function of the change point in Figure 2, it can be seen that there are at least two change points in the data. There many methods are used to detect multiple change points in the data. One of the basic methods used for this purpose is binary

Change Point Analysis in Earthquake Data 31

Fig. 2. The maximum point of log likelihood function given in Eq.(16) is shown as a red star

Change points and the occurrence rates of the earthquake throughout time axis 1ˆ 2ˆ 3ˆ 4ˆ 5ˆ 6ˆ

0.0073\* 0.0107\* 0.0136\* 0.0254\* 0.1161\* 0.1603\* 0.0284\*

Fig. 3. The maximum likelihood estimations of parameters for multiple change points in the

03-Feb-2002 N<sup>ˆ</sup> =429

<sup>ˆ</sup> <sup>5</sup>

16-Jul-2007 N<sup>ˆ</sup> =660

<sup>ˆ</sup> <sup>6</sup>

12-Jan-2010 N<sup>ˆ</sup> =806

ˆ

<sup>ˆ</sup> <sup>7</sup>

and the corresponding change point estimate is shown with a red circle.

10-May-1997 N<sup>ˆ</sup> =385

<sup>ˆ</sup> <sup>4</sup>

03-Nov-1928 N<sup>ˆ</sup> =74

<sup>ˆ</sup> <sup>2</sup>

earthquake data (\*in days).

1

0

09-Mar-1957 N<sup>ˆ</sup> =185

<sup>ˆ</sup> <sup>3</sup>

segmentation procedure (Chen & Gupta, 1997; Yang & Kuo, 2001). In the binary segmentation procedure, the data is divided into two homogeneous groups according to the estimated change point, and a change point is searched in each subdivided data until there is no change in the subdivided data. The results are given in Figure 3.

Fig. 1. The time intervals (days) between the earthquakes according to the earthquake occurrences.


Table 1. Estimated parameters for the earthquake data.

segmentation procedure (Chen & Gupta, 1997; Yang & Kuo, 2001). In the binary segmentation procedure, the data is divided into two homogeneous groups according to the estimated change point, and a change point is searched in each subdivided data until there

Fig. 1. The time intervals (days) between the earthquakes according to the earthquake

n = 818, T= 40185 days

Estimated change point ˆ N<sup>ˆ</sup> =429

4.1896804 (in years) 03-Feb-2002 0.1170628 (in days)

After the change point 2 ˆ

42.1426080 (in years)

occurrences.

Before the change point 1 ˆ

0.0116380 (in days)

Table 1. Estimated parameters for the earthquake data.

is no change in the subdivided data. The results are given in Figure 3.

Fig. 2. The maximum point of log likelihood function given in Eq.(16) is shown as a red star and the corresponding change point estimate is shown with a red circle.

Change points and the occurrence rates of the earthquake throughout time axis

Fig. 3. The maximum likelihood estimations of parameters for multiple change points in the earthquake data (\*in days).

Change Point Analysis in Earthquake Data 33

where the prior distributions are called uninformative priors. The joint posterior

112 12 012 0

p ( , , ) L( , , )p ( , )p ( )

With integrate Eq. (18) with respect to the parameters <sup>1</sup> , <sup>2</sup> , the marginal posterior

where (.) is gamma function. The marginal posterior distributions of <sup>1</sup> and 2 are not obtained analytically because of discontinuity in and the close form of the marginal

<sup>1</sup> (n N ) p( ) e <sup>d</sup>

<sup>1</sup> (N ) p( ) e <sup>d</sup>

where, for the possible values of , s ,s , ,s 12 n , and the variable N takes the values

Form Eq.(19), the Bayesian estimate ˆ of the change point in the earthquake data is found to be the date 03-February-2002 2002 which corresponds to the posterior mode. The posterior

 <sup>2</sup>

n N (T ) 1 2 2 N 0

T N 1 1 1 n N

0

T

T

Fig. 6. The posterior distribution of change point and the posterior mode.

 <sup>1</sup>

 <sup>1</sup> N nN

1 2

<sup>T</sup> (T ) ,

T

(N ) (n N ) <sup>1</sup> p() (T ) <sup>T</sup> (19)

(18)

1


N nN (T )

e e

1 2

posterior distributions of <sup>1</sup> and 2 are respectively proportional to

distribution of the parameters is can be written as:

distribution of change point is proportional to

distribution of the change point is given in Figure 6.

1,2, ,n .

Fig. 4. Multiple change points in the earthquake data.

Estimated change points are shown on the scatter plot of the time intervals between the earthquakes.

Fig. 5. The estimated occurrence rate of earthquakes in order of change points.

As can be seen in Figure 4 and Figure 5, the occurrence rate of the earthquakes increases slowly up until 03-February-2002, goes up sharply between 03- February-2002 and 12- January-2010, and tends to fall thereafter.

When using Bayesian method, the prior distributions of the parameters <sup>1</sup> , 2 and , are taken to be respectively,

$$\begin{aligned} \mathbf{p}\_0(\lambda\_1, \lambda\_2 \mid \tau) &= \frac{1}{\lambda\_1 \lambda\_2}, \quad \lambda\_1, \lambda\_2 > 0\\ \mathbf{p}\_0(\tau) &= \frac{1}{\mathbf{T}}, \qquad 0 < \tau < \mathbf{T} \end{aligned} \tag{17}$$

Estimated change points are shown on the scatter plot of the time intervals between the

Fig. 5. The estimated occurrence rate of earthquakes in order of change points.

0

As can be seen in Figure 4 and Figure 5, the occurrence rate of the earthquakes increases slowly up until 03-February-2002, goes up sharply between 03- February-2002 and 12-

When using Bayesian method, the prior distributions of the parameters <sup>1</sup> , 2 and , are

 

<sup>1</sup> p( , ) , , 0

012 1 2 1 2

(17)

<sup>1</sup> p ( ) , 0 T <sup>T</sup>

Fig. 4. Multiple change points in the earthquake data.

January-2010, and tends to fall thereafter.

taken to be respectively,

earthquakes.

where the prior distributions are called uninformative priors. The joint posterior distribution of the parameters is can be written as:

$$\begin{aligned} \mathbf{p}\_1(\boldsymbol{\lambda}\_1, \boldsymbol{\lambda}\_2, \tau) &\approx \mathcal{L}(\boldsymbol{\lambda}\_1, \boldsymbol{\lambda}\_2, \tau) \mathbf{p}\_0(\boldsymbol{\lambda}\_1, \boldsymbol{\lambda}\_2 \mid \tau) \mathbf{p}\_0(\tau) \\ &\approx \boldsymbol{\lambda}\_1^{\text{N}\_\text{r}} \mathbf{e}^{\boldsymbol{\lambda}\_1 \cdot \tau} \boldsymbol{\lambda}\_2^{\text{n-N}\_\text{r}} \mathbf{e}^{\boldsymbol{\lambda}\_2(\mathcal{T} \cdot \tau)} \frac{1}{\boldsymbol{\lambda}\_1 \boldsymbol{\lambda}\_2 \mathbf{T}} \end{aligned} \tag{18}$$

With integrate Eq. (18) with respect to the parameters <sup>1</sup> , <sup>2</sup> , the marginal posterior distribution of change point is proportional to

$$\mathbf{p}\_1(\mathbf{r}) \propto \frac{\Gamma(\mathbf{N}\_\tau)\Gamma(\mathbf{n} - \mathbf{N}\_\tau)}{\pi^{\mathbf{N}\_\tau}(\mathbf{T} - \mathbf{r})^{\mathbf{n} - \mathbf{N}\_\tau}} \left(\frac{1}{\mathbf{T}}\right) \tag{19}$$

where (.) is gamma function. The marginal posterior distributions of <sup>1</sup> and 2 are not obtained analytically because of discontinuity in and the close form of the marginal posterior distributions of <sup>1</sup> and 2 are respectively proportional to

$$\begin{aligned} \mathbf{p}\_1(\lambda\_1) &\propto \frac{1}{\mathbf{T}} \Big\|\_{0}^{\mathrm{T}\_\mathrm{r}} \mathbf{e}^{-\lambda\_1 \mathsf{r}} \frac{\Gamma(\mathbf{n} - \mathbf{N}\_\mathrm{r})}{(\mathbf{T} - \mathbf{r})^{\mathrm{n} - \mathbf{N}\_\mathrm{r}}} \mathrm{d}\mathsf{r} \,\mathrm{d}\mathsf{r} \\\\ \mathbf{p}\_1(\lambda\_2) &\propto \frac{1}{\mathbf{T}} \Big\|\_{0}^{\mathrm{n} - \mathbf{N}\_\mathrm{r}} \mathbf{e}^{-\lambda\_2 (\mathbf{T} - \mathbf{r})} \frac{\Gamma(\mathbf{N}\_\mathrm{r})}{\mathbf{r}^{\mathrm{N}\_\mathrm{r}}} \mathrm{d}\mathsf{r} \end{aligned}$$

where, for the possible values of , s ,s , ,s 12 n , and the variable N takes the values 1,2, ,n .

Form Eq.(19), the Bayesian estimate ˆ of the change point in the earthquake data is found to be the date 03-February-2002 2002 which corresponds to the posterior mode. The posterior distribution of the change point is given in Figure 6.

Fig. 6. The posterior distribution of change point and the posterior mode.

Change Point Analysis in Earthquake Data 35

The sequence of the earthquake magnitudes Y ,Y , ,Y 12 n can be assumed to be normally distributed random variables. The change point model in the mean of the normally

2

Y ,Y , ,Y ~ N( , ) (20)

<sup>n</sup> 2 2

. (21)

2

Y ,Y , ,Y ~ N( , )

1 2 1

12 n 2

where is an unknown change point in the mean of the sequence of the normally distributed random variables. To investigate whether there is a change point in the earthquake magnitudes, under the change point model, the likelihood function of the

v n 2 2

<sup>2</sup> i1 i2 2 n i 1 i 1

 2 2 1 2 2 i1 i2

i 1 i 1

<sup>1</sup> nL( , , ) (y ) (y ) 2 (22)

1 1 exp (y ) (y )

 

1 2 i 1 i 2 i 1 iv1

(2 ) 2

and the logarithm of the likelihood function in Eq. (21) is:

L( , , ) f(y , , ) f(y , , )

 

Fig. 8. Histogram of the earthquake magnitudes.

distributed random sequences is written,

observed values, y12 n ,y , , y ,

The Bayesian estimates of the parameters <sup>1</sup> and 2 would be the same as the maximum likelihood estimates because of uninformative prior distributions.

For multiple change points, the binary Bayesian segmentation procedure can be easily used. The procedure is employed as the mode of the posterior distribution of change point decreases considerably until the one before the posterior mode is found. The results are given in Figure 7.


Fig. 7. The Bayesian estimations of parameters for multiple change points in the earthquake data (\*in days).

When comparing the maximum likelihood estimates and the Bayesian estimates of the change points, we can see that the change points corresponding to dates such as 9-March-1957, 10-May-1997, 03-February-2002 and 12-January-2010 overlap. These dates are investigated in depth from many aspects, including overall world temperature, other disasters, and astronomical events.

When we look at the magnitudes of the earthquakes in the data, those of magnitude 6 to 6.9 number 291 (35.57%) and those of magnitude 7 to 7.9 number 258 (31.54%) (Table 2). The histogram of the magnitudes is given in Figure 6. Furthermore, using the Kolmogorov-Smirnov test, the distribution of the magnitudes is found to be almost normal (p=0.000).


Table 2. The number of earthquakes with respect to magnitudes.

The Bayesian estimates of the parameters <sup>1</sup> and 2 would be the same as the maximum

For multiple change points, the binary Bayesian segmentation procedure can be easily used. The procedure is employed as the mode of the posterior distribution of change point decreases considerably until the one before the posterior mode is found. The results are

Change points and the occurrence rates of the earthquake throughout time axis 1ˆ 2ˆ 3ˆ 4ˆ 5ˆ 6ˆ

0.0093\* 0.0025\* 0.0136\* 0.0254\* 0.1246\* 0.2244\* 0.0284\*

Fig. 7. The Bayesian estimations of parameters for multiple change points in the earthquake

When comparing the maximum likelihood estimates and the Bayesian estimates of the change points, we can see that the change points corresponding to dates such as 9-March-1957, 10-May-1997, 03-February-2002 and 12-January-2010 overlap. These dates are investigated in depth from many aspects, including overall world temperature, other

When we look at the magnitudes of the earthquakes in the data, those of magnitude 6 to 6.9 number 291 (35.57%) and those of magnitude 7 to 7.9 number 258 (31.54%) (Table 2). The histogram of the magnitudes is given in Figure 6. Furthermore, using the Kolmogorov-Smirnov test, the distribution of the magnitudes is found to be almost normal (p=0.000).

Magnitude 4-4.9 5-5.9 6-6.9 7-7.9 8-8.9 9

earthquakes 54 93 291 258 40 82

Table 2. The number of earthquakes with respect to magnitudes.

% 0.0660 0.1137 0.3557 0.3154 0.0489 0.1002

03-Feb-2002 N<sup>ˆ</sup> =429

<sup>ˆ</sup> <sup>5</sup>

09-Aug-2009 N<sup>ˆ</sup> =771

<sup>ˆ</sup> <sup>6</sup>

12-Jan-2010 N<sup>ˆ</sup> =806

ˆ

<sup>ˆ</sup> <sup>7</sup>

10-May-1997 N<sup>ˆ</sup> =385

<sup>ˆ</sup> <sup>4</sup>

likelihood estimates because of uninformative prior distributions.

09-Mar-1957 N<sup>ˆ</sup> =185

<sup>ˆ</sup> <sup>3</sup>

given in Figure 7.

21-Dec-1954 N<sup>ˆ</sup> =183

<sup>ˆ</sup> <sup>2</sup>

disasters, and astronomical events.

1

0

data (\*in days).

The number of

Fig. 8. Histogram of the earthquake magnitudes.

The sequence of the earthquake magnitudes Y ,Y , ,Y 12 n can be assumed to be normally distributed random variables. The change point model in the mean of the normally distributed random sequences is written,

$$\begin{aligned} \mathbf{Y}\_{1}, \mathbf{Y}\_{2}, \dots, \mathbf{Y}\_{\mathbf{v}} &\sim \mathbf{N}(\mu\_{1}, \sigma^{2}) \\ \mathbf{Y}\_{\mathbf{v}+1}, \mathbf{Y}\_{\mathbf{v}+2}, \dots, \mathbf{Y}\_{\mathbf{n}} &\sim \mathbf{N}(\mu\_{2}, \sigma^{2}) \end{aligned} \tag{20}$$

where is an unknown change point in the mean of the sequence of the normally distributed random variables. To investigate whether there is a change point in the earthquake magnitudes, under the change point model, the likelihood function of the observed values, y12 n ,y , , y ,

$$\begin{split} \operatorname{L}(\theta\_{1}, \theta\_{2}, \mathbf{v}) &= \prod\_{i=1}^{\mathrm{v}} \operatorname{f}(\mathbf{y}\_{i}, \boldsymbol{\mu}\_{1}, \sigma^{2}) \prod\_{i=\mathrm{v}+1}^{\mathrm{n}} \operatorname{f}(\mathbf{y}\_{i}, \boldsymbol{\mu}\_{2}, \sigma^{2}) \\ &\propto \frac{1}{(\sqrt{2\pi\sigma^{2}})^{n}} \exp\left\{ -\frac{1}{2\sigma^{2}} \Big( \sum\_{i=1}^{\mathrm{v}} (\mathbf{y}\_{i} - \boldsymbol{\mu}\_{1})^{2} - \sum\_{i=\mathrm{v}+1}^{\mathrm{n}} (\mathbf{y}\_{i} - \boldsymbol{\mu}\_{2})^{2} \Big) \right\} . \end{split} \tag{21}$$

and the logarithm of the likelihood function in Eq. (21) is:

$$\ell \ln \text{L} \{ \theta\_1, \theta\_2, \mathbf{v} \} \propto -\frac{1}{2\sigma^2} \left( \sum\_{i=1}^{\nu} (\mathbf{y}\_i - \boldsymbol{\mu}\_1)^2 - \sum\_{i=1}^{\nu} (\mathbf{y}\_i - \boldsymbol{\mu}\_2)^2 \right) \tag{22}$$

Change Point Analysis in Earthquake Data 37

The expected number of earthquakes of magnitude y is computed by the multiplying the probability of occurrence class of magnitude with the expected number of earthquakes over

4-4.9 0.0028 0.0246 0.0218 0.2257 4.4306

5-5.9 0.0302 0.1402 0.1100 1.1404 0.8769

6-6.9 0.1609 0.4238 0.2629 2.7248 0.3670

7-7.9 0.4588 0.7565 0.2977 3.0856 0.3241

8-8.9 0.7834 0.9432 0.1598 1.6566 0.6037

9 0.9526 1.0000 0.0474 0.4912 2.0360

Table 4. The estimation of recurrence periods for certain earthquakes after the change point

Expected number of earthquakes

30 days 0.8520 6.0634

3 months 2.5560 18.1902

6 months 5.1120 36.3805

1 years 10.3660 73.7717

2 years 20.7320 147.5434

3 years 31.0980 221.3151

4 years 41.4640 295.0868

5 years 51.8300 368.8585

10 years 103.6600 737.7171

20 years 207.3200 1475.4342

with respect to the change point 12-January-2010.

Table 5. Expected number of earthquakes and the expected total magnitudes of earthquakes

P( L y U ) Expected

Expected number of earthquakes (year)

Expected total magnitudes of earthquakes

Average recurrence period (year)

F (U) <sup>Y</sup> Expected

a certain period of time.

F (L) <sup>Y</sup> Expected

Class Magnitude (*y*)

12-January-2010.

Time

For any given value of , the maximum estimates of unknown parameters 1, 2, are

 <sup>i</sup> i 1 1 y y and n i i 1 2 y y n respectively. The estimates are substituted into Eq.(22), and with the maximized the log likelihood function given by Eq.(22), we can obtain the maximum likelihood estimate of the change point,

$$\hat{\mathbf{v}} = \arg\max\_{\mathbf{k}=1,\ldots,\mathbf{n}\cdot\mathbf{1}} \left[ \frac{1}{2\sigma^2} \Big( \mathbf{v}\,\mathbf{\overline{y}}\_1^2 - (\mathbf{n}-\mathbf{v})\,\mathbf{\overline{y}}\_2^2 \Big) \right] \tag{23}$$

From Eq. (23), the change point is estimated to be the date 28-February-2011 corresponding to =817. The estimated change point is close to end of the sequence such that it refers to no change point in the sequence. We can model both the number of the earthquake occurrences and the magnitudes of the earthquakes with compound Poisson process. Let {X ,t 0} <sup>t</sup> be a compound Poisson process as defined in Subsection 2.1

$$\mathbf{X}\_{\mathbf{t}} = \sum\_{i=1}^{N\_i} \mathbf{Y}\_{\mathbf{i}}$$

where Nt is the number of earthquakes in time interval (0,t] , which is a Poisson distributed random variable with parameter t, and {Y ,i 1,2, } <sup>i</sup> are the magnitudes of the earthquakes, which are normally distributed random variables with parameters i and 2. The mean and variance of the magnitudes of the earthquakes according to the estimated change points dates, (9-March-1957, 10-May-1997, 03-February-2002 and 12-January-2010) are given in Table 3.


Table 3. Means and variances of the magnitudes with respect to estimated change point (\*in days).

Using Table 3, we compute the recurrence periods for certain earthquakes of magnitude y, corresponding to the change point, 12-January-2010. The recurrence period of an earthquake with magnitude y is computed by

Reccurence period=1/(Expected number of earthquakes).

For any given value of , the maximum estimates of unknown parameters 1, 2, are

with the maximized the log likelihood function given by Eq.(22), we can obtain the

From Eq. (23), the change point is estimated to be the date 28-February-2011 corresponding to =817. The estimated change point is close to end of the sequence such that it refers to no change point in the sequence. We can model both the number of the earthquake occurrences and the magnitudes of the earthquakes with compound Poisson process. Let {X ,t 0} <sup>t</sup> be a

 - 

<sup>2</sup> 1 2 k 1, n 1

 Ni t i i 1 X Y

where Nt is the number of earthquakes in time interval (0,t] , which is a Poisson distributed random variable with parameter t, and {Y ,i 1,2, } <sup>i</sup> are the magnitudes of the earthquakes, which are normally distributed random variables with parameters i and 2. The mean and variance of the magnitudes of the earthquakes according to the estimated change points dates, (9-March-1957, 10-May-1997, 03-February-2002 and 12-January-2010)

> Between 9-Mar-1957 and 10-May-1997

earthquakes 0.0090\* 0.0136\* 0.0254\* 0.1300\* 0.0284\*

magnitudes 6.9622 6.7250 7.0795 6.4562 7.1167

magnitudes 0.9305 0.8821 0.5579 0.9845 1.2706

Table 3. Means and variances of the magnitudes with respect to estimated change point

Using Table 3, we compute the recurrence periods for certain earthquakes of magnitude y, corresponding to the change point, 12-January-2010. The recurrence period of an earthquake

Reccurence period=1/(Expected number of earthquakes).

respectively. The estimates are substituted into Eq.(22), and

<sup>1</sup> <sup>ˆ</sup> arg max y (n )y 2 (23)

Between 10-May-1997 and 03-Feb-2002

Between 03-Feb-2002 and 12-Jan-2010

After 12- Jan-2010

2 2

y

are given in Table 3.

Estimated Change points

Occurrence rate of the

Mean of the

Variance of the

with magnitude y is computed by

(\*in days).

y and

2

y

 n

n

i i 1

y

maximum likelihood estimate of the change point,

compound Poisson process as defined in Subsection 2.1

Up to 9- Mar-1957

 <sup>i</sup> i 1 1


The expected number of earthquakes of magnitude y is computed by the multiplying the probability of occurrence class of magnitude with the expected number of earthquakes over a certain period of time.

Table 4. The estimation of recurrence periods for certain earthquakes after the change point 12-January-2010.


Table 5. Expected number of earthquakes and the expected total magnitudes of earthquakes with respect to the change point 12-January-2010.


Table 6. Probability of earthquake occurrences after the change point 12-January-2010.

Change Point Analysis in Earthquake Data 39

Predicting earthquakes is of crucial importance these days. Many researchers have studied earlier attempts at earthquake detection. The change point analysis is used in both backward (off-line) and forward (on-line) statistical research. In this study, it is used with the backward approach in the worldwide earthquake data. The change points found in the worldwide earthquake data are useful in making reliable inferences and interpreting the results for further research. Each date found as the change point in the earthquake data should be carefully investigated with respect to other geographical, ecological, and

Akman V. E. & Raftery A. E. (1986). Asymptotic inference for a change-point Poisson

Aktaş, S.; Konşuk, H. & Yiğiter, A. (2009). Estimation of change point and compound

Amorese, D. (2007). Applying a Change-Point Detection Method on Frequency-magnitude Distributions, *Bulletin of the Seismological of America*, Vol.97, No.5, pp.1742-1749 Boudjellaba, H.; MacGibbon, B. & Sawyer, P. (2001). On exact inference for change in a

Chen J. & Gupta A. K. (1997). Testing and locating variance change points with application

Fotopoulos S. B. & Jandhyala V. K. (2001). Maximum likelihood estimation of a change-point

Kagan, Y. Y. (2010). Statistical Distributions of Earthquake Numbers: Consequence of Branching Process, *Geophysical Journal International*, Vol.180, No.3, pp.1313-1328 Leonard T. & Papasouliotis, O. (2001). A Poisson Model for Identifying Characteristic Size

Lindsey H. (June 2007). 'Natural' catastrophes increasing worldwide, 14.06.2007, Available from: http://www.wnd.com/news/article.asp?ARTICLE\_ID=40584/ Mandeville M. W. (June 2007). "Eight Charts Which Prove That Chandler's Wobble Causes

Parzen, E. (1962). *Stochastic Processes*, Holden-Day, ISBN 0-8162-6664-6, San Francisco, USA Raftery A. E. & Akman V. E. (1986). Bayesian analysis of a Poisson process with a change

Rotondi R. & Garavaglia, E. (2002). Statistical Analysis of the Completeness of a Seismic

Poisson process parameters for the earthquake data in Turkey, *Environmetrics*, Vol.

Poisson Sequence. *Communications in Statistics A: Theory and Methods,* Vol.30, No.3,

to stock prices, Journal of the American Statistical Association: JASA, Vol.92, No.

for exponentially distributed random variables, *Statistics & Probability Letters*, 51:

Effects in Frequency Data: Application to Frequency-size Distributions for Global Earthquakes, "Starquakes", and Fault Lengths, *Journal of Geophysical Research*,

Earthquakes, Volcanism, El Nino and Global Warming", 14.06.2007, Available from: http://www.michaelmandeville.com/polarmotion/spinaxis/vortex\_correlations2.

process, *The Annals of Statistics*, Vol.14, No.4, pp. 15831590.

**5. Conclusion** 

**6. References** 

geological events or structures in the world.

20, No.4, pp. 416-427

pp. 407434

423429

htm

438, pp. 739-747

Vol.106, No.B7, pp.13,413-13,484

point, *Biometrika*, Vol.73, pp.8589

Catalogue, *Natural Hazards*, Vol.25, No. 3, pp. 245-258
