Programs ETAS and RETAS

Estimation of parameter values in the ETAS model and diagnostic analysis of seismic activities

1. The ETAS model and the residual point process

1.1 The MLE computation

The FORTRAN program [etas.f] computes the maximum likelihood estimates (MLEs) of five parameters of the Epidemic Type Aftershock Sequence (ETAS) model for seismic activities (Ogata, 1988) from a set of data on the occurrence times and magnitudes of earthquakes. An example of data [work.etas] is given in the same directory. The executive file of the present FORTRAN source code is compiled by GNU-FORTRAN on the free UNIX Operating System gCygwinh on the Microsoft Windows XP.

The ETAS model is a point process model representing the activity of earthquakes of magnitude Mz and larger occurring in a certain region during a certain interval of time. Total number of such earthquakes is denoted by N. The seismic activity includes primary activity of constant occurrence rate  in time (Poisson process). Each earthquake (including aftershock of another earthquake) is followed by its aftershock activity, though only aftershocks of magnitude Mz and larger are included in the data. The aftershock activity is represented by the Omori-Utsu formula (modified Omori formula; see section 1 of Program AFT) in the time domain. The rate of aftershock occurrence at time t following the ith earthquake (time: ti, magnitude: M_i) is given by

n_i(t) = K exp[(M_i - Mz)] / (t - t_i + c)^p,  for  tt_i     (1)

where K, , c, and p are constants, which are common to all aftershock sequences in the region. The rate of occurrence of the whole earthquake series at time t becomes

(t)  =  + _i  n_i(t). @  (2)

The summation is done for all i satisfying t_it. Five parameters , K, c, , and p represent characteristics of seismic activity of the region. The data set should includes the occurrence times ti of earthquakes associate with Mi in the period [Tstart, Tend] such as the prototype data file [work.etas] given in this directory as an example.

1.2 The residual point process of the ETAS model

The FORTRAN program [retas.f] computes the following output for displaying the goodness-of-fit of the ETAS model to the data. The cumulative number of earthquakes at time t since t0 is given by the integration of (2) with respect to the time t,

(t) = (t-t0) + K_i exp[(M_i-Mz)]{c^(1-p) - (t-t_i+c)^(1-p)} / (p-1),   (3)

where the summation of i is taken for all data event. The output of the program is given in the file [work.res] which includes the column of {(t_i), I = 1, 2, c, N}. This file is used for displaying the cumulative curve and magnitude v.s. transformed time (t_i). For the users who have the free graphic software R, see Section 3 (Graphs) below. I wrote a R-module [r.retas] to display them, in addition to the module [r.seis] for displaying the cumulative curve and magnitude v.s. ordinary time ti. If the observed rate of occurrence is compared with the calculated one from the model, period of decreased or increased seismic activity (relative quiescence or activation) can be recognized. [r.seisetas] See section 3 (Graphs) below.
  
2. Programs and data

2.1 Exact and approximate algorithms

A FORTRAN program [etas.f] consists of two (exact and approximated) versions of the calculation algorithm for the maximization of likelihood. The approximated one (Ogata et al., 1993) to reduce the processing time for relatively large data sets (number of quakes N larger than several hundreds). The processing time is proportional to NxN for the exact version but proportional to N for the approximation version. The approximated version uses the subroutine [func9] which runs at one of the five levels 1, 2, 4, 8, and 16. The higher level means faster processing but lower accuracy. At levels 1 or 2, the results are very close to the exact solution by the exact version which uses the subroutine [func4]. The subroutines and the levels are selected by the control file [etas.open], in addition to other input variables such as the range of data components that restrict data and initial values of parameters.

The maximum number of earthquakes can be processed by ETAS and RETAS is 17,777. The data file may contain as many as 17,777 quakes based on the set dimensions of the ETAS program.

Starting from a given set of initial guess of the parameters that is given in the control file [etas.open], the program [etas.f] repeats calculations of the function values and its gradients at each step of parameter vector. At each cycle of iteration, the linearly searched step (lambda), negative log-likelihood value (-LL), and two estimates of square sum of gradients are shown. At each linear search of the likelihood computation, -LL and the 5 parameter values are shown. -LL decreases and tends to a final value. With the convergence of -LL, one of the square sum of gradients becomes nearly zero and the iteration is terminated. The time of the termination is displayed. In the present case, the following results are displayed, where AIC = -2 LL + 2x5. If -LL value and gradients do not seem to converge, we recommend you to use the last parameter estimates as the initial values in the control file [etas.open] and restart the execution.

2.2 Precursory and target intervals

The period for which the ETAS model parameters are computed is called target interval. The seismicity in this period may be affected by earthquakes which occurred before this period due to the long-lived nature of aftershock activity. To consider this effect, a time interval precursory to the target interval (called precursory interval) is chosen and aftershock activities following earthquakes in this period are taken into computation. When the ETAS parameters are estimated for two or more successive periods (usually divided at turning points of seismicity), it is better to set a precursory period before each target period (see Ogata, 1992, for instance).

         Prec. interval         Target interval        Prediction interval
   |--------------------|--------------------------|--------------------------|
@zts                 tstart                      zte                        ztend


2.3 Data file [work.etas]

