A simple and efficient numerical method for pricing discretely monitored early-exercise options

TL;DR

This paper introduces a simple, fast, and accurate quadrature-based numerical method for pricing discretely monitored early-exercise options within the Black-Scholes model, suitable for various structured products and options.

Contribution

It proposes a novel quadrature technique that uses elementary calculations and a fixed grid, achieving high convergence rate and efficiency for pricing discretely monitored options.

Findings

Convergence rate of $O(1/N^4)$

Computational complexity of $O(MN ext{log} N)$

Applicable to a wide range of options including barrier and Bermudan options

Abstract

We present a simple, fast, and accurate method for pricing a variety of discretely monitored options in the Black-Scholes framework, including autocallable structured products, single and double barrier options, and Bermudan options. The method is based on a quadrature technique, and it employs only elementary calculations and a fixed one-dimensional uniform grid. The convergence rate is $O(1/N^4)$ and the complexity is $O(MN\log N)$, where $N$ is the number of grid points and $M$ is the number of observation dates.

Tables

Table 1. Table 1. An autocallable structured product.

Observation date	Barrier level	Risk-free rate
0.2	3050	2
0.4	3100	2.1
0.6	3150	2.2
0.8	3200	2.3
1	3250	2.4

Table 2. Table 2. A double barrier option.

Observation date	Barrier level 1	Barrier level 2	Risk-free rate
0.25	2200	2800	1
0.50	2100	2900	1.1
0.75	2000	3000	1.2
1	1900	3100	1.3
1.25	1800	3200	1.2
1.50	1700	3300	1.3
1.75	1600	3400	1.4
2	$-$	$-$	1.5

Equations

$$V(t_m,S)$$

$$=e^{-r(t_{m+1}-t_m)}{\int }_0^{\infty }V(t_{m+1},y)f(y|S)dy,$$

$$dS(t)=\bigl{[}r(t)-q(t)\bigr{]}S(t)\,dt+\sigma(t)S(t)\,dW(t),$$

$$S(t)=S(t_{0})\exp\biggl{\{}\int_{t_{0}}^{t}\bigl{[}r(s)-q(s)-\tfrac{1}{2}\sigma^{2}(s)\bigr{]}\,ds+\int_{t_{0}}^{t}\sigma(s)\,dW(s)\biggr{\}},$$

$$S(t_{m})\sim S(t_{0})\exp\biggl{\{}\int_{t_{0}}^{t_{m}}\bigl{[}r(s)-q(s)-\tfrac{1}{2}\sigma^{2}(s)\bigr{]}\,ds+\Bigl{(}\int_{t_{0}}^{t_{m}}\sigma^{2}(s)\,ds\Bigr{)}^{1/2}\,Z\biggr{\}},$$

$$r_m={\int }_{t_{m-1}}^{t_m}\frac{r(s)}{\Delta t_m}ds,q_m={\int }_{t_{m-1}}^{t_m}\frac{q(s)}{\Delta t_m}ds,{\sigma }_m^2={\int }_{t_{m-1}}^{t_m}\frac{{\sigma }^2(s)}{\Delta t_m}ds,$$

$$\tilde{r}(t)=r_m,\tilde{q}(t)=q_m,\text{and}\tilde{\sigma }(t)={\sigma }_m,\text{for}{t}_{m-1}<t\le t_m.$$

$$d\tilde{S}(t)=\bigl{[}\tilde{r}(t)-\tilde{q}(t)\bigr{]}\tilde{S}(t)\,dt+\tilde{\sigma}(t)\tilde{S}(t)\,dW(t),$$

$${\int }_{t_0}^{t_m}r(s)ds$$

$${\int }_{t_0}^{t_m}q(s)ds$$

$${\int }_{t_0}^{t_m}{\sigma }^2(s)ds$$

$$r(t)=r_m,q(t)=q_m,\text{and}\sigma (t)={\sigma }_m,\text{for}{t}_{m-1}<t\le t_m.$$

$$V(t_M,S)=a_{M}S+b_M,\text{for}{K}_M^{-}<S<K_M^{+}.$$

$$1\le m\le M:$$

$$m=M:$$

$$1\le m\le M-1:$$

$$m=M:$$

$$a_M^{+}=b_M^{+}=0.$$

$$1\le m\le M-1:$$

$$m=M:$$

$$a_M^{-}=b_M^{-}=0.$$

$$1\le m\le M:$$

$$K_m^{+}=\infty ,a_m^{-}=-1,b_m^{-}=K;$$

$$m=M:$$

$$V_m(S)=V(t_m,S),m=0,1,\dots ,M,$$

$${\tau }_m=\frac{1}{2}{\sigma }_m^2\Delta t_m=\frac{1}{2}{\sigma }_m^2(t_m-t_{m-1}).$$

$$\rho_{m}(y,S)=\frac{1}{2\sqrt{\pi\tau_{m}}\,y}\exp\biggl{\{}-\frac{1}{4\tau_{m}}\Bigl{(}\log\frac{y}{S}-\frac{2}{\sigma_{m}^{2}}\bigl{[}r_{m}-q_{m}-\tfrac{1}{2}\sigma_{m}^{2}\bigr{]}\tau_{m}\Bigr{)}^{2}\biggr{\}}.$$

$$\displaystyle V_{m-1}(S)=e^{-2r_{m}\tau_{m}/\sigma_{m}^{2}}\mathbb{E}\bigl{[}V_{m}(\cdot)|S\bigr{]}=e^{-2r_{m}\tau_{m}/\sigma_{m}^{2}}\int_{0}^{\infty}V_{m}(y)\rho_{m}(y,S)\,dy$$

$$\displaystyle=\frac{e^{-2r_{m}\tau_{m}/\sigma_{m}^{2}}}{2\sqrt{\pi\tau_{m}}}\int_{0}^{\infty}\frac{1}{y}\,V_{m}(y)\exp\biggl{\{}-\frac{1}{4\tau_{m}}\Bigl{(}\log\frac{y}{S}-\frac{2}{\sigma_{m}^{2}}\bigl{[}r_{m}-q_{m}-\tfrac{1}{2}\sigma_{m}^{2}\bigr{]}\tau_{m}\Bigr{)}^{2}\biggr{\}}\,dy,$$

