# 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

1 2 |

### Arguments

`object` |
An object of class |

`param0` |
Initial guess for model parameters. A numeric vector of appropriate length (currently 8). See details. |

`bwd` |
Optional. Bandwidths for smoothness and integration
on the geographical region |

`nnp` |
Number of nearest neighbors for bandwidth calculations. An integer. |

`bwm` |
Minimum bandwidth. A positive numeric value. |

`verbose` |
Logical flag indicating whether to print progress reports. |

`plot.it` |
Logical flag indicating whether plot probabilities of each event being a background event on a map. |

`no.itr` |
An integer indicating the number of iterations for convergence of the iterative approach of declustering algorithm. See details. |

### 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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 | ```
# 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)
``` |