An Alternative Numerical Scheme to Approximate the Early Exercise Boundary of American Options

Abstract

This paper deals with a new numerical method for the approximation of the early exercise boundary in the American option pricing problem. In more detail, using the mean-value theorem for integrals, we provide a flexible algorithm that allows for reaching a more accurate numerical solution with fewer calculations rather than other previously described methods.

1. Introduction

Volterra integral equations have practical applications in many fields, including biology, medicine, finance and engineering. Although various methods are available to find analytical solutions for these integral equations, in most cases, finding a closed-form solution may be unfeasible. To this end, several numerical methods have been developed in recent years (see, e.g., [1,2]). In the present work, we deal specifically with the Volterra integral equation of the second kind, defined below

$$\varphi \left(t\right)=f(t,\varphi \left(t\right))+{\int }_0^{t}k(t,s)\psi (t,s,\varphi \left(t\right),\varphi \left(s\right))ds,$$

(1)

with $t\in I=[0,T]$. Such an equation in which the integrand depends both on $\varphi \left(t\right)$ and $\varphi \left(s\right)$ is also a nonstandard Volterra integral equation ([3]).

In the field of mathematical finance, Equation (1) and its solution are a matter of interest in order to determine the value of an American financial option. More generally, a financial option is a derivative contract allowing for the holder to buy or sell an underlying financial asset at a fixed price, namely, the strike price. If the holder can buy or sell the underlying asset only at the expiration date, the option is called European; conversely, an American financial option provides the holder with the possibility to also exercise its right before the expiration date. Depending on the value of the underlying asset, it may become optimal for the holder to exercise its right before the expiration time. Such an asset value is the optimal exercise price, and by considering the optimal exercise times during the option duration, we can achieve a collection of optimal exercise prices, namely, the optimal exercise boundary. From a mathematical perspective, the definition of optimal exercise boundary identifies a free boundary problem whose solution must be determined numerically. For instance, Kim [4] showed that the early exercise boundary of an American put option is the solution of the integral equation having the form (1).

Among others, the authors in [5,6], Barone-Adesi, Whaley [7], Bunch et al. [8], Aitsahlia and Lai [9,10] considered integral equations and/or the optimal stopping problem for pricing American options. The integral representation proposed by Kim gives a suggestive economic interpretation, but a closed-form solution for these kind of problems is not available. In particular, Kim’s equation ([4]) contains a double integral due to the presence of a cumulative normal distribution term. To avoid such a problem, the authors in [11,12] suggested a transformation of the original equation in order to reduce it to an one-dimensional integral equation and afterwards consider its numerical solution. The early exercise boundary shows some kind of singularity near expiration (see, e.g., [13,14]). Considering the numerical approximation problem of the early boundary exercise for the American put, there have been more efforts in recent years. The authors in [15] proposed a four-point extrapolation scheme. Bunch and Johnson [8] proposed a modified version of the two-point extrapolation scheme of Geske and Johnson. Broadie and Detemple [16] suggested a lower and upper bound approximation method. The authors in [17] showed the exponential function method combined with a three-point Richardson extrapolation. Starting from an integral equation representation for the early boundary exercise, the authors in [12] used an iteration method. The authors in [18] proposed trapezoidal formulas approximations combined with the Newton–Raphson iteration. In [19], the authors solved a variational inequality representation of the American put option pricing problem by applying the projected successive over-relaxation method. The authors in [20] suggested an analytical expression to the value of American put options and their optimal exercise boundary. The authors in [21] derived a local iterative numerical scheme based on a solution of the integral equation proposed by [13]. More recently, Nedaiasl et al. [22], starting from the cited one-dimensional reformulation of Kim’s integral equations, obtained the numerical solution by using a modified version of the Nyström method in order to take into account the singularity near expiry presented by the early exercise boundary. Besides the method of Zhu [20] and others involving semianalytical approximations that are not cited for the brevity of treatment, in the American option pricing literature, numerical methods aim to numerically solve the integral equations describing the early exercise boundary. We proceed in the latter direction.

More recently, some research studies applied the integral equation approach to deal with particular financial evaluations. For instance, the authors in [23] exploited the Mellin’s transform aiming to provide an integral representation for the barrier option price. In [24], the authors gave an integral equation representation for a two-free-boundaries problem arising in the American better-of option on two assets.

In the present work, using a very simplified approach, we straightforward solve the nonstandard Volterra integral equation of the form (1). In more detail, using the mean-value theorem for integrals, we provide a flexible algorithm that allows for reaching a more accurate numerical solution with fewer calculations rather than other previously described methods.

The paper is organized as follows. In Section 2, we propose our numerical method for nonstandard Volterra integral equations. In Section 3, such a method is applied to the case of the American put. In Section 4, we present some numerical results confirming the accuracy of the proposed method. Lastly, in Section 5 we show the conclusions.

2. A New Numerical Scheme for Nonstandard Volterra Integral Equations

Let us fix interval $I=[0,T]$ where $T>0$ and consider the following nonstandard Volterra integral equation:

$$\varphi \left(t\right)=f(t,\varphi \left(t\right))+{\int }_0^{t}k(t,s)\psi (t,s,\varphi (t),\varphi (s))ds,$$

(2)

with $t\in I=[0,T]$, where the kernel function k is continuous in $I\times I$, $\psi$ is continuous in $I\times I\times \mathbb{R}\times \mathbb{R}$.

Let n be a positive integer, and let us consider the following partition $\Gamma$ of the interval $[0,T]$ into n intervals of equal length $\Delta =T/n$:

$$0=t_0<t_1<t_2<\cdots <t_{n-1}<t_n=T.$$

(3)

Let us consider Equation (2). By means of the additive properties for integrals, we can rewrite for each $t_i$, $i=1,2,\cdots ,n$ Equation (2) in the following way:

$$\varphi \left(t_i\right)=f(t_i,\varphi \left(t_i\right))+\sum _{m=1}^i{\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds.$$

(4)

The following result holds.

Proposition 1. 

Let:

(a) 

Function $k(t,s)$ be continuous such that it never changes sign for all $(t,s)\in I\times I$. Let $L>0$ be a constant such that $\left|k\right(t,s\left)\right|\le L$ for each $(t,s)\in I\times I$.

(b) 

Function $\varphi \left(s\right)$ be continuous in I.

(c) 

Function $\psi (t,s,x,y)$ be continuous in $I\times I\times \mathbb{R}\times \mathbb{R}$.

Then, for each $m=1,2\cdots ,i$, the following approximation holds:

$$\begin{array}{c}{\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds=\hfill \\ =\frac{1}{2}\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\\ \hfill +\frac{1}{2}\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+{ϵ}_{m,i},\end{array}$$

(5)

where error $|{ϵ}_{m,i}|$ satisfies $|{ϵ}_{m,i}|\le h_m/\left(2n\right)$, with $h_m$ a constant, and ${ϵ}_{m,i}\to 0$ for $n\to \infty$.

Proof. 

Using the mean value theorem for integrals, it follows that

$${\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds=\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds,$$

where ${\xi }_m$ is, for each m, a number belonging to the interval $[t_{m-1},t_m]$. It follows

$$\begin{array}{c}{\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds=\frac{1}{2}\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\hfill \\ +\frac{1}{2}\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds-\frac{1}{2}\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\\ -\frac{1}{2}\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\frac{1}{2}\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\\ \hfill +\frac{1}{2}\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds.\end{array}$$

(6)

Let us rearrange the elements in (6) in the following way:

$$\begin{array}{c}{\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds=\frac{1}{2}\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\hfill \\ \hfill +\frac{1}{2}\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+{ϵ}_{m,i}^A+{ϵ}_{m,i}^B,\end{array}$$

(7)

where we have set

$${ϵ}_{m,i}^A=\frac{1}{2}\left[\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right))\right]{\int }_{t_{m-1}}^{t_m}k(t_i,s)ds$$

and

$${ϵ}_{m,i}^B=\frac{1}{2}\left[\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right))\right]{\int }_{t_{m-1}}^{t_m}k(t_i,s)ds.$$

Then, it follows that

$$\begin{array}{c}|{ϵ}_{m,i}^A|=\frac{1}{2}\left|\left[\psi (t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right))\right]{\int }_{t_{m-1}}^{t_m}k(t_i,s)ds\right|=\hfill \\ \frac{1}{2}\left|\psi \right(t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)\left)\right|·\left|{\int }_{t_{m-1}}^{t_m}k(t_i,s)ds\right|\le \\ \frac{1}{2}\left|\psi \right(t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)\left)\right|·{\int }_{t_{m-1}}^{t_m}\left|k(t_i,s)\right|ds\le \\ \hfill \frac{1}{2}\left|\psi \right(t_i,{\xi }_m,\varphi \left(t_i\right),\varphi \left({\xi }_m\right))-\psi (t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)\left)\right|·\frac{TL}{n}=\frac{h_{m,i}^A}{2n},\end{array}$$

