Next Article in Journal
Theme-Aware Semi-Supervised Image Aesthetic Quality Assessment
Next Article in Special Issue
A Nonlinear Multigrid Method for the Parameter Identification Problem of Partial Differential Equations with Constraints
Previous Article in Journal
Composite Indicator of the Organisational Information and Communication Technologies Infrastructure—A Novel Statistical Index Tool
Previous Article in Special Issue
Feature Reconstruction from Incomplete Tomographic Data without Detour
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Convergence of Inverse Volatility Problem Based on Degenerate Parabolic Equation

Department of Mathematics, Lanzhou Jiaotong University, Lanzhou 730070, China
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(15), 2608; https://doi.org/10.3390/math10152608
Submission received: 12 June 2022 / Revised: 21 July 2022 / Accepted: 22 July 2022 / Published: 26 July 2022
(This article belongs to the Special Issue Inverse Problems and Imaging: Theory and Applications)

Abstract

:
Based on the theoretical framework of the Black–Scholes model, the convergence of the inverse volatility problem based on the degenerate parabolic equation is studied. Being different from other inverse volatility problems in classical parabolic equations, we introduce some variable substitutions to convert the original problem into an inverse principal coefficient problem in a degenerate parabolic equation on a bounded area, from which an unknown volatility can be recovered and deficiencies caused by artificial truncation can be solved. Based on the optimal control framework, the problem is transformed into an optimization problem and the existence of the minimizer is established, and a rigorous mathematical proof is given for the convergence of the optimal solution. In the end, the gradient-type iteration method is applied to obtain the numerical solution of the inverse problem, and some numerical experiments are performed.

1. Introduction

With the continuous development of science and technology, a large number of inverse problems [1,2,3] have appeared in the fields of weather forecasting, material nondestructive testing, wave inverse scattering, biological imaging and option pricing. The inverse volatility problem has received much attention in recent years due to its significant impacts on the market value of options (see, e.g., [4,5,6,7,8,9,10,11]). In the real market, the volatility of the underlying asset price cannot be directly observed, but the market price of an option can be directly observed and can provide information about the volatility of the underlying asset price. The associated forward operators are based on solutions to the corresponding Black–Scholes/Dupire equations (see [12], and Section 2 of this paper). In this paper, we investigate the following problem:
Problem P: In terms of a stock option without paying dividends, it is matter of public knowledge that V ( s , t ; K , T ) for a call option satisfies the following Black–Scholes equation
V t + 1 2 σ 2 ( s ) s 2 2 V s 2 + s μ V s r V = 0 , ( s , t ) R + × ( 0 , T ) , V ( s , t ) = ( s K ) + = max ( 0 , s K ) , s R + .
Herein, the parameters s , K are the price of the underlying stock and the strike price. T , μ and r are, respectively, the date of expiry, the risk-neutral drift and the risk-free interest rate, which are assumed to be constants. The parameter σ ( s ) represents the volatility coefficient to be identified.
Given the following additional condition:
V ( s * , t * , K , T ) = V * ( K , T ) , K R + ,
in which s * and V * ( K , T ) indicate market price of the stock at time t * and market price of the option with strike K at a given expiry time T. The inverse problem is to identify a pair of functions ( u , a ) satisfying (1) and (2), respectively.
The problem P was first considered by Dupire in [12]. In his paper, Dupire applied the symmetric property of the transition probability density function to replace the option pricing inverse problem with an equation that contains parameters K, T and Dupire’s proposed formula for calculating the implied volatility. Although Dupire’s formula is seriously ill-posed, the solution of Dupire provided the foundation for later understanding of this problem. Time-dependent and space-related volatility has been studied in [4,5]. For space-dependent volatility, authors generate a broad class of non-Gaussian stochastic processes. By using the dual equation, the problem is transformed into a inverse coefficient problem with final observations and some new uniqueness and stability theorems are established. In [6,7,8], an optimal control framework is applied to determine implied volatility in option pricing. A rigorous mathematical analysis about the existence, necessary condition and convergence of the inverse problem is made. The authors propose a well-posed approximation algorithm to identify the coefficients. In recent years, the linearization method has been applied to the inverse volatility problem in option pricing. In [13,14,15], the linearization method is applied to transform the problem P into an inverse source problem in classical parabolic equation, from which an unknown volatility can be recovered. However, there are several deficiencies in these studies that need to be improved. One of the significant deficiencies is the assumption that the original problem P is in the unbounded region, and it is under this assumption many scholars conducted numerical simulations by applying the artificial truncation. There is a potential issue in this approach; that is, if the truncation interval is too large, it will increase the amount of calculation, and if the truncation interval is too small, it will increase the error. In practical applications, this approach may fail to control the volatility risk. There are a number of papers devoted to this problem, see, e.g., [16,17,18,19].
The purpose of this paper is to better control the volatility risk of financial markets with precision and to address the deficiencies caused by the artificial truncation in previous studies performed by other scholars. In order to achieve our objective, in this paper we transform the principal coefficient inversion problem P of classical parabolic equation into the principal coefficient inversion problem of the degenerate parabolic equation in bounded regions via the variable substitutions. For the study of the inverse problems of degenerate parabolic equation, we refer the reader to [20,21,22], and for other topics of degenerate parabolic equation, such as null controllability, see references [23,24,25]. In the present paper, we use the optimal control framework to provide a full analysis of the convergence for the degenerate problem. To the best of our knowledge, there are no other works investigating the theoretical studies on the convergence of inverse volatility problem for a degenerate parabolic equation.
The structure of this paper is arranged as follows: In Section 2, we propose inverse the principal coefficient problem of degenerate parabolic equations through some variable substitutions, and use the optimal control framework to transform problem P1 into an optimal control problem P2. In Section 3, we prove the existence of the minimum of the optimal control problem. We show that the approximate solution of the optimal control problem converges to the true solution of the original problem in Section 4. In Section 5, we use the gradient iteration method to perform some numerical experiments.

2. Problem Setting and Methodology

