Quasi-Monte Carlo for discontinuous integrands with singularities along the boundary of the unit cube

TL;DR

This paper improves the understanding of randomized quasi-Monte Carlo methods' accuracy when applied to discontinuous integrands with boundary singularities, especially relevant in financial engineering applications.

Contribution

It establishes new error bounds for randomized QMC on discontinuous functions with boundary singularities, including cases with boundary alignment to coordinate axes.

Findings

Expected error rate of O(n^{-1/2-1/(4d-2)+ε}) under mild conditions

Better error rates if discontinuity boundaries are aligned with axes

Error rate of O(n^{-1+ε}) for QMC-friendly discontinuities

Abstract

This paper studies randomized quasi-Monte Carlo (QMC) sampling for discontinuous integrands having singularities along the boundary of the unit cube $[0,1]^d$. Both discontinuities and singularities are extremely common in the pricing and hedging of financial derivatives and have a tremendous impact on the accuracy of QMC. It was previously known that the root mean square error of randomized QMC is only $o(n^{-1/2})$ for discontinuous functions with singularities. We find that under some mild conditions, randomized QMC yields an expected error of $O(n^{-1/2-1/(4d-2)+\epsilon})$ for arbitrarily small $\epsilon>0$. Moreover, one can get a better rate if the boundary of discontinuities is parallel to some coordinate axes. As a by-product, we find that the expected error rate attains $O(n^{-1+\epsilon})$ if the discontinuities are QMC-friendly, in the sense that all the discontinuity boundaries are parallel to coordinate axes. The results can be used to assess the QMC accuracy for some typical problems from financial engineering.

Equations

$$f(\mathbf{u})=g(\mathbf{u})\mathbb{I}\{\mathbf{u}\in \Omega \},$$

$$I(f)={\int }_{[0,1{]}^d}f(\mathbf{u})\mathrm{d}\mathbf{u}$$

$$\hat{I}(f)=\frac{1}{n}\sum _{i=1}^{n}f({\mathbf{u}}_i),$$

$$\tilde{f}(\mathbf{u})=\tilde{g}(\mathbf{u})\mathbb{I}\{\mathbf{u}\in \Omega \}.$$

$$\left|I(f)-\hat{I}(f)\right|\le \left|I(f)-I(\tilde{f})\right|+\left|I(\tilde{f})-\hat{I}(\tilde{f})\right|+\left|\hat{I}(\tilde{f})-\hat{I}(f)\right|.$$

$$\mathbb{E}\left[\left|\hat{I}(\tilde{f})-\hat{I}(f)\right|\right]\le \frac{1}{n}\sum _{i=1}^n\mathbb{E}\left[\left|f({\mathbf{u}}_i)-\tilde{f}({\mathbf{u}}_i)\right|\right]=I(|f-\tilde{f}|).$$

$$\mathbb{E}\left[|I(f)-\hat{I}(f)|\right]\le 2I(|f-\tilde{f}|)+\mathbb{E}\left[|I(\tilde{f})-\hat{I}(\tilde{f})|\right].$$

$$E=\prod_{i=1}^{d}\bigg{[}\frac{t_{i}}{b^{k_{i}}},\frac{t_{i}+1}{b^{k_{i}}}\bigg{)},$$

$${\tilde{a}}_{ijk}={\pi }_{j{a}_{ij1}{a}_{ij2}\dots a_{ijk-1}}(a_{ijk}).$$

$$\mathcal{M}(\partial \Omega )=\underset{ϵ\downarrow 0}{\lim }\frac{{\lambda }_d((\partial \Omega {)}_{ϵ})}{2ϵ},$$

$$\Omega ={\Omega }_u\times \prod _{i\notin u}[a_i,b_i),$$

$$\mathrm{Var}[\hat{I}(f)]\le c_{d,\Omega }{M}_g^2{n}^{-1-1/(2d-1)}(\log n{)}^{2d/(2d-1)},$$

$$M_g=\max \left(V_{\mathrm{HK}}(g),\underset{\mathbf{u}\in [0,1{]}^d}{\sup }\left|g(\mathbf{u})\right|\right).$$

$$\mathrm{Var}[\hat{I}(f)]\le c_{d,\Omega }{M}_g^2{n}^{-1-1/(2d_u-1)}(\log n{)}^{2d/(2d_u-1)}.$$

$$\left|{\partial }^{v}g(\mathbf{u})\right|\le B\prod _{i\in v}\min (u_i,1-u_i{)}^{-A_i-1}\prod _{i\notin v}\min (u_i,1-u_i{)}^{-A_i}$$

$$K(ϵ)=\{\mathbf{u}\in [0,1{]}^d|\prod _{1\le i\le d}\min (u_i,1-u_i)\ge ϵ\},$$

$$g(\mathbf{u})=g(\mathbf{c})+\sum _{v\ne \varnothing }{\int }_{[{\mathbf{c}}^v,{\mathbf{u}}^v]}{\partial }^{v}g({\mathbf{z}}^v:{\mathbf{c}}^{-v})\mathrm{d}{\mathbf{z}}^v,$$

$$g_{ϵ}(\mathbf{u})=g(\mathbf{c})+\sum _{v\ne \varnothing }{\int }_{[{\mathbf{c}}^v,{\mathbf{u}}^v]}{\partial }^{v}g({\mathbf{z}}^v:{\mathbf{c}}^{-v})\mathbb{I}\{{\mathbf{z}}^v:{\mathbf{c}}^{-v}\in K(ϵ)\}\mathrm{d}{\mathbf{z}}^v,$$

