QuantLib-Python: Heston Monte Carlo Valuation for Autocallable Memory Coupon Note

In the spirit of the previous post, I was woodshedding an implementation for valuing Autocallable Memory Coupon note by using libraries available in QuantLib-Python. These products are embedding a series of out-of-the-money barrier options and for this specific reason, it is important to capture implied volatility smile by using appropriate model. For this implementation example, Heston stochastic volatility model has been used. In addition to the actual Monte Carlo algorithm and path generator, I also implemented a simple method for calibrating Heston model to volatility surface by using SciPy optimization package. The complete program can be downloaded from my GitHub page.

Autocallable Notes

Autocallable notes are path-dependent structured products, which may be linked to any index. However, this product is usually linked to equity index or equity basket. If condition called autocall barrier will be reached on any coupon date, the product will be terminated immediately, after paying back coupon (usually fixed) and notional amount. However, if the underlying index will not reach autocall barrier on a coupon date, but reaches condition called coupon barrier, only fixed coupon will be paid on notional and the product will be kept alive until maturity (or until the next autocall event). Autocallables may also have a memory, which means that if the underlying index was below coupon barrier in the previous coupon date (and there was no coupon payment), but will be above coupon barrier in the next coupon date, the product will pay all cumulated (and unpaid) coupons on notional in the next coupon date. Autocallable products are (usually) not capital protected, which means that the investor is not guaranteed to receive the initial investment. If the index fixing on a final coupon date will be below condition called protection barrier, the final redemption amount will be calculated according to specific formula and investor will lose some of the invested capital. There is an example brochure and term sheet available for this product in my GitHub page.

Library Imports

import QuantLib as ql
import numpy as np
import scipy.optimize as opt

Heston Path Generator

Below is a simple (hard-coded) method for generating paths by using Heston process for a given set of QuantLib dates, which can be unevenly distributed. A couple of notes about Heston path generation process in general. QuantLib MultiPathGenerator has to be used for simulating paths for this specific process, because "under the hood", QuantLib simulates two correlated random variates: one for the asset and another for the volatility. As a result for "one simulation round", MultiPath object will be received from generator. This MultiPath object contains generated paths for the both asset and volatility (multiPath[0] = asset path, multiPath[1] = volatility path). So, for any further purposes (say, to use simulated asset path for something), one can just use generated asset path and ignore volatility path. However, volatility path is available, if there will be some specific purpose for obtaining the information on how volatility has been evolving over time.

# hard-coded generator for Heston process
def HestonPathGenerator(dates, dayCounter, process, nPaths):
    t = np.array([dayCounter.yearFraction(dates[0], d) for d in dates])
    nGridSteps = (t.shape[0] - 1) * 2
    sequenceGenerator = ql.UniformRandomSequenceGenerator(nGridSteps, ql.UniformRandomGenerator())
    gaussianSequenceGenerator = ql.GaussianRandomSequenceGenerator(sequenceGenerator)
    pathGenerator = ql.GaussianMultiPathGenerator(process, t, gaussianSequenceGenerator, False)
    paths = np.zeros(shape = (nPaths, t.shape[0]))
    
    for i in range(nPaths):
        multiPath = pathGenerator.next().value()
        paths[i,:] = np.array(list(multiPath[0]))
        
    # return array dimensions: [number of paths, number of items in t array]
    return paths

Heston Model Calibration

Below is a simple (hard-coded) method for calibrating Heston model into a given volatility surface. Inside this method, process, model and engine are being created. After this, calibration helpers for Heston model are being created by using given volatility surface data. Finally, calibrated model and process are being returned for any further use. The actual optimization workhorse will be given outside of this method. For this specific example program, SciPy's Differential Evolution solver is being used, in order to guarantee global minimization result.

# hard-coded calibrator for Heston model
def HestonModelCalibrator(valuationDate, calendar, spot, curveHandle, dividendHandle, 
    v0, kappa, theta, sigma, rho, expiration_dates, strikes, data, optimizer, bounds):
    
    # container for heston calibration helpers
    helpers = []
    
    # create Heston process, model and pricing engine
    # use given initial parameters for model
    process = ql.HestonProcess(curveHandle, dividendHandle, 
        ql.QuoteHandle(ql.SimpleQuote(spot)), v0, kappa, theta, sigma, rho)
    model = ql.HestonModel(process)
    engine = ql.AnalyticHestonEngine(model)
    
    # nested cost function for model optimization
    def CostFunction(x):
        parameters = ql.Array(list(x))        
        model.setParams(parameters)
        error = [helper.calibrationError() for helper in helpers]
        return np.sqrt(np.sum(np.abs(error)))

    # create Heston calibration helpers, set pricing engines
    for i in range(len(expiration_dates)):
        for j in range(len(strikes)):
            expiration = expiration_dates[i]
            days = expiration - valuationDate
            period = ql.Period(days, ql.Days)
            vol = data[i][j]
            strike = strikes[j]
            helper = ql.HestonModelHelper(period, calendar, spot, strike,
                ql.QuoteHandle(ql.SimpleQuote(vol)), curveHandle, dividendHandle)
            helper.setPricingEngine(engine)
            helpers.append(helper)
    
    # run optimization, return calibrated model and process
    optimizer(CostFunction, bounds)
    return process, model

