多變數 MC：我做錯了什麼？

我正在嘗試生成本文中介紹的多元 MC 結果A Simple Generalization of Kirk’s Approximation for Multi-Asset Spread Options by the Lie-Trotter Operator Splitting Method，作者 Chi-Fai Lo https://file.scirp.org/pdf/ JMF_2014050615380663.pdf。

我使用不同的方法進行了不同的嘗試，但最終總是產生相同的錯誤結果。例如，我生成表 5.1 中第一個元素的程式碼生成的13.7424 +/- 0.0045不是13.5763 ± 0.0089.

下面是我製作的 Python 程式碼

如果有人好心告訴我我做錯了什麼……

（有關資訊，我可以確認論文中的 EK 值是正確的，因為我可以找到 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

基於 Quantuple 評論（謝謝），我修復了許多錯誤，並提出了以下程式碼：

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

引用自：https://quant.stackexchange.com/questions/55636