Package {bayesianETAS}

Draws samples from the posterior distribution of the ETAS model

Description

This function implements the latent variable MCMC scheme from (Ross 2019) which draws samples from the Bayesian posterior distribution of the Epidemic Type Aftershock Sequence (ETAS) model. Either the temporal ETAS or the spatio-temporal (with a bivariate Gaussian trigger kernel) can be sampled.

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

and s() is a spatial triggering kernel. The public spatio-temporal fitting routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component. The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i at spatial location ( x_i,, y_i) aand 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

estimateETAS(
  ts,
  ms,
  M0,
  maxTime = NULL,
  xs = NULL,
  ys = NULL,
  area = NULL,
  sims = 5000,
  burnin = 1000,
  initval = NA,
  sigma2_prior_shape = 3,
  sigma2x_prior_rate = NULL,
  sigma2y_prior_rate = NULL,
  handle_ties = c("error", "jitter"),
  c_lower = NULL
)

Arguments

Details

In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \sigma^2_x and \sigma^2_y. The spatial background is homogeneous over the supplied area and no spatial edge correction is implemented.

The inverse-gamma priors for \sigma^2_x and \sigma^2_y use a shape/rate convention. If the spatial prior rates are not supplied, the defaults are scale-aware: 2 * (0.1 * diff(range(xs)))^2 for \sigma^2_x and 2 * (0.1 * diff(range(ys)))^2 for \sigma^2_y. Degenerate coordinate ranges require explicit prior rates.

The likelihood assumes strictly increasing event times. By default duplicate times are rejected. With handle_ties = "jitter", tied event times are separated deterministically by a small amount before fitting. The Omori c parameter is constrained to be greater than c_lower; this prevents finite-resolution catalogues with exact or near ties from driving c to zero.

Value

A data frame containing the posterior samples. Each row is a single sample, and the columns correspond to (mu, K, alpha, c, p) for the temporal model and (mu, K, alpha, c, p,sigma2x,sigma2y) for the spatio-temporal model

Author(s)

Gordon J Ross

References

Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)

Examples

#temporal model
ts <- c(1, 2, 3, 5, 8)
ms <- c(3.0, 3.1, 3.3, 3.2, 3.0)
samples <- estimateETAS(ts, ms, M0=3, maxTime=10, sims=20, burnin=5,
    initval=c(0.1, 0.1, 0.5, 0.2, 1.5))
apply(samples, 2, median) #point estimate

#spatio-temporal model
xs <- c(1, 2, 3, 4, 5)
ys <- c(1, 1.5, 2, 1.5, 1)
samples <- estimateETAS(ts, ms, M0=3, maxTime=10, xs=xs, ys=ys,
    area=25, sims=20, burnin=5,
    initval=c(0.1, 0.1, 0.5, 0.2, 1.5, 1, 1))
apply(samples, 2, median) #point estimate

Estimate the parameters of the ETAS model using maximum likelihood.

Description

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

and s() is a spatial triggering kernel. The public spatio-temporal fitting routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component. The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i at spatial location ( x_i,, y_i) aand 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

maxLikelihoodETAS(
  ts,
  ms,
  M0,
  maxTime = NULL,
  xs = NULL,
  ys = NULL,
  area = NULL,
  initval = NA,
  displayOutput = FALSE,
  zeromu = FALSE,
  handle_ties = c("error", "jitter"),
  c_lower = NULL
)

Arguments

Details

In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \sigma^2_x and \sigma^2_y. The spatial background is homogeneous over the supplied area and no spatial edge correction is implemented. This function estimates the parameters of the ETAS model (both temporal and spatio-temporal) using maximum likelihood. Parameters are optimized on constrained transformed scales from multiple starting values, but the result is still a numerical maximum of a non-convex likelihood.

The likelihood assumes strictly increasing event times. By default duplicate times are rejected. With handle_ties = "jitter", tied event times are separated deterministically by a small amount before fitting. The Omori c parameter is constrained to be greater than c_lower; this prevents finite-resolution catalogues with exact or near ties from driving c to zero.

Value

A list consisting of

Author(s)

Gordon J Ross

References

Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)

Examples

#temporal model
ts <- c(1, 2, 3, 5, 8)
ms <- c(3.0, 3.1, 3.3, 3.2, 3.0)
maxLikelihoodETAS(ts, ms, M0=3, maxTime=10)

#spatio-temporal model
xs <- c(1, 2, 3, 4, 5)
ys <- c(1, 1.5, 2, 1.5, 1)
maxLikelihoodETAS(ts, ms, M0=3, maxTime=10, xs=xs, ys=ys, area=25)

Simulates synthetic data from the ETAS model

Description

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

and s() is a spatial triggering kernel. The public spatio-temporal routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component. The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time t_i at spatial location ( x_i,, y_i) aand 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

simulateETAS(
  mu,
  K,
  alpha,
  c,
  p,
  beta,
  M0,
  maxTime,
  sigma2x = NULL,
  sigma2y = NULL,
  xmin = NULL,
  xmax = NULL,
  ymin = NULL,
  ymax = NULL
)

Arguments

Details

In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \sigma^2_x and \sigma^2_y. Spatial simulation uses a bounded-window discard convention: proposed spatial offspring outside the requested window are discarded and discarded offspring do not generate descendants. This is not an edge-corrected simulation of an infinite-plane ETAS process observed through a bounded window.

This function simulates sample data from the ETAS model (either temporal or spatio-temporal) over a particular time interval [0,T], For the spatio-temporal model we also specify the range of the x axis region [xmin,xmax] and the y axis region [ymin,ymax].

Value

A list consisting of

Author(s)

Gordon J Ross

References

Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)

Examples

set.seed(1)

#temporal model
beta <- 2.4; M0 <- 3; maxTime <- 20
catalog <- simulateETAS(0.1, 0.05, 1.0, 0.5, 2, beta, M0, maxTime)

#spatio-temporal model
sigma2x <- 0.2; sigma2y <- 0.2
xmin <- 0; xmax <- 10; ymin <- 0; ymax <- 10
catalog <- simulateETAS(0.1, 0.05, 1.0, 0.5, 2, beta, M0, maxTime,
   sigma2x, sigma2y, xmin, xmax, ymin, ymax)