Next Article in Journal
Geomechanical Properties of Deep-Sea Pore-Filled Methane Hydrate-Bearing Soils at Critical State Using DEM Analysis
Next Article in Special Issue
Pricing European Options under a Fuzzy Mixed Weighted Fractional Brownian Motion Model with Jumps
Previous Article in Journal
Sliding-Window TD-FrFT Algorithm for High-Precision Ranging of LFM Signals in the Presence of Impulse Noise
Previous Article in Special Issue
On Finite-Time Blow-Up Problem for Nonlinear Fractional Reaction Diffusion Equation: Analytical Results and Numerical Simulations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Option Pricing with Fractional Stochastic Volatilities and Jumps

1
School of Science, Xi’an University of Posts and Telecommunications, Xi’an 710121, China
2
School of Computer Science and Technology, Xi’an University of Posts and Telecommunications, Xi’an 710121, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(9), 680; https://doi.org/10.3390/fractalfract7090680
Submission received: 20 July 2023 / Revised: 30 August 2023 / Accepted: 8 September 2023 / Published: 11 September 2023

Abstract

:
Empirical studies suggest that asset price fluctuations exhibit “long memory”, “volatility smile”, “volatility clustering” and asset prices present “jump”. To fit the above empirical characteristics of the market, this paper proposes a fractional stochastic volatility jump-diffusion model by combining two fractional stochastic volatilities with mixed-exponential jumps. The characteristic function of the log-return is expressed in terms of the solution of two-dimensional fractional Riccati equations of which closed-form solution does not exist. To obtain the explicit characteristic function, we approximate the pricing model by a semimartingale and convert fractional Riccati equations into a classic PDE. By the multi-dimensional Feynman-Kac theorem and the affine structure of the approximate model, we obtain the solution of the PDE with which the explicit characteristic function and its cumulants are derived. Based on the derived characteristic function and Fourier cosine series expansion, we obtain approximate European options prices. By differential evolution algorithm, we calibrate our approximate model and its two nested models to S&P 500 index options and obtain optimal parameter estimates of these models. Numerical results demonstrate the pricing method is fast and accurate. Empirical results demonstrate our approximate model fits the market best among the three models.

1. Introduction

Many studies [1,2,3,4,5] suggest that asset price fluctuations exhibit “long memory”. In addition, recent empirical studies [6,7,8,9] show that the roughness of the volatility process is observed. Fractional stochastic volatility models driven by fractional Brownian motions with a Hurst index H ( 0 < H < 1 / 2 , H 1 / 2 ) have flourished in the financial field, which can fit “long memory” ( H > 1 / 2 ) and “the roughness of volatility” ( H < 1 / 2 ) well. However, almost all these works ignore jumps in asset prices. In fact, in a setting driven by standard Brownian motions, the studies in [10,11,12] by real data suggest the real asset price is not a continuous process but jumps occasionally. Moreover, extensive empirical studies [13,14,15] show that jumps and stochastic volatilities exist in the asset price. On the one hand, compared to the other jump distribution, the mixed-exponential jump [12] is more general in approximating asset return distributions; on the other hand, the double Heston model [16] exhibits good performance in fitting the long-term “volatility smile” and “volatility clustering”. Therefore, this paper considers a fractional adaptation of the double Heston mixed-exponential jump-diffusion model by replacing standard Brownian motions with fractional Brownian motions in two volatility processes.
Due to the introduction of fractional Brownian motions, the model is no longer Markovian, and the PDE discretization schemes are not easily applied to the fractional setting [17]. Monte Carlo simulation is usually used. By discretizing the stochastic integral representation of the process, Bennedsen et al. [18] developed a hybrid scheme for the fractional Bergomi model. Fukasawa and Hirano [19] refined the hybrid scheme in reference [20] by reducing and reusing random numbers. Mehrdoust et al. [20] developed an Euler scheme for the fractional Heston model in which asset price and volatility process are both driven by mixed fractional Brownian motions. Chronopoulou and Spiliopoulos [21] proposed a sequential Monte Carlo simulation for fractional stochastic volatility models. More Monte Carlo methods for fractional stochastic volatility models can be found in [22]. However, an enormous amount of simulation time makes the above methods time-consuming and are not suitable for practical use.
Given the characteristic function, the Fourier transform method is an effective method for pricing European-type options [23,24]. Carr and Madan [25] first applied fast Fourier transform algorithms to evaluate call options. Jackson et al. [26] proposed a Fourier space time-stepping method for pricing options with jump risk by representing the pricing problem as the solution of a PIDE. Lord et al. [27] introduced a novel convolution method for pricing options by taking the probability density as a function of the difference between two variables. Fang and Oosterlee [28] proposed an efficient Fourier-cosine series expansion for a class of stochastic volatility (jump-diffusion) model. More Fourier transform methods can be found in [29]. For the affine model within the semimartingale framework [30,31], the characteristic function of the log-return can be obtained by solving a Riccati equation. For a fractional stochastic volatility model with affine structure [32], the characteristic function is still associated with the solution of a Riccati equation, but the Riccati equation is replaced by a fractional Riccati equation [33].
Unlike the standard Riccati equation, there is no closed-form solution for the fractional Riccati equation. The Adams discretization scheme [34,35,36] is the classical numerical method to deal with fractional ODEs. However, the high complexity of the scheme makes the computation of option prices not feasible for most practitioners [37]. Furthermore, some analytical approximation methods [38,39,40] are employed to approximate the solution of the fractional Riccati equation. For example, Cang et al. [38] approximated the solution of the fractional Riccati equation by the homotopy analysis method. Based on the Laplace-Adomian-decomposition method, Jeng and Kilicman [40] provided a Padé approximation method to solve the fractional Riccati equation. Callegaro et al. [41] proposed a hybrid scheme to solve the fractional Riccati equation by combining power series expansion with Richardson Romberg extrapolation. Due to the fact that the obtained solution is usually a semi-analytical expression, analytical approximation methods are more efficient compared to the Adams discretization scheme.
However, by introducing jumps and two fractional Brownian motions into the double Heston model, the aforementioned Riccati equation will become more complex, and the model distributions may no longer be stable [42]. Recently, Wang and Guo [43] evaluated variance and volatility swaps under a similar model. Under the approximative fractional Brownian motion framework, they obtained pricing formulas of variance and volatility swaps by the forward characteristic function which was derived by the property of conditional expectation and governing PIDE. Different from the technology in [43], we evaluate options based on Fourier-cosine series expansion. We first approximate the pricing model by a semimartingale and then convert corresponding fractional Riccati equations into a classic PDE with which we obtain the characteristic function of the log-return.
The main goal of this paper is to provide an efficient method for solving the fractional Riccati equation in the characteristic function of the fractional stochastic volatility jump-diffusion model. Since long-memory models are capable of capturing long-range dependence in time series data, which is helpful in future predictions [44,45], this work focuses on the case with H > 1 / 2 . The main contributions of this paper are threefold. First, this paper proposes a fractional stochastic volatility jump-diffusion model with two fractional stochastic volatilities and mixed-exponential jumps in asset prices. Secondly, this paper provides a method for solving the fractional Riccati equation in the characteristic function of the log-return. Thirdly, this paper provides the calibration of the pricing model. The rest of the paper is organized as follows. Section 2 presents the pricing model and related fractional Riccati equations. Section 3 provides the model approximation and the derivation of the characteristic function. Section 4 dilates the pricing method for European options. Section 5 presents some numerical experiments. Section 6 provides the calibration of the pricing model and some empirical experiments. Section 7 concludes.

