In the last 20 years, the parameter inversion problem in

option pricing ﬁeld has been extensively studied by many

scholars, and the results of these studies have all relied on the

famous Black-Sholes model. An important parameter in the

Black-Sholes model is the volatility of the underlying asset

associated with the option, which has a signiﬁcant impact on

market value of the options, and as such many scholars and

practitioners in the ﬁnancial industry have focused intensively

on the volatility of an underlying asset in option pricing.

The derivation of the Black-Scholes partial differential

equations builds on the basic components of derivatives theory,

such as delta hedging and no arbitrage. One of the erroneous

assumptions of the Black-Sholes model is that volatility of

the underlying asset is a constant. Empirical research on

implied volatility shows that implied volatility depends on

strike prices. The value of a call option is obviously a function

of various parameters in the contract, such as strike price K

and expiration time T−t, where Tis the expiration time and

tis the current time. For our inverse problem, we will just use

u(s, t;K, T )for the option value.

Problem P1: Considering the option on the stock without

paying dividend, it is well-known that u(s, t;K, T )for a call

option satisﬁes the following Black-Sholes equation

∂u

∂t +LBS = 0,(s, t)∈R+×(0, T ),

u(s, T ) = (s−K)+= max(0, s −K), s ∈R+(1)

Here, sis the price of underlying stock, Kis the strike price,

Tis the time of expiry, and µand rare, respectively, the risk-

neutral drift and the risk-free interest rate which are assumed

to be constants. The Black-Sholes operator LBS is given by

LBS =1

2σ2(s)s2∂2u

∂s2+sµ ∂u

∂s −ru,

The parameter σ(s)is the volatility coefﬁcient to be identiﬁed.

We assume that

2σ2(s) = 1

2σ2

0+g(s),

where g(s)is small perturbation of constant σ0. Given the

following additional condition:

u(s∗,0, K, T ) = u∗(K, T ), K ∈R+,(2)

where s∗is market price of the stock at time t∗= 0 , and

u∗(K, T )indicates market price of the option with strike K

at a given expiry time T. The inverse problem is to determine

the functions uand σsatisfying (1.1) and (1.2), respectively.

The inverse volatility problem for the Black-Scholes equa-

tion has been discussed intensively in the literature. The

inverse problem was ﬁrst considered by Dupire in [4]. He

applied the symmetric property of the transition probability

density function to replace the option pricing inverse prob-

lem with an equation containing parameters K, T, which

has duality, and proposed Dupire’s formula for calculating

implied volatility. Although this formula is seriously ill-

posed, Dupire’s solution lays an important foundation for later

scholars to study this problem. In [5], the authors reduce

the identiﬁcation of volatility to an inverse parabolic prob-

lem with terminal observation and establish uniqueness and

stability results by using Carleman estimates. This approach

produces a nonlinear Fredholm integral equation in which the

approximated solution is obtained from solving the integral

equation iteratively. In [6], a time-dependent and a space-

dependent volatility have been studied, respectively. A class

of non-Gaussian stochastic processes has been generated in

the study of spatially correlated volatility. The problem is

transformed into a known inverse coefﬁcient problem with

ﬁnal observations and uniqueness and stability theorems are

established by using the dual equations. In [7], L.S. Jiang

used an optimal control framework to determine the implied

Determining the volatility in option pricing from degenerate parabolic

equation

YILIHAMUJIANG YIMAMU

Department of Mathematics, Lanzhou Jiaotong University, Lanzhou, Gansu, 730070,

PEOPLE’S REPUBLIC OF CHINA

Abstract: This contribution deals with the inverse volatility problem for a degenerate parabolic equation from

numerical perspective. Being different from other inverse volatility problem in classical parabolic equations,

the model in this paper is degenerate parabolic equation. Due to solve the deficiencies caused by artificial

truncation and control the volatility risk with precision, the linearization method and variable substitutions are

applied to transformed the inverse principal term coefficient problem for classical parabolic equation into the

inverse source problem for degenerate parabolic equation in bounded region. An iteration algorithm of

Landweber type is designed to obtain the numerical solution of the inverse problem. Some numerical

experiments are performed to validate that the proposed algorithm is robust and the unknown coefficient is

recovered quite well.

Keywords: Inverse volatility problem, Linearization method, Landweber iteration, Numerical experiments

Received: July 24, 2021. Revised: June 24, 2022. Accepted: July 28, 2022. Published: September 13, 2022.

1. Introduction

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

629

Volume 21, 2022

volatility, and carried out a rigorous mathematical analysis

