Week 8: Implementing Risk Forecasting

Version 5.0 - August 2026

Author
Affiliation

Jon Danielsson

London School of Economics

8 Week 8

Value at Risk (VaR) asks for a loss threshold. How large is the loss that will be exceeded with a chosen small probability? It is not the largest possible loss, and it says nothing about how far beyond the threshold a loss might go. Expected Shortfall (ES) answers that second question, asking how large the loss is on average once the threshold has been crossed.

This week calculates VaR and ES using historical simulation, extends historical simulation to a portfolio, and compares VaR estimates from historical simulation, EWMA and GARCH. The methods make different assumptions, so they do not have to give the same answer.

NoteBefore you start

Open the main course folder and check that Returns.csv, created in Week 2, is present. You will use quantiles from Week 3, EWMA from Week 7 and GARCH from Week 6. The arch library must be installed in the Anaconda Python selected in Positron.

TipBy the end of this session

You should be able to calculate historical VaR and ES, calculate portfolio historical VaR, produce EWMA and GARCH VaR forecasts, and explain why the three methods differ.

8.1 The plan for this week

  1. Calculate historical VaR and ES for one asset
  2. Calculate historical VaR and ES for a portfolio
  3. Forecast VaR with EWMA
  4. Forecast VaR with GARCH
  5. Compare the methods

8.2 Loading data and setting choices

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from arch import arch_model

Returns = pd.read_csv("Returns.csv", index_col="date", parse_dates=True)
y = Returns["JPM"].dropna()
portfolio_returns = Returns[["JPM", "MCD"]].dropna()
p = 0.01
value = 1000
estimation_window = 1000
lambda_ = 0.94
CautionUnits and signs used this week
  • Returns are one-day decimals: 1% is 0.01.
  • p = 0.01 means the lower 1% return tail, often called 99% VaR.
  • VaR and ES are reported as positive loss amounts in dollars.
  • A later VaR violation occurs when return < -VaR / value.
  • Volatility forecasts are daily unless explicitly annualised.
NoteCheck your result

y should contain at least 1,000 observations. The portfolio data should have two columns with the same dates and no missing values.

assert len(y) >= estimation_window
assert portfolio_returns.shape[1] == 2
print(y.index.min(), "to", y.index.max())
2000-01-03 00:00:00 to 2025-09-23 00:00:00

8.3 Historical VaR and ES

Historical simulation treats the latest observed returns as possible outcomes for tomorrow. It does not assume a normal distribution.

method="lower" makes the threshold an observed return rather than an interpolation between two of them, which is what historical simulation is supposed to deliver. numpy offers several finite-sample quantile conventions and interpolates by default. With 1,000 observations and \(p=0.01\) the tail holds about ten returns, with 500 about five, and with 250 about three, which is worth remembering when the estimation window is shortened.

estimation_returns = y.tail(estimation_window).to_numpy()
return_threshold = np.quantile(estimation_returns, p, method="lower")
tail_returns = estimation_returns[estimation_returns <= return_threshold]

hs_var = -value * return_threshold
hs_es = -value * tail_returns.mean()

print("Historical VaR: $", round(hs_var, 2))
print("Historical ES:  $", round(hs_es, 2))
Historical VaR: $ 45.61
Historical ES:  $ 57.63
Risk idea Formula Python name
Lower-tail return threshold \(q_p\) return_threshold
Positive VaR loss \(-Vq_p\) hs_var
Returns beyond VaR \(r \leq q_p\) tail_returns
Positive average tail loss \(-V E[r\mid r\leq q_p]\) hs_es
NoteCheck your result

Both numbers should be positive. ES should be at least as large as VaR because it averages losses from beyond the VaR threshold.

assert hs_var > 0
assert hs_es >= hs_var
print("tail observations used for ES:", len(tail_returns))
tail observations used for ES: 10

The histogram shows where the VaR threshold sits in the return distribution.

fig, ax = plt.subplots()
ax.hist(estimation_returns, bins=50, edgecolor="white")
ax.axvline(return_threshold, color="red", label="1% return threshold")
ax.set_title("JPM returns used for historical simulation")
ax.set_xlabel("Daily return")
ax.set_ylabel("Count")
ax.legend(frameon=False)
plt.show()

