etas-diagnostic-framework

A Machine Learning Diagnostic Framework for Nonstationary ETAS Aftershock Modeling

Giuseppe Petrillo and Luca Dal Zilio — Earth Observatory of Singapore / Nanyang Technological University

---

Overview

This repository provides the complete numerical framework associated with the study:

G. Petrillo and L. Dal Zilio, A Machine Learning Diagnostic Framework for Nonstationary ETAS Aftershock Modeling, submitted (2025).

The framework integrates three independent analytical components:

1. Stationary ETAS inference — maximum likelihood estimation of the five standard ETAS parameters from an earthquake catalog.
2. Nonstationary ETAS inference — penalized log-spline extension in which the background rate µ(t) and triggering productivity K₀(t) are allowed to vary continuously in time, with ABIC-based smoothing selection.
3. Model-agnostic anomaly detection — unsupervised Isolation Forest analysis applied to multivariate catalog features, serving as an independent diagnostic to assess whether statistically inferred parameter transients correspond to observable changes in seismicity properties.

The framework is applied to the aftershock sequence of the 2025 Mw 7.1 Dingri earthquake in southern Tibet (Yao et al. 2025), and validated on a family of semi-synthetic catalogs with injected perturbations of controlled amplitude. While developed for the Dingri sequence, all components are general and directly transferable to any well-recorded aftershock sequence.

---

Repository Structure

etas-diagnostic-framework/
│
├── catalogs/
│   ├── catalog.txt                   # Original Dingri aftershock catalog
│   ├── decimal_days_mL.txt           # Catalog in decimal days + local magnitude
│   ├── modified_decimal_days_mL.txt  # Catalog in decimal years + magnitude (ETAS input)
│   ├── fort.10                       # Time-dependent completeness Mc(t) estimates
│   ├── catalog_rate_bump2.txt        # Semi-synthetic catalog Sim1
│   ├── catalog_rate_bump3.txt        # Semi-synthetic catalog Sim2
│   ├── catalog_rate_bump4.txt        # Semi-synthetic catalog Sim3
│   ├── catalog_rate_bump5.txt        # Semi-synthetic catalog Sim4
│   └── catalog_rate_bump6.txt        # Semi-synthetic catalog Sim5
│
├── codes/
│   ├── figure1_catalog.py            # Figure 1 — catalog overview and Mc(t)
│   ├── figure2_stat_etas.py          # Figure 2 — stationary ETAS fit and residuals
│   ├── nonstat_etas.py               # Nonstationary ETAS inference (core engine)
│   ├── figure3_plot_nonstat_etas.py  # Figure 3 — nonstationary ETAS results
│   ├── figure4-5_iso_forest.py       # Figures 4–5 — Isolation Forest on real catalog
│   ├── simu_bump.py                  # Semi-synthetic catalog generation
│   └── figure6_isolation_bump.py     # Figure 6 — Isolation Forest on synthetic catalogs
│
├── figures/
│   └── [output figures — generated by the scripts above]
│
└── README.md

---

Scientific Background

2.1 Stationary ETAS Model

The temporal ETAS model (Ogata 1988, 1998) represents seismicity as a self-exciting point process. The conditional intensity function is:

$$\lambda(t \mid \mathcal{H}_t) = \mu + \sum_{j:, t_j < t} K \exp[\alpha(m_j - m_0)],(t - t_j + c)^{-p}$$

where:

Parameter	Description
µ	Stationary background seismicity rate [events/day]
K	Triggering productivity
α	Magnitude scaling of productivity
c	Omori-Utsu time offset [days]
p	Omori-Utsu temporal decay exponent
m₀	Reference magnitude (set equal to Mₜ)

The five parameters θ = (µ, K, α, c, p) are estimated by maximizing the log-likelihood:

$$\log \mathcal{L}(\theta) = \sum_{i=1}^{N} \log \lambda(t_i \mid \mathcal{H}_{t_i}) - \int_{T_0}^{T} \lambda(t \mid \mathcal{H}_t), dt$$

The integral decomposes analytically into a background term µ(T − T₀) and a triggering term with closed-form antiderivative:

$$I_j = \begin{cases} \dfrac{(T - t_j + c)^{1-p} - c^{1-p}}{1 - p} & p \neq 1 \[6pt] \log!\left(\dfrac{T - t_j + c}{c}\right) & p = 1 \end{cases}$$

Model selection uses the Akaike Information Criterion:

$$\mathrm{AIC} = 2k - 2\log\mathcal{L}, \quad k = 5$$

Model adequacy is assessed via the time-rescaling theorem (Zhuang and Ogata 2006). Transformed times

$$\tau_i = \int_{T_0}^{t_i} \lambda(t \mid \mathcal{H}_t), dt$$

should constitute a homogeneous Poisson process with unit rate if the model is correctly specified. The empirical cumulative count N(τ) is compared against the expected relation N(τ) = τ with 95% confidence band ±1.96 √τ.

---

2.2 Nonstationary ETAS Model

The nonstationary extension allows µ and K₀ to vary continuously in time while holding the decay parameters (c, α, p) constant:

$$\lambda(t \mid \mathcal{H}_t) = \mu(t) + \sum_{j:, t_j < t} K(t_j),\exp[\alpha(m_j - M_0)],(t - t_j + c)^{-p}$$

Time dependence is parameterized via multiplicative anomaly factors relative to fixed reference values (µ_ref, K₀_ref):

$$\mu(t) = \mu_{\mathrm{ref}}, q_\mu(t), \qquad K(t) = K_{0,\mathrm{ref}}, q_K(t)$$

where positivity is enforced through an exponential link applied to cubic spline functions:

$$q_\mu(t) = \exp!\bigl(s_\mu(t)\bigr), \qquad q_K(t) = \exp!\bigl(s_K(t)\bigr)$$

The splines s_µ(t) and s_K(t) are defined on a uniform set of N_knots = 6 interior knots spanning the observation window [T₀, T], with natural boundary conditions.

2.2.1 Penalized inference

To prevent overfitting, the penalized negative log-likelihood is minimized:

$$\mathcal{J}(\theta;, w) = -\log\mathcal{L}(\theta) + w, \mathcal{R}$$

where the roughness penalty is based on the integrated squared curvature of the log-anomaly functions:

$$\mathcal{R} = \int_{T_0}^{T} \Bigl[\bigl(s_\mu''(t)\bigr)^2 + \bigl(s_K''(t)\bigr)^2\Bigr], dt$$

The smoothing weight w is selected by minimizing the Akaike Bayesian Information Criterion (ABIC) over a candidate grid {w_ℓ}:

$$\mathrm{ABIC}(w) = 2\min_\theta \mathcal{J}(\theta;, w) + 2, n_{\mathrm{hyper}}$$

where n_hyper = 1 counts the single hyperparameter w. The full nonstationary model has k = 2 × N_knots + 3 = 15 free parameters (12 spline coefficients plus c, α, p), and the corresponding AIC is:

$$\mathrm{AIC}_{\mathrm{nonstat}} = 2k + 2\min_\theta\bigl(-\log\mathcal{L}\bigr), \quad k = 15$$

Optimization uses a bound-constrained quasi-Newton algorithm (L-BFGS-B). The log-likelihood is evaluated in vectorized form via a precomputed time-difference matrix, reducing the triggering sum to a single matrix-vector product with O(N²) complexity.

---

2.3 Model-Agnostic Anomaly Detection

The earthquake catalog is partitioned into non-overlapping temporal windows of width Δt = 0.5 days. For each window k, a 13-dimensional feature vector x_k is constructed:

Feature index	Description
1	Event count N_k
2	Seismicity rate N_k / Δt [events/day]
3	Mean magnitude
4	Maximum magnitude
5	Fraction of events with M ≥ 4.5
6	Mean focal depth
7	Median focal depth
8	Standard deviation of focal depth
9	Mean longitude
10	Standard deviation of longitude
11	Mean latitude
12	Standard deviation of latitude
13	Spatial dispersion proxy σ_lat × σ_lon

Each feature is standardized using robust Median Absolute Deviation (MAD) normalization:

