# Construction of Interval Shannon Wavelet and Its Application in Solving Nonlinear Black-Scholes Equation

Authors: Liwei Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Interval wavelet numerical method for nonlinear PDEs can improve the calculation precision compared with the common wavelet. A new interval Shannon wavelet is constructed with the general variational principle. Compared with the existing interval wavelet, both  the gradient and the smoothness near the boundary of the approximated function are taken into account. Using the new interval Shannon wavelet, a multiscale interpolation wavelet operator was constructed in this paper, which can transform the nonlinear partial differential equations into matrix differential equations; this can be solved by the coupling technique of the wavelet precise integration method (WPIM) and the variational iteration method (VIM). At last, the famous Black-Scholes model is taken as an example to test this new method. The numerical results show that this method can decrease the boundary effect greatly and improve the numerical precision in the whole definition domain compared with Yan’s method.


---

## Body

## 1. Introduction

The nonlinear Black-Scholes equations have been increasingly attracting interest over the last two decades, since they provide more accurate values by taking into account more realistic assumptions, such as the transaction costs, risks from an unprotected portfolio, large investor’s preferences, or illiquid markets, which may have an impact on the stock price, the volatility, the drift, and the option price itself [1]. In contrast with the traditional Black-Scholes equation, it is almost impossible to find the analytical solution of the nonlinear Black-Scholes equations. Most of the numerical schemes focus on the finite difference methods [2, 3]. As the matter of fact, the wavelet precise integration method (WPIM) is a simple and effective method for linear partial differential equations proposed by Mei et al. [4]. For linear steady structural dynamic systems, their numerical results at the integration points are almost equal to those of the exact solution in machine accuracy. Yan proposed a wavelet precise integration method for the nonlinear Black-Scholes model [5] recently by combining the variational iteration method (VIM) and homotopy perturbation method (HPM) [6].Since the definition domain of wavelet transformation is an infinite interval, the boundary effect would occur when being applied for resolving the engineering problems with bounded interval, for example, ordinary differential equations (ODEs). Consequently, it will decrease the precision and computational efficiency of the solution. Nevertheless, the boundary effect can be eliminated effectively by constructing an interval wavelet using numerical methods. There are several ways available to construct an interval wavelet. Mei et al. [7] proposed a general construction method of interval wavelet based on the restricted variational principle. According to this method, there exist some interval wavelets that have been constructed to solve PDEs in engineering [8–11].The main purpose of this paper is to deduce a novel interval wavelet based on the restricted variational principle and construct an interval wavelet precise integration method (WPIM) for the nonlinear Black-Scholes model. According to WPIM, the nonlinear differential equation should be transformed to a system of ordinary differential equations with the multiscales wavelet interpolation operator [12], and then the nonlinear PDEs become a system of nonlinear differential equations. The matrix differential equation (MDE) can be solved almost exactly by the coupling technique of the precise integration method and the variational iteration method (VIM) [13–17]. In contrast to the traditional finite difference approximation, the numerical results obtained by the PIM are close to the computer precision and also are free from the stiff problem in solving a set of simultaneous linear time-invariant ODEs.


## 2. Construction of Interval Shannon Wavelet Based on Restricted Variational Principle

Mei et al. [7] proposed a kind of construction method of the interval interpolation wavelet based on the restricted variational principle. It does not matter with the representation of the wavelet function. In this section, a new interval wavelet function with this method will be constructed. It takes into account the gradient and the smoothness of the approximated function. This is helpful to improve the calculation precision of the wavelet numerical method.The representation of Shannon wavelet [18–21] is based upon approximating the Dirac delta function as a band-limited function and is given by
(1)
ϕ
(
x
)
=
sin
(
π
x
)
π
x
.
The Shannon wavelet possesses many excellent numerical properties such as interpolating, relative sparse, and orthogonal properties. A perceived disadvantage of (1) is that it tends to zero quite slowly as |
x
|
→
∞. A direct consequence of this is that there are a large number of grid points will contribute to the derivatives calculation of approximated function. It is for this reason that Hoffman et al. [22] have suggested using the Shannon-Gabor wavelet as follows:
(2)
w
(
x
)
=
sin
(
π
x
)
π
x
exp
⁡
(
-
x
2
2
σ
2
)
,
σ
>
0
.
It has been proved that (2) can improve the localized and asymptotic behavior of the Shannon scaling function. σ is the width parameter (or called window size). However, the presence of the Gaussian window destroys the orthogonal properties possessed by the Shannon wavelet, effectively worsening the approximation to a Dirac delta function. So, the Shannon wavelet representation of the Dirac delta function is adopted in this paper, and it is shown that this representation ensures that the approach is identical to the weighted residual approach. The corresponding discrete formula of the Shannon scaling function is
(3)
w
(
x
-
x
i
)
=
sin
⁡
(
π
/
Δ
)
(
x
-
x
i
)
(
π
/
Δ
)
(
x
-
x
i
)
.
In the following, the construction of the interval Shannon wavelet is based on the general variational principle [7]. Consider one-dimensional problem u
(
x
), whose domain of definition Ω is [
m
,
n
]. Assuming that the number of discrete points is 2
j
+
1, j
∈
Z, the ith discrete point of variable x could be written as
(4)
x
i
=
m
+
n
-
m
2
j
·
i
,
i
∈
Z
.
If the Shannon wavelet with interpolation property is employed as the trial function N
j, that is,
(5)
N
j
=
w
j
(
x
-
x
i
)
=
sin
⁡
(
2
j
π
/
(
n
-
m
)
)
(
x
-
x
i
)
(
2
j
π
/
(
n
-
m
)
)
(
x
-
x
i
)
,
any element in the matrices K and G
T could be expressed, respectively, as
(6)
K
k
,
n
=
A
(
w
j
(
k
-
n
)
)
,
(7)
G
k
,
n
T
=
∫
Ω
N
~
k
T
L
1
(
N
n
)
d
Ω
,
where N
~
k is the Lagrange basic function.As the definition domain of wavelet transformation is a double infinite interval, the wavelet coefficients are large near the endpoints of the bounded signal, which increases the computational error. To obtain a high computational precision, the additional condition of scalar fonctionelle in boundaryΓ could be changed to
(8)
u
=
u
.
-
Thus
(9)
L
1
(
N
n
)
=
∑
i
N
i
u
-
i
,
where u
-
i is a known value. According to the general variational principle, we obtain
(10)
K
0
C
+
R
=
[
K
G
G
T
0
]
{
a
b
}
-
{
P
0
}
=
0
,
where
(11)
K
=
∫
Ω
N
T
A
(
N
)
d
Ω
,
G
T
=
∫
Ω
N
~
T
L
1
(
N
)
d
Ω
,
Q
=
∫
Ω
N
~
T
C
1
d
Ω
.
Correspondingly, (7) changed to the form
(12)
G
k
,
n
T
=
N
~
k
T
N
n
.
The interpolation properties of N
~ and N imply that a and b are the vectors consisting of the approximate values of u in discrete points x
i and extended discrete points outside the domain of [
m
,
n
], respectively. Typically, periodic, symmetric, and zero extension methods could be employed. Assuming that there are R
=
2
j discrete points in [
m
,
n
], that is, x
0
,
x
1
,
…
,
x
2
j, and L extended points in both outsides of the domain [
m
,
n
], respectively, that is, x
-
L
-
1
,
x
-
L
,
…
,
x
-
1 and x
R
+
1
,
x
R
+
2
,
…
,
x
R
+
L, the function u can be extended outside of the definition domain using the symmetry method as below:
(13)
u
(
x
-
1
)
=
2
u
(
x
0
)
-
u
(
x
1
)
,
u
(
x
-
2
)
=
2
u
(
x
0
)
-
u
(
x
2
)
,
aaaaaaai
⋮
u
(
x
-
1
-
i
)
=
2
u
(
x
0
)
-
u
(
x
1
+
i
)
,
aaaaaaaai
⋮
u
(
x
-
L
-
1
)
=
2
u
(
x
0
)
-
u
(
x
1
+
L
)
,
u
(
x
R
+
1
)
=
2
u
(
x
R
)
-
u
(
x
R
-
1
)
,
u
(
x
R
+
2
)
=
2
u
(
x
R
)
-
u
(
x
R
-
2
)
,
aaaaaaaai
⋮
u
(
x
R
+
i
)
=
2
u
(
x
R
)
-
u
(
x
R
-
i
)
,
aaaaaaaa
⋮
u
(
x
R
+
L
)
=
2
u
(
x
R
)
-
u
(
x
R
-
L
)
.
Substituting (13) into (10) results in
(14)
[
G
1
A
G
2
]
{
u
-
L
-
1
u
-
L
⋮
u
-
1
u
0
u
1
⋮
u
2
j
u
R
+
1
u
R
+
2
⋮
u
R
+
L
}
=
{
P
′
0
}
,
where
(15)
A
=


[


A
(
w
j
(
0
-
0
)
)
A
(
w
j
(
0
-
1
)
)
⋯
A
(
w
j
(
0
-
2
j
)
)
A
(
w
j
(
1
-
0
)
)
A
(
w
j
(
1
-
1
)
)
⋯
A
(
w
j
(
1
-
2
j
)
)
⋮
⋮
⋮
⋮
A
(
w
j
(
2
j
-
0
)
)
A
(
w
j
(
2
j
-
1
)
)
⋯
A
(
w
j
(
2
j
-
2
j
)
)


]

,
G
1
=


[


a
-
L
-
1
,
-
L
-
1
w
j
(
0
)
a
-
L
-
1
,
-
L
w
j
(
-
1
)
⋯
a
-
L
-
1
,
-
1
w
j
(
-
L
)
a
-
L
,
-
L
-
1
w
j
(
1
)
a
-
L
,
-
L
w
j
(
0
)
⋯
a
-
L
,
-
1
w
j
(
1
-
L
)
⋮
⋮
⋮
⋮
a
-
1
,
-
L
-
1
w
j
(
L
)
a
-
1
,
-
L
w
j
(
L
-
1
)
⋯
a
-
1
,
-
1
w
j
(
0
)


]

,
G
2
=
[
b
R
,
R
w
j
(
0
)
b
R
,
R
+
1
w
j
(
-
1
)
⋯
b
R
,
R
+
L
w
j
(
-
L
)
b
R
+
1
,
R
w
j
(
1
)
b
R
+
1
,
R
+
1
w
j
(
0
)
⋯
b
R
+
1
,
R
w
j
(
1
-
L
)
⋮
⋮
⋮
⋮
b
R
+
L
,
R
w
j
(
L
)
b
R
+
L
,
R
+
1
w
j
(
L
-
1
)
⋯
b
R
+
L
,
R
+
L
w
j
(
0
)
]
.
In the matrices G
1 and G
2, a
n
k and b
n
k could be calculated, respectively, as
(16)
a
n
k
=
l
j
k
1
(
x
j
n
)
,
b
n
k
=
l
j
k
2
(
x
j
n
)
,
l
j
,
k
1
=
∏
i
=
-
L
-
1
i
≠
k
-
1
x
-
x
j
,
i
x
j
,
k
-
x
j
,
i
,
l
j
,
k
2
=
∏
i
=
2
j
+
1
i
≠
k
2
j
+
1
+
L
x
-
x
j
,
i
x
j
,
k
-
x
j
,
i
.
According to G
T
b
=
0 and the relational expression (13), (14) could be condensed in the following way:
(17)
[
G
1
+
A
11
A
12
A
13
A
21
A
22
A
23
A
31
A
32
G
2
+
A
33
]
{
u
0
u
1
⋮
u
2
j
}
=
{
P
0
}
.
Therefore, the function u(x) can be represented as
(18)
u
(
x
)
=
∑
k
=
-
N
2
j
+
N
w
(
2
j
x
-
k
)
u
k
j
=
∑
k
=
-
N
-
1
w
(
2
j
x
-
k
)
(
2
f
0
-
u
-
k
j
)
+
∑
k
=
0
2
j
w
(
2
j
x
-
k
)
u
k
j
+
∑
k
=
2
j
+
1
2
j
+
N
w
(
2
j
x
-
k
)
(
2
u
2
j
-
u
2
j
+
1
-
k
)
=
-
∑
k
=
-
N
-
1
w
(
2
j
x
-
k
)
u
-
k
j
+
2
∑
k
=
-
N
-
1
w
(
2
j
x
-
k
)
u
0
+
∑
k
=
0
2
j
w
(
2
j
x
-
k
)
u
k
j
+
2
∑
k
=
2
j
+
1
2
j
+
N
w
(
2
j
x
-
k
)
u
2
j
-
∑
k
=
2
j
+
1
2
j
+
N
w
(
2
j
x
-
k
)
f
2
j
+
1
-
k
=
-
∑
k
=
1
N
w
(
2
j
x
+
k
)
u
-
k
j
+
2
∑
k
=
1
N
w
(
2
j
x
+
k
)
u
0
+
∑
k
=
0
2
j
w
(
2
j
x
-
k
)
u
k
j
+
2
∑
k
=
2
j
-
N
-
2
j
-
1
w
(
2
j
x
+
k
)
u
2
j
-
∑
k
=
2
j
+
1
2
j
+
N
w
(
2
j
x
-
k
)
u
2
j
+
1
-
k
,
where
(19)
∑
k
=
0
2
j
w
(
2
j
x
-
k
)
f
k
j
=
w
(
2
j
x
)
f
0
j
+
∑
k
=
1
N
w
(
2
j
x
-
k
)
f
k
j
+
∑
k
=
N
+
1
2
j
-
N
-
1
w
(
2
j
x
-
k
)
f
k
j
+
∑
k
=
2
j
-
N
2
j
-
1
w
(
2
j
x
-
k
)
f
k
j
+
w
(
2
j
x
-
2
j
)
f
2
j
j
.
Subsequently, as interval interpolation basic functions, the interval Shannon wavelet can be obtained from (18) as follows:
(20)
w
j
k
=
{
w
(
2
j
x
-
0
)
+
2
∑
k
=
1
N
w
(
2
j
x
+
k
)
,
k
=
0
,
w
(
2
j
x
-
k
)
-
w
(
2
j
x
+
k
)
,
k
=
1,2
,
…
,
N
,
w
(
2
j
x
-
k
)
,
k
=
N
+
1
,
…
,
i
i
i
i
2
j
-
N
-
1
,
w
(
2
j
x
-
k
)
-
w
(
2
j
x
-
2
j
+
1
+
k
)
,

k
=

2
j
-
N
,
…
,
2
j
-
1
,
w
(
2
j
x
-
2
j
)
+
2
∑
k
=
2
j
+
1
2
j
+
N
w
(
2
j
x
-
k
)
,
k
=
2
j
,
where N is the support domain of the wavelet function; that is, sup ϕ
=
[
-
N
,
N
].The interval wavelet constructed in [8] is different from (20), which is represented as follows:
(21)
w
j
k
=

{
w
(
2
j
x
-
k
)
+

∑
n
=
-
L
+
1
-
1
a
n
k
w
(
2
j
x
-
n
)
,
k
=
0,1
,
2
,
…
,
L
,
w
(
2
j
x
-
k
)
,
k
=
L
+
1
,
…
,
i
i
i
i
2
j
-
L
-
1
,
w
(
2
j
x
-
k
)
+

∑
n
=
2
j
+
1
2
j
+
L
-
1
b
n
k
w
(
2
j
x
-
n
)
,
k
=
2
j
-

L
,
…
,
2
j
,
where
(22)
a
n
k
=
∏
i
=
L
-
1
i
≠
k
-
1
x
j
,
n
-
x
j
,
i
x
j
,
k
-
x
j
,
i
,
b
n
k
=
∏
i
=
2
j
+
1
i
≠
k
2
j
+
1
+
L
x
j
,
n
-
x
j
,
i
x
j
,
k
-
x
j
,
i
,
x
j
,
k
=
k
x
max
⁡
-
x
min
⁡
2
j
,
k
∈
ℤ
,whereL is the amount of the external collocation points, the amount of discrete points in the definition domain is 2
j
+
1
(
j
∈
Z
), and [
x
min
⁡
,
x
max
⁡
] is the definition domain of the approximated function.It is obvious that (21) is just the Lagrange extrapolation representation. The parameter L appeared in (21) is relative to the smoothness and gradient of the function to be approximated. In theory, the increase of L can improve the approximate precision. But it was pointed out in [23] that the increase of L can also bring the increase of the condition number of the problem to be solved, and this can decrease the numerical precision greatly. In addition, the larger L can also bring more calculation.Table1 shows the maximum absolute errors of the numerical approximation to curves with different gradient near the boundary corresponding to different value of L. It illustrates that L is not only related to the smoothness, but also to the gradient of original function near the boundary. For the same L, the greater the original function’s gradient is, the bigger the error is. To the different gradient original functions in Table 1, the reasonable L should be L
≥
3, L
≥
4, L
≥
4, and L
≥
8, respectively. In contrast, the interval wavelet proposed in this paper is not sensitive to the gradient of the function to be approximated, and the calculation amount is the same as the interval wavelet in [7] with L
=
1. As the matter of fact, the proposed interval wavelet is constructed based on that the derivative function should be continuous, and the smaller L avoids the Gibbs phenomenon appeared in the Lagrange interpolator. This is the reason why it is not sensitive to the gradient value.Table 1
Errors comparison of the numerical approximation to curves with different gradients (j
=
10, r
=
3).
L
sin(x)(10−15)
sin(5x)(10−15)
sin(10x)(10−15)
sin(100x)(10−14)
Interval wavelet in [7]
1
3964.88397669
123954.735364
793825.449818
2336702.80902
2
5.60662627435
634.825525480
1759.70349403
599685.701097
3
0.22204460492
2.22044604925
46.9624339416
18207.2357191
4
0.22204460492
0.44408920985
0.44408920985
3343.21459405
5
0.27755575615
0.44408920985
0.44408920985
129.352084599
6
0.22204460492
0.44408920985
0.44408920985
18.0300219199
7
0.27755575615
0.44408920985
0.44408920985
0.88817841970
8
0.27755575615
0.44408920985
0.44408920985
0.88817841970
New interval wavelet
0.27755575613
0.27755575625
0.27755575615
0.02784584512


## 3. Black-Scholes Equation and Its Interval Wavelet Approximation

The famous Black-Scholes equation can be represented as(23)
∂
V
∂
t
=
r
V
-
1
2
σ
2
S
2
∂
2
V
∂
S
2
-
r
S
∂
V
∂
S
,
where S(t) denotes the underlying asset, t
∈
(
0
,
T
), T denotes the expiry date, σ is the volatility (measures the standard deviation of the returns), and r is the riskless interest rate. In general, the parameter σ is taken as a constant since the transaction cost is taken as zero. Obviously, the σ is not really a constant, and then we can obtain the nonlinear Black-Scholes equation as follows:
(24)
∂
V
∂
t
=
r
V
-
1
2
σ
~
2
(
t
,
S
,
∂
V
∂
S
,
∂
2
V
∂
S
2
)
S
2
∂
2
V
∂
S
2
-
r
S
∂
V
∂
S
,
where σ
~ denotes a nonconstant volatility. In order to solve the problem, it is necessary to perform a variable transformation as follows:
(25)
x
=
ln
⁡
(
S
K
)
,
τ
=
1
2
σ
2
(
T
-
t
)
,
u
(
x
,
τ
)
=
e
-
x
V
(
s
,
t
)
K
.
Substituting (25) into (24), the following equation can be obtained:
(26)
∂
u
∂
t
=
σ
~
2
σ
2
(
∂
2
u
∂
x
2
+
∂
u
∂
x
)
+
D
∂
u
∂
x
,
where
(27)
D
=
2
r
σ
2
,
x
∈
R
,
0
≤
τ
≤
T
~
=
σ
2
2
.
Initial condition is
(28)
u
(
x
,
0
)
=
(
1
-
e
-
x
)
+
for
x
∈
R
.
And the boundary conditions are
(29)
u
(
x
,
τ
)
=
0
,
as
x
⟶
-
∞
,
u
(
x
,
τ
)
~
1
-
e
-
D
τ
-
x
,
as
x
⟶
∞
.
According to the transformation relation (25), it is easy to understand that the point x
=
0 is corresponding to the strike price S
=
K. Obviously, the initial solution curve is smooth except around the position x
=
0, where a sharp steep wave happened. So, an adaptive numerical method is suited to this problem.According to the wavelet precision integration method, the nonlinear Black-Scholes equation can be approximately represented as(30)
d
u
j
(
x
i
,
t
)
d
t
=
D
∑
m
=
0
2
j
u
j
(
x
m
,
t
n
)
[
w
′
(
x
i
-
x
m
)
]
+
σ
~
2
σ
2
∑
m
=
0
2
j
u
j
(
x
m
,
t
n
)
iiiiiiiiiiiiiiiiiiiii
×
[
w
′
(
x
i
-
x
m
)
+
w
′′
(
x
i
-
x
m
)
]
,
where i
=
0
,
1,2
,
…
,
2
j. Let
(31)
U
j
n
=
(
u
j
(
x
0
,
t
n
)
,
u
j
(
x
1
,
t
n
)
,
…
,
u
j
(
x
2
j
,
t
n
)
)
T
,
F
=
diag
⁡
(
f
′
(
x
0
)
,
f
′
(
x
1
)
,
…
,
f
′
(
x
2
j
)
)
,
W
1
=
[
w
′
(
x
0
-
x
0
)
w
′
(
x
0
-
x
1
)
⋯
w
′
(
x
0
-
x
2
j
)
w
′
(
x
1
-
x
0
)
w
′
(
x
1
-
x
1
)
⋯
w
′
(
x
1
-
x
2
j
)
⋮
⋮
⋱
⋮
w
′
(
x
2
j
-
x
0
)
w
′
(
x
2
j
-
x
1
)
⋯
w
′
(
x
2
j
-
x
2
j
)
]
,
W
2
=
[
w
′′
(
x
0
-
x
0
)
w
′′
(
x
0
-
x
1
)
⋯
w
′′
(
x
0
-
x
2
j
)
w
′′
(
x
1
-
x
0
)
w
′′
(
x
1
-
x
1
)
⋯
w
′′
(
x
1
-
x
2
j
)
⋮
⋮
⋱
⋮
w
′′
(
x
2
j
-
x
0
)
w
′′
(
x
2
j
-
x
1
)
⋯
w
′′
(
x
2
j
-
x
2
j
)
]
.
Then the system of (30) can be expressed as the matrix format:
(32)
d
U
j
n
d
t
=

D
·
W
1
U
j
n
+
σ
~
2
σ
2
(
W
1
+
W
2
)
U
j
n

=
M
0
U
j
n
+
M
1
(
U
j
n
)
U
j
n
.It is obvious thatM
0 is a constant matrix according to (32). M
1(V
j)  V
j is the nonlinear term in (32), which can be linearized subsectionally in terms of first-order Taylor series expansion; that is, the nonlinear terms can be taken as linear ones in the time step (t
k,t
k
+
1), and then (32) can be rewritten as
(33)
d
d
t
U
j
=
M
0
U
j
+
r
0
+
r
1
(
t
-
t
k
)
,
where
(34)
r
0
=
M
1
U
j
,
r
1
=
d
(
M
1
U
j
)
d
t
=
d
M
1
d
t
U
j
+
M
1
d
U
j
d
t
.
Letting M
2
=
d
M
1
/
d
t
=
-
diag
⁡
(
d
U
j
/
d
t
)
W
1 and combining with (32), the expression of r
1 can be obtained as
(35)
r
1
=
M
2
U
j
+
M
1
(
(
M
0
+
M
1
)
U
j
)
.The precise time integration scheme for the linear ordinary differential equations (33) can be expressed as
(36)
U
j
k
+
1
=
T
×
[
V
j
k
+
M
0
-
1
(
r
0
+
M
0
-
1
r
1
)
]
-
M
0
-
1
(
r
0
+
M
0
-
1
r
1
+
r
1
τ
)
,
where
(37)
T
=
exp
⁡
(
M
0
·
τ
)
(
τ
is time step
)
,
and then solving (33) is changed into the calculation of matrix T. the exponential matrix T can be calculated accurately by PIM as follows:
(38)
T
(
t
)
=
exp
⁡
⁡
(
M
0
t
)
=
[
exp
⁡
⁡
(
M
0
·
t
2
N
)
]
2
N
.
Let Δ
t
=
τ
/
2
N, where N is a positive integer (usually take N
=
20, and then Δ
t
=
τ
/
1048576). Taking into account that t is a small time interval and Δt is a very smaller value, then
(39)
exp
⁡
⁡
(
M
0
·
Δ
t
)
=
I
+
T
a
=
I
+
M
0
Δ
t
+
(
M
0
Δ
t
)
2
[
I
+
(
M
0
Δ
t
)
/
3
+
(
M
0
Δ
t
)
2
/
12
]
2
which is the Taylor series expansion of exp(M
0Δt). In order to calculate the matrix T more accurately, it is necessary to factorize the matrix T as
(40)
T
(
t
)
=
[
exp
⁡
⁡
(
M
0
·
Δ
t
)
]
2
N
=
(
I
+
T
a
)
2
N
-
1
(
I
+
T
a
)
2
N
-
1
.
After doing N times of factorization as above, a more accurate solution of T can be obtained.The calculation of the exponent matrixT
(
i
h
) at different time steps is needed in solving nonlinear equations through iteration based on the precise integration method, and the algorithm of the matrix T
(
i
h
) presented here can obtain all the matrices at different time steps for once.The numerical solution of the nonlinear Black-Scholes equation is shown in Figure1, where the riskless rate r
=
0.1, the volatility σ
=
0.2, and strike price K
=
100. With the increase of the parameter t, the smoothness of the solution function near the strike price K becomes better and better, and the amount of the collocation points is decreasing correspondingly. This illustrates the performance of the adaptive wavelet numerical method.Evolution of the call option price with the parametert.
(a)
(b)
(c)
(d)
(e)
(f)Compared with the wavelet precision integration method employed in Yan’s method [5], the amount of the collocation points in the interval wavelet precision integration method proposed in this paper is decreased greatly, especially around the boundary of the definition domain. This illustrates that the interval wavelet precision integration method can eliminate the boundary effect completely. It is obviously that the decrease of the collocation points can improve the efficiency of the algorithm. The performance of Mei’s interval wavelet [4] relates to the external interpolation points (L), and the choice scheme of L was not given in [3]. It was pointed out in [23] that the bigger L usually introduces bigger condition number which can decrease the precision greatly. This can be illustrated in Figures 1(d) and 1(e). As L
=
1, the boundary effect can be eliminated effectively. On the contrary, larger L almost cannot restrict it.In addition, the interval wavelet can also improve the numerical precision compared with the common wavelet method. For convenience, we take the linear Black-Scholes equation as an example, which has exact analytical solution as follows:(41)
C
=
S
·
N
(
d
1
)
-
K
e
-
r
T
N
(
d
2
)
,
where
(42)
d
1
=
ln
⁡
⁡
(
S
/
K
)
+
(
r
+
(
1
/
2
)
σ
2
)
T
σ
T
,
d
2
=
d
1
-
σ
T
.
C is the call price, S is the underlying asset price, K is the strike price, r is the riskless rate, T is the maturity, σ is the volatility, and N
(
d
1
) express the normal distribution.All the comparisons in Table2 are made qualitatively by comparing the calculation precision in the same time step and space mesh grid size. The first measure of error e
1 is given by
(43)
e
1
=
∥
U
j
n
-
U
exact
n
∥
∞
which provides a measure of the accuracy of the solution near the boundary. The second measure of error e
2 is given by
(44)
e
2
=
1
2
j
+
1
∑
i
=
0
2
j
(
u
(
x
i
)
-
u
exact
(
x
i
)
)
2
.
In Table 2, the parameter L is taken as 1 in Mei’s method. The error e
1 of Mei’s method almost equals that of Yan’s method, and the error e
2 of Mei’s method is smaller evidently than Yan’s method. This indicates that the interval wavelet method proposed by Mei can eliminate the boundary effect to some extent as L
=
1, but it cannot decrease the maximum of the absolute error evidently. In addition, with the increasing of the parameters j and t, the errors of both Yan and Mei’s method become larger and larger evidently. This can be explained by the condition number theory [23]. On the contrary, the interval wavelet method proposed in this paper is not sensitive to the parameters j and t, and the precision is the best among all above methods.Table 2
The numerical precision comparison.
j
t
e
1
e
2
Interval WPIM
Mei’s method
Yan’s method
Interval WPIM
Mei’s method
Yan’s method
12
10
5.5367
×
10
-
6
8.15
79
×
10
-
5
8.1588
×
10
-
5
4.1037
×
10
-
6
1.0237
×
10
-
5
5.8298
×
10
-
5
100
5.7907
×
10
-
6
4.1
246
×
10
-
4
4.1158
×
10
-
4
4.2424
×
10
-
6
0.6785
×
10
-
4
2.9468
×
10
-
4
200
3.5309
×
10
-
6
9.1
429
×
10
-
4
9.1673
×
10
-
4
2.6232
×
10
-
6
1.2321
×
10
-
4
6.3692
×
10
-
4
15
10
5.5551
×
10
-
6
8.2
327
×
10
-
5
8.2118
×
10
-
5
4.1718
×
10
-
6
1.0624
×
10
-
5
5.9148
×
10
-
5
100
5.7760
×
10
-
6
4
.
1025
×
10
-
4
4.0932
×
10
-
4
4.2907
×
10
-
6
0.6031
×
10
-
4
2.9753
×
10
-
4
200
6.5910
×
10
-
6
0.3812
0.3971
3
.
92
907
×
10
-
6
0.0211
0.1195


## 4. Conclusion

The solution of the nonlinear Black-Scholes equation has an evident sharp steep wave near the boundary, and this can introduce evident boundary effect in the wavelet approximation such as Yan’s method. The interval wavelet approximation proposed in this paper makes sure that the derivative function near the boundary of the definition domain is continued, which is helpful to eliminate the boundary effect and improve the numerical precision in the whole definition domain. Compared with Mei’s interval wavelet method, the interval wavelet proposed in this paper is not sensitive to the amount of the levels and the time step and does not have to choose the external collocation points. So, it is a simple and robust tool for solving the Black-Scholes equation. In fact, we can construct various interval wavelet functions aiming at different problems in engineering based on the general variational principle.



---

*Source: 541023-2014-02-20.xml*