# Reconstruction of the Time-Dependent Volatility Function Using the Black–Scholes Model

Authors: Yuzi Jin; Jian Wang; Sangkwon Kim; Youngjin Heo; Changwoo Yoo; Youngrock Kim; Junseok Kim; Darae Jeong
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

We propose a simple and robust numerical algorithm to estimate a time-dependent volatility function from a set of market observations, using the Black–Scholes (BS) model. We employ a fully implicit finite difference method to solve the BS equation numerically. To define the time-dependent volatility function, we define a cost function that is the sum of the squared errors between the market values and the theoretical values obtained by the BS model using the time-dependent volatility function. To minimize the cost function, we employ the steepest descent method. However, in general, volatility functions for minimizing the cost function are nonunique. To resolve this problem, we propose a predictor-corrector technique. As the first step, we construct the volatility function as a constant. Then, in the next step, our algorithm follows the prediction step and correction step at half-backward time level. The constructed volatility function is continuous and piecewise linear with respect to the time variable. We demonstrate the ability of the proposed algorithm to reconstruct time-dependent volatility functions using manufactured volatility functions. We also present some numerical results for real market data using the proposed volatility function reconstruction algorithm.


---

## Body

## 1. Introduction

The accurate calibration of models using market option data is one of the most important problems in finance [1]. The reason is for accurate pricing and accordingly hedging strategy [2]. According to [3], the authors proposed the accurate computations for Greeks using the numerical solutions of the Black–Scholes partial differential equation. The well-known standard Black–Scholes (BS) model is not adequate for calibrating the market option data because it uses the constant volatility [4–6]. As an alternative to the BS equation with constant volatility [7], local volatility models were introduced to explain the volatility smiles or skews observed in the market. There has been much research carried out regarding the reconstruction of local volatility functions from market data [8, 9]. For example, in [10], radial basis functions were used to construct local volatility surfaces. In [1], the authors proposed a regularized optimization formulation, using spline kernels to ensure both accuracy and stability in the local volatility function calibration. In [9], the authors mentioned the calibration of a local volatility surface for European options using a nonparametric approach by employing a second-order Tikhonov regularization.The main goal of this study is to propose a new simple and robust numerical method for the construction of a time-dependent volatility function, using the BS partial differential equation with nonconstant volatility [11, 12]:(1)∂uS,t∂t+12σtS2∂2uS,t∂S2+rS∂uS,t∂S=ruS,t,for (S,t)∈R+×[0,T), where u(S,t) is the option value of the underlying price S and time t, σ(t) is the volatility function of time t, and r is the riskless interest rate. The final condition is the payoff function u(S,T)=Λ(S) at expiry T. Typically, the local volatility function is given as a function of an underlying asset and time, that is, σ(S,t) [13, 14]. However, for simplicity and robustness of the solution algorithm, we assume that the local volatility function only depends on the time and is a piecewise linear function with respect to the time variable.The outline of this paper is as follows. In Section2, we describe our numerical algorithm for constructing the time-dependent volatility function. In Section 3, numerical experiments are presented. Finally, conclusions are drawn in Section 4.


## 2. Numerical Algorithm

In this section, we present a new numerical algorithm for constructing the time-varying volatility function.


### 2.1. Numerical Solution

