# Efficient Simulation of $p$-Tempered $\alpha$-Stable OU Processes

Michael<EMAIL_ADDRESS>
University of North Carolina Charlotte
Piergiacomo<EMAIL_ADDRESS>333The views, opinions,
positions or strategies expressed in this article are those of the authors and
do not necessarily represent the views, opinions, positions or strategies of,
and should not be attributed to E.ON SE.
Quantitative Risk Management, E.ON SE
Brüsseler Platz 1, 45131 Essen, Germany and
Department of Mathematics and Statistics,
University of Helsinki P.O. Box 68 FI-00014 Finland

###### Abstract

We develop efficient methods for simulating processes of Ornstein-Uhlenbeck
type related to the class of $p$-tempered $\alpha$-stable
($\mathrm{TS}^{p}_{\alpha}$) distributions. Our results hold for both the
univariate and multivariate cases and we consider both the case where the
$\mathrm{TS}^{p}_{\alpha}$ distribution is the stationary law and where it is
the distribution of the background driving Lévy process (BDLP). In the latter
case, we also derive an explicit representation for the transition law as this
was previous known only in certain special cases and only for $p=1$ and
$\alpha\in[0,1)$. Simulation results suggest that our methods work well in
practice.

##  1 Introduction

Tempered stable processes of Ornstein Uhlenbeck type have been the subject of
much research in recent years. They combine two important directions. First,
they are examples of non-Gaussian processes of Ornstein Uhlenbeck type (OU
processes), which have a more intricate dependence structure than the more
commonly used Lévy processes and, in particular, they are mean reverting. In
financial applications this makes them natural models for various quantities
including stochastic volatility, stochastic interest rates, and commodity
prices, see, e.g., [1], [2], [26] and [5]. Second, they are based on the class
of tempered stable (TS) distributions, which has been gaining in prominence
over the last decade. These distributions can approximate the more common
Gaussian and stable distributions, but their tail behavior is more realistic,
which makes them useful for a variety of application areas. We are
particularly motivated by their use in the modeling of financial returns, see
[15], [31], and the references therein. TS distributions were first formalized
in the classic paper [24]. Since then, they have been extended in several
directions in [25], [20], and [8], see also the monograph [9].

There are two types of OU processes related to TS distributions. The first,
denoted TSOU processes, correspond to the case where the stationary
distribution is TS. The second, denoted OUTS processes, correspond to the case
where the background driving Lévy process (BDLP) has a TS distribution. The
study of the transition laws of TSOU processes has been primarily focused on
the fairly simple univariate class of so-called classical tempered stable
(CTS) distributions, see [35], [32], [19], [22], [28], and the references
therein. More detailed results for the special cases of gamma and inverse
Gaussian distributions are given in [34], [21], and [27]. Extensions to
general classes of TS distributions, including the multivariate case, are
given in [11] and [13]. Significantly less attention has been paid to the
class of OUTS processes. To the best of our knowledge, in this case, the
transition laws have only been studied in the case of gamma distributions in
[21] and [27], and in the case of CTS distributions with parameter
$\alpha\in[0,1)$ in [22] and [28]. Although, it should be noted that some
preliminary results about certain univariate TS distributions beyond CTS are
given in [22].

The purpose of the current paper is two-fold. First, we derive the transition
laws of OUTS processes for the class of $p$-tempered $\alpha$-stable
($\mathrm{TS}^{p}_{\alpha}$) distributions with any $\alpha\in(-\infty,2)$ and
$p>0$. These distributions form a large and flexible class of both univariate
and multivariate models and include CTS distributions as a special case. Our
results compliment those in [13] for the corresponding class of TSOU
processes. Second, we develop efficient methods for simulating from the
transition laws of both TSOU and OUTS processes based on
$\mathrm{TS}^{p}_{\alpha}$ distributions. These can then be used to simulate
the corresponding OU process on a finite grid. Our simulation methods extend
the ideas introduced in [28] for CTS distributions and work well in practice.
The main idea is based on showing that certain components of the transition
law can be represented as generalized gamma scale mixtures (GGSMs). For this
reason we develop the theory of such mixtures and give multiple approaches for
simulation.

The rest of the paper is organized as follows. In Section 2 we introduce GGSM
distributions and discuss several examples that are important for simulating
TSOU and OUTS processes. In Section 3 we recall the definition of
$\mathrm{TS}^{p}_{\alpha}$ distributions and give some properties. Then in
Section 4 we recall basic properties of OU processes and, in particular, we
give our results on the transition laws of OUTS processes. In Section 5 we
perform a series of numerical experiments to better understand the performance
of our simulation methods. Proofs are postponed to Section 6.

Before proceeding we introduce some notation. We write cdf, pdf, and pmf for
_cumulative distribution function_ , _probability density function_ , and
_probability mass function_ , respectively. For a distribution $F$ we write
$X\sim F$ to denote that $X$ is a random variable with distribution $F$ and we
write $X_{1},X_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}F$ to
denote that $X_{1},X_{2},\dots$ are independently and identically distributed
random variables with distribution $F$. For simplicity, instead of $F$, we
sometimes write the corresponding pdf or pmf. We write $U(0,1)$ to denote a
uniform distribution on $(0,1)$, $\delta_{a}$ to denote a point-mass at $a$,
and $1_{A}$ to denote the indicator function on $A$. We write $\vee$ and
$\wedge$ to denote the maximum and minimum, respectively, and we write
$\lfloor\cdot\rfloor$ to denote the floor function. Further, we use the
convention that $\sum_{n=1}^{0}$ is $0$. We write $\stackrel{{\scriptstyle
d}}{{=}}$ and $=:$ to denote equality in distribution and a defining equality,
respectively.

##  2 Generalized Gamma Scale Mixtures

In this section we introduce the class of generalized gamma scale mixture
(GGSM) distributions and discuss various properties and approaches for
simulation. We then show that the incomplete gamma (IGa) distribution, which
was introduced in [13] as an important component of the transition law of a
TSOU process, is a GGSM and use this fact to develop efficient simulation
techniques. Next, we introduce two new GGSM distributions, which play a
similar role in the study of OUTS processes. We note that the results in this
section may be of independent interest.

We begin by recalling that the generalized gamma distribution was introduced
in [30] and has a pdf given by

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

where $\gamma,p,\theta>0$ are parameters. We denote this distribution by
$\mathrm{GGa}(\gamma,p,\theta)$. When $p=1$ it reduces to the usual gamma
distribution, which we denote by $\mathrm{Ga}(\gamma,\theta)$. We can simulate
from $\mathrm{GGa}(\gamma,p,\theta)$ by using the fact that

$\displaystyle\mbox{if }X\sim\mathrm{Ga}(\gamma/p,1)\mbox{, then
}\left(X/\theta\right)^{1/p}\sim\mathrm{GGa}(\gamma,p,\theta),$ (1)

see [11] and [13] for more on this distribution.

A GGSM distribution has a pdf of the form

$\displaystyle
f(u)=\int_{0}^{\infty}g_{\gamma,p,\theta^{p}}(u)m(\theta)\mathrm{d}\theta,$
(2)

where $m$ is the pdf of the so-called mixing distribution, whose support is
contained in $[0,\infty)$. Note that, in (2), the parameter is $\theta^{p}$
and not just $\theta$. To simulate from a GGSM we can first simulate $Z\sim m$
and then, given $Z$, simulate $X\sim\mathrm{GGa}(\gamma,p,Z^{p})$.
Equivalently, using (1), we get the following Algorithm.

Algorithm GGSM1.
Step 1. Independently simulate $Y\sim\mathrm{Ga}(\gamma/p,1)$ and $Z\sim m$.
Step 2. Return $Y^{1/p}/Z$.

In practice, it is not always easy to simulate from $m$. Under general
assumptions, which always hold in the situations that are of interest to us,
we can set up a rejection sampling approach by using the following result,
which follows immediately from (2).

######  Lemma 2.1.

If the support of $m$ is lower bounded by some $a>0$ and if
$\int_{a}^{\infty}\theta^{\gamma}m(\theta)\mathrm{d}\theta<\infty$, then

$f(u)\leq Vg_{\gamma,p,a^{p}}(u),$

where
$V=a^{-\gamma}\int_{a}^{\infty}\theta^{\gamma}m(\theta)\mathrm{d}\theta$.

Let

$\varphi(u)=\frac{\int_{a}^{\infty}e^{-u(\theta^{p}-a^{p})}\theta^{\gamma}m(\theta)\mathrm{d}\theta}{\int_{a}^{\infty}\theta^{\gamma}m(\theta)\mathrm{d}\theta}.$

Algorithm GGSM2.
Step 1. Independently simulate $U\sim U(0,1)$ and
$Y\sim\mathrm{Ga}(\gamma/p,1)$.
Step 2. If $U\leq\varphi(Y/a^{p})$ return $Y^{1/p}/a$, otherwise go back to
step 1.

From standard results about rejection sampling algorithms, the probability of
rejection on a given iteration is given by $1/V$.

###  2.1 Incomplete Gamma Distribution

For $u>0$ and $\gamma>0$, let

$G_{\gamma}(u)=\frac{1}{\Gamma(\gamma)}\int_{0}^{u}x^{\gamma-1}e^{-x}\mathrm{d}x$

be the scaled lower incomplete gamma function. It is the cdf of the
$\mathrm{Ga}(\gamma,1)$ distribution. The incomplete gamma (IGa) distribution
has a pdf given by

$f_{\beta,\gamma,p,\eta}(u)=\frac{1}{K_{\beta,\gamma,p,\eta}}G_{\gamma}(u^{p}(\eta-1))e^{-u^{p}}u^{-1-\beta},\qquad
u>0,$

where $\gamma>0$, $p>0$, $\eta>1$, $\beta\in(-\infty,p\gamma)$, and
$K_{\beta,\gamma,p,\eta}$ is a normalizing constant given by

$K_{\beta,\gamma,p,\eta}=\frac{\Gamma(\gamma-\beta/p)}{p\Gamma(\gamma)}K^{*}_{\beta,\gamma,p,\eta}$

with

$\displaystyle
K^{*}_{\beta,\gamma,p,\eta}=\int_{1/\eta}^{1}(1-x)^{\gamma-1}x^{-\beta/p-1}\mathrm{d}x=\int_{1}^{\eta}(x-1)^{\gamma-1}x^{\beta/p-\gamma}\mathrm{d}x.$

We denote this distribution by $\mathrm{IGa}(\beta,\gamma,p,\eta)$. In [13] it
was shown that if $W\sim\mathrm{IGa}(\beta,\gamma,p,\eta)$ and
$\xi>\beta-p\gamma$, then

$\displaystyle\mathrm{E}\left[W^{\xi}\right]=\frac{K_{(\beta-\xi),\gamma,p,\eta}}{K_{\beta,\gamma,p,\eta}}.$
(3)

Further, Proposition 1 in that paper shows that

$\displaystyle
K^{*}_{\beta,\gamma,p,\eta}\sim\frac{(\eta-1)^{\gamma}}{\gamma}\mbox{ as
}\eta\downarrow 1.$ (4)

We now show that the IGa distribution is a GGSM. When $p=1$ and $\gamma=1$,
this was already observed in [28].

######  Lemma 2.2.

We have

$\displaystyle
f_{\beta,\gamma,p,\eta}(u)=\int_{1}^{\eta^{1/p}}g_{(p\gamma-\beta),p,\theta^{p}}(u)m_{\beta,\gamma,p,\eta}(\theta)\mathrm{d}\theta,$

where

$\displaystyle
m_{\beta,\gamma,p,\eta}(\theta)=\frac{p}{K^{*}_{\beta,\gamma,p,\eta}}(\theta^{p}-1)^{\gamma-1}\theta^{p+\beta-p\gamma-1},\
\ 1<\theta<\eta^{1/p}.$ (5)

Since IGa is a GGSM, we can use Algorithms GGSM1 and GGSM2 to simulate from
it. It is readily checked that if $Z\sim m_{\beta/p,\gamma,1,\eta}$, then
$Z^{1/p}\sim m_{\beta,\gamma,p,\eta}$. This implies that Algorithm GGSM1
reduces to the following.

Algorithm IGa1. Simulation from $\mathrm{IGa}(\beta,\gamma,p,\eta)$.
Step 1. Independently simulate $Y\sim\mathrm{Ga}(\gamma-\beta/p,1)$ and $Z\sim
m_{\beta/p,\gamma,1,\eta}$.
Step 2. Return $(Y/Z)^{1/p}$.

Next, we note that for the IGa distribution Algorithm GGSM2 reduces to the
algorithm introduced in [13]. It can be stated as follows. Let

$\displaystyle\varphi_{1}(u)$ $\displaystyle=$
$\displaystyle\frac{\int_{1}^{\eta^{1/p}}e^{-u(\theta^{p}-1)}(\theta^{p}-1)^{\gamma-1}\theta^{p-1}\mathrm{d}\theta}{\int_{1}^{\eta^{1/p}}(\theta^{p}-1)^{\gamma-1}\theta^{p-1}\mathrm{d}\theta}=\frac{\Gamma(\gamma+1)}{(\eta-1)^{\gamma}}G_{\gamma}((\eta-1)u)u^{-\gamma}.$

Algorithm IGa2. Simulation from $\mathrm{IGa}(\beta,\gamma,p,\eta)$.
Step 1. Independently simulate $U\sim U(0,1)$ and
$Y\sim\mathrm{Ga}(\gamma-\beta/p,1)$.
Step 2. If $U\leq\varphi_{1}(Y)$ return $Y^{1/p}$, otherwise go back to step
1.

In [13] it is shown that, on a given iteration, the probability of acceptance
is $1/V_{1}$, where $V_{1}=\frac{(\eta-1)^{\gamma}}{\gamma
K^{*}_{\beta,\gamma,p,\eta}}$. From (4) it follows that $1/V_{1}\to 1$ as
$\eta\downarrow 1$. As we will see, when simulating TSOU processes, we
typically take $\eta$ close to $1$.

In order to use Algorithm IGa1, we need a way to simulate from
$m_{\beta,\gamma,1,\eta}$. We will provide several algorithms for doing this.
First, we introduce the pdf

$\ell_{\delta,\eta}(\theta)=\frac{\delta+1}{\eta^{\delta+1}-1}\theta^{\delta}\quad
1<\theta<\eta,$

where $\delta\in\mathbb{R}$ and $\eta>1$ are parameters. Here and throughout
we interpret $\frac{\delta+1}{\eta^{\delta+1}-1}$ by its limiting value of
$1/\ln\eta$ when $\delta=-1$. It is easy to check that we can simulate from
this distribution as follows.

Algorithm $\ell$1. Simulation from $\ell_{\delta,\eta}$.
Step 1. Simulate $U\sim U(0,1)$.
Step 2. If $\delta=-1$ return $\eta^{U}$. Otherwise, return
$\left[1+U\,\left(\eta^{\delta+1}-1\right)\right]^{\frac{1}{\delta+1}}$.

Note that, when $\gamma=1$, we have $m_{\beta,1,1,\eta}=\ell_{\beta-1,\eta}$,
which leads to the following algorithm.

Algorithm M0. Simulation from $m_{\beta,1,1,\eta}$.
Step 1. Simulate $U\sim U(0,1)$.
Step 2. If $\beta=0$ return $\eta^{U}$. Otherwise, return
$\left(U(\eta^{\beta}-1)+1\right)^{1/\beta}$.

We now turn to the case $\gamma>1$. It is readily checked that

$m_{\beta,\gamma,1,\eta}(\theta)\leq
V_{1}^{*}\ell_{\beta-\gamma,\eta}(\theta),$

where

$V_{1}^{*}=\frac{1}{K^{*}_{\beta,\gamma,1,\eta}}(\eta-1)^{\gamma-1}\frac{\eta^{\beta-\gamma+1}-1}{\beta-\gamma+1}.$

Letting

$\varphi_{1}^{*}(y)=\left(\frac{y-1}{\eta-1}\right)^{\gamma-1}$

leads to the following algorithm.

Algorithm M1. Simulation from $m_{\beta,\gamma,1,\eta}$ with $\gamma>1$.
Step 1. Independently simulate $U\sim U(0,1)$ and
$Y\sim\ell_{\beta-\gamma,\eta}$.
Step 2. If $U_{2}\leq\varphi_{1}^{*}(Y)$ return $Y$, otherwise go back to step
1.

On a given iteration, the probability of acceptance is $1/V^{*}_{1}$. From (4)
it follows that $1/V^{*}_{1}\to 1/\gamma$ as $\eta\downarrow 1$. Thus, when
$\eta$ is close to $1$ this method works better for smaller values of
$\gamma>1$.

######  Remark 1.

Alternatively, we can note that for $\gamma>1$ and $\beta\neq 0$

$m_{\beta,\gamma,1,\eta}(\theta)\leq V_{1,1}^{*}\ell_{\beta-1,\eta}(\theta),$

where

$V_{1,1}^{*}=\frac{1}{K^{*}_{\beta,\gamma,1,\eta}}\frac{\eta^{\beta}-1}{\beta}.$

This can be used to develop another rejection sampling method. This method may
work better than Algorithm M1 for some choices of the parameters. However, (4)
implies that $1/V_{1,1}^{*}\to 0$ as $\eta\downarrow 1$. As such it will not
work well for the situation of interest.

For the remaining methods we only consider the case where
$\gamma\in\\{1,2,3,\dots\\}$ is an integer. In this case the binomial theorem
gives

$m_{\beta,\gamma,1,\eta}(\theta)=\frac{1}{K^{*}_{\beta,\gamma,1,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\theta^{-1-k+\beta},\qquad 1<y<\eta.$ (6)

Integrating from $1$ to $y\in(1,\eta)$ we get the cdf

$\displaystyle
M_{\beta,\gamma,1,\eta}(y)=\frac{1}{K^{*}_{\beta,\gamma,1,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\frac{y^{\beta-k}-1}{\beta-k},$ (7)

where, in the case $k=\beta$, we interpret $\frac{y^{-k+\beta}-1}{\beta-k}$ by
its limiting value of $\ln y$. Let $M_{\beta,\gamma,1,\eta}^{-1}$ be the
inverse function of $M_{\beta,\gamma,1,\eta}$. This can be calculated
numerically, which leads to the following algorithm.

Algorithm M2. Simulation from $m_{\beta,\gamma,1,\eta}$ with
$\gamma\in\\{1,2,3,\dots\\}$.
Step 1. Simulate $U\sim U(0,1)$.
Step 2. Return $M_{\beta,\gamma,1,\eta}^{-1}(U)$.

When $\gamma=1$ the inverse function has a simple form and this algorithm
reduces to Algorithm M0. Our last algorithm is based on the methodology in
[3]. The idea is to use the positive terms in (6) to obtain the bound

$\displaystyle m_{\beta,\gamma,1,\eta}(\theta)$ $\displaystyle\leq$
$\displaystyle\frac{1}{K^{*}_{\beta,\gamma,1,\eta}}\sum_{k=0}^{\lfloor(\gamma-1)/2\rfloor}{\gamma-1\choose
2k}\theta^{-1-2k+\beta}$ $\displaystyle=$ $\displaystyle
V^{*}_{2}\sum_{k=0}^{\lfloor(\gamma-1)/2\rfloor}p_{\beta,\gamma,\eta}(k)\,\ell_{\beta-2k-1,\eta}(\theta)=:V_{2}^{*}\tilde{m}_{\beta,\gamma,1,\eta}(\theta),\quad
1<\theta<\eta,$

where

$V_{2}^{*}=\frac{\sum_{k=0}^{\lfloor(\gamma-1)/2\rfloor}H^{*}_{2k,\beta,\gamma,\eta}}{K^{*}_{\beta,\gamma,1,\eta}},$
$p_{\beta,\gamma,\eta}(k)=\frac{H^{*}_{2k,\beta,\gamma,\eta}}{\sum_{k=0}^{\lfloor(\gamma-1)/2\rfloor}H^{*}_{2k,\beta,\gamma,1,\eta}},\
\ \ k=0,\dots,\lfloor(\gamma-1)/2\rfloor,$

and

$H^{*}_{k,\beta,\gamma,\eta}={\gamma-1\choose
k}\frac{\eta^{\beta-k}-1}{\beta-k},\ \ \ k=0,\dots,\gamma-1.$

Using the binomial theorem, it can be checked that
$K^{*}_{\beta,\gamma,1,\eta}=\sum_{k=0}^{\gamma-1}(-1)^{k}H^{*}_{k,\beta,\gamma,\eta}$,
which guarantees that $V_{2}^{*}\geq 1$. Clearly,
$\sum_{k=1}^{\lfloor(\gamma-1)/2\rfloor}p_{\beta,\gamma,\eta}(k)=1$ and thus
$p_{\beta,\gamma,\eta}$ is a valid pmf. It follows that
$\tilde{m}_{\beta,\gamma,1,\eta}$ is a mixture distribution and we can
simulate from it as follows.

Algorithm BD. Simulation from $\tilde{m}_{\beta,\gamma,1,\eta}$.
Step 1. Simulate $S\sim p_{\beta,\gamma,\eta}$.
Step 2. Simulate $Y\sim\ell_{\beta-2S-1,\eta}$ and return $Y$.

Now, letting

$\varphi_{2}^{*}(y)=\frac{(\theta-1)^{\gamma-1}\theta^{1-\gamma}}{\sum_{k=0}^{\lfloor(\gamma-1)/2\rfloor}{\gamma-1\choose
2k}\theta^{-2k}}$

leads to the following algorithm.

Algorithm M3. Simulation from $m_{\beta,\gamma,1,\eta}$ with
$\gamma\in\\{2,3,\dots\\}$.
Step 1. Independently simulate $U\sim U(0,1)$ and
$Y\sim\tilde{m}_{\beta,\gamma,1,\eta}$.
Step 2. If $U\leq\varphi_{2}^{*}(Y)$ return $Y$, otherwise go back to step 1.

On a given iteration, the probability of acceptance is $1/V^{*}_{2}$. From
l’Hôpital’s Rule and (4) it follows that $1/V^{*}_{2}\to 0$ as $\eta\downarrow
1$. Thus, this method does not work well when $\eta$ is close to $1$. However,
it may work well in other cases.

###  2.2 Incomplete Beta Gamma Mixture Distribution

We now introduce a distribution, which is important for the simulation of OUTS
processes. To the best of our knowledge this distribution has not been studied
previously. It has a pdf of the form

$\displaystyle
f^{\sharp}_{\beta,\gamma,p,\eta}(v)=\frac{1}{C_{\beta,\gamma,p,\eta}}v^{p\gamma-\beta-1}\int_{1}^{\eta}\theta^{p\gamma-\beta-1}e^{-v^{p}\theta^{p}}\int_{1/\theta}^{1}(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u\mathrm{d}\theta,\
\ v>0,$

where $p>0$, $\gamma>0$, and $\beta\in(-\infty,p\gamma)$ are parameters and
$C_{\beta,\gamma,p,\eta}$ is a normalizing constant. We call this an
incomplete beta gamma mixture (IBGM) distribution and we denote it by
$\mathrm{IBGM}(\beta,\gamma,p,\eta)$. It is readily seen that

$\displaystyle C_{\beta,\gamma,p,\eta}$ $\displaystyle=$ $\displaystyle
p^{-1}\Gamma(\gamma-\beta/p)C^{*}_{\beta,\gamma,p,\eta},$

where

$\displaystyle
C^{*}_{\beta,\gamma,p,\eta}=\int_{1}^{\eta}\theta^{-1}\int_{1/\theta}^{1}(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u\mathrm{d}\theta=\int_{1/\eta}^{1}\ln\left(\eta
u\right)(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u.$

It can be checked that

$f^{\sharp}_{\beta,\gamma,p,\eta}(v)=\int_{1}^{\eta}g_{p\gamma-\beta,p,\theta^{p}}(v)m^{\sharp}_{\beta,\gamma,p,\eta}(\theta)\mathrm{d}\theta$

where

$\displaystyle m^{\sharp}_{\beta,\gamma,p,\eta}(\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{C^{*}_{\beta,\gamma,p,\eta}}\theta^{-1}\int_{1/\theta}^{1}(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{1}{pC^{*}_{\beta,\gamma,p,\eta}}\theta^{-1}\int_{1/\theta^{p}}^{1}(1-u)^{\gamma-1}u^{-1-\beta/p}\mathrm{d}u,\
\ 1<\theta<\eta.$

Thus, this is a GGSM with mixing density $m^{\sharp}_{\beta,\gamma,p,\eta}$.
The presence of the incomplete beta function in the mixing density gives the
distribution its name. It is easily checked that if $X\sim
m^{\sharp}_{\beta/p,\gamma,1,\eta^{p}}$ then $X^{1/p}\sim
m^{\sharp}_{\beta,\gamma,p,\eta}$.

We only focus on the case where $\gamma\in\\{1,2,3,\dots\\}$ as the other
values are not relevant for simulating OUTS processes. In this case, the
binomial theorem gives

$\displaystyle
m^{\sharp}_{\beta,\gamma,p,\eta}(\theta)=\frac{1}{C^{*}_{\beta,\gamma,p,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\theta^{-1}\frac{1-\theta^{\beta-pk}}{pk-\beta},\quad
1<\theta<\eta,$ (8)

where we replace $\frac{1-\theta^{\beta-pk}}{pk-\beta}$ by its limiting value
of $\ln\theta$ when $\beta=pk$. Integrating shows that the cdf is

$\displaystyle
M^{\sharp}_{\beta,\gamma,p,\eta}(y)=\frac{1}{C^{*}_{\beta,\gamma,p,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\frac{(pk-\beta)\ln y-1+y^{\beta-pk}}{(pk-\beta)^{2}},\quad
1<y<\eta,$ (9)

where we replace $\frac{(pk-\beta)\log y-1+y^{\beta-pk}}{(pk-\beta)^{2}}$ by
its limiting value of $(\ln y)^{2}/2$ when $\beta=pk$. Similarly, we can check
that

$C^{*}_{\beta,\gamma,p,\eta}=\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\frac{(pk-\beta)\log\eta-1+\eta^{\beta-pk}}{(pk-\beta)^{2}}.$

######  Proposition 2.3.

If $W\sim\mathrm{IBGM}(\beta,\gamma,p,\eta)$ with
$\gamma\in\\{1,2,3,\dots\\}$, then for any $\xi>\beta-\gamma p$ we have

$\mathrm{E}[W^{\xi}]=\frac{\Gamma\left(\gamma+(\xi-\beta)/p\right)}{\Gamma\left(\gamma-\beta/p\right)C^{*}_{\beta,\gamma,p,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}\frac{(-1)^{k}}{pk-\beta}\left(\frac{1-\eta^{-\xi}}{\xi}+\frac{1-\eta^{\beta-
pk-\xi}}{\beta-pk-\xi}\right),$

where we replace
$\frac{1}{pk-\beta}\left(\frac{1-\eta^{-\xi}}{\xi}+\frac{1-\eta^{\beta-
pk-\xi}}{\beta-pk-\xi}\right)$ by
$\frac{1-\eta^{-\xi}(\xi\ln\eta+1)}{\xi^{2}}$ if $\beta=pk$.

We can use Algorithms GGSM1 and GGSM2 to simulate from the IBGM distribution.
In this case Algorithm GGSM1 reduces to the following.

Algorithm IBGM1. Simulation from $\mathrm{IBGM}(\beta,\gamma,p,\eta)$.
Step 1. Independently simulate $Y\sim\mathrm{Ga}(\gamma-\beta/p,1)$ and $Z\sim
m^{\sharp}_{\beta/p,\gamma,1,\eta^{p}}$.
Step 2. Return $\left(Y/Z\right)^{1/p}$.

To specialize Algorithm GGSM2 let
$\varphi_{2}(u)=p^{-1}e^{u}\varphi_{2,n}(u)/\varphi_{2,d}(u)$, where

$\displaystyle\varphi_{2,n}(u)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}\frac{(-1)^{k}}{pk-\beta}u^{\beta/p-\gamma}\int_{u}^{u\eta^{p}}e^{-\theta}\left(\theta^{\gamma-\beta/p-1}-\theta^{\gamma-k-1}u^{k-\beta/p}\right)\mathrm{d}\theta$

and

$\displaystyle\varphi_{2,d}(u)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}\frac{(-1)^{k}}{pk-\beta}\left(\frac{\eta^{p\gamma-\beta}-1}{p\gamma-\beta}-\frac{\eta^{p(\gamma-k)}-1}{p(\gamma-k)}\right).$

In the above, if $k=\beta/p$, we replace the summand in $\varphi_{2,n}$ by

$\displaystyle p^{-1}{\gamma-1\choose
k}(-1)^{k}\int_{1}^{\eta^{p}}e^{-u\theta}\theta^{\gamma-\beta/p-1}\ln\theta\mathrm{d}\theta$

and the summand in $\varphi_{2,d}$ by

$\displaystyle{\gamma-1\choose
k}(-1)^{k}\frac{\eta^{p\gamma-\beta}((p\gamma-\beta)\ln\eta-1)+1}{(p\gamma-\beta)^{2}}.$

We note that for $k\neq\beta/p$ the integral in $\varphi_{2,n}$ can be written
in terms of incomplete gamma functions. With this notation we can specialize
Algorithm GGSM2 as follows.

Algorithm IBGM2. Simulation from $\mathrm{IBGM}(\beta,\gamma,p,\eta)$.
Step 1. Independently simulate $U\sim U(0,1)$ and
$Y\sim\mathrm{Ga}(\gamma-\beta/p,1)$.
Step 2. If $U\leq\varphi_{2}(Y)$ return $Y^{1/p}$, otherwise go back to step
1.

In this case, the probability of rejection is $1/V_{2}$, where

$V_{2}=\frac{\int_{1}^{\eta}\theta^{p\gamma-\beta-1}\int_{1/\theta}^{1}(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u\mathrm{d}\theta}{\int_{1}^{\eta}\theta^{-1}\int_{1/\theta}^{1}(1-u^{p})^{\gamma-1}u^{-1-\beta}\mathrm{d}u\mathrm{d}\theta}.$

Applying L’Hôpital’s rule and Leibniz Rule shows that $V_{2}\to 1$ as
$\eta\downarrow 1$.

We now turn to the problem of simulation from
$m^{\sharp}_{\beta,\gamma,p,\eta}$. We begin with the important case when
$\gamma=1$. In this case, $m^{\sharp}_{\beta,\gamma,p,\eta}$ does not depend
on the parameter $p$. When $\gamma=1$ and $\beta=0$ we have

$m^{\sharp}_{0,1,p,\eta}(\theta)=\frac{2\log\theta}{\theta(\log\eta)^{2}},\quad
1<\theta<\eta.$

It is easily checked that we can use the following algorithm in this case.

Algorithm M♯0. Simulation from $m^{\sharp}_{0,1,p,\eta}$.
Step 1. Simulate $U\sim U(0,1)$ and set $X=\eta^{\sqrt{U}}$.
Step 2. Return $X$.

When $\gamma=1$ and $\beta\neq 0$ we get

$m^{\sharp}_{\beta,1,p,\eta}(\theta)=\frac{\beta}{\eta^{\beta}-\beta\log\eta-1}\frac{\theta^{\beta}-1}{\theta},\quad
1<\theta<\eta.$

Simulation from this distribution was studied in [22] and [28]. We follow the
approach given in [28]. The idea is to first observe that $X\sim
m^{\sharp}_{\beta,1,p,\eta}$ can be represented as $X\stackrel{{\scriptstyle
d}}{{=}}\eta^{W}$, where the pdf of $W$ is

$f_{W}(w)=\frac{\beta\log\eta}{\eta^{\beta}-\beta\log\eta-1}\left(\eta^{\beta
w}-1\right),\quad 0\leq w\leq 1.$

Since $f_{W}$ is monotone and convex in $[0,1]$, simulation can be done in a
fast and efficient way by using the decomposition method illustrated in
Section 4.3 of [6]. Let $0=w_{0}<w_{1}<\cdots<w_{L-1}<w_{L}=1$ for some
positive integer $L$ and define a partition of $(0,1]$ given by the disjoint
intervals $\mathcal{I}_{\ell}=(w_{\ell-1},w_{\ell}],\;\ell=1,\dots,L$. From
the definition of a convex function, it follows that for each
$\ell=1,2,\dots,L$

$f_{W}(w)\leq\frac{f_{W}(w_{\ell})-f_{W}(w_{\ell-1})}{w_{\ell}-w_{\ell-1}}\,(w-w_{\ell-1})+f_{W}(w_{\ell-1})=:g_{\ell}(w),\
\ w\in\mathcal{I}_{\ell}.$

Now let

$\displaystyle\bar{g}_{\ell}(w)=\frac{g_{\ell}(w)}{q_{\ell}},\quad
q_{\ell}=\int_{\mathcal{I}_{\ell}}g_{\ell}(w)\mathrm{d}w,\qquad
p_{\ell}=\frac{q_{\ell}}{V_{L}},\qquad V_{L}=\sum_{\ell=1}^{L}q_{\ell}$

and note that each $\bar{q}_{\ell}$ is a pdf and that

$\displaystyle p(\ell)=p_{\ell},\ \ \ \ell=1,2,\dots,L$ (10)

is a pmf. With this notation we have

$\displaystyle f_{W}(w)\leq V_{L}\,\bar{g}_{L}(w),$ (11)

where

$\bar{g}_{L}(w)=\frac{1}{V_{L}}\sum_{\ell=1}^{L}g_{\ell}(w)\mathds{1}_{\mathcal{I}_{\ell}}(w)=\sum_{\ell=1}^{L}p_{\ell}\,\bar{g}_{\ell}(w)\mathds{1}_{\mathcal{I}_{\ell}}(w)$

is a pdf. Using (11) we can set up a rejection sampling method for simulation
from $f_{W}$. However, as observed in [28], the probability of acceptance can
be made arbitrarily close to one when the intervals are of equal length and
$L$ is large enough. Thus, in this case, we can skip the rejection sampling
step and just use the approximation $f_{W}\approx\bar{g}_{L}$. Simulation
results in [28], suggest that this approximation is very fast and works very
well. It is almost exact. This leads to the following approximate simulation
method for $m^{\sharp}_{\beta,1,p,\eta}$ with $\beta\neq 0$.

Algorithm M♯-CS. Approximate simulation from $m^{\sharp}_{\beta,1,p,\eta}$
with $\beta\neq 0$.
Step 1. Simulate $\ell\sim p$, where $p$ is the pmf in (10).
Step 2. Simulate $W_{\ell}\sim\bar{g}_{\ell}$.
Step 3. Return $X=\eta^{W_{\ell}}$.

In the above, simulation from $\bar{g}_{\ell}$ is straightforward as the pdf
is a linear function and we can use the inverse transform method, as its cdf
has a simple form. In practice we used the the _random.triangular_ routine in
the _numpy_ package for _Python_. We now turn to the case $\gamma>1$. A simple
brute force approach is to numerically invert the formula for
$M^{\sharp}_{\beta,\gamma,p,\eta}$ as given in (9). Denoting this inverse
function by $(M^{\sharp}_{\beta,\gamma,p,\eta})^{-1}$, leads to the following
algorithm.

Algorithm M♯1. Simulation from $m^{\sharp}_{\beta,\gamma,p,\eta}$ with
$\gamma\in\\{1,2,3,\dots\\}$.
Step 1. Simulate $U\sim U(0,1)$.
Step 2. Return $(M^{\sharp}_{\beta,\gamma,p,\eta})^{-1}(U)$.

We also develop a rejection sampling algorithm, which follows from the fact
that for $\gamma>1$

$m^{\sharp}_{\beta,\gamma,p,\eta}(\theta)\leq\frac{C^{*}_{\beta,1,p,\eta}}{C^{*}_{\beta,\gamma,p,\eta}}m^{\sharp}_{\beta,1,p,\eta}(\theta).$

Let

$\varphi^{\sharp}_{1}(u)=\frac{\beta}{u^{\beta}-1}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}(-1)^{k}\frac{1-u^{\beta-pk}}{pk-\beta}$

where we interpret $\frac{\beta}{u^{\beta}-1}$ as $1/\ln u$ when $\beta=0$ and
$\frac{1-u^{\beta-pk}}{pk-\beta}$ as $\ln u$ when $\beta=pk$.

Algorithm M♯2. Simulation from $m^{\sharp}_{\beta,\gamma,p,\eta}$ with
$\gamma\in\\{2,3,\dots\\}$.
Step 1. Independently simulate $U\sim U(0,1)$ and $Y\sim
m^{\sharp}_{\beta,1,p,\eta}$.
Step 2. If $U\leq\varphi^{\sharp}_{1}(Y)$ return $Y$, otherwise go back to
step 1.

In this case the probability of acceptance on any given iteration is given by
$\frac{C^{*}_{\beta,\gamma,p,\eta}}{C^{*}_{\beta,1,p,\eta}}$. Applying
L’Hôpital’s rule and Leibniz Rule shows that this approaches $0$ as
$\eta\downarrow 1$. Nevertheless, in simulations we found that this methods
works well for choices of $\eta$ that are not too close to $1$.

We end this section by noting that one can derive another rejection sampling
algorithm by taking only the positive terms of (8) as was done for the IGa
law. Unfortunately, in this case, simulation from the law of the normalized
positive sum is not straightforward as it requires an additional rejection
sampling step. For this reason we do not to consider this approach here.

###  2.3 Difference Generalized Gamma Distribution

Let $F$ be the cdf of some distribution with support contained in $[0,\infty)$
and consider the function

$f(x)=\frac{F(\eta
x)-F(x)}{x\log\eta}=\frac{1}{x\log\eta}\int_{(x,x\eta]}\mathrm{d}F(v),\ \ x>0$

for some $\eta>0$. Is is readily checked that this is a pdf. Such pdf’s arise
in the study of the transition laws of OU processes, where the BDLP is
compound Poisson, see [33]. We are interested in the case where $F$ is the cdf
of the $\mathrm{GGa}(\gamma,p,1)$ distribution. In this case, the pdf becomes

$\displaystyle
h_{\gamma,p,\eta}(x)=\frac{1}{x\log\eta}\int_{x}^{x\eta}g_{\gamma,p,1}(u)\mathrm{d}u=\frac{1}{\log\eta}\int_{1}^{\eta}g_{\gamma,p,1}(\theta
x)\mathrm{d}\theta.$

We call this the Difference Generalized Gamma Distribution and denote it by
$\mathrm{DGGa}(\gamma,p,\eta)$. It is readily checked that

$\displaystyle
h_{\gamma,p,\eta}(x)=\int_{1}^{\eta}g_{\gamma,p,\theta^{p}}(x)m_{0,1,1,\eta}(\theta)\mathrm{d}\theta,$

where

$m_{0,1,1,\eta}(\theta)=\frac{1}{\log\eta}\theta^{-1}\ \ 1<\theta<\eta.$

is a special case of the pdf given in (5). We can use Algorithm M0 to simulate
from $m_{0,1,1,\eta}$ and we can combine this with Algorithm GGSM1 to simulate
from $\mathrm{DGGa}(\gamma,p,\eta)$.

##  3 $p$-Tempered $\alpha$-Stable Distributions

Fix $\alpha\in(-\infty,2)$ and $p>0$. A $p$-tempered $\alpha$-stable
($\mathrm{TS}^{p}_{\alpha}$) distribution $\mu$ on $\mathbb{R}^{d}$ has a
characteristic function of the form $\hat{\mu}(z)=e^{C_{\mu}(z)}$, where

$c_{\mu}(z)=i\langle
b,z\rangle+\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\left(e^{it\langle
x,z\rangle}-1-it\langle x,z\rangle 1_{[\alpha\geq
1]}\right)t^{-1-\alpha}e^{-t^{p}}\mathrm{d}tR(\mathrm{d}x),$

$b\in\mathbb{R}^{d}$, and $R$ is a finite Borel measure on $\mathbb{R}^{d}$
satisfying $R(\\{0\\})=0$ and

$\displaystyle\int_{|x|>2}|x|^{\alpha}R(\mathrm{d}x)<\infty$
$\displaystyle\mbox{ if }\alpha\in(0,2)\setminus\\{1\\}$
$\displaystyle\int_{|x|>2}|x|\log|x|R(\mathrm{d}x)<\infty$
$\displaystyle\mbox{ if }\alpha=1$
$\displaystyle\int_{|x|>2}\log|x|R(\mathrm{d}x)<\infty$ $\displaystyle\mbox{
if }\alpha=0.$

No additional assumptions on $R$ are needed when $\alpha<0$. We denote this
distribution by $\mathrm{TS}^{p}_{\alpha}(R,b)$. We call $b$ the shift and $R$
the Rosiński measure after the author of [24]. One can consider extensions to
certain cases where $R$ is not a finite measure (see [9]), but we will not do
so here. The class of $\mathrm{TS}^{p}_{\alpha}$ distributions with $p=1$ and
$\alpha\in(0,2)$ was introduce in [24] and the class with $p=2$ and
$\alpha\in[0,2)$ was introduced in [20]. The general class was introduced in
[8].

Every $p$-tempered $\alpha$-stable distribution is infinitely divisible and
the Lévy measure of $\mathrm{TS}^{p}_{\alpha}(R,b)$ is given by

$\displaystyle
M(B)=\int_{\mathbb{R}^{d}}\int_{0}^{\infty}1_{B}(xt)t^{-1-\alpha}e^{-t^{p}}\mathrm{d}tR(\mathrm{d}x),\
\ \ B\in\mathfrak{B}(\mathbb{R}^{d})$ (12)

where $\mathfrak{B}(\mathbb{R}^{d})$ denotes the Borel sets on
$\mathbb{R}^{d}$. Formulas for the cumulants of $\mathrm{TS}^{p}_{\alpha}$
distributions are given in Theorem 2.16 of [9]. For simplicity we only recall
the formulas in the one-dimensional ($d=1$) case. In this case for
distribution $\mathrm{TS}^{p}_{\alpha}(R,b)$ if
$\int_{|x|>1}|x|^{k}R(\mathrm{d}x)<\infty$, then the $k$th cumulant exists and
is given for $\alpha\in(-\infty,1)$ by

$\displaystyle\mathrm{c}_{k}=p^{-1}\Gamma\left(\frac{k-\alpha}{p}\right)\int_{\mathbb{R}}x^{k}R(\mathrm{d}x)+1_{[k=1]}b.$
(13)

If $\alpha\in[1,2)$, then this formula still holds for $k\geq 2$, but for
$k=1$ it is given by $\mathrm{c}_{1}=b$. For more on $p$-tempered
$\alpha$-stable distributions and their associate Lévy processes see [9] and
the references therein.

While we present our results for general Rosiński measures $R$, we are
especially interest in the class of so-called $p$-rapidly decreasing tempered
stable ($p$-RDTS) distributions, see [12]. These correspond to the case where
the dimension $d=1$ and the Rosiński measure is of the form
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$ for some
$c,\beta>0$. In this case, after a change of variables, we get

$c_{\mu}(z)=ibz+c\int_{0}^{\infty}\left(e^{itz}-1-itz1_{[\alpha\geq
1]}\right)t^{-1-\alpha}e^{-(\beta t)^{p}}\mathrm{d}t.$

If one understands these distributions, then one can easily extend to the
bilateral case, where
$R(\mathrm{d}x)=c_{-}\beta_{-}^{\alpha}\delta_{-1/\beta_{-}}(\mathrm{d}x)+c_{+}\beta_{+}^{\alpha}\delta_{1/\beta_{+}}(\mathrm{d}x)$
for some $c_{-},c_{+},\beta_{-},\beta_{+}>0$. When $p=1$ these are sometimes
called classical tempered stable (CTS) distributions.

##  4 OU Processes

Let $L=\\{L_{t}:t\geq 0\\}$ be a Lévy process on $\mathbb{R}^{d}$. Fix
$\lambda>0$ and define a process $Y=\\{Y_{t}:t\geq 0\\}$ as the strong
solution of the stochastic differential equation (_SDE_)

$\mathrm{d}Y_{t}=-\lambda Y_{t}\mathrm{d}t+\mathrm{d}L_{t},\qquad Y_{0}=Y\ \
a.s.$

This process can be written as

$Y_{t}=Y_{0}\,e^{-\lambda t}+\int_{0}^{t}e^{-\lambda\,(t-s)}\mathrm{d}L_{s}.$

We say that $Y$ is an OU process with parameter $\lambda$ and that $L$ is the
BDLP. We refer to the distribution of $L_{1}$ as the BDLP distribution. Every
OU process is a Markov process and, so long as $\mathrm{E}[\log(|L_{1}|\vee
e)]<\infty$, the process has a stationary (also sometimes called a limiting)
distribution. An OU process whose stationary distribution is
$\mathrm{TS}^{p}_{\alpha}$ is called a TSOU process and an OU process whose
BDLP distribution is $\mathrm{TS}^{p}_{\alpha}$ is called an OUTS process.

In the remainder of this section we study the transition laws of both TSOU and
OUTS processes. We begin by giving formulas for their cumulants. For
simplicity we focus on the one-dimensional ($d=1$) case. In [27] a simple
formula relating the cumulants of the transition law of an OU process and
those of the stationary law are provided. Specifically, it is shown there
that, if the stationary law has a finite $k$th cumulant, then so does the
transition law and, in this case, the $k$th cumulant of the conditional
distribution of $Y_{s+t}$ given $Y_{s}=y$ is given by

$\hat{c}_{k,t}=ye^{-\lambda t}1_{[k=1]}+(1-e^{-k\lambda t})c_{k},$

where $c_{k}$ is the $k$th cumulant of the stationary law. Thus, when $d=1$
and the stationary law is $\mathrm{TS}^{p}_{\alpha}(R,b)$, if
$\int_{|x|>1}|x|^{k}R(\mathrm{d}x)<\infty$, then

$\displaystyle\hat{c}_{k,t}=\left(ye^{-\lambda t}+(1-e^{-k\lambda
t})b\right)1_{[k=1]}+(1-e^{-k\lambda
t})p^{-1}\Gamma\left(\frac{k-\alpha}{p}\right)\int_{\mathbb{R}}x^{k}R(\mathrm{d}x).$
(14)

The one exception to this formula is that when $k=1$ and $\alpha\in[1,2)$ we
have $\hat{c}_{1,t}=\left(ye^{-\lambda t}+(1-e^{-k\lambda t})b\right)$. We now
turn to the cumulants of the transition law of an OUTS process. In this case,
combining Proposition 3.12 in [4] with Lemma 17.1 in [29] shows that, so long
as the $k$th cumulant of the BDLP distribution exists, the $k$th cumulant of
the conditional distribution of $Y_{s+t}$ given $Y_{s}=y$ exists and is given
by

$\check{c}_{k,t}=ye^{-\lambda t}1_{[k=1]}+\frac{1-e^{-k\lambda
t}}{k\lambda}c_{k},$

where $c_{k}$ is the $k$th cumulant of the BDLP distribution. Thus, when $d=1$
and the BDLP distribution is $\mathrm{TS}^{p}_{\alpha}(R,b)$, if
$\int_{|x|>1}|x|^{k}R(\mathrm{d}x)<\infty$, then

$\displaystyle\check{c}_{k,t}=\left(ye^{-\lambda t}+\frac{1-e^{-k\lambda
t}}{k\lambda}b\right)1_{[k=1]}+\frac{1-e^{-k\lambda t}}{k\lambda
p}\Gamma\left(\frac{k-\alpha}{p}\right)\int_{\mathbb{R}}x^{k}R(\mathrm{d}x).$
(15)

We must again modify this formula when $k=1$ and $\alpha\in[1,2)$. In this
case we have $\check{c}_{1,t}=\left(ye^{-\lambda t}+(1-e^{-k\lambda
t})(k\lambda)^{-1}b\right)$.

We now give explicit representations for the transition laws of both TSOU and
OUTS processes and discuss simulation. These are given in $d$-dimensions. We
begin with TSOU processes. Since only selfdecomposable distributions can serve
is stationary distributions of OU processes, we only consider the case when
$\alpha\in[0,2)$ as $\mathrm{TS}^{p}_{\alpha}$ distributions are not
selfdecomposable when $\alpha<0$, see Proposition 3.14 in [9]. The following
result is given in [13].

######  Theorem 4.1.

Let $Y=\\{Y_{t}:t\geq 0\\}$ be a TSOU process with parameter $\lambda>0$ and
stationary distribution $\mathrm{TS}^{p}_{\alpha}(R,b)$ with $p>0$,
$\alpha\in[0,2)$, and $R\neq 0$. Set $\gamma=1+\lfloor\alpha/p\rfloor$. If
$t>0$, then, given $Y_{s}=y$, we have

$\displaystyle Y_{s+t}\stackrel{{\scriptstyle d}}{{=}}e^{-\lambda
t}y+(1-e^{-\lambda t})b-\sum_{n=0}^{\gamma-1}b_{n}+X_{0}+e^{-\lambda
t}\sum_{n=1}^{\gamma-1}X_{n}+\sum_{j=n}^{N}V_{n}W_{n},$ (16)

where $b_{0},\dots,b_{\gamma-1}\in\mathbb{R}^{d}$ are constants and
$N,X_{0},X_{1},\dots,X_{\gamma-1}$, $V_{1},V_{2},\dots$, $W_{1},W_{2},\dots$
are independent random variables with:
1\. $X_{0}\sim\mathrm{TS}^{p}_{\alpha}(R_{0},0)$ with
$R_{0}(\mathrm{d}x)=(1-e^{-\alpha\lambda t})R(\mathrm{d}x)$,
2\. if $\gamma\geq 2$ then $X_{n}\sim\mathrm{TS}^{p}_{\alpha-np}(R_{n},0)$
with $R_{n}(\mathrm{d}x)=\frac{1}{n!}(1-e^{-p\lambda t})^{n}R(\mathrm{d}x)$
for $n=1,2,\dots,(\gamma-1)$,
3\. $V_{1},V_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}R^{1}$,
where $R^{1}(\mathrm{d}x)=R(\mathrm{d}x)/R(\mathbb{R}^{d})$,
4\.
$W_{1},W_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{IGa}(\alpha,\gamma,p,e^{p\lambda
t})$,
5\. $N$ has a Poisson distribution with mean $e^{-\alpha\lambda
t}R(\mathbb{R}^{d})K_{\alpha,\gamma,p,e^{p\lambda t}}$,
6.

$b_{0}=\left\\{\begin{array}[]{ll}e^{-\alpha\lambda
t}\int_{\mathbb{R}^{d}}xR(\mathrm{d}x)K_{(\alpha-1),\gamma,p,e^{p\lambda
t}}&\alpha\in[1,2)\\\ 0&\alpha\in[0,1)\end{array}\right.,$

and if $\gamma\geq 2$ then for $n=1,2,\dots,(\gamma-1)$

$b_{n}=\left\\{\begin{array}[]{ll}e^{-\lambda
t}\int_{\mathbb{R}^{d}}xR_{n}(\mathrm{d}x)p^{-1}\Gamma\left(\frac{1-\alpha+np}{p}\right)&1\leq\alpha<1+np\\\
0&\mbox{otherwise}\end{array}\right..$

Note that when $\alpha=0$ we have $\gamma=1$, $b_{0}=0$, and $X_{0}=0$ with
probability one, thus the transition law is essentially just compound Poisson.
Note further, that the $\mathrm{IGa}$ distribution needed in the theorem has
parameter $\eta=e^{p\lambda t}$. When simulating a TSOU process on a finite
grid, one typically takes a small time step $t>0$. Thus one often uses a value
of $\eta$ that is close to $1$. Next, we turn to OUTS processes. In this case
we can allow for any $\alpha\in(-\infty,2)$. To the best of our knowledge the
transition law has not been studied previously in this case, except for CTS
and closely related distributions and only in the one-dimensional case with
$\alpha\in[0,1)$.

######  Theorem 4.2.

Let $Y=\\{Y_{t}:t\geq 0\\}$ be an OUTS process with parameter $\lambda>0$ and
BDLP distribution $\mathrm{TS}^{p}_{\alpha}(\lambda R,\lambda b)$ with $p>0$,
$\alpha\in(-\infty,2)$, and $R\neq 0$. If $\alpha\in(0,2)$ set
$\gamma=1+\lfloor\alpha/p\rfloor$, otherwise set $\gamma=1$. If $t>0$, then,
given $Y_{s}=y$, we have

$\displaystyle Y_{s+t}\stackrel{{\scriptstyle d}}{{=}}e^{-\lambda
t}y+(1-e^{-\lambda t})b-\sum_{n=0}^{\gamma-1}b^{*}_{n}+e^{-\lambda
t}\sum_{n=0}^{\gamma-1}X_{n}+\sum_{n=1}^{N}V_{n}W_{n},$ (17)

where $b_{0},\dots,b_{\gamma-1}\in\mathbb{R}^{d}$ are constants and
$N,X_{0},X_{1},\dots,X_{\gamma-1}$, $V_{1},V_{2},\dots$, $W_{1},W_{2},\dots$
are independent random variables with:
1\. $X_{0}\sim\mathrm{TS}^{p}_{\alpha}(R^{*}_{0},0)$ with
$R^{*}_{0}(\mathrm{d}x)=\frac{e^{\alpha\lambda t}-1}{\alpha}R(\mathrm{d}x)$,
2\. if $\gamma\geq 2$ then $X_{n}\sim\mathrm{TS}^{p}_{\alpha-np}(R^{*}_{n},0)$
with $R^{*}_{n}(\mathrm{d}x)=\kappa_{\lambda,t,n}R(\mathrm{d}x)$ and
$\kappa_{\lambda,t,n}=\int^{1}_{e^{-\lambda
t}}\frac{(1-u^{p})^{n}}{n!}u^{-1-\alpha}\mathrm{d}u$, for
$n=1,2,\dots,(\gamma-1)$,
3\. $V_{1},V_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}R^{1}$,
where $R^{1}(\mathrm{d}x)=R(\mathrm{d}x)/R(\mathbb{R}^{d})$,
4\.
$W_{1},W_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{IBGM}(\alpha,\gamma,p,e^{\lambda
t})$,
5\. $N$ has a Poisson distribution with mean
$\frac{pC_{\alpha,\gamma,p,e^{\lambda t}}R(\mathbb{R}^{d})}{(\gamma-1)!}$,
6.

$b^{*}_{0}=\left\\{\begin{array}[]{ll}\frac{pC_{\alpha,\gamma,p,e^{\lambda
t}}R(\mathbb{R}^{d})}{(\gamma-1)!}\mathrm{E}[V_{1}]\mathrm{E}[W_{1}]&\alpha\in[1,2)\\\
0&\alpha<1\end{array}\right.,$

and if $\gamma\geq 2$ then for $n=1,2,\dots,(\gamma-1)$

$b^{*}_{n}=\left\\{\begin{array}[]{ll}e^{-\lambda
t}\int_{\mathbb{R}^{d}}xR^{*}_{n}(\mathrm{d}x)p^{-1}\Gamma\left(\frac{1-\alpha+np}{p}\right)&1\leq\alpha<1+np\\\
0&\mbox{otherwise}\end{array}\right..$

In the theorem and its proof, when $\alpha=0$, we interpret
$\frac{e^{\alpha\lambda t}-1}{\alpha}$ by its limiting value of $t\lambda$. We
can, of course, state the theorem for the case where the BDLP distribution is
$\mathrm{TS}^{p}_{\alpha}(R,b)$ instead of $\mathrm{TS}^{p}_{\alpha}(\lambda
R,\lambda b)$. However, the formulas would be a bit more complicated and we do
not do so here.

Note that Theorem 4.2 holds even if the OUTS process does not have a
stationary distribution. A stationary distribution exists if and only if

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

where $M$ is the Lévy measure of $\mathrm{TS}^{p}_{\alpha}(R,b)$. A simple
sufficient condition is

$\int_{|x|>2}|x|^{\epsilon}R(\mathrm{d}x)<\infty\mbox{ for some }\epsilon>0.$

Under our assumptions, this always holds for $\alpha\in(0,2)$, see [9].

While Theorem 4.2 holds for any $\alpha\in(-\infty,2)$, when $\alpha<0$ we can
get a significantly simpler representation as, in this case,
$\mathrm{TS}^{p}_{\alpha}$ distributions are simply compound Poisson (with
drift). In the one-dimensional case, a general representation of the
transition law of an OU process with a compound Poisson BDLP is given in [33].
Although we cannot use those results directly as we are in $d$-dimensions, our
results are related to the ones in that paper.

######  Theorem 4.3.

Let $Y=\\{Y_{t}:t\geq 0\\}$ be an OUTS process with parameter $\lambda>0$ and
BDLP distribution $\mathrm{TS}^{p}_{\alpha}(\lambda R,\lambda b)$ with $p>0$,
$\alpha\in(-\infty,0)$, and $R\neq 0$. If $t>0$, then, given $Y_{s}=y$, we
have

$\displaystyle Y_{s+t}\stackrel{{\scriptstyle d}}{{=}}e^{-\lambda
t}y+(1-e^{-\lambda t})b+\sum_{n=1}^{N}V_{n}W_{n},$ (18)

where $N$, $V_{1},V_{2},\dots$, $W_{1},W_{2},\dots$ are independent random
variables with:
1\. $V_{1},V_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}R^{1}$,
where $R^{1}(\mathrm{d}x)=R(\mathrm{d}x)/R(\mathbb{R}^{d})$,
2\.
$W_{1},W_{2},\dots\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{DGGa}(|\alpha|,p,e^{\lambda
t})$,
3\. $N$ has a Poisson distribution with mean $p^{-1}\lambda
t\Gamma(|\alpha|/p)R(\mathbb{R}^{d})$.

Our main goal in studying the transition laws is to use them to simulate the
corresponding TSOU or OUTS process on a finite grid. To do this, we need a way
to simulate from the transition law, or equivalently from the various
components of this law. We have already discussed the simulation of
$\mathrm{IGa}$, $\mathrm{IBGM}$, and $\mathrm{DGGa}$ distributions in Section
2. There is no one approach for simulating from $R^{1}$ as it can be,
essentially, any probability measure on $\mathbb{R}^{d}$. However, when
simulating specifically TSOU processes, there is a way to avoid simulating
from $R^{1}$. In this case one can directly simulate the product $V_{i}W_{i}$,
where $V_{i}\sim R^{1}$ and $W_{i}$ has an $\mathrm{IGa}$ distribution, see
[13].

The remaining components of the transition law are $\mathrm{TS}^{p}_{\alpha-
np}$ for $n=0,1,\dots,\gamma-1$. There are several approaches for simulating
from these distributions. First, one can use the inverse transform method,
which requires one to numerically invert the cdf. While this method can work
well, the fact that there is no closed formula for the cdfs of
$\mathrm{TS}^{p}_{\alpha}$ distributions makes this method impractical in many
cases. Second, under mild assumptions, one can use the rejection sampling
approach of [10]. However, this method requires one to numerically calculate
pdfs, which may also be computationally intensive. A third approach is to use
an approximate method based on truncating an infinite series representation. A
number of such representations appear in the literature, see [24], [25], [17],
or [20]. We note that several of the methods discussed here are easier to
implement in the univariate case. An approach for extending univariate
simulation methods of $\mathrm{TS}^{p}_{\alpha}$ random variables to the
multivariate case is given in [31]. Finally, we note that numerical methods
for simulation and the evaluation of pdfs and cdfs of certain classes of
symmetric $\mathrm{TS}^{p}_{\alpha}$ distribution can be found in the SymTS
package [14] for the statistical software R.

######  Remark 2.

We are particularly interested in the class of $p$-RDTS distributions, which
correspond to the case where the dimension $d=1$ and
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$ for some
$c,\beta>0$. In this case simulation of the various components of the
transition law is fairly simple. First, we have
$R^{1}(\mathrm{d}x)=\delta_{1/\beta}(\mathrm{d}x)$ and thus if $V\sim R^{1}$
then $V=1/\beta$ with probability $1$. Second, a simple method for simulating
from $\mathrm{TS}^{p}_{\alpha}(R,b)$ is given in [12] for the case $\alpha<1$
and $p>1$. Finally, when $p=1$, this class reduces to the class of CTS
distributions. Exact simulation methods for CTS distributions are well known
and can be found in, e.g., [7], [16], [18], and the references therein.

##  5 Numerical Experiments

In this section we illustrate and compare the performance and effectiveness of
the simulation algorithms discussed in this paper. All simulations were
conducted using _Python_ with a $64$-bit Intel Core i7-7500U CPU @270-290 GHz,
8GB. We first investigate the performance of the simulation methods for the
IGa distribution as described in Section 2.1 and then the simulation methods
for the IBGM distribution presented in Section 2.2. Finally, we focus on the
generation of TSOU and OUTS processes on a finite grid. To ensure that we are
simulating from the correct distributions, we compare the empirical moments to
the true moments. For simplicity, for the OU processes we compare the
cumulants instead of the moments. To see how close the empirical values are to
the true values, we consider the relative errors as given by

$\text{err}\;\%=\frac{\text{true value - estimated value}}{\text{true
value}}\times 100\%.$

###  5.1 Results for IGa

In this section we compare the performance of four methods for simulating from
an IGa distribution, which are discussed in Section 2.1. Three of them are new
and use Algorithm IGa1 in conjunction with an algorithm for simulating from
$m_{\beta,\gamma,p,\eta}$. We denote these by _ARGS_ , _Inverse_ , and _ARBD_
and for $\gamma\geq 2$ they use Algorithms M1, M2, and M3, respectively. When
$\gamma=1$, Algorithms M1 and M2 are no longer meaningful and Algorithm M2
reduces to Algorithm M0. For this reason, when $\gamma=1$ we use Algorithm M0
for all three methods. The fourth method uses Algorithm IGa2 and is denoted
_ARG_. It was introduced in [13].

Simulation using the _Inverse_ method when $\gamma\geq 2$ depends on the
numerical inversion of the cdf given in (7). This, in turn, depends on an
initial guess which can, of course, affect the final computation time. Instead
of blindly taking the middle term $(\eta^{1/p}+1)/2$, we choose the initial
guess equal to the random variate drawn from the corresponding distribution
with $\gamma=1$. In other words, we start with the value returned by Algorithm
M0. Calculating this value is fast and its impact on the overall computation
time is negligible.

As discussed in Section 2.1, without loss of generality we take $p=1$. For the
other parameters, we take $\beta=0.9$, $\eta\in\\{1.1,2\\}$, and
$\gamma\in\\{1,2,3,4,5,10\\}$. The choice of the $\eta$’s stems from the fact
that when simulating a TSOU process on a finite grid one often needs $\eta>1$
close to $1$. Table 1 shows the computation times for the four methods. We use
the method _ARGS_ as the baseline and for it all values are given in seconds,
while the values for the other methods are given as multiplicative factors
with respect to it. We can see that the new _ARGS_ method performed the
fastest, while the new _ARBD_ method performed the slowest. Further, the new
methods _ARGS_ and _Inverse_ performed significantly faster than the _ARG_
method of [13].

To ensure that the methods are simulating from the correct distributions,
Table 2 shows the comparison between the first four true moments computed in
(3) and the empirical moments estimated based on $R=5\times 10^{4}$ simulated
values. In the interest of space, we only present the results for $p=1$,
$\beta=0.9$, and $\eta=2$. In all cases the _err %_ is small suggesting that
all methods are simulating from the correct distributions.

|  | $\gamma=1$ | $\gamma=2$ | $\gamma=3$ | $\gamma=4$ | $\gamma=5$ | $\gamma=10$
---|---|---|---|---|---|---|---
Method | $R$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$
ARGS (sec) | $1000$ | $0.00020$ | $0.00040$ | $0.0030$ | $0.0080$ | $0.0049$ | $0.0100$ | $0.0080$ | $0.0170$ | $0.0120$ | $0.0269$ | $0.0289$ | $0.6902$
$10000$ | $0.00125$ | $0.00289$ | $0.0498$ | $0.0549$ | $0.0678$ | $0.0858$ | $0.0817$ | $0.1685$ | $0.1067$ | $0.2822$ | $0.2604$ | $6.8407$
$20000$ | $0.00247$ | $0.00598$ | $0.1054$ | $0.0858$ | $0.1197$ | $0.1685$ | $0.1576$ | $0.3311$ | $0.2055$ | $0.5586$ | $0.5037$ | $13.559$
$50000$ | $0.00669$ | $0.00740$ | $0.2042$ | $0.2105$ | $0.2942$ | $0.4488$ | $0.4019$ | $0.7660$ | $0.4957$ | $1.4192$ | $1.2825$ | $38.205$
Inverse | $1000$ | $1$ | $1$ | $2.3$ | $1.6$ | $1.9$ | $2.5$ | $1.9$ | $2.2$ | $4.6$ | $2.0$ | $1.7$ | $1.9$
$10000$ | $1$ | $1$ | $2.9$ | $1.4$ | $1.9$ | $2.5$ | $3.1$ | $3.1$ | $4.2$ | $4.0$ | $3.0$ | $1.6$
$20000$ | $1$ | $1$ | $2.4$ | $2.8$ | $1.7$ | $5.2$ | $3.9$ | $3.1$ | $3.5$ | $4.2$ | $2.5$ | $1.9$
$50000$ | $1$ | $1$ | $4.0$ | $2.2$ | $3.1$ | $4.9$ | $3.0$ | $4.2$ | $4.7$ | $4.3$ | $2.7$ | $3.5$
ARG | $1000$ | $4.5$ | $6.4$ | $17$ | $6.8$ | $12$ | $9.2$ | $7.5$ | $7.6$ | $5.9$ | $6.4$ | $3.0$ | $4.7$
$10000$ | $5.7$ | $6.0$ | $6.4$ | $10$ | $8.4$ | $10$ | $7.5$ | $7.6$ | $6.6$ | $8.3$ | $3.4$ | $4.8$
$20000$ | $6.1$ | $6.5$ | $4.6$ | $12$ | $9.4$ | $11$ | $7.8$ | $8.8$ | $6.5$ | $11$ | $3.5$ | $3.9$
$50000$ | $12$ | $11$ | $8.5$ | $12$ | $9.6$ | $15$ | $7.7$ | $22$ | $6.8$ | $20$ | $3.4$ | $3.5$
ARBD | $1000$ | $1$ | $1$ | $21$ | $18$ | $24$ | $26$ | $19$ | $20$ | $15$ | $15$ | $6.8$ | $7.2$
$10000$ | $1$ | $1$ | $11$ | $25$ | $18$ | $21$ | $18$ | $21$ | $20$ | $17$ | $8.8$ | $8.5$
$20000$ | $1$ | $1$ | $12$ | $28$ | $21$ | $23$ | $20$ | $22$ | $33$ | $28$ | $11$ | $11$
$50000$ | $1$ | $1$ | $13$ | $30$ | $22$ | $22$ | $44$ | $54$ | $48$ | $47$ | $12$ | $10$

Table 1: Results for IGa: computation times. Here we take
$(\eta_{1},\eta_{2})=(1.1,2)$ and $(p,\beta)=(1,0.9)$. For _ARGS_ the values
are in seconds, otherwise they are the multiplicative factors with respect to
the _ARGS_ method

| $m_{1}$ | $m_{2}$
---|---|---
$\gamma$ | True | ARGS | Inverse | ARG | ARBD | True | ARGS | Inverse | ARG | ARBD
$1$ | $0.070$ | $1.0\%$ | $1.0\%$ | $-0.5\%$ | $1.0\%$ | $0.055$ | $1.3\%$ | $1.3\%$ | $1.9\%$ | $1.3\%$
$2$ | $0.694$ | $0.5\%$ | $0.1\%$ | $-0.3\%$ | $-2.8\%$ | $0.946$ | $0.8\%$ | $-0.5\%$ | $-0.5\%$ | $-5.1\%$
$3$ | $1.261$ | $0.3\%$ | $0.0\%$ | $0.2\%$ | $0.0\%$ | $2.395$ | $0.7\%$ | $-0.3\%$ | $0.1\%$ | $-0.2\%$
$4$ | $1.804$ | $0.1\%$ | $0.2\%$ | $0.0\%$ | $0.8\%$ | $4.372$ | $0.4\%$ | $0.4\%$ | $0.3\%$ | $1.6\%$
$5$ | $2.336$ | $0.0\%$ | $-0.2\%$ | $0.2\%$ | $1.1\%$ | $6.866$ | $0.0\%$ | $-0.5\%$ | $0.5\%$ | $2.3\%$
$10$ | $4.916$ | $0.2\%$ | $0.6\%$ | $-0.5\%$ | $-0.5\%$ | $26.96$ | $0.3\%$ | $1.5\%$ | $-0.6\%$ | $-3.7\%$
| $m_{3}$ | $m_{4}$
$\gamma$ | True | ARGS | Inverse | ARG | ARBD | True | ARGS | Inverse | ARG | ARBD
$1$ | $0.088$ | $2.9\%$ | $2.9\%$ | $4.0\%$ | $2.9\%$ | $0.212$ | $-3.5\%$ | $-3.5\%$ | $-2.9\%$ | $-3.5\%$
$2$ | $1.962$ | $0.8\%$ | $-1.5\%$ | $0.1\%$ | $-7.4\%$ | $5.546$ | $-0.3\%$ | $-2.5\%$ | $1.6\%$ | $3.4\%$
$3$ | $6.152$ | $1.0\%$ | $-0.9\%$ | $-0.3\%$ | $0.1\%$ | $20.13$ | $1.1\%$ | $-1.5\%$ | $-1.2\%$ | $1.3\%$
$4$ | $13.40$ | $1.0\%$ | $0.5\%$ | $0.9\%$ | $2.5\%$ | $50.06$ | $1.9\%$ | $0.8\%$ | $1.9\%$ | $3.0\%$
$5$ | $24.47$ | $-0.1\%$ | $-0.8\%$ | $0.8\%$ | $3.6\%$ | $103.0$ | $-0.3\%$ | $-1.0\%$ | $1.1\%$ | $5.3\%$
$10$ | $163.3$ | $0.4\%$ | $1.7\%$ | $-1.8\%$ | $3.4\%$ | $1085$ | $0.4\%$ | $5.1\%$ | $-4.0\%$ | $6.1\%$

Table 2: Results for IGa: moment comparison. Here we take
$(p,\beta,\eta)=(1,0.9,2)$ and evaluate empirical moments based on $R=5\times
10^{4}$ simulated values. Column _True_ gives the true values of the moments,
while the other columns give the err %.

###  5.2 Results for IBGM

In this section we compare the performance of three methods for simulating
from an IBGM distribution, which are presented in Section 2.2. The first
method, denoted _Inverse_ combines Algorithm IBGM1 with Algorithm M♯1\. In
this case, we always chose our initial guess for the numerical inversion to be
the midpoint, $(1+\eta)/2$. The second, denoted _ARGS_ combines Algorithm
IBGM1 with Algorithm M♯-CS when $\gamma=1$ and with Algorithm M♯2 when
$\gamma\geq 2$. The third, denoted _GGSM_ , uses Algorithm IBGM2. When
implementing Algorithm M♯2, we use Algorithm M♯-CS in the first step to
generate an observation from $m^{\sharp}_{\beta,1,p,\eta}$. In all cases, when
we use Algorithm M♯-CS we take $L=2000$ equally spaced intervals. This
algorithm was introduced in [28] and, while it is an approximate algorithm, it
works very well and is almost exact.

Table 3 presents the computation times of the different methods for several
choices of the parameters. Here, we take _Inverse_ as the baseline. For it all
values are given in seconds, while the values for the other methods are given
as multiplicative factors with respect to it. When $\gamma=1$, _ARGS_ , which
uses Algorithm M♯-CS, performed the fastest. However, when $\gamma\geq 2$
there was a dichotomy. In this case _GGSM_ is always faster than _ARGS_ for
$\eta=1.1$ and slower than _ARGS_ for $\eta=2$. This is likely related to the
asymptotic results (as $\eta\downarrow 1$) for the probability of acceptance
in Algorithms IBGM2 and $M^{\sharp}2$, which are given in Section 2.2. Method
_Inverse_ tends to work better for larger values for $\gamma$. To summarize,
when $\gamma=1$ it is better to use _ARGS_ , when $\gamma$ is large it is
better to use _Inverse_ , and when $\gamma\geq 2$ is not too big, the
situation depends on the value of $\eta$. In the context of the simulation of
an OUTS process on a finite grid, a larger value of $\eta$ corresponds to a
grid of time points with larger time-steps, whereas $\eta$ approaches $1$ as
the grid gets finer. Thus, in this case, the selection of the fastest approach
depending on the granularity of the grid.

Table 4 gives the comparison between the true and the empirical moments based
on $5\times 10^{4}$ simulated values. We see that all of the methods seem to
be simulating from the correct distributions. This is especially important to
note in the case of _ARGS_ as this method is only approximate. We see that it
works well and that the err % is no worse than it is for the other methods.

|  | $\gamma=1$ | $\gamma=2$ | $\gamma=3$ | $\gamma=4$ | $\gamma=5$ | $\gamma=10$
---|---|---|---|---|---|---|---
Method | $R$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$ | $\eta_{1}$ | $\eta_{2}$
Inverse (sec) | $1000$ | $0.4853$ | $1.0034$ | $0.5923$ | $0.5855$ | $1.2915$ | $0.6171$ | $0.8489$ | $1.2302$ | $0.0385$ | $0.4398$ | $0.0204$ | $0.7247$
$10000$ | $2.2132$ | $2.6426$ | $3.5021$ | $1.5841$ | $8.9314$ | $2.6201$ | $2.0610$ | $3.4214$ | $0.2153$ | $3.6592$ | $0.1602$ | $4.5543$
$20000$ | $3.9827$ | $3.2911$ | $8.8375$ | $5.7192$ | $15.896$ | $6.4364$ | $3.6689$ | $8.9187$ | $0.4276$ | $12.247$ | $0.3198$ | $12.514$
$50000$ | $12.636$ | $9.1716$ | $27.139$ | $17.987$ | $45.205$ | $23.394$ | $8.253$ | $36.293$ | $1.0157$ | $48.394$ | $0.792$ | $21.124$
ARGS | $1000$ | $0.0082$ | $0.0030$ | $1.6$ | $1.5$ | $0.5$ | $1.2$ | $0.8$ | $0.7$ | $18$ | $2.1$ | $386$ | $1.8$
$10000$ | $0.0050$ | $0.0034$ | $1.8$ | $4.6$ | $0.7$ | $3.6$ | $3.2$ | $2.5$ | $33$ | $2.5$ | $478$ | $2.9$
$20000$ | $0.0048$ | $0.0061$ | $1.4$ | $3.2$ | $0.8$ | $2.9$ | $3.7$ | $1.9$ | $33$ | $1.5$ | $455$ | $2.1$
$50000$ | $0.0050$ | $0.0030$ | $1.1$ | $2.5$ | $0.7$ | $2.1$ | $4.1$ | $1.4$ | $35$ | $1.0$ | $449$ | $3.1$
GGSM | $1000$ | $0.1665$ | $0.0885$ | $0.2$ | $3.0$ | $0.1$ | $8.4$ | $0.2$ | $1.5$ | $6.8$ | $6.1$ | $37$ | $219$
$10000$ | $0.3880$ | $0.2597$ | $0.3$ | $18$ | $0.2$ | $17$ | $1.0$ | $15$ | $12$ | $7.0$ | $46$ | $331$
$20000$ | $0.3176$ | $0.3285$ | $0.2$ | $9.8$ | $0.2$ | $14$ | $1.1$ | $7.1$ | $12$ | $3.9$ | $46$ | $241$
$50000$ | $0.1822$ | $0.2909$ | $0.2$ | $7.2$ | $0.2$ | $11$ | $1.2$ | $16$ | $13$ | $13$ | $48$ | $351$

Table 3: Results for IBGM: computation times. Here we take
$(\eta_{1},\eta_{2})=(1.1,2)$ and $(p,\beta)=(1,0.9)$. For _Inverse_ the
values are in seconds, otherwise they are the multiplicative factors with
respect to the _Inverse_ method

| $m_{1}$ | $m_{2}$
---|---|---
$\gamma$ | True | Inverse | ARGS | GGSM | True | Inverse | ARG | ARBD
$1$ | $0.0630$ | $0.2\%$ | $0.1\%$ | $0.3\%$ | $0.0450$ | $-0.2\%$ | $0.2\%$ | $1.6\%$
$2$ | $0.6566$ | $-0.1\%$ | $0.1\%$ | $0.4\%$ | $0.8393$ | $0.6\%$ | $0.4\%$ | $0.5\%$
$3$ | $1.2146$ | $-0.3\%$ | $-1.0\%$ | $-0.1\%$ | $2.2096$ | $0.5\%$ | $-1.8\%$ | $0.1\%$
$4$ | $1.7556$ | $0.1\%$ | $-2.3\%$ | $-0.3\%$ | $4.1220$ | $0.2\%$ | $-4.5\%$ | $-0.9\%$
$5$ | $2.2865$ | $0.1\%$ | $-3.0\%$ | $0.0\%$ | $6.5618$ | $-0.1\%$ | $-1.7\%$ | $0.0\%$
$10$ | $4.8739$ | $0.1\%$ | $-0.7\%$ | $-3.0\%$ | $26.470$ | $0.2\%$ | $-0.9\%$ | $-3.2\%$
| $m_{3}$ | $m_{4}$
$\gamma$ | True | Inverse | ARGS | GGSM | True | Inverse | ARGS | GGSM
$1$ | $0.0630$ | $-2.5\%$ | $-2.5\%$ | $4.6\%$ | $0.1346$ | $4.7\%$ | $6.4\%$ | $-4.9\%$
$2$ | $1.6182$ | $-2.1\%$ | $1.3\%$ | $0.1\%$ | $4.2236$ | $-5.3\%$ | $2.8\%$ | $0.1\%$
$3$ | $5.4023$ | $0.5\%$ | $-2.0\%$ | $0.4\%$ | $16.724$ | $3.5\%$ | $-0.7\%$ | $0.8\%$
$4$ | $12.190$ | $-0.1\%$ | $-5.8\%$ | $-1.7\%$ | $43.721$ | $-1.0\%$ | $-5.1\%$ | $-3.2\%$
$5$ | $22.749$ | $-0.3\%$ | $-4.6\%$ | $-0.1\%$ | $92.820$ | $-0.6\%$ | $-2.7\%$ | $-0.3\%$
$10$ | $158.67$ | $0.2\%$ | $-0.7\%$ | $-3.0\%$ | $1041.7$ | $0.3\%$ | $1.6\%$ | $-0.7\%$

Table 4: Results for IBGM: moment comparison. Here we take
$(p,\beta,\eta)=(1,0.9,2)$ and evaluate empirical moments based on $R=5\times
10^{4}$ simulated values. Column _True_ gives the true values of the moments,
while the other columns give the err %.

###  5.3 Results for TSOU processes

Theorem 4.1 characterizes the transition laws of TSOU processes with $p>0$ and
$\alpha\in[0,2)$. This can be used to simulate such a process on a finite grid
of times. In this section we illustrate this approach by performing a series
of simulations. We focus on the important class of $p$-RDTS distributions,
which correspond to the case where the dimension $d=1$ and the stationary
distribution is $\mathrm{TS}^{p}_{\alpha}(R,b)$, where
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$ for some
$c,\beta>0$. This means that
$R^{1}(\mathrm{d}x)=\delta_{1/\beta}(\mathrm{d}x)$ and hence that each
$V_{i}=1/\beta$ with probability $1$. For simplicity we take $b=0$ and for
tractability, we take $\alpha<1$ and $p>1$. In this case $\gamma=1$, we can
simulate $X_{0}$ using the method given in [12], and we can simulate from the
required $\mathrm{IGa}$ distribution using the _Inverse_ method, which
combines Algorithm IGa1 with Algorithm M0.

Figure 1 displays the sample trajectories of TSOU processes with several
choices of the parameters. These were simulated using the time-step $t=1/365$
over an equally-spaced grid with $365$ points. It is well-known that $p$-RDTS
distributions with $\alpha=0.5$ and $p=1$ reduce to the well-known class of
inverse Gaussian distributions. In our simulations we take $\alpha=0.5$ and
$p>1$. Thus, these processes are generalizations of inverse Gaussian OU
processes. We note that the transition laws of inverse Gaussian OU processes
were studied in [34].

Next, we check the correctness and of our algorithm. We simulate $10^{5}$
observations from the stationary law with a time step of $t=0.1$. A comparison
of the true cumulants and the empirical cumulants is given in Table 5 for
several choices of the parameters. The values of the true cumulants are
evaluated using (14) with the appropriate choice of $R$ and the parameters. We
see that for all cumulants the err % is small, which suggests that the
algorithm stemming from Theorem 4.1 is simulating from the correct
distribution.

Figure 1: Sample trajectories of TSOU processes with parameter $\lambda=10$,
initial value $y_{0}=0$, and time step $t=1/365$. The stationary distribution
is $\mathrm{TS}^{p}_{\alpha}(R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=1$,
$\beta=1$, $\alpha=0.5$, and $p\in\\{1.5,2,2.5,3\\}$.

|  | $c_{X,1}(0,t)$ | $c_{X,2}(0,t)$ | $c_{X,3}(0,t)$ | $c_{X,4}(0,t)$
---|---|---|---|---|---
$p$ | $\alpha$ | true | err % | true | err % | true | err % | true | err %
1.5 | $0.1$ | $0.628$ | $-0.61\%$ | $0.521$ | $-0.53\%$ | $0.617$ | $-0.34\%$ | $0.936$ | $1.25\%$
$0.3$ | $0.800$ | $-0.03\%$ | $0.541$ | $-0.04\%$ | $0.590$ | $-1.07\%$ | $0.850$ | $-2.93\%$
$0.5$ | $1.129$ | $0.06\%$ | $0.576$ | $-0.47\%$ | $0.572$ | $-2.52\%$ | $0.779$ | $-4.01\%$
$0.7$ | $1.935$ | $-0.04\%$ | $0.632$ | $0.29\%$ | $0.562$ | $0.10\%$ | $0.721$ | $-1.58\%$
$0.9$ | $6.104$ | $-0.71\%$ | $0.720$ | $-0.48\%$ | $0.562$ | $-1.74\%$ | $0.674$ | $-1.14\%$
2 | $0.1$ | $0.622$ | $-0.69\%$ | $0.446$ | $-0.76\%$ | $0.421$ | $-0.56\%$ | $0.481$ | $-0.01\%$
$0.3$ | $0.805$ | $0.53\%$ | $0.481$ | $0.47\%$ | $0.423$ | $1.06\%$ | $0.464$ | $4.53\%$
$0.5$ | $1.146$ | $0.31\%$ | $0.530$ | $0.28\%$ | $0.431$ | $-1.11\%$ | $0.451$ | $-4.04\%$
$0.7$ | $1.966$ | $-0.04\%$ | $0.599$ | $0.15\%$ | $0.443$ | $-0.75\%$ | $0.442$ | $-4.39\%$
$0.9$ | $6.154$ | $0.63\%$ | $0.699$ | $-0.78\%$ | $0.463$ | $-3.14\%$ | $0.436$ | $-4.21\%$
2.5 | $0.1$ | $0.625$ | $0.69\%$ | $0.419$ | $0.36\%$ | $0.353$ | $-0.65\%$ | $0.349$ | $-1.95\%$
$0.3$ | $0.813$ | $0.35\%$ | $0.460$ | $0.56\%$ | $0.365$ | $0.32\%$ | $0.348$ | $-4.66\%$
$0.5$ | $1.161$ | $-0.03\%$ | $0.515$ | $-0.98\%$ | $0.380$ | $-2.13\%$ | $0.348$ | $-3.77\%$
$0.7$ | $1.988$ | $0.04\%$ | $0.590$ | $0.03\%$ | $0.400$ | $0.50\%$ | $0.351$ | $2.29\%$
$0.9$ | $6.185$ | $-0.34\%$ | $0.696$ | $0.77\%$ | $0.427$ | $-2.16\%$ | $0.357$ | $3.67\%$
3 | $0.1$ | $0.630$ | $-0.36\%$ | $0.409$ | $0.24\%$ | $0.323$ | $0.97\%$ | $0.294$ | $-1.64\%$
$0.3$ | $0.822$ | $0.25\%$ | $0.453$ | $0.03\%$ | $0.339$ | $-0.54\%$ | $0.298$ | $-3.60\%$
$0.5$ | $1.173$ | $-0.22\%$ | $0.511$ | $-0.49\%$ | $0.358$ | $-1.28\%$ | $0.304$ | $-2.44\%$
$0.7$ | $2.005$ | $-0.18\%$ | $0.589$ | $-0.85\%$ | $0.381$ | $-3.16\%$ | $0.311$ | $-3.13\%$
$0.9$ | $6.206$ | $-0.62\%$ | $0.699$ | $-2.17\%$ | $0.411$ | $-2.23\%$ | $0.321$ | $-4.23\%$

Table 5: Comparison of the first four true cumulants with their estimated
values obtained from $10^{5}$ simulations from the transition law of a TSOU
process with $\lambda=10$, initial value $y_{0}=0$, and time step $t=0.1$. The
stationary distribution is $\mathrm{TS}^{p}_{\alpha}(R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=1$,
$\beta=1$, and with several choices for $\alpha$ and $p$.

###  5.4 OUTS processes

We now turn to the simulation of OUTS process on a finite grid. We again focus
on the case of $p$-RDTS distributions and for simplicity we assume that the
shift $b=0$. Here, we are assuming that the dimension $d=1$ and that the BDLP
distribution is $\mathrm{TS}^{p}_{\alpha}(\lambda R,0)$, where $\lambda>0$ is
the parameter of the OUTS process and
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$ for some
$c,\beta>0$.

We begin with the case $\alpha\geq 0$, for which the transition law is
characterized in Theorem 4.2. For tractability, we again focus on the case
$\alpha\in[0,1)$ and $p>1$. Here $\gamma=1$ and we can simulate the term
$X_{0}$ using the approach given in [12]. To simulate from the $\mathrm{IBGM}$
distribution we use the _ARGS_ method discussed in Section 5.2. We again take
$L=2000$ equally spaced intervals, which leads to an efficient approximate
simulation method. Figure 3 shows sample trajectories of OUTS processes for
several choices of the parameters. In all cases we take $y_{0}=0$ as the
initial value, a time step of $t=1/365$, $\alpha=0.5$, and $p>1$. With this
choice for $\alpha$, we can think of the processes as extensions of OU process
with inverse Gaussian BDLP distributions. Next, to check the correctness of
the algorithm, we simulate $10^{5}$ observations from the transition law with
a time step of $t=0.1$ and several values for the parameters. We evaluate the
empirical cumulants and compare them to the true cumulants in Table 6. We can
see that err % is small. Here the true cumulants are evaluated using (15).

When $\alpha<0$ the transition law is given in Theorem 4.3. The formula is
very simple and essentially boils down to simulating from the $\mathrm{DGGa}$
distribution, which is easily done using the approach described in Section
2.3. For several choices of the parameters, plots of the sample trajectories
of these processes using a time step of $t=1/365$ are given in Table 7 and a
comparison of the empirical and true cumulants again using $10^{5}$
observations from the transition law is given in Figure 3. We can again see
that err % is small.

|  | $c_{X,1}(0,t)$ | $c_{X,2}(0,t)$ | $c_{X,3}(0,t)$ | $c_{X,4}(0,t)$
---|---|---|---|---|---
$p$ | $\alpha$ | true | err % | true | err % | true | err % | true | err %
1.5 | $0.1$ | $0.063$ | $1.0\%$ | $0.026$ | $0.8\%$ | $0.021$ | $-1.3\%$ | $0.023$ | $-3.6\%$
$0.3$ | $0.080$ | $0.1\%$ | $0.027$ | $1.1\%$ | $0.020$ | $2.6\%$ | $0.021$ | $3.8\%$
$0.5$ | $0.113$ | $0.8\%$ | $0.029$ | $2.0\%$ | $0.019$ | $3.7\%$ | $0.020$ | $4.1\%$
$0.7$ | $0.194$ | $0.1\%$ | $0.032$ | $0.5\%$ | $0.019$ | $0.6\%$ | $0.018$ | $0.3\%$
$0.9$ | $0.610$ | $-0.3\%$ | $0.036$ | $-2.0\%$ | $0.019$ | $-4.0\%$ | $0.017$ | $-4.6\%$
2 | $0.1$ | $0.062$ | $1.0\%$ | $0.022$ | $2.2\%$ | $0.014$ | $3.7\%$ | $0.012$ | $3.8\%$
$0.3$ | $0.081$ | $1.3\%$ | $0.024$ | $2.1\%$ | $0.014$ | $1.2\%$ | $0.012$ | $-3.8\%$
$0.5$ | $0.115$ | $0.4\%$ | $0.027$ | $0.2\%$ | $0.014$ | $-0.6\%$ | $0.011$ | $-0.7\%$
$0.7$ | $0.197$ | $-0.6\%$ | $0.030$ | $-1.2\%$ | $0.015$ | $0.3\%$ | $0.011$ | $4.2\%$
$0.9$ | $0.615$ | $-0.2\%$ | $0.035$ | $-1.3\%$ | $0.015$ | $-1.2\%$ | $0.011$ | $0.8\%$
2.5 | $0.1$ | $0.063$ | $1.4\%$ | $0.021$ | $1.7\%$ | $0.012$ | $1.5\%$ | $0.009$ | $0.3\%$
$0.3$ | $0.081$ | $0.8\%$ | $0.023$ | $1.4\%$ | $0.012$ | $-0.4\%$ | $0.009$ | $-4.0\%$
$0.5$ | $0.116$ | $0.3\%$ | $0.026$ | $0.9\%$ | $0.013$ | $1.1\%$ | $0.009$ | $1.6\%$
$0.7$ | $0.199$ | $0.4\%$ | $0.030$ | $2.6\%$ | $0.013$ | $3.8\%$ | $0.009$ | $4.5\%$
$0.9$ | $0.619$ | $-0.1\%$ | $0.035$ | $-0.2\%$ | $0.014$ | $-1.2\%$ | $0.009$ | $-2.2\%$
3 | $0.1$ | $0.063$ | $-0.7\%$ | $0.020$ | $-1.6\%$ | $0.011$ | $-2.6\%$ | $0.007$ | $-3.0\%$
$0.3$ | $0.082$ | $-0.5\%$ | $0.023$ | $-1.5\%$ | $0.011$ | $-2.6\%$ | $0.007$ | $-5.3\%$
$0.5$ | $0.117$ | $0.5\%$ | $0.026$ | $0.4\%$ | $0.012$ | $0.1\%$ | $0.008$ | $0.3\%$
$0.7$ | $0.201$ | $-0.1\%$ | $0.030$ | $0.0\%$ | $0.013$ | $0.7\%$ | $0.008$ | $2.2\%$
$0.9$ | $0.621$ | $0.0\%$ | $0.035$ | $0.8\%$ | $0.014$ | $0.9\%$ | $0.008$ | $0.3\%$

Table 6: Comparison of the first four true cumulants with their estimated
values obtained from $10^{5}$ simulations from the transition law of an OUTS
process with $\lambda=10$, initial value $y_{0}=0$, and time step $t=0.1$. The
BDLP distribution is $\mathrm{TS}^{p}_{\alpha}(\lambda R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=0.1$,
$\beta=1$, and with several choices for $\alpha>0$ and $p>1$. Figure 2: Sample
trajectories of OUTS processes with parameter $\lambda=10$, initial value
$y_{0}=0$, and time step $t=1/365$. The BDLP distribution is
$\mathrm{TS}^{p}_{\alpha}(\lambda R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=0.1$,
$\beta=1$, $\alpha=0.5$, and $p\in\\{1.5,2,2.5,3\\}$.

|  | $c_{X,1}(0,t)$ | $c_{X,2}(0,t)$ | $c_{X,3}(0,t)$ | $c_{X,4}(0,t)$
---|---|---|---|---|---
$p$ | $\alpha$ | true | err % | true | err % | true | err % | true | err %
1.5 | $-0.1$ | $0.053$ | $1.17\%$ | $0.026$ | $1.73\%$ | $0.022$ | $1.89\%$ | $0.026$ | $1.32\%$
$-0.3$ | $0.046$ | $-0.16\%$ | $0.026$ | $0.34\%$ | $0.023$ | $0.00\%$ | $0.029$ | $-3.25\%$
$-0.5$ | $0.042$ | $1.28\%$ | $0.026$ | $1.41\%$ | $0.025$ | $2.43\%$ | $0.033$ | $5.27\%$
$-0.7$ | $0.040$ | $-2.02\%$ | $0.027$ | $-1.41\%$ | $0.027$ | $-0.46\%$ | $0.037$ | $0.70\%$
$-0.9$ | $0.038$ | $-0.13\%$ | $0.028$ | $-0.43\%$ | $0.030$ | $-0.16\%$ | $0.042$ | $0.87\%$
2 | $-0.1$ | $0.051$ | $0.97\%$ | $0.021$ | $2.48\%$ | $0.014$ | $5.99\%$ | $0.013$ | $6.12\%$
$-0.3$ | $0.044$ | $0.01\%$ | $0.020$ | $0.26\%$ | $0.014$ | $2.02\%$ | $0.013$ | $6.27\%$
$-0.5$ | $0.039$ | $-0.09\%$ | $0.020$ | $0.70\%$ | $0.015$ | $2.44\%$ | $0.014$ | $4.09\%$
$-0.7$ | $0.035$ | $0.78\%$ | $0.019$ | $0.55\%$ | $0.015$ | $0.81\%$ | $0.015$ | $2.08\%$
$-0.9$ | $0.033$ | $0.67\%$ | $0.019$ | $1.36\%$ | $0.016$ | $0.74\%$ | $0.016$ | $-1.95\%$
2.5 | $-0.1$ | $0.051$ | $-0.38\%$ | $0.019$ | $-0.27\%$ | $0.012$ | $1.50\%$ | $0.009$ | $6.01\%$
$-0.3$ | $0.043$ | $1.55\%$ | $0.018$ | $4.21\%$ | $0.011$ | $7.91\%$ | $0.009$ | $6.40\%$
$-0.5$ | $0.038$ | $-1.85\%$ | $0.017$ | $-3.46\%$ | $0.011$ | $-5.58\%$ | $0.009$ | $-4.02\%$
$-0.7$ | $0.034$ | $-0.67\%$ | $0.017$ | $-0.71\%$ | $0.011$ | $-1.32\%$ | $0.009$ | $-2.60\%$
$-0.9$ | $0.031$ | $-0.54\%$ | $0.016$ | $1.02\%$ | $0.011$ | $3.21\%$ | $0.010$ | $6.10\%$
3 | $-0.1$ | $0.051$ | $-0.09\%$ | $0.019$ | $0.62\%$ | $0.010$ | $1.74\%$ | $0.007$ | $4.35\%$
$-0.3$ | $0.043$ | $0.13\%$ | $0.017$ | $0.56\%$ | $0.010$ | $1.97\%$ | $0.007$ | $5.30\%$
$-0.5$ | $0.037$ | $0.10\%$ | $0.016$ | $-1.69\%$ | $0.010$ | $-4.80\%$ | $0.007$ | $-4.52\%$
$-0.7$ | $0.033$ | $1.02\%$ | $0.015$ | $1.22\%$ | $0.010$ | $2.16\%$ | $0.007$ | $4.73\%$
$-0.9$ | $0.030$ | $-0.19\%$ | $0.015$ | $-1.00\%$ | $0.009$ | $-2.36\%$ | $0.007$ | $-3.78\%$

Table 7: Comparison of the first four true cumulants with their estimated
values obtained from $10^{5}$ simulations from the transition law of an OUTS
process with $\lambda=10$, initial value $y_{0}=0$, and time step $t=0.1$. The
BDLP distribution is $\mathrm{TS}^{p}_{\alpha}(\lambda R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=0.1$,
$\beta=1$, and with several choices for $\alpha<0$ and $p>1$.

Comparison of the first four true cumulants with their estimated values
obtained from $10^{5}$ simulations from the transition law of an OUTS process
with $\lambda=10$, initial value $y_{0}=0$, and time step $t=0.1$. The
stationary distribution is $\mathrm{TS}^{p}_{\alpha}(R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=1$,
$\beta=1$, and several choices for $\alpha$ and $p$.

Figure 3: Sample trajectories of OUTS processes with parameter $\lambda=10$,
initial value $y_{0}=0$, and time step $t=1/365$. The BDLP distribution is
$\mathrm{TS}^{p}_{\alpha}(\lambda R,0)$ with
$R(\mathrm{d}x)=c\beta^{\alpha}\delta_{1/\beta}(\mathrm{d}x)$, where $c=0.1$,
$\beta=1$, $\alpha=-0.5$, and $p\in\\{1.5,2,2.5,3\\}$.

##  6 Proofs

In this section we give the proofs.

###### Proof of Lemma 2.2.

Note that

$\displaystyle f_{\beta,\gamma,p,\eta}(u)$ $\displaystyle=$
$\displaystyle\frac{u^{-1-\beta}}{K_{\beta,\gamma,p,\eta}\Gamma(\gamma)}\int_{0}^{u^{p}(\eta-1)}x^{\gamma-1}e^{-x-u^{p}}\mathrm{d}x$
$\displaystyle=$
$\displaystyle\frac{u^{p\gamma-\beta-1}p}{K_{\beta,\gamma,p,\eta}\Gamma(\gamma)}\int_{1}^{\eta^{1/p}}(\theta^{p}-1)^{\gamma-1}e^{-\theta^{p}u^{p}}\theta^{p-1}\mathrm{d}\theta$
$\displaystyle=$
$\displaystyle\int_{1}^{\eta^{1/p}}\left(p\frac{\theta^{p\gamma-\beta}u^{p\gamma-\beta-1}}{\Gamma(\gamma-\beta/p)}e^{-(\theta
u)^{p}}\right)\left(\frac{p}{K^{*}_{\beta,\gamma,p,\eta}}(\theta^{p}-1)^{\gamma-1}\theta^{p+\beta-p\gamma-1}\right)\mathrm{d}\theta,$

where the second line follows by the change of variables
$\theta^{p}=(xu^{-p}+1)$. ∎

###### Proof of Proposition 2.3..

First, let $Y\sim\mathrm{GGa}(p\gamma-\beta,p,\theta^{p})$ and note that

$\mathrm{E}[Y^{\xi}]=\frac{\Gamma\left(\gamma+(\xi-\beta)/p\right)}{\Gamma\left(\gamma-\beta/p\right)}\theta^{-\xi}.$

From here, by a conditioning argument, we have

$\displaystyle\mathrm{E}[W^{\xi}]$
$\displaystyle=\frac{\Gamma\left(\gamma+\xi-\beta)/p\right)}{\Gamma\left(\gamma-\beta/p\right)}\int_{1}^{\eta}\theta^{-\xi}m^{\sharp}_{\beta,\gamma,p,\eta}(\theta)\mathrm{d}\theta$
$\displaystyle=\frac{\Gamma\left(\gamma+(\xi-\beta)/p\right)}{\Gamma\left(\gamma-\beta/p\right)C^{*}_{\beta,\gamma,p,\eta}}\sum_{k=0}^{\gamma-1}{\gamma-1\choose
k}\frac{(-1)^{k}}{pk-\beta}\int_{1}^{\eta}\theta^{-\xi-1}\left(1-\theta^{\beta-
pk}\right)\mathrm{d}\theta,$ (19)

where

$\int_{1}^{\eta}\theta^{-\xi-1}\left(1-\theta^{\beta-
pk}\right)\mathrm{d}\theta=\frac{1-\eta^{-\xi}}{\xi}+\frac{1-\eta^{\beta-
pk-\xi}}{\beta-pk-\xi}.$

Next, using the fact that

$\int_{1}^{\eta}\theta^{-\xi-1}\ln\theta\mathrm{d}\theta=\frac{1-\eta^{-\xi}(\xi\ln\eta+1)}{\xi^{2}},$

the result can be proved in a similar way if $\beta=pk$ for some
$k\in\\{0,1,2,\dots,\gamma-1\\}$. ∎

Theorem 4.2 is an immediately consequence of the following lemma.

######  Lemma 6.1.

In the context of Theorem 4.2, $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

$\displaystyle C_{t}(y,z)$ $\displaystyle=$ $\displaystyle ie^{-\lambda
t}\langle y,z\rangle+i\left(1-e^{-\lambda t}\right)\langle
b,z\rangle-\sum_{n=0}^{\gamma-1}i\langle b_{n},z\rangle$
$\displaystyle\quad+\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{\alpha}(e^{-\lambda
t}z,xv)e^{-v^{p}}v^{-1-\alpha}\mathrm{d}vR^{*}_{0}(\mathrm{d}x)$
$\displaystyle\quad+\sum_{n=1}^{\gamma-1}\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{\alpha-
np}(e^{-\lambda
t}z,xv)e^{-v^{p}}v^{np-1-\alpha}\mathrm{d}vR^{*}_{n}(\mathrm{d}x)$
$\displaystyle\quad+\frac{pC_{\alpha,\gamma,p,e^{\lambda
t}}R(\mathbb{R}^{d})}{(\gamma-1)!}\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{0}(z,xv)f^{\sharp}_{\alpha,\gamma,p,e^{\lambda
t}}(v)\mathrm{d}vR^{1}(\mathrm{d}x)$

and

$\displaystyle\psi_{\alpha}(z,x)=e^{i\langle z,x\rangle}-1-i\langle z,x\rangle
1_{[\alpha\geq 1]}.$

###### Proof.

Proposition 2.13 in [23] implies that

$\displaystyle C_{t}(y,z)=ie^{-\lambda t}\langle
y,z\rangle+i\left(1-e^{-\lambda t}\right)\langle
b,z\rangle+\lambda\int_{0}^{t}\int_{\mathbb{R}^{d}}\psi_{\alpha}(e^{-\lambda
s}z,x)M(\mathrm{d}x)\mathrm{d}s,$

where $M$ is the Lévy measue of $\mathrm{TS}^{p}_{\alpha}(R,b)$. Now using
(12) and the fact that $\psi_{\alpha}(az,x)=\psi_{\alpha}(z,ax)$ for any
$a\in\mathbb{R}$ gives

$\displaystyle\lambda\int_{0}^{t}\int_{\mathbb{R}^{d}}\psi_{\alpha}(e^{-\lambda
s}z,x)M(\mathrm{d}x)\mathrm{d}s$
$\displaystyle\qquad=\lambda\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\int_{0}^{t}\psi_{\alpha}(z,xue^{-\lambda
s})u^{-1-\alpha}e^{-u^{p}}\mathrm{d}s\mathrm{d}uR(\mathrm{d}x)$
$\displaystyle\qquad=\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\int_{ue^{-\lambda
t}}^{u}\psi_{\alpha}(z,xv)u^{-1-\alpha}e^{-u^{p}}v^{-1}\mathrm{d}v\mathrm{d}uR(\mathrm{d}x)$
$\displaystyle\qquad=\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{\alpha}(z,xv)v^{-1}\int_{v}^{ve^{\lambda
t}}u^{-1-\alpha}e^{-u^{p}}\mathrm{d}u\mathrm{d}vR(\mathrm{d}x).$ (20)

Note that for $v>0$ we have

$\displaystyle v^{-1}\int_{v}^{ve^{\lambda
t}}u^{-1-\alpha}e^{-u^{p}}\mathrm{d}u=\int_{1}^{e^{\lambda
t}}(uv)^{-1-\alpha}e^{-(vu)^{p}}\mathrm{d}u$
$\displaystyle\quad=e^{-v^{p}e^{\lambda tp}}\int_{1}^{e^{\lambda
t}}(uv)^{-1-\alpha}\mathrm{d}u+\int_{1}^{e^{\lambda
t}}\left(e^{-(vu)^{p}}-e^{-v^{p}e^{\lambda
tp}}\right)(uv)^{-1-\alpha}\mathrm{d}u$
$\displaystyle\quad=\frac{e^{-v^{p}e^{\lambda
tp}}}{\alpha}\left(1-{e^{-\alpha\lambda
t}}\right)v^{-1-\alpha}+e^{-v^{p}e^{\lambda
tp}}\sum_{n=1}^{\gamma-1}v^{np-1-\alpha}\int_{1}^{e^{\lambda
t}}\frac{(e^{\lambda tp}-u^{p})^{n}}{n!}u^{-1-\alpha}\mathrm{d}u$
$\displaystyle\qquad+\int_{1}^{e^{\lambda
t}}\left(e^{-(vu)^{p}}-e^{-v^{p}e^{\lambda
tp}}\sum_{n=0}^{\gamma-1}\frac{(e^{\lambda
tp}-u^{p})^{n}}{n!}v^{np}\right)(uv)^{-1-\alpha}\mathrm{d}u.$

Now applying Lemma 1 in [13]

$\displaystyle\int_{1}^{e^{\lambda t}}\left(e^{-(vu)^{p}}-e^{-v^{p}e^{\lambda
tp}}\sum_{n=0}^{\gamma-1}\frac{(e^{\lambda
tp}-u^{p})^{n}}{n!}v^{np}\right)(uv)^{-1-\alpha}\mathrm{d}u$
$\displaystyle\qquad=\frac{1}{(\gamma-1)!}\int_{1}^{e^{\lambda
t}}\int_{0}^{v^{p}(e^{\lambda
tp}-u^{p})}e^{-(x+v^{p}u^{p})}x^{\gamma-1}\mathrm{d}x(uv)^{-1-\alpha}\mathrm{d}u$
$\displaystyle\qquad=\frac{pv^{p\gamma-\alpha-1}}{(\gamma-1)!}\int_{1}^{e^{\lambda
t}}\int_{u}^{e^{\lambda
t}}e^{-v^{p}y^{p}}(y^{p}-u^{p})^{\gamma-1}y^{p-1}\mathrm{d}yu^{-1-\alpha}\mathrm{d}u$
$\displaystyle\qquad=\frac{pv^{p\gamma-\alpha-1}}{(\gamma-1)!}\int_{1}^{e^{\lambda
t}}y^{p-1}e^{-v^{p}y^{p}}\int_{1}^{y}(y^{p}-u^{p})^{\gamma-1}u^{-1-\alpha}\mathrm{d}u\mathrm{d}y$
$\displaystyle\qquad=\frac{pv^{p\gamma-\alpha-1}}{(\gamma-1)!}\int_{1}^{e^{\lambda
t}}y^{p\gamma-\alpha-1}e^{-v^{p}y^{p}}\int_{1/y}^{1}(1-s^{p})^{\gamma-1}s^{-1-\alpha}\mathrm{d}s\mathrm{d}y$
$\displaystyle\qquad=\frac{pC_{\alpha,\gamma,p,e^{\lambda
t}}}{(\gamma-1)!}f^{\sharp}_{\alpha,\gamma,p,e^{\lambda t}}(v),$

where the third line follows by the substitution $y^{p}=v^{-p}x+u^{p}$ and the
fifth by the substitution $s=u/y$. Now putting everything together and using
the readily checked facts that

$\int_{0}^{\infty}\psi_{\alpha}(z,xv)e^{-v^{p}e^{\lambda
tp}}v^{np-1-\alpha}\mathrm{d}v=e^{\lambda t(\alpha-
np)}\int_{0}^{\infty}\psi_{\alpha}(z,e^{-\lambda
t}xv)e^{-v^{p}}v^{np-1-\alpha}\mathrm{d}v$

and

$\int_{1}^{e^{\lambda t}}\frac{(e^{\lambda
tp}-u^{p})^{n}}{n!}u^{-1-\alpha}\mathrm{d}u=e^{\lambda
t(pn-\alpha)}\int^{1}_{e^{-\lambda
t}}\frac{(1-u^{p})^{n}}{n!}u^{-1-\alpha}\mathrm{d}u$

gives the result. ∎

###### Proof of Theorem 4.3.

Following the proof of Lemma 6.1 to (20) shows that the characteristic
function of the temporally homogenous transition function
$P_{t}(y,\mathrm{d}x)$ is given by $\int_{\mathbb{R}^{d}}e^{i\langle
x,z\rangle}P_{t}(y,\mathrm{d}x)=\exp\left\\{C_{t}(y,z)\right\\}$, where

$\displaystyle C_{t}(y,z)$ $\displaystyle=$ $\displaystyle ie^{-\lambda
t}\langle y,z\rangle+i\left(1-e^{-\lambda t}\right)\langle b,z\rangle$
$\displaystyle\quad+\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{\alpha}(z,xv)v^{-1}\int_{v}^{ve^{\lambda
t}}u^{-1-\alpha}e^{-u^{p}}\mathrm{d}u\mathrm{d}vR(\mathrm{d}x).$

Noting that the formula on the second line equals

$p^{-1}\lambda
t\Gamma(|\alpha|/p)R(\mathbb{R}^{d})\int_{\mathbb{R}^{d}}\int_{0}^{\infty}\psi_{\alpha}(z,xv)h_{|\alpha|,p,e^{\lambda
t}}(v)\mathrm{d}vR^{1}(\mathrm{d}x)$

gives the result. ∎

## References

  * [1] O.E. Barndorff-Nielsen and N. Shephard.  Non-Gaussian Ornstein-Uhlenbeck-based models and some of their uses in financial economics.  Journal of the Royal Statistical Society: Series B, 63(2):167–241, 2001.
  * [2] F.E. Benth, T. Meyer-Brandis, and J. Kallsen.  A non-Gaussian Ornstein-Uhlenbeck process for electricity spot price modeling and derivatives pricing.  Applied Mathematical Finance, 14(2): 153–169, 2007.
  * [3] A. Bignami and A. de Matteis.  A note on sampling from combination of distribution.  Journal of the Institute of Mathematics and its Applications, 8:80-81, 1971.
  * [4] R. Cont and P. Tankov.  Financial Modelling with Jump Processes.  Chapman and Hall, London, 2004.
  * [5] N. Cufaro Petroni and P. Sabino.  Fast pricing of energy derivatives with mean-reverting jump-diffusion processes.  Applied Mathematical Finance, 0(0):1–22, 2021.
  * [6] L. Devroye.  Non-Uniform Random Variate Generation.  Springer-Verlag, New York, 1986
  * [7] L. Devroye.  Random variate generation for exponentially and polynomially tilted stable distributions.  ACM Transactions on Modeling and Computer Simulation, 19(4): Article 18, 2009.
  * [8] M. Grabchak.  On a new class of tempered stable distributions: Moments and regular variation.  Journal of Applied Probability, 49(4):1015–1035, 2012.
  * [9] M. Grabchak.  Tempered stable distributions.  Springer, Cham, Switzerland, 2016.
  * [10] M. Grabchak.  Rejection sampling for tempered Lévy processes.  Statistics and Computing, 29(3):549–558, 2019.
  * [11] M. Grabchak.  On the simulation of general tempered stable Ornstein–Uhlenbeck processes.  Journal of Statistical Computation and Simulation, 90(6):1057–1081, 2020.
  * [12] M. Grabchak.  An exact method for simulating rapidly decreasing tempered stable distributions.  Statistics & Probability Letters, 170:109015, 2021.
  * [13] M. Grabchak.  On the transition laws of $p$-tempered $\alpha$-stable OU-processes.  Computational Statistics, 36:1415–1436, 2021.
  * [14] M. Grabchak and L. Cao. SymTS: Symmetric tempered stable distributions. Ver. 1.0, R Package. https://cran.r-project.org/web/packages/SymTS/index.html, 2017.
  * [15] M. Grabchak and G. Samorodnitsky.  Do financial returns have finite or infinite variance? A paradox and an explanation.  Quantitative Finance, 10(8):883–893, 2010.
  * [16] M. Hofert.  Sampling exponentially tilted stable distributions.  ACM Transactions on Modeling and Computer Simulation, 22(1): Article 3, 2011.
  * [17] J. Imai an R. Kawai.  On finite truncation of infinite shot noise series representation of tempered stable laws. Physica A, 390(23–24), 4411–4425, 2011.
  * [18] R. Kawai and H. Masuda.  On simulation of tempered stable random variates.  Journal of Computational and Applied Mathematics, 235(8):2873–2887, 2011.
  * [19] R. Kawai and H. Masuda.  Infinite variation tempered stable Ornstein–Uhlenbeck processes with discrete observations.  Communications in Statistics – Simulation and Computation, 41(1):125–139, 2012.
  * [20] Y.S. Kim M.L. Bianchi, S.T. Rachev and F.J. Fabozzi.  Tempered infinitely divisible distributions and processes.  Theory of Probability and Its Applications, 55(1):2–26, 2011.
  * [21] Y. Qu, A. Dassios, and H. Zhao.  Exact simulation of gamma-driven Ornstein–Uhlenbeck processes with finite and infinite activity jumps.  Journal of the Operational Research Society, 72(2):471–484, 2021.
  * [22] Y. Qu, A. Dassios, and H. Zhao.  Exact simulation of Ornstein–Uhlenbeck tempered stable processes.  Journal of Applied Probability, 58(2):347–371, 2021.
  * [23] A. Rocha-Arteaga and K. Sato (2019).  Topics in Infinitely Divisible Distributions and Lévy Processes, Revised Edition.  Springer, Cham.
  * [24] J. Rosiński.  Tempering stable proceses.  Stochastic Processes and their Applications, 117(6):677–707, 2007\.
  * [25] J. Rosiński and J.L. Sinclair.  Generalized tempered stable processes.  Banach Center Publications, 90:153–170, 2010.
  * [26] P. Sabino.  Exact simulation of variance gamma related OU proceses: Application to the pricing of energy derivatives.  Applied Mathematical Finance, 27(3):207–227, 2020.
  * [27] P. Sabino and N. Cufaro Petroni.  Gamma-related Ornstein–Uhlenbeck processes and their simulation*.  Journal of Statistical Computation and Simulation, 91(6):1108-1133, 2021.
  * [28] P. Sabino and N. Cufaro Petroni.  Fast simulation of tempered stable Ornstein-Uhlenbeck processes.  To appear in Computational Statistics.
  * [29] K. Sato.  Lévy Processes and Infinitely Divisible Distributions.  Cambridge University Press, Cambridge, 1999.
  * [30] E.W. Stacy. A generalization of the gamma distribution. The Annals of Mathematical Statistics, 33(3):1187–1192, 1962.
  * [31] Y. Xia and M. Grabchak. Estimation and simulation for multivariate tempered stable distributions. Journal of Statistical Computation and Simulation, 92(3):251-275, 2022
  * [32] S. Zhang.  Exact simulation of tempered stable Ornstein–Uhlenbeck proceses.  Journal of Statistical Computation and Simulation, 81(11):1533–1544, 2011.
  * [33] S. Zhang, Z. Sheng, and W. Deng. On the transition law of O-U compound Poisson processes. In 2011 Fourth International Conference on Information and Computing, pp. 260–263, 2011.
  * [34] S. Zhang and X. Zhang. Exact simulation of IG-OU processes. Methodology and Computing in Applied Probability, 10(3):337–355, 2008.
  * [35] S. Zhang and X. Zhang (2009).  On the transition law of tempered stable Ornstein-Uhlenbeck processes.  Journal of Applied Probability, 46(3):721–731.