Monte Carlo Valuation

The actual Autocallable valuation algorithm has been implemented in this part of the code. For the sake of being able to value this product also after its inception, valuation method takes past fixings as Python dictionary (key: date, value: fixing). Now, if one changes valuation date to any possible date after inception and before transaction maturity, the product will also be valued accordingly. It should be noted, that in such scheme it is crucial to provide all possible past fixings. Failure to do this, will lead to an exception. Algorithm implementation is a bit long and might look scary, but it is actually really straightforward if one is familiar enough with this specific product.

def AutoCallableNote(valuationDate, couponDates, strike, pastFixings, 
    autoCallBarrier, couponBarrier, protectionBarrier, hasMemory, finalRedemptionFormula, 
    coupon, notional, dayCounter, process, generator, nPaths, curve):    
    
    # immediate exit trigger for matured transaction
    if(valuationDate >= couponDates[-1]): return 0.0
    
    # immediate exit trigger for any past autocall event
    if(valuationDate >= couponDates[0]):
        if(max(pastFixings.values()) >= (autoCallBarrier * strike)): return 0.0

    # create date array for path generator
    # combine valuation date and all the remaining coupon dates
    dates = np.hstack((np.array([valuationDate]), couponDates[couponDates > valuationDate]))
    
    # generate paths for a given set of dates, exclude the current spot rate
    paths = generator(dates, dayCounter, process, nPaths)[:,1:]
    
    # identify the past coupon dates
    pastDates = couponDates[couponDates <= valuationDate]

    # conditionally, merge given past fixings from a given dictionary and generated paths
    if(pastDates.shape[0] > 0):
        pastFixingsArray = np.array([pastFixings[pastDate] for pastDate in pastDates])        
        pastFixingsArray = np.tile(pastFixingsArray, (paths.shape[0], 1))
        paths = np.hstack((pastFixingsArray, paths))
    
    # result accumulator
    global_pv = []
    expirationDate = couponDates[-1]
    hasMemory = int(hasMemory)
    
    # loop through all simulated paths
    for path in paths:
        payoffPV = 0.0
        unpaidCoupons = 0
        hasAutoCalled = False
        
        # loop through set of coupon dates and index ratios
        for date, index in zip(couponDates, (path / strike)):
            # if autocall event has been triggered, immediate exit from this path
            if(hasAutoCalled): break
            payoff = 0.0
                
            # payoff calculation at expiration
            if(date == expirationDate):
                # index is greater or equal to coupon barrier
                # pay 100% redemption, plus coupon, plus conditionally all unpaid coupons
                if(index >= couponBarrier):
                    payoff = notional * (1 + (coupon * (1 + unpaidCoupons * hasMemory)))
                # index is greater or equal to protection barrier and less than coupon barrier
                # pay 100% redemption, no coupon
                if((index >= protectionBarrier) & (index < couponBarrier)):
                    payoff = notional
                # index is less than protection barrier
                # pay redemption according to formula, no coupon
                if(index < protectionBarrier):
                    # note: calculate index value from index ratio
                    index = index * strike
                    payoff = notional * finalRedemptionFormula(index)
                
            # payoff calculation before expiration
            else:
                # index is greater or equal to autocall barrier
                # autocall will happen before expiration
                # pay 100% redemption, plus coupon, plus conditionally all unpaid coupons
                if(index >= autoCallBarrier):
                    payoff = notional * (1 + (coupon * (1 + unpaidCoupons * hasMemory)))
                    hasAutoCalled = True
                # index is greater or equal to coupon barrier and less than autocall barrier
                # autocall will not happen
                # pay coupon, plus conditionally all unpaid coupons
                if((index >= couponBarrier) & (index < autoCallBarrier)):
                    payoff = notional * (coupon * (1 + unpaidCoupons * hasMemory))
                    unpaidCoupons = 0
                # index is less than coupon barrier
                # autocall will not happen
                # no coupon payment, only accumulate unpaid coupons
                if(index < couponBarrier):
                    payoff = 0.0
                    unpaidCoupons += 1                    

            # conditionally, calculate PV for period payoff, add PV to local accumulator
            if(date > valuationDate):
                df = curveHandle.discount(date)
                payoffPV += payoff * df
            
        # add path PV to global accumulator
        global_pv.append(payoffPV)
        
    # return PV
    return np.mean(np.array(global_pv))

Main Program

Finally, in this part of the program we will actually use the stuff presented above. First, the usual QuantLib-related parameters are being created, as well as parameters for Autocallable product. Note, that since in this example we are valuing this product at inception, there is no need to provide any past fixings for valuation method (however, it is not forbidden either). After this, interest rate curve and dividend curve will be created, as well as volatility surface data for calibration purposes. Next, initial guesses for Heston parameters will be set and the model will then be calibrated by using dedicated method. Finally, calibrated process will be given to the actual valuation method and Monte Carlo valuation will be processed.

# general QuantLib-related parameters
valuationDate = ql.Date(20,11,2019)
ql.Settings.instance().evaluationDate = valuationDate
convention = ql.ModifiedFollowing
dayCounter = ql.Actual360()
calendar = ql.TARGET()

