Heston Stochastic Volatility Model with Euler Discretisation in C++

Up until this point we have priced all of our options under the assumption that the volatility, $\sigma$, of the underlying asset has been constant over the lifetime of the option. In reality financial markets do not behave this way. Assets exist under market regimes where their volatility can vary signficantly during different time periods. The 2007-2008 financial crisis and the May Flash Crash are good examples of periods of intense market volatility.

Thus a natural extension of the Black Scholes model is to consider a non-constant volatility. Steven Heston formulated a model that not only considered a time-dependent volatility, but also introduced a stochastic (i.e. non-deterministic) component as well. This is the famous Heston model for stochastic volatility.

In this article we will outline the mathematical model and use a discretisation technique known as Full Truncation Euler Discretisation, coupled with Monte Carlo simulation, in order to price a European vanilla call option with C++. As with the majority of the models implemented on QuantStart, the code is object-oriented, allowing us to "plug-in" other option types (such as Path-Dependent Asians) with minimal changes.

Mathematical Model

The Black Scholes model uses a stochastic differential equation with a geometric Brownian motion to model the dynamics of the asset path. It is given by:

\begin{eqnarray} d S_t = \mu S_t dt + \sigma S_t dW^S_t \end{eqnarray}

Where I am using the notation of the Wikipedia Heston Model article. $S_t$ is the price of the underlying asset at time $t$, $\mu$ is the (constant) drift of the asset, $\sigma$ is the (constant) volatility of the underlying and $dW^S_t$ is a Weiner process (i.e. a random walk).

The Heston model extends this by introducing a second stochastic differential equation to represent the "path" of the volatility of the underlying over the lifetime of the option. The SDE for the variance is given by a Cox-Ingersoll-Ross process:

\begin{eqnarray} d S_t &=& \mu S_t dt + \sqrt{\nu_t} S_t dW^S_t \\ d \nu_t &=& \kappa (\theta - \nu_t) dt + \xi \sqrt{\nu_t} dW^\nu_t \end{eqnarray}

Where:

$\mu$ is the drift of the asset
$\theta$ is the expected value of $\nu_t$, i.e. the long run average price variance
$\kappa$ is the rate of mean reversion of $\nu_t$ to the long run average $\theta$
$\xi$ is the "vol of vol", i.e. the variance of $\nu_t$

Note that none of the parameters have any time-dependence. Extensions of the Heston model generally allow the values to become piecewise constant.

In order for $\nu_t > 0$, the Feller condition must be satisfied:

\begin{eqnarray} 2 \kappa \theta > \xi^2 \end{eqnarray}

In addition, the model enforces that the two separate Weiner processes making up the randomness are in fact correlated, with instantaneous constant correlation $\rho$:

\begin{eqnarray} dW^S_t dW^\nu_t &=& \rho dt \end{eqnarray}

Euler Discretisation

Given that the SDE for the asset path is now dependent (in a temporal manner) upon the solution of the second volatility SDE, it is necessary to simulate the volatility process first and then utilise this "volatility path" in order to simulate the asset path. In the case of the original Black Scholes SDE it is possible to use Ito's Lemma to directly solve for $S_t$. However, we are unable to utilise that procedure here and must use a numerical approximation in order to obtain both paths. The method utilised is known as Euler Discretisation.

The volatility path will be discretised into constant-increment time steps of $\Delta t$, with the updated volatility, $\nu_{i+1}$ given as an explicit function of $\nu_i$:

\begin{eqnarray} \nu_{i+1} &=& \nu_i + \kappa (\theta - \nu_i) \Delta t + \xi \sqrt{\nu_i} \Delta W^\nu_{i+1} \end{eqnarray}

However, since this is a finite discretisation of a continuous process, it is possible to introduce discretisation errors where $\nu_{i+1}$ can become negative. This is not a "physical" situation and so is a direct consequence of the numerical approximation. In order to handle negative values, we need to modify the above formula to include methods of eliminating negative values for subsequent iterations of the volatility path. Thus we introduce three new functions $f_1, f_2, f_3$, which lead to three separate "schemes" for how to handle the negative volatility values:

