# 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 })


'''

• Please create an MVCE. This code is missing definitions. Commented Apr 17, 2019 at 17:43
• Thanks for the comment, I added the full code Commented Apr 17, 2019 at 18:10

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)))