# Accurate and Efficient Computations of the Greeks for Options Near Expiry Using the Black-Scholes Equations

Authors: Darae Jeong; Minhyun Yoo; Junseok Kim
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

We investigate the accurate computations for the Greeks using the numerical solutions of the Black-Scholes partial differential equation. In particular, we study the behaviors of the Greeks close to the maturity time and in the neighborhood around the strike price. The Black-Scholes equation is discretized using a nonuniform finite difference method. We propose a new adaptive time-stepping algorithm based on local truncation error. As a test problem for our numerical method, we consider a European cash-or-nothing call option. To show the effect of the adaptive stepping strategy, we calculate option price and its Greeks with various tolerances. Several numerical results confirm that the proposed method is fast, accurate, and practical in computing option price and the Greeks.


---

## Body

## 1. Introduction

In this paper, we investigate the accurate and efficient computations for the Greeks using the numerical solutions of the Black-Scholes (BS) partial differential equation (PDE) [1]. Let s
i and t denote the price of the underlying asset and time, respectively. Let s
=
(
s
1
,
s
2
,
…
,
s
n
) be the set of n underlying assets. Then, the value of an option u
(
s
,
t
) is governed by the following n-dimensional BS equation [2].For(
s
,
t
)
∈
R
+
d
×
[
0
,
T
), (1)
∂
u
s
,
t
∂
t
+
∑
i
=
1
d
r
s
i
∂
u
s
,
t
∂
s
i
+
1
2
∑
i
,
j
=
1
d
ρ
i
j
σ
i
σ
j
s
i
s
j
∂
2
u
s
,
t
∂
s
i
∂
s
j
-
r
u
s
,
t
=
0with a final condition u
(
s
,
T
)
=
Λ
(
s
), where r is the constant riskless interest rate, σ
i are volatility values of s
i, and ρ
i
j are the asset correlations between s
i and s
j. Λ
(
s
) is the payoff function at maturity T. There are three classical techniques which are the finite difference method (FDM) [3–13], the finite element method [14], and the finite volume method [15] for the numerical solutions of the BS PDE.Sensitivities of option price, the so-called Greeks of option values, are derivatives with respect to market variables or model parameters. In this paper, we will focus on the behaviors of the Greeks close to the maturity time and in the neighborhood of the strike price. The outline of the paper is as follows. In Section2, we present the numerical solution of the BS PDE. In Section 3, we show the results of the several numerical experiments and the conclusions are drawn in Section 4.


## 2. Numerical Solution

In this section, we present a numerical scheme and its solution for the one-dimensional BS equation. By introducingτ
=
T
-
t which means the time to expiry, the one-dimensional version of (1) becomes the following initial value problem: (2)
∂
u
x
,
τ
∂
τ
=
1
2
σ
x
2
∂
2
u
x
,
τ
∂
x
2
+
r
x
∂
u
∂
x
-
r
u
x
,
τ
,
for
x
,
τ
∈
Ω
×
0
,
Twith an initial condition u
(
x
,
0
) for x
∈
Ω
=
(
0
,
x
m
a
x
)
. Here, x
m
a
x is a sufficiently large asset price. The BS equation (2) is discretized on a grid defined by x
0
=
0 and x
i
+
1
=
x
i
+
h
i for i
=
0
,
…
,
N
x
-
1, where N
x is the number of grid intervals and h
i is the grid spacing (see Figure 1). We assume that x
N
x
=
x
m
a
x and the ghost point x
N
x
+
1
=
x
m
a
x
+
h
N
x
-
1.Figure 1
A nonuniform grid with spatial step sizesh
j
=
x
j
+
1
-
x
j.Letu
i
n
≈
u
(
x
i
,
n
Δ
τ
) be the numerical approximate solution, where Δ
τ
=
T
/
N
τ is the time step size and N
τ is the total number of time steps. By applying implicit Euler’s scheme to (2), we have (3)
u
i
n
+
1
-
u
i
n
Δ
τ
=
σ
2
x
i
2
2
2
u
i
-
1
n
+
1
h
i
-
1
h
i
-
1
+
h
i
-
2
u
i
n
+
1
h
i
-
1
h
i
+
2
u
i
+
1
n
+
1
h
i
h
i
-
1
+
h
i
+
r
x
i
-
h
i
u
i
-
1
n
+
1
h
i
-
1
h
i
-
1
+
h
i
+
h
i
-
h
i
-
1
u
i
n
+
1
h
i
-
1
h
i
+
h
i
-
1
u
i
+
1
n
+
1
h
i
h
i
-
1
+
h
i
-
r
u
i
n
+
1
,for i
=
1
,
…
,
N
x and n
=
0
,
…
,
N
τ
-
1. Then, we can rewrite (3) as (4)
r
x
i
h
i
-
σ
x
i
2
h
i
-
1
h
i
-
1
+
h
i
u
i
-
1
n
+
1
+
σ
x
i
2
-
r
x
i
h
i
-
h
i
-
1
h
i
-
1
h
i
+
r
+
1
Δ
τ
u
i
n
+
1
-
r
x
i
h
i
-
1
+
σ
x
i
2
h
i
h
i
-
1
+
h
i
u
i
+
1
n
+
1
=
u
i
n
Δ
τ
.In this paper, we restrict our attention to European call options. Since the option price at x
=
0 approaches zero, we impose the zero Dirichlet boundary condition as u
(
0
,
τ
)
=
0. Also, the option value at sufficiently large asset price is asymptotically linear. Therefore, we use the linear boundary condition [16–18] at x
=
x
m
a
x as ∂
2
u
/
∂
x
2
(
x
m
a
x
,
τ
)
=
0; that is, u
N
x
+
1
n
+
1
=
2
u
N
x
n
+
1
-
u
N
x
-
1
n
+
1. The linear system from (4) is solved by the Thomas algorithm [19].


### 2.1. Adaptive Time-Stepping Strategy

