remarkRemark hypothesisHypothesis claimClaim assumAssumption Efficient risk
estimation via nested multilevel QMC simulation Z. Xu, Z. He, and X. Wang

# Efficient risk estimation via nested multilevel quasi-Monte Carlo simulation
††thanks: Submitted to the editors DATE. This work was supported by the
National Key R&D Program of China (No. 2016QY02D0301) and the National Science
Foundation of China (No. 72071119, 12071154), and the Fundamental Research
Funds for the Central Universities (No. 2019MS106).

Zhenghang Xu Department of Mathematical Sciences, Tsinghua University, Beijing
100084, People’s Republic of China<EMAIL_ADDRESS><EMAIL_ADDRESS>Zhijian He Corresponding author. School of
Mathematics, South China University of Technology, Guangzhou 510641, People’s
Republic of China<EMAIL_ADDRESS>Xiaoqun Wang22footnotemark: 2

###### Abstract

We consider the problem of estimating the probability of a large loss from a
financial portfolio, where the future loss is expressed as a conditional
expectation. Since the conditional expectation is intractable in most cases,
one may resort to nested simulation. To reduce the complexity of nested
simulation, we present a method that combines multilevel Monte Carlo (MLMC)
and quasi-Monte Carlo (QMC). In the outer simulation, we use Monte Carlo to
generate financial scenarios. In the inner simulation, we use QMC to estimate
the portfolio loss in each scenario. We prove that using QMC can accelerate
the convergence rates in both the crude nested simulation and the multilevel
nested simulation. Under certain conditions, the complexity of MLMC can be
reduced to $O(\epsilon^{-2}(\log\epsilon)^{2})$ by incorporating QMC. On the
other hand, we find that MLMC encounters catastrophic coupling problem due to
the existence of indicator functions. To remedy this, we propose a smoothed
MLMC method which uses logistic sigmoid functions to approximate indicator
functions. Numerical results show that the optimal complexity
$O(\epsilon^{-2})$ is almost attained when using QMC methods in both MLMC and
smoothed MLMC, even in moderate high dimensions.

###### keywords:

nested simulation; quasi-Monte Carlo; multilevel Monte Carlo; risk estimation

65C05, 62P05

##  1 Introduction

We consider the problem of estimating

(1)
$\theta=\mathbb{P}(g(\omega)>c)=\mathbb{E}[\mathbb{I}\\{g(\omega)>c\\}]=\mathbb{E}[\mathbb{I}\\{\mathbb{E}[X|\omega]>c\\}]$

via simulation for a given constant $c$, where
$g(\omega):=\mathbb{E}[X|\omega]$. The inner expectation of the one-
dimensional random variable, $X$, is conditional on the value of the outer
multidimensional random variable $\omega$. This nested expectation appears,
for instance, when estimating the probability of a large loss from a financial
portfolio. If the portfolio consists of some complex financial derivatives
(such as path-dependent options and exotic options) or the underlying model is
complicated, the future loss $g(\omega)$ of the portfolio on a fixed period is
only known as a conditional expectation with respect to the risk factor
$\omega$, which does not have an analytical form of $\omega$. A typical
approach to estimate such an expectation of a function of a conditional
expectation is to use nested estimation. Specifically, nested simulation
refers to a two-level simulation procedure. In the outer level, one generates
a number of scenarios of $\omega$. Then, in the inner level, a number of
samples of $X$ are generated for each simulated $\omega$ to estimate the
conditional expectation $\mathbb{E}[X|\omega]$; see [2] and [14].

Nested simulation has been widely studied in the literature due to its broad
applicability, particularly in portfolio risk measurement. Gordy and Juneja
[14] considered and analyzed uniform nested simulation estimators which employ
a constant number of inner samples across all scenarios in the outer level.
They showed that to achieve a root mean squared error (RMSE) of $\epsilon$,
allocating proper computational effort in each level results in a total cost
of $O(\epsilon^{-3})$. Nested simulation can be made more efficient by
allocating computational effort nonuniformly across scenarios. Particularly,
Broadie et al. [2] proposed an adaptive procedure to allocate computational
effort to inner simulations. The resulting nonuniform nested simulation
estimator enjoys a reduced complexity of $O(\epsilon^{-5/2})$ under certain
conditions. Recently, Giles and Haji-Ali [10] proposed to use multilevel Monte
Carlo (MLMC) method in nested simulation. They showed that in the original
MLMC, the complexity is $O(\epsilon^{-5/2})$. By incorporating the adaptive
allocations procedure of [2], Giles and Haji-Ali [10] showed that the
complexity of MLMC can be reduced to $O(\epsilon^{-2}(\log\epsilon)^{2})$
under certain conditions. For some applications of nested MLMC, we refer to
[8, 9, 12, 13] and references therein.

MLMC is a sophisticated variance reduction technique introduced by Heinrich
[19] for parametric integration and by Giles [6] for the estimation of the
expectations arising from stochastic differential equations. Nowadays MLMC
method has been extended extensively. For a thorough review of MLMC methods,
we refer to [7]. On the other hand, quasi-Monte Carlo (QMC) and randomized QMC
(RQMC) methods are alternatives to improve the efficiency of traditional Monte
Carlo; see [5, 26, 28] for details. QMC methods have achieved great success in
finance applications, such as option pricing and hedging [25]. It is natural
to incorporate (R)QMC in the MLMC framework. Giles and Waterhouse [11] first
attempted to apply QMC method into multilevel path simulation in financial
problems. In recent years, multilevel QMC methods have received increasing
attention amongst researchers due to its broad applicability, particularly in
problems of partial differential equations with random coefficients, see e.g.,
[4, 22, 23] and in uncertainty quantification, see e.g., [3, 32].

In this paper, we focus on the combination of MLMC and QMC in nested
simulation. Specifically, in the outer simulation, we use Monte Carlo to
generate a number of scenarios of $\omega$. But, in the inner simulation, we
use QMC to estimate the portfolio loss in each scenario. Our work is closely
related to Goda et al. [13], who used nested multilevel RQMC method to deal
with the expected value of partial perfect information (EVPPI) problem. A
central problem of EVPPI is to estimate an expectation of the form
$\mathbb{E}[f(g(\omega))]$, where $g(\omega)=\mathbb{E}[X|\omega]$ and
$f(\cdot)$ is a continuous function. As Goda et al. [13] pointed out, the
multilevel RQMC estimator can achieve the optimal complexity of
$O(\epsilon^{-2})$ under some mild conditions. However, for the problem Eq. 1
considered in this paper, the performance function $f(\cdot)$ becomes an
indicator function $\mathbb{I}(\cdot>c)$, making it much harder than EVPPI for
MLMC algorithms [8]. Moreover, the antithetic MLMC estimator for EVPPI used in
[13] does not help to reduce the variance convergence rate in our setting
because the discontinuity in the indicator function violates the
differentiability requirements of the antithetic estimator [10]. As a result,
the efficiency of using the antithetic form in multilevel RQMC is subtle for
estimating Eq. 1. The gain of using RQMC in nested MLMC is also unclear. Can
the multilevel RQMC estimator achieve the optimal complexity of
$O(\epsilon^{-2})$ or the sub-optimal complexity of
$O(\epsilon^{-2}(\log\epsilon)^{2})$? It is well-known that the performance of
(R)QMC integration depends on the smoothness of the integrand and the
dimension of the problem [16, 17, 18, 35]. So it is natural to ask how these
factors affect the efficiency of the multilevel RQMC algorithm.

The main contribution of this work is to reduce the complexity of the MLMC
method for financial risk management by using QMC methods in the multilevel
scheme. We also discuss some considerations of using antithetic MLMC with
RQMC. In Section 2, we review QMC methods and nested simulation. The
complexity of the uniform nested simulation combined with QMC is analyzed.
After introducing the basic MLMC method, we develop a multilevel QMC procedure
for the problem Eq. 1. The effects of using the QMC method and the indicator
function on the nested simulation are discussed. In Section 3, we develop a
new smoothed coupling method, which aims to overcome the catastrophic coupling
caused by small differences between the “coarse” and “fine” estimates for the
conditional expectation. Some numerical experiments are performed in Section
4. Finally, we conclude this paper with some remarks in Section 5.

##  2 Nested simulation and Multilevel Monte Carlo

###  2.1 Nested simulation and quasi-Monte Carlo

Our goal is to estimate an expectation of a function of a conditional
expectation via nested simulation, with the problem Eq. 1 of interest. In the
uniform nested simulation, one takes

$\hat{\theta}_{n,m}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}\\{\hat{g}_{m}(\omega_{i})>c\\},$

and

(2) $\hat{g}_{m}(y)=\frac{1}{m}\sum_{j=1}^{m}X_{j}(y),$

where $\omega_{i}$ are independent and identically distributed (iid)
replications of $\omega$, and $X_{j}(y)$ are iid replications of $X$ given
$\omega=y$. The quadrature rule $\hat{g}_{m}(y)$ is used to estimate the
conditional expectation $g(y)=\mathbb{E}[X|\omega=y]$ in the inner simulation.

In this paper, we incorporate QMC methods within the inner simulation. Let
$\Omega$ be the possible scenarios of the random variable $\omega$. To fit MC
or QMC framework, we assume that given $\omega=y$, $X$ can be generated via
the mapping

(3) $X(y)=\psi(\bm{u};y),\ y\in\Omega,$

for some function $\psi$, where $\bm{u}\sim\mathbb{U}([0,1)^{d})$. As a
result, the inner estimator Eq. 2 is replaced by

(4) $\hat{g}_{m}(y)=\frac{1}{m}\sum_{j=1}^{m}\psi(\bm{u}_{j};y),$

and $g(y)=\int_{[0,1)^{d}}\psi(\bm{u};y)\mathrm{d}\bm{u}$. If
$\bm{u}_{1},\dots,\bm{u}_{m}$ are iid uniform points over $[0,1)^{d}$, the
approximation Eq. 4 refers to the MC method, which attains a probabilistic
convergence rate of $O(m^{-1/2})$. If $\bm{u}_{1},\dots,\bm{u}_{m}$ are QMC
points (known as low discrepancy points), which are deterministic points
chosen from $[0,1)^{d}$ and are more uniformly distributed than random points,
the approximation Eq. 4 refers to the QMC method.

The error of the QMC quadrature can be bounded by the Koksma-Hlawka inequality
(see [28])

(5) $|g(y)-\hat{g}_{m}(y)|\leq
V_{\mathrm{HK}}\big{(}\psi(\cdot;y)\big{)}D^{*}(\bm{u}_{1},\dots,\bm{u}_{m}),$

where $V_{\mathrm{HK}}\big{(}\psi(\cdot;y)\big{)}$ is the variation of the
integrand $\psi(\cdot;y)$ for given $y$ in the sense of Hardy and Krause which
measures the smoothness of $\psi(\cdot;y)$, and
$D^{*}(\bm{u}_{1},\dots,\bm{u}_{m})$ is the star discrepancy which measures
the uniformity of the points set $\\{\bm{u}_{1},\dots,\bm{u}_{m}\\}$. The
Koksma-Hlawka inequality Eq. 5 implies that for functions of finite variation,
the convergence rate of QMC approximation is determined by the factor
$D^{*}(\bm{u}_{1},\dots,\bm{u}_{m})$, which is of order $O(m^{-1}(\log
m)^{d})$ for low discrepancy points.

