On the computation of the extremal index for time series

TL;DR

This paper reviews the extremal index in dynamical systems, introduces a new estimator, and explores its applications to noisy systems and climate data, enhancing understanding of local and global system properties.

Contribution

It provides a comprehensive review of the extremal index in dynamical systems and introduces a novel estimator linked to existing statistical methods.

Findings

New estimator of the extremal index introduced

Relation established between the new and existing estimators

Applications demonstrated on climate data

Abstract

The extremal index is a quantity introduced in extreme value theory to measure the presence of clusters of exceedances. In the dynamical systems framework, it provides important information about the dynamics of the underlying systems. In this paper we provide a review of the meaning of the extremal index in dynamical systems. Depending on the observables used, this quantity can inform on local properties of attractors such as periodicity, stability and persistence in phase space, or on global properties such as the Lyapunov exponents. We also introduce a new estimator of the extremal index and shows its relation with those previously introduced in the statistical literature. We reserve a particular focus to the systems perturbed with noise as they are a good paradigm of many natural phenomena. Different kind of noises are investigated in the annealed and quenched situations. Applications to climate data are also presented.

Figures

Equations

$$H_n(x):=\text{{first time the iterate Tn(x) enters the open ball B(z,e−un)}},$$

$$n\mu (B(z,e^{-u_n}))\to \tau ,$$

$$\mu (x;H_n(x)>u_n)\to e^{-\theta \tau },$$

$$\{M_n(x)\le u_n\},\text{where}{M}_n=\max \{\varphi (x),\varphi (Tx),\dots ,\varphi (T^{n-1}x\},$$

$$\varphi (x)=-\log d(x,z).$$

$${\Omega }_n^{(k)}(z):=\{\varphi (x)>u_n,\underset{i=1,...,k}{\max }\varphi (T^{i}x)\le u_n,\varphi (T^{k+1}x)>u_n\},$$

$$q_k:=\underset{n\to \infty }{\lim }\frac{\mu ({\Omega }_n^{(k)}(z))}{\mu (B(z,e^{-u_n}))}.$$

$$\theta =1-\sum _{k=0}^{\infty }{q}_k.$$

$$q_{p-1}=\frac{1}{|D{T}^p(z)|}.$$

$${\theta }^{-1}=\underset{n\to \infty }{\lim }\sum _{j=1}^{\infty }j{\pi }_n(j),$$

$$A_n^{(k)}:=\{\varphi (x)>u_n\cap \underset{i=1,\dots ,k}{\max }\varphi (T^{i}x)\le u_n\}$$

$${\theta }_k^{(n)}:=\frac{\mu (A_n^{(k)})}{\mu (\varphi (x)>u_n)}$$

$$q_k={\theta }_k-{\theta }_{k+1},$$

$${\theta }_k=\underset{n\to \infty }{\lim }{\theta }_k^{(n)},$$

$$D_{u_n}^{(k)}=\underset{n\to \infty }{\lim }n\sum _{j=k+1}^{[n/l_n]-1}\mu (A_n^{(k)}\cap T^{-j}(A_n^{(k)}))=0,$$

$$\mu (M_n\le u_n)-e^{-\tau {\theta }_k^{(n)}}\to 0,$$

$$q_j={\theta }_j-{\theta }_{j+1}=0,\forall j\ge k.$$

$$\theta =\underset{k\to \infty }{\lim }{\theta }_k.$$

$$\theta =\underset{k\to \infty }{\lim }\left(1-\sum _{j=0}^k{q}_j\right)=\underset{k\to \infty }{\lim }{\theta }_{k+1}.$$

$$\theta =\underset{k\to \infty }{\lim }\underset{n\to \infty }{\lim }{\theta }_k^{(n)}.$$

$${\hat{\theta }}_m:=\frac{{\sum }_{i=0}^{N-1-m}1(\varphi (T^i{x}_0)>u\cap {\max }_{j=1,...,m}\varphi (T^{j+i}{x}_0)\le u)/(N-m)}{{\sum }_{i=0}^{N-1}1(\varphi (T^i{x}_0)\ge u)/N}.$$

$${\hat{\theta }}_{Su}=\frac{{\sum }_{i=1}^{N_c}(1-p)S_i+N+N_c-\sqrt{({\sum }_{i=1}^{N_c}(1-p)S_i+N+N_c{)}^2-8N_c{\sum }_{i=1}^{N_c}(1-p)S_i}}{2{\sum }_{i=1}^{N_c}(1-p)S_i},$$

$${\hat{q}}_j=\frac{{\sum }_{i=0}^{N-2-m}1(\varphi (T^i{x}_0)\ge u\cap {\max }_{l=1,\dots ,m}\varphi (T^{l+i}{x}_0)<u\cap \varphi (T^{i+j+1}{x}_0)\ge u)/(N-1-m)}{{\sum }_{i=0}^{N-1}1(\varphi (T^i{x}_0)\ge u)/N}.$$

$$f_{\omega }^n:x=f_{{\sigma }^{n-1}\omega }\circ \cdots \circ f_{\omega }.$$

$$\omega \to f_{\omega }\text{such that}{f}_{\omega }(x)=T(x)+\omega -\mathrm{mod}1.$$

$$\sigma \omega \to f_{\sigma \omega }\text{such that}{f}_{\sigma \omega }(x)=T(x)+\sigma \omega -\mathrm{mod}1.$$

$$⋮$$

$${\sigma }^k\omega \to f_{{\sigma }^k\omega }\text{such that}{f}_{{\sigma }^k\omega }(x)=T(x)+{\sigma }^k\omega -\mathrm{mod}1.$$

$$q_0=\underset{n\to \infty }{\lim }\frac{\int d\mathbb{G}(\xi ){\rho }_s(H_n(0)\cap f_{\xi }^{-1}{H}_n(0))}{{\rho }_s(H_n(0))}=\underset{n\to \infty }{\lim }\left[\frac{1}{2}\frac{{\rho }_s(H_n(0)\cap f_0^{-1}{H}_n(0))}{{\rho }_s(H_n(0))}+\frac{1}{2}\frac{{\rho }_s(H_n(0)\cap f_1^{-1}{H}_n(0))}{{\rho }_s(H_n(0))}\right]$$

$$q_k\ge \underset{n\to \infty }{\lim }\frac{\mathbb{G}({\xi }_1,\dots ,{\xi }_{k+1}){\rho }_s(H_n(0)\cap (f_{{\xi }_1}^{-1}{H}_n(0){)}^c\cap \cdots \cap (f_{{\xi }_k}^{-k}{H}_n(0){)}^c\cap f_{{\xi }_{k+1}}^{-k-1}{H}_n(0))}{{\rho }_s(H_n(0))}=\mathbb{G}({\xi }_1,\dots ,{\xi }_{k+1}))/2^k.$$

$$f_{{\xi }_m}\circ \cdots \circ f_{{\xi }_2}\circ f_{{\xi }_{m^{\prime }}^{\prime }}\circ \cdots \circ f_{{\xi }_2^{\prime }}\circ f_1(0)=f_{{\xi }_m}\circ \cdots \circ f_{{\xi }_2}(p/q)=0.$$

$$q_k=\frac{1}{2^{k+1}}\underset{n\to \infty }{\lim }\frac{{\rho }_s(H_n(0)\cap (f_1^{-1}{H}_n(0){)}^c\cap ...\cap (f_1^{-k}{H}_n(0){)}^c\cap f_0^{-k-1}{H}_n(0))}{{\rho }_s(H_n(0))}=\frac{1}{4^{k+1}}.$$

Full text

, P. Yiou111Laboratoire des Sciences du Climat et de l’Environnement, UMR 8212 CEA-CNRS-UVSQ, IPSL and Université Paris-Saclay, 91191 Gif-sur-Yvette, France. E-mail: .

1. 1 Introduction

2. 2 The deterministic case

3. 2.1 Statistical estimators of the extremal index

4. 2.2 Numerical estimates of the extremal index in deterministic systems

5. 3 The random case: non-stationary situations (quenched noise)

6. 3.1 Sequential systems

7. 3.2 Random fibred systems

8. 3.2.1 Random Lasota-Yorke maps

9. 3.2.2 Rotations

10. 4 The random case: stationary situations (annealed noise)

11. 4.1 I.i.d. random transformations

12. 4.1.1 Additive noise

13. 4.1.2 Discrete noise

14. 4.2 Moving target

15. 4.2.1 Discrete noise

16. 4.2.2 Continuous noise

17. 4.2.3 Observational noise

18. 5 The dynamical extremal index

19. 5.1 Additive noise

20. 5.2 Discrete noise

21. 6 Point processes: statistics of persistence

22. 6.1 The deterministic case

23. 6.2 The random case

24. 6.3 The rotational case

25. 6.4 Climate data

26. 7 Conclusions

27. 8 appendix

28. 9 Acknowledgement

1. Introduction

In the last ten years extreme value theory (EVT) has been successfully applied to the study of time series generated by dynamical systems. Illustrations can be found in the book [13] and the articles [16, 36, 35] for an exhaustive account on the formalism, the methodologies and several applications. In particular the extremal index (EI) — a quantity defined in the unit interval — has been used as a powerful statistical indicator to discriminate among different qualitative types of dynamical behaviors in the phase space of a few climate models and in real situations (an extensive overview can be found in [37]). In a series of papers [17, 7, 19, 37, 18], the EI was renamed as local persistence indicator, suitable to estimate the average cluster size of the trajectories within the neighborhood of a given point. The aim of this note is to give an overview on a few recent rigorous mathematical results which show that the EI is very sensitive to stochastic perturbations affecting the deterministic evolution of a system states. We will also discuss the influence of noise on the dynamical extremal index, which was recently introduced to characterize the local divergence of chaotic orbits[20].