where $h_{m,i}^A$ is a positive constant because the argument in modulus is a number for each m and for any fixed $t_i$.

As $n\to \infty$, it follows ${ϵ}_{m,i}^A\to 0$, and the same considerations apply to ${ϵ}_{m,i}^B$ with a positive constant $h_{m,i}^B$. If we consider the sum ${ϵ}_{m,i}={ϵ}_{m,i}^A+{ϵ}_{m,i}^B$, it follows ${ϵ}_{m,i}\to 0$, as $n\to \infty$.

In addition, it follows easily that $|{ϵ}_{m,i}|\le h_{m,i}/\left(2n\right)$, with $h_{m,i}=h_{m,i}^A+h_{m,i}^B$. □

Remark. 

If we consider Equation (4), it follows, through Proposition 1 for each fixed $t_i$, that

$$\begin{array}{c}\sum _{m=1}^i{\int }_{t_{m-1}}^{t_m}k(t_i,s)\psi (t_i,s,\varphi \left(t_i\right),\varphi \left(s\right))ds=\hfill \\ =\frac{1}{2}\sum _{m=1}^i\left[\psi \right(t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s),ds+\\ \hfill +\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s),ds]+{ϵ}_i,\end{array}$$

where

$${ϵ}_i=\sum _{m=1}^i{ϵ}_{m,i}.$$

From Proposition 1, it is easy to check that: $|{ϵ}_i|\le H/\left(2n\right)$, with $H={\sum }_{m=1}^i{h}_{m,i}$. It is straightforward to notice that ${ϵ}_i\to 0$ for $n\to \infty$.

We provide the following algorithm in order to find the numerical solution.

Step 1. 

Let n be a positive integer. Let us consider the following partition Γ of interval $[0,T]$ into n intervals of equal length $\Delta =T/n$:

$$0=t_0<t_1<t_2<\cdots <t_{n-1}<t_n=T.$$

If Proposition 1 is verified, and from the additive property for integrals, for each $t_i$, with $i=0,1,2,\cdots ,n$, we can rewrite Equation (2) in the following way:

$$\begin{array}{c}\hfill \varphi \left(t_i\right)=f(t_i,\varphi \left(t_i\right))+\frac{1}{2}\left[\sum _{m=1}^i\psi \right(t_i,t_{m-1},\varphi \left(t_i\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds+\\ \hfill +\psi (t_i,t_m,\varphi \left(t_i\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_i,s)ds]\end{array}$$

(8)

Step 2. 

We observe that

$$\varphi \left(t_0\right)=f(t_0,\varphi \left(t_0\right))+{\int }_{t_0}^{t_0}k(t_0,s)\psi (t_0,s,\varphi \left(t_0\right),\varphi \left(s\right))ds=f(t_0,\varphi \left(t_0\right)).$$

(9)

Then, for $i=0,\cdots ,n$ we consider the following nonlinear system.

$$\left\{\begin{array}{c}\varphi \left(t_0\right)=f(t_0,\varphi \left(t_0\right)),\hfill \\ \varphi \left(t_1\right)=f(t_1,\varphi \left(t_1\right))+\frac{1}{2}\left[\psi \right(t_1,t_0,\varphi \left(t_1\right),\varphi \left(t_0\right)){\int }_{t_0}^{t_1}k(t_1,s)ds+\hfill \\ \\ +\psi (t_1,t_1,\varphi \left(t_1\right),\varphi \left(t_1\right)){\int }_{t_0}^{t_1}k(t_1,s)ds],\hfill \\ \\ \varphi \left(t_2\right)=f(t_2,\varphi \left(t_2\right))+\frac{1}{2}\left[\sum _{m=1}^2\psi \right(t_2,t_{m-1},\varphi \left(t_2\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_2,s)ds+\hfill \\ \\ +\psi (t_2,t_m,\varphi \left(t_2\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_2,s)ds],\hfill \\ ⋮\hfill \\ \varphi \left(t_n\right)=f(t_n,\varphi \left(t_n\right))+\frac{1}{2}\left[\sum _{m=1}^n\psi \right(t_n,t_{m-1},\varphi \left(t_n\right),\varphi \left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}k(t_n,s)ds+\hfill \\ +\psi (t_n,t_{m-1},\varphi \left(t_n\right),\varphi \left(t_m\right)){\int }_{t_{m-1}}^{t_m}k(t_n,s)ds].\hfill \end{array}\right.$$

(10)

The above nonlinear system is solved with a numerical method that provides the approximate solution $\{\tilde{\varphi }\left(t_0\right),\tilde{\varphi }\left(t_2\right),\cdots ,\tilde{\varphi }\left(t_n\right)\}$.

3. Approaching the American Put Pricing Problem

Let us assume the asset price process $\left\{S\right(t),t\ge 0\}$ following the log-normal distribution of the form

$$d{S}_t=(r-\delta )S_{t}dt+\sigma S_{t}d{W}_t,$$

(11)

where $W_t$ is the standard Wiener process, r is the constant interest rate, and $\delta$ is the dividend yield. For the aim of the present work, we follow [25] assuming that the volatility term, $\sigma$, is constant. The latter is a simplifying hypothesis, and Equation (11) can be generalized considering nonconstant volatility, as in [26]. In this case, the boundary problem is represented by a multidimensional Volterra integral equation involving more complicated integrals (see, e.g., [27]). We leave the nonconstant volatility problem to future research works.

Let us consider the interval $I=[0,T]$. We denote with $\mathsf{\Phi }(·)$ the standard cumulative normal distribution function. Let us denote with $\mathcal{B}\left(t\right)$ the early exercise boundary of an American put option, where $t\in I$. Let P the price of an American put. The authors in [28] proved the existence and uniqueness of the pair $(P,\mathcal{B})$, and the continuity and monotonic behaviour of $\mathcal{B}$. In [4], the authors showed that $\mathcal{B}\left(t\right)$, as a function of time to expiration, satisfies the following weakly singular Volterra integral equation:

$$\begin{array}{c}\hfill K-\mathcal{B}\left(t\right)=K{e}^{-rt}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{K}\right)+(r-\delta -\frac{{\sigma }^2}{2})t}{\sigma \sqrt{t}})+\\ \hfill -\mathcal{B}\left(t\right)e^{-\delta t}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{K}\right)+(r-\delta +\frac{{\sigma }^2}{2})t}{\sigma \sqrt{t}})+\\ \hfill +Kr{\int }_0^t{e}^{-r(t-s)}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(s\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}})ds+\\ \hfill -\delta \mathcal{B}\left(t\right){\int }_0^t{e}^{-\delta (t-s)}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(s\right)}\right)+(r-\delta +\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}})ds,\end{array}$$

(12)

where K is the exercise price.

In Equation (12), let us consider the two integrals and in particular the two functions $k_1(t,s)=e^{-r(t-s)}$ and $k_2(t,s)=e^{-\delta (t-s)}$. It is easy to check that $k_1(t,s)$ and $k_2(t,s)$, for any $(t,s)\in I\times I$, satisfy Proposition 1 with $L=1$.

We apply the algorithm presented in Section 2. Let us consider the following result.

Proposition 2. 

Assume that the asset price, $S_t$, follows a log-normal distribution process of the form $d{S}_t=(r-\delta )S_{t}dt+\sigma S_{t}d{W}_t$, in which $W_t$ is the standard Wiener process.

Let $\mathcal{B}\left(t\right)$ be the solution of Integral Equation (12). Then, $\mathcal{B}\left(t\right)$ is a continuously differentiable function on $(0,T]$ and

• for $r\le \delta$:

  $$\underset{t\to 0}{lim}\mathcal{B}\left(t\right)=K;$$

  (13)
