Package 'ETAS'

Create an Earthquake Catalog

Description

Creates an object of class "catalog" representing an earthquake catalog dataset. An earthquake catalog is a chronologically ordered list of time, epicenter and magnitude of all recorded earthquakes in geographical region during a specific time period.

Usage

catalog(data, time.begin=NULL, study.start=NULL,
          study.end=NULL, study.length=NULL,
          lat.range=NULL, long.range=NULL,
          region.poly=NULL, mag.threshold=NULL,
          flatmap=TRUE, dist.unit = "degree", 
          roundoff=TRUE, tz="GMT")

Arguments

data	A data.frame containing date, time, latitude, longitude and magnitude of earthquakes.
time.begin	The beginning of time span of the catalog. A character string or an object that can be converted to date-time (calendar dates plus time to the nearest second) by as.POSIXlt. The default NULL sets it to the date-time of the first event.
study.start	The start of the study period. A character string or an object that can be converted to date-time by as.POSIXlt. If not specified (NULL), then time.begin is used.
study.end	The end of the study period. A character string or an object that can be converted to date-time by as.POSIXlt. The default NULL sets it to the date-time of the last event.
study.length	A single numeric value specifying the length of the study period in decimal days. Incompatible with study.end: either study.end or study.length can be specified, but not both.
lat.range	The latitude range of a rectangular study region. A numeric vector of size 2 giving (latmin, latmax). By default (NULL) the range of the latitudes of events is used.
long.range	The longitude range of a rectangular study region. A numeric vector of size 2 giving (longmin, longmax). By default (NULL) the range of the longitudes of events is used.
region.poly	Polygonal boundary of a non-rectangular study region. A list with components lat and long of equal length specifying the coordinates of the vertices of a polygonal study region. The vertices must be listed in anticlockwise order.
mag.threshold	The magnitude threshold of the catalog. A positive numeric value. The default (NULL) sets it to the minimum magnitude of all events.
flatmap	Logical flag indicating whether to transform the spherical coordinates $(long, lat)$ on the earth surface to flat map (planar) coordinates $(x, y)$ in order to approximate the great-circle distance on the sphere by the corresponding Euclidean distance on the flat map.
dist.unit	A character string specifying the unit of geographical coordinates and spatial distances between events. Options are "degree" (the default case) and "km".
roundoff	Logical flag indicating whether to account for round-off error in coordinates of epicenters.
tz	A character string specifying the time zone to be used for the date-time conversion in as.POSIXlt. The default "GMT" is the UTC (Universal Time, Coordinated).

Details

The data is required to have at least 5 columns with names date, time, lat, long and mag containing, respectively, the date, time, latitude, longitude and magnitude of each event in the catalog.

The geographical study region can be rectangular or polygonal:

• rectangular study region can be specified by lat.range and long.range which must be numeric vectors of length 2.

• polygonal study region can be specified by region.poly which contains coordinates of the vertices of the polygon. It must be either a list with components lat and long of equal length or a data.frame with columns lat and long. The vertices must be listed in anticlockwise order and no vertex should be repeated (i.e. do not repeat the first vertex).

The function inside.owin in the spatstat is used to indicate whether events lie inside the study region. Only events inside the study region and the study period (study.start, study.end) are considered as target events. Other events are assumed to be complementary events.

If the events in data are not chronologically sorted, then a warning will be produced and the events will be sorted in ascending order with respect to time of occurrence.

If flatmap=TRUE, longitude-latitude coordinates convert to flat map coordinates:

• if dist.unit="degree", then the Equirectangular projection

  $x = \cos(cnt.lat/180 \pi) (long - cnt.long)$

  and $y = lat - cnt.lat$ is used to obtain the flat map coordinates $(x, y)$ in degrees, where $cnt.lat$ and $cnt.long$ are, respectively, the latitude and longitude of the centroid of the geographical region.

• if dist.unit="km", then the projection

  $x = 111.32 \cos(lat/180 \pi) long$

  and $y = 110.547 lat$ is used where $x$ and $y$ are in (approximate) kilometers.