In this section, we introduce some variables substitutions to transform the problem P into a degenerate parabolic principal term coefficient inversion problem P1. First, we apply the Dupire’s solution to the problem (1) by setting
2 V K 2 = G ( s , t ; K , T ) ,
it is well-known that G satisfies the following equation (see Appendix A for detailed derivation):
G t + 1 2 σ 2 ( s ) s 2 2 G s 2 + r s G s r G = 0 , 0 s < , 0 t < T , G ( s , t ; K , T ) = δ ( K s ) , 0 s < .
Since δ ( x ) = δ ( x ) , it is easy to see that G ( s , t ; K , T ) , as a function of K , T is the fundamental solution of the conjugate problem of the above problem (see [12]), namely
G T + 1 2 2 G K 2 ( σ 2 ( K ) K 2 G ) r G K ( K G ) r G = 0 , 0 K < , 0 t < T , G ( s , t ; K , T ) = δ ( K s ) , 0 K < .
Substitute G ( s , t ; K , T ) into the above conjugate problem and integrate K twice on [ K , ] , if s is fixed, as K ,
V , K V K , σ 2 K 2 G , K G K , ( σ 2 K 2 G ) K 0 .
Therefore, the problem (1) is transformed into the Cauchy problem of the following Dupire equation (see [4,5,6,12]):
V T + 1 2 σ 2 ( K ) K 2 2 V K 2 r K V K = 0 , ( T , K ) ( t , ) × ( 0 , ) , V ( s , t ; K , T ) | T = 0 = ( s 0 K ) + , K R + .
From above equation, Dupire obtained the formulas for local volatility under the different strike price and expiry date, i.e.,
σ ( K ) = V T + r K V K 1 2 K 2 2 V K 2 .
However, this formula is extremely sensitive to changes in the data, and small changes in the options market quotes may cause great changes in the second-order partial derivative, so the formula is unstable.
For such a case, we take
x = K K + E , x ( 0 , 1 ) ; u ( x , τ ) = V ( s 0 , 0 , K , T ) K + E ; τ = T ,
where E is a non-zero constant, u satisfies the following problem P1 (see Appendix A for detailed derivation).
Poblem P1: Consider the following degenerate parabolic equation
u τ 1 2 σ 2 ( x ) x 2 ( 1 x ) 2 2 u x 2 + r x ( 1 x ) u x + r x u = 0 , ( x , τ ) Q = ( 0 , 1 ) × ( 0 , T ] , u ( x , 0 ) = ( s 0 x E 1 x ) + , 0 < x < 1 ,
with the additional data
u ( x , τ ) = w ( x , τ ) , 0 < x < 1 ,   T p τ T .
where p ( 0 , T ) is a relatively small constant. The goal is to determine the unknown space-dependent volatility coefficient σ ( x ) and u satisfying the (4) and (5). For simplicity, we denote
a ( x ) = 1 2 σ 2 ( x ) x 2 ( 1 x ) 2 ,
b ( x ) = r x ( 1 x ) ,
c ( x ) = r x ,
ϕ ( x ) = ( s 0 x E 1 x ) + .
where ϕ ( x ) ,   b ( x ) , and c ( x ) are the smooth functions given on ( 0 , 1 ) , and a ( x ) is the unknown term satisfying
a ( 0 ) = a ( 1 ) = 0 , a ( x ) > 0 ,   x ( 0 , 1 ) .
Equation (4) are the degenerate parabolic equations. Compared with the classical parabolic equations, the boundary conditions of degenerate equations may vanish on some degenerate part. In [3], whether or not boundary conditions should be given is determined by the sign of the Fichera’s function. For general degenerate parabolic equations with non-divergent form, the Fichera function can be defined as follows:
B ( x , t ) = x k ( x , t ) b ( x , t ) n i , i = 1 , 2 ,
where k ( x , t ) satisfies k ( 0 , t ) = k ( l , t ) = 0 , k ( x , t ) > 0 , ( x , t ) ( 0 , l ) × ( 0 , T ] , and n i are taken as 1, −1 when x = 0 , x = l , respectively. Specifically, if the sign of the Fichera function is non-negative, then boundary condition should not be given. For Equation (4), it can be easily seen that B ( x , t ) | x = 0 , 1 = 0 . Therefore, the boundary conditions on x = 0 and x = 1 should not be given either.

3. Optimal Control Framework

