Movatterモバイル変換

Type:	Package
Title:	Modeling Earthquake Data Using 'ETAS' Model
Description:	Fits the space-time Epidemic Type Aftershock Sequence ('ETAS') model to earthquake catalogs using a stochastic 'declustering' approach. The 'ETAS' model is a 'spatio-temporal' marked point process model and a special case of the 'Hawkes' process. The package is based on a Fortran program by 'Jiancang Zhuang' (available athttps://bemlar.ism.ac.jp/zhuang/software.html), which is modified and translated into C++ and C such that it can be called from R. Parallel computing with 'OpenMP' is possible on supported platforms.
Version:	0.7.2
Date:	2025-10-23
Author:	Abdollah Jalilian [aut, cre], Jiancang Zhuang [ctb]
Maintainer:	Abdollah Jalilian <stat4aj@gmail.com>
Depends:	R (≥ 3.5.0), stats, graphics, utils, maps
Imports:	lattice, goftest, spatstat.geom, spatstat.explore,spatstat.random, Rcpp, fields
LinkingTo:	Rcpp (≥ 1.0.0)
License:	GPL-2 |GPL-3 [expanded from: GPL (≥ 2)]
LazyData:	true
LazyLoad:	true
NeedsCompilation:	yes
URL:	https://github.com/jalilian/ETAS
BugReports:	https://github.com/jalilian/ETAS/issues
RoxygenNote:	7.3.3
Packaged:	2025-10-23 21:08:30 UTC; aj
Repository:	CRAN
Date/Publication:	2025-10-24 05:20:30 UTC

Internal ETAS functions

Description

Internal ETAS functions.

Usage

rates.inter(theta, object, bwd, lat.range=NULL, long.range=NULL,            mver=1, dimyx=NULL, plot.it=TRUE)lambdax(rt, rx, ry, theta, revents)etasfit(theta, revents, rpoly, tperiod, integ0, ihess, verbose,    ndiv, eps, cxxcode, nthreads, mver)decluster(theta, rbwd, revents, rpoly, tperiod, ndiv, cxxcode, mver)cxxfit(tht, revents, rpoly, tperiod, rinteg0, ihess, ndiv, eps, verbose, nthreads, mver)cxxdeclust(param, revents, rpoly, bwd, tperiod, ndiv, mver)cxxrates(param, revents, bwd, tperiod, gx, gy, mver)cxxtimetrans(theta, revents, rpoly, tperiod, integ0, ndiv, mver)cxxlambdtemp(tg, theta, revents, rpoly, tperiod, integ0, ndiv, mver)cxxlambspat(xg, yg, theta, revents, rpoly, tperiod, bwd, mver)cxxSmooth(x, y, bwd, gx, gy, expand)cxxstpoisstest(xrank, yrank, M)cxxstpoisstestMP(xrank, yrank, M, nthreads)cdeclust(theta, rbwd,revents, rpoly, tperiod)cfit(theta, rdata, ihess, rverbose)clambdax(rt, rx, ry, theta, revents)timetransform(fit)lambdaspatial(x, y, fit)lambdatemporal(t, fit)longlat2xy(long, lat, region.poly, dist.unit="degree")xy2longlat(x, y, region.poly, dist.unit="degree")poiss.test(object, which="joint", r=NULL, lambda=NULL, bwd=NULL, dimyx=NULL,   nsim=299, n.perm=1000, verbose=TRUE, nthreads=1)Smoothcatalog(object, type="spatial", bwd=NULL, bwm=NULL, nnp=NULL,   dimyx=NULL, convert=FALSE)allanfactor(object, K=100, nsim=500, cat.name=NULL)morisitaindex(object, K=11, bwd=NULL, dimyx=NULL, cat.name=NULL)decimalplaces(x) roundoffErr(x)ttrend(obj, tmax=NULL, tmin=NULL)simetas(param, bkgd, sim.start, sim.end=NULL, sim.length=NULL,         lat.range=NULL, long.range=NULL, region.poly=NULL,         mag.threshold=NULL, flatmap=TRUE, dist.unit="degree",         cat0=NULL, roundoff=FALSE, tz="GMT")## S3 method for class 'etas'print(x,...)## S3 method for class 'etas'plot(x, which="est", dimyx=NULL, ...)## S3 method for class 'catalog'print(x,...)## S3 method for class 'catalog'plot(x,...)

Details

These are usually not to be called by the user.

Create an Earthquake Catalog

Description

Creates an object of class"catalog" representingan earthquake catalog dataset. An earthquake catalog is achronologically ordered list of time, epicenter and magnitudeof all recorded earthquakes in geographical region duringa 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

Details

Thedata is required to have at least 5 columns with namesdate,time,lat,long andmagcontaining, respectively, the date, time, latitude, longitudeand magnitude of each event in the catalog.

The geographical study region can be rectangular or polygonal:

• rectangular study region can be specified bylat.rangeandlong.range which must be numeric vectors of length 2.

• polygonal study region can be specified byregion.polywhich contains coordinates of the vertices of the polygon. It mustbe either alist with componentslat andlongof equal length or adata.frame with columnslatandlong. The vertices must be listed inanticlockwise order and no vertex should be repeated(i.e. do not repeat the first vertex).

The functioninside.owin in thespatstatis 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 astarget events. Other events are assumed to becomplementary events.

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

Ifflatmap=TRUE, longitude-latitude coordinates convert toflat map coordinates:

• ifdist.unit="degree", then theEquirectangular projection

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

  andy = lat - cnt.latis used to obtain the flat map coordinates(x, y) indegrees, wherecnt.lat andcnt.long are,respectively, the latitude and longitude of the centroid of thegeographical region.

• ifdist.unit="km", then the projection

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

  andy = 110.547 latis used wherex andy are in (approximate) kilometers.

Value

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

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

References

Zhuang J (2012).Long-term Earthquake Forecasts Based on the Epidemic-type AftershockSequence (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

Details

The argumentsdates andstart must be ofappropriate format to be passed toas.POSIXlt function.

Value

A numeric vector of the same length asdates.

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

See Also

as.POSIXlt anddifftime for appropriate format of the datato 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 TypeAftershock Sequence (ETAS) model to a catalog of earthquakes(a spatio-temporal point pattern) and perform a stochasticdeclustering 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

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 atemporal marked point process model, and a special case of markedHawkes 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 notincluding t; that is

H_t=\{(t_i, x_i, y_i, m_i): t_i < t\}

\mu(x,y):

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

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

where\mu is an unknown constant andu(x,y)is an unknown function.

k(m):

is the expected number of events triggeredfrom an event of magnitudem given by

k(m) = A\exp(\alpha(m - m_0))

g(t):

is the p.d.f of the occurrencetimes 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 ofthe locations of the triggered events, considered to beeither 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 ifcxxcode = 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, backgroundseismicity\mu(x, y) and clustering seismicity\lambda(t, x, y|H_t) - \mu(x, y), whereeach earthquake event, whether it is a background event or generated byanother event, produces its own offspring according to the branching rulescontrolled byk(m),g(m) andf(x, y|m).

Background seismicity rateu(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 initialu_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 backgroundevent 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 kerneland appropriate choice of bandwidth (determined bybwd ornnpandbwm arguments), the background rateu_0(x, y)is updated. These steps are repeated until the estimates converge(stabilize).

Theno.itr argument specifies the maximum number of iterationsin the iterative simultaneous estimation and declustering algorithm.The estimates often converge in less than ten iterations. The relativeiteration convergence tolerance and the optimizationconvergence tolerance are, respectively, determined byrel.tol andeps arguments.The progress of the computations can be traced by settingtheverbose andplot.it arguments to beTRUE.

Ifcxxcode = TRUE, then the internal functionetasfituses the C++ code implemented using theRcpp package,which allows multi-thread parallel computing on multi-core processorswith OpenMP.The argumentnthreads in this case determinesthe number of threads in the parallel region of the code.Ifnthreads = 1 (the default case), then a serial version of theC++ code carries out the computations.

This version of the ETAS model assumes that the earthquake catalogis complete and the data are stationary in time. If the catalogis incomplete or there is non-stationarity (e.g. increasing or cyclictrend) in the time of events, then the results of this function arenot reliable.

Value

A list with components

param:

The ML estimates of model parameters.

bk:

An estimate of theu(x, y).

pb:

The probabilities of being background event.

opt:

The results of optimization: the value of the log-likelihoodfunction 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 aC port of the originalFortran code by Jiancang Zhuang, Yosihiko Ogata andtheir colleagues. Theetas function is intended to beused for small and medium-size earthquake catalogs.For large earthquake catalogs, due to time-consumingcomputations, it is highly recommended touse the parallelFortran code on a server machine.TheFortran code (implemented forparallel/non-parallel computing) can be obtained fromhttps://bemlar.ism.ac.jp/zhuang/software.html.

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

References

Ogata Y (1988).Statistical Models for Earthquake Occurrences and Residual Analysis forPoint 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.InCase 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 bythe 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 classetas to represent a fitted ETAS model.The output ofetas.

Details

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

The classetas has methods for the followingstandard generic functions:

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

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 thanor equal to 4.5 are included.

Usage

data(iran.quakes)

Format

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

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

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

• longLatitude of epicenter in decimal degrees

• latLatitude of epicenter in decimal degrees

• magMagnitude 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 earthquakesoccurred 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:

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

• timeOccurrence time in decimal days after the first earthquake

• longLatitude of epicenter in decimal degrees

• latLatitude of epicenter in decimal degrees

• magMagnitude of each earthquake by JMA (Japan Meteorological Agency)

• depthDepth of each earthquake

Source

Data are retrieved from the Italian Seismological Instrumental andParametric Data Base (ISIDE) by Marcello Chiodi and Giada Adelfio in theetasFLP package.

References

Marcello Chiodi and Giada Adelfio (2015). etasFLP: Mixed FLP and MLEstimation of ETAS Space-Time Point Processes. R package version1.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 till2007-12-29 in Japan and its vicinity (128-145E and 27-45N).Only earthquakes with magnitude greater than or equal to 4.5are included.

Usage

data(japan.quakes)

Format

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

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

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

• longLatitude of epicenter in decimal degrees

• latLatitude of epicenter in decimal degrees

• magMagnitude of each earthquake by JMA (Japan Meteorological Agency)

• depthDepth of each earthquake

Source

Data are retrieved from the Japan Meteorological Agency (JMA) byJiancang Zhuang and acoppany theFortran code athttps://bemlar.ism.ac.jp/zhuang/software.html.

References

Zhuang J, Ogata Y, Vere-Jones D (2002). Stochastic Declustering ofSpace-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-TimeBranching Processes for Earthquakes. InCase Studies in Spatial PointProcess 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 conditionalintensity function of the ETAS model at specified time and location.

Usage

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

Arguments

Details

For a givent,x andy, this functioncomputes

\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 Jalilianjalilian@razi.ac.ir

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.InCase Studies in Spatial Point Process Modeling,pp. 275–292. Springer Nature.doi:10.1007/0-387-31144-0_15.

See Also

etascatalog

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 andClustering Coefficient

Description

Functions to estimate the declustering probabilities, background seismicityrate 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

Details

The functionprobs 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})}

where1-p_j is the probability that eventjis a background event.

The functionrates returns kernel estimate of the backgroundseismicity 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 intensityfunction.

The argumentdimyx determines the rectangular discretizationgrid dimensions. If it is given, then it must be a numeric vectorof length 2 where the first componentdimyx[1] is thenumber of subdivisions in the y-direction (latitude) and thesecond componentdimyx[2] is the number of subdivisionsin the x-direction (longitude).

Value

Ifplot.it=TRUE, the function produces plots of thebackground seismicity and total spatial rate, clustering coefficientand conditional intensity function at the end of study period.

Ifplot.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 attimet=t_{\mathrm{start}}

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

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.InCase 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 bythe 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 astransformed times for a fitted ETAS model.

Usage

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

Arguments

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

forI_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 obtainwithh=1,h=1/\lambda andh=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 againstj and Q-Q plot ofU_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 Jalilianjalilian@razi.ac.ir

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 forPoint 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 andApplications in Model EvaluationJournal of the Royal Statistical Society: Series B (StatisticalMethodology),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)

Arguments

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 onhttps://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 earthquakecatalog dataset.

Author(s)

Abdollah Jalilianjalilian@razi.ac.ir

References

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

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

International Seismological Centre (2015).EHB Bulletin,https://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)