LetS be the value of the underlying asset price and τ=T-t be the time to expiry, then (1) can be given as the following initial-value problem:(2)∂uS,τ∂τ=12στS2∂2uS,τ∂S2+rS∂uS,τ∂S-ruS,τ,for (S,τ)∈Ω×(0,T] with an initial condition u(S,0)=Λ(S)forS∈Ω=(0,L), where the infinite domain is truncated to a finite computational domain [15]. Now, to solve (2) numerically, we apply a finite difference method (FDM). Let us denote the numerical approximation of the solution of (2) by uin≡u(Si,τn)=u(ih,nΔτ), for i=1,2,…,NS and n=0,1,…,Nτ. Here, h=L/(NS-1) and Δτ=T/Nτ are uniform spatial and temporal step sizes, respectively. NS is the number of grid points and Nτ is the number of time steps. Figure 1 illustrates the uniform grid with a spatial step size h. Furthermore, the variable volatility σn is defined similarly as σn≡σ(τn).Figure 1
Uniform grid with a spatial step sizeh.By applying the fully implicit-in-time and space-centered difference scheme to (2), we obtain that (3)uin+1-uinΔτ=σn+1Si22ui-1n+1-2uin+1+ui+1n+1h2+rSiui+1n+1-ui-1n+12h-ruin+1.We can rewrite (3) as(4)αiui-1n+1+βiuin+1+γiui+1n+1=bi,fori=2,…,NS,where αi=rSi/2h-σn+1Si2/2h2, βi=1/Δτ+σn+1Si2/h2+r, γi=-rSi/2h-σn+1Si2/2h2, and bi=uin/Δτ. For boundary conditions, we use the zero Dirichlet condition at S1, that is, u1n=0, and the linear condition at SNS, that is, uNS+1n=2uNSn-uNS-1n, for all n [16]. To solve the resulting discrete system (4), we apply the Thomas algorithm [17].


### 2.2. Algorithm of Volatility Construction

Next, we describe the proposed algorithm for constructing the time-dependent volatility using option prices. Suppose that we have a set of measurements{ωβα}, where ωβα is the market price of the options with the exercise time Tα for α=1,…,Mt and the exercise price Kβ for β=1,…,Mk. Here, we assume that T1≤⋯≤TMt and K1≤⋯≤KMk. Using the given data, we determine a volatility function σ(t) in the least-squares sense. That is, we minimize the following mean-square error:(5)Γασ=1Mk∑β=1MkuKβσ;S0,Tα-ωβα2χβαforα=1,…,Mt,where uKβ(σ;S0,Tα) is the numerical solution at S=S0 of 2 with the strike price Kβ at time Tα. Here, χβα is a characteristic function, which is equal to one if ωβα is available and otherwise zero. In this paper, we use the steepest descent method [18] to find the optimal value σ that minimizes the cost function Γα(σ). Now, we present the detailed process of the proposed algorithm. We use the notation Tα+1/2=(Tα+Tα+1)/2 for α=0,1,…, where T0=0.Step 1 (find the constant volatility function on[0,T1]).
By assumingσ(t)=y1 on [0,T1], we find that y1 minimizes the following cost function:(6)Γ1σ=1Mk∑β=1MkuKβσ;S0,T1-ωβ12χβ1.
After determining the volatility function, we denoteσ1(t)≔σ(t). Figure 2(a) illustrates the constant volatility function at this step.Figure 2
Schematic of procedure for time-dependent volatility function.
(a) (b) (c) (d) (e) (f)Step 2 (find the linear volatility function on[0,T2]).
Withy1 obtained inStep1, we define the linear volatility function σ(t) on [0,T2] as(7)σt=y2-y1T2-0.5T1t-T2+y2.Here, (7) represents a linear function passing through the point (T1/2,y1) and having the optimal value y2 that minimizes the following cost function:(8)Γ2σ=1Mk∑β=1MkuKβσ;S0,T2-ωβ22χβ2.Note that, by setting (7), we only estimate a single parameter y2 using the given value for y1. After determining the volatility function, we denote σ2(t):=σ(t). Figure 2(b) illustrates the linear volatility function σ(t) on [0,T2].Step 3 (find the piecewise linear volatility function on[0,T3]).
First, we assume that the linear volatility functionσ3(t) on [T3/2,T3] is given by(9)σ3t=y3-σ2T3/2T3-T3/2t-T3+y3, where the linear equation passes through the given point (T3/2,σ2(T3/2)) and unknown point (T3,y3) as shown in Figure 2(c). From the linear function σ3, we define the following piecewise linear volatility function on [0,T3]:(10)σt=σ2t,ift∈0,T3/2,σ3t,ift∈T3/2,T3. Then, we estimate the optimal value of y3 that minimizes the cost function as(11)Γ3σ=1Mk∑β=1MkuKβσ;S0,T3-ωβ32χβ3.Figure 2(d) presents a schematic of the piecewise linear volatility function on [0,T3].Step 4 (find the linear volatility function on[0,TMt]).
The following process is repeated fromα=4 to α=Mt. By using yα-1 as obtained in previous step, we set the linear volatility function σα(t) on [Tα-3/2,Tα] as(12)σαt=yα-σα-1Tα-3/2Tα-Tα-3/2t-Tα+yα. In general, by using the volatility function(13)σt=σ2tift∈0,T3/2,σjtift∈Tj-3/2,Tj-1/2for2<j<α,σαtift∈Tα-3/2,Tα, we determine the optimal parameter yα that minimizes the cost function as (14)Γασ=1Mk∑β=1MkuKβσ;S0,Tα-ωβα2χβα.Note that in the proposed algorithm we do not compute the implied volatility and do not use Dupire’s formula [12, 19–21]. Instead, we directly compute the time-dependent volatility function from option prices and the BS equation. Our reconstructed volatility function is piecewise linear. It is very robust to compute the volatility because we only need to find a single parameter at each steepest descent algorithm.