# Autocallable Memory Coupon Note
notional = 1000000.0
spot = 3550.0
strike = 3550.0
autoCallBarrier = 1.0
couponBarrier = 0.8
protectionBarrier = 0.6
finalRedemptionFormula = lambda indexAtMaturity: min(1.0, indexAtMaturity / strike)
coupon = 0.05
hasMemory = True

# coupon schedule for note
startDate = ql.Date(20,11,2019)
firstCouponDate = calendar.advance(startDate, ql.Period(1, ql.Years))
lastCouponDate = calendar.advance(startDate, ql.Period(7, ql.Years))
couponDates = np.array(list(ql.Schedule(firstCouponDate, lastCouponDate, ql.Period(ql.Annual), 
    calendar, ql.ModifiedFollowing, ql.ModifiedFollowing, ql.DateGeneration.Forward, False)))

# create past fixings into dictionary
pastFixings = {}
#pastFixings = { ql.Date(20,11,2020): 99.0, ql.Date(22,11,2021): 99.0 }

# create discounting curve and dividend curve, required for Heston model
curveHandle = ql.YieldTermStructureHandle(ql.FlatForward(valuationDate, 0.01, dayCounter))
dividendHandle = ql.YieldTermStructureHandle(ql.FlatForward(valuationDate, 0.0, dayCounter))

# Eurostoxx 50 volatility surface data
expiration_dates = [ql.Date(19,6,2020), ql.Date(18,12,2020), 
    ql.Date(18,6,2021), ql.Date(17,12,2021), ql.Date(17,6,2022),
    ql.Date(16,12,2022), ql.Date(15,12,2023), ql.Date(20,12,2024), 
    ql.Date(19,12,2025), ql.Date(18,12,2026)]

strikes = [3075, 3200, 3350, 3550, 3775, 3950, 4050]

data = [[0.1753, 0.1631, 0.1493, 0.132 , 0.116 , 0.108 , 0.1052],
       [0.1683, 0.1583, 0.147 , 0.1334, 0.1212, 0.1145, 0.1117],
       [0.1673, 0.1597, 0.1517, 0.1428, 0.1346, 0.129 , 0.1262],
       [0.1659, 0.1601, 0.1541, 0.1474, 0.1417, 0.1381, 0.1363],
       [0.1678, 0.1634, 0.1588, 0.1537, 0.1493, 0.1467, 0.1455],
       [0.1678, 0.1644, 0.1609, 0.1572, 0.1541, 0.1522, 0.1513],
       [0.1694, 0.1666, 0.1638, 0.1608, 0.1584, 0.1569, 0.1562],
       [0.1701, 0.168 , 0.166 , 0.164 , 0.1623, 0.1614, 0.161 ],
       [0.1715, 0.1698, 0.1682, 0.1667, 0.1654, 0.1648, 0.1645],
       [0.1724, 0.171 , 0.1697, 0.1684, 0.1675, 0.1671, 0.1669]]

# initial parameters for Heston model
theta = 0.01
kappa = 0.01
sigma = 0.01
rho = 0.01
v0 = 0.01

# bounds for model parameters (1=theta, 2=kappa, 3=sigma, 4=rho, 5=v0)
bounds = [(0.01, 1.0), (0.01, 10.0), (0.01, 1.0), (-1.0, 1.0), (0.01, 1.0)]

# calibrate Heston model, print calibrated parameters
calibrationResult = HestonModelCalibrator(valuationDate, calendar, spot, curveHandle, dividendHandle, 
        v0, kappa, theta, sigma, rho, expiration_dates, strikes, data, opt.differential_evolution, bounds)
print('calibrated Heston parameters', calibrationResult[1].params())

# monte carlo parameters
nPaths = 10000

# request and print PV
PV = AutoCallableNote(valuationDate, couponDates, strike, pastFixings, 
    autoCallBarrier, couponBarrier, protectionBarrier, hasMemory, finalRedemptionFormula, 
    coupon, notional, dayCounter, calibrationResult[0], HestonPathGenerator, nPaths, curveHandle)

print(PV)

Finally, thanks for reading this blog and have a pleasant wait for the coming Christmas.
-Mike

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

QuantLib-Python: Path Generator Method for Uncorrelated and Correlated 1-D Stochastic Processes

This Python program presents one compact method for simulating paths for the both uncorrelated and correlated stochastic processes.

Thanks for reading my blog.
-Mike

%config IPCompleter.greedy = True
from QuantLib import *
import numpy as Numpy
import matplotlib.pyplot as Matplotlib

