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()

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)