// import umontreal.iro.lecuyer.rng.*; import umontreal.iro.lecuyer.randvar.NormalGen; import umontreal.iro.lecuyer.rng.RandomStream; import cern.colt.matrix.DoubleMatrix2D; import cern.colt.matrix.impl.DenseDoubleMatrix2D; import cern.colt.matrix.linalg.*; /* Outils de base pour simuler un mouvement Brownien multivarie, en c dimensions, avec d instants d'observations. Apres avoir appele le constructeur, on doit appeler init */ public class MultivariateBM { public int c; // Nombre de processus (nombre d'actifs). public int d; // Nombre d'instants d'observation. public int s; // c * d. public double[] sgi; // sgi[i] est la valeur de \sigma_i, pour le processus i. // (Vecteur de dimension c.) public double[] tj; // tj[0]=0 et tj[j] est le j-ieme instant d'observation, // pour j=1,...,d. Donc tj doit etre de dimension d+1. public double[][] rho; // rho[i][j] est la correlation entre les accroissements // des processus i et j sur un intervalle de longueur 1. // (Matrice de dimension c par c.) public double[][] sigmacd; // Matrice de covariance du vecteur // W = (W_1(t_1),...,W_c(t_1), ..., W_1(t_d),...,W_c(t_d)) // dans l'enonce du devoir. public double[][] sigmac; // Les vecteurs tj et sgi, et la matrice rho, doivent etre crees de l'exterieur. public MultivariateBM (int c, int d, double[] tj, double[] sgi, double[][] rho) { if (sgi.length != c) throw new IllegalArgumentException ("Invalid length of sgi"); if (tj.length != d + 1) throw new IllegalArgumentException ("Invalid length of tj"); if (rho.length != c || rho[0].length != c) throw new IllegalArgumentException ("Invalid dimensions of rho"); this.c = c; this.d = d; s = c * d; sigmacd = new double[s][s]; sigmac = new double[c][c]; this.tj = tj; this.sgi = sgi; this.rho = rho; } // Initialise la matrice de covariance de taille (cd)*(cd), // dont l'element [i*c+k][j*c+l] est la covariance entre les processus // k et l observee aux temps i et j. // Le vecteur tj[] contient les instants d'observations et sgi[] les // \sigma_i pour les differents processus. public void initBigSigma () { for (int i = 0; i < d; ++i) { for (int j = 0; j < d; ++j) { for (int k = 0; k < c; ++k) { for (int l = 0; l < c; ++l) { int g = (i <= j ? i + 1 : j + 1); if (k == l) sigmacd [i * c + k][j * c + l] = g * sgi[k] * sgi[l] * (tj[i + 1] - tj[i]); else sigmacd[i * c + k][j * c + l] = g * rho[k][l] * sgi[k] * sgi[l] * (tj[i + 1] - tj[i]); } } } } } // Calcule la matrice de covariance de taille c * c, // dont l'element [i][j] est la covariance entre les accroissements des // processus i et j sur un intervalle de temps de longueur deltat. // Resultat stocke dans sigmac. public void initSmallSigma (double deltat) { for (int i = 0; i < c; ++i) for (int j = 0; j < c; ++j) sigmac [i][j] = rho[i][j] * sgi[i] * sgi[j] * deltat; } // Prend en entree une matrice sigma et retourne la matrice A telle que // sigma = A A' selon la decomposition de Cholesky. public double[][] decompCholesky (double[][] sigma) { CholeskyDecomposition sv = new CholeskyDecomposition (new DenseDoubleMatrix2D (sigma)); DoubleMatrix2D L = sv.getL (); return L.toArray (); } // Prend en entree une matrice sigma et retourne la matrice A telle que // sigma = A A' selon la decomposition ACP. public double[][] decompACP (double[][] sigma) { // L'objet SingularValueDecomposition permet de recuperer la matrice // des valeurs propres en ordre decroissant et celle des vecteurs propres de // sigma (pour une matrice symetrique et definie-positive seulement). SingularValueDecomposition sv = new SingularValueDecomposition (new DenseDoubleMatrix2D (sigma)); DoubleMatrix2D D = sv.getS (); // Calculer la racine carree des valeurs propres for (int i = 0; i < D.rows(); ++i) D.setQuick (i, i, Math.sqrt (D.getQuick (i, i))); DoubleMatrix2D P = sv.getV(); return P.zMult (D, null).toArray (); } // Generates a vector of independent standard normals. public void generateZ (int dim, double z[], RandomStream stream) { for (int i = 0; i < dim; ++i) z[i] = NormalGen.nextDouble (stream, 0.0, 1.0); } // Computes A * z and returns it in vector Az. public void computeAZ (int dim, double[][] A, double z[], double Az[]) { for (int i = 0; i < dim; ++i) { double sum = 0.0; for (int k = 0; k < dim; ++k) sum += A[i][k] * z[k]; Az[i] = sum; } } }