Calibration of time-dependent volatility for European options under the fractional Vasicek model

: In this paper, we calibrate the time-dependent volatility function for European options under the fractional Vasicek interest rate model. A fully implicit ﬁnite di ﬀ erence method is applied to solve the partial di ﬀ erential equation of option pricing numerically. To ﬁnd the volatility function, we minimize a cost function that is the sum of the squared errors between the theoretical prices and market prices with Tikhonov L 2 regularization and L 1 / 2 regularization respectively. Finally numerical experiments with simulated and real market data verify the e ﬃ ciency of the proposed methods.


Introduction
Black and Scholes [1] proved that option price satisfies a partial differential equation (PDE) under the strict assumptions and proposed the corresponding option pricing formula.However, the growing option market contradicts the assumptions in Black-Scholes model such as the logarithmic return distribution of underlying assets which usually shows the self similarity and long-term dependence in the real market [2][3][4][5][6][7].A common method is to introduce another random term, such as stochastic interest rate or stochastic volatility, into the model of underlying price.Vasicek [8] derived a general form for the term structure of interest rates and gave the bond pricing formula in the form of a stochastic integral representation.Heston [28] assumed that the volatility follows a Cox-Ingersoll-Ross (CIR) process and deduced a closed-form pricing formula of European options.He and Chen [34] modified the Heston model with a stochastic mean-reversion level and gave a closed-form pricing formula based on the dimensional reduction technique.And there are many other researches under different stochastic volatility models [31][32][33].
Moreover, fractional Brownian motion (FBM) can describe the dynamics of underlying assets which has long-term dependence.Mehrdoust and Najaf [9] obtained the explicit solution of fractional Black-Scholes model with weak payoff function.But Rogers [10] proved that it allows arbitrage opportunities.In order to solve this problem and consider the characteristics of long memory, Duncan et al. [11], Hu and Øksendal [12] developed a fractional white noise calculus and applied it to the market, which is modeled by Wick-It ô type of stochastic differential equations driven by fractional Brownian motion B H (t)(1/2 < H < 1).The corresponding It ô fractional Black-Scholes market has no arbitrage, and it is complete in contrast to the situation using the pathwise integration.Then Necula [13] obtained a fractional Black-Scholes formula for option price at t ∈ [0, T ] and a risk-neutral valuation theorem for the underlying which is driven by a fractional Brownian motion.Merton [14] gave the pricing formula of zero coupon bonds under the assumption that the corporate liabilities process obeys geometric Brownian motion and the interest rate changes randomly.Huang et al. [15] considered a model for complete and continuous market, and assumed that the process of asset price follows the geometric fractional Brownian motion and the interest rate process follows fractional Vasicek Interest rate model.The pricing formulas of European call option and put option were derived by using quasi martingale and partial differential equation method, and the parity formula was further obtained.
Volatility is a crucial parameter that affects option pricing, it is used to describe the price change of the underlying assets, and investors can use it to avoid some risks of assets price fluctuation in the future spot market.But it can not be directly observed in the market like the stock price.It is a considerable issue in finance to estimate volatility from the market price.Lagnado and Osher [18] put forward a technique for calibrating derivative security pricing models with respect to observed market prices, and estimated volatility from price observations by solving the inverse problem associated with the parabolic partial differential equation governing arbitrage-free derivative security prices.Chiarella et al. [19] suggested an improvement on the basis of Lagnado and Osher, avoiding the possibility of negative volatility, and gave the numerical algorithm of Euler-Lagrange equation.Jiang and Tao [20] applied an optimal control framework to determine implied volatility and made a rigorous mathematical analysis of this inverse problem.Ngnepieba [21] employed the unconstrained minimization algorithm of the quasi-Newton limited memory type to find the optimal volatility function, and computed the gradient via the adjoint method.He and Zhu [30] developed a two-step approach to calibrate the local volatility under the regime-switching models.Georgiev and Vulkov [22] presented a robust and fast numerical algorithm to reconstruct the implied volatility as a piecewise linear function of time.In this paper, we try to calibrate the time-dependent volatility of European options from a set of market observations under the fractional Vasicek model, using the Euler-Lagrange iterative method [16] and alternating direction method of multiplier (ADMM) [25], and give an empirical analysis of China's option market.Moreover, the proposed schemes also apply to the stochastic volatility models with mean reversion such as Heston model [28].
This paper is organized as follows: In Section 2, we present the the Vasicek fractional model for European option pricing and the corresponding inverse problem.In Section 3, we introduce the Tikhonov L 2 and L 1/2 regularization and analyse the stability of regularization model and the convergence of ADMM.The numerical verification is demonstrated in Section 4 with synthetic and real market data.Finally, some conclusions are stated in Section 5.

The mathematical formulation
In this section, we discuss the Vasicek fractional model for the European option.The dynamics of the model can be described as where b = b − λ a τ and λ is the market price of risk.Lemma 1. (European call option pricing formula) Suppose that the underlying asset price S and interest rate r follow Vasicek fractional model (2.1), then the price of a European call option with strike price K and maturity T can be expressed as where P (r, t; T ) is the zero coupon bond P(r, t; T ) = e −rB(t,T )−A(t,T ) , B(t, T ) = 1 a (1 − e −a(T −t) ), and N(•) denotes the standard normal cumulative distribution function.
In this paper, we only consider the European call option, the put option can be obtained by the putcall parity [15].If the volatility σ(t) is specified, the option price V(S 0 , r 0 , 0, K, T, σ) can be uniquely determined by Lemma 1. Assume that the option prices {V i j , i = 1, 2, ..., N, j = 1, 2, ..., M i } with maturities {T i } and the corresponding strike prices {K i j } are known, the inverse problem requires us to find a volatility function which makes the calibrated option value satisfy the following market quotes where V b i j and V a i j represent the bid and ask price in market respectively.Let H denote the measurable function space defined in [0, T max ], where T max is the longest maturity.To satisfy (2.5), we use the available data to solve the following nonlinear problem.
Calibration problem: Given a series of market option prices {V i j , i = 1, 2, ..., N, j = 1, 2, ..., M i }, find a suitable volatility function σ(t) ∈ H to minimize the following error function where , and the theoretical option value V(S 0 , r 0 , 0, K i j , T i , σ(t)) is obtained by Lemma 1.
However, there are not sufficient market observations to confirm the volatility σ(t) uniquely, and the value of Π(σ) is discontinuously dependent on the market data.Thus the calibration problem is ill-posed.In general, regularization is a practical technique to solve ill-posedness.

The regularization schemes for calibration problem
We employ two regularization schemes due to the ill-posedness of the volatility function minimizer.One is the Tikhonov

The Tikhonov L 2 regularization strategy
The Tikhonov L 2 regularization term we add is based on [16], where 1 2 , V i j is a set of option prices in market and α is the regularization parameter.So we have Now, rewrite the right hand side of (3.2) as an integral where δ(t) is the Dirac delta function.Suppose that We apply the Euler-Lagrange equations to L(t, σ, σ ) and obtain with δV i j δσ = ∂ ∂σ V(S 0 , r 0 , 0, K i j , T i , σ)δ(t).However, the above strategy calibrates volatility function at the fixed time t 0 = 0. Thus, the function σ(t) derived from this method does not necessarily apply to future time.We can calibrate the volatility from the historical data [16], that is where T cur is the current time.
Then apply L 2 regularization to minimize the following function Now let us rewrite J α (σ) as the integral form Similarly, the Euler-Lagrange equation for (3.8) is with initial conditions σ(t = 0) which can be determined from market data.To solve Eq (3.9), we need to evaluate the term ∂V(t; S 0 , r 0 , K i j , T i , σ)/∂σ.This variational derivative is defined for general arguments by where the perturbation δ(t) is a Dirac delta function [18].Such variational derivative can be calculated by solving a PDE (2.2) with an additional source term.Now consider the following partial differential operator where I denotes the identity operator.
The pricing function with a perturbed volatility satisfies a PDE similar to (2.2) of the form P σ+εδ V(t; S 0 , r 0 , K i j , T i , σ + εδ) = 0. (3.12) Then differentiate (3.12) with respect to ε and evaluate the value when ε = 0, we have The variational derivative satisfies homogeneous boundary and initial conditions [18], we can solve (3.13) for ∂V/∂σ numerically.

The L 1/2 regularization strategy
In this subsection, the L 1/2 regularization is considered, i.e., we calibrate volatility function by minimizing the following function where and ξ, η are regularization parameters.