2. The Pricing Model and Fractional Riccati Equations

Assume that Ω , F , { F t } 0 t T , P is a complete probability space, where the filtration { F t } 0 t T is continuous on the right, F 0 contains all P -null sets, and P is a risk-neutral probability measure. Assume that the asset price S t follows the stochastic partial differential system
d S t S t - = ( r - λ δ ) d t + j = 1 2 V j t d W j t S + d k = 1 N t ( ζ k - 1 ) , d V 1 t = k 1 θ 1 - V 1 t d t + σ 1 V 1 t d B 1 t H 1 , d V 2 t = k 2 θ 2 - V 2 t d t + σ 2 V 2 t d B 2 t H 2 ,
where k j , θ j , σ j   ( j = 1 ,   2 ) denote the mean reversion speed, mean reversion level, and volatility of variance processes V j t ( j = 1 ,   2 ) , respectively. Assume 2 k j θ j σ j 2   ( j = 1 ,   2 ) to guarantee processes V j t ( j = 1 ,   2 ) positive. For j = 1 ,   2 , W j t S are standard Brownian motions and B j t H j are fractional Brownian motions with Hurst parameters H j ( 1 / 2 ,   1 ) . The fractional Brownian motions B j t H j ( j = 1 ,   2 ) are both centered Gaussian processes with the covariance function E B j t H j B j s H j = 1 2 t 2 H j + s 2 H j - ( t - s ) 2 H j ( j = 1 ,   2 ) . Suppose S 0 = S ,   V 10 = V 1 , V 20 = V 2 .
N t is a Poisson process with constant intensity λ > 0 . Let δ = E [ ζ - 1 ] , where ζ = ( ζ k ) k 1 is a sequence of i.i.d. nonnegative random variables, such that Y = ln ζ has a mixed-exponential distribution with the following probability density
f Y ( y ) = p u k = 1 m p k η k e - η k y I y 0 + q d l = 1 n q l θ ^ l e θ ^ l y I y < 0 ,
where η k > 1   ( k = 1 , , m ) , θ ^ l > 0   ( l = 1 , , n ) , p u 0 , q d = 1 p u 0 ,   k = 1 m p k = 1 ,   l = 1 n q l = 1 ,   p k ,   q l ( ,   ) . To make f Y ( y ) be a probability density, a necessary condition is p 1 > 0 , q 1 > 0 , k = 1 m p k η k 0 , l = 1 n q l θ ^ l 0 , and a sufficient condition is k = 1 m p k η k 0 , l = 1 n q l θ ^ l 0 . By Formula (2), one has
δ = p u k = 1 m p k η k 1 η k - 1 + q d l = 1 n q l θ ^ l 1 θ ^ l + 1 - 1 .
For j = 1 ,   2 , suppose the process W j t S and B j t H j are independent of ζ and N t , while W j t S and B j t H j are correlated such that Cov d W 1 t S , d B 1 t H 1 = ρ 1 d t , Cov d W 2 t S ,   d B 2 t H 2 = ρ 2 d t .
According to Mandelbrot and Van Ness [46], B j t H j   ( j = 1 , 2 ) can be decomposed as follows
B j t H j = 1 Γ ( H + 1 / 2 ) 0 ( t s ) H 1 / 2 ( s ) H 1 / 2 d W j s v + 1 Γ ( H + 1 / 2 ) 0 t ( t s ) H 1 / 2 d W j s v , j = 1 ,   2 ,
where Γ ( ) is the gamma function, W j s v are both standard Brownian motions. Since the processes
0 ( t s ) H 1 / 2 d W j s v , j = 1 ,   2
are key for fitting long memory. We introduce these processes in the model (1) and have
d S t S t = ( r λ δ ) d t + j = 1 2 V j t d W j t S + d k = 1 N t ( ζ k 1 ) , d V 1 t = 1 Γ ( H + 1 / 2 ) 0 t ( t s ) H 1 / 2 k 1 θ 1 V 1 s d s + 1 Γ ( H + 1 / 2 ) 0 t ( t s ) H 1 / 2 σ 1 V 1 s d W 1 s v , d V 2 t = 1 Γ ( H + 1 / 2 ) 0 t ( t s ) H 1 / 2 k 2 θ 2 V 2 s d s + 1 Γ ( H + 1 / 2 ) 0 t ( t s ) H 1 / 2 σ 2 V 2 s d W 2 s v .
Let x ( t ) = ln ( S t / S ) . Since the process W j t S   ( j = 1 ,   2 ) and B j t H j   ( j = 1 ,   2 ) are independent of ζ and N t , based on the affine structure of the model (6), the characteristic function ϕ ( u , τ ) can be written as
ϕ ( u ,   τ ) = E exp i u λ τ δ + i u k = 1 N τ Y k exp g 1 ( u ,   τ ) + g 2 ( u ,   τ ) + j = 1 2 h j ( u ,   τ ) V j ,
where i denotes the imaginary unit, h j ( u ,   τ ) = 1 Γ ( H + 1 / 2 ) 0 τ ( τ s ) H 1 / 2 l j ( u ,   s ) d s   ( j = 1 ,   2 ) ,   g j ( u ,   τ ) = k j θ j 0 τ h j ( u ,   s ) d s   ( j = 1 ,   2 ) , l j ( u ,   τ )   ( j = 1 ,   2 ) are the solutions of the following two-dimensional fractional Riccati equations
D H + 1 / 2 l 1 ( u ,   τ ) = 1 / 2 ( u 2 i u ) + ( i u ρ 1 σ 1 k 1 ) l 1 ( u ,   τ ) + 1 / 2 σ 1 2 l 1 2 ( u ,   τ ) , D H + 1 / 2 l 2 ( u ,   τ ) = 1 / 2 ( u 2 i u ) + ( i u ρ 2 σ 2 k 2 ) l 2 ( u ,   τ ) + 1 / 2 σ 2 2 l 2 2 ( u ,   τ ) , I 1 / 2 H l 1 ( u ,   0 ) = I 1 / 2 H l 2 ( u ,   0 ) = 0 ,
where
D H + 1 / 2 l j ( u ,   τ ) = 1 Γ ( 1 / 2 H ) d d τ 0 τ ( τ s ) H 1 / 2 l j ( u ,   s ) d s ,   j = 1 ,   2 ,
I 1 / 2 H l j ( u ,   τ ) = 1 Γ ( 1 / 2 H ) 0 τ ( τ s ) H 1 / 2 l j ( u ,   s ) d s ,   j = 1 ,   2 .

