etasFLP (version 1.0.2)

etasclass: Mixed estimation of an ETAS model

Description

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.

Usage

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),
    declustering=TRUE,thinning=FALSE,
    flp=TRUE,ndeclust=5, onlytime=FALSE,is.backconstant=FALSE,
##### end of  main input arguments. 
##### Control and secondary arguments:
    description="", cat.back=NULL, back.smooth=1.0,
    sectoday=TRUE,usenlm=TRUE, method	="BFGS", compsqm=TRUE,
    epsmax=0.0001, iterlim=100, ntheta=100)

Arguments

cat.orig

An earthquake catalog, possibly an object of class eqcat, or however a data.frame with variables of names time, lat, long, z, magn1. No missing values are

magn.threshold

Threshold magnitude (only events with a magnitude at least magn.threshold will be used). Default value = 2.5.

magn.threshold.back

Threshold magnitude used to build the catalog cat.back for the first estimation of the background seismicity. Default value = magn.threshold+2.

mu

Parameter 1 ($\mu$) of the ETAS model: background general intensity; see details. Default value = 1.

k0

Parameter 2 ($\kappa_0$) of the ETAS model: measures the strength of the aftershock activity; see details. Default value = 1.

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 ($\alpha$) of the ETAS model; efficiency of an event of given magnitude in generating aftershocks; see details. Default value = 1.2.

gamma

Parameter 6 ($\gamma$) of the ETAS model; together with a is related to the efficiency of an event of given magnitude in generating aftershocks; see details. Default value = 0.5.

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.

params.ind

vector of 8 logical values: 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,

declustering

if 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

if thinning=TRUE a background catalog is obtained sampling from the original catalog with probabilities estimated during the iterations. Default value =FALSE.

flp

if flp=TRUE then background seismicity is estimated through Forward Likelihood Predictive (see details). Otherwise the Silverman rule is used. Default value =TRUE.

ndeclust

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

onlytime

if 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 value

is.backconstant

if TRUE then background seismicity is assumed to be homogeneous in space (and declustering, flp are set to FALSE). Default value = FALSE.

description

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

cat.back

external catalog used for the estimation of the background seismicity. Default value = NULL.

back.smooth

Controls the level of smoothing for the background seismicity (meaningful only if flp=FALSE). Default value = 1.

sectoday

if TRUE, then time variable of cat.orig is converted from seconds to days. Default value = TRUE.

usenlm

if 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 = TR

method

used if usenlm=FALSE: method used by optim. Default value = "BFGS".

compsqm

if TRUE, then standard errors are computed. Default value = TRUE.

epsmax

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

iterlim

maximum number of iterations in the maximum likelihood steps (used in nlm or optim). Default value = 100.

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.

Value

• 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 (sqm[i]=0 if params.ind[i]=FALSE and for the situation where hessian is not computed or near to singularity).
• AIC.iter
  AIC values at each iteration.
• hdef
  final bandwidth used for the kernel estimation of background spatial intensity (however estimated, with flp=TRUE or flp=FALSE).
• rho.weights
  Estimated probability for each event to be a background event ($\rho$).
• 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 rho.weights at each iteration.
• 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

Details

Estimates the components of an ETAS (Epidemic type aftershock sequences) 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

$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=\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.

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{\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$$

References

Chiodi, M. and Adelfio, G., (2011) Forward Likelihood-based predictive approach for space-time processes. Environmetrics, vol. 22 (6), pp. 749-757.

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.

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).

See Also

eqcat, plot.etasclass, summary.etasclass, profile.etasclass

Examples

data("italycatalog")
# load a sample catalog of the italian seismicity

etas.flp=etasclass(italycatalog,  magn.threshold = 3.1,  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 = 10e-04)

# execution of etasclass for events with minimum magnitude of 3.1. 
# 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)

summary(etas.flp)
# summary merhod for the etasclass object
>summary(etas.flp)

Call: 
 
etasclass(cat.orig = italycatalog, magn.threshold = 3.1, 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", sectodacatalog-search.htmly = TRUE, usenlm = TRUE, 
    epsmax = 0.001)

 
etas flp 
Execution started:                  2014-02-04 13:17:29 
Elapsed time of execution (hours)   0.358737 
Number of observations             1700 
Magnitude threshold                3.1 
Number of declustering iterations   7 
Kind of declustering                weighting 
sequence of AIC values for each iteration 
40444.81 39058.33 39100.61 39101.3 39101.45 39027.12 39025.27 
 
------------------------------------------------------- 
 
ETAS Parameters: 
            Estimates       std.err.
mu           0.299141       0.010177
k0           0.008847       0.002832
c            0.012747       0.002752
p            1.149504       0.020206
a            1.640902       0.070027
gamma        0.932499       0.094271
d            2.010186       0.384154
q            1.926670       0.089593
------------------------------------------------------- 

# plot results with maps of intensities

plot(etas.flp)