# method for simulating paths for the both uncorrelated and correlated processes
# arguments:
# process = QuantLib 1-dimensional stochastic process object or 
#           StochasticProcessArray (Array of correlated 1-D stochastic processes)
# timeGrid = QuantLib TimeGrid object
# n = number of paths
def GeneratePaths(process, timeGrid, n):
    
    # correlated processes, use GaussianMultiPathGenerator
    if(type(process) == StochasticProcessArray):
        times = []; [times.append(timeGrid[t]) for t in range(len(timeGrid))]        
        nGridSteps = (len(times) - 1) * process.size()
        sequenceGenerator = UniformRandomSequenceGenerator(nGridSteps, UniformRandomGenerator())
        gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
        pathGenerator = GaussianMultiPathGenerator(process, times, gaussianSequenceGenerator, False)        
        paths = Numpy.zeros(shape = (n, process.size(), len(timeGrid)))
        
        # loop through number of paths
        for i in range(n):
            # request multiPath, which contains the list of paths for each process
            multiPath = pathGenerator.next().value()
            # loop through number of processes
            for j in range(multiPath.assetNumber()):
                # request path, which contains the list of simulated prices for a process
                path = multiPath[j]
                # push prices to array
                paths[i, j, :] = Numpy.array([path[k] for k in range(len(path))])
        # resulting array dimension: n, process.size(), len(timeGrid)
        return paths

    # uncorrelated processes, use GaussianPathGenerator
    else:
        sequenceGenerator = UniformRandomSequenceGenerator(len(timeGrid), UniformRandomGenerator())
        gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
        maturity = timeGrid[len(timeGrid) - 1]
        pathGenerator = GaussianPathGenerator(process, maturity, len(timeGrid), gaussianSequenceGenerator, False)
        paths = Numpy.zeros(shape = (n, len(timeGrid)))   
        for i in range(n):
            path = pathGenerator.next().value()
            paths[i, :] = Numpy.array([path[j] for j in range(len(timeGrid))])
        # resulting array dimension: n, len(timeGrid)
        return paths

# create simulation-related parameters
today = Date(30, November, 2018)
maturity = 5.0
nSteps = int(maturity) * 365
# create regularly spaced QuantLib TimeGrid object
timeGrid = TimeGrid(maturity, nSteps)
nPaths = 25

# create HW1F model
reversionSpeed = 0.05
rateVolatility = 0.0099255
curve = RelinkableYieldTermStructureHandle(FlatForward(today, 0.01, Actual360()))
HW1F = HullWhiteProcess(curve, reversionSpeed, rateVolatility)
hw1f_paths = GeneratePaths(HW1F, timeGrid, nPaths)

# create GBM model
initialValue = 0.01
mue = 0.01
sigma = 0.0099255
GBM = GeometricBrownianMotionProcess(initialValue, mue, sigma)
gbm_paths = GeneratePaths(GBM, timeGrid, nPaths)

# plot uncorrelated paths
times = []; [times.append(timeGrid[t]) for t in range(len(timeGrid))]
Matplotlib.rcParams['figure.figsize'] = [12.0, 8.0]
f, subPlots = Matplotlib.subplots(2, sharex = True)
f.suptitle('Uncorrelated paths n=' + str(nPaths))
subPlots[0].set_title('HW1F')
subPlots[1].set_title('GBM')

for i in range(hw1f_paths.shape[0]):
    path = hw1f_paths[i, :] 
    subPlots[0].plot(times, path)

for i in range(gbm_paths.shape[0]):
    path = gbm_paths[i, :] 
    subPlots[1].plot(times, path)

# create correlated paths
rho = 1.0
correlation = [[1.0, rho], [rho, 1.0]]
processArray = StochasticProcessArray([HW1F, GBM], correlation)
correlated_paths = GeneratePaths(processArray, timeGrid, nPaths)

# plot correlated paths
f2, subPlots2 = Matplotlib.subplots(processArray.size(), sharex = True)
f2.suptitle('Correlated paths n=' + str(nPaths) + ', rho=' + str(rho))
subPlots2[0].set_title('HW1F')
subPlots2[1].set_title('GBM')

for i in range(nPaths):
    for j in range(processArray.size()):
        path = correlated_paths[i, j, :]
        subPlots2[j].plot(times, path)

QuantLib-Python: Simulating Paths for Correlated 1-D Stochastic Processes

This program, which is just an extension to my previous post, will create two correlated Geometric Brownian Motion processes, then request simulated paths from dedicated generator function and finally, plots all simulated paths to charts. For the two processes in this example program, correlation has been set to minus one and total of 20 paths has been requested for the both processes. Each path has maturity of one year, which has been divided into 90 time steps. This example can be generalized for higher dimensions by adjusting size of a given correlation matrix and corresponding lists of other parameters (nProcesses, names, spot, mue, sigma).

Thanks for reading my blog.
-Mike 

%config IPCompleter.greedy = True
from QuantLib import *
import numpy as Numpy
import matplotlib.pyplot as Matplotlib

# processArray = StochasticProcessArray (Array of correlated 1-D stochastic processes)
# timeGrid = TimeGrid object
def GenerateCorrelatedPaths(processArray, timeGrid, nPaths):
    times = []; [times.append(timeGrid[t]) for t in range(len(timeGrid))]
    generator = UniformRandomGenerator()
    nProcesses = processArray.size()
    nGridSteps = len(times) - 1 # deduct initial time (0.0)
    nSteps = nGridSteps * nProcesses
    sequenceGenerator = UniformRandomSequenceGenerator(nSteps, generator)
    gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
    multiPathGenerator = GaussianMultiPathGenerator(processArray, times, gaussianSequenceGenerator)
    paths = Numpy.zeros(shape = (nPaths, nProcesses, len(timeGrid)))

    # loop through number of paths
    for i in range(nPaths):
        # request multiPath, which contains the list of paths for each process
        multiPath = multiPathGenerator.next().value()
        # loop through number of processes
        for j in range(multiPath.assetNumber()):
            # request path, which contains the list of simulated prices for a process
            path = multiPath[j]
            # push prices to array
            paths[i, j, :] = Numpy.array([path[k] for k in range(len(path))])
    return paths

