Finite Difference Pricer using heap allocated vectors

[yfnaji]

Important Notes

This Pull Request is regarding issue https://github.com/avhz/RustQuant/issues/98 for the creation of a finite difference pricer.

I have made 2 PRs for this particular issue. Only one of them needs to be merged if the decision has been made to merge.

nalgebra vs heap allocated vectors

This PR uses heap allocated vectors to represent vectors and matrices with functions to handle the appropriate linear algebra calculations (+ matrix multiplication has been optimised for tridiagonal matrices). The other PR utilises the nalgebra package.

The way in which the objects and methods should be used are exactly the same for both PRs.

There are differences in speed performance between the two implementations.

An example with the following parameters:

initial_price: 3.22,
strike_price: 12.87,
risk_free_rate: 0.02,
volatility: 0.2,
evaluation_date: None,
expiration_date: today() + Duration::days(365),
time_steps: 1000,
price_steps: 100,
type_flag: TypeFlag::Put,
exercise_flag: ExerciseFlag::American

We get the following runtimes when we run one instance of both implementations

Method	vec (latest commit)	nalgebra branch
Explicit	14.912845ms	447.828518ms
Implicit	212.050909ms	548.617466ms
Crank-Nicolson	242.015886ms	1.034193406s

We can immediately see that the explicit method is about ~17 22~ 24 times faster using vectors (this PR) than using the nalgebra package ~and runtime for the Crank-Nicolson method is roughly halved~ and runtime for the Crank-Nicolson method is 4 times faster and implicit method is roughly halved (This is at least the case for the machine I am running the code on).

Contents of this PR

This PR implements the following features:

• A struct, FiniteDifferencePricer in the module instruments::options::finite_difference_pricer, where the attributes are the parameters of an option, namely:

  • initial_price: f64
  • strike_price: f64
  • risk_free_rate: f64
  • volatility: f64
  • evaluation_date: Option<Date>,
  • expiration_date: Date,
  • time_steps: u32,
  • price_steps: u32,
  • type_flag: TypeFlag,
  • exercise_flag: ExerciseFlag

  where TypeFlag and ExerciseFlag are imported from instruments::options::option

• A constructor for the above struct with the new() method which takes in the above defined attributes as arguments
• Three public methods implemented in the FiniteDifferencePricer for the three widely known finite difference methods, namely explicit(), implicit() and crank_nicolson() which returns the option price rounded to 2 decimal places

Demonstration by example

We can start off by defining the finite difference object by calling the new() constructor

use RustQuant::instruments::options::{
    option::{TypeFlag, ExerciseFlag},
    finite_difference_pricer::FiniteDifferencePricer,
};

let finite_difference_obj = FiniteDifferencePricer::new(
    15.55,                            // initial_price
    10.2,                             // strike_price
    0.1,                              // risk_free_rate
    0.025,                            // volatility
    None,                             // evaluation_date
    today() + Duration::days(1),      // expiration_date
    1000,                             // time_steps
    100,                              // price_steps
    TypeFlag::Call,                   // type_flag
    ExerciseFlag::American            // exercise_flag
);

Now that the object has been defined, we can run any of the 3 methods to provide a numerical approximation of the option price:

println!("Explicit method: {}", finite_difference_obj.explicit());
println!("Implicit method: {}", finite_difference_obj.implicit());
println!("Crank-Nicolson: {}", finite_difference_obj.crank_nicolson());

Running the above code yields the following results:

Explicit method: 5.35
Implicit method: 5.35
Crank-Nicolson: 5.35

Notes

1. Since I am using heap allocated vectors in this PR, I have had to manually calculate the inverse for the tridiagonal matrix for the implicit and Crank-Nicolson finite difference methods. I have done so by implementing the theorem outlined in this paper.

2. I have opted to solve this problem using a struct to define the option parameters and from there we can run each finite difference method from the object. Alternatively, I could amend the PR so that it can be purely functional.

3. The code utilises the Dirichlet condition for the boundaries of the stock price extremities.

If we define the following:

$$V^{m}_{n} - \text{The option price at time step } n \text{ and price step } m\text{ (with maximum } M\text{)}$$

$$T - \text{Maturity time}$$

$$K - \text{Strike price}$$

$$S_{M} - \text{Maximum price of the underlying asset}$$

$$\Delta t - \text{Time step size}$$

Then Dirichlet boundary conditions for the call options are:

$$V^{0}_{n} = 0$$

$$V^{M}_{n} = S_M - Ke^{-r(T-n\Delta t)}$$

and for the put options:

$$V^{0}_{n} =Ke^{-r(T-n\Delta t)}$$

$$V^{M}_{n} = 0$$

4. We can also implement the option use the Neumann conditions, perhaps in another pull request in the future.

~EDIT 1: Optimisation amendments~

~In the code there are two private methods that seemingly do the same thing, create_tridiagonal_matrix() and create_trimmed_tridiagonal_matrix(). The difference between the two is that the former returns the matrx representation whole i.e. inclusive of all outer zeros, whilst the former only stores the the values of the three diagonals on a row by row basis.~

~As such, two private methods for matrix multiplication had to be created to handle the two data structures, appropriately named matrix_multiply_vector() and trimmed_tridiagonal_matrix_multiply_vector().~

~The reason both representation exist is because for explicit methods, only the diagonal values of the matrix are required to implement the approximation, whereas for implicit methods we require the whole matrix in order to calculate its inverse (though, this can be amended in a future PR).~

~This change has fractionally increased speed performances for explicit() and crank_nicolson().~

EDIT 2: Condensing methods + single data structure for tridiagonal matrices

Previously two distinct representations of a tridiagonal matrix were implemented:

• Data structure 1: A vector containing vectors representing rows in the matrix
• Data structure 2: A vector containing vectors of only the tridiagonal elements in a row

An example for each case respectively:

Data structure 1:

[
  [2,1,0,0,0],
  [1,2,1,0,0],
  [0,1,2,1,0],
  [0,0,1,2,1],
  [0,0,0,1,2],
]

Data structure 2:

[
  [2,1],
  [1,2,1],
  [1,2,1],
  [1,2,1],
  [1,2],
]

This branch will now only utilise data structure 2 to represent tridiagonal matrices and operations.

[yfnaji]

Changes have been made as per suggestions from @avhz to this branch and the nalgebra branch

[avhz]

thanks for your work !