umontreal.iro.lecuyer.stochprocess.CIRProcess 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: CIRProcess
* 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 a CIR (Cox, Ingersoll, Ross) 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
* ν = 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);
}
}