$$z_k^{(j)} = \frac{0.6745,\bigl(x_k^{(j)} - \mathrm{median}(x^{(j)})\bigr)}{\mathrm{MAD}(x^{(j)})}$$

yielding a robustly standardized feature matrix Z ∈ ℝ^{K × 13}. Missing values (empty bins) are replaced by the column median prior to standardization.

Anomaly detection is performed using the Isolation Forest algorithm (Liu et al. 2008) with T = 500 trees. The anomaly score for observation z_k is:

$$s(z_k) = 2^{-\mathbb{E}[h(z_k)],/,c(K)}$$

where E[h(z_k)] is the mean path length over all trees and c(K) = 2H_{K-1} − 2(K−1)/K is the normalization constant based on the harmonic number H_{K-1}. Higher scores indicate greater anomalousness. If scikit-learn is unavailable, a fallback score is used:

$$s(z_k) = \sqrt{\sum_j \bigl(z_k^{(j)}\bigr)^2}$$

---

2.4 Semi-Synthetic Catalog Generation

Semi-synthetic catalogs are generated by injecting a controlled perturbation into a real catalog within a specified time window [t₁, t₂]. The additional events are drawn from a Poisson process with rate:

$$\lambda_{\mathrm{add}} = a \cdot \lambda_{\mathrm{obs}}$$

where λ_obs is the observed event rate within the window and a (amp) is the amplification factor. Injected events form a spatially compact Gaussian cluster centered at a shifted location relative to the real centroid:

$$\mathrm{lat}_{\mathrm{new}} \sim \mathcal{N}(\bar{\mathrm{lat}} + \Delta\mathrm{lat},, \sigma_{\mathrm{cluster}}), \qquad \mathrm{lon}_{\mathrm{new}} \sim \mathcal{N}(\bar{\mathrm{lon}} + \Delta\mathrm{lon},, \sigma_{\mathrm{cluster}})$$

Magnitude and depth of synthetic events are sampled from real events within the window with small additive jitter. Outside [t₁, t₂] the catalog is identical to the original.

---

Catalog Format

catalog.txt and catalog_rate_bump*.txt

Whitespace-separated, one event per line:

yyyy  mm  dd  hh  mm  ss.ss  lat  lon  depth  Ml  seqID

Lines beginning with # are treated as comments and ignored.

decimal_days_mL.txt and modified_decimal_days_mL.txt

Two-column format used as direct input to the ETAS fitting scripts:

time_decimal_years   magnitude

fort.10

Output of the time-dependent completeness estimation (CV stabilization method). The last column contains the Mc(t) values used to produce Figure 1.

---

Installation

All scripts require Python 3.8 or later. Dependencies can be installed with:

pip install numpy scipy matplotlib scikit-learn

scikit-learn is required for Isolation Forest. If unavailable, the scripts automatically fall back to a MAD-norm-based anomaly score.

---

Usage

All scripts are designed to be run from within the codes/ directory. Input file paths are relative to this location by default and can be overridden via command-line arguments or by editing the path variables at the top of each script.

---

Figure 1 — Catalog overview and time-dependent Mc(t)

python figure1_catalog.py

Reads catalogs/decimal_days_mL.txt and catalogs/fort.10. Produces a scatter plot of the earthquake catalog with magnitude color-coded by the viridis colormap, marker size scaled quadratically with magnitude, and the smoothed time-dependent completeness magnitude Mc(t) overlaid as a red curve.

Output: figure1.svg, figure1.eps

---

Figure 2 — Stationary ETAS fit and residuals

python figure2_stat_etas.py

Reads catalogs/modified_decimal_days_mL.txt. Converts times from decimal years to days, estimates the five ETAS parameters (µ, K, c, α, p) by L-BFGS-B maximization of the log-likelihood, reports AIC, and produces a two-panel figure:

• (a) Conditional intensity λ(t) on a logarithmic scale with the earthquake catalog overlaid.
• (b) Time-rescaling residuals: empirical N(τ) versus transformed time τ, with the expected relation N(τ) = τ and 95% confidence band.

Terminal output includes the optimization status, log-likelihood, AIC, and all estimated parameters.

Output: etas_intensity_and_residuals.svg, etas_intensity_and_residuals.eps

---

Nonstationary ETAS inference (core engine)

python nonstat_etas.py

Reads catalogs/modified_decimal_days_mL.txt. Runs the full nonstationary ETAS pipeline:

1. Loads and converts the catalog to days.
2. Precomputes the time-difference matrix dt[i,j] = t_i − t_j and the lower-triangular mask.
3. Constructs uniform cubic spline knots for s_µ(t) and s_K(t).
4. Performs ABIC grid search over the smoothing weight candidates W_GRID = [0.0, 0.05, 0.1].
5. Reports the optimal w, ABIC, and estimated parameters (c, α, p).
6. Saves all outputs to output_nonstat_etas/ and produces a two-panel summary figure.

Key settings (editable at the top of the script):

Mth       = 2.5          # Magnitude of completeness
N_knots_1 = 6            # Number of spline knots for mu(t)
N_knots_2 = 6            # Number of spline knots for K0(t)
W_GRID    = [0.0, 0.05, 0.1]   # ABIC smoothing grid
mu_ref    = 0.01         # Reference background rate [events/day]
K0_ref    = 0.02         # Reference productivity [events/day]

Output files written to output_nonstat_etas/:

File	Description
catalog_used.dat	Events used for fitting (decimal year, day, magnitude)
abic_grid_results.dat	ABIC values for each candidate smoothing weight
params_optimal_vector.dat	Full optimal parameter vector
fit_summary.txt	Human-readable summary of the fit
fitted_time_series.dat	Time series of λ(t), µ(t), K₀(t) on a 600-point grid
residuals_time_transform.dat	Transformed times τᵢ and cumulative event count

---

Figure 3 — Nonstationary ETAS results

python figure3_plot_nonstat_etas.py

Reads the output files from output_nonstat_etas/ produced by nonstat_etas.py. Generates a three-panel figure:

• (a) Conditional intensity λ(t) with catalog in background.
• (b) Time-rescaling residuals N(τ) versus τ.
• (c) Time-varying background rate µ(t) and productivity K₀(t) on a logarithmic scale.

Output: output_nonstat_etas/figure_multipanel_nonstat_etas.png, .svg, .eps, .pdf

---

Figures 4–5 — Isolation Forest on the real catalog

python "figure4-5_iso_forest.py" --file ../catalogs/catalog.txt --dt 0.5 --mc 2.5 --highlight 8 10

Arguments:

Argument	Default	Description
--file	../catalogs/catalog.txt	Path to input catalog
--dt	0.5	Temporal bin width [days]
--tmin	auto	Start of analysis window [days]
--tmax	auto	End of analysis window [days]
--mc	2.5	Magnitude of completeness
--highlight	8 10	Interval to highlight [days]
--seed	0	Random seed for Isolation Forest
--no_heatmap	—	Suppress heatmap figure
--out_prefix	anomaly	Prefix for output figure files

Output: anomaly_score.png/.svg/.eps (Figure 4) and anomaly_heatmap.png/.svg/.eps (Figure 5).

---

Semi-synthetic catalog generation

python simu_bump.py \
    --input ../catalogs/catalog.txt \
    --output ../catalogs/catalog_rate_bump2.txt \
    --t1 8 --t2 10 \
    --amp 1.0 \
    --new_lat_shift 0.28 --new_lon_shift 0.28 \
    --cluster_sigma 0.10 \
    --depth_shift -4.0 \
    --seed 0

Arguments:

Argument	Default	Description
--input	../catalogs/catalog.txt	Input catalog
--output	../output_nonstat_etas/catalog_cluster_test.txt	Output catalog
--t1	8.0	Start of perturbation window [days]
--t2	10.0	End of perturbation window [days]
--amp	1.5	Rate amplification factor relative to observed rate
--dt	0.01	Poisson time step [days]
--new_lat_shift	0.08	Spatial shift of injected cluster [degrees latitude]
--new_lon_shift	0.08	Spatial shift of injected cluster [degrees longitude]
--cluster_sigma	0.003	Standard deviation of injected spatial cluster [degrees]
--depth_shift	-2.0	Systematic depth shift [km]
--mag_jitter	0.05	Random magnitude perturbation σ
--depth_jitter	0.1	Random depth perturbation σ [km]
--seed	0	Random seed

