In this post, I started using a open-source numerical library, Math.NET Numerics. Don Syme wrote a blog on how to install and use it. Basically I use it in my code for two purposes:

- using the Normal class to draw normal random numbers. Internally the Normal class also uses the Box Muller transform to generate samples, but it uses the polar form of the transform, which is faster than the basic form I had been using in previous posts. Another good thing about the polar form is that it does not suffer the infinity issue I mentioned in a previous post.
- using some handy statistical functions, like mean, standard deviation, and inverse CDF, because I like to compute standard error and confidence interval.

Before we look into the code, here are the parameters I use:

Below is the code snippet to construct a function which generates Brownian Motion increments. You can see that I use System.Random as the random source, because it's a few times faster than all other uniform random sources provided by Math.NET Numerics. (I actually have a minor concern about System.Random. I might talk about it later in another post.)

**// Monte Carlo parameters****let T = 1.0 // time horizon in years****let M = 100000 // number of simulations****let N = 252 // number of time intervals****let p = 0.99 // confidence level**

**// underlying dynamics****let S0 = 50.0 // spot price****let sigma = 0.2 // annualized vol****let r = 0.01 // risk-free rate**Below is the code snippet to construct a function which generates Brownian Motion increments. You can see that I use System.Random as the random source, because it's a few times faster than all other uniform random sources provided by Math.NET Numerics. (I actually have a minor concern about System.Random. I might talk about it later in another post.)

**// risk factor dW**

**let rnd = new System.Random(1)**

**let get_dW dt =**

**let dW = Normal.WithMeanVariance(0.0, dt)**

**dW.RandomSource <- rnd**

**(fun () -> dW.Sample())**

The below function,

As you can see from the main function of the code, after generation of GBM paths, I extracted final prices out of all the paths and passed them to:

Finally, here are the results from the main function:

The entries in the above tuple are mean, standard deviation, standard error, and confidence interval. I also calculated theoretical mean and standard deviation as follows:

And below is the KDE plot:

///////////////////////////////////////////////////////////////////////////////////////////////////////

**generate_GBM_paths**, is the one that actually simulate all the entire paths. Here I use an array of arrays to keep all the paths simulated. Storing every complete simulated path is memory consuming, but it's convenient for me to do post-simulation analysis. I think this code snippet is a good example to show that F#'s built-in array processing functions are really expressive and handy, compared to other major general-purpose languages.**// GBM path generation based on the Euler method****let generate_GBM_paths S0 r sigma T M N =****let dt = T / (float N)****let drift = (r - 0.5 * (sigma**2.0)) * dt****let dW = get_dW dt****Array.init M (fun _ ->****Array.init N (fun _ -> drift + sigma * dW()) |>****Array.scan (+) 0.0 |>****Array.map (fun x -> S0 * exp(x)) )**As you can see from the main function of the code, after generation of GBM paths, I extracted final prices out of all the paths and passed them to:

- the statistical functions provided by Math.NET Numerics to compute
- mean, i.e., the forward price of the underlying
- standard deviation of the simulated final price distribution
- standard error of the final price simulation
- confidence interval
- Matlab to plot KDE (Kernel Density Estimation) to see if the simulated final price distribution looks like a lognormal distribution. I recently just learnt KDE from the book, "Data Analysis with Open Source Tools," written by Philipp K. Janert. In case you haven't heard about KDE, it's something like an advanced version of histograms. In a KDE, the graph will be smoothed by a density function and you don't need to find out what bin size you should use. If you don't like to run this KDE part, you could simply comment it out.

Finally, here are the results from the main function:

**Real: 00:00:07.541, CPU: 00:00:06.817, GC gen0: 201, gen1: 135, gen2: 2****val it : float * float * float * (float * float) =****(50.47657684, 10.21848436, 0.03231368481, (50.3933423, 50.55981138))**The entries in the above tuple are mean, standard deviation, standard error, and confidence interval. I also calculated theoretical mean and standard deviation as follows:

**val true_forward : float = 50.50250835****val true_std_dev : float = 10.20235347**And below is the KDE plot:

///////////////////////////////////////////////////////////////////////////////////////////////////////

**open MathNet.Numerics.Distributions**

**open MathNet.Numerics.Statistics**

**// Monte Carlo parameters**

**let T = 1.0 // time horizon in years**

**let M = 100000 // number of simulations**

**let N = 252 // number of time intervals**

**let p = 0.99 // confidence level**

**// underlying dynamics**

**let S0 = 50.0 // spot price**

**let sigma = 0.2 // annualized vol**

**let r = 0.01 // risk-free rate**

**// risk factor dW**

**let rnd = new System.Random(1)**

**let get_dW dt =**

**let dW = Normal.WithMeanVariance(0.0, dt)**

**dW.RandomSource <- rnd**

**(fun () -> dW.Sample())**

**// GBM path generation**

**let generate_GBM_paths S0 r sigma T M N =**

**let dt = T / (float N)**

**let drift = (r - 0.5 * (sigma**2.0)) * dt**

**let dW = get_dW dt**

**Array.init M (fun _ ->**

**Array.init N (fun _ -> drift + sigma * dW()) |>**

**Array.scan (+) 0.0 |>**

**Array.map (fun x -> S0 * exp(x)) )**

**// plot the simulated distribution using KDE in Matlab**

**let plot_KDE (prices:float array) =**

**let matlab = new MLApp.MLAppClass()**

**matlab.PutWorkspaceData("x", "base", prices)**

**matlab.Execute("[f,xi] = ksdensity(x);") |> ignore**

**matlab.Execute("plot(xi,f);") |> ignore**

**()**

**// main**

**let main () =**

**let final_prices = generate_GBM_paths S0 r sigma T M N |>**

**Array.map (fun p -> p.[p.Length - 1])**

**let forward = final_prices.Mean()**

**let std_dev = final_prices.StandardDeviation()**

**let std_error = std_dev / sqrt(float M)**

**let std_normal = Normal.WithMeanVariance(0.0,1.0)**

**let radius = std_normal.InverseCumulativeDistribution(0.5*(1.0+p))**

*** std_error**

**let conf_int = (forward - radius, forward + radius)**

**plot_KDE final_prices**

**(forward, std_dev, std_error, conf_int)**

**// run simulation**

**main()**

**// calculate theoretical mean and deviation**

**let true_forward = S0 * exp(r*T)**

**let v = sigma*sigma*T**

**let true_std_dev = sqrt((exp(v)-1.0)*exp(2.0*(log(S0)+r*T-0.5*v*T)+v))**

## No comments:

## Post a Comment