\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Ninth MSU-UAB Conference on Differential Equations and Computational
Simulations.
\emph{Electronic Journal of Differential Equations},
Conference 20 (2013),  pp. 79--91.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}\setcounter{page}{79}
\title[\hfilneg EJDE-2013/Conf/20/ \hfil Stabilized Adams type method]
{Stabilized Adams type method with a block extension for the
valuation of options}

\author[S. N. Jator, D. Y. Nyonna, A. D. Kerr \hfil EJDE-2013/conf/20 \hfilneg]
{Samuel N. Jator, Dong Y. Nyonna, Andrew D. Kerr}  % in alphabetical order

\address{Samuel N. Jator \newline
Department of Mathematics and Statistics,
Austin Peay State University,
Clarksville, TN 37044, USA}
\email{Jators@apsu.edu}

\address{Dong Y. Nyonna \newline
Department of Accounting, Finance, and Economics,  
Austin Peay State University, Clarksville,
TN 37044, USA}
\email{NyonnaD@apsu.edu}

\address{Andrew D. Kerr \newline
Department of Physics and Astronomy,  
Austin Peay State University, Clarksville,
Clarksville, TN 37044}
\email{akerr@my.apsu.edu}

\thanks{Published October 31, 2013.}
\subjclass[2000]{65L05, 65L06}
\keywords{Stabilized Adams method; extended block; options;
\hfill\break\indent Black-Scholes partial differential equation}

\begin{abstract}
 We construct a continuous stabilized Adams type method (CSAM) that
 is defined for all values of the independent variable on the range
 of interest. This continuous scheme has the ability to provide a
 continuous solution between all the grid points with a uniform
 accuracy comparable to that obtained at the grid points. Hence,
 discrete schemes which are recovered from the CSAM as by-products
 are combined to form a stabilized block Adams type method (SBAM).
 The SBAM is then extended on the entire interval and applied as a
 single block matrix equation for the valuation of options on a
 non-dividend-paying stock by solving a system resulting from the
 semi-discretization of the Black-Scholes model. The stability of the
 SBAM is discussed and the convergence of the block extension of the
 SBAM is given. A numerical example is given to show the accuracy of
 the method.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks

\section{Introduction}

The Black-Scholes option pricing model is one of the most celebrated
achievements in financial economics in the previous four decades. The
model gives the theoretical value of European style options on a
non-dividend-paying stock given the stock price, the strike price,
the volatility of the stock, the time to maturity, and the risk-free
rate of interest. However, since it is optimal to exercise early an
American put option on a non-dividend paying stock, the
Black-Scholes formula cannot be used.
Hull \cite{Hull}
argues that it is never optimal for an American call option on a
non-dividend-paying stock to be exercised early. Therefore the
Black-Scholes formula can be used to value American Style call
options on non-dividend-paying stocks.
In fact, no exact analytic
formula for valuing American put options on non-dividend paying
stocks exists. As a result, numerous numerical procedures are
utilized. A discussion of some of these numerical techniques is
found in Hull \cite{Hull}. In addition to that, several other
numerical procedures for solving the Black-Scholes model abound in
the literature (see Chawla et al. \cite{CAE} and Khaliq et al.
\cite{KVK}). Since there is the possibility of an early exercise,
Khaliq et al. \cite{KVK}) consider the pricing of an American put
option as a free boundary problem. In effect, the early exercise
feature of the American put option transforms the Black-Scholes
linear differential equation into a non-linear type. In order to do
away with the free and moving boundary, Khaliq et al. \cite{KVK} add
a small continuous penalty term to the Black-Scholes equation and
treat the nonlinear penalty term explicitly. They conclude that
their method maintains superior accuracy and stability properties
when compared to standard methods that are based on the Newton-type
iteration procedure in valuing American options.

Furthermore, Chawla et al. \cite{CAE} employ a technique based on
the Generalized Trapezoidal Formulas (GTF) and compare the
computational performance of the scheme obtained  with the
Crank-Nicolson scheme for the case of European option pricing. They
note that their GTF (1/3) scheme is superior to the Crank-Nicholson
scheme. While all these techniques try to accomplish the same goal
by solving the Black-Scholes differential equation for a particular
derivative security, they are applied only after transforming the
model to be forward in time. In this paper, we propose SBAM that is
$A$-stable and applied to solve the model in its original form
without transforming it  into a forward parabolic equation. Thus,
consider the Black-Scholes model
\begin{equation} \label{e1}
\frac{\partial V}{\partial t}
+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2} V}{\partial S^{2}}
+rS\frac{\partial V}{\partial S}-rV=0,
\end{equation}
subject to the initial/boundary conditions
\begin{gather*}
V (0, t) = X, \\
V (S, t) \to 0\quad\text{as }S \to \infty,\\
V (S, T) = \max (X-S, 0),
\end{gather*}
 where $V (S, t)$ denotes the value of the option, $\sigma$
the volatility of the underlying asset, $X$, the exercise price, $T$
the expiry, and $r$ the interest rate.

The method considered in this article involves the method of line
approach to solve \eqref{e1} in which we discretize the space derivatives
in such as way that the resulting system of ordinary differential
equations is stable  (see Lambert \cite{Lam1}, Ramos and Vigo-Aguiar
\cite{RA}, and Cash \cite{CA}). We then discretize time by using the
SBAM. In particular, we seek a solution in the strip (rectangle)
$[a, b]\times[c, d]$ by first discretizing the variable $S$ with
mesh spacings $\Delta S=1/M$,
\[
S_{m}=m\Delta S, \quad m=0, 1, \dots, M.
\]
We then define $v_{m}(t)\approx V(S_{m}, t)$,
$\mathbf{v}(t)=[v_0(t), v_{1}(t), \dots, v_{M-1}(t)]^{T}$, and
replace the partial derivatives
 $\frac{\partial^{2} V(S,t)}{\partial S^{2}}$ and