There are various constructions of QMC point sets in the literature [3], such
as digital nets and lattice rule point sets. In this paper, we use
$(t,d)$-sequences in base $b\geq 2$. In practice, one uses RQMC points for
ease of evaluating quadrature errors. RQMC retains the essential equi-
distribution structure, while allowing a statistical error estimation based on
independent replications. Moreover, RQMC is able to improve the rate of
convergence for smooth integrands, such as the Owen’s scrambling technique
[29]. For a survey on RQMC, we refer to [26]. It should be noted that, in both
MC and RQMC settings, $\mathbb{E}[\hat{g}_{m}(y)]=g(y)$.

######  Assumption 1.

There exist a constant $\eta\geq 1$ and a function $\sigma^{2}(\cdot)\geq 0$
such that for any $y\in\Omega$,

$\mathrm{Var}(\hat{g}_{m}(y))\leq\frac{\sigma^{2}(y)}{m^{\eta}}.$

In the MC setting, we can take $\eta=1$ and
$\sigma^{2}(y)=\mathrm{Var}(X|\omega=y)$. In the RQMC setting, we can expect a
larger value of $\eta$. Indeed, a digit scrambling of [29] applied to
$(t,d)$-sequences leads to a variance of $O(m^{-2}(\log m)^{2d})$ for
integrands of finite variation in the sense of Hardy and Krause, and a
variance of $O(m^{-3}(\log m)^{d-1})$ for smooth enough integrands [31].
Furthermore, for any square integrable integrand, the scrambled net variance
$\mathrm{Var}(\hat{g}_{m}(y))$ has a conservative upper bound
$M\mathrm{Var}(X|\omega=y)/m$, where the constant $M$ depends on $t,\ d$ and
$b$ [30]. As a result, it is reasonable to assume that $\eta\geq 1$ for RQMC.
If the function $\psi$ in Eq. 4 is sufficiently smooth, then $\eta\approx 2$
or even larger.

The next theorem is a generalization of Proposition 1 in [14].

######  Theorem 2.1.

Suppose that 1 is satisfied, and let $f(\cdot)$ be the probability density
function of $g(\omega)$. Assume the following:

  * •

The joint density $p_{m}(x,y)$ of $g(\omega)$ and
$m^{\eta/2}[\hat{g}_{m}(\omega)-g(\omega)]$ and partial derivatives
$(\partial/\partial x)p_{m}(x,y)$ and $(\partial^{2}/\partial
x^{2})p_{m}(x,y)$ exist for each $m$ and $(x,y)$.

  * •

For each $m\geq 1$, there exist functions $f_{i,m}(\cdot)$ such that

$\bigg{|}\frac{\partial^{i}}{\partial x^{i}}p_{m}(x,y)\bigg{|}\leq
f_{i,m}(y),i=0,1,2.$

In addition,

$\sup_{m\geq 1}\int\left|y\right|^{r}f_{i,m}(y)\mathrm{d}y<\infty,\ for\
i=0,1,2\ and\ 0\leq r\leq 4.$

Then the bias of the nested estimator $\hat{\theta}_{n,m}$ asymptotically
satisfies

(6)
$|\mathbb{E}[\hat{\theta}_{n,m}-\theta]|\leq\frac{|\Theta^{\prime}(c)|}{m^{\eta}}+O(m^{-3\eta/2}),$

where

$\Theta(c)=\frac{1}{2}f(c)\mathbb{E}[\sigma^{2}(\omega)|g(\omega)=c],$

and $\sigma^{2}(y)$ is given in 1.

######  Proof 2.2.

Let

$\theta_{m}:=\mathbb{E}(\hat{\theta}_{n,m})=\mathbb{P}(\hat{g}_{m}(\omega)>c).$

Note that

$\theta_{m}=\int_{\mathbb{R}}\int_{c-ym^{-\eta/2}}^{\infty}p_{m}(x,y)dxdy,$

so that

$\theta_{m}-\theta=\int_{\mathbb{R}}\int_{c-ym^{-\eta/2}}^{c}p_{m}(x,y)dxdy.$

Consider the Taylor series expansion of the density function $p_{m}(x,y)$ at
$x=c$,

$p_{m}(x,y)=p_{m}(c,y)+(x-c)\frac{\partial}{\partial
x}p_{m}(c,y)+\frac{(x-c)^{2}}{2}\frac{\partial^{2}}{\partial
x^{2}}p_{m}(\tilde{c},y),$

where $\tilde{c}$ is a number between $c$ and $x$. From this expansion and the
assumptions in the theorem, it follows that

(7)
$\theta_{m}-\theta=\int_{\mathbb{R}}\frac{y}{m^{\eta/2}}p_{m}(c,y)dy-\int_{\mathbb{R}}\frac{y^{2}}{2m^{\eta}}\frac{\partial}{\partial
x}p_{m}(c,y)dy+O(m^{-3\eta/2}).$

For the first term on the right hand side of Eq. 7, we have

$\displaystyle\int_{\mathbb{R}}\frac{y}{m^{\eta/2}}p_{m}(c,y)dy$
$\displaystyle=\frac{f(c)}{m^{\eta/2}}\int_{\mathbb{R}}y\cdot\frac{p_{m}(c,y)}{f(c)}dy$
$\displaystyle=\frac{f(c)}{m^{\eta/2}}\mathbb{E}[m^{\eta/2}(\hat{g}_{m}(\omega)-g(\omega))|g(\omega)=c].$

This term is zero because

$\mathbb{E}[m^{\eta/2}(\hat{g}_{m}(\omega)-g(\omega))|g(\omega)=c]=\mathbb{E}[\mathbb{E}[m^{\eta/2}(\hat{g}_{m}(\omega)-g(\omega))|g(\omega)=c,\omega]]=0.$

For the second term on the right hand side of Eq. 7, we have

$\displaystyle\int_{\mathbb{R}}\frac{y^{2}}{2m^{\eta}}\frac{\partial}{\partial
x}p_{m}(c,y)dy$
$\displaystyle=\frac{1}{2m^{\eta}}\frac{d}{dc}\int_{\mathbb{R}}y^{2}p_{m}(c,y)dy$
$\displaystyle=\frac{1}{2m^{\eta}}\frac{d}{dc}\bigg{[}f(c)\int_{\mathbb{R}}y^{2}\cdot\frac{p_{m}(c,y)}{f(c)}dy\bigg{]}$
$\displaystyle=\frac{1}{2m^{\eta}}\frac{d}{dc}\bigg{[}f(c)\mathbb{E}[m^{\eta}(\hat{g}_{m}(\omega)-g(\omega))^{2}|g(\omega)=c]\bigg{]}$
$\displaystyle=\frac{1}{2m^{\eta}}\frac{d}{dc}\bigg{[}f(c)\mathbb{E}[\mathbb{E}[m^{\eta}(\hat{g}_{m}(\omega)-g(\omega))^{2}|\omega]|g(\omega)=c]\bigg{]}$
$\displaystyle=\frac{1}{2m^{\eta}}\frac{d}{dc}\bigg{[}f(c)\mathbb{E}[\sigma^{2}(\omega)|g(\omega)=c]\bigg{]}$
$\displaystyle=\frac{\Theta^{\prime}(c)}{m^{\eta}}$

By (7), the proof is completed.

We now use Theorem 2.1 to analyze the mean squared error (MSE) of the nested
estimator. By Eq. 6, it is easy to see that

$\displaystyle\mathrm{Var}(\hat{\theta}_{n,m})$
$\displaystyle=\frac{1}{n}\mathrm{Var}(\mathbb{I}(\hat{g}_{m}(\omega_{1})>c))=\frac{\theta_{m}(1-\theta_{m})}{n}$
$\displaystyle=\frac{\theta(1-\theta)}{n}+\frac{(\theta_{m}-\theta)(1-\theta_{m})}{n}+\frac{\theta(\theta-\theta_{m})}{n}$
$\displaystyle\leq\frac{\theta(1-\theta)}{n}+O(n^{-1}m^{-\eta}).$

So we have

$\displaystyle\mathbb{E}[(\hat{\theta}_{n,m}-\theta)^{2}]$
$\displaystyle=\mathrm{Var}(\hat{\theta}_{n,m})+(\mathbb{E}[\hat{\theta}_{n,m}-\theta])^{2}$
$\displaystyle\leq\frac{\theta(1-\theta)}{n}+\frac{\Theta^{\prime}(c)^{2}}{m^{2\eta}}+O(m^{-5\eta/2})+O(n^{-1}m^{-\eta}).$

To achieve an RMSE of $\epsilon$, this suggests the optimal allocations
$n=O(\epsilon^{-2})$ and $m=O(\epsilon^{-1/\eta})$. The total computational
cost is then $O(\epsilon^{-2-1/\eta})$. In the MC setting where $\eta=1$, the
complexity is $O(\epsilon^{-3})$, see [14]. In RQMC setting, it is possible to
achieve $\eta\approx 2$, so that the complexity is around
$O(\epsilon^{-5/2})$. Broadie et al. [2] used adaptive allocations in the
inner simulation to reduce the nested MC down to $O(\epsilon^{-5/2})$. Giles
and Haji-Ali [10] showed that in the original MLMC, the complexity is also
$O(\epsilon^{-5/2})$. By incorporating the adaptive allocations idea of [2],
Giles and Haji-Ali [10] showed that the complexity of MLMC can be reduced to
$O(\epsilon^{-2}(\log\epsilon)^{2})$. We next focus on using the framework of
multilevel method to reduce the complexity of RQMC-based nested simulation.

###  2.2 MLMC estimators

We review briefly the basic idea of MLMC. Let
$P=\theta=\mathbb{I}\\{g(\omega)>c\\}$, which requires an infinite cost to
evaluate. Instead of dealing with $P$ directly, we consider a sequence of
random variables $P_{0},P_{1},\dots$ with increasing approximation accuracy to
$P$ but with increasing cost per sample. The $\ell$th level approximation of
$P$ is defined as $P_{\ell}=\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}$,
where $m_{\ell}=2^{\ell+\ell_{0}}$ in our setting and $\ell_{0}\geq 0$. Due to
the linearity of expectation, we have the following telescoping representation

(8)
$\mathbb{E}[P_{L}]=E[P_{0}]+\sum_{\ell=1}^{L}\mathbb{E}[P_{\ell}-P_{\ell-1}].$

Let $Y_{\ell}=P_{\ell}-P_{\ell-1}$ for $\ell\geq 1$ and $Y_{0}=P_{0}$, we
further have

(9) $\mathbb{E}[P_{L}]=\sum_{\ell=0}^{L}\mathbb{E}[Y_{\ell}].$

The MLMC method uses the above equality and estimates each term on the right
hand side of Eq. 9 independently. The resulting MLMC estimator is given by

$\hat{\theta}_{\mathrm{MLMC}}=\sum_{\ell=0}^{L}Z_{\ell}$

with

$Z_{\ell}=\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}Y_{\ell}^{(i)},$

