Option pricing in fractional Heston-type model

TL;DR

This paper develops discretization schemes for option pricing in a fractional Heston model with H>1/2, analyzes their convergence, and employs Malliavin calculus to handle discontinuous payoffs.

Contribution

It introduces new discretization methods for the fractional Heston model and derives convergence rates, also applying Malliavin calculus for discontinuous payoffs.

Findings

Discretization schemes converge as mesh size decreases.

Convergence rate of the expectation approximation is established.

Malliavin calculus provides an alternative expectation formula for discontinuous payoffs.

Abstract

In this paper, we consider option pricing in a framework of the fractional Heston-type model with $H>1/2$. As it is impossible to obtain an explicit formula for the expectation $\mathbb E f(S_T)$ in this case, where $S_T$ is the asset price at maturity time and $f$ is a payoff function, we provide a discretization schemes $\hat Y^n$ and $\hat S^n$ for volatility and price processes correspondingly and study convergence $\mathbb E f(\hat S^n_T) \to \mathbb E f(S_T)$ as the mesh of the partition tends to zero. The rate of convergence is calculated. As we allow $f$ to have discontinuities of the first kind which can cause errors in straightforward Monte-Carlo estimation of the expectation, we use Malliavin calculus techniques to provide an alternative formula for $\mathbb E f(S_T)$ with smooth functional under the expectation.

Tables

Table 1. Table 1. f ​ ( x ) = ( x − K ) + 𝑓 𝑥 superscript 𝑥 𝐾 f(x)=(x-K)^{+} , K = 1 𝐾 1 K=1 , σ ​ ( x ) = 0.5 ​ ( x + 0.01 ) 0.9 𝜎 𝑥 0.5 superscript 𝑥 0.01 0.9 \sigma(x)=0.5(x+0.01)^{0.9} , μ = 0.5 𝜇 0.5 \mu=0.5 , H = 0.7 𝐻 0.7 H=0.7

$n$	Mean	Standard deviation	Coefficient of variation	Min.	1st Qu.	Median	3rd Qu.	Max.
100	0.7019	0.05628101	0.0802	0.5171	0.6630	0.7006	0.7380	0.8989
500	0.7040	0.05476103	0.0778	0.5406	0.6655	0.7025	0.7406	0.9305
1000	0.7004	0.05459163	0.0779	0.5463	0.6625	0.6978	0.7375	0.9344

Table 2. Table 2. f ​ ( x ) = 𝟙 [ 0.5 , 1 ] ​ ( x ) 𝑓 𝑥 subscript 1 0.5 1 𝑥 f(x)=\mathbbm{1}_{[0.5,1]}(x) , σ ​ ( x ) = 0.5 ​ ( x + 0.01 ) 0.9 𝜎 𝑥 0.5 superscript 𝑥 0.01 0.9 \sigma(x)=0.5(x+0.01)^{0.9} , μ = 0.5 𝜇 0.5 \mu=0.5 , H = 0.7 𝐻 0.7 H=0.7

$n$	Mean	Standard deviation	Coefficient of variation	Min.	1st Qu.	Median	3rd Qu.	Max.
100	0.2126	0.01196734	0.0563	0.1790	0.2046	0.2131	0.2206	0.2518
500	0.2123	0.01266216	0.0596	0.1652	0.2037	0.2124	0.2206	0.2553
1000	0.2129	0.01272749	0.0598	0.1725	0.2042	0.2132	0.2210	0.2505

Table 3. Table 3. f ​ ( x ) = 𝟙 ( 0.5 , ∞ ) ​ ( x ) + 1 2 ​ ∑ k = 2 6 𝟙 ( 0.5 ​ k , ∞ ) ​ ( x ) 𝑓 𝑥 subscript 1 0.5 𝑥 1 2 superscript subscript 𝑘 2 6 subscript 1 0.5 𝑘 𝑥 f(x)=\mathbbm{1}_{(0.5,\infty)}(x)+\frac{1}{2}\sum_{k=2}^{6}\mathbbm{1}_{(0.5k,\infty)}(x) , σ ​ ( x ) = 0.5 ​ ( x + 0.01 ) 0.9 𝜎 𝑥 0.5 superscript 𝑥 0.01 0.9 \sigma(x)=0.5(x+0.01)^{0.9} , μ = 0.5 𝜇 0.5 \mu=0.5 , H = 0.7 𝐻 0.7 H=0.7

$n$	Mean	Standard deviation	Coefficient of variation	Min.	1st Qu.	Median	3rd Qu.	Max.
100	1.804	0.08973507	0.0497	1.476	1.748	1.803	1.864	2.066
500	1.806	0.08873267	0.0491	1.546	1.745	1.805	1.866	2.136
1000	1.806	0.09001699	0.0498	1.547	1.747	1.809	1.865	2.105

Equations

$$B_t=e^{\lambda t},t\in [0,T],$$

$$d{S}_t=\mu S_{t}dt+\sigma (Y_t)S_{t}d{W}_t,$$

$$d{Y}_t=\frac{1}{2}\left(\frac{\kappa }{Y_t}-\theta Y_t\right)dt+\frac{\nu }{2}d{B}_t^H,t\in [0,T],$$

$$d{X}_t=(\kappa -\theta X_t)dt+\nu \sqrt{X_t}d{B}_t^H,X_0=Y_0^2>0,$$

$$\sum _{k=1}^n\frac{X_{t_k}+X_{t_{k-1}}}{2}(B_{t_k}^H-B_{t_{k-1}}^H),$$

$${\int }_0^t\sqrt{X_s}d{B}_s^H={\int }_0^t{Y}_{s}d{B}_s^H$$

$$\mathbb{E}f({\hat{S}}_T^n)\to \mathbb{E}f(S_T),n\to \infty ,$$

$$K(t,s)=c_H{s}^{\frac{1}{2}-H}{\int }_s^t{u}^{H-\frac{1}{2}}(u-s{)}^{H-\frac{3}{2}}du\text{mathbbm}{1}_{s<t},$$

$$c_H={\left(\frac{H(2H-1)}{B(2-2H,H-\frac{1}{2})}\right)}^{1/2},$$

