This vignette does not use qrmtools, but shows how Value-at-Risk (VaR) can be fitted and predicted based on an underlying ARMA-GARCH process (which of course also concerns QRM in the wider sense).
library(rugarch)
library(qrmtools)
We consider an ARMA(1,1)-GARCH(1,1) process with \(t\) distributed innovations.
## Model specification (for simulation)
3 # d.o.f. of the standardized distribution of Z_t
nu <- list(mu = 0, # our mu (intercept)
fixed.p <-ar1 = 0.5, # our phi_1 (AR(1) parameter of mu_t)
ma1 = 0.3, # our theta_1 (MA(1) parameter of mu_t)
omega = 4, # our alpha_0 (intercept)
alpha1 = 0.4, # our alpha_1 (GARCH(1) parameter of sigma_t^2)
beta1 = 0.2, # our beta_1 (GARCH(1) parameter of sigma_t^2)
shape = nu) # d.o.f. nu for standardized t_nu innovations
c(1,1) # ARMA order
armaOrder <- c(1,1) # GARCH order
garchOrder <- list(model = "sGARCH", garchOrder = garchOrder)
varModel <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
spec <-fixed.pars = fixed.p, distribution.model = "std") # t standardized residuals
Simulate one path (for illustration purposes).
## Simulate (X_t)
1000 # sample size (= length of simulated paths)
n <- ugarchpath(spec, n.sim = n, m.sim = 1, rseed = 271) # n.sim length of simulated path; m.sim = number of paths
x <-## Note the difference:
## - ugarchpath(): simulate from a specified model
## - ugarchsim(): simulate from a fitted object
## Extract the resulting series
fitted(x) # simulated process X_t = mu_t + epsilon_t for epsilon_t = sigma_t * Z_t
X <- sigma(x) # volatilities sigma_t (conditional standard deviations)
sig <- x@path$residSim # unstandardized residuals epsilon_t = sigma_t * Z_t
eps <-## Note: There are no extraction methods for the unstandardized residuals epsilon_t
## for uGARCHpath objects (only for uGARCHfit objects; see below).
## Sanity checks (=> fitted() and sigma() grab out the right quantities)
stopifnot(all.equal(X, x@path$seriesSim, check.attributes = FALSE),
all.equal(sig, x@path$sigmaSim, check.attributes = FALSE))
As a sanity check, let’s plot the simulated path, conditional standard deviations and residuals.
## Plots
plot(X, type = "l", xlab = "t", ylab = expression(X[t]))
plot(sig, type = "h", xlab = "t", ylab = expression(sigma[t]))
plot(eps, type = "l", xlab = "t", ylab = expression(epsilon[t]))
Fit an ARMA-GARCH process to X
(with the correct, known orders here; one would normally fit processes of different orders and then decide).
## Fit an ARMA(1,1)-GARCH(1,1) model
ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
spec <-distribution.model = "std") # without fixed parameters here
ugarchfit(spec, data = X) # fit
fit <-
## Extract the resulting series
fitted(fit) # fitted hat{mu}_t (= hat{X}_t)
mu. <- sigma(fit) # fitted hat{sigma}_t
sig. <-
## Sanity checks (=> fitted() and sigma() grab out the right quantities)
stopifnot(all.equal(as.numeric(mu.), fit@fit$fitted.values),
all.equal(as.numeric(sig.), fit@fit$sigma))
Again let’s consider some sanity checks.
## Plot data X_t and fitted hat{mu}_t
plot(X, type = "l", xlab = "t",
ylab = expression("Data"~X[t]~"and fitted values"~hat(mu)[t]))
lines(as.numeric(mu.), col = adjustcolor("blue", alpha.f = 0.5))
legend("bottomright", bty = "n", lty = c(1,1),
col = c("black", adjustcolor("blue", alpha.f = 0.5)),
legend = c(expression(X[t]), expression(hat(mu)[t])))
## Plot the unstandardized residuals epsilon_t
as.numeric(residuals(fit))
resi <-stopifnot(all.equal(fit@fit$residuals, resi))
plot(resi, type = "l", xlab = "t", ylab = expression(epsilon[t])) # check residuals epsilon_t
## Q-Q plot of the standardized residuals Z_t against their specified t
## (t_nu with variance 1)
as.numeric(residuals(fit, standardize = TRUE))
Z <-stopifnot(all.equal(Z, fit@fit$z, check.attributes = FALSE),
all.equal(Z, as.numeric(resi/sig.)))
qq_plot(Z, FUN = function(p) sqrt((nu-2)/nu) * qt(p, df = nu),
main = substitute("Q-Q plot of ("*Z[t]*") against a standardized"~italic(t)[nu.],
list(nu. = round(nu, 2))))
Compute VaR estimates. Note that we could have also used the GPD-based estimators here.
## VaR confidence level we consider here
0.99
alpha <-
## Extract fitted VaR_alpha
as.numeric(quantile(fit, probs = alpha))
VaR. <-
## Build manually and compare the two
fit@fit$coef[["shape"]] # extract (fitted) d.o.f. nu
nu. <- as.numeric(mu. + sig. * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.)) # VaR_alpha computed manually
VaR.. <-stopifnot(all.equal(VaR.., VaR.))
## => quantile(<rugarch object>, probs = alpha) provides VaR_alpha = hat{mu}_t + hat{sigma}_t * q_Z(alpha)
Let’s backtest the VaR estimates.
## Note: VaRTest() is written for the lower tail (not sign-adjusted losses)
## (hence the complicated call here, requiring to refit the process to -X)
VaRTest(1-alpha, actual = -X,
btest <-VaR = quantile(ugarchfit(spec, data = -X), probs = 1-alpha))
$expected.exceed # number of expected exceedances = (1-alpha) * n btest
## [1] 10
$actual.exceed # actual exceedances btest
## [1] 12
## Unconditional test
$uc.H0 # corresponding null hypothesis btest
## [1] "Correct Exceedances"
$uc.Decision # test decision btest
## [1] "Fail to Reject H0"
## Conditional test
$cc.H0 # corresponding null hypothesis btest
## [1] "Correct Exceedances & Independent"
$cc.Decision # test decision btest
## [1] "Fail to Reject H0"
Now predict VaR.
## Predict from the fitted process
getspec(fit) # specification of the fitted process
fspec <-setfixed(fspec) <- as.list(coef(fit)) # set the parameters to the fitted ones
ceiling(n / 10) # number of steps to forecast (roll/iterate m-1 times forward with frequency 1)
m <- ugarchforecast(fspec, data = X, n.ahead = 1, n.roll = m-1, out.sample = m-1) # predict from the fitted process
pred <-
## Extract the resulting series
fitted(pred) # extract predicted X_t (= conditional mean mu_t; note: E[Z] = 0)
mu.predict <- sigma(pred) # extract predicted sigma_t
sig.predict <- as.numeric(quantile(pred, probs = alpha)) # corresponding predicted VaR_alpha
VaR.predict <-
## Checks
## Sanity checks
stopifnot(all.equal(mu.predict, pred@forecast$seriesFor, check.attributes = FALSE),
all.equal(sig.predict, pred@forecast$sigmaFor, check.attributes = FALSE)) # sanity check
## Build predicted VaR_alpha manually and compare the two
as.numeric(mu.predict + sig.predict * sqrt((nu.-2)/nu.) *
VaR.predict. <- qt(alpha, df = nu.)) # VaR_alpha computed manually
stopifnot(all.equal(VaR.predict., VaR.predict))
Simulate paths, estimate VaR for each simulated path (note that quantile()
can’t be used here so we have to construct VaR manually) and compute bootstrapped confidence intervals for \(\mathrm{VaR}_\alpha\).
## Simulate B paths
1000
B <-set.seed(271)
ugarchpath(fspec, n.sim = m, m.sim = B) # simulate future paths
X.sim.obj <-
## Compute simulated VaR_alpha and corresponding (simulated) confidence intervals
## Note: Each series is now an (m, B) matrix (each column is one path of length m)
fitted(X.sim.obj) # extract simulated X_t
X.sim <- sigma(X.sim.obj) # extract sigma_t
sig.sim <- X.sim.obj@path$residSim # extract epsilon_t
eps.sim <- (X.sim - eps.sim) + sig.sim * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.) # (m, B) matrix
VaR.sim <- apply(VaR.sim, 1, function(x) quantile(x, probs = c(0.025, 0.975))) VaR.CI <-
Finally, let’s display all results.
## Setup
range(X, # simulated path
yran <-# fitted conditional mean and VaR_alpha
mu., VaR., # predicted mean, VaR and CIs
mu.predict, VaR.predict, VaR.CI) max(abs(yran))
myran <- c(-myran, myran) # y-range for the plot
yran <- c(1, length(X) + m) # x-range for the plot
xran <-
## Simulated (original) data (X_t), fitted conditional mean mu_t and VaR_alpha
plot(X, type = "l", xlim = xran, ylim = yran, xlab = "Time t", ylab = "",
main = "Simulated ARMA-GARCH, fit, VaR, VaR predictions and CIs")
lines(as.numeric(mu.), col = adjustcolor("darkblue", alpha.f = 0.5)) # hat{\mu}_t
lines(VaR., col = "darkred") # estimated VaR_alpha
mtext(paste0("Expected exceed.: ",btest$expected.exceed," ",
"Actual exceed.: ",btest$actual.exceed," ",
"Test: ", btest$cc.Decision),
side = 4, adj = 0, line = 0.5, cex = 0.9) # label
## Predictions
length(X) + seq_len(m) # future time points
t. <-lines(t., mu.predict, col = "blue") # predicted process X_t (or mu_t)
lines(t., VaR.predict, col = "red") # predicted VaR_alpha
lines(t., VaR.CI[1,], col = "orange") # lower 95%-CI for VaR_alpha
lines(t., VaR.CI[2,], col = "orange") # upper 95%-CI for VaR_alpha
legend("bottomright", bty = "n", lty = rep(1, 6), lwd = 1.6,
col = c("black", adjustcolor("darkblue", alpha.f = 0.5), "blue",
"darkred", "red", "orange"),
legend = c(expression(X[t]), expression(hat(mu)[t]),
expression("Predicted"~mu[t]~"(or"~X[t]*")"),
substitute(widehat(VaR)[a], list(a = alpha)),
substitute("Predicted"~VaR[a], list(a = alpha)),
substitute("95%-CI for"~VaR[a], list(a = alpha))))