Comparison of Sobol’ sequences in financial applications

Following [5, 7, 22], we recall a digital method to construct QMC point sets $P$ and (infinite) sequences $S$. Sobol’ sequences are included in these classes. Let $\mathbb{F}_{2}:=\{0,1\}$ be the two-element field. We perform addition and multiplication over $\mathbb{F}_{2}$ (or modulo $2$).

The concept of digital nets can be extended to (infinite) sequences $S=\{\mathbf{x}_{0},\mathbf{x}_{1},\ldots\}\subset[0,1)^{s}$ for $\infty\times\infty$ generating matrices $C_{1},\ldots,C_{s}\in\mathbb{F}_{2}^{\infty\times\infty}$ and infinite expansions $n=\sum_{l=0}^{\infty}n_{l}2^{l}$ and $x_{n,i}:=\sum_{l=0}^{\infty}x_{n,i,l}2^{-l-1}$ that contain only a finite number of nonzero terms. The resulting sequence $S$ is called a digital sequence and the matrices $C_{1},\ldots,C_{s}$ are called the generating matrices of the digital sequence.

As a quality parameter of uniformity for a point set $P$, we recall the definition of the $t$-value. See [5, 7, 22] for details.

A point set $P$ is well distributed if the $t$-value is small. The integration error is bounded by $O(2^{t}(\log N)^{s-1}/N)$ for $N=2^{m}$ when $f$ is smooth. The factor $(\log N)^{s-1}$ is not negligible if $s$ is large, but QMC works well for high-dimensional numerical integration in finance possibly because $f(\mathbf{x})$ has low effective dimension. In the case of digital nets, the $t$-value can be easily calculated by some algorithms [6, 26].

Sobol’ [27] proposed a construction method for generating matrices $C_{1},\ldots,C_{s}\in\mathbb{F}_{2}^{\infty\times\infty}$ that have a good structure of $(t,m,s)$-nets. His sequences are now called Sobol’ sequences and are included in a subclass of generalized Niederreiter sequences [31, 32]. Recently, Faure and Lemieux [8] described the relationships between them in detail. From this viewpoint, Sobol’ sequences are formulated as follows:

1. 1.

   Let $p_{1}(x)=x\in\mathbb{F}_{2}[x]$ and $p_{i}(x)\in\mathbb{F}_{2}[x]$, $2\leq i\leq s$, be the $(i-1)$th primitive polynomials in a list of primitive polynomials that are sorted in non-decreasing order of degree, i.e., $p_{2}(x)=x+1,p_{3}(x)=x^{2}+x+1$, and so on. Let $e_{i}:=\deg(p_{i})$.

2. 2.

   For each $1\leq i\leq s$, set polynomials $g_{i,0}(x),\ldots,g_{i,{e_{i}-1}}(x)\in\mathbb{F}_{2}[x]$ such that

   $\displaystyle\deg g_{i,k}(x)=e_{i}-1-k$

   (3)

   for $0\leq k\leq e_{i}-1$, in advance. These polynomials are the parameters that can be freely chosen by the user and correspond one-to one to the so-called direction numbers (see Remark 2.4 for details).

3. 3.

   For $u=1,2,\ldots$, consider the formal power series expansion

   $\displaystyle\frac{g_{i,k}(x)}{p_{i}(x)^{u}}=\sum_{v=1}^{\infty}a^{(i)}(u,k,v)x^{-v}\in\mathbb{F}_{2}((x^{-1})).$

   (4)

4. 4.

   Define $C_{i}=(c_{j,v}^{(i)})_{j\geq 1,v\geq 1}\in\mathbb{F}_{2}^{\infty\times\infty}$ as $c_{j,v}^{(i)}=a^{(i)}(Q+1,k,v)\in\mathbb{F}_{2}$ for $1\leq i\leq s,j\geq 1,v\geq 1$, where

   $\displaystyle j-1=Qe_{i}+k,$

   (5)

   with integers $Q=Q(i,j)$ and $k=k(i,j)$ satisfying $0\leq k\leq e_{i}-1$. Note that each row of $C_{i}$ corresponds to each formal power series expansion in (4). Note that the conditions (3) and (5) correspond to the reordering of rows of $C_{i}$ so as to obtain non-singular upper triangular (NUT) matrices.

The first $2^{m}$ points $P$ can be viewed as a digital net generated by the upper-left $m\times m$ submatrices of $C_{1},\ldots,C_{s}\in\mathbb{F}_{2}^{\infty\times\infty}$. We can prove that $P$ is a $(t,m,s)$-net with the following properties (see [7, 22]):

