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 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 'Uncomment if you think you need this.
With Application.WorksheetFunction
On Error Resume Next 'Ignore error: VBA calls already with 1-dim array.
vtemp = .Transpose(vtemp)
On Error GoTo 0
n = UBound(vtemp)
j = n
Do While j > 0
  u = Rnd()
  k = Int(j * u + 1)
  temp = vtemp(j)
  vtemp(j) = vtemp(k)
  vtemp(k) = 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
'V0.3 PB 02-Nov-2024 by Bernd Plumhoff
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
Dim state As SystemState

With Application
Set state = New SystemState
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(Cholesky(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 = RandomShuffle(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(Cholesky(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