Value

An object of class "catalog" containing an earthquake catalog dataset.

Author(s)

Abdollah Jalilian [email protected]

References

Zhuang J (2012). Long-term Earthquake Forecasts Based on the Epidemic-type Aftershock Sequence (ETAS) Model for Short-term Clustering. Research in Geophysics, 2(1), 52–57. doi:10.4081/rg.2012.e8.

See Also

etas.

Examples

summary(iran.quakes)

  # creating a catalog with rectangular study region
  iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1985/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)

  # equivalently, specifying the length of the study period
  iran.cat2 <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1985/01/01", study.length=11322,
     lat.range=c(25, 42), long.range=c(42, 63),
     mag.threshold=4.5)

  print(iran.cat2)


  # specifying a polygonal geographical region
  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))
  # creating a catalog with polygonal study region
  japan.cat <- catalog(japan.quakes, time.begin="1966-01-01",
      study.start="1970-01-01", study.end="2010-01-01",
      region.poly=jpoly, mag.threshold=4.5)

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

---

Convert date-time data to numeric data in decimal days

Description

A function to convert date-time data to decimal days with respect to a date-time origin.

Usage

date2day(dates, start=NULL, tz="", ...)

Arguments

dates	A date-time or date object. Typically, it is a character vector containing date-time information.
start	A date-time or date object. Determines the origin of the conversion.
tz	Optional. Timezone specification to be used for the conversion.
...	Arguments to be passed to as.POSIXlt.

Details

The arguments dates and start must be of appropriate format to be passed to as.POSIXlt function.

Value

A numeric vector of the same length as dates.

Author(s)

Abdollah Jalilian [email protected]

See Also

as.POSIXlt and difftime for appropriate format of the data to be converted.

Examples

# date-time data of Iran's earthquakes between 1973/01/01 and 2016/01/01
  dt <- paste(iran.quakes$date, iran.quakes$time)
  # origin of the conversion
  start <- "1973/01/01 00:00:00"
  # time in days since 1973-01-01 (UTC)
  date2day(dt, start, tz="GMT")

---

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 = NULL, bwd = NULL, nnp = 5, bwm = 0.05,
      verbose = TRUE, plot.it = FALSE, ndiv = 1000, no.itr = 11,
      rel.tol=1e-03, eps = 1e-06, cxxcode = TRUE, nthreads = 1,
      mver = 1)

Arguments

object	An object of class "catalog" containing an earthquake catalog dataset.
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 win. A numeric vector which has the length of the number of events. If not supplied, the following arguments nnp and bwm determine bandwidths.
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.
ndiv	An integer indicating the number of knots on each side of the geographical region for integral approximation.
no.itr	An integer indicating the number of iterations for convergence of the iterative approach of simultaneous estimation and declustering algorithm. See details.
rel.tol	Relative iteration convergence tolerance of the iterative estimation approach.
eps	Optimization convergence tolerance in the Davidon-Fletch-Powell algorithm
cxxcode	Logical flag indicating whether to use the C++ code. The C++ code is slightly faster and allows parallel computing.
nthreads	An integer indicating number of threads in the parallel region of the C++ code
mver	An integer indicating which spatial probability density function for locations of triggered events should be use. The default mver=1 corresponds to the inverse power density and mver=2 corresponds to the Gaussian density.

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 Hawkes 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) = A\exp(\alpha(m - m_0))$

$g(t)$:

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

$g(t) = \frac{p-1}{c}(1 + \frac{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 (mver = 1)

$f(x, y|m) = \frac{q-1}{\pi \sigma(m))} (1 + \frac{x^2 + y^2}{\sigma(m)})^{-q}$

or the light tail Gaussian density (mver = 2, only can be used if cxxcode = TRUE)

$f(x,y|m)= \frac{1}{2\pi \sigma(m)}\exp(-\frac{x^2 + y^2}{2\sigma(m)})$

with

$\sigma(m) = D\exp(\gamma(m - m_0))$

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 $u_0(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_i}) - \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

