How to implement Girardi & Ergun's (2013) three-step multivariate GARCH estimation of CoVaR in R?

Question

I'm trying to calculate multivariate GARCH estimation of conditional value-at-risk, by adopting a three-step model from Girardi & Ergun (2013) paper entitled "Systemic risk measurement: Multivariate GARCH estimation of CoVaR". I'm using R and up until now I can't find any good references as to how to construct the model. For your reference, this is the three-step calculation from Girardi & Ergun (2013):

Step 1

enter image description here

Step 2

enter image description here

enter image description here

Step 3

enter image description here

I have tried to put some prompt into chatGPT, and this is what it gave me.

# Install and load required packages
install.packages("quantmod")
library(quantmod)
library(rmgarch)

# Step 1: Compute VaR of each stock using univariate model

# Define the stock symbols and date range
symbols <- c("YAHOO", "GOOG")
start_date <- "2018-01-01"
end_date <- "2022-12-31"

# Download stock price data
getSymbols(symbols, from = start_date, to = end_date)

# Extract adjusted close prices
prices <- merge(Ad(YAHOO), Ad(GOOG))
colnames(prices) <- symbols

# Compute daily returns
returns <- na.omit(ROC(prices))

# Function to compute VaR
computeVaR <- function(returns, alpha = 0.05) {
  VaR <- matrix(NA, ncol = ncol(returns))
  for (i in 1:ncol(returns)) {
    VaR[, i] <- quantile(returns[, i], alpha)
  }
  return(VaR)
}

# Compute VaR for each stock
VaR <- computeVaR(returns)

# Step 2: Estimate bivariate GARCH model using Engle's (2002) DCC specification

# Combine the returns into a multivariate matrix
data <- as.matrix(returns)

# Function to estimate DCC-GARCH model
estimateDCC <- function(data) {
  spec <- dccspec(uspec = multispec(replicate(2, ugarchspec(variance.model = list(model = "sGARCH")))),
                  dccOrder = c(1, 1), distribution = "mvnorm")
  fit <- dccfit(spec, data = data, out.sample = 0)
  return(fit)
}

# Estimate DCC-GARCH model
dcc_fit <- estimateDCC(data)

# Step 3: Calculate CoVar measurement for each stock using double integral

# Function to calculate CoVaR
calculateCoVaR <- function(returns, VaR, dcc_fit, alpha = 0.05) {
  CoVaR <- matrix(NA, ncol = ncol(returns))
  for (i in 1:ncol(returns)) {
    conditional_variance <- as.matrix(rmgarch::rcov(dcc_fit)[,,i])
    exceedances <- t(matrix(VaR[, i] > returns[, i], nrow = nrow(returns), ncol = ncol(returns)))
    CoVaR[, i] <- alpha + (1/alpha) * sum(conditional_variance * exceedances)
  }
  return(CoVaR)
}

# Calculate CoVaR for each stock
CoVaR <- calculateCoVaR(returns, VaR, dcc_fit)

# Print the results
print("VaR:")
print(VaR)
print("CoVaR:")
print(CoVaR)

I have tried run those code, somehow it work, but I'm not 100% convinced. I have tried to calculate step 1 in Excel and step 2 using Eviews, the result from step 2 is quite different. Calculation in step 2 was supposed to bring out the conditional variance series from the model GARCH(1,1), then calculate the VaR from there.

Now what I have in mind is:

1. Is there any way to prove that those output chatGPT gave me, is the similar calculation following Girardi & Ergun (2013)?
2. And if Q no. 1 is correct, I'm sure that calculation in step 3 is wrong, because from Girardi & Ergun (2013), I have to calculate double integral from joint PDF functions of 2 different stocks. How can I construct the function to calculate step 3?

P.S: I'm apologize in advanced if my question seems messed up, I'm a beginner in here and this is only my 3rd time asking around.