我的蒙特卡罗模拟不断给出错误的 Theta 值。其他估计选项希腊语是正确的，但 Theta 不正确

Question

我正在尝试运行蒙特卡罗模拟来帮助我在 python 中为希腊期权定价。我得到的大多数值都是正确的，但我不确定为什么我从来没有得到正确的 Theta 结果。另外，Theta 的估计值不断变化，但它永远不会围绕正确答案。

正确答案： {“Delta”：-0.3611，“Gamma”：0.0177，“Theta”：-0.745，“Vega”：26.483，“Rho”：-49.635}

我的代码给了我什么： {'德尔塔'：-0.36114866036456306，'伽玛'：0.01765540297071766，'西塔'：4.969082186799323，'织女星'：26.482584628353578，'罗'：-49.64926774 1554404}

import numpy as np
from scipy.stats import norm

def monte_carlo_option_greeks(option_type, S, K, T, r, sigma, steps, N, dS=0.01, dSigma=0.01, dT=0.01):
    """
    Calculate option Greeks using Monte Carlo simulation.

    Inputs:
    - option_type: 'call' or 'put'
    - S: Current stock price
    - K: Strike price
    - T: Time to maturity (in years)
    - r: Risk-free interest rate
    - sigma: Volatility
    - steps: Time steps for the simulation
    - N: Number of simulation trials
    - dS: Increment for calculating Delta (optional, default is 0.01)
    - dSigma: Increment for calculating Vega (optional, default is 0.01)
    - dT: Increment for calculating Theta (optional, default is 0.01)

    Returns a dictionary containing the estimated Greeks: {'Delta', 'Gamma', 'Theta', 'Vega', 'Rho'}
    """

    def geo_paths(S, T, r, sigma, steps, N):
        dt = T / steps
        ST = np.cumprod(np.exp((r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) *
                               np.random.normal(size=(steps, N))), axis=0) * S
        return ST

    def black_scholes_price(S, K, T, r, sigma):
        d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
        d2 = d1 - sigma * np.sqrt(T)
        if option_type == 'call':
            price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
        else:  # option_type == 'put'
            price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
        return price

    # Simulate underlying asset price paths
    paths = geo_paths(S, T, r, sigma, steps, N)

    # Calculate option payoffs
    if option_type == 'call':
        payoffs = np.maximum(paths[-1] - K, 0)
    else:  # option_type == 'put'
        payoffs = np.maximum(K - paths[-1], 0)

    # Discount to present value
    option_price = np.mean(payoffs) * np.exp(-r * T)

    # Calculate Greeks
    delta = (black_scholes_price(S + dS, K, T, r, sigma) - black_scholes_price(S - dS, K, T, r, sigma)) / (2 * dS)
    
    # Corrected gamma calculation
    gamma = (black_scholes_price(S + dS, K, T, r, sigma) - 2 * black_scholes_price(S, K, T, r, sigma) + black_scholes_price(S - dS, K, T, r, sigma)) / (dS**2)
    
    # Corrected theta calculation with negative sign for a European put option
    dT_theta = dT / 365  # Convert to a daily increment for theta calculation
    theta = -(black_scholes_price(S, K, T - dT_theta, r, sigma) - option_price) / dT

    vega = (black_scholes_price(S, K, T, r, sigma + dSigma) - black_scholes_price(S, K, T, r, sigma - dSigma)) / (2 * dSigma)
    rho = (black_scholes_price(S, K, T, r + dT, sigma) - black_scholes_price(S, K, T, r - dT, sigma)) / (2 * dT)

    greeks = {
        'Delta': delta,
        'Gamma': gamma,
        'Theta': theta,
        'Vega': vega,
        'Rho': rho
    }

    return greeks

# Example usage:
option_type = 'put'  # 'call' or 'put'
S = 50  # stock price S_{0}
K = 52  # strike
T = 2  # time to maturity (in years)
r = 0.05  # risk-free interest rate
sigma = 0.3  # annual volatility in decimal form
steps = 100  # time steps
N = 10000  # number of simulation trials

greeks = monte_carlo_option_greeks(option_type, S, K, T, r, sigma, steps, N)
print(greeks)

我最初获得的 Theta 值非常大，因此我使用了更大的增量（dT = 0.001 到 dT = 0.01）来获得更合理的数字。但我仍然不确定实际的公式。

Answer 1

如果计算 theta 作为

T

的偏移量，则使用

dT_theta = 0.01 / 365

的量，即仅

2.739726027e-05

年（因为这是

black_scholes_price

作为其参数

T

所采用的单位，其变化范围为theta 计算），即 T 的相对差异仅为 1.3698630135e-05，那么我的猜测是精度（浮点运算）开始成为问题，因为价格差异太小。这一假设也得到了重复运行

monte_carlo_option_greeks

时获得的巨大的 θ 方差的支持（而其他希腊人没有到非常小的方差）：

print(np.var([monte_carlo_option_greeks(
    option_type, S, K, T, r, sigma, steps, N
)['Theta'] for k in range(100)]))

打印

65.54872356722292

我还建议坚持使用中心方法（即除以变化参数的增量的两倍并在两个方向上偏移）来计算导数，就像您对其他希腊人所做的那样。

如果完成这些调整，即将 theta 计算替换为

theta  = -(black_scholes_price(S, K, T + dT, r, sigma) - black_scholes_price(S, K, T - dT, r, sigma)) / (dT * 2)

我得到了

-0.7453608249772259

的 theta 值，即预期结果。另外，重复计算表明现在 theta 的方差为零。

（另外，如果我没有记错的话，希腊字母的计算也应该考虑期权类型（看跌期权或看涨期权），而所提供的代码中似乎仍然缺少一些东西（即我对 theta 计算的建议适用于看跌期权情况） .不过，请注意，我对此并不确定，因为自从我上次涉足这类东西以来已经有一段时间了，并且懒得去仔细检查。我只是从这里得到了这种感觉。）