| extBQuantx | R Documentation |
Given posterior samples for the parameters of the continuous or discrete generalized Pareto distribution and scedasis function for a set of covariates,
return the posterior mean and 1-\alpha level credibility intervals of the extreme quantile for each value of the scedasis function
extBQuantx(
cx,
postsamp,
threshold,
n,
qlev,
k,
type = c("continuous", "discrete"),
confint = c("asymmetric", "symmetric"),
alpha = 0.05,
...
)
cx |
an |
postsamp |
a |
threshold |
threshold for the generalized Pareto model, corresponding to the |
n |
integer, number of observations in the full sample. Must be greater than |
qlev |
double indicating the percentile level at which the extreme quantile is estimated. Must be smaller than |
k |
integer, number of exceedances for the generalized Pareto (only used if |
type |
string indicating distribution types. Default: |
confint |
type of confidence interval. Default: |
alpha |
level for credibility interval. Default 0.05, giving 95% credibility intervals |
... |
additional arguments, for back-compatilibity |
a list with components
mQ posterior mean of the extreme quantile
ciQ vector of dimension 2 returning the \alpha/2 and 1-\alpha/2 empirical quantiles of the posterior distribution of the extreme quantile
## Not run:
# generate data
set.seed(1234)
n <- 500
samp <- evd::rfrechet(n,0,1:n,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,
control = list(maxit = 500))
# 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)
# conditional predictive density estimation
yg <- seq(0, 50, by = 2)
nyg <- length(yg)
# estimation of scedasis function
# setting
M <- 1e3
C <- 5
alpha <- 0.05
bw <- .5
nsim <- 5000
burn <- 1000
# create covariate
# in sample obs
n_in = n
# number of years ahead
nY = 1
n_out = 365 * nY
# total obs
n_tot = n_in + n_out
# total covariate (in+out sample period)
x <- seq(0, 1, length = n_tot)
# in sample grid dimension for covariate
ng_in <- 150
xg <- seq(0, x[n_in], length = ng_in)
# in+out of sample grid
xg <- c(xg, seq(x[n_in + 1], x[(n_tot)], length = ng_in))
# in+out sample grid dimension
nxg <- length(xg)
xg <- array(xg, c(nxg, 1))
# in sample observations
samp_in <- samp[1:n_in]
ssamp_in <- sort(samp_in, decreasing = TRUE, index = TRUE)
x_in <- x[1:n_in] # in sample covariate
xs <- x_in[ssamp_in$ix[1:k]] # in sample concomitant covariate
# estimate scedasis function over the in and out of sample period
res_stat <- apply(
xg,
1,
cpost_stat,
N = nsim - burn,
x = x_in,
xs = xs,
bw = bw,
k = k,
C = C
)
# conditional predictive posterior density
probq = 1 - 0.99
res_AQ <- extBQuantx(
cx = res_stat,
postsamp = proc$post_sample,
threshold = proc$t,
n = n,
qlev = probq,
k = k,
type = "continuous",
confint = "asymmetric")
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.