The numerical solutions of the BS PDE near maturity are very sensitive to the size of the time step used. In this paper, for the sake of efficiency and accuracy of the numerical solution, we consider an adaptive time-stepping strategy [20]. In this strategy, the time step is chosen by using criteria in accordance with a truncation error. Before we start, we consider the relation between the exact solution and the numerical approximation in terms of Δ
τ. We denote the exact solution for an advance from τ to τ
+
2
Δ
τ by u
(
x
,
τ
+
2
Δ
τ
) and the two approximate solutions by v (one step with 2
Δ
τ) and w (two Δ
τ steps). In this study, since we use the fully implicit scheme for time derivative, the numerical solution of (4) has a first-order accuracy with respect to time. For the numerical solution with one step by 2
Δ
τ, (5)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
2
Δ
τ
u
τ
x
,
τ
+
1
2
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
=
v
x
,
τ
+
2
Δ
τ
+
1
2
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
.And, for the numerical solutions with two steps by Δ
τ, (6)
u
x
,
τ
+
Δ
τ
=
u
x
,
τ
+
Δ
τ
u
τ
x
,
τ
+
1
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
,
(7)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
Δ
τ
+
Δ
τ
u
τ
x
,
τ
+
Δ
τ
+
1
2
Δ
τ
2
u
τ
τ
x
,
τ
+
Δ
τ
+
O
Δ
τ
3
.By adding (6) and (7), we obtain(8)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
Δ
τ
u
τ
x
,
τ
+
u
τ
x
,
τ
+
Δ
τ
+
Δ
τ
2
u
τ
τ
x
,
τ
+
u
τ
τ
x
,
τ
+
Δ
τ
2
+
O
Δ
τ
3
=
w
x
,
τ
+
2
Δ
τ
+
Δ
τ
2
u
τ
τ
x
,
τ
+
u
τ
τ
x
,
τ
+
Δ
τ
2
+
O
Δ
τ
3
.Let E be the difference between the two numerical estimates; that is, (9)
E
=
w
x
,
τ
+
2
Δ
τ
-
v
x
,
τ
+
2
Δ
τ
=
Δ
τ
2
ϕ
+
O
Δ
τ
3
,where ϕ means the constant value whose order of magnitude is u
τ
τ
(
x
,
τ
). As shown in Figure 2, there is a difference of two numerical solutions by one step with 2
Δ
τ and two steps with Δ
τ. In general, the numerical approximation by two steps with Δ
τ is more accurate than by one step with 2
Δ
τ.Figure 2
Schematic illustration of the numerical solutionsw
(
x
,
τ
) and v
(
x
,
τ
) by two different time step sizes in Euler scheme. Here, w
(
x
,
τ
) denotes the numerical solution by one step with 2
Δ
τ and v
(
τ
) denotes the numerical solution by two steps with Δ
τ.In this study, we apply the adaptive time step strategy which is based on the local truncation error. First, we set the maximum and minimum time step sizes to avoid using too large or small time step,Δ
τ
m
a
x and Δ
τ
m
i
n. Next, let Δ
τ
0 be an initial time step size. With a given numerical solution u
n and a time step Δ
τ
=
Δ
τ
0, we solve (4) twice to get u
n
+
2. Next, with the given numerical solution u
n and a twice larger time step Δ
τ
=
2
Δ
τ
0, we solve (4) to get v
n
+
2. We define the time step scaled error as (10)
E
t
r
=
u
n
+
2
-
v
n
+
2
2
,where u
n
+
2 and v
n
+
2 are the numerical solution with times Δ
τ
0 and 2
Δ
τ
0, respectively. If the error is below the given tolerance, then we set the (
n
+
2
)th numerical solution as u
n
+
2. Otherwise, we solve (4) using each time step Δ
τ
0
/
2 and then check the scaled error. This process repeats until the scaled error meets the given tolerance. In our strategy, the next time step size is automatically determined by the given tolerance t
o
l and the error E
t
r as Δ
τ
n
e
w
=
Δ
τ
0
×
t
o
l
/
E
t
r. If E
t
r
<
t
o
l, that is, t
o
l
/
E
t
r
>
1, the new time step size is larger than the old one. Otherwise, which is E
t
r
>
t
o
l, the new one is smaller than the old one.The adaptive time-stepping strategy can be summarized in Algorithm1.Algorithm 1:Algorithm of adaptive time-stepping strategy.
Require: Set the initial conditions u
0 and v
0, the expiry time T, the maximum and minimum
time stepsΔ
τ
m
a
x
=
T and Δ
τ
m
i
n
=
T
/
86400, tolerance tol, truncation error E
t
r, safety
factors
=
0.8, time τ
=
0, iteration number n
=
0.
While
τ
<
T
do
n
=
n
+
1
if
n
=
1
then
Set Δ
τ
0
=
t
o
l
×
Δ
τ
m
a
x.
else
Set Δ
τ
0
=
s
×
Δ
τ
0
×
t
o
l
/
E
t
r.
end if
if
τ
+
Δ
τ
0
>
T
then
Set Δ
τ
0
=
T
-
τ.
Solve (4) with time step Δ
τ
0 and set its approximation as u
n.
end if
Solve (4) twice with time step Δ
τ
0 and set its approximation as u
∗.
Solve (4) with time step 2
Δ
τ
0 and set its approximation as v
∗.
Calculate E
t
r
=
∑
x
i
∈
Ω
-
(
u
i
∗
-
v
i
∗
)
2
/
M where Ω
-
=
[
0.9
K
,
1.1
K
] and M is the total number of x
i
∈
Ω
-.
while
E
t
r
>
t
o
l
do
Set Δ
τ
0
=
Δ
τ
0
/
2.
if
Δ
τ
0
<
Δ
τ
m
i
n
then
Set Δ
τ
0
=
2
Δ
τ
0.
end if
Solve (4) twice with Δ
τ
0 and set its approximation as u
∗.
Solve (4) with 2
Δ
τ
0 and set its approximation as v
∗.
Calculate E
t
r
=
∑
x
i
∈
Ω
-
(
u
i
∗
-
v
i
∗
)
2
/
M
.
end while
Set τ
=
τ
+
Δ
τ
0.
Set u
n
=
u
∗ and v
n
=
u
∗.
end while


## 2.1. Adaptive Time-Stepping Strategy

