C++11 : modelling one-factor processes using functional programming paradigm

There was an interesting technical article on July 2017 Wilmott magazine written by Daniel Duffy and Avi Palley. This multi-page article was giving an overview on some "game changing" Boost libraries, which have been accepted as a part of C++11/14 standard, such as smart pointers, function wrappers, lambda expressions and tuples. The second part of this article series will be published in the next Wilmott magazine and it will present (according to editor) design and implementation for Monte Carlo option pricing framework. It is also an important point of mentioning, that (according to Amazon.com) Daniel Duffy is about to publish long-awaited second edition of his book on pricing derivatives using C++. The book was initially published somewhere back in 2004 and the landscape has changed quite dramatically since these days.

Within the last chapter of this article, functional programming paradigm was nicely applied for modelling one-factor stochastic differential equations, generally used in Finance. By applying more functional programming paradigm, the usually observed code bloating can be substantially reduced. As an example of such code bloat, I reviewed my own implementation for path generator, which models one-factor processes for Monte Carlo purposes. There is an abstract base class (OneFactorProcess) and implementations for GBM and Vasicek processes. Even there is nothing fundamentally wrong with this approach (class hierarchy), one may ask, whether there would be a bit more flexible ways to implement this kind of a scheme.

Within this post, I have been re-designing modelling part for one-factor processes by applying functional programming paradigm, as presented in that article. Reduction in code bloating is present, since there is currently only one single class for modelling different types of processes (before, there was a class hierarchy). Moreover, since the both functions for handling drift and diffusion terms will be constructed outside of this class, their construction process is now much more flexible than before.

The program

 

Implement the following program (two header files and one implementation file for tester) into a new project. For brevity reasons, I have re-designed only the part of the program, which models one-factor processes. Monte Carlo part has been implemented as free function in tester implementation file. First, one-factor process object (Process) will be created by using ProcessBuilder object (Builder Pattern). Within this example, I have implemented a builder for constructing Process object by using console. However, the flexibility in this program allows different types of builders to be implemented. As soon as Process object is created, "skeleton path" (std::vector) will be sent to Monte Carlo method (MonteCarloLite), along with all simulation-related attributes and Process object. As a result, this method will fill vector with the values from one simulation path for chosen parameters and applied stochastic process. Finally, a path will be printed back to console.

