Attribute VB_Name = "FINAN_DERIV_BS_SIMULATION_LIBR"


Option Explicit     'Requires that all variables to be declared explicitly.
Option Base 1       'The "Option Base" statement allows to specify 0 or 1 as the
                    'default first index of arrays.

'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : CC_PP_LC_SP_SIMULATION_FUNC
'DESCRIPTION   : Simulation of Optioned Portfolios: Return simulation for Covered
'Call, Protective Put, Long Call and Short Put portfolios and Kernel
'density estimation (non-parametric estimation of the probability
'distribution function).

'LIBRARY       : DERIVATIVES
'GROUP         : PORTFOLIO
'ID            : 001
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************

Function CC_PP_LC_SP_SIMULATION_FUNC(ByVal STRIKE As Double, _
ByVal EXPIRATION As Double, _
ByVal INITIAL_UNDERLYING_PRICE As Double, _
ByVal UNDERLYING_VOLATILITY As Double, _
ByVal RISK_FREE_RATE As Double, _
ByVal UNDERLYING_EXPECTED_RETURN As Double, _
Optional ByVal OPTION_WEIGHT As Double = 0.5, _
Optional ByVal NBINS As Long = 100, _
Optional ByVal nLOOPS As Long = 1000, _
Optional ByVal NTRIALS As Long = 1, _
Optional ByVal VERSION As Integer = 4, _
Optional ByVal OUTPUT As Integer = 1, _
Optional ByVal CND_TYPE As Integer = 0)

'Strike: 100
'Expiration: 1
'Initial Underlying Price: 100
'Underlying Volatility: 15%
'Riskfree Rate: 10%
'Expected Return of Underlying: 3%
'% of Portfolio which is optioned: 50%

Dim i As Long
Dim j As Long
Dim DATA_MATRIX As Variant
'% of Portfolio which is optioned

On Error GoTo ERROR_LABEL

'----------------------------------------------------------------------
DATA_MATRIX = MATRIX_RANDOM_NORMAL_FUNC(nLOOPS, NTRIALS, 0, _
              UNDERLYING_EXPECTED_RETURN, UNDERLYING_VOLATILITY, 0)
'----------------------------------------------------------------------
Select Case VERSION
'----------------------------------------------------------------------
Case 0 'Long Underlying
'----------------------------------------------------------------------
    
    For j = 1 To NTRIALS
        For i = 1 To nLOOPS
            DATA_MATRIX(i, j) = INITIAL_UNDERLYING_PRICE * _
                                Exp(DATA_MATRIX(i, j) * EXPIRATION)
            DATA_MATRIX(i, j) = Log(DATA_MATRIX(i, j) / INITIAL_UNDERLYING_PRICE)
        Next i
    Next j
'----------------------------------------------------------------------
Case 1 'Short Call + Underlying
'----------------------------------------------------------------------
    For j = 1 To NTRIALS
        For i = 1 To nLOOPS
            DATA_MATRIX(i, j) = INITIAL_UNDERLYING_PRICE * _
                                Exp(DATA_MATRIX(i, j) * EXPIRATION)
            
            DATA_MATRIX(i, j) = OPTION_WEIGHT * Log((-MAXIMUM_FUNC( _
                                DATA_MATRIX(i, j) - STRIKE, 0) + _
                                DATA_MATRIX(i, j) + EUROPEAN_CALL_OPTION_FUNC( _
                                INITIAL_UNDERLYING_PRICE, STRIKE, EXPIRATION, _
                                RISK_FREE_RATE, 0, UNDERLYING_VOLATILITY, _
                                CND_TYPE) * Exp(EXPIRATION * RISK_FREE_RATE)) / _
                                INITIAL_UNDERLYING_PRICE) + (1 - OPTION_WEIGHT) * _
                                Log(DATA_MATRIX(i, j) / INITIAL_UNDERLYING_PRICE)
        Next i
    Next j
