# Calculate implied volatility for options on stocks and futures with two models

I wrote a code for a quant finance job and they told me that, besides it worked, it was poorly written. I asked for a more detailed feedback but they did not send it to me. I paste it here with minimal info about it, so it will be difficult to connect it to the company.

The scope of the code is to calculate implied volatility for options on two different underlyings (stocks, futures) with two different models (Black and Scholes and another one, for which they gave me some publications).

They asked to write all the math functions from scratch and do not use third party libraries.

What do you suggest in order to improve?

import csv
from math import *

def cdf_of_normal(x): #reference Paul Wilmott on Quantitative Finance

a1 = 0.31938153
a2 = -0.356563782
a3 = 1.781477937
a4 = -1.821255978
a5 = 1.330274429

if x >= 0.0:
d = 1.0/(1.0+0.2316419*x)
N_x = 1.0 - 1.0/sqrt(2.0*pi)*exp(-x**2/2.0)*(a1*d + a2*d**2 + a3*d**3 + a4*d**4 + a5*d**5)
else:
N_x = 1.0 - cdf_of_normal(-x)

return(N_x)

def pdf_of_normal(x): #No sigma, just for Bachelier. Thomson 2016
return( 1.0/sqrt(2*pi)*exp(-0.5*x**2) )

def brent_dekker(a,b,f,epsilon,N_step,N_corr): #reference https://en.wikipedia.org/wiki/Brent%27s_method

if f(a)*f(b)<0:

if abs(f(a))<abs(f(b)):
return(brent_dekker(b,a,f,epsilon,N_step,N_corr))

else:
i=0
c=a
s=b
mflag=True
while(f(b)!=0.0 and f(s)!=0.0 and abs(b-a)>=epsilon and i<N_step):
i=i+1
if (f(a)!=f(c) and f(b)!=f(c)):
s=a*f(b)*f(c)/((f(a)-f(b))*(f(a)-f(c)))+b*f(a)*f(c)/((f(b)-f(a))*(f(b)-f(c)))+c*f(a)*f(b)/((f(c)-f(a))*(f(c)-f(b)))
else:
#secant method
s=b-f(b)*(b-a)/(f(b)-f(a))

if((s<(3.0*a+b)/4.0 or s>b) or (mflag==True and abs(s-b)>=abs(b-c)/2.0) or (mflag==False and abs(s-b)>=abs(c-d)/2.0) or (mflag==True and abs(b-c)<epsilon) or (mflag==False and abs(c-d)<epsilon)):
#bisection method
s=(a+b)/2.0
mflag=True
else:
mflag=False
d=c
c=b
if f(a)*f(s)<0.0:
b=s
else:
a=s
if abs(f(a))<abs(f(b)):
aux=a
a=b
b=aux

if i>=N_step: #it did not converge, it never happened
return(float('nan'))
else:
if(f(s)==0.0):
return s
else:
return b

else:
#I try to automate it knowing that volatility will be between 0 and infinity
#        return 'error'
if N_corr>0:

if a>b:
a=a*10.0
b=b/10.0
else:
a=a/10.0
b=b*10.0
#print a,b,N_corr,f(a),f(b) #for debug
N_corr = N_corr-1
return(brent_dekker(a,b,f,epsilon,N_step,N_corr))

else:
return(float('nan'))
#it happens. The problem is f(a) and f(b) remain constant.

#ABSTRACT CLASS FOR A PRICING MODEL
#I think these classes are useful if a person wants to play with implied volatility (iv) for a specific model with a diverse set of assets (call/put options on futures/stocks). In this way he can define a model, then update it and run the different methods.
#I don't know exactly how a normal day working in finance is so it might not be the smartest way to do it.

class pricing_model_iv(object):

def __init__(self, z):
self.z = z

def update(self):
raise NameError("Abstract Pricing Model")

def iv_call_stock(self):
raise NameError("Abstract Pricing Model")

def iv_put_stock(self):
raise NameError("Abstract Pricing Model")

def iv_call_future(self):
raise NameError("Abstract Pricing Model")

def iv_put_future(self):
raise NameError("Abstract Pricing Model")

#BLACK & SCHOLES PRICING MODEL

class bl_sch_iv(pricing_model_iv):

