Evaluation of integrals with fractional Brownian motion for different Hurst indices

Abstract

In this paper, we will evaluate integrals that define the conditional expectation, variance and characteristic function of stochastic processes with respect to fractional Brownian motion (fBm) for all relevant Hurst indices, i.e. $H\in (0,1)$. Particularly, the fractional Ornstein–Uhlenbeck (fOU) process gives rise to highly nontrivial integration formulas that need careful analysis when considering the whole range of Hurst indices. We will show that the classical technique of analytic continuation, from complex analysis, provides a way of extending the domain of validity of an integral from $H\in (1/2,1)$ to the larger domain $H\in (0,1)$. Numerical experiments for different Hurst indices confirm the robustness and efficiency of the integral formulations presented. Moreover, we provide accurate and highly efficient financial option pricing results for processes that are related to the fOU process, with the help of Fourier cosine expansions.

Keywords:

• Fractional Brownian motion
• fractional Ornstein–Uhlenbeck process
• conditional density function
• analytic continuation
• COS method option pricing

2020 AMS Subject Classifications:

• 60G22
• 65D30
• 91G60
• 91G20

1. Introduction

In this paper, we will evaluate the integrals that define the conditional expectation and variance of stochastic processes that are based on fractional Brownian motion (fBm). For this, we first define fBm, to which a Hurst index is associated, as a generalization of standard Brownian motion (sBm).

Let $(\mathrm{\Omega },\mathrm{\Sigma },\mathbb{P})$ be a complete probability space, ${\mathbf{B}}^H=\{B_t^H,t\ge 0\}$ be an fBm with Hurst index $H\in (0,1)$ defined on this probability space, and ${\mathcal{F}}_t$ be its natural filtration.

Fractional Brownian motion, ${\mathbf{B}}^H$, is then uniquely characterized by the following properties [22]:

1. ${\mathbf{B}}^H$ has stationary increments, which means $B_t^H-B_s^H\sim B_{t-s}^H$, for $0\le s\le t$;

2. $B_0^H=0$ and $\mathbb{E}\left[B_t^H\right]=0$, for $t\ge 0$;

3. $\mathbb{E}[B_t^H{]}^2=t^{2H}$, for $t\ge 0$;

4. $B_t^H$ has a Gaussian distribution, for t>0.

From the first three properties, it follows that the covariance function is given by (1) $\mathrm{C}\mathrm{o}\mathrm{v}(B_t^H,B_s^H)=\frac{1}{2}(|t{|}^{2H}+|s{|}^{2H}-|t-s{|}^{2H}).$(1) A special case of Equation (1(1) $\mathrm{C}\mathrm{o}\mathrm{v}(B_t^H,B_s^H)=\frac{1}{2}(|t{|}^{2H}+|s{|}^{2H}-|t-s{|}^{2H}).$(1) ) is sBm ($H=1/2$), with covariance function $\mathrm{C}\mathrm{o}\mathrm{v}(B_t^{\frac{1}{2}},B_s^{\frac{1}{2}})$ $=min(t,s)$. Particularly for fBm, and different from sBm, the increments are not independent and fBm is neither a Markov process nor a semi-martingale, for $H\ne 1/2$. As a consequence, fBm enables modelling of stochastic processes that are self-similar or exhibit long-range dependence on time. Properties of fBm have been presented, for example, in Refs [4] and [24].

Fractional Brownian motion forms the basis for a broad class of stochastic processes, like the fractional Ornstein–Uhlenbeck (fOU), the fractional Cox–Ingersoll–Ross (fCIR) and the geometric fOU (GfOU) processes. Such processes play an increasingly important role in diverse application fields, like hydrology [19,21], telecommunication [23,27], weather forecast [30,33] and epidemic disease (such as the Coronavirus) modelling [1]. FBm has also been used in finance for the modelling of commodities, weather (considered as an asset) and currencies in financial products, that are based on these assets, in Refs [2,3,6,29,36] among others, for bond markets in Ref. [14] and recently for rough volatility in Refs [15,18]. In these application areas, the computation of conditional expectations is required to compute the average behaviour under uncertainty and also to price options on these non-tradable assets in finance.

For Markov processes, the accurate computation of conditional expectations and variances has been given ample attention. In the diverse literature, we encounter, for such processes, approaches based on numerical integration, Fourier or Laplace inversion, Monte Carlo simulation and recently also neural network approximations. For non-Markov processes, to which fBm belongs, the literature is more scarce, but some significant books and papers have been presented in the past 20 years, like Refs [2,4,14,24]. The conditional characteristic function for the fOU process has been derived in Ref. [14] and was evaluated in Refs [14,34], where involved integration formulas had to be calculated numerically. These integration formulas are nontrivial when considering all Hurst indices $H\in (0,1)$, due to the occurrence of singularities in the integrand. Particularly, for the range $H\in (0,1/2)$, we encounter approximation formulas in the literature that are not easily understood from a numerical analysis point-of-view.

In this paper, we analyse the integration formulas for the expectation and variance for all relevant Hurst indices and aim to approximate them numerically on the basis of accurate (classical) approximations and numerical techniques. The technique of analytic continuation [31,32] is used in this paper to deal with the integration formulas for certain Hurst index values. Analytic continuation provides a way of extending the domain over which a complex function is defined. Here, it is used to extend the domain of validity of an integral which is defined on an interval with respect to the Hurst index, from $H\in (1/2,1)$, to the larger domain, $H\in (0,1)$. Here, we calculate the conditional expectation, variance and the conditional characteristic function of fBm and of fBm-driven stochastic differential equations (SDEs). The accuracy of our numerical results is verified by results obtained by Monte Carlo simulation.

As an application in finance, we discuss the option pricing problem where the asset price is modelled as a fractional stochastic process. Clearly, fBm is not a semi-martingale, but some assets, like those based on the weather or volatility, are non-tradable, so the well-known financial derivative pricing theory, based on the no-arbitrage and market completeness, cannot be used anyway. Valuation of such a derivative is then based on the expected discounted value approach under the real-world probability measure. With the characteristic function for the fOU process available, we can compute the conditional expectation by means of the COS method, a Fourier-based numerical integration technique, see Refs [9,10]. The COS method is based on an approximation of the appearing (conditional) probability density function by a Fourier cosine series expansion. The characteristic function of some of the fractional processes considered here is either unavailable or not easy to derive. However, by means of a transformation of variables, the COS method can still be employed.

This paper is organized as follows. Section 2 briefly introduces fBm and the related stochastic processes. In Section 3, we analyse the conditional variance of the fOU process and present a convenient numerical method to compute the appearing integrals. Option pricing with the COS method is described in Section 4. Section 5 presents numerical results, including conditional distributions and option pricing. Section 6 concludes the paper. Some examples and pointers for the analytic continuation technique are presented in Appendix 1.

2. Fractional Brownian motion and related processes

In this section, we discuss some well-known properties of fBm and also present several stochastic processes driven by fBm.

2.1. Fractional Brownian motion

Fractional Brownian motion ${\mathbf{B}}^H$ has the following three important properties.

1. Self-similarity. ${\mathbf{B}}^H$ is self-similar with Hurst index H, which means $B_{pt}^H$ and $p^H{B}_t^H$ have the same finite-dimensional distributions, for all p>0.

2. Regularity. Trajectories of ${\mathbf{B}}^H$ are locally Hölder continuous of any order strictly less than H, which means for every $ϵ>0$, there exists a (random) constant C such that $|B_t^H-B_s^H|\le C|t-s{|}^{H-ϵ}$.

3. Long-range dependence. The autocovariance function of the increment process, $\mathbf{Y}=\{Y_t,t\ge 0\}$, defined by $Y_t:=B_{t+1}^H-B_t^H$, reads, (2) $\begin{array}{rl}\mathrm{C}\mathrm{o}\mathrm{v}(Y_t,Y_s)& =\mathbb{E}\left[\right(B_{t+1}^H-B_t^H\left)\right(B_{s+1}^H-B_s^H\left)\right]\\ & =\frac{1}{2}\left[\right(t-s+1{)}^{2H}-2(t-s{)}^{2H}+(t-s-1{)}^{2H}],\end{array}$(2) where $0\le s\le t$. The function $\mathrm{C}\mathrm{o}\mathrm{v}(Y_t,Y_s)$ is presented, for various H, in Figure 1. $\mathrm{C}\mathrm{o}\mathrm{v}(Y_t,Y_s)=0$ for any t and s when $H=1/2$, since sBm has independent increments, whereas increments of fBm are not independent. The increments are positively correlated when $H>1/2$ and negatively correlated when $H<1/2$. Moreover, when $H>1/2$, $\mathrm{C}\mathrm{o}\mathrm{v}(Y_t,Y_s)$ decays very slowly with increasing time increment $(t-s)$, and therefore fBm exhibits long-range dependence when $H>1/2$.

Figure 1. Covariance of fBm increments for various H.

2.2. Fractional Ornstein–Uhlenbeck process

The fractional Ornstein–Uhlenbeck (fOU) process, $\mathbf{X}=\{X_t,t\ge 0\}$, is driven by fBm and is defined by the following SDE: (3) $\mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H,t\in [0,T],$(3) with an fBm ${\mathbf{B}}^H$ and given initial value $X_0$. The parameter μ represents the long-term mean value of the process, λ the speed of mean reversion, and σ the volatility. It has a unique path-wise solution as (4) $X_t=X_0{\mathrm{e}}^{-\lambda t}+\mu (1-{\mathrm{e}}^{-\lambda t})+\sigma {\int }_0^t{\mathrm{e}}^{-\lambda (t-s)}\mathrm{d}{B}_s^H.$(4) For $H=1/2$, using It$\hat{\mathrm{o}}$'s isometry property, the analytical solution of the Ornstein–Uhlenbeck process driven by standard Brownian motion reads: (5) $X_t\stackrel{d}{=}{X}_0{\mathrm{e}}^{-\lambda t}+\mu (1-{\mathrm{e}}^{-\lambda t})+\sigma \sqrt{\frac{1-{\mathrm{e}}^{-2\lambda t}}{2\lambda }}W,$(5) where W is a sample from the standard normal distribution $\mathcal{N}(0,1)$.

