Remarks on Lévy Process Simulation

Abstract

Algorithms for simulation of a Lévy process X(t) are discussed, with particular emphasis on two algorithms approximating jumps that are in some sense small. One is classical, defining small jumps as those of absolute value ${<}\varepsilon$. The other one appears to be new and relies on an completely monotone structure of the Lévy density n(x). One then truncates the representing measure of n(x) to [0, A], meaning that jumps of mean ${<}1/A$ are left out. In both algorithms, the large jump part is simulated as compound Poisson and the small jumps are approximated. The standard choice of such an approximation is normal with the same mean and variance, but we also consider gamma approximations in two variants, and show that in some cases these perform substantially better. Other algorithms are briefly surveyed and we sketch a new one for simulation of a tempered stable (CGMY) process with infinite variation.

Appendix

1.1 An A-R Scheme for Generation from v(t)/t in the TS Case

We need to generate a r.v. Z with density proportional to u(t)/t, $\beta<t<A$ where $u(t)=(t-\beta )^\alpha$. To this end, write $Z=\beta +Z_0$ where $Z_0$ has density proportional to $u(t+\beta )/(t+\beta )$ $=t^\alpha /(t+\beta )$ , $0<t<B$ where $B=A-\beta$. Here $Y=1/Z_0$ has density proportional to $1/y^{1+\alpha }/(1+\beta y)$, $1<1/B<y<\infty$, and can therefore be generated by A-R with either a Pareto$(\alpha )$ proposal and acceptance probability proportional to $1/(1+\beta y)$ (high for small y) or a Pareto$(1+\alpha )$ proposal and acceptance probability proportional to $y/(1+\beta y)$ (high for large y). As in Sect. 4, we use a mixture, corresponding to breaking the compound Poisson part in 2) above into two. So, let

$$\begin{aligned}\begin{gathered} \lambda _1\ =\ \int _{1/B}^{y_0}\frac{1}{y^{1+\alpha }(1+\beta y)}\,{\textrm{d}} y\,,\ \ \mu _1\ =\ \int _{1/B}^{y_0}\frac{1}{y^{1+\alpha }}\,{\textrm{d}} y\ =\ \frac{1}{\alpha }[B^\alpha -1/y_0^\alpha ]\,,\\ \lambda _2\ =\ \int _{y_0}^\infty \frac{1}{y^{1+\alpha }(1+\beta y)}\,{\textrm{d}} y\,,\ \ \mu _2\ =\ \int _{y_0}^\infty \frac{1}{y^{2+\alpha }}\,{\textrm{d}} y\ =\ \frac{1}{(1+\alpha )y_0^{1+\alpha }}\,,\\ C_1\, =\,\frac{\mu _1}{\lambda _1(1+\beta /B)}\,,\quad C_2\, =\,\frac{\mu _2}{\beta \lambda _2} \end{gathered}\end{aligned}$$

The target densities are then

$$f_1(y)=\frac{1}{\lambda _1y^{1+\alpha }(1+\beta y)},\ 1/B<y<y_0,\text {\ \ and\ \ } f_2(y)=\frac{1}{\lambda _2y^{1+\alpha }(1+\beta y)},\ y_0<y<\infty ,$$

and chosen with probabilities $\lambda _1/(\lambda _1+\lambda _2)$, resp. $\lambda _2/(\lambda _1+\lambda _2)$, and the proposals are

$$g_1(y)=\frac{1}{\mu _1y^{1+\alpha }},\ 1/B<y<y_0,\text {\ \ and\ \ } g_2(y)=\frac{1}{\mu _2y^{2+\alpha }},\ y_0<y<\infty ,$$

Then $f_1(y)\le C_1(y_0)g_1(y)$ and $f_2(y)\le C_2(y_0)g_2(y)$, and we may use A-R with acceptance probabilities

$$\frac{f_1(y)}{C_1(y_0)g_1(y)}=\frac{1}{1+\beta y},\quad \frac{f_2(y)}{C_2(y_0)g_2(y)}=\frac{y}{1+\beta y}$$

for r.v. generation from $f_1$, resp. $f_2$. This gives expected numbers $C_1(y_0)$, $C_2(y_0)$ of samplings from $g_1(y)$, resp. $g_2(y)$, and as measure $E(y_0)$ of the computational effort, we shall use the total number of these samplings, i.e.

$$\begin{aligned} E(y_0)\ =\ \lambda _1C_1+\lambda _2C_2\ =\ \frac{\mu _1}{1+\beta /B}+\frac{\mu _2}{\beta }\ =\ \frac{\beta \mu _1+(1+\beta /B)\mu _2}{\beta (1+\beta /B)} \end{aligned}$$

(18)

Proposition 4

The function $E(y_0)$, $1/B<y_0<\infty$, is minimized for $y_0=y_0^*=1/\beta +1/B$.

