1 Superpositions of Ornstein–Uhlenbeck processes
We consider the superpositions of Ornstein–Uhlenbeck (supOU) processes which are an important class of ambit stochastic processes [11] and play a fundamental role in finance and insurance, degradation modelling, diffusion, turbulence, astrophysics and transport, see [20, 52, 30, 51] and references therein. SupOU processes are proving to be invaluable tools for modeling complex phenomena with long-range dependence, intermittency, and non-Gaussian behaviors.
In finance, the supOU processes can be incorporated in different ways. The popular continuous time stochastic volatility model for financial assets
generalizes the classical Black–Scholes model by replacing the volatility parameter σ by the volatility process $\sigma (t)$ in the form of a superposition of positive OU processes, where $B(t)$ is the standard Brownian motion that is independent of $\sigma (t)$, see the fundamental work [8]. If $\sigma (t)$ is a supOU process with inverse gamma marginals, then log-returns have approximately the Student distribution, see [8, p. 170] for a general correspondence.
The fractal activity time geometric Brownian motion model introduced in [35] has the form
where $S(t)$ is the asset price and the fractal activity time $T(t)$ is a positive, nondecreasing stochastic process in the form of the linear spline of a discrete integrated supOU process with positive marginals that is independent of $B(t)$. If $T(t)$ is an integrated supOU process with inverse gamma marginals, then log-returns have exactly the Student distribution, see [27, 28] for other marginals.
For modeling high-frequency financial data, supOU processes were explored in [37] by demonstrating their efficacy in capturing the stochastic volatility of stock prices with improved accuracy over traditional models like GARCH. The authors leveraged the superposition property to model multiscale temporal dependencies, which proved particularly effective for intraday trading strategies. Furthermore, supOU processes were applied to stock indexes [36], zero-coupon bonds [61], and option pricing [48, 59].
SupOU processes are well suited for modeling migrating fish counts [62] and intricate hydrological phenomena, such as river discharge time series [65, 63] and water characteristics [64]. In the context of wind energy systems, the utility of graph supOU models in capturing the spatial and temporal variability of wind capacity factors across a European electricity network was demonstrated in [50]. In astrophysics, the supOU processes were applied to modelling the X-ray light curves generated by black holes [42].
For the main properties of supOU processes, we refer to several works by Barndorff-Nielsen and his collaborators [5, 9, 14, 6]. In particular, for every self-decomposable (SD) distribution, there is a stationary supOU process with this SD distribution as marginal distribution and quite flexible class of covariance structures including long-range dependence (LRD), that follows from the following theorem.
Theorem 1.
For any SD random variable V there is a double-sided Lévy process $Z(t)$, called backward driving Lévy process (BDLP), such that
where we assume $\mathbb{E}[1+\log |Z(1)|]\lt \infty $ and $\log \mathbb{E}{e^{izZ(t)}}=t\hspace{0.1667em}{\kappa _{Z}}(z)$, where
is the cumulant function of an infinitely divisible (ID) random variable $Z(1)$ with Lévy–Khintchine characteristic triplet $(b,{\sigma ^{2}},{\mu _{Z}})$, that is,
where
We note that Theorem 1 is a reformulation of Theorem 3.9.3 in [41]. Let us now give the definition of supOU processes.
Definition 1.
A supOU process $Y(t),\hspace{2.83862pt}t\in \mathbb{R}$, is a strictly stationary process defined by the stochastic integral
in the sense of the paper [53], where Λ is a homogeneous infinitely divisible independently scattered random measure (Lévy basis) on ${\mathbb{R}_{+}}\times \mathbb{R}$ with cumulant function
\[ \log \mathbb{E}{e^{iz\Lambda (A)}}={m_{c}}(A){\kappa _{Z}}(z)=(\pi \times \mathit{Leb})(A){\kappa _{Z}}(z),\hspace{1em}z\in \mathbb{R},\]
for any Borel set $A\in \mathcal{B}({\mathbb{R}_{+}}\times \mathbb{R})$, where the control measure ${m_{c}}=\pi \times \mathit{Leb}$ is the product of a measure π on $(0,\infty )$, such that
and the Lebesgue measure $\mathit{Leb}$ on $\mathbb{R}$, the Lévy measure ${\mu _{Z}}$ of the cumulant function ${\kappa _{Z}}(z)$ is a positive Radon measure on $\mathbb{R}\setminus \{0\}$ such that
The quadruple $(b,{\sigma ^{2}},{\mu _{Z}},\pi )$ is referred to as the characteristic quadruple and it completely determines the supOU process (1) with any given marginal SD distribution with Lévy triplet $({b_{Y}},{\sigma _{Y}^{2}},{\mu _{Y}})$, that is,
Recall that the ID random variable $V\stackrel{d}{=}Y(0)$ with characteristic function $\phi (z)=\exp ({\kappa _{Y}}(z))$ is self-decomposable if for every constant $c\in (0,1)$ there is a characteristic function ${\phi _{c}}(z)$ such that $\phi (z)=\phi (cz){\phi _{c}}(z)$, $z\in \mathbb{R}$. The ID random variable V is self-decomposable if and only if its Lévy measure has the form ${\mu _{Y}}(dx)=\frac{g(x)}{|x|}dx$, where $g(x)$ is increasing in $(-\infty ,0)$ and decreasing on $(0,\infty )$, see [58, 54]. We also note that ${\kappa _{Z}}(z)=z\frac{d}{dz}{\kappa _{Y}}(z)$.
(3)
\[ {\kappa _{Y}}(z)=iz{b_{Y}}-\frac{1}{2}{z^{2}}{\sigma _{Y}^{2}}+{\int _{\mathbb{R}\setminus \{0\}}}({e^{izx}}-1-izx{1_{[-1,1]}}(x)){\mu _{Y}}(dx).\]If $\mathbb{E}{Y^{2}}(t)\lt \infty $, the correlation function of the supOU process (1) is given by
which can exhibit either short-range dependence (SRD) or LRD depending on the measure π.
For explanation of the concept of the superposition, let us introduce the Lévy-driven Ornstein–Uhlenbeck process
where $Z(t)$ is a two-sided Lévy process with Lévy–Khintchine triplet $(b,{\sigma ^{2}},{\mu _{Z}})$.
The correlation function of the OU process $X(t)$ is of the form
provided $\mathbb{E}{X^{2}}(t)\lt \infty $.
Thus, if π is a probability measure, then the measure π provides a randomization of the rate parameter ξ of the Lévy-driven OU process $X(t)$ and hence the correlation structure (4) of the supOU process is more flexible and can exhibit LRD. The popular example of the measure π is the gamma distribution, some other examples can be found in [7]. Although the OU process $X(t)$ is a linear process, the supOU process is nonlinear in general, that is, $Y(t)$ cannot be represented in the form ${\textstyle\int _{-\infty }^{t}}a(t-s)dZ(s)$ for some function $a(\cdot )$, see [23, 5, 47].
In [38–40] self-decomposability was shown for the gamma distribution, the inverse gamma distribution, the inverse Gaussian distribution, the Student distribution, the hyperbolic cosine distribution, the Gumbel distribution and the Bessel distribution. Also, self-decomposability was shown for the generalized inverse Gaussian distribution in [34], for the generalized Gaussian distribution in [24], for the positive α-stable distribution in [9], and for the generalized Linnik distribution in [3].
Self-decomposability of the log-normal distribution follows from [18], where it is proven that the log-normal distribution belongs to the class of generalized gamma convolutions, which is a subclass of SD distributions. Finally, self-decomposability of the Rosenblatt distribution has been proven in [49].
Let us present the fundamental result which explains the structure of supOU processes.
Theorem 2.
Let $Y(t)$ be a supOU process with generating quadruple $(b,0,{\mu _{Z}},\pi )$, where the measure π is arbitrary and satisfies (2). Then there exists a modification of Λ such that for $A\in \mathcal{B}({\mathbb{R}_{+}}\times \mathbb{R})$ we have
where $-\infty \lt \cdots \lt {S_{-1}}\lt {S_{0}}\le 0\lt {S_{1}}\lt \cdots \lt \infty $ are the jump times of a two-sided Poisson process on $\mathbb{R}$ with some intensity θ, $\{{R_{k}},\hspace{0.1667em}k\in \mathbb{Z}\}$ are i.i.d. variables with the distribution π independent of $\{{S_{k}},\hspace{2.83862pt}k\in \mathbb{Z}\}$, and $\{{H_{k}},\hspace{0.1667em}k\in \mathbb{Z}\}$ are i.i.d. variables with the shot-height distribution F such that
and equivalently
where ${\mu _{Y}}$ is the Lévy measure for the supOU process $Y(t)$. Thus, the supOU process can be represented as
Proof.
From the Lévy–Itô decomposition, the Lévy basis Λ in (1) can be written as
\[ \Lambda (d\xi ,ds)=b\pi (d\xi )ds+{\int _{(0,1]}}x(\nu -{\mu _{Z}})(d\xi ,ds,dx)+{\int _{(1,\infty )}}x\nu (d\xi ,ds,dx),\]
where $b\in \mathbb{R}$ and ν is a Poisson random measure on ${\mathbb{R}_{+}}\times \mathbb{R}\times {\mathbb{R}_{+}}$ with intensity measure $\pi (d\xi )ds{\mu _{Z}}(dx)$. By the Poisson construction theorem, for the Poisson random measure ν there exists a sequence $\{({H_{k}},{R_{k}},{S_{k}})\}$ such that the equality (5) holds. Substituting the Lévy basis (5) into the formula (1), we obtain the representation (7). Direct calculation of the cumulant function of the process (7) shows that it coincides with the cumulant function of the process (1). We refer to [56] for further details, see also [10, Theorem 2.2] and [26]. □The representation of ν in the proof of Theorem 6 is most convenient for us among several series representations of random measures which are reviewed in [56].
The most famous example of the application of Theorem 2 is the gamma process, however, the relation (6) does not work for many other marginal distributions, see Section 2. We note that the equation (7) is called the Bondesson–Rosinski representation due to [16, 56], which is not unique, if we change the intensity θ, then the shot-height distribution $F(x)$ should be scaled. If ${R_{k}}=\lambda $ for all k, then the formula (7) becomes the representation of the usual OU process. Since ${S_{k}}$ are jump times of the Poisson process with intensity θ, we have
that is, ${S_{k}}-{S_{k-1}}$ has the exponential distribution with rate θ (implying mean $1/\theta )$. Thus, the Bondesson–Rosinski representation (7) has in average θ summands per any unit interval. We note that ${S_{k}}/\eta $ are jumps of the Poisson process with intensity $\theta /\eta $. In Figure 1 we show a typical sample path of the supOU process $Y(t)$, which consists of jumps with different heights and mixed-exponential decays with different rates.
Overall, we can see that the supOU processes form a wide class of stochastic processes and the representation (7) enables to independently specify the marginal distribution by the choice of the shot-height distribution and the correlation structure by the choice of the measure π.
We note that the Lévy measure ${\mu _{Y}}$ for the specified marginal distribution of the process $Y(t)$ can be derived from the one-to-one correspondence between the Lévy density and the characteristic function as stated in the following result.
Theorem 3.
Let ${\kappa _{Y}}(z)=\log \phi (z)$ be the cumulant function of an ID marginal distribution of the supOU process $Y(t)$ with Lévy triplet $({b_{Y}},0,{\mu _{Y}})$ and ${\textstyle\int _{|x|\gt 1}}\hspace{-0.1667em}{x^{2}}{\mu _{Y}}(dx)\hspace{-0.1667em}\lt \hspace{-0.1667em}\infty $. Then the Lévy measure ${\mu _{Y}}(dx)=m(x)dx$ satisfies the equation
If ${\textstyle\int _{-\infty }^{\infty }}|{(\log \phi (z))^{\prime\prime }}|dz\lt \infty $, then the Lévy measure ${\mu _{Y}}$ is absolutely continuous with respect to the Lebesgue measure.
Proof.
By differentiating twice the equality (3) with ${\sigma _{Y}^{2}}=0$, we obtain the equation (8), which can be interpreted as the Fourier transform of ${x^{2}}m(x)$. For its inversion, we can use the identity
where we refer to [19] for further details. □
\[ {\int _{B}}{x^{2}}m(x)dx=-\frac{1}{2\pi }v.p.{\int _{-\infty }^{\infty }}{(\log \phi (z))^{\prime\prime }}\left({\int _{B}}{e^{-ixz}}dx\right)dz,\]
where B is any Borel set. If ${\textstyle\int _{-\infty }^{\infty }}|{(\log \phi (z))^{\prime\prime }}|dz\lt \infty $, then
(9)
\[ m(x)=-\frac{1}{2\pi {x^{2}}}{\int _{-\infty }^{\infty }}{e^{-ixz}}{(\log \phi (z))^{\prime\prime }}dz,\]Let us now reformulate Theorem 3 for the supOU processes $Y(t)$ with positive marginals using the Laplace transform $\psi (s)=\mathbb{E}{e^{-sY(t)}}={\textstyle\int _{0}^{\infty }}{e^{-sx}}p(x)dx$, $s\gt 0$, where $p(x)$ is the density of the marginal distribution. Following [17, p. 15], the density $p(x)$ is infinitely divisible if and only if
where $a\ge 0$ is the left-extremity of the distribution. Taking the first and second derivatives of the above identity with $a=0$, we obtain
which can be interpreted as the Laplace transforms of $xm(x)$ and ${x^{2}}m(x)$, respectively. Therefore, the Lévy density $m(x)$ can be found by the inversion of either (10) or (11) using one of several methods: the Bromwich integral, the Post–Widder formula, or Fredholm equations, see [21] for details. In particular, the Lévy density for the log-normal distribution was numerically computed by the method of Fredholm equations in [18]. In our numerical examples, where the Lévy density is not known analytically, we use the method based on the Bromwich integral.
The following theorem establishes another useful relation for the OU-type processes, one that does do not involve the Lévy measure and is motivated by the theory of shot-noise processes [58].
Theorem 4.
Define a shot-noise process
\[ {Y_{h}}(t)={\sum \limits_{k=-\infty }^{\infty }}{H_{k}}h(t-{S_{k}}){1_{[0,\infty )}}(t-{S_{k}}),\]
where $h(\cdot )$ is an impulse response, ${S_{k}}$ are jumps of the Poisson process with intensity θ and ${H_{k}}$ are i.i.d. random variables with cdf $F(x)$, $x\in \mathbb{R}$. Then the characteristic function $\phi (z)$ of the marginal distribution of ${Y_{h}}(t)$ can be written as
Proof.
Since $\{{S_{k}}\}$ are the points of a Poisson process on $\mathbb{R}$ with intensity θ, the marked point process $\{({S_{k}},{H_{k}})\}$ is a Poisson random measure on $\mathbb{R}\times \mathbb{R}$ with intensity measure
The characteristic function of ${Y_{h}}(t)$ is
\[ \phi (z)=\mathbb{E}\big[{e^{iz{Y_{h}}(t)}}\big]=\mathbb{E}\Big[\exp \Big(iz\sum \limits_{k}{H_{k}}h(t-{S_{k}}){1_{[0,\infty )}}(t-{S_{k}})\Big)\Big].\]
By the Laplace functional of a Poisson random measure, we have
\[ \mathbb{E}\Big[\exp \Big(\sum f({S_{k}},{H_{k}})\Big)\Big]=\exp \Big(\int ({e^{f(s,x)}}-1)\hspace{0.1667em}\nu (ds,dx)\Big)\]
for any measurable function f. Taking $f(s,x)=izxh(t-s){1_{[0,\infty )}}(t-s)$, we obtain
\[ \log \phi (z)={\int _{\mathbb{R}}}{\int _{\mathbb{R}}}\big({e^{izxh(t-s){1_{[0,\infty )}}(t-s)}}-1\big)\hspace{0.1667em}\nu (ds,dx).\]
Since $\nu (ds,dx)=\theta \hspace{0.1667em}ds\hspace{0.1667em}dF(x)$ and ${1_{[0,\infty )}}(t-s)$ restricts $s\le t$, we have
By the change of variables, we obtain the required statement. □From Theorem 4 for the exponential impulse response $h(t)={e^{-Rt}}$, we have
1.1 Choice of the correlation function via the measure π
From the relation (4) we have
and this integral can be either finite or infinite. Hence, we will say that a supOU process exhibits LRD (long memory) if the integral in (12) is infinite, and we will say that it exhibits SRD (short memory) otherwise. Moreover, if π is regularly varying at zero, i.e.
for some $\alpha \gt 0$ and a slowly varying function $\ell (\cdot )$ at infinity, then
(13)
\[ \pi ((0,\xi ])\sim \ell ({\xi ^{-1}}){\xi ^{1+\alpha }},\hspace{1em}\hspace{2.5pt}\text{as}\hspace{2.5pt}\xi \to 0,\]
\[ r(\tau )\sim \frac{\eta (1+\alpha )}{\Gamma (\alpha )}\ell (\tau ){\tau ^{-\alpha }},\hspace{1em}\hspace{2.5pt}\text{as}\hspace{2.5pt}\tau \to \infty ,\]
see [30]. In particular, if $\alpha \in (0,1)$ in (13), then the supOU process exhibits LRD; see [26] and [31] for details.We consider two cases of the measure π in the present paper. If π is the Dirac measure at point λ, then the supOU process is actually a usual OU process. If the measure π is given by the gamma distribution $\Gamma (1+\alpha ,\beta )$ with density
\[ {p_{\Gamma (1+\alpha ,\beta )}}(\xi )=\frac{{\beta ^{1+\alpha }}}{\Gamma (1+\alpha )}{\xi ^{\alpha }}{e^{-\beta \xi }},\hspace{1em}\xi \gt 0,\hspace{3.33333pt}\alpha \gt 0,\hspace{3.33333pt}\beta \gt 0,\]
where $\Gamma (\cdot )$ is the gamma function, then (2) holds and the supOU process has the correlation function
In particular, for $\alpha \in (0,1)$ we obtain processes with LRD. More examples of possible choices of π and corresponding supOU processes are given in [7].1.2 Choice of the marginal distribution via the shot-height distribution $F(x)$
Due to (6), we can derive the shot-height distribution $F(x)$ via the Lévy measure if the Lévy measure of the marginal distribution is known. However, most of the known simulation algorithms for OU processes with GIG marginals are constructed using the autoregressive formulae of order 1 and are rather difficult, see [66]. We also note that we cannot extend the simulation algorithm for OU processes from [66] to the case of superpositions because a linear combination of OU processes does not have the desired fat-tail marginal distribution.
The gamma OU process is the most famous stochastic process which has a very simple simulation algorithm obtained in [16]. The gamma distribution $\Gamma (a,b)$ has the density $p(x)={b^{a}}{x^{a-1}}{e^{-bx}}/\Gamma (a)$ with mean $a/b$, variance $a/{b^{2}}$, the cumulant function ${\kappa _{Y}}(z)=-a\log (1-iz/b)$, the Lévy density $m(x)=a{x^{-1}}{e^{-bx}}$, and the Bondesson–Rosinski representation (7) holds with
The Dickman OU process is another process with simple simulation algorithm, see [32]. We note that the Dickman distribution $D(a)$ is the distribution of a random variable V that satisfies the distributional fixed-point equation $V\stackrel{d}{=}{U^{1/a}}(1+V)$, where $a\gt 0$, “$\hspace{0.1667em}\stackrel{d}{=}\hspace{0.1667em}$” denotes the equality in distribution, U is independent of V and has the uniform distribution on $(0,1]$. The Dickman distribution $D(a)$ has the density
\[ p(x)=\left\{\begin{array}{l@{\hskip10.0pt}l}\frac{{e^{-\gamma a}}}{\Gamma (a)}{x^{a-1}},\hspace{1em}& 0\lt x\le 1,\\ {} \frac{{e^{-\gamma a}}}{\Gamma (a)}{x^{a-1}}-a{x^{a-1}}{\textstyle\textstyle\int _{0}^{x-1}}\frac{p(u)}{{(1+u)^{a}}}du,\hspace{1em}& x\gt 1,\end{array}\right.\]
where $\gamma =-{\Gamma ^{\prime }}(1)\approx 0.5772$ is Euler’s constant, with mean a, variance $a/2$, the cumulant function ${\kappa _{Y}}(z)=a{\textstyle\int _{0}^{1}}{u^{-1}}({e^{izu}}-1)du$, the Lévy density $m(x)=a{x^{-1}}{\mathbf{1}_{(0,1]}}(x)$, and the Bondesson–Rosinski representation (7) holds with
The above two processes have marginals with short tails and their Lévy measures are such that $xm(x)$ is finite. The case of OU processes with fat-tail marginals is more difficult because $xm(x)$ is not finite at zero. This means that sample paths of $Z(t)$ have infinitely many jumps and exact simulation is impossible, see [57].
2 Main results
We will concentrate on simulation of supOU processes with fat-tail marginals such as the inverse gamma distribution and the inverse Gaussian distribution. We note that supOU processes with fat-tail marginals usually have the Lévy density $m(x)$ such that $x\hspace{0.1667em}m(x)$ is singular at zero that disables the direct use of (6).
We propose to adapt the Bondesson–Rosinski representation (7) for simulation of the supOU processes which may have infinitely many jumps. In numerical studies, we usually simulate a stochastic process at equidistant points ${t_{j}}=j\Delta $ with time step Δ and, therefore, the contribution of frequent small jumps between points ${t_{1}},{t_{2}},\dots $ can be disregarded that can be viewed as the truncation of the Lévy measure.
Indeed, consider the Lévy measure ${\mu _{Y}}$ on the interval $[0,\infty )$. This measure on the interval $[0,\epsilon ]$ corresponds to small jumps of the supOU process $Y(t)$. We approximate $Y(t)$ by a stochastic process ${Y_{\epsilon }}(t)$ with the Lévy measure ${\mu _{Y}}$ truncated to the interval $[\epsilon ,\infty )$. Specifically, we construct ${Y_{\epsilon }}(t)$ as follows.
For the given Lévy density $m(x)$ of the process $Y(t)$, from the relation (6), we obtain
which gives negative values near zero if $x\hspace{0.1667em}m(x)$ is singular at zero, and $F(x)$ is an increasing function with $F(x)\to 1$ as $x\to \infty $. Solving $F(\epsilon )=0$ with respect to θ for given $\epsilon \in (0,\infty )$, we get the solution ${\theta _{\epsilon }}=\epsilon \hspace{0.1667em}m(\epsilon )$. By truncating $F(x)$ to the interval $[\epsilon ,\infty )$, we define the function
which is a cumulative distribution function. Now we will state a main theorem which will be used for simulation.
Theorem 5.
Let $Y(t)$ be a supOU process with positive marginals such that $xm(x)$ is singular at zero. Define the process ${Y_{\epsilon }}(t)$ via the Bondesson–Rosinski representation (7), where ${H_{k}}$ has the distribution function ${F_{\epsilon }}(x)$ defined in (14). Then ${\kappa _{Y}}(z)-{\kappa _{{Y_{\epsilon }}}}(z)\to 0$ as $\epsilon \to 0$ and the correlation functions of $Y(t)$ and ${Y_{\epsilon }}(t)$ coincide.
Proof.
We note that ${\theta _{\epsilon }}$ is the intensity of the Poisson process with jumps at time moments ${S_{k}}$. By construction, the process ${Y_{\epsilon }}(t)$ is well defined. Moreover, we have that ${F_{\epsilon }}(x)\to F(x)$ and ${\kappa _{Y}}(z)-{\kappa _{{Y_{\epsilon }}}}(z)={\textstyle\int _{0}^{\varepsilon }}({e^{izx}}-1-izx{1_{[-1,1]}}(x))m(x)dx\to 0$ as $\epsilon \to 0$. The correlation structure is specified by the measure $\pi (d\xi )$, which is the same for both supOU processes $Y(t)$ and ${Y_{\epsilon }}(t)$. □
Since the probability density function (pdf) is the inverse Fourier transform of the exponent of the cumulant function, the pdf of the marginal distribution of ${Y_{\epsilon }}(t)$ converges to that of $Y(t)$. Moreover, the finite-dimensional distributions of the process ${Y_{\epsilon }}(t)$ converge to those of $Y(t)$ and, consequently, the statistical properties of sample paths of ${Y_{\epsilon }}(t)$ approximate those of $Y(t)$ for small ϵ. Since the supOU process ${Y_{\epsilon }}(t)$ can be viewed as a surrogate model of $Y(t)$, the sampling error for modelling the stochastic process $Y(t)$ on a finite interval is larger than the error due to the replacement of the model $Y(t)$ by ${Y_{\epsilon }}(t)$ with small ϵ, see Supplementary Materials in [45] for details. From a practical point of view, we can take $\epsilon =0.001$. For other Lévy processes, the convergence under truncation of the Lévy measure is investigated in [29].
In summary, the simulation algorithm for a supOU process with the given measure π and the Lévy density $m(x)$ of the given marginal distribution is as follows.
In the above algorithm, we replace the infinite sum in the Bondesson–Rosinski representation (7) by the sum in which k is such that ${S_{k}}\gt -{T_{\mathrm{WarmUp}}}$, where ${T_{\mathrm{WarmUp}}}$ is a constant such that the correlation function $r(t)$ at $t={T_{\mathrm{WarmUp}}}$ is close to zero, say, $r({T_{\mathrm{WarmUp}}})\leqq 0.1$. We also note that the above algorithm produces processes with positive marginals and the computational time is proportional to $(T+{T_{\mathrm{WarmUp}}}){\theta _{\epsilon }}/(\eta \Delta )$ and takes seconds in our further examples.
If the marginal distribution is defined on $\mathbb{R}$, then the Lévy density $m(x)$ is also defined on $\mathbb{R}$ and we construct the simulation algorithm for a process with nonpositive marginals as the difference of two positive processes,
where M is the mean-correction parameter, ${Y_{\epsilon ,1}}(t)$ and ${Y_{\epsilon ,2}}(t)$ are independent stochastic processes, ${Y_{\epsilon ,1}}(t)$ has the representation (7) with positive ${H_{k}}$ from the distribution (14) and ${Y_{\epsilon ,2}}(t)$ has the representation (7) with positive ${H_{k}}$ from the distribution (14) with $m(x)$ replaced by $m(-x)$. In this construction, the process ${Y_{\epsilon ,2}}(t)$ corresponds to the Lévy density $m(x)$ for $x\lt 0$.
In the following sections we will confirm the theoretical behavior of the proposed algorithm by simulation for a wide range of supOU processes.
2.1 SupOU processes with positive marginals
2.1.1 SupOU process with inverse gamma marginals
The density of the inverse gamma distribution, R$\Gamma (a,b)$ for short, is
${J_{|a|}}(\cdot )$ and ${Y_{|a|}}(\cdot )$ are Bessel functions of the first and second kind, respectively.
\[ p(x)=\frac{{b^{a}}}{\Gamma (a)}\frac{1}{{x^{a+1}}}{e^{-b/x}},\hspace{1em}x\gt 0,\hspace{3.33333pt}a\gt 0,\hspace{3.33333pt}b\gt 0,\]
which has mean $\frac{b}{a-1}$ if $a\gt 1$, variance $\frac{{b^{2}}}{{(a-1)^{2}}(a-2)}$ if $a\gt 2$, and the cumulant function
where ${K_{a}}(\cdot )$ is the modified Bessel function of the second kind. Following [8], the Lévy density of the supOU process with inverse gamma marginals is
where
(15)
\[ {g_{a}}(u)=\frac{2}{t{\pi ^{2}}\big({J_{|a|}^{2}}(\sqrt{u})+{Y_{|a|}^{2}}(\sqrt{u})\big)},\]Fig. 2.
Top: Realizations of the supOU process ${Y_{\epsilon }}(t)$ with R$\Gamma (3,2)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of R$\Gamma (3,2)$ and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations, the x-axis represents the lag between points
We note that the intensity ${\theta _{\epsilon }}$ depends on a and b. For example, for $a=4$ and $\epsilon =0.001$, we have ${\theta _{\epsilon }}=23.5$ for $b=2$ and ${\theta _{\epsilon }}=29.2$ for $b=3$. Thus, we can evaluate the cumulative distribution function ${F_{\epsilon }}(x)$ by (14) and simulate the process ${Y_{\epsilon }}(t)$ using the Bondesson–Rosinski representation (7) with intensity ${\theta _{\epsilon }}$, ${R_{k}}\sim \pi $ and ${H_{k}}\sim {F_{\epsilon }}$. The R code of simulation of ${Y_{\epsilon }}(t)$ is given in [45].
In Figure 2 we depict several realizations of the supOU process ${Y_{\epsilon }}(t)$ with inverse gamma marginals and their characteristics. We need the large interval [0,3000] to show the behavior of the supOU process which shows distinct patterns on short intervals. The measure $\pi =\Gamma (1+\alpha ,\alpha )$ implies the long-range dependence for $\alpha \in (0,1]$ and short-range dependence for $\alpha \gt 1$. We can see that the empirical density is close to the true density of the inverse gamma distribution and the empirical autocorrelation function (acf) is close the true acf. We note that the empirical estimators of parameters of realizations of the supOU process are usually very far from the true values if the supOU process is observed over short intervals.
2.1.2 SupOU process with inverse Gaussian marginals
The inverse Gaussian distribution, IG$(a,b)$ for short, has the density
\[ p(x)=\frac{a{e^{ab}}}{\sqrt{2\pi }}{x^{-3/2}}{e^{-({a^{2}}/x+{b^{2}}x)/2}},\hspace{1em}x\gt 0,\hspace{3.33333pt}a\gt 0,\hspace{3.33333pt}b\gt 0,\]
with mean $\frac{a}{b}$, variance $\frac{a}{{b^{3}}}$ and the cumulant function ${\kappa _{Y}}(z)=ab-a\sqrt{{b^{2}}-2iz}$. Following [8], the Lévy density of the supOU process with inverse Gaussian marginals is
Fig. 3.
Top: Realizations of the supOU process with IG$(2,2)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of IG$(2,2)$ and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 3 we depict several realizations of the supOU process with inverse Gaussian marginals and their characteristics.
2.1.3 SupOU process with generalized inverse Gaussian marginals
The properties of supOU processes with GIG marginals were studied in [8]. In [67], simulation algorithms of OU processes with GIG marginals were constructed using the autoregressive formulae of order 1. The GIG distribution has the density
The GIG distribution is essential for constructing the family of Generalized Hyperbolic Lévy processes [11], which are highly effective tools for modeling phenomena in finance and turbulence [8, 15].
\[ p(x)=\frac{{(c/b)^{a/2}}}{2{K_{a}}(\sqrt{cb})}{x^{a-1}}{e^{-(cx+b/x)/2}},\hspace{1em}x\gt 0,\hspace{3.33333pt}c\gt 0,\hspace{3.33333pt}b\gt 0,\hspace{3.33333pt}a\in \mathbb{R},\]
with mean $\sqrt{b/c}\frac{{K_{a+1}}(\sqrt{cb})}{{K_{a}}(\sqrt{cb})}$, variance $\frac{b}{c}\left[\frac{{K_{a+2}}(\sqrt{cb})}{{K_{a}}(\sqrt{cb})}-{\left(\frac{{K_{a+1}}(\sqrt{cb})}{{K_{a}}(\sqrt{cb})}\right)^{2}}\right]$ and the cumulant function
\[ {\kappa _{Y}}(z)=\log \left({\left(\frac{c}{c-2iz}\right)^{a/2}}\frac{{K_{a}}(\sqrt{b(c-2iz)})}{{K_{a}}(\sqrt{cb})}\right).\]
The following distributions are special cases of GIG$(a,b,c)$:
-
• GIG$(a,0,c)$ is the gamma distribution $\Gamma (a,c)$, $a\gt 0$,
-
• GIG$(-a,b,0)$ is the inverse gamma distribution R$\Gamma (a,b)$, $a\gt 0$,
-
• GIG$(-1/2,b,c)$ is the inverse Gaussian distribution IG$(\sqrt{b},\sqrt{c})$,
-
• GIG$(1/2,b,c)$ is the reciprocal inverse Gaussian distribution RIG$(\sqrt{b},\sqrt{c})$,
-
• GIG$(1,b,c)$ is the positive hyperbolic distribution PH$(b,c)$,
-
• GIG$(-1,b,c)$ is the reciprocal positive hyperbolic distribution RPH$(b,c)$.
Following [8], the Lévy density of the supOU process with GIG marginals is
\[ m(x)=\frac{1}{x}\left(\frac{1}{2}{\int _{0}^{\infty }}{e^{-\frac{xu}{2b}}}{g_{a}}(u)du+\max \{0,a\}\right){e^{-\frac{cx}{2}}},\hspace{1em}x\gt 0,\]
where ${g_{a}}(u)$ is defined in (15).Fig. 4.
Top: Realizations of the supOU process with GIG$(1.5,2,1)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of GIG$(1.5,2,1)$ and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 4 we depict several realizations of the supOU process with generalized inverse Gaussian marginals and their characteristics.
2.1.4 SupOU process with Bessel marginals
The Bessel distribution B(a) (which can be viewed as the generalized Mc Kay distribution) has the density
with infinite mean, the Laplace transform $\mathbb{E}{e^{-sV}}={(1+s-\sqrt{{s^{2}}+2s})^{a}}$ and the cumulant function ${\kappa _{Y}}(z)=a\log (1-iz-\sqrt{-2iz-{z^{2}}})$, where ${I_{a}}(\cdot )$ is the modified Bessel functions of the first kind. Following [12], the Lévy density of the supOU process with Bessel marginals is given by
which implies a small finite value of ${\theta _{\epsilon }}$ even if $\epsilon \to 0$ because $x\hspace{0.1667em}m(x)$ is not singular at zero.
2.1.5 SupOU process with Mittag-Leffler marginals
The Mittag-Leffler distribution ML$(a)$ has the cumulative distribution function
with infinite mean, the density $p(x)={x^{a-1}}{E_{a,a}}(-{x^{a}})$, the Laplace transform $\mathbb{E}{e^{-sV}}=1/(1+{s^{a}})$ and the cumulant function
\[ {\kappa _{Y}}(z)=-\log \left(1+{(-iz)^{a}}\right)=-\log \left(1+|z{|^{a}}{e^{-i\operatorname{sign}(z)\pi a/2}}\right),\]
where ${E_{a}}(z)={E_{a,1}}(z)$ and
\[ {E_{a,b}}(z)={\sum \limits_{k=0}^{\infty }}\frac{{z^{k}}}{\Gamma (b+ak)},\hspace{1em}z\in \mathbb{C},\]
is the two-parameter Mittag-Leffler function. Following [12, 44], the Lévy density of the supOU process with Mittag-Leffler marginals is given by
which implies a small finite value of ${\theta _{\epsilon }}$ even if $\epsilon \to 0$ because $x\hspace{0.1667em}m(x)$ is not singular at zero.Fig. 6.
Top: Realizations of the supOU process with ML$(0.9)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom: The true density of ML$(0.9)$ and empirical densities of realizations
In Figure 6 we depict several realizations of the supOU process with Mittag-Leffler marginals and their characteristics.
2.1.6 SupOU process with positive α-stable marginals
The positive α-stable distribution has the Laplace transform $\mathbb{E}{e^{-sV}}={e^{-{s^{a}}}}$, $s\ge 0$, the density
implying infinite mean, and the cumulant function
Following [6], the Lévy density of the supOU process with α-stable marginals is given by
Therefore, we obtain ${\theta _{\epsilon }}=-1/({\epsilon ^{a}}\Gamma (-a))$ and
that enables us to exactly simulate the random variables with the cumulative distribution function ${F_{\epsilon }}(x)$.
(16)
\[ {p_{a}}(x)=\frac{1}{\pi }{\sum \limits_{k=1}^{\infty }}{(-1)^{k+1}}\frac{\Gamma (ka+1)}{k!{x^{ak+1}}}\sin (ka\pi ),\hspace{1em}x\gt 0,\hspace{3.33333pt}a\in (0,1),\]Fig. 7.
Top: Realizations of the supOU process with positive 0.9-stable marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom: The true density of the positive 0.9-stable distribution and empirical densities of realizations
In Figure 7 we depict several realizations of the supOU process with positive 0.9-stable marginals and their characteristics.
2.1.7 SupOU process with tempered stable marginals
The tempered stable distribution TS$(a,b)$ has the cumulant function
the mean ${b^{a-1}}a$, the variance ${b^{a-2}}a(1-a)$ and the density $p(x)={p_{a}}(x){e^{{b^{a}}-bx}}$, where ${p_{a}}(x)$ is defined in (16), see [13]. Following [9], the Lévy density of the supOU process with tempered stable marginals is given by
Fig. 8.
Top: Realizations of the supOU process with TS$(0.4,1)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of TS$(0.4,1)$ and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 8 we depict several realizations of the supOU process with tempered stable marginals and their characteristics.
2.1.8 SupOU process with log-normal marginals
The log-normal distribution LN$(a,b)$ has the density
\[ p(x)=\hspace{3.33333pt}\frac{1}{xb\sqrt{2\pi }}\exp \left(-\frac{{\left(\log x-a\right)^{2}}}{2{b^{2}}}\right),\hspace{1em}x\gt 0,\hspace{3.33333pt}a\in \mathbb{R},\hspace{3.33333pt}b\gt 0.\]
with mean $\exp (a+\frac{{b^{2}}}{2})$ and variance $(\exp ({b^{2}})-1)\exp (2a+{b^{2}})$. The log-normal distribution is not determined by its moments. The characteristic function of LN$(a,b)$ is
\[ \phi (z)=\frac{{e^{\frac{{\pi ^{2}}}{8{b^{2}}}}}}{\sqrt{\pi }}{\int _{-\infty }^{\infty }}{e^{-z{e^{a+\sqrt{2}bx}}-\frac{i}{\sqrt{2}b}\pi x-{x^{2}}}}dx,\]
see [33]. Following [2], the characteristic function has a good approximation
\[ \phi (z)\approxeq \frac{\exp \left(-\frac{{W^{2}}(-iz{b^{2}}{e^{a}})+2W(-iz{b^{2}}{e^{a}})}{2{b^{2}}}\right)}{\sqrt{1+W(-iz{b^{2}}{e^{a}})}},\]
where $W(\cdot )$ is the Lambert W function which satisfies $W(x){e^{W(x)}}=x$, see [22].Following [18], there is no explicit formula for the Lévy density of the log-normal distribution. Thus, we directly compute the Lévy density $m(x)$ by the method based on the Bromwich integral, which is equivalent to formula (9) with
\[ {(\log \phi (z))^{\prime\prime }}\approxeq \frac{\left({W^{3}}(-iz{b^{2}}{e^{a}})+3{W^{2}}(-iz{b^{2}}{e^{a}})+(3+{b^{2}}/2)W(-iz{b^{2}}{e^{a}})+1+3{b^{2}}/2\right){W^{2}}(-iz{b^{2}}{e^{a}})}{{(1+W(-iz{b^{2}}{e^{a}}))^{4}}{z^{2}}{b^{4}}},\]
which is absolutely integrable because $W(x)\sim \log (x)$ as $x\to \infty $.Fig. 9.
Top: Realizations of the supOU process with LN$(0,1)$ marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of LN$(0,1)$ and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
2.2 SupOU processes with marginals on entire real line
2.2.1 SupOU process with hyperbolic cosine marginals
The hyperbolic cosine distribution has the density $p(x)=1/(\pi \cosh (x))$, $x\in \mathbb{R}$, with zero mean, variance ${\pi ^{2}}/4$ and the characteristic function $\phi (z)=1/\cosh (\pi z/2)$. Following [19], the Lévy density of the supOU process with hyperbolic cosine marginals is given by
which implies a small finite value of ${\theta _{\epsilon }}$ even if $\epsilon \to 0$ because $x\hspace{0.1667em}m(x)$ is not singular at zero.
Fig. 10.
Top: Realizations of the supOU process with hyperbolic cosine marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the hyperbolic cosine distribution and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 10 we depict several realizations of the supOU process with hyperbolic cosine marginals and their characteristics.
2.2.2 SupOU process with normal inverse Gaussian marginals
The supOU processes with NIG marginals were extensively studied in [4]. The NIG distribution has the density
\[ p(x)=\frac{a\delta {K_{1}}\left(a\sqrt{{\delta ^{2}}+{(x-c)^{2}}}\right)}{\pi \sqrt{{\delta ^{2}}+{(x-c)^{2}}}}\hspace{0.2778em}{e^{\delta \gamma +\beta (x-c)}},\hspace{1em}x\in \mathbb{R},\hspace{3.33333pt}a\gt 0,\hspace{3.33333pt}\gamma =\sqrt{{a^{2}}-{\beta ^{2}}},\]
where c is the location parameter, a is the tail parameter, β is the asymmetry parameter and δ is the scale parameter. The NIG distribution has semi-heavy tails, specifically, the tails of $p(x)$ behave like $|x{|^{-3/2}}{e^{-a|x|+\beta x}}$ as $x\to \pm \infty $. The NIG distribution has the mean $c+\delta \beta /\gamma $, the variance $\delta {a^{2}}/{\gamma ^{3}}$ and the cumulant function
Following [4, 6], the Lévy density of the supOU process with NIG marginals is given by
Consider the case of $c=0$, $\beta =0$ and $\delta =1$, which will be denoted as NIG$(a)$.Fig. 11.
Top: Realizations of the supOU process with NIG(0.6) marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the NIG distribution and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 11 we depict several realizations of the supOU process with NIG(0.6) marginals and their characteristics.
2.2.3 SupOU process with Student marginals
The Student distribution ST$(c,a,\nu )$ has the density
\[ p(x)=\frac{\Gamma ((\nu +1)/2)}{a\sqrt{\pi \nu }\Gamma (\nu /2)}{\left(1+\frac{1}{\nu }{\left(\frac{x-c}{a}\right)^{2}}\right)^{-(\nu +1)/2}},\hspace{1em}x\in \mathbb{R},\hspace{3.33333pt}\nu \gt 0,\hspace{3.33333pt}a\gt 0,\]
which has mean c if $\nu \gt 1$, variance $\frac{\nu a}{\nu -2}$ if $\nu \gt 2$ and the cumulant function
\[ {\kappa _{Y}}(z)=\log \left(\frac{{\left(\sqrt{\nu }a|z|\right)^{\nu /2}}{K_{\nu /2}}\left(\sqrt{\nu }a|z|\right)}{\Gamma (\nu /2)\hspace{3.33333pt}{2^{\nu /2-1}}\hspace{3.33333pt}}\right)+icz.\]
Following [35], the Lévy density of the supOU process with Student marginals is given by
\[ m(x)=\frac{a}{2|x|}{\int _{0}^{\infty }}{e^{-|x|\sqrt{u/\nu }}}{g_{\nu /2}}(u)du,\hspace{1em}x\in \mathbb{R},\]
where ${g_{a}}(u)$ is defined in (15), and has the property ${x^{2}}m(x)/a=\sqrt{\nu }/\pi +(1-\nu )|x|/4+o(x)$ as $x\to 0$.Fig. 12.
Top: Realizations of the supOU process with ST(0,1,2.5) marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the Student distribution and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 12 we depict several realizations of the supOU process with ST(0,1,2.5) marginals and their characteristics.
2.2.4 SupOU process with Cauchy marginals
The Cauchy distribution has the density
\[ p(x)=\frac{1}{\pi a}{\left(1+{\left(\frac{x-c}{a}\right)^{2}}\right)^{-1}},\hspace{1em}x\in \mathbb{R},\hspace{3.33333pt}a\gt 0,\]
which has infinite absolute moment of order 1 and the cumulant function ${\kappa _{Y}}(z)=a|z|+icz$. Following [19], the Lévy density of the supOU process with Cauchy marginals is given by
Fig. 13.
Top: Realizations of the supOU process with Cauchy marginals with $a=1,\hspace{3.33333pt}c=0$ on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom: The true density of the Cauchy distribution and empirical densities of realizations
In Figure 13 we depict several realizations of the supOU process with Cauchy marginals and their characteristics.
2.2.5 SupOU process with generalized Linnik marginals
The generalized Linnik distribution GL$(a,b)$ has the characteristic function
\[ \phi (z)=\frac{1}{{(1+|z{|^{a}})^{b}}},\hspace{1em}z\in \mathbb{R},\hspace{3.33333pt}a\in (0,2],\hspace{3.33333pt}b\gt 0.\]
The density does not have an explicit form and can be computed as
\[ p(x)=\frac{1}{2\pi }{\int _{-\infty }^{\infty }}\phi (z)\cos (xz)dz=\frac{1}{\pi }\mathrm{Im}{\int _{0}^{\infty }}\frac{{e^{-y|x|}}}{{(1+{e^{-i\pi a/2}}{y^{a}})^{b}}}dy,\hspace{1em}x\ne 0\hspace{3.33333pt}\mathrm{for}\hspace{3.33333pt}a\le 1,\]
which has zero mean if $a\gt 1$ and has the property
\[ p(x)=\frac{b}{2\pi }\sin (a\pi /2)\Gamma (1+a)|x{|^{-1-a}}+o(|x{|^{-1-a}})\hspace{1em}\mathrm{as}\hspace{3.33333pt}|x|\to \infty ,\hspace{3.33333pt}a\in (0,2),\]
see [3, 25] for details. The GL$(2,b)$ distribution is a special case, for example, for $a=2$ and $b=1$ we have the density $p(x)=\frac{1}{2}{e^{-|x|}}$ of the Laplace distribution with zero mean, variance 2 and the Lévy density $m(x)=|x{|^{-1}}{e^{-|x|}}$.The GL$(a,1)$ distribution is the Linnik distribution with the Lévy density
where ${p_{a}}(x)$ is defined in (16).
The Lévy density of the generalized Linnik distribution is given by
where
\[ f(x,u)=\left\{\begin{array}{l@{\hskip10.0pt}l}\frac{1}{\pi }{\textstyle\textstyle\sum _{k=1}^{\infty }}{(-1)^{k+1}}\frac{\Gamma (ak+1)}{k!}\frac{{u^{k}}}{{x^{ak+1}}}\sin (\pi ak/2),\hspace{1em}& a\in (0,1),\\ {} \frac{1}{\pi }{\textstyle\textstyle\sum _{k=1}^{\infty }}{(-1)^{k+1}}\frac{\Gamma (k/a+1)}{k!}\frac{{x^{k-1}}}{{u^{k/a}}}\sin (\pi k/2),\hspace{1em}& a\in (1,2),\end{array}\right.\]
see [43]. Alternatively, we can directly compute the Lévy density $m(x)$ by formula (9) with
which is absolutely integrable.Fig. 14.
Top: Realizations of the supOU process with generalized Linnik marginals with $a=1.5,\hspace{3.33333pt}b=1$ on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the generalized Linnik distribution and empirical densities of realizations. Bottom-Right: The specified acf (dotted) and empirical acf (solid line) of realizations
In Figure 14 we depict several realizations of the supOU process with generalized Linnik marginals and their characteristics.
2.2.6 SupOU process with generalized Gaussian marginals
The random variable V has the (one-parameter) generalized Gaussian distribution if its density is
\[ p(x)={c_{a}}\exp (-|x{|^{a}}/2),\hspace{1em}{c_{a}}={2^{-\frac{a+1}{a}}}a/\Gamma (1/a),\hspace{3.33333pt}x\in \mathbb{R},\hspace{3.33333pt}a\gt 0,\]
with zero mean and variance ${2^{2/a}}\Gamma (\frac{3}{a})/\Gamma (\frac{1}{a})$. Following [24], the characteristic function of V is given by
\[ \phi (z)=2{c_{a}}{\int _{0}^{\infty }}\cos (xz)\exp (-{x^{a}}/2)dx=\pi {c_{a}}\frac{a|z{|^{\frac{a}{a-1}}}}{|a-1|}{\int _{0}^{1}}{U_{a}}(x){e^{-|z{|^{\frac{a}{a-1}}}{U_{a}}(x)}}dx,\]
where
\[ z\in \mathbb{R},\hspace{1em}{U_{a}}(x)={\left(\frac{\sin (\pi xa/2)}{\cos (\pi x/2)}\right)^{\frac{a}{a-1}}}\frac{\cos (\pi x(a-1)/2)}{\cos (\pi x/2)}.\]
The random variable V is self-decomposable for $a\in (0,1)\cup \{2\}$.There is no explicit formula for the Lévy density of the generalized Gaussian distribution. Thus, we directly compute the Lévy density $m(x)$ by formula (9) since ${(\log \phi (z))^{\prime\prime }}$ is absolutely integrable, see [24, App. J].
Fig. 15.
Top: Realizations of the supOU process with GG(0.9) marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the generalized Gaussian and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 15 we depict several realizations of the supOU process with GG(0.9) marginals and their characteristics.
2.2.7 SupOU process with Gumbel marginals
The Gumbel distribution G$(c,a)$, which is also known as the extreme value distribution, has the cumulative distribution function
\[ {F_{\mathrm{G}}}(x)=\exp \left(-\exp \left(-\frac{x-c}{a}\right)\right),\hspace{1em}x\in \mathbb{R},\hspace{3.33333pt}c\in \mathbb{R},\hspace{3.33333pt}a\gt 0,\]
with mean $c+a\gamma $, where $\gamma \approx 0.5772$, variance $\frac{{\pi ^{2}}}{6}{a^{2}}$ and the characteristic function
Following [1], the Lévy density of the supOU process with Gumbel marginals is given by
Since the Lévy density produces the mean-shifted marginal distribution, we have to set the mean-correcting parameter $M=M(a,c)$.Fig. 16.
Top: Realizations of the supOU process with G(0,1) marginals on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the Gumbel and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 16 we depict several realizations of the supOU process with G(0,1) marginals and their characteristics.
2.2.8 SupOU process with Rosenblatt marginals
The Rosenblatt distribution was introduced in [55] and has zero mean, unit variance and the characteristic function
where ${S_{0}}$ is a small neighborhood of zero, ${\sigma _{a}}=\sqrt{(1-2a)(1-a)/2}$ and
(17)
\[ \phi (z)=\exp \left(\frac{1}{2}{\sum \limits_{k=2}^{\infty }}{(2iz{\sigma _{a}})^{k}}\frac{{c_{a,k}}}{k}\right),\hspace{1em}a\in [0,1/2],\hspace{3.33333pt}z\in {S_{0}},\]
\[ {c_{a,k}}={\int _{0}^{1}}\cdots {\int _{0}^{1}}|{x_{1}}-{x_{2}}{|^{-a}}|{x_{2}}-{x_{3}}{|^{-a}}\cdots |{x_{k-1}}-{x_{k}}{|^{-a}}|{x_{k}}-{x_{1}}{|^{-a}}d{x_{1}}\cdots d{x_{k}},\]
an explicit formula for the density is not available. The random variable V with the Rosenblatt distribution can be given as
where ${\varepsilon _{n}}$ are i.i.d. random variables with the standard normal distribution and ${\lambda _{a,1}},{\lambda _{a,2}},\dots $ are such that ${\textstyle\sum _{n=1}^{\infty }}{\lambda _{a,n}^{k}}={\sigma _{a}^{k}}{c_{a,k}}$ for all $k=2,3,\dots $. In particular, we have ${\textstyle\sum _{n=1}^{\infty }}{\lambda _{a,n}^{2}}=1/2$, ${\textstyle\sum _{n=1}^{\infty }}{\lambda _{a,n}}=\infty $ and the Laplace transform of V is given by
\[ {\phi _{LT}}(s)=\mathbb{E}({e^{-sV}})=\exp \left(-{\sum \limits_{n=1}^{\infty }}\left(\frac{1}{2}\log (1+2{\lambda _{a,n}}s)-{\lambda _{a,n}}s\right)\right),\hspace{1em}s\gt -\frac{1}{2{\lambda _{a,1}}}.\]
Following [60], the Lévy density of the supOU process with Rosenblatt marginals is given by
\[ m(x)=\frac{1}{2x}{\sum \limits_{n=1}^{\infty }}\exp \left(-\frac{x}{2{\lambda _{a,n}}}\right),\hspace{1em}x\gt 0,\]
and ${\lambda _{a,n}}$ can also be computed as eigenvalues of the integral operator ${\tilde{K}_{a}}:{L^{2}}(0,1)\to {L^{2}}(0,1)$ defined as
Following [60], the eigenvalues ${\lambda _{a,n}}$ admit the accurate approximation
\[ {\lambda _{a,n}}\approx {C_{a}}{n^{a-1}}\hspace{1em}\mathrm{for}\hspace{3.33333pt}n\gt 30,\hspace{1em}{C_{a}}=\frac{2{\sigma _{a}}}{{\pi ^{1-a}}}\Gamma (1-a)\sin (\pi a/2),\]
and should be computed numerically for small n. An accurate approximation of the eigenvalues ${\lambda _{a,n}}$ for all n is proposed in [46]. Since the Lévy density of the Rosenblatt distribution is defined for positive x like in cases of supOU processes with positive marginals, we have to set the mean-correcting parameter $M=M(a)$.Fig. 17.
Top: Realizations of the supOU process with Rosenblatt marginals with $a=0.3$ on the interval $[0,3000]$ with time step $\Delta =0.5$ for the measure $\pi =\Gamma (1+\alpha ,\alpha )$ and $\pi ={\delta _{\lambda }}$. Bottom-Left: The true density of the Rosenblatt and empirical densities of realizations. Bottom-Right: The true acf (dotted) and empirical acf (solid line) of realizations
In Figure 17 we depict several realizations of the supOU process with Rosenblatt marginals and their characteristics.
3 Conclusion
We proposed the universal simulation algorithm for supOU processes with specified marginal distributions and correlation functions that provides the complete viewpoint on the structure of these processes. Our algorithm is based on the truncation of the Lévy density in cases where $x\hspace{0.1667em}m(x)$ has a singularity at zero. We have applied the simulation algorithm for supOU processes with 16 marginal distributions and established a repository [45] containing R scripts to facilitate the use of this algorithm. Our extensive numerical study confirms that the empirical density of realizations is close to the specified marginal density and the empirical acf is close to the true acf if the second moment of the marginal distribution is finite and to the specified acf otherwise. If a process is required with marginals not included in the list of 18 distributions (see Table 1), the target marginal distribution can be approximated by one of the 18 considered distributions. The supOU process can then be simulated using this fitted marginal distribution.
