1 Introduction

Time series displaying temporally inhomogeneous behaviour have received strong interest in the recent statistical physics literature (Barabási 2005; Oliveira and Barabási 2005; Vasquez et al. 2006; Vazquez et al. 2007; Omi and Shinomoto 2011; Min, Goh, and Vazquez 2011; Karsai et al. 2011; Bagrow and Brockmann 2013). They have been observed in the context of earthquakes, sunspots, neuronal activity and human communication (see Karsai et al. 2012; Vajna, Tóth, and Kertész 2013; Meerschaert and Stoev 2008 for a list of references). Such time series exhibit high activity in some ‘bursty’ intervals, which alternate with other, quiet intervals. Although several mechanisms are plausible explanations for bursty behaviour (most prominently self-exciting point process by Hawkes (1971)), there seems to be one salient feature which very typically indicates the departure from temporal homogeneity: a heavy-tailed distribution of waiting times (Vasquez et al. 2006; Karsai et al. 2012; Vajna, Tóth, and Kertész 2013). As we show below in simulations, a simple renewal process with heavy-tailed waiting times can capture this type of dynamics. For many systems, the renewal property is appropriate; a simple test of the absence of correlations in a succession of waiting times can be undertaken by randomly reshuffling the waiting times (Karsai et al. 2012).

Often a magnitude, or mark can be assigned to each event in the renewal process, such as for earthquakes, solar flares or neuron voltages. The Peaks-Over-Threshold model (POT, see e.g. Coles 2001) applies a threshold to the magnitudes, and fits a Generalized Pareto distribution to the threshold exceedances. A commonly made assumption in POT models is that times between events are either fixed or light-tailed, and this entails that the threshold crossing times form a Poisson process (Hsing, Hüsler, and Leadbetter 1988). Then as one increases the threshold and thus decreases the threshold crossing probability \(p\), the Poisson process is rarefied, i.e. its intensity decreases linearly with \(p\) (see e.g. Beirlant et al. 2006).

As will be shown below, in the heavy-tailed waiting time scenario threshold crossing times form a fractional Poisson process (Laskin 2003; Meerschaert, Nane, and Vellaisamy 2011), which is a renewal process with Mittag-Leffler distributed waiting times. The family of Mittag-Leffler distributions nests the exponential distribution (Haubold, Mathai, and Saxena 2011), and hence the fractional Poisson process generalizes the standard Poisson process. Again as the threshold size increases and the threshold crossing probability \(p\) decreases, the fractional Poisson process is rarefied: The scale parameter of the Mittag-Leffler inter-arrival times of threshold crossing times increases, but superlinearly; see the Theorem below.

Maxima of events which occur according to a renewal process with heavy-tailed waiting times have been studied under the names “Continuous Time Random Maxima process” (CTRM) (Benson, Schumer, and Meerschaert 2007; Meerschaert and Stoev 2008; Hees and Scheffler 2016, 2017), “Max-Renewal process” (Silvestrov 2002; Silvestrov and Teugels 2004; Basrak and Špoljarić 2015), and “Shock process” (Esary and Marshall 1973; Shanthikumar and Sumita 1983, 1984, 1985; Anderson 1987; Gut and Hüsler 1999). The existing literature focuses on probabilistic results surrounding these models. In this work, however, we introduce a method of inference for this type of model, which is seemingly not available in the literature.