def __init__(self, V, S, E, time_to_exp, r, epsilon, nstep_conv):

self.V = V #option price
self.S = S #underlying price, either stock or future.
self.E = E #strike price
self.time_to_exp = time_to_exp #time to expiration
self.r = r #risk-free interest rate
self.epsilon = epsilon #precision
self.nstep_conv = nstep_conv #maximum number of steps to convergence

def update(self, V, S, E, time_to_exp, r, epsilon, nstep_conv):

self.V = V #option price
self.S = S #underlying price, either stock or future.
self.E = E #strike price
self.time_to_exp = time_to_exp #time to expiration
self.r = r #risk-free interest rate
self.epsilon = epsilon #precision
self.nstep_conv = nstep_conv #maximum number of steps to convergence

def iv_call_stock(self):
#Black & Scholes
def price_tozero(sigma):
d1=( log(self.S/self.E)+(self.r+sigma**2/2.0)*self.time_to_exp )/( sigma*sqrt(self.time_to_exp) )
d2=d1-sigma*sqrt(self.time_to_exp)
return(self.S*cdf_of_normal(d1)-self.E*exp(-self.r*self.time_to_exp)*cdf_of_normal(d2)-self.V)

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,10))

def iv_put_stock(self):
#Black & Scholes
def price_tozero(sigma):
d1=( log(self.S/self.E)+(self.r+sigma**2/2.0)*self.time_to_exp )/( sigma*sqrt(self.time_to_exp) )
d2=d1-sigma*sqrt(self.time_to_exp)
return(-self.S*cdf_of_normal(-d1)+self.E*exp(-self.r*self.time_to_exp)*cdf_of_normal(-d2)-self.V)

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

def iv_call_future(self):
#Black & Scholes + http://finance.bi.no/~bernt/gcc_prog/recipes/recipes/node9.html
def price_tozero(sigma):
d1=( log(self.S/self.E)+(self.r+sigma**2/2.0)*self.time_to_exp )/( sigma*sqrt(self.time_to_exp) )
d2=d1-sigma*sqrt(self.time_to_exp)
return((self.S*cdf_of_normal(d1)-self.E*cdf_of_normal(d2))*exp(-self.r*self.time_to_exp)-self.V)

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

def iv_put_future(self):
#Obtained using Put-Call Parity using the call future formula
#In the put-call parity relation I should discount also S because we are considering a future contract
def price_tozero(sigma):
d1=( log(self.S/self.E)+(self.r+sigma**2/2.0)*self.time_to_exp )/( sigma*sqrt(self.time_to_exp) )
d2=d1-sigma*sqrt(self.time_to_exp)
return((self.S*cdf_of_normal(d1)-self.E*cdf_of_normal(d2))*exp(-self.r*self.time_to_exp) + self.E*exp(-self.r*self.time_to_exp) - self.S*exp(-self.r*self.time_to_exp)  -self.V)

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

#BACHELIER PRICING MODEL

class bachl_iv(pricing_model_iv):

def __init__(self, V, S, E, time_to_exp, r, epsilon, nstep_conv):

self.V = V #option price
self.S = S #underlying price, either stock or future.
self.E = E #strike price
self.time_to_exp = time_to_exp #time to expiration
self.r = r #risk-free interest rate
self.epsilon = epsilon #precision
self.nstep_conv = nstep_conv #maximum number of steps to convergence

def update(self, V, S, E, time_to_exp, r, epsilon, nstep_conv):

self.V = V #option price
self.S = S #underlying price, either stock or future.
self.E = E #strike price
self.time_to_exp = time_to_exp #time to expiration
self.r = r #risk-free interest rate
self.epsilon = epsilon #precision
self.nstep_conv = nstep_conv #maximum number of steps to convergence

#Following 4.5.1 and 4.5.2 of Thomson 2016
#I converted the arithmetic compounding to continuous compounding, but not the normal to log-normal distribution of returns
#I thought the distribution of the returns was an important part of the Bachelier model, while continuous compounding is related to the interest rate that is given in the data

def iv_call_stock(self):
#From Thomson 2016, Eq. 31c, paragraph 4.2.3
def price_tozero(sigma):
d=(self.S - self.E*exp(-self.r*self.time_to_exp))/(self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp))
return( (self.S - self.E*exp(-self.r*self.time_to_exp))*cdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp)*pdf_of_normal(d) - self.V )

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

