How to calculate CoVaR and delta CoVaR using multivariate DCC-GJR-GARCH

Question

I'm trying to calculate conditional value-at-risk (CoVaR/CVaR) and delta conditional value-at-risk using multivariate DCC-GJR-GARCH, but I'm stuck. I'm using three step procedure, the first is calculate value-at-risk using univariate GJR-GARCH model, second step is build a joint distribution using DCC-GARCH, then I'm lost at step three where many of my references paper instruct to solve double integration numerically to calculate CoVaR.

This is the R code that I'm using.

library(tidyverse)
library(quantmod)
library(rugarch)
library(rmgarch)

# Download dataset
d1=getSymbols("BTC-USD",auto.assign=F)|>Ad()|>ROC()|>fortify.zoo()
d2=getSymbols("^GSPC",auto.assign=F)|>Ad()|>ROC()|>fortify.zoo()
colnames(d1)=c("Index","BTC")
colnames(d2)=c("Index","SP500")
d3=na.omit(merge(d1,d2,all=F))

# Step 1 (univariate GJR-GARCH model)
m1=ugarchspec(variance.model=list(model="gjrGARCH",garchOrder=c(1,1)),
              mean.model=list(armaOrder=c(0,0),include.mean=T),
              distribution.model="std")
m2=ugarchfit(spec=m1,data=d3$BTC)
m3=ugarchfit(spec=m1,data=d3$SP500)

d4=fitdist(distribution='std',x=d3$BTC)$pars[3]
d5=qdist(distribution='std',shape=d4,p=0.05)
d6=fitdist(distribution='std',x=d3$SP500)$pars[3]
d7=qdist(distribution='std',shape=d6,p=0.05)

d8=m2@fit$sigma*d5
d9=m3@fit$sigma*d7
d10=cbind(d3,d8,d9)
colnames(d10)[4]="BTCvol"
colnames(d10)[5]="SP500vol"

f1=ggplot(data=d10,aes(x=Index))+geom_line(aes(y=BTCvol,color="BTC"))+
  geom_line(aes(y=SP500vol,color="SP500"))+ylab("")+xlab("")+
  scale_color_manual(name="Volatility Series", values=c("BTC"="red","SP500"="darkblue"))+
  theme_light()+theme(legend.position="top")

# Step 2 (DCC-GARCH specification)
d11=d10[,c("BTC","SP500")]
m4=dccspec(uspec=multispec(replicate(2,m1)),dccOrder=c(1,1),distribution="mvt")
m5=dccfit(m4,data=d11)
d12=m5@model$sigma

That was the last line of code that I wrote, I'm stuck as I don't know what else to do. While my goal is to make something like this:

And something like this:

Now my question is that, is there any way to turn output from d12=m5@model$sigma into a time-varying value-at-risks? I have tried to contact the corresponding author but I got no luck :(