SB-ETAS: using simulation based inference for scalable, likelihood-free inference for the ETAS model of earthquake occurrences
Samuel Stockman Daniel J. Lawson Maximilian J. Werner


SB-ETAS: using simulation based inference for scalable, likelihood-free inference for the ETAS model of earthquake occurrences]SB-ETAS: using simulation based inference for scalable, likelihood-free inference for the ETAS model of earthquake occurrences 










[1]Samuel Stockmansam.stockman@bristol.ac.uk

1]Daniel J. Lawsondan.lawson@bristol.ac.uk


2]Maximilian J. Wernermax.werner@bristol.ac.uk


*[1]School of Mathematics, University of Bristol, Woodland Rd, Bristol, BS8 1UG, United Kingdom

[2]School of Earth Sciences
, University of Bristol, Wills Memorial Building, Bristol, BS8 1RL, United Kingdom






Performing Bayesian inference for the Epidemic-Type Aftershock Sequence (ETAS) model of earthquakes typically requires MCMC sampling using the likelihood function or estimating the latent branching structure. These tasks have computational complexity O(n^2) with the number of earthquakes and therefore do not scale well with new enhanced catalogs, which can now contain an order of 10^6 events. On the other hand, simulation from the ETAS model can be done more quickly O(n log n ). We present SB-ETAS: simulation-based inference for the ETAS model. This is an approximate Bayesian method which uses Sequential Neural Posterior Estimation (SNPE), a machine learning based algorithm for learning posterior distributions from simulations. SB-ETAS can successfully approximate ETAS posterior distributions on shorter catalogues where it is computationally feasible to compare with MCMC sampling. Furthermore, the scaling of SB-ETAS makes it feasible to fit to very large earthquake catalogs, such as one for Southern California dating back to 1932. SB-ETAS can find Bayesian estimates of ETAS parameters for this catalog in less than 10 hours on a standard laptop, which would have taken over 2 weeks using MCMC. Looking beyond the standard ETAS model, this simulation based framework would allow earthquake modellers to define and infer parameters for much more complex models that have intractable likelihood functions.

