Using $\kappa :=H-1/2$, and based on the derivation of the conditional characteristic function, it is proved in Ref. [14] that under the filtration ${\mathcal{F}}_s(s\ge 0)$, $X_t(t\ge s)$ is normally distributed, with (6) $\begin{array}{rl}\mathbb{E}[X_t\mid {\mathcal{F}}_s]& =X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+{\int }_0^s{\mathrm{\Psi }}_c(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(6) (7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) where $c\left(r\right)=\sigma {\mathrm{e}}^{-\lambda (t-r)}$; and, for $v\in (0,s)$, (8) ${\mathrm{\Psi }}_c(s,t,v)=\frac{\sin \left(\pi \kappa \right)}{\pi }{v}^{-\kappa }(s-v{)}^{-\kappa }{\int }_s^t\frac{r^{\kappa }(r-s{)}^{\kappa }}{r-v}c(r)\mathrm{d}r,$(8) while for $v\in \{0,s\}$, ${\mathrm{\Psi }}_c(s,t,v)=0$. For $\kappa \ne 0$, the norm of a function f is defined by the fractional Riemann–Liouville integral (9) $\Vert f{\Vert }_{\kappa ,T}^2=A_{\kappa }{\int }_0^T{z}^{-2\kappa }{\left({\int }_z^T\frac{r^{\kappa }f\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z,$(9) with (10) $A_{\kappa }=\frac{\pi \kappa (2\kappa +1)}{\mathrm{\Gamma }(1-2\kappa )\sin \left(\pi \kappa \right){\mathrm{\Gamma }}^2\left(\kappa \right)},$(10) where $\mathrm{\Gamma }(\cdot )$ is the gamma function. For $\kappa =0$, $\Vert f{\Vert }_{0,T}^2$ is set as the $L_2$ norm, i.e. (11) $\Vert f{\Vert }_{0,T}^2={\int }_0^T{f}^2\left(r\right)\mathrm{d}r.$(11) The accurate numerical evaluation of Equations (6(6) $\begin{array}{rl}\mathbb{E}[X_t\mid {\mathcal{F}}_s]& =X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+{\int }_0^s{\mathrm{\Psi }}_c(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(6) ) and (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ) is the focus of the present paper. These equations also form the basis for the conditional characteristic function.

Remark 2.1

Conditional distribution of fBm

FBm can be considered a special case of the fOU process, setting $\lambda =0$ and $\sigma =1$ in Equation (3(3) $\mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H,t\in [0,T],$(3) ). Then, the conditional distribution of the fBm also is a normal distribution, with (12) $\begin{array}{rl}\mathbb{E}[B_t^H\mid {\mathcal{F}}_s]& =B_s^H+{\int }_0^s\mathrm{\Psi }(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(12) (13) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[B_t^H\mid {\mathcal{F}}_s]& =\Vert {\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert \mathrm{\Psi }(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(13) where, for $v\in (0,s)$, (14) $\mathrm{\Psi }(s,t,v)=\frac{\sin \left(\pi \kappa \right)}{\pi }{v}^{-\kappa }(s-v{)}^{-\kappa }{\int }_s^t\frac{r^{\kappa }(r-s{)}^{\kappa }}{r-v}\mathrm{d}r,$(14) for $v\in \{0,s\}$, $\mathrm{\Psi }(s,t,v)=0$. See Ref. [14].

2.3. Stochastic processes related to the fOU process

There are important processes $Z_t$ that are related to the fOU process, i.e. they can be defined based on the fOU process, as follows, (15) $Z_t=g\left(X_t\right)$(15) where $X_t$ is the fOU process and the function $g(\cdot )$ is invertible. With $z=g\left(x\right)$, the cumulative distribution functions (CDFs) of Z and X are related, as follows, (16) $F_Z\left(z\right)\stackrel{def}{=}\mathbb{P}[Z\le z]=\mathbb{P}[X\le x]\stackrel{def}{=}{F}_X\left(x\right).$(16) By differentiation, we have the relation of their probability density functions (PDFs), (17) $f_Z\left(z\right)\stackrel{def}{=}\frac{\mathrm{d}{F}_Z\left(z\right)}{\mathrm{d}z}=\frac{\mathrm{d}{F}_X\left(x\right)}{\mathrm{d}x}\frac{\mathrm{d}x}{\mathrm{d}z}\stackrel{def}{=}{\left[g^{\mathrm{\prime }}\right(x\left)\right]}^{-1}{f}_X\left(x\right).$(17) So, because the PDF of the fOU process is available, the PDF of stochastic process Z is also available by means of Equation (17(17) $f_Z\left(z\right)\stackrel{def}{=}\frac{\mathrm{d}{F}_Z\left(z\right)}{\mathrm{d}z}=\frac{\mathrm{d}{F}_X\left(x\right)}{\mathrm{d}x}\frac{\mathrm{d}x}{\mathrm{d}z}\stackrel{def}{=}{\left[g^{\mathrm{\prime }}\right(x\left)\right]}^{-1}{f}_X\left(x\right).$(17) ). We are interested in the following three processes.

1. The geometric fractional Ornstein–Uhlenbeck (GfOU) process is defined as the process whose logarithm is the fOU process, (18) $\{\begin{array}{l}{Z}_t={\mathrm{e}}^{X_t},\\ \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H.\end{array}$(18)

2. The fractional Cox–Ingersoll–Ross (fCIR) process is an extension of the classical Cox–Ingersoll–Ross (CIR) process, which is used for short-term interest rates and stochastic volatility modelling [7,17]. Compared to the CIR process, the fCIR process has the capability of modelling ‘memory’ in financial data [20], like in rough volatility. It is shown in Ref. [13] that, for $\kappa >0$, an fCIR process with zero mean, i.e. $\mathrm{d}{Z}_t=-\lambda Z_t\mathrm{d}t+\sigma \sqrt{|Z_t|}\mathrm{d}{B}_t^H,$is the square of an fOU process before hitting zero for the first time, i.e. (19) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{4}{\sigma }^2{X}_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=-\frac{\lambda }{2}{X}_t\mathrm{d}t+\mathrm{d}{B}_t^H.}\end{array}$(19)

3. Equation (19(19) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{4}{\sigma }^2{X}_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=-\frac{\lambda }{2}{X}_t\mathrm{d}t+\mathrm{d}{B}_t^H.}\end{array}$(19) ) is an example of a polynomial process. Polynomial processes are a class of processes that are obtained by a polynomial map of the underlying process [11]. This class plays an increasingly important role in finance as it provides a tractable relationship between underlying factors and resulting prices. For example, a second-order polynomial process is employed for the valuation of electricity storage contracts in Refs [5,35]. Using the same increasing polynomial map construction method as in Ref. [35], we consider a polynomial process with a third-order polynomial map, i.e. (20) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{6}\delta X_t^3+\frac{1}{2}(1-\delta )X_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H,}\end{array}$(20) where $\delta \in [0,1]$.

3. Conditional distribution of the fOU process

We reviewed the conditional distribution of the fOU process, which is a conditional normal distribution with its expectation as in Equation (6(6) $\begin{array}{rl}\mathbb{E}[X_t\mid {\mathcal{F}}_s]& =X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+{\int }_0^s{\mathrm{\Psi }}_c(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(6) ) and variance in Equation (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ). In this section, we calculate the relevant integrals by which the conditional distributions of the fOU process and several related processes can be determined for all relevant Hurst indices.

With a discretization in the time interval $[0,s]$, i.e. $0=s_0\le s_1\le \cdots \le s_n=s$, and the corresponding discrete realizations of ${\mathbf{B}}^H$ [28], the conditional expectation can be calculated by (21) $\mathbb{E}[X_t\mid {\mathcal{F}}_s]\approx X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+\sum _{i=0}^{n-1}{\mathrm{\Psi }}_c(s,t,t_i)(B_{t_{i+1}}^H-B_{t_i}^H).$(21) The variance in Equation (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ) is expressed as the difference of two terms, each being the square of a norm that is not easily calculated when $\kappa \ne 0$. The expression involves a nontrivial integration with multiple variables and singularities in its integral function. To deal with this, we simplify the expression of the conditional variance for $\kappa \ne 0$ in Section 3.1 and obtain a solution for all $\kappa \in (-1/2,1/2)$ (i.e. $H\in (0,1)$) in Section 3.2, which is summarized in Algorithm 1.

3.1. Simplification of the conditional variance

For $\kappa \ne 0$, since both terms in Equation (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ) contain indicator functions, we can rewrite them as (22) $\begin{array}{rl}\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2& =A_{\kappa }{\int }_0^T{z}^{-2\kappa }{\left({\int }_z^T\frac{r^{\kappa }c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z\\ & =A_{\kappa }[{\int }_0^s{z}^{-2\kappa }{\left({\int }_s^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z\\ & +{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z]\end{array}$(22) and (23) $\begin{array}{rl}\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2& =A_{\kappa }{\int }_0^T{z}^{-2\kappa }{\left({\int }_z^T\frac{v^{\kappa }{\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right)}{(v-z{)}^{1-\kappa }}\mathrm{d}v\right)}^2\mathrm{d}z\\ & =A_{\kappa }{\int }_0^s{z}^{-2\kappa }{G}^2\left(z\right)\mathrm{d}z,\end{array}$(23) respectively, where (24) $\begin{array}{rl}G\left(z\right)& ={\int }_z^s\frac{v^{\kappa }{\mathrm{\Psi }}_c(s,t,v)}{(v-z{)}^{1-\kappa }}\mathrm{d}v\\ & =\frac{\sin \left(\pi \kappa \right)}{\pi }{\int }_z^s(v-z{)}^{\kappa -1}(s-v{)}^{-\kappa }\left({\int }_s^t\frac{r^{\kappa }(r-s{)}^{\kappa }}{r-v}c\right(r\left)\mathrm{d}r\right)\mathrm{d}v.\end{array}$(24) By a variable transformation and the new variable, $u\in (0,\mathrm{\infty })$, (25) $u=\frac{v-z}{s-v}⟹v=\frac{us+z}{u+1},\frac{\mathrm{d}v}{\mathrm{d}u}=\frac{s-z}{(1+u{)}^2},\frac{1}{v-z}=\frac{1+u}{u(s-z)},$(25) and interchanging the order of integration, we find (26) $\begin{array}{rl}G\left(z\right)& =\frac{\sin \left(\pi \kappa \right)}{\pi }{\int }_0^{\mathrm{\infty }}\frac{u^{\kappa -1}}{1+u}\left({\int }_s^t\frac{r^{\kappa }(r-s{)}^{\kappa }}{r-v}c\right(r\left)\mathrm{d}r\right)\mathrm{d}u\\ & =\frac{\sin \left(\pi \kappa \right)}{\pi }{\int }_s^t{r}^{\kappa }(r-s{)}^{\kappa -1}c(r){\int }_0^{\mathrm{\infty }}\frac{u^{\kappa -1}}{u+\frac{r-z}{r-s}}\mathrm{d}u\mathrm{d}r.\end{array}$(26) Here, we recognize the Mellin transform (see [26, Table 1.14.4])¹ (27) ${\int }_0^{\mathrm{\infty }}\frac{u^{\kappa -1}}{u+a}\mathrm{d}u=\frac{\pi }{\sin \left(\pi \kappa \right)}{a}^{\kappa -1},0<\kappa <1.$(27) This gives us, (28) $G\left(z\right)={\int }_s^t{r}^{\kappa }(r-z{)}^{\kappa -1}c(r)\mathrm{d}r.$(28) We have thus derived this relation for $0<\kappa <1/2$, see Equation (27(27) ${\int }_0^{\mathrm{\infty }}\frac{u^{\kappa -1}}{u+a}\mathrm{d}u=\frac{\pi }{\sin \left(\pi \kappa \right)}{a}^{\kappa -1},0<\kappa <1.$(27) ), however, we need this expression, or a similar result, also for $-1/2<\kappa <0$. For this, the classical technique of analytic continuation is employed, see also Appendix 1.

Substituting Equations (28(28) $G\left(z\right)={\int }_s^t{r}^{\kappa }(r-z{)}^{\kappa -1}c(r)\mathrm{d}r.$(28) ) into (23(23) $\begin{array}{rl}\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2& =A_{\kappa }{\int }_0^T{z}^{-2\kappa }{\left({\int }_z^T\frac{v^{\kappa }{\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right)}{(v-z{)}^{1-\kappa }}\mathrm{d}v\right)}^2\mathrm{d}z\\ & =A_{\kappa }{\int }_0^s{z}^{-2\kappa }{G}^2\left(z\right)\mathrm{d}z,\end{array}$(23) ) and comparing with Equation (22(22) $\begin{array}{rl}\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2& =A_{\kappa }{\int }_0^T{z}^{-2\kappa }{\left({\int }_z^T\frac{r^{\kappa }c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z\\ & =A_{\kappa }[{\int }_0^s{z}^{-2\kappa }{\left({\int }_s^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z\\ & +{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z]\end{array}$(22) ), we find that $\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2$ is part of $\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2$ and the expression of the conditional variance in Equation (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ) can be simplified, from two terms to one term only: (29) $\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]=A_{\kappa }{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z.$(29)

3.2. Calculation of the conditional variance

Next, we illustrate how to calculate Equation (29(29) $\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]=A_{\kappa }{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z.$(29) ) for the conditional variance with $\kappa \in (-1/2,0)\cup (0,1/2)$.

Knowing that $\mathrm{\Gamma }(\kappa +1)=\kappa \mathrm{\Gamma }\left(\kappa \right)\text{and}\frac{\pi \kappa }{\sin \left(\pi \kappa \right)}=\mathrm{\Gamma }(1-\kappa )\mathrm{\Gamma }(1+\kappa ),$the front factor $A_{\kappa }$ in Equation (29(29) $\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]=A_{\kappa }{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z.$(29) ) is given by (30) $A_{\kappa }=\frac{\kappa (1-4{\kappa }^2)\mathrm{\Gamma }(1-\kappa )}{\mathrm{\Gamma }(2-2\kappa )\mathrm{\Gamma }\left(\kappa \right)}.$(30) Equation (29(29) $\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]=A_{\kappa }{\int }_s^t{z}^{-2\kappa }{\left({\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r\right)}^2\mathrm{d}z.$(29) ) can now be written as, (31) $\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]=\frac{\mathrm{\Gamma }(1-\kappa )}{\mathrm{\Gamma }(2-2\kappa )\mathrm{\Gamma }(\kappa +1)}H(s,t,\kappa ),$(31) with (32) $\begin{array}{rl}H(s,t,\kappa )& =(1-4{\kappa }^2){\int }_s^t{z}^{-2\kappa }{h}^2(\kappa ,z)\mathrm{d}z,\end{array}$(32) (33) $\begin{array}{rl}h(\kappa ,z)& =\kappa {\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r,\end{array}$(33) and (34) $c\left(r\right)=\sigma {\mathrm{e}}^{-\lambda (t-r)}.$(34)

3.2.1. Calculation of $H(s,t,\kappa )$

A set of transformations is now proposed for the calculation of $H(s,t,\kappa )$ in Equation (32(32) $\begin{array}{rl}H(s,t,\kappa )& =(1-4{\kappa }^2){\int }_s^t{z}^{-2\kappa }{h}^2(\kappa ,z)\mathrm{d}z,\end{array}$(32) ):

1. First, we use the transformation $x=z^{1-2\kappa }$, and obtain (35) $H(s,t,\kappa )=(1+2\kappa ){\int }_{s^{1-2\kappa }}^{t^{1-2\kappa }}{h}^2(\kappa ,z)\mathrm{d}x;$(35)

2. Then, using $y=(x-s^{1-2\kappa })/(t^{1-2\kappa }-s^{1-2\kappa })$ yields (36) $H(s,t,\kappa )=(1+2\kappa )(t^{1-2\kappa }-s^{1-2\kappa }){\int }_0^1{h}^2(\kappa ,z)\mathrm{d}y;$(36)

3. The final transformation is given by $y=\frac{1}{2}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(-w)$, and we get (37) $H(s,t,\kappa )=\frac{(1+2\kappa )(t^{1-2\kappa }-s^{1-2\kappa })}{\sqrt{\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}{h}^2(\kappa ,z)\mathrm{d}w,$(37) where $\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(\cdot )$ is the complementary error function (38) $\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(w\right)=\frac{2}{\sqrt{\pi }}{\int }_w^{\mathrm{\infty }}{\mathrm{e}}^{-t^2}\mathrm{d}t,$(38) with the following properties (39) $\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(-\mathrm{\infty })=2,\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\mathrm{\infty }\right)=0,\frac{\mathrm{d}y}{\mathrm{d}w}=\frac{1}{\sqrt{\pi }}{\mathrm{e}}^{-w^2}.$(39)

The resulting integral in Equation (37(37) $H(s,t,\kappa )=\frac{(1+2\kappa )(t^{1-2\kappa }-s^{1-2\kappa })}{\sqrt{\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}{h}^2(\kappa ,z)\mathrm{d}w,$(37) ) is approximated by the Trapezoidal Rule, (40) ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}f\left(w\right)\mathrm{d}w\approx m\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-(mn{)}^2}f\left(mn\right),$(40) where m>0 is the step size. Because of the fast convergence of the series in (40(40) ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}f\left(w\right)\mathrm{d}w\approx m\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-(mn{)}^2}f\left(mn\right),$(40) ), the range of summation can be cut off to a smaller interval $[-a,a]$, a>0: (41) ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}f\left(w\right)\mathrm{d}w\approx m\sum _{|mn|\le a}{\mathrm{e}}^{-(mn{)}^2}f\left(mn\right).$(41) For details about the Trapezoidal Rule, see [16, §5.4]. In our case, $f\left(w\right)=h^2(\kappa ,z)$, with w a function of z.

3.2.2. Calculation of $h(\kappa ,z)$

Here we will derive convenient expressions for the $h(\kappa ,z)$ term in Equation (33(33) $\begin{array}{rl}h(\kappa ,z)& =\kappa {\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r,\end{array}$(33) ), which is the involved part remaining.

We expand $c\left(r\right)$ from Equation (34(34) $c\left(r\right)=\sigma {\mathrm{e}}^{-\lambda (t-r)}.$(34) ) in powers of r and get an expansion with fast convergence, i.e. (42) $c\left(r\right)=\sum _{n=0}^{\mathrm{\infty }}{c}_n{r}^n,c_n=\sigma {\mathrm{e}}^{-\lambda t}\frac{{\lambda }^n}{n!}.$(42) Then, we obtain (43) $h(\kappa ,z)=\sum _{n=0}^{\mathrm{\infty }}{c}_n{R}_n(\kappa ,z),R_n(\kappa ,z)=\kappa {\int }_z^t{r}^{\kappa +n}(r-z{)}^{\kappa -1}\mathrm{d}r.$(43) Initially, we consider $\kappa >0$, because for r = z, we need a convergence of the integral. Integrating by parts gives us, (44) $R_n(\kappa ,z)={\int }_z^t{r}^{\kappa +n}\mathrm{d}(r-z{)}^{\kappa }=t^{\kappa +n}(t-z{)}^{\kappa }-(\kappa +n){\int }_z^t{r}^{\kappa +n-1}(r-z{)}^{\kappa }\mathrm{d}r.$(44) In this representation, we can accept $\kappa >-1$, so integration by parts gives us again an analytic continuation² of $R_n(\kappa ,z)$ with respect to κ.

Using $(r-z{)}^{\kappa }=(r-z{)}^{\kappa -1}(r-z)$, a simple recursion formula follows (45) $R_n=\frac{\kappa t^{\kappa +n}(t-z{)}^{\kappa }+z(\kappa +n)R_{n-1}}{2\kappa +n},n\ge 1.$(45) The functions $R_n$ can be expressed in terms of Gauss hypergeometric functions. For details about these functions, see [25]. In particular, we have the integral representation and power series, (46) $\begin{array}{rl}\text{FG}\alpha \text{β}\gamma \zeta & =\frac{\mathrm{\Gamma }\left(\gamma \right)}{\mathrm{\Gamma }\left(\text{β}\right)\mathrm{\Gamma }(\gamma -\text{β})}{\int }_0^1{t}^{\text{β}-1}(1-t{)}^{\gamma -\text{β}-1}(1-\zeta t{)}^{-\alpha }\mathrm{d}t,\\ & =\sum _{n=0}^{\mathrm{\infty }}\frac{\left(\alpha {)}_n\right(\text{β}{)}_n}{n!(\gamma {)}_n}{\zeta }^n,(\alpha {)}_n=\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }\left(\alpha \right)}.\end{array}$(46) For the integral, we assume that $\mathrm{\Re }\left\{\gamma \right\}>\mathrm{\Re }\left\{\text{β}\right\}>0$, where $\mathrm{\Re }\{\cdot \}$ denotes the real part of the input argument. Furthermore, we assume that ζ is a complex number not in the interval $[1,\mathrm{\infty })$, and for the series, we need the conditions $\left|\zeta \right|<1$ and $\gamma \ne 0,-1,-2,\dots$.

Again, initially we assume that $\kappa >0$ and use the transformation $r=z/(1-u)$, to obtain (47) $R_n(\kappa ,z)=\kappa z^{2\kappa +n}{\int }_0^{1-z/t}{u}^{\kappa -1}(1-u{)}^{-n-2\kappa -1}\mathrm{d}u.$(47) In this integral, we take $u=(1-z/t)v$, which gives us, by invoking the integral in (46(46) $\begin{array}{rl}\text{FG}\alpha \text{β}\gamma \zeta & =\frac{\mathrm{\Gamma }\left(\gamma \right)}{\mathrm{\Gamma }\left(\text{β}\right)\mathrm{\Gamma }(\gamma -\text{β})}{\int }_0^1{t}^{\text{β}-1}(1-t{)}^{\gamma -\text{β}-1}(1-\zeta t{)}^{-\alpha }\mathrm{d}t,\\ & =\sum _{n=0}^{\mathrm{\infty }}\frac{\left(\alpha {)}_n\right(\text{β}{)}_n}{n!(\gamma {)}_n}{\zeta }^n,(\alpha {)}_n=\frac{\mathrm{\Gamma }(\alpha +n)}{\mathrm{\Gamma }\left(\alpha \right)}.\end{array}$(46) ), (48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) The representation in the second line of Equation  (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ) is valid for $\kappa >-1$, and it gives us again the analytic continuation to the domain of $h(\kappa ,z)$ initially defined for $\kappa >0$.

Remark 3.1

From Equation (43(43) $h(\kappa ,z)=\sum _{n=0}^{\mathrm{\infty }}{c}_n{R}_n(\kappa ,z),R_n(\kappa ,z)=\kappa {\int }_z^t{r}^{\kappa +n}(r-z{)}^{\kappa -1}\mathrm{d}r.$(43) ), we observe that we do not need to deal with a singularity at the endpoint z = s when s>0. However, this will be needed when s = 0, in which case the variable of integration, z in Equation (32(32) $\begin{array}{rl}H(s,t,\kappa )& =(1-4{\kappa }^2){\int }_s^t{z}^{-2\kappa }{h}^2(\kappa ,z)\mathrm{d}z,\end{array}$(32) ), will reach the value z = 0, where the ${}_2{F}_1$-function in Equation  (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ) is not defined at z = 0. For this case, a different form of $R_n(\kappa ,z)$ will be required, which is given by (49) $R_n(\kappa ,z)=\frac{\mathrm{\Gamma }(\kappa +1)\mathrm{\Gamma }(n+1+\kappa )}{2\cos \left(\pi \kappa \right)\mathrm{\Gamma }(n+1+2\kappa )}{z}^{2\kappa +n}+\frac{\kappa (t-z{)}^{\kappa }{t}^{n+\kappa }}{n+2\kappa }\text{FG}-\kappa -n11-n-2\kappa \frac{z}{t}.$(49) The derivation of this form is presented in Appendix 2. The ${}_2{F}_1$-function is not defined for $z\to t$ when $\kappa >0$. In that case, Equation  (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ) should be used.

4. Option pricing methodology

In this section, we will compute conditional expectations based on the calculation of the expectation and variance of the fOU process described in the previous sections.

Usually, financial option valuation starts from the risk-neutral valuation formula, which is the discounted expectation of the asset price process under a risk-neutral measure. Risk neutrality is related to complete markets, no arbitrage and the asset price process being a martingale, i.e. the expected fair value of the discounted price process at a future time point is the present asset price. Processes based on fractional Brownian motion are not (semi-) martingales, so the risk-neutral pricing formula does not apply in this case. For financial option pricing based on non-tradable assets, the real-world probability measure is often employed. Under the assumption that such a measure has been defined, option pricing can be cast in an integral form: (50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) where $\mathrm{\Delta }t:=T-t$, V denotes the option value, x the state variable at time t and y at time T, $f(y\mid T,{\mathcal{F}}_t)$ the transition probability density, and r the interest rate.

We focus on pricing European-style options where the asset price is modelled by means of the processes described in Section 2.3. For these processes, the transition probability, $f(y\mid T,{\mathcal{F}}_t)$ in Equation (50(50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) ), is the PDF of a normal distribution with expectation ${\mu }_y=\mathbb{E}[X_T\mid {\mathcal{F}}_t]$ and variance ${\sigma }_y^2=\mathrm{V}\mathrm{a}\mathrm{r}[X_T\mid {\mathcal{F}}_t]$, as in Equations (6(6) $\begin{array}{rl}\mathbb{E}[X_t\mid {\mathcal{F}}_s]& =X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+{\int }_0^s{\mathrm{\Psi }}_c(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(6) ) and (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ), respectively. It is therefore possible to obtain a closed-form solution of Equation (50(50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) ) for some specific processes. For example, if the asset price follows a GfOU process, see Equation (18(18) $\{\begin{array}{l}{Z}_t={\mathrm{e}}^{X_t},\\ \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H.\end{array}$(18) ), and the payoff reads (51) $V(T,y)=max[\eta ({\mathrm{e}}^y-K),0]\text{with}\eta =\{\begin{array}{ll}1& \text{for a call},\\ -1& \text{for a put},\end{array}$(51) with the strike price K, the closed-form solution of Equation (50(50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) ) for a call option can be derived as (52) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t+\frac{1}{2}{\sigma }_y^2+{\mu }_y}[1-\mathrm{\Phi }\left(\frac{\log K-{\sigma }_y^2-{\mu }_y}{{\sigma }_y}\right)]-K{\mathrm{e}}^{-r\mathrm{\Delta }t}[1-\mathrm{\Phi }\left(\frac{\log K-{\mu }_y}{{\sigma }_y}\right)],$(52) and for a put option, (53) $V(t,x)=K{\mathrm{e}}^{-r\mathrm{\Delta }t}\mathrm{\Phi }\left(\frac{\log K-{\mu }_y}{{\sigma }_y}\right)-{\mathrm{e}}^{-r\mathrm{\Delta }t+\frac{1}{2}{\sigma }_y^2+{\mu }_y}\mathrm{\Phi }\left(\frac{\log K-{\sigma }_y^2-{\mu }_y}{{\sigma }_y}\right).$(53) Here, $\mathrm{\Phi }(\cdot )$ denotes the CDF of the standard normal distribution.

For more general processes, e.g. polynomial processes, a closed-form option value expression is difficult or even impossible to obtain. It is, however, possible to employ the COS method in such cases.

4.1. The COS method

The COS method  [9,10] can be applied to compute expectations of (underlying) processes for which the characteristic function is available. For Markov processes with affine dynamics, the characteristic function can be derived by solving the Ricatti differential equations [8]. Affinity is not invariant under polynomial transformations [12] and, therefore, the characteristic function of a stochastic model which is based on a polynomial process, $Z_t=g\left(X_t\right)$, generally does not exist. However, by using a transformation $z=g\left(x\right)$ in the definition of the COS method, the state variables can be appropriately defined so that the characteristic function of the underlying process $X_t$, (54) $\varphi (u\mid T,{\mathcal{F}}_t)={\mathrm{e}}^{i{\mu }_{y}u-\frac{1}{2}{\sigma }_y^2{u}^2},$(54) can be employed within the pricing formula, with, again, ${\mu }_y=\mathbb{E}[X_T\mid {\mathcal{F}}_t]$ and ${\sigma }_y^2=\mathrm{V}\mathrm{a}\mathrm{r}[X_T\mid {\mathcal{F}}_t]$, as in Equations (6(6) $\begin{array}{rl}\mathbb{E}[X_t\mid {\mathcal{F}}_s]& =X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+{\int }_0^s{\mathrm{\Psi }}_c(s,t,v)\mathrm{d}{B}_v^H,\end{array}$(6) ) and (7(7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}[X_t\mid {\mathcal{F}}_s]& =\Vert c\left(r\right){\mathbf{1}}_{[s,t]}\left(r\right){\Vert }_{\kappa ,T}^2-\Vert {\mathrm{\Psi }}_c(s,t,v){\mathbf{1}}_{[0,s]}\left(v\right){\Vert }_{\kappa ,T}^2,\end{array}$(7) ).

The conditional probability density function in Equation (50(50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) ) is approximated by a truncated Fourier cosine expansion, which uses the characteristic function to recover the Fourier coefficients, as follows [9]: (55) $f(y\mid T,{\mathcal{F}}_t)\approx \frac{2}{d-b}\sum _{l=0}^{L-1}{}^{\mathrm{\prime }}\mathrm{\Re }\{\varphi (\frac{l\pi }{d-b}|T,{\mathcal{F}}_t)\cdot {\mathrm{e}}^{-il\pi \frac{b}{d-b}}\}\cos \left(l\pi \frac{y-b}{d-b}\right),$(55) where $[b,d]$ is the integration interval, L is the number of terms in the Fourier-cosine expansion and ${\sum }^{\prime }$ implies that the first term of the summation is multiplied by $\frac{1}{2}$.

By replacing the conditional density function $f(y\mid T,{\mathcal{F}}_t)$ in Equation (50(50) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t}\mathbb{E}\left[V\right(T,y)\mid {\mathcal{F}}_t]={\mathrm{e}}^{-r\mathrm{\Delta }t}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}V(T,y)f(y\mid T,{\mathcal{F}}_t)\mathrm{d}y,$(50) ) by its Fourier cosine expansion approximation Equation (55(55) $f(y\mid T,{\mathcal{F}}_t)\approx \frac{2}{d-b}\sum _{l=0}^{L-1}{}^{\mathrm{\prime }}\mathrm{\Re }\{\varphi (\frac{l\pi }{d-b}|T,{\mathcal{F}}_t)\cdot {\mathrm{e}}^{-il\pi \frac{b}{d-b}}\}\cos \left(l\pi \frac{y-b}{d-b}\right),$(55) ) and interchanging integration and summation, with the help of Fubini's theorem, the COS formula is obtained: (56) $V(t,x)\approx {\mathrm{e}}^{-r\mathrm{\Delta }t}\sum _{l=0}^{L-1}{}^{\mathrm{\prime }}\mathrm{\Re }\left\{\varphi (\frac{l\pi }{d-b}|T,{\mathcal{F}}_t){\mathrm{e}}^{il\pi \frac{-b}{d-b}}\right\}{\mathcal{V}}_l,$(56) where coefficients ${\mathcal{V}}_l$ are defined by: (57) ${\mathcal{V}}_l=\frac{2}{d-b}{\int }_b^{d}V(T,y)\cos \left(l\pi \frac{y-b}{d-b}\right)\mathrm{d}y,$(57) where, importantly, the state variables x and y can be any function of respectively the asset prices $Z_t$ and $Z_T$, e.g. in our case, $x=g^{-1}\left(Z_t\right)=X_t$ and $y=g^{-1}\left(Z_T\right)=X_T$, with $g(\cdot )$, invertible, as in Equation (15(15) $Z_t=g\left(X_t\right)$(15) ). So, we can work with the characteristic function of the fOU process $X_t$, but price options under the more involved process $Z_t$. A closed-form solution of the coefficients ${\mathcal{V}}_l$ is available for various payoff functions and several choices of the state variables x and y [9].

5. Numerical results

We show numerical results for the methods presented, starting with the calculation of the conditional variances of the fBm and fOU processes, followed by the computation of the probability density functions of the GfOU, fCIR and polynomial processes, and some option pricing results using the COS method. These results show the accuracy and efficiency of the proposed computational techniques.

5.1. Calculation of variance

First of all, we focus on the accurate representation of the conditional variance by means of the integrals that are derived in Section 3.

5.1.1. Choice of hyperparameters

As explained in the previous sections and also shown in Algorithm 1, there are three hyperparameters in the numerical method to compute the integrals in the approximation of the conditional variance for the fOU process, i.e. the step size m, summation range a for the Trapezoidal Rule in Equation (41(41) ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-w^2}f\left(w\right)\mathrm{d}w\approx m\sum _{|mn|\le a}{\mathrm{e}}^{-(mn{)}^2}f\left(mn\right).$(41) ), and the number of expansion terms N in Equation (43(43) $h(\kappa ,z)=\sum _{n=0}^{\mathrm{\infty }}{c}_n{R}_n(\kappa ,z),R_n(\kappa ,z)=\kappa {\int }_z^t{r}^{\kappa +n}(r-z{)}^{\kappa -1}\mathrm{d}r.$(43) ). To determine suitable values for these parameters, we perform the following experiment. For two values of H, and the value of N sufficiently large (i.e. N = 20), the computed results by varying m and a are presented in Figure 2(a). The figure shows that the summation range value a has a bigger impact than the step size m. The numerical approximation converges when a is larger than 3 and m is smaller than 0.6. We therefore set a = 5 and m = 0.5 in the subsequent numerical experiments. Next, we consider the effect of the number of expansion terms N on the computed results, also for the two values of H, see Figure 2(b). The results show that the numerical scheme converges rapidly as N increases, and setting N>10 already provides reliable computational results. We may define the stopping criterion by checking when $\left|c_n{R}_n/h\right|$ is less than ${10}^{-8}$, then stop the summation. The results presented are representative for a large range of model parameters.

Figure 2. Calculated variances with $\sigma =0.3$, $\lambda =0.5$, s = 0 and t = 5 for two values of H. (a) Step size m varies from 0.1 to 1.0, and summation range a from 1 to 5, with fixed N = 20. (b) N varies from 1 to 20, with m = 0.5 and a = 5 fixed. (a) m and a; (b) N.

5.1.2. Numerical performance

To analyse the accuracy of the numerical technique, we compare the conditional (s = 3) and unconditional (s = 0) variances obtained by our method with Monte Carlo simulation. FBm and fOU, with $\sigma =0.3$, $\lambda =0.5$, are considered, with different Hurst indices H. For the simulation of fBm, we use a Python package named ‘fbm’³.

Two sets of parameters are chosen, one with a step size of 0.01 and the number of Monte Carlo paths of ${10}^4$, and another with the step size halved and the number of paths multiplied by 10. We report the relative errors in Table 1, which are defined as (58) $ϵ=\frac{|{\overline{\sigma }}_1-{\overline{\sigma }}_2|}{{\overline{\sigma }}_1},$(58) where the standard deviation of the Monte Carlo results is ${\overline{\sigma }}_1$ and the computed standard deviation by numerical integration is ${\overline{\sigma }}_2$.

Table 1. Conditional (s = 3) and unconditional (s = 0) standard deviations of fBm and the fOU process, with $\sigma =0.3$, $\lambda =0.5$, t = s + 5 and Hurst index H from 0.1 to 0.9.

	fBm		fOU
H	s = 0	s = 3	s = 0	s = 3
0.1	1.1662 (0.92, 0.75)	0.9774 (0.54, 0.32)	0.2186 (0.80, 0.22)	0.2175 (0.82, 0.25)
0.2	1.3796 (0.59, 0.19)	1.2546 (0.41, 0.15)	0.2310 (0.72, 0.16)	0.2296 (0.42, 0.16)
0.3	1.6207 (0.44, 0.18)	1.5544 (0.60, 0.15)	0.2482 (0.46, 0.18)	0.2470 (0.40, 0.19)
0.4	1.9036 (0.65, 0.22)	1.8832 (0.43, 0.12)	0.2708 (0.66, 0.17)	0.2702 (0.47, 0.23)
0.5	2.2361 (0.52, 0.18)	2.2361 (0.63, 0.23)	0.2990 (0.61, 0.17)	0.2990 (0.76, 0.17)
0.6	2.6265 (0.99, 0.21)	2.5924 (0.68, 0.14)	0.3334 (0.80, 0.21)	0.3317 (0.62, 0.16)
0.7	3.0852 (0.55, 0.19)	2.9025 (0.55, 0.14)	0.3746 (0.53, 0.18)	0.3633 (0.53, 0.25)
0.8	3.6239 (0.31, 0.17)	3.0555 (0.52, 0.23)	0.4238 (0.67, 0.21)	0.3808 (0.54, 0.19)
0.9	4.2568 (0.41, 0.11)	2.7760 (0.50, 0.22)	0.4822 (0.57, 0.13)	0.3491 (0.61, 0.21)

Notes: The numbers in brackets are relative errors (%) compared with Monte Carlo simulations with two configurations, i.e. time step 0.01, number of paths ${10}^4$, and time step 0.005, number of paths ${10}^5$, respectively. Errors reported are the averages of 10 experiments.

As shown in Table 1, our numerical technique provides highly accurate results very close to the Monte Carlo simulation results, with average relative errors of less than $1\mathrm{\%}$. By using a smaller time step and simulating more paths within the Monte Carlo simulation, the differences are even smaller. The CPU times⁴ are also compared in Table 2. Clearly, a numerical integration technique costs much less computing time than a Monte Carlo sampling method, because Monte Carlo simulation requires the generation of many paths with a small time step to obtain reliable results.

Table 2. CPU times (s) using our computational technique and Monte Carlo simulations for the variances of fBm and the fOU process with three representative H indices (i.e. $H=0.3,0.5,0.7$).

	fBm		fOU
H	s = 0	s = 3	s = 0	s = 3
0.3	$1.1\times {10}^{-3}$, 1.30, 20.05	$2.4\times {10}^{-3}$, 2.13, 65.74	$1.1\times {10}^{-3}$, 1.47, 21.47	$2.4\times {10}^{-3}$, 2.37, 63.12
0.5	$1.2\times {10}^{-5}$, 0.44, 8.58	$6.7\times {10}^{-6}$, 0.58, 11.68	$1.1\times {10}^{-5}$, 0.44, 8.50	$7.7\times {10}^{-6}$, 0.54, 11.70
0.7	$1.1\times {10}^{-3}$, 1.30, 20.27	$2.4\times {10}^{-3}$, 2.12, 65.67	$1.1\times {10}^{-3}$, 1.54, 21.42	$2.6\times {10}^{-4}$, 2.15, 62.91

Notes: The triples are CPU times of our computational technique, Monte Carlo simulations with time step 0.01, number of paths ${10}^4$ and time step 0.005, number of paths ${10}^5$, respectively. Results reported are the averages of 10 experiments.

We further present the computed variances for different values of H, time points s and with increasing sizes of the time interval t−s in Figure 3 (fBm) and Figure 4 (fOU process with $\sigma =0.3$, $\lambda =0.5$). From these figures, we conclude that our computational technique converges very well for all relevant H indices, in a robust way. The curves for the different values of H are very accurate, as compared to the literature [14]. As the time interval t−s increases, the variance increases too. With a sufficiently large time interval t−s, large values of the Hurst parameter H give rise to higher variances. This makes sense as the Hurst index indicates the roughness of a process, with a higher value leading to smoother paths. Comparing the variances for s = 0 and s = 3, we see different results for the same time interval t−s, which confirms that the variance is not a function t−s, but of s and t separately.

Figure 3. Variances of the fBm process for different H and t values. (a) s = 0; (b) s = 3.

Figure 4. Variances of the fOU process for different H and t, and $\sigma =0.3$, $\lambda =0.5$. (a) s = 0; (b) s = 3.

5.2. PDF recovery of fOU-related processes

As described in Section 2.3, due to the relation between the PDFs of the fOU process and the processes in Equation (17(17) $f_Z\left(z\right)\stackrel{def}{=}\frac{\mathrm{d}{F}_Z\left(z\right)}{\mathrm{d}z}=\frac{\mathrm{d}{F}_X\left(x\right)}{\mathrm{d}x}\frac{\mathrm{d}x}{\mathrm{d}z}\stackrel{def}{=}{\left[g^{\mathrm{\prime }}\right(x\left)\right]}^{-1}{f}_X\left(x\right).$(17) ), we can compute the PDFs of the related processes by means of the PDF of the fOU process, which is normally distributed with the expectation in Equation (21(21) $\mathbb{E}[X_t\mid {\mathcal{F}}_s]\approx X_s{\mathrm{e}}^{-\lambda (t-s)}+\mu (1-{\mathrm{e}}^{-\lambda (t-s)})+\sum _{i=0}^{n-1}{\mathrm{\Psi }}_c(s,t,t_i)(B_{t_{i+1}}^H-B_{t_i}^H).$(21) ) and the variance computed as in Algorithm 1, and the transformation $g(\cdot )$ in Equation (15(15) $Z_t=g\left(X_t\right)$(15) ).

Figure 5 presents the PDFs of the GfOU, fCIR and polynomial process, as in Equations (18(18) $\{\begin{array}{l}{Z}_t={\mathrm{e}}^{X_t},\\ \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H.\end{array}$(18) ), (19(19) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{4}{\sigma }^2{X}_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=-\frac{\lambda }{2}{X}_t\mathrm{d}t+\mathrm{d}{B}_t^H.}\end{array}$(19) ) and (20(20) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{6}\delta X_t^3+\frac{1}{2}(1-\delta )X_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H,}\end{array}$(20) ), respectively. We calculate the PDFs on the condition of either initial value (s = 0, t = 3) or realized historical paths (s = 3, t = 6), and compare them with the results by Monte Carlo simulation. It is shown that the calculated PDFs resemble the histograms of Monte Carlo simulation very well, which is an indication of the accuracy of our technique. Moreover, by comparing the first and third columns of Figure 5, we can see the impact of the memory of the processes (the second column of Figure 5) on the PDFs.

Figure 5. PDFs of the related processes with H = 0.75, $\sigma =0.3$, $\lambda =0.5$, $Z_0=10$. The red curves are the PDFs obtained by our scheme, the blue histograms are Monte Carlo simulation results. From left to right, first column: PDFs at time t = 3 with given initial value; second column: realized paths until time s = 3; third column, PDFs at time t = 6 on the condition of corresponding path in the second column. From top to bottom, first row: the GfOU process as in Equation (18(18) $\{\begin{array}{l}{Z}_t={\mathrm{e}}^{X_t},\\ \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H.\end{array}$(18) ) with $\mu =\log \left(Z_0\right)$; second row: the fCIR process as in Equation (19(19) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{4}{\sigma }^2{X}_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=-\frac{\lambda }{2}{X}_t\mathrm{d}t+\mathrm{d}{B}_t^H.}\end{array}$(19) ); third row: the polynomial process as in Equation (20(20) $\{\begin{array}{l}{\displaystyle Z_t=\frac{1}{6}\delta X_t^3+\frac{1}{2}(1-\delta )X_t^2,}\\ {\displaystyle \mathrm{d}{X}_t=\lambda (\mu -X_t)\mathrm{d}t+\sigma \mathrm{d}{B}_t^H,}\end{array}$(20) ) with $\delta =0.8$ and $\mu =g^{-1}\left(Z_0\right)$, where $g\left(x\right)=\delta x^3/6+(1-\delta )x^2/2$.

5.3. Option pricing

The COS method enables us to price options if we have the characteristic function available for a function of the asset price process, like $X_t$ in $Z_t=g\left(X_t\right)$. When the asset prices are modelled as one of the fOU-related processes, the use of the COS method is based on the characteristic function of the fOU process $X_t$.

To test the performance of the COS method, we price a European call option under the GfOU process. This is just a reference example, since a closed-form solution is available in this case, as in Equation (52(52) $V(t,x)={\mathrm{e}}^{-r\mathrm{\Delta }t+\frac{1}{2}{\sigma }_y^2+{\mu }_y}[1-\mathrm{\Phi }\left(\frac{\log K-{\sigma }_y^2-{\mu }_y}{{\sigma }_y}\right)]-K{\mathrm{e}}^{-r\mathrm{\Delta }t}[1-\mathrm{\Phi }\left(\frac{\log K-{\mu }_y}{{\sigma }_y}\right)],$(52) ). The errors of the COS method are presented in Table 3 with the strike price K = 10, the interval of integration $[b,d]=[{\mu }_y-10\sqrt{{\sigma }_y^2},{\mu }_y+10\sqrt{{\sigma }_y^2}]$ and a varying number of cosine terms L, see Equations (55(55) $f(y\mid T,{\mathcal{F}}_t)\approx \frac{2}{d-b}\sum _{l=0}^{L-1}{}^{\mathrm{\prime }}\mathrm{\Re }\{\varphi (\frac{l\pi }{d-b}|T,{\mathcal{F}}_t)\cdot {\mathrm{e}}^{-il\pi \frac{b}{d-b}}\}\cos \left(l\pi \frac{y-b}{d-b}\right),$(55) ) and (56(56) $V(t,x)\approx {\mathrm{e}}^{-r\mathrm{\Delta }t}\sum _{l=0}^{L-1}{}^{\mathrm{\prime }}\mathrm{\Re }\left\{\varphi (\frac{l\pi }{d-b}|T,{\mathcal{F}}_t){\mathrm{e}}^{il\pi \frac{-b}{d-b}}\right\}{\mathcal{V}}_l,$(56) ). It is shown that the COS method converges rapidly, as with only L = 16 Fourier-cosine expansion terms, the errors are of the order ${10}^{-8}$.

Table 3. Errors of the COS method for a European call option under the GfOU process with K = 10, t = 0, T = 3, r = 0.1.

	H
L	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
4	$8.4\times {10}^{-02}$	$9.0\times {10}^{-02}$	$9.8\times {10}^{-02}$	$1.1\times {10}^{-01}$	$1.2\times {10}^{-01}$	$1.3\times {10}^{-01}$	$1.5\times {10}^{-01}$	$1.6\times {10}^{-01}$	$1.8\times {10}^{-01}$
8	$1.8\times {10}^{-03}$	$2.0\times {10}^{-03}$	$2.1\times {10}^{-03}$	$2.4\times {10}^{-03}$	$2.6\times {10}^{-03}$	$2.9\times {10}^{-03}$	$3.3\times {10}^{-03}$	$3.7\times {10}^{-03}$	$4.1\times {10}^{-03}$
16	$2.8\times {10}^{-08}$	$2.5\times {10}^{-08}$	$2.1\times {10}^{-08}$	$1.6\times {10}^{-08}$	$9.1\times {10}^{-09}$	$1.0\times {10}^{-09}$	$8.6\times {10}^{-09}$	$2.0\times {10}^{-08}$	$3.2\times {10}^{-08}$

We further present the prices of European call and put options under the GfOU process in Figure 6 and the polynomial process in Figure 7, respectively, for different Hurst indices H and strike prices K. Obviously, the Hurst index H has a significant impact on the option price, which indicates the importance of modelling the asset price as a stochastic process with a suitable H.

Figure 6. European option values under the GfOU process with different H and K, t = 0 and T = 3 fixed. (a) Call option; (b) Put option.

Figure 7. European option values under the polynomial process with different H and K, t = 0 and T = 3 fixed. (a) Call option; (b) Put option.

6. Conclusion

We have shown that, on the basis of several transformations of variables and the analytic continuation technique, integrals that define the conditional expectation, variance and characteristic function of stochastic processes with respect to fractional Brownian motion (fBm) can be reformulated for all relevant Hurst indices. Numerical experiments confirm the robustness and efficiency of the integral formulations and of the numerical technique regarding different values of the Hurst index. Moreover, by means of the COS method, we have presented accurate and highly efficient option pricing results for processes connected to the fOU process, on the basis of the conditional characteristic function of the fOU process.

The methodology presented can be used as an alternative to Monte Carlo simulation when appropriate, or as a validation technique for implementations of other numerical techniques.

Acknowledgements

N. M. Temme thanks CWI, Amsterdam, for general support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

F. Gao would like to thank the China Scholarship Council [CSC, grant number 202006280439] for the financial support. N. M. Temme acknowledges financial support from Spanish Ministry of Science and Innovation (Ministerio de Ciencia e Innovación), Spain, project MTM2012-11686.

Notes

1 http://dlmf.nist.gov/1.14.T4

2 For details on analytic continuation, we refer again to Appendix 1.

3 https://pypi.org/project/fbm/

4 The computer used has an Intel Core i7-8700K CPU. The code is written in Python 3.

Appendices

Appendix 1.

Analytic continuation

In complex analysis, analytic continuation is a classical technique to extend the domain of definition of a given analytic function. Analytic continuation often succeeds in defining further values of a function, for example, in a new region where an infinite series representation, in terms of which it is initially defined, becomes divergent.

As a standard form to explain analytic continuation for a class of integrals, we consider the function $F\left(x\right)$ defined by $F\left(x\right)=\frac{1}{\mathrm{\Gamma }\left(x\right)}{\int }_0^{\mathrm{\infty }}{t}^{x-1}{\mathrm{e}}^{-t}f\left(t\right)\mathrm{d}t,\mathrm{\Re }\left\{x\right\}>0,$where the function $f\left(t\right)$ is analytic in $[0,\mathrm{\infty })$, and in a complex domain around this interval. $F\left(x\right)$ is an analytic function in the domain $\mathrm{\Re }\left\{x\right\}>0$. When we assume that $f\left(0\right)\ne 0$, the integral does not converge if x = 0. The reciprocal gamma function in front of the integral will allow us to give an interpretation to $F\left(x\right)$ for $x\le 0$, but a different choice of a front factor may act in a similar way. By a standard method, we can change the representation of $F\left(x\right)$ into a different form from which we can conclude that $F\left(x\right)$ is analytic in a larger domain, in particular at x = 0.

We demonstrate the idea in detail by using the example $F\left(x\right)=\frac{1}{\mathrm{\Gamma }\left(x\right)}{\int }_0^1{t}^{x-1}\cos \left(at\right)\mathrm{d}t,\mathrm{\Re }\left\{x\right\}>0,$where $a\ne 0$ is some real or complex parameter. Here, the condition $\mathrm{\Re }\left\{x\right\}>0$ controls the convergence at the origin.

Integration by parts gives, by writing $t^{x-1}=\frac{1}{x}\frac{\mathrm{d}{t}^x}{\mathrm{d}x}$ and using $x\mathrm{\Gamma }\left(x\right)=\mathrm{\Gamma }(x+1)$, $F\left(x\right)=\frac{\cos a}{\mathrm{\Gamma }(x+1)}+\frac{a}{\mathrm{\Gamma }(x+1)}{\int }_0^1{t}^x\sin \left(at\right)\mathrm{d}t$The new integral converges inside the domain $\mathrm{\Re }\left\{x\right\}>-2$.

We have two representations of $F\left(x\right)$. In the overlapping domain, $\mathrm{\Re }\left\{x\right\}>0$, these representations give the same analytic function $F\left(x\right)$, in the new domain $\mathrm{\Re }\left\{x\right\}>-2$ the second representation is analytic, and the principle of analytic continuation gives that, hence, $F\left(x\right)$ can be analytically extended from all points in $\mathrm{\Re }\left\{x\right\}>0$ to all points in $\mathrm{\Re }\left\{x\right\}>-2$. The first representation cannot be used in that new domain.

The integration by parts technique can be continued, and it gives the result that $F\left(x\right)$ is analytic in a domain $\mathrm{\Re }\left\{x\right\}>-n$ for each nonnegative integer n.

As a bonus, we observe that by expanding the cosine function as a power series, we obtain for $\mathrm{\Re }\left\{x\right\}>0$ $\begin{array}{rl}F\left(x\right)& =\frac{1}{\mathrm{\Gamma }\left(x\right)}{\int }_0^1{t}^{x-1}\sum _{n=0}^{\mathrm{\infty }}(-1{)}^n\frac{(at{)}^{2n}}{\left(2n\right)!}\mathrm{d}t\\ & =\frac{1}{\mathrm{\Gamma }\left(x\right)}\sum _{n=0}^{\mathrm{\infty }}(-1{)}^n\frac{a^{2n}}{\left(2n\right)!}{\int }_0^1{t}^{x-1+2n}\mathrm{d}t\\ & =\sum _{n=0}^{\mathrm{\infty }}(-1{)}^n\frac{a^{2n}}{\mathrm{\Gamma }\left(x\right)(2n+x)\left(2n\right)!}.\end{array}$Because for $n=0,1,2,,\dots$ $\frac{1}{\mathrm{\Gamma }\left(x\right)(x+2n)}=\frac{\mathrm{\Gamma }(1-x)}{\mathrm{\Gamma }(1+x+2n)\mathrm{\Gamma }(1-x-2n)}\to \left(2n\right)!,$as $x\to -2n$, we see that the series gives a representation of $F\left(x\right)$ as an analytic function for all complex values of x.

For further details on analytic continuation, we refer to the references [31, Section 1.10(ii)] and [32, Chapter IV].

Appendix 2.

The case s = 0

In (32(32) $\begin{array}{rl}H(s,t,\kappa )& =(1-4{\kappa }^2){\int }_s^t{z}^{-2\kappa }{h}^2(\kappa ,z)\mathrm{d}z,\end{array}$(32) ), we see that we do not need to handle a singularity at the endpoint z = s when s>0. However, this will be needed when s = 0, in which case the variable of integration, z in (33(33) $\begin{array}{rl}h(\kappa ,z)& =\kappa {\int }_z^t\frac{r^{\kappa }c\left(r\right)}{(r-z{)}^{1-\kappa }}\mathrm{d}r,\end{array}$(33) ), will take the value z = 0.

To deal with this, first of all, we write the representation of $R_n(\kappa ,z)$ in (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ) in different form because the ${}_2{F}_1$-function is not defined at z = 0 (although the front factor $z^{2\kappa +n}$ will control the product of the two quantities).

We use (A1) $w_1\left(z\right)=A{w}_3\left(z\right)+B{w}_4\left(z\right)$(A1) of [25, Eqn. 15.10.21], with $w_3\left(z\right)$ and $w_4\left(z\right)$ in the first lines of Equations (15.10.13) and (15.10.14) of that reference. Here, (A2) $A=\frac{\mathrm{\Gamma }(\kappa +1)\mathrm{\Gamma }(-n-2\kappa )}{\mathrm{\Gamma }(-n-\kappa )}=\frac{\mathrm{\Gamma }(\kappa +1)\mathrm{\Gamma }(n+1+\kappa )}{2\cos \left(\pi \kappa \right)\mathrm{\Gamma }(n+1+2\kappa )},$(A2) and (A3) $B=\frac{\mathrm{\Gamma }(\kappa +1)\mathrm{\Gamma }(n+2\kappa )}{\mathrm{\Gamma }(n+2\kappa +1)\mathrm{\Gamma }\left(\kappa \right)}=\frac{\kappa }{n+2\kappa }.$(A3) For $w_3\left(z\right)$, we have an elementary form, i.e. (A4) $w_3\left(z\right)=\text{FG}n+2\kappa +1\kappa n+2\kappa +1\frac{z}{t}=(1-z/t{)}^{-\kappa },$(A4) and $w_4\left(z\right)$ is given by (A5) $w_4\left(z\right)={\left(\frac{z}{t}\right)}^{-n-2\kappa }\text{FG}-\kappa -n11-n-2\kappa \frac{z}{t}.$(A5) This gives us, for (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ), (A6) $R_n(\kappa ,z)=\frac{\mathrm{\Gamma }(\kappa +1)\mathrm{\Gamma }(n+1+\kappa )}{2\cos \left(\pi \kappa \right)\mathrm{\Gamma }(n+1+2\kappa )}{z}^{2\kappa +n}+\frac{\kappa (t-z{)}^{\kappa }{t}^{n+\kappa }}{n+2\kappa }\text{FG}-\kappa -n11-n-2\kappa \frac{z}{t}.$(A6) This ${}_2{F}_1$-function is not defined for $z\to t$ when $\kappa >0$. In that case Equation  (48(48) $\begin{array}{rl}{R}_n(\kappa ,z)& =\kappa z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }{\int }_0^1{v}^{\kappa -1}{[1-(1-\frac{z}{t})v]}^{-n-2\kappa -1}\mathrm{d}v\\ & =z^{2\kappa +n}{(1-\frac{z}{t})}^{\kappa }\text{FG}n+2\kappa +1\kappa \kappa +11-\frac{z}{t}.\end{array}$(48) ) should be used.