I implement a GARCH-DCC model in Python. My implementation is the following :

def garch_dcc_specification(
    self,
    eps_last: Optional[np.ndarray],
    cond_var_last: Optional[np.ndarray],
    q_last_t: Optional[np.ndarray],
) -> SpecResult:
    if eps_last is None:
        eps_last = np.zeros(self.n)

    if q_last_t is None:
        q_last_t = np.zeros((self.n, self.n))

    epsilon_square_last = np.array([eps_last_i ** 2 for eps_last_i in eps_last])
    # first, evaluate the garch cond. variance. 
    # (garch_alpha, beta and omega are 1D arrays)
    cond_var_t = (self.garch_omega
                  + self.garch_alpha * epsilon_square_last
                  + self.garch_beta * (cond_var_last if cond_var_last is not None else np.zeros(self.n)))

    d_t = np.diag([math.sqrt(v) for v in cond_var_t])

    if cond_var_last is not None:
        d_last_t = np.diag([math.sqrt(v) for v in cond_var_last])
        v_last_t = inv(d_last_t).dot(eps_last)
    else:
        v_last_t = np.zeros(self.n)

    # DCC specification for the conditional correlation
    # note: DO NOT DO v[t - 1].dot(v[t - 1].transpose()) : since v[t - 1] is a 1D array, result would be a number
    q_t = (self.dcc_r * (1 - self.dcc_alpha - self.dcc_beta)
                       + self.dcc_alpha * v_last_t.reshape(1, -1).transpose().dot(v_last_t.reshape(1, -1))
                       + self.dcc_beta * q_last_t)

    # standardize q to get a real correlation matrix 
    r_t = np.zeros((self.n, self.n))
    for i in range(self.n):
        for j in range(self.n):
            r_t[i][j] = q_t[i][j] / math.sqrt(q_t[i][i] * q_t[j][j])

    # transforms to a variance-covariance matrix by incorporing the cond variances 
    h_t = d_t.dot(r_t).dot(d_t)
    return GarchDccParams.SpecResult(cond_var = cond_var_t, q = q_t, h = h_t)

def generate_innovations(self, length: int) -> np.ndarray:
    innovations = np.zeros((length + 1, self.n))
    spec_res: List[GarchDccParams.SpecResult] = []
    for t in range(0, length + 1):
        spec_res.append(self.garch_dcc_specification(
            eps_last = innovations[t - 1] if t != 0 else None,
            cond_var_last = spec_res[t - 1].cond_var if t != 0 else None,
            q_last_t = spec_res[t - 1].q if t != 0 else None,
        ))
        innovations[t] = np.random.multivariate_normal(np.zeros(self.n), spec_res[t].h)

To check my implentation, I control that the generate_innovation() empirical pearson correlation coefficient (with np.corrcoef) is equal to the input self.dcc_r matrix correlation coefficient, which should be the unconditional correlation of the overrall generated innovation, if I understand correctly.

When running with constant variance in the GARCH (garch alpha and beta = 0), I get a good pearson coef coefficient that is equal to the one I set in the input self.dcc_r Though, when the conditional variance is moving (garch alpha and beta > 0), I don't get the same coefficient, I get always a lesser empirical correlation coefficient then expected in the DCC_R. For example, when running with 10000 points, and with an input dcc_r correlation coef. of 0.9, I get an empirical unconditional correlation, in my generated innovations, of around 0.75

PS: To simplificate, I set dcc_alpha and dcc_beta to 0 so only the dcc_r matrix is taken into account (we have a GARCH-CCC model instead of a GARCH-DCC, cond. correlation is always the same). The "problem" (if it is one) still occurs, still when GARCH alpha/beta > 0.

Is it normal ?