# Fitting and Predicting VaR based on an ARMA-GARCH Process

#### 2024-03-04

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)

## 1 Simulate (-log-return) data $$(X_t)$$ from an ARMA-GARCH process

We consider an ARMA(1,1)-GARCH(1,1) process with $$t$$ distributed innovations.

## Model specification (for simulation)
nu <- 3 # d.o.f. of the standardized distribution of Z_t
fixed.p <- list(mu = 0, # our mu (intercept)
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
armaOrder <- c(1,1) # ARMA order
garchOrder <- c(1,1) # GARCH order
varModel <- list(model = "sGARCH", garchOrder = garchOrder)
spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
fixed.pars = fixed.p, distribution.model = "std") # t standardized residuals

Simulate one path (for illustration purposes).

## Simulate (X_t)
n <- 1000 # sample size (= length of simulated paths)
x <- ugarchpath(spec, n.sim = n, m.sim = 1, rseed = 271) # n.sim length of simulated path; m.sim = number of paths
## Note the difference:
## - ugarchpath(): simulate from a specified model
## - ugarchsim():  simulate from a fitted object

## Extract the resulting series
X <- fitted(x) # simulated process X_t = mu_t + epsilon_t for epsilon_t = sigma_t * Z_t
sig <- sigma(x) # volatilities sigma_t (conditional standard deviations)
eps <- x@path$residSim # unstandardized residuals epsilon_t = sigma_t * Z_t ## 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])) ## 2 Fit an ARMA-GARCH model to the (simulated) data 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 spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder), distribution.model = "std") # without fixed parameters here fit <- ugarchfit(spec, data = X) # fit ## Extract the resulting series mu. <- fitted(fit) # fitted hat{mu}_t (= hat{X}_t) sig. <- sigma(fit) # fitted hat{sigma}_t ## 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 resi <- as.numeric(residuals(fit)) 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)
Z <- as.numeric(residuals(fit, standardize = TRUE))
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)))) ## 3 Calculate the VaR time series Compute VaR estimates. Note that we could have also used the GPD-based estimators here. ## VaR confidence level we consider here alpha <- 0.99 ## Extract fitted VaR_alpha VaR. <- as.numeric(quantile(fit, probs = alpha)) ## Build manually and compare the two nu. <- fit@fit$coef[["shape"]] # extract (fitted) d.o.f. nu
VaR.. <- as.numeric(mu. + sig. * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.)) # VaR_alpha computed manually
stopifnot(all.equal(VaR.., VaR.))
## => quantile(<rugarch object>, probs = alpha) provides VaR_alpha = hat{mu}_t + hat{sigma}_t * q_Z(alpha)

## 4 Backtest VaR estimates

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)
btest <- VaRTest(1-alpha, actual = -X,
VaR = quantile(ugarchfit(spec, data = -X), probs = 1-alpha))
btest$expected.exceed # number of expected exceedances = (1-alpha) * n ## [1] 10 btest$actual.exceed # actual exceedances
## [1] 12
## Unconditional test
btest$uc.H0 # corresponding null hypothesis ## [1] "Correct Exceedances" btest$uc.Decision # test decision
## [1] "Fail to Reject H0"
## Conditional test
btest$cc.H0 # corresponding null hypothesis ## [1] "Correct Exceedances & Independent" btest$cc.Decision # test decision
## [1] "Fail to Reject H0"

## 5 Predict VaR based on fitted model

Now predict VaR.

## Predict from the fitted process
fspec <- getspec(fit) # specification of the fitted process
setfixed(fspec) <- as.list(coef(fit)) # set the parameters to the fitted ones
m <- ceiling(n / 10) # number of steps to forecast (roll/iterate m-1 times forward with frequency 1)
pred <- ugarchforecast(fspec, data = X, n.ahead = 1, n.roll = m-1, out.sample = m-1) # predict from the fitted process

## Extract the resulting series
mu.predict <- fitted(pred) # extract predicted X_t (= conditional mean mu_t; note: E[Z] = 0)
sig.predict <- sigma(pred) # extract predicted sigma_t
VaR.predict <- as.numeric(quantile(pred, probs = alpha)) # corresponding predicted VaR_alpha

## 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
VaR.predict. <- as.numeric(mu.predict + sig.predict * sqrt((nu.-2)/nu.) *
qt(alpha, df = nu.)) # VaR_alpha computed manually
stopifnot(all.equal(VaR.predict., VaR.predict))

## 6 Simulate future trajectories of $$(X_t)$$ and compute corresponding VaRs

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
B <- 1000
set.seed(271)
X.sim.obj <- ugarchpath(fspec, n.sim = m, m.sim = B) # simulate future paths

## 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)
X.sim <- fitted(X.sim.obj) # extract simulated X_t
sig.sim <- sigma(X.sim.obj) # extract sigma_t
eps.sim <- X.sim.obj@path$residSim # extract epsilon_t VaR.sim <- (X.sim - eps.sim) + sig.sim * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.) # (m, B) matrix VaR.CI <- apply(VaR.sim, 1, function(x) quantile(x, probs = c(0.025, 0.975))) ## 7 Plot Finally, letâ€™s display all results. ## Setup yran <- range(X, # simulated path mu., VaR., # fitted conditional mean and VaR_alpha mu.predict, VaR.predict, VaR.CI) # predicted mean, VaR and CIs myran <- max(abs(yran)) yran <- c(-myran, myran) # y-range for the plot xran <- c(1, length(X) + m) # x-range for the plot ## 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
t. <- length(X) + seq_len(m) # future time points
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))))