ETAS : An R Package for Fitting the Space-Time ETAS Model to Earthquake Data

The epidemic-type aftershock sequence (ETAS) model is the most widely used statistical model to describe earthquake catalogs. ETAS is an R package for ﬁtting the space-time ETAS model to an earthquake catalog using the stochastic declustering approach introduced by Zhuang, Ogata, and Vere-Jones (2002). The package provides two classes and several functions to facilitate data preparation, model ﬁtting and some simple diagnostic checks. The present paper is a description of the package and illustrates modeling earthquake data using the space-time ETAS model.


Introduction
Since the establishment of seismograph networks for instrumental monitoring of earthquakes in the past decades, earthquake catalogs in computer-readable formats are regularly accessible from national and international seismological agencies or institutions (Di Giacomo, Harris, Villaseñor, Storchak, Engdahl, and Lee 2015).An earthquake catalog is a chronologically ordered list of historic and/or instrumental earthquake records in a particular geographical region S and during a specific time period.It often consists of time, coordinates of epicenter, magnitude and focal depth of all recorded earthquakes with magnitudes greater than or equal to a certain threshold m 0 that occurred inside or in the vicinity of S during the specified time period.For example, Table 1 shows parts of an earthquake catalog covering the study period from 1973-01-01 to 2016-01-01 and the rectangular geographical region 22 • -42 • N and 40 • -65 • E (Iran and surrounding areas) with the magnitude threshold m 0 = 4.0.This catalog is described in more detail in Section 6 and is used for illustration purposes in this paper.
From a statistical point of view, an earthquake catalog is the available data on earthquakes in a space-time window (Ogata 1999) model for the underlying earthquake process.Such a model provides a theoretical description of the seismic activities in the study region and makes it possible to estimate the probabilities of future events (Kagan 1991;Vere-Jones 2003).Therefore, fitting an appropriate statistical model to a given earthquake catalog is of great importance in the probabilistic assessment of seismic hazard.
Among different proposed models, the epidemic-type aftershock sequence (ETAS) model, introduced by Ogata (1988) and Ogata (1998), is the most widely used statistical model to describe earthquake catalogs (see e.g., Helmstetter and Sornette 2003;Ogata 2006a;Console, Jackson, and Kagan 2010;Lombardi and Marzocchi 2010;Zhuang 2011;Harte 2013;Zhuang 2012).The space-time version of the ETAS model is a semi-parametric model that describes the background and triggering seismic activities in a geographical region and can be used for earthquake declustering (Zhuang et al. 2002).However, estimation of the ETAS model parameters is computationally challenging (Veen and Schoenberg 2008;Schoenberg 2013;Adelfio and Chiodi 2015a).
In this article we introduce the R (R Core Team 2018b) package ETAS (Jalilian and Zhuang 2019) which is an R implementation (through C and C++ ports) of the original Fortran implementation etas8p (Zhuang, Ogata et al. 2017).The package is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=ETAS and GitHub at https://github.com/jalilian/ETASand can be distributed under the GPL-2 license (GNU General Public License).
The material in the remainder of the paper is organized as follows.The ETAS model is described in Section 2, and Section 3 introduces the simultaneous estimation of model parameters and the background seismicity rate using a stochastic declustering approach.Residual analysis and model assessment are briefly discussed in Section 4. In Section 5, the main functions in the ETAS package are discussed and Section 6 provides a worked example of fitting the ETAS model to earthquake data based on a data set included in the package.The ETAS package is compared with other available software packages in Section 7. Further possible extensions of the package are discussed in Section 8.

Temporal marked point processes
Given an earthquake catalog consisting of N events, let t i , x i , y i and m i , i = 1, . . ., N , denote, respectively, the time, longitude of epicenter, latitude of epicenter and magnitude of the ith event in the catalog.Then the catalog data (as for example exemplified in Table 2) can be regarded as a point pattern (t i , x i , y i , m i ) : i = 1, . . ., N on R + × R 2 × [m 0 , ∞) which is generated by a temporal marked point process X that governs the occurrence of earthquakes in time and space (Vere-Jones 1970;Ogata 1998).
The distribution of the temporal marked point process X is uniquely determined by its conditional intensity function λ(t, x, y, m|H t ) (see Daley and Vere-Jones 2003, Section 7.3).Roughly speaking, for each t > 0, (x, y) ∈ S, m ≥ m 0 and infinitesimal dt, dx, dy and dm, the conditional intensity function λ(t, x, y, m|H t ) satisfies where is the occurrence history of the earthquakes up to time t (Ogata 1998;Zhuang et al. 2002).
We assume that the time t is measured in decimal days with t = 0 corresponding to a specific time origin, say 1945-01-01, 00:00:00, and a study period of length T days started at time t start ≥ 0 is of interest.Only events inside the study region S and the study period [t start , t start + T ] are considered as target events.Other events with t < t start or (x, y) ∈ S, if any, are assumed to be complementary events (Zhuang, Ogata, and Vere-Jones 2006).
Complementary events are included in the history H t in order to take into account their effect on target events and to some extent resolve what is known as the boundary (edge) effect (Ogata 1998).Note that N is the total number of target and complementary events in the catalog.