$$B_t^H={\int }_0^{t}K(t,s)d{V}_s,t\in [0,T],$$

$$\mathbb{E}{W}_t{V}_t=\rho t,t\in [0,T],$$

$$\sigma (x)\le C_{\sigma }(1+|x{|}^q),x\in \mathbb{R};$$

$$|\sigma (x)-\sigma (y)|\le C_{\sigma }|x-y{|}^r,x,y\in \mathbb{R};$$

$${\sigma }^{\prime }(x)\le C_{\sigma }(1+|x{|}^{q^{\prime }})a.e.$$

$$\mathbb{E}\left[\exp \left\{{\int }_0^t\sigma (Y_s)d{W}_s\right\}\right]$$

$$f(x)\le C_f(1+x^p).$$

$$Y_t^p\le C_1+C_2\underset{s\in [0,T]}{\sup }|B_s^H{|}^p.$$

$$\underset{t\in [0,T]}{\sup }\mathbb{E}{Y}_t^p<\infty .$$

$$\kappa \exp \left\{\frac{\theta t}{2}\right\}\ge H(2H-1)(p+1){\int }_0^t\frac{{\nu }^2}{2}\exp \left\{\frac{\theta s}{2}\right\}|t-s{|}^{2H-2}ds.$$

$$\underset{t\in [0,T]}{\sup }\mathbb{E}\left[\frac{1}{Y_t^p}\right]<C.$$

$$p+1\le \frac{2\kappa }{{\nu }^{2}H{T}^{2H-1}}.$$

$$\mathbb{E}\left[{\int }_{t_0}^t\frac{1}{Y_u^{1+\beta }}du\right]\le \frac{4}{\kappa (1-\beta )}\mathbb{E}(Y_t^{1-\beta }-Y_{t_0}^{1-\beta })+\frac{2\theta }{\kappa }{\int }_{t_0}^t\mathbb{E}{Y}_u^{1-\beta }du<\infty .$$

$$\mathbb{E}|Y_t-Y_s{|}^{1+\beta }\le C|t-s{|}^{\beta }.$$

$$\kappa \exp \left\{\frac{\theta t}{2}\right\}\ge H(2H-1)(1+p){\int }_0^t\frac{{\nu }^2}{2}\exp \left\{\frac{\theta s}{2}\right\}|t-s{|}^{2H-2}ds,$$

$$\underset{t\in [0,T]}{\sup }\mathbb{E}\left[\frac{1}{Y_t^p}\right]<\infty .$$

$$(t-s{)}^{p-1}\mathbb{E}{\int }_s^t\frac{1}{Y_u^p}du<C(t-s{)}^p,$$

$$\mathbb{E}|Y_t-Y_s{|}^p\le C|t-s{|}^{pH}.$$

$$\mathbb{E}\exp \left\{x\underset{t\in [0,T]}{\sup }|Y_t{|}^{\varrho }\right\}<\infty .$$

$$S_t=S_0\exp \left\{\mu t+{\int }_0^t\sigma (Y_s)d{W}_s-\frac{1}{2}{\int }_0^t{\sigma }^2(Y_s)ds\right\}.$$

$$\mathbb{E}\left[\exp \left\{{\int }_0^t\sigma (Y_s)d{W}_s\right\}\right].$$

Full text

1. Introduction

Despite its undoubtedly significant historical and theoretical value, the classical Black-Scholes model does not explain numerous empirical phenomena that can be observed on real-life markets, such as implied volatility smile and skew. In order to overcome this issue, [17] and, later, [15] introduced stochastic volatility models that emerged into an essential subject of research activity in financial modeling nowadays.

To illustrate the range of existing models (without trying to list all possible references), we recall the approaches of [1], [4], [5], [8], [11], [18], [20], [29], and so on.

A separate class of stochastic volatility models are those based on fractional Brownian motion. They allow to reflect the so-called “memory phenomenon” of the market (for more detail on market models with memory see, for instance, [3, 12, 31]). In this context, we should also mention [7, 9, 10] and [6].

In the present paper, we consider option pricing in a framework of the fractional modification of the Heston-type model, namely a financial market with a finite maturity time $T$ that is composed of two assets:

(i) a risk-free bond (or bank account) $B=\{B_{t},~{}t\in[0,T]\}$, the dynamics of which is characterized by the formula

[TABLE]

where $\lambda\in\mathbb{R}^{+}$ represents the risk-free interest rate;

(ii) a risky asset $S=\{S_{t},~{}t\in[0,T]\}$, the evolution in time of which is given by the system of stochastic differential equations

[TABLE]

[TABLE]

with non-random initial values $S_{0},Y_{0}>0$, where the process $W=\{W_{t},~{}t\geq 0\}$ is a standard Wiener process, $\mu\in\mathbb{R}$, $\kappa,\theta,\nu>0$ are constants, $\sigma$: $[0,\infty)\to[0,\infty)$ is a function that satisfies some regularity properties and $B^{H}=\{B^{H}_{t},~{}t\in[0,T]\}$ is a fractional Brownian motion with the Hurst index $\frac{1}{2}<H<1$, which corresponds to the “long memory” case. $W$ and $B^{H}$ are assumed to be correlated.

The process $Y$ was extensively studied in [26, 27] and, for the case $\kappa=0$, in [25]. Note that, according to [28], the process $Y$ exists, is unique and has continuous paths until the first moment of zero hitting. Moreover, in Theorem 2 of [26] it was shown that in case of $\kappa>0$ and $H>\frac{1}{2}$ such process is strictly positive and never hits zero, therefore exists, is unique and continuous on the entire $[0,T]$.

Such choice of the volatility process can be explained by the fact that $Y$ can be interpreted as the square root of the fractional version of Cox-Ingersoll-Ross process. Indeed, according to [26], Theorem 1, the process $X=\{Y^{2}(t),~{}t\in[0,T]\}$ satisfies the stochastic differential equation of the form

[TABLE]

until the first moment of zero hitting, where the integral $\int_{0}^{t}\sqrt{X_{s}}dB_{s}^{H}$ is considered as the pathwise limit of the sums

