Implied Volatility Solver and Vega

Article Quant Q&A · Author: Joan Arau

Summary

The document presents a Python Newton-style procedure for estimating implied volatility from an option’s market price using Black–Scholes price and vega calculations. It describes a failure mode in which iterations become unstable when vega is very small, especially for options away from the money, and asks whether using a model for American options would help. The example also calculates option Greeks after estimating volatility, with a binomial engine available when a dividend yield is supplied.

The accepted response identifies a unit scaling issue: the stated vega is expressed per one-unit volatility change, while the iteration treats volatility in percentage-point terms, so the update is divided by 100. This addresses the reported overshooting attributed to the mismatch. The note does not establish that this adjustment makes the solver robust for every strike or market price. Newton iterations can still be unreliable when sensitivity is low or inputs are inconsistent, and the example’s model choice does not by itself resolve the numerical behavior.

Key ideas

The example estimates implied volatility by iteratively matching a model price to the market option price.
Newton updates depend on vega, so a small vega can make volatility estimates unstable.
The accepted response attributes the issue to a factor-of-100 mismatch in vega scaling.
Changing the option pricing model to a binomial model does not by itself address the solver’s scaling or convergence behavior.

Tags

Options Volatility Derivatives pricing Statistics

Full text

# Solving for Implied Volatility Vega gets stuck at 0 (Python)


# Solving for Implied Volatility Vega gets stuck at 0 (Python)












So my goal is to calculate option greeks with as few manual inputs as possible. I managed to get the IV for at the money options but then when I try further OTM strikes my results get completely screwed.

Im running this loop to estimate IV:

```
def find_vol(target_value, call_put, S, K, T, r):
MAX_ITERATIONS = 200
PRECISION = 0.01

sigma = 0.01
s = {}
for i in range(0, MAX_ITERATIONS):
    price = bs_price(call_put, S, K, T, r, sigma)
    vega = bs_vega(call_put, S, K, T, r, sigma)

    price = price
    diff = target_value - price  # our root
    print({'sigma':sigma,'vega':vega,'price':price, 'target_value':target_value,'diff':diff})

    if abs(diff) < PRECISION:
        return sigma
    sigma = sigma + diff/vega
    s[diff]=sigma

mindiff =min(list(s.keys()),key=abs)
sigma = s[mindiff]
# value wasn't found, return best guess so far
return sigma
```

My issue is that with certain options the loop will get stuck because Vega goes to 0 and then this line doesn't make sense anymore: sigma = sigma + diff/Vega when the problem happens it looks like this :

The issue is probably due to the fact that I am using BS model on American options which is clearly wrong but I am unsure if my results will improve when using Binomial. Suggestions would be much appreciated!

Full working Code: '''python

```
import datetime
from scipy.stats import norm
from math import exp, log,sqrt
import quantsbin.derivativepricing as qbdp

def find_vol(target_value, call_put, S, K, T, r):
    MAX_ITERATIONS = 200
    PRECISION = 0.01

    sigma = 0.5
    s = {}
    for i in range(0, MAX_ITERATIONS):
        price = bs_price(call_put, S, K, T, r, sigma)
        vega = bs_vega(call_put, S, K, T, r, sigma)

        price = price
        diff = target_value - price  # our root
        #print({'sigma':sigma,'vega':vega,'price':price,   'target_value':target_value,'diff':diff})

        if abs(diff) < PRECISION:
            return sigma

        sigma = sigma + diff/vega # f(x) / f'(x)
        s[diff]=sigma

    mindiff =min(list(s.keys()),key=abs)
    #print(mindiff)
    sigma = s[mindiff]
    #print(sigma)
    # value wasn't found, return best guess so far
    return sigma

n = norm.pdf
N = norm.cdf

def bs_price(cp_flag,S,K,T,r,v,q=0.0):
    d1 = (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    d2 = d1-v*sqrt(T)
    if cp_flag == 'c':
        price = S*exp(-q*T)*N(d1)-K*exp(-r*T)*N(d2)
    else:
        price = K*exp(-r*T)*N(-d2)-S*exp(-q*T)*N(-d1)
    return price

def bs_vega(cp_flag,S,K,T,r,v,q=0.0):
    d1 = (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    return S * sqrt(T)*n(d1)

def     calculate_iv_greeks(price_market,strike,expiration,date_calc,spot_underlying,rfr,kind,div_yield=0):
exp_y,exp_m,exp_d = expiration
calc_y, calc_m, calc_d = date_calc

expiration = datetime.date(exp_y, exp_m, exp_d)
date_calc = datetime.date(calc_y, calc_m, calc_d)

T = (expiration - date_calc).days / 365.
#T = (expiration - date_calc).days
r = ((rfr/100)/365)*(expiration - date_calc).days
iv = find_vol(price_market,kind,spot_underlying,strike,T,r)
if kind =='C': type = 'Call'
else: type = 'Put'

expiration = datetime.date(exp_y,exp_m,exp_d).strftime("%Y%m%d")
date_calc = datetime.date(calc_y, calc_m, calc_d).strftime("%Y%m%d")

equity_option1 = qbdp.EqOption(option_type=type, strike=strike, expiry_date=str(expiration))
if div_yield != 0:
    eq1_engine = equity_option1.engine(model='Binomial', pricing_date=str(date_calc), spot0=spot_underlying, rf_rate=rfr/100,
                                   volatility=iv, yield_div=div_yield)
    print('Engine:Binomial')
else:
    eq1_engine = equity_option1.engine(model='BSM', pricing_date=str(date_calc), spot0=spot_underlying, rf_rate=rfr/100,
                                   volatility=iv)
    print('Engine:BSM')
D,G,TH,V,R,P = eq1_engine.risk_parameters().values()
print(eq1_engine.risk_parameters())
return eq1_engine.valuation(),iv, D,TH

V_market = 5.85
K = 285
S = 289.85
cp = 'c'

greeks=calculate_iv_greeks(price_market=V_market,strike=K,expiration=. [2019,5,1],date_calc=[2019,4,17],
                   spot_underlying=S,rfr=2.41,kind=cp,div_yield=0)

print(greeks)

print('Broker values: ', {'IV':0.1217 ,'Delta':0.765 })
```

'''

## Answer by Charles Fox (score 1, accepted)

https://quant.stackexchange.com/a/45258

Your vega is off by a factor of 100. Change it to sigma = sigma + diff/vega/100.

```
import datetime
from scipy.stats import norm
from math import exp, log,sqrt
import quantsbin.derivativepricing as qbdp

def find_vol(target_value, call_put, S, K, T, r):
    MAX_ITERATIONS = 200
    PRECISION = 0.01

    sigma = 0.5
    s = {}
    for i in range(0, MAX_ITERATIONS):
        price = bs_price(call_put, S, K, T, r, sigma)
        vega = bs_vega(call_put, S, K, T, r, sigma)

        price = price
        diff = target_value - price  # our root
        #print({'sigma':sigma,'vega':vega,'price':price,   'target_value':target_value,'diff':diff})

        if abs(diff) < PRECISION:
            return sigma

        sigma = sigma + diff/vega/100 # f(x) / f'(x)
        s[diff]=sigma

    mindiff =min(list(s.keys()),key=abs)
    #print(mindiff)
    sigma = s[mindiff]
    #print(sigma)
    # value wasn't found, return best guess so far
    return sigma

n = norm.pdf
N = norm.cdf

def bs_price(cp_flag,S,K,T,r,v,q=0.0):
    d1 = (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    d2 = d1-v*sqrt(T)
    if cp_flag == 'c':
        price = S*exp(-q*T)*N(d1)-K*exp(-r*T)*N(d2)
    else:
        price = K*exp(-r*T)*N(-d2)-S*exp(-q*T)*N(-d1)
    return price

def bs_vega(cp_flag,S,K,T,r,v,q=0.0):
    d1 = (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    return S * sqrt(T)*n(d1)

def calculate_iv_greeks(price_market,strike,expiration,date_calc,spot_underlying,rfr,kind,div_yield=0):
    exp_y,exp_m,exp_d = expiration
    calc_y, calc_m, calc_d = date_calc

    expiration = datetime.date(exp_y, exp_m, exp_d)
    date_calc = datetime.date(calc_y, calc_m, calc_d)

    T = (expiration - date_calc).days / 365
    #T = (expiration - date_calc).days
    r = ((rfr/100)/365)*(expiration - date_calc).days
    iv = find_vol(price_market,kind,spot_underlying,strike,T,r)
    if kind =='C': type = 'Call'
    else: type = 'Put'

    expiration = datetime.date(exp_y,exp_m,exp_d).strftime("%Y%m%d")
    date_calc = datetime.date(calc_y, calc_m, calc_d).strftime("%Y%m%d")

    equity_option1 = qbdp.EqOption(option_type=type, strike=strike, expiry_date=str(expiration))
    if div_yield != 0:
        eq1_engine = equity_option1.engine(model='Binomial', pricing_date=str(date_calc), spot0=spot_underlying, rf_rate=rfr/100,
                                       volatility=iv, yield_div=div_yield)
        print('Engine:Binomial')
    else:
        eq1_engine = equity_option1.engine(model='BSM', pricing_date=str(date_calc), spot0=spot_underlying, rf_rate=rfr/100,
                                       volatility=iv)
        print('Engine:BSM')

    D,G,TH,V,R,P = eq1_engine.risk_parameters().values()
    print(eq1_engine.risk_parameters())
    return eq1_engine.valuation(),iv, D,TH

V_market = 5.85
K = 285
S = 289.85
cp = 'c'

expiration = datetime.date(2019, 5, 1)
date_calc = datetime.date(2019, 4, 17)

sig = find_vol(V_market, cp, S, K, T = (expiration - date_calc).days / 365, r=.0241)
print("Estimated implied vol: " + str(np.round(sig*100,2)))
```

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.

Newton Implied Volatility Solver and Vega Scaling

Summary

Key ideas

Tags