Abstract
In the statistical physics literature, heavy-tailed inter-arrival times are said to be the main signature of b burstyb dynamics. Such dynamics have been observed for financial time series, earthquakes, solar flares and neuron voltage spikes. This talk is concerned with the modelling of extremes of events when events follow such bursty dynamics. We assign i.i.d. magnitudes (marks) to the events in a heavy-tailed renewal process and apply the b Peaks Over Thresholdb method from Extreme Value Theory. Leveraging geometric stability, it can be shown that the threshold exceedance times asymptotically form a b fractional Poisson processb with Mittag-Leffler inter-arrival times. Finally, we discuss methods for inference on the tail and scale parameters of the bursty dynamics.
Earthquakes, Neurons, Phone calls
Karsai, M., Kaski, K., BarabC!si, A. L., & KertC)sz, J. (2012). Universal features of correlated bursty behaviour. Scientific Reports, 2. http://doi.org/10.1038/srep00397
Heavy-tailed renewal processes model intermittency
Power law inter-arrival times in the literature
- Resnick and Starica (1995): quiet periods between transmissions for a networked computer terminal, \(\beta \approx 0.6\)
- Mainardi et al. (2000): waiting times between trades of bond futures, \(\beta \approx 0.95\)
- BarabC!si (2005): \(\beta \approx 1\) for waiting time between emails, and for waiting time for a reply
- Berhondo et al. (2006): \(1 < \beta < 2\), for waiting times between coronal mass ejections (solar flares)
- Brockmann et al. (2006): \(\beta \approx 1.05\), waiting times for paper currency to travel beyond the limits of a metropolitan area. They also report a power law distribution in the ensuing travel distance
- Cheng et al. (1995): \(\beta \approx 1.6\) in the waiting time between earthquakes, and a similar value for waiting time between b starquakesb in neutron stars
- Lavergnat and GolC) (1998): \(\beta \approx 0.68\), waiting time between large raindrops
- Smethurst and Williams (2001): \(\beta \approx 1.4\), waiting times for a doctor appointment
- Sabatelli et al. (2002): \(\beta \approx 0.4\), for 19th century Irish stock prices, and \(\beta \approx 1.87\) for the modern Japanese Yen currencymarket.
Testing the renewal assumption: ACF
But: for renewal processes with \(1 < \beta < 2\) (finite mean, infinite variance) the ACF decays with exponent \(2-\beta\).
Quick fix: log-transform the waiting times.
Testing the renewal assumption: event trains
- Train of events given a time window \(\Delta t\):
a sequence of consecutive events where all inter-arrival times are \(< \Delta t\).
- Plot the distribution of train lengths \(E(\Delta t)\) for varying \(\Delta t\),
- Reshuffle the inter-event times and plot again (b Null-distributionb )
If the two plots are visually different, renewal assumption is false.
(Karsai et al 2012, Universal features of correlated bursty behaviour, Scientific Reports)
Renewal assumption b OKb
Renewal assumption b not OKb
Waiting time b Magnitude pairs
Clustering coefficients for a growing network, where nodes are preferentially attached to nodes with large clustering coefficients.
Bagrow, J. P., & Brockmann, D. (2013). Natural emergence of clusters and bursts in network evolution. Physical Review X, 3(2), 1b 6. http://doi.org/10.1103/PhysRevX.3.021016
Waiting time b Magnitude pairs
| Server requests |
filesize |
| Earthquakes |
magnitudes |
| Email sent |
email replied |
| Solar flares |
magnitudes |
| Time between rest for paper currency |
distance travelled |
| Time between raindrops |
raindrop size |
| Time between price jumps |
price jumps |
Research Question
How does one model magnitude and timing of extreme events, if the underlying dynamics are a heavy-tailed renewal process?
Example: Earthquakes in Coral Sea
CTRM-app
3 Inference via POT method (Peaks over Threshold)
Exceedance
Define the Exceedance as
\[X(\ell) = M(T_\ell) - \ell.\]
In the uncoupled case, it follows easily that Exceedance \(X(\ell)\) and Exceedance Time \(T(\ell)\) are independent.
- Theorem (POT, e.g.B Leadbetter et al.):
Asymptotically for large \(\ell\), \(X_\ell\) follows a generalized Pareto distribution:
\[X_\ell \sim GP(\xi, \bar \sigma)\]
Extract the red lines and fit GEV distributions to them!
POT method for magnitudes
For threshold \(\ell\) in \(\{J_{(1)} > J_{(2)} > \ldots > J_{(n)}\}\):
- extract number of exceedances \(k\).
- extract exceedances \(X(\ell)\)
- fit \({\rm GP}\) distribution to the \(X(\ell)\)
Then plot parameters estimates \(\hat \xi, \hat \sigma\) against \(k\).
Example: earthquake magnitudes shape parameterb &
b & and scale parameter: 
By inspection, pick \(\hat \xi = -0.1\), \(\hat \sigma = 0.6\).
- Bias-variance trade-off:
Small \(k\): few data (high variance) but asymptotically close to GP distribution (low bias)
High \(k\): many data (low variance) but asymptotically far away from GP distribution (high bias)
POT method for Exceedance Times
For threshold \(\ell\) in \(\{J_{(1)} > J_{(2)} > \ldots > J_{(n)}\}\):
- extract number of exceedances \(k\).
- extract exceedance times \(T(\ell)\)
- fit a distribution to the \(T(\ell)\)
Then plot the parameters of this distribution against \(k\).
But what is the distribution of \(T(\ell)\)??
\(T(\ell) =\) distances between blue circles.
Mittag-Leffler Function
- Definition:
\[E_\beta(z) = \sum_{k=0}^\infty \frac{z^k}{\Gamma(1+\beta k)}\]
where \(0 < \beta \le 1\).
\(E_\beta(z)\) nests the exponential function: \(E_1(z) = e^z\).
For \(\beta < 1\), \(E_\beta(-z^\beta)\) is \(\rm{RV}_\infty(-\beta)\).
Mittag-Leffler-1 Distribution
- Definition:
A r.v. \(T_\beta\) is Mittag-Leffler-1 distributed with tail parameter \(\beta \in (0,1]\) if
\[\mathbf P(T > t) = E_\beta(-t^\beta), \quad t > 0.\]
- Laplace transform: \[\int_0^\infty e^{-st} \mathbf P(T \in dt) =
\frac{1}{1+s^\beta}\]
- geometrically stable
- very heavy-tailed (infinite mean), but: \(T \stackrel{d}{\to} {\rm Exp(1)}\) as \(\beta \uparrow 1\)
- Introduce the scale parameter \(\sigma > 0\) via \(\sigma T \sim \rm{ML}(\beta,\sigma)\)
Mittag-Leffler-2 Distribution
- Definition:
A r.v. \(T_\beta\) is Mittag-Leffler-2 distributed with parameter \(\beta \in (0,1]\) if
\[\int_0^\infty e^{-st} \mathbf P(T \in dt) = E_\beta(-s).\]
We write \(\sigma T_\beta \sim {\rm ML2}(\beta, \sigma)\) for \(\sigma > 0\).
- Example:
Let \(D(u)\) be a stable subordinator with Laplace Transform \(\exp(-us^\beta)\).
Let \(E(t) = \inf\{u: D(u) > t\}\) (inverse stable subordinator).
Then
\[E(t) \sim {\rm ML2}(\beta, t^\beta).\]
Geometric stable distributions
- Definition:
A random variable \(Y\) is geometric stable if and only if its characteristic function is
\[\psi(t) = \frac{1}{1 + \sigma^\alpha |t|^\alpha
\omega_{\alpha, \beta}(t) - i\mu t}\]
where for \(\alpha \neq 1\)
\[\omega_{\alpha, \beta}(t)=
\begin{cases}
1 - i\beta {\rm sign}(t) \tan(\pi \alpha / 2), & \alpha \neq 1
\\
1 + i\beta(2/\pi) {\rm sign}(t) \log|t|, & \alpha = 1
\end{cases}\]
Write \(Y \sim {\rm GS}(\alpha, \beta, \sigma, \mu)\); (stability, skewness, scale, location).
(Parametrization by Samorodnitsky & Taqqu, see Kozubowski & Rachev, 1999)
- Proposition:
A r.v. \(Y\) on \(\mathbb R\) is GS (geometric stable) if for \(N_p \sim {\rm Geom}(p)\) there exist i.i.d. r.v.b s \(X_1, X_2, \ldots\), independent of \(N_p\), a norming sequence \(A_p\) and centering constants \(b_p\) such that
\[A_p \sum_{k=1}^{N_p} (X_k + b_p) \stackrel{d}{\to} Y \text{ as } p \downarrow 0.\]
Then \(X_1\) is weakly geometrically attracted to \(Y\), and the collection of such \(X_k\) is called the generalized domain of geometric attraction (GDOGA) of \(Y\).
- Proposition (Gnedenko 1972):
Let \(\beta \in (0,1]\). Then \(X \in \text{GDOGA}({\rm ML(\beta,1)})\) if and only if
\[1 - F_X(t) \in {\rm RV}_\infty(-\beta)\]
Exceedance Time = geometric sum
Recall:
- \(M(t)\) is CTRM process, and
- \(T(\ell) = \inf\{t: M(t) > \ell\}\) is exceedance time.
Moreover, \(N_{p_\ell} := \min\{k: J_k > \ell\}\) is geometric, independent of \(W_k\), and
\[T(\ell) = \sum_{k=1}^{N_{p_\ell}} W_k.\]
Scaling limit theorem for Exceedance Time
- Theorem:
Suppose the magnitudes \(J_k\) are supported on the interval \([x_L, x_R]\).
Suppose that \(\mathbf P(W_k > t) = L(t) t^{-\beta}\) for some slowly varying function \(L(t)\) and \(\beta \in (0,1)\).
Let \(\ell \in [x_L, x_R]\) be a threshold, and write \(p_\ell = 1 - F_J(\ell)\) for the probability that a magnitude exceeds threshold \(\ell\).
Then for some \(\sigma > 0\),
\[p_\ell^{1/\beta} T(\ell) \stackrel{d}{\to} {\rm ML}(\beta, \sigma),
\quad \ell \uparrow x_R\]
Side note
Gnedenko (1972) predates the Shock-process literature (1988+); Yet geometric stability does not appear in the Shock process literature.
It took the end of the Cold War and the invention of the internet to put the two things together!
Consequence
For large thresholds \(\ell\), we may approximate
\[T(\ell) \approx {\rm ML}(\beta, \sigma p_\ell^{-1/\beta}).\]
- Bias-variance trade-off:
Small k: few data (high variance) but many geometric summands, hence asymptotically close to M-L distribution (low bias)
High k: many data (low variance) but few geometric summands, hence far away from GP distribution (high bias)
Two-stage estimator for \(\beta\) and \(\sigma\)
CTRM-app
- Let \(p_\ell \approx |\{k: J_k > \ell\}| / n.\)
- For various thresholds \(\ell \in [x_L, x_R]\), fit an ML1 distribution to \(\{T_i(\ell)\}\) and get \[\{\hat \beta(\ell)\}_\ell, \{\hat \sigma(\ell)\}_\ell\]
- Plot \(\hat \beta(\ell)\) vs. \(\ell\), and pick \(\hat \beta\) where stable.
- Plot \(p_\ell^{1/ \hat \beta} \hat \sigma(\ell)\), and pick \(\hat \sigma\) where stable.
Stability plot for M-L tail parameter
Stability plot for M-L scale parameter
Mittag-Leffler estimation via log-moments
Let \(T \sim ML(\beta, \sigma)\).
- The Mittag-Leffler distribution has no moments if \(\beta \in (0,1)\); but the first two centered moments of its logarithmic transform are finite: \(\newcommand{\ex}{\mathbf E}\) \[\mu_{\log T} = \ex[\log T] = \log(\sigma) - \gamma,
\quad \sigma_{\log T} = Var[\log T] = \frac{\pi^2}{6} \left( \frac{2}{\beta^2} - 1 \right)\]
- Find parameters \(\hat \beta\) and \(\hat \sigma_T\) for which the first two centered moments match their empirical counterparts \[\hat \mu_{\log T} = \frac{1}{n} \sum_{j=1}^n \log(T_j), \quad
\hat \sigma_{\log T} = \frac{1}{n} \sum_{j=1}^n (\log(T_j) - \hat \mu_{\log T})^2\]
(Cahoy et al b 10)
ML-estimation via maximum-likelihood
Recent algorithm (Garrappa, 2015) calculates Mittag-Leffler Function by inversion of Laplace transform.
The MittagLeffleR package (Gill & Straka, 2017) implements this algorithm and calculates pdf, cdf and random variates for both types of Mittag-Leffler distributions.
But this is more costly than method of log-moments.
The case 1 < beta < 2
Suppose now that \(\mathbf P(W_k > t) = L(t) t^{-\beta}\), where \(1 < \beta < 2\) and \(L(t)\) is a slowly varying function.
Convergence to b homogeneousb dynamics takes very long:
\[\mathbf E[N(t)] \approx \frac{t}{\mathbf E[W_1]} + \frac{C}{\mathbf E[W_1]^3} t^{2-\beta}, \quad \text{large } t\]
Exceedance times
If the tail parameter satisfies \(1<\beta<2\), then \(\mathbf E[W_1] < \infty\), and in the limit as \(\ell \uparrow x_R\),
\[p_\ell T(\ell) = p_\ell \sum_{k=1}^{N_{p_\ell}} W_k
\to {\rm Exp}(\mu_W).\]
Hence we b couldb approximate
\[T(\ell) \sim {\rm Exp}(\mu_W/p_\ell)\]
But we expect this convergence to be very slow, as seen in previous slide!
Second order the approximation
\[
\begin{multline}
\sum_{k=1}^{N_{p_\ell}} W_k
= \sum_{k=1}^{N_{p_\ell}} (W_k - \mu_W) + N_{p_\ell} \mu_W
\\
= p_\ell^{-1/\beta} p_\ell^{1/\beta} \sum_{k=1}^{N_{p_\ell}} (W_k - \mu_W)
+ p_\ell^{-1}p_\ell N_{p_\ell} \mu_W
\approx p_\ell^{-1/\beta} \sigma Z^{1/\beta} X + p_\ell^{-1} \mu_W Z
\end{multline}
\]
where \(Z \sim {\rm Exp}(1)\) and \(X \sim S(\beta, 1, 1, 0)\). The latter r.v. is geometric stable (Kozubowski & Rachev, 1999), and so:
\[T(\ell) = \sum_{k=1}^{N_{p_\ell}} W_k
\approx {\rm GS}(\beta, 1, \sigma p_\ell^{-1/\beta}, p_\ell^{-1} \mu_W) \text{ for large } \ell. \]
Proposed estimator
- For varying thresholds \(\ell \in [x_L, x_R]\), extract i.i.d. \(\{T_i(\ell)\}\).
- For every \(\ell\), fit a \({\rm GS}(\beta, 1, \sigma, \mu)\). This yields \(\{\hat \beta(\ell)\}_\ell\), \(\{\hat \sigma(\ell)\}_\ell\) and \(\{\hat \mu(\ell)\}_\ell\).
- Plot \(\ell\) vs. [\(\hat \beta(\ell)\) and \(p_\ell \hat \mu(\ell)\)], and choose \(\hat \beta\) and \(\hat \mu\) where stable
- Plot \(\ell\) vs. \(p^{1/\beta} \hat \sigma(\ell)\) and choose \(\hat \sigma\) where stable.
Available estimators for \({\rm GS}\) distributions:
- fitting the characteristic function (Kozubowski 1999)
- ABC (approximate Bayesian computation, work in progress).