We review the marked renewal process in Section 2, and derive a scaling limit theorem for inter-exceedance times in Section 3. We give a statistical procedure to estimate model parameters via stability plots in Section 5, but to set the stage we first discuss inference for the Mittag-Leffler distribution in Section 4. (A simulation study of the effectiveness of our statistical procedure is given in the appendix.) Diagnostic plots for model criticism are discussed in Section 6. In Section 7, we discuss the memory property of the Mittag-Leffler distribution, and how it affects the predictive distribution for the time until the next threshold crossing event. Section 8 concludes. For all statistical computations we have used R (R Core Team 2018). All code and data used for the analysis in this article has been organized into an R package CTRE (https://github.com/UNSW-MATH/CTRE). The source code for the figures generated in this manuscript is available online at https://github.com/UNSW-MATH/bursty-POT.

2 Continuous Time Random Exceedances (CTRE)

As a model for extreme observations, we use a Marked Renewal Process (MRP):

Definition (MRP):

Let \((W,J), (W_1, J_1), (W_2, J_2), \ldots\) be i.i.d. pairs of random variables, where the \(W_k > 0\) are interpreted as the waiting times and \(J_k \in [x_L, x_R]\) as the event magnitudes (\(x_L \in [-\infty, +\infty), x_R \in (-\infty, +\infty]\)). If \(W\) and \(J\) are independent, the Marked Renewal Process is said to be uncoupled.

Note that the \(k\)-th magnitude \(J_k\) occurs at time \(T_k = W_1 + \ldots + W_k\). Based on an MRP, we define the Continuous Time Random Exceedance model (CTRE) as follows:

Definition (CTRE):

Given a threshold \(\ell \in (x_L, x_R)\), consider the stopping time \[\tau(\ell) := \min\{k: J_k > \ell\},\quad \ell \in (x_L, x_R).\] Define the pair of random variables \((X(\ell), T(\ell))\) via \[X(\ell) = J_{\tau(\ell)} - \ell, \quad T(\ell) = \sum_{k=1}^{\tau(\ell)} W_k.\] By restarting the MRP at \(\tau(\ell)\), inductively define the two i.i.d. sequences \(T(\ell,n)\) and \(X(\ell, n)\), \(n \in \mathbb N\), called the “interarrival times” and the “exceedances”, respectively. The pair sequence \((T(\ell, n), W(\ell, n))_{n \in \mathbb N}\) is called a Continuous Time Random Exceedance model (CTRE). If the underlying MRP is uncoupled, then the CTRE is also called uncoupled.

Figure 2.1: Exceedances (red) and Times until Exceedance (durations between blue crosses) for a given threshold \(\ell\) (dashed line).

In this article, we restrict ourselves to the uncoupled case, where the two sequences \(X(\ell, n)_{n \in \mathbb N}\) and \(T(\ell, n)_{n \in \mathbb N}\) are independent.¹ Figure 2.1 shows a simulated dataset in the top panel, where \(W\) has a stable distribution with tail parameter \(\beta =\) 0.8 (and skewness \(1\) and location \(0\)), and where \(J\) is from a standard Gumbel distribution. In the bottom panel, we plot a time series of solar flare intensities derived from a NASA dataset (Dennis et al. 1991)². Clearly, the simulated data exhibit long intervals without any events, whereas in the real-world dataset events appear continuously. The threshold exceedances, however, appear to have similar statistical behaviour in both models. Observations below a threshold are commonly discarded in Extreme Value Theory (POT approach); likewise, the CTRE model interprets these observations as noise and filters them out.

3 Scaling limit of Exceedance Times

In this section we state and prove the key theorem, see below. For an accessible introduction to regular variation and stable limit theorems, we recommend the book by Mark M Meerschaert and Sikorskii (2012).

Theorem:

Let the waiting times \(J_k\) be in the domain of attraction of a positively skewed sum-stable law with stability parameter \(0 < \beta < 1\); more precisely, \[\begin{align} \frac{W_1 + \ldots + W_n}{b(n)} \overset{d}{\longrightarrow} D, \quad n \to \infty \tag{3.1} \end{align}\] for a function \(b(n)\) which is regularly varying at \(\infty\) with parameter \(1/\beta\), and where \(\mathbf E[\exp(-sD)] = \exp(-s^\beta)\). Write \(p := \mathbf P(J > \ell)\). Then the weak convergence \[ \frac{T(\ell)} {b(1/p)} \to W_\beta \quad \text{ as } \quad \ell \uparrow x_R \] holds, where the Mittag-Leffler random variable \(W_\beta\) is defined on the positive real numbers via \[ \mathbf E[\exp(-sW_\beta)] = \frac{1}{1+s^\beta}. \]

For a scale parameter \(\sigma > 0\), we write \({\rm ML}(\beta, \sigma)\) for the distribution of \(\sigma W_\beta\). The Mittag-Leffler distribution with parameter \(\beta \in (0,1]\) is a heavy-tailed positive distribution for \(\beta < 1\), with infinite mean. However, as \(\beta \uparrow 1\), \({\rm ML}(\beta, \sigma)\) converges weakly to the exponential distribution \({\rm Exp}(\sigma)\). This means that although its moments are all infinite, the Mittag-Leffler distribution may (if \(\beta\) is close to 1) be indistinguishable from the exponential distribution, for the purposes of applied statistics. For a detailed reference on the Mittag-Leffler distribution, see e.g. Haubold, Mathai, and Saxena (2011), and for algorithms, see e.g. the R package MittagLeffleR (Gill and Straka 2017).

4 Inference for the Mittag-Leffler distribution

Figure 4.1: Pareto vs. Mittag-Leffler distribution, both with same tail exponent 0.8. Top row: solar flare magnitudes data, fitting nicely with a Pareto distribution. Bottom row: simulated Mittag-Leffler data. Left column: static QQ-Estimator plots, assessing a fit against the Pareto distribution. Centre column: dynamic QQ-Estimator plots. Right column: a usual QQ-Plot based on Mittag-Leffler quantiles.

Since the Mittag-Leffler distribution is heavy-tailed, many researchers would intuitively give the highest importance to the tail behaviour of the distribution, and estimate the exponent of the tail function with established methods such as the Hill estimator. The QQ-estimator (Kratz and Resnick 1996) is closely related to the Hill estimator and fits a least squares line through the logarithms of the ordered statistics (y-axis) and the corresponding quantiles of the exponential distribution (x-axis). The reciprocal slope is returned as the estimate of the tail exponent. For instance, the top-left panel in Figure 4.1 shows the (static) QQ-estimator for the magnitudes of the solar flare data, returning a good fit to a Pareto distribution with tail parameter 0.79. The dynamic QQ-estimator plot (top-center panel) plots the tail exponent estimate for the largest values at cutoffs varying up to the 5th order statistic, and the region of stability at 0.8 indicates a recommended estimate.

However, as described nicely by Resnick (1997), the less similar a heavy-tailed distribution is to a Pareto distribution, the less useful a Hill or QQ-Plot estimator becomes. The bottom-left panel of Figure 4.1, for instance, shows the static QQ-estimator plot for \(2000\) draws from the Mittag-Leffler distribution with tail parameter \(0.8\). The stretched exponential shape means that the Mittag-Leffler distribution has more probability near \(0\) than the Pareto distribution, severely biasing the estimator downwards (0.61). Even by looking at different cutoffs via the dynamic QQ-estimator plot (bottom-center panel) one hardly identifies \(0.8\) as a clear candidate for a tail parameter estimate.

QQ-estimators, and the closely related Hill estimator, are hence not suitable to detect a heavy-tailed Mittag-Leffler distribution. Moreover, QQ-estimator plots of exponentially distributed data are virtually indistinguishable from the bottom-left panel. This shows that Mittag-Leffler distributed data may look like exponentially distributed data if examined via a QQ-estimator.

Hence if there is some prior expectation that the data are drawn from the Mittag-Leffler distribution (as is the the case for threshold exceedance times), then we recommend avoiding QQ-estimators and Hill plots altogether, and instead examining QQ-plots directly, on a logarithmic scale (see Figure 4.1, right column). The scale parameter \(\sigma\) is irrelevant for a QQ-Plot. The tail parameter \(\beta\) can be estimated quickly via the log-moment method by Cahoy (2013), or via maximum likelihood. Both estimators are implemented in MittagLeffleR (Gill and Straka 2017). Since the exponential distribution is nested in the Mittag-Leffler family of distributions, a standard likelihood ratio test can be performed, with the exponential distribution as a null model against a Mittag-Leffler distribution as the alternative. As an example, the threshold crossing times for the solar flare dataset (Figure 6.2, right panel) yield a difference in deviance of \(\approx\) 324, which evaluates to a strongly significant \(\chi^2_1\) test statistic (p-value = \(10^{-72}\)) for the null hypothesis \(\beta = 1\).

5 Inference on Exceedance times

Figure 5.1: Stability plots for the tail and scale parameter of the Mittag-Leffler distribution of the Solar Flare dataset. Dotted horizontal lines are at \(\beta = 0.85\) and \(\sigma_0 = 3 \times 10^7\) seconds \(\approx 0.95\) years.

The Theorem in Section 3 implies that for a high threshold \(\ell\) we may approximate the distribution of \(T(\ell)\) with an \({\rm ML}(\beta, b(1/p))\) distribution, where \(b(\,)\) varies regularly at \(\infty\) with parameter \(1/\beta\). Building on the POT (Peaks Over Threshold) method, we propose the following estimation procedure for the distribution of \(T(\ell)\):

1. For a range of thresholds \(\ell\) near the largest order statistics, extract datasets of exceedance times \(\{T(\ell, i)\}_i\).

2. For each choice of threshold \(\ell\), fit a Mittag-Leffler distribution to the resulting dataset \(\{T(\ell, i)\}_i\). This results in the estimates \(\{\hat\beta(\ell)\}_\ell\) and \(\{\hat \sigma(\ell)\}_\ell\).

3. Plot \(\ell\) vs. \(\hat \beta(\ell)\). As \(\ell\) increases towards \(x_R\), \(\hat \beta(\ell)\) stabilizes around a constant \(\hat \beta\). Use \(\hat \beta\) as an estimate for the tail parameter \(\beta\) of the Mittag-Leffler distribution of exceedance times.

4. Approximate \(p \approx |\{k: J_k > \ell\}| / n\). Recall that \(b(c)\) is regularly varying with parameter \(1/\beta\), and hence has the representation \(b(c) = L(c) c^{1/\beta}\) for some slowly varying function \(L(c)\). Assuming that the variation of \(L(c)\) is negligible, we hence plot \(\ell\) vs. \(p^{1/\hat \beta} \hat \sigma(\ell)\). Again as \(\ell\) increases towards \(x_R\), \(p^{1/\hat \beta} \hat \sigma(\ell)\) is expected to stabilize around a constant \(\hat \sigma_0\). We then use \(p^{-1/\hat \beta} \hat \sigma_0\) as an estimate of the scale parameter of the Mittag-Leffler distribution of exceedance times for the level \(\ell\).

The above approach, though theoretically sound, benefits from the following practical adjustments (compare with Figure 5.1):

• We choose \(\ell\) from the order statistics, i.e. \(\ell\) is the \(k\)-th largest of the observations \(X_j\), where \(k\) runs from \(k_\text{min}, k_\text{min} + 1, \ldots, k_\text{max}\). The datasets are then of length \(k-1\).
• We use \(k\) rather than \(\ell\) for the horizontal axis of our plots.
• In Step 4, rather than plotting \(p^{1/\hat \beta} \hat \sigma(\ell)\) we plot \(k^{1/\hat \beta} \hat \sigma(\ell)\). This changes \(\hat \sigma_0\) by the multiplicative constant \(n^{1/\hat \beta}\), but has the advantage that \(\hat \sigma_0\) does not change if one pre-processes the data by removing all observations below a certain threshold.

The estimates \(\hat \beta\) and \(\hat \sigma_0\) give an estimate of the distribution of exceedance times, dependent on the threshold \(\ell\): \[\begin{align*} T(\ell) \sim {\rm ML}(\hat \beta, k^{-1/\hat \beta} \hat \sigma_0). \end{align*}\] For quick estimates of the Mittag-Leffler parameters we have used the method of log-transformed moments by Cahoy (2013). We have verified the validity of our estimation algorithm via simulations, see the appendix.

6 Checking Model Assumptions

Figure 6.1: Diagnostic plots for the solar flare data: auto-correlation function.

Figure 6.2: Diagnostic plots for the solar flare data: empirical copula and QQ Plot.

The CTRE model is based on three main assumptions, which are repeated below. For each assumption, we suggest one means of checking if it holds:

i.i.d.:

After removing the “noise observations” below the smallest threshold \(\ell_0\), the pair sequence \((T(\ell_0, i), X(\ell_0,i))\) is i.i.d. An indication if this is true is given by an auto-correlation plot for the logarithms (to ensure finite moments) of the two time series.

Uncoupled:

Each \(T(\ell, i)\) is independent of each \(X(\ell, i)\). We propose an empirical copula plot to check for any dependence.

\({\rm ML}(\beta, \sigma)\) distribution of \(T(\ell, i)\):

Apply a cutoff at the lowest threshold \(\ell_0\), extract the threshold crossing times, and create a QQ Plot for the Mittag-Leffler distribution. Use a log-Moment estimate of the tail parameter for the theoretical / population quantiles of the plot.

Figures 6.1 and 6.2 show the diagnostic plots for a minimum threshold chosen at the 200th order statistic. There is some residual autocorrelation for the sequence of threshold exceedance times that is not accounted for by the CTRE model. The fit with a Mittag-Leffler distribution (\(\beta = 0.8\)) is good, though there are signs that the power-law tail tapers off for very large inter-threshold crossing times. There is no apparent dependence between threshold exceedance times and event magnitudes seen in the copula plot.

7 Predicting the time of the next threshold crossing

According to Figure 5.1, for a threshold \(\ell\) at the \(k\)-th order statistic, the fitted threshold exceedance time distribution is \[ T_\ell \sim {\rm ML}(\beta, k^{-1/\beta} \sigma_0), \] where \(\beta = 0.85\) and \(\sigma_0 = 3.0 \times 10^7 {\rm sec}\). Unlike the exponential distribution, the Mittag-Leffler distribution is not memoryless, and the probability density of the time \(t\) until the next threshold crossing will depend on the time \(t_0\) elapsed since the last threshold crossing. This density equals \[ p(t|\beta, \sigma_0, \ell, t_0) = \frac{f(t + t_0 | \beta, k^{-1/\beta} \sigma_0)}{\mathbf P[T_\ell > t_0]} \] where \(f(\,\cdot\, | \beta, k^{1/\beta} \sigma_0)\) is the probability density of \({\rm ML}(\beta, k^{1/\beta} \sigma_0)\). The more time has passed without a threshold crossing, the more the probability distribution shifts towards larger values for the next crossing (see Figure 7.1, left panel). The hazard rate \[ h(t) = \frac{f(t| \beta, k^{1/\beta} \sigma_0))}{\int_t^\infty f(\tau| \beta, k^{-1/\beta} \sigma_0))\,d\tau} \] represents the risk of a threshold crossing per unit time, and is a decreasing function for the Mittag-Leffler distribution. The closer \(\beta\) is to \(1\), the more the hazard rate mimics that of an exponential distribution (a constant function, see Figure 7.1, right panel).

