Fit the space-time ETAS model to data

Description

A function to fit the space-time version of the Epidemic Type Aftershock Sequence (ETAS) model to a catalog of earthquakes (a spatio-temporal point pattern) and perform a stochastic declustering method.

Usage

 etas(object, param0, bwd = NULL, nnp = 5, bwm = 0.05,
      verbose = TRUE, plot.it = FALSE, no.itr = 11)

Arguments

Details

Ogata (1988) introduced the epidemic type aftershock sequence (ETAS) model based on Gutenberg-Richter law and modified Omori law. In its space-time representation (Ogata, 1998), the ETAS model is a temporal marked point process model, and a special case of marked Hawks process, with conditional intensity function

lambda(t, x, y | H_t ) = mu(x, y) + sum[t[i] < t] k(m[i]) g(t - t[i]) f(x - x[i], y - y[i]|m[i])

where

H_t:

is the observational history up to time t, but not including t; that is

H_t = {(t[i], x[i], y[i], m[i]): t[i] < t}

mu(x, y):

is the background intensity. Currently it is assumed to take the semi-parametric form

mu(x, y) = mu u(x, y)

where mu is an unknown constant and u(x, y) is an unknown function.

k(m):

is the expected number of events triggered from an event of magnitude m given by

k(m[i]) = A exp(alpha(m - m0))

g(t):

is the p.d.f of the occurrence times of the triggered events, taking the form

g(t) = ((p - 1)/c)(1 + t/c)^(-p)

f(x, y|m):

is the p.d.f of the locations of the triggered events, considered to be either the long tail inverse power density

f(x, y|m) = (q - 1)/(pi sigma(m)) (1 + (x^2 + y^2)/(sigma(m)))^(-q)

or the light tail Gaussian density (currently not implemented)

f(x,y|m)= exp(-(x^2 + y^2)/(2 sigma(m)))/(2 pi sigma(m))

with

sigma(m) = D exp( gamma (m - m0) )

The ETAS models classify seismicity into two components, background seismicity mu(x,y) and clustering seismicity lambda(t, x, y|H_t) - mu(x, y), where each earthquake event, whether it is a background event or generated by another event, produces its own offspring according to the branching rules controlled by k(m), g(m) and f(x, y|m).

Background seismicity rate u(x, y) and the model parameters

theta = (mu, A, c, alpha, p, D, q, gamma)

are estimated simultaneously using an iterative approach proposed in Zhuang et al. (2002). First, for an initial u0(x, y), the parameter vector theta is estimated by maximizing the log-likelihood function

l(theta) = sum_[i] lambda(t[i], x[i], y[i]|H_t) - int lambda(t, x, y|H_t) dx dy dt.

Then the procedure calculates the probability of being a background event for each event in the catalog by

mu(x[i], y[i])/lambda(t[i], x[i], y[i]|H_t[i]).

Using these probabilities and kernel smoothing method with Gaussian kernel and appropriate choice of bandwidth (determined by bwd or nnp and bwm arguments), the background rate u0(x, y) is updated. These steps are repeated enough times (determined by no.itr argument) such that the results converge.

This version of the ETAS model assumes that the earthquake catalog is complete and the data are stationary in time. If the catalog is incomplete or there is instationarity (e.g. increasing trend) in the time of events, then the results of this function are not reliable.

Value

A list with components

param:

The ML estimates of model parameters.

bk:

An estimate of the u(x, y).

pb:

The probabilities of being background event.

opt:

The results of optimization: the value of the log-likelihood function at the optimum point, its gradient at the optimum point and AIC of the model.

rates:

Pixel images of the estimated total intensity, background intensity, clustering intensity and conditional intensity.

Note

This function is based on a C port of the original Fortran code by Jiancang Zhuang, Yosihiko Ogata and their colleagues. The etas function is intended to be used for small and medium-size earthquake catalogs. For large earthquake catalogs, due to time-consuming computations, it is highly recommended to use the parallel Fortran code on a server machine. The Fortran code (implemented for parallel/non-parallel computing) can be obtained from http://bemlar.ism.ac.jp/zhuang/software.html.

Author(s)

Abdollah Jalilian jalilian@razi.ac.ir

References

Zhuang, J., Ogata, Y. and Vere-Jones, D. (2005). Diagnostic analysis of space-time branching processes for earthquakes. Lecture Note in Statistics: Case Studies in Spatial Point Process Models (Baddeley, A., Gregori, P., Mateu, J., Stoica, R. and Stoyan, D.), Springer-Verlag, New York, 185, 276–292.

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.

Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. Annals of the Institute of Statistical Mathematics, 50, 379–402.

Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes. Journal of American Statistical Association, 83, 9–27.

See Also

catalog for constructing data.