The numerical solutions of the BS PDE near maturity are very sensitive to the size of the time step used. In this paper, for the sake of efficiency and accuracy of the numerical solution, we consider an adaptive time-stepping strategy [20]. In this strategy, the time step is chosen by using criteria in accordance with a truncation error. Before we start, we consider the relation between the exact solution and the numerical approximation in terms of Δ
τ. We denote the exact solution for an advance from τ to τ
+
2
Δ
τ by u
(
x
,
τ
+
2
Δ
τ
) and the two approximate solutions by v (one step with 2
Δ
τ) and w (two Δ
τ steps). In this study, since we use the fully implicit scheme for time derivative, the numerical solution of (4) has a first-order accuracy with respect to time. For the numerical solution with one step by 2
Δ
τ, (5)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
2
Δ
τ
u
τ
x
,
τ
+
1
2
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
=
v
x
,
τ
+
2
Δ
τ
+
1
2
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
.And, for the numerical solutions with two steps by Δ
τ, (6)
u
x
,
τ
+
Δ
τ
=
u
x
,
τ
+
Δ
τ
u
τ
x
,
τ
+
1
2
Δ
τ
2
u
τ
τ
x
,
τ
+
O
Δ
τ
3
,
(7)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
Δ
τ
+
Δ
τ
u
τ
x
,
τ
+
Δ
τ
+
1
2
Δ
τ
2
u
τ
τ
x
,
τ
+
Δ
τ
+
O
Δ
τ
3
.By adding (6) and (7), we obtain(8)
u
x
,
τ
+
2
Δ
τ
=
u
x
,
τ
+
Δ
τ
u
τ
x
,
τ
+
u
τ
x
,
τ
+
Δ
τ
+
Δ
τ
2
u
τ
τ
x
,
τ
+
u
τ
τ
x
,
τ
+
Δ
τ
2
+
O
Δ
τ
3
=
w
x
,
τ
+
2
Δ
τ
+
Δ
τ
2
u
τ
τ
x
,
τ
+
u
τ
τ
x
,
τ
+
Δ
τ
2
+
O
Δ
τ
3
.Let E be the difference between the two numerical estimates; that is, (9)
E
=
w
x
,
τ
+
2
Δ
τ
-
v
x
,
τ
+
2
Δ
τ
=
Δ
τ
2
ϕ
+
O
Δ
τ
3
,where ϕ means the constant value whose order of magnitude is u
τ
τ
(
x
,
τ
). As shown in Figure 2, there is a difference of two numerical solutions by one step with 2
Δ
τ and two steps with Δ
τ. In general, the numerical approximation by two steps with Δ
τ is more accurate than by one step with 2
Δ
τ.Figure 2
Schematic illustration of the numerical solutionsw
(
x
,
τ
) and v
(
x
,
τ
) by two different time step sizes in Euler scheme. Here, w
(
x
,
τ
) denotes the numerical solution by one step with 2
Δ
τ and v
(
τ
) denotes the numerical solution by two steps with Δ
τ.In this study, we apply the adaptive time step strategy which is based on the local truncation error. First, we set the maximum and minimum time step sizes to avoid using too large or small time step,Δ
τ
m
a
x and Δ
τ
m
i
n. Next, let Δ
τ
0 be an initial time step size. With a given numerical solution u
n and a time step Δ
τ
=
Δ
τ
0, we solve (4) twice to get u
n
+
2. Next, with the given numerical solution u
n and a twice larger time step Δ
τ
=
2
Δ
τ
0, we solve (4) to get v
n
+
2. We define the time step scaled error as (10)
E
t
r
=
u
n
+
2
-
v
n
+
2
2
,where u
n
+
2 and v
n
+
2 are the numerical solution with times Δ
τ
0 and 2
Δ
τ
0, respectively. If the error is below the given tolerance, then we set the (
n
+
2
)th numerical solution as u
n
+
2. Otherwise, we solve (4) using each time step Δ
τ
0
/
2 and then check the scaled error. This process repeats until the scaled error meets the given tolerance. In our strategy, the next time step size is automatically determined by the given tolerance t
o
l and the error E
t
r as Δ
τ
n
e
w
=
Δ
τ
0
×
t
o
l
/
E
t
r. If E
t
r
<
t
o
l, that is, t
o
l
/
E
t
r
>
1, the new time step size is larger than the old one. Otherwise, which is E
t
r
>
t
o
l, the new one is smaller than the old one.The adaptive time-stepping strategy can be summarized in Algorithm1.Algorithm 1:Algorithm of adaptive time-stepping strategy.
Require: Set the initial conditions u
0 and v
0, the expiry time T, the maximum and minimum
time stepsΔ
τ
m
a
x
=
T and Δ
τ
m
i
n
=
T
/
86400, tolerance tol, truncation error E
t
r, safety
factors
=
0.8, time τ
=
0, iteration number n
=
0.
While
τ
<
T
do
n
=
n
+
1
if
n
=
1
then
Set Δ
τ
0
=
t
o
l
×
Δ
τ
m
a
x.
else
Set Δ
τ
0
=
s
×
Δ
τ
0
×
t
o
l
/
E
t
r.
end if
if
τ
+
Δ
τ
0
>
T
then
Set Δ
τ
0
=
T
-
τ.
Solve (4) with time step Δ
τ
0 and set its approximation as u
n.
end if
Solve (4) twice with time step Δ
τ
0 and set its approximation as u
∗.
Solve (4) with time step 2
Δ
τ
0 and set its approximation as v
∗.
Calculate E
t
r
=
∑
x
i
∈
Ω
-
(
u
i
∗
-
v
i
∗
)
2
/
M where Ω
-
=
[
0.9
K
,
1.1
K
] and M is the total number of x
i
∈
Ω
-.
while
E
t
r
>
t
o
l
do
Set Δ
τ
0
=
Δ
τ
0
/
2.
if
Δ
τ
0
<
Δ
τ
m
i
n
then
Set Δ
τ
0
=
2
Δ
τ
0.
end if
Solve (4) twice with Δ
τ
0 and set its approximation as u
∗.
Solve (4) with 2
Δ
τ
0 and set its approximation as v
∗.
Calculate E
t
r
=
∑
x
i
∈
Ω
-
(
u
i
∗
-
v
i
∗
)
2
/
M
.
end while
Set τ
=
τ
+
Δ
τ
0.
Set u
n
=
u
∗ and v
n
=
u
∗.
end while


## 3. Numerical Experiments

In this section, for numerical experiments, we consider a European cash-or-nothing option which pays an amountC at maturity if option is in-the-money state. The payoff function is given by (11)
u
x
,
0
=
C
if
x
>
K
,
0
o
t
h
e
r
w
i
s
e
,where K is the strike price and C denotes the return value. The closed-form solution [21] for the option is given as (12)
u
x
,
τ
=
C
e
-
r
τ
N
d
,
for
∀
x
,
τ
∈
0
,
L
×
0
,
T
,where d
=
[
ln
⁡
x
/
K
+
r
-
0.5
σ
2
τ
]
/
(
σ
τ
) and N
(
d
)
=
1
/
2
π
∫
-
∞
d
‍
exp
⁡
-
x
2
/
2
d
x is the cumulative distribution function for the standard normal distribution [1]. In the Appendix, MATLAB code for the closed-form solution is presented.Figure3 shows the cash-or-nothing option prices at τ
=
0, τ
=
1
/
365, and τ
=
2
/
365 on Ω
=
[
0,300
]. Here, the option prices at τ
=
1
/
365 and τ
=
2
/
365 are obtained by (12). And we use strike price K
=
100, cash C
=
100, the risk-free interest rate r
=
0.03, and volatility σ
=
0.3. As shown in Figure 3, the option price has dropped drastically for the first time step. Therefore, we need to take smaller time step sizes in early times since the solution rapidly changes.Figure 3
Valuesu
(
x
,
τ
) of cash-or-nothing option at τ
=
0, τ
=
1
/
365, and τ
=
2
/
365.In the following sections, unless otherwise specified, we use strike priceK
=
100, cash C
=
100, the risk-free interest rate r
=
0.03, the volatility σ
=
0.3, and L
=
250. All computations are performed using MATLAB version 8 [22].


### 3.1. Convergence Test

First, we present the performance of the numerical scheme with respect to uniform spatial and temporal step sizes. For measurement of accuracy of the numerical scheme, we compute the absolute errore
=
u
e
x
a
c
t
-
u, where u
e
x
a
c
t and u denote the exact and the numerical solutions for cash-or-nothing option, respectively. For consistent comparison of the accuracy, we evaluate the absolute error at x
=
100.5 that is the position near the predetermined strike price K
=
100. Figure 4 shows the spatial grid structure near K with respect to h. And the marked circle (x
=
100.5) denotes the point where the absolute error is measured.Figure 4
Construction of the spatial grid near the strike priceK
=
100 with respect to h. Here, the marked circle denotes the point where the absolute error is measured.In this study, we are concerned with the option value and its Greeks near the expiry. Therefore, we only consider the cash-or-nothing option at the day before the option is expired; that is,τ
=
1
/
365. Table 1 represents the absolute error between the numerical and exact solutions with respect to spatial and temporal step sizes.Table 1
Comparison of absolute error|
u
e
x
a
c
t
-
u
| between the numerical and exact solutions at x
=
100.5 and τ
=
1
/
365 with respect to h and Δ
τ.
h
Δ
τ
1
/
365
1
/
730
1
/
1460
1
/
2920
1
/
5840
1
/
11680
1
8.638266
5.403661
3.560908
2.618485
2.155975
1.929931
1/3
6.541975
3.457916
1.959400
1.288850
0.979618
0.831513
1
/
9
6.276309
3.262838
1.829017
1.187406
0.889680
0.746603
1
/
27
6.246320
3.241839
1.815166
1.176523
0.879980
0.737425
1
/
81
6.242982
3.239515
1.813635
1.175318
0.878906
0.736408As shown in Table1, the absolute error decreases as h and Δ
τ decrease. Also, we can see these results in Figure 5.Figure 5
Absolute errors between numerical and exact solutions onx
=
100.5 and τ
=
1
/
365 with respect to spatial and temporal step sizes.In the results, we confirm that the absolute error is more influenced by temporal step sizeΔ
τ than by spatial step size h. In particular, the numerical results with h
≤
1
/
3 are sufficiently accurate. Therefore, if we control time step size Δ
τ, we can obtain a more accurate numerical value. In the next section, we will present several numerical results with h
=
1
/
3 when we use the adaptive time-stepping strategy.