where $Y_{\ell}^{(1)},\dots,Y_{\ell}^{(N_{\ell})}$ are iid replications of
$Y_{\ell}$ for $\ell=0,\dots,L$. Denote the variance and the computational
cost of $Y_{\ell}$ by $V_{\ell}$ and $C_{\ell}$, respectively. The MSE of
$\hat{\theta}_{\mathrm{MLMC}}$ is then given by

$\mathbb{E}[(\hat{\theta}_{\mathrm{MLMC}}-\theta)^{2}]=\sum_{\ell=0}^{L}\frac{V_{\ell}}{N_{\ell}}+(\mathbb{E}[P_{L}-\theta])^{2}.$

The total cost of $\hat{\theta}_{\mathrm{MLMC}}$ is
$C=\sum_{\ell=0}^{L}N_{\ell}C_{\ell}$.

The pioneering work by Giles [6] established the following theorem for MLMC.

######  Theorem 2.3.

If there are constants $\alpha,\beta,\gamma,c_{1},c_{2},c_{3}$ such that
$\alpha\geq\min(\beta,\gamma)/2$ and

  * •

$|\mathbb{E}[P_{\ell}-\theta]|\leq c_{1}m_{\ell}^{-\alpha}$,

  * •

$V_{\ell}\leq c_{2}m_{\ell}^{-\beta}$,

  * •

$C_{\ell}\leq c_{3}m_{\ell}^{\gamma}$,

then there exists a constant $c_{4}>0$ such that for any $\epsilon<1/e$, MLMC
estimator $\hat{\theta}_{\mathrm{MLMC}}$ has an MSE bound
$\mathbb{E}[(\hat{\theta}_{\mathrm{MLMC}}-\theta)^{2}]\leq\epsilon^{2}$ with a
total computational cost $C$ bounded by

$C\leq\begin{cases}c_{4}\epsilon^{-2},&\beta>\gamma,\\\
c_{4}\epsilon^{-2}(\log\epsilon)^{2},&\beta=\gamma,\\\
c_{4}\epsilon^{-2-(\gamma-\beta)/\alpha},&\beta<\gamma.\end{cases}$

The constants $\alpha$ and $\beta$ in Theorem 2.3 describe the rates of bias
and variance decreasing, respectively, which are usually called weak
convergence and strong convergence. The constant $\gamma$ controls the
increasing of computing budget for each sample.

Now we develop our nested multilevel QMC estimator. Here we use the following
form for coupling the two consecutive levels:
$Y_{0}=P_{0}=\mathbb{I}\\{\hat{g}_{m_{0}}(\omega)>c\\}$, and for $\ell\geq 1$,

$Y_{\ell}=P_{\ell}-P_{\ell-1}=\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\mathbb{I}\\{\hat{g}_{m_{\ell-1}}(\omega)>c\\},$

where

$\hat{g}_{m_{\ell}}(y)=\frac{1}{m_{\ell}}\sum_{j=1}^{m_{\ell}}\psi(\bm{u}_{j};y),$

and $\psi$ is given by Eq. 3.

It is easy to see that $\mathbb{E}[Y_{\ell}]=\mathbb{E}[P_{\ell}-P_{\ell-1}]$
for $\ell\geq 1$. Also, $C_{\ell}=O(m_{\ell})$, i.e., $\gamma=1$. Theorem 2.1
suggests that $\alpha=\eta\geq 1/2$ in Theorem 2.3.

We next study the variance of $Y_{\ell}$. Note that for $\ell\geq 1$, we have

(10) $\mathrm{Var}(Y_{\ell})\leq
2\mathrm{Var}(P_{\ell}-P)+2\mathrm{Var}(P_{\ell-1}-P),$

since $\mathrm{Var}(A+B)\leq 2[\mathrm{Var}(A)+\mathrm{Var}(B)]$ for any
random variables $A,B$ with finite variances. It thus suffices to study the
decay rate of $\mathrm{Var}(P_{\ell}-P)$, which can be described by value of
the constant $\beta$.

######  Assumption 2.

Assume that the density of the random variable $|g(\omega)-c|/\sigma(\omega)$,
where $\sigma(\omega)$ is the square root of $\sigma^{2}(\omega)$ in 1,
denoted by $\rho(\cdot)$, exists. Moreover, we assume that there exist
constants $\rho_{0}>0$ and $x_{0}>0$ such that $\rho(x)\leq\rho_{0}$ for all
$x\in[0,x_{0}]$.

The next theorem is a generalization of Proposition 2.2 in [10].

######  Theorem 2.4.

Suppose that 1 and 2 are satisfied. Then

$\mathrm{Var}(P_{\ell}-P)\leq\mathbb{E}[(P_{\ell}-P)^{2}]=O(m_{\ell}^{-\eta/2}).$

######  Proof 2.5.

Let us start from

$\displaystyle\mathbb{E}\bigg{[}(\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\mathbb{I}\\{g(\omega)>c\\})^{2}\big{|}\omega\bigg{]}$
$\displaystyle=\mathbb{P}\bigg{[}|\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\mathbb{I}\\{g(\omega)>c\\}|=1\big{|}\omega\bigg{]}$
$\displaystyle\leq\mathbb{P}\bigg{[}|\hat{g}_{m_{\ell}}(\omega)-g(\omega)|\geq|g(\omega)-c|\big{|}\omega\bigg{]}.$

By Chebyshev’s inequality, we have

$\displaystyle\mathbb{P}\bigg{[}|\hat{g}_{m_{\ell}}(\omega)-g(\omega)|\geq|g(\omega)-c|\big{|}\omega\bigg{]}$
$\displaystyle\leq\min\bigg{(}1,|g(\omega)-c|^{-2}\mathrm{Var}(\hat{g}_{m_{\ell}}(\omega)\big{|}\omega)\bigg{)}$
$\displaystyle\leq\min\bigg{(}1,|g(\omega)-c|^{-2}\frac{\sigma^{2}(\omega)}{m_{\ell}^{\eta}}\bigg{)}.$

Taking expectation over $\omega$ yields

$\displaystyle\mathbb{E}\bigg{[}(\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\mathbb{I}\\{g(\omega)>c\\})^{2}\bigg{]}$
$\displaystyle\leq\int_{0}^{\infty}\min(1,x^{-2}m_{\ell}^{-\eta})\rho(x)dx$
$\displaystyle\leq\rho_{0}\int_{0}^{\infty}\min(1,x^{-2}{m_{\ell}^{-\eta}})dx+\min(1,x_{0}^{-2}m_{\ell}^{-\eta})$
$\displaystyle\leq 2\rho_{0}m^{-\eta/2}_{\ell}+x^{-2}_{0}m^{-\eta}_{\ell}.$

Hence, the variance is $O(m^{-\eta/2}_{\ell})$.

Based on Eq. 10 and Theorem 2.4, we have
$\mathrm{Var}(Y_{\ell})=O(m_{\ell}^{-\beta})$ with $\beta=\eta/2$. Since
$\gamma=1$, by Theorem 2.3, the total cost of MLQMC thus becomes

(11) $C=\begin{cases}O(\epsilon^{-2}),&\eta>2,\\\
O(\epsilon^{-2}(\log\epsilon)^{2}),&\eta=2,\\\
O(\epsilon^{-3/2-1/\eta}),&\eta\in[1,2).\end{cases}$

It is a common strategy to use an antithetic form for coupling the consecutive
levels, i.e.,

$Y_{\ell}=\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\frac{1}{2}\left(\mathbb{I}\\{\hat{g}_{m_{\ell-1}}^{(1)}(\omega)>c\\}+\mathbb{I}\\{\hat{g}_{m_{\ell-1}}^{(2)}(\omega)>c\\}\right),$

where

$\hat{g}_{m_{\ell-1}}^{(i)}(\omega)=\frac{1}{m_{\ell-1}}\sum_{j=1+(i-1)m_{\ell-1}}^{im_{\ell-1}}\psi(\bm{u}_{j};\omega),\
i=1,2.$

Antithetic sampling can sometimes reduce variance apparently, see [9] and
[13]. But it should be cautious to adapt antithetic sampling when the
underlying function is discontinuous. It is required that

$\mathbb{P}[\hat{g}_{m_{\ell-1}}^{(1)}(\omega)>c|\omega]=\mathbb{P}[\hat{g}_{m_{\ell-1}}^{(2)}(\omega)>c|\omega]$

in order to ensure the telescoping representation Eq. 8 under RQMC scheme.
This can be achieved if the first half RQMC points
$\bm{u}_{1},\dots,\bm{u}_{m_{\ell-1}}$ have the same joint distribution as
that of the second half RQMC points
$\bm{u}_{m_{\ell-1}+1},\dots,\bm{u}_{m_{\ell}}$. As pointed out by Goda et al.
[13], the well-known explicitly constructed $(t,d)$-sequences in base $b=2$
satisfy this condition if using Owen’s scrambling method [29]. However,
antithetic sampling does not change the variance convergence rate in this
problem because the discontinuity in the indicator function violates the
differentiability requirements of the antithetic estimator.

##  3 Smoothed MLMC method

Recall that

$Y_{\ell}=P_{\ell}-P_{\ell-1}=\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}-\mathbb{I}\\{\hat{g}_{m_{\ell-1}}(\omega)>c\\}.$

The “fine” estimator $\hat{g}_{m_{\ell}}(\omega)$ and the “coarse” estimator
$\hat{g}_{m_{\ell-1}}(\omega)$ conditional the same scenario $\omega$ have
small differences in most situations, especially for large $\ell$, making
$Y_{\ell}$ different from zero only in a tiny proportion of the scenarios.
Moraes et al. [27] called this phenomenon as catastrophic coupling.

This may lead to a higher kurtosis, which is defined by

$\kappa=\frac{\mathbb{E}[(Y_{\ell}-\mathbb{E}[Y_{\ell}])^{4}]}{(\mathrm{Var}(Y_{\ell}))^{2}}.$

This is a challenge for both MLMC and MLQMC. As shown in Hammouda et al. [1],
the sample variance of $Y_{\ell}^{(1)},\dots,Y_{\ell}^{(N_{\ell})}$ has a
standard error of

$\sigma_{\ell}:=\frac{\mathrm{Var}(Y_{\ell})}{\sqrt{N_{\ell}}}\sqrt{\kappa-1+\frac{2}{N_{\ell}-1}}.$

And the high kurtosis makes it challenging to estimate $V_{\ell}$ accurately
in deeper levels, since there are less samples in deeper levels under MLMC
scheme [1].

In fact, the high-kurtosis phenomenon is inevitable due to the indicator
function: $P_{\ell}$ takes values in $\\{0,1\\}$, then
$Y_{\ell}=P_{\ell}-P_{\ell-1}$ ($\ell\geq 1$) takes values in $\\{-1,0,1\\}$.
As a result, the iid replications of $Y_{\ell}$,
$Y^{(1)}_{\ell},\dots,Y^{(N_{\ell})}_{\ell}$ take values in $\\{-1,0,1\\}$ as
well. The sample mean of $Y^{(1)}_{\ell},\dots,Y^{(N_{\ell})}_{\ell}$, denoted
by $\bar{Y}_{\ell}$, is usually negligible because of small differences
between $P_{\ell}$ and $P_{\ell-1}$. We find that the sample kurtosis

