Python: Path Generator for Correlated Processes

One reader was interested to know, how to simulate correlated asset paths by using just Python libraries, without using QuantLib. This blog post is presenting the result of woodshedding this stuff. A couple of notes:

• GeneratePaths method can be used to simulate paths for a single process or multiple processes, based on a given (or non-existing) correlation matrix.
• Method returns Numpy array having dimensions based on the given number of processes, number of paths and number of time steps.
• Discretized stochastic processes for generator method to be simulated are defined as lambda methods. This approach makes it relatively easy to implement (almost) any desired single-factor model.

The program is enough commented and should be self-explainable. Thanks for reading this blog.
-Mike

import numpy as np
import matplotlib.pyplot as pl

# returns ndarray with the following dimensions: nProcesses, nPaths, nSteps
def GeneratePaths(spot, process, maturity, nSteps, nPaths, correlation = None):
    dt = maturity / nSteps
    
    # case: given correlation matrix, create paths for multiple correlated processes
    if (isinstance(correlation, np.ndarray)):
        nProcesses = process.shape[0]
        result = np.zeros(shape = (nProcesses, nPaths, nSteps))
        
        # loop through number of paths
        for i in range(nPaths):
            # create one set of correlated random variates for n processes
            choleskyMatrix = np.linalg.cholesky(correlation)
            e = np.random.normal(size = (nProcesses, nSteps))            
            paths = np.dot(choleskyMatrix, e)
            # loop through number of steps
            for j in range(nSteps):
                # loop through number of processes
                for k in range(nProcesses):
                    # first path value is always current spot price
                    if(j == 0):
                        result[k, i, j] = paths[k, j] = spot[k]
                    else:
                        # use SDE lambdas (inputs: previous spot, dt, current random variate)
                        result[k, i, j] = paths[k, j] = process[k](paths[k, j - 1], dt, paths[k, j])

    # case: no given correlation matrix, create paths for a single process
    else:
        result = np.zeros(shape = (1, nPaths, nSteps))
        # loop through number of paths
        for i in range(nPaths):
            # create one set of random variates for one process
            path = np.random.normal(size = nSteps)
            # first path value is always current spot price
            result[0, i, 0] = path[0] = spot
            # loop through number of steps
            for j in range(nSteps):
                if(j > 0):
                    # use SDE lambda (inputs: previous spot, dt, current random variate)
                    result[0, i, j] = path[j] = process(path[j - 1], dt, path[j])
    return result

# Geometric Brownian Motion parameters
r = 0.03
v = 0.25

# define lambda for process (inputs: spot, dt, random variate)
BrownianMotion = lambda s, dt, e: s + r * s * dt + v * s * np.sqrt(dt) * e   

# general simulation-related parameters
maturity = 1.0
nPaths = 10
nSteps = 250

# case: one process
SingleAssetPaths = GeneratePaths(100.0, BrownianMotion, maturity, nSteps, nPaths)
for i in range(nPaths):
    pl.plot(SingleAssetPaths[0, i, :])
pl.show()

# case: two correlated processes
matrix = np.array([[1.0, 0.999999], [0.999999, 1.0]])
spots = np.array([100.0, 100.0])
processes = np.array([BrownianMotion, BrownianMotion])
MultiAssetPaths = GeneratePaths(spots, processes, maturity, nSteps, nPaths, matrix)
f, subPlots = pl.subplots(processes.shape[0], sharex = True)
for i in range(processes.shape[0]): 
    for j in range(nPaths):
        subPlots[i].plot(MultiAssetPaths[i, j, :])
pl.show()

C++/Python: Creating Python Wrapper for C++ Class by Using SWIG

If the need sometimes arises, existing C++ libraries can be interfaced relatively easy to be used in Python by using SWIG wrapper. SWIG stands for Simplified Wrapper and Interface Generator and it is an open-source software tool used to connect computer programs or libraries written in C or C++ with scripting languages, such as Python. As one extremely motivating example, QuantLib C++ libraries have been made available for Python by using SWIG wrappers. Complete documentation for SWIG 3.0 can be found in here. The actual SWIG package can be downloaded in here. By taking a look at the documentation, one may understand quickly, that in this post only the tip of the iceberg will get scratched.