$\phi_i = \frac{\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 $u_0(x, y)$ is updated. These steps are repeated until the estimates converge (stabilize).

The no.itr argument specifies the maximum number of iterations in the iterative simultaneous estimation and declustering algorithm. The estimates often converge in less than ten iterations. The relative iteration convergence tolerance and the optimization convergence tolerance are, respectively, determined by rel.tol and eps arguments. The progress of the computations can be traced by setting the verbose and plot.it arguments to be TRUE.

If cxxcode = TRUE, then the internal function etasfit uses the C++ code implemented using the Rcpp package, which allows multi-thread parallel computing on multi-core processors with OpenMP. The argument nthreads in this case determines the number of threads in the parallel region of the code. If nthreads = 1 (the default case), then a serial version of the C++ code carries out the computations.

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 non-stationarity (e.g. increasing or cyclic 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 [email protected]

References

Ogata Y (1988). Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes. Journal of the American Statistical Association, 83(401), 9–27. doi:10.2307/2288914.

Ogata Y (1998). Space-time Point-process Models for Earthquake Occurrences. Annals of the Institute of Statistical Mathematics, 50(2), 379–402. doi:10.1023/a:1003403601725.

Zhuang J, Ogata Y, Vere-Jones D (2002). Stochastic Declustering of Space-Time Earthquake Occurrences. Journal of the American Statistical Association, 97(458), 369–380. doi:10.1198/016214502760046925.

Zhuang J, Ogata Y, Vere-Jones D (2006). Diagnostic Analysis of Space-Time Branching Processes for Earthquakes. In Case Studies in Spatial Point Process Modeling, pp. 275–292. Springer Nature. doi:10.1007/0-387-31144-0_15.

Zhuang J (2011). Next-day Earthquake Forecasts for the Japan Region Generated by the ETAS Model. Earth, Planets and Space, 63(3), 207–216. doi:10.5047/eps.2010.12.010.

See Also

catalog for constructing data. probs for estimated declustering probabilities. resid.etas for diagnostic plots.

Examples

# fitting the ETAS model to an Iranian catalog
  # preparing the catalog
  iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1986/01/01", study.end="2016/01/01",
     lat.range=c(26, 40), long.range=c(44, 63), mag.threshold=5)
  print(iran.cat)
  ## Not run: 
  plot(iran.cat)
## End(Not run)

  # setting initial parameter values
  param0 <- 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=param0)
## End(Not run)


  # fitting the ETAS model to an Italian catalog
  # preparing the catalog
  italy.cat <- catalog(italy.quakes, dist.unit="km")
  ## Not run: 
  plot(italy.cat)
## End(Not run)

  # setting initial parameter values
  mu <- 1
  k0 <- 0.005
  c <- 0.005
  alpha <- 1.05
  p <- 1.01
  D <- 1.1
  q <- 1.52
  gamma <- 0.6
  # reparametrization: transform k0 to A
  A <- pi * k0 / ((p - 1) * c^(p - 1) * (q - 1) * D^(q - 1))
  param0 <- c(mu, A, c, alpha, p, D, q, gamma)

  # fitting the model
  ## Not run: 
  nthreads <- parallel::detectCores()
  italy.fit <- etas(italy.cat, param0, nthreads=nthreads)
## End(Not run)


  # fitting the ETAS model to a Japanese catalog
  # setting the target polygonal study region
  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))
  # preparing the catalog
  japan.cat <- catalog(japan.quakes, study.start="1953-05-26",
      study.end="1990-01-08", region.poly=jpoly, mag.threshold=4.5)
  ## Not run: 
  plot(japan.cat)
## End(Not run)

  # setting initial parameter values
  param0 <- c(0.592844590, 0.204288231, 0.022692883, 1.495169224,
  1.109752319, 0.001175925, 1.860044210, 1.041549634)

  # fitting the model
  ## Not run: 
  nthreads <- parallel::detectCores()
  japan.fit <- etas(japan.cat, param0, nthreads=nthreads)
## End(Not run)

---