8.3.1 A function for historical VaR and ES

This function packages the same steps so that we can use them again for a portfolio. Run it unchanged and focus on what goes in and what comes out.

def historical_var_es(returns, probability=0.01, portfolio_value=1000):
    values = np.asarray(returns)
    threshold = np.quantile(values, probability, method="lower")
    tail = values[values <= threshold]
    var = -portfolio_value * threshold
    es = -portfolio_value * tail.mean()
    return {"VaR": var, "ES": es, "threshold": threshold}
historical_var_es(
    estimation_returns,
    probability=p,
    portfolio_value=value,
)
{'VaR': np.float64(45.6079093685133),
 'ES': np.float64(57.63329358114164),
 'threshold': np.float64(-0.0456079093685133)}

8.4 Portfolio historical VaR and ES

A portfolio return is a weighted sum of its asset returns. The example invests 10% in JPM and 90% in MCD. These are log returns, and a weighted sum of log returns is not exactly the portfolio’s log return. Over one day the difference is negligible, and the exact version would convert each return with np.exp(r) - 1 before applying the weights.

weights = np.array([0.10, 0.90])
portfolio_window = portfolio_returns.tail(estimation_window)
weighted_returns = portfolio_window.to_numpy() @ weights

portfolio_hs = historical_var_es(
    weighted_returns,
    probability=p,
    portfolio_value=value,
)
portfolio_hs
{'VaR': np.float64(28.368866684977093),
 'ES': np.float64(38.76966539469716),
 'threshold': np.float64(-0.028368866684977092)}
NoteCheck your result

The weights should sum to one. Portfolio ES should be at least as large as portfolio VaR. Do not add the two individual VaRs. Historical simulation keeps the way the two returns moved together on each day.

assert np.isclose(weights.sum(), 1)
assert portfolio_hs["ES"] >= portfolio_hs["VaR"]

8.5 EWMA VaR

EWMA assumes a normal one-day return but lets its variance change. Both it and the GARCH forecast below take the conditional mean to be zero, so tomorrow’s return is treated as \(N(0, \sigma_{t+1}^2)\). The average daily return is small beside the daily volatility, which is what makes that acceptable here. The update is

\[ \sigma_{t+1}^2 = \lambda\sigma_t^2 + (1-\lambda)y_t^2. \]

Formula Python name Meaning
\(\lambda\) lambda_ Weight retained from the previous variance
\(y_t\) observation Latest observed return
\(\sigma_{t+1}^2\) ewma_variance Next-day variance forecast
\(z_p\) stats.norm.ppf(p) Normal lower-tail quantile

The recursion needs a starting variance. It comes from the first 30 returns, as their average squared value so that the zero-mean assumption is used throughout, after which the remaining returns update it one at a time.

initial_days = 30
ewma_variance = np.mean(estimation_returns[:initial_days] ** 2)
for observation in estimation_returns[initial_days:]:
    ewma_variance = (
        lambda_ * ewma_variance
        + (1 - lambda_) * observation**2
    )

ewma_volatility = np.sqrt(ewma_variance)
ewma_var = -value * stats.norm.ppf(p) * ewma_volatility
print("EWMA VaR: $", round(ewma_var, 2))
EWMA VaR: $ 23.72
NoteCheck your result

The forecast variance, volatility and VaR should all be positive. The recursion runs from the oldest observation to the newest so that the newest return has the most recent effect.

assert ewma_variance > 0
assert ewma_volatility > 0
assert ewma_var > 0

8.6 GARCH VaR

GARCH uses the Week 6 variance equation. As in Week 6, the fit temporarily uses percentage returns for numerical stability. The forecast is then divided by 100**2 to return to decimal variance units.

8.6.1 A function for the GARCH variance forecast

def garch_variance_forecast(returns):
    scaled_returns = 100 * pd.Series(returns)
    model = arch_model(
        scaled_returns,
        mean="Zero",
        vol="GARCH",
        p=1,
        q=1,
        dist="normal",
        rescale=False,
    )
    result = model.fit(disp="off")
    if result.convergence_flag != 0:
        raise RuntimeError("The GARCH optimiser did not converge")
    forecast = result.forecast(horizon=1, reindex=False)
    decimal_variance = forecast.variance.iloc[-1, 0] / 100**2
    return decimal_variance, result
