Pricing equity basket option with Monte Carlo in VBA

Within the last three postings, I have been chewing the creation of correlated random numbers with Gaussian Copula model and some efficiency issues, when handling large matrices in VBA. This post is finally presenting some real-world application for these issues by pricing equity basket option with Monte Carlo method in VBA.

EQUITY BASKET OPTION
A basket option is an option where the payout is related to the cumulative performance of a specified basket of underlying assets. Typically examples include options on baskets of currencies and equities. In this option, the multifactor element is the multiple underlying assets incorporated into option structure. The payoff of single option is based on the value at expiry or exercise value of individual asset. In contrast, the payoff of basket option is based on the aggregate value of basket of assets.

There are analytical pricing formulas available for basket options, but in this example we will use Monte Carlo method for pricing. Pricing procedure for M assets and N simulations can be described with the following steps

1. Create matrix (NxM) of correlated normal random numbers by using Gaussian Copula
2. Simulate M asset paths by using standard Geometric Brownian Motion (GBM)
3. Sum all M simulated asset prices at expiration to get aggregate asset value at expiration
4. Calculate option payoff for aggregate asset value at expiration
5. Discount calculated payoff back to today
6. Repeat steps 2-5 N times
7. Calculate average for all discounted payoffs to get the option value today

For the structure presented in this example, the strike price is calculated as the sum of observed spot prices today.

PROGRAM DESIGN
We are aiming to have as flexible design as possible. This means, that first we try to identify the main components what might change in the future. The current design is presented in the picture below.





















We can identify at least four important components, which we might want to change in the future.

1. Payoff - even we are pricing multifactor option, the payoff itself is one-factor payoff since we are calculating the difference of aggregate asset value at expiration and the strike price. In the future, we might want to create implementations for other types of payoffs. As the other parts of the program are communicating with payoff interface, it is easy to change the payoff type without having any modifications for other parts of the program.
2. Copula - in the future, we might want to create implementations for other types of copula models.
3. Random - in the future, we might want to create implementations for other types of random number generators.
4. Option - this example is pricing equity basket option, but it is also an implementation of multi-factor option interface. In the future, we might want to create implementations for other types of multi-factor options, such as the best-of or worst-of basket options. 

With this flexibility in our hands, we are already able to create pricers for many different types of multi-factor options with this proposed design. There are also some other small functionalities, which we might also have isolated behind separate interfaces, like Stochastic differential equation used (GBM) and handling interest rates and discount factors.

EXCEL INTERFACE
Input parameters for the program are read from Excel worksheet. The screenshot of interface is shown in the picture below.

Important input parameters are marked with yellow color. The program is configured to read all required input parameters from the following named Excel ranges.

1. Correlation matrix - range("F5:N13") = "_correlation"
2. Spot prices - range("P5:P13") = "_spot"
3. Spot voltilities - range("Q5:Q13") = "_vol"
4. Discount curve - range("T5:U19") = "_curve"
5. Number of simulations - range("C18") = "_simulations"
6. Maturity - range("C20") = "_maturity"
7. Error reporting - range("C21") = "_status"
8. Results reporting - range("C25") = "_result"

Random number generator type (Mersenne Twister) has been hard-coded inside the main program. The pink area in the picture is the correlation matrix override region. Correlations in the yellow area are coming directly from bloomberg (like the other market data), but can be overriden with values from override matrix, if needed. This is useful feature to have, if we want to investigate the effect of correlation between the assets on option values.

Note, that the number of assets in the basket is not hard-coded in anywhere. There is not any particular limitation for having exactly those nine stocks in the basket. When the parameters are created for the actual program, Copula implementation is creating its random number matrix based on the dimensions of given correlation matrix and the number of simulations. If we define the named Excel range for correlation to be (3 x 3) matrix and the number of simulations to be 250 000, Copula creates (250 000 x 3) matrix of correlated random numbers. This means, that the actual simulation procedure is then having total of 250 000 paths for each asset. The only requirement is, that the number of rows and columns of correlation matrix must be correspond to the number of rows in spot and volatility ranges.

THE PROGRAM
Copy-paste the following program and break it up into different VBA modules and classes. Follow the instructions given. For external library needed (dll), check my post on Gaussian Copula. Remember also to create reference into Microsoft Scripting Runtime library. Also, insert one ActiveX command button into pricer worksheet, name it as "btnExecute" and create event as shown below. The program starts, as we press this button.

