3.1. Bias of TSB-Alg

An application of Theorem 2.1 implies that the bias of TSB-Alg equals

(6) \begin{equation}\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]-\mathbb{E}_{\boldsymbol0}[g(\chi_n)\Upsilon_{\boldsymbol\lambda}]=\mathbb{E}_{\boldsymbol0}\big[\Delta_n^g\big], \quad \text{where $\Delta_n^g\,:\!=\,(g(\chi)-g(\chi_n))\Upsilon_{\boldsymbol\lambda}$}.\end{equation}

The natural coupling $\big(\chi,\chi_n, Y_T^{{(+)}}, Y_T^{{(-)}}\big)$ in (6) is defined by $Y_T^{{(\pm)}}\,:\!=\,\sum_{k=1}^{\infty}\xi_k^{{(\pm)}}$ , $\xi_k\,:\!=\,\xi_k^{{(+)}}+\xi_k^{{(-)}}+ \eta_k$ for all $k\in\mathbb{N}$ , $\chi$ in (2), and $\chi_n$ in (5) with $Y_{L_n}\,:\!=\,\sum_{k=n+1}^\infty \xi_k$ , where, conditional on the stick-breaking process $\ell=(\ell_k)_{k\in\mathbb{N}}$ , the random variables $\{\xi_k^{{(\pm)}},\eta_k:k\in\mathbb{N}\}$ are independent and distributed as $\xi_k^{{(\pm)}}\sim F^{{(\pm)}}(\ell_k,\cdot)$ and $\eta_k\sim N\big(\ell_k b,\sigma^2 \ell_k\big)$ for $k\in\mathbb{N}$ .

The following result presents the decay of the strong error $\Delta_n^g$ for Lipschitz, locally Lipschitz, and two classes of barrier-type discontinuous payoffs that arise frequently in applications. Observe that, in all cases and under the corresponding mild assumptions, the pth moment of the strong error $\Delta_n^g$ decays exponentially fast in n with a rate $\vartheta>0$ that depends on the payoff g, the index $\beta_*$ defined in (21) below, and the degree p of the considered moment. In Proposition 3.1 and throughout the paper, the notation $f(n)=\mathcal{O}(g(n))$ as $n\to\infty$ for functions $f,g:\mathbb{N}\to(0,\infty)$ means $\limsup_{n\to\infty}f(n)/g(n)<\infty$ . Put differently, $f(n)=\mathcal{O}(g(n))$ is equivalent to f(n) being bounded above by $C_0 g(n)$ for some constant $C_0>0$ and all $n\in\mathbb{N}$ .

(Lipschitz.) Suppose $|g(x,y,t)-g(x,y',t')|\le K(|y-y'|+|t-t'|)$ for some K and all $(x,y,y',t,t')\in\mathbb{R}\times\mathbb{R}_+^2\times[0,T]^2$ . Then $\mathbb{E}_{{\boldsymbol0}}\big[\big|\Delta_n^g\big|^p\big]=\mathcal{O}\big(e^{-\vartheta_p n}\big)$ , where $\vartheta_p\in[\log(3/2),\log2]$ is in (23) below.

(Locally Lipschitz.) Let $|g(x,y,t)-g(x,y',t')|\le K(|y-y'|+|t-t'|)e^{\max\{y,y'\}}$ for some constant $K>0$ and all $(x,y,y',t,t')\in\mathbb{R}\times\mathbb{R}_+^2\times[0,T]^2$ . If $\lambda_+\geq q > 1$ and $\int_{[1,\infty)}e^{p(q-\lambda_+)x}\nu(\textrm{d} x)<\infty$ , then $\mathbb{E}_{{\boldsymbol0}}\big[\big|\Delta_n^g\big|^p\big]=\mathcal{O}\big(e^{-( \vartheta_{pr}/r)n}\big)$ , where $r\,:\!=\,(1-1/q)^{-1}>1$ and $\vartheta_{pr}\in[\log(3/2),\log2]$ is as in (23).

(Barrier type 1.) Suppose $g(\chi)=h(X_T)\cdot\mathbb{1}\{\overline{X}_T\le M\}$ for some $M>0$ and a measurable bounded function $h:\mathbb{R}\to\mathbb{R}$ . If $\mathbb{P}_{{\boldsymbol0}}(M<\overline{X}_T\le M+x)\le Kx$ for some $K>0$ and all $x\ge0$ , then for $\alpha_*\in (1,2]$ in (22) and any $\gamma\in(0,1)$ we have $\mathbb{E}_{{\boldsymbol0}}\big[\big|\Delta_n^g\big|^p\big]=\mathcal{O}\big(e^{-[\gamma\log(2)/(\gamma+\alpha_*)]n}\big)$ . Moreover, we may take $\gamma=1$ if any of the following hold: $\sigma^2>0$ or $\int_{(-1,1)}|x|\nu(\textrm{d} x)<\infty$ or Assumption 5.1 below.

(Barrier type 2.) Suppose $g(\chi)=h(X_T,\overline{X}_T)\cdot\mathbb{1}\{\tau_T\le s\}$ , where $s\in(0,T)$ , and h is measurable and bounded with $|h(x,y)-h(x,y')|\le K|y-y'|$ for some $K>0$ and all $(x,y,y')\in\mathbb{R}\times\mathbb{R}_+^2$ . If $\sigma^2>0$ or $\nu(\mathbb{R}\setminus\{0\})=\infty$ , then $\mathbb{E}_{{\boldsymbol0}}\big[\big|\Delta_n^g\big|^p\big]=\mathcal{O}\big(e^{-n/e}\big)$ .

3.2. Computational complexity and the central limit theorem for the Monte Carlo and multilevel Monte Carlo estimators

Consider the MC estimator

(7) \begin{equation}\hat\theta^{g,n}_{\textrm{MC}} \,:\!=\, \frac{1}{N}\sum_{i=1}^N \theta_i^{g,n},\end{equation}

where $\{\theta_i^{g,n}\}_{i\in\mathbb{N}}$ is the i.i.d. output of TSB-Alg with $\theta_1^{g,n}\overset{d}{=} g(\chi_n)\Upsilon_{\boldsymbol\lambda}$ (under $\mathbb{P}_{\boldsymbol0}$ ) and $n,N\in\mathbb{N}$ . The corresponding MLMC estimator is given by

(8) \begin{equation}\hat\theta^{g,n}_{\textrm{ML}}: = \sum_{k=1}^n\frac{1}{N_k}\sum_{i=1}^{N_k}D^g_{k,i},\end{equation}

where $\{D^g_{k,i}\}_{k,i\in\mathbb{N}}$ is an array of independent variables satisfying $D^g_{k,i}\overset{d}{=} (g(\chi_{k})-g(\chi_{k-1}))\Upsilon_{\boldsymbol\lambda}$ and $D^g_{1,i}\overset{d}{=} g(\chi_1)\Upsilon_{\boldsymbol\lambda}$ (under $\mathbb{P}_{\boldsymbol0}$ ), for $i\in\mathbb{N}$ , $k\ge 2$ , and $n,N_1,\ldots,N_n\in\mathbb{N}$ . Note that TSB-Alg can easily be adapted to sample the variable $D^g_{k,i}$ by drawing the ingredients for $(\chi_k,\Upsilon_{\boldsymbol\lambda})$ and computing $(\chi_{k-1},\chi_k,\Upsilon_{\boldsymbol\lambda})$ deterministically from the output; see [Reference González Cázares, Mijatović and Uribe Bravo18, Subsection 2.4] for further details. In the following, we refer to $\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]$ as the level variance of the MLMC estimator.

The computational complexity analysis of the MC and MLMC estimators is given in the next result (the usual notation $\lceil x\rceil\,:\!=\,\inf\{k\in\mathbb{N}:k\ge x\}$ , $x\in\mathbb{R}_+$ , is used for the ceiling function). In Proposition 3.2 and throughout the paper, the computational cost of an algorithm is measured as the total number of operations carried out by the algorithm. In particular, we assume that the following operations have computational costs uniformly bounded by some constant (measured, for instance, in units of time): simulation from the uniform law; simulation from the laws $F^{(\pm)}(t,\cdot)$ , $t>0$ ; evaluation of elementary mathematical operations such as addition, subtraction, multiplication, and division; and evaluation of elementary functions such as $\exp$ , $\log$ , $\sin$ , $\cos$ , $\tan$ , and $\arctan$ .

( MC.) Suppose $n=n(\varepsilon)$ and $N=\lceil2\varepsilon^{-2} \mathbb{V}_{\boldsymbol0}[g(\chi_n)\Upsilon_{\boldsymbol\lambda}]\rceil$ . Then the MC estimator $\hat\theta^{g,n}_{\textrm{MC}}$ is $\varepsilon$ -accurate, i.e. $\mathbb{E}_{\boldsymbol0}\big[\big|\hat\theta^{g,n}_{\textrm{MC}} -\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]\big|^2\big]\le\varepsilon^2$ , with expected cost $\mathcal{C}_{\textrm{MC}}(\varepsilon)\,:\!=\,c(n+1)N=\mathcal{O}\big(\varepsilon^{-2}\log(1/\varepsilon)\big)$ as $\varepsilon\to0$ .

(MLMC.) Suppose $n=n(\varepsilon)$ and set

(9) \begin{equation} N_k\,:\!=\,\Bigg\lceil2\varepsilon^{-2} \sqrt{\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]/k} \Bigg(\sum_{j=1}^{n}\sqrt{j\mathbb{V}_{\boldsymbol0}\big[D_{j,1}^g\big]}\Bigg) \Bigg\rceil, \quad k\in\{1,\ldots,n\}.\end{equation}

Then the MLMC estimator $\hat\theta^{g,n}_{\textrm{ML}}$ is $\varepsilon$ -accurate, and the corresponding expected cost equals

(10) \begin{equation}\mathcal{C}_{\textrm{ML}}(\varepsilon)\,:\!=\,2c\varepsilon^{-2} \Bigg(\sum_{k=1}^{n} \sqrt{k\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]}\Bigg)^2= \mathcal{O}\big(\varepsilon^{-2}\big)\quad \text{as $\varepsilon\to0$.}\end{equation}