garch_variance, garch_result = garch_variance_forecast(estimation_returns)
garch_volatility = np.sqrt(garch_variance)
garch_var = -value * stats.norm.ppf(p) * garch_volatility

print("GARCH VaR: $", round(garch_var, 2))
GARCH VaR: $ 30.78
NoteCheck your result

The GARCH variance and VaR should be positive. alpha[1] + beta[1] should be below one for this fitted series, as in Week 6.

garch_persistence = (
    garch_result.params["alpha[1]"]
    + garch_result.params["beta[1]"]
)
assert garch_variance > 0
assert garch_var > 0
assert garch_persistence < 1
print("GARCH persistence:", round(garch_persistence, 4))
GARCH persistence: 0.7739

8.7 Comparing the methods

Under a normal conditional distribution, ES has a closed form. The average of the tail beyond \(z_p\) is \(\varphi(z_p)/p\) standard deviations below zero, where \(\varphi\) is the normal density, so the loss is that multiple of the forecast volatility:

normal_es_multiple = stats.norm.pdf(stats.norm.ppf(p)) / p
ewma_es = value * normal_es_multiple * ewma_volatility
garch_es = value * normal_es_multiple * garch_volatility

risk_comparison = pd.DataFrame(
    {
        "VaR": [hs_var, ewma_var, garch_var],
        "ES": [hs_es, ewma_es, garch_es],
    },
    index=["Historical", "EWMA-normal", "GARCH-normal"],
)
risk_comparison.round(2)
VaR ES
Historical 45.61 57.63
EWMA-normal 23.72 27.18
GARCH-normal 30.78 35.27 
fig, ax = plt.subplots()
risk_comparison["VaR"].plot.bar(ax=ax, color=["tab:blue", "tab:orange", "tab:green"])
ax.set_title("One-day JPM VaR estimates")
ax.set_ylabel("Positive loss amount ($)")
ax.tick_params(axis="x", rotation=0)
plt.show()

Historical simulation uses the empirical tail and gives equal weight to every observation in its window. EWMA and GARCH assume a normal conditional return, but react to changing volatility. Different answers therefore reflect different assumptions, not automatically a coding error.

NoteCheck your result

Identify the largest VaR, then give a model-based reason it may be largest. Do not choose a “winner” from one forecast. Week 10 evaluates repeated forecasts.

8.8 Estimation-window choice

Historical simulation can change when the window changes because old crises enter or leave the sample.

window_results = {}
for days in [250, 500, 1000]:
    window_results[days] = historical_var_es(
        y.tail(days), probability=p, portfolio_value=value
    )["VaR"]

pd.Series(window_results, name="Historical VaR").rename_axis("Window days")
Window days
250     45.607909
500     46.015118
1000    45.607909
Name: Historical VaR, dtype: float64

8.9 Recap

  • VaR is a loss threshold, and ES is the average loss beyond that threshold.
  • VaR and ES are stored here as positive loss amounts.
  • Historical simulation uses observed returns without imposing a distribution.
  • Portfolio historical simulation preserves the observed dependence between assets.
  • EWMA and GARCH respond to changing volatility.
  • One forecast cannot establish which risk model is best. That requires backtesting.

8.10 Exercises

8.10.1 Apply the method

  1. Repeat the three JPM VaR calculations for AAPL and compare the results.
  2. Explain in words why historical ES is greater than historical VaR.
  3. Change the portfolio weights to 50% JPM and 50% MCD. Recalculate portfolio historical VaR and ES.

8.10.2 Compare results

  1. Compare historical VaR using 250, 500 and 1,000 observations. Explain why the results change.
  2. Change lambda_ from 0.94 to 0.90 and 0.97. Which EWMA estimate responds most strongly to recent returns?

8.10.3 Extend the analysis

  1. Calculate Student-t GARCH VaR. Fit with dist="t", take nu from the fitted parameters, and use the quantile of a Student-t standardised to variance one, stats.t.ppf(p, df=nu) * np.sqrt((nu - 2) / nu), in place of stats.norm.ppf(p). The scaling is needed because arch standardises its Student-t while stats.t does not, as in Week 6.
  2. Turn the three methods into functions and produce a table for every asset in Returns.csv.