[
    [
    April 29, 2024
==================





§ INTRODUCTION

The Epidemic Type Aftershock Sequence (ETAS) model <cit.> has been the most dominant way of modelling seismicity in both retrospective and fully prospective forecasting experiments <cit.> as well as in operational earthquake forecasting <cit.>.
The model characterises the successive triggering of earthquakes, making it effective for forecasting aftershock sequences. The parameters of the model are also used by seismologists to characterise physical phenomena of different tectonic regions such as structural heterogeneity, stress, and temperature <cit.>.

In recent years, an accelerated growth in the number of seismic sensors and machine learning algorithms for detecting the arrival times of earthquake phases <cit.>, has meant that earthquake catalogs have grown to a size that has overtaken what is computationally feasible to fit an ETAS model to. In California, a deployment of a dense network of seismic sensors over the last century combined with an active tectonic regime has resulted in a comprehensive dataset of earthquakes in the region <cit.>. Furthermore, in more specific areas of California, through machine learning based seismic phase picking and template matching, enhanced earthquake catalogs have been created which contain many small previously undetected earthquakes <cit.>. It is fair to assume that these datasets will only continue to grow in the future as past continuous data is reprocessed by these new detection techniques as well as being applied to future earthquakes. Determining whether increased data size results in improved earthquake forecasts is an important question for the seismological community. This work provides some methodology for answering that question by providing a simulation based estimation procedure for the ETAS model, with a scalability that make it feasible to infer ETAS parameters for such large datasets.

 

Most commonly, point estimates of ETAS parameters are found through maximum likelihood estimation (MLE). From these point estimates, forecasts can be issued by simulating multiple catalogs over the forecasting horizon. Forecast uncertainty is quantified by the distribution of simulations from MLE parameter values, however, this approach fails to quantify uncertainty contained in estimating the parameters themselves. Parameter uncertainty for MLE can be estimated using the Hessian of the likelihood <cit.>, which requires a very large sample size to be effective, or through multiple runs of the MLE procedure with different initial conditions <cit.>. 

Full characterisation of the parameter uncertainty is achieved with Bayesian inference, a procedure which returns the entire probability distribution over parameters conditioned on the observed data and updated from prior knowledge. This distribution over parameters, known as the posterior, does not have a closed form expression for the temporal and spatio-temporal ETAS model and so several approaches have used Markov Chain Monte Carlo (MCMC) to obtain samples from this distribution <cit.>. These approaches evaluate the likelihood of the ETAS model during the procedure, which is a an operation with quadratic complexity 𝒪(n^2), and therefore are only suitable for catalogs of around 10,000 events.

An alternative to MCMC is based on the Integrated Nested Laplace Approximation (INLA) method <cit.>, which approximates the marginal posterior distributions using latent Gaussian models. The implementation of INLA for the ETAS model, named , by <cit.> demonstrated a factor 10 speed-up over the MCMC method for a catalog of 3,500 events but they do not provide results on larger catalogs. Section <ref> provides a more in depth overview of the MCMC and INLA methods.

We propose SB-ETAS: an alternative method for Bayesian inference for the ETAS model. SB-ETAS is based on Simulation Based Inference (SBI), a family of approximate procedures which infer posterior distributions for parameters using simulations in place of the likelihood. Simulation from the ETAS model is 𝒪(n log n), making it a procedure which scales better than the evaluation of the likelihood, 𝒪(n^2). From a suite of SBI methods developed recently using neural density estimators <cit.>, we choose Sequential Neural Posterior Estimation (SNPE) since it does not involve any MCMC. SNPE trains a neural density estimator to approximate the posterior distribution from pairs of simulations and parameters. Section <ref> provides an overview of SNPE and other methods of SBI.

SBI is typically used in the case where the likelihood function is intractable, for example when the likelihood involves integrating over many unobserved latent variables. However, in our use case the likelihood function is tractable, presenting an alternative potential for SBI for problems where the simulator scales better than the likelihood. Looking beyond the example presented in this study, the likelihood-free inference procedure that we present, could be applied to other branching models for earthquakes. Such branching models need only be defined as a forward simulator and wouldn't be constrained to have a likelihood function. For example, a large problem in earthquake modeling is the time varying detection of earthquakes, particularly following very large magnitude events <cit.>. These undetected events bias models such as ETAS that do not account for them <cit.>. A branching model which deletes simulated events with a time varying detection rate, does not result in a tractable likelihood but could still be fit to data using our procedure.

Other branching process models for earthquakes that incorporate physical quantities <cit.>, on the other hand, do not have tractable likelihoods and currently use parameter values that are fit independent of the physical part of the model. This work provides a framework for estimating the parameters for such models, with the hope that better estimates would lead to improved forecasts.

The remainder of this paper is structured as follows: In section <ref> we give an overview of the ETAS model along with existing procedures for performing Bayesian inference; section <ref> gives an overview of SBI, following which we describe the details of SB-ETAS in section <ref>. We present empirical results based on synthetic earthquake catalogs in section <ref> and observational earthquake data from Southern California in section <ref>, before finishing with a discussion in section <ref>.


§ THE ETAS MODEL

The temporal Epidemic Type Aftershock Sequence (ETAS) model <cit.> is a marked Hawkes process <cit.> that describes the random times of earthquakes t_i along with their magnitudes m_i in the sequence x = {(t_1,m_1),(t_2,m_2),…,(t_n,m_n)  }∈ [0,T]^n ×ℳ^n ⊂ℝ_+^n ×ℳ^n. The natural ordering of events leads to the definition of a filtration ℋ_t known as the history of the process, which contains information on all events up to time t. Marked Hawkes processes are usually specified by their conditional intensity function <cit.>,

    λ(t,m|ℋ_t) =
       lim_Δ t, Δ m → 0𝔼[N([t,t+Δ t) × B(m,Δ m)|ℋ_t]/Δ t |B(m,Δ m)|.

Where the counting measure N(A) counts the events in the set A ⊂ℝ_+ ×ℳ and |B(m,Δ m)| is the Lebesgue measure of the ball B(m,Δ m) with radius Δ m.

The ETAS model typically has the form,

    λ(t,m|ℋ_t)  = (μ + ∑_i:t_i<t g(t-t_i,m_i) )f_GR(m),

where μ is a constant background rate of events and g(t,m) is a non-negative excitation kernel which describes how past events contribute to the likelihood of future events. For the ETAS model this triggering kernel factorises the contribution from the magnitude and the time,

    g(t,m) =  k(m; K, α) h(t;c,p)
       k(m; K, α)  = Ke^α(m-M_0) : m ≥ M_0
       h(t;c,p) =  c^p-1(p-1)(t+c)^-p: t ≥ 0

where the k(m; K, α) is known as the Utsu law of productivity <cit.> and h(t;c,p) is a power law known as the Omori decay <cit.>.
The magnitudes are said to be “unpredictable" since they do not depend on previous events and are distributed according to the Gutenberg-Richter law for magnitudes <cit.> with probability density f_GR(m) = β e^β(m-M_0) on the support m : m ≥ M_0.


 §.§ Branching Process Formulation

An alternative way of formulating the ETAS model is as a Poisson cluster process <cit.>. In this formulation a set of immigrants I are realisations of a Poisson process with rate μ. Each immigrant t_i ∈ I has a magnitude m_i with probability density f_GR and generates offspring S_i from an independent non-homogeneous Poisson process, with rate g(t-t_i,m_i). Each offspring t_j ∈ S_i also has magnitude m_j with probability density f_GR and generate offspring S_j of their own. This process is repeated over generations until a generation with no offspring in time interval [0,T] is produced. Although it is not observed in the data, this process is accompanied by latent branching variables B = {B_1, …, B_n } which define the branching structure of the process,

    B_i = 
            0 if  t_i ∈ I (i.e.  i  is a background event)
    
            j if  t_i ∈ S_j  (i.e.  i  is an offspring of   j)

Both the intensity function formulation as well as the branching formulation define different methods for simulating as well as inferring parameters for the Hawkes process. We now give a brief overview of some of these methods, however, we direct the reader to <cit.> for a more detailed review.


 §.§ Simulation

A simulation algorithm based on the conditional intensity function was proposed by <cit.>. This algorithm requires generating events sequentially using a thinning procedure. Simulating forward from an event t_i, the time to the next event τ is proposed from a Poisson process with rate λ(t_i|ℋ_t_i). The proposed event t_i+1 = t_i + τ is then rejected with probability 1- λ(t_i+τ|ℋ_t_i)/λ(t_i|ℋ_t_i). This procedure is then repeated from the newest simulated event until a proposed event falls outside a predetermined time window [0,T]. 

This procedure requires evaluating the intensity function λ(t|ℋ(t)) at least once for each one of n events that are simulated. Evaluating the intensity function requires a summation over all events before time t, thus giving this simulation procedure time complexity 𝒪(n^2).
 

Another algorithm which instead simulates using the branching process formulation of the ETAS model was proposed by <cit.>. This algorithm simulates events over generations G^(i), i=0,…, until no more events fall within the interval [0,T].


This procedure has time complexity 𝒪(n) for steps 1-5, since there is only a single pass over all events. Additional time constraint is added in step 6. where the whole set of events are sorted chronologically, which is at best 𝒪(nlog n).


 §.§ Bayesian Inference

Given we observe the sequence x_obs = {(t_1,m_1),(t_2,m_2),…,(t_n,m_n)  } in the interval [0,T], we are interested in the posterior probability p(θ|x_obs) for the parameters θ = (μ, K, α,c,p) of the ETAS model defined in (<ref>)-(<ref>), updated from some prior probability p(θ). The posterior distribution, expressed in Bayes rule,

    p(θ|x_obs) ∝ p(x_obs|θ)p(θ),

is known up to a constant of proportionality through the product of the prior p(θ) and the likelihood p(x_obs|θ), where,

    log p(x_obs|θ)    = ∑_i=1^n log[μ + ∑_j=1^i-1 h(t;c,p) k(m; K, α) ] 
       - μ T - ∑_i=1^nk(m; K, α)H(T-t_i; c, p).

Here, H(t; c,p) = ∫_0^t h(s;c,p)ds, denotes the integral of the Omori decay kernel.

<cit.> draw samples from the posterior p(θ|x_obs) through independent random walk Markov Chain Monte Carlo (MCMC) with Metropolis-Hastings rejection of proposed samples. This approach, however, can suffer from slow convergence due parameters of the ETAS model having high correlation <cit.>.

<cit.> developed an MCMC sampling scheme, , which conditions on the latent branching variables B = {B_1, …, B_n }. The scheme iteratively estimates the branching structure,

    ℙ(B_i^(k) = j| x_obs , θ^(k-1)) =
       μ^(k-1)/μ^(k-1) + ∑_j=1^i-1k(m_j)h(t_i-t_j) : j=0
    k(m_j)h(t_i-t_j)/μ^(k-1) + ∑_j=1^i-1k(m_j)h(t_i-t_j) : j = 1,2,…,i-1
 
and then samples parameters, θ^(k) = (μ^(k),K^(k), α^(k),c^(k),p^(k)), from the conditional likelihood,

    log p(x_obs|θ, B^(k)) = 
       |S_0^(k)|logμ - μ T
        + ∑_j=1^n (-k(m_j;K,α)H(T-t_j;c,p) 
       +  |S_j^(k)| log k(m_j;K,α)  + ∑_t_i ∈ S_j^(k)log h(t_i-t_j; c,p) )

Where |S_j| denotes the number of events that were triggered by the event at t_j.

By conditioning on the branching structure, the dependence between parameters (K,α) and (c,p) is reduced, decreasing the time it takes for the sampling scheme to converge. We can see from equation (<ref>) that estimating the branching structure from the data is a procedure that is 𝒪(n^2). Since for every event i = 1,…,n, to estimate its parent we must sum over j = 1,…,i-1. For truncated version of the time kernel h(t), this operation can be streamlined to 𝒪(n), however, this is not the case in general.
 

More recently <cit.> have constructed an approximate method of Bayesian inference for the ETAS model based on an Integrated Nested Laplace Approximation (INLA) implemented in the R-package  as well as a linear approximation of the likelihood. This approach expresses the log-likelihood as 3 terms,

    log p(x_o | θ) =
        -Λ_0(x_obs,θ) - ∑_i=1^n Λ_i(x_obs,θ) + ∑_i=1^n logλ(t_i|ℋ_t_i),

where, 

    Λ_0(x_obs,θ)    = ∫_0^Tμ dt, 
    Λ_i(x_obs,θ)    = ∑_h=1^B_i∫_b_h,i g(t-t_i,m_i)dt
       = ∑_h=1^B_iΛ_i(x_obs,θ,b_h,i).

Where b_1,i,…,b_B_i,i are chosen to partition the interval [0,T].
The log-likelihood is then linearly approximated with a first order Taylor expansion with respect to the posterior mode θ^*,

    log p(x_obs|θ ; θ^*) = -Λ_0(x_obs,θ,θ^*) 
       - ∑_i=1^n∑_h=1^B_iΛ_i(x_obs,θ,b_h,i; θ^*) + ∑_i=1^n logλ(x_obs,θ;θ^*) 
        = -exp{logΛ_0(x_obs,θ,θ^*)}
       - ∑_i=1^n∑_h=1^B_iexp{logΛ_i(x_obs,θ,b_h,i; θ^*)}
       + ∑_i=1^n logλ(x_obs,θ;θ^*),

where the notation, Λ, denotes the approximation of Λ and logΛ(   ;θ^*) denotes the first order Taylor expansion of logΛ about the point θ^*.

The posterior mode θ^* is found through a Quasi-Newton optimisation method and the final posterior densities are found using INLA, which approximates the marginal posteriors p(θ_i|x_obs) using a latent Gaussian model.

This approach speeds up computation of the posterior densities, since it only requires evaluation of the likelihood function during the search for the posterior mode. However, the approximation of the likelihood requires partitioning the space into a number of bins, which the authors recommend choosing as greater than 3 per observation. This results in the approximate likelihood having complexity 𝒪(n^2).


§ SIMULATION BASED INFERENCE

A family of Bayesian Inference methods have evolved from application settings in science, economics or engineering where stochastic models are used to describe complex phenomena. In this setting, the model may simulate data from a given set of input parameters, however, the likelihood of observing data given parameters is intractable. The task in this setting is to approximate the posterior p(θ|x_obs) ∝ p(x_obs|θ)p(θ), with the restriction that we cannot evaluate p(x|θ) but we have access to the likelihood implicitly through samples x_r ∼ p(x|θ_r) from a simulator, for r = 1,…,R and where θ_r ∼ p(θ). This approach is commonly referred to as Simulation Based Inference (SBI) or likelihood-free inference.

Until recently, the most predominant approach for SBI was Approximate Bayesian Computation <cit.>. In its simplest form, parameters are chosen from the prior θ_r ∼ p(θ), r=1,…,R, the simulator then generates samples x_r ∼ p(x|θ_r), r=1,…,R, and each sample is kept if it is within some tolerance ϵ of the observed data, i.e. d(x_r,x_obs) < ϵ. 

This approach, although exact when ϵ→ 0, is inefficient with the use of simulations. Sufficiently small ϵ requires simulating an impractical number of times, and this issue scales poorly with the dimension of x. In light of this an MCMC approach to ABC makes proposals for new simulator parameters θ_r ∼ q( |θ_r-1) using a Metropolis-Hastings kernel <cit.>. This leads to a far higher acceptance of proposed simulator parameters but still scales poorly with the dimension of x.

In order to cope with high dimensional simulator outputs x∈ℝ^n, summary statistics S(x)∈ℝ^d are chosen to reduce the dimension of the sample whilst still retaining as much information as possible. These are often chosen from domain knowledge or can be learnt as part of the inference procedure <cit.>. Summary statistics S(x) are then used in place of x in any of the described methods for SBI.


 §.§ Neural Density Estimation

Recently, SBI has seen more rapid development as a result of neural network based density estimators, which seek to approximate the density p(x) given samples of points x ∼ p(x). A popular method for neural density estimation is normalising flows <cit.>, in which a neural network parameterizes an invertible transformation x = g_ϕ(u), of a variable u from a simple base distribution p(u) into the target distribution of interest. In practice, the transformation is typically composed of a stack of invertible transformations, which allows it to learn the complex target density. The parameters of the transformation are trained through maximising the likelihood of observing p_g(x), which is given by the change of variables formula. Since x is expressed as a transformation of a simple distribution u ∼ p(u), samples from the learnt distribution p_g(x) can be generated by sampling from p(u) and passing the samples through the transformation. Neural density estimators may also be generalised to learn conditional densities p(x|y) by conditioning the transformation g_ϕ on the variable y <cit.>.

In the task of SBI, a neural density estimator can be trained on pairs of samples θ_r ∼ p(θ), x_r ∼ p(x|θ_r) to approximate either the likelihood p(x_obs|θ) or the posterior density p(θ|x_obs), from which posterior samples can be obtained. If the posterior density is estimated, in a procedure known as Neural Posterior Estimation (NPE) <cit.>, then samples can be drawn from the normalising flow. If the likelihood is estimated, known as Neural Likelihood estimation (NLE) <cit.>, then the approximate likelihood can be used in place of the true likelihood in a MCMC sampling algorithm to obtain posterior samples. These density estimation techniques generally outperform ABC-based methods since they are able to efficiently interpolate between different simulations (Figure <ref>, <cit.>). Other neural network methods exist for SBI such as ratio estimation <cit.> or score matching <cit.>, however, we direct the reader to <cit.> for a more comprehensive review of modern SBI.





§ SB-ETAS


We now present SB-ETAS, our scalable method for Bayesian inference for the ETAS model. The method avoids having to evaluate the computationally expensive likelihood function and instead leverages fast simulation from the ETAS branching process. The inference method uses Sequential Neural Posterior Estimation (SNPE) <cit.>, a modified version of NPE which performs inference over rounds. Each round, an estimate of the posterior proposes new samples for the simulator, a neural density estimator is trained on those samples and the estimated posterior is updated (Figure <ref>, Algorithm <ref>). SNPE was chosen over other methods of Neural-SBI, as it avoids the need to perform MCMC sampling, a slow procedure. Instead sampling from the posterior is fast since the approximate posterior is a normalising flow.




 §.§ Summary Statistics

What differentiates Hawkes process models from other simulator models used in Neural-SBI is that the output of the simulator x = (t_1,m_1),…,(t_n,m_n) itself has random dimension. For a specified time interval over which to simulate earthquakes [0,T], one particular parameter θ_1 will generate different numbers of events if simulated repeatedly. This is problematic for neural density estimators since even though they are successful over high dimensional data, they require a fixed dimensional input. For this reason, we fix the dimension of simulator output through calculating summary statistics of the data S(x). 

Works to perform ABC on the univariate Hawkes process with exponential decay kernel have found summary statistics that perform well in that setting. <cit.> use the histogram of inter-event times as summary statistics as well as the number of events. <cit.> extend these summary statistics by adding Ripley's K statistic <cit.>, which is a popular choice of summary statistic in spatial point processes <cit.>. Figure <ref> shows the performance of the ABC-MCMC method developed by <cit.>, using the aforementioned summary statistics. Using SNPE on the same summary statistics yields a more confident estimation of the “true” posterior which is found through MCMC sampling using the likelihood and requires far fewer simulations (10,000 versus 300,000).

The ETAS model is more complex than a univariate Hawkes process since it is both marked (i.e. it contains earthquake magnitudes) and contains a power law decay kernel which decays much more slowly than exponential making it harder to estimate <cit.>. For SB-ETAS we borrow similar summary statistics to <cit.>, namely S_1(x) = log( # events), S_2,…,S_4(x) =20th, 50th and 90th quantiles of the inter-event time histogram. Similar to <cit.>, we use another statistic S_5(x) = the ratio of the mean and median of the inter-event time histogram.


  §.§.§ Ripley's K Statistic

For the remaining summary statistics, we develop upon the introduction of Ripley's K statistic by <cit.>. For a univariate point process x = (t_1,…,t_n), Ripley's K statistic is <cit.>,

    K(x,w) = 
       1/λ𝔼(# of events within w of a random event ).

Here, λ is the unconditional rate of events in the time window [0,T].  An estimator for the K-statistic is derived in <cit.>,

    K̂(x,w) = T/n^2∑_i=1^n ∑_j≠ i𝕀(0< t_j-t_i ≤ w).

Despite containing a double-sum, computation of this estimator has complexity 𝒪(n) since {t_i}_i=1^n is an ordered sequence, i.e. (t_3 - t_1 < w) ⇒ (t_3 - t_2 < w). Calculation of Ripley K-statistic therefore satisfies the complexity requirement of our procedure if the number of windows w for which we evaluate K̂(x,w) is less than log n. In fact, our results suggest that less than 20 are required.

The use of Ripley's K-statistic for non-marked Hawkes data is motivated by <cit.>, where the authors show that second-order properties fully characterise a Hawkes process and can be used to estimate a non-parametric triggering kernel. <cit.> go on to give a recommendation for a binning strategy to estimate slow decay kernels such as a power law, using a combination of linear and log-scaling. It therefore seems reasonable to define S_6(x)… S_23(x)= K̂(x,w), where w scales logarithmically between [0,1] and linearly above 1.

We modify Ripley's K statistic to account for the particular interaction between marks and points in the ETAS model. Namely, the magnitude of an earthquake directly affects the clustering that occurs following it, expressed in the productivity relationship (<ref>). In light of this, a quantity of interest is a magnitude thresholded Ripley K-statistic,

    K_T(x,w,M_T) = 
       1/λ_T𝔼(# events within w of an event m_i≥ M_T).

Where λ_T is the unconditional rate of events above M_T. One can see that K_T(x,w,M_0) = K(x,w). We estimate K_T with,

    K̂_T(x,w,M_T) = T/ω^2∑_i:m_i≥ M_T∑_j≠ i𝕀(0< t_i-t_j ≤ w),

where ω is the number of events above magnitude threshold M_T. For general M_T, we lose the 𝒪(n) complexity that the previous statistic has, instead it is 𝒪(ω n). However, if the threshold is chosen to be large enough, evaluation of this estimator is fast. In our experiments, M_T is chosen to be (4.5,5,5.5,6), with w = (0.2,0.5,1,3). This defines the remaining statistics S_24(x),…,S_39(x).



§ EXPERIMENTS AND RESULTS

To evaluate the performance of SB-ETAS, we conduct inference experiments on a series of synthetic ETAS catalogs. On each simulated catalog we seek to obtain 5000 samples from the posterior distribution of ETAS parameters. The latent variable MCMC inference procedure, , will be used as a reference model in our experiments since it uses the ETAS likelihood without making any approximations. We compare samples from this exact method with samples from approximate methods,  and SB-ETAS.

Multiple catalogs are simulated from a fixed set of ETAS parameters, (μ,k,α,c,p) = (0.2,0.2,1.5,0.5,2) with magnitude of completeness M_0 = 3 and Gutenberg-Richter distribution parameter β = 2.4. Each new catalog is simulated in a time window [0,T], where T ∈ (10,20,30,40,50,60,70,80,90,100,250,500,1000) ×  10^3. This allows us to vary the number of events in each simulated catalog. 


Figure <ref> shows the runtime of each inference method as a function of the number of events in each catalog. Each method was run on a high-performance computing node with 8 2.4 GHz Intel E5-2680 v4 (Broadwell) CPUs, which is equivalent to what is commonly available on a standard laptop. On the catalogs with up to 100,000 events,  is the fastest inference method, around ten times quicker on a catalog of 20,000 events. However, the superior scaling of SB-ETAS allows it to be run on the catalog of ∼ 500,000 events, which was unfeasible for  given the same computational resources i.e. it exceeded a 336 hour time limit. The gradient of 2 for both  and  in log-log space confirm the 𝒪(n^2) time complexity of both inference methods. SB-ETAS, on the other hand, has gradient 2/3 which suggests that the theoretical 𝒪(nlog n) time complexity is a conservative upper-bound.

The prior distribution of parameters for each implementation are not identical since each has its own requirements. Priors are chosen to replicate the fixed implementation in the  package,


    μ   ∼Gamma(0.1,0.1)
    
        K,α,c    ∼Unif(0,10)
    
        p    ∼Unif(1,10).

The implementation of  uses a transformation K_b = K(p-1)/c, with prior K_b ∼LogN(-1,2.03) chosen by matching 1% and 99% quantiles with the  prior for K. SB-ETAS uses a μ∼Unif(0.05,0.3) prior in place of the gamma prior as well as enforcing a non-critical parameter region Kβ < β-α <cit.>. Both the uniform prior and the restriction on K and α stop unnecessarily long or infinite simulations.

Once samples are obtained from SB-ETAS and , we measure their (dis)similarity with samples from the exact method  using the Maximum Mean Discrepancy (MMD) <cit.> and the Classifier Two-Sample Test (C2ST) <cit.>. Figure <ref> shows the value of these performance metrics for samples from each of 15 rounds of simulations in SB-ETAS compared with the performance of . Since SB-ETAS involves random sampling in the procedure, we repeat it for 10 different initial seeds and plot a 95% confidence intervals. Across the 10 synthetic catalogs, SB-ETAS performs best in terms of MMD and  performs best in terms of C2ST. Figure <ref> shows samples from the T=60,000 catalog where the performance metrics rank the inference methods differently. Samples from  are overconfident with respect to the  samples, whereas SB-ETAS samples are more conservative. This phenomenon is shared across the samples from all the simulated catalogs and we speculate that it accounts for the metrics ranking the methods differently.




A common measure for the appropriateness of a prediction's uncertainty is the coverage property <cit.>. The coverage of an approximate posterior assesses the quality of its credible regions Ĉ_x_obs which satisfy,

    γ = 𝔼_p̂(θ|x_obs)(𝕀{θ̂∈Ĉ_x_obs})

An approximate posterior has perfect coverage if its operational coverage,

    b(x_obs) = 𝔼_p(θ|x_obs)(𝕀{θ̂∈Ĉ_x_obs})

is equal to the credibility level γ. The approximation is conservative if it has operational coverage b(x_obs) > γ and is overconfident if b(x_obs) < γ <cit.>. Expectations in equations (<ref>)-(<ref>) cannot be computed exactly and so are replaced with Monte Carlo averages, resulting in empirical coverage c(x_obs). Figure <ref> shows the empirical coverage for both SB-ETAS, averaged across the 10 initial seeds, along with  on the 10 synthetic catalogs.  consistently gives overconfident approximations, as the empirical coverage lies well below the credibility level. SB-ETAS has empirical coverage that indicates conservative estimates, but that is generally closer to the credibility level.





§ SCEDC CATALOG

We now continue our evaluation of the performance of SB-ETAS on some observational data from Southern California. The Southern California Seismic Network has produced an earthquake catalog for Southern California since 1932 <cit.>. This catalog contains many famous large earthquakes such as the 1952 M_W 7.5 Kern County, 1992 M_W 7.3 Landers, 1999
M_W 7.1 Hector Mine and M_W7.1 Ridgecrest sequences. The dataset contains N=43,537 events from  -  with earthquake magnitudes ≥ M_W 2.5 and can be downloaded from the Southern California Earthquake Data Center <https://service.scedc.caltech.edu/ftp/catalogs/SCSN/>. 

This size of catalog contains too many events to find ETAS posteriors using  (i.e. it would take longer than 2 weeks). Therefore we run only SB-ETAS and  on the entire catalog and validate their performance by comparing the compensator, Λ^*(t;θ) = ∫_0^t λ^*(s;θ)ds, with the observed cumulative number of events in the catalog N(t). Λ^*(t;θ) gives the expected number of events at time t, and therefore a model and its parameters are consistent with the observed data if Λ^*(t;θ) = N(t).



We generate 5,000 samples using SB-ETAS and  and use each sample to generate a compensator curve Λ^*(t;θ). We display 95% confidence intervals of these curves in Figure <ref>, along with a curve for the Maximum Likelihood Estimate (MLE). Consistent with the synthetic experiments, we find that SB-ETAS gives a conservative estimate of the cumulative number of events across the catalog, whereas  is overconfident and does not contain the observed number of events within its very narrow confidence interval. Both  and the MLE fit to the total observed number of events in the catalog, since this value, Λ^*(T), is a dominant term in each of their loss functions (the likelihood) during estimation.


We fix the α parameter equal to the β parameter of the Gutenberg-Richter law f_GR(m), a result that is consistent with other temporal only studies of Southern California <cit.>, and which reproduces Båth’s law for aftershocks <cit.>.

We were unable to successfully fix α for  and therefore use the 5 parameter implementation of ETAS. Posterior distributions are displayed in Figures <ref> and <ref> including α = β and free α implementations of the MLE. Although the modes of the marginal distributions do not match the MLE, the SB-ETAS posteriors contain the MLE parameters within their wider confidence ranges. Since  has much narrower confidence, although the distributions are relatively close to the MLE, the confidence ranges do not contain MLE parameters.







§ DISCUSSION AND CONCLUSION

The growing size of earthquake catalogs generated through machine learning based phase picking and an increased density of seismic networks has overtaken the ability of our current forecasting models to fit to such large datasets. Although improvements have been made to reduce the computational time of performing Bayesian inference for the ETAS model, first with  followed by , both of these approaches do not improve upon the scalability of inference. The implication of this is that as catalogs continue to grow in size, these methods become less feasible to use. On experiments where we give SB-ETAS,  and  access to the same 8 CPUs, only SB-ETAS could be used to fit a catalog of 500,000 events and was the fastest method for catalogs above 100,000 events. Both  and SB-ETAS are parallelized methods and would therefore see a reduction in runtime if given access to more CPUs. This is unlike  which is not parrallelized in its current implementation. It is also worth noting that although SB-ETAS and  were given the same CPUs,  required over 4 times the amount of memory than SB-ETAS with catalogs over 100,000 events, (Figure <ref>).



Comparing samples from both approximate methods to samples from  demonstrates the quality of each inference method. Our general finding is that  provides overconfident estimates of posteriors whereas SB-ETAS gives conservative estimates. Although it might seem reasonable to judge an approximate posterior by its closeness to the exact posterior, for practical use, overconfident estimates should be penalised more than under-confident ones. Bayesian inference for the ETAS model seeks to identify a range of parameter values which are then used to give confidence over a range of forecasts. However, failure to identify regions of the parameter space that give likely parameters would result in omission of a range of likely forecasts.

SB-ETAS provides an efficient procedure for performing Bayesian Inference for a temporal ETAS model on large earthquake catalogs. However, Operational earthquake forecasting requires spatial forecasts in addition to temporal forecasts. The spatio-temporal form of the ETAS model extends the temporal model used in the study by modeling earthquake spatial interactions with an isotropic Gaussian spatial triggering kernel <cit.>. This spatio-temporal ETAS model is also defined as a branching process and would therefore retain the 𝒪(nlog n) complexity of simulation. This study has illustrated that the Ripley K-statistic is an informative summary statistic for the triggering parameters of the temporal ETAS model. It seems fair to assume that the spatio-temporal Ripley K-statistic,

    K̂(x,w_t,w_s) = 
       A T/n^2∑_i=1^n ∑_j≠ i𝕀(0< t_j-t_i ≤ w_t)𝕀(||s_j-s_i||_2 ≤ w_s).

where A is the area of the study region, would be a reasonable choice for the spatio-temporal form of SB-ETAS. This statistic loses the 𝒪(n) efficiency that the purely temporal one benefits from. Instead <cit.> have developed a distributed procedure for calculating this statistic with 𝒪(nlog n) complexity that would retain the overall time complexity that SB-ETAS has.

SB-ETAS demonstrates that inference for parameters of a clustered temporal point process can be achieved without having to evaluate an expensive likelihood function. However, generally simulation based inference is used in cases where the likelihood function is intractable. For example, consider the same ETAS branching process used in this study, but instead events are deleted with a time varying probability h(t). The induced likelihood of this process,

    p(x|θ) ∝∫ p(x,x_u|θ) ∏_t_i ∈xh(t_i)∏_t_j ∈x_u(1-h(t_j))dx_u,

is intractable since it involves integrating over the set of unobserved events x_u <cit.>. This time-varying detection of events is a huge source of bias in earthquake models and current methods to deal with it typically estimate the true earthquake rate from the apparent earthquake rate assuming no contribution from undetected events <cit.>. There are biases from ignoring such triggering <cit.> which could be avoided by using this simulation based random deletion model with intractable likelihood Eq. <ref>.

Another application of this simulation based framework is to earthquake branching models that include complex physical dependencies. One possible example would be to calibrate the UCERF3-ETAS model, a unified model for fault rupture and ETAS earthquake clustering <cit.>. This model extends the standard ETAS model by explicitly modeling fault ruptures in California and includes a variable magnitude distribution which significantly affects the triggering probabilities of large earthquakes. This model is only defined as a simulator and uses ETAS parameters found independently to the joint ETAS and fault model.  In fact, <cit.> validate the models performance through a comparison of summary statistics from the outputs of the model. This validation could be extended to comprise part of the inference procedure for model parameters using the same simulation based framework as SB-ETAS.




























































§ MEMORY


