3. Model Approximation and the Derivation of the Characteristic Function

By introducing two perturbed parameters ε 1 ,   ε 2   ( 0 < ε 1 ,   ε 2 1 ) , B j t H j ( j = 1 ,   2 ) can be approximated as
B j t ε j , H j = 0 t ( t s + ε j ) H j 1 / 2 d W j s v , j = 1 ,   2 .
According to Thao [47], B j t ε j , H j ( j = 1 ,   2 ) is a semimartingale and B j t ε j , H j B j t H j   ( j = 1 ,   2 ) in L 2 ( Ω ) as ε 1 0 , ε 2 0 . Differentiate Formula (9), one has
d B j t ε j , H j = ( H j 1 / 2 ) ψ j t d t + ε j H j 1 / 2 d W ˜ j t v , j = 1 ,   2 ,
where
ψ j t = 0 t ( t s + ε j ) H j 3 / 2 d W j s v , j = 1 ,   2 .
Substituting Formula (10) into model (1) and rearranging the terms, model (6) can be approximated as follows
d S t ε 1 , ε 2 S t - ε 1 , ε 2 = ( r - λ δ ) d t + j = 1 2 V j t ε j d W j t S + d k = 1 N t ( ζ k - 1 ) , d V 1 t ε 1 = ( H 1 - 1 / 2 ) ψ 1 t σ 1 V 1 t ε 1 + k 1 θ 1 - V 1 t ε 1 d t + ε 1 H 1 - 1 / 2 σ 1 V 1 t ε 1 d W ˜ 1 t v , d V 2 t ε 2 = ( H 2 - 1 / 2 ) ψ 2 t σ 2 V 2 t ε 2 + k 2 θ 2 - V 2 t ε 2 d t + ε 2 H 2 - 1 / 2 σ 2 V 2 t ε 2 d W ˜ 2 t v .
For convenience, we call model (12) as the FDHestonMEM. Different from the model (6), the FDHestonMEM is a semimartingale, and Itô calculus can be applied to the stochastic differential equation in the FDHestonMEM, and thus the corresponding Riccati equations have the closed-form solution.
Remark. 
The model contains most existing models as special cases. For example, (1) the double Heston mixed-exponential jump-diffusion model driven by standard Brownian motions (DHestonMEM) by setting H 1 = H 2 = 1 / 2 ; (2) the double Heston model (DHeston) by setting H 1 = H 2 = 1 / 2 ,   λ = 0 .
Theorem 1.
The characteristic function of x ( t ) is given by
ϕ ( u ,   τ ) = exp i u ln ( S / K ) + A ( u ,   τ ) + j = 1 2 B j ( u ,   τ ) V j ,
where for j = 1 ,   2 ,
A ( u ,   τ ) = i u r τ + λ τ Λ + j = 1 2 k j θ j Δ j 2 ( D j + γ j ) τ 2 ln 1 G j e D j τ 1 G j ,
B j ( u ,   τ ) = j = 1 2 1 Δ j 2 ( D j + γ j ) 1 e D j τ 1 G j e D j τ ,
D j = γ j 2 4 a j b j ,   γ j = i u ρ j Δ j k j ,
a j = 1 2 Δ j 2 ,   b j = 1 2 u ( i + u ) ,   τ = T - t , Δ j = ε j H j 1 / 2 σ j ,
Λ = p u k = 1 m p k η k 1 η k i u + q d l = 1 n q l θ ^ l 1 θ ^ l + i u 1 i u δ .
Proof of Theorem 1.
By Formulas (2) and (3), direct calculation yields
E exp i u λ τ δ + i u j = 1 N τ Y j = λ τ p u k = 1 m p k η k 1 η k i u + q d l = 1 n q l θ ^ l 1 θ ^ l + i u 1 i u δ ,
Let ϕ ˜ ( ) be the characteristic function of the continuous part of the log-return in the FDHestonMEM. Applying the multi-dimensional Feynman-Kac theorem, we have the following PDE
0 = ϕ ˜ τ + r 1 2 ( V 1 + V 2 ) ϕ ˜ x + 1 2 ( V 1 + V 2 ) ϕ ˜ x x             + j = 1 2 k j ( θ j V j ) + ( H j 1 2 ) σ j ψ j V j ϕ ˜ V j + 1 2 σ j 2 ε j 2 H j 1 V j ϕ ˜ V j V j + + ρ j σ j V j ε j 2 H j 1 ϕ ˜ x V j , ϕ ˜ ( u ,   0 ) = e i u ln ( S / K ) .
According to the affine structure of the FDHestonMEM, we have
ϕ ˜ ( u ,   τ ) = exp i u ln ( S / K ) + A ˜ ( u ,   τ ) + j = 1 2 B j ( u ,   τ ) V j
with initial conditions A ˜ ( u ,   0 ) = 0 , B j ( u ,   0 ) = 0   ( j = 1 ,   2 ) .
Substituting Formula (16) into PDE (15) yields
0 = - A ˜ τ - B 1 τ V 1 - B 2 τ V 2 + r 1 2 ( V 1 + V 2 ) i u - 1 2 ( V 1 + V 2 ) u 2 + j = 1 2 k j ( θ j V j ) + ( H j 1 2 ) σ j ψ j V j B j + 1 2 σ j 2 ε j 2 H j 1 V j B j 2 + ρ j σ j V j ε j H j 1 / 2 B j i u .
Since ψ j = ψ j t   ( j = 1 ,   2 ) is a martingale, one has ψ j 0 = E ( ψ j t ) = 0   ( j = 1 ,   2 ) . From the above equation, one can obtain a system of ODEs as follows
A ˜ τ = i u r + k 1 θ 1 B 1 + k 2 θ 2 B 2 , B 1 τ = 1 2 σ 1 2 ε 1 2 H j 1 B 1 2 + ( ρ 1 σ 1 ε 1 H 1 1 / 2 k 1 ) B 1 + 1 2 i u ( i u 1 ) , B 2 τ = 1 2 σ 2 2 ε 2 2 H j 1 B 2 2 + ( ρ 2 σ 2 ε 2 H 2 1 / 2 k 2 ) B 2 + 1 2 i u ( i u 1 ) .
By solving ODEs (17), we can obtain
A ˜ ( u ,   τ ) = i u r τ + j = 1 2 k j θ j Δ j 2 ( D j + γ j ) τ 2 ln 1 G j e D j τ 1 G j ,
B j ( u ,   τ ) = j = 1 2 1 Δ j 2 ( D j + γ j ) 1 e D j τ 1 G j e D j τ ,   j = 1 ,   2 ,
D j = γ j 2 4 a j b j ,   γ j = i u ρ j Δ j k j ,   j = 1 ,   2 ,
a j = 1 2 Δ j 2 ,   b j = 1 2 u ( i + u ) ,   τ = T - t , Δ j = ε j H j 1 / 2 σ j ,   j = 1 ,   2 .
Combining Formulas (7), (14), and (16), Theorem 1 follows. □
Theorem 2.
Assume that the asset price S t is governed by the FDHestonMEM and c n = n ln ϕ u n u = 0 denotes the n-th cumulant of x ( t ) , one has
c 1 = ln ( S / K ) + j = 1 2 θ j V j 2 k j 1 e k j τ + τ { r θ 1 + θ 2 2 + λ [ p u k = 1 m p k η k q d l = 1 n q l θ ^ l ( p u k = 1 m p k η k 1 η k 1 + q d l = 1 n q l θ ^ l 1 θ ^ l + 1 1 ) ] } c 2 = λ τ p u k = 1 m 2 p k η k 2 + q d l = 1 n 2 q l θ ^ l 2 + j = 1 2 ( Π 0 j + Π 1 j e k j τ + Π 2 j e 2 k j τ ) ,
where for j = 1 ,   2 ,
Π 0 j = θ j τ + 1 k j ( θ j + V j θ j Δ j ρ j τ ) + 1 k j 2 Δ j ρ j ( 2 θ j V j ) + 1 4 θ j Δ j 2 τ Δ j 2 8 k j 3 ( 5 θ j 2 V j ) , Π 1 j = Δ j ρ j k j 2 ( V j 2 θ j ) + 1 k j ( 1 Δ j ρ j τ ) ( θ j V j ) + θ j Δ j 2 2 k j 3 + τ Δ j 2 ( θ j V j ) 2 k j 2 , Π 2 j = Δ j 2 8 k j 3 ( θ j 2 V j ) , Δ j = ε j H j 1 / 2 σ j .