Proposition 3.1 (with $p=1$ ) implies that the bias in (6) equals $\mathbb{E}_{\boldsymbol0}\big[\Delta_n^g\big]=\mathcal{O}\big(e^{-\vartheta_g n}\big)$ as $n\to\infty$ for some $\vartheta_g>0$ . Thus, the integer $n(\varepsilon)$ in Proposition 3.2 is finite for all payoffs g considered in Proposition 3.1, and moreover, $n(\varepsilon)=\mathcal{O}(\log(1/\varepsilon))$ as $\varepsilon\to0$ in all cases. In addition, by Remark 3.1(i) above, the variance of $\theta^{g,k}_1$ is bounded in $k\in\mathbb{N}$ :

$$\mathbb{V}_{\boldsymbol0}\big[\theta^{g,k}_1\big]\leq \mathbb{E}_{\boldsymbol0}\big[g(\chi_k)^2\Upsilon_{\boldsymbol\lambda}^2\big]\leq2\mathbb{E}_{\boldsymbol0}\big[g(\chi)^2\Upsilon_{\boldsymbol\lambda}^2\big]+2\mathbb{E}_{\boldsymbol0}\big[\big(\Delta_{k}^g\big)^2\big]\to2\mathbb{E}_{\boldsymbol0}\big[g(\chi)^2\Upsilon_{\boldsymbol\lambda}^2\big]<\infty\quad \text{as $k\to\infty$.}$$

Note that barrier-type payoffs g considered in Proposition 3.1 satisfy the second moment assumption, while in the Lipschitz or locally Lipschitz cases it is sufficient to assume additionally that $\lambda_+$ is either positive or strictly greater than one, respectively. Moreover, $\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]\le 2\mathbb{E}_{\boldsymbol0}\big[\big(\Delta_{k}^g\big)^2+\big(\Delta_{k-1}^g\big)^2\big]=\mathcal{O}\big(e^{-\vartheta_gk}\big)$ for $\vartheta_g>0$ bounded away from zero (see Remark 3.1(iii) above). These facts and the standard complexity analysis for MLMC (see e.g. [Reference González Cázares and Mijatović16, Appendix A] and the references therein) imply that the estimators $\hat\theta^{g,n}_{\textrm{MC}}$ and $\hat\theta^{g,n}_{\textrm{ML}}$ are $\varepsilon$ -accurate with the stated computational costs, proving Proposition 3.2.

We stress that the payoffs g in Proposition 3.2 include discontinuous payoffs in the supremum $\overline{X}_T$ (barrier type 1) or the time $\tau_T$ when this supremum is attained (barrier type 2), with the corresponding computational complexities of the MC and MLMC estimators given by $\mathcal{O}\big(\varepsilon^{-2}\log(1/\varepsilon)\big)$ and $\mathcal{O}\big(\varepsilon^{-2}\big)$ , respectively. This theoretical prediction matches the numerical performance of TSB-Alg for barrier options and the modified ulcer index; see Section 4.2 below.

For obtaining confidence intervals (CIs) for the MC and MLMC estimators $\hat\theta^{g,n}_{\textrm{MC}}$ and $\hat\theta^{g,n}_{\textrm{ML}}$ , a CLT can be very helpful. (The CIs derived in this paper do not account for model uncertainty or the uncertainty in the estimation or calibration of the parameters.) In fact, the CLT is necessary to construct a CI if the constants in the bounds on the bias in Proposition 3.1 are not explicitly known (e.g. for barrier-type-1 payoffs, the constant depends on the unknown value of the density of the supremum $\overline{X}_T$ at the barrier); see the discussion in [Reference González Cázares, Mijatović and Uribe Bravo18, Section 2.3]. Moreover, even if the bias can be controlled explicitly, the concentration inequalities typically lead to wider CIs than those based on the CLT; see [Reference Ben Alaya and Kebaier4, Reference Hoel and Krumscheid20]. The following result establishes the CLT for the MC and MLMC estimators valid for payoffs considered in Proposition 3.1.

Theorem 3.1. [CLT for TSB-Alg] Let g be any of the payoffs in Proposition 3.1, satisfying in addition $\mathbb{E}_{\boldsymbol0}\big[g(\chi)^2\Upsilon_{\boldsymbol\lambda}^2\big]<\infty$ . Let $\vartheta_g\in(0,\log2]$ be the rate satisfying $\mathbb{E}_{\boldsymbol0}\big[\big|\Delta_n^g\big|\big]=\mathcal{O}\big(e^{-\vartheta_gn}\big)$ , given in Proposition 3.1 and Remark 3.1(iii) above (with $p=1$ ). Fix $c_0>1/\vartheta_g$ , let $n=n(\varepsilon)\,:\!=\,\lceil c_0\log(1/\varepsilon)\rceil$ , and suppose N and $N_1,\ldots,N_n$ are as in Proposition 3.2. Then the MC and MLMC estimators $\hat\theta_{\textrm{MC}}^{g,n}$ and $\hat\theta_{\textrm{ML}}^{g,n}$ respectively satisfy the following CLTs (Z is a standard normal random variable):

(11) \begin{equation}\sqrt{2}\varepsilon^{-1}\big(\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}-\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]\big) \overset{d}{\to}Z,\qquad \text{and}\qquad \sqrt{2}\varepsilon^{-1}\big(\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}-\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]\big) \overset{d}{\to}Z,\quad \text{as }\varepsilon\to 0.\end{equation}

Theorem 3.1 works well in practice: in Figure 7 of Section 4.2 below we construct CIs (as a function of decreasing accuracy $\varepsilon$ ) for an MLMC estimator of a barrier option price under a tempered stable model. The rate $c_0$ can be taken arbitrarily close to $1/\vartheta_g$ , where $\vartheta_g$ is the corresponding geometric rate of decay of the error for the payoff g in Proposition 3.1 ( $\vartheta_g$ is bounded away from zero by Remark 3.1(iii) above), ensuring that the bias of the estimators vanishes in the limit.

By Lemma 5.2 below, the definition of the sample sizes N and $N_1,\ldots,N_n$ in Proposition 3.2 implies that the variances of the estimators $\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}$ and $\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}$ (under $\mathbb{P}_{\boldsymbol0}$ ) satisfy

\[\frac{\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}\big]}{\varepsilon^{2}/2}=\frac{\mathbb{V}_{\boldsymbol0}\big[\theta_1^{g,n(\varepsilon)}\big]}{\varepsilon^2 N/2}\to 1\quad \text{and} \quad \frac{\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}\big]}{\varepsilon^{2}/2}=\sum_{k=1}^{n(\varepsilon)}\frac{\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]}{\varepsilon^2 N_k/2}\to 1\quad \text{ as $\varepsilon\to0$}.\]

Hence, the CLT in (11) can also be expressed as follows:

\begin{align*}&\big(\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}-\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]\big)/\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}\big]^{1/2} \overset{d}{\to}Z \\\quad \text{and}\quad &\big(\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}-\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]\big)/\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}\big]^{1/2} \overset{d}{\to}Z,\quad \text{as }\varepsilon\to 0.\end{align*}

Since the variances $\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{MC}}^{g,n(\varepsilon)}\big]$ and $\mathbb{V}_{\boldsymbol0}\big[\hat\theta_{\textrm{ML}}^{g,n(\varepsilon)}\big]$ can be estimated from the sample, this is in fact how the CLT is often applied in practice. The proof of Theorem 3.1 is based on the CLT for triangular arrays and amounts to verifying Lindeberg’s condition for the estimators $\hat\theta_{\textrm{MC}}^{g,n}$ and $\hat\theta_{\textrm{ML}}^{g,n}$ ; see Section 5 below.

3.3. Unbiased estimators

It is known that when the MLMC complexity is optimal, it is possible to construct unbiased estimators without altering the optimal computational complexity $\mathcal{O}\big(\varepsilon^{-2}\big)$ as $\varepsilon\to 0$ . Such debiasing techniques are based around randomising the number of levels $\sup\{k\in\mathbb{N}:N_k>0\}$ and number of samples $(N_k)_{k\in\mathbb{N}}$ at each level of the variables $\{D^g_{k,i}\}_{k,i\in\mathbb{N}}$ in the MLMC estimator in (8); see e.g. [Reference McLeish30, Reference Rhee and Glynn34, Reference Vihola40]. More precisely, following [Reference Vihola40, Theorem 7], for any g in Proposition 3.1 and a parameter $N\in\mathbb{N}$ , we may construct an almost surely finite sequence of random integers $(N_k)_{k\in\mathbb{N}}$ (i.e. $\sup\{k\in\mathbb{N}:N_k>0\}<\infty$ a.s.) with explicit means $\mathbb{E}_{\boldsymbol0}[N_k]>0$ , such that the estimator

(12) \begin{equation}\hat{P}_N=:\sum_{k=1}^\infty\frac{1}{\mathbb{E}_{\boldsymbol0} [N_k]}\sum_{i=1}^{N_k}D_{k,i}^g\end{equation}

is unbiased, $\mathbb{E}_{\boldsymbol0}[\hat{P}_N]=\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]$ , and its variance $\mathbb{V}_{\boldsymbol0}[\hat{P}_N]$ and expected computational cost (under $\mathbb{P}_{\boldsymbol0}$ ) are proportional to $1/N$ and N, respectively, as $N\to\infty$ . The MC complexity analysis of the estimator $\hat{P}_N$ is then nearly identical to that of the classical MC estimator for exact simulation algorithms.

There are several parametric ways of constructing the random sequence $(N_k)_{k\in\mathbb{N}}$ (see e.g. [Reference Vihola40]), and it is also possible to optimise for the parameters appearing in these constructions as a function of the considered payoff g and other characteristics of X (see e.g. [Reference González Cázares, Mijatović and Uribe Bravo18, Section 2.5]). The details of such optimisations have been omitted in the present work in the interest of brevity, since they coincide with those found in [Reference González Cázares, Mijatović and Uribe Bravo18, Section 2.5 and Appendix A.3]. We comment that, as pointed out by the referee, in addition to achieving the optimal complexity, unbiased estimators (referred to as ‘single-term estimators’ in [Reference Rhee and Glynn34, Reference Vihola40]) also achieve the same asymptotic variance as the MLMC estimators (see [Reference Vihola40, Theorem 19]). Moreover, there are cases where unbiased estimators of a different form (referred to as ‘coupled-sum estimators’ in [Reference Rhee and Glynn34, Reference Vihola40]) achieve much smaller asymptotic variance than the MLMC (see Table 17 in the online supplement of [Reference Rhee and Glynn34]). We are grateful to the referee for this observation.

3.4. Comparisons

This subsection provides non-asymptotic and asymptotic performance comparisons of estimators based on TSB-Alg. The main aim is to develop rule-of-thumb principles guiding the user to the most effective estimator. In Subsection 3.4.1, for a given value of the accuracy $\varepsilon$ , we compare the computational complexity of the MC and MLMC estimators based on TSB-Alg. The MLMC estimator based on TSB-Alg is compared with the ones based on SB-Alg [Reference González Cázares, Mijatović and Uribe Bravo18] with rejection sampling (available only when the jumps of X are of finite variation) and SBG-Alg [Reference González Cázares and Mijatović16] in Subsections 3.4.2 and 3.4.3, respectively. In both cases we analyse the behaviour of the computational complexity in two regimes: (I) $\varepsilon$ tending to 0 and fixed time horizon T; (II) fixed $\varepsilon$ and time horizon T tending to 0 or $\infty$ .

Regime (II) is useful when there is a limited benefit to arbitrary accuracy in $\varepsilon$ but the constants may vary greatly with the time horizon T. For example, in option pricing, estimators with accuracy smaller than a basis point are of limited interest. For simplicity, in the remainder of this subsection the payoff g is assumed to be Lipschitz. However, analogous comparisons can be made for other payoffs under appropriate assumptions.

3.4.1. Comparison between the Monte Carlo and multilevel Monte Carlo estimators based on TSB-Alg

Recall first that both MC and MLMC estimators have the same bias, since the latter estimator is a telescoping sum of a sequence of the former estimators, controlled by $n(\varepsilon)$ given in Theorem 3.1 above.

Regime (I). Propositions 3.1 and 3.2 imply that the MLMC estimator outperforms the MC estimator as $\varepsilon\to 0$ . Moreover, since $\mathbb{V}_{\boldsymbol0}[g(\chi_n)\Upsilon_{\boldsymbol\lambda}]\to\mathbb{V}_{\boldsymbol0}[g(\chi)\Upsilon_{\boldsymbol\lambda}]$ and $\varepsilon^2\mathcal{C}_{\textrm{ML}}(\varepsilon)\to 2c\big(\sum_{k=1}^\infty(k\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big])^{1/2}\big)^2<\infty$ as $\varepsilon\to0$ , the MLMC estimator is preferable to the MC estimator for $\varepsilon > 0$ satisfying

\[n(\varepsilon)>\Bigg(\sum_{k=1}^\infty\sqrt{k\mathbb{V}_{\boldsymbol0}\big[D_{k,1}^g\big]}\Bigg)^2/\mathbb{V}_{\boldsymbol0}\big[g(\chi)\Upsilon_{\boldsymbol\lambda}\big].\]

Since the payoff g is Lipschitz, Proposition 3.1 implies that this condition is close to

\[\log(1/\varepsilon)>\vartheta_1\Bigg(\sum_{k=1}^\infty\sqrt{k\big(2^{-k}-e^{-2\vartheta_1k}\big)}\Bigg)^2,\]

where we recall that $\vartheta_1\in[\log(3/2),\log2]$ is defined in (23) below. In particular, the latter condition is equivalent to $\varepsilon<0.0915$ if $\vartheta_1=\log(3/2)$ , or $\varepsilon<5.06\times 10^{-5}$ if $\vartheta_1=\log2$ .

Regime (II). Assume $\varepsilon>0$ is fixed. In this case, the MC and MLMC estimators share the value of $n=n(\varepsilon)$ , which is $\mathcal{O}(\max\{1,\log T\})$ as either $T\to 0$ or $T\to\infty$ . Moreover, the variances $\mathbb{V}_{\boldsymbol0}[g(\chi_n)\Upsilon_{\boldsymbol\lambda}]$ (appearing in $\mathcal{C}_{\textrm{MC}}$ ) and $\mathbb{V}_{\boldsymbol0}\big[D^g_{k,1}\big]$ , $k\in\mathbb{N}$ (appearing in $\mathcal{C}_{\textrm{ML}}$ ; see Proposition 3.2 above) are all proportional to $\mathcal{O}\big(\big(T+T^2\big)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}\big)$ as either $T\to0$ or $T\to\infty$ . Therefore, by Proposition 3.2, the quotient $\mathcal{C}_{\textrm{MC}}/\mathcal{C}_{\textrm{ML}}$ is proportional to a constant as $T\to 0$ and a multiple of $\log T$ as $T\to\infty$ .

