“When your work speaks for itself, don’t interrupt.” [Henry J. Kaiser]
Abstract
If you need to generate correlated random numbers, the Iman Conover approach is a good method which is to be preferred to the Cholesky decomposition.
In 1982 Iman and Conover published their original article (external link!) “A distribution-free approach to inducing rank correlation among input variables”.
Rick Wicklin wrote in 2021 on (external link!) “Simulate correlated variables by using the Iman-Conover transformation”. His article includes a SAS implementation of the Iman Conover approach.
In 2005 Stephen J. Mildenhall published (external link!) “Correlation and Aggregate Loss Distributions with an Emphasis on the Iman-Conover Method”.
I implemented the example given in Mildenhall’s paper with Excel worksheet functions as well as with Excel / VBA:
The input matrix X:
The target correlation S:
The Cholesky decomposition C of S:
(Please compare to Cholesky)
The intermediate matrix M (constant values to equal Mildenhall’s data):
You can create similar data automatically with an array formula in A1:A20:
=NORMSINV(SEQUENCE(20,1,1,1)/21)/STDEVPA(NORMSINV(SEQUENCE(20,1,1,1)/21))
or with the array formula
=RandomShuffle($A$1:$A$20)
in cells B1:B20 (copy to columns C and D respectively).
Now you get the covariance matrix E:
And its Cholesky decomposition F:
The intermediate matrix T:
You can check the generated correlations:
Calculate the ranks of numbers in the columns of T:
Finally you get your result:
Appendix – RandomShuffle Code
Please read my Disclaimer.
Function RandomShuffle(vtemp As Variant) As Variant
Dim j As Long, k As Long, n As Long
Dim temp As Double, u As Double
Application.Volatile
With Application.WorksheetFunction
vtemp = .Transpose(.Transpose(vtemp))
n = UBound(vtemp, 1)
j = n
Do While j > 0
u = Rnd()
k = Int(j * u + 1)
temp = vtemp(j, 1)
vtemp(j, 1) = vtemp(k, 1)
vtemp(k, 1) = temp
j = j - 1
Loop
RandomShuffle = vtemp
End With
End Function
Appendix – ImanConover Code
A VBA solution to generate correlated numbers with the Iman Conover approach is given here. I have included this code into my sbGenerateTestData application, too.
Please read my Disclaimer.
Function ImanConover(rInputMatrix As Range, _
rTargetCorrelation As Range) As Variant
'Implements the Iman-Conover method to generate random
'number vectors with a given correlation.
'Algorithm as described in:
'Mildenham, November 27, 2005
'Correlation and Aggregate Loss Distributions With An
'Emphasis On The Iman-Conover Method
Dim vX As Variant 'Input matrix
Dim vS As Variant 'Target correlation matrix
Dim vC As Variant 'Cholesky decomposition of vS
Dim vM As Variant 'Intermediate matrix M
Dim vE As Variant 'Covariance matrix E
Dim vF As Variant 'Cholesky decomposition of vE
Dim vT As Variant 'Intermediate matrix T
Dim d As Double, dS As Double
Dim i As Long, j As Long, k As Long
Dim lRow As Long, lCol As Long
With Application.WorksheetFunction
vX = .Transpose(.Transpose(rInputMatrix))
lRow = rInputMatrix.Rows.Count
lCol = rInputMatrix.Columns.Count
'#############################################################################
'# Check inputs #
'#############################################################################
If lCol <> rTargetCorrelation.Columns.Count _
And rTargetCorrelation.Rows.Count <> rTargetCorrelation.Columns.Count Then
'Structure of target correlation matrix needs to fit input matrix
ImanConover = CVErr(xlErrNum)
Exit Function
End If
vS = .Transpose(.Transpose(rTargetCorrelation))
For i = 1 To lCol
If vS(i, i) <> 1# Then
'Target correlation matrix not 1 on diagonal
ImanConover = CVErr(xlErrValue)
Exit Function
End If
For j = 1 To i - 1
If vS(i, j) <> vS(j, i) Then
'Target correlation matrix not symmetric
ImanConover = CVErr(xlErrValue)
Exit Function
End If
Next j
Next i
vC = .Transpose(Cholesky2(vS))
'#############################################################################
'# Create intermediate matrix M #
'#############################################################################
ReDim vMV(1 To lRow) As Double
d = 0#
dS = 0#
For i = 1 To Int(lRow / 2)
vMV(i) = .NormSInv(i / (lRow + 1))
vMV(lRow - i + 1) = -vMV(i)
d = d + 2# * vMV(i) * vMV(i)
Next i
If lRow Mod 2 = 1 Then vMV((lRow + 1) / 2) = 0 'Just for clarity, it's already 0
d = Sqr(d / lRow)
For i = 1 To lRow
vMV(i) = vMV(i) / d
Next i
vM = vX
For i = 1 To lRow
vM(i, 1) = vMV(i)
Next i
Dim vMW As Variant
For i = 2 To lCol
vMW = RandomShuffle2(vMV)
For j = 1 To lRow
vM(j, i) = vMW(j)
Next j
Next i
'#############################################################################
'# Calculate covariance matrix E #
'#############################################################################
vE = vC
For i = 1 To lCol
vE(i, i) = .Covar(.Index(.Transpose(vM), i), .Index(.Transpose(vM), i))
For j = i + 1 To lCol
vE(i, j) = .Covar(.Index(.Transpose(vM), i), .Index(.Transpose(vM), j))
vE(j, i) = vE(i, j)
Next j
Next i
vF = .Transpose(Cholesky2(vE))
vT = .MMult(.MMult(vM, .MInverse(vF)), vC)
'#############################################################################
'# Compute ranks of matrix T #
'#############################################################################
Dim vRT As Variant, vR As Variant
vRT = vX
For j = 1 To lCol
vR = IndexX(lRow, vT, j)
For i = 1 To lRow
vRT(i, j) = vR(i)
Next i
vR = IndexX(lRow, vX, j)
For i = 1 To lRow
vX(i, j) = vX(vR(i), j)
Next i
Next j
'#############################################################################
'# Calculate result matrix Y #
'#############################################################################
Dim vY As Variant
vY = vX
For i = 1 To lRow
For j = 1 To lCol
vY(i, j) = vX(vRT(i, j), j)
Next j
Next i
ImanConover = vY
End With
End Function
Function IndexX(n As Long, arr As Variant, colNo As Long) As Variant
'Indexes an array arr[1..n], i.e., outputs the array indx[1..n] such
'that arr[indx[j]] is in ascending order for j = 1, 2, . . . ,n. The
'input quantities n and arr are not changed. Translated from [31].
Const m As Long = 7
Const NSTACK As Long = 50
Dim i As Long, indxt As Long, ir As Long, itemp As Long, j As Long
Dim k As Long, l As Long
Dim jstack As Long, istack(1 To NSTACK) As Long
Dim a As Double
ir = n
l = 1
ReDim indx(1 To n) As Long
For j = 1 To n
indx(j) = j
Next j
Do While 1
If (ir - l < m) Then
For j = l + 1 To ir
indxt = indx(j)
a = arr(indxt, colNo)
For i = j - 1 To l Step -1
If (arr(indx(i), colNo) <= a) Then Exit For
indx(i + 1) = indx(i)
Next i
indx(i + 1) = indxt
Next j
If (jstack = 0) Then Exit Do
ir = istack(jstack)
jstack = jstack - 1
l = istack(jstack)
jstack = jstack - 1
Else
k = (l + ir) / 2
itemp = indx(k)
indx(k) = indx(l + 1)
indx(l + 1) = itemp
If (arr(indx(l), colNo) > arr(indx(ir), colNo)) Then
itemp = indx(l)
indx(l) = indx(ir)
indx(ir) = itemp
End If
If (arr(indx(l + 1), colNo) > arr(indx(ir), colNo)) Then
itemp = indx(l + 1)
indx(l + 1) = indx(ir)
indx(ir) = itemp
End If
If (arr(indx(l), colNo) > arr(indx(l + 1), colNo)) Then
itemp = indx(l)
indx(l) = indx(l + 1)
indx(l + 1) = itemp
End If
i = l + 1
j = ir
indxt = indx(l + 1)
a = arr(indxt, colNo)
Do While 1
Do
i = i + 1
Loop While (arr(indx(i), colNo) < a)
Do
j = j - 1
Loop While (arr(indx(j), colNo) > a)
If (j < i) Then Exit Do
itemp = indx(i)
indx(i) = indx(j)
indx(j) = itemp
Loop
indx(l + 1) = indx(j)
indx(j) = indxt
jstack = jstack + 2
If (jstack > NSTACK) Then
'STACK too small in indexx
IndexX = CVErr(xlErrNum)
Exit Function
End If
If (ir - i + 1 >= j - l) Then
istack(jstack) = ir
istack(jstack - 1) = i
ir = j - 1
Else
istack(jstack) = j - 1
istack(jstack - 1) = l
l = i
End If
End If
Loop
IndexX = indx
End Function
Function Cholesky2(vA As Variant) As Variant
'Array version
Dim d As Double
Dim i As Long, j As Long, k As Long, n As Long
n = UBound(vA, 1)
If n <> UBound(vA, 2) Then
Cholesky2 = CVErr(xlErrRef)
Exit Function
End If
ReDim vR(1 To n, 1 To n) As Variant
For j = 1 To n
d = 0#
For k = 1 To j - 1
d = d + vR(j, k) * vR(j, k)
Next k
vR(j, j) = vA(j, j) - d
If vR(j, j) <= 0 Then
Cholesky2 = CVErr(xlErrNum)
Exit Function
End If
vR(j, j) = Sqr(vR(j, j))
For i = j + 1 To n
d = 0#
For k = 1 To j - 1
d = d + vR(i, k) * vR(j, k)
Next k
vR(i, j) = (vA(i, j) - d) / vR(j, j)
Next i
For i = 1 To j - 1
vR(i, j) = 0
Next i
Next j
Cholesky2 = vR
End Function
Download
Please read my Disclaimer.
mildenhall_example_on_iman_conover.xlsm [53 KB Excel file, open and use at your own risk]