Frequentist Likelihood-Based Inference for Time Series Extremes

Cheeseboro wind gust data

We use an example dataset to illustrate the use of the two main functions in lite, that is, flite and returnLevel. The 744 by 10 matrix cheeseboro, available from the exdex package, contains hourly maximum wind gusts (in miles per hour) recorded at the Cheeseboro weather station near Thousand Oaks, Southern California, USA during the month of January over the period 2000-2009. Column \(i\) contains the hourly values for year \(1999+i\). These data are sourced from the Cheeseboro page of the Remote Automated Weather Stations USA Climate Archive. The following plot shows the general behaviour of these data.

For estimating \((p_u, \sigma_u, \xi)\) we treat observations from different years as being in separate clusters. In the code below we do not need to set the argument cluster explicitly because when the data are supplied as a matrix the function flite sets cluster automatically to achieve this. A total of 42 observations are missing in cheeseboro. This causes no problems for inferences about \((p_u, \sigma_u, \xi)\). For making inferences about \(\theta\) using the \(K\)-gaps model the kgaps function splits the data further into sequences of non-missing values and constructs the sample \(K\)-gaps separately for each sequence.

Inferences for model parameters

We need to choose values for the tuning parameters threshold \(u\) and \(K\)-gaps run parameter \(K\). There is no definitive way to do this but the information matrix test developed in Süveges and Davison (2010) can help to make this choice. Based on the \(K\)-gaps section of the introductory vignette for the exdex package we choose \(u = 45\)mph and \(K = 3\).

library(lite)
cdata <- exdex::cheeseboro
# Each column of the matrix cdata corresponds to data from a different year
# flite() sets cluster automatically to correspond to column (year)
cfit <- flite(cdata, u = 45, k = 3)
summary(cfit)
#> 
#> Call:
#> flite(data = cdata, u = 45, k = 3)
#> 
#>          Estimate Std. Error
#> p[u]      0.02771   0.005988
#> sigma[u]  9.27400   2.071000
#> xi       -0.09368   0.084250
#> theta     0.24050   0.023360

summary(cfit, adjust = FALSE)
#> 
#> Call:
#> flite(data = cdata, u = 45, k = 3)
#> 
#>          Estimate Std. Error
#> p[u]      0.02771   0.001908
#> sigma[u]  9.27400   0.858100
#> xi       -0.09368   0.061040
#> theta     0.24050   0.023360

The effect of the adjustment of the log-likelihood for clustering is to increase the estimated standard errors for the parameters \(p_u\) and \((\sigma_u, \xi)\) in comparison to those obtained with no adjustment.

A plot method for objects inheriting from class "flite" enables us to inspect the log-likelihood functions for \(p_u\), \((\sigma_u, \xi)\) and \(\theta\). This first set of plots use the default vertical adjustment to the log-likelihoods for \(p_u\) and \((\sigma_u, \xi)\).

plot(cfit)

The line superimposed on the on the plot on the top right shows the boundary of the log-likelihood imposed by the constraint \(\xi > -\sigma_u / x_{(n)}\), where \(x_{(n)}\) is the largest threshold excess. Using the vertical adjustment ensures that this constraint is respected after adjustment. However, if we use one of the horizontal adjustments, cholesky or spectral, then the adjusted log-likelihood has non-zero density values outside of the correct support. This is why the vertical adjustment is the default.

plot(cfit, which = "gp", adj_type = "spectral")

Inferences for return levels

The \(m\)-year return level is the value that is exceeded with probability \(1/m\) in any given year. To infer this from time series data we need to take account of the (intended) sampling frequency of the data. In the context of the Cheeseboro data, a year is effectively one month because we are making inferences about extremal behaviour in January. If there are no missing values then there are \(31 \times 24 = 744\) hourly values per January. We denote this by \(n_y\) (number of observations per year).

The \(m\)-year return level is given by \[x_m = u+\frac{\sigma_u}{\xi} \left[ p_u^\xi \left\{ 1-\left( 1-\displaystyle \frac1m\right)^{1/n_{y}\theta} \right\}^{-\xi} - 1 \right].\]

If \(m n_y \theta\) is large then \[x_m \approx u+\frac{\sigma_u}{\xi} \left[ (m \, n_{y}\, \theta \, p_u)^\xi - 1 \right].\] The returnLevel function performs inferences about return levels. It uses the former expression for \(x_m\). By default it calculates two types of (95%) confidence intervals: (a) intervals based on approximate large-sample normality of the maximum likelihood estimator for return level, which are symmetric about the point estimate, and (b) profile likelihood-based intervals based on an (adjusted) log-likelihood.

rl <- returnLevel(cfit, m = 100, level = 0.95, ny = 31 * 24)
summary(rl)
#> 
#> Call:
#> returnLevel(x = cfit, m = 100, level = 0.95, ny = 31 * 24)
#> 
#>         Estimate Std. Error
#> m = 100 88.61986   9.600889

rl
#> 
#> Call:
#> returnLevel(x = cfit, m = 100, level = 0.95, ny = 31 * 24)
#> 
#> MLE and 95% confidence limits for the 100-year return level
#> 
#> Normal interval:
#>  lower     mle   upper  
#>  69.80   88.62  107.44  
#> Profile likelihood-based interval:
#>  lower     mle   upper  
#>  75.89   88.62  125.43

plot(rl)

References