Convergence analysis of quasi-Monte Carlo sampling for quantile and expected shortfall

TL;DR

This paper analyzes the convergence properties of quasi-Monte Carlo methods for estimating quantiles and expected shortfalls, providing theoretical error bounds and demonstrating improved convergence rates over traditional Monte Carlo techniques.

Contribution

It establishes convergence proofs and error bounds for QMC-based quantile and expected shortfall estimates, advancing the theoretical understanding of QMC in risk measurement.

Findings

QMC quantile estimates converge under mild conditions.

Deterministic error bound of $O(N^{-1/d})$ for QMC quantiles.

MSE of randomized QMC for expected shortfall can be $o(N^{-1})$ or better under certain conditions.

Abstract

Quantiles and expected shortfalls are usually used to measure risks of stochastic systems, which are often estimated by Monte Carlo methods. This paper focuses on the use of quasi-Monte Carlo (QMC) method, whose convergence rate is asymptotically better than Monte Carlo in the numerical integration. We first prove the convergence of QMC-based quantile estimates under very mild conditions, and then establish a deterministic error bound of $O(N^{-1/d})$ for the quantile estimates, where $d$ is the dimension of the QMC point sets used in the simulation and $N$ is the sample size. Under certain conditions, we show that the mean squared error (MSE) of the randomized QMC estimate for expected shortfall is $o(N^{-1})$. Moreover, under stronger conditions the MSE can be improved to $O(N^{-1-1/(2d-1)+\epsilon})$ for arbitrarily small $\epsilon>0$.

Equations

$$v:=F^{-1}(p)=\inf \{x\in \mathbb{R}|F(x)\ge p\},$$

$$c:=v-\frac{1}{p}\mathbb{E}[(v-X{)}^{+}],$$

$${\hat{F}}_N(x)=\frac{1}{N}\sum _{i=1}^{N}1\{X_i\le x\},$$

$${\hat{v}}_N$$

$${\hat{c}}_N={\hat{v}}_N-\frac{1}{pN}\sum _{i=1}^N({\hat{v}}_N-X_i{)}^{+}.$$

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

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

$$\left|{\hat{I}}_N(f)-I(f)\right|\le V_{\mathrm{HK}}(f)D_N^{\ast }(\mathcal{P}),$$

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

$$X=\varphi (\mathbf{u}),\mathbf{u}\sim \mathbb{U}([0,1{)}^d),$$

$${\hat{F}}_N(x)=\frac{1}{N}\sum _{i=1}^{N}1\{\varphi ({\mathbf{u}}_i)\le x\},$$

$${\hat{c}}_N={\hat{v}}_N-\frac{1}{pN}\sum _{i=1}^N({\hat{v}}_N-\varphi ({\mathbf{u}}_i){)}^{+}.$$

$$p\le {\hat{F}}_N({\hat{v}}_N)\le p+1/N$$

$$F(v-\eta )<p<F(v+\eta ).$$

$$F(v+ϵ)=\underset{i\to \infty }{\lim }{\hat{F}}_{n_i}(v+ϵ)\le p.$$

$$F(v-ϵ)=\underset{i\to \infty }{\lim }{\hat{F}}_{n_i}(v-ϵ)\ge p.$$

$${\Omega }_x:=\{\mathbf{u}\in [0,1{)}^d|\varphi (\mathbf{u})\le x\}$$

$$(A{)}_{ϵ}:=\{\mathbf{x}\in {\mathbb{R}}^d|{\Vert \mathbf{x}-\mathbf{y}\Vert }_2\le ϵ\mathrm{for}\mathrm{some}\mathbf{y}\in A\},$$

$${\lambda }_d((\partial \Omega {)}_{ϵ})\le g(ϵ)\text{ for all }ϵ>0.$$

$$\underset{x\in B(v,\delta )}{\sup }{\lambda }_d((\partial {\Omega }_x{)}_{ϵ})\le g(ϵ)$$

$$\left|{\hat{v}}_N-v\right|\le \frac{2g(\sqrt{d}{b}^{1+t/d}{N}^{-1/d})}{f_X(v)}$$

$$\left|{\hat{F}}_N(x)-F(x)\right|$$

