Simulating a DCC-GARCH(1,1) model
Given that you already have a given set of proper defined parameters for the DCC-GARCH model, the standardized residuals $\varepsilon_{t-1}$ are recovered from the univariate GARCH models and fed into the DCC structure to yield your simulated correlation dynamics. The simulation of a DCC-GARCH model can be defined in 3 overall steps:
Simulate univariate GARCH dynamics using a filter algorithm for pre-defined parameters and recover the standardized residuals. After that, construct a general function that simulates $d$ GARCH models with different sets of parameters (or the same set of parameters).
Use the standardized residuals from (1) to simulate the DCC dynamics and recover the correlation structure. From construction, we know that $\varepsilon_{t-1} \overset{iid}{\sim} N(0,1)$ and thus you can simulate the DCC dynamics and back out the correlation without the GARCH component.
Recover the covariance matrix via the univariate variances and the correlation structure using the following formula, $\Sigma_t=D_t^{1/2}R_t D_t^{1/2}$, with $D_t = diag(\sigma_{i,t}^2)$
Throughout the answer, the code will be written in R. Moreover, the code will be non-optimized.
1. Simulating multiple univariate GARCH models
Formally let us define the (demeaned) returns in a univariate setting,
\begin{equation}
r_t = \sigma_t \epsilon_t \iff \frac{r_t}{\sigma_t} = \epsilon_t,
\end{equation}
where $\epsilon_t \overset{iid}{\sim}N(0,1)$ are univariate standardized residuals. This further implies that $r_t \sim N(0, \sigma_t^2)$, which we will be using in the simulation. The variance dynamics of the univariate returns follows a GARCH(1,1) on the form:
\begin{equation}
\sigma^2_t = \omega + \alpha r_{t-1}^2 + \beta \sigma^2_{t-1}.
\end{equation}
The way to simulate a GARCH model we initialize the simulation by sampling $r_1 \sim N(0,\sigma^2_1)$, where $\sigma^2_1 = \frac{\omega}{1-\alpha - \beta}$ is the unconditional variance. After that, we recusively update the variance dynamics and simulate the returns using the rnorm() function.
The function is presented below under the name GARCHSim and has been extensively commented. For ease of implementation, we Construct a function (GARCHmultisim) that performs multiple simulations in one function call. Here, the parameters variable (params) will be a $d \times 3$ matrix with $d$ denoting the amount of GARCH models to simulate. The function outputs a list of two matrices with dimensions $sims \times d$, where one contains the standardized residuals and the other contains the variance.
2 + 3. Simulating the DCC structure using standardized residuals
Constructing the filter algorithm for the DCC structure is, more or less, the same procedure as with the GARCH model. We need to update $Q_t$ for each time-point and then reconstruct the correlation matrix following $R_t = \widetilde{Q}_t^{-1/2} Q_t \widetilde{Q}_t^{-1/2}$ with $\widetilde{Q}_t$ containing the diagonal elements of $Q_t$. Then we can reconstruct the covariance matrix, $\Sigma_t$, using the relation in point 3.
We further model $Q_t$ following a GARCH(1,1)-like structure where we impose weak-stationarity (ie. $a+b<1$) to ensure that $Q_t$ is positive definite and implying that $R_t$ is positive definite:
\begin{equation*}
Q_t = \bar{Q} (1-a-b) + a (\eta_{t-1}\eta_{t-1}^\intercal) + b Q_{t-1},
\end{equation*}
where $\bar{Q}$ is the unconditional covariance matrix of the multivariate standardized residuals $\eta_t$. We can recover both parameters following
\begin{align*}
\bar{Q} &= \frac{1}{T} \sum_{t=1}^{T} \eta_t \eta_t^\intercal\\
\eta_t &= D_t^{-\frac{1}{2}} r_t \sim N(0, R_t),
\end{align*}
which we will utilize in our simulation procedure. The function has been named DCCsim and can be found in the end of the code snippet.
Graphical illustration
I have provided a graphical illustration of a simulated bivariate DCC-GARCH(1,1) model with arbitrary parameters given by:
$$\omega_{1}=0.02, \quad \alpha_1 = 0.08, \quad \beta_1 = 0.89, \quad \omega_{2}=0.02, \quad \alpha_2 = 0.05, \quad \beta_2 = 0.94, \quad a=0.1, \quad b=0.89,$$
with $\omega_i, \: \alpha_i, \: \beta_i$ for $i=1,2$ being the parameters for the first and second GARCH process, which are also depicted below as the red and green one respectively.
Some comments
The correlation structure in the DCC model is dependent on the residuals of the GARCH models. Essentially, the univariate GARCH models do not capture the dependence structure and thus it will be contained within the residuals. However, the residuals follows a univariate standard Gaussian distribution (and are iid) thus having no correlation structure whatsoever. This implies that the unconditional correlation tends toward zero as $sims \rightarrow \infty$.
You can also simulate multivariate normal variates with dynamic correlations in the DCC model, where the correlations are estimated from real world stock prices (that is, you can use the mvrnorm function in R with a changing covariance matrix, based on different correlations), or you can try to use alternative distributions for the residuals.
The Constant Conditional Correlation (CCC) model is recovered by reconstructing the covariance matrix (DCC_COVAR) where you replace acor with the unconditional correlation.
Any extensions in the univariate GARCH models is valid due to the separation of the correlation- and volatility dynamics.
The above implementation provides a basic simulation framework. Consequently, there might be more efficient ways to construct realistic simulated correlations within the DCC structure or completely alternative ways to simulate this model. I hope this atleast provides some insight for you and future readers.
Code appendix
The code is presented below.
# Univariate GARCH simulation
GARCHSim <- function(omega, alpha, beta, sims) {
#sims (integer): Number of simulations performed.
## initialize the vector of simulated returns and variances
ret = numeric(sims)
sigma2 = numeric(sims)
## initialize the variance at time t = 1 with the unconditional variance.
sigma2[1] = omega/(1.0 - alpha - beta)
## sample the first observations
ret[1] = rnorm(1, mean = 0, sd = sqrt(sigma2[1]))
## loop over sims.
for (t in 2:sims) {
#update volatility
sigma2[t] = omega + alpha * ret[t - 1]^2 + beta * sigma2[t - 1]
# sample new observarion
ret[t] = rnorm(1, mean = 0, sd = sqrt(sigma2[t]))
}
## we return a list with three components: the sampled returns, the variance,
## and the standardized residuals.
lOut = list()
lOut[["ret"]] = ret
lOut[["sigma2"]] = sigma2
# gives standardized residuals for demeaned returns.
lOut[["stdres"]] = ret/sqrt(sigma2)
return(lOut)
}
# Multiple simulations of the GARCH model.
GARCHmultisim <- function(params, sims){
multisim_res <- matrix(0L, nrow=sims, ncol = nrow(params))
multisim_sigma <- matrix(0L, nrow=sims, ncol = nrow(params))
# rows in parameter matrix determines amount of garch simulations.
for(i in 1:nrow(params)){
# matrix() transforms it into a vector.
multisim_res[, i] <- matrix(GARCHSim(params[i,1], params[i,2], params[i,3], sims)[["stdres"]])
multisim_sigma[, i] <- matrix(GARCHSim(params[i,1], params[i,2], params[i,3], sims)[["sigma2"]])
}
lOut = list()
lOut[["res"]] = multisim_res
lOut[["sigma"]] = multisim_sigma
return(lOut)
}
# Simulation for the Dynamic Conditional Correlations (DCC).
DCCsim <- function(GARCH_params, a, b, sims){
stdres <- GARCHmultisim(GARCH_params, sims)[["res"]]
Gsigs <- GARCHmultisim(GARCH_params, sims)[["sigma"]]
mQ = cor(stdres)
iN = ncol(stdres)
iT = nrow(stdres)
# initialize the array for the correlations
aCor = array(0, dim = c(iN, iN, iT))
# initialize the array for the Q matrices
aQ = array(0, dim = c(iN, iN, iT))
## initialization at the unconditional cor
aCor[,, 1] = mQ
aQ[,,1] = mQ
# main loop
for (t in 2:iT) {
# update the Q matrix
aQ[,, t] = mQ * (1 - a - b) + a * t(stdres[t - 1, , drop = FALSE]) %*%
stdres[t - 1, , drop = FALSE] + b * aQ[,,t - 1]
## Compute the correlation as Q_tilde^{-1/2} Q Q_tilde^{-1/2}
aCor[,, t] = diag(sqrt(1/diag(aQ[,, t]))) %*% aQ[,, t] %*%
diag(sqrt(1/diag(aQ[,, t])))
}
#Getting univariate variance matrix D_t:
msigma <- array(0L, dim = c(ncol(Gsigs), ncol(Gsigs), nrow(Gsigs)))
for(i in 1:ncol(Gsigs)){
msigma[i, i, ] <- Gsigs[, i]
}
# constructing the covariance matrix
DCC_COVAR <- array(0, dim = c(ncol(Gsigs), ncol(Gsigs), nrow(Gsigs)))
DCC_returns <- matrix(0, nrow = nrow(Gsigs), ncol = ncol(Gsigs))
for(i in 1:sims){
DCC_COVAR[,,i] <- msigma[,,i]^0.5 %*% aCor[,,i] %*% msigma[,,i]^0.5
# Getting returns:
DCC_returns[i,] <- mvrnorm(1, rep(0,ncol(Gsigs)), DCC_COVAR[,,i])
}
lOut = list()
# correlations
lOut[["aCor"]] = aCor
# univariate variances
lOut[["Gsigs"]] = Gsigs
# covariances
lOut[["DCC_COVAR"]] = DCC_COVAR
# returns
lOut[["DCC_returns"]] = DCC_returns
return(lOut)
}
```
