#' Joint Quantile and Expected Shortfall Regression
#'
#' Estimates a joint linear regression model for the pair (VaR, ES):
#' \deqn{Q_\alpha(Y | Xq) = Xq'\beta_q}
#' \deqn{ES_\alpha(Y | Xe) = Xe'\beta_e}
#'
#' @param formula Formula: y ~ x1 + x2 ... | x1 + x2 ...
#' where the first part after the response variable specifies the quantile equation
#' and the second the expected shortfall part. If only one set of regressors is
#' provided it is used for both model specifications.
#' @param data data.frame that holds the variables
#' @param y Response vector
#' @param xq Explanatory variables for the quantile regression equation
#' @param xe Explanatory variables for the expected shortfall regression equation
#' @param alpha Probability level
#' @param g1 1, 2 (see \link{G1_fun}, \link{G1_prime_fun}), defaults to 1
#' @param g2 1, 2, 3, 4, 5 (see \link{G2_curly_fun}, \link{G2_fun}, \link{G2_prime_fun}). defaults to 2
#' @param early_stopping Stop the iterated local search if there is no improvement in early_stopping steps.
#' @param ... Further arguments (does not apply here)
#' @return An esreg object
#' @seealso \code{\link{vcov.esreg}} for covariance estimation
#' @examples
#' # Simulate data (DGP-(2) in the linked paper)
#' set.seed(0)
#' x <- rchisq(1000, df=1)
#' y <- -x + (1 + 0.5 * x) * rnorm(1000)
#'
#' # True quantile and expected shortfall regression parameters (for alpha=0.025)
#' alpha=0.025
#' true_pars <- c(-1.959964, -1.979982, -2.337803, -2.168901)
#'
#' # Estimate the model using the standard settings
#' fit <- esreg(y ~ x, alpha=alpha)
#'
#' # Compare the different variance-covariance estimators
#' cov1 <- vcov(object=fit, sparsity="iid", sigma_est="ind")
#' cov2 <- vcov(object=fit, sparsity="nid", sigma_est="scl_N")
#' cov3 <- vcov(object=fit, sparsity="nid", sigma_est="scl_sp")
#'
#' print("Comparison of the variance-covariance estimators")
#' print(cbind(Truth=true_pars,
#' Estimate=coef(fit),
#' SE_iid_ind=sqrt(diag(cov1)),
#' SE_nid_N=sqrt(diag(cov2)),
#' SE_nid_sp=sqrt(diag(cov3))))
#'
#' # Compares estimates using different G2 functions
#' fit1 <- esreg(y ~ x, alpha=alpha, g2=1)
#' fit2 <- esreg(y ~ x, alpha=alpha, g2=2)
#' fit3 <- esreg(y ~ x, alpha=alpha, g2=3)
#' fit4 <- esreg(y ~ x, alpha=alpha, g2=4)
#' fit5 <- esreg(y ~ x, alpha=alpha, g2=5)
#' fits <- sapply(list(fit1, fit2, fit3, fit4, fit5), coef)
#' colnames(fits) <- sapply(1:5, function(i) esreg:::.G_function_names(1, i)[2])
#' print("Comparison of the five G2 functions")
#' print(rbind(Truth=true_pars, t(fits)))
#'
#' # Usage of different covariates
#' x <- rchisq(1000, df=1)
#' noise <- rnorm(1000)
#' y <- -x + (1 + 0.5 * x) * rnorm(1000)
#' fit <- esreg(y ~ x | x + noise, alpha=0.025)
#' print("Using different covariates for VaR and ES")
#' print(summary(fit))
#'
#' @references \href{https://arxiv.org/abs/1704.02213}{A Joint Quantile and Expected Shortfall Regression Framework}
#' @rdname esreg
#' @export
esreg <- function(...) UseMethod('esreg')
#' @rdname esreg
#' @export
esreg.formula <- function(formula, data=parent.frame(), alpha, g1 = 2L, g2 = 1L,
early_stopping = 10, ...) {
chkDots(...)
# Prepare the formula
if(is.matrix(data)) {
data <- as.data.frame(data)
}
formula <- Formula::Formula(formula)
if (length(formula)[2] == 1) {
formula <- Formula::as.Formula(formula(formula),
formula(formula, lhs = 0, rhs = 1))
} else if (length(formula)[2] > 2) {
stop('You must not have more than two formula objects.')
}
mf <- stats::model.frame(formula = formula, data = data, na.action = stats::na.pass)
terms <- stats::terms(formula)
# Extract the data matrices
xq <- stats::model.matrix(formula, mf, rhs=1)
xe <- stats::model.matrix(formula, mf, rhs=2)
y <- stats::model.response(mf)
# Fit the model
fit <- esreg.fit(xq, xe, y, alpha, g1, g2, early_stopping)
fit$call <- match.call()
fit$terms <- terms
fit$formula <- formula
fit$data <- data
fit
}
#' @rdname esreg
#' @export
esreg.default <- function(xq, xe, y, alpha, g1 = 2L, g2 = 1L,
early_stopping = 10, ...) {
chkDots(...)
xq <- cbind("(Intercept)" = 1, xq)
xe <- cbind("(Intercept)" = 1, xe)
# Fit the model
fit <- esreg.fit(xq, xe, y, alpha, g1, g2, early_stopping)
fit$call <- match.call()
fit
}
#' @export
model.matrix.esreg <- function(object, ...) {
chkDots(...)
mm1 <- stats::model.matrix(object$formula, rhs=1, data=object$data)
mm2 <- stats::model.matrix(object$formula, rhs=2, data=object$data)
mm <- cbind(mm1, mm2)
attributes(mm)$assign <- c(attr(mm1, 'assign'), attr(mm2, 'assign'))
mm
}
#' @export
print.esreg <- function(x, digits = 4, ...) {
chkDots(...)
cat("Call:\n")
cat(deparse(x$call), "\n")
cat("\nQuantile Coefficients:\n")
print(format(x$coefficients_q, digits = digits), quote = FALSE)
cat("\nExpected Shortfall Coefficients:\n")
print(format(x$coefficients_e, digits = digits), quote = FALSE)
}
#' @export
summary.esreg <- function(object, ...) {
chkDots(...)
xq <- object$xq
xe <- object$xe
n <- nrow(xq)
k <- ncol(xq) + ncol(xe)
cov <- vcov.esreg(object, ...)
se <- sqrt(diag(cov))
tval <- stats::coef(object) / se
coef_mat <- cbind(
Estimate = stats::coef(object),
`Std. Error` = se,
`t value` = tval,
`Pr(>|t|)` = 2 * stats::pt(abs(tval), n - 2*k, lower.tail = FALSE)
)
structure(
c(object, list(cov = cov, coef_mat = coef_mat)),
class = "summary.esreg")
}
#' @export
print.summary.esreg <- function(x, ...) {
kq <- length(x$coefficients_q)
ke <- length(x$coefficients_e)
cat("Call:\n", paste0(deparse(x$call), sep = "\n", collapse = "\n"))
cat("\nQuantile Coefficients:\n")
stats::printCoefmat(x$coef_mat[1:kq,,drop=FALSE], signif.legend = FALSE)
cat("\nExpected Shortfall Coefficients:\n")
stats::printCoefmat(x$coef_mat[(kq+1):(kq+ke),,drop=FALSE])
}
#' @export
fitted.esreg <- function(object, ...) {
chkDots(...)
cbind(object$xq %*% object$coefficients_q,
object$xe %*% object$coefficients_e)
}
#' @export
residuals.esreg <- function(object, ...) {
chkDots(...)
cbind(object$y - object$xq %*% object$coefficients_q,
object$y - object$xe %*% object$coefficients_e)
}
#' @export
predict.esreg <- function(object, newdata=NULL, ...) {
chkDots(...)
if (is.null(newdata)) {
yhat <- fitted.esreg(object)
} else {
if (nrow(newdata) == 0)
stop('newdata is empty')
if (is.null(object$formula)) {
stop('The predict method is only supported using the formula interface')
} else {
mf <- stats::model.frame(object$formula, newdata)
xq <- stats::model.matrix(object$formula, mf, rhs=1)
xe <- stats::model.matrix(object$formula, mf, rhs=2)
yhat <- cbind(xq %*% object$coefficients_q,
xe %*% object$coefficients_e)
}
}
yhat
}
#' Estimating function
#'
#' This function matches the estfun function of the sandwich package and
#' returns the estimating functions for the fitted model.
#' It can for instance be used for an OPG estimator of the sigma matrix.
#' For esreg, the dimension of the estimating functions is n x (kq + ke).
#'
#' @param x An \link{esreg} object
#' @param ... Further arguments (does not apply here)
#' @export
estfun.esreg <- function(x, ...) {
chkDots(...)
# Extract elements from object
object <- x
y <- object$y
xq <- object$xq
xe <- object$xe
coefficients_q <- object$coefficients_q
coefficients_e <- object$coefficients_e
alpha <- object$alpha
n <- length(y)
kq <- ncol(xq)
ke <- ncol(xe)
# Transform the data and coefficients
if (object$g2 %in% c(1, 2, 3)) {
max_y <- max(y)
y <- y - max_y
coefficients_q[1] <- coefficients_q[1] - max_y
coefficients_e[1] <- coefficients_e[1] - max_y
}
# Precompute some quantities
xbq <- as.numeric(xq %*% coefficients_q)
xbe <- as.numeric(xe %*% coefficients_e)
uq <- as.numeric(y - xbq)
# Evaluate G1 / G2 functions
G1_prime_xq <- G_vec(z = xbq, g = "G1_prime", type = object$g1)
G2_xe <- G_vec(z = xbe, g = "G2", type = object$g2)
G2_prime_xe <- G_vec(z = xbe, g = "G2_prime", type = object$g2)
G2_prime_prime_xe <- G_vec(z = xbe, g = "G2_prime_prime", type = object$g2)
psi <- estimating_function_loop(
y = y, xq = xq, xe = xe, xbq = xbq, xbe = xbe, alpha = alpha,
G1_prime_xq = G1_prime_xq,
G2_xe = G2_xe, G2_prime_xe = G2_prime_xe)
psi
}
#' Covariance Estimation
#'
#' Estimate the variance-covariance matrix of the joint (VaR, ES) estimator
#'
#' @param object An \link{esreg} object
#' @param method For asymptotic use \link{vcovA}, for boot use \link{vcovB}
#' @param ... All possible values which can be passed to \link{vcovA} and \link{vcovB}
#' @export
vcov.esreg <- function(object, method='asymptotic', ...) {
if (method == 'asymptotic') {
cov <- vcovA(object, ...)
} else if(method == 'boot') {
cov <- vcovB(object, ...)
} else {
stop('method can be asymptotic or boot')
}
cov
}
#' Asymptotic Covariance Estimation
#'
#' Estimate the variance-covariance matrix of the joint (VaR, ES) estimator by the sandwich formula:
#' \deqn{\lambda^{-1} \Sigma \lambda^{-1}}
#' Several estimators are available for both matrices and the default options are selected to take into account
#' possible misspecifications in the underlying data.
#'
#' @param object An esreg object
#' @param sigma_est The estimator to be used for \eqn{\Sigma}, see \link{conditional_truncated_variance}
#' \itemize{
#' \item ind - Variance over all negative residuals
#' \item scl_N - Scaling with the normal distribution
#' \item scl_sp - Scaling with the kernel density function
#' }
#' @param sparsity The estimator to be used for the sparsity in \eqn{\Lambda}, see \link{density_quantile_function}
#' \itemize{
#' \item iid - Piecewise linear interpolation of the distribution
#' \item nid - Hendricks and Koenker sandwich
#' }
#' @param bandwidth_estimator The bandwidth estimator to be used for the iid and nid sparsity estimator, see \link{density_quantile_function}
#' \itemize{
#' \item Bofinger
#' \item Chamberlain
#' \item Hall-Sheather
#' }
#' @param misspec if TRUE, the estimator accounts for potential misspecification in the model
#' @export
vcovA <- function(object, sigma_est = 'scl_sp', sparsity = 'nid', misspec = TRUE, bandwidth_estimator = 'Hall-Sheather') {
lambda <- lambda_matrix(object = object, sparsity = sparsity,
bandwidth_estimator = bandwidth_estimator, misspec = misspec)
lambda_inverse <- solve(lambda)
sigma <- sigma_matrix(object = object, sigma_est = sigma_est, misspec = misspec)
n <- length(object$y)
cov <- 1/n * (lambda_inverse %*% sigma %*% lambda_inverse)
rownames(cov) <- colnames(cov) <- names(stats::coef(object))
cov
}
#' @title Bootstrap Covariance Estimation
#' @description Estimate the variance-covariance matrix of the joint (VaR, ES) estimator
#' using the bootstrap.
#' @param object An esreg object
#' @param bootstrap_method The bootstrap sampling scheme to be used
#' \itemize{
#' \item iid - The iid bootstrap of Efron (1979)
#' }
#' @param B The number of bootstrap iterations
#' @export
vcovB <-function(object, bootstrap_method='iid', B=1000) {
if (!(bootstrap_method %in% c("iid")))
stop("bootstrap_method can be iid")
if (B < 1000)
warning("The number of bootstrap iterations is small!")
# Draw the bootstrap indices
n <- length(object$y)
idx <- matrix(sample(1:n, size = n * B, replace = TRUE), nrow = n)
# Get the bootstrap parameter estimates
b <- sapply(seq_len(B), function(i) {
fitb <- esreg(object$y[idx[,i]] ~ object$xq[idx[,i], -1] | object$xe[idx[,i], -1],
alpha = object$alpha, g1 = object$g1, g2 = object$g2,
early_stopping = 0)
fitb$coefficients
})
# Compute the covariance
cov <- cov(t(b))
rownames(cov) <- colnames(cov) <- names(stats::coef(object))
cov
}
#' @export
esreg.fit <- function(xq, xe, y, alpha, g1, g2, early_stopping, ...) {
t0 <- Sys.time()
# Check input parameters and data
if (!(g1 %in% c(1, 2))) stop("G1 can be 1 or 2.")
if (!(g2 %in% c(1, 2, 3, 4, 5))) stop("G2 can be 1, 2, 3, 4 or 5.")
if ((alpha < 0) | (1 < alpha)) stop("alpha not in (0, 1)")
if (any(is.na(y)) | any(is.na(xq)) | any(is.na(xe)))
stop("Data contains NAs!")
if (!(all(is.finite(y)) & all(is.finite(xq)) & all(is.finite(xe))))
stop("Not all values are finite!")
if (length(y) != nrow(xq))
stop("Dimensions of y and xq do not match!")
if (length(y) != nrow(xe))
stop("Dimensions of y and xe do not match!")
# Transform the data
if (g2 %in% c(1, 2, 3)) {
max_y <- max(y)
y <- y - max_y
}
## Find starting values
# Match quantile and expected shortfall levels
e <- -stats::dnorm(stats::qnorm(alpha))/alpha
alpha_tilde <- stats::uniroot(function(x) stats::qnorm(x) - e, c(0, alpha))$root
# Fit the quantile regression
fit_qr_q <- suppressWarnings(quantreg::rq(y ~ xq - 1, tau = alpha))
fit_qr_e <- suppressWarnings(quantreg::rq(y ~ xe - 1, tau = alpha_tilde))
# Get standard errors for the perturbation
qr_sum_q <- suppressWarnings(summary(fit_qr_q, se="iid"))
qr_sum_e <- suppressWarnings(summary(fit_qr_e, se="iid"))
se <- c(qr_sum_q$coefficients[,2], qr_sum_e$coefficients[,2])
# Store the coefficient estimates
b0 <- c(fit_qr_q$coef, fit_qr_e$coef)
# If G2 is 1, 2, 3 ensure that x'be < 0 by moving the es intercept
kq <- ncol(xq)
ke <- ncol(xe)
if (g2 %in% c(1, 2, 3)) {
max_xe <- max(xe %*% b0[(kq+1):(kq+ke)])
b0[kq+1] <- b0[kq+1] - (max_xe + 0.1) * (max_xe >= -0.1)
}
# Set the target function
fun <- function(b) suppressWarnings(
esr_rho_lp(b = b, y = y, xq = xq, xe = xe, alpha = alpha, g1 = g1, g2 = g2)
)
## Optimize the model using iterated local search
# Initial optimization
fit <- try(stats::optim(par = b0, fn = fun, method = "Nelder-Mead",
control = list(maxit = 2000)), silent = TRUE)
# Counts the iterations without decrease of the loss
counter <- 0
while (counter < early_stopping) {
# Perturbe b
bt <- fit$par + stats::rnorm(kq+ke, sd=se)
# If G2 is 1, 2, 3 ensure that x'be < 0 by moving the es intercept
if (g2 %in% c(1, 2, 3)) {
max_xe <- max(xe %*% bt[(kq+1):(kq+ke)])
bt[kq+1] <- bt[kq+1] - (max_xe + 0.1) * (max_xe >= 0)
}
# Fit the model with the perturbed parameters
tmp_fit <- try(stats::optim(par = bt, fn = fun, method = "Nelder-Mead",
control = list(maxit = 2000)), silent = TRUE)
# Replace the fit if the new loss is smaller than the old. Otherwise increase the counter.
if (!inherits(tmp_fit, "try-error")) {
if (tmp_fit$value < fit$value) {
fit <- tmp_fit
counter <- 0
} else {
counter <- counter + 1
}
}
}
# Set names of the parameters
names(fit$par) <- c(paste0("bq_", 1:kq - 1), paste0("be_", 1:ke - 1))
# Undo the transformation
if (g2 %in% c(1, 2, 3)) {
y <- y + max_y
fit$par[1] <- fit$par[1] + max_y
fit$par[kq+1] <- fit$par[kq+1] + max_y
}
# Retun results
revtal <- list(
coefficients = fit$par,
coefficients_q = fit$par[1:kq],
coefficients_e = fit$par[(kq+1):(kq+ke)],
y = y,
xq = xq,
xe = xe,
alpha = alpha,
g1 = g1,
g2 = g2,
loss = fit$value,
time = Sys.time() - t0
)
class(revtal) <- 'esreg'
revtal
}
#' Lambda Matrix
#'
#' Estimate the lambda matrix.
#'
#' @inheritParams vcovA
#' @export
lambda_matrix <- function(object, sparsity, bandwidth_estimator, misspec) {
if(!(sparsity %in% c("iid", "nid")))
stop("sparsity can be iid or nid")
if(!(bandwidth_estimator %in% c("Bofinger", "Chamberlain", "Hall-Sheather")))
stop("bandwidth_estimator can be Bofinger, Chamberlain or Hall-Sheather")
# Extract elements from object
y <- object$y
xq <- object$xq
xe <- object$xe
coefficients_q <- object$coefficients_q
coefficients_e <- object$coefficients_e
alpha <- object$alpha
n <- length(y)
kq <- ncol(xq)
ke <- ncol(xe)
# Transform the data and coefficients
if (object$g2 %in% c(1, 2, 3)) {
max_y <- max(y)
y <- y - max_y
coefficients_q[1] <- coefficients_q[1] - max_y
coefficients_e[1] <- coefficients_e[1] - max_y
}
# Precompute some quantities
xbq <- as.numeric(xq %*% coefficients_q)
xbe <- as.numeric(xe %*% coefficients_e)
uq <- as.numeric(y - xbq)
# Evaluate G1 / G2 functions
G1_prime_xq <- G_vec(z = xbq, g = "G1_prime", type = object$g1)
G1_prime_prime_xq <- G_vec(z = xbq, g = "G1_prime_prime", type = object$g1)
G2_xe <- G_vec(z = xbe, g = "G2", type = object$g2)
G2_prime_xe <- G_vec(z = xbe, g = "G2_prime", type = object$g2)
G2_prime_prime_xe <- G_vec(z = xbe, g = "G2_prime_prime", type = object$g2)
# Check the methods in case of sample quantile / es
if ((kq == 1) & (ke == 1) & sparsity != "iid") {
warning("Changed sparsity estimation to iid!")
sparsity <- "iid"
}
# Density quantile function
dens <- density_quantile_function(y = y, x = xq, u = uq, alpha = object$alpha,
sparsity = sparsity, bandwidth_estimator = bandwidth_estimator)
# Conditional CDF evaluated at conditional quantile
cdf <- cdf_at_quantile(y = y, x = xq, q = xbq)
# Compute lambda
lambda <- lambda_matrix_loop(
xq = xq, xe = xe, xbq = xbq, xbe = xbe, alpha = alpha,
G1_prime_xq = G1_prime_xq, G1_prime_prime_xq = G1_prime_prime_xq,
G2_xe = G2_xe, G2_prime_xe = G2_prime_xe, G2_prime_prime_xe = G2_prime_prime_xe,
density = dens, cdf = cdf, include_misspecification_terms = misspec)
lambda
}
#' Sigma Matrix
#'
#' Estimate the sigma matrix.
#'
#' @inheritParams vcovA
#' @export
sigma_matrix <- function(object, sigma_est, misspec) {
if(!(sigma_est %in% c("ind", "scl_N", "scl_sp")))
stop("sigma_estimator can be ind, scl_N or scl_sp")
# Extract elements from object
y <- object$y
xq <- object$xq
xe <- object$xe
coefficients_q <- object$coefficients_q
coefficients_e <- object$coefficients_e
alpha <- object$alpha
n <- length(y)
kq <- ncol(xq)
ke <- ncol(xe)
# Transform the data and coefficients
if (object$g2 %in% c(1, 2, 3)) {
max_y <- max(y)
y <- y - max_y
coefficients_q[1] <- coefficients_q[1] - max_y
coefficients_e[1] <- coefficients_e[1] - max_y
}
# Precompute some quantities
xbq <- as.numeric(xq %*% coefficients_q)
xbe <- as.numeric(xe %*% coefficients_e)
uq <- as.numeric(y - xbq)
# Evaluate G1 / G2 functions
G1_prime_xq <- G_vec(z = xbq, g = "G1_prime", type = object$g1)
G2_xe <- G_vec(z = xbe, g = "G2", type = object$g2)
G2_prime_xe <- G_vec(z = xbe, g = "G2_prime", type = object$g2)
G2_prime_prime_xe <- G_vec(z = xbe, g = "G2_prime_prime", type = object$g2)
# Check the methods in case of sample quantile / es
if ((kq == 1) & (ke == 1) & sigma_est != "ind") {
warning("Changed conditional truncated variance estimation to ind!")
sigma_est <- "ind"
}
# Estimate the (conditional) truncated variance
cv <- conditional_truncated_variance(y = uq, x = xq, approach = sigma_est)
# Estimate the CDF at the quantile predictions
cdf <- cdf_at_quantile(y = y, x = xq, q = xbq)
# Compute sigma
sigma <- sigma_matrix_loop(
xq = xq, xe = xe, xbq = xbq, xbe = xbe, alpha = alpha,
G1_prime_xq = G1_prime_xq,
G2_xe = G2_xe, G2_prime_xe = G2_prime_xe,
conditional_variance = cv, cdf = cdf,
include_misspecification_terms = misspec)
sigma
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.