' VBA worksheet (name = pricer)
Option Explicit
'
Private Sub btnExecute_Click()
    '
    On Error GoTo errorhandler: Range("_status") = ""
    '
    ' create parameter wrapper for copula
    Dim p_copula As New Scripting.Dictionary
    '
    ' add parameters for copula wrapper
    Dim correlation() As Double: correlation = MatrixOperations.matrixFromRange(Range("_correlation"))
    p_copula.Add P_CORRELATION_MATRIX, correlation
    p_copula.Add P_SIMULATIONS, CLng(Range("_simulations").value)
    p_copula.Add P_GENERATOR_TYPE, New MersenneTwister ' HARD-CODED !!
    p_copula.Add P_TRANSFORM_TO_UNIFORM, CBool("FALSE") ' HARD-CODED !!
    '
    ' create copula implementation and feed wrapper into it
    Dim copula As ICopula: Set copula = New GaussianCopula
    copula.init p_copula
    '
    '
    ' create parameter wrapper for multifactor option
    Dim p_option As New Scripting.Dictionary
    '
    ' add parameters for multifactor option wrapper
    p_option.Add P_COPULA_MODEL, copula
    p_option.Add P_SIMULATIONS, CLng(Range("_simulations").value)
    p_option.Add P_MATURITY, CDbl(Range("_maturity"))
    '
    Dim curve() As Double: curve = MatrixOperations.matrixFromRange(Range("_curve"))
    p_option.Add P_ZERO_CURVE, curve
    '
    Dim spot() As Double: spot = MatrixOperations.matrixFromRange(Range("_spot"))
    p_option.Add P_SPOT, spot
    '
    Dim vol() As Double: vol = MatrixOperations.matrixFromRange(Range("_vol"))
    p_option.Add P_VOLATILITY, vol
    '
    ' calculate basket strike price (assuming equal weights)
    Dim i As Long, strike As Double
    For i = 1 To UBound(spot)
        strike = strike + spot(i, 1)
    Next i
    '
    ' create payoff implementation and feed strike into it
    Dim payoff As IOneFactorPayoff: Set payoff = New OneFactorCallPayoff
    payoff.init strike
    p_option.Add P_PAYOFF, payoff
    '
    '
    ' create equity basket option and feed wrapper into it
    Dim basketOption As IMultiFactorOption: Set basketOption = New EquityBasketOption
    basketOption.init p_option
    '
    ' receive and write calculation results into Excel worksheet
    Dim result() As Double: result = basketOption.getStatistics()
    MatrixOperations.matrixToRange result, Range("_result")
    Exit Sub
    '
errorhandler:
    Range("_status") = Err.Description
End Sub
'
'
'
'
' standard VBA module (name = Enumerators)
Option Explicit
'
Public Enum E_
    P_CORRELATION_MATRIX = 1
    P_SPOT = 2
    P_VOLATILITY = 3
    P_ZERO_CURVE = 4
    P_SIMULATIONS = 5
    P_GENERATOR_TYPE = 6
    P_MATURITY = 7
    P_COPULA_MODEL = 8
    P_TRANSFORM_TO_UNIFORM = 9
    P_PAYOFF = 10
    P_STRIKE = 11
End Enum
'
'
'
'
' standard VBA module (name = MatrixOperations)
Option Explicit
'
Public Function multiplication(ByRef m1() As Double, ByRef m2() As Double) As Double()
    '
    ' get matrix multiplication
    Dim result() As Double: ReDim result(1 To UBound(m1, 1), 1 To UBound(m2, 2))
    Dim i As Long, j As Long, k As Long
    '
    Dim r1 As Long: r1 = UBound(m1, 1)
    Dim c1 As Long: c1 = UBound(m1, 2)
    Dim r2 As Long: r2 = UBound(m2, 1)
    Dim c2 As Long: c2 = UBound(m2, 2)
    Dim v As Double
    '
    For i = 1 To r1
        For j = 1 To c2
            v = 0
            '
            For k = 1 To c1
                v = v + m1(i, k) * m2(k, j)
            Next k
            result(i, j) = v
        Next j
    Next i
    multiplication = result
