NAG G05 Copula wrapper
With wrapper, one can create correlated random numbers
- from any distribution (assuming appropriate inverse transformation method is provided)
- by using Gaussian or Student Copula
Client can create NAGCopula object with two different constructors, with or without inverse transformation. Also, NAGCopula model selection (Gaussian, Student) will be made, based on user-given value for degrees-of-freedom parameter. This parameter, being greater than zero, will automatically trigger the use of Student Copula model.
NAGCopula class, InverseTransformation class hierarchy and Client tester program is presented below. Program will first retrieve a matrix filled with correlated random numbers from NAGCopula object and then writing the content of that matrix into a given text file. Define new addresses for file outputs, if needed. Create a new C# console project, copyPaste the whole content into a new cs file and run the program.
Copula outputs
The results of one client program run is presented in the table below. For a given bi-variate test scheme, the both Copula models (Gaussian, Student) were used to create 5000 correlated random number pairs. The created random numbers were also transformed back to Standard Normal and Exponential distributions.
The Program
using System; using System.Collections.Generic; using System.Linq; using System.Text; using NagLibrary; using System.IO; namespace NAGCopulaNamespace { // CLIENT class Program { static void Main(string[] args) { try { double[,] covariance = CreateCovarianceMatrix(); int[] seed = new int[1]; int n = 5000; double[,] randomMatrix; // // create correlated uniform random numbers using gaussian copula NAGCopula gaussian = new NAGCopula(covariance, 1, 1, 2); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = gaussian.Create(n, seed); Print(@"C:\temp\GaussianCopulaNumbers.txt", randomMatrix); // // create correlated uniform random numbers using student copula NAGCopula student = new NAGCopula(covariance, 1, 1, 2, 3); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = student.Create(n, seed); Print(@"C:\temp\StudentCopulaNumbers.txt", randomMatrix); // // create correlated standard normal random numbers using gaussian copula NAGCopula gaussianTransformed = new NAGCopula(covariance, 1, 1, 2, new StandardNormal()); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = gaussianTransformed.Create(n, seed); Print(@"C:\temp\GaussianTransformedCopulaNumbers.txt", randomMatrix); // // create correlated standard normal random numbers using student copula NAGCopula studentTransformed = new NAGCopula(covariance, 1, 1, 2, new StandardNormal(), 3); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = studentTransformed.Create(n, seed); Print(@"C:\temp\StudentTransformedCopulaNumbers.txt", randomMatrix); // // create correlated exponential random numbers using gaussian copula NAGCopula gaussianTransformed2 = new NAGCopula(covariance, 1, 1, 2, new Exponential(0.078)); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = gaussianTransformed2.Create(n, seed); Print(@"C:\temp\GaussianTransformedCopulaNumbers2.txt", randomMatrix); // // create correlated exponential random numbers using student copula NAGCopula studentTransformed2 = new NAGCopula(covariance, 1, 1, 2, new Exponential(0.078), 3); seed[0] = Math.Abs(Guid.NewGuid().GetHashCode()); randomMatrix = studentTransformed2.Create(n, seed); Print(@"C:\temp\StudentTransformedCopulaNumbers2.txt", randomMatrix); } catch (Exception e) { Console.WriteLine(e.Message); } } static void Print(string filePathName, double[,] arr) { // write generated correlated random numbers into text file StreamWriter writer = new StreamWriter(filePathName); // for (int i = 0; i < arr.GetLength(0); i++) { for (int j = 0; j < arr.GetLength(1); j++) { writer.Write("{0};", arr[i, j].ToString().Replace(',','.')); } writer.WriteLine(); } writer.Close(); } static double[,] CreateCovarianceMatrix() { // hard-coded data for 2x2 co-variance matrix // correlation is (roughly) 0.75 double[,] c = new double[2, 2]; c[0, 0] = 0.00137081575700264; c[0, 1] = 0.00111385579983184; c[1, 0] = 0.00111385579983184; c[1, 1] = 0.00161530858115086; return c; } } // // // NAG COPULA WRAPPER public class NAGCopula { private InverseTransformation transformer; private double[,] covariance; private int genID; private int subGenID; private int mode; private int df; private int m; private int errorNumber; private double[] r; private double[,] result; // public NAGCopula(double[,] covariance, int genID, int subGenID, int mode, int df = 0) { // ctor : create correlated uniform random numbers // degrees-of-freedom parameter (df), being greater than zero, // will automatically trigger the use of student copula this.covariance = covariance; this.genID = genID; this.subGenID = subGenID; this.mode = mode; this.df = df; this.m = covariance.GetLength(1); r = new double[m * (m + 1) + (df > 0 ? 2 : 1)]; } public NAGCopula(double[,] covariance, int genID, int subGenID, int mode, InverseTransformation transformer, int df = 0) { // ctor : create correlated random numbers from distribution X // degrees-of-freedom parameter (df), being greater than zero, // will automatically trigger the use of student copula this.transformer = transformer; this.covariance = covariance; this.genID = genID; this.subGenID = subGenID; this.mode = mode; this.df = df; this.m = covariance.GetLength(1); r = new double[m * (m + 1) + (df > 0 ? 2 : 1)]; } public double[,] Create(int n, int[] seed) { result = new double[n, m]; // G05.G05State g05State = new G05.G05State(genID, subGenID, seed, out errorNumber); if (errorNumber != 0) throw new Exception("G05 state failure"); // if (this.df != 0) { // student copula G05.g05rc(mode, n, df, m, covariance, r, g05State, result, out errorNumber); } else { // gaussian copula G05.g05rd(mode, n, m, covariance, r, g05State, result, out errorNumber); } // if (errorNumber != 0) throw new Exception("G05 method failure"); if (transformer != null) transform(); return result; } private void transform() { // map uniform random numbers to distribution X for (int i = 0; i < result.GetLength(0); i++) { for (int j = 0; j < result.GetLength(1); j++) { result[i, j] = transformer.Inverse(result[i, j]); } } } } // // // INVERSE TRANSFORMATION LIBRARY // base class for all inverse transformation classes public abstract class InverseTransformation { public abstract double Inverse(double x); } // // mapping uniform random variate to exponential distribution public class Exponential : InverseTransformation { private double lambda; public Exponential(double lambda) { this.lambda = lambda; } public override double Inverse(double x) { return Math.Pow(-lambda, -1.0) * Math.Log(x); } } // // mapping uniform random variate to standard normal distribution public class StandardNormal : InverseTransformation { public override double Inverse(double x) { const double a1 = -39.6968302866538; const double a2 = 220.946098424521; const double a3 = -275.928510446969; const double a4 = 138.357751867269; const double a5 = -30.6647980661472; const double a6 = 2.50662827745924; // const double b1 = -54.4760987982241; const double b2 = 161.585836858041; const double b3 = -155.698979859887; const double b4 = 66.8013118877197; const double b5 = -13.2806815528857; // const double c1 = -7.78489400243029E-03; const double c2 = -0.322396458041136; const double c3 = -2.40075827716184; const double c4 = -2.54973253934373; const double c5 = 4.37466414146497; const double c6 = 2.93816398269878; // const double d1 = 7.78469570904146E-03; const double d2 = 0.32246712907004; const double d3 = 2.445134137143; const double d4 = 3.75440866190742; // const double p_low = 0.02425; const double p_high = 1.0 - p_low; double q = 0.0; double r = 0.0; double c = 0.0; // if ((x > 0.0) && (x < 1.0)) { if (x < p_low) { // rational approximation for lower region q = Math.Sqrt(-2.0 * Math.Log(x)); c = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } if ((x >= p_low) && (x <= p_high)) { // rational approximation for middle region q = x - 0.5; r = q * q; c = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1); } if (x > p_high) { // rational approximation for upper region q = Math.Sqrt(-2 * Math.Log(1 - x)); c = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } } else { // throw if given x value is out of bounds // (less or equal to 0.0, greater or equal to 1.0) throw new Exception("input parameter out of bounds"); } return c; } } }
Afterthoughts
It is quite easy to appreciate the numerical work performed by NAG G05 library. In the case you might be interested about how to create non-correlated random numbers with NAG, check out my blog post in here. If you like the stuff, you can get yourself more familiar with NAG libraries in here. So, that was all I wanted to share with you again. Thanks for reading my blog. Good night. Mike.
No comments:
Post a Comment