[TABLE]

as the mesh of the partition $0=t_{0}<t_{1}<...<t_{n}=t$ tends to zero.

Note that, due to Kolmogorov theorem, fractional Brownian motion $B^{H}$ has a modification with Hölder continuous paths up to order $H$. Hence, from the form of the equation (3), the process $Y$ also has a modification with trajectories that are Hölder-continuous up to order $H$. Therefore, in case of $H>\frac{1}{2}$, the sum of Hölder exponents of the integrator and integrand in the integral

[TABLE]

exceeds 1 and, due to [32], the corresponding integral exists as the pathwise limit of Riemann-Stieltjes integral sums.

It should be also mentioned that for the case $H<1/2$, the process $Y$ can hit zero and it is not clear whether the solution exists on the entire $[0,T]$ (see [27] for more detail). Therefore, we will concentrate on the case $H>1/2$. For more information on markets with rough volatility see, for example, [14] or [19].

An analogue of the model (2), (3) was considered in [6] with fractional Ornstein-Uhlenbeck process instead of $Y$. However, Ornstein-Uhlenbeck process can take negative values with positive probability which is a notable drawback for a stochastic volatility model.

Note that it is impossible to calculate $\mathbb{E}f(S_{T})$ (with $f$ being a payoff function) for option pricing analytically, so numerical methods should be used. Therefore it is required to provide a decent discretization scheme for $S_{T}$ and prove the convergence

[TABLE]

where $\hat{S}^{n}$ is a discretized version of the process $S$. Moreover, we allow $f$ to have discontinuities of the first kind which can cause errors in straightforward Monte-Carlo estimation of the expectation, so we provide an alternative formula with smooth functional under the expectation. In such framework, we also give the rate of convergence (4).

It should be mentioned that the market with risky asset defined by (2)–(3) is arbitrage-free, incomplete but admits minimal martingale measure (see Section 3). However, the expectations calculated with respect to the minimal martingale and objective measures differ only by non-random coefficient, therefore, for simplicity, we concentrate on expectation with respect to the objective measure. In order to model the volatility $Y$, we use the inverse Euler approximation scheme studied in [16].

The paper is organized as follows. In Section 2, we describe main assumptions concerning relation between the Wiener process and the fractional Brownian motion as well as volatility function $\sigma$ and payoff function $f$. In Section 3 several important properties of both price and volatility processes are presented and the arbitrage-free property is discussed. In Section 4 we apply the Malliavin calculus techniques, following [1] and [6], to obtain the formula for option price that does not contain discontinuities (which are allowed for the payoff function $f$). In Section 5, we study the rate of convergence of Monte-Carlo estimation of the option price $\mathbb{E}f(S_{T})$ based on inverse Euler approximation scheme for fractional CIR process presented in [16]. In Section 6, we give results of numerical simulations for different payoff functions $f$. Section 7 contains the proofs of all results of the paper. Appendix A is devoted to several well-known results from the Malliavin calculus used in this paper.

2. Model description and main assumptions

Consider the market with risk-free asset $B$ given by (1) and risky asset $S$, the dynamics of which is described by stochastic differential equations (2), (3).

Denote

[TABLE]

[TABLE]

where $B(\cdot,\cdot)$ is the Beta function. Then, according to [21], the process $B^{H}=\{B^{H}_{t},~{}t\in[0,T]\}$ given by

[TABLE]

where $V=\{V_{t},~{}t\in[0,T]\}$ is a Wiener process, is the fractional Brownian motion with Hurst parameter $H$.

The processes $W$ and $B^{H}$ from (2), (3) are assumed to be correlated and the form of the dependence is defined on the basis of representation (5) as follows.

Assumption 1**.**

The processes $W$ and $V$ from (2) and (5) correspondingly are correlated:

[TABLE]

with some constant $\rho\in[-1,1]$.

Remark 2.1**.**

Assumption 1 means that $W_{t}=\rho V_{t}+\sqrt{1-\rho^{2}}\tilde{V}_{t}$, $t\in[0,T]$, where $\tilde{V}$ is a Wiener process independent of $V$.

The function $\sigma$: $\mathbb{R}\to\mathbb{R}$ is assumed to satisfy the following conditions.

Assumption 2**.**

For some constant $C_{\sigma}>0$:

• (i)

there exists such $\sigma_{\min}>0$ that for all $x\in\mathbb{R}$: $\sigma(x)>\sigma_{\min}>0$;

• (ii)

$\sigma$* has moderate polyniomial growth, i.e. there is such $q\in(0,1)$ that*

[TABLE]

• (iii)

$\sigma$* is uniformly Hölder continuous, i.e. there is such $r\in(0,1]$ that*

[TABLE]

• (iv)

$\sigma$* is differentiable a.e. w.r.t. the Lebesgue measure on $\mathbb{R}$ and there exists such $q^{\prime}>0$ that*

[TABLE]

Remark 2.2**.**

1) Item (i) in Assumption 2 is required for theoretical calculations as we will divide on $\sigma$ in what follows.

2) Item (ii) is necessary to ensure the finiteness of expectations of the form

[TABLE]

in case if the Wiener process $W$ and the fractional Brownian motion $B^{H}$ from (2) and (3) are correlated (see Remark 3.3 for discussion). Note that in standard Heston model moment explosions may appear as well, see e.g. [2].

3) (ii) follows from (iii) in the case $r<1$, while in (iii) we also allow $r=1$.

In the framework above, we consider an option with a measurable payoff function $f:~{}\mathbb{R}^{+}\to\mathbb{R}^{+}$ depending on the value $S_{T}$ of the stock at maturity time $T$ which satisfies the following properties:

Assumption 3**.**

For some constant $C_{\sigma}>0$:

• (i)

$f$* is of polynomial growth, i.e. there are such $C_{f}>0$ and $p>0$ that*

[TABLE]

• (ii)

$f$* is locally Riemann integrable, possibly, having discontinuities of the first kind.*

Remark 2.3**.**

In what follows, we will denote $C$ any positive constant that does not depend on time variable or diameter of the partition and the exact value of which is not important. Note that $C$ may change from line to line (and even within one line).

