I am a supporter of
St. Joseph's hospice.
 If you find this site useful or if it helped you, consider a small donation to
St. Joseph's, please.

Information on
St. Joseph's

JustGiving - Sponsor me now!

 

TestCorrelatedNumbers

Option Explicit

Option Base 1

Sub sbGenCorrData()
'Generate correlated numbers, applying the Iman-Conover method.
'Reverse("moc.LiborPlus.www")
'Change history:
'Version Date        Who   Comment
'0.10    06-Apr-2013 Bernd Initial version
Const CMaxIter = 300

Dim dMinMaxAbsDiff As Double
Dim i As Long, j As Long, k As Long
Dim lCol As Long
Dim lRow As Long
Dim wsIn As Worksheet
Dim wsOut As Worksheet
Dim wsCorr As Worksheet
Dim vR As Variant
Dim state As SystemState

With Application.WorksheetFunction
Set state = New SystemState
Set wsIn = Sheets("Gen_Corr_Input_Series")
Set wsOut = Sheets("Gen_Corr_Output_Series")
Set wsCorr = Sheets("Gen_Corr_Target_Corr")
wsOut.Cells.ClearContents

lRow = wsIn.Range("A1").End(xlDown).Row
lCol = wsIn.Range("A1").End(xlToRight).Column
wsCorr.Cells(lCol + 2, 1).Resize(lCol * 2 + 6, lCol).ClearContents
dMinMaxAbsDiff = 1#

For k = 1 To CMaxIter
    wsOut.Range("A1").Resize(lRow, lCol) = ImanConover(wsIn.Range("A1").Resize(lRow, lCol), _
        wsCorr.Range("A1").Resize(lCol, lCol))
    wsCorr.Cells(lCol + 2, 1) = "Result correlation"
    For i = lCol + 3 To 2 * lCol + 2
        For j = 1 To lCol
            wsCorr.Cells(i, j).FormulaR1C1 = "=CORREL(INDEX(Gen_Corr_Output_Series!R1C1:R" & _
                lRow & "C" & lCol & ",,ROW()-" & lCol + 2 & ")," & _
                "INDEX(Gen_Corr_Output_Series!R1C1:R" & lRow & "C" & lCol & ",,COLUMN()))"
            wsCorr.Cells(lCol + 2 + i, j).FormulaR1C1 = "=R[" & -lCol - 2 & "]C" & " - R[" & _
                -2 * lCol - 4 & "]C"
        Next j
    Next i
    wsCorr.Cells(2 * lCol + 4, 1) = "Correlation difference"
    wsCorr.Cells(3 * lCol + 6, 1) = "Largest absolute error"
    wsCorr.Cells(3 * lCol + 7, 1).FormulaArray = "=MAX(ABS(R[" & -lCol - 2 & "]C:R[-3]C[" & _
        lCol - 1 & "]))"
    wsCorr.Calculate
    If wsCorr.Cells(3 * lCol + 7, 1) < dMinMaxAbsDiff Then
        dMinMaxAbsDiff = wsCorr.Cells(3 * lCol + 7, 1)
        vR = wsOut.Range("A1").Resize(lRow, lCol)
    End If
Next k
If wsCorr.Cells(3 * lCol + 7, 1) <> dMinMaxAbsDiff Then
    wsOut.Range("A1").Resize(lRow, lCol) = vR
    wsCorr.Calculate
End If
End With
End Sub

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

Function RandomShuffle2(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
n = UBound(vtemp, 1)
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
RandomShuffle2 = vtemp
End With
End Function
 

Sulprobil   Get it done   Contact   Disclaimer   Download