$$V_{\mathrm{HK}}(g_{ϵ})\le C_1{ϵ}^{-{\max }_i{A}_i-\eta }.$$

$$I(\left|f-f_{ϵ}\right|)\le C_2{ϵ}^{1-{\max }_i{A}_i-\eta }.$$

$$\underset{\mathbf{u}\in [0,1{]}^d}{\sup }\left|g_{ϵ}(\mathbf{u})\right|\le C_3{ϵ}^{-{\max }_i{A}_i-\eta }.$$

$$\left|g_{ϵ}(\mathbf{u})\right|$$

$$\le \left|g(\mathbf{c})\right|+B\sum _{v\ne \varnothing }{I}_v\prod _{i\notin v}\min (c_i,1-c_i{)}^{-A_i},$$

$$I_v:={\int }_{[0^v,1^v]}\prod _{i\in v}\min (z_i,1-z_i{)}^{-A_i-1}\mathbb{I}\{{\mathbf{z}}^v:{\mathbf{c}}^{-v}\in K(ϵ)\}\mathrm{d}{\mathbf{z}}^v.$$

$$I_v$$

$$=2{\int }_{[0^{\tilde{v}},1^{\tilde{v}}]}\prod _{i\in \tilde{v}}\min (z_i,1-z_i{)}^{-A_i-1}\left({\int }_{ϵ/e({\mathbf{z}}_{\tilde{v}})}^{1/2}{y}^{-A_{m(v)}-1}\mathrm{d}y\right)\mathrm{d}{\mathbf{z}}_{\tilde{v}}$$

$$\le 2{\int }_{[0^{\tilde{v}},1^{\tilde{v}}]}\prod _{i\in \tilde{v}}\min (z_i,1-z_i{)}^{-A_i-1}\frac{(ϵ/e({\mathbf{z}}_{\tilde{v}}){)}^{-A_{m(v)}}}{A_{m(v)}}\mathrm{d}{\mathbf{z}}_{\tilde{v}}$$

$$=2\frac{{ϵ}^{-A_{m(v)}}}{A_{m(v)}}\prod _{i\notin v}\min (c_i,1-c_i{)}^{A_{m(v)}}{\int }_{[0^{\tilde{v}},1^{\tilde{v}}]}\prod _{i\in \tilde{v}}\min (z_i,1-z_i{)}^{A_{m(v)}-A_i-1}\mathrm{d}{\mathbf{z}}_{\tilde{v}}$$

$$\le 2\frac{{ϵ}^{-A_{m(v)}}}{A_{m(v)}}\prod _{i\notin v}\min (c_i,1-c_i{)}^{A_{m(v)}}\prod _{i\in \tilde{v}}\frac{2}{A_{m(v)}-A_i}$$

$$=C_v{ϵ}^{-A_{m(v)}},$$

$$\left|g_{ϵ}(\mathbf{u})\right|\le \left|g(\mathbf{c})\right|+\tilde{B}{ϵ}^{-{\max }_i{A}_i}$$

$$\mathbb{E}\left[\left|I(f)-\hat{I}(f)\right|\right]=O(n^{-\gamma (1/2+1/(4d-2))}(\log n{)}^{\gamma d/(2d-1)}),$$

$$\mathbb{E}\left[\left|I(f)-\hat{I}(f)\right|\right]=O(n^{-\gamma (1/2+1/(4d_u-2))}(\log n{)}^{\gamma d/(2d_u-1)}).$$

Full text

This work was supported by the National Science Foundation of China under grant $71601189$.

1. Introduction

It is known that quasi-Monte Carlo (QMC) integration over the unit cube $[0,1]^{d}$ yields an asymptotic error rate of $O(n^{-1}(\log n)^{d})$ when the integrand has bounded variation in the sense of Hardy and Krause (BVHK); see [9] for details. In this paper we consider integrands that are discontinuous and have singularities along the boundary of the unit cube $[0,1]^{d}$. Such integrands cannot be BVHK because they are unbounded. Both discontinuities and singularities are extremely common in computational finance. Specifically, many problems arising from option pricing can be formulated as an integral over an unbounded domain $\mathbb{R}^{d}$ (see Glasserman [5] and references therein). A necessary first step in applying QMC methods to a practical integral formulated over $\mathbb{R}^{d}$ is to transform the integral into the unit cube $[0,1]^{d}$. The transformation may introduce singularities at the boundary. In addition, discontinuities appear in the pricing and hedging of financial derivatives (e.g., barrier options) and have a tremendous impact on the accuracy of QMC method [6, 16].

Formally, we are interested in integrands of the form

[TABLE]

where $\Omega\subset[0,1]^{d}$ and $g$ has singularities along the unit cube $[0,1]^{d}$. The integrand $f$ has a singularity at the boundary if $\Omega\cap[0,1]^{d}\neq\emptyset$. The QMC estimate of the integral

[TABLE]

is given by the average of $n$ samples

[TABLE]

with carefully chosen $\bm{u}_{i}\in[0,1]^{d}$. He and Wang [7] studied the convergence rate of RQMC for discontinuous functions of the form (1.1), but $g$ is assumed to be BVHK that excludes singularities. Under some mild assumptions on $\Omega$, they proved that the root mean square error of RQMC is $O(n^{-1/2-1/(4d-2)+\epsilon})$ for any $\epsilon>0$. If some discontinuity boundaries are parallel to some coordinate axes, the rate can be further improved to $O(n^{-1/2-1/(4d_{u}-2)+\epsilon})$, where $d_{u}$ denotes the so-called irregular dimension, that is the number of axes which are not parallel to the discontinuity boundaries. The results in He and Wang [7] cannot be applied to our setting because $g$ is not BVHK.

