Next Article in Journal
Particle Swarm Optimization Fractional Slope Entropy: A New Time Series Complexity Indicator for Bearing Fault Diagnosis
Previous Article in Journal
Variable Step Hybrid Block Method for the Approximation of Kepler Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simultaneous Identification of Volatility and Mean-Reverting Parameter for European Option under Fractional CKLS Model

School of Mathematics, Renmin University of China, Beijing 100872, China
*
Author to whom correspondence should be addressed.
Submission received: 21 May 2022 / Revised: 13 June 2022 / Accepted: 14 June 2022 / Published: 21 June 2022

Abstract

:
In this paper, we reconstruct the time-dependent volatility function of the underlying asset and the mean-reverting parameter γ of the interest rate for European options under the fractional Chan–Karolyi–Longstaff–Sanders (CKLS) stochastic interest rate model. Tikhonov regularization is used to solve the ill-posedness of the inverse problem. The existence and stability of the solution of the regularization problem are given. We employ the alternating direction method of multipliers (ADMM) to iteratively optimize the volatility function and the parameter γ . Finally, numerical simulations and the empirical analysis are presented to illustrate the efficiency of the proposed method.

1. Introduction

Black-Scholes (B-S) model [1] made a breakthrough in option pricing, which gives the analytic expression of European option based on the Geometric Brownian motion. However, B-S model can not accurately describe the growing option market due to the long-term dependence of the underlying logarithmic returns. Fractional Brownian motion(fBm) was employed to describe it by B H , H ( 0 , 1 ) which was suggested by Mandelbrot and Van Ness [2] as an integral representation of Brownian motion process. If H ( 1 / 2 , 1 ) , the increment of fBm has a positive correlation. Furthermore, it has a negative correlation when H ( 0 , 1 / 2 ) . Hu and Øskendal [3] applied Wick-It o ^ product to give the analytic expression of European option under the fBm framework. Necula [4] gave a risk-neutral valuation theorem when the underlying is driven by fBm B H , H ( 1 / 2 , 1 ) . Huang et al. [5] deduced the analytic formula of European option for the fractional Vasicek stochastic interest rate model.
The short-term riskless interest rate is one of the most basic and vital parameters in the financial market. A lot of models have been proposed to interpret its behavior such as Vasicek model [6], Cox–Ingersoll–Ross (CIR) model [7], Constant Elasticity of Variance (CEV) model [8] and so on. Chan et al. [9] popularized the above models to the Chan–Karolyi–Longstaff–Sanders (CKLS) model as
d r = a ( b r ) d t + ω r γ d B H
Note that Equation (1) transforms to Vasicek model and CIR model when γ = 0 and γ = 1 / 2 , respectively. Both models are widely used in market due to their analytical properties. However, for the general parameter γ except for γ = 0 (Vasicek), γ = 1 / 2 (CIR) and γ = 1 [10], the explict expression for the bond price is hard to obtain. Thus, some methods have been developed for the approximation to the closed form solution such as Monte Carlo simulations [11], the analytical approximation formula [12,13,14] and numerical solutions of the corresponding partial differential Equation (PDE) [15,16]. Mao et al. [17] considered the Euler–Maruyama (EM) scheme for the typical hybrid mean-reverting θ -process and the convergence of EM solution when γ [ 1 / 2 , 1 ] . Wu et al. [18] gave the analytical properties and the convergence of numerical solutions in probability with γ > 1 . As for the fractional CKLS model, Marie [19] discussed the local existence of solution in rough paths sense. Kubilius and Med z ˇ i u ¯ nas [20] found the sufficiently simple conditions when the solution of fractional CKLS for γ > 1 and H ( 1 / 2 , 1 ) is positive, and obtained the convergence rate of the backward Euler approximation method for solutions of fractional CKLS model.
Volatility describes the randomness of the underlying asset price, it can efficaciously reduce the risks caused by market prices fluctuations but it can not be observed directly in market. Thus, there are many researches estimating volatility function from market data [21,22,23]. On the other hand, the parameter γ controls the relationship between the volatility of interest changes and the level of rate. Therefore, the reconstruction of the volatility and parameter γ from market data is practical to financial market. Cristofol and Roques [24] showed that the nonlinear drift and diffusion coefficients are simultaneously uniquely determined by the observation of the expectation and variance of a one-dimensional I t o ^ diffusion process, and solved the corresponding inverse problem with PDE technics mainly by the strong parabolic maximum principle.
Regularization is a practical tool to treat the ill-posedness occurred in the inverse problem. Naumova and Pereverzyev [25] proposed a multi-penalty regularization with a component-wise penalization and exhibited the compensatory property of the regularized approximation. Hofmann et al. [26] used a two-parameter Tikhonov regularization with heuristic parameter choice rules to settle the ill-posed nonlinear inverse problem which is the simultaneous recovery of implied volatility and interest rate functions from market data. Georgiev and Vulkov [27] considered an algorithm of predictor-corrector type to recover the piecewise linear time-dependent volatility and risk-free interest rate functions. In this paper, we try to reconstruct the volatility function and parameter γ from market option prices by alternating direction method of multiplier (ADMM) and Euler–Lagrange iterative method [22], then present the empirical analysis of Shanghai Stock Exchange (SSE) 50ETF.
The rest of this paper is organized as follows. In Section 2, we give the PDE for European option pricing under the fractional CKLS stochastic interest rate model. Section 3 introduces the Tikhonov regularization and gives the existence and stability of the regularized problem, and presents the ADMM algorithm which is used to iteratively optimize the volatility and parameter γ . The numerical verification is demonstrated in Section 4 with synthetic and real market data. Finally, some conclusions are stated in Section 5.