### 3.2. Adaptive Time-Stepping Scheme

Now, we consider an adaptive time-stepping strategy to efficiently solve the BS PDE. In this strategy, the time step size is determined at everyτ
n along with the truncation error of the numerical solution as we described before.As the default parameters, we set the initial time step sizeΔ
τ
0
=
t
o
l
×
Δ
τ
m
a
x, where Δ
τ
m
a
x
=
1
/
365 and Δ
τ
m
i
n
=
Δ
τ
m
a
x
/
(
24
·
3600
). Here, Δ
τ
m
a
x
=
1
/
365 and Δ
τ
m
i
n
=
1
/
(
365
·
24
·
3600
) represent one day and one second, respectively. Also, we use the spatial step size h
=
1
/
3 according to results in Table 1.Table2 represents the root mean square error (RMSE) of the numerical solution which is measured in the area [
0.9
K
,
1.1
K
]. Here, RMSE is calculated by ∑
i
=
1
N
(
u
i
-
u
i
e
x
)
2
/
N, where u and u
e
x denote the numerical and exact solutions. In Table 2, min
⁡
Δ
τ
0 and max
⁡
Δ
τ
0 represent the minimum and maximum time step sizes which are used during the numerical iteration. And N
τ denotes the total number of iterations during one day, T
=
1
/
365. Note that the RMSE by only one time step size Δ
τ
=
1
/
365 is 2.20437414. For comparison, we evaluate the ratio of 2.20437414 to the RMSE by adaptive time-stepping strategy. As shown in Table 2, the numerical solution by the adaptive time-stepping method is about 40 times more accurate than that by the one time step.Table 2
RMSE by adaptive time stepping strategy with respect to different tolerance.
Tolerance
min
⁡
Δ
τ
0
max
⁡
Δ
τ
0
N
τ
RMSE
Ratio
1.0
E
-
2
5.94366
E
-
6
2.04576
E
-
4
46
0.321825625
6.89882
1.0
E
-
3
1.28176
E
-
6
4.77865
E
-
5
213
0.067703445
32.79326
1.0
E
-
4
5.47945
E
-
7
2.1621
E
-
5
480
0.060855408
36.48347
1.0
E
-
5
5.47945
E
-
8
4.82658
E
-
6
2150
0.056760465
39.11554
1.0
E
-
6
5.47945
E
-
9
2.15918
E
-
6
4801
0.052404827
42.36664
1.0
E
-
7
5.47945
E
-
10
6.82985
E
-
7
14214
0.054571982
40.68418
1.0
E
-
8
5.47945
E
-
11
2.16019
E
-
7
39776
0.055265093
40.17394Figure6 shows the time step size Δ
τ against time τ over one day, T
=
1
/
365, in adaptive time-stepping strategy with t
o
l
=
1.0
E
-
4. In Figure 6, early adaptive time step is very small because the errors between numerical and exact solutions are greatly generated during the total time T
=
1
/
365.Figure 6
Graph of time step sizeΔ
τ against time τ with tolerance 1.0
E
-
4.Next, we perform numerical tests around maturity and compare the numerical solution with the exact solution. Let us defineu
e
x and u as the analytic solution by (12) and numerical solution by adaptive time strategy, respectively. Figure 7(a) shows that the usage of time step as Δ
τ
=
1
/
365 near maturity makes inaccurate solution. Figure 7(b) illustrates u
-
u
e
x with respect to the changes of t
o
l.Figure 7
(a) Option pricesu
(
x
,
1
/
365
) as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors u
-
u
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)We note that as thet
o
l decreases, the corresponding numerical solution remains accuracy, however, for a longer time. Therefore, the tolerance, t
o
l, can be determined by considering speed and accuracy trade-off.


### 3.3. The Greeks

The Greeks are defined as changes in option value relative to changes in each independent variable. For example, Delta is the first derivative of the option value with respect to the underlying asset. Gamma is the second derivative of the value with respect to the asset. Theta is the first derivative of the option with respect to time. Vega is the option’s sensitivity to changes in the volatility. Rho is the first derivative of the option price with respect to interest rates. For more details about the Greeks, we refer the reader to [21].In general, the payoff structure of the option is very stiff because it has a discontinuity around the strike price. Therefore, there exists a difficulty to numerically measure an accurate value of the Greeks. However, as we use the proposed adaptive time-stepping strategy, we can control the accuracy of the Greeks.Next, we evaluate the option sensitivities of the cash-or-nothing option atx
=
100. In particular, because of the given grid points as shown in Figure 1, we use the linear interpolation to calculate the Greeks at x
=
100 in some numerical tests. For comparison, we also evaluate the Greeks by the exact formula using MATLAB code which is presented in the Appendix.


#### 3.3.1. Delta

Delta (Δ) is defined as the rate of change of the option value with respect to small changes in the underlying asset price. By differentiating (12), we have the following exact formula: (13)
Δ
e
x
=
∂
u
∂
x
=
C
e
-
r
τ
N
′
d
σ
x
τ
.Also, we can define the Delta by using the finite difference discretizations at x
i. For example, (14)
Δ
i
n
≔
∂
u
∂
x
i
n
=
-
h
i
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
+
h
i
-
h
i
-
1
u
i
n
h
i
-
1
h
i
+
h
i
-
1
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 8(a) shows the Delta by closed-form formula (13), numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. To represent the accuracy of the Delta with respect to t
o
l, we investigate the error between closed-form solution and numerical approximation by our proposed method in Figure 8(b). Through these results, we can see that the adaptive time step strategy reduces the numerical error generated around the strike price K
=
100.Figure 8
(a) DeltaΔ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Δ
-
Δ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table3 shows the Delta and its error by the adaptive time step strategy with respect to t
o
l. As shown in Table 3, the Delta by our proposed method is more accurate as t
o
l is smaller.Table 3
DeltaΔ and its error Δ
-
Δ
e
x at x
=
100 and T
=
1
/
365. Here, Delta by the exact formula is Δ
e
x
=
25.40376.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Δ
24.7890
25.5980
25.5713
25.5553
25.5425
25.5481
Δ
-
Δ
e
x
-
0.61472
0.19425
0.16751
0.15149
0.13876
0.14432


#### 3.3.2. Gamma

Gamma (Γ) is the rate of change of the Delta with respect to changes in the underlying asset price. In other words, Gamma is defined by the second partial derivative of the portfolio value with respect to underlying asset price. According to closed-form solution (12), we get the exact formula for Gamma as (15)
Γ
e
x
=
∂
2
u
∂
x
2
=
-
C
e
-
r
τ
d
+
σ
τ
N
′
d
σ
2
x
2
τ
.As a different approach, Gamma is described by applying numerical discretization at x
i as follows: (16)
Γ
i
n
≔
∂
2
u
∂
x
2
i
n
=
2
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
-
2
u
i
n
h
i
-
1
h
i
+
2
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 9 shows the Gamma at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. As shown in Figure 9, while the numerical solution by one time step Δ
τ
=
1
/
365 has a great deviation from the exact formula (15), numerical approximations by adaptive time step strategy and exact solution of Gamma are in good agreement. Also, from Figure 9(b), we observe that the accuracy of the Gamma can be controlled by t
o
l in our proposed method.Figure 9
(a) GammaΓ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Γ
-
Γ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table4 displays the numerical value by our proposed method and its corresponding errors of the Gamma at x
=
100. As tolerance is smaller, we can obtain a convergent value of Gamma.Table 4
GammaΓ and its error Γ
-
Γ
e
x at x
=
100 and T
=
1
/
365. Here, Γ
e
x
=
-
0.21170.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Γ
-
0.19796
-
0.20396
-
0.20378
-
0.20368
-
0.20359
-
0.20363
Γ
-
Γ
e
x
0.01374
0.00774
0.00792
0.00802
0.00811
0.00807