(12) $\displaystyle\hat{\kappa}$
$\displaystyle=\frac{\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}(Y_{\ell}^{(i)}-\bar{Y}_{\ell})^{4}}{(\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}(Y_{\ell}^{(i)}-\bar{Y}_{\ell})^{2})^{2}}$
(13)
$\displaystyle\approx\frac{\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}{Y^{(i)}_{\ell}}^{4}}{(\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}{Y^{(i)}_{\ell}}^{2})^{2}}=\frac{1}{\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}{Y^{(i)}_{\ell}}^{2}}\approx\frac{1}{\mathcal{S}^{2}_{\ell}},$

where
$\mathcal{S}^{2}_{\ell}=\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}(Y^{(i)}_{\ell}-\bar{Y}_{\ell})^{2}$
is the sample variance. So the kurtosis increases with the same rate as the
variance decreases.

It should be noted that this phenomenon is due to the structure of the
indicator function. Moreover, the indicator function eliminates further the
difference between $\hat{g}_{m_{\ell}}(\omega)$ and
$\hat{g}_{m_{\ell-1}}(\omega)$, then this aggravates the catastrophic
coupling. This phenomenon enlightens us to replace the indicator function with
a smoother function. To this end, we introduce the logistic sigmoid function

$S(x)=\frac{1}{1+\exp(-x)},\ x\in\mathbb{R}.$

Logistic sigmoid function is widely applied in traditional classification
task, acting as activation function in deep learning, see e.g. [15] and [24].
Based on this function, we construct a family of sigmoid-like functions

$S(x;k)=\frac{1}{1+\exp(-kx)},\ x\in\mathbb{R},$

where $k>0$. It is obvious that $S(x;k)$ has the following properties:

  * •

$S(x;k)$ takes value in $(0,1)$ for any $k$;

  * •

$S(x;k)$ is centrosymmetric about $S(0;k)=\frac{1}{2}$ for any $k$;

  * •

$S(x;k)$ converges to $\mathbb{I}\\{x>0\\}$ as $k\rightarrow\infty$ for any
fixed $x$.

Fig. 1 presents several examples of $S(x;k)$. It can be seen that the sigmoid-
like functions approximate the indicator function. When the value of $k$ gets
larger, the approximation of $S(x;k)$ to $\mathbb{I}(x>0)$ gets better.

Figure 1: $S(x;k)$ with different $k$s.

Now we are ready to develop a smoothed MLMC method. Define the $\ell$th level
approximation of $P=\mathbb{I}\\{g(\omega)>c\\}$ as

$\tilde{P}_{\ell}=S^{(k_{0},r)}_{\ell}(\hat{g}_{m_{\ell}}(\omega)-c),$

instead of $P_{\ell}=\mathbb{I}\\{\hat{g}_{m_{\ell}}(\omega)>c\\}$ used in
Section 2.2, where $S^{(k_{0},r)}_{\ell}(x)=S(x;k_{0}r^{\ell})$, $k_{0}>0$ is
a constant and $r>1$ controls the rate of $S^{(k_{0},r)}_{\ell}(x)$ converging
to $\mathbb{I}\\{x>0\\}$ as $\ell\to\infty$.

Similar to Eq. 8, we have the following telescoping representation

$\mathbb{E}[\tilde{P}_{L}]=E[\tilde{P}_{0}]+\sum_{\ell=1}^{L}\mathbb{E}[\tilde{P}_{\ell}-\tilde{P}_{\ell-1}].$

Let $\tilde{Y}_{\ell}=\tilde{P}_{\ell}-\tilde{P}_{\ell-1}$ for $\ell\geq 1$
and $\tilde{Y}_{0}=\tilde{P}_{0}$. It is obvious that

$\tilde{P}_{L}=\sum_{\ell=0}^{L}\tilde{Y}_{\ell},$

and

$\mathbb{E}[\tilde{P}_{L}]=\sum_{\ell=0}^{L}\mathbb{E}[\tilde{Y}_{\ell}].$

This yields a new MLMC estimator

$\tilde{\theta}_{\mathrm{MLMC}}=\sum_{\ell=0}^{L}\tilde{Z}_{\ell}$

with

$\tilde{Z}_{\ell}=\frac{1}{N_{\ell}}\sum_{j=1}^{N_{\ell}}\tilde{Y}_{\ell}^{(i)},$

where $\tilde{Y}_{\ell}^{(1)},\dots,\tilde{Y}_{\ell}^{(N_{\ell})}$ are iid
replications of $\tilde{Y}_{\ell}$ for $\ell=0,\dots,L$. Notice that
$\tilde{P}_{L}$ converges to $P$ as $L$ goes to infinity, just like $P_{L}$
does. Denoting the variance and cost of each sample of $\tilde{Y}_{\ell}$ by
$\tilde{V}_{\ell}$ and $\tilde{C}_{\ell}$, the MSE of
$\tilde{\theta}_{\mathrm{MLMC}}$ is given by

$\mathbb{E}[(\tilde{\theta}_{\mathrm{MLMC}}-\theta)^{2}]=\sum_{\ell=0}^{L}\frac{\tilde{V}_{\ell}}{N_{\ell}}+(\mathbb{E}[\tilde{P}_{L}-\theta])^{2}.$

The total cost of $\tilde{\theta}_{\mathrm{MLMC}}$ is
$\tilde{C}=\sum_{\ell=0}^{L}N_{\ell}\tilde{C}_{\ell}$. It is critical to
verify the strong convergence. To this end, we make the following assumption.

######  Assumption 3.

Let $\Pi_{\ell}(x)$ denote the the cumulative distribution function of the
random variable $\hat{g}_{m_{\ell}}(\omega)-c$. Assume that the densities of
these variables, denoted by $\pi_{\ell}(x)$ for $\ell=0,1,\dots$, exist in a
common neighborhood of the origin, say $(-\delta_{0},\delta_{0})$. Assume also
that there exists a maximum
$\pi^{*}_{\ell}=\sup_{x\in(-\delta_{0},\delta_{0})}\pi_{\ell}(x)<+\infty$ for
all $\ell$ and there is a constant $C_{\pi}$ such that
$\sup_{\ell}\pi_{\ell}^{*}<C_{\pi}$.

The first part of 3 is similar to the assumptions in Theorem 2.1, but here we
only require the densities to exist in a neighborhood of the origin. For the
later part of 3, heuristically, we have
$\lim_{\ell\rightarrow\infty}\pi_{\ell}(x)=\pi(x)$, where $\pi(x)$ is density
of $g(\omega)$. Then we can expect that $\pi_{\ell}^{*}$ can be controlled by
a constant $C_{\pi}$.

######  Theorem 3.1.

Suppose that the conditions of 1, 2 and 3 are satisfied. Then

$\mathrm{Var}(\tilde{P}_{\ell}-P)=O(\max(r^{-\ell},m_{\ell}^{-\eta/2})).$

######  Proof 3.2.

First, we have

$\displaystyle\mathbb{E}[(\tilde{P}_{\ell}-P_{\ell})^{2}]$
$\displaystyle=\int_{-\infty}^{-\delta_{0}}[S_{\ell}^{(k_{0},r)}(x)]^{2}d\Pi_{\ell}(x)+\int_{-\delta_{0}}^{0}[S_{\ell}^{(k_{0},r)}(x)]^{2}\pi_{\ell}(x)dx$
$\displaystyle+\int_{0}^{\delta_{0}}[S_{\ell}^{(k_{0},r)}(-x)]^{2}\pi_{\ell}(x)dx+\int_{\delta_{0}}^{+\infty}[S_{\ell}^{(k_{0},r)}(-x)]^{2}d\Pi_{\ell}(x)$
$\displaystyle\leq
2[S_{\ell}^{(k_{0},r)}(-\delta_{0})]^{2}+2\pi^{*}_{\ell}\int_{-\delta_{0}}^{0}[S_{\ell}^{(k_{0},r)}(x)]^{2}dx$
$\displaystyle=\frac{2}{[1+\exp(k_{0}r^{\ell}\delta_{0})]^{2}}+\bigg{[}\ln\frac{2}{1+e^{-k_{0}r^{\ell}\delta_{0}}}+\frac{1}{2}-\frac{1}{1+e^{-k_{0}r^{\ell}\delta_{0}}}\bigg{]}\cdot\frac{2C_{\pi}}{k_{0}r^{\ell}}$
$\displaystyle=O(r^{-\ell}).$

By Theorem 2.4, we have $\mathbb{E}[(P_{\ell}-P)^{2}]=O(m_{\ell}^{-\eta/2})$.
It then follows that

$\mathbb{E}[(\tilde{P}_{\ell}-P)^{2}]\leq
2\mathbb{E}[(\tilde{P}_{\ell}-P_{\ell})^{2}]+2\mathbb{E}[(P_{\ell}-P)^{2}]=O(\max(r^{-\ell},m_{\ell}^{-\eta/2})).$

Similarly, we can verify the weak convergence of the smoothed MLMC method.

######  Theorem 3.3.

Suppose that the conditions of Theorem 2.1 and 3 are satisfied. Then

(14)
$|\mathbb{E}[\tilde{P}_{\ell}-\theta]|=O(\max(r^{-\ell},m_{\ell}^{-\eta})).$

######  Proof 3.4.

Note that

(15)
$|\mathbb{E}[\tilde{P}_{\ell}-\theta]|\leq\mathbb{E}[|\tilde{P}_{\ell}-P_{\ell}|]+|\mathbb{E}[P_{\ell}-\theta]|.$

By Theorem 2.1, $|\mathbb{E}[P_{\ell}-\theta]|=O(m^{-\eta}_{\ell})$. Under 3,
we have

$\displaystyle\mathbb{E}[|\tilde{P}_{\ell}-P_{\ell}|]$
$\displaystyle=\int_{-\infty}^{-\delta_{0}}S_{\ell}^{(k_{0},r)}(x)d\Pi_{\ell}(x)+\int_{-\delta_{0}}^{0}S_{\ell}^{(k_{0},r)}(x)\pi_{\ell}(x)dx$
$\displaystyle+\int_{0}^{\delta_{0}}S_{\ell}^{(k_{0},r)}(-x)\pi_{\ell}(x)dx+\int_{\delta_{0}}^{+\infty}S_{\ell}^{(k_{0},r)}(-x)d\Pi_{\ell}(x)$
$\displaystyle\leq
2S_{\ell}^{(k_{0},r)}(-\delta_{0})+2\pi^{*}_{\ell}\int_{-\delta_{0}}^{0}S_{\ell}^{(k_{0},r)}(x)dx$
$\displaystyle=\frac{2}{1+\exp(k_{0}r^{\ell}\delta_{0})}+\ln\frac{2}{1+\exp(-k_{0}r^{\ell}\delta_{0})}\cdot\frac{2C_{\pi}}{k_{0}r^{\ell}}$
$\displaystyle=O(r^{-\ell}).$

The result Eq. 14 follows from Eq. 15.

