Recurrence analysis of extreme event like data

The identification of recurrences at various time scales in extreme event-like time series is challenging because of the rare occurrence of events which are separated by large temporal gaps. Most of the existing time series analysis techniques cannot be used to analyse extreme event-like time series in its unaltered form. The study of the system dynamics by reconstruction of the phase space using the standard delay embedding method is not directly applicable to event-like time series as it assumes a Euclidean notion of distance between states in the phase space. The edit distance method is a novel approach 5 that uses the point-process nature of events. We propose a modification of edit distance to analyze the dynamics of extreme event-like time series by incorporating a nonlinear function which takes into account the sparse distribution of extreme events and utilizes the physical significance of their temporal pattern. We apply the modified edit distance method to event-like data generated from point process as well as flood event series constructed from discharge data of the Mississippi River in USA, and compute their recurrence plots. From the recurrence analysis, we are able to quantify the deterministic properties of extreme 10 event-like data. We also show that there is a significant serial dependency in the flood time series by using the random shuffle surrogate method. Copyright statement. TEXT

. Climate change has already influenced river flood magnitudes (Bloschl et al., 2019), and has been related to increases in the intensity and frequency of heavy precipitation events aggregating flash flood and river flood risk (Donat et al., 2016;Kemter et al., 2020). Natural climate variability at different time scales may lead to flood-rich and flood-poor periods (Merz et al., 2018). In addition, human interventions in river systems and catchments also heavily influence flood magnitudes and frequencies (Hall et al., 2014). Since floods can massively affect life quality of our societies, it is desirable to understand the underlying dynamics and, thus, put forward precautionary measures to avert potential disasters. 25 The occurrence of extreme events is not random but rather a manifestation of complex dynamics, and such events tend to have long-term correlation. Comparing data of extreme events with other non-event data using standard methods (linear/nonlinear) is problematic, because the temporal sampling differs largely (time points of events vs. continuous sampling). Furthermore, linear methods such as the Fourier transform (Bloomfield, 2004) and wavelet analysis (Percival and Walden, 2007) are often insufficient to capture the full range of dynamics occurring due to the underlying nonlinearities (Marwan, 2019). Hence, 30 defining a principled nonlinear method is necessary for the analysis of extreme event time series, in particular when the correlation or coupling between several variables is investigated.
Out of the various approaches to study nonlinear dynamical systems (Bradley and Kantz, 2015), the reconstruction of the phase space using delay coordinates (Takens, 1981) is a widely used method that allows us to estimate dynamical invariants by constructing a topologically equivalent dynamical trajectory of the original (often high-dimensional and unknown) dynamics 35 from the measured (scalar) time series. In the delay coordinate approach, the distance between states in the phase space plays a pivotal role in describing the underlying dynamics of a system. After reconstructing the dynamical trajectory, we can extract further information about the dynamics of a system encoded in the evolution of the distances between the trajectories, e.g., through recurrence plots (Marwan et al., 2007), correlation dimension (Grassberger and Procaccia, 1983b), Kolmogorov entropy (Grassberger and Procaccia, 1983a), or Lyapunov exponents (Wolf et al., 1985). 40 Although there are powerful techniques based on phase space reconstruction of a wide range of nonlinear dynamical processes, they are not directly applicable to event-like time series. Extreme events such as flood, earthquakes, or solar flares are known to have long-term correlations (Jentsch et al., 2006). However, capturing the correlations of extreme events using such methods is difficult as the phase space reconstruction and the Euclidean distance for measuring the distances of states are not suitable for event-like time series because, by definition, extreme events are small in number and are separated by large 45 temporal gaps. In this case, it becomes necessary to define an appropriate distance measure that can help analyze the dynamics of extreme event-like time series. Event-like time series can be analyzed in their unaltered form by considering a time series of discrete events as being generated by a point process. Victor and Purpura (1997) presented a new distance metric called edit distance to calculate a distance between two spike trains (binary event sequences) as a measure of similarity. The method has also been adopted to 50 measure a distance between marked point processes to analyze foreign currencies (Suzuki et al., 2010) and irregularly sampled palaeoclimate data (Eroglu et al., 2016). Although the existing definition of edit distance is quite suitable to measure similarity between event series, it introduces a bias (discussed in Sec. 3), when there are large gaps in the data as in extreme event time series. Additionally, the method depends on multiple parameters and often it is difficult to associate a physical meaning to them. This further complicates the parameterization of the method in case of extreme events 55 In this study, we propose a modification of the edit distance metric for analyzing extreme event-like time series. The proposed extension allows to consider the shifting parameter of the edit distance metric in terms of a temporal delay which can be physically interpreted as a tolerance introduced to deal with the quasi-periodic nature of a real world extreme event time series. We demonstrate the efficacy of the proposed modified edit distance measure by employing recurrence plots and their quantification for characterizing the dynamics of flood time series from the Mississippi river in the United States. Flood 60 time series shows a complex time-varying behavior. Moreover, flood generation is often characterized by nonlinear catchment response to precipitation input or antecedent catchment state (Schröter et al., 2015), requiring methods able to deal with nonstationarity and nonlinearity. By using the random shuffle surrogate method, we show that there is a significant serial dependency in the flood events.
2 Edit Distance for event-like time series 65 Distance measurements between two data points play an important role for many time series analysis methods, for example, in recurrence plot analysis (RP) (Marwan et al., 2007), estimation of the maximum Lyapunov exponent (Rosenstein et al., 1993), scale-dependent correlations (Rodó and Rodríguez-Arias, 2006), data classification (Sakoe and Chiba, 1978), and correlation dimension estimation (Grassberger, 1983). In case of regularly sampled data, the Euclidean distance is often used. However, for event-like data where big gaps between events are common, this approach is not directly applicable. 70 Victor and Purpura (1997) introduced a method, called edit distance, to calculate a distance between spike trains. Later, Suzuki et al. (2010) extended this method to include changes in the amplitude of events. Ozken et al. (2015) suggested a novel interpolation scheme for irregular time series based on the edit distance, and Ozken et al. (2018) extended this approach to perform recurrence analysis for irregularly sampled data.
In order to apply the edit distance as a distance measure for, e.g., recurrence analysis, the whole time series is divided into 75 small, possibly overlapping, segments (windows), which should contain some data points (Fig. 1). To transform one segment to another, three elementary operations are required: (1) delete or (2) insert a data point, and (3) shifting a data point to a different For the edit distance we use the sum of the costs of the three operations. The cost function is where α and β are events in segments S a and S b occurring at times t a (α); C is the set containing all such pairs (α, β) from the two segments; I and J are the sets of indices of events in S a and S b ; | I |, | J | and | C | are the cardinalities of I, J and C; Λ s is the cost of deletion/insertion and Λ 0 the cost assigned for shifting events in time. The first summand in the cost function deals with deletion/insertion and the second summand (the summation) deals with shifting of the pairs (α, β) ∈ C.

85
The distance D is the minimum cost needed to transform the event sequence in S a to the event sequence in S b The operation costs Λ s and Λ 0 are defined as Λ s = const., and (2a) where M is the number of all events in the full time series. Originally, Λ s was chosen as a constant value 1 (Victor and Purpura, 1997). Later, Λ s was optimized in a range [0, 4] (Ozken et al., 2018). The units of the parameters are (time) −1 for Λ 0 and unit-less for Λ s .
The treatment of an event series as point process makes the edit distance measure a good starting point for defining a distance 95 between segments of an extreme event time series. However, the existing form of the edit distance has a linear dependency on the difference between the occurrences of events which is inappropriate for an extreme event time series, as the rare occurrence leads to large gaps between events. Also, as already mentioned, the existing method depends on a number of parameters.
Therefore, we suggest a modification of the cost function to address these two concerns.

100
The cost of transformation between related events (e.g. events belonging to the same climatological phenomenon in a climate event series) should be lower than that between independent, unrelated events. Hence, the shifting operation should be a more likely a choice for comparison between segments if the events in each segment are related or belong to the same phenomenon.
Consequently, deletion/insertion as a choice of operation for transformation tends to be associated with unrelated events. Now, we consider two event sequences S a and S b , where each of them has only one event at times t a and t b , respectively 105 ( Fig. 3). We want to transform segment S a to S b by using either deletion/insertion or shifting. The path for shifting the solitary spike has the cost Λ 0 ∆t where ∆t =| t a − t b |, i.e., it grows linearly with the distance between the two events. On the other hand, the cost to delete the event at t a from S a and insert it at time t b in S a for it to resemble S b has cost 2, as the single operation cost for deletion and insertion is each 1. So the shifting operation will take place only as long as the time difference between both events is smaller then ∆t < 2 Λ0 . The above condition has two limiting cases which need to be considered for an unbiased choice between shifting and deletion/insertion -when the cost of shifting is too low or when it is too high. The former case arises for an extreme event time series where the number of extreme events M total length of the time series, i.e., Λ 0 → 0, Eq. (2b). This assigns a low cost for the transformation due to a biased choice of shifting over deletion/insertion for largely separated events which may be independent.

115
The other limiting case, when the cost of shifting per unit time, Λ 0 , is moderately high, leads to a biased preference for the deletion/insertion operation over shifting depending on how high Λ 0 is, because the cost of shifting increases linearly with the distance between two events according to the definition in Eq. (1). In this case, related events separated by a relatively small gap may be considered as independent as the shifting cost may exceed the deletion/insertion cost because of higher Λ 0 . The following example of a real extreme event series illustrates how this biased preference can lead to erroneous results. In regions 120 with seasonal precipitation regimes, high precipitation events are more similar to those of the same season of another year compared to other seasons. However, the exact time of occurrence of the extreme precipitation events varies from year to year, i.e., even if the events in each segment are related, there might be a certain time delay between the events (Fig. 4). Now, a linear dependency of the shifting cost on the difference between times of events, Eq. (1), does not allow us to consider small delays between potentially similar events. In addition, a not so high value of Λ 0 could increase the cost of shifting for small delays 125 higher than the cost of deletion/insertion, implying a lower similarity between the segments. As the maximum cost for shifting is limited by the cost for deletion/insertion, we need to lower the cost of shifting to get a higher similarity between the segments for the above case. This modification is done by considering the temporal distance ∆t at which the maximum cost of shifting occurs as a delay between the two events.
Selecting a certain temporal delay is relevant for climate time series, such as precipitation, or hydrological variables, like 130 discharge, which are of an event-like nature, where the synchronization between extreme events (Malik et al., 2010(Malik et al., , 2012 6 https://doi.org/10.5194/npg-2020-41 Preprint. Discussion started: 31 October 2020 c Author(s) 2020. CC BY 4.0 License. cost Δt Figure 5. Linear cost function for shifting (black) and alternative, nonlinear cost function (orange). Boers et al., 2013Boers et al., , 2014Ozturk et al., 2018) at different geographical locations and recurrence of extreme events for the same location are of interest.
The above discussion illustrates the importance of an appropriate choice of the cost of shifting per unit time, Λ 0 , based on properties of the event series, such as time series length, total number of events, and event rate. Thus it is more reasonable to 135 optimize Λ 0 as opposed to Λ s . Victor and Purpura (1997) suggested generalizations of the cost assigned to finite translations, such as modifying the simple Λ 0 ∆t to more general functions, satisfying the triangle inequality.
In view of comparing two events under a predefined delay, we suggest to replace the linear cost function for shifting, Eq. (1), by a nonlinear cost function which allows a temporal tolerance and also ensures a smooth change from the cost function for shifting to the cost for deletion/insertion (Fig. 5). We propose to use the Logistic function (Cramer, 2002) as the cost for shifting: where τ is the chosen time interval between events marking the transition of the function from exponential growth (for closelyspaced events) to bounded exponential growth (for events separated by large time gaps) (Fig. 6a),

145
-Λ s is the maximum cost, and k is a parameter that affects the rate of the exponential growth, in this study we choose k = 1. For k = ∞, the Logistic function becomes a step function.
The Logistic function (Eq.(3)) is used for modeling the population growth in an area with limited carrying capacity. The sigmoid nature of the Logistic function is apt for representing the cost of shifting with low values for closely-spaced events 150 and which tends to saturate for events separated by a time delay greater than a certain τ .
A step function (e.g., Heaviside function) would be another choice for the cost function according to whether shifting is chosen below a certain delay and deletion/insertion is chosen above it (Fig. 6b): where Λ 0 is the cost for shifting, according to Eq. (2b), and τ is the given delay choice which decides between shifting and 155 deletion/insertion. However, such a cost function maintains a constant value irrespective of the event-spacing and switches to a different value only at the delay threshold. Therefore, the logistic function, as a smooth approximation of the step function, is a preferable choice.
The parameter optimization problem of the modified edit distance equation now becomes a problem of solving a single linear equation with two unknown variables -the coefficient related to the cost of shifting events in time and Λ s . Since, in our 160 method, we choose to optimize the former by using the Logistic function, we keep Λ s constant. This constant can be absorbed in the optimization function of the first term, and therefore, the coefficient of the second term related to the cost of addition/ deletion.
For a particular pair of events (α, β), we can arbitrarily set the maximum cost Λ s for shifting or deletion/insertion to be 1, because it is the only free cost parameter (we have no Λ k because of neglecting the amplitudes). Thus, the cost for transforming 165 one segment to another using the modified edit distance (mED) is defined as and the corresponding distance function to calculate the distance between two segments S a and S b is This new definition of cost depends only on the parameter τ , which can be interpreted in the sense of a delay between events.

170
Moreover, this definition holds the triangular inequality when τ satisfies a certain condition (see the Appendix for proof).
In the next section we use the modified version of edit distance to perform recurrence analysis of extreme event-like data.

Recurrence plots
A recurrence plot (RP) is a visualization of the recurrences of states of a dynamical system, capturing the essential features of the underlying dynamics into a single image. The quantification of the patterns in a RP by various measures of complexity 175 provides further (quantitative) insights into the system's dynamics (recurrence quantification analysis, RQA) (Marwan et al., 2007). RQA was designed to complement the nonlinearity measures such as the Lyapunov exponent, Kolmogorov entropy, information dimension, and correlation dimension (Kantz, 1994). RP-based techniques have been used in many real world problems from various disciplines. In the field of finance, Strozzi et al. (2002)  activity in the last 150 years. It has also been used to classify dynamical systems (Corso et al., 2018), or to detect regime changes Eroglu et al., 2016;Trauth et al., 2019) and has many applications in Earth science (Marwan et al., 2003;Chelidze and Matcharashvili, 2007), econophysics (Kyrtsou and Vorlow, 2005;Crowley and Schultz, 2010), physiology (Webber, Jr. and Zbilut, 1994;Zbilut et al., 2002;Marwan et al., 2002;Schinkel et al., 2009), and engineering (Gao et al., 2013;Oberst and Tuttle, 2018).

185
Consider a time series that encodes N measured states of a dynamical system {x i ∈ M } N i=1 on space M . A state of this system is said to be recurrent if it falls into a certain ε-neighborhood of another state. Given a distance function D : M × M → {0} ∪ R + , for a given trajectory x i , the recurrence matrix of the system is defined as: In a RP, when R i,j (ε) = 1 a point is plotted at (i, j), otherwise nothing is plotted. Different classes of dynamics result in 190 different patterns in their respective RPs (Marwan et al., 2007).
The RP contains a main diagonal line, called the line of identity (LOI), corresponding to the recurrence of a state with itself.
The RP is symmetrical about the LOI when D is a metric (e.g., a symmetrical norm). We are interested in the line structures in RP as they capture several aspects of the underlying dynamical behaviour of a system. For instance, long, continuous lines parallel to the LOI denote (pseudo-)periodic behavior, whereas short, discontinuous diagonal lines are indicative of a chaotic 195 system.
In order to incorporate mED in RP, we divide the time series into small segments and compute the distance between these segments; the time indices are the centre points of these segments. We can now define the RP in terms of the distance between the segments calculated by mED: One of the most important measures of RQA is the determinism (DET), based on the diagonal line structures in the RP. The diagonal lines indicate those time periods where two branches of the phase space trajectory evolve parallel to each other in the phase space. The frequency distribution P (l) of the lengths of the diagonal lines is directly connected to the dynamics of the system (Marwan et al., 2007). The ratio of the recurrence points which form diagonal lines to all points on the RP is the RPs of stochastic processes mainly contain single points, resulting in low DET values, where RPs of deterministic processes contain long diagonal lines, resulting in high DET values.
In this work, we are interested in the deterministic nature of flood event time series. For this purpose, we focus on the DET 210 measure.

Choice of window size
We divide the time series into segments/windows. If adjacent windows overlap, we call it overlapping sliding window, else, non-overlapping window (Fig. 1). Sliding window techniques are widely used for signal processing (Bastiaans, 1985), activity recognition processes (Dehghani et al., 2019). The window size should be selected properly, i.e. each window should contain 215 enough data points to be differentiable from similar movements.
Consider a time series of data x i ∈ R at times t i ∈ N. For generalization we consider a constant sampling rate, i.e., Each window consists of n (n ∈ N, n > 1) data points, so the window size is w = n∆T . For overlapping windows, a fraction of the data is shared between consecutive windows, denoted by L ∈ {1, 2, 3, .., n − 1} as the number of data points within the 220 overlapping range, where L = ∅ signifies the non-overlapping case. The overlapping range is OL = L∆T and in terms of percentage OL(%) = L n . The number of data points being shared for overlapping windows is p = w∆T − OL. Because mED focuses on event-like data, and finding an optimum window size is important. An inappropriate window size may lead to missing events in windows ("empty windows") due to the sparse occurrence of events. Here we propose some criteria for choosing the optimum window size:

225
-To avoid empty windows we try to fix the number of events (n) in each window. Here we choose the window size as 1 year, as most climate phenomena such as flood events exhibit annual periodic behaviour and therefore we get similar number of events each year.
-In case we need to choose a larger window size to have enough data points in each window, the number of segments decreases if the windows do not overlap which in turn reduces the dimension of the recurrence matrix. The underlying 230 dynamical behaviour of the system will not be completely captured by the resulting coarse-grained recurrence plot. Overlapping windows alleviate this problem by increasing the number of windows. In our work, we use a certain percentage of overlapping.
-An inappropriate window size can lead to the aliasing effect (Fig. 11). As a result different signal would be indistinguishable and we might loose important transitions.

6 Applications of modified edit distance
We apply mED first to generated data and then to real-world flood observations to understand how well our method can identify recurrences in extreme event-like time series.

Numerically generated event-series using Poisson process
The Poisson process is used as a natural model in numerous disciplines such as astronomy (Babu and Feigelson, 1996), biol-240 ogy (Othmer et al., 1988), ecology (Thomson, 1955), geology (Connor and Hill, 1995), or trends in flood occurrence (Swierczynski et al., 2013). It is not only used to model many real-world phenomena but also allows for a tractable mathematical analysis.
Consider N (t) to be a stochastic counting process which represents the number of events above some specified base level in the time period (0, t). Suppose the events occur above the base level at a constant rate λ > 0 (units of 1/time). So, the 245 probability that n events occur in the time between t and t + s is given by (Ross, 1997;Loaiciga and Mariño, 1991) P {N (t + s) − N (s) = n} = e −λt (λt) n n! (n = 0, 1, 2, . . .).
A Poisson process is a set of random events whose stochastic properties do not change with time (stationary) and every event is independent of other events, i.e., the waiting time between events is memory-less (Kampen, 2007).
Here, we study three cases of numerically generated event series which are motivated by the occurrence of natural events.

250
First, we test mED for a simple Poisson process. Each subsequent case tries to capture features of these real world phenomena by adding an element of and memory to the simple Poisson process. We compare the RPs and the structures in the RPs (by considering the RQA measure DET) derived from the standard edit distance and from the modified edit distance.
For ED (Eq. 1), the cost parameter for shifting is calculated according to Eq. (2a) and Λ s = 1, whereas Eq. (4) is used for mED, recurrence threshold ε = 10 percentile of the distance matrix for each case. We choose the upper bound for the range of 255 τ to be less than the mean inter-event time gap for the complete time series. The physical interpretation of this choice is that the temporal tolerance or delay in the arrival of events in a particular season (event cluster) should not only be less than the time period of the seasonal cycle but also relatively less than the length of the season.
It is expected that the quasi-periodic behaviour of the event series will lead to high DET values. In case of ED, the DET should be constantly high at all τ , as ED does not include the concept of time delay. On the other hand, DET computed using

Homogeneous Poisson process
The homogeneous Poisson process models the occurrence of random spaced, i.e. stochastic, events in time, where the average time between the events is known. It is used to model shot noise, radioactive decay, arrival of customers at a store, earthquakes, etc. Here we construct an event series (Fig. 7a) from Eq. (9) with λ = 1 5 . The RP of a realisation of such a homogeneous

265
Poisson process is shown in Fig. 8a,b. In a stochastic process, a recurrence of a randomly selected state occurs by chance, resulting in randomly distributed points in the RP. Accordingly, DET has low values (Fig. 8c).

Repeating Poisson process
We generate a small segment of a simple Poisson process and repeat the same segment after certain time gaps (Fig. 7b). Here, the Poisson process is generated for a segment of length 50 with mean rate of events λ = 1 5 . The gap between the event segments chunks make the diagonal lines discontinuous implying a quasi-periodic behaviour. Quasi-periodicity is often a characteristic of an extreme event time series such as flood events. RP using ED Fig.9b contains more short and discontinuous diagonal lines corresponds to less DET value. In Fig.9c, we find a certain range where the DET value calculated using mED is higher than 275 ED.
(a) (b) (c) Figure 9. RP of a repeating Poisson process (Fig.7b) using (a) mED, (b) ED; (c) comparison of DET for 300 realisations using ED (blue horizontal line) and using mED for varying τ in range 0 to 15.

Poisson process with periodical forcing
Quasi-periodicity in a time series can be expressed as a sum of harmonics with linearly independent periods. Here we construct a time series with superposition of a slow moving signal f (T 1 ) = sin(2πt/T 1 ) and a fast moving signal f (T 2 ) = sin(2πt/T 2 ), whereas the events are drawn by a Poisson process (Fig. 7c). We select the extreme events by using a certain threshold and 280 13 https://doi.org/10.5194/npg-2020-41 Preprint. Discussion started: 31 October 2020 c Author(s) 2020. CC BY 4.0 License. take the peak over the threshold (POT). Their occurrences and positions can be seen as the outcome of a point process (Coles, 2001). Such a time series can be used to mimic a streamflow time series which may exhibit periodic behaviour on several scales (annual, seasonal, decadal, etc.). Here the window size is equal to the time period of the slow moving signal. The RP is shown in Fig. 10. The periodic occurrences of the event chunks is visible in the RP by the line-wise accumulation of recurrence points with a constant vertical distance (which corresponds to the period). The stochastic nature of the "local" event pattern is 285 visible by the short diagonal lines of varying lengths. As mentioned earlier, when using an improper window size we can lose recurrence points and a smaller number of diagonal lines occurs due to the aliasing effect (Fig. 11). For a certain range of τ , DET for mED (Fig.8, Fig.9, and Fig.10) to be higher than the DET for ED, thus capturing the underlying periodic behaviour better.

Recurrence analysis of flood events 290
As a case study, we apply our method to the Mississippi river which has a rich history of flood events. Here we are motivated to study the recurring nature of flood events using recurrence plots. Recurrence plot analysis also helps to quantify the serial dependency of flood time series. (Wendi et al., 2019) used RQA to study the similarity of flood events.
We use the mean daily discharge data of the Mississippi River from the Clinton, Iowa station in the USA for the period 1874-2018. The data is obtained from the Global Runoff Data Centre (www.bafg.de/GRDC).

295
We construct the flood event series from the complete discharge time series as follows: We first apply a threshold (corresponding to the 99 th percentile) for each year on the discharge values which gives a few events (say 3/4) per year. If several successive days fall above the threshold, we pick only the maximum event from this cluster. Then we lower the threshold by 0.1 percentiles (the threshold typically lies between 99-90 percentile) and repeat the above step until we get the desired number of independent events. 300 Figure 11. RP of Poisson process with periodical forcing (Fig. 7c) using mED showing aliasing effect due to improper window size (note that this effect also occurs for the standard ED approach).
We apply mED for finding recurrence patterns for flood events. To this end, we choose our window size equal to the annual cycle (1 year) with 6 months and 9 months of overlapping and with the cost function of delay, τ = 15, 30, 45 and days. The recurrence threshold ε = 8 percentile of the distance matrix. The recurrence plot of flood events is shown in Fig. 12 and 13 .
We also compute the recurrence plot using ED and calculate DET Fig. 14.
The black points in the RP (Fig. 12, 13) denote the segments which have higher similarity and the white points imply less  Fig. 10a for events generated by Poisson process with periodical forcing. Diagonal structures are seen in the flood RP ( Fig. 12, 13) denoting an inherent serial dependency in the data as indicated by the DET value. Serial dependence is a property by virtue of which the future depends on the past. When diagonal lines are likely to appear in a RP, the current 310 neighbors tend to be neighbors in the near future, and thus, the corresponding system is said to have serial dependence.
To determine the statistical significance of the RP analysis, we develop the following statistical test. We use random shuffle surrogates (Scheinkman and LeBaron, 1989) statistical testing. Suzuki et al. (2010) used this method for finding serial dependency on foreign exchange tick data by quantifying the diagonal lines in the RP (DET measure).
We first set a null hypothesis, then we generate a set of random surrogates that preserve the null hypothesis property. After 315 that, we compare the test statistics of the original data with the surrogates data. We can reject the null hypothesis if the test statistics obtained from the surrogate data and the original data is out of the specified range, else we cannot reject the null hypothesis.
The null hypothesis is that there is no serial dependency in the data. Thus, we expect that after each random shuffling the information will be preserved. We create each random surrogate by randomly shuffling rows and columns of a recurrence plot 320 simultaneously. We can reject the null hypothesis if the test statistics of the original data K 0 (DET value) is out of the specified range from the distribution of those surrogate datasets K S (Theiler et al., 1992). The algorithm works as follows: 1. Calculate the distance matrix and compute the RP for a certain window size.
2. Reorder the row and columns randomly to create surrogate recurrence plots.
3. Calculate the DET of the surrogate recurrence plots. We assume the statistics obtained from the surrogates K s is normally distributed and we consider only one experimental data set. The measure of "significance" as defined by Theiler et al. (1992): follows a t-distribution, the number of surrogates are n = 50 with (n − 1) degrees of freedom. µ Ks and σ Ks are the mean and 330 standard deviation of the statistics of surrogate data. The value of ρ has to be compared to the value of the t-distribution that corresponds to the 99 th percentile and the mentioned degrees of freedom, i.e., t 0.01/2 (49) = 2.68. If ρ exceeds this value, the null-hypothesis has to be rejected. We measure ρ for the data generated from a homogeneous Poisson process, Eq. (9), and also for the flood event data. First, for the Poisson process the value of ρ = 1.91, < 2.68, so we cannot reject the null hypothesis, confirming the missing serial  For the flood event data, the distribution of the test statistics (DET) is far away from the value of the original data (Fig. 15).
The ρ value for the flood event RP with 6 months overlapping and τ = 15 days is 23.80 which is much larger than 2.68. Hence, we can reject the null hypothesis at the significance level of α 0 = 0.01, getting strong indications that there is serial dependency at high significance level in the occurrence of flood events.

340
In this paper, we propose a distance measure for recurrence based analysis of extreme event time series. The proposed measure is based on a modification of the edit distance measure proposed by Victor and Purpura (1997). We include the concept of time delay to incorporate the slight variation in the occurrence of recurring events of a real world time series due to changes in seasonal patterns by replacing the linear dependency of the cost of shifting events by a nonlinear dependency using the logistic 345 function. This is substantial improvement over the previous definition of the edit distance as used for the TACTS algorithm proposed by Ozken et al. (2015) as the optimization of Λ 0 is based on the temporal delay between events, which is of physical relevance to the study of extreme events and can be chosen according to the phenomenon being studied. The modified edit distance also reduces the number of independent parameters. We tested mED on prototypical event series generated by a point process (Poisson process). We found that for the quasi-periodic repeating Poisson process event series, the determinism of RP 350 computed using mED varies with temporal delay and is higher than the measured DET value of RP computed using ED for a certain range of delay. We applied the method to study recurrences in the flood events of the Mississippi river. Our analysis revealed deterministic patterns in the occurrence of flood events from RP. Finally, using the random shuffle surrogate method we have shown that the data of the occurrence of flood events has a statistically significant serial dependency.
In this work, we have only considered binary extreme event time series, and ignored the amplitude of events. Next, the mED 355 measure should be further modified by including the cost due to difference in amplitude for applications in real world time series where patterns of intensity of extremes might be of interest. Moreover, an elaborate method to optimize for a proper window size may be devised to capture recurrences of events at different time scales.

Appendix A
For deletions and insertions, the new definition inherits the triangle inequality for the original definition of edit distance by 360 Victor and Purpura (1997). Thus, we can show that the new definition of the edit distance preserves the triangle inequality as a total if a shift preserves the triangle inequality. For this sake, considering the following situation is sufficient: Suppose that there are three simple point process each of which has an event; the first point process has an event at u(u > 0); the second point process has an event at u + s(s > 0); the third point process has an event at u + s + t(t > 0). Therefore, the first and the second point processes have an inter-event interval of s; the second and the third point processes have an inter-event interval 365 of t; and the first and the third point processes have an inter-event interval of s + t.
We assume that τ ≤ s+t 2 . Then, the distance between the first and the second point processes is D(S 1 , S 2 ) = 1 1 + e −s+τ .
Likewise, the distance between the second and the third point processes is D(S 2 , S 3 ) = 1 1 + e −t+τ .

370
The distance between the first and the third point processes is D(S 1 , S 3 ) = 1 1 + e −s−t+τ .
Then the following chain of inequalities holds: