蒙特卡羅
多變數 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