\begin{eqnarray} \nu_{i+1} &=& f_1(\nu_i) + \kappa (\theta - f_2(\nu_i)) \Delta t + \xi \sqrt{f_3(\nu_i)} \Delta W^\nu_{i+1} \end{eqnarray}

The three separate schemes are listed below:

Scheme	$f_1$	$f_2$	$f_3$
Reflection	$\|x\|$	$\|x\|$	$\|x\|$
Partial Truncation	$x$	$x$	$x^{+}$
Full Truncation	$x$	$x^{+}$	$x^{+}$

Where $x^{+} = \max(x, 0).$

The literature tends to suggest that the Full Truncation method is the "best" and so this is what we will utilise here. The Full Truncation scheme discretisation equation for the volatility path will thus be given by:

\begin{eqnarray} \nu_{i+1} &=& \nu_i + \kappa (\theta - \nu_i^{+}) \Delta t + \xi \sqrt{\nu_i^{+}} \Delta W^\nu_{i+1} \end{eqnarray}

In order to simulate $\Delta W^\nu_{i+1}$, we can make use of the fact that since it is a Brownian motion, $W^\nu_{i+1} - W^\nu_{i}$ is normally distributed with variance $\Delta t$ and that the distribution of $W^\nu_{i+1} - W^\nu_{i}$ is independent of $i$. This means it can be replaced with $\sqrt{\Delta t} N(0,1)$, where $N(0,1)$ is a random draw from the standard normal distribution.

We will return to the question of how to calculate the $W^\nu_{i}$ terms in the next section. Assuming we have the ability to do so, we are able to simulate the price of the asset path with the following discretisation:

\begin{eqnarray} S_{i+1} = S_i \exp \left(\mu - \frac{1}{2} v_i^{+}\right) \Delta t + \sqrt{v_i^{+}} \sqrt{\Delta t} \Delta W^S_{i+1} \end{eqnarray}

As with the Full Truncation mechanism outlined above, the volatility term appearing in the asset SDE discretisation has also been truncated and so $\nu_i$ is replaced by $\nu_i^{+}$.

Correlated Asset Paths

The next major issue that we need to look at is how to generate the $W^\nu_{i}$ and $W^S_{i}$ terms for the volatility path and the asset path respectively, such that they remain correlated with correlation $\rho$, as prescribed via the mathematical model. If you recall, I wrote an article about to how to generate correlated asset paths in C++ some time ago. In that article I described how to use two std::vector<> structures, filled with uniform random draws, in order to create a new vector of uniform draws which was correlated to the first with correlation $\rho$. This is exactly what is necessary here.

Once we have two uniform random draw vectors, it is possible to use the StandardNormalDistribution class, outlined in the article on statistical distributions in C++ to create two new vectors containing standard normal random draws - exactly what we need for the volatility and asset path simulation!

Monte Carlo Algorithm

In order to price a European vanilla call option under the Heston stochastic volatility model, we will need to generate many asset paths and then calculate the risk-free discounted average pay-off. This will be our option price.

If you would like to see this method (risk-neutral pricing by Monte Carlo) applied on a simpler model (Black Scholes), take a look at this article on European Option Pricing with Monte Carlo in C++.

The algorithm that we will follow to calculate the full options price is as follows:

Choose number of asset simulations for Monte Carlo and number of intervals to discretise asset/volatility paths over
For each Monte Carlo simulation, generate two uniform random number vectors, with the second correlated to the first
Use the statistics distribution class to convert these vectors into two new vectors containing standard normal draws
For each time-step in the discretisation of the vol path, calculate the next volatility value from the normal draw vector
For each time-step in the discretisation of the asset path, calculate the next asset value from the vol path vector and normal draw vector
For each Monte Carlo simulation, store the pay-off of the European call option
Take the mean of these pay-offs and then discount via the risk-free rate to produce an option price, under risk-neutral pricing.

