Exponentially fitted block backward differentiation formulas for pricing options

Abstract A family of Exponentially Fitted Block Backward Differentiation Formulas (EFBBDFs) whose coefficients depend on a parameter and step-size is developed and implemented on the Black–Scholes partial differential equation (PDE) for the valuation of options on a non-dividend-paying stock. Specific EFBBDFs of order 2 and 4 are applied to solve the PDE after reducing it into a system of ordinary differential equations via the method of lines. The methods are shown to be superior to the well-known Crank–Nicolson method since they are -stable and do not exhibit oscillations usually triggered by discontinuities inherent in the payoff function of financial contracts. We confirmed the accuracy of the methods by initially applying them to a prototype example based on the one-dimensional time-dependent convection–diffusion equation with a known analytical solution. It is demonstrated that the American put can be exercised early by computing the hedging parameter “delta”, which specifies the condition for early exercise of the put option. Although the methods can be used to price all vanilla options, we elect to focus on the put due to its optimality.


PUBLIC INTEREST STATEMENT
It is well known that the celebrated Black-Scholes option pricing model provides the theoretical value of European style options and the American call option on a non-dividend-paying stock. Nevertheless, the exact solution cannot be used to value the American put option, which can be optimal when early exercise is elected. Hence, several numerical methods have been proposed in the literature. In this manuscript, a family of Exponentially Fitted Block Backward Differentiation Formulas whose coefficients depend on a parameter and step-size is developed and implemented on the Black-Scholes partial differential equation for the valuation of options on a non-dividend-paying stock. The methods are superior to the popular Crank-Nicolson method since they are stable and do not exhibit oscillations usually triggered by discontinuities inherent in the payoff function of financial contracts. It is demonstrated that the American put can be exercised early via delta, which specifies the condition for early exercise.

Introduction
To make informed decisions, financial analysts use models together with certain parameters to determine the theoretical fair values of assets. The computed values are then compared to the market prices of the assets in question. Both computational procedures and mathematical models abound in the financial economics literature for all kinds of assets including derivative securities. One such celebrated financial model in recent memory is the Black-Scholes (Merton, 1973) option pricing model. A closed form solution of Black-Scholes model does not exist for American style put options due to the possibility of an early exercise. Consequently, several analytical approximations and numerical procedures are employed for the valuation of American put options.
In addition to the procedures described in (Hull, 2015), a great deal of research exists on numerical methods for solving the Black-Scholes differential equation. For example, Khaliq et al. (2006) considered the pricing of an American put option as a free boundary problem and noted that the early exercise feature of the American put option transforms the Black-Scholes linear differential equation into a non-linear type. In another study, Huang et al. (1996) presented a method for valuing and hedging American style options. The study asserted that a "complicated path integral" implicitly defines the early exercise boundary of an American option. It employed a 'unified framework' that made use of an analytic formula and approximation method (Geske & Johnson, 1984). This combined framework was then used to price options on dividend paying stocks by estimating the early exercise boundary for a few points and used the Richardson extrapolation to approximate the entire boundary.
In this paper, we present EFBBDFs, motivated by the fact that the methods have good stability properties, such as A-stability and L-stability. They performed excellently when applied to the so called stiff differential equations. The traditional backward differentiation formulas are implicit and are generally applied in predictor-corrector mode in a forward in time manner, which is called a matching technique. The implementation is generally facilitated by the Newton's method or variable step techniques (Cash, 1980(Cash, , 1984Gear, 1971). It is known that the reduction of parabolic PDEs into a system of first order differential equations triggers a stiff system which can be efficiently solved by L-stable methods. Therefore, we are motivated to solve the Black-Scholes equation (a parabolic PDE) to derive methods that can efficiently solve the Black-Scholes equation via the method of lines (MOLs). In particular, the idea is to convert the Black-Scholes equation into a system of ODEs after which the system is solved using L-stable EFBBDFs that do not exhibit oscillations. It is well known that discontinuities in the payoff function of financial contracts lead to oscillations in the Crank-Nicolson scheme for both the option price and hedging parameters such as delta (Giles & Carter, 2006). This is not the case with the EFBBDFs presented in this paper and according to (Tangman et al., 2008) the commonly used Crank-Nicolson scheme does not exhibit L 0 -stability. Since the methods are L-stable and do not exhibit oscillations, the methods have the advantage of using larger time steps which can be very advantageous when the problem is solved on a wide time interval (Liao & Khaliq, 2009).
Specifically, we derive a family of EFBBDFs based on mixed basis and apply them in a block-byblock manner to solve the Black-Scholes equation. Methods involving mixed basis have also been discussed in the literature (Coleman & Duxbury, 2000;Ndukum et al., 2016;Nguyen et al., 2007). Akinfenwa et al. (2013) discussed extended Block Backward Differentiation formula which was extended to solve a system resulting from a semi-discretization of the Black-Scholes equation (Jator & Nyonna, 2014). We note that the methods are initially applied in a block-by-block fashion to a prototype example based on a forward in time one-dimensional time-dependent convectiondiffusion equation with a known analytical solution. In order to replicate the implementation of the EFBBDFs on the Black-Scholes equation which backward in time, we transform the equation into a forward in time equation. In what follows, we reduce the Black-Scholes equation into a system of ordinary differential equations via the method of lines.
The methods considered in this paper are facilitated by the method of lines approach (Cash, 1984;Lambert, 1991;Ramos & Vigo-Aguiar, 2007) which involves seeking a solution in the strip ½a; b� � ½c; d�, where a; b; c; d are real constants, by first discretizing the variable s with mesh spacings Δs ¼ 1=M, s m ¼ mΔS, m ¼ 0; 1; . . . ; M.
12ðΔsÞ ; @V MÀ 1 ðtÞ @s 12ðΔsÞ ; @ 2 V1ðtÞ We note that (2) or (3) and (4) can be written in the vector form as where f(t, V) = AV+ G, A is an M � M matrix arising from the semi-discretized system (2) or (3) and (4), and G is a vector of constants. The problem (5) is now a system of ordinary differential equations which is solved by EFBBDFs which are L-stable and hence can effectively solve stiff problems. The rest of the paper proceeds as follows. In Section 2, we construct the EFBBDFs and discuss the block formulation and its implementation. In Section 3, the analysis of the methods is given. Section 4 is devoted to numerical examples and the conclusion is given in Section 5.

