package com.imsl.test.example.stat;

import com.imsl.stat.*;
import com.imsl.math.*;

/**
 * <p>
 * Generates a pseudorandom multivariate sequence with user defined marginal
 * distributions.</p>
 * <p>
 * This example uses method <code>Random.nextGaussianCopula</code> to generate a
 * multivariate sequence
 * <code>GCdevt[k=0..nseq-1][j=0..nvar-1]</code> whose marginal distributions
 * are user-defined and imprinted with a user-specified correlation matrix
 * <code>CorrMtrxIn[i=0..nvar-1][j=0..nvar-1]</code> and then uses method
 * <code>Random.canonicalCorrelation</code> to extract from this multivariate
 * random sequence a canonical correlation matrix
 * <code>CorrMtrx[i=0..nvar-1][j=0..nvar-1]</code>.</p>
 *
<p>
 * 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 procedure is the inverse of the
 * first: an arbitrary 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.</p>
 *
<p>
 * This example demonstrates that the <code>nextGaussianCopula</code> method
 * correctly embeds 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.</p>
 *
 *
 * @see <a href="RandomEx2.java">Code</a>
 * @see <a href="doc-files/RandomEx2.out">Output</a>
 */

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:
                     */
                    switch (j) {
                        case 0:
                            //convert probs into ChiSquare(df=10) deviates:
                            GCdevt[k][j] = InvCdf.chi(GCvart[k][j], 10.);
                            break;
                        case 1:
                            //convert probs into F(dfn=15,dfd=10) deviates:
                            GCdevt[k][j] = InvCdf.F(GCvart[k][j], 15., 10.);
                            break;
                        default:
                            //convert probs into Normal(mean=0,variance=1) deviates:
                            GCdevt[k][j] = InvCdf.normal(GCvart[k][j]);
                            break;
                    }
                }
            }
            /*
             *  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);
        }
    }
}