GammaProcess.java

/*
 * Class:        GammaProcess
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre Tremblay, Jean-Sébastien Parent & Maxime Dion
 * @since        July 2003
 *
 *
 * 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.randvar.*;
import umontreal.ssj.probdist.GammaDist;
import umontreal.ssj.util.Num;

public class GammaProcess extends StochasticProcess {
   // Make sure that process is strictly increasing: when two successive steps
   // are equal, reset the new one to old*(1 + EPS)
   protected static final double EPS = 1.0e-15;
   protected boolean usesAnti = false;
 
   protected RandomStream stream; // used in the generatePath(), so must be kept
   protected GammaGen Ggen;
   protected double mu, nu;
   protected double muOverNu, mu2OverNu;
 
   // For precomputations for standard GP
   protected double[] mu2dtOverNu;
 
   protected void setLarger(double[] path, int left, int mid, int right) {
      /*
       * make sure that path[left] < path[mid] < path[right] by adding a small epsilon
       * to the last two
       */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      path[mid] = path[left] * (1.0 + myEps);
      if (path[mid] >= path[right])
         path[right] = path[mid] * (1.0 + myEps);
   }
 
   protected double setLarger(double[] path, int left, int right) {
      /*
       * make sure that path[left] < s < path[right] by adding a small epsilon to the
       * last two; return s.
       */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      double s = path[left] * (1.0 + myEps);
      if (s >= path[right])
         path[right] = s * (1.0 + myEps);
      return s;
   }
 
   protected double setLarger(double v) {
      /* return a number that is slightly larger than v */
      if (v >= 0)
         return v * (1.0 + EPS) + Num.DBL_EPSILON;
      else
         return v * (1.0 - EPS);
   }

   public GammaProcess(double s0, double mu, double nu, RandomStream stream) {
      this(s0, mu, nu, new GammaGen(stream, 1.0));
// Note that the parameters (the 1.0) supplied to
// the GammaGen constructors are meaningless.
// The 'real' parameters are supplied at time of generation
   }

   public GammaProcess(double s0, double mu, double nu, GammaGen Ggen) {
      this.mu = mu;
      this.nu = nu;
      this.x0 = s0;
      this.Ggen = Ggen;
      stream = Ggen.getStream();
   }
 
   public double nextObservation() {
      double s = path[observationIndex];
      double v = s;
      s += Ggen.nextDouble(stream, mu2dtOverNu[observationIndex], muOverNu);
      if (s <= v)
         s = setLarger(v);
      observationIndex++;
      observationCounter = observationIndex;
      path[observationIndex] = s;
      return s;
   }

   public double nextObservation(double nextT) {
      // This method is useful for generating variance gamma processes
      double s = path[observationIndex];
      double v = s;
      double previousT = t[observationIndex];
      observationIndex++;
      observationCounter = observationIndex;
      t[observationIndex] = nextT;
      double dt = nextT - previousT;
      s += Ggen.nextDouble(stream, mu2OverNu * dt, muOverNu);
      if (s <= v)
         s = setLarger(v);
      path[observationIndex] = s;
      return s;
   }

   public double[] generatePath() {
      double s = x0;
      double v;
      for (int i = 0; i < d; i++) {
         v = s;
         s += Ggen.nextDouble(stream, mu2dtOverNu[i], muOverNu);
         if (s <= v)
            s = setLarger(v);
         path[i + 1] = s;
      }
      observationIndex = d;
      observationCounter = d;
      return path;
   }

   public double[] generatePath(double[] uniform01) {
      double s = x0;
      double v;
      for (int i = 0; i < d; i++) {
         v = s;
         s += GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniform01[i]);
         if (s <= v)
            s = setLarger(v);
         path[i + 1] = s;
      }
      observationIndex = d;
      observationCounter = d;
      return path;
   }

   public void setParams(double s0, double mu, double nu) {
      this.x0 = s0;
      this.mu = mu;
      this.nu = nu;
      if (observationTimesSet)
         init(); // Otherwise no need to.
   }

   public double getMu() {
      return mu;
   }

   public double getNu() {
      return nu;
   }

   public void setStream(RandomStream stream) {
      this.stream = stream;
      Ggen.setStream(stream);
   }

   public RandomStream getStream() {
      return stream;
   }
 
   protected void init() {
      // init of StochasticProcess only sets path[0] = x0 if observation times are
      // set.
      super.init();
      muOverNu = mu / nu;
      mu2OverNu = mu * mu / nu;
      mu2dtOverNu = new double[d];
      if (observationTimesSet) {
         for (int i = 0; i < d; i++)
            mu2dtOverNu[i] = mu2OverNu * (t[i + 1] - t[i]);
      }
   }
 
}