4. The Pricing Method

Let ξ = ln S T ln K , z = ln ( S t / K ) , a European option price can be written by
V ( z , t ) = e r ( T t ) E [ α ( e ξ 1 ) + ] = e r ( T t ) K + α ( e ξ 1 ) + f ξ z d ξ
with strike price K and maturity T , where r is the risk-neutral interest rate, f ξ z denotes the probability density of ξ conditional on z , α = 1 for calls and α = 1 for puts.
According to the cosine expansion, one has
f ξ z = k = 0 F k ( z ) cos k π ξ m n m ,
where F k ( z ) = 2 n m m n ξ z cos k π ξ m n m d ξ . By choosing appropriate m ,   n , F k ( z ) can be approximated by F ˜ k ( z ) as follows
F ˜ k ( z ) = 2 n m + f ( ξ z ) cos k π ξ m n m d ξ = 2 n m ϕ k π n m e i k π m n m ,
where ( ) denotes real part. By replacing F k ( z ) by F ˜ k ( z ) in (20), one has
f ξ z 2 n m k = 0 ϕ k π n m e i k π m n m cos k π ( ξ m ) n m .
By choosing appropriate N , one can truncate the series (22) and approximate f ξ z by the following Fourier cosine series expansion (COS)
f ξ z 2 n m k = 0 N 1 ϕ k π n m e i k π m n m cos k π ( ξ m ) n m ,
where ϕ ( ) denotes the characteristic function of the log-return, and the characteristic function can be obtained by Theorem 1 for the FDHestonMEM.
Inserting Formula (23) in Formula (19) and interchanging integration and summation, one has
V ( z , t ) 2 n m e r ( T t ) K k = 0 N 1 ϕ k π n m e i k π z m n m U k ,
where U k = m n α ( e ξ 1 ) + cos k π ( ξ m ) n m d ξ . For put options, basic calculus gives the following analytical solutions
U k = 1 1 + k π ( n m ) 2 k π n m sin m k π n m cos m k π n m + e m + m n k π sin m k π n m   k 0 , e m m 1                                                                                                                                                                 k = 0 .
According to [28], one can choose the integration interval [ m ,   n ] as follows
[ m ,   n ] = c 1 L c 2 ,   c 1 + L c 2 ,
where L is a constant, c 1 = ln ϕ u u = 0 , c 2 = 2 ln ϕ u 2 u = 0 . For the FDHestonMEM, c 1 ,   c 2 can be obtained by Theorem 2.

5. Numerical Experiments

We use the FDHestonMEM, the characteristic function (13), and Formula (24) (COS-based method) to evaluate European puts. To test the effectiveness of the pricing method, we take the numerical integration method as a benchmark and set 64 points in Formula (24) and L = 10 in Formula (25). According to [12] and [16], model parameters are set as follows:
k 1 = 12 , θ 1 = 0.05 , σ 1 = 0.9 , ρ 1 = 0.5 , V 1 = 0.05 , k 2 = 16 , θ 2 = 0.03 , σ 2 = 0.9 , ρ 2 = 0.5 , V 2 = 0.02 ,
λ = 1 , m = 2 , n = 2 , p u = 0.4 , p 1 = 1.3 , q 1 = 1.2 , θ ^ 1 = 20 , η 1 = 50 , θ ^ 2 = 20 , η 2 = 50 , ε 1 = 0.02 , ε 2 = 0.02 ,
H 1 = 0.8 , H 2 = 0.7 . For option parameters, we set   r = 0.0165 ,   S = 100 and specify three maturities T =   1 / 6 (short-term), T = 1 / 3 (medium-term), T = 1 (long-term) and set K = 80 ,   90 ,   100 ,   110 ,   120 for each maturity. We compare the two methods for evaluating European puts under the FDHestonMEM. Table 1 reports the main results.
From Table 1, we can see that the option prices calculated by the two methods are very close. The maximum relative error of the COS-based method under T = 1 / 6 , T = 1 / 3 , and T = 1 is 0.1932%, 0.0063%, and 0.0130%, respectively, and all the relative error of the COS-based method is not more than 0.1932%. The COS-based method presents high accuracy for pricing medium-term and long-term options. In terms of computational speed, the COS-based method is faster than numerical integration. To calculate an option price, the COS-based method needs 0.0203 s, while numerical integration almost needs three times that of the COS-based method. Table 1 shows that the COS-based method proposed in the paper is fast and accurate.

6. Calibration

We explore the calibration of the FDHestonMEM in this section. We use the implied volatility mean error sum of squares (IVMSE) as the loss function
IVMSE = 1 N T × N K t = 1 N T k = 1 N K ( I V t k I V t k Θ ) 2 ,
where I V t k Θ and I V t k are the implied volatilities from the model and market, respectively. We calibrate the FDHestonMEM by solving the following optimization problem
Θ * = arg min IVMSE ,
where Θ * is the optimal parameter vector.
We used S&P 500 index options on March 2, 2020, as market data. We filter market data according to the criteria in [48]. The filtered options have strike prices ranging from 2975 to 3150 and maturities ranging from 45 days to 377 days. The price of the underlying is 3060.725. According to the T-Bill interest rate in 2020 in the United States, we set 3.1% as the risk-free interest rate. We calculate model prices by the COS-based method. We calculate the market and model implied volatilities by putting the obtained market and model prices into the BS Formula [49], respectively. By the differential evolution algorithm, we solve the optimization problem (23). We also calibrate the DHestonMEM and the DHeston for comparison. Table 2 lists calibrated results and associated IVMSE for the three models.
Table 2 shows that the FDHestonMEM, the DHestonMEM, and the DHeston are all well calibrated. The IVMSE of these models are 1.871 × 10−7, 4.611 × 10−6, and 1.354 × 10−5, respectively. It suggests that the FDHestonMEM fits market data best among the three models; the second best is the DHestonMEM and the worst is the DHeston.
Furthermore, we set r = 0.031 ,   S = 3060.725 and specify four maturities T = 45, 107, 258, and 377 days, and set K = 2975, 3000, 3025, 3050, 3075, 3100, 3125, and 3150 for each maturity. Based on the parameters estimated in Table 2, we calculate European puts prices under the FDHestonMEM, the DHestonMEM, and the DHeston, respectively. With the obtained prices, we calculate the implied volatilities of these models and plot the implied volatility curves on strike prices under maturity T = 45, 107, 258, and 377 days. Figure 1 compares implied volatilities generated from these models and the market.
Figure 1 shows that the implied volatility curve of the FDHestonMEM is closest to that of the market in all maturities cases; the second closest is the DHestonMEM, and the farthest is the DHeston. The results are consistent with that in Table 2. Figure 1 further confirms that the FDHestonMEM fits the market best among the three models.
Let moneyness = S / K . We divide put options into three categories according to moneyness: out-of-the-money (OTM) if S / K 1.03 , at-the-money (ATM) if S / K ∈ (0.97, 1.03), and in-the-money (ITM) if S / K 0.97 . We set S = 100 ,   r = 0.031 ,   ε 1 = 0.01 ,   ε 2 = 0.01 and specify three maturities T = 1/12, 1/2, 1 years and set K = 90, 95, 100, 105, 110 for each maturity. With calibrated parameters in Table 2, we examine the influence of two Hurst parameters on OTM, ATM, and ITM prices under T = 1/12, 1/2, 1 years, respectively. Table 3, Table 4 and Table 5 report main results.
From Table 3, Table 4 and Table 5we can see that for all maturities cases, increasing H 2 lowers OTM and ATM prices but raises ITM prices when fixing H 1 = 0.5 and H 1 = 0.7 . When fixing H 1 = 0.9 , increasing H 2 presents the same influence on OTM and ITM as fixing H 1 = 0.5 and H 1 = 0.7 . However, for ATMs, increasing H 2 first raises option prices and then lowers option prices, which exhibits an inverted “smile” shape. This phenomenon is more obvious for long-term options. For short-term options, increasing H 1 lowers OTM prices but raises ATM and ITM prices when fixing H 2 . For medium-term options and long-term options, increasing H 1 raises option prices in all strike price cases when fixing H 2 . Table 3, Table 4 and Table 5 show that the influence of two Hurst parameters on options prices depends on maturity and moneyness.