End Function
'
Public Function transpose(ByRef m() As Double)
    '
    ' get matrix transpose
    Dim nRows As Long: nRows = UBound(m, 1)
    Dim nCols As Long: nCols = UBound(m, 2)
    Dim result() As Double: ReDim result(1 To nCols, 1 To nRows)
    '
    Dim i As Long, j As Long
    For i = 1 To nRows
        For j = 1 To nCols
            result(j, i) = m(i, j)
        Next j
    Next i
    transpose = result
End Function
'
Public Function cholesky(ByRef c() As Double) As Double()
    '
    ' create cholesky decomposition, a lower triangular matrix
    ' d = decomposition, c = correlation matrix
    Dim s As Double
    Dim n As Long: n = UBound(c, 1)
    Dim m As Long: m = UBound(c, 2)
    Dim d() As Double: ReDim d(1 To n, 1 To m)
    '
    Dim i As Long, j As Long, k As Long
    For j = 1 To n
        s = 0
        For k = 1 To j - 1
            s = s + d(j, k) ^ 2
        Next k
        d(j, j) = (c(j, j) - s)
        If (d(j, j) <= 0) Then Exit For
        '
        d(j, j) = Sqr(d(j, j))
        '
        For i = (j + 1) To n
            s = 0
            For k = 1 To (j - 1)
                s = s + d(i, k) * d(j, k)
            Next k
            d(i, j) = (c(i, j) - s) / d(j, j)
        Next i
    Next j
    cholesky = d
End Function
'
Public Function matrixToRange(ByRef m() As Double, ByRef r As Range)
    '
    ' write matrix into Excel range
    r.Resize(UBound(m, 1), UBound(m, 2)) = m
End Function
'
Public Function matrixFromRange(ByRef r As Range) As Double()
    '
    ' read matrix from Excel range
    Dim v As Variant: v = r.Value2
    Dim r_ As Long: r_ = UBound(v, 1)
    Dim c_ As Long: c_ = UBound(v, 2)
    Dim m() As Double: ReDim m(1 To r_, 1 To c_)
    '
    ' transform variant to double
    Dim i As Long, j As Long
    For i = 1 To r_
        For j = 1 To c_
            m(i, j) = CDbl(v(i, j))
        Next j
    Next i
    matrixFromRange = m
End Function
'
'
'
'
' VBA class module (name = ICopula)
Option Explicit
'
Public Function init(ByRef parameters As Scripting.Dictionary)
    ' interface - no implementation
End Function
'
Public Function getMatrix() As Double()
    ' interface - no implementation
End Function
'
'
'
'
' VBA class module (name = GaussianCopula)
Option Explicit
'
Implements ICopula
'
Private n As Long ' number of simulations
Private transform As Boolean ' condition for uniform transformation
Private generator As IRandom ' random number generator implementation
'
Private c() As Double ' correlation matrix
Private d() As Double ' cholesky decomposition matrix
Private z() As Double ' independent normal random variables
Private x() As Double ' correlated normal random variables
'
Private Function ICopula_init(ByRef parameters As Scripting.Dictionary)
    '
    ' initialize class data and objects
    n = parameters(P_SIMULATIONS)
    transform = parameters(P_TRANSFORM_TO_UNIFORM)
    Set generator = parameters(P_GENERATOR_TYPE)
    c = parameters(P_CORRELATION_MATRIX)
End Function
'
Private Function ICopula_getMatrix() As Double()
    '
    ' create matrix of independent normal random numbers
    ReDim z(1 To n, 1 To UBound(c, 2))
    z = generator.getNormalRandomMatrix(UBound(z, 1), UBound(z, 2))
    '
    ' create cholesky decomposition
    d = MatrixOperations.cholesky(c)
    '
    ' create correlated normal random numbers
    z = MatrixOperations.transpose(z)
    x = MatrixOperations.multiplication(d, z)
    x = MatrixOperations.transpose(x)
    '
    ' transform correlated normal random numbers
    ' into correlated uniform random numbers
    If (transform) Then uniformTransformation
    ICopula_getMatrix = x
End Function
'
Private Function uniformTransformation()
    '
    ' map normal random number to uniform plane
    Dim nRows As Long: nRows = UBound(x, 1)
    Dim nCols As Long: nCols = UBound(x, 2)
    '
    Dim i As Long, j As Long
    For i = 1 To nRows
        For j = 1 To nCols
            x(i, j) = WorksheetFunction.NormSDist(x(i, j))
        Next j
    Next i