# create two 1-D stochastic processes
process = []
nProcesses = 2
correlation = -1.0
names = ['equity_1', 'equity_2']
spot = [100.0, 100.0]
mue = [0.01, 0.01]
sigma = [0.10, 0.10]
[process.append(GeometricBrownianMotionProcess(spot[i], mue[i], sigma[i])) for i in range(nProcesses)]
matrix = [[1.0, correlation], [correlation, 1.0]]

# create timegrid object and define number of paths
maturity = 1.0
nSteps = 90
timeGrid = TimeGrid(maturity, nSteps)
nPaths = 20

# create StochasticProcessArray object
# (array of correlated 1-D stochastic processes)
processArray = StochasticProcessArray(process, matrix)
# request simulated correlated paths for all processes
# result array dimensions: nPaths, nProcesses, len(timeGrid)
paths = GenerateCorrelatedPaths(processArray, timeGrid, nPaths)

# plot paths
f, subPlots = Matplotlib.subplots(nProcesses, sharex = True)
f.suptitle('Path simulations rho=' + str(correlation) + ', n=' + str(nPaths))

for i in range(nPaths):
    for j in range(nProcesses):
        subPlots[j].set_title(names[j])
        path = paths[i, j, :]
        subPlots[j].plot(timeGrid, path)

QuantLib-Python: Simulating Paths for 1-D Stochastic Processes

This simple Python program will create two 1-dimensional stochastic process objects (Hull-White 1-Factor and Geometric Brownian Motion), then request simulated paths from dedicated generator function and finally, plots all simulated paths to charts. The original C++ implementation can be found in here.

Thanks for reading my blog.
-Mike

%config IPCompleter.greedy = True
from QuantLib import *
import numpy as Numpy
import matplotlib.pyplot as Matplotlib

# process = QuantLib 1-dimensional stochastic process object
def GeneratePaths(process, maturity, nPaths, nSteps):
    generator = UniformRandomGenerator()
    sequenceGenerator = UniformRandomSequenceGenerator(nSteps, generator)
    gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
    paths = Numpy.zeros(shape = ((nPaths), nSteps + 1))
    pathGenerator = GaussianPathGenerator(process, maturity, nSteps, gaussianSequenceGenerator, False)
    for i in range(nPaths):
        path = pathGenerator.next().value()
        paths[i, :] = Numpy.array([path[j] for j in range(nSteps + 1)])
    return paths


# general parameters and objects
tradeDate = Date(23, November, 2018)
Settings_instance().evaluationDate = tradeDate
dayCounter = Actual360()
calendar = UnitedStates()
settlementDate = calendar.advance(tradeDate, 2, Days)

# common simulation-related parameters for all processes
maturity = 3.0
nPaths = 50
nSteps = int(maturity * 365)
timeGrid = Numpy.linspace(0.0, maturity, nSteps + 1)

# create HW1F model, request paths from generator
reversionSpeed = 0.05
rateVolatility = 0.0099255
r = QuoteHandle(SimpleQuote(0.01))
curve = RelinkableYieldTermStructureHandle(FlatForward(settlementDate, r, dayCounter))
HW1F = HullWhiteProcess(curve, reversionSpeed, rateVolatility)
hw1f_paths = GeneratePaths(HW1F, maturity, nPaths, nSteps)

# create GBM model, request paths from generator
initialValue = 0.01
mue = 0.01
sigma = 0.0099255
GBM = GeometricBrownianMotionProcess(initialValue, mue, sigma)
gbm_paths = GeneratePaths(GBM, maturity, nPaths, nSteps)

# plot all paths for the both processes
f, subPlots = Matplotlib.subplots(2, sharex = True)
Matplotlib.rcParams['figure.figsize'] = [16.0, 10.0]
f.suptitle('Path simulations n=' + str(nPaths))
subPlots[0].set_title('Hull-White 1-Factor')
subPlots[1].set_title('Geometric Brownian Motion')

for i in range(hw1f_paths.shape[0]):
    path = hw1f_paths[i, :] 
    subPlots[0].plot(timeGrid, path)

for i in range(gbm_paths.shape[0]):
    path = gbm_paths[i, :] 
    subPlots[1].plot(timeGrid, path)

Excel : implementing multithreading using C# and Excel-DNA

Multithreading possibilities for VBA has been chewed in this post. My conclusion back then was, that there is no multi-threading possibilities for VBA. One reader has then been commenting this as follows :

"Correct to say there is no multithreading in VBA however, you can very easily create your multithreading functionality in a c# dll and expose the interface via Interops for VBA to consume. Then all you need to do is to pass in the data into the function and walla, you have multithreading. Interops is a very powerful feature and VBA fully supports it."

For the sake of completeness, I will present one possible program to implement such multithreading scheme for Excel, by using Parallel class for the actual multi-threading part and Excel-DNA for interfacing the program back to Excel (using Excel only as Input-Output platform). Excel interfacing part used in this post, has already been completely chewed within this post. Note, that by using instructions given in here or here, one can relatively easy interface this multithreading program back to VBA, instead of Excel.

The program

