PRODUCING SETS OF CORRELATED NORMAL RANDOM NUMBERS USING MATRIX MATHEMATICS

The previous example allowed us to see the effects of a single type of correlation between entities and to understand the effects on simulation results that correlation exerts. However, the example was limited in scope because the correlation coefficient was provided, and there was only a single correlation coefficient assumed. When financial simulations become more developed, they may require more robust methods of generating correlated normal random numbers.

Using matrix mathematics is an efficient means for creating multiple correlated normal random numbers. In the next Model Builder, we will examine historical data on fictitious companies, perform calculations to understand how those companies are correlated, and then set up a system for generating correlated random numbers for eventual use in a broader financial simulation.

MODEL BUILDER 3.3: Advanced Correlation Concepts through the Lens of Corporate Performance

1. Create a new workbook and save it as MB_3.3_User.xls.
2. For this Model Builder we will need a large set of raw data. Open MB_3.3_Complete.xls from the website, copy the range B3:G505, and paste the copied range in your workbook's first sheet in the same range (B3:G505). We need only one sheet for this Model Builder, so you can delete the other sheets and name the sheet you are working in “MB 3.3” to tie with the complete version.
3. The values that we copied over could represent the observable asset values in a company over time. The first real calculation we should do to understand the correlation is to see how each entity's asset value moves each period. To do this, enter the following formula in I7:

   =C7/C6-1

   Copy and paste this formula over the range I7:M505. Also, for labeling purposes copy over the names of the entities in I3:M3.

4. We now need to set up a system for analyzing the correlation of each entity's asset value changes to one another. This can be done on the sheet or in code. Since it is not too difficult to implement on the sheet using formulas we will use such a method. The first part of this is administrative. Enter 1 to 5 going down the sheet, increasing by 1 in O7:O11. Enter the same series of 1 to 5 going across the sheet, increasing by 1 in Q3:U3. Next to the numbers in O7:O11, enter the names of each entity in P7:P11. For example, P7 should contain TrustedFinCo1, P8 should contain TrustedFinCo2, etc. Similarly in Q6:U6 enter the same names (e.g. Q6 should contain TrustedFinCo1, R6 should contain TrustedFinCo2, etc.).
5. We will now use the CORREL function to calculate the correlations. Recall that the mathematics behind the CORREL function were detailed in Model Builder 3.1. The implementation of the CORREL function in this case is a little bit challenging because we have to compute correlations for every combination of entities. This can be done using the OFFSET function within the CORREL function. The OFFSET function allows us to “move” the range references within the CORREL function as we drag it across. Enter the following function in Q7:

   FIGURE 3.9 We can advance our correlation knowledge by doing multiple correlations.

   

   =CORREL(OFFSET($H$7,0,Q$3):OFFSET($H$505,0,Q$3),OFFSET($H$7,0,$O7):

   OFFSET($H$505,0,$O7))

   Let's examine this formula, since it is complicated. The CORREL function as the outermost function will return the correlation coefficient of two ranges of data. To obtain each range of data we use the OFFSET function. The key to this function is that is allows us to reference fixed locations H7 and H505 as the starting and end points of our data from top to bottom, and then select which range of data we want to use in the calculations, depending on the values in Q3:U3 and O7:O11. Recall that we created a two-dimensional range in step 4 and entered values to reference locations in that range. The OFFSET function uses those values as the formula is dragged. Implement this formula by copying and pasting it over the range Q7:U11. Figure 3.9 shows what the sheet should like so far.

6. After completing the correlation matrix in step 5, take a moment to look at the resulting numbers. You may start to see a pattern. Notice that a diagonal of 1s appears where the CORREL function calculated the correlation coefficient of the entity against itself? Now imagine drawing a line down that diagonal. If one were to do this we would have two triangles: an upper triangle and a lower triangle. What's interesting is that the numbers in the upper triangle mirror those in the lower triangle. For example, the value in Q8 is the same as in R8 and the value in Q9 is the same as in S9. The technical term for what we have created here is a symmetric, positive-definite matrix. Look at Figure 3.10 for further clarity.
7. Once we have our symmetric, positive-definite matrix, we then need to find a matrix where U^TU = C, where U^T signifies the transpose of the matrix U. One method of solving such an equation is using Cholesky decomposition. The best way to implement Cholesky decomposition in Excel is by creating a custom function in code. In the same workbook you have been working in go to the Visual Basic editor and insert a new module. Enter the following code to start the function, declare variables, assign values, and dimension a dynamic range:

   Public Function Cholesky(WorkRange As Range)
   
   Dim WorkArray As Variant

   Dim AggVar As Double
   Dim i As Integer, j As Integer, k As Integer
   
   WorkArray = WorkRange
   rCount = WorkRange.Rows.count
   cCount = WorkRange.Columns.count
   
   ReDim TransArray(1 To rCount, 1 To cCount) As Double

   FIGURE 3.10 A symmetric, positive-definite matrix has unique characteristics.

   

