APPROXIMATE PRICING OF DERIVATIVES UNDER FRACTIONAL STOCHASTIC VOLATILITY MODEL

Abstract

This paper examines the issue of derivative pricing within the framework of a fractional stochastic volatility model. We present a deterministic partial differential equation system to derive an approximate expression for the derivative price. The proposed approach allows for the stochastic volatility to be expressed as a composition of deterministic functions of time and a fractional Ornstein–Uhlenbeck process. We apply this method to the European option pricing under the fractional Stein–Stein volatility model, demonstrating its feasibility and reliability through numerical simulations. Our numerical simulations also illustrate the impact of the parameters in the fractional stochastic volatility model on the option price.

1. Introduction

The Black–Scholes–Merton pricing formula has long been regarded as a fundamental tool for derivatives analysis [5, 20]. However, with the continuous fluctuation of volatilities observed in financial markets, the stochastic volatility models have been shown to provide a better description of market behaviour than the model with constant volatility. Stein and Stein [25] examined the stock price distribution under a diffusion process with a stochastic volatility parameter, known as the Stein–Stein model. Chernov et al. [7] evaluated the effectiveness of various volatility specifications, such as multiple stochastic volatility factors and jump components, in the appropriate modelling of equity return distributions. Furthermore, Johnson and Shanno [18] used the Monte Carlo method to price a European call option with the stochastic variance, while Wiggins [26] numerically solved the call option valuation problem under a general continuous stochastic process for return volatility.

For continuous sampling, Neuberger [22] proposed a nonparametric approach to study Delta hedging strategies based on variance swaps under a log contract. This method applies to arbitrary stochastic volatility processes, eliminating the need to assume a particular stochastic volatility model. In terms of discrete sampling, Elliot et al. [10] solved the pricing problem of swaps using probabilistic methods and partial differential equation methods. Additionally, Sepp [24] analysed the effect of discrete sampling on the valuation of options on the realized variance in the Heston stochastic volatility model. Short-term options on the realized variance can be priced by the semi-analytical Fourier transform methods, while Zhu and Lian [27] proposed a closed-form exact solution for pricing variance swaps under Heston’s two-factor stochastic volatility model based on the partial differential equation system. Further, Rujivan and Zhu [23] developed a simplified analytical approach and explored the relationship between the parameter space and effectiveness.

While all the aforementioned stochastic volatility models are driven by Brownian motion, empirical studies showed the presence of long-range correlation in the returns of stocks in financial markets [3, 17, 19]. To address this phenomenon, Mandelbrot and Van Ness [19] proposed fractional Brownian motion as a process based on the path integral form of standard Brownian motion. Decreusefond and Ustunel [8] used the stochastic calculus of variations to develop stochastic analysis theory for the functionals of fractional Brownian motions, while Duncan et al. [9] defined the multiple and iterated integrals of a fractional Brownian motion and provided various properties of these integrals. Elliott and Van Der Hoek [11] presented an extended framework for fractional Brownian motion, enabling processes with all indices to be considered under the same probability measure. Biagini et al. [4] introduced the theory of stochastic integration for fractional Brownian motion based on white-noise theory and differentiation (see, for example, [1]). As an application, Necula [21] generalized the risk-neutral valuation pricing formula in the framework of fractional Wick-type integrals [9]. Hu and Øksendal [16] proved that the fractional Black–Scholes market has no arbitrage if using the stochastic integration developed by Duncan et al. [9], contrary to the situation when the pathwise integration is used.

Gatheral et al. [13] showed that the fractional stochastic volatility (fSV) models have an excellent fit to financial time series data. Bayer et al. [2] showed how the rough fractional stochastic volatility model can be used to price claims, and they found that the rough Bergomi model with fewer parameters fits the volatility of the S&P 500 index markedly better than conventional Markovian stochastic volatility models. Cheridito et al. [6] proposed the fractional Ornstein–Uhlenbeck process [12] as a model for stochastic volatility, by proving the existence of a stationary solution to the Langevin equation [14] with fractional white noise. Garnier and Sølna [12] analysed the case where the stationary stochastic volatility model is constructed by a fractional Ornstein–Uhlenbeck process. However, some classical models (for example, the Stein–Stein model and the Heston model) introducing fractional Brownian motion have more general requirements for the volatility form. In this paper, we propose an approximate pricing method for the fractional stochastic volatility model by solving stochastic partial differential equations with variable coefficients, where the volatility is constructed as a deterministic function of time and the fractional Ornstein–Uhlenbeck process.

The paper is organized as follows. In Section 2, we introduce some basic background on the fractional Brownian motion and the fractional Ornstein–Uhlenbeck process. Our main results are presented in Section 3. We derive the approximate pricing formula, and prove that the approximation error can be limited. In Section 4, we calculate the price of the European option under the fractional Stein–Stein model [15] as an example to illustrate the feasibility and operability of the method. Numerical simulations are presented in Section 5 as well as the conclusions in Section 6.

2. Fractional Brownian motion and Ornstein–Uhlenbeck process

The fractional Brownian motion with Hurst parameter [19] $H\in (0,1) $ is a zero-mean Gaussian process $(B_{t}^{H} )_{t \in R } $ with covariance

$$ \begin{align*} E[ B_{t}^{H} B_{s}^{H}]=\frac{\sigma_{H}^{2}}{2}(|t|^{2 H}+|s|^{2 H}-|t-s|^{2 H}), \end{align*} $$

where

