PHILADELPHIA REFLECTIONS
Musings of a Philadelphia Physician who has served the community for six decades

Return to Home

Related Topics

Using XML to Retrieve Financial Data
It's often desirable to get live financial data and everyone knows. XML is the thing to use but actually writing programs that work takes a bit of trouble. Plus, once you've got the data you need to display it.

Black-Scholes-Merton in Python

The provenance is numerous sources including Back, Benninga, Hull and Wilmott but I wrote the majority of this myself while I was at MIT in the Master of Finance program.

# Black Scholes Merton for Python (from Excel VBA)
# George Fisher MIT Spring 2012
#
# I created BSM routines for Excel first and then converted to Python
# I drew upon published work by Back, Benninga, Hull and Wilmott but the majority is my own original work
#
#      d1, d2
#      N  (en/phi) std normal CDF
#      N' (nprime) std normal PDF
#
#      Binary Options
#      Euro Call & Put
#      Greeks
#      Implied Volatility

#      Put/Call Parity

#      American_Call_Dividend
#      American_Put_Binomial

#
# Also includes
#      risk-neutral prob
#      forward prices & rates
#      CAGR
#      randn/randn_ssdt
#      convert discrete to continuous interest rate

# Interest is
#       risk-free rate
#       domestic risk-free rate for currencies

# Yield is
#       dividend yield for stock
#       lease rate for commodities
#       foreign currency risk-free rate for currencies

# Sigma is the standard deviation of the underlying asset

# Time is a year fraction: for 3-months ... Time = 3/12

# Stock is S_0

# Exercise is K

# => Interest, Yield, Sigma, Time are all annual numbers
# => Time = 0 is the value at maturity
#        most of the functions accomodate this
#        for some, it's infinity or otherwise meaningless
# => Sigma = 0 is also accomodated in most functions

##
##   Utilities
##   ---------
##

# N: the standard-normal CDF

def en(x):
    return phi(x)
    

def phi(x):
    import math
    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # Save the sign of x
    sign = 1
    if x < 0:
        sign = -1
    x = abs(x)/math.sqrt(2.0)

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)

    return 0.5*(1.0 + sign*y)


# N': the first derivative of N(x) ... the standard normal PDF

def nprime(x):
    import math
    return math.exp(-0.5 * x * x) / math.sqrt(2 * 3.1415926)

# Random Normal (epsilon)

def RandN():
    # produces a standard normal random variable epsilon
    import random
    return random.gauss(0,1)

def RandN_ssdt(ssdt):
    # produces a standard normal random variable epsilon times sigma*sqrt(deltaT)
    import random
    return random.gauss(0,ssdt)

# binomial tree risk-neutral probability (Hull 7th edition Ch 19 P 409)

def RiskNeutralProb(Interest, Yield, sigma, deltaT):
    import math

    u = math.exp( sigma * math.sqrt(deltaT))
    d = math.exp(-sigma * math.sqrt(deltaT))
    
    a = math.exp((Interest - Yield) * deltaT)
    
    numerator = a - d
    denominator = u - d
    
    return numerator / denominator

# Call & Put prices derived from put-call parity
#                                ---------------

def CallParity(Stock, Exercise, Time, Interest, Yield, Put_price):

    import math
    return Put_price + Stock * math.exp(-Yield * Time) - Exercise * math.exp(-Interest * Time)

def PutParity(Stock, Exercise, Time, Interest, Yield, Call_price):

    import math
    return Call_price + Exercise * math.exp(-Interest * Time) - Stock * math.exp(-Yield * Time)

# forward price

def ForwardPrice(Spot, Time, Interest, Yield):

    import math
    return Spot * math.exp((Interest - Yield) * Time)

# forward rate from Time1 to Time2 (discrete compounding)

def ForwardRate(SpotInterest1, Time1, SpotInterest2, Time2):

    numerator   = (1 + SpotInterest2) ** Time2
    denominator = (1 + SpotInterest1) ** Time1
    
    return ((numerator / denominator) ** (1 / (Time2 - Time1))) - 1

# CAGR

def CAGRd(Starting_value, Ending_Value, Number_of_years):

    # discrete CAGR

    return ((Ending_Value / Starting_value) ** (1 / Number_of_years)) - 1

