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

Keywords:

• Block backward differentiation formula
• exponential fitting
• oscillations
• options
• Black–Scholes partial differential equation

Jel codes:

• C63
• G10

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.

1. 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, 1984; Gear, 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-by-block 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 convection–diffusion 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.

1.1. Black–Scholes model and stock variable discretization

Consider the Black–Scholes model for the American put option

(1) $\frac{\mathrm{\partial }V}{\mathrm{\partial }t}+L(V)=0,$(1)

subject to the initial/boundary conditions

$V(0,t)=X$, $V(s,t)=0$ if $s>0$, $V(s,T)=max(X-s,0)$,

where $L=\frac{1}{2}{\sigma }^2{s}^2\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{s}^2}+rs\frac{\mathrm{\partial }}{\mathrm{\partial }s}-r$ is the differential operator, $V(s,t)$ denotes the value of the put option, $\sigma$, $X$, $T$, $r$ denote the volatility of the underlying asset, exercise price, option expiration date, and interest rate, respectively.

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]\times [c,d]$, where $a,b,c,d$ are real constants, by first discretizing the variable $s$ with mesh spacings $\mathrm{\Delta }s=1/M$, $s_m=m\mathrm{\Delta }S$, $m=0,1,\dots ,M$.

We then define

$V_m(t)\approx V(s_m,t),V(t)=[V_0(t),V_1(t),V_2(t),\dots ,V_M(t){]}^T,\mathrm{a}\mathrm{n}\mathrm{d}$

replace the partial derivatives $\frac{{\mathrm{\partial }}^{2}V(s,t)}{\mathrm{\partial }{s}^2}$ and $\frac{\mathrm{\partial }V(s,t)}{\mathrm{\partial }s}$ occurring in (1) by central difference approximations of order 2 or 4. Thus, the central differences of order 2 are given by

