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.

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 and denote the price of the underlying asset and time, respectively. Let be the set of underlying assets. Then, the value of an option is governed by the following -dimensional BS equation [2].

For , with a final condition , where is the constant riskless interest rate, are volatility values of , and are the asset correlations between and . is the payoff function at maturity . There are three classical techniques which are the finite difference method (FDM) [313], 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 Section 2, 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 which means the time to expiry, the one-dimensional version of (1) becomes the following initial value problem: with an initial condition for Here, is a sufficiently large asset price. The BS equation (2) is discretized on a grid defined by and for , where is the number of grid intervals and is the grid spacing (see Figure 1). We assume that and the ghost point .

Let be the numerical approximate solution, where is the time step size and is the total number of time steps. By applying implicit Euler’s scheme to (2), we have for and . Then, we can rewrite (3) as In this paper, we restrict our attention to European call options. Since the option price at approaches zero, we impose the zero Dirichlet boundary condition as . Also, the option value at sufficiently large asset price is asymptotically linear. Therefore, we use the linear boundary condition [1618] at as ; that is, . 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 by and the two approximate solutions by (one step with ) and (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 , And, for the numerical solutions with two steps by , By adding (6) and (7), we obtainLet be the difference between the two numerical estimates; that is, where means the constant value whose order of magnitude is . As shown in Figure 2, there is a difference of two numerical solutions by one step with and two steps with . In general, the numerical approximation by two steps with is more accurate than by one step 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, and . Next, let be an initial time step size. With a given numerical solution and a time step , we solve (4) twice to get . Next, with the given numerical solution and a twice larger time step , we solve (4) to get . We define the time step scaled error as where and are the numerical solution with times and , respectively. If the error is below the given tolerance, then we set the th numerical solution as . Otherwise, we solve (4) using each time step 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 and the error as . If , that is, , the new time step size is larger than the old one. Otherwise, which is , the new one is smaller than the old one.

The adaptive time-stepping strategy can be summarized in Algorithm 1.

Require: Set the initial conditions and , the expiry time , the maximum and minimum
time steps and , tolerance tol, truncation error , safety
factor , time , iteration number .
 While do
  
  if then
   Set .
  else
   Set .
  end if
  if then
   Set .
   Solve (4) with time step and set its approximation as .
  end if
  Solve (4) twice with time step and set its approximation as .
  Solve (4) with time step and set its approximation as .
  Calculate where and is the total number of .
  while do
   Set .
   if then
    Set .
   end if
   Solve (4) twice with and set its approximation as .
   Solve (4) with and set its approximation as .
   Calculate
  end while
  Set .
  Set and .
end while

3. Numerical Experiments

In this section, for numerical experiments, we consider a European cash-or-nothing option which pays an amount at maturity if option is in-the-money state. The payoff function is given by where is the strike price and denotes the return value. The closed-form solution [21] for the option is given as where and is the cumulative distribution function for the standard normal distribution [1]. In the Appendix, MATLAB code for the closed-form solution is presented.

Figure 3 shows the cash-or-nothing option prices at , , and on . Here, the option prices at and are obtained by (12). And we use strike price , cash , the risk-free interest rate , and volatility . 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.

In the following sections, unless otherwise specified, we use strike price , cash , the risk-free interest rate , the volatility , and . All computations are performed using MATLAB version [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 error , where and 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 that is the position near the predetermined strike price . Figure 4 shows the spatial grid structure near with respect to . And 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, . Table 1 represents the absolute error between the numerical and exact solutions with respect to spatial and temporal step sizes.

As shown in Table 1, the absolute error decreases as and decrease. Also, we can see these results in Figure 5.

In the results, we confirm that the absolute error is more influenced by temporal step size than by spatial step size . In particular, the numerical results with 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 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 along with the truncation error of the numerical solution as we described before.

As the default parameters, we set the initial time step size , where and . Here, and represent one day and one second, respectively. Also, we use the spatial step size according to results in Table 1.

Table 2 represents the root mean square error (RMSE) of the numerical solution which is measured in the area . Here, RMSE is calculated by , where and denote the numerical and exact solutions. In Table 2, and represent the minimum and maximum time step sizes which are used during the numerical iteration. And denotes the total number of iterations during one day, . Note that the RMSE by only one time step size is . For comparison, we evaluate the ratio of to the RMSE by adaptive time-stepping strategy. As shown in Table 2, the numerical solution by the adaptive time-stepping method is about times more accurate than that by the one time step.

Figure 6 shows the time step size against time over one day, , in adaptive time-stepping strategy with . 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 .

Next, we perform numerical tests around maturity and compare the numerical solution with the exact solution. Let us define and 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 near maturity makes inaccurate solution. Figure 7(b) illustrates with respect to the changes of .

We note that as the decreases, the corresponding numerical solution remains accuracy, however, for a longer time. Therefore, the tolerance, , 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 at . In particular, because of the given grid points as shown in Figure 1, we use the linear interpolation to calculate the Greeks at 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: Also, we can define the Delta by using the finite difference discretizations at . For example, Figure 8(a) shows the Delta by closed-form formula (13), numerical approximation for only one time step size , and adaptive time step strategy. To represent the accuracy of the Delta with respect to , 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 .

Table 3 shows the Delta and its error by the adaptive time step strategy with respect to . As shown in Table 3, the Delta by our proposed method is more accurate as is smaller.

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 As a different approach, Gamma is described by applying numerical discretization at as follows: Figure 9 shows the Gamma at time as a function of the asset price by closed-form formula, numerical approximation for only one time step size , and adaptive time step strategy. As shown in Figure 9, while the numerical solution by one time step 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 in our proposed method.

Table 4 displays the numerical value by our proposed method and its corresponding errors of the Gamma at . As tolerance is smaller, we can obtain a convergent value of Gamma.

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 Also, the numerical Theta is described as Here, we calculate the Theta by using the central difference approximation as where .

In Figure 10(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 is smaller (see Figure 10(b) and Table 5).

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, We can write the Vega by using the discretizations as where .

As shown in Figure 11(a), we obtain a reasonable approximation for the Vega when we use the adaptive time step strategy. Moreover, the difference of and decays as smaller tolerance level is determined.

The results in Table 6 show that our proposed method could be improved taking a lower tolerance parameter .

3.3.5. Rho

Rho () is the rate of change of the option value with respect to small changes in the riskless interest rate . The exact Rho formula on cash-or-nothing option is derived from (12) as Rho is evaluated by numerical discretization as follows: where .

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.

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. Table 8 shows RMSE, , 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 (). However, we can see that when tolerance is below , 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.

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 . 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.

Appendix

MATLAB Code

The MATLAB [22] script is composed of cash-or-nothing option value and its Greeks by the closed-form solutions of Black-Scholes equation (see Code 1).

clear all; close all; clc; clf;
K = 100; L = 2.5*K; T = 1/365; sigma = 0.3; r = 0.03;
cash = 100; Nx = L+1; x = linspace(0,L,Nx); h = x()-x();
payoff(1:Nx) = 0.0; payoff(x >= K) = cash;
%% Price %%
d1 = (log(x/K) + (r+sigma2/2)*T)/(sigma*sqrt(T));
d2 = d1 - (sigma*sqrt(T));
exact = cash*exp(-r*T).*normcdf(d1);
figure(); hold on; grid on;
plot(x,payoff,,k-,,x,exact,,r-,); axis([0 L -5 1.05*cash]);
%% Delta %%
del = cash*exp(-r*T).*normpdf(d1)./(sigma*sqrt(T)*x);
figure(); hold on; grid on;
plot(del,,r-,,,LineWidth,,1); axis([0 L -5 1.1*max(del)]);
%% Gamma %%
gam = -cash*exp(-r*T) * (d1.*normpdf(d2))./((sigma*x).2*T);
figure(); hold on; grid on; plot(gam,,r-,,,LineWidth,,1);
axis([0 L -1.1*max(gam) 1.1*max(gam)]);
%% theta %%
the = cash*exp(-r*T)*(r*normcdf(d2) + normpdf(d2)...
.* (d1/(2*T) - r/(sigma*sqrt(T))));
figure(); hold on; grid on; plot(the,,r-,,,LineWidth,,1);
axis([0 L -1.1*max(the) 1.1*max(the)]);
%% Vega %%
veg = -cash*exp(-r*T)*(normpdf(d2)).*d1/sigma;
figure(); hold on; grid on; plot(veg,,r-,,,LineWidth,,1);
axis([0 L -1.1*max(veg) 1.1*max(veg)]);
%% Rho %%
rho = cash*exp(-r*T)*((-T*normcdf(d2)) + sqrt(T)/sigma*(normpdf(d2)));
figure(); hold on; grid on; plot(rho,,r-,,,LineWidth,,1);
axis([0 L 3*min(rho) 1.1*max(rho)]);

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

The first author (Darae Jeong) was supported by a Korea University grant. The corresponding author (Junseok Kim) was supported by a subproject of the project Research for Applications of Mathematical Principles (no. C21501) and supported by the National Institute for Mathematical Sciences (NIMS).