This unsurprising program is simply processing a task (simulate N amount of normally distributed random numbers) shared for M processors. Amount of simulated numbers and number of processors used is decided by the client (Excel). Finally, the program will print the aggregated time consumed for this task back to Excel. Needless to say, the effect of multithreading can be demonstrated by changing the number of processors in Excel worksheet.

Copy-Paste the following program presented below into a new C# project. Also, follow carefully all instructions given in this post.

Thanks for reading this blog.
-Mike

// ExcelInterface.cs
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Threading.Tasks;
using System.Windows.Forms;
using ExcelDna.Integration;

namespace ExcelInterface {
    // type definition, delegate for exporting data from processor to queue
    using Processor = Action<int, DataMessenger>;
    public delegate void DataMessenger(List<double> result);

    public static class ExcelInterface {
        // method available in excel
        public static void execute() {
            try {
                // create Excel application object
                dynamic Excel = ExcelDnaUtil.Application;
                //
                // get parameters from worksheet named ranges
                int n = (int)Excel.Range["_n"].Value2;
                int nProcessors = (int)Excel.Range["_nProcessors"].Value2;
                int nPerProcessor = n / nProcessors;
                // create data thread-safe queue and random processors
                RandomQueue messageQueue = new RandomQueue();
                List<Processor> processors = new List<Processor>();
                for (int i = 0; i < nProcessors; ++i) {
                    processors.Add(Factory.CreateProcessor());
                }
                // host timer, execute all processors in parallel
                Stopwatch timer = new Stopwatch();
                timer.Start();
                Parallel.For(0, nProcessors, i => {
                    processors[i].Invoke(nPerProcessor, messageQueue.Enqueue); 
                });
                timer.Stop();
                //
                // print processing time to excel
                double timeElapsed = timer.Elapsed.TotalSeconds;
                Excel.Range["_timer"] = timeElapsed;
            }
            catch (Exception e) {
                MessageBox.Show(e.Message);
            }
        }
    }
}

// RandomQueue.cs
using System;
using System.Collections.Concurrent;
using System.Collections.Generic;

namespace ExcelInterface {
    // wrapper class for ConcurrentQueue data structure
    public class RandomQueue {
        private ConcurrentQueue<List<double>> randomQueue;
        public RandomQueue() {
            randomQueue = new ConcurrentQueue<List<double>>();
        }
        // add
        public void Enqueue(List<double> data) {
            randomQueue.Enqueue(data);
        }
        // remove
        public List<double> Dequeue() {
            List<double> data = null;
            bool hasValue = false;
            while (!hasValue) hasValue = randomQueue.TryDequeue(out data);
            return data;
        }
    }
}

// Factory.cs
using System;
using System.Collections.Generic;

namespace ExcelInterface {
    // type definition
    using Processor = Action<int, DataMessenger>;
    
    public static class Factory {
        // factory method for creating processor for random generator
        public static Processor CreateProcessor() {
            //
            // create standard normal variable generator
            // using action delegate and lambda expression
            Processor process =
            (int n, DataMessenger dataMessenger) => {
                Random generator = new Random(Guid.NewGuid().GetHashCode());
                List<double> normalRandomList = new List<double>();
                // loop through n
                for (int i = 0; i < n; ++i) {
                    double uniformRandomSum = 0.0;
                    // create standard normal random variable using CLT
                    for (int j = 0; j < 12; ++j) {
                        double uniformRandom = 0.0;
                        while (true) {
                            uniformRandom = generator.NextDouble();
                            if (uniformRandom > 0.0) break;
                        }
                        uniformRandomSum += uniformRandom;
                    }
                    normalRandomList.Add(uniformRandomSum - 6.0);
                }
                // send list of normal random variables to queue
                dataMessenger(normalRandomList);
            };
            return process;
        }
    }
}

C++11 : Multi-Threaded PathGenerator using PPL

 

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

C++11 : Multi-Threaded Random Generator using PPL

After some woodshedding with C++11 random and publishing distributional random generator in my blog a couple of weeks ago, the next question was, how to make this thing to happen a bit faster? Since I already have been snooping PPL (Parallel Patterns Library), I realized that Parallel Algorithms in that library would be an excellent candidate for this task. There is a great and comprehensive internet book available on this topic. After getting familiar with some of the issues presented in these sources, one should be able to parallelize algorithms succesfully.

There are some very specific issues, related to random generator seeding. As I started to parallelize my random number generation (creating several generators and processing these using parallel algorithm), I realized that my random generators were not working as expected. After some further woodshedding I came to conclusion, that (in a multi-threaded setting) each generator has to be seeded independently. Finally, concerning my current generator implementation, I made a decision that a client has to provide explicit seed values for all generators. Based on my hands-on experience with a system used for counterparty credit risk simulations (widely ranked as state-of-the-art numerical software), I dare to say that this approach is well justified. Now, the question of what kind of tools will be used for creating the actual seed values, has been left to a client judgement. The rabbit hole of these profound issues related to random generators have been chewed in this blog. Do not enter, if you are not absolutely ready to take the red pill.

ROUGH DIAGNOSTICS

I 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. The average running time for this sample was 5187 milliseconds (sequential) and 2317 milliseconds (parallel). So, parallel scheme will be completed about two times faster. The actual CPU usages during simulation process for the both schemes are shown in the screenshot below. When using parallel algorithm, computer is able to utilize all four cores in my computer.

