Calculates maximum likelihood estimates of the extremal index \(\theta\)based on the \(K\)-gaps model for threshold inter-exceedances times ofSuveges and Davison (2010).
kgaps(data,u, k=1, inc_cens=TRUE)A numeric vector or numeric matrix of raw data. Ifdata is a matrix then the log-likelihood is constructed as the sum of (independent) contributions from different columns. A common situation is where each column relates to a different year.
Ifdata contains missing values thensplit_by_NAs is used to divide the data further into sequences of non-missing values, stored in different columns in a matrix. Again, the log-likelihood is constructed as a sum of contributions from different columns.
A numeric scalar. Extreme value threshold applied to data.
A non-negative numeric scalar. Run parameter \(K\), as defined inSuveges and Davison (2010). Threshold inter-exceedances times that arenot larger thank units are assigned to the same cluster, resultingin a \(K\)-gap equal to zero. Specifically, the \(K\)-gap \(S\)corresponding to an inter-exceedance time of \(T\) is given by\(S = \max(T - K, 0)\). In practice, \(k\) shouldbe no smaller than 1, because when \(k\) is less than 1 the estimateof \(\theta\) is always equal to 1.
A logical scalar indicating whether or not to includecontributions from right-censored inter-exceedance times, relating to thefirst and last observations. It is known that these times are greaterthan or equal to the time observed. See Attalides (2015) for details.Ifdata has multiple columns then there will be right-censoredfirst and last inter-exceedance times for each column.
An object (a list) of classc("kgaps", "exdex") containing
thetaThe maximum likelihood estimate (MLE) of \(\theta\).
seThe estimated standard error of the MLE, calculated using an algebraic expression for the observed information. Ifk = 0 thense is returned as0.
se_expThe estimated standard error of the MLE, calculated using an algebraic expression for the expected information. If the estimate of \(\theta\) is 0 or 1 thense_exp isNA.
ssThe list of summary statistics returned fromkgaps_stat.
k, u, inc_censThe input values ofk,u andinc_cens.
max_loglikThe value of the log-likelihood at the MLE.
callThe call tokgaps.
Ifinc_cens = FALSE then the maximum likelihood estimate of the extremal index \(\theta\) under the \(K\)-gaps model of Suveges and Davison (2010) is calculated.
Ifinc_cens = TRUE then information from right-censored first and last inter-exceedance times is also included in the likelihood to be maximized, following Attalides (2015). The form of the log-likelihood is given in theDetails section ofkgaps_stat.
It is possible that the estimate of \(\theta\) is equal to 1, and also possible that it is equal to 0.kgaps_stat explains the respective properties of the data that cause these events to occur.
Suveges, M. and Davison, A. C. (2010) Model misspecification in peaks over threshold analysis,Annals of Applied Statistics,4(1), 203-221.doi:10.1214/09-AOAS292
Attalides, N. (2015) Threshold-based extreme value modelling, PhD thesis, University College London.https://discovery.ucl.ac.uk/1471121/1/Nicolas_Attalides_Thesis.pdf
kgaps_confint to estimate confidence intervals for \(\theta\).
kgaps_methods for S3 methods for"kgaps" objects.
kgaps_imt for the information matrix test, which may be used to inform the choice of the pair (u, k).
choose_uk for a diagnostic plot based onkgaps_imt.
kgaps_stat for the calculation of sufficient statistics for the \(K\)-gaps model.
kgaps_post in therevdbayes package for Bayesian inference about \(\theta\) using the \(K\)-gaps model.
### S&P 500 indexu<-quantile(sp500, probs=0.60)theta<-kgaps(sp500,u)theta#>#> Call:#> kgaps(data = sp500, u = u)#>#> Estimate of the extremal index theta:#> theta#> 0.6953summary(theta)#>#> Call:#> kgaps(data = sp500, u = u)#>#> Estimate Std. Error#> theta 0.6953 0.007234coef(theta)#> theta#> 0.6953391nobs(theta)#> [1] 2901vcov(theta)#> theta#> theta 5.232578e-05logLik(theta)#> 'log Lik.' -3811.894 (df=1)### Newlyn sea surgesu<-quantile(newlyn, probs=0.60)theta<-kgaps(newlyn,u, k=2)theta#>#> Call:#> kgaps(data = newlyn, u = u, k = 2)#>#> Estimate of the extremal index theta:#> theta#> 0.1758summary(theta)#>#> Call:#> kgaps(data = newlyn, u = u, k = 2)#>#> Estimate Std. Error#> theta 0.1758 0.009211### Cheeseboro wind guststheta<-kgaps(cheeseboro,45, k=3)theta#>#> Call:#> kgaps(data = cheeseboro, u = 45, k = 3)#>#> Estimate of the extremal index theta:#> theta#> 0.2405summary(theta)#>#> Call:#> kgaps(data = cheeseboro, u = 45, k = 3)#>#> Estimate Std. Error#> theta 0.2405 0.02336