I want to simulate some exchange rates with a DCC GARCH. I know the package rmgarch but I want to code the simulation my self. The following are the main equations of the model:

$r_t = a_t$

$a_t = H^{1/2}_{t} Z_t$

$H_t = D_t R_t D_t$

$R_t = Q^{*-1}_t Q_t Q^{*-1}_t$

$Q_t = (1-a-b)\bar{Q} + a\epsilon_{t-1}\epsilon'_{t-1} + b Q_{t-1}$

My question is: if I all ready have the univariate GARCH parameters for $D_t$, $a$ and $b$ and $\bar{Q}$ how do I get $\epsilon_{t-1}$ to perform the simulations? and what are the steps to perform it.

  • 1
    $\begingroup$ you set the first errors term to zero (their mean) and then derive the other ones by fitting data to your model (recursively for each $t$ - by going backward from the origin) $\endgroup$
    – Malick
    Jul 14, 2018 at 14:39
  • $\begingroup$ @Malick I don't think is like that because with what data will i compare it if it's a simulation? $\endgroup$ Jul 14, 2018 at 16:46
  • $\begingroup$ @Malick but isn't $Z_t$ the variable that follows a N(0,1)? $\endgroup$ Jul 16, 2018 at 16:25
  • $\begingroup$ Sorry Alejandro, I have deleted my answer as it was wrong, I'll try to come up with a proper answer. $\endgroup$
    – Malick
    Jul 18, 2018 at 10:58
  • $\begingroup$ @Malick thanks anyway if you manage to come up with the right answer i would appreciated $\endgroup$ Jul 18, 2018 at 15:17

1 Answer 1


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:

  1. 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).

  2. 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.

  3. 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)

# 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


# 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
  # returns
  lOut[["DCC_returns"]] = DCC_returns

  • $\begingroup$ How is it done for the CCC-GARCH model? @Pleb $\endgroup$
    – Adrian
    Jun 21, 2023 at 16:56
  • $\begingroup$ @Adrian As specified under my "some comments" section in my answer, you recover the CCC model by setting acor equal to the unconditional correlation. This can be done by setting $a=0$ and $b=0$ in the DCC model, so the unconditional correlation is recovered from the factorization $R = \widetilde{Q}^{-1/2} Q \widetilde{Q}^{-1/2}$. $\endgroup$
    – Pleb
    Jun 21, 2023 at 18:15
  • $\begingroup$ Are the returns simple returns or log returns? @Pleb $\endgroup$
    – Adrian
    Aug 9, 2023 at 11:42
  • $\begingroup$ @Adrian purely from a mathematical perspective, there is no specification on whether it's log- or simple-returns (I have never seen one). However, we mainly (almost always) use log-returns when dealing with GARCH models. Therefore, I would treat it as such. $\endgroup$
    – Pleb
    Aug 11, 2023 at 6:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.