GammaProcessPCABridge.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-Sébastien 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.*;

public class GammaProcessPCABridge extends GammaProcessPCA {
   protected BrownianMotionBridge BMBridge;
   protected double mu2OverNu, mu2dTOverNu;
   protected double[] bMu2dtOverNuL, // For precomputations for G Bridge
         bMu2dtOverNuR;
   protected int[] wIndexList;

   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);
   }

   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;
            }
         }
      }
   }
 
}