Class of Fitted ETAS Models

Description

A class etas to represent a fitted ETAS model. The output of etas.

Details

An object of class etas represents an ETAS model that has been fitted to a spatio-temporal point pattern (catalog) of earthquakes. It is the output of the model fitter, etas.

The class etas has methods for the following standard generic functions:

generic	method	description
print	print.etas	print details

Author(s)

Abdollah Jalilian [email protected]

See Also

etas,

Examples

# fitting the ETAS model to an Iranian catalog

  data(iran.quakes)
  summary(iran.quakes)

    # fitting the ETAS model to an Iranian catalog
  # preparing the catalog
  iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1986/01/01", study.end="2016/01/01",
     lat.range=c(26, 40), long.range=c(44, 63), mag.threshold=5)
  print(iran.cat)
  ## Not run: 
  plot(iran.cat)
## End(Not run)

  # setting initial parameter values
  param0 <- 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=param0)
## End(Not run)

---

An Iranian Earthquake Catalog

Description

A data frame with 5970 rows and 5 columns giving occurrence date, time, longitude, latitude and magnitude of shallow earthquakes (depth < 100 km) occurred since 1973-01-01 till 2016-01-01 in Iran and its vicinity (40-65E and 22-42N). Only earthquakes with magnitude greater than or equal to 4.5 are included.

Usage

data(iran.quakes)

Format

An object of class "data.frame" containing the following columns:

• date Occurrence date in the format "yyyy-mm-dd"

• time Occurrence time (UTC) in the format "hh:mm:ss"

• long Latitude of epicenter in decimal degrees

• lat Latitude of epicenter in decimal degrees

• mag Magnitude in body-wave magnitude scale (mb)

Source

The ANSS Comprehensive Catalog (ComCat): https://earthquake.usgs.gov/earthquakes/search/

Examples

summary(iran.quakes)

  gregion <- list(lat = c(26, 25, 29, 38, 35), long = c(52, 59, 58, 45, 43))
  # creat an earthquake catalog
  iran.cat <- catalog(iran.quakes, study.start = "1991/01/01",
     study.end = "2011/01/01", region.poly = gregion, mag.threshold = 4.5)


  ## Not run: 
  plot(iran.cat)
  iran.fit <- etas(iran.cat)
## End(Not run)

  zagros <- list(lat = c(27, 26, 29, 29, 35, 33),
     long = c(52, 58, 58, 54, 48, 46))
  iran.cat <- catalog(iran.quakes, study.start = "1991/01/01",
     study.end = "2011/01/01", region.poly = zagros, mag.threshold = 4)

---

An Italian Earthquake Catalog

Description

A data frame with 2158 rows and 6 columns giving occurrence date, time, longitude, latitude, magnitude and depth of earthquakes occurred since 2005-04-16 till 2013-11-01 in Italy and its vicinity (6.15-19E and 35-48N) with magnitude greater than or equal to 3.

Usage

data(italy.quakes)

Format

An object of class "data.frame" containing the following columns:

• date Occurrence date in the format "yyyy/mm/dd"

• time Occurrence time in decimal days after the first earthquake

• long Latitude of epicenter in decimal degrees

• lat Latitude of epicenter in decimal degrees

• mag Magnitude of each earthquake by JMA (Japan Meteorological Agency)

• depth Depth of each earthquake

Source

Data are retrieved from the Italian Seismological Instrumental and Parametric Data Base (ISIDE) available at http://iside.rm.ingv.it by Marcello Chiodi and Giada Adelfio in the etasFLP package.

References

Marcello Chiodi and Giada Adelfio (2015). etasFLP: Mixed FLP and ML Estimation of ETAS Space-Time Point Processes. R package version 1.3.0. https://CRAN.R-project.org/package=etasFLP.

Examples

