On approximations of subordinators in $L^p$ and the simulation of tempered stable distributions

Abstract

Subordinators are infinitely divisible distributions on the positive half-line. They are often used as mixing distributions in Poisson mixtures. We show that appropriately scaled Poisson mixtures can approximate the mixing subordinator and we derive a rate of convergence in $L^p$ for each $p\in [1,\infty ]$. This includes the Kolmogorov and Wasserstein metrics as special cases. As an application, we develop an approach for approximate simulation of the underlying subordinator. In the interest of generality, we present our results in the context of more general mixtures, specifically those that can be represented as differences of randomly stopped Lévy processes. Particular focus is given to the case where the subordinator belongs to the class of tempered stable distributions.

1 Introduction

Subordinators are a wide class of positive random variables. They correspond to the class of infinitely divisible distributions on the positive half-line. See Sato (1999) or Steutel and van Harn (2004) for many properties and examples. Subordinators have many uses and applications. They are often used as mixing distributions in the context of Poisson mixtures, which are common models for over-dispersed count data. We show that appropriately scaled Poisson mixtures can approximate the mixing subordinator and we derive a rate of convergence in $L^p$ for each $p\in [1,\infty ]$. When $p=\infty$ this is the Kolmogorov metric and when $p=1$ it is the Kantorovich or Wasserstein metric. As an application, we develop a methodology for approximate simulation of the underlying subordinator. In the interest of generality, we present our theoretical results in the context of more general mixtures. Specifically, those that can be represented as differences of randomly stopped Lévy processes.

We are particularly interested by the case where the mixing distribution is tempered stable (TS). TS distributions form a large class of models, which modify the tails of infinite variance stable distributions to make them lighter. This leads to more realistic models for many application areas. We are motivated by applications to finance, where TS subordinators are common models for stochastic volatility (Barndorff-Nielsen and Shephard 2001, 2003) and their bilateral versions are often used to model financial returns (Küchler and Tappe 2013, 2014). Tweedie distributions introduced in Tweedie (1984) are, perhaps, the earliest class of TS distributions to be studied. More general classes were introduced in Rosiński (2007) and Rosiński and Sinclair (2010). See also the monograph (Grabchak 2016) and the references therein. Poisson mixtures, where the mixing distribution is TS are called discrete TS distributions, see Grabchak (2018) and Grabchak (2022). In the special case where the mixing distribution belongs to the Tweedie class, they are called Poisson-Tweedie. See Baccini et al. (2016) for a recent overview.

Formally, a Poisson mixture can be described as follows. Let X be a subordinator and, independent of X, let $\{Z(t):t\ge 0\}$ be a Poisson process with rate 1. The Poisson mixture is then given by Z(X). For a survey of Poisson mixtures see Karlis and Xekalaki (2005) or Sect. VI.7 of Steutel and van Harn (2004). We consider a more general bilateral situation, where $X_+,X_-$ are independent subordinators and, independent of these, $Z^+=\{Z^+(t):t\ge 0\}$ and $Z^-=\{Z^-(t):t\ge 0\}$ are independent Lévy processes. Later we will allow these processes to be quite general, but, for the moment, assume that they are independent Poisson processes each with rate 1 and set

$$\begin{aligned} a Z^+(X^+/a)- a Z^-(X^-/a). \end{aligned}$$

(1)

This gives a random variable whose support is concentrated on the lattice $a\mathbb {Z}=\{0,\pm a,\pm 2a,\pm 3a,\dots \}$.

Distributions of this type are important in the context of high frequency finance. In practice, trades on an exchange cannot be for arbitrary amounts. There is a smallest price increment called the ‘tick size,’ and all prices changes must be in the form of an integer times the tick size. Thus, if a is the tick size, then price changes must take values on the lattice $a\mathbb {Z}$. For this reason, distributions of the form (1) are often used to model high frequency price changes, see Barndorff-Nielsen et al. (2012) or Koopman et al. (2017). On the other hand, less frequent price changes are typically modeled using continuous distributions, often bilateral subordinators. The idea is that over larger time frames, the tick size does not matter (so long as it is relatively small). We explain the transition from discrete to continuous models, by showing that

$$\begin{aligned} a Z^+(X^+/a)- a Z^-(X^-/a) {\mathop {\rightarrow }\limits ^{d}}X^+-X^- \text{ as } a\downarrow 0. \end{aligned}$$

(2)

Further, this result gives an approximate simulation method for simulating subordinators. The ability to simulate subordinators is particularly important as these distributions form the building blocks from which other distributions, including many multivariate distributions, can be simulated, see Xia and Grabchak (2022).

The purpose of this paper is two-folds. First, we prove that (2) holds and we derive a rate of convergence in $L^p$ for each $p\in [1,\infty ]$. In the interest of generality, we do not require $Z^+$ and $Z^-$ to be Poisson processes, but allow them to be more general Lévy processes. Second, we use this convergence to justify an approximate simulation method. We show that one can approximately simulate $X^+-X^-$ by simulating $a Z^+(X^+/a)- a Z^-(X^-/a)$ for a relatively small a and we derive approaches for simulating the latter. This is particularly important for TS distributions, where, with the exception of some relatively simple special cases (Hofert 2011; Grabchak 2021a), no exact simulation methods are known.

The rest of this paper is organized as follows. In Sect. 2 we recall some basic facts about infinitely divisible distributions and subordinators and in Sect. 3 we give our main theoretical results. In Sect. 4 we discuss Poisson mixtures more broadly and give our methodology for approximate simulation of the mixing distribution. In Sect. 5 we give a small scale simulation study. Proofs are postponed to Sect. 6. In Sect. 7 we consider extension of our results. These include an application to normal variance models, which are popular in financial applications.

Before proceeding we introduce some notation. We use the standard abbreviations: iid for independent and identically distributed, cdf for cumulative distribution function, pdf for probability density function, and pmf for probability mass function. We write $\mathbb {R}$ to denote the set of real numbers and we write ${\mathfrak {B}}(\mathbb {R})$ to denote the class of Borel sets on $\mathbb {R}$. We write $\mathbb {C}$ to denote the set of complex numbers and, for $z\in \mathbb {C}$, we write $\Re z$ and $\Im z$ to denote the real and imaginary parts of z, respectively. We write $\wedge$ and $\vee$ to denote the minimum and maximum, respectively. We write $\textrm{Pois}(\lambda )$, U(a, b), N(0, 1), and $\delta _a$ to denote, respectively, the Poison distribution with mean $\lambda$, the uniform distribution on (a, b), the standard normal distribution, and the point mass at a. Further, we write $\textrm{Ga}(\beta ,\zeta )$ to denote a gamma distribution with pdf

$$\begin{aligned} f_{\beta ,\zeta }(x)=\frac{\zeta ^{\beta }}{\Gamma (\beta )}x^{\beta -1}e^{-\zeta x},\ x>0, \end{aligned}$$

where $\beta ,\zeta >0$ are parameters. For a set A we write $1_A$ to denote the indicator function on A and we write ${:=}$ to denote a defining equality. If $\mu$ is a probability distribution, we write $X\sim \mu$ to denote that X is a random variable with distribution $\mu$ and we write $X_1,X_2,\dots {\mathop {\sim }\limits ^{\textrm{iid}}}\mu$ to denote that $X_1,X_2,\dots$ are iid random variables with common distribution $\mu$.

2 Background

In this paper we focus on infinitely divisible distributions with finite variation and no drift. A distribution $\mu$ of this type has a characteristic function of the form

$$\begin{aligned} {\hat{\mu }}(s) =e^{C_\mu (s)}, \ \ s\in \mathbb {R}, \end{aligned}$$

where

$$\begin{aligned} C_\mu (s)=\int _{-\infty }^\infty \left( e^{ixs}-1\right) M(\textrm{d}x), \ \ s\in \mathbb {R} \end{aligned}$$

is the cumulant generating function and M is the Lévy measure satisfying

$$\begin{aligned} M(\{0\}) = 0 \text{ and } \int _{-\infty }^\infty \left( |x|\wedge 1\right) M(\textrm{d}x)<\infty . \end{aligned}$$

(3)

We write $\mu =\textrm{ID}_0(M)$ to denote this distribution. For future reference, we note that

$$\begin{aligned} \Re C_\mu (s)=\int _{-\infty }^\infty \left( \cos (xs)-1\right) M(\textrm{d}x)\le 0, \ \ s\in \mathbb {R}. \end{aligned}$$

(4)

Associated with every infinitely divisible distribution $\mu$ is a Lévy process $\{Z(t):t\ge 0\}$, where $Z(1)\sim \mu$. This process has independent and stationary increments. An important example is when $M(\textrm{d}x) = \lambda \delta _1(\textrm{d}x)$. In this case, $\textrm{ID}_0(M)$ is a Poisson distribution with mean $\lambda$ and the associated Lévy process is a Poisson process with rate $\lambda$. See Sato (1999) for many results about infinitely divisible distributions and their associated Lévy processes.

When $M((-\infty ,0))=0$, the distribution $\mu =\textrm{ID}_0(M)$ is called a subordinator. By a slight abuse of terminology, we also use this term to describe any random variable with this distribution. In this case, $\mu ((-\infty ,0))=0$ and the associated Lévy process is increasing. When $M((0,\infty ))=0$, the distribution $\mu =\textrm{ID}_0(M)$ satisfies $\mu ((0,\infty ))=0$ and the associated Lévy process is decreasing. In general, under our assumptions, we can decompose the Lévy measure M into two parts:

$$\begin{aligned} M^+(B) = M(B\cap (0,\infty )),\ \ B\in {\mathfrak {B}}(\mathbb {R}) \end{aligned}$$

(5)

and

$$\begin{aligned} M^-(B) = M((-B)\cap (-\infty ,0)), \ \ B\in {\mathfrak {B}}(\mathbb {R}). \end{aligned}$$

(6)

In this case, if $X^-\sim \textrm{ID}_0(M^-)$ and $X^+\sim \textrm{ID}_0(M^+)$ are independent and $X=X^+-X^-$, then $X\sim \textrm{ID}_0(M)$ and $X^-,X^+$ are subordinators. Since X is the difference of two independent subordinators, we sometimes call it (or equivalently its distribution) a bilateral subordinator.

