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.
import com.imsl.stat.*;
import com.imsl.math.*;
import com.imsl.math.PrintMatrix;
public class RandomEx2 {
static Random IMSLRandom() {
Random r = new Random();
r.setSeed(123457);
r.setMultiplier(16807);
return r;
}
public static void main(String args[]) throws com.imsl.IMSLException
{
double CorrMtrxIn[][] = {
{ 1., -0.9486832980505138, 0.8164965809277261},
{-0.9486832980505138, 1., -0.6454972243679028},
{ 0.8164965809277261, -0.6454972243679028, 1.}
};
int nvar = 3;
System.out.println("Random Example 2:");
System.out.println("");
for (int i = 0; i < nvar; i++) {
for (int j = 0; j < i; j++) {
System.out.println("CorrMtrxIn["+i+"]["
+j+"] = " + CorrMtrxIn[i][j]);
}
}
PrintMatrixFormat pmf = new PrintMatrixFormat();
pmf.setNumberFormat(new java.text.DecimalFormat("0.000000000"));
new PrintMatrix("Input Correlation Matrix: ").print(pmf, CorrMtrxIn);
System.out.println("Correlation Matrices calculated from");
System.out.println(" Gaussian Copula imprinted multivariate sequence:");
System.out.println("");
// Compute the Cholesky factorization of CorrMtrxIn
Cholesky CholMtrx = new Cholesky(CorrMtrxIn);
for (int kmax = 500; kmax < 1000000; kmax *= 10) {
System.out.println("# vectors in multivariate sequence: "+ kmax);
double GCvart[][] = new double[kmax][];
double GCdevt[][] = new double[kmax][nvar];
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.);
}
else if (j == 1) {
//convert probs into F(dfn=15,dfd=10) deviates:
GCdevt[k][j] = InvCdf.F(GCvart[k][j], 15., 10.);
}
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. - (CorrMtrx[i][j]/CorrMtrxIn[i][j]));
System.out.println("CorrMtrx["+i+"]["+j+"] = " + CorrMtrx[i][j]
+ "; relerr = " + relerr);
}
}
new PrintMatrix("Correlation Matrix: ").print(pmf, CorrMtrx);
}
}
}
Random Example 2:
CorrMtrxIn[1][0] = -0.9486832980505138
CorrMtrxIn[2][0] = 0.8164965809277261
CorrMtrxIn[2][1] = -0.6454972243679028
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.9502956556814602; relerr = 0.0016995741721812507
CorrMtrx[2][0] = 0.8052605146732508; relerr = 0.013761314519784795
CorrMtrx[2][1] = -0.6402027401666432; relerr = 0.008202179655294572
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.9486118836203709; relerr = 7.527741901824925E-5
CorrMtrx[2][0] = 0.815532446740145; relerr = 0.001180818401573247
CorrMtrx[2][1] = -0.6462558369400558; relerr = 0.001175237543268981
Correlation Matrix:
0 1 2
0 1.000000000 -0.948611884 0.815532447
1 -0.948611884 1.000000000 -0.646255837
2 0.815532447 -0.646255837 1.000000000
# vectors in multivariate sequence: 50000
CorrMtrx[1][0] = -0.9483147908254509; relerr = 3.8844072180899136E-4
CorrMtrx[2][0] = 0.8178271707817083; relerr = 0.0016296330995904107
CorrMtrx[2][1] = -0.6466691282330974; relerr = 0.0018155056613018417
Correlation Matrix:
0 1 2
0 1.000000000 -0.948314791 0.817827171
1 -0.948314791 1.000000000 -0.646669128
2 0.817827171 -0.646669128 1.000000000
# vectors in multivariate sequence: 500000
CorrMtrx[1][0] = -0.9486520691467761; relerr = 3.291815488037919E-5
CorrMtrx[2][0] = 0.8165491846305497; relerr = 6.442611524937192E-5
CorrMtrx[2][1] = -0.6459060179737274; relerr = 6.333003309577645E-4
Correlation Matrix:
0 1 2
0 1.000000000 -0.948652069 0.816549185
1 -0.948652069 1.000000000 -0.645906018
2 0.816549185 -0.645906018 1.000000000
Link to Java source.