# creat an earthquake catalog
  italy.cat <- catalog(italy.quakes, dist.unit="km")

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

  # set initial parameter values
  mu <- 1
  k0 <- 0.005
  c <- 0.005
  alpha <- 1.05
  p <- 1.01
  D <- 1.1
  q <- 1.52
  gamma <- 0.6
  # reparametrization: transform k0 to A
  A <- pi * k0 / ((p - 1) * c^(p - 1) * (q - 1) * D^(q - 1))
  param0 <- c(mu, A, c, alpha, p, D, q, gamma)

  ## Not run: 
  italy.fit <- etas(italy.cat, param0)
## End(Not run)

---

A Japanese Earthquake Catalog

Description

A data frame with 13724 rows and 6 columns giving occurrence date, time, longitude, latitude, magnitude and depth of shallow earthquakes (depth < 100 km) occurred since 1926-01-08 till 2007-12-29 in Japan and its vicinity (128-145E and 27-45N). Only earthquakes with magnitude greater than or equal to 4.5 are included.

Usage

data(japan.quakes)

Format

An object of class "data.frame" containing the following columns:

• date Occurrence date in the format "yyyy/mm/dd"

• time Occurrence time in the format "hh:mm:ss.ss"

• long Latitude of epicenter in decimal degrees

• lat Latitude of epicenter in decimal degrees

• mag Magnitude of each earthquake by JMA (Japan Meteorological Agency)

• depth Depth of each earthquake

Source

Data are retrieved from the Japan Meteorological Agency (JMA) by Jiancang Zhuang and acoppany the Fortran code at http://bemlar.ism.ac.jp/zhuang/software.html.

References

Zhuang J, Ogata Y, Vere-Jones D (2002). Stochastic Declustering of Space-Time Earthquake Occurrences. Journal of the American Statistical Association, 97(458), 369–380. doi:10.1198/016214502760046925.

Zhuang J, Ogata Y, Vere-Jones D (2006). Diagnostic Analysis of Space-Time Branching Processes for Earthquakes. In Case Studies in Spatial Point Process Modeling, pp. 275–292. Springer Nature. doi:10.1007/0-387-31144-0_15.

Examples

# set the target polygonal study region
  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))

  # creat an earthquake catalog
  japan.cat <- catalog(japan.quakes, study.start="1953-05-26",
      study.end="1990-01-08", region.poly=jpoly, mag.threshold=4.5)

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

  param0 <- c(0.592844590, 0.204288231, 0.022692883, 1.495169224,
  1.109752319, 0.001175925, 1.860044210, 1.041549634)

  ## Not run: 
  nthreads <- parallel::detectCores()
  japan.fit <- etas(japan.cat, param0, nthreads=nthreads)
## End(Not run)

---

Clustering Part of Conditional Intensity Function of the ETAS Model

Description

A function to compute the clustering part of the conditional intensity function of the ETAS model at specified time and location.

Usage

lambda(t, x, y, param, object)

Arguments

t	A numeric value. The time that the conditional intensity is to be computed at.
x	A numeric value. The x-coordinate of the location that the conditional intensity is to be computed at.
y	A numeric value. The y-coordinate of the location that the conditional intensity is to be computed at.
param	Vector of model parameters.
object	An object of class "catalog" containing an earthquake catalog dataset.

Details

For a given $t$, $x$ and $y$, this function computes

$\sum_{t_i < t} k(m_i)g(t - t_i)f(x - x_i, y - y_i|m_i).$

Value

A numeric value.

Author(s)

Abdollah Jalilian [email protected]

References

Zhuang J, Ogata Y, Vere-Jones D (2002). Stochastic Declustering of Space-Time Earthquake Occurrences. Journal of the American Statistical Association, 97(458), 369–380. doi:10.1198/016214502760046925.

Zhuang J, Ogata Y, Vere-Jones D (2006). Diagnostic Analysis of Space-Time Branching Processes for Earthquakes. In Case Studies in Spatial Point Process Modeling, pp. 275–292. Springer Nature. doi:10.1007/0-387-31144-0_15.

See Also

etas catalog

Examples

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)

  param <- c(0.46, 0.23, 0.022, 2.8, 1.12, 0.012, 2.4, 0.35)

  ## Not run: 
  lambda(15706, 40.12, 34.5, param, iran.cat)