2. Problem Formulation

In this section, we consider the fractional CKLS model and derive the PDE for pricing the European option. Let S and r be the underlying asset price and stochastic interest rate, respectively. Their dynamics can be expressed as follows
d S = r S d t + σ ( t ) S d B H 1 d r = a ( b r ) d t + ω r γ d B H 2
where H is the Hurst parameter and 1 / 2 < H < 1 , B H j = { B H j ( t ) , t 0 , j = 1 , 2 } are two correlated fractional Brown motions with C o v ( d B H 1 , d B H 2 ) = ρ ( d t ) 2 H , ( | ρ | < 1 ) . σ ( t ) is the volatility. a , b , ω are constants. γ implies the mean-reverting γ -process.
Firstly we discuss the PDE of zero coupon pricing. Assume that we have a zero coupon bond whose value is given by P = P ( r , t ) . Construct an imaginary portfolio Π 1 consisting of one bond P 1 that matures at time T 1 and a number of ( ) of the bond P 2 which expires at T 2 as follows
Π 1 = P 1 P 2
Using the fractional I t o ^ ’s Lemma, we get the following formula in a small time interval d t
d Π 1 = P 1 t d t + P 1 r d r + H ω 2 r 2 γ t 2 H 1 2 P 1 r 2 d t P 2 t d t + P 2 r d r + H ω 2 r 2 γ t 2 H 1 2 P 2 r 2 d t .
To eliminate the random term in (3), let
= P 1 / r P 2 / r .
The no-arbitrage condition says that
E ( d Π 1 ) = r Π 1 d t
which implies that
P 1 t + H ω 2 r 2 γ t 2 H 1 2 P 1 r 2 r P 1 / P 1 r = P 2 t + H ω 2 r 2 γ t 2 H 1 2 P 2 r 2 r P 2 / P 2 r .
Note that
P t + H ω 2 r 2 γ t 2 H 1 2 P r 2 r P / P r = f ( r , t )
from (5), and f ( r , t ) can be set as the following form
f ( r , t ) = a ( b r ) λ ω r γ
where λ is the introducing market price of risk. Then
P t + H ω 2 r 2 γ t 2 H 1 2 P r 2 r P + a ( b r ) + λ ω r γ P r = 0
with the condition as
P ( r , T ) = 1 .
Lemma 1.
Suppose that the underlying asset S and the interest rate r follow (2), then the European call option V ( S , r , t ) with strike K and maturity T satisfies the PDE
V t + H σ 2 ( t ) S 2 t 2 H 1 2 V S 2 + 2 H ρ σ ( t ) ω r γ s t 2 H 1 2 V S r + H ω 2 r 2 γ t 2 H 1 2 V r 2 + a ( b r ) + λ ω r γ V r + r S V S r V = 0
conditional to
V ( S , r , T ) = ( S K ) + .
Proof. 
Construct an imaginary portfolio Π 2 with one option V ( S , r , t ) , a number of ( 1 ) of the underlying asset S and ( 2 ) number of the zero coupon bond P ( r , t ) as follows
Π 2 = V 1 S 2 P .
Using the fractional I t o ^ ’s Lemma in a small time d t , we have
d Π 2 = V t + H σ 2 ( t ) S 2 t 2 H 1 2 V S r + H ω 2 r 2 γ t 2 H 1 2 V r 2 d t + V S 1 d S + V r 2 P r d r + 2 P t + H ω 2 r 2 γ t 2 H 1 2 P r 2 d t
To make this portfolio riskless over the time interval d t , we choose
1 = V S , 2 = V / r P / r
in (4), then combine the Formula (7)
d Π 2 = V t + H σ 2 ( t ) S 2 t 2 H 1 2 V S r + H ω 2 r 2 γ t 2 H 1 2 V r 2 d t V / r P / r r P a ( b r ) + λ ω r γ P r d t .
Using the no-arbitrage pricing principle (4) gives
V t + H σ 2 ( t ) S 2 t 2 H 1 2 V S 2 + 2 H ρ σ ( t ) ω r γ s t 2 H 1 2 V S r + H ω 2 r 2 γ t 2 H 1 2 V r 2 + a ( b r ) + λ ω r γ V r + r S V S r V = 0
Direct problem: To solve the option pricing problem (9) and (10), we set the boundary conditions as in [28]:
V ( 0 , r , t ) = 0 , S = 0 , V S ( S , r , t ) = 1 , S = S m a x , V r ( S , r , t ) = 0 , r = r m a x .
The boundary condition at r = 0 can be obtained by setting r = 0 in Equation (9) which leads to the dynamical boundary
V t + H σ 2 ( t ) S 2 t 2 H 1 2 V S 2 + a b V r = 0
We reverse the time as τ : = T t which means the residual time to maturity. Equation (9) becomes a forward PDE, the terminal condition (10) transforms to an initial term, and the conditions (14) and (15) remain the same. Thus, we get that
V τ = H σ 2 ( τ ) S 2 ( T τ ) 2 H 1 2 V S 2 + 2 H ρ σ ( τ ) ω r γ s ( T τ ) 2 H 1 2 V S r + H ω 2 r 2 γ ( T τ ) 2 H 1 2 V r 2 + a ( b r ) + λ ω r γ V r + r S V S r V
with the condition
V ( τ = 0 ) = ( S K ) + .
Then the European call option price can be uniquely determined by the initial and boundary problem as soon as the volatility and parameter γ are accessible. Furthermore, we apply an implicit finite difference method to numerically solve it.
Calibration problem: Assume that the option prices { V ^ i j , i = 1 , 2 , , N , j = 1 , 2 , , M i } with maturities { T i } and strike prices { K i j } in market are given. Market calibration requires us to find a volatility function σ ( τ ) and a parameter γ which make the recalculated option price V ( S 0 , r 0 , τ 0 , K i j , T i , σ ( τ ) , γ ) satisfy
V i j b V ( S 0 , r 0 , τ 0 , K i j , T i , σ , γ ) V i j a
where V i j b and V i j a represent the bid and ask price in market, respectively. To accord with (18), we employ the available data to minimize the following error function
Γ ( σ , γ ) = i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , γ ) V ^ i j ] 2
where V ^ i j = ( V i j a + V i j b ) / 2 , and V ( S 0 , r 0 , τ 0 , K i j , T i , σ , γ ) is the numerical solution of (16) and (17) with the recovered volatility σ ( τ ) and parameter γ .
The option data in market is not sufficient to determine the volatility σ ( τ ) and γ uniquely, and the minimum value of the optimization problem (19) depends on market data discontinuously, so the calibration problem is ill-posed. Regularization technique is practical to solve the ill-posedness generally.

