Extremes of Dependent and/or Nonstationary Sequences

Extremes of Dependent and/or Nonstationary Sequences

 Much of the theory applied thus far assumes independence of the data, which may not be the case when looking at extreme values because of the tendency for extreme conditions to persist over several observations. The most natural generalization of a sequence of independent random variables is to a stationary series, which is realistic for many physical processes. Here the variables may be mutually dependent, but the stochastic properties are homogeneous over time (see Coles (2001) (b) Ch. 5). Extreme value theory still holds, without any modification, for a wide class of stationary processes; for example, for a Gaussian autoregressive moving average process. With modification, the theory can be extended to an even broader class of stationary processes.

Parameter Variation

 It is possible to allow parameters of the extreme value distributions to vary as a function of time or other covariates. In doing so, it is possible to account for some nonstationarity sequences. One could, for example, allow the location parameter, , of the GEV(, , ) distribution to vary cyclically with time by replacing by (t) = + sin(2 t/365.25) + 2cos(2 t/365.25). When allowing the scale parameter to vary, it is important to ensure that (t)0, for all t. Often a link function that only yields positive output is employed. The {\em log} link function is available for this purpose as an option with extRemes. For example, the model (x) = exp( + x) can be employed using the default linear representation log(x)= + x by checking the appropriate Link button. While it is also possible to allow the shape parameter to vary, it is generally difficult to estimate this parameter with precision; so it is unrealistic to allow this parameter to vary as a smooth function. One alternative is to allow it to vary on a larger scale (e.g., fit a different distribution for each season) if enough data are available (see, for example, Coles (2001) (b) section 6.1).
Example 2: Fort Collins Precipitation (annual cycle)
It is also possible to include a seasonal trend in the model; either within the model parameters or within the threshold. Here, we shall include an annual cycle in the scale parameter. To do this, we first need to create a few new columns in the data. First, we require an indicator variable that is 1 whenever the precipitation exceeds 0.395 inches, and 0 otherwise. Using extRemes: There should now be a new column called Prec.ind0.395 in the Fort Collins precipitation data matrix, Fort$data. Next, we need to add columns that will account for annual cycles. Specifically, we want to add columns that give sin(2 t/365.25) and cos(2 t/365.25), where t is simply the obs column found in Fort$data (i.e., t = 1,...,36524). Using extRemes: There should now be two new columns in Fort$data with the names obs.sin3659 and obs.cos3659 Now, we are ready to incorporate a seasonal cycle into some of the parameters of the Poisson-GP model for the Fort Collins precipitation data. We begin by fitting the Poisson rate parameter () as a function of time. Specifically, we want to find

log(t) = 0 + 1sin(2 t/365.25) + 2cos(2 t/365.25) = 0 + 1. obs.sin365 + 2. obs.cos365. (7.1)


Results from fitting the Poisson rate parameter with an annual cycle (Eq. (7.1)) are 0 -3.72 (0.037), 1 0.22 (0.046) and 2 -0.85 (0.049). Note also that the likelihood-ratio against the null model (Example 1 above) is about 355 with associated p-value 0, which indicates that the addition of an annual cycle is significant. Next, we fit the GPD with the same annual cycle as a covariate in the scale parameter. That is, the scale parameter is modeled by

log(t)= + sin(2 t/365.25) + 2cos(2 t/365.25). (7.2)

MLE parameter estimates for the scale parameter from Eq. (7.2) are 0 -1.24 (0.053), 1 0.09 (0.048) and 2 -0.30 (0.069), and for the shape parameter 0.18 (0.037). The negative log-likelihood value is about 73, and the likelihood-ratio test between this fit and that of section~\ref{subsec:FitGPD} Example 2 is about 24 (associated p-value nearly zero) indicating that inclusion of the annual cycle is significant.

Back to Top

Nonconstant Thresholds

 In addition to varying parameters of the GPD to account for dependencies, it is also possible to vary the threshold. For some, such as engineers, interest may be only in the absolute maximum event, but others, such as climatologists, may be interested in modeling exceedances not only of the absolute maximum, but also in exceedances during a lower point in the cycle.