def iv_put_stock(self):
#From IV_CALL_STOCK method plus put-call parity (Paul Wilmott on Quantitative Finance, it is model independent.) => P = C + E*e^(-r*time_to_exp) - S
def price_tozero(sigma):
d=(self.S - self.E*exp(-self.r*self.time_to_exp))/(self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp))
return( (self.S - self.E*exp(-self.r*self.time_to_exp))*cdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp)*pdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp) - self.S - self.V )

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

#For the options on futures, I tried to replicate the scheme I used in Black and Scholes of discounting S together with E.
#This means that the first term of the price pass from:
#    (self.S - self.E*exp(-self.r*self.time_to_exp))*cdf_of_normal(d)
#to:
#    (self.S - self.E)*exp(-self.r*self.time_to_exp)*cdf_of_normal(d)

def iv_call_future(self):
def price_tozero(sigma):
d=(self.S - self.E*exp(-self.r*self.time_to_exp))/(self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp))
return( (self.S - self.E)*exp(-self.r*self.time_to_exp)*cdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp)*pdf_of_normal(d) - self.V )

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

#In the put-call parity relation I should discount also S because we are considering a future contract

def iv_put_future(self):
def price_tozero(sigma):
d=(self.S - self.E*exp(-self.r*self.time_to_exp))/(self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp))
return( (self.S - self.E)*exp(-self.r*self.time_to_exp)*cdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp)*sigma*sqrt(self.time_to_exp)*pdf_of_normal(d) + self.E*exp(-self.r*self.time_to_exp) - self.S*exp(-self.r*self.time_to_exp) - self.V )

return(brent_dekker(1.0,99.0,price_tozero,self.epsilon,self.nstep_conv,3))

#NOW I START THE CODE

blsc = bl_sch_iv(1,1,1,1,1,1,1)
bach = bachl_iv(1,1,1,1,1,1,1)

precision = 1.e-8
max_n = 10**6 #max number of steps to convergence for the Brent Dekker zero finding algorithm

ids_u_F=set()
ids_u_S=set()
ids_o_C=set()
ids_o_P=set()
ids_m_Ba=set()
ids_m_BS=set()

sick=set()

with open('input.csv','rb') as csvfile, open('output.csv','wb') as csvout:

csvfile.seek(0) #rewind
writing = csv.writer(csvout, delimiter=',')
writing.writerow(['ID','Spot','Strike','Risk-Free Rate','Years to Expiry','Option Type','Model Type','Implied Volatility','Market Price'])
# I did this just in case the next time you give me a CSV without the header
#print ', '.join(row)
ID=row[0]
underlying_type=row[1]
underlying=float(row[2])
risk_free_rate=float(row[3])*(-1.) #all the interest rates are negative values (in the CSV file). I think the - sign is given for granted. It's usually a yearly interest rate
days_to_expiry=float(row[4])/365.0 #Converting from days to expiry into years to expiry
strike=float(row[5])
option_type=row[6]
model_type=row[7]
market_price=float(row[8])

if model_type=='BlackScholes':
else:
if underlying_type=='Stock':
else:
if option_type=='Call':
else:

if model_type=='BlackScholes':
blsc.update(market_price, underlying, strike, days_to_expiry, risk_free_rate, precision, max_n)
if underlying_type=='Stock':
if option_type=='Call':
iv=blsc.iv_call_stock()
else:
iv=blsc.iv_put_stock()
else:
if option_type=='Call':
iv=blsc.iv_call_future()
else:
iv=blsc.iv_put_future()
else:
bach.update(market_price, underlying, strike, days_to_expiry, risk_free_rate, precision, max_n)
if underlying_type=='Stock':
if option_type=='Call':
iv=bach.iv_call_stock()
else:
iv=bach.iv_put_stock()
else:
if option_type=='Call':
iv=bach.iv_call_future()
else:
iv=bach.iv_put_future()

#Writing the csv file
if underlying_type=='Stock':
writing.writerow([row[0],row[2],row[5],row[3],str(days_to_expiry),row[6],row[7],str(iv),row[8]])
else:#Spot price for futures
writing.writerow([row[0],str(underlying*exp(-risk_free_rate*days_to_expiry)),row[5],row[3],str(days_to_expiry),row[6],row[7],str(iv),row[8]])

