Double digital option pricing with C++ via Monte Carlo methods

Double digital option pricing with C++ via Monte Carlo methods

This article will discuss the pricing of a double digital option using Monte Carlo methods. So far we've seen how to do this for vanilla calls and puts as well as digital calls and puts. As before, we will modify the code from the article on vanilla calls and puts pricing, by writing a double digital pay-off function.

Double Digital Options

Double digital options are similar to standard digitals, with the exception that they possess two "strike" prices, $K_L$ and $K_U$, with $K_L < K_U$. These represent the Lower and Upper bounds. If the spot at expiry, $S(T)$, of the underlying falls within the two strike values then the pay-off returns unity, else it returns zero.

\begin{eqnarray*} f(S(T)) = \begin{cases} 1, & \text{if }S(T) \in [K_L, K_U] \\ 0, & \text{if }S(T) \in (-\infty, K_L) \cup (K_U, \infty) \end{cases} \end{eqnarray*}

C++ Implementation

Following on from our approach to pricing digital options via Monte Carlo we will now modify the code presented in that article to handle double digital options. Since we do not distinguish between calls and puts we only need to write a single pay-off function. The Heaviside function has been removed and replaced with the payoff_double_digital function. Note that it takes three parameters: the two strike values and the spot price of the underlying asset at expiry. I've also replaced the payoff_sum line in the Monte Carlo function to make use of this pay-off function. Here is the full listing:

#define _USE_MATH_DEFINES

#include <cmath>
#include <iostream>

double gaussian_box_muller() {
  double x = 0.0;
  double y = 0.0;
  double euclid_sq = 0.0;

  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);
}

double payoff_double_digital(const double& K_L, const double& K_U, const double& x) {
  if (x >= K_L && x <= K_U) {
    return 1.0;
  } else {
    return 0.0;
  }
}

double monte_carlo_double_digital_price(const int& num_sims, const double& S, const double& K_L, const double& K_U, 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 += payoff_double_digital(K_L, K_U, 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_L = 100.0;  // Lower strike price                                                                                                                                                                                             
  double K_U = 120.0;  // Upper 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 double digital values                                                                                                                                                                                         
  double dd = monte_carlo_double_digital_price(num_sims, S, K_L, K_U, 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 << "Lower Strike:    " << K_L << std::endl;
  std::cout << "Upper Strike:    " << K_U << std::endl;
  std::cout << "Risk-Free Rate:  " << r << std::endl;
  std::cout << "Volatility:      " << v << std::endl;
  std::cout << "Maturity:        " << T << std::endl;

  std::cout << "Options Price:   " << dd << std::endl;

  return 0;
}

As before, we can see the Box-Muller function as well as the pay-off function needed to price the double digital option by the Monte Carlo method.

// The pay-off function for the double digital. It takes two strike values
// as well as the spot price of the underlying at expiry. It simply returns
// unity if the spot is within the strikes (inclusive) and zero otherwise.
double payoff_double_digital(const double& K_L, const double& K_U, const double& x) {
  if (x >= K_L && x <= K_U) {
    return 1.0;
  } else {
    return 0.0;
  }
}

And the modification to the monte_carlo_double_digital_price:

..
payoff_sum += payoff_double_digital(K_L, K_U, S_cur);
..

The output of the program is given as follows:

Number of Paths: 10000000
Underlying:      100
Lower Strike:    100
Upper Strike:    120
Risk-Free Rate:  0.05
Volatility:      0.2
Maturity:        1
Options Price:   0.32009

Once again, I will emphasise that we are continuing to duplicate the code among these articles, which is not best practice for production deployment! However, by providing the full listing above you can get started straight away pricing double digital options - if this is your goal. Note that in future articles we will create pay-off hierarchies using inheritance and polymorphism, which will allow us to resolve these issues of maintainability and duplication.