Example 2: Using Copulas to imprint and extract correlation information

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);
        }
    }
}

Output

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.