# Monte Carlo Simulation

`Function European_Call_bMCBS( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    N As Long, _    M As Long) As Double​'European Call Option by Monte Carlo with Black-Scholes'Clewlow/Strickland: Implementing Derivatives Models, page 85​Dim dt As DoubleDim nudt As Double, sigsdt As DoubleDim lnS As Double, lnSt As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St As Double, CT As DoubleDim SD As Double, SE As DoubleDim i As Long, j As Long​'Precompute constants​dt = T / Nnudt = (r - div - 0.5 * sig ^ 2#) * dtsigsdt = sig * Sqr(dt)lnS = Log(S)​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        lnSt = lnS        For i = 1 To N            'For each time step        '        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample        lnSt = lnSt + nudt + sigsdt * e 'Evolve the stock price            Next i        St = Exp(lnSt)    CT = St - K: If CT < 0# Then CT = 0#    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Call_bMCBS = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Call_bMCBS_AVR( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    N As Long, _    M As Long) As Double​'European Call Option by Monte Carlo with Black-Scholes and Antithetic Variance Reduction'Clewlow/Strickland: Implementing Derivatives Models, page 89​Dim dt As DoubleDim nudt As Double, sigsdt As DoubleDim lnS As Double, LnSt1 As Double, LnSt2 As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St1 As Double, St2 As DoubleDim CT As Double, CT1 As Double, CT2 As DoubleDim SD As Double, SE As DoubleDim i As Long, j As Long​'Precompute constants​dt = T / Nnudt = (r - div - 0.5 * sig ^ 2#) * dtsigsdt = sig * Sqr(dt)lnS = Log(S)​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        LnSt1 = lnS    LnSt2 = lnS        For i = 1 To N            'For each time step        '        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample        LnSt1 = LnSt1 + nudt + sigsdt * e 'Evolve the stock price        LnSt2 = LnSt2 + nudt + sigsdt * -e 'Evolve the stock price            Next i        St1 = Exp(LnSt1)    St2 = Exp(LnSt2)    CT1 = St1 - K: If CT1 < 0# Then CT1 = 0#    CT2 = St2 - K: If CT2 < 0# Then CT2 = 0#    CT = (CT1 + CT2) / 2#    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Call_bMCBS_AVR = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Call_bMCBS_DCV( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    N As Long, _    M As Long) As Double​'European Call Option by Monte Carlo with Black-Scholes and Delta-based Control Variate'Clewlow/Strickland: Implementing Derivatives Models, page 97​Dim dt As DoubleDim nudt As Double, sigsdt As Double, erddt As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St As Double, Stn As DoubleDim CT As DoubleDim SD As Double, SE As DoubleDim beta1 As Double, t1 As DoubleDim cv As DoubleDim delta1 As DoubleDim i As Long, j As Long​'Dim eps As Variant'eps = Array(0.5391, 0.1692, 0.4583, 0.4164, 0.9223, 0.1597, 0.4011, -0.2382, 1.2152, -0.5506)​'Precompute constants​dt = T / Nnudt = (r - div - 0.5 * sig ^ 2#) * dtsigsdt = sig * Sqr(dt)erddt = Exp((r - div) * dt)beta1 = -1#'Debug.Print "dt = " & dt & ", nudt = " & nudt & ", sigsdt = " & sigsdt & ", erddt = " & erddt & ", beta1 = " & beta1​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        St = S    cv = 0#        For i = 1 To N            'For each time step                t1 = CDbl(i - 1) * dt        delta1 = Black_Scholes_Delta(St, t1, K, T, sig, r, div)'        Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv'        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample'        e = eps(i - 1)        Stn = St * Exp(nudt + sigsdt * e)        cv = cv + delta1 * (Stn - St * erddt) '* exp(T - (t1 + dt)) 'inflation factor *exp can be omitted        St = Stn            Next i    '    Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv    CT = St - K: If CT < 0# Then CT = 0#    CT = CT + beta1 * cv    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT'    Debug.Print "CT = " & CT & ", CT2 = " & CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Call_bMCBS_DCV = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Call_bMCBS_ADCV( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    N As Long, _    M As Long) As Double​'European Call Option by Monte Carlo with Black-Scholes and Antithetic and Delta-based Control Variates'Clewlow/Strickland: Implementing Derivatives Models, page 100​Dim dt As DoubleDim nudt As Double, sigsdt As Double, erddt As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St1 As Double, Stn1 As DoubleDim St2 As Double, Stn2 As DoubleDim CT As Double, CT1 As Double, CT2 As DoubleDim SD As Double, SE As DoubleDim beta1 As Double, t1 As DoubleDim cv1 As Double, cv2 As DoubleDim delta1 As Double, delta2 As DoubleDim i As Long, j As Long​'Dim eps As Variant'eps = Array(0.5391, 0.1692, 0.4583, 0.4164, 0.9223, 0.1597, 0.4011, -0.2382, 1.2152, -0.5506)​'Precompute constants​dt = T / Nnudt = (r - div - 0.5 * sig ^ 2#) * dtsigsdt = sig * Sqr(dt)erddt = Exp((r - div) * dt)beta1 = -1#'Debug.Print "dt = " & dt & ", nudt = " & nudt & ", sigsdt = " & sigsdt & ", erddt = " & erddt & ", beta1 = " & beta1​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        St1 = S    St2 = S    cv1 = 0#    cv2 = 0#        For i = 1 To N            'For each time step                t1 = CDbl(i - 1) * dt        delta1 = Black_Scholes_Delta(St1, t1, K, T, sig, r, div)        delta2 = Black_Scholes_Delta(St2, t1, K, T, sig, r, div)'        Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv'        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample'        e = eps(i - 1)        Stn1 = St1 * Exp(nudt + sigsdt * e)        Stn2 = St2 * Exp(nudt + sigsdt * -e)        cv1 = cv1 + delta1 * (Stn1 - St1 * erddt) '* exp(T - (t1 + dt)) 'inflation factor *exp can be omitted        cv2 = cv2 + delta2 * (Stn2 - St2 * erddt) '* exp(T - (t1 + dt)) 'inflation factor *exp can be omitted        St1 = Stn1        St2 = Stn2            Next i    '    Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv    CT1 = St1 - K: If CT1 < 0# Then CT1 = 0#    CT2 = St2 - K: If CT2 < 0# Then CT2 = 0#    CT = (CT1 + beta1 * cv1 + CT2 + beta1 * cv2) / 2#    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT'    Debug.Print "CT = " & CT & ", CT2 = " & CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Call_bMCBS_ADCV = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Call_bMCBS_ADGCV( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    N As Long, _    M As Long) As Double​'European Call Option by Monte Carlo with Black-Scholes and Antithetic, Delta-'and Gamma-based Control Variates'Clewlow/Strickland: Implementing Derivatives Models, page 102​Dim dt As DoubleDim nudt As Double, sigsdt As DoubleDim erddt As Double, egamma As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St1 As Double, Stn1 As DoubleDim St2 As Double, Stn2 As DoubleDim CT As Double, CT1 As Double, CT2 As DoubleDim SD As Double, SE As DoubleDim beta1 As Double, beta2 As Double, t1 As DoubleDim cv1 As Double, cv2 As DoubleDim delta1 As Double, delta2 As DoubleDim gamma1 As Double, gamma2 As DoubleDim i As Long, j As Long​'Dim eps As Variant'eps = Array(0.5391, 0.1692, 0.4583, 0.4164, 0.9223, 0.1597, 0.4011, -0.2382, 1.2152, -0.5506)​'Precompute constants​dt = T / Nnudt = (r - div - 0.5 * sig ^ 2#) * dtsigsdt = sig * Sqr(dt)erddt = Exp((r - div) * dt)egamma = Exp((2# * (r - div) + sig * sig) * dt) - 2# * erddt + 1#beta1 = -1#beta2 = -0.5'Debug.Print "dt = " & dt & ", nudt = " & nudt & ", sigsdt = " & sigsdt & ", erddt = " & erddt & ", beta1 = " & beta1​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        St1 = S    St2 = S    cv1 = 0#    cv2 = 0#        For i = 1 To N            'For each time step                'Compute hedge sensitivities        t1 = CDbl(i - 1) * dt        delta1 = Black_Scholes_Delta(St1, t1, K, T, sig, r, div)        delta2 = Black_Scholes_Delta(St2, t1, K, T, sig, r, div)        gamma1 = Black_Scholes_Gamma(St1, t1, K, T, sig, r, div)        gamma2 = Black_Scholes_Gamma(St2, t1, K, T, sig, r, div)                'Evolve asset prices'        Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv'        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample'        e = eps(i - 1)        Stn1 = St1 * Exp(nudt + sigsdt * e)        Stn2 = St2 * Exp(nudt + sigsdt * -e)                'Accumulate control variates        cv1 = cv1 + delta1 * (Stn1 - St1 * erddt) + _                    delta2 * (Stn2 - St2 * erddt)        cv2 = cv2 + gamma1 * ((Stn1 - St1) ^ 2 - St1 ^ 2 * egamma) + _                    gamma2 * ((Stn2 - St2) ^ 2 - St2 ^ 2 * egamma)        St1 = Stn1        St2 = Stn2            Next i    '    Debug.Print "e = " & e & ", St = " & St & ", delta1 = " & delta1 & ", cv = " & cv    CT1 = St1 - K: If CT1 < 0# Then CT1 = 0#    CT2 = St2 - K: If CT2 < 0# Then CT2 = 0#    CT = (CT1 + CT2 + beta1 * cv1 + beta2 * cv2) / 2#    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT'    Debug.Print "CT = " & CT & ", CT2 = " & CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Call_bMCBS_ADGCV = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Spread_bMCBS( _    K As Double, _    T As Double, _    S1 As Double, _    S2 As Double, _    sig1 As Double, _    sig2 As Double, _    div1 As Double, _    div2 As Double, _    rho As Double, _    r As Double, _    N As Long, _    M As Long) As Double​'European Spread Option by Monte Carlo with Black-Scholes'Clewlow/Strickland: Implementing Derivatives Models, page 110​Dim dt As DoubleDim nu1dt As Double, nu2dt As DoubleDim sig1sdt As Double, sig2sdt As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e1 As Double, e2 As DoubleDim z1 As Double, z2 As DoubleDim St1 As Double, St2 As DoubleDim CT As DoubleDim SD As Double, SE As DoubleDim srho As DoubleDim i As Long, j As Long​'Precompute constants​dt = T / Nnu1dt = (r - div1 - 0.5 * sig1 ^ 2#) * dtnu2dt = (r - div2 - 0.5 * sig2 ^ 2#) * dtsig1sdt = sig1 * Sqr(dt)sig2sdt = sig2 * Sqr(dt)srho = Sqr(1# - rho ^ 2#)​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        St1 = S1    St2 = S2        For i = 1 To N            'For each time step                e1 = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e2 = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample'        e1 = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _'            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample'        e2 = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _'            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample        z1 = e1        z2 = rho * e1 + srho * e2        St1 = St1 * Exp(nu1dt + sig1sdt * z1)        St2 = St2 * Exp(nu2dt + sig2sdt * z2)            Next i        CT = St1 - St2 - K: If CT < 0# Then CT = 0#    sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT    Next j​SD = Sqr(sum_CT * sum_CT / M - sum_CT2) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Spread_bMCBS = sum_CT / M * Exp(-r * T)​End Function`
`Function European_Down_and_Out_Call_bMCBS( _    K As Double, _    T As Double, _    S As Double, _    sig As Double, _    r As Double, _    div As Double, _    H As Long, _    N As Long, _    M As Long) As Double​'European Down and Out Call Option by Monte Carlo with Black-Scholes'Clewlow/Strickland: Implementing Derivatives Models, page 116​Dim dt As DoubleDim nudt As Double, sigsdt As DoubleDim lnS As Double, lnSt As DoubleDim sum_CT As Double, sum_CT2 As DoubleDim e As DoubleDim St As Double, CT As DoubleDim SD As Double, SE As DoubleDim i As Long, j As LongDim bBarrierCrossed As Boolean​#If MODULE_TEST Then    e = Rnd(-1)    Randomize 1 'maybe even additional parameter to function#End If​'Precompute constants​dt = T / Nnudt = (r - div - sig * sig / 2#) * dtsigsdt = sig * Sqr(dt)​sum_CT = 0#sum_CT2 = 0#​For j = 1 To M​    'For each simulation        St = S    bBarrierCrossed = False        For i = 1 To N            'For each time step        '        e = Application.WorksheetFunction.NormSInv(Rnd()) 'Standard normal sample        e = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + _            Rnd() + Rnd() + Rnd() + Rnd() - 6# 'Standard normal sample        St = St * Exp(nudt + sigsdt * e) 'Evolve the stock price        If St <= H Then            bBarrierCrossed = True            Exit For        End If            Next i        If bBarrierCrossed Then        CT = 0#    Else        CT = St - K        If CT < 0# Then CT = 0#    End If        sum_CT = sum_CT + CT    sum_CT2 = sum_CT + CT * CT    Next j​SD = Sqr(Abs(sum_CT * sum_CT / M - sum_CT2)) * Exp(-2# * r * T) / (M - 1) 'Standard deviationSE = SD / Sqr(M) 'Standard errorEuropean_Down_and_Out_Call_bMCBS = sum_CT / M * Exp(-r * T)​End Function`