The study of event time series is of general interest in data analysis and modeling, because of their ubiquitous nature in almost all scientific fields, such as investigating financial transactions, customer interactions, life-threatening cardiac events, system failures, or natural phenomena. Event series can be single, discrete events, binary events or events with different amplitude, e.g., events extracted from data with heavy tail distributions, short-term extreme events, or anomalies in time series. In neuroscience, event series are also called “spike trains” [1]. A time series is generally be denoted by a set of ordered pairs {(t_i, x_i)} of time t_i with t_i+1 > t_i and corresponding data value x_i; and with sampling index i (usually constant sampling time t_i+1 − t_i = const., i.e., equidistant time axis). An event series, instead, is considered as a series of event times, defined by the associated specific time or timestamp of the single events, finally resulting in a set of event time points {t_i}. As events could also have some amplitude, a definition as an event time series as a tuple of time and event strength {(t_i, x_i)} is also possible. Because the events usually do not occur at regular intervals, such event time series are usually irregularly sampled t_i+1 − t_i ≠ const. The alternative is using a regularly sampled, discretized time axis with binary (or amplitude) values at those points of time where the event happens (this is similar to categorical data, another class of discrete data, but not necessarily representing separated single events). However, this approach is usually limited and not appropriate for many research questions, because the timing of events often does not fit the sampling points and, even more important, the time series can be filled with many zeros. Standard time series analysis tools have their limits when analysing such data.

Examples of event data are time series of extreme events, which are of specific interest because of their technical and medical importance, and their potential of serious societal impacts: Extreme rainfall (flush floods) and river floods are of strong concern because they are increasing worldwide due to the global warming [2, 3]. Extreme loading conditions are considered and modeled in material sciences to monitor and predict serious failures, e.g., on bridges caused by extreme traffic or on airplane structures due to sudden stress or birdstrikes [4, 5]. Ventricular tachyarrhythmias are life-threatening cardiac arrhythmias, usually analyzed by investigating the beat-to-beat intervals of the heart, expressed by a series of events [described as the R-waves in an electrocardiogram (ECG)], themselves [6]. Examples of natural, rare event time series are sequences of landslides. Such landslide events cause serious damages and can be triggered by specific weather phenomena, like atmospheric rivers or El Niño/ Southern Oscillation [7, 8]. They are increasingly affecting urban settlements, because of the spreading of cities and climate change [9]. Another example is brain activity which is controlled by the firing of neurons. For example, the coherence of neuron firings can cause seizures [10]. The investigation of extreme events in dynamical systems is an important subject in statistics and statistical physics. It covers many research questions, from the emergence of extremes to predicting extreme events [11].

The research questions related to event series are often the same as for other kind of data, e.g., comparing different time series, classifying the dynamics of the process behind, identifying regime changes, or use as the base for simulations and predictions. The growing availability of data and computational facilities in almost all scientific disciplines has significantly advanced data science in general. Several approaches have been introduced that allow to study different research questions related of event data [12–14]. Among them are probabilistic methods based on large deviation and extreme value theory (parametric, semi-parametric approaches, and multivariate extensions), pattern-based prediction algorithms and BDE modeling [15–17], as well as modern learning based approaches for predicting extreme events [18, 19]. Another class of methods are based on the property of recurrences of states. In general, recurrence based methods provide versatile approaches for classifying data, identification of regime transitions, but also for unveiling interrelationships, synchronization, and causal links between different dynamical systems [20–23]. Due to its broad usability, simplicity, and growing number of software allowing recurrence analysis [24], this method is attracting more and more interest and becoming increasingly popular [23, 25]. By modifying the definition of what a recurrence is, it is, in general, possible to adapt recurrence analysis to be usable for analysing discrete and event-like data [26–28]. In the following, I describe briefly the basics of recurrence analysis, its extension to work with event data, and the related challenges and future perspectives.

A recurrence plot (RP) is the graphical representation of the recurrence matrix, which is simply representing all pair-wise time combinations (i, j) of a data sequence which have similar values or states x→i:

with a similarity measure d(·, ·) and the Heaviside function Θ(·) which sets R_i,j = 1 if the similarity value d(·, ·) falls below the threshold ε [29]. For dynamical systems with continuous change of the state variables, i.e., x→i∈ℝm (with m the dimension of the system), the similarity between states is often defined by the Euclidean norm d(xi→,xj→)=‖x→i-x→j‖ [29]. For discrete data of regular sampling (equidistant time instances), e.g., categorical data, the recurrence matrix R can be simply defined by the Kronecker delta function R_i,j = δ(x_i, x_j), which is one if both arguments are identical [26, 30–32]. This approach works well for discrete data, such as categorical data or symbolic sequences, with applications, e.g., in life science to detect atrial fibrillation or congestive heart failure [33, 34], to measure synchronization in an epileptic brain [30], or in engineering to optimize manufacturing networks [35]. This concept is easily extendable for bivariate analysis. Cross-RPs, CRi,jx,y=Θ(ε-d(xi→,yj→)), and joint-RPs, JRi,jx,y=Θ(ε-d(xi→,xj→))∘Θ(ε-d(yi→,yj→)), are two basic concepts for measuring different aspects of synchronization [29]. To modify, the cross-RP for discrete data, we can simply use the Kronecker delta CR_i,j = δ(x_i, y_j) [36]. Joint-RP even allows us to measure the synchronization between different types of data, such as discrete and continuous data [37], where