## End(Not run)

---

Declustering Probabilities, Background Seismicity Rate and Clustering Coefficient

Description

Functions to estimate the declustering probabilities, background seismicity rate and clustering (triggering) coefficient for a fitted ETAS model.

Usage

probs(fit)
  rates(fit, lat.range = NULL, long.range = NULL,
        dimyx=NULL, plot.it=TRUE)

Arguments

fit	A fitted ETAS model. An object of class "etas".
lat.range	Latitude range of the rectangular grid. A numeric vector of length 2.
long.range	Longitude range of the rectangular grid. A numeric vector of length 2.
dimyx	Dimensions of the rectangular discretization grid for the geographical study region. A numeric vector of length 2.
plot.it	Logical flag indicating whether to plot the rates or return them as pixel images.

Details

The function probs returns estimates of the declustering probabilities

$p_j = 1 - \frac{\mu(x_j, y_j)}{lambda(t_j, x_j, y_j|H_{t_j})}$

where $1-p_j$ is the probability that event $j$ is a background event.

The function rates returns kernel estimate of the background seismicity rate $\mu(x,y)$ and the clustering (triggering) coefficient

$\omega(x,y)=1-\frac{\mu(x,y)}{\Lambda(x,y)}$

where $\Lambda(x,y)$ is the total spatial intensity function.

The argument dimyx determines the rectangular discretization grid dimensions. If it is given, then it must be a numeric vector of length 2 where the first component dimyx[1] is the number of subdivisions in the y-direction (latitude) and the second component dimyx[2] is the number of subdivisions in the x-direction (longitude).

Value

If plot.it=TRUE, the function produces plots of the background seismicity and total spatial rate, clustering coefficient and conditional intensity function at the end of study period.

If plot.it=FALSE, it returns a list with components

• bkgd the estimated background siesmicity rate

• total the estimated total spatial rate

• clust the estimated clustering coefficient

• lamb the estimated conditional intensity function at time $t=t_{\mathrm{start}}$

Author(s)

Abdollah Jalilian [email protected]

References

Zhuang J, Ogata Y, Vere-Jones D (2002). Stochastic Declustering of Space-Time Earthquake Occurrences. Journal of the American Statistical Association, 97(458), 369–380. doi:10.1198/016214502760046925.

Zhuang J, Ogata Y, Vere-Jones D (2006). Diagnostic Analysis of Space-Time Branching Processes for Earthquakes. In Case Studies in Spatial Point Process Modeling, pp. 275–292. Springer Nature. doi:10.1007/0-387-31144-0_15.

Zhuang J (2011). Next-day Earthquake Forecasts for the Japan Region Generated by the ETAS Model. Earth, Planets and Space, 63(3), 207–216. doi:10.5047/eps.2010.12.010.

See Also

etas

Examples

# 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 and
  ## Not run: 
  iran.fit <- etas(iran.cat, param0=param01)
## End(Not run)

  # estimating the declustering probabilities
  ## Not run: 
  pr <- probs(iran.fit)
  plot(iran.cat$longlat.coord[,1:2], cex=2 * (1 - pr$prob))
## End(Not run)

  # estimating the  background seismicity rate and clustering coefficient
  ## Not run: 
  rates(iran.fit, dimyx=c(100, 125))
  iran.rates <- rates(iran.fit, dimyx=c(200, 250), plot.it=FALSE)
  summary(iran.rates$background)
## End(Not run)

---

Residuals Analysis and Diagnostics Plots

Description

A function to compute and plot spatial and temporal residuals as well as transformed times for a fitted ETAS model.

Usage

resid.etas(fit, type="raw", n.temp=1000, dimyx=NULL)

Arguments

fit	A fitted ETAS model. An object of class "etas".
type	A character string specifying the type residuals to be computed. Options are "raw" (the default case), "reciprocal" and "pearson".
n.temp	An integer specifying the number of partition points for temporal residuals.
dimyx	Dimensions of the discretization for the smoothed spatial residuals. A numeric vector of length 2.

