One of the core financial applications of derivatives pricing theory is to be able to manage risk via a liquid options market. Such a market provides the capability for firms and individuals to tailor their risk exposure depending upon their hedging or speculation requirements. In order to effectively assess such risks, it is necessary to calculate the sensitivity of an options price to the factors that affect it, such as the underlying asset price, volatility and time to option expiry.
If we assume the existence of an analytical formula for an option price (such as in the case of European vanilla call/put options on a single asset) then it is possible to differentiate the call price with respect to its parameters in order to generate these sensitivities. The common sensitivities of interest include:
- Delta - Derivative of an option with respect to (w.r.t.) the spot price, $\frac{\partial C}{\partial S}$
- Gamma - Second derivative of an option w.r.t. the spot price, $\frac{\partial^2 C}{\partial S^2}$
- Vega - Derivative of an option w.r.t. the underlying volatility, $\frac{\partial C}{\partial \sigma}$
- Theta - (Negative) derivative of an option w.r.t. the time to expiry, $\frac{\partial C}{\partial t}$
- Rho - Derivative of an option w.r.t. the interest rate, $\frac{\partial C}{\partial \rho}$
Since all of the sensitivities are commonly denoted by letters of the Greek alphabet (except Vega!) they have come to be known colloquially as "the Greeks".
In this chapter we will calculate the Greeks using three separate methods. In the first instance we will utilise formula derived directly from the analytic formulae for European vanilla call and put options on a single asset. This will provide us with a baseline to determine the accuracy of subsequent numerical methods.
The first numerical approach utilised will be based on a Finite Difference Method (FDM) and the original analytical formulae. The second numerical method will use a combination of the FDM technique and Monte Carlo for pricing. The latter approach is readily applicable to a wider range of contingent claims as it is not dependent upon the existence of an analytic solution.
Analytic Formulae
The formulae of the Greeks for a European vanilla call and put option on a single asset are tabulated below:
| Calls | Puts | |
|---|---|---|
| Delta, $\frac{\partial C}{\partial S}$ | $N(d_1)$ | $N(d_1) - 1$ |
| Gamma, $\frac{\partial^{2} C}{\partial S^{2}}$ | $\frac{N'(d_1)}{S\sigma\sqrt{T - t}}\,$ | |
| Vega, $\frac{\partial C}{\partial \sigma}$ | $S N'(d_1) \sqrt{T-t}$ | |
| Theta, $\frac{\partial C}{\partial t}$ | $-\frac{S N'(d_1) \sigma}{2 \sqrt{T - t}} - rKe^{-r(T - t)}N(d_2)$ | $-\frac{S N'(d_1) \sigma}{2 \sqrt{T - t}} + rKe^{-r(T - t)}N(-d_2)$ |
| Rho, $\frac{\partial C}{\partial r}$ | $K(T - t)e^{-r(T - t)}N( d_2)$ | $-K(T - t)e^{-r(T - t)}N(-d_2)$ |
Where $N$ is the cumulative distribution function of the standard normal distribution, $N'$ is the probability density function of the standard normal distribution, $d_1 = (log(S/K)+(r+\frac{\sigma^2}{2})T)/(\sigma \sqrt{T})$ and $d_2 = d_1 ? \sigma \sqrt{T}$.
The following listing implements these methods in C++ making use of the formulae for the probability functions given in the article on European Option Pricing. It also includes a basic Monte Carlo pricer, which is taken from this article. Here is the listing for black_scholes.h
#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
#include <algorithm> // Needed for the "max" function
// =================
// ANALYTIC FORMULAE
// =================
// Standard normal probability density function
double norm_pdf(const double x) {
return (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x);
}
// An approximation to the cumulative distribution function
// for the standard normal distribution
// Note: This is a recursive function
double norm_cdf(const double x) {
double k = 1.0/(1.0 + 0.2316419*x);
double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k))));
if (x >= 0.0) {
return (1.0 - (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum);
} else {
return 1.0 - norm_cdf(-x);
}
}
// This calculates d_j, for j in {1,2}. This term appears in the closed
// form solution for the European call or put price
double d_j(const int j, const double S, const double K, const double r, const double v, const double T) {
return (log(S/K) + (r + (pow(-1,j-1))*0.5*v*v)*T)/(v*(pow(T,0.5)));
}
// Calculate the European vanilla call price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double call_price(const double S, const double K, const double r, const double v, const double T) {
return S * norm_cdf(d_j(1, S, K, r, v, T))-K*exp(-r*T) * norm_cdf(d_j(2, S, K, r, v, T));
}
// Calculate the European vanilla call Delta
double call_delta(const double S, const double K, const double r, const double v, const double T) {
return norm_cdf(d_j(1, S, K, r, v, T));
}
// Calculate the European vanilla call Gamma
double call_gamma(const double S, const double K, const double r, const double v, const double T) {
return norm_pdf(d_j(1, S, K, r, v, T))/(S*v*sqrt(T));
}
// Calculate the European vanilla call Vega
double call_vega(const double S, const double K, const double r, const double v, const double T) {
return S*norm_pdf(d_j(1, S, K, r, v, T))*sqrt(T);
}
// Calculate the European vanilla call Theta
double call_theta(const double S, const double K, const double r, const double v, const double T) {
return -(S*norm_pdf(d_j(1, S, K, r, v, T))*v)/(2*sqrt(T))
- r*K*exp(-r*T)*norm_cdf(d_j(2, S, K, r, v, T));
}
// Calculate the European vanilla call Rho
double call_rho(const double S, const double K, const double r, const double v, const double T) {
return K*T*exp(-r*T)*norm_cdf(d_j(2, S, K, r, v, T));
}
// Calculate the European vanilla put price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double put_price(const double S, const double K, const double r, const double v, const double T) {
return -S*norm_cdf(-d_j(1, S, K, r, v, T))+K*exp(-r*T) * norm_cdf(-d_j(2, S, K, r, v, T));
}
// Calculate the European vanilla put Delta
double put_delta(const double S, const double K, const double r, const double v, const double T) {
return norm_cdf(d_j(1, S, K, r, v, T)) - 1;
}
// Calculate the European vanilla put Gamma
double put_gamma(const double S, const double K, const double r, const double v, const double T) {
return call_gamma(S, K, r, v, T); // Identical to call by put-call parity
}
// Calculate the European vanilla put Vega
double put_vega(const double S, const double K, const double r, const double v, const double T) {
return call_vega(S, K, r, v, T); // Identical to call by put-call parity
}
// Calculate the European vanilla put Theta
double put_theta(const double S, const double K, const double r, const double v, const double T) {
return -(S*norm_pdf(d_j(1, S, K, r, v, T))*v)/(2*sqrt(T))
+ r*K*exp(-r*T)*norm_cdf(-d_j(2, S, K, r, v, T));
}
// Calculate the European vanilla put Rho
double put_rho(const double S, const double K, const double r, const double v, const double T) {
return -T*K*exp(-r*T)*norm_cdf(-d_j(2, S, K, r, v, T));
}
// ===========
// MONTE CARLO
// ===========
// A simple implementation of the Box-Muller algorithm, used to generate
// gaussian random numbers - necessary for the Monte Carlo method below
// Note that C++11 actually provides std::normal_distribution<> in
// the <random> library, which can be used instead of this function
double gaussian_box_muller() {
double x = 0.0;
double y = 0.0;
double euclid_sq = 0.0;
// Continue generating two uniform random variables
// until the square of their "euclidean distance"
// is less than unity
do {
x = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
y = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
euclid_sq = x*x + y*y;
} while (euclid_sq >= 1.0);
return x*sqrt(-2*log(euclid_sq)/euclid_sq);
}
Here is the listing for the main.cpp file making use of the above header:
#include "black_scholes.h"
int main(int argc, char **argv) {
// First we create the parameter list
double S = 100.0; // Option price
double K = 100.0; // Strike price
double r = 0.05; // Risk-free rate (5%)
double v = 0.2; // Volatility of the underlying (20%)
double T = 1.0; // One year until expiry
// Then we calculate the call/put values and the Greeks
double call = call_price(S, K, r, v, T);
double call_delta_v = call_delta(S, K, r, v, T);
double call_gamma_v = call_gamma(S, K, r, v, T);
double call_vega_v = call_vega(S, K, r, v, T);
double call_theta_v = call_theta(S, K, r, v, T);
double call_rho_v = call_rho(S, K, r, v, T);
double put = put_price(S, K, r, v, T);
double put_delta_v = put_delta(S, K, r, v, T);
double put_gamma_v = put_gamma(S, K, r, v, T);
double put_vega_v = put_vega(S, K, r, v, T);
double put_theta_v = put_theta(S, K, r, v, T);
double put_rho_v = put_rho(S, K, r, v, T);
// Finally we output the parameters and prices
std::cout << "Underlying: " << S << std::endl;
std::cout << "Strike: " << K << std::endl;
std::cout << "Risk-Free Rate: " << r << std::endl;
std::cout << "Volatility: " << v << std::endl;
std::cout << "Maturity: " << T << std::endl << std::endl;
std::cout << "Call Price: " << call << std::endl;
std::cout << "Call Delta: " << call_delta_v << std::endl;
std::cout << "Call Gamma: " << call_gamma_v << std::endl;
std::cout << "Call Vega: " << call_vega_v << std::endl;
std::cout << "Call Theta: " << call_theta_v << std::endl;
std::cout << "Call Rho: " << call_rho_v << std::endl << std::endl;
std::cout << "Put Price: " << put << std::endl;
std::cout << "Put Delta: " << put_delta_v << std::endl;
std::cout << "Put Gamma: " << put_gamma_v << std::endl;
std::cout << "Put Vega: " << put_vega_v << std::endl;
std::cout << "Put Theta: " << put_theta_v << std::endl;
std::cout << "Put Rho: " << put_rho_v << std::endl;
return 0;
}
The output of this program is as follows:
Underlying: 100
Strike: 100
Risk-Free Rate: 0.05
Volatility: 0.2
Maturity: 1
Call Price: 10.4506
Call Delta: 0.636831
Call Gamma: 0.018762
Call Vega: 37.524
Call Theta: -6.41403
Call Rho: 53.2325
Put Price: 5.57352
Put Delta: -0.363169
Put Gamma: 0.018762
Put Vega: 37.524
Put Theta: -1.65788
Put Rho: -41.8905
We are now going to compare the analytical prices with those derived from a Finite Difference Method.
Finite Difference Method
FDM is widely used in derivatives pricing (as well as engineering/physics in general) to solve partial differential equations (PDE). I have written before about using FDM to solve the Black-Scholes equation via the Explicit Euler Method. In this section we are going to apply the same technique, namely the discretisation of the partial derivatives, to create a simple approximation to the Greek sensitivities with which we can compare to the analytical solution.
The essence of the method is that we will approximate the partial derivative representing the particular sensitivity of interest. To do this we make use of the analytical formulae for the Black-Scholes prices of the call and puts. These formulae are covered in this article.
As an example, let's assume we want to calculate the Delta of a call option. The Delta is given by $\frac{\partial C}{\partial S} (S, T, \sigma, r, K)$. If we calculate two call prices, one at spot $S$ and the other at spot $S+\Delta S$, subtract the prices and divide by $\Delta S$, we have a forward difference approximation to the derivative:
\begin{eqnarray} \frac{\partial C}{\partial S} \approx \frac{C(S+\Delta S, T, \sigma, r, K) - C(S,T,\sigma,r,K)}{\Delta S} \end{eqnarray}Each of the additional first order sensitivities (Vega, Rho and Theta) can be calculated in this manner by simply incrementing the correct parameter dimension. Gamma on the other hand is a second order derivative and so must be approximated in a different way. The usual approach in FDM is to use a central difference approximation to produce the following formula:
\begin{eqnarray} \frac{\partial^2 C}{\partial S^2} \approx \frac{C(S+\Delta S, T, \sigma, r, K) - 2 C(S,T,\sigma,r,K) + C(S-\Delta S, T, \sigma, r, K)}{(\Delta S)^2} \end{eqnarray}At this stage we will keep the code procedural as we wish to emphasise the mathematical formulae. We are now able to implement the FDM numerical approximations in C++. For the sake of brevity we will restrict ourselves to the calculation of the call Delta and Gamma, as the remaining sensitivities are similar:
#include "black_scholes.h"
// This uses the forward difference approximation to calculate the Delta of a call option
double call_delta_fdm(const double S, const double K, const double r, const double v, const double T, const double delta_S) {
return (call_price(S + delta_S, K, r, v, T) - call_price(S, K, r, v, T))/delta_S;
}
// This uses the centred difference approximation to calculate the Gamma of a call option
double call_gamma_fdm(const double S, const double K, const double r, const double v, const double T, const double delta_S) {
return (call_price(S + delta_S, K, r, v, T) - 2*call_price(S, K, r, v, T) + call_price(S - delta_S, K, r, v, T))/(delta_S*delta_S);
}
int main(int argc, char **argv) {
// First we create the parameter list
double S = 100.0; // Option price
double delta_S = 0.001; // Option price increment
double K = 100.0; // Strike price
double r = 0.05; // Risk-free rate (5%)
double v = 0.2; // Volatility of the underlying (20%)
double T = 1.0; // One year until expiry
// Then we calculate the Delta and the Gamma for the call
double call_delta_f = call_delta_fdm(S, K, r, v, T, delta_S);
double call_gamma_f = call_gamma_fdm(S, K, r, v, T, delta_S);
// Finally we output the parameters and greeks
std::cout << "Underlying: " << S << std::endl;
std::cout << "Delta underlying: " << delta_S << std::endl;
std::cout << "Strike: " << K << std::endl;
std::cout << "Risk-Free Rate: " << r << std::endl;
std::cout << "Volatility: " << v << std::endl;
std::cout << "Maturity: " << T << std::endl << std::endl;
std::cout << "Call Delta: " << call_delta_f << std::endl;
std::cout << "Call Gamma: " << call_gamma_f << std::endl;
}
The output of this program is as follows:
Underlying: 100
Delta underlying: 0.001
Strike: 100
Risk-Free Rate: 0.05
Volatility: 0.2
Maturity: 1
Call Delta: 0.636845
Call Gamma: 0.0187633
Monte Carlo Method
The final method of calculating the Greeks is to use a combination of the FDM and Monte Carlo. The overall method is the same as above, with the exception that we will replace the analytical prices of the call/puts in the Finite Difference approximation and use a Monte Carlo engine instead to calculate the prices. This method is significantly more versatile as it can be extended to many differing types of contingent claim prices.
It is extremely important to re-use the random draws from the initial Monte Carlo simulation in calculating the incremented/decremented parameter paths. Thus, when calculating $C(S,K,r,v,T)$ and $C(S+\Delta S,K,r,v,T)$ it is necessary to use the same random draws. Further, it is significantly more optimal as there is no need to calculate additional asset paths for each incremented parameter instance.
The following code calculates the Monte Carlo price for the Delta and the Gamma, making use of separate Monte Carlo prices for each instance. The essence of the Monte Carlo method is to calculate three separate stock paths, all based on the same Gaussian draws. Each of these draws will represent an increment (or not) to the asset path parameter, $S$. Once those paths have been generated a price for each can be calculated. This price is then substituted into the FDM derivative approximations, in exactly the same way as with the analytical formulae.
The final prices are passed to the monte_carlo_call_price function by reference and the function itself is void. This provides a simple mechanism for returning multiple doubles from the function. An alternative would be to pass a vector of values by reference, to be set.
It is straightforward to modify the code to calculate the Vega, Rho or Theta (based on the Delta). A more 'production ready' implementation would utilise an object-oriented framework. However, for the purposes of this article an OOP-based implementation would likely obscure the algorithm:
#include "black_scholes.h"
// Pricing a European vanilla call option with a Monte Carlo method
// Create three separate paths, each with either an increment, non-
// increment or decrement based on delta_S, the stock path parameter
void monte_carlo_call_price(const int num_sims,
const double S, const double K, const double r,
const double v, const double T, const double delta_S,
double& price_Sp, double& price_S, double& price_Sm) {
// Since we wish to use the same Gaussian random draws for each path, it is
// necessary to create three separated adjusted stock paths for each
// increment/decrement of the asset
double Sp_adjust = (S+delta_S) * exp(T*(r-0.5*v*v));
double S_adjust = S * exp(T*(r-0.5*v*v));
double Sm_adjust = (S-delta_S) * exp(T*(r-0.5*v*v));
// These will store all three 'current' prices as the Monte Carlo
// algorithm is carried out
double Sp_cur = 0.0;
double S_cur = 0.0;
double Sm_cur = 0.0;
// There are three separate pay-off sums for the final price
double payoff_sum_p = 0.0;
double payoff_sum = 0.0;
double payoff_sum_m = 0.0;
// Loop over the number of simulations
for (int i=0; i<num_sims; i++) {
double gauss_bm = gaussian_box_muller(); // Random gaussian draw
// Adjust three stock paths
double expgauss = exp(sqrt(v*v*T)*gauss_bm); // Precalculate
Sp_cur = Sp_adjust * expgauss;
S_cur = S_adjust * expgauss;
Sm_cur = Sm_adjust * expgauss;
// Calculate the continual pay-off sum for each increment/decrement
payoff_sum_p += std::max(Sp_cur - K, 0.0);
payoff_sum += std::max(S_cur - K, 0.0);
payoff_sum_m += std::max(Sm_cur - K, 0.0);
}
// There are three separate prices
price_Sp = (payoff_sum_p / static_cast<double>(num_sims)) * exp(-r*T);
price_S = (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
price_Sm = (payoff_sum_m / static_cast<double>(num_sims)) * exp(-r*T);
}
double call_delta_mc(const int num_sims, const double S, const double K, const double r, const double v, const double T, const double delta_S) {
// These values will be populated via the monte_carlo_call_price function.
// They represent the incremented Sp (S+delta_S), non-incremented S (S) and
// decremented Sm (S-delta_S) prices.
double price_Sp = 0.0;
double price_S = 0.0;
double price_Sm = 0.0;
// Call the Monte Carlo pricer for each of the three stock paths
// (We only need two for the Delta)
monte_carlo_call_price(num_sims, S, K, r, v, T, delta_S, price_Sp, price_S, price_Sm);
return (price_Sp - price_S)/delta_S;
}
double call_gamma_mc(const int num_sims, const double S, const double K, const double r, const double v, const double T, const double delta_S) {
// These values will be populated via the monte_carlo_call_price function.
// They represent the incremented Sp (S+delta_S), non-incremented S (S) and
// decremented Sm (S-delta_S) prices.
double price_Sp = 0.0;
double price_S = 0.0;
double price_Sm = 0.0;
// Call the Monte Carlo pricer for each of the three stock paths
// (We need all three for the Gamma)
monte_carlo_call_price(num_sims, S, K, r, v, T, delta_S, price_Sp, price_S, price_Sm);
return (price_Sp - 2*price_S + price_Sm)/(delta_S*delta_S);
}
int main(int argc, char **argv) {
// First we create the parameter list
double S = 100.0; // Option price
double delta_S = 0.001; // Option price increment
double K = 100.0; // Strike price
double r = 0.05; // Risk-free rate (5%)
double v = 0.2; // Volatility of the underlying (20%)
double T = 1.0; // One year until expiry
int num_sims = 100000000; // Number of simulations to carry out for Monte Carlo
// Then we calculate the Delta and the Gamma for the call
double call_delta_m = call_delta_mc(num_sims, S, K, r, v, T, delta_S);
double call_gamma_m = call_gamma_mc(num_sims, S, K, r, v, T, delta_S);
// Finally we output the parameters and greeks
std::cout << "Number of sims: " << num_sims << std::endl;
std::cout << "Underlying: " << S << std::endl;
std::cout << "Delta underlying: " << delta_S << std::endl;
std::cout << "Strike: " << K << std::endl;
std::cout << "Risk-Free Rate: " << r << std::endl;
std::cout << "Volatility: " << v << std::endl;
std::cout << "Maturity: " << T << std::endl << std::endl;
std::cout << "Call Delta: " << call_delta_m << std::endl;
std::cout << "Call Gamma: " << call_gamma_m << std::endl;
}
Here is the output of the program:
Number of sims: 100000000
Underlying: 100
Delta underlying: 0.001
Strike: 100
Risk-Free Rate: 0.05
Volatility: 0.2
Maturity: 1
Call Delta: 0.636894
Call Gamma: 0.0188733
Here is a summary table with the calculation of the Delta and Gamma for a European vanilla call option across all of the methods listed above:
| Delta | Gamma | |
|---|---|---|
| Analytic | 0.636831 | 0.018762 |
| FDM/Analytic | 0.636845 | 0.0187633 |
| FDM/MC | 0.636894 | 0.0188733 |
The values are extremely similar, even for the Monte Carlo based approach, albeit at the computational expense of generating the random draws. In later articles it will be shown how the FDM/MC approach can be applied to other contingent claims where an analytical solution is not forthcoming.