## 2.1. Numerical Solution

LetS be the value of the underlying asset price and τ=T-t be the time to expiry, then (1) can be given as the following initial-value problem:(2)∂uS,τ∂τ=12στS2∂2uS,τ∂S2+rS∂uS,τ∂S-ruS,τ,for (S,τ)∈Ω×(0,T] with an initial condition u(S,0)=Λ(S)forS∈Ω=(0,L), where the infinite domain is truncated to a finite computational domain [15]. Now, to solve (2) numerically, we apply a finite difference method (FDM). Let us denote the numerical approximation of the solution of (2) by uin≡u(Si,τn)=u(ih,nΔτ), for i=1,2,…,NS and n=0,1,…,Nτ. Here, h=L/(NS-1) and Δτ=T/Nτ are uniform spatial and temporal step sizes, respectively. NS is the number of grid points and Nτ is the number of time steps. Figure 1 illustrates the uniform grid with a spatial step size h. Furthermore, the variable volatility σn is defined similarly as σn≡σ(τn).Figure 1
Uniform grid with a spatial step sizeh.By applying the fully implicit-in-time and space-centered difference scheme to (2), we obtain that (3)uin+1-uinΔτ=σn+1Si22ui-1n+1-2uin+1+ui+1n+1h2+rSiui+1n+1-ui-1n+12h-ruin+1.We can rewrite (3) as(4)αiui-1n+1+βiuin+1+γiui+1n+1=bi,fori=2,…,NS,where αi=rSi/2h-σn+1Si2/2h2, βi=1/Δτ+σn+1Si2/h2+r, γi=-rSi/2h-σn+1Si2/2h2, and bi=uin/Δτ. For boundary conditions, we use the zero Dirichlet condition at S1, that is, u1n=0, and the linear condition at SNS, that is, uNS+1n=2uNSn-uNS-1n, for all n [16]. To solve the resulting discrete system (4), we apply the Thomas algorithm [17].


## 2.2. Algorithm of Volatility Construction

Next, we describe the proposed algorithm for constructing the time-dependent volatility using option prices. Suppose that we have a set of measurements{ωβα}, where ωβα is the market price of the options with the exercise time Tα for α=1,…,Mt and the exercise price Kβ for β=1,…,Mk. Here, we assume that T1≤⋯≤TMt and K1≤⋯≤KMk. Using the given data, we determine a volatility function σ(t) in the least-squares sense. That is, we minimize the following mean-square error:(5)Γασ=1Mk∑β=1MkuKβσ;S0,Tα-ωβα2χβαforα=1,…,Mt,where uKβ(σ;S0,Tα) is the numerical solution at S=S0 of 2 with the strike price Kβ at time Tα. Here, χβα is a characteristic function, which is equal to one if ωβα is available and otherwise zero. In this paper, we use the steepest descent method [18] to find the optimal value σ that minimizes the cost function Γα(σ). Now, we present the detailed process of the proposed algorithm. We use the notation Tα+1/2=(Tα+Tα+1)/2 for α=0,1,…, where T0=0.Step 1 (find the constant volatility function on[0,T1]).
By assumingσ(t)=y1 on [0,T1], we find that y1 minimizes the following cost function:(6)Γ1σ=1Mk∑β=1MkuKβσ;S0,T1-ωβ12χβ1.
After determining the volatility function, we denoteσ1(t)≔σ(t). Figure 2(a) illustrates the constant volatility function at this step.Figure 2
Schematic of procedure for time-dependent volatility function.
(a) (b) (c) (d) (e) (f)Step 2 (find the linear volatility function on[0,T2]).
Withy1 obtained inStep1, we define the linear volatility function σ(t) on [0,T2] as(7)σt=y2-y1T2-0.5T1t-T2+y2.Here, (7) represents a linear function passing through the point (T1/2,y1) and having the optimal value y2 that minimizes the following cost function:(8)Γ2σ=1Mk∑β=1MkuKβσ;S0,T2-ωβ22χβ2.Note that, by setting (7), we only estimate a single parameter y2 using the given value for y1. After determining the volatility function, we denote σ2(t):=σ(t). Figure 2(b) illustrates the linear volatility function σ(t) on [0,T2].Step 3 (find the piecewise linear volatility function on[0,T3]).
First, we assume that the linear volatility functionσ3(t) on [T3/2,T3] is given by(9)σ3t=y3-σ2T3/2T3-T3/2t-T3+y3, where the linear equation passes through the given point (T3/2,σ2(T3/2)) and unknown point (T3,y3) as shown in Figure 2(c). From the linear function σ3, we define the following piecewise linear volatility function on [0,T3]:(10)σt=σ2t,ift∈0,T3/2,σ3t,ift∈T3/2,T3. Then, we estimate the optimal value of y3 that minimizes the cost function as(11)Γ3σ=1Mk∑β=1MkuKβσ;S0,T3-ωβ32χβ3.Figure 2(d) presents a schematic of the piecewise linear volatility function on [0,T3].Step 4 (find the linear volatility function on[0,TMt]).
The following process is repeated fromα=4 to α=Mt. By using yα-1 as obtained in previous step, we set the linear volatility function σα(t) on [Tα-3/2,Tα] as(12)σαt=yα-σα-1Tα-3/2Tα-Tα-3/2t-Tα+yα. In general, by using the volatility function(13)σt=σ2tift∈0,T3/2,σjtift∈Tj-3/2,Tj-1/2for2<j<α,σαtift∈Tα-3/2,Tα, we determine the optimal parameter yα that minimizes the cost function as (14)Γασ=1Mk∑β=1MkuKβσ;S0,Tα-ωβα2χβα.Note that in the proposed algorithm we do not compute the implied volatility and do not use Dupire’s formula [12, 19–21]. Instead, we directly compute the time-dependent volatility function from option prices and the BS equation. Our reconstructed volatility function is piecewise linear. It is very robust to compute the volatility because we only need to find a single parameter at each steepest descent algorithm.


## 3. Numerical Experiments

In this section, we demonstrate the performance of the proposed time-dependent volatility construction algorithm by numerical experiments with manufactured volatility functions and option data from a real market. All computations were performed on a 2.7 GHz Intel PC with 8 GB of RAM loaded with MATLAB 2016a [22].


### 3.1. Convergence Test

As the first numerical test, we perform a convergence test for the numerical scheme (4) with respect to h and Δτ. Let us consider a payoff function Λ(S)=max⁡(S-100,0) with r=0.015, S0=100, T=1, and L=400. We set a complicated volatility function, σ(t)=0.1cos⁡(4πt)-0.1t+0.2 (see Figure 3), to check the convergence of the employed numerical scheme.Figure 3
Volatility function,σ(t)=0.1cos⁡(4πt)-0.1t+0.2.Table1 shows the convergence of European call option prices at (S0,T) as we refine h and Δτ. From this point on, we will use h=1 and Δτ=1/360, because these values are sufficiently accurate for the present work, and we will use r=0.015 and L=400 unless otherwise specified.Table 1
Convergence test for European call option values.
Case h = 4 h = 2 h = 1 h = 0.5 Δτ=1/90 7.397038 7.432546 7.441332 7.443523 Δτ=1/180 7.389061 7.424523 7.433299 7.435488 Δτ=1/360 7.385044 7.420485 7.429255 7.431443 Δτ=1/720 7.383028 7.418458 7.427226 7.429413


