European vanilla option pricing with C++ via Monte Carlo methods

European vanilla option pricing with C++ via Monte Carlo methods

In the previous article on using C++ to price a European option with analytic solutions we were able to take the closed-form solution of the Black-Scholes equation for a European vanilla call or put and provide a price.

This is possible because the boundary conditions generated by the pay-off function of the European vanilla option allow us to easily calculate a closed-form solution. Many option pay-off functions lead to boundary conditions which are much harder to solve analytically and some are impossible to solve this way. Thus in these situations we need to rely on numerical approximation.

In this article we will price the same European vanilla option with a very basic Monte Carlo solver in C++ and then compare our numerical values to the analytical case. We won't be concentrating on an extremely efficient or optimised implementation at this stage. Right now I just want to show you how to get up and running to give you an understanding of how risk neutral pricing works numerically. We will also see how our values differ from those generated analytically.

Risk Neutral Pricing of a European Vanilla Option

The first stage in implementation is to briefly discuss how risk neutral pricing works for a vanilla call or put option. We haven't yet discussed on QuantStart about risk neutral pricing in depth, so in the meantime, it might help to have a look at the article on Geometric Brownian Motion, which is the underlying model of stock price evolution that we will be using. It is given by the following stochastic differential equation:

\begin{eqnarray*} dS(t) = \mu S(t) dt + \sigma S(t) dB(t) \end{eqnarray*}

where $S$ is the asset price, $\mu$ is the drift of the stock, $\sigma$ is the volatility of the stock and $B$ is a Brownian motion (or Wiener process).

You can think of $dB$ as being a normally distributed random variable with zero mean and variance $dt$. We are going to use the fact that the price of a European vanilla option is given as the discounted expectation of the option pay-off of the final spot price (at maturity $T$):

\begin{eqnarray*} e^{-rT} \mathbb{E}(f(S(T))) \end{eqnarray*}

This expectation is taken under the appropriate risk-neutral measure, which sets the drift $\mu$ equal to the risk-free rate $r$:

\begin{eqnarray*} dS(t) = r S(t) dt + \sigma S(t) dB(t) \end{eqnarray*}

We can make use of Ito's Lemma to give us:

\begin{eqnarray*} d \log S(t) = \left(r - \frac{1}{2} \sigma^2 \right) dt + \sigma dB(t) \end{eqnarray*}

This is a constant coefficient SDE and therefore the solution is given by:

\begin{eqnarray*} \log S(t) = \log S(0) + \left(r - \frac{1}{2} \sigma^2 \right) T + \sigma \sqrt{T} N(0,1) \end{eqnarray*}

Here I've used the fact that since $B(t)$ is a Brownian motion, it has the distribution as a normal distribution with variance $T$ and mean zero so it can be written as $B(t) = \sqrt{T} N(0,1)$.

Rewriting the above equation in terms of $S(t)$ by taking the exponential gives:

\begin{eqnarray*} S(t) = S(0) e^{ (r - \frac{1}{2} \sigma^2) T + \sigma \sqrt{T} N(0,1) } \end{eqnarray*}

Using the risk-neutral pricing method above leads to an expression for the option price as follows:

\begin{eqnarray*} e^{-rT} \mathbb{E} (f(S(0) e^{ (r - \frac{1}{2} \sigma^2) T + \sigma \sqrt{T} N(0,1)} )) \end{eqnarray*}

The key to the Monte Carlo method is to make use of the law of large numbers in order to approximate the expectation. Thus the essence of the method is to compute many draws from the normal distribution $N(0,1)$, as the variable $x$, and then compute the pay-off via the following formula:

\begin{eqnarray*} f(S(0) e^{ (r - \frac{1}{2} \sigma^2) T + \sigma \sqrt{T} x} ) \end{eqnarray*}

In this case the value of $f$ is the pay-off for a call or put. By averaging the sum of these pay-offs and then taking the risk-free discount we obtain the approximate price for the option.

C++ Implementation