The space-time ETAS model
The ETAS model is a multiple type branching process and a special case of the marked Hawkes process (see Daley and Vere-Jones 2003, pp. 202-205). where is the probability density function (PDF) of the magnitude of an event and (Ogata 1998) λ θ (t, x, y|H t ) = ũ(x, y) with the following components: ũ(x, y) is the background seismicity rate and is assumed to be independent of time.The semi-parametric form ũ(x, y) = µu(x, y), where µ > 0 and u(x, y) is a smooth function on S, is usually considered for the background seismicity rate.
κ A,α (m i ) is the expected number of triggered events (aftershocks) generated from an event of magnitude m i and it can be expressed as where A > 0 and α > 0 are unknown parameters.
g c,p (t − t i ) is the PDF of the occurrence time of a triggered event generated from an event of magnitude m i occurring at time t i .It is assumed that the probability distribution of the time until the occurrence of a triggered event is a function of the time lag t − t i from its direct main shock and is independent of m i .Based on the modified Omori's law, the following inverse power density is considered where c > 0 and p > 1 are unknown parameters.
f D,γ,q (x − x i , y − y i ; m i ) is the PDF of the occurrence location of a triggered event generated from an event of magnitude m i occurring at the location (x i , y i ).It is assumed that the probability distribution of the location of a triggered event depends on the magnitude and the location of its direct main shock.The radially symmetric density function , is usually considered for this distribution, where D > 0, γ > 0 and q > 1 are unknown parameters.
The ETAS model assumes that each event is either a background (spontaneous) event or triggered by a previous event (Ogata 1998).The background events are generated by a Poisson process with intensity ũ(x, y), which is stationary in time.The relaxing coefficient coefficient µ in ũ(x, y) = µu(x, y) is introduced in order to fasten the convergence of the model fitting algorithm (see Algorithm 2 below).Previous events, whether they are background or triggered events, generate further events according to a non-stationary Poisson process with intensity function (Zhuang 2011) i:t i <t κ A,α (m i )g c,p (t − t i )f D,γ,q (x − x i , y − y i ; m i ).
The expected number of triggered events generated from a typical event, regardless of its magnitude, is If Aβ/(β − α) < 1, then the stability (existence of a stationary version) of the model is guaranteed (Daley and Vere-Jones 2003, p. 204).
Under the stability condition, the expected number of events in an interval of length T is proportional to T .Thus, the expected number of events tend to infinity as T → ∞.
Theoretically, the total spatial intensity function is defined as In practice, the total spatial intensity function is approximated by Λ(x, y) ≈ µu(x, y) + 1 T i:t i <T κ A,α (m i )f D,γ,q (x − x i , y − y i ; m i ) and used to obtain the clustering (triggering) coefficient ω(x, y) = 1 − u(x, y) Λ(x, y) , which quantifies the clustering effect relative to the total spatial intensity at any location (x, y) ∈ S (Zhuang et al. 2006).Zhuang et al. (2002) proposed a probabilistic approach for declustering an earthquake catalog to obtain the background seismicity rate.They suggested to use

Stochastic declustering
as the probability that event j is triggered by event i.Consequently, is the probability that event j is triggered by a previous event and 1 − p j = µu(x j , y j ) λ θ (t j , x j , y j |H t j ) is the probability that event j is a background event.

