The Option Greek formulae express the change in the option price with respect to a parameter change taking as fixed all the other inputs. (Haug explores multiple parameter changes at once.) One significant use of Greek measures is to calibrate risk exposure. A market-making financial institution with a portfolio of options, for instance, would want a snap shot of its exposure to asset price, interest rates, dividend fluctuations. It would try to establish impacts of volatility and time decay. In the formulae below, the Greeks merely evaluate change to only one input at a time. In reality, we might expect a conflagration of changes in interest rates and stock prices etc.
No discussion of Black Scholes is complete without estimating the Greeks. For a more exhaustive list of Greek formulae that incorporate dividends please go to wikipedia. To understand the intuition behind Greeks and parameter sensitivities please view this playlist. Parameters are represented by the notation S, K, T, r, q and Sigma. See VBA code below. See Martin Haugh notes for formulae and clear exposition of same. Parameter Sensitivity can be estimated more precisely using the Black Scholes Greeks. The Greeks are important because we rely upon these to assist in hedging. The VBA code below can be used to estimate the effect of altering parameter inputs on the value of options. These are important for hedging and grasping the importance of variables changing. I have adapted here the coding style set by the late Simon Benninga. For a deep dive into estimating the Greeks using both analytical and numerically estimation please check out this playlist graphing included. Below we implement Analytical Greeks in VBA. Also see https://github.com/jrvarma/MSExcel-Black-Scholes
'S = Asset Price
'K = Option Exercise
'T = Option Maturity
'r = Risk free rate
'q = Dividend yield
'Sigma = Asset price return volatility
'Calculating d1
Function d1One(S, K, T, r, q, Sigma)
d1One = (Log(S / K) + (r - q) * T) / (Sigma * Sqr(T)) + 0.5 * Sigma * Sqr(T)
End Function
'Calculating d2
Function d2Two(S, K, T, r, q, Sigma)
d2Two = d1One(S, K, T, r, q, Sigma) - Sigma * Sqr(T)
End Function
'The normal probability density, N'(x), is defined as pdf
Function pdf(x)
pdf = Exp(-x ^ 2 / 2) / (Sqr(2 * Application.Pi()))
End Function
'The Black Schole Merton Call Function
' Application.NormSDist() calls up the built in excel function to estimate Normal Cumulative Probability
Function BSMCall(S, K, T, r, q, Sigma)
BSMCall = S * Exp(-q * T) * Application.NormSDist(d1One(S, K, T, r, q, Sigma)) - K * Exp(-T * r) * _
Application.NormSDist(d2Two(S, K, T, r, q, Sigma))
End Function
'The put is estimated using Put Call Parity
Function BSMPut(S, K, T, r, q, Sigma)
BSMPut = BSMCall(S, K, T, r, q, Sigma) + K * Exp(-r * T) - S * Exp(-q * T)
End Function
'The Delta for the Call
Function DeltaC(S, K, T, r, q, Sigma)
DeltaC = Exp(-q * T) * Application.NormSDist(d1One(S, K, T, r, q, Sigma))
End Function
'The Delta for the Put
Function DeltaP(S, K, T, r, q, Sigma)
DeltaP = -Exp(-q * T) * Application.NormSDist(-d1One(S, K, T, r, q, Sigma))
End Function
'The Gamma for the Call or Put
Function GammaCP(S, K, T, r, q, Sigma)
GammaCP = ((pdf(d1One(S, K, T, r, q, Sigma))) * Exp(-q * T)) / (S * Sigma * Sqr(T))
End Function
'The Vega for the Call or Put
Function Vega(S, K, T, r, q, Sigma)
Vega = S * Sqr(T) * pdf(d1One(S, K, T, r, q, Sigma)) * Exp(-q * T)
End Function
'The Theta for the Call
Function ThetaC(S, K, T, r, q, Sigma)
ThetaC = -S * pdf(d1One(S, K, T, r, q, Sigma)) * Sigma * Exp(-q * T) / (2 * Sqr(T)) + q * S * Application.NormSDist(d1One(S, K, T, r, q, Sigma)) * Exp(-q * T) - r * K * Exp(-r * T) * _
Application.NormSDist(d2Two(S, K, T, r, q, Sigma))
End Function
'The Theta for the Put
Function ThetaP(S, K, T, r, q, Sigma)
ThetaP = -S * pdf(d1One(S, K, T, r, q, Sigma)) * Sigma * Exp(-q * T) / (2 * Sqr(T)) - q * S * Application.NormSDist(-d1One(S, K, T, r, q, Sigma)) * _
Exp(-q * T) + r * K * Exp(-r * T) * Application.NormSDist(-d2Two(S, K, T, r, q, Sigma))
End Function
'The Rho for the Call
Function RhoC(S, K, T, r, q, Sigma)
RhoC = K * T * Exp(-r * T) * Application.NormSDist(d2Two(S, K, T, r, q, Sigma))
End Function
'The Rho for the Put
Function RhoP(S, K, T, r, q, Sigma)
RhoP = -K * T * Exp(-r * T) * Application.NormSDist(-d2Two(S, K, T, r, q, Sigma))
End Function
//Fabrice Rouah and Brian Byrne
#include <vector>
#include <stdio.h>
#include <math.h>
#include <iostream>
using namespace std;
// N(0,1) density
double f(double x) {
double pi = 4.0*atan(1.0);
return exp(-x*x*0.5) / sqrt(2 * pi);
}
// Boole's Rule
double Boole(double StartPoint, double EndPoint, int n) {
vector<double> X(n + 1, 0.0);
vector<double> Y(n + 1, 0.0);
double delta_x = (EndPoint - StartPoint) / double(n);
for (int i = 0; i <= n; i++) {
X[i] = StartPoint + i*delta_x;
Y[i] = f(X[i]);
}
double sum = 0;
for (int t = 0; t <= (n - 1) / 4; t++) {
int ind = 4 * t;
sum += (1 / 45.0)*(14 * Y[ind] + 64 * Y[ind + 1] + 24 * Y[ind + 2] + 64 * Y[ind + 3] + 14 * Y[ind + 4])*delta_x;
}
return sum;
}
// N(0,1) cdf by Boole's Rule
double N(double x) {
return Boole(-10.0, x, 240);
}
// Black-Scholes Call Price
double BSPrice(double S, double K, double T, double r, double q, double v, char OpType)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
double call = S* exp(-q*T)*N(d) - exp(-r*T)*K*N(d - v*sqrt(T));
if (OpType == 'C')
return call;
else
return call - S* exp(-q*T) + K*exp(-r*T);
}
// Black-Scholes Delta
double BSDelta(double S, double K, double T, double r, double q, double v, char OpType)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
if (OpType == 'C')
return exp(-q*T)*N(d);
else
return exp(-q*T)*(N(d) - 1);
}
// Black-Scholes Gamma
double BSGamma(double S, double K, double T, double r, double q, double v)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
return f(d) * exp(-q*T) / S / v / sqrt(T);
}
// Black-Scholes Vega
double BSVega(double S, double K, double T, double r, double q, double v)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
return S*f(d) * exp(-q*T) *sqrt(T);
}
// Black-Scholes Rho
double BSRho(double S, double K, double T, double r, double q, double v, char OpType)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
if (OpType == 'C')
return T*K*exp(-r*T)*N(d - v*sqrt(T));
else
return -T*K*exp(-r*T)*N(v*sqrt(T) - d);
}
// Black-Scholes Theta
double BSTheta(double S, double K, double T, double r, double q, double v, char OpType)
{
double d = (log(S / K) + T*(r - q + 0.5*v*v)) / (v*sqrt(T));
if (OpType == 'C')
return (-S*f(d) * exp(-q*T) *v / 2 / sqrt(T)) + q*exp(-q*T)*S*N(d) - r*K*exp(-r*T)*N(d - v*sqrt(T));
else
return (-S*f(d) * exp(-q*T) *v / 2 / sqrt(T)) - q*exp(-q*T)*S*N(-d) + r*K*exp(-r*T)*N(v*sqrt(T) - d);
}
int main()
{
double S = 100.0; // Stock Price
double K = 100.0; // Strike Price
double T = 1; // Years to maturity
double r = 0.05; // Risk free interest rate
double q = 0.0;
double v = 0.20; // Yearly volatility
char OpType = 'P'; // 'C'all or 'P'ut
cout << "Black Scholes Price " << BSPrice(S, K, T, r, q, v, OpType) << endl;
cout << "Black Scholes Delta " << BSDelta(S, K, T, r, q, v, OpType) << endl;
cout << "Black Scholes Gamma " << BSGamma(S, K, T, r, q, v) << endl;
cout << "Black Scholes Vega " << BSVega(S, K, T, r, q, v) << endl;
cout << "Black Scholes Rho " << BSRho(S, K, T, r, q, v, OpType) << endl;
cout << "Black Scholes Theta " << BSTheta(S, K, T, r, q, v, OpType)
<< " per year, or " << BSTheta(S, K, T, r, q, v, OpType) / 365.0 << " per day." << endl;
system("PAUSE");
};
The statistical software package R is one of the most promising tools for rapid prototyping of financial applications. Here we will introduce programs or functions written in R code to estimate the Black Scholes model, their Greeks (or parameter sensitivities) and Implied Volatility. The code is of a very high standard and follows Espen Haug. Haug E.G., The Complete Guide to Option Pricing Formulas. The R code is written to achieve high efficiency and optimization.
Robert McDonald, Professor of Finance at the Kellogg School of Management, Northwestern University and author of "Derivatives Markets" (Pearson, 2013) has put together the derivmkts package for R. Here I used the template for estimating Greeks as presented in Robert's derivmkts package and in his book's namesake.
library("derivmkts")
s=100; k=100; v=0.2; r=0.05; tt=1; d=0;
greeks(bscall(s, k, v, r, tt, d), complete=FALSE, long=FALSE, initcaps=TRUE)
greeks2(bscall, list(s=s, k=k, v=v, r=r, tt=tt, d=d))
greeks2(bscall, list(s=s, k=k, v=v, r=r, tt=tt, d=d))[c('Delta', 'Gamma'), ]
bsopt(s, k, v, r, tt, d)
bsopt(s, c(90, 100, 110), v, r, tt, d)
bsopt(s, c(90, 100, 110), v, r, tt, d)[['Call']][c('Delta', 'Gamma'), ]
## plot Greeks for calls and puts for 500 different stock prices
k <- 100; v <- 0.20; r <- 0.05; tt <- 1; d <- 0
S <- seq(10, 200, by=5)
Call <- greeks(bscall(S, k, v, r, tt, d))
Put <- greeks(bsput(S, k, v, r, tt, d))
y <- list(Call=Call, Put=Put)
par(mfrow=c(4, 4), mar=c(2, 2, 2, 2)) ## create a 4x4 plot
for (i in names(y)) {
for (j in rownames(y[[i]])) { ## loop over greeks
plot(S, y[[i]][j, ], main=paste(i, j), ylab=j, type='l')
}
}
Options greeks are mathematical derivatives of the option price with respect to inputs; see McDonald (2013), Chapters 12 and 13, for a discussion of the greeks for vanilla options. Greeks for vanilla and barrier options can be computed using the greeks function, which is a wrapper for any pricing function that returns the option price and which uses the default naming of inputs. In the Derivmkts package ,Prefessor McDonald has two alternative functions that return Greeks:
The bsopt function by default produces prices and Greeks for European calls and puts.
The greeks2 function takes as arguments the name of the pricing function and then inputs as a list
The Python code below was obtained from Davis Edwards github. We stripped back the full code to focus on estimating the Black Scholes Greeks and we follow the Generalized Black Scholes approach referenced in Espen Haug.
# This document demonstrates a Python implementation of some option models described in books written by Davis
# Edwards: "Energy Trading and Investing", "Risk Management in Trading", "Energy Investing Demystified".
# The gbs approach follows Espen Haug
# for backward compatability with Python 2.7
from __future__ import division
# import necessary libaries
import math
import numpy as np
from scipy.stats import norm
from scipy.stats import mvn
# The primary class for calculating Generalized Black Scholes option prices and deltas
# Inputs: option_type = "p" or "c", fs = price of underlying, x = strike, t = time to expiration, r = risk free rate
# b = cost of carry, v = implied volatility
# Outputs: value, delta, gamma, theta, vega, rho
def _gbs(option_type, fs, x, t, r, b, v):
# Create preliminary calculations
t__sqrt = math.sqrt(t)
d1 = (math.log(fs / x) + (b + (v * v) / 2) * t) / (v * t__sqrt)
d2 = d1 - v * t__sqrt
if option_type == "c":
# it's a call
# _debug(" Call Option")
value = fs * math.exp((b - r) * t) * norm.cdf(d1) - x * math.exp(-r * t) * norm.cdf(d2)
delta = math.exp((b - r) * t) * norm.cdf(d1)
gamma = math.exp((b - r) * t) * norm.pdf(d1) / (fs * v * t__sqrt)
theta = -(fs * v * math.exp((b - r) * t) * norm.pdf(d1)) / (2 * t__sqrt) - (b - r) * fs * math.exp(
(b - r) * t) * norm.cdf(d1) - r * x * math.exp(-r * t) * norm.cdf(d2)
vega = math.exp((b - r) * t) * fs * t__sqrt * norm.pdf(d1)
rho = x * t * math.exp(-r * t) * norm.cdf(d2)
else:
# it's a put
# _debug(" Put Option")
value = x * math.exp(-r * t) * norm.cdf(-d2) - (fs * math.exp((b - r) * t) * norm.cdf(-d1))
delta = -math.exp((b - r) * t) * norm.cdf(-d1)
gamma = math.exp((b - r) * t) * norm.pdf(d1) / (fs * v * t__sqrt)
theta = -(fs * v * math.exp((b - r) * t) * norm.pdf(d1)) / (2 * t__sqrt) + (b - r) * fs * math.exp(
(b - r) * t) * norm.cdf(-d1) + r * x * math.exp(-r * t) * norm.cdf(-d2)
vega = math.exp((b - r) * t) * fs * t__sqrt * norm.pdf(d1)
rho = -x * t * math.exp(-r * t) * norm.cdf(-d2)
return value, delta, gamma, theta, vega, rho
# This is the public interface for European Options
# Each call does a little bit of processing and then calls the calculations located in the _gbs module
# Inputs:
# option_type = "p" or "c"
# fs = price of underlying
# x = strike
# t = time to expiration
# v = implied volatility
# r = risk free rate
# q = dividend payment
# b = cost of carry
# Outputs:
# value = price of the option
# delta = first derivative of value with respect to price of underlying
# gamma = second derivative of value w.r.t price of underlying
# theta = first derivative of value w.r.t. time to expiration
# vega = first derivative of value w.r.t. implied volatility
# rho = first derivative of value w.r.t. risk free rates
# Black Scholes: stock Options (no dividend yield)
def black_scholes(option_type, fs, x, t, r, v):
b = r
return _gbs(option_type, fs, x, t, r, b, v)
# Merton Model: Stocks Index, stocks with a continuous dividend yields
def merton(option_type, fs, x, t, r, q, v):
b = r - q
return _gbs(option_type, fs, x, t, r, b, v)
def black_76(option_type, fs, x, t, r, v):
b = 0
return _gbs(option_type, fs, x, t, r, b, v)
# FX Options
def garman_kohlhagen(option_type, fs, x, t, r, rf, v):
b = r - rf
return _gbs(option_type, fs, x, t, r, b, v)
#black_scholes('c', fs=100, x=100, t=1, r=0.05, v=0.20)
# merton(option_type, fs, x, t, r, q, v):
# b = r - q
merton('c', fs=100, x=100, t=1, r=0.05, q = 0, v=0.20)
# c = black_scholes('c', fs=100, x=100, t=1, r=0.05, v=0.20)
# print( c )