Runs of extremes of observables on dynamical systems and applications.

The impact of successive extreme weather events on populations has become a significant topic of discussion in climate literature. Recent notable examples include the 2019 global heat-waves, which had severe impacts in Australia, Europe, and the United States, and the 2022 heat-wave in western Europe resulting in a nationally declared emergency in the United Kingdom. Other examples include the 2021-22 floods in Eastern Australia, the 2021 flood in the United Kingdom and Western Europe and the 2017 flood in Texas, all attributed to prolonged heavy rainfall. Understanding the returns of successive extreme weather events like these is crucial for preparing and mitigating disastrous effects. For instance, fire retardant materials can be used to prevent wildfires during heat-waves, and earlier evacuation of populations from high-risk flood zones can be initiated prior to episodes of prolonged rainfall.

Accurately predicting the likelihood of weather anomalies, such as heat-waves, cold spells, or prolonged periods of rainfall, is a challenging problem. In this paper, we explore the use of extreme value theory (EVT) to estimate the probability of successive exceedances of a threshold of a time-series of observations on a variety of uniformly hyperbolic dynamical systems. The study of chaotic dynamical systems provides insights into the statistical behavior of complex systems, such as climate models. We test some of the predictions we derive from uniformly hyperbolic systems on climate data and climate models.

Provided the time-series of an observable $\phi$, namely $(\phi\circ T^{j})_{j\geq 0}$, follow the assumptions outlined in [32], an extreme value law holds for the maxima of the time-series, given by $M_{j}=\max\{\phi\dots,(\phi\circ T^{j})\}$, which guarantees a limiting probability distribution for the sequence of maxima $(M_{j})$. The traditional measure of successive extremes is the extremal index, a parameter which appears in the limiting distribution of the extreme value law and measures the expected number of over-threshold exceedances observed in a short block of the time-series.

For appropriate scaling choices, the exact limiting distribution is the generalized extreme value distribution (GEV), $G(\xi,\mu,\sigma)$, with shape $\xi$, location $\mu$ and scale $\sigma$ parameters. In practice, the block maximum approach is commonly employed which involves dividing the time-series into blocks of a fixed length and performing some likelihood fitting of the parameters of the GEV using the sequence of block maxima. By nature of the block maximum approach, the extremal index is normalized to 1 so that all information on the successive properties of extremes of the data becomes lost. However, a careful choice of observable can preserve some of the properties of successive extremes in this setting.

In climate applications, successive extremes are often thought of in two ways: a high average in a short run, for example, a high daily temperature average over 7 days; or successive over-threshold values occurring in a short run, for example, a high daily temperature observed every day for 7 days. Motivated by [21] and these definitions, we investigate the time-averages of observables over a time-window of integer length $k$, $Y(x)=\frac{1}{k}\sum_{j=0}^{k-1}\phi(T^{j}x)$, as well the exceedance function $Y(x)=\min\{\phi(x),\ldots,\phi(T^{k-1}x)\}$.

This paper’s focus is on whether the scaling constants in the time-series of $k$-exceedances or time-averages can be derived from that of the original time-series. We establish scaling constants in some interesting scenarios for hyperbolic systems. Specifically, we consider an observable $\phi:X\rightarrow\mathbb{R}$ on an ergodic dynamical system $(T,X,\mu)$, which is maximized on an invariant repelling set $S$. The prolonged bouts of extreme weather events are associated with, for example, a weather system being in a quasi-invariant state of extreme rainfall or temperature. This is broadly modeled by the time-series of an observable maximized on an invariant set in a chaotic dynamical system.

We begin our investigation by establishing a general lemma relating the parameters of the GEV for the $k$-exceedances and time-averages coming from any system, with a Fréchet or Weibull limiting distribution, satisfying the conditions of [32]. This general result is then used to establish more concrete theorems for climate-relevant observables taken on hyperbolic systems. We end with a detailed illustration of the applications of these lemmas and theorems to numerically model successive extremes of time-series data coming from: certain hyperbolic dynamical systems; temperature data simulated from the general circulation model, PUMA; and real-world temperature and precipitation data taken from weather stations throughout Germany. Finally, we illustrate that our scaling results can be applied to obtain more accurate numerical estimates of the GEV parameters for maxima of $k$-exceedance and time-averaged functionals of increasing $k$ window lengths than those of traditional maximum likelihood estimation.

