Asian option pricing with C++ via Monte Carlo Methods

In this article I'm going to discuss how to price a certain type of *Exotic* option known as a **Path-Dependent Asian** in C++ using Monte Carlo Methods. It is considered "exotic" in the sense that the pay-off is a function of the underlying asset at multiple points throughout its lifetime, rather than just the value at expiry. An Asian option actually utilises the mean of the underlying asset price sampled at appropriate intervals as the basis for its pay-off, which is where the "path-dependency" of the asset comes from. The name actually arises because they were first devised in 1987 in Tokyo as options on crude oil futures.

There are two types of Asian option that we will be pricing. They are the *arithmetic Asian* and the *geometric Asian*. They differ only in how the mean of the underlying values is calculated. We will be studying discrete Asian options in this article, as opposed to the theoretical continuous Asian options. The first step will be to break down the components of the program and construct a set of classes that represent the various aspects of the pricer. Before that however, I will brief explain how Asian options work (and are priced by Monte Carlo).

An Asian option is a type of *exotic* option. Unlike a vanilla European option where the price of the option is dependent upon the price of the underlying asset at expiry, an Asian option pay-off is a function of multiple points up to and including the price at expiry. Thus it is *path-dependent* as the price relies on knowing how the underlying behaved at certain points *before* expiry. Asian options in particular base their price off the *mean average* price of these sampled points. To simplify this article we will consider $N$ equally distributed sample points beginning at time $t=0$ and ending at maturity, $T$.

Unlike in the vanilla European option Monte Carlo case, where we only needed to generate multiple spot values at expiry, we now need to generate multiple spot *paths*, each sampled at the correct points. Thus instead of providing a `double`

value representing spot to our option, we now need to provide a `std::vector<double>`

(i.e. a vector of double values), each element of which represents a sample of the spot price on a particular path. We will still be modelling our asset price path via a Geometric Brownian Motion (GBM), and we will create each path by adding the correct drift and variance at each step in order to maintain the properties of GBM.

In order to calculate the arithmetic mean $A$ of the spot prices we use the following formula:

\begin{eqnarray*} A(0,T) = \frac{1}{N} \sum^{N}_{i=1} S(t_i) \end{eqnarray*}Similarly for the geometric mean:

\begin{eqnarray*} A(0,T) = \exp\left(\frac{1}{N} \sum^{N}_{i=1} \log(S(t_i)) \right) \end{eqnarray*}These two formulae will then form the basis of our `pay_off_price`

method, which is attached to our `AsianOption`

class.

In order to increase maintainability we will separate the components of the Asian options pricer. As mentioned in the class on option pay-off hierarchies we are able to create an *abstract base class* called `PayOff`

, which defines an interface that all subsequent inherited pay-off classes will implement. The major benefit of this approach is that we can encapsulate multiple various types of pay-off functionality without the need to modify the remaining classes, such as our AsianOption class (to be discussed below). We make use of the `operator()`

to turn our pay-off classes into a *functor* (i.e. a function object). This means that we can "call" the object just as we would a function. "Calling" a `PayOff`

object has the effect of calculating the pay-off and returning it.

Here is the declaration for the `PayOff`

base class:

