Title: | Bayesian Estimation of the ETAS Model for Earthquake Occurrences |
Version: | 1.0.3 |
Date: | 2017-01-16 |
Depends: | R (≥ 2.15.0) |
Author: | Gordon J. Ross |
Maintainer: | Gordon J. Ross <gordon@gordonjross.co.uk> |
Description: | The Epidemic Type Aftershock Sequence (ETAS) model is one of the best-performing methods for modeling and forecasting earthquake occurrences. This package implements Bayesian estimation routines to draw samples from the full posterior distribution of the model parameters, given an earthquake catalog. The paper on which this package is based is Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2016), available from the below URL. |
URL: | http://www.gordonjross.co.uk/bayesianetas.pdf |
License: | GPL-3 |
RoxygenNote: | 5.0.1 |
NeedsCompilation: | yes |
Packaged: | 2017-01-17 05:18:50 UTC; gordonross |
Repository: | CRAN |
Date/Publication: | 2017-01-17 08:36:26 |
Bayesian estimation of the ETAS model for earthquake occurrences
Description
Bayesian estimation of the ETAS model for earthquake occurrences
Author(s)
Gordon J Ross gordon@gordonjross.co.uk
References
Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2016), available from http://www.gordonjross.co.uk/bayesianetas.pdf
Estimate the parameters of the ETAS model using maximum likelihood.
Description
The Epidemic Type Aftershock Sequence (ETAS) model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2016). The temporal ETAS model is a point process where the probability of an earthquake occurring at time t
depends on the previous seismicity H_t
, and is defined by the conditional intensity function:
\lambda(t|H_t) = \mu + \sum_{t[i] < t} \kappa(m[i]|K,\alpha) h(t[i]|c,p)
where
\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}
and
h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}
where the summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i
and having magnitude m_i
. The quantity M_0
denotes the magnitude of completeness of the catalog, so that m_i \geq M_0
for all i. The temporal ETAS model has 5 parameters: \mu
controls the background rate of seismicity, K
and \alpha
determine the productivity (average number of aftershocks) of an earthquake with magnitude m
, and c
and p
are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law p(m_i) = \beta e^{\beta(m_i-M_0)}
.
This function estimates the parameters of the ETAS model using maximum likelihood
Usage
maxLikelihoodETAS(ts, magnitudes, M0, T, initval = NA, displayOutput = TRUE)
Arguments
ts |
Vector containing the earthquake times |
magnitudes |
Vector containing the earthquake magnitudes |
M0 |
Magnitude of completeness. |
T |
Length of the time window [0,T] the catalog was observed over. If not specified, will be taken as the time of the last earthquake. |
initval |
Initial value at which to start the estimation. A vector, with elements (mu, K, alpha, c, p) |
displayOutput |
If TRUE then prints the out the likelihood during model fitting. |
Value
A list consisting of
params |
A vector containing the estimated parameters, in the order (mu,K,alpha,c,p,beta) |
loglik |
The corresponding loglikelihood |
Author(s)
Gordon J Ross
References
Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2016), available from http://www.gordonjross.co.uk/bayesianetas.pdf
Examples
## Not run:
beta <- 2.4; M0 <- 3; T <- 500
catalog <- simulateETAS(0.2, 0.2, 1.5, 0.5, 2, beta, M0, T)
maxLikelihoodETAS(catalog$ts, catalog$magnitudes, M0, 500)
## End(Not run)
Draws samples from the posterior distribution of the ETAS model
Description
This function implements the latent variable MCMC scheme from (Ross 2016) which draws samples from the Bayesian posterior distribution of the Epidemic Type Aftershock Sequence (ETAS) model.
The ETAS model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2016). The temporal ETAS model is a point process where the probability of an earthquake occurring at time t
depends on the previous seismicity H_t
, and is defined by the conditional intensity function:
\lambda(t|H_t) = \mu + \sum_{t[i] < t} \kappa(m[i]|K,\alpha) h(t[i]|c,p)
where
\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}
and
h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}
where the summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i
and having magnitude m_i
. The quantity M_0
denotes the magnitude of completeness of the catalog, so that m_i \geq M_0
for all i. The temporal ETAS model has 5 parameters: \mu
controls the background rate of seismicity, K
and \alpha
determine the productivity (average number of aftershocks) of an earthquake with magnitude m
, and c
and p
are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law p(m_i) = \beta e^{\beta(m_i-M_0)}
.
Usage
sampleETASposterior(ts, magnitudes, M0, T = NA, initval = NA,
approx = FALSE, sims = 5000, burnin = 500)
Arguments
ts |
Vector containing the earthquake times |
magnitudes |
Vector containing the earthquake magnitudes |
M0 |
Magnitude of completeness. |
T |
Length of the time window [0,T] the catalog was observed over. If not specified, will be taken as the time of the last earthquake. |
initval |
Initial value at which to start the estimation. If specified, should be a vector, with elements (mu, K, alpha, c, p). If unspecified, the sampler will be initialized at the maximum likelihood estimate of the model parameters |
approx |
If TRUE then will approximate the true posterior using the infinite time approximation discussed in (Ross 2016) |
sims |
Number of posterior samples to draw |
burnin |
Number of burnin samples |
Value
A matrix containing the posterior samples. Each row is a single sample, and the columns correspond to (mu, K, alpha, c, p)
Author(s)
Gordon J Ross
References
Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2016), available from http://www.gordonjross.co.uk/bayesianetas.pdf
Examples
## Not run:
beta <- 2.4; M0 <- 3; T <- 500
catalog <- simulateETAS(0.2, 0.2, 1.5, 0.5, 2, beta, M0, T)
sampleETASposterior(catalog$ts, catalog$magnitudes, M0, T, sims=5000)
## End(Not run)
Simulates synthetic data from the ETAS model
Description
This function simulates sample data from the ETAS model over a particular interval [0,T].
The Epidemic Type Aftershock Sequence (ETAS) model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2016). The temporal ETAS model is a point process where the probability of an earthquake occurring at time t
depends on the previous seismicity H_t
, and is defined by the conditional intensity function:
\lambda(t|H_t) = \mu + \sum_{t[i] < t} \kappa(m[i]|K,\alpha) h(t[i]|c,p)
where
\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}
and
h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}
where the summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i
and having magnitude m_i
. The quantity M_0
denotes the magnitude of completeness of the catalog, so that m_i \geq M_0
for all i. The temporal ETAS model has 5 parameters: \mu
controls the background rate of seismicity, K
and \alpha
determine the productivity (average number of aftershocks) of an earthquake with magnitude m
, and c
and p
are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law p(m_i) = \beta e^{\beta(m_i-M_0)}
.
This function simulates sample data from the ETAS model over a particular interval [0,T].
Usage
simulateETAS(mu, K, alpha, c, p, beta, M0, T, displayOutput = TRUE)
Arguments
mu |
Parameter of the ETAS model as described above. |
K |
Parameter of the ETAS model as described above. |
alpha |
Parameter of the ETAS model as described above. |
c |
Parameter of the ETAS model as described above. |
p |
Parameter of the ETAS model as described above. |
beta |
Parameter of the Gutenberg-Richter law used to generate earthquake magnitudes. |
M0 |
Magnitude of completeness. |
T |
Length of the time window [0,T] to simulate the catalog over. |
displayOutput |
If TRUE then prints the number of earthquakes simulated so far. |
Value
A list consisting of
ts |
The simulated earthquake times |
magnitudes |
The simulated earthquake magnitudes |
branching |
The simulated branching structure, where branching[i] is the index of the earthquake that triggered earthquake i, or 0 if earthquake i is a background event |
Author(s)
Gordon J Ross
References
Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2016), available from http://www.gordonjross.co.uk/bayesianetas.pdf
Examples
## Not run:
beta <- 2.4; M0 <- 3
simulateETAS(0.2, 0.2, 1.5, 0.5, 2, beta, M0, T=500, displayOutput=FALSE)
## End(Not run)
Simulates event times from an inhomogenous Poisson process on [0,T]
Description
Simulates event times from an inhomogenous Poisson process on [0,T]
Usage
simulateNHPP(targetfn, maxintensity, T = Inf)
Arguments
targetfn |
A first order function defining the process intensity |
maxintensity |
The maximum values of targetfn |
T |
Length of the interval [0,T] on which to simulate the process |
Value
The simulated event times
Author(s)
Gordon J Ross
Examples
simulateNHPP(function(x) {sin(x)+1}, 2, 100)
simulateNHPP(function(x) {x^2}, 100, 10)