$$ \begin{align*} \sigma_{H}^{2} &=\frac{1}{\Gamma(H+1/2)^{2}}\bigg[\int_{0}^{\infty}((1+s)^{H-1/2}-s^{H-1/2})^{2} d s+\frac{1}{2 H}\bigg] \\ &=\frac{1}{\Gamma(2 H+1) \sin (\pi H)}. \end{align*} $$

We use the following moving-average stochastic integral representation of the fractional Brownian motion [16]:

$$ \begin{align*} B_{t}^{H}=\frac{1}{\Gamma(H+1/2)} \int_{R}(t-s)_{+}^{H-1/2}-(-s)_{+}^{H-1/2}\,d B_{s}, \end{align*} $$

where $(B_{t} )_{t \in R} $ is a standard Brownian motion. The filtration $\mathcal {F}_{t}$ generated by $B_{t}^{H}$ is also the one generated by $B_{t}.$

The fractional Brownian motion is self-similar, that is, $ (B_{ct}^{H},t \in R )$ and $ (c^{H} B_{t}^{H}, t \in R ) $ have the same probability law for all $c>0$ . Compared with Brownian motion, it displays a long-range dependence and positive correlation properties when $1/2<H<1$ , and it displays negative correlation property when $0<H<1/2$ . This special property of fractional Brownian motion allows it to describe path-dependent models.

The fractional Ornstein–Uhlenbeck process [12]

$$ \begin{align*} Z_{t}^{H}=\int_{-\infty}^{t} e^{-a(t-s)}\,d B_{s}^{H}=B_{t}^{H}-a \int_{-\infty}^{t} e^{-a(t-s)} B_{s}^{H}\,d s \end{align*} $$

is a zero-mean, stationary Gaussian process, with variance

$$ \begin{align*} \sigma_{\mathrm{ou}}^{2}=E[(Z_{t}^{H})^{2}]=\tfrac{1}{2} a^{-2 H} \Gamma(2 H+1) \sigma_{H}^{2}, \end{align*} $$

and covariance

$$ \begin{align*} E[Z_{t}^{H} Z_{t+s}^{H}]&=\sigma_{\text {ou }}^{2} \frac{1}{\Gamma(2 H+1)}\bigg[\frac{1}{2} \int_{R} e^{-|v|}|a s+v|^{2 H}\,d v-|a s|^{2 H}\bigg]\\ &=\sigma_{\text {ou }}^{2} \frac{2 \sin (\pi H)}{\pi} \int_{0}^{\infty} \cos (asx ) \frac{x^{1-2 H}}{1+x^{2}} \, d x. \end{align*} $$

Note that $Z_{t}^{H}$ is neither a martingale nor a Markov process. The fractional Ornstein–Uhlenbeck process also has the following representation of a moving-average integral:

$$ \begin{align*} Z_{t}^{H}=\int_{-\infty}^{t} \mathcal{K}(t-s)\, d B_{s}, \end{align*} $$

where

$$ \begin{align*} \mathcal{K}(t)=\frac{1}{\Gamma(H+1/2)}\bigg[t^{H-1/2}-a \int_{0}^{t}(t-s)^{H-1/2} e^{-a s} d s\bigg]. \end{align*} $$

3. Main results

Suppose that the market is self-financing and $Z_{t}^{H} $ is adapted to the Brownian motion $B^{\prime }_{t}$ . Here, $B_{t}$ and $B^{\prime }_{t}$ are two standard Brownian motions with correlation coefficient $\rho $ . In this section, we consider an option pricing problem, when the dynamics of the underlying asset is driven by the following stochastic differential equation:

(3.1) $$ \begin{align} \begin{cases} d X_{t}=\mu X_{t} \,d t+v_{t}X_{t} \,d B_{t}, \\ v_{t}=\bar{v}(t ) +F(\gamma Z_{t}^{H} ) , \end{cases} \end{align} $$

where F and $\bar {v}$ are smooth, positive-valued functions, bounded away from zero, with bounded derivatives. Here, $\mu $ and $\gamma $ are two constants.

Our objective is to calculate the price of the following derivative:

(3.2) $$ \begin{align} W_{t} : =E[g(X_{T}) \mid \mathcal{F}_{t}]. \end{align} $$

For notational simplicity, we introduce the operator

$$ \begin{align*} \mathcal{L}_{\bar{v}(t )} =\partial_{t} +\mu x\partial_{x}+\tfrac{1}{2} \bar{v}(t )^{2}x^{2}\partial^{2}_{xx}. \end{align*} $$

The following theorem gives an approximate expression of derivative price when $\gamma $ is small.

Theorem 3.1. If the underlying asset and the derivative follow the dynamics given by equations (3.1) and (3.2), we approximate the price of the derivative as follows:

$$ \begin{align*} W_{t}=M(t,X_{t} )+ O(\gamma^{2} ), \end{align*} $$

where

$$ \begin{align*} M(t,X_{t} )&=M_{1} (t,X_{t} )\\ &\quad +a(t,X_{t} )\gamma\bar{v}(t )\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )+a(t,X_{t} )\gamma \rho M_{2} ( t,X_{t} )\\ &\quad +\gamma\phi_{t}M_{3} (t,X_{t} ) +\gamma M_{4} (t,X_{t} )+\gamma \rho M_{5} (t,X_{t} ), \end{align*} $$

and $\phi _{t}=E[\int _{t}^{T} Z_{s}^{H}\,d s \mid \mathcal {F}_{t}]$ . The other elements are deterministic, and can be solved by the following partial differential equation system:

(3.3) $$ \begin{align} \begin{cases} \mathcal{L}_{\bar{v}(t )}M_{1} (t,x ) =0 ,\\ \mathcal{L}_{\bar{v}(t )}M_{2} (t,x )=-\bar{v}(t )^{2} (x \partial_{x}(x^{2}\partial_{xx}^{2} ) ) M_{1}(t,x )\theta _{t,T},\\ \mathcal{L}_{\bar{v}(t )}M_{3} (t,x ) =-(x^{2}\partial_{xx}^{2} )M_{1}(t,x )[a(t,x ){\bar{v}}' (t )+\bar{v} (t )\mathcal{L}_{\bar{v}(t )}a (t,x ) ] ,\\ \mathcal{L}_{\bar{v}(t )}M_{4} (t,x )=- \bar{v}(t ) (x\partial_{x}) M_{3} (t,x )\theta_{t,T} ,\\ \mathcal{L}_{\bar{v}(t )}M_{5} (t,x )=- M_{2} (t,x ) \mathcal{L}_{\bar{v}(t )}a (t,x ) ,\\ (1-a(t,x ) )\bar{v}(t ) (x^{2}\partial^{2}_{xx}) M_{1} (t,x )-M_{3} (t,x )=0,\\ M_{1} (T,x )=g(x),\\ M_{2} (T,x )=M_{3} (T,x )=M_{4} (T,x )=M_{5} (T,x )=0, \end{cases} \end{align} $$

where $ \theta _{t, T}=\int _{0}^{T-t} \mathcal {K}(v) \,d v$ .

Proof. For the smooth function $M_{1} (t,x )$ , we have by Itô’s formula [1],

$$ \begin{align*} dM_{1} (t,X_{t} )&=\bigg(\gamma\bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{2} g^{\gamma}(Z_{t}^{H} ) }{2} \bigg) ( x^{2}\partial^{2}_{xx}) M_{1}(t,X_{t} )\,dt\\ &\quad+v_{t} (x\partial_{x}) M_{1} (t,X_{t} )\,dB_{t}, \end{align*} $$

and

$$ \begin{align*} d(\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) ) &= (x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )d\phi_{t}+\phi_{t}d[(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) ]\\ &\quad +d\phi_{t}d((x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) )\\ &=-Z_{t}^{H}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )dt+(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )d\psi_{t}\\ &\quad +\phi_{t}v_{t} (x \partial_{x}(x^{2}\partial_{xx}^{2} ) )M_{1}(t,X_{t} )\,dB_{t}\\ &\quad +\phi_{t}\bigg(\gamma\bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{2} g^{\gamma}(Z_{t}^{H} ) }{2} \bigg) (x^{2} \partial^{2}_{xx}(x^{2}\partial_{xx}^{2} ) )M_{1}(t,X_{t} )\,dt \\ &\quad +v_{t}(x \partial_{x}(x^{2}\partial_{xx}^{2} ) ) M_{1}(t,X_{t} ) d{\langle} \phi_{t},B_{t} \rangle, \end{align*} $$

where

$$ \begin{align*} \psi_{t}=E\bigg[\int_{0}^{T} Z_{s}^{H} \, d s \mid \mathcal{F}_{t}\bigg], \quad g^{\gamma}(y )=2\bar{v}(t )\frac{F(\gamma y) -\gamma y}{\gamma^{2}}+\frac{F(\gamma y)^{2}}{\gamma^{2}}. \end{align*} $$

Notice that $v(t)$ is a function of t, so we cannot eliminate the first term in $dM_{1} (t,X_{t} )$ using $d(\phi _{t}(x^{2}\partial _{xx}^{2} )M_{1}(t,X_{t} ) )$ only. This problem is solved by introducing the function $a(t,x)$ in equation (3.3). We have $\langle \phi _t, B_t\rangle =\rho \langle \psi _t, B^{\prime }_t\rangle $ , and therefore,

$$ \begin{align*} &d(M_{1} (t,X_{t} ) +a(t,X_{t} ) \gamma\bar{v} (t )\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) ) \\[2pt]&\quad=\bigg(\gamma (1-a(t,X_{t} ) )\bar{v}(t )Z_{t}^{H} + \frac{\gamma^{2} g^{\gamma}(Z_{t}^{H} ) }{2} \bigg) ( x^{2}\partial^{2}_{xx}) M_{1} (t,X_{t} )\,dt\\[2pt]&\qquad+a(t,X_{t} )\phi_{t}\bigg(\gamma^{2}\bar{v}(t )^{2} Z_{t}^{H}+\frac{\gamma^{3}\bar{v}(t ) g^{\gamma}(Z_{t}^{H} )}{2} \bigg) (x^{2} \partial^{2}_{xx}(x^{2}\partial_{xx}^{2} ) )M_{1}(t,X_{t} )\,dt \\[2pt]&\qquad+a(t,X_{t} )\gamma\rho\bar{v}(t )v_{t}X_{t}(\partial_{x}(x^{2}\partial_{xx}^{2} ) ) M_{1}(t,X_{t} ) d \langle \psi_{t},B^{\prime}_{t} \rangle\\[2pt]&\qquad+\gamma(\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) )da(t,X_{t} )\bar{v}(t )+M_{t}^{(1 ) }, \end{align*} $$

where $M_{t}^{(1 ) }$ is a martingale satisfying

$$ \begin{align*} dM_{t}^{(1 ) }&=a(t,X_{t} ) v_{t} X_{t}\partial_{x}M_{1} (t,X_{t} )\,dB_{t}\\[2pt]&\quad+a(t,X_{t} )\gamma\bar{v}(t )(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )\,d\psi_{t}\\[2pt]&\quad+a(t,X_{t} )\gamma\bar{v}(t )\phi_{t}v_{t}X_{t} (\partial_{x}(x^{2}\partial_{xx}^{2} ) )M_{1}(t,X_{t} )\,dB_{t}. \end{align*} $$

