tag:blogger.com,1999:blog-50632076327966001852024-02-19T14:31:27.574-08:00Programming CradleChao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.comBlogger15125tag:blogger.com,1999:blog-5063207632796600185.post-57704848633169514032014-03-26T00:24:00.000-07:002014-03-26T00:41:17.826-07:00F#: Capture real time bars from InteractiveBrokersLet'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:<br />
<br />
<ul>
<li>TRADES: actual trade prices;</li>
<li>BID: bid prices;</li>
<li>ASK: ask prices;</li>
<li>MIDPOINT: mid prices.</li>
</ul>
<br />
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.<br />
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.<br />
<br />
let contract = tws1.createContract()<br />
contract.secType <- "STK"<br />
contract.symbol <- "INFY"<br />
contract.exchange <- "NSE"<br />
contract.currency <- "INR"<br />
<br />
Then, we use the reqRealTimeBarsEx method to subscribe to real time bar events. This method takes 5 parameters:<br />
<br />
<ul>
<li>ticker ID and contract: the same as in the previous two blog posts;</li>
<li>bar size: the time length of each bar. Currently it must be 5, as only 5-second bars are supported;</li>
<li>price type: TRADES, BID, ASK, or MIDPOINT; </li>
<li>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.</li>
</ul>
<br />
The following line instructs the broker that we wish to listen for real time bar events for the stock Infosys:<br />
<br />
tws1.reqRealTimeBarsEx(1, contract, 5, "TRADES", 0)<br />
<br />
As for to unsubscribe, we just need invoke the cancelRealTimeBars method with the ticker ID for which we want to cancel subscription.<br />
<br />
tws1.cancelRealTimeBars(1)<br />
<br />
Finally, below is a simple demo program, which captures real time bar events for 10 minutes and simply prints them out.<br />
<br />
---------------------------------------------------------------------------------------<br />
open AxTWSLib<br />
open TWSLib<br />
open System<br />
open System.Drawing<br />
open System.Windows.Forms<br />
<br />
type RealTimeBarEvent =<br />
AxTWSLib._DTwsEvents_realtimeBarEvent<br />
<br />
type RealTimeBarEventHandler =<br />
AxTWSLib._DTwsEvents_realtimeBarEventHandler<br />
<br />
let now() = DateTime.Now<br />
<br />
let formatTimeString (t:TimeSpan) =<br />
System.String.Format("{0:D2}:{1:D2}:{2:D2}",<br />
t.Hours, t.Minutes, t.Seconds)<br />
<br />
let processRealTimeBarEvent _ (e:RealTimeBarEvent) =<br />
// ticker ID<br />
printf "id=%A," e.tickerId<br />
// time: the start of the bar interval<br />
let t = TimeSpan.FromSeconds(float e.time)<br />
printf "time=%s," (formatTimeString t)<br />
// open<br />
printf "o=%A," e.``open``<br />
// high<br />
printf "h=%A," e.high<br />
// low<br />
printf "l=%A," e.low<br />
// close<br />
printf "c=%A," e.close<br />
// wap<br />
printf "WAP=%A," e.wAP<br />
// volume<br />
printf "v=%A," e.volume<br />
// count<br />
printfn "cnt=%A" e.count<br />
<br />
[<EntryPoint; STAThread>]<br />
let main _ =<br />
printfn "start time: %A" (now())<br />
let form1 = new Form(Text="Dummy Form")<br />
let tws1 = new AxTws()<br />
tws1.BeginInit()<br />
form1.Controls.Add(tws1)<br />
tws1.EndInit()<br />
<br />
tws1.connect("", 7496, 1)<br />
printfn "server version = %d" tws1.serverVersion<br />
<br />
if tws1.serverVersion > 0 then<br />
tws1.realtimeBar.AddHandler(<br />
new RealTimeBarEventHandler(processRealTimeBarEvent))<br />
<br />
let contract = tws1.createContract()<br />
contract.secType <- "STK"<br />
contract.symbol <- "INFY"<br />
contract.exchange <- "NSE"<br />
contract.currency <- "INR"<br />
tws1.reqRealTimeBarsEx(1, contract, 5, "TRADES", 0)<br />
<br />
let timerWorkflow = async {<br />
do! Async.Sleep (60000 * 10)<br />
Application.Exit()<br />
}<br />
Async.Start timerWorkflow<br />
<br />
Application.Run()<br />
tws1.cancelRealTimeBars(1)<br />
tws1.disconnect()<br />
printfn "Diconnected!"<br />
<br />
printfn "end time: %A" (now())<br />
0<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com1tag:blogger.com,1999:blog-5063207632796600185.post-2929563072707748152014-03-04T16:42:00.000-08:002014-03-04T16:43:22.260-08:00F#: Capture market depth data from InteractiveBrokersIn the <a href="http://programmingcradle.blogspot.sg/2014/01/f-capture-tick-price-data-from.html" target="_blank">last blog post</a>, 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.<br />
<br />
Before we talk about how to process market depth events to maintain an <a href="http://en.wikipedia.org/wiki/Order_book_(trading)" target="_blank">order book</a>, 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.<br />
<br />
<table border="1" cellpadding="0" cellspacing="0" style="border-collapse: collapse; text-align: center; width: 320px;">
<colgroup><col span="5" style="width: 48pt;" width="64"></col>
</colgroup><tbody>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15.0pt; width: 48pt;" width="64"></td>
<td class="xl64" colspan="2" style="width: 96pt;" width="128">bid</td>
<td class="xl64" colspan="2" style="width: 96pt;" width="128">ask</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15.0pt;">position</td>
<td class="xl64">price</td>
<td class="xl64">size</td>
<td class="xl64">price</td>
<td class="xl64">size</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td class="xl65">3748.70</td>
<td class="xl64">14</td>
<td class="xl65">3749.00</td>
<td class="xl64">71</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
2</div>
</td>
<td class="xl65">3748.50</td>
<td class="xl64">290</td>
<td class="xl65">3749.40</td>
<td class="xl64">10</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
3</div>
</td>
<td class="xl65">3748.45</td>
<td class="xl64">20</td>
<td class="xl65">3749.45</td>
<td class="xl64">1</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
4</div>
</td>
<td class="xl65">3748.15</td>
<td class="xl64">25</td>
<td class="xl65">3749.85</td>
<td class="xl64">50</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
5</div>
</td>
<td class="xl65">3748.10</td>
<td class="xl64">125</td>
<td class="xl65">3749.90</td>
<td class="xl64">15</td>
</tr>
</tbody></table>
<div style="text-align: center;">
<br /></div>
<div style="text-align: center;">
<div style="text-align: left;">
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.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> let contract = tws1.createContract()</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.secType <- "STK"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.symbol <- "INFY"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.exchange <- "NSE"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.currency <- "INR"</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
Then, we use the <span style="font-family: Verdana, sans-serif;">reqMktDepthEx</span> 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.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> tws1.reqMktDepthEx(1, contract, 5)</span><br />
<br /></div>
<div style="text-align: left;">
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 <span style="font-family: Verdana, sans-serif;">reqMktDepthEx</span> 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:</div>
<div style="text-align: left;">
</div>
<ul>
<li style="text-align: left;">id: ticker ID, which is specified by us via the <span style="font-family: Verdana, sans-serif;">reqMktDepthEx</span> method;</li>
<li style="text-align: left;">side: 0 = ask, 1 = bid;</li>
<li style="text-align: left;">operation: 0 = insert, 1 = update, 2 = delete;</li>
<li style="text-align: left;">position: indicate which position in the order book to update;</li>
<li style="text-align: left;">price: bid/ask price;</li>
<li style="text-align: left;">size: bid/ask size.</li>
</ul>
<br />
<div style="text-align: left;">
For example, the snapshot example above was constructed by the following initial batch of events:</div>
<div style="text-align: left;">
<br /></div>
<table border="1" cellpadding="0" cellspacing="0" style="border-collapse: collapse; width: 528px;">
<colgroup><col style="mso-width-alt: 5266; mso-width-source: userset; width: 108pt;" width="144"></col>
<col span="6" style="width: 48pt;" width="64"></col>
</colgroup><tbody>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15.0pt; width: 108pt;" width="144"><div style="text-align: center;">
event sequence</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
id</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
side</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
operation</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
position</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
price</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
size</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.7</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
14</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
2</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.5</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
290</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
3</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
2</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.45</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
20</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
4</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.15</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
25</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
5</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
4</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
125</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
6</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3749</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
71</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
7</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3749.4</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
10</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
8</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
2</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3749.45</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
9</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3749.85</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
50</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
10</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
4</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3749.9</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
15</div>
</td>
</tr>
</tbody></table>
<br />
<div style="text-align: left;">
From the 11th event onward, we receive increment changes. For example, the 11th event was as follows:<br />
<br />
<table border="1" cellpadding="0" cellspacing="0" style="border-collapse: collapse; width: 528px;">
<colgroup><col style="mso-width-alt: 5266; mso-width-source: userset; width: 108pt;" width="144"></col>
<col span="6" style="width: 48pt;" width="64"></col>
</colgroup><tbody>
<tr height="20" style="height: 15.0pt;">
<td height="20" style="height: 15.0pt; width: 108pt;" width="144"><div style="text-align: center;">
event sequence</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
id</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
side</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
operation</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
position</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
price</div>
</td>
<td style="width: 48pt;" width="64"><div style="text-align: center;">
size</div>
</td>
</tr>
<tr height="20" style="height: 15.0pt;">
<td class="xl65" height="20" style="height: 15pt; text-align: right;"><div style="text-align: center;">
11</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
1</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
0</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
3748.7</div>
</td>
<td style="text-align: right;"><div style="text-align: center;">
17</div>
</td>
</tr>
</tbody></table>
</div>
<div style="text-align: left;">
</div>
<div style="text-align: left;">
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.<br />
<br />
As for to unsubscribe, we just need invoke the <span style="font-family: Verdana, sans-serif;">cancelMktDepth</span> method with the ticker ID for which we want to cancel the subscription.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> tws1.cancelMktDepth(1)</span><br />
<br />
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.<br />
<br />
----------------------------------------------------------------------------------------------<br />
<span style="font-family: Verdana, sans-serif;">open AxTWSLib</span><br />
<span style="font-family: Verdana, sans-serif;">open TWSLib</span><br />
<span style="font-family: Verdana, sans-serif;">open System</span><br />
<span style="font-family: Verdana, sans-serif;">open System.Drawing</span><br />
<span style="font-family: Verdana, sans-serif;">open System.Windows.Forms</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">type MktDepthEvent =</span><br />
<span style="font-family: Verdana, sans-serif;"> AxTWSLib._DTwsEvents_updateMktDepthEvent</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">type MktDepthEventHandler =</span><br />
<span style="font-family: Verdana, sans-serif;"> AxTWSLib._DTwsEvents_updateMktDepthEventHandler</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">let now() = DateTime.Now</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">let processMktDepthEvent _ (e:MktDepthEvent) =</span><br />
<span style="font-family: Verdana, sans-serif;"> // event timestamp </span><br />
<span style="font-family: Verdana, sans-serif;"> printf "%A," (now())</span><br />
<span style="font-family: Verdana, sans-serif;"> // event ID, which should match the subscription ID</span><br />
<span style="font-family: Verdana, sans-serif;"> printf "id=%A," e.id</span><br />
<span style="font-family: Verdana, sans-serif;"> // side, 0=ask, 1=buy</span><br />
<span style="font-family: Verdana, sans-serif;"> printf "side=%A," e.side</span><br />
<span style="font-family: Verdana, sans-serif;"> // operation, 0=insert, 1=update, 2=delete</span><br />
<span style="font-family: Verdana, sans-serif;"> printf "op=%A," e.operation</span><br />
<span style="font-family: Verdana, sans-serif;"> // position</span><br />
<span style="font-family: Verdana, sans-serif;"> printf "pos=%A," e.position</span><br />
<span style="font-family: Verdana, sans-serif;"> // price</span><br />
<span style="font-family: Verdana, sans-serif;"> printf "price=%A," e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> // size</span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "size=%A" e.size</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">[<EntryPoint; STAThread>]</span><br />
<span style="font-family: Verdana, sans-serif;">let main _ = </span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "start time: %A" (now())</span><br />
<span style="font-family: Verdana, sans-serif;"> let form1 = new Form(Text="Dummy Form")</span><br />
<span style="font-family: Verdana, sans-serif;"> let tws1 = new AxTws()</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.BeginInit()</span><br />
<span style="font-family: Verdana, sans-serif;"> form1.Controls.Add(tws1)</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.EndInit()</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> tws1.connect("", 7496, 1)</span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "server version = %d" tws1.serverVersion</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> if tws1.serverVersion > 0 then</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.updateMktDepth.AddHandler(</span><br />
<span style="font-family: Verdana, sans-serif;"> new MktDepthEventHandler(processMktDepthEvent))</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> let contract = tws1.createContract()</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.secType <- "STK"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.symbol <- "INFY"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.exchange <- "NSE"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.currency <- "INR"</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.reqMktDepthEx(1, contract, 5)</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> let timerWorkflow = async { </span><br />
<span style="font-family: Verdana, sans-serif;"> do! Async.Sleep 10000</span><br />
<span style="font-family: Verdana, sans-serif;"> Application.Exit() </span><br />
<span style="font-family: Verdana, sans-serif;"> }</span><br />
<span style="font-family: Verdana, sans-serif;"> Async.Start timerWorkflow</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> Application.Run()</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.cancelMktDepth(1)</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.disconnect()</span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "Diconnected!"</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> printfn "end time: %A" (now())</span><br />
<span style="font-family: Verdana, sans-serif;"> 0</span><br />
<div>
<br /></div>
<br />
<br /></div>
</div>
<div style="text-align: center;">
<br /></div>
<div>
</div>
<div style="text-align: center;">
<br /></div>
<div style="-webkit-text-stroke-width: 0px; color: black; font-family: 'Times New Roman'; font-size: medium; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; margin: 0px; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px;">
<div style="text-align: center;">
<span style="font-family: Verdana, sans-serif;"> </span></div>
</div>
Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-21950779438107770242014-01-21T04:51:00.000-08:002014-01-21T04:52:53.842-08:00F#: Capture tick price data from InteractiveBrokersIn 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.<br />
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:<br />
<br />
<span style="font-family: Verdana, sans-serif;"> let contract = tws1.createContract()</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.secType <- "STK"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.symbol <- "SPY"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.exchange <- "SMART"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.currency <- "USD"</span><br />
<br />
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.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> tws1.reqMktDataEx(1, contract, "", 0)</span><br />
<br />
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.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> tws1.cancelMktData(1)</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
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:<br />
<br />
<span style="font-family: Verdana, sans-serif;"> tws1.tickPrice.AddHandler(</span><br />
<span style="font-family: Verdana, sans-serif;"> new TickPriceEventHandler(processTickPriceEvent))</span><br />
<br />
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.<br />
<br />
<span style="font-family: Verdana, sans-serif;"> let timerWorkflow = async { </span><br />
<span style="font-family: Verdana, sans-serif;"> do! Async.Sleep 10000</span><br />
<span style="font-family: Verdana, sans-serif;"> Application.Exit() </span><br />
<span style="font-family: Verdana, sans-serif;"> }</span><br />
<span style="font-family: Verdana, sans-serif;"> Async.Start timerWorkflow </span><br />
<br />
OK, now that we have covered all the major steps, below is the complete code of the program.<br />
<br />
/////////////////////////////////////////////////////////////////////////////////////////////////////////<br />
<span style="font-family: Verdana, sans-serif;">open AxTWSLib</span><br />
<span style="font-family: Verdana, sans-serif;">open TWSLib</span><br />
<span style="font-family: Verdana, sans-serif;">open System</span><br />
<span style="font-family: Verdana, sans-serif;">open System.Drawing</span><br />
<span style="font-family: Verdana, sans-serif;">open System.Windows.Forms</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">type TickPriceEvent = </span><br />
<span style="font-family: Verdana, sans-serif;"> AxTWSLib._DTwsEvents_tickPriceEvent</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">type TickPriceEventHandler =</span><br />
<span style="font-family: Verdana, sans-serif;"> AxTWSLib._DTwsEvents_tickPriceEventHandler</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">let now() = DateTime.Now</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">let processTickPriceEvent _ (e:TickPriceEvent) =</span><br />
<span style="font-family: Verdana, sans-serif;"> match e.tickType with</span><br />
<span style="font-family: Verdana, sans-serif;"> | 1 -> printfn "%A|bid=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | 2 -> printfn "%A|ask=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | 4 -> printfn "%A|last=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | 6 -> printfn "%A|high=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | 7 -> printfn "%A|low=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | 9 -> printfn "%A|close=%A" (now()) e.price</span><br />
<span style="font-family: Verdana, sans-serif;"> | x -> printfn "%A|UNKNOW_TICK_TYPE(%d)=%A" (now()) x e.price</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;">[<EntryPoint; STAThread>]</span><br />
<span style="font-family: Verdana, sans-serif;">let main _ = </span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "start time: %A" (now())</span><br />
<span style="font-family: Verdana, sans-serif;"> let form1 = new Form(Text="Dummy Form")</span><br />
<span style="font-family: Verdana, sans-serif;"> let tws1 = new AxTws()</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.BeginInit()</span><br />
<span style="font-family: Verdana, sans-serif;"> form1.Controls.Add(tws1)</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.EndInit()</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> tws1.connect("", 7496, 1)</span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "server version = %d" tws1.serverVersion</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> if tws1.serverVersion > 0 then</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.tickPrice.AddHandler(</span><br />
<span style="font-family: Verdana, sans-serif;"> new TickPriceEventHandler(processTickPriceEvent))</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> let contract = tws1.createContract()</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.secType <- "STK"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.symbol <- "SPY"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.exchange <- "SMART"</span><br />
<span style="font-family: Verdana, sans-serif;"> contract.currency <- "USD"</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.reqMktDataEx(1, contract, "", 0)</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> let timerWorkflow = async { </span><br />
<span style="font-family: Verdana, sans-serif;"> do! Async.Sleep 10000</span><br />
<span style="font-family: Verdana, sans-serif;"> Application.Exit() </span><br />
<span style="font-family: Verdana, sans-serif;"> }</span><br />
<span style="font-family: Verdana, sans-serif;"> Async.Start timerWorkflow</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> Application.Run()</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.cancelMktData(1)</span><br />
<span style="font-family: Verdana, sans-serif;"> tws1.disconnect()</span><br />
<span style="font-family: Verdana, sans-serif;"> printfn "Diconnected!"</span><br />
<span style="font-family: Verdana, sans-serif;"><br /></span>
<span style="font-family: Verdana, sans-serif;"> printfn "end time: %A" (now())</span><br />
<span style="font-family: Verdana, sans-serif;"></span><br />
<span style="font-family: Verdana, sans-serif;"> 0</span><br />
<div>
<br /></div>
<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com2tag:blogger.com,1999:blog-5063207632796600185.post-83085357579363262172014-01-05T22:15:00.000-08:002014-01-10T02:51:05.684-08:00F#: Connect to InteractiveBrokers via ActiveXIn this blog post we show how to write a F# console program which runs on Windows and connect to <a href="https://www.interactivebrokers.com/en/main.php" target="_blank">InteractiveBrokers</a> (IB). There are a few different <a href="https://www.interactivebrokers.com/en/?f=%252Fen%252Fsoftware%252Fibapi.php" target="_blank">options of programmatic connectivity</a> 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:<br />
<br />
<div style="text-align: center;">
<b>F# code <=> Tws.ocx <=> TWS <=> IB server.</b></div>
<br />
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.<br />
<br />
<b>Install TWS and API:</b><br />
<br />
<ul>
<li>Trade Workstation: please download and install a standalone version from <a href="https://www.interactivebrokers.com/en/index.php?f=tws&p=software" target="_blank">here</a>. You need to configure TWS to accept ActiveX clients by doing the following:</li>
<ul>
<li>Select Edit->Global Configuration</li>
<li>Select API->Settings</li>
<li>Check "Enable ActiveX and Socket Clients" and note down the port number specified at "Socket Port". By default, the port number is 7496.</li>
</ul>
<li>API: can be downloaded from <a href="http://interactivebrokers.github.io/#" target="_blank">here</a>, which contains Tws.ocx.</li>
</ul>
<br />
<b>Generate wrapper DLLs to expose Tws.ocx to .NET:</b><br />
Assuming that the API has been installed at C:\TWS_API and you have Visual Studio loaded on your machine, you can then run <a href="http://msdn.microsoft.com/en-us/library/8ccdh774(v=vs.110).aspx" target="_blank">Aximp.exe</a> as follows:<br />
<br />
<span style="font-family: Courier New, Courier, monospace;">C:\TWS_API\bin\ActiveX>aximp Tws.ocx</span><br />
<br />
And you shall see the following 2 DLLs generated:<br />
<br />
<span style="font-family: Courier New, Courier, monospace;">Generated Assembly: C:\TWS_API\bin\ActiveX\TWSLib.dll</span><br />
<span style="font-family: Courier New, Courier, monospace;">Generated Assembly: C:\TWS_API\bin\ActiveX\AxTWSLib.dll</span><br />
<br />
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.)<br />
<br />
<b>Write F# code to connect and disconnect:</b><br />
Below you can see the code of our console program. There are some remarks about the code I would like to make:<br />
<br />
<ul>
<li>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.</li>
<li>To connect to TWS, we invoke the connect method, which takes 3 parameters: IP address, port number, client ID. </li>
<ul>
<li>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.</li>
<li>port number: 7496 by default, unless you change it in the Global Configuration page of TWS.</li>
<li>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. </li>
</ul>
<li>Once TWS receives a connection request, a dialog pops up. Click on yes to accept the connection. </li>
<li>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. </li>
</ul>
<br />
<br />
------------------------------------------------------------------------<br />
<span style="font-family: Courier New, Courier, monospace;">open AxTWSLib</span><br />
<span style="font-family: Courier New, Courier, monospace;">open System</span><br />
<span style="font-family: Courier New, Courier, monospace;">open System.Drawing</span><br />
<span style="font-family: Courier New, Courier, monospace;">open System.Windows.Forms</span><br />
<span style="font-family: Courier New, Courier, monospace;"><br /></span>
<span style="font-family: Courier New, Courier, monospace;">[<EntryPoint; STAThread>]</span><br />
<span style="font-family: Courier New, Courier, monospace;">let main _ = </span><br />
<span style="font-family: Courier New, Courier, monospace;"> // initialize TWS ActiveX Control tws1</span><br />
<span style="font-family: Courier New, Courier, monospace;"> let form1 = new Form(Text="Dummy Form")</span><br />
<span style="font-family: Courier New, Courier, monospace;"> let tws1 = new AxTws()</span><br />
<span style="font-family: Courier New, Courier, monospace;"> tws1.BeginInit()</span><br />
<span style="font-family: Courier New, Courier, monospace;"> form1.Controls.Add(tws1)</span><br />
<span style="font-family: Courier New, Courier, monospace;"> tws1.EndInit()</span><br />
<span style="font-family: Courier New, Courier, monospace;"><br /></span>
<span style="font-family: Courier New, Courier, monospace;"> // connect to local TWS</span><br />
<span style="font-family: Courier New, Courier, monospace;"> tws1.connect("", 7496, 1)</span><br />
<span style="font-family: Courier New, Courier, monospace;"> printfn "server version = %d" tws1.serverVersion</span><br />
<span style="font-family: Courier New, Courier, monospace;"><br /></span>
<span style="font-family: Courier New, Courier, monospace;"> // disconnect from TWS</span><br />
<span style="font-family: Courier New, Courier, monospace;"> if tws1.serverVersion > 0 then </span><br />
<span style="font-family: Courier New, Courier, monospace;"> tws1.disconnect()</span><br />
<span style="font-family: Courier New, Courier, monospace;"><br /></span>
<span style="font-family: Courier New, Courier, monospace;"> 0</span><br />
------------------------------------------------------------------------<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com4tag:blogger.com,1999:blog-5063207632796600185.post-16735454280849094482013-01-25T10:27:00.001-08:002013-05-29T00:32:01.899-07:00F#: contributing a chapter to an upcoming book, "F# Deep Dives"Thanks to <a href="http://tomasp.net/blog/manning-deep-dives.aspx">Tomas Petricek's invitation</a>, 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 <b>4:</b> Numerical computing in financial domain. The first draft of the chapter is now available online via <a href="http://manning.com/petricek2/">Manning's early access program</a>.<br />
<div>
<br /></div>
<div>
Currently the chapter covers the following:</div>
<div>
<div class="Head1">
<span lang="EN-US"><span style="font-family: inherit;">1 Introducing financial derivatives<o:p></o:p></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">2 Using probability functions of Math.NET</span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">2.1 Configuring F# Interactive (profiling and floating point formatting)<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">2.2 Setting up Math.NET Numerics<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">2.3 Introducing Random Variables, Expectation
and Variance<o:p></o:p></span></span></span></div>
<div class="Head1">
<span style="font-family: inherit;"><span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;">2.4 </span></span><span lang="EN-US">Generating Normal</span><span lang="EN-US">ly Distributed Samples</span></span></div>
<div class="SidebarHead">
<span lang="EN-US"><span style="font-family: inherit;"><o:p></o:p></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3 Geometric Brownian Motion and Monte Carlo
Estimate<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3.1 Modeling stock prices using geometric
Brownian motion<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3.2 Payoff Function, Discounted Payoff, and
Monte Carlo estimate<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3.3 Analyzing variance of Monte Carlo estimates<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3.4 Pricing Path-dependent options (Asian options and barrier options)<o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;">3.5 Reducing Variance by antithetic
variates </span><span style="font-family: Calibri, sans-serif;"><o:p></o:p></span></span></span></div>
<div class="Head1">
<span lang="EN-US"><span lang="EN-US" style="font-size: 11pt; line-height: 115%;"><span style="font-family: inherit;"><br /></span></span></span></div>
</div>
Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com2tag:blogger.com,1999:blog-5063207632796600185.post-77446783854566226402012-09-07T05:41:00.001-07:002012-09-07T05:44:28.341-07:00F#: K-S test on final prices of GBM pathsIn the <a href="http://programmingcradle.blogspot.sg/2012/09/f-simulate-entire-gbm-path.html">previous post</a>, 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 <a href="http://connect.microsoft.com/VisualStudio/feedback/details/634761/system-random-serious-bug">raised this issue</a> 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. <br /><br />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., <strike>close my eyes to this issue</strike> do <a href="http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test">Kolmogrov-Smirnov test</a> (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.<br /><br />Basically the F# code is pretty much the same as in the previous post except the following: <br /><br /><ul><br /><li>I add a run_stat_test function to run the statistical test and generate a Q-Q plot in Matlab.</li><br /><li>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.</li><br /></ul><br /><strong>// Monte Carlo parameters<br />let T = 1.0 // time horizon in years<br />let M = 100000 // number of simulations<br />let N = 1 // number of time intervals</strong><br /><br />You can see that I take the logarithms of all final prices and test them against a normal distribution with mean = <strong>(log S0)+(r-0.5*sigma*sigma)*T</strong> and standard deviation = <strong>sigma*sqrt(T)</strong> using KS test. And the output from the test is as follows:<br /><br /><strong>h = false, p = 0.7243765247</strong><br /><br />And here is the Q-Q plot:<br /><br /><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9VQPXPEQ-7CCDjZZCDwtCdekIenXnj15rHpRcTnMV-OXoek5OHTDQnti7R1CSc09D13zcdeThfFNI1ErqRFj2GrhlrIu46MKupCEhVIUcpbxhKZf5cZ269qcti-mJ10ttO40bvKJI-0b7/s1600-h/untitled%25255B5%25255D.png"><img alt="untitled" border="0" height="476" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhkiqTi61p62eb7Jvgrog4IwJyooPbol2EOXDia9VwjbTsTdw28AxGAsF5YULcosLVm5EX8mFe-H-veWMaixX-2uINnzYpqTXjUjHnjVyMY6hWBTZB08t7DNtgpxnqw2dl1idE08Wz-KwmX/?imgmax=800" style="background-image: none; border: 0px currentColor; display: inline; padding-left: 0px; padding-right: 0px; padding-top: 0px;" title="untitled" width="570" /></a><br /><br />p-value is about 0.72, which means it’s quite likely that the simulated prices were indeed sampled from the normal distribution, and <strong>h=false</strong> 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.<br /><br />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.<br /><br />//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br /><br />open MathNet.Numerics.Distributions<br />open MathNet.Numerics.Statistics<br /><br />// Monte Carlo parameters<br />let T = 1.0 // time horizon in years<br />let M = 100000 // number of simulations<br />let N = 1 // number of time intervals<br />let p = 0.99 // confidence level<br /><br />// underlying dynamics<br />let S0 = 50.0 // spot price<br />let sigma = 0.2 // annualized vol<br />let r = 0.01 // risk-free rate<br /><br />// risk factor dW<br />let rnd = new System.Random(1)<br />let get_dW dt =<br /> let dW = Normal.WithMeanVariance(0.0, dt)<br /> dW.RandomSource <- rnd<br /> (fun () -> dW.Sample())<br /><br />// GBM path generation<br />let generate_GBM_paths S0 r sigma T M N =<br /> let dt = T / (float N)<br /> let drift = (r - 0.5 * (sigma**2.0)) * dt<br /> let dW = get_dW dt<br /> Array.init M (fun _ -> <br /> Array.init N (fun _ -> drift + sigma * dW()) |><br /> Array.scan (+) 0.0 |><br /> Array.map (fun x -> S0 * exp(x)) )<br /><br />let run_stat_test (prices:float array) =<br /> let m = (log S0)+(r-0.5*sigma*sigma)*T<br /> let s = sigma*sqrt(T)<br /> let log_prices = Array.map log prices<br /> let matlab = new MLApp.MLAppClass()<br /> matlab.PutWorkspaceData("log_prices", "base", log_prices)<br /> matlab.PutWorkspaceData("m", "base", m)<br /> matlab.PutWorkspaceData("s", "base", s)<br /> matlab.Execute("pd = ProbDistUnivParam('normal',[m s])") |> ignore<br /> matlab.Execute("[h,p] = kstest(log_prices,pd,0.05,0)") |> ignore<br /> matlab.Execute("qqplot(log_prices,pd)") |> ignore<br /> let h = matlab.GetVariable("h","base")<br /> let p = matlab.GetVariable("p","base")<br /> (h,p)<br /><br />// main<br />let main () =<br /> let final_prices = generate_GBM_paths S0 r sigma T M N |> <br /> Array.map (fun p -> p.[p.Length - 1])<br /> let forward = final_prices.Mean()<br /> let std_dev = final_prices.StandardDeviation()<br /> let std_error = std_dev / sqrt(float M)<br /> let std_normal = Normal.WithMeanVariance(0.0,1.0)<br /> let radius = std_normal.InverseCumulativeDistribution(0.5*(1.0+p)) <br /> * std_error<br /> let conf_int = (forward - radius, forward + radius)<br /> let h,p = run_stat_test final_prices<br /> printfn "h = %A, p = %A" h p<br /> (forward, std_dev, std_error, conf_int)<br /><br />// run simulation<br />main()<br /><br />// calculate theoretical mean and deviation<br />let true_forward = S0 * exp(r*T)<br />let v = sigma*sigma*T<br />let true_std_dev = sqrt((exp(v)-1.0)*exp(2.0*(log(S0)+r*T-0.5*v*T)+v))Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-85421947819609375562012-09-02T23:30:00.000-07:002012-09-03T02:01:08.569-07:00F#: Simulate entire GBM pathIn 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.<br />
<br />
In this post, I started using a open-source numerical library, <a href="http://numerics.mathdotnet.com/">Math.NET Numerics</a>. Don Syme wrote <a href="http://blogs.msdn.com/b/dsyme/archive/2012/07/06/getting-started-with-math-net-and-f-programming.aspx">a blog</a> on how to install and use it. Basically I use it in my code for two purposes:<br />
<br />
<ul>
<li>using the Normal class to draw normal random numbers. Internally the Normal class also uses the <a href="http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform">Box Muller transform</a> 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 <a href="http://programmingcradle.blogspot.sg/2012/08/f-how-could-monte-calro-return.html">infinity issue</a> I mentioned in a previous post. </li>
<li>using some handy statistical functions, like mean, standard deviation, and inverse CDF, because I like to compute standard error and confidence interval. </li>
</ul>
<div>
Before we look into the code, here are the parameters I use:<br />
<br />
<br />
<b>// Monte Carlo parameters</b><br />
<b>let T = 1.0 // time horizon in years</b><br />
<b>let M = 100000 // number of simulations</b><br />
<b>let N = 252 // number of time intervals</b><br />
<b>let p = 0.99 // confidence level</b><br />
<b><br /></b>
<b>// underlying dynamics</b><br />
<b>let S0 = 50.0 // spot price</b><br />
<b>let sigma = 0.2 // annualized vol</b><br />
<b>let r = 0.01 // risk-free rate</b><br />
<br />
<br />
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.)</div>
<div>
<br /></div>
<div>
<div>
<span style="font-family: inherit;"><b>// risk factor dW</b></span></div>
<div>
<span style="font-family: inherit;"><b>let rnd = new System.Random(1)</b></span></div>
<div>
<span style="font-family: inherit;"><b>let get_dW dt =</b></span></div>
<div>
<span style="font-family: inherit;"><b> let dW = Normal.WithMeanVariance(0.0, dt)</b></span></div>
<div>
<span style="font-family: inherit;"><b> dW.RandomSource <- rnd</b></span></div>
<div>
<span style="font-family: inherit;"><b> (fun () -> dW.Sample())</b></span></div>
</div>
<div>
<br /></div>
<div>
The below function, <b>generate_GBM_paths</b>, 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.<br />
<br />
<br />
<b>// GBM path generation based on the Euler method</b><br />
<b>let generate_GBM_paths S0 r sigma T M N =</b><br />
<b> let dt = T / (float N)</b><br />
<b> let drift = (r - 0.5 * (sigma**2.0)) * dt</b><br />
<b> let dW = get_dW dt</b><br />
<b> Array.init M (fun _ -> </b><br />
<b> Array.init N (fun _ -> drift + sigma * dW()) |></b><br />
<b> Array.scan (+) 0.0 |></b><br />
<b> Array.map (fun x -> S0 * exp(x)) )</b><br />
<br />
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:<br />
<br />
<ul>
<li>the statistical functions provided by Math.NET Numerics to compute </li>
<ul>
<li>mean, i.e., the forward price of the underlying</li>
<li>standard deviation of the simulated final price distribution</li>
<li>standard error of the final price simulation</li>
<li>confidence interval</li>
</ul>
<li>Matlab to plot KDE (<a href="http://en.wikipedia.org/wiki/Kernel_density_estimation">Kernel Density Estimation</a>) 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.</li>
</ul>
<br />
Finally, here are the results from the main function:<br />
<br />
<br />
<b>Real: 00:00:07.541, CPU: 00:00:06.817, GC gen0: 201, gen1: 135, gen2: 2</b><br />
<b>val it : float * float * float * (float * float) =</b><br />
<b> (50.47657684, 10.21848436, 0.03231368481, (50.3933423, 50.55981138))</b><br />
<br />
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:<br />
<br />
<b>val true_forward : float = 50.50250835</b><br />
<b>val true_std_dev : float = 10.20235347</b><br />
<br />
And below is the KDE plot:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh0IeG96GAbnccAZmzduXsJqpvWBfhWqaqGL9fxMIhb8WxfqMNjlC776Xdec2Yh7kxoGVumE8XffJNRPGvYA3qCk2dgkB28r-yR5nU6tkdr2vtjlDs73KZWHczM0XY4YkBPy16nAromuCv4/s1600/untitled.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="257" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh0IeG96GAbnccAZmzduXsJqpvWBfhWqaqGL9fxMIhb8WxfqMNjlC776Xdec2Yh7kxoGVumE8XffJNRPGvYA3qCk2dgkB28r-yR5nU6tkdr2vtjlDs73KZWHczM0XY4YkBPy16nAromuCv4/s640/untitled.png" width="640" /></a></div>
<br />
///////////////////////////////////////////////////////////////////////////////////////////////////////<br />
<br /></div>
<br />
<b>open MathNet.Numerics.Distributions</b><br />
<b>open MathNet.Numerics.Statistics</b><br />
<b><br /></b>
<b>// Monte Carlo parameters</b><br />
<b>let T = 1.0 // time horizon in years</b><br />
<b>let M = 100000 // number of simulations</b><br />
<b>let N = 252 // number of time intervals</b><br />
<b>let p = 0.99 // confidence level</b><br />
<b><br /></b>
<b>// underlying dynamics</b><br />
<b>let S0 = 50.0 // spot price</b><br />
<b>let sigma = 0.2 // annualized vol</b><br />
<b>let r = 0.01 // risk-free rate</b><br />
<b><br /></b>
<b>// risk factor dW</b><br />
<b>let rnd = new System.Random(1)</b><br />
<b>let get_dW dt =</b><br />
<b> let dW = Normal.WithMeanVariance(0.0, dt)</b><br />
<b> dW.RandomSource <- rnd</b><br />
<b> (fun () -> dW.Sample())</b><br />
<b><br /></b>
<b>// GBM path generation</b><br />
<b>let generate_GBM_paths S0 r sigma T M N =</b><br />
<b> let dt = T / (float N)</b><br />
<b> let drift = (r - 0.5 * (sigma**2.0)) * dt</b><br />
<b> let dW = get_dW dt</b><br />
<b> Array.init M (fun _ -> </b><br />
<b> Array.init N (fun _ -> drift + sigma * dW()) |></b><br />
<b> Array.scan (+) 0.0 |></b><br />
<b> Array.map (fun x -> S0 * exp(x)) )</b><br />
<b><br /></b>
<b>// plot the simulated distribution using KDE in Matlab</b><br />
<b>let plot_KDE (prices:float array) =</b><br />
<b> let matlab = new MLApp.MLAppClass()</b><br />
<b> matlab.PutWorkspaceData("x", "base", prices)</b><br />
<b> matlab.Execute("[f,xi] = ksdensity(x);") |> ignore</b><br />
<b> matlab.Execute("plot(xi,f);") |> ignore</b><br />
<b> ()</b><br />
<b><br /></b>
<b>// main</b><br />
<b>let main () =</b><br />
<b> let final_prices = generate_GBM_paths S0 r sigma T M N |> </b><br />
<b> Array.map (fun p -> p.[p.Length - 1])</b><br />
<b> let forward = final_prices.Mean()</b><br />
<b> let std_dev = final_prices.StandardDeviation()</b><br />
<b> let std_error = std_dev / sqrt(float M)</b><br />
<b> let std_normal = Normal.WithMeanVariance(0.0,1.0)</b><br />
<b> let radius = std_normal.InverseCumulativeDistribution(0.5*(1.0+p)) </b><br />
<b> * std_error</b><br />
<b> let conf_int = (forward - radius, forward + radius)</b><br />
<b> plot_KDE final_prices</b><br />
<b> (forward, std_dev, std_error, conf_int)</b><br />
<b><br /></b>
<b>// run simulation</b><br />
<b>main()</b><br />
<b><br /></b>
<b>// calculate theoretical mean and deviation</b><br />
<b>let true_forward = S0 * exp(r*T)</b><br />
<b>let v = sigma*sigma*T</b><br />
<b>let true_std_dev = sqrt((exp(v)-1.0)*exp(2.0*(log(S0)+r*T-0.5*v*T)+v))</b><br />
<div>
<br /></div>
Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-1272451605838406422012-08-24T00:07:00.000-07:002012-08-24T05:04:34.571-07:00F#: Running Monte Carlo Simulation on GPUMulti-threading on CPU is a bit old fashioned. Let's run Monte Carlo Simulation <strike>with style</strike> on GPU. Here I'm using <a href="http://research.microsoft.com/en-us/projects/accelerator/">the Accelerator API</a> from Microsoft Research. Tomas Petricek has a nice and concise <a href="http://tomasp.net/blog/accelerator-intro.aspx">introduction</a> about the API. In fact, this post mainly relies on the introduction. If you never use the API, it's a good idea to take a look at the article to pick up some basic concepts on how the API works.<br />
<br />
Here are some points I like to mention about the sample code below:<br />
<ul>
<li>Currently Accelerator does not support random number generation. Since I want to time only the execution done on GPU, I put the code for normal random number generation outside the main pricing function (i.e., <strong>European_Option_Price</strong>), as what I did in my previous posts on Monte Carlo. Inside the the pricing function, I try to put every computation on GPU.</li>
<li>As I still use the Box Muller transform, to avoid the issue of generating a normal random number of infinity, I apply the simple transform (i.e., Y = 1-U[0,1)) mentioned in <a href="http://programmingcradle.blogspot.sg/2012/08/f-how-could-monte-calro-return.html">the blog post</a>. Well, if you want to try, you can choose not to apply the transform and set the seed to 388, and then see what your GPU will return. I tried that for fun.</li>
<li>Why I set grid_size = 3163 is because the number of simulations will be the square of grid_size, which is slightly over the simulation size (10 millions) I used in my previous posts on Monte Carlo. Easier for me to compare.</li>
<li>Obviously, using the API reduces F#'s expressiveness, but with some custom syntactic sugars, like overloaded operators, I believe that it can be mitigated. However, I don't bother to do that at this moment.</li>
</ul>
After running the code, the output is as follows:<br />
<br />
<strong>option_price = 10.117122</strong><br />
<strong><div>
CPU time = 234ms</div>
<div>
Absolute time = 433ms</div>
</strong><br />
Wow, in terms of absolute time, this GPU version is almost 4 times faster than my previous champion,<a href="http://programmingcradle.blogspot.sg/2012/08/f-use-find-and-replace-to-speed-up.html"> the Array version</a>. Though we fix the seed to 1, the reason why the option price given by the GPU version is different from the price by the Array version is mainly because of the simple transform I mentioned above. If you don't apply the transform, you will see both versions produce "almost" the same option price (some very minor difference is still there, due to slightly different number of simulations).<br />
<br />
My laptop is 4 years old and equiped with only a <strike>lousy</strike> Intel 965 Express chipset, because I seldom play PC games and hence I used to look at mainly CPU speed and memory when purchasing a laptop. Now, given this impressive performance result, next time when I buy a new laptop, I will definitely choose one with a high-end NVIDIA graphics chipset. <br />
<br />
Why NVIDIA? because I like to try the following:<br />
<ul>
<li>Actually the latest release of Accelerator contains a CUDA target.</li>
<li>NVIDIA claims their <a href="http://developer.nvidia.com/cuda/curand">cuRAND</a> library can generate random numbers at light speed. </li>
</ul>
<span style="font-family: "Courier New", Courier, monospace;"></span><br />
<span style="font-family: "Courier New", Courier, monospace;">////////////////////////////////////////////////////////////</span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
open System</div>
<div>
open System.Diagnostics</div>
<div>
open Microsoft.ParallelArrays</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">type PA = Microsoft.ParallelArrays.ParallelArrays</span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
type FPA = Microsoft.ParallelArrays.FloatParallelArray</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">let grid_size = 3163</span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
let shape = [| grid_size; grid_size |]</div>
<div>
let dxTarget = new DX9Target()</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">let rnd = new Random(1) </span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
let next_double_pair () = </div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;">(1.0 - rnd.NextDouble(), 1.0 - rnd.NextDouble())</span></div>
<span style="font-family: "Courier New", Courier, monospace;"><div>
//let rnd = new Random(388) </div>
<div>
//let next_double_pair () = (rnd.NextDouble(), rnd.NextDouble())</div>
<div>
let box_muller_transform (u1,u2) = </div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;"> sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)</span></div>
<span style="font-family: "Courier New", Courier, monospace;"><div>
let generator _ _ = </div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;"> next_double_pair() |> box_muller_transform |> (float32)</span></div>
<span style="font-family: "Courier New", Courier, monospace;"><div>
let random_normals = Array2D.init grid_size grid_size generator</div>
</span><span style="font-family: "Courier New", Courier, monospace;">let European_Option_Price K S0 r vol T =</span><span style="font-family: "Courier New", Courier, monospace;"><div>
let FPA_zero = new FPA(float32 0.0, shape)</div>
<div>
let FPA_N = new FPA(float32 (grid_size * grid_size), [|1|])</div>
<div>
let FPA_discount = new FPA(float32 (exp(-r*T)), [|1|])</div>
<div>
let FPA_K = new FPA(float32 K,shape)</div>
<div>
let FPA_S0 = new FPA(float32 S0,shape)</div>
<div>
let FPA_mean = </div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;"><span style="font-family: Times New Roman;"> </span>new FPA(float32 ((r - 0.5*(vol**2.0))*T), shape)</span></div>
<span style="font-family: "Courier New", Courier, monospace;"><div>
let FPA_sigma = new FPA(float32 (vol*sqrt(T)), shape)</div>
<div>
let FPA_exp = new FPA(float32 (exp 1.0), shape)</div>
<div>
let FPA_normals = new FPA(random_normals)</div>
<div>
let FPA_S_T = </div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;"> FPA_S0 * </span></div>
<span style="font-family: "Courier New", Courier, monospace;"> PA.Pow(FPA_exp, </span><br />
<span style="font-family: "Courier New", Courier, monospace;"> (FPA_mean + FPA_sigma * FPA_normals))</span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
let FPA_payouts = PA.Max(FPA_S_T - FPA_K, FPA_zero)</div>
<div>
let FPA_price = FPA_discount * PA.Sum(FPA_payouts) / FPA_N </div>
<div>
let result = dxTarget.ToArray1D(FPA_price)</div>
<div>
(float) result.[0]</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">let S0 = 50.0</span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
let K = 40.0</div>
<div>
let r = 0.01</div>
<div>
let vol = 0.2</div>
<div>
let T = 0.25</div>
<div>
let price() = </div>
<div>
let price = European_Option_Price K S0 r vol T</div>
<div>
printfn "option_price = %f" price</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">let time f = </span><br />
<span style="font-family: "Courier New", Courier, monospace;"><div>
let proc = Process.GetCurrentProcess()</div>
<div>
let cpu_time_stamp = proc.TotalProcessorTime</div>
<div>
let timer = new Stopwatch()</div>
<div>
timer.Start()</div>
<div>
try</div>
<div>
f()</div>
<div>
finally</div>
<div>
let cpu_time = (proc.TotalProcessorTime -</div>
</span><div>
<span style="font-family: "Courier New", Courier, monospace;"> cpu_time_stamp).TotalMilliseconds</span></div>
<span style="font-family: "Courier New", Courier, monospace;"><div>
printfn "CPU time = %dms" (int64 cpu_time)</div>
<div>
printfn "Absolute time = %dms" timer.ElapsedMilliseconds</div>
</span><br />
<span style="font-family: "Courier New", Courier, monospace;">time price</span><br />
<div>
</div>
<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-81128221601290948262012-08-20T20:20:00.001-07:002012-08-25T05:33:10.882-07:00F#: Infinity and Infinitesimal-nessThese days many programming languages have the notion of infinity defined explicitly with their float type. F# (or .Net actually) is no exception. In contrast to infinity, <strong><a href="http://msdn.microsoft.com/en-us/library/system.double.epsilon">Double.Epsilon</a></strong> is the smallest positive float in F#, i.e., anything between <strong>Double.Epsilon</strong> and 0.0 is viewed as 0.0 in F#. Let's use F# Interactive to do some simple empirical studies on both infinity and <strong>Double.Epsilon</strong> in F#.<br />
<br />
First of all, F# has a built-in float value of <strong><a href="http://msdn.microsoft.com/en-us/library/ee353754">infinity</a></strong>, which is equivalent to <strong><a href="http://msdn.microsoft.com/en-us/library/system.double.positiveinfinity.aspx">System.Double.PositiveInfinity</a></strong> in .Net.<br />
<br />
<strong>> infinity;;<br />val it : float = infinity<br />> System.Double.PositiveInfinity;;<br />val it : float = infinity</strong><br />
<br />
Well, negative infinity is also there.<br />
<br />
<strong>> -infinity;;<br />val it : float = -infinity<br />> System.Double.NegativeInfinity;;<br />val it : float = -infinity</strong><br />
<br />
Then let's see how small <strong>Double.Epsilon</strong> is:<br />
<br />
<strong> > System.Double.Epsilon;;<br />val it : float = 4.940656458e-324</strong><br />
<br />
Not sure how you feel about this, but it looks "scary" small to me. By the way, we also have <strong>Double.MaxValue</strong> and <strong>Double.MinValue</strong>.<br />
<br />
<strong>> System.Double.MaxValue;;<br />val it : float = 1.797693135e+308<br />> infinity > System.Double.MaxValue;;<br />val it : bool = true<br />> System.Double.MinValue;;<br />val it : float = -1.797693135e+308<br />> System.Double.MinValue > -infinity;;<br />val it : bool = true</strong><br />
<br />
In financial mathematics, log(x) and exp(x) are used excessively, so let's take a look at how these 2 functions operate with infinity and <strong>Double.Epsilon</strong>. <br />
<br />
<strong>> exp(infinity);;<br />val it : float = infinity<br />> log(infinity);;<br />val it : float = infinity</strong><br />
<strong>> log(0.0);;<br />val it : float = -infinity</strong><br />
<br />
We know that log(x) goes to negative infinity when x goes to 0.0. What if we evaluate log(x) at x = <strong>Double.Epsilon</strong>? Can F# properly handle this seemingly extreme case?<br />
<br />
<strong>> log(System.Double.Epsilon);;<br />val it : float = -744.4400719<br />> exp(log(System.Double.Epsilon));;<br />val it : float = 4.940656458e-324<br />> exp(log(System.Double.Epsilon)) = System.Double.Epsilon;;<br />val it : bool = true</strong><br />
<br />
OK, <strong>log(Double.Epsilon) </strong>is only around -744.44, not too bad relative to the infinitesimal-ness of Double.Epsilon, and if you substitute this value into exp(x) you basically recover <strong>Double.Epsilon </strong>as shown above. If you run similar tests in Matlab, Matlab will return you pretty much the same results. <br />
<br />
The above example actually implies a few things about the resolution of log(x) and exp(x) in F#:<br />
<ul>
<li>Other than negative infinity, you won't see that log(x) could ever return any numerical value lower than <strong>-744.4400719</strong>, i.e., <strong>log(Double.Epsilon)</strong> </li>
<li>For any x smaller than <strong>-744.4400719,</strong> i.e., <strong>log(Double.Epsilon)</strong>, exp(x) is evaluated as 0.0.</li>
</ul>
Besides exp(x) and log(x), you might also try how well F# handles sin(x) and cos(x) when x is extremely big or small. You might see something unexpected.<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-56054533569426209582012-08-16T22:25:00.000-07:002012-08-18T06:39:19.291-07:00F#: How could the Monte Calro pricer return option_price = infinity?If you leave the random seed blank in <a href="http://programmingcradle.blogspot.sg/2012/08/f-use-find-and-replace-to-speed-up.html">the Monte Carlo code</a> and let <strong>System.Random</strong> pick up a time-dependent seed for you (or you deliberately set the seed to 388), you might see a funny pricing result like this:<br />
<br />
<strong>option_price = Infinity<br />CPU time = 2636ms<br />Absolute time = 2435ms</strong><br />
<br />
Why you see this funny "infinity" price is because:<br />
<ul>
<li><strong>System.Random.NextDouble()</strong> could generate 0.0, though the probability is extremely low. As such, it's possible to generate a tuple of <strong>(u1,u2)</strong> where <strong>u1 = 0.0</strong> and<strong> (u2 < 0.25 or u2 > 0.75)</strong>. This kind of <strong>(u1,u2)</strong> drives <a href="http://en.wikipedia.org/wiki/Box-Muller_transform">the Box Muller transform</a> to infinity. Then, not surprisingly, we can see that a single sample of infinite payout screws up the <a href="http://en.wikipedia.org/wiki/Law_of_Large_Numbers">Law of Large Numbers</a>, as the integrability condition is not satisfied.</li>
<li>The notion of infinity is defined in F# (or actually via <strong>Double.PositiveInfinity</strong> and <strong>Double.NegativeInfinity</strong> in .Net). F# knows how to handle <strong>log(0.0)</strong> and <strong>exp(infinity)</strong> so that F# elegantly tells you "Sir, the option is literally priceless." rather than returns you a garbage float number. Nevertheless, elegance without correctness means nothing, as infinity is a wrong answer to our option pricing problem .</li>
</ul>
Well, there's nothing wrong for <strong>System.Random</strong> to generate a value of 0.0, because this is simply <a href="http://msdn.microsoft.com/en-us/library/system.random.nextdouble.aspx">how it's designed</a>. The <strong>NextDouble()</strong> method is meant to approximate <strong>U[0,1)</strong>, but the possibility of getting a 0.0 value is not desirable here. One workaround is to take a simple transform: <strong>Y = 1.0 - U[0,1)</strong>. Then <strong>Y</strong> follows <strong>U(0,1]</strong>. So the code snippet we use to draw random numbers can be modified as follows:<br />
<br />
<strong>let N = 10000000<br />let rnd = new Random(388) <br />let next_double () = 1.0 - rnd.NextDouble()<br />let randomPairs = Array.init N (fun i -> (next_double(), next_double()))<br />let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)<br />let randomNormals = randomPairs |> Array.map box_muller_transform</strong><br />
<br />
You might argue that mathematically <strong>P(U=0)</strong> is zero in either <strong>U(0,1)</strong> or <strong>U[0,1)</strong> or <strong>U(0,1] </strong>or <b>U[0,1]</b>. When it comes to implementation of a pseudo-random number generator, at least in the case of <strong>System.Random</strong>, <strong>U[0,1)</strong> usually implies <strong>P(U=0)>0</strong>.<br />
<br />
Perhaps you like to know how I find out the seed of 388. For sure, I don't pick up this seed out of thin air. Here is the F# code I use to search for a seed that leads to infinity in our case:<br />
<br />
<strong>///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////</strong><br />
<strong>open System<br />open System.Diagnostics</strong><br />
<br />
<strong>let search_for_0 seed_min seed_max draw_max =<br /> let found u1 u2 = if (u1 = 0.0 && (u2 < 0.25 || u2 > 0.75)) then true else false</strong><br />
<strong> let rec seed_loop seed =<br /> let rnd = new Random(seed)<br /> <br /> let rec draw_loop draw =<br /> let (u1,u2) = (rnd.NextDouble(),rnd.NextDouble())<br /> if found u1 u2 || draw = draw_max then (u1,u2, draw)<br /> else draw_loop (draw + 1)<br /> <br /> let (u1,u2, draw) = draw_loop 1<br /> if found u1 u2 || seed = seed_max then (u1,u2, draw, seed)<br /> else seed_loop (seed + 1)<br /> <br /> let (u1,u2, draw, seed) = seed_loop seed_min<br /> printfn "seed_min = %d, seed_max = %d, draw_max = %d" seed_min seed_max draw_max<br /> if found u1 u2 then </strong><br />
<strong> printfn "found u1 = 0.0 and u2 = %f at draw %d of seed %d" u2 draw seed<br /> else printfn "couldn't find 0.0 within the range specified."</strong><br />
<br />
<strong>let time () = <br /> let proc = Process.GetCurrentProcess()<br /> let cpu_time_stamp = proc.TotalProcessorTime<br /> let timer = new Stopwatch()<br /> timer.Start()<br /> try<br /> search_for_0 1 1000 10000000<br /> finally<br /> let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds<br /> printfn "CPU time = %dms" (int64 cpu_time)<br /> printfn "Absolute time = %dms" timer.ElapsedMilliseconds</strong><br />
<strong>time()</strong><br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-68920632010738224772012-08-15T19:46:00.000-07:002012-08-16T01:29:22.132-07:00F#: Monte Carlo Simulation using Tail RecursionPreviously we used <a href="http://programmingcradle.blogspot.sg/2012/08/f-monte-carlo-pricer-for-european-calls.html">lists</a> or <a href="http://programmingcradle.blogspot.sg/2012/08/f-use-find-and-replace-to-speed-up.html">arrays</a> to implement Monte Carlo simulation, where in the beginning of the code we draw a fixed number of random numbers and keep all of them in memory throughout the execution of the code. Needless to say it's very memory-consuming. Indeed, if I increase the number of simulation runs to 100 millions, the code (either the List or Array version) throws <strong>System.OutOfMemoryException</strong> on my laptop.<br />
<br />
We in fact don't really need to keep all the random numbers in memory, if the option price is the only thing we want to calculate, which is the case in the simple option pricing problem we are handling. At the end of the day, what really matters and hence we have to keep track of is the sum of all simulated option payouts rather than individual simulated payouts. In other words, every time when we generate a random number to simulate an option payout, there's no harm to throw away the random number as long as we have added the simulated payout to a partial sum value we keep track of.<br />
<br />
That's where tail recursion can help. We want to use tail recursion with the accumulator pattern to achieve constant space complexity. In the code below, the sum parameter of the loop function keeps track of the partial sum and plays the role of accumulator.<br />
<br />
Well, here is the performance results for 10 million simulation runs on my laptop:<br />
<table border="0" cellpadding="0" cellspacing="0" class="MsoNormalTable" style="border-collapse: collapse; margin: auto auto auto 4.65pt; mso-padding-alt: 0cm 5.4pt 0cm 5.4pt; mso-yfti-tbllook: 1184; width: 533px;">
<tbody>
<tr style="height: 15pt; mso-yfti-firstrow: yes; mso-yfti-irow: 0;"><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0) windowtext windowtext; border-style: solid none solid solid; border-width: 1pt 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">version<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">option price<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">CPU time (ms)<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext windowtext windowtext rgb(0, 0, 0); border-style: solid solid solid none; border-width: 1pt 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">absolute time (ms)<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 1;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">List<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">5,335<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">5,467<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 2;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Array<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">1,606<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">1,639<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 3; mso-yfti-lastrow: yes;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext windowtext; border-style: none none solid solid; border-width: 0px 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Tail Recursion</span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">2,730</span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext windowtext rgb(0, 0, 0); border-style: none solid solid none; border-width: 0px 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">2,732</span></span></div>
</td></tr>
</tbody>
</table>
<br />
The performance of the tail-recursive version is quite OK relative to List and Array. Anyhow, the following results show the true power of tail recursion when it comes to 100 million runs.<br />
<br />
<table border="0" cellpadding="0" cellspacing="0" class="MsoNormalTable" style="border-collapse: collapse; margin: auto auto auto 4.65pt; mso-padding-alt: 0cm 5.4pt 0cm 5.4pt; mso-yfti-tbllook: 1184; width: 533px;">
<tbody>
<tr style="height: 15pt; mso-yfti-firstrow: yes; mso-yfti-irow: 0;"><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0) windowtext windowtext; border-style: solid none solid solid; border-width: 1pt 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">version<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">option price<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">CPU time (ms)<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: windowtext windowtext windowtext rgb(0, 0, 0); border-style: solid solid solid none; border-width: 1pt 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">absolute time (ms)<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 1;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">List<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;"><span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span><br />
<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;"><span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span><br />
<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 2;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Array<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;"><span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span><br />
<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;"><span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span><br />
<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;"><span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">out of memory</span></span><br />
<o:p></o:p></span></span></div>
</td></tr>
<tr style="height: 15pt; mso-yfti-irow: 3; mso-yfti-lastrow: yes;"><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext windowtext; border-style: none none solid solid; border-width: 0px 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Tail Recursion</span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118305<o:p></o:p></span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">26,800</span></span></div>
</td><td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext windowtext rgb(0, 0, 0); border-style: none solid solid none; border-width: 0px 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">26,898</span></span></div>
</td></tr>
</tbody>
</table>
<br />
Perhaps memory consumption is a more important factor to monitor than CPU time or absolute time. If I manually bring up Windows Task Manager, I can see that the memory use of the tail-recursive code remains steadily at the level of 1.4MB while the memory use of either the List or Array version keeps growing up to 1.8GB until it eventually throws out the out-of-memory exception and crashes.<br />
<br />
<strong>///////////////////////////////////////////////////////////////////////////////<br />open System<br />open System.Diagnostics</strong><br />
<br />
<strong>let N = 10000000<br />let rnd = new Random(1)</strong><br />
<br />
<strong>let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)</strong><br />
<br />
<strong>let call_payout K S = max (S-K) 0.0</strong><br />
<br />
<strong>let European_Option_Price (payout:float->float) S0 r vol T =<br /> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)<br /> let rec loop n sum =<br /> if n = N then sum<br /> else let new_sample = (rnd.NextDouble(),rnd.NextDouble()) <br /> |> box_muller_transform|> S_T |> payout<br /> loop (n+1) (sum+new_sample)<br /> (loop 0 0.0) / (float N) |> (*) (exp(-r*T))</strong><br />
<br />
<strong>let S0 = 50.0<br />let K = 40.0<br />let r = 0.01<br />let vol = 0.2<br />let T = 0.25<br />let price() = <br /> let price = European_Option_Price (call_payout K) S0 r vol T<br /> printfn "option_price = %f" price</strong><br />
<br />
<strong>let time f = <br /> let proc = Process.GetCurrentProcess()<br /> let cpu_time_stamp = proc.TotalProcessorTime<br /> let timer = new Stopwatch()<br /> timer.Start()<br /> try<br /> f()<br /> finally<br /> let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds<br /> printfn "CPU time = %dms" (int64 cpu_time)<br /> printfn "Absolute time = %dms" timer.ElapsedMilliseconds</strong><br />
<br />
<strong>time price</strong>Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com2tag:blogger.com,1999:blog-5063207632796600185.post-36031236030186639702012-08-14T20:03:00.001-07:002012-08-16T00:06:20.431-07:00F#: speed up Monte Carlo Simulation by reducing creation of a new listIn <a href="http://programmingcradle.blogspot.sg/2012/08/f-use-find-and-replace-to-speed-up.html">the previous post</a>, I learnt that creation of a new list could be much more costly than that of a new array. As such, I try to identify where in <a href="http://programmingcradle.blogspot.sg/2012/08/f-monte-carlo-pricer-for-european-calls.html">my original "List" version</a> I can reduce creation of a new list.<br />
<br />
Here is the original <b>European_Option_Price</b> function in the code:<br />
<br />
<b>let European_Option_Price (payout:float->float) S0 r vol T =</b><br />
<b> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)</b><br />
<b> randomNormals |> List.map S_T |> List.map payout |> List.average </b><br />
<b> |> (*) (exp(-r*T))</b><br />
<br />
I modify the third line by using an anonymous function to wrap up the two function calls to <b>S_T</b> and <b>payout</b>:<br />
<br />
<b>let European_Option_Price (payout:float->float) S0 r vol T =</b><br />
<b> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)</b><br />
<b> randomNormals |> List.map (fun x -> S_T x |> payout) |> List.average </b><br />
<b> |> (*) (exp(-r*T))</b><br />
<br />
The performance results are as follows:<br />
<br />
<table border="0" cellpadding="0" cellspacing="0" class="MsoNormalTable" style="border-collapse: collapse; margin: auto auto auto 4.65pt; mso-padding-alt: 0cm 5.4pt 0cm 5.4pt; mso-yfti-tbllook: 1184; width: 533px;">
<tbody>
<tr style="height: 15pt; mso-yfti-firstrow: yes; mso-yfti-irow: 0;">
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0) windowtext windowtext; border-style: solid none solid solid; border-width: 1pt 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">version</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">option
price</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">CPU
time (ms)</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext windowtext windowtext rgb(0, 0, 0); border-style: solid solid solid none; border-width: 1pt 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">absolute
time (ms)</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
</tr>
<tr style="height: 15pt; mso-yfti-irow: 1;">
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">List
(original)</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">10.118465</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">5,335</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">5,467</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
</tr>
<tr style="height: 15pt; mso-yfti-irow: 2; mso-yfti-lastrow: yes;">
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext windowtext; border-style: none none solid solid; border-width: 0px 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">List
(modified)</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">10.118465</span><span style="font-family: "Times New Roman","serif"; font-size: 12pt; mso-fareast-font-family: "Times New Roman";"><o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">3,993</span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext windowtext rgb(0, 0, 0); border-style: none solid solid none; border-width: 0px 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="font-family: Calibri;"><span style="color: black; font-size: 12pt; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;">4,035</span></span></div>
</td>
</tr>
</tbody></table>
<br />
Obviously, the improvement is significant (>25%), although it's not as big as that of "<a href="http://programmingcradle.blogspot.sg/2012/08/f-use-find-and-replace-to-speed-up.html">Find and Replace</a>". Nevertheless, if you like to stick with List in Monte Carlo simulation, perhaps you can consider creating less lists in your F# functions. <br />
<br />
By adding the anonymous function, we actually introduce some overhead because now we call <strong>S_T</strong> and <strong>payout</strong> indirectly and hence have more states to maintain on the call stack. If you like to avoid the overhead, you could choose to move the payout call to the end of the S_T function rather than using an anonymous function. <br />
<br />
<strong>let European_Option_Price (payout:float->float) S0 r vol T =<br /> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x) |> payout<br /> randomNormals |> List.map S_T |> List.average <br /> |> (*) (exp(-r*T))</strong><br />
<br />
You might see a bit improvment by saving the overhead incurred by the use of an anonymous function. However, I also tried that, but the improvement seems to be immaterial in this case so I don't bother to adopt it at this moment. <br />
<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-37688266115816238512012-08-13T21:16:00.000-07:002012-08-16T20:00:45.909-07:00F#: use "Find and Replace" to speed up Monte Carlo SimulationIn the <a href="http://programmingcradle.blogspot.sg/2012/08/f-monte-carlo-pricer-for-european-calls.html">previous post</a>, I implemented a Monte Carlo pricer using List in F#. If I would like to improve the performance of the pricer, what could I do? One may suggest the following solutions:<br />
<br />
<ul>
<li>(Computer Science Solution 1) Given that nowadays almost every machine (or even a mobile phone) is a multicore and there is no dependency between simulation runs, let's run them in parallel. </li>
<li>(Computer Science Solution 2) Profile the code, identify the bottleneck, figure out asymptotic time complexity, and then redesign the steps in the algorithm accordingly.</li>
<li>(Stochastic Process Solution) Apply variance reduction techniques so that I can achieve the same level of simulation accuracy with less runs.</li>
<li>(Numerical Solution?) pre-compute as much as possible, like <strong>(r - 0.5*(vol**2.0))*T</strong> or <strong>2.0*Math.PI</strong>...etc, so that I don't have to repeatedly compute the same constants.</li>
<li>(hmm... Finance Solution perhaps) Simply use money to resolve this issue, i.e., buy another machine with faster processors.</li>
</ul>
Well, all solutions above sound promising, but I actually took another even simpler approach. I just use the "Find and Replace" function in the Visual Studio code editor to replace all the "List" occurrences in my code with "Array" (and I tried "Seq" also). And then measure the performances of the 3 different versions: the original "List" version and the two new versions: "Array" and "Seq".<br />
<br />
The results are as follows:<br />
<table border="0" cellpadding="0" cellspacing="0" class="MsoNormalTable" style="border-collapse: collapse; margin: auto auto auto 4.65pt; mso-padding-alt: 0cm 5.4pt 0cm 5.4pt; mso-yfti-tbllook: 1184; width: 533px;">
<tbody>
<tr style="height: 15pt; mso-yfti-firstrow: yes; mso-yfti-irow: 0;">
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0) windowtext windowtext; border-style: solid none solid solid; border-width: 1pt 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">version<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">option price<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext rgb(0, 0, 0); border-style: solid none; border-width: 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">CPU time (ms)<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: windowtext windowtext windowtext rgb(0, 0, 0); border-style: solid solid solid none; border-width: 1pt 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; mso-border-top-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">absolute time (ms)<o:p></o:p></span></span></div>
</td>
</tr>
<tr style="height: 15pt; mso-yfti-irow: 1;">
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">List<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">5,335<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">5,467<o:p></o:p></span></span></div>
</td>
</tr>
<tr style="height: 15pt; mso-yfti-irow: 2;">
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none none solid; border-width: 0px 0px 0px 1pt; height: 15pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Array<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border: rgb(0, 0, 0); height: 15pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">1,606<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext rgb(0, 0, 0) rgb(0, 0, 0); border-style: none solid none none; border-width: 0px 1pt 0px 0px; height: 15pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">1,639<o:p></o:p></span></span></div>
</td>
</tr>
<tr style="height: 15pt; mso-yfti-irow: 3; mso-yfti-lastrow: yes;">
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext windowtext; border-style: none none solid solid; border-width: 0px 0px 1pt 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-left-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">Sequence<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">10.118465<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) rgb(0, 0, 0) windowtext; border-style: none none solid; border-width: 0px 0px 1pt; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">7,612<o:p></o:p></span></span></div>
</td>
<td nowrap="" style="background-color: transparent; border-color: rgb(0, 0, 0) windowtext windowtext rgb(0, 0, 0); border-style: none solid solid none; border-width: 0px 1pt 1pt 0px; height: 15pt; mso-border-bottom-alt: solid windowtext .5pt; mso-border-right-alt: solid windowtext .5pt; padding: 0cm 5.4pt; width: 100pt;" valign="bottom" width="133"><div align="center" class="MsoNormal" style="line-height: normal; margin: 0cm 0cm 0pt; text-align: center;">
<span style="color: black; mso-ascii-font-family: Calibri; mso-bidi-font-family: Calibri; mso-fareast-font-family: "Times New Roman"; mso-hansi-font-family: Calibri;"><span style="font-family: Calibri;">7,631<o:p></o:p></span></span></div>
</td>
</tr>
</tbody></table>
( Please note that optimization is turned on while compiling the 3 versions. The code of the Array version is attached at the bottome of this post.)<br />
<br />
All the 3 option prices are the same because they use the same seed, 1, to generate random numbers. As you can see, the performance improvement is obviously worth the effort of "Find and Replace". Now the new "Array" version can finish the same computations within less than 1/3 of the time needed by the original "List" version. <br />
<br />
The fact that Array is the fastest and Sequence is the slowest is not surprising to me. And I don't bother to find out why Sequence took so much time relatively at this moment, because the way I implemented the code doesn't really need its main strengths: lazy evaluation and ability to model infinite series. However, what I found intriguing is why Array is so fast relative to List. <br />
<br />
What I learnt from some F# introductory books is that Array is efficient for random access, which is also where List tends to perform horribly given its singly-linked nature. Though, in our Monte Carlo computations, we walk through a list or an array in a sequential manner only. As such, I don't think it's the efficiency of random access to cause the performance difference I'm looking at.<br />
<br />
After running a few simple tests on List and Array using F# Interactive as described below, it seems to me that in our case what really makes the difference is creation of a data structure rather than access to individual elements in the data structure. As you can see from below, it took more than 2 seconds to initialize a list of 10,000,000 floats, while the same action with an array took less than 0.1 seconds. Therefore, in terms of creation or allocation, arrays perform much faster than lists. In terms of iterating over the 10,000,000 floats, although the result below shows that the list was even slightly faster than the array, the difference is immaterial as far as our case is concerned.<br />
<br />
Please note that by this post I am not saying that arrays are always superior than lists. At least we should remember that lists have execellent constant time complexity in stack-style operations, which I don't use in my Monte Carlo code.<br />
<br />
Based on the finding that creation of a new list can be costly, I will be updating the code of the "List" version accordingly in next post.<br />
<br />
<br />//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br />
//// Here are the simple tests on Array and List using F# Interactive<br />
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br />
<u>Firstly, I turn on the timing function:</u><br />
<strong>> #time;;</strong><br />
<strong>--> Timing now on</strong><br />
<br />
<u>Secondly, I create a list of 10,000,000 floats and sum up using the fold function:</u><br />
<strong>> let X = List.init 10000000 (fun i -> float i);;<br />Real: 00:00:02.118, CPU: 00:00:02.246, GC gen0: 38, gen1: 20, gen2: 2</strong><br />
<br />
<strong>val X : float list =<br /> [0.0; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; 9.0; 10.0; 11.0; 12.0; 13.0;<br /> 14.0; 15.0; 16.0; 17.0; 18.0; 19.0; 20.0; 21.0; 22.0; 23.0; 24.0; 25.0;<br /> 26.0; 27.0; 28.0; 29.0; 30.0; 31.0; 32.0; 33.0; 34.0; 35.0; 36.0; 37.0;<br /> 38.0; 39.0; 40.0; 41.0; 42.0; 43.0; 44.0; 45.0; 46.0; 47.0; 48.0; 49.0;<br /> 50.0; 51.0; 52.0; 53.0; 54.0; 55.0; 56.0; 57.0; 58.0; 59.0; 60.0; 61.0;<br /> 62.0; 63.0; 64.0; 65.0; 66.0; 67.0; 68.0; 69.0; 70.0; 71.0; 72.0; 73.0;<br /> 74.0; 75.0; 76.0; 77.0; 78.0; 79.0; 80.0; 81.0; 82.0; 83.0; 84.0; 85.0;<br /> 86.0; 87.0; 88.0; 89.0; 90.0; 91.0; 92.0; 93.0; 94.0; 95.0; 96.0; 97.0;<br /> 98.0; 99.0; ...]</strong><br />
<br />
<strong>> let sum_X = List.fold (fun s x -> s + x) 0.0 X;;<br />Real: 00:00:00.062, CPU: 00:00:00.062, GC gen0: 0, gen1: 0, gen2: 0</strong><br />
<br />
<strong>val sum_X : float = 4.9999995e+13</strong><br />
<br />
<u>Finally, I do the same things with Array:</u><br />
<strong>> let Y = Array.init 10000000 (fun i -> float i);;<br />Real: 00:00:00.076, CPU: 00:00:00.078, GC gen0: 0, gen1: 0, gen2: 0</strong><br />
<br />
<strong>val Y : float [] =<br /> [|0.0; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; 9.0; 10.0; 11.0; 12.0; 13.0;<br /> 14.0; 15.0; 16.0; 17.0; 18.0; 19.0; 20.0; 21.0; 22.0; 23.0; 24.0; 25.0;<br /> 26.0; 27.0; 28.0; 29.0; 30.0; 31.0; 32.0; 33.0; 34.0; 35.0; 36.0; 37.0;<br /> 38.0; 39.0; 40.0; 41.0; 42.0; 43.0; 44.0; 45.0; 46.0; 47.0; 48.0; 49.0;<br /> 50.0; 51.0; 52.0; 53.0; 54.0; 55.0; 56.0; 57.0; 58.0; 59.0; 60.0; 61.0;<br /> 62.0; 63.0; 64.0; 65.0; 66.0; 67.0; 68.0; 69.0; 70.0; 71.0; 72.0; 73.0;<br /> 74.0; 75.0; 76.0; 77.0; 78.0; 79.0; 80.0; 81.0; 82.0; 83.0; 84.0; 85.0;<br /> 86.0; 87.0; 88.0; 89.0; 90.0; 91.0; 92.0; 93.0; 94.0; 95.0; 96.0; 97.0;<br /> 98.0; 99.0; ...|]</strong><br />
<br />
<strong>> let sum_Y = Array.fold (fun s x -> s + x) 0.0 Y;;<br />Real: 00:00:00.107, CPU: 00:00:00.109, GC gen0: 0, gen1: 0, gen2: 0</strong><br />
<br />
<strong>val sum_Y : float = 4.9999995e+13</strong><br />
<br />
<br />
<br />//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br />
//// The Array version of Monte Calro Simulation<br />//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br />
<strong>open System<br />open System.Diagnostics</strong><br />
<br />
<strong>let N = 10000000<br />let rnd = new Random(1) <br />let randomPairs = Array.init N (fun i -> (rnd.NextDouble(), rnd.NextDouble()))<br />let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)<br />let randomNormals = randomPairs |> Array.map box_muller_transform</strong><br />
<br />
<strong>let call_payout K S = max (S-K) 0.0</strong><br />
<br />
<strong>let European_Option_Price (payout:float->float) S0 r vol T =<br /> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)<br /> randomNormals |> Array.map S_T |> Array.map payout |> Array.average <br /> |> (*) (exp(-r*T))</strong><br />
<br />
<strong>let S0 = 50.0<br />let K = 40.0<br />let r = 0.01<br />let vol = 0.2<br />let T = 0.25<br />let price() = <br /> let price = European_Option_Price (call_payout K) S0 r vol T<br /> printfn "option_price = %f" price</strong><br />
<br />
<strong>let time f = <br /> let proc = Process.GetCurrentProcess()<br /> let cpu_time_stamp = proc.TotalProcessorTime<br /> let timer = new Stopwatch()<br /> timer.Start()<br /> try<br /> f()<br /> finally<br /> let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds<br /> printfn "CPU time = %dms" (int64 cpu_time)<br /> printfn "Absolute time = %dms" timer.ElapsedMilliseconds</strong><br />
<br />
<strong>time price</strong><br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com2tag:blogger.com,1999:blog-5063207632796600185.post-25613930671283034122012-08-09T03:55:00.001-07:002012-08-09T03:55:09.791-07:00F#: a Monte Carlo Pricer for European Calls and PutsHere we implement Monte Carlo simulation for pricing European Calls and Puts in a Black-Scholes world. Although we already have closed-form solutions, I think it's a good, simple but non-trivial example to show how Monte Carlo simulation can be nicely and concisely done in F#, thanks to the following functional programming techniques:<br />
<br />
<ul>
<li>List processing</li>
<li>Higher-order function</li>
<li>Function currying</li>
<li>Function pipelining</li>
</ul>
<br />
There are a few things to note about the code example:<br />
<br />
<ul>
<li>Normal random numbers are generated using Box-Muller transform.</li>
<li>The number of simulation runs is 10,000,000.</li>
<li><a href="http://programmingcradle.blogspot.sg/2012/08/f-profiling-how-to-measure-cpu-time-and.html">The time function</a> is used to measure CPU time and absolute time.</li>
<li>The seed to initiate System.Random is deliberately fixed as 1 so that I can get the same series of normal random numbers every time I run the code. If you don't want this, please remember to remove the seed.</li>
<li>The underlying dynamics and European call option details are as follows:</li>
<ul>
<li>S0 (spot price) = 50.0</li>
<li>vol (annualized volatility) = 20%</li>
<li>risk free rate = 1%</li>
<li>strike = 40.0</li>
<li>expiration (in years) = 0.25</li>
</ul>
<li>As a sanity check, the closed-form
theoretical price based on the parameters above is 10.11842.</li>
</ul>
After running the code on my 4-year-old laptop, I got the following output :<br />
<br />
<b>option_price = 10.118465</b><br />
<b>CPU time = 5709ms</b><br />
<b>Absolute time = 7115ms</b><br />
<br />
It took about 7 seconds to run 10 million simulations. And the simulated price agrees with the theoretical price up to 4 decimal places.<br />
<br />
Finally, here is the code:<br />
<b>///////////////////////////////////////////////////////////////////////////////</b><br />
<br />
<b>open System</b><br />
<b>open System.Diagnostics</b><br />
<b><br /></b>
<b>let N = 10000000</b><br />
<b>let rnd = new Random(1)</b><br />
<b>let randomPairs = List.init N (fun i -> (rnd.NextDouble(), rnd.NextDouble()))</b><br />
<b>let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)</b><br />
<b>let randomNormals = randomPairs |> List.map box_muller_transform</b><br />
<b><br /></b>
<b>let call_payout K S = max (S-K) 0.0</b><br />
<b><br /></b>
<b>let European_Option_Price (payout:float->float) S0 r vol T =</b><br />
<b> let S_T x = S0 * exp((r - 0.5*(vol**2.0))*T + vol*sqrt(T)*x)</b><br />
<b> randomNormals |> List.map S_T |> List.map payout |> List.average </b><br />
<b> |> (*) (exp(-r*T))</b><br />
<b><br /></b>
<b>let S0 = 50.0</b><br />
<b>let K = 40.0</b><br />
<b>let r = 0.01</b><br />
<b>let vol = 0.2</b><br />
<b>let T = 0.25</b><br />
<b>let price() = </b><br />
<b> let price = European_Option_Price (call_payout K) S0 r vol T</b><br />
<b> printfn "option_price = %f" price</b><br />
<b><br /></b>
<b>let time f = </b><br />
<b> let proc = Process.GetCurrentProcess()</b><br />
<b> let cpu_time_stamp = proc.TotalProcessorTime</b><br />
<b> let timer = new Stopwatch()</b><br />
<b> timer.Start()</b><br />
<b> try</b><br />
<b> f()</b><br />
<b> finally</b><br />
<b> let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds</b><br />
<b> printfn "CPU time = %dms" (int64 cpu_time)</b><br />
<b> printfn "Absolute time = %dms" timer.ElapsedMilliseconds</b><br />
<b><br /></b>
<b>time price</b><br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com0tag:blogger.com,1999:blog-5063207632796600185.post-62164691781395839922012-08-06T22:05:00.000-07:002012-08-16T00:48:13.862-07:00F# Profiling: How to measure CPU time and absolute time?Jon Harrop wrote a nice book, "F# for Scientists." In Chapter 8 of the book, he explained how to time the execution of a function in F#. Here I simply summarized what he said in the book and I revised his F# code examples to fit my needs.<br />
<br />
Basically I'm interested in measuring the following two kinds of time:<br />
<br />
<ul>
<li>Absolute time (or wall-clock time): the time elapsed in the real world.</li>
<li>CPU time: the amount of time that CPUs have spent running the function.</li>
</ul>
<div>
<span style="font-family: inherit;">The .Net class <b>System.Diagnostics.Stopwatch</b> can be used to measure absolute time. And the <b>TotalProcessorTime</b> property of <b>System.Diagnostics.Process</b> returns a <b>TimeSpan</b> object, which indicates the amount of CPU time that the associated process has spent running the F# function.</span></div>
<div>
<span style="font-family: inherit;"><b><br /></b></span></div>
<div>
<div>
<b>///////////////////////////////////////////////////////////</b></div>
<div>
<br />
<br />
<b>open System.Diagnostics</b><br />
<b><br /></b>
<b>let time f = </b><br />
<b> let proc = Process.GetCurrentProcess()</b><br />
<b> let cpu_time_stamp = proc.TotalProcessorTime</b><br />
<b> let timer = new Stopwatch()</b><br />
<b> timer.Start()</b><br />
<b> try</b><br />
<b> f()</b><br />
<b> finally</b><br />
<b> let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds</b><br />
<b> printfn "CPU time = %dms" (int64 cpu_time)</b><br />
<b> printfn "Absolute time = %dms" timer.ElapsedMilliseconds</b><br />
<b><br /></b>
<b>let rec loop n f x =</b><br />
<b> if n > 0 then</b><br />
<b> f x |> ignore</b><br />
<b> loop (n-1) f x</b><br />
<b><br /></b>
<b>time (fun () -> loop 1000000 List.sum [1..100])</b></div>
</div>
<div>
<br /></div>
<div>
</div>
Chao-Jen Chenhttp://www.blogger.com/profile/17567134700618641403noreply@blogger.com2