Monte-Carlo Engine, seed parameter implementation and new numerical schemes for SDEs

[yfnaji]

This Pull Request was created to address issue 244.

Latest update

With these latest changes, the emphasis was on implementing the Monte-Carlo only once - so that we would only need to update the code in one place should we need to amend anything regarding the algorithm.

The following changes have been made:

• The euler_maruyama() method has been replaced with monte_carlo()
• The type of scheme (Euler-Maruyama, Milstein etc) can be chosen via the StochasticScheme enum as a parameter for StochasticConfig
• Monte-Carlo engine centralised in a private (public crate-wide) function run_monte_carlo() in monte_carlo.rs
• Jump diffusion processes have a simplified Monte-Carlo implementation (no need to reimplement the algorithm)
• The Monte-Carlo simulation for fractional stochastic processes has been centralised (no need to implement the same algorithm for each fractional process) and simplified
• Fractional stochastic process functionalities have been moved to the new fractional_process.rs module

Note: Jump/Fractional processes require a specific implementation of run_monte_carlo(). Although this works fine, there might be some issues when user's construct custom jump/fractional processes (like in examples/examples/custom_process.rs) as default implementation will overlook jump/fractional properties. Added to future works section below.

Benchmark tests

The following benchmark tests were run using the Euler-Maruyama scheme with the following configuration:

x_0 = 1.0
t_0 = 0.0
t_n = 1.0
n_steps = 1000

Standard Stochastic Process

GeometricBrownianMotion with

mu = 1.0
sigma = 0.2
m_path = 100

Branch	Run 1	Run 2	Run 3
monte-carlo	72.253916ms	77.608711ms	70.91479ms
main	110.545326ms	90.941425ms	101.169815ms

This branch on average is 27% faster than main with Geometric Brownian Motion.

Jump Diffusion

MertonJumpDiffusion with

mu = 1.0
sigma = 0.2
lambda = 1
m = 0.0
v = 0.1
m_path = 100

Branch	Run 1	Run 2	Run 3
monte-carlo	192.875529ms	185.308736ms	195.708128ms
main	171.309267ms	185.378698ms	159.632454ms

This branch is on average 11% slower than main with Merton Jump Diffusion.
We should consider whether generalisation of the Monte-Carlo method is worth this cost.
Of course, we can try to improve performance in a follow up (added to Future work section at the bottom).

Fractional Process

FractionalBrownianMotion with

hurst = 0.1
method = FractionalProcessGeneratorMethod::Cholesky
m_path = 10

Branch	Run 1	Run 2	Run 3
monte-carlo	98.090812801s	105.315690374s	85.730874942s
main	119.905169119s	89.912333031s	97.421714688s

The performance on a Fraction Brownian motion process is roughly the same.

Centralize Monte-Carlo functionality

Previously, the monte-carlo method needed to be implemented for each numerical method individually. Now there exists a private monte_carlo() function in the process.rs module which can be utilised for each numerical method defined in the StochasticProcess trait, removing the need for re-implementing the Monte-Carlo method every time.

Consolidate seeded and "unseeded" numerical methods

The Euler Maruyama method had two implementations, one where a seed can be defined and another where the seed is "randomly" generated, seedable_euler_maruyama and euler_maruyama respectively.

A seed: Option<u64> parameter is introduced into StochasticProcessConfig to give the user the ability to fix outcomes, other can be entered as None, removing the need for two implementations of the same numerical methods.

Milstein's scheme & Strang Splitting

The Milstein's scheme and Strang Splitting are numerical methods for SDEs that have been implemented in process.rs as part of this PR (and utilize monte_carlo() mentioned before). As such, price_monte_carlo() has been amended slightly to accomodate the new numerical methods.

Their methodologies are outlined below.

Define the SDE

$$ dX = A(X, t)dt + B(X, t) dW_t $$

where $W_t\sim N\left(0, t\right)$ is the Wiener process. We define $t_n = t_0 + n\Delta t$ where $\Delta t = \frac{t_N-t_0}{M}$.

Milstein's Scheme

The Milstein scheme at time $t_n\in\left[t_0, t_N\right]$ is defined as

$$ X_{n+1} = X_n + A(X_n, t_n)\Delta t + B(X_n, t_n)\Delta W_t + \frac{1}{2}B(X_n, t_n)\frac{\partial}{\partial X}B\left(X_n, t_n\right)\left(\Delta W_n^2 - \Delta t\right) $$

with the initial condition $X_0 = x_0$.

Strang Splitter

The Strang Splitter method can be split into three steps:

1. Calculate the drift term over half a time step:

$$ X_{n+1/2} = X_n + \frac{1}{2}a\left(X_n, t\right)\Delta t $$

2. Calculate the Diffusion term:

$$ X_{n+1/2} =X_{n+1/2} + b\left(X_{n+1/2}, t+\frac{1}{2}\Delta t\right)\Delta W_t $$

3. Complete the full drift step:

$$ X_{n+1} = X_{n+1/2} + \frac{1}{2}a\left(X_{n+1/2}, t+\Delta t\right)\Delta t $$

All appropriate unit tests have been added.

Future work

• Generalise monte_carlo() to handle stochastic volatility
• Bring Heston model to life in order to test monte carlo methods on stochastic volatility models
• Enable monte_carlo() for custom fractional processes
• Improve performance of jump diffusion processes
• Enable user's to utilise the appropriate monte carlo method when creating customised jump/fractional processes

Although these points are relevant, they will not be addressed as part of this PR to avoid scope screep.

[avhz]

I will hopefully review and merge this PR this weekend :) Thanks for this one !

[avhz]

Sorry it's taken me a while to have a look at this.
Just a couple of things:

• Why not put the StochasticMethod into the config ?
• Did you benchmark it ?

[yfnaji]

Hey @avhz

I have moved StochasticMethod (now called StochasticScheme) into the config + added some benchmark tests to the description.

I have also included other changes that I have noted at the top of the description.

[avhz]

Looks good :)

Just one thing, I would change the method names to a verb, something like process.generate() rather than process.monte_carlo(), we can reserve Monte Carlo for pricing and other estimation, rather than path simulation.

[yfnaji]

Thanks for having a look @avhz ! Agreed with your comment, all references to Monte-Carlo (method names, file names and documentation) have been amended 👍

Tested changes locally.

[avhz]

thanks for your patience !