Simultaneous Calibration of European Option Volatility and Fractional Order under the Time Fractional Vasicek Model

Abstract

In this paper, we recover the European option volatility function $\sigma \left(t\right)$ of the underlying asset and the fractional order $\alpha$ of the time fractional derivatives under the time fractional Vasicek model. To address the ill-posed nature of the inverse problem, we employ Tikhonov regularization. The Alternating Direction Multiplier Method (ADMM) is utilized for the simultaneous recovery of the parameter $\alpha$ and the volatility function $\sigma \left(t\right)$. In addition, the existence of a solution to the minimization problem has been demonstrated. Finally, the effectiveness of the proposed approach is verified through numerical simulation and empirical analysis.

1. Introduction

The option pricing problem is an important issue in the financial field. Black and Scholes [1] introduced the Black–Scholes (BS) formula, assuming a constant risk-free rate r and that the underlying asset S satisfies the following equation

$$\frac{dS}{S}=\mu dt+\sigma d{W}_t,$$

where $\mu$ is the expected return rate, $\sigma$ is the volatility, and $d{W}_t$ is the standard Brown motion. Under these conditions, the BS formula is derived as

$$\mathrm{Call} \mathrm{option} \mathrm{price}C=SN\left(d_1\right)-K{e}^{-r(T-t)}N\left(d_2\right),$$

$$\mathrm{Put} \mathrm{option} \mathrm{price}P=K{e}^{-r(T-t)}N(-d_2)-SN(-d_1),$$

where

$$d_1=\frac{ln(S/K)+(r+{\sigma }^2/2)(T-t)}{\sigma \sqrt{T-t}},$$

$$d_2=\frac{ln(S/K)+(r-{\sigma }^2/2)(T-t)}{\sigma \sqrt{T-t}},$$

$N(·)$ denotes the cumulative standard normal distribution function, and K is the exercise price. However, in real financial markets, interest rates change over time. Therefore, in the option pricing process, the impact of interest rate uncertainty on option prices must be considered. There are various interest rate models, including the Constant Elasticity of Variance (CEV) model [2], Cox–Ingersoll–Ross (CIR) model [3], Vasicek model [4], etc. Specifically, our focus in this article is on the Vasicek model.

Fractional derivatives are widely used in fractal theory [5], diffusion theory [6], signal processing [7], and financial theory ([8,9,10]) due to their non-locality advantages, i.e., the current state is influenced not only by the past instantaneous state but also by the state of the past period of time, which is more in line with the options market. Fractional derivative have two main forms: the Riemann–Liouville derivative [11] and the Caputo derivative [12]. In our discussion, we are mainly concerned with fractional derivatives in the Caputo sense. Since most partial differential equations lack analytic solutions, discretization of the equations in time and space is necessary. The commonly used methods for discrete Caputo time fractional derivatives include the Gr$\ddot{\mathrm{u}}$nwald–Letnikov (GL) formula and L1 formula. Specifically, the GL formula has a first-order accuracy in time direction after discretization, while the L1 formula has a (2-$\alpha$)-order accuracy in time direction. Subsequently, Gao et al. [13] introduced the L1-2 formula with (3-$\alpha$)-order accuracy, Alikhanov [14] proposed the L2-$1_{\sigma }$ formula with (3-$\alpha$)-order accuracy, and Cao et al. [15] improved the accuracy of the discretization to the (3-$\alpha$)-order Caputo derivative. In recent studies, Cao et al. [16] and Mokhtari et al. [17] improved the accuracy to (4-$\alpha$)-order. It is worth noting that the process of improving accuracy involves a higher-order interpolation of the function under the integral of the Caputo derivative. Using different interpolation functions results in different levels of precision.

Nourian et al. [8] solved the time fractional Black–Scholes equation by introducing a new set of wavelet functions. Roul [9] constructed a high-order computational scheme for European option pricing; this method is (3-$\alpha$)-order in the time direction. Chen et al. [10] introduced a novel operator splitting method to address American options within the framework of the time-fractional Black–Scholes model. Li et al. [18] presented a Newton-linearized Galerkin finite element method for solving time-oriented nonsmooth nonlinear time-fraction parabolic problems. The optimal error estimate of L2 norm is obtained without limitation of time step of spatial mesh size. The theoretical results are verified by numerical experiments. In [19], the convergence of a class of fully implicit time-stepping Galerkin finite element method for one-dimensional nonlinear subdiffusion equations is proved by using the generalized discrete Gronwall inequality. Based on the convergence of the direct time step scheme, a simple method to prove the convergence of the fast time-stepping Galerkin finite element method is presented. Yuan et al. [20] proposed a linearized fast high-order time-stepping scheme to solve the spatiotemporal fractional Schrodinger equation. A fast L2-$1_{\sigma }$ formula was used to approximate the time of non-uniform mesh. The Fourier spectrum method is used to discretize the space. The result of unconditional convergence is given. Zhang and Jiang [21] considered the two-dimensional nonlinear time-space fractional Schrodinger equation and applied the second-order fractional backward difference formula in the time direction and the Fourier spectrum method in the space direction to solve the model numerically. By using the generalized discrete Gronwall inequality and the space-time error splitting argument, the convergence of the fast time-step numerical method is simply proved without applying the Courant–Friedrichs–Lewy (CFL) condition.