'----------------------------------------------------------------------
Case 2 'Long Put + Underlying
'----------------------------------------------------------------------
    For j = 1 To NTRIALS
        For i = 1 To nLOOPS
            DATA_MATRIX(i, j) = INITIAL_UNDERLYING_PRICE * _
                                Exp(DATA_MATRIX(i, j) * EXPIRATION)
            
            DATA_MATRIX(i, j) = OPTION_WEIGHT * Log((MAXIMUM_FUNC(STRIKE - _
                                DATA_MATRIX(i, j), 0) + DATA_MATRIX(i, j)) / _
                                (INITIAL_UNDERLYING_PRICE + _
                                EUROPEAN_PUT_OPTION_FUNC(INITIAL_UNDERLYING_PRICE, _
                                STRIKE, EXPIRATION, RISK_FREE_RATE, 0, _
                                UNDERLYING_VOLATILITY, CND_TYPE))) + _
                                (1 - OPTION_WEIGHT) * Log(DATA_MATRIX(i, j) / _
                                INITIAL_UNDERLYING_PRICE)
        Next i
    Next j
'----------------------------------------------------------------------
Case 3 'Short Put + Underlying
'----------------------------------------------------------------------
    For j = 1 To NTRIALS
        For i = 1 To nLOOPS
            DATA_MATRIX(i, j) = INITIAL_UNDERLYING_PRICE * _
                                Exp(DATA_MATRIX(i, j) * EXPIRATION)
            DATA_MATRIX(i, j) = OPTION_WEIGHT * Log((-MAXIMUM_FUNC(STRIKE - _
                                DATA_MATRIX(i, j), 0) + DATA_MATRIX(i, j) + _
                                EUROPEAN_PUT_OPTION_FUNC(INITIAL_UNDERLYING_PRICE, _
                                STRIKE, EXPIRATION, RISK_FREE_RATE, 0, _
                                UNDERLYING_VOLATILITY, CND_TYPE) * Exp(EXPIRATION * _
                                RISK_FREE_RATE)) / INITIAL_UNDERLYING_PRICE) _
                                + (1 - OPTION_WEIGHT) * Log(DATA_MATRIX(i, j) / _
                                INITIAL_UNDERLYING_PRICE)
        Next i
    Next j
'----------------------------------------------------------------------
Case Else 'Long Call + Underlying
'----------------------------------------------------------------------
    For j = 1 To NTRIALS
        For i = 1 To nLOOPS
            DATA_MATRIX(i, j) = INITIAL_UNDERLYING_PRICE * _
                                Exp(DATA_MATRIX(i, j) * EXPIRATION)
            DATA_MATRIX(i, j) = OPTION_WEIGHT * Log((MAXIMUM_FUNC( _
                                DATA_MATRIX(i, j) - STRIKE, 0) + _
                                DATA_MATRIX(i, j)) / (INITIAL_UNDERLYING_PRICE _
                                + EUROPEAN_CALL_OPTION_FUNC( _
                                INITIAL_UNDERLYING_PRICE, STRIKE, _
                                EXPIRATION, RISK_FREE_RATE, 0, _
                                UNDERLYING_VOLATILITY, CND_TYPE))) + _
                                (1 - OPTION_WEIGHT) * _
                                Log(DATA_MATRIX(i, j) / INITIAL_UNDERLYING_PRICE)
        Next i
    Next j
'----------------------------------------------------------------------
End Select
'----------------------------------------------------------------------

'----------------------------------------------------------------------
Select Case OUTPUT
'----------------------------------------------------------------------
Case 0
'----------------------------------------------------------------------
    CC_PP_LC_SP_SIMULATION_FUNC = DATA_MATRIX
'----------------------------------------------------------------------
Case 1
'----------------------------------------------------------------------
    CC_PP_LC_SP_SIMULATION_FUNC = _
        ADVANCED_MOMENTS_FUNC(DATA_MATRIX, 0, 0, 0.05, 0)
'----------------------------------------------------------------------
Case Else
'----------------------------------------------------------------------
    'Probability distributions for Returns of Optioned Portfolios
    CC_PP_LC_SP_SIMULATION_FUNC = _
        HISTOGRAM_KERNEL_DIST_FUNC(DATA_MATRIX, NBINS) 'X & PDF
'----------------------------------------------------------------------
End Select
'----------------------------------------------------------------------
Exit Function
ERROR_LABEL:
CC_PP_LC_SP_SIMULATION_FUNC = Err.number
End Function