The data file for ETAS should have a a sequential file with a 9 columns-format of the sequential number, longitude, latitude, magnitude, time from the mainshock in days, depth, year, month, and day, as is given by the prototype data file [work.etas]. The time is usually measured in days. Time data t(i) in the 5th column should be given at least to the fourth place of decimals (in double precision), since the differences between them are important. The first row of the data is the title. Any copy of [work.etas] is saved by the name [*.etas] where * is an alphabetical word.

2.4 Control file [etas.open]

Furthermore, the file [etas.open] controls the data reading such as [work.etas] and other input variables such as the range of data components that restrict data and initial values of parameters, ztstart, ztend, which represent the beginning and the end of the interval for which data are prepared, respectively. Any copy of [etas.open] is saved by the name [etas.*] where * is an alphabetical word, corresponding to the data [*.etas].

The first row consists of two numbers; either 4 or 9 for the selection of the subroutine [func4] or [func9], respectively, and then one of 1, 2, 4, 8, or 16 when [func9] is used; when [func4] is used any number should be set.

The second row consists of three numbers; zts, zte and tstart in the above illustrated variables when the file is applied to the program [etas.f]. When the file is applied to the program [retas.f] for the diagnostic outputs, the last number can be ztend (>= zte) that is the last time of the period of available data. The example shows the first 6.5days of aftershock period since the mainshock due to [work.etas], where the target interval is [0.03, 6.5]days and the data of first 0.03days period is only used as the history of the ETAS model. 

The third row consists of two numbers; threshold magnitude Mth and reference magnitude Mz. Here I have taken the magnitude of the mainshock as the reference magnitude, but you may take Mz = Mth for general seismicity. Note that the MLE of the parameter K differs depending on the reference parameter. 

The fourth row provides the initial estimates of the five parameters , K, c,  and p when this is applied the program [etas.f], and this is the last row when we use for the program [etas.f]. If you set = 0.0 in the initial value, this variable remains 0.0 throughout the calculations. Also, if you set p = 1.0 (the original Omori formula), this variable remains 1.0 throughout the calculations. When this control file is used for the program [retas.f] to make the diagnostic outputs, these values in this row have to be the MLEs estimated by the program [etas.f]. 

The fifth row consists of three numbers; zts, zte and tarst. This row is only used when the file is applied to the program [retas.f] for the diagnostic outputs, the last number must be zte that is the end of the target period. 

The sixth row is the same as the fourth row, only used for the program [retas.f], and the seventh row of the value of the negative of the maximum likelihood value corresponding to the MLEs, used just for the userf note to compare with other cases or models.

3. Graphs

The R language modules [r.seis], [r.seisetas] and [r.retas] are for the graphs of cumulative number and magnitude of earthquakes against the ordinary time and transformed time, respectively. The transformation is made so that the theoretical cumulative frequency curve becomes a straight line. These use the files [work.etas] and [work.res], respectively. The last graph includes printed number with certain time values for vertical dotted lines showing target interval and MLEs in control file [etas.open]. The outputs of the graphs in postscript format, named as seis.ps, seisetas.ps, and retas.ps are given in the same directory.

The R language modules [r.retas] and [r.seisetas] shows the theoretical curve of cumulative frequency (red color) which is a straight line by definition of the transformed time, and the cumulative frequency curves of the data (red color). 

The graphical free software R is also installed in the other directory based on the Microsoft Windows XP (RGui). To command in the R console, write 
> source(er.*f) 
where * correspond above R-module, then the graphical window should appear and draw appropriate figures.

4. Remarks

4.1 Application to aftershock sequences

When the ETAS model is applied to an aftershock sequence, the main shock should always be included in the data. Since the primary seismicity (backgroud seismicity)  is usually absent in the aftershock sequence data, it may be preferable to use a model with=0. To constrain  to be 0, set 0.0 as the initial value for m in the control file [etas.open].

4.2 Abnormal termination

If the data are very different from the proposed model, an overflow error may occur or the solution may not converge within the fixed number of iterations (30N, N is the number of parameters). The results in this case are often incorrect and not recommended to adopt as the estimates.

4.3 Test data set

The test data made from the JMA earthquake catalogs named as [work.etas] is attached, which is same as [miyagiN.etas]. This is the aftershock data of 26th July 2003 earthquake of M6.2 at the northern Miyagi-Ken, northern Japan. There is another dataset [fukuokaW.etas], which is the aftershock data of 16th August 2005 earthquake of M7.2 at the offshore of western Fukuoka-Ken, Kyushu, Japan.  

References

Ogata, Y., Statistical models for earthquake occurrences and residual analysis for point processes, J. Amer. Statis. Assoc., 83, No.401, 9-27, 1988.

Ogata, Y., Statistical model for standard seismicity and detection of anomalies by residual analysis for point process, Tectonophysics, 169, 1-16, 1989.

Ogata, Y., Detection of precursory relative quiescence before great earthquakes through a statistical model, J. Geophys. Res., 97, 19845-19871, 1992.

Ogata, Y., R.S. Matsu'ura, and K. Katsura, Fast likelihood computation of Epidemic Type Aftershock Sequence model, Geophys. Res. Lett., 20, 2143-2146, 1993.

Utsu, T., Y. Ogata, and R. S. Matsu'ura, The centenary of the Omori formula for a decay law of aftershock activity, J. Phys. Earth, 43, 1-33, 1995.