is the Hadamard product of the RP of the discrete system x_i and the RP of the continuous system y→i.

This concept reaches its limits when considering event time series which consist of rare events and many zeros between them, or, even more limiting, consist only of the events {t_i} or have strong non-equidistant time instances (t_i+1 − t_i ≠ const.). For this kind of data, the similarity measure d(·, ·) has to be replaced by a specific metric which measures the coincidence of event sequences. Several measures (event metrics) are available, mainly developed in neuroscience [13]. A widespread measure would be the event synchronization, which allows varying delays between events to be considered as coinciding [38]. This measure is successfully applied for investigating, e.g., the spatio-temporal relationships between extreme rainfall events [39]. Another candidate is the edit distance, an extension of the Levenshtein distance [40, 41]. The distance is calculated by the minimum cost needed to modify one event sequence into another with a limited set of operations (Figure 1). Edit distance is a metric and has been successfully integrated with recurrence analysis [27, 28].

The edit distance measure is the (minimum) sum of the costs of the transform operations addition, deletion, and shifting applied to modify a sequence Si={t1(i),t2(i),…,tNi(i)} of N_i events (with events at time points ta(i)) into sequence Sj={t1(j),t2(j),…,tNj(j)} (with events at time points tb(j)):

where a and b are indices of the events in segments S_i and S_j; N_i and N_j the number of events in segments S_i and S_j, respectively; N_{(i, j)} the number of events in S_i and S_j to be shifted, which form the set C; λ_s is the cost of deletion/ insertion, and λ₀ the cost assigned for shifting events in time. Thus, the first summand corresponds to deletion and insertion operations and the second summand to the shifting of the events (Figure 1). Extensions of this cost function include considering costs for amplitude changes or to modify the shifting term by a continuous cost function allowing a more intuitive interpretation in terms of a delay [27, 28]. To apply the edit distance for recurrence analysis, the event series has to be divided into sequences S_i defined by a time window of length T_w. The shifting of this interval can be with smaller steps s < T_w resulting in overlapping time intervals. In order to get reliable costs d(S_i, S_j), the resulting sequences S_i and S_j should have at least one event (i.e., should not be empty).

This edit distance measure has been used as a simple synchronization measure between event series to study the stimulus responses in neuron spike trains [41], as a similarity measure between extreme rainfall data to reconstruct climate networks [42], and to create regularly sampled time series from non-regularly sampled time series (TACTS approach) [43, 44]. It was also used as a distance measure for computing RPs directly from event data (Figure 2), e.g., to study stock exchange data [27], flood events [28], or to allow calculation of RPs directly from irregularly sampled palaeoclimate data [46, 47]. The integration of the edit distance metric into the RP definition, Equation (1), provides all the applications of recurrence based time series analysis for the specific data of event time series.

Despite the recent advances in recurrence analysis of event time series, there are still several challenges.

Event time series can have missing data which are not easy to be detectable. For example, data on landslide events is mainly available at sites where they affect infrastructure [48], but their statistical analysis with respect to, e.g., climate change would require reliable event series [49]. Missing or sparse data can, therefore, bias the results of any analysis, and is subject of research in time series analysis in general, including interpolation, modeling, or advanced data reconstruction methods [50–52], but mainly not applicable for event data.

The process behind the analyzed study object could be non-stationary (e.g., life-threatening cardiac arrhythmias or seizures [6, 10]), meaning that the statistical properties of the event series may change over time (such as the distribution of events could change over time – events may be sparse, meaning that there could be some periods of time without events), which can make it difficult to apply the event based recurrence analysis (e.g., using edit distance). In particular, if the time interval defined by length T_w is too small, many sequences S_i could be empty, resulting in non-defined costs d(S_i, S_j). The selection of the time interval length T_w is, thus, crucial. For simple periodically recurring events, the choice might be easy, but its selection if multiple time scales are present is not straightforward [28].

The number of events in an interval can also change due to sampling issues, as it is a common problem in palaeoclimate research, where the deposition rate in sediments is varying over time, thus, leading to palaeoclimate time series of non-equidistant sampling in general [53–55]. Event based metrics, such as event synchronization, event coincidence analysis, or edit distance cost depend on the number of events in the interval and produce different types of biases which impact the results of the quantitative analysis and call for correction schemes [55, 56].

In general, the comparison of event time series with continuous time series is very challenging. Such problems occur, e.g., in climate research when studying recurring pattern of special weather phenomena (e.g., atmospheric rivers) or extreme events (such as heavy precipitation or river floods) with respect to large scale climate phenomena, such as El Niño/ Southern Oscillation or North Atlantic Oscillation [7, 57]. The RP approach offers a promising way by modifying Equation (2) to Kodama et al. [37]. However, event series often consist of much less events than the number of sampling points of the continuous time series, resulting in RPs of rather different length and making it impossible to directly apply Equation (4). An approach to finally match the event based RP with those of the continuous data would be required.