We partially eliminate the first term in $dM_{1} (t,X_{t} )$ . To eliminate it completely, we introduce $M_{2} (t,x )$ and $M_{3} (t,x )$ . Applying [12, Lemma A.1], notice that

$$ \begin{align*}d\langle\psi_t, B^{\prime}_t\rangle =\theta_{t, T}\,d t.\end{align*} $$

Thus, we can write

$$ \begin{align*} &d(M_{1} (t,X_{t} )+a(t,X_{t} )\gamma\bar{v}(t )\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )) \\[2pt]&\qquad+d(a(t,X_{t} )\gamma \rho M_{2} (t,X_{t} )+\gamma\phi_{t}M_{3} (t,X_{t} ) ) \\[2pt]&\quad=\gamma \rho M_{2} (t,X_{t} )da(t,X_{t} )+\gamma v_{t}(x\partial_{x}) M_{3} (t,X_{t} )\theta_{t, T}\,d t \\[2pt]&\qquad+dR_{t}^{(1) }+dM_{t}^{(1 ) }+dM_{t}^{(2 ) }, \end{align*} $$

where

$$ \begin{align*} dM_{t}^{(2) }&=\gamma(\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) )\bar{v} (t )v_{t}(x\partial_{x}) a(t,X_{t} )\,dB_{t}\\[2pt]&\quad+a(t,X_{t} )\gamma \rho v_{t} (x\partial_{x}) M_{2} (t,X_{t} )\,dB_{t}\\[2pt]&\quad+\gamma M_{3} (t,X_{t} )\,d\psi_{t}\\[2pt]&\quad+\gamma v_{t}\phi_{t}(x\partial_{x})M_{3} (t,X_{t} )\,dB_{t}, \end{align*} $$

and

$$ \begin{align*} dR_{t}^{( 1) }&=\frac{\gamma^{2} g^{\gamma}(Z_{t}^{H} )}{2} ( x^{2}\partial^{2}_{xx}) M_{1} ( t,X_{t} )\,dt\\&\quad+a ( t,X_{t} )\phi_{t}\bigg(\gamma^{2}\bar{v} ( t )^{2} Z_{t}^{H}+\frac{\gamma^{3}\bar{v} ( t ) g^{\gamma}(Z_{t}^{H} )}{2} \bigg) (x^{2} \partial^{2}_{xx}(x^{2}\partial_{xx}^{2} ) )M_{1} ( t,X_{t} )\, dt \\&\quad+a ( t,X_{t} )\gamma^{2}\rho\bar{v} ( t )Z_{t}^{H} (x \partial_{x}(x^{2}\partial_{xx}^{2} ) ) M_{1} ( t,X_{t} ) \theta_{t, T} \,d t \\&\quad+\bigg(a ( t,X_{t} )\gamma^{2} \rho \bar{v} ( t ) Z_{t}^{H}+\frac{a ( t,X_{t} )\gamma^{3} \rho g^{\gamma}(Z_{t}^{H} )}{2} \bigg) (x^{2} \partial^{2}_{xx}) M_{2} ( t,X_{t} )\,dt\\&\quad+\frac{1}{2}(\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} ) )\bar{v} ( t )\bigg(\gamma^{2}\bar{v} ( t ) Z_{t}^{H}+\frac{\gamma^{3} g^{\gamma}(Z_{t}^{H} )}{2} \bigg)(x^{2} \partial^{2}_{xx}) a ( t,X_{t} )\,dt\\&\quad+\phi_{t}\bigg(\gamma^{2}\bar{v} ( t ) Z_{t}^{H}+\frac{\gamma^{3} g^{\gamma}(Z_{t}^{H} )}{2} \bigg) ( x^{2}\partial^{2}_{xx}) M_{3} ( t,X_{t} )\,dt. \end{align*} $$

In summary,

$$ \begin{align*} &d(M_{1} (t,X_{t} )+a(t,X_{t} )\gamma\bar{v}(t )\phi_{t}(x^{2}\partial_{xx}^{2} )M_{1}(t,X_{t} )) \\&\qquad+d(a(t,X_{t} )\gamma \rho M_{2} (t,X_{t} )+\gamma\phi_{t}M_{3} (t,X_{t} )) \\&\qquad+d(\gamma M_{4} (t,X_{t} )+\gamma \rho M_{5} (t,X_{t} ) ) \\&\quad=dR_{t}^{(1) }+dR_{t}^{(2) }+dM_{t}^{(1 ) }+dM_{t}^{(2 ) }+dM_{t}^{(3 ) }, \end{align*} $$

where

$$ \begin{align*} dR_{t}^{(2) }&=M_{2} (t,X_{t} )\bigg(\gamma^{2} \rho \bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{3} \rho g^{\gamma}(Z_{t}^{H} )}{2} \bigg)\,dt\\&\quad+\gamma^{2}\rho Z_{t}^{H}(x\partial_{x}) M_{3} (t,X_{t} )\theta_{t, T} \,d t \\&\quad+\bigg(\gamma^{2}\bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{3} g^{\gamma}(Z_{t}^{H} )}{2} \bigg) (x^{2}\partial^{2}_{xx}) M_{4} ( t,X_{t} )\,dt\\&\quad+\bigg(\gamma^{2} \rho\bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{3}\rho g^{\gamma}(Z_{t}^{H} )}{2} \bigg) (x^{2}\partial^{2}_{xx}) M_{5} (t,X_{t} )\,dt,\\dM_{t}^{(3 ) }&=v_{t}X_{t}M_{2} (t,X_{t} )dB_{t}+\gamma v_{t} (x\partial_{x})M_{4} (t,X_{t} )\,dB_{t}\\&\quad+\gamma\rho v_{t} (x\partial_{x})M_{5} (t,X_{t} )\,dB_{t}. \end{align*} $$

Therefore,

$$ \begin{align*} dM(t,X_{t} )=dM_{t}^{(1 ) }+dM_{t}^{(2 ) }+dM_{t}^{(3 ) }+dR_{t}^{(1) }+dR_{t}^{(2) }, \end{align*} $$

where $M_{t}^{(1 ) }$ , $M_{t}^{(2 ) }$ , $M_{t}^{(3 ) }$ are martingales. Let

$$ \begin{align*} M_{t}=M_{t}^{(1 ) }+M_{t}^{(2 ) }+M_{t}^{(3 ) },\quad R_{t}=R_{t}^{(1) }+R_{t}^{(2) }. \end{align*} $$

In addition, according to equation (3.3), we have $M(T,X_{T} )=g(X_{T} ) $ and

$$ \begin{align*} W_{t}&=E[g(X_{T}) \mid \mathcal{F}_{t}]\\ &=E[M (T,X_{T} ) \mid \mathcal{F}_{t}]\\ &=M (t,X_{t} )+E[ M_{T}-M_{t} \mid \mathcal{F}_{t} ]+ E\bigg[ \int_{t}^{T} dR_{s} \mid \mathcal{F}_{t} \bigg]\\ &=M (t,X_{t} )+ E\bigg[ \int_{t}^{T} dR_{s} \mid \mathcal{F}_{t} \bigg]. \end{align*} $$

Note that $g^\gamma (y)$ is bounded uniformly in $\gamma $ by

$$ \begin{align*} |g^\gamma(y)| \leq(\|\bar{v}\|_{\infty}\|F^{\prime \prime}\|_{\infty}+\|F^{\prime}\|_{\infty}^2) y^2. \end{align*} $$

This completes the proof of Theorem 3.1, since $E[ \int _{t}^{T} dR_{s} \mid \mathcal {F}_{t} ]$ is of order $\gamma ^{2}$ .

4. European option pricing under the fractional Stein–Stein volatility model

In this section, we calculate the approximate price of a European call option with a strike price of K under the fractional Stein–Stein volatility model as an example, that is,

(4.1) $$ \begin{align} \begin{cases} d X_{t}=\mu X_{t} \,d t+ | v_{t} | X_{t} \,d B_{t}, \\ d v_{t}=\beta(\alpha-v_{t}) d t+\gamma \,d B_{t}^{H}, \end{cases} \end{align} $$

(4.2) $$ \begin{align} W_{t}=E[(X_{T}-K)^{+} \mid \mathcal{F}_{t}], \end{align} $$

where $X_{t}$ is risky asset price process. Here, $\mu $ is the drift rate of the risk asset price process. Since the volatility process is a mean-reverting process, $v_{t}$ tends towards a long-term value $\alpha $ with rate $\beta $ . Here, $\gamma $ is a constant and $ B_{t}^{H}$ is a fractional Brownian motion with Hurst parameter $H>1 / 2$ .

Lemma 4.1. The equation

(4.3) $$ \begin{align} dv_t=\beta (\alpha- v_t)\,dt+\gamma \,dB_{t}^{H} \end{align} $$

has a unique solution of the form

$$ \begin{align*} v_{t}&=e^{-\beta t} v_{0}+\beta \alpha \int_{0}^{t} e^{\beta(s-t)} \,d s+\gamma B_{t}^{H}-\beta \int_{0}^{t} \gamma e^{\beta(s-t)} B_{s}^{H} \,d s\\ &=\bar{v}(t ) +\gamma Z_{t}^{H}. \end{align*} $$

Proof. By using Itô’s formula under fractional Brownian motion [9, Theorem 4.3], proof of Lemma 4.1 can be obtained directly.

Proposition 4.2. If the underlying asset and the derivative follow the dynamics given by equations (4.1) and (3.2) instead of equations (3.1) and (3.2), then Theorem 3.1 still holds.

Proof. Applying Lemma 4.1 and noting that $| v_{t} |^2=v_{t}^2$ , we replace $v_{t}$ in the proof of Theorem 3.1 with $| v_{t} |$ , such as

$$ \begin{align*} dM_{1} (t,X_{t} )&=\bigg(\gamma\bar{v}(t ) Z_{t}^{H}+\frac{\gamma^{2} g^{\gamma}(Z_{t}^{H} ) }{2} \bigg) ( x^{2}\partial^{2}_{xx}) M_{1}(t,X_{t} )\,dt\\ &\quad+| v_{t} | (x\partial_{x}) M_{1} (t,X_{t} )\,dB_{t}. \end{align*} $$

Then, the proof of Proposition 4.2 is complete, since $M_{t}$ is still a martingale and $E[ \int _{t}^{T} dR_{s} \mid \mathcal {F}_{t} ]$ is still of order $\gamma ^{2}$ .

If the underlying asset and the derivative follow the dynamics given by equations (4.1) and (4.2), applying Proposition 4.2, our approximate pricing method proceeds as follows.

Step 1. We solve $M_{1} (t,x )$ by

