#' Exceedance Residuals Backtest
#'
#' Tests whether the mean of the exceedance residuals, respectively the
#' mean of the standardized exceedance residuals is zero.
#'
#' @inheritParams parameter_definition
#' @param B Number of bootstrap iterations
#' @return Returns a list with the following components:
#' * pvalue_twosided_simple
#' * pvalue_onesided_simple
#' * pvalue_twosided_standardized
#' * pvalue_onesided_standardized
#' @examples
#' data(risk_forecasts)
#' r <- risk_forecasts$r
#' q <- risk_forecasts$q
#' e <- risk_forecasts$e
#' s <- risk_forecasts$s
#' er_backtest(r = r, q = q, e = e, s = s)
#' @references McNeil & Frey (2000) \doi{10.1016/S0927-5398(00)00012-8}
#' @export
#' @md
er_backtest <- function(r, q, e, s = NULL, B = 1000) {
er_backtest_fun <- function(x) {
set.seed(1)
boot_x <- matrix(sample(x, size = length(x) * B, replace = TRUE), nrow = B)
f <- function(x) mean(x) / stats::sd(x) * sqrt(length(x))
t0 <- f(x)
t <- apply(boot_x, 1, f)
t <- t[is.finite(t) & !is.na(t)]
c(pv_2s = mean(abs(t - mean(t)) >= abs(t0)),
pv_1s = mean(t - mean(t) <= t0))
}
# Get the backtest p-values
er_simple <- er_backtest_fun((r - e)[r <= q])
if (!is.null(s)) {
er_standardized <- er_backtest_fun(((r - e) / s)[r <= q])
} else {
er_standardized <- c(NA, NA)
}
# Return results
ret <- list(
pvalue_twosided_simple = unname(er_simple[1]),
pvalue_onesided_simple = unname(er_simple[2]),
pvalue_twosided_standardized = unname(er_standardized[1]),
pvalue_onesided_standardized = unname(er_standardized[2])
)
ret
}
#' Conditional Calibration Backtest
#'
#' The simple and general conditional calibration backtests of Nolde & Ziegel (2007).
#'
#' @inheritParams parameter_definition
#' @param hommel If TRUE, use Hommels correction,
#' otherwise use the classical Bonferroni correction.
#' @return Returns a list with the following components:
#' * pvalue_twosided_simple
#' * pvalue_onesided_simple
#' * pvalue_twosided_general
#' * pvalue_onesided_general
#' @examples
#' data(risk_forecasts)
#' r <- risk_forecasts$r
#' q <- risk_forecasts$q
#' e <- risk_forecasts$e
#' s <- risk_forecasts$s
#' cc_backtest(r = r, q = q, e = e, s = s, alpha = 0.025)
#' @references Nolde & Ziegel (2007) \doi{10.1214/17-AOAS1041}
#' @export
#' @md
cc_backtest <- function(r, q, e, s=NULL, alpha, hommel=TRUE) {
# Sample length
n <- length(r)
# Identification function matrix (n x 2)
V <- cbind(alpha - (r <= q),
e - q + (r <= q) * (q - r) / alpha)
# Test functions; n x q x 2
H1 <- aperm(replicate(n, diag(2))) # q = 2
if (!is.null(s)) {
H2 <- array(cbind((q - e) / alpha / s, 1 / s), dim = c(n, 1, 2)) # q = 1
H3 <- array(NA, dim = c(n, 4, 2)) # q = 4
H3[,, 1] <- cbind(1, abs(q), 0, 0)
H3[,, 2] <- cbind(0, 0, 1, 1 / s)
}
# n x q matrices of h * V
hV1 <- apply(H1, 2, function(x) rowSums(x * V))
if (!is.null(s)) {
hV2 <- apply(H2, 2, function(x) rowSums(x * V))
hV3 <- apply(H3, 2, function(x) rowSums(x * V))
}
# Estimates of the covariance matrix of h * V
omega1 <- crossprod(hV1) / n
if (!is.null(s)) {
omega2 <- crossprod(hV2) / n
omega3 <- crossprod(hV3) / n
}
# Test statistics and p-values
t1 <- as.numeric(n * colMeans(hV1) %*% solve(omega1) %*% colMeans(hV1))
t3 <- sqrt(n) * diag(omega1)^(-1/2) * colMeans(hV1)
p1 <- 1 - stats::pchisq(t1, 2)
p3 <- ifelse(hommel,
min(2 * sum(1 / (1:2)) * min(sort(1 - stats::pnorm(t3)) / (1:2)), 1),
min(2 * min(1 - stats::pnorm(t3)), 1))
if (!is.null(s)) {
t2 <- as.numeric(n * colMeans(hV2) %*% solve(omega2) %*% colMeans(hV2))
t4 <- sqrt(n) * diag(omega3)^(-1/2) * colMeans(hV3)
p2 <- 1 - stats::pchisq(t2, 1)
p4 <- ifelse(hommel,
min(4 * sum(1 / (1:4)) * min(sort(1 - stats::pnorm(t4)) / (1:4)), 1),
min(4 * min(1 - stats::pnorm(t4)), 1))
} else {
p2 <- NA
p4 <- NA
}
# Return results
ret <- list(
pvalue_twosided_simple = p1,
pvalue_onesided_simple = p3,
pvalue_twosided_general = p2,
pvalue_onesided_general = p4
)
ret
}
#' Expected Shortfall Regression Backtest
#'
#' This function implements multiple expected shortfall regression (esreg)
#' based backtests.
#' Using the `version` argument, the following backtests are available:
#' 1. ("Strict ESR") Regresses the returns on the expected shortfall forecasts
#' and tests the ES coefficients for (0, 1).
#' 1. ("Auxiliary ESR") Regresses the returns on the quantile and the expected shortfall forecasts
#' and tests the ES coefficients for (0, 1).
#' 1. ("Strict Intercept") Tests whether the expected shortfall of the forecast error r - e is zero.
#'
#' @inheritParams parameter_definition
#' @param version Version of the backtest to be used
#' @param cov_config a list with three components: sparsity, sigma_est, and misspec, see \link[esreg]{vcovA}
#' @return Returns a list with the following components:
#' * pvalue_two_sided_asymptotic
#' * pvalue_one_sided_asymptotic (for version 3)
#' * pvalue_two_sided_bootstrap
#' * pvalue_one_sided_bootstrap (for version 3)
#' @examples
#' data(risk_forecasts)
#' r <- risk_forecasts$r
#' q <- risk_forecasts$q
#' e <- risk_forecasts$e
#' esr_backtest(r = r, q = q, e = e, alpha = 0.025, version = 1)
#' @references Bayer & Dimitriadis (2020) \doi{10.1093/jjfinec/nbaa013}
#' @export
#' @md
esr_backtest <- function(r, q, e, alpha, version, B = 0,
cov_config=list(sparsity='nid', sigma_est='scl_sp', misspec=TRUE)) {
if (missing(q) & version %in% c(2)) {
stop('You need to supply VaR forecast `q` for backtest version ', version)
}
if (missing(q)) {
data <- data.frame(r = r, e = e)
} else {
data <- data.frame(r = r, q = q, e = e)
}
# Set the details for the selected version of the backtest
if (version == 1) {
model <- r ~ e
h0 <- c(NA, NA, 0, 1)
one_sided <- FALSE
} else if (version == 2) {
model <- r ~ q | e
h0 <- c(NA, NA, 0, 1)
one_sided <- FALSE
} else if (version == 3) {
model <- I(r - e) ~ e | 1
h0 <- c(NA, NA, 0)
one_sided <- TRUE
} else {
stop('This is a non-supported backtest version!')
}
# Fit the model
fit0 <- esreg::esreg(model, data = data, alpha = alpha, g1 = 2, g2 = 1)
cov0 <- esreg::vcovA(fit0,
sparsity = cov_config$sparsity,
sigma_est = cov_config$sigma_est,
misspec = cov_config$misspec)
s0 <- fit0$coefficients - h0
mask <- !is.na(h0)
# Compute the asymptotic test statistic and p-value
if (version %in% c(1, 2)) {
t0 <- as.numeric(s0[mask] %*% solve(cov0[mask, mask]) %*% s0[mask])
pv0_1s <- NA
pv0_2s <- 1 - stats::pchisq(t0, sum(mask))
} else if (version %in% c(3)) {
t0 <- as.numeric(s0[mask] / sqrt(cov0[mask, mask]))
pv0_1s <- stats::pnorm(t0)
pv0_2s <- 2 * (1 - stats::pnorm(abs(t0)))
}
# Compute the bootstrap p-values
if (B > 0) {
n <- length(r)
idx <- matrix(sample(1:n, n*B, replace=TRUE), nrow=n)
bs_estimates <- apply(idx, 2, function(id) {
tryCatch({
fitb <- esreg::esreg(model, data = data[id,], alpha = alpha, g1 = 2, g2 = 1, early_stopping = 0)
sb <- fitb$coefficients - fit0$coefficients
covb <- esreg::vcovA(fitb,
sparsity = cov_config$sparsity,
sigma_est = cov_config$sigma_est,
misspec = cov_config$misspec)
list(sb = sb, covb = covb)
}, error=function(e) NA)
})
idx_na <- is.na(bs_estimates)
share_na <- mean(idx_na)
if (share_na >= 0.05) stop('More than 5% of the bootstrap replications failed!')
bs_estimates <- bs_estimates[!idx_na]
if (version %in% c(1, 2)) {
tb <- sapply(bs_estimates, function(x) {
as.numeric(x$sb[mask] %*% solve(x$covb[mask, mask]) %*% x$sb[mask])
})
tb <- tb[!is.na(tb)]
pvb_2s <- mean(tb >= t0)
pvb_1s <- NA
} else if (version %in% c(3)) {
tb <- sapply(bs_estimates, function(x) {
x$sb[mask] / sqrt(x$covb[mask, mask])
})
tb <- tb[!is.na(tb)]
pvb_2s <- mean(abs(t0) <= abs(tb))
pvb_1s <- mean(tb <= t0)
}
} else {
pvb_2s <- NA
pvb_1s <- NA
}
# Return results
ret <- list(
pvalue_twosided_asymptotic = pv0_2s,
pvalue_twosided_bootstrap = pvb_2s
)
if (version %in% c(3)) {
ret['pvalue_onesided_asymptotic'] <- pv0_1s
ret['pvalue_onesided_bootstrap'] <- pvb_1s
}
ret
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.