La descomposición Cholesky reduce la volatilidad de los procesos de Wiener simulados / movimientos brownianos

Question

Estoy tratando de simular $n$ correlacionado movimientos geométricos brownianos (GBM) dada una matriz de correlación especificada $\Sigma$ al seguir este procedimiento que utiliza la descomposición Cholesky.

Sin embargo, cuando implemento el código en Python uno de los elementos altamente correlacionados ( $\rho=0.99$ ) Los Procesos de Wiener casi no tienen volatilidad, como se puede ver en el gráfico de abajo en la línea roja.

Yo esperaría una correlación de $\rho=1$ para implicar que los dos procesos serían idénticos y para $\rho=-1$ yo esperaría un "espejo" alrededor de 0.

El gráfico muestra la realización de $W$ a través de $t$ .

Nota: En mi intento trato de evitar el uso de bucles.

# PARAMETERS
n = 2                                              # number of assets
T = 5                                              # number of years
N = 252                                            # number of steps pr. year (T)
r = np.array([0.03, 0.03])                         # drift-rate
q = np.array([0.0, 0.0])                           # divivend-rate
S0 = np.array([1.0, 1.0])                          # initial value
sigma = np.array([0.2, 0.2])                       # diffusion-coefficients 
corr_mat = np.array([[1.0, 0.99], [0.99, 1.0]])    # correlation-matrix

# Calculate step size
dt = 1 / (T * N)

# Vector of time index
t = np.array(range(N*T+1)).reshape((N * T + 1, 1)) * dt

# == PERFORM SIMULATION ==
# Draw normal distributed random variables
eps = np.random.normal(loc=0, scale=np.sqrt(dt), size=(N * T, n))

# Use Cholesky decomposition to get the matrix R satisfying: corr_mat = L x L^*
L = np.linalg.cholesky(corr_mat)

# Transform into realizations of correlated Wiener processes
W = eps @ L

# Calculate realizations of each stock as the realization of a GBM
# S_t = S0 * exp{ ((r-q) - ½sigma)*dt + sigma * W}
S = np.exp(
    (((r - q) - np.power(sigma, 2) / 2) * dt).T + W * sigma.T
)
S = np.vstack([np.ones(S.shape[1]), S])
S = S0 * S.cumprod(axis=0)

Answer 1

Arreglar la siguiente parte

# Transform into realizations of correlated Wiener processes
W = eps @ L.T

Originalmente, debería ser L @ z donde z es un vector normal multivariado. Como cada fila de eps es un vector de este tipo, eps @ L.T es correcto.

Obtendrá el resultado deseado como se indica a continuación: