On the Simulation of General Tempered Stable Ornstein-Uhlenbeck Processes

TL;DR

This paper provides an explicit transition law for tempered stable Ornstein-Uhlenbeck processes and introduces an exact simulation method applicable to broad classes of univariate and multivariate distributions.

Contribution

It offers a novel explicit representation for the transition law and a rejection sampling algorithm for exact simulation of these processes.

Findings

Explicit transition law derived for tempered stable Ornstein-Uhlenbeck processes

Rejection sampling algorithm developed for exact simulation

Results encompass many previously studied cases as special instances

Abstract

We give an explicit representation for the transition law of a tempered stable Ornstein-Uhlenbeck process and use it to develop a rejection sampling algorithm for exact simulation of increments from this process. Our results apply to general classes of both univariate and multivariate tempered stable distributions and contain a number of previously studied results as special cases.

Figures

Equations

$$g(x;\alpha ,p,\delta )$$

$$(p-\alpha ){\delta }^{-(p-\alpha )}{x}^{p-\alpha -1},0<x\le \delta$$

$$\alpha {\delta }^{\alpha }{x}^{-1-\alpha },x>\delta .$$

$$h=\frac{\alpha {\delta }^p}{\alpha {\delta }^p+p-\alpha }.$$

$$Y=U_1^{1/(p-\alpha )}\delta 1_{[U_2<h]}+U_1^{-1/\alpha }\delta 1_{[U_2\ge h]}\sim \mathrm{MLL}(\alpha ,p,\delta ).$$

$$\frac{{\theta }^{\beta }}{\Gamma (\beta )}{u}^{\beta -1}{e}^{-\theta u},u>0,$$

$$\frac{p{\theta }^{\beta /p}}{\Gamma (\beta /p)}{u}^{\beta -1}{e}^{-\theta u^p},u>0,$$

$$C_{\mu }(z)=-\frac{1}{2}⟨z,Az⟩+i⟨b,z⟩+{\int }_{{\mathbb{R}}^d}\left(e^{i⟨z,x⟩}-1-i⟨z,x⟩h(x)\right)M(\mathrm{d}x).$$

$$M(\{0\})=0\text{mbox}and{\int }_{{\mathbb{R}}^d}(|x{|}^2\wedge 1)M(\mathrm{d}x)<\infty .$$

$${\int }_{{\mathbb{R}}^d}\left|e^{i⟨z,x⟩}-1-i⟨z,x⟩h(x)\right|M(\mathrm{d}x)<\infty$$

$$M(B)={\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }{1}_B(u\xi )k(\xi ,u)u^{-1}\mathrm{d}u\sigma (\mathrm{d}\xi ),B\in \mathfrak{B}({\mathbb{R}}^d).$$

$$X\stackrel{d}{=}\frac{1}{c}X+Y_c.$$

$$L_{\alpha }(B)={\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }{1}_B(u\xi )u^{-1-\alpha }q(\xi ,u)\mathrm{d}u\sigma (\mathrm{d}\xi ),B\in \mathfrak{B}({\mathbb{R}}^d),$$

$$\displaystyle h_{\alpha}(x)=\left\{\begin{array}[]{ll}0&\alpha\in(0,1)\\ 1_{[|x|\leq 1]}&\alpha=1\\ 1&\alpha\in(1,2)\end{array}\right.,$$

$$\exp \left\{i⟨b,z⟩+{\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }{\psi }_{\alpha }(z,u\xi )q(\xi ,u)u^{-1-\alpha }\mathrm{d}u\sigma (\mathrm{d}\xi )\right\},z\in {\mathbb{R}}^d,$$

$${\psi }_{\alpha }(z,x)=e^{i⟨z,x⟩}-1-i⟨z,x⟩h_{\alpha }(x)\text{mbox}forz,x\in {\mathbb{R}}^d.$$

$$\underset{u\to 0}{\lim }q(\xi ,u)=1\text{mbox}and\underset{u\to \infty }{\lim }q(\xi ,u)=0.$$

$$\mathrm{d}{Y}_t=-\lambda Y_t\mathrm{d}t+\mathrm{d}{Z}_t,$$

$$Y_t=e^{-\lambda t}{Y}_0+{\int }_0^t{e}^{-\lambda (t-s)}\mathrm{d}{Z}_s.$$

$${\int }_{|x|>2}\log |x|M(\mathrm{d}x)<\infty ,$$

$$c_{\alpha}=b-\left\{\begin{array}[]{ll}\int_{|x|>1}\frac{x}{|x|}M(\mathrm{d}x)&\mbox{if }\alpha=1\\ 0&\mbox{otherwise}\end{array}\right.,$$

$$M(B)={\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }{1}_B(u\xi )\mathrm{d}{\rho }_{\xi }(u)\sigma (\mathrm{d}\xi ),B\in \mathfrak{B}({\mathbb{R}}^d),$$

$$C_t(y,z)$$

$$d_{\alpha }$$

$${\tilde{d}}_1={\int }_{1<|x|\le e^{\lambda t}}x\left(|x{|}^{-1}-e^{-\lambda t}\right)M(\mathrm{d}x)-\left(1-e^{-\lambda t}\right){\int }_{|x|>1}\frac{x}{|x|}M(\mathrm{d}x).$$

$$K={\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }\left(q(\xi ,u)-q(\xi ,u{e}^{\lambda t})\right)u^{-1-\alpha }\mathrm{d}u\sigma (\mathrm{d}\xi ).$$

$$Y_{s+t}\stackrel{d}{=}{e}^{-\lambda t}y+d_{\alpha }+X_0$$

$$Y_{s+t}\stackrel{d}{=}{e}^{-\lambda t}y+d_{\alpha }+X_0+\sum _{j=1}^N{X}_j-\psi .$$

$$\psi =e^{-\alpha \lambda t}{\int }_{{\mathbb{S}}^{d-1}}{\int }_0^{\infty }{h}_{\alpha }(u\xi )\left(q(\xi ,u)-q(\xi ,u{e}^{\lambda t})\right)u^{-\alpha }\mathrm{d}u\xi \sigma (\mathrm{d}\xi )$$

$$H(\mathrm{d}\xi ,\mathrm{d}u)=\frac{1}{K}\left(q(\xi ,u)-q(\xi ,u{e}^{\lambda t})\right)u^{-1-\alpha }\mathrm{d}u\sigma (\mathrm{d}\xi ),u>0,\xi \in {\mathbb{S}}^{d-1},$$

$${\kappa }_{\xi }=\frac{1}{{\int }_0^{\infty }\left(q(\xi ,u)-q(\xi ,u{e}^{\lambda t})\right)u^{-1-\alpha }\mathrm{d}u},\xi \in {\mathbb{S}}^{d-1},$$

Full text

Keywords: Tempered stable distributions, Ornstein-Uhlenbeck processes, rejection sampling, selfdecomposability

1 Introduction

Non-Gaussian processes of Ornstein-Uhlenbeck-type (henceforth OU-processes) form a rich and flexible class of stochastic models. They are the continuous time analogues of autoregressive AR(1) processes and are important in the study of selfdecomposable distributions. Further, they are mean reverting, which makes them useful for many applications. We are particularly motivated by their applications to mathematical finance, where they have been used to model stochastic volatility [3] [4], stochastic interest rates [24] [19], and commodity prices [5] [12]. Here we focus on the class of TSOU-processes, which are OU-processes with tempered stable limiting distributions. These distributions are obtained by modifying the tails of infinite variance stable distributions to make them lighter, which leads to models that are more realistic for a variety of application areas. In particular, discussions of financial applications can be found in [18], [25], [14], and the references therein. Some of the earliest tempered stable distributions were introduced in the seminal paper Tweedie (1984) [32] and are called Tweedie distributions. A general framework was developed in Rosiński (2007) [27]. This was further generalized, in several directions, in [28], [6], and [13]. A survey, along with a historical overview and many references can be found in [14].

In this paper, we give an explicit representation for the transition law of a TSOU-process and use it to develop a simulation method based on rejection sampling. This method extends the rejection sampling technique for simulating tempered stable Lévy processes, which was introduced in [16]. Similar representations for the transition laws of TSOU-processes with Tweedie (and closely related) limiting distributions are given in [34], [33], [21], [22], and [7]. Most of these are special cases of our result.

The rest of this paper is organized as follows. In Section 2 we discuss two distributions that will be useful for simulating TSOU-processes and in Section 3 we recall the definition of the class of tempered stable distributions. Then, in Section 4 we formally introduce TSOU-processes, give explicit representations of their transition laws, and introduce our methodology for simulation. In Section 5 we specialize our results to the important class of $p$-tempered $\alpha$-stable distributions. A small simulation study is given in Section 6. Proofs are postponed to Section 7.

Before proceeding, we introduce some notation. Let $\mathbb{R}^{d}$ be the space of $d$-dimensional column vectors of real numbers equipped with the usual inner product $\langle\cdot,\cdot\rangle$ and the usual norm $|\cdot|$. Let $\mathbb{S}^{d-1}=\{x\in\mathbb{R}^{d}:|x|=1\}$ denote the unit sphere in $\mathbb{R}^{d}$. Let $\mathfrak{B}(\mathbb{R}^{d})$ and $\mathfrak{B}(\mathbb{S}^{d-1})$ denote the Borel sets in $\mathbb{R}^{d}$ and $\mathbb{S}^{d-1}$, respectively. If $a,b\in\mathbb{R}$, we write $a\vee b$ and $a\wedge b$ to denote, respectively, the maximum and the minimum of $a$ and $b$. We write $U(0,1)$ to denote the uniform distribution on $[0,1]$, we write $\delta_{x}$ to denote a point mass at $x$, and we write $1_{A}$ to denote the indicator function on set $A$. If $\mu$ is a probability measure on $\mathbb{R}^{d}$, we write $X\sim\mu$ to denote that $X$ is an $\mathbb{R}^{d}$-valued random variable with distribution $\mu$. We write $\stackrel{{\scriptstyle d}}{{=}}$ to denote equality in distribution and $\stackrel{{\scriptstyle w}}{{\rightarrow}}$ to denote weak convergence.

2 Auxiliary Distributions

In this section we briefly introduce two distributions, which will be important for simulating increments from TSOU-processes. These are the modified log-Laplace distribution and the generalized gamma distribution.

We begin with the modified log-Laplace distribution. Consider a distribution with probability density function (pdf) given by

[TABLE]

where $p>\alpha>0$ and $\delta>0$ are parameters. We call this the modified log-Laplace distribution and denote it by $\mathrm{MLL}(\alpha,p,\delta)$. This distribution is a mixture of the beta distribution with pdf

[TABLE]

and the Pareto distribution with pdf

[TABLE]

The mixing parameter is

[TABLE]

Using this mixture representation, we can simulate from the distribution as follows. If $U_{1}$ and $U_{2}$ are independent and identically distributed (iid) random variables from $U(0,1)$, then

[TABLE]

The reason that we call this the modified log-Laplace distribution is because the log-Laplace distribution is also a mixture of a beta distribution and a Pareto distribution, but with a different mixing parameter. The modified log-Laplace distribution is a log-Laplace distribution only when $\delta=1$. For more on the log-Laplace distribution see [23].

We now turn to the generalized gamma distribution. Recall that the pdf of a gamma distribution is given by

[TABLE]

where $\beta,\theta>0$ are parameters. We denote this distribution by $\mathrm{Ga}(\beta,\theta)$. The generalized gamma distribution was introduced in [30]. It has a pdf of the form

[TABLE]

where $\beta,p,\theta>0$ are parameters. We denote this distribution by $\mathrm{GGa}(\beta,p,\theta)$. It is readily checked that if $X\sim\mathrm{Ga}(\beta/p,\theta)$, then $X^{1/p}\sim\mathrm{GGa}(\beta,p,\theta)$. Thus, to simulate from a generalized gamma distribution it suffices to know how to simulate from a gamma distribution. Approaches for doing this are well-known, see e.g. [2] and the references therein.

3 Tempered Stable Distributions

An infinitely divisible distribution $\mu$ on $\mathbb{R}^{d}$ is a probability measure with a characteristic function of the form $\hat{\mu}(z)=\exp\{C_{\mu}(z)\}$, where, for $z\in\mathbb{R}^{d}$,

[TABLE]

Here, $A$ is a symmetric nonnegative-definite $d\times d$-dimensional matrix called the Gaussian part, $b\in\mathbb{R}^{d}$ is called the shift, and $M$ is a Borel measure, called the Lévy measure, which satisfies

[TABLE]

The function $h:\mathbb{R}^{d}\mapsto\mathbb{R}$, which we call the $h$-function, can be any Borel function satisfying

[TABLE]

for all $z\in\mathbb{R}^{d}$. For a fixed $h$-function, $A$, $M$, and $b$ uniquely determine the distribution $\mu$, and we write $\mu=\mathrm{ID}_{d}(A,M,b)_{h}$. The choice of $h$ does not affect $A$ and $M$, but different choices of $h$ require different values for $b$, see Section 8 in [29].

Associated with every infinitely divisible distribution $\mu=\mathrm{ID}_{d}(A,M,b)_{h}$ is a stochastic process $\{X_{t}:t\geq 0\}$, which is called a Lévy process. This processes is stochastically continuous and has independent and stationary increments. The characteristic function of $X_{t}$ is $\left(\hat{\mu}(z)\right)^{t}$. It follows that, for each $t\geq 0$, $X_{t}\sim\mathrm{ID}_{d}(tA,tM,tb)_{h}$. For more on infinitely divisible distributions and their associated Lévy processes see [29].

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

[TABLE]

Here $\sigma$ is a finite Borel measure on $\mathbb{S}^{d-1}$ and $k:\mathbb{S}^{d-1}\times(0,\infty)\mapsto[0,\infty)$ is a Borel function, which, for each fixed $\xi\in\mathbb{S}^{d-1}$, is nonincreasing in $u$. A more intuitive characterization is as follows. A distribution $\mu$ is selfdecomposable if and only if there is an $X\sim\mu$ such that, for all $c>1$, there exists a random variable $Y_{c}$, independent of $X$, satisfying

[TABLE]

Following [28], we define a tempered stable distribution on $\mathbb{R}^{d}$ as an infinitely divisible distribution with no Gaussian part and a Lévy measure of the form

[TABLE]

where $\alpha\in(0,2)$, $q:\mathbb{S}^{d-1}\times(0,\infty)\mapsto[0,\infty)$ is a Borel function, and $\sigma$ is a finite Borel measure on $\mathbb{S}^{d-1}$. The function $q$ is called the tempering function. Throughout this paper we assume that $q$ satisfies the following two assumptions:

• A1.

$0\leq q(\xi,u)\leq 1$ for all $u>0$ and $\xi\in\mathbb{S}^{d-1}$, and

• A2.

$q(\xi,u)$ is nonincreasing in $u$ for each fixed $\xi\in\mathbb{S}^{d-1}$.

Assumption A1 ensures that the tempering function is bounded, while Assumption A2 ensures that the corresponding distribution is selfdecomposable. For tempered stable distributions we use the $h$-function given by

[TABLE]

and we denote the distribution $\mathrm{ID}_{d}(0,L_{\alpha},b)_{h_{\alpha}}$ by $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$. When $q(\xi,u)\equiv 1$, the distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ is $\alpha$-stable and we sometimes write $\mathrm{S}_{\alpha,d}(\sigma,b)$ in this case. The characteristic function of $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ is given by

[TABLE]

where

[TABLE]

Remark 1**.**

The name “tempering function” can be explained as follows. It is often assumed that $q$ satisfies the additional assumption that for each $\xi\in\mathbb{S}^{d-1}$

[TABLE]

In this case, the tempered stable distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ has a distribution function that looks like that of the stable distribution $\mathrm{S}_{\alpha,d}(\sigma,b)$ in some central region, but with lighter tails. In this sense, $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ “tempers” the tails of $\mathrm{S}_{\alpha,d}(\sigma,b)$. While this is part of the motivation for defining tempered stable distributions, we will not require (3) to hold in this paper.

4 TSOU-Processes

We begin this section by recalling the definition of a process of Ornstein-Uhlenbeck-type (OU-process). Toward this end, let $Z=\{Z_{t}:t\geq 0\}$ be a Lévy process with $Z_{1}\sim\mathrm{ID}_{d}(A,M,b)_{h}$ and define the process $\{Y_{t}:t\geq 0\}$ by the stochastic differential equation

[TABLE]

where $\lambda>0$ is a parameter. This has a strong solution of the form

[TABLE]

In this case $Y=\{Y_{t}:t\geq 0\}$ is called an OU-process with parameter $\lambda$ and $Z$ is called the background driving Lévy process (BDLP). The process $Y$ is a Markov process and, so long as

[TABLE]

it has a limiting distribution. This distribution is necessarily selfdecomposable. Further, every selfdecomposable distribution is the limiting distribution of some OU-process. For details see [29] or [26].

Since $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ is selfdecomposable, it is the limiting distribution of some OU-process, which we call a TSOU-process. We now characterize the BDLP and the transition law of this process.

Theorem 1**.**

Consider the distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ and denote its Lévy measure by $L_{\alpha}$. Let $Y=\{Y_{t}:t\geq 0\}$ be an OU-process with parameter $\lambda>0$ and BDLP $Z=\{Z_{t}:t\geq 0\}$. Assume that $Z_{1}\sim\mathrm{ID}_{d}(0,\lambda M,\lambda c_{\alpha})_{h_{\alpha}}$ with

[TABLE]

and

[TABLE]

where $\rho_{\xi}(u)=-u^{-\alpha}q(\xi,u)$. In this case, $Y$ is a Markov process with temporally homogenous transition function $P_{t}(y,\mathrm{d}x)$ having characteristic function $\int_{\mathbb{R}^{d}}e^{i\langle x,z\rangle}P_{t}(y,\mathrm{d}x)=\exp\left\{C_{t}(y,z)\right\}$, where

[TABLE]

and

[TABLE]

Further, for any $y\in\mathbb{R}^{d}$, $P_{t}(y,\cdot)\stackrel{{\scriptstyle w}}{{\rightarrow}}\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ as $t\to\infty$.

In the above $\rho_{\xi}$ is a nondecreasing function and the integral with respect to $\mathrm{d}\rho_{\xi}(u)$ is the Stieltjes integral. We now give an explicit representation for an increment of a TSOU-process, under an additional assumption.

Theorem 2**.**

Let $\{Y_{t}:t\geq 0\}$ be a TSOU-process with parameter $\lambda>0$ and limiting distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$. Fix $t>0$ and let

[TABLE]

If $K=0$, then, given $Y_{s}=y$, we have

[TABLE]

and if $0<K<\infty$, then, given $Y_{s}=y$, we have

[TABLE]

In the above

[TABLE]

*and $N,X_{0},X_{1},X_{2},X_{3},\dots$ are independent random variables with:

1. $X_{0}\sim\mathrm{TS}_{\alpha,d}(\sigma_{0},q,0)$ where $\sigma_{0}(\mathrm{d}\xi)=(1-e^{-\alpha\lambda t})\sigma(\mathrm{d}\xi)$,

2. $X_{1},X_{2},X_{3},\dots$ are iid random variables such that $X_{i}\stackrel{{\scriptstyle d}}{{=}}\xi W$, where $\xi\in\mathbb{S}^{d-1}$, $W>0$, and the joint distribution of $\xi$ and $W$ is*

[TABLE]

3. $N$ has a Poisson distribution with mean $Ke^{-\alpha\lambda t}$.

Remark 2**.**

1. If $q(\xi,u)\equiv 1$, then $\mathrm{TS}_{\alpha,d}(\sigma,q,b)=\mathrm{S}_{\alpha,d}(\sigma,b)$. In this case $K=0$ and hence the transition law for an OU-process with an $\alpha$-stable limiting distribution can be represented as in (6). 2. If $\alpha\in(0,1)$ then $h_{\alpha}\equiv 0$ and hence $\psi=0$.

Theorem 2 provides a simple recipe for simulating an increment from a TSOU-process when $K<\infty$. The main ingredients are the ability to simulate from three distributions: the Poisson distribution, the $\mathrm{TS}_{\alpha,d}(\sigma_{0},q,0)$ distribution, and distribution $H$. The problem of simulating from a Poisson distribution is well-studied, see e.g. [1]. Under mild conditions, a rejection sampling technique for simulating from $\mathrm{TS}_{\alpha,d}(\sigma_{0},q,0)$ is given in [16]. To simulate from $H$ we first define the quantities

[TABLE]

the probability measure on $\mathbb{S}^{d-1}$ given by

[TABLE]

and the family of pdfs on $(0,\infty)$ given by

[TABLE]

It is readily checked that

[TABLE]

Hence, one can simulate $(\xi,W)$ from $H$ by first simulating $\xi$ from $\sigma_{1}$ and then $W$ from the distribution with pdf $f_{\xi}$. Thus, the problem reduces to that of simulating from $\sigma_{1}$ and from the pdf $f_{\xi}$.

We begin by discussing simulation from $\sigma_{1}$. One of the most important situations is when the support of $\sigma_{1}$ is finite. This always holds, in particular, when the dimension $d=1$. In this case, the problem reduces to simulating from a multinomial distribution. Another important situation is when $\sigma_{1}$ follows a uniform distribution. This problem is well-studied, see e.g. [31]. For a discussion of simulation from a variety of other standard distributions on $\mathbb{S}^{d-1}$, see the monograph [20]. While no method works in general, one can often set up an approximate simulation method by first approximating $\sigma_{1}$ by a distribution with a finite support, see Lemma 1 in [9].

We now discuss simulation from $f_{\xi}$ in two important situations with $K<\infty$. The first, which we call “hard truncation,” covers a useful but fairly specific situation, while the second, which we call “Class F,” is very general.

4.1 Hard Truncation

Consider the distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$, where $q(\xi,u)=1_{[0\leq u\leq\gamma_{\xi}]}$. Assume that $\gamma_{\xi}\geq 0$ for each $\xi\in\mathbb{S}^{d-1}$, that the function $\xi\mapsto\gamma_{\xi}$ is measurable, and that $\int_{\mathbb{S}^{d-1}}\gamma_{\xi}^{-\alpha}\sigma(\mathrm{d}\xi)<\infty.$ Such distributions appear in certain limit theorems, see [10] and [15]. In this case, it is readily checked that $K<\infty$ and that the pdf $f_{\xi}$, as given in (8), can be written as

[TABLE]

The corresponding cumulative distribution function (cdf) can be written as

[TABLE]

and, for $y\in(0,1)$,

[TABLE]

Hence, we can use the inverse transform method to simulate $W$ from $f_{\xi}$ by first simulating $U\sim U(0,1)$ and then taking $W=F_{\xi}^{-1}(U)$.

4.2 Class F

Consider the distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ and assume that there exists an $\epsilon>0$ and Borel functions $p(\cdot),M_{\bullet}:\mathbb{S}^{d-1}\mapsto[0,\infty)$ with

[TABLE]

such that for any $u,v\in(0,\epsilon)$ and $\sigma$-a.e. $\xi$

[TABLE]

In this case, we say that $q$ belongs to Class F and write $q\in\mathrm{CF}_{\alpha,d}(\epsilon,M_{\bullet},p(\cdot),\sigma)$. For simplicity, we will sometimes also say that the distribution $\mu=\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ belongs to Class F and write $\mu\in\mathrm{CF}_{\alpha,d}(\epsilon,M_{\bullet},p(\cdot),\sigma)$. It is easily verified that, in this case, $K<\infty$.

It is not always easy to check if a given tempering function belongs to Class F. A sufficient condition is that $q(\xi,u)=q_{1}(\xi,u^{p(\xi)})$, where, for $\sigma$-a.e. $\xi$, the function $q_{1}(\xi,\cdot)$ is uniformly Lipschitz continuous on a neighborhood of zero. By this, we mean that there are constants $M^{\prime},\epsilon^{\prime}>0$ such that for any $u,v\in(0,\epsilon^{\prime})$ and $\sigma$-a.e. $\xi\in\mathbb{S}^{d-1}$

[TABLE]

In this case $q\in\mathrm{CF}_{\alpha,d}(\epsilon,M_{\bullet},p(\cdot),\sigma)$ with $\epsilon=\epsilon^{\prime}$ and $M_{\bullet}\equiv M^{\prime}$. In particular, by the mean value theorem, (9) holds whenever $\left|\frac{\partial}{\partial u}q_{1}(\xi,u)\right|\leq M^{\prime}$ for all $u\in(0,\epsilon^{\prime})$ and $\sigma$-a.e. $\xi\in\mathbb{S}^{d-1}$.

We now develop a rejection sampling approach for simulating from $f_{\xi}$ in this case. The approach is based on the following result.

Proposition 1**.**

If $q\in\mathrm{CF}_{\alpha,d}(\epsilon,M_{\bullet},p(\cdot),\sigma)$ and $f_{\xi}$ is as in (8), then

[TABLE]

where $g_{1}$ is the pdf of $\mathrm{MLL}(\alpha,p(\xi),\delta_{0})$, $\delta_{0}=\epsilon e^{-\lambda t}$, and

[TABLE]

Let $\varphi_{1}(u)=f_{\xi}(u)/(V_{1}g_{1}(u))$ and note that

[TABLE]

With this notation, we get the following rejection sampling algorithm for simulating from $f_{\xi}$ for a fixed $\xi$.

**Algorithm 1.

**1. Independently simulate $U\sim U(0,1)$ and $Y\sim\mathrm{MLL}(\alpha,p(\xi),\delta_{0})$.

2. If $U\leq\varphi_{1}(Y)$ return $Y$, otherwise go back to step 1.

From general facts about rejection sampling algorithms, we know that, on a given iteration, the probability of accepting the observation is $1/V_{1}$ and the number of iterations until we accept follows a geometric distribution with mean $V_{1}$. Thus, the algorithm is more efficient when $V_{1}$ is small.

5 $p$-Tempered $\alpha$-Stable Distributions

In the previous sections, we considered tempered stable distributions with very general tempering functions. However, it is often convenient to work with families of tempering functions that have additional structure. One such family, which is commonly used, corresponds to the case where

[TABLE]

Here $p>0$ and $\bar{Q}=\{Q_{\xi}:\xi\in\mathbb{S}^{d-1}\}$ is a measurable family of probability measures on $(0,\infty)$. For fixed $p>0$ and $\alpha\in(0,2)$, we refer to the class of tempered stable distributions with such tempering functions as $p$-tempered $\alpha$-stable. For $p=1$, these models were introduced in [27], for $p=2$ they were introduced in [6], and the general case was introduced in [13]. See also the recent monograph [14] for an overview.

Now, consider the distribution $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$, where $q$ is as in (10) and define the Borel measures

[TABLE]

and

[TABLE]

When we know $R$, we can calculate $Q$ by using the formula

[TABLE]

It can be shown that, in this case, the Lévy measure as given by (1) can be written as

[TABLE]

and the measure $\sigma$ can be written as

[TABLE]

see Chapter 3 in [14]. Further, for fixed $\alpha\in(0,2)$ and $p>0$, the measure $R$ uniquely determines both $q$ and $\sigma$. For this reason, we will generally write $\mathrm{TS}^{p}_{\alpha,d}(R,b)$ instead of $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$. We call $R$ the Rosiński measure of the distribution. A Borel measure $R$ on $\mathbb{R}^{d}$ is the Rosiński measure of some $p$-tempered $\alpha$-stable distribution if and only if

[TABLE]

Remark 3**.**

One can define a more general class of distributions with Lévy measures of the form (14), but where $R$ does not satisfy (16). In [14] distributions that satisfy (16) are called proper $p$-tempered $\alpha$-stable.

We now characterize when we can use the representation given in Theorem 2 and when the distributions belong to class F.

Proposition 2**.**

*Let $\mu=\mathrm{TS}^{p}_{\alpha,d}(R,b)$ and let $K$ be as in (5).

1. If $p\leq\alpha$ or $R$ is an infinite measure, then $K=\infty$.

2. If $p>\alpha$ and $R$ is a finite measure, then*

[TABLE]

and the result of Theorem 2 holds. If, in addition,

[TABLE]

then $\mu\in\mathrm{CF}_{\alpha,d}(\epsilon,M_{\bullet},p,\sigma)$ with $M_{\xi}=\int_{(0,\infty)}sQ_{\xi}(\mathrm{d}s)$ and any $\epsilon>0$.

Remark 4**.**

For more on the finiteness of $M_{\xi}$, as given above, see Lemma 1 in Section 7 below. When $p\leq\alpha$, we have $K=\infty$ and the result of Theorem 2 does not hold. In this case a different framework is needed. Building on the work of [22] for certain types of Tweedie distributions, we will develop such a framework in a future work.

We now turn to the problem of simulating increments from the transition law of a TSOU-process with limiting measure $\mu=\mathrm{TS}^{p}_{\alpha,d}(R,b)$. For the remainder of this section assume that $p>\alpha$, that $0<R(\mathbb{R}^{d})<\infty$, and that (17) holds. Thus $\mu$ belongs to Class F and we can use Theorem 2. This requires us to simulate from the distributions $\mathrm{TS}_{\alpha,d}(\sigma_{0},q,0)$ and $H(\mathrm{d}\xi,\mathrm{d}u)$. It is easy to see that $\mathrm{TS}_{\alpha,d}(\sigma_{0},q,0)=\mathrm{TS}^{p}_{\alpha,d}(R_{0},b)$, where $R_{0}(\mathrm{d}x)=(1-e^{-\alpha\lambda t})R(\mathrm{d}x)$. So long as the support of $\sigma$ contains $d$ linearly independent vectors, we can use Algorithm 1 in [16] to simulate from this distribution, see Proposition 2 in that paper. To simulate from $H(\mathrm{d}\xi,\mathrm{d}u)$ we must simulate from $\sigma_{1}$ and $f_{\xi}$. We now derive more explicit formulas for this case.

Proposition 3**.**

If $\mu=\mathrm{TS}^{p}_{\alpha,d}(R,b)$ has $\bar{Q}=\{Q_{\xi}:\xi\in\mathbb{S}^{d-1}\}$, then

[TABLE]

and

[TABLE]

In light of Proposition 2, we can simulate from $f_{\xi}$ by using Algorithm 1 with $\delta_{0}=\epsilon e^{-\lambda t}=1$, $Y\sim\mathrm{MLL}(\alpha,p,1)$, and

[TABLE]

In this case,

[TABLE]

In some cases we can improve on Algorithm 1. An issue with the MLL distribution is that it has heavy tails, which can lead to many rejections when the tails of $f_{\xi}$ are lighter. However, when the support of $Q_{\xi}$ is lower bounded, we can use the pdf of the generalized gamma distribution, which has lighter tails. The method is based on the following result.

Proposition 4**.**

Let $\zeta=\sup\{c>0:Q_{\xi}((0,c))=0\}$. If $\zeta>0$ and $\int_{(0,\infty)}sQ_{\xi}(\mathrm{d}s)<\infty$, then

[TABLE]

where $g_{2}$ is the pdf of a $\mathrm{GGa}(p-\alpha,p,\zeta)$ distribution and

[TABLE]

Now, let $\varphi_{2}(u)=f_{\xi}(u)/(V_{2}g_{2}(u))$ and note that

[TABLE]

With this notation, when $\zeta>0$, we get the following rejection sampling algorithm for simulating from $f_{\xi}$ for a fixed $\xi\in\mathbb{S}^{d-1}$.

**Algorithm 2.

**1. Independently simulate $U\sim U(0,1)$ and $Y\sim\mathrm{GGa}(p-\alpha,p,\zeta)$.

2. If $U\leq\varphi_{2}(Y)$ return $Y$, otherwise go back to step 1.

Algorithm 2 works better than Algorithm 1 when $V_{2}<V_{1}$. Since

[TABLE]

Algorithm 2 is always better when $\zeta>\left(\Gamma(2-\alpha/p)\frac{\alpha}{p}\right)^{\frac{1}{1-\alpha/p}}$.

Example. In [21] a version of Algorithm 2 was derived for one-dimensional TSOU-processes with certain types of Tweedie limiting distributions. Specifically, the distributions considered are of the form $\mathrm{TS}^{p}_{\alpha,d}(R,0)$ where $d=1$, $p=1$, $\alpha\in(0,1)$, and $R(\mathrm{d}x)=a\zeta^{\alpha}\delta_{1/\zeta}(\mathrm{d}x)$, where $a,\zeta>0$. These correspond to $\sigma(\mathrm{d}\xi)=a\delta_{1}(\mathrm{d}\xi)$ and

[TABLE]

where $Q_{1}(\mathrm{d}r)=\delta_{\zeta}(\mathrm{d}r)$. In this case, the trial distribution is $\mathrm{GGa}(1-\alpha,1,\zeta)$ and

[TABLE]

It follows that our Algorithm 2 reduces to Algorithm 3 in [21]. We note that, in that paper, there appears to be a typo in the formula for $\varphi_{2}$, which they denote by $g_{2,\Delta}$. It should be as given above.

6 Simulation Study

In this section we perform a small-scale simulation study to see how well our methodology works in practice. We focus on a family of tempered stable distributions for which the transition law had not been previously derived.

6.1 Power Tempered Stable Distributions

A distribution is said to be power tempered stable if it is $\mathrm{TS}^{p}_{\alpha,d}(R,b)$ with $d=1$, $p=1$, $\alpha\in(0,1)$, and

[TABLE]

where

[TABLE]

and $\ell,c>0$ are parameters. We denote such distributions by $\mathrm{PT}_{\alpha}(\ell,c)$. These models were introduced in [14] as a class of tempered stable distributions with a finite mean, but still fairly heavy tails. In fact, if $Y\sim\mathrm{PT}_{\alpha}(\ell,c)$, then, for $\beta\geq 0$,

[TABLE]

Thus, $\ell$ controls how heavy the tails are. Methods for evaluating the pdfs and related quantities for these distribution are available in the SymTS package [17] for the statistical software R. We now give some useful facts.

Proposition 5**.**

For a $\mathrm{PT}_{\alpha}(\ell,c)$ distribution with $\alpha\in(0,1)$ and $\ell,c>0$ we have

[TABLE]

Further,

[TABLE]

[TABLE]

and for $\xi\in\mathbb{S}^{0}=\{-1,1\}$

[TABLE]

where

[TABLE]

For these distributions, an exact simulation technique, based on rejection sampling, was developed in [16]. Specifically, to simulate from $\mathrm{PT}_{\alpha}(\ell,c)$ we can use Algorithm 1 in [16] with, in the notation of that paper,

[TABLE]

see Proposition 2 in [16]. Detailed simulations can be found in that paper.

Here, we are interested in simulating increments from a TSOU-process with limiting distribution $\mathrm{PT}_{\alpha}(\ell,c)$. Proposition 2 implies that this distribution belongs to Class F and that we can use the representation given in Theorem 2. To use this representation, we need to simulate from $\mathrm{TS}^{1}_{\alpha,1}(R_{0},0)$ and $H(\mathrm{d}\xi,\mathrm{d}u)$. It is readily checked that, in this case, $\mathrm{TS}^{1}_{\alpha,1}(R_{0},0)=\mathrm{PT}_{\alpha}(\ell,\left(1-e^{-\alpha\lambda t}\right)c)$ and, thus, that we can use the methodology described above to simulate this component. Further, to implement our methodology, we can use the following.

Proposition 6**.**

For a $\mathrm{PT}_{\alpha}(\ell,c)$ distribution with $\alpha\in(0,1)$ and $\ell,c>0$ we have

[TABLE]

and

[TABLE]

Further, for $\xi\in\mathbb{S}^{0}=\{-1,1\}$ and $u>0$ we have

[TABLE]

and for $u>0$ we have

[TABLE]

Note that none of the quantities discussed in Proposition 6 depend on the parameter $c$. Further, both $f_{\xi}$ and $\varphi_{1}$ do not depend on $\xi$. Now, we can simulate from $H$ as follows. First, we simulate $\xi$ from $\sigma_{1}$ by taking

[TABLE]

Next, we simulate $W$ from $f_{\xi}$. We do this using Algorithm 1, in which we take $Y\sim\mathrm{MLL}(\alpha,1,1)$ and $\varphi_{1}$ as given above.

6.2 Univariate Simulation Results

We now present the results of our simulation study, where we focused on simulating increments from a TSOU-process with parameter $\lambda$ and limiting distribution $\mathrm{PT}_{\alpha}(\ell,c)$. For simplicity, throughout this section, we fix the parameters $\lambda=1$ and $c=1$. First, we consider simulation from $H$. Here we choose the parameter values $\alpha=0.55$, $\ell=1$, and $t=0.1$. Figure 1 gives a plot of the pdf $f_{1}$ (solid line) overlaid with a plot of the pdf of $\mathrm{MLL}(0.55,1,1)$ (dashed line). We then simulate $n=10^{7}$ observations from $\mathrm{MLL}(0.55,1,1)$, and, using Algorithm 1, we decide which observations should be rejected. We numerically calculated that $V_{1}=12.8837$ and thus we expect to get $n/V_{1}=776175.4$ observations. We actually obtained $776528$ observations. We then assigned to each of these a positive or a negative sign based on (22). A plot of the kernel density estimate (KDE) of the pdf of the resulting observations is given in Figure 2. This is overlaid (dashed line) with the pdf of the true density of $X=\xi W$.

We now turn to the simulation of TSOU-processes. Figure 3 presents four simulated paths of the TSOU-process with limiting distribution $\mathrm{PT}_{\alpha}(\ell,1)$ for different choices of the parameters. The paths go from time $T=0$ to time $T=100$ in increments of $t=0.1$. Thus, each path consists of $1000$ increments. For simplicity, we start each path at [math], which is the mean of the limiting distribution. To check that we really get the correct limiting distribution, we simulated $3000$ paths of the process for two choices of the parameters. We then took the final observation (at time $T=100$) from each of these. These are independent random variables from $\mathrm{PT}_{\alpha}(\ell,1)$. In Figure 4 we plot the KDE of their densities (solid line) with the true pdf of $\mathrm{PT}_{\alpha}(\ell,1)$ overlaid (dashed line).

6.3 Multivariate Simulation Results

In this section we give simulations for the multivariate case. We begin by introducing a multivariate extension of the power tempered stable distribution with a finite measure $\sigma$. Let $d\geq 2$ be the dimension, let $s_{1},s_{2},\dots,s_{k}\in\mathbb{S}^{d-1}$ with $s_{i}\neq s_{j}$ and $s_{i}\neq-s_{j}$ for $i\neq j$, and let

[TABLE]

where $X_{1},X_{2},\dots,X_{k}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{PT}_{\alpha}(\ell,c)$. It is not difficult to show that $X\sim\mathrm{TS}^{1}_{\alpha,d}(R,0)$ with

[TABLE]

where

[TABLE]

and $S=\{\pm s_{1},\pm s_{2},\dots,\pm s_{k}\}$. We denote this distribution by $\mathrm{PT}_{\alpha,d,k}(\ell,c,S)$. Further, we can show that, in this case,

[TABLE]

and that, for $\xi\in S$, $Q_{\xi}$ is given by (18) and $f_{\xi}$ is given by (19).

We are interested in simulating a multivariate TSOU-process with parameter $\lambda>0$ and limiting distribution $\mathrm{PT}_{\alpha,d,k}(\ell,c,S)$. Proposition 2 implies that this distribution belongs to Class F and that we can use the representation of the transition law given in (7). To use this representation, we need to simulate from $\mathrm{TS}^{1}_{\alpha,d}(R_{0},0)$ and $H(\mathrm{d}\xi,\mathrm{d}u)$. It is readily checked that $\mathrm{TS}^{1}_{\alpha,d}(R_{0},0)=\mathrm{PT}_{\alpha,d,k}(\ell,\left(1-e^{-\alpha\lambda t}\right)c,S)$ and (since we already know how to simulate from $\mathrm{PT}_{\alpha}(\ell,\left(1-e^{-\alpha\lambda t}\right)c)$) we can use (23) to simulate from this distribution. To simulate from $H$, we first simulate $\xi$ uniformly from the finite set $S$, then we simulate $W$ from $f_{\xi}$. For simplicity, in our simulations, we take $d=2$, $k=3$, and we fix the parameters $\lambda=1$, $c=1$, $\ell=1$, and $\alpha=0.55$. We simulated two TSOU-processes with different choices of $S$. As in the univariate case, the paths go from time $T=0$ to time $T=100$ in increments of $t=0.1$. For simplicity we start each path at $\left({0\atop 0}\right)$, which is the mean of the limiting distribution. The plots of the $x$ and $y$ coordinates of the two processes are given in Figure 5.

7 Proofs

We begin with several lemmas.

Lemma 1**.**

Let $Q$ be given by (11) and $R$ by (12). For any $\gamma>0$, we have

[TABLE]

if and only if

[TABLE]

Proof.

The fact that

[TABLE]

gives the result. ∎

Lemma 2**.**

If $0\leq a<b<\infty$ then

[TABLE]

Proof.

The result follows from the fact that $e^{-a}-e^{-b}=\int_{a}^{b}e^{-u}\mathrm{d}u$. ∎

Lemma 3**.**

For any $\kappa>0$ and $0<\alpha<p$

[TABLE]

Proof.

We have

[TABLE]

where the second line follows by (3.38) in [14]. ∎

Proof of Theorem 1.

From Lemma 17.1 in [29] it follows that $Y$ is a Markov process with temporally homogenous transition function $P_{t}(y,\mathrm{d}x)$ on $\mathbb{R}^{d}$ having a characteristic function given by

[TABLE]

We begin with the case $\alpha\neq 1$. In this case it is easily checked that, for $a\in\mathbb{R}$, we have $\psi_{\alpha}(az,x)=\psi_{\alpha}(z,ax)$. It follows that

[TABLE]

where the third line follows by the substitution $u=e^{-\lambda s}r$ and the fifth line by the fact that

[TABLE]

When $\alpha=1$ and $a\in\mathbb{R}$, we have

[TABLE]

From here, proceeding as before, we get

[TABLE]

Now, note that the third term equals

[TABLE]

We now show the convergence to $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$. It suffices to show that the characteristic function of the transition law approaches the characteristic function of $\mathrm{TS}_{\alpha,d}(\sigma,q,b)$ as $t\to\infty$. First, since $|q(\xi,u)|\leq 1$ for each $\xi\in\mathbb{S}^{d-1}$ and $u>0$, we have

[TABLE]

where the final equality follows from the fact that the integral is finite, which itself follows by standard facts about the function $\psi_{\alpha}$, see (26.4) in [8]. Now, note that $d_{\alpha}\to b$ as $t\to\infty$. This is immediate for $\alpha\neq 1$ and holds for $\alpha=1$ by dominated convergence. From here it follows that

[TABLE]

which is the characteristic function of $\mathrm{TS}_{\alpha}(\sigma,q,b)$ as required. ∎

Proof of Theorem 2.

We only prove (7) as the proof of (6) is similar. It suffices to check that the characteristic function of the random variable on the right side of (7) is given by $\exp\left\{C_{t}(y,z)\right\}$, where $C_{t}(y,z)$ is as in Theorem 1. This follows immediately from the fact that the Lévy measure of $\mathrm{TS}_{\alpha,d}(\sigma,q_{0},0)$ is $(1-e^{-\alpha\lambda t})L_{\alpha}(\mathrm{d}x)$ and $\sum_{j=1}^{N}X_{j}$ is a compound Poisson random variable, whose characteristic function is given in e.g. Proposition 3.4 of [11]. ∎

Proof of Proposition 1.

If $u\in(0,\epsilon e^{-\lambda t})$, then

[TABLE]

On the other hand, if $u\geq\epsilon e^{-\lambda t}$, then the fact that $q(\xi,u)\leq 1$ implies that

[TABLE]

and the result follows. ∎

Proof of Proposition 2.

First, note that by change of variables

[TABLE]

Thus, if $R$ is an infinite measure, then $K=\infty$. Next, if $\alpha\geq p$, then Lemma 2 implies that

[TABLE]

and we again have $K=\infty$. Henceforth, assume that $R$ is a finite measure and $p>\alpha$. In this case, Lemma 3 implies that

[TABLE]

which gives the form of $K$. Now assume that (17) holds. We again use Lemma 2, which gives, for any $u,v>0$,

[TABLE]

From here the result follows by Lemma 1. ∎

Proof of Proposition 3..

Note that

[TABLE]

From here the formula for $\kappa_{\xi}$ follows from Lemma 3. The rest of the formulas are immediate. ∎

Proof of Proposition 4.

Lemma 2 implies that

[TABLE]

which gives the result. ∎

Proof of Propositions 5 and 6.

The results in Proposition 5 follow immediately from (13), (15), and (11). We now turn to Propositions 6. Let $\phi(x)=(\alpha+\ell+1)^{-1}\left(\frac{x}{x+1}\right)^{\alpha+\ell+1}$ for $x>0$ and note that $\phi^{\prime}(x)=x^{\alpha+\ell}(x+1)^{-\alpha-\ell-2}$. It follows that

[TABLE]

Using this fact we can easily put the formulas into the required form. ∎