The first rigorous computations of the EI for expanding and uniformly hyperbolic systems (from now on named chaotic systems), showed a dichotomy for the value of the EI, which is equal to $1$ everywhere, except on periodic points. [24] give a detailed account of this matter. On the other hand, the EI exhibits local variations which have been also connected to the local fractal dimensions. By assuming that the data we are interested in are chaotic, one could wonder about the wider variation of the extremal index, offering a spectrum of values beyond the previous dichotomy. In a recent paper [8] we commented about a similar effect concerning the computation of the local dimensions. Indeed, while it is well known for a large class of systems that the local dimensions are all equal to the so-called information dimension $D_{1}$ with probability one, large deviation theory estimates the likelihood of deviations from this value at finite resolution. In this perspective, the spread of the experimentally observed values of the local dimensions can be thought of as originating from the multifractal structure of the invariant measure, which in turn is revealed by the non-constant value of generalized dimensions.

Although finite resolution affects the estimation of the extremal index, it is unlikely that it enjoys some large deviation property. This is due to the previous dichotomy, which prevents the existence of a smooth spectrum of values for the EI. Yet, the computation of the EI is very sensitive to randomness and disturbances in the measuring process. The main object of our note is to show how the EI is affected by the presence of noise. Our starting point is to assume for the EI a general formula obtained by Keller and Liverani [33] in the context of stationary dynamical systems enjoying the so-called spectral gap property. For this reason we call it the spectral formula. It encompasses all the other rigorous formulae previously obtained for the computation of the extremal index. We will show that when the spectral formula holds, one can obtain a generalized version of the O’Brien formula, which inspired several numerical algorithms and statistical estimators (see Eq. (3.2.4) in [13]). This is discussed in Proposition 2 and Remark 2, which constitute two main contributions of this paper.

The first sections are devoted to the computation of the EI for the standard observable given by the distribution of the first visit of the trajectory in a decreasing net of balls around a given point, which we qualify as the target set. The system will be perturbed according to three different classes of noises: sequential, quenched and annealed. In particular, we will focus on the annealed class since it gives stationary random processes particularly suitable for the application of the spectral formula. We will show on simple but non trivial examples that noises with discrete distributions tend to make the EI less than one, while noises with continuous distributions usually make the extremal index equal to one. The latter case applies also to deterministic dynamics with the target set affected by disturbances with smooth densities, which corresponds to common physical situations.

The final section deals with a different extremal index we introduced recently with the name of Dynamical Extremal Index and that explores the mutual distances between the coordinates of a point moving in a suitable higher dimensional direct product space. This new index captures the rate of phase space contraction and is naturally related to positive Lyapunov exponent(s). We will show again how this index is sensitive to the different kinds of perturbations described above.

We finally analyze how the random perturbations affect the statistics of the number of visits in the target regions. The expectation of the limiting law is the reciprocal of the extremal index and the differences in the value of the EI reflect in different types of Poisson compound distributions.

2. The deterministic case

Our intent is to provide a critical discussion of the application of the EI to time series. Therefore, we start defining it for particular random processes. Let us therefore consider a discrete dynamical systems defined by a map $T$ acting on a smooth manifold $X$ and preserving a Borel measure $\mu$. Suppose $z$ is a point of $X$ and $B(z,r)$ an open ball around $z$ and of radius $r$. Given any other point $x\in X$ let us consider the random variables $H_{n}(x)$ given by

[TABLE]

where the boundary level $u_{n}$ is defined by asking that

[TABLE]

for some positive number $\tau$. Notice that by the stationarity (or invariance) of the measure $\mu$, the quantity $\mu(B(z,e^{-u_{n}}))$ gives the probability that any iterate of the map $T$ be in the open balls $B(z,e^{-u_{n}})$. Moreover as soon as the measure is not atomic, the measure of a ball varies continuously with the radius and this allows us to explicit $u_{n}$ as a function of $\tau$ and $n$.

We say that we have an extreme value law for the process $H_{n}$ if

[TABLE]

where $0\leq\theta\leq 1$ is called the extremal index.

We remind that the events $H_{n}>u_{n}$ are equivalent to the following ones:

[TABLE]

and the observable $\phi$ is defined as

[TABLE]

It is remarkable that for a large class of systems whose transfer operator admits a spectral gap, the extremal index can be explicitly computed. This is done by G. Keller [32], who applied the perturbative theory developed in [33]. This spectral approach has the advantage to give a formula for the EI which holds for general target sets when their measure goes to zero. We will use it in the last section of this note when balls are replaced by tubular neighborhoods.

Let us first define the event

[TABLE]

namely the set of points in $B(z,e^{-u_{n}})$ whose first $k$ iterates are outside $B(z,e^{-u_{n}})$ and whose $(k+1)$-th iterate falls again in $B(z,e^{-u_{n}})$.

Suppose that the following limit exists:

[TABLE]

Then the perturbative theory gives an estimate of the EI:

[TABLE]

This formula allows us to reproduce a few existing results whenever the target point $z$ is of a particular type. For instance, for one-dimensional expanding systems with strong mixing properties and preserving a measure absolutely continuous with respect to Lebesgue, the extremal index is $1$ everywhere but in periodic points. When we have a periodic point $z$ of (minimal) period $p,$ only the $q_{p-1}$ term is nonzero and it is of the form [32, 24]

[TABLE]

The dependence of the EI on the periodicity of the target point $z$ means that the sojourns and returns of the point into the ball $B(z,e^{-u_{n}})$, which constitute the clusters of exceedances, keep memory either of the past orbit and of the topological structure of the target point.

This can be made more precise by defining a cluster size distribution $\pi_{n}(j)$ as the probability of having $j$ returns into the ball up to a rescaled time $n/k_{n}$, for a suitable sequence $k_{n}=o(n)$. Section 3.2.1 in [13] provides more details on this derivation; a more formal definition is also given in section 7. It can be therefore proved that (see Sec. 3.3.3 in [13]):

[TABLE]

which is interpreted by saying that the EI is equal to the inverse of the average cluster size.222Abadi et al. [1] built a dynamically generated stochastic processes with an extremal index for which that equality does not hold. They considered observable functions maximised at least two points of the phase space, where one of them is an indifferent periodic point and another one is either a repelling periodic point or a non periodic point. We will not consider these kind of observables in this paper. In the applications to dynamical systems, the periodicity translates into the fact that the clustering $\pi=\lim_{n\rightarrow\infty}\pi_{n}$ is actually a geometric distribution of parameter $\theta\in(0,1]$, i.e., $\pi_{k}=\theta(1-\theta)^{k-1}$, for every $k\in\mathbb{N}_{0}$ [30]. We will return to this matter in section 7.

The interpretation of the EI given in Eq. (8) has been used in applications to time series climate data and it was also emphasized the local character of such an indicator and its strong correlation with the fractal local dimension of the invariant measure [17, 7, 19].

Let us now consider for $k\in\mathbb{N}_{0}$ the event (compare with Eq. (6)):

[TABLE]

and the quantity:

[TABLE]

where we set $\theta_{0}^{(n)}=1,\forall n$. By introducing the event $\phi(T^{k+1}\leq u_{n})$ and its complement, and passing to the limit for $n\rightarrow\infty$, we get

[TABLE]

where

[TABLE]

when the limit exists.

Notice that as soon as one of the two sequences $\{q_{k}\}_{k\in\mathbb{N}},\{\theta_{k}\}_{k\in\mathbb{N}}$ is known, the other one is determined as well.

In some circumstances, the limit in Eq. (14) gives exactly the extremal index. To explain this point, we need to be more precise about the assumptions we make. The spectral approach briefly sketched above applies to a large class of dynamical systems admitting a spectral gap for the transfer (Perron-Fröbenius) operator. This means that the systems have exponential decay of correlations for smooth enough observables, usually of bounded variation type. This is not sufficient to establish the existence of the limit in Eq. (7). Computations of that limit under various circumstances are given in [32, 33, 4, 8]. The more standard approach to EVT, which consists in adapting the classical Leadbetter theory for i.i.d. random variables to stationary dependent processes [34], requires milder mixing conditions that allow to get asymptotic independence for the process. The theory of O’Brien [38] was developed for strictly stationary processes verifying asymptotic independence, which is a strong probabilistic mixing condition like the $\phi$-mixing (see [6] for the definition). We stress that asymptotic independence is very difficult to check in practice. Recently it has been shown that weaker conditions are enough. For our purposes, it is sufficient to remind condition $D(u_{n})$ in section 2.3 of [13] successively improved in [24]333This condition holds for instance when the invariant measure $\mu$ is mixing with decay of correlation fast enough; sometimes a rate of decay as $n^{-2}$ is sufficient. Asymptotic independence is not enough to get the convergence of the limit distribution of the maxima. One needs to control short returns, which could be in particular affected by clustering. A powerful condition which takes care of that is:

[TABLE]

where $l_{n}=o(n)$ is a sequence slowly diverging to infinity and verifying the $D(u_{n})$ conditions (see Remark 4.1.2 in [13] for the details). We call the integer $k$ the clustering order. We now summarize in the following proposition a few facts that are useful for our paper and which refer to the stationary process $\{\phi\circ T^{k}\}_{k\in\mathbb{N}}$, on the probability space $(X,\mu)$:

Proposition 1**.**

• •

Suppose the sequence $u_{n}$ verifies Eq. (2) and condition $D(u_{n})$ holds. Moreover suppose that $\liminf_{n\rightarrow\infty}\mu(M_{n}\leq u_{n})>0$ and condition in Eq. (15) holds with clustering order $k$. Then

[TABLE]

as shown by **[10]**, Proposition 1.1.

• •

If the limit in Eq. (14) exists, then the extremal index is given by $\theta_{k}$, as shown by Corollary 1.3 in **[10]** or Corollary 4.1.7 in **[13]**. This is sometimes called the O’Brien formula.