TS distributions form an important class of infinitely divisible distributions. They correspond to the case, where the Lévy measure is of the form

$$\begin{aligned} M(\text {d}x)&= \eta _- g(x) |x|^{-1-\alpha } 1_{[x<0]}\text {d}x\nonumber \\ &\quad +\eta _+ g(x) x^{-1-\alpha } 1_{[x>0]}\text {d}x, \end{aligned}$$

(7)

where $\eta _{\pm }\ge 0$ with $\eta _++\eta _->0$, $\alpha \in (0,1)$, and $g:\mathbb {R}\mapsto [0,\infty )$ is a bounded Borel function satisfying

$$\begin{aligned} \lim _{x\rightarrow 0}g(x) = 1. \end{aligned}$$

We call g the tempering function and we denote the distribution $\mu =\textrm{ID}_0(M)$ by $\mu =\textrm{TS}_{\alpha }(g,\eta _-,\eta _+)$. We note that one can consider the more general case where $\alpha \in (0,2)$. However, when $\alpha \in [1,2)$, the corresponding distributions do not have finite variation and are, thus, beyond the scope of this paper.

When $g\equiv 1$, the distribution $\textrm{TS}_{\alpha }(g,\eta _-,\eta _+)$ reduces to an infinite variance stable distribution, see Samorodnitsky and Taqqu (1994). It is often assumed that $\lim _{x\rightarrow \pm \infty }g(x)=0$, which leads to distributions that are similar to stable distributions in some central region, but with lighter (i.e., tempered) tails. This justifies the name and makes these distributions useful for many application areas, see Grabchak and Samorodnitsky (2010) for a discussion in the context of finance. Many facts about TS distributions can be found in Rosiński (2007), Rosiński and Sinclair (2010), Grabchak (2016), and the references therein.

3 Main results

Our main theoretical results are concerned with randomly stopped Lévy processes. For this reason, we generally work with two infinitely divisible distributions: $\mu =\textrm{ID}_0(M)$ and $\nu =\textrm{ID}_0(L)$. For us, $\nu$ is the distribution of the Lévy process, while $\mu$ is the distribution that governs the random time when the process is stopped. It is not the distribution of this random time, but it will be decomposed into two parts. The first will be the distribution of the time when we stop the positive jumps and second will be the distribution of the time when we stop the negative jumps. Throughout, we take

$$\begin{aligned} X\sim \mu =\textrm{ID}_0(M) \text{ and } Z\sim \nu =\textrm{ID}_0(L) \end{aligned}$$

to denote generic random variables from distributions $\mu$ and $\nu$, respectively. At this point, we make no additional assumptions on $\mu$, but for distribution $\nu$, we assume, throughout, that

$$\begin{aligned} \zeta _1{:=}\int _{-\infty }^\infty |z|L(\textrm{d}z)<\infty \text{ and } \gamma {:=}\int _{-\infty }^\infty z L(\textrm{d}z)>0. \end{aligned}$$

It can be shown that $\gamma =\textrm{E}[Z]$.

We decompose M into $M^+$ and $M^-$ as in (5) and (6). Let $X^-\sim \textrm{ID}_0(M^-)$ and $X^+\sim \textrm{ID}_0(M^+)$ be independent random variables and note that $X{\mathop {=}\limits ^{d}}X^+-X^-$. Independent of these, let $Z^+=\{Z^+(t):t\ge 0\}$ and $Z^-=\{Z^-(t):t\ge 0\}$ be independent Lévy processes with $Z^+(1),Z^-(1)\sim \nu$. For $a>0$ set

$$\begin{aligned} Y_a = a Z^+(X^+/(a\gamma ))- a Z^-(X^-/(a\gamma )) \end{aligned}$$

and let $\mu _a$ be the distribution of $Y_a$. A simple conditioning argument shows that the characteristic function of $Y_a$ is given by

$$\begin{aligned} & {\hat{\mu }}_a(s) = \textrm{E}\left[ e^{isY_a}\right] = \exp \left\{ \int _{-\infty }^\infty \left( e^{\frac{|x|}{a\gamma } C_\nu \left( sa\frac{x}{|x|}\right) }-1\right) M(\textrm{d}x)\right\} , \nonumber \\ & s\in \mathbb {R}, \end{aligned}$$

(8)

where $C_\nu$ is the cumulant generating function of $\nu$.

Proposition 1

We have

$$\begin{aligned} Y_a{\mathop {\rightarrow }\limits ^{d}}X \ \text{ as } \ a\downarrow 0. \end{aligned}$$

In the case where $Z^-$ and $Z^+$ are Poisson processes, versions of this result can be found in Klebanov and Slámova (2013), Grabchak (2018), and Proposition 6.5 of Steutel and van Harn (2004). Our goal is to characterize the rate of convergence in $L^p$. Let F be the cdf of distribution $\mu$ and let $F_a$ be the cdf of distribution $\mu _a$. Recall that, for $p\in [1,\infty )$, the $L^p$ distance between these distributions is given by

$$\begin{aligned} \Vert F-F_a\Vert _p = \left( \int _{-\infty }^\infty \left| F(x)-F_a(x)\right| ^p\textrm{d}x\right) ^{1/p} \end{aligned}$$

and for $p=\infty$ it is given by

$$\begin{aligned} \Vert F-F_a\Vert _\infty = \sup _{x\in \mathbb {R}} \left| F(x)-F_a(x)\right| . \end{aligned}$$

Since $\left| F(x)-F_a(x)\right| \le 1$ for each $x\in \mathbb {R}$, it is readily checked that for $1\le p_1\le p_2<\infty$ we have

$$\begin{aligned} \Vert F-F_a\Vert _{p_2}^{p_2} \le \Vert F-F_a\Vert ^{p_1}_{p_1}. \end{aligned}$$

In this context, the $L^\infty$ distance is sometimes called the Kolmogorov metric and the $L^1$ distance is sometimes called the Kantorovich or the Wasserstein metric. Many results about these and related distances can be found in Gibbs and Su (2002), Bobkov (2016), Arras and Houdré (2019), and the references therein. We just note the following equivalent formulations of $L^1$ distance. We have

$$\begin{aligned} & \Vert F-F_a\Vert _1 = \inf \{\textrm{E}\left[ |X-Y_a|\right] \}\\ & = \sup _h \left\{ \left| \textrm{E}\left[ h(X)\right] -\textrm{E}\left[ h(Y_a)\right] \right| \right\} , \end{aligned}$$

where the infimum is taken over all possible joint distributions of X and $Y_a$ and the supremum is taken over all continuous functions h satisfying the Lipschitz condition $|h(x)-h(y)|\le |x-y|$ for $x,y\in \mathbb {R}$.

Before stating our main results, we introduce several quantities. We do not, in general, assume these to be finite. Let

$$\begin{aligned} & r_0 {:=} \int _{\mathbb {R}} |{\hat{\mu }}(s)|\textrm{d}s, \ m_1 {:=} \int _{-\infty }^\infty |x| M(\textrm{d}x),\ m_2 \\ & {:=} \int _{-\infty }^\infty x^2 M(\textrm{d}x), \end{aligned}$$

and let

$$\begin{aligned} \zeta _2 {:=} \int _{-\infty }^\infty z^2 L(\textrm{d}z). \end{aligned}$$

It is well-known that if $r_0<\infty$, then $\mu$ has a pdf f satisfying

$$\begin{aligned} {\mathop {\mathrm{ess\,sup}}\limits _{x\in \mathbb {R}}}\, f(x)\le \frac{r_0}{2\pi }. \end{aligned}$$

Further, by Corollary 25.8 in Sato (1999) we have

$$\begin{aligned} m_i<\infty \text{ if } \text{ and } \text{ only } \text{ if } \textrm{E}[|X|^i]<\infty , \ \ i=1,2 \end{aligned}$$

and it is readily checked that, when they are finite, $m_1\ge \textrm{E}[|X|]$, $m_2=\textrm{Var}(X)$, and $\zeta _2=\textrm{Var}(Z)$. We now give our main result.

Theorem 1

Assume that $r_0,m_1,\zeta _1,\zeta _2$ are all finite. For any $a>0$,

$$\begin{aligned} \left\| F_{a} - F \right\| _\infty \le \sqrt{a}\left( \frac{1 }{\pi }e^{0.5m_1\zeta _2/\gamma } m_1\frac{\zeta _2}{2\gamma }+ \frac{12}{\pi ^2}\right) r_0= \mathcal {O}(a^{1/2}) \end{aligned}$$

and for any $p\in [2,\infty )$ and any $a>0$

$$\begin{aligned} & \left\| F_{a} - F \right\| _p \le \left( \sqrt{a} e^{0.5m_1\zeta _2/\gamma } m_1\frac{\zeta _2}{4\gamma } r_0^{1-1/p} \right. \\ & \quad \left. + a^{1/(2p)}4 (p-1)\right) = \mathcal {O}\left( a^{1/(2p)}\right) . \end{aligned}$$

If, in addition, $m_2<\infty$, then for any $p\in [1,\infty )$ and any $a>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p\le & \left\| F_{a} - F \right\| _1^{1/p} \le C_a^{1/p} a^{1/(2p)}= \mathcal {O}\left( a^{1/(2p)}\right) , \end{aligned}$$

where

$$\begin{aligned} & C_a = \left( \left( \sqrt{a} +1\right) e^{0.5m_1\zeta _2/\gamma } m_1+ \frac{\zeta _1}{\gamma } m_2 + 2\sqrt{a} m_1 + e^{0.5 m_1 \zeta _2/\gamma } m_1^2\frac{\zeta _1}{\gamma } \right) \\ & \frac{\zeta _2}{2\gamma } r_0^{1/2} + 4 \pi . \end{aligned}$$

It is well-known that the $L^1$ and $L^\infty$ distances provide bounds on a number of other distances, including the Lévy, discrepancy, and Prokhorov metrics, see Gibbs and Su (2002). Thus, our bounds give rates of convergence for these distances as well. We now give a better rate of convergence for the case where $\mu$ has a TS distribution.

