Copula-Tail VaR Engine

Mo Minoneshan

2026

GitHub Demo

Overview

This engine estimates portfolio tail risk by modelling each asset's marginal loss tail with Extreme Value Theory (EVT) and the joint dependence between assets with a Student-t copula. It then computes Value-at-Risk (VaR) and Expected Shortfall (ES) by Monte Carlo simulation of the fitted joint distribution.

The design principle is separation of concerns: the marginal behaviour of each asset (how fat its own tail is) is modelled independently of the dependence structure (how assets crash together). This decoupling — the defining property of a copula — lets us fit heavy tails and tail co-movement that a single joint-Gaussian model cannot represent.

Headline result

On synthetic ground-truth data the engine reduces the absolute 99% Expected Shortfall error by ≈30% (measured 30.43% in the reference demo) relative to a bootstrap historical simulation baseline.

0Executive Summary

Gaussian models understate joint extreme losses because they have thin tails and zero tail dependence. By separating marginals from dependence, this engine captures both fat tails (via EVT/GPD) and tail co-movement (via the Student-t copula) that drive real portfolio drawdowns.

Marginals: Peaks-Over-Threshold EVT — an empirical body below a high quantile spliced to a Generalised Pareto Distribution (GPD) tail above it.
Dependence: Student-t copula fitted by pseudo maximum likelihood, with degrees of freedom $\nu$ optimised and a correlation matrix projected to the nearest positive-definite correlation matrix.
Aggregation: Monte Carlo — sample the copula, invert each marginal, weight into portfolio losses, then read off VaR and ES.
Benchmark: a bootstrap historical simulator provides an apples-to-apples baseline; a large-sample simulation from the known data-generating process provides ground truth.

1Motivation & Problem

In crises, correlations rise and tails fatten simultaneously. Historical simulation — resampling past returns — is simple but has two weaknesses in the tail: it can only ever reproduce losses it has already seen, and with a few hundred observations the 99% quantile is estimated from a mere handful of points, so ES is noisy and typically biased.

A Student-t copula encodes tail dependence — the tendency for assets to hit their extremes together — which a Gaussian copula cannot, and EVT extrapolates smoothly beyond the observed sample using a theoretically-motivated tail shape. The combination yields materially more accurate high-quantile risk estimates.

2Marginal Tails (EVT)

Each asset's losses are modelled by MarginalEVTModel. We work in loss space (losses $= -\text{returns}$) so that the right tail is the region of interest.

2.1Peaks-Over-Threshold

The Pickands–Balkema–de Haan theorem states that for a high threshold $u$, the distribution of excesses $Y = X - u$ given $X > u$ converges to a Generalised Pareto Distribution (GPD):

$$ G_{\xi,\beta}(y) = 1 - \left(1 + \frac{\xi y}{\beta}\right)^{-1/\xi}, \quad y \ge 0,\ \beta > 0. $$

Here $\xi$ is the shape (tail-index) parameter and $\beta$ the scale. $\xi > 0$ gives a heavy, power-law tail — the regime relevant for financial losses. The threshold $u$ is chosen as an empirical quantile (default: the 97th percentile of losses).

2.2Hybrid Body + GPD Tail

Rather than force a parametric form on the whole distribution, the model uses the empirical distribution for the body (below $u$) and the fitted GPD for the tail (above $u$). Writing $p_u = P(X \le u)$, the composite CDF is:

$$ F(x) = \begin{cases} \widehat{F}_{\text{emp}}(x), & x \le u, \\[4pt] p_u + (1 - p_u)\, G_{\xi,\beta}(x - u), & x > u. \end{cases} $$

The inverse (quantile) function ppf mirrors this split, interpolating the empirical body for $u \le p_u$ and inverting the GPD for the tail. This hybrid CDF is exactly the probability-integral transform used to feed the copula.

model = MarginalEVTModel(threshold_quantile=0.97, min_exceedances=25)
model.fit(losses)                 # empirical body + GPD tail
u = model.to_uniform(losses)      # PIT -> pseudo-uniforms in (0, 1)
x = model.ppf(u)                  # inverse transform back to loss space

2.3Fallback Strategy

Robustness by design

If the number of exceedances above $u$ falls below min_exceedances (default 25), fitting a GPD would be unstable. The model automatically reverts to a purely empirical distribution for that asset, so the pipeline degrades gracefully on short or thin-tailed histories instead of producing an unreliable parametric fit.