Parameter calibration is a challenging problem in the financial field. Different from the forward problem of solving option prices with known model parameters, the parameter calibration problem refers to solving the model parameters based on the price data observed in the market. The volatility is one of the most important parameters in the option pricing, which is an indicator to measure the uncertainty of the return rate of assets and is used to reflect the risk level of financial assets. However, the volatility cannot be derived directly from market data and needs to be calibrated. This is a typical inverse problem, and it is usually ill-posed because the calibration problem may have no solution or have infinite solutions, or there may be solutions that are discontinuously dependent on observation data; many scholars have conducted research in this area. Dupire [22] proposed a local volatility model for option pricing in 1994, calibrated the volatility, and obtained an expression for volatility, which is called the Dupire formula. Xu and Jia [23] studied the Tikhonov regularization method for the volatility calibration problem under the jump diffusion model and used an iterative algorithm to solve the problem to obtain the volatility. Jiang et al. [24] reconstructed the implicit local volatility function using the optimal control framework and presented numerical experiments. Based on the Legendre pseudo-spectral method of space discretization and the finite difference method of time discretization, Li and Zhou [25] gave numerical solutions for the optimal control problem of time-fractional diffusion equations. Zhao and Xu [26] jointly recovered the mean-reversion parameter $\gamma$ and the volatility function under the fractional Chan–Karolyi–Longstaff–Sanders (CKLS) stochastic interest rate model. At the same time, numerical experiments are presented that show the fractional CKLS model ($\gamma$ = 0.8451) is more suitable for SSE 50ETF than the fractional Vasicek model. Yimamu and Deng [27] investigated the inverse problem of identifying volatility in option pricing. They accomplished this problem by transforming the parabolic equation defined on an unbounded region into a degenerate parabolic equation on a bounded region through a variable substitution. The study establishes that the optimal approximate solution converges effectively to the true solution of the original problem. In our paper, we use the alternating direction multiplier method to simultaneously reconstruct fractional order $\alpha$ and time-dependent volatility function $\sigma \left(t\right)$. Numerical examples and empirical analysis demonstrate the effectiveness of this method. Our research aims to reconstruct the volatility $\sigma$ and the fractional order $\alpha$ in the time-fractional Vasicek model. This unique perspective, to the best of our knowledge, has not been thoroughly explored in previous studies. We aim to provide readers with a novel perspective, offering a potential reference for future studies dealing with time-fractional direct problems, and offer a valuable contribution to the field by addressing the challenges associated with the simultaneous reconstruction of the volatility and the fractional order.

The remainder of this article is organized as follows. Section 2 gives the pricing formula of European put options under the time fractional Vasicek model. Section 3 gives the existence of the solution to the minimization problem and introduces the Tikhonov regularization method and ADMM algorithm for solving the volatility $\sigma \left(t\right)$ and the fractional order $\alpha$. Section 4 combines numerical examples and empirical analysis to obtain the stability of the algorithm. Finally, Section 5 provides conclusions.

2. Pricing Formula under the Time Fractional Vasicek Model

In this section, we present some fundamental definitions for our exploration and the pricing formula of European put options under the Caputo time fractional Vasicek model, laying the foundation for the subsequent exploration and reconstruction of the volatility $\sigma$ and the parameter $\alpha$.

Definition 1

(Kolmogorov [28]). Assume that the system state $Y_t$ satisfies the stochastic differential equation

$$d{Y}_t=a(Y_t,t)dt+b(Y_t,t)d{W}_t,$$

where W(t) is a Wiener process (also called Brownian motion), then the conditional mathematical expectation $u(y,t)=E\left(f\left(Y_T\right|Y_t=y)\right)$ is the solution of the Cauchy problem (terminal value problem) of the following backward parabolic equation

$$\frac{\partial u}{\partial t}+a(y,t)\frac{\partial u}{\partial y}+\frac{1}{2}{b}^2(y,t)\frac{{\partial }^{2}u}{\partial y^2}=0.$$

$$u(y,T)=f\left(y\right)(y\in \mathbb{R}).$$

Definition 2

(Feynman–Kac [28]). Let v be the conditional mathematical expectation

$$v(y,t)=E\left(f\left(Y_T\right)e^{{\int }_t^{T}g(Y_s,s)ds}|Y_t=y\right)=E_{y,t}\left(f\left(Y_T\right)e^{{\int }_t^{T}g(Y_s,s)ds}\right),$$