Proof

In (20), $\beta$ and B as well as term $B^\alpha /\alpha$ in $\mu _1$ do not depend on $y_0$, so we are left with the minimization of $-\beta /\alpha /y_0^\alpha +(1+\beta /B)/(1+\alpha )/y_0^{1+\alpha }$. The derivative is $1/y_0^{1+\alpha }/(1+\beta /B)-1/y_0^{2+\alpha }$ which changes sign from negative to positive at $y_0^*$. From this the result follows. $\square$

1.2 An A-R Scheme for Spectrally Positive Infinite Variation TS Processes

Let $f_0$ be the density of a strictly $\alpha$-stable distribution $S_\alpha (\sigma ,1,0)$ distribution with $\sigma =\bigl (-\delta \Gamma (-\alpha )\cos (\pi \alpha /2)\bigr )^{1/\alpha }$. The goal is to generate a r.v. X from the density $f(x)=$ $\exp \{-\beta x-\psi \}f_0(x)$ in the case $\alpha >1$ where f and $f_0$ have support on the whole of ${\mathbb {R}}$; here $\psi =\delta \Gamma (-\alpha )\beta ^\alpha$. We use that $f_0(x)$ has the asymptotics [29, p. 100]

$$\begin{aligned} f_0(x)\ \sim \ \frac{c_1}{|x|^\ell }\exp \bigl \{-c_2|x|^\eta \bigr \}\ \ \text {as }x\rightarrow -\infty \end{aligned}$$

(19)

for suitable (explicit) constants) $c_1,c_2$ and $\ell =\alpha /(2\alpha -2)$, $\eta =\alpha /(\alpha -1)$.

For initialization of the algorithm:

(1) Select $-A<0$ and compute $p=\int _{-A}^\infty f(x)\,{\textrm{d}} x$

(2) Select $c_3<c_2$ and find $c_4<\infty$ such that

$$\begin{aligned} h(x)\ =\ \frac{{\textrm{e}}^{\beta |x|}f_0(x)}{|x|^{\eta -1}\exp \bigl \{-c_3|x|^\eta \bigr \}}\ \le \ c_4\text {\ \ \ for all }x<-A. \end{aligned}$$

The algorithm is then as follows:

(3) Generate I as Bernoulli(p).

(4) If $I=1$, generate $X\in (-A,\infty )$ having density f(x)/p, $-A<x<\infty$, by A-R with proposal $Z_0$ having a strictly $\alpha$-stable distribution $S_\alpha (\sigma ,1,0)$ conditioned to $(-A,\infty )$ and acceptance probability ${\textrm{e}}^{-b(z+A)}$ when $Z_0=z$

(5) If $I=0$, generate $X\in (-\infty ,-A)$ with density $\widetilde{f}(x)=$ ${\textrm{e}}^{b|x|-\psi }f_0(x)/(1-p)$, $=\infty<x<-A$, by an A-R scheme defined as follows. As proposal, take a r.v. $Z_1$ distributed as $-Z_2$ given $Z_2>A$ where $Z_2>0$ is Weibull with ${\mathbb {P}}(Z_2>z)=$ ${\textrm{e}}^{-c_3z^\eta }$. If $Z_1=x$, accept w.p. $c_4h(x)$.

(8) return X.

Explanation: Step (2) is possibly because (21), $c_2>c_3$ and $\eta >1$ imply $h(x)\rightarrow 0$ as $x\rightarrow -\infty$. In (5), the proposal density is $g(x)={\mathbb {P}}\bigl (X_2\in {\textrm{d}} |x|\,\big |\,Z_2>A\bigr )=c_3\eta |x|^{\eta -1}{\textrm{e}}^{-c_3z^\eta }/{\textrm{e}}^{-c_3A^\eta }$. Thus the ratio of the target density to the proposal density is

$$\begin{aligned} \frac{\widetilde{f}(x)}{g(x)}\ =\ c_5h(x)\quad \text {where}\quad c_5=\frac{\exp \{-\psi -c_3A^\eta \}}{c_3\eta (1-p)} \end{aligned}$$

Hence $\widetilde{f}(x)/g(x)\le c_0h(x)$ where $c_0=c_4c_5$, and acceptance w.p. $\widetilde{f}(x)/c_0/g(x)=c_4h(x)$ will produce the correct result. The conditioned sampling of proposals in (6) and (7) is straightforward by sampled by sampling a $S_\alpha (\sigma ,1,0)$, resp. Weibull, r.v. until the conditioning requirement is met. Available software, say Matlab or Nolan’s stable package (see the Preface to [29]) accounts for computing $f_0(x)$ and generating the $S_\alpha (\sigma ,1,0)$ r.v.’s. The Weibulls can be generated by inversion.