6
\$\begingroup\$

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)  
\$\endgroup\$
  • \$\begingroup\$ not a full answer but my immediate thought was that the mass of "magic numbers" makes the code impenetrable. \$\endgroup\$ – MD-Tech Jun 25 '18 at 8:46
  • \$\begingroup\$ 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. \$\endgroup\$ – Bob Jansen Jun 25 '18 at 9:00
  • 4
    \$\begingroup\$ @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. \$\endgroup\$ – 409_Conflict Jun 25 '18 at 9:38
  • \$\begingroup\$ @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. \$\endgroup\$ – Mast Jun 25 '18 at 10:05
  • 1
    \$\begingroup\$ @spec3 Can you edit your question to include such information (possibly with links to the mentionned videos) into the question? \$\endgroup\$ – 409_Conflict Jun 25 '18 at 13:26
2
\$\begingroup\$
  • 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.

\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – 200_success Jun 26 '18 at 4:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.