library(xts)
library(tsmarch)
library(tsgarch)
library(MASS)
source("common/functions.r", chdir = TRUE)24 Portfolio risk
Portfolio risk is not the sum of individual asset risks. Correlations between assets drive the diversification benefit, and these correlations change over time. While the VaR methods so far applied to a single return series, practical risk management requires understanding how individual asset risks combine and contribute to overall portfolio risk.
We start with three standard VaR methods. The parametric and Monte Carlo approaches use the sample covariance as a benchmark, historical simulation works directly from the empirical returns, and an EWMA covariance variant follows. We then show how replacing the sample covariance with DCC covariance forecasts from Chapter 22 captures time-varying correlations.
24.1 Data and libraries
A multivariate normal generator for the Monte Carlo method, and tsmarch for the DCC forecast.
import numpy as np
from scipy import stats
import sys
sys.path.insert(0, 'common')
from functions import ProcessRawData, EmpiricalESinclude("common/functions.jl");data = ProcessRawData()
Return = data$Return
dates = Return$dateWe work with a two-asset portfolio of JPM and Intel:
assets = c("JPM", "INTC")
y = as.matrix(Return[, assets])
y_xts = xts(y, order.by = dates)data_py = ProcessRawData()
Return_py = data_py['Return']
assets_py = ['JPM', 'INTC']
y_py = Return_py[assets_py].valuesdata_jl = ProcessRawData();
Return_jl = data_jl["Return"];
assets_jl = ["JPM", "INTC"];
y_jl = Matrix(Return_jl[:, Symbol.(assets_jl)]);24.2 Portfolio returns
For a portfolio with weight vector \(\weights\), the portfolio return at time \(t\) is:
\[\CompoundReturns_{\Portfolio,t} = \weights^\top \CompoundReturns_t = \sum_{i=1}^{\NumberAssets} \weight_i \CompoundReturns_{i,t}\]
where \(\CompoundReturns_t\) is the vector of individual asset returns and \(\NumberAssets\) is the number of assets. This is exact for simple returns. For log returns it is a first-order approximation accurate for daily data.
We use equal weights for illustration:
w = c(0.5, 0.5)
names(w) = assets
# Portfolio returns
r_portfolio = y %*% w
head(r_portfolio) [,1]
2 0.005217571
3 0.056833446
4 0.002702257
5 -0.039573955
6 0.021066306
7 0.006579230w_py = np.array([0.5, 0.5])
r_portfolio_py = y_py @ w_py
print(r_portfolio_py[:6])[ 0.00521757 0.05683345 0.00270226 -0.03957396 0.02106631 0.00657923]w_jl = [0.5, 0.5];
r_portfolio_jl = y_jl * w_jl;
println(r_portfolio_jl[1:6])Union{Missing, Float64}[0.005217571194306103, 0.05683344608209895, 0.002702256936126801, -0.03957395509423889, 0.021066305570529575, 0.006579229777937012]# Compare individual and portfolio volatility
cat("JPM annualised volatility: ", round(sd(y[, "JPM"]) * sqrt(252) * 100, 1), "%\n")
cat("INTC annualised volatility: ", round(sd(y[, "INTC"]) * sqrt(252) * 100, 1), "%\n")
cat("Portfolio annualised volatility:", round(sd(r_portfolio) * sqrt(252) * 100, 1), "%\n")
cat("Correlation: ", round(cor(y[, 1], y[, 2]), 3), "\n")JPM annualised volatility: 36.3 %
INTC annualised volatility: 31.4 %
Portfolio annualised volatility: 29.1 %
Correlation: 0.478print(f"JPM annualised volatility: {np.std(y_py[:, 0], ddof=1) * np.sqrt(252) * 100:.1f}%")
print(f"INTC annualised volatility: {np.std(y_py[:, 1], ddof=1) * np.sqrt(252) * 100:.1f}%")
print(f"Portfolio annualised volatility: {np.std(r_portfolio_py, ddof=1) * np.sqrt(252) * 100:.1f}%")
print(f"Correlation: {np.corrcoef(y_py[:, 0], y_py[:, 1])[0, 1]:.3f}")JPM annualised volatility: 36.3%
INTC annualised volatility: 31.4%
Portfolio annualised volatility: 29.1%
Correlation: 0.478using Printf, Statistics
@printf("JPM annualised volatility: %.1f%%\n", std(y_jl[:, 1]) * sqrt(252) * 100)
@printf("INTC annualised volatility: %.1f%%\n", std(y_jl[:, 2]) * sqrt(252) * 100)
@printf("Portfolio annualised volatility: %.1f%%\n", std(r_portfolio_jl) * sqrt(252) * 100)
@printf("Correlation: %.3f\n", cor(y_jl[:, 1], y_jl[:, 2]))JPM annualised volatility: 36.3%
INTC annualised volatility: 31.4%
Portfolio annualised volatility: 29.1%
Correlation: 0.478In this example, the portfolio volatility is lower than the weighted average of individual volatilities due to diversification.
24.3 Portfolio VaR and ES: three approaches
The parametric and Monte Carlo methods begin with a static sample covariance. Historical simulation uses the empirical portfolio returns directly. Section 24.3.4 replaces the sample covariance with an EWMA covariance, and Section 24.4 with DCC forecasts, to capture time-varying correlations.
The parametric method (also called variance-covariance) assumes returns are multivariate normal. Historical simulation uses the empirical distribution of portfolio returns instead, with no distributional assumption. Monte Carlo simulates from the estimated distribution to build the same empirical quantile.
All three use the same parameters:
par = list()
par$probability = 0.05 # 5% VaR
par$value = 1000 # Portfolio value
par$WE = 500 # Estimation windowpar_py = {}
par_py['probability'] = 0.05
par_py['value'] = 1000
par_py['WE'] = 500par_jl = Dict{String, Any}();
par_jl["probability"] = 0.05;
par_jl["value"] = 1000;
par_jl["WE"] = 500;24.3.1 Parametric VaR and ES
The parametric approach uses the covariance matrix to compute portfolio variance:
\[\Vol_\Portfolio^2 = \weights^\top \CovMatrix \weights\]
For normally distributed returns with zero mean, VaR as a positive loss magnitude is:
\[\VaR(\probability) = -\NormalQuantile(\probability)\,\Vol_\Portfolio \times \PortfolioValue\]
where \(\NormalQuantile(\probability)\) is the normal quantile at probability \(\probability\). We assume zero-mean returns, which is reasonable for daily data where the mean is small relative to volatility.
ParametricVaR = function(w, Sigma, p = 0.05, value = 1000) {
sigma_p = sqrt(t(w) %*% Sigma %*% w)
VaR = -qnorm(p) * sigma_p * value
ES = (dnorm(qnorm(p)) / p) * sigma_p * value
return(list(VaR = as.numeric(VaR), ES = as.numeric(ES), sigma = as.numeric(sigma_p)))
}def ParametricVaR(w, Sigma, p=0.05, value=1000):
sigma_p = np.sqrt(w @ Sigma @ w)
VaR = -stats.norm.ppf(p) * sigma_p * value
ES = (stats.norm.pdf(stats.norm.ppf(p)) / p) * sigma_p * value
return {'VaR': VaR, 'ES': ES, 'sigma': sigma_p}using Distributions
function ParametricVaR(w, Sigma, p=0.05, value=1000)
sigma_p = sqrt(w' * Sigma * w)
VaR = -quantile(Normal(), p) * sigma_p * value
ES = (pdf(Normal(), quantile(Normal(), p)) / p) * sigma_p * value
return Dict("VaR" => VaR, "ES" => ES, "sigma" => sigma_p)
endUsing the sample covariance matrix:
Sigma_sample = cov(tail(y, par$WE))
res_param = ParametricVaR(w, Sigma_sample, par$probability, par$value)
cat("Parametric VaR:", round(res_param$VaR, 2), "\n")
cat("Parametric ES: ", round(res_param$ES, 2), "\n")Parametric VaR: 26.95
Parametric ES: 33.8Sigma_sample_py = np.cov(y_py[-par_py['WE']:], rowvar=False)
res_param_py = ParametricVaR(w_py, Sigma_sample_py, par_py['probability'], par_py['value'])
print(f"Parametric VaR: {res_param_py['VaR']:.2f}")
print(f"Parametric ES: {res_param_py['ES']:.2f}")Parametric VaR: 26.95
Parametric ES: 33.80using Statistics, Printf
Sigma_sample_jl = cov(y_jl[end-par_jl["WE"]+1:end, :]);
res_param_jl = ParametricVaR(w_jl, Sigma_sample_jl, par_jl["probability"], par_jl["value"]);
@printf("Parametric VaR: %.2f\n", res_param_jl["VaR"])
@printf("Parametric ES: %.2f\n", res_param_jl["ES"])Parametric VaR: 26.95
Parametric ES: 33.80The parametric figures are the benchmark the other three methods are measured against. Everything about them follows from the covariance matrix and the normal quantile, so they inherit whatever that assumption gets wrong.
24.3.2 Historical simulation
Historical simulation computes portfolio returns for each historical day and takes the empirical quantile, using the rank convention set out in Section 23.5. Here \(\probability \EstWindow = 0.05 \times 500 = 25\), so ceiling and truncation give the same rank.
HistoricalVaR = function(y, w, p = 0.05, value = 1000, WE = 500) {
y_window = tail(y, WE)
r_portfolio = y_window %*% w
r_sorted = sort(r_portfolio)
idx = ceiling(p * WE)
VaR = -r_sorted[idx] * value
ES = -EmpiricalES(r_sorted, p) * value
return(list(VaR = VaR, ES = ES))
}
res_hs = HistoricalVaR(y, w, par$probability, par$value, par$WE)
cat("Historical VaR:", round(res_hs$VaR, 2), "\n")
cat("Historical ES: ", round(res_hs$ES, 2), "\n")Historical VaR: 25.89
Historical ES: 35.91def HistoricalVaR(y, w, p=0.05, value=1000, WE=500):
y_window = y[-WE:]
r_portfolio = y_window @ w
r_sorted = np.sort(r_portfolio)
idx = int(np.ceil(p * WE))
VaR = -r_sorted[idx - 1] * value # -1 for 0-based indexing
ES = -EmpiricalES(r_sorted, p) * value
return {'VaR': VaR, 'ES': ES}
res_hs_py = HistoricalVaR(y_py, w_py, par_py['probability'], par_py['value'], par_py['WE'])
print(f"Historical VaR: {res_hs_py['VaR']:.2f}")
print(f"Historical ES: {res_hs_py['ES']:.2f}")Historical VaR: 25.89
Historical ES: 35.91using Statistics, Printf
function HistoricalVaR(y, w, p=0.05, value=1000, WE=500)
y_window = y[end-WE+1:end, :]
r_portfolio = y_window * w
r_sorted = sort(r_portfolio)
idx = ceil(Int, p * WE)
VaR = -r_sorted[idx] * value
ES = -EmpiricalES(r_sorted, p) * value
return Dict("VaR" => VaR, "ES" => ES)
end
res_hs_jl = HistoricalVaR(y_jl, w_jl, par_jl["probability"], par_jl["value"], par_jl["WE"]);
@printf("Historical VaR: %.2f\n", res_hs_jl["VaR"])
@printf("Historical ES: %.2f\n", res_hs_jl["ES"])Historical VaR: 25.89
Historical ES: 35.91Historical VaR and ES are deterministic given the data and window, so all three languages return the same values.
The comparison with the parametric figures is worth reading now rather than later. Historical simulation gives the lower VaR of the two and the higher ES. The normal assumption puts the 5% quantile slightly further out than the data do, and then understates how bad things are beyond it. That is the shape of a fat tail, and it is exactly what a single VaR number cannot show.
24.3.3 Monte Carlo VaR and ES
Monte Carlo simulates from the multivariate normal distribution using the estimated covariance matrix, drawing on each language’s own multivariate-normal generator.
MonteCarloVaR = function(w, Sigma, p = 0.05, value = 1000, nsim = 10000) {
Y_sim = mvrnorm(nsim, mu = rep(0, length(w)), Sigma = Sigma)
# Portfolio returns
r_portfolio = Y_sim %*% w
r_sorted = sort(r_portfolio)
idx = ceiling(p * nsim)
VaR = -r_sorted[idx] * value
ES = -EmpiricalES(r_sorted, p) * value
return(list(VaR = VaR, ES = ES))
}
set.seed(42)
res_mc = MonteCarloVaR(w, Sigma_sample, par$probability, par$value)
cat("Monte Carlo VaR:", round(res_mc$VaR, 2), "\n")
cat("Monte Carlo ES: ", round(res_mc$ES, 2), "\n")Monte Carlo VaR: 27.22
Monte Carlo ES: 34.03def MonteCarloVaR(w, Sigma, p=0.05, value=1000, nsim=10000, seed=42):
rng = np.random.default_rng(seed)
Y_sim = rng.multivariate_normal(np.zeros(len(w)), Sigma, size=nsim)
r_portfolio = Y_sim @ w
r_sorted = np.sort(r_portfolio)
idx = int(np.ceil(p * nsim))
VaR = -r_sorted[idx - 1] * value
ES = -EmpiricalES(r_sorted, p) * value
return {'VaR': VaR, 'ES': ES}
res_mc_py = MonteCarloVaR(w_py, Sigma_sample_py, par_py['probability'], par_py['value'])
print(f"Monte Carlo VaR: {res_mc_py['VaR']:.2f}")
print(f"Monte Carlo ES: {res_mc_py['ES']:.2f}")Monte Carlo VaR: 26.82
Monte Carlo ES: 33.37rand(MvNormal(...), nsim) returns assets × simulations, the transpose of R and Python’s orientation, so the draws are transposed before use.
using Random, Distributions, Statistics, Printf
function MonteCarloVaR(w, Sigma, p=0.05, value=1000, nsim=10000; seed=42)
Random.seed!(seed)
draws = rand(MvNormal(zeros(length(w)), Sigma), nsim)'
r_portfolio = draws * w
r_sorted = sort(r_portfolio)
idx = ceil(Int, p * nsim)
VaR = -r_sorted[idx] * value
ES = -EmpiricalES(r_sorted, p) * value
return Dict("VaR" => VaR, "ES" => ES)
end
res_mc_jl = MonteCarloVaR(w_jl, Sigma_sample_jl, par_jl["probability"], par_jl["value"]);
@printf("Monte Carlo VaR: %.2f\n", res_mc_jl["VaR"])
@printf("Monte Carlo ES: %.2f\n", res_mc_jl["ES"])Monte Carlo VaR: 26.92
Monte Carlo ES: 33.84Monte Carlo lands close to the parametric figure, as it must — both draw on the same covariance matrix and the same normality assumption. The gap between them is simulation error and nothing else, which makes it a useful check that the simulation is doing what we think.
24.3.4 EWMA covariance
The sample covariance weights every observation in the window equally. The EWMA recursion from Chapter 22 instead gives more weight to recent observations, producing a next-day forecast \(\CovMatrix_{\SampleSize+1|\SampleSize}\). The recursion is initialised with the full-sample covariance (Sigma_ewma_initial) — a separate quantity from the 500-observation Sigma_sample used above — playing the role of \(\CovMatrix_1\), and run forward through the last observation:
\[\CovMatrix_{t+1} = \EWMAdecay \times \CovMatrix_t + (1-\EWMAdecay) \times \CompoundReturns_t \times \CompoundReturns_t'\]
with \(\EWMAdecay = 0.94\). Running the recursion through \(t = \SampleSize\) (the last observation) gives the one-step-ahead forecast \(\CovMatrix_{\SampleSize+1|\SampleSize}\), which feeds the same parametric VaR formula.
Sigma_ewma_initial = cov(y)
lambda = 0.94
S = Sigma_ewma_initial
for (t in 1:nrow(y)) {
S = lambda * S + (1 - lambda) * y[t, ] %*% t(y[t, ])
}
Sigma_ewma = S
res_ewma = ParametricVaR(w, Sigma_ewma, par$probability, par$value)
cat("EWMA VaR:", round(res_ewma$VaR, 2), "\n")
cat("EWMA ES: ", round(res_ewma$ES, 2), "\n")EWMA VaR: 24.74
EWMA ES: 31.02Sigma_ewma_initial_py = np.cov(y_py, rowvar=False)
lam = 0.94
S_py = Sigma_ewma_initial_py.copy()
for t in range(len(y_py)):
S_py = lam * S_py + (1 - lam) * np.outer(y_py[t, :], y_py[t, :])
Sigma_ewma_py = S_py
res_ewma_py = ParametricVaR(w_py, Sigma_ewma_py, par_py['probability'], par_py['value'])
print(f"EWMA VaR: {res_ewma_py['VaR']:.2f}")
print(f"EWMA ES: {res_ewma_py['ES']:.2f}")EWMA VaR: 24.74
EWMA ES: 31.02using Statistics, Printf
Sigma_ewma_initial_jl = cov(y_jl);
lambda_jl = 0.94;
S_jl = copy(Sigma_ewma_initial_jl);
for t in 1:size(y_jl, 1)
global S_jl = lambda_jl * S_jl + (1 - lambda_jl) * y_jl[t, :] * y_jl[t, :]'
end
Sigma_ewma_jl = S_jl;
res_ewma_jl = ParametricVaR(w_jl, Sigma_ewma_jl, par_jl["probability"], par_jl["value"]);
@printf("EWMA VaR: %.2f\n", res_ewma_jl["VaR"])
@printf("EWMA ES: %.2f\n", res_ewma_jl["ES"])EWMA VaR: 24.74
EWMA ES: 31.02EWMA returns the lowest VaR of the four. The sample covariance averages the whole 500-day window, while EWMA weights the recent past most heavily, so the difference says that the last few weeks were calmer than the window as a whole. Whether that is prudent or complacent depends on whether the calm holds.
24.3.5 Comparison
comparison = data.frame(
Method = c("Parametric", "Historical", "Monte Carlo", "EWMA"),
VaR = round(c(res_param$VaR, res_hs$VaR, res_mc$VaR, res_ewma$VaR), 2),
ES = round(c(res_param$ES, res_hs$ES, res_mc$ES, res_ewma$ES), 2)
)
comparison Method VaR ES
1 Parametric 26.95 33.80
2 Historical 25.89 35.91
3 Monte Carlo 27.22 34.03
4 EWMA 24.74 31.02import pandas as pd
comparison_py = pd.DataFrame({
'Method': ['Parametric', 'Historical', 'Monte Carlo', 'EWMA'],
'VaR': [round(res_param_py['VaR'], 2), round(res_hs_py['VaR'], 2), round(res_mc_py['VaR'], 2), round(res_ewma_py['VaR'], 2)],
'ES': [round(res_param_py['ES'], 2), round(res_hs_py['ES'], 2), round(res_mc_py['ES'], 2), round(res_ewma_py['ES'], 2)]
})
print(comparison_py.to_string(index=False)) Method VaR ES
Parametric 26.95 33.80
Historical 25.89 35.91
Monte Carlo 26.82 33.37
EWMA 24.74 31.02using DataFrames
comparison_jl = DataFrame(
Method = ["Parametric", "Historical", "Monte Carlo", "EWMA"],
VaR = round.([res_param_jl["VaR"], res_hs_jl["VaR"], res_mc_jl["VaR"], res_ewma_jl["VaR"]], digits=2),
ES = round.([res_param_jl["ES"], res_hs_jl["ES"], res_mc_jl["ES"], res_ewma_jl["ES"]], digits=2)
);
println(comparison_jl)4×3 DataFrame
Row │ Method VaR ES
│ String Float64 Float64
─────┼───────────────────────────────
1 │ Parametric 26.95 33.8
2 │ Historical 25.89 35.91
3 │ Monte Carlo 26.92 33.84
4 │ EWMA 24.74 31.02The parametric and Monte Carlo results are similar because both assume multivariate normality. Historical simulation can differ when returns are non-normal. EWMA and the DCC forecast below both respond more quickly to new observations than the sample covariance, and DCC additionally allows correlations themselves to change.
24.4 Using DCC covariances
The DCC model differs from the sample covariance in letting correlations change over time, which the sample covariance cannot do. That is not the same as being better in every case.
Standard DCC buys its tractability with restrictions. A single pair of parameters governs the correlation dynamics of every asset pair, so all pairs adapt at the same speed, and the update responds to the outer product of the standardised shocks, so a shock vector and its sign-reversed counterpart move correlations identically.
Estimation also gets harder as the cross-section grows. Each likelihood evaluation factorises an \(\NumberAssets \times \NumberAssets\) matrix at every date, and the correlation target is singular once the sample is short relative to the number of assets.
Whether the added dynamics repay those restrictions is an empirical question, settled out of sample rather than by assumption. We use the DCC model from Chapter 22, which is R-only here for the reasons given in that chapter.
# Fit DCC model
garch_models = lapply(1:ncol(y_xts), function(i) {
spec = garch_modelspec(y_xts[, i], model = "garch", order = c(1, 1), constant = FALSE)
estimate(spec, keep_tmb = TRUE)
})
names(garch_models) = colnames(y_xts)
multi_garch = to_multi_estimate(garch_models)
dcc_spec = dcc_modelspec(multi_garch, dynamics = "dcc", dcc_order = c(1, 1), distribution = "mvn")
dcc_res = estimate(dcc_spec)Attaching package: ‘data.table’
The following objects are masked from ‘package:reshape2’:
dcast, melt
The following objects are masked from ‘package:lubridate’:
hour, isoweek, mday, minute, month, quarter, second, wday, week,
yday, year
The following objects are masked from ‘package:xts’:
first, last
The following objects are masked from ‘package:zoo’:
yearmon, yearqtrExtract the one-step-ahead conditional covariance forecast:
# Forecast the conditional covariance one step ahead
dcc_pred = predict(dcc_res, h = 1, seed = 42)
Sigma_dcc = tscov(dcc_pred, distribution = FALSE)[, , 1]
Sigma_dcc [,1] [,2]
[1,] 0.0001292289 0.0001102352
[2,] 0.0001102352 0.0004799305res_dcc = ParametricVaR(w, Sigma_dcc, par$probability, par$value)
cat("DCC-based VaR:", round(res_dcc$VaR, 2), "\n")
cat("DCC-based ES: ", round(res_dcc$ES, 2), "\n")DCC-based VaR: 23.69
DCC-based ES: 29.71The DCC forecast gives the lowest VaR of any method here, below even EWMA. Both weight recent data, but DCC also lets the JPM–INTC correlation itself move, and it evidently reads that correlation as currently below its long-run level. A lower correlation means more diversification and a smaller portfolio loss for the same two volatilities.
24.5 Risk decomposition
How much does each asset contribute to total portfolio risk? The Euler decomposition provides an answer that sums to total VaR.
For parametric VaR, the component VaR for asset \(i\) is:
\[\CompVaR_i = \weight_i \frac{(\CovMatrix \weights)_i}{\Vol_\Portfolio^2} \VaR(\probability)\]
where \((\CovMatrix \weights)_i\) is the \(i\)-th element of the vector \(\CovMatrix \weights\). “Component VaR” is distinct from “Conditional VaR” (another name for Expected Shortfall). We decompose on the sample covariance Sigma_sample here, in all three languages. The DCC-based decomposition below stays R-only, consistent with the DCC forecast itself.
ComponentVaR = function(w, Sigma, p = 0.05, value = 1000) {
sigma_p = as.numeric(sqrt(t(w) %*% Sigma %*% w))
VaR_p = -qnorm(p) * sigma_p * value
# Marginal VaR per unit of normalised weight, not per dollar of exposure.
# Component VaR is the same either way: w_i * marginal_per_weight
# equals value * w_i * marginal_per_dollar.
Sigma_w = Sigma %*% w
marginal_per_weight = as.numeric(Sigma_w) / sigma_p * (-qnorm(p)) * value
# Component VaR
component = w * marginal_per_weight
return(list(
total_VaR = VaR_p,
component_VaR = component,
pct_contribution = component / VaR_p * 100
))
}def ComponentVaR(w, Sigma, p=0.05, value=1000):
sigma_p = np.sqrt(w @ Sigma @ w)
VaR_p = -stats.norm.ppf(p) * sigma_p * value
# Marginal VaR per unit of normalised weight, not per dollar of exposure.
Sigma_w = Sigma @ w
marginal_per_weight = Sigma_w / sigma_p * (-stats.norm.ppf(p)) * value
component = w * marginal_per_weight
return {
'total_VaR': VaR_p,
'component_VaR': component,
'pct_contribution': component / VaR_p * 100
}using Distributions
function ComponentVaR(w, Sigma, p=0.05, value=1000)
sigma_p = sqrt(w' * Sigma * w)
VaR_p = -quantile(Normal(), p) * sigma_p * value
Sigma_w = Sigma * w
# Marginal VaR per unit of normalised weight, not per dollar of exposure.
marginal_per_weight = Sigma_w ./ sigma_p .* (-quantile(Normal(), p)) .* value
component = w .* marginal_per_weight
return Dict(
"total_VaR" => VaR_p,
"component_VaR" => component,
"pct_contribution" => component ./ VaR_p .* 100
)
enddecomp = ComponentVaR(w, Sigma_sample, par$probability, par$value)
cat("Total VaR:", round(decomp$total_VaR, 2), "\n\n")
cat("Component VaR:\n")
for (i in 1:length(assets)) {
cat(" ", assets[i], ":", round(decomp$component_VaR[i], 2),
"(", round(decomp$pct_contribution[i], 1), "%)\n")
}
cat("\nSum of components:", round(sum(decomp$component_VaR), 2), "\n")Total VaR: 26.95
Component VaR:
JPM : 10.57 ( 39.2 %)
INTC : 16.38 ( 60.8 %)
Sum of components: 26.95decomp_py = ComponentVaR(w_py, Sigma_sample_py, par_py['probability'], par_py['value'])
print(f"Total VaR: {decomp_py['total_VaR']:.2f}\n")
print("Component VaR:")
for i, a in enumerate(assets_py):
print(f" {a}: {decomp_py['component_VaR'][i]:.2f} ({decomp_py['pct_contribution'][i]:.1f}%)")
print(f"\nSum of components: {sum(decomp_py['component_VaR']):.2f}")Total VaR: 26.95
Component VaR:
JPM: 10.57 (39.2%)
INTC: 16.38 (60.8%)
Sum of components: 26.95using Printf
decomp_jl = ComponentVaR(w_jl, Sigma_sample_jl, par_jl["probability"], par_jl["value"]);
@printf("Total VaR: %.2f\n\n", decomp_jl["total_VaR"])
println("Component VaR:")
for i in 1:length(assets_jl)
@printf(" %s: %.2f (%.1f%%)\n", assets_jl[i], decomp_jl["component_VaR"][i], decomp_jl["pct_contribution"][i])
end
@printf("\nSum of components: %.2f\n", sum(decomp_jl["component_VaR"]))Total VaR: 26.95
Component VaR:
JPM: 10.57 (39.2%)
INTC: 16.38 (60.8%)
Sum of components: 26.95The components sum to total VaR in all three languages. The percentage contributions show each asset’s share of portfolio risk.
24.5.1 DCC-based decomposition (R only)
Since the DCC forecast is R-only, so is its risk decomposition:
decomp_dcc = ComponentVaR(w, Sigma_dcc, par$probability, par$value)
cat("Total VaR (DCC):", round(decomp_dcc$total_VaR, 2), "\n\n")
cat("Component VaR (DCC):\n")
for (i in 1:length(assets)) {
cat(" ", assets[i], ":", round(decomp_dcc$component_VaR[i], 2),
"(", round(decomp_dcc$pct_contribution[i], 1), "%)\n")
}
cat("\nSum of components:", round(sum(decomp_dcc$component_VaR), 2), "\n")Total VaR (DCC): 23.69
Component VaR (DCC):
JPM : 6.84 ( 28.9 %)
INTC : 16.85 ( 71.1 %)
Sum of components: 23.69The two decompositions disagree about how much. On the sample covariance INTC carries about 61% of the risk, and on the DCC forecast about 71%. Equal money was never equal risk, and which covariance we use changes the answer by ten percentage points. A risk budget built on one of these numbers would look wrong under the other.
24.6 Summary
Portfolio VaR extends single-asset methods to account for correlations, and the four methods above trade off differently. The parametric method reads VaR off the portfolio variance \(\weights^\top\CovMatrix \weights\) and is the cheapest to compute, at the price of assuming normality. Historical simulation drops that assumption and takes the empirical quantile of past portfolio returns, so it inherits whatever tail the window holds and adapts slowly. Monte Carlo simulates from the estimated distribution, which accommodates payoffs the other two cannot handle, and costs the most to run. EWMA keeps the parametric formula but replaces the sample covariance with one that weights recent observations more heavily, so it adapts faster to changing conditions.
The DCC model, R-only in this notebook, provides time-varying covariance forecasts that capture changing correlations. Risk decomposition shows where the danger sits in the portfolio. That matters when you set risk budgets or size positions.