We will now present a C++ implementation of this algorithm using a mixture of new code and prior classes written for other articles on the site. I will only reproduce new code in this article, but I will link to specific articles where the other classes can be found in their entirety.

C++ Implementation

We are going to take an object-oriented approach and break the calculation domain into various re-usable classes. In particular we will split the calculation into the following objects:

PayOff - This class represents an option pay-off object. We have discussed it at length on QuantStart.
Option - This class holds the parameters associated with the term sheet of the European option, as well as the risk-free rate. It requires a PayOff instance.
StandardNormalDistribution - This class allows us to create standard normal random draw values from a uniform distribution or random draws.
CorrelatedSND - This class takes two standard normal random draws and correlates the second with the first by a correlation factor $\rho$.
HestonEuler - This class accepts Heston model parameters and then performs a Full Truncation of the Heston model, generating both a volatility path and a subequent asset path.

We will now discuss the classes individually.

PayOff Class

I won't dwell on the PayOff class in any great detail within this article. It is described fully in the article on Path-Dependent Asians in C++. The PayOff class is a functor and as such is callable.

Option Class

The Option class is straightforward. It simply contains a set of public members for the option term sheet parameters (strike $K$, time to maturity $T$) as well as the (constant) risk-free rate $r$. The class also takes a pointer to a PayOff object, making it straightforward to "swap out" another pay-off (such as that for a Put option).

The listing for option.h follows:

#ifndef __OPTION_H
#define __OPTION_H

#include "payoff.h"

class Option {
 public:
  PayOff* pay_off;
  double K;
  double r;
  double T;

  Option(double _K, double _r, 
         double _T, PayOff* _pay_off);

  virtual ~Option();
};

#endif

The listing for option.cpp follows:

#ifndef __OPTION_CPP
#define __OPTION_CPP

#include "option.h"

Option::Option(double _K, double _r, 
               double _T, PayOff* _pay_off) : 
  K(_K), r(_r), T(_T), pay_off(_pay_off) {}

Option::~Option() {}

#endif

As can be seen from the above listings, the class doesn't do much beyond storing some data members and exposing them.

Statistics and CorrelatedSND Classes

The StandardNormalDistribution and CorrelatedSND classes are described in detail within the articles on statistical distributions in C++ and generating correlated asset paths in C++. Those articles contain the full listings for each of the classes, as well as the mathematics behind their implementation.

HestonEuler Class

The HestonEuler class is designed to accept the parameters of the Heston Model - in this case $\kappa$, $\theta$, $\xi$ and $\rho$ - and then calculate both the volatility and asset price paths. As such there are private data members for these parameters, as well as a pointer member representing the option itself. There are two calculation methods designed to accept the normal draw vectors and produce the respective volatility or asset spot paths.

The listing for heston_mc.h follows:

#ifndef __HESTON_MC_H
#define __HESTON_MC_H

#include <cmath>
#include <vector>
#include "option.h"

// The HestonEuler class stores the necessary information
// for creating the volatility and spot paths based on the
// Heston Stochastic Volatility model. 
class HestonEuler {
 private:
  Option* pOption;
  double kappa;
  double theta;
  double xi;
  double rho;

 public:
  HestonEuler(Option* _pOption, 
              double _kappa, double _theta, 
              double _xi, double _rho);
  virtual ~HestonEuler();

  // Calculate the volatility path
  void calc_vol_path(const std::vector<double>& vol_draws, 
                     std::vector<double>& vol_path);

  // Calculate the asset price path
  void calc_spot_path(const std::vector<double>& spot_draws, 
                      const std::vector<double>& vol_path, 
                      std::vector<double>& spot_path);
};

#endif

The listing for heston_mc.cpp follows:

#ifndef __HESTON_MC_CPP
#define __HESTON_MC_CPP

#include "heston_mc.h"

// HestonEuler
// ===========

HestonEuler::HestonEuler(Option* _pOption,
                         double _kappa, double _theta,
                         double _xi, double _rho) :
  pOption(_pOption), kappa(_kappa), theta(_theta), xi(_xi), rho(_rho) {}

HestonEuler::~HestonEuler() {}

void HestonEuler::calc_vol_path(const std::vector<double>& vol_draws, 
                                std::vector<double>& vol_path) { 
  size_t vec_size = vol_draws.size();
  double dt = pOption->T/static_cast<double>(vec_size);

  // Iterate through the correlated random draws vector and
  // use the 'Full Truncation' scheme to create the volatility path
  for (int i=1; i<vec_size; i++) {
    double v_max = std::max(vol_path[i-1], 0.0);
    vol_path[i] = vol_path[i-1] + kappa * dt * (theta - v_max) +
      xi * sqrt(v_max * dt) * vol_draws[i-1];
  }
}

void HestonEuler::calc_spot_path(const std::vector<double>& spot_draws, 
                                 const std::vector<double>& vol_path, 
                                 std::vector<double>& spot_path) {
  size_t vec_size = spot_draws.size();
  double dt = pOption->T/static_cast<double>(vec_size);

  // Create the spot price path making use of the volatility
  // path. Uses a similar Euler Truncation method to the vol path.
  for (int i=1; i<vec_size; i++) {
    double v_max = std::max(vol_path[i-1], 0.0);
    spot_path[i] = spot_path[i-1] * exp( (pOption->r - 0.5*v_max)*dt + 
        sqrt(v_max*dt)*spot_draws[i-1]);
  }
}

#endif

The calc_vol_path method takes references to a const vector of normal draws and a vector to store the volatility path. It calculates the $\Delta t$ value (as dt), based on the option maturity time. Then, the stochastic simulation of the volatility path is carried out by means of the Full Truncation Euler Discretisation, outlined in the mathematical treatment above. Notice that $\nu_i^{+}$ is precalculated, for efficiency reasons.

The calc_spot_path method is similar to the calc_vol_path method, with the exception that it accepts another vector, vol_path that contains the volatility path values at each time increment. The risk-free rate $r$ is obtained from the option pointer and, once again, $\nu_i^{+}$ is precalculated. Note that all vectors are passed by reference in order to reduce unnecessary copying.

Main Program

This is where it all comes together. There are two components to this listing: The generate_normal_correlation_paths function and the main function. The former is designed to handle the "boilerplate" code of generating the necessary uniform random draw vectors and then utilising the CorrelatedSND object to produce correlated standard normal distribution random draw vectors.

I wanted to keep this entire example of the Heston model tractable, so I have simply used the C++ built-in rand function to produce the uniform standard draws. However, in a production environment a Mersenne Twister uniform number generator (or something even more sophisticated) would be used to produce high-quality pseudo-random numbers. The output of the function is to calculate the values for the spot_normals and cor_normals vectors, which are used by the asset spot path and the volatility path respectively.

The main function begins by defining the parameters of the simulation, including the Monte Carlo values and those necessary for the option and Heston model. The actual parameter values are those give in the paper by Broadie and Kaya^[1]. The next task is to create the pointers to the PayOff and Option classes, as well as the HestonEuler instance itself.

After declaration of the various vectors used to hold the path values, a basic Monte Carlo loop is created. For each asset simulation, the new correlated values are generated, leading to the calculation of the vol path and the asset spot path. The option pay-off is calculated for each path and added to the total, which is then subsequently averaged and discounted via the risk-free rate. The option price is output to the terminal and finally the pointers are deleted.

Here is the listing for main.cpp:

#include <iostream>

#include "payoff.h"
#include "option.h"
#include "correlated_snd.h"
#include "heston_mc.h"

void generate_normal_correlation_paths(double rho, 
    std::vector<double>& spot_normals, std::vector<double>& cor_normals) {
  unsigned vals = spot_normals.size();

  // Create the Standard Normal Distribution and random draw vectors
  StandardNormalDistribution snd;
  std::vector<double> snd_uniform_draws(vals, 0.0);

  // Simple random number generation method based on RAND
  for (int i=0; i<snd_uniform_draws.size(); i++) {
    snd_uniform_draws[i] = rand() / static_cast<double>(RAND_MAX);
  }

  // Create standard normal random draws
  snd.random_draws(snd_uniform_draws, spot_normals);

  // Create the correlated standard normal distribution
  CorrelatedSND csnd(rho, &spot_normals);
  std::vector<double> csnd_uniform_draws(vals, 0.0);

  // Uniform generation for the correlated SND
  for (int i=0; i<csnd_uniform_draws.size(); i++) {
    csnd_uniform_draws[i] = rand() / static_cast<double>(RAND_MAX);
  }

  // Now create the -correlated- standard normal draw series
  csnd.random_draws(csnd_uniform_draws, cor_normals);
}

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 = 1000;  // Number of intervals for the asset path to be sampled 

  double S_0 = 100.0;    // Initial spot price
  double K = 100.0;      // Strike price
  double r = 0.0319;     // Risk-free rate
  double v_0 = 0.010201; // Initial volatility 
  double T = 1.00;       // One year until expiry

  double rho = -0.7;     // Correlation of asset and volatility
  double kappa = 6.21;   // Mean-reversion rate
  double theta = 0.019;  // Long run average volatility
  double xi = 0.61;      // "Vol of vol"

  // Create the PayOff, Option and Heston objects
  PayOff* pPayOffCall = new PayOffCall(K);
  Option* pOption = new Option(K, r, T, pPayOffCall);
  HestonEuler hest_euler(pOption, kappa, theta, xi, rho);

  // Create the spot and vol initial normal and price paths
  std::vector<double> spot_draws(num_intervals, 0.0);  // Vector of initial spot normal draws
  std::vector<double> vol_draws(num_intervals, 0.0);   // Vector of initial correlated vol normal draws
  std::vector<double> spot_prices(num_intervals, S_0);  // Vector of initial spot prices
  std::vector<double> vol_prices(num_intervals, v_0);   // Vector of initial vol prices

  // Monte Carlo options pricing
  double payoff_sum = 0.0;
  for (unsigned i=0; i<num_sims; i++) {
    std::cout << "Calculating path " << i+1 << " of " << num_sims << std::endl; 
    generate_normal_correlation_paths(rho, spot_draws, vol_draws);
    hest_euler.calc_vol_path(vol_draws, vol_prices);
    hest_euler.calc_spot_path(spot_draws, vol_prices, spot_prices);
    payoff_sum += pOption->pay_off->operator()(spot_prices[num_intervals-1]);
  }
  double option_price = (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);
  std::cout << "Option Price: " << option_price << std::endl;

  // Free memory
  delete pOption;
  delete pPayOffCall;

  return 0;
}

For completeness, I have included the makefile utilised on my MacBook Air, running Mac OSX 10.7.4:

heston: main.cpp heston_mc.o correlated_snd.o statistics.o option.o payoff.o
    clang++ -o heston main.cpp heston_mc.o correlated_snd.o statistics.o option.o payoff.o -arch x86_64

heston_mc.o: heston_mc.cpp option.o
    clang++ -c heston_mc.cpp option.o -arch x86_64

correlated_snd.o: correlated_snd.cpp statistics.o
    clang++ -c correlated_snd.cpp statistics.o -arch x86_64

statistics.o: statistics.cpp 
    clang++ -c statistics.cpp -arch x86_64

option.o: option.cpp payoff.o
    clang++ -c option.cpp payoff.o -arch x86_64 

payoff.o: payoff.cpp
    clang++ -c payoff.cpp -arch x86_64

Here is the output of the program:

..
..
Calculating path 99997 of 100000
Calculating path 99998 of 100000
Calculating path 99999 of 100000
Calculating path 100000 of 100000
Option Price: 6.81982

The exact option price is 6.8061, as reported by Broadie and Kaya^[1]. It can be made somewhat more accurate by increasing the number of asset paths and discretisation intervals.

There are a few extensions that could be made at this stage. One is to allow the various "schemes" to be implemented, rather than hard-coded as above. Another is to introduce time-dependence into the parameters. The next step after creating a model of this type is to actually calibrate to a set of market data such that the parameters may be determined. That will be the subject of another article!