Hostname: page-component-586b7cd67f-l7hp2 Total loading time: 0 Render date: 2024-11-21T06:44:54.290Z Has data issue: false hasContentIssue false

A remark on exact simulation of tempered stable Ornstein–Uhlenbeck processes

Published online by Cambridge University Press:  02 May 2024

Takuji Arai*
Affiliation:
Keio University
Yuto Imai*
Affiliation:
Nishogakusha University
*
*Postal address: Department of Economics, Keio University, 2-15-45 Mita, Minato-ku, Tokyo, 108-8345, Japan. Email: arai@econ.keio.ac.jp
**Postal address: Faculty of International Politics and Economics, Nishogakusha University, 6-16 Sanbancho, Chiyoda-ku, Tokyo, 102-8336, Japan. Email: y-imai@nishogakusha-u.ac.jp
Rights & Permissions [Opens in a new window]

Abstract

Qu, Dassios, and Zhao (2021) suggested an exact simulation method for tempered stable Ornstein–Uhlenbeck processes, but their algorithms contain some errors. This short note aims to correct their algorithms and conduct some numerical experiments.

Type
Letter to the Editor
Copyright
© The Author(s), 2024. Published by Cambridge University Press on behalf of Applied Probability Trust

1. Introduction

A stochastic process $Z=\{Z_t\}_{t\geq0}$ is said to be a tempered stable (TS) subordinator if it is a driftless subordinator with the Lévy measure $\nu(\mathrm{d} y)=\theta y^{-\alpha-1}{\mathrm{e}}^{-\beta y}\,\mathrm{d} y$ , $y>0$ , where $\alpha\in(0,1)$ and $\beta,\theta\in{\mathbb R}^+$ . In this case, we call the distribution of $Z_1$ a tempered stable distribution with parameters $\alpha,\beta,\theta$ , and denote it by $\mathop{\mathrm{TS}}\nolimits(\alpha,\beta,\theta)$ . In addition, a process $X=\{X_t\}_{t\geq0}$ is said to be a TS-based Ornstein–Uhlenbeck (OU-TS) process if it is a solution to the following stochastic differential equation:

(1) \begin{align} \mathrm{d} X_t = -\delta X_t\,\mathrm{d} t + \rho\,\mathrm{d} Z_t, \qquad X_0>0,\end{align}

where $\delta>0$ and $\rho>0$ . For any $t\geq0$ and $\tau>0$ , we have

\[ X_{t+\tau} = {\mathrm{e}}^{-\delta \tau}X_t + \rho\int_t^{t+\tau}{\mathrm{e}}^{-\delta(t+\tau-s)}\,\mathrm{d} Z_s.\]

Qu et al. [Reference Qu, Dassios and Zhao2] suggested an exact simulation algorithm for $X_{t+\tau}$ given $X_t$ . In addition, they separately gave another algorithm available only for the case of $\alpha=\frac{1}{2}$ . Here we correct the two algorithms and introduce the results of some numerical experiments.

2. Mathematical background and algorithms

The infinitesimal operator ${\mathcal A}$ of X is given by

\begin{align*} {\mathcal A}f(x,t) = \frac{\partial f}{\partial t} - \delta x\frac{\partial f}{\partial x} + \int_0^\infty\{f(x+\rho y,t)-f(x,t)\}\,\nu(\mathrm{d} y),\end{align*}

where $f\colon{\mathbb R}^+\times{\mathbb R}^+\to{\mathbb R}$ is differentiable on x and t. We can derive this by applying [Reference Applebaum1, (6.36)] to the stochastic differential equation (1). [Reference Qu, Dassios and Zhao2, (3.2)] gave a representation of ${\mathcal A}$ incorrectly as follows:

\[ {\mathcal A}f(x,t) = \frac{\partial f}{\partial t} - \delta x\frac{\partial f}{\partial x} + \rho\int_0^\infty\{f(x+y,t)-f(x,t)\}\,\nu(\mathrm{d} y).\]

Thus, all the subsequent arguments in [Reference Qu, Dassios and Zhao2] must be corrected, but this error does not affect the case of $\rho=1$ . Now, we fix $t\geq0$ and $\tau>0$ , and define a process $Y=\{Y_s\}_{t\leq s\leq t+\tau}$ as

\[ Y_s \;:\!=\; \exp\bigg\{{-}X_s\kappa{\mathrm{e}}^{\delta s} + \int_0^s\Phi(\rho\kappa{\mathrm{e}}^{\delta u})\,\mathrm{d} u\bigg\},\]

where $\kappa\in{\mathbb R}$ , and $\Phi$ is the Laplace exponent of Z, i.e. $\Phi(u)\;:\!=\;\int_0^\infty(1-{\mathrm{e}}^{-uy})\,\nu(\mathrm{d} y)$ . When ${\mathcal A} f(x,t)=0$ , the process $f(X_t,t)$ is a martingale. From this point of view, we can see that Y is a martingale. For any $\eta\in{\mathbb R}^+$ , taking $\kappa=\eta{\mathrm{e}}^{-\delta(t+\tau)}$ , we obtain

(2) \begin{align} {\mathbb E}\big[{\mathrm{e}}^{-\eta X_{t+\tau}}\mid X_t\big] = \exp\bigg\{{-}\eta X_t{\mathrm{e}}^{-\delta\tau} - \int_{\rho\eta{\mathrm{e}}^{-\delta\tau}}^{\rho\eta}\frac{\Phi(z)}{\delta z}\,\mathrm{d} z\bigg\}.\end{align}

From the view of [Reference Qu, Dassios and Zhao2, (3.7)–(3.9)], (2) implies

(3) \begin{align} {\mathbb E}[{\mathrm{e}}^{-\eta X_{t+\tau}} \mid X_t] & = \exp\bigg\{{-}\eta wX_t - \frac{\rho^\alpha\theta(1-w^\alpha)}{\alpha\delta} \int_0^\infty(1-{\mathrm{e}}^{-\eta s})s^{-\alpha-1}{\mathrm{e}}^{-({\beta}/{w\rho})s}\,\mathrm{d} s\bigg\} \nonumber \\ & \quad \times \exp\bigg\{{-}\frac{\theta\beta^\alpha\Gamma(1-\alpha)D_w}{\alpha\delta} \int_0^\infty(1-{\mathrm{e}}^{-\eta s})\nonumber \\ & \quad \times\int_1^{1/w}\frac{(({\beta}/{\rho})u)^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1} {\mathrm{e}}^{-({\beta}/{\rho})us}f_V(u)\,\mathrm{d} u\,\mathrm{d} s\bigg\},\end{align}

where $w\;:\!=\;{\mathrm{e}}^{-\delta\tau}$ , $\Gamma(\cdot)$ is the Gamma function, and

\[ D_w \;:\!=\; \int_1^{1/w}(u^{\alpha-1}-u^{-1})\,\mathrm{d} u = \frac{w^{-\alpha}-1}{\alpha} + \ln w, \quad f_V(u) \;:\!=\; \frac{u^{\alpha-1}-u^{-1}}{D_w}, \qquad u\in[1,{1}/{w}].\]

Equation (3) can be obtained by replacing $\theta$ and $\beta$ in [Reference Qu, Dassios and Zhao2, Theorem 3.1] with $\rho^{\alpha-1}\theta$ and $\beta/\rho$ , respectively. Thus, [Reference Qu, Dassios and Zhao2, Algorithm 3.2] can be corrected by replacing all $\theta$ s and $\beta$ s appearing there in the same way. As for the correction of [Reference Qu, Dassios and Zhao2, Algorithm 3.4], we have only to change the distribution of $\widetilde{IG}$ into

$$\mathop{\mathrm{IG}}\nolimits\bigg(\frac{2\rho}{c\delta}(\sqrt{w}-w),\frac{4\rho}{\delta^2}(1-\sqrt{w})^2\bigg),$$

where $\mathop{\mathrm{IG}}\nolimits(\mu,\lambda)$ denotes the IG distribution with the mean parameter $\mu$ and the rate parameter $\lambda$ . Furthermore, [Reference Qu, Dassios and Zhao2, Proposition 6.1] can be corrected with the same replacements of $\theta$ and $\beta$ as above, but additionally, the function $h(\cdot)$ needs to be replaced by $h(\cdot/\rho)$ .

Remark 1. Denoting the solution to (1) by $X^{\rho,X_0}$ with emphasis on $\rho>0$ and the initial value $X_0>0$ , we have $X^{\rho,X_0}=\rho X^{1,X_0/\rho}$ for any $\rho>0$ and $X_0>0$ . Thus, the algorithm with $\rho=1$ can be generalized to any case of $\rho>0$ .

3. Numerical results

As can be seen in [Reference Qu, Dassios and Zhao2, Tables 1 and 2], even using the original algorithms in [Reference Qu, Dassios and Zhao2] the errors are kept small enough as long as the means are computed. Thus, we compute the second moments instead and compare the results of the original and corrected algorithms.

Here, we execute simulations for an OU-TS process with $\alpha=0.25$ . For the other parameters, we set $\delta=0.2$ , $\beta=0.5$ , $\theta=0.25$ and vary the value of $\rho$ as 0.5, 1, 2, 5. We set $X_0=10.0$ and simulate $X_{0.5}$ ; next, we simulate $X_1$ using the value of $X_{0.5}$ , which is repeated until we simulate $X_5$ . We carried out the simulation one million times, and compared their mean square with the second moment ${\mathbb E}[X_5^2]$ , where we can calculate ${\mathbb E}[X_5^2]$ by using [Reference Qu, Dassios and Zhao2, (3.5)] and (2) as follows:

\[ {\mathbb E}[X_5^2] = \bigg\{wX_0 + \frac{\rho\theta}{\delta\beta^{1-\alpha}}(1-w)\Gamma(1-\alpha)\bigg\}^2 + \frac{\rho^2\theta}{2\delta\beta^{2-\alpha}}(1-w^2)\Gamma(2-\alpha),\]

where $w={\mathrm{e}}^{-5\delta}$ . The algorithms were coded in MATLAB (R2022b).

The simulation results are given in Table 1. Note that “Estim” in the third column represents the mean square of one million simulation results, and “Diff” and “Error” in the fourth and fifth columns are defined as $\mathrm{Diff} \;:\!=\; \mathrm{Estim} - {\mathbb E}[X_5^2]$ and $\mathrm{Error} \;:\!=\; ({\mathrm{Diff}}/{{\mathbb E}[X_5^2]})\times100$ , respectively. The last three columns display the results for the original algorithm, [Reference Qu, Dassios and Zhao2, Algorithm 3.2]. As seen in Table 1, the errors of the corrected algorithm are small enough regardless of the value of $\rho$ , but for the original algorithm this is not the case. Furthermore, similar results were obtained for the cases of OU-TS with $\alpha=0.75$ and OU-IG (inverse Gaussian), but they are not tabulated.

Table 1. OU-TS process with $\alpha=0.25$ .

Funding information

Takuji Arai gratefully acknowledges the financial support of the MEXT Grants in Aid for Scientific Research (C) Nos. 18K03422 and 22K03419. Yuto Imai also gratefully acknowledges the financial support of the MEXT Grant in Aid for Young Research No. 21K13327.

Competing interests

There were no competing interests to declare which arose during the preparation or publication process of this article.

References

Applebaum, D. (2009). Lévy Processes and Stochastic Calculus. Cambridge University Press.CrossRefGoogle Scholar
Qu, Y., Dassios, A. and Zhao, H. (2021). Exact simulation of Ornstein–Uhlenbeck tempered stable processes. J. Appl. Prob. 58, 347371.CrossRefGoogle Scholar
Figure 0

Table 1. OU-TS process with $\alpha=0.25$.