Theorem 2

Let $\mu =\textrm{TS}_{\alpha }(g,\eta _-,\eta _+)$ and assume that $m_1,\zeta _1,\zeta _2$ are all finite. For any $a>0$,

$$\begin{aligned} \left\| F_{a} - F \right\| _\infty \le \mathcal {O}\left( a^{\frac{1}{2-\alpha }} \right) \end{aligned}$$

and for any $p\in [2,\infty )$ and any $a>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p \le \mathcal {O}\left( a^{\frac{1}{p(2-\alpha )}} \right) . \end{aligned}$$

If, in addition, $m_2<\infty$, then for any $p\in [1,\infty )$ and any $a>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p \le \left\| F_{a} - F \right\| _1^{1/p} \le \mathcal {O}\left( a^{\frac{1}{p(2-\alpha )}} \right) . \end{aligned}$$

In the proof we will see that we always have $r_0<\infty$ in this case.

4 Poisson mixtures and simulation

Our main results imply that $\mu _a \approx \mu$ when a is small, which suggests that we can approximately simulate from $\mu$ by simulating from $\mu _a$. In practice, simulation from $\mu _a$ is often easier than simulation from $\mu$. This is particularly true when the Lévy processes $Z^+,Z^-$ are Poisson processes. In this case, the distribution of $\mu _a$ is a scaled Poisson mixture. A general review of Poisson mixtures can be found in Karlis and Xekalaki (2005) or Sect. VI.7 of Steutel and van Harn (2004). Our discussion primarily extends results about discrete TS distributions given in Grabchak (2022).

Throughout this section we assume that $M((-\infty ,0])=0$, which means that $\mu =\textrm{ID}_0(M)$ is a subordinator. Further, we take $L(\textrm{d}z) = \delta _1(\textrm{d}z)$. This means that $\nu =\textrm{ID}_0(L)=\textrm{Pois}(1)$ and thus that the Lévy process $Z^+=\{Z^+(t):t\ge 0\}$ with $Z^+(1)\sim \nu$ is a Poisson process with rate 1 and that $\gamma =1$. Let $X\sim \textrm{ID}_0(M)$ be independent of $Z^+$ and, for any $a>0$, set

$$\begin{aligned} Y_a = a Z^+(X/a). \end{aligned}$$

As before, the distribution of $Y_a$ is denoted by $\mu _a$. In this case, we sometimes write $\mu _a=\textrm{PM}(M,a)$. Note that the support of this distribution is concentrated on $\{0,a,2a,3a,\dots \}$.

Specializing (8) to the current situation shows that the characteristic function of $Y_a$ is given by

$$\begin{aligned} {\hat{\mu }}_a(z) = \textrm{E}\left[ e^{izY_a}\right] = \exp \left\{ \int _{0}^\infty \left( e^{x (e^{iza}-1)/a}-1\right) M(\textrm{d}x)\right\} , \ \ \ z\in \mathbb {R}. \end{aligned}$$

Next, we introduce several useful quantities. These are closely related to the canonical sequence of discrete infinitely divisible distributions discussed in Steutel and van Harn (2004). For $a>0$, let

$$\begin{aligned} \ell _k^{(a)} = \int _{0}^{\infty } e^{-x/a} (x/a)^{k} M(\textrm{d}x), \quad k = 1,2, \ldots \end{aligned}$$

and let

$$\begin{aligned} \ell _+^{(a)} = \sum _{k=1}^{\infty } \frac{\ell _k^{(a)}}{k!} = \int _{0}^{\infty } (1-e^{-x/a}) M(\textrm{d}x). \end{aligned}$$

We now characterize the distribution $\mu _a$ in terms of these quantities.

Proposition 2

• 1. The characteristic function of $\mu _a$ is given by

  $$\begin{aligned} & {\hat{\mu }}_a(s) = \exp \left\{ \sum _{k=1}^{\infty } \frac{\ell _k^{(a)}}{k!}\left( e^{isa k} -1\right) \right\} \\ & = e^{-\ell _+^{(a)}} \exp \left\{ \sum _{k=1}^{\infty } \frac{\ell _k^{(a)}}{k!}e^{isa k} \right\} , \ \ s \in \mathbb {R}. \end{aligned}$$

• 2. We have $\mu _a=\textrm{ID}_0(M^{(a)})$, where the Lévy measure $M^{(a)}$ is concentrated on $\{a,2a,\ldots \}$ with

  $$\begin{aligned} M^{(a)}(\{a k\}) = \frac{1}{k!} \ell _k^{(a)}, \quad k = 1,2,\ldots , \end{aligned}$$

  and $M^{(a)} ((0, \infty )) = \ell _+^{(a)} < \infty$.

• 3. Let $p_k^{(a)} = P(Y_a = ak)$ for $k = 0, 1, 2, \ldots$. We have $p_0 = e^{- \ell _+^{(a)}}$ and

  $$\begin{aligned} p_k^{(a)} = \frac{1}{k} \sum _{j=0}^{k-1} \frac{\ell _{k-j}^{(a)}}{(k-j-1)!} p_j^{(a)}, \quad k = 1,2,\ldots . \end{aligned}$$

• 5. Let $M_1^{(a)}$ be the probability measure on $\{a,2a, \ldots \}$ with

  $$\begin{aligned} M_1^{(a)}(\{ak\}) = \frac{M^{(a)}(\{ak\})}{M^{(a)}((0, \infty ))} = \frac{\ell _k^{(a)}}{k! \ell _+^{(a)}}, \quad k = 1,2,\ldots . \end{aligned}$$

  If $N \sim \textrm{Pois}\left( \ell _+^{(a)}\right)$ and $W_1, W_2, \ldots {\mathop {\sim }\limits ^{iid}} M_1^{(a)}$ are independent of N, then

  $$\begin{aligned} Y_a{\mathop {=}\limits ^{d}}\sum _{i=1}^{N} W_i. \end{aligned}$$

Here, $M^{(a)}$ is a discretization of M. Parts 3 and 4 can be used to simulate from $\mu _a$. First we describe the algorithm suggested by Part 3. It is a version of the inverse transform method.

Algorithm 1: Simulation from $\mu _a=\textrm{PM}(M,a)$.

1. Simulate $U \sim U(0,1)$.

2. Evaluate $Y = \min \left\{ m=0,1,2, \ldots | \sum _{k=0}^m p^{(a)}_k>U\right\}$ and return aY.

This method often works well. However, for heavier tailed distributions, we may need to calculate $p^{(a)}_k$ for a large value of k, which can be inefficient. An alternate simulation method is suggested by Part 4.

Algorithm 2: Simulation from $\mu _a=\textrm{PM}(M,a)$.

1. Simulate $N \sim \textrm{Pois}(\ell _+^{(a)})$.

2. Simulate $W_1, W_2, \ldots , W_N {\mathop {\sim }\limits ^{iid}} M_1^{(a)}$.

3. Return $\sum _{i=1}^{N} W_i$.

To use Algorithm 2 we need a way to simulate from $M_1^{(a)}$. Following ideas in Grabchak (2022) for discrete TS distributions, we can represent $M_1^{(a)}$ as a mixture of scaled zero-truncated Poisson distributions. Let $\textrm{Pois}_+(\lambda )$ denote the zero truncated Poisson distribution with parameter $\lambda$ and pmf

$$\begin{aligned} p_{\lambda }(k) = \frac{\lambda ^{k} e^{- \lambda }}{k! (1-e^{-\lambda })}, \ \ k=1,2,\dots . \end{aligned}$$

We have

$$\begin{aligned} M_1^{(a)}(\{ak\})&= \frac{1}{\ell _+^{(a)}} \int _{0}^{\infty } \frac{(x/a)^k e^{-x/a}}{k! (1-e^{-x/a})} (1-e^{-x/a}) M(\textrm{d}x) \\&= \int _{0}^{\infty } p_{x/a}(k) M_2^{(a)}(\textrm{d}x), \end{aligned}$$

where

$$\begin{aligned} M_2^{(a)}(\textrm{d}x)= \frac{1}{\ell _+^{(a)}} (1-e^{-x/a}) M(\textrm{d}x) \end{aligned}$$

is a probability measure. Thus, we can simulate from $M_1^{(a)}$ as follows.

Algorithm 3: Simulation from $M_1^{(a)}$.

1. Simulate $\lambda \sim M_2^{(a)}$.

2. Simulate $W \sim \textrm{Pois}_+(\lambda )$ and return aW.

In this paper, to simulate from $\textrm{Pois}_+(\lambda )$, we use the rztpois method in the actuar package for the statistical software R, see Dutang et al. (2008). Simulation from $M_2^{(a)}$ depends on the structure of M. For certain discrete TS distributions and $a=1$ two methods were given in Grabchak (2022). We now extend them to the case of general $a>0$. For TS distributions, M is of the form (7) and, to ensure $M((-\infty ,0])=0$, we take $\eta _-=0$. It follows that, in this case,

$$\begin{aligned} M_2^{(a)}(\textrm{d}x)= \frac{\eta _+}{\ell _+^{(a)}} (1-e^{-x/a}) g(x)x^{-1-\alpha }1_{[x>0]}\textrm{d}x. \end{aligned}$$

Let

$$\begin{aligned} M_3^{(a)}(\textrm{d}x)= \frac{1}{K_a} (1-e^{-x}) g(xa)x^{-1-\alpha } 1_{[x>0]} \textrm{d}x, \end{aligned}$$

where $K_a=\int _0^\infty (1-e^{-x}) g(xa)x^{-1-\alpha } \text {d}x= a^\alpha \ell _+^{(a)}/\eta _+$ is a normalizing constant. It is easily checked that, if $Y\sim M_3^{(a)}$, then $aY\sim M_2^{(a)}$. Two rejection sampling methods for simulating from $M_3^{(a)}$ are given in Grabchak (2022). Toward this end, let $C_g=\operatorname {ess\,sup}_{x>0} g(x)$ and set