• •

If condition Eq. (15) holds for some particular $k^{\prime}\in\mathbb{N}$, it also holds for all $k\geq k^{\prime}$ and therefore the limit in Eq. (14) for all $k\geq k^{\prime}$ gives again the extremal index $\theta$, as shown by **[24]** and **[13]**, Remark 4.1.10.

When the second item in the Proposition holds for a process with clustering order $k$ we have

[TABLE]

The natural question is therefore to ask when condition (15) is verified. For a large class of dynamical systems, in particular verifying the assumptions of theorem 4.2.7 in [13], it can be shown that condition (15) holds with $k=0$ if the target point $z$ is not periodic, and with $k=p$ if $z$ is periodic of prime period $p$ (see Proposition 4.2.13 in [13]). In particular, for one-dimensional expanding systems with strongly mixing properties and preserving a measure absolutely continuous with respect to Lebesgue, the extremal index is $1-q_{p-1}=1-\frac{1}{|DT^{p}(z)|}$ in periodic points $z$ of (minimal) period $p$. For higher dimensional systems, condition (15) gives a precise value for $\theta$ (see for instance [9]).

Remark 1**.**

For the preceding example around periodic points, the spectral and the standard approaches both give $q_{p-1}$ and $\theta_{p}$ different from zero. Yet, the spectral approach shows that all $q_{j}$ are zero for $j\neq p-1$. This is therefore coherent with Eq. (16) for which $\theta_{j}=\theta_{j+1}$ for $j\geq p$. Let us now consider $q_{p-1}=\theta_{p-1}-\theta_{p}$. The quantity $\theta_{p-1}$ is well defined as the limit of $\theta_{p-1}^{(n)},$ when $n$ goes to infinity and it is equal to $1$. This implies that the condition (15) is violated, because otherwise $\theta_{p}$ should be equal to $1$ as well. The fact that all the $\theta_{j}$, for $j\leq p-1$ are equal to $1$ is not surprising: it means that the full conditional measure of points in the ball $B(z,e^{-u_{n}})$ are outside it when iterated up to $p-1$ times.

In the assumptions of Proposition 1, the EI is equal to $\theta_{k}$ when the limit in Eq. (14) defining this quantity exists and $k$ being the clustering order. We now show that this is a particular case of a more general result.

Proposition 2**.**

Let us suppose that the limits in Eq. (7) defining the quantities $q_{k}$ exist for any $k\geq 0$. Then the extremal index $\theta$ is given by

[TABLE]

Proof.

We noticed that whenever the sequence $q_{k}$ exists, the same happens for the sequence $\theta_{k}$. Then by a simple telescopic trick we get

[TABLE]

∎

Remark 2**.**

Under the assumptions of the preceding proposition, we got that

[TABLE]

We stress that the limits in Eq. (14) or Eq. (7) have been computed up to now and in the framework of dynamical systems only around periodic points, or around the diagonal in coupled systems (see [15]) where at most only one of them is different from zero. We will present later on further examples of explicit computations of the $q_{k}$. In particular section 4.1.2 will give an example of a stationary random dynamical system for which all the $q_{k}$ are different from zero and from each other, which implies that condition (15) is violated at all orders. In this respect our Eq. (18) is a generalization of the O’Brien formula.

2.1. Statistical estimators of the extremal index

We saw in the preceding section that the value $\theta_{p}$ provides the extremal index if the target point $z$ is of (minimal) period $p$. It is therefore not surprising that estimators based on the O’Brien formula have been proposed for the computation of the EI. In particular the following formula has been introduced to compute the extremal index associated to time series in [5] (formula 10.26): we start with a trajectory $(x_{0},Tx_{0},\dots,T^{N-1}x_{0})$ and for a high enough threshold $u$ we compute:

[TABLE]

If we are in presence of a periodic point of period $m$ for the observable $\phi(x)=-\log d(z,x)$, then $\hat{\theta}_{m}$ is a good approximation of the EI.

We notice that when $m=1$, Eq. (19) gives the likelihood estimator of Süveges $\hat{\theta}_{Su}$ [39], and the obtained result is an indicator of the persistence of the system close to a state $z$ (see for example [17]). We remind that for a high quantile $p$ of the time series, the Süveges estimator is given by:

[TABLE]

where $N$ is the number of exceedances over $\tau$, $N_{c}$ is the number of clusters of two or more exceedances, and $S_{i}$ is the size of the cluster $i$. This estimator is presented in [39]. A Matlab code and numerical details related to this computation in the context of dynamical systems are available in [13].

We warn that in presence of periodicity of order, say $p$, the estimator of Eq. (19) is reliable when it is strictly less than $1$, because all the $\theta_{m}$, for $m<p$ are equal to $1$. This observation applies as well to the computation of the $q_{k}$: they are all zero up to $q_{p-1}$.

Ferreira [22] reviews several numerical algorithms to compute the EI (including those shown before). In this paper we choose to work with the $q_{k}$ and we list a few features of this approach.

• •

In the presence of periodicity, we do not need to check condition Eq. (15); the clustering order, if any, is obtained automatically by progressing in the evaluation of the $q_{j}$. We will show in sections 4.1.2 and 4.2 examples of stationary processes constructed with random perturbation where several $q_{k}$ are different from zero.

• •

Eq. (8) works as well when the ball $B(z,e^{-u_{n}})$ is replaced by any other measurable sets whose measure goes to zero when $n\rightarrow\infty$. This comes to modify suitably the observable $\phi$. A concrete example will be given in section 5.

• •

Whenever one of the $q_{j}$ is strictly positive, we get an EI strictly less than $1$. This will be extensively used in sections 4.1.2, 4.2, 5.2. Notice that it is equivalent to find a certain $\theta_{k}<1$ for $k\geq 1$.

In practice, one cannot compute all the $q_{j}$ values and must stop at an order $m$. We stress that the sequence $\{q_{j}\}$ is rapidly decreasing to [math] and the quantity $\theta_{m}=1-q_{0}-\cdots-q_{m-1}$, where $m$ should be taken not too big, should give a fair estimate of $\theta$. The approximate value $\hat{q_{j}}$ for $q_{j}$ is obtained by considering a trajectory of length $N$ as

[TABLE]

2.2. Numerical estimates of the extremal index in deterministic systems

We present below the computations for different maps and various target points $z$ and show the difference between the approximate estimate $\hat{\theta}_{m}$ of $\theta_{m}$ (we performed our computations at the order $m=5$) and Süveges estimate (with order [math]). We find that $\hat{\theta}_{m}$ obtains values close to the predictions of Eq. (9), while the other estimate gives $1$ for points of period larger than $1$.

[TABLE]

We observe in figure 1 that for a target point of period $k$, $\hat{\theta}_{m}$ is equal to 1 for $m<k$, is equal to $\theta$ for $m\geq k$, due to the fact that $\hat{q}_{k-1}$ is non zero. For $m\neq k$, we also have $\hat{q}_{m-1}=0$, as expected. Although the Süveges estimator is not in principle suitable to detect periodic points with period strictly larger than $1$, it is very often used to analyze time series since it is particularly simple to implement. But there is a more interesting reason. If the process is not periodic, for instance generated by random dynamical systems like those presented in the next sections, whenever the Süveges estimator is strictly less than one, the same happens for the EI, thanks to our Proposition 2.

Instead, as we saw in the table above, the difference with real situations of larger periodicity could be important. In such cases, the estimate $\hat{\theta}_{m}$ performs systematically better. We will use this estimator in the numerical computations of this paper.

To evaluate the differences between the two estimates in high dimensional datasets, we test them on the North Atlantic daily sea-level pressure data described in [17]. For each atmospheric state $z$, the EI index associated to the observable $-\log\text{dist}(z,.)$ was evaluated using the Süveges estimate and called the inverse persistence of $z$. The two estimations are shown in Figure 3

Estimates using $\hat{\theta}_{5}$ for these data are systematically lower than those with the order [math] method (in average $0.065$ less), due to the contribution of some $\hat{q}_{k}$ for $k>0$ (see the empirical distribution in figure 2). In Figure 3 we present a scatter plot of the daily values of $\theta$ obtained with the two estimators. We see that there is a strong linear relation between the two estimates with an offset of about 0.1 days*-1*. The estimation of the cross-correlation coefficient, namely the zeroth lag of the normalized covariance function, yields 0.94, meaning that the information contained in the two estimators is the same except for a restricted set of states $z$. Moreover, if we remove the time averaged values $\langle\hat{\theta}_{Su}\rangle$ of the Süveges estimator and $\langle\hat{\theta}_{5}\rangle$ of the new estimator from the time series, we can compare the distribution of $\hat{\theta}_{Su}-\langle\hat{\theta}_{Su}\rangle$ with that of $\hat{\theta}_{5}-\langle\hat{\theta}_{5}\rangle$. For this comparison, we have used a two sided Kolmogorov-Smirnov test and verified the null hypothesis that the estimators are from the same continuous distribution at the 5$\%$ confidence level.

The lowest differences between the two estimators are found for patterns close to the average field $z$. This corresponds to anomalies close to 0 hPa in Figure (3b). Instead, large differences between estimators correspond to a peculiar $z$ pattern, consisting of a deep low pressure anomaly over the north of the domain and an anticyclonic anomaly over the southern part of the domain. This pattern resembles those observed in [17] for the minima of the local dimensions $d(z)$. This analysis confirms that the two estimators are different even beyond the simple shift in the values. The patterns where the two estimators mostly diverge can be thought as the pre-images of higher order periodic points of the attractor underlying the mid-latitude atmospheric circulation. In this sense, the combined use of both $\hat{\theta}_{Su}$ and $\hat{\theta}_{5}$ estimators could be useful to detect these special points of the dynamics in climate data and other natural phenomena.

3. The random case: non-stationary situations (quenched noise)

The dynamics and the consequent detection of the extremal index could be affected in different ways. We begin with the dynamics by distinguishing first two cases:

3.1. Sequential systems

These systems are defined by concatenating maps chosen in some set, usually in the close neighborhood of a given map. As a probability measure one usually takes some ambient measure like Lebesgue, which makes the process given by the sequences of concatenations non-stationary. The extremal value theory and the extremal index must be redefined (see [25], in particular Eq. (2.8)). In section 4.5 of the aforementioned paper, we showed an example of a sequential system modeled on maps chosen in the neighborhood $\mathcal{U}$ of a given $\beta$-transformation $T_{\beta}$. According to suitable choices of $\mathcal{U}$ it is possible to show either that the EI of the unperturbed map and of the sequential system are the same, or the two differ, in particular when the elements of the concatenation are far enough from $T_{\beta}$. In this case the EI is simply $1$. Sequential systems are a model for non-autonomous physical systems that do not admit a natural stationary probability measure. Statistical limit theorems can be proved as well in this context provided the observables are suitably centered (see [12, 29] for more details). The previous example with $\beta$-transformations suggests that when the maps are far enough from each other, the EI converges to $1$.

To the best of our knowledge, there are not numerical investigations of extreme value behaviors for sequential systems. We propose one of them in the Appendix.

3.2. Random fibred systems

The sequential systems are in some sense too general, since there are no real prescriptions on the choice of the maps. A much more realistic class of non-autonomous systems is given by random fibred systems. They are constructed by taking a driving map $\sigma$ that preserves a probability measure $\nu$ on the measurable space $\Omega$, and which codes a family of transformations $f_{\omega},$ for $\omega\in\Omega$ on the fiber $X$ via the composition rule

[TABLE]

The driving system encodes the external influences on the system of interest: it acts deterministically in the choice of the evolution transformations on the fibers. As L. Arnold wrote at the beginning of his monograph [2]:

Imagine a mechanism which at each discrete time $n$ tosses a (possibly complicated, many-sided) coin to randomly select a mapping $\phi_{n}$ by which a given point $x_{n}$ is moved to $x_{n+1}=\phi(x_{n}).$ The selection mechanism is permitted at time $n$ to remember the choices made prior to $n$, and even to foresee the future. The only assumption made is that the same mechanism is used at each step. This scenario, called product of random mapping, is one of the prototypes of a random dynamical systems.

The mechanism used at each step is just the driving system. It is interesting to observe for the objectives of this paper, that random fibred systems have been used to analyze the transport phenomena in non-autonomous dynamical systems, such as geophysical flows (see the enlightening review article [27]). In order to perform statistics on these systems, we notice that a family of sample measures $\mu_{\omega}$ lives on the fiber, which verifies the quasi-invariant equation $(f_{\omega})_{*}\mu_{\omega}=\mu_{\sigma\omega}$, where $(f_{\omega})_{*}$ is the push-forward of the measure. These sample measures will be taken as the probability measures that describe the statistical properties along the fiber and they do not give rise to stationary processes. The relation between the sample measures and $\nu$ is that the measure $\mu:=\int\mu_{\omega}d\nu(\omega)$ is preserved by the skew (deterministic) transformation $F(x,\omega)=(f_{\omega}(x),\sigma\omega),$ acting on the product space $X\times\Omega.$ The extremal index can be suitably defined and it could be different from $1$ (see Corollary 5.1 and Theorem 5.3 in [25] dealing with random subshift). Transparent examples are constructed in the following subsections.

3.2.1. Random Lasota-Yorke maps

[23]. Take $\Omega=I^{\mathbb{Z}},$ where $I=(1,\dots,m)$ is a finite alphabet with $m$ letters. We associate to each letter a piecewise expanding map of the interval $f_{\omega}$ satisfying some smoothness and distortion standard assumptions. The map $\sigma$ is therefore the bilateral shift and $\nu$ any ergodic shift-invariant non-atomic probability measure, for instance the Bernoulli measure with weights $p_{1},\dots,p_{m}$. The sequence of concatenations is given by $f^{k}_{\omega}=f_{\omega_{k}}\circ\cdots\circ f_{\omega_{1}}$, with $\omega_{j}\in I$ and $j=1,\dots,k$ are the first $k$ symbols of the word $\omega$.444Notice that this is equivalent to Eq. (22) by defining the map $\omega\rightarrow f_{(\omega)^{\prime}}$, where $(\omega)^{\prime}$ is the $\omega_{1}$ coordinate of $\omega$. In this setting it can proved that the sample measures are equivalent to the Lebesgue measure $m$ and for $m$-almost any target point $z\in X$ and $\nu$ almost any realization $\omega\in\Omega$, the EI is equal to $1$ (see [23]). This corresponds to a quenched result since it depends on the choice of the sequences $\omega$.

3.2.2. Rotations

Another more physical quenched example is constructed in the following way. Let us take as $\Omega$ the unit circle $\mathbb{S}^{1}$ and as $\sigma$ the irrational rotation: $\sigma(\omega)=\omega+\alpha-\mathrm{mod}1$, $\omega\in\mathbb{S}^{1}$ with $\alpha\in\mathbb{R}$. Then we define $T(x)=3x-\mathrm{mod}1$ and make the correspondence:

[TABLE]

and so on, by composing after that as: $f_{\sigma^{k}\omega}\circ\cdots f_{\omega}$.

This last example was implemented numerically using the estimate $\hat{\theta}_{5}$ for different trajectories and choices of $\alpha$, either irrational or rational. We found an extremal index equal or very close to $1$ in all cases (see table 2).

Remark 3**.**

In the preceding two examples we considered the extremal index, but we did not indicate any method to compute it. We are in fact outside the stationary regime which provided us with the various formulae described in the previous sections. Actually a definition of the EI in the non-stationary case is given in formula (2.8) in [25]. This formula relies also on a suitable definition of the boundary level $u_{n}$ given in Eq. (2.2) in the aforementioned paper. The above formula (2.8) is a slight modification of the O’Brien formula and it could be reasonably recovered numerically by using Eq. (19) or the Süveges formula when we are not in presence of clustering. Although the EI was computed rigorously for the Random Lasota-Yorke maps as indicated in section 3.2.1, we do not dispose of a rigourous argument for the rotations and the numerical results we presented have been obtained with the $\hat{\theta}_{5}$ estimate. We observe that we could assume as a definition of the extremal index in these non-stationary situations the reciprocal of the expectation given by the distribution of the number of visits (see Eq. (46) in section 6). We will show in section 6.3 that the value for the EI computed in that way is the same as that provided by the $\hat{\theta}_{5}$ estimate.

[TABLE]

4. The random case: stationary situations (annealed noise)

The preceding two situations dealt with non-stationary processes. We now focus on stationary processes given (i) by i.i.d. randomly chosen transformations, and (ii) by i.i.d. moving target sets.

4.1. I.i.d. random transformations

We now suppose that the maps $f_{\omega}$ are no longer driven by the measure-preserving map $\sigma$, but are chosen in an i.i.d. manner.

4.1.1. Additive noise

A common way to build such a process is to choose a map $f$ and to construct the family of maps $f_{\xi}=f+\epsilon\xi$, where $\xi$ a random variable sampled from some distribution $\mathbb{G}$. This distribution can be the uniform distribution on some small ball of radius $\epsilon$ around [math], where $\epsilon$ is the intensity of the noise. The iteration of the single map $f$ is now replaced by the concatenation $f_{\xi_{n}}\circ\cdots\circ f_{\xi_{1}}$, where the $\xi_{k}$ are i.i.d. random variables with distribution $\mathbb{G}$. If $T$ has good expanding or hyperbolic properties, it is possible to show the existence of the so-called stationary measure $\rho_{s}$, verifying for any real bounded function $g$: $\int gd\rho_{s}=\int\int g\circ T_{\xi}d\rho_{s}d\mathbb{G},$ (see [13] Chap. 7, for a general introduction to the matter). If we now take the probability product measure $\mathbb{Q}:=\mathbb{G}^{\mathbb{N}}\times\rho_{s},$ any process of type $\{g(f_{\xi_{n}}\circ\cdots\circ f_{\xi_{1}})\}_{n\in\mathbb{N}}$ is $\mathbb{Q}$-stationary. The extremal index is obtained again by applying Eq. (8) where the sets $\Omega_{n}^{(k)}(z)$ and $\mu(B(z,e^{-u_{n}}))$ are now weighted with the measure $\mathbb{Q}$. The expectation being taken with respect to $\mathbb{G}$ too, makes this an annealed type random perturbation, while the fibred perturbation described above is a quenched one, the expectations depending on the realization $\omega.$ In the article [4] we rigorously proved that for a large class of piecewise expanding maps of the interval perturbed with additive noise with the uniform distribution $\mathbb{G}$, the extremal index is $1$ as soon as $\epsilon$ becomes positive. The same happens for other kinds of maps and smooth distributions, as we showed in [14] almost numerically. We believe that this mostly happens for the systems perturbed with i.i.d. noise, since the latter destroys all the periodic points. Nevertheless it is possible to construct annealed examples giving an EI less than one. This will be done in the next section.

4.1.2. Discrete noise

Let us consider:

• •

Two maps $f_{0}=2x-\mathrm{mod}1$ and $f_{1}=2x+b-\mathrm{mod}1$, $0<b<1.$

• •

A fixed target set around [math]: $H_{n}(0):=B(0,e^{-u_{n}}),$ the ball of radius $e^{-u_{n}}$ around [math]. Notice that [math] is a fixed point of $f_{0}$, but not of $f_{1}$.

