Description Usage Arguments Details Value Note Author(s) References See Also Examples

`etasclass`

is the main function of the package `etasFLP`

.

Performs the estimation of the components of the ETAS (Epidemic Type Aftershock Sequence) model for the description of the seismicity in a space-time region. Background seismicity is estimated non-parametrically, while triggered seismicity is estimated by MLE. In particular also the bandwidth for a kernel smoothing can be estimated through the Forward Likelihood Predictive (FLP) approach . For each event the probability of being a background event or a triggered one is estimated.

An ETAS with up to 8 parameters can be estimated, with several options and different methods.

Returns an `etasclass`

object, for which `plot`

, `summary`

, `print`

and `profile`

methods are defined.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 | ```
etasclass(cat.orig,
magn.threshold=2.5, magn.threshold.back=magn.threshold+2,
mu=1,k0=1,c=0.5,p=1.01, a=1.2,gamma=.5,d=1.,q=1.5, params.ind=replicate(8,TRUE),
hdef=c(1,1),
declustering=TRUE,thinning=FALSE,
flp=TRUE, m1=NULL, ndeclust=5, onlytime=FALSE,is.backconstant=FALSE,
w=replicate(nrow(cat.orig),1),
##### end of main input arguments.
##### Control and secondary arguments:
description="", cat.back=NULL, back.smooth=1.0,
sectoday=TRUE,longlat.to.km=TRUE,
usenlm=TRUE, method ="BFGS", compsqm=TRUE,
epsmax=0.0001, iterlim=100, ntheta=100)
``` |

`cat.orig` |
An earthquake catalog, possibly an object of class |