(2) $\left\{\begin{array}{c}\frac{\mathrm{\partial }{V}_m(t)}{\mathrm{\partial }s}=\frac{V_{m+1}(t)-V_{m-1}(t)}{2\mathrm{\Delta }s},\\ \frac{{\mathrm{\partial }}^2{V}_m(t)}{\mathrm{\partial }{s}^2}=\frac{V_{m+1}(t)-2V_m(t)+V_{m-1}(t)}{{(\mathrm{\Delta }s)}^2},m=1,\dots ,M-1.\end{array}\right.$(2)

and in the spirit of (Oosterlee et al., 2005) the central differences of order 4 are given by

(3) $\left\{\begin{array}{c}\frac{\mathrm{\partial }{V}_m(t)}{\mathrm{\partial }s}=\frac{-V_{m+2}(t)+8V_{m+1}(t)-8V_{m-1}(t)+V_{m-2}(t)}{12(\mathrm{\Delta }s)},\\ \frac{{\mathrm{\partial }}^2{V}_m(t)}{\mathrm{\partial }{s}^2}=\frac{-V_{m+2}(t)+16V_{m+1}(t)-30V_m,t)+16V_{m-1}(t)-V_{m-2}(t)}{12(\mathrm{\Delta }s{)}^2}.\end{array}\right.$(3)

The following backward differences are required to find the derivative points for $V_1(t)$ and $V_{M-1}(t)$.

(4) $\left\{\begin{array}{c}\frac{\mathrm{\partial }{V}_1(t)}{\mathrm{\partial }s}=\frac{-3V_0(t)-10V_1(t)+18V_2(t)-6V_3(t)+V_4(t)}{12(\mathrm{\Delta }s)},\\ \frac{\mathrm{\partial }{V}_{M-1}(t)}{\mathrm{\partial }s}=\frac{-3V_M(t)-10V_{M-1}(t)+18V_{M-2}(t)-6V_{M-3}(t)+V_{M-4}(t)}{12(\mathrm{\Delta }s)},\\ \frac{{\mathrm{\partial }}^2{V}_1(t)}{\mathrm{\partial }{s}^2}=\frac{10V_0(t)-15V_1(t)-4V_2(t)+14V_3(t)-6V_4(t)+V_5(t)}{12(\mathrm{\Delta }s{)}^2},\\ \frac{{\mathrm{\partial }}^2{V}_{M-1}(t)}{\mathrm{\partial }{s}^2}=\frac{10V_M(t)-15V_{M-1}(t)-4V_{M-2}(t)+14V_{M-3}(t)-6V_{M-4}(t)+V_{M-5}(t)}{12(\mathrm{\Delta }s{)}^2}.\end{array}\right.$(4)

We note that (2) or (3) and (4) can be written in the vector form as

(5) $\frac{dV(t)}{dt}=f(t,V),$(5)

where f($t$, V) = AV+ G, A is an $M\times 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.

2. 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 $⟨1,t,\dots ,t^{k-1},e^{\omega t}⟩$. In order to solve (5) for a chosen step size $h$ and step number $k$, we define the EFBBDFs on the partition $\{t_n=t_0+nh,h=\frac{t_N-t_0}{N},n=0,1,\dots ,N\}$ in which the step $[t_n,V_n]\mapsto [t_{n+k},V_{n+k}]$ is defined by combining the following main and additional methods.

(6) $\left\{\begin{array}{c}\sum _{j=0}^k{Q}_j{v}_{n+j}=h{\mathrm{{\rm Y}}}_k{f}_{n+k},\\ \sum _{j=0}^{k-1}{Q}_{i,j}{v}_{n+j}=h{\mathrm{{\rm Y}}}_{i,k}{f}_{n+k}+h{\mathrm{{\rm Y}}}_{i,i}{f}_{n+i},i=1,\dots ,k-1,\end{array}\right.$(6)

where $Q_j$, $Q_{i,j}$, ${\mathrm{{\rm Y}}}_j$, ${\mathrm{{\rm Y}}}_{i,k}$, ${\mathrm{{\rm Y}}}_{i,i}$ are coefficients. We note that $V_{n+j}$ is the numerical approximation to the analytical solution $V(t_{n+j})$, $f_{n+j}=f(t_{n+j},V_{n+j}),j=0,\dots ,k$. 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 $\mathrm{\Phi }(t)$ on the interval $[t_n,t_n+kh]$ as

(7) $\mathrm{\Phi }(t)=\left[\begin{array}{cccccc}1& t& \hspace{0.17em}& \dots & t^{k-1}& e^{\omega t}\end{array}\right]\left[\begin{array}{c}{\rho }_0\\ {\rho }_1\\ ⋮\\ {\rho }_{k-1}\\ {\rho }_k\end{array}\right],$(7)

where $\omega$ is an appropriate parameter and ${\rho }_0,{\rho }_1,\dots ,{\rho }_k$ are coefficients to be uniquely determined by solving the system of equations obtain by imposing the following conditions:

(8) $\mathrm{\Phi }(t_{n+j})=V_{n+j},j=0,\dots ,k-1,$(8)

(9) $\mathrm{\Phi }{ }^{\mathrm{\prime }}(t_{n+k})=f_{n+k}.$(9)

The continuous approximation $\mathrm{\Phi }(t)$ is then constructed by substituting the values of ${\rho }_j$ into Equation (7)(7) $\mathrm{\Phi }(t)=\left[\begin{array}{cccccc}1& t& \hspace{0.17em}& \dots & t^{k-1}& e^{\omega t}\end{array}\right]\left[\begin{array}{c}{\rho }_0\\ {\rho }_1\\ ⋮\\ {\rho }_{k-1}\\ {\rho }_k\end{array}\right],$(7) , which is then evaluated at $t=t_{n+j},j=0,1,2,\dots ,k$ to obtain the method (6). In what follows, we give two particular methods.

2.1. Particular methods

Case k = 2: When k = 2, the EFBBDF is given by

(10) $\left\{\begin{array}{c}{Q}_2{V}_{n+2}+Q_1{V}_{n+1}+Q_0{V}_n=h{\mathrm{{\rm Y}}}_2{f}_{n+2},\\ Q_{1,1}{V}_{n+1}+Q_{1,0}{V}_n=h({\mathrm{{\rm Y}}}_{1,1}{f}_{n+1}+{\mathrm{{\rm Y}}}_{1,2}{f}_{n+2}),\end{array}\right.$(10)

where $q=\omega h$, $Q_2=1-e^q+e^{2q}q$, $Q_1=-1+e^{2q}-2e^{2q}q$, $Q_0=e^q-e^{2q}+e^{2q}q$,${\mathrm{{\rm Y}}}_2=1-2e^q+e^{2q}$,

$Q_{1,1}=e^{q}q-e^{2q}q$, $Q_{1,0}=-e^{q}q+e^{2q}q$, ${\mathrm{{\rm Y}}}_{1,1}=1-e^q+e^{2q}q$, ${\mathrm{{\rm Y}}}_{1,2}=-1+e^q-e^{q}q$.

Case k = 4: When k = 4, the EFBBDF is given by

(11) $\left\{\begin{array}{c}{Q}_4{V}_{n+4}+Q_3{V}_{n+3}+Q_2{V}_{n+2}+Q_1{V}_{n+1}+Q_0{V}_n=h{\mathrm{{\rm Y}}}_4{f}_{n+4},\\ Q_{1,3}{V}_{n+3}+Q_{1,2}{V}_{n+2}+Q_{1,1}{V}_{n+1}+Q_{1,0}{V}_n=h({\mathrm{{\rm Y}}}_{1,1}{f}_{n+1}+{\mathrm{{\rm Y}}}_{1,4}{f}_{n+4}),\\ Q_{2,3}{V}_{n+3}+Q_{2,2}{V}_{n+2}+Q_{2,1}{V}_{n+1}+Q_{2,0}{V}_n=h({\mathrm{{\rm Y}}}_{2,2}{f}_{n+2}+{\mathrm{{\rm Y}}}_{2,4}{f}_{n+4}),\\ Q_{3,3}{V}_{n+3}+Q_{3,2}{V}_{n+2}+Q_{3,1}{V}_{n+1}+Q_{3,0}{V}_n=h({\mathrm{{\rm Y}}}_{3,3}{f}_{n+3}+{\mathrm{{\rm Y}}}_{3,4}{f}_{n+4}),\end{array}\right.$(11)

where $q=\omega h$,

$Q_4=11-42e^q+57e^{2q}-26e^{3q}+6e^{4q}q$, $Q_3=-18+64e^q-72e^{2q}+26e^{4q}+24e^{4q}q$, $Q_2=9-24e^q+72e^{3q}-57e^{4q}-24e^{4q}q$, $Q_1=-2+24e^{2q}-64e^{3q}+42e^{4q}q$,

$Q_0=e^q(2-9e^q+18e^{2q}+e^{3q}(-11+6q)$, ${\mathrm{{\rm Y}}}_4=6(-1+e^q{)}^4$,

$Q_{1,3}=21+12e^q-33e^{2q}+52e^{q}q+2e^{4q}q$, $Q_{1,2}=-60+27e^q+33e^{3q}-114e^{q}q-12e^{4q}q$, $Q_{1,1}=39-27e^{2q}-12e^{3q}+84e^{q}q+6e^{4q}q$, $Q_{1,0}=-39e^q+60e^{2q}-21e^{3q}-22e^{q}q+4e^{4q}q$, ${\mathrm{{\rm Y}}}_{1,1}=-22+84e^q-114e^{2q}+52e^{3q}-12e^{4q}q$, ${\mathrm{{\rm Y}}}_{1,4}=4-12e^{2q}+2e^{3q}+6e^q+12e^{q}q$,

$Q_{2,3}=8-40e^q+32e^{2q}-26e^{2q}q+2e^{4q}q$, $Q_{2,2}=-4+36e^q-32e^{3q}+57e^{2q}+3e^{4q}q$,

$Q_{2,1}=-4-36e^{2q}+40e^{3q}-42e^{2q}q-6e^{4q}$, $Q_{2,0}=4e^q-8e^{3q}+e^{4q}q+4e^{2q}+11e^{2q}q$,

${\mathrm{{\rm Y}}}_{2,2}=11-42e^q+57e^{2q}-26e^{3q}+6e^{4q}q$, ${\mathrm{{\rm Y}}}_{2,4}=1-6e^q+2e^{3q}+3e^{2q}-6e^{2q}q$,

$Q_{3,3}=23-76e^q+53e^{2q}-52e^{3q}q+22e^{4q}q$, $Q_{3,2}=-28+81e^q-53e^{3q}+114e^{3q}-36e^{4q}q$,

$Q_{3,1}=5-81e^{2q}+76e^{3q}-84e^{3q}q+18e^{4q}q$, $Q_{3,0}=-5e^q+28e^{2q}-23e^{3q}+44e^{3q}q-4e^{4q}q$,

${\mathrm{{\rm Y}}}_{3,3}=22-84e^q+114e^{2q}-52e^{3q}+12e^{4q}q$, ${\mathrm{{\rm Y}}}_{3,4}=-4+18e^q-36e^{2q}+22e^{3q}-12e^{3q}q$.

2.2. Block formulation and Implementation

The methods given in (6) can be expressed in block form as

(12) ${\mathrm{\Psi }}_1{Y}_{\kappa +1}={\mathrm{\Psi }}_0{Y}_{\kappa }+h({\phi }_0{F}_{\kappa }+{\phi }_1{F}_{\kappa +1}),$(12)

where

$Y_{\kappa +1}=(V_{n+1},\dots ,V_{n+k}{)}^{\mathrm{T}},$

$Y_{\kappa }=(V_{n-k+1},\dots ,V_{n-1},V_n{)}^{\mathrm{T}},$

$F_{\kappa +1}=(f_{n+1},\dots ,f_{n+k}{)}^{\mathrm{T}},$

$F_{\kappa }=(f_{n-k+1},\dots ,f_{n-1},f_n{)}^{\mathrm{T}},$

$\kappa =0,\dots ,\varsigma ,n=0,k,2k,\dots ,N=\varsigma k$, $\varsigma$ is the number of blocks in the interval $[t_0,T]$, ${\mathrm{\Psi }}_0,{\mathrm{\Psi }}_1,{\phi }_0,$ and ${\phi }_1$ square matrices of dimension $k$ whose elements are the corresponding coefficients of the methods in (6).

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

${\pi }_M=\{c=s_0<s_1<\dots <s_M=d,s_m=m\mathrm{\Delta }s\},$

$\mathrm{\Delta }s=\frac{d-c}{M}$ is a constant step-size of the partition of ${\pi }_M$, $m=1,2,\dots ,M$, $M$ is a positive integer and $m$ the grid index. Let the partition

${\pi }_N:\{a=t_0<t_1<\dots <t_N=b,t_n=n\mathrm{\Delta }t\},$

$h=\mathrm{\Delta }t=\frac{b-a}{N}$ is a constant step-size of the partition of ${\pi }_N$, $n=1,2,\dots ,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 ${\pi }_N$, ${\pi }_M$, $t_0,s_0,h,N,M,\mathrm{\Delta }s$

2: On ${\pi }_M$, discretize (1) using (2) or (3) to obtain (5)

3: Then, on ${\pi }_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,\dots ,N-k$, obtain all solutions on ${\pi }_N$

6: end procedure

3. Analysis of methods

3.1. Order and local truncation error

The local truncation error (LTE) of (12) is defined as

(13) $L[Z_{\kappa +1};h]={\mathrm{\Psi }}_1{Z}_{\kappa +1}-({\mathrm{\Psi }}_0{Z}_{\kappa }+h({\phi }_0{\bar{F}}_{\kappa }+{\phi }_1{\bar{F}}_{\kappa +1})),$(13)

where

$Z_{\kappa +1}=(z(t_{n+1}),\dots ,z(t_{n+k}){)}^{\mathrm{T}},$

$Z_{\kappa }=(z(t_{n-k+1}),\dots ,z(t_n){)}^{\mathrm{T}},$

${\bar{F}}_{\kappa +1}=(f(t_{n+1},z(t_{n+1})),\dots ,f(t_{n+k},z(t_{n+k})){)}^{\mathrm{T}},$

${\bar{F}}_{\kappa }=(f(t_{n-k+1},z(t_{n-k+1})),\dots ,f(t_n,z(t_n)){)}^{\mathrm{T}}.$

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

(14) $L[Z_{\kappa +1};h]=K_{0}z(t_n)+K_{1}hz{ }^{\mathrm{\prime }}(t_n)+\dots +K_q{h}^q{z}^{(q)}(t_n)+\dots ,$(14)

where the constant coefficients $K_q$, $q=0,1,\dots$ are column vectors of size $k$ since the entries of the matrices ${\mathrm{\Psi }}_0,{\mathrm{\Psi }}_1,{\phi }_0,$ and ${\phi }_1$ are given by the coefficients of the methods in (6).

Definition 3.1. Let $K_0=K_1=\dots =K_p=0$, $K_{p+1}\ne 0$, then the method (12) has order $p\ge 1$ provided there exists a constant $K_{p+1}$ such that the LTE satisfies

$L[Z_{\kappa +1};h]=K_{p+1}{h}^{p+1}+O(h^{p+2}).$

Remark 3.2. In Table 1Table 2, the order and error truncation terms are displayed.

Table 1. Orders and principal local truncation error terms for (10) and (11) generated by evaluating (7) at $t=t_{n+j}$

j	$p$	$K_{p+1}$
1	2	$\frac{-5}{12}(wV{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}(t_n)-V{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}(t_n))$
2	2	$\frac{2}{9}(wV{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}(t_n)-V{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}(t_n))$
1	4	$\frac{29}{390}(w{V}^{(iv)}(t_n)-V^{(v)}(t_n))$
2	4	$\frac{31}{90}(w{V}^{(iv)}(t_n)-V^{(v)}(t_n))$
3	4	$\frac{-111}{1970}(w{V}^{(iv)}(t_n)-V^{(v)}(t_n))$
4	4	$\frac{12}{125}(w{V}^{(iv)}(t_n)-V^{(v)}(t_n))$

Table 2. Optimal exercise price for Examples 4.2, 4.3, and 4.4

Example	$M$	$N$	$t$	$s_P$	$P$
4.2	$20$	$4$	$0.25$	$\mathrm{\$}6$	$\mathrm{\$}14.83$
4.3	$200$	$12$	$0.9167$	$\mathrm{\$}8$	$\mathrm{\$}101.45$
4.4	$40$	$12$	$0.9167$	$\mathrm{\$}0.12$	$\mathrm{\$}0.68$

3.2. Stability

The Linear stability regions is obtained by applying (12) to the test equation $V{ }^{\mathrm{\prime }}=\lambda V$ to give

(15) $Y_{\kappa +1}=M(W;q)Y_{\mu },W=\lambda h,q=wh,$(15)

where the stability matrix $M(W;q))$ is given by

$M(W;q)=({\mathrm{\Psi }}_1-W{\phi }_1{)}^{-1}({\mathrm{\Psi }}_0+W{\phi }_0).$

Definition 3.3. For a fixed q, the method (12) is $A$-stable if for all $W\in {\mathbb{C}}^{-}$, $M(W;q)$ has a dominant eigenvalue $W_{max}$ such that

$|W_{max}|\le 1.$

In particular, its region of absolute stability contains the left half-plane $\{W\in \mathbb{C}|Re(W)<0\}$.

Definition 3.4. For a fixed q, the method (12) is $L$-stable if it is $A$-stability and in addition, $W_{max}\to 0$ as $Re(W)\to -\mathrm{\infty }.$

Specifically, for k = 2, the spectral radius $W_{max}$ of the matrix $M(W;q)$ is a rational function of $W$ given by

$W_{max}=|\frac{e^q(\mu -e^q\mu +q(-1+e^q(1+\mu )))}{e^{2q}(q-\mu )(-1+\mu )-{\mu }^2-e^q(q-\mu )(-1+2\mu )}|,$

and for k = 4,

$W_{max}=|\frac{{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_2}|,$

where

$\begin{array}{rl}& {\mathrm{\Gamma }}_1=e^q(-(-1+e^q)\mu (2(3+3\mu +{\mu }^2)-e^q(12+18\mu +7{\mu }^2)+e^{2q}(6+12\mu +11{\mu }^2))+q(-2(3+\\ & 3\mu +{\mu }^2)-6e^{2q}(3+5\mu +3{\mu }^2)+3e^q(6+8\mu +3{\mu }^2)+e^{3q}(6+12\mu +11{\mu }^2+6{\mu }^3))),\end{array}$

and

$\begin{array}{rl}& {\mathrm{\Gamma }}_2=-6{\mu }^4-6e^{3q}(q-\mu )(-3+7\mu -7{\mu }^2+4{\mu }^3)+e^{4q}(q-\mu )(-6+12\mu -11{\mu }^2+6{\mu }^3)+\\ & 3e^{2q}(q-\mu )(-6+16\mu -19{\mu }^2+12{\mu }^3)-2e^q(q-\mu )(-3+9\mu -13{\mu }^2+12{\mu }^3).\end{array}$

Remark 3.5. In Figure 1, the stability regions are displayed.

Figure 1. The stability regions for the EFBBDF2 (10) and EFBBDF4 (11) plotted in the ($\mu$, q)-plane.

4. Numerical examples

The following acronyms are used in the Figures:

• EFBBDF2 is the Exponentially Fitted Block Backward Differentiation Formula of order 2

• EFBBDF4 is the Exponentially Fitted Block Backward Differentiation Formula of order 4

• EFBBDF2P is the Exponentially Fitted Block Backward Differentiation Formula of order 2 for the put option

• EFBBDF4P is the Exponentially Fitted Block Backward Differentiation Formula of order 4 for the put option

We begin this section by initially apply the EFBBDFs to a prototype example based on the one-dimensional time-dependent convection–diffusion equation with a known analytical solution.

Example 4.1. We consider the following one dimensional convection diffusion equation.

$\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=\frac{1}{2}\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+\frac{1}{2}\frac{\mathrm{\partial }u}{\mathrm{\partial }x},x\in (0,1),t\in (0,1],$

$u(x,0)=e^x,x\in [0,1],$

$u(0,t)=e^t,t\in (0,1],$

$u(1,t)=e^{1+t},t\in (0,1].$

The exact solution $u(x,t)=e^{t+x}$.

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 ($\omega =0$), ($\omega =1$), and the optimal parameter for ($\omega =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.

Figure 2. Errors for the EFBBDF2 and Crank–Nicolson method

Figure 3. Errors for the EFBBDF4

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-\tau$, $\tau$ is a dimensionless variable and $V(s,t)=V(s,T-\tau )=P(s,\tau )$ to yield the following corresponding equation.

$\frac{\mathrm{\partial }P}{\mathrm{\partial }\tau }=\frac{1}{2}{\sigma }^2{s}^2\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{s}^2}+rs\frac{\mathrm{\partial }P}{\mathrm{\partial }s}-rP,$

subject to the final condition: $P(s,0)=max(X-s,0)$,

and boundary conditions

$P(0,\tau )=X,$

${lim}_{s\to \mathrm{\infty }}P(s,\tau )=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).

In order to compute the put option for this example, we use the standard notations to denote $X=21$, s=20, $r=0.10$, $\sigma =0.30$, $\mathrm{\Delta }s=2$, $\mathrm{\Delta }t=0.0833$, and $T=0.3333$.

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). In order to compute the put option, we use the standard notations to denote $X=110$, $s=100$, $r=0.06$, $\sigma =0.40$, $\mathrm{\Delta }s=2$, $\mathrm{\Delta }t=0.0833$, and $T=1.0$.

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$, $\sigma =0.10$, $\mathrm{\Delta }s=0.0405$, $\mathrm{\Delta }t=0.0833$, and $T=1.0$, $0\le s\le 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.Figures 5 and 6

Figure 4. The values of the put option as a function of stock price and time computed using the EFBBDF2 and EFBBDF4 for $M=20$ and $N=4$. The Exact Solution is included for comparison. The value of the put at the stock price of $s=20$ is $\mathrm{\$}1.55$

Figure 5. The value of the put option as a function of the underlying asset, $s$, at $t=0$ computed (a) for $\mathrm{\Delta }s=1/8$ with two time step of $\mathrm{\Delta }t=0.1667$ using the Crank–Nicolson Method and the EFBBDFP2 and (b) computed with four time step for $\mathrm{\Delta }t=0.0833$ using the Crank–Nicolson Method and the EFBBDFP4.

Figure 6. The comparison of $\mathrm{\Delta }$ computed using (a) EFBBDF2 showing no oscillations and (b) Crank Nicolson showing oscillations for $\mathrm{\Delta }s=0.125$, $\mathrm{\Delta }t=0.1667$

Figure 7. The values of the put option as a function of stock price and time computed using the EFBBDF2 and EFBBDF4 for $M=200$ and $N=12$. Exact Solution is included for comparison. The value of the put at the stock price of $s=100$ is $\mathrm{\$}17.98$

Figure 8. The value of the put option as a function of the underlying asset, $s$, at $t=0$ computed (a) for $\mathrm{\Delta }s=2$ with two time step of $\mathrm{\Delta }t=0.5$ using the Crank–Nicolson Method and the EFBBDFP2 and (b) computed with four time step for $\mathrm{\Delta }t=0.25$ using the Crank–Nicolson Method and the EFBBDFP4.

Figure 9. The comparison of methods: (a) Showing no oscillations in $\mathrm{\Delta }$ generating using EFBBDF2 for $\mathrm{\Delta }s=0.125$, $\mathrm{\Delta }t=0.1667$. (b) Showing oscillations in $\mathrm{\Delta }$ generating using the Crank–Nicolson method.

Figure 10. The values of the put option as a function of stock price and time computed using EFBBDF2 and EFBBDF4 for $M=40$ and $N=12$. The Exact Solution is included for comparison. The value of the put at the stock price of $s=0.81$ is $\mathrm{\$}0.0176$

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.

Figure 11. The value of the put option as a function of the underlying asset, $s$, at $t=0$ computed (a) for $\mathrm{\Delta }s=0.0125$ with two time step of $\mathrm{\Delta }t=0.5$ using the Crank–Nicolson Method and the EFBBDFP2 and (b) computed with four time step for $\mathrm{\Delta }t=0.25$ using the Crank–Nicolson Method and the EFBBDFP4.

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, $\mathrm{\Delta }\approx \frac{\mathrm{\partial }P}{\mathrm{\partial }s}$, where $-1<\mathrm{\Delta }<0$. When a put option is deep in-the-money, the strike price, $X$ is far greater than the stock price $s$, $\mathrm{\Delta }\to -1$. Similarly, when a put option is far out-of-the-money, the strike price, $X$ is far less than the stock price, $s$, $\mathrm{\Delta }\to 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.

Figure 12. The comparison of methods: (a) Showing no oscillations in $\mathrm{\Delta }$ generating using EFBBDF2 for $\mathrm{\Delta }s=0.125$, $\mathrm{\Delta }t=0.1667$. (b) Showing oscillations in $\mathrm{\Delta }$ generating using the Crank–Nicolson method.

In Figure 13, early exercise for the American put is discussed in which the EFBBDF2 was used to generate $\mathrm{\Delta }$. 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 $\mathrm{\Delta }=\frac{\mathrm{\partial }P}{\mathrm{\partial }s}=-1$. We note that the optimal exercise price corresponds to the $s=s_P$ value for which the time is minimum and $\mathrm{\Delta }=-1$. In Figure 13 we also show the variation of delta with stock price for the put option on a non-dividend-paying stock on the entire grid for Examples 4.2, 4.3, and 4.4.

Figure 13. Delta for the entire grid computed from EFBBDF2: (a) $M=20$ and $N=4$ for Example, 4.2; (b) $M=200$ and $N=12$ for Example, 4.3; (c) $M=40$ and $N=12$ for Example, 4.4.

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

Acknowledgements

The authors are grateful to the anonymous referees for their constructive suggestions that tremendously improved the quality of the manuscript.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

S. N. Jator

S.N. Jator is the Chair and Professor of Mathematics, Department of Mathematics and Statistics, Austin Peay State University, Clarksville, Tennessee, USA. His research interests include numerical analysis, ordinary and partial differential equations, and computational finance.

R. K. Sahi

R. K. Sahi is the Graduate Coordinator for Mathematical Finance and Professor of Mathematics, Department of Mathematics and Statistics, Austin Peay State University, Clarksville, Tennessee, USA. Her research interest is in risk and interest rate models.

M. I. Akinyemi

M. I. Akinyemi is a Lecturer in the Department of Mathematics and Statistics at the University of Lagos, Lagos, Nigeria. Her research interest is in mathematical finance.

D. Nyonna

D. Nyonna is a Professor of Finance and Chair of the Department of Accounting, Finance, and Economics, Austin Peay State University, Clarksville, Tennessee, USA. His research interest is in financial derivatives.