基于R语言实现Newton-Raphson算法计算隐含波动率出错求助

已完成的实现 #

Black-Scholes期权定价函数 #

BS = function(Flag,St, K, D, r, Ti, sigma) {
  d1  =  (log(St/K) + (r - D + (sigma^2)/2)*Ti) / (sigma*sqrt(Ti))
  d2  =  d1- sigma*sqrt(Ti)
  
  if(Flag == "call") price = St*exp(-D*(Ti)) * pnorm(d1)  - K*exp(-r*Ti)*pnorm(d2)
  if(Flag != "call") price = K*exp(-r*Ti)*pnorm(-d2)-St*exp(-D*Ti)*pnorm(-d1)
  return(price)}
# 测试调用
BS("call",St=505.15, K=500, D=0, r=0.033, Ti=33/250, sigma=0.2)
# 输出：[1] 18.48827

Vega（波动率一阶导数）函数 #

vega_BS = function(St, K, D, r, Ti, sigma){
  d1  =  (log(St/K) + (r - D + (sigma^2)/2)*Ti) / (sigma*sqrt(Ti))
  vega =  St * dnorm(d1) * sqrt(Ti)
  return(round(vega,4))
}

基于uniroot的隐含波动率计算（正确结果） #

通过uniroot求解得到正确隐含波动率0.394：

sig_implied = function(St, K, D,r, Ti,sigma,Market) {
  root_find = function(sigma){
    BS("call",St, K,D,r, Ti, sigma) - Market}
  round(uniroot(root_find, c(0,1))$root,3)
}
Market = 32.4
sig_implied(St=505.15, K=500, r=0.033,D=0, Ti=33/250,sigma=0.2,Market=Market)
# 输出：[1] 0.394

牛顿-拉夫逊算法实现问题 #

第一版错误实现及结果 #

以下代码运行后得到结果22.73685，与预期0.394偏差极大：

ImpliedVolNewton = function(Market,Flag, St, K, Ti, r, D,sigma, tol=0.0001, maxiter = 100) {

  s = 0.3
  not_converged = Ti
  vega = vega_BS(St, K, D, r, Ti, sigma)
  i = 1
  while (not_converged & (i < maxiter)) {
    err = (Market - BS(Flag,St, K, D, r, Ti, sigma) ) 
    s  =  s + err/vega
    not_converged = (abs(err/vega) > tol)
    i = i + 1
  }
  s }

# 调用测试
ImpliedVolNewton(Market=32.4,"call",St=505.15, K=500, Ti=33/250, r=0.033, D=0,sigma=0.2,tol=0.0001)
# 输出：[1] 22.73685

更新版代码仍无法运行 #

修改后的代码依旧无法得到正确结果：

implied_volatility = function(Market,Flag,St,K,Ti,r,D,sigma,tol=0.0001,max_iterations=100){
  sigma0 = sqrt(abs(log(St/K)+r*Ti)*(2/Ti))
  for(i in max_iterations){
    diff = BS(Flag,St,K,Ti,r,D,sigma)-Market
    if(abs(diff)<tol){
      break
    }
    Sigma = sigma0 -diff/vega_BS(St,K,r,D,Ti,sigma)
  }
  return(Sigma)
}
# 调用测试
implied_volatility (Market=32.4,"call",St=505.15, K=500, Ti=33/250, r=0.033, D=0,sigma=0.2,tol=0.0001)

错误原因分析及修正 #

第一版代码核心错误 #

1. 未迭代更新波动率：始终使用初始传入的sigma计算BS价格和Vega，而非每次迭代后的s，导致修正方向完全错误。
2. 收敛条件逻辑混乱：not_converged初始值设为Ti（数值），且收敛条件判断的是修正步长而非价格误差，不符合迭代收敛要求。
3. Vega未动态更新：Vega是波动率的函数，每次迭代都需用当前波动率重新计算，而非提前固定。

更新版代码核心错误 #

1. 循环逻辑错误：for(i in max_iterations)仅执行一次，正确写法应为for(i in 1:max_iterations)。
2. 未迭代更新波动率：始终用初始sigma计算价格和Vega，未将每次迭代得到的Sigma赋值给下一次计算的参数。
3. 参数顺序错误：调用vega_BS时r和D位置颠倒，导致Vega计算错误。

修正后的牛顿-拉夫逊算法代码 #

ImpliedVolNewton_Fixed = function(Market, Flag, St, K, Ti, r, D, init_sigma=0.2, tol=0.0001, maxiter=100) {
  sigma = init_sigma
  for(i in 1:maxiter) {
    bs_price = BS(Flag, St, K, D, r, Ti, sigma)
    diff = bs_price - Market
    if(abs(diff) < tol) break
    vega = vega_BS(St, K, D, r, Ti, sigma)
    # 牛顿迭代公式：sigma_new = sigma_old - f(sigma_old)/f'(sigma_old)
    sigma = sigma - diff / vega
    # 限制波动率在合理范围
    sigma = max(0.0001, min(sigma, 2))
  }
  return(round(sigma, 3))
}

# 测试调用
ImpliedVolNewton_Fixed(Market=32.4, "call", St=505.15, K=500, Ti=33/250, r=0.033, D=0, init_sigma=0.2)
# 输出：[1] 0.394

内容的提问来源于stack exchange，提问作者Homer Jay Simpson