umontreal.ssj.stochprocess.GammaProcessPCABridge 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
Stochastic Simulation in Java
/*
* Class: GammaProcessPCABridge
* Description:
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Universite de Montreal
* Organization: DIRO, Universite de Montreal
* @author Jean-Sebastien Parent and Maxime Dion
* @since july 2008
*
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*/
package umontreal.ssj.stochprocess;
import umontreal.ssj.rng.*;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;
/**
* Same as @ref GammaProcessPCA, but the generated uniforms correspond to a
* bridge transformation of the @ref BrownianMotionPCA instead of a
* sequential transformation.
*
*
*/
public class GammaProcessPCABridge extends GammaProcessPCA {
protected BrownianMotionBridge BMBridge;
protected double mu2OverNu,
mu2dTOverNu;
protected double[] bMu2dtOverNuL, // For precomputations for G Bridge
bMu2dtOverNuR;
protected int[] wIndexList;
/**
* Constructs a new `GammaProcessPCABridge` with parameters @f$\mu=
* \mathtt{mu}@f$, @f$\nu= \mathtt{nu}@f$ and initial value @f$S(t_0)
* = \mathtt{s0}@f$. The same stream is used to generate the gamma and
* beta random numbers. All these numbers are generated by inversion in
* the following order: the first uniform random number generated is
* used for the gamma and the other @f$d-1@f$ for the beta’s.
*/
public GammaProcessPCABridge (double s0, double mu, double nu,
RandomStream stream) {
super (s0, mu, nu, stream);
this.BMBridge = new BrownianMotionBridge(0.0, 0.0, Math.sqrt(nu), stream);
}
public double[] generatePath (double[] uniform01) {
// uniformsV[] of size d+1, but element 0 never used.
double[] uniformsV = new double[d+1];
// generate BrownianMotion PCA path
double[] BMPCApath = BMPCA.generatePath(uniform01);
int oldIndexL;
int newIndex;
int oldIndexR;
double temp, y;
// Transform BMPCA points to uniforms using an inverse bridge.
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = BMBridge.wIndexList[j];
newIndex = BMBridge.wIndexList[j + 1];
oldIndexR = BMBridge.wIndexList[j + 2];
temp = BMPCApath[newIndex] - BMPCApath[oldIndexL];
temp -= (BMPCApath[oldIndexR] - BMPCApath[oldIndexL]) * BMBridge.wMuDt[newIndex];
temp /= BMBridge.wSqrtDt[newIndex];
uniformsV[newIndex] = NormalDist.cdf01( temp );
}
double dT = BMPCA.t[d] - BMPCA.t[0];
uniformsV[d] = NormalDist.cdf01( ( BMPCApath[d] - BMPCApath[0] - BMPCA.mu*dT )/
( BMPCA.sigma * Math.sqrt(dT) ) );
// Reconstruct the bridge for the GammaProcess from the Brownian uniforms.
path[0] = x0;
path[d] = x0 + GammaDist.inverseF(mu2dTOverNu, muOverNu, 10, uniformsV[d]);
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = wIndexList[j];
newIndex = wIndexList[j + 1];
oldIndexR = wIndexList[j + 2];
y = BetaDist.inverseF(bMu2dtOverNuL[newIndex], bMu2dtOverNuR[newIndex], 8, uniformsV[newIndex]);
path[newIndex] = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL]) * y;
}
observationIndex = d;
observationCounter = d;
return path;
}
public double[] generatePath() {
double[] u = new double[d];
for(int i =0; i < d; i++)
u[i] = stream.nextDouble();
return generatePath(u);
}
public void setParams (double s0, double mu, double nu) {
super.setParams(s0, mu, nu);
BMBridge.setParams(0.0, 0.0, Math.sqrt(nu));
BMPCA.setParams(0.0, 0.0, Math.sqrt(nu));
}
public void setObservationTimes (double[] t, int d) {
super.setObservationTimes(t, d);
BMBridge.setObservationTimes(t, d);
}
/**
* Returns the inner @ref BrownianMotionPCA.
*/
public BrownianMotionPCA getBMPCA() {
return BMPCA;
}
protected void init() {
super.init();
if (observationTimesSet) {
// Quantities for gamma bridge process
bMu2dtOverNuL = new double[d+1];
bMu2dtOverNuR = new double[d+1];
wIndexList = new int[3*d];
int[] ptIndex = new int[d+1];
int indexCounter = 0;
int newIndex, oldLeft, oldRight;
ptIndex[0] = 0;
ptIndex[1] = d;
mu2OverNu = mu * mu / nu;
mu2dTOverNu = mu2OverNu * (t[d] - t[0]);
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]));
bMu2dtOverNuL[newIndex] = mu * mu
* (t[newIndex] - t[ptIndex[oldLeft]]) / nu;
bMu2dtOverNuR[newIndex] = mu * mu
* (t[ptIndex[oldRight]] - t[newIndex]) / nu;
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.
bMu2dtOverNuL[ptIndex[k-1]+1] = mu * mu
* (t[ptIndex[k-1]+1] - t[ptIndex[k-1]]) / nu;
bMu2dtOverNuR[ptIndex[k-1]+1] = mu * mu
* (t[ptIndex[k]] - t[ptIndex[k-1]+1]) / nu;
wIndexList[indexCounter] = ptIndex[k]-2;
wIndexList[indexCounter+1] = ptIndex[k]-1;
wIndexList[indexCounter+2] = ptIndex[k];
indexCounter += 3;
}
}
}
}
}