Details

The function computes the temporal residuals

$R^{temp}(I_j, h) = \sum_{i=1}^{N} \delta_i 1[t_i \in I_j] h(t_i) \lambda^{temp}(t_i|H_{t_i}) - \int_{I_j} h(t)\lambda^{temp}(t|H_t) d t$

for $I_j=((j-1)T/n.temp, jT/n.temp]$, $j=1,...,n.temp$, and the (smoothed version of) spatial residuals

$R^{spat}(B_j, h) = h(\tilde{x}_i, \tilde{y}_i) \lambda^{spat}(\tilde{x}_i, \tilde{y}_i)(\tilde{\delta}_i - \tilde{w}_i)$

for a Berman-Turner quadrature scheme with quadrature points $(\tilde{x}_i, \tilde{y}_i)$ and quadrature weights $\tilde{w}_i$, $i=1,...,n.spat$. Raw, reciprocal and Pearson residuals obtain with $h=1$, $h=1/\lambda$ and $h=1/\sqrt{\lambda}$, respectively.

In addition, the function computes transformed times

$\tau_j=\int_{0}^{t_j} \lambda^{temp}(t|H_t) d t$

and

$U_j = 1 - \exp(-(t_j - t_{j-1}))$

Value

The function produces plots of temporal and smoothed spatial residuals, transformed times $\tau_j$ against $j$ and Q-Q plot of $U_j$.

It also returns a list with components

• tau the transformed times

• U related quantities with the transformed times

• tres the temporal residuals

• sres the smoothed spatial residuals

Author(s)

Abdollah Jalilian [email protected]

References

Baddeley A, Rubak E, Turner R (2015). Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press, London. https://www.routledge.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/p/book/9781482210200.

Baddeley A, Turner R (2000). Practical Maximum Pseudolikelihood for Spatial Point Patterns. Australian & New Zealand Journal of Statistics, 42(3), 283–322. doi:10.1111/1467-842X.00128.

Baddeley A, Turner R, Moller J, Hazelton M (2005). Residual Analysis for Spatial Point Processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(5), 617–666. doi:10.1111/j.1467-9868.2005.00519.x.

Ogata Y (1988). Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes. Journal of the American Statistical Association, 83(401), 9–27. doi:10.2307/2288914.

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.

See Also

etas

Examples

iran.cat <- catalog(iran.quakes, time.begin="1973/01/01",
     study.start="1986/01/01", study.end="2016/01/01",
     lat.range=c(26, 40), long.range=c(44, 63), mag.threshold=5)
  print(iran.cat)
  ## Not run: 
  plot(iran.cat)
## End(Not run)

  # setting initial parameter values
  param0 <- 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=param0)

  # diagnostic plots
  resid.etas(iran.fit)
## End(Not run)

---

Retrieving data from the ISC Bulletin

Description

Searches the International Seismological Centre (ISC) bulletin for recorded eqrthqukes ocurred during a specified time period in a specified geographical region (rectangular or circular) with a specified magnitude type and magnitude threshould.

Usage

search.isc(start.year=1900, start.month=1, start.day=01,
             end.year=2018, end.month=12, end.day=31,
             searchshape="RECT", lat.bot=NULL, lat.top=NULL,
             long.left=NULL, long.right=NULL, lat.ctr=NULL, 
             long.ctr=NULL, radius=NULL, dist.units="deg", 
             dep.min=0, dep.max=100, nulldep=TRUE,
             mag.min=4, mag.max=NULL, mag.type='MB', 
             mag.agcy=NULL, mirror=TRUE)

Arguments