# Convert TO continuous compounding FROM discrete

def r_continuous(r_discrete, compounding_periods_per_year):

    import math
    m = compounding_periods_per_year
    return m * math.log(1 + r_discrete / m)

# Convert TO discrete compounding FROM continuous
#
# t_discrete = m * (exp(r_continuous / m) - 1)
#
# where m is the number of compounding periods per year
#
def r_discrete(r_continuous, compounding_periods_per_year):

  import math
  m = compounding_periods_per_year
  return m * (math.exp(r_continuous / m) - 1)

# --------------------------------------------------------------------------------

##
##   Black Scholes
##   -------------
##

def dOne(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    return (math.log(Stock / Exercise) + (Interest - Yield + 0.5 * sigma * sigma) * Time) / (sigma * math.sqrt(Time))

def dTwo(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    return (math.log(Stock / Exercise) + (Interest - Yield - 0.5 * sigma * sigma) * Time) / (sigma * math.sqrt(Time))

#
# Binary Options
#

# Digital: Cash or Nothing

def CashCall(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time < 0.000000005:
        if Stock >= Exercise:
            return 1
        else:
            return 0
    
    if sigma == 0:
        sigma = 0.0000000001
    
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    Nd2 = phi(d2_)
    
    return math.exp(-Interest * Time) * Nd2

def CashPut(Stock, Exercise, Time, Interest, Yield, sigma):

    if Time < 0.000000005:
        if Stock >= Exercise:
            return 0
        else:
            return 1
    
    if sigma == 0:
        sigma = 0.0000000001
    
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    Nminusd2 = phi(-d2_)
    
    return math.exp(-Interest * Time) * Nminusd2

# Asset or Nothing

def AssetCall(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time < 0.000000005:
        if Stock >= Exercise:
            return Stock
        else:
            return 0   
    
    if sigma == 0:
        sigma = 0.0000000001
    
    if Exercise < 0.000000005:
        Nd1 = 1
    else:
        d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
        Nd1 = phi(d1_)
    
    return Stock * math.exp(-Yield * Time) * Nd1

def AssetPut(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time < 0.000000005:
        if Stock >= Exercise:
            return 0
        else:
            return Stock
    
    if sigma == 0:
        sigma = 0.0000000001
    
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    Nminusd1 = phi(-d1_)
    
    return Stock * math.exp(-Yield * Time) * Nminusd1

#
# European Call and Put
# ---------------------
#

def EuroCall(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time == 0:
        return max(0, Stock - Exercise)
    
    if sigma == 0:
        return max(0, math.exp(-Yield * Time) * Stock - math.exp(-Interest * Time) * Exercise)
    
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)

    return Stock * math.exp(-Yield * Time) * phi(d1_) - Exercise * math.exp(-Interest * Time) * phi(d2_)

def EuroPut(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time == 0:
        return max(0, Exercise - Stock)
    
    if sigma == 0:
        return max(0, math.exp(-Interest * Time) * Exercise - math.exp(-Yield * Time) * Stock)
        
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)

    return Exercise * math.exp(-Interest * Time) * phi(-d2_) - Stock * math.exp(-Yield * Time) * phi(-d1_)


#
# American Put
# ------------
# Per Kerry Back Chapt5.bas
#

def American_Put_Binomial(S0, K, r, sigma, q, T, N):
    import math
    """
    #
    # Inputs are S0 = initial stock price
    #            K = strike price
    #            r = risk-free rate
    #            sigma = volatility
    #            q = dividend yield
    #            T = time to maturity
    #            N = number of time periods
    #
    """
    dt = T / N                                   # length of time period
    u = math.exp(sigma * math.sqrt(dt))          # size of up step
    d = 1 / u                                    # size of down step
    pu = (math.exp((r - q) * dt) - d) / (u - d)  # probability of up step
    dpu = math.exp(-r * dt) * pu                 # one-period discount x prob of up step
    dpd = math.exp(-r * dt) * (1 - pu)           # one-period discount x prob of down step
    u2 = u * u
    S = S0 * d ** N                              # stock price at bottom node at last date
    PutV[0] = max(K - S, 0)                      # put value at bottom node at last date
    for j in range(1,N+1):
        S = S * u2
        PutV[j] = max(K - S, 0)

    for i in range(N - 1, 0, -1):                # back up in time to date 0
        S = S0 * d ** i                          # stock price at bottom node at date i
        PutV[0] = max(K - S, dpd * PutV(0) + dpu * PutV(1))
        for j in range(1,i+1):                   # step up over nodes at date i
            S = S * u2
            PutV[j] = max(K - S, dpd * PutV(j) + dpu * PutV(j + 1))

    return PutV[0]                               # put value at bottom node at date 0

#
# Greeks from Hull (Edition 7) Chapter 17 p378
# --------------------------------------------
#

# per the Black Scholes PDE for a portfolio of options
# on a single non-dividend-paying underlying stock
#
# Theta + Delta * S * r + Gamma * 0.5 * sigma**2 * S**2 = r * Portfolio_Value

# Per Hull: for large option portfolios, usually created by banks in the
# course of buying and selling OTC options to clients, the portfolio is
# Delta hedged every day and Gamma/Vega hedged as needed
#
#             Delta      Gamma      Vega
# Portfolio   P_delta    P_gamma    P_vega
# Option1     w1*1_delta w1*1_gamma w1*1_vega
# Option2     w2*2_delta w2*2_gamma w2*2_vega
#
# Set the columns equal to zero and solve the simultaneous equations

# Most OTC options are sold close to the money; high gamma and vega
# as (if) the price of the underlying move away gamma and vega decline

# Delta
# -----
#
# If a bank sells a call to a client (short a call)
#   ... it hedges itself with a synthetic long call
#
# Synthetic long call = Delta * Stock_price - bond
# ie., borrow the money to buy Delta shares of the stock
#

def DeltaCall(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time == 0:
        if Stock > Exercise:
            return 1
        else:
            return 0
    
    if sigma == 0:
        sigma = 0.0000000001

    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)

    return math.exp(-Yield * Time) * phi(d1_)

def DeltaPut(Stock, Exercise, Time, Interest, Yield, sigma):

    import math
    if Time == 0:
        if Stock < Exercise:
            return -1
        else:
            return 0
    
    if sigma == 0:
        sigma = 0.0000000001
    
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)

    return math.exp(-Yield * Time) * (phi(d1_) - 1)


#
# Gamma the convexity
# -----
#

def OptionGamma(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)

    OptionGamma = nprime(d1_) * math.exp(-Yield * Time) _
        / (Stock * sigma * math.sqrt(Time))

#
# Theta the decay in the value of an option/portfolio of options as time passes
# -----
#
# divide by 365 for "per calendar day"; 252 for "per trading day"
#
# In a delta-neutral portfolio, Theta is a proxy for Gamma
#

def ThetaCall(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    Dim d2_
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nd1_
    Nd1_ = phi(d1_)
    Dim Nd2_
    Nd2_ = phi(d2_)
    
    ThetaCall = -Stock * nprime(d1_) * sigma * math.exp(-Yield * Time) / (2 * math.sqrt(Time)) _
        + Yield * Stock * Nd1_ * math.exp(-Yield * Time) _
        - Interest * Exercise * math.exp(-Interest * Time) * Nd2_

def ThetaPut(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    Dim d2_
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nminusd1_
    Nminusd1_ = phi(-d1_)
    Dim Nminusd2_
    Nminusd2_ = phi(-d2_)
    
    ThetaPut = -Stock * nprime(d1_) * sigma * math.exp(-Yield * Time) / (2 * math.sqrt(Time)) _
        - Yield * Stock * Nminusd1_ * math.exp(-Yield * Time) _
        + Interest * Exercise * math.exp(-Interest * Time) * Nminusd2_

#
# Vega the sensitivity to changes in the volatility of the underlying
# ----
#
def Vega(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Vega = Stock * math.sqrt(Time) * nprime(d1_) * math.exp(-Yield * Time)

#
# Rho the sensitivity to changes in the interest rate
# ---
#

#
# Note the various Rho calculations see Hull 7th Edition Ch 17 P378
#

def RhoFuturesCall(Stock, Exercise, Time, Interest, Yield, sigma):

    RhoFuturesCall = -EuroCall(Stock, Exercise, Time, Interest, Yield, sigma) * Time

def RhoFuturesPut(Stock, Exercise, Time, Interest, Yield, sigma):

    RhoFuturesPut = -EuroPut(Stock, Exercise, Time, Interest, Yield, sigma) * Time

#
# The Rho corresponding to the domestic interest rate is RhoCall/Put, below
#                              foreign  interest rate is RhoFXCall/Put, shown here
#
def RhoFXCall(Stock, Exercise, Time, Interest, Yield, sigma):

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nd1_
    Nd1_ = phi(d1_)
    
    RhoFXCall = -Time * math.exp(-Yield * Time) * Stock * Nd1_

def RhoFXPut(Stock, Exercise, Time, Interest, Yield, sigma):

    Dim d1_
    d1_ = dOne(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nminusd1_
    Nminusd1_ = phi(-d1_)
    
    RhoFXPut = Time * math.exp(-Yield * Time) * Stock * Nminusd1_

#
# "Standard" Rhos
#

def RhoCall(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d2_
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nd2_
    Nd2_ = phi(d2_)
    
    RhoCall = Exercise * Time * math.exp(-Interest * Time) * Nd2_

def RhoPut(Stock, Exercise, Time, Interest, Yield, sigma):

    If sigma = 0 Then
        sigma = 0.0000000001
    End If

    Dim d2_
    d2_ = dTwo(Stock, Exercise, Time, Interest, Yield, sigma)
    
    Dim Nminusd2_
    Nminusd2_ = phi(-d2_)
    
    RhoPut = -Exercise * Time * math.exp(-Interest * Time) * Nminusd2_

#
# Since Bennigna and Back produce identical numbers
# and MATLAB produced numbers that are +/- 2%, I'm
# inclined to go with these numbers
#

#
# Implied Volatility from Benningna
# ---------------------------------
#
def EuroCallVol(Stock, Exercise, Time, Interest, Yield, Call_price):

    Dim High, Low
    High = 2
    Low = 0
    Do While (High - Low) > 0.000001
    If EuroCall(Stock, Exercise, Time, Interest, Yield, (High + Low) / 2) > _
        Call_price Then
             High = (High + Low) / 2
             Else: Low = (High + Low) / 2
    End If
    Loop
    EuroCallVol = (High + Low) / 2

def EuroPutVol(Stock, Exercise, Time, Interest, Yield, Put_price):

    Dim High, Low
    High = 2
    Low = 0
    Do While (High - Low) > 0.000001
    If EuroPut(Stock, Exercise, Time, Interest, Yield, (High + Low) / 2) > _
        Put_price Then
             High = (High + Low) / 2
             Else: Low = (High + Low) / 2
    End If
    Loop
    EuroPutVol = (High + Low) / 2

#
# Implied Volatility from Kerry Back p64
# Chapt3.bas Newton Raphson technique
# Answer IDENTICAL to Bennigna (EuroCallVol)
#
def Black_Scholes_Call(S, K, r, sigma, q, T):

    Black_Scholes_Call = EuroCall(S, K, T, r, q, sigma)

def Black_Scholes_Call_Implied_Vol(S, K, r, q, T, CallPrice):
#
# Inputs are S = initial stock price
#            K = strike price
#            r = risk-free rate
#            q = dividend yield
#            T = time to maturity
#            CallPrice = call price
#
Dim tol, lower, flower, upper, fupper, guess, fguess
If CallPrice < math.exp(-q * T) * S - math.exp(-r * T) * K Then
    MsgBox ("Option price violates the arbitrage bound.")
    Exit Function
End If
tol = 10 ^ -6
lower = 0
flower = Black_Scholes_Call(S, K, r, lower, q, T) - CallPrice
upper = 1
fupper = Black_Scholes_Call(S, K, r, upper, q, T) - CallPrice
Do While fupper < 0                   # double upper until it is an upper bound
    upper = 2 * upper
    fupper = Black_Scholes_Call(S, K, r, upper, q, T) - CallPrice
Loop
guess = 0.5 * lower + 0.5 * upper
fguess = Black_Scholes_Call(S, K, r, guess, q, T) - CallPrice
Do While upper - lower > tol               # until root is bracketed within tol
    If fupper * fguess < 0 Then            # root is between guess and upper
        lower = guess                      # make guess the new lower bound
        flower = fguess
        guess = 0.5 * lower + 0.5 * upper  # new guess = bi-section
        fguess = Black_Scholes_Call(S, K, r, guess, q, T) - CallPrice
    Else                                   # root is between lower and guess
        upper = guess                      # make guess the new upper bound
        fupper = fguess
        guess = 0.5 * lower + 0.5 * upper  # new guess = bi-section
        fguess = Black_Scholes_Call(S, K, r, guess, q, T) - CallPrice
    End If
Loop
Black_Scholes_Call_Implied_Vol = guess

#
# Implied Volatility from Wilmott Into Ch 8 p192 Newton Raphson***NOT DEBUGGED***
#
def ImpVolCall(Stock, Exercise, Time, Interest, Yield, Call_price):

    Volatility = 0.2
    epsilon = 0.0001
    dv = epsilon + 1
    
    While Abs(dv) > epsilon
        PriceError = EuroCall(Stock, Exercise, Time, Interest, Yield, Volatility) - Call_price
        dv = PriceError / Vega(Stock, Exercise, Time, Interest, Yield, Volatility)
        Volatility = Volatility - dv
    Wend
    
    ImpVolCall = Volatility



#
# from Kerry Back Chapt8.bas ... need Python's "BiNormalProb"
#
def American_Call_Dividend(S, K, r, sigma, Div, TDiv, TCall):
    import math
    """
#
# Inputs are S = initial stock price
#            K = strike price
#            r = risk-free rate
#            sigma = volatility
#            Div = cash dividend
#            TDiv = time until dividend payment
#            TCall = time until option matures >= TDiv
#
"""
    LessDiv = S - math.exp(-r * TDiv) * Div               # stock value excluding dividend
    If Div / K <= 1 - math.exp(-r * (TCall - TDiv)):      # early exercise cannot be optimal
        return Black_Scholes_Call(LessDiv, K, r, sigma, 0, TCall)
    #
    # Now we find an upper bound for the bisection.
    #
    upper = K
    while upper + Div - K < Black_Scholes_Call(upper, K, r, sigma, 0, TCall - TDiv):
        upper = 2 * upper
    #
    # Now we use bisection to compute Zstar = LessDivStar.
    #
    tol = 10 ** -6
    lower = 0
    flower = Div - K
    fupper = upper + Div - K - Black_Scholes_Call(upper, K, r, sigma, 0, TCall - TDiv)
    guess = 0.5 * lower + 0.5 * upper
    fguess = guess + Div - K - Black_Scholes_Call(guess, K, r, sigma, 0, TCall - TDiv)
    
    while upper - lower > tol:
        if fupper * fguess < 0:
            lower = guess
            flower = fguess
            guess = 0.5 * lower + 0.5 * upper
            fguess = guess + Div - K - Black_Scholes_Call(guess, K, r, sigma, 0, TCall - TDiv)
        else:
            upper = guess
            fupper = fguess
            guess = 0.5 * lower + 0.5 * upper
            fguess = guess + Div - K - Black_Scholes_Call(guess, K, r, sigma, 0, TCall - TDiv)

    LessDivStar = guess
    #
    # Now we calculate the probabilities and the option value.
    #
    d1 = (math.log(LessDiv / LessDivStar) + (r + sigma ** 2 / 2) * TDiv) / (sigma * math.sqrt(TDiv))
    d2 = d1 - sigma * math.sqrt(TDiv)
    d1prime = (math.log(LessDiv / K) + (r + sigma ** 2 / 2) * TCall) / (sigma * math.sqrt(TCall))
    d2prime = d1prime - sigma * math.sqrt(TCall)
    rho = -math.sqrt(TDiv / TCall)
    N1 = phi(d1)
    N2 = phi(d2)
    M1 = BiNormalProb(-d1, d1prime, rho)
    M2 = BiNormalProb(-d2, d2prime, rho)
    return LessDiv * N1 + math.exp(-r * TDiv) * (Div - K) * N2 + LessDiv * M1 - math.exp(-r * TCall) * K * M2


My thanks to Encode / Decode HTML Entities

(2396)

Please Let Us Know What You Think


(HTML tags provide better formatting)

Because of robot spam we ask you to confirm your comment: we will send you an email containing a link to click. We apologize for this inconvenience but this ensures the quality of the comments. (Your email will not be displayed.)
Thank you.