为何该Python代码计算的欧式看涨期权值与BS公式结果不符？

首先，先修正你提供的Python代码中的语法错误（导入模块的语句需要分开写），正确的代码如下：

import math
import numpy as np

S0 = 100.
K = 100.
T = 1.0
r = 0.05
sigma = 0.2

M = 100
dt = T / M
I = 500000

S = np.zeros((M + 1, I))
S[0] = S0

for t in range(1, M + 1):
    z = np.random.standard_normal(I)
    S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)

C0 = math.exp(-r * T) * np.sum(np.maximum(S[-1] - K, 0)) / I
print("European Option Value is ", C0)

接下来解答你的疑问：你提到模拟结果约为10.45，但认为Black-Scholes（BS）公式的结果应为10.09，这其实是你对BS公式的计算出现了错误——正确的BS公式计算结果应该约为10.46，和蒙特卡洛模拟的结果几乎一致。

详细验证BS公式计算过程 #

欧式看涨期权的BS公式为：
$$C_0 = S_0 N(d_1) - K e^{-rT} N(d_2)$$
其中：
$$d_1 = \frac{\ln(S_0/K) + (r + \sigma^2/2)T}{\sigma \sqrt{T}}$$
$$d_2 = d_1 - \sigma \sqrt{T}$$

代入你的参数：

• $S_0=100, K=100, T=1, r=0.05, \sigma=0.2$

1. 计算$\ln(S_0/K) = \ln(1) = 0$
2. 分子部分：$0 + (0.05 + 0.5*0.2^2)*1 = 0.05 + 0.02 = 0.07$
3. $d_1 = 0.07/(0.2*1) = 0.35$，查标准正态分布表得$N(0.35)≈0.6368$
4. $d_2 = 0.35 - 0.2*1 = 0.15$，查标准正态分布表得$N(0.15)≈0.5596$
5. 折现因子$e^{-rT} = e^{-0.05*1}≈0.9512$
6. 代入公式：$C_0 = 1000.6368 - 1000.9512*0.5596 ≈ 63.68 - 53.22 = 10.46$

为什么你的模拟结果和正确BS值接近？ #

你的蒙特卡洛代码逻辑是完全正确的：

• 采用了几何布朗运动的欧拉离散化方法模拟股票价格路径
• 对到期日的股票价格与执行价的差值取最大值后折现，再取平均值得到期权价值

当模拟次数$I=500000$时，结果已经非常接近BS解析解，微小的差异（10.45 vs 10.46）是蒙特卡洛模拟固有的随机性误差，随着模拟次数进一步增加，结果会逐渐收敛到10.46左右。

你可能的BS计算错误点 #

你得到10.09的错误结果，大概率是以下某一步出错：

• 误将连续复利折现因子$e^{-rT}$当成了单利折现$1/(1+rT)$
• 计算$d_1$时误用了$(r - \sigma^2/2)$而非$(r + \sigma^2/2)$
• 查标准正态分布表时取错了$N(d_1)$或$N(d_2)$的值
• 参数输入错误（比如误将$\sigma$设为0.18而非0.2）

内容的提问来源于stack exchange，提问作者M Thomas