In conclusion, the MLMC estimator is preferable to the MC estimator when either $\varepsilon$ is small or T is large. Otherwise, when $\varepsilon$ is not small and T is small, the two estimators have similar performance.

3.4.2. Comparison with SB-Alg

In the special case when the jumps of X have finite variation (or equivalently, $\int_{(-1,1)}|x|\nu(\textrm{d} x)<\infty$ ), the increments $X_t$ can be simulated under $\mathbb{P}_{\boldsymbol\lambda}$ using rejection sampling (see [Reference Grabchak19, Reference Kawai and Masuda23]), making SB-Alg [Reference González Cázares, Mijatović and Uribe Bravo18] applicable to sample $\chi_n$ (see (5) for its definition) under $\mathbb{P}_{\boldsymbol\lambda}$ . The rejection sampling is performed for each of the increments of the subordinators

\[\widetilde Y^{{(\pm)}}_t \,:\!=\, \pm Y^{{(\pm)}}_t + d_\pm t,\quad \text{where}\quad d_+ \,:\!=\, \int_{(0,1)}x\nu(\textrm{d} x)\quad \text{and}\quad d_- \,:\!=\, \int_{({-}1,0)}|x|\nu(\textrm{d} x),\]

and the processes $Y^{{(\pm)}}$ are as in Theorem 2.1. The algorithm proposes samples under $\mathbb{P}_{\boldsymbol0}$ and rejects independently with probability $\exp\big({-}\lambda_\pm\widetilde Y_t^{{(\pm)}}\big)$ . Let ${\boldsymbol\lambda}_+ \,:\!=\, (\lambda_+,0)$ and ${\boldsymbol\lambda}_-\,:\!=\,(0,\lambda_-)$ ; then the expected number of proposals required for each sample equals $\exp\big(\gamma^{{(\pm)}}_{\boldsymbol\lambda} t\big)=1/\mathbb{E}_{\boldsymbol0}\big[\exp\big({-}\lambda_\pm\widetilde Y_t^{{(\pm)}}\big)\big]$ , where we define

(13) \begin{equation}\gamma^{{(\pm)}}_{\boldsymbol\lambda}\,:\!=\, \lambda_\pm d_\pm-\mu_{{\boldsymbol\lambda}_\pm}= \int_{\mathbb{R}_\pm}\big(1-e^{-\lambda_{\pm}|x|}\big)\nu(\textrm{d} x) \in [0,\infty).\end{equation}

(Note that $\mu_{p{\boldsymbol\lambda}}-p \mu_{\boldsymbol\lambda} = p \big(\gamma^{{(+)}}_{{\boldsymbol\lambda}}+\gamma^{{(-)}}_{{\boldsymbol\lambda}}\big) - \big(\gamma^{{(+)}}_{p{\boldsymbol\lambda}}+\gamma^{{(-)}}_{p{\boldsymbol\lambda}}\big)$ ; see (4).)

We need the following elementary lemma to analyse the computational complexity of SB-Alg with rejection sampling.

Lemma 3.1. (a) Let $\ell$ be a stick-breaking process on [0,Reference Andersen and Lipton1]; then for any $n\in\mathbb{N}$ we have

(14) \begin{equation} 0\le n + \int_0^1\frac{1}{x}\big(e^{cx}-1\big)\textrm{d} x - \sum_{k=1}^n \mathbb{E}\big[e^{c\ell_k}\big] \le 2^{-n}\int_0^1\frac{1}{x}\big(e^{cx}-1\big)\textrm{d} x.\end{equation}

(b) We have $c^{-1}e^{-c}\big(1+c^2\big)\int_0^1x^{-1}\big(e^{cx}-1\big)\textrm{d} x\to1$ as either $c\to 0$ or $c\to\infty$ .

Assume that the simulation of the increments $Y^{{(\pm)}}_t$ under $\mathbb{P}_{\boldsymbol0}$ has constant cost independent of the time horizon t (we also assume without loss of generality that the simulation of uniform random variables and the evaluation of operators such as sums and products and of the exponential functions have constant cost). Since SB-Alg requires the rejection sampling to be carried out over random stick-breaking lengths, the expected cost of SB-Alg with rejection sampling is, by (14) in Lemma 3.1, asymptotically proportional to

(15) \begin{equation}\sum_{k=1}^n \mathbb{E}\Big[e^{\gamma^{{(+)}}_{\boldsymbol\lambda}\ell_k} + e^{\gamma^{{(-)}}_{\boldsymbol\lambda}\ell_k}\Big]= 2n + \big(1+\mathcal{O}\big(2^{-n}\big)\big)\int_0^1\frac{1}{x} \Big(e^{\gamma_{\boldsymbol\lambda}^{{(+)}} Tx} + e^{\gamma_{\boldsymbol\lambda}^{{(-)}} Tx} - 2\Big)\textrm{d} x,\quad \text{as }n\to\infty.\end{equation}

In fact, by Lemma 3.1(a), the absolute value of the term $\mathcal{O}\big(2^{-n}\big)$ is bounded by $2^{-n}$ . Moreover, by Lemma 3.1(b), the integral in (15) may be replaced with an explicit expression

(16) \begin{equation}\Gamma_{{\boldsymbol\lambda}}(T)\,:\!=\,\frac{\gamma_{\boldsymbol\lambda}^{{(+)}} T}{1+\big(\gamma_{\boldsymbol\lambda}^{{(+)}} T\big)^2} e^{\gamma_{\boldsymbol\lambda}^{{(+)}} T}+\frac{\gamma_{\boldsymbol\lambda}^{{(-)}} T}{1+\big(\gamma_{\boldsymbol\lambda}^{{(-)}} T\big)^2} e^{\gamma_{\boldsymbol\lambda}^{{(-)}} T},\end{equation}

as T tends to either 0 or $\infty$ .

Table 2 shows how SB-Alg with rejection sampling compares to TSB-Alg above.

Table 2. Level variance and expected cost of SB-Alg and TSB-Alg, up to multiplicative constants that do not depend on the time horizon T. The bounds on the level variances (defined in the first paragraph of Subsection 3.2 above) follow from [Reference González Cázares, Mijatović and Uribe Bravo18, Theorem 2] for SB-Alg and Proposition 3.1 for TSB-Alg.

Regime (I). By Table 2, we can deduce that the MC and MLMC estimators of both algorithms have the complexities $\mathcal{O}\big(\varepsilon^{-2}\log(1/\varepsilon)\big)$ and $\mathcal{O}\big(\varepsilon^{-2}\big)$ as $\varepsilon\to0$ , respectively, for all the payoffs considered in Proposition 3.1.

Regime (II). Assume $\varepsilon$ is fixed. The biases of the two algorithms agree and equal $\mathcal{O}\big(e^{-\vartheta_1n}\big(\sqrt{T}+T\big)\big)$ , implying $n=\log\big(\big(\sqrt{T}+T\big)/\varepsilon\big)/\vartheta_1+\mathcal{O}(1)$ (with $\vartheta_1$ defined in (23)). The level variances of SB-Alg and TSB-Alg are $\mathcal{O}\big(e^{-\log(2)n}\big(T+T^2\big)\big)$ and $\mathcal{O}\big(e^{-\log(2)n}\big(T+T^2\big)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}\big)$ , with costs $\mathcal{O}(n+\Gamma_{{\boldsymbol\lambda}}(T))$ and $\mathcal{O}(n)$ , respectively. Thus, by (16), the ratios of the level variance and cost converge to 1 as $T\to 0$ , so the ratio of the complexities of both algorithms converges to 1, implying that one should use TSB-Alg, as its performance in this regime is no worse than that of the other algorithm, but its implementation is simpler. For moderate or large values of T, by [Reference González Cázares and Mijatović16, Equation (A.3)], the computational complexity of the MLMC estimator based on TSB-Alg is proportional to $\varepsilon^{-2}e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}(T+T^2)$ and that of SB-Alg is proportional to $\varepsilon^{-2}(1+\Gamma_{{\boldsymbol\lambda}}(T))(T+T^2)$ , where, in both cases, the proportionality constants do not depend on the time horizon T. Since both constants are exponential in T, for large T we need only compare $\max\big\{\gamma_{\boldsymbol\lambda}^{{(+)}},\gamma_{\boldsymbol\lambda}^{{(-)}}\big\}$ with $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ . Indeed, TSB-Alg has a smaller complexity than SB-Alg for all sufficiently large T if and only if $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}<\max\big\{\gamma_{\boldsymbol\lambda}^{{(+)}},\gamma_{\boldsymbol\lambda}^{{(-)}}\big\}$ . In Subsection 4.3.1 below, we provide an explicit criterion for the tempered stable process in terms of the parameters; see Figure 8.