$$V_{m-1}(S)=\{\begin{array}{ll}{a}_{m-1}^{-}S+b_{m-1}^{-},& S\le K_{m-1}^{-}\\ a_{m-1}^{+}S+b_{m-1}^{+},& S\ge K_{m-1}^{+}\end{array}.$$

$$V_m^a(S,K,ϵ)=e^{-2q_m{\tau }_m/{\sigma }_m^2}SN(ϵd_1),V_m^b(S,K,ϵ)=e^{-2r_m{\tau }_m/{\sigma }_m^2}N(ϵd_2),$$

$$d_{1}=\frac{1}{\sqrt{2\tau_{m}}}\Bigl{(}\log\frac{S}{K}+\frac{2}{\sigma_{m}^{2}}\bigl{[}r_{m}-q_{m}+\tfrac{1}{2}\sigma_{m}^{2}\bigr{]}\tau_{m}\Bigr{)},\qquad d_{2}=d_{1}-\sqrt{2\tau_{m}}.$$

$$V_{\mathrm{smile}}=V_{\mathrm{no}\mathrm{smile}}-\frac{\partial V_{\mathrm{vanilla}}}{\partial \sigma }\frac{\partial \sigma }{\partial K}.$$

$$V_M(S)=\{\begin{array}{ll}{a}_M^{-}S+b_M^{-},& S\le K_M^{-}\\ a_{M}S+b_M,& K_M^{-}<S<K_M^{+}\\ a_M^{+}S+b_M^{+},& S\ge K_M^{+}\end{array}.$$

$${\tilde{V}}_{M-1}(S)$$

$$\displaystyle\qquad{}+a_{M}\bigl{[}V_{M}^{a}(S,K_{M}^{-},1)-V_{M}^{a}(S,K_{M}^{+},1)\bigr{]}$$

Full text

Disclaimer: The statements made and opinions expressed here are solely our own, and do not reflect the views of China Mechants Bank or its employees and affiliates.

1. Introduction

Exotic options are commonly traded throughout the world. Many popular exotic options are path-dependent and have early-exercise features. These options can often be priced using analytical formulas if they are continuously monitored (e.g. barrier options). In practice, however, most path-dependent exotic options are discretely monitored [4], in which case they need to be priced using numerical techniques. Due to the complicated structures of these options, traditional pricing models based on Monte-Carlo simulations and finite difference methods are often too time-consuming to be useful in practical situations. More recent pricing methods based on advanced mathematical techniques, on the other hand, tend to be more efficient (e.g. [14, 12, 10, 11]), but for many financial institutions, these methods are often too difficult to understand and to properly implement. To strike a balance between model performance and practical utility, we propose a new quadrature-based method that is much faster and more accurate than the traditional Monte-Carlo and PDE methods, yet at the same time is easy to understand and to implement. We will first give a brief review of the types of products considered, as well as the quadrature-based pricing model which is the foundation of our work. Then we will explain our method and provide numerical examples.

1.1. Autocallable Structured Products

Autocallable structured products belong to a class of exotic options with early-exercise features. Many different types of autocallable products have been created and traded in financial markets, and they have become increasingly popular in recent years. We refer to the Appendix of [9] for a description of the main features of various autocallable products.

We will consider a very common autocallable product with discrete observation dates. At each observation date there is a pre-specified barrier level. If the price of the underlying asset is greater (less) than or equal to the barrier level (depending on the terms of the product), the option is exercised and a pre-specified fixed-rate return is paid. If the asset price is below (above) the barriers at all observation dates, the option is never exercised and the investor receives a negative return at maturity. In addition, autocallable products may have a knock-in feature. In this case, if the option is never exercised, the negative return the investor receives depends on whether the asset price at maturity reaches a pre-specified knock-in level. While the value of a continuously monitored autocallable product has a simple closed-form solution, the value of a discretely monitored autocallable product cannot be calculated easily. In the discrete case, there exist analytical solutions in terms of multiple integrals, cf. [25, 26, 9]. The numerical calculation of these integrals, however, can become prohibitive if the number of observation dates exceeds five. In practice, discretely monitored autocallable products are commonly priced using Monte-Carlo simulations. This method is straightforward, but convergence is usually slow and acceleration techniques such as variance reduction are often needed (cf. [2, 13]). Another popular method for pricing discretely monitored autocallable products is to solve the governing Black-Scholes partial differential equation (PDE) using finite difference method (cf. [9]). Assuming a second-order central difference approximation in space, the overall convergence rate of a typical finite difference based pricing method is $O((\delta x)^{2}+\delta t)$ if the explicit forward Euler method is used in time, and $O((\delta x)^{2}+(\delta t)^{2})$ if the implicit Crank-Nicolson method is used. Since two-dimensional grids are needed for finite difference methods, computational complexities are at least of order $1/(\delta x\delta t)$. In addition, since the payoffs of autocallable products are discontinuous (in asset price), additional care (such as smoothing of payoff functions) must be taken to ensure the accuracy of any finite difference approximations.

Remark 1.1**.**

Some autocallable products may have payoffs at maturity that are of the same type as that of European vanilla options. These products can be effectively viewed as combinations of autocallable products and barrier options (see below).

1.2. Discrete Barrier Options

Barrier options are among the most popular types of exotic options. A barrier option may be activated (knock-in option) or deactivated (knock-out option) when the price of the underlying asset crosses certain barrier levels. Barrier options may be discretely or continuously monitored. A single barrier option has one barrier at each observation date, while a double barrier option has two barriers at each observation date. The final payoff of a barrier option (if it is active at maturity) may be of the same type as that of a vanilla option or that of a digital option. A special type of barrier option has a constant amount of cash as the final payoff, and such an option is called a touch option. Most barrier options have time-independent barrier levels, but options with time-dependent barrier levels have also been studied [7].

Similar to the case of autocallable structured products, there are closed-form solutions for discrete barrier options in terms of multiple integrals [23], but such solutions are often difficult to evaluate directly. In practice, discrete barrier options can be priced using Monte-Carlo simulations or standard binomial tree methods, but these methods are usually slow [5]. Other methods that have been proposed to price discrete barrier options include continuity correction approximations [4, 27], Wiener-Hopf methods [14], adaptive mesh methods [1], Hilbert transform methods [12], finite element methods [16], Fourier-cosine series expansion methods [10], and quadrature methods [3, 6]. These methods, while useful in certain contexts, have not been as widely used as the traditional Monte-Carlo and finite difference methods, usually due to their complexity.

1.3. Bermudan Options

Bermudan options are discrete versions of American options. A Bermudan option can be exercised at any of the prescribed observation dates, and the payoff is of the same type as that of a vanilla option. Similar to discrete barrier options, Bermudan options can be priced using Monte-Carlo simulations [17], Hilbert transform methods [11], Fourier-cosine series expansion methods [10], and quadrature methods [22, 21, 6].

1.4. Overview of the Quadrature Method

Among the various methods proposed to price discretely monitored options, the quadrature method is particularly appealing because of its high efficiency and accuracy. The method has been applied to discrete barrier options and Bermudan options [3, 22, 21, 6]. The main idea is to solve for option values at each observation date via backward induction in time. The risk-neutral valuation formula is expressed as a single integral, which is then evaluated numerically to produce the option price. Specifically, let $V$ denote the value of the option, $S$ the value of the underlying asset, $r$ the risk-neutral interest rate, $t_{1},\dotsc,t_{M}$ the observation dates, and $\mathbb{E}$ the risk-neutral expectation. If the underlying asset $S$ does not trigger an early exercise at $t_{m}$, we have

[TABLE]

where $f(y|S)$ is a probability density function whose form depends on the model of the underlying asset. If $S$ triggers an early exercise at $t_{m}$, on the other hand, the option price $V(t_{m},S)$ would be equal to a prescribed value. The integral above can be calculated using FFT [21, 22] or Fast Gauss Transform [6]. Since, however, $V(t_{m+1},y)$ is discontinuous (in $y$) for autocallable products and barrier options, and non-differentiable (in $y$) for Bermudan options, care must be taken to ensure the accuracy of the numerical evaluation of the integral. While several shifted or nonuniform grids have been designed in previous studies to address this difficulty [21, 6], the problem becomes particularly challenging when multiple discontinuities are present at each observation point, for instance in the case of double barrier options with time-dependent barrier levels.

1.5. Motivation of Our Work

Although a good number of advanced techniques have been proposed to improve pricing models’ accuracy and efficiency, in most practical situations, simple methods that are more cost effective are usually preferred over their more sophisticated counterparts. The reason for this is twofold. First of all, when data quality is not high enough, sophisticated models are not necessarily beneficial. For instance, interest rates and volatilities are crucial components of nearly every pricing model, but they need to be estimated from available market data. If the estimated parameters contain large errors, which is not uncommon in products sold in emerging markets, any advantages gained from the use of sophisticated models may be (more than) offset by these errors, making the simpler models more attractive. Secondly, implementations of pricing models usually involve staff members from multiple business departments, and the resulting products often need active maintenance and updates. As a result, models that are too complicated in nature may hinder effective business communications, which increases maintenance costs and operational risks. In view of these considerations, it is not difficult to see why traditional Monte-Carlo and PDE methods are still among the most popular methods in the valuation of discretely monitored options, even though their computational costs are already high enough to adversely impact their applicability in business.

In view of the practical concerns mentioned above, we propose a new quadrature method to price the aforementioned discretely monitored options in the Black-Scholes framework. The convergence rate of our method is $O(1/N^{4})$ and the complexity is $O(MN\log N)$, where $N$ is the number of grid points and $M$ is the number of observation dates. The performance of our method is on par with previous quadrature-based methods such as the CONV method [21], but it is more straightforward, and is better suited for products with multiple discontinuities. Our method differs from other quadrature methods mainly in three aspects. First, we work with probability density functions directly instead of using characteristic functions or Toeplitz matrices. Secondly, we use only a fixed one-dimensional uniform grid to compute all integrals. Thirdly, we utilize explicit Black-Scholes formulas to improve the accuracy of the calculations. Due to these novel modifications, our method is very easy to implement, and is capable of handling sophisticated products such as double-barrier options with time-dependent barrier levels.

1.6. Organization of the Paper

The rest of the paper is organized as follows. Section 2 specifies the class of (discrete) option pricing problems that our quadrature method is intended to solve, and Section 3 presents the main recursion formula for our method. After detailing the implementation of our method in Section 4–5, we summarize the algorithm in Section 6 and then present numerical examples in Section 7. The Appendix collects a few useful theoretical results, which lay the foundation of a class of (discrete) option pricing algorithms (including the one described here) but which, to our best knowledge, do not seem to have been rigorously proved or even properly formulated in the literature. We supply the proofs here in the hope that they would be useful to interested readers.

2. Basic Assumptions

We assume that the price of the underlying asset $S(t)$ satisfies the following stochastic differential equation in the risk-neutral measure:

[TABLE]

where $r(t)$ is the risk-neutral interest rate, $q(t)$ is the yield rate, $\sigma(t)$ is the volatility, and $W(t)$ is the Wiener process.

In practice, interest rates are always time dependent. Yield rates for FX products are simply foreign interest rates, and for other types of products they may be implied from futures prices. Thus yield rates are usually time dependent as well. Implied volatilities are time dependent, whereas historical volatilities can often be taken as constant.

The solution to (2.1) is

[TABLE]

where $t_{0}$ is the present date. We consider a discretely monitored option with observation dates $t_{1},\dotsc,t_{M}$, where the last observation date $t_{M}$ is the maturity date. It follows from (2.2) that each $S(t_{m})$ for $1\leq m\leq M$ has the lognormal distribution

[TABLE]

where $Z$ denotes the standard normal distribution. Now define

[TABLE]

for $1\leq m\leq M$ where $\Delta t_{m}=t_{m}-t_{m-1}$, and define piecewise constant functions

[TABLE]

For the process

[TABLE]

since

[TABLE]

it follows that $\tilde{S}(t_{m})$ has the same distribution as $S(t_{m})$ in (2.3) for each $m$. Since the value of the option depends only on probability distributions of the asset price at observation dates, the option value remains the same if we replace the process $S$ by the process $\tilde{S}$. In other words, we may safely assume that $r(t),\ q(t)$, and $\sigma(t)$ are piecewise constant functions. Thus in what follows, we shall assume

[TABLE]

Consider now a general class of discretely monitored options with barriers. Since the sum of a knock-in barrier option and a knock-out barrier option with the same observation dates and barrier levels is a vanilla option (or a digital option if the barrier options are digital), to study the pricing of these discretely monitored options, it suffices to consider a knock-out barrier option which ceases to exist when barrier levels are crossed. To this end, assume that

1. (A)

The option has two strike prices $K_{m}^{-},\,K_{m}^{+}\in[0,\infty]$, with $K_{m}^{-}\leq K_{m}^{+}$, at each observation date $t_{m},\ m=1,2,\dotsc,M$. 2. (B)

The option is exercised if $S\leq K_{m}^{-}$ or $S\geq K_{m}^{+}$ at some $t_{m}$, and the payoffs are given by $a_{m}^{-}S+b_{m}^{-}$ (if $S\leq K_{m}^{-}$) and $a_{m}^{+}S+b_{m}^{+}$ (if $S\geq K_{m}^{+}$), respectively, for some $a_{m}^{\pm},\,b_{m}^{\pm}\in\mathbb{R}$. 3. (C)

The final payoff at maturity is

[TABLE]

These assumptions are general enough to cover a wide class of discretely monitored options, such as the ones mentioned in the introduction. For instance, common up-and-out autocallable products would have

[TABLE]

Down-and-out put barrier options would have

[TABLE]

Double barrier knock-out call options would have

[TABLE]

Bermudan put options with strike $K$ would have

[TABLE]

We will give a proof of the uniqueness of $K_{m}^{\pm}$ for Bermudan options in the Appendix.

To summarize, our basic assumptions are

1. (1)

The underlying asset price $S$ follows a geometric Brownian motion with piecewise constant interest rates, yield rates, and volatilities. 2. (2)

There are finitely many observation points, and two exercise levels (possibly $\infty$) at each observation point. If $S$ is above the upper exercise level or below the lower exercise level at any observation point, the option is exercised and the payoff is a linear function in $S$. 3. (3)

At maturity, if $S$ is between the two exercise levels, a payoff is incurred which is also a linear function in $S$.

3. Outline of the method

Let $V(t,S)$ denote the value of the option (as a function of asset price $S$) at any time $t$, and let

[TABLE]

denote the value of the option at the observation dates. Our goal is to find $V_{0}(S(t_{0}))$, and our strategy is to use backward induction in time. Since $V_{M}(S)$ is piecewise linear in $S$, $V_{M-1}(S)$ has a simple explicit expression. For each $m=M-1,\dotsc,1$, we write $V_{m-1}(S)$ as the risk-neutral expectation of $V_{m}(S)$ for $K_{m-1}^{-}<S<K_{m-1}^{+}$, and as $a_{m-1}^{\pm}S+b_{m-1}^{\pm}$ otherwise. The expectation is given by an explicit integral and is calculated numerically. The core of the quadrature method is the calculation of $M-1$ expectation integrals step-by-step. Let

[TABLE]

For each $1\leq m\leq M-1$, note that $S(t_{m})$ has a lognormal distribution as in (2.3). The relevant probability density functions are known to be [19]

[TABLE]

For simplicity of notations we define $K_{0}^{+}=\infty$ and $K_{0}^{-}=0$. By the fundamental theorems of asset pricing, we have the risk-neutral pricing formula [24]:

[TABLE]

for $K_{m-1}^{-}<S<K_{m-1}^{+}$ and $1\leq m\leq M-1$. By Assumption (B) from Section 2, we also have

[TABLE]

To further study the formulas (3.2) and (3.3), we first recall some classical results on the pricing of binary options.

Lemma 3.1**.**

Let $K>0$, and let $\chi_{A}$ denote the characteristic function of a set $A$. Consider an option with no early-exercise features.

1. (1)

If the option has payoff $\hat{V}_{m}(y)=\chi_{[K,\infty)}y$, then $\hat{V}_{m-1}(S)=V_{m}^{a}(S,K,1)$. 2. (2)

If $\hat{V}_{m}(y)=\chi_{(0,K]}y$, then $\hat{V}_{m-1}(S)=V_{m}^{a}(S,K,-1)$. 3. (3)

If $\hat{V}_{m}(y)=\chi_{[K,\infty)}$, then $\hat{V}_{m-1}(S)=V_{m}^{b}(S,K,1)$. 4. (4)

If $\hat{V}_{m}(y)=\chi_{(0,K]}$, then $\hat{V}_{m-1}(S)=V_{m}^{b}(S,K,-1)$.

The functions $V_{m}^{a}$ and $V_{m}^{b}$ are defined as

[TABLE]

where $N$ is the cumulative normal distribution function, and

[TABLE]

Proof.

By definition $V_{m}^{a}$ is the value of an asset-or-nothing option, and $V_{m}^{b}$ is the value of a cash-or-nothing option. The valuation formulas are just standard results for binary options [18]. ∎

Remark 3.2**.**

The standard Black-Scholes formulas in Lemma 3.1 ignore the possible effects of volatility smiles. If such effects need to be taken into account, one may amend the definitions of $V_{m}^{a}$ and $V_{m}^{b}$ (as given in the Lemma) by incorporating suitable vega-induced correction terms [15]. For instance, the value of a cash-or-nothing call option in the presence of volatility smiles would become

[TABLE]

We can use Lemma 3.1 to obtain an explicit formula for the value of $V_{M-1}(S)$. By Assumption (B)–(C) from Section 2, we have

[TABLE]

Without loss of generality we may assume $0<K_{M}^{\pm}<\infty$, since otherwise we may choose some arbitrary $0<K_{M}^{\pm}<\infty$ and set $a_{M}^{\pm}=a_{M},\ b_{M}^{\pm}=b_{M}$.

Proposition 3.3**.**

The value of the option at $t_{M-1}$ is given by

[TABLE]

for $K_{M-1}^{-}<S<K_{M-1}^{+}$.

Proof.

Clearly, the option from $t_{M-1}$ to $t_{M}$ is equivalent to a linear combination of binary options consisting of two put options with strike $K_{M}^{-}$, two call options with strike $K_{M}^{-}$ and four call options with strike $K_{M}^{+}$. The result then follows from Lemma 3.1. ∎

With the aid of (3.2) and Proposition 3.3, we may write the main recursion of our quadrature method as follows.

Proposition 3.4**.**

Let $\tilde{V}_{M}=V_{M}$ be defined in (3.4), and $\tilde{V}_{m-1}$ be given by the following recursion formula:

[TABLE]

for $1\leq m\leq M$. Then we have $\tilde{V}_{m}(S)=V_{m}(S)$ for $K_{m}^{-}<S<K_{m}^{+}$ and $0\leq m\leq M$. In particular, $\tilde{V}_{0}(S(t_{0}))=V_{0}(S(t_{0}))$.

Proof.

This is merely the classical recursion formula for quadrature methods specialized to the Black-Scholes model. To prove the formula, we only need to show

[TABLE]

for all $0\leq m\leq M$. By assumption (3.6) is true for $m=M$. Assume now (3.6) holds for some $1\leq m\leq M$. Substituting the equation into (3.2), applying Lemma 3.1, comparing the result with (3.5) and using (3.3), we observe that (3.6) holds for $m-1$. The result then follows from induction. ∎

Remark 3.5**.**

Recursion formula (3.5) lies at the heart of our quadrature method and distinguishes our method from other quadrature methods, which are primarily based on (3.2) or one of its many variants. The significance of the formula (3.5) lies in the fact that it makes explicit use of Black-Scholes formulas to separate the expectation integral \mathbb{E}\bigl{[}V_{m}(\cdot)|S\bigr{]} into a “quadrature part”

[TABLE]

and an “early-exercise part”

[TABLE]

Since the function $\tilde{V}_{m}(S)$ is smooth for $S\in(K_{m}^{-},K_{m}^{+})$ (in fact for all $S\in(0,\infty)$, as we will show below), the integral $F_{m-1}(S)$ can be evaluated accurately and efficiently using a high-order quadrature method such as Simpson’s rule. In contrast, the integrand $V_{m}(S)$ in the original recursion formula (3.2) is discontinuous on $(0,\infty)$ (in either $V_{m}$ itself or in its first derivative); this makes the accurate evaluation of the expectation integral a difficult and challenging task. Although (3.5) applies specifically to the Black-Scholes model, the same idea can be used for other asset price models, as long as a suitable analytical formula (exact or approximate) can be found for the probability density function $\rho_{m}(y,S)$ and early-exercise part $E_{m-1}(S)$.

With $\tilde{V}_{M-1}(S)$ given in Proposition 3.3, Proposition 3.4 implies that we may apply (3.5) successively to obtain $\tilde{V}_{M-2}(S),\tilde{V}_{M-3}(S),\dotsc,\tilde{V}_{0}(S)$. The value of the option is equal to $\tilde{V}_{0}(S(t_{0}))$.

4. Details of Implementation

In (3.5), $\tilde{V}_{m-1}(S)$ is written as a sum of explicit functions and an integral, namely

[TABLE]

where

[TABLE]

We may truncate the integral by replacing its upper and lower bounds by

[TABLE]

respectively, where $C>1$ is a suitable constant. In practice, the choice

[TABLE]

where $\sigma_{0}=\max_{1\leq m\leq M}\sigma_{m}$, is sufficient to reduce the truncation errors to round-off level. Heuristically this is clear from (2.3), which suggests that the chance that $S(t_{m})$ move outside the range $(S(t_{0})/C,S(t_{0})C)$ is negligibly small. The rigorous derivation of the error bounds can be obtained using a recursive argument, as will be explained in the Appendix.

Now we consider the truncated integral

[TABLE]

If $K_{m}^{-}\geq S(t_{0})C$ or $K_{m}^{+}\leq S(t_{0})/C$, then by convention the integral is zero. Thus in what follows we shall assume $K_{m}^{-}<S(t_{0})C$ and $K_{m}^{+}>S(t_{0})/C$. Let

[TABLE]

and denote

[TABLE]

The truncated integral can be rewritten as

[TABLE]

One can show by differentiating (4.2) that $\tilde{F}_{m}$, and thus $\tilde{V}_{m}$ and $u_{m}$, are smooth functions in $x$. This means we can compute the integrals efficiently using a high-order quadrature such as Simpson’s rule.

In general, $B_{m}^{\pm}$ are different for different values of $m$, so they cannot all be placed on one grid. Now we choose a uniform grid $\mathbf{x}=\{x_{1},x_{2},\dotsc,x_{N}\}$, where $x_{1}=-\log C$ and $x_{N}=\log C$. Let

[TABLE]

For each $m$, let

[TABLE]

where by definition $p_{m}^{-}\geq 1$ and $p_{m}^{+}<N$. Since we will use Simpson’s rule which requires an odd number of grid points, we define

[TABLE]

and rewrite (4.2) as

[TABLE]

For each $2\leq m\leq M-1$, we will compute $\tilde{F}_{m-1}(S(t_{0})e^{x})$ for all

[TABLE]

where

[TABLE]

For $m=1$ we only need to compute $\tilde{F}_{m-1}(S(t_{0})e^{x})$ for $x=0$, since the value of the option is given by $\tilde{V}_{0}(S(t_{0}))$.

4.1. Computation of the first integral in (4.3)

To compute the first integral in (4.3) using Simpson’s rule, we let

[TABLE]

The integral is discretized as

[TABLE]

since Simpson’s rule is of order 4 [8]. Note that $U_{m}(i)$ is known from the previous step (or by Proposition 3.3 for $m=M-1$) for all $i=1,2,\dotsc,N$. For $x\in\{B_{m-1}^{\pm},\xi_{m-1}^{\pm},0\}$, the sum (4.4) can be computed directly with complexity $O(N)$. For all grid points $x\in\{x_{1},x_{2},\dotsc,x_{N}\}$, on the other hand, the sum (4.4) can be computed altogether using FFT with complexity $O(N\log N)$. This latter fact is crucial to the efficient implementation of our quadrature method and is a consequence of the following simple observation.

Proposition 4.1**.**

Define $(2N-1)$-periodic grid functions $\hat{z},\ \hat{U}_{m}$, and $\hat{F}_{m}$ by

[TABLE]

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the discrete Fourier transform and the inverse discrete Fourier transform of size $2N-1$, respectively. Then

[TABLE]

for all $1\leq j\leq N$. The above discrete Fourier transforms and inverse discrete Fourier transform can be calculated using FFT, and the total computational complexity is $O(N\log N)$.

Proof.

We consider the discrete convolution

[TABLE]

for $j\in\mathbb{Z}$. Note that by definition,

[TABLE]

for all $1-N\leq j-i\leq N-1$. Thus for $1\leq j\leq N$, we have

[TABLE]

Therefore (4.4) with $x=x_{j}$ can be written as

[TABLE]

The discrete convolution $G_{m}$ can be calculated using FFT as

[TABLE]

with a complexity of $O(N\log N)$ [8]. ∎

Remark 4.2**.**

Our method differs from other well-known FFT-based methods (such as [21, 22]) in that we express the discrete quadrature rule (4.4) directly in terms of discrete Fourier transforms, instead of applying continuous Fourier transform to the integral and then discretizing the Fourier integrals (in other words, we have exchanged the order of Fourier transform and discretization). The direct application of the discrete Fourier transform (to the discrete quadrature rule) not only eliminates the need for artificially-introduced damping factors, which are required for the existence of the continuous Fourier transforms, but also eliminates the need for additional specially-designed computational grids which are required to satisfy Nyquist relations. This enables us to carry out the main recursion (3.5) on a fixed uniform grid, without any additional artificial parameters.

4.2. Computation of the last two integrals in (4.3)

The last two integrals in (4.3) are calculated in similar ways using Simpson’s rule. First, note that we may use Proposition 3.3 to calculate $\tilde{V}_{M-1}(S(t_{0})e^{x})$ for $x\in\{B_{M-1}^{\pm},\xi_{M-1}^{\pm}\}$. Generally all four points are needed if the option has two barriers, and only two are needed if the option has one barrier. For each $2\leq m\leq M-1$, assume $u_{m}(x)=\tilde{V}_{m}(S(t_{0})e^{x})$ has been calculated for $x\in\{B_{m}^{\pm},\xi_{m}^{\pm}\}$. The last two integrals in (4.3) are calculated using Simpson’s rule as follows:

[TABLE]

5. Finding Optimal Exercise Prices for Bermudan options

Unlike autocallable products and barrier options, Bermudan options do not have pre-specified exercise levels. Instead, one needs to solve for $K_{m}^{\pm}$ from the equations

[TABLE]

for call options and

[TABLE]

for put options, where $\tilde{V}_{m}$ is determined by (3.5). For simplicity we assume the yield rates $q_{m}\geq 0$, which is almost always the case in practice. We will demonstrate how to find $K_{m}^{-}$, as the same procedure applies to $K_{m}^{+}$. Let

[TABLE]

If $p=1$ there is no early exercise, so $K_{m}^{-}=0$. Otherwise we have $S(t_{0})e^{x_{p-1}}\leq K_{m}^{-}<S(t_{0})e^{x_{p}}$ by Corollary 8.3. The value of $K_{m}^{-}$ can be found using classical root-finding methods such as the bisecting method or the secant method. Note that the bisecting method is guaranteed to converge by Corollary 8.3, and it takes $O(\log N)$ steps to reduce the error of the approximate root to an order of $O(h^{4})$. Since the cost for calculating $\tilde{V}_{m}$ at one point using (4.1) is $O(N)$, the total cost for finding the optimal exercise price is $O(N\log N)$. The secant method is superlinear and converges faster than the bisecting method, though its error estimates are not as straightforward.

6. Summary of the Algorithm

We summarize our algorithm as follows:

1:  Define the functions $V_{m}^{a},V_{m}^{b}$ as in Lemma 3.1

2:  Define the function $\tilde{V}_{M-1}$ as in Proposition 3.3

3:  if option style is Bermudan then

4:     Calculate $K_{M-1}^{\pm}$ as in Section 5

5:  end if

6:  Calculate $p_{M-1}^{\pm}$ and $p_{0}$ to find the bounds of integration in (4.3)

7:  Use Proposition 3.3 to compute $\tilde{V}_{M-1}(S(t_{0})e^{x})$ for $x\in\{B_{M-1}^{\pm},\xi_{M-1}^{\pm}\}$, and assign their values to $v_{1}^{\pm},v_{2}^{\pm}$ respectively

8:  Define a vector $\mathbf{S}$ as $S(i)\leftarrow S(t_{0})e^{x_{i}}$ for $i=1,2,\dotsc,N$

9:  Define a vector $\mathbf{y}$ as $y(i)\leftarrow\tilde{V}_{M-1}(S(i))$ for $i=1,2,\dotsc,N$

10:  for $m=M-1$ downto 2 do

11:     Let $\hat{z}$ and $\hat{U}_{m}$ be as defined in Proposition 4.1

12:     \hat{F}_{m}\leftarrow\mathcal{F}^{-1}\Bigl{\{}\mathcal{F}\bigl{(}w_{m}(\hat{z})\bigr{)}\mathcal{F}(\hat{U}_{m})\Bigr{\}}

13:     Define (or redefine) the vector $\mathbf{Y_{1}}$ as $Y_{1}(j)\leftarrow\frac{h}{3}\hat{F}_{m}(j+N)$ for $j=1,2,\dotsc,N$

14:     Use (4.6), (4.7), and the values of $v_{1}^{\pm},v_{2}^{\pm},y(p_{m}^{-}),y(p_{m}^{+}+p_{0})$ to compute the last two integrals in (4.3) at $x_{i}$ for $i=1,2,\dotsc,N$, and assign their values to $\mathbf{Y_{2}}$ and $\mathbf{Y_{3}}$

15:     if option style is Bermudan then

16:        Calculate $K_{m-1}^{\pm}$ as in Section 5

17:     end if

18:     Calculate $p_{m-1}^{\pm}$ and $p_{0}$ to find the bounds of integration in (4.3)

19:     Compute $\tilde{V}_{m-1}(S(t_{0})e^{x})$ for $x\in\{B_{m-1}^{\pm},\xi_{m-1}^{\pm}\}$ using (4.1) (where the first integral in (4.3) is computed using (4.4), and the last two integrals using (4.6), (4.7), and the existing values of $v_{1}^{\pm},v_{2}^{\pm},y(p_{m}^{-}),y(p_{m}^{+}+p_{0})$), and assign their values to $v_{1}^{\pm},v_{2}^{\pm}$ respectively

20:     

[TABLE]

as in (4.1) (note that $\mathbf{y}$ now stores $\tilde{V}_{m-1}(S(i))$ for $i=1,2,\dotsc,N$)

21:  end for

22:  Compute $\tilde{V}_{0}(S(t_{0}))$ using (4.1), where the first integral in (4.3) is computed using (4.4), and the last two integrals using (4.6), (4.7), and the existing values of $v_{1}^{\pm},v_{2}^{\pm},y(p_{1}^{-}),y(p_{1}^{+}+p_{0})$

Since the computational complexity of each step of the loop is $O(N\log N)$, the total complexity is $O(MN\log N)$.

Remark 6.1**.**

While other quadrature methods typically employ multiple uniform grids or specially-designed (nonuniform or shifted) grids, our method utilizes only a fixed one-dimensional uniform grid, which not only eliminates the need for complicated inter-grid data transfer procedures, but also eliminates the need for special subroutines that are often required to interpolate data across discontinuities. This makes our method particularly easy to implement.

7. Numerical Examples

We will demonstrate the accuracy and efficiency of the proposed method using two examples, in which the value of an autocallable structured product and that of a double barrier option with time-dependent barriers are found.

7.1. Example 1: Autocallable Structured Product

We consider a knock-out autocallable structured product maturing in one year. The price of the underlying asset is 3000, the nominal amount is 1, and the volatility is 20%. The observation dates (in years from now), barrier levels, and risk-free rates (in %) are given below in Table 1.

If the asset price reaches or goes above the barrier level at some observation date $t$, the investor receives a payment of $4\%\times t$. If the asset price is below the barrier at every observation date, the investor will have to pay a premium of $1\%$. The relative errors of the computed option values with varying grid sizes are shown below in Figure 1, where the exact option value is taken to be the one computed on the grid of size 70,001.

As a comparison, the relative errors of the option values computed using Monte-Carlo simulations with antithetic variates technique are shown in Figure 1. As is clear from the figures, the error of the option value computed using the proposed method is well within $10^{-5}$ with just 501 points in the grid, and drops very quickly as the grid size increases. In contrast, it takes more than ten million paths for Monte-Carlo simulations to reduce the error of the computed option value to within $10^{-3}$, and the error decays very slowly as the number of paths increases.

Figure 2 below shows the CPU time used by the proposed method to price the autocallable structured product, where the code is developed in Matlab and is run on a personal computer.

As is clear from the figure, the CPU time required by the proposed method is well within 0.01 seconds, and it increases approximately linearly as grid size increases. It is difficult to compare the speed of the proposed method with that of Monte-Carlo simulations directly, since the CPU time required by the latter depends largely on specific implementations. Nevertheless, the typical CPU time consumed by a Monte-Carlo simulation with tens of millions of paths ranges from tens of seconds to a few minutes.

7.2. Example 2: Double Barrier Option

As another example, consider a knock-out double barrier put option with time-dependent barrier levels. The price of the underlying asset is 2500, the strike price is 2600, the nominal amount is 1, and the volatility is 25%. The option matures in two years. The observation dates (in years from now), barrier levels, and risk-free rates (in %) are given below in Table 2.

If the asset price falls below barrier level 1 or rises above barrier level 2 at any observation date, the option ceases to exist. If the option is still valid at maturity, the payoff is the same as that of a vanilla put option. The relative errors of the computed option values with varying grid sizes are shown below in Figure 3, where the exact option value is taken to be the one computed on the grid of size 50,001.

As a comparison, the relative errors of the option values computed using Monte-Carlo simulations with antithetic variates technique are shown in Figure 3. As is clear from the figures, the error of the option value computed using the proposed method is within $10^{-5}$ with just 701 points in the grid, and drops very quickly as the grid size increases. In contrast, it takes more than ten million paths for Monte-Carlo simulations to reduce the error of the computed option value to within $10^{-3}$, and the error decays very slowly as the number of paths increases.

Figure 4 below shows the CPU time used by the proposed method to price the double barrier option, where the code is developed in Matlab and is run on a personal computer.

As is clear from the figure, the CPU time required by the proposed method is within 0.01 seconds, and it increases approximately linearly as grid size increases. In contrast, the typical CPU time consumed by a Monte-Carlo simulation with tens of millions of paths ranges from tens of seconds to a few minutes.

Remark 7.1**.**

Although the designed order of the proposed method is 4, in the above numerical examples, a lower order of convergence (close to 3) is actually observed for the grid sizes considered. It is interesting to note that this apparent “loss” of order of accuracy is not a defect of our method; rather, it is a manifestation of the subtle influences that barrier levels can have on option pricing algorithms. These influences can be understood from two perspectives. First, in the above examples, the barrier levels $K_{m}^{\pm}$ are close to the spot price $S(t_{0})$. This gives rise to a relatively small integration domain $[B_{m}^{-},B_{m}^{+}]$ compared with the entire computational domain $[{-}\log C,\log C]$ (recall that

[TABLE]

), which means that the set of grid points that are available for the discrete quadrature rule (4.4) represents only a relatively small fraction of the set of grid points introduced on the entire computational domain. Secondly, in the above examples the option prices $V_{m+1}(S)$ contain discontinuities at each observation date $t_{m+1}$. These discontinuities necessarily show up in the form of large gradients in the (smooth) functions $\tilde{V}_{m}(S)$ (via the expectation integrals \mathbb{E}\bigl{[}V_{m+1}(\cdot)|S\bigr{]}), which means that the discrete quadrature rule (4.4) is being applied to fast-varying functions with only a relatively small number of grid points, leaving the integrands only marginally resolved and hence explaining the degeneracy observed in the convergence rate. If the barrier levels $K_{m}^{\pm}$ are pushed farther away from the spot price $S(t_{0})$, so that the option becomes increasingly like a vanilla option, then the discrete quadrature rule (4.4) effectively applies to slow-varying functions with a relatively large number of grid points, which improves the resolution of the integrands and hence the convergence rate (to close to 4). Despite these caveats on convergence rate, we emphasize that our method is capable of pricing a sophisticated discretely monitored option and obtaining five to six significant digits within a fraction of a second, while at the same time being very easy to understand and to implement. Thus the (relatively technical) issue of convergence rate should pose no real concerns in practice.

Remark 7.2**.**

Although a relative error of the order $10^{-5}$ or $10^{-6}$ may not always seem necessary for option pricing problems considered in the real financial world, this extra accuracy is actually needed in the calculation of the Greeks, which are typically approximated by finite difference formulas and which are much more sensitive to numerical errors incurred in the calculation of option prices.

8. Appendix

8.1. Estimate of Truncation Errors

To estimate the truncation error for the integral in (4.1), we first introduce

Lemma 8.1**.**

Let

[TABLE]

Then

[TABLE]

Proof.

Clearly, by assumption (cf. (3.4)),

[TABLE]

Now assume

[TABLE]

for some $1\leq m\leq M$. By (3.2) and Proposition 3.4, we have

[TABLE]

for $K_{m-1}^{-}<S<K_{m-1}^{+}$. By definition (3.1), it is clear that $\rho_{m}(y,S)=\rho_{m}(y/S,1)/S$. Thus a simple change of variable $z=y/S$ yields

[TABLE]

The integral $\int_{0}^{\infty}z\rho_{m}(z,1)\,dz$ is the expectation of the lognormal distribution, which is simply $\exp\{2(r_{m}-q_{m})\tau_{m}/\sigma_{m}^{2}\}$ [19]. Thus we have (observe $R\leq 0$ and $Q\leq 0$)

[TABLE]

for all $S\in(K_{m-1}^{-},K_{m-1}^{+})$. By (3.3) and Proposition 3.4, we then deduce

[TABLE]

for all $S\in(0,\infty)$. The result then follows from induction. ∎

The possibly infinite integral in (3.5) is approximated by a finite integral. To be specific, let $C>1,\ S_{0}=S(t_{0})$, and $\tilde{G}_{M-1}(S)=\tilde{V}_{M-1}(S)$. We consider $\tilde{G}_{m-1}$ ($1\leq m\leq M-1$) defined recursively by

[TABLE]

A direct calculation shows that the errors $\tilde{R}_{m}=\tilde{V}_{m}-\tilde{G}_{m}$ satisfy the recursion

[TABLE]

We use the operator notation

[TABLE]

to write the recursion of $\tilde{R}_{m}$ as

[TABLE]

for $1\leq m\leq M-1$. Note also that $\tilde{R}_{M-1}(S)=0$.

Now we consider $\tilde{Q}_{m}$ defined by the recursion

[TABLE]

and $\tilde{Q}_{M-1}=0$. It is easy to show using Lemma 8.1 and induction that $|\tilde{R}_{m}|\leq\tilde{Q}_{m}$ for all $0\leq m\leq M-1$. We can also apply the recursion formula (8.1) to get the expansion

[TABLE]

Since $\mathcal{T}_{m}$ is the risk-neutral expectation operator discounted from $t_{m}$ to $t_{m-1}$, $\mathcal{T}_{1}\circ\dotsb\circ\mathcal{T}_{m}(f)$ is simply the value of a European option whose payoff at $t_{m}$ is given by $f$. Thus $\tilde{Q}_{0}$ is equal to the value of a sum of $2(M-1)$ binary call options with strike $S_{0}C$ and $2(M-1)$ binary put options with strike $S_{0}/C$. As a result,

[TABLE]

where $C_{a},P_{a},C_{b},P_{b}$ denote values of asset-or-nothing call, asset-or-nothing put, cash-or-nothing call, and cash-or-nothing put respectively. According to Lemma 3.1, the values of these binary options are given by

[TABLE]

where

[TABLE]

and $\bar{r}_{m},\bar{q}_{m},\bar{\sigma}_{m}$ represent the time-weighted averages of $\{r_{n},q_{n},\sigma_{n}\}_{n=1}^{m}$ respectively. Generally, in practice, the absolute values of annual interest and yield rates will not exceed $50\%$, $t_{M}-t_{0}$ will not exceed $10$ years, and $A$ will not exceed $1$. For general autocallable structured products, $B$ will not exceed $t_{M}-t_{0}$, and for Bermudan options $B$ will not exceed $K$, which is not much larger than $S_{0}$. We may also assume $M\leq 120$, which corresponds to products that are not too frequently monitored, say monthly (for more frequently monitored products, such as daily monitored products, continuity correction methods [4, 5] are usually more appropriate). Let $\sigma_{0}=\max_{1\leq m\leq M}\sigma_{m}$. If we choose

[TABLE]

we can make sure that

[TABLE]

A crude estimate using (8.2) then shows that the error bound $\tilde{Q}_{0}(S_{0})$ does not exceed $10^{-15}(S_{0}+1)$. This means the relative truncation error is negligible for all practical purposes.

8.2. Analysis of $K_{m}^{\pm}$ for Bermudan Options

The proper application of Proposition 3.4 requires the uniqueness of the exercise prices $K_{m}^{\pm}$, which we now establish for Bermudan options.

To begin with, observe that the risk-neutral pricing formulas (3.2)–(3.3) applied to Bermudan options can be written in an alternative form as

[TABLE]

where

[TABLE]

and

[TABLE]

Since, by definition,

[TABLE]

(8.3) and (8.4) define $V_{m}$ and $\tilde{V}_{m}$ recursively for all $0\leq m\leq M-1$. It is easy to see that $V_{m}(S)\geq\tilde{V}_{m}(S)>0$ for all $0\leq m\leq M-1$ and $S\in(0,\infty)$.

Proposition 8.2**.**

Assume $q_{m}\geq 0$ for all $1\leq m\leq M$. For a Bermudan call option with strike $K$, the equation $\tilde{V}_{m}(K_{m}^{+})=K_{m}^{+}-K$ has at most one finite solution, and

[TABLE]

For a Bermudan put option with strike $K$, the equation $\tilde{V}_{m}(K_{m}^{-})=K-K_{m}^{-}$ has at most one finite solution, and

[TABLE]

Proof.

Clearly $Q=0$ since $q_{m}\geq 0$. The proofs for call and put options are very similar, and we will present the argument for put options only which proceeds by induction. Since $\tilde{V}_{M-1}$ is the value of a European vanilla put option and $q_{m}\geq 0$, its delta is between $-1$ and 0 [18]. Thus

[TABLE]

for all $S_{2}>S_{1}>0$. Now suppose $\tilde{V}_{M-1}(S_{3})=K-S_{3}$ and $\tilde{V}_{M-1}(S_{4})=K-S_{4}$ for some $S_{4}>S_{3}>0$. This means

[TABLE]

which is a contradiction. So the proposition is true for $m=M-1$.

Suppose now the results hold for some $2\leq m\leq M-1$, which implies that the function $S-K+\tilde{V}_{m}(S)$ is strictly increasing. Since $\tilde{V}_{m}(S)>0$, for sufficiently large $S$ we have $S-K+\tilde{V}_{m}(S)>0$, or $\tilde{V}_{m}(S)>K-S$. This means that

[TABLE]

is well-defined and satisfies $K_{m}^{-}<\infty$. Consider now (cf. (8.3))

[TABLE]

and any $S_{2}>S_{1}>0$. If $S_{2}<K_{m}^{-}$ (and hence $S_{1}<K_{m}^{-}$), we have

[TABLE]

by (8.5). If $S_{1}\geq K_{m}^{-}$ (and hence $S_{2}\geq K_{m}^{-}$), we have

[TABLE]

by (8.5) and inductive hypothesis. If $S_{1}<K_{m}^{-}\leq S_{2}$, we have

[TABLE]

by (8.5), inductive hypothesis, and the definition of $K_{m}^{-}$. In conclusion, we have shown that

[TABLE]

With the aid of (8.4) and (8.6), we write

[TABLE]

where

[TABLE]

As a result,

[TABLE]

since by elementary properties of lognormal distributions [19],

[TABLE]

Now suppose $\tilde{V}_{m-1}(S_{3})=K-S_{3}$ and $\tilde{V}_{m-1}(S_{4})=K-S_{4}$. This means

[TABLE]

so by (8.7) we must have $S_{3}=S_{4}$. The proposition then follows from induction. ∎

Corollary 8.3**.**

Assume a Bermudan put option with strike $K$ has an optimal early-exercise level $K_{m}^{-}>0$ for some $1\leq m\leq M-1$. Then we have $\tilde{V}_{m}(S)>K-S$ for $S>K_{m}^{-}$ and $\tilde{V}_{m}(S)<K-S$ for $S<K_{m}^{-}$. Similarly, for a Bermudan call option we have $\tilde{V}_{m}(S)<S-K$ for $S>K_{m}^{+}$ and $\tilde{V}_{m}(S)>S-K$ for $S<K_{m}^{+}$.

Proof.

We will present the proof for put options only as the argument for call options is similar. It follows from Proposition 8.2 that the function $\tilde{V}_{m}(S)+S-K$ is increasing in $S$. Since $\tilde{V}_{m}(K_{m}^{-})+K_{m}^{-}-K=0$, we have $\tilde{V}_{m}(S)>K-S$ for $S>K_{m}^{-}$ and $\tilde{V}_{m}(S)<K-S$ for $S<K_{m}^{-}$. ∎

Acknowledgements

We are very grateful to Qingshuo Song for his valuable insights and helpful suggestions, as well as to Zhenan Sui for her careful reading and editing of the manuscript.