我用50ETF期权，1月28日收盘数据，来计算一下，上面的为解析解，vega = 0.19922296046043902
下面是数值解，vega = 0.320760687225449,这个与期权软件上显示的vega值一致。 为什么两个会相差这么大，解析解的公式有误么？

def calculate_original_vega(
s: float,
k: float,
r: float,
t: float,
v: float,
d1: float = 0.0
) -> float:
"""Calculate option vega"""
if v <= 0:
return 0

if not d1:
    d1: float = calculate_d1(s, k, r, t, v)

vega: float = s * pdf(d1) * sqrt(t)

return vega

对于解析解，我直接用的C:\vnstudio\Lib\site-packages\vnpy_optionmaster\pricing，这个目录下的源文件，并不是我学习的时候手打出来的

提供下你的数值解计算函数代码？

from math import log, pow, sqrt, exp
from scipy import stats

cdf = stats.norm.cdf  # 累积概率分布


def calculate_d1(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float
) -> float:
    """Calculate option D1 value"""
    d1: float = (log(s / k) + (r + 0.5 * pow(v, 2)) * t) / (v * sqrt(t))
    return d1


def calculate_price(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float,
    cp: int
) -> float:
    """Calculate option price"""
    # Return option space value if volatility not positive
    if v <= 0:
        return max(0, cp * (s - k))

    d1: float = calculate_d1(s, k, r, t, v)
    d2: float = d1 - v * sqrt(t)

    price: float = cp * (s * cdf(cp * d1) - k * cdf(cp * d2) * exp(-r * t))
    return price


def calculate_delta(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float,
    cp: int,
    dx: float = 1/1000
) -> float:
    """"""
    y_up = calculate_price(s * (1+dx), k, r, t, v, cp)
    y_down = calculate_price(s * (1-dx), k, r, t, v, cp)
    y_change = y_up - y_down
    x_change = s * dx * 2
    delta = y_change / x_change
    return delta


def calculate_vega(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float,
    cp: int,
    dx: float = 1/1000
) -> float:
    """"""
    y_up = calculate_price(s, k, r, t, v * (1+dx), cp)
    y_down = calculate_price(s, k, r, t, v * (1-dx), cp)
    y_change = y_up - y_down
    x_change = v * dx * 2
    vega = y_change / x_change
    return vega


def calculate_theta(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float,
    cp: int,
    dx: float = 1/1000
) -> float:
    """"""
    y_up = calculate_price(s, k, r, t, v, cp)
    y_down = calculate_price(s, k, r, t * (1-dx), v, cp)
    y_change = y_up - y_down
    x_change = t * dx
    theta = y_change / x_change
    return theta


def calculate_gamma(
    s: float,
    k: float,
    r: float,
    t: float,
    v: float,
    cp: int,
    dx: float = 1/1000
) -> float:
    """"""
    y_up = calculate_delta(s * (1+dx), k, r, t, v, cp)
    y_down = calculate_delta(s * (1-dx), k, r, t, v, cp)
    y_change = y_up - y_down
    x_change = s * dx * 2
    gamma = y_change / x_change

    return gamma