- The paper introduces ST-Score, a likelihood-free sequential method that jointly estimates when a spatio-temporal point process changes and which continuous spatial region is affected using localized weighted Hyvärinen scores.
- The method provides false-alarm, detection-delay, and localization guarantees, including an asymptotic Jaccard bound, while experiments report lower delays and higher localization accuracy than CUSUM-based baselines on Hawkes-process data.
- ST-Score operates directly on event streams and successfully identifies earthquake swarms and wildfire damage regions, although neighborhood size, optimization epochs, score-model accuracy, and Hawkes-process assumptions remain important limitations.
The paper proposes ST-Score, a likelihood-free sequential procedure for detecting changes in spatio-temporal point processes (STPPs) that jointly estimates the change time and localizes the affected spatial region in continuous space-time (2602.04798). The work addresses a gap between classical quickest change detection, which provides delay and false-alarm guarantees but treats space as fixed or ancillary, and multi-stream detection-isolation methods that assume a finite set of post-change scenarios. In the spatio-temporal setting considered here, the post-change behavior belongs to an infinite-dimensional family of alternatives, and the latent object of interest is an intensity evolving over continuous space while observations arrive as sparse event times and locations.
Events are generated by an STPP with conditional intensity λ(x∣Ht(x)) over X=[0,T)×S. At an unknown time τ, the intensity shifts from λ0 to λ1 within an unknown spatial region Ω, leaving the complement unchanged. The detection task is cast as testing this composite alternative against the null of no change. The setup covers homogeneous Poisson processes and self-exciting Hawkes processes with base-intensity changes; notably, even when only the base rate changes, self-excitation propagates pre-change history into post-change intensities, blurring both the onset and spatial extent of the change. The authors also note the framework extends to marked point processes by augmenting the domain with a mark space.
Classic CUSUM is inadequate here for two reasons: the post-change likelihood requires knowledge of Ω, which must be inferred jointly with detection, and evaluating region-specific log-likelihoods over continuous space is computationally prohibitive under unknown dynamics.
Methodology
The core statistic replaces likelihood ratios with a localized anomaly measure. For each observed event x, evaluated on its transformed inter-arrival representation (t(x)−tn,s(x)), regime-specific scores ψi(x) are computed via a weighted Hyvärinen score:
X=[0,T)×S0
where X=[0,T)×S1 approximates X=[0,T)×S2. The weighting function downweights boundary regions because the one-sided support of inter-arrival times violates standard integration-by-parts regularity conditions; the paper relies on weighted score matching results to handle truncation. The detection statistic aggregates the score difference X=[0,T)×S3 over candidate change windows and regions:
X=[0,T)×S4
with stopping when X=[0,T)×S5; the maximizers X=[0,T)×S6 at stopping provide joint estimates.
A key technical obstacle is that the true conditional density involves an integral of the intensity over all of X=[0,T)×S7, which is infeasible when X=[0,T)×S8 is unknown. The authors resolve this with a localized score approximation: each event's density is conditioned only on its X=[0,T)×S9 neighborhood τ0 of radius τ1, yielding a locally stable approximation that avoids mixing regimes across distant regions. Score models are trained by denoising score matching on offline reference datasets from each regime, using an LSTM-plus-feedforward architecture. An online variant handles the absence of post-change reference data by initializing τ2 from pre-change data and updating it via gradient descent on the denoising loss restricted to the currently estimated change region.
Optimization of the nested maximization proceeds via alternating I-steps (region update) and O-steps (change-point update). A lemma shows the region subproblem's solution set is characterized by level sets of τ3: any set containing all positive-anomaly locations and excluding negative ones is optimal; the practical estimator unions τ4-balls around positive locations minus negative ones. The O-step reduces to a finite maximization since the counting measure jumps only at event times.
Theoretical guarantees
Three guarantees are established, assuming exact score models (τ5), which the authors acknowledge requires sufficiently expressive model classes and ample training data.
False alarms: the average run length satisfies τ6, so threshold calibration scales logarithmically in the target ARL. The proof constructs exponential martingales over Voronoi-cell-generated candidate regions and applies optional sampling.
Drift structure: the expected per-event score difference equals signed weighted Fisher divergences, τ7 and τ8, confirming properness of the scoring rule. For Hawkes processes, closed-form drift expressions follow: for homogeneous Poisson processes τ9, and for exponential-kernel Hawkes processes the drift is discounted by λ00. This last result formalizes an intuitive phenomenon—stronger self-excitation (larger branching ratio λ01) dilutes detectable evidence because endogenous bursts confound exogenous change signatures—and is corroborated empirically.
Detection delay and localization: under law-of-large-numbers-type concentration assumptions on the counting process, the expected detection delay obeys
λ02
which improves as more post-change events arrive—the delay is governed by the event rate rather than wall-clock time alone. The Jaccard index of the estimated region is asymptotically lower bounded by a min of dilation/erosion ratios involving λ03:
λ04
The irreducible localization error stems from neighborhoods straddling the boundary of λ05, where mixed-regime events disrupt the drift argument. Notably, these concentration assumptions hold trivially for Poisson processes but require stationarity and light-tailedness conditions for Hawkes processes, so the delay guarantee inherits those caveats.
Numerical results
Synthetic experiments use spatio-temporal Hawkes processes on a unit box with λ06 and a square change region, comparing against CUSUM, SCUSUM, PP-CUSUM, and MinCUSUM with uniform discretization. Across all ARL levels and branching ratios λ07, ST-Score achieves lower EDD and higher Jaccard index than all baselines, operating directly on raw event streams without binning or manual channel partitioning. Qualitatively, CUSUM flags the entire domain as changed, and MinCUSUM variants produce coarse rectangular regions misaligned with the truth, whereas ST-Score traces irregular and disconnected change-region contours without prior shape information. An ablation on λ08 reveals a tradeoff: small radii yield sparse, uninformative neighborhoods while large radii contaminate scores with mixed-regime events; λ09 balances both metrics in the synthetic setting. Performance is robust to the pre-change base rate λ10 spanning two orders of magnitude—an observation the authors explicitly flag as contrasting with their EDD upper bound, whose dependence on event rate is not borne out empirically. Increasing optimization epochs λ11 from 1 to 5 improves both metrics, and runtime stays below one second per trial throughout.
Two case studies exercise the online variant without post-change reference data. On the JMA earthquake catalog for the Noto Peninsula (2018–2024), the detection statistic rises sharply in late 2021 and mid-2022 during the documented seismic swarm preceding the 2023 λ12 6.5 and 2024 λ13 7.6 events, and the estimated change region converges onto the main-shock epicentral area after swarm onset—consistent with swarm activity serving as a precursor signal. On CAL FIRE damage-inspection data for the January 2025 Palisades wildfire, the statistic drifts upward only at ignition: no precursor was detected, which the authors attribute to the abrupt, possibly human-initiated nature of wildfire ignition. The method still correctly localizes the affected area post-ignition despite unrelated hotspot activity earlier in the record.
Limitations and open questions
Several limitations are stated plainly in the paper. The theory assumes exact score models, and the drift analysis in closed form covers only translation-invariant (e.g., exponential or Gaussian) excitation kernels; extending parametric drift characterizations to broader kernel families remains open. The delay bound rests on concentration assumptions that are verified for Poisson processes but only conditionally established for Hawkes processes. Hyperparameters λ14 and λ15 materially affect performance—the ablations show non-monotone sensitivity—yet no principled selection rule or automated calibration procedure is provided, and linking them theoretically to detection performance is left as an open problem. Finally, the empirical robustness of EDD to event frequency stands in apparent tension with the derived bound, which the authors do not resolve.
Conclusion
This paper formulates sequential change-point detection for STPPs as a joint temporal-spatial inference problem over continuous domains and solves it with a score-based CUSUM statistic built from localized, conditionally weighted Hyvärinen scores trained via denoising score matching. The alternating optimization makes the nested maximization tractable, the theoretical analysis delivers false-alarm, delay, and Jaccard-index guarantees, and experiments on synthetic Hawkes data plus earthquake and wildfire catalogs demonstrate accurate region localization where binned or discretized baselines fail. The main open issues concern kernel-general drift analysis, principled hyperparameter selection, and reconciling the empirical insensitivity of delay to event rate with the theoretical scaling.