Week 10: Backtesting Risk Models

Version 5.0 - August 2026

Author

Affiliation

Jon Danielsson

London School of Economics

10 Week 10

A risk model produces a forecast before the return is observed. Backtesting keeps that order intact, compares each forecast with the return that followed, and records whether the loss exceeded VaR.

If a 1% VaR model is correctly calibrated, violations should occur on about one day in a hundred. This week compares historical simulation, EWMA and GARCH over many forecast days, then applies a Bernoulli coverage test to their violation counts.

NoteBefore you start

Open the main course folder and check that Returns.csv, created in Week 2, is present. You will use historical simulation and EWMA from Week 8 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 have three aligned VaR forecast series, a violation count and rate for each model, and a coverage-test result that you can explain in words.

10.1 The plan for this week

1. Reserve an initial estimation period
2. Produce historical, EWMA and GARCH VaR forecasts
3. Record violations using one sign convention
4. Compare violation rates
5. Test whether each rate is consistent with 1%

10.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)

p = 0.01
lambda_ = 0.94
value = 1
T = 5000
estimation_days = 1000
y = Returns["MCD"].dropna().tail(T)
y_values = y.to_numpy()

CautionUnits, signs and timing

• Returns are one-day decimals.
• VaR is stored as a positive loss on a portfolio worth value = 1.
• Therefore a violation is exactly return < -VaR.
• The forecast at position t may use observations only up to t - 1.
• GARCH volatility is a daily decimal standard deviation.

NoteCheck your result

There should be 5,000 MCD returns. The first 1,000 form the initial estimation period, leaving 4,000 forecasts to evaluate.

assert len(y) == T
testing_days = T - estimation_days
print("estimation observations:", estimation_days)
print("testing observations:", testing_days)
print(y.index.min(), "to", y.index.max())

estimation observations: 1000
testing observations: 4000
2005-11-07 00:00:00 to 2025-09-23 00:00:00

10.3 Storing aligned forecasts

Each row refers to one date. Forecast columns begin with missing values because the models need earlier data before they can make a forecast.

VaR = pd.DataFrame(index=y.index)
VaR["return"] = y_values
VaR["HS"] = np.nan
VaR["EWMA"] = np.nan
VaR["GARCH"] = np.nan

10.4 Rolling historical simulation

At each forecast date, historical simulation uses the preceding 1,000 returns. The slice ends at t, so it does not include the return being forecast.

hs_forecasts = np.full(T, np.nan)

for t in range(estimation_days, T):
    window = y_values[t - estimation_days : t]
    threshold = np.quantile(window, p, method="lower")
    hs_forecasts[t] = -value * threshold

VaR["HS"] = hs_forecasts

NoteCheck your result

The first forecast is at position 1,000. Exactly 1,000 earlier observations feed it, and the realised return at position 1,000 is not in that window.

assert np.isnan(hs_forecasts[:estimation_days]).all()
assert np.isfinite(hs_forecasts[estimation_days:]).all()
assert (hs_forecasts[estimation_days:] > 0).all()

10.5 EWMA

The one-step-ahead EWMA variance is

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

The subscript matters. The variance forecast for date t uses return t - 1, not the still-unknown return at t.

The recursion starts from the average squared return over the first 30 observations, placed at position 30 rather than at position 0, so that no forecast is made from returns that had not yet occurred. As in Week 8, the starting value squares the returns rather than taking their sample variance, which keeps the zero-mean assumption of the recursion.

initial_days = 30
ewma_variance = np.full(T, np.nan)
ewma_variance[initial_days] = np.mean(y_values[:initial_days] ** 2)

for t in range(initial_days + 1, T):
    ewma_variance[t] = (
        lambda_ * ewma_variance[t - 1]
        + (1 - lambda_) * y_values[t - 1] ** 2
    )

ewma_forecasts = -value * stats.norm.ppf(p) * np.sqrt(ewma_variance)
ewma_forecasts[:estimation_days] = np.nan
VaR["EWMA"] = ewma_forecasts

NoteCheck your result

The 30-observation starting value is burned in long before evaluation begins at position 1,000. All 4,000 evaluated EWMA forecasts should be finite and positive.

assert np.isfinite(ewma_forecasts[estimation_days:]).all()
assert (ewma_forecasts[estimation_days:] > 0).all()

10.6 Fixed-parameter recursive GARCH

We estimate one normal GARCH(1,1) from the initial 1,000 observations. During the testing period its parameters remain fixed, while its conditional variance is updated after each newly observed return. This cleanly separates initial parameter estimation from out-of-sample recursive forecasting.

The variance recursion is

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

Formula	Python name	Meaning
\(\omega\)	omega	Baseline variance contribution
\(\alpha\)	alpha	Response to the latest squared return
\(\beta\)	beta	Carry-over from the previous variance
\(\sigma_t^2\)	current_variance	Forecast variance for date t

