Forecasting volatility and Value at Risk

Today let's learn how to put a number on the most a position can lose in one day, and how to check that the number holds up.

Say you hold a position worth 1,000,000 in a fund that tracks the DAX, the German stock index. Your risk manager asks a simple question: how much can this position lose tomorrow?

Nobody can answer that exactly. But there is a useful way to answer it: on 99 days out of 100, the loss stays below a certain amount. That amount is called the Value at Risk, or VaR.

The size of that amount depends on how volatile the market is going to be tomorrow. Volatility is how widely the daily returns spread around their average, so a high-volatility market gives a larger VaR and a low-volatility market gives a smaller one. That means we first need a forecast of volatility, and then a way to turn it into a VaR.

And a VaR that nobody has checked is only a claim. So the last stage is to test it on days that were not used to fit the model.

That gives three stages, and each one uses the result of the one before it.

Each of the three boxes above is run on real DAX returns.

The DAX returns and the fitted GARCH(1,1) model

We start with the data. The EuStockMarkets dataset that comes with R holds the daily closing level of four European stock indices from 1991 to 1998, and the DAX is one of them.

What we model is not the closing level but the daily log return in percent: 100 times the change in the log of the closing level from one day to the next. For daily moves it is very close to the ordinary percentage change. That gives 1,859 returns.

We then split them in time order. The first 1,000 returns are used to fit the model. The last 859 are the held-out days: the fit does not use them, so they are the days we test the VaR on later.

Press Run.

RInteractive R
# Build the DAX log returns in percent and hold out the last 859 days dax <- EuStockMarkets[, "DAX"] ret <- diff(log(as.numeric(dax))) * 100 fit_ret <- ret[1:1000] test_ret <- ret[1001:1859] length(ret) #> [1] 1859 round(c(mean = mean(ret), sd = sd(ret)), 4) #> mean sd #> 0.0652 1.0301

  

The mean daily return is 0.0652% and the standard deviation is 1.0301%. That standard deviation is one number for the whole period. But volatility does not stay at one level, and the model we fit next describes how it changes from day to day.

The model is a GARCH(1,1). It says the variance of a day's return depends on what happened the day before. That variance is called the conditional variance, because it is the variance of the return given everything up to the previous day. Here is the model as a formula:

\[ \sigma_t^2 = \omega + \alpha\,\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2 \]

Read it from left to right:

  • \(\sigma_t^2\) is the conditional variance of day t.
  • \(\varepsilon_{t-1}\) is the demeaned return of the day before: that day's return minus the average return.
  • \(\sigma_{t-1}^2\) is the conditional variance of the day before.
  • \(\omega\), \(\alpha\) and \(\beta\) are the three coefficients the fit estimates. \(\omega\) is a constant, \(\alpha\) sets how strongly yesterday's squared demeaned return feeds in, and \(\beta\) sets how much of yesterday's variance carries over.

The function garch() in the tseries package fits this model. It has no equation for the mean, so it expects a series whose mean is zero. So we subtract the fit-sample mean from every return. We subtract the same number from the held-out returns too, because on those days the mean has to come from the fit sample as well.

RInteractive R
# Fit the GARCH(1,1) to the first 1,000 days and read its three coefficients suppressMessages(library(tseries)) mu <- mean(fit_ret) e_fit <- fit_ret - mu e_test <- test_ret - mu round(mu, 4) #> [1] 0.0214 fit <- garch(e_fit, order = c(1, 1), trace = FALSE) round(cbind(estimate = coef(fit), std_error = sqrt(diag(vcov(fit)))), 4) #> estimate std_error #> a0 0.1126 0.0255 #> a1 0.0550 0.0170 #> b1 0.8264 0.0394

  

tseries labels the three coefficients a0, a1 and b1, which are \(\omega\), \(\alpha\) and \(\beta\). The estimates are 0.1126, 0.0550 and 0.8264, and each one is more than 3 standard errors away from zero.

Two more numbers follow from the coefficients, and both matter later. The persistence is \(\alpha + \beta\): the closer it is to 1, the longer a jump in variance takes to fade. The long-run variance is \(\omega / (1 - \alpha - \beta)\): the average variance the model implies over a long stretch of days.

RInteractive R
# Name the coefficients and derive the persistence and the long-run variance omega <- coef(fit)[["a0"]] alpha <- coef(fit)[["a1"]] beta <- coef(fit)[["b1"]] persistence <- alpha + beta long_run <- omega / (1 - persistence) round(c(persistence = persistence, long_run_variance = long_run, long_run_sd = sqrt(long_run)), 4) #> persistence long_run_variance long_run_sd #> 0.8813 0.9493 0.9743

  

The persistence is 0.8813 and the long-run variance is 0.9493, which is a long-run standard deviation of 0.9743% a day.

So the fitted model gives one thing for every day: a conditional variance. That is the quantity we forecast.