Parameter estimation and model fitting
The ETAS model, defined through (1) and ( 2), is a semi-parametric model with Euclidean parameters β and θ = (µ, A, α, c, p, D, γ, q) and an infinite dimensional parameter u(x, y), (x, y) ∈ S (Zhuang et al. 2006).The log-likelihood function of the model is given by (Ogata 1998) λ β,θ (t, x, y, m|H t )dxdydtdm, where δ i = 1 if event i is a target event and δ i = 0 otherwise.Due to separability of β and θ in (1) and the fact that ∞ m 0 ν β (m)dm = 1, the log-likelihood function can be written as where where tstart+T tstart S λ θ (t, x, y|H t )dxdydt = T µ S u(x, y)dxdy and Thus, by maximizing l 1 (β|H T ), the maximum likelihood estimate (MLE) of β is , where N ′ = N i=1 δ i , while to obtain the MLE of θ, an estimate of u(x, y), numerical approximations of the integral terms over S in (6) and an appropriate numerical maximization method for l 2 (θ|H T ) are required.Ogata (1998) suggested a radial partitioning of the geographical region S in order to efficiently approximate the integral term S f D,γ,q (x − x i , y − y i ; m i )dxdy = S (i)   f D,γ,q (x, y; m i )dxdy, where S (i) = (x − x i , y − y i ) : (x, y) ∈ S , i = 1, . . ., N.
In this method, for each event i, S (i) is partitioned into n v subregions S (i) 1 , . . ., S (i) nv by dividing the boundary of S (i) into n v knots and using radial segments to connect the origin, (0, 0), and each of these knots.Then the integral on each subregion is transformed from the Cartesian coordinates (x, y) to the polar coordinates (r, ϑ) and approximated using the length of the radial segments and the angles between them (see Figure 1).As follows in (10), the integral of u(x, y) over S is approximated similarly.Zhuang et al. (2002) introduced an iterative approach to simultaneous estimation of the smoothed function u(x, y), (x, y) ∈ S, and the parameter vector θ under the stochastic declustering framework.Given an initial estimate for u(u, v) and using the above mentioned integral approximation, the MLE of θ is obtained by minimizing ξ(θ) = −l 2 (θ|H T ) using the Davidon-Fletcher-Powell method, which is a gradient-based nonlinear optimization algorithm (Ogata 1998).Let ∇ξ(θ) denote the gradient of ξ(θ).Then Algorithm 1 summarizes the Davidon-Fletcher-Powell (DFP) procedure.As k → ∞ in Algorithm 1, θ k → θ * , where θ * is a local minimum of ξ(θ), and the matrix H k converges to the inverse of the Hessian matrix of ξ(θ) at θ = θ * ; i.e., 10: ⊲ convergence criteria 13: return θ k and H k (Fletcher and Powell 1963).Thus, when the algorithm converges, H k is an estimate (approximation) for the inverse of the Fisher information matrix I( θ k ), where Therefore, by asymptotic properties (T → ∞) of maximum likelihood estimators (Schoenberg 2005;Rathbun 1996), I( θ k ) = H k is an estimate of the asymptotic covariance matrix of θ.
Once the estimates of u(x, y) and θ are available, the declustering probabilities p j can be estimated from (3) and ( 4) and the variable bandwidth kernel estimator provides a better estimate for u(x, y), where is the bivariate isotropic Gaussian kernel function with bandwidth h.The bandwidth h j is considered to be where h min is a minimum threshold bandwidth value and r(j, n p ) denotes the distance between the location of event j and its n p th nearest neighbor.The n p th nearest neighbor distance r(j, n p ) can be very small or even zero because events may overlap at the same location.Thus, h min is used to prevent very small or even zero h i values (Zhuang et al. 2002).Since spatial point patterns of locations of earthquake epicenters are clustered, using a variable bandwidth in the kernel estimator ( 7) is more preferable than a fixed bandwidth.A fixed bandwidth Algorithm 2 Fitting the ETAS model using stochastic declustering method.1: Inputs: θ 0 ⊲ initial guess for the model parameters u 0 (x i , y i ), i = 1, . . ., N ⊲ initial background rate 2: H 0 ← I ⊲ initial guess for the inverse of the Hessian matrix 3: k ← 0 4: repeat 5: use Algorithm 1 with inputs θ k , u k (x i , y i ) and H k and get θ k+1 and H k+1 6: 11: ⊲ relative change in the log-likelihood function 12: until e 1 < ε and e 2 < ε and e 3 < ε ⊲ convergence criteria 13: return θ k , u k (x, y) and H k produces under-smoothing in areas with sparse observations and over-smoothing in areas with dense observations.The bandwidth ( 9) is adaptive in the sense that r(j, n p ), and hence h j , is small in areas with dense observations and large in areas with sparse observations.As suggested by Zhuang (2011), n p = 3, . . ., 6 seems reasonable and h min should be equivalent to the order of earthquake location error (due to location inaccuracy or rounded coordinates at a certain precision), say h min = 0.05 • equivalent to 5.56 kilometers on the earth surface.(1 − p j ) S (j)   ϕ(x, y; h j )dxdy (10) and since ( 8) is a radially symmetric function, integral terms in the right hand side of (10) can be approximated by transforming to the polar coordinate system and using the radial partitioning of the geographical regions S (j) as described above.Replacing the initial estimate of u(x, y) by the new estimate from (7), estimating θ and u(x, y) is repeated until the estimates converge (Zhuang et al. 2006).This iterative approach is summarized in Algorithm 2.
Finally, similar to (7), the total spatial intensity function can be estimated by which along with u(x, y) yields the estimate ω(x, y) = 1 − u(x, y)/ Λ(x, y) for the clustering coefficient.
The so-called transformed times for target events j ∈ J = {i : δ i = 1} are closely related to raw temporal residuals R temp ([t start , t j ]; 1) (Zhuang 2006).If the fitted model is close enough to the true underlying earthquake process, then {τ j : j ∈ J} follows a standard (unit rate) Poisson process and, consequently, U j = 1 − exp − (τ j − τ j−1 ) , j ∈ J, are independent and identically distributed uniform random variables on (0, 1) (Ogata 1988;Rathbun 1996).Thus, plotting τ j against j and a Q-Q plot of U j yield simple diagnostic checks for the temporal adequacy of the fitted model: If the model is good enough then points of the first plot approximately lie on the line y = x and the second plot shows agreement between the empirical quantiles of U j and the theoretical quantiles of a U (0, 1) distribution.Furthermore, the Kolmogorov-Smirnov test for goodness-of-fit can be employed to determine whether random variables U 2 , . . ., U N follow a U (0, 1) distribution.

The ETAS package
Assume that the catalog data are imported into R as a 'data.frame'object (see the R Data Import/Export manual for information about importing and working with different data formats in R; R Core Team 2018a).This 'data.frame'object is required to have at least 5 columns with (exact or partial match of) names date, time, lat, long and mag containing, respectively, the date, time, latitude, longitude and magnitude of each event in the catalog.
The date and time columns should be in appropriate formats such that paste(date, time) can be converted to date-time (calendar dates plus time to the nearest second) by the R function as.POSIXlt.For information on calendar dates and times in R, see the help on "DateTimeClasses" obtained with: R> help("DateTimeClasses")