$\frac{\partial V(S, t)}{\partial S}$
occurring in \eqref{e1} by central difference approximations to obtain
\begin{gather*}
\frac{\partial^{2} V(S_{m}, t)}{\partial S^{2}}
=[v(S_{m+1},t)-2v(S_{m},t)+v(S_{m-1},t)]/(\Delta S)^{2};
\\
\frac{\partial V(S_{m}, ~t)}{\partial S}
=[v(S_{m+1},t)-v(S_{m-1},t)]/(\Delta S),\quad m=0, 1, \dots, M-1.
\end{gather*}
Problem \eqref{e1} then leads to the resulting semi-discrete
problem
\begin{align*}
\frac{d v_{i}(t)}{dt}
&=-\frac{1}{2}\sigma^{2}S_{i}^{2}[v_{i+1}(t)-2v_{i}(t)+v_{i-1}(t)]/(\Delta
S)^{2}\\
&\quad -rS_{i}[v_{i+1}(t)-v_{i-1}(t)]/(\Delta S)+r v_{i}(t)=0,
\end{align*}
 which can be written in the form
\begin{equation} \label{e2}
\frac{d \mathbf{v}(t)}{d t} =\mathbf{f}(t,\mathbf{v}),
\end{equation}
where
$\mathbf{f}(t,\mathbf{v})=\textbf{A v}+ \mathbf{Q}$ and
$\mathbf{A}$ is  an
$M-1 \times M-1$ matrix arising from the central difference
approximations to the derivatives of $S$ and $\mathbf{Q}$ is a vector
of given constants. Problem \eqref{e2} is now a system of ordinary
differential equations which can be solved by the SBAM.

We will assume the scalar form of \eqref{e2} for notational simplification
and will return to the system at the implementation stage in Section
5. We note that the Adams Moulton is one of the most popular methods
available for solving \eqref{e2}. The $k$-Adams Method is given by
\begin{equation} \label{e3}
v_{n+k}-v_{n+k-1}=h\sum _{j=0}^{k}\beta_jf_{n+j}
\end{equation}
where $\beta_{j}$ are constants. 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(t_{n+j}))$, $j=0,\dots ,k$.

For non-stiff problems, \eqref{e3} performs excellently, while for stiff
for problems, \eqref{e3} is restricted by the step-size. For instance, for
$k=4$, \eqref{e3} gives the standard $4-$step Adams-Moulton Method which is
of order 5 and has a stability interval of $[-1.84, 0]$. The method
\eqref{e3} is implemented in a step-by-step fashion in which on the
partition $\pi _N$, an approximation is obtained at $t_n$ only
after an approximation at $t_{n-1}$ has been computed, where
\[
\pi _N:a=t_0<t_{1}<\dots <t_N=b, \quad
t_n=t_{n-1}+h ,\quad n = 1, \dots , N ,
\]
$h=\frac{b-a}{N}$ is the constant step-size of the partition of
$\pi_N$, $N$ is a positive integer, and $n$ is the grid index.

Our objective is to propose a $4$-step stabilized Adams type method
given by
\begin{equation} \label{e4}
v_{n+4}-v_{n+3}=\frac{h}{48}(f_n-2f_{n+1}-4f_{n+2}+34
f_{n+3}+19f_{n+4})
\end{equation}
 with an extended stability interval of $[-6, 0]$.
The method has order 4, since we are more interested
in the stability than the order as the system resulting from the
semi-discretization of \eqref{e1} could be stiff.
The discretization of \eqref{e2} using only
\eqref{e4} leads to an indeterminate, hence we are required to
look for additional methods to complete the system as discussed in
Section 2. The method \eqref{e4} is a by-product of the CSAM constructed in
Section 2. The construction of the CSAM is enhanced by perturbing
the differential by adding a perturbation term involving the shifted
Chebyshev's polynomial $\top^{*}(t)$ of degree 4.

In what follows, we explain how $\top^{*}(t)$ is obtained. We recall
that the trigonometric definition of the Chebychev's polynomials of
degree $m$ are
\[
\top(\xi)=\cos \{m \arccos \xi \},\quad
\xi \in [-1, 1],\quad  m=1,2,3 \dots
\]
which can also be defined by the recurrence relation
$T_0(\xi)=1$, $T_{1}(\xi)=\xi$,
$T_{n+1}(\xi)=2 \xi
T_n(\xi)-T_{n-1}(\xi)$. We then obtain the shifted Chebyshev's
polynomials $\top^{*}(t)$ by transforming the Chebychev's
polynomials $\top(\xi)$ defined on $[-1, 1]$ to the interval
$[t_n, t_{n+4}]$. This is done by using the linear transformation
$ t=d_{1}\xi+d_{2}$, $\xi \in [-1, 1]$, $t \in [t_n,t_{n+4}]$,
where $d_{1}$ and $d_{2}$ are chosen such that $t=t_n$
when $\xi=-1$ and $t=t_{n+4}$ when $\xi=1$ (see  Johnson and Riess
\cite{JL}). It is easily shown that after a simple algebraic
computation $\xi=\frac{2(t-t_n)}{t_{n+4}-t_n}-1$. Thus, we
define the shifted Chebyshev's polynomials of degree $4$ on the
interval $[t_n, t_{n+4}]$ as
\[
\top^{*}(t)=\top(\xi)=\cos [m \cos^{-1}(\frac{2(t-t_n)}{t_{n+4}-t_n}-1)]
\]
which is used to generate $\top^{*}(t_{n+j})$,
$j=0,\dots, 4$ involved with the perturbation term $\eta$ in \eqref{e7}.

