BrownianMotionBridge.java

/*
 * Class:        BrownianMotionBridge
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * 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 BrownianMotionBridge extends BrownianMotion {
   protected int bridgeCounter = -1; // Before 1st observ
 
   // For precomputations for B Bridge
   protected double[] wMuDt, wSqrtDt;
   protected int[] wIndexList, ptIndex;

   public BrownianMotionBridge(double x0, double mu, double sigma, RandomStream stream) {
      super(x0, mu, sigma, stream);
   }

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