## Sunday, January 26, 2014

### 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
'
' 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
'
' 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

1. Hi Mikael,

This is just what I'm looking for so thanks. I'm going to extend some features as I need, but I just made one change that you might find useful.

As I intend to move my final solution around across various file systems etc the "C:\temp' folder is a bit of a pain. You can make the program just look in the folder you have the workbook. I have shown the modified code below - you remove the path in the Declare statement and then you need to change the directory just once before any calls are made.

Private Declare Function nextMT Lib "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

ChDir ActiveWorkbook.Path

For i = 1 To nRows
For j = 1 To nCols
r(i, j) = InverseCumulativeNormal(nextMT())
Next j
Next i
'
IRandom_getNormalRandomMatrix = r
End Function

2. Glad to hear, that you have found this to be useful for your needs and thanks for the tip. Usually, for any "VBA production code", I am handling all external parameters (those nasty changing constants) like described in the following blog post:

http://mikejuniperhill.blogspot.fi/2013/05/configuring-your-constants-from-global.html

3. Hi Michael,
Would you mind sharing the example Excel spreadsheet with me (gkchang@gmail.com)?
Thanks so much,
George