Owen [12] considered functions singular around any or all of the $2^{d}$ corners of $[0,1]^{d}$ and got some error rates that can be as good as $O(n^{-1+\epsilon})$ if the singular function obeys a strict enough growth rate. Owen [13] found the convergence rate of RQMC for integrands with point singularities with unknown locations. More recently, Basu and Owen [2] considered functions on the square $[0,1]^{2}$ that may are singular along a diagonal in the square. A key strategy in [12], [13] and [2] is to employ another function that has finite variation to approximate the singular function. The approximation has low variation. Motivated by these works, we use a low variation approximation $\tilde{g}$ to $g$, resulting in an approximation of $f$, given by

[TABLE]

Then using triangle inequality gives

[TABLE]

Suppose that $\bm{u}_{1},\dots,\bm{u}_{n}$ in (1.2) are RQMC points where each $\bm{u}_{i}\sim\mathbb{U}([0,1]^{d})$ individually; then

[TABLE]

As a result,

[TABLE]

To get an expected error bound, it suffices to bound the approximation error $I{(|f-\tilde{f}|)}$ and the RQMC integration error for the function (1.3). An upper bound of the approximation error can be obtained similarly as in Owen [12], which requires a growth condition on $g$. For the later, we will follow the analysis in [7] since $\tilde{g}$ is BVHK.

This paper finds some rates of convergence for RQMC integration of the function (1.1). Suppose that $g$ obeys a strict enough growth rate. We find that the expected error in RQMC is $O(n^{-1/2-1/(4d-2)+\epsilon})$. Moreover, one can get a better rate $O(n^{-1/2-1/(4d_{u}-2)+\epsilon})$ if the boundary of $\Omega$ is parallel to some coordinate axes. These results are similar to those found in He and Wang [7]. As a by-product, the expected error rate attains $O(n^{-1+\epsilon})$ if the discontinuities involved in (1.1) are QMC-friendly (they are parallel to coordinate axes as discussed in [15]). Our theoretical results can explain why QMC integration can be effective for some problems with both discontinuities and singularities in financial engineering. They also reveal the effects of discontinuities and singularities on QMC accuracy.

This paper is organized as follows. Section 2 gives the background on $(t,m,d)$-nets, $(t,d)$-sequences and the randomization technique proposed by [10]. The singular function $g$ is assumed to satisfy the growth condition. Some results from [7] are also reviewed. Convergence results of the expected error in RQMC for the function (1.1) are formally stated and proved in Section 3. Section 4 presents some examples arising from computational finance in which the growth condition is satisfied with arbitrarily small positive rates. Section 5 concludes this paper.

2. Background

2.1. Digital nets and sequences

Throughout this paper, we work with scrambled nets and sequences following the framework of He and Wang [7]. The integer $b\geq 2$ serves as a base. To begin with, we define an elementary interval in base $b$.

Definition 2.1**.**

An elementary interval in base $b$ is a subset of $[0,1)^{d}$ of the form

[TABLE]

where $k_{i},t_{i}\in\mathbb{N}$ with $t_{i}<b^{k_{i}}$ for $i=1,\dots,d$.

Definition 2.2**.**

Let $t$ and $m$ be nonnegative integers with $t\leq m$. A finite sequence $\bm{u}_{1},\dots,\bm{u}_{b^{m}}\in[0,1)^{d}$ is a $(t,m,d)$-net in base $b$ if every elementary interval in base $b$ of volume $b^{t-m}$ contains exactly $b^{t}$ points of the sequence.

Definition 2.3**.**

Let $t$ be a nonnegative integer. An infinite sequence $(\bm{u}_{i})_{i\geq 1}$ with $\bm{u}_{i}\in[0,1)^{d}$ is a $(t,d)$-sequence in base $b$ if for all $k\geq 0$ and $m\geq t$ the finite sequence $\bm{u}_{kb^{m}+1},\dots,\bm{u}_{(k+1)b^{m}}$ is a $(t,m,d)$-net in base $b$.

2.2. Scrambling

Owen [10] applied a scrambling scheme on the nets that retains the net property. Let $\bm{u}_{1},\dots,\bm{u}_{n}$ be a $(t,m,d)$-net or the first $n$ elements of a $(t,d)$-sequence in base $b$ where $\bm{u}_{i}=(u_{i}^{1},\dots,u_{i}^{d})$. We may write the components of $\bm{u}_{i}$ in their base $b$ expansion $u_{i}^{j}=\sum_{k=1}^{\infty}a_{ijk}b^{-k},$ where $a_{ijk}\in\{0,\dots,b-1\}$ for all $i,j,k$. The scrambled version of $\bm{u}_{1},\dots,\bm{u}_{n}$ is a sequence $\tilde{\bm{u}}_{1},\dots,\tilde{\bm{u}}_{n}$ with $\tilde{\bm{u}}_{i}=(\tilde{u}_{i}^{1},\dots,\tilde{u}_{i}^{d})$ written as $\tilde{u}_{i}^{j}=\sum_{k=1}^{\infty}\tilde{a}_{ijk}b^{-k},$ where $\tilde{a}_{ijk}$ are defined in terms of random permutations of the $a_{ijk}$. The permutation applied to $a_{ijk}$ depends on the values of $a_{ijh}$ for $h=1,\dots,k-1$. Specifically, $\tilde{a}_{ij1}=\pi_{j}(a_{ij1}),\ \tilde{a}_{ij2}=\pi_{ja_{ij1}}(a_{ij2}),\ \tilde{a}_{ij3}=\pi_{ja_{ij1}a_{ij2}}(a_{ij3})$, and in general

