The POT Package: An R package to model peaks over threshold

The POT package aims to provide operational tools to analyze POT. This package relies on the EVT to model the tail of any continuous distribution. Tail modelling, in particular POT modelling, is of great importance for many financial and environmental applications.

The POT package was first committed to the CRAN in April 2005 and is still in active development.The main motivation was to provide practical tools for probabilistic modelling of high flood flows. However, the strength of the EVT is that results do not depend on the process to be modelled. Thus, one can use the POT package to analyze precipitations, floods, financial times series, earthquakes and so on...

Features

The POT package can perform univariate and bivariate extreme value analysis; first order Markov chains can also be considered. For instance, the (univariate) GPD is currently fitted using 18 estimators. These estimators rely on three different techniques:

Likelihood maximization: MLE, LME, MPLE
Moment Approaches: MOM, PWM, MED
Distance Minimization: MDPD and MGF estimators.

Contrary to the univariate case, there is no finite parametrisation to model bivariate exceedances over thresholds. The POT packages allows 6 parametrisation for the bivariate GPD: the logistic, negative logistic and mixed models - with their respective asymmetric counterparts.

Lastly, first order Markov chains can be fitted using the bivariate GPD for the joint distribution of two consecutive observations.

Screen Shots

The POT Package in a Few Lines

In this section, we explicit some of the most useful function of the package. However, for a full description, users may want to have a look to the package vignette and the html help of the package.

GPD Computations:

##Simulate a sample from a GPD(0,1,0.2):
x <- rgpd(100, 0, 1, 0.2) [1] 0.123422302 0.297063966 ...
##Evaluate density at x=3 and probability of non-exceedance:
dgpd(3, 0, 1, 0.2); pgpd(3, 0, 1, 0.2) [1] 0.05960464
[1] 0.9046326
##Compute the quantile with non-exceedance probability 0.95:
qgpd(0.95, 0, 1, 0.2) [1] 4.102821
##What about the bivariate case? Just the same

y <- rbvgpd(100, model = "alog", alpha = 0.2,
	     asCoef1 = 0.8, asCoef2 = 0.2, mar1 = c(0, 1,
	     0.2), mar2 = c(10, 1, 0.5))

[,1] [,2]
[1,] 1.3671957344 10.32043
[2,] 0.3074876998 10.16365
[3,] 0.2277643365 10.64219
[4,] 1.3870434019 10.57433
[5,] 1.3678122283 11.74367
[6,] 1.7164821847 10.56805
[7,] 0.3058534945 11.29733
...
##Evaluate the probability to not exceed (5,14)

pbvgpd(c(3,15), model = "alog", alpha = 0.2,
    asCoef1 = 0.8, asCoef2 = 0.2, mar1 = c(0, 1,
    0.2), mar2 = c(10, 1, 0.5))

[1] 0.8450499

GPD Fitting

##Maximum likelihood estimate (threshold = 0):
mle <- fitgpd(x, 0) Estimator: MLE
Deviance: 115.4406
AIC: 119.4406

Varying Threshold: FALSE

Threshold Call: 0
Number Above: 50
Proportion Above: 1

Estimates
scale shape
1.24924 -0.06813
...
##Probability Weighted Moments:
pwu <- fitgpd(x, 0, "pwmu") Estimator: PWMU

Varying Threshold: FALSE

Threshold Call: 0
Number Above: 100
Proportion Above: 1

Estimates
scale shape
0.9714 0.1837
...
##Maximum Goodness-of-Fit estimators:
adr <- fitgpd(x, 0, "mgf", stat = "ADR") Estimator: MGF
Statistic: ADR

Varying Threshold: FALSE

Threshold Call: 0
Number Above: 100
Proportion Above: 1

Estimates
scale shape
0.9438 0.1938

Optimization Information
Convergence: successful
Function Evaluations: 21
Gradient Evaluations: 10
##Specifying a known parameter:
fitgpd(x, 0, "mple", shape = 0.2) Estimator: MPLE

Penalized Deviance: 231.0240
Penalized AIC: 233.0240

Varying Threshold: FALSE

Threshold Call: 0
Number Above: 100
Proportion Above: 1

Estimates
scale
0.955
...
##Specifying starting values for numerical optimizations:
fitgpd(x, 0, "mdpd", start = list(scale = 1, shape = 0.2)) Estimator: MDPD

Varying Threshold: FALSE

Threshold Call: 0
Number Above: 100
Proportion Above: 1

Estimates
scale shape
0.9693 0.1809

Optimization Information
Convergence: successful
Function Evaluations: 8
Gradient Evaluations: 5
##Fit a bivariate GPD with a logistic dependence:
log <- fitbvgpd(y, c(0,10), "log") Call: fitbvgpd(data = y, threshold = c(0, 10), model = "log")
Estimator: MLE
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = 0.063
Deviance: 489.5613
AIC: 499.5613

Marginal Threshold: 0 10
Marginal Number Above: 100 100
Marginal Proportion Above: 1 1
...

Plots

##Generic function for the univariate and bivariate cases:
plot(mle);

plot(log)

##Return level plots:
retlev(mle, npy = 2);

retlev(log)

##Probability-Probability and Q-Q plots:
pp(mle); qq(mle)

##Plot the density:
dens(mle)

##Plot the Pickands' dependence function:
pickdep(log)

##Spectral density plot:
specdens(log)

##Profile Likelihood (quantiles):
confint(mle, prob = 0.95)

##Profile Likelihood (parameters):
confint(mle, "scale"); confint(mle, "shape")

Manuals

We have written a package vignette to help new users. This user's guide is a part of the package - just run vignette(POT) once the package is loaded.
For a quick overview, one can have a look at the R News (7):1 (2OO7) article or by running the univariate and bivariate demos - i.e. demo(POT) and demo(bvPOT).

Contribute to the Project

If you are interested in joining the project, you must first create an R-forge account. Then, just go here.

Contact

Any suggestions, feature requests, bugs: select the appropriate tracker
Author: Mathieu Ribatet; Maintainer: Christophe Dutang.