#### 3.3.3. Theta

Theta (Θ) is the rate of change of the option value with respect to changes in the time to maturity, whose exact formula on cash-or-nothing option is given as (17)
Θ
e
x
=
∂
u
∂
t
=
-
∂
u
∂
τ
=
C
e
-
r
τ
r
N
d
+
N
′
d
d
+
σ
τ
2
τ
-
r
σ
τ
. Also, the numerical Theta is described as (18)
Θ
i
n
≔
∂
u
∂
τ
i
n
=
u
i
n
+
1
-
u
i
n
-
1
2
Δ
τ
. Here, we calculate the Theta by using the central difference approximation as (19)
Θ
i
=
u
x
i
,
T
+
Δ
τ
end
-
u
x
i
,
T
-
Δ
τ
end
2
Δ
τ
end
,where Δ
τ
end
=
0.01
T.In Figure10(a), we represent Θ by each one of the methods. As a result, numerical results by the adaptive time algorithm are superior to the other ones. Also, we can see that numerical results of Θ converge as t
o
l is smaller (see Figure 10(b) and Table 5).Table 5
ThetaΘ and its error Θ
-
Θ
e
x at x
=
100 and T
=
1
/
365. Here, Θ
e
x
=
20.54956.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Θ
19.41235
19.96587
19.88387
19.94125
19.93973
19.93678
Θ
-
Θ
e
x
-
1.13721
-
0.58369
-
0.66570
-
0.60832
-
0.60984
-
0.61279Figure 10
(a) ThetaΘ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Θ
-
Θ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


#### 3.3.4. Vega

Vega (ν) is the rate of the option value with respect to small changes in the volatility of the underlying asset. By the closed-form solution of cash-or-nothing option, (20)
ν
e
x
=
∂
u
∂
σ
=
C
e
-
r
τ
d
σ
+
τ
N
′
d
. We can write the Vega by using the discretizations as (21)
ν
i
n
≔
∂
u
∂
σ
i
n
=
u
i
n
σ
+
Δ
σ
-
u
i
n
σ
-
Δ
σ
2
Δ
σ
,where Δ
σ
=
0.01.As shown in Figure11(a), we obtain a reasonable approximation for the Vega ν when we use the adaptive time step strategy. Moreover, the difference of ν and ν
e
x decays as smaller tolerance level t
o
l is determined.Figure 11
(a) Vegaν at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ν
-
ν
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)The results in Table6 show that our proposed method could be improved taking a lower tolerance parameter t
o
l.Table 6
Vegaν and its error ν
-
ν
e
x at x
=
100 and T
=
1
/
365. Here, ν
e
x
=
-
1.73998.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ν
-
1.72129
-
1.72138
-
1.72165
-
1.71663
-
1.72054
-
1.72260
ν
-
ν
e
x
0.01870
0.01860
0.01833
0.02335
0.01945
0.01739


#### 3.3.5. Rho

Rho (ρ) is the rate of change of the option value with respect to small changes in the riskless interest rate r. The exact Rho formula on cash-or-nothing option is derived from (12) as (22)
ρ
e
x
=
∂
u
∂
r
=
C
e
-
r
τ
-
τ
N
d
+
τ
σ
N
′
d
. Rho is evaluated by numerical discretization as follows: (23)
ρ
i
n
≔
∂
u
∂
r
i
n
=
u
i
n
r
+
Δ
r
-
u
i
n
r
-
Δ
r
2
Δ
r
,where Δ
r
=
0.001.Similar to previous tests, we obtain good results for Rho (ρ) by the adaptive time step algorithm as shown in Figure 12 and Table 7. Particularly, when we compare results by our proposed method with one time step simulation, we see the superiority in terms of accuracy.Table 7
Rhoρ and its error ρ
-
ρ
e
x at x
=
100 and T
=
1
/
365. Here, ρ
e
x
=
6.82325.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ρ
7.02722
6.77944
6.78185
6.78329
6.78617
6.78449
ρ
-
ρ
e
x
0.20394
-
0.04384
-
0.04143
-
0.03999
-
0.03711
-
0.03879Figure 12
(a) Rhoρ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ρ
-
ρ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


### 3.4. Comparison with Adaptive and Uniform Time Step

In this section, we investigate the difference of numerical solutions with adaptive and uniform time step size. For comparison, we find the numerical results by uniform time step, which have similar RMSE to the results by the adaptive time step method. Table8 shows RMSE, N
τ, and CPU time with adaptive and uniform time step methods. As shown in Table 8, the uniform time step method is more efficient than the adaptive time step method at large tolerance (t
o
l
≥
1.0
E
-
5). However, we can see that when tolerance t
o
l is below 1.0
E
-
6, the results by adaptive time step are better than those by uniform time step. Furthermore, we obtain more accurate results by adaptive time method than by uniform method.Table 8
Comparison with adaptive and uniform time step methods.
tol
RMSE
N
τ
CPU time
Adaptive
1.0
E
-
2
0.32183
46
0.02790
Uniform
0.33442
9
0.00415
Adaptive
1.0
E
-
3
0.06770
213
0.11714
Uniform
0.06774
207
0.02237
Adaptive
1.0
E
-
4
0.06086
480
0.15395
Uniform
0.06086
476
0.04640
Adaptive
1.0
E
-
5
0.05676
2150
1.04035
Uniform
0.05676
2144
0.20096
Adaptive
1.0
E
-
6
0.05240
4801
1.41838
Uniform
0.05559
10000000
902.92370
Adaptive
1.0
E
-
7
0.05457
14214
4.24171
Uniform
0.05559
10000000
902.92370
Adaptive
1.0
E
-
8
0.05527
39776
12.89830
Uniform
0.05559
10000000
902.92370


## 3.1. Convergence Test

First, we present the performance of the numerical scheme with respect to uniform spatial and temporal step sizes. For measurement of accuracy of the numerical scheme, we compute the absolute errore
=
u
e
x
a
c
t
-
u, where u
e
x
a
c
t and u denote the exact and the numerical solutions for cash-or-nothing option, respectively. For consistent comparison of the accuracy, we evaluate the absolute error at x
=
100.5 that is the position near the predetermined strike price K
=
100. Figure 4 shows the spatial grid structure near K with respect to h. And the marked circle (x
=
100.5) denotes the point where the absolute error is measured.Figure 4
Construction of the spatial grid near the strike priceK
=
100 with respect to h. Here, the marked circle denotes the point where the absolute error is measured.In this study, we are concerned with the option value and its Greeks near the expiry. Therefore, we only consider the cash-or-nothing option at the day before the option is expired; that is,τ
=
1
/
365. Table 1 represents the absolute error between the numerical and exact solutions with respect to spatial and temporal step sizes.Table 1
Comparison of absolute error|
u
e
x
a
c
t
-
u
| between the numerical and exact solutions at x
=
100.5 and τ
=
1
/
365 with respect to h and Δ
τ.
h
Δ
τ
1
/
365
1
/
730
1
/
1460
1
/
2920
1
/
5840
1
/
11680
1
8.638266
5.403661
3.560908
2.618485
2.155975
1.929931
1/3
6.541975
3.457916
1.959400
1.288850
0.979618
0.831513
1
/
9
6.276309
3.262838
1.829017
1.187406
0.889680
0.746603
1
/
27
6.246320
3.241839
1.815166
1.176523
0.879980
0.737425
1
/
81
6.242982
3.239515
1.813635
1.175318
0.878906
0.736408As shown in Table1, the absolute error decreases as h and Δ
τ decrease. Also, we can see these results in Figure 5.Figure 5
Absolute errors between numerical and exact solutions onx
=
100.5 and τ
=
1
/
365 with respect to spatial and temporal step sizes.In the results, we confirm that the absolute error is more influenced by temporal step sizeΔ
τ than by spatial step size h. In particular, the numerical results with h
≤
1
/
3 are sufficiently accurate. Therefore, if we control time step size Δ
τ, we can obtain a more accurate numerical value. In the next section, we will present several numerical results with h
=
1
/
3 when we use the adaptive time-stepping strategy.