• for $r>\delta$:

  $$\underset{t\to 0}{lim}\mathcal{B}\left(t\right)=\frac{r}{\delta }K.$$

  (14)

Proof. 

See [4]. □

For our purposes, let us rewrite Equation (12) in contract form:

$$\begin{array}{c}\hfill \mathcal{B}\left(t\right)=K-f(t,\mathcal{B}\left(t\right))-rK{\int }_0^t{e}^{-r(t-s)}\mathsf{\Phi }(-d_2(t,s,\mathcal{B}\left(t\right),\mathcal{B}\left(s\right)))ds+\\ \hfill +\delta \mathcal{B}\left(t\right){\int }_0^t{e}^{-\delta (t-s)}\mathsf{\Phi }(-d_1(t,s,\mathcal{B}\left(t\right),\mathcal{B}\left(s\right)))ds,\end{array}$$

(15)

where

• $f(t,\mathcal{B}\left(t\right))=K{e}^{-rt}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{K}\right)+(r-\delta -\frac{{\sigma }^2}{2})t}{\sigma \sqrt{t}})-\mathcal{B}\left(t\right)e^{-\delta t}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{K}\right)+(r-\delta +\frac{{\sigma }^2}{2})t}{\sigma \sqrt{t}})$;
• $d_1(t,s,\mathcal{B}\left(t\right),\mathcal{B}\left(s\right))=\frac{log\left(\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(s\right)}\right)+(r-\delta +\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}}$;
• $d_2(t,s,\mathcal{B}\left(t\right),\mathcal{B}\left(s\right))=\frac{log\left(\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(s\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}}.$

Let n be a positive integer. Let us consider the following partition $\Gamma$ of the interval $[0,T]$ into n intervals of equal length $\Delta =T/n$:

$$0=t_0<t_1<t_2<\cdots <t_{n-1}<t_n=T.$$

For each $t_i$, with $i=0,1,2,\cdots ,n$, using the additive property of integral, let us rewrite Equation (15) in this way:

$$\begin{array}{c}\mathcal{B}\left(t_i\right)=K-f(t_i,\mathcal{B}\left(t_i\right))+\sum _{m=1}^i[-rK{\int }_{t_{m-1}}^{t_m}{e}^{-r(t_i-s)}\mathsf{\Phi }(-d_2(t_i,s,\mathcal{B}\left(t_i\right),\mathcal{B}\left(s\right)))ds+\hfill \\ \hfill +\delta \mathcal{B}\left(t\right){\int }_{t_{m-1}}^{t_m}{e}^{-\delta (t_i-s)}\mathsf{\Phi }(-d_1(t_i,s,\mathcal{B}\left(t_i\right),\mathcal{B}\left(s\right)))ds].\end{array}$$

(16)

Proposition 3. 

For the two integrals appearing in Equation (16), the following equalities hold:

$$\begin{array}{c}\hfill {\int }_{t_{m-1}}^{t_m}{e}^{-r(t_i-s)}\mathsf{\Phi }(-d_2(t_i,s,\mathcal{B}\left(t_i\right),\mathcal{B}\left(s\right)))ds=\\ \hfill =\frac{1}{2}\mathsf{\Phi }(-d_2(t_i,t_{m-1},\mathcal{B}\left(t_i\right),\mathcal{B}\left(t_{m-1}\right)){\int }_{t_{m-1}}^{t_m}{e}^{-r(t_i-s)}ds+\\ \hfill +\frac{1}{2}\mathsf{\Phi }(-d_2(t_i,t_m,\mathcal{B}\left(t_i\right),\mathcal{B}\left(t_m\right)){\int }_{t_m-1}^{t_m}{e}^{-r(t_i-s)}ds+{ϵ}_{m,i}^A\end{array}$$

(17)

and

$$\begin{array}{c}\hfill {\int }_{t_{m-1}}^{t_m}{e}^{-\delta (t_i-s)}\mathsf{\Phi }(-d_1(t_i,s,\mathcal{B}\left(t_i\right),\mathcal{B}\left(s\right)))ds=\\ \hfill =\frac{1}{2}\mathsf{\Phi }(-d_1(t_i,t_{m-1},\mathcal{B}\left(t_i\right),\mathcal{B}\left(t_{m-1}\right))){\int }_{t_{m-1}}^{{\zeta }_m}{e}^{-\delta (t_i-s)}ds+\\ \hfill +\frac{1}{2}\mathsf{\Phi }(-d_1(t_i,t_m,\mathcal{B}\left(t_i\right),\mathcal{B}\left(t_m\right))){\int }_{t_{m-1}}^{t_m}{e}^{-\delta (t_i-s)}ds+{ϵ}_{m,i}^B,\end{array}$$

(18)

where ${ϵ}_{m,i}^A$ and ${ϵ}_{m,i}^B$ are the errors as defined in Proposition 1.

Proof. 

The result follows immediately from the continuity of $\mathsf{\Phi }$ and $\mathcal{B}$ and by means of Proposition 1. □

Equation (12) represents a nonstandard Volterra integral equation of the second kind having the form (1). We solve it by using of the method described in Section 2. By applying this method, when $s=t_i$, the indeterminate form

$$\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t_i\right)}{\mathcal{B}\left(t_i\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t_i-t_i)}{\sigma \sqrt{t_i-t_i}})=\mathsf{\Phi }\left(\frac{0}{0}\right)$$