3Dependence (t-Copula)

A copula $C$ couples uniform margins into a joint distribution. By Sklar's theorem any joint CDF $H$ decomposes as $H(x_1,\dots,x_d) = C\big(F_1(x_1),\dots,F_d(x_d)\big)$. The Student-t copula produces symmetric tail dependence that strengthens as the degrees of freedom $\nu$ decrease.

3.1PIT & Pseudo-Uniforms

Each asset's losses are mapped to $(0,1)$ via its own fitted hybrid CDF (the probability integral transform). Stacking these columns yields a matrix of pseudo-uniform observations $\{u_{t,j}\}$ on which the copula is estimated — dependence is thus learned free of the marginals.

3.2Likelihood & Degrees of Freedom

The copula density is evaluated on the t-quantiles $z_{t,j} = t_\nu^{-1}(u_{t,j})$. The Student-t copula log-density for a $d$-dimensional observation with correlation $R$ and DoF $\nu$ is:

$$ \ln c(u) = \ln\!\frac{\Gamma\!\left(\frac{\nu+d}{2}\right)\Gamma\!\left(\frac{\nu}{2}\right)^{d-1}} {\Gamma\!\left(\frac{\nu+1}{2}\right)^{d}\,|R|^{1/2}} - \frac{\nu+d}{2}\ln\!\left(1 + \frac{z^\top R^{-1} z}{\nu}\right) + \sum_{j=1}^{d}\frac{\nu+1}{2}\ln\!\left(1 + \frac{z_j^2}{\nu}\right). $$

Fitting proceeds by an alternating scheme: initialise $R$ from the (Gaussian-score) correlation, optimise $\nu$ on a bounded interval $[2.1, 30]$ with scipy.optimize.minimize_scalar, re-estimate $R$ from the t-scores, and re-optimise $\nu$. The lower bound 2.1 keeps the variance finite; the upper bound approaches the Gaussian limit.

copula = TCopula(df_bounds=(2.1, 30.0), regularization=1e-6)
copula.fit(uniforms_matrix)       # pseudo-MLE for R and nu
print(copula.params.df)           # fitted degrees of freedom
print(copula.log_likelihood)      # goodness-of-fit diagnostic

3.3Correlation PD Projection

Sample correlation matrices can be indefinite or singular, especially as the asset count grows. Before use, the matrix is projected to the nearest valid correlation matrix: symmetrise, clip eigenvalues to a small floor $\varepsilon$, reconstruct, and renormalise to unit diagonal. This guarantees a Cholesky factor exists for sampling.

$$ R \;\leftarrow\; D^{-1/2}\,\big(Q \max(\Lambda, \varepsilon) Q^\top\big)\,D^{-1/2}, $$

where $Q\Lambda Q^\top$ is the eigendecomposition of the symmetrised matrix and $D$ its diagonal.

3.4Sampling

To draw from the fitted copula we use the standard multivariate-t construction: correlate standard normals with the Cholesky factor $L$ of $R$, divide by $\sqrt{\chi^2_\nu/\nu}$, and map through the univariate t-CDF to return to uniforms.

$$ z = \frac{L\,\varepsilon}{\sqrt{s/\nu}},\quad \varepsilon \sim \mathcal N(0, I_d),\ s \sim \chi^2_\nu,\qquad u = t_\nu(z). $$

4Risk Measures

For a portfolio loss $L$ at confidence level $\alpha$:

4.1Value-at-Risk

$$ \mathrm{VaR}_\alpha(L) = \inf\{\ell : P(L \le \ell) \ge \alpha\}. $$

In code this is the empirical $\alpha$-quantile of the simulated loss vector, np.quantile(losses, alpha).

4.2Expected Shortfall

ES (a.k.a. Conditional VaR) is the average loss given we are beyond VaR — a coherent risk measure that, unlike VaR, is sensitive to how bad the tail actually is:

$$ \mathrm{ES}_\alpha(L) = \mathbb{E}\!\left[L \mid L \ge \mathrm{VaR}_\alpha(L)\right]. $$

4.3Monte Carlo Aggregation

The engine draws $N$ copula samples (default 120k in the demo), inverts each asset's marginal ppf, forms weighted portfolio losses $L = \sum_j w_j\,\ell_j$, and computes VaR and ES on that vector. Simulated paths are cached and reused unless a larger draw or an explicit refresh is requested.