Derivation of the methods
For notational simplification, we derive the continuous form of the EFBBDF for the scalar form of (5) since the vector form can be solved using the same methods with obvious notational modifications. Motivated by the use of mixed basis functions for the derivation of a collocation method considered by Coleman and Duxbury (2000) and Nguyen et al. (2007), we construct the continuous method using mixed basis functions that belong to the linear space h1; t; . . . ; t kÀ 1 ; e ωt i. In order to solve (5) for a chosen step size h and step number k, we define the EFBBDFs on the partition ft n ¼ . . . ; Ng in which the step ½t n ; V n �7 !½t nþk ; V nþk � is defined by combining the following main and additional methods.
We note that the EFBBDFs are provided by the continuous form and in order to obtain the continuous form, we assume that the solution is given by the function ΦðtÞ on the interval ½t n ; t n þ kh� as where ω is an appropriate parameter and ρ 0 ; ρ 1 ; . . . ; ρ k are coefficients to be uniquely determined by solving the system of equations obtain by imposing the following conditions: The continuous approximation ΦðtÞ is then constructed by substituting the values of ρ j into Equation (7), which is then evaluated at t ¼ t nþj ; j ¼ 0; 1; 2; . . . ; k to obtain the method (6). In what follows, we give two particular methods.

Implementation
Using Mathematica 11, the EFBBDF is implemented in a block-by-block fashion to solve (5) using NSolve[] for linear problems and FindRoot[] for non-linear problems. We recall that the system (5) is solved on the partition N is a constant step-size of the partition of π N , n ¼ 1; 2; . . . ; N, N is a positive integer and n the grid index.
We summarize the process in following algorithm: Algorithm 1 Block-by-Block Algorithm 1: procedure Enter π N , π M , t 0 ; s 0 ; h; N; M; Δs 2: On π M , discretize (1) using (2) or (3) to obtain (5) 3: Then, on π N , discretize (5) using (6) and generate the members of the first block and the variables to be determined for n ¼ 0, say System 4: Solve½System; variables� to obtain the discrete solutions in the first block 5: For n ¼ k; 2k; . . . ; N À k, obtain all solutions on π N 6: end procedure

Order and local truncation error
The local truncation error (LTE) of (12) is defined as Assuming that zðtÞ is a sufficiently differentiable function, the terms in (13) can expanded by Taylor series about the point t n and obtain the expression of the LTE as where the constant coefficients K q , q ¼ 0; 1; . . . are column vectors of size k since the entries of the matrices Ψ 0 ; Ψ 1 ; φ 0 ; and φ 1 are given by the coefficients of the methods in (6).

Numerical examples
The following acronyms are used in the Figures: We begin this section by initially apply the EFBBDFs to a prototype example based on the onedimensional time-dependent convection-diffusion equation with a known analytical solution.
Example 4.1. We consider the following one dimensional convection diffusion equation. This prototype convection-diffusion equation resembles the Black-Scholes equation and is solved using the Crank-Nicolson method and EFBBDFs (k ¼ 2; 4) for a space step size of 0:01 and time step size of 0:25. The results for the errors between the calculated and exact solutions for parameters (ω ¼ 0), (ω ¼ 1), and the optimal parameter for (ω ¼ w n ) are displayed in Figures 2 and  3. It is obvious from Figure 2 that the results obtained using the EFBBDF2 are better than those from the Crank-Nicholson method.
In order to solve the rest of the numerical examples, we note that since (1) is backward in time, we transform (1) into a forward in time equation in the spirit of (Chawla et al., 2003) by letting t ¼ T À τ, τ is a dimensionless variable and Vðs; tÞ ¼ Vðs; T À τÞ ¼ Pðs; τÞ to yield the following corresponding equation.