The rest of the paper proceeds as follows. In Section 2, we
construct the CSAM and use it to produce the SBAM in Section 3. The
block extension of the SBAM is given in Section 4.  In Section 5,
the implementation of the method to solve the Black-Scholes model is
given. Section 6 is devoted to a numerical example and the
conclusion of the paper is given in Section 7.

\section{Construction of CSAM, and formulation of SBAM}

In this section, we develop a $4$-step SAM for \eqref{e1} on the interval
from $t_n$ to $t_{n+4} = t_n + 4h$, where h is the chosen
step-length. We assume that the solution on the interval
$[t_n, t_{n+4}]$ is initially locally approximated by the polynomial
\begin{equation}\label{e5}
U_{k}(t)=\alpha _{k-1}(t)v_{n+k-1}+h\sum_{j=0}^{k}\beta
_{j}(t)f_{n+j}
\end{equation}
 where $k=4$, $\alpha_3(t)$ and $\beta_{j}(t)$ are
continuous coefficients.

We assume that $v_{n+j} =U_{k}(t_n +jh)$ is the numerical
approximation to the analytical solution $v(t_{n+j})$,
$v'_{n+j} =U_{k}'(t_n +jh)$ is an approximation to
$v'(t_{n+j})$. We also note that
$f_{n+j}=f(t_{n+j}, v_{n+j})$, $j=0,\dots,k$. The
continuous method \eqref{e5} is piecewise continuous on $[a,b]$ and defined
for all $t\in [ a, b]$. That is,  $U_{k}(t)$ is defined
such that $U_{k}(t)=v(t)+O (h^{k+1}))$,
$t \in (t_n, t_{n+k})$. Since $k=4$, the polynomials
$\{U_0(t), U_4(t), \dots, U_{N-4}(t)\}$, then define piecewise
polynomials $U(t)$
which is also continuous on $[a, b]$. Hence, \eqref{e5} has the ability to
provide a continuous solution on $[a,b]$ with a uniform accuracy
comparable to that obtained at the grid points (see \cite{Onum1})
and can also be used to produce additional discrete methods (see
Onumanyi et al. \cite{Onum2}). In what follows, we state the theorem
that facilitates the construction of the CSAM \eqref{e5}.

\begin{theorem} \label{thm2.1}
Let the following conditions be satisfied
\begin{gather} \label{e6}
U_{k}(t_{n+j})=v_{n+k-1},\\
\label{e7}
U_{k}'(t_{n+k})=f_{n+j} +
\eta \top^{*}(t_{n+j}),\quad j=0,\dots,k\,.
\end{gather}
Then, the continuous representations \eqref{e2} and \eqref{e3}
are equivalent to 
\begin{equation} \label{e8}
U_{k}(t)=\Phi_{k}^{T}\left(W_{k}^{-1}\right)^{T}P_{k}(t),
\end{equation}
where
\begin{gather*}
W_{k}=\begin{pmatrix}
P_0(t_{n+3}) & \cdots & P_4(t_{n+3})& 0 \\
P_0'(t_n) & \cdots & P_4'(t_n)& \top_0\\
P_0'(t_{n+1}) & \cdots & P_4'(t_{n+1})& \top_{1} \\
P_0'(t_{n+2}) & \cdots & P_4'(t_{n+2})& \top_{2}\\
P_0'(t_{n+3}) & \cdots & P_4'(t_{n+3})& \top_3\\
P_0'(t_{n+4}) & \cdots & P_4'(t_{n+4})& \top_4\\
\end{pmatrix},
\\
\Phi_{k}=(y_{n+3}, f_n, f_{n+1}, f_{n+2},\dots, f_{n+k})^{T},
\\
P_{k}(t)=(P_0(t),P_{1}(t), \dots, P_{k}(t), 0)^{T}.
\end{gather*}
\end{theorem}

The proof of the above theorem is the same as in Jator et al. \cite{JASF}
with only minor notational modifications.
Note that $T$ denotes the
transpose and $P_{j}(t)=t^{j}, j=0, \dots, k$ are basis functions,
and the perturbation term involves $\eta$ as a parameter and
$\top^{*}(t_{n+j})$ are obtained from $\top^{*}(t)$ which is the
shifted Chebychev's polynomial of degree $4$.

\begin{remark} \label{rmk2.2} \rm
The continuous method \eqref{e5} is obtained by solving a system of
6 equations resulting from conditions \eqref{e6} and \eqref{e7} given in
Theorem \ref{thm2.1}. It is also used to obtain method \eqref{e3} and three
additional methods.
\end{remark}

In particular,  method \eqref{e5} is evaluated at $\{t=t_n, t_{n+1},
t_{n+2}, t_{n+4}\}$ to produce the following methods.
\begin{equation} \label{e9}
\begin{gathered}
v_n-v_{n+3}=\frac{h}{16} (-5 f_n-22 f_{n+1}-12 f_{n+2}-10
f_{n+3}+f_{n+4})\\
v_{n+1}-v_{n+3}=\frac{h}{12} (f_n-8f_{n+1}1-10 f_{n+2}-8 f_{n+3}+f_{n+4})\\
v_{n+2}-v_{n+3}=\frac{h}{48}( f_n-2f_{n+1}-20f_{n+2}-30
f_{n+3}+3f_{n+4})\\
v_{n+4}-v_{n+3}=\frac{h}{48}(f_n-2f_{n+1}-4f_{n+2}+34
f_{n+3}+19f_{n+4})
\end{gathered}
\end{equation}