(4.4) $$ \begin{align} \begin{cases} \mathcal{L}_{\bar{v}(t )}M_{1} (t,x ) =0, \\ M_{1} (T,x )=(x-K)^{+}. \end{cases} \end{align} $$

Let

$$ \begin{align*} \begin{cases} u=M_{1}(t,x ) ,\\ y=xe^{\mu(T-t ) }, \end{cases} \end{align*} $$

then equation (4.4) becomes

(4.5) $$ \begin{align} \begin{cases} \partial_{t}u +\frac{1}{2} \bar{v}(t )^{2}y^{2}\partial^{2}_{yy}u =0 ,\\ u| _{t=T}=(y-K)^{+}. \end{cases} \end{align} $$

Likewise, let $\tau =\int _{0}^{t} \bar {v}(s )^{2} \,ds$ , then equation (4.5) becomes

$$ \begin{align*} \begin{cases} \partial_{\tau}u +\frac{1}{2} y^{2}\partial^{2}_{yy}u=0, \\ u|_{\tau =\hat{T}}=(y-K)^{+}, \end{cases} \end{align*} $$

where $\hat {T}=\int _{0}^{T} \bar {v}(s )^{2}\,ds$ . Applying the Black–Scholes formula (see, [5, 20]),

$$ \begin{align*} M_{1}(t,x )&=u(y,\tau) \\ &=yN(\hat{d_{1}})-KN(\hat{d_{2}}) \\ &=xe^{\mu(T-t ) }N(\hat{d_{1}})-KN(\hat{d_{2}}), \end{align*} $$

where

$$ \begin{align*} \hat{d_{1}}=\frac{\ln{(y/K) }+ (\hat{T}-\tau )/2 }{\sqrt{ \hat{T}-\tau } } =\frac{\ln{(x/K) }+\mu (T-t ) +(1/2) \int_{t}^{T}\bar{v}(s )^{2}\,ds }{\sqrt{ \int_{t}^{T}\bar{v}(s )^{2}\,ds } }, \end{align*} $$

and

$$ \begin{align*} \hat{d_{2}}=\hat{d_{1}}-\sqrt{ \hat{T}-\tau }=\hat{d_{1}}-\sqrt{ \int_{t}^{T}\bar{v}(s )^{2}\,ds } , \quad N(x )=\frac{1}{\sqrt{2\pi } }\int_{-\infty }^{x} e^{-s^{2}/2 }\,ds. \end{align*} $$

Moreover,

$$ \begin{align*} \partial_{x}M_{1}(t,x )&=e^{\mu(T-t ) }N(\hat{d_{1}})+\frac{e^{\mu(T-t ) }e^{-\hat{d_{1}}^{2}/2 }}{\sqrt{2\pi \int_{t}^{T}\bar{v}(s )^{2}\,ds }}-\frac{Ke^{-\hat{d_{2}}^{2}/2 }}{x\sqrt{2\pi \int_{t}^{T}\bar{v}(s )^{2}\,ds }} ,\\\partial^{2}_{xx}M_{1}(t,x )&=\frac{e^{\mu(T-t ) }e^{-\hat{d_{1}}^{2}/2 }}{x\sqrt{2\pi \int_{t}^{T}\bar{v}(s )^{2}\,ds }}-\frac{e^{\mu(T-t ) }e^{-\hat{d_{1}}^{2}/2 }\hat{d_{1}}}{x\sqrt{2\pi }\int_{t}^{T}\bar{v}(s )^{2}\,ds } \\& \quad +\frac{Ke^{-\hat{d_{2}}^{2}/2 }\hat{d_{2}}}{x^{2}\sqrt{2\pi }\int_{t}^{T}\bar{v}(s )^{2}\,ds } +\frac{Ke^{-\hat{d_{2}}^{2}/2 }}{x^{2}\sqrt{2\pi \int_{t}^{T}\bar{v}(s )^{2}\,ds} }. \end{align*} $$

Step 2. We solve $M_{2} (t,x )$ by

(4.6) $$ \begin{align} \begin{cases} \mathcal{L}_{\bar{v}(t )}M_{2} (t,x ) =-\bar{v}(t )^{2} (x \partial_{x}(x^{2}\partial_{xx}^{2} ) ) M_{1}(t,x )\theta _{t,T} =\overline{M_{1}}(t,x ) , \\ M_{2} (T,x )=0. \end{cases} \end{align} $$

Let

$$ \begin{align*} \begin{cases} z=\text{ln}x, \\ \tau=T-t, \end{cases} \end{align*} $$

then equation (4.6) becomes

(4.7) $$ \begin{align} \begin{cases} \partial_{\tau}M_{2} -\dfrac{\bar{v}^{2}}{2} \partial^{2}_{zz}M_{2}-\bigg(\mu-\dfrac{\bar{v}^{2}}{2} \bigg)\partial_{z}M_{2} =\overline{M_{1}}, \\[8pt] M_{2}|_{\tau=0}=0. \end{cases} \end{align} $$

Let

$$ \begin{align*} M_{2}=ue^{\alpha \tau +\beta z}, \quad \zeta=\int_{0}^{\tau} \bar{v}^{2}\,ds, \end{align*} $$

where

$$ \begin{align*} \alpha=\bigg(\frac{1}{2}+\frac{\mu}{\bar{v}^{2}} \bigg) \bigg(\frac{3\mu}{2}-\frac{\bar{v}^{2}}{4} \bigg), \quad \beta=\frac{1}{2}+\frac{\mu}{\bar{v}^{2}}. \end{align*} $$

Then equation (4.7) becomes

