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.
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]:for , where is the option value of the underlying price and time , is the volatility function of time , and is the riskless interest rate. The final condition is the payoff function at expiry . Typically, the local volatility function is given as a function of an underlying asset and time, that is, [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 Section 2, 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
Let be the value of the underlying asset price and be the time to expiry, then (1) can be given as the following initial-value problem:for with an initial condition , 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 , for and . Here, and are uniform spatial and temporal step sizes, respectively. is the number of grid points and is the number of time steps. Figure 1 illustrates the uniform grid with a spatial step size . Furthermore, the variable volatility is defined similarly as .
By applying the fully implicit-in-time and space-centered difference scheme to (2), we obtain that We can rewrite (3) aswhere , , , and . For boundary conditions, we use the zero Dirichlet condition at , that is, , and the linear condition at , that is, , for all [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 for and the exercise price for . Here, we assume that and . Using the given data, we determine a volatility function in the least-squares sense. That is, we minimize the following mean-square error:where is the numerical solution at of 2 with the strike price at time . 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 for , where .
Step 1 (find the constant volatility function on ). By assuming on , we find that minimizes the following cost function:After determining the volatility function, we denote . Figure 2(a) illustrates the constant volatility function at this step.
(a)
(b)
(c)
(d)
(e)
(f)
Step 2 (find the linear volatility function on ). With obtained in Step 1, we define the linear volatility function on asHere, (7) represents a linear function passing through the point and having the optimal value that minimizes the following cost function:Note that, by setting (7), we only estimate a single parameter using the given value for . After determining the volatility function, we denote . Figure 2(b) illustrates the linear volatility function on .
Step 3 (find the piecewise linear volatility function on ). First, we assume that the linear volatility function on is given by where the linear equation passes through the given point and unknown point as shown in Figure 2(c). From the linear function , we define the following piecewise linear volatility function on : Then, we estimate the optimal value of that minimizes the cost function asFigure 2(d) presents a schematic of the piecewise linear volatility function on .
Step 4 (find the linear volatility function on ). The following process is repeated from to . By using as obtained in previous step, we set the linear volatility function on as In general, by using the volatility function we determine the optimal parameter that minimizes the cost function as
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 and . Let us consider a payoff function with , , , and . We set a complicated volatility function, (see Figure 3), to check the convergence of the employed numerical scheme.
Table 1 shows the convergence of European call option prices at as we refine and . From this point on, we will use and , because these values are sufficiently accurate for the present work, and we will use and unless otherwise specified.
3.2. Nonuniqueness of the Volatility Function
We consider a numerical test to demonstrate the nonuniqueness of a linear volatility function. Let be a given volatility function on a time interval . Using this given volatility function, we generate European call prices by solving (4). The other parameters used are and for . The payoff functions are set as for . Table 2 presents the European call option prices that are generated using the given conditions.
Now, we will determine the following linear volatility function which minimizes the cost function:where and are constants to be computed. Figure 4(a) shows the cost function on a parameter domain, . 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 . 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.
(a)
(b)
(c)
3.3. Volatility Construction Using Manufactured Data
Next, we consider the following nonlinear manufactured volatility function:First, we obtain reference values for the call option that is based on the given volatility function (16) by solving (4) with for and for . The European call option prices generated by the volatility function (16) are represented in Table 3.
In Figure 5, we see the procedure for constructing the volatility function up until . Here, we denote the temporal volatility function at each time level 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 for . From this figure, we can see that our proposed method satisfactorily recovers the exact volatility function.
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:If we apply the SPL volatility function (17) to approximate the two exact volatility functions and , 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 should be smaller than the exact value to minimize the cost function. Likewise, is larger than and is again smaller than . 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.
(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 Figure 7. Table 4 shows the value of .
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 for and the maturity times are , , , and . At this time, the current value of the KOSPI 200 index was and interest rate was .
By the proposed algorithm, we construct the time-dependent volatility function 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 , , and .
(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.
Disclosure
An earlier version of the manuscript was presented as a poster at SIAM 2017, Spring Conference, Seoul National University, Seoul, Korea, June 23-24, 2017.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author (Yuzi Jin) was supported by the University-Level Major Social Science Project of Jilin Institute of Chemical Technology, no. 20, 2016. The author (Youngrock Kim) was supported by the Basic Science Research Program of the NRF (Korea) under Grant no. 2015R1D1A1A01059643/2.