7. Conclusions

This paper proposes a fractional stochastic volatility jump-diffusion model driven by two fractional Brownian motions. The characteristic function of the log-return is expressed in terms of the solution of two-dimensional fractional Riccati equations. By introducing two perturbed parameters, we approximate our pricing model by a semi-martingale and convert the fractional Riccati equations into a classic PDE. By the affine structure of the approximate model, we obtain the solution of the PDE with which we derive the characteristic function and its two cumulants. Based on the derived characteristic function and Fourier cosine series expansion, we obtain European option prices under the approximate model. Numerical results show that the proposed pricing method is fast and accurate. Our pricing method can be extended to stochastic rate models, such as the fractional stochastic volatility stochastic rate (jump-diffusion) model. Based on the differential evolution algorithm and S&P 500 index data, we calibrate the approximate model and its two nested models. Empirical analysis shows that our model is better than the double Heston mixed-exponential jump-diffusion model and the double Heston model in fitting the market. We also examine the influence of two Hurst parameters on option prices. Results show that the influence of Hurst parameters on options prices depends on maturity and moneyness.
The evaluation of path-dependent options under the fractional stochastic volatility jump-diffusion model is not trivial because the payoff of the type of option depends on the changing process of the asset price during the option’s lifetime. This is left for further research.

Author Contributions

