Tools to simulate Markov chains with the Array-RQMC method. More...

Classes
class	ArrayOfComparableChains
	This class provides tools to simulate an array of MarkovChainComparable objects with the array-RQMC method of [139], [141] . More...

class	ArrayOfDoubleChains
	Similar to ArrayOfComparableChains, except that instead of working with $n$ clones of a MarkovChain, we use a single MarkovChainDouble object for all the chains. More...

class	MarkovChain
	This class defines a generic Markov chain and provides basic tools to simulate it for a given number of steps or until it stops, and to recover the performance measure. More...

class	MarkovChainComparable
	A subclass of `MarkovChain` for which there is a total ordering between the states, induced by the implementation of the umontreal.ssj.util.MultiDimComparable interface. More...

class	MarkovChainDouble
	A special kind of Markov chain whose state space is a subset of the real numbers. More...

Detailed Description

Tools to simulate Markov chains with the Array-RQMC method.

This package provides facilities designed specifically to simulate discrete-time Markov chains (DTMC) with the Array-RQMC method [139], [141], [143] . With this method, several realizations of the Markov chain are simulated in parallel. At each step, the chains are sorted (in some way) based on the value of their current state, then all the chains are simulated for one more step using an RQMC point set. We now give more details on the setting.

A DTMC can be written as $\{X_i, i\in I\}$ where $I=\{0,1,2,…\}$, $X_i \in \mathcal{S}$ represents the state at step $i$, and the state space $\mathcal{S}$ is arbitrary. Typically, $\mathcal{S} \subseteq \mathbb{R}^\ell$ for some $\ell\geq 1$. We assume that the state evolves according to the stochastic recurrence

\[ X_0 = x_0, \qquad\mbox{ and } X_j = \varphi(X_{j-1},\mathbf{U}_j) \mbox{ for } j\ge 1, \]

where the $\mathbf{U}_j$ are independent random variables uniformly distributed over $[0,1)^d$ for some integer $d \geq 1$ (it is often 1, but can be larger). A performance mesure $Y_\tau$ is defined over this sequence as

\[ Y = Y_{\tau} = \sum_{j=1}^\tau c_j(X_j), \]

where $c_j$ is a cost (or revenue) function for step $j$ and $\tau$ is either a constant of a random stopping time. Often, $c_j$ does not depend on $j$. The goal is to estimate $\mu=\mathbb E[Y_{\tau}]$. It is not rare that $c_j(\cdot)=0$ for all $j<\tau$, i.e., that the cost or revenue occurs only at the end.

The base class MarkovChain offers methods to simulate one or more copies of the Markov chain one step at a time, or over several steps, collect the realizations of $Y$ in statistical probes, make several independent replications of that, etc. The chains can be simulated by Monte Carlo, RQMC, or Array-RQMC.

To use these tools, one must define a subclass of MarkovChain and implement its three abstract methods: initialState() which resets the chain to its initial state $x_0$; nextStep(stream) which advances the chain by one step from the current state using a random stream (it represents function $\varphi(\cdot)$); and getPerformance() which returns the performance realization for the chain, i.e., the value taken by $Y_{\tau}$, assuming that the chain has reached its stopping time $$.

If the chains have to be sorted as in the Array-RQMC method, one must implement the MarkovChainComparable interface, unless the chain is a subclass of MarkovChainDouble,
in which case the state is just a real number and the chains are sorted by increasing order of their state in a trivial way. For a direct subclass of MarkovChain, other methods are then needed. See the examples below for more details.

The classes ArrayOfComparableChains and ArrayOfDoubleChains work with multiple Markov chains in parallel. They implement Array-RQMC. The chains can be sorted using the method ArrayOfComparableChains.sortChains.

Examples

The following examples demonstrate how to use this package.

Remarks: In the following example, the state is just a real number, so it could be implemented more simply as a MarkovChainDouble.
We should put more elaborate examples, in which the state is multivariate.

The class displayed in Listing Brownian shows a simple implementation of a umontreal.ssj.markovchainrqmc.MarkovChainComparable. It represents a Brownian motion over the real line. The starting position x0 as well as the time step dt are given in the constructor. Each step represents a move which is represented by the addition of a normal variable of mean $0$ and variance dt to the current position. The performance mesure here is just the positive distance between the current position and the initial position, but it could be anything else.

