MC multivariante: ¿qué estoy haciendo mal?

Question

Estoy tratando de generar los resultados multivariados de MC presentados en este documento Una simple generalización de la aproximación de Kirk para las opciones de diferencial multiactivo mediante el método de división del operador Lie-Trotter por Chi-Fai Lo https://file.scirp.org/pdf/JMF_2014050615380663.pdf .

He hecho diferentes intentos utilizando diferentes enfoques pero siempre termino produciendo el mismo resultado incorrecto. Por ejemplo, mi código para generar el primer elemento de la tabla 5.1 produce 13.7424 +/- 0.0045 en lugar de 13.5763 ± 0.0089 .

Aquí abajo está el código Python que hice

Si alguien tiene la amabilidad de decirme qué estoy haciendo mal ....

(para información puedo confirmar que los valores de EK en el papel son correctos porque pude encontrar los valores proporcionados por Kirk)

import numpy as np
import scipy.stats

nb_simuls = 10_000_000

# parameters from https://file.scirp.org/pdf/JMF_2014050615380663.pdf Table 1 / cell 1
S1, S2, S3 = 50, 60, 150
v1, v2, v3 = 0.3, 0.3, 0.3

rho_12, rho_23, rho_13 = 0.40, 0.20, 0.80

K = 30.0
T = 0.25
r = 0.05

# from spot to forward values
F1 = S1 * np.exp(r*T)
F2 = S2 * np.exp(r*T)
F3 = S3 * np.exp(r*T)

# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3]) * np.sqrt(T)

# determine covariance matrix from correls and volatilities

correl = np.asarray([[1,     rho_12,   rho_13],
                     [rho_12,    1,     rho_23],
                     [rho_13,   rho_23,   1]])

cov = np.diag(vols).dot(correl).dot(np.diag(vols))

# simulate prices
returns = 1 + np.random.default_rng().multivariate_normal((0, 0, 0), cov, nb_simuls)

# determine exercise and value
values = []
for i in range(nb_simuls):
    v = max(0, F3 * returns[i, 2] - F1 * returns[i, 0] - F2 * returns[i, 1] - K)
    values.append(v)

v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"The option value is: {v} +/- {se}")
# The option value is: 13.742449851577431 +/- 0.004475778801668813

Answer 1

Basándome en los comentarios de Quantuple (gracias), he corregido muchos errores y me ha salido el siguiente código:

import numpy as np
import scipy.stats

nb_simuls = 5_000_000

# parameters from https://file.scirp.org/pdf/JMF_2014050615380663.pdf Table 1 / cell 1
S1, S2, S3 = 50, 60, 150
v1, v2, v3 = 0.3, 0.3, 0.3

rho_12, rho_23, rho_13 = 0.40, 0.20, 0.80

K = 30.0
T = 0.25
r = 0.05

# from spot to forward values
F1 = S1 * np.exp(r*T)
F2 = S2 * np.exp(r*T)
F3 = S3 * np.exp(r*T)

# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3])

# determine covariance matrix from correls and volatilities

correl = np.asarray([[1,     rho_12,   rho_13],
                     [rho_12,    1,     rho_23],
                     [rho_13,   rho_23,   1]])

# simulate prices
returns = np.random.default_rng().multivariate_normal((0, 0, 0), correl, nb_simuls)

# determine exercise and value
values = []
for i in range(nb_simuls):
    F3_exp = F3 * np.exp(-0.5 * v3 ** 2 * T + v3 * np.sqrt(T) * returns[i, 2])
    F2_exp = F2 * np.exp(-0.5 * v2 ** 2 * T + v2 * np.sqrt(T) * returns[i, 1])
    F1_exp = F1 * np.exp(-0.5 * v1 ** 2 * T + v1 * np.sqrt(T) * returns[i, 0])

    v = max(0, F3_exp - F1_exp - F2_exp - K)

    values.append(v)

v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"Method 1 (using correlation matrix):   The option value is: {v} +/- {se}")
# The option value is: 13.742449851577431 +/- 0.004475778801668813

# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3]) * np.sqrt(T)

# determine covariance matrix from correls and volatilities

correl = np.asarray([[1,     rho_12,   rho_13],
                     [rho_12,    1,     rho_23],
                     [rho_13,   rho_23,   1]])

cov = np.diag(vols).dot(correl).dot(np.diag(vols))

# simulate prices
returns = np.random.default_rng().multivariate_normal((0, 0, 0), cov, nb_simuls)

# determine exercise and value
values = []
for i in range(nb_simuls):
    F3_exp = F3 * np.exp(-0.5 * np.diag(cov)[2] + returns[i, 2])
    F2_exp = F2 * np.exp(-0.5 * np.diag(cov)[1] + returns[i, 1])
    F1_exp = F1 * np.exp(-0.5 * np.diag(cov)[0] + returns[i, 0])

    v = max(0, F3_exp - F1_exp - F2_exp - K)
    values.append(v)

v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"Method 2 (using covariance matrix):   The option value is: {v} +/- {se}")

# Method 1 (using correlation matrix): The option value is: 13.5738 +/- 0.0066
# Method 2 (using covariance matrix):  The option value is: 13.5803 +/- 0.0066