Before we consider the regularization of problem P1, we discuss the well-posedness of direct problem and give some basic definitions, lemmas and estimations.
Definition 1.
Define H a 1 ( 0 , 1 ) to be the closure of C 0 ( 0 , 1 ) under the following norm:
u H a 1 ( 0 , 1 ) 2 = 0 1 a ( x ) | u | 2 + | u | 2 d x , u H a 1 ( 0 , 1 ) .
Definition 2.
A function u H 1 [ 0 , T ] ; L 2 ( 0 , 1 ) L 2 ( 0 , T ) ; H a 1 ( 0 , 1 ) is called the weak solution of (4), if u satisfies
u ( x , 0 ) = ϕ ( x ) , x ( 0 , 1 ) ,
and the following integration identity
0 1 u τ ψ d x + 0 1 a u x ψ x d x + 0 1 a x u x ψ d x + 0 1 b u x ψ d x + 0 1 c u ψ d x = 0 ,
holds for a.e., τ ( 0 , T ] .
Theorem 1.
For any given f L ( Q ) , φ L ( 0 , 1 ) , there exists a unique weak solution to (4) and satisfies the following estimate:
max 0 τ T 0 1 u 2 d x + 0 T 0 1 a u 2 d x d τ C 0 1 ϕ 2 d x .
Proof. 
Firstly, multiplying u on both sides of (4) and integrating over Q τ = ( 0 , 1 ) × ( 0 , τ ) , we have
0 τ 0 1 u τ u d x d τ 0 τ 0 1 a ( x ) u x x u d x d τ + 0 τ 0 1 b ( x ) u x u d x d τ + 0 τ 0 1 c ( x ) u 2 d x d τ = 0 .
Integration by parts, we obtain
1 2 0 1 u 2 d x 0 τ 0 1 ( b a x ) u u x d x d τ + 0 τ 0 1 a u x 2 d x d τ + 0 τ 0 1 c ( x ) u 2 d x d τ = 1 2 0 1 ϕ 2 d x .
Integration the second term from the left side of the equality (9) by parts, according to the Fichera condition and the Grönwall’s inequality, one has
max 0 τ T 0 1 u 2 d x + 0 T 0 1 a u 2 d x d τ C 0 1 ϕ 2 d x ,
where C is independent of T. □
Since the problem P1 is ill-posed, that is, the solution of the problem does not continuously depend on the observed data. Therefore, we use the optimal control framework to convert the problem P1 into the following optimal control problem P2:
Problem P2: Find a ¯ ( x ) A , such that
J ( a ¯ ) = min a A J ( a ) ,
where
J ( a ) = 1 2 p T p T 0 1 u ( x , τ ; a ) w ( x , τ ) 2 d x d τ + α 2 0 1 a 2 d x ,
A = a ( x ) | a ( 0 ) = a ( 1 ) = 0 ,   a ( x ) > 0 ,   a ( x ) H 1 ( 0 , 1 ) H 0 1 ( 0 , 1 ) ,
u ( x , τ ; a ) is the solution of (4) for a given coefficient a ( x ) A , α is the regularization parameter.
We are now going to show the existence of minimizers to the problem (10). Firstly, we assert that the functional J ( a ) is of some continuous property in A with the following sense.
Lemma 1.
For any sequence { a n } in A which converges to some { a } A in L 1 ( 0 , 1 ) as n , we have
lim n 0 1 u ( x , τ ; a n ) w ( x , τ ) 2 d x = 0 1 u ( x , τ ; a ) w ( x , τ ) 2 d x .
Proof. 
Step1: Setting a = a n and selecting the test function as u ( a n ) ( · , τ ) in (7), one can derive the following equality by integrating on τ .
1 2 u ( x , τ ; a n ) L 2 ( 0 , 1 ) 2 + 0 τ 0 1 a u ( x , τ ; a n ) 2 d x d τ + 0 τ 0 1 ( b + a x ) u ( x , τ ; a n ) x u ( x , τ ; a n ) d x d τ + 0 τ 0 1 c u 2 ( x , τ ; a n ) d x d τ = 1 2 ϕ L 2 ( 0 , 1 ) 2 ,
for any τ ( 0 , T ] . By the Fichera condition, Young’s inequality and the smoothness of a , b , we obtain
u ( x , τ ; a n ) L 2 ( 0 , 1 ) 2 + 0 τ 0 1 a u ( x , τ ; a n ) 2 d x d τ C ϕ L 2 ( 0 , 1 ) 2 .
Since the sequence { u ( a n ) } is uniformly bounded in the space L 2 ( 0 , T ) ; H a 1 ( 0 , 1 ) , we can extract a subsequence denoted as { u ( a n ) } , such that
u ( x , τ ; a n ) u * ( x , τ ) L 2 ( ( 0 , T ) ; H a 1 ( 0 , 1 ) ) L 2 ( ( 0 , T ) ; L 2 ( 0 , 1 ) ) .
Step 2: u * ( x , τ ) = u ( a ) ( x , τ ) .
By taking a = a n in (7), multiplying both sides by a function η ( τ ) C 1 [ 0 , T ] with η ( T ) = 0 , and integrating with respect to τ , we have
0 T 0 1 u ( a n ) τ ψ η ( τ ) d x d τ + 0 T 0 1 a u ( a n ) x ψ x η ( τ ) d x d τ + 0 T 0 1 ( a x + b ) u ( a n ) x ψ η ( τ ) d x d τ + 0 T 0 1 c u ( a n ) ψ η ( τ ) d x d τ = 0 .
Integration by part for the first term on the left side of (17), letting n , in (17), and using (16), we obtain
0 1 ϕ ψ η ( 0 ) d x 0 T 0 1 u * ψ η ( τ ) τ d x d τ + 0 T 0 1 a u x * ψ x η ( τ ) d x d τ + 0 T 0 1 ( a x + b ) u x * ψ η ( τ ) d x d τ + 0 T 0 1 c u * ψ η ( τ ) d x d τ = 0 .
Noticing that (18) is valid for any η ( τ ) C 1 [ 0 , T ] with η ( T ) = 0 , we have
0 1 u τ * ψ d x + 0 1 a u x * ψ x d x + 0 1 a x u x * ψ d x + 0 1 b u x * ψ d x + 0 T 0 1 c u * ψ d x d τ = 0 , ψ L 2 ( 0 , 1 ) H a 1 ( 0 , 1 ) ,
and u * ( x , 0 ) = ϕ ( x ) . Therefore, u * = u ( x , τ ; a ) by the definition of u ( a ) .
Step 3: Prove u ( x , τ ; a n ) w ( x , τ ) L 2 ( 0 , 1 ) u ( x , τ ; a ) w ( x , τ ) L 2 ( 0 , 1 ) as n .
We rewrite (7) for a = a n in the form
0 1 u ( a n ) w τ ψ d x + 0 1 a n u ( a n ) w ψ d x + 0 1 b + ( a n ) x u ( a n ) w ψ d x + 0 1 c ( u ( a n ) w ) ψ d x = 0 1 a n w ψ d x 0 1 b + ( a n ) x w ψ d x .
Taking ψ = u ( a n ) w in (20), we have
1 2 d d τ u ( a n ) w L 2 ( 0 , 1 ) 2 + 0 1 a n u ( a n ) w 2 d x + 0 1 1 2 b + ( a n ) x x + c u ( a n ) w 2 d x = 0 1 a n w u ( a n ) w d x 0 1 b + ( a n ) x w u ( a n ) w d x .
Similar relations hold for u ( a ) , namely
1 2 d d τ u ( a ) w L 2 ( 0 , 1 ) 2 + 0 1 a u ( a ) w 2 d x + 0 1 1 2 b + a x x + c u ( a ) w 2 d x = 0 1 a w u ( a ) w d x 0 1 b + a x w u ( a ) w d x .
Subtracting (22) from (21), one has
1 2 d d τ u ( a n ) u ( a ) L 2 ( 0 , 1 ) 2 + 0 1 a n | u ( a n ) w | 2 a | u ( a ) w | 2 d x + 0 1 1 2 b + a x x u ( a ) w 2 1 2 b + ( a n ) x x u ( a n ) w 2 d x = 0 1 a w u ( a ) w a n w u ( a n ) w d x + 0 1 b + a x w u ( a ) w b + ( a n ) x w u ( a n ) w d x 0 1 d d τ u ( a ) w u ( a n ) u ( a ) d x .
Taking ψ = u ( a n ) u ( a ) in (7), we have
0 1 u ( a ) τ u ( a n ) u ( a ) d x = 0 1 a u ( a ) u ( a n ) u ( a ) d x 0 1 ( a x + b ) u ( a ) u ( a n ) u ( a ) d x .
Similarly, by taking ψ = u ( a ) w we have
0 1 u ( a n ) u ( a ) τ u ( a ) w d x = 0 1 a u ( a n ) u ( a ) u ( a ) w d x 0 1 ( a x + b ) u ( a n ) u ( a ) u ( a ) w d x .
Substituting (24) and (25) into (23), after some manipulations, we derive
1 2 d d τ u ( a n ) u ( a ) L 2 ( 0 , 1 ) 2 + 0 1 a n | u ( a n ) w | 2 a | u ( a ) w | 2 d x + 0 1 1 2 b + a x x u ( a ) w 2 1 2 b + ( a n ) x x u ( a n ) w 2 d x + 0 1 a n w u ( a n ) w a w u ( a ) w d x + 0 1 b + ( a n ) x w u ( a n ) w b + a x w u ( a ) w d x = 0 1 a u ( a ) u ( a n ) u ( a ) d x + 0 1 ( a x + b ) u ( a ) u ( a n ) u ( a ) d x + 0 1 a u ( a n ) u ( a ) u ( a ) w d x + 0 1 ( a x + b ) u ( a n ) u ( a ) u ( a ) w d x = A n + B n .
Integrating over the interval ( 0 , τ ) for any τ T , we obtain
1 2 u ( a n ) u ( a ) L 2 ( 0 , 1 ) 2 0 T A n + B n d τ .
By the convergence of a n and weak convergence of u ( a n ) , one can easily obtain
0 T A n + B n d τ 0 , a s n .
Combining (27) and (28) we have
max τ [ 0 , T ] u ( x , τ ; a n ) u ( x , τ ; a ) L 2 ( 0 , 1 ) 2 0 , a s n .
Besides, by the H ÷ l d e r inequality, there holds
0 1 u ( x , τ ; a n ) w ( x , τ ) 2 d x 0 1 u ( x , τ ; a ) w ( x , τ ) 2 d x 0 1 u ( x , τ ; a n ) u ( x , τ ; a ) · u ( x , τ ; a n ) + u ( x , τ ; a ) 2 w d x u ( x , τ ; a n ) u ( x , τ ; a ) L 2 ( 0 , 1 ) · u ( x , τ ; a n ) + u ( x , τ ; a ) 2 w L 2 ( 0 , 1 ) .
From (16), (29) and (30) we obtain
lim n 0 1 u ( x , τ ; a n ) w ( x , τ ) 2 d x = 0 1 u ( x , τ ; a ) w ( x , τ ) 2 d x .
This completes the proof of Lemma 1. □
Theorem 2.
There exists a minimizer a ¯ A of J ( a ) , i.e.,
J ( a ¯ ) = min a A J ( a ) .
Proof. 
It is obvious that J ( a ) is non-negative and thus J ( a ) has the greatest lower bound inf a A J ( a ) . Let { a n } be a minimizing sequence, i.e.,
inf a A J ( a ) J ( a n ) inf a A J ( a ) + 1 n , n = 1 , 2 , .
Notice that J ( a n ) C , we deduce
a n L 2 ( 0 , 1 ) C ,
where C is independent of n. From the boundedness of { a n } and (31), we also have
a n H 1 ( 0 , 1 ) C .
So we can extract a subsequence, still denoted by { a n } , such that
a n ( x ) a ¯ ( x ) H 1 ( 0 , 1 ) a s n .
By using the Sobolev’s imbedding theorem (see [1]) we obtain
a n ( x ) a ¯ ( x ) L 1 ( 0 , 1 ) 0 a s n .
It can be easily seen that { a n ( x ) } A . So we get as n that
a n ( x ) a ¯ ( x ) A ,
in L 1 ( 0 , 1 ) .
Moreover, from (35) we have
0 l a ¯ 2 d x = lim n 0 l a n · a ¯ d x lim n 0 l a n 2 d x · 0 l a ¯ 2 d x .
From Lemma 1 and the convergence of { a n } , there exists a subsequence of { a n } , still denoted by { a n } , such that
lim n 0 1 u ( x , τ ; a n ) w ( x , τ ) 2 d x = 0 1 u ( x , τ ; a ¯ ) w ( x , τ ) 2 d x .
By (35), (36) and (37), we obtain
J ( a ¯ ) = lim n J ( a n ) = min a A J ( a ) ,
Hence,
J ( a ¯ ) = min a A J ( a ) .

4. Convergence

In this section, we consider whether the solution of the control problem (10) converges to the solution of the original problem (4). What needs to be declared is that the optimal control problem P2 is non-convex, so a unique solution may not be obtained. In fact, the optimization technique is a classic tool for obtaining “general solutions” for inverse problems that do not have a unique solution. However, we can obtain the necessary conditions for the optimal solution, and through the necessary conditions it can be proved that the optimal control problem has local uniqueness when T is relatively small. We do not elaborate here because the techniques are similar to [6]. We are more interested in its convergence. We assume the true solution a * ( x ) is achievable, namely, there exists a * ( x ) H 1 ( 0 , 1 ) H 0 1 ( 0 , 1 ) , such that
u ( x , t ; a * ) = w ( x , τ ) .
Observation data tends to have errors, so it can be recorded as w δ , assuming that the error level satisfies:
w w δ L 2 ( ( 0 , 1 ) × ( T p , T ) ) p δ .
We define a forward operator u ( a ) :
u ( a ) : A H 1 ( ( 0 , T ) ; L 2 ( 0 , 1 ) ) L 2 ( ( 0 , T ) ; H a 1 ( 0 , 1 ) ) ,
u ( a ) ( x , t ) = u ( x , t ; a ( x ) ) ,
where u ( x , t ; a ( x ) ) is the solution of Equation (10) corresponding to a A . For any a A and ρ A , the G a ^ t e a u x derivative u ( a ) ρ satisfies the homogeneous initial conditions, and for any φ L 2 ( 0 , 1 ) H a 1 ( 0 , 1 ) , the following equality holds:
0 1 ( u ( a ) ρ ) τ φ d x + 0 1 a ( u ( a ) ρ ) x φ x d x + 0 1 a x ( u ( a ) ρ ) x φ d x + 0 1 b ( u ( a ) ρ ) x φ d x + 0 1 c ( u ( a ) ρ ) φ d x = 0 1 u ( a ) x ( ρ φ ) x d x .
Lemma 2.
For any a A and ρ A , such that a + ρ A , we have the following equality
0 1 ( R ( a ) ) τ φ d x + 0 1 a ( R ( a ) ) x φ x d x + 0 1 a x ( R ( a ) ) x φ d x + 0 1 b ( R ( a ) ) x φ d x + 0 1 c R ( a ) φ d x = 0 1 u ( a ) u ( a + ρ ) x ( ρ φ ) x d x ,
where
R ( a ) = u ( a + ρ ) u ( a ) u ( a ) ρ , R ( a ) H 1 ( ( 0 , T ) ; H 0 1 ( 0 , 1 ) ) ,
holds for any φ L 2 ( 0 , 1 ) H a 1 ( 0 , 1 ) .
Proof. 
Notice that u ( a + ρ ) satisfies
0 1 u ( a + ρ ) τ φ d x + 0 1 a φ x + ( ρ φ ) x u ( a + ρ ) x d x + 0 1 a x u ( a + ρ ) x φ d x + 0 1 b u ( a + ρ ) x φ d x + 0 1 c u ( a + ρ ) φ d x = 0 .
Subtracting (7) from (45), and replacing ψ in Equation (7) with φ , we obtain
0 1 u ( a + ρ ) u ( a ) τ φ d x + 0 1 a u ( a + ρ ) u ( a ) x φ x d x + 0 1 a x u ( a + ρ ) u ( a ) x φ d x + 0 1 b u ( a + ρ ) u ( a ) x φ d x + 0 1 c u ( a + ρ ) u ( a ) φ d x = 0 1 u ( a + ρ ) x ( ρ φ ) x d x .
Combining (46) and (43), (44) can be obtained.
In order to get the convergence of the optimal solution, some source conditions need to be added to the direct problem. We introduce the following linear operator F ( a ) :
F ( a ) : L 2 ( ( 0 , T ) ; L 2 ( 0 , 1 ) ) L 2 ( 0 , 1 ) ,
F ( a ) ϕ = 1 p T p T u ( a ) x x ϕ d τ , ϕ L 2 ( ( 0 , T ) ; L 2 ( 0 , 1 ) ) ,
where u ( a ) is the solution of Equation (10). From (43), it is easy to see that for ρ A and φ L 2 ( 0 , 1 ) H a 1 ( 0 , 1 ) , the following equation holds:
F ( a ) φ , ρ = 1 p T p T 0 1 u ( a ) x x ρ φ d x d τ = 1 p T p T 0 1 u ( a ) x ( ρ φ ) x d x d τ = 1 p T p T 0 1 ( u ( a ) ρ ) τ φ + a ( u ( a ) ρ ) x φ x d x d τ + 1 p T p T 0 1 [ a x ( u ( a ) ρ ) x φ + b ( u ( a ) ρ ) x φ + c u ( a ) ρ ] d x d τ ,
where · , · L 2 ( 0 , 1 ) represents the inner product in the L 2 ( 0 , 1 ) . Notice that is a linear operator, now we define * as its conjugate operator such as:
* w , φ L 2 ( 0 , 1 ) = w , φ L 2 ( 0 , 1 ) , w H 1 ( 0 , 1 ) , φ H 1 ( 0 , 1 ) .
If φ H 0 1 ( 0 , 1 ) , then * = . Here, we only use the weak form of ( * ) . □
Theorem 3.
Assume that there exists a function φ
φ H 0 1 ( ( T p , T ) ; H 0 1 ( 0 , 1 ) ) L 2 ( ( T p , T ) ; H 2 ( 0 , 1 ) ) ,
such that in the weak sense the following source condition holds
F ( a * ) φ = a x x * ,
where F ( a * ) is defined by (47). Moreover, as α δ , here are the following estimates:
a α δ a * L 2 ( 0 , 1 ) 2 C δ ,
and
1 p T p T 0 1 u ( a α δ ) u ( a * ) 2 d x d τ C δ 2 ,
where a α δ is the minimum when w replaces w δ in (11). u ( a α δ ) is the solution of variational problem (7) when a = a α δ , and C is independent of δ and α.
Proof. 
Because a α δ is the minimum of (11), then
J p ( a α δ ) J p ( a * ) ,
that is
1 2 p T p T 0 1 u ( a α δ ) w δ 2 d x d τ + α 2 0 1 ( a α δ ) x 2 d x 1 2 δ 2 + α 2 0 1 a x * 2 d x .
From (53), one can deduce
1 2 p T p T 0 1 u ( a α δ ) w δ 2 d x d τ + α 2 0 1 ( a α δ ) x a x * 2 d x 1 2 δ 2 + α 2 0 1 a x * 2 d x α 2 0 1 ( a α δ ) x 2 d x + α 2 0 1 ( a α δ ) x a x * 2 d x 1 2 δ 2 + α 0 1 a x * a * a α δ x d x 1 2 δ 2 + α a x * , ( a * a α δ ) x ,
Combining Equations (48) and (50), for the last term in (54), we have
a x * , ( a * a α δ ) x = a x x * , a * a α δ = F ( a * ) φ , a * a α δ = 1 p T p T 0 1 u ( a * ) x a * a α δ φ x d x d τ = 1 p T p T 0 1 [ u ( a * ) a * a α δ τ φ + a * u ( a * ) a * a α δ x φ x + a x * u ( a * ) a * a α δ x φ ] d x d τ + 1 p T p T 0 1 [ b u ( a * ) a * a α δ x φ + c ( u ( a * ) ( a * a α δ ) φ ] d x d τ .
Let
R α δ : = u ( a α δ ) u ( a * ) u ( a * ) ( a α δ a * ) ,
one can deduce
α a x * , ( a * a α δ ) x = α p T p T 0 1 ( R α δ ) τ φ + a * ( R α δ ) x φ x + a x * ( R α δ ) x φ + b ( R α δ ) x φ d x d τ α p T p T 0 1 u ( a α δ ) u ( a * ) τ φ d x d τ α p T p T 0 1 a * u ( a α δ ) u ( a * ) x φ x d x d τ α p T p T 0 1 a x * u ( a α δ ) u ( a * ) x φ d x d τ α p T p T 0 1 b [ u ( a α δ ) u ( a * ) ] x φ d x d τ α p T p T 0 1 c [ u ( a α δ ) u ( a * ) ] φ d x d τ = I 1 + I 2 + I 3 + I 4 + I 5 + I 6 .
The main idea now is to use the left side of inequality (54) to control I 1 I 6 .
For I 1 , we can deduce from (44)
I 1 = α p T p T 0 1 u ( a * ) u ( a α δ ) x a α δ a * φ x d x d τ .
Integration by part, we obtain
| I 1 | = α p | T p T { u ( a * ) u ( a α δ ) d ( a α δ a * ) φ d x | x = 0 x = 1 0 1 u ( a * ) u ( a α δ ) ( a α δ a * ) φ x x d x } d τ | = α p T p T 0 1 u ( a * ) u ( a α δ ) ( a α δ a * ) φ x x d x d τ .
By the triangle inequality, we have
| I 1 | α p T p T 0 1 u ( a * ) w δ ( a α δ a * ) φ x x d x d τ + α p T p T 0 1 w δ u ( a α δ ) ( a α δ a * ) φ x x d x d τ .
Then from Equations (39), (40) and (58) and using Young’s inequality, we have
| I 1 | 1 2 p T p T 0 1 u ( a * ) w δ 2 d x d τ + C α 2 T p T 0 1 ( a α δ a * ) φ x x 2 d x d τ + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 T p T 0 1 ( a α δ a * ) φ x x 2 d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
For I 2 , from (39), (40), and using integration by parts, trigonometric inequality and Young inequality, we can obtain
| I 2 | = α p T p T 0 1 u ( a α δ ) u ( a * ) τ φ d x d τ α p T p T 0 1 u ( a α δ ) u ( a * ) τ φ d x d τ α p T p T 0 1 u ( a α δ ) w δ τ φ d x d τ + α p T p T 0 1 w δ u ( a * ) τ φ d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
Similarly, for I 3 I 6 , there are
| I 3 | = α p T p T 0 1 a * u ( a α δ ) u ( a * ) x φ x d x d τ = α p T p T a * φ x u ( a α δ ) u ( a * ) | x = 0 x = 1 0 1 u ( a α δ ) u ( a * ) ( a * φ x ) x d x d τ = α p T p T 0 1 u ( a α δ ) u ( a * ) ( a * φ x ) x d x d τ α p T p T 0 1 u ( a α δ ) w δ ( a * φ x ) x d x d τ + α p T p T 0 1 w δ u ( a * ) ( a * φ x ) x d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
| I 4 | = α p T p T 0 1 a x * u ( a α δ ) u ( a * ) x φ d x d τ = α p T p T a x * φ u ( a α δ ) u ( a * ) | x = 0 x = 1 0 1 u ( a α δ ) u ( a * ) ( a x * φ ) x d x d τ = α p T p T 0 1 u ( a α δ ) u ( a * ) ( a x * φ ) x d x d τ α p T p T 0 1 u ( a α δ ) w δ ( a x * φ ) x d x d τ + α p T p T 0 1 w δ u ( a * ) ( a x * φ ) x d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
| I 5 | = α p T p T 0 1 b u ( a α δ ) u ( a * ) x φ d x d τ = α p T p T b φ u ( a α δ ) u ( a * ) | x = 0 x = 1 0 1 u ( a α δ ) u ( a * ) ( b φ ) x d x d τ = α p T p T 0 1 u ( a α δ ) u ( a * ) ( b φ ) x d x d τ α p T p T 0 1 u ( a α δ ) w δ ( b φ ) x d x d τ + α p T p T 0 1 w δ u ( a * ) ( b φ ) x d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
| I 6 | = α p T p T 0 1 c [ u ( a α δ ) u ( a * ) ] φ d x d τ α p T p T 0 1 | u ( a α δ ) u ( a * ) | | c φ | d x d τ α p T p T 0 1 | u ( a α δ ) w δ | | c φ | d x d τ + α p T p T 0 1 | w δ u ( a * ) | | c φ | d x d τ 1 2 δ 2 + 1 24 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 .
Combining (54), (55), and (58) to (64), we obtain
1 2 p T p T 0 1 u ( a α δ ) w δ 2 d x d τ + α 2 0 1 ( a α δ ) x a x * 2 d x 1 2 δ 2 + α a x * , a * a α δ x 1 2 δ 2 + i = 1 6 | I i | 7 2 δ 2 + 1 4 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + C α 2 ,
by (65), one has
1 4 p T p T 0 1 w δ u ( a α δ ) 2 d x d τ + α 2 0 1 ( a α δ ) x a x * 2 d x 3 δ 2 + C α 2 .
Choosing α δ , we obtain
1 p T p T 0 1 u ( a α δ ) w δ 2 d x d τ + α 0 1 ( a α δ ) x a x * 2 d x C δ 2 .
Then by (67) and using P o i n c a r e ´ inequality, we can also obtain (51).
This completes the proof of Theorem 3. □

5. Numerical Experiments

Based on the optimal control theory, we theoretically analyze the optimal control problem P2. In this section, we will use the synthetic data to reconstruct the implied volatility by gradient-type iteration method, and give some numerical examples. To achieve that, we apply the finite difference method to solve the forward problem (4). Since the explicit difference scheme is conditionally stable, the Crank–Nicolson scheme, which is absolutely stable, is employed to obtain the numerical solution. To construct more general testing functions a and u, we add a source function f ( x , τ ) to the right-hand side of Equation (4), i.e., we consider the more general equation:
u τ a ( x ) 2 u x 2 + b ( x ) u x + c ( x ) u = f ( x , τ ) , ( x , τ ) Q = ( 0 , 1 ) × ( 0 , T ] , u ( x , 0 ) = ϕ ( x ) , x ( 0 , 1 ) .
Assume that the domain Q ¯ = [ 0 , 1 ] × [ 0 , T ] is divided into M × N mesh with the spacial step size h = 1 M and the time step size t = T τ , respectively. Grid points ( x j , τ n ) are defined by
x j = j h , j = 0 , 1 , 2 , , M ,
τ n = n t , t = 1 , 2 , , N ,
in which M and N are two integers. The notation u j n is used for the finite difference approximation of u ( j h , n t ) .
The Crank–Nicolson scheme leads to the following difference equation for (68):
u j n + 1 u j n t a j u j + 1 n + 1 2 u j n + 1 + u j 1 n + 1 2 h 2 + u j + 1 n 2 u j n + u j 1 n 2 h 2
+ b j u j + 1 n + 1 u j 1 n + 1 2 h + u j + 1 n u j 1 n 2 h + c j ( u j n + 1 + u j 1 n + 1 2 ) = f j n ,
for 1 j M 1 , and 1 n N .
Note that on the left boundary x = 0 , Equation (68) degenerates into the following ordinary differential equation:
u 0 n + 1 u 0 n t = f 0 n , n = 1 , 2 , , N .
Likewise, on the right boundary x = 1 , we can obtain the following discretized equation:
u M n + 1 u M n t = f M n c j ( u M n + 1 + u M 1 n + 1 2 ) , n = 1 , 2 , , N .
The initial condition can be discretized by
u j 0 = ϕ j , j = 0 , 1 , 2 , , M .
This method is stable in the maximum-norm without any restriction on the mesh parameters t and h, and the truncation error is O ( t 2 + h 2 ) .
For the inverse problem, we consider using the gradient-type iteration method to obtain the numerical solution. The key of this method is to use the G a ^ t e a u x derivative. Before discussing the algorithm, we need to provide the following lemmas:
Lemma 3.
Let a be the solution of the optimal control problem (10). Then there exists a triple of functions ( E , Y ; a ) satisfying the following system:
E τ = a E x x b E x c E + ϕ , ( x , τ ) Q = ( 0 , 1 ) × ( 0 , T ] , E ( x , 0 ) = 0 , x ( 0 , 1 ) ,
Y τ a Y x x b Y x + c Y = 0 , ( x , τ ) Q = ( 0 , 1 ) × ( 0 , T ] , Y ( x , τ ) = u ( x , τ ; a ) w ( x , τ ) , ( x , τ ) ( 0 , 1 ) × [ T p , T ] ,
and
1 p T p T 0 1 Y E x x ( h a ) d x d τ + α 0 1 a ( h a ) d x ,
for any h A .
This lemma represents the necessary conditions for the minimum of optimal control problem, and the proof is similar to [6].
Lemma 4.
The G a ^ t e a u x derivative of J ( a ) in the direction of ρ ( x ) can be expressed as:
J ( a ) ρ = 1 p T p T 0 1 ρ Y E x x d x d τ + α 0 1 a ρ d x ,
where E ( x , τ ; a ) and Y ( x , τ ; a ) are solutions of (69)/(70) at the given coefficient a ( x ) A .
The procedure for the stable reconstruction of the solutions u and a can be stated as follows:
  • Step 1. Choose the initial iteration value a ( x ) = a 0 ( x ) . For simplicity, we choose a 0 ( x ) = c o n s t a n t .
  • Step 2. Solve the initial boundary value problem (4), and obtain the solution u 0 corresponding to a = a 0 .
  • Step 3. Solve the adjoint problem of (69) to get Y 0 , where Y ( x , τ ) = u ( x , τ ; a ) w ( x , τ ) , ( x , τ ) ( 0 , 1 ) × [ T p , T ] .
  • Step 4. Compute the G a ^ t e a u x derivative J ( a ) ψ j = m j , j = 0 , 1 , 2 , , M ,
    m j = 1 p T p T 0 1 ψ j ( x ) Y 0 E x x d x d τ + α 0 1 a ψ j ( x ) d x ,
    where ψ j is a function determined by the following grid:
    ψ 0 ( x ) = x 1 x h x 0 x x 1 , 0 o t h e r ;
    ψ j ( x ) = x x j 1 h x j 1 x x j , x j + 1 x h x j x x j + 1 , 0 o t h e r ;
    ψ M ( x ) = x x M h x M 1 x x M , 0 o t h e r .
Then, the iterative value from step j to step j + 1 can be defined as:
C 0 ( x ) = j = 0 M m j ψ j ( x ) .
  • Step 5. Compute the norm of C 0 ( x ) at step j:
    e = h j = 0 M m j 2 ( x ) 1 2 .
  • Step 6. Select arbitrarily small positive number ε as the error bound, and whether to cease or iterate is determined by the following steps:
Step 6.1. Let k = 1 .
Step 6.2. Compute e r r o r : = J ( a 0 + k C 0 ( x ) ) J ( a 0 ) + 1 2 k e j 2 .
Step 6.3. If e r r 0 , then we execute the Step 6.4; otherwise, let k = θ k , we execute the Step 6.2, where θ is an adjustment parameter.
Step 6.4. Take a 1 ( x ) = a 0 ( x ) + k C 0 ( x ) . If k C 0 ( x ) ε , then iteration should cease; otherwise, let j = j + 1 continue to execute the Step 2. We take a 1 ( x ) as the initial value of the new iteration and follow the above steps until the stopping rule is met.
We present two numerical experiments to test the stability of our algorithm. Since a ( x ) = 1 2 σ 2 ( x ) x 2 ( 1 x ) 2 is the unknown term, it can be directly seen that all but σ 2 ( x ) is known in a ( x ) , we only give the unknown term σ 2 ( x ) in our following examples. For simplicity, in all experiments, the values of some basic parameters are
r = 0 , T = 1 , p = 0.1 , h = t = 0.01 .
We use the symbol ϵ to denote the stopping parameter in the iteration procedure, i.e.,
ϵ = u ( x , τ ; a ) w δ ( x , τ ) L 2 ( Q ) ,
and the symbols ε 1 and ε 2 to denote the absolute and relative L 2 -norm error between the exact solution a ( x ) to be identified and the numerically reconstructed solution a ˜ ( x ) , i.e.,
ε 1 = a ˜ ( x ) a ( x ) L 2 ( Q ) ,
ε 2 = a ˜ ( x ) a ( x ) L 2 ( Q ) a ( x ) L 2 ( Q ) .
Example 1.
In the first numerical experiment, we take
u ( x , τ ) = x ( 1 x ) e τ , ( x , τ ) Q ,
σ 2 ( x ) = 1 x 2 , x ( 0 , 1 ) ,
f ( x , τ ) = x ( 1 x ) e τ [ 1 + ( 1 x 2 ) x ( 1 x ) ] .
The additional condition w ( x , τ ) is given by
w ( x , τ ) = x ( 1 x ) e τ , ( x , τ ) ( 0 , 1 ) × [ 0.9 , 1 ] .
The noisy data w δ ( x , τ ) is generated in the following form:
w δ ( x , τ ) = w ( x , τ ) [ 1 + δ × r a n d o m ( x , t ) ] ,
where the parameter δ is the noisy level.
We present the results of the reconstructed non-linear function σ 2 ( x ) after 100 iterations in Figure 1, where the figure is obtained by k = 100 , 200 , and k = 300 , respectively. It can be seen from this figure that the main shape of the unknown function is recovered well. The initial guess of the unknown term is chosen as σ 2 ( x ) = 1 . In terms of the number of iteration steps alone, it is obvious that the more iterations are better than those with fewer iterations. However, Figure 1 shows that the difference between them is not significant. It can be observed from Figure 1 that the result of k = 300 is only slightly better than that of k = 100 , not necessarily indicating that the more iterations, the more obvious the effect.
We also consider the case of noisy input data to test the stability of our algorithm. The reconstruction results of σ 2 ( x ) from the noisy data w δ ( x , τ ) are shown in Figure 2, where the noise level δ is taken as 0.03, 0.05 and 0.08, respectively. In fact, the numerical solution can achieve a satisfactory accuracy for only about 28 iterations (the L 2 error ε 2 is less than 7.2 × 10 3 ). Satisfactory approximation is obtained even under the case of noisy data. The iteration number j, the stopping parameter ϵ , the accuracy error ε 1 and relative L 2 -norm error ε 2 for various noisy level δ are given in Table 1. We can see from this figure that after 12 iterations (denoted by j) the stopping parameter ϵ is 6.1 × 10 3 which is much less than δ .
Remark 1.
Compared with the exact case, we cannot execute endless iterations for the noisy case. Normally, we have to stop the iteration at a right time, or the data will overflow. This is the cause why the iteration number for the noisy case is much less than the exact case. In fact, if we do not stop the iteration in time after reaching the stopping point, then the residual ϵ may be smaller, but the error of solution ε 1 will become more and more large.
Example 2.
In the second numerical experiment, we consider the piecewise function
u ( x , τ ) = τ sin ( π x ) , ( x , τ ) Q ,
σ 2 ( x ) = sin ( 2 π x ) , 0 < x 1 2 , 4 x 2 , 1 2 < x < 1 ,
f ( x , τ ) = sin ( π x ) [ 1 + π 2 2 σ 2 ( x ) x 2 ( 1 x ) 2 ] .
The additional condition w ( x , τ ) is given by
w ( x , τ ) = τ sin ( π x ) , ( x , τ ) ( 0 , 1 ) × [ 0.9 , 1 ] .
The noisy data w δ ( x , τ ) is generated in the following form:
w δ ( x , τ ) = w ( x , τ ) [ 1 + δ × r a n d o m ( x , t ) ] ,
where the parameter δ is the noisy level.
The initial guess of source function is taken as σ 2 ( x ) = 0 . In this case, we also present the results of the reconstructed unknown function σ 2 ( x ) after 250 iterations in Figure 3, where the figure is obtained by k = 250 , 500 and k = 1000 , respectively. One can see from Figure 3 that although the main shape of the unknown function is reconstructed after 250 iterations, the latter iterative process converges rather slowly. Generally speaking, it is quite difficult to reconstruct the information of unknown functions near the boundary of parabolic equations. From these two examples, one can see that the unknown function in the middle part of the region is well recovered, and the main error appears near the boundary. In summary, our algorithm is satisfactory.
Analogously, we also consider the noisy case, where the noisy levels are δ = 0.001 and δ = 0.01 . The corresponding numerical result is shown in Figure 4. Finally, some computation parameters j , ϵ , ε 1 and ε 2 for the exact case and noisy one are given in Table 2.

6. Conclusions

In this paper, the inverse problem of identifying volatility in option pricing is considered, and the parabolic equation on the unbounded region is transformed into a degenerate parabolic equation on the bounded region through some variable substitutions. At the same time, based on the optimal control framework, it is proved that the optimal approximate solution converges to the true solution of the original problem. The method in this paper can provide a certain reference for the study of inverse problems in option pricing. At the same time, we use the gradient iteration method to obtain the numerical solution of inverse problem, and present some numerical experiments. In this article, we are only using synthetic data, not actual market data. Our goal is simply to support our aforementioned approach and reconstruct the unknown function through stable algorithms. In future work, we will continue to study our problems by applying fast algorithms with the real market data. Based on the rigorous theoretical analysis, the results obtained in this paper are interesting and useful, and people who engage in risk management or volatility trading may apply the method to reconstruct the implied volatility.

Author Contributions

Formal analysis, Y.Y.; supervision, Z.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China (Nos. 11061018, 11261029), Youth Foundation of Lanzhou Jiaotong University (2011028), Long Yuan Young Creative Talents Support Program (No. 252003), and the Joint Funds of NSF of Gansu Province of China (No. 1212RJZA043).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1

We give a brief derivation of P1 from the following equation:
V T + 1 2 σ 2 ( K ) K 2 2 V K 2 r K V K = 0 , ( T , K ) ( t , ) × ( 0 , ) , V ( s , t ; K , T ) | T = 0 = ( s 0 K ) + , K R + .
Let
x = K K + E ; u = V K + E .
Then we have
1 x = E K + E ; x ( 1 x ) = K E ( K + E ) 2 ; K = x E 1 x ,
x 2 ( 1 x ) 2 = K 2 E 2 ( K + E ) 4 ; V = ( K + E ) u ,
and
V T = V u u τ τ T = ( K + E ) u τ ,
V K = V u u x x K + u = ( K + E ) u x E ( K + E ) 2 + u = u x E K + E + u ,
2 V K 2 = ( u x E K + E ) x + u x x K = E K + E 2 u x 2 E ( K + E ) 2 u x + E ( K + E ) 2 u x = E K + E 2 u x 2 .
Substituting these formulas into (A1), we obtain
( K + E ) u τ + 1 2 σ 2 ( K ) K 2 E 2 ( K + E ) 3 2 u x 2 r K ( u x E K + E + u ) = 0 .
Divide both sides by K + E , one has
u τ + 1 2 σ 2 ( K ) K 2 E 2 ( K + E ) 4 2 u x 2 r K E ( K + E ) 2 u x r K K + E u = 0 ,
finally, Equation (A1) transforms into the following degenerate parabolic equation:
u τ 1 2 σ 2 ( x ) x 2 ( 1 x ) 2 2 u x 2 + r x ( 1 x ) u x + r x u = 0 , ( x , τ ) Q = ( 0 , 1 ) × ( 0 , T ] , u ( x , 0 ) = ( s 0 x E 1 x ) + , 0 < x < 1 ,

Appendix A.2

We present an idea of Dupire [12] and derive the partial differential equation dual to the Black–Scholes Equation (1). We consider 1 θ of the following spread option
V 1 ( · ; K , θ ) = 1 θ ( V ( · ; K + θ ) V ( · ; K ) ) .
We know this finite difference satisfies (1) and at the final time T is equal to 0 on ( 0 , K ) , 1 on ( K + θ , ) and linear on ( K , K + θ ) . By the maximum principle, one can conclude that 1 < V 1 ( . ; K , θ ) < 0 on R + × ( t * , T ) . The available theory for parabolic equations ([1,5]) implies that there exists a limit when θ 0 . This limit solves (1) with the terminal data given by the negative Heavyside step function (see [12]). We continue repeating this argument for 1 θ 2 butterflies
V 2 ( · ; K , θ ) = 1 θ 2 ( V ( · ; K θ ) 2 V ( · ; K ) + V ( · ; K + θ ) ) .
Then, we build that there exists a second derivative with respect to K
2 V ( s , t ; K , T ) ) K 2 = G ( s , t ; K , T ) ,
which also solves (1), and
G ( s , t ; K , T ) = δ ( s K ) .
Here, δ ( s K ) is the Dirac delta function concentrated at K. Because G ( s , t ; K , T ) is the fundamental solution to (1) with the final data at t = T , it also satisfies the dual equation (see [12]) with respect to K and T
G T = 1 2 2 K 2 ( K 2 σ 2 ( K ) G ) μ K ( K G ) r G .
By using the definition of G ( s , t ; K , T ) and integrating this equation twice from K to ∞, the second term in the right-hand side can be integrated by parts as follows
K η γ ( γ 2 V γ 2 ) d γ d η = K η 2 V η 2 d η = K V K + K V η d η = K V K V ,
where we have employed the following behaviour at infinity
V , K V K , σ 2 K 2 G , K G K , ( σ 2 K 2 G ) K 0 K .
Finally, this generates the following dual equation for V ( · ; K , T )
V T + 1 2 σ 2 ( K ) K 2 2 V K 2 μ K V K + ( μ r ) V = 0 .
In our paper, we assume μ = r , then we obtain
V T + 1 2 σ 2 ( K ) K 2 2 V K 2 r K V K = 0 , ( T , K ) ( t , ) × ( 0 , ) , V ( s , t ; K , T ) | T = 0 = ( s 0 K ) + .

References

  1. Isakov, V. Inverse Problems for Partial Differential Equations; Springer: New York, NY, USA, 1998. [Google Scholar]
  2. Evans, L.C. Partial Differential Equations; Intersxcience Publishers: Providence, RI, USA, 1964. [Google Scholar]
  3. Oleinik, O. Second-Order Equations with Nonnegative Characteristic Form; Springer Science & Business Media: New York, NY, USA, 2012. [Google Scholar]
  4. Bouchouev, I.; Isakov, V. Uniqueness, stability and numerical methods for the inverse problem that arises in financial markets. Inverse Probl. 1995, 13, R95. [Google Scholar] [CrossRef]
  5. Bouchouev, I.; Isakov, V. The inverse problem of option pricing. Inverse Probl. 1997, 13, 7–11. [Google Scholar] [CrossRef]
  6. Jiang, L.S.; Tao, Y.S. Identifying the volatility of underlying assets from option prices. Inverse Probl. 2001, 17, 137–155. [Google Scholar]
  7. Jiang, L.S.; Chen, Q.H.; Wang, L.J.; Zhang, J.E. A new well-posed algorism to recover implied local volatilit. Quant. Financ. 2003, 3, 451–457. [Google Scholar] [CrossRef]
  8. Jiang, L.S.; Bian, B.J. Identifying the principal coefficient of parabolic equations with non-divergent form. J. Phys. Conf. Ser. 2005, 12, 58–65. [Google Scholar] [CrossRef]
  9. Egger, H.; Engl, H.W. Iikhonov regularization applied to the inverse problem of option pricing: Convergence analysis and rates. Inverse Probl. 2005, 21, 58–65. [Google Scholar] [CrossRef] [Green Version]
  10. Ota, Y.; Kaji, S. Reconstruction of local volatility for the binary option model. Inverse-Ill-Posed Probl. 2016, 24, 727–741. [Google Scholar] [CrossRef]
  11. Xiao, L.; Chen, Z. Taxation problems in the dual model with capital injections. Acta Math. Sci. 2016, 10, 142–149. [Google Scholar]
  12. Dupire, B. Pricing with a Smile. Risk 1994, 10, 18–20. [Google Scholar]
  13. Deng, Z.C.; Yang, L. An Inverse Volatility Problem of Financial Products Linked with Gold Price. Iran. Math. Soc. 2019, 45, 1243–1267. [Google Scholar] [CrossRef]
  14. Isakov, V. Recovery of time dependent volatility coefficient by linearization. Evol. Equ. Control Theory 2014, 3, 119–134. [Google Scholar] [CrossRef]
  15. Deng, Z.C.; Yang, L.; Isakov, V. Recovery of time-dependent volatility in option pricing model. Inverse Probl. 2016, 32, 142–149. [Google Scholar] [CrossRef] [Green Version]
  16. Kangro, R.; Nicolaides, R. Far field boundary conditions for Black-Scholes equations. SIAM J. Numer. Anal. 2000, 38, 1357–1368. [Google Scholar] [CrossRef]
  17. Valkov, R. Fitted finite volume method for a generalized Black-Scholes equation transformed on finite interval. Numer. Algorithms 2014, 65, 195–220. [Google Scholar] [CrossRef] [Green Version]
  18. Windcliff, H.; Forsyth, P.A.; Vetzal, K.R. Analysis of the stability of the linear boundary condition for the Black-Scholes equation. J. Comput. Financ. 2004, 8, 65–92. [Google Scholar] [CrossRef] [Green Version]
  19. Slavi, G.; Vulkov, L.G. Numerical solution of the right boundary condition inverse problem for the Black-Scholes equation. AIP Conf. Proc. 2017, 1910, 30008. [Google Scholar]
  20. Deng, Z.C.; Yang, L. Multi-parameters identification problem for a degenerate parabolic equation. J. Comput. Appl. Math. 2020, 366, 112422. [Google Scholar]
  21. Kamynin, V.L. Inverse problem of determining the absorption coefficient in a degenerate parabolic equation in the class of L2-functions. J. Math. Sci. 2020, 250, 322–336. [Google Scholar] [CrossRef]
  22. Kamynin, V.L. The inverse problem of simultaneous determination of the two time-dependent lower coefficients in a nondivergent parabolic equation in the plane. Math. Notes 2020, 107, 93–104. [Google Scholar] [CrossRef]
  23. Cannarsa, P.; Fragnelli, G. Null controllability of semilinear weakly degenerate parabolic equations in bounded domains. Electron. J. Differ. Equ. 2006, 107, 1–20. [Google Scholar]
  24. Cannarsa, P.; Fragnelli, G.; Rocchetti, D. Null controllability of degenerate parabolic operators with drift. Netw. Heterog. Media 2017, 2, 695–715. [Google Scholar] [CrossRef]
  25. Fragnelli, G. Null controllability of degenerate parabolic equations in non divergence form via Carleman estimates. Discret. Contin. Dyn. Syst.-S 2013, 3, 687–701. [Google Scholar] [CrossRef]