In conclusion, when X has jumps of finite variation, it is preferable to use TSB-Alg if T is small or $e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}<1+\Gamma_{\boldsymbol\lambda}(T)$ . Moreover, this typically holds if the Blumenthal–Getoor index of X is larger than $\log_2(3/2)<0.6$ ; see Subsection 4.3.1 and Figure 8 below for details.

3.4.3. Comparison with SBG-Alg

Given any cutoff level $\kappa\in(0,1]$ , the algorithm SBG-Alg approximates the Lévy process X (under $\mathbb{P}_{\boldsymbol\lambda}$ with the generating triplet $\big(\sigma^2, \nu_{\boldsymbol\lambda}, b_{\boldsymbol\lambda}\big)$ ) with a jump-diffusion $X^{(\kappa)}$ and samples exactly the vector $\big(X^{(\kappa)}_T,\overline{X}^{(\kappa)}_T,\tau_T\big(X^{(\kappa)}\big)\big)$ ; see [Reference González Cázares and Mijatović16] for details. The bias of SBG-Alg is $\mathcal{O}\big(\min\{\sqrt{T}\kappa^{1-\beta_*/2},\kappa\} \big(1+\log^+\big(T\kappa^{-\beta_*}\big)\big)\big)$ with corresponding expected cost $\mathcal{O}\big(T\kappa^{-\beta_*}\big)$ . Thus, when parametrised by computational effort n, the bias of SBG-Alg is $\mathcal{O}\big(T^{1/\beta_*}n^{-1/\beta_*}\log n\big)$ while the bias of TSB-Alg is $\mathcal{O}\big(e^{-\vartheta_1 n}\big(\sqrt{T}+T\big)\big)$ .

Regime (I). The complexities of the MC and MLMC estimators based on SBG-Alg are, by [Reference González Cázares and Mijatović16, Theorem 22], equal to $\mathcal{O}\big(\varepsilon^{-2-\beta_*}\big)$ and $\mathcal{O}\big(\varepsilon^{-\max\{2,2\beta_*\}}\big(1+\mathbb{1}_{\{\beta_*=1\}}\log(\varepsilon)^2\big)\big)$ , respectively, where $\beta_*$ , defined in (21) below, is slightly larger than the Blumenthal–Getoor index of X. Since, by Subsection 3.2 above, the complexities of the MC and MLMC estimators based on TSB-Alg are $\mathcal{O}\big(\varepsilon^{-2}\log(1/\varepsilon)\big)$ and $\mathcal{O}\big(\varepsilon^{-2}\big)$ , respectively, the complexity of TSB-Alg is always dominated by that of SBG-Alg as $\varepsilon\to0$ .

Regime (II). Suppose $\varepsilon>0$ is fixed. Then, as in Subsection 3.4.2 above, the computational complexity of the MLMC estimator based on TSB-Alg is $\mathcal{O}\big(\varepsilon^{-2}(T+T^2)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}\big)$ , where the multiplicative constant does not depend on the time horizon T. By [Reference González Cázares and Mijatović16, Theorems 3 and 22] and [Reference González Cázares and Mijatović16, Equation (A.3)], the complexity of the MLMC estimator based on SBG-Alg is $\mathcal{O}\big(\varepsilon^{-2}\big(C_1 T+C_2T^2\big)\big)$ if $\beta_*< 1$ , $\mathcal{O}\big(\varepsilon^{-2}\big(C_1T+C_2T^2\log^2(1/\varepsilon)\big)\big)$ if $\beta_*= 1$ , and $\mathcal{O}\big(\varepsilon^{-2}\big(C_1T+C_2T^2\varepsilon^{-2(\beta_*-1)}\big)\big)$ otherwise, where

\begin{align*}C_1 &\,:\!=\, e^{r\beta_*}/\big(1-e^{r(\beta_*/2-1)}\big)^{2},\\C_2 &\,:\!=\, e^{r\beta_*}/\big(1-e^{r(\beta_*-1)}\big)^{2}\cdot\mathbb{1}_{\{\beta_*\ne 1\}}+(e/2)^2\cdot\mathbb{1}_{\{\beta_*= 1\}}, \text{ and}\\r &\,:\!=\, (2/|\beta_*-1|)\log(1+|\beta_*-1|/\beta_*)\cdot\mathbb{1}_{\{\beta_*\ne 1\}}+ 2\cdot\mathbb{1}_{\{\beta_*=1\}}. \end{align*}

None of the other multiplicative constants in these bounds depend on the time horizon T. Thus TSB-Alg outperforms SBG-Alg if and only if we are in one of the following cases:

• • $\beta_*< 1$ and $(1+T)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T} <C_1+C_2T$ ,

• • $\beta_*=1$ and $(1+T)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T} <C_1+C_2T\log^2(1/\varepsilon)$ ,

• • $\beta_*>1$ and $(1+T)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T} <C_1+C_2T\varepsilon^{-2(\beta_*-1)}$ .

Note that the constant $C_2$ is unbounded for $\beta_*$ close to 1, favouring TSB-Alg.

In conclusion, TSB-Alg is simpler than SBG-Alg [Reference González Cázares and Mijatović16] and outperforms it asymptotically as $\varepsilon\to0$ . Moreover, TSB-Alg performs better than SBG-Alg for a fixed accuracy $\varepsilon>0$ if either (I) $\beta_*<1$ and the time horizon $T\ll 1$ satisfies the inequality $T<\log(C_1)/(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})$ or (II) $\beta_*\geq 1$ and T is not large. In Subsection 4.3.2 below, we apply this criterion to tempered stable processes; see Figure 9 for the case (I) with $\beta_*< 1$ and Figure 10 for the case (II) with $\beta_*\ge 1$ .

3.5. Variance reduction via exponential martingales