## 3.2. Adaptive Time-Stepping Scheme

Now, we consider an adaptive time-stepping strategy to efficiently solve the BS PDE. In this strategy, the time step size is determined at everyτ
n along with the truncation error of the numerical solution as we described before.As the default parameters, we set the initial time step sizeΔ
τ
0
=
t
o
l
×
Δ
τ
m
a
x, where Δ
τ
m
a
x
=
1
/
365 and Δ
τ
m
i
n
=
Δ
τ
m
a
x
/
(
24
·
3600
). Here, Δ
τ
m
a
x
=
1
/
365 and Δ
τ
m
i
n
=
1
/
(
365
·
24
·
3600
) represent one day and one second, respectively. Also, we use the spatial step size h
=
1
/
3 according to results in Table 1.Table2 represents the root mean square error (RMSE) of the numerical solution which is measured in the area [
0.9
K
,
1.1
K
]. Here, RMSE is calculated by ∑
i
=
1
N
(
u
i
-
u
i
e
x
)
2
/
N, where u and u
e
x denote the numerical and exact solutions. In Table 2, min
⁡
Δ
τ
0 and max
⁡
Δ
τ
0 represent the minimum and maximum time step sizes which are used during the numerical iteration. And N
τ denotes the total number of iterations during one day, T
=
1
/
365. Note that the RMSE by only one time step size Δ
τ
=
1
/
365 is 2.20437414. For comparison, we evaluate the ratio of 2.20437414 to the RMSE by adaptive time-stepping strategy. As shown in Table 2, the numerical solution by the adaptive time-stepping method is about 40 times more accurate than that by the one time step.Table 2
RMSE by adaptive time stepping strategy with respect to different tolerance.
Tolerance
min
⁡
Δ
τ
0
max
⁡
Δ
τ
0
N
τ
RMSE
Ratio
1.0
E
-
2
5.94366
E
-
6
2.04576
E
-
4
46
0.321825625
6.89882
1.0
E
-
3
1.28176
E
-
6
4.77865
E
-
5
213
0.067703445
32.79326
1.0
E
-
4
5.47945
E
-
7
2.1621
E
-
5
480
0.060855408
36.48347
1.0
E
-
5
5.47945
E
-
8
4.82658
E
-
6
2150
0.056760465
39.11554
1.0
E
-
6
5.47945
E
-
9
2.15918
E
-
6
4801
0.052404827
42.36664
1.0
E
-
7
5.47945
E
-
10
6.82985
E
-
7
14214
0.054571982
40.68418
1.0
E
-
8
5.47945
E
-
11
2.16019
E
-
7
39776
0.055265093
40.17394Figure6 shows the time step size Δ
τ against time τ over one day, T
=
1
/
365, in adaptive time-stepping strategy with t
o
l
=
1.0
E
-
4. In Figure 6, early adaptive time step is very small because the errors between numerical and exact solutions are greatly generated during the total time T
=
1
/
365.Figure 6
Graph of time step sizeΔ
τ against time τ with tolerance 1.0
E
-
4.Next, we perform numerical tests around maturity and compare the numerical solution with the exact solution. Let us defineu
e
x and u as the analytic solution by (12) and numerical solution by adaptive time strategy, respectively. Figure 7(a) shows that the usage of time step as Δ
τ
=
1
/
365 near maturity makes inaccurate solution. Figure 7(b) illustrates u
-
u
e
x with respect to the changes of t
o
l.Figure 7
(a) Option pricesu
(
x
,
1
/
365
) as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors u
-
u
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)We note that as thet
o
l decreases, the corresponding numerical solution remains accuracy, however, for a longer time. Therefore, the tolerance, t
o
l, can be determined by considering speed and accuracy trade-off.


## 3.3. The Greeks

The Greeks are defined as changes in option value relative to changes in each independent variable. For example, Delta is the first derivative of the option value with respect to the underlying asset. Gamma is the second derivative of the value with respect to the asset. Theta is the first derivative of the option with respect to time. Vega is the option’s sensitivity to changes in the volatility. Rho is the first derivative of the option price with respect to interest rates. For more details about the Greeks, we refer the reader to [21].In general, the payoff structure of the option is very stiff because it has a discontinuity around the strike price. Therefore, there exists a difficulty to numerically measure an accurate value of the Greeks. However, as we use the proposed adaptive time-stepping strategy, we can control the accuracy of the Greeks.Next, we evaluate the option sensitivities of the cash-or-nothing option atx
=
100. In particular, because of the given grid points as shown in Figure 1, we use the linear interpolation to calculate the Greeks at x
=
100 in some numerical tests. For comparison, we also evaluate the Greeks by the exact formula using MATLAB code which is presented in the Appendix.


### 3.3.1. Delta

Delta (Δ) is defined as the rate of change of the option value with respect to small changes in the underlying asset price. By differentiating (12), we have the following exact formula: (13)
Δ
e
x
=
∂
u
∂
x
=
C
e
-
r
τ
N
′
d
σ
x
τ
.Also, we can define the Delta by using the finite difference discretizations at x
i. For example, (14)
Δ
i
n
≔
∂
u
∂
x
i
n
=
-
h
i
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
+
h
i
-
h
i
-
1
u
i
n
h
i
-
1
h
i
+
h
i
-
1
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 8(a) shows the Delta by closed-form formula (13), numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. To represent the accuracy of the Delta with respect to t
o
l, we investigate the error between closed-form solution and numerical approximation by our proposed method in Figure 8(b). Through these results, we can see that the adaptive time step strategy reduces the numerical error generated around the strike price K
=
100.Figure 8
(a) DeltaΔ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Δ
-
Δ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table3 shows the Delta and its error by the adaptive time step strategy with respect to t
o
l. As shown in Table 3, the Delta by our proposed method is more accurate as t
o
l is smaller.Table 3
DeltaΔ and its error Δ
-
Δ
e
x at x
=
100 and T
=
1
/
365. Here, Delta by the exact formula is Δ
e
x
=
25.40376.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Δ
24.7890
25.5980
25.5713
25.5553
25.5425
25.5481
Δ
-
Δ
e
x
-
0.61472
0.19425
0.16751
0.15149
0.13876
0.14432


### 3.3.2. Gamma