$$\begin{aligned} \varphi _1(x) = C_g^{-1}\frac{g(xa)(1-e^{-x})}{x\wedge 1}, \ \ x>0. \end{aligned}$$

Algorithm 4. Simulation from $M_2^{(a)}$ for TS distributions.

1. Independently simulate $U,U_1,U_2{\mathop {\sim }\limits ^{\textrm{iid}}}U(0,1)$ and set $V=U_1^{1/(1-\alpha )}U_2^{-1/\alpha }$.

2. If $U\le \varphi _1(V)$ return aV, otherwise go back to step 1.

On a given iteration, the probability of acceptance is $K_a \alpha (1-\alpha )/C_g$. In some cases, we can increase the acceptance rate by using a different trial distribution. The second method works when there exist $p,C,\zeta >0$ such that, for Lebesgue a.e. x,

$$\begin{aligned} g(x) \le C e^{-\zeta x^p}. \end{aligned}$$

(9)

In this case let

$$\begin{aligned} \varphi _2(x) = \frac{g(x)(1-e^{-x})}{Cxe^{-\zeta x^p}},\ x>0. \end{aligned}$$

We get the following rejection sampling algorithm, where, on a given iteration, the probability of acceptance is $\frac{K_a p\zeta ^\frac{1-\alpha }{p}}{C\Gamma \left( \frac{1-\alpha }{p}\right) }$.

Algorithm 5. Simulation from $M_2^{(a)}$ for TS distributions when (9) holds.

1. Independently simulate $U\sim U(0,1)$ and $V\sim \textrm{Ga}((1-\alpha )/p,\zeta )$.

2. If $U\le \varphi _2(V^{1/p})$ return $aV^{1/p}$, otherwise go back to step 1.

Note that, in this section, the assumption that $M((-\infty ,0])=0$ is without loss of generality. When it does not hold, we can decompose M into $M^+$ and $M^-$ as in (5) and (6). We can then independently simulate $Y_{a,+}\sim \textrm{PM}(M^+,a)$ and $Y_{a,-}\sim \textrm{PM}(M^-,a)$ and take

$$\begin{aligned} Y_a = Y_{a,+}-Y_{a,-}. \end{aligned}$$

5 Simulations

In this section we give a small simulation study to illustrate our methodology for approximate simulation from a subordinator $\mu$ by simulating from $\mu _a$ with some small $a>0$. We focus on two classes of TS distributions: the class of classical tempered stable (CTS) distributions and the class of TS distributions with power tempering (PT distributions). CTS distributions are, arguably, the best known and most commonly used TS distributions, see Küchler and Tappe (2013) for a survey. While exact simulation methods for CTS distributions are well-known, see e.g. Kawai and Masuda (2011) or Hofert (2011), we use them to illustrate our methodology because they are so popular. On the other hand, the only known simulation methods for PT distributions are either approximate or computationally intensive, see Grabchak (2019). For a discussion of PT distributions and their applications, see Grabchak (2016) or Grabchak (2021). In the symmetric case, methods to numerically evaluate pdfs, cdfs, and quantiles of CTS and PT distributions are given in the SymTS package for the statistical software R.

We begin with the class of CTS distributions. Distribution $\textrm{TS}_{\alpha }(\eta _-,\eta _+,g)$ is CTS if

$$\begin{aligned} g(x) = e^{-|x|/\ell _-}1_{[x\le 0]} + e^{-x/\ell _+} 1_{[x\ge 0]} \end{aligned}$$

for some $\ell _-,\ell _+>0$. It is often convenient to take $\eta _-=c_-\ell _-^\alpha$ and $\eta _+=c_+\ell _+^\alpha$, where $c_-,c_+\ge 0$ with $c_-+c_+>0$ are parameters. In this case we denote this distribution by CTS$_\alpha (c_-,\ell _-,c_+,\ell _+)$. When $c_-=0$ we get an important subset of the class of Tweedie distributions, see Tweedie (1984).

Our methodology for approximate simulation from $\mu =$CTS$_\alpha (c_-,\ell _-,c_+,\ell _+)$ is as follows. First choose $a>0$ small, then independently simulate $Y_{a,+}\sim \textrm{PM}(M^+,a)$ and $Y_{a,-}\sim \textrm{PM}(M^-,a)$, where $M^+(\textrm{d}x) = c_+\ell _+^\alpha e^{-x/\ell _+}x^{-1-\alpha } 1_{[x>0]}\textrm{d}x$ and $M^-(\textrm{d}x) = c_-\ell _-^\alpha e^{-x/\ell _-}x^{-1-\alpha } 1_{[x>0]}\textrm{d}x$. The random variable $Y_a = Y_{a,+}-Y_{a,-}$ gives an approximate simulation from $\mu$. We denote the distribution of $Y_a$ by $\mu _a$. We can simulate $Y_{a,+}$ and $Y_{a,-}$ using Algorithm 1 or Algorithm 2. In this paper, we use Algorithm 2, where we simulate from $M_1^{(a)}$ using Algorithm 3. To simulate from $M^{(a)}_2$ we can use Algorithm 4 or Algorithm 5. We select the one that has the larger probability of acceptance. Acceptance probabilities for several choices of the parameters are given in Table 1.

Table 1 Acceptance probabilities for Algorithms 4 and 5 for several choices of $\alpha$ and a

We now illustrate the performance of our methodology. For simplicity, we focus on the symmetric case. In this study we choose $\alpha \in \{0.25, 0.5, 0.75\}$, $c_-=c_+ = 1$, and $\ell _-=\ell _+ = 0.5$. We take $a \in \{0.5, 0.1, 0.01, 0.0001\}$ to demonstrate how $\mu _a$ approaches $\mu$. The results of our simulations are given in Fig. 1 and Table 2.

First, for each choice of the parameters and each value of a, we simulate 5000 observations. In Fig. 1 we give two types of diagnostic plots. First, we compare kernel density estimates (KDE) of our observations from $\mu _a$ with the density of $\mu$. We can see that, as a decreases, the KDE approaches the density of $\mu$. Second, we present qq-plots and again we see that, as a decreases, the qq-plots approach the diagonal line.

Next, we perform formal goodness-of-fit testing using the Kolmogorov-Smirnov (KS) and Cramer-Von Mises (CVM) tests. Here, for each choice the parameters and of a, we simulate 100 samples each comprised of 5000 observations. For each sample we evaluate the p-values of the KS and the CVM tests. We then average these together over the 100 samples. Recall that under the null hypothesis, the p-value follows a U(0, 1) distribution. Thus, when averaging over many runs, we expect the average of the p-values to be close to 0.5. The results are given in Table 2. We can see that the means of the p-values get closer to 0.5 as a decreases. Taken together, the goodness-of-fit results given in Fig. 1 and Table 2 suggest that the approximation works well for $a=10^{-4}$.

Fig. 1

Diagnostic plots for CTS. In the plots on the left, the solid black line is the pdf of the CTS distribution. Overlaid are KDEs based on observations from $\mu _a$ for several choices of a. The plots on the right give qq-plots based on samples from $\mu _a$ for several choices of a. The solid line is the reference line $y=x$

Table 2 These are p-values for KS and CVM tests for different values of $\alpha$ and a in approximating simulations from CTS and PT distributions

We now turn to PT distributions. Distribution $\textrm{TS}_{\alpha }(\eta _-,\eta _+,g)$ is PT if

$$\begin{aligned} g(x)= & 1_{[x\le 0]} \int _{0}^{\infty } e^{-|x|u} (1+u)^{-\alpha -\ell _{-}-2} u^{\ell _-} \textrm{d}u \\ & \quad + 1_{[x\ge 0]} \int _{0}^{\infty } e^{-xu} (1+u)^{-\alpha -\ell _+-2} u^{\ell _+} \textrm{d}u, \end{aligned}$$

where $\ell _+,\ell _->0$. It is often convenient to take $\eta _-=c_-(\alpha +\ell _-)(\alpha +\ell _-+1)$ and $\eta _+=c_+(\alpha +\ell _+)(\alpha +\ell _++1)$, where $c_-,c_+\ge 0$ with $c_-+c_+>0$ are parameters. In this case, we denote this distribution by PT$_\alpha (c_-,\ell _-,c_+,\ell _+)$. These distributions have finite means, but heavier tails than CTS distributions. If $Y\sim$PT$_\alpha (c_-,\ell _-,c_+,\ell _+)$ with $c_-,c_+>0$, then

$$\begin{aligned} \textrm{E}[|Y|^\beta ]<\infty \text{ if } \text{ and } \text{ only } \text{ if } \beta <1+\alpha +\min \{\ell _-,\ell _+\}. \end{aligned}$$

For more on these distributions, see Grabchak (2016) or Grabchak (2021).

Our methodology for approximate simulation from $\mu =$PT$_\alpha (c_-,\ell _-,c_+,\ell _+)$ is as follows. First choose $a>0$ small, then independently simulate $Y_{a,+}\sim \textrm{PM}(M^+,a)$ and $Y_{a,-}\sim \textrm{PM}(M^-,a)$, where $M^-(\textrm{d}x) = \eta _- g(-x) x^{-1-\alpha }1_{[x>0]}\textrm{d}x$ and $M^+(\textrm{d}x) = \eta _+ g(x)x^{-1-\alpha } 1_{[x>0]}\textrm{d}x$. The random variable $Y_a = Y_{a,+}-Y_{a,-}$ gives an approximate simulation from $\mu$ and we denote the distribution of $Y_a$ by $\mu _a$. We can simulate $Y_{a,+}$ and $Y_{a,-}$ using Algorithm 1 or Algorithm 2. In this paper, we use Algorithm 2, where we simulate from $M_1^{(a)}$ using Algorithm 3 and from $M^{(a)}_2$ using Algorithm 4. Note that we cannot use Algorithm 5 since (9) does not hold in this case.

We now illustrate the performance of our methodology for PT distributions. We again focus on the symmetric case. In this study we choose $\alpha \in \{0.25, 0.5, 0.75\}$, $c_-=c_+ = 1$, and $\ell _-=\ell _+ = 1$. We take $a\in \{0.5, 0.1, 0.01, 0.0001\}$ to demonstrate how $\mu _a$ converges to $\mu$. The results of our simulations are given in Fig. 2 and Table 2.

Again, for each choice of the parameters and each value of a, we simulate 5000 observations. In Fig. 2 we give two types of diagnostic plots. First, we compare the KDE of our observations from $\mu _a$ with the density of $\mu$. We can see that, as a decreases, the KDE approaches the density of $\mu$. Second, we present qq-plots and again we see that, as a decreases, the qq-plots approach the diagonal line.

Next, we perform formal goodness-of-fit testing. For each choice of the parameters and of a, we simulate 100 samples each comprised of 5000 observations. For each sample we evaluate the p-values of the KS and the CVM tests. We then average these together over the 100 samples. Again, when averaging over many runs, we expect the average of the p-values to be close to 0.5 under the null hypothesis. The results are given in Table 2. We can see that the means of the p-values get closer to 0.5 as a decreases. Taken together, the goodness-of-fit results given in Fig. 2 and Table 2 suggest that the approximation works well for $a=10^{-4}$.

Fig. 2

Diagnostic plots for PT. In the plots on the left, the solid black line is the pdf of the PT distribution. Overlaid are KDEs based on observations from $\mu _a$ for several choices of a. The plots on the right give qq-plots based on samples from $\mu _a$ for several choices of a. The solid line is the reference line $y=x$

6 Proofs

In this section we prove out theoretical results. We begin with several lemmas.

Lemma 1

Fix $z\in \mathbb {C}$.

1. We have $\left| e^{z}-1\right| \le |z|e^{|z|}$.

2. If $\Re z\le 0$, then $\left| e^{z}-1\right| \le 2\wedge |z|$.

3. If $\Im z=0$, then for any $n\in \{0,1,2,\dots \}$

$$\begin{aligned} \left| e^{iz}-\sum _{k=0}^n \frac{(iz)^k}{k!}\right| \le \min \left\{ \frac{|z|^{n+1}}{(n+1)!}, \frac{2|z|^n}{n!}\right\} . \end{aligned}$$

4. If $z_1,z_2\in \mathbb {C}$ with $\Re z_1,\Re z_2\le 0$, then

$$\begin{aligned} \left| e^{z_1}-e^{z_2}\right| \le 2\wedge |z_1-z_2|. \end{aligned}$$

Proof

We have

$$\begin{aligned} & \left| e^{z}-1\right| \le \sum _{n=1}^\infty \frac{|z|^n}{n!} = |z|\sum _{n=0}^\infty \frac{|z|^{n}}{(n+1)! }\\ & \le |z|\sum _{n=0}^\infty \frac{|z|^{n}}{n! }=|z|e^{|z|}. \end{aligned}$$

Next, let $z=-a+ib$, where $a=-\Re z$ and $b=\Im z$ and assume that $-a=\Re z\le 0$. We have

$$\begin{aligned} \left| e^{z}-1\right| = \left| z\int _0^1 e^{tz}\textrm{d}t\right| \le |z| \int _0^1 e^{-at}\textrm{d}t \le |z| \end{aligned}$$

and, by the triangle inequality,

$$\begin{aligned} \left| e^{z}-1\right| \le |e^{-a}| |e^{ib}| + 1 \le 2. \end{aligned}$$

The third part can be found in Sect. 26 of Billingsley (1995). We now turn to the fourth part. Without loss of generality assume that $\Re z_1\le \Re z_2\le 0$ and note that

$$\begin{aligned} \left| e^{z_1}-e^{z_2}\right| \le \left| 1-e^{z_1-z_2}\right| \le 2\wedge |z_1-z_2|, \end{aligned}$$

where the final inequality follows by the second part. $\square$

Lemma 2

For any $s\in \mathbb {R}$ we have

$$\begin{aligned} \left| C_\nu (s) \right| \le \zeta _1|s| \text{ and } \left| C_\nu (s) -is\gamma \right| \le \frac{1}{2}s^2 \zeta _2. \end{aligned}$$

Proof

By Lemma 1

$$\begin{aligned} & \left| C_\nu (s) \right| \le \int _{-\infty }^\infty \left| e^{isz}-1\right| L(\textrm{d}z)\\ & \le |s| \int _{-\infty }^\infty \left| z\right| L(\textrm{d}z) = |s|\zeta _1. \end{aligned}$$

Next, from the definition of $C_\nu$ and $\gamma$, we have

$$\begin{aligned} & \left| C_\nu (s) -is\gamma \right| \le \int _{-\infty }^\infty \left| e^{isz}-1-isz\right| L(\textrm{d}z)\\ & \le \frac{1}{2}s^2 \int _{-\infty }^\infty z^2L(\textrm{d}z), \end{aligned}$$

where the final inequality follows by Lemma 1. $\square$

Proof of Proposition 1

Our proof is based on showing the convergence of characteristic functions. By l’Hôpital’s rule

$$\begin{aligned} & \lim _{a\downarrow 0} \frac{|x|}{a\gamma } (e^{isazx/|x|}-1) \\ & = |x| \lim _{a\downarrow 0} \frac{\cos (sazx/|x|)-1+i\sin (sazx/|x|)}{a\gamma } = iszx/\gamma . \end{aligned}$$

By Lemma 1

$$\begin{aligned} \int _{-\infty }^\infty \frac{|x|}{a\gamma } \left| e^{isazx/|x|}-1\right| L(\textrm{d}z) \le |x||s|\int _{-\infty }^\infty |z|L(\textrm{d}z)/\gamma <\infty . \end{aligned}$$

Thus, by dominated convergence

$$\begin{aligned} \lim _{a\downarrow 0} \frac{|x|}{a\gamma } C_\nu (sax/|x|)= & \lim _{a\downarrow 0} \int _{-\infty }^\infty \frac{|x|}{a\gamma } \left( e^{isazx/|x|}-1\right) L(\textrm{d}z) \\= & ixs\int _{-\infty }^\infty zL(\textrm{d}z) /\gamma =ixs. \end{aligned}$$

From here

$$\begin{aligned} \lim _{a\downarrow 0}{\hat{\mu }}_a(z) = {\hat{\mu }}(z), \ \ \ z\in \mathbb {R} \end{aligned}$$

by another application of dominated convergence. We just note that by Lemmas 1 and 2

$$\begin{aligned} \left| e^{\frac{|x|}{a\gamma } C_\nu (sax/|x|) }-1\right|\le & 2\wedge \frac{|x|}{a\gamma } \left| C_\nu (sax/|x|)\right| \\\le & 2\wedge |x||s|\zeta _1/\gamma \\\le & \left( 2+|s|\zeta _1/\gamma \right) (1\wedge |x|). \end{aligned}$$

From here (3) guarantees that dominated convergence holds. $\square$

We now recall several versions of Esseen’s smoothing lemma, see Bobkov (2016) for a survey of such results.

Lemma 3

For any $T>0$, we have

$$\begin{aligned} \left\| F_{a} - F \right\| _\infty \le \frac{1}{\pi }\int _{-T}^T \left| \frac{{\hat{\mu }}_{a}(s) - {\hat{\mu }}(s)}{s} \right| \textrm{d}s+ \frac{12 r_0}{\pi ^2 T} \end{aligned}$$

(10)

and for any $T>0$ and any $p\in [1,\infty )$

$$\begin{aligned} \left\| F_{a} - F \right\| _p^p\le & \left\| F_{a} - F \right\| _1 \le \left( \int _{-T}^T\left| \frac{{\hat{\mu }}_{a}(s) - {\hat{\mu }}(s)}{s}\right| ^2\textrm{d}s\right) ^{1/2} \nonumber \\ & \quad + \left( \int _{-T}^T \left| \frac{\hat{\mu }'_{a}(s) - \hat{\mu }'(s)}{s}\right| ^2\textrm{d}s\right) ^{1/2} \nonumber \\ & \quad + \left( \int _{-T}^T \left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s^2}\right| ^2\textrm{d}s\right) ^{1/2} +\frac{4\pi }{T}. \end{aligned}$$

(11)

Further, for any $p\in [2,\infty )$, if $q = p/(p-1)$, then for any $T>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p \le \frac{1}{2}\left( \int _{-T}^T\left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s}\right| ^q\textrm{d}s\right) ^{1/q} + \frac{4(p-1)}{T^{1/p}}. \end{aligned}$$

(12)

Proof

Here, (10) is a version of the classical Esseen’s smoothing lemma, see, e.g., page 538 in Feller (1971) or Theorem 1.5.2 in Ibragimov and Linnik (1971). We just use the fact that the essential supremum of the pdf of $\mu$ is upper bounded by $r_0/(2\pi )$. The formula in (11) is essentially Corollary 8.3 in Bobkov (2016), which it itself a version of Theorem 1.5.4 in Ibragimov and Linnik (1971). The constant $4\pi$ is given in Ibragimov and Linnik (1971). To get the bound in (11), we apply the quotient rule and Minkowski’s inequality to the bound in Bobkov (2016). The formula in (12) is given in Corollary 7.2 of Bobkov (2016). It is a simple application of the classical Hausdorff-Young inequality, as given in, e.g., Proposition 2.2.16 of Grafakos (2004). $\square$

Lemma 4

1. For any $s\in \mathbb {R}$, we have

$$\begin{aligned} \left| \hat{\mu }_{a}(s) - \hat{\mu }(s)\right| \le a |\hat{\mu }(s)| s^2 e^{0.5 as^2 m_1\zeta _2/\gamma } m_1\zeta _2\frac{1}{2\gamma } \end{aligned}$$

and

$$\begin{aligned} & \left| \hat{\mu }'_{a}(s) - \hat{\mu }'(s)\right| \le a|\hat{\mu }(s)| \zeta _2\\ & \left( \frac{\zeta _1}{2\gamma ^2} s^2 m_2 + \frac{\left| s \right| }{\gamma } m_1 + s^2 e^{0.5 as^2 m_1 \zeta _2/\gamma } m_1^2 \frac{\zeta _1}{2\gamma ^2}\right) . \end{aligned}$$

2. For any $a>0$, $T>0$, and $q\ge 1$, we have

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s}\right| ^q\textrm{d}s\right) ^{1/q} \nonumber \\ & \le a T e^{0.5aT^2 m_1\zeta _2/\gamma } m_1\frac{\zeta _2}{2\gamma } r_0^{1/q}, \end{aligned}$$

(13)

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s^2}\right| ^q\textrm{d}s\right) ^{1/q} \\ & \le a e^{0.5aT^2 m_1\zeta _2/\gamma } m_1\frac{\zeta _2}{2\gamma } r_0^{1/q}, \end{aligned}$$

and

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }'_{a}(s) - \hat{\mu }'(s)}{s}\right| ^2\textrm{d}s\right) ^{1/2} \\ & \le a \frac{\zeta _2 }{\gamma } \left( \frac{\zeta _1}{2\gamma } T m_2 + m_1 + T e^{0.5aT^2 m_1 \zeta _2/\gamma } m_1^2\frac{\zeta _1}{2\gamma }\right) r_0^{1/2}. \end{aligned}$$

Clearly, for the bounds in the above to be meaningful, we need the quantities $r_0,m_1,m_2,\zeta _1$, and $\zeta _2$ be finite, where appropriate. For ease of presentation, we do not include this assumption in the statement of the lemma. This should not cause any confusion as the bounds are trivial when these quantities are infinite. Our proof of the first part is based on Lemma 1. A different approach for getting related bounds is used in the proof of Theorem 4.1 in Arras and Houdré (2019). However, that approach seems to lead to sub-optimal bounds in this case.

Proof

By the first part of Lemma 1 and the fact that $|e^{ixs}|=1$

$$\begin{aligned} & \left| \hat{\mu }_{a}(s) - \hat{\mu }(s)\right| \\ & = |\hat{\mu }(s)|\left| \exp \left\{ \int _{-\infty }^\infty \left( e^{C_\nu (sax/|x|)|x|/(a\gamma )}-e^{isx}\right) M(\textrm{d}x) \right\} -1\right| \\ & \quad \le |\hat{\mu }(s)| \exp \left\{ \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )-ixs}-1\right| M(\textrm{d}x) \right\} \\ & \quad \times \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )-ixs}-1\right| M(\textrm{d}x). \end{aligned}$$

Combining Lemmas 1 and 2 with (4) gives

$$\begin{aligned} & \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma ) -ixs}-1\right| M(\textrm{d}x)\nonumber \\ & \quad \le \int _{-\infty }^\infty \left| C_\nu (sax/|x|)|x|/(a\gamma ) -ixs\right| M(\textrm{d}x)\nonumber \\ & \quad = \frac{1}{a\gamma }\int _{-\infty }^\infty \left| C_\nu (sax/|x|) -isa\gamma x/|x|\right| |x| M(\textrm{d}x)\nonumber \\ & \quad \le \frac{a}{2\gamma }s^2\zeta _2\int _{-\infty }^\infty |x| M(\textrm{d}x) = \frac{a}{2\gamma } s^2 \zeta _2m_1. \end{aligned}$$

(14)

From here the first inequality follows. Turning to the second inequality, let

$$\begin{aligned} A_a(s) = \frac{1}{\gamma }\int _{-\infty }^\infty e^{C_\nu (sax/|x|)|x|/(a\gamma ) }C'_\nu (sax/|x|)xM(\textrm{d}x) \end{aligned}$$

and

$$\begin{aligned} A(s) = i\int _{-\infty }^\infty e^{ixs}xM(\textrm{d}x). \end{aligned}$$

By dominated convergence, we have $\hat{\mu }_a'(s) = A_a(s)\hat{\mu }_a(s)$ and $\hat{\mu }'(s) = A(s)\hat{\mu }(s)$. It follows that

$$\begin{aligned} & \left| \hat{\mu }'_{a}(s) - \hat{\mu }'(s)\right| \le \left| A_{a}(s) - A(s)\right| |\hat{\mu }(s)|\\ & +|A_a(s)|\left| \hat{\mu }_{a}(s) - \hat{\mu }(s)\right| . \end{aligned}$$

By dominated convergence,

$$\begin{aligned} C'_\nu (s) = i\int _{-\infty }^\infty e^{izs} zL(\textrm{d}z), \end{aligned}$$

and thus, for any $s\in \mathbb {R}$

$$\begin{aligned} |C'_\nu (s)| \le \int _{-\infty }^\infty |z|L(\textrm{d}z)=\zeta _1 \end{aligned}$$

and

$$\begin{aligned} & \left| A_a(s) \right| \le \frac{\zeta _1}{\gamma }\int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma ) }\right| |x|M(\textrm{d}x) \\ & \le \frac{1}{\gamma } m_1 \zeta _1, \end{aligned}$$

where the last inequality follows from (4). Next, using Lemmas 1 and 2 and (4) gives

$$\begin{aligned} & \left| A_a(s) -A(s) \right| \le \frac{1}{\gamma }\int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )}\frac{1}{i}C'_\nu (sax/|x|)\right. \\ & \quad \left. -\gamma e^{ixs}\right| |x|M(\textrm{d}x) \\ & \le \frac{1}{\gamma }\int _{-\infty }^\infty \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )}e^{izsax/|x|}\right. \\ & \left. - e^{ixs}\right| |z|L(\textrm{d}z)|x|M(\textrm{d}x) \\ & = \frac{1}{\gamma }\int _{-\infty }^\infty \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )+izsax/|x|-ixs}\right. \\ & \left. -1\right| |z|L(\textrm{d}z)|x|M(\textrm{d}x) \\ & \le \frac{1}{\gamma }\ \int _{-\infty }^\infty \int _{-\infty }^\infty \left| C_\nu (sax/|x|)\frac{|x|}{a\gamma }\right. \\ & \left. +izsa\frac{x}{|x|}-ixs\right| |z|L(\textrm{d}z) |x|M(\textrm{d}x) \\ & \le \frac{1}{a\gamma ^2} \int _{-\infty }^\infty \left| C_\nu (sax/|x|) -i\gamma as\frac{x}{|x|}\right| |x|^2M(\textrm{d}x)\\ & \int _{-\infty }^\infty |z|L(\textrm{d}z)\\ & \qquad + \frac{a\left| s \right| }{\gamma }\ \int _{-\infty }^\infty |x|M(\textrm{d}x)\int _{-\infty }^\infty |z|^2L(\textrm{d}z) \\ & \le \frac{a}{2\gamma ^2} s^2 m_2 \zeta _1\zeta _2 + \frac{a\left| s \right| }{\gamma } m_1\zeta _2. \end{aligned}$$

Putting everything together gives the first part. The second part follows immediately from the first, Minkowski’s inequality, and the fact that $|\hat{\mu }(s)|\le 1$ for each $s\in \mathbb {R}$. $\square$

Proof of Theorem 1

The proof follows by combining Lemmas 3 and 4. We illustrate the approach for the first bound only, as the rest can be verified in a similar way. By (10) for any $T>0$

$$\begin{aligned} & \left\| F_{a} - F \right\| _\infty \le \frac{1}{\pi }\int _{-T}^T \left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s} \right| \textrm{d}s+ \frac{12 r_0}{\pi ^2 T}\\ & \le \frac{1}{\pi }a T e^{0.5aT^2 m_1\zeta _2/\gamma } m_1\frac{\zeta _2}{2\gamma } r_0 + \frac{12 r_0}{\pi ^2 T}, \end{aligned}$$

where the second line follows by (13) with $q=1$. Taking $T=a^{-1/2}$ gives the result. $\square$

Proof of Theorem 2

First, fix $\delta \in (0,1)$ and note that, since $\lim _{x\rightarrow 0}g(x)=1$, there exists $\beta >0$ such that for any $x\in (-\beta ,\beta )$ we have $g(x)\ge 1-\delta$. It follows that for $|s|>1$

$$\begin{aligned} & |\hat{\mu }(s)| = \exp \left\{ - \int _{-\infty }^\infty \left( 1-\cos (|sx|)\right) M(\textrm{d}x)\right\} \\ & = \exp \left\{ -\int _0^\infty \left( 1-\cos (|s|x)\right) \left( \eta _+ g(x)\right. \right. \\ & \quad \left. \left. +\eta _- g(-x)\right) x^{-1-\alpha }\textrm{d}x\right\} \\ & = \exp \left\{ -|s|^\alpha \int _0^\infty \left( 1-\cos (x)\right) \left( \eta _+ g(x/|s|)\right. \right. \\ & \quad \left. \left. +\eta _- g(-x/|s|)\right) x^{-1-\alpha }\textrm{d}x\right\} \\ & \le \exp \left\{ -|s|^\alpha \int _0^\beta \left( 1-\cos (x)\right) \left( \eta _+ g(x/|s|)\right. \right. \\ & \quad \left. \left. +\eta _- g(-x/|s|)\right) x^{-1-\alpha }\textrm{d}x\right\} \\ & \le \exp \left\{ -|s|^\alpha \left( \eta _++\eta _-\right) (1-\delta ) \int _0^\beta \left( 1-\cos (x)\right) x^{-1-\alpha }\textrm{d}x\right\} \\ & = \exp \left\{ -|s|^\alpha A\right\} , \end{aligned}$$

where

$$\begin{aligned} A = (1-\delta )\left( \eta _++\eta _-\right) \int _0^\beta \left( 1-\cos (x)\right) x^{-1-\alpha }\textrm{d}x>0. \end{aligned}$$

Clearly $|\hat{\mu }(s)|$ is integrable and, hence, $r_0<\infty$ in this case.

Lemma 4 implies that for any $q\ge 1$

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }_a(s)-\hat{\mu }(s)}{s}\right| ^q\textrm{d}s\right) ^{1/q} \\ & \le am_1\frac{\zeta _2}{2\gamma }\left( \int _{-T}^T |s|^q |\hat{\mu }(s)|^q e^{0.5qas^2m_1\zeta _2/\gamma } \textrm{d}s\right) ^{1/q} \\ & \le am_1 \frac{\zeta _2}{\gamma }\left( \left( \int _{0}^1 s^q e^{qa m_1 \zeta _2/\gamma } \textrm{d}s \right) ^{1/q} \right. \\ & \ \left. + \left( \int _{1}^T s^qe^{-qs^\alpha A} e^{0.5qAs^{\alpha }s^{2-\alpha }am_1\zeta _2/(\gamma A)} \textrm{d}s\right) ^{1/q} \right) \\ & \le am_1 \frac{\zeta _2}{\gamma } \left( e^{a m_1 \zeta _2/\gamma } \right. \\ & \ + \left. \left( \int _{1}^\infty s^qe^{-qs^\alpha A } e^{0.5qAs^{\alpha }T^{2-\alpha }am_1\zeta _2/(\gamma A)} \textrm{d}s \right) ^{1/q}\right) , \end{aligned}$$

where the second line follows by Minkowski’s inequality and symmetry. Taking $q=1$ and combining this with Lemma 3 gives

$$\begin{aligned} & \Vert F-F_a\Vert _\infty \le \inf _{T>0}\left( \frac{m_1 }{\pi }a e^{a m_1 \zeta _2/\gamma } \frac{\zeta _2}{\gamma } \right. \\ & \quad \left. + \frac{m_1}{\pi }a \frac{\zeta _2}{\gamma } \int _{1}^\infty s e^{-s^\alpha A } e^{0.5As^{\alpha }T^{2-\alpha }am_1\zeta _2/(\gamma A)} \textrm{d}s +\frac{12 r_0}{\pi ^2 T}\right) \\ & \le a\frac{m_1}{\pi } e^{a m_1 \zeta _2/\gamma }\frac{\zeta _2}{\gamma } + a\frac{m_1}{\pi } \frac{\zeta _2}{\gamma } \int _{1}^\infty s e^{-s^{\alpha }A/2} \textrm{d}s\\ & \quad + a^{1/(2-\alpha )} \left( \frac{m_1\zeta _2}{A\gamma }\right) ^{1/(2-\alpha )} \frac{12 r_0}{\pi ^2} = \mathcal {O}\left( a^{1/(2-\alpha )}\right) , \end{aligned}$$

where we take $T = a^{-1/(2-\alpha )} \left( \frac{A\gamma }{m_1\zeta _2}\right) ^{1/(2-\alpha )}$ in the final inequality. The other parts of the theorem can be shown in a similar way. $\square$

Proof of Proposition 2

The first part follows from the fact that

$$\begin{aligned} & \sum _{k=1}^{\infty } \frac{\ell _k^{(a)}}{k!}\left( e^{iza k} -1\right) = \int _{0}^{\infty } e^{-x/a} \sum _{k=1}^{\infty } \frac{1}{k!}\left( e^{iza k} -1\right) (x/a)^{k} M(\textrm{d}x)\\ & = \int _{0}^{\infty } e^{-x/a} \sum _{k=1}^{\infty } \frac{1}{k!}\left( (e^{iza}x/a)^k -(x/a)^{k} \right) M(\textrm{d}x)\\ & = \int _{0}^{\infty } e^{-x/a} \left( e^{e^{iza}x/a} -e^{x/a} \right) M(\textrm{d}x) \\ & = \int _{0}^{\infty } \left( e^{(e^{iza}-1)x/a} - 1\right) M(\textrm{d}x). \end{aligned}$$

The second part follows immediately from the first. We now turn to the third part. We have

$$\begin{aligned} p_k^{(a)}= & P(Y_a = ak) = P(Z(X/a) = k) = \frac{1}{k!}\textrm{E}\left[ e^{-X/a}(X/a)^k\right] \\= & \frac{a^{-k}}{k!}\textrm{E}\left[ e^{-X/a}X^k\right] = \frac{a^{-k}}{k!} (-1)^k \psi ^{(k)}(1/a). \end{aligned}$$

where $\psi (s) = \textrm{E}[e^{-sX}]$ for $s\ge 0$ is the Laplace transform of X and $\psi ^{(k)}$ is the kth derivative of $\psi$. Formal substitution shows that

$$\begin{aligned} \psi (s) = \hat{\mu }(is)= \exp \left\{ \int _{0}^\infty \left( e^{-xs}-1\right) M(\textrm{d}x) \right\} . \end{aligned}$$

From here we see that $p_0^{(a)}=\psi (1/a)=e^{-\ell _+^{(a)}}$. Next, setting $f(s) = \log \psi (s)$ we note that the kth derivative of f is given by

$$\begin{aligned} f^{(k)}(s)= (-1)^k\int _{0}^\infty e^{-xs} x^k M(\textrm{d}x), \ \ k=1,2,\dots . \end{aligned}$$

It follows that $f^{(k)}(1/a) =(-a)^k \ell _k^{(a)}$. From standard results about derivatives of exponential functions, we get

$$\begin{aligned} \psi ^{(k)}(s) = \sum _{j=1}^{k-1} \left( {\begin{array}{c}k-1\\ j\end{array}}\right) \psi ^{(j)}(s)f^{(k-j)}(s), \ \ s>0. \end{aligned}$$

see, e.g., Lemma 1 in Grabchak (2022). It follows that

$$\begin{aligned} & p_k^{(a)} = \frac{a^{-k}}{k!} (-1)^k \psi ^{(k)}(1/a)\\ & = \frac{1}{k!} (-a)^{-k} \sum _{j=1}^{k-1} \left( {\begin{array}{c}k-1\\ j\end{array}}\right) \psi ^{(j)}(1/a)f^{(k-j)}(1/a)\\ & = \frac{1}{k!} (-a)^{-k} \sum _{j=1}^{k-1} \left( {\begin{array}{c}k-1\\ j\end{array}}\right) p_j^{(a)}(-a)^j j! (-a)^{k-j} \ell _{k-j}^{(a)}\\ & = \frac{1}{k} \sum _{j=1}^{k-1} \frac{1}{(k- j-1)!} p_j^{(a)} \ell _{k-j}^{(a)}, \end{aligned}$$

as required. The fourth part follows by a standard conditioning argument. $\square$

7 Extensions

In this section we consider some extensions of our work. These results are not complete, but aim to suggest directions for future research.

7.1 Infinite $r_0$

The results in Theorem 1 assume that $r_0=\int _{-\infty }^\infty |\hat{\mu }(s)|\textrm{d}s<\infty$, which implies that $\mu$ has a bounded density. Here, we briefly discuss what happens if we allow $r_0=\infty$. In this case, we no longer need $|\hat{\mu }(s)|$ in the bounds in Lemma 4. While we can apply $|\hat{\mu }(s)|\le 1$, we can get better bounds as follows. By (4), the fourth part of Lemma 1, and (14), we have for $s\in \mathbb {R}$

$$\begin{aligned} & \left| \hat{\mu }_{a}(s) - \hat{\mu }(s)\right| \\ & =\left| e^{\int _{-\infty }^\infty \left( e^{C_\nu (sax/|x|)|x|/(a\gamma )}-1\right) M(\textrm{d}x)}-e^{\int _{-\infty }^\infty \left( e^{isx}-1\right) M(\textrm{d}x)}\right| \\ & \le \int _{-\infty }^\infty \left| e^{C_\nu (sax/|x|)|x|/(a\gamma )}- e^{isx}\right| M(\textrm{d}x)\le \frac{a}{2\gamma } s^2 \zeta _2 m_1. \end{aligned}$$

From here, by arguments as in the proof of Lemma 4, we get for $s\in \mathbb {R}$

$$\begin{aligned} \left| \hat{\mu }'_{a}(s) - \hat{\mu }'(s)\right| \le a \zeta _2 \left( \frac{1}{2\gamma ^2} s^2 \zeta _1m_2 + \frac{\left| s \right| }{\gamma } m_1 + s^2 m_1^2\zeta _1\frac{1}{2\gamma ^2}\right) . \end{aligned}$$

It follows that for any $a>0$, $T>0$, and $q\ge 1$, we have

$$\begin{aligned} \left( \int _{-T}^T\left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s}\right| ^q\textrm{d}s\right) ^{1/q} \le aT^{1+1/q} \frac{2^{1/q-1}}{\gamma (q+1)^{1/q} } \zeta _2 m_1, \end{aligned}$$

(15)

$$\begin{aligned} \left( \int _{-T}^T\left| \frac{\hat{\mu }_{a}(s) - \hat{\mu }(s)}{s^2}\right| ^q\textrm{d}s\right) ^{1/q} \le aT^{1/q} 2^{1/q-1}\zeta _2 m_1, \end{aligned}$$

and

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }'_{a}(s) - \hat{\mu }'(s)}{s}\right| ^2\textrm{d}s\right) ^{1/2} \\ & \le a \sqrt{2} \zeta _2 \left( \frac{T^{3/2}}{2\sqrt{3} \gamma ^2} \zeta _1m_2 + \frac{T^{1/2}}{\gamma } m_1 + m_1^2\zeta _1\frac{T^{3/2}}{2\sqrt{3}\gamma ^2}\right) . \end{aligned}$$

Now, arguments as in the proof of Theorem 1 give the following.

Proposition 3

Assume that $m_1,\zeta _1,\zeta _2$ are all finite. For any $p\in [2,\infty )$ and any $a>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p \le \mathcal {O}\left( a^{1/(2p)}\right) . \end{aligned}$$

If, in addition, $m_2<\infty$, then for any $p\in [1,\infty )$ and any $a>0$

$$\begin{aligned} \left\| F_{a} - F \right\| _p\le & \mathcal {O}\left( a^{2/(5p)}\right) . \end{aligned}$$

Here, for $p\in [1,2)$, we get a slightly worse bound than in Theorem 1. We do not know if this due to an actually slower rate of convergence or if it is an artifact of the proof. We cannot give a bound for the $L^\infty$ metric without assuming $r_0<\infty$, as (10) requires this assumption. However, we can get a bound for the closely related Lévy metric, which metrizes weak convergence. Specifically, combining (15) with Theorem 3.2 in Bobkov (2016) leads to the rate $\mathcal {O}\left( a^{1/3}\log (1/a)\right)$.

7.2 Normal variance mixtures