### 3.2. Nonuniqueness of the Volatility Function

We consider a numerical test to demonstrate the nonuniqueness of a linear volatility function. Letσ(t)=-0.2t+0.3 be a given volatility function on a time interval [0,1]. Using this given volatility function, we generate European call prices by solving (4). The other parameters used are S0=100 and Kβ=92+2β for β=1,…,8. The payoff functions are set as Λβ(S)=max⁡(S-Kβ,0) for β=1,…,8. Table 2 presents the European call option prices that are generated using the given conditions.Table 2
Generated European call option prices for different strike pricesKβ.
β 1 2 3 4 5 6 7 8 Kβ 94 96 98 100 102 104 106 108 ωβ 12.22 11.07 10.00 9.00 8.08 7.23 6.45 5.74Now, we will determine the following linear volatility function which minimizes the cost function:(15)σt=σ1-σ0t+σ0,where σ0=σ(0) and σ1=σ(1) are constants to be computed. Figure 4(a) shows the cost function Γ(σ)=∑β=18uKβσ;100,1-ωβ2/8 on a parameter domain, 0.1≤σ0,σ1≤0.3. Figure 4(b) shows contour lines of the cost function Γ(σ). As shown in Figure 4(b), we obtain the different points (star markers) having the local minimum of the cost function Γ(σ) by applying the steepest descent method with different initial guesses (from A to G). In Figure 4(c), we represent these estimated linear volatility functions (15), satisfying Γ(σ)≈0. This means that the optimal parameter values for (15) are nonunique. Therefore, in the first step of the construction of the volatility function, it is reasonable to assume a constant volatility function.Figure 4
(a) Cost functionΓ(σ)=∑β=18uKβσ;100,1-ωβ2/8 on a parameter domain, 0.1≤σ0,σ1≤0.3. (b) Contour lines of the cost function and behavior of eight different points. (c) Estimated volatility functions for corresponding points.
(a) (b) (c)


### 3.3. Volatility Construction Using Manufactured Data

Next, we consider the following nonlinear manufactured volatility function:(16)σt=0.33t.First, we obtain reference values for the call option that is based on the given volatility function (16) by solving (4) with Tα=90α/360 for α=1,2,3,4 and Kβ=92+2β for β=1,2,…,8. The European call option prices generated by the volatility function (16) are represented in Table 3.Table 3
European call option prices generated by the volatility functionσ(t)=0.3/3t.
Kβ 94 96 98 100 102 104 106 108 T1=90/360 8.87 7.60 6.45 5.42 4.51 3.72 3.05 2.47 T2=180/360 10.27 9.06 7.95 6.94 6.03 5.20 4.47 3.82 T3=270/360 11.08 9.89 8.79 7.78 6.86 6.02 5.26 4.58 T4=360/360 11.61 10.43 9.33 8.32 7.39 6.54 5.76 5.06In Figure5, we see the procedure for constructing the volatility function up until t=T4. Here, we denote the temporal volatility function σ(t) at each time level t=Ti by solid lines and star markers. Also, we have inserted small figures to illustrate the cost function Γ(σ) against the volatility and highlighted the optimal point that minimizes Γ(σ) at t=Ti for i=1,2,3,4. From this figure, we can see that our proposed method satisfactorily recovers the exact volatility function.Figure 5
Estimated time-dependent volatility by our algorithm and the exact volatility function (16). Here, each small figure shows the cost function Γ(σ) as a function of σ and highlights the optimal value that minimizes Γ(σ) at time level t=Tα for α=1,2,3,4.


### 3.4. Basic Mechanism of the Proposed Algorithm

