Extreme value theory for strictly stationary sequences has been extensively studied, initiated in the works of Watson (1954), Berman (1964), Loynes (1965), and continued by Leadbetter (1974, 1983) and O’Brien (1987) amongst others. One of the key findings in this line of research is that unlike in independent and identically distributed sequences where extreme values tend to occur in isolation, stationary sequences possess an intrinsic potential for clustering of extremes, i.e., several extreme values may be observed in quick succession. Understanding the extremal clustering characteristics of a random process is critical in many applications where a run of extreme values may have serious consequences. For example, if a sequence consists of daily temperatures at some fixed location then a cluster of extremes may correspond to a heatwave.
The extent to which extremal clustering may occur is naturally measured for stationary sequences by a parameter known as the extremal index. Suppose that the stationary sequence of random variables has marginal distribution function and let . In the special case where the random variables in the sequence are independent then, for any , , whereas in general, provided a suitable long range dependence restriction is satisfied then, for large , . Here, and in what follows, is a sequence of real numbers converging to in such a way that as .
Leadbetter (1983) showed that is the limiting mean cluster size, and Hsing et al. (1988) further developed this result by considering the point process of exceedance times above a large threshold. In particular, Hsing et al. (1988) found that the point process of normalized exceedance times converges weakly as to a compound Poisson process on . The compound nature of the Poisson process arises from the fact that the cluster positions follow a homogeneous Poisson process, each of which is marked by an independent realization from the so-called cluster size distribution, the mean of which is . The cluster size distribution is defined by , for , where
and . The stated result about the limiting mean cluster size can then be expressed as In a more general non-stationary setting, one may expect the cluster size distribution to vary with time. Although we will not consider a point process approach here, the results presented in the following sections implicitly confirm such an expectation.
Another characterization of that links it to the extremal clustering properties of a stationary sequence can be found in O’Brien (1987). Defining , O’Brien (1987) found that the distribution function of satisfies
for , and provided the limit exists, as . This result clarifies that smaller values of are indicative of a larger degree of extremal clustering, since the conditional probability in (
are indicative of a larger degree of extremal clustering, since the conditional probability in (1.3) is small when an exceedance of a large threshold is likely to soon be followed by another exceedance.
Early attempts at estimating were based on associating with the limiting mean cluster size. Different methods for identifying clusters gave rise to to different estimators, well known examples being the runs and blocks estimators (Smith and Weissman, 1994). For the runs estimator, a cluster is identified as being initialized when a large threshold is exceeded and ends when a fixed number, known as the run length, of non-exceedances occur. The extremal index is then estimated by the ratio of the number of identified clusters to the total number of exceedances. A difficulty that arises when using this estimator is its sensitivity to the choice of run length (Hsing, 1991).
The problem of cluster identification was carefully studied by Ferro and Segers (2003) who considered the distribution of the time between two exceedances of a large threshold. They found that the limiting distribution of appropriately normalized interexceedance times converges to a distribution that is indexed by . In particular, for a given threshold , they define the random variable , and find that as converges in distribution to a mixture of a point mass at zero and an exponential distribution with mean
converges in distribution to a mixture of a point mass at zero and an exponential distribution with mean. Thus, by computing theoretical moments of this limiting distribution and comparing them with their empirical counterparts, they construct their so-called intervals estimator.
Motivated by the fact that many real world processes are non-stationary, in this paper we are led to investigate the effect on the extremal clustering properties of a random sequence when the stationarity assumption is dropped. Previous statistical works that invoke the concept of extremal clustering in a non-stationary sequence include Süveges (2007) and Coles et al. (1994). Süveges (2007) uses the likelihood function introduced by Ferro and Segers (2003) for the extremal index together with smoothing methods to capture non-stationarity in a time series of temperature measurements, whereas in a similar application, Coles et al. (1994) use a Markov model together with simulation techniques to estimate the extremal index within different months.
use a Markov model together with simulation techniques to estimate the extremal index within different months.
Early works that develop the extreme value theory for non-stationary sequences with emphasis on the asymptotic distribution of sample maxima include Hüsler (1983, 1986). Hüsler (1983) focuses on the case of sequences with a common marginal distribution and Hüsler (1986) considers the more general case where the margins may differ. While neither papers discuss extremal clustering or statistical inference, Hüsler (1986) discusses the difficulties of generalising the extremal index to non-stationary sequences.
In this paper we see that the notion of a time varying extremal index may be defined for a class of non-stationary sequences in a way that naturally extends the theory for stationary sequences. In the simple case of an identically distributed but not necessarily stationary sequence with common marginal distribution function , we find that under similar assumptions as in O’Brien (1987) that
Thus, we find in this setting, that the limiting distribution of the sample maximum at large thresholds is characterized by a parameter , (provided the limit exists). By analogy with equation (1.3), we see from (1.5) that may be regarded as the average of local extremal indices, each of which measures the potential for extremal clustering at a particular location in the sequence. We will see that these local extremal indices can be estimated using the intervals estimator of Ferro and Segers (2003) which is a function of the observed times between consecutive exceedances of a high threshold. In the case where the sequence is stationary, so that all terms in the summation (1.5) are equal, the formula for reduces to in (1.3).
The result stated in equations (1.4) and (1.5) is reminiscent of the main result in de Haan (2015). In de Haan (2015) , independent but non-identically distributed random variables
with proportional tails are considered, a regime known as heteroscedastic extremes (see also
, independent but non-identically distributed random variables with proportional tails are considered, a regime known as heteroscedastic extremes (see alsoEinmahl et al. (2016)). In such a setting, de Haan (2015) finds that the asymptotic distribution of appropriately normalized sample maxima is characterized by the Cesàro mean of the tail proportionality constants. We deduce a more general version of this result in Section 5 where we allow for dependence in the sequence.
The structure of the paper is as follows. Section 2 defines the notation and assumed mixing condition used throughout the paper and states the main theoretical results regarding the asymptotic distribution of the sample maxima and normalized interexceedance times in the case of a common marginal distribution. Section 3 discusses approaches to parameter estimation using the result from Section 2 on the distribution of the interexceedance times. Section 4 considers the estimation problem for two simple non-stationary Markov sequences with periodic dependence structures. Section 5 discusses extensions of the main results to sequences with different marginal distributions and Section 6 gives the proofs of the main theoretical results.
2 Theoretical results
2.1 Notation and definitions
Unless otherwise stated, we will assume that all random variables in the sequence have common marginal distribution with upper endpoint , though we do not assume stationarity. We discuss the consequence of dropping the common marginal assumption in Section 5. We write . In addition to the definitions for and given in the introduction, we define where is an arbitrary set of positive integers. We denote the number of elements in the set by . We also refer to a set of consecutive integers as an interval. We denote asymptotic equivalence of the real valued functions and by , which means that . For , the -algebra generated by the events is denoted by where
In all results in this section, we use the standard technique of block-clipping, see for example Section 10.2.1 in Beirlant et al. (2004), to divide observations up in to sets of alternating sizes of and where and . Specifically, we define intervals and by
for , where .
We assume the asymptotic independence (AIM) mixing condition of O’Brien (1987) which restricts long range dependence and gives approximate independence of the random variables and , .
(AIM()). The sequence is said to have asymptotic independence of maxima relative to the sequence (AIM()) if there exists a sequence of positive integers with such that for any two intervals and separated by we have
where the maximum is taken over all positive integers and such that , and .
Remark: In the following results we abbreviate the mixing coefficents in Definition 2.1 to .
2.2 Asymptotic distribution of
Before giving the main result on the asymptotic distribution of , we restate Lemma 3.1. from O’Brien (1987) amended to allow for non-stationarity.
Remarks: If we can find a sequence such that the conditions in Definition 2.1 hold, then we can automatically find a sequence such that (2.2) holds, for example, by taking . Moreover, if is a sequence of positive integers such that and , then Lemma 2.1 holds with replaced by and replaced by In the special case where the sequence is stationary, so that all terms in the product in (2.5) are equal, we have that . To see this, suppose that as Then, and similarly and so the stated result is a consequence of the fact that . By the stationarity of the sequence, it follows that if and are any two subsets of of lengths and respectively, then In Theorem 2.1 below, we will require a weak version of this to hold in the non-stationary setting, and this is stated in equation (2.6). This condition, loosely speaking, says that we are much more likely to observe an exceedance of a large threshold within a large set of observations than within a small set. One way that such a condition may fail to hold is if the dependence in the sequence is too strong. For example, in the contrived case of complete dependence where for we have that regardless of the sizes of and . Of course, such a situation is excluded by the mixing condition in Definition 2.1, however it would seem that Definition 2.1 does not imply (2.6) in the general non-stationary case. An obvious sufficient condition for (2.6) to be satisfied is that the terms in the product of (2.5) are equal, or at least asymptotically equivalent, so that the same argument as in the stationary case holds.
We may construct non-stationary sequences with periodicity in the dependence structure that satisfy this property. For example, suppose that and for , let where independent of and is periodic sequence, i.e., there is an such that for all . In such a case, and are identically distributed when is a multiple of . Non-stationary Markov models of this type are considered further in Section 4.
We can now state the main result.
2.3 Interexceedance times
Ferro and Segers (2003) provided a method for estimating the extremal index of a stationary sequence without the need for identifying independent clusters of extremes. This was achieved by considering the distribution of the time between two exceedances of a threshold . Specifically, for a fixed they define the random variable by
and investigate the distribution of as approaches . A simple normalization of gives rise to a non-degenerate limiting distribution. In particular, Ferro and Segers (2003) show that
as , converges in distribution to a mixture of a point mass at zero (with probability ) and an exponential random variable with mean (with probability ). The mixture arises from the fact that the interexceedance times can be classified in to two categories: within cluster and between cluster times. The mass at zero stems from the fact that the within cluster times, which tend to be small relative to the between cluster times, are dominated by the factor
). The mixture arises from the fact that the interexceedance times can be classified in to two categories: within cluster and between cluster times. The mass at zero stems from the fact that the within cluster times, which tend to be small relative to the between cluster times, are dominated by the factorwhich goes to zero.
In the stationary case, conditioning on the event in equation (2.10) defining may be replaced with for any (then must be replaced by ) without affecting the distribution of . In the non-stationary case we must allow for the possibility that the distribution of interexceedance times may vary with location in the sequence. Thus we consider for each and threshold , the random variable defined by
We find that the distribution of converges as to a mixture of a mass at zero (with probability ) and an exponential random variable with mean (with probability where is the extremal index at time point (defined below). As in Ferro and Segers (2003), a slightly stronger mixing condition is required to derive this result than was needed for Theorem 2.1. We define the mixing coefficients
where the supremum is over all with and .
We also use the notation and . The result on the limiting distribution of the interexceedance times is now given.
Let be a sequence of random variables with common marginal distribution and suppose that . Suppose (2.6) holds and and as . If there is a sequence of positive integers such that as for all then for
In this section we consider moment and likelihood estimators for the and based on the limiting distribution of normalized interexceedance times given in Theorem 2.2. We first show that the intervals estimator of Ferro and Segers (2003) may be used to estimate the and then consider likelihood based estimation along the lines of Süveges (2007). For simplicity, in the likelihood formulation, we shall assume a finite number of parameters to infer, say with
3.2 Moment based estimators
The limiting result in Theorem 2.2 implies that the first two moments of satisfy and as . If we assume that the threshold is chosen to be suitably large so that the terms can be dropped, we can solve these two equations to get the unknown parameters in terms of the moments of the normalized interexceedance times as
A complication that arises in the non-stationary settting is that, since is defined via a probability conditional on the event if does not exceed the threshold then we will have no interexceedance times relevant to estimating . This problem doesn’t arise in the stationary case where every interexceedance time may be used to estimate the extremal index .
In order to estimate the , it is natural to assume that they are structured in some way, e.g., periodic or piecewise constant. Such a situation may arise, for example, through a periodicity in the dependence structure of the sequence. In this way, we obtain a set of interexceedance times that can be used for estimating Equation suggests the estimator
Note that, while equation (3.1) also suggests an estimator for , this is based only on the interexceedances relevant to estimating and also requires an estimate of . One possibility is to obtain such estimates and take the mean of these as the estimate of . However, this estimator need not respect the relation , a consequence of the fact that we dropped the terms when solving the first two moment equations for for and . In the examples that we consider later, we estimate using the mean of the estimates for the parameters.
which motivates consideration of the positive integer valued random variable defined by
where and and we may identify with In a similar manner to Ferro and Segers (2003), we find that and , so that upon simplification we find that
A Taylor expansion of the right hand side of equation (3.5) about gives
so that the first order bias of is .
On the other hand, since
this motivates the estimator
whose first order bias in zero. This estimator forms the key component of the intervals estimator of Ferro and Segers (2003), which we can use to estimate .
3.3 Maximum likelihood estimation
Theorem 2.2 also allows for the construction of the likelihood function for the vector of unknown parameters. This is an attractive approach due to the modelling possibilities that become available, e.g., parameterizing , smoothing spline interpolation methods or the use of Bernstein polynomials.
However, as discussed in
also allows for the construction of the likelihood function for the vector of unknown parameters. This is an attractive approach due to the modelling possibilities that become available, e.g., parameterizingas a function of
, smoothing spline interpolation methods or the use of Bernstein polynomials. However, as discussed inFerro and Segers (2003) in the stationary case, problems arise with maximum likelihood estimation due to the uncertainty in how to assign the interexceedance times to which component of the limiting mixture distribution. Since the asymptotically valid likelihood is used as an approximation at some subaymptotic threshold , all observed normalized interexceedance times are strictly positive. Assigning all interexceedance times to the exponential part of the limiting mixture means that they are all being classified as between cluster times. This is tantamount to exceedances of a large threshold occuring in isolation, and so the maximum likelihood estimator based on this (typically misspecified) likelihood converges to 1 in distribution regardless of the true underlying value of .
A partial solution to this problem was proposed in Süveges (2007), for a certain class of sequences having the property that the limiting extremal clusters consist only of runs of consecutive exceedances. This contrasts with the more general situation where within clusters, some observations may fall below the threshold. For a stationary sequence, the technical condition required to be satisfied for the model of consecutive exceedances of the clusters is due to Chernick et al. (1991). This condition states that , where , i.e., the probability of again exceeding the threshold after dropping below it within a cluster falls to zero (faster that ) as . For such a sequence, the conditional probability characterization of in O’Brien (1987) reduces to . A diagnostic plot was suggested by Süveges (2007) to assess the plausibility of this assumption from a given realization of a stationary process and it is found that maximum likelihood estimation outperforms the intervals estimator, in the sense of lower root mean squared error, for sequences satisfying this assumption.
If we were to make a similar assumption in the non-stationary case, e.g., for each , so that the consecutive exceedances model for clusters is accurate, we obtain the likelihood function as
where are the interexceedance times for estimating . The full log-likelihood is then
where and in practice must be replaced with an estimate. Unlike in the stationary case, the likelihood equations don’t have a closed form solution, essentially due to the dependence of on all the , however equation (3.11) is easily optimized numerically provided is not too large. If is large, it is more natural to parameterize in terms of a small number of parameters which we may estimate by maximum likelihood or consider non-parametric estimation along the lines of Einmahl et al. (2016).
We may generalise this idea and assign all interexceedance times below some value to the zero part of the mixture so that the corresponding expression for becomes
In the stationary setting the technical condition required is as in which case (Chernick et al., 1991). Selection of an appropriate value of is equivalent to the selection of the run length for the runs estimator, and this problem is considered in Süveges and Davison (2010). However, in a non-stationary setting, where the clustering characteristics of the sequence may change in time, the appropriate value of may also be time varying, so that should be replaced with in equation (3.12). If too small a value is selected for then some of the interexceedance times may be wrongly assigned to the exponential component of the likelihood leading to an overestimate of whereas if is selected to be too large then we tend to underestimate .
In this section we consider two simple examples of non-stationary Markov sequences with a periodic dependence structure and common marginal distributions. In the first example, no limiting extremal clustering occurs at any position in the sequence, so that for each , in contrast to the second example where for each .
We note that for sufficiently well behaved stationary Markov processes, mixing conditions much stronger than those considered in Section 2 hold (Athreya and Pantula, 1986). For example, for the autoregressive sequence , , with independent and identically distributed as , Theorems 1 and 2 from Athreya and Pantula (1986) give that that the mixing conditions of Section 2 hold for any sequence such that , , for any . One would naturally expect that making small changes to the dependence structure of a stationary Markov process, such as letting in the previous example vary along the sequence, would not significantly alter the mixing properties of the sequence. For the two examples that follow this is indeed the case, see for example Bradley (2005) Theorem 3.3 and Davydov (1973) Theorem 4.
4.2 Normal autoregressive model
Stationary sequences where each is a standard normal random variable, are extensively studied in Chapter 4 of Leadbetter et al. (1983). It is shown there that if the lag autocorrelation satisfies as , a condition originally due to Berman (1964), then the extremal index of the sequence equals one and thus no limiting extremal clustering occurs. For example, the autoregressive process , where , has extremal index one provided since the lag correlation is . This is a special case of a more general result that says that a stationary asymptotically independent Markov sequence has an extremal index of one (Smith, 1992). We say that the sequence with common marginal distribution function is asymptotically independent at lag if as , and asymptotically independent if it is asymptotically independent at all lags.
Non-stationary normal sequences are considered in Chapter 6 of Leadbetter et al. (1983) , although the emphasis is on the limiting distribution of sample maximum rather than clustering.
Here, we consider a simple non-stationary autoregressive model is the standard normal distribution function, and thus conclude that
, although the emphasis is on the limiting distribution of sample maximum rather than clustering. Here, we consider a simple non-stationary autoregressive modelwith so that for each , and specify a periodic lag one correlation function for . We note that as and the sequence is asymptotically independent. Applying Theorem 6.3.4 of Leadbetter et al. (1983), and comparing the non-stationary sequence to an independent standard normal sequence, we deduce that as where
is the standard normal distribution function, and thus conclude thatand for
We simulated 1000 realizations of this sequence of length and for each realisation, estimated and for a range of high thresholds, using both the intervals estimator and maximum likelihood with in equation (3.12) equal to zero and one. We then repeated this procedure for sequences of length . We found that maximum likelihood with performed best as measured by root mean squared error in , at all but one combination of threshold and sample size. The one exception was for a sample size of and the threshold equal to the 0.99 theoretical
quantile, for which the expected total number of exceedances is 100, giving approximately 14 interexceedance times to estimate each
and the threshold equal to the 0.99 theoretical quantile, for which the expected total number of exceedances is 100, giving approximately 14 interexceedance times to estimate each. However, in this case the root mean squared error for when using the intervals estimator and maximum likelihood with only differed in the third decimal place. Maximum likelihood with typically performed slighty poorer than the intervals estimator. This is not suprising as we will assign several interexceedance times to the zero component of the likelihood whereas asymptotically, they all belong to the exponential component.
Table 1 shows the 0.025 and 0.975 quantiles for the sampling distributions of the maximum likelhood estimators with In the table, corresponds to the threshold that there is probability of exceeding at each time point i.e., . We note that although the true values of each is 1, so that no extremal clustering occurs in the limit as , clustering may occur at subasymptotic levels. Moreover, there will tend to be more subasymptotic clustering in the sequence at positions with a corresponding larger lag one autocorrelation (larger ). This point has been thoroughly discussed in the context of stationary sequences and estimation of the extremal index, Ancona-Navarrete and Tawn (2000), Eastoe and Tawn (2012), and leads to the notion of a subasymptotic or threshold based extremal index. In the non-stationary setting we may define a threshold based version of . Specifically, for a fixed threshold and positive integer , we define
which provides a measure of the potential for clustering above the threshold at position . The probability in equation (4.1 ) can be evaluated by integration of the
joint distribution of
) can be evaluated by integration of the joint distribution of, which in this case is multivariate normal, normalised by . We performed this calculation for equal to and , for various combinations of and using the technique for evaluating multivariate normal probabilities given in Genz (1992) and implemented in the R package mvtnorm. The results are presented in Table 2. The reason for considering equal to 3 and 6 is that these positions have respectively the highest and lowest associated lag one autocorrelations, namely and , and consequently correspond to the positions where there is greatest and least subasymptotic clustering potential. This feature is clearly seen in Table 2 which shows that for each combination of and .
Another feature of the subasymptotic extremal indices that is clear from Table 2, is that there is less extremal clustering as the threshold increases, i.e, when , so that is a monotone function of for fixed . In fact, for the normal autoregressive model we have, for fixed , as for each and we can use results from Ledford and Tawn (1996, 1997) to make this more precise. First, note that
from which it follows
By the asymptotic independence of the sequence, all conditional probabilites in (4.3) tend to zero as , which verifies that as for fixed . Applying the results of Ledford and Tawn (1996, 1997) for the joint tails of bivariate normal random vectors, we obtain how the conditional probabilities in (4.3) tend toward zero as
where and is the correlation between and . Consequently, we have for ,