8. The main code for the Cholesky transformation is difficult to explain in a top to bottom order since it is written to work through a matrix with many conditional characteristics. The best way to explain this is to graphically depict a matrix, the location that we are working on, and the code that affects that location. It is recommended that you review the whole code at the end of this chapter or from the MB_3.3_Complete.xls file on the website, since pieces of code will be explained individually.

   The easiest place to start in explaining this function is on the matrix in the upper left corner at row 1, column 1. When the function runs, it skips a section of code and doesn't start calculating until it reaches the following section:

   TransArray(j, j) = WorkArray(j, j)-AggVar
   
   If TransArray(j, j) <=  0 Then
       Exit For
   Else
       TransArray(j, j) = Sqr(TransArray(j, j))
   End If

   At this point it looks to the array on the sheet entered as the parameter and only for the diagonals (row 1, column 1 is not only the starting point, but also a diagonal) and subtracts from it an aggregated variable. The aggregated variable does not exist for row 1, column 1, so the transform array (i.e., TransArray(j, j)) is simply the same value as the referenced array from the sheet. If that value is less than or equal to 0 then we move on to the next row; otherwise we take the square root of it. In this specific example the referenced array at row 1, column 1 has a value of 1, meaning the square root is 1.

9. With the diagonal done the code runs until it hits the following section:

   For i = j + 1 To rCount
     AggVar = 0
   
     For k = 1 To j-1
         AggVar = AggVar + TransArray(j, k) * TransArray(i, k)
     Next k
   
   TransArray(i, j) = (WorkArray(i, j)-AggVar) / TransArray(j, j)
    Next i

   For row 1, column 1 there will be no aggregate variable, so the focus of the code is on filling out TransArray(i, j). The counter variable i will loop through rows 2 to 5 (i.e., effectively a vertical fill of column 1, below the diagonal point we calculated in step 8. Each of these points takes the value in the sheet array and divides it by the transform array's diagonal for the column under consideration.

10. The outermost j loop then continues, which spatially can be thought of as moving to row 2. This time through the code though there are a few differences than when we moved through the first j loop. The initial difference that we encounter is a part of code that we skipped over earlier:

    For j = 1 To rCount
        AggVar = 0
        For k = 1 To j-1
          AggVar = AggVar + TransArray(j, k) ^ 2
        Next k

    Here the k loop will be active and the aggregate variable will take the value of the squared sum of all the points to the left of the diagonal for the row the j loop is working on. See Figure 3.11 for a visual depiction. Looking back at step 8, we can see that when we calculate the diagonal value it is the value from the referenced sheet array, less the sum of the squared values that are to the left of the diagonal of concern. Keep in mind that the final value for the diagonal is the square root of that calculated value, unless it is less than or equal to 0.

    FIGURE 3.11 The Cholesky function transforms the symmetric, positive-definite matrix.

    

11. Another part that is different for the second j loop is when we move down from the diagonal (as in step 9). As you can see in the code referenced at the beginning of step 9, we have an aggregate variable again. As with the first aggregate variable we encountered it calculates values from points to the left of the diagonal of concern; however, this time it takes the value on the transform array immediately to the left of the value being calculated and multiplies that value on the left by the value above the value to the left. This process will be done for every column to the left of the diagonal and for every row below the diagonal.
12. The process continues for each diagonal, looping through steps 8 through 11 until the bottom half of the triangle is complete. Once this is done the final piece of code returns the array to the sheet.

    Cholesky = TransArray
    End Function

13. The final code for Cholesky, in the proper order, should read:

    Public Function Cholesky(WorkRange As Range)
    Dim WorkArray As Variant
    Dim AggVar As Double
    Dim i As Integer, j As Integer, k As Integer
    
    WorkArray = WorkRange
    rCount = WorkRange.Rows.count
    cCount = WorkRange.Columns.count
    
    ReDim TransArray(1 To rCount, 1 To cCount) As Double
    
    For j = 1 To rCount
    
        AggVar = 0
    
        For k = 1 To j-1
            AggVar = AggVar + TransArray(j, k) ^ 2
        Next k
    
        TransArray(j, j) = WorkArray(j, j)-AggVar
    
        If TransArray(j, j) <= 0 Then
            Exit For
        Else
            TransArray(j, j) = Sqr(TransArray(j, j))
        End If
    
        For i = j + 1 To rCount
            AggVar = 0
    
            For k = 1 To j-1
                AggVar = AggVar + TransArray(j, k) * TransArray(i, k)
            Next k

    TransArray(i, j) = (WorkArray(i, j)-AggVar)
                / TransArray(j, j)
        Next i
    Next j
    
    Cholesky = TransArray
    End Function

14. Note that this custom function is built with an array in mind, and when entered on the sheet, Excel must know that this is an array function. Excel identifies functions as array functions when a user presses Control-Shift-Enter after the formula involving the function is complete. Select W7:AA11, enter the formula “= Cholesky(Q7:U11)”, and press Control-Shift-Enter. The sheet should look like Figure 3.11 so far.
15. Even with all of this work we still have more to do. So far we can think of ourselves has having worked on the “correlated” part of creating correlated normal random variables. We should now create the “normal random variables” part. This is much easier given what was learned in Chapter 2. Enter the following formula in AC7:

    =NORMSINV(RAND())

    Copy and paste this to AC7:AG11.

16. The final step is to multiply the entire matrix of normal random variables created in step 15 by the transposed matrix derived from the Cholesky function. This is simplified by using the MMULT and TRANSPOSE functions native to Excel and shown with the formula that should be entered in AI7:AM11:

    =MMULT(AC7:AG11,TRANSPOSE(W7:AA11))

    Note that this is also an array formula and must be entered with Control-Shift-Enter, exactly as the custom Cholesky function must be used. The last section for this Model Builder should look like Figure 3.12.

FIGURE 3.12 Correlated normal random numbers created from asset values and Cholesky decomposition.