of the inverse problem, proving that the approximate optimal

solution converges to an appropriate solution to the original

problem. Similar results are derived in [15]. In [8]-[9], the

inverse problem of identifying the principal coefﬁcient is

investigated when the solution is known, and a well-posed ap-

proximation algorithm to identify the coefﬁcient is proposed.

The existence, uniqueness and stability of such solutions are

proved. The Tikhonov regularization method has always been

an important tool for solving ill-posed problems. In [9], a new

two-dimensional numerical differentiation method is proposed

through Tikhonov regularization. Convergence analysis and

numerical examples are given. The authors of [10] studied

the stable identiﬁcation problem of the local volatility surface

σ(S, t)in the Black-Scholes/Dupire equation from the market

price of European options. The stability and convergence of the

approximation obtained by Tikhonov regularization. In case of

a known term-structure of volatility, based on the assumption

that the volatility is constant in time σ(S, t) = σ(S), the con-

vergence rate under simple smoothness and decay conditions

on the true volatility is proved. In recent years, linearization

techniques have been applied to the inverse problem of option

pricing. In [11,12,14], linearization techniques are applied to

transform the problem into an inverse source problem, from

which unknown volatility can be recovered. A stable numerical

solution to the inverse problem is obtained by using the

integral equation method and the Landweber iteration method.

Both the theoretical analysis and the numerical examples

demonstrate the effectiveness of the proposed method.

It is worth mentioning that the aforementioned scholars

and their research have made outstanding contributions to

the inverse volatility problem in option pricing. However,

there are some deﬁciencies in these studies that need to be

improved. One of the signiﬁcant deﬁciencies is to consider

that the original problem is in the unbounded region, so

many scholars conduct numerical simulations by artiﬁcial

truncation. There is a potential trouble in this approach, that

is, if we truncate the interval 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 will fail to precisely control the volatility

risk. In order to solve this deﬁciency, our main objective

in this paper is to introduce a new variable substitutionand

the linearization method that transforms original problem

into the degenerate problem. Since the problem is ill-

posedness, we adapt tikhonov regularization and design the

landweber-type iteration to solve the invese valatility problem.

We shall assume that the option price premium u(·,·;K, T )

satisﬁes the equation dual to the Black-Sholes equation (1)

with respect to the strike price Kand expiry time T:

∂u

∂T −1

2K2σ2(K)∂2u

∂K2+µK ∂u

∂K + (r−µ)u= 0.(3)

The equation (3) was found by Dupire in [4].

In the lognormal variables

y= ln K

s∗, τ =T,

a(y) = σ(s∗ey), U(y, τ ) = u(s∗ey, τ ),(4)

The inverse problem P1 transforms into

Uτ−LU = 0, τ > 0,

U(y, 0) = s∗(1 −ey)+, y ∈R, (5)

where

LU =1

2a2(y)Uyy −1

2a2(y, τ) + µUy−(r−µ)U,

with the additional market data

U(y, T ) = U(y), y ∈R, (6)

where

U(y) = u∗(s∗ey, T ).(7)

The goal is to recover the unknown space-dependent volatil-

ity coefﬁcient a(y)from market data U(y).

If 1

2a2(y) = 1

2σ2

0+g(y),

we have

U=V0+V+v,

where V0is the solution of (5) when a=σ0and vis

quadratically small with respect functions g. The principal

linear term Vsatisﬁes

Vτ−AV =α0(y, τ)g(y),

V(y, 0) = 0, y ∈R, (8)

where

AV =1

2σ2

0Vyy −1

2σ2

0+µVy−(r−µ)V,

with the additional ﬁnal data

V(y, T ) = V(y) = U(y)−V0(y, T ), y ∈R, (9)

in which

α0(y, τ) = s∗1

σ0√2πτ e

−y2

2τ σ2

+c′y+dτ ,

c′=1

2+µ

σ2

, d =−1

2σ2

0σ2

2+µ2

+µ−r,

is known.

The recovery of gin (5) and (6) is a linear inverse source

problem. However, this is a matter of unbounded areas,

we proposed some new variable substitutions here, so that

the above problem is transformed into linear inverse source

degenerate parabolic problem on bounded areas.

Taking

x= arctan y, W =V. (10)

2. Preliminary Knowledge

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

630

Volume 21, 2022

Problem P: Consider the following parabolic equation:

Wτ−BW =α(x, τ)g(x),(x, τ)∈Q=−π

2,π

2×(0, T ],

W(x, 0) = 0,−π

2< x < π

(11)

where

BW =1

2σ2

0cos4xWxx

−σ2

0sinxcos3x+1

2σ2

0+µcos2xWx

−(r−µ)W,

with the additional ﬁnal data

V(y, T ) = W(tan x, T ) = w(x), x ∈−π

2,π

2,(12)

in which

α(x, τ) = s∗1

σ0√2πτ e

−(tan x)2

2τ σ2

+c′∗tan x+d∗τ,

c′=1

2+µ

σ2

, d =−1

2σ2

0σ2

2+µ2

+µ−r,

For the sake of analysis, we take

a(x) = σ2

0cos4x,

b(x) = σ2

0sin xcos3x+1

2σ2

0+µcos2x,

c= (r−µ),

W(x, 0) = ϕ(x),

f(x, τ) = α(x, τ)g(x),

where a(x), b(x), c, ϕ(x)and α(x, τ)are given smooth func-

tion on (−π

2,π

2)which satisﬁes

a(−π

2) = a(π

2) = 0, a(x)>0, x ∈(−π

2,π

2),(13)

b(−π

2) = b(π

2) = 0,

and g(x)is an unknown source term in (11). We shall

determine the functions Wand gsatisfying (11) and (12).

There are many tools such as optimal control frameworks

that can analyze the uniqueness and stability of solutions to

such inverse problems. In this paper, we solve the inverse

source problem from a numerical perspective. Since the in-

verse problem P is ill-posedness, the Tikhonov regularization

method based on the L2gradient norm should be adopted. We

consider the following linear system:

T x =y, x ∈X, y ∈Y,

where Xand Yare Hilbert spaces, and T:X→Yis

a bounded linear operator. Then Tikhonov functional can be

written as follows:

J(x) := kT x −yk2+αkxk2, x ∈X.

Obviously, the minimal element of the functional described

above is equivalent to the solution of the following equation

αx +T∗T x =T∗y,

which is

x= (αI +T∗T)−1T∗y.

However, this method is not suitable for solving the problem

of this article. There are two main difﬁculties. First, we do

not know the speciﬁc form of (αI +T∗T)−1. In fact, we can

write the speciﬁc form of (αI +T∗T)−1only for individual

operators, for example, Tis a matrix. Second, a second

derivative term of the unknown function appears in the form

of Euler’s equation, which makes the numerical simulation

process very complicated.

Therefore, we use the iterative method to solve the inverse

problem P. In this article, we particularly use the Landweber

iteration method to get the numerical results.

Deﬁne the linear operator as shown below:

K:L2(−π

2,π

2)→H1

a(−π

2,π

2),(14)

Kg =W(·, T ) = w(x),(15)

where Wis the solution of equation (11) under the following

initial value condition:

ϕ(x)≡0, x ∈(−π

2,π

2).(16)

Noticed that (15) can also be written as the following form:

g= (I−αK∗K)g+αK∗w. (17)

for α > 0,K∗is the adjoint operator of K, so the following

iterative format can be used to solve equation (17):

g= 0,

gm= (I−αK∗K)gm−1+αK∗w, m = 1,2,3,···.

(18)

It is easy to verify that equation (18) is the fastest descent

method to solve the following equation, where αis the step

size.

φ(g) = 1

2kKg −wk.(19)

Besides, there is the following lemma for the conjugate

adjoint K∗.

Lemma 1. For any given h(x)∈H1

a(−π

2,π

2),let ω(x, 0) =

K∗h, which is

K∗:H1

a(−π

2,π

2)→L2(−π

2,π

2),

K∗h=ω(x, 0),

then ωsatisﬁes the following parabolic equation:







−ωτ−(aω)xx −(bω)x+cω

=α(x, τ)h(x),(x, τ)∈Q,

ω(x, T ) = 0.

The proof of lemma 1 is similar to the references(see [15]).

3. Landweber Iterations

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

631

Volume 21, 2022

We have already known the Kis a linear mapping, so

problem (11)-(12) can be transformed into the ﬁrst type of

operator equation:

Kg =w−Hϕ. (20)

Noticing the (20) can be rewritten as the following form:

g= (I−αK∗K)g+αK∗(w−Hϕ).(21)

so the following iterative format can be used to solve equation

(21):







g= 0,

gm= (I−αK∗K)gm−1+αK∗(w−Hϕ),

m= 1,2,3,···.

(22)

equation (22) is the fastest descent method to solve the

following equation, where αis the step size.

φ(g) = 1

2kKg −(w−Hϕ)k.(23)

from the deﬁnition of Hand (22), we have:

gm=gm−1−αK∗(Kgm−1−(w−Hϕ))