Gamma (Γ) is the rate of change of the Delta with respect to changes in the underlying asset price. In other words, Gamma is defined by the second partial derivative of the portfolio value with respect to underlying asset price. According to closed-form solution (12), we get the exact formula for Gamma as (15)
Γ
e
x
=
∂
2
u
∂
x
2
=
-
C
e
-
r
τ
d
+
σ
τ
N
′
d
σ
2
x
2
τ
.As a different approach, Gamma is described by applying numerical discretization at x
i as follows: (16)
Γ
i
n
≔
∂
2
u
∂
x
2
i
n
=
2
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
-
2
u
i
n
h
i
-
1
h
i
+
2
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 9 shows the Gamma at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. As shown in Figure 9, while the numerical solution by one time step Δ
τ
=
1
/
365 has a great deviation from the exact formula (15), numerical approximations by adaptive time step strategy and exact solution of Gamma are in good agreement. Also, from Figure 9(b), we observe that the accuracy of the Gamma can be controlled by t
o
l in our proposed method.Figure 9
(a) GammaΓ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Γ
-
Γ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table4 displays the numerical value by our proposed method and its corresponding errors of the Gamma at x
=
100. As tolerance is smaller, we can obtain a convergent value of Gamma.Table 4
GammaΓ and its error Γ
-
Γ
e
x at x
=
100 and T
=
1
/
365. Here, Γ
e
x
=
-
0.21170.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Γ
-
0.19796
-
0.20396
-
0.20378
-
0.20368
-
0.20359
-
0.20363
Γ
-
Γ
e
x
0.01374
0.00774
0.00792
0.00802
0.00811
0.00807


### 3.3.3. Theta

Theta (Θ) is the rate of change of the option value with respect to changes in the time to maturity, whose exact formula on cash-or-nothing option is given as (17)
Θ
e
x
=
∂
u
∂
t
=
-
∂
u
∂
τ
=
C
e
-
r
τ
r
N
d
+
N
′
d
d
+
σ
τ
2
τ
-
r
σ
τ
. Also, the numerical Theta is described as (18)
Θ
i
n
≔
∂
u
∂
τ
i
n
=
u
i
n
+
1
-
u
i
n
-
1
2
Δ
τ
. Here, we calculate the Theta by using the central difference approximation as (19)
Θ
i
=
u
x
i
,
T
+
Δ
τ
end
-
u
x
i
,
T
-
Δ
τ
end
2
Δ
τ
end
,where Δ
τ
end
=
0.01
T.In Figure10(a), we represent Θ by each one of the methods. As a result, numerical results by the adaptive time algorithm are superior to the other ones. Also, we can see that numerical results of Θ converge as t
o
l is smaller (see Figure 10(b) and Table 5).Table 5
ThetaΘ and its error Θ
-
Θ
e
x at x
=
100 and T
=
1
/
365. Here, Θ
e
x
=
20.54956.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Θ
19.41235
19.96587
19.88387
19.94125
19.93973
19.93678
Θ
-
Θ
e
x
-
1.13721
-
0.58369
-
0.66570
-
0.60832
-
0.60984
-
0.61279Figure 10
(a) ThetaΘ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Θ
-
Θ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


### 3.3.4. Vega

Vega (ν) is the rate of the option value with respect to small changes in the volatility of the underlying asset. By the closed-form solution of cash-or-nothing option, (20)
ν
e
x
=
∂
u
∂
σ
=
C
e
-
r
τ
d
σ
+
τ
N
′
d
. We can write the Vega by using the discretizations as (21)
ν
i
n
≔
∂
u
∂
σ
i
n
=
u
i
n
σ
+
Δ
σ
-
u
i
n
σ
-
Δ
σ
2
Δ
σ
,where Δ
σ
=
0.01.As shown in Figure11(a), we obtain a reasonable approximation for the Vega ν when we use the adaptive time step strategy. Moreover, the difference of ν and ν
e
x decays as smaller tolerance level t
o
l is determined.Figure 11
(a) Vegaν at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ν
-
ν
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)The results in Table6 show that our proposed method could be improved taking a lower tolerance parameter t
o
l.Table 6
Vegaν and its error ν
-
ν
e
x at x
=
100 and T
=
1
/
365. Here, ν
e
x
=
-
1.73998.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ν
-
1.72129
-
1.72138
-
1.72165
-
1.71663
-
1.72054
-
1.72260
ν
-
ν
e
x
0.01870
0.01860
0.01833
0.02335
0.01945
0.01739


### 3.3.5. Rho

Rho (ρ) is the rate of change of the option value with respect to small changes in the riskless interest rate r. The exact Rho formula on cash-or-nothing option is derived from (12) as (22)
ρ
e
x
=
∂
u
∂
r
=
C
e
-
r
τ
-
τ
N
d
+
τ
σ
N
′
d
. Rho is evaluated by numerical discretization as follows: (23)
ρ
i
n
≔
∂
u
∂
r
i
n
=
u
i
n
r
+
Δ
r
-
u
i
n
r
-
Δ
r
2
Δ
r
,where Δ
r
=
0.001.Similar to previous tests, we obtain good results for Rho (ρ) by the adaptive time step algorithm as shown in Figure 12 and Table 7. Particularly, when we compare results by our proposed method with one time step simulation, we see the superiority in terms of accuracy.Table 7
Rhoρ and its error ρ
-
ρ
e
x at x
=
100 and T
=
1
/
365. Here, ρ
e
x
=
6.82325.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ρ
7.02722
6.77944
6.78185
6.78329
6.78617
6.78449
ρ
-
ρ
e
x
0.20394
-
0.04384
-
0.04143
-
0.03999
-
0.03711
-
0.03879Figure 12
(a) Rhoρ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ρ
-
ρ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


## 3.3.1. Delta

Delta (Δ) is defined as the rate of change of the option value with respect to small changes in the underlying asset price. By differentiating (12), we have the following exact formula: (13)
Δ
e
x
=
∂
u
∂
x
=
C
e
-
r
τ
N
′
d
σ
x
τ
.Also, we can define the Delta by using the finite difference discretizations at x
i. For example, (14)
Δ
i
n
≔
∂
u
∂
x
i
n
=
-
h
i
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
+
h
i
-
h
i
-
1
u
i
n
h
i
-
1
h
i
+
h
i
-
1
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 8(a) shows the Delta by closed-form formula (13), numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. To represent the accuracy of the Delta with respect to t
o
l, we investigate the error between closed-form solution and numerical approximation by our proposed method in Figure 8(b). Through these results, we can see that the adaptive time step strategy reduces the numerical error generated around the strike price K
=
100.Figure 8
(a) DeltaΔ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Δ
-
Δ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table3 shows the Delta and its error by the adaptive time step strategy with respect to t
o
l. As shown in Table 3, the Delta by our proposed method is more accurate as t
o
l is smaller.Table 3
DeltaΔ and its error Δ
-
Δ
e
x at x
=
100 and T
=
1
/
365. Here, Delta by the exact formula is Δ
e
x
=
25.40376.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Δ
24.7890
25.5980
25.5713
25.5553
25.5425
25.5481
Δ
-
Δ
e
x
-
0.61472
0.19425
0.16751
0.15149
0.13876
0.14432


## 3.3.2. Gamma

Gamma (Γ) is the rate of change of the Delta with respect to changes in the underlying asset price. In other words, Gamma is defined by the second partial derivative of the portfolio value with respect to underlying asset price. According to closed-form solution (12), we get the exact formula for Gamma as (15)
Γ
e
x
=
∂
2
u
∂
x
2
=
-
C
e
-
r
τ
d
+
σ
τ
N
′
d
σ
2
x
2
τ
.As a different approach, Gamma is described by applying numerical discretization at x
i as follows: (16)
Γ
i
n
≔
∂
2
u
∂
x
2
i
n
=
2
u
i
-
1
n
h
i
-
1
h
i
-
1
+
h
i
-
2
u
i
n
h
i
-
1
h
i
+
2
u
i
+
1
n
h
i
h
i
-
1
+
h
i
.Figure 9 shows the Gamma at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. As shown in Figure 9, while the numerical solution by one time step Δ
τ
=
1
/
365 has a great deviation from the exact formula (15), numerical approximations by adaptive time step strategy and exact solution of Gamma are in good agreement. Also, from Figure 9(b), we observe that the accuracy of the Gamma can be controlled by t
o
l in our proposed method.Figure 9
(a) GammaΓ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Γ
-
Γ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)Table4 displays the numerical value by our proposed method and its corresponding errors of the Gamma at x
=
100. As tolerance is smaller, we can obtain a convergent value of Gamma.Table 4
GammaΓ and its error Γ
-
Γ
e
x at x
=
100 and T
=
1
/
365. Here, Γ
e
x
=
-
0.21170.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Γ
-
0.19796
-
0.20396
-
0.20378
-
0.20368
-
0.20359
-
0.20363
Γ
-
Γ
e
x
0.01374
0.00774
0.00792
0.00802
0.00811
0.00807