[TABLE]

Each permutation $\pi_{\bullet}$ is uniformly distributed over the $b!$ permutations of $\{0,\dots,b-1\}$, and the permutations are mutually independent.

2.3. Convergence results from He and Wang [7]

He and Wang [7] considered integrands of the form $f(\bm{u})=g(\bm{u})\mathbb{I}\{\bm{u}\in\Omega\}$, where $g$ is BVHK and the boundary of $\Omega$ admits a $(d-1)$-dimensional Minkowski content defined below.

Definition 2.4**.**

For a set $\Omega\subset[0,1]^{d}$, define

[TABLE]

where $(A)_{\epsilon}:=\{x+y|x\in A,\|y\|\leq\epsilon\}$, and $\|\cdot\|$ denotes the usual Euclidean norm. If $\mathcal{M}(\partial\Omega)$ exists and finite, then $\partial\Omega$ is said to admit a $(d-1)$-dimensional Minkowski content.

In the terminology of geometry, $\mathcal{M}(\partial\Omega)$ is known as the surface area of the set $\Omega$. The Minkowski content has a clear intuitive basis, compared to the Hausdorff measure [8] that provides an alternative to quantify the surface area. We should note that the Minkowski content coincides with the Hausdorff measure, up to a constant factor, in regular cases. It is known that the boundary of any convex set in $[0,1]^{d}$ has a $(d-1)$-dimensional Minkowski content. In this case, $\mathcal{M}(\partial\Omega)\leq 2d$ since the surface area of a convex set in $[0,1]^{d}$ is bounded by the surface area of the unit cube $[0,1]^{d}$, which is $2d$. More generally, Ambrosio et al. [1] found that if $\Omega$ has a Lipschitz boundary, then $\partial\Omega$ admits a ($d-1$)-dimensional Minkowski content. In their terminology, a set $\Omega$ is said to have a Lipschitz boundary if for every boundary point $a$ there exists a neighborhood $A$ of $a$, a rotation $R$ in $\mathbb{R}^{d}$ and a Lipschitz function $f:\mathbb{R}^{d-1}\to\mathbb{R}$ such that $R(\Omega\cap A)=\left\{(x,y)\in(\mathbb{R}^{d-1}\times\mathbb{R})\cap R(A)|y\geq f(x)\right\}$. In other words, $\Omega\cap A$ is the epigraph of a Lipschitz function.

He and Wang [7] showed that a faster convergence rate of RQMC can be achieved if the set $\Omega$ has some regularity. They studied partially axis-parallel sets as defined below. For a positive integer $d$, denote $1{:}d=\{1,\dots,d\}$. For a set $u\subset{1{:}d}$, denote the cardinality of $u$ as $\left|u\right|$ and $-u=1{:}d\backslash u$.

Definition 2.5**.**

A set $\Omega$ is said to be a partially axis-parallel set with irregular dimension $d_{u}=\left|u\right|$ if

[TABLE]

where $u\subset 1{:}d$, $d_{u}<d$, $0\leq a_{i}<b_{i}\leq 1$ for $i\notin u$, and $\Omega_{u}$ is a Lebesgue measurable subset of $\prod_{i\in u}[0,1)$. The quantity $d_{u}$ counts the number of axes which are not parallel to the boundaries of $\Omega$.

Denote $V_{\mathrm{HK}}(g)$ as the variation of the function $g$ in the sense of Hardy and Krause. See [11] for an outline of the variation. The following proposition summarizes the convergence results found in He and Wang [7].

Proposition 2.6**.**

Suppose that $f(\bm{u})=g(\bm{u})\mathbb{I}\{\bm{u}\in\Omega\}$, where $g\in L^{2}[0,1]^{d}$ satisfies $V_{\mathrm{HK}}(g)<\infty$. Assume that the sequence $\bm{u}_{1},\dots,\bm{u}_{n}$ in (1.2) is a scrambled $(t,d)$-sequence in base $b\geq 2$. If $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content, then for all sufficiently large $n$,

[TABLE]

where $c_{d,\Omega}$ depends only on $\Omega$ and $d$, and

[TABLE]

If $\Omega$ is a partially axis-parallel set with irregular dimension $d_{u}$ defined by (2.3), where $\partial\Omega_{u}$ admits a $(d_{u}-1)$-dimensional Minkowski content, then for all sufficiently large $n$,

[TABLE]

Proof.

The first part (2.4) has been proven in Theorem 3.5 of [7], and the second part (2.6) has been proven in Theorem 3.6 of [7]. Here we show the implicit constants in the upper bounds of the scrambled net variances because they are useful in the following proofs. ∎

Functions of BVHK must necessarily be bounded. So $M_{g}$ given by (2.5) is finite. However, for many practical problems arising from computational finance, $g$ has singularities on the surface of the unit cube $[0,1]^{d}$. In this cases, $g$ is unbounded, and hence $g$ has infinite variation. The conditions in Proposition 2.6 are thus not satisfied. Before establishing the convergence rate of scrambled net errors for singular integrands, we suppose that $g$ satisfies the growth condition as studied in Owen [12].

2.4. Growth Condition

For a set $v\subseteq 1{:}d$, $\partial^{v}g$ denotes the mixed partial derivative of $g$ taken once with respect to components with indices in $v$. Following Owen [12], we first introduce a growth condition for $g$ on $(0,1)^{d}$ that may become singular at the boundary of $[0,1]^{d}$ as shown in some integrands in the valuation of options with unbounded payoffs (see Section 4 for some examples).