In this section, we introduce the basic mechanism of the proposed algorithm. For comparison, we consider a simple piecewise linear (SPL) volatility function:(17)σSPLt=σ1t=y1ift∈0,T1,σαt=yα-yα-1Tα-Tα-1t-Tα+yαift∈Tα-1,Tαforα=2,…,Mt.If we apply the SPL volatility function (17) to approximate the two exact volatility functions σ(t)=0.3/3t and σ(t)=0.1cos⁡(4πt)-0.1t+0.2, then we obtain the results (dashed lines) illustrated in Figure 6. In Figure 6(a), the volatility is constant in the first interval, and in the second interval the value of y2 should be smaller than the exact value σ(T2) to minimize the cost function. Likewise, y3 is larger than σ(T3) and y4 is again smaller than σ(T4). Figure 6(b) shows similar behavior; that is, an oscillatory solution is generated. However, the proposed algorithm for constructing the volatility function is essentially a predictor-corrector technique. It dampens spurious oscillations and generates satisfactory results. As shown in Figure 6, the results (thin solid lines) using the proposed algorithm are accurate and robust.Figure 6
Construction of volatility by the SPL and proposed algorithms when the exact volatility is (a)σ(t)=0.3/3t and (b) σ(t)=0.1cos⁡(4πt)-0.1t+0.2.
(a) (b)


### 3.5. Calculate for the Whole Time at Once

Until now, we have found an optimized volatility for each section of time. But from now on, we will not calculate the time for each section but the volatility for all of the time. The comparison between the two results is shown in Figure7. Table 4 shows the value of Γ(σ).Table 4
The value ofΓ(σ).
Proposed volatility New proposed volatility Γ ( σ ) 1.4796 e - 04 6.1036 e - 06Figure 7
Construction of volatility by proposed algorithm and new proposed algorithm when the exact volatility isσ(t)=0.3/3t.


### 3.6. Volatility Construction from KOSPI 200 Data

Now, we construct the volatility function by the proposed algorithm with the KOSPI 200 index call option data [23, 24] on 29 July 2016. Table 5 shows real market data for the KOSPI 200 index call option price with respect to the strike and maturity. Here, the strike indices are Kβ=245+2.5(β-1) for β=1,2,…,8 and the maturity times are Δτ=1/365, T1=13Δτ, T2=41Δτ, and T3=76Δτ. At this time, the current value of the KOSPI 200 index was S0=251.48 and interest rate was r=0.0136.Table 5
KOSPI 200 index call option price on 29 July 2016 with respect to the strike and maturity.
Kβ 245.0 247.5 250.0 252.5 255.0 257.5 260.0 262.5 T1=13Δτ 7.34 5.29 3.39 1.93 0.92 0.37 0.11 0.04 T2=41Δτ 8.85 6.90 5.28 3.82 2.56 1.70 1.05 0.62 T3=76Δτ 9.34 8.41 6.95 4.72 4.24 3.21 2.36 1.74By the proposed algorithm, we construct the time-dependent volatility functionσ(t) from the KOSPI 200 market data as shown in Figure 8(d). Figures 8(a)–8(c) represent the numerical option values given by the proposed volatility function and the real market data for the KOSPI 200 index at t=T1, T2, and T3.Figure 8
Comparison of KOSPI 200 index call option price (circle-marked solid line) with numerical option prices (dashed line) at (a)T=13/365, (b) 41/365, and (c) 76/365. (d) Constructed time-dependent volatility function.
(a) (b) (c) (d)


## 3.1. Convergence Test

As the first numerical test, we perform a convergence test for the numerical scheme (4) with respect to h and Δτ. Let us consider a payoff function Λ(S)=max⁡(S-100,0) with r=0.015, S0=100, T=1, and L=400. We set a complicated volatility function, σ(t)=0.1cos⁡(4πt)-0.1t+0.2 (see Figure 3), to check the convergence of the employed numerical scheme.Figure 3
Volatility function,σ(t)=0.1cos⁡(4πt)-0.1t+0.2.Table1 shows the convergence of European call option prices at (S0,T) as we refine h and Δτ. From this point on, we will use h=1 and Δτ=1/360, because these values are sufficiently accurate for the present work, and we will use r=0.015 and L=400 unless otherwise specified.Table 1
Convergence test for European call option values.
Case h = 4 h = 2 h = 1 h = 0.5 Δτ=1/90 7.397038 7.432546 7.441332 7.443523 Δτ=1/180 7.389061 7.424523 7.433299 7.435488 Δτ=1/360 7.385044 7.420485 7.429255 7.431443 Δτ=1/720 7.383028 7.418458 7.427226 7.429413