then it is the solution of the terminal value problem of the following backward parabolic equation

$$\frac{\partial v}{\partial t}+a(y,t)\frac{\partial v}{\partial y}+\frac{1}{2}{b}^2(y,t)\frac{{\partial }^{2}v}{\partial y^2}+g(y,t)v=0,(y\in \mathbb{R},t\in [0,T)),$$

$$v(y,T)=f\left(y\right)(y\in \mathbb{R}).$$

Definition 3.

The Caputo time fractional derivative is defined as follows

$${}_c{}^{0}D_t^{\alpha }V\left(t\right)=\left\{\begin{array}{cc}\hfill & \frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{V^{\prime }\left(t^{\prime }\right)}{{(t-t^{\prime })}^{\alpha }}d{t}^{\prime },0<\alpha <1,\hfill \\ \hfill & V^{\prime }\left(t\right),\alpha =1.\hfill \end{array}\right.$$

The pricing problem of the European option under the Vasicek model is considered. Suppose the asset price $S_t$ adheres to the following stochastic process

$$d{S}_t=r_t{S}_{t}dt+{\sigma }_t{S}_{t}d{W}_t^S$$

(1)

and the interest rate process $r_t$ satisfies

$$d{r}_t=(a+b{r}_t)dt+\omega d{W}_t^r,$$

(2)

where ${\sigma }_t$ is the volatility, $W_t^S$ and $W_t^r$ are two correlated Brown motions with $Cov(d{W}_t^s,d{W}_t^r)=\rho dt$ and $\rho \in (-1,1)$, and $a,b,\omega$ are constants.

Zero coupon bond is the carrier of interest rate, and we first consider the pricing of zero coupon bond. Let $P_t=P(r,t)$ represent the value of zero-coupon bonds at time t. When $r=r_t$ is a stochastic process, the value of zero coupon bonds $P_t=P(r_t,t)$ is

$$P_t=E\left(e^{-{\int }_t^{T}r\left(s\right)ds}|r\left(t\right)=r_t\right).$$

Therefore, from the Kolmogorov theorem and the Feynman–Kac formula, we can obtain the backward parabolic equation that the price of the zero coupon bond satisfies under the Vasicek model

$$\frac{\partial P}{\partial t}+\frac{{\sigma }^2}{2}\frac{{\partial }^{2}P}{\partial r^2}+(a+br)\frac{\partial P}{\partial r}-rP=0,(r\in \mathbb{R},t\in [0,T)),$$

(3)

$$P(r,T)=1.$$

Next, the pricing formula of the European put option under the Vasicek model is derived.

Lemma 1.

Assuming that the stochastic interest rate r and the underlying asset price S satisfy the Formulas (1) and (2), then the price of the European put option $V(r,S,t)$ with maturity T and strike K satisfies the following equation

$$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2\left(t\right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\rho \sigma \left(t\right)\omega S\frac{{\partial }^{2}V}{\partial S\partial r}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}+(a+br)\frac{\partial V}{\partial r}+rS\frac{\partial V}{\partial S}-rV=0,$$

(4)

with the condition

$$V(r,S,T)={(K-S)}^{+}.$$

Proof. 

Consider a portfolio $\Pi$, which consists of $(-{\Delta }_1)$ units of the underlying asset S, $(-{\Delta }_2)$ units of the zero coupon bond $P(r,t)$, and one option $V(S,r,t)$, then the value of the portfolio at time t is

$$\Pi =V-{\Delta }_{1}S-{\Delta }_{2}P.$$

Applying $It\widehat{o}$’s Lemma, we obtain

$$\begin{array}{cc}d\Pi \hfill & =\left(\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2\left(t\right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}+S\sigma \left(t\right)\omega \rho \frac{{\partial }^{2}V}{\partial S\partial r}\right)dt+\left(\frac{\partial V}{\partial S}-{\Delta }_1\right)dS\hfill \\ \hfill & +\left(\frac{\partial V}{\partial r}-{\Delta }_2\frac{\partial P}{\partial r}\right)dr+{\Delta }_2\left(\frac{\partial P}{\partial t}+\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}P}{\partial r^2}\right)dt.\hfill \end{array}$$

(5)

To eliminate risk in this portfolio, we take

$${\Delta }_1=\frac{\partial V}{\partial S},{\Delta }_2=\frac{\partial V/\partial r}{\partial P/\partial r},$$

substitute it into (5), then combine Formula (3) to obtain

$$d\Pi =\left(\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2\left(t\right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}+S\sigma \left(t\right)\omega \rho \frac{{\partial }^{2}V}{\partial S\partial r}\right)dt-\frac{\partial V/\partial r}{\partial P/\partial r}\left[rP-(a+br)\frac{\partial P}{\partial r}\right]dt.$$

The no-arbitrage principle gives us