A simple implementation of MarkovChainComparable [Brownian]

package markovchainrqmc;
import umontreal.ssj.markovchainrqmc.MarkovChainComparable;
import umontreal.ssj.rng.RandomStream;
import umontreal.ssj.probdist.NormalDist;
class Brownian extends MarkovChainComparable  {
   final double x0;                   // Initial position.
   final double dt, sqrtDt;           // Time interval between observations.
   double x;                           // Position.              
   public Brownian (double x0, double dt) {
       this.x0 = x0;
       this.dt = dt;
       if (dt < 0)
          throw new IllegalArgumentException("dt must be positive");
       sqrtDt = Math.sqrt (dt);        // Just for faster computation
       stateDim = 1;                   // Dimension of state.
       initialState();
   }
   // Sets initial position
   public void initialState () {
      x = x0;                          
   }
   // Simulates the next step.
   public void nextStep (RandomStream stream) {
      x += sqrtDt * NormalDist.inverseF01 (stream.nextDouble());  
   }
   // Returns performance mesure.
   public double getPerformance () {
      return Math.abs(x-x0);
   }
   
   // Compares value of x between two chains.
   public int compareTo (MarkovChainComparable m, int i) {
      if (!(m instanceof Brownian))     
         throw new IllegalArgumentException("Can't compare a " +
           "Brownian Markov chain with other types of Markov chains."); 
      switch(i) {
         case 0: 
            double mx = ((Brownian)m).x;
            return (x>mx ? 1 : (x<mx ? -1 : 0));
         default: 
            throw new AssertionError("Invalid state index" );
      }
   }
}

The program displayed in Listing BrownianTest shows different ways to use the Markov chain.

1- How to simulate the trajectory and print the state of the chain at each step and the performance at the end.

2- How to simulate using Monte Carlo to get an unbiased estimator of the expectation of the performance and an estimation of its variance. If stream is a umontreal.ssj.hups.PointSetIterator, use umontreal.ssj.markovchainrqmc.MarkovChain.simulRunsWithSubstreams instead of umontreal.ssj.markovchainrqmc.MarkovChain.simulRuns. The umontreal.ssj.stat.Tally is a statistical collector; see package umontreal.ssj.stat for how to use it.

3- Same as 2 but with randomized quasi-Monte Carlo. Basically, it takes a umontreal.ssj.hups.PointSet where the dimension of the points is the number of steps and the number of points is the number of trajectories. The umontreal.ssj.hups.PointSetRandomization must be compatible with the point set. See package umontreal.ssj.hups more information on these classes.

4- Same as 2 but with the array-RQMC method of [139] . The ArrayOfComparableChains is used to simulate chains in parallel. It uses a umontreal.ssj.hups.PointSetRandomization to randomize the point sets and a umontreal.ssj.util.MultiDimSort to sort the chains. Here, as the chain is one-dimensional, the sort used is a umontreal.ssj.util.OneDimSort. It is important to call umontreal.ssj.markovchainrqmc.ArrayOfComparableChains.makeCopies(int) in order to set the number of chains. See package umontreal.ssj.util for more information on sorts.

5- How to simulate the trajectories with array-RQMC and do something with the chains at each step. The Do something with mc comment should be replaced by anything, using the MarkovChain mc. For example to store or print the state x of each chain for a later use.

Tests using a MarkovChainComparable [BrownianTest]

package markovchainrqmc;
import umontreal.ssj.markovchainrqmc.*;
import umontreal.ssj.hups.*;
import umontreal.ssj.rng.*;
import umontreal.ssj.stat.Tally;
import umontreal.ssj.util.*;
import umontreal.ssj.util.sort.*;
public class BrownianTest {
   Brownian brownian;
   
   public BrownianTest (double x0, double dt) {
      brownian = new Brownian(x0,dt);
      tests();
   }
   
