如何基于导数定义在R中用蒙特卡洛模拟计算期权一阶Greeks值

实现思路 #

一阶Greeks本质是期权价格对对应参数的一阶偏导数，用有限差分法（即你老师提到的导数定义实现）可以直接适配你已经写好的蒙特卡洛、二叉树定价函数，完全不需要依赖Black-Scholes模型的解析公式。
首先你需要把已有的定价逻辑封装成独立函数，输入参数至少包含：标的价格S、行权价K、剩余期限T、无风险利率r、波动率σ、分红率q、期权类型（看涨/看跌），输出为对应的期权价格V，函数结构参考：

# 蒙特卡洛定价函数，替换内部逻辑为你自己写的代码即可
mc_option_price <- function(S, K, T, r, sigma, q, type = "call") {
  # 你的蒙特卡洛定价逻辑
  # 最后返回期权价格
}
# 二叉树定价函数
bt_option_price <- function(S, K, T, r, sigma, q, type = "call") {
  # 你的二叉树定价逻辑
  # 最后返回期权价格
}

各一阶Greeks的有限差分计算规则 #

h取极小值即可，通常取1e-4到1e-6量级，用中心差分比单边差分精度更高：

• Delta：期权价格对标的价格的偏导数，计算式为 (V(S+h, K, T, r, sigma, q) - V(S-h, K, T, r, sigma, q)) / (2*h)，h可取标的价格的0.01%
• Vega：期权价格对波动率的偏导数，计算式为 (V(S, K, T, r, sigma+h, q) - V(S, K, T, r, sigma-h, q)) / (2*h)，h可取1e-4
• Theta：期权价格对剩余期限的偏导数，取负值代表时间价值衰减，计算式为 -(V(S, K, T+h, r, sigma, q) - V(S, K, T-h, r, sigma, q)) / (2*h)，h可取1/365（1天的年化时长）
• Rho：期权价格对无风险利率的偏导数，计算式为 (V(S, K, T, r+h, sigma, q) - V(S, K, T, r-h, sigma, q)) / (2*h)，h可取1e-4
• Lambda：期权杠杆率，等于Delta乘以标的价格再除以期权价格，计算式为 delta * S / V
• Epsilon：期权价格对分红率的偏导数，计算式为 (V(S, K, T, r, sigma, q+h) - V(S, K, T, r, sigma, q-h)) / (2*h)，h可取1e-4

R完整实现示例 #

# 示例蒙特卡洛定价函数，你可以替换成自己的定价逻辑
mc_option_price <- function(S, K, T, r, sigma, q, type = "call") {
  set.seed(123) # 固定随机种子保证差分计算时路径一致，避免随机误差干扰
  n_sim <- 10000
  ST <- S * exp((r - q - 0.5*sigma^2)*T + sigma*sqrt(T)*rnorm(n_sim))
  payoff <- ifelse(type == "call", pmax(ST - K, 0), pmax(K - ST, 0))
  price <- exp(-r*T) * mean(payoff)
  return(price)
}

# 封装Greeks计算函数，切换定价模型只需要替换price_fun参数即可
calc_greeks <- function(S, K, T, r, sigma, q, type = "call", price_fun = mc_option_price) {
  # 基准期权价格
  V0 <- price_fun(S, K, T, r, sigma, q, type)
  
  # 计算Delta
  h_S <- 1e-4 * S
  V_plus_S <- price_fun(S + h_S, K, T, r, sigma, q, type)
  V_minus_S <- price_fun(S - h_S, K, T, r, sigma, q, type)
  delta <- (V_plus_S - V_minus_S) / (2 * h_S)
  
  # 计算Vega
  h_sigma <- 1e-4
  V_plus_sigma <- price_fun(S, K, T, r, sigma + h_sigma, q, type)
  V_minus_sigma <- price_fun(S, K, T, r, sigma - h_sigma, q, type)
  vega <- (V_plus_sigma - V_minus_sigma) / (2 * h_sigma)
  
  # 计算Theta
  h_T <- 1/365
  V_plus_T <- price_fun(S, K, T + h_T, r, sigma, q, type)
  V_minus_T <- price_fun(S, K, T - h_T, r, sigma, q, type)
  theta <- -(V_plus_T - V_minus_T) / (2 * h_T)
  
  # 计算Rho
  h_r <- 1e-4
  V_plus_r <- price_fun(S, K, T, r + h_r, sigma, q, type)
  V_minus_r <- price_fun(S, K, T, r - h_r, sigma, q, type)
  rho <- (V_plus_r - V_minus_r) / (2 * h_r)
  
  # 计算Lambda
  lambda <- delta * S / V0
  
  # 计算Epsilon
  h_q <- 1e-4
  V_plus_q <- price_fun(S, K, T, r, sigma, q + h_q, type)
  V_minus_q <- price_fun(S, K, T, r, sigma, q - h_q, type)
  epsilon <- (V_plus_q - V_minus_q) / (2 * h_q)
  
  return(list(
    price = V0,
    delta = delta,
    vega = vega,
    theta = theta,
    rho = rho,
    lambda = lambda,
    epsilon = epsilon
  ))
}

# 测试调用
calc_greeks(S = 100, K = 100, T = 1, r = 0.03, sigma = 0.2, q = 0.01, type = "call")

注意事项 #

• 用蒙特卡洛定价计算Greeks时，必须固定随机种子，保证计算f(x+h)和f(x-h)时使用同一批随机路径，否则随机误差会覆盖差分的有效信号，结果完全失真
• 二叉树计算时不需要处理随机数问题，但要注意h的取值不要太小，避免因二叉树的离散截断误差导致结果异常
• 要切换为二叉树模型计算Greeks，只需要把calc_greeks的price_fun参数替换为你自己的二叉树定价函数即可，其他逻辑不需要修改

内容的提问来源于stack exchange，提问作者Gaia Mancassola