/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */


package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * .
 * 
 * This class represents a CIR process
 * as in {@link CIRProcess}, but
 * the process is generated using the simple Euler scheme
 * 

 * 
 * 

 * X(tj) - X(tj-1) = α(b - X(tj-1))(tj - tj-1) + σ((t_j - t_j-1)X(t_j-1))1/2 Zj
 * 

 * where 
 * Zj∼N(0, 1). This is a good approximation only for small
 * time intervals 
 * tj - tj-1.
 * 
 */
public class CIRProcessEuler extends StochasticProcess  {
    protected NormalGen    gen;
    protected double       alpha,
                           beta,
                           sigma;
    // Precomputed values
    protected double[]     alphadt,
                           sigmasqrdt;



   /**
    * Constructs a new CIRProcessEuler with parameters
    * α = alpha, b, σ = sigma and initial value
    * 
    * X(t0) = x0. The normal variates Zj will be
    * generated by inversion using the stream stream.
    * 
    */
   public CIRProcessEuler (double x0, double alpha, double b, double sigma,
                           RandomStream stream)  {
      this (x0, alpha, b, sigma, new NormalGen (stream, new NormalDist()));
   }


   /**
    * The normal variate generator gen is specified directly
    * instead of specifying the stream.
    *  gen can use another method than inversion.
    * 
    */
   public CIRProcessEuler (double x0, double alpha, double b, double sigma,
                           NormalGen gen)  {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      this.gen   = gen;
   }


   public double nextObservation() {
      double xOld = path[observationIndex];
      double x;
      x = xOld + (beta - xOld) * alphadt[observationIndex]
           + sigmasqrdt[observationIndex] * Math.sqrt(xOld) * gen.nextDouble();
      observationIndex++;
      if (x >= 0.0)
         path[observationIndex] = x;
      else
         path[observationIndex] = 0.;
      return x;
   }


   /**
    * Generates and returns the next observation at time tj+1 =
    *  nextTime, using the previous observation time tj defined earlier
    * (either by this method or by setObservationTimes),
    * as well as the value of the previous observation X(tj).
    * Warning: This method will reset the observations time tj+1
    * for this process to nextTime. The user must make sure that
    * the tj+1 supplied is 
    *  >= tj.
    * 
    */
   public double nextObservation (double nextTime)  {
      double previousTime = t[observationIndex];
      double xOld = path[observationIndex];
      observationIndex++;
      t[observationIndex] = nextTime;
      double dt = nextTime - previousTime;
      double x = xOld + alpha * (beta - xOld) * dt
           + sigma * Math.sqrt (dt*xOld) * gen.nextDouble();
      if (x >= 0.0)
         path[observationIndex] = x;
      else
         path[observationIndex] = 0.;
      return x;
   }


   /**
    * Generates an observation of the process in dt time units,
    * assuming that the process has value x at the current time.
    * Uses the process parameters specified in the constructor.
    * Note that this method does not affect the sample path of the process
    * stored internally (if any).
    * 
    * 
    */
   public double nextObservation (double x, double dt)  {
      x = x + alpha * (beta - x) * dt +
          sigma * Math.sqrt (dt*x) * gen.nextDouble();
      if (x >= 0.0)
         return x;
      return 0.0;
   }


   public double[] generatePath() {
      double x;
      double xOld = x0;
      for (int j = 0; j < d; j++) {
          x = xOld + (beta - xOld) * alphadt[j]
              + sigmasqrdt[j] * Math.sqrt(xOld) * gen.nextDouble();
          if (x < 0.0)
              x = 0.0;
          path[j + 1] = x;
          xOld = x;
      }
      observationIndex = d;
      return path;
   }

   public double[] generatePath (RandomStream stream) {
        gen.setStream (stream);
        return generatePath();
   }


   /**
    * Resets the parameters 
    * X(t0) = x0, α = alpha,
    *  b = b and σ = sigma of the process.
    * Warning: This method will recompute some quantities stored internally,
    * which may be slow if called too frequently.
    * 
    */
   public void setParams (double x0, double alpha, double b, double sigma)  {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      if (observationTimesSet) init(); // Otherwise not needed.
   }


   /**
    * Resets the random stream of the normal generator to stream.
    * 
    */
   public void setStream (RandomStream stream)  {
      gen.setStream (stream);
   }


   /**
    * Returns the random stream of the normal generator.
    * 
    */
   public RandomStream getStream()  {
      return gen.getStream ();
   }


   /**
    * Returns the value of α.
    * 
    */
   public double getAlpha()  { return alpha; }


   /**
    * Returns the value of b.
    * 
    */
   public double getB()  { return beta; }


   /**
    * Returns the value of σ.
    * 
    */
   public double getSigma()  { return sigma; }


   /**
    * Returns the normal random variate generator used.
    * The RandomStream used for that generator can be changed via
    * getGen().setStream(stream), for example.
    * 
    */
   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();
       alphadt = new double[d];
       sigmasqrdt = new double[d];
       double dt;
       for (int j = 0; j < d; j++) {
           dt = t[j+1] - t[j];
           alphadt[j]      = alpha * (dt);
           sigmasqrdt[j]   = sigma * Math.sqrt (dt);
       }
    }

}