$$\mathbb{E}(d\Pi )=r\Pi dt,$$

then

$$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2\left(t\right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\rho \sigma \left(t\right)S\frac{{\partial }^{2}V}{\partial S\partial r}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}+(a+br)\frac{\partial V}{\partial r}+rS\frac{\partial V}{\partial S}-rV=0.$$

   □

To ensure the completeness of the narrative, we first introduce the direct problem.

Replacing the term in Equation (4) that involves the integer-order derivative with the Caputo fractional-order derivative, denoted as $\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}$, where order $\alpha \in (0,1)$, we have

$$\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2\left(t\right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\rho \sigma \omega S\frac{{\partial }^{2}V}{\partial S\partial r}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}+(a+br)\frac{\partial V}{\partial r}+rS\frac{\partial V}{\partial S}-rV=0,$$

(6)

the fractional order derivative $\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}$ is defined as

$$\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}\frac{\partial }{\partial t}{\int }_t^T\frac{V(S,t^{\prime })-V(S,T)}{{(t^{\prime }-t)}^{\alpha }}d{t}^{\prime },0<\alpha <1.$$

Then, we invert the time variable as $\tau :=T-t$ and the relationship between the Caputo derivative ${}_c{}^{0}D_{\tau }^{\alpha }V$ and the time fractional derivative $\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}$ given by article [29], that is, $\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}={-}_0^c{D}_{\tau }^{\alpha }V$.

In this way, we obtain the pricing formula of European put options under the Caputo time fractional Vasicek model

$$\begin{array}{cc}{}_c{}^{0}D_{\tau }^{\alpha }V\hfill & =\frac{1}{2}{\sigma }^2\left(\tau \right)S^2\frac{{\partial }^{2}V}{\partial S^2}+\rho \sigma \left(\tau \right)\omega S\frac{{\partial }^{2}V}{\partial S\partial r}+\frac{1}{2}{\omega }^2\frac{{\partial }^{2}V}{\partial r^2}\hfill \\ & +(a+br)\frac{\partial V}{\partial r}+rS\frac{\partial V}{\partial S}-rV,\hspace{1em}\hspace{1em} \hfill \end{array}$$

(7)

with the conditions

$$V(r,S,0)=max(K-S,0),$$

$$\underset{S\to 0}{lim}V(r,S,\tau )=K,\underset{S\to +\infty }{lim}V(r,S,\tau )=0,$$

$$\underset{r\to 0}{lim}V(r,S,\tau )=\underset{r\to +\infty }{lim}V(r,S,\tau )=max(K-S,0).$$

Direct problem: For the given initial boundary problem, the price of the European put option can be precisely computed as long as the volatility and parameter $\alpha$ are known. Subsequently, we can use the implicit finite difference method to numerically address this problem.

Inverse problem: Let $V_{ij}$ denote the market price of an option with strike price $K_{ij}(j=1,2,\cdots ,M_i)$ and expiration date $T_i(i=1,2,\cdots ,N)$, where $K_1\le K_2\le \cdots \le K_{M_i}$ and $T_1\le T_2\le \cdots \le T_N$. The inverse problem is a calibration problem, which can be formulated as finding a time-dependent volatility function and parameter $\alpha$ such that the predicted price $V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))$ falls between the corresponding bid price $V_{ij}^b$ and ask price $V_{ij}^a$, that is

$$V_{ij}^b\le V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))\le V_{ij}^a.$$

From the above formula, the calibration problem is transformed to solve such an optimization problem:

Minimize the following error function associated with the given set of market option prices

$$\Gamma (\alpha ,\sigma )=\sum _{i=1}^N\sum _{j=1}^{M_i}{[V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))-\widehat{V_{ij}}]}^2,$$

(8)

where $\widehat{V_{ij}}=(V_{ij}^a+V_{ij}^b)/2$ is the mean of the bid and ask prices.

The calibration problem can be formulated as finding the parameters that minimize the error function $\Gamma (\alpha ,\sigma )$. However, the parameter $\alpha$ and the volatility function $\sigma \left(t\right)$ are discontinuously dependent on the market price, which means that small disturbances in the market price may lead to large changes in the error function $\Gamma (\alpha ,\sigma )$. This is an ill-posed problem, so we use regularization methods to solve it.

3. Regularization Method

In this section, we introduce the Tikhonov regularization [30,31] to make the problem well posed. We calibrate fractional order $\alpha$ and volatility function $\sigma \left(t\right)$ by minimizing the following optimization problem

$$\underset{\alpha ,\sigma }{min}{\left|\right|(\alpha ,\sigma )\left|\right|}_{*}^2+\lambda \sum _{i=1}^N\sum _{j=1}^{M_i}{[V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))-\widehat{V_{ij}}]}^2,$$

(9)

