trjhtr

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) 

• 
  

  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.