Data preparation with function catalog
As the first step, the 'data.frame'object that contains the earthquake data should be passed to the user accessible catalog function in order to specify the data as an earthquake catalog.
The function catalog creates an object of class 'catalog' to represent an earthquake catalog.An object of this class contains the earthquake data as well as the intended study period, geographical region and the magnitude threshold.The magnitude threshold can be specified by the argument mag.threshold.If not specified, the minimum magnitude value in the catalog will be used as the magnitude threshold.• time.begin is the beginning of time span of the catalog.By default, it is set to the time of the first event in the catalog.
• study.start is the start of the study period.If not specified, then time.begin will be used instead.
• study.end is the end of the study period.By default, it is set to the time of the last even in the catalog.
• study.length is the length of the study period (T ) in decimal days.
Arguments time.begin,study.startand study.endmust be character strings or objects that can be converted to date-time by the as.POSIXlt function.Either study.endor study.lengthcan be specified, but not both.The occurrence times are converted to elapsed decimal days since time.begin.If the events in the catalog are not chronologically sorted, then the catalog function will produce a warning and sort the events in ascending order with respect to time of occurrence.
In ETAS, a geographical region is defined as a set of point coordinates in a rectangular or irregular polygonal two-dimensional study region (see Figure 3).For a rectangular geographical region, the latitude and longitude ranges should be specified by lat.range and long.rangearguments.A polygonal study region can be specified by the region.polyargument which contains coordinates of the vertices of the polygon.The region.polyargument 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., the first vertex is not repeated, see Figure 3).(study.start, study.end)are considered as target events and other events are complementary events.
By default (flatmap = TRUE) the spherical coordinates (long, lat) on the earth surface are transformed 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.The default case dist.unit= "degree" uses the Equirectangular projection (see Snyder 1997, pp. 5-8) x = cos(cnt.lat/180× pi) × (long − cnt.long) and y = lat − cnt.lat,where here x and y are in (approximate) kilometers.Note that changing dist.unitaffects the unit of u(x, y) and D parameters (see Table 3).
Specific print and plot methods associated with 'catalog' objects are designed to allow users to see a summary of the specified earthquake catalog.

Fitting the model with function etas
The function etas fits the ETAS model to a catalog of earthquakes.The function takes an object of class 'catalog' and calls the internal functions decluster and etasfit for simultaneous estimation of the background seismicity rate u(x, y) and the parameter vector θ using the iterative approach discussed in Section 3.An initial guess for the parameter vector θ = (µ, A, α, c, p, D, γ, q) can be provided by the param0 argument.Thus, param0 needs to be a numeric vector of length 8.If param0 is not specified, as suggested by Ogata (1998), the default values µ = N/(4T |S|), A = 0.01, c = 0.01, α = 1, p = 1.3, D = 0.01, q = 2 and γ = 1 are used as components of param0.Note that this is a crude initial estimate and does not guarantee the convergence of Algorithm 2. The algorithm is relatively sensitive to the choice of the initial estimate of θ and, if possible, providing more informative initial estimate is highly recommended.
The bandwidth for each event in the kernel estimators ( 7) and ( 11) of the background seismicity rate and total spatial intensity can be specified by the bwd argument, which needs to be a numeric vector of length N .If bwd is not specified, the nnp and bwm arguments determine, respectively, the number of nearest neighbors n p and the minimum threshold bandwidth h min and ( 9) is used to obtain the bandwidth for each event.The values nnp = 5 and bwm = 0.05 are used as default values for nnp and bwm (see Zhuang 2011).The argument ndiv, with default value ndiv = 1000, sets the number of knots on each side of the geographical region S in radial partitioning of S for integral approximation as described in Section 3.
The no.itr argument specifies the maximum number of iterations in the iterative estimation method (Algorithm 2).The estimates often converge in less than ten iterations.If no.itr is not specified, no.itr = 11 will be used.Moreover, the relative iteration convergence tolerance ε in Algorithm 2 and the optimization convergence tolerance ǫ in Algorithm 1 are, respectively, determined by the rel.tol and eps arguments.The progress of the computations can be traced by setting the the verbose and plot.itarguments to TRUE.For example, the following output snippet shows the value of the target function ξ(θ) = −l 2 (θ|H T ) and its gradient ∇ξ( θ) for the current estimate of θ and the value of the optimal coefficient ζ of the line search procedure in one iteration of Algorithm 1.The etas function returns an object of class 'etas'.Specific print and plot methods associated with 'etas' objects are designed to allow users to see a summary of the fitted model.An object of class 'etas' can be passed to the rates function to obtain and plot estimates of the total spatial intensity, background seismicity rate and clustering coefficient from the fitted model over a grid on S. The probs and resid.etasfunctions take an object of class 'etas' and return estimates of declustering probabilities p j and model residuals and diagnostic plots, respectively.
Lack of completeness and time-non-stationarity of earthquake catalogs produce some problems in the statistical analysis of the catalog data.This version of the ETAS model assumes that the earthquake catalog is complete and the data are stationary in time (Zhuang 2011).If the catalog is incomplete or there is non-stationarity (e.g., increasing or seasonal trend) in the time of events, then the results of this function are not reliable.2011).Thus, parallel computing of the computations of the log-likelihood function (5) reduces the computation time for large earthquake catalogs.
By default (cxxcode = TRUE) the C++ code is used and the argument nthreads 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.As pointed out in the Writing R Extensions manual (R Core Team 2018c): The performance of OpenMP varies substantially between platforms.The Windows implementation has substantial overheads, so is only beneficial if quite substantial tasks are run in parallel.
Thus, in particular on Windows, the parallel computing (nthreads > 1) should only be considered for catalogs with a large number of events (say more than 2000).Also, the number of threads nthreads should be carefully specified.The detectCores function in the parallel package (R Core Team 2018b) can be consulted to find out the overall number of available threads on a given machine.For example, on a laptop with an Intel Core i7 processor we get R> parallel::detectCores() [1] 8 However, resource usage and limitations should be considered when setting nthreads.