End Function
'
'
'
'
' VBA class module (name = IRandom)
Option Explicit
'
Public Function getNormalRandomMatrix( _
    ByVal nRows As Long, _
    ByVal nCols As Long) As Double()
    '
    ' interface - no implementation
    ' takes in two parameters (number of rows and columns) and
    ' returns double() filled with normal random variates
End Function
'
'
'
'
' VBA class module (name = MersenneTwister)
Option Explicit
'
Implements IRandom
'
Private Declare Function nextMT Lib "C:\temp\mt19937.dll" Alias "genrand" () As Double
'
Private Function IRandom_getNormalRandomMatrix( _
    ByVal nRows As Long, _
    ByVal nCols As Long) As Double()
    '
    ' retrieve NxM matrix with normal random numbers
    Dim r() As Double: ReDim r(1 To nRows, 1 To nCols)
    Dim i As Long, j As Long
    For i = 1 To nRows
        For j = 1 To nCols
            r(i, j) = InverseCumulativeNormal(nextMT())
        Next j
    Next i
    '
    IRandom_getNormalRandomMatrix = r
End Function
'
Public Function InverseCumulativeNormal(ByVal p As Double) As Double
    '
    ' Define coefficients in rational approximations
    Const a1 = -39.6968302866538
    Const a2 = 220.946098424521
    Const a3 = -275.928510446969
    Const a4 = 138.357751867269
    Const a5 = -30.6647980661472
    Const a6 = 2.50662827745924
    '
    Const b1 = -54.4760987982241
    Const b2 = 161.585836858041
    Const b3 = -155.698979859887
    Const b4 = 66.8013118877197
    Const b5 = -13.2806815528857
    '
    Const c1 = -7.78489400243029E-03
    Const c2 = -0.322396458041136
    Const c3 = -2.40075827716184
    Const c4 = -2.54973253934373
    Const c5 = 4.37466414146497
    Const c6 = 2.93816398269878
    '
    Const d1 = 7.78469570904146E-03
    Const d2 = 0.32246712907004
    Const d3 = 2.445134137143
    Const d4 = 3.75440866190742
    '
    'Define break-points
    Const p_low = 0.02425
    Const p_high = 1 - p_low
    '
    'Define work variables
    Dim q As Double, r As Double
    '
    'If argument out of bounds, raise error
    If p <= 0 Or p >= 1 Then Err.Raise 5
    '
    If p < p_low Then
        '
        'Rational approximation for lower region
        q = Sqr(-2 * Log(p))
        InverseCumulativeNormal = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
        ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
        '
    ElseIf p <= p_high Then
        'Rational approximation for lower region
        q = p - 0.5
        r = q * q
        InverseCumulativeNormal = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / _
        (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1)
        '
    ElseIf p < 1 Then
        'Rational approximation for upper region
        q = Sqr(-2 * Log(1 - p))
        InverseCumulativeNormal = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
        ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    End If
End Function
'
'
'
'
' VBA class module (name = IOneFactorPayoff)
Option Explicit
'
Public Function init(ByVal x As Double)
    ' interface
End Function
'
Public Function getPayoff(ByVal s As Double) As Double
    ' interface
End Function
'
'
'
'
' VBA class module (name = OneFactorCallPayoff)
Option Explicit
'
Implements IOneFactorPayoff
'
Private x_ As Double
'
Private Function IOneFactorPayoff_getPayoff(ByVal s As Double) As Double
    IOneFactorPayoff_getPayoff = 0
    If (s > x_) Then IOneFactorPayoff_getPayoff = s - x_
End Function
'
Private Function IOneFactorPayoff_init(ByVal x As Double)
    x_ = x
End Function
'
'
'
'
' VBA class module (name = IMultiFactorOption)
Option Explicit
'
Public Function init(ByRef wrapper As Scripting.Dictionary)
    ' interface
End Function
'
Public Function getStatistics() As Double()
    ' interface
End Function
'
'
'
'
' VBA class module (name = EquityBasketOption)
Option Explicit
'
Private Declare Function GetTickCount Lib "kernel32.dll" () As Long
'
Implements IMultiFactorOption
'
Private payoff As IOneFactorPayoff
Private copula As ICopula
'
Private spot() As Double
Private vol() As Double
Private statistics() As Double
Private random() As Double
Private curve() As Double
'
Private maturity As Double
Private n As Long
Private startTime As Long
'
Private Function IMultiFactorOption_init(ByRef wrapper As Scripting.IDictionary)
    '
    ' extract objects and class member data from wrapper
    Set copula = wrapper(P_COPULA_MODEL)
    Set payoff = wrapper(P_PAYOFF)
    spot = wrapper(P_SPOT)
    vol = wrapper(P_VOLATILITY)
    curve = wrapper(P_ZERO_CURVE)
    n = wrapper(P_SIMULATIONS)
    maturity = wrapper(P_MATURITY)
