Estimator of small exceedance probabilities and large return periods using censored GPD-MLE
Computes estimates of a small exceedance probability P(X>q) or large return period 1/P(X>q) using the GPD-ML estimator adapted for right censoring.
cProbGPD(data, censored, gamma1, sigma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnGPD(data, censored, gamma1, sigma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of n observations. |
censored |
A logical vector of length n indicating if an observation is censored. |
gamma1 |
Vector of n-1 estimates for the EVI obtained from |
sigma1 |
Vector of n-1 estimates for σ_1 obtained from |
q |
The used large quantile (we estimate P(X>q) or 1/P(X>q) for q large). |
plot |
Logical indicating if the estimates should be plotted as a function of k, default is |
add |
Logical indicating if the estimates should be added to an existing plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The probability is estimated as
\hat{P}(X>q)=(1-km) \times (1+ \hat{γ}_1/a_{k,n} \times (q-Z_{n-k,n}))^{-1/\hat{γ}_1}
with Z_{i,n} the i-th order statistic of the data, \hat{γ}_1 the generalised Hill estimator adapted for right censoring and km the Kaplan-Meier estimator for the CDF evaluated in Z_{n-k,n}. The value a is defined as
a_{k,n} = \hat{σ}_1 / \hat{p}_k
with \hat{σ}_1 the ML estimate for σ_1 and \hat{p}_k the proportion of the k largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter k. |
P |
Vector of the corresponding probability estimates, only returned for |
R |
Vector of the corresponding estimates for the return period, only returned for |
q |
The used large quantile. |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
# Set seed set.seed(29072016) # Pareto random sample X <- rpareto(500, shape=2) # Censoring variable Y <- rpareto(500, shape=1) # Observed sample Z <- pmin(X, Y) # Censoring indicator censored <- (X>Y) # GPD-MLE estimator adapted for right censoring cpot <- cGPDmle(Z, censored=censored, plot=TRUE) # Exceedance probability q <- 10 cProbGPD(Z, gamma1=cpot$gamma1, sigma1=cpot$sigma1, censored=censored, q=q, plot=TRUE) # Return period cReturnGPD(Z, gamma1=cpot$gamma1, sigma1=cpot$sigma1, censored=censored, q=q, plot=TRUE)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.