Definition 2.7**.**

A function $g$ defined on $(0,1)^{d}$ is said to satisfy the boundary growth condition if

[TABLE]

holds for some $A_{i}>0$, some $B<\infty$ and all $v\subseteq 1{:}d$.

The boundary growth condition is the second growth condition described in Owen [12]. Owen [13] and Basu and Owen [2] studied other types of growth conditions for point singularities and singularities along a diagonal in the square, respectively. Large values of $A_{i}$ correspond to more severe singularities. When $\max_{i}A_{i}\geq 1$ the upper bound for $\left|g\right|$ is not even integrable. When $\max_{i}A_{i}<1/2$, then $f^{2}$ is integrable and Monte Carlo sampling has a root mean square error of $O(n^{-1/2})$. We use a region to avoid the singularities as

[TABLE]

for small $\epsilon>0$. We now define an extension $g_{\epsilon}$ of $g$ from $K(\epsilon)$ to $[0,1]^{d}$ such that $g_{\epsilon}(\bm{u})=g(\bm{u})$ for $\bm{u}\in K(\epsilon)$.

Definition 2.8**.**

A set $K\subset[0,1]^{d}$ is said to be Sobol’ extensible with anchor $\bm{c}$ if for every $\bm{u}\in K$ the rectangle $\prod_{i=1}^{d}[\min(u_{i},c_{i}),\max(u_{i},c_{i})]\subset K$.

It is easy to see that $K(\epsilon)$ is Sobol’ extensible with anchor $\bm{c}=(1/2,\dots,1/2)$. So one may write

[TABLE]

and then the desired low variation approximation of $g$ is given by

[TABLE]

where $\bm{z}^{v}{:}\bm{c}^{-v}$ denotes the point $\bm{y}\in[0,1]^{d}$ with $y_{j}=z_{j}$ for $j\in v$ and $y_{j}=c_{j}$ for $j\notin v$.

3. Expected errors in RQMC

Proposition 3.1**.**

If $g$ satisfies the boundary growth condition (2.7), then for any $\eta>0$ there exists $C_{1}<\infty$ such that

[TABLE]

If there is a unique maximum among $A_{1},\dots,A_{d}$, then (3.1) holds with $\eta=0$.

Proof.

See the proof of Theorem 5.5 in [12]. ∎

Proposition 3.2**.**

Let $f_{\epsilon}(\bm{u})=g_{\epsilon}(\bm{u})\mathbb{I}\{\bm{u}\in\Omega\}$, where $g_{\epsilon}$ is given by (2.10). If $g$ satisfies the boundary growth condition (2.7) with $\max_{i}A_{i}<1$, then for any $\eta\in(0,1-\max_{i}A_{i})$, there exists $C_{2}<\infty$ such that

[TABLE]

If there is a unique maximum among $A_{1},\dots,A_{d}$, then (3.2) holds with $\eta=0$.

Proof.

From the proof of Theorem 5.5 in [12] which is based on a result in [14], we have $I(\left|g-g_{\epsilon}\right|)\leq C_{2}\epsilon^{1-\max_{i}A_{i}-\eta}$. The upper bound (3.2) then follows from $I(\left|f-f_{\epsilon}\right|)=I(\left|g-g_{\epsilon}\right|\mathbb{I}\{\bm{u}\in\Omega\})\leq I(\left|g-g_{\epsilon}\right|)$. ∎

Proposition 3.3**.**

If $g$ satisfies the boundary growth condition (2.7), then for any $\eta>0$ there exists $C_{3}<\infty$ such that

[TABLE]

If there is a unique maximum among $A_{1},\dots,A_{d}$, then (3.3) holds with $\eta=0$.

Proof.

The procedure is similar to the proof of Theorem 5.5 in [12]. Combining (2.10) with the boundary growth condition (2.7), we have

[TABLE]

where

[TABLE]

We first assume that $A_{1},\dots,A_{d}$ are distinct positive numbers. Let $m(v)=\arg\max_{i\in v}A_{i}$ and $\tilde{v}=v-\{m(v)\}$. Let $e(\bm{z}_{\tilde{v}})=\prod_{i\in\tilde{v}}\min(z_{i},1-z_{i})\prod_{i\notin v}\min(c_{i},1-c_{i})$. Then

[TABLE]

where $C_{v}$ is a finite constant. It then follows from (3.4) that

[TABLE]

for some finite $\tilde{B}$.

If $A_{j}=A_{k}<\max_{i}A_{i}$ for some $j\neq k$, then we increase some of the $A_{i}$ so that $A_{1},\dots,A_{d}$ are distinct, while leaving $\max_{i}A_{i}$ unchanged. Then (3.5) also holds if there is a unique maximum among $A_{1},\dots,A_{d}$. We thus have (3.3) with $\eta=0$ due to $\left|g(\bm{c})\right|<\infty$. If there are two or more maximums among $A_{1},\dots,A_{d}$, then these maximums can be increased to distinct values, while raising $\max_{i}A_{i}$ by no more than $\eta$. ∎

Theorem 3.4**.**

Suppose that $f$ is given by (1.1), where $g$ satisfies the boundary growth condition (2.7) with $\max_{i}A_{i}<1$. Assume that the sequence $\bm{u}_{1},\dots,\bm{u}_{n}$ in (1.2) is a scrambled $(t,d)$-sequence in base $b\geq 2$. If $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content, then for any $\eta\in(0,1-\max_{i}A_{i})$,

