Gaussian Copula implementation revisited

In my previous posting on matrix data structure, I compared the efficiency of processing matrix structures with different array data structures. The message was clear: two-dimensional arrays are the most efficient data structure when operating with large matrix operations. I have modified my implementation for Gaussian Copula in a such way, that it is using two-dimensional arrays instead of Matrix class.

Second version

The modified program is presented below. For required input data (Excel ranges), library references and external dll file, follow the instructions given on Gaussian Copula posting.

' standard VBA module (name = Enumerators)
Option Explicit
'
Public Enum E_
    P_SIMULATIONS = 1
    P_GENERATOR_TYPE = 2
    P_CORRELATION_MATRIX = 3
    P_TRANSFORM_TO_UNIFORM = 4
End Enum
'
'
'
' standard VBA module (name = MainProgram)
Option Explicit
'
Public Sub tester()
    '
    ' create correlation matrix object
    Dim correlation() As Double
    correlation = MatrixOperations.matrixFromRange(Sheets("Sheet1").Range("_correlation"))
    '
    ' create parameters for copula
    Dim parameters As New Scripting.Dictionary
    parameters.Add P_SIMULATIONS, CLng(Sheets("Sheet1").Range("_simulations").value)
    parameters.Add P_GENERATOR_TYPE, New MersenneTwister
    parameters.Add P_CORRELATION_MATRIX, correlation
    parameters.Add P_TRANSFORM_TO_UNIFORM, CBool(Sheets("Sheet1").Range("_transform").value)
    '
    ' create copula implementation
    Dim copula As ICopula: Set copula = New GaussianCopula
    copula.init parameters
    '
    ' get results from copula and write these into Excel
    Dim result() As Double: result = copula.getMatrix
    MatrixOperations.matrixToRange result, Sheets("Sheet1").Range("_dataDump")
    '
    ' release objects
    Set copula = Nothing
    Set parameters = Nothing
End Sub
'
'
'
' 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
'
' interface for copula model
'
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
'

Afterthoughts

In this version, Matrix class was replaced with plain two-dimensional array. Separate standard VBA module has been set to perform matrix operations (multiplication, transpose, cholesky decomposition, data input/output). Everything else is unchanged. With this small modification, we have reached reduction in processing times, when simulating correlated random numbers with Gaussian Copula model.

-Mike