This research was motivated in part by the interesting paper [21] in which the authors analyzed the time and spatial averages of daily temperature anomalies over a large number of days (heat-waves or cold spells) via a large deviation approach. Large deviation theory was used to estimate probabilities associated to heat-wave or cold-spell occurrences. The averaging period ranged roughly from 10 to 40 days. The main testing ground was the PUMA climate model of mid-atlantic latitudes, 30 degrees to 60 degrees. The important question they addressed is how to estimate the probability of very rare events, such as long runs of extreme weather, from a limited time-series of recordings. Their predictions were favorably compared to those estimated by standard extreme value theory (Peaks over threshold model leading to a generalized Pareto distribution (GPD)) applied to long simulations of the PUMA model. A related approach based on a large deviation algorithm combined with an importance sampling technique was given in [38].

Our dynamical models and numerical results suggest that in certain settings we do not need to average over long time windows to ensure the applicability of extreme value theory and in some sense we may estimate rare prolonged anomalies from the original time-series of data. The advantage of our method is that while data on, for example, prolonged spells of rainfall is sparse, the simple scalings we find enable us to use all available daily maxima data to estimate the probabilities of prolonged rain of a certain duration.

Suppose $\phi:X\to\mathbb{R}$ is an observable on an ergodic dynamical system $(T,X,\mu)$. Extreme value theory (EVT) considers distributional limits of the maxima process

under suitable scaling constants $a_{n}(M_{n}-b_{n})$.

P. Collet [3] introduced techniques from EVT to establish hitting and return time statistics to ‘generic’ non-recurrent points in some one-dimensional non-uniformly hyperbolic dynamical systems. He considered the observable $\phi(x)=-\ln d(x,x_{0})$, where $x_{0}$ is non-recurrent, in particular not periodic. In this setting EVT is related to hitting and return time statistics since $M_{n}$ increases as orbits make closer approaches to $x_{0}$. Following Collet’s work there have been many papers developing EVT for hyperbolic dynamical systems. The theory is well-developed for $\phi$ maximized at a point $x_{0}$ and a function of distance, as in Collet’s case of $\phi(x)=-\ln d(x,x_{0})$. The paper [2] extended some of this work to the more general setting of observables maximized on submanifolds in phase space, while  [15] considered observables maximized on Cantor sets. In a series of influential papers Freitas et al [16, 18] adapted techniques from [32] to the context of deterministic dynamical systems and elucidated the close relation between extremes, hitting time statistics and Poisson processes. For these and other developments we refer to the book [33] for more details.

More generally suppose $(X_{n})$ is a stationary process with probability distribution function $F_{X}(u):=\mu(X\leq u).$ When determining the distributional limit of

it is well-known that to find the appropriate scaling, given $\tau\in\mathbb{R}$, one should define $u_{n}(\tau)$ as a sequence satisfying $n\mu(X_{0}>u_{n}(\tau))\to\tau$, as $n\to\infty$. We say that $(X_{n})$ satisfies an extreme value law if

for some $\theta\in(0,1]$. The number $\theta$ is called the extremal index and $\frac{1}{\theta}$ roughly measures the average cluster size of exceedances given that one exceedance has occurred. When $(X_{n})$ is iid and has a regularly varying tail it can be shown that under the scaling $u_{n}$ there exists an extreme value law (EVL) and $\theta=1$. If $(X_{n})$ is only assumed stationary rather than iid then the existence of an EVL with $\theta=1$ has been shown provided two conditions hold: (1) $D(u_{n})$ (a mixing condition) and; (2) $D^{\prime}(u_{n})$ (a non-recurrence condition). Thus in the case of an observable on a dynamical system if the time series of observations $X_{n}=\phi\circ T^{n}$ satisfy $D(u_{n})$ and $D^{\prime}(u_{n})$ (or some variation thereof) then an EVL holds with $\theta=1$.

We will use two conditions $\operatorname{\textrm{Д}}(u_{n})$ and $\operatorname{\textrm{Д}}^{{}^{\prime}}_{q}(u_{n})$, adapted to the dynamical setting, introduced in the work [18] that imply a non-degenerate EVL in the case the EI $\theta\not=1$ and also give the value of $\theta$.

We let $(T,X,\mu)$ be an ergodic dynamical system. Here $T$ is a uniformly hyperbolic measure-preserving map of a Riemannian metric space $(X,\mu)$ equipped with a probability measure $\mu$ which is equivalent to Lebesgue measure. We let $\phi:X\to\mathbb{R}$ be an observable which is a function of Euclidean distance and maximized on a set $S$. In this section we consider the extreme value theory of $k$-exceedances and time-averages of the observable. More precisely we consider the EVT of the derived time-series $(Y\circ T^{i})$ where $Y=\min\{\phi,\phi\circ T,\ldots,\phi\circ T^{k-1}\}$ in the case of $k$-exceedances or $Y=\frac{1}{k}\sum_{j=0}^{k-1}\phi\circ T^{j}$ for time-averages.

