BrownianMotionPCA.java

/*
 * Class:        BrownianMotionPCA
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.stochprocess;
 
import umontreal.ssj.rng.*;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.*;

public class BrownianMotionPCA extends BrownianMotion {
 
   protected double[][] sigmaCov; // Matrice de covariance du vecteur des observ.
   // sigmaCov [i][j] = Cov[X(t_{i+1}),X(t_{j+1})].
   protected double[][] A; // sigmaCov = AA' (PCA decomposition).
   protected double[] z; // vector of standard normals.
   protected double[] sortedEigenvalues;
   protected boolean isDecompPCA;

   public BrownianMotionPCA(double x0, double mu, double sigma, RandomStream stream) {
      super(x0, mu, sigma, stream);
      isDecompPCA = false;
   }

   public BrownianMotionPCA(double x0, double mu, double sigma, NormalGen gen) {
      super(x0, mu, sigma, gen);
      isDecompPCA = false;
   }
 
   public double nextObservation() {
      throw new UnsupportedOperationException("nextObservation() is not defined in BrownianMotionPCA");
   }
 
   public void setParams(double x0, double mu, double sigma) {
      super.setParams(x0, mu, sigma);
      isDecompPCA = true; // the setParams method defined in BrownianMotion calls init() which does
                          // the PCA decompostion
   }
 
   public double[] generatePath() {
      if (!isDecompPCA) {
         init();
      } // if the decomposition is not done, do it...
      for (int j = 0; j < d; j++)
         z[j] = gen.nextDouble();
      for (int j = 0; j < d; j++) {
         double sum = 0.0;
         for (int k = 0; k < d; k++)
            sum += A[j][k] * z[k];
         path[j + 1] = x0 + mu * t[j + 1] + sum;
      }
      observationIndex = d;
      observationCounter = d;
      return path;
   }
 
   public double[] generatePath(double[] uniform01) {
      if (!isDecompPCA) {
         init();
      } // if the decomposition is not done, do it...
      for (int j = 0; j < d; j++)
         z[j] = NormalDist.inverseF01(uniform01[j]);
      for (int j = 0; j < d; j++) {
         double sum = 0.0;
         for (int k = 0; k < d; k++)
            sum += A[j][k] * z[k];
         path[j + 1] = x0 + mu * t[j + 1] + sum;
      }
      observationIndex = d;
      observationCounter = d;
      return path;
   }
 
   public double[][] decompPCA(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(); // diagonal
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); i++) {
         sortedEigenvalues[i] = D.getQuick(i, i);
         D.setQuick(i, i, Math.sqrt(D.getQuick(i, i)));
      }
      DoubleMatrix2D P = sv.getV(); // right factor matrix
      return P.zMult(D, null).toArray();
   }

   public double[] getSortedEigenvalues() {
      return sortedEigenvalues;
   }
 
   protected void init() {
      super.init();
      z = new double[d];
      sortedEigenvalues = new double[d];
 
      // Initialize sigmaCov, based on the observation times.
      sigmaCov = new double[d][d];
      for (int i = 0; i < d; i++) {
         for (int j = i; j < d; j++) {
            sigmaCov[i][j] = sigma * sigma * t[i + 1];
            sigmaCov[j][i] = sigmaCov[i][j];
         }
      }
      A = decompPCA(sigmaCov);
      isDecompPCA = true;
   }
}