class PayOff { public: PayOff(); // Default (no parameter) constructor virtual ~PayOff() {}; // Virtual destructor // Overloaded () operator, turns the PayOff into an abstract function object virtual double operator() (const double& S) const = 0; };

The second class will represent many aspects of the exotic path-dependent Asian option. This class will be called `AsianOption`

. It will once again be an abstract base class (which means that it contains at least one *pure virtual function*). Since there are many types of Asian option - arithmetic, geometric, continuous, discrete - we will once again make use of the inheritance hierarchy. In particular, we will override the `pay_off_price`

function, which determines how the mean pay-off is to be calculated, once the appropriate `PayOff`

object has been provided.

Here is the declaration for the `AsianOption`

base class:

class AsianOption { protected: PayOff* pay_off; // Pay-off class (in this instance call or put) public: AsianOption(PayOff* _pay_off); virtual ~AsianOption() {}; // Pure virtual pay-off operator (this will determine arithmetic or geometric) virtual double pay_off_price(const std::vector<double>& spot_prices) const = 0; };

In C++ there is always a trade-off between simplicity and extensibility. More often than not a program must become more complicated if it is to be extendable elsewhere. Thus we have the choice of creating a separate object to handle the path generation used by the Asian option or write it procedurally. In this instance I have elected to use a procedural approach because I don't feel that delving into the complexities of random number generation classes has a beneficial effect on learning how Asian options are priced *at this stage*. In later articles we will encapsulate the random number and path generation, but right now we will **keep it simple**.

The first class we will consider is the `PayOff`

. As stated above this is an abstract base class and so can never be instantiated directly. Instead it provides an *interface* through which all inherited classes will be bound to. The destructor is `virtual`

to avoid memory leaks when the inherited and base classes are destroyed. We'll take a look at the full listing for the `payoff.h`

header file and then step through the important sections:

#ifndef __PAY_OFF_HPP #define __PAY_OFF_HPP #include <algorithm> // This is needed for the std::max comparison function, used in the pay-off calculations class PayOff { public: PayOff(); // Default (no parameter) constructor virtual ~PayOff() {}; // Virtual destructor // Overloaded () operator, turns the PayOff into an abstract function object virtual double operator() (const double& S) const = 0; }; class PayOffCall : public PayOff { private: double K; // Strike price public: PayOffCall(const double& K_); virtual ~PayOffCall() {}; // Virtual function is now over-ridden (not pure-virtual anymore) virtual double operator() (const double& S) const; }; class PayOffPut : public PayOff { private: double K; // Strike public: PayOffPut(const double& K_); virtual ~PayOffPut() {}; virtual double operator() (const double& S) const; }; #endif

The declaration for the `PayOff`

class is straightforward except for the overload of `operator()`

. The syntax says that this method cannot be implemented within this class (`= 0`

) and that it should be overridden by an inherited class (`virtual`

). It is also a `const`

method as it does not modify anything. It simply takes a spot price $S$ and returns the option pay-off, for a given strike, $K$:

// Overloaded () operator, turns the PayOff into an abstract function object virtual double operator() (const double& S) const = 0;

Beyond the `PayOff`

class we implement the `PayOffCall`

and `PayOffPut`

classes which implement call and put functionality. We could also introduce digital and double digital pay-off classes here, but for the purposes of demonstrating Asian options, I will leave that as an exercise! Notice now that in each of the declarations of the `operator()`

the `= 0`

pure virtual declaration has been removed. This states that the methods should be implemented by these classes:

// Virtual function is now over-ridden (not pure-virtual anymore) virtual double operator() (const double& S) const;

That sums up the header file for the `PayOff`

class hierarchy. We will now look at the source file (below) and then step through the important sections:

#ifndef __PAY_OFF_CPP #define __PAY_OFF_CPP #include "payoff.h" PayOff::PayOff() {} // ========== // PayOffCall // ========== // Constructor with single strike parameter PayOffCall::PayOffCall(const double& _K) { K = _K; } // Over-ridden operator() method, which turns PayOffCall into a function object double PayOffCall::operator() (const double& S) const { return std::max(S-K, 0.0); // Standard European call pay-off } // ========= // PayOffPut // ========= // Constructor with single strike parameter PayOffPut::PayOffPut(const double& _K) { K = _K; } // Over-ridden operator() method, which turns PayOffPut into a function object double PayOffPut::operator() (const double& S) const { return std::max(K-S, 0.0); // Standard European put pay-off } #endif

The only major point of note within the source file is the implementation of the call and put `operator()`

methods. They make use of the `std::max`

function found within the `<algorithm>`

standard library:

// Over-ridden operator() method, which turns PayOffCall into a function object double PayOffCall::operator() (const double& S) const { return std::max(S-K, 0.0); // Standard European call pay-off } .. // Over-ridden operator() method, which turns PayOffPut into a function object double PayOffPut::operator() (const double& S) const { return std::max(K-S, 0.0); // Standard European put pay-off }

The concludes the `PayOff`

class hierarchy. It is quite straightforward and we're not using any real advanced features beyond abstract base classes and pure virtual functions. However, we have written quite a lot of code to encapsulate something seemingly as simple as a simple pay-off function. The benefits of such an approach will become clear later on, when or if we decide we wish to create more complex pay-offs as we will not need to modify our `AsianOption`

class to do so.

To generate the paths necessary for an Asian option we will use a procedural approach. Instead of encapsulating the random number generator and path generator in a set of objects, we will make use of two functions, one of which generates the random Gaussian numbers (via the Box-Muller method) and the other which takes these numbers and generates sampled Geometric Brownian Motion asset paths for use in the `AsianOption`

object. The article on European option pricing via Monte Carlo explains the concept of risk-neutral pricing, Monte Carlo techniques and the Box-Muller method. I have included the full function in the listing below for completeness.

The second function, `calc_path_spot_prices`

is more relevant. It takes a reference to a vector of doubles (`std::vector<double>&`

) and updates it to reflect a random spot path based on a set of random Gaussian increments of the asset price. I will show the full listing for the header file which contains these functions, then I will step through `calc_path_spot_prices`

in detail:

#ifndef __PATH_GENERATION_H #define __PATH_GENERATION_H #include <vector> #include <cmath> // For random Gaussian generation // Note that there are many ways of doing this, but we will // be using the Box-Muller method for demonstration purposes 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); } // This provides a vector containing sampled points of a // Geometric Brownian Motion stock price path void calc_path_spot_prices(std::vector<double>& spot_prices, // Vector of spot prices to be filled in const double& r, // Risk free interest rate (constant) const double& v, // Volatility of underlying (constant) const double& T) { // Expiry // Since the drift and volatility of the asset are constant // we will precalculate as much as possible for maximum efficiency double dt = T/static_cast<double>(spot_prices.size()); double drift = exp(dt*(r-0.5*v*v)); double vol = sqrt(v*v*dt); for (int i=1; i<spot_prices.size(); i++) { double gauss_bm = gaussian_box_muller(); spot_prices[i] = spot_prices[i-1] * drift * exp(vol*gauss_bm); } } #endif

Turning our attention to `calc_path_spot_prices`

we can see that the function requires a vector of doubles, the risk-free rate $r$, the volatility of the underlying $\sigma$ and the time at expiry, $T$:

void calc_path_spot_prices(std::vector<double>& spot_prices, // Vector of spot prices to be filled in const double& r, // Risk free interest rate (constant) const double& v, // Volatility of underlying (constant) const double& T) { // Expiry

Since we are dealing with constant increments of time for our path sampling frequency we need to calculate this increment. These increments are always identical so in actual fact it can be pre-calculated outside of the loop for the spot price path generation. Similarly, via the properties of Geometric Brownian Motion, we know that we can increment the drift and variance of the asset in a manner which can be pre-computed. The only difference between this increment and the European case is that we are replacing $T$ with $dt$ for each subsequent increment of the path. See the European option pricing article for a comparison:

double dt = T/static_cast<double>(spot_prices.size()); double drift = exp(dt*(r-0.5*v*v)); double vol = sqrt(v*v*dt);

The final part of the function calculates the new spot prices by iterating over the `spot_price`

vector and adding the drift and variance to each piece. We are using the arithmetic of logarithms here, and thus can *multiply* by our drift and variance terms, since it is the *log of the asset price* that is subject to normally distributed increments in Geometric Brownian Motion. Notice that the loop runs from $i=1$, not $i=0$. This is because the `spot_price`

vector is already pre-filled with $S$, the initial spot, elsewhere in the program:

for (int i=1; i<spot_prices.size(); i++) { double gauss_bm = gaussian_box_muller(); spot_prices[i] = spot_prices[i-1] * drift * exp(vol*gauss_bm); }

That concludes the path-generation header file. Now we'll take a look at the `AsianOption`

classes themselves.

The final component of our program (besides the `main`

file of course!) is the Asian option inheritance hierarchy. We wish to price multiple types of Asian option, including geometric Asian options and arithmetic Asian options. One way to achieve this is to have separate methods on an `AsianOption`

class. However this comes at the price of having to continually add more methods if we wish to make more granular changes to our `AsianOption`

class. In a production environment this would become unwieldy. Instead we can use the abstract base class approach and generate an abstract `AsianOption`

class with a pure virtual method for `pay_off_price`

. This method is implemented in subclasses and determines how the averaging procedure for the asset prices over the asset path lifetime will occur. Two publicly-inherited subclasses `AsianOptionArithmetic`

and `AsianOptionGeometric`

implement this method.

Let's take a look at the full listing for the header declaration file and then we'll step through the interesting sections:

#ifndef __ASIAN_H #define __ASIAN_H #include <vector> #include "payoff.h" class AsianOption { protected: PayOff* pay_off; // Pay-off class (in this instance call or put) public: AsianOption(PayOff* _pay_off); virtual ~AsianOption() {}; // Pure virtual pay-off operator (this will determine arithmetic or geometric) virtual double pay_off_price(const std::vector<double>& spot_prices) const = 0; }; class AsianOptionArithmetic : public AsianOption { public: AsianOptionArithmetic(PayOff* _pay_off); virtual ~AsianOptionArithmetic() {}; // Override the pure virtual function to produce arithmetic Asian Options virtual double pay_off_price(const std::vector<double>& spot_prices) const; }; class AsianOptionGeometric : public AsianOption { public: AsianOptionGeometric(PayOff* _pay_off); virtual ~AsianOptionGeometric() {}; // Overide the pure virtual function to produce geometric Asian Options virtual double pay_off_price(const std::vector<double>& spot_prices) const; }; #endif

The first thing to notice about the abstract `AsianOption`

class is that it has a pointer to a `PayOff`

class as a protected member:

protected: PayOff* pay_off; // Pay-off class (in this instance call or put)

This is an example of *polymorphism*. The object will not know what type of `PayOff`

class will be passed in. It could be a `PayOffCall`

or a `PayOffPut`

. Thus we can use a pointer to the `PayOff`

abstract class to represent "storage" of this as-yet-unknown `PayOff`

class. Note also that the constructor takes this as its only parameter:

public: AsianOption(PayOff* _pay_off);

Finally we have the declaration for the pure virtual `pay_off_price`

method. This takes a vector of spot prices, but notice that it is passed as a reference to `const`

, so this vector will not be modified within the method:

// Pure virtual pay-off operator (this will determine arithmetic or geometric) virtual double pay_off_price(const std::vector<double>& spot_prices) const = 0;

The listings for the `AsianOptionArithmetic`

and `AsianOptionGeometric`

classes are analogous to those in the `PayOff`

hierarchy, with the exception that their constructors take a pointer to a `PayOff`

object. That concludes the header file.

The source file essentially implements the two `pay_off_price`

methods for the inherited subclasses of `AsianOption`

:

#ifndef __ASIAN_CPP #define __ASIAN_CPP #include <numeric> // Necessary for std::accumulate #include <cmath> // For log/exp functions #include "asian.h" // ===================== // AsianOptionArithmetic // ===================== AsianOption::AsianOption(PayOff* _pay_off) : pay_off(_pay_off) {} // ===================== // AsianOptionArithmetic // ===================== AsianOptionArithmetic::AsianOptionArithmetic(PayOff* _pay_off) : AsianOption(_pay_off) {} // Arithmetic mean pay-off price double AsianOptionArithmetic::pay_off_price(const std::vector& spot_prices) const { unsigned num_times = spot_prices.size(); double sum = std::accumulate(spot_prices.begin(), spot_prices.end(), 0); double arith_mean = sum / static_cast<double>(num_times); return (*pay_off)(arith_mean); } // ==================== // AsianOptionGeometric // ==================== AsianOptionGeometric::AsianOptionGeometric(PayOff* _pay_off) : AsianOption(_pay_off) {} // Geometric mean pay-off price double AsianOptionGeometric::pay_off_price(const std::vector & spot_prices) const { unsigned num_times = spot_prices.size(); double log_sum = 0.0; for (int i=0; i<spot_prices.size(); i++) { log_sum += log(spot_prices[i]); } double geom_mean = exp(log_sum / static_cast<double>(num_times) ); return (*pay_off)(geom_mean); } #endif

Let's take a look at the arithmetic option version. First of all we determine the number of sample points via the size of the `spot_price`

vector. Then we use the `std::accumulate`

algorithm and iterator syntax to sum the spot values in the vector. Finally we take the arithmetic mean of those values and use pointer dereferencing to call the `operator()`

for the `PayOff`

object. For this program it will provide a call or a put pay-off function for the average of the spot prices:

double AsianOptionArithmetic::pay_off_price(const std::vector& spot_prices) const { unsigned num_times = spot_prices.size(); double sum = std::accumulate(spot_prices.begin(), spot_prices.end(), 0); double arith_mean = sum / static_cast<double>(num_times); return (*pay_off)(arith_mean); }

The geometric Asian is similar. We once again determine the number of spot prices. Then we loop over the spot prices, summing the logarithm of each of them and adding it to the grand total. Then we take the geometric mean of these values and finally use pointer dereferencing once again to determine the correct call/put pay-off value:

double AsianOptionGeometric::pay_off_price(const std::vector& spot_prices) const { unsigned num_times = spot_prices.size(); double log_sum = 0.0; for (int i=0; i<spot_prices.size(); i++) { log_sum += log(spot_prices[i]); } double geom_mean = exp(log_sum / static_cast<double>(num_times) ); return (*pay_off)(geom_mean); }

Note here what the `AsianOption`

classes *do not* require. Firstly, they don't require information about the underlying (i.e. vol). They also don't require time to expiry or the interest rate. Thus we are really trying to encapsulate the *term sheet* of the option in this object, i.e. all of the parameters that would appear on the contract when the option is made. However, for simplicity we have neglected to include the actual *sample times*, which would also be written on the contract. This instead is moved to the path-generator. However, later code will amend this, particularly as we can re-use the `AsianOption`

objects in more sophisticated programs, where interest rates and volatility are subject to stochastic models.

The final part of the program is the `main.cpp`

file. It brings all of the previous components together to produce an output for the option price based on some default parameters. The full listing is below:

#include <iostream> #include "payoff.h" #include "asian.h" #include "path_generate.h" int main(int argc, char **argv) { // First we create the parameter list // Note that you could easily modify this code to input the parameters // either from the command line or via a file unsigned num_sims = 100000; // Number of simulated asset paths unsigned num_intervals = 250; // Number of intervals for the asset path to be sampled double S = 30.0; // Option price double K = 29.0; // Strike price double r = 0.08; // Risk-free rate (8%) double v = 0.3; // Volatility of the underlying (30%) double T = 1.00; // One year until expiry std::vector<double> spot_prices(num_intervals, S); // The vector of spot prices // Create the PayOff objects PayOff* pay_off_call = new PayOffCall(K); // Create the AsianOption objects AsianOptionArithmetic asian(pay_off_call); // Update the spot price vector with correct // spot price paths at constant intervals double payoff_sum = 0.0; for (int i=0; i<num_sims; i++) { calc_path_spot_prices(spot_prices, r, v, T); payoff_sum += asian.pay_off_price(spot_prices); } double discount_payoff_avg = (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T); delete pay_off_call; // Finally we output the parameters and prices std::cout << "Number of Paths: " << num_sims << std::endl; std::cout << "Number of Ints: " << num_intervals << 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 << "Asian Price: " << discount_payoff_avg << std::endl; return 0; }

The first interesting aspect of the `main.cpp`

program is that we have now added a `unsigned num_intervals = 250;`

line which determines how frequently the spot price will be sampled in the Asian option. As stated above, this would usually be incorporated into the option *term sheet*, but I have included it here instead to help make the pricer easier to understand without too much object communication overhead. We have also created the vector of spot prices, which are pre-filled with the default spot price, $S$:

unsigned num_intervals = 250; // Number of intervals for the asset path to be sampled .. std::vector<double> spot_prices(num_intervals, S); // The vector of spot prices

Then we create a `PayOffCall`

object and assign it to a `PayOff`

pointer. This allows us to leverage *polymorphism* and pass that object through to the `AsianOption`

, without the option needing to know the actual type of `PayOff`

. *Note that whenever we use the new operator, we must make use of the corresponding delete operator.*

PayOff* pay_off_call = new PayOffCall(K); .. delete pay_off_call;

The next step is to create the `AsianOptionArithmetic`

object. We could have as easily chosen the `AsianOptionGeometric`

and the program would be trivial to modify to do so. It takes in the pointer to the `PayOff`

as its lone constructor argument:

AsianOptionArithmetic asian(pay_off_call);

Then we create a loop for the total number of path simulations. In the loop we recalculate a new spot price path and then add that pay-off to a running sum of all pay-offs. The final step is to discount the average of this pay-off via the risk-free rate across the lifetime of the option ($T-0 = T$). This discounted price is then the final price of the option, with the above parameters:

double payoff_sum = 0.0; for (int i=0; i<num_sims; i++) { calc_path_spot_prices(spot_prices, r, v, T); payoff_sum += asian.pay_off_price(spot_prices); } double discount_payoff_avg = (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);

The final stage of the main program is to output the parameters and the options price:

Number of Paths: 100000 Number of Ints: 250 Underlying: 30 Strike: 29 Risk-Free Rate: 0.08 Volatility: 0.3 Maturity: 1 Asian Price: 2.85425

The benefits of such an object-oriented approach are now clear. We can easily add more `PayOff`

or `AsianOption`

classes without needing to extensively modify any of the remaining code. These objects are said to have a *separation of concerns*, which is exactly what is needed for large-scale software projects.

That is not to say that the program could not be improved! One obvious task is to determine how to incorporate the number of stock samples - and interval spacing - within the `AsianOption`

hierarchy. Another is to encapsulate the random number and path generation into its own "engine" that can be used in other option pricers.

As always, if you have any difficulty following the above, please feel free to email me at **mike@quantstart.com**.

You'll get instant access to a free 10-part email course packed with hints and tips to help you get started in quantitative trading!

Every week I'll send you a wrap of all activity on QuantStart so you'll never miss a post again.

Real, actionable quant trading tips with no nonsense.