• •

  Each one-dimensional (1D) projection is a $(0,m,1)$-net, which means that 1D projections have already been optimized, that is, each $t$-value is $0$.

• •

  The $t$-value is $\leq\sum_{i=1}^{s}(e_{i}-1)$ for any $m$, which means that the initial dimensions have already been optimized.

• •

  For any low-dimensional projection, the $t$-value is $\leq\sum(e_{i}-1)$, where $\sum$ is taken over the corresponding projections.

The condition (3) can be described as

for $0\leq k\leq e_{i}-1$, and a good selection of lower terms makes us obtain $t$-values smaller than those of the above upper bounds.

Assume that under the risk-neutral measure the asset price $S_{t}$ follows the Black–Scholes model (i.e., geometric Brownian motion):

where $r$ is the risk-free interest rate, $\sigma$ is the volatility, $B_{t}$ is a standard Brownian motion. The problem of pricing an Asian call option on the discrete arithmetic average is formulated as follows: the payoff function is given by $\max(0,\frac{1}{s}\sum_{i=1}^{s}S_{t_{i}}-K)$, where $K$ is the strike price at maturity $T$, and a time interval $[0,T]$ is discretized at equally spaced times $t_{i}=i\Delta t$ for $i=1,\ldots,s$, where $\Delta t=T/s$. Then, the value of the option at time 0 is given by

The analytical solution to (7) is given by $S_{t}=S_{0}\exp((r-\sigma^{2}/2)t+\sigma B_{t})$, so it is sufficient to simulate sample paths of Brownian motion. The standard construction of Brownian motion is to generate $B_{t_{i}}$ sequentially in time: given $B_{0}=0$,

where $Z_{1},\ldots,Z_{s}\sim N(0,1)$ are i.i.d. standard normally distributed random variables. The standard construction (9) can be written as

where $A$ is an $s\times s$ lower triangular matrix. Thus, the expectation (8) can be written as

where $\Phi:(0,1)\to\mathbb{R}$ denotes the cumulative distribution function of the standard normal distribution and $(w_{1},\ldots,w_{s})^{\top}:=A(\Phi^{-1}(x_{1}),\ldots,\Phi^{-1}(x_{s}))^{\top}$ for $\mathbf{x}=(x_{1},\ldots,x_{s})\in(0,1)^{s}$. We use the following parameters: $T=1,r=0.1,\sigma=0.2,S_{0}=100,K=100$, which were used in [34].

First, we consider the case of dimension $s=360$. Figure 1 shows a summary of the standard error $\mbox{stderr}(\bar{Q})$ in $\log_{2}$ scale for $m=1,\ldots,20$. In our experiments, we applied QMC methods based on Sobol’ and Niederreiter sequences (a)–(d) and crude MC methods using random number sequences from Mersenne Twister [19]. For this, we observed that Niederreiter (NUT) and Sobol’ (JoeKuo08) are more effective than the others, particularly for $m=10,\ldots,18$, which are often used in practice. Additionally, our result seems to agree with the frequency of $t$-values for 2D projections in Table 1. Thus, it is inferred that the Asian option using the standard construction (9) has low effective dimension in the superposition sense but high effective dimension in the truncation sense.

Further, we consider the higher dimensional case $s=1024$. Figure 2 shows a summary of $\mbox{stderr}(\bar{Q})$ in $\log_{2}$ scale. In the standard construction (9), Niederreiter (NUT) and Sobol’ (JoeKuo08) also provide better results than Sobol’ (JoeKuo03). Further, we recall dimension reduction techniques, such as the principal component analysis (PCA) construction [1] for generating Brownian motion $B_{t}$, which enhance the efficiency of QMC methods. Here, the sampled Brownian motion $(B_{t_{1}},\ldots,B_{t_{s}})^{\top}$ is normally distributed with mean $\bf{0}$ and covariance matrix $C=(\min(t_{i},t_{j}))_{i,j=1}^{s}$, i.e., $(B_{t_{1}},\ldots,B_{t_{s}})^{{\top}}\sim N(\mathbf{0},C)$. Generally, we obtain the equivalent paths of Brownian motion:

provided we apply the change of variables $\mathbf{x}=A\mathbf{z}$ with $AA^{\top}=C$. The matrix $A$ in the standard construction (10) is the Cholesky matrix of $C$, i.e., $AA^{\top}=C$. Conversely, the PCA construction is a method used to choose $A=[\sqrt{\lambda_{1}}\mathbf{v}_{1},\cdots,\sqrt{\lambda_{s}}\mathbf{v}_{s}]$, where $\lambda_{1}\geq\cdots\geq\lambda_{s}$ are the eigenvalues and $\mathbf{v}_{1},\ldots,\mathbf{v}_{s}$ are the corresponding unit-length eigenvectors of $C$. In our experiment, the Niederreiter (NUT) and Sobol’ sequences with PCA outperform those with the standard construction, but those with PCA have exactly the same convergence rates. Our result implies that PCA transforms the integrand so as to have low effective dimension in the truncation sense, that is, the important variables are concentrated in the first few dimensions (e.g., $\leq 2$ or $3$). However, for Sobol’ sequences, the lower terms of $g_{i,k}(x)$ in (6) have almost no choice and are fixed in these first dimensions because the degree $e_{i}$ in (6) is sufficiently small. Thus, it seems to be difficult to expect further improvement for dimension reduction techniques as a result of changing the direction numbers for Sobol’ sequences.

We also tested the Brownian bridge (BB) construction [20] as another dimension reduction technique and observed that there is no difference among Niederreiter (NUT) and Sobol’ sequences for the convergence rates, which are better than the standard construction but worse than PCA, so we omitted the results.

Assume that the asset price $S_{t}$ follows the Black–Scholes model (7). Papageorgiou [24] considered the following digital option:

where $(x)_{+}^{0}$ is equal to 1 if $x>0$ and is $0$ otherwise, $x\in\mathbb{R}$. He indicated that effective dimension reduction techniques perform worse than the standard construction (9). Wang and Tan [36] and Wang [33] found that if the paths are generated by the standard construction, then the discontinuities of the payoff function of the sum of the indicator functions are aligned with the coordinate axes, so good performance is expected, but BB and PCA do not have this type of discontinuity. Thus, the standard construction is a good choice in this QMC setting. Figure 3 shows a summary of $\mbox{stderr}(\bar{Q})$ in $\log_{2}$ scale. We used the parameters $s=128,T=1,r=0.045,\sigma=0.3,S_{0}=100$ from [24]. Indeed, PCA is worse than the standard construction. Niederreiter (NUT) and Sobol’ (JoeKuo08) with good 2D projections are useful for such a problem. Note that this example is very simple and the value (11) can be calculated analytically.

Following [34, 35], under the risk-neutral measure, we consider a European-style basket call option on the arithmetic average over $s$ assets $S_{t}^{(1)},\ldots,S_{t}^{(s)}$, and assume that each asset satisfies

for a mean return parameter $r$ and volatility parameters $\sigma^{(i)}$. Assume $B_{t}^{(1)},\ldots,B_{t}^{(s)}$ are correlated Brownian motions with correlations $\rho_{ij}$, and the terminal pay off at $T$ is given by $\max(0,\frac{1}{s}\sum_{i=1}^{s}S_{T}^{(i)}-K)$. For this, we compute the price of the basket option:

Note that $s$ is the number of assets, not the number of discretization steps. The solutions to (12) are given by $S_{t}^{(i)}=S_{0}^{(i)}\exp((r-(\sigma^{(i)})^{2}/2)t+\sigma^{(i)}B_{t}^{(i)})$. Here, the random vector $(B_{T}^{(1)},\ldots,B_{T}^{(s)})^{\top}$ is normally distributed with mean $\mathbf{0}$ and covariance matrix $C=(\rho_{ij}T)_{i,j=1}^{s}$. Let $(Z_{1},\ldots,Z_{s})^{\top}\sim N(\mathbf{0},I_{s})$. The standard construction for generating Brownian motion is $(B_{T}^{(1)},\ldots,B_{T}^{(s)})^{\top}=A(Z_{1},\dots,Z_{s})^{\top}$, where $A$ is the Cholesky matrix of $C$. By contrast, the PCA chooses $A=[\sqrt{\lambda_{1}}\mathbf{v}_{1},\cdots,\sqrt{\lambda_{s}}\mathbf{v}_{s}]$, where $\lambda_{1}\geq\cdots\geq\lambda_{s}$ are the eigenvalues and $\mathbf{v}_{1},\ldots,\mathbf{v}_{s}$ are the corresponding unit-length eigenvectors of $C$. Expectation (13) is expressed as