As an Amazon Associate, I earn from qualifying purchases. For more practice problems like this, see Schaum’s Outline of Statistics, 6th Edition.

1. What the student is asking (in plain language)

The student wants to compute CoVaR – the value‑at‑risk of institution i conditional on institution j being in distress – using a bivariate DCC‑GARCH model fitted in R.

The three‑step recipe is

  1. Univariate GARCH → obtain the one‑step‑ahead VaR for each series.
  2. Bivariate DCC‑GARCH → obtain the conditional joint density (f_t(x_i,x_j)) of the two returns for the same horizon.
  3. Solve

[ \Pr!\bigl(X_t^{i}\le \text{CoVaR}^{i|j}{\beta,t}, \;X_t^{j}\le \text{VaR}^{j}{\alpha,t}\bigr)=\alpha\beta ]

for the unknown (\text{CoVaR}^{i j}_{\beta,t}).

The difficulty is:

  • How do we turn the DCC‑GARCH output into a bivariate density?
  • Once we have the density, how do we evaluate the double integral and find the CoVaR value?

Below is a complete worked solution that shows every step, provides the required R code, and explains the underlying mathematics.

2. Step‑by‑step solution

2.1 Model set‑up and notation

  • Let (\mathbf{X}t=(X{i,t},X_{j,t})^{\prime}) be the (zero‑mean) return vector at time (t).
  • The DCC‑GARCH model writes the conditional covariance matrix as

[ \mathbf{H}_t = \mathbf{D}_t \,\mathbf{R}_t \,\mathbf{D}_t, ]

where

  • (\mathbf{D}t = \operatorname{diag}(\sigma{i,t},\sigma_{j,t})) with (\sigma_{i,t}^2,\sigma_{j,t}^2) the univariate GARCH(1,1) variances, and
  • (\mathbf{R}_t) is the conditional correlation matrix generated by the DCC dynamics.

Conditional distribution.
The most common assumption (and the one used by rmgarch) is that, conditional on information (\mathcal{F}_{t-1}),

[ \mathbf{X}t \mid \mathcal{F}{t-1} \;\sim\; \mathcal{N}\bigl(\mathbf{0},\mathbf{H}_t\bigr) ]

or a multivariate‑Student‑(t) with (\nu) degrees of freedom.
Hence the joint density is simply the multivariate normal (or t) density evaluated at (\mathbf{x}=(x_i,x_j)) with covariance (\mathbf{H}_t).

2.2 From DCC parameters to the conditional density

R packages

Task Package / function
Univariate GARCH (step 1) rugarch::ugarchfit
Bivariate DCC‑GARCH (step 2) rmgarch::dccfit
Evaluate a bivariate normal density mvtnorm::dmvnorm
Evaluate a bivariate t density mvtnorm::dmvt
Evaluate a bivariate CDF (the double integral) mvtnorm::pmvnorm (normal) or mvtnorm::pmvt (t)

2.2.1 Fit the two univariate GARCH models

library(rugarch)

# Example: daily returns of two banks
ret_i <- bank_i$ret   # vector of returns for institution i
ret_j <- bank_j$ret   # vector of returns for institution j

# GARCH(1,1) specification (Gaussian innovations)
spec_i <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                    variance.model = list(garchOrder = c(1,1)),
                    distribution.model = "norm")
spec_j <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
                    variance.model = list(garchOrder = c(1,1)),
                    distribution.model = "norm")

fit_i <- ugarchfit(spec = spec_i, data = ret_i)
fit_j <- ugarchfit(spec = spec_j, data = ret_j)

From these fits we later obtain the one‑step‑ahead conditional volatilities
(\hat\sigma_{i,t+1},\hat\sigma_{j,t+1}) that are needed for (\mathbf{D}_{t+1}).

2.2.2 Fit the DCC‑GARCH model

library(rmgarch)

# Combine the two series in a matrix (columns = assets)
data_mat <- cbind(ret_i, ret_j)

# Univariate specs for the DCC (must be the same as above)
uspec <- multispec(replicate(2, spec_i))   # both series use the same spec object

# DCC specification (Gaussian)
dccspec <- dccspec(uspec = uspec,
                   dccOrder = c(1,1),
                   distribution = "mvnorm")   # or "mvt" for Student‑t

dccfit <- dccfit(dccspec, data = data_mat)

The object dccfit contains, for any forecast horizon (h) (usually (h=1)):

  • Conditional variances sigma(dccfit) → a matrix with (\hat\sigma_{i,t+1}^2, \hat\sigma_{j,t+1}^2).
  • Conditional correlation matrix rcor(dccfit) → a 2×2 matrix (\hat{\mathbf R}_{t+1}).

From these we build the conditional covariance matrix:

# One‑step‑ahead forecasts
sigma_f   <- sigma(dccfit)               # 2 × (n) matrix, last column = t+1
R_f       <- rcor(dccfit)                # 2 × 2 × (n) array, last slice = t+1

