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

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 timestepping algorithm based on local truncation error. As a test problem for our numerical method, we consider a European cashor-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.


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 s = ( 1 ,  2 , . . .,   ) be the set of  underlying assets.Then, the value of an option (s, ) is governed by the following -dimensional BS equation [ with a final condition (s, ) = Λ(s), where  is the constant riskless interest rate,   are volatility values of   , and   are the asset correlations between   and   .Λ(s) is the payoff function at maturity .There are three classical techniques which are the finite difference method (FDM) [3][4][5][6][7][8][9][10][11][12][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 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.

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 (, 0) for  ∈ Ω = (0,  max ).Here,  max is a sufficiently large asset price.The BS equation ( 2) is discretized on a grid defined by  0 = 0 and  +1 =   + ℎ  for  = 0, . . .,   − 1, where   is the number of grid intervals and ℎ  is the grid spacing (see Figure 1).We assume that    =  max and the ghost point    +1 =  max + ℎ   −1 .
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  = 1, . . .,   and  = 0, . . .,   − 1.Then, we can rewrite (3) as In this paper, we restrict our attention to European call options.Since the option price at  = 0 approaches zero, we impose the zero Dirichlet boundary condition as (0, ) = 0. Also, the option value at sufficiently large asset price is asymptotically linear.Therefore, we use the linear boundary condition [16][17][18] at  =  max as The linear system from ( 4) is solved by the Thomas algorithm [19].

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 (,  + 2Δ) and the two approximate solutions by V (one step with 2Δ) 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 2Δ,  (,  + 2Δ) =  (, ) + 2Δ  (, ) Let  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 2Δ and two steps with Δ.In general, the numerical approximation by two steps with Δ is more accurate than by one step with 2Δ.
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, Δ max and Δ min .Next, let Δ 0 be an initial time step size.With a given numerical solution    a time step Δ = Δ 0 , we solve (4) twice to get  +2 .Next, with the given numerical solution   and a twice larger time step Δ = 2Δ 0 , we solve (4) to get V +2 .We define the time step scaled error as where  +2 and V +2 are the numerical solution with times Δ 0 and 2Δ 0 , respectively.If the error is below the given tolerance, then we set the (+2)th numerical solution as  +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 tol and the error  tr as Δ new = Δ 0 × tol/ tr .If  tr < tol, that is, tol/ tr > 1, the new time step size is larger than the old one.Otherwise, which is  tr > tol, the new one is smaller than the old one.The adaptive time-stepping strategy can be summarized in Algorithm 1.

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  = [ln(/) + ( − 0.5 2 )]/(√) and )d 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  = 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  = 100, cash  = 100, the risk-free interest rate  = 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.
In the following sections, unless otherwise specified, we use strike price  = 100, cash  = 100, the risk-free interest rate  = 0.03, the volatility  = 0.3, and  = 250.All computations are performed using MATLAB version 8 [22].

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 e = | exact − |, where  exact 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  = 100.5 that is the position near the predetermined strike price  = 100.Figure 4 shows the spatial grid structure near  with respect to ℎ.And the marked circle ( = 100.5)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.
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 ℎ ≤ 1/3 are Require: Set the initial conditions  0 and V 0 , the expiry time , the maximum and minimum time steps Δ max =  and Δ min = /86400, tolerance tol, truncation error  tr , safety factor  = 0.8, time  = 0, iteration number  = 0.
While  <  do 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 ℎ = 1/3 when we use the adaptive time-stepping strategy.

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.
Table 2    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.
Figure 6 shows the time step size Δ against time  over one day,  = 1/365, in adaptive time-stepping strategy with tol = 1.0 − 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  = 1/365.Next, we perform numerical tests around maturity and compare the numerical solution with the exact solution.Let us define  ex 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 Δ = 1/365 near We note that as the tol decreases, the corresponding numerical solution remains accuracy, however, for a longer time.Therefore, the tolerance, tol, can be determined by considering speed and accuracy trade-off.

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-ornothing option at  = 100.In particular, because of the given grid points as shown in Figure 1, we use the linear interpolation to calculate the Greeks at  = 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.asset price.By differentiating (12), we have the following exact formula:

Delta. Delta (Δ) is defined as the rate of change of the option value with respect to small changes in the underlying
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 Δ = 1/365, and adaptive time step strategy.To represent the accuracy of the Delta with respect to tol, 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  = 100.Table 3 shows the Delta and its error by the adaptive time step strategy with respect to tol.As shown in Table 3, the Delta by our proposed method is more accurate as tol is smaller.

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  = 1/365 as a function of the asset price  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 tol in our proposed method.Table 4 displays the numerical value by our proposed method and its corresponding errors of the Gamma at  = 100.As tolerance is smaller, we can obtain a convergent value of Gamma.

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 Δ end = 0.01.
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 tol is smaller (see Figure 10(b) and Table 5).

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 Δ = 0.01.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 ] ex decays as smaller tolerance level tol is determined.
The results in Table 6 show that our proposed method could be improved taking a lower tolerance parameter tol.
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 Δ = 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.

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.

Figure 5 :
Figure 5: Absolute errors between numerical and exact solutions on  = 100.5 and  = 1/365 with respect to spatial and temporal step sizes.

Figure 8 :
Figure 8: (a) Delta Δ at time  = 1/365 as a function of the asset price  by closed-form formula, numerical approximation for only one time step size Δ = 1/365, and adaptive time step strategy.(b) Errors Δ − Δ ex as a function of the asset price  by adaptive time-stepping strategy with respect to tol.

Figure 9 :Figure 10 :
Figure 9: (a) Gamma Γ at time  = 1/365 as a function of the asset price  by closed-form formula, numerical approximation for only one time step size Δ = 1/365, and adaptive time step strategy.(b) Errors Γ − Γ ex as a function of the asset price  by adaptive time-stepping strategy with respect to tol.

Figure 11 :Figure 12 :
Figure 11: (a) Vega ] at time  = 1/365 as a function of the asset price  by closed-form formula, numerical approximation for only one time step size Δ = 1/365, and adaptive time step strategy.(b) Errors ] − ] ex as a function of the asset price  by adaptive time-stepping strategy with respect to tol.
2  (s, ) andFigure 2: Schematic illustration of the numerical solutions (, ) and V(, ) by two different time step sizes in Euler scheme.Here, (, ) denotes the numerical solution by one step with 2Δ and V() denotes the numerical solution by two steps with Δ.
represents the root mean square error (RMSE) of the numerical solution which is measured in the area [0.9, 1.1].Here, RMSE is calculated by √ ∑  =1 (  −  ex  ) 2 /, where  and  ex denote the numerical and exact solutions.In Table2, min Δ 0 and max Δ 0 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,  = 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.

Table 2 :
RMSE by adaptive time stepping strategy with respect to different tolerance.

Table 7 :
Rho  and its error  −  ex at  = 100 and  = 1/365.Here,  ex = 6.82325.timestepmethod.Table8shows RMSE,   , and CPU time with adaptive and uniform time step methods.As shown in Table

Table 8 :
Comparison with adaptive and uniform time step methods.