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)$ .

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$ .