3. Problem Regularization

In this section, we consider the Tikhonov regularization method and derive the existence and stability of the solution of the regularized problem.
To obtain the stable solution, we calibrate volatility function and γ by solving the following minimization problem
min σ , g ( γ ) i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , g ( γ ) ) V ^ i j ] 2 + α | | ( σ , g ( γ ) ) | | * 2
where | | ( σ , g ( γ ) ) | | * = ( | | σ | | 2 2 + | | σ | | 2 2 + | | g ( γ ) | | 2 2 ) 1 / 2 and g ( γ ) = r 0 γ , V ^ i j is the market data and α > 0 is the regularization parameter.

3.1. Existence and Stability of the Solutions of the Minimization Problem

Now we discuss the existence and stability results of solutions for the minimization problem. Denote W = [ 0 , T ] × ( 1 / 2 , ) and
F α ( σ , g ( γ ) , v ) = i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , g ( γ ) ) v i j ] 2 + α | | ( σ , g ( γ ) ) | | * 2
for any vector v = ( v 11 , , v 1 M 1 , , v N M N ) ( 0 , ) i = 1 N M i of option prices. In what follows, f a , v ( σ , g ( γ ) ) stands for a function L 2 ( W ) ( σ , g ( γ ) ) F α ( σ , g ( γ ) , v ) [ 0 , ) .
Theorem 1
(Existence). Let α , V ^ be fixed and the same as in problem (20). Put
m = inf ( σ , g ( γ ) ) L 2 ( W ) f α , V ^ ( σ , g ( γ ) ) .
Assume that for some c ( m , ) the set M α , V ^ ( c ) = f α , V ^ 1 ( [ 0 , c ] ) L 2 ( W ) is bounded and f α , V ^ is weakly lower semi-continuous. Then problem (20) has at least one solution which belongs to M α , V ^ ( c ) .
Proof. 
There exists a minimizing sequence ( σ k , g ( γ ) k ) M α , V ^ ( c ) since c > m . L 2 ( W ) is a Hilbert space. Every Hilbert space is a reflexive Banach space. In the reflexive Banach space every bounded sequence contains a weakly convergent subsequence. Therefore, ( σ k , g ( γ ) k ) contains subsequence ( σ k n , g ( γ ) k n ) weakly tending to some ( σ ˜ , g ˜ ( γ ) ) . By weak lower semi-continuity of f α , V ^ ,
f α , V ^ ( σ ˜ , g ˜ ( γ ) ) inf min n f α , V ^ ( σ k n , g ( γ ) k n ) = m ,
and by the same weak lower semi-continuity M α , V ^ ( c ) = f α , V ^ 1 ( [ 0 , c ] ) is weakly closed. Therefore, ( σ ˜ , g ˜ ( γ ) ) is the minimizer contained in M α , V ^ ( c ) . □
The following theorem is a direct consequence of Theorem 1.
Theorem 2
(Stability). Let α , V ^ , m and f α , V ^ be the same as in Theorem 1 and let assumptions stated there hold. Let further α k ( 0 , ) , V k ( 0 , ) i = 1 N M i , m k [ 0 , ) , c k ( 0 , ) , ( σ ˜ k , g ˜ ( γ ) ) L 2 ( W ) be such that:
(1) α k k α , V k k V ^ ;
(2) k N + , m k = inf ( σ , g ( γ ) ) L 2 ( W ) f α k , V k ( σ , g ( γ ) ) ;
(3) k N + , c k > m k , M α k , V k ( c k ) = f α k , V k 1 ( [ 0 , c k ] ) L 2 ( W ) is bounded and f α k , V k is weakly lower semi-continuous;
(4) k N + , ( σ ˜ k , g ˜ ( γ ) k ) M α k , V ^ ( c k ) arg min ( σ , g ( γ ) ) L 2 ( W ) f α k , V k ( σ , g ( γ ) ) .
Finally, assume that ( σ ˜ k , g ˜ ( γ ) k ) is bounded as well. Then one can extract weakly convergent subsequence ( σ ˜ k n , g ˜ ( γ ) k n ) such that its limit is the solution of problem (20).