Conceptualization, S.Z.; formal analysis, S.Z. and H.X.; methodology, S.Z. and H.Y.; validation, H.Y. and H.X.; writing–original draft, S.Z.; writing–review and editing, S.Z. and H.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation (No. 11601420) and the Natural Science Foundation of Shaanxi Province (No. 2020JM-577).

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the editors and anonymous reviewers for their detailed and valuable suggestions and comments, which improved the quality of the original manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Comte, F.; Renault, E. Long memory in continuous-time stochastic volatility models. Math. Financ. 1998, 8, 291–323. [Google Scholar] [CrossRef]
  2. Hassani, H.; Yarmohammadi, M.; Mashhad, L.M. Uncovering hidden insights with long-memory process detection: An in-depth overview. Risks 2023, 11, 113. [Google Scholar] [CrossRef]
  3. Han, Y.; Li, Z.; Liu, C. Option pricing under the fractional stochastic volatility model. Anziam, J. 2021, 63, 123–142. [Google Scholar]
  4. Najafi, A.; Mehrdoust, F. Conditional expectation strategy under the long memory Heston stochastic volatility model. Commun. Stat.-Simul. Comput. 2023, in press. [Google Scholar] [CrossRef]
  5. Alhagyan, M.; Yassen, M.F. Incorporating stochastic volatility and long memory into geometric Brownian motion model to forecast performance of Standard and Poor’s 500 index. AIMS Math. 2023, 8, 18581–18595. [Google Scholar] [CrossRef]
  6. Gatheral, J.; Jaisson, T.; Rosenbaum, M. Volatility is rough. Quant. Financ. 2018, 18, 933–949. [Google Scholar] [CrossRef]
  7. Livieri, G.; Mouti, S.; Pallavicini, A.; Rosenbaum, M. Rough volatility: Evidence from option prices. IISE Trans. 2018, 50, 767–776. [Google Scholar] [CrossRef]
  8. Fukasawa, M. Volatility has to be rough. Quant. Financ. 2021, 21, 1–8. [Google Scholar] [CrossRef]
  9. Brandi, G.; Di Matteo, T. Multiscaling and rough volatility: An empirical investigation. Int. Rev. Financ. Anal. 2022, 84, 102324. [Google Scholar] [CrossRef]
  10. Bates, D. Jumps and stochastic volatility: Exchange rate processes implicit in Deutsche mark options. Rev. Financ. Stud. 1996, 9, 69–107. [Google Scholar] [CrossRef]
  11. Kou, S.G. A jump-diffusion model for option pricing. Manag. Sci. 2002, 48, 1086–1101. [Google Scholar] [CrossRef]
  12. Cai, N.; Kou, S.G. Option pricing under a mixed-exponential jump diffusion model. Manag. Sci. 2011, 57, 2067–2081. [Google Scholar] [CrossRef]
  13. Coqueret, G.; Tavin, B. An investigation of model risk in a market with jumps and stochastic volatility. Eur. J. Oper. Res. 2016, 253, 648–658. [Google Scholar] [CrossRef]
  14. Jin, X.; Hong, Y. Jump-diffusion volatility models for variance swaps: An empirical performance analysis. Int. Rev. Financ. Anal. 2023, 87, 102606. [Google Scholar] [CrossRef]
  15. Bates, D.S. Empirical option pricing models. Annu. Rev. Financ. Econ. 2022, 14, 369–389. [Google Scholar] [CrossRef]
  16. Christoffersen, P.; Heston, S.; Jacobs, K. The shape and term structure of the index option smirk: Why multifactor stochastic volatility models work so well. Manag. Sci. 2009, 55, 1914–1932. [Google Scholar] [CrossRef]
  17. Bayer, C.; Friz, P.K.; Gulisashvili, A.; Horvath, B.; Stemper, B. Short-time near-the-money skew in rough fractional volatility models. Quant. Financ. 2019, 19, 779–798. [Google Scholar] [CrossRef]
  18. Bennedsen, M.; Lunde, A.; Pakkanen, M.S. Hybrid scheme for Brownian semistationary processes. Financ. Stoch. 2017, 21, 931–965. [Google Scholar] [CrossRef]
  19. Fukasawa, M.; Hirano, A. Refinement by reducing and reusing random numbers of the hybrid scheme for Brownian semistationary processes. Quant. Financ. 2021, 21, 1127–1146. [Google Scholar] [CrossRef]
  20. Mehrdoust, F.; Najafi, A.R.; Fallah, S.; Samimi, O. Mixed fractional Heston model and the pricing of American options. J. Comput. Appl. Math. 2018, 330, 141–154. [Google Scholar] [CrossRef]
  21. Chronopoulou, A.; Spiliopoulos, K. Sequential Monte Carlo for fractional stochastic volatility models. Quant. Financ. 2018, 18, 507–517. [Google Scholar] [CrossRef]
  22. Ma, J.T.; Wu, H.F. A fast algorithm for simulation of rough volatility models. Quant. Financ. 2022, 22, 447–462. [Google Scholar] [CrossRef]
  23. Zhang, S.M.; Sun, Y.D. Forward starting options pricing with double stochastic volatility, stochastic interest rates and double jumps. J. Comput. Appl. Math. 2017, 325, 34–41. [Google Scholar] [CrossRef]
  24. Zhang, S.M.; Geng, J.H. Fourier-cosine method for pricing forward starting options with stochastic volatility and jumps. Commun. Stat.-Theory Methods 2017, 46, 9995–10004. [Google Scholar] [CrossRef]
  25. Carr, P.P.; Madan, D.B. Option valuation using the fast Fourier transform. J. Comput. Financ. 2001, 2, 61–73. [Google Scholar] [CrossRef]
  26. Jackson, K.R.; Jaimungal, S.; Surkov, V. Fourier space time-stepping for option pricing with Lévy models. J. Comput. Financ. 2008, 12, 1–29. [Google Scholar] [CrossRef]
  27. Lord, R.; Fang, F.; Bervoets, G.; Oosterlee, C.W. A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes. SIAM J. Sci. Comput. 2008, 30, 1678–1705. [Google Scholar] [CrossRef]
  28. Fang, F.; Oosterlee, C.W. A novel pricing method for European options based on Fourier cosine series expansions. SIAM J. Sci. Comput. 2008, 31, 826–848. [Google Scholar] [CrossRef]
  29. Kienitz, J.; Wetterau, D. Financial Modeling: Theory, Implementation and Practice (with Matlab Source); John Wiley & Sons Ltd.: Chichester, UK, 2012. [Google Scholar]
  30. Heston, S.L. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Rev. Financ. Stud. 1993, 6, 327–343. [Google Scholar] [CrossRef]
  31. Duffie, D.; Pan, J.; Singleton, K. Transform analysis and asset pricing for affine jump-diffusions. Econometrica 2000, 68, 1343–1376. [Google Scholar] [CrossRef]
  32. Comte, F.; Coutin, L.; Renault, E. Affine fractional stochastic volatility models. Ann. Financ. 2012, 8, 337–378. [Google Scholar] [CrossRef]
  33. Euch, O.E.; Rosenbaum, M. The characteristic function of rough Heston models. Math. Financ. 2019, 29, 3–38. [Google Scholar] [CrossRef]
  34. Dumitru, B.; Kai, D.; Enrico, S. Fractional Calculus: Models and Numerical Methods; World Scientific: Singapore, 2012. [Google Scholar]
  35. Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52. [Google Scholar] [CrossRef]
  36. Li, C.; Tao, C. On the fractional Adams method. Comput. Math. Appl. 2009, 58, 1573–1588. [Google Scholar] [CrossRef]
  37. Jeng, S.W.; Kilicman, A. Fractional Riccati equation and its applications to rough Heston model using numerical methods. Symmetry 2020, 12, 959. [Google Scholar] [CrossRef]
  38. Cang, J.; Tan, Y.; Xu, H.; Liao, S.J. Series solutions of non-linear Riccati differential equations with fractional order. Chaos Solitons Fractals 2009, 40, 1–9. [Google Scholar] [CrossRef]
  39. Jaber, E.A.; El Euch, O. Multifactor approximation of rough volatility models. SIAM J. Financ. Math. 2019, 10, 309–349. [Google Scholar] [CrossRef]
  40. Jeng, S.W.; Kilicman, A. Series expansion and fourth-order global Pade approximation for a rough Heston solution. Mathematics 2020, 8, 1968. [Google Scholar] [CrossRef]
  41. Callegaro, G.; Grasselli, M.; Paèes, G. Fast hybrid schemes for fractional Riccati equations (rough is not so tough). Math. Oper. Res. 2021, 46, 221–254. [Google Scholar] [CrossRef]
  42. Wang, L.; Xia, W.X. Power-type derivatives for rough volatility with jumps. J. Futures Mark. 2022, 42, 1369–1406. [Google Scholar] [CrossRef]
  43. Wang, K.; Guo, X.X. Valuations of variance and volatility swaps under double Heston jump-difusion model with approximative fractional stochastic volatility. Comput. Econ. 2023, 1–31. [Google Scholar] [CrossRef]
  44. Alhagyan, M. The effects of incorporating memory and stochastic volatility into GBM to forecast exchange rates of Euro. Alex. Eng. J. 2022, 61, 9601–9608. [Google Scholar] [CrossRef]
  45. Hassler, U.; Pohle, M.O. Forecasting under long memory. J. Financ. Econom. 2023, 21, 742–778. [Google Scholar] [CrossRef]
  46. Mandelbrot, B.B.; Van Ness, J.W. Fractional Brownian motions, fractional noises and applications. SIAM Rev. 1968, 10, 422–437. [Google Scholar] [CrossRef]
  47. Thao, T.H. An approximate approach to fractional analysis for finance. Nonlinear Anal.-Real World Appl. 2006, 7, 124–132. [Google Scholar] [CrossRef]
  48. Bakshi, G.; Cao, C.; Chen, Z. Empirical performance of alternative option pricing models. J. Comput. Financ. 1997, 52, 2003–2049. [Google Scholar] [CrossRef]
  49. Black, F.; Scholes, M. The pricing of options and corporate liabilities. J. Polit. Econ. 1973, 81, 637–654. [Google Scholar] [CrossRef]