For a given $\eta$, it is sufficient to take $r=2^{\eta/2}$, and then
$O(r^{-\ell})\leq O(m_{\ell}^{-\eta/2})$, since $m_{\ell}=O(2^{\ell})$. As a
result, the strong convergence of the smoothed method is the same as the
original ML(Q)MC methods while the weak convergence may slow down to
$O(m^{-\eta/2}_{\ell})$ because of the smoothed method. When $\eta=2$, we can
take $r=2$ such that there are constants $c_{2},c_{3}$ such that
$\tilde{V}_{\ell}\leq c_{2}m_{\ell}^{-1}$ and $\tilde{C}_{\ell}\leq
c_{3}m_{\ell}^{-1}$, which means $\beta=\gamma=1$, then the cost is in the
second regime in Theorem 2.3, where weak convergence makes little difference.
So $r=2$ is sufficient for this situation to get the expected improvement.

The next theorem can be obtained immediately from Theorem 2.3.

######  Theorem 3.5.

Suppose that 1 is satisfied with $\eta=2$ and $r=2$. Then for any
$\epsilon<1/e$, the smoothed MLMC estimator $\tilde{\theta}_{\mathrm{MLMC}}$
has an MSE bound
$\mathbb{E}[(\tilde{\theta}_{\mathrm{MLMC}}-\theta)^{2}]\leq\epsilon^{2}$ with
a total computational cost of $O(\epsilon^{-2}(\log\epsilon)^{2})$.

Theorem 3.1 and Theorem 3.3 guarantee that the smoothed MLMC method converges
to the true value under a certain rate. However, it should be noticed that
$r=2$ is sufficient but unnecessary for the smoothed MLMC method. In fact,
when the dimension is relatively high, we would like to take $r$ a relatively
small value. On the one hand, the efficiency of QMC method is affected by high
dimensionality leading an $\eta$ value smaller than 2. As a result, a smaller
$r$ can also match to $O(m_{\ell}^{-\eta/2})$ rate. On the other hand, if the
value of $r$ is smaller, the coupling between the “fine” estimator
$\tilde{P}_{\ell}$ and the “coarse” estimator $\tilde{P}_{\ell-1}$ is better
and then this leads to a smaller variance. By this way, we can get a better
$\beta$.

In addition, the logistic sigmoid function $S_{\ell}^{(k_{0},r)}(x)$ with a
smaller $r$ can maintain milder derivative in the neighborhood of the origin
which relieves the catastrophic coupling. Antithetic sampling can also benefit
from the differentiability and mild derivative of the logistic sigmoid
functions with small $r$. If $r$ is taken too large, $S_{\ell}^{(k_{0},r)}(x)$
converges to indicator function rapidly, making the smoothed MLMC method
inefficient. So in our numerical studies, we take $r=2$ in one-dimensional
problems, while take $r=\sqrt{2}$ in multiple-dimensional problems. We take
$k_{0}=8$ in
$\tilde{P}_{\ell}=S^{(k_{0},r)}_{\ell}(\hat{g}_{m_{\ell}}(\omega)-c)$, which
guarantees the sigmoid functions not differing from the indicator function too
much at the beginning.

Finally, we define kurtosis and variance factor (KVF) of level $\ell$ as the
product of kurtosis and variance in level $\ell$,

$\textrm{KVF}=\kappa V_{\ell}.$

KVF can be estimated by

$\widehat{\textrm{KVF}}=\hat{\kappa}\mathcal{S}^{2}_{\ell},$

where $\hat{\kappa}$ is defined by Eq. 12, and
$\mathcal{S}^{2}_{\ell}=\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}(Y^{(i)}_{\ell}-\bar{Y}_{\ell})^{2}$.
This factor describes the relationship between the rates of variance
decreasing and kurtosis increasing, which can reveal the effectiveness of a
method dealing with the high-kurtosis phenomenon. For the original MLMC and
MLQMC methods, the KVF is very close to 1 for mildly large levels as explained
in Eq. 13. Therefore, if a method can achieve a smaller KVF, then with a same
$V_{\ell}$, the kurtosis is not so large, which means the high-kurtosis
phenomenon is relieved and the estimation of $V_{\ell}$ can be more accurate.

##  4 Numerical study

In the numerical study we consider a portfolio that consists of European
options written on $d$ stocks whose price dynamics follow the Black-Scholes
model. For simplicity, we assume that the stock returns are the same, denoted
by $\mu$, while risk-free interest rate is $\mu_{0}$. Price dynamics of the
stocks $\bm{S}_{t}=(S^{1}_{t},\dots,S^{d}_{t})$ evolve according to

$\frac{dS^{i}_{t}}{S_{t}^{i}}=\mu^{\prime}dt+\sum_{j=1}^{d}\sigma_{ij}dW^{i}_{t},\
i=1,\dots,d,$

where $\mu^{\prime}=\mu$ under the real-world probability measure, and
$\mu^{\prime}=\mu_{0}$ under the risk-neutral probability measure. Here,
$\bm{W}_{t}=(W^{1}_{t},\dots,W^{d}_{t})$ is a $d$-dimensional standard
Brownian motion which represents $d$ risk factors in the model. We have,

$S^{i}_{t}=S^{i}_{0}\exp\\{(\mu^{\prime}-\frac{1}{2}\sum_{j=1}^{d}\sigma^{2}_{ij})t+\sum_{j=1}^{d}\sigma_{ij}W_{t}^{j}\\},\
i=1,\dots,d,$

where $\bm{S}_{0}=(S^{1}_{0},\dots,S^{d}_{0})$ are the initial prices of
stocks.

We assume that the maturities of all the European options in the portfolio are
the same, denoted by $T$. We want to measure the portfolio risk at a future
time $\tau$ ($\tau<T$). In the simulation, we first simulate the random
variable $\omega=\bm{S}_{\tau}=(S^{1}_{\tau},\dots,S^{d}_{\tau})$ under real-
world probability measure, which denotes the prices of stocks at the risk
horizon $\tau$ as the outer sample. And then simulate
$\bm{S}_{T}=(S^{1}_{T},\dots,S^{d}_{T})$ under the risk-neutral probability
measure which is the prices of stocks at maturity $T$ given $\omega$ as the
inner samples. Denote by $V_{0}=\sum_{i=1}^{d}v^{i}_{0}$ the initial value of
the portfolio, where each $v^{i}_{0}$ is known by the Black-Scholes formula
[21]. Then the portfolio value change is

$g(\omega):=V_{0}-\mathbb{E}[V_{T}(\bm{S}_{T})|\omega],$

where $V_{T}(\bm{S}_{T})$ is the discounted payoff of the portfolio at time
$T$ which is a known function of $\bm{S}_{T}$.

Our target is to estimate the loss probability
$\theta=\mathbb{P}[g(\omega)>c]$ for a given threshold $c$. We take
$m_{\ell}=32\times 2^{\ell}$ in all experiments and $k_{0}=8$ for the smoothed
methods. We compare MLMC, MLQMC, and smoothed MLQMC (denoted by SMLQMC) in
each example.

###  4.1 Single asset

Firstly, we consider a portfolio consists of a single put option, i.e., $d=1$.
This example was studied in [2].

In the simulation, the outer random variable is generated according to

$\omega=S_{\tau}=S_{0}\exp\\{(\mu-\sigma^{2}/2)\tau+\sigma\sqrt{\tau}Z\\},$

where the real-valued risk factor $Z$ is a standard normal random variable.
The portfolio value change is

$g(\omega)=v_{0}-\mathbb{E}[e^{-\mu_{0}(T-\tau)}(K-S_{T}(\omega,W))^{+}|\omega],$

where the expectation is taken over the random variable $W$, which is a
standard normal random variable independent with $Z$, and $S_{T}(\omega,W)$ is
given by

$S_{T}(\omega,W)=\omega\exp\\{(\mu_{0}-\sigma^{2}/2)(T-\tau)+\sigma\sqrt{T-\tau}W\\}.$

Now let $X=v_{0}-e^{-\mu_{0}(T-\tau)}(K-S_{T}(\omega,W))^{+}$. For $\omega=y$,
$X$ can be generated via

$X(y)=\psi(u;y)=v_{0}-e^{-\mu_{0}(T-\tau)}(K-y\exp\\{(\mu_{0}-\sigma^{2}/2)(T-\tau)+\sigma\sqrt{T-\tau}\Phi^{-1}(u)\\})^{+}.$

Note that $\psi(u;y)$ is a decreasing function with respect to $u$. So the
variation of Hardy and Krause for $\psi(u;y)$ can be computed easily, i.e.,
for any $y$,

$V_{\mathrm{HK}}(\psi(\cdot;y))=e^{-\mu_{0}(T-\tau)}K.$

As a result, when using scrambled $(t,1)$-sequence in base $b=2$ in the inner
simulation, we have

$\mathrm{Var}(\hat{g}_{m}(y))\leq\frac{C}{m^{2}},$

where $C$ is a constant that does not depend on $y$, and the sample size has
the form $m=2^{\ell}$. 1 is thus satisfied with $\eta=2$ and $\sigma(y)\equiv
C.$

Note that $g(y)$ is a strictly increasing continuous function. So the
cumulative distribution function of $g(\omega)$ can be computed easily,
namely,

$\displaystyle\mathbb{P}[g(\omega)\leq t]$
$\displaystyle=\mathbb{P}[\omega\leq g^{-1}(t)]$
$\displaystyle=\mathbb{P}[S_{0}\exp\\{(\mu-\sigma^{2}/2)\tau+\sigma\sqrt{\tau}Z\\}\leq
g^{-1}(t)]$
$\displaystyle=\Phi([\log(g^{-1}(t)/S_{0})-(\mu-\sigma^{2}/2)\tau]/(\sigma\sqrt{\tau})).$

The density of $g(\omega)$ is then given by

$f(t)=\frac{\phi([\log(g^{-1}(t)/S_{0})-(\mu-\sigma^{2}/2)\tau]/(\sigma\sqrt{\tau}))}{\sigma\sqrt{\tau}g^{-1}(t)g^{\prime}(g^{-1}(t))}.$

It is easy to see that $f(0)<\infty$ so that 2 is satisfied. The MLQMC thus
gives a complexity of $O(\epsilon^{-2}(\log\epsilon)^{2})$ by Eq. 11.

Figure 2: Comparison of MLMC, MLQMC and SMLQMC for the single asset case.

The parameters in this case are set as follows: $S_{0}=100$, $\mu=8\%$,
$\mu_{0}=3\%$, the annualized volatility $\sigma=20\%$, the strike of the put
option $K=95$, the maturity $T=0.25$ years (i.e., three months), and the risk
horizon $\tau=1/52$ years (i.e., one week). With these parameters, the initial
value of the put option is $v_{0}=1.669$ by the Black–Scholes formula. The
inner expectation $g(\omega)$ can be computed explicitly by Black-Scholes
formula. We choose $c=0.476887$ so that the loss probability $\theta=0.3$. For
this one-dimensional problem, we take $r=2$ in the smoothed MLQMC. The results
on testing the weak convergence and the strong convergence in Fig. 2 are based
on 500,000 outer samples in each level.

Fig. 2 displays the weak convergence of the three methods, from which we can
see that both crude MLQMC and smoothed MLQMC method result in larger $\alpha$
than the MLMC method. This implies that QMC methods enjoy a better weak
convergence and converge to the true value faster.