3.2. ADMM Algorithm for the Regularized Problem

The augmented Lagrangian for the above optimization problem (20) is
min σ , g ( γ ) l μ ( σ , g ( γ ) ) : = β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , γ ) V ^ i j ] 2 + | | σ | | 2 2 + | | σ | | 2 2 + | | g ( γ ) | | 2 2 + i = 1 N j = 1 M i μ i j [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , γ ) V ^ i j ]
where μ i j > 0 is the Lagrange multiplier and β = 2 / α is the penalty parameter.
We can divide the minimization problem (21) into two sub-problems as
σ k + 1 = arg min σ l μ k ( σ , g ( γ ) k ) , g ( γ ) k + 1 = arg min g ( γ ) l μ k ( σ k + 1 , g ( γ ) ) .
The iteration process of μ i j is
μ i j k + 1 = μ i j k + η V ( S 0 , r 0 , τ 0 , K i j , T i , σ k + 1 , g ( γ ) k + 1 ) V ^ i j
where η > 0 is the step size. The Lagrange multiplier update commonly uses a step size equal to the penalty parameter β [29]. In order to improve the convergence in practice and reduce the dependence on the initial selection of the penalty parameter, a simple scheme of varying penalty parameter such as (3.13) in [29] usually works well. For simplicity, we choose the same step size as the penalty parameter here.
For the σ sub-problem
σ k + 1 = arg min σ { | | σ | | 2 2 + | | σ | | 2 2 + β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , g ( γ ) k ) V ^ i j ] 2 + i = 1 N j = 1 M i μ i j k [ V ( S 0 , r 0 , τ 0 , K i j , T i , σ , g ( γ ) k ) V ^ i j ] }
Rewrite the right hand side of (24) as an integral that is
σ k + 1 = arg min σ 0 T { β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] 2 δ ( τ τ 0 ) + σ 2 ( τ ) + σ ( τ ) 2 + i = 1 N j = 1 M i μ i j k [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] δ ( τ τ 0 ) } d τ
where δ ( τ τ 0 ) is the Dirac delta function.
Since the above scheme only consider the appointed time τ 0 = T t 0 , the recovered volatility σ ( τ ) can not be well applied in the future [22]. So we can employ the historical data to get the calibrated results as
σ k + 1 = arg min σ T T c u r T { β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] 2 σ 2 ( τ ) + σ ( τ ) 2 + i = 1 N j = 1 M i μ i j k [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] } d τ
where T c u r represents the current time.
Let
L ( τ , σ , σ ) = σ 2 ( τ ) + σ ( τ ) 2 + β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] 2 + i = 1 N j = 1 M i μ i j k [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ]
and apply the Euler–Lagrange equation to L ( τ , σ , σ ) , then
σ σ = β 2 i = 1 N j = 1 M i V σ ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) [ V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) V ^ i j ] + 1 2 i = 1 N j = 1 M i μ i j k V σ ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) .
The term V ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) / σ need to be calculated and its definition is
V σ ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) = d d ε V ( S 0 , r 0 , τ , K i j , T i , σ + ε δ , g ( γ ) k ) ε = 0 .
We can evaluate the variational derivative by solving a PDE with a similar form to (16). Now consider the partial differential operator as
P σ = τ + H σ 2 ( τ ) S 2 ( T τ ) 2 H 1 2 S 2 + 2 H ρ σ ( t ) ω r γ S ( T τ ) 2 H 1 2 S r + H ω 2 r 2 γ ( T τ ) 2 H 1 2 r 2 + [ a ( b r ) + λ ω r γ ] r + r S S r I .
where I is the identity operator.
The option price V with a disturbed volatility satisfies the following PDE
P σ + ε δ V ( S 0 , r 0 , τ , K i j , T i , σ + ε δ , g ( γ ) k ) = 0 .
where δ is the Dirac delta function. Then differentiating the formula (31) with respect to ε and estimating its value at ε = 0 , we get
P σ V σ ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) = 2 H σ ( τ ) S 2 ( T τ ) 2 H 1 2 S 2 2 H ρ ω r γ S ( T τ ) 2 H 1 2 V S r ( S 0 , r 0 , τ , K i j , T i , σ , g ( γ ) k ) .
To solve Equation (28) numerically, we can introduce a ’false’ time parameter θ to make the iteration smoother [30] as follows
σ n k σ n k 1 Δ θ σ n + 1 k 2 σ n k + σ n 1 k ( Δ t ) 2 + σ n k + z n k 1 = 0
where
z n k = β 2 i = 1 N j = 1 M i V σ ( S 0 , r 0 , τ n , K i j , T i , σ k , g ( γ ) k ) [ V ( S 0 , r 0 , τ n , K i j , T i , σ k , g ( γ ) k ) V ^ i j ] + 1 2 i = 1 N j = 1 M i μ i j k V σ ( S 0 , r 0 , τ n , K i j , T i , σ k , g ( γ ) k ) .
For the sub-problem of g ( γ )
g ( γ ) k + 1 = arg min g ( γ ) { g ( γ ) 2 + β 2 i = 1 N j = 1 M i [ V ( S 0 , r 0 , τ n , K i j , T i , σ k + 1 , g ( γ ) k ) V ^ i j ] 2 + i = 1 N j = 1 M i μ i j k [ V ( S 0 , r 0 , τ n , K i j , T i , σ k + 1 , g ( γ ) k ) V ^ i j ] }
then minimize the objective function (34) by Particle Swarm optimization (PSO) algorithm.
As for the convergence of the solutions generated by ADMM algorithm, we can similarly deduce it according to Theorem 1 in [23].

