Related Topics
George (3)
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.
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
Originally published: Saturday, December 15, 2012; most-recently modified: Tuesday, May 14, 2019