\subsection{SBAM}

Method \eqref{e9} is used to form the block SBAM as follows:
\begin{equation} \label{e10}
A_{1} Y_{\mu+1}= A_0 Y_{\mu} + h[B_{1}
F_{\mu+1} + B_0F_{\mu}]
\end{equation}
where
\begin{gather*}
Y_{\mu+1}=(v_{n+1}, v_{n+2}, v_{n+3}, v_{n+4})^{T}, \quad
Y_{\mu}=(v_{n-3}, v_{n-2}, v_{n-1}, v_n)^{T},\\
F_{\mu}=(f_{n+1}, f_{n+2}, f_{n+3}, f_{n+4})^{T},\quad
F_{\mu-1}=(f_{n-3}, f_{n-2}, f_{n-1}, f_n)^{T}
\end{gather*}
for $\mu=0, \dots \Gamma$, where $\Gamma=N/4$ is the
number of blocks and $n=0, 4, \dots, N-4$,
and $A_{i}$, $B_{i}$, $i=0,1$ are $4$ by $4$ matrices
whose entries are given by the coefficients of \eqref{e9}.
In particular,
\begin{gather*}
A_0=\begin{pmatrix}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & -1 \\
 0 & 0 & 0 & 0
\end{pmatrix},\quad
A_{1}=\begin{pmatrix}
  1 & 0 & -1 & 0 \\
 0 & 1 & -1 & 0 \\
 0 & 0 & -1 & 0 \\
 0 & 0 & -1 & 1
\end{pmatrix},
\\
B_0=\begin{pmatrix}
  0 & 0 & 0 & \frac{1}{12} \\
 0 & 0 & 0 & \frac{1}{48} \\
 0 & 0 & 0 & -\frac{5}{16} \\
 0 & 0 & 0 & \frac{1}{48}
\end{pmatrix},\quad
B_{1}=\begin{pmatrix}
   -\frac{2 }{3} & -\frac{5 }{6} & -\frac{2 }{3} & \frac{1}{12} \\
 -\frac{1}{24} & -\frac{5 }{12} & -\frac{5}{8} & \frac{1}{16} \\
 -\frac{11 }{8} & -\frac{3}{4} & -\frac{5 }{8} & \frac{1}{16} \\
 -\frac{1}{24} & -\frac{1}{12} & \frac{17 }{24} & \frac{19}{48}
\end{pmatrix}.
\end{gather*}

\subsection*{Consistency of the block SBAM}
The consistency of \eqref{e10} is
accomplished by rewriting it as
\begin{equation} \label{e11}
Y_{\mu+1}= A_{1} ^{-1}A_0 Y_{\mu} + hA_{1} ^{-1}[B_{1}
F_{\mu+1} + B_0F_{\mu}]
\end{equation}
 and define the local truncation error of \eqref{e11} as