#pragma once
#include <functional>
// OneFactorSDE.h
namespace FunctionalOneFactorProcessExampleNamespace
{
 // blueprint for a function, which models drift or diffusion
 using Function = std::function<double(double s, double t)>;
 // wrapper for drift and diffusion function components
 using Functions = std::tuple<Function, Function>;
 //
 // class for modeling one-factor processes by employing drift and diffusion functions
 class Process
 {
 public:
  Process(Functions functions, double initialCondition)
   : driftFunction(std::get<0>(functions)), diffusionFunction(std::get<1>(functions)),
   initialCondition(initialCondition) { }
  double drift(double s, double t) { return this->driftFunction(s, t); }
  double diffusion(double s, double t) { return this->diffusionFunction(s, t); }
  double InitialCondition() const { return this->initialCondition; }
  //
 private:
  double initialCondition; // spot price for a process
  Function driftFunction; // function for modelling process drift
  Function diffusionFunction;  // function for modelling process diffusion
 };
}
//
//
//
// ProcessFactory.h
#pragma once
#include <iostream>
#include <string>
#include <memory>
#include "OneFactorSDE.h"
//
namespace FunctionalOneFactorProcessExampleNamespace
{
 // abstract base class for all process builders
 class ProcessBuilder
 {
 public:
  // return process object which is wrapped inside shared pointer
  // let implementation classes decide, how and from where to built process object
  virtual std::shared_ptr<Process> Build() = 0;
 };
 //
 // specific implementation for console process builder
 // process type and corresponding parameters will be requested from a client by using console
 class ConsoleProcessBuilder : public ProcessBuilder
 {
 public:
  std::shared_ptr<Process> Build() override
  {
   Functions functions;
   double initialCondition = 0.0;
   std::string consoleSelection;
   std::cout << "Select process [1 = GBM, 2 = Vasicek] > ";
   // if conversion cannot be performed, stoi will throw invalid argument exception
   std::getline(std::cin, consoleSelection);
   int processID = std::stoi(consoleSelection);
   //
   switch (processID)
   {
   // GBM process
   case 1:
   {
    // receive client inputs
    std::cout << "spot price > ";
    std::getline(std::cin, consoleSelection);
    initialCondition = std::stod(consoleSelection);
    //
    std::cout << "risk-free rate > ";
    std::getline(std::cin, consoleSelection);
    double r = std::stod(consoleSelection);
    //
    std::cout << "volatility > ";
    std::getline(std::cin, consoleSelection);
    double v = std::stod(consoleSelection);
    //
    // build drift and diffusion functions for GBM process object
    auto driftFunction = [r](double S, double T){ return r * S; };
    auto diffusionFunction = [v](double S, double T){ return v * S; };
    // wrap drift and diffusion functions into tuple
    functions = std::make_tuple(driftFunction, diffusionFunction);
    break;
   }
   case 2:
   {
    // receive client inputs
    std::cout << "spot price > ";
    std::getline(std::cin, consoleSelection);
    initialCondition = std::stod(consoleSelection);
    //
    std::cout << "reversion > ";
    std::getline(std::cin, consoleSelection);
    double reversion = std::stod(consoleSelection);
    //
    std::cout << "long-term rate > ";
    std::getline(std::cin, consoleSelection);
    double longTermRate = std::stod(consoleSelection);
    //
    std::cout << "rate volatility > ";
    std::getline(std::cin, consoleSelection);
    double v = std::stod(consoleSelection);
    //
    // build drift and diffusion functions for Vasicek process object
    auto driftFunction = [reversion, longTermRate](double S, double T)
     { return reversion * (longTermRate - S); };
    auto diffusionFunction = [v](double S, double T){ return v; };
    // wrap drift and diffusion functions into tuple
    functions = std::make_tuple(driftFunction, diffusionFunction);
    break;
   }
   default:
    // if selected process is not configured, program will throw invalid argument exception
    throw std::invalid_argument("invalid process ID");
    break;
   }
   // build and return constructed process object for a client
   // wrapped into shared pointer
   return std::shared_ptr<Process>(new Process(functions, initialCondition));
  }
 };
}
//
//
//
// Tester.cpp
#include <vector>
#include <random>
#include <algorithm>
#include <chrono>
#include "ProcessFactory.h"
namespace MJ = FunctionalOneFactorProcessExampleNamespace;
//
void MonteCarloLite(std::vector<double>& path, 
 std::shared_ptr<MJ::Process>& process, const double maturity);
//
int main()
{
 try
 {
  // create process object by using console builder
  MJ::ConsoleProcessBuilder builder;
  std::shared_ptr<MJ::Process> process = builder.Build();
  //
  // create simulation-related attributes
  int nSteps = 100;
  double timeToMaturity = 1.25;
  // create, process and print one path
  std::vector<double> path(nSteps);
  MonteCarloLite(path, process, timeToMaturity);
  std::for_each(path.begin(), path.end(), 
   [](double v) { std::cout << v << std::endl; });
 }
 catch (std::exception e)
 {
  std::cout << e.what() << std::endl;
 }
 return 0;
}
//
void MonteCarloLite(std::vector<double>& path, 
 std::shared_ptr<MJ::Process>& process, const double maturity)
{
 // lambda method for seeding uniform random generator
 std::function<unsigned long(void)> seeder =
  [](void) -> unsigned long { return static_cast<unsigned long>
  (std::chrono::steady_clock::now().time_since_epoch().count()); };
 //
 // create uniform generator, distribution and random generator function
 std::mt19937 uniformGenerator(seeder());
 std::normal_distribution<double> distribution;
 std::function<double(double)> randomGenerator;
 //
 // lambda method for processing standard normal random numbers
 randomGenerator = [uniformGenerator, distribution](double x) mutable -> double
 {
  x = distribution(uniformGenerator);
  return x;
 };
 //
 // create vector of standard normal random numbers
 // use lambda method created earlier
 std::transform(path.begin(), path.end(), path.begin(), randomGenerator);
 //
 double dt = maturity / (path.size() - 1);
 double dw = 0.0;
 double s = (*process).InitialCondition();
 double t = 0.0;
 path[0] = s; // 1st path element is always the current spot price
 //
 // transform random number vector into a path containing asset prices
 for (auto it = path.begin() + 1; it != path.end(); ++it)
 {
  t += dt;
  dw = (*it) * std::sqrt(dt);
  (*it) = s + (*process).drift(s, t) * dt + (*process).diffusion(s, t) * dw;
  s = (*it);
 }
}

As always, thanks a lot for reading this blog.
-Mike