This article will discuss the pricing of a digital call (and put) option using Monte Carlo methods. We've already seen how to do this for vanilla calls and puts. We will modify the code from the previous article to handle the pay-off function for digital options, which makes use of the Heaviside function.
Digital Options
Digital options are similar to vanilla options. They differ only in the fact that the pay-off at expiry, $f(T)$, only has two values. In this case $f(T) \in \{0,1\}$. In particular, the pay-off function is given as the Heaviside function with $S(T)-K$, the difference between the spot at expiry and the strike, as the argument:
\begin{eqnarray*} f(S(T)) = H(S(T)-K) \end{eqnarray*}where $H(x) = 1$ when $x\geq 0$ and $H(x)=0$ otherwise.
Thus it is apparent where the digital arises in the name!
C++ Implementation
Given that we have already considered the basic Monte Carlo approach in the article on pricing European vanilla calls and puts with Monte Carlo, I will only discuss the modifications to the code. However, I've still presented the full listing, which will give you everything you need to implement a basic digital option pricer.
#include <cmath>
#include <iostream>
// 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);
}
// This is the Heaviside step function, named after English
// mathematician Oliver Heaviside. It returns unity when val
// is greater than or equal to zero and returns zero otherwise
double heaviside(const double& val) {
if (val >= 0) {
return 1.0;
} else {
return 0.0;
}
}
// Pricing a digital call option with a Monte Carlo method
double monte_carlo_digital_call_price(const int& num_sims, const double& S, const double& K, const double& r, const double& v, const double& T) {
double S_adjust = S * exp(T*(r-0.5*v*v));
double S_cur = 0.0;
double payoff_sum = 0.0;
for (int i=0; i<num_sims; i++) {
double gauss_bm = gaussian_box_muller();
S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
payoff_sum += heaviside(S_cur - K);
}
return (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
}
// Pricing a digital put option with a Monte Carlo method
double monte_carlo_digital_put_price(const int& num_sims, const double& S, const double& K, const double& r, const double& v, const double& T) {
double S_adjust = S * exp(T*(r-0.5*v*v));
double S_cur = 0.0;
double payoff_sum = 0.0;
for (int i=0; i<num_sims; i++) {
double gauss_bm = gaussian_box_muller();
S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
payoff_sum += heaviside(K - S_cur);
}
return (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
}
int main(int argc, char **argv) {
// First we create the parameter list
int num_sims = 10000000; // Number of simulated asset paths
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 via Monte Carlo
double call = monte_carlo_digital_call_price(num_sims, S, K, r, v, T);
double put = monte_carlo_digital_put_price(num_sims, S, K, r, v, T);
// Finally we output the parameters and prices
std::cout << "Number of Paths: " << num_sims << std::endl;
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::cout << "Call Price: " << call << std::endl;
std::cout << "Put Price: " << put << std::endl;
return 0;
}
As before, we can see the Box-Muller function as well as the functions to price the options by the Monte Carlo method. However, I have added the Heaviside function, necessary for the pay-off for the digital options. In addition, I have removed the <algorithm> library from the include list as we do not need the max function for these pay-offs. The last modification alters the payoff_sum line to use the Heaviside function.
// This is the Heaviside step function, named after English
// mathematician Oliver Heaviside. It returns unity when val
// is greater than or equal to zero and returns zero otherwise
double heaviside(const double& val) {
if (val >= 0) {
return 1.0;
} else {
return 0.0;
}
}
And the modification to the monte_carlo_digital_call_price:
..
payoff_sum += heaviside(S_cur - K);
..
The output of the program is given as follows:
Number of Paths: 10000000
Underlying: 100
Strike: 100
Risk-Free Rate: 0.05
Volatility: 0.2
Maturity: 1
Call Price: 0.532424
Put Price: 0.418726
I will reiterate here that there is significant code duplication between this article and the article for vanilla calls and puts. I wanted to give you a full listing for digital options (and will do the same for all further variations of options) in order for those who haven't read the previous articles yet! Later on we will consider solutions to the code duplication issue by creating our own objects for different option pay-offs.