Examples

  # fitting the ETAS model to an Iranian catalog

  summary(iran.quakes)

  # preparing the catalog
  iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1996/01/01", study.end="2016/01/01",
     lat.range=c(25, 42), long.range=c(42, 63), mag.threshold=4.5)

  print(iran.cat)
  ## Not run: 
  plot(iran.cat)
## End(Not run)

  # initial parameters values
  param01 <- c(0.46, 0.23, 0.022, 2.8, 1.12, 0.012, 2.4, 0.35)

  # fitting the model
  ## Not run: 
  iran.fit <- etas(iran.cat, param0=param01, no.itr=5)
## End(Not run)



  # fitting the ETAS model to a Japanese catalog

  # preparing the catalog
  jpoly <- list(long=c(134.0, 137.9, 143.1, 144.9, 147.8,
      137.8, 137.4, 135.1, 130.6), lat=c(31.9, 33.0, 33.2,
      35.2, 41.3, 44.2, 40.2, 38.0, 35.4))
  jap.cat <- catalog(jap.quakes, time.begin="1966-01-01",
      study.start="1970-01-01", study.end="2010-01-01",
      region.poly=jpoly, mag.threshold=4.5)

  ## Not run: 
  plot(jap.cat)
## End(Not run)

  # initial parameters values
  param00 <- c(0.6, 0.2, 0.02, 1.5, 1.1, 0.0012, 1.8, 1.0)

  # fitting the model
  ## Not run: 
  jap.fit <- etas(jap.cat, param0=param00)
## End(Not run)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(ETAS)
Loading required package: maps

 # maps v3.1: updated 'world': all lakes moved to separate new #
 # 'lakes' database. Type '?world' or 'news(package="maps")'.  #


> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ETAS/etas.Rd_%03d_medium.png", width=480, height=480)
> ### Name: etas
> ### Title: Fit the space-time ETAS model to data
> ### Aliases: etas
> ### Keywords: spatial math earthquake modeling
> 
> ### ** Examples
> 
>   # fitting the ETAS model to an Iranian catalog
> 
>   summary(iran.quakes)
         date               time           long            lat       
 2014-08-18:  44   00:04:50.70:   2   Min.   :40.01   Min.   :22.13  
 2011-10-23:  42   04:46:33.00:   2   1st Qu.:48.34   1st Qu.:28.32  
 2013-04-09:  31   07:02:27.00:   2   Median :51.88   Median :32.13  
 2013-05-11:  27   00:00:19.42:   1   Mean   :52.20   Mean   :32.80  
 2011-10-24:  19   00:00:41.44:   1   3rd Qu.:56.39   3rd Qu.:37.52  
 1990-06-21:  18   00:00:46.52:   1   Max.   :65.00   Max.   :41.98  
 (Other)   :5789   (Other)    :5961                                  
      mag       
 Min.   :4.000  
 1st Qu.:4.200  
 Median :4.400  
 Mean   :4.468  
 3rd Qu.:4.700  
 Max.   :6.200  
                
> 
>   # preparing the catalog
>   iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
+      study.start="1996/01/01", study.end="2016/01/01",
+      lat.range=c(25, 42), long.range=c(42, 63), mag.threshold=4.5)
> 
>   print(iran.cat)
earthquake catalog:
  time begin 1973-01-01 
  study period: 1996-01-01  to  2016-01-01 (T = 7305 days)
geographical region:
    rectangular = [ 42 , 63 ] x [ 25 25 ]
threshold magnitude: 4.5
number of events:
  total events 2959 : 1240 target events,  1719 complementary events
  ( 67 events outside geographical region, 1652 events outside study period)>   ## Not run: 
> ##D   plot(iran.cat)
> ## End(Not run)
> 
>   # initial parameters values
>   param01 <- c(0.46, 0.23, 0.022, 2.8, 1.12, 0.012, 2.4, 0.35)
> 
>   # fitting the model
>   ## Not run: 
> ##D   iran.fit <- etas(iran.cat, param0=param01, no.itr=5)
> ## End(Not run)
> 
> 
> 
>   # fitting the ETAS model to a Japanese catalog
> 
>   # preparing the catalog
>   jpoly <- list(long=c(134.0, 137.9, 143.1, 144.9, 147.8,
+       137.8, 137.4, 135.1, 130.6), lat=c(31.9, 33.0, 33.2,
+       35.2, 41.3, 44.2, 40.2, 38.0, 35.4))
>   jap.cat <- catalog(jap.quakes, time.begin="1966-01-01",
+       study.start="1970-01-01", study.end="2010-01-01",
+       region.poly=jpoly, mag.threshold=4.5)
> 
>   ## Not run: 
> ##D   plot(jap.cat)
> ## End(Not run)
> 
>   # initial parameters values
>   param00 <- c(0.6, 0.2, 0.02, 1.5, 1.1, 0.0012, 1.8, 1.0)
> 
>   # fitting the model
>   ## Not run: 
> ##D   jap.fit <- etas(jap.cat, param0=param00)
> ## End(Not run)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>