Figure 7.1: Left: Conditional distribution of time until next threshold crossing, depending on elapsed time \(t_0\) since last crossing (\(\beta = 0.8\), \(\sigma_0 = 1\)). Right: Hazard rate depending on tail parameter \(\beta\).

It is beyond the scope of the current paper to incorporate parameter uncertainty into our predictive distribution for the next threshold crossing; however, methods as described by Scarrott and Macdonald (2012) and Lee, Fan, and Sisson (2015) are likely to extend to our setting.

8 Discussion & Conclusion

We have extended the POT (Peaks over Threshold) model, a mainstay of extreme value theory, to “bursty” time series, which have been studied intensively in statistical physics. Burstiness is characterized by power-law waiting times between events, and we have shown that the Mittag-Leffler distribution arises naturally as a scaling limit for the inter-exceedance times of high thresholds. Moreover, we have derived the following non-linear scaling behaviour: \(\sigma \sim p^{-1/\beta}\), where \(\sigma\) is the scale parameter of the distribution of threshold exceedance times, \(p\) is the fraction of magnitudes above the threshold, and \(\beta\) the exponent of the power law. This “anomalous” scaling behaviour in the bursty setting entails two phenomena: i) a heavy tail of the interarrival time distribution of threshold crossings (long rests), and ii) a high propensity for more threshold crossing events immediately after each threshold crossing event (bursts). The Mittag-Leffler distribution captures both phenomena, due to its heavy tail as well as its stretched exponential (peaked) asymptotics for small times. It generalizes the exponential distribution, and in the solar flare data example, this generalization is warranted, because the likelihood-ratio test is strongly significant.