4. Computational Experiments

4.1. Numerical Simulation

Here we consider the numerical examples. Let K = 64 , 68 , 72 , 76 , 80 , T = 0.5 , 1 , then S m a x = 2 K m a x = 160 , r m a x = 0.25 , the initial value S 0 = 62 , r 0 = 0.025 , the parameters in the interest rate model are
a = 0.2 , b = 0.05 , ω = 0.3 , ρ = 0.4 , λ = 0.2 , H = 0.7 .
Example 1.
Reconstruct the volatility function σ ( τ ) in L 2 [ 0 , 1 ] and γ > 1 / 2 by
σ ( τ ) = 0.2 0.1 l n ( 1.5 + 3 τ ) , γ = 1.3214 .
The values of regularization parameter and Lagrange multiplier are β = 0.7 , μ 0 = 0.5 , η = 0.7 . The iteration stop criterion is m a x | σ d ( τ ) σ d 1 ( τ ) | < 2.5 × 10 6 and the maximum number of iterations is 200, where σ d ( τ ) presents the recovered volatility in the dth iteration. Table 1 gives the true European call option price by Lemma 1 with (35). Figure 1 shows the comparison between the recovered and exact volatility. The recovered value of parameter γ is 1.327. A common measure A E (absolute error) is used to evaluate the volatility of the iterative algorithm for γ in Table 2 when ξ = 0 . It is clear that our algorithm is able to satisfactorily reconstruct volatility function and γ.
To illustrate the stability of our algorithm, we add the noise with intensity ξ of 1%, 3% and 5% to the exact volatility function (35)
σ ξ = [ 0.2 0.1 l n ( 1.5 + 3 τ ) ] + ξ × r a n d ( 1 , 1 )
where r a n d ( 1 , 1 ) presents the random number uniformly distributed in the interval (−1,1). Using the noised option prices generated by the disturbed volatility σ ξ , we can calibrate the volatility σ ( τ ) and parameter γ again. Figure 2 plots the reconstructed volatility from the noised option prices along with the exact volatility. The recovered parameter γ is 1.171, 1.028 and 0.9512 for the noised cases of 1%, 3% and 5% respectively. Table 2 gives the RMSE and maximum absolute error of recovered volatility σ ( τ ) and the absolute errors results of γ. Furthermore, the situation of ξ = 0 represents the calibration results without noise. We can see that the recovered results are not very accurate, but they can illustrate the stability of our scheme since they are closely located around the exact values.
Example 2.
Recover the σ ( τ ) in L 2 [ 0 , 1 ] and γ > 1 / 2 by
σ ( τ ) = 0.1 + 0.5 ( τ 0.5 ) 2 , γ = 0.8 .
This volatility function shows a ’smile’ curve. We choose the parameters as β = 0.4 , μ 0 = 0.49 , η = 0.4 .
Figure 3 plots the recovered volatility function against the exact one. Table 2 presents the calibrated results for volatility. Root mean square error(RMSE) is used to judge the accuracy of reconstructed volatility σ ( τ ) with exact volatility σ ^ as
R M S E = 1 n i = 1 n ( σ i σ ^ i ) 2
The recovered parameter γ is 0.8020 and its A E is shown in Table 3. We can see that the proposed scheme can effectively reconstruct the volatility function in the error of 10 4 . The absolute error of γ is in the 10 3 . Similarly, we add a little disturbance with the intensity ξ of 1%, 3% and 5% to the exact volatility function (36). Figure 4 plots the reconstructed volatility from the noised option prices along with the noised volatility and the exact one. The recovered parameter γ is 0.811, 0.754 and 0.712 for the noised cases of 1%, 3% and 5% respectively. Table 3 gives the RMSE and maximum absolute error of recovered volatility σ ( τ ) and the absolute errors results of γ. As expected, the reconstructed volatility and γ from the noised option prices no longer fit very well to the true values due to the disturbance. However, the numerical results also imply the stability of our method as the reconstructed values are close to the exact ones.

