FINAL DESTINATION
The circle has been closed. This post is kind of an aggregation, based on the last four posts published on generating random numbers.
Initially, I started just with a simple template class for distributional random generator,
then continued with a path generator using any one-factor stochastic process and
finally, ended up with a multi-threaded distributional random generation scheme using
Parallel algorithms. This final post (hopefully) is opening up my larger goal :
to be able to generate asset price paths for any one-factor process, using multi-threading scheme.
GROUNDHOG DAY
Again, I have tested the both sequential (for_each) and parallel (parallel_for_each)
schemes by using four generators, 10000 paths and 250 time steps for a
single run. After this, I repeated this run for 250 times. Conclusion :
- The average
running time for this sample was 17116 milliseconds (sequential) and 8209
milliseconds (parallel). So, parallel scheme will be completed about
two times faster.
- The actual CPU usage profiles during the simulation processes are behaving exactly as reported in this post.
- I also analyzed processed asset price paths for parallel scheme, just to
be absolutely sure there are no path duplicates (random number
generation would not be independent). Based on my analysis made in
Excel, all processed asset price paths are different and there
are no duplicates.
Presented scheme for path generator is again fulfilling my two
initial requirements : faster creation of asset price paths following any one-factor process and independency of random generators.
RandomGenerator.h
The basic functionality of this template class has not been changed, except for construction part : second constructor is allowing a client to give any probability distribution for uniform generator from outside of this class. Even there is actually no need for having this kind of optionality in real-life (most of the stuff in Monte Carlo is randomized by using standard normal distribution), I decided to implement this optionality for the sake of completeness.
#pragma once
#include <algorithm>
#include <functional>
#include <vector>
#include <random>
#include <memory>
//
namespace MikeJuniperhillRandomGeneratorTemplate
{
template <typename Generator = std::mt19937, typename Distribution = std::normal_distribution<double>>
/// <summary>
/// Template class for creating random number paths using mt19937 as default uniform
/// random generator and Standard Normal as default probability distribution.
/// </summary>
class RandomGenerator
{
public:
/// <summary>
/// Constructor with explicit seed value
/// </summary>
RandomGenerator(unsigned long seed)
{
// construct function for processing distributional random number
randomGenerator = [this](double x)-> double
{
x = distribution(uniformGenerator);
return x;
};
// seed generator once
uniformGenerator.seed(seed);
}
/// <summary>
/// Constructor for explicit seed value and client-given probability distribution.
/// </summary>
RandomGenerator(unsigned long seed, const Distribution& distribution)
// constructor delegation
: RandomGenerator(seed)
{
// assign client-given probability distribution
this->distribution = distribution;
}
/// <summary>
/// Fill a given vector reference with distributional random numbers
/// </summary>
void operator()(std::vector<double>& v) const
{
std::transform(v.begin(), v.end(), v.begin(), randomGenerator);
}
private:
std::function<double(double)> randomGenerator;
Generator uniformGenerator;
Distribution distribution;
};
}
//
OneFactorProcess.h
I decided to tag drift and diffusion functions with
const declaration, since these functions should not modify the internal state of class data members.
#pragma once
//
namespace MikeJuniperhillOneFactorProcessLibrary
{
/// <summary>
/// Abstract base class for all one-factor processes for customizing
/// drift and diffusion functions for different types of processes.
/// </summary>
class OneFactorProcess
{
public:
virtual double drift(double x, double t) const = 0;
virtual double diffusion(double x, double t) const = 0;
};
//
/// <summary>
/// Implementation for Vasicek short-rate model.
/// </summary>
class Vasicek : public OneFactorProcess
{
public:
Vasicek(double meanReversion, double longTermRate, double rateVolatility)
: meanReversion(meanReversion), longTermRate(longTermRate), rateVolatility(rateVolatility) { }
//
double drift(double x, double t) const override { return meanReversion * (longTermRate - x); }
double diffusion(double x, double t) const override { return rateVolatility; }
private:
double meanReversion;
double longTermRate;
double rateVolatility;
};
}
//
PathGenerator.h
As in the case with RandomGenerator, the basic functionality of this template class has not been changed either,
except for construction part : second constructor is allowing a client
to give any probability distribution to be delivered for distributional random generator.
#pragma once
//
#include "RandomGenerator.h"
#include "OneFactorProcess.h"
namespace MJRandom = MikeJuniperhillRandomGeneratorTemplate;
namespace MJProcess = MikeJuniperhillOneFactorProcessLibrary;
//
namespace MikeJuniperhillPathGenerator
{
template <typename Generator = std::mt19937, typename Distribution = std::normal_distribution<double>>
class PathGenerator
{
public:
/// <summary>
/// Constructor for PathGenerator template class.
/// </summary>
PathGenerator(double spot, double maturity, unsigned long seed,
const std::shared_ptr<MJProcess::OneFactorProcess>& process)
: spot(spot), maturity(maturity), process(process)
{
// create random generator
generator = std::unique_ptr<MJRandom::RandomGenerator<Generator, Distribution>>
(new MJRandom::RandomGenerator<Generator, Distribution>(seed));
}
/// <summary>
/// Constructor for PathGenerator template class, with a client-given probability distribution
/// </summary>
PathGenerator(double spot, double maturity, unsigned long seed,
const std::shared_ptr<MJProcess::OneFactorProcess>& process, const Distribution& distribution)
: spot(spot), maturity(maturity), process(process)
{
// create random generator with client-given probability distribution
generator = std::unique_ptr<MJRandom::RandomGenerator<Generator, Distribution>>
(new MJRandom::RandomGenerator<Generator, Distribution>(seed, distribution));
}
/// <summary>
/// Fill a given vector reference with asset prices, following a given stochastic process.
/// </summary>
void operator()(std::vector<double>& v) const
{
// transform initialized vector into a path containing random numbers
(*generator)(v);
//
double dt = maturity / (v.size() - 1);
double dw = 0.0;
double s = spot;
double t = 0.0;
v[0] = s; // 1st path element will always be the current spot price
//
// transform distributional random number vector into a path containing
// asset prices from a given stochastic one-factor process
for (auto it = v.begin() + 1; it != v.end(); ++it)
{
t += dt;
dw = (*it) * std::sqrt(dt);
(*it) = s + (*process).drift(s, t) * dt + (*process).diffusion(s, t) * dw;
s = (*it);
}
}
private:
double spot;
double maturity;
std::shared_ptr<MJProcess::OneFactorProcess> process;
std::unique_ptr<MJRandom::RandomGenerator<Generator, Distribution>> generator;
};
}
//
Tester.cpp
Tester program is closely tracking the program presented in
previous post. For the sake of additional clarity, I have used new type definitions in order to improve code readability and get rid of some lengthy variable names. The program is again using simple factory method for creating PathGenerator (function wrapped in shared pointer). In this program, OneFactorProcess implementation is created and delivered for factory method for processing. Finally, there is a method for printing processed paths to console for testing purposes.
#include <iostream>
#include <chrono>
#include <ppl.h>
#include <concurrent_vector.h>
#include "PathGenerator.h"
namespace MJGenerator = MikeJuniperhillPathGenerator;
//
// type definitions
using Path = std::vector<double>;
using Paths = concurrency::concurrent_vector<Path>;
using Process = std::shared_ptr<MJProcess::OneFactorProcess>;
using Processor = std::function<void(void)>;
using PathGenerator = std::shared_ptr<Processor>;
//
// thread-safe container for storing asset price paths, processed by path generators
Paths paths;
//
// printer for generated asset price paths
void Printer()
{
std::for_each(paths.begin(), paths.end(),
[](Path path) -> void
{
std::for_each(path.begin(), path.end(),
[](double s) -> void
{
std::cout << s << ",";
});
std::cout << std::endl;
});
}
//
// factory method :
// return path-generating function as function wrapper
// input arguments are common for all types of generators
PathGenerator Factory(double spot, double maturity, int nPaths,
int nSteps, unsigned long seed, const Process& process, Paths& paths)
{
// create function for processing one-factor paths
auto generator = [=, &process, &paths]() -> void
{
MJGenerator::PathGenerator<> oneFactorProcess(spot, maturity, seed, process);
Path path(nSteps);
for (auto i = 0; i != nPaths; ++i)
{
oneFactorProcess(path);
paths.push_back(path);
}
};
// return generator function as function wrapper
return PathGenerator(new Processor(generator));
}
//
int main()
{
// create vasicek process
double longTermRate = 0.05;
double meanReversion = 0.2;
double rateVolatility = 0.0075;
Process vasicek = Process(new MJProcess::Vasicek(meanReversion, longTermRate, rateVolatility));
//
// define parameters and seed values for path generators
int nGenerators = 4;
int nPaths = 100;
int nSteps = (250 + 1);
std::vector<unsigned long> seed = { 10322854, 65947, 387528, 772399573 };
//
// use factory method for creating path generators
double spot = 0.0095;
double maturity = 3.0;
std::vector<PathGenerator> generators;
for (auto i = 0; i < nGenerators; i++) generators.push_back(
Factory(spot, maturity, nPaths, nSteps, seed[i], vasicek, paths));
//
// parallel processing
auto start = std::chrono::steady_clock::now();
concurrency::parallel_for_each(generators.begin(), generators.end(),
[](PathGenerator pg) -> void { (*pg)(); });
auto end = std::chrono::steady_clock::now();
auto timeElapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
std::cout << timeElapsed.count() << std::endl;
//
// print paths
Printer();
return 0;
}
//
Finally, thanks again for reading this blog.
-Mike