In this subsection we give rates of convergence for approximating normal variance mixtures. A normal variance mixture is the distribution of a random variable of the form $\sqrt{X} W$, where $W\sim N(0,1)$ and X is a positive random variable independent of W. Such distributions often arise in financial application, see, e.g., Sabino (2023), Bianchi and Tassinari (2024), and the references therein.

Throughout this section, let $W\sim N(0,1)$ and let $Z=\{Z_t:t\ge 0\}$ be a Lévy process with $Z_1\sim \nu =\textrm{ID}_0(L)$. Assume that L is symmetric, i.e., that $L(B)=L(-B)$ for each $B\in {\mathfrak {B}}(\mathbb {R})$. This implies that $\nu$ is also symmetric and that its cumulant generating function $C_\nu$ is real. Assume further that

$$\begin{aligned} \int _{-\infty }^\infty \left( |z|\vee z^2\right) L(\textrm{d}z)<\infty \end{aligned}$$

(16)

and set

$$\begin{aligned} \gamma ^* = \int _{-\infty }^\infty z^2 L(\textrm{d}z). \end{aligned}$$

Let $X\sim \mu =\textrm{ID}_0(M)$ with $M((-\infty ,0])=0$ and

$$\begin{aligned} m_1 {:=} \int _0^\infty x M(\textrm{d}x)<\infty . \end{aligned}$$

Assume that X is independent of W and Z when appropriate. Let $\mu ^*$ be the distribution of $\sqrt{X} W$ and let $\mu _a^*$ be the distribution of $\sqrt{\frac{a}{\gamma ^*}} Z_{X/a}$. By a conditioning argument, it is easily checked that for any $a>0$ the characteristic function of $\mu ^*_{a}$ is

$$\begin{aligned} \hat{\mu }^*_a(s) = \exp \left\{ \int _0^\infty \left( e^{C_\nu (s\sqrt{a/\gamma ^*})x/a}-1\right) M(\textrm{d}x)\right\} \quad s\in \mathbb {R} \end{aligned}$$

and the characteristic function of $\mu ^*$ is

$$\begin{aligned} \hat{\mu }^*(s) = \exp \left\{ \int _0^\infty \left( e^{-s^2x/2}-1\right) M(\textrm{d}x)\right\} , \quad s\in \mathbb {R}. \end{aligned}$$

Proposition 4

We have

$$\begin{aligned} \sqrt{\frac{a}{\gamma ^*}} Z_{X/a} {\mathop {\rightarrow }\limits ^{d}}\sqrt{X} W \text{ as } a\downarrow 0. \end{aligned}$$

Proof

First note that the fact that L is symmetric and satisfies (16) implies that

$$\begin{aligned} & C_\nu (s{\sqrt{a/\gamma ^{*}}}) = \int _{-\infty }^{\infty }\left( e^{izs {\sqrt{a/\gamma ^{*}}}}-1\right) L(\textrm{d}z) \nonumber \\ & = \int _{-\infty }^{\infty }\left( e^{izs{\sqrt{a/\gamma ^{*}}} }-1-izs {\sqrt{\frac{a}{\gamma ^{*}}}} \right) L(\textrm{d}z). \end{aligned}$$

(17)

Now, by Lemma 1

$$\begin{aligned} \left| \frac{1}{a} \left( e^{izs{\sqrt{a/\gamma ^{*}}}}-1-izs {\sqrt{\frac{a}{\gamma ^{*}}}} \right) \right| \le \frac{z^2s^2 }{2 \gamma ^*} \end{aligned}$$

(18)

and

$$\begin{aligned} & \frac{1}{a} \left| e^{izs\sqrt{a/\gamma ^*}}-1-izs \sqrt{\frac{ a}{\gamma ^*}} + \frac{a}{2 \gamma ^*}z^2s^2 \right| \\ & \le \frac{\sqrt{a}}{6} \left| zs\right| ^3 (\gamma ^*)^{-3/2} \rightarrow 0 \text{ as } a\downarrow 0. \end{aligned}$$

It follows that

$$\begin{aligned} & \frac{1}{a} \left( e^{izs\sqrt{a/\gamma ^*}}-1-izs \sqrt{\frac{ a}{\gamma ^{*}}} \right) \\ & = -\frac{z^{2}s^{2}}{2\gamma ^{*}} +\frac{1}{a} \left( e^{izs\sqrt{a/\gamma ^*}}-1-izs \sqrt{\frac{a}{\gamma ^*}} + \frac{az^{2}s^{2}}{2\gamma ^{*}} \right) \end{aligned}$$

and by dominated convergence

$$\begin{aligned} & \lim _{a\downarrow 0} \frac{1}{a} C_{\nu }(s\sqrt{a/\gamma ^*}) \\ & = \lim _{a\downarrow 0} \frac{1}{a} \int _{-\infty }^{\infty }\left( e^{izs\sqrt{a/\gamma ^*}}-1-izs \sqrt{\frac{a}{\gamma ^*}}\right) L(\textrm{d}z)\\ & = -\frac{1}{2 \gamma ^{*}}s^{2} \int _{-\infty }^{\infty } z^{2} L(\textrm{d}z) = -\frac{1}{2}s^{2}. \end{aligned}$$

From here we get $\lim _{a\downarrow 0}\hat{\mu }^*_a(s) = \hat{\mu }^*(s)$ for each $s\in \mathbb {R}$ using another application of dominated convergence. To see that it holds note that by (4), Lemma 1, and (18), we have

$$\begin{aligned} & \left| e^{C_\nu (s\sqrt{a/\gamma ^*})x/a}-1\right| \le \left| C_\nu (s\sqrt{a/\gamma ^*})\frac{x}{a}\right| \\ & \le \frac{s^2 }{2 \gamma ^*}\int _{-\infty }^\infty z^2 L(\textrm{d}z) |x|, \end{aligned}$$

which is integrable with respect to M. $\square$

Henceforth, assume that

$$\begin{aligned} \zeta _3 {:=} \int _{-\infty }^\infty |z|^3 L(\textrm{d}z)<\infty \text{ and } r_0 {:=} \int _{-\infty }^\infty \left| \hat{\mu }^*(s)\right| \textrm{d}s <\infty . \end{aligned}$$

Theorem 3

Let $F_a$ be the cdf of $\mu _a^*$ and let F be the cdf of $\mu ^*$. Assume that $r_0,m_1,\zeta _3$ are all finite. For any $a>0$,

$$\begin{aligned} & \left\| F_{a} - F \right\| _\infty \le a^{1/6}\left( \frac{1}{6\pi } \left( \frac{1}{\sqrt{\gamma ^*}}\right) ^3 e^{\left( \gamma ^*\right) ^{-3/2} \zeta _3m_1/6} \zeta _3m_1 + \frac{12}{\pi ^2} \right) \\ & r_0= \mathcal {O}\left( a^{1/6}\right) \end{aligned}$$

and for any $p\in [2,\infty )$ and any $a>0$

$$\begin{aligned} & \left\| F_{a} - F \right\| _p \le \left( \frac{ a^{1/6}}{12} \left( \gamma ^*\right) ^{-3/2} e^{\left( \gamma ^*\right) ^{-3/2} \zeta _3m_1/6} \zeta _3m_1 r_0^{1-1/p}\right. \\ & \left. + 4(p-1) a^{1/(6p)} \right) \\ & = \mathcal {O}\left( a^{1/(6p)}\right) . \end{aligned}$$

For technical reasons, we are not able to get a rate of convergence for $p\in [1,2)$.

Proof

By arguments similar to those in the proof of Lemma 4, we have

$$\begin{aligned} & \left| \hat{\mu }^*_{a}(s) - \hat{\mu }^*(s)\right| \le |\hat{\mu }^*(s)| \exp \left\{ \int _{0}^\infty \left| e^{C_\nu (s\sqrt{a/\gamma ^*})x/a}\right. \right. \\ & \left. \left. -e^{-s^2x/2}\right| M(\textrm{d}x) \right\} \\ & \quad \times \int _{0}^\infty \left| e^{C_\nu (s\sqrt{a/\gamma ^*})x/a}-e^{-s^2x/2}\right| M(\textrm{d}x). \end{aligned}$$

By the fourth part of Lemma 1 and (17), we have

$$\begin{aligned} & \int _{0}^\infty \left| e^{C_\nu (s\sqrt{a/\gamma ^*})x/a}-e^{-s^2x/2}\right| M(\textrm{d}x) \\ & \quad \le \int _{0}^\infty \left| C_\nu (s\sqrt{a/\gamma ^*}) \frac{x}{a}+ \frac{s^2x \gamma ^*}{2\gamma ^*} \right| M(\textrm{d}x) \\ & \quad \le \int _{0}^\infty \frac{x}{a} \int _{-\infty }^\infty \left| e^{izs\sqrt{a/\gamma ^*} }-1-izs \sqrt{\frac{ a}{\gamma ^*}} + \frac{s^2z^2a}{2\gamma ^*} \right| \\ & L(\textrm{d}z) M(\textrm{d}x)\\ & \quad \le \frac{\sqrt{a}}{6}\left| \frac{s}{\sqrt{\gamma ^*}}\right| ^3\int _{-\infty }^\infty |z|^3 L(\textrm{d}z)\int _{0}^\infty xM(\textrm{d}x) \\ & = \frac{\sqrt{a}}{6} \left| \frac{s}{\sqrt{\gamma ^*}}\right| ^3 \zeta _3m_1. \end{aligned}$$

It follows that, for any $a>0$, $T>0$, and $q\ge 1$ we have

$$\begin{aligned} & \left( \int _{-T}^T\left| \frac{\hat{\mu }^*_{a}(s) - \hat{\mu }^*(s)}{s}\right| ^q\textrm{d}s\right) ^{1/q} \\ & \le \frac{\sqrt{a}}{6} T^2 \left( \frac{1}{\sqrt{\gamma ^*}}\right) ^3 e^{\frac{\sqrt{a}}{6} \left( \frac{T}{\sqrt{\gamma ^*}}\right) ^3 \zeta _3m_1} \zeta _3m_1 r_0^{1/q}. \end{aligned}$$

From here we proceed as in the proof of Theorem 1. $\square$