3. Model properties

3.1. Properties of stochastic volatility process

In what follows we will require an auxiliary result, presented in Corollary 2.2 of [27].

Theorem 3.1**.**

For all $H\in(0,1)$, $T>0$ and $p>0$ there are such non-random constants $C_{1}=C_{1}(T,p,Y_{0},\kappa,\theta)>0$ and $C_{2}=(T,p,\theta,\sigma)>0$ that for all $t\in[0,T]$:

[TABLE]

Furthermore,

[TABLE]

The next result is crucial for obtaining discrete approximation scheme for the process $Y$ and was presented in [16].

Theorem 3.2**.**

Let $p>0$ and $\kappa$, $\theta$, $\nu$ and $T$ are such that for all $t\in[0,T]$:

[TABLE]

Then there is such constant $C=C(T,Y_{0},\theta)$ that

[TABLE]

Remark 3.1**.**

Condition (6) is satisfied if, for example,

[TABLE]

See Remarks 3.1 and 3.2 in [16] for discussion.

Note that condition (6) involves $T$ and does not guarantee the existence of the inverse moments on whole $\mathbb{R}_{+}$. However, the following result concerning the integrated inverse moments of the volatility process $Y$ holds true.

Theorem 3.3**.**

Let $\beta\in\left(0,\min\{1,\frac{\kappa}{\nu^{2}HT^{2H-1}}\}\right)$. Then, for all $0\leq t_{0}<t\leq T$:

[TABLE]

Theorem 3.4**.**

Let $\beta\in\left(0,\min\{1,\frac{\kappa}{\nu^{2}HT^{2H-1}}\}\right)$. Then, there is such $C=C(\kappa,\theta,\nu,T,\beta)>0$ that for any $0\leq s<t\leq T$:

[TABLE]

Remark 3.2**.**

Let $p>1$ and for all $t\in[0,T]$:

[TABLE]

i.e., due to Theorem 3.2,

[TABLE]

Proceeding just as in proof of Theorem 3.4 and taking into account that

[TABLE]

we can easily obtain that

[TABLE]

3.2. Properties of the price process

Now let us consider several properties of the price process $S$ defined by the stochastic differential equation (2).

Theorem 3.5**.**

For any $x>0$ and $\varrho\in[0,2)$:

[TABLE]

• 2.

Equation (2) has a unique solution of the form

[TABLE]

Remark 3.3**.**

As it was mentioned in Remark 2.2, presence of function $\sigma$ in (2), the choice of which is restricted by Assumption 2, is required to ensure finiteness of the moments of the form

[TABLE]

Note that Assumption 2, (i) and (ii), does not allow $\sigma$ to be linear function, i.e. we do not consider straigthforward modification of the Heston model of the form

[TABLE]

[TABLE]

where $\mu\in\mathbb{R}$, $\kappa$, $\theta$, $\nu$, $\sigma>0$ are constants.

However, in case of independent $W$ and $B^{H}$, i.e. when $\rho=0$ in Assumption 1, it is easy to see (e.g. by conditioning on $Y$ and solving the conditioned equation) that equation (9) has a unique solution of the form

[TABLE]

Moreover, $\mathbb{E}S_{t}<\infty$ for all $t\in[0,T]$, because the process $\tilde{S}$, such that

[TABLE]

is a non-negative local martingale and, therefore, a supermartingale.

3.3. Arbitrage-free property and incompleteness

For the market (1)–(3), we can obtain the following result which is similar to the one in [6], Theorem 4.

Theorem 3.6**.**

Let the function $\sigma$ satisfy Assumption 2. Then the market (1)–(3) has the following properties.

• $(i)$

It is arbitrage-free and incomplete.

• $(ii)$

Any probability measure $\mathbb{Q}$ such that

[TABLE]

where $\eta_{i}$, $i=1,2$, are non-anticipative, bounded and satisfy the condition

[TABLE]

is a martingale measure.

• $(iii)$

Taking $\eta_{1}=\rho\frac{\lambda-\mu}{\sigma(Y_{s})}$ and $\eta_{2}=\sqrt{1-\rho^{2}}\frac{\lambda-\mu}{\sigma(Y_{s})}$, we get the minimal martingale measure.

4. Option pricing in fractional Heston model

In this section, we will use the tools of Malliavin calculus to obtain the formula that can be used for computation of

[TABLE]

Consider two-dimensional Wiener process $(V,\tilde{V})$, where $V$ is given in Volterra representation (5) and $\tilde{V}$ is defined in Remark 2.1. Denote $(D^{V},D^{\tilde{V}})$ the stochastic derivative with respect to the two-dimensional Wiener process $(V,\tilde{V})$ and recall $K$ is the kernel from representation (5). Denote also

[TABLE]

Lemma 4.1**.**

• (i)

The stochastic derivatives of the fBm $B^{H}$ are equal to

[TABLE]

• (ii)

The stochastic derivatives of the volatility process $Y$ are

[TABLE]

where $h(s):=\frac{1}{2}\left(\frac{\kappa}{Y_{s}^{2}}+\theta\right)$.

• (iii)

The stochastic derivatives of $X$ are equal to

[TABLE]

Denote

[TABLE]

and consider a random variable

[TABLE]

Note that, due to Assumption 2, (i), $Z_{T}$ is correctly defined.

Theorem 4.1**.**

Under Assumptions 2 and 3, the option price $\mathbb{E}f(S_{T})=\mathbb{E}g(X_{T})$ can be represented as

[TABLE]

or, alternatively,

[TABLE]

5. Inverse Euler approximation scheme for the volatility and price processes

Let $0=t^{n}_{0}<t^{n}_{1}<...<t^{n}_{n}=T$ be an equidistant partition of the interval $[0,T]$, $t^{n}_{i}=\frac{iT}{n}$, $\Delta_{n}:=\frac{1}{n}$, $\Delta B^{H}_{k+1}:=B^{H}_{t^{n}_{k+1}}-B^{H}_{t^{n}_{k}}$ and consider the approximation scheme of the form

[TABLE]

with linear interpolation between the points of the partition.

Note that approximations given by (15) are strictly positive and it is easy to verify that in points of partition they satisfy the following difference equation:

[TABLE]

Approximations of the form (15) were presented and studied in [16]. We give the result concerning the convergence rate of these approximations (for more detail, see Theorem 4.2 in [16]).

Theorem 5.1**.**

Let $\xi\in(0,1)$, $p\geq 2$, $\Delta_{n}<1-\xi$ and parameters $\theta,\kappa,\nu>0$ are such that for all $t\in[0,T]$:

[TABLE]

Then there is such $C=C(T,H,p,Y_{0},\theta,\kappa,\nu,\xi)>0$ that

[TABLE]

Remark 5.1**.**

Condition (17) is a sufficient condition for finiteness of the inverse moments of $Y$ of order $3p$, namely for

[TABLE]

Three approximations of the volatility process $Y$ trajectories given by the formula (15) with $T=1$, $\kappa=1$, $\theta=1$, $\nu=0.14$, $Y_{0}=1$, $H=0.7$ and $\Delta_{n}=0.0001$ are presented on Fig. 1.

For the sake of simplicity, instead of linear interpolation between the points of the partition, we put $\hat{Y}^{n}_{t}=\hat{Y}^{n}_{t^{n}_{k}}$ for $t\in[\hat{Y}^{n}_{t^{n}_{k}},\hat{Y}^{n}_{t^{n}_{k+1}})$. It should be noted that in this case speed of convergence of approximations remains the same as in Theorem 5.1 due to Remark 3.2 because

[TABLE]

Denote

[TABLE]

where $X_{0}:=\log S_{0}$, and consider the discretized process

[TABLE]

where $\Delta W_{j}=W_{t_{j+1}^{n}}-W_{t_{j}^{n}}$.

Before going to the main theorem of the paper, let us prove several auxiliary results.

Theorem 5.2**.**

Let $p\geq 1$. Then, for all $H\in(0,1)$:

[TABLE]

Remark 5.2**.**

Note that approximations (15) (see Fig. 2) are correctly defined for $H<1/2$ and Theorem 5.2 holds for an arbitrary Hurst parameter as well. However, for $H<1/2$ behaviour of $\hat{Y}^{n}$ as $n\to\infty$ remains obscure.

Corollary 5.1**.**

Approximating processes $\hat{Y}^{n}$ have bounded exponential moments, i.e. for any $x>0$ and $\varrho<2$:

[TABLE]

Remark 5.3**.**

From Theorem 5.2, Corollary 5.1 and Assumption 2 (ii), using the same argument as in the proof of Theorem 3.5, it is easy to verify that for any $m\in\mathbb{Z}$:

[TABLE]

Theorem 5.3**.**

Let $n\geq 2$ and conditions of Theorem 5.1 hold for $p=4$. Then, under Assumption 2, there exists a constant $C$ such that

[TABLE]

[TABLE]

where $\hat{Z}^{n}_{T}:=\int_{0}^{T}\sigma^{-1}(\hat{Y}^{n}_{u})d\tilde{V}_{u}$.

Lemma 5.1**.**

Let $n\geq 2$ and conditions of Theorem 5.1 hold for $p=32$. Then, under Assumptions 2 and 3, there is such $C_{F}>0$ that

[TABLE]

Theorem 5.4**.**

Let $n\geq 2$ and conditions of Theorem 5.1 for $p=32$ hold. Then, under Assumptions 2 and 3,

[TABLE]

6. Simulations

In this section, we use the discretization scheme studied previously to estimate option price for several payoff functions $f$. In all simulations we use $T=1$, $\kappa=1$, $\theta=1$ and $\nu=0.14$ to make sure that for all $H\in(1/2,1)$ the following condition is satisfied for $p=32$:

[TABLE]

which is sufficient for Theorem 5.4 to hold true. For simplicity, we also consider everywhere the case $\mu=0.5$, $\rho=0$ and $\sigma=0.5\left(x+0.01\right)^{0.9}$.

In Tables 1–3 we present descriptive statistics of Monte-Carlo estimations of $\mathbb{E}\left[\frac{F(\hat{S}^{n}_{T})}{\hat{S}^{n}_{T}}\left(1+\frac{\hat{Z}^{n}_{T}}{T}\right)\right]$ (and, therefore, $\mathbb{E}f(S_{T})$) for different functions $f$ and different partition sizes $\Delta_{n}$. On Fig. 3, (a)–(c), the data is visualized in a form of box-and-whisker plots. In each case, 1000 Monte-Carlo estimates of option price, calculated from samples of 1000 trials each, were analyzed. All calculations were performed in R using package somebm to generate trajectories of Wiener process and fractional Brownian motion.

As we can see, simulations show relatively small coefficient of variation in all cases. Nota that increasing partition size does not lead to any significant changes in standard deviation of the estimates.

7. Proofs

Proof of Theorem 3.3. Denote $\alpha:=1-\beta$ and let $\varepsilon>0$ be fixed. By applying the chain rule, we obtain:

[TABLE]

It is clear from (3) that the process $Y=\{Y_{t},~{}t\in[0,T]\}$ has trajectories that are $\delta$-Hölder-continuous for any $\delta\in(0,H)$, so the process

[TABLE]

also has Hölder-continuous trajectories up to the order $H$. Therefore, the sum of Hölder exponents of the integrator and integrand in the integral w.r.t. fractional Brownian motion in (20) exceeds 1. In this case this integral is the pathwise limit of Riemann-Stieltjes integral sums (see, for example, [32]), coincides with the pathwise Stratonovich integral and, by applying Theorem A.1, we can rewrite (20) as follows:

[TABLE]

where $D^{H}_{s}$ is the Malliavin derivative operator w.r.t. $B^{H}$ and $\int_{t_{0}}^{t}\frac{\nu\alpha}{2(Y(s)+\varepsilon)^{1-\alpha}}\delta B_{s}^{H}$ is the corresponding Skorokhod integral.

Note that

[TABLE]

From this, it is easy to verify that

[TABLE]

so

[TABLE]

Taking into account (21) and (22), we can rewrite (20) in the following form:

[TABLE]

where $\varphi(u,s):=H(2H-1)|u-s|^{2H-2}$.

Note that

[TABLE]

It is easy to verify that

[TABLE]

so

[TABLE]

Hence, if $\alpha\in\left(\max\{0,1-\frac{\kappa}{\nu^{2}HT^{2H-1}}\},1\right)$, i.e. when $0\geq-(1-\alpha)\frac{\nu^{2}}{2}HT^{2H-1}\geq-\frac{\kappa}{2}$,

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

Therefore, taking into account upper bounds (24), (25) and (26), it is obvious from (23) that

[TABLE]

or

[TABLE]

Since the expectation of the Skorokhod integral is zero, by letting $\varepsilon\to 0$ we obtain that

[TABLE]

Finiteness of the right-hand side of (27) follows from Theorem 3.1. $\square$

Proof of Theorem 3.4. From (3), Hölder’s and Jensen’s inequalities it is clear that

[TABLE]

where

[TABLE]

Note that form of $C_{3}$ follows from the fact that $B^{H}_{t}-B^{H}_{s}\sim\mathcal{N}(0,|t-s|^{2H})$ (see, for example, [30])

From Theorem 3.1 it is obvious that

[TABLE]

Finally, from Theorem 3.3,

[TABLE]

where $C_{1}=\tilde{C}_{1}\mathbb{E}\int_{s}^{t}\frac{1}{Y_{u}^{1+\beta}}du<C_{1}(t-s)^{\beta}$.

The statement of the Theorem now follows from (28), (29) and (30) as well as the fact that from condition $\beta<1$ it is easy to verify that for any $H\in(1/2,1)$:

[TABLE]

$\square$

Proof of Theorem 3.5. 1. From Theorem 3.1, for all $\varrho\in[0,2)$:

[TABLE]

and, due to [13], for all $x>0$ and $\varrho\in[0,2)$:

[TABLE]

Hence,

[TABLE]

2. In order to show that the representation (8) indeed holds, it is sufficient to prove that the integrals $\int_{0}^{t}\sigma(Y_{s})dW_{s}$ and $\int_{0}^{t}\sigma(Y_{s})S_{s}dW_{s}$ are well-defined, while the form of the representation can be obtained straightforwardly.

Note that (see, for example, [24]) for all $p>0$

[TABLE]

so, due to item $(ii)$ from Assumption 2 and Theorem 3.1,

[TABLE]

and the integral $\int_{0}^{t}\sigma(Y_{s})dW_{s}$ is well-defined.

Now consider the integral $\int_{0}^{t}\sigma(Y_{s})S_{s}dW_{s}$. As

[TABLE]

it is sufficient to check two conditions:

[TABLE]

Using Theorem 3.1 and Assumption 2, $(ii)$, it is easy to verify that

[TABLE]

Moreover, from (7), for any $x>0$:

[TABLE]

hence, for all $n\in\mathbb{Z}$, by putting $x:=\frac{4n^{2}}{2}$, we obtain the Novikov’s condition for the process $-2n\sigma(Y_{t})$, $t\in~{}[0,T]$.

Consequently,

[TABLE]

and so

[TABLE]

due to (33).

Therefore, from (31), (32) and (35),

[TABLE]

and so the integral $\int_{0}^{t}\sigma(Y_{s})S_{s}dW_{s}$ is well-defined. $\square$

Proof of Theorem 3.6. The proof is similar to the proof of Theorem 4 in [6].

Proof of Lemma 4.1. Item $(i)$ can be found in [6]. In particular, $D^{\tilde{V}}_{u}Y_{t}=0$ in $(ii)$ follows from independence of $Y$ and $V$.

Applying stochastic derivative operator to both parts of the integral form of (3), we get

[TABLE]

Application of the chain rule with the function $F(x)=1/x$ can be justified by the same argument as in Remark 10 of [6], since $F$ is locally Lipschitz on $(0,\infty)$.

According to [26], Theorem 2, $Y$ does not hit zero a.s. Therefore $h$ is well defined a.s., and (36) means that for a fixed $u$, the process $\{Z_{t},t\in[0,t]\}$ defined by $Z_{t}:=D^{B}_{u}Y_{t}$ satisfies a random linear integral equation of the form

[TABLE]

This is a Volterra equation, and its solution is given by

[TABLE]

Note that $K$ is differentiable in the first argument ($\frac{\partial}{dt}K(t,s)$ is well defined for $t>s$), so (38) can be checked by substituting in (37) and taking derivatives of both sides.

Both derivatives in $(iii)$ are obtained by direct differentiation following the Malliavin derivative rules, see e.g. [23], Proposition 3.4. Since $Y$ is independent of $\tilde{V}$,

[TABLE]

To find $D_{u}^{V}X_{t}$, we note that

[TABLE]

$\square$

Proof of Theorem 4.1. The result can be obtained by following the proof of Lemma 11 in [6], taking into account Lemma 4.1 and relation (35). $\square$

Proof of Theorem 5.2. First, note that for any fixed $n$ and $k=0,1,...,n$:

[TABLE]

By continuing calculations above recurrently and taking into account that $\hat{Y}^{n}_{t^{n}_{0}}=Y_{0}$, it is easy to see that there is such constant $C_{n}$ that

[TABLE]

Moreover, for any fixed $N$ there is such constant $C_{N}$ that

[TABLE]

Let us prove that there is such $C>0$ (which does not depend on $n$) that

[TABLE]

From calculations above, it will be enough to show that, for some $N\geq 1$,

[TABLE]

Let $n>2(8\theta)^{p}T^{p-1}$ be fixed. Consider the last moment of staying above level $Y_{0}/2$, i.e.

[TABLE]

Let us prove that for any point of the partition $t^{n}_{k}$, $k=1,...,n$, the following inequality holds:

[TABLE]

In order to do that, we will separately consider cases $t^{n}_{k}\leq t^{n}_{\tau_{1}}$ and $t^{n}_{k}>t^{n}_{\tau_{1}}$.

Step 1. Assume that $t^{n}_{k}\leq t^{n}_{\tau_{1}}$. Then, due to representation (16),

[TABLE]

Note that for all $t^{n}_{k}\leq t^{n}_{\tau_{1}}$:

[TABLE]

Moreover, from Jensen’s inequality,

[TABLE]

Finally,

[TABLE]

Hence, for all $t^{n}_{k}\leq t^{n}_{\tau_{1}}$:

[TABLE]

Step 2. Assume that $\tau_{1}\neq n$, i.e. there are points of partition on the interval $(t^{n}_{\tau_{1}},T]$. From definition of $\tau_{1}$, $\hat{Y}^{n}_{t^{n}_{\tau_{1}}}\geq\frac{Y_{0}}{2}$ and for all points of the partition $t^{n}_{k}$ such that $t^{n}_{k}\in(t^{n}_{\tau_{1}},T]$:

[TABLE]

Let $t^{n}_{k}\in(t^{n}_{\tau_{1}},T]$ be fixed and denote

[TABLE]

It is obvious that $t^{n}_{\tau_{1}}<t^{n}_{\tau_{2}^{k}}\leq t^{n}_{k}$ and $\hat{Y}^{n}_{t^{n}_{\tau_{2}^{k}}}<\frac{Y_{0}}{2}$, and

[TABLE]

In addition, if $t^{n}_{\tau_{2}^{k}}=t^{n}_{k}$,

[TABLE]

and if $t^{n}_{\tau_{2}^{k}}<t^{n}_{k}$,

[TABLE]

From definition of $\tau_{2}^{k}$, for all points of the partition $t^{n}_{l}\in(t^{n}_{\tau_{2}^{k}},t^{n}_{k}]$ it holds that $\hat{Y}^{n}_{t^{n}_{k}}\geq\frac{Y_{0}}{2}$, so

[TABLE]

Furthermore,

[TABLE]

and

[TABLE]

Hence,

[TABLE]

Finally, from (41) and (42),

[TABLE]

Therefore, (40) indeed holds for any point $t^{n}_{k}$ of the partition.

Step 3. As $n>2(8\theta)^{p}T^{p-1}$,

[TABLE]

therefore, as, due to (40),

[TABLE]

we have

[TABLE]

Using the discrete version of the Grönwall’s lemma, we obtain:

[TABLE]

i.e., taking into account that the right-hand side does not depend on $n$ and remarks in the beginning of the proof, there is such $C>0$ that

[TABLE]

Now the claim of the Theorem follows from the fact that the right-hand side of (43) does not depend on $n$ and that (see, for example, [24])

[TABLE]

$\square$

Proof of Corollary 5.1. From (43) it follows that there is such $C>0$ that

[TABLE]

The rest of the proof is similar to Theorem 3.5, 1. $\square$

Proof of Theorem 5.3. We shall proceed as in proof of Lemma 14, [6].

Using Hölder’s inequality, we write:

[TABLE]

From Assumption 2 (iii), Jensen’s inequality and Theorem 5.1,

[TABLE]

Moreover, Assumption 2, (ii) and (iii), implies that

[TABLE]

From Theorem 5.1,

[TABLE]

and, from Theorems 3.1 and 5.2,

[TABLE]

Therefore, taking into account bounds above, there is such constant $C>0$ that

[TABLE]

Now, let us prove (19). Taking into account Assumption 2 (i),

[TABLE]

so, from Assumption 2 (iii),

[TABLE]

$\square$

Proof of Lemma 5.1. It is clear that

[TABLE]

Now we shall estimate the right-hand side of (44) term by term.

[TABLE]

From Assumption 3 (i), both $f$ and $F$ are of polynomial growth, therefore, due to (35),

[TABLE]

Furthermore, using sequentially the inequalities

[TABLE]

and Hölder’s inequality, we obtain that

[TABLE]

Next, from (34) and Remark 5.3 it follows that

[TABLE]

so, using this together with Hölder and Burkholder-Davis-Gundy inequalities, we continue the chain as follows:

[TABLE]

By applying Assumption 2, (ii) and (iii),

[TABLE]

and

[TABLE]

From Theorems 3.1 and 5.2,

[TABLE]

and, according from Theorem 5.1,

[TABLE]

hence

[TABLE]

Now, let us move to the second term of the right-hand side of (44).

[TABLE]

Due to Remark 5.3,

[TABLE]

and, from Assumption 3 (i),

[TABLE]

According to (35) and Remark 5.3,

[TABLE]

so

[TABLE]

To get the final result, we can proceed just as in the upper bound for the first term in the right-hand side of (44). Thus

[TABLE]

Relations (45) and (46) together with (44) complete the proof. $\square$

Proof of Theorem 5.4. According to Theorem 4.1,

[TABLE]

According to Theorem 3.5, Assumption 3 (i) and the Cauchy-Schwartz inequality, $\mathbb{E}\left(\frac{F(S_{T})}{S_{T}}\right)^{2}<\infty$. Next,

[TABLE]

The proof now follows from Theorem 5.3 and Lemma 5.1. $\square$

Appendix A Necessary results from Malliavin Calculus

In this section, we recall several main definitions and results related to Malliavin calculus. For more detail, we refer to [22].

Let $B^{H}=\{B^{H}_{t},~{}t\in[0,T]\}$ be a fractional Brownian motion with $H\in[1/2,1)$ on the standard probability space $\{\Omega,\mathcal{F},\mathbb{F}=\{\mathcal{F}_{t}\}_{t\in[0,T]},\mathbb{P}\}$, where $\Omega=C([0,T],\mathbb{R})$, i.e. a centered Gaussian process that starts in zero and has a covariance function of the form

[TABLE]

Note that the covariance function of the fractional Brownian motion has the form

[TABLE]

where $\varphi(\tau,u):=H(2H-1)|u-\tau|^{2H-2}$.

On the set of all step functions on $[0,T]$, define an inner product that acts as follows for the indicator functions:

[TABLE]

Denote $\mathcal{H}$ the Hilbert space that is the closure of the space of all step functions on $[0,T]$ with respect to $\langle\cdot,\cdot\rangle_{\mathcal{H}}$.

Remark A.1**.**