I also tested all processed random paths for the both schemes, 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 paths for the both schemes are different and there are no duplicates. Presented scheme is hereby fulfilling my two requirements : faster creation of random paths and independency of random generators.

RandomGenerator.h

Current implementation for distributional random path generator is presented below. Note, that the actual seeder creation takes place outside of this class.

#pragma once
#include <algorithm>
#include <functional>
#include <vector>
#include <random>
//
namespace MikeJuniperhillRandomGeneratorTemplate
{
 template <typename Generator = std::mt19937, typename Distribution = std::normal_distribution<double>>
 /// <summary>
 /// Template class for creating random number paths using
 /// Mersenne Twister as default uniform random generator and
 /// Standard Normal (0.0, 1.0) as default probability distribution.
 /// </summary> 
 class RandomGenerator
 {
 public:
  /// <summary>
  /// set explicit seed value
  /// </summary>
  RandomGenerator(unsigned long seed)
  {
   // construct lambda method for processing standard normal random number
   randomGenerator = [this](double x)-> double
   {
    x = distribution(uniformGenerator);
    return x;
   };
   // seed generator
   uniformGenerator.seed(seed);
  }
  /// <summary>
  /// fill vector reference with distributional random numbers
  /// </summary> 
  void operator()(std::vector<double>& v)
  {
   std::transform(v.begin(), v.end(), v.begin(), randomGenerator);
  }
 private:
  std::function<double(double)> randomGenerator;
  Generator uniformGenerator;
  Distribution distribution;
 };
}
//

tester.cpp

First, all simulation parameters and explicit seed values for all generators are defined. Simple factory method is then used to create the actual generator function as lambda function. Next, factory method is returning generator function wrapped inside a shared pointer. Main program is then storing each created generator function into vector container and finally processing all generators (sequentially or parallel). Each generator is storing generated paths into thread-safe nested vector (std::vector inside concurrency::concurrent_vector).

#include <iostream>
#include <memory>
#include <chrono>
#include <ppl.h>
#include <concurrent_vector.h>
#include "RandomGenerator.h"
namespace MJ = MikeJuniperhillRandomGeneratorTemplate;
//
//
// thread-safe nested vector for storing random paths processed by generators
concurrency::concurrent_vector<std::vector<double>> paths;
//
// prototype for generator function
using Generator = std::function<void(void)> ;
//
// factory method :
// return generator function as function wrapper
std::shared_ptr<Generator> Factory(int nPaths, int nSteps, unsigned long seed, 
 concurrency::concurrent_vector<std::vector<double>>& paths)
{
 // create generator function
 auto generator = [=, &paths]()->void
 {
  // step 1 : create distributional random generator
  MJ::RandomGenerator<> randomGenerator(seed);
  // step 2 : create container for a random path
  std::vector<double> path(nSteps);
  // step 3 : create path using distributional generator
  // and import created path into thread-safe container
  for (auto i = 0; i != nPaths; ++i)
  {
   randomGenerator(path);
   paths.push_back(path);
  }
 };
 // return generator function as function wrapper
 return std::shared_ptr<Generator>(new Generator(generator));
}
//
int main()
{
 // define parameters and seed values for generator functions
 int nGenerators = 4;
 int nPaths = 10000;
 int nSteps = 250;
 std::vector<unsigned long> seed = { 10322854, 65947, 387528, 772399573 };
 //
 // use factory method for creating generator function
 std::vector<std::shared_ptr<Generator>> generators;
 for (auto i = 0; i < nGenerators; i++) generators.push_back(Factory(nPaths, nSteps, seed[i], paths));
 //
 // parallel processing
 auto start = std::chrono::steady_clock::now();
 concurrency::parallel_for_each(generators.begin(), generators.end(),
  [](std::shared_ptr<Generator> ppg) -> void { (*ppg)(); });
 //
 //// uncomment for sequential processing
 //std::for_each(generators.begin(), generators.end(),
 // [](std::shared_ptr<Generator> ppg) -> void { (*ppg)(); });
 //
 auto end = std::chrono::steady_clock::now();
 auto timeElapsed = (end - start);
 std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(timeElapsed).count() << std::endl;
 return 0;
}
//

Finally, thanks for reading my blog. I would like to wish Happy New Year for everybody.
-Mike

XLW : Interfacing C++11 PathGenerator to Excel

In my previous post on generating paths for different types of one-factor processes, I was writing all processed paths sequentially into separate CSV files. Later, these files could be opened in Excel for any further use. This approach is far from being an efficient way to interface C++ program to Excel. While being possible to create in-house tool for this task, the actual implementation part requires relatively deep understanding and hands-on experience on Excel/C API. The amount of issues to be learned in order to produce truly generic and reliable tool, is far from being something which could be internalized in very short period of time. There are limited amount of books available on this topic, but the most recommended book is the one written by Steve Dalton.

This time, I wanted to present kind of "industry standard" way to accomplish this task using "easy-to-use" Excel/C++ interfacing tool, which has been there already since 2002. XLW is an application which wraps Excel/C API into a simple C++ interface, which can then be used to customize Excel with user-defined worksheet functions. For any newcomer into this issue, it is highly recommended first to watch instructional "how-to-use" clips from project homepage. The complete workflow of using XLW wrapper is presented there from downloading to debugging. Also, AdaptiveRisk blog is presenting extremely useful stuff, well enough to get started.

