umontreal.iro.lecuyer.stochprocess.OrnsteinUhlenbeckProcess Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: OrnsteinUhlenbeckProcess
* Description:
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* 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
GPL licence site.
*/
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
*
*
*
* 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
*
*
*
* 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);
}
}