Fig. 2 shows the comparison of strong convergence, which plays a central role
in MLMC methods. It is obvious that QMC methods work much better than MLMC,
yielding the slopes of the lines lower than $-1$. Thus, QMC methods make the
constant $\beta$ exceed $\gamma=1$, which means the variance $V_{\ell}$
decreases faster than the cost $C_{\ell}$ increases. And the smoothed MLQMC
method shares the same rate with the crude MLQMC method but achieves an even
smaller variance. As a result, the costs in each level are reduced
accordingly, which is confirmed by Fig. 2.

As we stated earlier, QMC methods are able to accelerate the strong
convergence, leading to $\beta>\gamma=1$ in Theorem 2.3. Then the total cost
falls in the third regime in Theorem 2.3, which is the optimal complexity
bound $O(\epsilon^{-2})$ for MLMC methods. The numerical results support this
in Fig. 2. Both QMC methods are of the best complexity, while the smoothed
MLQMC method reduces further the computational burden.

Fig. 2 displays the tendency of kurtosis. As we explained, the kurtosis
increases with the variance decreasing as $\ell$ gets larger. This phenomenon
is inevitable due to the structure of the indicator function. But combined
with the smoothed method, the increasing of kurtosis slows down obviously,
which means the smoothed MLQMC method takes effect on high-kurtosis
phenomenon.

More visually, in Fig. 2, we can see the KVFs of the MLMC and the crude MLQMC
methods are close to 1 in all the cases, while the smoothed MLQMC method
yields much smaller KVF. This implies that the smoothed method slows down the
rate of kurtosis increasing. In this sense, we can see the effect of smoothed
MLQMC on overcoming the high-kurtosis phenomenon and catastrophic coupling.

###  4.2 Multiple assets

Now we consider a portfolio consists of $d$ European call options, which was
studied in [20]. Given the outer sample
$\omega=\bm{S}_{\tau}=(S^{1}_{\tau},\dots,S^{d}_{\tau})$ which denotes the
prices of stocks at the risk horizon $\tau$ under real-world measure, and
$\bm{K}=(K^{1},\dots,K^{d})$, the strike prices for call options, the
portfolio value change is

$g(\omega)=V_{0}-\mathbb{E}[e^{-\mu_{0}(T-\tau)}\sum_{i=1}^{d}(S^{i}_{T}(\omega^{i},\bm{W})-K^{i})^{+}|\omega],$

where $\omega^{i}$ is the $i$th element of the vector $\omega$ and the
expectation is taken over the random variable $\bm{W}=(W^{1},\dots,W^{d})\sim
N(\bm{0},\bm{I}_{d})$, and samples of

$S^{i}_{T}(\omega^{i},\bm{W})=\omega^{i}\exp\\{(\mu_{0}-\frac{1}{2}\sum_{j=1}^{d}\sigma^{2}_{ij})(T-\tau)+\sum_{j=1}^{d}\sigma_{ij}\sqrt{T-\tau}W^{j}\\}$

are simulated under the risk-neutral measure.

To handle the effect of high dimensionality, we combine the gradient PCA
(GPCA) method proposed by Xiao and Wang [34] within SMLQMC, which is
abbreviated as GMLQMC hereafter. GPCA is a general method to reduce the
effective dimension of functions for improving the efficiency of QMC method.
Here we apply GPCA in the inner simulation. We also examine the antithetic
sampling of GMLQMC (abbreviated as AMLQMC). As we know, using an antithetic
estimator does not change the variance convergence rate in the crude MLMC
because of the indicator function. It is of interest to test whether
antithetic sampling can get more benefit in the smoothed method.

The parameters in our experiments are set as follows:
$S^{1}_{0}=\dots=S^{d}_{0}=100$, $\mu=8\%$, $\mu_{0}=5\%$, the strikes
$K^{1}=\dots=K^{d}=95$, the maturity $T=0.1$, and the risk horizon
$\tau=0.02$. Without loss of generality, we let $\Sigma=(\sigma_{ij})$ be a
sub-triangular matrix satisfying $C=\Sigma\Sigma^{T}$ which corresponds with
Cholesky decomposition of $C$, where $C_{ij}=0.3\cdot 0.98^{|i-j|}$. The
different decompositions of $C$ do not influence the efficiency of MC method,
but take effect under QMC scheme, while the Cholesky decomposition is a common
and standard method in QMC simulation. The threshold $c=20\%V_{0}$. We take
$r=\sqrt{2}$ in the smoothed MLQMC methods.

(a) $d=4$ (b) $d=16$ (c) $d=32$

Figure 3: Estimations of $|\mathbb{E}{Y_{\ell}}|$ or
$|\mathbb{E}{\tilde{Y}_{\ell}}|$.

Firstly, we perform the convergence tests for estimating $\alpha$ and $\beta$
in Theorem 2.3. The results are based on 300,000 outer samples in each level.
Fig. 3 shows the behaviors of the absolute values of the expectations of
$Y_{\ell}$ for crude methods and $\tilde{Y}_{\ell}$ for smoothed methods. It
can be seen that QMC methods have smaller values but with more volatility.
This is caused by the catastrophic coupling mentioned in Section 2, because
there is only a tiny proportion of samples different from zero. The high
sensitivity makes it difficult to predict the value of $\alpha$, even with
such a size of outer samples.

Figure 4: Estimations of $|\mathbb{E}{Y_{\ell}}|$ for MLMC.

Not only MLQMC methods are faced with this difficulty, so the MLMC method is.
Fig. 4 shows the results of the MLMC method in deep levels. Under the same
magnitude, both MLMC and MLQMC have large volatility. But it can be observed
in Fig. 3 that the rates of the MLQMC methods are not worse than the MLMC
method and the value of $\alpha$ does not matter here actually for the strong
convergence being good enough. Nevertheless, it should be noted that the
smoothed MLQMC methods are affected less by the catastrophic coupling and look
more stable.

(a) $d=4$ (b) $d=16$ (c) $d=32$

Figure 5: Estimations of $\mathrm{Var}(Y_{\ell})$ or
$\mathrm{Var}(\tilde{Y}_{\ell})$.

Fig. 5 shows the empirical variances of $Y_{\ell}$ and $\tilde{Y}_{\ell}$ for
different levels, which can be used to predict the strong convergence rate
$\beta$ by the usual linear regression. For the plain MLMC, we observe
$\beta\approx 0.5$ for all dimensions. For the crude MLQMC method, the
dimension $d$ has an impact on $\beta$, and $\beta\approx 1$ for moderately
large $d=32$. When combined with GPCA method in the MLQMC method (i.e., the
GMLQMC method), we observe a larger $\beta$, usually $1.1$. When antithetic
sampling method is applied in smoothed MLQMC method, $\beta$ is further
improved to $1.2$. These insights suggest that antithetic sampling can benefit
from the smooth coupling, achieving a larger $\beta$. The strong convergence
gets apparent improvement with QMC methods, in which $\beta=\gamma$ even
$\beta>\gamma$. The total cost then is $O(\epsilon^{-2}(\log\epsilon)^{2})$
even $O(\epsilon^{-2})$.

(a) $d=4$ (b) $d=16$ (c) $d=32$

Figure 6: Tests of KVF.

Similar to the case of $d=1$, the kurtosis also becomes larger as the level
$\ell$ increases. Here we focus on KVF to take both variance and kurtosis into
consideration; see Fig. 6 for the results. It should be noticed that the curve
of AMLQMC showed is 4 times KVF of AMLQMC actually, in order to balance the
influence bought by the antithetic sampling. This is because the antithetic
sampling can reduce the variance and kurtosis by 1/2 for the plain MLMC and
crude MLQMC without changing the convergence rate, then it is fairer to
compare 4 times KVF of the method with antithetic sampling. The MLMC and the
MLQMC methods without the smooth coupling have KVFs close to 1. On the other
hand, the smoothed methods achieve faster strong convergence and slower
kurtosis increasing, and so smaller values of KVF are observed. It can be seen
that the KVF gets smaller in deeper levels, which means the smoothed MLQMC
methods make the kurtosis increase slowly. In that sense, the smoothed MLQMC
helps to overcome the catastrophic coupling and make the algorithm more
efficient.

(a) $d=4$ (b) $d=16$ (c) $d=32$

Figure 7: Tests of total cost.

Fig. 7 shows the total computation cost of each method. We can see that the
complexity bound of MLQMC methods is $O(\epsilon^{-2})$ or
$O(\epsilon^{-2}\log(\epsilon)^{2})$ as expected. The smooth coupling indeed
helps to reduce the total cost.

(a) $d=32$ (b) $d=64$

Figure 8: Tests of dimension effects.

As shown in Fig. 5 and Fig. 7, the crude MLQMC method is good enough, almost
attaining the optimal complexity, even for the large dimension $d=32$. The
effect of the dimension $d$ on the performance of the crude MLQMC method seems
minor. A possible explanation is that the integrand in the inner simulation
may enjoy low effective dimension, which is friendly to QMC method [33]. It is
clear that the choice of the matrix $\Sigma$ has an impact on the effective
dimension of the integrand.

Now we examine another risk factors structure, with $\Sigma$ satisfying
$\tilde{C}=\Sigma\Sigma^{T}$, where $\tilde{C}_{ij}=0.3(d-|i-j|)/d$. Fig. 8
shows the decays of the variances of $Y_{\ell}$ and $\tilde{Y}_{\ell}$ with
the new matrix $\Sigma$ for $d=32$ and $d=64$. For these cases, the strong
convergence rate $\beta$ of the crude MLQMC methods is not much better than
the MLMC method. But when it is combined with the dimension reduction
technique GPCA, the strong convergence is improved dramatically, in which
$\beta$ exceeds $1$. We thus conclude that the dimension of the integrand in
the inner simulation has a significant impact on the performance of MLQMC
methods. As the dimension gets higher, plain QMC methods may render a small
$\eta<2$ in 1. When combined with dimension reduction methods, such as the
GPCA method, the efficiency of QMC methods can be reclaimed, achieving a
better strong convergence. This highlights the importance of using dimension
reduction methods in nested MLQMC.

##  5 Conclusion

In this paper, we have incorporated randomized QMC methods into MLMC to deal
with financial risk estimation via nested simulation. We have proved that the
new MLQMC estimator can achieve better complexity bound under some
assumptions. At the meantime, we discussed the catastrophic coupling
phenomenon caused by the character of RQMC points and discontinuity of the
indicator function in the problem. Then we developed a new smoothed MLMC
method and we proved that the smoothed MLMC method still reserves the
advantages of MLMC without extra requirement. The smoothed method can also
take advantage of antithetic sampling, which does not work for the original
method. The superiority of the (smoothed) MLQMC methods over the original MLMC
method has been empirically shown in numerical studies.

Further improvements for the MLQMC methods can be expected in the follow
respects. On the one hand, we only applied RQMC methods in the inner
simulation. If RQMC methods are used not only in the inner sampling but also
in the outer sampling, the computational cost may be further reduced. And the
smooth coupling estimator is more suitable for RQMC methods and an improved
efficiency can be expected. On the other hand, we take a fixed number of inner
samples in an identical level in this paper. If we are able to allocate the
sample size according to the outer samples adaptively under the QMC scheme,
like the way in the Monte Carlo scheme [2, 10], a further improvement can be
also expected.

## Acknowledgments

