Chapter 17

Introduction to Monte Carlo and Simulation Methods

17.1 Introduction

The Monte Carlo method is a numerical technique that allows scientists to analyze various natural phenomena and compute complicated quantities by means of repeated generation of random numbers. For example, if you make several thousand tosses of a fair (i.e., balanced) coin, then you may notice that the long-run ratio of a count of heads to the total number of tosses is approaching one half. If the limiting ratio is not close to one half, then you may conclude that the coin is not balanced. Similarly, researchers can compute other more complicated quantities, although in reality nobody tosses coins and throws dice. Instead, computer simulations are used. First, a computer experiment that produces some quantitative information about a random phenomenon of interest needs to be designed. After that, a computer is used to perform independent repeated random simulations and then to calculate averages of results. Such averages are used to approximate the quantity of interest. In particular, the two fundamental theorems, namely, the Law of Large Numbers (LLN) and the Central Limit Theorem (CLT), allow researchers to construct a valid approximation of the quantity of interest and estimate the approximation error.

The Monte Carlo method was coined in the 1940s by John von Neumann, Stanislaw Ulam, and Nicholas Metropolis, while they were working on nuclear weapons project at the Los Alamos National Laboratory. It was named after a famous casino in Monte Carlo, where Ulam’s uncle often gambled away his money. The range of problems that can be analyzed by Monte Carlo methods is vast. Many quantities such as the probability of default of a financial company, distribution of galaxies in the universe, and characteristics of a nuclear reactor can be computed. Monte Carlo methods are especially useful for simulating complex systems with coupled degrees of freedom and with significant uncertainty in inputs, such as the calculation of risk in business. Monte Carlo methods are used to evaluate multidimensional integrals and to numerically solve large systems of equations. The success and popularity of Monte Carlo is partly explained by the enormous progress of computers that are so powerful and inexpensive these days.

The Monte Carlo method is a very popular computational tool in financial economics. Many typical problems such as the optimal allocation of financial assets, pricing of derivative contracts, and evaluation of business risk can be numerically solved by means of repeated generation of possible market scenarios. For example, we can calculate the no-arbitrage value of an option written on the asset by repeatedly generating sample paths of the underlying asset price process. The simplest possible model is the binomial price model. At each time step, the price may go up or down by a constant factor. Every possible scenario can be represented by a random walk on the binomial tree. An estimate of the option price is computed by averaging payoff values calculated for independent realizations of such a random walk.

17.1.1 The “Hit-or-Miss” Method

A typical application of the Monte Carlo method is the computation of areas and volumes of objects of complex shape and geometry. This type of Monte Carlo computation is based on two ideas: geometric probability and the frequentist definition of probability. The latter tells us that the likelihood of some event E can be calculated as a long-run ratio of the number of successful trials to the total number of trials. By a success we mean here a trial where the event E occurs. Count the number of times, NE(n), the event E occurs in the first n trials that are performed independently and under identical conditions. The probability ℙ(E) is approximated as

$\mathbb{P}\left(E\right)\approx \frac{N_E\left(n\right)}{n}.$

Now let us link the notion of probability to the area of a figure (or the volume of a multidimensional domain in the general case). Consider a planar figure D contained completely within a unit square S = [0, 1]2. Select a point at random uniformly on the square. This can be done by sampling two independent identically distributed (i.i.d.) Cartesian coordinates uniformly in the interval [0, 1]. According to the principle of geometric probability, the likelihood that such a randomly chosen point belongs to D is equal to the ratio of the area of D, denoted |D|, to the area of the unit square (which equals one). Choose at random n points independently and uniformly on the square. Let ND(n) points fall on the figure D. On the one hand, the ratio ND(n)/n is approximately equal to the probability that a random point chosen uniformly on the square lies on D. On the other hand, this probability equals the ratio of areas of D and S. Thus, we obtain the following approximation of the area of D:

$\left|D\right|\approx \frac{N_D\left(n\right)}{n}\cdot \left|S\right|.$

17.1.2 The Law of Large Numbers

The Monte Carlo method is based on two fundamental laws, namely, the Law of Large Numbers (LLN) and the Central Limit Theorem (CLT). There are several versions of the LLN. Here are some of them.

Theorem 17.1

(Borel’s LLN). Suppose that an experiment with an uncertain outcome is repeated a large number of times independently and under identical conditions. Then, for any event E we have

$\mathbb{P}\left(E\right)=\underset{n\to \infty }{\lim }\frac{N_E\left(n\right)}{n},$

where NE(n) is the number of times the event E occurs in the first n trials.

Borel’s LLN provides us with a theoretical foundation for the frequintist interpretation of probability. On the other hand, the probability of an event E can be expressed as a mathematical expectation of an indicator function of $E:\mathbb{P}\left(E\right)=\text{E}\left[{\mathbb{I}}_E\right]$. Thus, Borel’s LLN is a special case of a general LLN provided just below.

Let {Xk}k≥1 be a sequence of i.i.d. random variables with common finite mean μX, which are all defined on some probability space (Ω, ℱ, ℙ). Let ${\overline{X}}_n=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{X}_k}$ denote the arithmetic average of the first n variables.

Theorem 17.2

(Chebyshev’s Weak LLN). It is true that ${\overline{X}}_n\stackrel{p}{\to }{\mu }_X,asn\to \infty$. That is, $\forall \epsilon >0\mathbb{P}\left(\left|{\overline{X}}_n-{\mu }_X\right|>\epsilon \right)\to 0,asn\to \infty$.

Theorem 17.3

(Kolmogorov’s Strong LLN). It is true that ${\overline{X}}_n\stackrel{a.s.}{\to }{\mu }_X$, as n → ∞. That is, $\mathbb{P}\left(\left\{\omega |{\overline{X}}_n\left(\omega \right)\to {\mu }_{X}asn\to \infty \right\}\right)=1$.

The LLN provides us with a recipe for estimation of quantities of interest. For a given unknown Q, we first construct a random variable X so that Q = μX = E[X], i.e., we find the probabilistic representation of the quantity Q. The random variable X is called an estimator of the quantity Q. It is an unbiased estimator of Q, meaning that E[X] = Q. Let X1, X2, . . . , Xn be i.i.d. variates all having the same probability distribution as that of X. Following the LLN, construct a sample estimator of Q, ${\overline{X}}_n=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{X}_k}$, which converges almost surely to Q, as n → ∞. That is, the quantity Q can be approximated by an average of n independent sample values:

$Q\approx {\overline{x}}_n=\frac{1}{n}{\displaystyle \sum _{k=1}^n{x}_k},$

where xk = Xk(ω), 1 ≤ k ≤ n, for some outcome ω ∊ Ω. To obtain independent sample values {xk}k≥1, it suffices to construct a sequence of statistically independent numbers (or almost statistically independent pseudorandom numbers) sampled from the target distribution of X as follows. First, we sample from the uniform distribution Unif(0, 1). After that, the independent sample values (also called realizations or draws) xk, 1 ≤ k ≤ n, can be generated from a sequence of uniform random numbers by using a special transformation algorithm. One example of such algorithms is the inverse cumulative distribution function (CDF) method, which is based on the representation $X\stackrel{d}{=}{F}_X^{-1}\left(U\right)$, where $F_X^{-1}$ is the (generalized) inverse of the CDF of X and U ~ Unif(0, 1).

In what follows, we will distinguish the sample mean estimator ${\overline{X}}_n$, which is a random quantity equal to an average of n i.i.d. variates, and a sample mean estimate ${\overline{x}}_n$, which is a nonrandom quantity equal to an average of n statistically independent realizations of the estimator X. While the former is important for the theoretical analysis (e.g., to construct a confidence interval), the latter is used in practice to approximate Q.

17.1.3 Approximation Error and Confidence Interval

The next goal is to estimate the approximation error. The Central Limit Theorem (CLT) provides a solution to this problem.

Theorem 17.4

(The CLT). Consider a sequence {Xk}k≥1 of i.i.d. variates with finite common variance ${\sigma }_X^2$ and expected value μX. Then,

$\frac{{\overline{X}}_n-{\mu }_X}{{\sigma }_X/\sqrt{n}}\underrightarrow{d}Norm\left(0,1\right),asn\to \infty .$

That is, $\mathbb{P}\left(\frac{{\overline{X}}_{n-\mu X}}{{\sigma }_X/\sqrt{n}}\le z\right)\to \mathcal{N}\left(z\right)$, as n → ∞, for all z ∊ ℝ.

With the help of the CLT, we can construct a confidence interval for the mathematical expectation μX. We have that

$\mathbb{P}\left(\left|\frac{{\overline{X}}_n-{\mu }_X}{{\sigma }_X/\sqrt{n}}\right|\le z\right)\to \mathbb{P}\left(\left|Z\right|\le z\right)=2\mathcal{N}\left(z\right)-1\text{for}\text{all}\text{z}\in ℝ\text{.}$

We wish to make this probability as close to one as possible. Let us fix a confidence level 1 − α ∊ (0, 1) with α ≪ 1. Solve the equation

$2\mathcal{N}\left(z\right)-1=1-\alpha \iff \mathcal{N}\left(z\right)=1-\frac{\alpha }{2}$

for z. Since $\mathcal{N}$ is a strictly monotone function of z, the solution is unique. It is a so-called (1 − α/2)-quantile of Norm(0, 1) denoted by zα/2. Table 17.1 contains commonly used normal quantiles.

For large values of n, the distribution of ${\overline{X}}_n$ is approximately normal. By the CLT we obtain

$\left|{\overline{X}}_n-{\mu }_X\right|\le \frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}}\text{with}\text{probability}\approx 1-\alpha .$

Hence, $\mathbb{P}\left(({\overline{X}}_n-\frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}},{\overline{X}}_n+\frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}})∍{\mu }_X\right)\approx 1-\alpha$. Replacing the average ${\overline{X}}_n$ by its sample value ${\overline{x}}_n$ gives us a confidence interval for μX:

$\left({\overline{x}}_n-\frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}},{\overline{x}}_n+\frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}}\right)∍{\mu }_X\text{with}\text{the}\text{confidence}\text{level}\text{of}\left(1-\alpha \right).$

Typically, the variance ${\sigma }_X^2$ is unknown. However, it can be approximated by the sample variance:

${\sigma }_X^2\approx s_n^2:=\frac{1}{n}{\displaystyle \sum _{k=1}^n{\left(x_k-{\overline{x}}_n\right)}^2}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{x}_K^2-{\overline{x}}_n^2}.$

As is seen from the formula of the confidence interval, the relative error $\frac{z_{\alpha /2}{\sigma }_X}{\sqrt{n}{\mu }_X}$ is of order $\mathcal{O}\left(n^{-0.5}\right)$, as n → ∞. This fact points out the main drawback of the Monte Carlo method—a slow rate of convergence. For example, to decrease the error by a factor of 10, the number of sample values needs to be increased to 100 times. The quantity $\frac{s_n/\sqrt{n}}{{\overline{x}}_n}$ is often referred to as an accuracy measure for the sample estimate ${\overline{x}}_n$.

17.1.4 Parallel Monte Carlo Methods

One of the main advantages of the Monte Carlo method is the ease of its parallelization. In the simplest case, independent CPUs can calculate partial expectations and then a head CPU computes the final average and the confidence interval. Consider a computational cluster, where CPUs are numbered from 1 through ℓ. Let CPU # i generate ni independent draws, $\left\{x_1^i,...,x_{n_i}^i\right\}$. The partial averages

${\overline{x}}^{\left(i\right)}=\frac{x_1^i+x_2^i+\cdots +x_{n_i}^i}{n_i}\text{and}{\overline{y}}^{\left(i\right)}=\frac{{\left(x_1^i\right)}^2+{\left(x_2^i\right)}^2+\cdots +{\left(x_{n_i}^i\right)}^2}{n_i}$

are calculated and then their values are sent to the head CPU. The total number of draws is n = n1 + n2 + · · · + nk. To construct a confidence interval for the mean value, we only need to compute the sample mean ${\overline{x}}_n$ and the sample variance $s_n^2$ as follows:

$\begin{array}{l}{\overline{x}}_n=\frac{1}{n}{\displaystyle \sum _{i=1}^{\ell }{n}_i{\overline{x}}^{\left(i\right)}}=\frac{1}{n}{\displaystyle \sum _{i=1}^{\ell }{\displaystyle \sum _{j=1}^{n_i}{x}_j^i}}\hfill \\ s_n^2=\frac{1}{n}{\displaystyle \sum _{i=1}^{\ell }{n}_i{\overline{y}}^{\left(i\right)}-{\overline{x}}_n^2}=\frac{1}{n}{\displaystyle \sum _{i=1}^{\ell }{\displaystyle \sum _{j=1}^{n_i}{\left(x_j^i\right)}^2-{\overline{x}}_n^2}}.\hfill \end{array}$

Note that the numbers ni, i = 1, 2, . . . , ℓ, can be chosen so that all CPUs will finish their jobs almost simultaneously. If the CPUs we deal with have a similar performance, then we can simply set n1 = n2 = · · · = nℓ.

17.1.5 One Monte Carlo Application: Numerical Integration

Suppose that we wish to evaluate a definite integral of an integrable function g : ℝ → ℝ on a finite interval $\left[a,b\right]:I={\displaystyle {\int }_a^{b}g\left(x\right)\text{d}x}$. Let us consider several Monte Carlo methods of approximating I.

The “hit-or-miss” method is based on the geometric interpretation of a definite integral. Suppose that the function g is nonnegative and bounded: 0 ≤ g ≤ M. If a point with Cartesian coordinates (X, Y) is chosen uniformly on the rectangle R = [a, b] × [0, M], then the probability of the event {Y ≤ g(X)} (i.e., the event that (X, Y) belongs to the domain D bounded by the graph of g and the lines x = a, x = b, and y = 0) is equal to the ratio of two areas, $\frac{\left|D\right|}{\left|R\right|}$. Since |D| = I and |R| = (b − a)M, we obtain the following approximation:

$I\approx \left(b-a\right)M\frac{N_D\left(n\right)}{n},$

where n statistically independent random points (xk, yk), 1 ≤ k ≤ n, are sampled uniformly on the rectangle R, and ND(n)/n is the fraction of points belonging to D. This result can be viewed as an application of the LLN. Let us introduce an indicator function ${\mathbb{I}}_D\left(X,Y\right)$, which equals 1 if (X, Y) ∊ D and equals 0 otherwise. Then, the expected value of this indicator is $\text{E}\left[{\mathbb{I}}_D\left(X,Y\right)\right]=\mathbb{P}\left(Y\le g\left(X\right)\right)=\frac{|D|}{|R|}$. Therefore, $\text{E}\left[\left|R\right|\right].{\mathbb{I}}_D\left(X,Y\right)]=\left|D\right|=I$ and applying the LLN gives

$\frac{1}{n}{\displaystyle \sum _{k=1}^n\left|R\right|}{\mathbb{I}}_D\left(x_k,y_k\right)=\left(b-a\right)M\frac{N_D\left(n\right)}{n}\to I,asn\to \infty .$

For the sample mean method, the representation of the integral I in the form of a mathematical expectation w.r.t. a uniform probability distribution is used:

$I=\left(b-a\right){\displaystyle {\int }_a^{b}g\left(x\right)}\frac{1}{b-a}\text{d}x=\text{E}\left[\left(b-a\right)g\left(X\right)\right],$

where X ~ Unif(a, b). Therefore, by the LLN we have $I\approx \frac{1}{n}{\displaystyle \sum _{k=1}^n\left(b-a\right)g\left(x_k\right)}$, where {xk}k≥1 are independent draws from Unif(a, b). In comparison with the “hit-or-miss” method, g is only required to be integrable on the interval (a, b).

The weighted method generalizes the previous approach. Suppose that there exists a probability density function (PDF) f whose support is (a, b), where (a, b) can be a finite, semi-infinite, or infinite interval. Then, the integral I can be represented in the form of a mathematical expectation w.r.t. the PDF f:

$I={\displaystyle {\int }_a^b\frac{g\left(x\right)}{f\left(x\right)}f\left(x\right)}\text{d}x=\text{E}\left[\frac{g\left(X\right)}{f\left(X\right)}\right],\text{where}X~f,$

provided that the ratio $\frac{g\left(x\right)}{f\left(x\right)}$ is defined for all x ∊ (a, b) and is an integrable function of x.

The random variable $\frac{g\left(X\right)}{f\left(X\right)}$ is the estimator of I. Thus, we approximate $I\approx \frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{g\left(x_k\right)}{f\left(x_k\right)}}$ where {xk}k≥1 are independent draws from the PDF f.

Clearly, there are many PDFs that can be used in the weighted method. The importance sampling principle, which is discussed in Section 17.5, explains how to chose f. The optimal PDF f that minimizes the variance of the random estimator $\frac{g\left(X\right)}{f\left(X\right)}$, where X ~ f, i.e., the PDF that solves the minimization problem

$\text{Var}\left(\frac{g\left(X\right)}{f\left(X\right)}\right)\to \underset{f}{\min },$

is a function proportional to |g|:

$f\left(x\right)=\frac{\left|g\left(x\right)\right|}{{\displaystyle {\int }_a^b\left|g\left(t\right)\right|\text{d}t}},x\in \left(a,b\right).$

17.2 Generation of Uniformly Distributed Random Numbers

A central part of every stochastic simulation algorithm is a generator of random numbers, which produces a sequence of statistically independent samples (or draws) from a given distribution. As we demonstrate in the following sections, the sampling from a nonuniform probability distribution can be reduced to the sampling from the uniform distribution on (0, 1), denoted Unif(0, 1), by applying certain procedures such as transformation methods and acceptance-rejection techniques. Our main priority is to have available a reliable method of sampling from Unif(0, 1). Therefore, we first concentrate on methods of generating statistically independent (pseudo-)random numbers uniformly distributed on (0, 1).

The process of obtaining truly random numbers can be quite complicated. There exist different types of physical generators of random numbers. The simplest ones are balanced coins and dice, and even a playing roulette. More advanced generators are based on the use of physical phenomena such as thermal noise. However, there are common drawbacks of such hardware generators such as their slow speed, lack of portability, possible nonuniformness of numbers generated, and impossibility of reproducing the same sequence of draws. Algorithmic (software) generators of random numbers solve most of these issues although they create new ones. Numbers generated by such algorithms only mimic truly random numbers, hence the generated draws are called pseudo-random numbers (PRNs). However, the statistical properties of software generators are put to scrutiny so we can trust the numbers obtained.

17.2.1 Uniform Probability Distributions

As is well known, the distribution of a continuous random variable X can be characterized by its PDF, denoted fX, or a cumulative distribution function (CDF), denoted FX. For a continuous random variable U uniformly distributed on (0, 1),

