Black-Scholes Model in Python returning negative Call price-CodePudding

The Black-Scholes Code returns negative call price which is wrong. Can you please help me figure out the mistake.

    ## define two functions, d1 and d2 in Black-Scholes model
def d1(S,K1,T,r,sigma):
    return(log(S/K1) (r (sigma**2)/2.0)* (T) )/sigma*sqrt(T)
def d2(S,K1,T,r,sigma):
    return d1(S,K1,T,r,sigma)-sigma*sqrt(T)

## define the call options price function
def bs_call(S,K1,T,r,sigma):
  S = float(S)
  K1 = float(K1)
  T = float(T)
  r = float(r)
  sigma = float(sigma)
  return S*norm.cdf(d1(S,K1,T,r,sigma))-K1*exp(-r*T)*norm.cdf(d1(S,K1,T,r,sigma)-sigma*np.sqrt(T))

bsc = bs_call(125.5,187.5,0.70136,.005896996,.24263)
print(bsc)

The output of the print statement is:

-0.24854522041832539

Thanks.

CodePudding user response：

Do you know what you expect to be returned for those input parameters?

Anyway, one issue I can see is in your definition of d1:

def d1(S,K1,T,r,sigma):
    return(log(S/K1) (r (sigma**2)/2.0)* (T) )/sigma*sqrt(T)

the sigma*sqrt(T) term is being resolved as sqrt(T) * (log...)/sigma

Put some brackets around this term, or rewrite the function as something like:

def d1(S, K1, T, r, sigma):
     a = np.log(S / K1)
     b = T * (r   sigma ** 2 / 2)
     c = sigma * np.sqrt(T)
     return (a   b) / c

Using this returns a positive number (0.296164 to 6 s.f.). Generally speaking, you can probably help yourself out by not trying to one line everything and breaking it into logical chunks.

CodePudding user response：

Everything looks good except that you are missing a pair of parentheses in the d1() function.

Please see notes in comments in the code below:

        ## define two functions, d1 and d2 in Black-Scholes model
        def d1(S,K1,T,r,sigma):
            '''
            # Your original calculation was missing parentheses around sigma*sqrt(T)
            return(log(S/K1) (r (sigma**2)/2.0)* (T) )/sigma*sqrt(T)
            '''
            #This will give the correct result
            return(log(S/K1) (r (sigma**2)/2.0)* (T) )/  (sigma*sqrt(T))
        def d2(S,K1,T,r,sigma):
            return d1(S,K1,T,r,sigma)-sigma*sqrt(T)

        ## define the call options price function
        def bs_call(S,K1,T,r,sigma):
            S = float(S)
            K1 = float(K1)
            T = float(T)
            r = float(r)
            sigma = float(sigma)
            
            from scipy.stats import norm
            import numpy as np
            '''
            # For the value of a call option for a non-dividend-paying underlying stock, Wikipedia says:
            #   P(S, 0) = N(d1)*S - N(d2)*K*exp(-r*T)
            # where:
            #   d1 = (log(S/K)   (r   T/2 * sigma**2)) / (sigma * sqrt(T))
            #   d2 = D1 - sigma * sqrt(T)
            #
            # You have written the price as:
            #   S*N(d1) - K*exp(-r*T)*N(d1 - sigma*sqrt(T))
            # which is equivalent to:
            #   N(d1)*S - N(d2)*K*exp(-r*T)
            #
            # This looks fine. Note that you should be able to simplify the logic a bit as shown below.

            return S*norm.cdf(d1(S,K1,T,r,sigma))-K1*exp(-r*T)*norm.cdf(d1(S,K1,T,r,sigma)-sigma*np.sqrt(T))
            '''
            original_calculation = S*norm.cdf(d1(S,K1,T,r,sigma))-K1*exp(-r*T)*norm.cdf(d1(S,K1,T,r,sigma)-sigma*np.sqrt(T))
            
            alternative1 = S*norm.cdf(d1(S,K1,T,r,sigma))-K1*exp(-r*T)*norm.cdf(d2(S,K1,T,r,sigma))
            
            d1_value = d1(S,K1,T,r,sigma)
            d2_value = d1_value - sigma*np.sqrt(T)
            alternative2 = S*norm.cdf(d1_value)-K1*exp(-r*T)*norm.cdf(d2_value)
            
            print(original_calculation)
            print(alternative1)
            print(alternative2)
            
            return alternative2

        bsc = bs_call(125.5,187.5,0.70136,.005896996,.24263)
        print(bsc)

The above code prints:

0.29616408169941755
0.29616408169941755
0.29616408169941755
0.29616408169941755