When we introduced the CTRE model, we assumed that all events are i.i.d. This assumption is likely sufficient but not necessary for our limit theorem to hold. Moreover, any data below a (minimum) threshold \(\ell_0\) is discarded for CTREs, and hence need not satisfy the i.i.d. assumption. For the purposes of statistical inference, we merely require that the inter-threshold-crossing times are i.i.d.

The bursty CTRE approach to model “non-Poissonian” threshold crossing times should be contrasted with the (now standard) approach of clusters of extremes, see e.g. Ferro and Segers (2003). In this approach, i.i.d. event sequences of magnitudes are generalized to stationary sequences of event magnitudes (subject to a mixing condition). The two approaches are fundamentally different: A clustering model assumes that each event belongs to one particular (latent) group of events. For bursts, however, the aim is to identify an underlying scale-free pattern in the event dynamics, which is often characteristic of complex systems. It is an interesting open problem to develop quality criteria, based e.g. on measures of surprise (Lee, Fan, and Sisson 2015), which guide an applied statistician in the choice between a clustering and a CTRE approach for a particular problem. Moreover, we believe it may be possible to unify the two approaches by considering CTREs based on MRPs with a stationary, rather than i.i.d., sequence of magnitudes.

Finally, a purely scale-free pattern for event times may be too rigid an assumption for some bursty time series, because often the heavy-tailed character of the inter-arrival time distribution does not hold at all time scales; rather, it applies at short and intermediate time scales, and is truncated (or tempered, reverting to an exponential distribution) at very long time scales (see e.g. Mark M. Meerschaert, Roy, and Shao 2012; and Aban, Meerschaert, and Panorska 2006). In such situations, a “tempered” Mittag-Leffler distribution may provide a more realistic fit, which we aim to introduce in follow-up work.