$f_U\left(x\right)={\mathbb{I}}_{\left(0,1\right)}\left(x\right)=\{\begin{array}{l}1\text{if}0<x<1,\hfill \\ 0\text{otherwise},\hfill \end{array}{F}_U\left(X\right)={\displaystyle {\int }_{-\infty }^x=\{\begin{array}{l}0x\le 0,\hfill \\ x0<x<1,\hfill \\ 11\le x.\hfill \end{array}}$

The mathematical expectation and variance of U are, respectively,

$\text{E}\left[U\right]={\displaystyle {\int }_{-\infty }^{\infty }xfu\left(x\right)}\text{d}x={\displaystyle {\int }_0^{1}x\text{d}x}=\frac{1}{2},\text{Var}\left(U\right)=\text{E}\left[U^2\right]-{\left(\text{E}\left[U\right]\right)}^2=\frac{1}{3}-\frac{1}{4}=\frac{1}{12}.$

The continuous uniform distribution on (a, b) with a < b, denoted Unif(a, b), reduces to the case with the interval (0, 1) as follows. Let X ~ Unif(a, b) and U ~ Unif(0, 1). Then, we have

$X\stackrel{d}{=}a+\left(b-a\right)U,F_X\left(x\right)=F_U\left(\frac{x-a}{b-a}\right),f_X\left(x\right)=\frac{1}{b-a}{f}_U\left(\frac{x-a}{b-a}\right)=\frac{1}{b-a}{\mathbb{I}}_{\left(a,b\right)}\left(x\right).$

Now consider the multidimensional case. Suppose that a point X is chosen at random in a domain D ⊂ ℝm with a finite volume |D| so that X is equally likely to lie anywhere in D. Then, the random vector X = [X1, X2, . . . , Xm]⊤ is said to have a multivariate uniform distribution in D; its multivariate PDF is

$f_X\left(\text{x}\right)=\frac{1}{|D|}{\mathbb{I}}_D\left(\text{x}\right)=\{\begin{array}{l}\frac{1}{|D|}\text{x}\in \text{D,}\hfill \\ 0\text{otherwise,}\hfill \end{array}\text{x}\text{=}{\left[x_1{x}_2\cdots x_m\right]}^{⊺}\in D.$

In the case of a hyperparallelepiped D = (a1, b1) × (a2, b2) × · · · × (am, bm) with ai < bi, 1 ≤ i ≤ m, we have

$f\text{x}\left(\text{x}\right)={\displaystyle \prod _{i=1}^{m}f{\text{x}}_i}\left(x_i\right)={\displaystyle \prod _{i=1}^m\frac{1}{b_i-a_i}}{\mathbb{I}}_{\left(a_i,b_i\right)}\left(x_i\right).$

Since the multivariate PDF fX is a product of m univariate PDFs $f_{X_i}=\frac{1}{b_i-a_i}{\mathbb{I}}_{\left(a_i,b_i\right)}$, the entries of the vector X are independent uniformly distributed random variables. Therefore, the simulation of a vector uniformly distributed in a hyperparallelepiped reduces to sampling from Unif(0, 1).

17.2.2 Linear Congruential Generator

Many algorithms for generating pseudo-random numbers (PRNs for short) uniformly distributed in [0, 1] have the form of an iterative rule, ut+1 = F (ut), t = 0, 1, . . ., where both the range and the domain of the function F are [0, 1] and the initial seed u0 ∊ [0, 1] is given. Suppose that a sequence of PRNs {ut}t≥0 is generated by such a rule. Combine the numbers in pairs to obtain points (ut, ut−1) ∊ [0, 1]2, t = 1, 3, 5, . . . . On the one hand, these points are situated on the curve y = F(x). On the other hand, they have to be uniformly distributed in the unit square [0, 1]2 so the points will fill out the square without leaving gaps. Therefore, the plot of F should provide a sufficiently dense filling of the square. One example of a function that has such a property is the mapping y = {Mx}, where {x} denotes the fractional part of x, i.e., $\left\{x\right\}=x-\lfloor x\rfloor$, where $\lfloor \cdot \rfloor$ is the floor function, and M is a large positive number called a multiplier. The plot of y = {Mx} consists of parallel segments and the distance between them goes to zero as M → ∞ (see Figure 17.2).

Proposition 17.5

(Voitishek and Mikhaĭlov (2006)). Consider the transformation y = {Mx + a} where {x} denotes the fractional part of x, and M ∊ Z and a ∊ ℝ are positive constants.

1. If U ~ Unif(0, 1), then {MU + a} ~ Unif(0, 1).
2. Let U0 ~ Unif(0, 1) and the sequence {Uk}k≥1 be generated from the rule Uk+1 = {MUk}. Then, Uk ~ Unif(0, 1) and Corr(U0, Uk) = Corr(Un, Un+k) = M−k for all n ≥ 0 and k ≥ 0.

Proof. First, we prove the fact that X := {MU} ~ Unif(0, 1). By the definition of the fractional part function, we have X ∊ [0, 1). Moreover, X = 0 iff MU is an integer what happens with probability zero. Hence ℙ(X ∊ (0, 1)) = 1. For x ∊ (0, 1), we have

$\mathbb{P}\left(X\le x\right)={\displaystyle \sum _{k=0}^{M-1}\mathbb{P}\left(k\le MU\le k+x\right)}={\displaystyle \sum _{k=0}^{M-1}\mathbb{P}}\left(\frac{k}{M}\le U\le \frac{k}{M}+\frac{x}{M}\right)={\displaystyle \sum _{k=0}^{M-1}\frac{x}{M}=x.}$

Here, we use the fact that 0 ≤ MU ≤ M and {z} ∊ [0, x] iff $z\in {\cup }_{k\in \mathbb{Z}}\left[k,k+x\right]$. Thus, the CDF of X is FX(x) = x for x ∊ (0, 1). Therefore, X ~ Unif(0, 1).

Second, let us prove that {MU + a} ~ Unif(0, 1). Note that

$\left\{MU+a\right\}=\left\{\left\{MU\right\}+\left\{a\right\}\right\}\stackrel{d}{=}\left\{U+b\right\},$

where b = {a} ∊ [0, 1). Show that ℙ({U + b} ≤ x) = x for x ∊ (0, 1). Consider two cases. (x ≤ b): Since {U + b} ≤ x iff 1 ≤ U + b ≤ 1 + x, we have

$\mathbb{P}\left(\left\{U+b\right\}\le x\right)=\mathbb{P}(\underset{\in \left(0,1\right)}{\underbrace{1-b}}\le U\le \underset{\in \left(0,1\right)}{\underbrace{1+x-b}})=\left(1+x-b\right)-\left(1-b\right)=x.$

(b < x): Since {U + b} ≥ x iff x ≤ U + b ≤ 1, we have

$\mathbb{P}\left(\left\{U+b\right\}\ge x\right)=\mathbb{P}(\underset{\in \left(0,1\right)}{\underbrace{x-b}}\le U\le \underset{\in \left(0,1\right)}{\underbrace{1-b}})=\left(1-b\right)-\left(x-b\right)=1-x.$

Therefore, ℙ({U + b} ≤ x) = 1 − ℙ({U + b} ≥ x) = 1 − (1 − x) = x.

Finally, let us prove the last assertion. By induction, all Uk ~ Unif(0, 1), k = 0, 1, 2, . . . . Clearly, for any k ≥ 1, the pairs (Un, Un+k), n ≥ 0, all have the same distribution. Hence, we only need to prove that Corr(U0, Uk) = M−k for all k ≥ 0. Denote rk = Corr(U0, Uk). Let us show that $r_k=\frac{1}{M}{r}_{k-1}\text{for}k\ge 1$. Since r0 = Corr(U0, U0) = 1, the assertion will be proved by induction. A uniform random variable MU ~ Unif(0, M) can be expressed as a sum of its fractional and integer parts: $MU=\lfloor MU\rfloor +\left\{MU\right\}$. Clearly, $\lfloor MU\rfloor ~Unif\left\{0,1,...,M-1\right\}\left(\text{since}\mathbb{P}\left(\lfloor MU\rfloor \right)=k\right)=\mathbb{P}\left(MU\in [k,k+1)\right)=\frac{1}{M}$ for 0 ≤ k ≤ M − 1) and {MU} ~ Unif(0, 1) (as proved above). For any k = 0, 1, . . . , M − 1 and any x ∊ (0, 1), we have

$\mathbb{P}\left(\lfloor MU\rfloor =k;\left\{MU\right\}\le x\right)=\mathbb{P}\left(MU\in \left[k;k+x\right]\right)=\frac{x}{M}=\mathbb{P}\left(\lfloor MU\rfloor =k\right)\mathbb{P}\left(\left\{MU\right\}\le x\right).$

Therefore, the fractional and integer parts are independent random variables. It is also true that $\text{E}\left[\lfloor MU\rfloor \right]{=}^{\left(M-1\right)}/2$ and $\text{Var}\left(\lfloor MU\rfloor \right){=}^{\left(M^2-1\right)}/\sqrt{12}$. Moreover, we have

$\begin{array}{cc}\hfill \frac{U-\text{E}\left[U\right]}{\sqrt{\text{Var}\left(U\right)}}=& \frac{U-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{12}$}\right.}=\frac{MU-\raisebox{1ex}{$M$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\raisebox{1ex}{$M$}\!\left/ \!\raisebox{-1ex}{$\sqrt{12}$}\right.}\hfill \\ \hfill =& \frac{\lfloor MU\rfloor +\left\{MU\right\}-\raisebox{1ex}{$\left(M-1\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\raisebox{1ex}{$M$}\!\left/ \!\raisebox{-1ex}{$\sqrt{12}$}\right.}\hfill \\ \hfill =& \sqrt{\frac{M^2-1}{M^2}}\left(\frac{\lfloor MU\rfloor -\raisebox{1ex}{$\left(M-1\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\sqrt{\raisebox{1ex}{$\left(M^2-1\right)$}\!\left/ \!\raisebox{-1ex}{$12$}\right.}}\right)+\frac{1}{M}\left(\frac{\left\{MU\right\}-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{12}$}\right.}\right).\hfill \end{array}$

Now we are ready to calculate rk:

$r_k=\text{Corr}\left(U_0,U_k\right)=\sqrt{\frac{M^2-1}{M^2}}\text{Corr}\left(\lfloor M{U}_0\rfloor ,U_k\right)+\frac{1}{M}\text{Corr}\left(\left\{M{U}_0\right\},U_k\right).$

We proved that $\lfloor M{U}_0\rfloor$ and {MU0} = U1 are independent random variables. Moreover, since Uk = {MUk−1} = {M {MUk−2}} = · · · = g(U1) for some g : ℝ → ℝ, the random variable Uk as a function of U1 is independent of $\lfloor M{U}_0\rfloor$ as well. Thus, $\text{Corr}\left(\lfloor M{U}_0\rfloor ,U_k\right)=0$ and

$\begin{array}{ll}\hfill r_k=& \text{Corr}\left(U_0,U_k\right)=\frac{1}{M}\text{Corr}\left(\left\{M{U}_0\right\},U_k\right)=\frac{1}{M}\text{Corr}\left(U_1,U_k\right)\hfill \\ \hfill =& \frac{1}{M}\text{Corr}\left(U_0,U_{k-1}\right)=\frac{1}{M}{r}_{k-1}.\hfill \end{array}$

According to Proposition 17.5, a sequence {Ut, t = 0, 1, . . .} of random numbers uniformly distributed in (0, 1) can be generated from one random number U0 ~ Unif(0, 1) by using the rule Ut = {MUt−1} for t ≥ 1. We shall call this method a multiplicative method. Although the numbers obtained are dependent, the correlation between any two members of the sequence is negligible and rapidly goes to zero as the distance between the numbers in the sequence increases. Thus, the transformation y = {Mx + a} with a suitable choice of M and a can be used to generate pseudo-random numbers having good statistical properties such as uniformity and independence of samples. A drawback of the multiplicative method is that multidimensional tuples formed from the sequence {Ut}t≥0 lie on a family of multidimensional hyperplanes in the unit hypercube. For example, the points (U0, U1), (U2, U3), . . . lie on a family of M parallel lines.

The linear congruential generator (LCG) of PRNs is one of the oldest and simplest methods. It was proposed by Lehmer (1951). This algorithm is based on the mapping y = {Mx + a} but to avoid round-off errors the generating rule is written in a different form. As an example, let us consider the function y = {Mx} and suppose that x is a ratio of two integers: $x=\frac{s}{m}$. Then,

$\left\{M\frac{s}{m}\right\}=\frac{Ms}{m}-\lfloor \frac{Ms}{m}\rfloor =\frac{Ms-\lfloor \frac{Ms}{m}\rfloor m}{m}=\frac{Ms\mathrm{mod}m}{m}.$

The operation (s mod m) returns the remainder of an integer s after division by m:

$s\mathrm{mod}m:=s-\lfloor s/m\rfloor \cdot m$

This operation is also called the reduction of s by modulo m; the result is called the residue of s modulo m.

The LCG works as follows. First, a sequence {st, t = 0, 1, . . .} of integers is generated. The initial number s0 (called the initial seed) is chosen by the user and all subsequent numbers are calculated from the rule

$s_t=\left(M{s}_{t-1}+a\right)\mathrm{mod}m,t=1,2,...\text{ (17.1)}$

By construction, all integers st lie between 0 and m − 1. Then, the pseudo-random numbers in [0, 1) are obtained by dividing these integers by m:

$u_t=\frac{s_t}{m},t=0,1,2,...\text{ (17.2)}$

The LCG is equivalent to the recurrence uk = {Mut−1 + a}, where all ut are of the form $\frac{s}{m}$ with integers s and m, but it avoids round-off errors. Here the multiplier M, the increment a, and the modulus m are integers so that 1 ≤ M < m, and 0 ≤ a < m. The initial seed s0 is an integer from {0, 1, . . . , m − 1}. If a = 0, then the generator is called a multiplicative congruential generator (MCG). It operates on multiplicative group of integers modulo m. The generating rule of the MCG is

$s_t=M{s}_{t-1}\mathrm{mod}m,u_t=\frac{s_t}{m},t=1,2,...\text{ (17.3)}$

To guarantee that numbers produced from the rule (17.3) are all nonzero, the initial seed s0 has to be nonzero. Otherwise, all st = 0, t ≥ 1, if s0 = 0.

Example 17.1.

Construct PRNs generated by the LCG with M = 11, m = 16, a = 5, and s0 = 0.

Solution. The LCG rule is $u_t=\frac{s_t}{16}$, where st = (11st−1 + 5) mod 16. We first obtain the sequence of integers st, t ≥ 0:

$\begin{array}{c}{S}_1=\left(11\cdot 0+5\right)\mathrm{mod} 16=5\mathrm{mod} 16=5,s_2=\left(11\cdot 5+5\right)\mathrm{mod} 16=60\mathrm{mod} 16=12,\\ s_3=\left(11\cdot 12+5\right)\mathrm{mod} 16=137\mathrm{mod} 16=9,s_4=\left(11\cdot 9+5\right)\mathrm{mod} 16=104\mathrm{mod} 16=8,\\ s_5=\left(11\cdot 8+5\right)\mathrm{mod} 16=93\mathrm{mod} 16=13,s_6=\left(11\cdot 13+5\right)\mathrm{mod} 16=148\mathrm{mod} 16=4,\\ s_7=\left(11\cdot 4+5\right)\mathrm{mod} 16=49\mathrm{mod} 16=1,s_8=\left(11\cdot 1+5\right)\mathrm{mod} 16=16\mathrm{mod} 16=0,\\ s_9=\left(11\cdot 0+5\right)\mathrm{mod} 16=5\mathrm{mod} 16=5,...\end{array}$

As is seen, the sequence obtained is periodic. The numbers repeat themselves after eight steps. The PRNs $u_t=\frac{s_t}{m},t\ge 0$, are

$0,\frac{5}{16},\frac{12}{16},\frac{9}{16},\frac{8}{16},\frac{13}{16},\frac{4}{16},\frac{1}{16},0,\frac{5}{16},...$

The quality of the LCG depends on the choice of parameters M, a, and m. Often the modulus m is chosen as a prime number. Then, all calculations are done in the finite field Zm. Preferred moduli are Mersenne primes of the form 2r − 1, e.g., 231 − 1 = 2,147,483,647. Sometimes, m is a power of 2 since calculations can be done faster by exploiting the binary structure of computer arithmetic. Here are some choices of m and M for the MCG:

• m = 231 − 1 and M = 16807 (Park and Miller (1988));
• m = 240 and M = 517 (Ermakov and Mikhaĭlov (1982));
• m = 2128 and M = 5100109 mod 2128 (Dyadkin and Kenneth (2000)).

Since the set {0, 1, . . . , m − 1} is finite, the sequence generated by an LCG is periodic. That is, st+r = st holds for all t and some integer r > 0 called a period of the sequence. Let ℓ be the least possible period, which is called the length of period. Because there are only m possible different values of st, the maximum possible period of an LCG is m (or m − 1 if a = 0).

Proposition 17.6

(Ermakov (1975)). The maximum possible length of period for the sequence st = Mst−1 mod 2p, t ≥ 1, with p ≥ 3 is ℓ = 2n−2. The maximum period length is achieved if M ≡ 3 mod 8 or M ≡ 5 mod 8 and the seed s0 is odd.

There are many ways of improving a classical linear congruential generator. Let us consider some of them (a more comprehensive review of modern PRNGs can be found in L’Ecuyer (2012)). To increase the maximum period, one can combine two (or more) LCG methods as follows. Let {ut(1)}t≥1 and {ut(2)}t≥1 be the outputs of two LCGs. A new sequence {ut}t≥1 is given by

$u_t:=\left\{u_t^{\left(1\right)}+u_t^{\left(2\right)}\right\},t\ge 0.$

Another generalization of the classical LCG is a multiple recursive generator defined by:

$s_t=\left(M_1{s}_{t-1}+\cdots +M_k{s}_{t-k}\right)\mathrm{mod}m,u_t=\frac{s_t}{m},t\ge k,$

where M1, . . . , Mk are multipliers selected from $\mathcal{S}:=\left\{0,1,...,m-1\right\}$; k ≥ 2 is the order of recursion, and (s0, s1, . . . , sk−1) is the seed sequence with $s_i\in \mathcal{S}$ for 0 ≤ i ≤ k − 1. The maximum period length for this generator is mk − 1.

One of important properties of the MCG is the ease of its parallelization. Suppose that a sequence {st}t≥0 is generated by (17.3). To share this sequence among several independent processors, we split it into several disjoint subsequences. This can be achevied by using different seeds situated far apart along the original sequence. The seeds s0(j), j ≥ 0, that start respective subsequences of length K are generated by using the leaping-frog generator:

$s_0^{\left(j\right)}=A{s}_0^{\left(j-1\right)}\mathrm{mod}m,\text{where}{s}_0^{\left(0\right)}\equiv s_0\text{and}A\equiv M^K\mathrm{mod}m.$

Once the seeds are generated, the original generator (17.3) is used to obtain disjoint subsequences starting from these seeds:

$s_t^{\left(j\right)}=M{s}_{t-1}^{\left(j\right)}\mathrm{mod}m,u_t^{\left(j\right)}=\frac{s_t^{\left(j\right)}}{m},t=1,2,...,K,j\ge 0.$

The maximum number of processors that can be served by such a leaping-frog generator cannot exceed the ratio $\frac{\ell }{K}$ of the length of period to the length of an individual subsequence.

17.3 Generation of Nonuniformly Distributed Random Numbers

Suppose that we have a good PRN generator for the uniform distribution Unif(0, 1). However, our ultimate goal is to be able to sample from any given probability distribution. This can be achieved by transforming uniform PRNs into nonuniform random numbers. We are interested in general transformation methods that work for large classes of distributions. The transformation methods should be fast and efficient and should not use too much memory. In the sequel we are going to consider three main groups of sampling algorithms, namely, inversion methods, composition methods, and acceptance-rejection methods.

17.3.1 Transformations of Random Variables

Recall that a (univariate) random variable defined on a probability space (Ω, ℱ, P) is a measurable function from Ω to ℝ. The support of X is defined as the smallest closed set $\mathcal{S}$ whose compliment ${\mathcal{S}}^c$ has probability zero: $\mathbb{P}\left(X\in {\mathcal{S}}^c\right)=0$. The cumulative distribution function (CDF) of a random variable X is a function FX from ℝ to [0, 1] that is defined by FX(x) = ℙ(X ≤ x) for x ∊ ℝ. A CDF F satisfies the following properties:

1. F is a nondecreasing, right-continuous (i.e., F (x+) = F (x)) function;
2. F (−∞) = 0 and F (+∞) = 1.

A random variable X (and its CDF FX) is said to be

Discrete, if there exists a nonnegative function pX called a probability mass function (PMF for short) with a countable support such that

$F_X\left(x\right)={\displaystyle \sum _{y\le x:p_X\left(y\right)\ne 0}{p}_X\left(y\right)}\text{for}x\in ℝ.$

The CDF FX is a piecewise constant function with jumps.

(Absolutely) continuous, if there exists a nonnegative integrable function fX called a probability density function (PDF for short) such that

$F_X\left(x\right)={\displaystyle {\int }_{-\infty }^x{f}_X\left(x\right)}\text{d}x\text{for}x\in ℝ.$

The CDF FX is a continuous function (without jumps).

Since FX(∞) = 1, a PMF p satisfies ${\displaystyle {\sum }_{x\in \mathcal{S}}p\left(x\right)=1}$, and a PDF f ≥ 0 satisfies ${\displaystyle {\int }_{-\infty }^{\infty }f\left(x\right)\text{d}x=1}$. For a discrete random variable, the support is the smallest countable collection of points $\mathcal{S}$ so that $\mathbb{P}\left(X\in \mathcal{S}\right)=1$. It is defined by $\mathcal{S}=\left\{x\in ℝ:p_X\left(x\right)\ne 0\right\}$. Typically, the support $\mathcal{S}$ of a univariate continuous probability distribution is an interval of the real line. Since for a continuous random variable X the mass probability ℙ(X = x) is zero for any x ∊ ℝ, the support may be chosen to be an open interval.

Often, one random variable can be expressed as a function of another variate, say Y = f(X). Then sample values of Y can be obtained by applying the mapping f to sample values of X: yi = f(xi), i ≥ 1. Let us consider several useful examples.

A linear (or affine) mapping f(x) = α + βx. The density of Y = α + βX is given by $f_Y\left(x\right)=\frac{1}{\beta }{f}_X\left(\frac{x-\alpha }{\beta }\right)$. A normal random X ~ Norm(μ, σ2) can be expressed in terms of a standard normal variate Z ~ Norm(0, 1) as X = μ + σZ. Similarly, a uniform variable X ~ Unif(a, b) with a < b is a function of U ~ Unif(0, 1) given by X = a + (b − a)U.

A power mapping f(x) = xc for x > 0. Let X be a nonnegative random variable. The density of Y = Xc is then $f_Y\left(x\right)=\frac{1}{c}{x}^{\frac{1}{c}-1}f\left(x^{\frac{1}{c}}\right)$. For example, the density of Uc, where U ~ Unif(0, 1), is $f\left(x\right)=\frac{1}{c}{x}^{\frac{1}{c}-1}{\mathbb{I}}_{\left(0,1\right)}\left(x\right)$. A special case is a reciprocal of a nonzero variate: $Y=\frac{1}{X}$. The PDF is $f_Y\left(x\right)=\frac{1}{x^2}f\left(\frac{1}{x}\right)$.

An exponential mapping f(x) = exp(x). The PDF of Y = eX is $f_Y\left(x\right)=\frac{1}{x}f\left(\ln x\right)$ for x > 0. For example, the log-normal random variable is an exponential of a normal variate: Y = eμ+σZ, where Z ~ Norm(0, 1).

17.3.2 Inversion Method

The inversion method of generating nonuniform random numbers is based on the analytical or numerical inversion of the CDF. This method is most preferable if we wish to use it with low-discrepancy (quasi-random) sequences of numbers. Also it is compatible with some control variate techniques such as the antithetic variate method. The inversion method can be applied to both discrete and continuous distributions; however, its efficiency (in comparison with other methods) depends on how fast the inverse CDF can be evaluated.

17.3.2.1 Inverse Distribution Function

Consider a random variable X that has a continuous and strictly increasing on its support CDF F. Thus, the inverse function F−1 is well-defined and is also a strictly increasing function. Let us find the distribution function of Y := F(X). By using the definitions of a CDF and an inverse function, we obtain

$F_Y\left(y\right)=\mathbb{P}\left(F\left(X\right)\le y\right)=\mathbb{P}\left(F^{-1}\left(F\left(X\right)\right)\le F^{-1}\left(y\right)\right)=\mathbb{P}\left(X\le F^{-1}\left(y\right)\right)=F\left(F^{-1}\left(y\right)\right)=y$

for all y ∊ (0, 1). As is seen, the function FY is the CDF of a continuous random variable uniformly distributed on (0, 1). That is, F (X) ~ Unif(0, 1). This result provides us with a simple algorithm for sampling from continuous strictly increasing CDFs. Algorithm 17.1 is very simple and transparent. However, it relies on the knowledge of the inverse CDF F−1 in closed form (or the ability to compute F−1 in efficient manner). Figure 17.3 illustrates the method.

Example 17.2.

Using the inverse CDF method, find generating formulae for

1. (a) the uniform distribution Unif(a, b), a < b;
2. (b) the exponential distribution Exp(λ), λ > 0;
3. (c) the power distribution with the PDF $f\left(x\right)=c{x}^{c-1}{\mathbb{I}}_{\left(0,1\right)}\left(x\right),c>0$;
4. (d) the Weibull distribution with the PDF $f\left(x\right)=\frac{\alpha }{\beta }{x}^{\alpha -1}{\text{e}}^{-x^{\alpha }/\beta }{\mathbb{I}}_{{ℝ}_{+}}\left(x\right),\alpha ,\beta >0$.

Solution.

1. (a) The CDF of Unif(a, b) is $F\left(x\right)=\frac{x-a}{b-a}$ for x ∊ (a, b). Solve $\frac{x-a}{b-a}=u$ with u ∊ (0, 1) for x to find the inverse CDF: F−1(u) = a + (b − a)u. So the generating formula is X = a + (b − a)U.
2. (b) The CDF of the exponential distribution with rate λ > 0 is F(x) = 1 − e−x, x > 0. Solve 1 − e−λx = u for x to obtain $x=-\frac{1}{\lambda }\ln \left(1-u\right)$. If U ~ Unif(0, 1), then 1 − U ~ Unif(0, 1). So the generating formula for X ~ Exp(λ) simplifies: $X=-\frac{\ln U}{\lambda }$.

3. (c) Integrate $f\left(x\right)=c{x}^{c-1}{\mathbb{I}}_{\left(0,1\right)}\left(x\right)$ on (0, x) for 0 < x < 1, to obtain

   $F\left(x\right)={\displaystyle {\int }_0^{x}f\left(x\right)}\text{d}x=x^c.$

   Thus, F−1(u) = u1/c and X = U1/c.

4. (d) The CDF is

   $F\left(x\right)={\displaystyle {\int }_0^x\frac{\alpha }{\beta }{x}^{\alpha -1}{\text{e}}^{-x^{\alpha }/\beta }\text{d}x=1-}{\text{e}}^{-x^{\alpha }/\beta }\text{for}x>0.$

   Its inverse is F−1(u) = (−βln(1 − u))1/α. Thus, we obtain

   $X=F^{-1}\left(1-U\right)={\left(-\beta \ln U\right)}^{1/\alpha }.$

To generalize the inverse CDF sampling method to the case of any probability distribution, we define the generalized inverse CDF F−1 : (0; 1) → ℝ

$F^{-1}\left(u\right)=\inf \left\{x\in ℝ:u\le F\left(x\right)\right\}.\text{ (17.4)}$

To justify this formula, let us consider two special cases. First, suppose that a CDF F of a random variable X has a jump discontinuity at x0, i.e., F(x0−) < F(x0). Therefore, ℙ(X = x0) = F(x0) − F(x0−) ≠ 0. The formula in (17.4) gives that

$\begin{array}{l}{F}^{-1}\left(u\right)=x_0\text{for}F\left(x_0\text{−}\right)\le u\le F\left(x_0\right),\hfill \\ F^{-1}\left(u\right)<x_0\text{for}u<F\left(x_0\text{−}\right),\hfill \\ F^{-1}\left(u\right)>x_0\text{for}u>F\left(x_0\right).\hfill \end{array}$

Thus, the random variable F−1(U) has a nonzero mass probability at x0 and

$\mathbb{P}\left(F^{-1}\left(U\right)=x_0\right)=F\left(x_0\right)-F\left(x_0\text{−}\right)$

as expected. Now, suppose that F has a flat section on [x0, x1] with x0 < x1. There exists u0 ∊ [0, 1] such that F(x) = u0 for x ∊ (x0, x1), F (x) ≤ u0 for x ≤ x0, and F(x) ≥ u0 for x ≥ x1. In this case, we have ℙ(x0 < X < x1) = F (x1−)−F (x0) = 0. Then, the generalized inverse has a jump at u0: F−1(u0−) ≤ x0 and F−1(u0) ≥ x1. Thus, with probability zero the random variable F−1(U) takes on a value in (x0, x1) as expected. Now, as we can see, the same sampling formula X = F−1(U), U ~ Unif(0, 1), with a generalized inverse CDF in (17.4), works for any CDF F , including those having a jump discontinuity or a flat section.

Example 17.3.

Using the inverse CDF method, find a generating formula for

1. (a) a Bernoulli random variable $X=\{\begin{array}{cc}0& \text{with}\text{probability}1-p,\\ 1& \text{with}\text{probability}p;\end{array}$

2. (b) a random variable with the PMF

   $f\left(x\right)=p_1{\mathbb{I}}_{\left\{x_1\right\}}\left(x\right)+p_2{\mathbb{I}}_{\left\{x_2\right\}}\left(x\right)+p_3{\mathbb{I}}_{\left\{x_3\right\}}\left(x\right),$

   where pi, i = 1, 2, 3, are positive probabilities so that p1 + p2 + p3 = 1.

Solution.

1. (a) Sample U ~ Unif(0, 1). If U < p, then set X = 1; otherwise set X = 0. Verify that X has the Bernoulli distribution: ℙ(X = 1) = ℙ(U < p) = p and ℙ(X = 0) = ℙ(U ≥ p) = 1 − p.

2. (b) Suppose that x1 < x2 < x3. The CDF of X ~ f and the generalized inverse CDF are, respectively,

   $F_X\left(x\right)=\{\begin{array}{l}0\text{if}x<x_1,\hfill \\ p_1\text{if}{x}_1\le x<x_2,\hfill \\ p_1+p_2\text{if}{x}_2\le x<x_3,\hfill \\ 1\text{if}{x}_3\le x,\hfill \end{array}{F}_X^{-1}\left(u\right)=\{\begin{array}{l}{x}_1\text{if}u\le p_1,\hfill \\ x_2\text{if}{p}_1<u\le p_1+p_2,\hfill \\ x_3\text{if}{p}_1+p_2<u,\hfill \end{array}$

   for x ∊ ℝ and u ∊ (0, 1). Sample U ~ Unif(0, 1) and set

   $X=F_X^{-1}\left(U\right)=\{\begin{array}{l}{x}_1\text{if}U\le p_1,\hfill \\ x_3\text{if}U>p_1+p_2,\hfill \\ x_2\text{otherwise}.\hfill \end{array}$

Example 17.4.

Justify the following method of sampling from the discrete uniform distribution on a set of N distinct numbers x1, x2, . . . , xN:

1. (i) generate U ~ Unif(0, 1);
2. (ii) set X = xK, where $K=\lfloor N\cdot U+1\rfloor$.

Solution. We need to show that $\mathbb{P}\left(X=x_k\right)=\frac{1}{N}$ for any k = 1, 2, . . . , N. Indeed,

$\begin{array}{ll}\hfill \mathbb{P}(X=x_k)=& \mathbb{P}\left(\lfloor N\cdot U+1\rfloor =k\right)=\mathbb{P}\left(k\le N\cdot U+1<k+1\right)\hfill \\ \hfill =& \mathbb{P}\left(\frac{k-1}{N}\le U<\frac{k}{N}\right)=\frac{k}{N}-\frac{k-1}{N}=\frac{1}{N}.\hfill \end{array}$

17.3.2.2 The Chop-Down Search Method

Consider a general discrete random variable with a countable support $\mathcal{S}={\left\{x_j\right\}}_{j\ge 1}$ and mass probabilities {pj}j≥1, where pj = ℙ(X = xj) > 0 and Σj≥1pj = 1. The CDF F of such a discrete probability distribution is a piecewise-constant function. Let us assume that the mass points xj are sorted in increasing order: x1 < x2 < . . .. Then the CDF is given by $F\left(x\right)={\displaystyle \sum _{j:x_j\le x}{p}_j}$. Hence the generalized inverse CDF is calculated as follows:

$F^{-1}\left(u\right)=\inf \left\{x_k\in S:u\le {\displaystyle \sum _{j=1}^k{p}_j}\right\}=\left\{x_k\in \mathcal{S}:{\displaystyle \sum _{j=1}^{k-1}{p}_j}<u\le {\displaystyle \sum _{j=1}^k{p}_j}\right\}\text{ (17.5)}$

The requirement that mass points xj are sorted in increasing order is not necessary for the application of the inversion method. We can consider any arrangement for {xj}j≥1, since the sampling of X is equivalent to the sampling of a random index $K\in \mathbb{N}$ with probabilities ℙ(K = j) = pj, j ≥ 1. First, generate U ~ Unif(0, 1). Second, find the index K ≥ 1 such that

${\displaystyle \sum _{j=1}^{K-1}{p}_j}<U\le {\displaystyle \sum _{j=1}^K{p}_j}.\text{ (17.6)}$

Finally, set X = xK. Note that the probability of the event that U satisfies the above double inequality is exactly pK. One of possible implementations of this approach is the chop-down search (CDS) algorithm (see Algorithm 17.2).

The number of cycles in the chop-down search method is equal to the expected value ${\displaystyle {\sum }_{j\ge 1}j{p}_j}$. Indeed, if U ∊ (0, p1], then the algorithm stops after one cycle, and this happens with probability p1 = ℙ(U ∊ (0, p1]). If U ∊ (p1, p1 + p2], then the algorithm stops after two cycles, and this happens with probability p2 = ℙ(U ∊ (p1, p1 +p2]), and so on. Let cU be the computational cost of the generation of one draw from Unif(0, 1) and cI the computational cost of one cycle of the method. Then, the total cost is

$c_U+c_I{\displaystyle \sum _{j\ge 1}j{p}_j}.$

Example 17.5.

Find the computational cost of the CDS method for

1. (a) the geometric distribution Geom(p) with ℙ(X = j) = (1 − p)j−1p, j ≥ 1, 0 < p < 1;
2. (b) the Poisson distribution Pois(λ) with $\mathbb{P}\left(X=j\right)=\frac{{\lambda }^j}{j!}{\text{e}}^{-\lambda },j\ge 0,\lambda >0$.

Solution. Compute the mathematical expectation $\epsilon ={\displaystyle {\sum }_{j=1}^{\infty }j{p}_j}$ for both distributions. The total cost is then cU + CIε.

$\left(\text{a}\right)\epsilon ={\displaystyle \sum _{j=1}^{\infty }j{\left(1-p\right)}^{j-1}p=\frac{1}{p};\left(\text{b}\right)}\epsilon ={\displaystyle \sum _{j=0}^{\infty }\left(j+1\right)}\frac{{\lambda }^j}{j!}{\text{e}}^{-\lambda }=\lambda +1.$

The computational cost can be reduced by rearranging elements of {(xj, pj)}j≥1. Consider an arrangement {ji}i≥1 of the integers 1, 2, 3, . . .. Then, the sequence ${\left\{\left(x_{j_i},p_{j_i}\right)\right\}}_{i\ge 1}$ defines another discrete probability distribution equivalent to the original one.

Proposition 17.7.

The computational cost of the chop-down search method attains its minimal value iff the mass probabilities are arranged in the decreasing order

$p_1\ge p_2\ge p_3\ge \cdots \text{ (17.7)}$

Proof. Suppose that there exists another arrangement of the mass probabilities, {pj′}j≥1, for which (17.7) is violated and the expected value ${\epsilon }^{\prime }={\displaystyle {\sum }_{j\ge 1}j{p}_j^{\prime }}$ is minimal. There are two indices k and l, k < l, so that pl′ < pk′. Let us construct another arrangement ${\left\{{p^{″}}_j\right\}}_{j\ge 1}$, which is obtained from {pj′}j≥1 by swapping pl′ and pk′, i.e., pl″ = pk′, pk″ = pl′, and pj″ = pj′ for $j\notin \left\{k,l\right\}$. Let ${\epsilon }^{″}={\displaystyle {\sum }_{j\ge 1}j{p^{″}}_j}$. We have

${\epsilon }^{\prime }-{\epsilon }^{″}=l{p}_l^{\prime }+k{p}_k^{\prime }-l{p}_l^{″}-k{p}_k^{″}=\left(l-k\right)\left(p_l^{\prime }-p_k^{\prime }\right)>0,$

since l < k and pl′ < pk′. Hence ε″ < ε′. We arrive at a contradiction.

According to Proposition 17.7, the mass probabilities {pj}j≥1 should be arranged in decreasing order before applying the chop-down search method. Another way to speed up calculations is to use recurrence relations for mass probabilities. This allows us to reduce the parameter cI—the cost of one cycle of the CDS method. For example, the mass probabilities of a geometric random variable X ~ Geom(p) satisfy

$\mathbb{P}\left(X=j+1\right)=\mathbb{P}\left(X=j\right)\cdot \left(1-p\right),j=1,2,...$

For a Poisson random variable X ~ Pois(λ), we have

$\mathbb{P}\left(X=j\right)=\mathbb{P}\left(X=j-1\right)\cdot \frac{\lambda }{j},j=1,2,...$

17.3.2.3 The Binomial Search Method

The formula (17.5) of the generalized inverse CDF requires the cumulative probabilities ${\displaystyle {\sum }_{j=1}^k{p}_j,k\ge 1}$. Therefore, the inversion method can be speeded up if these cumulative probabilities are precalculated in advance and stored in memory. After that, to find K that satisfy (17.6), we can employ a fast search procedure such as the binary search method. Suppose that a random variable X takes on one of N distinct values x1 < x2 < · · · < xN with the respective mass probabilities p1, p2, . . . , pN , i.e., ℙ(X = xj) = pj, 1 ≤ j ≤ N. Note that a random variable with a countably infinite support can be reduced to the case with N mass points by suitable truncation of the support such that the total probability of removed mass points is very small. Suppose that N = 2r, then r = log2 N cycles are required to find K. For general N, the computation cost is proportional to $\lceil {\log }_{2}N\rceil$, where $\lceil \cdot \rceil$ denotes the ceiling function.

17.3.3 Composition Methods

It is a well-known fact that a linear combination (a mixture) of CDFs, PDFs, or PMFs with positive weights summing up to one is again a CDF, PDF, or PMF, respectively. The composition method aims to represent the probability distribution of interest as a mixture of simpler-organized distributions. The sampling from a mixture distribution is a two-step procedure. First, one of the distributions which appear in the composition is selected at random; second, a sample is drawn from the distribution selected (e.g., by using an inversion algorithm). In comparison with the inverse CDF method that requires one draw from Unif(0, 1), a composition method needs at least two uniform random numbers. However, the composition method allows us to express a probability distribution in a simpler form that allows for simplifying the sampling algorithm.

17.3.3.1 Mixture of PDFs

Consider a continuous random variable X with PDF f. Suppose that f can be represented as a linear combination of m PDFs f1, . . . , fm with m positive weights w1, . . . , wm so that w1 + · · · + wm = 1:

$f\left(x\right)={\displaystyle \sum _{j=1}^m{w}_j{f}_j}\left(x\right),x\in ℝ.$

The PDF f is called a mixture PDF. The support of f is a union of supports of fj, 1 ≤ j ≤ m. If all pairwise intersections of supports of fj are empty sets, then such a mixture is called a stratification.

Proof of Algorithm 17.4. Let us find the distribution function of X generated by the algorithm. Applying the total probability law gives

$\begin{array}{ll}\hfill \mathbb{P}\left(X\le x\right)=& {\displaystyle \sum _{j=1}^m\mathbb{P}}\left(X\le x;K=j\right)={\displaystyle \sum _{j=1}^m\mathbb{P}}\left(K=j\right)\mathbb{P}\left(X\le x|K=j\right)={\displaystyle \sum _{j=1}^m{w}_j}{\displaystyle {\int }_{-\infty }^x{f}_j\left(x\right)}\text{d}x\hfill \\ \hfill =& {\displaystyle {\int }_{-\infty }^x{\displaystyle \sum _{j=1}^m{w}_j}{f}_j\left(x\right)}\text{d}x={\displaystyle {\int }_{-\infty }^{x}f\left(x\right)\text{d}x=F_X\left(x\right)}\hfill \end{array}$

for all x ∊ ℝ.

Example 17.6.

Develop a stratification method for sampling from the PDF

$f\left(x\right)=\{\begin{array}{l}\frac{2}{3}x\text{if}0<x\le 1,\hfill \\ \frac{2}{3}\text{if}1<x<2,\hfill \\ 0\text{otherwise}\text{.}\hfill \end{array}$

Solution. As is seen, f is a piecewise linear function; it is constant on [1, 2] and linear on [0, 1]. Introduce two new PDFs: f1(x) ∝ for x ∊ (0, 1] and f2(x) ∝ 1 for x ∊ (1, 2). The second function is a PDF for Unif(1, 2). Hence, $f_2\left(x\right)={\mathbb{I}}_{\left(1,2\right)}\left(x\right)$. Calculate a normalizing constant for f1 to obtain

$f_1\left(x\right)=\frac{x}{{\displaystyle {\int }_0^{1}x\text{d}x}}{\mathbb{I}}_{(0,1]}\left(x\right)=2x{\mathbb{I}}_{(0,1]}\left(x\right).$

Now, the PDF f(x) can be decomposed as follows:

$\begin{array}{ll}\hfill f\left(x\right)=& \frac{2}{3}x{\mathbb{I}}_{(0,1]}\left(x\right)+\frac{2}{3}{\mathbb{I}}_{\left(1,2\right)}\left(x\right)=\frac{1}{3}\cdot \left(2x{\mathbb{I}}_{(0,1]}\left(x\right)\right)+\frac{2}{3}\cdot {\mathbb{I}}_{\left(1,2\right)}\left(x\right)\hfill \\ \hfill =& \frac{1}{3}{f}_1\left(x\right)+\frac{2}{3}{f}_2\left(x\right):=w_1{f}_1\left(x\right)+w_2{f}_2\left(x\right).\hfill \end{array}$

Sampling from f2 is easy (see Example 17.3): X = 1 + U with U ~ Unif(0, 1). To sample from f1, we apply the inverse CDF method. First, find the CDF:

$F_1\left(x\right)={\displaystyle {\int }_0^x{f}_1\left(x\right)\text{d}x=x^2,x\in \left[0,1\right].}$

The inverse CDF is $F_1^{-1}\left(u\right)=\sqrt{u},u\in \left[0,1\right]$. As a result, we obtain the following sampling algorithm:

1. (1) Generate two independent uniform random numbers U1, U2 ← Unif(0, 1).
2. (2) Sample K ∊ {1, 2} as follows. If $U_1<\frac{1}{3}$, then K = 1, else K = 2.

3. (3) Generate X from fK as follows:

   1. (i) if K = 1, then set X ← 1 + U2;
   2. (ii) if K = 2, then set $X\leftarrow \sqrt{U_2}$.

In general, the stratification method can be used with any partition of the support of a probability distribution. Suppose that A is the support of a PDF f. Let

$A={\displaystyle \underset{j=1}{\overset{m}{\cup }}{A}_k,\text{where}}i\ne j\Rightarrow A_i\cap A_j=\varphi .$

Then, the PDF f admits the following representation:

$f\left(x\right)=f\left(x\right){\mathbb{I}}_A\left(x\right)=f\left(x\right){\displaystyle \sum _{j=1}^m{\mathbb{I}}_{A_j}}\left(x\right)={\displaystyle \sum _{j=1}^m{w}_j{f}_j\left(x\right)},$

where $w_j={\displaystyle {\int }_{A_j}f\left(x\right)\text{d}x}$ and $f_j\left(x\right)=\frac{1}{{\omega }_j}f\left(x\right){\mathbb{I}}_{A_j}\left(x\right),1\le j\le m,x\in ℝ$.

17.3.3.2 Randomized Gamma Distributions

Probability distributions whose PDFs contain special functions such as Bessel and other hypergeometric functions present a real challenge for sampling random numbers. By replacing a special function by its integral or series representation in terms of simpler functions, the PDF can be expressed as a mixture of more regular densities with known sampling algorithms.

For example, consider the noncentral chi-square distribution with κ > 0 degrees of freedom and noncentrality parameter λ > 0. Its PDF is

$f\left(x;\kappa ,\lambda \right)=\frac{1}{2}{\text{e}}^{-\left(x+\lambda \right)/2}{\left(\frac{x}{\lambda }\right)}^{\frac{\kappa }{4}-\frac{1}{2}}{I}_{\frac{\kappa }{2}-1}\left(\sqrt{\lambda x}\right),x>0.\text{ (17.8)}$

The modified Bessel function Iμ of the first kind (of order μ) admits the following series expansion:

$I_{\mu }\left(x\right)={\left(\frac{x}{2}\right)}^{\mu }{\displaystyle \sum _{j=0}^{\infty }\frac{{\left(x^2/4\right)}^j}{j!\text{Γ}\left(\mu +j+1\right)}.}$

By using this expansion for $I_{\frac{\kappa }{2}-1}$ in (17.8), we can represent the noncentral chi-square PDF as a mixture of gamma densities with Poisson weights:

$\begin{array}{l}f\left(x;\kappa ,\lambda \right)=\frac{1}{2}{\text{e}}^{-\left(x+\lambda \right)/2}{\left(\frac{x}{\lambda }\right)}^{\frac{\kappa }{4}-\frac{1}{2}}{\left(\frac{\sqrt{\lambda x}}{2}\right)}^{\frac{\kappa }{2}-1}{\displaystyle \sum _{j=0}^{\infty }\frac{{\left(\lambda x/4\right)}^j}{j!\text{Γ}\left(\frac{\kappa }{2}+j\right)}}\hfill \\ ={\displaystyle \sum _{j=0}^{\infty }\underset{=p_j\left(\text{poisson}\text{prob}.\right)}{\underbrace{{\text{e}}^{-\lambda /2\frac{{\left(\lambda /2\right)}^j}{j!}}}}\underset{=f_j\left(x\right)\left(\text{gamma}\text{density}\right)}{\underbrace{\frac{1}{2}{\left(\frac{x}{2}\right)}^{\kappa /2+j-1}\frac{{\text{e}}^{-x/2}}{\text{Γ}\left(\frac{\kappa }{2}+j\right)}.}}}\hfill \end{array}\text{ (17.9)}$

Recall that a random variable X is said to be gamma-distributed with shape parameter α and scale parameter θ, denoted by Gamma(α, θ), if its PDF is

$f_X\left(x\right)=\frac{\theta }{\text{Γ}\left(\alpha \right)}{\left(\theta x\right)}^{\alpha -1}{\text{e}}^{-\theta x},x>0.\text{ (17.10)}$

In (17.9), the PDF f(x; κ, λ) is expressed as a mixture ${\displaystyle {\sum }_{j=0}^{\infty }{p}_j{f}_j\left(x\right)}$, where {pj}j≥0 are mass probabilities of the Poisson distribution with intensity $\frac{\lambda }{2}$ and fj is a gamma PDF of the form (17.10) with parameters $\alpha =\frac{\kappa }{2}+j$ and $\theta =\frac{1}{2}$ for all j = 0, 1, 2, . . . . Such a mixture probability distribution is called a randomized gamma distribution, denoted $Gamma\left(Y_1+\frac{\kappa }{2},\frac{1}{2}\right)$, where $Y_1~Pois\left(\frac{\lambda }{2}\right)$.

In general, we can consider a mixture gamma distribution Gamma(Y + ν + 1, θ) with parameters ν > −1 and θ > 0, where the randomizer Y is a discrete random variable taking its values in the set of nonnegative integers with probabilities ℙ(Y = j) = pj, j = 0, 1, 2, . . . . The PDF f of such a randomized gamma distribution admits the form of a series expansion:

$f\left(y\right)={\displaystyle \sum _{j=0}^{\infty }{p}_j}\frac{\theta }{\text{Γ}\left(v+j+1\right)}{\left(\theta y\right)}^{v+j}{\text{e}}^{-\theta \text{y}}.$

Let us consider three choices for the randomizer Y. The resulting probability distributions are called the randomized gamma distribution of the first, second, and third types, respectively.

Let Y1 ~ Pois(α) be a Poisson random variable with mean α > 0. The randomized gamma distribution of the first type is Gamma(Y1 + ν + 1, θ) with the PDF

$f_1\left(y\right)=\theta {\left(\frac{\theta y}{\alpha }\right)}^{v/2}{\text{e}}^{-\alpha -\theta y}{I}_v\left(2\sqrt{\alpha \theta y}\right),y>0.\text{ (17.11)}$

So the noncentral chi-square distribution with parameters κ > 0 and λ > 0 is the randomized gamma distribution of the first type with $\nu =\frac{\kappa }{2}-1,\theta =\frac{1}{2}$, and $\alpha =\frac{\lambda }{2}$.

A discrete random variable Y2 is said to have the Bessel probability distribution, denoted Bes(ν, b), with parameters ν > −1 and b > 0 if

$\mathbb{P}\left(Y_2=j\right)=\frac{{\left(b/2\right)}^{2j+v}}{I_v\left(b\right)j!\text{Γ}\left(j+v+1\right)},j=0,1,2,...\text{ (17.12)}$

This distribution is related to many other distributions, where the modified Bessel function I is involved in the density, including the squared Bessel bridge distribution. The randomized gamma distribution of the second type is a mixture distribution Gamma(Y1 + 2Y2 + ν + 1; θ), where $Y_1~Pois\left(\frac{a+b}{4\theta }\right)$ and $Y_2~Bes\left(\nu ,\frac{\sqrt{ab}}{2\theta }\right)$ are independent Poisson and Bessel random variables, respectively. For any positive numbers θ, a, b, and ν > −1, the PDF is

$f_2\left(y\right)=\theta {\text{e}}^{-\theta y-\left(a+b\right)/\left(4\theta \right)}\frac{I_v\left(\sqrt{ay}\right)I_v\left(\sqrt{by}\right)}{I_v\left(\sqrt{ab}/\left(2\theta \right)\right)},y>0.\text{ (17.13)}$

A discrete random variate Y3 is said to follow an incomplete gamma probability distribution, which we simply denote by IΓ(ν, λ) with parameters λ > 0 and ν > 0, if

$\mathbb{P}\left(Y_3=j\right)=\frac{{\text{e}}^{-\lambda }{\lambda }^{j+v}}{\text{Γ}\left(j+v+1\right)}\frac{\text{Γ}\left(v\right)}{\gamma \left(v,\lambda \right)},j=0,1,2,...,\text{ (17.14)}$

where $\gamma \left(a,x\right):={\displaystyle {\int }_0^x{t}^{a-1}}\text{e}-t\text{d}t$ is the lower incomplete gamma function. Note that if ν is a nonnegative integer, then the distribution of Y3 is simply a truncated and shifted Poisson distribution thanks to the property

$\frac{\gamma \left(m,a\right)}{\text{Γ}\left(m\right)}=1-\left(1+x+\cdots +\frac{x^{m-1}}{\left(m-1\right)!}\right){\text{e}}^{-x},m=0,1,2,...$

We call a mixture probability distribution Gamma(Y3 + 1, θ), Y3 ~ IΓ(ν, λ), the randomized gamma distribution of the third type. The PDF is

$f_3\left(y\right)=\frac{\theta \text{Γ}\left(v\right)}{\gamma \left(v,\lambda \right)}{\left(\frac{\theta y}{\lambda }\right)}^{-v/2}{\text{e}}^{-\lambda -\theta y}{I}_v\left(\sqrt{4\lambda \theta y}\right),y>0.\text{ (17.15)}$

As we will see in the next chapter, randomized gamma distributions play a significant role in simulation of the so-called constant elasticity of variance diffusion model (see also Makarov and Glew (2010)).

17.3.3.3 The Alias Method by Walker

Let us study the case with discrete random variables. A mixture of PMFs is defined in the same manner as a mixture of PDFs. Consider m PMFs pj(x) and m weights wj > 0, 1 ≤ j ≤ m, such that w1 + w2 + · · · + wm = 1. The function p defined by $p\left(x\right)={\displaystyle {\sum }_{j=1}^m{w}_j{p}_j\left(x\right)}$, x ∊ ℝ, is also a PMF called the mixture of the PMFs pj, 1 ≤ j ≤ m. To sample from p, we can apply Algorithm 17.4, where the sampling from the PDF fK is replaced by sampling from the PMF pK.

The alias method proposed by Walker (1977) allows us to represent any discrete probability distribution with m mass points $\mathcal{S}:=\left\{x_1,x_2,...,x_m\right\}$ as an equally weighted mixture of m two-point distributions. That is, there exist m two-point PMFs

$p_j\left(x\right)=q_j{\mathbb{I}}_{\left\{x_j\right\}}\left(x\right)+\left(1-q_j\right){\mathbb{I}}_{\left\{a_j\right\}}\left(x\right)\text{with}{a}_j\in \mathcal{S}\text{and}{q}_j\in \left[0,1\right],1\le j\le m,$

such that

$p\left(x\right)=\frac{1}{m}{\displaystyle \sum _{j=1}^m{p}_j\left(x\right)}\text{for}x\in ℝ.\text{ (17.16)}$

Since all weights in (17.16) are equal to $\frac{1}{m}$, Algorithm 17.4 is simplified and we obtain Algorithm 17.5.

To obtain such a decomposition of the PMF p, we need to construct two lists, namely, the list of probabilities Q = (q1, q2, . . . , qm) and the list of aliases A = (a1, a2, . . . , am). This can be achieved by using the “leveling the histogram” procedure, which is described below. During this procedure the original histogram {wj = pj : 1 ≤ j ≤ m} is transformed into an equally weighed histogram $\left\{w_j=\frac{1}{m}:1\le j\le m\right\}$; the lists Q and A are generated in the course of this process.

Step 1: Start with wj = pj, qj = 1, and aj = xj for all j = 1, 2, . . . , m.

Step 2: Find two indices ℓ, u ∊ {1, 2, . . . , m} (ℓ and u stand for “lower” and “upper,” respectively) such that

$w_{\ell }=\underset{1\le j\le m}{\min }\left\{w_j\right\}\text{and}{w}_u=\underset{1\le j\le m}{\max }\left\{w_j\right\}.$

Step 3: If $w_{\ell }=w_u=\frac{1}{m}$, then the histogram is levelled and the algorithm is stopped; otherwise (i.e., $w_{\ell }<\frac{1}{m}<w_u$) we proceed with Step 4

Step 4: Set qℓ ← N wℓ and aℓ ← xu. Change $w_u\leftarrow w_u-\left(\frac{1}{m}-w_{\ell }\right)$ and $w_{\ell }\leftarrow \frac{1}{m}$. Note that wu can become less than $\frac{1}{m}$. Go back to Step 2.

As is seen, the alias method requires at most m − 1 iterations until all columns of the histogram (i.e., the weights wj, 1 ≤ j ≤ m) become the same height.

Example 17.7.

Apply the alias method to the probability distribution with

$p\left(x\right)=0.1{\mathbb{I}}_{\left\{1\right\}}\left(x\right)+0.2{\mathbb{I}}_{\left\{2\right\}}\left(x\right)+0.3{\mathbb{I}}_{\left\{3\right\}}\left(x\right)+0.4{\mathbb{I}}_{\left\{4\right\}}\left(x\right).$

Solution. This is a four-point distribution with the list of mass points (x1, x2, x3, x4) = (1, 2, 3, 4) and the list of mass probabilities (p1, p2, p3, p4) = (0.1, 0.2, 0.3, 0.4). Initially, we set (a1, a2, a3, a4) = (1, 2, 3, 4) and (w1, w2, w3, w4) = (0.1, 0.2, 0.3, 0.4); all qi are equal to 1.

1. (i) The smallest weight is w1 = 0.1; the largest is w4 = 0.4. Set

   $\begin{array}{l}{w}_4\leftarrow w_4-\left(0.25-w_1\right)=0.4-\left(0.25-0.1\right)=0.25,q_1\leftarrow 4w_1=0.4,\hfill \\ w_1\leftarrow 0.25,a_1\leftarrow x_4=4.\hfill \end{array}$

2. (ii) Now, (w1, w2, w3, w4) = (0.25, 0.2, 0.3, 0.25). The smallest weight is w2 = 0.2; the largest is w3 = 0.3. Set

   $\begin{array}{l}{w}_3\leftarrow w_3-\left(0.25-w_1\right)=0.3-\left(0.25-0.2\right)=0.25,q_2\leftarrow 4w_2=0.8,\hfill \\ w_2\leftarrow 0.25,a_1\leftarrow x_3=3.\hfill \end{array}$

3. (iii) Now, (w1, w2, w3, w4) = (0.25, 0.25, 0.25, 0.25). All weights are equal to 0.25. Stop the algorithm.

As a result, we obtained the following two lists:

$\left(a_1,a_2,a_3,a_4\right)=\left(4,3,3,4\right)\text{and}\left(q_1,q_2,q_3,q_4\right)=\left(0.4,0.8,1,1\right).$

17.3.4 Acceptance-Rejection Methods

To generate realizations from some probability distribution, an acceptance-rejection method makes use of realizations of another random variable whose probability distribution is similar to the target one. The distribution from which the independent samples are generated is called a proposal distribution. Each sample can be accepted or rejected. The realizations being accepted have the target probability distribution. The computational cost of such a method is proportional to the average number of draws generated before one is accepted. Since the number of trials (draws) before the first success (acceptance of a draw) follows the geometric distribution, the average number of draws generated before an acceptance occurs is a reciprocal of the probability of the acceptance of a proposal draw.

Example 17.8.

Propose an acceptance-rejection method of sampling a point uniformly distributed in the circle C = {(x, y) : x2 + y2 ≤ 1}.

Solution. The circle C is contained in the square S = [−1, 1]2. To obtain a draw from the uniform distribution Unif(C), proceed as follows.

1. (1) Sample a random point (X, Y) uniformly distributed in S as follows: X = 2U1 − 1 and Y = 2U2 − 1, where U1 and U2 are i.i.d. Unif(0, 1)-distributed random variables.
2. (2) Accept the point if (X, Y) ∊ C, i.e., X2 + Y2 ≤ 1. Otherwise, the point is rejected and we return to (1).

The formal justification of this example and of the general acceptance-rejection method is based on Propositions 17.8–17.10, which follow.

Proposition 17.8.

Let a random vector X be uniformly distributed in a domain D ⊂ ℝd of a finite (d-dimensional) volume |D| < ∞. Let Ω be a subdomain of D. The distribution of X conditional on X ∊ Ω is uniform in Ω.

Proof. Let dΩ be an arbitrary subdomain of Ω. Then,

$\mathbb{P}\left(X\in \text{d}\Omega |X\in \Omega \right)=\frac{\mathbb{P}\left(X\in \text{d}\Omega ;X\in \Omega \right)}{\mathbb{P}\left(X\in \Omega \right)}=\frac{\mathbb{P}\left(X\in \text{d}\Omega \right)}{\mathbb{P}\left(X\in \Omega \right)}=\frac{|\text{d}\Omega |/|D|}{|\Omega |/|D|}=\frac{|\text{d}\Omega |}{|\Omega |}.$

Thus, the assertion is proved.

In Example 17.8 we deal with a case covered by Proposition 17.8. Proposed points are sampled uniformly on the square S. A point (X, Y) is accepted if it lies in the circle C. According to Proposition 17.8, the probability distribution of (X, Y) conditional on (X, Y) ∊ C is uniform in C. As a result, accepted points are uniformly distributed in C.

The simplest acceptance-rejection algorithm is the so-called Neumann method. Consider a bounded PDF f ≤ M with a support contained in a finite interval [a, b]. The plot of f on [a, b] is contained in the rectangle [a, b] × [0, M]. To sample from f we proceed as follows. First, sample (X, Y) uniformly in the rectangle [a, b] × [0, M]. This point is accepted if Y ≤ f(X) and rejected otherwise. According to Proposition 17.8, the distribution of accepted points is uniform in the region bounded by the plot of y = f(x) and the x-axis. Moreover, we can show that the distribution of the x-coordinate of an accepted point has the PDF f. That is, X conditional on Y ≤ f(X) has the distribution with the PDF f. This result will be proved in Proposition 17.9 in more general setting.

Definition 17.1.

Consider a nonnegative integrable function g with support D ⊆ ℝd, that is, g(x) ≥ 0 for x ∊ D, g(x) = 0 for x ∉ D, and $I_D\left(g\right):={\displaystyle {\int }_{{ℝ}^n}g\left(x\right)dx}={\displaystyle {\int }_{D}g\left(x\right)dx<\infty }$ hold. The region

$B_D\left(g\right):=\left\{\left(\text{x},y\right):\text{x}\in D,0<y<g\left(\text{x}\right)\right\}\subseteq {ℝ}^{d+1}$

is called the body of the function g.

Proposition 17.9.

Suppose that the (d + 1)-dimensional point (X, Y), with X ∊ ℝd and Y ∊ ℝ, is uniformly distributed in the body BD(g) of an integrable function g : D → [0, ∞) defined on D ⊆ ℝd. Then, the random vector X is distributed in D with the PDF proportional to g:

$p\text{x}\left(\text{x}\right)=\frac{1}{I_D\left(g\right)}g\left(\text{x}\right),\text{x}\in D.$

Proof. The joint PDF of (X, Y) is

$f_{\text{X},\text{Y}}\left(\text{x},y\right)=\frac{1}{|B_D\left(g\right)|}{\mathbb{I}}_{B_D\left(g\right)}\left(\text{x},y\right)=\frac{1}{|B_D\left(g\right)|}{\mathbb{I}}_D\left(\text{x}\right){\mathbb{I}}_{\left(0,g\left(\text{x}\right)\right)}\left(y\right).$

The marginal PDF of X is

$p_{\text{X}}\left(\text{x}\right)={\displaystyle {\int }_{-\infty }^{\infty }{f}_{\text{X},Y}\left(\text{x,}y\right)}\text{d}y={\displaystyle {\int }_0^{g\left(\text{x}\right)}\frac{1}{|B_D\left(g\right)|}{\mathbb{I}}_D\left(\text{x}\right)}\text{d}y=\frac{1}{I_D\left(g\right)}g\left(\text{x}\right){\mathbb{I}}_D\left(\text{x}\right),$

since ID(g) = |BD(g)|.

In the Neumann method the proposed sample values are drawn from the uniform PDF $p\left(x\right)=\frac{1}{b-a}{\mathbb{I}}_{\left(a,b\right)}\left(x\right)$. This choice is explained by the fact that the PDF p majorizes the target PDF f up to a multiplicative constant: f(x) ≤ C p(x) with C = M(b − a) (provided that f(x) ≤ M for all x ∊ [a, b]). A draw X is accepted if Y ≤ f(X) where Y ~ Unif(0, M). Therefore, the Neumann method can be generalized to the case with an arbitrary PDF f as long as we can find a majorizing function for f.

In what follows, we will require the following proposition that explains how to sample a point uniformly distributed in a body of a nonnegative function.

Proposition 17.10.

Let p be a multivariate PDF with support D. Suppose that X ~ p and (Y |X = x) ~ Unif(0, C p(x)) for some constant C > 0. Then, the point (X, Y) is uniformly distributed in the domain BD(C p).

Proof. The joint PDF of (X, Y) is

$\begin{array}{c}{f}_{\text{X,}Y}\left(\text{x,y}\right)=f_{\text{X}}\left(\text{x}\right)f_{Y|\text{X}}\left(y|\text{x}\right)=p\left(\text{x}\right)\frac{1}{Cp\left(\text{x}\right)}{\mathbb{I}}_{\left(0,Cp\left(\text{x}\right)\right)}\left(y\right)\\ =\frac{1}{C}{\mathbb{I}}_D\left(\text{x}\right){\mathbb{I}}_{\left(0,Cp\left(\text{x}\right)\right)}\left(y\right)=\frac{1}{|B_D\left(Cp\right)|}{\mathbb{I}}_{B_D\left(Cp\right)}\left(\text{x}\right),\end{array}$

since $|B_D\left(C_P\right)|={\displaystyle {\int }_D{C}_p\left(x\right)dx=C{\displaystyle {\int }_{D}p\left(x\right)dx=C.}}$.

Consider an n-variate PDF f with support D ⊆ ℝn. Suppose that there exists another PDF p called a proposal PDF and a constant C > 0 such that f(x) ≤ C p(x) for all x ∊ ℝn. Often p is chosen to be a simple function like a piecewise-linear function so that sampling from p is feasible.

The acceptance-rejection method can be visualized as choosing a subsequence of draws from a sequence of i.i.d. realizations from the PDF p in such a way that the resulting subsequence consists of i.i.d. realizations from the target PDF f.

The proof of the acceptance-rejection method is based on Propositions 17.8–17.10. The steps of Algorithm 17.6 can be reformulated as follows. First, sample two random variables X ~ p and Y ~ Unif(0, C p(X)). As is proved in Proposition 17.10, the point (X, Y) is uniformly distributed in BD(C p). If (X, Y) ∊ BD(f), then this point is accepted. In accordance with Proposition 17.8, the accepted point is uniformly distributed in BD(f). Finally, as follows from Proposition 17.9, the X-coordinate is distributed with the PDF f.

While justifying Algorithm 17.6, we did not use the fact that f is a normalized PDF. The acceptance-rejection method is often applied to complicated multivariate densities only known up to a multiplicative constant. Thus, the following generalization of Algorithm 17.6 is quite useful in dealing with such cases. Consider the sampling from a PDF proportional to some nonnegative integrable function f. Suppose that there exists another integrable function g so that it majorizes f, i.e., f(x) ≤ g(x) for all x. Let f and g have the same support D. The sampling algorithm is as follows.

A proposed draw X is accepted if the point (X, Y) being sampled uniformly in the body of g belongs to the body of f. Therefore, the probability of accepting X equals the ratio of the volumes of BD(f) and BD(g):

$\mathbb{P}\left(\text{Accept}\right)=\frac{|B_D\left(f\right)|}{|B_D\left(g\right)|}=\frac{{\displaystyle {\int }_{D}f\left(\text{x}\right)\text{dx}}}{{\displaystyle {\int }_{D}g\left(\text{x}\right)\text{dx}}}.$

To maximize this probability, we need to choose g as close to f as possible (see Figure 17.5). The average number of trials per one accepted draw (the computational cost of the acceptance-rejection method) is

$\text{Cost}=\mathbb{P}{\left(\text{Accept}\right)}^{-1}=\frac{{\displaystyle {\int }_{D}g\left(\text{x}\right)\text{dx}}}{{\displaystyle {\int }_{D}f\left(\text{x}\right)\text{dx}}}.$

Example 17.9.

Develop an acceptance-rejection method for the PDF

$f\left(x\right)=\frac{2\text{arcin}x}{\pi -2},0<x<1.$

Solution. By using the property that arcsin $x\le \frac{\pi }{2}x\text{for}0<x<1$, we obtain that

$f\left(x\right)\le g\left(x\right):=\frac{\pi }{\pi -2}x\text{for}0<x<1.$

The ratio of f and g is $\frac{\arcsin x}{\left(\pi /2\right)x}$. The proposal PDF p ∝ g is given by $p\left(x\right)=2x{\mathbb{I}}_{\left(0,1\right)}\left(x\right)$. The sampling algorithm is as follows.

1. (1) Generate two i.i.d. U1, U2 ~ Unif(0, 1).
2. (2) Sample X ~ p by using the inverse CDF method: $X=\sqrt{U_1}$.
3. (3) Accept X if $U_2\le \frac{\arcsin X}{\left(\pi /2\right)X}$. Otherwise, return to (1).

Example 17.10.

Develop an acceptance-rejection method for the standard normal distribution using the double-exponential sampling distribution as a proposal one. Find the computational cost.

Solution. The standard normal PDF n is proportional to $f\left(x\right)={\text{e}}^{-x^2/2}$. We have the following upper bound for f:

$\exp \left(-\frac{x^2}{2}\right)=\exp \left(-\frac{x^2-2|x|+1}{2}-|x|+\frac{1}{2}\right)=\exp \left(-\frac{{\left(|x|-1\right)}^2}{2}\right)\sqrt{\text{e}}{\text{e}}^{-|x|}\le \sqrt{\text{e}}{\text{e}}^{-|x|}.$

So the majorizing function is $g\left(x\right)=\sqrt{\text{e}}{\text{e}}^{-|x|}$. Therefore, the proposal probability distribution is the double-exponential distribution with the PDF $p\left(x\right)=\frac{1}{2}{\text{e}}^{-|x|}$, which can be expressed as a mixture:

$p\left(x\right)=\frac{1}{2}{\text{e}}^{-x}{\mathbb{I}}_{[0,\infty )}\left(x\right)+\frac{1}{2}{\text{e}}^{-|x|}{\mathbb{I}}_{[-\infty ,0)}\left(x\right).$

To sample from p, we first obtain a draw from Exp(1) and then assign a random sign to it:

$X=\{\begin{array}{l}Y\text{with}\text{probability}\frac{1}{2},\hfill \\ -Y\text{with}\text{probability}\frac{1}{2},\hfill \end{array}\text{where}Y~E\text{xp}\left(1\right).$

As a result, we obtain the following algorithm.

1. (1) Generate three i.i.d. U1, U2, U1 ~ Unif(0, 1).
2. (2) Sample X ~ p by using the composition method: $X=\mathrm{sgn} \left(U_1-\frac{1}{2}\right)\ln U_2$.
3. (3) Accept X if $U_3\le \frac{f\left(X\right)}{g\left(X\right)}=\exp \left(-\frac{{\left(|X|-1\right)}^2}{2}\right)$. Otherwise, return to (1).

The probability of acceptance PA is

$P_A=\frac{{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-x^2/2}\text{d}x}}{\sqrt{\text{e}}{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-|x|}\text{d}x}}=\frac{\sqrt{2\pi }}{2\sqrt{\text{e}}}=\sqrt{\frac{\pi }{2\text{e}}}\cong \mathrm{0.7602.}$

Therefore, the computational cost is $\text{E}\left[\#\text{trials}\text{per}\text{acceptance}\right]=\frac{1}{P_A}\cong 1.3155$.

17.3.5 Multivariate Sampling

17.3.5.1 Sampling by Conditioning

A d-variate joint PDF fX of the random vector X = [X1, X2, . . . , Xd]⊤ can be represented as a product of univariate conditional densities:

$f_{\text{X}}\left(\text{X}\right)=f_{X_1}\left(x_1\right)f_{X_2|X_1}\left(x_2|x_1\right)\cdots f_{X_d|X_1,...,X_{d-1}}\left(x_d|x_1,...,x_{d-1}\right),$

where x = [x1, x2, . . . , xd]⊤. The sampling procedure is as follows:

Step 1. Generate $X_1~f_{X_1}$.

Step 2. Generate X2 conditional on X1 from $f_{X_2|X_1}$.

⋮

Step d. Generate Xd conditional on X1, . . . , Xd−1 from $f_{X_d|X_{1,...,}{X}_{d-1}}$.

This method is simplified if the components Xj, 1 ≤ j ≤ d, are independent random variables. The joint PDF is then a product of marginal PDFs:

$f_{\text{X}}\left(\text{X}\right)={\displaystyle \prod _{j=1}^d{f}_{X_j}\left(x_j\right)}.$

Example 17.11.

Construct a sampling algorithm for a random vector X uniformly distributed in a hyperparallelepiped $D={\displaystyle {\prod }_{j=1}^d\left(a_j,b_j\right),a_j}<b_j,1\le j\le d$.

Solution. The joint PDF is a product of d marginal uniform densities:

$\begin{array}{l}{f}_{\text{X}}\left(X\right)=\frac{1}{|D|}{\mathbb{I}}_D\left(\text{X}\right)=\frac{1}{\left(b_1-a_1\right)\cdots \left(b_d-a_d\right)}{\mathbb{I}}_{\left(a_1,b_1\right)}\left(x_1\right)\cdots {\mathbb{I}}_{\left(a_d,b_d\right)}\left(x_d\right)\hfill \\ ={\displaystyle \prod _{j=1}^d\frac{1}{b_j-a_j}{\mathbb{I}}_{\left(a_j,b_j\right)}\left(x_j\right)}\equiv {\displaystyle \prod _{j=1}^d{f}_{X_j}\left(x_j\right)}.\hfill \end{array}$

Therefore, the vector X is formed of d i.i.d. uniformly distributed random variables:

$X_j=a_j+\left(b_j-a_j\right)U_j,U_j~Unif\left(0,1\right)1\le j\le d.$

17.3.5.2 The Box–Müller method

A pair of independent standard normal random variables Z1 and Z2 can be generated from two independent Unif(0, 1)-distributed random variables by using the following steps.

1. (1) Define the random variables R and Θ implicitly by

   $Z_1=R\cos \Theta \text{and}{\text{Z}}_{\text{2}}=R\sin \Theta .\text{ (17.17)}$

2. (2) One can show that R and Θ are independent random variables. Moreover, they can be simulated by the following formulae:

   $R=\sqrt{-2\ln U_1}\text{and}\Theta =2\pi U_2,\text{ (17.18)}$

   where U1 and U2 are independent Unif(0, 1)-distributed random variables.

3. (3) Therefore, Z1 and Z2 can be expressed in terms of U1 and U2 as follows:

   $Z_1=\sqrt{-2\ln U_1}\cos \left(2\pi U_2\right)\text{and}{Z}_2=\sqrt{-2\ln U_1}\sin \left(2\pi U_2\right).\text{ (17.19)}$

   To justify the Box–Müller method we apply the following theorem.

Theorem 17.11

(Bivariate Transformation Theorem, e.g., See Gut (2009)). Consider X and Y —jointly continuous random variables, and a one-to-one bivariate continuously differentiable transformation defined on the support of (X, Y) by u = g(x, y) and v = h(x, y). The joint PDF of U := g(X, Y) and V := h(X, Y) is $f_{U,V}\left(u,\upsilon \right)=\frac{1}{|J\left(x,y\right)|}{f}_{X,Y}\left(x,y\right)$, where (x, y) is a unique solution to $\{\begin{array}{l}g\left(x,y\right)=u\hfill \\ h\left(x,y\right)=v\hfill \end{array}$ and J(x, y) is the Jacobian determinant of the transformation defined by

$J\left(x,y\right)=\text{det}\left[\begin{array}{ll}\frac{\partial g}{\partial x}\left(x,y\right)\hfill & \frac{\partial g}{\partial y}\left(x,y\right)\hfill \\ \frac{\partial h}{\partial x}\left(x,y\right)\hfill & \frac{\partial h}{\partial y}\left(x,y\right)\hfill \end{array}\right].$

Since Z1 and Z2 are independent standard normal random variables, their joint PDF is

$f_{Z_1,Z_1}\left(z_1,z_2\right)=n\left(z_1\right)n\left(z_2\right)=\frac{1}{2\pi }{\text{e}}^{-\frac{1}{2}\left(z_1^2+z_2^2\right)}.$

The bivariate transformation theorem allows us to obtain a joint PDF of the pair (R, Θ). The Jacobian determinant of the transformation in (17.17) is equal to r. Thus,

$f_{Z_1,Z_2}\left(r\cos \theta ,r\sin \theta \right)=\frac{1}{r}{f}_{R,\Theta }\left(r,\theta \right)\Rightarrow f_{R,\Theta }\left(r,\theta \right)=r{\text{e}}^{-r^2/2}\frac{1}{2\pi },$

for r > 0 and 0 < θ < 2π. The joint PDF of R and Θ is a product of the marginal PDFs $f_R\left(r\right)=r{\text{e}}^{-r^2/2}{\mathbb{I}}_{[0,\infty )}\left(r\right)$ and $f_{\Theta }\left(\theta \right)=\frac{1}{2\pi }{\mathbb{I}}_{[0,2\pi )}\left(\theta \right)$. Therefore, R and Θ are independent random variable. Let us apply the inverse CDF method to generate R:

$F_R\left(r\right)=1-{\text{e}}^{-r^2/2}\Rightarrow R=F_R^{-1}\left(1-U_1\right)=\sqrt{-2\ln U_1,}$

where U1 ~ Unif(0, 1). Moreover, Θ ~ Unif(0, 2π), hence Θ = 2πU2 with U2 ~ Unif(0, 1). So, (17.18) is proved and (17.19) follows.

17.3.5.3 Simulation of Multivariate Normals

Consider a multivariate normal vector X = [X1, X2, . . . , Xd]⊤ ~ Normd(μ, Σ). If the covariance matrix Σ is a diagonal matrix diag $({\sigma }_1^2,.\; .\; .,{\sigma }_d^2)$, then Xj, 1 ≤ j ≤ d, are all independent normals which can be expressed in terms of independent standard normal variables as follows:

$X_j={\mu }_j+{\sigma }_j{Z}_j,Z_j~Norm\left(0,1\right),1\le j\le d.$

Independent standard normals can be generated by the Box–Müller method, by the acceptance-rejection method, or by the inverse CDF method. In the latter case we set $Z={\mathcal{N}}^{-1}\left(U\right)$ with U ~ Unif(0, 1), where the inverse normal CDF N−1(x) can be calculated numerically. One interesting application of standard normal random variables is the sampling of an isotropic vector in d dimensions.

1. (1) Sample i.i.d. Z1, . . . , Zd ~ Norm(0, 1)
2. (2) Define the vector X = [X1, X2, . . . , Xd]⊤ by Xj = Zj/Rd, 1 ≤ j ≤ d, where $R_d^2=Z_1^2+·\; ·\; ·+Z_d^2$. [Note that $R_d^2$ is a chi-square random variable with d degrees of freedom.]
3. (3) As a result, X is uniformly distributed on a unit d-dimensional sphere.

Suppose that the covariance matrix Σ has nonzero off-diagonal elements. Consider the following two general methods of sampling X:

1. (1) using the Cholesky factorization of the covariance matrix;
2. (2) using the conditional normal distribution.

Sampling by the Cholesky Factorization.

Let L be a lower-triangular matrix from the Cholesky factorization of Σ, i.e., Σ = L L⊤.

Let Z be a d-dimensional vector formed by i.i.d. standard normals Zj ~ Norm(0, 1), j = 1, 2, . . . , d. Then, we set

$X:=\mu +LZ~Nor{m}_d\left(\mu ,\Sigma \right).$

Conditional Normal.

Let us split the vector X into two parts:

$X=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right],\text{where}{X}_1\in {ℝ}^m\text{and}{X}_2\in {ℝ}^{d-m}$

for some 1 ≤ m < d. Split also the vector μ and matrix Σ to represent them in block form:

$\mu =\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right]\text{and}\Sigma =\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right],$

where μ1 ∊ ℝm and μ2 ∊ ℝd−m are vectors; Σ11 ∊ ℝm×m, Σ12 ∊ ℝm×(d−m), Σ21 ∊ ℝ(d−m)×m, and Σ22 ∊ ℝ(d−m)×(d−m) are matrices. Then, the conditional distribution of X1 given the value of X2 is normal:

$X_1|\left\{X_2={\text{x}}_{\text{2}}\right\}~Nor{m}_m\left({\mu }_1+{\Sigma }_{12}{\Sigma }_{22}^{-1}\left({\text{X}}_{\text{2}}-{\mu }_2\right),{\Sigma }_{11}-{\Sigma }_{12}{\Sigma }_{22}^{-1}{\Sigma }_{21}\right).\text{ (17.20)}$

Example 17.12.

Construct two methods of sampling from the trivariate normal distribution

$Nor{m}_3\left(\left[\begin{array}{c}3\\ 2\\ 4\end{array}\right],\left[\begin{array}{lll}9\hfill & 0\hfill & 0\hfill \\ 0\hfill & 4\hfill & 2\hfill \\ 0\hfill & 2\hfill & 3\hfill \end{array}\right]\right),$

using the Cholesky factorization and the conditional sampling approach, respectively.

Solution.

1. Find the Cholesky factorization of the covariance matrix. Let us solve the matrix equation Σ = L L⊤ to find L:

   $\left[\begin{array}{lll}9\hfill & 0\hfill & 0\hfill \\ 0\hfill & 4\hfill & 2\hfill \\ 0\hfill & 2\hfill & 3\hfill \end{array}\right]=\left[\begin{array}{lll}{\ell }_{11}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {\ell }_{21}\hfill & {\ell }_{22}\hfill & \hfill 0\hfill \\ {\ell }_{31}\hfill & {\ell }_{32}\hfill & {\ell }_{33}\hfill \end{array}\right]\left[\begin{array}{lll}{\ell }_{11}\hfill & {\ell }_{21}\hfill & {\ell }_{31}\hfill \\ \hfill 0\hfill & {\ell }_{22}\hfill & {\ell }_{32}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & {\ell }_{33}\hfill \end{array}\right]\Rightarrow L=\left[\begin{array}{lll}3\hfill & 0\hfill & \hfill 0\hfill \\ 0\hfill & 2\hfill & \hfill 0\hfill \\ 0\hfill & 1\hfill & \sqrt{2}\hfill \end{array}\right].$

   Thus, we obtain the following sampling formulae:

   $\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\end{array}\right]=\left[\begin{array}{c}3\\ 2\\ 4\end{array}\right]+\left[\begin{array}{lll}3\hfill & 0\hfill & \hfill 0\hfill \\ 0\hfill & 2\hfill & \hfill 0\hfill \\ 0\hfill & 1\hfill & \sqrt{2}\hfill \end{array}\right]\left[\begin{array}{c}{Z}_1\\ Z_2\\ Z_3\end{array}\right]\Rightarrow \{\begin{array}{l}{X}_1=3+3Z_1,\hfill \\ X_2=2+2Z_2,\hfill \\ X_3=4+Z_2+\sqrt{2}{Z}_3,\hfill \end{array}$

   where Z1, Z2, Z3 are i.i.d. standard normals.

2. First, sample X1 ~ Norm(3, 9): X1 = 3 + 3Z1. Second, sample X2 conditional on X1. Recall that two normal variables are independent iff they are uncorrelated. Hence X2 is independent of X1 since Cov(X1, X2) = 0. Therefore, $\left(X_2|X_1\right)\stackrel{d}{=}{X}_2~Norm\left(2,4\right)$ and X2 = 2 + 2Z2. Third, sample X3 conditional on X1 and X2. Again, X3 and X1 are uncorrelated, hence $\left(X_3|X_1,X_2\right)\stackrel{d}{=}\left(X_3|X_2\right)$. We have

   $\left[\begin{array}{c}{X}_3\\ X_2\end{array}\right]~Nor{m}_2\left(\left[\begin{array}{c}4\\ 2\end{array}\right],\left[\begin{array}{cc}3& 2\\ 2& 4\end{array}\right]\right)\Rightarrow X_3|X_2~Norm\left(3+\frac{X_2}{2},2\right).$

Thus, we have $X_3=3+\frac{X_2}{2}+\sqrt{2}{Z}_3=4+Z_2+\sqrt{2}{Z}_3$.

17.4 Simulation of Random Processes

A typical problem that requires simulation of sample paths of a stochastic process {X(t)}t≥0 is the estimation of a mathematical expectation of the form

$\text{E}\left[g\left(\left\{X\left(t\right):0\le t\le T\right\}\right)\right]\text{ (17.21)}$

with some function g of an X-path. There are several possible cases.

1. The function g depends on a discretely monitored skeleton of the process X:

   $g=g\left(X\left(t_1\right),X\left(t_2\right)...,X\left(t_m\right)\right),0\le t_1<t_2...<t_m\le T.$

   One special case is when g = g(X(T)). For example, the estimation of E[g(X(T))] is required to price a European-style option.

2. The function g depends on path-dependent quantities such as the running maximum/minimum of the process and the first passage time. It may be possible to sample such path-dependent quantities directly from their distributions rather than calculate them from a sample path.
3. The function g depends on a full sample path of process X on [0, T]. Since it may be not feasible to generate a complete sample path of a continuous-time process (unless we deal with a Poisson process or a similar process with piecewise paths that changes at a finite number of time points), such a full path can only be obtained by applying an interpolation algorithm to a path skeleton.

So our goal is to sample a path skeleton

$X\left(t_1\right),X\left(t_2\right),...,X\left(t_m\right)\text{for}0\le t_1<t_2<\cdots <t_m\le T.$

The skeleton can be generated from its exact multivariate distribution. In this case, the problem (17.21) can be reduced to the estimation of a multivariate integral of the form

${\displaystyle {\int }_{R^m}g\left(x_1,x_2,...,x_m\right)}{f}_{X\left(t_1\right),X\left(t_2\right),...,X\left(t_m\right)}\left(x_1,x_2,...,x_m\right)\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_m$

where f is a joint PDF of X(t1), X(t2), . . . , X(tm). Another approach is to sample an approximation path by applying some discretization scheme. Note that Brownian motion and other Gaussian processes as well as some jump processes can be sampled precisely form their path distributions. General diffusions can be simulated approximately by using, for example, the Euler approximation scheme.

17.4.1 Simulation of Brownian Processes

17.4.1.1 Sequential Sampling

The sequential sampling of Brownian motion (BM) and geometric Brownian motion is based on the property that Brownian increments on nonoverlapping intervals are independent. Consider a scaled Brownian motion with drift, $W_{x_0}^{\left(\mu ,\sigma \right)}\left(t\right):=x_0+\mu t+\sigma W\left(t\right)$. The standard BM is recovered from the process $W_{x_0}^{\left(\mu ,\sigma \right)}$ if we take x0 = 0, μ = 0, and σ = 1. Suppose that the process $W_{x_0}^{\left(\mu ,\sigma \right)}$ is to be sampled at a set of time points 0 = t0 < t1 < t2 < · · · < tm. For all j ≥ 1, we have

$W_{x_0}^{\left(\mu ,\sigma \right)}\left(t_j\right)=x_0+\mu t_j+\sigma W\left(t_j\right)=W_{x_0}^{\left(\mu ,\sigma \right)}\left(t_{j-1}\right)+\mu \left(t_j-t_{j-1}\right)+\sigma \left(W\left(t_j\right)-W\left(t_{j-1}\right)\right).$

Since the increment W(tj) − W(tj−1) ~ Norm(0, tj − tj−1) is independent of $W_{x_0}^{\left(\mu ,\sigma \right)}\left(t_{j-1}\right)$, we obtain the following simple algorithm.

The sample path of a geometric Brownian motion $S\left(t\right)=S_0{\text{e}}^{\mu t+\sigma W\left(t\right)}={\text{e}}^{\ln S_0+\mu t+\sigma W\left(t\right)}$ can be obtained by taking the exponential function of a sample path of the scaled BM with drift $W_{s_0}^{\left(\mu ,\sigma \right)}\left(t\right)$ that starts at x0 = ln S0, i.e., $S\left(t\right)=\exp \left(W_{\ln S_0}^{\left(\mu ,\sigma \right)}\left(t\right)\right)$.

Now we consider is a multidimensional BM $W\left(t\right)={\left[W_1\left(t\right),W_2\left(t\right),...,W_d\left(t\right)\right]}^{{\rm T}}$. Each component of W(t) is a standard Brownian motion. Suppose that the processes Wj, 1 ≤ j ≤ d, are correlated. For 1 ≤ i, j ≤ d, the correlation coefficient between Wi(t) and Wj(t) is

${\rho }_{ij}=\text{Corr}\left(W_i\left(t\right),W_j\left(t\right)\right)=\frac{\text{E}\left[W_i\left(t\right)W_j\left(t\right)\right]-\text{E}\left[W_i\left(t\right)\right]\text{E}\left[W_j\left(t\right)\right]}{\sqrt{\text{E}\left[W_i^2\left(t\right)\right]E\left[W_j^2\left(t\right)\right]}}=\frac{\text{E}\left[W_i\left(t\right)W_j\left(t\right)\right]}{t}.$

Let R = [ρij]1≤i,j≤d be the correlation matrix. R is a positive definite matrix with ones on the main diagonal. If we deal with independent Brownian motions then R = I. Apply the Cholesky factorization to find a lower triangular matrix L so that R = L L⊤. For example, for the two-dimensional case we have

$R=\left[\begin{array}{cc}1& {\rho }_{12}\\ {\rho }_{12}& 1\end{array}\right]=L{L}^{⊺}\Rightarrow L=\left[\begin{array}{cc}1& 0\\ {\rho }_{12}& \sqrt{1-{\rho }_{12}}\end{array}\right].$

Algorithm 17.9 allows us to obtain a realization of W at time points t0, t1, . . . , tm with 0 = t0 < t1 < · · · < tm.

17.4.1.2 Bridge Sampling

Previously, we derived the probability distribution of Brownian motion pinned at the endpoints of a time interval. Recall that Brownian motion conditional on W (0) = a and W (T) = b is called a Brownian bridge from a to b on [0, T]. There exist several applications of the Brownian bridge. First, the bridge distribution can be used to refine a sample skeleton. Second, it can be used as an alternative to the sequential simulation method for sampling a Brownian trajectory.

Suppose that a standard BM is sampled at m time moments 0 = t0 < t1 < · · · < tm = T and we wish to sample W (s) at some additional time moment s ∊ (0, T) conditional on these values. Let s ∊ (tj, tj+1) for some j ∊ {0, 1, . . . , m−1}. It follows from the Markov property of BM that

$\left(W\left(s\right)|\left\{W\left(t_i\right)=x_i,0\le i\le m\right\}\right)\stackrel{d}{=}\left(W\left(s\right)|\left\{W\left(t_j\right)=x_{j}W\left(t_{j+1}\right)=x_{j+1}\right\}\right).$

Hence W (s) can be sampled from the distribution of a Brownian bridge on [ti, tj+1]. By applying this procedure, a sample path can be refined without re-sampling its values at t1, t2, . . . , tm.

Consider the so-called dyadic partition of the time interval [0, T] with m = 2k points $t_j=\frac{j}{m}T$, where j = 0, 1, . . . , m and k ≥ 1. Let a realization of Brownian motion be sampled as follows:

$\begin{array}{ll}\text{Step}\text{1:}\hfill & \text{sample}W\left(t_m\right)\text{conditional}\text{on}W\left(t_0\right)=0,\hfill \\ \text{Step}\text{2:}\hfill & \text{sample}W\left(t_{m/2}\right)\text{conditional}\text{on}W\left(t_0\right),W\left(t_m\right),\hfill \\ \text{Step}3:\hfill & \text{sample}W\left(t_{m/4}\right)\text{conditional}\text{on}W\left(t_0\right),W\left(t_{m/2}\right),\hfill \\ \text{Step}4:\hfill & \text{sample}W\left(t_{3m/4}\right)\text{conditional}\text{on}W\left(t_{m/2}\right),W\left(t_m\right),\hfill \\ ⋮\hfill & \hfill \\ \text{Step}m:\hfill & \text{sample}W\left(t_{m-1}\right)\text{conditional}\text{on}W\left(t_{m-2}\right),W\left(t_m\right).\hfill \end{array}$

In other words, first we sample W(tm) and after that for each tj, 1 ≤ j ≤ m − 1, W(tj) is sampled conditionally on W (tℓ) and W(tk) previously generated, where the indices ℓ and k satisfy 0 ≤ ℓ < j < k ≤ m and $j=\frac{\ell +k}{2}$. As a result, a trajectory of BM is sampled at the time points in the following order of generation:

$\underset{}{\underbrace{t_m}},\underset{}{\underbrace{t_{m/2}}},\underset{\underset{}{\underbrace{t_2,t_6,t_{10},...,t_{m-2}}},}{\underbrace{t_{m/4},t_{3m/4}}},\underset{\underset{}{\underbrace{t_1,t_3,...,t_{m-1}}}}{\underbrace{t_{m/8},t_{3m/8},t_{5m/8},t_{7m/8},...}}\text{ (17.22)}$

Bridge sampling with m = 8 time points is illustrated in Figure 17.6.

The bridge sampling algorithm is useful in pricing path-dependent financial instruments. Being applied with (randomized) low-discrepancy numbers, i.e., when the (randomized) quasi-Monte Carlo method is used, it allows us to reduce the variance of a path-dependent estimator. Another advantage is that the bridge sampling algorithm can be easily parallelized.

17.4.2 Simulation of Gaussian Processes

In accordance with the definition of a Gaussian process {X(t)}t≥0, for any sequence of time points 0 < t1 < t2 < · · · < tm, the vector X = [X(t1), X(t2), . . . , X(tm)]⊤ has a multivariate normal distribution Normm(μm, Σm) with

${\mu }_m=\left[\begin{array}{c}{m}_X\left(t_1\right)\\ m_X\left(t_2\right)\\ ⋮\\ m_X\left(t_m\right)\end{array}\right]\text{and}{\Sigma }_m=\left[\begin{array}{llll}{c}_X\left(t_1,t_1\right)\hfill & c_X\left(t_1,t_2\right)\hfill & \cdots \hfill & c_X\left(t_1,t_m\right)\hfill \\ c_X\left(t_2,t_1\right)\hfill & c_X\left(t_2,t_2\right)\hfill & \cdots \hfill & c_X\left(t_2,t_m\right)\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \ddots \hfill & \hfill ⋮\hfill \\ c_X\left(t_m,t_1\right)\hfill & c_X\left(t_m,t_2\right)\hfill & \cdots \hfill & c_X\left(t_m,t_m\right)\hfill \end{array}\right],$

where mX(t) = E[X(t)] and cX(t, s) = Cov(X(t), X(s)) are, respectively, the mean and covariance functions. Therefore, the realization of X at time points t1, t2, . . . , tm can be constructed by sampling from the m-variate normal distribution Normm(μm, Σm) as follows.

1. (1) Apply the Cholesky factorization to find a lower-triangular matrix Lm so that

   ${\Sigma }_m=L_m{L}_m^{⊺}.$

2. (2) Sample m i.i.d. standard normals Z1, Z2, . . . , Zm.
3. (3) Set X = μm + Lm Z where Z = [Z1, Z2, . . . , Zm]⊤.

This generic algorithm can be applied to any Gaussian process including those listed below.

Brownian Motion $W_{x_0}^{\left(\mu ,\sigma \right)}\left(t\right):=x_0+\mu t+\sigma W\left(t\right)$ with m(t) = x0 + μt and c(t, s) = σ2 (tΛs) for s, t ≥ 0.

Itô Processes $X\left(t\right):={\displaystyle {\int }_0^t\mu \left(s\right)\text{d}s+{\displaystyle {\int }_0^t\sigma \left(t\right)\text{d}W\left(\mathcal{S}\right)}}$ with $m\left(t\right)={\displaystyle {\int }_0^t\mu \left(u\right)\text{d}u}$ and $c\left(t,\mathcal{S}\right)={\displaystyle {\int }_0^{t\wedge \mathcal{S}}{\sigma }^2\left(u\right)\text{d}u\text{for}\mathcal{S},t\ge 0}$.

Fractional Brownian Motion W(H)(t) with c(t, s) = (t2H + s2H − |t − s|2H)/2 and m(t) = 0. Here, H ∊ [0, 1] is the so-called Hurst parameter. Note that W(1/2) is a standard BM with c(t, s) = (t + s − |t − s|)/2 = t Λ s for s, t ≥ 0.

Standard Brownian Bridge from a to b on [0, T], denoted ${\left\{B_{\left[0,T\right]}^{a,b}\left(t\right)\right\}}_{0\le t\le T}$, with $m\left(t\right)=\frac{a\left(T-t\right)+bt}{T}$ and $c\left(t,s\right)=t\wedge s\left(T-t\vee s\right)\text{for}0\le s,t\le T$.

17.4.3 Diffusion Processes: Exact Simulation Methods

A diffusion process {X(t)}t≥0 is a solution to an initial value problem for a stochastic differential equation (SDE):

$\{\begin{array}{l}\text{d}X\left(t\right)=\mu \left(t,X\left(t\right)\right)\text{d}t+\sigma \left(t,X\left(t\right)\right)\text{d}W\left(t\right),t\ge 0,\hfill \\ X\left(0\right)=X_0.\hfill \end{array}\text{ (17.23)}$

The process can also be written in the integral form

$X\left(t\right)=X_0+{\displaystyle {\int }_0^t\mu \left(s,X\left(s\right)\right)\text{d}s+}{\displaystyle {\int }_0^t\sigma \left(s,X\left(s\right)\right)\text{d}W\left(s\right).}\text{ (17.24)}$

The functions μ and σ are called the drift coefficient and the diffusion coefficient, respectively.

The problem (17.23) (or (17.24)) can be solved analytically or numerically. One approach is to represent the process X as an explicit function of the underlying Brownian motion W :

$X\left(t\right)=f\left(t,\left\{W\left(s\right):0\le s\le t\right\}\right).\text{ (17.25)}$

Such an explicit representation is called a strong solution to (17.23). As a result of (17.25), a sample path of the X-process is obtained by transforming a Brownian trajectory. Another approach is to find the transition PDF for X by solving the Kolmogorov equation. The strong solution and/or the transition density can be used to generate a realization of the diffusion X from its exact finite-dimensional distribution. As usual, our goal is to generate a path skeleton for an arbitrary sequence of time points. Alternatively, being unable to analytically solve (17.23), one can apply a numerical scheme to find an approximate realization of the diffusion process. The Euler scheme, which is the simplest and most popular simulation method, is considered in Section 17.4.4.

17.4.3.1 The Stochastic Calculus Approach

Most of the SDEs for which we can find a strong solution are SDEs with linear (w.r.t. the space variable) drift and diffusion coefficients. Let us consider some examples of such SDEs and derive transition probability distributions of their solutions.

1. Geometric Brownian motion is a solution to

   $\text{d}X\left(t\right)=\mu X\left(t\right)\text{d}t+\sigma X\left(t\right)\text{d}W\left(t\right),X\left(0\right)=X_0.$

   Applying the Itô formula gives $X\left(t\right)=X_0{\text{e}}^{\left(\mu -{\sigma }^2/2\right)t+\sigma W\left(t\right)}$. To model the transition X(s) → X(t) for 0 ≤ s < t, we use the following representation:

   $X\left(t\right)=X\left(s\right){\text{e}}^{\left(\mu -{\sigma }^2/2\right)\left(t-s\right)+\sigma \left(W\left(t\right)-W\left(s\right)\right)}\stackrel{d}{=}X\left(s\right){\text{e}}^{\left(\mu -{\sigma }^2/2\right)\left(t-s\right)+\sigma \sqrt{t-s}Z},$

   where Z ~ Norm(0, 1) is independent of X(s).

2. The solution to

   $\text{d}X\left(t\right)=\mu \left(t\right)\text{d}t+\sigma \left(t\right)\text{d}W\left(t\right),X\left(0\right)=X_0,$

   is a Gaussian process given by $X\left(t\right)=X_0+{\displaystyle {\int }_0^t\mu \left(u\right)\text{d}u}+{\displaystyle {\int }_0^t\sigma \left(u\right)\text{d}W\left(u\right)}$. By representing X(t) in terms of X(s) for 0 ≤ s < t and using the fact that the Itô integral ${\displaystyle {\int }_{\mathcal{S}}^t\sigma \left(u\right)}\text{d}W\left(u\right)$ is normally distributed with mean 0 and variance ${\displaystyle {\int }_{\mathcal{S}}^t{\sigma }^2\left(u\right)\text{d}u}$, we obtain that the increments of X are normal:

   $X\left(t\right)-X\left(s\right)~Norm\left({\displaystyle {\int }_s^t\mu \left(u\right)\text{d}u,}{\displaystyle {\int }_s^t{\sigma }^2\left(u\right)\text{d}u}\right).$

   Since the increments of X are independent, we are able to model the transition X(s) → X(t) for all 0 ≤ s < t.

3. The Ornstein–Uhlenbeck process is a solution to the SDE with a constant diffusion coefficient and linear drift:

   $\text{d}X\left(t\right)=\alpha \left(b-X\left(t\right)\right)\text{d}t+\sigma \text{d}W\left(t\right),X\left(0\right)=X_0.$

   The strong solution is

   $X\left(t\right)={\text{e}}^{-\alpha t}{X}_0+\alpha b{\displaystyle {\int }_0^t{\text{e}}^{-\alpha \left(t-u\right)}\text{d}u+\sigma }{\displaystyle {\int }_0^t{\text{e}}^{-\alpha \left(t-u\right)}\text{d}W\left(u\right)}.$

   The conditional distribution of X(t) given X(s) for 0 ≤ s < t is normal:

   $\left(X\left(t\right)|X\left(s\right)\right)~Norm\left({\text{e}}^{-\alpha \left(t-s\right)}X\left(s\right)+\alpha b{\displaystyle {\int }_s^t{\text{e}}^{-\alpha \left(t-u\right)}\text{d}u,{\sigma }^2}{\displaystyle {\int }_s^t{\text{e}}^{-2\alpha \left(t-u\right)}\text{d}u}\right).$

17.4.3.2 The PDF Approach

Consider a Markov stochastic process {X(t)}t≥0 starting at X0 so that its transition PDF p given by p(s, t; y, x) dx = ℙ(X(t) ∊ dx | X(s) = y), 0 ≤ s < t, x, y ∊ ℝ, is known in closed form. Moreover, suppose that it is feasible to sample from the PDF p(s, t; y, ·). The general simulation algorithm is as follows.

For example, due to the time- and space-homogeneity property, the Brownian motion transition PDF p(s, t; y, x) reduces to a two-variable function p0 given by

$p_0\left(t;x\right)=\frac{1}{\sqrt{2\pi t}}\exp \left(-\frac{x^2}{2t}\right)$

as follows: p(s, t; y, x) = p0(t − s; x − y). The latter solves the PDE ${\partial }_t{p}_0=\frac{1}{2}{\partial }_{xx}{p}_0$. To generate (W(tj)|{W(tj−1) = y}) ~ p(tj−1, tj; y, x) = p0(tj − tj−1; x − y), we first sample Z ~ Norm(0, 1) and then set $W\left(t_j\right)\leftarrow y+\sqrt{t_j-t_{j-1}}Z$.

To sample from a transition PDF, we can use a whole arsenal of sampling techniques such as the inversion method, composition approach, and acceptance-rejection technique. The main example considered in the current section is the exact simulation of a family of Bessel diffusions. We start with a squared Bessel (SQB) process and show that its transition PDF reduces to a randomized gamma distribution. Moreover, as demonstrated in Chapter 16, other processes such as the Cox–Ingersolll–Ross (CIR) process and the constant elasticity of variance (CEV) diffusion model can be obtained from the SQB process by means of a scale and time transformation and change of variable.

Let us consider a λ0-dimensional squared Bessel (SQB) process {X(t) ∊ ℝ+}t≥0 obeying the stochastic differential equation (SDE)

$dX\left(t\right)={\lambda }_0\text{d}t+v\sqrt{X\left(t\right)}\text{d}W\left(t\right),\text{ (17.26)}$

with constant parameters λ0 and ν > 0. For simplicity of presentation, we assume here that ν = 2. The process X is a time-homogeneous Markov process and its transition PDF is given by (16.15).

The transition PDF (16.15) looks very similar to the PDF (17.8) of the noncentral chi-square distribution. In fact, for the case when μ ≥ 0 or μ ∊ (−1, 0) and x = 0 is a reflecting boundary (so there is no absorption at the origin), the transition PDF p(t; y, x) of the SQB process reduces to the PDF f(x; κ, λ) in (17.8) as follows:

$p\left(t;y,x\right)=\frac{1}{t}f\left(\frac{x}{t};\kappa =2\mu +2,\lambda =\frac{y}{t}\right).$

Following the composition method, the noncentral chi-square PDF (17.8) can be represented as a randomized gamma distribution (17.10). The value of X(t) conditional on X(s) = y has the randomized gamma distribution

$Gamma\left(Y+\mu +1,2\left(t-s\right)\right),\text{where}Y~Pois\left(\frac{y}{2\left(t-s\right)}\right),\text{ (17.27)}$

for 0 ≤ s < t, y > 0. Therefore, we have the following sampling algorithm for the SQB process without absorption.

Now consider the case when μ < 0 and x = 0 is an exit or a killing boundary. The stochastic process X admits absorption at the origin. It can be shown that the transition PDF p given by (16.15) with $\tilde{\mu }=\left|\mu \right|,\mu <0$ does not integrate to one. Let us define the probability Psv of surviving before time t and the probability Pab of absorption before time t for the SQB process starting at x0 > 0:

$P_{s\upsilon }\left(x_0;t\right)={\displaystyle {\int }_0^{\infty }p\left(t;x_0,x\right)\text{d}x>}0\text{and}{P}_{ab}\left(x_0;t\right)=1-P_{s\upsilon }\left(x_0;t\right)>0.$

The probabilities of surviving and absorption of the SQB process before time t are

$P_{s\upsilon }\left(x;t\right)=\mathbb{P}\left\{T_0>t\right\}=\frac{\gamma \left(|\mu |,\frac{x}{2t}\right)}{\Gamma \left(|\mu |\right)}\text{and}{P}_{ab}\left(x;t\right)=\mathbb{P}\left\{T_0\le t\right\}=\frac{\Gamma \left(|\mu |,\frac{x}{2t}\right)}{\Gamma \left(|\mu |\right)},$

respectively, where τ0 is the first hitting time (FHT) at zero. Here,

$\gamma \left(a,x\right)={\displaystyle {\int }_0^x{t}^{a-1}}{\text{e}}^{-t}\text{d}t\text{and}\Gamma \left(a,x\right):=\Gamma \left(a\right)-\gamma \left(a,x\right)$

are, respectively, the lower and upper incomplete gamma functions.

Observe that the actual transition probability distribution is then a mixture of continuous and discrete probability distributions with the following generalized density:

$f\left(X\left(s\right)\to X\left(t\right)\right)=P_{s\upsilon }\left(X\left(s\right);t-s\right)\cdot \left(\frac{p\left(t-s;X\left(s\right),X\left(t\right)\right)}{P_{s\upsilon }\left(X\left(s\right);t-s\right)}\right)+P_{ab}\left(X\left(s\right);t-s\right)\cdot \delta \left(X\left(t\right)\right),$

for 0 ≤ s < t. Here, δ denotes the Dirac delta function that can be viewed as a generalized density of the discrete distribution with an only mass point at zero. With the probability Pab, the process is absorbed at zero. With the additional probability Psv, the process survives. The normalized transition PDF of the SQB process conditioned on the survival of the process before time t is

$\frac{p\left(t;x_0,x\right)}{P_{s\upsilon }\left(x_0;t\right)}=\frac{\Gamma \left(|\mu |\right)}{\gamma \left(|\mu |,\frac{x_0}{2t}\right)}{\left(\frac{x}{x_0}\right)}^{\frac{\mu }{2}}\frac{{\text{e}}^{-\left(x+x_0\right)/\left(2t\right)}}{2t}{I}_{|\mu |}\left(\frac{\sqrt{x{x}_0}}{t}\right).\text{ (17.28)}$

As is seen, the function on the right-hand side of (17.28) reduces to the form of (17.15) with ν = |μ|, λ = x0/(2t), and θ = 2t. Thus, the above normalized transition PDF follows the randomized gamma distribution of the third kind, Gamma(Y + 1, 2t), where Y ~ IΓ(|μ|, x0/(2t)). As a result, we obtain the following sampling algorithm that returns a sample path [X(t1), X(t2), . . . , X(tm)] and an approximation ${\tilde{T}}_0\in \left\{t_1,...,t_m,\infty \right\}$ of the FHT τ0.

17.4.4 Diffusion Processes: Approximation Schemes

Approximation schemes can be used whenever an exact method is not available or not feasible. In particular, approximation methods are efficient for the numerical solution of multidimensional stochastic differential equations. Computational schemes for SDEs are based on the same ideas as the numerical methods for solving deterministic differential equations. A discrete-time approximation solution is calculated on a time grid with small time steps. To illustrate the main idea, let us first consider a Cauchy problem for a system of ordinary differential equations written in a vector form:

$\frac{\text{dx}\left(t\right)}{\text{d}t}=\text{a}\left(t,\text{x}\left(t\right)\right),t\in \left[0.T\right];\text{x}\left(0\right)={\text{x}}_{\text{0}}.$

The solution admits the following integral representation on a small time interval [t, t + h]:

$\text{x}\left(t+h\right)=\text{x}\left(t\right)+{\displaystyle {\int }_t^{t+h}\text{a}}\left(s,\text{x}\left(s\right)\right)\text{d}s.\text{ (17.29)}$

By applying a rectangle quadrature rule to the integral in (17.29), we derive the so-called Euler approximation scheme:

$\text{x}\left(t+h\right)\approx \text{x}\left(t\right)+\text{a}\left(t,\text{x}\left(t\right)\right)h.$

A discrete-time numerical solution ${\left\{X_k^h\right\}}_{k=0,1,2,...}$ on the time grid {tk := kh}k=0,1,2,... is given by

${\text{X}}_0^h={\text{X}}_{\text{0}},{\text{X}}_{k+1}^h={\text{X}}_k^h+\text{a}\left(t_k,{\text{X}}_k^h\right)h,k=0,1,2,...$

The numerical solution approximates the genuine solution: $X_k^h\approx X\left(t_k\right),k\ge 1$; it converges to the exact one as the time step h goes to zero. Applying a similar approach to stochastic differential equations, we can construct a discrete-time approximate realization (a skeleton) of a sample path. Since we deal with stochastic processes, there are different types of convergence of the numerical solution to the exact one.

17.4.4.1 Types of Convergence

Consider a continuous-time d-dimensional stochastic process {X(t)}t∊[0,T] and its discrete-time approximation ${\left\{X_k^h\right\}}_{0\le k\le m}$ defined on a time grid 0 = t0 < t1 < · · · < tm = T with the maximum step size h = max{tk − tk−1 : k = 1, 2, . . . , m}. The simplest time grid is an equally spaced grid with tk = kh with $h=\frac{T}{m}$ and k = 0, 1, . . . , m. Let us analyze the convergence of the approximate solution Xh to the genuine solution X, as h → ∞, in the Euclidian vector norm ${||x||}_2=\sqrt{x_1^2+...+x_d^2}$, where x = [x1, x2, . . . , xd]⊤ ∊ ℝd. We say that the approximation Xh has:

• the strong order of convergence α > 0 if there exists C > 0 so that E

  $\text{E}\left[||X_k^h-X\left(t_k\right)||\right]\le C{h}^{\alpha },\forall k=0,1,2,...,m;$

• the mean-square order of convergence β > 0 if there exists C > 0 so that

  ${\left(\text{E}\left[{||X_k^h-X\left(t_k\right)||}^2\right]\right)}^{1/2}\le C{h}^{\beta },\forall k=0,1,2,...,m;$

• the weak order of convergence γ > 0 if for any real-valued function f selected from some large class of functions (usually, f is a sufficiently smooth function) $\exists C_f>0$ so that

  $|\text{E}\left[f\left(X_m^h\right)-f\left(X\left(T\right)\right)\right]|\le C_f{h}^{\gamma }.$

17.4.4.2 The Euler Scheme

Consider a multivariate SDE

$\text{d}X\left(t\right)=\text{a}\left(t,X\left(t\right)\right)\text{d}t+\text{b}\left(t,X\left(t\right)\right)\text{dW}\left(t\right),\text{ (17.30)}$

where W(t) is an d-dimensional standard Brownian motion with independent components, and a: ℝ × ℝd → ℝd and b: ℝ × ℝd → ℝd×d are, respectively, the drift and diffusion coefficient functions. The SDE (17.30) can be written in the integral form on [t, t + h], t ≥ 0, h > 0:

$X\left(t+h\right)=X\left(t\right)+{\displaystyle {\int }_t^{t+h}\text{a}}\left(s,X\left(s\right)\right)\text{d}s+{\displaystyle {\int }_t^{t+h}\text{b}\left(s,X\left(s\right)\right)\text{d}W}\left(s\right).\text{ (17.31)}$

Application of the rectangular approximation formula to each integral in (17.31) gives

$X\left(t+h\right)\approx X\left(t\right)+\text{a}\left(t,X\left(t\right)\right)h+\text{b}\left(t,X\left(t\right)\right)\left(W\left(t+h\right)-W\left(t\right)\right).\text{ (17.32)}$

So, the Euler discrete-time approximation on a time grid 0 = t0 < t1 < · · · < tm = T is given by

$X_0^h=X_0,X_{k+1}^h=X_k^h+\text{a}\left(t_k,X_k^h\right)h+\text{b}\left(t_k,X_k^h\right)\sqrt{t_k-t_{k-1}}{Z}_k,0\le k\le m-1,\text{ (17.33)}$

where {Zk}k≥0 is a sequence of i.i.d. multivariate Normd(0, I)-distributed vectors.

Theorem 17.12

The Euler scheme has strong order $\frac{1}{2}$ and weak order 1. Moreover, the weak error admits the following expansion:

$\text{E}\left[f\left(X_m^h\right)-f\left(X\left(T\right)\right)\right]=C_{f}h+\mathcal{O}\left(h^2\right).$

For a proof, see Kloeden and Platen (2011).

17.4.4.3 Extrapolation

The process of extrapolation uses two approximations computed from the same formula but with different step sizes to obtain higher-order approximation. The Richardson extrapolation method allows us to achieve second-order accuracy for a first-order scheme. For the Euler method with a constant time step h, we have

$\text{E}\left[f\left(X_m^h\right)\right]=\text{E}\left[f\left(X\left(T\right)\right)\right]+C_{f}h+\mathcal{O}\left(h^2\right).$

for some constant Cf that depends on f. Suppose that the number of steps m is even. Apply the Euler method with steps $h=\frac{T}{m}\text{and}2h=\frac{T}{m/2}$ to obtain approximations $X_m^h$ and $X_{m/2}^{2h}$ of the time-T realization X(T), respectively. By combining

$\text{E}\left[f\left(X_m^h\right)\right]=\text{E}\left[f\left(X\left(T\right)\right)\right]+C_{f}h+\mathcal{O}\left(h^2\right)\text{and}\text{E}\left[f\left(X_{m/2}^{2h}\right)\right]=\text{E}\left[f\left(X\left(T\right)\right)\right]+C_{f}2h+\mathcal{O}\left(h^2\right),$

we can eliminate the leading term of the error expansion:

$\text{E}\left[2f\left(X_m^h\right)-f\left(X_{m/2}^{2h}\right)\right]=\text{E}\left[f\left(X\left(T\right)\right)\right]+\mathcal{O}\left(h^2\right).$

The error of the combined estimate is of (weak) order 2. The variance of the combined estimator is

$\text{Var}\left(2f\left(X_m^h\right)-f(X_{m/2}^{2h})\right)=4\text{Var}\left(f(X_m^h)\right)+\text{Var}\left(f(X_{m/2}^{2h})\right)-4\text{Cov}\left(f\left(X_m^h\right),f(X_{m/2}^{2h})\right).$

By making $f\left(X_m^h\right)$ and $f(X_{m/2}^{2h})$ positively correlated, we can reduce the variance. This can be achieved by using consistent Brownian increments in simulating paths of Xh and X2h. Suppose that $\left(X_m^h\right)$ is constructed from m i.i.d. normal vectors

$\sqrt{h}{Z}_0,...,\sqrt{h}{Z}_{m-1}~Nor{m}_d\left(0,\sqrt{h}I\right).$

The same normal increments can be used for sampling X2h. For example, to construct $X_{m/2}^{2h}$ we use the normally distributed vectors

$\sqrt{h}\left(Z_0+Z_1\right),...,\sqrt{h}\left(Z_{m-2}+Z_{m-1}\right)~Nor{m}_d\left(0,\sqrt{2h}I\right).$

Here, the property that a sum of two normal random variables is again normally distributed is applied.

17.4.4.4 Error Analysis

Suppose that we wish to approximate some quantity Q by using a biased estimator Yh where h is a discretization parameter approaching zero. For example,

$Q:=\text{E}\left[f\left(X\left(T\right)\right)\right]\approx \text{E}\left[f\left(X_m^h\right)\right],$

where Xh is the Euler approximation of a diffusion X and $h=\frac{T}{m}$ is a time step size. Introduce the approximation error ε(h) = Q − E[Yh]. Suppose that

$\epsilon \left(h\right)\approx C_1{h}^{\beta }\text{ (17.34)}$

holds for some positive constants C1 and β. Taking the logarithm of both parts of (17.34) gives

$\log \epsilon \left(h\right)\approx \log C_1+\beta \log h.$

As is seen from the above equation, we can use the linear regression method to calculate C1 and β. First, we compute the sample estimate of Q for a decreasing sequence of values of h:

$Q\approx {\overline{y}}_n^{h_k}=\frac{1}{n}{\displaystyle \sum _{j=1}^n{y}_j^{h_k},k=1,2,...,}$

where $y_j^{h_k}$ are i.i.d. samples of $Y^{h_k}$ and h1 > k2 > h3 > . . .. For example, we set hk = 2−k. Assuming that the number of draws, n, is sufficiently large so that the statistical error is negligible in comparison with the approximation error, we obtain

$\overline{\epsilon }\left(h_k\right):=Q-{\overline{y}}_n^{h_k}\approx Q-\text{E}\left[Y^{h_k}\right]\approx C_1{h}^{\beta }\Rightarrow \log \overline{\epsilon }\left(h_k\right)\approx \log C_1+\beta \log h_k.$

Now, we plot $\log \overline{\epsilon }\left(h_k\right)$ versus log hk for all hk and then perform a linear regression. As a result, the slope of the regression line gives us the order of approximation β.

In fact, the error $\overline{\epsilon }\left(h_k\right)=Q-{\overline{y}}_h^h$ includes two components: the approximation bias and the statistical (Monte Carlo) error. To optimize the method, it is helpful to separate these two errors:

$Q-{\overline{y}}_n^h=\underset{\approx C_1{h}^{\beta }}{\underbrace{Q-\text{E}\left[Y^h\right]}}+\text{E}\left[Y^h\right]-{\overline{y}}_n^h.$

Define the mean-square (statistical) error as

$\text{MSE}\left(h,n\right)=\text{E}\left[{\left(Q-{\overline{Y}}_n^h\right)}^2\right],$

where ${\overline{Y}}_n^h=\frac{1}{n}{\displaystyle {\sum }_{j=1}^n{Y}_j^h}$ is the sample estimator. The value of MSE(h, n) is

$\begin{array}{cc}\hfill \text{E}\left[{\left(Q-{\overline{Y}}_n^h\right)}^2\right]=& \text{E}\left[{\left(\left(Q-\text{E}\left[Y^h\right]\right)+\left(\text{E}\left[Y^h\right]-{\overline{Y}}_n^h\right)\right)}^2\right]\hfill \\ \hfill =& \text{E}\left[{\left(Q-\text{E}\left[Y^h\right]\right)}^2\right]+2\text{E}\left[\left(Q-\text{E}\left[Y^h\right]\right)\left(\text{E}\left[Y^h\right]-{\overline{Y}}_n^h\right)\right]+\text{E}\left[{\left(\text{E}\left[Y^h\right]-{\overline{Y}}_n^h\right)}^2\right]\hfill \\ \hfill \approx & {\left(C_1{h}^{\beta }\right)}^2+2\left(Q-\text{E}\left[Y^h\right]\right)\text{E}\left[\text{E}[Y^h]-{\overline{Y}}_n^h\right]+\text{Var}\left({\overline{Y}}_n^h\right)\hfill \\ \hfill =& {\left(C_1{h}^{\beta }\right)}^2+\text{Var}\left({\overline{Y}}_n^h\right)=C_1^2{h}^{2\beta }+\frac{1}{n}\text{Var}\left(Y^h\right)\approx C_1^2{h}^{2\beta }+\text{Var}\left(Y^{0+}\right)\hfill \\ \hfill =& C_1^2{h}^{2\beta }+\frac{C_2}{n}.\hfill \end{array}$

Here the constant C2 = Var(Y0+) is defined as the limiting value of the variance Var(Yh), as $h\searrow 0$. To find the optimal values of the sample volume n and the discretization parameter h, we will minimize the computational cost for a given level of error: MSE = ε2. Clearly, the computational cost (i.e., the number of operations) is directly proportional to m and n. The number of steps $m=\frac{T}{h}$ is inversely proportional to h. Hence, the computational cost (or the runtime) is given by C3n/h, where C3 is another positive constant. Let us minimize the computational cost under the constraint that the mean-square error is fixed and equal to ε. The optimization problem takes the following form:

$\{\begin{array}{l}\frac{C_{3}n}{h}\to \underset{n,h}{\min },\hfill \\ \text{MSE}\left(h,n\right)=C_1^2{h}^{2\beta }+\frac{C_2}{n}={\in }^2.\hfill \end{array}\text{ (17.35)}$

The problem (17.35) is easy to solve. As a result, we obtain that the computational cost is proportional to ε−2−1/β. Alternatively, we can solve the problem of minimizing MSE subject to a fixed computational cost s (see Exercise 17.23).

17.4.5 Simulation of Processes with Jumps

17.4.5.1 Poisson Processes

Recall that a Poisson process is a continuous-time stochastic process {Nλ(t)}t≥0 with independent, stationary, Poisson distributed increments that starts at zero. In other words,

1. (a) Nλ(0) = 0;
2. (b) for all m and 0 ≤ t0 ≤ t1 ≤ · · · ≤ tm, the increments Nλ(tk) − Nλ(tk−1), 1 ≤ k ≤ m, are independent random variables;
3. (c) Nλ(t) − Nλ(s) ~ Pois(λ(t − s)) for 0 ≤ s < t.

The parameter λ is called the intensity of the process.

Every realization of a Poisson process is a step function that starts at the origin. It stays at each level k ≥ 0 for a random time period and then jumps to the next level k + 1. The occurrence time Tk is the time when the process jumps from k − 1 to k for k ≥ 1. Set τ1 = T1 and τk = Tk − Tk−1 for k ≥ 2. These variables {τk}k≥1 are called durations. We can express {Tk}k≥1 in terms of {τk}k≥1 as

$T_m={\displaystyle \sum _{k=1}^m{T}_k.}$

Thus, a Poisson process can be defined via durations as follows:

$N_{\lambda }\left(t\right)=\sup \left\{m:{\displaystyle \sum _{k=1}^m{T}_k\le t}\right\}.\text{ (17.36)}$

A realization of a Poisson process can be generated from the sample values of the durations. The probability distribution of the random variables {τk}k≥1 is given by the following result.

Proposition 17.13.

Consider a Poisson process with occurrence times Tk, k ≥ 1, and durations τ1 = T1, τk = Tk − Tk−1, k ≥ 2. Then,

1. (a) τk, k ≥ 1, are jointly independent, Exp(λ)-distributed random variables;
2. (b) Tk ~ Gamma(k, λ), k ≥ 1.

For a complete proof of Propositions 17.13 and 17.14, see Gut (2009).

By using Proposition 17.13 and Equation 17.36, we come up with the following sampling algorithm of the Poisson process on [0, T].

Another algorithm for sampling a Poisson process is based on conditioning on the number of occurrences in a time interval. As it turns out, the joint distribution of the occurrence times conditional on the number of occurrences is the same as that of the order statistics of a sample from a uniform distribution.

Proposition 17.14.

The joint density of T1, T2, . . . , Tm conditional on Nλ(T) = m is

$f_{T_1,...,T_m|N_{\lambda }\left(T\right)=m}\left(t_1,...,t_m\right)=\{\begin{array}{l}\frac{m!}{T^m}for0<t_2<\cdots <t_n<T,\hfill \\ 0otherwise.\hfill \end{array}$

In other words,

$\left(T_1,T_2,...,T_m|N_{\lambda }\left(T\right)=m\right)\stackrel{d}{=}\left(U_{\left(1\right)},U_{\left(2\right)},...,U_{\left(m\right)}\right),$

where U(1), U(2), . . . , U(m) are the order statistics defined by sorting m i.i.d. Unif(0, T)-distributed random variables in increasing order.

Algorithms 17.14 and 17.15 allow us to sample a Poisson path on [0, T]. If we want to continue sampling on (T, T + T ′], T ′ > 0, then the sample path on [0, T] can be reused thanks to the following property.

Proposition 17.15.

If {N(t)}t≥0 is a Poisson process, then so are

1. {N(t + s) − N(s)}t≥0 for every fixed s > 0;
2. {N(t+Tk)−N(Tk)}t≥0 for every fixed k ≥ 1, where Tk is the time of the kth occurrence of the original Poisson process N.

This result, along with the property of independence of Poisson increments, allows us to simulate a Poisson process individually on disjoint intervals. Suppose we have generated the process Nλ on [0, T]. To continue the sample path on (T, T + T ′], we first generate another Poisson process ${\left\{{\tilde{N}}_{\lambda }\left(t\right)\right\}}_{t\in \left[0,T^{\prime }\right]}$ independently of {Nλ(t)}t∊[0,T]. Second, we set $N_{\lambda }\left(t\right)=N_{\lambda }\left(T\right)+{\tilde{N}}_{\lambda }\left(t-T\right)$ for all t ∊ (T, T + T′].

The jumps of a Poisson process all have size one. A compound Poisson process is defined in such a way that the size of each of its jumps is random (see Section 16.2). Let Y1, Y2, . . . be i.i.d. random variables which are all independent of the Poisson process Nλ(t). A compound Poisson process with jump sizes Yk, k ≥ 1, is

$X\left(t\right)={\displaystyle \sum _{k=1}^{N_{\lambda }\left(t\right)}{Y}_k,t\ge 0.}\text{ (17.37)}$

The simulation of a compound Poisson process on [0, T] is straightforward. Applying Algorithm 17.14 or 17.15 gives us a sequence of occurrence times Tk, 1 ≤ k ≤ m, where m = Nλ(t) ~ Pois(λT). After that we set

$X\left(t\right)={\displaystyle \sum _{k=1}^m{\mathbb{I}}_{\left\{T_k\le t\right\}}{Y}_k.}$

17.4.5.2 Subordinated Processes

The variance gamma (VG) process is a three-parameter generalization of the Brownian motion model for the dynamics of the logarithm of the stock price. It is obtained by evaluating a scaled Brownian motion with drift at a random time given by a gamma process (see Madan and Seneta (1990)). The gamma process G(t; μ, υ) with mean rate μ > 0 and variance rate υ > 0 is a random process with independent gamma increments over nonoverlapping intervals of time.

The VG process X(t; σ, υ, θ) is defined in terms of the scaled Brownian motion with drift, B(t) ≡ W(θ,σ)(t) = θt + σW (t), and the gamma process with unit mean rate, denoted G(t) ≡ G(t; 1, υ), as

$X\left(t;\sigma ,\upsilon ,\theta \right):=B\left(G\left(t\right)\right),t\ge 0.$

The PDF of the VG process at time t can be expressed as a normal density function conditional on the realization of the gamma time change. The risk-neutral process for the asset price is given by

$S\left(t\right):=S_0\exp \left(\left(r-q-\omega \right)t+X\left(t;\sigma ,\upsilon ,\theta \right)\right),t\ge 0,\text{ (17.38)}$

where r and q are, respectively, the risk-neutral interest rate and dividend yield. The parameter ω = ln(1 − θυ − σ2υ/2)/υ is chosen so that the discounted asset price process e−(r−q)tS(t) is a true martingale.

Modelling of the variance gamma process relies on sampling from the normal and gamma probability distributions. One needs first to generate the gamma process and then to sample the Brownian motion conditional on the obtained values of the stochastic time process. Note that both the gamma process and Brownian motion are random processes with independent stationary increments. Thus, to sample a path of the variance gamma process at a discrete sequence 0 = t0 < t1 < t2 < · · · < tN, it is sufficient to generate the gamma increment G(ti) − G(ti−1) and then the Brownian increment B(gi) − B(gi−1) conditional on G(ti) = gi and G(ti−1) = gi−1 for all i = 1, 2, . . . , N. The sample values of G(ti) and X(ti) = B(G(ti)) can then be obtained by calculating cumulative sums of the respective increments. The increments of the gamma process G with mean rate one and variance rate θ and the Brownian process B with parameters θ and σ can be simulated as stated below:

G(t2) − G(t1) has the Gamma((t2 − t1)/θ,θ) distribution for any 0 < t1 < t2;

B(g2) − B(g1) has the Norm(μ(g2 − g1), σ2(g2 − g1)) distribution for any 0 < g1 < g2.

By repeating the above procedure n times and taking the cumulative sums of the increments, we can obtain the values of the variance gamma process at a discrete sequence of time moments. Algorithm 17.16 is used for sampling paths of the variance gamma process.

17.5 Variance Reduction Methods

Consider the evaluation of a mathematical expectation E[h(X)], where h is a real-valued function of a random variable X having a PDF f. We call H ≡ h(X) a direct Monte Carlo estimator as opposed to an estimator with some variance reduction techniques embedded. The direct Monte Carlo sample estimator of E[h(X)] is

${\overline{H}}_n=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}h\left(X_i\right)},$

where {Xi}i≥1 are i.i.d. random variables with the common PDF f. A direct sample estimate is the following average with statistically independent realizations x1, x2, . . . , xn of X:

${\overline{h}}_n=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}h\left(x_i\right)}.$

A large variance of the estimator h(X) results in slow convergence of the sample estimate to E[h(X)]. Any modification of the direct Monte Carlo method that results in a decrease of the variance is called a variance reduction method. One example of variance reduction techniques is the importance sampling method, which is discussed below. The goal of this section is to summarize various techniques used to improve the direct Monte Carlo method. Since a modification of the original estimator may result in an increase in computing time, we compare methods on the basis of their computational costs.

17.5.1 Numerical Integration by a Direct Monte Carlo Method

Consider a multidimensional integral $I\left(g\right)={\displaystyle {\int }_{{ℝ}^d}g\left(x\right)dx}$. If the number of dimensions d is small (say, d ≤ 3) then deterministic quadrature rules can be successfully applied to evaluate I(g). For larger dimensions, it is more beneficial to use stochastic methods. As is pointed out in the beginning of this chapter, to apply the Monte Carlo method (MCM), the quantity of interest, say Q, needs to be represented in the form of a mathematical expectation of a random variable called the estimator of Q. Select a d-variate PDF f such that f(x) ≠ 0 if g(x) ≠ 0 for all x ∊ ℝd. If the integrand g has an integrable singularity at some point x0, then the PDF f should also be singular at x0 so that ${\lim }_{x\to x_0}\frac{g\left(x\right)}{f\left(x\right)}$ exists and is finite. Moreover, it is reasonable to select f having the same support as that of g. If the integral is taken on a manifold, e.g., a sphere, then the support of the PDF f has to be the same manifold. Rewrite the integral I(g) as follows:

$I\left(g\right)={\displaystyle {\int }_{{ℝ}^d}g\left(\text{x}\right)}\text{dx}={\displaystyle {\int }_{{ℝ}^d}\frac{g\left(\text{x}\right)}{f\left(\text{x}\right)}f\left(\text{x}\right)}d\text{x}=\text{E}\left[h\left(X\right)\right],\text{ (17.39)}$

where X ~ f; $h\left(x\right):=\frac{g\left(x\right)}{f\left(x\right)}$ if g(x) ≠ 0 and g(x) ≠ ∞, h(x) := 0 if g(x) = 0, and $h\left(x_0\right):={\lim }_{x\to x_0}\frac{g\left(x\right)}{f\left(x\right)}\text{if}g\left(x_0\right)=\infty$. The integral I(g) = E[h(X)] is estimated by the Monte Carlo method as follows.

1. (1) Generate n independent sample values ${\left\{x_i\right\}}_{i=1}^n$ drawn from the PDF f.

2. (2) Construct a sample estimate of the integral I(g):

   $I\left(g\right)\approx {\overline{h}}_n=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}h\left({\text{x}}_i\right)}.$

3. (3) For 0 < α < 1, construct an asymptotically valid 100(1 − α)% confidence interval for I(g):

   $\left({\overline{h}}_n-\frac{z_{\alpha /2}{s}_n}{\sqrt{n}},{\overline{h}}_n+\frac{z_{\alpha /2}{s}_n}{\sqrt{n}}\right)∍I\left(g\right),$

   where ${\mathcal{S}}_n^2=\frac{1}{n}{{\displaystyle {\sum }_{i=1}^{n}h{\left(x_i\right)}^2-\overline{h}}}_n^2$ is a sample variance for h(X), and zα/2 is a (1 − α/2)-quantile of the standard normal distribution.

As is seen, the Monte Carlo method allows us to simultaneously construct an approximation of I(g) and an upper bound for the error:

$|I\left(g\right)-{\overline{h}}_n|\le \frac{z_{\alpha /2}\sqrt{\text{Var}\left(h\left(X\right)\right)}}{\sqrt{n}},\text{ (17.40)}$

which is valid with confidence (1 − α)100%. The variance Var(h(X)) can be approximated by the sample variance $S_n^2$. As is seen from (17.40), the approximation error is of order $\mathcal{O}\left(n^{-1/2}\right)$, as n → ∞. Since the number of sample values n is equal to the number of times the integrand g is evaluated, the Monte Carlo method can be compared with deterministic quadrature rules for which the approximation error is known in terms of n. For example, the midpoint quadrature rule provides the approximation error of order $\mathcal{O}\left(n^{-2/d}\right)$, as n → ∞. In contrast to the Monte Carlo method, the error of a quadrature rule depends the dimensionality d. The Monte Carlo method outperforms the midpoint rule for problems with d > 4. In general, the Monte Carlo method outperforms deterministic methods if d ≥ 10. For lower dimensionality (say, d ≤ 3), it is better to use a deterministic quadrature rule. For the case with 3 < d < 10, a combination of stochastic and deterministic methods may be a more efficient approach to computing the integral. The Monte Carlo method is said to beat the curse of dimensionality since its performance is not considerably affected by the dimensionality of the problem. Another advantage of the Monte Carlo method is that it can be applied to integrals of nonsmooth functions, whereas deterministic quadrature rules require the smoothness of higher-order derivatives of integrands.

The computational time required to calculate the estimate ${\overline{h}}_n$ is a product of the number of samples, n, and the time tH required to calculate one sample value of H := h(X). Let δ be a given error level. Using (17.40), we obtain

$\frac{z_{\alpha /2}\sqrt{\text{Var}\left(h\left(X\right)\right)}}{\sqrt{n}}=\delta \iff n=\frac{z_{\alpha /2}^2\text{Var}\left(h\left(X\right)\right)}{{\delta }^2}.$

Therefore, the total computational time tH · n is proportional to $\frac{t_H\text{Var}\left(h\left(X\right)\right)}{{\delta }^2}$. The product tH ·Var(h(X)) is called the computational cost of the stochastic estimator h(X). Clearly, we can accelerate the Monte Carlo method by reducing the variance of the estimator. Another approach is to parallelize computations.

17.5.2 Importance Sampling Method

Clearly, there are many choices for the PDF f in (17.39) that fit for the estimation of I(g) by the Monte Carlo method. One obvious requirement on the density f is that it should include all singularities of the integrand to guarantee that the variance of the random estimator is finite. For example, let us estimate the value of the integral

$I={\displaystyle {\int }_0^1\frac{g\left(x\right)}{\sqrt{x}}\text{d}x,\text{where}g\in C\left[0,1\right]}\text{and}g\left(0\right)\ne 0\text{ (17.41)}$

by the Monte Carlo method. Consider the following two choices for the PDF f.

Choice 1: Let $f\left(x\right)={\mathbb{I}}_{\left(0,1\right)}\left(x\right)$, i.e., X ~ Unif(0, 1). It is easy to verify that the variance of the random estimator $h\left(X\right):=\frac{g\left(X\right)}{\sqrt{X}}$ is infinite.

Choice 2: Let $f\left(x\right)\propto \frac{1}{\sqrt{x}}$ for x ∊ (0, 1). Since f integrates to one, we set $f\left(x\right)=\frac{1}{2\sqrt{x}}{\mathbb{I}}_{\left(0,1\right)}\left(x\right)$. Then, the random estimator

$h\left(X\right)=\frac{g\left(X\right)}{\sqrt{X}f\left(X\right)}=\frac{g\left(X\right)}{2}$

is bounded and hence has a finite variance.

Example 17.13.

Construct a Monte Carlo estimator of the integral

$I={\displaystyle {\int }_0^{\infty }\cdots {\displaystyle {\int }_0^{\infty }\frac{{\text{e}}^{-x_1-x_2-\cdots -x_8}}{{\left(1+x_1\cdots x_8\right)}^2}\text{d}{x}_1\cdots \text{d}{x}_8.}}$

Solution. Notice that the integrand is a product of the function ${\text{e}}^{-x_1-x_2-\cdots -x_8}$, which is a product of exponential PDFs, and the positive, bounded function (1 + x1 · · · x8)−2. Let us select the PDF

$f\left(x_1,x_2,...,x_8\right)={\text{e}}^{-x_1-x_2-\cdots -x_8}={\text{e}}^{-x_1}{\text{e}}^{-x_2}\cdots {\text{e}}^{-x_8},x_1,x_2,...,x_8>0.$

That is, the entries of the random vector X are independent exponentially distributed variables Xj ~ Exp(1), 1 ≤ j ≤ 8. Set the random estimator to be

$h\left(X\right)=\frac{1}{{\left(1+X_1{X}_2\cdots X_8\right)}^2}.$

To find an upper bound of Var(h(X)), use the following property (see Exercise 17.24). Suppose that Y is a bounded random variable such that 0 ≤ m1 ≤ Y ≤ m2 for some constants m1 and m2. Then,

$\text{Var}\left(Y\right)\le \frac{{\left(m_2-m_1\right)}^2}{4}.$

Applying this property to Y = h(X) ∊ [0, 1], obtain that $\text{Var}\left(h\left(X\right)\right)\le \frac{1}{4}$.

As is seen from the above examples, it is reasonable to select the PDF f as close to the integrand g as possible. This suggestion is confirmed by the importance sampling principle, which is presented just below. Since the computational cost is proportional to the variance of the estimator, the optimal PDF f solves the following optimization problem:

$\text{Var}\left(\frac{g\left(X\right)}{f\left(X\right)}\right)={\displaystyle {\int }_{{ℝ}^d}\frac{g^2\left(\text{x}\right)}{f\left(\text{x}\right)}\text{dx}-I^2\left(g\right)\to \underset{f}{\min }}.\text{ (17.42)}$

Let us prove that the variance attains its minimum value if f ∝ |g|. The proof is based on the Cauchy–Schwartz–Bunyakovsky inequality

${\left({\displaystyle {\int }_{{ℝ}^d}\left|u\left(\text{x}\right)\upsilon \left(\text{x}\right)\right|\text{dx}}\right)}^2\le {\displaystyle {\int }_{{ℝ}^d}{u}^2\left(\text{x}\right)\text{dx}}\cdot {\displaystyle {\int }_{{ℝ}^d}{\upsilon }^2\left(\text{x}\right)\text{dx}},\text{ (17.43)}$

for u, υ ∊ L2(ℝd) (here L2 is the set of square-integrable functions). Let us set $u\left(x\right)=\frac{g\left(x\right)}{\sqrt{f\left(x\right)}}$ and $\upsilon \left(x\right)=\sqrt{f\left(x\right)}$ to obtain

$\begin{array}{c}{\left({\displaystyle {\int }_{{ℝ}^d}|g\left(\text{x}\right)|\text{dx}}\right)}^2={\left({\displaystyle {\int }_{{ℝ}^d}\left|\frac{g\left(\text{x}\right)}{\sqrt{f\left(\text{x}\right)}}\sqrt{f\left(\text{x}\right)}\right|\text{dx}}\right)}^2\\ \le {\displaystyle {\int }_{{ℝ}^d}\frac{g^2\left(\text{x}\right)}{f\left(\text{x}\right)}\text{dx}}\cdot \underset{=1}{\underbrace{{\displaystyle {\int }_{{ℝ}^d}f\left(\text{x}\right)\text{dx}}}}={\displaystyle {\int }_{{ℝ}^d}\frac{g^2\left(\text{x}\right)}{f\left(\text{x}\right)}\text{dx}}.\end{array}$

Therefore, we obtain a lower bound of the variance

$\text{Var}\left(\frac{g\left(X\right)}{f\left(X\right)}\right)\ge {\left({\displaystyle {\int }_{{ℝ}^d}|g\left(\text{x}\right)|\text{dx}}\right)}^2-I^2\left(g\right).$

Now, it remains to show that the variance attains the lower bound when the PDF is proportional to |g|, that is, when $f\left(x\right)=\frac{1}{c}|g\left(x\right)|$ with $c={\displaystyle {\int }_{{ℝ}^d}|g\left(x\right)|\text{d}x}$. Indeed, for $f=\frac{1}{c}|g|$, we have

${\displaystyle {\int }_{{ℝ}^d}\frac{g^2\left(\text{x}\right)}{f\left(\text{x}\right)}\text{dx}}={\displaystyle {\int }_{{ℝ}^d}\frac{c{|g\left(\text{x}\right)|}^2}{|g\left(\text{x}\right)|}\text{dx}}=c{\displaystyle {\int }_{{ℝ}^d}|g\left(\text{x}\right)|\text{dx}}={\left({\displaystyle {\int }_{{ℝ}^d}|g\left(\text{x}\right)|\text{dx}}\right)}^2.$

In practice, it may be impossible to use the “best” PDF since the numerical evaluation of the normalizing constant c = I(|g|) is equivalent to the original problem. Alternatively, we can use an acceptance-rejection method that does not require the normalizing constant to be calculated. If sampling from f ∝ |g| is not feasible or is computationally expensive, then one can use any approximation PDF that is close to the optimal PDF. For example, a piecewise approximation of |g| can be used to construct the sampling density.

Example 17.14.

Approximate the integral $I={\displaystyle {\int }_0^{\pi /2}\sin x\text{d}x=1}$ by the Monte Carlo method where the sampling PDF is

1. (a) $f_1\left(x\right)=\frac{2}{\pi }{\mathbb{I}}_{\left(0,\pi /2\right)}\left(x\right)$ (a constant approximation of g);
2. (b) $f_2\left(x\right)=\frac{8x}{{\pi }^2}{\mathbb{I}}_{\left(0,\pi /2\right)}\left(x\right)$ (a linear approximation of g).

Solution. Let us compare the variances ${\sigma }_i^2=\text{Var}(H_i)$ of the random estimators $H_i=\frac{g\left(X_i\right)}{f_i\left(X_i\right)}$ with Xi ~ fi for i = 1, 2. We have

$\begin{array}{l}{\sigma }_1^2={\displaystyle {\int }_0^{\pi /2}\frac{{\sin }^{2}x}{2/\pi }}\text{d}x-1=\frac{{\pi }^2}{8}-1\cong 0.2337;\hfill \\ {\sigma }_2^2={\displaystyle {\int }_0^{\pi /2}\frac{{\sin }^{2}x}{8x/\pi }}\text{d}x-1\cong \mathrm{0.0168.}\hfill \end{array}$

By using the PDF f2 instead of f1, we can reduce the computational cost by approximately $\frac{{\sigma }_1^2}{{\sigma }_2^2}\cong 14$ times.

17.5.3 Change of Probability Measure

Let us consider another PDF $\widehat{f}$ such that $\widehat{f}\left(x\right)>0$ for all x ∊ ℝ with f(x)h(x) ≠ 0. Apply the change of measure method to obtain:

$Q=\text{E}\left[h\left(X\right)\right]=\text{E}\left[h(\widehat{X})\frac{f\left(\widehat{X}\right)}{\widehat{f}\left(\widehat{X}\right)}\right]\text{ (17.44)}$

where the last mathematical expectation in (17.44) is relative to $\widehat{X}~\widehat{f}$. The weight function $f/\widehat{f}$ is called the likelihood function. According to the importance sampling method, the optimal PDF $\widehat{f}$ is proportional to the product f ·|g|. However, sampling from $\widehat{f}\propto f\cdot \left|g\right|$ may not be feasible; therefore, one may use one of the following simple methods to “improve” the sampling density.

1. (a) Shifting the PDF: ${\widehat{f}}_c\left(x\right):=f\left(x-c\right)$ for some c ∊ ℝ. The point c may be chosen in accordance with the maximum principle so that ${\widehat{f}}_c$ and f ·|g| attain their maximum at the same point. For example, consider a normal density f which reaches its maximum value at the mean μ. In this case, set c = ν − μ, where ν = arg max{f(x) · |g(x)|}.
2. (b) Reshaping the PDF: ${\widehat{f}}_c\left(x\right):=\frac{1}{c}f\left(x/c\right)$ for some c ∊ ℝ.

17.5.4 Control Variate Method

The main idea of the control variate method is to represent the unknown quantity as a sum of two parts. One part can be calculated analytically and the other part is to be estimated by the Monte Carlo method. Such a splitting is expected to decrease the variance. For example, consider the evaluation of $I\left(g\right)={\displaystyle {\int }_{{ℝ}^d}g\left(x\right)\text{d}x}$. Suppose that there exists another function g0, which is close to g and for which the integral $I\left(g_0\right)={\displaystyle {\int }_{{ℝ}^d}{g}_0\left(x\right)\text{d}x}$ can be calculated analytically. Write I(g) = I(g0) + I(g − g0). To approximate the integral I(g − g0) we apply the Monte Carlo method with a PDF f:

$I\left(g\right)=I\left(g_0\right)+\text{E}\left[\frac{g\left(X\right)-g_0\left(X\right)}{f\left(X\right)}\right],X~f.\text{ (17.45)}$

If g0 is close to g, then the difference g − g0 is close to zero. So we expect that the Monte Carlo estimator of I(g − g0) has a smaller variance than that of the original integral I(g). Rewrite (17.45) as follows:

$I\left(g\right)=\text{E}\left[\frac{g\left(X\right)}{f\left(X\right)}-\left(\frac{g_0\left(X\right)}{f\left(X\right)}-I\left(g_0\right)\right)\right]=\text{E}\left[Y-\left(Z-I\left(g_0\right)\right)\right],\text{ (17.46)}$

where $Y:=\frac{g\left(X\right)}{f\left(X\right)}$ and $Z:=\frac{g_0\left(X\right)}{f\left(X\right)}$ are, respectively, unbiased estimators of I(g) and I(g0).

Let us apply this method to the estimation of an arbitrary mathematical expectation Q = E[Y] of some random variable Y . Suppose that there exists another random variable Z with known mathematical expectation μZ such that Y and Z can be sampled simultaneously. Construct a new parametric family of random variables:

$Y\left(b\right):=Y-b\left(Z-{\mu }_Z\right),b\in ℝ.$

Clearly, the mathematical expectation of Y(b) is equal to Q for all b ∊ ℝ:

$\text{E}\left[Y\left(b\right)\right]=\text{E}\left[Y\right]-b\left[Z-{\mu }_Z\right]=\text{E}\left[Y\right]=Q.$

Thus, for every b ∊ ℝ, Y (b) is an unbiased estimator of Q. The random variable Z is called a control variate; Y (b), b ∊ ℝ, is called a controlled estimator. The sample estimate ${\overline{y}}_n\left(b\right)$ is constructed as usual.

Proposition 17.16.

The optimal value of b chosen so that Var(Y (b)) reaches its minimum value is given by

$b_{opt}=\frac{\text{Cov}\left(Y,Z\right)}{\text{Var}\left(Z\right)}.\text{ (17.47)}$

The minimum variance of the controlled estimator is

$\text{Var}\left(Y\left(b_{opt}\right)\right)=\text{Var}\left(Y\right)\cdot \left(1-\text{Corr}{\left(Y,Z\right)}^2\right).$

Proof. The variance of the controlled estimator is a quadratic function of b:

$\text{Var}\left(Y\left(b\right)\right)=\text{Var}\left(Z\right)b^2-2\text{Cov}\left(Y,Z\right)b+\text{Var}\left(Y\right).$

Differentiate it with respect to b and equate the derivative obtained to zero:

$\frac{\text{d}\text{Var}\left(Y\left(b\right)\right)}{\text{d}b}=2\text{Var}\left(Z\right)b-2\text{Cov}\left(Y,Z\right)=0\Rightarrow b=\frac{\text{Cov}\left(Y,Z\right)}{\text{Var}\left(Z\right)}.$

Since the second derivative is positive, the variance attains its minimum value at the point bopt given by (17.47). The variance of the controlled estimator Y (b) for b = bopt is

$\begin{array}{cc}\hfill \text{Var}\left(Y\left(b_{\text{opt}}\right)\right)=& \text{Var}\left(Z\right)\frac{{\text{Cov}}^2\left(Y,Z\right)}{{\text{Var}}^2\left(Z\right)}-2\text{Cov}\left(Y,Z\right)\frac{\text{Cov}\left(Y,Z\right)}{\text{Var}\left(Z\right)}+\text{Var}\left(Y\right)\hfill \\ \hfill =& \text{Var}\left(Y\right)-\frac{{\text{Cov}}^2\left(Y,Z\right)}{\text{Var}\left(Z\right)}-\text{Var}\left(Y\right)\cdot \left(1-\frac{{\text{Cov}}^2\left(Y,Z\right)}{\text{Var}\left(Y\right)\text{Var}\left(Z\right)}\right)\hfill \\ \hfill =& \text{Var}\left(Y\right)\cdot \left(1-\text{Corr}{\left(Y,Z\right)}^2\right).\hfill \end{array}$

Thus, Var(Y (bopt)) ≤ Var(Y) since 0 ≤ Corr(Y, Z)2 ≤ 1.

In practice, Var(Z) and/or Cov(Y, Z) are unknown in closed form. So, the optimal value of b can be approximated by the ratio of the sample covariance of Y and Z to the sample variance of Z. The variance reduction factor

$\frac{\text{Var}\left(Y\right)}{\text{Var}\left(Y\left(b_{\text{opt}}\right)\right)}=\frac{1}{1-{\rho }_{YZ}^2},\text{ (17.48)}$

increases very rapidly as |ρY Z| → 1, where ρY Z = Corr(Y, Z). For example, $\frac{1}{1-{\rho }^2}=\frac{4}{3}\approx 1.3$ for $\left|\rho \right|=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, but if |ρ| = 0.99 then $\frac{1}{1-{\rho }^2}\approx 50$. This observation implies that a very high degree of correlation between the original estimator and a control variate is required to yield a substantial variance reduction.

If a nonoptimal value of b is used, then the controlled estimator may have a larger variance than the original one, as is demonstrated in the following example. Let Z be a control variate for Y so that $\text{Cov}\left(Y,Z\right)>\frac{1}{2}\text{Var}\left(Z\right)$. The controlled estimator Y (1) = Y − (Z − μZ) has a smaller variance than Y :

$\text{Var}\left(Y\left(1\right)\right)=\text{Var}\left(Y\right)-2\text{Cov}\left(Y,Z\right)+\text{Var}\left(Z\right)<\text{Var}\left(Y\right)-\text{Var}\left(Z\right)+\text{Var}\left(Z\right)=\text{Var}\left(Y\right).$

However, the variance of Y (−1) = Y + (Z − μZ) (i.e., replace Z by −Z) is larger than that of Y :

$\text{Var}\left(Y\left(-1\right)\right)=\text{Var}\left(Y\right)+2\text{Cov}\left(Y,Z\right)+\text{Var}\left(Z\right)>\text{Var}\left(Y\right)+\text{Var}\left(Z\right)+\text{Var}\left(Z\right)>\text{Var}\left(Y\right).$

Example 17.15.

Apply the control variate method to the integral

$I={\displaystyle {\int }_{{\left(0,1\right)}^d}\left({\text{e}}^{x_1+x_2+\cdots +x_d}-1\right)}\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d={\left(\text{e}-1\right)}^d-1.$

Assume that a uniform sampling distribution is used in the direct Monte Carlo method. Calculate the variance reduction factor for the optimal control variate when d = 1, 5, 10.

Solution. The PDF of the uniform distribution in the unit hypercube (0, 1)d is $f\left(x\right)={\mathbb{I}}_{{\left(0,1\right)}^d}\left(x\right)$. The direct estimator of the integral is Y = eU1+U2+···+Ud −1, where U1, U2, . . . , Ud are i.i.d. Unif(0, 1)-distributed random variables. Applying Taylor’s formula to the integrand gives us the following approximation: ${\text{e}}^{x_1+x_2+\cdots +x_d}-1\approx x_1+x_2+\cdots +x_d$. Let us use Z = U1 + U2 + · · · + Ud as a control variate. The expected value of Z is

$\text{E}\left[Z\right]={\displaystyle {\int }_{{\left[0,1\right]}^d}\left(x_1+x_2+\cdots +x_d\right)}\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d=d{\displaystyle {\int }_0^{1}x\text{d}x=\frac{d}{2}.}$

So, the controlled estimator is

$Y\left(b\right)=Y-b\left(Z-\text{E}\left[Z\right]\right)={\text{e}}^{U_1+U_2+\cdots +U_d}-1-b\left(U_1+U_2+\cdots +U_d-d/2\right).$

To find the optimal value of b and the respective variance reduction, we calculate Cov(Y, Z), Var(Y), and Var(Z):

$\begin{array}{cc}\hfill \text{E}\left[YZ\right]=& \text{E}\left[({\text{e}}^{U_1+U_2+\cdots +U_d}-1)\left(U_1+U_2+\cdots +U_d\right)\right]\hfill \\ \hfill =& d{\displaystyle {\int }_0^{1}x}{\text{e}}^x\text{d}x{\left({\displaystyle {\int }_0^1{\text{e}}^y\text{d}y}\right)}^{d-1}-d{\displaystyle {\int }_0^{1}x}\text{d}x=d{\left(\text{e}-1\right)}^{d-1}-d/2,\hfill \\ \hfill \text{Cov}\left(Y,Z\right)=& \text{E}\left[YZ\right]-\text{E}\left[Y\right]\text{E}\left[Z\right]=d{\left(\text{e}-1\right)}^{d-1}-\frac{d}{2}-\left({\left(\text{e}-1\right)}^d-1\right)\cdot \frac{d}{2}\hfill \\ \hfill =& d{\left(\text{e}-1\right)}^{d-1}\cdot \frac{3-\text{e}}{2},\hfill \\ \hfill \text{Var}\left(Y\right)=& \text{E}\left[Y^2\right]-\text{E}{\left[Y\right]}^2={\left({\displaystyle {\int }_0^1{\text{e}}^{2x}\text{d}x}\right)}^d-2{\left({\displaystyle {\int }_0^1{\text{e}}^x\text{d}x}\right)}^d+1-{\left({\left(\text{e}-1\right)}^d-1\right)}^2\hfill \\ \hfill =& {\left(\frac{{\text{e}}^2-1}{2}\right)}^d-2{\left(\text{e}-1\right)}^d+1-{\left(\text{e}-1\right)}^{2d}+2{\left(\text{e}-1\right)}^d-1\hfill \\ \hfill =& {\left(\frac{{\text{e}}^2-1}{2}\right)}^d-{\left(\text{e}-1\right)}^{2d},\hfill \\ \hfill \text{Var}\left(Z\right)=& \text{Var}\left(U_1+U_2+\cdots +U_d\right)=d\text{Var}\left(U_1\right)=\frac{d}{12}.\hfill \end{array}$

Therefore, the optimal value of b is $b_{\text{opt}}=\frac{\text{Cov}\left(Y,Z\right)}{\text{Var}\left(Z\right)}=6{\left(\text{e}-1\right)}^{d-1}\left(3-\text{e}\right)$. The results for d = 1, 5, 10 are summarized in Table 17.7. As is seen from the table, the efficiency of this control variate method is decreasing with the growth of d.

The above example suggests a universal approach to the construction of a control variate. Consider the estimation of Q = E[Y] with Y = h(X), X ~ f. Suppose that h: ℝ → ℝ has continuous derivatives and it can be approximated by the Taylor series expansion about x0:

$h\left(x\right)\approx {\displaystyle \sum _{j=1}^k\frac{h^{\left(j\right)}\left(x_0\right)}{j!}}{\left(x-x_0\right)}^j.$

Suppose that we are able to exactly calculate the moments E [(X − x0)j], 1 ≤ j ≤ k. Then, the random variable $Z:={\displaystyle {\sum }_{j=1}^k\frac{h^{\left(j\right)}\left(x_0\right)}{j!}}{\left(X-x_0\right)}^j$ can be used as a control variate. The point x0 should be selected in a way such that | Corr(Y, Z)| is maximized.

17.5.5 Antithetic Variate

The antithetic variate method attempts to reduce the variance by exploiting a symmetry of a probability distribution. Consider the estimation of Q = E[h(U)], where U ~ Unif(0, 1). For example, $Q={\displaystyle {\int }_a^{b}g\left(x\right)dx}=\text{E}\left[h\left(U\right)\right]$ with h(U) = (b − a)g(a + (b − a)U). Since U and 1 − U have the same distribution, the estimators h(U) and h(1 − U) are equally distributed as well. The antithetic estimator is defined as $Y_{\text{anti}}=\frac{1}{2}\left(h\left(U\right)+h\left(1-U\right)\right)$. Compare the variances of the direct estimator Y = h(U) and the antithetic estimator Yanti:

$\begin{array}{c}\text{Var}\left(Y_{\text{anti}}\right)=\frac{1}{4}\left(\text{Var}\left(h\left(U\right)\right)+2\text{Cov}\left(h\left(U\right),h\left(1-U\right)\right)+\text{Var}\left(h\left(1-U\right)\right)\right)\\ =\frac{1}{2}\text{Var}\left(h\left(U\right)\right)+\frac{1}{2}\text{Cov}\left(h\left(U\right),h\left(1-U\right)\right).\end{array}$

If h(U) and h(1 − U) are negatively correlated, then $\text{Var}\left(Y_{\text{anti}}\right)\le \frac{1}{2}\text{Var}\left(Y\right)$. Since every realization of Yanti requires two evaluations of h, the computational time is doubled when we apply the antithetic variate method. However, we can still expect a reduction in the computational cost if Var(Yanti) is less than Var(Y) by half. This is the case if h is a monotone function, as is proved in the following proposition.

Proposition 17.17.

Let h be a monotone function defined on [0, 1]; let Cov(h(U), h(1−U)), where U ~ Unif(0, 1), be finite. Then, Cov(h(U), h(1 − U)) ≤ 0.

Proof. Without loss of generality, assume that the function h is nondecreasing. We need to show that

$\text{E}\left[h\left(U\right)h\left(1-U\right)\right]\le \text{E}{\left[h\left(U\right)\right]}^2.$

Since h is nondecreasing on [0, 1], the value Q lies between the values of h at the endpoints on the unit interval, i.e., h(0) ≤ Q ≤ h(1). Consider the function $f\left(y\right):={\displaystyle {\int }_0^{y}h\left(1-x\right)\text{d}x-Qy}$ defined on the unit interval, [0, 1]. The function f is zero at the endpoints, f(0) = f(1) = 0. Its derivative, f′(y) = h(1 − y) − Q, is a nonincreasing function. Since f′(0) = h(1) − Q ≥ 0 and f′(1) = h(0) − Q ≤ 0, the function f is nonnegative everywhere on [0, 1]. Therefore, ${\displaystyle {\int }_0^{1}f\left(y\right)h^{\prime }\left(y\right)\text{d}y\ge 0}$. Integration by parts gives

${\displaystyle {\int }_0^{1}f\left(y\right)h^{\prime }}\left(y\right)\text{d}y=f\left(y\right)h\left(y\right){|}_0^1-{\displaystyle {\int }_0^1{f}^{\prime }\left(y\right)h}\left(y\right)\text{d}y=-{\displaystyle {\int }_0^1{f}^{\prime }\left(y\right)h}\left(y\right)\text{d}y.$

Therefore, we have

${\displaystyle {\int }_0^1{f}^{\prime }\left(y\right)h}\left(y\right)\text{d}y={\displaystyle {\int }_0^1\left(h\left(y\right)h\left(1-y\right)-h\left(y\right)Q\right)}\text{d}y={\displaystyle {\int }_0^{1}h\left(y\right)}h\left(1-y\right)\text{d}y-Q^2\le 0.$

Hence, ${\displaystyle {\int }_0^{1}h\left(y\right)h\left(1-y\right)\text{d}y\le Q^2}$.

In summary, the antithetic variate method attempts to reduce the variance by introducing negative correlation between pairs of realizations. Here are some examples of antithetic pairs.

• U and 1 − U are both Unif(0, 1)-distributed. Corr(U, 1 − U) = −1.
• Let F be a CDF and F−1 its generalized inverse. The random variables F−1(U) and F−1(1 − U) have the same CDF F , and they are antithetic to each other (F and F−1 are both monotone functions).
• Z and −Z are both standard normal variables, and Corr(Z, −Z) = −1.
• Z and 2μ − Z, where Z ~ Norm(μ, σ2), form an antithetic pair.

17.5.6 Conditional Sampling

Consider the estimation of E[h(X)] where h is a function of a random variable X. Suppose that there exists another random variable V correlated with X so that the conditional expectation $\widehat{h}\left(V\right)=\text{E}\left[h\left(X\right)|V\right]$ can be calculated exactly. Applying the double expectation formula gives us

$\text{E}\left[h\left(X\right)\right]=\text{E}\left[\text{E}\left[h\left(X\right)|V\right]\right]=\text{E}\left[\widehat{h}\left(V\right)\right].$

As is well known from a standard course on probability theory, for any pair of square-integrable random variables (Y, V), we have Var(Y) = E[Var(Y | V)] + Var(E[Y | V]). Setting Y = h(X) gives that

$\text{Var}\left(h\left(X\right)\right)=\text{E}\left[\text{Var}\left(h\left(X\right)|V\right)\right]+\text{Var}\left(\text{E}\left[h\left(X\right)|V\right]\right)\ge \text{Var}\left(\text{E}\left[h\left(X\right)|V\right]\right).$

Therefore, using the estimator $\widehat{h}\left(V\right)=\text{E}\left[h\left(X\right)|V\right]$ instead of h(X) always leads to a variance reduction.

For this algorithm to be of practical use, the following conditions must be met.

1. (a) It is easy to generate V .
2. (b) E[h(X) | V = v] is readily computable for all v.
3. (c) E[Var(h(X) | V)] is large relative to Var(E[h(X) | V]).

Example 17.16.

Let {Xi}i≥1 be i.i.d. random variables with a common CDF F; let R be a positive-integer-valued random variable independent of all Xi. Apply the conditional sampling method to estimate

$q\left(x\right)=\mathbb{P}\left(S_R\le x\right)=\text{E}\left[{\mathbb{I}}_{\left\{S_R\le x\right\}}\right],\text{where}{S}_R={\displaystyle \sum _{i=1}^R{X}_i}.$

Solution. We can try to improve the direct Monte Carlo estimator by isolating X1:

$\begin{array}{c}\mathbb{P}\left(S_R\le x|R=r\right)=\mathbb{P}\left(X_1\le x-{\displaystyle \sum _{i=2}^r{X}_i}\right)\\ =\text{E}\left[\mathbb{P}\left(X_1\le x-{\displaystyle \sum _{i=2}^r{X}_i|X_2,X_3,...}\right)\right]=\text{E}\left[F\left(x-{\displaystyle \sum _{i=2}^r{X}_i}\right)\right],\end{array}$

where the mathematical expectation is relative to X2, X3, . . . . Therefore,

$q\left(x\right)=\text{E}\left[\mathbb{P}\left(S_R\le x|R\right)\right]=\text{E}\left[\text{E}\left[F\left(x-{\displaystyle \sum _{i=2}^R{X}_i}\right)\right]\right],$

where the external expectation is relative to R and the internal expectation is taken with respect to X2, X3, . . .. This method is efficient if ℙ(R = 1) is large.

17.5.7 Stratified Sampling

Stratified sampling is a special case of conditional sampling. Consider again the estimation of Q = E[h(X)]. We construct V such that Q = E[E[h(X) | V]] by stratifying the support of X as follows. Introduce mutually disjoint subsets A1, A2, . . . , Ak so that ℙ(X ∊ [Ui≥1Ai) = 1. Define a discrete random variable V ∊ {1, 2, . . . , k} so that V = i if X ∊ Ai, 1 ≤ i ≤ k. Set pi = ℙ(V = i) = ℙ(X ∊ Ai), 1 ≤ i ≤ k. The quantity of interest Q is given as a double expectation:

$Q=\text{E}\left[h\left(X\right)\right]={\displaystyle \sum _{i=1}^k{p}_i}\text{E}\left[h\left(X\right)|X\in A_i\right]=\text{E}\left[\text{E}\left[h\left(X\right)|V\right]\right].$

The stratified sample estimator of Q is

${\overline{H}}_n^{\text{start}}={\displaystyle \sum _{i=1}^k\frac{p_i}{n_i}{\displaystyle \sum _{j=1}^{n_i}h\left(X_{i,j}\right)},}$

where $X_{i,1,...,}{X}_{i,n_i}$ are i.i.d. samples of X conditional on V = i, for 1 ≤ i ≤ k. The total sample size is n = n1 +n2 +· · ·+nk. A crucial assumption is that sampling of X conditional on V is possible. The variance of the stratified estimator is

$\text{Var}\left({\overline{H}}_n^{\text{start}}\right)={\displaystyle \sum _{i=1}^k\frac{p_i^2}{n_i^2}{\displaystyle \sum _{j=1}^{n_i}{\sigma }_1^2}}={\displaystyle \sum _{i=1}^k\frac{p_i^2{\sigma }_i^2}{n_i}},\text{ (17.49)}$

where ${\sigma }_i^2=\text{Var}(h(x)|V=i)$.

The most important question is how to select optimal allocation. The following two results provide a simple method and an optimal method for choosing the sample sizes ni, 1 ≤ i ≤ k, respectively.

Proposition 17.18.

Let the sample sizes ni be proportional to the probabilities pi, i.e., ni = n · pi, for i = 1, 2, . . . , k. Then, the variance of the stratified sample estimator ${\overline{H}}_n^{strat}$ does not exceed that of the direct sample estimator ${\overline{H}}_n^{direct}$:

$\text{Var}\left({\displaystyle \sum _{i=1}^k\frac{p_i}{n_i}{\displaystyle \sum _{j=1}^{n_i}h}\left(X_{i,j}\right)}\right)\le \text{Var}\left(\frac{1}{n}{\displaystyle \sum _{i=1}^k{\displaystyle \sum _{j=1}^{n_i}h}\left(X_{i,j}\right)}\right),$

where n = n1 + n2 + · · · + nk.

Proof.

$\begin{array}{c}n\cdot \text{Var}\left({\overline{H}}_n^{\text{strat}}\right)=n\cdot {\displaystyle \sum _{i=1}^k\frac{p_i^2{\sigma }_i^2}{n_i}={\displaystyle \sum _{i=1}^k{p}_i{\sigma }_i^2\left(n{p}_i=n_i\Rightarrow n{p}_i^2/n_i=p_i\right)}}\\ =\text{E}\left[\text{Var}\left(h\left(X\right)|V\right)\right]\le \text{Var}\left(h\left(X\right)\right)=n\cdot \text{Var}\left({\overline{H}}_n^{\text{direct}}\right).\end{array}$

Proposition 17.19.

The optimal allocations are $n_i^{*}=n\frac{p_i{\sigma }_i}{{\displaystyle {\sum }_{j=1}^k{p}_j{\sigma }_j}}$, which give

$\text{Var}\left({\overline{H}}_n^{\text{start}}\right)=\frac{1}{n}{\left({\displaystyle \sum _{i=1}^k{p}_i{\sigma }_i}\right)}^2.$

Proof. The proof is left as an exercise for the reader (see Exercise 17.35).

The standard deviations σi, 1 ≤ i ≤ k, are usually unknown. In practice, one may estimate them from “pilot” runs and then estimate the optimal sample sizes $n_1^{*},n_2^{*},.\; .\; .,n_k^{*}$.

A typical problem where the stratified sampling method is efficient is numerical integration. Consider an integral of some function g on D ⊆ ℝd. Introduce a partition of D into k disjoint subdomains:

$D={\displaystyle \underset{i=1}{\overset{k}{\cup }}{D}_i},D_i\cap D_j=\varphi \text{for}i\ne j.$

Apply the stratified sampling method to the integral to obtain

$\begin{array}{c}{\displaystyle {\int }_{D}g\left(\text{x}\right)}\text{dx}={\displaystyle \sum _{i=1}^k{\displaystyle {\int }_{D_i}g\left(\text{x}\right)}}\text{dx}\text{=}{\displaystyle \sum _{i=1}^k{\displaystyle {\int }_{D_i}h\left(\text{x}\right)}}{f}_X\left(\text{x}\right)\text{dx}\left(\text{where}h\left(\text{x}\right):=\frac{g\left(\text{x}\right)}{f_X\left(\text{x}\right)}\right)\\ ={\displaystyle \sum _{i=1}^k{p}_i{\displaystyle {\int }_{D_i}h\left(\text{x}\right)}}{f}_i\left(\text{x}\right)\text{dx}\text{=E}\left[\text{E}\left[h\left(X\right)|V\right]\right],\end{array}$

where V = i with probability $p_i={\displaystyle {\int }_{D_i}f\left(x\right)\text{d}x}$ for 1 ≤ i ≤ k, and $f_i\left(x\right):=\frac{f_X\left(x\right)}{p_i}{\mathbb{I}}_{D_i}\left(x\right)$ is the conditional density of X given V = i for 1 ≤ i ≤ k.

Example 17.17.

Construct a direct sample estimator and a stratified sample estimator with two equal subintervals for the integral $I={\displaystyle {\int }_0^1{\text{e}}^x}\text{d}x$. Use 10 sample points for both estimators. Compare the variances of the estimators obtained.

Solution. Let us construct an estimator using a uniform sampling distribution, i.e., I = E[eU] with U ~ Unif(0, 1). The direct Monte Carlo estimator with n = 10 sample points is

${\overline{H}}_{10}^{\text{direct}}=\frac{1}{10}\left({\text{e}}^{U_1}+\cdots +{\text{e}}^{U_{10}}\right),$

where U1, . . . , U10 are i.i.d. Unif(0, 1)-distributed random variables. The variance of the direct MC estimator is

$\text{Var}\left({\overline{H}}_{10}^{\text{direct}}\right)=\frac{1}{10}\text{Var}\left({\text{e}}^{U_1}\right)=\frac{1}{10}\left({\displaystyle {\int }_0^1{\text{e}}^{2x}\text{d}x-I^2}\right)=\frac{4\text{e}-{\text{e}}^2-3}{20}\cong \mathrm{0.2421.}$

Divide the integration interval into two subintervals (k = 2): (0, 1) = (0,1/2) ∪ [1/2, 1). Now for each interval we need to find the conditional distribution of U and the probability, as is discussed above. The distribution of U ~ Unif(0, 1) conditional on {U ∊ (0, 1/2)} (or {U ∊ [1/2, 1)}) is uniform:

$\left(U|\left\{U\in \left(0,\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)\right\}\right)\stackrel{d}{=}\frac{U}{2}~Unif\left(0,\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right),\left(U|\left\{U\in \left[\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,1\right)\right\}\right)\stackrel{d}{=}\frac{U+1}{2}~Unif\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,1\right).$

The probabilities are p1 = ℙ(U ∊ (0, 1/2)) = 1/2 and p2 = ℙ(U ∊ [1/2, 1)) = 1/2. Let n1 and n2 be sample sizes so that n1 + n2 = 10. The stratified sample estimator is

${\overline{H}}_{10}^{\text{start}}=\frac{1}{2n_1}\left({\text{e}}^{U_1/2}+\cdots +{\text{e}}^{U_{n_1}/2}\right)+\frac{1}{2n_2}\left({\text{e}}^{\left(V_1+1\right)/2}+\cdots +{\text{e}}^{\left(V_{n_2}+1\right)/2}\right),$

where $U_1,...,U_{n_1}$ and $V_1,...,V_{n_2}$ are i.i.d. U(0, 1)-distributed random variables. The variance of the stratified estimator is

$\begin{array}{c}\text{Var}\left({\overline{H}}_{10}^{\text{start}}\right)=\frac{1}{4n_1}\text{Var}\left({\text{e}}^{U_1/2}\right)+\frac{\text{e}}{4n_2}\text{Var}\left({\text{e}}^{V_1/2}\right)\\ =\left(\frac{1}{4n_1}+\frac{\text{e}}{4n_2}\right)\text{Var}\left({\text{e}}^{U_1/2}\right)=\left(\frac{1}{4n_1}+\frac{\text{e}}{4n_2}\right)\left({\displaystyle {\int }_0^1{\text{e}}^x\text{d}x-}{\left({\displaystyle {\int }_0^1{\text{e}}^{x/2}\text{d}x}\right)}^2\right)\\ =\left(\frac{1}{4n_1}+\frac{\text{e}}{4n_2}\right)\left(8\sqrt{\text{e}}-3\text{e}-5\right)\cong \frac{0.008731}{n_1}+\frac{0.02373}{n_2}.\end{array}$

If $n_1=n_2=5$, then Var $\left({\overline{H}}_{10}^{\text{strat}}\right)\cong 0.006493<0.02421\cong \text{Var}\left({\overline{H}}_{10}^{\text{direct}}\right)$.

17.6 Exercises

1. Exercise 17.1. In the “hit-or-miss” method, the value of π is estimated as follows. A point is sampled uniformly in the square [−1, 1] × [−1, 1], i.e., sample two independent Cartesian coordinates X, Y ~ Unif(−1, 1). A point (X, Y) is accepted if it lies inside the circle inscribed in the square, i.e., X2 + Y2 < 1. The experiment is repeated N times, and the number of accepted points, NH, is recorded. The ratio $\frac{N_H}{N}$ converges to a limiting value whose expression involves π.
   1. (a) Construct a random estimator of π based on the experiment described above.
   2. (b) Construct a sample estimate of π with N trials. Express it in terms of the ratio $\frac{N_H}{N}$.
   3. (c) Construct a 99% confidence interval for π. Find the number of trials, N, required to obtain the confidence interval of length less than 10−3.

2. Exercise 17.2. Let U ~ Unif(0, 1) admit the following binary representation:

   $U={\left(0.B_1{B}_2...B_k...\right)}_2={\displaystyle \sum _{k=1}^{\infty }{B}_k\cdot 2^{-k}}.$

   Show that the digits Bk, k = 1, 2, . . ., are independent Bernoulli random variables uniformly distributed in {0, 1}.

3. Exercise 17.3. Propose an algorithm for simulating the occurrence of two dependent events A and B that uses only one uniform random variable U ~ Unif(0, 1). Assume that the probabilities ℙ(A) = pA, ℙ(B) = pB, and ℙ(A ∪ B) = q are given.
4. Exercise 17.4. Using the inverse CDF method, find generating formulae for the following PDFs:
   1. (a) $f\left(x\right)=\frac{3e^{-x}}{\left(e^{-x}+3\right)\sqrt{9-e^{-2x}}},-\ln 3<x<+\infty$;
   2. (b) $f\left(x\right)=\frac{2\cos x}{{\sin }^{3}x},\frac{\pi }{4}\le x\le \frac{\pi }{2}$;
   3. (c) $f\left(x\right)=\frac{8x}{{\left(x+1\right)}^3},0\le x\le 1$;
   4. (d) $f\left(x\right)=x\sin x^2,0\le x\le \sqrt{\pi }$;
   5. (e) $f\left(x\right)=\frac{3}{2{\left(2x+1\right)}^{3/2}},0\le x\le 4$;
   6. (f) $f\left(x\right)=\frac{1}{{\cos }^{2}x},0\le x\le \frac{\pi }{4}$.
5. Exercise 17.5. Using the inverse CDF method, find generating formulae for the Rayleigh distribution with the CDF F (x) = 1 − e−2x(x−b), x > b.
6. Exercise 17.6. Let X be an absolutely continuous random variable such that the inverse function of its distribution function F is well-defined. Consider the problem of sampling X conditional on a < X < b, with F (a) < F (b).
   1. (a) Find the CDF for such a conditional distribution.
   2. (b) Show that the generating formula X = F−1(F (a) + (F(b) − F(a))U) with U ~ Unif(0, 1) produces the desired conditional distribution.

7. Exercise 17.7. Justify the following sampling formula for the geometric random variable X with parameter p ∊ (0, 1) (i.e., ℙ(X = k) = (1 − p)k−1p for k = 1, 2, . . .):

   $X=\lfloor \frac{\ln \left(1-U\right)}{\ln \left(1-p\right)}\rfloor +1,U\in Unif\left(0,1\right).$

   [Hint: Evaluate the probability ℙ(X = k), k = 1, 2, . . ..]

8. Exercise 17.8. Consider a random variable with the piecewise-constant PDF

   $f\left(x\right)=\{\begin{array}{l}\frac{1}{2},0<x\le 1,\hfill \\ \frac{1}{8},1<x\le 3,\hfill \\ \frac{1}{12},3<x<6.\hfill \end{array}$

   Design a simulation algorithm using the inverse-transform method.

9. Exercise 17.9. Consider a random variable with a PDF f, which is proportional to the piecewise function

   $g\left(x\right)=\{\begin{array}{l}\frac{x}{3},0<x\le 2,\hfill \\ \frac{1}{3},2<x\le 4,\hfill \\ 1-\frac{x}{6},4<x<6.\hfill \end{array}$

   Find the constant c such that f(x) = cg(x) is a probability density. Design a simulation algorithm using the composition method.

10. Exercise 17.10. Present an acceptance-rejection method for a random variable on (0, 1) with the PDF of the form

    $f\left(x\right)=g\left(x\right)/x^a,0<g\left(x\right)\le M,0<a<1.$

11. Exercise 17.11. Construct the Neumann method (and estimate its computational cost) for sampling a random variable X whose density is proportional to

    1. (a) $f\left(u\right)=3-\sqrt[3]{2u},0\le u\le 1$;
    2. (b) f(u) = −u2 + 2 u + 3 , 0 < u < 3;
    3. (c) f(u) = u5/3 (1 − u)3/2, 0 ≤ u ≤ 1;
    4. (d) f(u) = u5/3 e−u, u > 0.

    For each case, compute the probability of acceptance.

12. Exercise 17.12. Consider the triangular probability distribution with the PDF

    $f\left(x\right)=\frac{x}{2}{\mathbb{I}}_{\left(0,1\right]}\left(x\right)+\frac{4-x}{6}{\mathbb{I}}_{\left(1,4\right)}\left(x\right).$

    1. (a) Obtain the CDF F and then its inverse F−1. Describe the inverse CDF sampling method.
    2. (b) Develop the decomposition method using the strata {(0, 1], (1, 4)}.

    3. (c) Determine which of the following majorizing functions provides the largest acceptance probability in the acceptance-rejection method:

       $g_1\left(x\right)=\frac{x}{2},g_2\left(x\right)=\frac{4-x}{6},g_3\left(x\right)=\frac{1}{2},x\in \left(0,4\right).$

       Find the acceptance probability for each case.

13. Exercise 17.13. Let X be a positive random variable with the probability density function $f\left(x\right)=\sqrt{\frac{2}{\pi {\sigma }^2}}{e}^{-x^2/\left(2{\sigma }^2\right)}$, i.e., X = |Z| where Z ~ Norm(0, σ2). Develop an acceptance-rejection algorithm for generating a random variable from the PDF f using the exponential distribution Exp(λ) as the majorizing distribution. Which λ gives the largest acceptance probability (give the answer in terms of σ)?

14. Exercise 17.14. Consider a bivariate distribution with the joint density

    $f\left(x,y\right)=\{\begin{array}{l}cxy\text{if}x>0,y<1,y-x>0,\hfill \\ 0\text{otherwise}.\hfill \end{array}$

    1. (a) Find the constant c.
    2. (b) Represent the joint density as a product of marginal and conditional univariate densities.
    3. (c) Construct the exact simulation algorithm based on the inverse CDF method. Consider the following approach: first model X and then Y.

15. Exercise 17.15. The density of a random variable X is represented in the integral form:

    $f_X\left(u\right)=c{\displaystyle {\int }_1^{\infty }{\upsilon }^{-c}{\text{e}}^{-u\upsilon }}\text{d}\upsilon ,u>0,c>0.$

    1. (a) Show that fX is a density function.

    2. (b) For the joint density

       $f_{X,Y}\left(u,\upsilon \right)=c{\upsilon }^{-c}{\text{e}}^{-u\upsilon },u>0,\upsilon >1,$

       find a representation of the form fX,Y (u, v) = fY (υ)fX|Y (u|v).

    3. (c) Based on the above representation of the marginal density

       $f_X\left(u\right)={\displaystyle {\int }_{-\infty }^{\infty }{f}_Y\left(\upsilon \right)}{f}_{X|Y}\left(u|\upsilon \right)\text{d}\upsilon ,$

       construct an algorithm for generating X, first by sampling Y from fY (υ), and then by sampling from the conditional density fX|Y (u|v).

16. Exercise 17.16. Demonstrate how three independent tosses of a balanced coin can be modelled by two rolls of a balanced die (with six faces). In other words, develop an exact algorithm for sampling a vector [X1, X2, X3] formed by three independent discrete random variables uniformly distributed in {0, 1} by using two independent discrete random variables Y1 and Y2 uniformly distributed in {1, 2, 3, 4, 5, 6}.
17. Exercise 17.17. Develop an algorithm that allows you to “replace” a fair roulette wheel (with 37 pockets numbered from 0 to 36) by a balanced die (with six faces labelled from 1 to 6). In other words, develop an exact algorithm for sampling a discrete random variable X uniformly distributed in {0, 1, 2, . . . , 36} by using a sequence (Yk)k≥0 of i.i.d. discrete random variables uniformly distributed in {1, 2, 3, 4, 5, 6}.

18. Exercise 17.81. Consider acceptance-rejection sampling from a standard normal distribution.

    1. (a) Find the range for λ > 0 and c > 0 so that the function

       $g\left(x\right)=\{\begin{array}{l}1\text{if}\left|x\right|<c,\hfill \\ {\text{e}}^{-\lambda \left(\left|x\right|-c\right)}\text{if}\left|x\right|\ge c,\hfill \end{array}$

       majorizes $e^{-x^2/2}\text{on}\left(-\infty ,\infty \right)$.

    2. (b) Calculate the acceptance probability (as a function of c and λ).
    3. (c) Let us set $c=1/\sqrt{2}$ and $\lambda =2c=\sqrt{2}$. Show that this choice of parameters maximizes the acceptance probability equal to $\frac{\sqrt{\pi }}{2}$.

    4. (d) Show that the proposal CDF G is given by

       $G\left(x\right)=\{\begin{array}{l}\frac{1}{4}{\text{e}}^{\sqrt{2x}+1}\text{if}x<-\frac{\sqrt{2}}{2},\hfill \\ \frac{\sqrt{2x}+2}{4}\text{if}\left|x\right|\le \frac{\sqrt{2}}{2},\hfill \\ 1-\frac{1}{4}{\text{e}}^{-\sqrt{2x}+1}\text{if}x>\frac{\sqrt{2}}{2}.\hfill \end{array}$

    5. (e) Develop computational algorithms

       1. (i) for sampling from the inverse of G;
       2. (ii) for sampling from the standard normal distribution by using the acceptance-rejection method with the proposal CDF G.

19. Exercise 17.19. Construct two methods of sampling from the trivariate normal distribution

    $Nor{m}_3\left(\left[\begin{array}{c}1\\ 0\\ -1\end{array}\right],\left[\begin{array}{ccc}1& 1& 1\\ 1& 4& 1\\ 1& 1& 9\end{array}\right]\right).$

    using the Cholesky factorization and the conditional sampling approach, respectively.

20. Exercise 17.20. By applying the conditioning formula for a multivariate normal distribution, find the conditional (bridge) distribution of a Brownian motion with drift

    $B\left(t\right)=\mu t+\sigma W\left(t\right),t\ge 0.$

    In particular, show that for every 0 ≤ s < u < t, B(u) conditional on values of B(s) and B(t) has a normal distribution. Find the mean and variance of this conditional distribution.

21. Exercise 17.21. Consider a time-homogeneous random process {X(t)}t≥0 with the transition PDF p(t; y, x) defined by

    $p\left(t;y,x\right)\text{d}x=\mathbb{P}\left(X\left(t+s\right)\in \left[x,x+\text{d}x\right)|X\left(s\right)=y\right)\text{for}t,s>0$

    with infinitesimally small dx. Using the notion of conditional probability, show that for every 0 ≤ t1 < t < t2 the bridge density b(t; x|t1, t2; x1, x2) of X(t) conditional on X(t1) = x1 and X(t2) = x2 is given by

    $b\left(t;x|t_1,t_2;x_1,x_2\right)=\frac{p\left(t-t_1;x_1,x\right)p\left(t_1-t;x,x_2\right)}{p\left(t_2-t_1;x_1,x_2\right)}.$

22. Exercise 17.22. To simplify the Euler scheme, the Brownian increments ΔW = W (t + h) − W (t) can be replaced by other random variables $\widehat{\Delta W}$ with moments up to order 5 that are within $\mathcal{O}(h^3)$ of those of ΔW. Find the moments of the three point distribution

    $P(\widehat{△W}=\pm \sqrt{3h})=\frac{1}{6},P(\widehat{△W}=0)=\frac{2}{3}$

    and compare them with those of ΔW.

23. Exercise 17.23. Solve the problem of minimizing the mean-square error subject to a computational budget s as follows:

    $\frac{C_2}{n}+C_1^2{h}^{2\beta }\to \min \text{subject}\text{to}\frac{n{C}_3}{h}=s,$

    where h is the discretization parameter, n is the sample size, and C1, C2, C3 are positive constants. Find the optimal values of h and n. Find the order (w.r.t. s) of the optimal MSE.

24. Exercise 17.24. Let a random variable Y satisfy 0 ≤ m1 ≤ Y ≤ m2 < ∞ a.s. Show that the variance of Y satisfies

    $\text{Var}\left(Y\right)\le {\left(m_2-m_1\right)}^2/4.$

    [Hint: Consider $\text{E}\left[{\left(Y-\frac{m_1+m_2}{2}\right)}^2\right]$.]

25. Exercise 17.25. Construct an importance sampling Monte Carlo method for evaluating the integral

    $I={\displaystyle {\int }_0^{\pi }\cdots {\displaystyle {\int }_0^{\pi }\cos }\left(1+\sin \left(u_1...u_6\right)\right)}\sin u_1...\sin u_6\text{d}{u}_1\cdots \text{d}{u}_6.$

    Choose a nonconstant density function f so that it is close enough to the integrand and, on the other hand, can be easily sampled from.

26. Exercise 17.26. Suppose that the integral $I={\displaystyle {\int }_0^1{\text{e}}^x}$ dx is evaluated by the Monte Carlo method using the density f(x) ∞ (1 + cx), c > 0. Find the optimal value of c so that the variance of the standard estimator eX/f(X), X ~ f, attains its minimum value.

27. Exercise 17.27. Show the advantage of using the control variate method for evaluating the integral

    $I={\displaystyle {\int }_0^1\cdots {\displaystyle {\int }_0^1\sqrt{1+u_1}}}...\sqrt{1+u_8}\text{d}{u}_1\cdots \text{d}{u}_8.$

    Use the first two terms of a Taylor series of $\sqrt{1+u}$ to construct the control variate.

28. Exercise 17.28. Let the integral $I={\displaystyle {\int }_0^1\sin \left(\frac{\pi x}{2}\right)}$ dx be evaluated using the Monte Carlo method.

    1. (a) Suppose that a combination of stratified random sampling and an antithetic variate method of the form

       $Y=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(f\left(\raisebox{1ex}{$U$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)+f\left(1-\raisebox{1ex}{$U$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)\right),U~Unif\left(0,1\right),$

       where f(x) = sin(πx/2), is applied to evaluate I = E[Y]. What is the efficiency gain of this estimator as compared to the direct Monte Carlo estimator X = f(U), U ~ Unif(0, 1).

    2. (b) How large a sample size do you need if you use the direct Monte Carlo or the antithetic stratified method of (a), respectively, in order to estimate the above integral, correct to four decimal places (i.e., the error does not exceed 10−4/2) with a confidence level of 95%?

29. Exercise 17.29. Suppose the integral

    $I={\displaystyle {\int }_0^2{x}^m\text{d}x\text{for}m=2,3,...}$

    is approximated by the Monte Carlo method.

    1. (a) Find the variance of the direct Monte Carlo estimator X = 2f(U), U ~ Unif(0, 2), where f(x) = xm.
    2. (b) Suppose that a combination of stratified random sampling and an antithetic variate method of the form Y = f(U) + f(2 − U), U ~ Unif(0, 2), is applied to evaluate the integral.
       1. (i) Show that I = E[Y].
       2. (ii) Find the variance Var(Y).
       3. (iii) What is the efficiency gain of this estimator as compared to the direct estimator?
       4. (iv) How does the efficiency gain behave with the increase of m?
    3. (c) Suppose Z = U is used as a control variate. Give the formula of the controlled estimator. Find the variance reduction factor for the optimal controlled estimator in comparison with the direct Monte Carlo estimator. How does the efficiency gain behave with the increase of m? [Note: You are not required to compute the optimal value of the control variate parameter.]
    4. (d) How large a sample size do you need if you use the direct Monte Carlo or the control variate method of (c), respectively, to approximate the above integral, correct to two decimal places (i.e., the absolute error does not exceed 10−2/2) with a confidence level of 95%?
30. Exercise 17.30. Let R ~ Geom(p), and X1, X2, . . . be i.i.d. Exp(λ) random variables. Find the distribution of $S_R={\displaystyle {\sum }_{i=1}^R{X}_i}$.

31. Exercise 17.31. Construct an importance sampling Monte Carlo method for evaluating the integral

    $I={\displaystyle {\int }_0^1\cdots {\displaystyle {\int }_0^1\ln }\left(3+u_1{u}_2\cdots u_{11}\right)}{u}_1^2{u}_2^2\cdots u_{11}^2\text{d}{u}_1\text{d}{u}_2\cdots \text{d}{u}_{11}.$

    Choose the density function f so that it is close enough to the integrand and, on the other hand, can be easily sampled from. A constant density is not acceptable!

32. Exercise 17.32. Show the advantage of using the control variate method for evaluating the integral

    $I={\displaystyle {\int }_0^{\pi /2}{\displaystyle {\int }_0^{\pi /2}\cdots {\displaystyle {\int }_0^{\pi ,2}\sin }}}{u}_1\sin u_2\cdots \sin u_5\text{d}{u}_1\text{d}{u}_2\cdots \text{d}{u}_5.$

    Use u1 u2 · · · u5 to construct a control variate.

33. Exercise 17.33. Prove that for any pair of random variables (U, V),

    $\text{Var}\left(U\right)=\text{E}\left[\text{Var}\left(U|V\right)\right]+\text{Var}\left(E\left[U|V\right]\right)$

    [Hint: Use the fact that E[U2] = E[E[U2|V]] and Var(U) = E[U2] − (E[U])2.]

34. Exercise 17.34. Let X1, X2, . . . , Xn be independent random variables with expected values E[Xi] = μi and nonzero variances. Consider the following estimator of E[Y]:

    $Z=Y+{\displaystyle \sum _{i=1}^n{c}_i}\left(X_i-{\mu }_i\right).$

    Show that the values of c1, c2, . . . , cn that minimize Var(Z) are

    $c_i=-\frac{\text{Cov}\left(Y,X_i\right)}{\text{Var}\left(X_i\right)},i=1,2,...,n.$

35. Exercise 17.35. Show that the solution to the minimization problem

    $\underset{n_1,...,n_k}{\min }{\displaystyle \sum _{i=1}^k\frac{p_i^2{\sigma }_i^2}{n_i}},\text{such}\text{that}{n}_1+n_2+\cdots +n_k=N,$

    is given by

    $n_i=N\frac{p_i{\sigma }_i}{{\Sigma }_{j=1}^k{p}_j{\sigma }_j},i=1,2,...,k.$

    [Hint: Use Lagrange multipliers.]

36. Exercise 17.36. Show that the solution of

    $\arg \underset{g}{\min }{\text{Var}}_g\left(h\left(X\right)\frac{f\left(X\right)}{g\left(X\right)}\right),X~g,$

    where g is a PDF, is given by

    $g\left(x\right)=\frac{\left|h\left(x\right)\right|f\left(x\right)}{{\displaystyle {\int }_{-\infty }^{\infty }\left|h\left(x\right)\right|f\left(x\right)\text{d}x}}.$

References

I. G. Dyadkin and G. H. Kenneth. A study of 128-bit multipliers for congruential pseudo-random number generators. Computer Physics Communications, 125(13):239–258, 2000.

S. M. Ermakov. Metod Monte-Karlo i smezhnye voprosy [The Monte Carlo Method and Related Questions]. Izdat. “Nauka,” Moscow, 1975. Second edition, augmented, Teoriya Veroyatnostei i Matematicheskaya Statistika. [Monographs in Probability and Mathematical Statistics].

S. M. Ermakov and G. A. Mikhaĭlov. Statisticheskoe modelirovanie [Statistical Modelling]. “Nauka”, Moscow, second edition, 1982.

A. Gut. An Intermediate Course in Probability. Springer Texts in Statistics. Springer, New York, second edition, 2009.

P. E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Applications of Mathematics. Springer, 2011.

P. L’Ecuyer. Random number generation. In Handbook of Computational Statistics, pages 35–71. Springer, 2012.

D. H. Lehmer. Mathematical methods in large-scale computing units. In Proceedings of a Second Symposium on Large-Scale Digital Calculating Machinery, 1949, pages 141–146, Cambridge, MA, 1951. Harvard University Press.

D. B. Madan and E. Seneta. The Variance Gamma model for share market returns. Journal of Business, 63(4):511–524, 1990.

R. N. Makarov and D. Glew. Exact simulation of Bessel diffusions. Monte Carlo Methods and Applications, 16(3):283–306, 2010.

S. K. Park and K. W. Miller. Random number generators: good ones are hard to find. Commun. ACM, 31(10):1192–1201, October 1988.

A. V. Voitishek and G. A. Mikhaĭlov. Numerical Statistical Modeling: Monte Carlo Methods. Akademiya, Moscow, 2006.

A. J. Walker. An efficient method for generating discrete random variables with general distributions. ACM Transactions on Mathematical Software, 3(3):253–256, 1977.