\begin{equation} \label{e12}
{L[z(t);h}]= Z_{\mu+1}-( A_{1} ^{-1}A_0
Z_{\mu} + hA_{1} ^{-1}[B_{1} \overline{F}_{\mu+1} +
B_0\overline{F}_{\mu}])
\end{equation}
where
\begin{gather*}
Z_{\mu+1}=(( v(t_{n+1}), v(t_{n+2}), v(t_{n+3}), v(t_{n+4}))^{T},
\\
\overline{F}_{\mu+1}=(( f(t_{n+1}, v(t_{n+2})), f(t_{n+3},
v(t_{n+4})), f(t_{n+1}, v(t_{n+2})), f(t_{n+3}, v(t_{n+4})))^{T},
\\
Z_{\mu}=(v(t_{n-3}), v(t_{n-2}), v(t_{n-1}), v(t_n))^{T},
\\
\overline{F}_{\mu}=(( f(t_{n-3}, v(t_{n-3})), f(t_{n-2},
v(t_{n-2})), f(t_{n-1}, v(t_{n-1})), f(t_n, v(t_n)))^{T},
\end{gather*}
and ${L[z(x);h}]=({L}_{1}[z(x);h], {L}_{2}[z(x);h],
{L}_3[z(x);h], {L}_4[z(x);h])^{T}$ is a linear difference
operator.

Thus, assuming that the arbitrary function $z(t)$ is the exact
solution and is sufficiently differentiable; we can expand the terms
in (12) as a Taylor series about the point $t_n$, to obtain the
expression for the local truncation error (LTE) as
\begin{equation} \label{e13}
L[z(t);h]=O(h^{5}).
\end{equation}
Thus, the method has order four.

\subsection*{Stability}
The Linear-stability regions are obtained by
applying \eqref{e10} to the test equation $y'=\lambda y$ to give
\begin{equation} \label{e14}
Y_{\mu+1}= M(q) Y_{\mu}, \quad q=\lambda h,
\end{equation}
where the stability matrix is
\[
M(q)=(A_{1}-qB_{1})^{-1}(A_0+qB_0)
=\frac{12 + 24 q + 20 q^2 + 8 q^3 + q^4}{12 - 24 q + 20 q^2 - 8 q^3 +
q^4}\,.
\]


\begin{definition} \label{def2.3} \rm
The block method \eqref{e10} is said to be $A$-stability if for all
$q \in  \mathbb C^{-}$, $M(q)$ has a dominant eigenvalue
$q_{\rm max}$ such that $|q_{\rm max}|\leq 1$; moreover, since $q_{\rm max}$ is
a rational function, the real part of the zeros of $q_{\rm max}$ must be
negative, while the real part of the poles of $q_{\rm max}$ must be
positive.
\end{definition}


A simple calculation gives the zeros as $\{ -1 - i,  -1 + i, -3 -
\sqrt{3}, -3 + \sqrt{3}\}$ and the poles as $\{1 - i,  1 + i,  3 -
\sqrt{3}, 3 + \sqrt{3}\}$ (see Figure \ref{fig1}).

\begin{figure}[ht]
\begin{center}
% \includegraphics[width=0.5\textwidth]{fig1} %{SBAM}
\setlength{\unitlength}{.8mm}
\begin{picture}(100,45)(0,0)
\scriptsize
\put(0,20){\line(1,0){100}}
\put(99.7,19.2){$\rightarrow$}
\put(90,15){Real axes}
\put(50,0){\line(0,1){40}}
\put(49.1,40){$\uparrow$}
\put(40,44){Imaginary axes}
\put(1,18.8){$\square$}
\put(36,18.8){$\square$}
\put(40,28.8){$\square$}
\put(40,8.8){$\square$}
\put(62.7,19){+}
\put(96,19){+}
\put(60,29){+}
\put(60,9){+}
\linethickness{0.05mm}
\multiput(51,2)(1,0){49}{\line(0,1){36}}
\end{picture} 
\end{center}
\caption{Stability region for the SBAM confirming A-stability}
\label{fig1}
\end{figure}


\begin{remark} \label{rmk2.5} \rm
Method \eqref{e10} can be used to solve initial value systems in a
block by block fashion as in Jator \cite{JASF}. In this paper, the
approach in \cite{JASF} is modified by first generating the blocks
on the entire interval and then combining these blocks to form a
single block matrix equation which is solved to provide the global
solution of \eqref{e2}.
\end{remark}

\section{Block extension and implementation}

\subsection*{Order Convergence}
 Since the block extension of \eqref{e10} gives
a global method, we discuss the convergence of the block extension
of \eqref{e10}.

\begin{theorem} \label{thm3.1}
Let $Y$, $Z$ be solution vectors formed by extending \eqref{e10} from the
interval $[t_0, t_4]$, to the intervals $[t_4, t_8], \dots,
[t_{N-4}, t_N]$, and $E=Z-Y$, where $Y$ is interpreted as an
approximation of the solution vector for the system formed from the
block extension of \eqref{e10} whose exact solution is $Z$. If
$e_{i}=|v(t_{i})-v_{i}|$, where the exact solution $v(t)$ is several
times differentiable on $[a,b]$ and if $\Vert E\Vert =\Vert
Z-Y\Vert$, then, the SBAM is a fourth-order convergent method. That
is $\Vert E\Vert =O(h^4)$.
\end{theorem}

\begin{proof}
Let the matrices obtained from the block extension of \eqref{e10} be
defined as follows:
\begin{gather*}
A=\begin{pmatrix}
 1 & 0 & -1 & 0 &0 & 0 & 0 & 0 &0 \cdots&0 \\
0 & 1 & -1& 0 &0 & 0 & 0 & 0 & 0 \cdots&0\\
0 & 0 & -1 & 0 & 0& 0 & 0 & 0 & 0 \cdots&0\\
0 & 0 & -1 & 1 & & 0 & 0 & 0 & 0 \cdots&0\\
0 & 0 & 0 & 0 &1& 0 & -1 & 0 & 0\cdots& 0\\
0 & 0 & 0 & 0 &0& 1 & -1 & 0 & 0 \cdots&0\\
0 & 0 & 0 & 1 & 0& 0 & -1 & 0 & 0\cdots&0 \\
0 & 0 & 0 & 0 & 0& 0 & -1 & 1 & 0\cdots&0 \\
\vdots &&&&&&\ddots&&&\vdots\\
0 & 0 & 0 & 0 \cdots& 0&0& 1 & 0 & -1 & 0 \\
0 & 0 & 0 & 0 \cdots& 0&0& 0 & 1 & -1 & 0 \\
0 & 0 & 0 & 0 \cdots& 0& 1& 0 & 0& -1 & 0 \\
0 & 0 & 0 & 0 \cdots&0  &0& 0 & 0 & -1 & 1%
\end{pmatrix},
\\
B=h\begin{pmatrix}
 -\frac{2 }{3} & -\frac{5 }{6} & -\frac{2 }{3} & \frac{1}{12} &0 & 0 & 0 & 0 &0 \cdots&0 \\
-\frac{1}{24} & -\frac{5 }{12} & -\frac{5}{8} & \frac{1}{16} & 0 & 0 & 0 & 0 &0 \cdots&0 \\
-\frac{11 }{8} & -\frac{3}{4} & -\frac{5 }{8} & \frac{1}{16}& 0 & 0 & 0 & 0 &0 \cdots&0 \\
 -\frac{1}{24} & -\frac{1}{12} & \frac{17 }{24} & \frac{19}{48}& 0 & 0 & 0 & 0 &0 \cdots&0 \\
0 & 0 & 0 & \frac{1}{12} &-\frac{2 }{3} & -\frac{5 }{6} & -\frac{2 }{3} & \frac{1}{12}  & 0\cdots& 0\\
0 & 0 & 0 & \frac{1}{48} &-\frac{1}{24} & -\frac{5 }{12} & -\frac{5}{8} & \frac{1}{16}   & 0\cdots& 0\\
0 & 0 & 0 & -\frac{5}{16} &-\frac{11 }{8} & -\frac{3}{4} & -\frac{5 }{8} & \frac{1}{16}  & 0\cdots& 0\\
0 & 0 & 0 & \frac{1}{48} & -\frac{1}{24} & -\frac{1}{12} & \frac{17 }{24} & \frac{19}{48}  & 0\cdots& 0\\
\vdots &&&&&&\ddots&&&\vdots\\
0 & 0 & 0 & 0 \cdots& 0&\frac{1}{12}& -\frac{2 }{3} & -\frac{5 }{6} & -\frac{2 }{3} & \frac{1}{12} \\
0 & 0 & 0 & 0 \cdots&0& \frac{1}{48} & -\frac{1}{24} & -\frac{5 }{12} & -\frac{5}{8} & \frac{1}{16} \\
0 & 0 & 0 & 0 \cdots&0& -\frac{5}{16} & -\frac{11 }{8} & -\frac{3}{4} & -\frac{5 }{8} & \frac{1}{16} \\
0 & 0 & 0 & 0 \cdots&0&  \frac{1}{48} &  -\frac{1}{24} & -\frac{1}{12} & \frac{17 }{24} & \frac{19}{48} \\
\end{pmatrix},
\\
C=(\frac{h}{12}f_0, \frac{h}{48}f_0, -y_0-\frac{5h}{16}f_0,
\frac{h}{48}f_0,0,\dots .,0)^{T}\,.
\end{gather*}
The local truncation errors are
\begin{equation} \label{e15}
\begin{gathered}
\tau_{i+1}=-\frac{13}{180}h^{5}y^{(5)}(t_{i}+\theta_{i})+O(h^{6}),\\
\tau_{i+2}=-\frac{13}{360}h^{5}y^{(5)}(t_{i}+\theta_{i})+O(h^{6}),\\
\tau_{i+3}=-\frac{1}{40}h^{5}y^{(5)}(t_{i}+\theta_{i})+O(h^{6}),\\
\tau_{i+4}=-\frac{17}{360}h^{5}y^{(5)}(t_{i}+\theta_{i})+O(h^{6}),
\end{gathered}
\end{equation}
for $i=0,4,\dots,N-4$, $|\theta_{i}|\leq 1$.
We further define the  vectors
\begin{gather*}
Z=(v(t_{1}),\dots ,v(t_N))^{T},\quad
Y=(v_{1},\dots ,v_N)^{T}, \\
F=(f_{1},\dots ,f_N)^{T},\quad
L(h)=(\tau _{1},\dots ,\tau _N)^{T},
\end{gather*}
 where $L(h)$ is the local truncation error.

Let $E=Y-Z=(e _{1},\dots ,e _N)^{T}$.
The exact form of the system given by the block extension
of \eqref{e10} is
\begin{equation} \label{e16}
AZ-BF(Z)+C+L(h)=0,
\end{equation}
and the approximate form of the system is
\begin{equation} \label{e17}
AY-BF(Y)+C=0,
\end{equation}
where $Y$ is the approximation of the solution vector $Z$.
Subtracting \eqref{e16} from \eqref{e17}, we obtain
\begin{equation} \label{e18}
AE-BF(Y)+BF(Z)=L(h)\,.
\end{equation}
Using the mean-value theorem, we write \eqref{e18} as
\[
(A- BJ)E=L(h),
\]
 where the Jacobian matrix is
\[
J=\begin{pmatrix}
\frac{\partial f_{1}}{\partial v_{1}} & \cdots
& \frac{\partial f_{1}}{\partial v_N} \\
\vdots &  & \vdots \\
\frac{\partial f_N}{\partial v_{1}} & \cdots &
\frac{\partial f_N}{\partial v_N}
\end{pmatrix}.
\]
Let $M = -BJ$ be a matrix of dimension $N$. We have
\begin{equation} \label{e19}
(A+M)E=L(h),
\end{equation}
and for sufficiently small $h$, the matrix $A+M$ is a monotone
and thus invertible (see Jain and Aziz \cite{JAN} and Jator and Li
\cite{JAL}). Therefore,
\begin{equation} \label{e20}
(A+M)^{-1}=D=(d_{i,j})\geq0,\quad
\sum_{j=1}^{N}d_{i,j}=O(h^{-1}).
\end{equation}
If $\|E\|$ is the norm of maximum global error and from
\eqref{e19}, $E=(A+M)^{-1}L(h)$, using \eqref{e20} and the truncation
error vector $L(h)$, it follows that
$\|E\|=O(h^4)$.
Therefore, the TDM is an fourth-order convergent method.
\end{proof}

\subsection{Implementation}

Recall that the semi-discretization of \eqref{e1} is initially performed on
the partition
\[
 \pi _{M}:\{c=S_0<S_{1}<\dots <S_{M}=d,\quad
S_{m}=S_{m-1}+\Delta S\},
\]
where $\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. The resulting system of ODEs \eqref{e2} is
then solved on the partition $\pi _N$. We emphasize the block
extension of \eqref{e10} lead to a single matrix of finite difference
equations, which is solved to provide all the solutions of \eqref{e2} on
the entire grid given by the rectangle $[a, b] \times [c, d]$.
\smallskip

\noindent\textbf{Step 1:}
Use the block extension of \eqref{e10} to
generated from the rectangles $[t_0, t_4]\times [c, d]$, to the rectangle
$[t_4, t_8]\times [c, d], \dots, [t_{N-4}, t_N]\times [c,d]$.
\smallskip

\noindent\textbf{Step 2:}
 Solve the system obtained in step 1, noting that
$v_{m}(t)\approx V(S_{m},~ t)$ and $\mathbf{v}(t_n)=\mathbf{v}_n$, to obtain
$\mathbf{v}_n=[V_{0, n}, V_{1, n}, \dots, V_{M-1, n}]^{T}$,
$ n=1, 2, \dots, N$.
\smallskip

\noindent\textbf{Step 3:}
The solution of \eqref{e1} is approximated by the
solutions in step 2 as  $\mathbf{v}(t_n)=\mathbf{v}_n$, where
$\mathbf{v}(t_n)=[V(S_0,t_n), V(S_{1},t_n), \dots,
V(S_{M-1},t_n)]^{T}$, $n=1, 2, \dots, N$.

\section{Numeriacl examples}
Our method was implemented Mathematica 8.0
enhanced by the feature NSolve[\;].

\begin{example}\label{examp1} \rm
Consider a five-month European call and put options on
a non-dividend-paying stock when the stock price is \$50, the strike
price is \$50, the risk-free interest rate is 10\% per annum, and
the volatility is 40\% per annum. This example is taken from Hull
\cite{Hull}.
\end{example}

To compute the call and put options we use the standard notation to
denote $X=50$, $S=50$, $r=0.10$, $\sigma=0.40$, and $T=0.4167$.
The theoretical solutions for the prices of the European call and
put options are given in Hull \cite{Hull} as follows.
\begin{gather*}
c=S N(d_{1})-Xe^{-r(T-t)}N(d_{2}),\\
p=Xe^{-r(T-t)}N(-d_{2}-S N(-d_{1}),
\end{gather*}
 where
\begin{gather*}
d_{1}=\frac{\ln(S/X)+(r+\sigma^{2}/2)(T-t)}{\sigma\sqrt{T-t}},\\
d_{2}=\frac{\ln(S/X)+(r-\sigma^{2}/2)(T-t)}{\sigma\sqrt{T-t}}
=d_{1}-\sigma\sqrt{T-t},
\end{gather*}
and  $N(x)$ is the cumulative probability distribution
function for the standard normal variable.

We note that Example \ref{examp1} was chosen to illustrate that the SBAM does
not exhibit oscillatory behavior near the exercise price which occurs
in well known methods such as the Crank-Nicolson method
(see \cite{CAE} and \cite{KVK}).
In Figures \ref{fig2} and \ref{fig3}, we plot the values of the call
and put options with their corresponding exact solutions at expiration
 versus $S$ respectively, while in Figures \ref{fig4} and \ref{fig5}
 we plot
the the call and the put options with their corresponding exact solutions
versus $S$ and $t$. In all cases it is obvious that the SBAM gives an
oscillation-free solution as the graphs are smooth through out.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2} % {CurveCall}
\end{center}
 \caption{Call curves (C) at $t=0$ for Example \ref{examp1},
where ECALL is the exact solution and SBAMC is the call option for
the current method}
\label{fig2}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3} %{CurvePut}
\end{center}
\caption{Put curves (V) at $t=0$ for Example \ref{examp1},
where EPUT is the exact solution and SBAMP is the put option for the
current method}
\label{fig3}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.46\textwidth]{fig4a} \quad % SBAMP3D
\includegraphics[width=0.46\textwidth]{fig4b} % EPUT3D
\end{center}
\caption{Approximate and exact solutions for the put option for
Example \ref{examp1} with $h=1/12$, $\Delta S=20$}
\label{fig4}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.46\textwidth]{fig5a} \quad % SBAMC3D
\includegraphics[width=0.46\textwidth]{fig5b} % ECALL3D
\end{center}
\caption{Approximate and exact solutions for the call option for
Example \ref{examp1} with $h=1/12$, $\Delta S=20$}
\label{fig5}
\end{figure}