4.2. Empirical Analysis

In this part, we try to recover the volatility and parameter γ from the market data of Shanghai Stock Exchange(SSE) 50ETF on 9 September 2021 [23]. The current price of the underlying asset is S 0 = 3.204 , the times to expiration are T 1 = 13 / 250 , T 2 = 48 / 250 , T 3 = 104 / 250 , T 4 = 195 / 250 , the initial interest rate r 0 = 2.323 % is selected according to Shanghai Inter Bank offered Rate (SHIBOR) 1M at 9 September. The regularization and ADMM parameters are chosen as
β = 0.6 , μ 0 = 0.1 , η = 0.6 .
Other parameters in our model are chosen by PSO algorithm as
a = 0.2322 , b = 0.0202 , ω = 0.3164 , ρ = 0.2503 , λ = 0.3113 , H = 0.5972 .
The market data shown in Table 4 is subdivided by cubic spline interpolation to improve the calibration accuracy. The RMSE between the recovered option prices and market prices is used to evaluate the precision of reconstruction results. Table 5 presents the RMSE and maximum absolute error of recalculated option prices, and gives the value of recovered parameter γ . Figure 5 and Figure 6 specifically report the calibration results. Compared with Table 6 in the empirical analysis of [23], we can see that the recalculated option prices with γ = 0.8451 are more accurate than the Vasicek model. It is clear that the recovered prices at T = T 1 are lower than the market value from Figure 6 due to the interpolation error.
Through the numerical simulation with two typical kinds of volatility functions and the empirical analysis for SSE 50ETF, we conclude that our method is effective and stable to recover the volatility function σ ( τ ) and parameter γ simultaneously. Furthermore, it has the potential to be applied in real market.

5. Conclusions

In this paper, we discuss the Tikhonov regularization for the simultaneous recovery of volatility function σ ( τ ) and mean-reverting parameter γ from a finite set of market observations by the fractional CKLS model. Firstly we introduce the direct and inverse problems in option market, respectively. Then Tikhonov regularization and ADMM algorithm are applied to solve the ill-posed calibration problem. Finally, several examples are given to demonstrate the accuracy and robustness of our method. Furthermore, empirical analysis implies that the fractional CKLS model with γ = 0.8451 is more suitable than fractional Vasicek model for SSE 50ETF.

Author Contributions

J.Z.: conceptualization, methodology, software, validation, writing—original draft preparation, writing—review and editing; Z.X.: conceptualization, methodology, funding acquisition, supervision, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by the National Natural Science Foundation of China (Grant No. 12071479).

Acknowledgments