where ${\left|\right|(\alpha ,\sigma )\left|\right|}_{*}={\left(\right|\left|\alpha \right||}_2^2+\left|\right|{\sigma }^{\prime }{\left|\right|}_2^2{)}^{1/2}$, $\left|\right|·{\left|\right|}_2$ represents $L^2$-norm, $\lambda >0$ is the regularization parameter, and $\widehat{V_{ij}}$ is the real market option price.

3.1. Existence of Solutions to Optimization Problems

In the following theorem, we establish the existence of a solution to the minimization problem. Let $M=(0,1)\times [0,T]$ and

$${\Phi }_{\lambda }(\alpha ,\sigma ,\widehat{V})=\left|\right|(\alpha ,\sigma ){\left|\right|}_{*}^2+\lambda \sum _{i=1}^N\sum _{j=1}^{M_i}{[V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))-\widehat{V_{ij}}]}^2$$

and we define ${\varphi }_{\lambda ,\widehat{V}}(\alpha ,\sigma ):(\alpha ,\sigma )\mapsto {\Phi }_{\lambda }(\alpha ,\sigma ,\widehat{V})$, where $(\alpha ,\sigma )\in L^2\left(M\right)$.

Theorem 1.

Fix λ, $\widehat{V}$ and define $q=\underset{(\alpha ,\sigma )\in L^2\left(M\right)}{inf}{\varphi }_{\lambda ,\widehat{V}}(\alpha ,\sigma )$, assume that for any $z\in (q,+\infty )$, the set $Q_{\lambda ,\widehat{V}}\left(z\right)={\varphi }_{\lambda ,\widehat{V}}^{-1}(\alpha ,\sigma )\subset L^2\left(M\right)$ is bounded and ${\varphi }_{\lambda ,\widehat{V}}(\alpha ,\sigma )$ is weakly lower semicontinuous. Then, the minimization problem (9) has a solution belonging to $Q_{\lambda ,\widehat{V}}\left(z\right)$.

Proof. 

Because $z>q$, then there must be a minimization sequence $({\alpha }_n,{\sigma }_n)$. $L^2\left(M\right)$ is a Hilbert space, that is, it is a reflexive Banach space, so we can deduce that any bounded sequence in $L^2\left(M\right)$ has a weakly convergent subsequence. Therefore, $({\alpha }_n,{\sigma }_n)$ contains subsequence $({\alpha }_{n_k},{\sigma }_{n_k})$, which weakly converges to $(\overline{\alpha },\overline{\sigma })$. And because of the weakly lower semi-continuity of ${\varphi }_{\lambda ,\widehat{V}}(\alpha ,\sigma )$,

$${\varphi }_{\lambda ,\widehat{V}}(\overline{\alpha },\overline{\sigma })\le inf\underset{k\to \infty }{min}{\varphi }_{\lambda ,\widehat{V}}({\alpha }_{n_k},{\sigma }_{n_k})=q$$

and $Q_{\lambda ,\widehat{V}}\left(z\right)={\varphi }_{\lambda ,\widehat{V}}^{-1}(\alpha ,\sigma )$ is weakly closed. Therefore, $(\overline{\alpha },\overline{\sigma })$ is a solution of the minimization problem belonging to $Q_{\lambda ,\widehat{V}}\left(z\right)$.    □

3.2. ADMM Algorithm

Now, we consider the augmented Lagrangian function

$$\begin{array}{cc}\hfill min{l}_{\mu }(\alpha ,\sigma ):& ={\left|\right|(\alpha ,\sigma )\left|\right|}_{*}^2+\lambda \sum _{i=1}^N\sum _{j=1}^{M_i}{[V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))-\widehat{V_{ij}}]}^2\hfill \\ \hfill & +\sum _{i=1}^N\sum _{j=1}^{M_i}{\mu }_{ij}[V({\tau }_0,r_0,S_0,T_i,K_{ij},\alpha ,\sigma \left(\tau \right))-\widehat{V_{ij}}],\hfill \end{array}$$

where ${\mu }_{ij}>0$ is the Lagrange multiplier.

The ADMM algorithm consists of staring from $({\sigma }^0,{\alpha }^0,{\mu }_{ij}^0)$ to generate inductively a sequence $({\sigma }^k,{\alpha }^k,{\mu }_{ij}^k)$ as follows:

−

Step 1: minimization with repect to $\sigma$:

$${\sigma }^{k+1}:=\underset{\sigma }{argmin}{l}_{{\mu }^k}({\alpha }^k,\sigma ),$$

−

Step 2: minimization with repect to $\alpha$:

$${\alpha }^{k+1}:=\underset{\alpha }{argmin}{l}_{{\mu }^k}(\alpha ,{\sigma }^{k+1}),$$

−

Step 3: update the Lagrange multiplier:

$${\mu }_{ij}^{k+1}:={\mu }_{ij}^k+\beta [V({\tau }_0,r_0,S_0,T_i,K_{ij},{\alpha }^{k+1},{\sigma }^{k+1}\left(\tau \right))-\widehat{V_{ij}}],$$