Figure 1. Comparison of the implied volatilities generated from the FDHestonMEM, the DHestonMEM, the DHeston, and the market.
Figure 1. Comparison of the implied volatilities generated from the FDHestonMEM, the DHestonMEM, the DHeston, and the market.
Fractalfract 07 00680 g001
Table 1. Comparison of the COS-based method and numerical integration for evaluating European put options under the FDHestonMEM.
Table 1. Comparison of the COS-based method and numerical integration for evaluating European put options under the FDHestonMEM.
TKCOS-BasedNumerical IntegrationRelative Error (%)
1/6800.15490.15520.1932
901.13641.13720.0703
1004.49494.49550.0133
11011.096111.09770.0144
12019.960519.96110.0030
1/3800.67400.67400.0000
902.47052.47050.0000
1006.34386.34410.0047
11012.539512.54030.0063
12020.579120.58030.0058
1803.06903.06940.0130
906.17466.17480.0032
10010.678610.67890.0028
11016.515216.51580.0036
12023.492023.49270.0029
CPU(s) 0.02030.0590
Table 2. Calibration results for the FDHestonMEM, the DHestonMEM, and the DHeston.
Table 2. Calibration results for the FDHestonMEM, the DHestonMEM, and the DHeston.
ModelCalibrated ParametersIVMSE
FDHestonMEM k 1 θ 1 σ 1 V 1 ρ 1 1.871 × 10−7
2.29960.00301.99640.1362−0.8256
k 2 θ 2 σ 2 V 2 ρ 2
6.55780.05250.36690.0003−0.8398
λ p u p 1 q 1 H 1
0.08350.06131.3330−0.27270.5010
η 1 η 2 θ ^ 1 θ ^ 2 H 2
19.50245.44903.721244.20410.5615
DHestonMEM k 1 θ 1 σ 1 V 1 ρ 1 4.611 × 10−6
9.72000.05630.93440.1268−0.9990
k 2 θ 2 σ 2 V 2 ρ 2
0.04830.00011.99980.0218−0.8584
λ p u p 1 q 1
0.03400.0001−1.4998−1.4999
η 1 η 2 θ ^ 1 θ ^ 2
49.998736.66690.00072.0791
DHeston k 1 θ 1 σ 1 V 1 ρ 1 1.354 × 10−5
5.58840.00012.00000.1296−0.5705
k 2 θ 2 σ 2 V 2 ρ 2
49.99980.05450.78460.0001−0.9990
Table 3. The influence of Hurst parameters on short-term European puts prices under the FDHestonMEM with T = 1/12 years.
Table 3. The influence of Hurst parameters on short-term European puts prices under the FDHestonMEM with T = 1/12 years.
H1H2K = 90K = 95K = 100K = 105K = 110
0.50.511.51087.45573.99591.52740.3555
0.711.50767.44983.98471.53360.3752
0.911.50877.44753.98021.53380.3829
0.70.511.30697.44164.32622.12190.8242
0.711.31577.43754.32412.12480.8319
0.911.31427.43594.32312.12570.8349
0.90.511.24277.37024.37852.29981.0494
0.711.23707.36704.37862.30421.0556
0.911.23517.36564.37832.30521.0581
Table 4. The influence of Hurst parameters on medium-term European puts prices under the FDHestonMEM with T = 1/2 years.
Table 4. The influence of Hurst parameters on medium-term European puts prices under the FDHestonMEM with T = 1/2 years.
H1H2K = 90K = 95K = 100K = 105K = 110
0.50.515.397811.82418.69286.07354.0090
0.715.289011.73948.66286.12284.1395
0.915.238711.69518.63996.13174.1787
0.70.516.166812.84599.90037.38185.2935
0.716.100812.79099.87297.39455.3530
0.916.073712.76209.85447.38785.3707
0.90.516.224813.103810.36858.02566.0706
0.716.189313.082510.36898.05246.1243
0.916.171913.069710.36388.05776.1404
Table 5. The influence of Hurst parameters on long-term European puts prices under the FDHestonMEM with T = 1 years.
Table 5. The influence of Hurst parameters on long-term European puts prices under the FDHestonMEM with T = 1 years.
H1H2K = 90K = 95K = 100K = 105K = 110
0.50.518.246014.988512.07979.53647.3693
0.718.123214.911112.06509.59877.5113
0.918.063114.866812.04579.61017.5555
0.70.519.144616.024813.214910.72648.5628
0.719.051315.959813.191710.75638.6512
0.919.005415.923813.171610.75728.6756
0.90.519.446816.486213.828211.47179.4085
0.719.391816.458513.832411.51079.4860
0.919.363616.439213.825711.51759.5084
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, S.; Yong, H.; Xiao, H. Option Pricing with Fractional Stochastic Volatilities and Jumps. Fractal Fract. 2023, 7, 680. https://doi.org/10.3390/fractalfract7090680

AMA Style

Zhang S, Yong H, Xiao H. Option Pricing with Fractional Stochastic Volatilities and Jumps. Fractal and Fractional. 2023; 7(9):680. https://doi.org/10.3390/fractalfract7090680

Chicago/Turabian Style

Zhang, Sumei, Hongquan Yong, and Haiyang Xiao. 2023. "Option Pricing with Fractional Stochastic Volatilities and Jumps" Fractal and Fractional 7, no. 9: 680. https://doi.org/10.3390/fractalfract7090680

APA Style

Zhang, S., Yong, H., & Xiao, H. (2023). Option Pricing with Fractional Stochastic Volatilities and Jumps. Fractal and Fractional, 7(9), 680. https://doi.org/10.3390/fractalfract7090680

Article Metrics

Back to TopTop