MultivariateBrownianMotion.java

/*
 * Class:        MultivariateBrownianMotion
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre L'Écuyer, Jean-Sébastien Parent & Clément Teule
 * @since        2008
 
 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.
 
 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 
 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stochprocess;
 
import umontreal.ssj.rng.*;
import umontreal.ssj.util.DMatrix;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.CholeskyDecomposition;

public class MultivariateBrownianMotion extends MultivariateStochasticProcess {
   protected NormalGen gen;
   protected double[] mu;
   protected double[] sigma;
   protected double[][] corrZ; // Correlation matrix for Z.
   protected DoubleMatrix2D covZ; // Covariance Matrix
   protected DoubleMatrix2D covZCholDecomp; // Matrix obtained using the cholesky decomposition on covZ
   protected CholeskyDecomposition decomp;
   protected boolean covZiSCholDecomp;
 
   // Precomputed values for standard BM
   protected double[] dt, sqrdt;
 
   protected MultivariateBrownianMotion() {
   };

   public MultivariateBrownianMotion(int c, double[] x0, double[] mu, double[] sigma, double[][] corrZ,
         RandomStream stream) {
      setParams(c, x0, mu, sigma, corrZ);
      this.gen = new NormalGen(stream, new NormalDist());
   }

   public MultivariateBrownianMotion(int c, double[] x0, double[] mu, double[] sigma, double[][] corrZ, NormalGen gen) {
      setParams(c, x0, mu, sigma, corrZ);
      this.gen = gen;
   }

   public void nextObservationVector(double[] obs) {
      if (!covZiSCholDecomp) 
         // the cholesky factorisation must be done to use the matrix covZCholDecomp
         initCovZCholDecomp();
      double z;
      for (int i = 0; i < c; i++) {
         z = 0.0;
         for (int k = 0; k < c; k++)
            z += covZCholDecomp.getQuick(i, k) * gen.nextDouble();
         obs[i] = path[c * observationIndex + i] + mu[i] * dt[observationIndex] + sqrdt[observationIndex] * z;
         path[c * (observationIndex + 1) + i] = obs[i];
      }
      observationIndex++;
   }

   public double[] nextObservationVector() {
      double[] obs = new double[c];
      nextObservationVector(obs);
      return obs;
   }

   public double[] nextObservationVector(double nextTime, double[] obs) {
      if (!covZiSCholDecomp) // the cholesky factorisation must be done to use the matrix covZCholDecomp
         initCovZCholDecomp();
      t[observationIndex + 1] = nextTime;
      double z;
      for (int i = 0; i < c; i++) {
         z = 0.0;
         for (int k = 0; k < c; k++)
            z += covZCholDecomp.getQuick(i, k) * gen.nextDouble();
         obs[i] = path[c * observationIndex + i] + mu[i] * (t[observationIndex + 1] - t[observationIndex])
               + Math.sqrt(t[observationIndex + 1] - t[observationIndex]) * z;
         path[c * (observationIndex + 1) + i] = obs[i];
      }
      observationIndex++;
      return obs;
   }

   public double[] nextObservationVector(double x[], double dt) {
      double[] obs = new double[c];
      double z;
      if (!covZiSCholDecomp) // the cholesky factorisation must be done to use the matrix covZCholDecomp
         initCovZCholDecomp();
 
      for (int i = 0; i < c; i++) {
         z = 0.0;
         for (int k = 0; k < c; k++)
            z += covZCholDecomp.getQuick(i, k) * gen.nextDouble();
         obs[i] = x[i] + mu[i] * dt + Math.sqrt(dt) * z;
      }
      observationIndex++;
      return obs;
   }
 
   public double[] generatePath() {
      double[] u = new double[c * d];
      for (int i = 0; i < c * d; i++)
         u[i] = gen.nextDouble();
      return generatePath(u);
   }

   public double[] generatePath(double[] uniform01) {
      if (!covZiSCholDecomp) { 
         // the cholesky factorisation must be done to use the matrix covZCholDecomp
         initCovZCholDecomp();
      }
      double z;
      double[] QMCpointsPart = new double[c];
      for (int i = 0; i < c; i++)
         path[i] = x0[i];
      for (int j = 0; j < d; j++) {
         for (int i = 0; i < c; i++)
            QMCpointsPart[i] = uniform01[j * c + i];
         for (int i = 0; i < c; i++) {
            z = 0.0;
            for (int k = 0; k < c; k++)
               z += covZCholDecomp.getQuick(i, k) * QMCpointsPart[k];
            path[c * (j + 1) + i] = path[c * j + i] + mu[i] * dt[j] + sqrdt[j] * z;
         }
      }
      observationIndex = d;
      return path;
   }
 
   public double[] generatePath(RandomStream stream) {
      gen.setStream(stream);
      return generatePath();
   }

   public void setParams(int c, double x0[], double mu[], double sigma[], double corrZ[][]) {
      this.c = c;
      this.x0 = x0;
      this.mu = mu;
      this.sigma = sigma;
      this.corrZ = corrZ;
 
      if (x0.length < c)
         throw new IllegalArgumentException(
               "x0 dimension :  " + x0.length + " is smaller than the process dimension : " + c);
      if (mu.length < c)
         throw new IllegalArgumentException(
               "mu dimension :  " + mu.length + " is smaller than the process dimension : " + c);
      if (sigma.length < c)
         throw new IllegalArgumentException(
               "sigma dimension :  " + sigma.length + " is smaller than the process dimension : " + c);
      if ((corrZ.length < c) || (corrZ[0].length < c))
         throw new IllegalArgumentException("corrZ dimensions :  " + corrZ.length + "x" + corrZ[0].length
               + " are smaller than the process dimension : " + c);
      covZ = new DenseDoubleMatrix2D(c, c);
      initCovZ();
      covZiSCholDecomp = false;
 
      if (observationTimesSet)
         init(); // Otherwise not needed.
   }

   public void setParams(double x0[], double mu[], double sigma[]) {
      this.x0 = x0;
      this.mu = mu;
      this.sigma = sigma;
      covZ = new DenseDoubleMatrix2D(c, c);
      initCovZ();
      covZiSCholDecomp = false;
 
      if (observationTimesSet)
         init(); // Otherwise not needed.
   }

   public void setStream(RandomStream stream) {
      gen.setStream(stream);
   }

   public RandomStream getStream() {
      return gen.getStream();
   }

   public NormalGen getGen() {
      return gen;
   }
 
   // This is called by setObservationTimes to precompute constants
   // in order to speed up the path generation.
   protected void init() {
      super.init();
      dt = new double[d + 1];
      sqrdt = new double[d + 1];
      for (int j = 0; j < d; j++) {
         dt[j] = t[j + 1] - t[j];
         sqrdt[j] = Math.sqrt(dt[j]);
      }
   }
 
   protected void initCovZCholDecomp() {
      covZCholDecomp = new DenseDoubleMatrix2D(c, c);
      covZCholDecomp = DMatrix.CholeskyDecompose(covZ);
      covZiSCholDecomp = true;
   }
 
   protected void initCovZ() {
      for (int i = 0; i < c; i++)
         for (int j = 0; j < c; j++)
            covZ.setQuick(i, j, corrZ[i][j] * sigma[i] * sigma[j]);
   }

   public double[] getMu() {
      return mu;
   }
 
}