Copula-Based Portfolio Risk Analysis
jan kowalski
2026-03-23
1. Data & Monthly Loss Rates
get_monthly_losses <- function(ticker) {
raw <- getSymbols(ticker, src = "yahoo", from = "2021-03-01", to = "2026-03-01",
auto.assign = FALSE)[, 6]
monthly <- as.numeric(Cl(to.monthly(raw)))
round(-diff(monthly) / head(monthly, -1) * 100, 4)
}
nvda_loss <- get_monthly_losses("NVDA")
aramco_loss <- get_monthly_losses("2222.SR")
byd_loss <- get_monthly_losses("BYDDF")
losses <- data.frame(NVDA = nvda_loss, Aramco = aramco_loss, BYD = byd_loss) |>
na.omit()
cat("Observations:", nrow(losses), "\n")## Observations: 59summary(losses)## NVDA Aramco BYD
## Min. :-36.344 Min. :-12.3839 Min. :-36.922
## 1st Qu.:-14.755 1st Qu.: -2.8377 1st Qu.: -6.975
## Median : -5.620 Median : 0.3610 Median : -2.511
## Mean : -5.539 Mean : -0.1617 Mean : -1.554
## 3rd Qu.: 4.827 3rd Qu.: 2.9672 3rd Qu.: 6.653
## Max. : 32.027 Max. : 7.7175 Max. : 19.6102. Marginal Distribution Fitting
fit_best <- function(x, name) {
families <- list(NO = NO, TF = TF, GU = GU, LO = LO)
fits <- lapply(families, \(f) try(gamlss(x ~ 1, family = f, trace = FALSE), silent = TRUE))
aics <- sapply(fits, \(m) if (inherits(m, "gamlss")) AIC(m) else Inf)
best_name <- names(which.min(aics))
cat(name, "-> best family:", best_name, "| AIC:", round(min(aics), 2), "\n")
fits[[best_name]]
}
fit_nvda <- fit_best(losses$NVDA, "NVDA")## NVDA -> best family: NO | AIC: 488.55fit_aramco <- fit_best(losses$Aramco, "Aramco")## Aramco -> best family: GU | AIC: 339.47fit_byd <- fit_best(losses$BYD, "BYD")## BYD -> best family: GU | AIC: 456.533. PIT — Parametric Uniform Transforms
get_pit <- function(model, x) {
fam <- model$family[1]
pfun <- get(paste0("p", fam))
mu <- fitted(model, "mu")
sigma <- fitted(model, "sigma")
if ("nu" %in% model$parameters) {
nu <- fitted(model, "nu")
return(pfun(x, mu = mu, sigma = sigma, nu = nu))
}
pfun(x, mu = mu, sigma = sigma)
}
pit_nvda <- get_pit(fit_nvda, losses$NVDA)
pit_aramco <- get_pit(fit_aramco, losses$Aramco)
pit_byd <- get_pit(fit_byd, losses$BYD)4. Pseudo-Observations (Empirical Uniform Transform)
# pobs is the standard non-parametric route for copula estimation
U <- pobs(as.matrix(losses))
colnames(U) <- c("NVDA", "Aramco", "BYD")
par(mfrow = c(1, 3))
for (i in 1:3) {
hist(U[, i], breaks = 10, freq = FALSE, col = "steelblue", border = "white",
main = colnames(U)[i], xlab = "u", ylim = c(0, 2))
abline(h = 1, col = "red", lty = 2, lwd = 2)
}
par(mfrow = c(1, 1))5. C-Vine Construction (NVDA as root)
Following the PCC decomposition: \[f(x_1,x_2,x_3) = f_1 \cdot f_2 \cdot f_3 \cdot c_{1,2}(F_1,F_2) \cdot c_{1,3}(F_1,F_3) \cdot c_{2,3|1}(F(x_2|x_1), F(x_3|x_1))\]
# Tree 1: unconditional pairs
bc12 <- BiCopSelect(U[,"NVDA"], U[,"Aramco"], familyset = NA, selectioncrit = "AIC",
indeptest = TRUE, level = 0.05)
bc13 <- BiCopSelect(U[,"NVDA"], U[,"BYD"], familyset = NA, selectioncrit = "AIC",
indeptest = TRUE, level = 0.05)
# H-functions: F(Aramco|NVDA), F(BYD|NVDA)
h_aramco_nvda <- pmax(pmin(BiCopHfunc2(U[,"NVDA"], U[,"Aramco"], bc12), 1-1e-10), 1e-10)
h_byd_nvda <- pmax(pmin(BiCopHfunc2(U[,"NVDA"], U[,"BYD"], bc13), 1-1e-10), 1e-10)
# Tree 2: conditional pair
bc23g1 <- BiCopSelect(h_aramco_nvda, h_byd_nvda, familyset = NA, selectioncrit = "AIC",
indeptest = TRUE, level = 0.05)
cat("=== C-VINE STRUCTURE (root: NVDA) ===\n")## === C-VINE STRUCTURE (root: NVDA) ===cat(sprintf("c_12 NVDA–Aramco: %-20s tau=%+.4f\n", bc12$familyname, bc12$tau))## c_12 NVDA–Aramco: Independence tau=+0.0000cat(sprintf("c_13 NVDA–BYD: %-20s tau=%+.4f\n", bc13$familyname, bc13$tau))## c_13 NVDA–BYD: Independence tau=+0.0000cat(sprintf("c_23|1 Aramco–BYD|NVDA: %-20s tau=%+.4f\n", bc23g1$familyname, bc23g1$tau))## c_23|1 Aramco–BYD|NVDA: t tau=+0.99106. Copula Visualization
par(mfrow = c(1, 3))
plot(U[,"NVDA"], U[,"Aramco"], pch = 19, col = rgb(.2,.4,.8,.6), cex = .8,
xlab = "NVDA", ylab = "Aramco",
main = sprintf("c12: %s | τ=%.3f", bc12$familyname, bc12$tau)); grid()
plot(U[,"NVDA"], U[,"BYD"], pch = 19, col = rgb(.8,.3,.2,.6), cex = .8,
xlab = "NVDA", ylab = "BYD",
main = sprintf("c13: %s | τ=%.3f", bc13$familyname, bc13$tau)); grid()
plot(h_aramco_nvda, h_byd_nvda, pch = 19, col = rgb(.2,.7,.3,.6), cex = .8,
xlab = "F(Aramco|NVDA)", ylab = "F(BYD|NVDA)",
main = sprintf("c23|1: %s | τ=%.3f", bc23g1$familyname, bc23g1$tau)); grid()
par(mfrow = c(1, 1))7. Simulation (100 000 paths)
set.seed(42)
N <- 100000
# Simulate from C-vine
sim <- BiCopSim(N, bc12) # (u1, u2)
u1 <- sim[, 1]; u2 <- sim[, 2]
# Conditional: u3 from c13 and c23|1
u3_given_u1 <- BiCopHinv2(u1, runif(N), bc13) # F^{-1}(v|u1)
u3 <- BiCopHinv2(u1,
BiCopHinv2(h_aramco_nvda[sample(nrow(U), N, replace=TRUE)],
runif(N), bc23g1),
bc13)
sim_U <- cbind(u1, u2, u3)
colnames(sim_U) <- c("NVDA", "Aramco", "BYD")
# Back-transform to loss scale using empirical quantiles
q_nvda <- quantile(losses$NVDA, probs = u1, type = 1)
q_aramco <- quantile(losses$Aramco, probs = u2, type = 1)
q_byd <- quantile(losses$BYD, probs = u3, type = 1)
portfolio_loss <- (q_nvda + q_aramco + q_byd) / 3 # equal-weight
VaR_95 <- quantile(portfolio_loss, 0.95)
ES_95 <- mean(portfolio_loss[portfolio_loss > VaR_95])
cat(sprintf("VaR (95%%): %.4f%%\nES (95%%): %.4f%%\n", VaR_95, ES_95))## VaR (95%): 11.4206%
## ES (95%): 14.1451%hist(portfolio_loss, breaks = 100, freq = FALSE, col = "lightblue", border = "white",
main = "Simulated Portfolio Loss Distribution",
xlab = "Monthly Loss (%)")
abline(v = VaR_95, col = "red", lwd = 2, lty = 2)
abline(v = ES_95, col = "darkred", lwd = 2, lty = 1)
legend("topright", legend = c(sprintf("VaR 95%% = %.2f%%", VaR_95),
sprintf("ES 95%% = %.2f%%", ES_95)),
col = c("red", "darkred"), lty = c(2, 1), lwd = 2, bty = "n")