Black-Scholes

Please see the attached code for an implementation of a standard Black-Scholes model for options pricing and risk management. Includes functions for valuation of first, second, and third order Greeks.

Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
import numpy as np
from math import log, sqrt, erf, exp
import scipy.stats

# Library for valuation of European call and put options using the
# standard Black-Scholes model

# S is the current stock price
# K is the strike price
# r is the risk-free rate
# sigma is the annualized volatility
# T is the time to expiration
# q is the continuous dividend yield
# call is a flag indicating whether the option is a call (1) or put (0)

def d1(S,K,r,sigma,T,q):
num = log(S/K)+(r-q+0.5*sigma**2)*T
den = sigma*sqrt(T)
return num/den

def d2(S,K,r,sigma,T,q):
num = log(S/K)+(r-q-0.5*sigma**2)*T
den = sigma*sqrt(T)
return num/den

# Cumulative distribution function for the standard normal distribution
def phi(x):
return (1.0 + erf(x / sqrt(2.0))) / 2.0

# Current value of the option given market parameters
def V(S,K,r,sigma,T,q,call):
if call == 1:
return S*exp(-q*T)*phi(d1(S,K,r,sigma,T,q)) - K*exp(-r*T)*phi(d2(S,K,r,sigma,T,q))
else:
return exp(-r*T)*K*phi(-d2(S,K,r,sigma,T,q)) - S*exp(-q*T)*phi(-d1(S,K,r,sigma,T,q))

# Delta is the instantaneous rate of change of the option value with respect
# to the price of the underlying
def delta(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
if call == 1:
return exp(-q*T)*scipy.stats.norm(0, 1).cdf(d1_new)
else:
return -exp(-q*T)*scipy.stats.norm(0, 1).cdf(-d1_new)

# Gamma is the instantaneous rate of change of delta with respect to the price
# of the underlying
def gamma(S,K,r,sigma,T,q,call):
# Same for both call and put
d1_new = d1(S,K,r,sigma,T,q)
return exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)/(S*sigma*sqrt(T))

# Vega measures the instantaneous rate of change of the option value with
# respect to volatility.
def vega(S,K,r,sigma,T,q,call):
# Same for both call and put
d2_new = d2(S,K,r,sigma,T,q)
return K*exp(-r*T)*scipy.stats.norm(0, 1).pdf(d2_new)*sqrt(T)/100