## 3.2. Nonuniqueness of the Volatility Function

We consider a numerical test to demonstrate the nonuniqueness of a linear volatility function. Letσ(t)=-0.2t+0.3 be a given volatility function on a time interval [0,1]. Using this given volatility function, we generate European call prices by solving (4). The other parameters used are S0=100 and Kβ=92+2β for β=1,…,8. The payoff functions are set as Λβ(S)=max⁡(S-Kβ,0) for β=1,…,8. Table 2 presents the European call option prices that are generated using the given conditions.Table 2
Generated European call option prices for different strike pricesKβ.
β 1 2 3 4 5 6 7 8 Kβ 94 96 98 100 102 104 106 108 ωβ 12.22 11.07 10.00 9.00 8.08 7.23 6.45 5.74Now, we will determine the following linear volatility function which minimizes the cost function:(15)σt=σ1-σ0t+σ0,where σ0=σ(0) and σ1=σ(1) are constants to be computed. Figure 4(a) shows the cost function Γ(σ)=∑β=18uKβσ;100,1-ωβ2/8 on a parameter domain, 0.1≤σ0,σ1≤0.3. Figure 4(b) shows contour lines of the cost function Γ(σ). As shown in Figure 4(b), we obtain the different points (star markers) having the local minimum of the cost function Γ(σ) by applying the steepest descent method with different initial guesses (from A to G). In Figure 4(c), we represent these estimated linear volatility functions (15), satisfying Γ(σ)≈0. This means that the optimal parameter values for (15) are nonunique. Therefore, in the first step of the construction of the volatility function, it is reasonable to assume a constant volatility function.Figure 4
(a) Cost functionΓ(σ)=∑β=18uKβσ;100,1-ωβ2/8 on a parameter domain, 0.1≤σ0,σ1≤0.3. (b) Contour lines of the cost function and behavior of eight different points. (c) Estimated volatility functions for corresponding points.
(a) (b) (c)


## 3.3. Volatility Construction Using Manufactured Data

Next, we consider the following nonlinear manufactured volatility function:(16)σt=0.33t.First, we obtain reference values for the call option that is based on the given volatility function (16) by solving (4) with Tα=90α/360 for α=1,2,3,4 and Kβ=92+2β for β=1,2,…,8. The European call option prices generated by the volatility function (16) are represented in Table 3.Table 3
European call option prices generated by the volatility functionσ(t)=0.3/3t.
Kβ 94 96 98 100 102 104 106 108 T1=90/360 8.87 7.60 6.45 5.42 4.51 3.72 3.05 2.47 T2=180/360 10.27 9.06 7.95 6.94 6.03 5.20 4.47 3.82 T3=270/360 11.08 9.89 8.79 7.78 6.86 6.02 5.26 4.58 T4=360/360 11.61 10.43 9.33 8.32 7.39 6.54 5.76 5.06In Figure5, we see the procedure for constructing the volatility function up until t=T4. Here, we denote the temporal volatility function σ(t) at each time level t=Ti by solid lines and star markers. Also, we have inserted small figures to illustrate the cost function Γ(σ) against the volatility and highlighted the optimal point that minimizes Γ(σ) at t=Ti for i=1,2,3,4. From this figure, we can see that our proposed method satisfactorily recovers the exact volatility function.Figure 5
Estimated time-dependent volatility by our algorithm and the exact volatility function (16). Here, each small figure shows the cost function Γ(σ) as a function of σ and highlights the optimal value that minimizes Γ(σ) at time level t=Tα for α=1,2,3,4.


## 3.4. Basic Mechanism of the Proposed Algorithm

