1. Introduction
The Black–Scholes–Merton pricing formula has long been regarded as a fundamental tool for derivatives analysis [Reference Black and Scholes5, Reference Merton20]. 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 [Reference Stein and Stein25] examined the stock price distribution under a diffusion process with a stochastic volatility parameter, known as the Stein–Stein model. Chernov et al. [Reference Chernov, Gallant, Ghysels and Tauchen7] 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 [Reference Johnson and Shanno18] used the Monte Carlo method to price a European call option with the stochastic variance, while Wiggins [Reference Wiggins26] numerically solved the call option valuation problem under a general continuous stochastic process for return volatility.
For continuous sampling, Neuberger [Reference Neuberger22] 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. [Reference Elliott, Siu and Chan10] solved the pricing problem of swaps using probabilistic methods and partial differential equation methods. Additionally, Sepp [Reference Sepp24] 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 [Reference Zhu and Lian27] 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 [Reference Rujivan and Zhu23] 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 [Reference Beben and Orłowski3, Reference Huang and Yang17, Reference Mandelbrot and Van Ness19]. To address this phenomenon, Mandelbrot and Van Ness [Reference Mandelbrot and Van Ness19] proposed fractional Brownian motion as a process based on the path integral form of standard Brownian motion. Decreusefond and Ustunel [Reference Decreusefond and Ustunel8] used the stochastic calculus of variations to develop stochastic analysis theory for the functionals of fractional Brownian motions, while Duncan et al. [Reference Duncan, Hu and Pasik-Duncan9] defined the multiple and iterated integrals of a fractional Brownian motion and provided various properties of these integrals. Elliott and Van Der Hoek [Reference Elliott and Van Der Hoek11] presented an extended framework for fractional Brownian motion, enabling processes with all indices to be considered under the same probability measure. Biagini et al. [Reference Biagini, Øksendal, Sulem and Wallner4] introduced the theory of stochastic integration for fractional Brownian motion based on white-noise theory and differentiation (see, for example, [Reference Alòs, Mazet and Nualart1]). As an application, Necula [Reference Necula21] generalized the risk-neutral valuation pricing formula in the framework of fractional Wick-type integrals [Reference Duncan, Hu and Pasik-Duncan9]. Hu and Øksendal [Reference Hu and Øksendal16] proved that the fractional Black–Scholes market has no arbitrage if using the stochastic integration developed by Duncan et al. [Reference Duncan, Hu and Pasik-Duncan9], contrary to the situation when the pathwise integration is used.
Gatheral et al. [Reference Gatheral, Jaisson and Rosenbaum13] showed that the fractional stochastic volatility (fSV) models have an excellent fit to financial time series data. Bayer et al. [Reference Bayer, Friz and Gatheral2] 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. [Reference Cheridito, Kawaguchi and Maejima6] proposed the fractional Ornstein–Uhlenbeck process [Reference Garnier and Sølna12] as a model for stochastic volatility, by proving the existence of a stationary solution to the Langevin equation [Reference Gillespie14] with fractional white noise. Garnier and Sølna [Reference Garnier and Sølna12] 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 [Reference Gulisashvili, Viens and Zhang15] 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 [Reference Mandelbrot and Van Ness19]
$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 [Reference Hu and Øksendal16]:
$$ \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 [Reference Garnier and Sølna12]
$$ \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:
$$ \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:
$$ \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:
$$ \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 [Reference Alòs, Mazet and Nualart1],
$$ \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 [Reference Garnier and Sølna12, 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,
$$ \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} $$
$$ \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
$$ \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 [Reference Duncan, Hu and Pasik-Duncan9, 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
$$ \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
$$ \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, [Reference Black and Scholes5, Reference Merton20]),
$$ \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
$$ \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
$$ \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
$$ \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
$$ \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
$$ \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
$$ \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$
.
Option prices with different H,
$\alpha $
and K.
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.
Option prices with different H,
$\alpha $
and K.
Simulated
$\phi_{0}=E[\int_{0}^{T} Z_{s}^{H}\,d s \mid {\mathcal{F}}_{0}]$
with different H and
$\alpha $
.
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.
$H=0.9$
. Option prices with different T,
$\alpha $
and K.
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.
$\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.
Acknowledgement
This work is partially supported by the National Natural Science Foundation of China (grant no. 11871244).