Figure 2. Errors for the EFBBDF2 and Crank-Nicolson method.
@P @τ ¼ 1 2 σ 2 s 2 @ 2 P @s 2 þ rs @P @s À rP; subject to the final condition: Pðs; 0Þ ¼ maxðX À s; 0Þ, and boundary conditions Pð0; τÞ ¼ X; lim s!1 Pðs; τÞ ¼ 0: Example 4.2. Consider an American put option on a non-dividend-paying stock that has four months to maturity. The exercise price is 21, USD the stock price is 20, USD the risk-free rate of interest is 10 % per annum, and the volatility is 30 % per annum. This example is taken from Hull (2015).

Example 4.3.
Consider an American put option on a non-dividend-paying stock that has one year to maturity. The exercise price is 110, USD the stock price is 100, USD the risk-free rate of interest is 6 % per annum, and the volatility is 40 % per annum. This example is taken from Carr and Hirsa (2003).
Example 4.4. Consider an American put option to sell a Swiss franc for dollars has a strike price of 0.80 USD and a time to maturity of one year. The Swiss franc's volatility is 10%, the dollar interest rate is 6%, the Swiss franc interest rate is 3%, and the current exchange rate is 0.81. USD This example is taken from Hull (2015). In order to compute the put option, we use the standard notations to denote X ¼ 0:80, s ¼ 0:81, r ¼ 0:03, σ ¼ 0:10, Δs ¼ 0:0405, Δt ¼ 0:0833, and T ¼ 1:0, 0 � s � 1:62.  In Figures 4, Figures 7 and Figures 10, we present the results obtained using EFBBDF2 and EFBBDF4 as well as the analytical solution. It is observed in all cases that the two methods accurately approximate the analytical solution.
In Figures 5, Figures 8, and Figures 11, we present a comparison of EFBBDFP2, EFBBDFP4, and the Crank-Nicolson Method. It is well documented in the literature that discontinuity in the payoff function can pose challenges for numerical schemes when they are used to price financial contracts (Le Floch, 2014). The Crank-Nicolson method is very popular for pricing options, but oscillates near the corner on expiry at the exercise price, since it is not L-stable as stated in (Chawla et al., 2003;Tangman et al., 2008). Our methods are L-stable and provide consistently superior approximations to the option than the Crank-Nicolson method, especially, near the corner at the exercise price.     In Figures 6, Figures 9, and Figures 12, we show the variations of delta with stock price for a put option on a non-dividend-paying stock at t ¼ 0. In the spirit of Hull (2015), Greeks are used by traders to measure the risk in an option position so that the risks are acceptable. Several traders use more hedging procedures that involve calculating measures such as delta. The delta of a put option is the first derivative of the put price P with respect to the stock price, s. That is, Δ � @P @s , where À 1<Δ<0. When a put option is deep in-the-money, the strike price, X is far greater than the stock price s, Δ ! À 1. Similarly, when a put option is far out-ofthe-money, the strike price, X is far less than the stock price, s, Δ ! 0. It is well-known that the application of the Crank-Nicolson method to generate Greeks such as delta can distort the Greeks due to oscillations in the scheme (Le Floch, 2014). Therefore, we examine the performances of our methods when applied to compute delta facilitated by the central difference method order 2 given in (1). In Figures 6, Figures 9, and Figures 12 the results obtained using the EFBBDF2 are superior to those obtained from the Crank-Nicolson method, since the EFBBDF2 does not exhibit spurious oscillations in delta.
In Figure 13, early exercise for the American put is discussed in which the EFBBDF2 was used to generate Δ. We compute the option value at each time and stock as well as determine the s value which is used to decide whether early exercise is optimal. In the spirit of Wilmott, Howison, and Dewynne (2009) at each time t there is a particular value of s which indicates the boundary between the "hold" and "exercise" the option regions. Since s depends on time, there could be several values of s for which Δ ¼ @P @s ¼ À 1. We note that the optimal exercise price corresponds to the s ¼ s P value for which the time is minimum and Δ ¼ À 1. In Figure 13 we also show the variation of delta with stock price for the put option on a nondividend-paying stock on the entire grid for Examples 4.2, 4.3, and 4.4.

Conclusion
We have derived and implemented L-stable EFBBDFs that do not produce spurious oscillations in both the option price and hedging parameters such as delta when applied to solve the Black-Scholes PDE for the valuation of options on a non-dividend-paying stock. The characteristics of the methods are discussed and their accuracy confirmed by initially applying them to a prototype example based on the one-dimensional time-dependent convection-diffusion equation with a known analytical solution. It is demonstrated that the American put can be exercised early by computing delta facilitated by incorporating an additional equation based on the central difference method which specifies the condition for early exercise leading to the optimality of the American put. Our future research will be based on extending the methods to solve problems involving multi-assets as well as nonlinear Black-Scholes models.