[TABLE]

where $\gamma=1-\max_{i}A_{i}-\eta$. If $\Omega$ is a partially axis-parallel set with irregular dimension $d_{u}$ defined by (2.3), where $\partial\Omega_{u}$ admits a $(d_{u}-1)$-dimensional Minkowski content, then

[TABLE]

If there is a unique maximum among $A_{1},\dots,A_{d}$, then (3.6) and (3.7) hold with $\gamma=1-\max_{i}A_{i}$.

Proof.

Using the triangle inequality and the unbiasedness of the estimate $\hat{I}(f_{\epsilon})$, we have

[TABLE]

Proposition 3.2 gives $I(\left|f-f_{\epsilon}\right|)=O(\epsilon^{\gamma})$, where $\gamma=1-\max_{i}A_{i}-\eta$. For the function $f_{\epsilon}(\bm{u})=g_{\epsilon}(\bm{u})\mathbb{I}\{\bm{u}\in\Omega\}$, it follows from Propositions 2.6, 3.1 and 3.3 that

[TABLE]

Consequently,

[TABLE]

Taking $\epsilon\propto n^{-1/2-1/(4d-2)}(\log n)^{d/(2d-1)}$ establishes (3.6). The rate (3.7) can be proved in the same way. ∎

From Theorem 3.4, the rates for discontinuous integrands with singularities are not faster than those in Proposition 2.6. RQMC is asymptotically superior to Monte Carlo when $A_{i}<1/(2d)$ for all $i$. For some applications in computational finance (see Section 4 for some examples), it is possible that $g$ obeys the growth condition (2.7) with arbitrarily small positive $A_{i}$ for all $i$. The associated rates are presented in the following corollary, which are asymptotically superior to plain Monte Carlo sampling. In this case, the singularities may be regarded as QMC-friendly singularities because they deliver the best possible rate in our setting.

Corollary 3.5**.**

Suppose that $f$ is given by (1.1), where $g$ satisfies the boundary growth condition with arbitrarily small positive $A_{i}$ for all $i$. Assume that the sequence $\bm{u}_{1},\dots,\bm{u}_{n}$ in (1.2) is a scrambled $(t,d)$-sequence in base $b\geq 2$. If $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content, then

[TABLE]

for arbitrary small $\epsilon>0$. If $\Omega$ is a partially axis-parallel set with irregular dimension $d_{u}$ defined by (2.3), where $\partial\Omega_{u}$ admits a $(d_{u}-1)$-dimensional Minkowski content, then

[TABLE]

4. Examples from computational finance

Let $S(t)$ denote the underlying price dynamics at time $t$ under the risk-neutral measure. In a simulation framework, it is common that the prices are simulated at discrete times $t_{1},\dots,t_{d}$ satisfying $0=t_{0}<t_{1}<\dots<t_{d}=T$, where $T$ is the maturity of the financial derivative of interest. Without loss of generality, we assume that the discrete times are evenly spaced, i.e., $t_{i}=i\Delta t$, where $\Delta t=T/d$. For simplicity, denote $S_{i}=S(t_{i})$, and $\bm{S}=(S_{1},\dots,S_{d})^{\top}$. Under the risk-neutral measure, the price and the sensitivities of the financial derivative can be expressed as an expectation $I=\mathbb{E}\left[f(\bm{S})\right]$ for a real function $f$ over $\mathbb{R}^{d}$. To translate the problem into QMC setting, we suppose that $S_{i}$ can be expressed as a function of $\bm{u}\sim\mathbb{U}([0,1]^{d})$, denoted by $S_{i}(\bm{u})$, after some appropriate transformations. Let $\bm{S}(\bm{u})=(S_{1}(\bm{u}),\dots,S_{d}(\bm{u}))^{\top}$. We thus have

[TABLE]

After the transformations, the integrand $f\circ\bm{S}$ is often unbounded at the boundary of the unit cube.

Many functions in the pricing and hedging of financial derivatives involve indicator functions, which can be expressed in the form

[TABLE]

where $g$ and $h$ are usually smooth functions over $\mathbb{R}^{d}$ (see [6]). For pricing financial options, the factor $g$ determines the magnitude of the payoff and $h(\bm{S})>0$ gives the payout condition. For calculating Greeks by the pathwise method, the target function often involves an indicator function as in (4.1) even though the underlying payoff function is continuous.

We assume that under the risk-neutral measure the asset follows the geometric Brownian motion

[TABLE]

where $r$ is the risk-free interest rate, $\sigma$ is the volatility and $B(t)$ is the standard Brownian motion. Under this assumption, the solution of (4.2) is analytically available

[TABLE]

where $S_{0}$ is the initial price of the asset. Let $\bm{x}:=(B(t_{1}),\dots,B(t_{d}))^{\,\mathrm{T}}$. We have $\bm{x}\sim N(\bm{0},\bm{\Sigma})$, where the entries of $\bm{\Sigma}$ are given by $\Sigma_{ij}=\Delta t\min(i,j)$.

Note that $\bm{\Sigma}$ is positive definite. Let $\bm{A}$ be a generating matrix satisfying $\bm{A}\bm{A}^{\top}=\bm{\Sigma}$. Let $\Phi$ be the cumulative distribution function of the standard normal distribution. Using the transformation $\bm{x}=\bm{A}\Phi^{-1}(\bm{u})$, it follows from (4.3) that

[TABLE]