Let's now fire up your favourite C++ development environment and get started on the implementation. As with the previous article, we won't be utilising any object oriented or generic programming techniques right now. At this stage the goal is to generate understanding of a basic C++ Monte Carlo implementation. I've written out the program in full and then below I'll explain how each component works subsequently:

#include <algorithm>    // Needed for the "max" function
#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);
}

// Pricing a European vanilla call option with a Monte Carlo method
double monte_carlo_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 += std::max(S_cur - K, 0.0);
  }

  return (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
}

// Pricing a European vanilla put option with a Monte Carlo method
double monte_carlo_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 += std::max(K - S_cur, 0.0);
  }

  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_call_price(num_sims, S, K, r, v, T);
  double put = monte_carlo_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;
}

Let's now look at the program in a step-by-step fashion.

The first section imports the algorithm, iostream and cmath libraries. The algorithm library provides access to the max comparison function. This allows us to use the std::cout command to output code. In addition we now have access to the C mathematical functions such as exp, pow, log and sqrt:

#include <algorithm>
#include <cmath>
#include <iostream>

The first function to implement is the Box-Muller algorithm. As stated above, it is designed to convert two uniform random variables into a standard Gaussian random variable. Box-Muller is a good choice for a random number generator if your compiler doesn't support the C++11 standard. However, if you do have a compiler that supports it, you can make use of the std::normal_distribution<> template class found in the <random> library. In later articles, we will introduce this as a means of avoiding the need to "roll our own code":

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

The next two functions are used to price a European vanilla call or put via the Monte Carlo method. We've outlined above how the theory works, so let's study the implementation itself. The comments will give you the best overview. I've neglected to include the function for the put, as it is almost identical to the call and only differs in the nature of the pay-off. In fact, this is a hint telling us that our code is possibly repeating itself (a violation of the Do-Not-Repeat-Yourself, DRY, principle):

// Pricing a European vanilla call option with a Monte Carlo method
double monte_carlo_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));   // The adjustment to the spot price   
  double S_cur = 0.0;   // Our current asset price ("spot")
  double payoff_sum = 0.0;  // Holds the sum of all of the final option pay-offs

  for (int i=0; i<num_sims; i++) {
    double gauss_bm = gaussian_box_muller();   // Generate a Gaussian random number via Box-Muller
    S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);   // Adjust the spot price via the Brownian motion final distribution
    payoff_sum += std::max(S_cur - K, 0.0);   // Take the option pay-off, then add it to the rest of the pay-off values
  }

  // Average the pay-off sum via the number of paths and then 
  // discount the risk-free rate from the price
  return (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
}

The main function is extremely similar to that found in European vanilla option pricing with C++ and analytic formulae. The exception is the addition of the num_sims variable which stores the number of calculated asset paths. The larger this value, the more accurate the option price will be. Hence an inevitable tradeoff between execution time and price accuracy exists. Unfortunately, this is not a problem limited to this article! Another modification is that the call and put values are now derived from the Monte Carlo function calls:

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_call_price(num_sims, S, K, r, v, T);
  double put = monte_carlo_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;
}

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:      10.4553
Put Price:       5.57388

We can compare this with the output from the analytical formulae generated in European vanilla option pricing with C++ and analytic formulae, which are given below for convenience. This serves two purposes. Firstly, we can see that the values are almost identical, which provides an extra level of validation that the code in both instances is pricing correctly. Secondly, we can compare accuracy:

Call Price:      10.4506
Put Price:       5.57352

As can be seen the prices are relatively accurate for $10^7$ simulation paths. On my MacBook Air this program took a few seconds to calculate. Adding a couple of orders of magnitude on to the number of paths would quickly make the program prohibitive to execute. Thus we have run into the first problem - that of optimising execution speed. The second is that of making the above code maintainable. Notice that there is a significant duplication of code in both the monte_carlo_call_price and monte_carlo_put_price functions. Later articles will discuss these issues in depth and provide solutions to them.