| Title: | Functions from "Reinsurance: Actuarial and Statistical Aspects" |
|---|---|
| Description: | Functions from the book "Reinsurance: Actuarial and Statistical Aspects" (2017) by Hansjoerg Albrecher, Jan Beirlant and Jef Teugels <https://www.wiley.com/en-us/Reinsurance%3A+Actuarial+and+Statistical+Aspects-p-9780470772683>. |
| Authors: | Tom Reynkens [aut, cre]  ,
Roel Verbelen [aut] (R code for Mixed Erlang distribution,
<https://orcid.org/0000-0002-2347-9240>),
Anastasios Bardoutsos [ctb] (Original R code for cEPD estimator),
Dries Cornilly [ctb] (Original R code for EVT estimators for truncated
data),
Yuri Goegebeur [ctb] (Original S-Plus code for basic EVT estimators),
Klaus Herrmann [ctb] (Original R code for GPD estimator) |
| Maintainer: | Tom Reynkens <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.0.15 |
| Built: | 2025-01-29 06:25:43 UTC |
| Source: | https://github.com/treynkens/reins |
Density, distribution function, quantile function and random generation for the Burr distribution (type XII).
dburr(x, alpha, rho, eta = 1, log = FALSE) pburr(x, alpha, rho, eta = 1, lower.tail = TRUE, log.p = FALSE) qburr(p, alpha, rho, eta = 1, lower.tail = TRUE, log.p = FALSE) rburr(n, alpha, rho, eta = 1)dburr(x, alpha, rho, eta = 1, log = FALSE) pburr(x, alpha, rho, eta = 1, lower.tail = TRUE, log.p = FALSE) qburr(p, alpha, rho, eta = 1, lower.tail = TRUE, log.p = FALSE) rburr(n, alpha, rho, eta = 1)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
alpha |
The |
rho |
The |
eta |
The |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the Burr distribution is equal to
for all and otherwise. We need that , and .
Beirlant et al. (2004) uses parameters which correspond to , and .
dburr gives the density function evaluated in , pburr the CDF evaluated in and qburr the quantile function evaluated in . The length of the result is equal to the length of or .
rburr returns a random sample of length .
Tom Reynkens.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dburr(x, alpha=2, rho=-1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pburr(x, alpha=2, rho=-1), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dburr(x, alpha=2, rho=-1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pburr(x, alpha=2, rho=-1), xlab="x", ylab="CDF", type="l")
Computes the EPD estimates adapted for right censored data.
cEPD(data, censored, rho = -1, beta = NULL, logk = FALSE, plot = FALSE, add = FALSE, main = "EPD estimates of the EVI", ...)cEPD(data, censored, rho = -1, beta = NULL, logk = FALSE, plot = FALSE, add = FALSE, main = "EPD estimates of the EVI", ...)
data |
Vector of |
censored |
A logical vector of length |
rho |
A parameter for the |
beta |
Parameter for EPD ( |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The function EPD uses which is equal to .
This estimator is only suitable for right censored data.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma1 |
Vector of the corresponding estimates for the |
kappa1 |
Vector of the corresponding MLE estimates for the |
beta |
Vector of estimates for (or values of) the |
Delta |
Difference between |
Tom Reynkens based on R code from Anastasios Bardoutsos.
Beirlant, J., Bardoutsos, A., de Wet, T. and Gijbels, I. (2016). "Bias Reduced Tail Estimation for Censored Pareto Type Distributions." Statistics & Probability Letters, 109, 78–88.
Fraga Alves, M.I. , Gomes, M.I. and de Haan, L. (2003). "A New Class of Semi-parametric Estimators of the Second Order Parameter." Portugaliae Mathematica, 60, 193–214.
# 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) # EPD estimator adapted for right censoring cepd <- cEPD(Z, censored=censored, plot=TRUE)# 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) # EPD estimator adapted for right censoring cepd <- cEPD(Z, censored=censored, plot=TRUE)
Exponential QQ-plot adapted for right censored data.
cExpQQ(data, censored, plot = TRUE, main = "Exponential QQ-plot", ...)cExpQQ(data, censored, plot = TRUE, main = "Exponential QQ-plot", ...)
data |
Vector of |
censored |
A logical vector of length |
plot |
Logical indicating if the quantiles should be plotted in an exponential QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The exponential QQ-plot adapted for right censoring is given by
for
with the -th order statistic of the data and the Kaplan-Meier estimator for the CDF.
Hence, it has the same empirical quantiles as an ordinary exponential QQ-plot but replaces the theoretical quantiles by .
This QQ-plot is only suitable for right censored data.
In Beirlant et al. (2007), only a Pareto QQ-plot adapted for right-censored data is proposed. This QQ-plot is constructed using the same ideas, but is not described in the paper.
A list with following components:
eqq.the |
Vector of the theoretical quantiles, see Details. |
eqq.emp |
Vector of the empirical quantiles from the data. |
Tom Reynkens
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
ExpQQ, cLognormalQQ, cParetoQQ, cWeibullQQ, KaplanMeier
# 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) # Exponential QQ-plot adapted for right censoring cExpQQ(Z, censored=censored)# 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) # Exponential QQ-plot adapted for right censoring cExpQQ(Z, censored=censored)
Computes the generalised Hill estimates adapted for right censored data.
cgenHill(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Generalised Hill estimates of the EVI", ...)cgenHill(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Generalised Hill estimates of the EVI", ...)
data |
Vector of |
censored |
A logical vector of length |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The generalised Hill estimator adapted for right censored data is equal to the ordinary generalised Hill estimator divided by the proportion of the largest observations that is non-censored.
This estimator is only suitable for right censored data.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma1 |
Vector of the corresponding generalised Hill estimates. |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
genHill, cHill, cProbGH, cQuantGH
# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE)# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE)
Computes ML estimates of fitting GPD to peaks over a threshold adapted for right censoring.
cGPDmle(data, censored, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...) cPOT(data, censored, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...)cGPDmle(data, censored, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...) cPOT(data, censored, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...)
data |
Vector of |
censored |
A logical vector of length |
start |
Vector of length 2 containing the starting values for the optimisation. The first element
is the starting value for the estimator of |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The GPD-MLE estimator for the EVI adapted for right censored data is equal to the ordinary GPD-MLE estimator for the EVI divided by the proportion of the largest observations that is non-censored. The estimates for are the ordinary GPD-MLE estimates for .
This estimator is only suitable for right censored data.
cPOT is the same function but with a different name for compatibility with POT.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma1 |
Vector of the corresponding MLE estimates for the |
sigma1 |
Vector of the corresponding MLE estimates for the |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
GPDmle, cProbGPD, cQuantGPD, cEPD
# 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-ML estimator adapted for right censoring cpot <- cGPDmle(Z, censored=censored, plot=TRUE)# 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-ML estimator adapted for right censoring cpot <- cGPDmle(Z, censored=censored, plot=TRUE)
Computes the Hill estimator for positive extreme value indices, adapted for right censoring, as a function of the tail parameter (Beirlant et al., 2007).
Optionally, these estimates are plotted as a function of .
cHill(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Hill estimates of the EVI", ...)cHill(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Hill estimates of the EVI", ...)
data |
Vector of |
censored |
A logical vector of length |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Hill estimator adapted for right censored data is equal to the ordinary Hill estimator divided by the proportion of the largest observations that is non-censored.
This estimator is only suitable for right censored data, use icHill for interval censored data.
See Section 4.3.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma1 |
Vector of the corresponding Hill estimates. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
Hill, icHill, cParetoQQ, cProb, cQuant
# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE)# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE)
Log-normal QQ-plot adapted for right censored data.
cLognormalQQ(data, censored, plot = TRUE, main = "Log-normal QQ-plot", ...)cLognormalQQ(data, censored, plot = TRUE, main = "Log-normal QQ-plot", ...)
data |
Vector of |
censored |
A logical vector of length |
plot |
Logical indicating if the quantiles should be plotted in a log-normal QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The log-normal QQ-plot adapted for right censoring is given by
for
with the -th order statistic of the data, the quantile function of the standard normal distribution and the Kaplan-Meier estimator for the CDF.
Hence, it has the same empirical quantiles as an ordinary log-normal QQ-plot but replaces the theoretical quantiles by .
This QQ-plot is only suitable for right censored data.
In Beirlant et al. (2007), only a Pareto QQ-plot adapted for right-censored data is proposed. This QQ-plot is constructed using the same ideas, but is not described in the paper.
A list with following components:
lqq.the |
Vector of the theoretical quantiles, see Details. |
lqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
LognormalQQ, cExpQQ, cParetoQQ, cWeibullQQ, KaplanMeier
# 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) # Log-normal QQ-plot adapted for right censoring cLognormalQQ(Z, censored=censored)# 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) # Log-normal QQ-plot adapted for right censoring cLognormalQQ(Z, censored=censored)
Computes the Method of Moment estimates adapted for right censored data.
cMoment(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Moment estimates of the EVI", ...)cMoment(data, censored, logk = FALSE, plot = FALSE, add = FALSE, main = "Moment estimates of the EVI", ...)
data |
Vector of |
censored |
A logical vector of length |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The moment estimator adapted for right censored data is equal to the ordinary moment estimator divided by the proportion of the largest observations that is non-censored.
This estimator is only suitable for right censored data.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma1 |
Vector of the corresponding moment estimates. |
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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE)# 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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE)
Pareto QQ-plot adapted for right censored data.
cParetoQQ(data, censored, plot = TRUE, main = "Pareto QQ-plot", ...)cParetoQQ(data, censored, plot = TRUE, main = "Pareto QQ-plot", ...)
data |
Vector of |
censored |
A logical vector of length |
plot |
Logical indicating if the quantiles should be plotted in a Pareto QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Pareto QQ-plot adapted for right censoring is given by
for
with the -th order statistic of the data and the Kaplan-Meier estimator for the CDF.
Hence, it has the same empirical quantiles as an ordinary Pareto QQ-plot but replaces the theoretical quantiles by .
This QQ-plot is only suitable for right censored data, use icParetoQQ for interval censored data.
A list with following components:
pqq.the |
Vector of the theoretical quantiles, see Details. |
pqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
ParetoQQ, icParetoQQ, cExpQQ, cLognormalQQ, cWeibullQQ, cHill, KaplanMeier
# 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) # Pareto QQ-plot adapted for right censoring cParetoQQ(Z, censored=censored)# 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) # Pareto QQ-plot adapted for right censoring cParetoQQ(Z, censored=censored)
Computes estimates of a small exceedance probability or large return period using the estimates for the EVI obtained from the Hill estimator adapted for right censoring.
cProb(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturn(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)cProb(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturn(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with the -th order statistic of the data,
the Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in .
A list with following components:
k |
Vector of the values of the tail parameter |
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
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
cHill, cQuant, Prob, KaplanMeier
# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProb(Z, censored=censored, gamma1=chill$gamma1, q=q, plot=TRUE) # Return period cReturn(Z, censored=censored, gamma1=chill$gamma1, q=q, plot=TRUE)# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProb(Z, censored=censored, gamma1=chill$gamma1, q=q, plot=TRUE) # Return period cReturn(Z, censored=censored, gamma1=chill$gamma1, q=q, plot=TRUE)
Computes estimates of a small exceedance probability or large return period using the parameters from the EPD fit adapted for right censoring.
cProbEPD(data, censored, gamma1, kappa1, beta, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnEPD(data, censored, gamma1, kappa1, beta, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)cProbEPD(data, censored, gamma1, kappa1, beta, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnEPD(data, censored, gamma1, kappa1, beta, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
kappa1 |
Vector of |
beta |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with
the CDF of the EPD with estimated parameters , and
and the Kaplan-Meier estimator for the CDF evaluated in (the -th largest data point).
A list with following components:
k |
Vector of the values of the tail parameter |
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.
Beirlant, J., Bardoutsos, A., de Wet, T. and Gijbels, I. (2016). "Bias Reduced Tail Estimation for Censored Pareto Type Distributions." Statistics & Probability Letters, 109, 78–88.
cEPD, ProbEPD, Prob, KaplanMeier
# 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) # EPD estimator adapted for right censoring cepd <- cEPD(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbEPD(Z, censored=censored, gamma1=cepd$gamma1, kappa1=cepd$kappa1, beta=cepd$beta, q=q, plot=TRUE) # Return period cReturnEPD(Z, censored=censored, gamma1=cepd$gamma1, kappa1=cepd$kappa1, beta=cepd$beta, q=q, plot=TRUE)# 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) # EPD estimator adapted for right censoring cepd <- cEPD(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbEPD(Z, censored=censored, gamma1=cepd$gamma1, kappa1=cepd$kappa1, beta=cepd$beta, q=q, plot=TRUE) # Return period cReturnEPD(Z, censored=censored, gamma1=cepd$gamma1, kappa1=cepd$kappa1, beta=cepd$beta, q=q, plot=TRUE)
Computes estimates of a small exceedance probability or large return period using the estimates for the EVI obtained from the generalised Hill estimator adapted for right censoring.
cProbGH(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnGH(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)cProbGH(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnGH(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with the -th order statistic of the data, the generalised Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ordinary Hill estimator
and the proportion of the largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
cQuantGH, cgenHill, ProbGH, cProbMOM, KaplanMeier
# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbGH(Z, censored=censored, gamma1=cghill$gamma1, q=q, plot=TRUE) # Return period cReturnGH(Z, censored=censored, gamma1=cghill$gamma1, q=q, plot=TRUE)# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbGH(Z, censored=censored, gamma1=cghill$gamma1, q=q, plot=TRUE) # Return period cReturnGH(Z, censored=censored, gamma1=cghill$gamma1, q=q, plot=TRUE)
Computes estimates of a small exceedance probability or large return period 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", ...)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 |
censored |
A logical vector of length |
gamma1 |
Vector of |
sigma1 |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with the -th order statistic of the data, the generalised Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ML estimate for
and the proportion of the largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
cQuantGPD, cGPDmle, ProbGPD, Prob, KaplanMeier
# 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)# 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)
Computes estimates of a small exceedance probability or large return period using the Method of Moments estimates for the EVI adapted for right censoring.
cProbMOM(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnMOM(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)cProbMOM(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) cReturnMOM(data, censored, gamma1, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with the -th order statistic of the data, the MOM estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ordinary Hill estimator
and the proportion of the largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
cQuantMOM, cMoment, ProbMOM, Prob, KaplanMeier
# 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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbMOM(Z, censored=censored, gamma1=cmom$gamma1, q=q, plot=TRUE) # Return period cReturnMOM(Z, censored=censored, gamma1=cmom$gamma1, q=q, plot=TRUE)# 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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE) # Small exceedance probability q <- 10 cProbMOM(Z, censored=censored, gamma1=cmom$gamma1, q=q, plot=TRUE) # Return period cReturnMOM(Z, censored=censored, gamma1=cmom$gamma1, q=q, plot=TRUE)
Computes estimates of large quantiles using the estimates for the EVI obtained from the Hill estimator adapted for right censoring.
cQuant(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)cQuant(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 quantile is estimated as
with the -th order statistic of the data,
the Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in .
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
cHill, cProb, Quant, KaplanMeier
# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuant(Z, gamma1=chill$gamma, censored=censored, p=p, plot=TRUE)# 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) # Hill estimator adapted for right censoring chill <- cHill(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuant(Z, gamma1=chill$gamma, censored=censored, p=p, plot=TRUE)
Computes estimates of large quantiles using the estimates for the EVI obtained from the generalised Hill estimator adapted for right censoring.
cQuantGH(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)cQuantGH(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 quantile is estimated as
with the -th order statistic of the data, the generalised Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ordinary Hill estimator
and the proportion of the largest observations that is non-censored, and
with
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
cProbGH, cgenHill, QuantGH, Quant, KaplanMeier
# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuantGH(Z, gamma1=cghill$gamma, censored=censored, p=p, plot=TRUE)# 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) # Generalised Hill estimator adapted for right censoring cghill <- cgenHill(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuantGH(Z, gamma1=cghill$gamma, censored=censored, p=p, plot=TRUE)
Computes estimates of large quantiles using the estimates for the EVI obtained from the GPD-ML estimator adapted for right censoring.
cQuantGPD(data, censored, gamma1, sigma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)cQuantGPD(data, censored, gamma1, sigma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
sigma1 |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 quantile is estimated as
ith the -th order statistic of the data, the generalised Hill estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ML estimate for
and the proportion of the largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
cProbGPD, cGPDmle, QuantGPD, Quant, KaplanMeier
# 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) # Large quantile p <- 10^(-4) cQuantGPD(Z, gamma1=cpot$gamma1, sigma1=cpot$sigma1, censored=censored, p=p, plot=TRUE)# 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) # Large quantile p <- 10^(-4) cQuantGPD(Z, gamma1=cpot$gamma1, sigma1=cpot$sigma1, censored=censored, p=p, plot=TRUE)
Computes estimates of large quantiles using the estimates for the EVI obtained from the MOM estimator adapted for right censoring.
cQuantMOM(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)cQuantMOM(data, censored, gamma1, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
censored |
A logical vector of length |
gamma1 |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 quantile is estimated as
ith the -th order statistic of the data, the MOM estimator adapted for right censoring and the Kaplan-Meier estimator for the CDF evaluated in . The value is defined as
with the ordinary Hill estimator
and the proportion of the largest observations that is non-censored.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens
Einmahl, J.H.J., Fils-Villetard, A. and Guillou, A. (2008). "Statistics of Extremes Under Random Censoring." Bernoulli, 14, 207–227.
cProbMOM, cMoment, QuantMOM, Quant, KaplanMeier
# 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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuantMOM(Z, censored=censored, gamma1=cmom$gamma1, p=p, plot=TRUE)# 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) # Moment estimator adapted for right censoring cmom <- cMoment(Z, censored=censored, plot=TRUE) # Large quantile p <- 10^(-4) cQuantMOM(Z, censored=censored, gamma1=cmom$gamma1, p=p, plot=TRUE)
Hill-type estimator for the conditional Extreme Value Index (EVI) adapted for censored data.
crHill(x, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"), logk = FALSE, plot = FALSE, add = FALSE, main = "", ...)crHill(x, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"), logk = FALSE, plot = FALSE, add = FALSE, main = "", ...)
x |
Value of the conditioning variable |
Xtilde |
Vector of length |
Ytilde |
Vector of length |
censored |
A logical vector of length |
h |
Bandwidth of the non-parametric estimator. |
kernel |
Kernel of the non-parametric estimator. One of |
logk |
Logical indicating if the Hill-type estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
This is a Hill-type estimator of the EVI of given . The estimator uses the censored sample , for , where and are censored at the same time. We assume that and the censoring variable are conditionally independent given .
See Section 4.4.3 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding Hill-type estimates. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Conditional Pareto QQ-plot crParetoQQ(x=1, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=2) # Plot Hill-type estimates crHill(x=1, Xtilde, Ytilde, censored, h=2, plot=TRUE)# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Conditional Pareto QQ-plot crParetoQQ(x=1, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=2) # Plot Hill-type estimates crHill(x=1, Xtilde, Ytilde, censored, h=2, plot=TRUE)
Conditional Pareto QQ-plot adapted for right censored data.
crParetoQQ(x, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"), plot = TRUE, add = FALSE, main = "Pareto QQ-plot", type = "p", ...)crParetoQQ(x, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"), plot = TRUE, add = FALSE, main = "Pareto QQ-plot", type = "p", ...)
x |
Value of the conditioning variable |
Xtilde |
Vector of length |
Ytilde |
Vector of length |
censored |
A logical vector of length |
h |
Bandwidth of the non-parametric estimator for the conditional survival function ( |
kernel |
Kernel of the non-parametric estimator for the conditional survival function ( |
plot |
Logical indicating if the quantiles should be plotted in a Pareto QQ-plot, default is |
add |
Logical indicating if the quantiles should be added to an existing plot, default is |
main |
Title for the plot, default is |
type |
Type of the plot, default is |
... |
Additional arguments for the |
We construct a Pareto QQ-plot for conditional on using the censored sample , for , where and are censored at the same time. We assume that and the censoring variable are conditionally independent given .
The conditional Pareto QQ-plot adapted for right censoring is given by
for
with the -th order statistic of the censored data and the non-parametric estimator for the conditional CDF of Akritas and Van Keilegom (2003), see crSurv.
See Section 4.4.3 in Albrecher et al. (2017) for more details.
A list with following components:
pqq.the |
Vector of the theoretical quantiles, see Details. |
pqq.emp |
Vector of the empirical quantiles from the log-transformed |
Tom Reynkens
Akritas, M.G. and Van Keilegom, I. (2003). "Estimation of Bivariate and Marginal Distributions With Censored Data." Journal of the Royal Statistical Society: Series B, 65, 457–471.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Conditional Pareto QQ-plot crParetoQQ(x=1, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=2) # Plot Hill-type estimates crHill(x=1, Xtilde, Ytilde, censored, h=2, plot=TRUE)# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Conditional Pareto QQ-plot crParetoQQ(x=1, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=2) # Plot Hill-type estimates crHill(x=1, Xtilde, Ytilde, censored, h=2, plot=TRUE)
Non-parametric estimator of the conditional survival function of given for censored data, see Akritas and Van Keilegom (2003).
crSurv(x, y, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"))crSurv(x, y, Xtilde, Ytilde, censored, h, kernel = c("biweight", "normal", "uniform", "triangular", "epanechnikov"))
x |
The value of the conditioning variable |
y |
The value(s) of the variable |
Xtilde |
Vector of length |
Ytilde |
Vector of length |
censored |
A logical vector of length |
h |
Bandwidth of the non-parametric estimator. |
kernel |
Kernel of the non-parametric estimator. One of |
We estimate the conditional survival function
using the censored sample , for , where and are censored at the same time. We assume that and the censoring variable are conditionally independent given .
The estimator is given by
where when is censored and 0 otherwise. The weights are given by
when and 0 otherwise.
See Section 4.4.3 in Albrecher et al. (2017) for more details.
Estimates for as described above.
Tom Reynkens
Akritas, M.G. and Van Keilegom, I. (2003). "Estimation of Bivariate and Marginal Distributions With Censored Data." Journal of the Royal Statistical Society: Series B, 65, 457–471.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Plot estimates of the conditional survival function x <- 5 y <- seq(0, 5, 1/100) plot(y, crSurv(x, y, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=5), type="l", xlab="y", ylab="Conditional survival function")# Set seed set.seed(29072016) # Pareto random sample Y <- rpareto(200, shape=2) # Censoring variable C <- rpareto(200, shape=1) # Observed (censored) sample of variable Y Ytilde <- pmin(Y, C) # Censoring indicator censored <- (Y>C) # Conditioning variable X <- seq(1, 10, length.out=length(Y)) # Observed (censored) sample of conditioning variable Xtilde <- X Xtilde[censored] <- X[censored] - runif(sum(censored), 0, 1) # Plot estimates of the conditional survival function x <- 5 y <- seq(0, 5, 1/100) plot(y, crSurv(x, y, Xtilde=Xtilde, Ytilde=Ytilde, censored=censored, h=5), type="l", xlab="y", ylab="Conditional survival function")
Compute Conditional Tail Expectation (CTE) of the fitted spliced distribution.
CTE(p, splicefit) ES(p, splicefit)CTE(p, splicefit) ES(p, splicefit)
p |
The probability associated with the CTE (we estimate |
splicefit |
A |
The Conditional Tail Expectation is defined as
where is the premium of the excess-loss insurance with retention u.
If the CDF is continuous in , we have with
the Tail Value-at-Risk.
See Reynkens et al. (2017) and Section 4.6 of Albrecher et al. (2017) for more details.
The ES function is the same function as CTE but is deprecated.
Vector with the CTE corresponding to each element of .
Tom Reynkens with R code from Roel Verbelen for the mixed Erlang quantiles.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
qSplice, ExcessSplice, SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) p <- seq(0.01, 0.99, 0.01) # Plot of CTE plot(p, CTE(p, splicefit), type="l", xlab="p", ylab=bquote(CTE[1-p])) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) p <- seq(0.01, 0.99, 0.01) # Plot of CTE plot(p, CTE(p, splicefit), type="l", xlab="p", ylab=bquote(CTE[1-p])) ## End(Not run)
Weibull QQ-plot adapted for right censored data.
cWeibullQQ(data, censored, plot = TRUE, main = "Weibull QQ-plot", ...)cWeibullQQ(data, censored, plot = TRUE, main = "Weibull QQ-plot", ...)
data |
Vector of |
censored |
A logical vector of length |
plot |
Logical indicating if the quantiles should be plotted in a Weibull QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Weibull QQ-plot adapted for right censoring is given by
for
with the -th order statistic of the data and the Kaplan-Meier estimator for the CDF.
Hence, it has the same empirical quantiles as an ordinary Weibull QQ-plot but replaces the theoretical quantiles by .
This QQ-plot is only suitable for right censored data.
In Beirlant et al. (2007), only a Pareto QQ-plot adapted for right-censored data is proposed. This QQ-plot is constructed using the same ideas, but is not described in the paper.
A list with following components:
wqq.the |
Vector of the theoretical quantiles, see Details. |
wqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens
Beirlant, J., Guillou, A., Dierckx, G. and Fils-Villetard, A. (2007). "Estimation of the Extreme Value Index and Extreme Quantiles Under Random Censoring." Extremes, 10, 151–174.
WeibullQQ, cExpQQ, cLognormalQQ, cParetoQQ, KaplanMeier
# 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) # Weibull QQ-plot adapted for right censoring cWeibullQQ(Z, censored=censored)# 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) # Weibull QQ-plot adapted for right censoring cWeibullQQ(Z, censored=censored)
Fit the Extended Pareto Distribution (GPD) to the exceedances (peaks) over a threshold. Optionally, these estimates are plotted as a function of .
EPD(data, rho = -1, start = NULL, direct = FALSE, warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "EPD estimates of the EVI", ...)EPD(data, rho = -1, start = NULL, direct = FALSE, warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "EPD estimates of the EVI", ...)
data |
Vector of |
rho |
A parameter for the |
start |
Vector of length 2 containing the starting values for the optimisation. The first element
is the starting value for the estimator of |
direct |
Logical indicating if the parameters are obtained by directly maximising the log-likelihood function, see Details. Default is |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
We fit the Extended Pareto distribution to the relative excesses over a threshold (X/u).
The EPD has distribution function
with and .
The parameters are determined using MLE and there are two possible approaches:
maximise the log-likelihood directly (direct=TRUE) or follow the approach detailed in
Beirlant et al. (2009) (direct=FALSE). The latter approach uses the score functions of the log-likelihood.
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding estimates for the |
kappa |
Vector of the corresponding MLE estimates for the |
tau |
Vector of the corresponding estimates for the |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Joossens, E. and Segers, J. (2009). "Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions." Journal of Statistical Planning and Inference, 139, 2800–2815.
Fraga Alves, M.I. , Gomes, M.I. and de Haan, L. (2003). "A New Class of Semi-parametric Estimators of the Second Order Parameter." Portugaliae Mathematica, 60, 193–214.
data(secura) # EPD estimates for the EVI epd <- EPD(secura$size, plot=TRUE) # Compute return periods ReturnEPD(secura$size, 10^10, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE)data(secura) # EPD estimates for the EVI epd <- EPD(secura$size, plot=TRUE) # Compute return periods ReturnEPD(secura$size, 10^10, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE)
Fit the Extended Pareto Distribution (EPD) to data using Maximum Likelihood Estimation (MLE).
EPDfit(data, tau, start = c(0.1, 1), warnings = FALSE)EPDfit(data, tau, start = c(0.1, 1), warnings = FALSE)
data |
Vector of |
tau |
Value for the |
start |
Vector of length 2 containing the starting values for the optimisation. The first element
is the starting value for the estimator of |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A vector with the MLE estimate for the parameter of the EPD as the first component and the MLE estimate for the parameter of the EPD as the second component.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Joossens, E. and Segers, J. (2009). "Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions." Journal of Statistical Planning and Inference, 139, 2800–2815.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Fit EPD to last 500 observations res <- EPDfit(SOAdata/sort(soa$size)[500], tau=-1)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Fit EPD to last 500 observations res <- EPDfit(SOAdata/sort(soa$size)[500], tau=-1)
Create an S3 object using an EVT (Extreme Value Theory) fit.
EVTfit(gamma, endpoint = NULL, sigma = NULL)EVTfit(gamma, endpoint = NULL, sigma = NULL)
gamma |
Vector of estimates for |
endpoint |
Vector of endpoints (with the same length as |
sigma |
Vector of scale estimates for the GPD (with the same length as |
See Reynkens et al. (2017) and Section 4.3 of Albrecher et al. (2017) for details.
An S3 object containing the above input arguments.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)
Estimate premiums of excess-loss reinsurance with retention and limit using EPD estimates.
ExcessEPD(data, gamma, kappa, tau, R, L = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)ExcessEPD(data, gamma, kappa, tau, R, L = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)
data |
Vector of |
gamma |
Vector of |
kappa |
Vector of |
tau |
Vector of |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance, default is |
warnings |
Logical indicating if warnings are displayed, default is |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
We need that , the -th largest observation.
If this is not the case, we return NA for the premium. A warning will be issued in
that case if warnings=TRUE.
The premium for the excess-loss insurance with retention and limit is given by
where is the premium of the excess-loss insurance with retention . When , the premium is equal to .
We estimate by
with and the estimates for the parameters of the EPD.
See Section 4.6 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
premium |
The corresponding estimates for the premium. |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
data(secura) # EPD estimator epd <- EPD(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessEPD(secura$size, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, R=R, ylim=c(0,2*10^4))data(secura) # EPD estimator epd <- EPD(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessEPD(secura$size, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, R=R, ylim=c(0,2*10^4))
Estimate premiums of excess-loss reinsurance with retention and limit using GPD-MLE estimates.
ExcessGPD(data, gamma, sigma, R, L = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)ExcessGPD(data, gamma, sigma, R, L = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)
data |
Vector of |
gamma |
Vector of |
sigma |
Vector of |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance, default is |
warnings |
Logical indicating if warnings are displayed, default is |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
We need that , the -th largest observation.
If this is not the case, we return NA for the premium. A warning will be issued in
that case if warnings=TRUE. One should then use global fits: ExcessSplice.
The premium for the excess-loss insurance with retention and limit is given by
where is the premium of the excess-loss insurance with retention . When , the premium is equal to .
We estimate by
with and the estimates for the parameters of the GPD.
See Section 4.6 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
premium |
The corresponding estimates for the premium. |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
data(secura) # GPDmle estimator mle <- GPDmle(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessGPD(secura$size, gamma=mle$gamma, sigma=mle$sigma, R=R, ylim=c(0,2*10^4))data(secura) # GPDmle estimator mle <- GPDmle(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessGPD(secura$size, gamma=mle$gamma, sigma=mle$sigma, R=R, ylim=c(0,2*10^4))
Estimate premiums of excess-loss reinsurance with retention and limit using a (truncated) Pareto model.
ExcessPareto(data, gamma, R, L = Inf, endpoint = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...) ExcessHill(data, gamma, R, L = Inf, endpoint = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)ExcessPareto(data, gamma, R, L = Inf, endpoint = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...) ExcessHill(data, gamma, R, L = Inf, endpoint = Inf, warnings = TRUE, plot = TRUE, add = FALSE, main = "Estimates for premium of excess-loss insurance", ...)
data |
Vector of |
gamma |
Vector of |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance, default is |
endpoint |
Endpoint for the truncated Pareto distribution. When |
warnings |
Logical indicating if warnings are displayed, default is |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
We need that , the -th largest observation.
If this is not the case, we return NA for the premium. A warning will be issued in
that case if warnings=TRUE. One should then use global fits: ExcessSplice.
The premium for the excess-loss insurance with retention and limit is given by
where is the premium of the excess-loss insurance with retention . When , the premium is equal to .
We estimate (for the untruncated Pareto distribution) by
with the Hill estimator.
The ExcessHill function is the same function but with a different name for compatibility with old versions of the package.
See Section 4.6 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
premium |
The corresponding estimates for the premium. |
R |
The retention level of the (re-)insurance. |
L |
The limit of the (re-)insurance. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Hill, ExcessEPD, ExcessGPD, ExcessSplice
data(secura) # Hill estimator H <- Hill(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessPareto(secura$size, H$gamma, R=R)data(secura) # Hill estimator H <- Hill(secura$size) # Premium of excess-loss insurance with retention R R <- 10^7 ExcessPareto(secura$size, H$gamma, R=R)
Estimate premiums of excess-loss reinsurance with retention and limit using fitted spliced distribution.
ExcessSplice(R, L=Inf, splicefit)ExcessSplice(R, L=Inf, splicefit)
R |
The retention level of the (re-)insurance or a vector of retention levels for the (re-)insurance. |
L |
The limit for the (re-)insurance or a vector of limits for the (re-)insurance, default is |
splicefit |
A |
The premium for the excess-loss insurance with retention and limit is given by
where is the premium of the excess-loss insurance with retention . When , the premium is equal to .
See Reynkens et al. (2017) and Section 4.6 of Albrecher et al. (2017) for more details.
An estimate for the premium is returned (for every value of R).
Tom Reynkens with R code from Roel Verbelen for the estimates for the excess-loss premiums using the mixed Erlang distribution.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.8) # Excess-loss premium ExcessSplice(R=2, splicefit=splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.8) # Excess-loss premium ExcessSplice(R=2, splicefit=splicefit) ## End(Not run)
Computes the empirical quantiles of a data vector and the theoretical quantiles of the standard exponential distribution. These quantiles are then plotted in an exponential QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
ExpQQ(data, plot = TRUE, main = "Exponential QQ-plot", ...)ExpQQ(data, plot = TRUE, main = "Exponential QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the quantiles should be plotted in an Exponential QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The exponential QQ-plot is defined as
for
with the -th order statistic of the data.
Note that the mean excess plot is the derivative plot of the Exponential QQ-plot.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
eqq.the |
Vector of the theoretical quantiles from a standard exponential distribution. |
eqq.emp |
Vector of the empirical quantiles from the data. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
MeanExcess, LognormalQQ, ParetoQQ, WeibullQQ
data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])
Density, distribution function, quantile function and random generation for the Extended Pareto Distribution (EPD).
depd(x, gamma, kappa, tau = -1, log = FALSE) pepd(x, gamma, kappa, tau = -1, lower.tail = TRUE, log.p = FALSE) qepd(p, gamma, kappa, tau = -1, lower.tail = TRUE, log.p = FALSE) repd(n, gamma, kappa, tau = -1)depd(x, gamma, kappa, tau = -1, log = FALSE) pepd(x, gamma, kappa, tau = -1, lower.tail = TRUE, log.p = FALSE) qepd(p, gamma, kappa, tau = -1, lower.tail = TRUE, log.p = FALSE) repd(n, gamma, kappa, tau = -1)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
gamma |
The |
kappa |
The |
tau |
The |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the EPD is equal to
for all and otherwise.
Note that an EPD random variable with and is GPD distributed with , and .
depd gives the density function evaluated in , pepd the CDF evaluated in and qepd the quantile function evaluated in . The length of the result is equal to the length of or .
repd returns a random sample of length .
Tom Reynkens.
Beirlant, J., Joossens, E. and Segers, J. (2009). "Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions." Journal of Statistical Planning and Inference, 139, 2800–2815.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, depd(x, gamma=1/2, kappa=1, tau=-1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pepd(x, gamma=1/2, kappa=1, tau=-1), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, depd(x, gamma=1/2, kappa=1, tau=-1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pepd(x, gamma=1/2, kappa=1, tau=-1), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the Fréchet distribution (inverse Weibull distribution).
dfrechet(x, shape, loc = 0, scale = 1, log = FALSE) pfrechet(x, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) qfrechet(p, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) rfrechet(n, shape, loc = 0, scale = 1)dfrechet(x, shape, loc = 0, scale = 1, log = FALSE) pfrechet(x, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) qfrechet(p, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) rfrechet(n, shape, loc = 0, scale = 1)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
shape |
Shape parameter of the Fréchet distribution. |
loc |
Location parameter of the Fréchet distribution, default is 0. |
scale |
Scale parameter of the Fréchet distribution, default is 1. |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the Fréchet distribution is equal to
for all and otherwise. Both shape and scale need to be strictly positive.
dfrechet gives the density function evaluated in , pfrechet the CDF evaluated in and qfrechet the quantile function evaluated in . The length of the result is equal to the length of or .
rfrechet returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(1,10,0.01) plot(x, dfrechet(x, shape=2), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1,10,0.01) plot(x, pfrechet(x, shape=2), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(1,10,0.01) plot(x, dfrechet(x, shape=2), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1,10,0.01) plot(x, pfrechet(x, shape=2), xlab="x", ylab="CDF", type="l")
Computes the generalised Hill estimator for real extreme value indices as a function of the tail parameter .
Optionally, these estimates are plotted as a function of .
genHill(data, gamma, logk = FALSE, plot = FALSE, add = FALSE, main = "Generalised Hill estimates of the EVI", ...)genHill(data, gamma, logk = FALSE, plot = FALSE, add = FALSE, main = "Generalised Hill estimates of the EVI", ...)
data |
Vector of |
gamma |
Vector of |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 generalised Hill estimator is an estimator for the slope of the last points of the generalised QQ-plot:
with the UH scores and the Hill estimates.
This is analogous to the (ordinary) Hill estimator which is the estimator of the slope of the last points of the Pareto QQ-plot when using constrained least squares.
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding generalised Hill estimates. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Beirlant, J., Vynckier, P. and Teugels, J.L. (1996). "Excess Function and Estimation of the Extreme-value Index". Bernoulli, 2, 293–318.
data(soa) # Hill estimator H <- Hill(soa$size, plot=FALSE) # Moment estimator M <- Moment(soa$size) # Generalised Hill estimator gH <- genHill(soa$size, gamma=H$gamma) # Plot estimates plot(H$k[1:5000], M$gamma[1:5000], xlab="k", ylab=expression(gamma), type="l", ylim=c(0.2,0.5)) lines(H$k[1:5000], gH$gamma[1:5000], lty=2) legend("topright", c("Moment", "Generalised Hill"), lty=1:2)data(soa) # Hill estimator H <- Hill(soa$size, plot=FALSE) # Moment estimator M <- Moment(soa$size) # Generalised Hill estimator gH <- genHill(soa$size, gamma=H$gamma) # Plot estimates plot(H$k[1:5000], M$gamma[1:5000], xlab="k", ylab=expression(gamma), type="l", ylim=c(0.2,0.5)) lines(H$k[1:5000], gH$gamma[1:5000], lty=2) legend("topright", c("Moment", "Generalised Hill"), lty=1:2)
Computes the empirical quantiles of the UH scores of a data vector and the theoretical quantiles of the standard exponential distribution. These quantiles are then plotted in a generalised QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
genQQ(data, gamma, plot = TRUE, main = "Generalised QQ-plot", ...) generalizedQQ(data, gamma, plot = TRUE, main = "Generalised QQ-plot", ...)genQQ(data, gamma, plot = TRUE, main = "Generalised QQ-plot", ...) generalizedQQ(data, gamma, plot = TRUE, main = "Generalised QQ-plot", ...)
data |
Vector of |
gamma |
Vector of |
plot |
Logical indicating if the quantiles should be plotted in a generalised QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The generalizedQQ function is the same function but with a different name for compatibility with the old S-Plus code.
The UH scores are defined as with the Hill estimates, but other positive estimates for the EVI can also be used. The appropriate positive estimates for the EVI need to be specified in gamma. The generalised QQ-plot then plots
for .
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
gqq.the |
Vector of the theoretical quantiles from a standard exponential distribution. |
gqq.emp |
Vector of the empirical quantiles from the logarithm of the UH scores. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Beirlant, J., Vynckier, P. and Teugels, J.L. (1996). "Excess Function and Estimation of the Extreme-value Index." Bernoulli, 2, 293–318.
data(soa) # Compute Hill estimator H <- Hill(soa$size[1:5000], plot=FALSE)$gamma # Generalised QQ-plot genQQ(soa$size[1:5000], gamma=H)data(soa) # Compute Hill estimator H <- Hill(soa$size[1:5000], plot=FALSE)$gamma # Generalised QQ-plot genQQ(soa$size[1:5000], gamma=H)
Density, distribution function, quantile function and random generation for the Generalised Pareto Distribution (GPD).
dgpd(x, gamma, mu = 0, sigma, log = FALSE) pgpd(x, gamma, mu = 0, sigma, lower.tail = TRUE, log.p = FALSE) qgpd(p, gamma, mu = 0, sigma, lower.tail = TRUE, log.p = FALSE) rgpd(n, gamma, mu = 0, sigma)dgpd(x, gamma, mu = 0, sigma, log = FALSE) pgpd(x, gamma, mu = 0, sigma, lower.tail = TRUE, log.p = FALSE) qgpd(p, gamma, mu = 0, sigma, lower.tail = TRUE, log.p = FALSE) rgpd(n, gamma, mu = 0, sigma)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
gamma |
The |
mu |
The |
sigma |
The |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the GPD for is equal to
for all and otherwise.
When , the CDF is given by
for all and otherwise.
dgpd gives the density function evaluated in , pgpd the CDF evaluated in and qgpd the quantile function evaluated in . The length of the result is equal to the length of or .
rgpd returns a random sample of length .
Tom Reynkens.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
tGPD, Pareto, EPD, Distributions
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dgpd(x, gamma=1/2, sigma=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pgpd(x, gamma=1/2, sigma=5), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dgpd(x, gamma=1/2, sigma=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, pgpd(x, gamma=1/2, sigma=5), xlab="x", ylab="CDF", type="l")
Fit the Generalised Pareto Distribution (GPD) to data using Maximum Likelihood Estimation (MLE).
GPDfit(data, start = c(0.1, 1), warnings = FALSE)GPDfit(data, start = c(0.1, 1), warnings = FALSE)
data |
Vector of |
start |
Vector of length 2 containing the starting values for the optimisation. The first element
is the starting value for the estimator of |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
See Section 4.2.2 in Albrecher et al. (2017) for more details.
A vector with the MLE estimate for the parameter of the GPD as the first component and the MLE estimate for the parameter of the GPD as the second component.
Tom Reynkens based on S-Plus code from Yuri Goegebeur and R code from Klaus Herrmann.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Fit GPD to last 500 observations res <- GPDfit(SOAdata-sort(soa$size)[500])data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Fit GPD to last 500 observations res <- GPDfit(SOAdata-sort(soa$size)[500])
Fit the Generalised Pareto Distribution (GPD) to the exceedances (peaks) over a threshold using Maximum Likelihood Estimation (MLE). Optionally, these estimates are plotted as a function of .
GPDmle(data, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...) POT(data, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...)GPDmle(data, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...) POT(data, start = c(0.1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "POT estimates of the EVI", ...)
data |
Vector of |
start |
Vector of length 2 containing the starting values for the optimisation. The first element
is the starting value for the estimator of |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The POT function is the same function but with a different name for compatibility with the old S-Plus code.
For each value of k, we look at the exceedances over the th largest observation:
for , with the th largest observation and the sample size. The GPD is then fitted to these k exceedances using MLE which yields estimates for the parameters of the GPD: and .
See Section 4.2.2 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding MLE estimates for the |
sigma |
Vector of the corresponding MLE estimates for the |
Tom Reynkens based on S-Plus code from Yuri Goegebeur and R code from Klaus Herrmann.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Plot GPD-ML estimates as a function of k GPDmle(SOAdata, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Plot GPD-ML estimates as a function of k GPDmle(SOAdata, plot=TRUE)
Residual plot to check GPD fit for peaks over a threshold.
GPDresiduals(data, t, gamma, sigma, plot = TRUE, main = "GPD residual plot", ...)GPDresiduals(data, t, gamma, sigma, plot = TRUE, main = "GPD residual plot", ...)
data |
Vector of |
t |
The used threshold. |
gamma |
Estimate for the EVI obtained from |
sigma |
Estimate for |
plot |
Logical indicating if the residuals should be plotted, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
Consider the POT values and the transformed variable
when and
otherwise.
We can assess the goodness-of-fit of the GPD when modelling POT values by
constructing an exponential QQ-plot of the transformed variable since is standard exponentially distributed if follows the GPD.
See Section 4.2.2 in Albrecher et al. (2017) for more details.
A list with following components:
res.the |
Vector of the theoretical quantiles from a standard exponential distribution. |
res.emp |
Vector of the empirical quantiles of |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Plot POT-MLE estimates as a function of k pot <- GPDmle(SOAdata, plot=TRUE) # Residual plot k <- 200 GPDresiduals(SOAdata, sort(SOAdata)[length(SOAdata)-k], pot$gamma[k], pot$sigma[k])data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Plot POT-MLE estimates as a function of k pot <- GPDmle(SOAdata, plot=TRUE) # Residual plot k <- 200 GPDresiduals(SOAdata, sort(SOAdata)[length(SOAdata)-k], pot$gamma[k], pot$sigma[k])
Computes the Hill estimator for positive extreme value indices (Hill, 1975) as a function of the tail parameter .
Optionally, these estimates are plotted as a function of .
Hill(data, k = TRUE, logk = FALSE, plot = FALSE, add = FALSE, main = "Hill estimates of the EVI", ...)Hill(data, k = TRUE, logk = FALSE, plot = FALSE, add = FALSE, main = "Hill estimates of the EVI", ...)
data |
Vector of |
k |
Logical indicating if the Hill estimates are plotted as a function of the tail parameter |
logk |
Logical indicating if the Hill estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 Hill estimator can be seen as the estimator of slope in the upper right corner ( last points) of the Pareto QQ-plot when using constrained least squares (the regression line has to pass through the point ). It is given by
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding Hill estimates. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Hill, B. M. (1975). "A simple general approach to inference about the tail of a distribution." Annals of Statistics, 3, 1163–1173.
data(norwegianfire) # Plot Hill estimates as a function of k Hill(norwegianfire$size[norwegianfire$year==76],plot=TRUE)data(norwegianfire) # Plot Hill estimates as a function of k Hill(norwegianfire$size[norwegianfire$year==76],plot=TRUE)
Computes bias-reduced ML estimates of gamma based on the quantile view.
Hill.2oQV(data, start = c(1,1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of the EVI", ...)Hill.2oQV(data, start = c(1,1,1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of the EVI", ...)
data |
Vector of |
start |
A vector of length 3 containing starting values for the first numerical optimisation (see Details). The elements
are the starting values for the estimators of |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the ML estimates for the EVI for each value of |
b |
Vector of the ML estimates for the parameter |
beta |
Vector of the ML estimates for the parameter |
Tom Reynkens based on S-Plus code from Yuri Goegebeur and R code from Klaus Herrmann.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Dierckx, G., Goegebeur Y. and Matthys, G. (1999). "Tail Index Estimation and an Exponential Regression Model." Extremes, 2, 177–200.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(norwegianfire) # Plot bias-reduced MLE (QV) as a function of k Hill.2oQV(norwegianfire$size[norwegianfire$year==76],plot=TRUE)data(norwegianfire) # Plot bias-reduced MLE (QV) as a function of k Hill.2oQV(norwegianfire$size[norwegianfire$year==76],plot=TRUE)
Select optimal threshold for the Hill estimator by minimising the Asymptotic Mean Squared Error (AMSE).
Hill.kopt(data, start = c(1, 1, 1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "AMSE plot", ...)Hill.kopt(data, start = c(1, 1, 1), warnings = FALSE, logk = FALSE, plot = FALSE, add = FALSE, main = "AMSE plot", ...)
data |
Vector of |
start |
A vector of length 3 containing starting values for the first numerical optimisation (see |
warnings |
Logical indicating if possible warnings from the optimisation function are shown, default is |
logk |
Logical indicating if the AMSE values are plotted as a function of |
plot |
Logical indicating if the AMSE values should be plotted as a function of |
add |
Logical indicating if the optimal value for |
main |
Title for the plot, default is |
... |
Additional arguments for the |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
AMSE |
Vector of the AMSE values for each value of |
kopt |
Optimal value of |
gammaopt |
Optimal value of the Hill estimator corresponding to minimal |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Beirlant J., Vynckier, P. and Teugels, J. (1996). "Tail Index Estimation, Pareto Quantile Plots, and Regression Diagnostics." Journal of the American Statistical Association, 91, 1659–1667.
data(norwegianfire) # Plot Hill estimator as a function of k Hill(norwegianfire$size[norwegianfire$year==76],plot=TRUE) # Add optimal value of k Hill.kopt(norwegianfire$size[norwegianfire$year==76],add=TRUE)data(norwegianfire) # Plot Hill estimator as a function of k Hill(norwegianfire$size[norwegianfire$year==76],plot=TRUE) # Add optimal value of k Hill.kopt(norwegianfire$size[norwegianfire$year==76],add=TRUE)
Computes the Hill estimator for positive extreme value indices, adapted for interval censoring, as a function of the tail parameter . Optionally, these estimates are plotted as a function of .
icHill(L, U, censored, trunclower = 0, truncupper = Inf, logk = FALSE, plot = TRUE, add = FALSE, main = "Hill estimates of the EVI", ...)icHill(L, U, censored, trunclower = 0, truncupper = Inf, logk = FALSE, plot = TRUE, add = FALSE, main = "Hill estimates of the EVI", ...)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
trunclower |
Lower truncation point. Default is 0. |
truncupper |
Upper truncation point. Default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
This estimator is given by
where is the Turnbull estimator for the CDF.
More specifically, we use the values for where
is the empirical quantile function corresponding to the Turnbull estimator.
We then denote
with
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u.
If the interval package is installed, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Use Hill for non-censored data or cHill for right censored data.
See Section 4.3 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding Hill estimates. |
X |
Vector of thresholds |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
cHill, Hill, MeanExcess_TB, icParetoQQ, Turnbull, icfit
# 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) # Right boundary U <- Z U[censored] <- Inf # Hill estimator adapted for interval censoring icHill(Z, U, censored, ylim=c(0,1)) # Hill estimator adapted for right censoring cHill(Z, censored, lty=2, add=TRUE) # True value of gamma abline(h=1/2, lty=3, col="blue") # Legend legend("topright", c("icHill", "cHill"), lty=1:2)# 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) # Right boundary U <- Z U[censored] <- Inf # Hill estimator adapted for interval censoring icHill(Z, U, censored, ylim=c(0,1)) # Hill estimator adapted for right censoring cHill(Z, censored, lty=2, add=TRUE) # True value of gamma abline(h=1/2, lty=3, col="blue") # Legend legend("topright", c("icHill", "cHill"), lty=1:2)
Pareto QQ-plot adapted for interval censored data using the Turnbull estimator.
icParetoQQ(L, U = L, censored, trunclower = 0, truncupper = Inf, plot = TRUE, main = "Pareto QQ-plot", ...)icParetoQQ(L, U = L, censored, trunclower = 0, truncupper = Inf, plot = TRUE, main = "Pareto QQ-plot", ...)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
trunclower |
Lower truncation point. Default is 0. |
truncupper |
Upper truncation point. Default is |
plot |
Logical indicating if the quantiles should be plotted in a Pareto QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Pareto QQ-plot adapted for interval censoring is given by
for where is the Turnbull estimator for the CDF and with the empirical quantile function corresponding to the Turnbull estimator.
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u.
If the interval package is installed, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Use ParetoQQ for non-censored data or cParetoQQ for right censored data.
See Section 4.3 in Albrecher et al. (2017) for more details.
A list with following components:
pqq.the |
Vector of the theoretical quantiles, see Details. |
pqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
cParetoQQ, ParetoQQ, icHill, Turnbull, icfit
# 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) # Right boundary U <- Z U[censored] <- Inf # Pareto QQ-plot adapted for interval censoring icParetoQQ(Z, U, censored) # Pareto QQ-plot adapted for right censoring cParetoQQ(Z, censored)# 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) # Right boundary U <- Z U[censored] <- Inf # Pareto QQ-plot adapted for interval censoring icParetoQQ(Z, U, censored) # Pareto QQ-plot adapted for right censoring cParetoQQ(Z, censored)
Computes the Kaplan-Meier estimator for the survival function of right censored data.
KaplanMeier(x, data, censored, conf.type="plain", conf.int = 0.95)KaplanMeier(x, data, censored, conf.type="plain", conf.int = 0.95)
x |
Vector with points to evaluate the estimator in. |
data |
Vector of |
censored |
Vector of |
conf.type |
Type of confidence interval, see |
conf.int |
Confidence level of the two-sided confidence interval, see |
We consider the random right censoring model where one observes
where is the variable of interest and is the censoring variable.
This function is merely a wrapper for survfit.formula from survival.
This estimator is only suitable for right censored data. When the data are interval censored, one can use the Turnbull estimator implemented in Turnbull.
A list with following components:
surv |
A vector of length |
fit |
The output from the call to |
Tom Reynkens
Kaplan, E. L. and Meier, P. (1958). "Nonparametric Estimation from Incomplete Observations." Journal of the American Statistical Association, 53, 457–481.
data <- c(1, 2.5, 3, 4, 5.5, 6, 7.5, 8.25, 9, 10.5) censored <- c(0, 1, 0, 0, 1, 0, 1, 1, 0, 1) x <- seq(0, 12, 0.1) # Kaplan-Meier estimator plot(x, KaplanMeier(x, data, censored)$surv, type="s", ylab="Kaplan-Meier estimator")data <- c(1, 2.5, 3, 4, 5.5, 6, 7.5, 8.25, 9, 10.5) censored <- c(0, 1, 0, 0, 1, 0, 1, 1, 0, 1) x <- seq(0, 12, 0.1) # Kaplan-Meier estimator plot(x, KaplanMeier(x, data, censored)$surv, type="s", ylab="Kaplan-Meier estimator")
Computes the empirical quantiles of the log-transform of a data vector and the theoretical quantiles of the standard normal distribution. These quantiles are then plotted in a log-normal QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
LognormalQQ(data, plot = TRUE, main = "Log-normal QQ-plot", ...)LognormalQQ(data, plot = TRUE, main = "Log-normal QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the quantiles should be plotted in a log-normal QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
By definition, a log-transformed log-normal random variable is normally distributed. We can thus obtain a log-normal QQ-plot from a normal QQ-plot by replacing the empirical quantiles of the data vector by the empirical quantiles from the log-transformed data. We hence plot
for where is the standard normal CDF.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
lnqq.the |
Vector of the theoretical quantiles from a standard normal distribution. |
lnqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(norwegianfire) # Log-normal QQ-plot for Norwegian Fire Insurance data for claims in 1976. LognormalQQ(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Log-normal QQ-plot for Norwegian Fire Insurance data for claims in 1976. LognormalQQ(norwegianfire$size[norwegianfire$year==76])
Computes the derivative plot of the log-normal QQ-plot. These values can be plotted as a function of the data or as a function of the tail parameter .
LognormalQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of log-normal QQ-plot", ...)LognormalQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of log-normal QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the derivative values should be plotted, default is |
k |
Logical indicating if the derivative values are plotted as a function of the tail parameter |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The derivative plot of a log-normal QQ-plot is
or
with the Hill estimates and
Here is , the standard normal PDF and the standard normal CDF.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
xval |
Vector of the x-values of the plot ( |
yval |
Vector of the derivative values. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
LognormalQQ, Hill, MeanExcess, ParetoQQ_der, WeibullQQ_der
data(norwegianfire) # Log-normal QQ-plot for Norwegian Fire Insurance data for claims in 1976. LognormalQQ(norwegianfire$size[norwegianfire$year==76]) # Derivate plot LognormalQQ_der(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Log-normal QQ-plot for Norwegian Fire Insurance data for claims in 1976. LognormalQQ(norwegianfire$size[norwegianfire$year==76]) # Derivate plot LognormalQQ_der(norwegianfire$size[norwegianfire$year==76])
Computes the Least Squares (LS) estimates of the EVI based on the last observations of the generalised QQ-plot.
LStail(data, rho = -1, lambda = 0.5, logk = FALSE, plot = FALSE, add = FALSE, main = "LS estimates of the EVI", ...) TSfraction(data, rho = -1, lambda = 0.5, logk = FALSE, plot = FALSE, add = FALSE, main = "LS estimates of the EVI", ...)LStail(data, rho = -1, lambda = 0.5, logk = FALSE, plot = FALSE, add = FALSE, main = "LS estimates of the EVI", ...) TSfraction(data, rho = -1, lambda = 0.5, logk = FALSE, plot = FALSE, add = FALSE, main = "LS estimates of the EVI", ...)
data |
Vector of |
rho |
Estimate for |
lambda |
Parameter used in the method of Beirlant et al. (2002), only used when |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
We estimate (EVI) and using least squares on the following regression model (Beirlant et al., 2005): with and , where is the Hill estimator with threshold .
See Section 5.8 of Beirlant et al. (2004) for more details.
The function TSfraction is included for compatibility with the old S-Plus code.
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding LS estimates for the EVI. |
b |
Vector of the corresponding LS estimates for b. |
rho |
Vector of the estimates for |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Beirlant, J., Dierckx, G. and Guillou, A. (2005). "Estimation of the Extreme Value Index and Regression on Generalized Quantile Plots." Bernoulli, 11, 949–970.
Beirlant, J., Dierckx, G., Guillou, A. and Starica, C. (2002). "On Exponential Representations of Log-spacing of Extreme Order Statistics." Extremes, 5, 157–180.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # LS tail estimator LStail(soa$size, plot=TRUE, ylim=c(0,0.5))data(soa) # LS tail estimator LStail(soa$size, plot=TRUE, ylim=c(0,0.5))
Computes the mean excess values for a vector of observations. These mean excess values can then be plotted as a function of the data or as a function of the tail parameter .
MeanExcess(data, plot = TRUE, k = FALSE, main = "Mean excess plot", ...)MeanExcess(data, plot = TRUE, k = FALSE, main = "Mean excess plot", ...)
data |
Vector of |
plot |
Logical indicating if the mean excess values should be plotted in a mean excess plot, default is |
k |
Logical indicating if the mean excess scores are plotted as a function of the tail parameter |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The mean excess plot is
or
with
Note that the mean excess plot is the derivative plot of the Exponential QQ-plot.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
X |
Vector of the order statistics |
e |
Vector of the mean excess values corresponding to the tail parameters in |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
ExpQQ, LognormalQQ_der, ParetoQQ_der, WeibullQQ_der
data(norwegianfire) # Mean excess plots for Norwegian Fire Insurance data for claims in 1976. # Mean excess values as a function of k MeanExcess(norwegianfire$size[norwegianfire$year==76], k=TRUE) # Mean excess values as a function of the data MeanExcess(norwegianfire$size[norwegianfire$year==76], k=FALSE)data(norwegianfire) # Mean excess plots for Norwegian Fire Insurance data for claims in 1976. # Mean excess values as a function of k MeanExcess(norwegianfire$size[norwegianfire$year==76], k=TRUE) # Mean excess values as a function of the data MeanExcess(norwegianfire$size[norwegianfire$year==76], k=FALSE)
Computes mean excess values using the Turnbull estimator. These mean excess values can then be plotted as a function of the empirical quantiles (computed using the Turnbull estimator) or as a function of the tail parameter .
MeanExcess_TB(L, U = L, censored, trunclower = 0, truncupper = Inf, plot = TRUE, k = FALSE, intervalpkg = TRUE, main = "Mean excess plot", ...)MeanExcess_TB(L, U = L, censored, trunclower = 0, truncupper = Inf, plot = TRUE, k = FALSE, intervalpkg = TRUE, main = "Mean excess plot", ...)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
trunclower |
Lower truncation point, default is 0. |
truncupper |
Upper truncation point, default is |
plot |
Logical indicating if the mean excess values should be plotted in a mean excess plot, default is |
k |
Logical indicating if the mean excess values are plotted as a function of the tail parameter |
intervalpkg |
Logical indicating if the Turnbull estimator is computed using the implementation in the interval package if this package is installed. Default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The mean excess values are given by
where is the Turnbull estimator for the CDF.
More specifically, we use the values for where
is the empirical quantile function corresponding to the Turnbull estimator.
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u.
If the interval package is installed and intervalpkg=TRUE, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Use MeanExcess for non-censored data.
See Section 4.3 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
X |
Vector of the empirical quantiles, computed using the Turnbull estimator, corresponding to |
e |
Vector of the mean excess values corresponding to the tail parameters in |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
# 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) # Right boundary U <- Z U[censored] <- Inf # Mean excess plot MeanExcess_TB(Z, U, censored, k=FALSE)# 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) # Right boundary U <- Z U[censored] <- Inf # Mean excess plot MeanExcess_TB(Z, U, censored, k=FALSE)
Create an S3 object using a Mixed Erlang (ME) fit.
MEfit(p, shape, theta, M, M_initial = NULL)MEfit(p, shape, theta, M, M_initial = NULL)
p |
Vector of mixing weights, denoted by |
shape |
Vector of shape parameters |
theta |
Scale parameter |
M |
Number of mixture components. |
M_initial |
Initial value provided for |
The rate parameter used in Albrecher et al. (2017) is equal to .
See Reynkens et al. (2017) and Section 4.3 of Albrecher et al. (2017) for more details
An S3 object which contains the input arguments in a list.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)
Compute the moment estimates for real extreme value indices as a function of the tail parameter . Optionally, these estimates are plotted as a function of .
Moment(data, logk = FALSE, plot = FALSE, add = FALSE, main = "Moment estimates of the EVI", ...)Moment(data, logk = FALSE, plot = FALSE, add = FALSE, main = "Moment estimates of the EVI", ...)
data |
Vector of |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 moment estimator for the EVI is introduced by Dekkers et al. (1989) and is a generalisation of the Hill estimator.
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding moment estimates. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Dekkers, A.L.M, Einmahl, J.H.J. and de Haan, L. (1989). "A Moment Estimator for the Index of an Extreme-value Distribution." Annals of Statistics, 17, 1833–1855.
data(soa) # Hill estimator H <- Hill(soa$size, plot=FALSE) # Moment estimator M <- Moment(soa$size) # Generalised Hill estimator gH <- genHill(soa$size, gamma=H$gamma) # Plot estimates plot(H$k[1:5000], M$gamma[1:5000], xlab="k", ylab=expression(gamma), type="l", ylim=c(0.2,0.5)) lines(H$k[1:5000], gH$gamma[1:5000], lty=2) legend("topright", c("Moment", "Generalised Hill"), lty=1:2)data(soa) # Hill estimator H <- Hill(soa$size, plot=FALSE) # Moment estimator M <- Moment(soa$size) # Generalised Hill estimator gH <- genHill(soa$size, gamma=H$gamma) # Plot estimates plot(H$k[1:5000], M$gamma[1:5000], xlab="k", ylab=expression(gamma), type="l", ylim=c(0.2,0.5)) lines(H$k[1:5000], gH$gamma[1:5000], lty=2) legend("topright", c("Moment", "Generalised Hill"), lty=1:2)
Fire insurance claims for a Norwegian insurance company for the period 1972 to 1992 as studied in Beirlant et al. (1996).
A priority of 500 units was in force.
data("norwegianfire")data("norwegianfire")
A data frame with 9181 observations on the following 2 variables:
sizeSize of fire insurance claim (in 1000 NOK).
yearYear of claim occurence (expressed as instead of ).
Beirlant, J., Teugels, J. L. and Vynckier, P. (1996). Practical Analysis of Extreme Values, Leuven University Press.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])
Density, distribution function, quantile function and random generation for the Pareto distribution (type I).
dpareto(x, shape, scale = 1, log = FALSE) ppareto(x, shape, scale = 1, lower.tail = TRUE, log.p = FALSE) qpareto(p, shape, scale = 1, lower.tail = TRUE, log.p = FALSE) rpareto(n, shape, scale = 1)dpareto(x, shape, scale = 1, log = FALSE) ppareto(x, shape, scale = 1, lower.tail = TRUE, log.p = FALSE) qpareto(p, shape, scale = 1, lower.tail = TRUE, log.p = FALSE) rpareto(n, shape, scale = 1)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
shape |
The shape parameter of the Pareto distribution, a strictly positive number. |
scale |
The scale parameter of the Pareto distribution, a strictly positive number. Its default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the Pareto distribution is equal to
for all and otherwise. Both shape and scale need to be strictly positive.
dpareto gives the density function evaluated in , ppareto the CDF evaluated in and qpareto the quantile function evaluated in . The length of the result is equal to the length of or .
rpareto returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(1, 10, 0.01) plot(x, dpareto(x, shape=2), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1, 10, 0.01) plot(x, ppareto(x, shape=2), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(1, 10, 0.01) plot(x, dpareto(x, shape=2), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1, 10, 0.01) plot(x, ppareto(x, shape=2), xlab="x", ylab="CDF", type="l")
Computes the empirical quantiles of the log-transform of a data vector and the theoretical quantiles of the standard exponential distribution. These quantiles are then plotted in a Pareto QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
ParetoQQ(data, plot = TRUE, main = "Pareto QQ-plot", ...)ParetoQQ(data, plot = TRUE, main = "Pareto QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the quantiles should be plotted in a Pareto QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
It can be easily seen that a log-transformed Pareto random variable is exponentially distributed. We can hence obtain a Pareto QQ-plot from an exponential QQ-plot by replacing the empirical quantiles from the data vector by the empirical quantiles from the log-transformed data. We hence plot
for
with the -th order statistic of the data.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
pqq.the |
Vector of the theoretical quantiles from a standard exponential distribution. |
pqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
ParetoQQ_der, ExpQQ, genQQ, LognormalQQ, WeibullQQ
data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Exponential QQ-plot for Norwegian Fire Insurance data for claims in 1976. ExpQQ(norwegianfire$size[norwegianfire$year==76]) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76])
Computes the derivative plot of the Pareto QQ-plot. These values can be plotted as a function of the data or as a function of the tail parameter .
ParetoQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of Pareto QQ-plot", ...)ParetoQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of Pareto QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the derivative values should be plotted, default is |
k |
Logical indicating if the derivative values are plotted as a function of the tail parameter |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The derivative plot of a Pareto QQ-plot is
or
with the Hill estimates.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
xval |
Vector of the x-values of the plot ( |
yval |
Vector of the derivative values |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
ParetoQQ, Hill, MeanExcess, LognormalQQ_der, WeibullQQ_der
data(norwegianfire) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76]) # Derivate plot ParetoQQ_der(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Pareto QQ-plot for Norwegian Fire Insurance data for claims in 1976. ParetoQQ(norwegianfire$size[norwegianfire$year==76]) # Derivate plot ParetoQQ_der(norwegianfire$size[norwegianfire$year==76])
Compute approximations of the CDF using the normal approximation, normal-power approximation, shifted Gamma approximation or normal approximation to the shifted Gamma distribution.
pClas(x, mean = 0, variance = 1, skewness = NULL, method = c("normal", "normal-power", "shifted Gamma", "shifted Gamma normal"), lower.tail = TRUE, log.p = FALSE)pClas(x, mean = 0, variance = 1, skewness = NULL, method = c("normal", "normal-power", "shifted Gamma", "shifted Gamma normal"), lower.tail = TRUE, log.p = FALSE)
x |
Vector of points to approximate the CDF in. |
mean |
Mean of the distribution, default is 0. |
variance |
Variance of the distribution, default is 1. |
skewness |
Skewness coefficient of the distribution, this argument is not used for the normal approximation. Default is |
method |
Approximation method to use, one of |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The normal approximation for the CDF of the r.v. is defined as
where and are the mean and variance of , respectively.
This approximation can be improved when the skewness parameter
is available. The normal-power approximation of the CDF is then given by
for and .
The shifted Gamma approximation uses the approximation
Here, we need that .
The normal approximation to the shifted Gamma distribution approximates the CDF of as
for . We need again that .
See Section 6.2 of Albrecher et al. (2017) for more details.
Vector of estimates for the probabilities .
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Classical approximations p1 <- pClas(x, mean(X), var(X)) p2 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="normal-power") p3 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="shifted Gamma") p4 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="shifted Gamma normal") # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3, col="green") lines(x, p3, lty=4) lines(x, p4, lty=5, col="blue") legend("bottomright", c("True CDF", "normal approximation", "normal-power approximation", "shifted Gamma approximation", "shifted Gamma normal approximation"), lty=1:5, col=c("red", "black", "green", "black", "blue"), lwd=2)# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Classical approximations p1 <- pClas(x, mean(X), var(X)) p2 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="normal-power") p3 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="shifted Gamma") p4 <- pClas(x, mean(X), var(X), mean((X-mean(X))^3)/sd(X)^3, method="shifted Gamma normal") # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3, col="green") lines(x, p3, lty=4) lines(x, p4, lty=5, col="blue") legend("bottomright", c("True CDF", "normal approximation", "normal-power approximation", "shifted Gamma approximation", "shifted Gamma normal approximation"), lty=1:5, col=c("red", "black", "green", "black", "blue"), lwd=2)
Edgeworth approximation of the CDF using the first four moments.
pEdge(x, moments = c(0, 1, 0, 3), raw = TRUE, lower.tail = TRUE, log.p = FALSE)pEdge(x, moments = c(0, 1, 0, 3), raw = TRUE, lower.tail = TRUE, log.p = FALSE)
x |
Vector of points to approximate the CDF in. |
moments |
The first four raw moments if |
raw |
When |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
Denote the standard normal PDF and CDF respectively by and .
Let be the first moment, the variance, the third central moment and the fourth central moment of the random variable .
The corresponding cumulants are given by , , and .
Now consider the random variable , which has cumulants
0, 1, and .
The Edgeworth approximation for the CDF of () is given by
with , , and .
See Section 6.2 of Albrecher et al. (2017) for more details.
Vector of estimates for the probabilities .
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Cheah, P.K., Fraser, D.A.S. and Reid, N. (1993). "Some Alternatives to Edgeworth." The Canadian Journal of Statistics, 21(2), 131–138.
# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Empirical moments moments = c(mean(X), mean(X^2), mean(X^3), mean(X^4)) # Gram-Charlier approximation p1 <- pGC(x, moments) # Edgeworth approximation p2 <- pEdge(x, moments) # Normal approximation p3 <- pClas(x, mean(X), var(X)) # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3) lines(x, p3, lty=4, col="blue") legend("bottomright", c("True CDF", "GC approximation", "Edgeworth approximation", "Normal approximation"), col=c("red", "black", "black", "blue"), lty=1:4, lwd=2)# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Empirical moments moments = c(mean(X), mean(X^2), mean(X^3), mean(X^4)) # Gram-Charlier approximation p1 <- pGC(x, moments) # Edgeworth approximation p2 <- pEdge(x, moments) # Normal approximation p3 <- pClas(x, mean(X), var(X)) # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3) lines(x, p3, lty=4, col="blue") legend("bottomright", c("True CDF", "GC approximation", "Edgeworth approximation", "Normal approximation"), col=c("red", "black", "black", "blue"), lty=1:4, lwd=2)
Gram-Charlier approximation of the CDF using the first four moments.
pGC(x, moments = c(0, 1, 0, 3), raw = TRUE, lower.tail = TRUE, log.p = FALSE)pGC(x, moments = c(0, 1, 0, 3), raw = TRUE, lower.tail = TRUE, log.p = FALSE)
x |
Vector of points to approximate the CDF in. |
moments |
The first four raw moments if |
raw |
When |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
Denote the standard normal PDF and CDF respectively by and .
Let be the first moment, the variance, the third central moment and the fourth central moment of the random variable .
The corresponding cumulants are given by , , and .
Now consider the random variable , which has cumulants
0, 1, and .
The Gram-Charlier approximation for the CDF of () is given by
with , and .
See Section 6.2 of Albrecher et al. (2017) for more details.
Vector of estimates for the probabilities .
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Cheah, P.K., Fraser, D.A.S. and Reid, N. (1993). "Some Alternatives to Edgeworth." The Canadian Journal of Statistics, 21(2), 131–138.
# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Empirical moments moments = c(mean(X), mean(X^2), mean(X^3), mean(X^4)) # Gram-Charlier approximation p1 <- pGC(x, moments) # Edgeworth approximation p2 <- pEdge(x, moments) # Normal approximation p3 <- pClas(x, mean(X), var(X)) # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3) lines(x, p3, lty=4, col="blue") legend("bottomright", c("True CDF", "GC approximation", "Edgeworth approximation", "Normal approximation"), col=c("red", "black", "black", "blue"), lty=1:4, lwd=2)# Chi-squared sample X <- rchisq(1000, 2) x <- seq(0, 10, 0.01) # Empirical moments moments = c(mean(X), mean(X^2), mean(X^3), mean(X^4)) # Gram-Charlier approximation p1 <- pGC(x, moments) # Edgeworth approximation p2 <- pEdge(x, moments) # Normal approximation p3 <- pClas(x, mean(X), var(X)) # True probabilities p <- pchisq(x, 2) # Plot true and estimated probabilities plot(x, p, type="l", ylab="F(x)", ylim=c(0,1), col="red") lines(x, p1, lty=2) lines(x, p2, lty=3) lines(x, p3, lty=4, col="blue") legend("bottomright", c("True CDF", "GC approximation", "Edgeworth approximation", "Normal approximation"), col=c("red", "black", "black", "blue"), lty=1:4, lwd=2)
Compute estimates of a small exceedance probability or large return period using the approach of Weissman (1978).
Prob(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) Return(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...) Weissman.p(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) Weissman.r(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)Prob(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) Return(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...) Weissman.p(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) Weissman.r(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
gamma |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
Weissman.p and Weissman.r are the same functions as Prob and Return but with a different name for compatibility with the old S-Plus code.
A list with following components:
k |
Vector of the values of the tail parameter |
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 based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Weissman, I. (1978). "Estimation of Parameters and Large Quantiles Based on the k Largest Observations." Journal of the American Statistical Association, 73, 812–815.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Exceedance probability q <- 10^6 # Weissman estimator Prob(SOAdata,gamma=H$gamma,q=q,plot=TRUE) # Return period q <- 10^6 # Weissman estimator Return(SOAdata,gamma=H$gamma,q=q,plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Exceedance probability q <- 10^6 # Weissman estimator Prob(SOAdata,gamma=H$gamma,q=q,plot=TRUE) # Return period q <- 10^6 # Weissman estimator Return(SOAdata,gamma=H$gamma,q=q,plot=TRUE)
Computes estimates of a small exceedance probability or large return period using the parameters from the EPD fit.
ProbEPD(data, q, gamma, kappa, tau, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnEPD(data, q, gamma, kappa, tau, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)ProbEPD(data, q, gamma, kappa, tau, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnEPD(data, q, gamma, kappa, tau, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
q |
The used large quantile (we estimate |
gamma |
Vector of |
kappa |
Vector of |
tau |
Vector of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
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
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Joossens, E. and Segers, J. (2009). "Second-Order Refined Peaks-Over-Threshold Modelling for Heavy-Tailed Distributions." Journal of Statistical Planning and Inference, 139, 2800–2815.
data(secura) # EPD estimates for the EVI epd <- EPD(secura$size, plot=TRUE) # Compute exceedance probabilities q <- 10^7 ProbEPD(secura$size, q=q, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE) # Compute return periods ReturnEPD(secura$size, q=q, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE, ylim=c(0,10^4))data(secura) # EPD estimates for the EVI epd <- EPD(secura$size, plot=TRUE) # Compute exceedance probabilities q <- 10^7 ProbEPD(secura$size, q=q, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE) # Compute return periods ReturnEPD(secura$size, q=q, gamma=epd$gamma, kappa=epd$kappa, tau=epd$tau, plot=TRUE, ylim=c(0,10^4))
Computes estimates of a small exceedance probability or large return period using the generalised Hill estimates for the EVI.
ProbGH(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnGH(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)ProbGH(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnGH(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
gamma |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Beirlant, J., Vynckier, P. and Teugels, J.L. (1996). "Excess Function and Estimation of the Extreme-value Index". Bernoulli, 2, 293–318.
QuantGH, genHill, ProbMOM, Prob
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Generalised Hill estimator gH <- genHill(SOAdata, H$gamma) # Exceedance probability q <- 10^7 ProbGH(SOAdata, gamma=gH$gamma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnGH(SOAdata, gamma=gH$gamma, q=q, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Generalised Hill estimator gH <- genHill(SOAdata, H$gamma) # Exceedance probability q <- 10^7 ProbGH(SOAdata, gamma=gH$gamma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnGH(SOAdata, gamma=gH$gamma, q=q, plot=TRUE)
Computes estimates of a small exceedance probability or large return period using the GPD fit for the peaks over a threshold.
ProbGPD(data, gamma, sigma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnGPD(data, gamma, sigma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)ProbGPD(data, gamma, sigma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnGPD(data, gamma, sigma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
gamma |
Vector of |
sigma |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # GPD-ML estimator pot <- GPDmle(SOAdata) # Exceedance probability q <- 10^7 ProbGPD(SOAdata, gamma=pot$gamma, sigma=pot$sigma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnGPD(SOAdata, gamma=pot$gamma, sigma=pot$sigma, q=q, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # GPD-ML estimator pot <- GPDmle(SOAdata) # Exceedance probability q <- 10^7 ProbGPD(SOAdata, gamma=pot$gamma, sigma=pot$sigma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnGPD(SOAdata, gamma=pot$gamma, sigma=pot$sigma, q=q, plot=TRUE)
Computes estimates of a small exceedance probability or large return period using the Method of Moments estimates for the EVI.
ProbMOM(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnMOM(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)ProbMOM(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...) ReturnMOM(data, gamma, q, plot = FALSE, add = FALSE, main = "Estimates of large return period", ...)
data |
Vector of |
gamma |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
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.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Dekkers, A.L.M, Einmahl, J.H.J. and de Haan, L. (1989). "A Moment Estimator for the Index of an Extreme-value Distribution." Annals of Statistics, 17, 1833–1855.
QuantMOM, Moment, ProbGH, Prob
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # MOM estimator M <- Moment(SOAdata) # Exceedance probability q <- 10^7 ProbMOM(SOAdata, gamma=M$gamma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnMOM(SOAdata, gamma=M$gamma, q=q, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # MOM estimator M <- Moment(SOAdata) # Exceedance probability q <- 10^7 ProbMOM(SOAdata, gamma=M$gamma, q=q, plot=TRUE) # Return period q <- 10^7 ReturnMOM(SOAdata, gamma=M$gamma, q=q, plot=TRUE)
Estimator of small tail probability in the regression case where is constant and the regression modelling is thus only solely placed on the scale parameter.
ProbReg(Z, A, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)ProbReg(Z, A, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)
Z |
Vector of |
A |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 estimator is defined as
with the Hill estimator. Here, it is assumed that we have equidistant covariates .
See Section 4.4.1 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
P |
Vector of the corresponding probability estimates. |
q |
The used large quantile. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)
Compute estimates of an extreme quantile using the approach of Weissman (1978).
Quant(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) Weissman.q(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)Quant(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) Weissman.q(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
add |
Logical indicating if the estimates as a function of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
Weissman.q is the same function but with a different name for compatibility with the old S-Plus code.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Weissman, I. (1978). "Estimation of Parameters and Large Quantiles Based on the k Largest Observations." Journal of the American Statistical Association, 73, 812–815.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Bias-reduced estimator (QV) H_QV <- Hill.2oQV(SOAdata) # Exceedance probability p <- 10^(-5) # Weissman estimator Quant(SOAdata, gamma=H$gamma, p=p, plot=TRUE) # Second order Weissman estimator (QV) Quant.2oQV(SOAdata, gamma=H_QV$gamma, beta=H_QV$beta, b=H_QV$b, p=p, add=TRUE, lty=2)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Bias-reduced estimator (QV) H_QV <- Hill.2oQV(SOAdata) # Exceedance probability p <- 10^(-5) # Weissman estimator Quant(SOAdata, gamma=H$gamma, p=p, plot=TRUE) # Second order Weissman estimator (QV) Quant.2oQV(SOAdata, gamma=H_QV$gamma, beta=H_QV$beta, b=H_QV$b, p=p, add=TRUE, lty=2)
Compute second order refined Weissman estimator of extreme quantiles using the quantile view.
Quant.2oQV(data, gamma, b, beta, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) Weissman.q.2oQV(data, gamma, b, beta, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)Quant.2oQV(data, gamma, b, beta, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) Weissman.q.2oQV(data, gamma, b, beta, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
b |
Vector of |
beta |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.1 of Albrecher et al. (2017) for more details.
Weissman.q.2oQV is the same function but with a different name for compatibility with the old S-Plus code.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens based on S-Plus code from Yuri Goegebeur.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Bias-reduced estimator (QV) H_QV <- Hill.2oQV(SOAdata) # Exceedance probability p <- 10^(-5) # Weissman estimator Quant(SOAdata, gamma=H$gamma, p=p, plot=TRUE) # Second order Weissman estimator (QV) Quant.2oQV(SOAdata, gamma=H_QV$gamma, beta=H_QV$beta, b=H_QV$b, p=p, add=TRUE, lty=2)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Bias-reduced estimator (QV) H_QV <- Hill.2oQV(SOAdata) # Exceedance probability p <- 10^(-5) # Weissman estimator Quant(SOAdata, gamma=H$gamma, p=p, plot=TRUE) # Second order Weissman estimator (QV) Quant.2oQV(SOAdata, gamma=H_QV$gamma, beta=H_QV$beta, b=H_QV$b, p=p, add=TRUE, lty=2)
Compute estimates of an extreme quantile using generalised Hill estimates of the EVI.
QuantGH(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)QuantGH(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Beirlant, J., Vynckier, P. and Teugels, J.L. (1996). "Excess Function and Estimation of the Extreme-value Index". Bernoulli, 2, 293–318.
ProbGH, genHill, QuantMOM, Quant
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Generalised Hill estimator gH <- genHill(SOAdata, H$gamma) # Large quantile p <- 10^(-5) QuantGH(SOAdata, p=p, gamma=gH$gamma, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # Hill estimator H <- Hill(SOAdata) # Generalised Hill estimator gH <- genHill(SOAdata, H$gamma) # Large quantile p <- 10^(-5) QuantGH(SOAdata, p=p, gamma=gH$gamma, plot=TRUE)
Computes estimates of an extreme quantile using the GPD fit for the peaks over a threshold.
QuantGPD(data, gamma, sigma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)QuantGPD(data, gamma, sigma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
sigma |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # GPD-ML estimator pot <- GPDmle(SOAdata) # Large quantile p <- 10^(-5) QuantGPD(SOAdata, p=p, gamma=pot$gamma, sigma=pot$sigma, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # GPD-ML estimator pot <- GPDmle(SOAdata) # Large quantile p <- 10^(-5) QuantGPD(SOAdata, p=p, gamma=pot$gamma, sigma=pot$sigma, plot=TRUE)
Compute estimates of an extreme quantile using the Method of Moments estimates of the EVI.
QuantMOM(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)QuantMOM(data, gamma, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
See Section 4.2.2 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
Dekkers, A.L.M, Einmahl, J.H.J. and de Haan, L. (1989). "A Moment Estimator for the Index of an Extreme-value Distribution." Annals of Statistics, 17, 1833–1855.
ProbMOM, Moment, QuantGH, Quant
data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # MOM estimator M <- Moment(SOAdata) # Large quantile p <- 10^(-5) QuantMOM(SOAdata, p=p, gamma=M$gamma, plot=TRUE)data(soa) # Look at last 500 observations of SOA data SOAdata <- sort(soa$size)[length(soa$size)-(0:499)] # MOM estimator M <- Moment(SOAdata) # Large quantile p <- 10^(-5) QuantMOM(SOAdata, p=p, gamma=M$gamma, plot=TRUE)
Estimator of extreme quantile in the regression case where is constant and the regression modelling is thus only solely placed on the scale parameter.
QuantReg(Z, A, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)QuantReg(Z, A, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
Z |
Vector of |
A |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 estimator is defined as
with the Hill estimator. Here, it is assumed that we have equidistant covariates .
See Section 4.4.1 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)
Computes the estimator for the scale parameter as described in Beirlant et al. (2016).
Scale(data, gamma = NULL, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)Scale(data, gamma = NULL, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)
data |
Vector of |
gamma |
Vector of |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 scale estimates are computed based on the following model for the CDF:
, where is the scale parameter:
where are the Hill estimates.
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
A |
Vector of the corresponding scale estimates. |
C |
Vector of the corresponding estimates for |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Schoutens, W., De Spiegeleer, J., Reynkens, T. and Herrmann, K. (2016). "Hunting for Black Swans in the European Banking Sector Using Extreme Value Analysis." In: Jan Kallsen and Antonis Papapantoleon (eds.), Advanced Modelling in Mathematical Finance, Springer International Publishing, Switzerland, pp. 147–166.
data(secura) # Hill estimator H <- Hill(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l")data(secura) # Hill estimator H <- Hill(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l")
Computes the bias-reduced estimator for the scale parameter using the second-order Hill estimator.
Scale.2o(data, gamma, b, beta, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)Scale.2o(data, gamma, b, beta, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)
data |
Vector of |
gamma |
Vector of |
b |
Vector of |
beta |
Vector of |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 scale estimates are computed based on the following model for the CDF:
, where is the scale parameter.
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
A |
Vector of the corresponding scale estimates. |
C |
Vector of the corresponding estimates for |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Schoutens, W., De Spiegeleer, J., Reynkens, T. and Herrmann, K. (2016). "Hunting for Black Swans in the European Banking Sector Using Extreme Value Analysis." In: Jan Kallsen and Antonis Papapantoleon (eds.), Advanced Modelling in Mathematical Finance, Springer International Publishing, Switzerland, pp. 147–166.
data(secura) # Hill estimator H <- Hill(secura$size) # Bias-reduced Hill estimator H2o <- Hill.2oQV(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Bias-reduced scale estimator S2o <- Scale.2o(secura$size, gamma=H2o$gamma, b=H2o$b, beta=H2o$beta, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l") lines(S2o$k,log(S2o$A), lty=2)data(secura) # Hill estimator H <- Hill(secura$size) # Bias-reduced Hill estimator H2o <- Hill.2oQV(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Bias-reduced scale estimator S2o <- Scale.2o(secura$size, gamma=H2o$gamma, b=H2o$b, beta=H2o$beta, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l") lines(S2o$k,log(S2o$A), lty=2)
Computes the bias-reduced estimator for the scale parameter using the EPD estimator (Beirlant et al., 2016).
ScaleEPD(data, gamma, kappa, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)ScaleEPD(data, gamma, kappa, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of scale parameter", ...)
data |
Vector of |
gamma |
Vector of |
kappa |
Vector of |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 scale estimates are computed based on the following model for the CDF:
, where is the scale parameter.
Using the EPD approach we replace by .
See Section 4.2.1 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
A |
Vector of the corresponding scale estimates. |
C |
Vector of the corresponding estimates for |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Schoutens, W., De Spiegeleer, J., Reynkens, T. and Herrmann, K. (2016). "Hunting for Black Swans in the European Banking Sector Using Extreme Value Analysis." In: Jan Kallsen and Antonis Papapantoleon (eds.), Advanced Modelling in Mathematical Finance, Springer International Publishing, Switzerland, pp. 147–166.
data(secura) # Hill estimator H <- Hill(secura$size) # EPD estimator epd <- EPD(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Bias-reduced scale estimator Sepd <- ScaleEPD(secura$size, gamma=epd$gamma, kappa=epd$kappa, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l") lines(Sepd$k,log(Sepd$A), lty=2)data(secura) # Hill estimator H <- Hill(secura$size) # EPD estimator epd <- EPD(secura$size) # Scale estimator S <- Scale(secura$size, gamma=H$gamma, plot=FALSE) # Bias-reduced scale estimator Sepd <- ScaleEPD(secura$size, gamma=epd$gamma, kappa=epd$kappa, plot=FALSE) # Plot logarithm of scale plot(S$k,log(S$A), xlab="k", ylab="log(Scale)", type="l") lines(Sepd$k,log(Sepd$A), lty=2)
Estimator of the scale parameter in the regression case where is constant and the regression modelling is thus placed solely on the scale parameter.
ScaleReg(s, Z, kernel = c("normal", "uniform", "triangular", "epanechnikov", "biweight"), h, plot = TRUE, add = FALSE, main = "Estimates of scale parameter", ...)ScaleReg(s, Z, kernel = c("normal", "uniform", "triangular", "epanechnikov", "biweight"), h, plot = TRUE, add = FALSE, main = "Estimates of scale parameter", ...)
s |
Point to evaluate the scale estimator in. |
Z |
Vector of |
kernel |
The kernel used in the estimator. One of |
h |
The bandwidth used in the kernel function. |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 scale estimator is computed as
with the kernel function and the bandwidth.
Here, it is assumed that we have equidistant covariates .
See Section 4.4.1 in Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
A |
Vector of the corresponding scale estimates. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
ProbReg, QuantReg, scale, Hill
data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)data(norwegianfire) Z <- norwegianfire$size[norwegianfire$year==76] i <- 100 n <- length(Z) # Scale estimator in i/n A <- ScaleReg(i/n, Z, h=0.5, kernel = "epanechnikov")$A # Small exceedance probability q <- 10^6 ProbReg(Z, A, q, plot=TRUE) # Large quantile p <- 10^(-5) QuantReg(Z, A, p, plot=TRUE)
Automobile claims from 1988 to 2001, gathered from several European insurance companies, exceeding 1 200 000 Euro. Note that the data were, among others, corrected for inflation.
data("secura")data("secura")
A data frame with 371 observations on the following 2 variables:
yearYear of claim occurence.
sizeSize of automobile insurance claim (in EUR).
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(secura) # Exponential QQ-plot of Secura data ExpQQ(secura$size) # Pareto QQ-plot of Secura data ParetoQQ(secura$size) # Mean excess plot of Secura data (function of k) MeanExcess(secura$size, k=TRUE) # Mean excess plot of Secura data (function of order statistics) MeanExcess(secura$size, k=FALSE)data(secura) # Exponential QQ-plot of Secura data ExpQQ(secura$size) # Pareto QQ-plot of Secura data ParetoQQ(secura$size) # Mean excess plot of Secura data (function of k) MeanExcess(secura$size, k=TRUE) # Mean excess plot of Secura data (function of order statistics) MeanExcess(secura$size, k=FALSE)
The SOA Group Medical Insurance Large Claims Database records, among others, all the claim amounts exceeding 25,000 USD in the year 1991.
data("soa")data("soa")
A data frame with 75789 observations on the following variable:
sizeClaim size (in USD).
Grazier, K. L. and G'Sell Associates (1997). Group Medical Insurance Large Claims Database Collection and Analysis. SOA Monograph M-HB97-1, Society of Actuaries, Schaumburg.
Society of Actuaries, https://www.soa.org/resources/experience-studies/2000-2004/91-92-group-medical-claims/.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
data(soa) # Histogram of log-claim amount hist(log(soa$size),breaks=seq(10,16,0.2),xlab="log(Claim size)") # Exponential QQ-plot of claim amount ExpQQ(soa$size) # Mean excess plot of claim amount (function of k) MeanExcess(soa$size, k=TRUE) # Mean excess plot of claim amount (function of order statistics) MeanExcess(soa$size, k=FALSE)data(soa) # Histogram of log-claim amount hist(log(soa$size),breaks=seq(10,16,0.2),xlab="log(Claim size)") # Exponential QQ-plot of claim amount ExpQQ(soa$size) # Mean excess plot of claim amount (function of k) MeanExcess(soa$size, k=TRUE) # Mean excess plot of claim amount (function of order statistics) MeanExcess(soa$size, k=FALSE)
Density, distribution function, quantile function and random generation for the fitted spliced distribution.
dSplice(x, splicefit, log = FALSE) pSplice(x, splicefit, lower.tail = TRUE, log.p = FALSE) qSplice(p, splicefit, lower.tail = TRUE, log.p = FALSE) rSplice(n, splicefit)dSplice(x, splicefit, log = FALSE) pSplice(x, splicefit, lower.tail = TRUE, log.p = FALSE) qSplice(p, splicefit, lower.tail = TRUE, log.p = FALSE) rSplice(n, splicefit)
x |
Vector of points to evaluate the CDF or PDF in. |
p |
Vector of probabilities. |
n |
Number of observations. |
splicefit |
A |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
See Reynkens et al. (2017) and Section 4.3 in Albrecher et al. (2017) for details.
dSplice gives the density function evaluated in , pSplice the CDF evaluated in and qSplice the quantile function evaluated in . The length of the result is equal to the length of or .
rSplice returns a random sample of length .
Tom Reynkens with R code from Roel Verbelen for the mixed Erlang PDF, CDF and quantiles.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758.
VaR, SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD,
SpliceECDF, SpliceLL, SplicePP
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") p <- seq(0, 1, 0.01) # Plot of splicing quantiles plot(p, qSplice(p, splicefit), type="l", xlab="p", ylab="Q(p)") # Plot of VaR plot(p, VaR(p, splicefit), type="l", xlab="p", ylab=bquote(VaR[p])) # Random sample from spliced distribution x <- rSplice(1000, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") p <- seq(0, 1, 0.01) # Plot of splicing quantiles plot(p, qSplice(p, splicefit), type="l", xlab="p", ylab="Q(p)") # Plot of VaR plot(p, VaR(p, splicefit), type="l", xlab="p", ylab=bquote(VaR[p])) # Random sample from spliced distribution x <- rSplice(1000, splicefit) ## End(Not run)
This function plots the fitted survival function of the spliced distribution together with the
empirical survival function (determined using the Empirical CDF (ECDF)). Moreover, confidence bands are added.
SpliceECDF(x, X, splicefit, alpha = 0.05, ...)SpliceECDF(x, X, splicefit, alpha = 0.05, ...)
x |
Vector of points to plot the functions at. |
X |
Data used for fitting the distribution. |
splicefit |
A |
alpha |
|
... |
Additional arguments for the |
Use SpliceTB for censored data.
Confidence bands are determined using the Dvoretzky-Kiefer-Wolfowitz inequality (Massart, 1990).
See Reynkens et al. (2017) and Section 4.3.1 in Albrecher et al. (2017) for more details.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Massart, P. (1990). The Tight Constant in the Dvoretzky-Kiefer-Wolfowitz Inequality. Annals of Probability, 18, 1269–1283.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758.
SpliceTB, pSplice, ecdf, SpliceFitPareto, SpliceFitGPD, SpliceLL, SplicePP, SpliceQQ
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
Create an S3 object using ME-Pa or ME-GPD splicing fit obtained from SpliceFitPareto, SpliceFiticPareto or
SpliceFitGPD.
SpliceFit(const, trunclower, t, type, MEfit, EVTfit, loglik = NULL, IC = NULL)SpliceFit(const, trunclower, t, type, MEfit, EVTfit, loglik = NULL, IC = NULL)
const |
Vector of splicing constants or a single splicing constant. |
trunclower |
Lower truncation point. |
t |
Vector of splicing points or a single splicing point. |
type |
Vector of types of the distributions: |
MEfit |
|
EVTfit |
|
loglik |
Log-likelihood of the fitted model. When |
IC |
Information criteria of the fitted model. When |
See Reynkens et al. (2017) and Section 4.3 in Albrecher et al. (2017) for details.
An S3 object containing the above input arguments and values for , the splicing weights.
These splicing weights are equal to
when and
when , where is the length of const.
A summary method is available.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
MEfit, EVTfit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)# Create MEfit object mefit <- MEfit(p=c(0.65,0.35), shape=c(39,58), theta=16.19, M=2) # Create EVTfit object evtfit <- EVTfit(gamma=c(0.76,0.64), endpoint=c(39096, Inf)) # Create SpliceFit object splicefit <- SpliceFit(const=c(0.5,0.996), trunclower=0, t=c(1020,39096), type=c("ME","TPa","Pa"), MEfit=mefit, EVTfit=evtfit) # Show summary summary(splicefit)
Fit spliced distribution of a mixed Erlang distribution and a Generalised Pareto Distribution (GPD). The parameters of the GPD are determined using the POT-MLE approach.
SpliceFitGPD(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf)SpliceFitGPD(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf)
X |
Data used for fitting the distribution. |
const |
The probability of the quantile where the ME distribution will be spliced with the GPD distribution. Default is |
tsplice |
The point where the ME distribution will be spliced with the GPD distribution. Default is |
M |
Initial number of Erlang mixtures, default is 3. This number can change when determining an optimal mixed Erlang fit using an information criterion. |
s |
Vector of spread factors for the EM algorithm, default is |
trunclower |
Lower truncation point. Default is 0. |
ncores |
Number of cores to use when determining an optimal mixed Erlang fit using an information criterion.
When |
criterium |
Information criterion used to select the number of components of the ME fit and |
reduceM |
Logical indicating if M should be reduced based on the information criterion, default is |
eps |
Covergence threshold used in the EM algorithm (ME part). Default is |
beta_tol |
Threshold for the mixing weights below which the corresponding shape parameter vector is considered neglectable (ME part). Default is |
maxiter |
Maximum number of iterations in a single EM algorithm execution (ME part). Default is |
See Reynkens et al. (2017), Section 4.3.1 of Albrecher et al. (2017) and Verbelen et al. (2015) for details. The code follows the notation of the latter. Initial values follow from Verbelen et al. (2016).
A SpliceFit object.
Tom Reynkens with R code from Roel Verbelen for fitting the mixed Erlang distribution.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758.
Verbelen, R., Antonio, K. and Claeskens, G. (2016). "Multivariate Mixtures of Erlangs for Density Estimation Under Censoring." Lifetime Data Analysis, 22, 429–455.
SpliceFitPareto, SpliceFiticPareto, Splice,
GPDfit
## Not run: # GPD random sample X <- rgpd(1000, gamma = 0.5, sigma = 2) # Splice ME and GPD splicefit <- SpliceFitGPD(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # GPD random sample X <- rgpd(1000, gamma = 0.5, sigma = 2) # Splice ME and GPD splicefit <- SpliceFitGPD(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
Fit spliced distribution of a mixed Erlang distribution and a Pareto distribution adapted for interval censoring and truncation.
SpliceFiticPareto(L, U, censored, tsplice, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, ncores = NULL, criterium = c("BIC", "AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, cpp = FALSE)SpliceFiticPareto(L, U, censored, tsplice, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, ncores = NULL, criterium = c("BIC", "AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf, cpp = FALSE)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
tsplice |
The splicing point |
M |
Initial number of Erlang mixtures, default is 3. This number can change when determining an optimal mixed Erlang fit using an information criterion. |
s |
Vector of spread factors for the EM algorithm, default is |
trunclower |
Lower truncation point. Default is 0. |
truncupper |
Upper truncation point. Default is |
ncores |
Number of cores to use when determining an optimal mixed Erlang fit using an information criterion.
When |
criterium |
Information criterion used to select the number of components of the ME fit and |
reduceM |
Logical indicating if M should be reduced based on the information criterion, default is |
eps |
Covergence threshold used in the EM algorithm. Default is |
beta_tol |
Threshold for the mixing weights below which the corresponding shape parameter vector is considered neglectable (ME part). Default is |
maxiter |
Maximum number of iterations in a single EM algorithm execution. Default is |
cpp |
Use |
See Reynkens et al. (2017), Section 4.3.2 of Albrecher et al. (2017) and Verbelen et al. (2015) for details. The code follows the notation of the latter. Initial values follow from Verbelen et al. (2016).
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u.
A SpliceFit object.
Tom Reynkens based on R code from Roel Verbelen for fitting the mixed Erlang distribution (without splicing).
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758.
Verbelen, R., Antonio, K. and Claeskens, G. (2016). "Multivariate Mixtures of Erlangs for Density Estimation Under Censoring." Lifetime Data Analysis, 22, 429–455.
SpliceFitPareto, SpliceFitGPD, Splice
## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)
Fit spliced distribution of a mixed Erlang distribution and Pareto distribution(s). The shape parameter(s) of the Pareto distribution(s) is determined using the Hill estimator.
SpliceFitPareto(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, EVTtruncation = FALSE, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf) SpliceFitHill(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, EVTtruncation = FALSE, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf)SpliceFitPareto(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, EVTtruncation = FALSE, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf) SpliceFitHill(X, const = NULL, tsplice = NULL, M = 3, s = 1:10, trunclower = 0, truncupper = Inf, EVTtruncation = FALSE, ncores = NULL, criterium = c("BIC","AIC"), reduceM = TRUE, eps = 10^(-3), beta_tol = 10^(-5), maxiter = Inf)
X |
Data used for fitting the distribution. |
const |
Vector of length |
tsplice |
Vector of length |
M |
Initial number of Erlang mixtures, default is 3. This number can change when determining an optimal mixed Erlang fit using an information criterion. |
s |
Vector of spread factors for the EM algorithm, default is |
trunclower |
Lower truncation point. Default is 0. |
truncupper |
Upper truncation point. Default is |
EVTtruncation |
Logical indicating if the |
ncores |
Number of cores to use when determining an optimal mixed Erlang fit using an information criterion.
When |
criterium |
Information criterion used to select the number of components of the ME fit and |
reduceM |
Logical indicating if M should be reduced based on the information criterion, default is |
eps |
Covergence threshold used in the EM algorithm (ME part). Default is |
beta_tol |
Threshold for the mixing weights below which the corresponding shape parameter vector is considered neglectable (ME part). Default is |
maxiter |
Maximum number of iterations in a single EM algorithm execution (ME part). Default is |
See Reynkens et al. (2017), Section 4.3.1 of Albrecher et al. (2017) and Verbelen et al. (2015) for details. The code follows the notation of the latter. Initial values follow from Verbelen et al. (2016).
The SpliceFitHill function is the same function but with a different name for compatibility with old versions of the package.
Use SpliceFiticPareto when censoring is present.
A SpliceFit object.
Tom Reynkens with R code from Roel Verbelen for fitting the mixed Erlang distribution.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758.
Verbelen, R., Antonio, K. and Claeskens, G. (2016). "Multivariate Mixtures of Erlangs for Density Estimation Under Censoring." Lifetime Data Analysis, 22, 429–455.
SpliceFiticPareto, SpliceFitGPD, Splice,
Hill, trHill
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
This function plots the logarithm of the empirical survival function (determined using the Empirical CDF (ECDF)) versus the logarithm of the data. Moreover, the logarithm of the fitted survival function of the spliced distribution is added.
SpliceLL(x = sort(X), X, splicefit, plot = TRUE, main = "Splicing LL-plot", ...)SpliceLL(x = sort(X), X, splicefit, plot = TRUE, main = "Splicing LL-plot", ...)
x |
Vector of points to plot the fitted survival function at. By default we plot it at the data points. |
X |
Data used for fitting the distribution. |
splicefit |
A |
plot |
Logical indicating if the splicing LL-plot should be made, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The LL-plot consists of the points
for with the length of the data, the -th smallest observation
and the empirical distribution function.
Then, the line
with the fitted spliced distribution function, is added.
Use SpliceLL_TB for censored data.
See Reynkens et al. (2017) and Section 4.3.1 in Albrecher et al. (2017) for more details.
A list with following components:
logX |
Vector of the logarithms of the sorted data. |
sll.the |
Vector of the theoretical log-probabilities |
logx |
Vector of the logarithms of the points to plot the fitted survival function at. |
sll.emp |
Vector of the empirical log-probabilities |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceLL_TB, pSplice, ecdf, SpliceFitPareto, SpliceFitGPD, SpliceECDF, SplicePP, SpliceQQ
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
This function plots the logarithm of the Turnbull survival function (which is suitable for interval censored data) versus the logarithm of the data. Moreover, the logarithm of the fitted survival function of the spliced distribution is added.
SpliceLL_TB(x = sort(L), L, U = L, censored, splicefit, plot = TRUE, main = "Splicing LL-plot", ...)SpliceLL_TB(x = sort(L), L, U = L, censored, splicefit, plot = TRUE, main = "Splicing LL-plot", ...)
x |
Vector of points to plot the fitted survival function at. By default we plot it at the points |
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
splicefit |
A |
plot |
Logical indicating if the splicing LL-plot should be made, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The LL-plot consists of the points
for with the length of the data, the -th smallest observation
and the Turnbull estimator for the distribution function.
Then, the line
with the fitted spliced distribution function, is added.
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u. The limits trunclower and truncupper are obtained from the SpliceFit object.
If the interval package is installed, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Use SpliceLL for non-censored data.
See Reynkens et al. (2017) and Section 4.3.2 in Albrecher et al. (2017) for more details.
A list with following components:
logX |
Vector of the logarithms of the sorted left boundaries of the intervals. |
sll.the |
Vector of the theoretical log-probabilities |
logx |
Vector of the logarithms of the points to plot the fitted survival function at. |
sll.emp |
Vector of the empirical log-probabilities |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceLL, pSplice, Turnbull, icfit, SpliceFiticPareto, SpliceTB, SplicePP_TB, SpliceQQ_TB
## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)
This function plots the fitted survival function of the spliced distribution versus the empirical survival function (determined using the Empirical CDF (ECDF)).
SplicePP(X, splicefit, x = sort(X), log = FALSE, plot = TRUE, main = "Splicing PP-plot", ...)SplicePP(X, splicefit, x = sort(X), log = FALSE, plot = TRUE, main = "Splicing PP-plot", ...)
X |
Data used for fitting the distribution. |
splicefit |
A |
x |
Vector of points to plot the functions at. By default we plot them at the data points. |
log |
Logical indicating if minus the logarithms of the survival probabilities are plotted versus each other, default is |
plot |
Logical indicating if the splicing PP-plot should be made, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The PP-plot consists of the points
for with the length of the data, the -th smallest observation,
the empirical distribution function and the fitted spliced distribution function.
The minus-log version of the PP-plot consists of
Use SplicePP_TB for censored data.
See Reynkens et al. (2017) and Section 4.3.1 in Albrecher et al. (2017) for more details.
A list with following components:
spp.the |
Vector of the theoretical probabilities |
spp.emp |
Vector of the empirical probabilities |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SplicePP_TB, pSplice, ecdf, SpliceFitPareto, SpliceFitGPD, SpliceECDF, SpliceLL, SpliceQQ
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
This function plots the fitted survival function of the spliced distribution versus the Turnbull survival function (which is suitable for interval censored data).
SplicePP_TB(L, U = L, censored, splicefit, x = NULL, log = FALSE, plot = TRUE, main = "Splicing PP-plot", ...)SplicePP_TB(L, U = L, censored, splicefit, x = NULL, log = FALSE, plot = TRUE, main = "Splicing PP-plot", ...)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
splicefit |
A |
x |
Vector of points to plot the functions at. When |
log |
Logical indicating if minus the logarithms of the survival probabilities are plotted versus each other, default is |
plot |
Logical indicating if the splicing PP-plot should be made, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The PP-plot consists of the points
for with the length of the data, where , is the quantile function obtained using the Turnbull estimator,
the Turnbull estimator for the distribution function and the fitted spliced distribution function.
The minus-log version of the PP-plot consists of
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u. The limits trunclower and truncupper are obtained from the SpliceFit object.
If the interval package is installed, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Use SplicePP for non-censored data.
See Reynkens et al. (2017) and Section 4.3.2 in Albrecher et al. (2017) for more details.
A list with following components:
spp.the |
Vector of the theoretical probabilities |
spp.emp |
Vector of the empirical probabilities |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SplicePP, pSplice, Turnbull, icfit, SpliceFiticPareto, SpliceTB, SpliceLL_TB, SpliceQQ_TB
## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)
Computes the empirical quantiles of a data vector and the theoretical quantiles of the fitted spliced distribution. These quantiles are then plotted in a splicing QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
SpliceQQ(X, splicefit, p = NULL, plot = TRUE, main = "Splicing QQ-plot", ...)SpliceQQ(X, splicefit, p = NULL, plot = TRUE, main = "Splicing QQ-plot", ...)
X |
Vector of |
splicefit |
A |
p |
Vector of probabilities used in the QQ-plot. If |
plot |
Logical indicating if the quantiles should be plotted in a splicing QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
This QQ-plot is given by
for where is the quantile function of the fitted splicing model and is the empirical quantile function and .
See Reynkens et al. (2017) and Section 4.3.1 in Albrecher et al. (2017) for more details.
A list with following components:
sqq.the |
Vector of the theoretical quantiles of the fitted spliced distribution. |
sqq.emp |
Vector of the empirical quantiles from the data. |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceQQ_TB, qSplice, SpliceFitPareto, SpliceFitGPD, SpliceECDF, SpliceLL, SplicePP
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) x <- seq(0, 20, 0.01) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and empirical survival function SpliceECDF(x, X, splicefit) # Log-log plot with empirical survival function and fitted survival function SpliceLL(x, X, splicefit) # PP-plot of empirical survival function and fitted survival function SplicePP(X, splicefit) # PP-plot of empirical survival function and # fitted survival function with log-scales SplicePP(X, splicefit, log=TRUE) # Splicing QQ-plot SpliceQQ(X, splicefit) ## End(Not run)
This function plots the fitted quantile function of the spliced distribution versus quantiles based on the Turnbull survival function (which is suitable for interval censored data).
SpliceQQ_TB(L, U = L, censored, splicefit, p = NULL, plot = TRUE, main = "Splicing QQ-plot", ...)SpliceQQ_TB(L, U = L, censored, splicefit, p = NULL, plot = TRUE, main = "Splicing QQ-plot", ...)
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
splicefit |
A |
p |
Vector of probabilities used in the QQ-plot. If |
plot |
Logical indicating if the quantiles should be plotted in a splicing QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
This QQ-plot is given by
for where is the quantile function of the fitted splicing model, the quantile function obtained using the Turnbull estimator and .
If the interval package is installed, the icfit function is used to compute the Turnbull estimator. Otherwise, survfit.formula from survival is used.
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u. The limits trunclower and truncupper are obtained from the SpliceFit object.
Use SpliceQQ for non-censored data.
See Reynkens et al. (2017) and Section 4.3.2 in Albrecher et al. (2017) for more details.
A list with following components:
sqq.the |
Vector of the theoretical quantiles of the fitted spliced distribution. |
sqq.emp |
Vector of the empirical quantiles from the data (based on the Turnbull estimator). |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceQQ, qSplice, Turnbull, icfit, SpliceFiticPareto, SpliceTB, SplicePP_TB, SpliceLL_TB
## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)
This function plots the fitted survival function of the spliced distribution together with the
Turnbull survival function (which is suitable for interval censored data). Moreover, confidence intervals are added.
SpliceTB(x = sort(L), L, U = L, censored, splicefit, alpha = 0.05, ...)SpliceTB(x = sort(L), L, U = L, censored, splicefit, alpha = 0.05, ...)
x |
Vector of points to plot the functions at. By default we plot it at the points |
L |
Vector of length |
U |
Vector of length |
censored |
A logical vector of length |
splicefit |
A |
alpha |
|
... |
Additional arguments for the |
Right censored data should be entered as L=l and U=truncupper, and left censored data should be entered as L=trunclower and U=u. The limits trunclower and truncupper are obtained from the SpliceFit object.
Use SpliceECDF for non-censored data.
See Reynkens et al. (2017) and Section 4.3.2 in Albrecher et al. (2017) for more details.
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
SpliceECDF, pSplice, Turnbull, SpliceFiticPareto, SpliceLL_TB, SplicePP_TB, SpliceQQ_TB
## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)## Not run: # 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) # Right boundary U <- Z U[censored] <- Inf # Splice ME and Pareto splicefit <- SpliceFiticPareto(L=Z, U=U, censored=censored, tsplice=quantile(Z,0.9)) x <- seq(0,20,0.1) # Plot of spliced CDF plot(x, pSplice(x, splicefit), type="l", xlab="x", ylab="F(x)") # Plot of spliced PDF plot(x, dSplice(x, splicefit), type="l", xlab="x", ylab="f(x)") # Fitted survival function and Turnbull survival function SpliceTB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # Log-log plot with Turnbull survival function and fitted survival function SpliceLL_TB(x, L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and fitted survival function SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit) # PP-plot of Turnbull survival function and # fitted survival function with log-scales SplicePP_TB(L=Z, U=U, censored=censored, splicefit=splicefit, log=TRUE) # QQ-plot using Turnbull survival function and fitted survival function SpliceQQ_TB(L=Z, U=U, censored=censored, splicefit=splicefit) ## End(Not run)
Non-parametric estimators of the stable tail dependence function (STDF): and .
stdf(x, k, X, alpha = 0.5) stdf2(x, k, X)stdf(x, k, X, alpha = 0.5) stdf2(x, k, X)
x |
A |
k |
Value of the tail index |
X |
A data matrix of dimensions |
alpha |
The parameter |
The stable tail dependence function in can be estimated by
with
where is the rank of among the observations in the -th dimension:
This estimator is implemented in stdf.
The second estimator is given by
where is the -th smallest observation in the -th dimension.
This estimator is implemented in stdf2.
See Section 4.5 of Beirlant et al. (2016) for more details.
stdf returns the estimate and stdf2 returns the estimate .
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant J., Goegebeur Y., Segers, J. and Teugels, J. (2004). Statistics of Extremes: Theory and Applications, Wiley Series in Probability, Wiley, Chichester.
# Generate data matrix X <- cbind(rpareto(100,2), rpareto(100,3)) # Tail index k <- 20 # Point to evaluate the STDF in x <- c(2,3) # First estimate stdf(x, k, X) # Second estimate stdf2(x, k, X)# Generate data matrix X <- cbind(rpareto(100,2), rpareto(100,3)) # Tail index k <- 20 # Point to evaluate the STDF in x <- c(2,3) # First estimate stdf(x, k, X) # Second estimate stdf2(x, k, X)
Density, distribution function, quantile function and random generation for the truncated Burr distribution (type XII).
dtburr(x, alpha, rho, eta = 1, endpoint = Inf, log = FALSE) ptburr(x, alpha, rho, eta = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtburr(p, alpha, rho, eta = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtburr(n, alpha, rho, eta = 1, endpoint = Inf)dtburr(x, alpha, rho, eta = 1, endpoint = Inf, log = FALSE) ptburr(x, alpha, rho, eta = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtburr(p, alpha, rho, eta = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtburr(n, alpha, rho, eta = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
alpha |
The |
rho |
The |
eta |
The |
endpoint |
Endpoint of the truncated Burr distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated Burr distribution is equal to
for where is the CDF of the ordinary Burr distribution and is the endpoint (truncation point) of the truncated Burr distribution.
dtburr gives the density function evaluated in , ptburr the CDF evaluated in and qtburr the quantile function evaluated in . The length of the result is equal to the length of or .
rtburr returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtburr(x, alpha=2, rho=-1, endpoint=9), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptburr(x, alpha=2, rho=-1, endpoint=9), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtburr(x, alpha=2, rho=-1, endpoint=9), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptburr(x, alpha=2, rho=-1, endpoint=9), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the truncated exponential distribution.
dtexp(x, rate = 1, endpoint = Inf, log = FALSE) ptexp(x, rate = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtexp(p, rate = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtexp(n, rate = 1, endpoint = Inf)dtexp(x, rate = 1, endpoint = Inf, log = FALSE) ptexp(x, rate = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtexp(p, rate = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtexp(n, rate = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
rate |
The rate parameter for the exponential distribution, default is 1. |
endpoint |
Endpoint of the truncated exponential distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated exponential distribution is equal to
for where is the CDF of the ordinary exponential distribution and is the endpoint (truncation point) of the truncated exponential distribution.
dtexp gives the density function evaluated in , ptexp the CDF evaluated in and qtexp the quantile function evaluated in . The length of the result is equal to the length of or .
rtexp returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtexp(x, rate = 2, endpoint=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptexp(x, rate = 2, endpoint=5), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtexp(x, rate = 2, endpoint=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptexp(x, rate = 2, endpoint=5), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the truncated Fréchet distribution.
dtfrechet(x, shape, loc = 0, scale = 1, endpoint = Inf, log = FALSE) ptfrechet(x, shape, loc = 0, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtfrechet(p, shape, loc = 0, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtfrechet(n, shape, loc = 0, scale = 1, endpoint = Inf)dtfrechet(x, shape, loc = 0, scale = 1, endpoint = Inf, log = FALSE) ptfrechet(x, shape, loc = 0, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtfrechet(p, shape, loc = 0, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtfrechet(n, shape, loc = 0, scale = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
shape |
Shape parameter of the Fréchet distribution. |
loc |
Location parameter of the Fréchet distribution, default is 0. |
scale |
Scale parameter of the Fréchet distribution, default is 1. |
endpoint |
Endpoint of the truncated Fréchet distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated Fréchet distribution is equal to
for where is the CDF of an ordinary Fréchet distribution and is the endpoint (truncation point) of the truncated Fréchet distribution.
dtfrechet gives the density function evaluated in , ptfrechet the CDF evaluated in and qtfrechet the quantile function evaluated in . The length of the result is equal to the length of or .
rtfrechet returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(1, 10, 0.01) plot(x, dtfrechet(x, shape=2, endpoint=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1, 10, 0.01) plot(x, ptfrechet(x, shape=2, endpoint=5), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(1, 10, 0.01) plot(x, dtfrechet(x, shape=2, endpoint=5), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(1, 10, 0.01) plot(x, ptfrechet(x, shape=2, endpoint=5), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the truncated Generalised Pareto Distribution (GPD).
dtgpd(x, gamma, mu = 0, sigma, endpoint = Inf, log = FALSE) ptgpd(x, gamma, mu = 0, sigma, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtgpd(p, gamma, mu = 0, sigma, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtgpd(n, gamma, mu = 0, sigma, endpoint = Inf)dtgpd(x, gamma, mu = 0, sigma, endpoint = Inf, log = FALSE) ptgpd(x, gamma, mu = 0, sigma, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtgpd(p, gamma, mu = 0, sigma, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtgpd(n, gamma, mu = 0, sigma, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
gamma |
The |
mu |
The |
sigma |
The |
endpoint |
Endpoint of the truncated GPD. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated GPD is equal to
for where is the CDF of the ordinary GPD and is the endpoint (truncation point) of the truncated GPD.
dtgpd gives the density function evaluated in , ptgpd the CDF evaluated in and qtgpd the quantile function evaluated in . The length of the result is equal to the length of or .
rtgpd returns a random sample of length .
Tom Reynkens
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtgpd(x, gamma=1/2, sigma=5, endpoint=8), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptgpd(x, gamma=1/2, sigma=5, endpoint=8), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtgpd(x, gamma=1/2, sigma=5, endpoint=8), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptgpd(x, gamma=1/2, sigma=5, endpoint=8), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the truncated log-normal distribution.
dtlnorm(x, meanlog = 0, sdlog = 1, endpoint = Inf, log = FALSE) ptlnorm(x, meanlog = 0, sdlog = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtlnorm(p, meanlog = 0, sdlog = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtlnorm(n, meanlog = 0, sdlog = 1, endpoint = Inf)dtlnorm(x, meanlog = 0, sdlog = 1, endpoint = Inf, log = FALSE) ptlnorm(x, meanlog = 0, sdlog = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtlnorm(p, meanlog = 0, sdlog = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtlnorm(n, meanlog = 0, sdlog = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
meanlog |
Mean of the distribution on the log scale, default is 0. |
sdlog |
Standard deviation of the distribution on the log scale, default is 1. |
endpoint |
Endpoint of the truncated log-normal distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated log-normal distribution is equal to
for where is the CDF of the ordinary log-normal distribution and is the endpoint (truncation point) of the truncated log-normal distribution.
dtlnorm gives the density function evaluated in , ptlnorm the CDF evaluated in and qtlnorm the quantile function evaluated in . The length of the result is equal to the length of or .
rtlnorm returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtlnorm(x, endpoint=9), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptlnorm(x, endpoint=9), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtlnorm(x, endpoint=9), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptlnorm(x, endpoint=9), xlab="x", ylab="CDF", type="l")
Density, distribution function, quantile function and random generation for the truncated Pareto distribution.
dtpareto(x, shape, scale = 1, endpoint = Inf, log = FALSE) ptpareto(x, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtpareto(p, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtpareto(n, shape, scale = 1, endpoint = Inf)dtpareto(x, shape, scale = 1, endpoint = Inf, log = FALSE) ptpareto(x, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtpareto(p, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtpareto(n, shape, scale = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
shape |
The shape parameter of the truncated Pareto distribution, a strictly positive number. |
scale |
The scale parameter of the truncated Pareto distribution, a strictly positive number. Its default value is |
endpoint |
Endpoint of the truncated Pareto distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated Pareto distribution is equal to
for where is the CDF of an ordinary Pareto distribution and is the endpoint (truncation point) of the truncated Pareto distribution.
dtpareto gives the density function evaluated in , ptpareto the CDF evaluated in and qtpareto the quantile function evaluated in . The length of the result is equal to the length of or .
rtpareto returns a random sample of length .
Tom Reynkens
# Plot of the PDF x = seq(1,10,0.01) plot(x, dtpareto(x, shape=2, endpoint=10), xlab="x", ylab="PDF", type="l") # Plot of the CDF x = seq(1,10,0.01) plot(x, ptpareto(x, shape=2, endpoint=10), xlab="x", ylab="CDF", type="l")# Plot of the PDF x = seq(1,10,0.01) plot(x, dtpareto(x, shape=2, endpoint=10), xlab="x", ylab="PDF", type="l") # Plot of the CDF x = seq(1,10,0.01) plot(x, ptpareto(x, shape=2, endpoint=10), xlab="x", ylab="CDF", type="l")
Estimates of truncation odds of the truncated probability mass under the untruncated distribution using truncated Hill.
trDT(data, r = 1, gamma, plot = FALSE, add = FALSE, main = "Estimates of DT", ...)trDT(data, r = 1, gamma, plot = FALSE, add = FALSE, main = "Estimates of DT", ...)
data |
Vector of |
r |
Trimming parameter, default is |
gamma |
Vector of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The truncation odds is defined as
with the upper truncation point and the CDF of the untruncated distribution (e.g. Pareto distribution).
We estimate this truncation odds as
with .
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
DT |
Vector of the corresponding estimates for the truncation odds |
Tom Reynkens based on R code of Dries Cornilly.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
trHill, trEndpoint, trQuant, trDTMLE
# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dt <- trDT(X, gamma=trh$gamma, plot=TRUE, ylim=c(0,0.05))# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dt <- trDT(X, gamma=trh$gamma, plot=TRUE, ylim=c(0,0.05))
Estimates of truncation odds of the truncated probability mass under the untruncated distribution using truncated MLE.
trDTMLE(data, gamma, tau, plot = FALSE, add = FALSE, main = "Estimates of DT", ...)trDTMLE(data, gamma, tau, plot = FALSE, add = FALSE, main = "Estimates of DT", ...)
data |
Vector of |
gamma |
Vector of |
tau |
Vector of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The truncation odds is defined as
with the upper truncation point and the CDF of the untruncated distribution (e.g. GPD).
We estimate this truncation odds as
with .
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
DT |
Vector of the corresponding estimates for the truncation odds |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trMLE, trEndpointMLE, trProbMLE, trQuantMLE, trTestMLE, trDT
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=TRUE, ylim=c(0,0.05))# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=TRUE, ylim=c(0,0.05))
Estimator of endpoint using truncated Hill estimates.
trEndpoint(data, r = 1, gamma, plot = FALSE, add = FALSE, main = "Estimates of endpoint", ...)trEndpoint(data, r = 1, gamma, plot = FALSE, add = FALSE, main = "Estimates of endpoint", ...)
data |
Vector of |
r |
Trimming parameter, default is |
gamma |
Vector of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The endpoint is estimated as
with the Hill estimates adapted for truncation.
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Tk |
Vector of the corresponding estimates for the endpoint |
Tom Reynkens based on R code of Dries Cornilly.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Endpoint trEndpoint(X, gamma=trh$gamma, plot=TRUE, ylim=c(8,12)) abline(h=qpareto(0.99, shape=shape), lty=2)# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Endpoint trEndpoint(X, gamma=trh$gamma, plot=TRUE, ylim=c(8,12)) abline(h=qpareto(0.99, shape=shape), lty=2)
Estimator of endpoint using truncated ML estimates.
trEndpointMLE(data, gamma, tau, plot = FALSE, add = FALSE, main = "Estimates of endpoint", ...)trEndpointMLE(data, gamma, tau, plot = FALSE, add = FALSE, main = "Estimates of endpoint", ...)
data |
Vector of |
gamma |
Vector of |
tau |
Vector of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The endpoint is estimated as
with and the truncated ML estimates for and .
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Tk |
Vector of the corresponding estimates for the endpoint |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trMLE, trDTMLE, trProbMLE, trQuantMLE, trTestMLE, trEndpoint
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Endpoint trEndpointMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=TRUE, ylim=c(0,50)) abline(h=qgpd(0.99, gamma=gamma, sigma=sigma), lty=2)# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Endpoint trEndpointMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=TRUE, ylim=c(0,50)) abline(h=qgpd(0.99, gamma=gamma, sigma=sigma), lty=2)
Computes the Hill estimator for positive extreme value indices, adapted for upper truncation, as a function of the tail parameter (Aban et al. 2006; Beirlant et al., 2016). Optionally, these estimates are plotted as a function of .
trHill(data, r = 1, tol = 1e-08, maxiter = 100, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of the EVI", ...)trHill(data, r = 1, tol = 1e-08, maxiter = 100, logk = FALSE, plot = FALSE, add = FALSE, main = "Estimates of the EVI", ...)
data |
Vector of |
r |
Trimming parameter, default is |
tol |
Numerical tolerance for stopping criterion used in Newton-Raphson iterations, default is |
maxiter |
Maximum number of Newton-Raphson iterations, default is |
logk |
Logical indicating if the estimates are plotted as a function of |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The truncated Hill estimator is the MLE for under the truncated Pareto distribution.
To estimate the EVI using the truncated Hill estimator an equation needs to be solved. Beirlant et al. (2016) propose to use Newton-Raphson iterations to solve this equation. We take the trimmed Hill estimates as starting values for this algorithm. The trimmed Hill estimator is defined as
for
and is a basic extension of the Hill estimator for upper truncated data (the ordinary Hill estimator is obtained for ).
The equation that needs to be solved is
with .
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding estimates for |
H |
Vector of corresponding trimmed Hill estimates. |
Tom Reynkens based on R code of Dries Cornilly.
Aban, I.B., Meerschaert, M.M. and Panorska, A.K. (2006). "Parameter Estimation for the Truncated Pareto Distribution." Journal of the American Statistical Association, 101, 270–277.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
Hill, trDT, trEndpoint, trProb, trQuant, trMLE
# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2))# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2))
Computes the ML estimator for the extreme value index, adapted for upper truncation, as a function of the tail parameter (Beirlant et al., 2017). Optionally, these estimates are plotted as a function of .
trMLE(data, start = c(1, 1), eps = 10^(-10), plot = TRUE, add = FALSE, main = "Estimates for EVI", ...)trMLE(data, start = c(1, 1), eps = 10^(-10), plot = TRUE, add = FALSE, main = "Estimates for EVI", ...)
data |
Vector of |
start |
Starting values for |
eps |
Numerical tolerance, see Details. By default it is equal to |
plot |
Logical indicating if the estimates of |
add |
Logical indicating if the estimates of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
We compute the MLE for the and parameters of the truncated GPD.
For numerical reasons, we compute the MLE for and transform this estimate to .
The log-likelihood is given by
with .
In order to meet the restrictions and for , we require the estimates of these quantities to be larger than the numerical tolerance value eps.
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
gamma |
Vector of the corresponding estimates for |
tau |
Vector of the corresponding estimates for |
sigma |
Vector of the corresponding estimates for |
conv |
Convergence indicator of |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trDTMLE, trEndpointMLE, trProbMLE, trQuantMLE, trTestMLE, trHill, GPDmle
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2))# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2))
Extension of the Pareto QQ-plot as described in Beirlant et al. (2016).
trParetoQQ(data, r = 1, DT, kstar = NULL, plot = TRUE, main = "TPa QQ-plot", ...)trParetoQQ(data, r = 1, DT, kstar = NULL, plot = TRUE, main = "TPa QQ-plot", ...)
data |
Vector of |
r |
Trimming parameter, default is |
DT |
Vector of |
kstar |
Value for |
plot |
Logical indicating if the quantiles should be plotted in a Pareto QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Pareto QQ-plot for truncated data plots
for .
The value for can be be given by the user or can be determined automatically.
In the latter case, we use the that maximises the absolute value of the correlation between and for and .
When taking , one obtains the ordinary Pareto QQ-plot.
Note that the definition here differs slightly from the one in Beirlant et al. (2016). We plot the empirical and theoretical quantiles the other way around and therefore have to add a minus (before the log).
See Beirlant et al. (2016) for more details.
A list with following components:
pqq.the |
Vector of theoretical quantiles |
pqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
kstar |
Optimal value for |
DTstar |
Estimate of |
Tom Reynkens.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
# Endpoint of truncated Pareto distribution endpoint <- qpareto(0.99, shape=2) # Generate sample from truncated Pareto distribution X <- rtpareto(1000, shape=2, endpoint=endpoint) # Ordinary Pareto QQ-plot ParetoQQ(X) # Truncated Hill estimates gamma <- trHill(X)$gamma # Estimates for truncation odds dt <- trDT(X, gamma=gamma)$DT # Truncated Pareto QQ-plot trParetoQQ(X, DT=dt)# Endpoint of truncated Pareto distribution endpoint <- qpareto(0.99, shape=2) # Generate sample from truncated Pareto distribution X <- rtpareto(1000, shape=2, endpoint=endpoint) # Ordinary Pareto QQ-plot ParetoQQ(X) # Truncated Hill estimates gamma <- trHill(X)$gamma # Estimates for truncation odds dt <- trDT(X, gamma=gamma)$DT # Truncated Pareto QQ-plot trParetoQQ(X, DT=dt)
Computes estimates of a small exceedance probability using the estimates for the EVI obtained from the Hill estimator adapted for upper truncation.
trProb(data, r = 1, gamma, q, warnings = TRUE, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)trProb(data, r = 1, gamma, q, warnings = TRUE, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)
data |
Vector of |
r |
Trimming parameter, default is |
gamma |
Vector of |
q |
The used large quantile (we estimate |
warnings |
Logical indicating if warnings are shown, default is |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with and the Hill estimates adapted for truncation.
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
P |
Vector of the corresponding probability estimates. |
q |
The used large quantile. |
Tom Reynkens based on R code of Dries Cornilly.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
trHill, trQuant, Prob, trProbMLE
# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Small probability trProb(X, gamma=trh$gamma, q=8, plot=TRUE)# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Small probability trProb(X, gamma=trh$gamma, q=8, plot=TRUE)
Computes estimates of a small exceedance probability using the estimates for the EVI obtained from the ML estimator adapted for upper truncation.
trProbMLE(data, gamma, tau, DT, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)trProbMLE(data, gamma, tau, DT, q, plot = FALSE, add = FALSE, main = "Estimates of small exceedance probability", ...)
data |
Vector of |
gamma |
Vector of |
tau |
Vector of |
DT |
Vector of |
q |
The used large quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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
with and the ML estimates adapted for truncation and the estimates for the truncation odds.
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
P |
Vector of the corresponding probability estimates. |
q |
The used large quantile. |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trMLE, trDTMLE, trQuantMLE, trEndpointMLE, trTestMLE, trProb, Prob
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=FALSE) # Small exceedance probability trProbMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, q=26, ylim=c(0,0.005))# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=FALSE) # Small exceedance probability trProbMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, q=26, ylim=c(0,0.005))
trQuant computes estimates of large quantiles of the truncated distribution using the estimates for the EVI obtained from the Hill estimator adapted for upper truncation. trQuantW computes estimates of large quantiles of the parent distribution which is unobserved.
trQuant(data, r = 1, rough = TRUE, gamma, DT, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) trQuantW(data, gamma, DT, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)trQuant(data, r = 1, rough = TRUE, gamma, DT, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...) trQuantW(data, gamma, DT, p, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
r |
Trimming parameter, default is |
rough |
Logical indicating if rough truncation is present, default is |
gamma |
Vector of |
DT |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
We observe the truncated r.v. where is the truncation point and the untruncated r.v.
Under rough truncation, the quantiles for are estimated using
with the Hill estimates adapted for truncation and the estimates for the truncation odds.
Under light truncation, the quantiles are estimated using the Weissman estimator with the Hill estimates replaced by the truncated Hill estimates:
To decide between light and rough truncation, one can use the test implemented in trTest.
The quantiles for are estimated using
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens based on R code of Dries Cornilly.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
trHill, trDT, trProb, trEndpoint, trTest, Quant, trQuantMLE
# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dt <- trDT(X, gamma=trh$gamma, plot=TRUE, ylim=c(0,2)) # Large quantile p <- 10^(-5) # Truncated distribution trQuant(X, gamma=trh$gamma, DT=dt$DT, p=p, plot=TRUE) # Original distribution trQuantW(X, gamma=trh$gamma, DT=dt$DT, p=p, plot=TRUE, ylim=c(0,1000))# Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Truncated Hill estimator trh <- trHill(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dt <- trDT(X, gamma=trh$gamma, plot=TRUE, ylim=c(0,2)) # Large quantile p <- 10^(-5) # Truncated distribution trQuant(X, gamma=trh$gamma, DT=dt$DT, p=p, plot=TRUE) # Original distribution trQuantW(X, gamma=trh$gamma, DT=dt$DT, p=p, plot=TRUE, ylim=c(0,1000))
This function computes estimates of large quantiles of the truncated distribution using the ML estimates adapted for upper truncation. Moreover, estimates of large quantiles of the original distribution , which is unobserved, are also computed.
trQuantMLE(data, gamma, tau, DT, p, Y = FALSE, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)trQuantMLE(data, gamma, tau, DT, p, Y = FALSE, plot = FALSE, add = FALSE, main = "Estimates of extreme quantile", ...)
data |
Vector of |
gamma |
Vector of |
tau |
Vector of |
DT |
Vector of |
p |
The exceedance probability of the quantile (we estimate |
Y |
Logical indicating if quantiles from the truncated distribution ( |
plot |
Logical indicating if the estimates should be plotted as a function of |
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 |
We observe the truncated r.v. where is the truncation point and the untruncated r.v.
Under rough truncation, the quantiles for are estimated using
with and the ML estimates adapted for truncation and the estimates for the truncation odds.
The quantiles for are estimated using
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
Q |
Vector of the corresponding quantile estimates. |
p |
The used exceedance probability. |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trMLE, trDTMLE, trProbMLE, trEndpointMLE, trTestMLE, trQuant, Quant
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=FALSE) # Large quantile of X trQuantMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, p=0.005, ylim=c(15,30)) # Large quantile of Y trQuantMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, p=0.005, ylim=c(0,300), Y=TRUE)# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Truncation odds dtmle <- trDTMLE(X, gamma=trmle$gamma, tau=trmle$tau, plot=FALSE) # Large quantile of X trQuantMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, p=0.005, ylim=c(15,30)) # Large quantile of Y trQuantMLE(X, gamma=trmle$gamma, tau=trmle$tau, DT=dtmle$DT, plot=TRUE, p=0.005, ylim=c(0,300), Y=TRUE)
Test between non-truncated Pareto-type tails (light truncation) and truncated Pareto-type tails (rough truncation).
trTest(data, alpha = 0.05, plot = TRUE, main = "Test for truncation", ...)trTest(data, alpha = 0.05, plot = TRUE, main = "Test for truncation", ...)
data |
Vector of |
alpha |
The used significance level, default is |
plot |
Logical indicating if the P-values should be plotted as a function of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
We want to test
has non-truncated Pareto tails vs.
has truncated Pareto tails.
Let
with the -th order statistic.
The test statistic is then
which is asymptotically standard normally distributed.
We reject on level if
where is the quantile of a standard normal distribution.
The corresponding P-value is thus given by
with the CDF of a standard normal distribution.
See Beirlant et al. (2016) or Section 4.2.3 of Albrecher et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
testVal |
Corresponding test values. |
critVal |
Critical value used for the test, i.e. |
Pval |
Corresponding P-values. |
Reject |
Logical vector indicating if the null hypothesis is rejected for a certain value of |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Beirlant, J., Fraga Alves, M.I. and Gomes, M.I. (2016). "Tail fitting for Truncated and Non-truncated Pareto-type Distributions." Extremes, 19, 429–462.
# Sample from truncated Pareto distribution. # truncated at 95% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.95, shape=shape)) # Test for truncation trTest(X) # Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Test for truncation trTest(X)# Sample from truncated Pareto distribution. # truncated at 95% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.95, shape=shape)) # Test for truncation trTest(X) # Sample from truncated Pareto distribution. # truncated at 99% quantile shape <- 2 X <- rtpareto(n=1000, shape=shape, endpoint=qpareto(0.99, shape=shape)) # Test for truncation trTest(X)
Test between non-truncated GPD tails (light truncation) and truncated GPD tails (rough truncation).
trTestMLE(data, gamma, tau, alpha = 0.05, plot = TRUE, main = "Test for truncation", ...)trTestMLE(data, gamma, tau, alpha = 0.05, plot = TRUE, main = "Test for truncation", ...)
data |
Vector of |
gamma |
Vector of |
tau |
Vector of |
alpha |
The used significance level, default is |
plot |
Logical indicating if the P-values should be plotted as a function of |
main |
Title for the plot, default is |
... |
Additional arguments for the |
We want to test
has non-truncated GPD tails vs.
has truncated GPD tails.
Let and be the truncated MLE estimates for and .
The test statistic is then
which is asymptotically standard exponentially distributed.
We reject on level if
. The corresponding P-value is given by
.
See Beirlant et al. (2017) for more details.
A list with following components:
k |
Vector of the values of the tail parameter |
testVal |
Corresponding test values. |
critVal |
Critical value used for the test, i.e. |
Pval |
Corresponding P-values. |
Reject |
Logical vector indicating if the null hypothesis is rejected for a certain value of |
Tom Reynkens.
Beirlant, J., Fraga Alves, M. I. and Reynkens, T. (2017). "Fitting Tails Affected by Truncation". Electronic Journal of Statistics, 11(1), 2026–2065.
trMLE, trDTMLE, trProbMLE, trEndpointMLE, trTestMLE, trTest
# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Test for truncation trTestMLE(X, gamma=trmle$gamma, tau=trmle$tau)# Sample from GPD truncated at 99% quantile gamma <- 0.5 sigma <- 1.5 X <- rtgpd(n=250, gamma=gamma, sigma=sigma, endpoint=qgpd(0.99, gamma=gamma, sigma=sigma)) # Truncated ML estimator trmle <- trMLE(X, plot=TRUE, ylim=c(0,2)) # Test for truncation trTestMLE(X, gamma=trmle$gamma, tau=trmle$tau)
Computes the Turnbull estimator for the survival function of interval censored data.
Turnbull(x, L, R, censored, trunclower = 0, truncupper = Inf, conf.type = "plain", conf.int = 0.95)Turnbull(x, L, R, censored, trunclower = 0, truncupper = Inf, conf.type = "plain", conf.int = 0.95)
x |
Vector with points to evaluate the estimator in. |
L |
Vector of length |
R |
Vector of length |
censored |
Vector of |
trunclower |
Lower truncation point, default is 0. |
truncupper |
Upper truncation point, default is |
conf.type |
Type of confidence interval, see |
conf.int |
Confidence level of the two-sided confidence interval, see |
We consider the random interval censoring model where one observes
and where the variable of interest lies between and .
Right censored data should be entered as L=l and R=truncupper, and right censored data should be entered as L=trunclower and R=r.
This function calls survfit.formula from survival.
See Section 4.3.2 in Albrecher et al. (2017) for more details.
A list with following components:
surv |
A vector of length |
fit |
The output from the call to |
Tom Reynkens
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Turnbull, B. W. (1974). "Nonparametric Estimation of a Survivorship Function with Doubly Censored Data." Journal of the American Statistical Association, 69, 169–173.
Turnbull, B. W. (1976). "The Empirical Distribution Function with Arbitrarily Grouped, Censored and Truncated Data." Journal of the Royal Statistical Society: Series B (Methodological), 38, 290–295.
L <- 1:10 R <- c(1, 2.5, 3, 4, 5.5, 6, 7.5, 8.25, 9, 10.5) censored <- c(0, 1, 0, 0, 1, 0, 1, 1, 0, 1) x <- seq(0, 12, 0.1) # Turnbull estimator plot(x, Turnbull(x, L, R, censored)$cdf, type="s", ylab="Turnbull estimator")L <- 1:10 R <- c(1, 2.5, 3, 4, 5.5, 6, 7.5, 8.25, 9, 10.5) censored <- c(0, 1, 0, 0, 1, 0, 1, 1, 0, 1) x <- seq(0, 12, 0.1) # Turnbull estimator plot(x, Turnbull(x, L, R, censored)$cdf, type="s", ylab="Turnbull estimator")
Density, distribution function, quantile function and random generation for the truncated Weibull distribution.
dtweibull(x, shape, scale = 1, endpoint = Inf, log = FALSE) ptweibull(x, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtweibull(p, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtweibull(n, shape, scale = 1, endpoint = Inf)dtweibull(x, shape, scale = 1, endpoint = Inf, log = FALSE) ptweibull(x, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) qtweibull(p, shape, scale = 1, endpoint = Inf, lower.tail = TRUE, log.p = FALSE) rtweibull(n, shape, scale = 1, endpoint = Inf)
x |
Vector of quantiles. |
p |
Vector of probabilities. |
n |
Number of observations. |
shape |
The shape parameter of the Weibull distribution, a strictly positive number. |
scale |
The scale parameter of the Weibull distribution, a strictly positive number, default is 1. |
endpoint |
Endpoint of the truncated Weibull distribution. The default value is |
log |
Logical indicating if the densities are given as |
lower.tail |
Logical indicating if the probabilities are of the form |
log.p |
Logical indicating if the probabilities are given as |
The Cumulative Distribution Function (CDF) of the truncated Weibull distribution is equal to
for where is the CDF of the ordinary Weibull distribution and is the endpoint (truncation point) of the truncated Weibull distribution.
dtweibull gives the density function evaluated in , ptweibull the CDF evaluated in and qtweibull the quantile function evaluated in . The length of the result is equal to the length of or .
rtweibull returns a random sample of length .
Tom Reynkens.
# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtweibull(x, shape=2, scale=0.5, endpoint=1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptweibull(x, shape=2, scale=0.5, endpoint=1), xlab="x", ylab="CDF", type="l")# Plot of the PDF x <- seq(0, 10, 0.01) plot(x, dtweibull(x, shape=2, scale=0.5, endpoint=1), xlab="x", ylab="PDF", type="l") # Plot of the CDF x <- seq(0, 10, 0.01) plot(x, ptweibull(x, shape=2, scale=0.5, endpoint=1), xlab="x", ylab="CDF", type="l")
Compute Value-at-Risk () of the fitted spliced distribution.
VaR(p, splicefit)VaR(p, splicefit)
p |
The exceedance probability (we estimate |
splicefit |
A |
See Reynkens et al. (2017) and Section 4.6 of Albrecher et al. (2017) for details.
Note that VaR(p, splicefit) corresponds to qSplice(p, splicefit, lower.tail = FALSE).
Vector of quantiles .
Tom Reynkens with R code from Roel Verbelen for the mixed Erlang quantiles.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
Reynkens, T., Verbelen, R., Beirlant, J. and Antonio, K. (2017). "Modelling Censored Losses Using Splicing: a Global Fit Strategy With Mixed Erlang and Extreme Value Distributions". Insurance: Mathematics and Economics, 77, 65–77.
Verbelen, R., Gong, L., Antonio, K., Badescu, A. and Lin, S. (2015). "Fitting Mixtures of Erlangs to Censored and Truncated Data Using the EM Algorithm." Astin Bulletin, 45, 729–758
qSplice, CTE, SpliceFit, SpliceFitPareto, SpliceFiticPareto, SpliceFitGPD
## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) p <- seq(0,1,0.01) # Plot of quantiles plot(p, qSplice(p, splicefit), type="l", xlab="p", ylab="Q(p)") # Plot of VaR plot(p, VaR(p, splicefit), type="l", xlab="p", ylab=bquote(VaR[1-p])) ## End(Not run)## Not run: # Pareto random sample X <- rpareto(1000, shape = 2) # Splice ME and Pareto splicefit <- SpliceFitPareto(X, 0.6) p <- seq(0,1,0.01) # Plot of quantiles plot(p, qSplice(p, splicefit), type="l", xlab="p", ylab="Q(p)") # Plot of VaR plot(p, VaR(p, splicefit), type="l", xlab="p", ylab=bquote(VaR[1-p])) ## End(Not run)
Computes the empirical quantiles of the log-transform of a data vector and the theoretical quantiles of the standard Weibull distribution. These quantiles are then plotted in a Weibull QQ-plot with the theoretical quantiles on the -axis and the empirical quantiles on the -axis.
WeibullQQ(data, plot = TRUE, main = "Weibull QQ-plot", ...)WeibullQQ(data, plot = TRUE, main = "Weibull QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the quantiles should be plotted in a Weibull QQ-plot, default is |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The Weibull QQ-plot is given by
for with the -th order statistic of the data.
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
wqq.the |
Vector of the theoretical quantiles from a standard Weibull distribution. |
wqq.emp |
Vector of the empirical quantiles from the log-transformed data. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
WeibullQQ_der, ExpQQ, LognormalQQ, ParetoQQ
data(norwegianfire) # Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ(norwegianfire$size[norwegianfire$year==76]) # Derivative of Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ_der(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ(norwegianfire$size[norwegianfire$year==76]) # Derivative of Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ_der(norwegianfire$size[norwegianfire$year==76])
Computes the derivative plot of the Weibull QQ-plot. These values can be plotted as a function of the data or as a function of the tail parameter .
WeibullQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of Weibull QQ-plot", ...)WeibullQQ_der(data, k = FALSE, plot = TRUE, main = "Derivative plot of Weibull QQ-plot", ...)
data |
Vector of |
plot |
Logical indicating if the derivative values should be plotted, default is |
k |
Logical indicating if the derivative values are plotted as a function of the tail parameter |
main |
Title for the plot, default is |
... |
Additional arguments for the |
The derivative plot of a Weibull QQ-plot is
or
with the Hill estimates and
See Section 4.1 of Albrecher et al. (2017) for more details.
A list with following components:
xval |
Vector of the x-values of the plot ( |
yval |
Vector of the derivative values. |
Tom Reynkens.
Albrecher, H., Beirlant, J. and Teugels, J. (2017). Reinsurance: Actuarial and Statistical Aspects, Wiley, Chichester.
WeibullQQ, Hill, MeanExcess, LognormalQQ_der, ParetoQQ_der
data(norwegianfire) # Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ(norwegianfire$size[norwegianfire$year==76]) # Derivative of Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ_der(norwegianfire$size[norwegianfire$year==76])data(norwegianfire) # Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ(norwegianfire$size[norwegianfire$year==76]) # Derivative of Weibull QQ-plot for Norwegian Fire Insurance data for claims in 1976. WeibullQQ_der(norwegianfire$size[norwegianfire$year==76])