F#: How could the Monte Calro pricer return option_price = infinity?

If you leave the random seed blank in the Monte Carlo code and let 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

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 .

Well, there's nothing wrong for 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

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

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 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
        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."

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
time()