It follows from Subsection 3.4 that the performance of TSB-Alg deteriorates if the expectation $\mathbb{E}_{\boldsymbol0} \big[\Upsilon_{\boldsymbol\lambda}^2\big] = \exp((\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T)$ is large, making the variance of the estimator large. This problem is mitigated by using the control variates method, a variance reduction technique, based on exponential $\mathbb{P}_{\boldsymbol0}$ -martingales, at (almost) no additional computational cost.

Suppose we are interested in estimating $\mathbb{E}_{\boldsymbol\lambda}[\zeta]=\mathbb{E}_{\boldsymbol0}[\zeta\Upsilon_{\boldsymbol\lambda}]$ , where $\zeta$ is either $g(\chi_n)$ (MC) or $g(\chi_n)-g(\chi_{n-1})$ (MLMC). We propose to draw samples of $\tilde\zeta$ under $\mathbb{P}_{\boldsymbol0}$ , instead of $\zeta\Upsilon_{\boldsymbol\lambda}$ , where

\[\tilde\zeta\,:\!=\,\zeta\Upsilon_{\boldsymbol\lambda}- w_0 (\Upsilon_{\boldsymbol\lambda}-1)- w_+ (\Upsilon_{{\boldsymbol\lambda}_+}-1)- w_- (\Upsilon_{{\boldsymbol\lambda}_-}-1),\]

and $w_0,w_+,w_-\in\mathbb{R}$ are constants to be determined (recall ${\boldsymbol\lambda}_+=(\lambda_+,0)$ and ${\boldsymbol\lambda}_-=(0,\lambda_-)$ ). Clearly $\mathbb{E}_{\boldsymbol0}[\tilde\zeta]=\mathbb{E}_{\boldsymbol0}[\zeta\Upsilon_{\boldsymbol\lambda}]$ since the variables $\Upsilon_{\boldsymbol\lambda}$ , $\Upsilon_{{\boldsymbol\lambda}_+}$ , and $\Upsilon_{{\boldsymbol\lambda}_-}$ have unit mean under $\mathbb{P}_{\boldsymbol0}$ . We choose $w_0$ , $w_+$ , and $w_-$ to minimise the variance of $\tilde\zeta$ , by minimising the quadratic form [Reference Glasserman15, Section 4.1.2] of $\mathbb{V}_{\boldsymbol0}[\tilde\zeta]$ :

\[\left(\begin{array}{c}-1\\\omega_0\\\omega_+\\\omega_-\end{array}\right)^{\top}\left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\mathbb{V}_{\boldsymbol0}[\zeta\Upsilon_{\boldsymbol\lambda}]& \textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{\boldsymbol\lambda})& \textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& \textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})\\\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{\boldsymbol\lambda})& \mathbb{V}_{\boldsymbol0}[\Upsilon_{\boldsymbol\lambda}]& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})\\\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& \mathbb{V}_{\boldsymbol0}[\Upsilon_{{\boldsymbol\lambda}_+}]& \textrm{cov}_{\boldsymbol0}(\Upsilon_{{\boldsymbol\lambda}_+},\Upsilon_{{\boldsymbol\lambda}_-})\\\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{{\boldsymbol\lambda}_+},\Upsilon_{{\boldsymbol\lambda}_-})& \mathbb{V}_{\boldsymbol0}[\Upsilon_{{\boldsymbol\lambda}_-}]\end{array}\right)\left(\begin{array}{c}-1\\\omega_0\\\omega_+\\\omega_-\end{array}\right).\]

The solution, in terms of the covariance matrix (under $\mathbb{P}_{\boldsymbol0}$ ) of the vector $(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+},\Upsilon_{{\boldsymbol\lambda}_-})$ , is

\[\left(\begin{array}{c}w_0\\w_+\\w_-\end{array}\right)=\left(\begin{array}{c@{\quad}c@{\quad}c}2\mathbb{V}_{\boldsymbol0}[\Upsilon_{\boldsymbol\lambda}]& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})\\\textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})& 2\mathbb{V}_{\boldsymbol0}[\Upsilon_{{\boldsymbol\lambda}_+}]& \textrm{cov}_{\boldsymbol0}(\Upsilon_{{\boldsymbol\lambda}_+},\Upsilon_{{\boldsymbol\lambda}_-})\\\textrm{cov}_{\boldsymbol0}(\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})& \textrm{cov}_{\boldsymbol0}(\Upsilon_{{\boldsymbol\lambda}_+},\Upsilon_{{\boldsymbol\lambda}_-})& 2\mathbb{V}_{\boldsymbol0}[\Upsilon_{{\boldsymbol\lambda}_-}]\end{array}\right)^{-1}\left(\begin{array}{c}\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{\boldsymbol\lambda})\\\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_+})\\\textrm{cov}_{\boldsymbol0}(\zeta\Upsilon_{\boldsymbol\lambda},\Upsilon_{{\boldsymbol\lambda}_-})\end{array}\right).\]

In practice, these covariances are estimated from the same samples that were drawn to estimate $\mathbb{E}_{\boldsymbol0}[\zeta\Upsilon_{\boldsymbol\lambda}]$ . The additional cost is (nearly) negligible as all the variables in the exponential martingales are by-products of TSB-Alg. It is difficult to establish theoretical guarantees for the improvement of $\tilde \zeta$ over $\zeta\Upsilon_{\boldsymbol\lambda}$ . However, since most of the variance of the estimator based on $\zeta\Upsilon_{\boldsymbol\lambda}$ comes from $\Upsilon_{\boldsymbol\lambda}$ , the proposal $\tilde \zeta$ typically performs well in applications; see e.g. the CIs in Figures 3 and 4 for a tempered stable process.

4.1. Description of the model

In the present section we apply our results to the tempered stable process X. More precisely, given ${\boldsymbol\lambda}\in\mathbb{R}_+^2\setminus\{{\boldsymbol0}\}$ , the tempered stable Lévy measure $\nu_{\boldsymbol\lambda}$ specifies the Lévy measure $\nu$ via (1):

(17) \begin{equation}\frac{\nu(\textrm{d} x)}{\textrm{d} x}=\frac{c_+}{x^{\alpha_++1}} \cdot\mathbb{1}_{(0,\infty)}(x) + \frac{c_-}{|x|^{\alpha_-+1}} \cdot\mathbb{1}_{(-\infty,0)}(x),\end{equation}

where $c_\pm\ge 0$ and $\alpha_\pm\in(0,2)$ . The drift b is given by the tempered stable drift $b_{\boldsymbol\lambda}\in\mathbb{R}$ via (1), and the constant $\mu_{\boldsymbol\lambda}$ is given in (4). Both b and $\mu_{\boldsymbol\lambda}$ can be computed using (18) and (19) below. TSB-Alg samples from the distribution functions $F^{{(\pm)}}(t,x)=\mathbb{P}_{\boldsymbol0}\big(Y_t^{{(\pm)}}\le x\big)$ , where $Y^{{(\pm)}}$ are the spectrally one-sided stable processes defined in Theorem 2.1, using the Chambers–Mallows–Stuck algorithm [Reference Chambers, Mallows and Stuck9]. In the appendix we include a version of this algorithm, given explicitly in terms of the drift b and the parameters in the Lévy measure $\nu$ ; see Algorithm 2 below.

Next, we provide explicit formulae for b and $\mu_{\boldsymbol\lambda}$ in terms of special functions. We begin by expressing b in terms of $b_{\boldsymbol\lambda}$ (see (1) above):

(18) \begin{equation}b=b_{\boldsymbol\lambda}-c_+B_{\alpha_+,\lambda_+}+c_-B_{\alpha_-,\lambda_-},\quad \text{where}\quad B_{a,r} \,:\!=\, \int_0^1 \big(e^{-rx}-1\big)x^{-a}\textrm{d} x.\end{equation}

We have $B_{a,0}=0$ for any $a\ge0$ and, for $r>0$ ,

\[B_{a,r}=\sum_{n=1}^\infty \frac{({-}r)^n}{n!(n-a-1)}=\begin{cases}\big(e^{-r}-1+r^{a-1}\gamma(2-a,r)\big)/(1-a),&a\in(0,2)\setminus\{1\},\\-\gamma-\log r-E_1(r),&a=1,\end{cases}\]

where $\gamma(a,r)=\int_0^r e^{-x}x^{a-1}\textrm{d} x=\sum_{n=0}^\infty ({-}1)^{n}r^{n+a}/(n!(n+a))$ , $a>0$ , is the lower incomplete gamma function, $E_1(r)=\int_r^\infty e^{-x}x^{-1}\textrm{d} x$ , $r>0$ , is the exponential integral, and $\gamma=0.57721566\ldots$ is the Euler–Mascheroni constant. Similarly, to compute $\mu_{\boldsymbol\lambda}$ from (4), note that

(19) \begin{equation}\mu_{\boldsymbol\lambda} =c_+C_{\alpha_+,\lambda_+}+c_-C_{\alpha_-,\lambda_-},\quad \text{where}\quad C_{a,r}\,:\!=\, \int_0^\infty \big(e^{-rx}-1+rx\cdot \mathbb{1}_{(0,1)}(x)\big)x^{-a-1}\textrm{d} x.\end{equation}

Clearly, $C_{a,0}=0$ for any $a\ge0$ and, for $r>0$ ,

\begin{align*}C_{a,r}&=-\frac{1}{a}\big(e^{-r}+r(1+B_{a,r})\big)+r^{a}\Gamma({-}a,r),\end{align*}

where $\Gamma(a,r)=\int_r^\infty e^{-x}x^{a-1}dx$ is the upper incomplete gamma function. When $a<1$ , this simplifies to $C_{a,r}=r^a \Gamma({-}a) + r/(1-a)$ , where $\Gamma(\cdot)$ is the gamma function.