Numerical algorithm for L 1/2
To calibrate the time-dependent volatility σ (t) from market data by L 1/2 regularization, we employ ADMM algorithm to obtain the convergent solution.
The augmented Lagrangian for the nonconvex optimization problem in (3.1) is min σ,y,z where µ > 0 is a regularization parameter.
Then we solve the augmented Lagrangian (3.15).Concretely, the procedure is as follows However, the algorithm (3.16) can not guarantee sufficient reduction for y k+1 and z k+1 .Based on [25], we have (3.17) For the σ sub-problems, The corresponding Euler-Lagrange equation for (3.18) is Since the sub-problems for y and z are nonconvex, nonsmooth and discontinuous, Xu et al. [24] proposed a reweighted iterative algorithm to solve this problem, which transforms L 1/.2 regularization into L 1 term where ε is a sufficient small constant to avoid the case |y k i | = 0. Then the following threshold iterative Algorithm 1 gives the local minimizer of y by applying accelerated approximate gradient method [27], where L denotes the Lipschitz constant.

Convergence analysis
Now we discuss the convergence of the above algorithm [25,26].

Computational experiments
In this section, we use numerical simulation and empirical analysis to verify the effectiveness of the proposed methods.
Table 2 compares the reconstructed and real volatility when the Hurst index H = 0.5, 0.6, 0.7, 0.8, 0.9 under L 2 and L 1/2 regularization.A common measure RMS E is used to evaluate the accuracy of reconstructed volatility σ(t) compared with exact volatility σ(t), and it is defined as The European call option prices generated by Lemma 1 and volatility function (4.1) are represented in Table 3. Figure 1 and Figure 2 show the calibration results.We can see that both schemes can satisfactorily recover the exact volatility function, and L 1/2 regularization is slightly more accurate than L 2 regularization.The results indicate that there is a long-term correlation from the volatility itself long term correlation when Hurst index is higher than 0.5, and it is consistent with the existing research on the recurrence interval in finance [8].

Stability test
To emphasize the stability of both methods, we introduce a random variable x which follows the standard normal distribution [29,30], and the volatility function with noise is Using the noised volatility, we can obtain the corresponding noised option prices in market, which will be employed to reconstruct volatility function with the proposed algorithms.The iteration stop criterion is max |σ d (t) − σ d−1 (t)| < 10 −7 and the maximum number of iterations is 200, where σ d (t) is the volatility function reconstructed in the dth iteration.Figure 3 compares the reconstructed and real volatility under L 2 and L 1/2 schemes when Hurst parameter is H = 0.7.Table 4 gives the comparison of stability results between L 2 and L 1/2 schemes.We can conclude that the performance of Hurst index H in the range of 0.6 to 0.8 is better than that of standard Brownian motion with H = 0.5, and L 1/2 regularization has slightly smaller deviation than L 2 scheme.And both methods are stable since the recovered time-dependent volatility resembles the exact function though adding a little disturbance.In order to measure the quality of the calibrated volatility by different Hurst parameters, we use root mean squared error (RMSE) between the recalculated option price by model and market price.The accuracy of calibrated volatility has reported in Table 6 and Figures 4-6.The results indicate that H = 0.6 has a much smaller RMSE than H = 0.5 and performs better for far-month options.
It is clear that the fractional model such as H = 0.6 performs better than the general geometric Brownian model where H = 0.5 due to the long-term dependence of fractional Brownian motion.In China's stock market dominated by individual investors, the high-frequency and short-term trading behavior is relatively common, which results in the short-term and medium-term factors occupying the dominant position in the market [35].Thus the predicted price corresponding to the longest maturity T = T 4 significantly differ from the actual market price, and we can see that the predicted price of far-month option is higher than the actual price for European call option.In addition, for short-term contracts of the current month T = T 1 , the pricing error of out-of-the-money option is higher than the in-the-money, which is mainly caused by extrapolation and fitting [36].Table 5.The SSE 50ETF option price at 9th September, 2021.

Conclusions
In this paper, we consider the Tikhonov L 2 and L 1/2 regularization for the construction of a timedependent volatility function from a finite set of market observations by the fractional Vasicek model.A fully implicit finite difference method is applied to solve the direct problem numerically.Several examples are given to demonstrate the accuracy and robustness of our proposed algorithms.The calibration results indicate that the Hurst index H ranging from 0.6 to 0.8 performs better than H = 0.5, and L 1/2 regularization has slightly smaller deviation but longer running time than L 2 scheme.

Figure 3 .
Figure 3.The calibration results by noised data.

Table 2 .
The numerical results of different Hurst index.

Table 4 .
The numerical results of stability test.

Table 6 .
The numerical results of different Hurst index.