\begin{example}\label{examp2} \rm
As our second test example, we solve the  stiff parabolic equation
(see Cash\cite{CA})
\begin{gather*}
\frac{\partial V}{\partial t} = \kappa \frac{\partial^{2} V}{\partial S^{2}}\\
u(0,t)=u(1,t), \quad u(S,0)=\sin\pi S+\sin \omega \pi S, \quad \omega\gg 1.
\end{gather*}
The exact solution
$V(S, t)=e^{-\pi^{2} \kappa t}\sin\pi S+e^{-\omega^{2} \pi^{2}
\kappa t}\sin\omega\pi S$.
\end{example}

 Cash \cite{CA} notes that as $\omega$ increases, equations of the
type given in Example \ref{examp2} exhibit characteristics similar to model
stiff equations. Hence, the methods such as the Crank-Nicolson
method which are not $L_0$-stable are expected to perform poorly.
However, we found that the SBAM is not $L_0$-stable and performs
well when applied to this problem. In Table \ref{table1}, we display the results
for $\kappa=1$ and a range of values for $\omega$.

\begin{table}[ht]
\begin{center}
\begin{tabular}{cc}
\hline 
$\omega$& SBAM \\
\hline
1&  $3.23\times10^{-6}$\\
2&  $1.85\times10^{-5}$\\
10& $1.62\times10^{-6}$\\
20& $1.62\times10^{-6}$\\
30& $1.62\times10^{-6}$\\
40& $1.62\times10^{-6}$\\
\hline\\
\end{tabular}
\end{center}
\caption{Errors of the method for Example
\ref{examp2} at $t=1$ and $\Delta S=1/10$, $h=1/12$}
\label{table1}
\end{table}