#Count of nans
if isnan(iv):
if model_type=='BlackScholes':
else:
if underlying_type=='Stock':
else:
if option_type=='Call':
else:

#It returns how many options have NaN volatility. It also allows to study them

print len(bad_ids), ' out of 65535 that is the ', len(bad_ids)/65535.0*100.0, '% \n'


• not a full answer but my immediate thought was that the mass of "magic numbers" makes the code impenetrable. Commented Jun 25, 2018 at 8:46
• Hi, I'm a mod at Quant.Se and migrated from there, I think this is a first. From what I know of this Stack this would be on-topic here but please correct me if I'm wrong. Commented Jun 25, 2018 at 9:00
• @BobJansen Seems like an on-topic question here. A bit more context would be welcome though, for those of us not used to the technical jargon of quantitative finance. Commented Jun 25, 2018 at 9:38
• @BobJansen Looks fine here. If you have any questions left after reading the guide, feel free to find us in chat. There's usually plenty of people around to answer questions should anything be unclear.
– Mast
Commented Jun 25, 2018 at 10:05
• @spec3 Can you edit your question to include such information (possibly with links to the mentionned videos) into the question? Commented Jun 25, 2018 at 13:26

• Give the operators some breathing space. For example, the

    s=a*f(b)*f(c)/((f(a)-f(b))*(f(a)-f(c)))+b*f(a)*f(c)/((f(b)-f(a))*(f(b)-f(c)))+c*f(a)*f(b)/((f(c)-f(a))*(f(c)-f(b)))


line is a textbook example of unreadable. Consider at least

    s = a*f(b)*f(c)/((f(a)-f(b))*(f(a)-f(c))) +
b*f(a)*f(c)/((f(b)-f(a))*(f(b)-f(c))) +
c*f(a)*f(b)/((f(c)-f(a))*(f(c)-f(b)))

• Another example of unreadability is the

        if((s<(3.0*a+b)/4.0 or s>b) or (mflag==True and abs(s-b)>=abs(b-c)/2.0) or (mflag==False and abs(s-b)>=abs(c-d)/2.0) or (mflag==True and abs(b-c)<epsilon) or (mflag==False and abs(c-d)<epsilon)):


line which is 211 character long, and seriously hides the logics. Consider

        if (s < (3.0*a+b)/4.0 or s > b) or
(mflag==True and abs(s-b)>=abs(b-c)/2.0) or
(mflag==False and abs(s-b)>=abs(c-d)/2.0) or
(mflag==True and abs(b-c)<epsilon) or
(mflag==False and abs(c-d)<epsilon):


As a sidenote, you seem to assume that b is greater than (3.0*a+b)/4.0 (the reference code only says that s is not between (3.0*a+b)/4.0 and b). I don't see why it is necessarily the case. If it is so, the pythonic way to express it is

            not (3.0*a+b)/4.0 < s < b

• Comparing floating point values for equality is practically guaranteed to fail. Comparing them for inequality doesn't help much either.

• Comments like #inverse quadratic interpolation are clear indication that the commented code needs to be factored out into a function. Consider

        if (f(a)!=f(c) and f(b)!=f(c)):
s = inverse_quadratic_interpolation(f, a, b, c)
else:
s = secant_interpolation(f, a, b)

• A pythonic way to swap two values is

        a, b = b, a

• Given the f(a) * f(b) < 0 bracketing constraint, the Brent-Dekker method (just like any other decent solver) guarantees convergence. I don't understand the rationale of passing max number of iterations.

As a sidenote, selecting a Brent-Dekker solver is rather arbitrary, and has nothing to do with the problem domain. Leave the selection of solver open.

• It is not a solver's responsibility to deal with a violation of the bracketing constraint. It is the caller's job.

• All the price_to_zero functions are too similar to warrant 4 copies in bl_sch_iv, and 4 more in bachl_iv. Some parameter tweaking is all you need.

• Statements that span multiple lines need a backslash before the newline to indicate that a continuation line follows. Alternatively, an unclosed parenthesis will cause a statement to not be terminated. Commented Jun 26, 2018 at 4:48