/* *
 * 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 (Cox, Ingersoll, Ross) process
 *  
 * {X(t) : t >= 0}, sampled at times 
 * 0 = t0 < t1 <  ...  < td.
 * This process obeys the stochastic differential equation
 * 

 * 
 * 

 * dX(t) = α(b - X(t))dt + σ(X(t))1/2 dB(t)
 * 

 * with initial condition X(0) = x0,
 * where α, b and σ are positive constants,
 * and 
 * {B(t), t >= 0} is a standard Brownian motion
 * (with drift 0 and volatility 1).
 * This process is mean-reverting in the sense that it always tends to
 * drift toward its general mean b.
 * The process is generated using the sequential
 * technique
 * 

 * 
 * 

 * X(tj) = {σ2(1-e-α(tj-tj-1))/4α}χ′2ν({4αe-α(tj-tj-1)X(tj-1)}/{σ2(1-e-α(tj-tj-1))}),
 * 

 * where 
 * ν = 4bα/σ2, and 
 * χ′ 2ν(λ) is a noncentral
 * chi-square random variable with ν degrees of freedom and noncentrality
 * parameter λ.
 * 
 */
public class CIRProcess extends StochasticProcess  {
    private RandomStream stream;
    protected ChiSquareNoncentralGen gen;
    protected double alpha,
                     beta,
                     sigma,
                     nu;     // number of degrees of freedom
    // Precomputed values
    protected double[] parc,
                       parlam;



   /**
    * Constructs a new CIRProcess with parameters
    * 
    * α = alpha, b, 
    * σ = sigma and initial value
    * 
    * X(t0) = x0. The noncentral chi-square variates
    *   
    * χ′2ν(λ) will be
    * generated by inversion using the stream stream.
    * 
    */
   public CIRProcess (double x0, double alpha, double b, double sigma,
                      RandomStream stream)  {
      this (x0, alpha, b, sigma, new ChiSquareNoncentralGen (stream, null));
    }


   /**
    * The noncentral chi-square variate generator  gen
    *  is specified directly instead of specifying the stream.
    *  gen can use a method other than inversion.
    * 
    */
   public CIRProcess (double x0, double alpha, double b, double sigma,
                      ChiSquareNoncentralGen gen)  {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      nu = 4.0*b*alpha/(sigma*sigma);
      this.gen   = gen;
      stream = gen.getStream();
    }


   public double nextObservation() {
      double xOld = path[observationIndex];
      double lambda = xOld * parlam[observationIndex];
      double x;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = parc[observationIndex] *
             ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      observationIndex++;
      path[observationIndex] = x;
      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 = nextObservation (xOld, dt);
      path[observationIndex] = x;
      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)  {
      double c = -Math.expm1(-alpha * dt) * sigma * sigma / (4.0*alpha);
      double lambda = x * Math.exp(-alpha * dt) / c;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = c * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = c * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = c * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      return x;
    }


   public double[] generatePath() {
      double xOld = x0;
      double x, lambda;
      int j;

      if (gen.getClass() == ChiSquareNoncentralPoisGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else if (gen.getClass() == ChiSquareNoncentralGamGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
            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;
      nu = 4.0*b*alpha/(sigma*sigma);
      if (observationTimesSet) init(); // Otherwise not needed.
    }


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


   /**
    * Returns the random stream of the noncentral chi-square 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 noncentral chi-square random variate generator used.
    * The RandomStream used for that generator can be changed via
    * getGen().setStream(stream), for example.
    * 
    */
   public ChiSquareNoncentralGen getGen()  { return gen; }
 

   protected void initArrays(int d) {
      double dt, c;
      for (int j = 0; j < d; j++) {
         dt = t[j+1] - t[j];
         c = -Math.expm1(-alpha * dt)*sigma*sigma/(4.0*alpha);
         parc[j] = c;
         parlam[j] = Math.exp(-alpha * dt) / c;
      }
   }

   // This is called by setObservationTimes to precompute constants
   // in order to speed up the path generation.
   protected void init() {
       super.init();
       parc = new double[d];
       parlam = new double[d];
       initArrays(d);
    }

}