**Code for the
space-time ETAS model and stochastic declustering **

Source code and
example dataset:

System requirement:

OS:
Linux

Software: gfortran (or Intel fortran), mpich2

The code here is in a
parallel version. Suppose that you have already installed mpich2. (Not yet?!
Please go to http://www.mcs.anl.gov/research/projects/mpich2/ to have a look.)

(If
you need a non-parallel version, please use this
file. This program is coded by Dr. Qi Wang by modifying the parallel
version. However, I cannot provide much technical support for this serial
version)

After downloading sd.tar.gz,
please uncompress the file in a fold by

*$ tar -zxvf
sd.tar.gz *

and compile it by

*$ mpif77 *.f -o
etas8p*

and run the example
dataset by

*$ export
GFORTRAN_UNBUFFERED_ALL='y'*

*$ mpiexec -n 8
etas8p < jap.in *

There is 1 example
dataset:

jap.in: answers to the
input prompts.

alljpM45.etas: JMA
catalogue of shallow eqs. M>=4.5.

jap.out: outputs to
stdout.

Note: the polygon must
be input in the counterclockwise order.

The outputs of the
program will be

para: MLE of paramters
for each step (mu, A, c, alpha, p,
D, q,

gamma).

rates.dat: first row: m n

first column of Row 2 to Row last: m x n matrix of total rate

second column of Row 2 to Row Last: m x n matrix of

background rate
divided by mu

third column of Row 2 to Row Last: m x n matrix of

intensity rate at the end of time span of dataset.

probs.dat: first column : EQ sequence n

second column: background probability for this event.

third column: bandwidth

fourth column: background rate

pmatr.dat
: each row: j, i, rho_{ij} (probabability that j is trigger by i)

!Note:

2.
In this etas8p program, every earthquake is
considered as a point source, with its hypocenter or epicenter being the
initiating point of the rapture.

3.
Before you run the program, a careful
preparation of the input data is important. The input file that contains the
earthquake catalog should satisfy the following conditions.

(a)
The occurrence times of events in the input
catalog cannot be reversed or overlapped.

It frequently happens
that two or more events are not in the right chronical order in the catalog. In
such a case, etas8p simply stops and gives a warning. The user is required to
correct the input data file by putting these events in the correct chronic
order. Etas8p may fail to check overlapping times. When two or more events have
overlapping occurrence times, it is difficult to get the likelihood
optimization in etas8p to converge, mostly yielding c=0. The user can avoid
such a difficulty by adding some very small times to occurrence times of the
second, third, c events in the group of time-overlapping events and
putting them in a chronically ascending order.

(b) The occurrence
locations of events in the input catalog cannot be overlapped. Overlapping in
locations always happens within an aftershock sequence, when the records of
earthquake locations are rounded at a fixed digit and when there are many
events occurring in a small area. If such situations exist, optimizer usually
gives D=0. To avoid that, the user is suggested that he bring back the rounding
errors to the locations. If a location is rounded at the second digit, for
example, (138.45, 34.60), two random numbers that are uniformly distributed in
(-0.005, 0.005) can be added to the longitude 138.45 and the latitude 34.60.

4.
When setting up the input control parameters,
please note the followings:

(a) The catalog must be
complete above the magnitude threshold in the target polygon region, and there
should be no increasing trend in the occurrence rate of events in the catalog.
To find this you can simply draw some plots of (1) the magnitude of each event
vs the number of events before it, (2) the magnitude of each event vs number of
events with locations of less longitudes, and (3) the magnitude of each event
vs number of events with locations of less latitudes. The user can also add
some rounding errors, to each magnitude as well. For a complete catalog, the
magnitude distribution looks homogeneous along the horizontal axis.
@

(b) There are 2 times to
set up in the input. The begin time is for background seismicity calculation,
and the start time is the time window for fitting the model.

(c)
If some typical initial values of the ETAS
parameters are given, for example, a set of general values of (mu, A, c, alfa,
p, d, q, gamma) can be 0.5, 0.2, 0.02, 1.5, 1.1, 0.001, 1.5, 1.0, the result
should converge to the same.

(d) The target polygon must be input in the
counter-clockwise order.

5. On the output:

(a)
The gradient should be close to zero,
otherwise the fitting is divergent.

(b) The p value should not
be less than 1.0, with 1.1 is a most typical value.

If you find this
program useful in your research, please refer to the following articles in your
publications.

1.
__Zhuang J.__ (2006) *Second-order residual analysis of spatiotemporal
point processes and applications in model evaluation.* **Journal of the
Royal Statistical Society: Series B (Statistical Methodology)**, 68 (4), 635-653.
doi: 10.1111/j.1467-9868.2006.00559.x.

2.
__Zhuang J.__, Ogata Y. and Vere-Jones D. (2004). *Analyzing earthquake
clustering features by using stochastic reconstruction*. **Journal of
Geophysical Research**, 109, No. B5, B05301, doi:10.1029/2003JB002879.

3.
__Zhuang J.__, Ogata Y. and Vere-Jones D. (2002). *Stochastic declustering
of space-time earthquake occurrences*. **Journal of the American
Statistical Association**, 97: 369-380.

*updated on 09/June/2017*.