Code for the space-time ETAS model and stochastic declustering

Source code and example dataset:

sd.tar.gz

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

background rate divided by mu

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

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:

1. The ETAS model describes the triggering relationship between two earthquakes. It does not require that the triggering event has a larger or smaller magnitude than the triggered events. Thus, this model does not specify aftershocks, mainshocks or foreshocks.

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, … 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.

(c) By some experiments, the rounding of magnitudes may influence the convergence of the optimizer when the number of earthquakes in the catalog is small. This needs further confirmation.

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.