'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : ANTITHETIC_VARIABLES_CALL_OPTION_SIMULATION_FUNC
'DESCRIPTION   : Black Scholes Call price by Monte Carlo with Antithetic variables
'LIBRARY       : DERIVATIVES
'GROUP         : BS_MONTE_CARLO
'ID            : 001
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************

Function ANTITHETIC_VARIABLES_CALL_OPTION_SIMULATION_FUNC(ByVal SPOT As Double, _
ByVal STRIKE As Double, _
ByVal EXPIRATION As Double, _
ByVal RATE As Double, _
ByVal SIGMA As Double, _
Optional ByVal nLOOPS As Long = 1000, _
Optional ByVal OUTPUT As Integer = 0)

Dim i As Long  'paths counter

Dim FIRST_VAL As Double   ' stock values using z
Dim SECOND_VAL As Double   ' stock values using -z

Dim NORM_VAL As Double    'standard normal value
Dim FIRST_OPT_VAL As Double    'Black Scholes with z
Dim SECOND_OPT_VAL As Double    'Black Scholes with -z

Dim TEMP_MATRIX As Variant
Dim RND_NORMAL_MATRIX As Variant

On Error GoTo ERROR_LABEL

'Beginning the summation
FIRST_OPT_VAL = 0
SECOND_OPT_VAL = 0

If OUTPUT <> 0 Then: ReDim TEMP_MATRIX(1 To nLOOPS, 1 To 2)

RND_NORMAL_MATRIX = MATRIX_RANDOM_NORMAL_FUNC(nLOOPS, 1, 0, 0, 1, 0)

For i = 1 To nLOOPS 'Loop over M paths
    NORM_VAL = RND_NORMAL_MATRIX(i, 1)
    FIRST_VAL = SPOT * Exp((RATE - ((SIGMA ^ 2) / 2)) * _
                EXPIRATION + SIGMA * (EXPIRATION ^ 0.5) * NORM_VAL)
    SECOND_VAL = SPOT * Exp((RATE - ((SIGMA ^ 2) / 2)) * _
                EXPIRATION - SIGMA * (EXPIRATION ^ 0.5) * NORM_VAL)
    FIRST_OPT_VAL = FIRST_OPT_VAL + MAXIMUM_FUNC(FIRST_VAL - STRIKE, 0)
    SECOND_OPT_VAL = SECOND_OPT_VAL + MAXIMUM_FUNC(SECOND_VAL - STRIKE, 0)
    
    If OUTPUT <> 0 Then
        TEMP_MATRIX(i, 1) = FIRST_OPT_VAL
        TEMP_MATRIX(i, 2) = SECOND_OPT_VAL
    End If
Next i

Select Case OUTPUT
    Case 0
        ANTITHETIC_VARIABLES_CALL_OPTION_SIMULATION_FUNC = Exp(-RATE * EXPIRATION) * (FIRST_OPT_VAL + _
                SECOND_OPT_VAL) / (2 * nLOOPS)
    Case Else
        ANTITHETIC_VARIABLES_CALL_OPTION_SIMULATION_FUNC = TEMP_MATRIX
End Select

Exit Function
ERROR_LABEL:
ANTITHETIC_VARIABLES_CALL_OPTION_SIMULATION_FUNC = Err.number
End Function


'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : SDE_CALL_OPTION_SIMULATION_FUNC
'DESCRIPTION   : Black Scholes Call price by Monte Carlo with Antithetic
'variables and SDE discretization
'LIBRARY       : DERIVATIVES
'GROUP         : BS_MONTE_CARLO
'ID            : 002
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************


Function SDE_CALL_OPTION_SIMULATION_FUNC(ByVal SPOT As Double, _
ByVal STRIKE As Double, _
ByVal EXPIRATION As Double, _
ByVal RATE As Double, _
ByVal SIGMA As Double, _
Optional ByVal nLOOPS As Long = 1000, _
Optional ByVal nSTEPS As Long = 20, _
Optional ByVal OUTPUT As Integer = 0)

Dim i As Long  'paths counter
Dim j As Long  'time periods counter

Dim FIRST_VAL As Double   ' stock values using z
Dim SECOND_VAL As Double   ' stock values using -z

Dim NORM_VAL As Double   'standard normal value
Dim DELTA_VAL As Double  'time step

