I implement a GARCH-DCC model in Python, for number of asset = 2. 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 ?