evmix (version 1.0)

hpd: Hybrid Pareto Extreme Value Mixture Model

Description

Density, cumulative distribution function, quantile function and random number generation for the hybrid Pareto extreme value mixture model. The parameters are the normal mean nmean and standard deviation nsd and GPD shape xi.

Usage

dhpd(x, nmean = 0, nsd = 1, xi = 0, log = FALSE)

  phpd(q, nmean = 0, nsd = 1, xi = 0, lower.tail = TRUE)

  qhpd(p, nmean = 0, nsd = 1, xi = 0, lower.tail = TRUE)

  rhpd(n = 1, nmean = 0, nsd = 1, xi = 0)

Arguments

x

quantile

nmean

normal mean

nsd

normal standard deviation (non-negative)

xi

shape parameter

log

logical, if TRUE then log density

q

quantile

lower.tail

logical, if FALSE then upper tail probabilities

p

cumulative probability

n

sample size (non-negative integer)

Value

• dhpd gives the density, phpd gives the cumulative distribution function, qhpd gives the quantile function and rhpd gives a random sample.

Details

Extreme value mixture model combining normal distribution for the bulk below the threshold and GPD for upper tail which is continuous in its zeroth and first derivative at the threshold, but with one important difference to all the other mixture models. The hybrid Pareto does not include the usual tail fraction phiu scaling, i.e. so the GPD is not treated as a conditional model for the exceedances. The unscaled GPD is simply spliced with the normal truncated at the threshold, with no rescaling to account for the proportion above the threshold being applied. The parameters have to adjust for the lack of tail fraction scaling. The two continuity constraints lead to the threshold u and GPD scale sigmau being replaced by a function of the normal mean, standard deviation and GPD shape parameters. The continuity constraint on its first derivative at the threshold means that $h'(u) = g'(u)$. Then the Lambert W function is used for replacing the threshold u and GPD scale sigmau in terms of the normal mean, standard deviation and GPD shape xi. The cumulative distribution function defined upto the threshold $x \le u$, given by: $$F(x) = H(x) / r$$ and above the threshold $x > u$: $$F(x) = (H(u) + G(x)) / r$$ where $H(x)$ and $G(X)$ are the normal and conditional GPD cumulative distribution functions (i.e. pnorm(x, nmean, nsd), pgpd(x, u, sigmau, xi)). The normalisation constant $r$ ensures a proper density and is given byr = 1 + pnorm(u, mean = nmean, sd = nsd), i.e. the 1 comes from integration of the unscaled GPD and the second term is from the usual normal component. The continuity constraint on the density at the threshold means that $h(u) = g(u)$ where $h(x)$ and $g(x)$ are the normal and unscaled GPD density functions (i.e. dnorm(u, nmean, nsd) and dgpd(u, u, sigmau, xi)). See gpd for details of GPD upper tail component and dnorm for details of normal bulk component.

References

http://en.wikipedia.org/wiki/Normal_distribution http://en.wikipedia.org/wiki/Generalized_Pareto_distribution Scarrott, C.J. and MacDonald, A. (2012). A review of extreme value threshold estimation and uncertainty quantification. REVSTAT - Statistical Journal 10(1), 33-59. Available from http://www.ine.pt/revstat/pdf/rs120102.pdf Carreau, J. and Y. Bengio (2008). A hybrid Pareto model for asymmetric fat-tailed data: the univariate case. Extremes 12 (1), 53-76.

See Also

gpd, dnorm, dnormgpd and dnormgpdcon. The condmixt package written by one of the original authors of the hybrid Pareto model (Carreau and Bengio, 2008) also has similar functions for the hybrid Pareto hpareto and mixture of hybrid Paretos hparetomixt, which are more flexible as they also permit the model to be truncated at zero. Other hpd: fhpd, lhpd, nlhpd

Examples

par(mfrow = c(2, 2))
xx = seq(-5, 20, 0.01)
f1 = dhpd(xx, 0,1, 0.4)
plot(xx, f1, type = "l")
abline(v = 0.4942921)

# three tail behaviours
plot(xx, phpd(xx), type = "l")
lines(xx, phpd(xx, xi = 0.3), col = "red")
lines(xx, phpd(xx, xi = -0.3), col = "blue")
legend("bottomright", paste("xi =",c(0, 0.3, -0.3)),
  col=c("black", "red", "blue"), lty = 1)

sim = rhpd(10000, nmean = 0, nsd = 1.5, xi = 0.2)
hist(sim, freq = FALSE, 100, xlim = c(-5, 20), ylim = c(0, 0.2))
lines(xx, dhpd(xx, 0, 1.5, 0.2), col = "blue")

plot(xx, dhpd(xx, nmean = 0, nsd = 1.5, xi = 0), type = "l")
lines(xx, dhpd(xx, nmean = 0, nsd = 1.5, xi = 0.2), col = "red")
lines(xx, dhpd(xx, nmean = 0, nsd = 1.5, xi = -0.2), col = "blue")
legend("topright", c("xi = 0", "xi = 0.2", "xi = -0.2"),
  col=c("black", "red", "blue"), lty = 1)