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)):
#inverse quadratic interpolation
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
bad_ids=[]
bad_ids_u_F=set()
bad_ids_u_S=set()
bad_ids_o_C=set()
bad_ids_o_P=set()
bad_ids_m_Ba=set()
bad_ids_m_BS=set()
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:
has_header = csv.Sniffer().has_header(csvfile.read(10)) #I check the header existence with a little sample
csvfile.seek(0) #rewind
reading = csv.reader(csvfile, delimiter=',')
writing = csv.writer(csvout, delimiter=',')
writing.writerow(['ID','Spot','Strike','Risk-Free Rate','Years to Expiry','Option Type','Model Type','Implied Volatility','Market Price'])
if has_header:
next(reading) #skip header
# I did this just in case the next time you give me a CSV without the header
for row in reading:
#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':
ids_m_BS.add(ID)
else:
ids_m_Ba.add(ID)
if underlying_type=='Stock':
ids_u_S.add(ID)
else:
ids_u_F.add(ID)
if option_type=='Call':
ids_o_C.add(ID)
else:
ids_o_P.add(ID)
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):
bad_ids.append(ID)
if model_type=='BlackScholes':
bad_ids_m_BS.add(ID)
else:
bad_ids_m_Ba.add(ID)
if underlying_type=='Stock':
bad_ids_u_S.add(ID)
else:
bad_ids_u_F.add(ID)
if option_type=='Call':
bad_ids_o_C.add(ID)
else:
bad_ids_o_P.add(ID)
#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'
print 'BS-CALL-STOCK: ', len(bad_ids_m_BS & bad_ids_u_S & bad_ids_o_C)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_BS & bad_ids_u_S & bad_ids_o_C)
print 'BS-PUT-STOCK: ', len(bad_ids_m_BS & bad_ids_u_S & bad_ids_o_P)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_BS & bad_ids_u_S & bad_ids_o_P)
print 'BS-CALL-FUTURE: ', len(bad_ids_m_BS & bad_ids_u_F & bad_ids_o_C)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_BS & bad_ids_u_F & bad_ids_o_C)
print 'BS-PUT-FUTURE: ', len(bad_ids_m_BS & bad_ids_u_F & bad_ids_o_P)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_BS & bad_ids_u_F & bad_ids_o_P)
print 'Ba-CALL-STOCK: ', len(bad_ids_m_Ba & bad_ids_u_S & bad_ids_o_C)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_Ba & bad_ids_u_S & bad_ids_o_C)
print 'Ba-PUT-STOCK: ', len(bad_ids_m_Ba & bad_ids_u_S & bad_ids_o_P)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_Ba & bad_ids_u_S & bad_ids_o_P)
print 'Ba-CALL-FUTURE: ', len(bad_ids_m_Ba & bad_ids_u_F & bad_ids_o_C)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_Ba & bad_ids_u_F & bad_ids_o_C)
print 'Ba-PUT-FUTURE: ', len(bad_ids_m_Ba & bad_ids_u_F & bad_ids_o_P)*100.0/len(bad_ids), '% real: ', len(bad_ids_m_Ba & bad_ids_u_F & bad_ids_o_P)