(4.8) $$ \begin{align} \begin{cases} \partial_{\zeta}u - \frac{1}{2}\partial^{2}_{zz}u=\widetilde {M_{1}}(\zeta,z ), \\ u|_{\zeta=0}=0, \end{cases} \end{align} $$

where

$$ \begin{align*} \widetilde {M_{1}}(\zeta,z )=\frac{\overline{M_{1}}}{e^{\alpha \tau +\beta z}\,\bar{v}^{2}}. \end{align*} $$

We get the solution of equation (4.8) as follows:

$$ \begin{align*} u(\zeta,z ) =\int_{0}^{\zeta}\int_{z-(\zeta-m )/2 }^{z+(\zeta-m )/2}\widetilde {M_{1}}(m,n )\,dn\,dm, \end{align*} $$

$$ \begin{align*} M_{2}=e^{\alpha \tau +\beta z}\int_{0}^{\zeta}\int_{z-(\zeta-m )/2 }^{z+(\zeta-m )/2}\widetilde {M_{1}}(m,n )\,dn\,dm. \end{align*} $$

Step 3. We solve $M_{3} (t,x)$ and $a(t,x)$ by

(4.9) $$ \begin{align} \begin{cases} \mathcal{L}_{\bar{v}(t )}M_{3} (t,x ) =-(x^{2}\partial_{xx}^{2} )M_{1}(t,x )\big[a(t,x ){\bar{v}}' (t )+\bar{v} (t )\mathcal{L}_{\bar{v}(t )}a (t,x ) \big] , \\ (1-a(t,x ) )\bar{v}(t ) (x^{2}\partial^{2}_{xx}) M_{1} (t,x )-M_{3} (t,x )=0,\\ M_{3} (T,x )=0. \end{cases} \end{align} $$

Let

$$ \begin{align*} \tilde{a}(t,x ) =a(t,x )-1,\quad z=\text{ln}x, \end{align*} $$

then equation (4.9) is translated into

(4.10) $$ \begin{align} \begin{cases} \big[\bar{v}(t )\partial_{t}f(t,x )+\mu \bar{v}(t )\partial_{z}f(t,x )+\frac{1}{2}\bar{v}(t )^{3}\partial^{2}_{zz}f(t,x )-\frac{1}{2}\bar{v}(t )^{3}\partial_{z}f(t,x ) \big]\\ \times \tilde{a}(t,x )+\bar{v}(t )^{3}\partial_{z}f(t,x ) \partial_{z}\tilde{a}(t,x )={\bar{v}}' (t )f(t,x ) ,\\ M_{3} (t,x )=\tilde{a}(t,x )\bar{v}(t ) f(t,x ),\\ \tilde{a}(T,x )=0. \end{cases} \end{align} $$

The solutions of equation (4.10) are as follows:

$$ \begin{align*} \tilde{a}=e^{-\int_{0}^{z}m(t,s)\,ds }\int_{t}^{T} \int_{0}^{z} n(\tau ,s )\, ds \, d\tau ,\quad M_{3} (t,x )=\tilde{a}\bar{v}(t ) f(t,x ), \end{align*} $$

where

$$ \begin{align*} m(t,z ) &=\frac{\partial_{t}f}{\bar{v}^{2}\partial_{z}f}+\frac{\mu}{\bar{v}^{2}}+\frac{\partial^{2}_{zz}f}{2\partial_{z}f}-\frac{1}{2}, \\ n(t,z )&=-\partial_{t} (q(t,z ) e^{\int_{0}^{z}m(t,s )\,ds } ) ,\\ q(t,z )&=\frac{{\bar{v}}'f}{\bar{v}^{3}\partial_{z}f}. \end{align*} $$

Step 4. Using a similar approach as in Step 2, we solve the $M_{4} (t,x )$ and $M_{5} (t,x )$ by

(4.11) $$ \begin{align} \begin{cases} \mathcal{L}_{\bar{v}(t )}M_{4} (t,x ) =- \bar{v}(t ) (x\partial_{x}) M_{3} (t,x )\theta_{t,T} =\overline{M_{3}}, \\ \mathcal{L}_{\bar{v}(t )}M_{5} (t,x )=- M_{2} (t,x ) \mathcal{L}_{\bar{v}(t )}a (t,x ) =\overline{M_{2}},\\ M_{4} (T,x )=M_{5} (T,x )=0. \end{cases} \end{align} $$

We get the solution of the above equation as follows:

$$ \begin{align*} M_{4}=e^{\alpha \tau +\beta z}\int_{0}^{\zeta}\int_{z-(\zeta-m )/2 }^{z+(\zeta-m )/2}\widetilde {M_{3}}(m,n )\,dn\,dm, \end{align*} $$

$$ \begin{align*} M_{5}=e^{\alpha \tau +\beta z}\int_{0}^{\zeta}\int_{z-(\zeta-m )/2 }^{z+(\zeta-m )/2}\widetilde {M_{2}}(m,n )\,dn\, dm, \end{align*} $$

where

$$ \begin{align*} \widetilde {M_{3}}(\zeta,z )=\frac{\overline{M_{3}}}{e^{\alpha \tau +\beta z}\,\bar{v}^{2}}, \quad \widetilde {M_{2}}(\zeta,z )=\frac{\overline{M_{2}}}{e^{\alpha \tau +\beta z}\,\bar{v}^{2}}. \end{align*} $$

5. Numerical simulations

In this section, we compare the fractional Stein–Stein volatility model with different H. Taking European options as an example and applying the Proposition 4.2, we illustrate and analyse the properties of the model with different volatilities, maturities and strike prices. Subsequently, we fix other parameters and adjust $\gamma $ to illustrate the reliability of the asymptotic analysis.

To simplify the analysis, we set $t=0$ , $\mu =0$ and $\rho =0$ . In the following numerical examples, $H=0.5, \ 0.7, \ 0.9$ , respectively, and $X_0=50, \ \beta =0.5$ . Notice that the volatility process is an Ornstein–Uhlenbeck process when $H=0.5$ .

Table 1 Option prices with different H, $\alpha $ and K.

Figure 1 Option prices with different H, $\alpha $ and K.

First, we let $\gamma =0.1$ , $T=1$ and show the impact of K, $\alpha $ and H (see Table 1, Figures 1 and 2). For Figure 1, when $\alpha $ is small, the effect of H is weak. When $\alpha $ takes a larger value, the option prices under stochastic volatility models with different H reflect significant differences. Compared with the case where the volatility process is an Ornstein–Uhlenbeck process, when $\alpha =2.5$ , the option prices under stochastic volatility models with $H=0.7$ and $H=0.9$ are lower and higher, respectively. This indicates that option prices under this model are not positively or negatively correlated with the Hurst parameter. According to Lemma 4.1 and Proposition 4.2, $\alpha $ directly affects $\bar {v}(t )$ , and H affects

$$ \begin{align*}\phi_{t}=E\bigg[\int_{t}^{T} Z_{s}^{H}\,d s \mid \mathcal{F}_{t}\bigg]\end{align*} $$

by directly affecting $Z_{t}^{H}$ . When the other parameters except H are fixed, the solution of the partial differential equation system of equation (3.3) is fixed. Thus, H only affects $M (t,X_{t} )$ through $\phi _{t}$ . See Table 2 for the relationship between H and $\phi _{0}$ in this simulation. When $H=0.7$ , $\phi _{0}$ obtained from this simulation takes a higher value compared with the other two cases in Figure 1, which leads to lower approximation results. This result also shows that the method requires us to better work out the conditional expectation $\phi _{t}$ . Unlike the case with $\alpha =2.5$ , when $\alpha =0.5$ , we observe that options with different strike prices are influenced by K in different ways. More specifically, options with $K<50$ are more affected by changes in K, that is, the curve is steeper, while options with $K>50$ are less affected. In Figure 2, notice that the change of option price caused by the change of strike price narrows with the increase of $\alpha $ . Moreover, it is by no means the case that larger $\alpha $ leads to higher option prices. When the other parameters except $\alpha $ are fixed, the solution of the complex partial differential equation system of equation (3.3) is affected by $\bar {v}(t )$ . The solution and $\bar {v}(t )$ together lead to the complex result in Figure 2.

Second, we let $\gamma =0.1$ and show the impact of T (see Table 3 and Figure 3). In most cases (except when $H=0.9$ , $\alpha =2.5$ , $K=30,35$ ), the option prices with the same strike price increase as time to maturity T increases. As can be seen from the data in Table 3, option prices with different parameters have different sensitivities to T. In-the-money options are less sensitive to T compared to out-of-the-money options. Options with higher T and higher mean-reversion level $\alpha $ are less sensitive to K when other parameters are fixed. For short-term maturity option cases (when $T=0.25$ ), they are most sensitive to K and least affected by H.

Figure 2 Option prices with different H, $\alpha $ and K.

Table 2 Simulated $\phi_{0}=E[\int_{0}^{T} Z_{s}^{H}\,d s \mid {\mathcal{F}}_{0}]$ with different H and $\alpha $ .

Table 3 Option prices with different H, $\alpha $ , T and K.

Finally, we let $H=0.9$ and show the impact of $\gamma $ (see Table 4 and Figure 4). How $\gamma $ affects the option price depends on the value of $\alpha $ . For this part of the numerical simulation results, when $\alpha =0.5$ , an increasing $\gamma $ generally leads to a decrease in the option price. However, when $\alpha =2.5$ , an increasing $\gamma $ generally leads to an increase in the option price. Moreover, as the $\gamma $ decreases from 1 to 0.001, the option prices calculated by the approximation method in this paper gradually converge to a relatively stable level. This result is consistent with Proposition 4.2 we have obtained and demonstrates the reliability of the approximation method.

Figure 3 $H=0.9$ . Option prices with different T, $\alpha $ and K.

Table 4 Option prices with different $\alpha $ and K.

It is important to note that the requirement for $\gamma $ for the option price to reach a relatively stable level is related to the value of the $\alpha $ . When $\alpha =0.5$ , $\gamma =0.1$ is sufficient to make the option price reach a relatively stable level. However, when $\alpha =2.5$ , $\gamma =0.01$ is required to achieve this goal.

Figure 4 $\alpha =0.5$ . Option prices with different $\gamma $ and K.

6. Conclusions

In this paper, we investigate the problem of pricing derivatives under a fractional stochastic volatility model. We obtain a method for approximating the prices of derivatives where the stochastic volatility can be composed of deterministic functions of time and the fractional Ornstein–Uhlenbeck process. Some fractional stochastic volatility models can be generalized to this type of problem. As an example, we give an approximate pricing expression and numerical simulation of a European option under the fractional Stein–Stein model. Numerical simulation results demonstrate the impact of the parameters in the fractional stochastic volatility model on the option price. By numerical simulation, we also show that the price of the option can reach a relatively stable level as $\gamma $ decreases, which is consistent with the main results we have obtained and demonstrates the reliability of the approximation method.