The five semi-synthetic catalogs used in the paper (Sim1–Sim5) correspond to the following perturbation parameters:

Catalog	amp	Δlat	Δlon	σ_cluster	Δd (km)
Original	0.0	0.00	0.00	—	0.0
Sim1 (catalog_rate_bump2.txt)	1.0	0.28	0.28	0.10	−4.0
Sim2 (catalog_rate_bump3.txt)	1.3	0.35	0.35	0.08	−5.0
Sim3 (catalog_rate_bump4.txt)	1.6	0.42	0.42	0.06	−6.0
Sim4 (catalog_rate_bump5.txt)	2.0	0.50	0.50	0.04	−7.0
Sim5 (catalog_rate_bump6.txt)	2.5	0.60	0.60	0.03	−8.0

---

Figure 6 — Isolation Forest sensitivity validation

python figure6_isolation_bump.py --dt 0.5 --mc 2.5 --highlight 8 10

Automatically loads the original catalog and the five semi-synthetic realizations (Sim1–Sim5) from catalogs/. Applies the identical anomaly detection pipeline to each and produces a single comparison figure showing anomaly score time series for all catalogs, with the 8–10 day window highlighted.

Output: anomaly_compare.png/.svg/.eps

---

Computational Notes

• All ETAS fitting routines have O(N²) complexity in the number of catalog events. For the Dingri catalog (~1300 events above Mc = 2.5 within the analysis window), stationary fitting completes in seconds; nonstationary fitting with the ABIC grid search takes on the order of minutes.
• The nonstationary ETAS script precomputes the full N × N time-difference matrix dt_mat once before optimization. For catalogs with N > 5000 events, memory usage may become significant and the loop-based fallback should be considered.
• The Isolation Forest is run with 500 trees (n_estimators=500) and contamination set to "auto". Results are deterministic for a fixed random seed.

---

Dependencies

Package	Version (tested)	Purpose
numpy	≥ 1.22	Array operations, catalog I/O
scipy	≥ 1.8	Optimization (L-BFGS-B), cubic splines, integration
matplotlib	≥ 3.5	All figures
scikit-learn	≥ 1.0	Isolation Forest anomaly detection

---

Citation

If you use this code or framework in your research, please cite:

@article{petrillo2025etas,
  author  = {Petrillo, Giuseppe and Dal Zilio, Luca},
  title   = {A Machine Learning Diagnostic Framework for Nonstationary {ETAS} Aftershock Modeling},
  journal = {submitted},
  year    = {2025}
}

---

Funding

G.P. and L.D.Z. were supported by the Earth Observatory of Singapore (EOS), the Singapore Ministry of Education Tier 3b project "Investigating Volcano and Earthquake Science and Technology (InVEST)" (Award No. MOE-MOET32021-0002), and the Nanyang Assistant Professorship (NAP) (Award No. 025244-00001).

---

References

• Liu, F.T., Ting, K.M., Zhou, Z.-H. (2008). Isolation forest. Proceedings of the 8th IEEE ICDM, 413–422.
• Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401), 9–27.
• Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. Annals of the Institute of Statistical Mathematics, 50(2), 379–402.
• Utsu, T., Ogata, Y., Matsu'ura, R.S. (1995). The centenary of the Omori formula for a decay law of aftershock activity. Journal of Physics of the Earth, 43(1), 1–33.
• Yao, J., Yao, D., Chen, F., Zhi, M., Sun, L., Wang, D. (2025). A preliminary catalog of early aftershocks following the 7 January 2025 Ms 6.8 Dingri, Xizang earthquake. Journal of Earth Science, 36(2), 856–860.
• Zhuang, J., Ogata, Y. (2006). Properties of the probability distribution associated with the largest event in an earthquake cluster. Physical Review E, 73, 046134.