As discussed in Subsection 3.4 (see also Table 2 above), the performance of TSB-Alg deteriorates for large values of the constant $\mu_{2{\boldsymbol\lambda}}-2\mu_{{\boldsymbol\lambda}}$ . As a consequence of the formulae above, it is easy to see that this constant is proportional to $\lambda_+^{\alpha_+}/(2-\alpha_+) + \lambda_-^{\alpha_-}/(2-\alpha_-)$ as either $\lambda_\pm\to\infty$ or $\alpha_\pm\to 2$ ; see Figure 1.

Figure 1. The picture shows the map $l\mapsto \mu_{\boldsymbol\lambda}$ , where $c_+=c_-=1$ , $\alpha_+=\alpha_-=\alpha$ , ${\boldsymbol\lambda}=(l,l)$ , and $\alpha\in\{0.2,\ldots,1.8\}$ .

It is natural to expect the performance of TSB-Alg to deteriorate as $\lambda_\pm\to\infty$ or $\alpha_\pm\to 2$ , because the variance of the Radon–Nikodym derivative $\Upsilon_{\boldsymbol\lambda}$ tends to $\infty$ , making the variance of all the estimators very large.

4.2. Monte Carlo and multilevel Monte Carlo estimators for tempered stable processes

Let X denote the tempered stable process with generating triplet $\big(\sigma^2,\nu_{\boldsymbol\lambda},b_{\boldsymbol\lambda}\big)$ given in Subsection 4.1. Then $S_t=S_0e^{X_t}$ models the risky asset for some initial value $S_0>0$ and all $t\ge 0$ . Observe that $\overline{S}_t=S_0 e^{\overline{X}_t}$ and that the time at which S attains its supremum on [0,t] is also $\tau_t$ . Since the tails $\nu(\mathbb{R}\setminus({-}x,x))$ of $\nu$ decay polynomially as $x\to\infty$ (see (17)), $S_t$ and $\overline{S}_t$ satisfy the following moment conditions under $\mathbb{P}_{\boldsymbol\lambda}$ : $\mathbb{E}_{\boldsymbol\lambda}\big[S_t^r\big]<\infty$ (resp. $\mathbb{E}_{\boldsymbol\lambda}\big[\overline{S}_t^r\big]<\infty$ ) if and only if $r\in[-\lambda_-,\lambda_+]$ (resp. $r\le \lambda_+$ ). In the present subsection we apply TSB-Alg to estimate $\mathbb{E}_{\boldsymbol\lambda} [g(\chi)]$ using the MC and MLMC estimators in (7) and (8), respectively, for payoffs g that arise in applications and calibrated/estimated values of tempered stable model parameters.

Figure 2. Solid lines indicate the estimated value of $\mathbb{E}_{\boldsymbol\lambda} [g(\chi_n)]$ for the USD/JPY currency pair and the function $g(\chi)= \max\{S_T-K,0\}\mathbb{1}{\{\overline{S}_T \le M\}}$ using $N=10^6$ samples with $T\in\big\{\frac{30}{365},\frac{90}{365}\big\}$ , $K=95$ , and $M=102$ . Dotted lines (the symmetric bands around each solid line) indicate the CI of the estimates with confidence level 95%. These CIs are very tight, making them hard to discern in the plot.

Figure 3. Solid lines indicate the estimated value of the ulcer index (UI) $100\sqrt{\mathbb{E}_{\boldsymbol\lambda} [g(\chi_n)]}$ using $N=10^7$ samples, where $g(\chi)= (S_T/\overline{S}_T-1)^2$ . Dotted lines (the bands around each solid line) indicate the 95% CI of the estimates.

4.2.1. Monte Carlo estimators

Consider the estimator in (7) for the following payoff functions g: (I) the payoff of the up-and-out barrier call option $g(\chi)=\max\{S_T-K,0\}\mathbb{1}{\{\overline{S}_T\le M\}}$ and (II) the function $g(\chi)=(S_T/\overline{S}_T-1)^2$ associated to the ulcer index (UI) given by $100\sqrt{\mathbb{E}_{\boldsymbol\lambda} g(\chi)}$ (see [Reference Downey11]). We plot the estimated value along with a band representing the empirically constructed CIs of level 95% via Theorem 3.1; see Figures 2 and 3. We set the initial value $S_0=100$ , assume $b_{\boldsymbol\lambda}=0$ , and consider multiple time horizons T (given as fractions of a calendar year): we take $T\in\big\{\frac{30}{365},\frac{90}{365}\big\}$ in (I) and $T\in\big\{\frac{14}{365},\frac{28}{365}\big\}$ in (II). These time horizons are common in their respective applications, (I) option pricing and (II) real-world risk assessment (see e.g. [Reference Andersen and Lipton1, Reference Carr, Geman, Madan and Yor7]). The values of the other parameters, given in Table 3 below, are either calibrated or estimated: in (I) we are pricing an option and thus use calibrated parameters with respect to the risk-neutral measure obtained in [Reference Andersen and Lipton1, Table 3] for the USD/JPY foreign exchange (FX) rate, and in (II) we are assessing risks under the real-world probability measure and hence use the maximum likelihood estimates in [Reference Carr, Geman, Madan and Yor7, Table 2] for various time series of historic asset returns. In both (I) and (II), it is assumed that $\alpha_\pm=\alpha$ ; (II) additionally assumes that $c_+=c_-$ (i.e., X is a CGMY process). The stocks MCD, BIX, and SOX in (II) were chosen with $\alpha>1$ to stress-test the performance of TSB-Alg when $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ is big, forcing the variance of the estimator to be large; see the discussion at the end of Subsection 4.1 above and Figure 1.

Table 3. Calibrated/estimated parameters of the tempered stable model. The first two sets of parameters were calibrated in [Reference Andersen and Lipton1, Table 3] based on FX option data. The bottom three sets of parameters were estimated using the maximum likelihood estimators in [Reference Carr, Geman, Madan and Yor7, Table 2], based on a time series of equity returns.

The visible increase in the variance of the estimates for the SOX (see Figure 3) as the maturity increases from 14 to 28 days is a consequence of the unusually large value of $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ ; see Table 3 above. If $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ is not large, we may take long time horizons and relatively few samples to obtain an accurate MC estimate. In fact this appears to be the case for calibrated risk-neutral parameters (we were unable to find calibrated parameter values in the literature with a large value of $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ ). However, if $\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ is large, then for any moderately large time horizon, an accurate MC estimate would require a large number of samples. In such cases, the control variates method from Subsection 3.5 becomes very useful.

To illustrate the added value of the control variates method from Subsection 3.5, we apply it in the setting of Figure 3, where we observed the widest CIs. (Recall that the CIs derived in this paper do not account for model uncertainty or the uncertainty in the estimation or calibration of the parameters.) Figure 4 displays the resulting estimators and CIs, showing that this method is beneficial in the case where the variance of the Radon–Nikodym derivative $\Upsilon_{\boldsymbol\lambda}$ is large, i.e., when $(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T$ is large. The CIs for the SOX asset at $n=12$ shrank by a factor of $4.23$ ; in other words, the variance became $5.58\%$ of its original value.

Figure 4. Solid lines indicate the estimated value of the ulcer index (UI) $100\sqrt{\mathbb{E}_{\boldsymbol\lambda} g(\chi_n)}$ using $N=10^7$ samples and the control variates method of Subsection 3.5, where $g(\chi)= (S_T/\overline{S}_T-1)^2$ . Dotted lines (the bands around each solid line) indicate the CI of the estimates with confidence level 95%.

Figure 5. The top two panels show the bias and level variance decay as a function of k for the payoff $g(\chi)=\max\{S_T-K,0\}\mathbb{1}{\{\overline{S}_T \le M\}}$ with $S_0=100$ , $T=90/365$ , $K=95$ , $M=102$ and the parameter values for the USD/JPY FX rate (v2) in Table 3. The theoretical predictions (dashed) are based on Proposition 3.1 for barrier-type-1 payoffs. The bottom two panels show the corresponding values of the complexities and parameters $n,N_1,\ldots,N_n$ associated to the precision levels $\varepsilon\in\{2^{-9},2^{-10},\ldots,2^{-19}\}$ .

4.2.2. Multilevel Monte Carlo estimators

We will consider the MLMC estimator in (8) with parameters $n,N_1,\ldots,N_n\in\mathbb{N}$ given by [Reference González Cázares and Mijatović16, Equations (A.1)–(A.2)]. In this example we chose (I) the payoff $g(\chi)=\max\{S_T-95,0\}\mathbb{1}{\{\overline{S}_T \le 102\}}$ from Subsection 4.2.1 above and (II) the payoff $g(\chi)=(S_T/\overline{S}_T-1)^2\mathbb{1}_{\{\tau_T<T/2\}}$ , associated to the modified ulcer index (MUI) $100\sqrt{\mathbb{E}_{\boldsymbol\lambda}[g(\chi)]}$ , a risk measure which weighs trends more heavily than short-time fluctuations (see [Reference González Cázares and Mijatović16] and the references therein).

The payoff in (I) is that of a barrier up-and-out option, so it is natural to use the risk-neutral parameters for the USD/JPY FX rate (see (v2) in Table 3) over the time horizon $T=90/365$ . For (II) we take the parameter values of the MCD stock in Table 3 with $T=28/365$ . In both cases we set $S_0=100$ . Figure 5 (resp. Figure 6) shows the decay of the bias and level variance, the corresponding value of the constants $n,N_1,\ldots,N_n$ , and the growth of the complexity for the first (resp. second) payoff.

Figure 6. The top two panels show the bias and level variance decay as a function of k for the payoff $g(\chi)=(S_T/\overline{S}_T-1)^2\mathbb{1}_{\{\tau_T<T/2\}}$ with $S_0=100$ , $T=28/365$ and the parameter values for MCD in Table 3. The theoretical predictions (dashed) are based on Proposition 3.1 for barrier-type-2 payoffs. The bottom two panels show the corresponding values of the complexities and parameters $n,N_1,\ldots,N_n$ associated to the precision levels $\varepsilon\in\big\{2^{-9},2^{-10},\ldots,2^{-13}\big\}$ .

In Figure 7 we plot the estimator $\hat\theta^{g,n}_{\textrm{MC}}$ in Theorem 3.1 for (I), parametrised by $\varepsilon\to0$ . To further illustrate the CLT in Theorem 3.1, the figure also shows the CIs of confidence level 95% constructed using the CLT.

Figure 7. The solid line indicates the estimated value of the expectation $\mathbb{E}_{\boldsymbol\lambda} g(\chi_n)$ for the payoff $g(\chi)=\max\{S_T-K,0\}\mathbb{1}{\{\overline{S}_T \le M\}}$ with $S_0=100$ , $T=90/365$ , $K=95$ , $M=102$ and the parameter values for the USD/JPY FX rate (v2) in Table 3. We use the MLMC estimator $\hat\theta^{g,n}_{\textrm{ML}}$ in (8), and the CIs (dotted lines) are constructed using Theorem 3.1 with confidence level 95%.

4.3. Comparing TSB-Alg with existing algorithms for tempered stable processes

In this subsection we take the analysis from Subsection 3.4 and apply it to the tempered stable case.

4.3.1. Comparison with SB-Alg

Recall from Subsection 3.4.2 that SB-Alg is only applicable when $\alpha_\pm<1$ , and under Regime (II), SB-Alg is preferable over TSB-Alg for all sufficiently large T if and only if $\max\big\{\gamma_{\boldsymbol\lambda}^{{(+)}},\gamma_{\boldsymbol\lambda}^{{(-)}}\big\}\le\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ , where $\gamma_{\boldsymbol\lambda}^{{(\pm)}} = -c_\pm \lambda_\pm^{\alpha_\pm} \Gamma({-}\alpha_\pm)\ge 0$ is defined in (13). By the formulae in Subsection 4.1, it is easily seen that $\max\big\{\gamma_{\boldsymbol\lambda}^{{(+)}},\gamma_{\boldsymbol\lambda}^{{(-)}}\big\}\le\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda}$ is equivalent to

\begin{align*}&\min\{c_+\lambda_+^{\alpha_+}\Gamma({-}\alpha_+), c_-\lambda_-^{\alpha_-}\Gamma({-}\alpha_-)\}\\ &\qquad \ge c_+\lambda_+^{\alpha_+}(2-2^{\alpha_+})\Gamma({-}\alpha_+) +c_-\lambda_-^{\alpha_-}(2-2^{\alpha_-})\Gamma({-}\alpha_-).\end{align*}

Assuming that $\alpha_\pm=\alpha$ , the inequality simplifies to $\alpha\le\phi(\varrho)$ , where we define $\phi(x)\,:\!=\,\log_2\big(1+\frac{x}{1+x}\big)$ and $\varrho \,:\!=\, \min\{c_+\lambda_+^{\alpha},c_-\lambda_-^{\alpha}\} /\max\{c_+\lambda_+^{\alpha},c_-\lambda_-^{\alpha}\}$ . In particular, a symmetric Lévy measure yields $\varrho = 1$ and $\phi(1)=\log_2(3/2)=0.58496\ldots$ , and a one-sided Lévy measure gives $\varrho=0$ and $\phi(0)=0$ .

Figure 8. The picture shows the map $\varrho\mapsto\phi(\varrho)$ , $\varrho\in[0,1]$ . Assuming $\alpha_\pm=\alpha$ and defining $\varrho \,:\!=\, \min\{c_+\lambda_+^{\alpha},c_-\lambda_-^{\alpha}\} /\max\{c_+\lambda_+^{\alpha},c_-\lambda_-^{\alpha}\}$ , TSB-Alg is preferable to SB-Alg when $(\varrho,\alpha)$ lies in the shaded region.

4.3.2. Comparison with SBG-Alg

Recall from Subsection 3.4.3 that TSB-Alg is preferable to SBG-Alg when $\max\{\alpha_+,\alpha_-\}\ge 1$ (as it is equivalent to $\beta_\ast\ge1$ ). On the other hand, if $\alpha_\pm<1$ (or equivalently, $\beta_\ast<1$ ), TSB-Alg outperforms SB-Alg if and only if $(1+T)e^{(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T}\le C_1+C_2T$ . For large enough T, SB-Alg will outperform TSB-Alg; however, it is generally hard to determine when this happens. In Figure 9 we illustrate the region of parameters $(T,\alpha)$ where TSB-Alg is preferable, assuming $\alpha_\pm=\alpha\in(0,1)$ and all other parameters are as in the USD/JPY (v1) currency pair in Table 3.

Figure 9. The shaded region is the set of points $(T,\alpha)$ where TSB-Alg is preferable to SBG-Alg assuming $\alpha_\pm=\alpha\in(0,1)$ and all other parameters are as in the USD/JPY (v1) currency pair in Table 3.

Figure 10. The shaded region is the set of points $(T,\varepsilon)$ where TSB-Alg is preferable to SBG-Alg when $\alpha_\pm=\alpha\in[1,2)$ and all other parameters are as in the MCD stock in Table 3.

4.4. Concluding remarks

TSB-Alg is an easily implementable algorithm for which optimal MLMC (and even unbiased) estimators exist. TSB-Alg combines the best of both worlds: it is applicable to all tempered stable processes (as is SBG-Alg in [Reference González Cázares and Mijatović16]), while preserving the geometric convergence of SB-Alg in [Reference González Cázares, Mijatović and Uribe Bravo18]. The only downside of TSB-Alg is the enlarged variance by the factor $\exp((\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T)$ . This factor is typically small (see the discussion in Subsection 4.2.1 above), and when it is large, easily implementable variance reduction techniques exist (see details in Subsection 3.5 above). These facts favour the use of the MLMC estimator based on TSB-Alg over its competitors when $(\mu_{2{\boldsymbol\lambda}}-2\mu_{\boldsymbol\lambda})T$ is not large (for more concrete rules of thumb, see the concluding paragraphs of Subsections 3.4.1, 3.4.2, and 3.4.3 above).

In practice, in the implementation of the MC and MLMC estimators of TSB-Alg, the leading constants of the bias and variance decay are hard to compute and often overestimate the error. The general practice in this situation (see e.g. [Reference Giles13, Section 2]) is to numerically estimate these constants along with the integers $n(\varepsilon)$ , N, and $(N_k)_{k\in\mathbb{N}}$ in Proposition 3.2. Such estimation requires few simulations for the first few levels and some extrapolation but typically performs well in practice. This is particularly true in our setting as the MLMC estimator is optimal; see [Reference Giles13, Section 3] for a detailed discussion and a generic MATLAB implementation.