cppinterface.h

For this project, I have extracted a new XLW xll template and opened corresponding solution in my Visual Studio Express 2013. The content of this header file in my current project is presented below. I have method declarations for two functions, which are both returning matrix object (XLW data structure). It is my intention to use my previous PathGenerator project, in order to fill these matrix objects with paths using desired one-factor processes.

#ifndef TEST_H
#define TEST_H
//
#include "xlw/MyContainers.h"
#include <xlw/CellMatrix.h>
#include <xlw/DoubleOrNothing.h>
#include <xlw/ArgList.h>
#include <xlw/XlfServices.h>
//
using namespace xlw;
//
//<xlw:libraryname=XLWPathGenerator
//
// method for requesting vasicek paths
MyMatrix // return matrix of random paths following Vasicek SDE
//<xlw:volatile
GetPaths_Vasicek(double t, // time to maturity
double r, // current short rate
double longTermRate, // long-term average rate
double meanReversion, // mean reversion speed
double rateVolatility // rate volatility
);
//
// method for requesting GBM paths
MyMatrix // return matrix of random paths following Geometric Brownian Motion SDE
//<xlw:volatile
GetPaths_BrownianMotion(double t, // time to maturity
double s, // current spot rate
double rate, // risk-free rate
double volatility // volatility
);
//
#endif

Commenting may seem a bit strange first, but the following screenshot containing Excel function argument input box may help to catch the point.

source.cpp

Implementations for two methods declared in header file are presented below. Information concerning number of time steps (for a path) and number of paths to be created are extracted from matrix dimensions using XlfServices object. After this, desired OneFactorProcess and PathGenerator objects are created. Finally, PathGenerator object is used to process a path, which will be imported into resulting matrix object (paths) and returned for the client (Excel).

#include <cppinterface.h>
#include "PathGenerator.h"
#pragma warning (disable : 4996)
//
MyMatrix GetPaths_Vasicek(double t, double r, double longTermRate, double meanReversion, double rateVolatility)
{
 // request dimensions for calling matrix
 const unsigned int nPaths = XlfServices.Information.GetCallingCell().columns();
 const unsigned int nSteps = XlfServices.Information.GetCallingCell().rows();
 // create container for all processed paths
 MyMatrix paths(nSteps, nPaths);
 // create container for a single path to be processed
 MyArray path(nSteps);
 //
 // create vasicek process and path generator
 std::shared_ptr<MJProcess::OneFactorProcess> vasicek =
  std::shared_ptr<MJProcess::Vasicek>(new MJProcess::Vasicek(meanReversion, longTermRate, rateVolatility));
 PathGenerator<> shortRateProcess(r, t, vasicek);
 //
 // process paths using path generator
 for (unsigned int i = 0; i != nPaths; ++i)
 {
  shortRateProcess(path);
  // import processed path into paths container
  for (unsigned j = 0; j != nSteps; ++j)
  {
   paths[j][i] = path[j];
  }
 }
 return paths;
}
//
MyMatrix GetPaths_BrownianMotion(double t, double s, double rate, double volatility)
{
 // request dimensions for calling matrix
 const unsigned int nPaths = XlfServices.Information.GetCallingCell().columns();
 const unsigned int nSteps = XlfServices.Information.GetCallingCell().rows();
 // create container for all processed paths
 MyMatrix paths(nSteps, nPaths);
 // create container for a single path to be processed
 MyArray path(nSteps);
 //
 // create geometric brownian motion process and path generator
 std::shared_ptr<MJProcess::OneFactorProcess> brownianMotion =
  std::shared_ptr<MJProcess::GBM>(new MJProcess::GBM(rate, volatility));
 PathGenerator<> equityPriceProcess(s, t, brownianMotion);
 //
 // process paths using path generator
 for (unsigned int i = 0; i != nPaths; ++i)
 {
  equityPriceProcess(path);
  // import processed path into paths container
  for (unsigned j = 0; j != nSteps; ++j)
  {
   paths[j][i] = path[j];
  }
 }
 return paths;
}
//

In order to get this thing up and running, header file for PathGenerator has to be included. I have set the current XLL project as startup project. As a side, I have opened my PathGenerator project (containing header files for RandomGenerator, OneFactorProcess and PathGenerator). Since this side project is still unaccessible, it has to be linked to my current XLL project : Project - Properties - Configuration Properties - C/C++ - General - Additional Include Directories (Browse folder containing source files for the project to be linked). After completing these steps and building this project succesfully, I am finally ready to test the provided functionality in Excel.

Excel

After opening a new Excel, I need to drag-and-drop (or open from Excel) newly created xll template from my \\Projects\XLWTester\Debug folder to Excel. In my Excel (screenshot below), I have two boxes for input parameters (Vasicek, Brownian Motion) and two ranges for resulting one-factor process paths (36 steps, 15 paths). As soon as I hit F9 button, my one-factor paths will be re-created. Finally, it should be noted that the functions are both array formulas.

Finally, thanks for reading this blog. Pleasant waiting for a new year for everybody.
-Mike