# Theta, or time decay, measures the instantaneous rate of change of the option
# value with respect to the passage of time.
def theta(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
first_term = exp(-q*T)*S*sigma*scipy.stats.norm(0, 1).pdf(d1_new)/(2*sqrt(T))
second_term = r*exp(-r*T)*K*phi(d2_new)
third_term = q*S*exp(-q*T)*phi(d1_new)
if call == 1:
return (-first_term - second_term + third_term)/365
else:
return (-first_term + second_term - third_term)/365

# Rho measures sensitivity to the interest rate. It is the instantaneous rate
# of change of the option value with respect to the risk-free interest rate.
def rho(S,K,r,sigma,T,q,call):
d2_new = d2(S,K,r,sigma,T,q)
if call == 1:
return K*T*exp(-r*T)*phi(d2_new)/100
else:
return -K*T*exp(-r*T)*phi(-d2_new)/100

# Higher order Greeks

# Charm, or delta decay, measures the instantaneous rate of change of delta
# with the passage of time
def charm(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
first_term = q*exp(-q*T)*phi(d1_new)
second_term = exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)*(2*(r-q)*T - d2_new*sigma*sqrt(T))/(2*T*sigma*sqrt(T))
if call == 1:
return first_term - second_term
else:
return - first_term - second_term

# Speed is the rate of change of gamma with respect to the price of the
# underlying
def speed(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
gamma_new = gamma(S,K,r,sigma,T,q,call)
return (-gamma_new/S)*(d1_new/(sigma*sqrt(T))+1)

# Zomma measures the rate of change of gamma with respect to changes in volatility
def zomma(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
gamma_new = gamma(S,K,r,sigma,T,q,call)
return gamma_new*(d1_new*d2_new - 1)/sigma

# Vomma measures second-order exposure to volatility. Vomma is the instantaneous
# rate of change of vega with respect to volatility
def vomma(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
vega_new = vega(S,K,r,sigma,T,q,call)
return vega_new*d1_new*d2_new/sigma

def initialize(context):
context.stock = sid(24)

def handle_data(context, data):
pass

This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.
5 responses

Improvements to the Black-Scholes library. We now include a function to compute the implied volatility of an underlying asset given observed option prices. Under Black-Scholes, there is a one-to-one mapping of volatility to option price given a fixed set of other market parameters. The implied volatility is the volatility that would make the price given by the model equal to the observed price.

The implied volatility can be viewed as investors' expectation of the future volatility of an asset. When investing in stocks, one may use the implied vol to measure the riskiness of an asset in the future and adjust a position accordingly.

Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
import numpy as np
from math import sqrt, erf, exp
from math import log as logg
import scipy.stats
from scipy.optimize import minimize, newton

# Library for valuation of European call and put options using the
# standard Black-Scholes model

# S is the current stock price
# K is the strike price
# r is the risk-free rate
# sigma is the annualized volatility
# T is the time to expiration
# q is the continuous dividend yield
# call is a flag indicating whether the option is a call (1) or put (0)

def d1(S,K,r,sigma,T,q):
num = logg(S/K)+(r-q+0.5*sigma**2)*T
den = sigma*sqrt(T)
return num/den

def d2(S,K,r,sigma,T,q):
num = logg(S/K)+(r-q-0.5*sigma**2)*T
den = sigma*sqrt(T)
return num/den

# Cumulative distribution function for the standard normal distribution
def phi(x):
return (1.0 + erf(x / sqrt(2.0))) / 2.0

# Current value of the option given market parameters
def V(S,K,r,sigma,T,q,call):
if call == 1:
return S*exp(-q*T)*phi(d1(S,K,r,sigma,T,q)) - K*exp(-r*T)*phi(d2(S,K,r,sigma,T,q))
else:
return exp(-r*T)*K*phi(-d2(S,K,r,sigma,T,q)) - S*exp(-q*T)*phi(-d1(S,K,r,sigma,T,q))

# Delta is the instantaneous rate of change of the option value with respect
# to the price of the underlying
def delta(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
if call == 1:
return exp(-q*T)*scipy.stats.norm(0, 1).cdf(d1_new)
else:
return -exp(-q*T)*scipy.stats.norm(0, 1).cdf(-d1_new)

# Gamma is the instantaneous rate of change of delta with respect to the price
# of the underlying
def gamma(S,K,r,sigma,T,q,call):
# Same for both call and put
d1_new = d1(S,K,r,sigma,T,q)
return exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)/(S*sigma*sqrt(T))

# Vega measures the instantaneous rate of change of the option value with
# respect to volatility.
def vega(S,K,r,sigma,T,q,call):
# Same for both call and put
d2_new = d2(S,K,r,sigma,T,q)
return K*exp(-r*T)*scipy.stats.norm(0, 1).pdf(d2_new)*sqrt(T)/100

# Theta, or time decay, measures the instantaneous rate of change of the option
# value with respect to the passage of time.
def theta(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
first_term = exp(-q*T)*S*sigma*scipy.stats.norm(0, 1).pdf(d1_new)/(2*sqrt(T))
second_term = r*exp(-r*T)*K*phi(d2_new)
third_term = q*S*exp(-q*T)*phi(d1_new)
if call == 1:
return (-first_term - second_term + third_term)/365
else:
return (-first_term + second_term - third_term)/365

# Rho measures sensitivity to the interest rate. It is the instantaneous rate
# of change of the option value with respect to the risk-free interest rate.
def rho(S,K,r,sigma,T,q,call):
d2_new = d2(S,K,r,sigma,T,q)
if call == 1:
return K*T*exp(-r*T)*phi(d2_new)/100
else:
return -K*T*exp(-r*T)*phi(-d2_new)/100

# Higher order Greeks

# Charm, or delta decay, measures the instantaneous rate of change of delta
# with the passage of time
def charm(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
first_term = q*exp(-q*T)*phi(d1_new)
second_term = exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)*(2*(r-q)*T - d2_new*sigma*sqrt(T))/(2*T*sigma*sqrt(T))
if call == 1:
return first_term - second_term
else:
return - first_term - second_term

# Speed is the rate of change of gamma with respect to the price of the
# underlying
def speed(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
gamma_new = gamma(S,K,r,sigma,T,q,call)
return (-gamma_new/S)*(d1_new/(sigma*sqrt(T))+1)

# Zomma measures the rate of change of gamma with respect to changes in volatility
def zomma(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
gamma_new = gamma(S,K,r,sigma,T,q,call)
return gamma_new*(d1_new*d2_new - 1)/sigma

# Vomma measures second-order exposure to volatility. Vomma is the instantaneous
# rate of change of vega with respect to volatility
def vomma(S,K,r,sigma,T,q,call):
d1_new = d1(S,K,r,sigma,T,q)
d2_new = d2(S,K,r,sigma,T,q)
vega_new = vega(S,K,r,sigma,T,q,call)
return vega_new*d1_new*d2_new/sigma

# Compute the implied volatility of the underlying instrument given the observed
# option prices in the market. We find a zero using the Newton-Rhapson method.
# The convergence rate of the Newton-Rhapson method is quadratic, meaning that
# if the function is well-behaved the actual error in the estimated zero is
# approximately the square of the requested tolerance.

def compute_implied_volatility(S,K,r,T,q,call,price):
sigma = 0.5
zero = newton(implied_vol_obj, sigma, fprime = None, args=(S,K,r,T,q,call,price), tol = 1e-25)
return zero

def implied_vol_obj(sigma,S,K,r,T,q,call,price):
return abs(V(S,K,r,sigma,T,q,call) - price)

def initialize(context):
context.stock = sid(24)

def handle_data(context, data):
# Example to show how to compute implied volatility
S = 51
K = 50
r = .05
sigma = 0.8
T = 0.167
q = 0
call = 0
price = V(S,K,r,sigma,T,q,call)
implied_vol = compute_implied_volatility(S,K,r,T,q,call,price)
log.info('The implied volatility is %s' % implied_vol)

This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.

hopefully we can get some historical option data or surface data to use these.

Andrew - is there anything more you would like to see on Quantopian in terms of option pricing models, algorithms for constructing implied volatility surfaces, etc.? I'd like to be able to improve the convergence rate of the optimization routines used to compute implied volatility.

I think you have most of the basics, another 2nd order greek would be vanna (rate of change of vega with respect to spot)

I'm hoping we'll one day get an interpolated vol surface something like
(secid, daysToMaturity, delta, impliedVol) for each date for the major etfs

also, realizedVol calculation that uses log of price returns like http://www.volx.us/VolFormula.htm, but also includes dividend and corp. action adjustments
as a simple transform like mavg(days)

@ Ryan how come The implied volatility is always 0.8 ? on the print out? can this algo different IV for each individual stocks thanks... ;)