The indices $\xi_{1},\dots,\xi_{n}$ in the concatenation $f_{\xi_{n}}\circ\cdots\circ f_{\xi_{1}}$ are i.i.d. Bernoulli random variables with distribution, for instance, $\mathbb{G}=p\delta_{0}+(1-p)\delta_{1},$ with $p=\frac{1}{2}$. The spectral theory developed in [4] applies as well to this case, giving an absolutely continuous stationary measure $\rho_{s}$. Proposition 5.3 in [4] gives for $q_{0}$:

[TABLE]

Since [math] is not a fixed point of $f_{1}$, a standard argument gives [math] to the limit in the square bracket on the right hand side. Instead by the mean value theorem applied to the term in the left square bracket we get in the limit of large $n$: $q_{0}=\frac{1}{4}$.

Remark 4**.**

The preceding example could be interpreted as a perturbation of the map $f_{0}=2x-\mathrm{mod}1$ around its fixed point [math]. Instead we change the point of view and consider it as a perturbation of the map $f_{1}=2x+b-\mathrm{mod}1$ around [math]. Supposing the map $f_{0}$ is chosen with probability $p_{0}$, we see that when $p_{0}=0$ (absence of perturbation), the $EI=1$, since [math] is not a fixed point for $f_{1}$, but it jumps to $EI\leq 1-0.5\ p_{0}$ as soon as $p_{0}>0$ and regardless of the value of $b$. This behavior is specular of what happened for the additive noise described above: there it was enough to switch on the noise, no matter of its magnitude, to make the $EI$ equal to one. Here, the $EI$ changes in term of the probability of appearance of the perturbed map, no matter of its topological distance from the unperturbed one, the value of $|b|$ in this case. However, it is easy to make the $EI$ depending on the distance between the maps. For instance, let us take $f_{0}$ as $f_{0}=(2+j)x-\mathrm{mod}1$, $j\in\mathbb{N}$, which could be seen as a strong perturbation. By repeating the argument above we get that the $EI\leq 1-p_{0}\ \frac{1}{2+j}$.

The preceding computation of $q_{0}$ is valid for all values of $b\in(0,1)$. The question we may ask now is whether other $q_{k}$ are non zero. By a similar reasoning as for the computation of $q_{0}$, this problem is equivalent to ask for the existence of a concatenation $f_{\xi_{k+1}}\circ\cdots\circ f_{\xi_{1}}$ returning the point [math] to itself after $k+1$ iterations and not before. Indeed, in that case we have that

[TABLE]

If the random orbit will not return to [math], a standard argument gives the value [math] for all $q_{k}$ when $k\neq 0$.

To investigate the existence of such a concatenation, let us distinguish between the cases when $b$ is irrational or rational.

If $b$ is irrational and $\xi_{1}\neq 0$, any concatenation of any length maps [math] into an irrational point, and so the orbit cannot come back to [math] after having left it. For this reason, all the $q_{k}$ but $q_{0}$ are equal to [math].

We investigate the case where $0<b=p/q<1$ is rational, where $p,q\in\mathbb{N}$ are mutually prime.

We start by showing by induction that for $n\geq 1$ and for all $i\in[2^{n},2^{n+1}-1]$, the random orbit starting from $p/q$ can attain the points $\frac{ip}{q}-\mathrm{mod}1$ (and only them) in exactly $n$ steps.

This proposition is true at rank $n=1$, since $f_{0}(p/q)=2p/q$, $f_{1}(p/q)=3p/q-\mathrm{mod}1$. Suppose now that it holds at rank $n$. Then we have that

• •

$f_{0}(2^{n}p/q-mod1)=2^{n+1}p/q-\mathrm{mod}1$

• •

$f_{1}(2^{n}p/q-mod1)=(2^{n+1}+1)p/q-\mathrm{mod}1$

⋮

• •

$f_{0}((2^{n+1}-1)p/q-mod1)=(2^{n+2}-2)p/q-\mathrm{mod}1$

• •

$f_{1}((2^{n+1}-1)p/q-mod1)=(2^{n+2}-1)p/q-\mathrm{mod}1,$

and the proposition is true at rank $n+1$. This result implies that any point of the form $jp/q-\mathrm{mod}1$, with $j\in\mathbb{N},j\geq 2$ can be attained by the random orbit starting from the point $p/q$.

Let us consider $f_{\xi_{m}}\circ\cdots\circ f_{\xi_{2}}\circ f_{1}$ the concatenation of minimal length $m\geq 2$ such that $f_{\xi_{m}}\circ\cdots\circ f_{\xi_{2}}\circ f_{1}(0)=qp/q-\mathrm{mod}1=0$. Note that the first map involved in this concatenation is $f_{1}$, so that the first iterate gives $p/q$. By the preceding argument, there exists such a concatenation and as we consider the smallest of the concatenations having this property, the random orbit leaving [math] does not come back to [math] before $m$ iterations. As discussed earlier, this implies a non zero value of $q_{m-1}$.

In fact, there exists infinitely many non zero $q_{k}$: consider the concatenation of minimal length $m^{\prime}\geq 2,$ $f_{\xi^{\prime}_{m^{\prime}}}\circ\cdots\circ f_{\xi^{\prime}_{2}}\circ f_{1}$ such that $f_{\xi^{\prime}_{m^{\prime}}}\circ\cdots\circ f_{\xi^{\prime}_{2}}\circ f_{1}(0)=(q+1)p/q-\mathrm{mod}1=p/q-\mathrm{mod}1$. By the argument in the induction, we have that either $m^{\prime}=m$ or $m^{\prime}=m+1$. Then we see that

[TABLE]

This proves the existence of a concatenation returning the point [math] to itself after exactly either $2m$ or $2m+1$ iterations (and not before). Again we have proved that either $q_{2m}$ or $q_{2m-1}$ is non zero. Applying the same reasoning, we can prove that infinitely many $q_{k}$ have positive values.

An interesting situation happens when $b=1/2$. In this case, for all $k$, [math] can come back to itself after exactly $k+1$ iterations, and this happens only for the concatenation $f_{1}\circ\cdots\circ f_{1}\circ f_{0}$, where the map $f_{1}$ is applied succesively $k$ times. This sequence having probability $(1-p)^{k}p=1/2^{k+1}$, it is easily seen (by the very definition of $q_{k}$ in Eq. (7) and analogously to the proof for the computation of $q_{0}$), that the $q_{k}$ are non zero and are equal to:

[TABLE]

In this situation, $\theta$ can be computed analytically and we find:

[TABLE]

We computed in table 3 some estimates $\hat{q}_{k}$ and we find a very good agreement between theoretical and numerical estimates.

[TABLE]

4.2. Moving target

Another source of disturbance comes from the uncertainty of fixing the target ball for the hitting times of the orbit. A more general theory of random transformations and moving balls will be presented elsewhere. Here we simply consider the case of a deterministic dynamics and a target set shifted by a random displacement and we assume an annealed approach.

4.2.1. Discrete noise

Let us take a point $z\in X$; then we consider $\sigma$ the one-sided shift on $m$ symbols and set $\Omega=\{1,\cdots,m\}^{\mathbb{Z}}$. To each symbols, $j,j=1,\dots,m,$ there corresponds a $z(j)\in[z-\epsilon,z+\epsilon]$. As a probability $\mathbb{G}(\omega),\omega=(\omega_{0},\omega_{1},\dots)\in\Omega$, we put the Bernoulli measure of weights $p_{i}=\frac{1}{m},i=1,\dots,m$. It is invariant under the shift. Note that $\epsilon$ is fixed and is the strength of the uncertainty. At each temporal step the map $T$ moves $x$ and the target point moves as well around $z$ under the action of $\sigma$. We now construct the direct product $S(x,\omega):=(T(x),\sigma\omega)$. The probability measure $\mathbb{P}:=\mu\times\mathbb{G}$ is invariant under $S$. As in section 2 we are interested in the distribution of the random variable

[TABLE]

where $\omega^{\prime}=\omega_{0}$ is the first coordinate of $\omega$. We define the usual observable:

[TABLE]

and notice that for all $k\geq 1:$

[TABLE]

The distribution of $H^{r}_{n}(x)$ is given by $\mathbb{P}\left(M_{n}\leq u_{n}\right),$ where

[TABLE]

and the boundary levels $u_{n}$ verify

[TABLE]

for some positive number $\tau$. We now choose a map $T$ on the unit interval admitting a spectral gap for the transfer operator and preserving an absolutely continuous invariant measure $\mu$. That operator acts, for instance, on the space of bounded variation function on $[0,1]$, which can be enlarged to a Banach space by adding to the total variation, $|\cdot|_{TV}$, the $L^{1}$ norm with respect to the Lebesgue measure $m$, $||\cdot||_{L^{1}(m)}$. We endow the space of functions $g(x,\omega)$ defined on $X\times\Omega$ with the Banach norm defined by:

[TABLE]

