我正在尝试运行蒙特卡罗模拟来帮助我在 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)来获得更合理的数字。但我仍然不确定实际的公式。
如果计算 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 计算的建议适用于看跌期权情况) .不过,请注意,我对此并不确定,因为自从我上次涉足这类东西以来已经有一段时间了,并且懒得去仔细检查。我只是从这里得到了这种感觉。)