etasclass is the main function of the package etasFLP.Performs the estimation of the components of the ETAS (Epidemic type aftershock sequences) 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 approach (FLP). 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.
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)eqcat, or however a data.frame with variables of names time, lat, long, z, magn1. No missing values are magn.threshold will be used). Default value = 2.5.cat.back for
the first estimation of the background seismicity. Default value = magn.threshold+2.a is related to the efficiency of an event of given magnitude
in generating aftershocks; see details. Default value = 0.5.params.ind[i] = TRUE means that the i-th parameter must be estimated. params.ind[i] = FALSE means that the i-th parameter is fixed to its input value (the order of parametrs is: mu, x,y bandwidths used in the kernel estimator of background seismicity. Default value = 1,1.TRUE the catalog is iteratively declustered to optimally estimate the background intensity (through thinning, if thinning=TRUE, or through weighting if thinning=FALSE). Default value = TRUE.thinning=TRUE a background catalog is obtained sampling from the original catalog with probabilities estimated during the iterations. Default value =FALSE.flp=TRUE then background seismicity is estimated through Forward Likelihood Predictive (see details). Otherwise the Silverman rule is used. Default value =TRUE.flp=TRUE. Indicates the range of points used for the FLP steps. See details. If missing it is set to nrow(cat)/2.TRUE then a time process is fitted to data , regardless to space location (in this case is.backconstant is set to TRUE and declustering, flp are set to FALSE). Default valueTRUE then background seismicity is assumed to be homogeneous in space (and declustering, flp are set to FALSE). Default value = FALSE.NULL.flp=FALSE). Default value = 1.TRUE, then time variable of cat.orig is converted from seconds to days. Default value = TRUE.TRUE, then long and lat variables of cat.orig are treated as geographical coordinates and converted to kilometers. Default value = TRUE.TRUE, then nlm function (gauss-newton method) is used in the maximum likelihood steps; if FALSE, then optim function is used (with method =method ). Default value = TRusenlm=FALSE: method used by optim. Default value = "BFGS".TRUE, then standard errors are computed. Default value = TRUE.nlm or optim). Default value = 100.etasclass.The main items of the output are:
sqm[i]=0 if params.ind[i]=FALSE and for the situation where hessian is not computed or near to singularity).flp=TRUE or flp=FALSE).rho.weights at each iteration.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
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
$f(x,y)$ is estimated through a weighted kernel gaussian estimator; if Weights (computed only if
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 By default The eight parameters are internally ordered in this way: for example a 7 parameters model can be fitted with if if if If $$FLP_{k_1,k_2}(\hat{\boldsymbol{\psi}})\equiv\sum_{k=k_1}^{n-1}\delta_{k,{k+1}}(\hat{\boldsymbol{\psi}}(H_{t_k});
H_{t_{k+1}})$$ with $k_1=\frac{n}{2},k_2=n-1$ and where $\delta_{k,k+1}(\hat{\boldsymbol{\psi}}(H_{t_k});
H_{t_{k+1}})$ is
the predictive information
of the first $k$ observations on the $k+1$-th observation, and is so defined: $$\delta_{k,k+1}(\hat{\boldsymbol{\psi}}(H_{t_k});
H_{t_{k+1}})\equiv \log L(\hat{\boldsymbol{\psi}}(H_{t_k}); H_{t_{k+1}} )-\log L(\hat{\boldsymbol{\psi}}(H_{t_k});H_{t_k})$$ where $H_k$ is the history of the process until time $t_k$ and
$\hat{\boldsymbol{\psi}}(H_{t_k})$ is an estimate based only on history until the $k-th$ observation.
In the ML step, the vector of parameter $\theta=(\mu, \kappa_0, c , p, \alpha, \gamma, d, q)$ is estimated maximizing the sample log-likelihood given by:
$$\log L(\boldsymbol{\theta}; H_{t_n}) = \sum_{i=1}^{n}
\log \lambda(x_i,y_i,t_i; \boldsymbol{\theta})-
\int_{T_0}^{T_{max}} \int \int_{\Omega_{(x,y)}}\,
\lambda(x,y,t;\boldsymbol{\theta})\,d x \, d y \,d t$$flp is set to TRUE
then the bandwidth is estimated through a FLP step.declustering=TRUE) are given by the estimated probabilities of being a background event; for the i-th event this is given by
$\rho_i=\frac{\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 ($\kappa_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.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. params.ind=c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), and so a full 8 parameters ETAS model will be estimated.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);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); 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);is.backconstant=TRUE a process (space-time or time) with a constant background intensity $\mu$ is fitted;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).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):
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.
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).
eqcat, plot.etasclass, summary.etasclass, profile.etasclassdata("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)Run the code above in your browser using DataLab