ZH would like to appreciate Prof. Micheal B. Giles and Dr. Abdul-Lateef Haji-
Ali for useful discussions and comments at an early stage of this research.

## References

  * [1] C. Ben Hammouda, N. Ben Rached, and R. Tempone, Importance sampling for a robust and efficient multilevel Monte Carlo estimator for stochastic reaction networks, Stat. Comput., 30 (2020), pp. 1665–1689.
  * [2] M. Broadie, Y. Du, and C. C. Moallemi, Efficient risk estimation via nested sequential simulation, Manag. Sci., 57 (2011), pp. 1172–1194, https://www.jstor.org/stable/25835766.
  * [3] J. Dick, R. N. Gantner, Q. T. Le Gia, and C. Schwab, Multilevel higher-order quasi-Monte Carlo Bayesian estimation, Math. Models Methods Appl. Sci., 27 (2017), pp. 953–995, https://doi.org/10.1142/S021820251750021X.
  * [4] J. Dick, F. Y. Kuo, Q. T. L. Gia, and C. Schwab, Multilevel higher order QMC Petrov-Galerkin discretization for affine parametric operator equations, SIAM J. Numer. Anal., 54 (2016), pp. 2541–2568, https://doi.org/10.1137/16M1078690.
  * [5] J. Dick, F. Y. Kuo, and I. H. Sloan, High-dimensional integration: The quasi-Monte Carlo way, Acta Numer., 22 (2013), pp. 133–288, https://doi.org/10.1017/S0962492913000044.
  * [6] M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res., 56 (2008), pp. 607–617, https://www.jstor.org/stable/25147215.
  * [7] M. B. Giles, Multilevel Monte Carlo methods, Acta Numer., 24 (2015), pp. 259–328, https://doi.org/10.1017/S096249291500001X.
  * [8] M. B. Giles, MLMC for nested expectations, in Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan, Springer, 2018, pp. 425–442, https://doi.org/10.1007%2F978-3-319-72456-0_20.
  * [9] M. B. Giles and T. Goda, Decision-making under uncertainty: using MLMC for efficient estimation of EVPPI, Stat. Comput., 29 (2019), pp. 739–751, https://doi.org/10.1007/s11222-018-9835-1.
  * [10] M. B. Giles and A.-L. Haji-Ali, Multilevel nested simulation for efficient risk estimation, SIAM/ASA J. Uncertain. Quantif., 7 (2019), pp. 497–525, https://doi.org/10.1137/18M1173186.
  * [11] M. B. Giles and B. J. Waterhouse, Multilevel quasi-Monte Carlo path simulation, Advanced Financial Modelling, (2009), pp. 165–181.
  * [12] T. Goda, T. Hironaka, and T. Iwamoto, Multilevel Monte Carlo estimation of expected information gains, Stoch. Anal. Appl., 38 (2020), pp. 581–600, https://doi.org/10.1080/07362994.2019.1705168.
  * [13] T. Goda, D. Murakami, K. Tanaka, and K. Sato, Decision-theoretic sensitivity analysis for reservoir development under uncertainty using multilevel quasi-Monte Carlo methods, Comput. Geosci., 22 (2018), pp. 1009–1020, https://doi.org/10.1007/s10596-018-9735-7.
  * [14] M. B. Gordy and S. Juneja, Nested simulation in portfolio risk measurement, Manag. Sci., 56 (2010), pp. 1833–1848, https://www.jstor.org/stable/40864742.
  * [15] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data mining, Inference, and Prediction, New York: Springer, 2nd ed., 2016.
  * [16] Z. He, Quasi-Monte Carlo for discontinuous integrands with singularities along the boundary of the unit cube, Math. Comp., 87 (2018), pp. 2857–2870, https://doi.org/10.1090/mcom/3324.
  * [17] Z. He, On the error rate of conditional quasi-Monte Carlo for discontinuous functions, SIAM J. Numer. Anal., 57 (2019), pp. 854–874, https://doi.org/10.1137/18M118270X.
  * [18] Z. He and X. Wang, On the convergence rate of randomized quasi-Monte Carlo for discontinuous functions, SIAM J. Numer. Anal., 53 (2015), pp. 2488–2503, https://doi.org/10.1137/15M1007963.
  * [19] S. Heinrich, Monte Carlo complexity of global solution of integral equations, J. Complexity, 14 (1998), pp. 151–175, https://doi.org/10.1006/jcom.1998.0471.
  * [20] L. J. Hong, S. Juneja, and G. Liu, Kernel smoothing for nested estimation with application to portfolio risk measurement, Oper. Res., 65 (2017), pp. 657–673.
  * [21] J. C. Hull, Options, Futures, and Other Derivatives, Pearson, 10th ed., 2017.
  * [22] F. Y. Kuo, R. Scheichl, C. Schwab, I. H. Sloan, and E. Ullmann, Multilevel quasi-Monte Carlo methods for lognormal diffusion problems, Math. Comp., 86 (2017), pp. 2827–2860, https://doi.org/10.1090/mcom/3207.
  * [23] F. Y. Kuo, C. Schwab, and I. H. Sloan, Multi-level quasi-Monte Carlo finite element methods for a class of elliptic PDEs with random coefficients, Found. Comput. Math., 15 (2015), pp. 411–449, https://doi.org/10.1007/s10208-014-9237-5.
  * [24] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature (London), 521 (2015), pp. 436–444, https://doi.org/10.1038/nature14539.
  * [25] P. L’Ecuyer, Quasi-Monte Carlo methods with applications in finance, Finance Stoch., 13 (2009), pp. 307–349, https://doi.org/10.1007/s00780-009-0095-y.
  * [26] P. L’Ecuyer and C. Lemieux, Recent advances in randomized quasi-Monte Carlo methods, in Modeling Uncertainty, Springer, 2002, pp. 419–474, https://doi.org/10.1007/0-306-48102-2_20.
  * [27] A. Moraes, R. Tempone, and P. Vilanova, Multilevel hybrid Chernoff tau-leap, BIT, 56 (2016), pp. 189–239, https://doi.org/10.1007/s10543-015-0556-y.
  * [28] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, Society for Industrial and Applied Mathematics, 1992.
  * [29] A. B. Owen, Randomly permuted $(t,m,s)$-nets and $(t,s)$-sequences. in: Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing (H. Niederreiter and P. J.-S. Shiue, Eds.), J. Complexity, (1995), pp. 299–317, https://doi.org/10.1007/978-1-4612-2552-2_19.
  * [30] A. B. Owen, Scrambled net variance for integrals of smooth functions, Ann. Statist., 25 (1997), pp. 1541–1562, https://www.jstor.org/stable/2959062.
  * [31] A. B. Owen, Scrambling Sobol’ and Niederreiter-Xing points, J. Complexity, 14 (1998), pp. 259–328, https://doi.org/10.1006/jcom.1998.0487.
  * [32] R. Scheichl, A. M. Stuart, and A. L. Teckentrup, Quasi-Monte Carlo and multilevel Monte Carlo methods for computing posterior expectations in elliptic inverse problems, SIAM/ASA J. Uncertain. Quantif., 5 (2017), pp. 493–518, https://doi.org/10.1137/16M1061692.
  * [33] X. Wang, On the effects of dimension reduction techniques on some high-dimensional problems in finance, Oper. Res., 54 (2006), pp. 1063–1078.
  * [34] Y. Xiao and X. Wang, Enhancing quasi-Monte Carlo simulation by minimizing effective dimension for derivative pricing, Comput. Econom., 54 (2017), pp. 343–366, https://doi.org/10.1007/s10614-017-9732-2.
  * [35] F. Xie, Z. He, and X. Wang, An importance sampling-based smoothing approach for quasi-Monte Carlo simulation of discrete barrier options, European J. Oper. Res., 274 (2019), pp. 759–772, https://doi.org/10.1016/j.ejor.2018.10.030.


# Jet-medium photons as a probe of parton dynamics ††thanks: Talk presented at
the XXIXth International Conference on Ultra-relativistic Nucleus- Nucleus
Collisions (Quark Matter 2022)

Rouzbeh Modarresi Yazdia Speaker Shuzhe Shib Charles Galea Sangyong Jeona a
Department of Physics, McGill University, 3600 University Street, Montreal,
QC, Canada H3A 2T8
bCenter for Nuclear Theory, Department of Physics and Astronomy, Stony Brook
University, Stony Brook, New York 11794–3800, USA

###### Abstract

Photons resulting from jet-medium interactions offer the opportunity of
studying the evolving quark distribution in a heavy ion collision. The spectra
of jet-medium photons is presented within the jetscape framework for two
different energy loss models, martini and cujet. Jet-medium photons can
contribute significantly to the spectrum of direct photons in the intermediate
$p_{T}$ range.

_Introduction_ — Photons, real and virtual, are an important probe of the
quark gluon plasma. They are produced at all stages of the evolution of the
plasma and have a mean-free-path larger than the size of the created medium.
As such, their study and measurement can be a great tool for quantifying the
properties of QGP as well as jet-energy loss. Here we present a recent
calculation of photons produced from jet interactions with the QGP using the
jetscape [1] framework. To that comprehensive suite, we have added cujet [2]
as a low virtuality energy loss module, thus allowing for a systematic
comparison with martini [3]. What follows is a brief description of the
formalisms, a discussion of jet-medium photons and finally the results of the
simulation, a discussion, and outlook on future work.

_Low Virtuality Energy Loss_ — cujet implements the dglv [4, 5] inelastic
parton splitting rates, computed to leading order in the opacity expansion
[2]. The calculation of the gluon emission rates assumes that the hard probe
was created at some given time and position inside the plasma and then
evolved, scattering from a dynamical medium [2] and radiated a gluon. The
calculation of dglv rates is performed in the eikonal limit, with the jet and
the radiated gluon assumed to be distinguishable from the strongly interacting
medium around them. The Landau–Pomeranchuk–Migdal (LPM) effect manifests
itself as a phase factor. Soft, collinear divergences are regulated by a gluon
plasmon mass and the light quarks are considered massless [5]. It is also
assumed that the radiated gluon does not modify the direction of travel of the
incoming parton. In cujet collisional energy loss occurs via scattering with
the QGP medium. This is done using the Thoma–Gyulassy model where the HTL
gluon propagator includes a natural IR regulator [6]. Finally, cujet allows
for the running of the strong coupling, $\alpha_{\mathrm{s}}$, according to
the one-loop pQCD expression.

martini implements the AMY-McGill formalism for radiative and elastic
scattering energy loss on-shell energetic partons in a strongly interacting
medium [3]. It solves a rate equation with gain and loss terms for the time-
evolving parton distribution. As in cujet, radiative energy loss is strictly
collinear, with elastic scattering channels providing the momentum broadening
via space-like exchanges with the medium. Unlike the dglv model, the original
AMY radiative rates, which we use in this work, do not have an explicit time
dependence [3]. The reaction rates calculated with AMY are to all orders of
the opacity expansion and account fully for the LPM effect. Collisional energy
loss is also implemented, again using the gluon HTL propagator. Other than
gluon radiation and collisional energy loss, martini also includes $g\to
q+\bar{q}$ radiative channel as well as conversion processes $q\to g$ and
$g\to q$. The strong coupling in radiative and elastic channels is allowed to
run, using the one-loop expression for $\alpha_{\mathrm{s}}$.

