DCC-garch simulation not stationary?

Question

I am trying to replicate some outcomes from Becker et al.(2015). To recover table 2, I need to simulate returns from the DCC-Garch model. The parameters and unconditional values(matrices) have been estimated from the sample:

DCC: alpha, beta = 0.005559, 0.990274
     Q_bar = 
     matrix([[1.04272768, 0.24154633, 0.3195213 , 0.22292329, 0.2680927 ],
             [0.24154633, 0.9626568 , 0.26479132, 0.24719453, 0.32599102],
             [0.3195213 , 0.26479132, 1.1465517 , 0.33818774, 0.33139967],
             [0.22292329, 0.24719453, 0.33818774, 1.02975989, 0.27146692],
             [0.2680927 , 0.32599102, 0.33139967, 0.27146692, 1.03497847]])

Garch: omega = [1.71*1e-5, 2.71*1e-6, 9.14*1e-7, 5.84*1e-6, 1.37*1e-6]
       alpha = [0.1058, 0.0346, 0.0312, 0.0589, 0.0172]
       beta = [0.8344, 0.9600, 0.9663, 0.9311, 0.9798]

#The parameters are estimated from 5 time series, i.e., the return rate of 5 stocks.

I wrote the following codes to simulate return rates; however, the data looks weird because they don't seem to be stationary, e.g. some return rates are beyond 1e+13. I made sure that alpha + beta <1 so that the stationary condition was satisfied, and I am totally confused.

#----generate simulated data from DCC
#700 returns and 1000 repeated times
#Since it is too slow, just simu for 100 times
N = 5
K_simu = 100
T_simu = 700
mu_i = r_it.mean()
R_simu_t = np.zeros((N, N, T_simu, K_simu))
z_simu_t = np.zeros((N, T_simu, K_simu))
r_simu_t = np.zeros((N, T_simu, K_simu))
Q_simu_t = np.zeros((N, N, T_simu, K_simu))
#here I don't define sigma^2
sigma_simu_sq_t = np.zeros((N, T_simu, K_simu))
Sigma_simu_t = np.zeros((N, N, T_simu, K_simu))
#alpha, beta = 0.005559, 0.990274
DCCalpha, DCCbeta = dccResDF['alpha'][0], dccResDF['beta'][0]
# Garchomega, Garchalpha, Garchbeta = garchResDF['omega'], garchResDF['alpha'], garchResDF['beta']
Garchomega = np.asarray([1.71*1e-5, 2.71*1e-6, 9.14*1e-7, 5.84*1e-6, 1.37*1e-6])
Garchalpha = np.asarray([0.1058, 0.0346, 0.0312, 0.0589, 0.0172])
Garchbeta = np.asarray([0.8344, 0.9600, 0.9663, 0.9311, 0.9798])
#Get unconditional values
Q_bar = np.asmatrix(np.zeros((N,N)))
for t in range(T):
    Q_bar = Q_bar + np.matmul(z_it[t,:].T, z_it[t,:])    
Q_bar = Q_bar/(T-1)
DQsqrt = np.diag(np.sqrt(np.diag(Q_bar)) ** -1)
R_0 = np.matmul(np.matmul(DQsqrt, Q_bar), DQsqrt)
D_0 = np.diag(np.sqrt(sigma_sq_it[0, :]))
Sigma_t_0 = np.matmul(np.matmul(D_0, R_0), D_0)
#Recuisively simulate z and r
for k in range(K_simu):
    Q_simu_t[:, :, 0, k] = Q_bar
    R_simu_t[:, :, 0, k] = R_0
    rand = np.random.randn(5)
    z_simu_t[:, 0, k] = np.matmul(rand, np.sqrt(R_simu_t[:, :, 0, k ]))
    sigma_simu_sq_t[:, 0, k] = np.diag(Sigma_t_0)
    D_0 = np.diag(np.sqrt(sigma_simu_sq_t[:, 0, k]))
    r_simu_t[:, 0, k] = mu_i + np.matmul(D_0, z_simu_t[:, 0, k])
    for t in range(T_simu - 1):
        Q_simu_t[:, :, t+1, k] = Q_bar*(1 - DCCalpha - DCCbeta)\
        + DCCalpha*np.outer(z_simu_t[:, t, k], z_simu_t[:, t, k])\
        + DCCbeta*Q_simu_t[:, :, t, k]
        DQsqrt = np.diag(np.sqrt(np.diag(Q_simu_t[:, :, t+1, k])) ** -1)
        R_simu_t[:, :, t+1, k] = np.matmul(np.matmul(DQsqrt, Q_simu_t[:, :, t+1, k]), DQsqrt)
        z_simu_t[:, t+1, k] = np.matmul(np.random.randn(5), np.sqrt(R_simu_t[:, :, t+1, k]))
        for i in range(N):
            sigma_simu_sq_t[i, t+1, k] = Garchomega[i] \
            + Garchalpha[i] * (r_simu_t[i, t, k])**2 \
            + Garchbeta[i] * sigma_simu_sq_t[i, t, k]
        D = np.diag(np.sqrt(sigma_simu_sq_t[:, t+1, k]))
        r_simu_t[:, t+1, k] = mu_i + np.matmul(D, z_simu_t[:, t+1, k])
        Sigma_simu_t[:, :, t+1, k] = np.matmul(np.matmul(D, R_simu_t[:, :, t+1, k]), D)