If $H=1/2$, $\mathcal{H}$ coincides with $L_{2}([0,T])$.

The mapping $\mathbbm{1}_{[0,t]}\to B^{H}_{t}$ can be extended to a linear isometry from $\mathcal{H}$ onto a closed subspace $\mathcal{H}_{1}$ of $L^{2}(\Omega,\mathcal{F},\mathcal{P})$ associated with $B^{H}$. We will denote this isometry by $\phi\to B^{H}_{\phi}$. In this case, for all $\phi,\psi\in\mathcal{H}$:

[TABLE]

Denote by $C^{\infty}_{p}(\mathbb{R}^{n})$ the set of all infinitely differentiable functions with the derivatives of at most polynomial growth at infinity.

Definition A.1**.**

Random variables $\xi$ of the form

[TABLE]

where $h\in C^{\infty}_{p}(\mathbb{R}^{n})$, $\phi_{1},...,\phi_{n}\in\mathcal{H}$, $n\geq 1$, are called smooth.

Denote $\mathcal{S}$ the set of all smooth random variables.

Definition A.2**.**

Let $\xi\in\mathcal{S}$. The stochastic or Malliavin derivative of a smooth random variable $\xi$ of the form (47) is the $\mathcal{H}$-valued random variable given by

[TABLE]

Remark A.2**.**

If $\phi_{i}=\mathbbm{1}_{[0,t_{i}]}$, $t_{i}\in[0,T]$, $i=1,...,n$, then $B^{H}_{\mathbbm{1}_{[0,t_{i}]}}=B^{H}_{t_{i}}$ and the real-valued random variable of the form

[TABLE]

is called the stochastic derivative of $\xi$ at time $t$.

According to Proposition 1.2.1 from [22], $D$ as an operator from the subset of $L^{p}(\Omega)$ to $L^{p}(\Omega,\mathcal{H})$ is closable for any $p\geq 1$ and we shall use the same notation $D$ for the closure.

Definition A.3**.**

Let $p\geq 1$. The domain $\mathbb{D}^{1,p}$ of $D$ is the closure of the class of smooth random variables $\mathcal{S}$ with respect to the norm

[TABLE]

Remark A.3**.**

For $p=2$, the space $\mathbb{D}^{1,2}$ is the Hilbert space with respect to the inner product

[TABLE]

Proposition A.1**.**

([22], Proposition 1.2.3) Let $F$: $\mathbb{R}^{m}\to\mathbb{R}$ be a continuously differentiable function with bounded partial derivatives, and fix $p\geq 1$. Suppose that $\xi=(\xi_{1},...,\xi_{m})$ is a random vector whose components belong to the space $\mathbb{D}^{1,p}$. Then $F(\xi)\in\mathbb{D}^{1,p}$ and

[TABLE]

Remark A.4**.**

In what follows, we will consider the case $p=2$.

Definition A.4**.**

The divergence or Skorokhod operator $\delta$ is the adjoint of the operator $D$, i.e. an undounded operator on $L^{2}(\Omega,\mathcal{H})$ with values in $L^{2}(\Omega)$ such that:

• (i)

the domain of $\delta$, denoted by $\mathsf{Dom~{}}\delta$, is the set of $\mathcal{H}$-valued square integrable random variables $\zeta\in L^{2}(\Omega,\mathcal{H})$ such that for all $\xi\in\mathbb{D}^{1,2}$:

[TABLE]

where $C_{\zeta}$ is some constant depending on $\zeta$;

• (ii)

if $\zeta$ belongs to $\mathsf{Dom~{}}\delta$, then $\delta(\zeta)$ is the element of $L^{2}(\Omega)$ characterized by

[TABLE]

for any $\xi\in\mathbb{D}^{1,2}$.

The Skorokhod operator $\delta$ is closed.

Remark A.5**.**

Let $V=\{V_{t},~{}t\in[0,T]\}$ be the Wiener process, $\mathcal{H}_{V}=L_{2}([0,T])$ be the associated Hilbert space (see Remark A.1) and $\delta_{V}$ be the corresponding divergence operator. In this case, the elements of $\mathsf{Dom}~{}\delta_{V}\subset L^{2}([0,T]\times\Omega)$ are square-integrable processes, and the divergence $\delta_{V}(\zeta)$ is called the Skorokhod stochastic integral of the process $\zeta$ with respect to $V$ and is denoted as follows:

[TABLE]

According to [22], Section 1.3.2, the Skorokhod integral is correctly defined for all elements of the space $\mathbb{L}^{1,2}=L^{2}([0,T],\mathbb{D}^{1,2})$ with the norm $\lVert\cdot\rVert_{\mathbb{L}^{1,2}}$ such that

[TABLE]

Remark A.6**.**

Let $B^{H}$ be a fractional Brownian motion with $H>1/2$. Similarly to the Wiener process case, we shall call the corresponding divergence $\delta_{H}(\zeta)$ the Skorokhod stochastic integral with respect to fractional Brownian motion and shall denote it as

[TABLE]

In what follows, we shall use the definition of pathwise stochastic integral with respect to fractional Brownian motion proposed in [32] and denote it by $\int_{0}^{T}\zeta_{t}dB^{H}_{t}$. There is a useful result that connects stochastic and Skorokhod integrals, which is given below.

Let $H>1/2$ and

[TABLE]

Theorem A.1** ([22], Proposition 5.2.1).**

Let $\zeta=\{\zeta_{t},~{}t\in[0,T]\}$ be a stochastic process in the space $\mathbb{D}^{1,2}(|\mathcal{H}|)$ with Hölder continuous trajectories up to the order $H$ and $D_{s}^{H}$ be the Malliavin derivative operator with respect to $B^{H}$. Suppose that a.s.

[TABLE]

Then $\zeta$ is Stratonovich integrable and

[TABLE]

Acknowledgments

The first author acknowledges that the present research is carried through within the frame and support of the ToppForsk project nr. 274410 of the Research Council of Norway with title STORM: Stochastics for Time-Space Risk Models. The second author was partially supported by the grant 346300 for IMPAN from the Simons Foundation and the matching 2015-2019 Polish MNiSW fund.