=gm−1−αK∗(Wm−1−(·, T )−w),(24)

where Wm−1is the solution of (11)-(12), when g=gm−1.

From equation (24) and the deﬁned K∗, the following adjoint

equation is introduced:







−ωτ−(aω)xx −(bω)x+cω

=W(x, T )−w(x),(x, τ)∈Q,

ω(x, T ) = 0.

(25)

Assuming that the true solution w(x)is available, that is,

there is a g(x)∈L2(−π

2,π

2)such that

W(x, T ;g) = w(x),

and the noise of the observation data has an upper bound δ,

that is,

kwδ−wkL2(−π

2,π

2)≤δ.

In summary, the calculation steps of the iteration format can

be stated as follows:

Step one: Choose an initial iterative function g=g(x). The

initial function can be selected arbitrarily, for the convenience

of calculation. We generally choose g(x) = 0, x ∈(−π

2,π

2);

Step two: W0(x, τ)is obtained by solving the initial

boundary value problem (11), where g=g(x);

Step three: Solve the equation (25) to get ω0(x, τ), where

W(x, T ) = W0(x, T );

Step four: Let g1(x) = g(x)−αω0(x, T ), where α≥0,

and let W1(x, τ)be the solution of (11) when g=g1(x);

Step ﬁve: Choose an arbitrarily small normal number εas

the error limit, calculate kW1(x, T )−w(x)kand compare the

size with ε, if:

kW1(x, T )−w(x)k< ε,

then terminate the iteration, and take g=g1(x)at this time.

Otherwise, continue to execute step three, and let g1(x)be

the new initial value of the iteration and continue to execute

the inductive criterion until the iteration meets the termination

condition.

Normally, if the input data is accurate, the more iterations

of the Landweber iterative method, the higher the accuracy

of the output data. However, in the case of noise, there will

be errors in the initial iteration process. This calculation error

will initially decrease as the number of iterations increases,

but when a certain threshold is reached, it will increase rapidly

as the number of iterations increases. Therefore, the iteration

must be terminated at the appropriate time. In other words, in

order to balance accuracy and stability, a compromise solution

must be found, that is, a suitable parameter must be selected

so that the iteration format is both accurate and stable method.

If x∈(K∗K)r(X), r ∈N, then the following error

estimation formula can be obtained:



xN(δ),δ −x



≤CM 1

2r+1 δ2r

2r+1 ,

where Mis the boundary of (K∗K)rx. Therefore, in this case

it is different from Tikhonov’s regularization method.

We would like to give some numerical examples to test the

validity of the proposed methods in section 2. The simulated

data are generated by using the standard ﬁnite difference

method to solve the direct problems (11) and (12) under some

appropriate boundary conditions. We use T= 1, for numerical

convenience, xinterval (−π

2,π

2)is divided into 100 equal

interval and on xaxis we show number of an interval. It is

same for the case of τ.

Example 1. Take

