VBA5 - Monte Carlo Integration

To find the area under a curve, one can use integral calculus. If the curve has no close form, such as the normal curve, then the area can not be derived analytically. However, with today's computer technology, one can use Monte Carlo Integration to achieve such task. The area under a distribution is also known as probability. In this example, we want to compute the area under standard normal probability distribution from 0 to z. We will set z equal to 1.8.

The Monte Carlo Integration procedure is as follow:

1. Identify the range of X and Y coordinates where the random number will be placed.
(The random numbers that we use follow a uniform distribution.) Ex: In this case, the minimum coordinate of X and Y are both zero. The maximum coordinates of X are given by the user. The maximum coordinate of Y occurs where the X coordinate is zero. In this case, it is 1 over the square root of 2 pi, 0.39894.

2. Compute the area of the rectangle in question (defined by red dots in Figure 1) using the X and Y range.

3. Run the random process. All the random numbers (X,Y) will land within the rectangle. Count how many points land below the curve.

4. Divide the sum of points below the curve by the number of iterations to get the proportion of points below the curve to the total number of points.

5. Multiply this proportion by the area to get the probability. This probability is shown on cell D4 of Figure 2. The probability is 46.37%. Compare that result to 46.41%, the result that can be found in the z Table in virtually any statistics textbook.

The second program utilizes the random numbers generated from a normal distribution to derive the probability. The probability under 1,000,000 iterations using this method is 46.43%. In term of computer time required, this method is 3 to 5 times faster than the Monte Carlo Integration method. However, in order to use the latter method, one must be capable of generating the random numbers which associated to a specific distribution.

Option Explicit
Option Base 1

'***********************************************************************
'*     Obtain area under the normal distribution curve (probability)   *
'*                          using Monte Carlo Integration Method                  *
'***********************************************************************

Sub Integrate_Normal()
    Dim n As Variant
    Dim mean As Variant
    Dim sd As Variant
    Dim x As Variant
    Dim MaxY As Single, z As Single, f_of_z As Single
    Dim Proportion As Single, Area As Single, Prob As Single, y As Single
    Dim i As Long, count As Long

    Randomize
    n = Range("c4").Value
    mean = Range("c5").Value
    sd = Range("c6").Value
    x = Range("c7").Value
    count = 0
    MaxY = 1 / Sqr(2 * Application.Pi())

    For i = 1 To n
        z = Rnd * (x - mean) / sd
        y = Rnd * MaxY
        f_of_z = 1 / Sqr(2 * Application.Pi()) * Exp(-z ^ 2 / 2)
        If (y < f_of_z) Then count = count + 1
    Next i

    Proportion = count / n
    Area = MaxY * (x - mean) / sd
    Prob = Proportion * Area
    If (Prob > 1) Then Prob = 1
    Cells(4, 4) = Prob

End Sub

'***********************************************************************
'*     Obtain area under the normal distribution curve (probability)   *
'*      using random numbers from Standard Normal Distribution   *
'***********************************************************************

Sub NormNo()
    Dim n As Variant
    Dim mean As Variant
    Dim sd As Variant
    Dim MaxX As Variant
    Dim MinX As Variant
    Dim i As Long, count As Long
    Dim y As Double

    n = Range("c12").Value
    mean = Range("c13").Value
    sd = Range("c14").Value
    MaxX = Range("c15").Value
    MinX = mean
    Randomize
    count = 0

    For i = 1 To n
        y = gauss * sd + mean
        If y >= MinX And y <= MaxX Then count = count + 1
    Next i

    Cells(12, 4) = count / n

End Sub

'***********************************************************************
'*   Return random numbers from Standard Normal Distribution   *
'***********************************************************************

Function gauss()
    Dim fac As Double, r As Double, V1 As Double, V2 As Double

10      V1 = 2 * Rnd - 1
    V2 = 2 * Rnd - 1
    r = V1 ^ 2 + V2 ^ 2
    If (r >= 1) Then GoTo 10
    fac = Sqr(-2 * Log(r) / r)
    gauss = V2 * fac

End Function