- Create Cholesky Decomposition matrix A of input correlation matrix
- Generate a vector of independent normal random variables Z
- Compute a vector of correlated normal random variables by using Cholesky matrix X = AZ
- Convert X back to uniform plane [0,1] to get the matrix U containing correlated uniform random variables
- Use matrix U and Inverse Transform Sampling to generate correlated random variables for any marginal distributions.
Main program flow
Director (main program) creates parameter wrapper (Dictionary data structure) for Copula. Director creates correlation matrix for Copula as Matrix object (custom data structure). Director creates random number generator implementation for Copula. Director sets the number of simulations and uniform transformation condition for Copula. Director creates Copula implementation as Gaussian Copula and initializes the model with required data and objects (init).
Copula model gets independent normal random numbers as Matrix object from Random implementation (aggregated in Copula). Copula uses Cholesky decomposition for creating correlated normal random numbers. Copula transforms simulated correlated normal random numbers into uniform plane (optional). Finally, director gets correlated random numbers from Copula as Matrix object.
Generator for random numbers
Readers might be aware, that random number quality produced by Excel Rand function has been reported to be insufficient. Also, the use of Excel NormsInv function is irritatingly slow. For this implementation, the more efficient tools have been employed.
The source for an algorithm implementation of Mersenne Twister was found in this page. Download a zip file from the page and look for mt19937.dll. Save this dll file into C:\temp folder. The example program presented below has been configured so, that it searches that dll file from that folder. The source code for an algorithm for computing the inverse normal cumulative distribution function was acquired from Peter Acklam's web page.
Example program
The following program is a direct implementation for the UML presented above. Director is the main program in VBA (tester). At this point, we should reference the library for Dictionary data structure (VB editor - Tools - References - Microsoft Scripting Runtime).
CopyPaste enumerators into new VBA standard module (name = Enumerators)
Option Explicit ' Public Enum E_ P_SIMULATIONS = 1 P_GENERATOR_TYPE = 2 P_CORRELATION_MATRIX = 3 P_TRANSFORM_TO_UNIFORM = 4 End Enum '
Next, we create the actual tester program (director). CopyPaste the following code into standard VBA module (name = MainProgram). There are some source data read from Excel worksheet Sheet1 in the program presented below. Set correlation matrix into Excel and give a name for that range ("_correlation"). Similarly, give range names for number of simulations ("_simulations") and uniform transform condition ("_transform"). Also, set up a range name for output ("_dataDump").
Option Explicit ' Public Sub tester() ' ' create correlation matrix object Dim correlation As New Matrix correlation.matrixFromRange (Sheets("Sheet1").Range("_correlation")) ' ' create parameters for copula Dim parameters As New Scripting.Dictionary parameters.Add P_SIMULATIONS, CLng(Sheets("Sheet1").Range("_simulations").value) parameters.Add P_GENERATOR_TYPE, New MersenneTwister parameters.Add P_CORRELATION_MATRIX, correlation parameters.Add P_TRANSFORM_TO_UNIFORM, CBool(Sheets("Sheet1").Range("_transform").value) ' ' create copula implementation Dim copula As ICopula: Set copula = New GaussianCopula copula.init parameters ' ' get results from copula and write these into Excel Dim result As Matrix: Set result = copula.getMatrix result.matrixToRange Sheets("Sheet1").Range("_dataDump") ' ' release objects Set result = Nothing Set copula = Nothing Set parameters = Nothing Set correlation = Nothing End Sub '
After this, we create ICopula interface. CopyPaste the following code into a new VBA class module (name = ICopula).
Option Explicit ' ' interface for copula model ' Public Function init(ByRef parameters As Scripting.Dictionary) ' interface - no implementation End Function ' Public Function getMatrix() As Matrix ' interface - no implementation End Function '
CopyPaste the following ICopula implementation into a new VBA class module (name = GaussianCopula).
Option Explicit ' Implements ICopula ' Private n As Long ' number of simulations Private transform As Boolean ' condition for uniform transformation Private generator As IRandom ' random number generator implementation ' Private c As Matrix ' correlation matrix Private d As Matrix ' cholesky decomposition matrix Private z As Matrix ' independent normal random variables Private x As Matrix ' correlated normal random variables ' Private Function ICopula_init(ByRef parameters As Scripting.Dictionary) ' ' initialize class data and objects n = parameters(P_SIMULATIONS) transform = parameters(P_TRANSFORM_TO_UNIFORM) Set generator = parameters(P_GENERATOR_TYPE) Set c = parameters(P_CORRELATION_MATRIX) End Function ' Private Function ICopula_getMatrix() As Matrix ' ' create matrix of independent normal random numbers Set z = New Matrix: z.init n, c.get_c Set z = generator.getNormalRandomMatrix(z.get_r, z.get_c) ' ' create cholesky decomposition Set d = New Matrix: Set d = d.cholesky(c) ' ' create correlated normal random numbers z.transpose Set x = New Matrix: Set x = x.multiplication(d, z) x.transpose ' ' transform correlated normal random numbers ' into correlated uniform random numbers If (transform) Then uniformTransformation Set ICopula_getMatrix = x End Function ' Private Function uniformTransformation() ' ' map normal random number to uniform plane Dim nRows As Long: nRows = x.get_r Dim nCols As Long: nCols = x.get_c ' Dim i As Long, j As Long For i = 1 To nRows For j = 1 To nCols x.push i, j, WorksheetFunction.NormSDist(x.at(i, j)) Next j Next i End Function '
Then, we create interface for random number generator. CopyPaste the following code into a new VBA class module (name = IRandom).
Option Explicit ' Public Function getNormalRandomMatrix( _ ByVal nRows As Long, _ ByVal nCols As Long) As Matrix ' ' interface - no implementation ' takes in two parameters (number of rows and columns) and ' returns matrix object filled with normal random variates End Function '
Next, we create implementation for IRandom. CopyPaste the following code into a new VBA class module (name = MersenneTwister).
Option Explicit ' Implements IRandom ' Private Declare Function nextMT Lib "C:\temp\mt19937.dll" Alias "genrand" () As Double ' Private Function IRandom_getNormalRandomMatrix( _ ByVal nRows As Long, _ ByVal nCols As Long) As Matrix ' ' retrieve NxM matrix with normal random numbers Dim r As Matrix: Set r = New Matrix: r.init nRows, nCols Dim i As Long, j As Long For i = 1 To nRows For j = 1 To nCols r.push i, j, InverseCumulativeNormal(nextMT()) Next j Next i ' Set IRandom_getNormalRandomMatrix = r End Function ' Public Function InverseCumulativeNormal(ByVal p As Double) As Double ' ' Define coefficients in rational approximations Const a1 = -39.6968302866538 Const a2 = 220.946098424521 Const a3 = -275.928510446969 Const a4 = 138.357751867269 Const a5 = -30.6647980661472 Const a6 = 2.50662827745924 ' Const b1 = -54.4760987982241 Const b2 = 161.585836858041 Const b3 = -155.698979859887 Const b4 = 66.8013118877197 Const b5 = -13.2806815528857 ' Const c1 = -7.78489400243029E-03 Const c2 = -0.322396458041136 Const c3 = -2.40075827716184 Const c4 = -2.54973253934373 Const c5 = 4.37466414146497 Const c6 = 2.93816398269878 ' Const d1 = 7.78469570904146E-03 Const d2 = 0.32246712907004 Const d3 = 2.445134137143 Const d4 = 3.75440866190742 ' 'Define break-points Const p_low = 0.02425 Const p_high = 1 - p_low ' 'Define work variables Dim q As Double, r As Double ' 'If argument out of bounds, raise error If p <= 0 Or p >= 1 Then Err.Raise 5 ' If p < p_low Then ' 'Rational approximation for lower region q = Sqr(-2 * Log(p)) InverseCumulativeNormal = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _ ((((d1 * q + d2) * q + d3) * q + d4) * q + 1) ' ElseIf p <= p_high Then 'Rational approximation for lower region q = p - 0.5 r = q * q InverseCumulativeNormal = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / _ (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1) ' ElseIf p < 1 Then 'Rational approximation for upper region q = Sqr(-2 * Log(1 - p)) InverseCumulativeNormal = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _ ((((d1 * q + d2) * q + d3) * q + d4) * q + 1) End If End Function '
Finally, we create our custom data structure class, which is used extensively in this design. CopyPaste the following code into a new VBA class module (name = Matrix).
Option Explicit ' ' general matrix data structure for double data type ' variant array of double arrays, works like Excel Private outer() As Variant Private r_ As Long Private c_ As Long ' Public Function init(ByVal r As Long, ByVal c As Long) ' ' create new matrix object r_ = r c_ = c ReDim outer(1 To r_) ' Dim i As Long For i = 1 To r_ Dim inner() As Double: ReDim inner(1 To c_) outer(i) = inner Next i End Function ' Public Function multiplication(ByRef m1 As Matrix, ByRef m2 As Matrix) As Matrix ' ' get matrix multiplication from two external matrix objects ' return a new matrix object Dim result As New Matrix: result.init m1.get_r, m2.get_c Dim i As Long, j As Long, k As Long ' Dim r1 As Long: r1 = m1.get_r Dim c1 As Long: c1 = m1.get_c Dim r2 As Long: r2 = m2.get_r Dim c2 As Long: c2 = m2.get_c Dim v As Double ' For i = 1 To r1 For j = 1 To c2 v = 0 ' For k = 1 To c1 v = v + m1.at(i, k) * m2.at(k, j) Next k result.push i, j, v Next j Next i Set multiplication = result End Function ' Public Function clone(ByRef m As Matrix) ' ' clone external matrix object into Me object - no return matrix Dim nRows As Long: nRows = m.get_r Dim nCols As Long: nCols = m.get_c Me.set_r nRows Me.set_c nCols Me.setOuter m.getOuter End Function ' Public Function transpose() ' ' transpose matrix (Me) - no return matrix Dim nRows As Long: nRows = Me.get_r Dim nCols As Long: nCols = Me.get_c Dim m As New Matrix: m.init nCols, nRows ' Dim i As Long, j As Long For i = 1 To nRows For j = 1 To nCols 'result.push j, i, m.at(i, j) m.push j, i, Me.at(i, j) Next j Next i Me.clone m End Function ' Public Function cholesky(ByRef c As Matrix) As Matrix ' ' create cholesky decomposition, a lower triangular matrix ' d = decomposition, c = correlation matrix ' return a new matrix object Dim s As Double Dim n As Long: n = c.get_r Dim m As Long: m = c.get_c Dim d As New Matrix: d.init n, m ' Dim i As Long, j As Long, k As Long For j = 1 To n s = 0 For k = 1 To j - 1 s = s + d.at(j, k) ^ 2 Next k d.push j, j, c.at(j, j) - s If d.at(j, j) <= 0 Then Exit For d.push j, j, Sqr(d.at(j, j)) ' For i = j + 1 To n s = 0 For k = 1 To j - 1 s = s + d.at(i, k) * d.at(j, k) Next k d.push i, j, (c.at(i, j) - s) / d.at(j, j) Next i Next j Set cholesky = d End Function ' Public Function matrixToRange(ByRef r As Range) ' ' write matrix content (Me) into Excel range Dim nRows As Long: nRows = Me.get_r Dim nCols As Long: nCols = Me.get_c r.ClearContents ' Dim i As Long, j As Long For i = 1 To nRows For j = 1 To nCols r(i, j) = Me.at(i, j) Next j Next i End Function ' Public Function matrixFromRange(ByRef r As Range) ' ' initialize and copy matrix (Me) from Excel range Dim m As Variant: m = r.Value2 r_ = UBound(m, 1) c_ = UBound(m, 2) Me.init r_, c_ ' Dim i As Long, j As Long For i = 1 To r_ For j = 1 To c_ Me.push i, j, m(i, j) Next j Next i End Function ' Public Function push(ByVal r As Long, ByVal c As Long, ByVal value As Double) ' set value for matrix item outer(r)(c) = value End Function ' Public Function at(ByVal r As Long, ByVal c As Long) As Double ' get value from matrix at = outer(r)(c) End Function ' Public Function getOuter() As Variant ' get variant array getOuter = outer End Function ' Public Function setOuter(ByRef v As Variant) ' set variant array outer = v End Function ' Public Function get_r() As Variant ' get number of matrix rows get_r = r_ End Function ' Public Function set_r(ByVal r As Long) ' set number of matrix rows r_ = r End Function ' Public Function get_c() As Long ' get number of matrix columns get_c = c_ End Function ' Public Function set_c(ByVal c As Long) ' set number of matrix columns c_ = c End Function '
After creating all the previous components in VBA, the program is ready for the use.
Results
Correlated random numbers has been simulated for bi-variate case. The following scatter graph shows the results for 1000 simulated correlated normal random numbers (rho = 0.76). Note, that in Copula implementation, class member transform is FALSE and hereby, uniform transformation is not performed.
Setting class member transform to be TRUE, performs uniform transformation and the results are plotted within the following scatter graph.
So, this Copula implementation is leaving an option for its user to receive correlated random numbers as normal, or receive these numbers mapped into uniform plane.
The latter scheme might be useful, if we need to simulate correlated random numbers from any other distributions. For this task, we can then use inverse transform sampling. As an example, the following scatter chart is showing the results for 1000 simulated correlated exponential random numbers (rho = 0.76, lambda = 0.085).
Correlated normal random numbers were first mapped into uniform plane (setting class member transform to be TRUE) and then transformed to exponentially distributed correlated random numbers with the inverse cumulative distribution function.
Afterthoughts
Presented design for VBA is relatively easy to implement. Also, it has a lot of flexibility. Say, we would like to create implementation for Student Copula, we would just have to create new implementation from ICopula interface. Doing this does not have any changes to be made into existing program design, since we are programming to an interface, not to an implementation. The same applies for generating random numbers. It is now easy to create your own generator for any existing Copula design, if you some day manage to create something better than MT algorithm. A small performance penalty is paid with all the function calls made using Matrix object (which is technically just wrapping arrays into a manageable class).
Creating numerical algorithms and solutions is extremely interesting and rewarding. However, there are a lot of time spent for testing and debugging and still, after all efforts, you might feel like walking on thin ice sometimes. Self-made algorithms could be sometimes unreliably, unstable or not accurate enough and for that reason, doing "anything serious" with Copulas (such as creating production pricing tools) I would recommend to check numerical libraries with "proven industrial strength", such as NAG by Numerical Algorithms Group. Along with the existing libraries for C++/C#, NAG has also Fortran library, which can be easily used also with VBA.
Needless to say, of course, that we get the same results with just a few lines of code with the tools like Matlab or Python. However, the interest of this blog is on the development side of programs. My view is that, as a developer, implementing such algorithms are excellent way to learn in order to become a better developer. However, this is only my personal view on the matter.
Anyway, that's all I wanted to share this time. Again, I hope that this post could help you to solve some of your programming problems with VBA. Happy New Year for everybody!
-Mike
Hello sir,
ReplyDeleteNice to have a detailed code like this. But I am just a undergrad student , quite difficult for me to understand the whole code. Can you please mail me the excel file. Or can we have a skype session , to talk on certain queries of mine. I opted CQF but struggling with VBA. EMAIL- anash4444@gmail.com
I have sent a working example for you. Work with that one. Good luck.
Deleteujjwal bhateja could u pls send me the excel file?
Deletesonia.goretti2710@gmail.com
DeleteCan you also send the example file to ates_h@ymail.com as well. Thnx in advance
ReplyDeleteI have sent a working example also for you. Good luck.
DeleteI can get the example. Thank you very much for sharing knowledge.
ReplyDeleteEmail: jvr777m@gmail.com
Can i also pls get an example? kav.ghai@gmail.com
ReplyDeleteCan I also get an example sheet please.
ReplyDeletesunnykalra1087@gmail.com
Thanks a lot
Can i also please get an example?
ReplyDeleteluongbui@gmx.de
thanks alot
If it is possible, kindly please share the example excel file with me?
ReplyDeletelemonade22112@yahoo.co.kr
Many thanks in advance.
Can I get an example too please.
ReplyDeletemymacros@hotmail.com
Many thanks for sharing
Very nice example. Can you send the sample to g123456789l@gmail.com
ReplyDeleteSir, Nice explanation.. can you please send me a sample code?
ReplyDeleteJoining the other ones: would you be able to share the excel and send it to schuebe1971@yahoo.de? Many thanks.
ReplyDeleteDear All, as soon as I am back in Finland, I will send this Excel for all of you who have requested that. Thanks for your interest. -Mike
ReplyDeleteMore specifically at 26.6.2017
DeleteWould mine you give and Excel example to me, uta.kholiq@gmail.com
DeleteI have sent Excel for those who requested that. Sorry for the delay.
ReplyDeleteHi Mikael,
Deletethat is wonderful. Would you be so kind sharing the spreadsheet with me too?
May email: samsung1492@outlook.de
Thank you so much for sharing your experiences!
Hello sir, can you send me please the Excel file on :
Deletemohamed.whibi@outlook.com
Hi Mikael,
ReplyDeleteany chance to send me the excel file and /or C++ related script to makyou84@yahoo.com. I can't find anything about simulation of correlated stock returns based on 5 years past history of daily return with capitalizaion weights added to. Not sure you have something helpful on that as well. Many thanks in advance, for the nice blog too.
Hi Mikael,
ReplyDeleteMany thanks for this piece, that's very useful. Would you please be willing to share the spreadsheet as well? clement.larrue@edhec.com .
Great program. How can I get the file with the code? Thank you very much for giving knowledge. jvrico@yahoo.com
ReplyDeleteHi Mikael, this is fantastic! Would you mind sharing the spreadsheet with me as well? My E-Mail address is bekkers@live.co.za. Thank you so much for sharing your Knowledge and experience!
ReplyDeleteHi Mikael,
ReplyDeletethat is wonderful. Would you be so kind sharing the spreadsheet with me too?
May email: samsung1492@outlook.de
Thank you so much for sharing your experiences!
Hi Mikael, I'd be eternally grateful if you could please share the spreadsheet with me too? Uksurfer92@gmail.com thank you for sharing the knowledge, truly
ReplyDeleteHello Mikael: Would you be kind enough to share the Excel workbook with me as well. If you have an excel workbook with the student copula, I would greatly appreciate if you could share it with me.
ReplyDeletemy e-mail: mfritzner@hotmail.com
Thanks,
Fritz
Hi Mikael,
ReplyDeletewould you be so kind sharing the spreadsheet with me too?
My E-Mail: stevan@kanjo.eu
Could You pls sent the excel file to me too?
ReplyDeleteMail: sonia.goretti2710@gmail.com
Hi Michael,
ReplyDeleteCan you share the example Excel spreadsheet with me (gkchang@gmail.com)!? TIA!
Thanks a bunch,
George
Hi Mikael,
ReplyDeletewould you be so kind sharing the spreadsheet with me.. thanks..
My E-Mail: sgainsbourg2@gmail.com
Hi Mikael,
ReplyDeleteCould you plz send me the excel file?
My email: thang2379@gmail.com
Many thanks!
Hi Mike!
ReplyDeleteCould you please send me the excel as well?
Thank you: bidain16@bu.edu
Hey Mike
ReplyDeleteCould I ask you to send me the file as well? I've been reading up on CDOs and would like a better understanding of how to price them. I've read the theory, starting with the basic one factor model and making my way to using Copulas to price them. I would really appreciate seeing the theory put to use in an excel file.
Thanks.
sinpan@gmail.com
Hi Mikael,
ReplyDeleteCould you please send me the excel as well? Thank you.
christoph.hofstetter.1990@gmail.com
Hi Mikael,
ReplyDeleteI believe your excel file would open up my understanding on this topic, Could you please send me the excel as well? hluish@gmail.com
Hi Mikael,
ReplyDeleteMay I know if you could send me the excel? Thank you.
marcooo13@gmail.com
Hi Mikeel, It is fantastic! Could you please send me the excel with your VBA code? Thank you very much.My email's thangdh73@gmail.com
ReplyDeleteHi Mikael - could I have a copy also - many thanks peter.urbani@gmail.com
ReplyDelete