The authors thank the referees for their comments and detailed suggestions. These have significantly improved the presentation of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Black, F.; Scholes, M. The pricing of option and corporate liabilities. J. Polit. Econ. 1973, 81, 637–654. [Google Scholar] [CrossRef] [Green Version]
  2. Mandelbrot, B.B.; Ness, J. Fractional brownian motions, fractional noises and applications. SIAM Rev. 1968, 10, 422–437. [Google Scholar] [CrossRef]
  3. Hu, Y.Z.; Øksendal, B. Fractional white noise calculus and applications to finance. Infin. Dimens. Anal. Quantum Probab. Relat. Top. 2003, 6, 1–32. [Google Scholar] [CrossRef]
  4. Necula, C. Option pricing in a fractional Brownian motion environment. Adv. Econ. Financ.-Res.—Dofin Work. Pap. Ser. 2008, 2, 259–273. [Google Scholar] [CrossRef]
  5. Huang, W.L.; Tao, X.X.; Li, S.H. Pricing Formulae for European Options under the Fractional Vasicek Interest Rate Model. Acta Math. Sin. 2012, 55, 219–230. [Google Scholar]
  6. Vasicek, O.A. An Equilibrium Characterization of the Term Structure. J. Financ. Econ. 1977, 5, 177–188. [Google Scholar] [CrossRef]
  7. Cox, J.C.; Ingersoll, J.E.; Ross, S.A. A theory of the term structure of interest rates. Econometrica 1985, 53, 385–407. [Google Scholar] [CrossRef]
  8. Cox, J.C. The constant elasticity of variance option pricing model. J. Portf. Manag. 1996, 23, 15–17. [Google Scholar] [CrossRef]
  9. Chan, K.C.; Karolyi, G.A.; Longstaff, F.A.; Sanders, A.B. An empirical comparison of alternative models of the short term interest rate. J. Financ. 1992, 47, 1209–1227. [Google Scholar] [CrossRef]
  10. Brennan, M.J.; Schwartz, E.S. Analyzing convertible bonds. J. Financ. Quant. Anal. 1980, 15, 907–929. [Google Scholar] [CrossRef]
  11. Glasserman, P. Monte Carlo Methods in Financial Engineering; Springer: Berlin, Germany, 2013. [Google Scholar]
  12. Choi, Y.; Wirjanto, T.S. An analytic approximation formula for pricing zero-coupon bonds. Financ. Res. Lett. 2007, 4, 116–126. [Google Scholar] [CrossRef]
  13. Stehlíková, B. A simple analytic approximation formula for the bond price in the Chan-Karolyi-Longstaff-Sanders model. Int. J. Numer. Anal. Model. Ser. B 2013, 4, 224–234. [Google Scholar]
  14. Stehlíková, B.; Ševcovič, D. Approximate formulae for pricing zero-coupon bonds and their asymptotic analysis. Int. J. Numer. Anal. Model. 2009, 6, 274–283. [Google Scholar]
  15. Nowman, K.B.; Sorwar, G. Pricing UK and US securities within the CKLS model Further results. Int. Rev. Financ. Anal. 1999, 8, 235–245. [Google Scholar] [CrossRef]
  16. Nowman, K.B.; Sorwar, G. Derivative prices from interest rate models: Results for Canada, Hong Kong, and United States. Int. Rev. Financ. Anal. 2005, 14, 428–438. [Google Scholar] [CrossRef]
  17. Mao, X.R.; Truman, A.; Yuan, C. Euler–Maruyama approximations in mean-reverting stochastic volatility model under regime-switching. J. Appl. Math. Stoch. Anal. 2006, 7, 1–20. [Google Scholar] [CrossRef]
  18. Wu, F.K.; Mao, X.R.; Chen, K. A highly sensitive mean-reverting process in finance and the Euler–Maruyama approximations. J. Math. Anal. Appl. 2008, 348, 540–554. [Google Scholar] [CrossRef] [Green Version]
  19. Marie, N. A generalized mean-reverting equation and applications. ESAIM Probab. Stat. 2014, 18, 799–828. [Google Scholar] [CrossRef] [Green Version]
  20. Kubilius, K.; Medžiūnas, A. Positive solutions of the fractional SDEs with fon-lipschitz diffusion coefficient. Mathematics 2021, 9, 18. [Google Scholar] [CrossRef]
  21. Jiang, L.S.; Tao, Y.S. Identifying the volatility of underlying assets from option prices. Inverse Probl. 2001, 17, 137–155. [Google Scholar]
  22. Xu, Z.L.; Jia, X.Y. The calibration of volatility for option pricing models with jump diffusion processes. Appl. Anal. 2019, 98, 810–827. [Google Scholar] [CrossRef]
  23. Zhao, J.J.; Xu, Z.L. Calibration of time-dependent volatility for European options under the fractional Vasicek model. AIMS Math. 2022, 7, 11053–11069. [Google Scholar] [CrossRef]
  24. Cristofol, M.; Roques, L. Simultaneous determination of the drift and diffusion coefficients in stochastic differential equations. Inverse Probl. 2017, 33, 095006. [Google Scholar] [CrossRef] [Green Version]
  25. Naumova, V.; Pereverzyev, S.V. Multi-penalty regularization with a component-wise penalization. Inverse Probl. 2013, 29, 075002. [Google Scholar] [CrossRef]
  26. Hofmann, C.; Hofmann, B.; Pichler, A. Simultaneous identification of volatility and interest rate functions-a two-parameter regularization approach. Electron. Trans. Numer. Anal. 2019, 51, 99–117. [Google Scholar] [CrossRef]
  27. Georgiev, S.G.; Vulkov, L.G. Simultaneous identification of time-dependent volatility and interest rate for European options. Electron. Trans. Numer. Anal. 2020, 2333, 1–8. [Google Scholar] [CrossRef]
  28. Haentjens, T.; Hout, K. ADI finite difference schemes for the Heston-Hull-White PDE. J. Comput. Financ. 2012, 16, 83–110. [Google Scholar] [CrossRef]
  29. Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; Eckstein, J. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Now Found. Trends 2010, 3, 1–122. [Google Scholar]
  30. Lagnado, R.; Osher, S. A technique for calibrating derivative security pricing models: Numerical solution of an inverse problem. J. Comput. Financ. 1997, 1, 13–25. [Google Scholar] [CrossRef]