In this paper $S$ is taken to be either a generic non-recurrent point $x_{0}\in X$; a fixed or periodic point $x_{0}$; or a line or curve, invariant or not. The case of a periodic orbit of minimal period $q$ reduces to the case of a fixed point for $T^{q}$, so we only discuss the case of fixed points. We will consider three classes of uniformly hyperbolic dynamical system, namely: uniformly expanding maps of the interval; hyperbolic toral automorphisms; and coupled expanding maps. We aim not for the greatest generality of dynamical system, as the statements and proofs soon become tiresomely technical. We wish rather for simple illustrative models from chaotic dynamics which hopefully shed light on general principles and the behavior of complex physical models. The extreme value theory for observables $\phi:X\to\mathbb{R}$ maximized on such sets $S$ for these systems, among others, was investigated in [2]. A motivation of [2] was to extend the theory of EVT for dynamical systems beyond that developed for observables maximized at points to more relevant observables to applications. One goal of this paper is to extend this investigation to observables such as $k$-exceedances and time-averages. Another goal is to allow predictions of the GEV for $k$-exceedances and time averages from the GEV of the original time series, both in dynamical models and climate data.

We now compare the GEV scaling constants for an observable $\phi$ for the maximal process $M_{n}=\max\{\phi,\phi\circ T,\ldots,\phi\circ T^{n-1}\}$ with those of the functional $Y=\min\{\phi,\phi\circ T,\ldots,\phi\circ T^{k-1}\}$ or $Y=\frac{1}{k}\sum_{j=0}^{k-1}\phi\circ T^{j}$ with corresponding maximal process $B_{n}=\max_{1\leq j\leq n-1}Y_{j}$.

We first consider the scaling in the case of a Fréchet distribution, a setting in which the scaling constants $b_{n}$ in a linear scaling are zero and then the Weibull. The constant $g(k,T)$ in the lemma below is a constant depending upon $k$ and the map$T$. One of the goals of this paper is to calculate $g(k,T)$ in some simple examples.

Proof:

Assume $\phi$ has a Fréchet distribution, then $b_{n}=0$. To obtain linear scaling for the sequence $u_{n}(\tau)$ we write $t=\rho(\tau)=\tau^{-\xi_{1}}$ where $\xi_{1}>0$ is the shape parameter for the limiting distribution of $M_{n}$. Thus $u_{n}(\tau)=a_{n}t$ and $w_{n}(\tau)=g(k,T)a_{n}\rho(\tau)$. For simplicity of exposition we define $\alpha_{n}=g(k,T)a_{n}$. We obtain

So for fixed large $n$,

and

Note that if $H(t)$ is a generalized EVT distribution with parameters $(\mu,\sigma,\xi)$ then $H(\gamma t)$, $\gamma>0$, has the same shape parameter $\xi_{\gamma}=\xi$ while $\sigma_{\gamma}=\frac{\sigma}{\gamma}$, $\mu_{\gamma}=\frac{\mu}{\gamma}$.

Thus if $G^{\theta_{1}}(t/a_{n})$ has parameters $\xi_{1}$, $\mu_{1}$ and $\sigma_{1}$ then $G^{\theta_{1}}((t/a_{n})\cdot(\frac{a_{n}}{\alpha_{n}}))$ has GEV parameters $\xi^{{}^{\prime}}_{2}=\xi_{1}$, $\sigma^{{}^{\prime}}_{2}=g(k,T)\sigma_{1}$ and $\mu_{2}^{{}^{\prime}}=g(k,T)\mu_{1}$.

Furthermore it is known [4, Theorem 5.2, Page 96] that if $G(t)$ has GEV parameters $\xi$, $\mu$ and $\sigma$ then $G^{\theta}$ has GEV parameters $\xi_{\theta}=\xi$, $\mu_{\theta}=\mu-\frac{\sigma}{\xi}(1-\theta^{\xi})$ and $\sigma_{\theta}=\sigma\theta^{\xi}$.

Thus $[G^{\theta_{1}}((t/a_{n})\cdot(\frac{a_{n}}{\alpha_{n}})]^{\frac{\theta_{2}}{\theta_{1}}}]$ has GEV parameters $\xi_{2}=\xi_{1}$, $\sigma_{2}=g(k,T)\sigma_{1}(\frac{\theta_{2}}{\theta_{1}})^{\xi_{1}}$ and $\mu_{2}=g(k,T)\mu_{1}-\frac{g(k,T)\sigma_{1}}{\xi_{1}}(1-(\theta_{2}/\theta_{1})^{\xi_{1}})$.