$$=\left|\sum _{k\in \mathcal{T}}\frac{\#(E_k\cap {\Omega }_x)}{N}-\sum _{k\in \mathcal{T}}{\lambda }_d(E_k\cap {\Omega }_x)\right|$$

$$\le \max \left\{\sum _{k\in \mathcal{T}}\frac{\#(E_k\cap {\Omega }_x)}{N},\sum _{k\in \mathcal{T}}{\lambda }_d(E_k\cap {\Omega }_x)\right\}$$

$$\le \frac{b^t\left|\mathcal{T}\right|}{N}.$$

$$\sqrt{d}{b}^{-\lfloor (m-t)/d\rfloor }<\sqrt{d}{b}^{1+t/d}{N}^{-1/d}=:r(N).$$

$$\left|\mathcal{T}\right|\le \frac{{\lambda }_d((\partial {\Omega }_x{)}_r)}{b^{t-m}}\le b^{-t}Ng(r(N)).$$

$$\underset{x\in B(v,\delta )}{\sup }\left|{\hat{F}}_N(x)-F(x)\right|\le g(r(N)).$$

$${\hat{F}}_N(v+ϵ)\ge F(v+ϵ)-g(r(N)).$$

$$F(v+ϵ)-F(v)=f_X(\xi )ϵ$$

$${\hat{F}}_N(v+ϵ)\ge F(v)+f_X(\xi )ϵ-g(r(N))\ge \alpha +f_X(v)ϵ/2-g(r(N))=\alpha .$$

$$\left|{\hat{v}}_N-v\right|\le \frac{10d^{3/2}{b}^{1+t/d}}{f_X(v)}{N}^{-1/d}.$$

Full text

This work was supported by the National Science Foundation of China under Grant No. 71601189 and the National Key R&D Program of China under Grant No. 2016QY02D0301.

1. Introduction

Many application areas use quantiles or expected shortfalls to measure risks of stochastic systems. For instance, in the financial industry, a quantile (known as value-at-risk) plays an important role for quantifying and managing portfolio risk. On the other hand, expected shortfall (known as conditional value-at-risk) may provide incentives for risk managers to take into account tail risks beyond quantile. We refer to [8] for a review on the two measures. This paper focuses on estimating quantiles and expected shortfalls via simulation-based methods. Monte Carlo (MC) is a natural method to estimate them. However, the MC approach is often criticized for time-consuming, since value-at-risk estimation is often relevant to rare events simulation. That usually calls for a large number of runs to get accurate estimation. To address this issue, various variance reduction techniques are employed to increase the accuracy of MC. Importance sampling (IS) is a promising variance reduction technique for value-at-risk estimation (see, e.g., [5, 6]).

Beyond the use of MC, Avramidis and Wilson [1] proposed correlation-induction techniques to improve quantile estimation based on Latin hypercube sampling (LHS). They showed that the correlation-inducted LHS estimator is asymptotically normal and unbiased with smaller variance than that of the crude MC. Subsequently, Jin et al. [9] modified the correlation-inducted LHS estimator of [1] and proposed a new quantile estimator based on an indirect means of realizing full stratification of [16] that reuses samples. They showed that the error probability for the stratified quantile estimator is zero for sufficiently large, but finite sample size $N$. Moreover, in some special cases, the convergence rate is $O(N^{-1})$, as opposed to the crude MC rate $O(N^{-1/2})$. However, the stratified quantile estimator requires sample sizes that grow exponentially with the dimension of the problem.

Quasi-Monte Carlo (QMC) methods are deterministic versions of the MC methods, and have an asymptotically faster convergence rate than MC as shown in the field of numerical integration. It is straightforward to use QMC methods for estimating quantiles and expected shortfalls. Particularly, Papageorgiou and Paskov [18] observed from empirical studies that QMC methods provide a highly efficient alternative to MC for quantile calculation. Jin and Zhang [10] aimed at smoothing QMC estimators via Fourier transformation so that the faster convergence rate of QMC methods can be reclaimed. To the best of our knowledge, the convergence and the rates of convergence for plain QMC in estimating quantile and expected shortfall are still unclear.

In this paper, we focus on the use of QMC and randomized QMC (RQMC) for estimating quantile and expected shortfall. We first prove the convergence of QMC-based quantile estimates, and establish some useful error bounds for assessing the error rate. We then provide an error bound for the expected shortfall estimate, and find that the efficiency of the expected shortfall estimate is strongly tied to the efficiency of (R)QMC quadrature for a specific discontinuous function and a specific function with kinks. Under mild conditions, we show that the mean squared error (MSE) of the RQMC-based expected shortfall estimate is $o(N^{-1})$, which is asymptotically better than plain MC and LHS. Moreover, under stronger conditions the MSE can be improved to $O(N^{-1-1/(2d-1)+\epsilon})$ for arbitrarily small $\epsilon>0$, where $d$ is the dimension of the problem.

The rest of this paper is organized as follows. In Section 2, we introduce some background on quantile estimation and some preliminary results on QMC methods. In Section 3, we study the convergence and the convergence rate of QMC-based quantile estimate. In Section 4, we study the MSE of QMC-based expected shortfall estimate. In Section 5, we perform a numerical study for stochastic network models for which our theoretical results can be applied. Section 6 concludes this paper.

2. Preliminaries

Let $X$ be a real-valued random variable of interest with a cumulative distribution function (CDF) $F(x)$. For instance, $X$ is the loss or profit of a portfolio over a given holding period. We are interested in the left tail of the distribution of $X$. For a fixed $p\in(0,1)$, the quantity

[TABLE]

is called the $p$’th quantile of $X$ (or the value-at-risk of $X$ in the context of risk management). The expected shortfall of $X$ is defined as

[TABLE]

where $x^{+}:=\max\{x,0\}$. The expected shortfall is also known as the tail conditional expectation or conditional value-at-risk. Assume that the variable $X$ can be simulated easily. Our goal is to estimate the quantile $v$ and the expected shortfall $c$ by means of simulation.

In the MC setting, the CDF $F(x)$ of $X$ can be estimated by the empirical CDF

[TABLE]

where $X_{i}$’s are independent and identically distributed random replications of $X$. The quantile $v$ is then estimated by

[TABLE]

Let $X_{(i)}$ be the $i$th-order statistic of $X_{1},\dots,X_{N}$. It is easy to see that $\hat{v}_{N}=X_{(\lceil pN\rceil)}$, where $\lceil x\rceil$ denotes the smallest integer no less than $x$. The corresponding estimate of the expected shortfall $c$ is given by

[TABLE]

Serfling [20] showed that $\hat{v}_{N}\to v$ with probability 1 (w.p.1) as $N\to\infty$ under very mild assumptions. If $X$ has a continuous density $f_{X}(\cdot)$ in a neighborhood of $v$ and $f_{X}(v)>0$, Serfling [20] further showed that $\hat{v}_{N}$ is asymptotically normally distributed. For the expected shortfall estimate (2.5), Trindade et al. [22] found that under certain conditions, $\hat{c}_{N}\to c$ w.p.1 as $N\to\infty$, and $\hat{c}_{N}$ is asymptotically normally distributed.

MC is often criticized for its slow convergence. QMC has the potential to improve the convergence rate. We now turn to the regime of QMC in estimating quantiles and expected shortfalls. To start with, let’s consider the problem of estimating an integral over the unit cube $[0,1)^{d}$

[TABLE]

QMC quadrature rule takes the average

[TABLE]

where $\bm{u}_{1},\dots,\bm{u}_{N}$ are carefully chosen points in $[0,1)^{d}$. The Koksma-Hlawka inequality gives a deterministic error bound for the quadrature rule (2.6)

[TABLE]

where $\mathcal{P}:=\{\bm{u}_{1},\dots,\bm{u}_{N}\}$, $V_{\mathrm{HK}}(f)$ is the variation of $f(\bm{u})$ in the sense of Hardy and Krause, and $D^{*}_{N}(\mathcal{P})$ is the star-discrepancy of points in $\mathcal{P}$; see [13] for details. There are many ways to construct point sets such that $D^{*}_{N}(\mathcal{P})=O(N^{-1}(\log N)^{d})$. As a result, the QMC error is $O(N^{-1}(\log N)^{d})$ for integrands with bounded variation in the sense of Hardy and Krause (BVHK). In this paper, we restrict our attention to $(t,d)$-sequences or $(t,m,d)$-nets in base $b\geq 2$ (see the definitions below).

Definition 2.1**.**

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

[TABLE]

where $k_{j}\in{\rm Nature}$, $t_{j}\in{\rm Nature}$ with $t_{j}<b^{k_{j}}$ for $j=1,\dots,d$.

The elementary interval (2.8) is a hyperrectangle of volume $b^{-\sum_{j=1}^{d}k_{j}}$. For given $k_{j}$, the unit cube $[0,1)^{d}$ is partitioned into $b^{\sum_{j=1}^{d}k_{j}}$ elementary intervals of the form (2.8).

Definition 2.2**.**

Let $t$ and $m$ be nonnegative integers with $t\leq m$. A point set of $b^{m}$ points $\bm{u}_{1},...,\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 point set.

Definition 2.3**.**

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

For QMC, it is important to obtain an estimate of the quadrature error $|\hat{I}_{N}(f)-I(f)|$. But both the variation and the star discrepancy in the upper bound (2.7) are very hard to compute, and the upper bound is restricted to functions of finite variation. Instead, one can randomize the points $\bm{u}_{1},\dots,\bm{u}_{N}$ and treat the random version of the quadrature $\hat{I}_{N}(f)$ in (2.6) as an RQMC quadrature rule. Usually, the randomized points are uniformly distributed over $[0,1)^{d}$, and the low discrepancy property of the points is preserved under the randomization (see [11] and Chapter 13 of the monograph [3] for a survey of various RQMC methods). In this paper, we focus on the scrambling technique proposed by [14] to randomize $(t,d)$-sequences or $(t,m,d)$-nets.

QMC methods are designed to sample $d$-dimensional vectors that are uniformly distributed on the unit cube $[0,1)^{d}$. To fit into the setting of quantile estimation, one needs to know the mechanism of sampling the target variable $X$ via standard uniform distributed variables. In what follows, we assume that the target variable $X$ can be generated by

[TABLE]

where the function $\phi:[0,1)^{d}\to\mathbb{R}$ is easily computed. For practical problems, it may be easy to obtain the mapping $\phi$ by using the multivariate inverse transformation proposed by Rosenblatt [19]. In the QMC setting, we shall rewrite the empirical CDF (2.3) as

[TABLE]

where $\bm{u}_{1},\dots,\bm{u}_{N}$ are QMC or RQMC points. The QMC estimate of the quantile $v$ is then obtained by the formula (2.4). The expected shortfall estimate (2.5) is then replaced by

[TABLE]

Note that the empirical CDF (2.10) can be viewed as a QMC quadrature rule $\hat{I}_{N}(f)$ for the indicator function $f(\bm{u})=\bm{1}\{\phi(\bm{u})\leq x\}$. However, the upper bound (2.7) of the error does not provide useful information for this case because the variation $V_{\mathrm{HK}}(f)$ is usually infinite for discontinuous integrands [17].

3. Convergence analysis for QMC quantile estimation

In this section, we first show the convergence of QMC estimates for quantile estimation under very mild conditions. Then we give a deterministic error bound for QMC estimates under some relatively stronger conditions.

Theorem 3.1**.**

Let $X$ be a random variable with CDF $F(x)$. The empirical CDF $\hat{F}_{N}(x)$ given by (2.10) is based on QMC points. Assume that

1. (i)

$v=F^{-1}(p)$* is the unique solution $x$ of $F(x-)\leq p\leq F(x)$, and* 2. (ii)

$\lim_{N\to\infty}\hat{F}_{N}(x)=F(x)$* for all $x\in\mathbb{R}$.*

Then $\hat{v}_{N}=\hat{F}_{N}^{-1}(p)\to v$ as $N\to\infty$.

Proof.

By the definition of $\hat{v}_{N}$ in (2.4), we have

[TABLE]

for all $N\geq 1$. By the uniqueness condition (i) and the definition of $v$ in (2.1), we find that for any $\eta>0$,

[TABLE]

Assume that $\hat{v}_{N}$ does not converge to $v$. Then there exists an $\epsilon>0$ and an infinite sequence of positive integers $n_{i}$ with $\lim_{i\to\infty}n_{i}=\infty$ such that $|\hat{v}_{n_{i}}-v|\geq\epsilon$ for all $i$. If $\hat{v}_{n_{i}}\geq v+\epsilon$, then $\hat{F}_{n_{i}}(v+\epsilon)\leq\hat{F}_{n_{i}}(\hat{v}_{n_{i}})\leq p+1/n_{i}.$ By condition (ii),

[TABLE]

This leads to a contradiction because $F(v+\epsilon)>p$ by using (3.1).

On the other hand, if $\hat{v}_{n_{i}}\leq v-\epsilon$, then $\hat{F}_{n_{i}}(v-\epsilon)\geq\hat{F}_{n_{i}}(\hat{v}_{n_{i}})\geq p$. By condition (2),

[TABLE]

That also leads to a contradiction because $F(v-\epsilon)<p$. As a result, $\hat{v}_{N}$ converges to $v$ as $N$ goes to infinity. ∎

The uniqueness condition (i) is also the minimal requirement for establishing the strong consistency of the associated MC estimate [20, p. 75]. Condition (ii) implies that the empirical CDF $\hat{F}_{N}(x)$ converges to the true CDF $F(x)$ for all $x\in\mathbb{R}$ in the QMC setting. Recall that the random variable $X$ can be generated via the mapping (2.9). This calls for the Jordan measurability of the set

[TABLE]

for all $x\in\mathbb{R}$.

Corollary 3.2**.**

Suppose that the point set $\{\bm{u}_{1},\dots,\bm{u}_{N}\}$ used in (2.10) is the first $N$ points of a $(t,d)$-net in base $b\geq 2$. If $v=F^{-1}(p)$ is the unique solution $x$ of $F(x-)\leq p\leq F(x)$ and $\Omega_{x}$ is Jordan measurable for all $x\in\mathbb{R}$, then $\hat{v}_{N}\to v$ as $N\to\infty$.

Proof.

Let $f(\bm{u})=\bm{1}\{\phi(\bm{u})\leq x\}=\bm{1}\{\bm{u}\in\Omega_{x}\}$. Since $\Omega_{x}$ is Jordan measurable, $f(\bm{u})$ is Riemann integrable. Note that $\hat{I}_{N}(f)=\hat{F}_{N}(x)$ and $I(f)=F(x)$. It is known that the QMC quadrature $\hat{I}_{N}(f)$ converges to $I(f)$ for all Riemann-integrable functions $f$ (see, e.g., [13]). This implies that $\lim_{N\to\infty}\hat{F}_{N}(x)=F(x)$ for all $x\in\mathbb{R}$. Applying Theorem 3.1 completes the proof. ∎

Note that Riemann integrability of $\phi(\bm{u})$ may not lead to Riemann integrability of the indicator function $1\{\phi(\bm{u})\leq x\}$. Chen et al. [2] gave such an example by using Thomae’s function. By Lebesgue’s theorem (see [12]), $\bm{1}\{\phi(\bm{u})\leq x\}$ is Riemann integrable (or equivalently, $\Omega_{x}$ is Jordan measurable) iff $\lambda_{d}(\partial\Omega_{x})=0$, where $\lambda_{d}(\cdot)$ is the Lebesgue measure in $\mathbb{R}^{d}$.

Definition 3.3**.**

For a set $A\subset\mathbb{R}^{d}$, the outer parallel body of $A$ at distance $\epsilon$ is defined as

[TABLE]

where $\left\|\cdot\right\|_{2}$ denotes the Euclidean norm. When $A=\varnothing$, we use a convention that $(A)_{\epsilon}=\varnothing$ for any $\epsilon>0$.

Let $g=g(\epsilon)$ a positive nondecreasing function defined for all $\epsilon>0$ and satisfying $\lim_{\epsilon\to 0_{+}}g(\epsilon)=0$. Then we let $\mathcal{M}_{g}$ be the family of all Lebesgue-measurable $\Omega\subset[0,1]^{d}$ for which

[TABLE]

Every $\Omega\in\mathcal{M}_{g}$ is actually Jordan measurable. Conversely, every Jordan measurable subset of $[0,1)^{d}$ belongs to $\mathcal{M}_{g}$ for a suitable function $g$ (see [13, pp. 168-169]). To establish an error bound of quantile estimate, we need a stronger condition that $\lambda_{d}((\partial\Omega_{x})_{\epsilon})$ has a common upper bound $g(\epsilon)$ for $x$ in a neighborhood of $v$. For $\delta>0$, denote $B(x,\delta):=\{t\in\mathbb{R}||t-x|\leq\delta\}$ as a $\delta$-neighborhood of $x$.

Assumption 3.4**.**

Assume that $X$ has a density $f_{X}(x)$ in a neighborhood of $v$ and $f_{X}(x)$ is positive and continuous at $v$.

Assumption 3.5**.**

Assume that there exist a positive nondecreasing function $g(\epsilon)$ satisfying $\lim_{\epsilon\to 0_{+}}g(\epsilon)=0$ and $\delta,\epsilon_{0}>0$ such that

[TABLE]

for any $\epsilon\leq\epsilon_{0}$.

Theorem 3.6**.**

Suppose that the point set $\{\bm{u}_{1},\dots,\bm{u}_{N}\}$ used in (2.10) is a $(t,m,d)$-net in base $b\geq 2$, where $N=b^{m}$. If Assumptions 3.4 and 3.5 are satisfied, then

[TABLE]

for $N$ large enough.

Proof.

By Definition 2.2, there exist $K=b^{m-t}$ disjoint elementary intervals $E_{1},\dots,E_{K}$ with volume $b^{t-m}$ such that all $E_{k}$ contain exactly $b^{t}$ points of the $(t,m,d)$-net in base $b$. Let $\mathcal{T}:=\{k=1,\dots,K|E_{k}\cap\partial\Omega_{x}\neq\varnothing\}$, and denote $\#(A)$ as the number of the points of the $(t,m,d)$-net contained in the set $A$. By the fairness of the elementary intervals, we have

[TABLE]

Similarly to the proof of Lemma 4.1 in [7], one can choose $E_{k}$ with length as small as possible. By doing so, the length of $E_{k}$ is no larger than

[TABLE]

Let’s assume that $r(N)<\epsilon_{0}$ by taking large enough $N$. By Assumption 3.5,

[TABLE]

Therefore, we have

[TABLE]

By Assumption 3.4, there exists $\delta^{\prime}>0$ such that $f_{X}(x)\geq f_{X}(v)/2$ for all $x\in B(v,\delta^{\prime})$. Let $\epsilon=2g(r(N))/f_{X}(v)$. Since $r(N)\to 0$ as $N\to\infty$, there exists a $N_{0}(\delta,\delta^{\prime},\epsilon_{0})$ such that $r(N)\leq\epsilon_{0}$ and $\epsilon<\min(\delta,\delta^{\prime})$ for any $N\geq N_{0}$. By (3.2), we have

[TABLE]

By the mean value theorem, we obtain

[TABLE]

for some $\xi\in[v,v+\epsilon]\subset B(v,\delta^{\prime})$. Since $f_{X}(\xi)\geq f_{X}(v)/2$ and $F(v)=\alpha$,

[TABLE]

We therefore have $\hat{v}_{N}\leq v+\epsilon$. Conversely, we can prove in a similar manner that $\hat{v}_{N}\geq v-\epsilon$. Consequently, $|\hat{v}_{N}-v|\leq\epsilon=2g(r(N))/f_{X}(v)$ all $N\geq N_{0}$. ∎

If the set $\Omega_{x}$ is convex for all $x\in B(v,\delta)$, then $\Omega_{x}\in\mathcal{M}_{g}$ with $g(\epsilon)=5d\epsilon$ for $\epsilon$ small enough. This is because the volume $\lambda_{d}((\partial\Omega)_{\epsilon})$ for any convex set $\Omega\subset[0,1)^{d}$ is bounded by that of the case $\Omega=[0,1)^{d}$, which is no larger than $5d\epsilon$ for $\epsilon$ small enough (see the proof of Lemma 3.7). By Theorem 3.6, the deterministic error bound for the QMC-based quantile estimate becomes

[TABLE]

The result may be extended to pseudo-convex sets (see [23]). However, the convex conditions on $\Omega_{x}$ may be restrictive for practical problems. We next show that under the Lipschitz continuity condition on $\phi$, $g(\epsilon)=\kappa\epsilon$ for some constant $\kappa>0$. The same rate $O(N^{-1/d})$ also applies for this case.

Lemma 3.7**.**

Suppose that $\phi(\bm{u})$ is Lipschitz continuous over $[0,1)^{d}$ with modulus $L>0$. If Assumption 3.4 is satisfied, there exist $\epsilon_{0},\delta>0$ such that

[TABLE]

for any $\epsilon\leq\epsilon_{0}$.

Proof.

Let $S_{1}$ be the boundary of the unit cube $[0,1)^{d}$, and let $S_{2}=\{\bm{u}\in[0,1)^{d}|\phi(\bm{u})=x\}$. Note that $\partial\Omega_{x}\subset S_{1}\cup S_{2}$. As a result,

[TABLE]

Note that $\lambda_{d}((S_{1})_{\epsilon})\leq 2[(1+2\epsilon)^{d}-1]=4d\epsilon+O(\epsilon^{2})$. So there exists $\epsilon^{\prime}>0$ such that $\lambda_{d}((S_{1})_{\epsilon})\leq 5d\epsilon$ for any $\epsilon\leq\epsilon^{\prime}$. For any $\bm{v}\in(S_{2})_{\epsilon}\backslash(S_{1})_{\epsilon}$, there exists $\bm{w}\in S_{2}$ such that $||\bm{v}-\bm{w}||\leq\epsilon$. Since $\phi(\cdot)$ is Lipschitz, $|\phi(\bm{v})-\phi(\bm{w})|=|\phi(\bm{v})-x|\leq L||\bm{v}-\bm{w}||\leq L\epsilon$. Therefore,

[TABLE]

By Assumption 3.4, there exists $\delta>0$ such that $f_{X}(x)\leq(3/2)f_{X}(v)$ for all $x\in B(v,2\delta)$. Let $\epsilon_{0}=\min(\epsilon^{\prime},\delta/L)$. Assume that $\epsilon\leq\epsilon_{0}$ and $x\in B(v,\delta)$. By the mean value theorem, $\lambda_{d}(S_{3})=F(x+L\epsilon)-F(x-L\epsilon)=2L\epsilon f_{X}(\xi)$ for some $\xi\in B(x,L\epsilon)\subset B(v,2\delta)$. Using $f_{X}(\xi)\leq(3/2)f_{X}(v)$ gives $\lambda_{d}((\partial\Omega_{x})_{\epsilon})\leq 5d\epsilon+\lambda_{d}(S_{3})\leq(5d+3f_{X}(v)L)\epsilon$ for all $\epsilon\leq\epsilon_{0}$ and all $x\in B(v,\delta)$.

∎

Theorem 3.8**.**

Suppose that the point set $\{\bm{u}_{1},\dots,\bm{u}_{N}\}$ used in (2.10) is a $(t,m,d)$-net in base $b\geq 2$, where $N=b^{m}$. If Assumptions 3.4 is satisfied and $\phi(\bm{u})$ is Lipschitz continuous over $[0,1)^{d}$ with modulus $L>0$, then

[TABLE]

for $N$ large enough.

Proof.

By Lemma 3.7 and Theorem 3.6, there exists $N_{0}>0$ such that for $N\geq N_{0}$,

[TABLE]

∎

Assumption 3.4 is typically used in establishing the asymptotic normality of $\hat{v}_{N}$ in the MC setting. The Lipschitz continuity condition on $\phi$ can be easily verified for some applications; see Section 5 for the greater detail.

4. Convergence analysis for RQMC expected shortfall estimation

In this section, we study the MSE of the expected shortfall estimate $\hat{c}_{N}$ when using RQMC. Define $K(x)=\mathbb{E}[(x-X)^{+}]$, which is estimated by

[TABLE]

where $\bm{u}_{i}$ are given in (2.10). Note that $\hat{K}_{N}(x)$ can be viewed as a QMC quadrature rule $\hat{I}_{N}(f)$ and $K(x)=I(f)$ for the kink function $f(\bm{u})=(x-\phi(\bm{u}))^{+}$. Also, $c=v-K(v)/p$ and $\hat{c}_{N}=\hat{v}_{N}-\hat{K}_{N}(\hat{v}_{N})/p$. The following lemma gives a relationship between the quantile estimation error and the expected shortfall estimation error.

Lemma 4.1**.**

If Assumption 3.4 is satisfied, then

[TABLE]

where

[TABLE]

Proof.

Let $B_{N}=\hat{c}_{N}-c-[K(v)-\hat{K}_{N}(v)]/p$. By Equation (12) in [21], we find that

[TABLE]

Under Assumption 3.4, we have

[TABLE]

which completes the proof. ∎

Theorem 4.2**.**

Suppose that the point set $\{\bm{u}_{1},\dots,\bm{u}_{N}\}$ used in (2.10) and (2.11) is a scrambled $(t,m,d)$-net in base $b\geq 2$, where $N=b^{m}$. Suppose additionally that Assumptions 3.4 and 3.5 are satisfied and $(v-\phi(\bm{u}))^{+}\in L^{2}([0,1)^{d})$. Then for $N$ large enough,

[TABLE]

where

[TABLE]

and $g(\cdot)$ is given in Assumption 3.5. Particularly, $\mathbb{E}[(\hat{c}_{N}-c)^{2}]=o(1/N)$. If $\phi(\bm{u})$ is of BVHK and $g(\epsilon)=\kappa\epsilon$ for some constant $\kappa>0$, then

[TABLE]

for arbitrarily small $\epsilon>0$.

Proof.

It is known that a scrambled $(t,m,d)$-net is a $(t,m,d)$-net w.p.1 (see [14]). By combining Theorem 3.6 and Lemma 4.1, we have

[TABLE]

For any square-integrable integrands, as shown in [15], the scrambled net variance with sample size $N$ is $o(1/N)$. This implies $\mathrm{Var}[\hat{K}_{N}(v)]=o(1/N)$ and $\mathrm{Var}[\hat{F}_{N}(v)]=o(1/N)$, leading to $\mathbb{E}[(\hat{c}_{N}-c)^{2}]=o(1/N)$.

Since $g(\epsilon)=\kappa\epsilon$ for some constant $\kappa>0$, it is easy to see that $\partial\Omega_{v}$ admits ($d-1$)-dimensional Minkowski content (see [7]). By Theorem 3.5 in [7], we have $\mathrm{Var}[\hat{K}_{N}(v)]=O(N^{-1-1/(2d-1)+\epsilon})$ for arbitrarily small $\epsilon>0$. By Theorem 4.4 in [7], we have $\mathrm{Var}[\hat{F}_{N}(v)]=O(N^{-1-1/d})$. Consequently, $\mathbb{E}[(\hat{c}_{N}-c)^{2}]=O(N^{-1-1/(2d-1)+\epsilon})$.

∎

Remark 4.3**.**

The deterministic error bound for the quantile estimate established in Theorem 3.6 plays an important role in studying the MSE of the expected shortfall estimate. The convergence result in Theorem 3.1 does not help to bound the MSE. Observed from (4.2), the accuracy of the expected shortfall estimate depends strongly on the RQMC integration of the kink function $(v-\phi(\bm{u}))^{+}$. This implies that if the RQMC quadrature rule yields a faster rate of convergence for the kink function, one can expect a better performance of the expected shortfall estimate. However, if $\phi(\bm{u})$ is not of BVHK, it may be hard to predict the MSE rate for the function $(v-\phi(\bm{u}))^{+}$ unless using the worst-case rate $o(1/N)$.

5. Numerical Examples

A stochastic activity network (SAN) models the time to compute a project having activities with random durations and precedence constraints. Figure 1 shows an instance of SAN with $d=15$ activities, which correspond to the edges in the network. Dong and Nakayama [4] studied this model with LHS. Let $Y_{i}$ denote the time to complete the activity $i$. Assume that the activity durations $Y_{i}$ are independent exponential random variables $\mathrm{Exp}(\lambda_{i})$, i.e., the density of $Y_{i}$ is given by $p_{i}(x)=\lambda_{i}\exp(-\lambda_{i}x)\bm{1}\{x\geq 0\}$, where $\lambda_{i}>0$. The network in Figure 1 has $q=10$ paths form nodes $s$ to $t$, denoted by $B_{1},\dots,B_{q}$. Specially, $B_{1}=\{1,4,11,15\}$, $B_{2}=\{1,4,12\}$, $B_{3}=\{2,5,11,15\}$, $B_{4}=\{2,5,12\}$, $B_{5}=\{2,6,13\}$, $B_{6}=\{2,7,14\}$, $B_{7}=\{3,8,11,15\}$, $B_{8}=\{3,8,12\},$ $B_{9}=\{3,9,15\},\ B_{10}=\{3,10,14\}$. The time to complete the project can be modeled by the random variable

[TABLE]

We are interested in estimating the quantile of $X$. To simulate the model using QMC, we generate $Y_{j}=-(1/\lambda_{j})\log u_{j}$ for $j=1,\dots,d$. As a result, $X$ can be expressed a function of $\bm{u}=(u_{1},\dots,u_{d})$, denoted by $\phi(\bm{u})$. It should be noted that $\phi(\bm{u})$ is not Lipschitz continuous. So Theorem 3.8 cannot be applied directly. To circumvent this, we rewrite the set $\Omega_{x}$ as $\Omega_{x}=\cap_{i=1,\dots,q}A_{i}$, where

[TABLE]

Let $\lambda_{\mathrm{max}}=\max_{j=1,\dots,d}\lambda_{j}$, and let $\phi_{i}(\bm{u})=\prod_{j\in B_{i}}u_{j}^{\lambda_{\mathrm{max}}/\lambda_{j}}$. Then $A_{i}=\{\bm{u}\in[0,1)^{d}|\phi_{i}(\bm{u})\geq e^{-\lambda_{\mathrm{max}}x}\}$. It is easy to see that $\phi_{i}$ is Lipschitz continuous over $[0,1]^{d}$ because $\lambda_{\mathrm{max}}/\lambda_{j}\geq 1$ for all $j$. So by Lemma 3.7 and using $(\partial\Omega_{x})_{\epsilon}\subset\cup_{i=1,\dots,q}(\partial A_{i})_{\epsilon}$, the conditions in Theorem 3.6 are satisfied with $g(\epsilon)=\kappa\epsilon$ for some constant $\kappa>0$. The QMC error for the quantile estimation is $O(N^{-1/d})$.

We now study the MSE of the expected shortfall estimate $\hat{c}_{N}$ when using RQMC. By Theorem 4.4 in [7], we have $\mathrm{Var}[\hat{F}_{N}(v)]=O(N^{-1-1/d})$ since $g(\epsilon)=\kappa\epsilon$. Using (4.2) gives

[TABLE]

However, since $\phi(\bm{u})$ is not of BVHK, the rate $O(N^{-1-1/(2d-1)+\epsilon})$ established in Theorem 3.5 of [7] cannot be applied for the integrand $(v-\phi(\bm{u}))^{+}$. Instead, using the worst-case rate $\mathrm{Var}[\hat{K}_{N}(v)]=o(1/N)$ arrives at $\mathbb{E}[(\hat{c}_{N}-c)^{2}]=o(1/N)$. This confirms that RQMC performs asymptotically better than MC and LHS for expected shortfall estimation. The MSE rate $o(1/N)$ may be too conservative when $d$ is small.

Figure 2 shows the numerical results for the SAN model in Figure 1 with $\lambda_{i}=1/2$ for $i\leq 8$ and $\lambda_{i}=1$ for $i>8$. In the numerical experiments, we use Sobol’ points as inputs for QMC-based estimates and scrambled Sobol’ points for RQMC-based estimates. The MSEs in right panel of Figure 2 are computed based on 100 independent repetitions. Estimation of the errors requires knowing the true value of the quantity being estimated. Here we use the MC method with a very large sample size (say, $N=10^{9}$) to obtain accurate estimates of $v$ and $c$ and treat them as the true values. We consider the case $p=0.1$ for which the true values are $v=2.5446,\ c=2.1596.$ The empirical evidence shows convergence rates of (R)QMC beyond the crude MC rate of $N^{-1/2}$. Particularly, RQMC yields lower MSEs than MC for both the quantile and the expected shortfall estimations.

6. Conclusion

In this paper, we proved the convergence of QMC-based quantile estimates under very mild assumptions. More importantly, we proved that the QMC error is bounded from above by $N^{-1/d}$. The error rate $O(N^{-1/d})$ is worse than the usual MC rate $O(N^{-1/2})$ for $d>2$. But this error rate is critical to establish considerable MSE rates of RQMC for expected shortfall estimation. It is possible to obtain a faster error rate for RQMC-based quantile estimates as suggested by the numerical study. Owen [15] showed that scrambled net quadrature rules can yield an MSE of $o(N^{-1})$ for square-integrable functions. We conjecture that RQMC-based quantile estimation can also lead to an MSE of $o(N^{-1})$ under some technical conditions. We leave this problem open for future research.