Dim FIRST_OPT_VAL As Double    'Black Scholes with z
Dim SECOND_OPT_VAL As Double    'Black Scholes with -z

Dim TEMP_MATRIX As Variant
Dim RND_NORMAL_MATRIX As Variant

On Error GoTo ERROR_LABEL

'Beginning the summation
FIRST_OPT_VAL = 0
SECOND_OPT_VAL = 0
DELTA_VAL = EXPIRATION / nSTEPS

If OUTPUT <> 0 Then: ReDim TEMP_MATRIX(1 To nLOOPS, 1 To nSTEPS, 1 To 2)
RND_NORMAL_MATRIX = MATRIX_RANDOM_NORMAL_FUNC(nLOOPS, nSTEPS, 0, 0, 1, 0)

For i = 1 To nLOOPS 'Loop over M paths
    FIRST_VAL = SPOT
    SECOND_VAL = SPOT
    For j = 1 To nSTEPS 'Loop over N time steps
        NORM_VAL = RND_NORMAL_MATRIX(i, j)
        FIRST_VAL = FIRST_VAL + FIRST_VAL * (RATE * DELTA_VAL + _
                SIGMA * (DELTA_VAL ^ 0.5) * NORM_VAL)
        SECOND_VAL = SECOND_VAL + SECOND_VAL * (RATE * DELTA_VAL - _
                SIGMA * (DELTA_VAL ^ 0.5) * NORM_VAL)
        
        If OUTPUT <> 0 Then
            TEMP_MATRIX(i, j, 1) = FIRST_VAL
            TEMP_MATRIX(i, j, 2) = SECOND_VAL
        End If
    Next j
FIRST_OPT_VAL = FIRST_OPT_VAL + MAXIMUM_FUNC(FIRST_VAL - STRIKE, 0)
SECOND_OPT_VAL = SECOND_OPT_VAL + MAXIMUM_FUNC(SECOND_VAL - STRIKE, 0)
Next i

Select Case OUTPUT
    Case 0
        SDE_CALL_OPTION_SIMULATION_FUNC = Exp(-RATE * EXPIRATION) * (FIRST_OPT_VAL + _
            SECOND_OPT_VAL) / (2 * nLOOPS)
    Case Else
        SDE_CALL_OPTION_SIMULATION_FUNC = TEMP_MATRIX
End Select

Exit Function
ERROR_LABEL:
    SDE_CALL_OPTION_SIMULATION_FUNC = Err.number
End Function


'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : EUROPEAN_OPTION_SIMULATION_FUNC
'DESCRIPTION   : EUROPEAN OPTION MC SIMULATION
'LIBRARY       : DERIVATIVES
'GROUP         : BS_MONTE_CARLO
'ID            : 003
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************

Function EUROPEAN_OPTION_SIMULATION_FUNC(ByVal nLOOPS As Long, _
ByVal SPOT As Double, _
ByVal STRIKE As Double, _
ByVal EXPIRATION As Double, _
ByVal RATE As Double, _
ByVal SIGMA As Double, _
Optional ByVal OPTION_FLAG As Integer = 1)

'OPTION_FLAG = 1 --> CALL OPTION
'OPTION_FLAG = -1 --> PUT OPTION

Dim i As Long
Dim PAYOFF As Double
Dim TEMP_RND As Double
    
Dim TEMP_SAMPLE As Double
    
Dim CUMUL_PAYOFF As Double
Dim CUMUL_PAYOFF_SQR As Double
    
Dim TEMP_MATRIX As Variant

On Error GoTo ERROR_LABEL
    
If OPTION_FLAG <> 1 Then OPTION_FLAG = -1
'the default is a call option, unless OPTION_FLAG = -1 to indicate put option

