ExtremeRisks (version 0.0.5)

predDens: Predictive posterior density of peak-over-threshold models

Description

Given posterior samples for the parameters of the continuous or discrete generalized Pareto distribution, return the predictive posterior density of a peak above an intermediate or extreme threshold using the threshold stability property.

Usage

predDens(
  x,
  postsamp,
  threshold,
  type = c("continuous", "discrete"),
  extrapolation = FALSE,
  p,
  k,
  n
)

Value

a vector of length r of posterior predictive density values associated to x

Arguments

x

vector of length r containing the points at which to evaluate the density

postsamp

a m by 2 matrix with columns containing the posterior samples of scale and shape parameters of the generalized Pareto distribution, respectively

threshold

threshold for the generalized Pareto model, corresponding to the \(n-k\)th order statistic of the sample

type

data type, either "continuous" or "discrete" for count data.

extrapolation

logical; if TRUE, extrapolate using the threshold stability property

p

scalar tail probability for the extrapolation. Must be smaller than \(k/n\) (only used if extrapolation=TRUE)

k

integer, number of exceedances for the generalized Pareto (only used if extrapolation=TRUE)

n

integer, number of observations in the full sample. Must be greater than k (only used if extrapolation=TRUE)

Examples

if (FALSE) {
# generate data
set.seed(1234)
n <- 500
samp <- evd::rfrechet(n,0,3,4)
# set effective sample size and threshold
k <- 50
threshold <- sort(samp,decreasing = TRUE)[k+1]
# preliminary mle estimates of scale and shape parameters
mlest <- evd::fpot(samp, threshold)
# empirical bayes procedure
proc <- estPOT(
  samp,
  k = k,
  pn = c(0.01, 0.005),
  type = "continuous",
  method = "bayesian",
  prior = "empirical",
  start = as.list(mlest$estimate),
  sig0 = 0.1)
# predictive density estimation
yg <- seq(0, 50, by = 2)
nyg <- length(yg)
predDens_int <- predDens(
  yg,
  proc$post_sample,
  proc$t,
  "continuous",
  extrapolation = FALSE)
predDens_ext <- predDens(
  yg,
  proc$post_sample,
  proc$t,
  "continuous",
  extrapolation = TRUE,
  p = 0.001,
  k = k,
  n = n)
# plot
plot(
 x = yg,
 y = predDens_int,
 type = "l",
 lwd = 2,
 col = "dodgerblue",
 ylab = "",
 main = "Predictive posterior density")
lines(
 x = yg,
 y = predDens_ext,
 lwd = 2,
 col = "violet")
legend(
  "topright",
  legend = c("Intermediate threshold", "Extreme threshold"),
  lwd = 2,
  col = c("dodgerblue", "violet"))
  }