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

[yfnaji]

This Pull Request was created to address issue 244.

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 and jump diffusion models
• Bring Heston model to life in order to test monte carlo methods on stochastic volatility models
• Implement Longstaff-Schwartz to calculate American options

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 ?