End Function
'
Private Function IMultiFactorOption_getStatistics() As Double()
    '
    ' record start time
    startTime = GetTickCount()
    '
    ' get matrix of correlated normal random variates
    random = copula.getMatrix()
    '
    ' execute monte carlo and return results
    executeMonteCarloProcess
    IMultiFactorOption_getStatistics = statistics
End Function
'
Private Function executeMonteCarloProcess()
    '
    ' the actual monte carlo process
    Dim terminalValue() As Double: ReDim terminalValue(1 To n) ' container for all calculated payoffs
    Dim i As Long, j As Long, m As Long: m = UBound(spot)
    Dim r As Double: r = getRate(maturity) ' rate for risk-neutral drift
    Dim SDE As Double, drift As Double, noise As Double ' SDE process and its sub-components
    Dim value As Double
    '
    ' process random numbers
    For i = 1 To n
        value = 0
        '
        ' SDE process for m assets
        For j = 1 To m
            '
            ' since we have european option, we can create SDE process directly
            ' into maturity instead of simulating small time increments
            drift = (r - 0.5 * vol(j, 1) * vol(j, 1))
            noise = vol(j, 1) * Sqr(maturity) * random(i, j)
            SDE = spot(j, 1) * Exp(drift * maturity + noise)
            value = value + SDE
        Next j
        '
        ' push calculated basket payoffs into container
        terminalValue(i) = payoff.getPayoff(value)
    Next i
    '
    ' average and discount payoffs to get basket option value today
    Dim df As Double: df = getDF(maturity) ' zero-coupon bond price for discounting expectation
    Dim sum As Double, sumSq As Double
    For i = 1 To n
        value = terminalValue(i) * df
        sum = sum + value
        sumSq = sumSq + (value * value)
    Next i
    '
    ' HARD-CODED !!
    ' push statistics into result container (average, sample standard error, time elapsed)
    ReDim statistics(1 To 3, 1 To 1)
    statistics(1, 1) = (sum / n)
    statistics(2, 1) = Sqr((sumSq / n) - ((sum / n) * (sum / n))) / Sqr(n)
    statistics(3, 1) = ((GetTickCount - startTime) / 1000)
End Function
'
Private Function getDF(ByVal t As Double) As Double
    '
    ' calculate continuous-time discount factor
    Dim r As Double: r = linearInterpolation(t)
    getDF = Exp(-r * t)
End Function
'
Private Function getRate(ByVal t As Double) As Double
    '
    ' get zero-coupon bond yield
    getRate = linearInterpolation(t)
End Function
'
Private Function linearInterpolation(ByVal t As Double) As Double
    '
    ' get linear interpolation
    Dim n As Integer: n = UBound(curve, 1)
    Dim v As Double
    Dim i As Long
    '
    ' checking for boundaries
    If (t < curve(1, 1)) Then v = curve(1, 2): GoTo exitPoint
    If (t > curve(UBound(curve, 1), 1)) Then v = curve(UBound(curve, 1), 2): GoTo exitPoint
    '
    For i = 1 To (n - 1)
        If (t >= curve(i, 1)) And (t <= curve(i + 1, 1)) Then
            v = curve(i, 2) + (curve(i + 1, 2) - curve(i, 2)) * ((t - curve(i, 1)) / (curve(i + 1, 1) - curve(i, 1)))
            Exit For
        End If
    Next i
    '
exitPoint:
    linearInterpolation = v
End Function
'

PRICING RUN
In the testing run, equity basket has nine different GBP-nominated shares. Strike price has been calculated as the sum of current spot prices to be 12308.4 GBP. With the current correlation structure between the assets and 1.2 million simulated sample paths for each asset in the basket, the option value has been calculated to be approximately 631.8 GBP.

When changing all correlations to be exactly one with correlation override matrix, we get the value of 856.73 GBP, which should be approximately the same value as the sum of values of individually priced equity options. Plain Black & Scholes pricing model is giving the value of 856.32 GBP for such basket, so we are close enough.

Thanks for reading!

-Mike