# finding Black Scholes implied volatility

I am trying to imply BS volatilities using Yahoo Finance! API (MSFT options) and Reuters zero curves. I am ignoring dividends for the moment. Note that I've also cleant the options' data (e.g., I make sure that last trade day of the option was the day on which I am valuing, and to get a nice smile I've only kept OTM calls and OTM puts under the assumption that OTM options are chea

def f_call(vol, S, r, K, dt, q):
d1 = 1 / vol * np.sqrt(dt) * (np.log(S/K) + (r - q + vol**2 / 2) * dt)
d2 = d1 - vol * np.sqrt(dt)
return np.exp(-q*dt) * S * stats.norm.cdf(d1) - np.exp(-r*dt) * K * stats.norm.cdf(d2)

def f_put(vol, S, r, K, dt, q):
d1 = 1 / vol * np.sqrt(dt) * (np.log(S/K) + (r - q + vol**2 / 2) * dt)
d2 = d1 - vol * np.sqrt(dt)
return stats.norm.cdf(-d2) * K * np.exp(-r*dt) - stats.norm.cdf(-d1) * S * np.exp(-q*dt)

vol_guess = self.vol_yh[i]
print(i)
params = {
'S' : self.S,
'K' : self.K[i],
'dt' : self.dt[i],
'r' : self.r[i],
'q' : self.q[i],
}
option_type = self.type[i]
if option_type == 'C':
f = self.f_call
else:
f = self.f_put

bounds = optimize.Bounds(lb=0, ub=100)
e = lambda s: np.sqrt(np.mean((f(s, **params) - market_price[i]) ** 2))
# e = lambda s: f(s, **params) - market_price[i]
res = optimize.minimize(fun=e, x0=(vol_guess,), method='SLSQP', bounds=bounds)
res = res.x
# res = optimize.newton(func=e, x0=vol_guess)
# print(res)
print(f'Yahoo volatility: {self.vol_yh[i]}')
print(f'Model volatility: {res}')
self.vol[i] = res


In this example, params = {'S': 268.75, 'K': 80.0, 'dt': 0.28055555555555556, 'r': 0.015503659068441186, 'q': 0.0}, where dt is 'ACT/365' and 'r' is found by interpolating linearly on dt.

For reference, I am posting also a snapshot of the zero rates ('deltaTime' is also ACT/365). These come from RIC ticker USDSROISZ=, which are zero rates boostrapped from USDSROIS (SOFR OIS par rates), which I think is good enough to get sensible results. As you can see the rate for my option is between 5 and 6.

deltaTime  zeroRates
0    0.002778   0.007940
1    0.005556   0.008077
2    0.019444   0.009127
3    0.083333   0.011614
4    0.169444   0.013341
5    0.255556   0.015036
6    0.508333   0.019158
7    0.758333   0.022491
8    1.013889   0.024985
9    1.269444   0.026642
10   1.522222   0.027639
11   1.775000   0.028254
12   2.030556   0.028567
13   2.286111   0.028716


However, there is something terribly wrong going on since the results I am getting are way different than the quoted Yahoo Finance! volatility.

Yahoo volatility: 0.7812521875
Model volatility: 8.279821273049492e-12


Interestingly, if I move toward the money, I get sensible results:

params = {'S': 268.75, 'K': 250.0, 'dt': 0.28055555555555556, 'r': 0.015503659068441186, 'q': 0.0}
Yahoo volatility: 0.3357610369873046
Model volatility: 0.3448098337803358


And finally, if I run the whole continuum of options for that particular expiry, I get the following volatility smile plots:  I understand that I can't match it perfectly, also dividends are playing a part (which I need to add), Yahoo could also have a different model, however I don't understand why the deep OTM and ITM options are SO different. Note that I've also cleant the options' data (e.g., I make sure that last trade day of the option was the day on which I am valuing, and to get a nice smile I've only kept OTM calls and OTM puts under the assumption that OTM options are cheaper and therefore more liquid). But, in my opinion, none of this can explain the huge differences in the deep moneyness. What is happening is that the minimization hits the lower bound so it's very close to 0, but then again, why would the minimization want to move to negative volatility?

EDIT:

Fixed the typo in d1, and I now get a much more sensible smile, however I have a jump ATM. It might be because I use puts on the right and calls on the left, the OTM calls are actually overpriced because dividends are not taken into account I believe. • I once had a similar problem. Then it was the normal distribution function in the package I used that was ill-behaved in the tails. It gave non-smooth derivatives which gave the equation solving algorithm a hard time. I like to remember it that I substituted the function in its tails for a term from its taylor expansion and that solved it. Jun 9 at 13:09
• Ah, you approximated the normal CDF? Jun 9 at 13:29
• yes, because the built-in package didn't produce precise enough results near 0 and 1, the approximations are in the wikipedia Jun 9 at 13:57

I think you get it wrong your definition of $$d_1$$. Shouldn't it be $$d_1 = 1 / (\text{vol} * \text{dt}) * (\log(S/K) + (r - q + \text{vol}^2 / 2) * dt)$$ and not $$d_1 = 1 / \text{vol} * \text{dt} * (\log(S/K) + (r - q + \text{vol}^2 / 2) * dt)$$ ?