## 3.3.3. Theta

Theta (Θ) is the rate of change of the option value with respect to changes in the time to maturity, whose exact formula on cash-or-nothing option is given as (17)
Θ
e
x
=
∂
u
∂
t
=
-
∂
u
∂
τ
=
C
e
-
r
τ
r
N
d
+
N
′
d
d
+
σ
τ
2
τ
-
r
σ
τ
. Also, the numerical Theta is described as (18)
Θ
i
n
≔
∂
u
∂
τ
i
n
=
u
i
n
+
1
-
u
i
n
-
1
2
Δ
τ
. Here, we calculate the Theta by using the central difference approximation as (19)
Θ
i
=
u
x
i
,
T
+
Δ
τ
end
-
u
x
i
,
T
-
Δ
τ
end
2
Δ
τ
end
,where Δ
τ
end
=
0.01
T.In Figure10(a), we represent Θ by each one of the methods. As a result, numerical results by the adaptive time algorithm are superior to the other ones. Also, we can see that numerical results of Θ converge as t
o
l is smaller (see Figure 10(b) and Table 5).Table 5
ThetaΘ and its error Θ
-
Θ
e
x at x
=
100 and T
=
1
/
365. Here, Θ
e
x
=
20.54956.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
Θ
19.41235
19.96587
19.88387
19.94125
19.93973
19.93678
Θ
-
Θ
e
x
-
1.13721
-
0.58369
-
0.66570
-
0.60832
-
0.60984
-
0.61279Figure 10
(a) ThetaΘ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors Θ
-
Θ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


## 3.3.4. Vega

Vega (ν) is the rate of the option value with respect to small changes in the volatility of the underlying asset. By the closed-form solution of cash-or-nothing option, (20)
ν
e
x
=
∂
u
∂
σ
=
C
e
-
r
τ
d
σ
+
τ
N
′
d
. We can write the Vega by using the discretizations as (21)
ν
i
n
≔
∂
u
∂
σ
i
n
=
u
i
n
σ
+
Δ
σ
-
u
i
n
σ
-
Δ
σ
2
Δ
σ
,where Δ
σ
=
0.01.As shown in Figure11(a), we obtain a reasonable approximation for the Vega ν when we use the adaptive time step strategy. Moreover, the difference of ν and ν
e
x decays as smaller tolerance level t
o
l is determined.Figure 11
(a) Vegaν at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ν
-
ν
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)The results in Table6 show that our proposed method could be improved taking a lower tolerance parameter t
o
l.Table 6
Vegaν and its error ν
-
ν
e
x at x
=
100 and T
=
1
/
365. Here, ν
e
x
=
-
1.73998.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ν
-
1.72129
-
1.72138
-
1.72165
-
1.71663
-
1.72054
-
1.72260
ν
-
ν
e
x
0.01870
0.01860
0.01833
0.02335
0.01945
0.01739


## 3.3.5. Rho

Rho (ρ) is the rate of change of the option value with respect to small changes in the riskless interest rate r. The exact Rho formula on cash-or-nothing option is derived from (12) as (22)
ρ
e
x
=
∂
u
∂
r
=
C
e
-
r
τ
-
τ
N
d
+
τ
σ
N
′
d
. Rho is evaluated by numerical discretization as follows: (23)
ρ
i
n
≔
∂
u
∂
r
i
n
=
u
i
n
r
+
Δ
r
-
u
i
n
r
-
Δ
r
2
Δ
r
,where Δ
r
=
0.001.Similar to previous tests, we obtain good results for Rho (ρ) by the adaptive time step algorithm as shown in Figure 12 and Table 7. Particularly, when we compare results by our proposed method with one time step simulation, we see the superiority in terms of accuracy.Table 7
Rhoρ and its error ρ
-
ρ
e
x at x
=
100 and T
=
1
/
365. Here, ρ
e
x
=
6.82325.
tol
1.0
E
-
2
1.0
E
-
3
1.0
E
-
4
1.0
E
-
5
1.0
E
-
6
1.0
E
-
7
ρ
7.02722
6.77944
6.78185
6.78329
6.78617
6.78449
ρ
-
ρ
e
x
0.20394
-
0.04384
-
0.04143
-
0.03999
-
0.03711
-
0.03879Figure 12
(a) Rhoρ at time T
=
1
/
365 as a function of the asset price x by closed-form formula, numerical approximation for only one time step size Δ
τ
=
1
/
365, and adaptive time step strategy. (b) Errors ρ
-
ρ
e
x as a function of the asset price x by adaptive time-stepping strategy with respect to t
o
l.
(a)
(b)


## 3.4. Comparison with Adaptive and Uniform Time Step

In this section, we investigate the difference of numerical solutions with adaptive and uniform time step size. For comparison, we find the numerical results by uniform time step, which have similar RMSE to the results by the adaptive time step method. Table8 shows RMSE, N
τ, and CPU time with adaptive and uniform time step methods. As shown in Table 8, the uniform time step method is more efficient than the adaptive time step method at large tolerance (t
o
l
≥
1.0
E
-
5). However, we can see that when tolerance t
o
l is below 1.0
E
-
6, the results by adaptive time step are better than those by uniform time step. Furthermore, we obtain more accurate results by adaptive time method than by uniform method.Table 8
Comparison with adaptive and uniform time step methods.
tol
RMSE
N
τ
CPU time
Adaptive
1.0
E
-
2
0.32183
46
0.02790
Uniform
0.33442
9
0.00415
Adaptive
1.0
E
-
3
0.06770
213
0.11714
Uniform
0.06774
207
0.02237
Adaptive
1.0
E
-
4
0.06086
480
0.15395
Uniform
0.06086
476
0.04640
Adaptive
1.0
E
-
5
0.05676
2150
1.04035
Uniform
0.05676
2144
0.20096
Adaptive
1.0
E
-
6
0.05240
4801
1.41838
Uniform
0.05559
10000000
902.92370
Adaptive
1.0
E
-
7
0.05457
14214
4.24171
Uniform
0.05559
10000000
902.92370
Adaptive
1.0
E
-
8
0.05527
39776
12.89830
Uniform
0.05559
10000000
902.92370


## 4. Conclusion

In this paper, we presented the accurate numerical method for the Greeks close to the maturity time and in the neighborhood around the strike price. To reduce the difference between the numerical solution and the exact solution of the Greeks, we proposed an adaptive time step size strategy which is based on the truncation error. As a test problem, we considered a European cash-or-nothing call option. In order to show the effect of the adaptive time step, we had numerical tests for Delta, Gamma, Vega, Rho, and Theta on a variety of tolerances. Also, we compared the numerical results with the numerical solution by only one evolution having time stepΔ
τ
=
T. As a result, we obtained that the proposed method is accurate and practical in computing option price and its Greeks close to the maturity time.



---

*Source: 1586786-2016-04-04.xml*