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

#### 2018-01-25

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(qrmtools) # for qq_plot()
library(rugarch)

## 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.

## Setup to simulate from
nu <- 3 # d.o.f. of the standardized distribution of Z_t
fixed.p <- list(mu = 0, ar1 = 0.5, ma1 = 0.3, omega = 4, alpha1 = 0.4, beta1 = 0.2,
shape = nu)
armaOrder <- c(1,1)
garchOrder <- c(1,1)
varModel <- list(model = "sGARCH", garchOrder = garchOrder)
spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder),
fixed.pars = fixed.p, distribution.model = "std")

Simulate one path (for illustration purposes).

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

## 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) # volatility sigma_t (conditional standard deviations)
eps <- x@path$residSim # unstandardized residuals epsilon_t = sigma_t * Z_t ## Sanity checks (=> fitted() and sigma() grab out the right quantities) ## Note: There is no such function for epsilon_t 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("Simulated process"~X[t])) plot(sig, type = "l", xlab = "t", ylab = expression("Conditional standard deviation"~sigma[t])) plot(eps, type = "l", xlab = "t", ylab = expression("Residuals"~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). ## Fit an ARMA(1,1)-GARCH(1,1) model varModel <- list(model = "sGARCH", garchOrder = garchOrder) spec <- ugarchspec(varModel, mean.model = list(armaOrder = armaOrder), distribution.model = "std") fit <- ugarchfit(spec, data = X) # components fit, model ## 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. ## Data X_t vs 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 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

## Plot the standardized residuals Z_t
Z <- fit@fit$z stopifnot(all.equal(Z, as.numeric(resi/sig))) qq_plot(Z, FUN = function(p) sqrt((nu-2)/nu) * qt(p, df = nu)) ## 3 Calculate the VaR time series Compute VaR estimates. Note that we could have also used the GPD-based estimators here. ## VaR estimates and check alpha <- 0.99 # VaR confidence level we consider here VaR <- as.numeric(quantile(fit, probs = alpha)) # a vector (since fit is a rugarch object) 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 Backtesting via randomness check

Letâ€™s backtest the VaR estimates.

## Backtest VaR_0.99
btest <- VaRTest(alpha, actual = X, VaR = VaR, conf.level = 0.95)
btest$expected.exceed # 0.99 * n ## [1] 990 btest$actual.exceed
## [1] 988
all.equal(sig.predict, pred@forecast$sigmaFor, check.attributes = FALSE)) # sanity check ## Check VaR nu. <- pred@model$fixed.pars$shape # extract (fitted) d.o.f. nu VaR. <- as.numeric(mu.predict + sig.predict * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.)) # VaR_alpha computed manually stopifnot(all.equal(VaR., 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 X.boot <- ugarchpath(fspec, n.sim = m, m.sim = B, rseed = 271) # simulate future paths; components path, model, seed ## Bootstrap VaR ## Note: Each series is now an (m, B) matrix (each column is one path of length m) X.t.boot <- fitted(X.boot) # extract simulated X_t sig.t.boot <- sigma(X.boot) # extract sigma_t eps.t.boot <- X.boot@path$residSim # extract epsilon_t
VaR.boot <- (X.t.boot - eps.t.boot) + sig.t.boot * sqrt((nu.-2)/nu.) * qt(alpha, df = nu.) # (m, B) matrix
## => Bootstrapped VaR_alpha computed manually

## Compute bootstrapped two-sided 95%-confidence intervals for VaR
VaR.CI <- apply(VaR.boot, 1, function(x) quantile(x, probs = c(0.025, 0.975)))

## 7 Plot

Finally, letâ€™s display all results.

## Setup
yran <- range(X, mu, VaR, mu.predict, VaR.predict, VaR.CI)
myran <- max(abs(yran))
yran <- c(-myran, myran)
xran <- c(1, length(X) + m)

## Simulated data (X_t) and estimated conditional mean and VaR_alpha
plot(X, type = "l", xlim = xran, ylim = yran, xlab = "Time t", ylab = "",
main = "Data, fitted ARMA-GARCH process, VaR, VaR predictions and VaR CIs",
sub = paste0("Expected exceedances: ", btest$expected.exceed, " Actual exceedances: ", btest$actual.exceed, "   Test decision: ", btest\$uc.Decision))
lines(as.numeric(mu), col = adjustcolor("darkblue", alpha.f = 0.5)) # hat{\mu}_t
lines(VaR, col = "darkred") # estimated VaR_alpha

## 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))))