Comparing the two forms of $G^{\theta_{1}}(\frac{t}{a_{n}})$, namely $e^{-\theta_{1}(\frac{t}{a_{n}})^{-\frac{1}{\xi_{1}}}}$ and $\exp\bigg{\{}-\theta_{1}\bigg{[}1+\xi_{1}\bigg{(}\frac{t-\mu_{1}}{\sigma_{1}}\bigg{)}\bigg{]}^{\frac{-1}{\xi_{1}}}\bigg{\}}$, we see that the relation $1=\frac{\xi_{1}\mu_{1}}{\sigma_{1}}$ holds. Hence,

This concludes the proof in the case of a Fréchet distribution.

Now suppose that $\phi$ has a Weibull distribution, a setting in which the scaling constants $b_{n}$ may be taken as $C$, the supremum of the range of $\phi$. To obtain a linear scaling for the sequence $u_{n}(\tau)$ we write $t=\rho(\tau)=-\tau^{-\xi_{1}}$ where $\xi_{1}<0$ is the shape parameter for the limiting distribution of $M_{n}$. Thus $u_{n}(\tau)=a_{n}t+C$ and $w_{n}(\tau)=\alpha_{n}t+C$ where $\alpha_{n}=g(k,T)a_{n}$ by assumption. We obtain

We now put the distribution into the standard form of the GEV. So for fixed large $n$, let $z=a_{n}t+C$, so that

and similarly

In $G^{\theta_{1}}$ we have

which gives the relations $1-\frac{\xi_{1}\mu_{1}}{\sigma_{1}}=\frac{C}{a_{n}}$ and $\frac{\xi_{1}}{\sigma_{1}}=\frac{-1}{a_{n}}$. Similarly if $\alpha_{n}=g(k,T)a_{n}$ replaces $a_{n}$ in $G^{\theta_{1}}$ we have, as the shape parameter $\xi_{1}$ does not change, $1-\frac{\xi_{1}\mu_{2}^{{}^{\prime}}}{\sigma_{2}^{{}^{\prime}}}=\frac{C}{g(k,T)a_{n}}$ and $\frac{\xi_{1}}{\sigma_{2}^{{}^{\prime}}}=\frac{-1}{g(k,T)a_{n}}$. We obtain $\sigma_{2}^{{}^{\prime}}=g(k,T)\sigma_{1}$ and $\mu_{2}^{{}^{\prime}}=\frac{\sigma_{1}}{\xi_{1}}(g(k,T)-1)+\mu_{1}$.

Now we account for the transformation $G^{\theta_{1}}\to G^{\theta_{2}}$ by considering the parameters in $[G^{\theta_{1}}]^{\frac{\theta_{2}}{\theta_{1}}}$. As in the Fréchet case we use the fact that if $G(t)$ has GEV parameters $\xi$, $\mu$ and $\sigma$ then $G^{\theta}$ has GEV parameters $\xi_{\theta}=\xi$, $\mu_{\theta}=\mu-\frac{\sigma}{\xi}(1-\theta^{\xi})$ and $\sigma_{\theta}=\sigma\theta^{\xi}$.

Hence the parameters for $B_{n}$ are $\sigma_{2}=g(k,T)\sigma_{1}(\frac{\theta_{2}}{\theta_{1}})^{\xi_{1}}$ and $\mu_{2}=\frac{\sigma_{1}}{\xi_{1}}(g(k,T)-1)+\mu_{1}-\frac{g(k,T)\sigma_{1}}{\xi_{1}}(1-(\frac{\theta_{2}}{\theta_{1}})^{\xi_{1}})$.

We assume that $S$ is a repelling fixed point $x_{0}$ in the case of (A) or (B) or repelling invariant hyperplane of synchrony in the case of $(C)$. Let

where $\phi(x)=d(x,S)^{-\alpha}$ ($\alpha>0$) in the Fréchet case or $\phi(x)=C-d(x,S)^{-\alpha}$ ($\alpha<0$) in the Weibull case . Recall $M_{n}(x):=\max\{\phi(x),\phi(Tx),\ldots,\phi(T^{n-1}x)\}$ and $B_{n}=\max\{Y(x),\ldots,Y(T^{n-1}(x))\}$.

We first consider the Fréchet case. Consider $Y(x)=\min\{\phi(x),\ldots,\phi(T^{k-1}x)\}$, where $\phi(x)=d(x,x_{0})^{-\alpha}$ and $\alpha>0$. Suppose further that $x_{0}$ is non-recurrent. Then there exists a $\delta>0$ such that at least one of the iterates $T^{j}x$, $j=0,\ldots,k-1$ satisfies $d(T^{j}x,x_{0})>\delta$ and hence $Y(T^{n}x)$ is uniformly bounded. Thus the process $Y(T^{n}(x))$ is in the domain of attraction of a Weibull distribution rather than a Fréchet distribution. In general the form of the Weibull distribution cannot be discerned readily from $\phi$, except for some special cases. We illustrate with an example.