where $\beta$ is the step size.

For step 1

Similar to Reference [31], introduce a ‘false’ time parameter $\theta$ and a function $\widehat{\sigma }(t,\theta )$. Starting with initial guess $\widehat{\sigma }(t,0)$, solve the following parabolic equation

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial \widehat{\sigma }}{\partial \theta }& =\frac{{\partial }^2\widehat{\sigma }}{\partial t^2}-\lambda \sum _{i=1}^N\sum _{j=1}^{M_i}\frac{\partial V}{\partial \sigma }(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right))[V(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right))-\widehat{V_{ij}}]\hfill \\ \hfill & -\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^{M_i}{\mu }_{ij}\frac{\partial V}{\partial \sigma }(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right)).\hfill \end{array}\end{array}$$

(10)

From the necessary conditions for the existence of extreme values, we can obtain that if $\widehat{\sigma }(t,\theta )$ tends to a steady-state solution as $\theta \to \infty$, then this solution will satisfy the Euler Lagrange equation for the functional $l_{\mu }$.

In order to solve the above equation, we need to calculate $\frac{\partial V}{\partial \sigma }(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right))$ by its definition, define $\delta (·)$ as the Dirac delta function, then, calculate the following expression

$$\frac{\partial V}{\partial \sigma }(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right))={[\frac{d}{dϵ}V(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma \left(\tau \right)+ϵ\delta )]}_{ϵ=0}.$$

Solving the variational derivative is similar to the solution approach for Equation (7).

First, we define an operator as follows

$$G_{\sigma }={-}_0^c{D}_{\tau }^{\alpha }+\frac{1}{2}\sigma \left(\tau \right)S^2\frac{{\partial }^2}{\partial S^2}+\rho \sigma \left(\tau \right)\omega S\frac{{\partial }^2}{\partial S\partial r}+\frac{1}{2}{\omega }^2\frac{{\partial }^2}{\partial r^2}+(a+br)\frac{\partial }{\partial r}+rS\frac{\partial }{\partial S}-rI,$$

where I is the identity operator.

Therefore, when the volatility function is disturbed, the expected pricing function satisfies the following equation

$$G_{\sigma +ϵ\delta }V(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma +ϵ\delta )=0.$$

(11)

Then, differential (11) with $ϵ$ and evaluating for $ϵ=0$, we have

$$\begin{array}{cc}\hfill G_{\sigma }\frac{\partial V}{\partial \sigma }(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma )& =-\sigma \left(\tau \right)S^2\frac{{\partial }^{2}V}{\partial S^2}(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma )\hfill \\ \hfill & -\rho \omega S\frac{{\partial }^{2}V}{\partial S\partial r}(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,\sigma ).\hfill \end{array}$$

The variational derivative adheres to both the homogeneous boundary condition and the initial condition.

Finally, we iteratively solve Equation (10)

$$\begin{array}{cc}\frac{{\sigma }_n^k-{\sigma }_n^{k-1}}{\Delta \theta }=\hfill & \frac{{\sigma }_{n+1}^k-2{\sigma }_n^k+{\sigma }_{n-1}^k}{{(\Delta t)}^2}-\frac{1}{2}{\displaystyle \sum _{i=1}^N\sum _{j=1}^{M_i}}{\mu }_{ij}^k\frac{\partial V}{\partial \sigma }({\tau }_n,r_0,S_0,T_i,K_{ij},{\alpha }^k,{\sigma }^k\left({\tau }_n\right))\hfill \\ \hfill & -\lambda {\displaystyle \sum _{i=1}^N\sum _{j=1}^{M_i}}\frac{\partial V}{\partial \sigma }({\tau }_n,r_0,S_0,T_i,K_{ij},{\alpha }^k,{\sigma }^k\left({\tau }_n\right))[V({\tau }_n,r_0,S_0,T_i,K_{ij},{\alpha }^k,{\sigma }^k\left({\tau }_n\right))-\widehat{V_{ij}}].\hfill \end{array}$$

For step 2

$$\begin{array}{cc}\hfill {\alpha }^{k+1}& =\underset{\alpha }{argmin}{\left\{\right|\left|\alpha \right||}_2^2+\lambda \sum _{i=1}^N\sum _{j=1}^{M_i}{[V(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,{\sigma }^{k+1})-\widehat{V_{ij}}]}^2\hfill \\ \hfill & +\sum _{i=1}^N\sum _{j=1}^{M_i}{\mu }_{ij}[V(\tau ,r_0,S_0,T_i,K_{ij},{\alpha }^k,{\sigma }^{k+1})-\widehat{V_{ij}}]\},\hfill \end{array}$$