In this section, we introduce the basic mechanism of the proposed algorithm. For comparison, we consider a simple piecewise linear (SPL) volatility function:(17)σSPLt=σ1t=y1ift∈0,T1,σαt=yα-yα-1Tα-Tα-1t-Tα+yαift∈Tα-1,Tαforα=2,…,Mt.If we apply the SPL volatility function (17) to approximate the two exact volatility functions σ(t)=0.3/3t and σ(t)=0.1cos⁡(4πt)-0.1t+0.2, then we obtain the results (dashed lines) illustrated in Figure 6. In Figure 6(a), the volatility is constant in the first interval, and in the second interval the value of y2 should be smaller than the exact value σ(T2) to minimize the cost function. Likewise, y3 is larger than σ(T3) and y4 is again smaller than σ(T4). Figure 6(b) shows similar behavior; that is, an oscillatory solution is generated. However, the proposed algorithm for constructing the volatility function is essentially a predictor-corrector technique. It dampens spurious oscillations and generates satisfactory results. As shown in Figure 6, the results (thin solid lines) using the proposed algorithm are accurate and robust.Figure 6
Construction of volatility by the SPL and proposed algorithms when the exact volatility is (a)σ(t)=0.3/3t and (b) σ(t)=0.1cos⁡(4πt)-0.1t+0.2.
(a) (b)


## 3.5. Calculate for the Whole Time at Once

Until now, we have found an optimized volatility for each section of time. But from now on, we will not calculate the time for each section but the volatility for all of the time. The comparison between the two results is shown in Figure7. Table 4 shows the value of Γ(σ).Table 4
The value ofΓ(σ).
Proposed volatility New proposed volatility Γ ( σ ) 1.4796 e - 04 6.1036 e - 06Figure 7
Construction of volatility by proposed algorithm and new proposed algorithm when the exact volatility isσ(t)=0.3/3t.


## 3.6. Volatility Construction from KOSPI 200 Data

Now, we construct the volatility function by the proposed algorithm with the KOSPI 200 index call option data [23, 24] on 29 July 2016. Table 5 shows real market data for the KOSPI 200 index call option price with respect to the strike and maturity. Here, the strike indices are Kβ=245+2.5(β-1) for β=1,2,…,8 and the maturity times are Δτ=1/365, T1=13Δτ, T2=41Δτ, and T3=76Δτ. At this time, the current value of the KOSPI 200 index was S0=251.48 and interest rate was r=0.0136.Table 5
KOSPI 200 index call option price on 29 July 2016 with respect to the strike and maturity.
Kβ 245.0 247.5 250.0 252.5 255.0 257.5 260.0 262.5 T1=13Δτ 7.34 5.29 3.39 1.93 0.92 0.37 0.11 0.04 T2=41Δτ 8.85 6.90 5.28 3.82 2.56 1.70 1.05 0.62 T3=76Δτ 9.34 8.41 6.95 4.72 4.24 3.21 2.36 1.74By the proposed algorithm, we construct the time-dependent volatility functionσ(t) from the KOSPI 200 market data as shown in Figure 8(d). Figures 8(a)–8(c) represent the numerical option values given by the proposed volatility function and the real market data for the KOSPI 200 index at t=T1, T2, and T3.Figure 8
Comparison of KOSPI 200 index call option price (circle-marked solid line) with numerical option prices (dashed line) at (a)T=13/365, (b) 41/365, and (c) 76/365. (d) Constructed time-dependent volatility function.
(a) (b) (c) (d)


## 4. Conclusions

In this paper, we considered a simple and robust numerical algorithm for the construction of a time-dependent volatility function from a finite set of market observations by the Black–Scholes model. In order to numerically solve the BS equation, we applied a fully implicit finite difference method. By assuming a continuous and piecewise linear function with respect to time, we constructed a volatility function at each time level. The construction algorithm for the volatility function is based on the minimization of a cost function, defined by the sum of the squared errors between market values and theoretical values obtained from the BS model using the time-dependent volatility function. Here, we used a steepest descent method to minimize the cost function. However, in general, volatility functions for minimizing the cost function are nonunique. To resolve this problem, we proposed a predictor-corrector technique. As the first step, we construct the volatility function as a constant. Then, in the next step, our algorithm follows the prediction step and correction step at half-backward time level. Using several examples, we demonstrated the accuracy and robustness of our proposed algorithm.



---

*Source: 3093708-2018-05-07.xml*