CUMUL_PAYOFF = 0#
CUMUL_PAYOFF_SQR = 0#
CUMUL_PAYOFF = 0
For i = 1 To nLOOPS
    TEMP_RND = RANDOM_NORMAL_FUNC(0, 1, 0)
    TEMP_SAMPLE = SPOT * Exp((RATE - 0.5 * SIGMA ^ 2) * EXPIRATION _
                            + SIGMA * Sqr(EXPIRATION) * TEMP_RND)
                            'Geometric Brownian Motion
    PAYOFF = MAXIMUM_FUNC(0#, OPTION_FLAG * TEMP_SAMPLE - OPTION_FLAG * STRIKE)
    CUMUL_PAYOFF = CUMUL_PAYOFF + PAYOFF
    CUMUL_PAYOFF_SQR = CUMUL_PAYOFF_SQR + PAYOFF ^ 2
Next i
      
ReDim TEMP_MATRIX(1 To 2, 1 To 2)
      
TEMP_MATRIX(1, 1) = "OPT VALUE"
TEMP_MATRIX(2, 1) = "STD. ERROR"
      
TEMP_MATRIX(1, 2) = Exp(-RATE * EXPIRATION) * CUMUL_PAYOFF / nLOOPS
TEMP_MATRIX(2, 2) = (Sqr((CUMUL_PAYOFF_SQR - CUMUL_PAYOFF ^ 2 / nLOOPS) / _
                    nLOOPS - 1)) / Sqr(nLOOPS)
      
EUROPEAN_OPTION_SIMULATION_FUNC = TEMP_MATRIX

Exit Function
ERROR_LABEL:
EUROPEAN_OPTION_SIMULATION_FUNC = Err.number
End Function

'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : EUROPEAN_CALL_CORPUT_SEQUENCE_SIMULATION_FUNC
'DESCRIPTION   : Corput sequence technique for European Call Option Simulation
'LIBRARY       : DERIVATIVES
'GROUP         : BS_MONTE_CARLO
'ID            : 004
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************

Function EUROPEAN_CALL_CORPUT_SEQUENCE_SIMULATION_FUNC(ByVal SPOT As Double, _
ByVal STRIKE As Double, _
ByVal EXPIRATION As Double, _
ByVal RATE As Double, _
ByVal DIVD As Double, _
ByVal SIGMA As Double, _
Optional ByVal nLOOPS As Long = 1000)
     
Dim i As Long
     
Dim TEMP_SUM As Double
Dim ATEMP_VAL As Double
Dim BTEMP_VAL As Double
     
Dim INV_VALUE As Double
Dim SIM_VALUE As Double
Dim PROB_VALUE As Double
     
On Error GoTo ERROR_LABEL
       
ATEMP_VAL = (RATE - DIVD - ((SIGMA * SIGMA) / 2)) * EXPIRATION
BTEMP_VAL = SIGMA * Sqr(EXPIRATION)
     
TEMP_SUM = 0
     
'   Generating the Quasi-Random vector and the option calculus to
'   feed the simulation
            
For i = 1 To nLOOPS
    PROB_VALUE = CORPUT_SEQUENCE_NUMBER_FUNC(i, 2)
    INV_VALUE = NORMSINV_FUNC(PROB_VALUE, 0, 1, 1)
    SIM_VALUE = SPOT * Exp(ATEMP_VAL + (BTEMP_VAL * INV_VALUE))
    If SIM_VALUE >= STRIKE Then: TEMP_SUM = TEMP_SUM + (SIM_VALUE - STRIKE)
Next i
            
EUROPEAN_CALL_CORPUT_SEQUENCE_SIMULATION_FUNC = (TEMP_SUM / nLOOPS) * Exp(-RATE * EXPIRATION)

Exit Function
ERROR_LABEL:
EUROPEAN_CALL_CORPUT_SEQUENCE_SIMULATION_FUNC = Err.number
End Function

'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : BOX_MULLER_CALL_OPTION_SIMULATION_FUNC

'DESCRIPTION   : Uses Box-Muller transformation control variate uses
'(STRIKE-CONTROL_VAL) antithetic variates
'Option Pricing via Monte Carlo Method using control and antithetic
'variates Vanilla Call

'LIBRARY       : DERIVATIVES
'GROUP         : BS_MONTE_CARLO
'ID            : 005
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 01/29/2009
'************************************************************************************
'************************************************************************************

Function BOX_MULLER_CALL_OPTION_SIMULATION_FUNC(ByVal SPOT As Double, _
ByVal STRIKE As Double, _
ByVal RATE As Double, _
ByVal SIGMA As Double, _
ByVal EXPIRATION As Double, _
Optional ByVal CONTROL_VAL As Double = 5, _
Optional ByVal MIN_LOOPS As Double = 100, _
Optional ByVal MAX_LOOPS As Double = 200, _
Optional ByVal DELTA_LOOPS As Double = 5, _
Optional ByVal RANDOMIZE_INT As Integer = 2, _
Optional ByVal CND_TYPE As Integer = 0)

Dim i As Long
Dim NSIZE As Long

Dim j As Double
Dim NROWS As Double

Dim A_VAL As Double
Dim B_VAL As Double
Dim C_VAL As Double
Dim D_VAL As Double
Dim E_VAL As Double
Dim F_VAL As Double
Dim G_VAL As Double
Dim H_VAL As Double
Dim I_VAL As Double
Dim J_VAL As Double
Dim K_VAL As Double

Dim PI_VAL As Double
Dim TEMP_SUM As Double

Dim TEMP_VAR As Double
Dim TEMP_VAL As Double

Dim FIRST_VAL As Double
Dim SECOND_VAL As Double

Dim OPTION_VAL As Double

Dim FIRST_RND_VAL As Double
Dim SECOND_RND_VAL As Double

Dim TEMP_MULT As Double
Dim TEMP_FACTOR As Double

Dim tolerance As Double

On Error GoTo ERROR_LABEL

PI_VAL = 3.14159265358979
NSIZE = CInt((MAX_LOOPS - MIN_LOOPS) / DELTA_LOOPS)

Rnd (-1)
Randomize (RANDOMIZE_INT)
tolerance = 0.000000000000001

'/////////////////////////Randomize [Number]\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
'The optional number argument is a Variant or any valid
'numeric expression.

'Randomize uses number to initialize the Rnd function's random-number
'generator, giving it a new seed value. If you omit number, the value
'returned by the system timer is used as the new seed value.

'If Randomize is not used, the Rnd function (with no arguments) uses
'the same number as a seed the first time it is called, and thereafter
'uses the last generated number as a seed value.

'Note: To repeat sequences of random numbers, call Rnd with a negative
'argument immediately before using Randomize with a numeric argument.
'Using Randomize with the same value for number does not repeat the
'previous sequence.
'///////////////////////////////////////////////////////////////////////

'Black-Scholes formula
A_VAL = (Log(SPOT / STRIKE) + (RATE + 0.5 * SIGMA * SIGMA) * _
    EXPIRATION) / (SIGMA * Sqr(EXPIRATION))
B_VAL = A_VAL - SIGMA * Sqr(EXPIRATION)
C_VAL = CND_FUNC(A_VAL, CND_TYPE)
D_VAL = CND_FUNC(B_VAL, CND_TYPE)

TEMP_VAL = SPOT * C_VAL - STRIKE * Exp(-RATE * EXPIRATION) * D_VAL

'Black-Scholes formula for the control variate
A_VAL = (Log(SPOT / (STRIKE - CONTROL_VAL)) + (RATE + 0.5 * SIGMA * SIGMA) * _
    EXPIRATION) / (SIGMA * Sqr(EXPIRATION))
B_VAL = A_VAL - SIGMA * Sqr(EXPIRATION)
C_VAL = CND_FUNC(A_VAL, CND_TYPE)
D_VAL = CND_FUNC(B_VAL, CND_TYPE)
K_VAL = SPOT * C_VAL - (STRIKE - CONTROL_VAL) * Exp(-RATE * EXPIRATION) * D_VAL

ReDim TEMP_MATRIX(0 To NSIZE, 1 To 4)

TEMP_MATRIX(0, 1) = "TRIALS"
TEMP_MATRIX(0, 2) = "ESTIMATE"
TEMP_MATRIX(0, 3) = "95% CONFIDENCE"
TEMP_MATRIX(0, 4) = "INTERVAL"

For i = 1 To NSIZE

    If i = 1 Then
      j = MIN_LOOPS
    Else: j = MIN_LOOPS + DELTA_LOOPS * (i - 1)
    End If
    If (j < 4) Or (j > 30000) Then: j = 4

    If (j Mod 2 = 1) Then: j = j + 1
    TEMP_MATRIX(i, 1) = j
    NROWS = j / 2

    TEMP_SUM = 0
    TEMP_VAR = 0
    C_VAL = SPOT * Exp((RATE - 0.5 * SIGMA * SIGMA) * EXPIRATION)
    D_VAL = Sqr(EXPIRATION)

    For j = 1 To NROWS
        FIRST_RND_VAL = Rnd
        SECOND_RND_VAL = Rnd
        'Box-Muller
        If FIRST_RND_VAL < tolerance Then FIRST_RND_VAL = tolerance _
        'choose smallest machine number > 0
        
        FIRST_VAL = Sqr(-2 * Log(FIRST_RND_VAL)) * Cos(2 * PI_VAL * SECOND_RND_VAL)
        SECOND_VAL = Sqr(-2 * Log(FIRST_RND_VAL)) * Sin(2 * PI_VAL * SECOND_RND_VAL)
        'simulate price in EXPIRATION
        A_VAL = C_VAL * Exp(SIGMA * D_VAL * FIRST_VAL)
        B_VAL = C_VAL * Exp(SIGMA * D_VAL * SECOND_VAL)
        'Caculate payoff
        If A_VAL > (STRIKE - CONTROL_VAL) Then _
            E_VAL = A_VAL - (STRIKE - CONTROL_VAL) Else E_VAL = 0
        If B_VAL > (STRIKE - CONTROL_VAL) Then _
            F_VAL = B_VAL - (STRIKE - CONTROL_VAL) Else F_VAL = 0
        If A_VAL > STRIKE Then A_VAL = A_VAL - STRIKE Else A_VAL = 0
        If B_VAL > STRIKE Then B_VAL = B_VAL - STRIKE Else B_VAL = 0
     
        TEMP_SUM = TEMP_SUM + (A_VAL - E_VAL) + (B_VAL - F_VAL)
        TEMP_VAR = TEMP_VAR + (A_VAL - E_VAL) * (A_VAL - E_VAL) + _
                (B_VAL - F_VAL) * (B_VAL - F_VAL)
        
        'Simulate price in EXPIRATION with antithetic variate
        G_VAL = C_VAL * Exp(SIGMA * D_VAL * -FIRST_VAL)
        H_VAL = C_VAL * Exp(SIGMA * D_VAL * -SECOND_VAL)
        
        'Calculate payoff
        If G_VAL > (STRIKE - CONTROL_VAL) Then _
            I_VAL = G_VAL - (STRIKE - CONTROL_VAL) Else I_VAL = 0
        If H_VAL > (STRIKE - CONTROL_VAL) Then _
            J_VAL = H_VAL - (STRIKE - CONTROL_VAL) Else J_VAL = 0
        If G_VAL > STRIKE Then G_VAL = G_VAL - STRIKE Else G_VAL = 0
        If H_VAL > STRIKE Then H_VAL = H_VAL - STRIKE Else H_VAL = 0
              
              
        TEMP_SUM = TEMP_SUM + (G_VAL - I_VAL) + (H_VAL - J_VAL)
        TEMP_VAR = TEMP_VAR + (G_VAL - I_VAL) * (G_VAL - I_VAL) + _
            (H_VAL - J_VAL) * (H_VAL - J_VAL) _
            + 2 * (G_VAL - I_VAL) * (A_VAL - E_VAL) + 2 * (H_VAL - J_VAL) * _
            (B_VAL - F_VAL)
            
    Next j

    TEMP_MULT = TEMP_SUM / (2 * 2 * NROWS)
    OPTION_VAL = K_VAL + Exp(-RATE * EXPIRATION) * TEMP_MULT

    TEMP_FACTOR = (TEMP_VAR / (4 * 2 * NROWS) - TEMP_MULT * TEMP_MULT) _
    / (2 * NROWS - 1)

    TEMP_FACTOR = Sqr(Exp(-RATE * EXPIRATION * 2) * TEMP_FACTOR)

    TEMP_MATRIX(i, 2) = OPTION_VAL
    TEMP_MATRIX(i, 3) = OPTION_VAL - 1.96 * TEMP_FACTOR
    TEMP_MATRIX(i, 4) = OPTION_VAL + 1.96 * TEMP_FACTOR
Next i

BOX_MULLER_CALL_OPTION_SIMULATION_FUNC = TEMP_MATRIX

Exit Function
ERROR_LABEL:
BOX_MULLER_CALL_OPTION_SIMULATION_FUNC = Err.number
End Function