Finally, the uncertainty of the timing of events (timing jitter) is strongly affecting any measure of coincidence. It is expected that timing jitter is a common problem in measuring real-world event series. This challenge might be addressed by evaluating the sensitivity of the results on the jitter using specific models.

The extension of the edit distance can also take amplitude variations into account. However, this mixes two different aspects of the data: the temporal pattern of event sequences and amplitude differences. The optimal choice of the corresponding parameters might be less clear then, but have to be used to balance between these aspects.

The perspective future methodical developments will consider several important challenges to study interesting research questions related to (discrete) event data.

For recurrence analysis of event data, so far only the edit distance metric has been applied. It would be important to test and compare also other measures, such as Needleman-Wunsch distance, event synchronization, event coincidence analysis, or ARI-SPIKE-distance [13, 38, 56, 58]. Specific discrete data might call for distance metrics considering amplitude differences, e.g., edit distance or longest common subsequence [59].

Data with missing events is a general problem. Different strategies might be considered to solve this challenge, including correction and gap filling schemes [51, 55]. Correction schemes are also important for data with non-stationarities (varying sparsity of events). Such correction schemes needs further development to be applicable in a more general way (e.g., independent of the event distribution) and be more computationally efficient.

Events can exhibit some kind of temporal dependencies, meaning that the likelihood of an event may depend on the occurrence of previous events. The RQA measures could be used to study temporal dependencies from event series [28]. In general, diagonal lines in a RP represent the tendency that current neighbors in phase space will remain to be neighbors in the near future, thus represent serial dependence. The RQA measure determinism is quantifying the fraction of recurrences forming such diagonal lines and can, thus, be used as an indicator of serial dependence.

The classification of dynamical processes by event series based on duration, frequency of events or their characteristics (e.g., shape), will be another interesting application which will also involve machine learning approaches. The combination of machine learning with recurrence analysis is currently a strongly developing field with applications mainly in classification and prediction, using RPs and RQA measures as inputs in machine learning workflows [23]. A typical example is to convert time series into images by using the RP approach which are finally fed into the machine learning workflow for classification [60]. RPs of event series can be used in a similar manner for such kind of classification tasks. Other characteristics of event series (like serial dependence) would be accessible to machine learning approaches by the RQA measures [61–63].

The detection of interdependencies or synchronization of (sub-)systems represented by different kinds of data (e.g., event data with continuous time series) is an important methodical challenge. New approaches based on RPs seem to be promising, including the concept of joint-RPs [37] and the comparison of the probability of recurrences [64]. The advantage is the comparison by the recurrence structure, which would allow comparing time series of different kinds (e.g., event series vs. continuous data). It includes further developments to finally match the size of event based RPs with those of the continuous data, e.g., considering coarse-graining, interpolation, or specific window selections schemes (for event-sequence based metrics like edit distance) [28].

RP based analysis can be used to infer coupling directions or even causal links between different systems [65, 66]. Thus, the next step would be to test this approach for its potential on causality testing even for event data.

RPs also allow to identify patterns or regularities in challenging data, such as event series, including the estimation of the power spectral density of event series [67]. The most obvious way to derive a spectrum from a RP is to use the probability of recurrence after lag τ, which is simply the density of recurrence points along the diagonals (with distance τ from the main diagonal). This probability of recurrence is related to the auto-correlation [68, 69]. Using the edit distance measure, the temporal dependency structure within the event series can be visualized and quantified with this approach. Finally, the power spectrum can then be estimated from this probability of recurrence, either by applying the Fourier transform or any other advanced decomposition [67, 69].

The uncertainty of the timing of events (timing jitter) needs to be considered in the analysis, leading to new concepts such as Monte Carlo based ensemble approaches or Bayesian approaches. A recently proposed concept combines a Bayesian approach with RPs to derive a RP which explicitly represents the uncertainties of the timing of data points [70]. The resulting recurrence matrix contains the probabilities of recurrences instead of the binary information of recurrences. The recurrence quantification of such matrix is still subject of future research.

Although various distance measures for event based RP computation are available, the already applied one, edit distance, provides already a bunch of interesting directions for future research. For example, the choice of an optimal window length T_w or the different cost parameters λ. Including the cost for amplitude differences require an optimal choice of the corresponding parameters, which would need some systematic studies to provide some guidance to balance between the differences in the temporal and spatial domain.

The recurrence analysis as a concept is rather novel approach, with a lot of interesting and powerful developments and extensions in the last two decades [23]. It is also a promising concept for studying different aspects related to (discrete) event time series, where other methods have their limits.

Julia code to reproduce Figure 2 is available at: Zenodo, https://doi.org/10.5281/zenodo.7467886.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

The author confirms being the sole contributor of this work and has approved it for publication.