9 Validating our inference method on simulated data

To test our inference method via stability plots, we have simulated 10000 independent waiting time and magnitude pairs \((W_k, J_k)\) (see upper panel in Figure 2.1. In order to have exact analytical values available for \(\beta\) and \(\sigma_0\), a distribution for \(W_k\) needs to be chosen for which \(b(n)\) from (3.1) is known. If we choose \(W_k \stackrel{d}{=} D\), where \(D\) is as in (3.1), then due to the stability property we have the equality of distribution \(W_1 + \ldots + W_n \stackrel{d}{=} b(n) D\), for \(b(n) = n^{1/\beta}\). Using the parametrisation of Samorodnitsky and Taqqu (1994), a few lines of calculation (see e.g. the vignette on parametrisation in Gill and Straka (2017)) show that \(D\) must have the stable distribution \(S_\beta(\cos(\pi \beta/2)^{1/\beta}, +1, 0)\), which is implemented in the R package stabledist by Wuertz, Maechler, and members. (2016).

Figure 9.1: Tail and scale estimates for simulated data, with waiting times drawn from the stable distribution \(S_\beta(\cos(\pi \beta/2)^{1/\beta}, +1, 0)\) with \(\beta = 0.8\). Dashed lines are 95% confidence intervals, dotted lines are the known theoretical values (\(0.8\) and \(10000^{1/0.8}\)).

By the Theorem, the distribution of \(T(\ell)\) is approximately \[ {\rm ML}(\beta, p^{-1/\beta}) = {\rm ML}(\beta, k^{-1/\beta} n^{1/\beta}), \] which means \(\sigma_0 = n^{1/\beta}\). The distribution of \(J_k\) is irrelevant for the inference on \(\beta\) and \(\sigma_0\) (we have chosen unit exponential random variables). Figure 9.1 displays plots of \(\hat \beta(\ell)\) and \(\hat \sigma(\ell)\) vs. \(k\); recall that \(k\) is the index of the order statistics of \(J_k\) at which the threshold \(\ell\) is placed. Dotted lines show 95% confidence intervals, which are derived from the asymptotic normality of the log-moments estimators (Cahoy 2013) and the \(\delta\)-method (Gill and Straka 2017). The dashed lines show the actual values of \(\beta\) resp. \(\sigma_0\), showing that our inference method identifies the parameters correctly.