Also, it should be noted at this point, that this post is assuming all operations taking place in Linux environment. Corresponding instructions for Windows users can be found in SWIG documentation. This blog post is presenting, how to interface custom C++ random generator class to be used in Python. The original (template) class implementation can be found in here.

Header and implementation files for simple C++ random number generator class are presented below.

Header file (RG.h)

#include <algorithm>
#include <functional>
#include <random>

class Generator {
public:
 Generator(unsigned long seed);
 void operator()(std::vector<double>& v);
private:
 std::function<double(double)> randomGenerator;
 std::normal_distribution<double> distribution;
 std::mt19937 uniformGenerator;
};

Implementation file (RG.cpp)

#include "RG.h"

Generator::Generator(unsigned long seed) {
  // construct lambda method for processing standard normal random number
  randomGenerator = [this](double x)-> double {
   x = distribution(uniformGenerator);
   return x;
  };
  uniformGenerator.seed(seed);
}

// fill client-given vector with random variates from standard normal distribtuion
void Generator::operator()(std::vector<double>& v) {
  std::transform(v.begin(), v.end(), v.begin(), randomGenerator);
}

The class functionality is simple and straightforward. First, create a class implementation by giving desired seed value in constructor. In constructor, C++ function object will be created. This function object will ultimately create our normally distributed random variates. Client will request a new set of random variates by giving a reference for STL vector by using parenthesis operator, which has been overloaded here for syntactic reasons only.

SWIG interface file (RG.i)

A simple SWIG interface for our class can be built by simply wrapping the header file as follows.

%module RG
%{
#include "RG.h"
%}

%include "std_vector.i"
%template(DoubleVector) std::vector<double>;

%include "RG.h"

Python does not know anything about std::vector. The std_vector.i library provides support for C++ STL vector class. Then, all we need to do is to instantiate different versions of vector for all the specific types we would like to use in Python. For our example class, we want to use vector consisting of double data types. Using this library involves the use of the %template directive. In Python code, C++ vector will be used by using alias DoubleVector.

SWIG build

All in all we have now three files in our folder: RG.h, RG.cpp and RG.i. With terminal opened in the folder consisting the previous three files, the following three commands will create Python wrapper for our C++ class.

$ swig -c++ -python RG.i
$ g++ -fpic -c -I/home/mikejuniperhill/anaconda3/include/python3.7m RG.cpp RG_wrap.cxx
$ g++ -shared RG.o RG_wrap.o -o _RG.so

In the second command, a correct parameter for -I argument can be found by using the following command in terminal.

$ python3-config --cflags

After these three steps, our C++ program can be accessed and used in Python program. The following simple example program is requesting a vector of standard normal random numbers from our wrapped C++ class and uses those for calculating Monte Carlo PV for one European call option.

# import RG module
import SwigTester.RG as RG
import numpy as np

# set seed, create generator object
seed = 0
generator = RG.Generator(seed)

# set vector size, initialize vector
steps = 1000
vector = RG.DoubleVector(steps)

# request generator to fill vector with standard normal variates 
generator(vector)
# copy content to numpy array
e = np.array(vector)

# use variates to value one european call option
# set parameters
s0 = 100.0
x = 100.0
t = 1.25
v = 0.25
r = 0.01

# use vectorization in numpy, pre-calculate components
drift = (r - 0.5 * v * v) * t
diffusion = v * np.sqrt(t)
df = np.exp(-r * t)

# get option value as discounted average of all terminal payoffs
c0 = np.mean(np.maximum(s0 * np.exp(drift + diffusion * e) - x, 0.0)) * df
print(c0)

Note, that I imported SwigTester.RG, because I created Python wrapper module directly into this specific folder.

As a SWIG newcomer, one may expect to face all kinds of cryptic installation and other infuriating compatibility issues. As always, Google is the best friend along this cruel and unusual phase. But, have faith - and you'll be rewarded with all wonders and joys of SWIG wrapper.

All relevant files can be found directly from my GitHub repository. Thanks for reading my blog.
-Mike

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