Example: Fort Collins Precipitation Data
As in example 1 of this section, it will be necessary to create a vector from the R prompt that will be used as the nonconstant threshold. There are many ways to decide upon a threshold for these data. One could have a single threshold, similar to example 1, or one might use a trigonometric function to vary the threshold for each month. The latter will be employed here.
mths Fort$data[,"month"]
u.fortcollins 0.475+5*(-0.03*cos(2*pi*mths/12)) Figure 7.1 shows a plot of the Fort Collins precipitation data with both the previously used constant threshold of 0.4 inches and the above cyclical threshold. The following R commands created the plot in Figure 7.1.
prec Fort$data[,"Prec"]
plot( mths, Fort$data[,"Prec"], xlab="Month", ylab="precipitation (inches)", xaxt="n")
axis(1, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), at=1:12)
abline( h=0.4)
lines( mths[order(mths)], u.fortcollins[order(mths)], col="blue") Fitting data to a point process model using u.fortcollins to fit a nonconstant (seasonal) threshold gives parameter estimates: 1.40 inches (0.043 inches), 0.53 inches (0.034 inches) and 0.16 (0.040); and associated negative log-likelihood of about -619.64. The ideal model would be based on a nonconstant threshold, but it is also possible to include annual cycles in the parameters; compare estimates to those found when including a seasonal cycle in the scale parameter from section PP. Inspection of the diagnostic plots (Figure 7.2) suggests that the model assumptions seem reasonable. For different cycles in the threshold with higher peaks in the summer months resulted in rather poor fits suggesting that too much data is lost, so the lower thresholds are necessary.


Figure 7.1: Fort Collins, C.O. precipitation data with constant threshold of 0.4 inches (solid black line) and nonconstant (cyclic) threshold (solid blue line). Note that although the varying threshold appears to vary smoothly on a daily basis, the threshold used in the example is constant for each month.



Figure 7.2: Probability and quantile plots for fitting data to a point process model to the Fort Collins, C.O. precipitation (inches) data with a seasonal cycle incorporated into the threshold.


Back to Top

Declustering

 Clustering of extremes can introduce dependence in the data that subsequently invalidates the log-likelihood associated with the GPD for independent data. The most widely adopted method for dealing with this problem is declustering, which filters the dependent observations to obtain a set of threshold excesses that are approximately independent. Specifically, some empirical rule is used to define clusters of exceedances, maximums within each cluster are identified and cluster maxima are fit to the GPD; assuming independence among cluster maxima. One simple way to determine clusters is commonly known as runs declustering. First, specify a threshold and define clusters to be wherever there are consecutive exceedances of this threshold. Once a certain number of observations, the run length, call it r, falls below the threshold, the cluster is terminated. There are issues regarding how large both the threshold and r should be, and improper choices can lead to either bias or large variance. Therefore, the sensitivity of results should be checked for different choices of threshold and r. See Coles (2001) (b) Ch. 5 for more on this method and Ch. 9 for some alternatives to declustering. extRemes provides for declustering the data using runs declustering, but in practice declustering is a more involved process that should be executed by the user, and is not supported by extRemes itself. The general procedure for declustering data with the toolkit is as follows.
