This example uses method Random.NextGaussianCopula
to generate a multivariate sequence GCdevt[k=0..nseq-1][j=0..nvar-1]
whose marginal distributions are user-defined and imprinted with a user-specified correlation matrix CorrMtrxIn[i=0..nvar-1][j=0..nvar-1]
and then uses method Random.CanonicalCorrelation
to extract from this multivariate random sequence a canonical correlation matrix CorrMtrx[i=0..nvar-1][j=0..nvar-1]
.
This example illustrates two useful copula related procedures. The first procedure generates a random multivariate sequence with arbitrary user-defined marginal deviates whose dependence is specified by a user-defined correlation matrix. The second procedue is the inverse of the first: an arbitary multivariate deviate input sequence is first mapped to a corresponding sequence of empirically derived variates, i.e. cumulative distribution function values representing the probability that each random variable has a value less than or equal to the input deviate. The variates are then inverted, using the inverse Normal(0,1) function, to N(0,1) deviates; and finally, a canonical covariance matrix is extracted from the multivariate N(0,1) sequence using the standard sum of products.
This example demonstrates that the nextGaussianCopula
method correctly imbeds the user-defined correlation information into an arbitrary marginal distribution sequence by extracting the canonical correlation from these sequences and showing that they differ from the original correlation matrix by a small relative error, which generally decreases as the number of multivariate sequence vectors increases.
using System;
using Imsl.Math;
using Imsl.Stat;
public class RandomEx2
{
internal static Imsl.Stat.Random IMSLRandom()
{
Imsl.Stat.Random r = new Imsl.Stat.Random(123457);
r.Multiplier = 16807;
return r;
}
public static void Main(string[] args)
{
double[,] CorrMtrxIn = new double[,]{
{1.0, - 0.9486832980505138, 0.8164965809277261},
{- 0.9486832980505138, 1.0, -0.6454972243679028},
{0.8164965809277261, -0.6454972243679028, 1.0}};
int nvar = 3;
Console.WriteLine("Random Example 2:");
Console.WriteLine();
for (int i = 0; i < nvar; i++)
{
for (int j = 0; j < i; j++)
{
Console.WriteLine("CorrMtrxIn[" + i + "," + j + "] = " + CorrMtrxIn[i, j]);
}
}
PrintMatrixFormat pmf = new PrintMatrixFormat();
pmf.NumberFormat = "0.000000000";
new PrintMatrix("Input Correlation Matrix: ").Print(pmf, CorrMtrxIn);
Console.WriteLine("Correlation Matrices calculated from");
Console.WriteLine(" Gaussian Copula imprinted multivariate sequence:");
Console.WriteLine();
// Compute the Cholesky factorization of CorrMtrxIn
Cholesky CholMtrx = new Cholesky(CorrMtrxIn);
for (int kmax = 500; kmax < 1000000; kmax *= 10)
{
Console.WriteLine("# vectors in multivariate sequence: " + kmax);
double[][] GCvart = new double[kmax][];
double[][] GCdevt = new double[kmax][];
for (int i2 = 0; i2 < kmax; i2++)
{
GCdevt[i2] = new double[nvar];
}
Imsl.Stat.Random r = IMSLRandom();
for (int k = 0; k < kmax; k++)
{
GCvart[k] = r.NextGaussianCopula(CholMtrx); //probs
for (int j = 0; j < nvar; j++)
{
/*
* invert Gaussian Copula probabilities to deviates using
* variable-specific inversions: j = 0: Chi Square; 1: F;
* 2: Normal(0,1); will end up with deviate sequences ready
* for mapping to canonical correlation matrix:
*/
if (j == 0)
{
//convert probs into ChiSquare(df=10) deviates:
GCdevt[k][j] = InvCdf.Chi(GCvart[k][j], 10.0);
}
else if (j == 1)
{
//convert probs into F(dfn=15,dfd=10) deviates:
GCdevt[k][j] = InvCdf.F(GCvart[k][j], 15.0, 10.0);
}
else
{
//convert probs into Normal(mean=0,variance=1) deviates:
GCdevt[k][j] = InvCdf.Normal(GCvart[k][j]);
}
}
}
/*
* extract Canonical Correlation matrix from arbitrarily distributed
* deviate sequences GCdevt[k=0..kmax-1][j=0..nvar-1] which have been
* imprinted with CorrMtrxIn[i=1..nvar][j=1..nvar] above:
*/
double[][] CorrMtrx = r.CanonicalCorrelation(GCdevt);
double relerr;
for (int i = 0; i < nvar; i++)
{
for (int j = 0; j < i; j++)
{
relerr = Math.Abs(1.0 -
(CorrMtrx[i][j] / CorrMtrxIn[i, j]));
Console.WriteLine("CorrMtrx[" + i + "][" + j + "] = " +
CorrMtrx[i][j] + "; relerr = " + relerr);
}
}
new PrintMatrix("Correlation Matrix: ").Print(pmf, CorrMtrx);
}
}
}
Random Example 2:
CorrMtrxIn[1,0] = -0.948683298050514
CorrMtrxIn[2,0] = 0.816496580927726
CorrMtrxIn[2,1] = -0.645497224367903
Input Correlation Matrix:
0 1 2
0 1.000000000 -0.948683298 0.816496581
1 -0.948683298 1.000000000 -0.645497224
2 0.816496581 -0.645497224 1.000000000
Correlation Matrices calculated from
Gaussian Copula imprinted multivariate sequence:
# vectors in multivariate sequence: 500
CorrMtrx[1][0] = -0.95029565568146; relerr = 0.00169957417218125
CorrMtrx[2][0] = 0.805260514673251; relerr = 0.0137613145197848
CorrMtrx[2][1] = -0.640202740166643; relerr = 0.00820217965529457
Correlation Matrix:
0 1 2
0 1.000000000 -0.950295656 0.805260515
1 -0.950295656 1.000000000 -0.640202740
2 0.805260515 -0.640202740 1.000000000
# vectors in multivariate sequence: 5000
CorrMtrx[1][0] = -0.948611734364931; relerr = 7.54347480657058E-05
CorrMtrx[2][0] = 0.815532446740145; relerr = 0.00118081840157325
CorrMtrx[2][1] = -0.646255361981671; relerr = 0.00117450174090306
Correlation Matrix:
0 1 2
0 1.000000000 -0.948611734 0.815532447
1 -0.948611734 1.000000000 -0.646255362
2 0.815532447 -0.646255362 1.000000000
# vectors in multivariate sequence: 50000
CorrMtrx[1][0] = -0.948314953115713; relerr = 0.00038826965285188
CorrMtrx[2][0] = 0.817827046711005; relerr = 0.00162948114463246
CorrMtrx[2][1] = -0.646669093465958; relerr = 0.00181545180028175
Correlation Matrix:
0 1 2
0 1.000000000 -0.948314953 0.817827047
1 -0.948314953 1.000000000 -0.646669093
2 0.817827047 -0.646669093 1.000000000
# vectors in multivariate sequence: 500000
CorrMtrx[1][0] = -0.94873295551488; relerr = 5.23435634081082E-05
CorrMtrx[2][0] = 0.81728795761655; relerr = 0.000969234540976194
CorrMtrx[2][1] = -0.646929884520875; relerr = 0.00221946756529401
Correlation Matrix:
0 1 2
0 1.000000000 -0.948732956 0.817287958
1 -0.948732956 1.000000000 -0.646929885
2 0.817287958 -0.646929885 1.000000000
Link to C# source.