As in Weeks 6 and 8, estimation temporarily uses percentage returns for numerical stability. omega and conditional variance are converted back to decimal units.

initial_returns = 100 * y.iloc[:estimation_days]
garch_model = arch_model(
    initial_returns,
    mean="Zero",
    vol="GARCH",
    p=1,
    q=1,
    dist="normal",
    rescale=False,
)
garch_fit = garch_model.fit(disp="off")
if garch_fit.convergence_flag != 0:
    raise RuntimeError("The initial GARCH fit did not converge")

omega = garch_fit.params["omega"] / 100**2
alpha = garch_fit.params["alpha[1]"]
beta = garch_fit.params["beta[1]"]
current_variance = (garch_fit.conditional_volatility.iloc[-1] / 100) ** 2

print("omega:", omega)
print("alpha:", round(alpha, 4))
print("beta:", round(beta, 4))
print("persistence:", round(alpha + beta, 4))

omega: 2.755514595534986e-06
alpha: 0.0624
beta: 0.9253
persistence: 0.9877

garch_forecasts = np.full(T, np.nan)

for t in range(estimation_days, T):
    current_variance = (
        omega
        + alpha * y_values[t - 1] ** 2
        + beta * current_variance
    )
    garch_forecasts[t] = (
        -value * stats.norm.ppf(p) * np.sqrt(current_variance)
    )

VaR["GARCH"] = garch_forecasts

NoteCheck your result

alpha + beta should be below one for this fit. All evaluated GARCH variances and VaRs should be positive. The return at t is used only after its forecast has been made.

assert alpha >= 0 and beta >= 0
assert alpha + beta < 1
assert np.isfinite(garch_forecasts[estimation_days:]).all()
assert (garch_forecasts[estimation_days:] > 0).all()

10.7 Comparing the forecast paths

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(VaR.index, VaR["return"], color="grey", linewidth=0.5, label="Return")
ax.plot(VaR.index, -VaR["HS"], linewidth=1, label="-HS VaR")
ax.plot(VaR.index, -VaR["EWMA"], linewidth=1, label="-EWMA VaR")
ax.plot(VaR.index, -VaR["GARCH"], linewidth=1, label="-GARCH VaR")
ax.set_title("MCD returns and one-day VaR thresholds")
ax.set_xlabel("Date")
ax.set_ylabel("Return / threshold")
ax.legend(frameon=False)
plt.show()

Historical simulation changes when observations enter or leave its window. EWMA and GARCH change recursively and react more directly to recent squared returns.

Two things are being compared here beyond the three variance rules. EWMA and GARCH convert a volatility forecast into a VaR through stats.norm.ppf(p), so each carries a normal-tail assumption that historical simulation does not. A rejection below may be a failure of that assumption rather than of the variance forecast. And historical simulation refreshes its estimation window every day, while the GARCH parameters are fixed after the initial 1,000 observations. What is evaluated is therefore three complete forecasting schemes, including how each updates what it knows, rather than three volatility models in isolation.

10.8 Violations

A violation occurs when the realised return is more negative than the negative VaR threshold. Because value = 1, the rule is return < -VaR.

evaluation = VaR.iloc[estimation_days:].copy()
violations = pd.DataFrame(index=evaluation.index)

for method in ["HS", "EWMA", "GARCH"]:
    violations[method] = (
        evaluation["return"] < -evaluation[method]
    ).astype(int)

violation_summary = pd.DataFrame(
    {
        "Violations": violations.sum(),
        "Rate": violations.mean(),
        "Rate / expected rate": violations.mean() / p,
    }
)
violation_summary.round(4)

	Violations	Rate	Rate / expected rate
HS	38	0.0095	0.950
EWMA	78	0.0195	1.950
GARCH	51	0.0127	1.275

NoteCheck your result

With 4,000 forecast days and p = 0.01, about 40 violations are expected. A ratio above one means too many violations and underestimated risk, and a ratio below one means too few violations and overestimated risk.

The count is only part of the picture. Plotting the dates on which each method was violated shows whether they arrive evenly or in bursts:

fig, ax = plt.subplots(figsize=(9, 3))
for position, method in enumerate(["HS", "EWMA", "GARCH"]):
    violation_dates = violations.index[violations[method] == 1]
    ax.plot(
        violation_dates,
        np.full(len(violation_dates), position),
        "|",
        markersize=12,
    )
ax.set_yticks([0, 1, 2])
ax.set_yticklabels(["HS", "EWMA", "GARCH"])
ax.set_title("Dates of VaR violations")
ax.set_xlabel("Date")
fig.tight_layout()
plt.show()

Violations that cluster in a few turbulent stretches are worse than the same number spread evenly, because they arrive when losses are already accumulating. The test below does not see that difference.