`magn.threshold` |
Threshold magnitude (only events with a magnitude at least |

`magn.threshold.back` |
Threshold magnitude used to build the catalog |

*Values for the 8 parameters of the ETAS model* (starting values or fixed values according to `params.ind`

):

`mu` |
Parameter 1 ( |

`k0` |
Parameter 2 ( |

`c` |
Parameter 3 of the ETAS model; a shift parameter of the Omori law for temporal decay rate of aftershocks; see details. Default value = 0.5. |

`p` |
Parameter 4 of the ETAS model; the exponent of the Omori law for temporal decay rate of aftershocks; see details. Default value = 1.01. |

`a` |
Parameter 5 ( |

`gamma` |
Parameter 6 ( |

`d` |
Parameter 7 of the ETAS model; parameter related to the spatial influence of the mainshock; see details. Default value = 1. |

`q` |
Parameter 8 of the ETAS model; parameter related to the spatial influence of the mainshock; see details. Default value = 1.5. |

*End of model pararameter input*

`params.ind` |
vector of 8 logical values: |

*Flags for the kind of declustering and smoothing*:

`hdef` |
Starting values for the |

`declustering` |
if |

`thinning` |
if |

`flp` |
if |

`m1` |
Used only if |

`ndeclust` |
maximum number of iterations for the general declustering procedure. Default=5. |

`onlytime` |
if |

`is.backconstant` |
if |

`w` |
initial weights |

*Other control parameters*:

`description` |
a description string used for the output. Default value = "". |

`cat.back` |
external catalog used for the estimation of the background seismicity.
Default value = |

`back.smooth` |
Controls the level of smoothing for the background seismicity (meaningful only if |

`sectoday` |
if |

`longlat.to.km` |
if |

`usenlm` |
if |

`method` |
used if |

`compsqm` |
if |

`epsmax` |
maximum allowed difference between estimates in subsequent iterations (default = 0.0001). |

`iterlim` |
maximum number of iterations in the maximum likelihood steps (used in |

`ntheta` |
number of subdivisions of the round angle, used in the approximation of the integral involved in the likelihood computation of the ETAS model. Default value = 100. |

Estimates the components of an ETAS (Epidemic type aftershock sequence) model for the description of the seismicity of a space-time region. Background seismicity is estimated nonparametrically, while triggered seismicity is estimated by MLE.

The bandwidth of the kernel density estimator is estimated through the
Forward Likelihood Predictive approach (FLP),
(theoretical reference on Adelfio and Chiodi, 2013) if `flp`

is set to `TRUE`

.
Otherwise the bandwidth is estimated trough Silverman's rule.
FLP steps for the estimation of nonparametric background component is alternated with the Maximum Likelihood step for the estimation
of parametric components (only if `declustering=TRUE`

).
For each event the probability of being a background event or a triggered one is estimated,
according to a declustering procedure in a way similar to the proposal of Zhuang, Ogata, and Vere-Jones (2002).

The ETAS model for conditional space time intensity *lambda(x,y,t)* is given by:

*
lambda(x,y,t)=mu*f(x,y)+
sum_(t_j<t)(k0\ e^(a-gamma) \ (m_j-m_0))/(t-t_j +c)^p
[ ((x-x_j)^2+(y-y_j)^2)/(e^(gamma \ (m_j-m_0))+d) ]^(-q)*

*f(x,y)* is estimated through a weighted kernel gaussian estimator; if `flp`

is set to `TRUE`

then the bandwidth is estimated through a FLP step.

Weights (computed only if `declustering=TRUE`

) are given by the estimated probabilities of being a background event; for the i-th event this is given by
*rho_i=(mu f(x_i,y_i))/(lambda(x_i,y_i,t_i))*.
The weights *rho_i* are updated after a whole iteration.

`mu`

(*mu*) measures the background general intensity (which is assumed temporally homogeneous);

`k0`

(*k_0*) is a scale parameter related to the importance of the induced seismicity;

`c`

and `p`

are the characteristic parameters of the seismic
activity of the given region; `c`

is a shift parameter while `p`

, which characterizes the pattern of seismicity, is the exponent parameter of the modified Omori law for temporal decay rate of aftershocks;

`a`

(*alpha*) and `gamma`

(*gamma*) measure the efficiency of an event of given magnitude in generating aftershock sequences;

`d`

and `q`

are two parameters related to the spatial influence of the mainshocks.

Many kinds of ETAS models can be estimated, managing some control input arguments.
The eight ETAS parameters can be fixed to some input value, or can be estimated, according to `params.ind`

:
if `params.ind`

[i]=FALSE the i-th parameter is kept fixed to its input value, otherwise, if `params.ind[i]`

` = TRUE`

, the i-th parameter is estimated and the input value is used as a starting value.

By default `params.ind=c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE)`

, and so a full 8 parameters ETAS model will be estimated.

The eight parameters are internally ordered in this way: `params`

= (`mu`

, `k0`

, `c`

, `p`

, `a`

, `gamma`

, `d`

, `q`

); for example
a model with a fixed value `p=1`

(and `params.ind`

[4] = FALSE) can be estimated and compared with the model where `p`

is estimated (`params.ind`

[4]=TRUE);

for example a 7 parameters model can be fitted with `gamma=0`

and `params.ind[6]`

=`FALSE`

,
so that input must be in this case: `params.ind=c(TRUE,TRUE,TRUE,TRUE,TRUE,FALSE,TRUE,TRUE)`

;

if `onlytime=TRUE`

a time process is fitted to data (with a maximum of 5 parameters), regardless to space location (however the input catalog `cat.orig`

must contain three columns named `long`

, `lat, z`

);

if `is.backconstant=TRUE`

a process (space-time or time) with a constant background intensity *mu* is fitted;

if `mu`

is fixed to a very low value a process with very low background intensity is fitted, that is with only clustered intensity (useful to fit a model to a single cluster of events).

If `flp=TRUE`

the bandwidth for the kernel estimation of the background intensity is evaluated maximizing the
Forward Likelihood Predictive (FLP) quantity, given by
(Chiodi, Adelfio, 2011; Adelfio, Chiodi, 2013):

*FLP_{k_1,k_2}(\hat{\boldsymbol{ψ}})\equiv∑_{k=k_1}^{n-1}δ_{k,{k+1}}(\hat{\boldsymbol{ψ}}(H_{t_k});
H_{t_{k+1}})*

with *k_1=\frac{n}{2},k_2=n-1* and where *(d_(k,k+1))(psi(Ht);H(t+1))* is
the *predictive information*
of the first *k* observations on the *k+1*-th observation, and is so defined:

*d_(k,k+1)=log L(psi^(H_k);H_(k+1)-log L(psi^(H_k);H_k)*

where *H_k* is the history of the process until time *t_k* and
*psi^(H_{t_k})* is an estimate based only on history until the *k-th* observation.

In the ML step, the vector of parameter *θ=(μ, κ_0, c , p, α, γ, d, q)* is estimated maximizing the sample log-likelihood given by:

*\log L(\boldsymbol{θ}; H_{t_n}) = ∑_{i=1}^{n}
\log λ(x_i,y_i,t_i; \boldsymbol{θ})-
\int_{T_0}^{T_{max}} \int \int_{Ω_{(x,y)}}\,
λ(x,y,t;\boldsymbol{θ})\,d x \, d y \,d t
*

returns an object of class `etasclass`

.

The main items of the output are:

`this.call` |
reports the exact call of the function |

`params.ind` |
indicates which parameters have been estimated (see details) |

`params` |
ML estimates of the ETAS parameters. |

`sqm` |
Estimates of standard errors of the ML estimates of the ETAS parameters ( |

`AIC.iter` |
AIC values at each iteration. |

`hdef` |
final bandwidth used for the kernel estimation of background spatial intensity (however estimated, with |

`rho.weights` |
Estimated probability for each event to be a background event ( |

`time.res` |
rescaled time residuals (for time processes only). |

`params.iter` |
A matrix with estimates values at each iteration. |

`sqm.iter` |
A matrix with the estimates of the standard errors at each iteration. |

`rho.weights.iter` |
A matrix with the values of |

`l` |
A vector with estimated intensities, corresponding to observed points |

`summary`

,
`print`

and
`plot`

methods are defined for an object of class `etasclass`

to obtain main output.

A `profile`

method (`profile.etasclass`

) is also defined to make approximate inference on a single parameter

In this first version the x-y space region, where the point process is defined, is a rectangle embedding the catalog values.

The optimization algorithm depends on the choice of initial values. Some default guess choice is performed
inside the function for parameters without input starting values. If convergence problem are experienced, a useful strategy can be to start with an high magnitude threshold value *m0* (that is, with a smaller catalog with bigger earthquakes), and then using this first output as starting guess for a running with a lower magnitude threshold value *m0*.
In this trial executions avoid declustering (`declustering=FALSE`

) or at least use a small value of `ndeclust`

; small values of `iterlim`

and `ntheta`

can speed first executions.

Quicker executions are obtained using smaller values of `iterlim`

and `ntheta`

in the input.

Also a first execution with `is.backconstant = TRUE`

, to fit a first approximation model with constant background, can be useful.

Some other useful information can be obtained estimating a pure time process, that can give a good guess at least for some parameters, like *mu, k0, a, c, p*.

Input times are expected in days, and so final intensities are expected number of events per day. If input values are in seconds, then set `sectoday=TRUE`

Marcello Chiodi, Giada Adelfio

Adelfio G, Chiodi M (2015). Alternated Estimation in Semi-Parametric Space-Time
Branching-Type Point Processes with Application to Seismic Catalogs. *Stochastic Environmental Research and Risk Assessment*, **29**(2), 443-450. doi:10.1007/s00477-014-0873-8.

Adelfio G, Chiodi M (2015). FLP Estimation of Semi-Parametric Models for Space-Time
Point Processes and Diagnostic Tools. *Spatial Statistics*, **14**(B), 119-132. doi:10.1016/j.spasta.2015.06.004.

Adelfio, G. and Chiodi, M. (2013) Mixed estimation technique in semi-parametric space-time point processes for earthquake description.
*Proceedings of the 28th International Workshop on Statistical Modelling 8-13 July, 2013, Palermo* (Muggeo V.M.R., Capursi V., Boscaino G., Lovison G., editors). Vol. **1** pp.65-70.

Chiodi, M. and Adelfio, G., (2011) Forward Likelihood-based predictive approach for space-time processes. *Environmetrics*, vol. **22** (6), pp. 749-757. DOI:10.1002/env.1121.

Chiodi, M. and Adelfio, G., (2017) Mixed Non-Parametric and Parametric Estimation Techniques in R Package etasFLP for Earthquakes' Description. *Journal of Statistical Software*, vol. **76** (3), pp. 1-28.
DOI: 10.18637/jss.v076.i03.

Zhuang, J., Ogata, Y. and Vere-Jones, D.
Stochastic declustering of space-time earthquake occurrences.
*Journal of the American Statistical Association*,
**97**, 369–379 (2002). DOI:10.1198/016214502760046925.

`eqcat`

, `plot.etasclass`

, `summary.etasclass`

, `profile.etasclass`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | ```
## Not run:
data("italycatalog")
# load a sample catalog of the italian seismicity
etas.flp=etasclass(italycatalog,
magn.threshold = 3.0, magn.threshold.back = 3.5,
k0 = 0.005, c = 0.005, p = 1.01, a = 1.05, gamma = 0.6, q = 1.52, d = 1.1,
params.ind = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
declustering = TRUE, thinning = FALSE, flp = TRUE, ndeclust = 15,
onlytime = FALSE, is.backconstant = FALSE,
description = "etas flp",sectoday = TRUE, usenlm = TRUE, epsmax = 0.001)
# execution of etasclass for events with minimum magnitude of 3.0.
# The events with magnitude at least 3.5 are used to build a first approximation
# for the background intensity function
# (magn.threshold.back=3.5)
# print method for the etasclass object
etas.flp
Call:
etasclass(cat.orig = italycatalog, magn.threshold = 3, magn.threshold.back = 3.5,
k0 = 0.005, c = 0.005, p = 1.01, a = 1.05, gamma = 0.6, d = 1.1,
q = 1.52, params.ind = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE), declustering = TRUE, thinning = FALSE, flp = TRUE,
ndeclust = 15, onlytime = FALSE, is.backconstant = FALSE,
description = "etas flp", sectoday = TRUE, usenlm = TRUE,
epsmax = 0.001)
etas flp
Number of observations 2158
ETAS Parameters:
mu k0 c p a gamma d q
0.355850 0.008373 0.009404 1.121630 1.509371 0.857945 1.915139 1.836391
# plot results with maps of intensities and diagnostic tools
plot(etas.flp)
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.