_Jet-Medium Photons_ — Two dominant mechanisms of jet-medium photon emission
are jet bremsstrahlung and jet-photon conversion [7]. Jet bremsstrahlung has a
similar structure as the jet gluon bremsstrahlung channel: the jet receives
kicks from the medium and emits a photon. The differential rate for this
process in martini is given by

$\displaystyle\frac{\mathrm{d}\Gamma^{\textsc{amy}}_{q\to
q\gamma}}{\mathrm{d}z}(p,z)=\,\frac{e_{f}^{2}\alpha_{\mathrm{em}}P_{q\to
q\gamma}(z)}{[2p\,z(1{-}z)]^{2}}$ $\displaystyle[1-f_{q}((1-z)p)]$
$\displaystyle\times\int\,\frac{\mathrm{d}^{2}\mathbf{k}_{\perp}}{(2\pi)^{2}}~{}\mathrm{Re}\left[2\,\mathbf{k}_{\perp}\cdot\mathbf{g}_{(z,p)}(\mathbf{k}_{\perp})\right],$
(1)

where $e_{f}$ is the fractional charge of quark(anti-quark) and $z\equiv k/p$
is the momentum fraction carried away by the photon.
$\mathbf{g}_{(z,p)}(\mathbf{k}_{\perp})$ is the solution to the integral
equation[8]

$2\,\mathbf{k}_{\perp}=i\,\delta
E(z,p,\mathbf{k}_{\perp})\mathbf{g}(\mathbf{k}_{\perp})+\int\frac{\mathrm{d}^{2}\mathbf{q}_{\perp}}{(2\pi)^{2}}~{}{C}(\mathbf{q}_{\perp})\Big{[}\mathbf{g}(\mathbf{k}_{\perp})-\mathbf{g}(\mathbf{k}_{\perp}-z\,\mathbf{q}_{\perp})\Big{]}\
$ (2)

where $q_{\perp}$ and $k_{\perp}$ are the exchanged and photon transverse
momenta, respectively, $C(\mathbf{q}_{\bot})$ is a scattering kernel, and
$\delta\,E$ is the energy difference between the initial and final states[8].

Conversion photons arise from soft momentum exchanges of $q(\bar{q})$ jets
with the QGP medium. At leading order in the strong coupling, the process can
go through either QCD Compton scattering or $q\bar{q}$ annihilation. Both
processes include a Mandelstam-t channel where a soft quark is exchanged. A
general property of t-channel processes is the preference of the outgoing
particle to be in the same direction as the incoming owing to $1/t$ term in
the matrix element and the matrix element is maximum when the relative angle
is very nearly zero. As such, we can perform an approximation where the
entirety of the contribution to the conversion photon rate stems from this
collinear region and write the rate as

$\frac{\mathrm{d}\Gamma^{\gamma\mathrm{-Conv.}}}{\mathrm{d}p\,\mathrm{d}k}(p,k,T)=e^{2}_{f}\frac{2\pi\,\alpha\,\alpha_{s}}{3}\frac{T}{p}\left[\frac{1}{2}\ln{\frac{2pT}{m^{2}_{q}}-0.36149}\right]\,\delta\left(p-k\right).$
(3)

where the collinearity is made explicit by the delta function.

Figure 1: Calculation of (a) charged hadron and (b) jet $R_{\mathrm{AA}}$ for
martini and cujet for $\mathrm{Pb+Pb}$ collisions at $\sqrt{s}=2.76$ ATeV.
Data from [9, 10, 11, 12, 13].

_Results_ — Our calculations were performed using the jetscape framework
simulating a $\mathrm{Pb+Pb}$ collision at $\sqrt{s}=2.76$ ATeV. Two sets of
simulations were made using cujet and martini used as the low-virtuality
energy loss modules with all other parameters held fixed. pythia was used to
generate the hard scattering event, performing the initial state shower as
well as handling multi-parton interactions. High virtuality partons were then
handed over to matter[14, 15], to be further evolved down in virtuality. For
the $\mathrm{pp}$ spectra used in the $R_{\mathrm{AA}}$ calculation and the
tuning of the other modules in jetscape Ref. [1] was used. The hydrodynamic
background was provided by a vishnu $\left(2+1\right)$D viscous hydrodynamic
simulation with temperature dependent specific shear and bulk viscosities
[16]. Finally fastJet [17] was used for the jet-clustering.

Fig. 1 (a) and (b) show the results of the simulation for charged hadrons and
jets, respectively, for the $0$-$5\%$ centrality class. We see satisfactory
agreement between the two energy loss modules and the data. Despite the
theoretical differences between them, the two energy loss models are nearly
indistinguishable. We now turn our attention to photons. Fig. 2 provides a
channel-by-channel calculation of photon yield at mid-rapidity for the same
colliding system as before (at $20$-$40\%$ centrality) and compared to data
from ALICE. The prompt, pre-equilibrium and thermal photons are taken from
Ref. [18].

Figure 2: Channel break down of the direct photon yield at midrapidity of
$20$-$40\%$ collisions of $\mathrm{Pb+Pb}$ at $\sqrt{s}=2.76$ ATeV. Data from
[19].

It is clear from Fig. 2 that the significance of medium sources of photons
(pre-equilibrium and thermal) goes down while that of the prompt photons goes
up as we consider larger and larger values of photon $p_{T}$. The jet-medium
channels contribute $\approx 30\%$ at intermediate values of photon $p_{T}$,
peaking at $\approx 7$ GeV. Fig. 3 makes the effect of including jet-medium
photons even more stark. Including photons from jet-medium interactions brings
the overall curve into much better agreement with the data, particularly for
the intermediate $p^{\gamma}_{T}\in\left[5,8\right]$ GeV.

Figure 3:  Transverse momentum spectrum of direct photons, with and without
jet-medium contribution. Jet-medium photons are calculated using martini.
Error bars and shaded regions are statistical and systematic uncertainties,
respectively.

To compare the equivalent photon yield from cujet we require a cujet
calculation of bremsstrahlung photons , which is currently being completed.
Thus we restrict the comparison to conversion photons in Fig. 4 for three
centrality classes. The more peripheral collisions result in a smaller medium
and lower temperatures. Given the temperature dependence of the rates of
conversion photons, it is natural to see the $30$-$40\%$ curves to be nearly
identical. The $0$-$5\%$ centrality class shows larger difference between the
two modules: conversion photon yield from martini is nearly $40\%$ larger than
that of cujet. The direct proportionality of conversion photons to the
underlying $q,\bar{q}$ distribution, then, indicates a difference in the two
modules and how they modify the evolving $q/\bar{q}$ distribution.

Figure 4:  Comparison of conversion photon spectra from cujet and martini, for
3 centrality classes of $\mathrm{Pb+Pb}$ at $\sqrt{s}=2.76$ ATeV. See text for
details.

_Conclusion & Outlook_ — Jet energy loss and suppression of charged hadron
yield relative to the proton-proton baseline is an important signal of the
creation of QGP. In this work we presented the results of the implementation
of cujet into jetscape and a comparison of charged hadron and jet
$R_{\mathrm{AA}}$ from cujet and martini. We further presented the first
dynamic calculation of jet-medium photons in a realistic plasma and a dynamic
initial jet distribution with martini. The inclusion of jet-medium photons in
Fig. 3 made a clear contribution to the total photon yield, bringing the
theoretical curve into better agreement with the experimental measurement. We
saw that the two jet-medium channels, bremsstrahlung and conversion photons,
together contribute nearly $\approx 30\%$ to the total photon yield for
intermediate values of photon transverse momentum. The comparison of cujet and
martini in photon yield is currently limited to conversion photons only, since
the former does not contain bremsstrahlung photons. We found that the
difference in conversion photon yield from the two modules can be significant:
$40\%$ in central ($0$-$5\%$) collisions which further motivates the usage of
jet-medium photons as clean probes of energy loss models.

## References

  * [1] A. Kumar et al.  JETSCAPE framework: $p+p$ results.  Phys. Rev. C, 102(5):054906, 2020.
  * [2] Jiechen Xu et al.  Azimuthal jet flavor tomography with CUJET2.0 of nuclear collisions at RHIC and LHC.  JHEP, 08:063, 2014.
  * [3] Bjoern Schenke et al.  MARTINI: An Event generator for relativistic heavy-ion collisions.  Phys. Rev. C, 80:054913, 2009.
  * [4] M. Gyulassy et al.  Reaction operator approach to non-Abelian energy loss.  Nucl. Phys. B, 594:371–419, 2001.
  * [5] Magdalena Djordjevic et al.  Heavy quark radiative energy loss in QCD matter.  Nucl. Phys. A, 733:265–298, 2004.
  * [6] Markus H. Thoma et al.  Quark damping and energy loss in the high temperature qcd.  Nuclear Physics B, 351(3):491–506, 1991.
  * [7] Rainer J. Fries et al.  High-energy photons from passage of jets through quark gluon plasma.  Phys. Rev. Lett., 90:132301, 2003.
  * [8] Peter B. Arnold et al.  Photon emission from ultrarelativistic plasmas.  JHEP, 11:057, 2001.
  * [9] Betty Abelev et al.  Centrality Dependence of Charged Particle Production at Large Transverse Momentum in Pb–Pb Collisions at $\sqrt{s_{\rm{NN}}}=2.76$ TeV.  Phys. Lett. B, 720:52–62, 2013.
  * [10] Georges Aad et al.  Measurement of charged-particle spectra in Pb+Pb collisions at $\sqrt{{s}_{\mathsf{{NN}}}}=2.76$ TeV with the ATLAS detector at the LHC.  JHEP, 09:050, 2015.
  * [11] Serguei Chatrchyan et al.  Study of high-pT charged particle suppression in PbPb compared to $pp$ collisions at $\sqrt{s_{NN}}=2.76$ TeV.  Eur. Phys. J. C, 72:1945, 2012.
  * [12] Georges Aad et al.  Measurements of the Nuclear Modification Factor for Jets in Pb+Pb Collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV with the ATLAS Detector.  Phys. Rev. Lett., 114(7):072302, 2015.
  * [13] Vardan Khachatryan et al.  Measurement of inclusive jet cross sections in $pp$ and PbPb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV.  Phys. Rev. C, 96(1):015202, 2017.
  * [14] Abhijit Majumder.  Incorporating Space-Time Within Medium-Modified Jet Event Generators.  Phys. Rev. C, 88:014909, 2013.
  * [15] Shanshan Cao et al.  Nuclear modification of leading hadrons and jets within a virtuality ordered parton shower.  Phys. Rev. C, 101(2):024903, 2020.
  * [16] Jonah E. Bernhard et al.  Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma.  Nature Phys., 15(11):1113–1117, 2019.
  * [17] Matteo Cacciari et al.  FastJet User Manual.  Eur. Phys. J. C, 72:1896, 2012.
  * [18] Charles Gale et al.  Multimessenger heavy-ion collision physics.  Phys. Rev. C, 105(1):014909, 2022.
  * [19] Jaroslav Adam et al.  Direct photon production in Pb-Pb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV.  Phys. Lett. B, 754:235–248, 2016.