inst/examples.R

In ExtremeRisks: Extreme Risk Measures

library(ExtremeRisks)

#### fitdGPD----
rm(list=ls())
fitdGPD(rpois(1000,2))

#### estPOT----
rm(list=ls())
# 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)

#### extBQuant----
rm(list=ls())
# generate data
set.seed(1234)
n <- 500
samp <- evd::rfrechet(500,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)
# extreme quantile corresponding to a return period of 100
extBQuant(
  proc$t,proc$post_sample,
  k,
  n,
  100,
  0.05,
  type = "continuous")


#### predDens----
rm(list=ls())
# 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"))



#### predQuant----
rm(list=ls())
# 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)
predQuant_int <- predQuant(
  0.5,
  proc$post_sample,
  proc$t,
  proc$t + 0.01,
  50,
  "continuous",
  extrapolation = FALSE)
predDens_ext <- predDens(
  yg,
  proc$post_sample,
  proc$t,
  "continuous",
  extrapolation = TRUE,
  p = 0.001,
  k = k,
  n = n)
predQuant_ext <- predQuant(
  0.5,
  proc$post_sample,
  proc$t,
  proc$t + 0.01,
  100,
  "continuous",
  extrapolation = TRUE,
  p = 0.005,
  k = k,
  n = n)


#### scedastic.test----
# generate data
rm(list=ls())
# 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
# test on covariate effect
test <- scedastic.test(
  cbind(samp, x[1:n]),
  k,
  M,
  array(xg[1:ng_in], c(ng_in, 1)),
  ng_in,
  TRUE,
  C,
  0.05
)


#### cpost_stat----
rm(list=ls())
# 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
)




#### extBQuantx----
rm(list=ls())
# generate data# 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")




#### predDensx----
rm(list=ls())
# 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
yg <- seq(500, 5000, by = 100)
nyg = length(yg)
# intermediate threshold
predDens_intx <- predDensx(
  x = yg,
  postsamp = proc$post_sample,
  scedasis = res_stat,
  threshold = proc$t,
  "continuous",
  extrapolation = FALSE)
# extreme threshold
predDens_extx <- predDensx(
  x = yg,
  postsamp = proc$post_sample,
  scedasis = res_stat,
  threshold = proc$t,
  "continuous",
  extrapolation = TRUE,
  p = 0.001,
  k = k,
  n = n)
# plot intermediate and extreme density conditional on a specific value of scedasis function
# disclaimer: to speed up the procedure, we specify a coarse grid
plot(
 x = predDens_intx$x,
 y = predDens_intx$preddens[,20],
 type = "l",
 lwd = 2,
 col="dodgerblue",
 ylab = "",
 xlab = "yg",
 main = "Conditional predictive posterior density")
lines(
  x = predDens_extx$x,
  y = predDens_extx$preddens[,20],
  lwd = 2,
  col = "violet")
legend("topright",
  legend = c("Intermediate threshold","Extreme threshold"),
  lwd = 2,
  col = c("dodgerblue", "violet"))






#### quantF----
rm(list=ls())
# 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
yg <- seq(500, 5000, by = 100)
nyg = length(yg)
# intermediate threshold
predDens_intx <- predDensx(
  x = yg,
  postsamp = proc$post_sample,
  scedasis = res_stat,
  threshold = proc$t,
  "continuous",
  extrapolation = FALSE)
# extreme threshold
predDens_extx <- predDensx(
  x = yg,
  postsamp = proc$post_sample,
  scedasis = res_stat,
  threshold = proc$t,
  "continuous",
  extrapolation = TRUE,
  p = 0.001,
  k = k,
  n = n)
  # predictive conditional quantiles
predQuant_intxL <- predQuant_intxU <- predQuant_extxL <- predQuant_extxU <- numeric()
for (i in 1:nxg){
  predQuant_intxU[i] <- apply(
    array(predDens_intx$preddens[,i], c(1, nyg)),
    1,
    quantF,
    x = yg,
    p = c(0.975)
  )
  predQuant_intxL[i] <- apply(
    array(predDens_intx$preddens[,i], c(1, nyg)),
    1,
    quantF,
    x = yg,
    p = c(0.025)
  )
  predQuant_extxU[i] <- apply(
    array(predDens_extx$preddens[,i], c(1, nyg)),
    1,
    quantF,
    x = yg,
    p = c(0.975)
  )
  predQuant_extxL[i] <- apply(
    array(predDens_extx$preddens[,i], c(1, nyg)),
    1,
    quantF,
    x = yg,
    p = c(0.025)
  )
}





#### testTailHomo----
rm(list=ls())
# generate two samples of data
set.seed(1234)
y1 <- evd::rgpd(500, 0, 1, 0.2)
y2 <- evd::rgpd(700, 0, 2, 0.7)
y <- list(y1 = y1, y2 = y2)
# set effective sample size
k <- 50
# perform test
test <- testTailHomo(y,k)





#### plotBayes----
rm(list=ls())
# 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)
plotBayes(
  proc$post_sample[,1],
  mlest$estimate,
  param = "scale")
plotBayes(
  proc$post_sample[,2],
  param = "shape")