umontreal.iro.lecuyer.stochprocess.BrownianMotion Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer, in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal. It provides facilities for generating uniform and nonuniform random variates, computing different measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both events and processes.
The newest version!


/*
 * Class:        BrownianMotion
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Université de Montréal
 * Organization: DIRO, Université de Montréal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   GPL licence site.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * This class represents a Brownian motion process 
 * {X(t) : t >= 0},
 * sampled at times 
 * 0 = t₀ < t₁ < ^... < t_d.
 * This process obeys the stochastic differential equation
 * 
 * 
 * 
 * dX(t) = μdt + σdB(t),
 * 

 * with initial condition X(0) = x₀,
 * where μ and σ are the drift and volatility parameters,
 * and 
 * {B(t), t >= 0} is a standard Brownian motion
 * (with drift 0 and volatility 1).
 * This process has stationary and independent increments over disjoint
 * time intervals (it is a Lévy process) and the increment over an interval
 * of length t is normally distributed with mean μt and variance 
 * σ²t.
 * 
 * 
 * In this class, this process is generated using the sequential (or random walk)
 * technique:  X(0) = x₀ and
 * 
 * 

 * 
 * X(t_j) - X(t_j-1) = μ(t_j - t_j-1) + σ(t_j - t_j-1)^1/2Z_j
 * 

 * where 
 * Z_j∼N(0, 1).
 * 
 */
public class BrownianMotion extends StochasticProcess  {
    protected NormalGen    gen;
    protected double       mu,
                           sigma;
    // Precomputed values for standard BM
    protected double[]     mudt,
                           sigmasqrdt;



   /**
    * Constructs a new BrownianMotion with
    * parameters μ = mu, σ = sigma and initial value
    * 
    * X(t₀) = x0.
    * The normal variates Z_j in will be
    * generated by inversion using stream.
    * 
    */
   public BrownianMotion (double x0, double mu, double sigma,
                          RandomStream stream) {
        this (x0, mu, sigma, new NormalGen (stream));
    }


   /**
    * Constructs a new BrownianMotion with parameters μ =
    * mu, σ = sigma and initial value 
    * X(t₀) = x0.
    * Here, the normal variate generator
    * {@link umontreal.iro.lecuyer.randvar.NormalGen NormalGen} is specified
    * directly instead of specifying the stream and using inversion.
    * The normal generator gen can use another method than inversion.
    * 
    */
   public BrownianMotion (double x0, double mu, double sigma, NormalGen gen) {
        this.mu    = mu;
        this.sigma = sigma;
        this.x0    = x0;
        this.gen   = gen;
    }
 

   public double nextObservation() {
        double x = path[observationIndex];
        x += mudt[observationIndex]
             + sigmasqrdt[observationIndex] * gen.nextDouble();
        observationIndex++;
        path[observationIndex] = x;
        return x;
    }


   /**
    * Generates and returns the next observation at time t_j+1 =
    *  nextTime. It uses the previous observation time t_j defined earlier
    * (either by this method or by setObservationTimes),
    * as well as the value of the previous observation X(t_j).
    * Warning: This method will reset the observations time t_j+1
    * for this process to nextTime. The user must make sure that
    * the t_j+1 supplied is 
    *  >= t_j.
    * 
    */
   public double nextObservation (double nextTime)  {
        // This method is useful for generating variance gamma processes
        double x = path[observationIndex];
        double previousTime = t[observationIndex];
        observationIndex++;
        t[observationIndex] = nextTime;
        double dt = nextTime - previousTime;
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        path[observationIndex] = x;
        return x;
    }


   /**
    * Generates an observation of the process in dt time units,
    * assuming that the process has value x at the current time.
    * Uses the process parameters specified in the constructor.
    * Note that this method does not affect the sample path of the process
    * stored internally (if any).
    * 
    */
   public double nextObservation (double x, double dt)  {
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        return x;
    }


   public double[] generatePath() {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * gen.nextDouble();
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * Same as generatePath(), but a vector of uniform random numbers
    * must be provided to the method.  These uniform random numbers are used
    * to generate the path.
    * 
    */
   public double[] generatePath (double[] uniform01)  {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * NormalDist.inverseF01(uniform01[j]);
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   public double[] generatePath (RandomStream stream) {
        gen.setStream (stream);
        return generatePath();
    }

   /**
    * Resets the parameters 
    * X(t₀) = x0, 
    * μ = mu and
    * 
    * σ = sigma of the process.
    * Warning: This method will recompute some quantities stored internally,
    * which may be slow if called too frequently.
    * 
    */
   public void setParams (double x0, double mu, double sigma)  {
        this.x0    = x0;
        this.mu    = mu;
        if (sigma <= 0)
           throw new IllegalArgumentException ("sigma <= 0");
        this.sigma = sigma;
        if (observationTimesSet) init(); // Otherwise not needed.
    }


   /**
    * Resets the random stream of the normal generator to stream.
    * 
    */
   public void setStream (RandomStream stream)  { gen.setStream (stream); }


   /**
    * Returns the random stream of the normal generator.
    * 
    */
   public RandomStream getStream()  { return gen.getStream (); }


   /**
    * Returns the value of μ.
    * 
    */
   public double getMu()  { return mu; }


   /**
    * Returns the value of σ.
    * 
    */
   public double getSigma()  { return sigma; }


   /**
    * Returns the normal random variate generator used.
    * The {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    * used by that generator can be changed via
    * getGen().setStream(stream), for example.
    * 
    */
   public NormalGen getGen()  { return gen; }
 

    // This is called by setObservationTimes to precompute constants
    // in order to speed up the path generation.
   protected void init() {
        super.init();
        mudt       = new double[d];
        sigmasqrdt = new double[d];
        for (int j = 0; j < d; j++) {
            double dt     = t[j+1] - t[j];
            mudt[j]       = mu * dt;
            sigmasqrdt[j] = sigma * Math.sqrt (dt);
        }
     }

}