To verify the boundary growth condition, we need partial derivatives of ${g\circ\bm{S}}$ of order up to the dimension of the unit cube. The multivariate Faa di Bruno formula from [4] gives an arbitrary mixed partial derivative of $g\circ\bm{S}$ in terms of partial derivatives of $g$ and $S_{i}$. Basu and Owen [3] also used the formula to study the variation of some composition functions. The formula requires that the needed derivatives exist. Let $\bm{\lambda}=(\lambda_{1},\dots,\lambda_{d})$ be a vector of nonnegative integers. Define $\left|\bm{\lambda}\right|=\sum_{i=1}^{d}\lambda_{i}$. Denote $g_{\bm{\lambda}}$ as the derivative of $g$ taken $\lambda_{i}$ times with respect to the $i$th component. It follows from Basu and Owen [3] that for $\emptyset\neq v\subseteq 1{:}d$,

[TABLE]

where

[TABLE]

The following lemma is a result of Owen [12]. We prove it here also.

Lemma 4.1**.**

Suppose that $S_{i}$ is given by (4.4); then for any $v\subseteq 1{:}d$ and $i\in 1{:}d$,

[TABLE]

holds for arbitrarily small $A_{j}>0$ and $C_{i}<\infty$.

Proof.

It follows from (4.4) that

[TABLE]

Note that $\Phi^{-1}(\epsilon)=-\sqrt{-2\log(\epsilon)}+o(1)$ and $\Phi^{-1}(1-\epsilon)=\sqrt{-2\log(\epsilon)}+o(1)$ as $\epsilon\downarrow 0$. This leads to $\exp(a\Phi^{-1}(u_{j}))=O(\min(u_{j},1-u_{j})^{-A_{j}/2})$ for any $A_{j}>0$ and an arbitrary $a\in\mathbb{R}$. Denote $\phi(x)=(2\pi)^{-1/2}\exp(-x^{2}/2)$ as the probability density of the standard normal distribution. We find that

[TABLE]

for any $A_{j}>0$. The inequality (4.6) is thus obtained. ∎

Since the function $S_{i}$ admits the boundary growth condition for arbitrarily small $A_{j}>0$, Owen [12] showed that $S_{i}$ can be integrated with error $O(n^{-1+\epsilon})$ by the Halton sequence. However the results of Owen [12] cannot be applied to our target function (4.1) because it is discontinuous. To apply Theorem 3.4, we need to verify the boundary growth condition for the composition $g\circ\bm{S}$. Combining (4.5) and (4.6), we have

[TABLE]

for some finite $B_{1}$, arbitrarily small $A_{j}>0$ and $\emptyset\neq v\subseteq 1{:}d$. Therefore, the function $g\circ\bm{S}$ satisfies the growth condition (2.7) as long as

[TABLE]

holds for all $\left|\bm{\lambda}\right|\leq\left|v\right|$, $\tilde{A}_{j}>0$ and $B_{2}<\infty$. This may be verified for a broad range of functions since (4.6) admits that

[TABLE]

holds for arbitrarily small $A_{j}>0$. In our applications, $g$ is rather simple so that $g_{\bm{\lambda}}$ is available. As illustrative examples, we next show that the growth condition (2.7) can be satisfied with arbitrarily small growth rates.

Example 1**.**

The discounted payoff of an arithmetic Asian option is

[TABLE]

where $S_{A}=(1/d)\sum_{j=1}^{d}S_{j}$ and $K$ is the strike price.

Example 2**.**

The pathwise estimate of the delta of an arithmetic Asian option is

[TABLE]

The delta of an option is the sensitivity with respect to the initial price $S_{0}$ of the underlying asset.

Example 3**.**

An estimate of the gamma of an arithmetic Asian option is

[TABLE]

which results from applying the pathwise method first and then the likelihood ration method (see [5]). The gamma is the second derivative with respect to the initial price $S_{0}$ of the underlying asset.

Example 4**.**

The pathwise estimate of the rho of an arithmetic Asian option is

[TABLE]

where

[TABLE]

The rho of an option is the sensitivity with respect to the risk-free interest rate $r$.

Example 5**.**

The pathwise estimate of the theta of an arithmetic Asian option is

[TABLE]

where

[TABLE]

The theta of an option is the sensitivity with respect to the maturity of the option $T$.

Example 6**.**

The pathwise estimate of the vega of an arithmetic Asian option is

[TABLE]

in which

[TABLE]

The vega of an option is the sensitivity with respect to the volatility $\sigma$.

Theorem 4.2**.**

Suppose that $f$ is one of the functions (4.10)–(4.15), where $S_{i}$ is given by (4.4). Letting $f$ be expressed as the form (4.1), then $g\circ\bm{S}$ satisfies the boundary growth condition (2.7) with arbitrarily small $A_{i}>0$ for all $i$.

Proof.

For the functions (4.10)–(4.15), $g(\bm{S})$ is a linear combination of some components $S_{i}$ and $\log(S_{i})S_{i^{\prime}}$ for $i,i^{\prime}\in 1{:}d$. It suffices to verify that these components satisfy (4.8) because the linear combination then also satisfies (4.8).

Consider $g(\bm{S})=S_{i}$ for any $i\in 1{:}d$. We have $\left|g_{\bm{\lambda}}(\bm{S})\right|\leq 1$ for any $1\leq\left|\bm{\lambda}\right|\leq\left|v\right|$. For $\left|\bm{\lambda}\right|=0$, $\left|g_{\bm{\lambda}}(\bm{S})\right|=\left|g(\bm{S})\right|=S_{i}=O(\prod_{j=1}^{d}\min(u_{j},1-u_{j})^{-\tilde{A}_{j}})$ due to (4.6), for arbitrarily small $\tilde{A}_{j}>0$. In this case, $g_{\bm{\lambda}}$ satisfies (4.8) with arbitrarily small growth rates.

