F#: Capture real time bars from InteractiveBrokers

Let's continue our journey of capturing data from InteractiveBrokers. In this blog, we discuss how to capture real time bars. A bar here means a summary of prices over the last 5-second interval. Currently you can choose one of the following 4 types of prices that occurred in the interval to generate bars:

• TRADES: actual trade prices;
• BID: bid prices;
• ASK: ask prices;
• MIDPOINT: mid prices.

For the type of prices you chose, the information in a bar includes OHLC (open-high-low-close), volume, WAP (weighted average price), and count (the number of trades that occurred, which is applicable to TRADES only). Compared to ad hoc data like tick prices or market depth, real time bars are sent over to you on a regular basis -- you receive an update every 5 seconds.
The code example here is to capture real time bars on TRADES. As in the previous two blog posts, the first thing we need to do is to specify the stock we are interested in. Again, here I chose Infosys listed in Mumbai.

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "INFY"
        contract.exchange <- "NSE"
        contract.currency <- "INR"

Then, we use the reqRealTimeBarsEx method to subscribe to real time bar events. This method takes 5 parameters:

• ticker ID and contract: the same as in the previous two blog posts;
• bar size: the time length of each bar. Currently it must be 5, as only 5-second bars are supported;
• price type: TRADES, BID, ASK, or MIDPOINT; 
• useRTH: if a bar's interval fully or partially falls outside RTH (regular trading hours), this flag controls whether data outside RTH should be included in the bar. 0 means to exclude data outside RTH, and 1 means otherwise.

The following line instructs the broker that we wish to listen for real time bar events for the stock Infosys:

        tws1.reqRealTimeBarsEx(1, contract, 5, "TRADES", 0)

As for to unsubscribe, we just need invoke the cancelRealTimeBars method with the ticker ID for which we want to cancel subscription.

        tws1.cancelRealTimeBars(1)

Finally, below is a simple demo program, which captures real time bar events for 10 minutes and simply prints them out.

---------------------------------------------------------------------------------------
open AxTWSLib
open TWSLib
open System
open System.Drawing
open System.Windows.Forms

type RealTimeBarEvent =
    AxTWSLib._DTwsEvents_realtimeBarEvent

type RealTimeBarEventHandler =
    AxTWSLib._DTwsEvents_realtimeBarEventHandler

let now() = DateTime.Now

let formatTimeString (t:TimeSpan) =
    System.String.Format("{0:D2}:{1:D2}:{2:D2}",
            t.Hours, t.Minutes, t.Seconds)

let processRealTimeBarEvent _ (e:RealTimeBarEvent) =
    // ticker ID
    printf "id=%A," e.tickerId
    // time: the start of the bar interval
    let t = TimeSpan.FromSeconds(float e.time)
    printf "time=%s," (formatTimeString t)
    // open
    printf "o=%A," e.``open``
    // high
    printf "h=%A," e.high
    // low
    printf "l=%A," e.low
    // close
    printf "c=%A," e.close
    // wap
    printf "WAP=%A," e.wAP
    // volume
    printf "v=%A," e.volume
    // count
    printfn "cnt=%A" e.count

[<EntryPoint; STAThread>]
let main _ =
    printfn "start time: %A" (now())
    let form1 = new Form(Text="Dummy Form")
    let tws1 = new AxTws()
    tws1.BeginInit()
    form1.Controls.Add(tws1)
    tws1.EndInit()

    tws1.connect("", 7496, 1)
    printfn "server version = %d" tws1.serverVersion

    if tws1.serverVersion > 0 then
        tws1.realtimeBar.AddHandler(
            new RealTimeBarEventHandler(processRealTimeBarEvent))

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "INFY"
        contract.exchange <- "NSE"
        contract.currency <- "INR"
        tws1.reqRealTimeBarsEx(1, contract, 5, "TRADES", 0)

        let timerWorkflow = async {
            do! Async.Sleep (60000 * 10)
            Application.Exit()
            }
        Async.Start timerWorkflow

        Application.Run()
        tws1.cancelRealTimeBars(1)
        tws1.disconnect()
        printfn "Diconnected!"

    printfn "end time: %A" (now())
    0

F#: Capture market depth data from InteractiveBrokers

In the last blog post, we discussed how to capture real-time tick data from InteractiveBrokers. Now let's see how we can capture real-time market depth data from InteractiveBrokers. The example code in this post is very similar to the one in the previous post. As far as the code is concerned, the main difference is just that now we subscribe to and process market depth events instead of tick price events. However, to interpret market depth events is slightly more complicated than to interpret tick data events, since an order book is basically a table with multiple entries.

Before we talk about how to process market depth events to maintain an order book, let's first take a look at an example below, which is a snapshot of the order book for the stock "Infosys" listed in Mumbai. The snapshot was taken at some time on 24th Feb 2014. You can see that the depth of this order book is 5, i.e., it shows the best 5 entries for both bid side and ask side, so there are 10 entries in total. Moreover, all entries are sorted by price, so that the entries for bid are in descending order and those for ask are in ascending order. By observing the market depth data, it's not difficult to see that the first bid/ask entry of an order book should be always the same as the tick price data. As such, subscribing to tick price events is equivalent to subscribing to market depth events with depth 1.

	bid		ask
position	price	size	price	size
1	3748.70	14	3749.00	71
2	3748.50	290	3749.40	10
3	3748.45	20	3749.45	1
4	3748.15	25	3749.85	50
5	3748.10	125	3749.90	15

So, how do we subscribe to market depth events? It's almost as easy as to subscribe to tick price events. The following code snippet creates an contract object specifying the stock we are interested in. In this case, it's Infosys listed in Mumbai.

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "INFY"
        contract.exchange <- "NSE"
        contract.currency <- "INR"

Then, we use the reqMktDepthEx method, which takes 3 parameters: ticker ID, contract, and market depth. The following line instructs the broker that we wish to listen for market depth events up to depth 5.

        tws1.reqMktDepthEx(1, contract, 5)

Please note that though one can specify an arbitrarily large number for market depth, the effective market depth will still be subject to the exchange and your data subscription with InteractiveBrokers. More often than not, market depth data are not free, but there are some exceptions, like NSE in Mumbai. Once we have invoked the reqMktDepthEx method, we immediately receive the initial batch of market depth events, which are meant to initialize a complete snapshot of the order book. After the initial batch, we will receive subsequent market depth events on an ad hoc basis. These subsequent events represent increment changes to the order book. Basically, a market depth event consists of the following properties:

• id: ticker ID, which is specified by us via the reqMktDepthEx method;
• side: 0 = ask, 1 = bid;
• operation: 0 = insert, 1 = update, 2 = delete;
• position: indicate which position in the order book to update;
• price: bid/ask price;
• size: bid/ask  size.

For example, the snapshot example above was constructed by the following initial batch of events:

event sequence	id	side	operation	position	price	size
1	1	1	0	0	3748.7	14
2	1	1	0	1	3748.5	290
3	1	1	0	2	3748.45	20
4	1	1	0	3	3748.15	25
5	1	1	0	4	3748.1	125
6	1	0	0	0	3749	71
7	1	0	0	1	3749.4	10
8	1	0	0	2	3749.45	1
9	1	0	0	3	3749.85	50
10	1	0	0	4	3749.9	15

From the 11th event onward, we receive increment changes. For example, the 11th event was as follows:

event sequence	id	side	operation	position	price	size
11	1	1	1	0	3748.7	17

   

Based on this event, we should update the first bid entry of the order book with price 3748.7 and size 17. So I hope this is clear enough on how to interpret market depth events.

As for to unsubscribe, we just need invoke the cancelMktDepth method with the ticker ID for which we want to cancel the subscription.

        tws1.cancelMktDepth(1)

Finally, below is a simple demo program, which captures market depth events for 10 minutes and simply prints them out. This program does not present an order book visually.

----------------------------------------------------------------------------------------------
open AxTWSLib
open TWSLib
open System
open System.Drawing
open System.Windows.Forms

type MktDepthEvent =
    AxTWSLib._DTwsEvents_updateMktDepthEvent

type MktDepthEventHandler =
    AxTWSLib._DTwsEvents_updateMktDepthEventHandler

let now() = DateTime.Now

let processMktDepthEvent _ (e:MktDepthEvent) =
    // event timestamp 
    printf "%A," (now())
    // event ID, which should match the subscription ID
    printf "id=%A," e.id
    // side, 0=ask, 1=buy
    printf "side=%A," e.side
    // operation, 0=insert, 1=update, 2=delete
    printf "op=%A," e.operation
    // position
    printf "pos=%A," e.position
    // price
    printf "price=%A," e.price
    // size
    printfn "size=%A" e.size

[<EntryPoint; STAThread>]
let main _ = 
    printfn "start time: %A" (now())
    let form1 = new Form(Text="Dummy Form")
    let tws1 = new AxTws()
    tws1.BeginInit()
    form1.Controls.Add(tws1)
    tws1.EndInit()

    tws1.connect("", 7496, 1)
    printfn "server version = %d" tws1.serverVersion

    if tws1.serverVersion > 0 then
        tws1.updateMktDepth.AddHandler(
            new MktDepthEventHandler(processMktDepthEvent))

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "INFY"
        contract.exchange <- "NSE"
        contract.currency <- "INR"
        tws1.reqMktDepthEx(1, contract, 5)

        let timerWorkflow = async { 
            do! Async.Sleep 10000
            Application.Exit() 
            }
        Async.Start timerWorkflow

        Application.Run()
        tws1.cancelMktDepth(1)
        tws1.disconnect()
        printfn "Diconnected!"

    printfn "end time: %A" (now())
    0

   

F#: Capture tick price data from InteractiveBrokers

In the last blog post, we saw how to set up connectivity to InteractiveBrokers. Let's move on to see how we can capture real-time tick price data from InteractiveBrokers. Before I show the complete example code, let me first go over some key steps.
The first thing is to specify for which asset you want to capture tick prices. We do that by filling in an IContract object. In the code snippet below, tws1 represents the connection to InteractiveBrokers and we use the factory method, createContract(), to create a blank IContract object. If the asset we are interested in is S&P 500 ETF Trust, we can fill in the IContract object as follows:

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "SPY"
        contract.exchange <- "SMART"
        contract.currency <- "USD"

Then we can submit the IContract object via the reqMktDataEx method so as to subscribe to prices information of S&P 500 ETF Trust. The reqMktDataEx method takes 4 parameters. Here we need to focus on only the first 2: ticker ID and contract. Ticker ID is an integer, which we can choose any number as long as we assign a unique number to each asset we subscribe to. Below is how we invoke the reqMktDataEx method to submit the IContract object we created.

        tws1.reqMktDataEx(1, contract, "", 0)

So now we let ticker ID 1 represent S&P 500 ETF Trust. If later on we no longer want to listen to the prices of S&P 500 ETF Trust, we can pass ticker ID 1 to the cancelMktData method to stop receiving any further updates.

        tws1.cancelMktData(1)

After subscribing via the reqMktDataEx method, any tick price updates will be sent to us via events. As such, we write the event handler, processTickPriceEvent, to process tick price events. The handler is simply just to distinguish different types of tick prices and print them out. (As for the code of the handler, please refer to the complete code at the end of this blog post.) Then we register the handler with tick price events as follows:

        tws1.tickPrice.AddHandler(
            new TickPriceEventHandler(processTickPriceEvent))

After registering with the events, we have to invoke Application.Run() to enter the event loop to start capturing tick prices. My intention here is to capture tick prices for 10 seconds and then quit. To end the event loop after 10 seconds, I set up the following asynchronous workflow.

        let timerWorkflow = async { 
            do! Async.Sleep 10000
            Application.Exit() 
            }
        Async.Start timerWorkflow 

OK, now that we have covered all the major steps, below is the complete code of the program.

/////////////////////////////////////////////////////////////////////////////////////////////////////////
open AxTWSLib
open TWSLib
open System
open System.Drawing
open System.Windows.Forms

type TickPriceEvent = 
    AxTWSLib._DTwsEvents_tickPriceEvent

type TickPriceEventHandler =
    AxTWSLib._DTwsEvents_tickPriceEventHandler

let now() = DateTime.Now

let processTickPriceEvent _ (e:TickPriceEvent) =
    match e.tickType with
    | 1 -> printfn "%A|bid=%A" (now()) e.price
    | 2 -> printfn "%A|ask=%A" (now()) e.price
    | 4 -> printfn "%A|last=%A" (now()) e.price
    | 6 -> printfn "%A|high=%A" (now()) e.price
    | 7 -> printfn "%A|low=%A" (now()) e.price
    | 9 -> printfn "%A|close=%A" (now()) e.price
    | x -> printfn "%A|UNKNOW_TICK_TYPE(%d)=%A" (now()) x e.price

[<EntryPoint; STAThread>]
let main _ = 
    printfn "start time: %A" (now())
    let form1 = new Form(Text="Dummy Form")
    let tws1 = new AxTws()
    tws1.BeginInit()
    form1.Controls.Add(tws1)
    tws1.EndInit()

    tws1.connect("", 7496, 1)
    printfn "server version = %d" tws1.serverVersion

    if tws1.serverVersion > 0 then
        tws1.tickPrice.AddHandler(
            new TickPriceEventHandler(processTickPriceEvent))

        let contract = tws1.createContract()
        contract.secType <- "STK"
        contract.symbol <- "SPY"
        contract.exchange <- "SMART"
        contract.currency <- "USD"
        tws1.reqMktDataEx(1, contract, "", 0)

        let timerWorkflow = async { 
            do! Async.Sleep 10000
            Application.Exit() 
            }
        Async.Start timerWorkflow

        Application.Run()
        tws1.cancelMktData(1)
        tws1.disconnect()
        printfn "Diconnected!"

    printfn "end time: %A" (now())

    0

F#: Connect to InteractiveBrokers via ActiveX

In this blog post we show how to write a F# console program which runs on Windows and connect to InteractiveBrokers (IB). There are a few different options of programmatic connectivity provided by IB. Here we use the ActiveX control (Tws.ocx) and Trade Workstation (TWS), so the end-to-end communication between our console program and IB goes like this:

F# code <=> Tws.ocx <=> TWS <=> IB server.

Before we begin, you need to have an account with IB, which can be either a real account or a paper account. Once you have the account ready, here is the step-by-step guide to the console program.

Install TWS and API:

• Trade Workstation: please download and install a standalone version from here. You need to configure TWS to accept ActiveX clients by doing the following:
• Select Edit->Global Configuration
• Select API->Settings
• Check "Enable ActiveX and Socket Clients" and note down the port number specified at "Socket Port". By default, the port number is 7496.
• API: can be downloaded from here, which contains Tws.ocx.

Generate wrapper DLLs to expose Tws.ocx to .NET:
Assuming that the API has been installed at C:\TWS_API and you have Visual Studio loaded on your machine, you can then run Aximp.exe as follows:

C:\TWS_API\bin\ActiveX>aximp Tws.ocx

And you shall see the following 2 DLLs generated:

Generated Assembly: C:\TWS_API\bin\ActiveX\TWSLib.dll
Generated Assembly: C:\TWS_API\bin\ActiveX\AxTWSLib.dll

To make use of the ActiveX control, we need to add references to these 2 DLLs in the Visual Studio project. (Note that as the only thing we are trying to do here is connecting to and disconnecting from TWS, we need only AxTWSLib.dll. However, in following few blog posts, we will try to capture market data and submit trade orders, which will requires both DLLs.)

Write F# code to connect and disconnect:
Below you can see the code of our console program. There are some remarks about the code I would like to make:

• Tws.ocx is an ActiveX control which is meant to be added to a GUI container. To get the control initialized properly, a dummy Windows Form is created to host it, as you can see at the beginning of the main function. As our intention here is a console program, we don't display the form.
• To connect to TWS, we invoke the connect method, which takes 3 parameters: IP address, port number, client ID. 
• IP address: to indicate where is the machine running TWS. Since I run TWS and the program on the same machine, there's no need to provide an IP address.
• port number: 7496 by default, unless you change it in the Global Configuration page of TWS.
• client ID: an integer used to identify this client connection. If you want to connect multiple clients to the same TWS instance at the same time, each client should be assigned a unique ID. 
• Once TWS receives a connection request, a dialog pops up. Click on yes to accept the connection. 
• To see whether the connection is successful, we check the serverVersion property of the control. If it's successful, we shall a positive serverVersion. Otherwise, we receive a zero. 

------------------------------------------------------------------------
open AxTWSLib
open System
open System.Drawing
open System.Windows.Forms

[<EntryPoint; STAThread>]
let main _ = 
    // initialize TWS ActiveX Control tws1
    let form1 = new Form(Text="Dummy Form")
    let tws1 = new AxTws()
    tws1.BeginInit()
    form1.Controls.Add(tws1)
    tws1.EndInit()

    // connect to local TWS
    tws1.connect("", 7496, 1)
    printfn "server version = %d" tws1.serverVersion

    // disconnect from TWS
    if tws1.serverVersion > 0 then 
        tws1.disconnect()

    0
------------------------------------------------------------------------

F#: contributing a chapter to an upcoming book, "F# Deep Dives"

