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

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