Consider $g(\bm{S})=\log(S_{i})S_{i^{\prime}}$ for any $i\neq i^{\prime}$. We have $g_{\bm{\lambda}}(\bm{S})=0$ if $\lambda_{k}>0$ for some $k\notin\{i,i^{\prime}\}$ or $\lambda_{i^{\prime}}>1$. So it reduces to consider three cases:

1. (1)

$1\leq\lambda_{i}\leq\left|v\right|,\ \lambda_{k}=0$ for any $k\neq i$; 2. (2)

$0\leq\lambda_{i}\leq\left|v\right|-1,\ \lambda_{i^{\prime}}=1,\ \lambda_{k}=0$ for $k\notin\{i,i^{\prime}\}$; and 3. (3)

all $\lambda_{k}=0$.

For Case (1), we have

[TABLE]

where $c(1)=1$ and $c(\lambda_{i})=(\lambda_{i}-1)!$ for $\lambda_{i}>1$. For Case (2), we have

[TABLE]

For Case (3),

[TABLE]

From the proof of Lemma 4.1, we have

[TABLE]

and

[TABLE]

for arbitrarily small $\tilde{A}_{j}>0$. So $g_{\bm{\lambda}}$ satisfies (4.8) with arbitrarily small growth rates.

Consider $g(\bm{S})=\log(S_{i})S_{i}$ for any $i\in 1{:}d$. If $\lambda_{k}=0$ for all $k\neq i$, we have

[TABLE]

If $\lambda_{k}>0$ for some $k\neq i$, then $g_{\bm{\lambda}}(\bm{S})=0$. In this case, $g_{\bm{\lambda}}$ satisfies (4.8) with arbitrarily small growth rates.

Based on the reasoning above, it follows from (4.7) that for the functions (4.10)–(4.15), $g\circ\bm{S}$ satisfies the boundary growth condition with arbitrarily small growth rates. ∎

Note that the statement in Theorem 4.2 holds for any decomposition of $\bm{\Sigma}=\bm{A}\bm{A}^{\top}$. To handle discontinuities, Wang and Tan [16] proposed the orthogonal transformation (OT) method to make the discontinuities QMC friendly, in the sense that all the discontinuity boundaries are parallel to coordinate axes. The OT method delivers a special matrix $\bm{A}$ satisfying $\bm{A}\bm{A}^{\top}=\bm{\Sigma}$ to generate the path (4.4). To illustrate its effects, let us consider the function

[TABLE]

where $S_{G}=\prod_{i=1}^{d}S_{i}^{1/d}$ is the geometric average of the prices. For this function, applying the OT method arrives at

[TABLE]

for some constant $\kappa$ (see [16] for determining the matrix $\bm{A}$). As a result, discontinuities occur only on the axis-parallel hyperplane $u_{1}=\kappa$, which are QMC-friendly. The function (4.16) is then transformed to $g(\bm{S}(\bm{u}))\mathbb{I}\{\bm{u}\in\Omega\}$, where $\Omega=\{\bm{u}\in[0,1]^{d}|u_{1}>\kappa\}$. Note that the irregular dimension $d_{u}$ of the set $\Omega$ is one. Corollary 3.5 admits that the expected error rate of RQMC for the transformed function $g(\bm{S}(\bm{u}))\mathbb{I}\{\bm{u}\in\Omega\}$ is $O(n^{-1+\epsilon})$ if $g\circ\bm{S}$ has the same kind of singularities as examined in Theorem 4.2. This suggests that making the discontinuities of the function (4.16) QMC friendly by the OT method can improve the efficiency of QMC greatly. For the functions of the form $g(\bm{S}(\bm{u}))\mathbb{I}\{S_{A}>K\}$ in the examples above, Wang and Tan [16] suggested that using the obtained matrix $\bm{A}$ for the function (4.16) can still be effective since $S_{G}$ is a good substitute for $S_{A}$. The usefulness of this strategy was illustrated by several numerical examples in [16, 6].

5. Conclusion

We find that for discontinuous functions with singularities along the boundary of the unit cube $[0,1]^{d}$, RQMC has an expected error of $O(n^{-\gamma(1/2+1/(4d-2))+\epsilon})$ for $\gamma=1-\max_{i}A_{i}\in(0,1)$ depending on the growth rates $A_{i}$. The convergence rate $O(n^{-\gamma(1/2+1/(4d-2))+\epsilon})$ is a bit disappointing for large values of $A_{i}$. However, the error rate can be as good as $O(n^{-(1/2+1/(4d-2))+\epsilon})$ for some problems from computational finance in which the growth rates are arbitrarily small. In these cases, it seems that the singularities have insignificant impact on QMC accuracy, compared to the rate for discontinuous integrands (without singularities) found in He and Wang [7]. We also show theoretically the benefits of making discontinuities QMC-friendly, which have been shown empirically in various numerical examples of Wang and Sloan [15] and Wang and Tan [16].

For singular functions (even discontinuous) satisfying the growth condition with arbitrarily small growth rates, QMC can lead to improved accuracy. It would be interesting to know how generally the problems from financial engineering fit into this setting, beyond those under the Gaussian model discussed in Section 4.

Acknowledgments

The author thanks two anonymous referees for helpful suggestions on improving this paper. The author also thanks Professor Art B. Owen and Kinjal Basu for sharing their work [2] with him.