this optimization problem can be tackled by employing the Particle Swarm Optimization (PSO) algorithm (Algorithm 1). Alternatively, a direct exploration of $\alpha$ from 0 to 1 can also be conducted to solve the problem.

Algorithm 1: Particle Swarm Optimization (PSO)

Equations for the velocity and the position updates:

$$\mathrm{Velocity}\mathrm{update}:v_i(t+1)=w·v_i\left(t\right)+c_1·r_1·({\mathbf{p}}_i-{\alpha }_i\left(t\right))+c_2·r_2·(\mathbf{pg}-{\alpha }_i\left(t\right)),$$

(12)

$$\mathrm{Position}\mathrm{update}:{\alpha }_i(t+1)={\alpha }_i\left(t\right)+v_i(t+1),$$

(13)

where w is the inertia weight, $c_1$ and $c_2$ are acceleration constants typically set to $0.2$, $r_1$ and $r_2$ are random values in the range [0, 1], ${\mathbf{p}}_i$ represents the historical best position of particle ${\alpha }_i$, and $\mathbf{pg}$ is the global best position among all particles.

For step 3

Directly update the Lagrange multiplier ${\mu }_{ij}$.

4. Numerical Experiments

4.1. Numerical Simulation

In this part, we study the effectiveness of the reconstruction algorithm for two examples where the volatility $\sigma \left(t\right)$ and the fractional order $\alpha$ are known. In our numerical experiments, we utilized the parameters mentioned in [32] for consistency; let $K=[10,12,14,16],T=[0.25,0.5,0.75,1],$ $[S_{min},S_{max}]=[0.25,40],[r_{min},r_{max}]=[0.002,1.2]$, the initial value $S_0=8,r_0=0.25$, and other parameters

$$a=0.001,\omega =0.3,b=-0.2,\rho =0.4.$$

Example 1.

Considering volatility function $\sigma \left(\tau \right)=0.2-0.1log(1.5+3\tau )$, $\alpha =0.7$, the Lagrangian multiplier and the update step size of ${\mu }_{ij}$ are respectively taken as $\lambda =5,{\mu }_{ij}^0=0.3,\beta =5$.

We will use the European put option price calculated using the exact volatility $\sigma \left(\tau \right)$ and parameter $\alpha$ as the market price of the option. Stop iterating when $|{\sigma }_d\left(\tau \right)-{\sigma }_{d-1}\left(\tau \right)|<1\times {10}^{-6}$, where ${\sigma }_d\left(\tau \right)$ represents the volatility generated by the d-th iteration, and the total number of iterations is not greater than 100. Figure 1 presents a graph of exact volatility versus reconstructed volatility, and the recover value of parameter $\alpha$ is 0.7. We add noise with intensity $\delta$ of 0.01, 0.03, and 0.05 to the exact volatility to verify the stability

$${\sigma }^{\delta }=0.2-0.1log(1.5+3\tau )+\delta \times rand(-1,1).$$

We first use perturbed volatility to calculate option prices and then use the resulting option prices to reconstruct fractional $\alpha$ and volatility $\sigma \left(\tau \right)$. Figure 2 shows the comparison of real volatility and reconstructed volatility under different noise intensities.

Table 1. Numerical results under different intensity disturbances.

	$max|{\mathit{\sigma }}_{\mathit{d}}-{\mathit{\sigma }}_{\mathit{d}-1}|$	RMSE	$max|\mathit{\sigma }-\widehat{\mathit{\sigma }}|$	$\mathit{\alpha }$	$max|\mathit{V}-\widehat{\mathit{V}}|$
$\delta =0$	9.4106 × ${10}^{-7}$	9.4772 × ${10}^{-4}$	1.836 × ${10}^{-3}$	0.700	1.07 × ${10}^{-3}$
$\delta =0.01$	9.8816 × ${10}^{-7}$	0.0026	6.565 × ${10}^{-3}$	0.700	1.249 × ${10}^{-3}$
$\delta =0.03$	9.9998 × ${10}^{-7}$	0.0052	1.9694 × ${10}^{-2}$	0.697	9.465 × ${10}^{-3}$
$\delta =0.05$	1.0000 × ${10}^{-6}$	0.0090	3.2824 × ${10}^{-2}$	0.696	1.309 × ${10}^{-2}$

gives the root mean square error (RMSE), the maximum volatility error, the reconstructed fractional $\alpha$, and the maximum option price error, where RMSE can be calculated by

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}\sum _{i=1}^n{({\sigma }_i-\widehat{{\sigma }_i})}^2}.$$

Example 2.

In this example, we examine the volatility function represented as $\sigma \left(\tau \right)=0.5{(\tau -0.5)}^2+0.1$, $\alpha =0.37$; the Lagrangian multiplier and the update step size of ${\mu }_{ij}$ are respectively taken as

$$\lambda =1.5,{\mu }_{ij}^0=0.2,\beta =1,$$