start.year	A single numeric value specifying the start year of the time period.
start.month	A single numeric value (1 to 12) specifying the start month of the time period.
start.day	A single numeric value (1 to 31) specifying the start day of the time period.
end.year	A single numeric value specifying the end year of the time period.
end.month	A single numeric value (1 to 12) specifying the end month of the time period.
end.day	A single numeric value (1 to 31) specifying the end day of the time period.
searchshape	A character string specifying the shape of the search region. Rectangular search with searchshape="RECT" (default) and circular search searchshape="CIRC" are possible.
lat.bot	A single numeric value (-90 to 90) specifying the bottom latitude of the rectangular search region (if searchshape="RECT").
lat.top	A single numeric value (-90 to 90) specifying the top latitude of the rectangular search region (if searchshape="RECT").
long.left	A single numeric value (-180 to 180) specifying the left longitude of the rectangular search region (if searchshape="RECT").
long.right	A single numeric value (-180 to 180) specifying the right longitude of the rectangular search region (if searchshape="RECT").
lat.ctr	A single numeric value (-90 to 90) specifying the central latitude of the circular search region (if searchshape="CIRC").
long.ctr	A single numeric value (-180 to 180) specifying the central longitude of the circular search region (if searchshape="CIRC").
radius	A single numeric value specifying the search radius of the circular search region (if searchshape="CIRC").
dist.units	A character string specifying the unit of the search radius: degrees with "deg" and kilometers with "km" are the available options.
dep.min	A single numeric value specifying the minumum depth (km).
dep.max	A single numeric value specifying the maximum depth (km).
nulldep	Logical flag indicating whether to include events with unknown depths.
mag.min	A single numeric value specifying the minumum magnitude.
mag.max	A single numeric value specifying the maximum magnitude.
mag.type	A character string specifying the magnitude type. Due to the large number of magnitude types and the various notation used, the selected type will search all magnitudes of that type. E.g. the default "MB" will match all body-wave magnitudes: mb, mB, MB, mb1mx, etc.
mag.agcy	A character string specifying the magnitude author. "ISC", "NEIC", "GCMT/HRVD", "JMA", amongst others are considered reliable magnitude authors. The default uses "Any", which will match any author that has a magnitude for both prime and secondary hypocentres included in an event. Selecting "prime author" will limit the search by magnitudes directly associated with the prime hypocentre (usually the same author as the prime hypocentre) - i.e. for events with an ISC authored prime hypocentre, selecting 'Prime author', will limit events to those with ISC determined magnitudes.
mirror	Logical flag indicating whether to use the ISC mirror server hosted at IRIS DMC (http://isc-mirror.iris.washington.edu) as it may be faster.

Details

Any use of data from the ISC by academic or commercial organisations, as well as individuals should always be cited. The correct format for citations may be found on http://www.isc.ac.uk/citations/.

The ISC is named as a valid data centre for citations within American Geophysical Union (AGU) publications. As such, please follow the AGU guidelines when referencing ISC data in one of their journals. The ISC may be cited as both the institutional author of the Bulletin and the source from which the data were retrieved.

Value

An object of class "data.frame" which can be passed to the "catalog" function in order to create an earthquake catalog dataset.

Author(s)

Abdollah Jalilian [email protected]

References

International Seismological Centre (2015). On-line Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent., Thatcham, United Kingdom.

International Seismological Centre (2015). Reference Event Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent., Thatcham, United Kingdom.

International Seismological Centre (2015). EHB Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent., Thatcham, United Kingdom.

See Also

catalog.

Examples

# rectangular search
  ## Not run: 
  mydata <- search.isc(start.year=1980, start.month=1, start.day=01,
                   end.year=2017, end.month=11, end.day=11, 
                   lat.bot=33, lat.top=37,
                   long.left=44, long.right=48,
                   mag.min=3.5, mag.type='MB')  
  mycatalog <- catalog(mydata, study.start = "1990-01-01")
  plot(mycatalog)
  
## End(Not run)
  # circular search
  ## Not run: 
  mydata2 <- search.isc(start.year=2017, start.month=11, start.day=12,
                   end.year=2018, end.month=3, end.day=01,
                   searchshape="CIRC", 
                   lat.ctr=34.905, long.ctr=45.956,
                   radius=150, dist.units="km",
                   mag.min=3.5, mag.type='ML') 
  mycatalog2 <- catalog(mydata2)
  plot(mycatalog2)
  
## End(Not run)

---