engine = TailRiskEngine(threshold_quantile=0.97)
engine.fit(returns, weights=weights)
result = engine.expected_shortfall(alpha=0.99, n_paths=120_000, random_state=42)
print(result.var, result.expected_shortfall)

5Architecture & Data Flow

The pipeline is a directed flow from raw returns to risk numbers, with the historical baseline branching off the same inputs:

raw returns
   │
   ▼
data prep (drop NaNs, normalise weights)
   │
   ▼
marginal EVT fit  ──▶  pseudo-uniforms (PIT)  ──▶  t-copula fit
                                                        │
                                                        ▼
                                              sample copula uniforms
                                                        │
                                                        ▼
                                              EVT inverse transforms
                                                        │
                                                        ▼
                                    weight & aggregate → portfolio losses
                                                        │
                            ┌───────────────────────────┴───────────┐
                            ▼                                        ▼
                     model VaR / ES                       historical bootstrap ES

Package layout:

src/copula_tail_var/
    data.py         # ingestion + synthetic t-copula/GPD generators
    evt.py          # MarginalEVTModel — POT hybrid body + GPD tail
    copula.py       # TCopula — pseudo-MLE fit, PD projection, sampling
    engine.py       # TailRiskEngine — orchestration, simulate, ES
    historical.py   # BootstrapHistoricalSimulator — baseline
    metrics.py      # value_at_risk, expected_shortfall, error helpers
examples/run_demo.py
tests/ (test_evt, test_copula, test_engine)

6Implementation

End-to-end programmatic workflow, from returns to a benchmarked ES comparison:

import numpy as np
import pandas as pd
from copula_tail_var import TailRiskEngine

returns = pd.read_csv("daily_returns.csv", index_col=0)
weights = np.array([0.4, 0.35, 0.25])

engine = TailRiskEngine(threshold_quantile=0.97)
engine.fit(returns, weights=weights)

result = engine.expected_shortfall(alpha=0.99, n_paths=120_000, random_state=42)
print(f"99% VaR: {result.var:.6f}")
print(f"99% ES:  {result.expected_shortfall:.6f}")

# Compare to a bootstrap historical baseline (and known truth if available)
summary = engine.compare_to_historical(alpha=0.99, true_es=0.159)
print(summary["error_reduction"])   # relative ES-error reduction

Key design decisions baked into the code:

Weights are normalised to sum to one and default to equal weighting when omitted.
Every stochastic routine accepts random_state, making runs and tests fully reproducible.
Simulated paths are cached; requesting fewer paths slices the cache, more paths (or refresh=True) triggers a fresh draw.

7Synthetic Data & Benchmark

Because real tail data is scarce, the engine ships a generator that samples returns from a known t-copula with GPD-tailed marginals (generate_evt_copula_dataset). This gives an exact data-generating process, so a very large simulation from the true parameters yields a ground-truth ES against which both the model and the historical baseline can be scored.

from copula_tail_var.data import MarginalSpec, generate_evt_copula_dataset

corr = np.array([[1.0, 0.7, 0.6],
                 [0.7, 1.0, 0.65],
                 [0.6, 0.65, 1.0]])
specs = [
    MarginalSpec(mu=0.0,   sigma=0.019, threshold=0.030, gpd_shape=0.38, gpd_scale=0.026),
    MarginalSpec(mu=0.0,   sigma=0.022, threshold=0.033, gpd_shape=0.40, gpd_scale=0.028),
    MarginalSpec(mu=-0.001,sigma=0.018, threshold=0.028, gpd_shape=0.34, gpd_scale=0.024),
]
returns = generate_evt_copula_dataset(n_obs=600, corr=corr, copula_df=3.2,
                                      marginals=specs, random_state=42)

8Results

Reference demo (examples/run_demo.py): trained on 600 synthetic observations, 120k model simulations, 40k bootstrap resamples, with ground truth from a ~1M-path simulation of the true process.

Reference demo output

Target confidence level: 99%
True ES (benchmark):  0.158991
Model ES:             0.216317  (VaR: 0.119178)
Historical ES:        0.241389  (VaR: 0.094927)
Model absolute error:      0.057326
Historical absolute error: 0.082399
Relative error reduction vs historical: 30.43%

The copula/EVT estimate sits closer to the ground-truth ES than the historical baseline, cutting absolute error by 30.43%.

9Validation

The pytest suite exercises each layer of the stack:

EVT: monotonicity and inverse consistency — ppf and cdf round-trip, and the fitted tail behaves correctly at the threshold.
Copula: the fitted t-copula log-likelihood improves over a naive initial guess, and sampled uniforms recover the target dependence.
Engine: an end-to-end regression test asserts a $\ge 30\%$ ES-error reduction versus historical simulation on controlled synthetic data.

pytest -q

For real deployments, VaR exceptions would additionally be backtested with Kupiec's unconditional-coverage and Christoffersen's independence tests.

10Limitations

Assumptions & caveats

EVT threshold selection and the copula degrees of freedom are sensitive to the sample; estimates can be unstable for short histories. The Student-t copula imposes symmetric tail dependence, so it cannot distinguish downside-only clustering. Marginals are fitted independently per asset (returns are treated as i.i.d. through time — no volatility clustering / GARCH filter), and the reported 30% figure is measured on a controlled synthetic process, not a specific live portfolio.

AAPI Reference

class TailRiskEngine:
    def __init__(self, threshold_quantile: float = 0.97, min_exceedances: int = 25): ...
    def fit(self, returns: pd.DataFrame, weights: np.ndarray | None = None) -> "TailRiskEngine": ...
    def simulate(self, n_paths: int = 20000, random_state: int | None = None,
                 refresh: bool = False) -> np.ndarray: ...
    def expected_shortfall(self, alpha: float = 0.99, n_paths: int = 20000,
                           random_state: int | None = None) -> TailRiskResult: ...
    def compare_to_historical(self, alpha: float = 0.99, n_paths: int = 20000,
                              n_bootstrap: int = 20000, random_state: int | None = None,
                              true_es: float | None = None) -> dict: ...

class MarginalEVTModel:
    def __init__(self, threshold_quantile: float = 0.97, min_exceedances: int = 25): ...
    def fit(self, losses: np.ndarray) -> "MarginalEVTModel": ...
    def cdf(self, values: np.ndarray) -> np.ndarray: ...
    def ppf(self, uniforms: np.ndarray) -> np.ndarray: ...
    def to_uniform(self, losses: np.ndarray) -> np.ndarray: ...

class TCopula:
    def __init__(self, df_bounds=(2.1, 30.0), regularization: float = 1e-6): ...
    def fit(self, uniforms: np.ndarray) -> "TCopula": ...
    def sample(self, n_samples: int, random_state: int | None = None) -> np.ndarray: ...
    def log_density(self, uniforms: np.ndarray) -> float: ...

class BootstrapHistoricalSimulator:
    def __init__(self, returns: pd.DataFrame, weights: np.ndarray | None = None): ...
    def expected_shortfall(self, alpha: float, n_paths: int = 10000,
                           random_state: int | None = None) -> HistoricalESResult: ...

# metrics
def value_at_risk(losses: np.ndarray, alpha: float) -> float: ...
def expected_shortfall(losses: np.ndarray, alpha: float) -> float: ...
def absolute_error(estimate: float, truth: float) -> float: ...
def relative_improvement(baseline_error: float, new_error: float) -> float: ...

BGlossary

Tail risk — risk of extreme losses in the far end of the return distribution.
VaR — the loss quantile not exceeded with probability $\alpha$.
Expected Shortfall (ES / CVaR) — mean loss conditional on exceeding VaR; a coherent risk measure.
Copula — a function linking uniform margins into a joint distribution (Sklar's theorem).
Student-t copula — copula with symmetric tail dependence controlled by DoF $\nu$.
Tail dependence — probability that one asset is extreme given another is extreme.
EVT — Extreme Value Theory; statistics of maxima / threshold exceedances.
POT — Peaks-Over-Threshold; models excesses above a high threshold.
GPD — Generalised Pareto Distribution; limiting law for POT excesses.
Shape $\xi$ — GPD tail index; $\xi>0$ is heavy-tailed.
PIT — Probability Integral Transform; maps data to $(0,1)$ via its CDF.
Pseudo-MLE — likelihood maximisation on PIT pseudo-observations.
Monte Carlo — estimation by repeated random sampling.
Bootstrap — resampling observed data with replacement.

CReferences

Balkema, A. & de Haan, L. (1974). Residual life time at great age.
Pickands, J. (1975). Statistical inference using extreme order statistics.
Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges.
McNeil, Frey & Embrechts (2015). Quantitative Risk Management.
Artzner, Delbaen, Eber & Heath (1999). Coherent Measures of Risk.
Kupiec (1995); Christoffersen (1998). VaR backtesting.