the volatility curve exhibits a smile-shaped pattern.

Figure 3 shows an image of exact volatility versus reconstructed volatility. We can obtain the precise volatility and reconstructed volatility under different noise intensities from Figure 4. Finally,

Table 2. Numerical results under different intensity disturbances.

	$max|{\mathit{\sigma }}_{\mathit{d}}-{\mathit{\sigma }}_{\mathit{d}-1}|$	$\mathbf{RMSE}$	$max|\mathit{\sigma }-\widehat{\mathit{\sigma }}|$	$\mathit{\alpha }$	$max|\mathit{V}-\widehat{\mathit{V}}|$
$\delta =0$	9.9770 × ${10}^{-7}$	0.0030	4.621 × ${10}^{-3}$	0.364	2.855 × ${10}^{-3}$
$\delta =0.01$	9.9975 × ${10}^{-7}$	0.0044	7.820 × ${10}^{-3}$	0.367	5.162 × ${10}^{-3}$
$\delta =0.03$	9.9999 × ${10}^{-7}$	0.0073	1.9694 × ${10}^{-2}$	0.372	1.9478 × ${10}^{-2}$
$\delta =0.05$	4.8906 × ${10}^{-7}$	0.0156	3.2824 × ${10}^{-2}$	0.359	3.4620 × ${10}^{-2}$

shows the RMSE, the absolute error of volatility, and the reconstructed fractional $\alpha$ under different noise intensities.

4.2. Empirical Analysis

In this part, we utilize the proposed algorithm to recover the volatility function and fractional order based on the Shanghai Stock Exchange (SSE) 50ETF data on 12 December 2023.

Table 3. The option price of the SSE 50ETF on 12 December 2023.

K	$\mathit{V}(\mathit{K},{\mathit{T}}_1)$	$\mathit{V}(\mathit{K},{\mathit{T}}_2)$	$\mathit{V}(\mathit{K},{\mathit{T}}_3)$	$\mathit{V}(\mathit{K},{\mathit{T}}_4)$
$2.20$	0.0021	0.0099	0.0255	0.0411
$2.25$	0.0060	0.0189	0.0391	0.0569
$2.30$	0.0177	0.0353	0.0578	0.0759
$2.35$	0.0421	0.0588	0.0816	0.0986
$2.40$	0.0796	0.0902	0.1107	0.1241
$2.45$	0.1241	0.1290	0.1442	0.1542
$2.50$	0.1712	0.1714	0.1827	0.1885
$2.55$	0.2206	0.2182	0.2238	0.2256
$2.60$	0.2695	0.2651	0.2668	0.2648
$2.65$	0.3191	0.3145	0.3140	0.3092

shows real market price for SSE 50ETF put option price with respect to the maturity and strike. The maturity times are $T_1=11/250,T_2=31/250,T_3=76/250,$ $T_4=141/250$ and the strike prices are $K_i=2.20+(i-1)0.05$ for i = 1, 2, ..., 10. The current value of the underlying asset is $S_0=2.337$, and the interest rate $r_0=2.43\%$. The Lagrangian multiplier and the update step size of ${\mu }_{ij}$ are respectively chosen as

$$\lambda =20,{\mu }_{ij}^0=0.2,\beta =10$$

and other parameters are taken as:

$$a=0.004,\omega =0.1,b=-0.2,\rho =0.8.$$

Figure 5 and Figure 6 depict the comparison between the actual SSE 50ETF market prices and the reconstructed option prices at various time points, complementing the visual analysis, and

Table 4. The error outcomes following the calibration of SSE 50ETF options.

	$max|{\mathit{\sigma }}_{\mathit{d}}-{\mathit{\sigma }}_{\mathit{d}-1}|$	$\mathit{\alpha }$	$max|\mathit{V}-\widehat{\mathit{V}}|$
$\mathrm{Empirical}\mathrm{results}$	9.9351 × ${10}^{-7}$	0.120	1.810 × ${10}^{-2}$

provides information on multiple parameters obtained from the experimental results, and we can obtain $\alpha =0.120$, consistent with the recent trends in the SSE 50ETF.

5. Conclusions

In this article, we reconstruct both the volatility $\sigma$ and the fractional order $\alpha$ under the time fractional Vasicke model. Since the calibration problem is ill-posed, we used Tikhonov regularization and the ADMM algorithm for an iterative solution. At the same time, the existence of the solution to the minimization problem is given. And the effectiveness and robustness of the method can be obtained from numerical examples and empirical analysis. From empirical analysis, it can be concluded that $\alpha =0.120$ is more suitable for the current SSE 50ETF. Furthermore, it serves as an inspiration for readers to employ the derived fractional order in solving time-fractional forward problems, given its relevance to the current Chinese options market. It is worth noting that the possibility of negative interest rates could be addressed by considering alternative models such as the CIR model or CEV model.