One can show that with respect to this norm the transfer operator $\mathcal{L}_{S}$ for the direct product $S$ has a spectral gap on the largest eigenvalue $1$. We introduce now the perturbed operator $\tilde{\mathcal{L}}_{S}(g(x,\omega)=\mathcal{L}_{S}(g(x,\omega){\bf 1}_{B^{c}(z(\omega^{\prime}),e^{-u_{n}})}(x))$, where $g\in\mathcal{B}$ is a function defined on the Banach space just introduced. We first notice that

[TABLE]

where $h$ is the density of $\mu$, which, with our assumptions, is bounded from below by the constant $C$. Since $||(\tilde{\mathcal{L}}_{S}-\mathcal{L}_{S})h||_{L^{1}(\mathbb{G}\times m)}\leq Ce^{-u_{n}}||h||_{\mathcal{B}}$, we can apply the perturbation theory mentioned in section 2, (see [32], [15], [13] ch. 7), and show that $\mathbb{P}\left(M_{n}\leq u_{n}\right)$ converges to the Gumbel law $e^{-\theta\tau}$. The extremal index $\theta$ is given by the adaptation to the actual case of Eq. (7), with the $q_{0}$ which reads:

[TABLE]

We now give a simple example for which $q_{0}>0$, which implies that the extremal index is strictly less than $1$. Take as $T$ the map $T(x)=2x-\text{mod}1$, and an alphabet of $4$ letters: $\{0,1,2,3\}$ with equal weights $1/4$. Moreover $\mu$ is the Lebesgue measure. Then we set the associations:

• •

$0\rightarrow z_{0}$

• •

$1\rightarrow z_{1}$

• •

$2\rightarrow z_{2}$

• •

$3\rightarrow z_{3}$,

where $z_{i},i=0,1,2,3$ are points in the unit interval verifying the following assumptions:

• •

$T(z_{1})=T(z_{2})=z_{0};\ T(z_{0})=z_{3}.$

• •

$T(z_{3})\neq z_{i},i=0,1,2,3.$

The numerator of $q_{0}$ is:

[TABLE]

[TABLE]

[TABLE]

where $[\omega_{0},\omega_{-1}]$ denotes a cylinder with fixed coordinates $\omega_{0}$ and $\omega_{-1}$.

For the denominator we get the value:

[TABLE]

By using the same arguments as in Eq. (23), we see that for all the cylinders in the numerator different from the three explicitly given above, the integrals give zero in the limit of large $n$. For the first of the three cases above (the others being similar), we get

[TABLE]

The numerator therefore contributes with $\frac{3}{2}(\frac{1}{4})^{2}\mu\left(B(z_{0},n)\right)$. The denominator is (by the Lebesgue translation invariance): $\mu\left(B(z_{0},n)\right)$. Therefore we get $q_{0}=\frac{3}{2}(\frac{1}{4})^{2}=0.09375.$

To get $q_{1}$, we need to compute the following quantity:

[TABLE]

We start, as before, by summing the integral in the numerator over all the possible cylinders of length 3. Among those cylinders, it is easy to check that only 6 of them contribute to the mass, namely $([3,1,1],[3,2,1],[3,3,1],[3,1,2],[3,2,2]$ and $[3,3,2])$. For the integral associated to the first cylinder, we have (the five others are similar):

[TABLE]

[TABLE]

By summing over the six cylinders and dividing by the denominator, we get the result:

[TABLE]

We tested this result numerically with different sets of target points following the aforementioned assumptions. We find good agreements of $\hat{q}_{0}$ with the theoretical value of $q_{0}=0.09375$ and of $\hat{q}_{1}$ with the theoretical value of $q_{1}=0.0234375$ with a precision of order $10^{-3}$ (see table 4).

In the case where none of the points $z_{i}$ is the $k^{th}$ iterate of another (in particular none of the points $z_{i}$ is $k$-periodic), it can be shown that $q_{k-1}$ is equal to [math] for $k>2$. This is indeed what we find in our numerical simulations, which were performed up to the order $5$.

[TABLE]

This example is in some sense atypical since the four points $z_{0},\dots,z_{3}$ could be far away from each other. For instance, $z_{1}$ and $z_{2}$ are the two predecessors of $z_{0}$ and therefore they are on the opposite sides of $\frac{1}{2}$. By constraining the points to be in a small neighborhood of a given privileged center, the previous effects should be absent and the EI should be one, or very close to it. However our example shows that in the presence of moving target, it is not the periodicity which makes the EI eventually less than one. This gives another concrete example where the series in the spectral formula has at least two terms different from zero.

4.2.2. Continuous noise

We claimed above that when the center of the target ball can take any value in a small neighborhood of a given point $z_{0}$, the EI collapses to $1$, even for periodic points $z_{0}$. To model this more physically realistic situation, let us fix $\epsilon>0$ and consider the set $Z_{\epsilon}=[z_{0}-\epsilon,z_{0}+\epsilon]$. We define a map $f$ acting on $Z_{\epsilon}$ with an associated invariant probability measure $\nu$ that drives the dynamics of the target point $z\in Z_{\epsilon}$. We suppose that $\nu$ is not atomic. The observable considered is now $\phi(x,z)=-\log|x-z|$ on the product space $\{X\times Z_{\epsilon},\mu\times\nu\}$. By similar arguments to the ones we described for the discrete perturbation of the target point, we can show the existence of an extreme value law for the process $\phi\circ(T^{k},f^{k})$, with an EI given by Eq. (7), with

[TABLE]

where $B(y,n)$ denotes a ball around $y$ of radius $e^{-u_{n}}$. We have:

Proposition 3**.**

Suppose that for all $k\in\mathbb{N}$, $\nu(\{z\in Z_{\epsilon}|T^{k+1}(z)=f^{k+1}(z)\})=0$ and that $\mu$ is absolutely continuous with respect to Lebesgue with a bounded density $h$ such that $h\geq\iota>0$. Then the extremal index is $1$.

Proof.

We compute the term $q_{k}$ given by Eq. (28) and denote

[TABLE]

where $r_{n}=e^{-u_{n}}$. We can write the numerator in Eq. (28) as a sum of integrals over $Z_{1}^{n,k}$ and its complementary set $Z_{2}^{n,k}=Z_{\epsilon}\backslash Z_{1}^{n,k}$.

Since for $z\in Z_{1}^{n,k}$, $z$ and $z^{*}$ are at a distance larger than $r_{n},$ for all $z^{*}\in T^{-k-1}B(f^{k+1}(z),n)$, we have that $B(z,n)\cap T^{-k-1}B(f^{k+1}(z),n)=\emptyset$, and the integral over $Z_{1}^{n,k}$ is zero for all $n$.

It remains now to treat the integral over $Z_{2}^{n,k}$. We have that:

[TABLE]

With our assumptions on $h$, the fraction in the last term is finite. Let us consider the limit set $Z_{2}^{k}=\cap_{n}Z_{2}^{n,k}=\{z\in Z_{\epsilon}|T^{k+1}(z)=f^{k+1}(z)\}$. We have that $\lim_{n\to\infty}\nu(Z_{2}^{n,k})=\nu(Z_{2}^{k})=0$, by hypothesis and so $q_{k}=0$ for all $k$ and therefore $\theta=1.$ ∎

Remark 5**.**

In practice, the requirement that $\nu(\{z\in Z_{\epsilon}|T^{k+1}(z)=f^{k+1}(z)\})=0$ for all $k$ is satisfied for a large class of situations. We take two maps $T$ and $f$ intersecting at a countable number of points $Z_{2}^{1}$. We suppose that the preimage of any point by the applications $T$ and $f$ is a countable set. Then the set $Z_{2}^{k}=T^{-k-1}Z_{2}^{1}\cup f^{-k-1}Z_{2}^{1}$ is countable as well for all $k\in\mathbb{N}$. Since $\nu$ is non-atomic, we have that $\nu(\{z\in Z_{\epsilon}|T^{k+1}(z)=f^{k+1}(z)\})=0$ is verified and the proposition holds.

4.2.3. Observational noise

Instead of considering a perturbation on the target set driven by the map $f$, we consider the center of the target ball as a random variable uniformly distributed in the neighborhood of size $\epsilon$ of the point $z_{0}$. This situation models data with uncertainty or disturbances in their detection, and is equivalent to the case observational noise which we considered in [21] in the framework of extreme value theory. In this approach, the dynamics of a point $x$ is given by $T^{k}x+\xi_{k}$, where $\xi_{k}$ is a random variable uniformly distributed in a ball centered at [math] of radius $\epsilon$. The process is therefore given by $-\log(dist(T^{k}x+\xi_{k},z_{0}))$. When this is transposed in the moving target case, it becomes $-\log(dist(T^{k}x,z_{0}+\xi_{k}))$ and the two approaches define the same process, which give the same stationary distribution. In [21] we proved rigorously in some cases and numerically in others, that for a large class of chaotic maps, an extreme value law holds with an extremal index equal to $1$. By equivalence of the two processes, an extreme value law also holds in the moving target scenario and the EI is equal to $1$.

We tested numerically the scenario of a target point following a uniform distribution in an interval of length $\epsilon$ centered in a point $z_{0}$, for several maps of the circle (Gauss map, $3x-\mathrm{mod}1$, rotation of the circle). For all of them, we find that for periodic $z_{0}$, the extremal index which is less than $1$ when unperturbed, converges to $1$ as the noise intensity $\epsilon$ increases (see an example for the $3x-\mathrm{mod}1$ map in figure 4-A). When $z_{0}$ is not periodic (we chose it at random on the circle), the EI is equal to $1$ for all $\epsilon$ (see again figure 4-B).

5. The dynamical extremal index

Up to now the target set was a ball around a point. In the attempt to describe the synchronisation of coupled map lattices, Faranda et al. [15] introduced the neighborhood of the diagonal in the $n$-dimensional hypercube and defined accordingly the extremal index related to the first time the maps on the lattices become close together. Faranda et al. [15] then used that approach to define a new type of observable, first in [20] and successively generalized in [8], which brought to a different interpretation of the extremal index, in terms of Lyapunov exponents instead of periodic behaviors.

Let us consider the $k-$fold ($k>1$) direct product $(X,\mu,T)^{\bigotimes k}$ with the direct product map $T_{k}=T\circ\cdots\circ T$ acting on the product space $X^{k}$ and the product measure $\mu_{k}=\mu\circ\cdots\circ\mu$.

Let us now define the observable on $X^{k}$

[TABLE]

where each $x_{i}\in X$. We also write $\overline{x}_{k}:=(x_{1},x_{2},\dots,x_{k})$ and $T_{k}(\overline{x}_{k})=(Tx_{1},\dots,Tx_{k})$. As we explained in [20], if we put $M_{n}(\overline{x}_{k})=\max\{\phi(\overline{x}_{k}),\dots,\phi(T^{n-1}_{k}(\overline{x}_{k})\}$ and look at the distribution of the maximum $\mu_{k}(M_{n}\leq u_{n})$, this distribution is non-degenerate and converges to the Gumbel law for $n\to\infty$ $e^{-\theta_{k}\tau}$, provided we can find a sequence $u_{n}\to\infty$ and verifying $\mu_{k}(\phi>u_{n})\to\frac{\tau}{n}$, where $\tau$ is a positive number. The quantity $\theta_{k}$ is our new extremal index and it will be described later on. Notice first that

[TABLE]

Our first motivation to investigate the observable in Eq. (30) was the fact that the integral on the right hand side of Eq. (31) scales like $e^{-u_{n}D_{k}(k-1)}$, where $D_{k}$ denotes the generalized dimension of order $k$ of the measure $\mu$. We now describe $\theta_{k}$. We first define:

[TABLE]

By using the spectral technique in [20], and the analytical results of [15], it is possible to show that:

[TABLE]

The quantity $\theta_{k}$ is the dynamical extremal index (DEI) appearing at the exponent of the Gumbel law.

For $C^{2}$ expanding maps of the interval, which preserve an absolutely continuous invariant measure $\mu=hdx$ with strictly positive density $h$ of bounded variation, it is possible to compute the right hand side of Eq. (33) and get:

[TABLE]

All the $k-1$ integrals above factorize, and depend on the parameter $x_{1}$. Therefore they can be treated as in the proof of Proposition 5.5 in [15], yielding the rigorous result:

Proposition 4**.**

Suppose that: the map $T$ belongs to $C^{2}$; it preserves an absolutely continuous invariant measure $\mu=hdx$, with strictly positive density $h$ of bounded variation; it verifies conditions $P1-P5$ and $P8$ in [15]555These conditions essentially ensure that the transfer operator associated with the map $T$ has a spectral gap and that the density $h$ has finite oscillation in the neighborhood of the diagonal.. Then

[TABLE]

This formula uses the translational invariance of the Lebesgue measure: we refer to sections II-B and II-C in [20] for analogous extensions to more general invariant measures and to SRB measures for attractors. As remarked in [20], whenever the density does not vary too much, or alternatively the derivative (or the determinant of the Jacobian in higher dimensions) are almost constant, we expect a scaling of the kind: $\theta_{k}\sim 1-e^{-(k-1)h_{m}},$ where $h_{m}$ is the metric entropy (the sum of the positive Lyapunov exponents).

In section D of [20] and for the case $k=2$, we replaced the iterations of a single map with concatenations of i.i.d. maps chosen with additive noise, see above, section 4.1. Although the dynamical extremal index is related to the Lyapunov exponents, it is also influenced by the fact that the set $\Delta=\cup_{i=2}^{k}\{x_{i}=x_{1}\}$ is invariant in the deterministic situation. By looking at Eq. (33), we see that we estimate the proportion of the neighborhood of the invariant set $\Delta$ returning to itself. As argued in [20], that estimate gives information on the rate of backward volume contraction in the unstable direction. Since the noise generally destructs these invariants sets, we expect the extremal index be equal to $1$ or quickly approaching $1$ for all $k$ when the noise increases. The situation depends again on the characteristics of the noise.

5.1. Additive noise

The prediction at the end of the preceding section is confirmed by concatenating maps perturbed with additive noise and smooth probability density function (see the numerical computations shown in figure 9). As in section 4.1.2 the perturbation is of annealed type. For each map, trajectories of $10^{7}$ points were simulated and the $0.999-$quantile of the distribution of the observable $\phi$ was selected as a threshold. As only the $q_{0}$ term is non-zero, we choose the Süveges estimate for our computations, to be assured that no error coming from $\hat{q}_{k}$ terms of higher order are added.

5.2. Discrete noise

A different scenario happens if we look at sequences of finitely maps chosen in an i.i.d. way according to a Bernoulli process as we did in section 4.1.2. We choose now the two maps $f_{0}(x)=3x-\mathrm{mod}1$ and $f_{1}(x)=3x+b-\mathrm{mod}1$ and we put the same distribution $\mathbb{G}$ defined in 4.1.2. A combination of Eqs. (23) and (33) gives for the term $q_{0}$ entering the infinite sum defining the DEI:

[TABLE]

where: $\xi^{k}$ is the vector $\xi^{k}=(\xi_{1},\dots,\xi_{k})$, $f_{\xi^{k}}^{-1}=f_{\xi_{1}}^{-1}\cdots f_{\xi_{k}}^{-1}$ and $\rho_{k}$ is the stationary measure constructed as in [4]. We now split the integral over the cylinders $\xi^{k,0}=[\xi_{1}=0,\dots,\xi_{k}=0],\xi^{k,1}=[\xi_{1}=1,\dots,\xi_{k}=1]$ and their complement. Therefore we have:

[TABLE]

[TABLE]

The two fractions in Eq. (37) are equal to the $q_{0}^{(\text{unp})}$ of the unperturbed systems, which are the same by the particular choices of the maps and by Eq. (35).

All the other terms in Eq. (38) are zero, we now explain why. Let us in fact consider a vector $\tilde{\xi}^{k}$ different from $(0,0,\dots,0)$ and $(1,1,\dots,1)$. Suppose that the first coordinate $\tilde{\xi}^{k}_{1}$ of $\tilde{\xi}^{k}$ is zero and let $\tilde{\xi}^{k}_{j}\neq 0$ for some $1<j\leq k$. If $x_{k}\in\Delta_{n}^{k}$, and since $f_{0}x_{1}=3x_{1}-\mathrm{mod}1$, $f_{1}x_{j}=3x_{j}+b-\mathrm{mod}1$, we have that $d(f_{0}x_{1},f_{1}x_{i})\rightarrow_{n\to\infty}\min\{b,1-b\}$ as $d(x_{1},x_{i})<e^{-u_{n}}$.

Therefore when $n$ is large enough, $\forall x_{k}\in\Delta_{n}^{k}$, $f^{\tilde{\xi}^{k}}x_{k}:=f^{\tilde{\xi}_{k}^{k}}\cdots f^{\tilde{\xi}_{1}^{k}}x_{k}\notin\Delta_{n}^{k}$ and $\Delta_{n}^{k}\cap f^{-1}_{\tilde{\xi}^{k}}\Delta_{n}^{k}=\emptyset$. In conclusion we have

[TABLE]

where $p$ is the weight of the Bernoulli measure. This last formula has been tested numerically with the maps $f_{0}$ and $f_{1}$, for different values of $b$ and with $p=0.5$. Good agreement is found for different values of $b$ (see figure 8). We use trajectories of $5.10^{7}$ points and the $0.995-$ quantile of the observable distribution as a threshold. Estimations are made using the Süveges estimate, for the same reasons as described earlier.

6. Point processes: statistics of persistence

We said in section 2 that the extremal index is the inverse of the average cluster size when we look at exceedances in the shrinking neighborhood of a given target point $z$. To be more precise let us start with the deterministic case.

6.1. The deterministic case

We consider the following counting function

[TABLE]

where the radius $r_{n}$ goes to [math] when $n$ tends to infinity. We are interested in the distribution

[TABLE]

for $n\rightarrow\infty.$ It has been proved (see for instance [30, 24]) that when $z$ is not a periodic point, $\mu(N_{n}(t)=k)$ converge to the Poisson distribution $\frac{t^{k}e^{-t}}{k!}$, while for a periodic point of minimal period $p$ we have the Polyà-Aeppli distribution

[TABLE]

where $\theta$ is the extremal index.

We now invoke a very recent theory developed by [31], which allows us to put the statistics of the number of visits in a wider context and to make explicit the connection with the EI, which is the main object of this note.

We first replace the ball $B(z,r_{n})$ with a sequence of sets $B_{n}(\Gamma)$ whose intersection over $n\in\mathbb{N}$ gives the zero measure set $\Gamma$. For instance in section 5 such sets are the neighborhoods $\Delta^{k}_{n}$. We keep in mind this example for the next application. We then define the cluster size distribution $\pi_{l}$ as

[TABLE]

Under a few smoothness assumptions for the map $T$, which are verified for the majority of what we called chaotic systems in this note, it is possible to prove ([31], Th. 1) that

[TABLE]

as $n\rightarrow\infty$, where $\nu$ is the compound Poisson distribution for the parameter $t\pi_{l}$. This distribution can be recovered by its generating function $g_{\nu}(z)=\exp{\int_{0}^{\infty}(z^{x}-1)dH(x)}$, where $H$ is the counting measure on $\mathbb{N}$ defined as $H=\sum_{l}t\pi_{l}\delta_{l}$, with $\delta_{l}$ the point mass at $l$. The Polyà-Aeppli distribution is a particular case of compound Poisson distribution with the $\pi_{l}$ having a geometric distribution $\pi_{l}=\theta(1-\theta)^{l}$. The interesting feature of the theory in [31] is that the $\pi_{l}$ can be expressed in terms of more accessible quantities, in particular of the distribution function $\hat{\alpha}_{l}$ of the $(l-1)$ return into the set $B_{n}(\Gamma)$ for $n$ tending to infinity, namely $\pi_{l}=\frac{\hat{\alpha}_{l}-2\hat{\alpha}_{l+1}+\hat{\alpha}_{l+2}}{\hat{\alpha}_{1}-\hat{\alpha}_{2}},$ where

[TABLE]

being $\tau_{B_{n}(\Gamma)}^{\ell}(x)=\tau_{B_{n}(\Gamma)}^{\ell-1}(x)+\tau_{B_{n}(\Gamma)}(T^{\tau_{U}^{\ell-1}}(x)),$ with $\tau_{B_{n}(\Gamma)}^{1}=\tau_{B_{n}(\Gamma)}=\min\{j\geq 1:T^{j}(x)\in B_{n}(\Gamma),x\in B_{n}(\Gamma)\},$ the first return time. Moreover if $\theta$ denotes the extremal index (as defined in section 5 for general target sets), we have the identity:

[TABLE]

This result is a generalization of Eq. (10) stated in section 2 and which targeted neighborhoods of periodic points.

We now apply Eq. (46) to the dynamical extremal index introduced in section 5 in the case of direct product of two maps. For the class of maps verifying Proposition 35 we have the general result for the $\hat{\alpha}_{l}$ [31]:

[TABLE]

This formula was also previously obtained by [11]. In order to make rigorous computations, we consider a Markov map of the interval $T$ for which the density $h$ is piecewise constant (see for instance [28]):

[TABLE]

The density reads

[TABLE]

The extremal index is given by the integral in Eq. (35) which has already been computed in [8] and gives $\theta\approx 0.5926$. The quantities $\hat{\alpha}_{l}$ can also be computed explicitly. We defer to [26] for the details, and obtain:

[TABLE]

The interesting point is that the $\pi_{k}$ do not follow a geometric distribution. For the statistics of the number of visits, we get a compound Poisson distribution which is not Polyà-Aeppli. To illustrate that, we show the distributions $\mu(N_{n}(t))=k$ for $t=50$ and different values of $k$ for the real process associated to the markov map $T$, and the distribution obtained by supposing a Polyà-Aeppli with parameter equal to $0.5926$. In this context, the sets $B_{n}(\Gamma)$ are the tubular neighborhoods of the diagonal $\Delta_{n}^{2}$ defined in Eq. (32). For our numerical simulation, we fix a small set $\Delta_{n}^{2}$ by taking $u_{n}$ in Eq. (32) equal to the $0.99-$quantile of the distribution of the observable computed along a pre-runned trajectory of $10^{4}$ points. We obtain that for a small enough neighborhood of the diagonal, $\mu_{2}(\Delta_{n}^{2})\approx(h_{1}^{2}+h_{2}^{2}+h_{3}^{2})\sqrt{2}\exp{(-u_{n})}$.

6.2. The random case

There are only few results on the statistics of the number of visits in the case of random systems. In the forthcoming paper [23], we show that the quenched Lasota-Yorke example studied in section 3.2.1 gives a Poisson distribution for the visits around almost all target point $z$. A Polyà-Aeppli distribution is seen to hold for some particular random subshits.

We are not aware of a rigorous treatment of these point processes for stationary (annealed) random systems. We believe that the theory developed in [31] could be applied as well to those systems. In particular we expect to get the extremal index as the inverse of the average cluster size. To corroborate this claim, we computed the statistics of the number of visits for two examples: the case in section 4.1.2 for the composition of two maps with $b=1/2$, for which there are infinitely many non zero $q_{k}$, and the moving target example in section 4.2. In both cases we found that the real distribution coincide with Polyà-Aeppli with parameter given by $t=50$ and the extremal index of respectively is $2/3$ and $1-\frac{3}{2}(\frac{1}{4})^{2}-6(\frac{1}{4})^{3}(\frac{1}{2})^{2}\approx 0.8828$. These results are stable against different values of $t$.

6.3. The rotational case

We now compute the statistics of the number of visits for the process described in section 3.2.2. The statistics are computed in the same way as for the preceding examples and the theoretical value of 1 for the extremal index suggests that the number of visits is purely Poissonian. This is indeed what is observed in figure 10. Again the results are stable against different values of $t$.

6.4. Climate data

We study now the statistics of the number of visit for the atmospheric data presented in section 2.1. The target set considered in the analysis is a ball centered at a point $z$ that corresponds to the mean sea-level pressure recorded on Jan 5th 1948. The radius $r=e^{-\tau}$, is fixed using as $\tau$ the $0.95-$quantile of the observable distribution. We considered the number of visits in an interval of length 365 days, that corresponds to the annual cycle666This time scale is relevant for the underlying dynamics of the atmospheric circulation. so that the parameter $t$ in Eq. (40) is taken equal to $t=\lfloor 365\mu(B(z,r))\rfloor$. For the target set considered, we get $t=18$, that is close to the mean value that we find for the number of visits of $18,13$. In figure 11, we compare the obtained empirical distribution with the Polyà-Aeppli distributions of parameters $t$ and the $EI$ obtained using different estimates: Süveges, $\hat{\theta}_{5}$ and $\hat{\theta}_{10}$. The distribution using the Süveges estimate seems to fit the empirical data better than the other two. In fact if we increase the order of the estimate $\hat{\theta}_{m}$, the value it attributes to the EI decreases and the associated Polyà-Aeppli distribution flattens. We believe that this is due to finite effects; indeed we observe that if the fixed threshold $u$ in formula 19 is too small, the numerator approaches 0 as $m$ increases, yielding bad estimates for the EI. Qualitatively similar results are found for different atmospheric states and different quantiles considered (although taking higher quantiles gives too few data for the statistics), suggesting that a Polyà-Aeppli distribution for the number of visits holds for this type of high dimensional observational data. This analysis, applied systematically to a state of a complex system, can provide important information on the recurrence properties of the target state, namely how many visit to expect in a given time interval. In this sense, the extremal index only provides an averaged information. The drawback is in the number of data required to provide a reliable estimates of the visits distribution (see Figure 11). This problem could be however overcome by averaging the distribution of visits for class of events (several states $z$) instead of considering just one state.

7. Conclusions

In this paper we have reviewed the meaning of the extremal index in dynamical systems and provided numerical and theoretical strategies for its estimation. We have detailed how, depending on the nature (stochastic or deterministic) of the underlying system the estimators may provide different insights on the dynamics of the system. Our main findings can be summarized as follows:

• •

As soon as a deterministic system is perturbed with noise having a smooth distribution, the EI becomes equal to $1$. With discrete (point masses) distributions, the EI could be $1$ for almost choices of the target and of the realisations, (see sections 3.2.1, 3.2.2, 5.1), but it could also be less than $1$, (see sections 4.1.2, 5.2).

• •

What we described in the previous item happens either with annealed and quenched stochastic perturbations. In the annealed case with discrete probability distributions, the value of the EI could depend either on the choice of the masses and on the closeness of the maps (see Remark 4). This suggests that in the situations where the EI is generically $1$, if one finds a lesser value, this could be the effect of a discrete perturbation.

• •

The situation where the dynamics is deterministic, but the target set is known with some probability with a continuous distribution function is interesting, as discussed in section 4.2.2. This describes physical situations where the localisation and observation of the extreme event are perturbed. No matter of the intensity of those perturbations, the extremal index is everywhere equal to one.

• •

The statistics of the number of visits (persistence) in the shrinking target set gives rise to point processes of Poisson type. We have shown that the geometric nature of the target set produces different compound Poisson distributions. Moreover these compound Poisson distributions persist even when the system is randomly perturbed. Indeed the fact that the extremal index is the inverse of the average cluster size seems to be a robust property of chaotic systems in the deterministic and random settings.

• •

In presence of periodicity one expects asymptotic distributions of Polyà-Aeppli type. Nevertheless for target sets composed by more general invariant sets, different compound Poisson distribution could emerge up as we showed above. This suggest that for the statistics of the number of visits around sets with a more complicated geometry, and this could be relevant in applications, the computation of the EI is probably not enough, and one should go to higher moments to get the real asymptotic distribution.

8. appendix

As mentioned in section 3, results proving the existence of extreme value laws for non-stationary sequential systems are lacking. Nevertheless, the methods of estimation described in the paper can still be used to evaluate a quantity $\theta$ that would correspond to the extremal index in the eventuality that an extreme value law holds for this kind of systems, which by the way was the case in [25]. We could also alternatively define the extremal index by computing the statistics of the number of visit and evaluating the expectation in Eq. (46). We choose this second option and consider the motion given by the concatenation

[TABLE]

where the probability law of the $\xi_{i}$ changes over time. In particular, we consider the 10 maps

[TABLE]

where $b_{i}$ is the $i^{th}$ component of a vector $\bar{b}$ of size 10, with entries equally spaced between [math] and $1/2$. We consider sequences of time intervals $[\tau k+1,(k+1)\tau]$, with $\tau=10$, for $k=0,1,2,\dots$ in which the weights associated to $\xi_{i}$ are equal to $p_{i}^{k}$. For every $\tau$ iterations, the weight associated to each $\xi_{i}$ changes randomly, with the only constraint that they sum to 1. $\hat{\theta}_{5}$ is computed considering a trajectory of $5.10^{7}$ points, and a threshold $u$ corresponding to the $0.995-$quantile of the observable distribution. Using the same trajectory, we computed the empirical distribution of the number of visits in a ball centered at the origin and of radius $r=e^{-u}$ in intervals of time of length $\lfloor 2tr\rfloor$, with $t=50$. This is indeed $t$ times the Lebesgue measure (which is the invariant measure associated to all the maps $f_{i}$) of a ball of radius $r$ centered in the origin. Choosing the same trajectory for the computation of $\theta$ and for the statistics of visits is of crucial importance here, because different probability laws imply large variations of $\theta$ depending on the trajectory considered. In figure 12, we observe again a perfect agreement between the empirical distribution of the number of visits and the Polyà-Aeppli distribution of parameters $t=50$ and $\hat{\theta}_{5}$, which is equal to 0.92 for the trajectory presented in the figure 12. Although different trajectories give variations for $\theta$, this agreement is stable against different trajectories and different values of $t$ and $\tau$, suggesting that a geometric Poisson distribution for the number of visits given by $\theta$ is a universal feature, even for non-stationary scenarii.

9. Acknowledgement

We would like to thank Jorge Freitas, Michele Gianfelice, Nicolai Haydn and Giorgio Mantica for several discussions related to different parts of this work. PY was supported by ERC grant No. 338965-A2C2.