Figure 1. Comparison between the recovered and exact volatility for Example 1.
Figure 1. Comparison between the recovered and exact volatility for Example 1.
Fractalfract 06 00344 g001
Figure 2. Recovered volatility function against exact one with noise.
Figure 2. Recovered volatility function against exact one with noise.
Fractalfract 06 00344 g002
Figure 3. Comparison between the recovered and exact volatility for Example 2.
Figure 3. Comparison between the recovered and exact volatility for Example 2.
Fractalfract 06 00344 g003
Figure 4. Recovered volatility function against exact one for Example 2.
Figure 4. Recovered volatility function against exact one for Example 2.
Fractalfract 06 00344 g004
Figure 5. Left: Market prices; Right: Recovered prices.
Figure 5. Left: Market prices; Right: Recovered prices.
Fractalfract 06 00344 g005
Figure 6. Comparison of SSE 50ETF option price with recovered prices.
Figure 6. Comparison of SSE 50ETF option price with recovered prices.
Fractalfract 06 00344 g006
Table 1. The exact option prices for Example 1.
Table 1. The exact option prices for Example 1.
K6468727680
T 1 = 0.5 1.40320.40940.09170.01650.0025
T 2 = 1.0 2.51491.14650.45550.16000.0506
Table 2. The numerical results of recovered σ ( τ ) and γ .
Table 2. The numerical results of recovered σ ( τ ) and γ .
max | σ d σ d 1 | RMSE max | σ σ ^ | γ AE
ξ = 0 1.9305 × 10 7 5.7762 × 10 4 1.0623 × 10 3 1.3270.0056
ξ = 1 % 3.8220 × 10 7 3.0564 × 10 3 9.3685 × 10 3 1.1710.1504
ξ = 3 % 7.7861 × 10 7 8.1093 × 10 3 2.8106 × 10 3 1.0280.2934
ξ = 5 % 1.2144 × 10 6 1.1219 × 10 2 4.6803 × 10 2 0.95120.3702
Table 3. The numerical results of recovered σ ( τ ) and γ .
Table 3. The numerical results of recovered σ ( τ ) and γ .
max | σ d σ d 1 | RMSE max | σ σ ^ | γ AE
ξ = 0 7.7267 × 10 7 7.7984 × 10 4 1.6680 × 10 3 0.8020.002
ξ = 1 % 7.6843 × 10 7 9.3685 × 10 3 2.1668 × 10 3 0.8110.011
ξ = 3 % 7.5211 × 10 7 5.2394 × 10 3 2.8106 × 10 3 0.7540.046
ξ = 5 % 7.4837 × 10 7 8.4450 × 10 3 4.6843 × 10 2 0.7120.088
Table 4. The SSE 50ETF option price at 9 September 2021.
Table 4. The SSE 50ETF option price at 9 September 2021.
Strike Price K V K , T 1 V K , T 2 V K , T 3 V K , T 4
2.850.39280.40130.41790.4401
2.900.34340.35480.37580.4019
2.950.29660.30820.33470.3659
3.000.24700.26470.29840.3334
3.100.15230.18370.22690.2695
3.200.07260.11920.16710.2167
3.300.02750.07070.12100.1727
3.400.00840.03930.08600.1343
3.500.00300.02070.05900.1038
3.600.00140.01070.04220.0787
3.700.00080.00590.02900.0591
Table 5. The calibration results for SSE 50ETF option.
Table 5. The calibration results for SSE 50ETF option.
max | σ d σ d 1 | RMSE max | σ σ ^ | γ
Empirical results4.99 × 10 7 0.01510.04020.8451
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhao, J.; Xu, Z. Simultaneous Identification of Volatility and Mean-Reverting Parameter for European Option under Fractional CKLS Model. Fractal Fract. 2022, 6, 344. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract6070344

AMA Style

Zhao J, Xu Z. Simultaneous Identification of Volatility and Mean-Reverting Parameter for European Option under Fractional CKLS Model. Fractal and Fractional. 2022; 6(7):344. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract6070344

Chicago/Turabian Style

Zhao, Jiajia, and Zuoliang Xu. 2022. "Simultaneous Identification of Volatility and Mean-Reverting Parameter for European Option under Fractional CKLS Model" Fractal and Fractional 6, no. 7: 344. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract6070344

Article Metrics

Back to TopTop