# Extract the t+1 values
sigma_t1  <- sigma_f[, ncol(sigma_f)]               # vector (σ_i, σ_j)
R_t1      <- R_f[,, ncol(sigma_f)]                   # 2×2 matrix

# Build H_{t+1}
D_t1      <- diag(sigma_t1)                         # diag(σ_i, σ_j)
H_t1      <- D_t1 %*% R_t1 %*% D_t1                 # conditional covariance

Now we have (\mathbf{H}_{t+1}), the scale matrix of the conditional bivariate normal (or t) distribution.

2.2.3 The joint density

If the DCC model was estimated with distribution = "mvnorm":

[ f_{t+1}(x_i,x_j) \;=\; \text{dmvnorm}\bigl(\mathbf{x};\; \mathbf{0},\; \mathbf{H}_{t+1}\bigr). ]

If distribution = "mvt" with (\nu) degrees of freedom (available via model$shape in the fitted object):

[ f_{t+1}(x_i,x_j) \;=\; \text{dmvt}\bigl(\mathbf{x};\; \mathbf{0},\; \mathbf{H}_{t+1},\;\nu\bigr). ]

Thus the density is known analytically; we only need to evaluate it at points ((x_i,x_j)).

2.3 Computing the joint probability required for CoVaR

The definition we need to satisfy is

[ \Pr!\bigl( X_i \le c,\; X_j \le v_j \bigr) = \alpha\beta, ]

where

  • (v_j = \text{VaR}^{j}_{\alpha,t+1}) (already obtained from step 1),
  • (c = \text{CoVaR}^{i j}_{\beta,t+1}) is the unknown we are solving for.

The left‑hand side is the bivariate CDF evaluated at ((c, v_j)):

[ F_{t+1}(c, v_j) \;=\; \int_{-\infty}^{c}\int_{-\infty}^{v_j} f_{t+1}(x_i,x_j)\,dx_j\,dx_i . ]

The mvtnorm package provides exact CDF calculators, so we never have to code a double integral ourselves.

library(mvtnorm)

# VaR of institution j (step 1)
alpha <- 0.05               # example confidence level
VaR_j  <- quantile(ret_j, probs = alpha, type = 7)   # empirical quantile
# or, if you want the model‑based VaR:
# VaR_j <- qnorm(alpha, mean = 0, sd = sigma_t1[2])   # Gaussian
# VaR_j <- qt(alpha, df = nu, mean = 0, sd = sigma_t1[2])  # t‑dist

# Function that returns the joint CDF for a candidate c (CoVaR)
joint_cdf <- function(c) {
  lower <- c(-Inf, -Inf)                 # integration limits
  upper <- c(c, VaR_j)                   # (c , VaR_j)
  if (dist == "norm") {
    pmvnorm(lower = lower, upper = upper,
            mean = c(0,0), sigma = H_t1)  # Gaussian CDF
  } else {
    pmvt(lower = lower, upper = upper,
          df = nu, delta = c(0,0), sigma = H_t1)  # t‑CDF
  }
}

joint_cdf(c) returns the probability (\Pr(X_i\le c, X_j\le v_j)).

2.4 Solving for CoVaR

We need the root of

[ g(c) \;=\; \text{joint_cdf}(c) \;-\; \alpha\beta \;=\;0 . ]

Because the CDF is monotone increasing in (c), a simple uniroot (bisection) works.

beta  <- 0.05                # second confidence level
target <- alpha * beta

# Choose a reasonable search interval.
# Lower bound: a very low quantile (e.g., -5σ). Upper bound: 0 (or a positive number).
lower_bound <- -10 * sigma_t1[1]   # -10 standard deviations
upper_bound <-   0                 # CoVaR will usually be below zero for losses

# Root-finding
CoVaR_i_j <- uniroot(function(c) joint_cdf(c) - target,
                     lower = lower_bound,
                     upper = upper_bound,
                     tol = 1e-8)$root
CoVaR_i_j

CoVaR_i_j is the one‑step‑ahead CoVaR of institution i conditional on institution j being at its (\alpha)-VaR level.

2.5 Closed‑form shortcut (Gaussian case)

If the conditional joint distribution is exactly Gaussian, the double integral can be expressed analytically, yielding a simple formula:

[ \text{CoVaR}^{i|j}{\beta,t} \;=\; \mu{i,t}

  • \sigma_{i,t}\, \Phi^{-1}!\bigl(\beta \mid \Phi^{-1}(\alpha) \rho_{t}\bigr), ]

where

  • (\Phi) = standard normal CDF,
  • (\Phi^{-1}(\cdot \mid \cdot)) = conditional normal quantile,
  • (\rho_t) = conditional correlation (\hat R_{t,ij}).

In R this becomes:

```r rho_t <- R_t1[1,2] z_alpha <- qnorm(alpha) z_beta <- qnorm(beta)

Conditional mean of X_i given X_j = VaR_j

cond_mean <- rho_t * sigma_t1[1] / sigma_t1[2] * VaR_j cond_sd <- sigma

Original question: Implementation of CoVaR (a systemic risk measure) in R on Cross Validated (Stats Stack Exchange), licensed CC BY-SA.