library(reshape2)
library(ggplot2)
library(zoo)
library(tsgarch)
library(xts)
source("common/functions.r", chdir = TRUE)23 Implementing risk forecasts
The volatility models of the previous chapters do not by themselves tell us what we stand to lose. This chapter turns them into Value-at-Risk and Expected Shortfall, by four methods — historical simulation, EWMA, Gaussian GARCH and GARCH-t. We build each as a function in functions.r so the backtesting chapter can reuse it.
23.1 Risk measures: VaR and Expected Shortfall
Value-at-Risk (VaR) and Expected Shortfall (ES) are standard risk measures used throughout the financial industry:
VaR answers the question: “What is the loss threshold that we expect to be exceeded with probability \(\probability\) over the next day?” For example, a 5% daily VaR of £100,000 means there is a 5% chance of losing £100,000 or more tomorrow.
ES, also known as Conditional VaR, asks a second question. “What is the average loss in the worst \(\probability\) of outcomes?” For a continuous profit and loss distribution that equals the conditional mean beyond the VaR cutoff. In general it is the average of the quantiles from zero to \(\probability\), which is the convention Risk_HS implements below, and the distinction starts to matter as soon as ES is computed on an empirical distribution.
These measures are commonly used for:
- Regulatory capital requirements (Basel framework, Basel Committee on Banking Supervision 2019)
- Internal risk limits and position sizing
They also serve as inputs for portfolio optimisation, risk budgeting, stress testing and scenario analysis, though these applications are beyond the scope of this notebook.
Both measures are read off the distribution of future returns. Forecasting that distribution is the whole problem.
23.2 Data and libraries
The GARCH packages, and EmpiricalES from the shared function files.
import numpy as np
import pandas as pd
from scipy import stats
from arch import arch_model
import sys
sys.path.insert(0, 'common')
from functions import ProcessRawData, EmpiricalESusing CSV, DataFrames, Statistics, Distributions, Printf
using ARCHModels
include("common/functions.jl");data = ProcessRawData()
sp500 = data$sp500
y_r = sp500$ydata_py = ProcessRawData()
sp500_py = data_py['sp500']
y_py = sp500_py['y'].valuesdata_jl = ProcessRawData();
sp500_jl = data_jl["sp500"];
y_jl = Float64.(sp500_jl.y);23.3 Specification
Backtesting needs a small set of inputs — the probability, the portfolio value, the sample size, the estimation and testing windows and the EWMA decay. We use the techniques discussed in Section 16.5 and keep track of these inputs in a list that we call par for parameters.
One convention to fix before we start. Every risk function below turns a return-scale figure into currency by multiplying by the portfolio value, -quantile * value. Our returns are compound (log) returns, and for those this is a first-order approximation rather than an identity. A log return \(\CompoundReturns\) moves a position worth \(\PortfolioValue\) by \[\PortfolioValue\left(\exp(\CompoundReturns) - 1\right) \approx \PortfolioValue \CompoundReturns,\] since \(\exp(\CompoundReturns) - 1 \approx \CompoundReturns\) for small \(\CompoundReturns\). At daily magnitudes the gap is minor — a 2.33% loss on 1,000 is 23.3 under the linear map against 23.0 exactly, about 1% of the reported figure — but it grows with the size of the move, so it is larger at longer horizons and deeper in the tail. We keep the linear map throughout, because it is the convention in general use and it keeps the risk measure linear in the position, and we flag it here so the approximation is visible rather than buried. Simple returns need no such correction. There the multiplication is exact.
par = list()
par$probability = 0.05 # probability level (0.05 = 5% VaR)
par$value = 1000 # portfolio value in currency units
par$WE = 2000 # estimation window size (number of observations)
par$T = length(y_r) # total sample size
par$WT = par$T-par$WE # testing window size (T - WE)
par$lambda = 0.94 # EWMA decay; 0.94 is the RiskMetrics daily value
cat("The backtesting parameters are:",
"\n\tprobability = ",par$probability,
"\n\tvalue = ",par$value,
"\n\tWE = ",par$WE,
"\n\tT = ",par$T,
"\n\tWT = ",par$WT,
"\n\tlambda = ",par$lambda,
"\n"
)The backtesting parameters are:
probability = 0.05
value = 1000
WE = 2000
T = 8938
WT = 6938
lambda = 0.94par_py = {
'probability': 0.05, # probability level (0.05 = 5% VaR)
'value': 1000, # portfolio value in currency units
'WE': 2000, # estimation window size (number of observations)
'T': len(y_py), # total sample size
'lambda': 0.94 # EWMA decay; 0.94 is the RiskMetrics daily value
}
par_py['WT'] = par_py['T'] - par_py['WE'] # testing window size
print(f"The backtesting parameters are:")
print(f"\tprobability = {par_py['probability']}")
print(f"\tvalue = {par_py['value']}")
print(f"\tWE = {par_py['WE']}")
print(f"\tT = {par_py['T']}")
print(f"\tWT = {par_py['WT']}")
print(f"\tlambda = {par_py['lambda']}")The backtesting parameters are:
probability = 0.05
value = 1000
WE = 2000
T = 8938
WT = 6938
lambda = 0.94par_jl = Dict(
"probability" => 0.05, # probability level (0.05 = 5% VaR)
"value" => 1000, # portfolio value in currency units
"WE" => 2000, # estimation window size (number of observations)
"T" => length(y_jl), # total sample size
"lambda" => 0.94 # EWMA decay; 0.94 is the RiskMetrics daily value
);
par_jl["WT"] = par_jl["T"] - par_jl["WE"]; # testing window size
println("The backtesting parameters are:")
println("\tprobability = ", par_jl["probability"])
println("\tvalue = ", par_jl["value"])
println("\tWE = ", par_jl["WE"])
println("\tT = ", par_jl["T"])
println("\tWT = ", par_jl["WT"])
println("\tlambda = ", par_jl["lambda"])The backtesting parameters are:
probability = 0.05
value = 1000
WE = 2000
T = 8938
WT = 6938
lambda = 0.9423.4 Functions
Functions let us avoid repeating identical code when we test different parameters. They also allow us to reuse the same implementations in the backtesting chapter (Chapter 26).
We label these functions as Risk_* where * refers to any of the methods we are using, like HS or EWMA. We also put these functions into functions.r so we can use them later.
Each function takes three arguments. These are the returns y, the parameter list par and Print. If Print = TRUE, the function reports the forecast VaR and ES, and for the volatility-based methods the forecast volatility as well. Historical simulation has no forecast volatility, so it reports only VaR and ES.
23.5 Historical simulation — Risk_HS
Historical simulation is a method that assumes that history repeats itself. It uses the empirical distribution of the returns to find the \(\probability\)-th quantile.
We sort the values using sort():
Risk estimation, however, uses only the most recent \(\EstWindow\) observations. The figure below sorts just that estimation window rather than the full sample.
23.5.1 Risk_HS
The Risk_HS function sorts the most recent returns and takes the empirical quantile. The rank of the observation we want is \(k = \lceil \probability \EstWindow \rceil\), the ceiling rather than the truncation, so that the index is well defined whenever \(\probability \EstWindow\) is not a whole number. With this chapter’s values, \(0.05 \times 2000 = 100\), the two conventions agree. Quantile conventions differ between packages and textbooks, so a reported figure should say which one it uses.
ES needs a second decision that VaR does not. VaR reads off a single observation, but ES averages the tail, so we have to say how much probability mass the average covers. Averaging the worst \(k\) observations covers \(k / \EstWindow\), which is more than \(\probability\) whenever \(\probability \EstWindow\) is not a whole number. EmpiricalES() instead covers exactly \(\probability\). Every observation strictly inside the tail gets its full weight \(1/\EstWindow\), and the cutoff observation gets only the fraction needed to reach \(\probability\), \[\ES(\probability) = -\frac{1}{\EstWindow \probability}\left[\sum_{i=1}^{m} \CompoundReturns_{(i)} + (\EstWindow\probability - m)\,\CompoundReturns_{(m+1)}\right], \qquad m = \lfloor \EstWindow \probability \rfloor .\] When \(\EstWindow \probability\) is a whole number the fractional weight is zero and this is just the mean of the worst \(\EstWindow\probability\) returns, which is the case in this chapter. When \(\EstWindow \probability < 1\) the sum is empty and the fraction cancels the denominator, leaving the single worst return — the whole tail sits inside one observation, and no average can do better than report it.
23.5.1.1 ES on a discrete distribution
It would be reasonable to think this is a rounding detail. It is not. On a discrete distribution the two ways of defining ES — averaging the quantiles that make up the worst \(\probability\), or taking the conditional mean of everything at or beyond the \(\VaR\) cutoff — can give different answers, and only the first is coherent.
A small constructed example shows it. Take two positions, each of which loses 60 with probability 4% and nothing otherwise, and suppose their outcomes are independent. Look at \(\probability = 5\%\).
Each position on its own has a 4% chance of any loss at all, which is less than 5%, so its 5% quantile sits at zero and \(\VaR(5\%) = 0\). Now hold both, equally weighted. The portfolio loses 60 with probability \(0.04^2 = 0.16\%\), loses 30 if exactly one position fails, with probability \(2 \times 0.04 \times 0.96 = 7.68\%\), and otherwise nothing. The cumulative probability of losing at least 30 is 7.84%, which is above 5%, so the portfolio’s \(\VaR(5\%) = 30\). Diversifying has produced a portfolio VaR of 30 out of two positions whose VaRs were both zero — the subadditivity failure of VaR, in three lines of arithmetic.
We now compute ES on the same numbers, under both definitions.
L = 60; q = 0.04; p_ex = 0.05
# Each position: outcome and its probability, worst first
v1 = c(-L, 0); pr1 = c(q, 1 - q)
# Equally weighted portfolio of two independent copies
v2 = c(-L, -L/2, 0); pr2 = c(q^2, 2*q*(1-q), (1-q)^2)
# Quantile-average ES: take exactly p of probability mass, worst first
ES_qa = function(v, pr, p) {
take = pmin(pr, pmax(0, p - c(0, head(cumsum(pr), -1))))
-sum(v * take) / p
}
# Conditional-expectation ES: everything at or beyond the VaR cutoff
ES_ce = function(v, pr, p) {
k = which(cumsum(pr) >= p)[1]
-sum(v[1:k] * pr[1:k]) / sum(pr[1:k])
}
cat(sprintf("quantile-average : position %.3f portfolio %.3f subadditive: %s\n",
ES_qa(v1, pr1, p_ex), ES_qa(v2, pr2, p_ex),
ES_qa(v2, pr2, p_ex) <= ES_qa(v1, pr1, p_ex)))
cat(sprintf("conditional mean : position %.3f portfolio %.3f subadditive: %s\n",
ES_ce(v1, pr1, p_ex), ES_ce(v2, pr2, p_ex),
ES_ce(v2, pr2, p_ex) <= ES_ce(v1, pr1, p_ex)))quantile-average : position 48.000 portfolio 30.960 subadditive: TRUE
conditional mean : position 2.400 portfolio 30.612 subadditive: FALSEL, q, p_ex = 60, 0.04, 0.05
v1 = np.array([-L, 0.0]); pr1 = np.array([q, 1 - q])
v2 = np.array([-L, -L/2, 0.0]); pr2 = np.array([q**2, 2*q*(1-q), (1-q)**2])
def ES_qa(v, pr, p):
before = np.concatenate(([0.0], np.cumsum(pr)[:-1]))
take = np.minimum(pr, np.maximum(0.0, p - before))
return -np.sum(v * take) / p
def ES_ce(v, pr, p):
k = int(np.argmax(np.cumsum(pr) >= p))
return -np.sum(v[:k+1] * pr[:k+1]) / np.sum(pr[:k+1])
print(f"quantile-average : position {ES_qa(v1,pr1,p_ex):.3f} portfolio {ES_qa(v2,pr2,p_ex):.3f}"
f" subadditive: {ES_qa(v2,pr2,p_ex) <= ES_qa(v1,pr1,p_ex)}")
print(f"conditional mean : position {ES_ce(v1,pr1,p_ex):.3f} portfolio {ES_ce(v2,pr2,p_ex):.3f}"
f" subadditive: {ES_ce(v2,pr2,p_ex) <= ES_ce(v1,pr1,p_ex)}")quantile-average : position 48.000 portfolio 30.960 subadditive: True
conditional mean : position 2.400 portfolio 30.612 subadditive: FalseL, q, p_ex = 60, 0.04, 0.05
v1 = [-L, 0.0]; pr1 = [q, 1 - q]
v2 = [-L, -L/2, 0.0]; pr2 = [q^2, 2*q*(1-q), (1-q)^2]
function ES_qa(v, pr, p)
before = vcat(0.0, cumsum(pr)[1:end-1])
take = min.(pr, max.(0.0, p .- before))
return -sum(v .* take) / p
end
function ES_ce(v, pr, p)
k = findfirst(cumsum(pr) .>= p)
return -sum(v[1:k] .* pr[1:k]) / sum(pr[1:k])
end
@printf("quantile-average : position %.3f portfolio %.3f subadditive: %s\n",
ES_qa(v1,pr1,p_ex), ES_qa(v2,pr2,p_ex), ES_qa(v2,pr2,p_ex) <= ES_qa(v1,pr1,p_ex))
@printf("conditional mean : position %.3f portfolio %.3f subadditive: %s\n",
ES_ce(v1,pr1,p_ex), ES_ce(v2,pr2,p_ex), ES_ce(v2,pr2,p_ex) <= ES_ce(v1,pr1,p_ex))quantile-average : position 48.000 portfolio 30.960 subadditive: true
conditional mean : position 2.400 portfolio 30.612 subadditive: falseThe quantile-average definition behaves. The portfolio figure comes in below the position figure, so diversification cannot increase measured risk. The conditional-mean definition does not. Because each position’s \(\VaR\) cutoff is zero, conditioning on “at or beyond the cutoff” sweeps in the entire distribution and returns a tiny number, while the portfolio’s cutoff sits at a genuine loss. The comparison then runs the wrong way and subadditivity fails.
The discrepancy comes from probability mass sitting exactly at the cutoff. The conditional-mean form has to take all of it or none. The quantile-average form takes only as much as is needed to reach \(\probability\). That is precisely the fractional weight in EmpiricalES(), which is why the code carries it. A finite sample of returns is a discrete distribution, so this is not a pathological case we can wave away — it is the case we are always in.
Risk_HS = function(y, par, Print = FALSE) {
ys = sort(tail(y, par$WE))
k = ceiling(par$probability * par$WE)
if (k < 1 || k > par$WE) stop("probability * WE must give a rank between 1 and WE")
VaR = -ys[k] * par$value
ES = -EmpiricalES(ys, par$probability) * par$value
if (Print) {
cat(sprintf("The HS %s%% VaR is $%.1f, while the ES is $%.1f\n",
par$probability * 100, VaR, ES))
}
return(list(VaR = VaR, ES = ES, method = "HS", par = par))
}def Risk_HS(y, par, Print=False):
ys = np.sort(y[-par['WE']:])
k = int(np.ceil(par['probability'] * par['WE']))
if k < 1 or k > par['WE']:
raise ValueError("probability * WE must give a rank between 1 and WE")
VaR = -ys[k - 1] * par['value'] # -1 for 0-based indexing
ES = -EmpiricalES(ys, par['probability']) * par['value']
if Print:
print(f"The HS {int(par['probability']*100)}% VaR is ${VaR:.1f}, while the ES is ${ES:.1f}")
return {'VaR': VaR, 'ES': ES, 'method': 'HS'}function Risk_HS(y, par; Print=false)
ys = sort(y[end-par["WE"]+1:end])
k = ceil(Int, par["probability"] * par["WE"])
(1 <= k <= par["WE"]) || error("probability * WE must give a rank between 1 and WE")
VaR = -ys[k] * par["value"]
ES = -EmpiricalES(ys, par["probability"]) * par["value"]
if Print
println(@sprintf("The HS %d%% VaR is \$%.1f, while the ES is \$%.1f",
Int(par["probability"]*100), VaR, ES))
end
return Dict("VaR" => VaR, "ES" => ES, "method" => "HS")
end;23.6 EWMA — Risk_EWMA
Historical simulation gives equal weight to all observations in the estimation window and is slow to adapt when market conditions change. The EWMA approach addresses this by giving more weight to recent observations, allowing for faster adaptation to changing volatility.
The Risk_EWMA function uses exponentially weighted moving averages to forecast volatility. It assumes returns follow a normal distribution with time-varying volatility.
23.6.1 Normal distribution functions
Since EWMA and GARCH methods assume normally distributed returns, we need functions to calculate VaR and ES from the normal distribution:
NormalVaR = function(p, sd = 1) return(-qnorm(p, sd = sd))
NormalES = function(p, sd = 1) return(sd * dnorm(qnorm(p)) / p)def NormalVaR(p, sd=1):
return -stats.norm.ppf(p, scale=sd)
def NormalES(p, sd=1):
return sd * stats.norm.pdf(stats.norm.ppf(p)) / pfunction NormalVaR(p; sd=1)
return -quantile(Normal(0, sd), p)
end
function NormalES(p; sd=1)
return sd * pdf(Normal(), quantile(Normal(), p)) / p
end;These functions set the conditional mean to zero. At daily horizons the mean is well inside estimation error, so little is lost by leaving it out, and the same assumption carries through the Student-t versions below.
23.6.2 Risk_EWMA
The recursion below squares the raw returns rather than deviations from the mean, which is the RiskMetrics convention (J.P. Morgan 1994). The two GARCH functions de-mean the window first, so the treatments differ, but at daily frequency the numerical effect is negligible.
Risk_EWMA = function(y, par, Print = FALSE) {
y = tail(y, par$WE)
sigma2 = var(y)
for (i in 2:length(y)) {
sigma2 = par$lambda * sigma2 + (1 - par$lambda) * y[i - 1]^2
}
sigma2 = par$lambda * sigma2 + (1 - par$lambda) * y[length(y)]^2
sigma = sqrt(sigma2)
VaR = NormalVaR(p = par$probability, sd = sigma) * par$value
ES = NormalES(p = par$probability, sd = sigma) * par$value
if (Print) {
cat(sprintf("The EWMA forecast volatility is %.4f\n", sigma))
cat(sprintf("The EWMA %s%% VaR is $%.1f, while the ES is $%.1f\n",
par$probability * 100, VaR, ES))
}
return(list(VaR = VaR, ES = ES, method = "EWMA", sigma = sigma))
}def Risk_EWMA(y, par, Print=False):
y = y[-par['WE']:]
sigma2 = np.var(y, ddof=1) # ddof=1 for sample variance
for i in range(1, len(y)):
sigma2 = par['lambda'] * sigma2 + (1 - par['lambda']) * y[i-1]**2
sigma2 = par['lambda'] * sigma2 + (1 - par['lambda']) * y[-1]**2
sigma = np.sqrt(sigma2)
VaR = NormalVaR(par['probability'], sd=sigma) * par['value']
ES = NormalES(par['probability'], sd=sigma) * par['value']
if Print:
print(f"The EWMA forecast volatility is {sigma:.4f}")
print(f"The EWMA {int(par['probability']*100)}% VaR is ${VaR:.1f}, while the ES is ${ES:.1f}")
return {'VaR': VaR, 'ES': ES, 'method': 'EWMA', 'sigma': sigma}function Risk_EWMA(y, par; Print=false)
y = y[end-par["WE"]+1:end]
sigma2 = var(y)
for i in 2:length(y)
sigma2 = par["lambda"] * sigma2 + (1 - par["lambda"]) * y[i-1]^2
end
sigma2 = par["lambda"] * sigma2 + (1 - par["lambda"]) * y[end]^2
sigma = sqrt(sigma2)
VaR = NormalVaR(par["probability"], sd=sigma) * par["value"]
ES = NormalES(par["probability"], sd=sigma) * par["value"]
if Print
println(@sprintf("The EWMA forecast volatility is %.4f", sigma))
println(@sprintf("The EWMA %d%% VaR is \$%.1f, while the ES is \$%.1f",
Int(par["probability"]*100), VaR, ES))
end
return Dict("VaR" => VaR, "ES" => ES, "method" => "EWMA", "sigma" => sigma)
end;23.7 Gaussian GARCH — Risk_nGARCH
The GARCH approach models volatility clustering by estimating how current volatility depends on past volatility and past squared returns. For details on GARCH estimation, see Chapter 18. The chapter 5 slides call this model Normal GARCH, which is the same specification under a different name.
23.7.1 Risk_nGARCH
The Risk_nGARCH function estimates a GARCH(1,1) model to capture volatility clustering in returns. It uses maximum likelihood estimation and assumes normally distributed innovations.
Risk_nGARCH = function(y, par, Print = FALSE) {
ys = tail(y, par$WE)
ys = ys - mean(ys)
ys_xts = xts(ys, order.by = seq(as.Date("2000-01-01"), by = "day", length.out = length(ys)))
spec = garch_modelspec(ys_xts, model = "garch", order = c(1, 1), constant = FALSE)
res = estimate(spec)
pred = predict(res, h = 1)
sigma = as.numeric(pred$sigma)
params = coef(res)
VaR = NormalVaR(p = par$probability, sd = sigma) * par$value
ES = NormalES(p = par$probability, sd = sigma) * par$value
if (Print) {
cat(sprintf("The GARCH forecast volatility is %.4f\n", sigma))
cat(sprintf("The GARCH %s%% VaR is $%.1f, while the ES is $%.1f\n",
par$probability * 100, VaR, ES))
}
return(list(VaR = VaR, ES = ES, method = "nGARCH", parameters = params, sigma = sigma))
}def Risk_nGARCH(y, par, Print=False):
ys = y[-par['WE']:]
ys = ys - np.mean(ys)
# Scale for numerical stability (arch package works better with percentage returns)
scale = 100
ys_scaled = ys * scale
model = arch_model(ys_scaled, vol='Garch', p=1, q=1, mean='Zero', dist='normal')
res = model.fit(disp='off')
forecast = res.forecast(horizon=1)
sigma = np.sqrt(forecast.variance.values[-1, 0]) / scale
VaR = NormalVaR(par['probability'], sd=sigma) * par['value']
ES = NormalES(par['probability'], sd=sigma) * par['value']
if Print:
print(f"The GARCH forecast volatility is {sigma:.4f}")
print(f"The GARCH {int(par['probability']*100)}% VaR is ${VaR:.1f}, while the ES is ${ES:.1f}")
return {'VaR': VaR, 'ES': ES, 'method': 'nGARCH', 'sigma': sigma}function Risk_nGARCH(y, par; Print=false)
ys = y[end-par["WE"]+1:end]
ys = ys .- mean(ys)
# Fit GARCH(1,1) with normal distribution
model = fit(GARCH{1,1}, ys; meanspec=NoIntercept())
# Get one-step ahead forecast
fc = predict(model, :variance, 1)
sigma = sqrt(fc[1])
VaR = NormalVaR(par["probability"], sd=sigma) * par["value"]
ES = NormalES(par["probability"], sd=sigma) * par["value"]
if Print
println(@sprintf("The GARCH forecast volatility is %.4f", sigma))
println(@sprintf("The GARCH %d%% VaR is \$%.1f, while the ES is \$%.1f",
Int(par["probability"]*100), VaR, ES))
end
return Dict("VaR" => VaR, "ES" => ES, "method" => "nGARCH", "sigma" => sigma)
end;GARCH models capture volatility clustering, but they assume normally distributed innovations. This means returns are conditionally normal (given the volatility at time t) but unconditionally exhibit fat tails due to the time-varying volatility.
23.8 GARCH-t — Risk_tGARCH
The Risk_tGARCH function extends the GARCH model by using a Student-t distribution for the innovations instead of a normal distribution. This can better capture the heavy tails commonly observed in financial returns. The function keeps the name Risk_tGARCH throughout the pipeline, while the prose follows the slides and calls the model GARCH-t.
23.8.1 Student-t distribution functions
Whether the Student-t distribution is standardised to unit variance changes the quantile. A raw Student-t variate has variance \(\DOF/(\DOF-2)\), so most GARCH libraries standardise it to unit variance, which requires dividing the raw quantile by \(\sqrt{\DOF/(\DOF-2)}\), where \(\DOF\) is the degrees of freedom parameter. That division is what the code below performs.
That normalisation carries a restriction. The variance \(\DOF/(\DOF-2)\) only exists for \(\DOF > 2\), and at \(\DOF = 2\) exactly the divisor is infinite, so an unguarded implementation returns a VaR of zero, reporting an undefined variance as an absence of risk.
ES needs a second condition. The tail mean of a raw Student-t is finite only for \(\DOF > 1\), so the unstandardised branch needs that even though it skips the normalisation. Both functions check their own condition and stop rather than return a number that looks usable. Fitted GARCH packages normally constrain \(\DOF\) to an admissible range, so this guards against a hand-supplied value, not against the estimates below.
tVaR = function(p, df, sd = 1, Standardized = TRUE) {
if (Standardized && df <= 2) stop("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
Scale = 1
if (Standardized) Scale = sqrt(df / (df - 2))
return(-sd * qt(p = p, df = df) / Scale)
}
tES = function(p, df, sd = 1, Standardized = TRUE) {
if (Standardized && df <= 2) stop("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
if (df <= 1) stop("Student-t ES needs df > 1: the tail mean is infinite at or below 1")
Scale = 1
if (Standardized) Scale = sqrt(df / (df - 2))
ES = (sd * dt(qt(p = p, df = df), df = df) *
((df + (qt(p = p, df = df))^2) / (df - 1)) / p) / Scale
return(ES)
}def tVaR(p, df, sd=1, Standardized=True):
if Standardized and df <= 2:
raise ValueError("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
Scale = 1
if Standardized:
Scale = np.sqrt(df / (df - 2))
return -sd * stats.t.ppf(p, df) / Scale
def tES(p, df, sd=1, Standardized=True):
if Standardized and df <= 2:
raise ValueError("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
if df <= 1:
raise ValueError("Student-t ES needs df > 1: the tail mean is infinite at or below 1")
Scale = 1
if Standardized:
Scale = np.sqrt(df / (df - 2))
q = stats.t.ppf(p, df)
ES = (sd * stats.t.pdf(q, df) * ((df + q**2) / (df - 1)) / p) / Scale
return ESfunction tVaR(p, df; sd=1, Standardized=true)
(!Standardized || df > 2) || error("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
Scale = 1
if Standardized
Scale = sqrt(df / (df - 2))
end
return -sd * quantile(TDist(df), p) / Scale
end
function tES(p, df; sd=1, Standardized=true)
(!Standardized || df > 2) || error("Standardized Student-t needs df > 2: the variance is undefined at or below 2")
df > 1 || error("Student-t ES needs df > 1: the tail mean is infinite at or below 1")
Scale = 1
if Standardized
Scale = sqrt(df / (df - 2))
end
q = quantile(TDist(df), p)
ES = (sd * pdf(TDist(df), q) * ((df + q^2) / (df - 1)) / p) / Scale
return ES
end;23.8.2 How much do fat tails move the numbers?
Before fitting anything, we can see what the degrees of freedom actually buy. Both functions return a multiplier on \(\Vol \PortfolioValue\), so we can tabulate them directly and compare against the normal. Everything is at \(\probability = 1\%\) and unit variance, so only tail shape differs.
p_mult = 0.01
dfs = c(4, 5, 6, 30)
tab = data.frame(
Distribution = c("Normal", sprintf("Student-t(%d)", dfs)),
VaR = c(NormalVaR(p_mult), sapply(dfs, function(d) tVaR(p_mult, d))),
ES = c(NormalES(p_mult), sapply(dfs, function(d) tES(p_mult, d)))
)
tab$`ES/VaR` = tab$ES / tab$VaR
tab$`VaR vs norm` = sprintf("%+.1f%%", 100 * (tab$VaR / tab$VaR[1] - 1))
tab$`ES vs norm` = sprintf("%+.1f%%", 100 * (tab$ES / tab$ES[1] - 1))
print(format(tab, digits = 4), row.names = FALSE) Distribution VaR ES ES/VaR VaR vs norm ES vs norm
Normal 2.326 2.665 1.146 +0.0% +0.0%
Student-t(4) 2.649 3.692 1.393 +13.9% +38.5%
Student-t(5) 2.606 3.449 1.323 +12.0% +29.4%
Student-t(6) 2.566 3.293 1.283 +10.3% +23.5%
Student-t(30) 2.374 2.768 1.166 +2.0% +3.9%p_mult = 0.01
dfs = [4, 5, 6, 30]
tab = pd.DataFrame({
'Distribution': ['Normal'] + [f'Student-t({d})' for d in dfs],
'VaR': [NormalVaR(p_mult)] + [tVaR(p_mult, d) for d in dfs],
'ES': [NormalES(p_mult)] + [tES(p_mult, d) for d in dfs],
})
tab['ES/VaR'] = tab['ES'] / tab['VaR']
tab['VaR vs norm'] = [f"{100*(v/tab['VaR'][0]-1):+.1f}%" for v in tab['VaR']]
tab['ES vs norm'] = [f"{100*(v/tab['ES'][0]-1):+.1f}%" for v in tab['ES']]
print(tab.to_string(index=False, float_format=lambda x: f"{x:.4f}")) Distribution VaR ES ES/VaR VaR vs norm ES vs norm
Normal 2.3263 2.6652 1.1457 +0.0% +0.0%
Student-t(4) 2.6495 3.6915 1.3933 +13.9% +38.5%
Student-t(5) 2.6065 3.4488 1.3232 +12.0% +29.4%
Student-t(6) 2.5660 3.2925 1.2832 +10.3% +23.5%
Student-t(30) 2.3739 2.7685 1.1662 +2.0% +3.9%p_mult = 0.01
dfs = [4, 5, 6, 30]
names_jl = vcat("Normal", ["Student-t($d)" for d in dfs])
VaRs = vcat(NormalVaR(p_mult), [tVaR(p_mult, d) for d in dfs])
ESs = vcat(NormalES(p_mult), [tES(p_mult, d) for d in dfs])
@printf("%-15s %8s %8s %8s %12s %12s\n", "Distribution", "VaR", "ES", "ES/VaR", "VaR vs norm", "ES vs norm")
for i in eachindex(names_jl)
@printf("%-15s %8.4f %8.4f %8.4f %11.1f%% %11.1f%%\n", names_jl[i], VaRs[i], ESs[i],
ESs[i]/VaRs[i], 100*(VaRs[i]/VaRs[1]-1), 100*(ESs[i]/ESs[1]-1))
endDistribution VaR ES ES/VaR VaR vs norm ES vs norm
Normal 2.3263 2.6652 1.1457 0.0% 0.0%
Student-t(4) 2.6495 3.6915 1.3933 13.9% 38.5%
Student-t(5) 2.6065 3.4488 1.3232 12.0% 29.4%
Student-t(6) 2.5660 3.2925 1.2832 10.3% 23.5%
Student-t(30) 2.3739 2.7685 1.1662 2.0% 3.9%Two things stand out. The ES multiplier rises much faster than the VaR multiplier as \(\DOF\) falls. At \(\DOF = 4\), VaR is about 14% above the normal figure while ES is about 39% above. That is what it means for ES to average over the whole tail rather than read one point on it — extra tail mass beyond the quantile lifts ES without moving VaR nearly as much. Anyone comparing a normal and a Student-t model on VaR alone is looking at the least sensitive part of the comparison.
Second, \(\DOF = 30\) is already close to the normal, at 2% on VaR. The distinction matters at the single-digit degrees of freedom that daily equity returns actually produce, and stops mattering quickly above that. ### Risk_tGARCH
Risk_tGARCH = function(y, par, Print = FALSE) {
ys = tail(y, par$WE)
ys = ys - mean(ys)
ys_xts = xts(ys, order.by = seq(as.Date("2000-01-01"), by = "day", length.out = length(ys)))
spec = garch_modelspec(ys_xts, model = "garch", order = c(1, 1),
constant = FALSE, distribution = "std")
res = estimate(spec)
pred = predict(res, h = 1)
sigma = as.numeric(pred$sigma)
params = coef(res)
nu = params["shape"]
VaR = tVaR(p = par$probability, sd = sigma, df = nu, Standardized = TRUE) * par$value
ES = tES(p = par$probability, sd = sigma, df = nu, Standardized = TRUE) * par$value
if (Print) {
cat(sprintf("The GARCH-t forecast volatility is %.4f\n", sigma))
cat(sprintf("The GARCH-t degrees of freedom is %.2f\n", nu))
cat(sprintf("The GARCH-t %s%% VaR is $%.1f, while the ES is $%.1f\n",
par$probability * 100, VaR, ES))
}
return(list(VaR = VaR, ES = ES, method = "tGARCH", parameters = params, sigma = sigma, nu = nu))
}def Risk_tGARCH(y, par, Print=False):
ys = y[-par['WE']:]
ys = ys - np.mean(ys)
# Scale for numerical stability
scale = 100
ys_scaled = ys * scale
model = arch_model(ys_scaled, vol='Garch', p=1, q=1, mean='Zero', dist='t')
res = model.fit(disp='off')
forecast = res.forecast(horizon=1)
sigma = np.sqrt(forecast.variance.values[-1, 0]) / scale
nu = res.params['nu']
VaR = tVaR(par['probability'], nu, sd=sigma, Standardized=True) * par['value']
ES = tES(par['probability'], nu, sd=sigma, Standardized=True) * par['value']
if Print:
print(f"The GARCH-t forecast volatility is {sigma:.4f}")
print(f"The GARCH-t degrees of freedom is {nu:.2f}")
print(f"The GARCH-t {int(par['probability']*100)}% VaR is ${VaR:.1f}, while the ES is ${ES:.1f}")
return {'VaR': VaR, 'ES': ES, 'method': 'tGARCH', 'sigma': sigma, 'nu': nu}function Risk_tGARCH(y, par; Print=false)
ys = y[end-par["WE"]+1:end]
ys = ys .- mean(ys)
# Fit GARCH(1,1) with Student-t distribution
model = fit(GARCH{1,1}, ys; meanspec=NoIntercept(), dist=StdT)
# Get one-step ahead forecast
fc = predict(model, :variance, 1)
sigma = sqrt(fc[1])
# Extract degrees of freedom from the model
nu = model.dist.coefs[1]
VaR = tVaR(par["probability"], nu, sd=sigma, Standardized=true) * par["value"]
ES = tES(par["probability"], nu, sd=sigma, Standardized=true) * par["value"]
if Print
println(@sprintf("The GARCH-t forecast volatility is %.4f", sigma))
println(@sprintf("The GARCH-t degrees of freedom is %.2f", nu))
println(@sprintf("The GARCH-t %d%% VaR is \$%.1f, while the ES is \$%.1f",
Int(par["probability"]*100), VaR, ES))
end
return Dict("VaR" => VaR, "ES" => ES, "method" => "tGARCH", "sigma" => sigma, "nu" => nu)
end;23.9 Application and method comparison
We now apply all four risk forecasting methods to the S&P 500 data and compare their VaR and ES estimates.
res_r = list()
res_r$HS = Risk_HS(y = y_r, par = par, Print = TRUE)
res_r$EWMA = Risk_EWMA(y = y_r, par = par, Print = TRUE)
res_r$nGARCH = Risk_nGARCH(y = y_r, par = par, Print = TRUE)
res_r$tGARCH = Risk_tGARCH(y = y_r, par = par, Print = TRUE)The HS 5% VaR is $18.3, while the ES is $30.5
The EWMA forecast volatility is 0.0098
The EWMA 5% VaR is $16.1, while the ES is $20.1
The GARCH forecast volatility is 0.0074
The GARCH 5% VaR is $12.2, while the ES is $15.3
The GARCH-t forecast volatility is 0.0074
The GARCH-t degrees of freedom is 5.71
The GARCH-t 5% VaR is $11.7, while the ES is $16.5res_py = {}
res_py['HS'] = Risk_HS(y_py, par_py, Print=True)
res_py['EWMA'] = Risk_EWMA(y_py, par_py, Print=True)
res_py['nGARCH'] = Risk_nGARCH(y_py, par_py, Print=True)
res_py['tGARCH'] = Risk_tGARCH(y_py, par_py, Print=True)The HS 5% VaR is $18.3, while the ES is $30.5
The EWMA forecast volatility is 0.0098
The EWMA 5% VaR is $16.1, while the ES is $20.1
The GARCH forecast volatility is 0.0074
The GARCH 5% VaR is $12.2, while the ES is $15.3
The GARCH-t forecast volatility is 0.0074
The GARCH-t degrees of freedom is 5.76
The GARCH-t 5% VaR is $11.7, while the ES is $16.4res_jl = Dict{String, Any}();
res_jl["HS"] = Risk_HS(y_jl, par_jl, Print=true);
res_jl["EWMA"] = Risk_EWMA(y_jl, par_jl, Print=true);
res_jl["nGARCH"] = Risk_nGARCH(y_jl, par_jl, Print=true);
res_jl["tGARCH"] = Risk_tGARCH(y_jl, par_jl, Print=true);The HS 5% VaR is $18.3, while the ES is $30.5
The EWMA forecast volatility is 0.0098
The EWMA 5% VaR is $16.1, while the ES is $20.1
The GARCH forecast volatility is 0.0074
The GARCH 5% VaR is $12.2, while the ES is $15.3
The GARCH-t forecast volatility is 0.0074
The GARCH-t degrees of freedom is 5.69
The GARCH-t 5% VaR is $11.7, while the ES is $16.5The methods do not tell the same story. Each one treats recent history and tail risk differently, so the forecasts diverge.
Historical simulation reflects only what happened inside the estimation window. It imposes no distributional assumption and inherits whatever tail the window happens to contain, which is its strength and its weakness. It also moves slowly. A large loss keeps its full weight for as long as it stays in the window and then drops out all at once, so the estimate can sit high long after markets have calmed. That is why it gives the largest VaR here.
EWMA responds much faster, because the exponential weights let the most recent days dominate. The cost is a normal distribution and no mean reversion, so the forecast follows the recent past without ever being pulled back towards a long-run level.
Gaussian GARCH adds that mean reversion and estimates the volatility dynamics from the data instead of fixing the decay in advance. It still assumes normal innovations, so the shape of the tail comes entirely from the volatility path.
GARCH-t replaces those normal innovations with a Student-t. This raises VaR far out in the tail, at 1% and beyond, but at this chapter’s 5% level it gives a slightly lower VaR than Gaussian GARCH. The standardised Student-t and normal densities cross, and closer in the Student-t quantile is the smaller of the two. The fat tails show up instead in ES, which averages over the whole tail and is duly higher for GARCH-t than for Gaussian GARCH.
Formal backtesting (Chapter 26, Chapter 27) provides empirical guidance on which methods perform best for a given dataset.
23.10 Cross-language comparison
To verify consistency across implementations, we compare the VaR and ES estimates from all three languages.
23.10.1 VaR comparison
First, do the four methods agree with each other? One row per method, one column per language.
Method VaR_Julia VaR_Python VaR_R
2 HS 18.32 18.32 18.32
1 EWMA 16.06 16.06 16.06
3 nGARCH 12.23 12.22 12.23
4 tGARCH 11.72 11.72 11.7223.10.2 ES comparison
Second, does the tail behave differently from the quantile? Each ES sits above its own VaR, as it must, since it averages the losses beyond that threshold, but by how much varies with the method.
Method ES_Julia ES_Python ES_R
2 HS 30.53 30.53 30.53
1 EWMA 20.14 20.14 20.14
3 nGARCH 15.34 15.33 15.33
4 tGARCH 16.48 16.43 16.4623.10.3 Full comparison table
Third, the two together, which is how a risk report would present them. Reading across a row shows how far a method’s tail average sits from its threshold.
Method VaR_Julia VaR_Python VaR_R ES_Julia ES_Python ES_R
2 HS 18.32 18.32 18.32 30.53 30.53 30.53
1 EWMA 16.06 16.06 16.06 20.14 20.14 20.14
3 nGARCH 12.23 12.22 12.23 15.34 15.33 15.33
4 tGARCH 11.72 11.72 11.72 16.48 16.43 16.46The HS and EWMA methods produce identical results across all three languages since they use the same algorithms. The Gaussian GARCH and GARCH-t results show minor differences due to:
- Different optimisation algorithms in each package
- Different starting values for parameter estimation
- Numerical precision differences in the maximum likelihood estimation
For practical risk work, these differences are small enough to ignore.
23.11 Further directions
The four methods implemented above cover the standard toolkit for daily risk forecasting. Several extensions address limitations that arise in practice, particularly at extreme quantiles and when combining methods. This section outlines the main directions. The implementations are beyond the scope of this chapter.
23.11.1 Limitations at extreme quantiles
The methods above work well for moderate probability levels like 5% VaR. At more extreme quantiles (1% or 0.1%), additional challenges arise:
Historical simulation relies on having enough observations in the tail region. For 1% VaR with 1,000 observations, only about 10 data points determine the estimate. For 0.1% VaR, we need 1,000 observations just to expect a single exceedance.
Parametric methods are sensitive to distributional assumptions precisely where we have least data to validate them. The difference between normal and Student-t tails becomes substantial at extreme quantiles.
Backtesting becomes difficult because violations are rare. A 1% VaR should be exceeded about 2-3 times per year with daily data. Statistical tests have low power with so few observations.
23.11.2 Filtered historical simulation
Filtered historical simulation (FHS, Barone-Adesi et al. 1999) addresses these limitations by combining GARCH’s volatility dynamics with historical simulation’s empirical tail shape:
- Fit a GARCH model to estimate conditional volatility \(\Vol_t\)
- Compute standardised residuals \(\StdNormal_t = (\CompoundReturns_t - \Mean_t) / \Vol_t\), where \(\Mean_t\) is the conditional mean (often set to zero for daily returns)
- Resample from the empirical distribution of \(\StdNormal_t\)
- Scale by the current volatility forecast to get simulated returns
This approach captures volatility dynamics through GARCH while preserving the empirical shape of the residual distribution. It adapts to current market conditions (unlike pure HS) without imposing a parametric distribution on the tails (unlike pure GARCH).
The main implementation choice is the GARCH specification. The residual distribution inherits whatever tail properties remain after volatility scaling.
23.11.3 Extreme value theory
For very extreme quantiles, extreme value theory (EVT, Embrechts et al. 1997) provides tools tailored to tail estimation. Under weak regularity conditions, exceedances above a sufficiently high threshold are approximately distributed as a generalised Pareto distribution (GPD), and the approximation improves as the threshold rises. Threshold choice and diagnostic checks are an integral part of the method.
The peaks-over-threshold approach:
- Choose a high threshold \(\Threshold\) (e.g., the 95th percentile of losses)
- Fit a GPD to observations exceeding \(\Threshold\)
- Use the fitted distribution to extrapolate to more extreme quantiles
EVT is most useful when:
- Sample sizes are large (1,000+ observations)
- The target quantile is very extreme (0.1% or beyond)
- Unconditional risk estimates are acceptable
For conditional risk with moderate quantiles, GARCH-based methods are usually more practical.
23.11.4 Practical guidance
Sample size does most of the work when estimating tail risk. An elaborate method that yields unreliable tail estimates is often worse than a simpler method whose wider confidence intervals are at least honest about the uncertainty.
Compare estimates across methods and parameter choices, because a large discrepancy between them is itself a measure of how uncertain the estimate is. The application matters as well, since stress testing may warrant conservative assumptions where daily trading limits favour responsive methods. Whichever method is chosen, tail estimates carry substantial sampling error and that should inform how they are used.
23.12 Exercises
Create a generic function for obtaining VaR and ES from GARCH models, where the distribution (normal, t or skew t), the lag order and the apARCH option are arguments to a single function.
Combine all risk functions discussed above into a single function,
Risk = function(y,par,Model,Print), that computes ES and VaR based on a set of inputs. This function should support the HS, EWMA, Gaussian GARCH and GARCH-t models and be robust enough to handle improper inputs.
Barone-Adesi, Giovanni, Kostas Giannopoulos, and Les Vosper. 1999. “VaR Without Correlations for Portfolios of Derivative Securities.” Journal of Futures Markets 19 (5): 583–602.
Basel Committee on Banking Supervision. 2019. Minimum Capital Requirements for Market Risk. BCBS 457. Bank for International Settlements.
Embrechts, P., C. Klüppelberg, and T. Mikosch. 1997. Modelling Extremal Events for Insurance and Finance (Applications of Mathematics). Berlin: Springer Verlag.
J.P. Morgan. 1994. RiskMetrics Technical Document. Second. J.P. Morgan.