Thanks to Tomas Petricek's invitation, I am honored by being able to write a chapter for an upcoming book - F# Deep Dives. The chapter I am responsible for is Chapter 4: Numerical computing in financial domain. The first draft of the chapter is now available online via Manning's early access program.

Currently the chapter covers the following:

1 Introducing financial derivatives

2 Using probability functions of Math.NET

2.1 Configuring F# Interactive (profiling and floating point formatting)

2.2 Setting up Math.NET Numerics

2.3 Introducing Random Variables, Expectation and Variance

2.4 Generating Normally Distributed Samples

3 Geometric Brownian Motion and Monte Carlo Estimate

3.1 Modeling stock prices using geometric Brownian motion

3.2 Payoff Function, Discounted Payoff, and Monte Carlo estimate

3.3 Analyzing variance of Monte Carlo estimates

3.4 Pricing Path-dependent options (Asian options and barrier options)

3.5 Reducing Variance by antithetic variates

F#: K-S test on final prices of GBM paths

In the previous post, I mentioned that I chose System.Random, because it’s faster than any other random number generators available in Math.NET Numerics. And I also said that I have a minor concern about System.Random, because someone found a problem in the implementation of System.Random and raised this issue to the .NET team in early 2011. And the team admits it’s a problem. However, this issue is still there with .NET Framework 4.5. And it seems to me that the .NET team is not planning to take any immediate action. 

As I know little about the algorithm of System.Random, I’m not sure if it’s worth the efforts to pick up the details and quantify the impact. As such, I decided to take an easy way out, i.e., close my eyes to this issue  do Kolmogrov-Smirnov test (KS test) on only the final prices distribution. If the simulated distribution is statistically indistinguishable from the true, theoretic distribution in terms of shape, I will be fine with System.Random.  At the end of the day, I think what really matters to me is whether System.Random can ultimately give me the correct shape of the final prices distribution I want.

Basically the F# code is pretty much the same as in the previous post except the following:

• I add a run_stat_test function to run the statistical test and generate a Q-Q plot in Matlab.


• I change the Monte Carlo parameter N to 1 so that I simulate only one data point, i.e., final price, for each simulated path.

// Monte Carlo parameters
let T = 1.0             // time horizon in years
let M = 100000          // number of simulations
let N = 1               // number of time intervals

You can see that I take the logarithms of all final prices and test them against a normal distribution with mean = (log S0)+(r-0.5*sigma*sigma)*T and standard deviation = sigma*sqrt(T) using KS test. And the output from the test is as follows:

h = false, p = 0.7243765247

And here is the Q-Q plot:



p-value is about 0.72, which means it’s quite likely that the simulated prices were indeed sampled from the normal distribution, and h=false means that we can’t reject the null hypothesis. Moreover, the Q-Q plot does not show any obvious outliers off the diagonal line. Therefore, it’s ok for me to continue to use System.Random.

If any expert happens to know about the issue of System.Random and is willing to shed some light on it, it will be appreciated.

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

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 = 1               // 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)) )

let run_stat_test (prices:float array) =
    let m = (log S0)+(r-0.5*sigma*sigma)*T
    let s = sigma*sqrt(T)
    let log_prices = Array.map log prices
    let matlab = new MLApp.MLAppClass()
    matlab.PutWorkspaceData("log_prices", "base", log_prices)
    matlab.PutWorkspaceData("m", "base", m)
    matlab.PutWorkspaceData("s", "base", s)
    matlab.Execute("pd = ProbDistUnivParam('normal',[m s])") |> ignore
    matlab.Execute("[h,p] = kstest(log_prices,pd,0.05,0)") |> ignore
    matlab.Execute("qqplot(log_prices,pd)") |> ignore
    let h = matlab.GetVariable("h","base")
    let p = matlab.GetVariable("p","base")
    (h,p)

// 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)
    let h,p = run_stat_test final_prices
    printfn "h = %A, p = %A" h p
    (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))

F#: Simulate entire GBM path

In all the Monte Carlo simulation codes I have done in previous posts, I simulate only the underlying prices at expiry, because European Calls and Puts are dependent on only the final price at expiry and we knew that the final prices follow a lognormal distribution, for which there's a nice analytic expression and hence we can easily simulate. However, if we move further out of the safe zone of European options, it's quite often that we need to simulate entire GBM (Geometric Brownian Motion) paths. Before we really get into dealing with path-dependent options, let's first take a look at my implementation for generating GBM paths based on the Euler method.

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:


// 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, 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))