To test for convergence, Example \ref{examp2} was solved for various values
of $h=\Delta S$ and the results for the global maximum absolute
errors ($\mathrm{Err} = \max |v_{m}(t_n)- V(S_{m}, t_n)|)$ are reproduced
in Table \ref{table2}. We also give the rate of convergence (ROC)
which is calculated using the formula
$\mathrm{ROC}= \log_{2}(\mathrm{Err}^{2h}/\mathrm{Err}^{h})$, $\mathrm{Err}^{h}$
is the error obtained
using the step size $h$. In general, the ROC shows that the order of
the method is slightly greater than 2. This is expected since the
central difference method used for the spatial discretization
is of order 2 and hence affects the convergence of the SBAM
which is of order 4 with respect to the time variable.
In Figure \ref{fig6}, the solutions obtained using the SBAM are plotted
versus $S$ and $t$ and are in good agreement with the plots given by
the exact solution. Furthermore, In Figure \ref{fig7}, the errors (residuals)
at all the grid points are generally small when the residuals are
plotted versus $S$ and $t$.

\begin{table}[ht]
\begin{center}
\begin{tabular}{ccc}
\hline
$N=M$& Err & ROC \\ \hline
4  &$1.15\times10^{-1}$& \\
8  &$1.92\times10^{-2}$&2.6\\
16 &$1.62\times10^{-3}$&2.5\\
32 &$6.82\times10^{-4}$&2.4\\
64 &$1.48\times10^{-4}$&2.2\\
128&$3.67\times10^{-5}$&2.0\\
256&$9.23\times10^{-6}$&2.0\\
\hline \\
\end{tabular}
\end{center}
\caption{ Err and ROC for Example \ref{examp2}}
\label{table2}
\end{table}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.46\textwidth]{fig6a} \quad % SBAMstiff
\includegraphics[width=0.46\textwidth]{fig6b} % EXACTstiff
\end{center}
\caption{Approximate and exact solutions for Example \ref{examp2} with
$h=1/128$, $\Delta S=1/128$}
\label{fig6}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig7} % SBAMstiffR}
\end{center}
\caption{Residuals for $S=t=1/128$ for Example \ref{examp2}}
\label{fig7}
\end{figure}