Figure 1. Numerical solution of unknown function σ 2 ( x ) with various iterations for Example 1.
Figure 1. Numerical solution of unknown function σ 2 ( x ) with various iterations for Example 1.
Mathematics 10 02608 g001
Figure 2. Numerical solution of unknown function σ 2 ( x ) for Example 1 with noisy input data.
Figure 2. Numerical solution of unknown function σ 2 ( x ) for Example 1 with noisy input data.
Mathematics 10 02608 g002
Figure 3. Numerical solution of unknown function σ 2 ( x ) with various iterations for Example 1.
Figure 3. Numerical solution of unknown function σ 2 ( x ) with various iterations for Example 1.
Mathematics 10 02608 g003
Figure 4. Numerical solution of unknown function σ 2 ( x ) for Example 2 with noisy input data.
Figure 4. Numerical solution of unknown function σ 2 ( x ) for Example 2 with noisy input data.
Mathematics 10 02608 g004
Table 1. The values of j , ϵ , ε 1 and ε 2 with various noisy level δ for Example 1.
Table 1. The values of j , ϵ , ε 1 and ε 2 with various noisy level δ for Example 1.
δ = 0 δ = 0.03 δ = 0.05 δ = 0.08
j100282012
ϵ 1.22 × 10 7 1.18 × 10 3 3.06 × 10 3 6.1 × 10 3
ε 1 2.34 × 10 3 3.0 × 10 2 7.56 × 10 2 9.4 × 10 2
ε 2 4.76 × 10 4 7.2 × 10 3 1.08 × 10 2 5.91 × 10 2
Table 2. The values of j , ϵ , ε 1 and ε 2 with various noise levels δ for Example 2.
Table 2. The values of j , ϵ , ε 1 and ε 2 with various noise levels δ for Example 2.
j ϵ ε 1 ε 2
δ = 0 2502.38 × 10 6 4.52 × 10 3 8.08 × 10 3
δ = 0.001 1203.18 × 10 5 4 × 10 3 1.01 × 10 2
δ = 0.01 801.81 × 10 4 3.96 × 10 2 6.89 × 10 2
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yimamu, Y.; Deng, Z. Convergence of Inverse Volatility Problem Based on Degenerate Parabolic Equation. Mathematics 2022, 10, 2608. https://doi.org/10.3390/math10152608

AMA Style

Yimamu Y, Deng Z. Convergence of Inverse Volatility Problem Based on Degenerate Parabolic Equation. Mathematics. 2022; 10(15):2608. https://doi.org/10.3390/math10152608

Chicago/Turabian Style

Yimamu, Yilihamujiang, and Zuicha Deng. 2022. "Convergence of Inverse Volatility Problem Based on Degenerate Parabolic Equation" Mathematics 10, no. 15: 2608. https://doi.org/10.3390/math10152608

APA Style

Yimamu, Y., & Deng, Z. (2022). Convergence of Inverse Volatility Problem Based on Degenerate Parabolic Equation. Mathematics, 10(15), 2608. https://doi.org/10.3390/math10152608

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop