Nothing
#' Calculation of one-step ahead forecasts of Value at Risk and Expected Shortfall (parametric and
#' semiparametric)
#'
#'One-step ahead forecasts of Value at Risk and Expected Shortfall for a selection of short-memory
#'and long-memory parametric as well as semiparametric GARCH-type models are
#'computed.
#'
#' @importFrom stats 'dt' 'qt' 'sd' 'arima'
#' @importFrom utils 'tail'
#'
#' @export
#'
#' @param x a vector containing the price series.
#' @param a.v confidence level for calculating VaR; is set to \code{0.99} by default.
#' @param a.e confidence level for calculating ES; is set to \code{0.975} by default.
#' @param model model to be estimated. Options are 'sGARCH', 'eGARCH', 'apARCH',
#' 'lGARCH', 'fiGARCH' and 'filGARCH'; is set to \code{'sGARCH'} by default.
#' @param garchOrder orders to be estimated; c(1, 1), i.e. p = q = 1, is the
#' default.
#' @param distr distribution to use for the innovations of the respective GARCH model;
#' is set to \code{'std'} by default
#' @param n.out size of out-sample; is set to \code{250} by default.
#' @param smooth a character object; defines the data-driven smoothing approach
#' for the estimation of the nonparametric scale function; for
#' \code{smooth = 'lpr'}, the scale function is obtained from the logarithm
#' of the squared centralized returns by means of the \code{msmooth()} function or
#' \code{tsmoothlm()} function if \code{model} is set to \code{'sGARCH'}, \code{'eGARCH'},
#' \code{'apARCH'} and \code{lGARCH'} or \code{'fiGARCH} and \code{'filGARCH'},
#' respectively; is set to \code{smooth = 'none'} by default.
#' @param ... depending on the choice of \code{model}, further arguments can be
#' passed to either \code{smoots::msmooth()} or to \code{tsmoothlm()};
#' if no other arguments are given, the default settings
#' are used for both functions with the exception of \code{p = 3}.
#'
#' @details
#' Let \eqn{Y_t} be a (demeaned) return series. The semiparametric extension of
#' the GARCH(p,q) model (Bollerslev, 1986) is called a Semi-GARCH model
#' (Feng, 2004) and is defined by
#'
#' \deqn{Y_t = s(x_t)\sigma_t \eta_t,}
#' with \eqn{\eta_t \sim IID(0,1)} and
#'
#' \deqn{\sigma^2_t = \alpha_0 + \sum_{i=1}^p \alpha_i Y^2_{t-i}
#' + \sum_{j=1}^q \beta_j \sigma^2_{t-j},}
#'
#' where \eqn{\sigma_t > 0} are the conditional standard deviations, \eqn{s(x_t) > 0} is
#' a nonparametric scale function with \eqn{x_t} being the rescaled observation
#' time points on the interval [0,1] and \eqn{\alpha_i}
#' and \eqn{\beta_j} are non-negative real valued coefficients, except for
#' \eqn{\alpha_0}, which must satisfy \eqn{\alpha_0 > 0}. Furthermore, it is assumed that
#' Var\eqn{(\sigma_t \eta_t) = 1}. In this functions, different
#' short-memory and long-memory GARCH-type models are selectable for the parametric part of the
#' model. More specifically, the standard GARCH (Bollerslev, 1986), the
#' Log-GARCH (Pantula, 1986; Geweke, 1986; Milhoj, 1988), the eGARCH
#' (Nelson, 1991), the APARCH (Ding et al., 1993), the FIGARCH (Baillie et al., 1996)
#' and the FI-Log-GARCH (Feng et al., 2020) model are implemented. For more
#' information on the representations of the last three models mentioned, we
#' refer the reader to the corresponding references listed in the references
#' section.
#'
#' While the innovations \eqn{\eta_t} must be i.i.d. (independent and identically
#' distributed) with zero-mean and unit-variance and while any
#' distribution that satisfies these conditions is suitable, the standardized
#' t-distribution is selected for the estimation of the models and computation
#' of the Value at Risk (VaR) as well as the Expected Shortfall (ES) within
#' this function.
#'
#' For a given level \eqn{\alpha \in (0, 1)},
#'
#' \deqn{VaR(\alpha) = inf \{z \in R: F_L(z) \geq \alpha\}}
#'
#' defines the VaR at level alpha. In this definition, \eqn{L} is the
#' loss variable (making a loss is denoted as a positive value, whereas gains
#' are negative values) and \eqn{F_L} is its cumulative distribution function.
#' Explained differently, \eqn{VaR(\alpha)} is the \eqn{\alpha}-quantile of the loss
#' distribution.
#'
#' The ES for a level \eqn{\alpha}, however, is given by
#'
#' \deqn{ES(\alpha) = (1 / (1 - \alpha)) \int_{\alpha}^1 VaR(u)du,}
#'
#' i.e. it is the expected loss in case \eqn{VaR(\alpha)} is exceeded. More
#' information on these risk measures can be found on pp. 64-72 in McNeil et
#' al. (2015).
#'
#' To apply the function, a numeric vector that contains the price series that
#' is to be analyzed ordered from past to present must be passed to the
#' argument \code{x}. Furthermore, the user can set different levels of alpha
#' for the VaR and the ES via the arguments \code{a.v} and \code{a.e},
#' respectively. A parametric short-memory or long-memory GARCH-type model can be
#' selected by means of \code{model}, which only accepts a single-element character vector
#' as input. At the time of the release of package version 1.0.0, a standard
#' GARCH ('sGARCH'), a Log-GARCH ('lGARCH'), an eGARCH ('eGARCH'), an APARCH
#' ('apARCH'), a FIGARCH ('fiGARCH') and a FI-Log-GARCH ('filGARCH') model can be selected,
#' each with conditional t-distribution. By default, a standard GARCH model is applied.
#' The orders of the GARCH-type models can be defined with \code{garchOrder},
#' which is a numeric vector with two elements. Its first element is the
#' ARCH order p, whereas the GARCH order q can be adjusted via the second
#' element. If no adjustments are made, the orders p = q = 1 are selected. The
#' number of out-sample observations is set via the argument \code{n.out}. If n
#' is the total number of observations of the
#' whole price series, the model is estimated for the first n - n.out
#' observations (in-sample), while the VaR and the ES are obtained for the last
#' n.out observations (out-sample) based on the estimated model for the
#' in-sample. Moreover, the data-driven estimation method of the underlying
#' scale function can be adjusted via the argument \code{smooth}. If
#' \code{smooth = 'lpr'} is selected, the scale function is obtained by
#' applying an iterative plug-in algorithm logarithm of the squared
#' centralized returns. Depending on the setting of \code{model} an algorithm
#' proposed by Feng, Gries and Fritz (2020) or by Letmathe, Feng and Uhde
#' (2021) is employed. In the former case, the function \code{msmooth()} of the
#' \code{smoots} package is applied and for the latter the \code{tsmoothlm()}
#' function of the \code{esemifar} package is used. An ellipsis \code{...} is
#' implemented to allow for additional arguments for \code{msmooth()} and
#' \code{tsmoothlm()}.
#'
#' NOTE:
#'
#' This function makes use of the \code{arima()} function of the stats package,
#' of the \code{fracdiff()} function of the \code{fracdiff} package, of the
#' \code{ugarchspec()} and \code{ugarchfit()} functions of the \code{rugarch}
#' package, of the \code{msmooth()} function of the \code{smoots} package
#' and of the \code{esemifar()} function of the \code{esemifar} for estimation.
#' Moreover, Log-GARCH and FI-Log-GARCH models in the parametric part of the
#' complete models are estimated via their ARMA and FARIMA representations,
#' respectively, and must therefore satisfy \eqn{p \geq q}.
#'
#' @return
#' This function returns a list with the following elements.
#' \describe{
#' \item{model}{selected model for estimation}
#' \item{mean}{the estimated mean of the in-sample returns}
#' \item{model.fit}{estimated model parameters for the parametric part of the
#' in-sample}
#' \item{np.est}{the estimation results for the nonparametric part of the
#' in-sample model}
#' \item{ret.in}{in-sample return series}
#' \item{ret.out}{out-sample return series}
#' \item{sig.in}{estimated in-sample total volatility}
#' \item{sig.fc}{out-sample forecasts of the total volatility}
#' \item{scale}{the estimated nonparametric scale function values for the
#' in-sample}
#' \item{scale.fc}{the scale function forecast for the out-sample}
#' \item{VaR.e}{out-sample forecasts of the (1-\code{a.e})100\% VaR}
#' \item{VaR.v}{out-sample forecasts of the (1-\code{a.v})100\% VaR}
#' \item{ES}{out-sample forecasts of the (1-\code{a.e})100\% ES}
#' \item{dfree}{estimated degrees of freedom for the standardized returns}
#' \item{a.v}{coverage level for the 99 \% VaR}
#' \item{a.e}{coverage level for 97.5 \% VaR}
#' \item{garchOrder}{the orders p and q of the implemented GARCH-type model}
#' }
#'
#' @author
#'\itemize{
#'\item Sebastian Letmathe (Scientific Employee) (Department of Economics,
#'Paderborn University) \cr
#'\item Dominik Schulz (Scientific Employee) (Department of Economics,
#'Paderborn University), \cr
#'}
#'
#' @references
#' Baillie, R. T., Bollerslev, T., & Mikkelsen, H. O. (1996). Fractionally
#' integrated generalized autoregressive conditional heteroskedasticity.
#' In: Journal of Econometrics, 74.1, pp. 3-30.
#'
#' Bollerslev, T. (1986) Generalized autoregressive conditional
#' heteroskedasticity. In: Journal of Econometrics 31.3, pp. 307-327.
#'
#' Ding, Z., Granger, C.W., and Engle, R.F. (1993). A long memory property
#' of stock market returns and a new model. In: Journal of Empirical Finance
#' 1.1, pp. 83-106.
#'
#' Feng, Y. (2004). Simultaneously modeling conditional heteroskedasticity and
#' scale change. In: Econometric Theory 20.3, pp. 563-596.
#'
#' Feng, Y., Beran, J., Letmathe, S., & Ghosh, S. (2020). Fractionally
#' integrated Log-GARCH with application to value at risk and expected
#' shortfall (No. 137). Paderborn University, CIE Center for
#' International Economics.
#'
#' Pantula, S.G. (1986). Modeling the persistence of conditional variances:
#' a comment. In: Econometric Reviews 5, pp. 79-97.
#'
#' Geweke, J. (1986). Comment on: Modelling the persistence of conditional
#' variances. In: Econometric Reviews 5, pp. 57-61.
#'
#' Letmathe, S., Feng, Y., & Uhde, A. (2021). Semiparametric GARCH models with
#' long memory applied to Value at Risk and Expected Shortfall (No. 141).
#' Paderborn University, CIE Center for International Economics.
#'
#' McNeil, A.J., Frey, R., and Embrechts, P. (2015). Quantitative risk
#' management: concepts, techniques and tools - revised edition. Princeton
#' University Press.
#'
#' Milhoj, A. (1988). A Multiplicative parametrization of ARCH models.
#' Universitetets Statistiske Institut.
#'
#' Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A
#' new approach. In: Econometrica: Journal of the Econometric Society, 347-370.
#'
#'
#' @examples
#'
#'# Example for Walmart Inc. (WMT)
#' prices = WMT$price.close
#'
#' # forecasting VaR and ES
#' results = varcast(prices, model = 'sGARCH', n.out = 250)
#' ret.out = results$ret.out
#' n.out = length(ret.out)
#' VaR97.5 = results$VaR.e
#' VaR99 = results$VaR.v
#' ES = results$ES
#'
#'# plotting VaR at 99% coverage
#' matplot(1:n.out, cbind(-ret.out, VaR99),
#' type = 'hl',
#' xlab = 'number of out-of-sample obs.', ylab = 'losses, VaR and ES',
#' main = '99% VaR (red) for the WMT return series')
#'
#'# plotting VaR at 97.5% coverage and corresponding ES
#' matplot(1:n.out, cbind(-ret.out, ES, VaR97.5),
#' type = 'hll',
#' xlab = 'number of out-of-sample obs.', ylab = 'losses, VaR and ES',
#' main = '97.5% VaR (green) and ES (red) for the WMT return series')
#'
varcast <- function(x, a.v = 0.99, a.e = 0.975,
model = c("sGARCH", "lGARCH", "eGARCH", "apARCH", "fiGARCH",
"filGARCH"),
garchOrder = c(1, 1),
distr = c("norm", "std"),
n.out = 250,
smooth = c("none", "lpr"), ...) {
if (length(x) <= 1 || !all(!is.na(x)) || !is.numeric(x)) {
stop("A numeric vector of length > 1 and without NAs must be passed to",
" 'x'.")
}
if (length(a.v) != 1 || is.na(a.v) || !is.numeric(a.v) || a.v <= 0 || a.v >= 1) {
stop("A single numeric value that satisfies >0 and <1 must be passed to",
" 'a.v'")
}
if (length(a.e) != 1 || is.na(a.e) || !is.numeric(a.e) || a.e <= 0 || a.e >= 1) {
stop("A single numeric value that satisfies >0 and <1 must be passed to",
" 'a.e'")
}
if (!(length(model) %in% c(1, 6)) || !all(!is.na(model)) || !is.character(model)) {
stop("A single character value must be passed to 'model'.")
}
if (all(model == c("sGARCH", "lGARCH", "eGARCH", "apARCH", "fiGARCH", "filGARCH"))) {
model <- "sGARCH"
}
if (length(garchOrder) != 2 || !all(!is.na(garchOrder)) || !is.numeric(garchOrder) ||
(garchOrder[[1]] == 0 && garchOrder[[2]] == 0)) {
stop("A vector of length 2 must be passed to 'garchOrder' giving the",
" first and second order of the GARCH-type model.")
}
garchOrder <- floor(garchOrder)
if (length(n.out) != 1 || is.na(n.out) || !is.numeric(n.out) || n.out <= 0) {
stop("A single positive integer must be passed to 'n.out'.")
}
n.out <- floor(n.out)
if (!(length(distr) %in% c(1, 2)) || !all(!is.na(distr)) ||
!is.character(distr)) {
stop("Input argument 'distr' must be a single character ",
"value. Valid choices are 'norm' and 'std'.")
}
if (all(distr == c("norm", "std"))) distr <- "std"
if (length(distr) != 1 || !(distr) %in% c("norm", "std")) {
stop("Input argument 'distr' must be a single character ",
"value. Valid choices are 'norm' and 'std'.")
}
if (!(length(smooth) %in% c(1, 2)) || !all(!is.na(smooth)) ||
!is.character(smooth)) {
stop("Input argument 'smooth' must be a single character ",
"value. Valid choices are 'none' and 'lpr'.")
}
if (all(smooth == c("none", "lpr"))) smooth <- "none"
if (length(smooth) != 1 || !(smooth) %in% c("none", "lpr")) {
stop("Input argument 'smooth' must be a single character ",
"value. Valid choices are 'none' and 'lpr'.")
}
if (smooth == "lpr" && model %in% c("sGARCH", "lGARCH", "eGARCH", "apARCH")) {
smooth.method = "smoots"
}
if (smooth == "lpr" && model %in% c("fiGARCH", "filGARCH")) {
smooth.method = "esemifar"
}
dots <- list(...)
ret <- diff(log(x))
n.ret <- length(ret)
n.in <- n.ret - n.out
ret.in.t <- ret[1:n.in]
mean.ret.in <- mean(ret.in.t)
ret.in <- ret.in.t - mean.ret.in
ret.out.t <- ret[(n.in + 1):n.ret]
ret.out <- ret.out.t - mean.ret.in
p.true <- p <- garchOrder[1]
q.true <- q <- garchOrder[2]
if (p.true == 0) p <- 1
if (q.true == 0) q <- 1
l <- max(p, q)
if (is.null(dots[["p"]])) {
dots[["p"]] <- 3
}
dots[["y"]] <- log(ret.in^2)
switch(smooth,
none = {
np.est <- NA
zeta.in <- ret.in
zeta.out <- ret.out
res.in <- log(ret.in^2)
res.out <- log(ret.out^2)
sxt <- 1
sfc <- 1
mule <- mean(res.in)
Csig <- 1
},
lpr = {
switch(smooth.method,
smoots = {
if (is.null(dots[["alg"]])) {
dots[["alg"]] <- "A"
}
np.est <- suppressMessages(do.call(what = smoots::msmooth,
args = dots))
res.in <- np.est[["res"]]
res.out <- log(ret.out^2) - np.est[["ye"]][n.in]
mule <- -log(mean(exp(res.in)))
sxt <- exp(0.5 * (np.est[["ye"]] - mule))
sfc <- sxt[n.in]
zeta.in <- ret.in / sxt
zeta.out <- ret.out / sfc
},
esemifar = {
np.est <- suppressMessages(do.call(what = esemifar::tsmoothlm,
args = dots))
res.in <- np.est[["res"]]
res.out <- log(ret.out^2) - np.est[["ye"]][n.in]
Csig <- sd(ret.in / exp(0.5 * np.est[["ye"]]))
sxt <- exp(0.5 * np.est[["ye"]]) * Csig
sfc <- sxt[n.in]
zeta.in <- ret.in / sxt
zeta.out <- ret.out / sfc
})
}
)
zeta.n <- tail(zeta.in, n = l)
zeta.fc <- c(zeta.n, zeta.out[1:n.out])
if (model %in% c("sGARCH", "eGARCH", "apARCH", "fiGARCH")) {
spec <- rugarch::ugarchspec(variance.model =
list(model = model,
garchOrder = c(p.true, q.true)),
mean.model = list(armaOrder = c(0, 0),
include.mean = FALSE),
distribution.model = distr)
model_fit <- rugarch::ugarchfit(spec = spec, data = c(zeta.in, zeta.out),
out.sample = n.out,
fit.control = list(trunclag = min(n.in, 1000)))
if (distr == "std") dfree <- model_fit@fit$coef[["shape"]]
if (distr == "norm") dfree <- NA
c.sig <- as.numeric(rugarch::sigma(model_fit))
sig.in <- c.sig * sxt
sig.fc <- rugarch::ugarchforecast(model_fit, n.ahead = 1, n.roll = n.out - 1)
sig.fc <- as.numeric(rugarch::sigma(sig.fc)) * sfc
}
switch(model,
lGARCH = {
if (p.true < q.true) {
stop("p >= q must be satisfied for the estimation of a Log-GARCH
model ", "via its ARMA representation.")
}
model_fit <- arima(res.in,
order = c(p.true, 0, q.true),
include.mean = FALSE)
y.cent <- c(res.in, res.out) - mule * isTRUE(smooth == "none")
mulz <- -log(mean(exp(model_fit[["residuals"]])))
ar <- model_fit$model$phi
ma <- -model_fit$model$theta
k <- 50
d <- 0
coef.all <- arfilt(ar, ma, d, k)
pre0 <- c()
add0 <- 0
if (n.in < k) {
add0 <- k - n.in
pre0 <- rep(0, add0)
}
y.fc.out <- (1:n.out) * 0
y.cent <- c(pre0, y.cent)
for (i in (add0 + n.in + 1):(add0 + n.ret)) {
y.fc.out[i - (add0 + n.in)] <- coef.all %*%
y.cent[(i - 1):(i - k)]
}
sig.fc <- exp(0.5 * (y.fc.out + mule - mulz))
sig.fc.in <- exp(0.5 * (y.cent[1:n.in] - model_fit[["residuals"]]
+ mule - mulz))
sig.fc <- sig.fc * sfc
sig.in <- sig.fc.in * sxt
ret.sd <- zeta.in / sig.in
if (distr == "std") {
dfree <- as.numeric(rugarch::fitdist("std", ret.sd)$pars[[3]])
}
if (distr == "norm") {
dfree <- NA
}
},
filGARCH = {
if (p.true < q.true) {
stop("p >= q must be satisfied for the estimation of a Log-GARCH
model ", "via its ARMA representation.")
}
model_fit <- suppressWarnings(fracdiff::fracdiff(res.in,
nar = p.true,
nma = q.true))
ar <- model_fit[["ar"]]
ma <- model_fit[["ma"]]
d <- model_fit[["d"]]
k <- 50
coef.all <- arfilt(ar, ma, d, k)
pre0 <- c()
add0 <- 0
y.cent <- c(res.in, res.out)
y.cent <- y.cent - mean(res.in)
if (n.in < k) {
add0 <- k - n.in
pre0 <- rep(0, add0)
}
y.fc.out <- (1:n.out) * 0
y.cent <- c(pre0, y.cent)
for (i in (add0 + n.in + 1):(add0 + n.ret)) {
y.fc.out[i - (add0 + n.in)] <- coef.all %*% y.cent[(i - 1):(i - k)]
}
sig.fc <- exp(0.5 * (y.fc.out + mean(res.in)))
y.fc.in <- model_fit[["fitted"]]
sig.fc.in <- exp(0.5 * y.fc.in)
ret.sd <- zeta.in / (sig.fc.in * sxt)
if (distr == "std") {
dfree <- as.numeric(rugarch::fitdist("std", ret.sd)$pars[3])
}
if (distr == "norm") {
dfree <- NA
}
if (smooth == "lpr") {
Ch <- sd(zeta.in / sig.fc.in)
}
if (smooth == "none") {
Ch <- sd(ret.sd)
}
sig.fc <- sig.fc * sfc * Ch
sig.in <- sig.fc.in * sxt * Ch
}
)
switch(distr,
norm = {
sdev <- 1
ES0t <- 1
arglst <- list(p = c(a.e, a.v), distr = distr)
},
std = {
sdev <- sqrt(dfree / (dfree - 2))
ES0t <- (dfree + quant(a.e, df = dfree, distr = distr)^2) / (dfree - 1)
arglst <- list(p = c(a.e, a.v), df = dfree, distr = distr)
}
)
quants <- do.call(what = quant, args = arglst)
q.e <- quants[1]
q.a <- quants[2]
VaR.fc <- -mean.ret.in + sig.fc * q.a / sdev
VaR.e.fc <- -mean.ret.in + sig.fc * q.e / sdev
arglst[[1]] <- q.e; names(arglst)[1] <- "x"
d.es <- do.call(what = dens, args = arglst)
ES0 <- d.es / (1 - a.e) * ES0t / sdev
Esfc <- -mean.ret.in + sig.fc * ES0
if (smooth == "lpr") {
model <- paste0("Semi-", model)
}
results <- list(model = model, mean = mean.ret.in,
model.fit = model_fit, ret.in = ret.in.t,
ret.out = ret.out.t, sig.in = sig.in,
sig.fc = sig.fc, scale = sxt, scale.fc = sfc, VaR.e = VaR.e.fc,
VaR.v = VaR.fc, ES = Esfc, dfree = dfree, a.v = 1 - a.v,
a.e = 1 - a.e, garchOrder = garchOrder, np.est = np.est)
class(results) <- "ufRisk"
attr(results, "function") <- "varcast"
results
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.