\subsection*{Conclusion}
An SBAM with a block extension for solving the the Black-Scholes
differential equation has been proposed.
The method is shown to give oscillation-free solutions as
 illustrated in Figures \ref{fig2}--\ref{fig5}.
The method is also shown to be convergent
and capable of solving stiff problems (see Table \ref{table2}.


\begin{thebibliography}{00}

\bibitem{CA} J. R. Cash;
 \emph{Two new finite difference schemes for parabolic equations},
SIAM J. Numer. Anal., 21 (1984) 433-446.

\bibitem{CAE} M. M. Chawla, M. A. Al-Zanaidia, D. J. Evans;
\emph{Generalized Trapezoidal Formulas for the Black-Scholes
Equation of Option Pricing}, International Journal of Computer
Mathematics, 80 (2003) 1521-1526.

\bibitem{Hull} J. C. Hull;
 \emph{Options, Futures, and Other Derivatives, Fourth Edition},
Prentice Hall, Upper Saddle River, NJ, 2000.

\bibitem{JAN} M. K. Jain, T. Aziz;
\emph{Cubic spline solution of two-point boundary
Value with signifigant first derivatives}, Comp. Meth. Appl. Mech.
Eng., 39 (1983) 83-91.

\bibitem{JASF} S.N. Jator, S. Swindle, R. French;
\emph{Trigonometrically fitted block Numerov type method for
 $y''=f(x,y,y') $}, Numer Algor., (2012): DOI
10.1007/s11075-012-9562-1.

\bibitem{JAL} S. N. Jator, J. Li;
\emph{An algorithm for second order initial
and boundary value problems with an automatic error estimate based
on a third derivative method}, Numerical Algorithms, 59 (2012)
333-346.

\bibitem{JL} L. W. Johnson, R. D. Riess;
 \emph{Numerical Analysis} (second edition),
 Addison-Wesley, Massachusetts, USA, 1982.

\bibitem{KVK} A. Q. M. Khaliq, D. A. Voss, S. H. K. Kazmi;
\emph{A linearly implicit predictor–corrector scheme for pricing
American options using a penalty method approach}, Journal of
Banking \& Finance, 30 (2006) 489-502.

\bibitem{Lam1} J. D. Lambert;
 \emph{Numerical methods for ordinary differential
systems,} John Wiley, New York, 1991.

\bibitem{Onum2} P. Onumanyi, U. W.  Sirisena, S. N. Jator;
\emph{Continuous finite difference approximations for solving differential
equations}, Inter. J. Compt. Maths. 72 (1999) 15-27.

\bibitem{Onum1} P. Onumanyi, D. O. Awoyemi, S. N. Jator,  U. w.  Sirisena;
\emph{New linear mutlistep methods with continuous coefficients for first
order initial value problems}, J. Nig. Math. Soc. 13 (1994) 37-51.

\bibitem{RA} H. Ramos, J. Vigo-Aguiar;
\emph{A fourth-order Runge-Kutta method based on BDF-type
Chebyschev approximations}, J. Comp. Appl. Math. 204 (2007) 124-136.

\end{thebibliography}

\end{document}