(19)

arises. Observing that

$$\underset{s\to t^{-}}{lim}\mathsf{\Phi }(-\frac{log\left(\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(s\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}})=\frac{1}{2},$$

(20)

we can define the continuous function

$$\mathsf{\Phi }(-\frac{log\left(\frac{u\left(t\right)}{u\left(s\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}})=\left\{\begin{array}{c}\mathsf{\Phi }(-\frac{log\left(\frac{u\left(t\right)}{u\left(s\right)}\right)+(r-\delta -\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}})\iff 0\le s<t\hfill \\ \frac{1}{2}\iff t=s.\hfill \end{array}\right.$$

The same considerations apply to

$$\mathsf{\Phi }(-\frac{log\left(\frac{u\left(t\right)}{u\left(s\right)}\right)+(r-\delta +\frac{{\sigma }^2}{2})(t-s)}{\sigma \sqrt{t-s}}).$$

Using aforementioned notations, from the numerical approximation of the early boundary exercise, the price of the American put may be found by means of the following formula:

$$\begin{array}{c}\hfill P(t,S)=f(t,S)+rK{\int }_0^t{e}^{-r(t-s)}\mathsf{\Phi }(-d_2(t,s,S,\mathcal{B}\left(s\right)))ds+\\ \hfill -\delta S{\int }_0^t{e}^{-\delta (t-s)}\mathsf{\Phi }(-d_1(t,s,S,\mathcal{B}\left(s\right)))ds.\end{array}$$

(21)

In Equation (21), S represents the current price of the underlying asset (see [4]). The authors in [29] showed that the value of American options can be written as the sum of the corresponding European option price and the early exercise premium. Equation (21) lends itself to this interpretation being $f(t,S)$ the price of an European put.

4. Numerical Results

In this section, we present some numerical results showing the accuracy of the approximations computed by means of the pricing model presented in Section 2 and Section 3. We denote by T the time expressed in year, and the interest rate, dividend yield, and volatility are expressed on an annual basis. The results were obtained using MATLAB on a MacBook Pro with a 2.6 GHz Intel Core i7 processor with 16 GB RAM. We considered $n=110$ steps. The required computational time was about 9 s. The nonlinear system of equations of the form (10), arising in the problem discussed in Section 3, was solved applying the Broyden’s method, a quasi-Newtonian method.

The obtained approximated values of $\mathcal{B}\left(t\right)$ are plotted in Figure 1 for several values of $\sigma$. We chose $S=K=100$, $T=2$, $r=0.05$, $\delta =0$ and $\sigma =0.2,0.40,0.60$.

For several values of $\delta$, the obtained values of $\mathcal{B}\left(t\right)$ are, instead, plotted in Figure 2. We chose $S=K=100$, $T=2$, $r=0.05$, $\sigma =0.2$ and $\delta =0.05,0.10,0.15$.

The numerical approximation of $\mathcal{B}\left(t\right)$ gained from our method was applied to (21); consequently, the trapezoidal numerical integration rule was engaged to calculate the price of an American put. The outcome values were compared with the benchmarks ones, stemming from a binomial tree with n = 10,000 time steps.

Table 1. Numerical approximations of American put price: S as in the first column, $T=3$, $\sigma =0.2$, $r=0.08$, $K=100$, $\delta =0.08$, and $n=32$. For each value of S, the resulting approximation errors are displayed in the second rows within brackets. Some data in Table 1 were extracted from [22]; in addition, underlined values represent approximated prices to which the lowest error corresponds.

S	C-R-R	G-J	B-J	B-D	Ju	K-J-K	K-K	N-B-R	D-M-M-V
80	22.2050	22.2079	22.7106	22.1985	22.2084	22.1942	22.1900	22.2048	22.2075
	(-)	(2.9 × 10−3)	(5.1 × 10−1)	(6.5 × 10−3)	(3.4 × 10−3)	(1.1 × 10−2)	(1.5 × 10−2)	(2 × 10−4)	(2.5 × 10−3)
90	16.2071	16.1639	16.5205	16.1986	16.2106	16.1999	16.1960	16.2068	16.2096
	(-)	(4.3 × 10−2)	(3.6 × 10−1)	(5.9 × 10−2)	(7.2 × 10−3)	(1.1 × 10−2)	(3.9 × 10−3)	(1.1 × 10−4)	(2.5 × 10−3)
100	11.7037	11.7053	11.8106	11.6988	11.7066	11.6991	11.6958	11.7037	11.7061
	(-)	(1.6 × 10−3)	(1.1 × 10−1)	(4.9 × 10−3)	(2.9 × 10−3)	(4.9 × 10−3)	(7.9 × 10−3)	(1.0 × 10−5)	(4.4 × 10−3)
110	8.3671	8.3886	8.4072	8.3630	8.3695	8.3638	8.3613	8.3669	8.3689
	(-)	(2.1 × 10−2)	(4.0 × 10−2)	(4.1 × 10−3)	(2.4 × 10−3)	(3.3 × 10−3)	(5.8 × 10−3)	(2.0 × 10−4)	(1.8 × 10−3)
120	5.9299	5.9435	5.9310	5.9261	5.9323	5.9278	5.9258	5.9298	5.9314
	(-)	(1.4 × 10−2)	(1.1 × 10−3)	(3.8 × 10−3)	(2.4 × 10−3)	(2.1 × 10−3)	(4.1 × 10−3)	(1.0 × 10−4)	(1.5 × 10−3)

displays values of the American put achieved through benchmarks methods. All the considered methods are able to produce an approximation of the boundary. In particular, the Cox, Ross and Rubenstein model ([30]), namely, C-R-R in the first column of Table 1, with n = 10,000 time steps (C-R-R), is considered. The second column contains the American put prices corresponding to different approximated values of $\mathcal{B}\left(t\right)$ using the four-point extrapolation scheme of [15], namely, G-J. The third column of Table 1 reports the modified two-point method of [8], namely, B-J. The fourth column embeds values stemming from the lower and upper bound approximation of [16], indicated as B-D. The fifth column, namely, Ju, shows the exponential functions method combined with a three-point Richardson extrapolation proposed by [17]. In addition, the iteration method proposed by [12] was considered and it is identified as K-J-K. Values within the seventh column in Table 1 represents the trapezoidal formulas approximations of [18] followed by the Newton–Raphson iteration, namely, K-K. N-B-R refers to the more recently product integration method for the approximation of the early boundary in the American option pricing problem exposed in [22]. Lastly, the last column in Table 1, indicated as D-M-M-V, represents the mean-value theorem approach that we proposed in this work.

In Table 1, we considered S as in first column, $T=3$, $\sigma =0.2$, $r=0.08$, $K=100$ and $\delta =0.08$ and $n=32$. The benchmark values in the first column stem from a binomial tree (C-R-R model) with n = 10,000 time steps. For ease of reading, we underlined the best approximations.

The proposed numerical approach employs about $0.66$ s to compute $n=32$ points. It is worth underlining that approximations attained with our method are coherent or even better compared to the similar stemming from benchmark methods. The only exception holds for approximations obtained by [22].

Table 2. Comparison between N-B-R and D-M-M-V: S as in the first column, $T=3$, $\sigma =0.2$, $r=0.08$, $K=100$ and $\delta =0.08$. The benchmark value was obtained considering the C-R-R model with n = 10,000 time steps. For each value of S, the resulting approximation errors are displayed in the second rows within brackets.

S	C-R-R	N-B-R	D-M-M-V
80	22.2050	22.2048	22.2051
	(-)	(2.0 × 10−4)	(1.0 × 10−4)
90	16.2071	16.2068	16.2072
	(-)	(1.1 × 10−4)	(1.0 × 10−4)
100	11.7037	11.7037	11.7040
	(-)	(1.0 × 10−5)	(3.0 × 10−4)
110	8.3671	8.3669	8.3672
	(-)	(2.0 × 10−4)	(1.0 × 10−4)
120	5.9299	5.9299	5.9299
	(-)	(1.0 × 10−4)	(1.0 × 10−4)

shows results from our proposal compared to the ones [22] considering, for D-M-M-V method, $n=110$ points.

As Table 2 shows, the approximations obtained with our method were coherent with the ones obtained by applying the method in [22]. The resulting errors in the second rows for each value of S were very small and comparable with that relative to the N-B-R method. In particular, [22] numerically solved the one-dimensional reformulation of Kim’s integral equations ([11]) using a modified version of the Nyström method. Such a procedure allows for taking into account the singularity close to expiry presented by the early exercise boundary. In our method, instead, we directly dealt with Equation (15). In addition, our method may be easily extended to integral equations coming from more complicated dynamics.

In [22], taking into account $n=32$ points, an error of range of ${10}^{-4}$ was obtained in about 14 s. Their calculations were performed on a PC with a 4.00 Intel Core i7 GHz processor with 16 GB RAM (see Figures 3–6 as reported by [22]). With our algorithm, also considering $n=110$ points, we obtained an error of the same range in about 9 s. Moreover, our calculations were performed using a machine with a less powerful processor, a MacBook Pro with a 2.6 GHz Intel Core i7 processor with 16 GB of RAM. To test the convergence of our algorithm, we report in

Table 3. Convergence analysis in approximating the American put price. The benchmark value, equal to $7.5009$, stems from a binomial tree (C-R-R model) with n = 10,000 steps. We considered $S=K=100$, $T=1$, $r=\delta =0.08$, $\sigma =0.2$.

n	D-D-M-V	Error
10	7.5059	0.0050
20	7.5028	0.0019
30	7.5020	0.0011
40	7.5017	8.0 × 10−4
90	7.5012	3.0 × 10−4
110	7.5012	1.0 × 10−4

some values of the American put corresponding to as many as the approximated values of $\mathcal{B}\left(t\right)$. We considered: $S=K=100$, $T=1$, $r=0.08$, $\sigma =0.2$ and $\delta =0.08$. In Column 1, we report the time steps. The errors were computed as the difference between our values (D-D-M-V) with the ones obtained considering a binomial tree (C-R-R) model with n = 10,000 used as a benchmark. The error size of order ${10}^{-4}$ obtained when increasing the value of n proves the quality and the precision of the proposed method.

5. Conclusions

The present work proposed a new numerical method for the approximation of the early exercise boundary of an American option. American options are financial contracts allowing for the holder to buy or sell an underlying financial asset before the contract expiration. Thus, the evaluation of an American option requires to consider that the holder may have convenience to interrupt the contract at any time depending on the value of the underlying asset; that is, a collection of optimal exercise prices, one for each optimal exercise time, must be determined. Representing such a boundary by means of a nonstandard Volterra integral equation, we propose a numerical scheme in order to approximate the solution of the free boundary problem concerning the American option. Our proposal allows for managing well-known numerical problems arising in the presence of a singularity close to expiry and of a double integral. First, we provided theoretical instruments that were at the core of the proposed method. Consequently, we compared existing pricing models with results stemming from our method, showing the accuracy of the latter. The direct solution of the integral equation suggests that our method can be applied to more complicated cases coming up in finance.

Future research works will concern the use of an asset price process characterized by a nonconstant volatility term. In particular, we aim to expand our proposal to the case of the Heston model for American options, considering the associated multidimensional Volterra integral equation for the boundary problem.