**System.Random**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:

**option_price = Infinity**

CPU time = 2636ms

Absolute time = 2435ms

CPU time = 2636ms

Absolute time = 2435ms

Why you see this funny "infinity" price is because:

**System.Random.NextDouble()**could generate 0.0, though the probability is extremely low. As such, it's possible to generate a tuple of**(u1,u2)**where**u1 = 0.0**and**(u2 < 0.25 or u2 > 0.75)**. This kind of**(u1,u2)**drives the Box Muller transform to infinity. Then, not surprisingly, we can see that a single sample of infinite payout screws up the Law of Large Numbers, as the integrability condition is not satisfied.- The notion of infinity is defined in F# (or actually via
**Double.PositiveInfinity**and**Double.NegativeInfinity**in .Net). F# knows how to handle**log(0.0)**and**exp(infinity)**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 .

**System.Random**to generate a value of 0.0, because this is simply how it's designed. The

**NextDouble()**method is meant to approximate

**U[0,1)**, but the possibility of getting a 0.0 value is not desirable here. One workaround is to take a simple transform:

**Y = 1.0 - U[0,1)**. Then

**Y**follows

**U(0,1]**. So the code snippet we use to draw random numbers can be modified as follows:

**let N = 10000000**

let rnd = new Random(388)

let next_double () = 1.0 - rnd.NextDouble()

let randomPairs = Array.init N (fun i -> (next_double(), next_double()))

let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)

let randomNormals = randomPairs |> Array.map box_muller_transform

let rnd = new Random(388)

let next_double () = 1.0 - rnd.NextDouble()

let randomPairs = Array.init N (fun i -> (next_double(), next_double()))

let box_muller_transform (u1,u2) = sqrt(-2.0*log(u1)) * cos(2.0*Math.PI*u2)

let randomNormals = randomPairs |> Array.map box_muller_transform

You might argue that mathematically

**P(U=0)**is zero in either

**U(0,1)**or

**U[0,1)**or

**U(0,1]**or

**U[0,1]**. When it comes to implementation of a pseudo-random number generator, at least in the case of

**System.Random**,

**U[0,1)**usually implies

**P(U=0)>0**.

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:

**///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////**

**open System**

open System.Diagnostics

open System.Diagnostics

**let search_for_0 seed_min seed_max draw_max =**

let found u1 u2 = if (u1 = 0.0 && (u2 < 0.25 || u2 > 0.75)) then true else false

let found u1 u2 = if (u1 = 0.0 && (u2 < 0.25 || u2 > 0.75)) then true else false

**let rec seed_loop seed =**

let rnd = new Random(seed)

let rec draw_loop draw =

let (u1,u2) = (rnd.NextDouble(),rnd.NextDouble())

if found u1 u2 || draw = draw_max then (u1,u2, draw)

else draw_loop (draw + 1)

let (u1,u2, draw) = draw_loop 1

if found u1 u2 || seed = seed_max then (u1,u2, draw, seed)

else seed_loop (seed + 1)

let (u1,u2, draw, seed) = seed_loop seed_min

printfn "seed_min = %d, seed_max = %d, draw_max = %d" seed_min seed_max draw_max

if found u1 u2 then

let rnd = new Random(seed)

let rec draw_loop draw =

let (u1,u2) = (rnd.NextDouble(),rnd.NextDouble())

if found u1 u2 || draw = draw_max then (u1,u2, draw)

else draw_loop (draw + 1)

let (u1,u2, draw) = draw_loop 1

if found u1 u2 || seed = seed_max then (u1,u2, draw, seed)

else seed_loop (seed + 1)

let (u1,u2, draw, seed) = seed_loop seed_min

printfn "seed_min = %d, seed_max = %d, draw_max = %d" seed_min seed_max draw_max

if found u1 u2 then

**printfn "found u1 = 0.0 and u2 = %f at draw %d of seed %d" u2 draw seed**

else printfn "couldn't find 0.0 within the range specified."

else printfn "couldn't find 0.0 within the range specified."

**let time () =**

let proc = Process.GetCurrentProcess()

let cpu_time_stamp = proc.TotalProcessorTime

let timer = new Stopwatch()

timer.Start()

try

search_for_0 1 1000 10000000

finally

let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds

printfn "CPU time = %dms" (int64 cpu_time)

printfn "Absolute time = %dms" timer.ElapsedMilliseconds

let proc = Process.GetCurrentProcess()

let cpu_time_stamp = proc.TotalProcessorTime

let timer = new Stopwatch()

timer.Start()

try

search_for_0 1 1000 10000000

finally

let cpu_time = (proc.TotalProcessorTime-cpu_time_stamp).TotalMilliseconds

printfn "CPU time = %dms" (int64 cpu_time)

printfn "Absolute time = %dms" timer.ElapsedMilliseconds

**time()**

## No comments:

## Post a Comment