g(x) = (cos x, x ∈[−π

2,π

2],

0, others.

and

α(x, τ) = 1 + 1

2τcos4x−τ(sin2xcos2x+ sin xcos x),

r=µ= 0.5,

the numerical results are shown in Fig. 1, where the iteration

number k= 1000 . It can be seen from this ﬁgure that the

main shape of unknown functions is recovered well.

0 50 100 150 200 250

x (-0.5pi~0.5pi)

-0.2

0.2

0.4

0.6

0.8

1.2

g(x)

k=1000

exact solution

4. Numerical Experiments

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

632

Volume 21, 2022

Fig. 1. Numerical solution of source function g(x)for Example 1, where

k= 1000

We also consider the case of noisy input data to test the

stability of our algorithm. The noisy data are generated in the

following form:

w(x) = Wδ(x, T ) = W(x, T )[1+δ×random(x)], x ∈[−π

2,π

2],

with δ= 0.001,0.01,0.08. The reconstruction result is

displayed in Fig. 2 in which satisfactory approximation is

obtained under the case of noisy data as well. As δ=

0.001,0.01,0.08, the corresponding iteration numbers are k=

1500,1200 and k= 1000, respectively. Since the observation

data contains error, to obtain stable numerical results, we shall

cease the iteration at some suitable time.

0 50 100 150 200 250

x (-0.5pi~0.5pi)

-0.2

0.2

0.4

0.6

0.8

1.2

g(x)

=0.01

=0.08

=0.1

exact solution

Fig. 2. Numerical solution of source function g(x)for Example 1 with

noisy input data.

Example 2. In the second numerical experiment, we take

T= 1, α(x, τ) = x2+6τcos4x+4xτ(sin xcos3x+cos2x),

r=µ= 0.5,

g(x) = (x2, x ∈[−π

2,π

2],

0, others.

the numerical results are shown in Fig. 3, where the iteration

number k= 800 .

0 50 100 150 200 250

x (-0.5pi~0.5pi)

0.5

1.5

2.5

g(x)

k=800

exact solution

Fig. 3. Numerical solution of source function g(x)for Example 2,where

k= 800.

Generally speaking, it is not easy to reconstruct the infor-

mation of unknown function near the boundary of parabolic

equations. From this ﬁgure, we can see the main error which

appears near the boundary is very small. Analogously, we

also consider the noisy case, where the noisy levels are

same to those in Example 1, i.e., δ= 0.001,0.01,0.08.The

corresponding numerical result is displayed in Fig. 4. One

can see that for the noisy case, our algorithm is still stable

and the unknown function is reconstructed very well. As

δ= 0.001,0.01,0.08,the corresponding iteration numbers are

k= 980,820 and 500, respectively.

0 50 100 150 200 250

x (-0.5pi~0.5pi)

-0.5

0.5

1.5

2.5

g(x)

=0.001

=0.01

=0.08

exact solution

Fig. 4. Numerical solution of source function g(x)for Example 2 with

noisy input data.

We investigate the inverse volatility problem from numerical

perspective. We apply the linearization method and variable

substitutions to transform the inverse principal term coefﬁcient

problem for classical parabolic equation into the inverse source

problem for a degenerate parabolic equation. We design an

iteration of Landweber-type to obtain the numerical solution

of the inverse problem and present several experiments to show

that the proposed algorithm is robust.

This work was supported by National Natural Science

Foundation of China (Nos.11061018, 11261029).

[1] V. Isakov. Inverse problems for partial differential equations, Springer,

New York, 1998.

[2] Evans, Lawrence C. Partial differential equations. Intersxcience Publish-

ers, 1964.

[3] Oleˆınik, O. A, Radkevich E V, Fife P C. Second order equations with

nonnegative characteristic form[M]. American Mathematical So, 1973.

[4] Dupire B. Pricing with a Smile[J]. risk, 1994.

[5] I. Bouchouev, V. Isakov. Uniqueness, stability and numerical methods for

the inverse problem that arises in ﬁnancial markets[J]. Inverse Problems,

1999, 15(3):R95.

[6] I. Bouchouev, V. Isakov. The inverse problem of option pricing[J].

Inverse Problem, 1997, (13): 7-11.

[7] L.S. Jiang, Y.S. Tao. Identifying the volatility of underlying assets from

option prices[J]. Inverse Problems, 2001, (17): 137-155.

[8] L.S. Jiang, Q.H. Chen, L.J. Wang, J.E. Zhang. A new well-posed

algorism to recover implied local volatility[J]. Quantitative Finance,

2003, (3): 451-457.

5. Conclusion

Acknowledgment

References

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

633

Volume 21, 2022

[9] L. S. Jiang and B. J. Bian. Identifying the principal coefﬁcient of

parabolic equations with non-divergent form. Journal of Physics: Con-

ference Series 12 (2005) 58-65.

[10] Egger H , Engl H W . Tikhonov regularization applied to the inverse

problem of option pricing: convergence analysis and rates[J]. Inverse

Problems, 2005, 21(3).

[11] Z.C. Deng, L. Yang. An Inverse Volatility Problem of Financial Products

Linked with Gold Price[J]. Iranian Mathematical Society, 2019, 45(4):

1243-1267.

[12] Isakov V. Recovery of time dependent volatility coefﬁcient by lineariza-

tion[J]. Evol. Equ. Control Theory, 2014, (3): 119-134.

[13] Yasushi Ota, Shunsuke Kaji. Reconstruction of local volatility for the

binary option model[J]. Inverse Ill-Posed Problem, 2016, 24 (6): 727-

741.

[14] Z.C. Deng, Y.C. Hon and V Isakov. Recovery of time-dependent

volatility in option pricing model[J]. Inverse Problems, 2016, (32):30.

[15] L. Xiao, Z. Chen. Taxation problems in the dual model with capital

injections[J]. Acta Mathematica Scientia, 2016.

Creative Commons Attribution License 4.0

(Attribution 4.0 International, CC BY 4.0)

This article is published under the terms of the Creative

Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en_US

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.73

Yilihamujiang Yimamu

E-ISSN: 2224-2880

634

Volume 21, 2022