/* *
 * 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 an Ornstein-Uhlenbeck 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 + σ 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) = e-α(tj-tj-1)X(tj-1) + b(1 - e-α(tj-tj-1)) + σ(1 - e^-2α(t_j - t_j-1)2α)1/2 Zj
 * 

 * where 
 * Zj∼N(0, 1). The time intervals 
 * tj - tj-1
 * can be arbitrarily large.
 * 
 */
public class OrnsteinUhlenbeckProcess extends StochasticProcess  {
    protected NormalGen    gen;
    protected double       alpha,
                           beta,
                           sigma;
    // Precomputed values 
    protected double[]     badt,
                           alphadt,
                           sigmasqrdt;



   /**
    * Constructs a new OrnsteinUhlenbeckProcess 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 OrnsteinUhlenbeckProcess (double x0, double alpha, double b,
                                    double sigma, RandomStream stream) {
        this (x0, alpha, b, sigma, new NormalGen (stream));
    }


   /**
    * Here, the normal variate generator is specified directly
    * instead of specifying the stream.
    * The normal generator gen can use another method than inversion.
    * 
    */
   public OrnsteinUhlenbeckProcess (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 = badt[observationIndex] + xOld * alphadt[observationIndex]
                   + sigmasqrdt[observationIndex] * gen.nextDouble();
        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 tem = Math.exp(-alpha * dt);
        double tem1 = -Math.expm1(-alpha * dt);
        double x = tem*xOld + beta*tem1 + sigma *
            Math.sqrt(tem1*(1.0 + tem)/(2.0*alpha)) * gen.nextDouble();
        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 tem = Math.exp(-alpha * dt);
        double tem1 = -Math.expm1(-alpha * dt);
        x = tem*x + beta*tem1 + sigma *
            Math.sqrt(tem1*(1.0 + tem)/(2.0*alpha)) * gen.nextDouble();
        return x;
    }


   public double[] generatePath() {
        double x;
        double xOld = x0;
        for (int j = 0; j < d; j++) {
            x = badt[j] + xOld * alphadt[j] + sigmasqrdt[j] * gen.nextDouble();
            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; }
 

   protected void initArrays(int d) {
      double dt, tem, tem1;
      for (int j = 0; j < d; j++) {
         dt = t[j+1] - t[j];
         tem = Math.exp(-alpha * dt);
         tem1 = -Math.expm1(-alpha * dt);
         badt[j] = beta*tem1;
         alphadt[j] = tem;
         sigmasqrdt[j] = sigma * Math.sqrt (tem1*(1.0 + tem)/(2.0*alpha));
      }
   }

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

}