GammaProcessPCASymmetricalBridge.java

/*
 * Class:        GammaProcessPCASymmetricalBridge
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @authors      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 GammaProcessPCASymmetricalBridge extends GammaProcessPCABridge {

   public GammaProcessPCASymmetricalBridge(double s0, double mu, double nu, RandomStream stream) {
      super(s0, mu, 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.
      // Here we have to hope that the bridge is implemented in the
      // same order for the Brownian and Gamma processes.
 
      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 = BetaSymmetricalDist.inverseF(bMu2dtOverNuL[newIndex], uniformsV[newIndex]);
 
         path[newIndex] = path[oldIndexL] + (path[oldIndexR] - path[oldIndexL]) * y;
      }
      // resetStartProcess();
      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);
   }
 
   // code taken from GammaSymmetricalBridge to check time is power of 2,
   // as it is required for the symmetrical bridge.
   protected void init() {
      super.init();
      if (observationTimesSet) {
 
         /* Testing to make sure number of observations n = 2^k */
         int x = d;
         int y = 1;
         while (x > 1) {
            x = x / 2;
            y = y * 2;
         }
         if (y != d)
            throw new IllegalArgumentException(
                  "GammaSymmetricalBridgeProcess:" + "Number 'n' of observation times is not a power of 2");
      }
   }
}