以程式方式推導隱含波動率

我正在做一個項目來計算使用Python. 我正在使用 Black-Scholes 模型，我可以通過插入隱含波動率的給定值來獲得準確的結果。我通常從我擁有的教科書或optionsprofitcalculator中獲得 IV 。

我怎樣才能自己推導出IV？到目前為止，我發現了兩種方法：

1. 這個答案說IV可以通過使用BSM迭代得出。我也在其他答案中讀到過這個，但我認為我沒有正確理解它。似乎它說要插入猜測並進行調整，直到找到為該選項提供正確值的 IV。但是要辨識正確的 IV，我需要已經有了答案，比如教科書或 OPC 中的答案，對吧？在不檢查其他來源的情況下，您將如何以程式方式實現這一點？
2. 這個問題說您可以從無模型變異數交換中得出 IV。與迭代方法相比，它是否更適合我的項目？如果是這樣，我將如何以程式方式在實踐中實現它？

近似歐式期權的隱含波動率可以通過幾種方式完成——這只是其中之一。下面是一個使用Newton Raphson的 python 實現。您可以使用該implied_volatility函式找到近似的隱含波動率。然後，您可以通過將輸出插入回option_price函式來檢查它。

import numpy as np
from scipy.stats import norm

"""
right = 'C' or 'P'
s = Spot Price
k =  Strike Price
t = Days to expiration
rfr = Risk-free Rate
sigma = volatility
div = Annual dividend rate. Defaulted to zero.
price = Known option price. Needed for implied_volatility function
"""


def d_one(s, k, t, rfr, sigma, div=0):
   """d1 calculation"""
   d_1 = (np.log(s / k) +
          (rfr - div + sigma ** 2 / 2) * t) / (sigma * np.sqrt(t))
   return d_1


def d_two(s, k, t, rfr, sigma, div=0):
   """d2 calculation"""
   d_2 = d_one(s, k, t, rfr, sigma, div) - sigma * np.sqrt(t)
   return d_2


def nd_one(right, s, k, t, rfr, sigma, div=0):
   """nd1 calculation"""
   if right == 'C':
       nd_1 = norm.cdf(d_one(s, k, t, rfr, sigma, div), 0, 1)
   elif right == 'P':
       nd_1 = norm.cdf(-d_one(s, k, t, rfr, sigma, div), 0, 1)
   return nd_1


def nd_two(right, s, k, t, rfr, sigma, div=0):
   """nd2 calculation"""
   if right == 'C':
       nd_2 = norm.cdf(d_two(s, k, t, rfr, sigma, div), 0, 1)
   elif right == 'P':
       nd_2 = norm.cdf(-d_two(s, k, t, rfr, sigma, div), 0, 1)
   return nd_2


def option_price(right, s, k, t, rfr, sigma, div=0):
   """option price"""
   right = right.upper()
   t /= 365
   if right == 'C':
       price = (s * np.exp(-div * t) *
                nd_one(right, s, k, t, rfr, sigma, div)
                - k * np.exp(-rfr * t) *
                nd_two(right, s, k, t, rfr, sigma, div))
   elif right == 'P':
       price = (k * np.exp(-rfr * t) *
                nd_two(right, s, k, t, rfr, sigma, div)
                - s * np.exp(-div * t) *
                nd_one(right, s, k, t, rfr, sigma, div))
   return price


def option_vega(s, k, t, rfr, sigma, div=0):
   """option vega"""
   t /= 365
   vega = (.01 * s * np.exp(-div * t) * np.sqrt(t)
           * norm.pdf(d_one(s, k, t, rfr, sigma, div)))
   return vega


def implied_volatility(right, s, k, t, rfr, price, div=0):
   """implied volatility approximation"""
   epsilon = 0.00000001
   sigma = 1.0

   def newton_raphson(right, s, k, t, rfr, sigma, price, epsilon, div=0):
       diff = option_price(right, s, k, t, rfr, sigma, div) - price
       while diff > epsilon:
           sigma = (sigma -
                    (option_price(right, s, k, t, rfr, sigma, div) - price) /
                    (option_vega(s, k, t, rfr, sigma, div) * 100))
           diff = np.abs(
                   option_price(right, s, k, t, rfr, sigma, div) - price)
       return sigma

   iv = newton_raphson(right, s, k, t, rfr, sigma, price, epsilon, div)
   return iv

引用自：https://quant.stackexchange.com/questions/49669