/* *
* 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.*;
/**
* Represents a Brownian motion process
* {X(t) : t >= 0}
* sampled using the bridge sampling technique
* (see for example).
* This technique generates first the value X(td) at the last observation time,
* then the value at time td/2 (or the nearest integer),
* then the values at time td/4 and at time t3d/4
* (or the nearest integers), and so on.
* If the process has already been sampled at times ti < tk but not
* in between, the next sampling point in that interval will be
* tj where
* j = floor((i + k)/2).
* For example, if the sampling times used are
* {
* t0, t1, t2, t3, t4, t5},
* then the observations are generated in the following order:
* X(t5), X(t2), X(t1), X(t3), X(t4).
*
*
* Warning:
* Both the generatePath and the nextObservation methods from
* {@link umontreal.iro.lecuyer.stochprocess.BrownianMotion BrownianMotion} are
* modified to use the bridge method.
*
* In the case of nextObservation, the user should understand
* that the observations returned are not ordered chronologically.
* However they will be once an entire path is generated and the observations
* are read from the internal array (referenced by the getPath method)
* that contains them.
*
*
* The method nextObservation(double nextTime) differs from that of
* the class
* {@link umontreal.iro.lecuyer.stochprocess.BrownianMotion BrownianMotion}
* in that nextTime represents
* the next observation time of the Brownian bridge.
* However, the ti supplied must still be non-decreasing with i.
*
*
* Note also that, if the path is not entirely generated before being read
* from this array, there will be ``pollution'' from the previous path generated,
* and the observations will not represent a sample path of this process.
*
*/
public class BrownianMotionBridge extends BrownianMotion {
protected int bridgeCounter = -1; // Before 1st observ
// For precomputations for B Bridge
protected double[] wMuDt,
wSqrtDt;
protected int[] wIndexList,
ptIndex;
/**
* Constructs a new BrownianMotionBridge with
* parameters
* μ = mu ,
* σ = sigma and initial value
*
* X(t0) = x0 .
* The normal variates will be
* generated by inversion using the
* {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} stream.
*
*/
public BrownianMotionBridge (double x0, double mu, double sigma,
RandomStream stream) {
super (x0, mu, sigma, stream);
}
/**
* Constructs a new BrownianMotionBridge with
* parameters
* μ = mu ,
* σ = sigma and initial value
*
* X(t0) = x0 .
* The normal variates will be
* generated by the
* {@link umontreal.iro.lecuyer.randvar.NormalGen NormalGen} gen.
*
*
*/
public BrownianMotionBridge (double x0, double mu, double sigma,
NormalGen gen) {
super (x0, mu, sigma, gen);
}
public double nextObservation() {
double x;
if (bridgeCounter == -1) {
x = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();
bridgeCounter = 0;
observationIndex = d;
} else {
int j = bridgeCounter*3;
int oldIndexL = wIndexList[j];
int newIndex = wIndexList[j + 1];
int oldIndexR = wIndexList[j + 2];
x = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL])
* wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();
bridgeCounter++;
observationIndex = newIndex;
}
observationCounter = bridgeCounter + 1;
path[observationIndex] = x;
return x;
}
public double nextObservation (double nextTime) {
double x;
if (bridgeCounter == -1) {
t[d] = nextTime;
wMuDt[0] = 0.0; // The end point of the Wiener process
// w/ Brownian bridge has expectation = 0
wSqrtDt[0] = sigma * Math.sqrt(t[d] - t[0]);
// = sigma*sqrt(Dt) of end point
x = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();
bridgeCounter = 0;
observationIndex = d;
} else {
int j = bridgeCounter*3;
int oldIndexL = wIndexList[j];
int newIndex = wIndexList[j + 1];
int oldIndexR = wIndexList[j + 2];
t[newIndex] = nextTime;
double dtRL = t[oldIndexR] - t[oldIndexL];
if (dtRL != 0.0) {
wMuDt[newIndex] = (t[newIndex]-t[oldIndexL]) / dtRL;
} else {
wMuDt[newIndex] = 0.0;
}
wSqrtDt[newIndex] = sigma * Math.sqrt (
wMuDt[newIndex] * (t[oldIndexR] - t[newIndex]));
x = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL])
* wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();
bridgeCounter++;
observationIndex = newIndex;
}
observationCounter = bridgeCounter + 1;
path[observationIndex] = x;
return x;
}
public double[] generatePath() {
// Generation of Brownian bridge process
int oldIndexL, oldIndexR, newIndex;
path[d] = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = wIndexList[j];
newIndex = wIndexList[j + 1];
oldIndexR = wIndexList[j + 2];
path[newIndex] = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL])
* wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();
}
// resetStartProcess();
observationIndex = d;
observationCounter = d;
return path;
}
public double[] generatePath (double[] uniform01){
int oldIndexL, oldIndexR, newIndex;
path[d] = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * NormalDist.inverseF01(uniform01[0]);
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = wIndexList[j];
newIndex = wIndexList[j + 1];
oldIndexR = wIndexList[j + 2];
path[newIndex] = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL])
* wMuDt[newIndex] + wSqrtDt[newIndex] * NormalDist.inverseF01(uniform01[1 + j/3]);
}
// resetStartProcess();
observationIndex = d;
observationCounter = d;
return path;
}
public void resetStartProcess() {
observationIndex = 0;
observationCounter = 0;
bridgeCounter = -1;
}
protected void init() {
super.init();
/* For Brownian Bridge */
// Quantities for Brownian Bridge process
wMuDt = new double[d + 1];
wSqrtDt = new double[d + 1];
wIndexList = new int[3 * (d)];
ptIndex = new int[d + 1];
double tem = 0;
int indexCounter = 0;
int newIndex, oldLeft, oldRight;
ptIndex[0] = 0;
ptIndex[1] = d;
wMuDt[0] = 0.0; // The end point of the Wiener process
// w/ Brownian bridge has expectation = 0
if (t[d] < t[0])
throw new IllegalStateException(" t[d] < t[0]");
wSqrtDt[0] = sigma * Math.sqrt(t[d] - t[0]);
// = sigma*sqrt(Dt) of end point
for (int powOfTwo = 1; powOfTwo <= d / 2; powOfTwo *= 2) {
/* Make room in the indexing array "ptIndex" */
for (int j = powOfTwo; j >= 1; j--) {
ptIndex[2*j] = ptIndex[j];
}
/* Insert new indices and Calculate constants */
for (int j = 1; j <= powOfTwo; j++) {
oldLeft = 2 * j - 2;
oldRight = 2 * j;
newIndex = (int) (0.5 * (ptIndex[oldLeft] + ptIndex[oldRight]));
wMuDt[newIndex] = (t[newIndex] - t[ptIndex[oldLeft]]) /
(t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);
tem = (t[newIndex] - t[ptIndex[oldLeft]]) *
(t[ptIndex[oldRight]] - t[newIndex])
/ (t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);
// Test for NaN (z != z); 0/0 gives a NaN
if (tem < 0 || tem != tem) {
System.out.printf ("t[newIndex] - t[ptIndex[oldLeft]] = %g%n", t[newIndex] - t[ptIndex[oldLeft]]);
System.out.printf ("t[ptIndex[oldRight]] - t[newIndex] = %g%n", t[ptIndex[oldRight]] - t[newIndex]);
System.out.printf ("t[ptIndex[oldRight]] - t[ptIndex[oldLeft]] = %g%n", t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);
System.out.printf ("t[ptIndex[oldRight]] = %g%n", t[ptIndex[oldRight]]);
System.out.printf ("t[ptIndex[oldLeft]] = %g%n", t[ptIndex[oldLeft]]);
throw new IllegalStateException(" tem < 0 or NaN");
}
wSqrtDt[newIndex] = sigma * Math.sqrt (tem);
ptIndex[oldLeft + 1] = newIndex;
wIndexList[indexCounter] = ptIndex[oldLeft];
wIndexList[indexCounter + 1] = newIndex;
wIndexList[indexCounter + 2] = ptIndex[oldRight];
indexCounter += 3;
}
}
/* Check if there are holes remaining and fill them */
for (int k = 1; k < d; k++) {
if (ptIndex[k - 1] + 1 < ptIndex[k]) {
// there is a hole between (k-1) and k.
wMuDt[ptIndex[k - 1] + 1] = (t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]) /
(t[ptIndex[k]] - t[ptIndex[k - 1]]);
tem = (t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]) *
(t[ptIndex[k]] - t[ptIndex[k - 1] + 1])
/ (t[ptIndex[k]] - t[ptIndex[k - 1]]);
// Test for NaN (z != z); 0/0 gives a NaN
if (tem < 0 || tem != tem) {
System.out.printf ("t[ptIndex[k-1]+1] - t[ptIndex[k-1]] = %g%n", t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]);
System.out.printf ("t[ptIndex[k]] - t[ptIndex[k-1]+1] = %g%n", t[ptIndex[k]] - t[ptIndex[k - 1] + 1]);
System.out.printf ("t[ptIndex[k]] - t[ptIndex[k-1]] = %g%n", t[ptIndex[k]] - t[ptIndex[k - 1]]);
System.out.printf ("t[ptIndex[k]] = %20.16g%n", t[ptIndex[k]]);
System.out.printf ("t[ptIndex[k-1]] = %20.16g%n", t[ptIndex[k - 1]]);
throw new IllegalStateException(" tem < 0 or NaN");
}
wSqrtDt[ptIndex[k - 1] + 1] = sigma * Math.sqrt (tem);
wIndexList[indexCounter] = ptIndex[k] - 2;
wIndexList[indexCounter + 1] = ptIndex[k] - 1;
wIndexList[indexCounter + 2] = ptIndex[k];
indexCounter += 3;
}
}
}
}
Ads help maintain this website.