   public void tests() {
      RandomStream stream = new MRG32k3a();
      
      //1- Print trajectory and performance for 20 steps of the chain.
      System.out.println("1- Print trajectory and performance");
      brownian.initialState ();
      System.out.println("step =  0,   position =  " + brownian.x);
      for (int step = 1; step < 21; ++step){
          brownian.nextStep (stream);
          System.out.printf("step = %2d,   position = %8.5f%n", step, brownian.x);
      }
      System.out.printf("Performance = %8.5f%n", brownian.getPerformance());
      
      //2- Monte Carlo, 100 replications, 2^12 trajectories, 20 steps 
      Tally performance = new Tally("Performance of Brownian Motion");
      Tally replicates = new Tally("Replicates of Brownian Motion");
      for(int i=0; i<100; ++i){
         brownian.simulRuns(4096,20,stream,performance);
         stream.resetNextSubstream();
         replicates.add(performance.average());
      }
      System.out.println("\n 2- MC, 100 reps, 2^12 trajetories, 20 steps ");
      System.out.println(replicates.report());
      
      //3- RQMC, 100 replications, 2^12 trajectories, 20 steps 
      PointSetRandomization rand = new RandomShift(stream);
      PointSet p = new SobolSequence(12,31,20);
      brownian.simulRQMC(p,100,20,rand,replicates);
      System.out.println("\n 3- RQMC, 100 reps, 2^12 trajectories, 20 steps");
      System.out.println(replicates.report());
      
      //4- Array-RQMC, 100 replications, 2^12 trajectories, 20 steps
      PointSet p2 = new SobolSequence(12,31,1);
      MultiDimSort sort = new OneDimSort(0);
      ArrayOfComparableChains array = 
                              new ArrayOfComparableChains(brownian, rand, sort);
      array.makeCopies(4096);
      array.simulReplicatesArrayRQMC (p2, rand, sort, 0, 100, 20, replicates);
      System.out.println("\n 4- Array-RQMC, 100 replications" + 
                                              ", 2^12 trajectories, 20 steps");
      System.out.println(replicates.report());
      
      //5- Array-RQMC Simulate and use 2^12 trajectories, 20 steps 
      array.makeCopies(4096);
      array.initialStates();
      for (int step = 1; step < 21; ++step){
          array.sortChains();
          array.simulOneStepArrayRQMC (p2);
          for(MarkovChainComparable mc: array.getChains()){
             //Do something with mc
          }
      }
   }
   public static void main (String[] args) {
      double x0 = 0.0;
      double dt = 0.1;       
      BrownianTest testBrownian = new BrownianTest(x0,dt);
   }
}

The output of this program is displayed in Listing BrownianTestOutput. For this example, the variance of the estimator with RQMC is 6.25 times less than MC, and 388 times less with array-RQMC compared to MC.

Output of BrownianTest.java [BrownianTestOutput]

1- Print trajectory and performance
step =  0,   position =  0.0
step =  1,   position = -0.36070
step =  2,   position = -0.50990
step =  3,   position = -0.66743
step =  4,   position = -0.37085
step =  5,   position = -0.61330
step =  6,   position = -0.58680
step =  7,   position = -0.60205
step =  8,   position = -0.71916
step =  9,   position = -1.06654
step = 10,   position = -0.84739
step = 11,   position = -0.78714
step = 12,   position = -0.85904
step = 13,   position = -1.00137
step = 14,   position = -1.22434
step = 15,   position = -1.13596
step = 16,   position = -0.72304
step = 17,   position = -0.88980
step = 18,   position = -1.46628
step = 19,   position = -0.88737
step = 20,   position = -1.22407
Performance =  1.22407
 2- MC, 100 reps, 2^12 trajetories, 20 steps 
REPORT on Tally stat. collector ==> Replicates of Brownian Motion
    num. obs.     min         max       average    standard dev.
       100       1.098       1.166       1.130       0.013
 3- RQMC, 100 reps, 2^12 trajectories, 20 steps
REPORT on Tally stat. collector ==> Replicates of Brownian Motion
    num. obs.     min         max       average    standard dev.
       100       1.115       1.140       1.127      5.2E-3
 4- Array-RQMC, 100 replications, 2^12 trajectories, 20 steps
REPORT on Tally stat. collector ==> Replicates of Brownian Motion
    num. obs.     min         max       average    standard dev.
       100       1.126       1.130       1.128      6.6E-4

Classes

Detailed Description

Examples