10.9 Bernoulli coverage test

The hypotheses are:

• Null hypothesis: the true violation probability equals p = 0.01.
• Alternative hypothesis: the true violation probability differs from p.

The function below compares the Bernoulli log-likelihood under those two hypotheses. This is the standard likelihood-ratio test of unconditional coverage, usually named after Kupiec. It uses nothing but the total number of violations, so it cannot detect clustering, and its chi-squared calibration assumes that violations are independent from day to day under the null. A model whose violations arrive in bursts can therefore pass it while breaking the assumption on which its p-value rests.

10.9.1 A function for the coverage test

def bernoulli_coverage_test(probability, violation_series):
    number_of_days = len(violation_series)
    number_of_violations = int(violation_series.sum())
    observed_probability = number_of_violations / number_of_days

    loglik_null = stats.binom.logpmf(
        number_of_violations,
        number_of_days,
        probability,
    )
    loglik_alternative = stats.binom.logpmf(
        number_of_violations,
        number_of_days,
        observed_probability,
    )
    statistic = 2 * (loglik_alternative - loglik_null)
    p_value = stats.chi2.sf(statistic, df=1)
    return {
        "Observed rate": observed_probability,
        "Test statistic": statistic,
        "p-value": p_value,
    }

coverage_results = pd.DataFrame(
    {
        method: bernoulli_coverage_test(p, violations[method])
        for method in ["HS", "EWMA", "GARCH"]
    }
).T
coverage_results.round(4)

	Observed rate	Test statistic	p-value
HS	0.0095	0.1027	0.7486
EWMA	0.0195	28.5472	0.0000
GARCH	0.0127	2.8111	0.0936

At the 5% significance level, reject the null when the p-value is below 0.05. Not rejecting does not prove the model is correct. This test considers only the overall violation rate, not whether violations arrive in clusters.

coverage_results["Decision at 5%"] = np.where(
    coverage_results["p-value"] < 0.05,
    "Reject equal coverage",
    "Do not reject equal coverage",
)
coverage_results

	Observed rate	Test statistic	p-value	Decision at 5%
HS	0.00950	0.102720	7.485901e-01	Do not reject equal coverage
EWMA	0.01950	28.547201	9.144194e-08	Reject equal coverage
GARCH	0.01275	2.811094	9.361451e-02	Do not reject equal coverage

NoteCheck your result

For each model, report its observed rate, p-value and decision. Then distinguish the statistical statement (“reject” or “do not reject”) from a broader claim that the model is good or bad.

10.10 Repeated GARCH re-estimation

The core GARCH backtest estimates parameters once and updates variance recursively. An alternative is to re-estimate omega, alpha and beta on a moving window before every forecast. That allows parameters themselves to change, but it is much slower and introduces more estimation variation. It is an advanced model choice, not a required step for understanding backtesting.

10.11 Recap

• A forecast must use only information available before the realised return.
• Historical simulation uses a rolling window, while EWMA and fixed-parameter GARCH update recursively.
• VaR is a positive loss and the violation rule here is return < -VaR.
• A correct 1% VaR should produce approximately 1% violations over a long sample.
• The Bernoulli coverage test asks whether the overall violation probability equals 1%. It does not test violation independence, although its p-value assumes it.
• EWMA and GARCH VaR combine a volatility forecast with a normal-tail assumption, and a rejection does not say which of the two failed.

10.12 Exercises

10.12.1 Apply the method

1. Repeat the backtest for JPM and compare its violation table with MCD’s.
2. For each MCD method, explain whether the violation ratio indicates too much or too little forecast risk.
3. State the null and alternative hypotheses of the coverage test in words.

10.12.2 Compare results

4. Repeat historical simulation with 500 estimation days. Explain what changes and why.
5. Repeat EWMA with lambda_ = 0.90 and 0.97. Compare violation rates and explain the responsiveness trade-off.

10.12.3 Extend the analysis

6. Re-estimate GARCH on a rolling window every 20 forecast days, holding the parameters fixed in between, and compare the result with the fixed-parameter recursion. Re-estimating before every one of the 4,000 forecasts is the same idea at two hundred times the cost, and some of those fits will not converge. Record the runtime as well as the violation results.
7. Test whether violations are independent through time. Count the four transitions between consecutive days, \(n_{00}\), \(n_{01}\), \(n_{10}\) and \(n_{11}\), where \(n_{01}\) is the number of days that follow a quiet day with a violation. Under independence the probability of a violation is the same after a quiet day as after a violation, so compare the log-likelihood of one common probability against the log-likelihood of two separate ones, exactly as the coverage test compares one probability against another. The statistic is again chi-squared with one degree of freedom. Explain why correct overall coverage can coexist with clustered violations.