Data example
The following example shows how to fit the space-time ETAS model to an earthquake catalog using the R package ETAS.Three sample earthquake catalog data sets iran.quakes,italy.quakesand japan.quakesare included within the package.They are objects of class 'data.frame'and can be directly accessed by loading the package.

R> library("ETAS")
The iran.quakes data set contains data of an earthquake catalog of Iran and surrounding areas (the rectangular geographical region 22 • -42 For the purpose of illustration, we consider the study period from 1991-01-01 to 2011-01-01 (t start = 6574 and T = 7305 days), the magnitude threshold m 0 = 4.5 and a polygonal geographical region with the vertices at (26, 52), (25, 59), (29, 58), (38, 45) and (35, 43).Thus, all events in iran.quakes with magnitude greater than or equal to 4.5 that occurred between 1991-01-01 and 2011-01-01, and inside the specified geographical region, are target events.Events with magnitude greater than or equal to 4.5 that occurred before 1991-01-01 or outside the specified geographical region act as complementary events.Note that the start date 1991-01-01 of the study period is t start = 6574 days after the beginning date 1973-01-01 of the time span of the catalog and that the length of the study period is T = 7305 days.

Time difference of 7305 days
To create a catalog with the specified study period, geographical region and magnitude threshold, we use the following commands.

R> plot(iran.cat)
It produces a plot consisting of six subplots: a plot to show the spatial distribution of events and the selected region, three plots to see how latitude, longitude and magnitude of events change over time and two plots to visually inspect the completeness and time stationarity of the catalog.Let N m denote the number of events in the catalog with magnitude greater than or equal to m.The Gutenberg-Richter law states that log 10 (N m ) = a − bm for some a and b, which is equivalent to assuming that m follows an exponential distribution (Ogata 1998).Departure from linearity in the plot of log 10 (N m ) versus m is an indication of the catalog incompleteness.In addition, the minimum m value for which this linear relation holds can be used as the magnitude threshold.Time stationarity of the catalog, on the other hand, can be assessed by inspection of the plot of N t versus t, where N t is the number of events in the catalog that occurred up to t.If the earthquake process is stationary in time, then we expect N t = λ 0 t for some λ 0 > 0.
The resulting plot is shown in Figure 4.The plot of log 10 (N m ) versus m shows the linear relation expected from the Gutenberg-Richter law.The plot of N t versus t is linear during the considered study period, but a nonlinear trend is evident over the time span 1973-01-01 to 2011-01-01 (13879 days).Nevertheless, we assume that the specified catalog is complete and stationary in time.
To fit the ETAS model to the created catalog iran.cat,we first provide an initial guess for the model parameters.
0 4000 8000 12000 4.5 5.0 5.5 6.0 The output iran.fit is an object of class 'etas'.The maximum likelihood estimates of the model parameters, their estimated asymptotic standard errors, the value of the log-likelihood function at its maximum point and the Akaike information criterion (AIC) of the model can be seen by printing iran.fit.The estimated asymptotic standard errors of c and specially D and q are much larger than their corresponding estimates.This may be caused by insufficient data (small catalog) or model inadequacy, as discussed in what follows.

R>
Plot of iran.fitdisplays the estimated parameters in each iteration.

R> plot(iran.fit)
Figure 5 shows the plot.It can be seen that the estimates converge after four iterations.Estimates of the background seismicity rate u(x, y), total spatial intensity Λ(x, y), clustering coefficient ω(x, y) and conditional intensity function λ θ (x, y, t|H t ) at t = t start + T can be obtain by passing iran.fit to the rates function.

R> rates(iran.fit)
This function computes and plots the estimates on a rectangular grid over the geographical region (see Figure 6).The grid span and dimension can be changed by arguments lat.range, long.range and dimyx.
In addition, the function probs returns estimates of the declustering probabilities p j .

R> pr <-probs(iran.fit) R> summary(pr$prob)
Min. 1st Qu.Median Mean 3rd Qu.Max.0.00000 0.02299 0.09680 0.36346 0.89058 1.00000 These probabilities may be used to indicate the most likely background and triggered target events (see Figure 7).+ pch = 3, col = 2) R> map("world", add = TRUE, col = "grey") R> legend ("bottomleft", c("background", "triggered") The first-order residual analysis for a fitted ETAS model can be performed by using the resid.etasfunction.The resid.etas function computes the temporal residuals R temp (I j ; h), where I j = t start + (j − 1)T /n temp , t start + jT /n temp , j = 1, . . ., n temp .The argument n.temp with default value n.temp = 1000 specifies the number of partition points n temp .For spatial residuals, the spatstat function quadscheme is used to create a Berman-Turner quadrature over S and compute R spat (B i ; h), as described in Section 4. The resid.etas function by default (type = "raw") returns raw temporal residuals, raw spatial residuals, transformed times τ j and their related quantities U j .It also produces plots of the temporal residuals, smoothed spatial residuals, transformed times τ j against j and the Q-Q plot of U j .Figure 8 shows the resulting plots for the fitted ETAS model to the Iranian catalog.It can be seen that points on the plots of transformed times τ j and the Q-Q plot of U j lie approximately on the y = x line.Likewise, the Kolmogorov-Smirnov test for goodness-of-fit indicates that U j follows a U (0, 1) distribution.

R> ks.test(iran.res$U, punif)
One-sample Kolmogorov-Smirnov test data: iran.res$UD = 0.031565, p-value = 0.4938 alternative hypothesis: two-sided Therefore, the fitted model is able to adequately describe the temporal variation in the data.However, a systematic northwest-southeast trend is evident in the raw spatial residuals.The large (inflated) estimated asymptotic standard errors of D and q perhaps are another sign that the fitted model is not capable of describing the complex spatial clustering of earthquakes in this region.In fact, the geographical region considered here contains the Zagros fold-andthrust belt, which due to the ongoing collision between the Arabian Plate and the Eurasian Plate, is the most seismically active intra-continental belts on Earth (Talebian and Jackson 2004).Hence, the ETAS model with fixed parameters may not be suitable to model the occurrence of earthquakes in such an active tectonic region.Although this earthquake catalog provides as an example for the ETAS package, finding an appropriate model for it is well beyond the scope of this paper.It is possible that improved versions of the ETAS model with spatially varying parameters (Ogata 2011) may be more suitable for this region.

Other software packages
The temporal ETAS model can be fitted to the occurrence time and magnitude of earthquakes using the Fortran programs SASeis2006 (Ogata 2006b) and etas_solve (Kasahara, Yagi, and Enescu 2016) and the R packages PtProcess (Harte 2010) and SAPP (The Institute of Statistical Mathematics 2016).To fit the space-time ETAS model to the occurrence time, location of epicenter and magnitude of earthquakes, the R package etasFLP (Chiodi and Adelfio 2017) and the original Fortran code by Jiancang Zhuang (Zhuang et al. 2017) are publicly available.

Comparison with the etasFLP package
The R package etasFLP is an alternative to ETAS and fits both the temporal and space-time ETAS models.The etasFLP package uses the forward likelihood predictive (FLP) method (Chiodi and Adelfio 2011) to estimate the bandwidth of the kernel density estimator (7) for  the background seismicity rate (Adelfio andChiodi 2013, 2015a).Also, instead of (2), the conditional intensity function in the etasFLP package is parameterized as The package comes with a sample catalog of 2158 Italian earthquakes of magnitude greater than or equal to 3.0 that occurred between 2005 and 2013, extracted from the Italian Seismological Instrumental and Parametric Data Base (ISIDE) available at http://iside.rm.ingv.it/.
R> k0 <-0.005R> c <-0.005 R> p <-1.01 R> D <-1.1 R> q <-1.52 R> A <-pi * k0 / ((p -1) * c^(p -1) * (q -1) * D^(q -1)) R> param0 <-c(mu = 1, A, c, alpha = 1.05, p, D, q, gamma = 0.6)Although the estimates obtained with both packages etasFLP and ETAS are to some extent similar (particularly the estimates of α, p and q), there are some differences between them (particularly the estimates of µ and D) due to the different bandwidth selection methods for kernel estimator of the background seismicity rate, different optimization algorithms and different transformations of (long, lat) to (x, y) coordinates.The p value indicates an adequate temporal fit to the catalog.

Comparison with the Fortran program
The original Fortran program for fitting the space-time ETAS model to earthquake data is available at http://bemlar.ism.ac.jp/zhuang/software/sd.tar.gz(accessed on November 3, 2018).The main program etas8p calls several subroutines for integral approximation by radial partitioning of the geographical region (frint and polyint), computing the likelihood function (xlamb and xint), declustering computations (bkgdcalc) and Davidon-Fletcher-Powell optimization (linear and davidn).The program estimates the model parameters and background seismicity rate using the stochastic declustering approach as in Algorithm 2, except that it does not check the convergence criteria in the repeat loop and uses a fixed number of iterations (eleven times).Moreover, in each iteration, the Davidon-Fletcher-Powell optimization described in Algorithm 1 (subroutine davidn) starts with the identity matrix, H 0 = I, as the initial guess for the inverse of Hessian matrix while in the ETAS package, the estimate of inverse Hessian matrix in each iteration is used as initial guess for the next iteration.These differences generally result in a faster convergence of ETAS.After convergence, etas8p calls two subroutines (outprob and outrates) to report declustering probabilities p j and the values of the estimated background seismicity rate on a rectangular grid.-log likelihood = 0.153109518614449D+05 aic = 30637.9-----x -----0.74203D+000.40714D+00 0.17209D+00 0.12876D+01 0.10740D+01 0.42823D-01 0.13966D+01 0.10330D+01 *** gradient *** -0.71766D-03 0.16731D-05 0.36575D-03 -0.53808D-03 0.14690D-03 -0.11811D-02 0.16704D-03 0.13040D-03 mle = 0.55061E+00 0.16576E+00 0.29615E-01 0.16579E+01 0.11534E+01 0.18338E-02 0.19505E+01 0.10670E+01 The program, in addition, returns estimates of the background seismicity rate, the total spatial rate and the conditional intensity function at t = t start + T on a 401 × 401 grid over 31.8 • -44.3 • N and 134 • -145 • E (rates.datfile) and the background seismicity rate, declustering probability and bandwidth for each target event (file probs.dat).Figure 11 shows the estimated background seismicity rate, imported to and plotted in R.
For comparison, the data set of the Japanese earthquakes is included in ETAS as japan.quakes.

R>
The maximum absolute difference between the estimates of the model parameters obtained with the etas8p program and ETAS package is smaller than 10 −3 .The absolute difference between the values of the log-likelihood function and AIC are less than 10 −2 .The estimate of the background seismicity rate, obtained with the command (elapsed execution time 5.7 minutes) R> rates (japan.fit, long.range = c(134, 145), lat.range = c(31.8, 44.3), + dimyx = c(401, 401)) and shown in the top left panel of Figure 13, is also in total agreement with Figure 11.Thus, as expected, results obtained with the etas8p program and ETAS are practically the same.
For completeness, we perform residual analysis of the fitted ETAS model.The resulting diagnostic plots are shown in Figure 14.Plot of transformed times τ j and the Q-Q plot of U j as well as the Kolmogorov-Smirnov test for goodness-of-fit indicate that the temporal variations in the data are sufficiently explained by the fitted model.

Conclusion
This paper introduces the R package ETAS that fits the epidemic type aftershock sequence (ETAS) model to an earthquake catalog.Given data on an earthquake catalog, the package provides two classes and several functions to facilitate data preparation and model fitting and offers some simple residual analysis and diagnostic checks for model assessment.
The ETAS package is based on a Fortran program by Jiancang Zhuang and colleagues.The results obtained with the ETAS package and the Fortran program are practically the same, but the ETAS package achieves a faster execution time by using the estimated inverse Hessian matrix at each iteration as an initial guess for the next iteration, and using modified convergence criteria.The R package etasFLP is an alternative to ETAS.These packages produce different, although comparable, results because they use different bandwidth selection methods for the kernel estimator of the background seismicity rate, different optimization algorithms and different flat map transformations.The etasFLP package is faster than ETAS with serial computing, but for large catalogs (more than 2000 events) the ETAS package with parallel computing can be as fast as etasFLP.
As of version 0.3, the ETAS package uses parallel computing with OpenMP to allow for faster computations, particularly in case of large earthquake catalogs.The package can be extended in several ways.For example, different modifications of the space-time ETAS model have been recently proposed (Ogata 2011;Harte 2014;Nicolis, Chiodi, and Adelfio 2015).Moreover, the ETAS model has been generalized to introduce nonstationarity into the ETAS model (Kumazawa and Ogata 2014).These modifications and generalizations could be implemented and integrated into the package.Simulation of the space-time ETAS model and more comprehensive residual analyses and model diagnostics could also be included in the package.

Figure 1 :
Figure 1: Radial partitioning of the geographical region S for integral approximation.

Figure 3 :
Figure 3: Specifying the rectangular or polygonal two-dimensional study region in the catalog function using lat.range and long.rangeor region.polyarguments.

Figure 4 :
Figure 4: Location of epicenters (top left panel), logarithm of frequency by magnitude (bottom left panel), cumulative frequency over time (bottom middle panel) and latitude, longitude and magnitude against time (right panels) of 5970 earthquakes with magnitude greater than or equal to 4.5 occurred between 1973-01-01 and 2016-01-01 in Iran and its vicinity (22 • -42 • N and 40 • -65 • E), extracted from the ANSS global comprehensive catalog (available at http://earthquake.usgs.gov/earthquakes/search/). The plots are created by calling plot on an object of class 'catalog'; an earthquake catalog created by the catalog function.

Figure 6 :
Figure 6: Plots of estimates of the background seismicity rate, total spatial intensity, clustering (triggering) coefficient and conditional intensity at t = t start + T for the fitted ETAS model to the Iranian catalog.The plots are created by calling the rates function on an object of class 'etas'.

Figure 8 :
Figure 8: Diagnostic plots for the fitted ETAS model to the Iranian catalog: the temporal residuals (top left), smoothed spatial residuals (top right), transformed times τ i against i (bottom left) and the Q-Q plot of U i (bottom right).The plots are created by calling the resid.etasfunction on an object of class 'etas'.

Figure 9 :
Figure 9: Location of epicenters (top left panel), logarithm of frequency by magnitude (bottom left panel), cumulative frequency over time (bottom middle panel) and latitude, longitude and magnitude against time (right panels) of 2158 earthquakes with magnitude greater than or equal to 3.0 that occurred between 2005-04-16 and 2013-11-01 in Italy and its vicinity (35 • -48 • N and 6.15 • -19 • E), extracted from the Italian Seismological Instrumental and Parametric Data Base (ISIDE) available at http://iside.rm.ingv.it/.The catalog is also available in the R package etasFLP.

Figure 10 :
Figure 10: Diagnostic plots for the fitted ETAS model to the Italian catalog in the etasFLP package using the ETAS package: the temporal residuals (top left), smoothed spatial residuals (top right), transformed times τ i against i (bottom left) and the Q-Q plot of U i (bottom right).

Figure 11 :
Figure 11: Estimate of the background seismicity rate for the Japanese catalog using the etas8p program.

Figure 12 :
Figure 12: Location of epicenters (top left panel), logarithm of frequency by magnitude (bottom left panel), cumulative frequency over time (bottom middle panel) and latitude, longitude and magnitude against time (right panels) of 13724 earthquakes with magnitude greater than or equal to 4.5 that occurred between 1926-01-08 and 2007-12-29 in Japan and its vicinity (27 • -45 • N and 128 • -145 • E), accompanying the Fortran program by Jiancang Zhuang (Zhuang et al. 2017).

Figure 13 :
Figure 13: Plots of estimates of the background seismicity rate, total spatial intensity, clustering (triggering) coefficient and conditional intensity at t = t start + T for the fitted ETAS model to the Japanese catalog accompanying the Fortran program etas8p using the ETAS package.

Figure 14 :
Figure 14: Diagnostic plots for the fitted ETAS model to the sample Japanese catalog accompanying the Fortran program etas8p using the ETAS package: the temporal residuals (top left), smoothed spatial residuals (top right), transformed times τ i against i (bottom left) and the Q-Q plot of U i (bottom right).

Table 1 :
. Statistical analysis can then be used to find a suitable 49.34A catalog of earthquakes for Iran and surrounding areas from the beginning of 1973 to the end of 2015.

Table 2 :
The conditional intensity function of Example of catalog data.the space-time ETAS model is given by λ β,θ (t, x, y, m|H t ) = ν β (m)λ θ (t, x, y|H t ), Specifying the study period using time.begin,study.start,study.endand study.lengtharguments in the catalog function.
The arguments time.begin,study.start,study.endandstudy.lengthmaybeused to determine the study period (see Figure2):

Table 3 :
Units of the ETAS model parameters.Dimensionless parameters are denoted by -.
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.Setting dist.unit= "km" performs the projection x = 111.32× cos(lat/180 × pi) × long and y = 110.547× lat, The catalog can be printed by simply typing its name.
The maximum number of iterations for the Algorithm 2 can be specified by the argument ndeclust.
• -48 • N and 6.15 • -19 • E), extracted from the Italian Seismological Instrumental and Parametric Data Base (ISIDE) available at http://iside.rm.ingv.it/.The catalog is also available in the R package etasFLP.The data set italy.quakes in the ETAS package contains data of the same Italian catalog.For compatibility with the etasFLP package, we use the argument dist.unit= "km" in the catalog function to create the earthquake catalog.
The catalog and initial parameters are passed to the etas function to fit the space-time ETAS model.the etasFLP package and 31 minutes for the ETAS package.Using parallel computing in the ETAS package, the elapsed execution time reduces to 15 minutes with nthreads = 4 and 10 minutes with nthreads = 8.
(Adelfio and Chiodi 2015bS model by the ETAS package is 46805.03which is smaller than the one by etasFLP (48414.86).Figure10shows the diagnostic plots for the fitted model by the ETAS package.In ks.test(italy.res$U,punif):tiesshouldnot be present for the Kolmogorov-Smirnov testExcept for a few extreme temporal residuals, the plot of the transformed times, the Q-Q plot of the U j values and the Kolmogorov-Smirnov test for goodness-of-fit indicate a good temporal fit.No systematic pattern is evident in the spatial residuals, but they are not symmetric around zero and a few small clusters of extreme residuals on the image of spatial residuals may suggest room for improvements in the fitted PDF of spatial distribution of triggered events, f D,γ,q , or in estimating the background seismicity rate.The etasFLP package also computes and visualizes transformed times τ i , background and triggered seismicity rates, total spatial intensity as well as different types of spatial residuals(Adelfio and Chiodi 2015b).The following commands compute U j values and perform the Kolmogorov-Smirnov test for temporal goodness-of-fit of the fitted ETAS model to the Italian catalog using the etasFLP package.
(The MPI Forum 2015)ses parallel computing, implemented by the message passing interface (MPI) standard(The MPI Forum 2015).The program's tarball sd.tar.gzcontains a sample data set (alljpM45.etasfile) of 13724 shallow (depth < 100 km) earthquakes that occurred between 1926-01-08 and 2007-12-29 in Japan and its vicinity (27 • -45 • N and 128 • -145 • E) with magnitudes greater than or equal to 4.5.As instructed on http://bemlar.ism.ac.The file jap.in contains all program inputs including the name of the data file, the beginning and length of the study period, the magnitude threshold, the coordinates of vertices of the boundary of the target geographical region, the initial values for the model parameters, the span and the dimension of the rectangular grid for estimation of the background seismicity rate and values of n p and h min for computing the bandwidth h j for each event according to (9).The program can be run on a computer with 8 threads processor byThe execution time takes 422.3 minutes on the same Linux laptop as before (2.4 GHz Intel Core i7-4700MQ CPU) with GNU Compiler Collection, Version 5.4.0 (https://gcc.gnu.org/), and MPICH, Version 3.2 (https://www.mpich.org/).Here is the final output of the program.