Example: Phoenix Minimum Temperature
 To decluster the Phoenix minimum temperature (see Poisson-GP Example 3) data using the toolkit (runs declustering), do the following. It is also possible to plot the data with vertical lines at the cluster breaks by clicking on the Plot data checkbox. Here, however, (as is often the case) the amount of data and relatively large number of clusters creates a messy, illegible plot. Therefore, leave this box unchecked for this example. A message will be displayed on the main toolkit window that 84 clusters were found and that the declustered data were assigned to MinT.neg.u-70r1dcbyYear. This column has been added to the original data matrix using this name (where u-70 corresponds to the threshold of -70 and r1 corresponds to r being 1). Other information given includes two estimates of the extremal index. The first estimate is a simple estimate that is calculated after declustering is performed; referred to in the display as being estimated from runs declustering. Namely, the estimate is = nc/N, where nc is the estimated number of clusters and N is the total number of exceedances over the threshold, u. The second estimate is more complicated, but is made prior to declustering the data, and is called the intervals estimator (Ferro and Segers (2003)). Please see appendix for the definition of this estimate. Other information given in the main toolkit dialog is a suggested run length based on the procedure of Ferro and Segers (2003) of r=11, but this number should be disregarded here because we are declustering by year. The procedure for determining the "best" run length employed with this software does not account for covariates when declustering. It is important to decluster by year here because we do not want values from August of one year to be clustered with values from July of the following year. If it were determined unnecessary to decluster by year, then r=11 would still apply for declustering without taking into account the year. Note that because this process reduces the number of data points, values below the threshold have been "filled in" so that the declustered data will have the correct dimensions in order to be added to the original data matrix. Specifically, every point not found to be a cluster maxima is converted to be the minimum of the data and the threshold--i.e., min(x,u). These filled-in values will not affect any POT analyses (using the same or higher threshold) because they are less than the threshold, and subsequently discarded. The original positions of the cluster maxima are preserved so that any covariates will not require further transformations. The optional use of the Decluster by feature ensures that, in this case, values from one year will not be clustered with values from another year. The next step is to fit the declustered data to a GPD. One detail to be carefull about, in general, is that the number of points per year (npy) may be different once the data have been declustered. This will not affect parameter estimates for the GPD, but can affect subsequent calculations such as return levels, which are usually expressed on an annual scale. See Coles (2001) (b) Ch. 5 for an adjustment to the return level that accounts for the extremal index. Results of fitting the GPD to these data are shown in Table 7.1. It is difficult to compare the models using the log-likelihoods here, but there does not appear to be much variability in parameter estimates from one model to the other suggesting that declustering is not important for these data. Particularly the 100-year return level estimates. Each estimate is within the 95\% confidence bounds of every other estimate.

Declustering None r = 1 r = 2 r = 5
3.91 (0.303)4.16 (0.501)4.21 (0.540)4.42 (0.620)
-0.25 (0.049) -0.24 (0.079) -0.24 (0.086) -0.25 (0.097)
100-yr r.l. 59.20 58.67 58.45 58.39
1 0.57 0.53 0.44
- 0.20 0.20 0.20
Table 7.1: Results of fitting data to the GPD to (negative) minimum temperature (degrees Fahrenheit) using a threshold of 73 degrees at Phoenix Sky Harbor airport with: no declustering, runs declustering with run length r=1 (150 clusters), runs declustering with r=2 (138 clusters) and runs declustering with r=5 (115 clusters). Here, is the extremal index estimated by (nc)/n, and is the extremal index estimated as in Ferro and Segers (2003).


These minimum temperature data for Phoenix, A.Z. (Figure 6.3), clearly have an upward trend over time, and possibly a varying standard deviation from one year to the next. The blue line in Figure 6.3 is the least squares fit (using all data, and not just those points above a threshold) t = 0 + 1.(t-1948), which has a significant positive slope. It is, therefore, of interest to also investigate incorporation of a trend into the GPD model. We will also make use here of the Poisson-GP model (or see appendix for more details). The fitted values for the Poisson rate parameter model log(t) = 0 + 1t, with t=1,...,43 (Year - 47) are:

0 -2.38 (0.160)

1 -0.04 (0.008).

The likelihood ratio against the null model is approximately 26 with associated p-value near zero indicating that the inclusion of a temporal trend is statistically significant. The parameter estimates from fitting the GPD with log() = + t, t as above, are given by:

0 1.66 (0.161)

1 -0.02 (0.008)

-0.24 (0.069)

The likelihood ratio between this model and the model without any temporal trend is about 4.23, which is (slightly) greater than the 1, 0.95 critical value of 3.84, and the associated p-value of about 0.040; indicating that the inclusion of a temporal trend in the scale parameter is statistically significant, if only slightly. Of course, the apparent trend in Figure 6.3 is relatively small, but nevertheless apparent so that it is reasonable to include such a trend in the Poisson-GP model.



9Note: because of the naming convention used by extRemes the trigonometric transformations with periods of 365 days cannot exist simultaneously with periods of, for example, 365.25 days. By default, and in order to prevent accidental deletion of data, extRemes will not allow a transformation if there is already a data column with the same name. In the present example, if a period of 365 is desired, the new names would also be obs.sin365 and obs.cos365; so both of these columns must be removed (e.g, using the Scrubber function under File) before invoking this transformation.

Back to Top

Back to Table of Contents