Welcome back! In this post, I’ll show how classes from QuantLib’s Monte Carlo framework can be utilized to model the path an asset’s price can take over a period of time, such as that depicted in the one-year Intel (INTC) stock chart below. Why would we want a model of asset prices? There are a number of good reasons. One is that a model of asset price dynamics is essential to the valuation of derivatives, such as equity and index options. Secondly, such a model is a powerful tool for risk management. Simulated asset prices can be used to create a range of what-if scenarios with which to calculate a portfolio’s aggregate market risk exposure as measured by metrics such as Value at Risk (VaR). And thirdly, simulated prices that conform to historical asset return parameters (e.g. annualized mean and standard deviation) can be employed as market data for back-testing trading strategies.

So now that we appreciate the why, let’s look at the how. To start, let’s briefly address the theoretical foundations of asset price dynamics- the s*tochastic process*. Informally, a stochastic process is a function of one or more time varying parameters where at least one of the parameters is non-deterministic; its values correspond to a sequence of independent random variables drawn from a selected probability distribution. More precisely, the particular stochastic process that describes the evolution of stock prices is termed Geometric Brownian Motion (GBM).

Geometric Brownian Motion can be formulated as a Stochastic Differential Equation (SDE) of the form:

where *S* is the stock price at time* t*, μ (*mu)* represents the constant drift or trend (i.e. annual return) of the process and σ (*sigma) *represents the amount of random variation around the trend (i.e. annualized standard deviation of log returns). Intuitively, *mu* can be viewed as the ‘signal’, while *sigma* is the ‘noise’ of the GBM stochastic process. What this equation tells us is that the change in a stock’s price over a small discrete time increment (*dt*) is a function of the stock’s return (*mu*) and the stock’s volatility (*sigma*), where the volatility is scaled by the output of a Wiener process (*dWt*). The Wiener process essentially provides random numbers in accordance with a given (usually Gaussian) probability distribution. For more information on Geometric Brownian Motion, I encourage you to check out the Wikipedia entry here.

### Normally Distributed Model of Asset Returns

With that brief exposition of the theory of Geometric Brownian Motion behind us, let’s proceed with the QuantLib C++ code. First, I’ll present code that implements the classical model, in which stock returns are assumed to be normally distributed, even though empirical evidence demonstrates that asset returns are, in fact, ‘fat tailed’. This means that, in reality, there are more stock market ‘crashes’ and ‘rallies’ than a strictly normal model of stock returns would suggest. Shortly, I will show how the classical model can be extended to explicitly account for fat-tailed returns.

Here is the QuantLib C++ code to generate a single stock price path under the assumption of normally distributed asset returns:

```
#include <ql/quantlib.hpp>
using namespace QuantLib;
BOOST_AUTO_TEST_CASE(testGeometricBrownianMotion) {
Real startingPrice = 20.16; //closing price for INTC on 12/7/2012
Real mu = .2312; //INTC one year historical annual return
Volatility sigma = 0.2116; //INTC one year historical volatility
Size timeSteps = 255; //trading days in a year (U.S.)
Time length = 1; //one year
//instantiate Geometric Brownian Motion (GBM) stochastic process
const boost::shared_ptr<StochasticProcess>& gbm =
boost::shared_ptr<StochasticProcess> (new GeometricBrownianMotionProcess(startingPrice,
mu, sigma));
//generate a sequence of normally distributed random numbers from a
//uniform distribution using Box-Muller transformation
BigInteger seed = SeedGenerator::instance().get();
typedef BoxMullerGaussianRng MersenneBoxMuller;
MersenneTwisterUniformRng mersenneRng(seed);
MersenneBoxMuller boxMullerRng(mersenneRng);
RandomSequenceGenerator<MersenneBoxMuller> gsg(timeSteps, boxMullerRng);
//generate simulated path of stock price using GBM stochastic process
PathGenerator<RandomSequenceGenerator<MersenneBoxMuller> > gbmPathGenerator(gbm, length,
timeSteps, gsg, false);
const Path& samplePath = gbmPathGenerator.next().value;
//calculate simulated sample returns using C++11 lambda expression
boost::function<Real, (Real, Real)> calcLogReturns = [](Real x, Real y) {return std::log(y/x);};
std::vector<Real> logReturns;
Path::iterator samplePathBegin = samplePath.begin();
Path::iterator samplePathEnd = samplePath.end();
Path::iterator endMinusOne = std::prev(samplePathEnd);
Path::iterator beginPlusOne = std::next(samplePathBegin);
std::transform(samplePathBegin, endMinusOne, beginPlusOne,
std::back_inserter(logReturns), calcLogReturns);
//calculate some general statistics
GeneralStatistics statistics;
//returns statistics
statistics.addSequence(logReturns.begin(), logReturns.end());
std::cout << boost::format("Standard deviation of simulated returns (Normal):
%.4f") % (statistics.standardDeviation() * std::sqrt(255)) << std::endl;
//price statistics
statistics.reset();
statistics.addSequence(samplePath.begin(), samplePath.end());
std::cout << boost::format("Price statistics: mean=%.2f, min=%.2f, max=%.2f") %
statistics.mean() % statistics.min() % statistics.max() << std::endl;
//write simulated path to a file for charting with gnuplot
std::ofstream gbmFile;
gbmFile.open("/tmp/gbm.dat",
std::ios::out);
for (Size i = 0; i < timeSteps; ++i) {
gbmFile << boost::format("%d %.4f") % i % samplePath.at(i) << std::endl;
}
gbmFile.close();
/* gnuplot script to chart stock price path
set key bottom center
set key bottom box
set xlabel "Time Step (Days)"
set ylabel "Stock Price"
plot "/tmp/gbm.dat" using 1:2 w lines t "INTC simulated stock price"
*/
}
}
```

When executed, this code will produce output similar to the following, though the statistics will vary with each run of the program:

Standard deviation of simulated returns (Normal): 0.2105

Price statistics: mean=24.99, min=19.87, max=29.02

The resulting simulated stock price data, when charted, will look something like this:

A few comments:

The code produces a time series of 255 daily closing prices (one for each trading day in the U.S. equity markets), using QuantLib’s PathGenerator class, which is a C++ template parameterized by a RandomSequenceGenerator. The RandomSequenceGenerator requires a source of randomness, which is provided by way of QuantLib’s MersenneTwisterUniformRng class. The PathGenerator also requires a StochasticProcess as its first constructor argument to which a (shared) pointer to a GeometricBrownianMotionProcess instance is passed. The starting price for the process is the closing price of Intel as of 12/7/2012. The *mu* parameter is bound to INTC’s annualized return, while the *sigma* parameter is equated to INTC’s annualized historical volatility for the period from 12/7/2012 through 12/7/2013.

The Box-Muller transformation is employed to convert a sequence of uniformly distributed random numbers to a sequence of normally distributed random numbers. You can read more about the Box-Muller transformation here.

QuantLib’s GeneralStatistics class is applied to collect some process statistics in order to measure how accurately the simulated asset path approximates the INTC time series. By design, the simulated asset path’s annualized return and volatility should resemble the *mu* and *sigma* process parameters, respectively. For some applications, such as back-testing path-dependent trading strategies or evaluating certain risk management scenarios, we might also like the simulated minimum, maximum and mean stock prices to fall within a desired range.

Towards the end of the listing, the data points of the simulated asset price path are written to a file so that they can be charted with gnuplot, an open-source charting package. The gnuplot script is included at the end of the source code listing as a comment.

### Leptokurtic Model of Asset Returns

Now let’s revisit a concept I mentioned earlier regarding the distribution of asset returns. It’s well known that asset returns are ‘fat-tailed’, which is an intuitive way of stating that daily returns are observed in practice to fall more often in the lower and upper tails of their frequency distribution than would be predicted by the normal probability distribution function (pdf). As such, asset returns can be characterized as *leptokurtic*, meaning that their frequency distribution exhibits *excess* *kurtosis *compared to the normal distribution, which has a kurtosis value of 3.

How can we improve our model of asset prices to incorporate the leptokurtic nature of asset returns that we see in the real-world? Essentially, we would like the simulated asset price paths to exhibit more days with higher volatility compared to a price path generated in accordance with the classical, normally distributed model of asset returns, while realizing the same annualized volatility as the normal model. To achieve this we need to sample our random numbers from a non-normal distribution with excess kurtosis. One such distribution often used for this purpose is the Student’s t distribution, which belongs to a family of stable distributions known as Cauchy distributions, which are symmetrical about the mean but with ‘fat tails’. The shape of the Student’s t distribution is determined by a single parameter,* v*, which is defined as *degrees of freedom. *The greater the number of degrees of freedom, the more the Student’s t distribution resembles the normal distribution, the smaller the ‘fatter’ the tails.

So let’s make the necessary changes to the code above to generate random numbers distributed in accordance with the Student’s t distribution. This is accomplished with QuantLib’s InverseCumulativeRsg class, which is a template parameterized by a RandomSequenceGenerator and an inverse cumulative distribution function. The implementation of the Student’s t distribution is provided by the Boost Math library class boost::math::students_t_distribution, which is instantiated with five degrees of freedom, so the shape of the distribution has the desired ‘fat tails’, as illustrated by the light blue density in the image below. Note that there is more mass in the tails of the light blue density than the black, standard normal density.

The revised source code listing is as follows:

```
#include <ql/quantlib.hpp>
#include <boost/math/distributions/students_t.hpp>
Real studentTInverse(boost::math::students_t_distribution d, const Real& p) {
return quantile(d,p)
}
using namespace QuantLib;
BOOST_AUTO_TEST_CASE(testGeometricBrownianMotionStudentT) {
Real startingPrice = 20.16; //closing price for INTC on 12/7/2012
Real mu = .2312; //INTC one year historical annual return
Volatility sigma = 0.2116; //INTC one year historical volatility
Volatility scaledSigma = std::sqrt(sigma * sigma * 3/5); //scaled by reciprocal of Student T variance (v/(v-2))
Size timeSteps = 255; //trading days in a year (U.S.)
Time length = 1; //one year
//instantiate Geometric Brownian Motion (GBM) stochastic process
const boost::shared_ptr<StochasticProcess>& gbm =
boost::shared_ptr<StochasticProcess> (new GeometricBrownianMotionProcess(startingPrice,
mu, scaledSigma));
//random sequence generator uses Mersenne Twiseter to generate uniformly distributed
//pseudo-random numbers
BigInteger seed = SeedGenerator::instance().get();
MersenneTwisterUniformRng mersenneRng(seed);
RandomSequenceGenerator<MersenneTwisterUniformRng> rsg(timeSteps, mersenneRng);
//instantiate Student T distribution from Boost math library
boost::math::students_t_distribution<> studentT(5); //5 degrees of freedom - want fat tails!
boost::function<Real (Real)> icd = boost::bind(studentTInverse, studentT, _1);
//sample random numbers from the Student T distribution
InverseCumulativeRsg<RandomSequenceGenerator<MersenneTwisterUniformRng>,
boost::function<Real (Real)> > invCumRsg(rsg, icd);
//generates a single path
PathGenerator<InverseCumulativeRsg<RandomSequenceGenerator<MersenneTwisterUniformRng>,
boost::function<Real (Real)> > > gbmPathGenerator(gbm, length, timeSteps, invCumRsg, false);
const Path& samplePath = gbmPathGenerator.next().value;
//calculate simulated sample returns using C++11 lambda expression
boost::function<Real, (Real, Real)> calcLogReturns = [](Real x, Real y) {return std::log(y/x);};
std::vector<Real> logReturns;
Path::iterator samplePathBegin = samplePath.begin();
Path::iterator samplePathEnd = samplePath.end();
Path::iterator endMinusOne = std::prev(samplePathEnd);
Path::iterator beginPlusOne = std::next(samplePathBegin);
std::transform(samplePathBegin, endMinusOne, beginPlusOne,
std::back_inserter(logReturns), calcLogReturns);
//calculate some general statistics
GeneralStatistics statistics;
//returns statistics
statistics.addSequence(logReturns.begin(), logReturns.end());
std::cout << boost::format("Standard deviation of simulated returns (Student-T):
%.4f") % (statistics.standardDeviation() * std::sqrt(255)) << std::endl;
//price statistics
statistics.reset();
statistics.addSequence(samplePath.begin(), samplePath.end());
std::cout << boost::format("Price statistics: mean=%.2f, min=%.2f, max=%.2f") %
statistics.mean() % statistics.min() % statistics.max() << std::endl;
//write simulated path to a file for charting with gnuplot
std::ofstream gbmFile;
gbmFile.open("/tmp/gbm-student.dat",
std::ios::out);
for (Size i = 0; i < timeSteps; ++i) {
gbmFile << boost::format("%d %.4f") % i % samplePath.at(i) << std::endl;
}
gbmFile.close();
/* gnuplot script to chart stock price path
set key bottom center
set key bottom box
set xlabel "Time Step (Days)"
set ylabel "Stock Price"
plot "/tmp/gbm.dat" using 1:2 w lines
t "INTC Normal", "/tmp/gbm-student.dat" using 1:2 w lines t "INTC Student-T"
*/
}
}
```

The test case, when run, will produce output such as the following, though again, the statistics will vary from run to run of the program:

Standard deviation of simulated returns (Student-T): 0.2112

Price statistics: mean=24.94, min=19.44, max=30.24

As depicted in the chart below, the `testGeometricBrownianMotionStudentT`

test case generates a time series that features more relatively large up and down moves than that of the classic, normal model. However, the realized volatility of the Student’s t time series is not significantly different than the normal:

There is much more to the modeling of asset price dynamics than I can present in this brief post, but I hope that I’ve imparted the fundamental concepts and explained how you can accomplish quite a lot with very little effort using QuantLib.

Also, for your convenience, I have now made all of the source code from my Introducing QuantLib series available on GitHub at https://github.com/mhittesdorf/allthingsfintech. I encourage you to download and experiment with the code yourself.

Thanks for reading my blog and, as always, please feel free to leave your comments and questions below. Have fun with QuantLib!

How does

Real studentTInverse(boost::math::students_t_distribution d, const Real& p) {

return quantile(d,p)

}

compile when you have “using namespace QuantLib” after you’ve elected to return a Real type.

Mick, your posts are very well detailed and helpful. Thanks for putting this up.

Very helpful, thanks.

How to compile this under Ubuntu?

All libs are installed, but compilation fails.

A Makefile would be appreciated.

All the source code and cmake build instructions can be found at https://github.com/mhittesdorf/allthingsfintech.

Is it possible to feed this back into one of QL’s option pricing models?

Conceptually, yes. The underlying value and volatility of the generated time series at a given point in time could be used to calculate the value of an option at that moment in time.

Really helpful explanation and code. I’m playing with it, but I’m stuck when I want to create sample paths of different number of days. What if I want to create a sample path for 100 days or for 500 days? Which variables should I change then? Should I only change the timeSteps variable or also the length variable?