An efficient convergent lattice method for Asian option pricing with superlinear complexity

Ling Lu; Wei Xu; Zhehui Qian

An Asian option (so-called as this type of option was originally traded in the Asian markets, particularly in Tokyo (Zhang 2003)) is a derivative security whose payoff is based on the average of the prices of the underlying asset over certain periods prior to maturity. Since its payoff depends on the arithmetic average price of the underlying asset, it is therefore useful for hedging transactions whose cost is related to the average price of the underlying asset. In this case the resulting option price does not have a simple closed-form solution under the standard Black–Scholes assumptions. Thus, the evaluation of arithmetic Asian options has become a subfield of its own in the discipline of computational finance. Several solutions have been proposed for the problem, including analytical approximations, Monte Carlo simulation, lattice methods and partial differential equations (PDEs).

In Zhang (2003), Fu et al (1999), Ju (2002) and Fusai (2004), Asian options were approximated by closed-form solutions under various assumptions, but some of these assumptions are too strong, and convergence in these methods is not guaranteed in most cases. Thus, numerical methods are popular in practice.

One popular numerical approach is the Monte Carlo method, surveyed in Glasserman (2004), but it is expensive and suffers from the fact that options may be exercised early, ie, American style. In lattice methods, an auxiliary state vector containing the exact values of the path-dependent features is introduced. However, the number of possible values of the arithmetic average grows exponentially with the number of time steps. It is thus not feasible to track every possible average in the auxiliary vector. Ritchken et al (1993) and Hull and White (1993) proposed lattice methods containing a grid that covers the range of possible averages in the state vector, and interpolating between the nodes of the grid when solving backward through the tree to find the initial value of the claim. A similar set of contracts was studied by Barraquand and Pudet (1996) through a lattice method called the forward shooting grid (FSG). In fact, the Hull–White binomial tree method is a special case of the FSG method.

Compared with the Monte Carlo method, the lattice methods are much faster and handle the early exercise situation more easily, but, due to the interpolation error at each time step, the cumulative effect of the error must be carefully monitored to guarantee convergence. Forsyth et al (2002) demonstrated that the FSG method with linear interpolation and a log-linear allocation converges proportionally to $N^{1.5}$ with a complexity of $O(N^{3.5})$ , where $N$ is the number of time steps. Hsu and Lyuu (2007) proposed an $O(N^{2})$ -time-convergent lattice algorithm, where the optimal number of states for each node of the lattice is given by the Lagrange multiplier. As for the traditional two-dimensional PDE method, Forsyth et al (2002) proposed an $O(N^{3})$ -time algorithm, while the one-dimensional (1D) PDE method (Dubois and Lelièvre 2004) can reduce the complexity to $O(N^{2})$ . Večeř’s (2001) 1D PDE method has a simple form and is easy to implement. Zhang (2001) developed a semi-analytical PDE method consisting of two parts: a closed-form formula and a 1D PDE correction term. However, the running time of these 1D PDE methods increases dynamically when the volatility is high or the maturity is long (Dubois and Lelièvre 2004; Večeř 2001; Zhang 2001). Another drawback of the 1D PDE methods is that they cannot be applied to American-style Asian options.

The willow tree method was first proposed by Curran (2001). Subsequently, Xu et al (2013) proposed a new sampling strategy for the willow tree method to make it more accurate and stable for option pricing, and applied this new strategy to European, American, Asian and American-moving-average options to verify its accuracy and stability. In Xu et al (2013), a generic implementation of the Asian option pricing was studied in order to compare the performance of the new sampling strategy with Curran’s strategy. In this implementation, an auxiliary state vector is introduced that contains all the exact values of the path-dependent feature. The computational cost of this implementation increases exponentially with the number of time steps. Therefore, some technique for reducing the computational complexity from exponential to polynomial has to be explored for Asian option pricing. Due to the high cost of the generic implementation, the number of time intervals is set to be less than ten in Xu et al (2013). In order to reduce the amount of bias in the constructed willow tree, the interval between two observation times was divided into $k$ subintervals. The transition probability matrix of these two observation times can be calculated as the product of the $k$ transition probability matrixes corresponding to $k$ subintervals. This subinterval division technique is not necessary when the number of observation times is sufficiently large, say $N=50$ or $N=100$ .

In this paper, we propose a converged $O(N^{1.5})$ -time numerical algorithm based on the willow tree model (Xu et al 2013) for Asian options pricing. For simplicity, only the continuously monitored Asian options are considered, but our proposed method can be easily extended to the fixed discrete monitoring case by setting each time step as the monitoring time. For example, if only twelve monitoring times are required between the initial time and maturity, we just need to set the number of time steps on the willow tree to be $N=12$ . Then, the willow tree method approximates the possible asset values at each monitoring time step. When we analyze the convergence of our method, we make the same assumptions as in Forsyth et al (2002). It turns out that our algorithm is guaranteed to converge to the true solution of the continuous average model. Furthermore, our algorithm can handle high volatilities or long maturities without a spike in computational time. Our numerical experiments in Section 5 will support all these claims.

Assume that the underlying asset follows a lognormal random walk under the risk-neutral probability measure

\mathrm{d}S=rS\,\mathrm{d}t+\sigma S\,\mathrm{d}Z,

where $r$ is the risk-free interest rate, $\sigma$ is the volatility and $\mathrm{d}Z$ is the increment of a Wiener process. For the continuous arithmetic Asian options, the average is defined as

A=\frac{1}{T}\int_{0}^{\mathrm{T}}S_{\tau}\,\mathrm{d}\tau,

which is equivalent to

\mathrm{d}A=\frac{(S-A)\,\mathrm{d}t}{t}.

Then, define the value of a European Asian option as a function of $S$ , $A$ and $t$ , ie, $V^{\mathrm{c}}\equiv V^{\mathrm{c}}(S,A,t)$ , where $V^{\mathrm{c}}$ satisfies

\frac{\partial V^{\mathrm{c}}}{\partial t}+rS\frac{\partial V^{\mathrm{c}}}{% \partial S}+\frac{\sigma^{2}S^{2}}{2}\frac{\partial^{2}V^{\mathrm{c}}}{% \partial S^{2}}+\frac{(S-A)}{t}\frac{\partial V^{\mathrm{c}}}{\partial A}-rV^{% \mathrm{c}}=0.

(1.1)

In the lattice method, the discrete arithmetic average Asian option is considered, rather than the continuous one. Let the observed discrete average be given by the asset price at discrete observation times $0\equiv t_{0},t_{1},\dots,t_{N}\equiv T$ , with corresponding asset prices $\smash{S^{0},S^{1},\dots,S^{N}}$ . The average price can be evaluated as

A^{n}=\frac{\sum_{i=0}^{n}S^{i}}{n+1},

with $A^{0}=S^{0}$ and $n=0,1,\dots,N$ . In other words, $A^{n+1}$ can be rewritten recursively as

A^{n+1}=A^{n}+\frac{S^{n+1}-A^{n}}{n+2}.

Similarly to the continuous case, we can define the value of an Asian option under discrete averaging as $V(S,A,t)$ . Thus, at times other than the observation dates, the Asian option value satisfies the usual Black–Scholes equation, ie,

V_{t}+\tfrac{1}{2}\sigma^{2}S^{2}V_{SS}+rSV_{S}-rV=0.

(1.2)

On the observation date, due to the no-arbitrage assumption, this implies that

V(S,A^{n+1},t_{n}^{+})=V(S,A^{n},t_{n}^{-}),

(1.3)

where $t_{n}^{+}$ ( $t_{n}^{-}$ ) is the time immediately after (before) the observation date $t_{n}$ , and

A^{n+1}=A^{n}+\frac{S-A^{n}}{n+2}.

(1.4)

Equations (1.2)–(1.4) represent a set of infinitely many one-dimensional PDEs on the domain $0\leq S\leq\infty$ and $0\leq A\leq\infty$ , where the average value $A$ appears as a parameter. The terminal condition is defined by the payoff option function at time $T$ , that is, for the European Asian call option,

V(S,A,T)=\max(A-K,0).

These one-dimensional PDEs (for a given value of $A$ ) depend only on the asset price $S$ . At each observation date $t_{n}$ , these one-dimensional problems exchange information based on conditions (1.3) and (1.4). As the number of observation dates goes to infinity, ie, the time interval between two observation dates approaches zero, the solution of (1.2)–(1.4) converges to the solution of (1.1). In fact, if the time interval $[0,T]$ is divided uniformly into $N$ subintervals, the discrete average solution converges to the continuous average model at a rate of $O(N^{-1})$ (Forsyth et al 2002).

The paper is structured as follows. A brief introduction of the willow tree method for Brownian motion is described in Section 2. Then, we apply the willow tree method to Asian option pricing using a complexity-reduction technique in Section 3. The convergence and error analysis of the willow tree method are discussed in Section 4. Numerical experiments are presented in Section 5 to show the efficiency, accuracy and stability of our proposed methods compared with two binomial tree methods from Hull and White (1993) and Hsu and Lyuu (2007), two PDE methods for European Asian options from Večeř (2001) and Zhang (2001), a PDE method for American Asian options from d’Halluin et al (2005) and the Monte Carlo method. Finally, we give concluding remarks in Section 6.

2 Willow tree for Brownian motion

Figure 1: The willow tree lattice with five space nodes and four time steps.

The willow tree method was first proposed by Curran (2001). Assuming that the underlying asset price follows a Brownian motion, a discrete Markov chain is constructed to simulate the Brownian motion. It was proved by Ho (2000) that as the number of time steps goes to infinity, the whole discrete Markov process converges to a Brownian motion. The willow tree in Figure 1 has five spatial points to represent the asset price at each time step, and four time steps in total from time $t=0$ to time $t=3$ . A transition probability matrix is defined to describe the transition probability from a node at time $t_{n}$ to another node at time $t_{n+1}$ . As shown in Figure 1, the number of spatial nodes is constant at each time step. In fact, the willow tree method described in Curran (2001) and Xu et al (2013) approximates the Brownian motion itself, independent of any financial model. In other words, it can be adapted to other models in addition to geometric Brownian motion.

Here, we summarize how to generate the asset prices under the geometric Brownian motion. Readers may refer to Curran (2001), Xu et al (2013) and Xu and Yin (2014) for more details about willow tree methods for different financial models.

In order to approximate the possible asset prices at time $t_{n}$ , pairs $\{(z_{i},q_{i})\}$ are generated for $i=1,2,\dots,m$ , where $m$ is the number of nodes at each time step, $z_{i}$ is the discrete value of the standard normal distribution and $q_{i}$ is the probability $P(z_{i-1}\leq Z\leq z_{i})$ , where $z_{0}=-\infty$ . Then, one of the $m$ possible values of the asset price $S_{i}^{n}$ can be evaluated as

S^{n}_{i}=S_{0}\exp((r-\tfrac{1}{2}\sigma^{2})t_{n}+\sigma\sqrt{t_{n}}z_{i}).

(2.1)

Given the discrete approximate standard normal distribution, the computational cost to construct each tree node $S_{i}^{n}$ is constant, based on (2.1). Curran (2001) gave a strategy to generate $\{(z_{i},q_{i})\}$ as

z_{i}=N^{-1}\bigg(\frac{i-\frac{1}{2}}{m}\bigg)\quad\text{and}\quad q_{i}=% \frac{1}{m},

where $N(\cdot)$ is the cumulative density function for the standard normal distribution.

When $m$ is small, Curran’s proposed strategy leads to a big bias when pricing deeply out-of-the-money options. To overcome this drawback, Xu et al (2013) proposed a sequence of $\{q_{i}\}$ where

q_{i}=\frac{\tilde{q}_{i}}{\sum_{i=1}^{m}\tilde{q}_{i}}\quad\text{and}\quad% \tilde{q}_{i}=\frac{(i-\tfrac{1}{2})^{\gamma}}{m},\quad i=1,2,\dots,\tfrac{1}{% 2}m,~{}0\leq\gamma\leq 1.

The sequence of $\{z_{i}\}$ is the solution of the following nonlinear least squares problem:

\min_{z_{i}}\bigg[\sum_{i=1}^{m}q_{i}z_{i}^{4}-3\bigg]^{2}\quad\text{such that~{}}\sum_{i=1}^{m}q_{i}z_{i}=0,\quad\sum_{i=1}^{m}q_{i}z_{i}^{2}=1,\quad Z% _{i-1}\leq z_{i}\leq Z_{i},

(2.2)

where $\smash{Z_{i}=N^{-1}(\sum_{j=1}^{i}q_{j})}$ for $i=1,2,\dots,m-1$ , $Z_{0}=-\infty$ and $Z_{m}=\infty$ . As shown in Xu et al (2013), the new strategy can price options with a greater precision than Curran’s strategy. Thus, we adopt this new strategy in our willow tree construction for option pricing (1.2)–(1.4). Table 1 lists the pairs $\{(z_{i},q_{i})\}$ for $m=30$ computed by Curran’s sampling strategy and Xu’s strategy with $\gamma=0.6$ . Due to the symmetry of the $\{(z_{i},q_{i})\}$ , we only record the first fifteen values; the other half can be mirrored.

Table 1: List of {zi,qi} from Curran’s strategy, and from Xu’s strategy with γ=0.6.

	Xu	Curran

$?$	$z_{i}$	$q_{i}$	$z_{i}$	$q_{i}$
01	2.8821	0.0069	2.2692	0.0333
02	2.2081	0.0134	1.6449	0.0333
03	1.8896	0.0182	1.383	0.0333
04	1.6486	0.0222	1.1918	0.0333
05	1.4491	0.0259	1.0364	0.0333
06	1.2749	0.0292	0.9027	0.0333
07	1.1175	0.0323	0.7835	0.0333
08	0.9717	0.0351	0.6745	0.0333
09	0.8341	0.0379	0.5730	0.0333
10	0.7020	0.0405	0.4770	0.0333
11	0.5737	0.0430	0.3853	0.0333
12	0.4474	0.0454	0.2967	0.0333
13	0.3216	0.0478	0.2104	0.0333
14	0.1948	0.0500	0.1257	0.0333
15	0.0655	0.0522	0.0418	0.0333

To estimate the transition probability matrix between two continuous-time steps, a constrained linear programming problem is solved (Curran 2001; Xu et al 2013) so that the discrete Markov process approximated by the willow tree converges to a Brownian motion. Assume the underlying asset price satisfies a geometric Brownian motion. The asset price of the $i$ th node at time $t_{i}$ on the willow tree can be estimated as in (2.1). The option can be priced by backward deduction with the transition probability matrixes from the constrained linear programming problem. The cost of determining the transition probability matrixes is high, but it is independent of the parameters in any financial model and can be computed offline. In other words, given the number of nodes at each time step, $m$ , and a number of time steps, $N$ , the whole willow tree can be constructed in advance and stored in a database as a specific approximation of the Brownian motion. When pricing an option, the appropriate willow tree can be extracted from the database to approximate a financial stochastic process under a Brownian motion, so the online computational time is small. For example, in Section 5, building up a willow tree with $m=30$ and $N=400$ takes about thirty-five seconds, but it can be reused 140 times to price different Asian options in our experiments. In fact, this willow tree can be reused to price more options. Another advantage of the willow tree method is that the total number of nodes increases linearly with the number of time steps, $N$ , rather than quadratically as in the binomial tree method. Thus, it is much faster than the binomial tree method, especially when $N$ is large. In the following section, the willow tree is applied to price a European fixed-strike Asian option via an interpolation technique to approximate the arithmetic average on each path, so the total complexity can be $O(N^{1.5})$ .

3 Willow Tree Method for Asian Option Pricing

In this section, we apply the willow tree method to price a European Asian option with our interpolation scheme. Rather than a state vector containing all the exact arithmetic averages, as in Xu et al (2013), we employ a vector to record possible ranges of the average and approximate the option value by interpolation. We first introduce a basic $O(N^{2})$ -time willow tree algorithm with an interpolation scheme to price the Asian option. Then, we propose some complexity-reduction techniques to further reduce the complexity to $O(N^{1.5})$ .

3.1 Basic algorithm

As described in Section 2, the number of possible asset prices at each time step on the willow tree is constant. Let $\Delta t$ be the discrete time step. The price, $S_{i}^{n}$ , of the $i$ th underlying asset at time $t_{n}=n\Delta t$ can be evaluated, for $i=1,2,\dots,m$ , as

S_{i}^{n}=S^{0}\exp((r-\tfrac{1}{2}\sigma^{2})t_{n}+\sigma\sqrt{t_{n}}z_{i}).

(3.1)

As shown in Table 1, $\{z_{i}\}$ approximates the standard normal distribution. It is obvious that $\smash{S_{m}^{n}}$ is the lowest possible asset price at time $t_{n}$ , while $S_{1}^{n}$ is the highest. Thus, the range of the average asset price from $t_{0}$ to $t_{n}$ can be denoted by $[A_{\text{min}}^{n},A_{\text{max}}^{n}]$ , where

A_{\text{min}}^{n}=\frac{1}{n+1}\bigg(S^{0}+\sum_{q=1}^{n}S_{m}^{q}\bigg)\quad% \text{and}\quad A_{\text{max}}^{n}=\frac{1}{n+1}\bigg(S^{0}+\sum_{q=1}^{n}S_{1% }^{q}\bigg).

(3.2)

Given a constant step size $h$ , we can find two integers, $\smash{k_{\text{min}}^{n}}$ and $\smash{k_{\text{max}}^{n}}$ , such that $\smash{k_{\text{min}}^{n}}$ is the largest integer satisfying $\smash{S^{0}\exp(k_{\text{min}}^{n}h)\leq A_{\text{min}}^{n}}$ , while $\smash{k_{\text{max}}^{n}}$ is the smallest integer satisfying $\smash{A_{\text{max}}^{n}\leq S^{0}\exp(k_{\text{max}}^{n}h)}$ . Then, a possible average price from time $t_{0}$ to $t_{n}$ can be evaluated as

A_{k}^{n}=S^{0}\mathrm{e}^{kh},

(3.3)

where $k$ is an integer belonging to $\smash{[k_{\text{min}}^{n},k_{\text{max}}^{n}]}$ . From (3.3), we see that $A_{k}^{n}$ , the average price from $t_{0}$ to $t_{n}$ , does not depend on the current asset price, $S_{i}^{n}$ , at $t_{n}$ . In other words, although there are $m$ possible asset prices at time $t_{n}$ on the willow tree, they share the same set of average prices $\smash{A_{k}^{n}}$ , $\smash{k\in[k_{\text{min}}^{n},k_{\text{max}}^{n}]}$ , for each node.

Denote by $V(S_{i}^{n},A_{k}^{n},t_{n})$ the computed value of the Asian option with underlying asset price $\smash{S_{i}^{n}}$ and arithmetic average $\smash{A_{k}^{n}}$ at time $t_{n}$ . For simplicity, we let $V_{ik}^{n}\equiv V(S_{i}^{n},A_{k}^{n},t_{n})$ for $n=1,\dots,N-1$ . Next, we show how to price a European Asian call option through backward deduction based on the willow tree.

First, at maturity $t_{N}\equiv N\Delta t=T$ , the average price is $\smash{A_{k}^{N}=S^{0}\mathrm{e}^{kh}}$ , where $\smash{k\in[k_{\text{min}}^{N},k_{\text{max}}^{N}]}$ . The European Asian call option value at time $T$ , $\smash{V_{ik}^{N}}$ , is computed as

V_{ik}^{N}=\max(A_{k}^{N}-K,0)\quad\text{for~{}}i=1,\dots,m.

(3.4)

Note that the option value at time $T$ does not depend on $\smash{S_{i}^{N}}$ , the underlying asset price at $T$ .

Second, we consider the option value at $\smash{t_{N-1}}$ , ie, $(N-1)\Delta t$ . Similar to the situation at time $T$ , there are $m$ possible asset prices $\smash{S_{i}^{N-1}}$ . At each $\smash{S_{i}^{N-1}}$ , we assume there are $\smash{k_{\text{max}}^{N-1}-k_{\text{min}}^{N-1}}$ possible paths from $t_{0}$ to $\smash{t_{N-1}}$ . The average price, $\smash{A_{k}^{N-1}}$ , for each path is determined by (3.3). Our goal is to evaluate the option value on $\smash{S_{i}^{N-1}}$ at time $\smash{t_{N-1}}$ with an average price $\smash{A_{k}^{N-1}}$ . Suppose the asset price moves from $\smash{S_{i}^{N-1}}$ at time $\smash{t_{N-1}}$ to $\smash{S_{j}^{N}}$ at $\smash{t_{N}}$ . Then the average price from $t_{0}$ to $t_{N}$ , with $\smash{S_{j}^{N}}$ the underlying asset price at time $\smash{t_{N}}$ , can be calculated as

\bar{A}_{jk}^{N}=A_{k}^{N-1}+\frac{(S_{j}^{N}-A_{k}^{N-1})}{N+1}.

Since $\smash{A_{\text{min}}^{N}\leq\bar{A}_{jk}^{N}\leq A_{\text{max}}^{N}}$ , there exist two indexes $k^{+}$ and $k^{-}$ such that $\smash{A_{k^{-}}^{N}}\leq\smash{\bar{A}_{jk}^{N}}\leq\smash{A_{k^{+}}^{N}}$ . The corresponding option value of $\smash{\bar{A}_{jk}^{N}}$ at time $T$ can be determined as a linear interpolation of the option values of $\smash{V_{jk^{-}}^{N}}$ and $\smash{V_{jk^{+}}^{N}}$ , ie,

\bar{V}_{jk}^{N}=\alpha_{jk}^{N}V_{jk^{+}}^{N}+(1-\alpha_{jk}^{N})V_{jk^{-}}^{% N},\quad\text{where~{}}\alpha_{jk}^{N}=\frac{\bar{A}_{jk}^{N}-A_{k^{-}}^{N}}{A% _{k^{+}}^{N}-A_{k^{-}}^{N}}.

After considering other possible asset prices at time $\smash{t_{N}}$ , the option value at time $\smash{t_{N-1}}$ with asset price $\smash{S_{i}^{N-1}}$ and average price $\smash{A_{k}^{N-1}}$ can be evaluated as

V_{ik}^{N-1}=\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{N-1}\bar{V}_{jk}^{N}% \quad\text{for~{}}i=1,\dots,m\text{~{}and~{}}k\in[k_{\text{min}}^{N-1},k_{% \text{max}}^{N-1}],

(3.5)

where $\smash{p_{ij}^{N-1}}$ is the transition probability from price $\smash{S_{i}^{N-1}}$ at time $\smash{t_{N-1}}$ to node $\smash{S_{j}^{N}}$ at time $\smash{t_{N}}$ . Based on the option values $\smash{V_{ik}^{N-1}}$ at $\smash{t_{N-1}}$ , we can evaluate the option value at $t_{i}$ ( $i=N-2,\dots,1,0$ ) by backward deduction. The following algorithm gives details of the deduction.

Algorithm 3.1.

Assume the initial asset price $S^{0}$ , strike price $K$ , time to maturity $T$ and risk-free interest rate $r$ are known. Given the number of time steps, $N$ , the number of possible asset prices at each time step, $m$ , and the average price step size, $h$ , the European Asian call option value, $V$ , can be computed by the following steps.

(1) $\Delta t=T/N$ ;

(2) Compute $A_{k}^{N}=S^{0}\mathrm{e}^{kh}$ , $k\in[k_{\text{min}}^{N},k_{\text{max}}^{N}]$ ;

(3) $V_{ik}^{N}=\max(A_{k}^{N}-K,0)$ , $i=1,\dots,m$ ;

(4) for $n = N - 1 : - 1 : 1$

for $k_{\text{max}}^{n}$

Compute $A_{k}^{n}=S^{0}\mathrm{e}^{kh}$ ;

for $}m$

$V_{ik}^{n}=0$ ;

for $}m$

$\bar{A}_{jk}^{n+1}=A_{k}^{n}+\dfrac{(S_{j}^{n+1}-A_{k}^{n})}{n+2}$ ;

$\bar{V}_{jk}^{n+1}=\alpha_{jk}^{n+1}V_{jk^{+}}^{n+1}+(1-\alpha_{jk}^{n+1})V_{% jk^{-}}^{n+1}$ , $\alpha_{jk}^{n+1}=\dfrac{\bar{A}_{jk}^{n+1}-A_{k^{-}}^{n+1}}{A_{k^{+}}^{n+1}-A% _{k^{-}}^{n+1}}$ ;

$V_{ik}^{n}=V_{ik}^{n}+p_{ij}^{n}\bar{V}_{jk}^{n+1}$ ;

end % end for $j$

$V_{ik}^{n}=\mathrm{e}^{-r\Delta t}V_{ik}^{n}$ ;

end % end for $i$

end % end for $k$

end % end for $n$

(5) $V=0$ ;

(6) for $}m$

$\bar{A}_{jk}^{1}=\dfrac{S^{0}+S_{j}^{1}}{2}$ ;

$\bar{V}_{jk}^{1}=\alpha_{j}^{1}V_{jk^{+}}^{1}+(1-\alpha_{j}^{1})V_{jk^{-}}^{1}$ ;

$V=V+q_{j}\bar{V}_{jk}^{1}$ ;

end % end for $j$

(7) $V=\mathrm{e}^{-r\Delta t}V$ ;

Note that at time $t_{0}$ the asset price is known as $S^{0}$ , and $q_{i}$ in step 6 is the transition probability from asset price $S^{0}$ at time $t_{0}$ to asset price $S_{j}^{1}$ at time $t_{1}$ .

The error analysis of Algorithm 3.1 is presented in Section 4. First, we give the complexity analysis. Based on (3.2), the upper bound of the average price at time $t_{n}$ can be computed as

	$\displaystyle A_{\text{max}}^{n}$	$\displaystyle=\frac{1}{(n+1)}\sum_{l=0}^{n}\exp((r-\tfrac{1}{2}\sigma^{2})t_{l% }+\sigma\sqrt{t_{l}}z_{1})$
		$\displaystyle=\frac{1}{(n+1)}\sum_{l=0}^{n}\exp(C_{1}\Delta tl+C_{2}\sqrt{% \Delta t}\sqrt{l})$
		$\displaystyle\leq\frac{1}{(n+1)}\sum_{l=0}^{n}\exp(C_{1}\Delta tl+C_{2}\sqrt{% \Delta t}\sqrt{n})$
		$\displaystyle=\frac{1}{(n+1)}\exp(C_{2}\sqrt{\Delta t}\sqrt{n})\frac{1-\exp(C_% {1}\Delta t(n+1))}{1-\exp(C_{1}\Delta t)}$
		$\displaystyle\approx\exp(C_{1}\Delta tn+C_{2}\sqrt{\Delta t}\sqrt{n}),$

where $C_{1}=(r-\tfrac{1}{2}\sigma^{2})$ and $C_{2}=\sigma z_{1}$ . According to the symmetry of the willow tree method (Xu et al 2013), ie, $z_{1}=-z_{m}$ , the lower bound of the average price can be evaluated as

	$\displaystyle A_{\text{min}}^{n}$	$\displaystyle=\frac{1}{(n+1)}\sum_{l=0}^{n}\exp((r-\tfrac{1}{2}\sigma^{2})t_{l% }-\sigma\sqrt{t_{l}}z_{1})$
		$\displaystyle\geq\frac{1}{(n+1)}\sum_{l=0}^{n}\exp(C_{1}\Delta tl-C_{2}\sqrt{% \Delta t}\sqrt{n})$
		$\displaystyle\approx\exp(C_{1}\Delta tn-C_{2}\sqrt{\Delta t}\sqrt{n}).$

Let $h=C\Delta t$ , where $C$ is a constant. Then integers $\smash{k_{\text{min}}^{n}}$ and $\smash{k_{\text{max}}^{n}}$ can be determined by

k_{\text{max}}^{n}=\bigg\lfloor\frac{C_{1}}{C}n+\frac{C_{2}}{C}\frac{\sqrt{n}}% {\sqrt{\Delta t}}\bigg\rfloor\quad\text{and}\quad k_{\text{min}}^{n}=\bigg% \lceil\frac{C_{1}}{C}n-\frac{C_{2}}{C}\frac{\sqrt{n}}{\sqrt{\Delta t}}\bigg\rceil,

(3.6)

where $\lfloor a\rfloor$ is the largest integer less than $a$ , while $\lceil a\rceil$ is the smallest integer larger than $a$ . This implies that there are $\smash{(k_{\text{max}}^{n}-k_{\text{min}}^{n})}$ average prices for each asset price $S_{j}^{n}$ at time $t_{n}$ . For each average price $A_{k}^{n}$ , we can compute an option value $\smash{V_{ik}^{n}}$ using (3.5), so the complexity at time $t_{n}$ is $\smash{O(\sqrt{n}/\sqrt{\Delta t})}$ , ie, $\smash{O(\sqrt{N}\sqrt{n})}$ for $\Delta t=T/N$ . Therefore, the total complexity of computing the Asian option by Algorithm 3.1 is $\smash{O(\sqrt{N}\sum_{n=1}^{N}\sqrt{n})=O(N^{2})}$ .

3.2 Complexity-reduction technique

In Algorithm 3.1, after the range of the average price from $t_{0}$ to $t_{n}$ , $[A_{\text{min}}^{n},A_{\text{max}}^{n}]$ , is computed by (3.2), the number of possible average prices is evaluated by (3.6) to cover the whole range. Obviously, as the discrete time $t_{n}$ is close to the maturity time $T$ , the range $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ is wide, so the number of possible average prices is large.

In this subsection, we will predetermine the total number of possible average prices on the whole willow tree as a constant, so the complexity can be reduced from $\smash{O(N^{2})}$ to $\smash{O(N^{1.5})}$ . In order to handle the trade-off between efficiency and accuracy, we adopt two techniques to increase the efficiency of the willow tree method without losing accuracy. First, we prune some branches that can be valued by a simple formula rather than backward deduction, in order to further reduce the complexity and improve accuracy without introducing interpolation errors. Moreover, the range of average prices also shrinks. Second, according to the total number of average prices, the corresponding number of average prices, $k_{j}^{n}$ , at each asset price $S_{j}^{n}$ is determined by minimizing the total interpolation error. Then, the reduced range is divided equally into $\smash{k_{j}^{n}}$ subintervals.

Given an average price $\smash{A_{k}^{n}}$ at time $t_{n}$ , it is obvious that the sum of all asset prices from $t_{0}$ to $t_{n}$ is $(n+1)A_{k}^{n}$ . If $\smash{(n+1)A_{k}^{n}\geq(N+1)K}$ , where $K$ is the strike price, then the option value $\smash{V(S_{j}^{n},A_{k}^{n},t_{n})}$ can be evaluated by using a simple formula, instead of backward deduction. The following “easy price theorem” was first proved in Aingworth et al (2000) under the binomial tree scheme. Here, we demonstrate that it also holds under the willow tree scheme.

Theorem 3.2 (Easy price theorem).

Suppose the total price associated with state $\smash{(S^{n}_{j},A_{k}^{n},t_{n})}$ is $(N+1)K+\epsilon$ , ie, $(n+1)A_{k}^{n}\geq(N+1)K$ for some $\epsilon\geq 0$ . Then the option value $\smash{V(S^{n}_{j},A_{k}^{n},t_{n})}$ can be evaluated as

V(S^{n}_{j},A_{k}^{n},t_{n})=\begin{cases}\displaystyle\frac{\epsilon+(N-n)S^{% n}_{j}}{N+1}&\text{for~{}}r=0,\\ \\ \displaystyle\frac{1}{N+1}\mathrm{e}^{-r(N-n)\Delta t}\bigg(\epsilon+S^{n}_{j}% \mathrm{e}^{r\Delta t}\frac{1-\mathrm{e}^{(N-n)r\Delta t}}{1-\mathrm{e}^{r% \Delta t}}\bigg)&\text{for~{}}r>0.\end{cases}

(3.7)

Proof.

If $r\geq 0$ , the expected value of $S^{n+1}$ given the asset price $\smash{S_{j}^{n}}$ based on the willow tree backward deduction can be evaluated as

	$\displaystyle\mathbb{E}(S^{n+1}\mid S_{j}^{n})$	$\displaystyle=\sum_{i=1}^{m}p_{ji}^{n}S_{i}^{n+1}=\sum_{i=1}^{m}p_{ji}^{n}S_{j% }^{n}\mathrm{e}^{\Delta X_{ij}}$
		$\displaystyle=S_{j}^{n}\sum_{i=1}^{m}p_{ji}^{n}(1+\Delta X_{ij}+\tfrac{1}{2}% \Delta X_{ij}^{2}+O(\Delta t^{2}))$
		$\displaystyle=S_{j}^{n}\bigg[\sum_{i=1}^{m}p_{ji}^{n}+\sum_{i=1}^{m}p_{ji}^{n}% \Delta X_{ij}+\hbox{\small$\displaystyle\frac{1}{2}$}\sum_{i=1}^{m}p_{ji}^{n}% \Delta X_{ij}^{2}+O(\Delta t^{2})\bigg]$
		$\displaystyle=S_{j}^{n}[1+(r-\tfrac{1}{2}\sigma^{2})\Delta t+\tfrac{1}{2}% \sigma^{2}\Delta t+O(\Delta t^{2})]$
		$\displaystyle=S_{j}^{n}(1+r\Delta t+O(\Delta t^{2}))$
		$\displaystyle=S_{j}^{n}\mathrm{e}^{r\Delta t},$

where $\Delta X_{ji}$ is defined as $\ln S_{j}^{n+1}-\ln S_{i}^{n}$ in (4.3) with the properties in (4.5). Similarly, we can prove

	$\displaystyle\mathbb{E}(S^{n+2}\mid S_{j}^{n})$	$\displaystyle=\sum_{i=1}^{m}\sum_{l=1}^{m}p_{ji}^{n}p_{il}^{n+1}S_{l}^{n+2}=% \sum_{i=1}^{m}p_{ji}^{n}\bigg[\sum_{l=1}^{m}p_{il}^{n+1}S_{l}^{n+2}\bigg]$
		$\displaystyle=\sum_{i=1}^{m}p_{ji}^{n}\mathbb{E}(S^{n+2}\mid S_{i}^{n+1})=\sum% _{i=1}^{m}p_{ji}^{n}S_{i}^{n+1}\mathrm{e}^{r\Delta t}$
		$\displaystyle=\mathrm{e}^{r\Delta t}\bigg[\sum_{i=1}^{m}p_{ji}^{n}S_{i}^{n+1}\bigg]$
		$\displaystyle=\mathrm{e}^{r\Delta t}\mathbb{E}(S^{n+1}\mid S_{j}^{n})$
		$\displaystyle=\mathrm{e}^{r\Delta t}S_{j}^{n}\mathrm{e}^{r\Delta t}$
		$\displaystyle=S_{j}^{n}\mathrm{e}^{2r\Delta t}.$

So, the expected value of the asset price summation from $\smash{S^{n+1}}$ to $\smash{S^{N}}$ given $\smash{S^{n}_{j}}$ is in the form

	$\displaystyle\mathbb{E}(S^{n+1}+\cdots+S^{N}\mid S_{j}^{n})$	$\displaystyle=S_{j}^{n}(\mathrm{e}^{r\Delta t}+\mathrm{e}^{2r\Delta t}+\cdots+% \mathrm{e}^{(N-n)r\Delta t})$
		$\displaystyle=S_{j}^{n}\mathrm{e}^{r\Delta t}\frac{1-\mathrm{e}^{(N-n)r\Delta t% }}{1-\mathrm{e}^{r\Delta t}}.$

The option value $\smash{V(S_{j}^{n},A_{k}^{n},t_{n})}$ therefore equals

$\displaystyle V(S_{j}^{n},A_{k}^{n},t_{n})$	$\displaystyle=\mathrm{e}^{-r(N-n)\Delta t}\mathbb{E}[\max(A(N)-K,0)\mid S_{j}^% {n}]$
	$\displaystyle=\mathrm{e}^{-r(N-n)\Delta t}\mathbb{E}\bigg[\max\bigg(\frac{(N+1% )K+\epsilon+\sum_{l=n+1}^{N}S^{l}}{N+1}-K,0\bigg)\biggm\|S_{j}^{n}\bigg]$
	$\displaystyle=\mathrm{e}^{-r(N-n)\Delta t}\frac{1}{N+1}\bigg(\epsilon+S^{n}% \mathrm{e}^{r\Delta t}\frac{1-\mathrm{e}^{(N-n)r\Delta t}}{1-\mathrm{e}^{r% \Delta t}}\bigg).$	(3.8)

The case for $r=0$ is similar. ∎

Theorem 3.2 shows that, when $A_{k}^{n}\geq(N+1)K/(n+1)$ , the option value $\smash{V(S_{j}^{n},A_{k}^{n},t_{n})}$ can be evaluated accurately by (3.7) without introducing any interpolation errors. In other words, we only need to calculate the option values for the average prices, $\smash{A_{k}^{n}}$ , that are less than $(N+1)K/(n+1)$ through the backward deduction and interpolation in Algorithm 3.1. Thus, the range for the average price is reduced to $\smash{[A_{\text{min}}^{n},\min\{(N+1)K/(n+1),A_{\text{max}}^{n}\}]}$ , on which an optimal $\smash{k_{j}^{n}}$ is determined to minimize the total interpolation error. This technique improves efficiency by pruning unnecessary average price branches, and accuracy by not resorting to interpolation on these branches.

Another complexity-reduction technique is to determine the corresponding number of average prices, $\smash{k_{j}^{n}}$ , at each asset price $S_{j}^{n}$ , in order to minimize the total interpolation error. First, an average number, $k_{\mathrm{a}}$ , of average prices on each possible price $S_{j}^{n}$ is given, so that the total number of average prices on the whole willow tree can be evaluated as $Nmk_{\mathrm{a}}$ . In other words, the number of paths considered in pricing is no greater than $Nmk_{\mathrm{a}}$ . At the discrete time $t_{n}$ with asset price $\smash{S_{j}^{n}}$ , the corresponding number of average prices is determined by solving the following nonlinear constrained optimization problem to minimize the total interpolation error:

\min_{k^{n}_{j}}\sum_{n=1}^{N}\,\sum_{i=1}^{m}\,\sum_{j=1}^{m}p_{ij}^{n-1}% \frac{1}{(n+1)^{s}(k_{j}^{n})^{s}}\quad\text{such that~{}}\sum_{n=1}^{N}\,\sum_{j=1}^{m}k_{j}^{n}=Nmk_{\mathrm{a}},

(3.9)

where $s$ is the number of points used in the polynomial interpolation. The derivation of (3.9) is described in Section 4 after we analyze the truncation error and interpolation error of Algorithm 3.1. Solving (3.9) using Lagrangian multipliers, given the transition probability $p_{ij}^{n-1}$ , we can obtain a formula for every $k_{j}^{n}$ , that is

k_{j}^{n}=Nmk_{\mathrm{a}}\frac{[\sum_{i=1}^{m}p_{ij}^{n-1}(n+1)^{-s}]^{1/(s+1% )}}{\sum_{n=1}^{N}\,\sum_{j=1}^{m}[\sum_{i=1}^{m}p_{ij}^{n-1}(n+1)^{-s}]^{1/(s% +1)}}\quad\text{for~{}}j=1,\dots,m\text{~{}and~{}}n=1,\dots,N.

(3.10)

Thus, the range $\smash{[A_{\text{min}}^{n},\min\{(N+1)K/(n+1),A_{\text{max}}^{n}\}]}$ at $\smash{S_{j}^{n}}$ can be divided into $k_{j}^{n}$ subintervals, and the end point of each interval represents a possible average price $\smash{A_{k}^{n}}$ for $\smash{k=1,\dots,k_{j}^{n}+1}$ in Algorithm 3.1. Since the total number of average prices is $Nmk_{\mathrm{a}}$ , the complexity of pricing the Asian option by the willow tree is reduced to $O(N^{1.5})$ , since $k_{\mathrm{a}}$ is taken to be $O(\sqrt{N})$ to guarantee the convergence.

4 Convergence Analysis

In this section, we analyze the convergence of the willow tree method and its truncation and interpolation errors when pricing the European Asian call option through Algorithm 3.1 with and without complexity-reduction techniques.

Before analyzing the error from Algorithm 3.1, we first discuss the convergence of the willow tree method in Curran (2001) and Xu et al (2013) for solving the Black–Scholes equation (1.2). For simplicity, let $X\equiv\ln S$ . Then (1.2) can be rewritten as

\tilde{V}_{t}+\tfrac{1}{2}\sigma^{2}\tilde{V}_{XX}+(r-\tfrac{1}{2}\sigma^{2})% \tilde{V}_{X}-r\tilde{V}=0.

(4.1)

The willow tree method can be applied to solve (4.1) by following the process in Curran (2001) and Xu et al (2013):

\tilde{V}_{i}^{n}\equiv\tilde{V}(X_{i}^{n},t_{n})=\mathrm{e}^{-r\Delta t}\sum_% {j=1}^{n}p_{ij}^{n}\tilde{V}_{j}^{n+1},

(4.2)

where $\smash{\tilde{V}_{i}^{n}}$ is the option value at asset price $\smash{X_{i}^{n}}$ on time $t_{n}$ .

Theorem 4.1.

The option value $\tilde{V}(X_{0},t_{0})$ computed by the Willow tree model through (4.2) converges to the solution of (4.1) as $\Delta t\rightarrow 0$ . Furthermore, the global truncation error of solving (4.1) by the willow tree method is $O(\Delta t)$ .

Proof.

Consider a time interval $[t_{n},t_{n+1}]$ . We can define $\Delta X_{ji}$ as follows:

\Delta X_{ji}=\ln S_{j}^{n+1}-\ln S_{i}^{n}=(r-\tfrac{1}{2}\sigma^{2})\Delta t% +\sigma(\sqrt{t_{n+1}}z_{j}-\sqrt{t_{n}}z_{i}).

(4.3)

As is well known, the willow tree based on the process $\{Y_{t_{n}}\equiv\sqrt{t_{n}}z_{i},~{}n=1,2,\dots,N,~{}i=1,\dots,m\}$ converges to a Brownian motion as $N\rightarrow\infty$ and $m\rightarrow\infty$ . Thus, the process $\smash{\{Y_{t_{n}}\}}$ and the corresponding transition probability $\smash{p_{ij}^{n}}$ have the following five properties:

$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}$	$\displaystyle=1,$	(4.4a)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}(\sqrt{t_{n+1}}z_{j}-\sqrt{t_{n}}z_{i})$	$\displaystyle=0,$	(4.4b)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}(\sqrt{t_{n+1}}z_{j}-\sqrt{t_{n}}z_{i})^{2}$	$\displaystyle=\Delta t,$	(4.4c)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}(\sqrt{t_{n+1}}z_{j}-\sqrt{t_{n}}z_{i})^{3}$	$\displaystyle=0,$	(4.4d)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}(\sqrt{t_{n+1}}z_{j}-\sqrt{t_{n}}z_{i})^{4}$	$\displaystyle=3\Delta t^{2}.$	(4.4e)

Thus, from (4.4), we can obtain

$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}\Delta X_{ji}$	$\displaystyle=(r-\tfrac{1}{2}\sigma^{2})\Delta t,$	(4.5a)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}\Delta X_{ji}^{2}$	$\displaystyle=(r-\tfrac{1}{2}\sigma^{2})^{2}\Delta t^{2}+\sigma^{2}\Delta t,$	(4.5b)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}\Delta X_{ji}^{3}$	$\displaystyle=(r-\tfrac{1}{2}\sigma^{2})^{3}\Delta t^{3}+3\sigma^{2}(r-\tfrac{% 1}{2}\sigma^{2})\Delta t^{2},$	(4.5c)
$\displaystyle\sum_{j=1}^{m}p_{ij}^{n}\Delta X_{ji}^{4}$	$\displaystyle=(r-\tfrac{1}{2}\sigma^{2})^{4}\Delta t^{4}+6\sigma^{2}(r-\tfrac{% 1}{2}\sigma^{2})^{2}\Delta t^{3}+3\sigma^{4}\Delta t^{2}.$	(4.5d)

In order to prove that the deduction (4.2) satisfies the Black–Scholes equation (4.1), we first expand $\smash{\tilde{V}_{j}^{n+1}}$ at point $\smash{(X_{i}^{n},t_{n})}$ and $\smash{\mathrm{e}^{-r\Delta t}}$ at point 0 in the form

$\displaystyle\tilde{V}_{j}^{n+1}$	$\displaystyle=\tilde{V}(X_{j},t_{n+1})$
	$\displaystyle=\tilde{V}(X_{i}+\Delta X_{ji},t_{n}+\Delta t)$
	$\displaystyle=\tilde{V}_{i}^{n}+\tilde{V}_{t}\Delta t+\tfrac{1}{2}\tilde{V}_{% tt}\Delta t^{2}+\tilde{V}_{X}\Delta X_{ji}+\tfrac{1}{2}\tilde{V}_{XX}\Delta X_% {ji}^{2}$
	$\displaystyle\qquad+\tfrac{1}{6}\tilde{V}_{XXX}\Delta X_{ji}^{3}+\tfrac{1}{24}% \tilde{V}_{XXXX}\Delta X_{ji}^{4}+\tilde{V}_{Xt}\Delta X_{ji}\Delta t$
	$\displaystyle\qquad+\tfrac{1}{2}\tilde{V}_{XXt}\Delta X_{ji}^{2}\Delta t+O(% \Delta t^{3}),$	(4.6)

and

\mathrm{e}^{-r\Delta t}=1-r\Delta t+\tfrac{1}{2}r^{2}\Delta t^{2}+O(\Delta t^{% 3}).

Then, substituting the above expansions into (4.2), we have

$\displaystyle\tilde{V}_{i}^{n}$	$\displaystyle=\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}\tilde{V}_{j}^{n+1}$
	$\displaystyle=\bigg(1-r\Delta t+\frac{r^{2}\Delta t^{2}}{2}\bigg)$
	$\displaystyle\quad\times\sum_{j=1}^{m}p_{ij}^{n}\bigg(\tilde{V}_{i}^{n}+\tilde% {V}_{t}\Delta t+\frac{\tilde{V}_{tt}\Delta t^{2}}{2}+\tilde{V}_{X}\Delta X_{ji% }+\frac{\tilde{V}_{XX}\Delta X_{ji}^{2}}{2}+\frac{\tilde{V}_{XXX}\Delta X_{ji}% ^{3}}{6}$
	$\displaystyle\qquad\qquad\qquad\qquad+\frac{\tilde{V}_{XXXX}\Delta X_{ji}^{4}}% {24}+\tilde{V}_{Xt}\Delta X_{ji}\Delta t+\frac{\tilde{V}_{XXt}\Delta X_{ji}^{2% }\Delta t}{2}\bigg)$
	$\displaystyle\qquad+O(\Delta t^{3})$
	$\displaystyle=\tilde{V}_{i}^{n}-r\tilde{V}_{i}^{n}\Delta t+\tilde{V}_{t}\Delta t% -r\tilde{V}_{t}\Delta t^{2}+\frac{\tilde{V}_{tt}\Delta t^{2}}{2}+\frac{r^{2}% \tilde{V}\Delta t^{2}}{2}$
	$\displaystyle\qquad+(\tilde{V}_{X}+\tilde{V}_{Xt}\Delta t-r\tilde{V}_{X}\Delta t% )\bigg(r-\frac{\sigma^{2}}{2}\bigg)\Delta t$
	$\displaystyle\qquad+\bigg(\frac{\tilde{V}_{XX}}{2}+\frac{\tilde{V}_{XXt}\Delta t% }{2}-\frac{r\tilde{V}_{XX}\Delta t}{2}\bigg)\bigg(\bigg(r-\frac{\sigma^{2}}{2}% \bigg)^{\!2}\Delta t^{2}+\sigma^{2}\Delta t\bigg)$
	$\displaystyle\qquad+\frac{\tilde{V}_{XXX}3\sigma^{2}}{6}\bigg(r-\frac{\sigma^{% 2}}{2}\bigg)\Delta t^{2}+\frac{\tilde{V}_{XXXX}3\sigma^{4}\Delta t^{2}}{24}+O(% \Delta t^{3}).$	(4.7)

Simplifying (4), we obtain

\tilde{V}_{t}+\tfrac{1}{2}\sigma^{2}\tilde{V}_{XX}+(r-\tfrac{1}{2}\sigma^{2})% \tilde{V}_{X}-r\tilde{V}+\varepsilon=0,

where

	$\displaystyle\varepsilon$	$\displaystyle=(-r\tilde{V}_{t}+\tfrac{1}{2}\tilde{V}_{tt}+\tfrac{1}{2}r^{2}% \tilde{V}+(\tilde{V}_{Xt}-r\tilde{V}_{X})(r-\tfrac{1}{2}\sigma^{2})+\tfrac{1}{% 2}\tilde{V}_{XX}(r-\tfrac{1}{2}\sigma^{2})^{2}$
		$\displaystyle\qquad+\tfrac{1}{2}\sigma^{2}(\tilde{V}_{XXt}-r\tilde{V}_{XX})+% \tfrac{1}{2}\sigma^{2}\tilde{V}_{XXX}(r-\tfrac{1}{2}\sigma^{2})+\tfrac{1}{8}% \sigma^{4}\tilde{V}_{XXXX})\Delta t+O(\Delta t^{2})$
		$\displaystyle=O(\Delta t).$

So, $\varepsilon\rightarrow 0$ when $\Delta t\rightarrow 0$ . ∎

Theorem 4.1 gives the convergence of the willow tree deduction (4.2). Moreover, it shows that the global truncation error of (4.2) is $O(\Delta t)$ . Now, we are ready to consider the interpolation errors in the willow tree scheme when pricing the Asian option by Algorithm 3.1.

Assume that $W(S,A,t)$ is the exact Asian option value. Then, for an asset price $\smash{S_{j}^{n}}$ and average price $\smash{A_{k}^{n}}$ at time $t_{n}$ , we have $\smash{W_{jk}^{n}\equiv W(S_{j}^{n},A_{k}^{n},t_{n})}$ . Given an average price $\smash{\bar{A}_{jk}^{n+1}}$ as in Algorithm 3.1, the exact option value, $\smash{W_{j\bar{k}}^{n+1}}\equiv\smash{W(S_{j}^{n+1},\bar{A}_{jk}^{n+1},t_{n+1% })}$ , is approximated by a linear interpolation, ie,

\bar{W}_{jk}^{n+1}=\alpha_{jk}^{n+1}W_{jk^{+}}^{n+1}+(1-\alpha_{jk}^{n+1})W_{% jk^{-}}^{n+1},

(4.8)

where

\alpha_{jk}^{n+1}=\frac{\bar{A}_{jk}^{n+1}-A_{k^{-}}^{n+1}}{A_{k^{+}}^{n+1}-A_% {k^{-}}^{n+1}}.

The interpolation error of $\smash{W_{j\bar{k}}^{n}}$ is

\beta_{jk}^{n+1}\equiv|W_{j\bar{k}}^{n+1}-\bar{W}_{jk}^{n+1}|=\frac{(A_{k^{+}}% ^{n+1}-\bar{A}_{jk}^{n+1})(A_{k^{+}}^{n+1}-A_{k^{-}}^{n+1})}{2}\frac{\partial^% {2}W_{j\bar{k}}^{n+1}(\eta)}{\partial A^{2}}.

(4.9)

Assume that there exists a positive constant, $M$ , independent of $\Delta t$ , such that $\smash{|\partial^{2}W_{j\bar{k}}^{n+1}(\eta)/\partial A^{2}|\leq M}$ for any $n, j$ and $\smash{\eta\in[A_{k^{-}}^{n+1},A_{k^{+}}^{n+1}]}$ , and that $\smash{\hat{A}}$ is a constant such that the effect of interpolation errors induced at $\smash{A_{k}^{n}>\hat{A}}$ is negligible at $\smash{S^{0}}$ (Forsyth et al 2002). Then (4.9) can be bounded as

|\beta_{jk}^{n+1}|\leq M\hat{A}^{2}(1-\mathrm{e}^{-h})^{2},

since

|A_{k^{+}}^{n+1}-\bar{A}_{jk}^{n+1}|\leq A_{k^{+}}^{n+1}(1-\mathrm{e}^{-h})% \leq\hat{A}(1-\mathrm{e}^{-h}).

Thus, applying (3.5) to $\smash{W_{ik}^{n}}$ , we have

W_{ik}^{n}=\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}[\alpha_{jk}^{n+1}W_% {jk^{+}}^{n+1}+(1-\alpha_{jk}^{n+1})W_{jk^{-}}^{n+1}]+\text{i-err}+\text{t-err},

(4.10)

where

\text{i-err}\equiv\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}\beta_{jk}^{n% +1}

represents the interpolation error in this deduction step, and t-err represents the local truncation error. Let $\smash{E_{ik}^{n}=W_{ik}^{n}-V_{ik}^{n}}$ . Then we can obtain

E_{ik}^{n}=\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}[\alpha_{jk}^{n+1}E_% {jk^{+}}^{n+1}+(1-\alpha_{jk}^{n+1})E_{jk^{-}}^{n+1}]+\text{i-err}+\text{t-err}.

(4.11)

Define the maximum error at time $t_{n}$ as $\|E^{n}\|=\max_{i,k}|E_{ik}^{n}|$ . It follows from (4.11) that

	$\displaystyle\\|E^{n}\\|$	$\displaystyle\leq\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}[\alpha_{jk}^{% n+1}+(1-\alpha_{jk}^{n+1})]\\|E^{n+1}\\|+\text{i-err}+\text{t-err}$
		$\displaystyle\leq\\|E^{n+1}\\|+\text{i-err}+\text{t-err}.$		(4.12)

After $N$ steps of backward deduction, the option price error at initial time $t_{0}$ , $\smash{\|E^{0}\|}$ , can be separated into $\smash{\|E^{0}\|=\|E^{0}\|^{\mathrm{I}}+\|E^{0}\|^{\mathrm{T}}}$ , where $\|E^{0}\|^{\mathrm{I}}$ is the global interpolation error and $\|E^{0}\|^{\mathrm{T}}$ is the global truncation error. Theorem 4.1 demonstrates that $\|E^{0}\|^{\mathrm{T}}=O(\Delta t)$ , while the total interpolation error $\|E^{0}\|^{\mathrm{I}}$ can be bounded as

\|E^{0}\|^{\mathrm{I}}\leq\mathrm{e}^{-r\Delta t}\sum_{n=0}^{N-1}\,\sum_{j=1}^% {m}p_{ij}^{n}\max_{k}|\beta_{jk}^{n+1}|\leq NM\hat{A}^{2}(1-\mathrm{e}^{-h})^{% 2}=O\bigg(\frac{h^{2}}{\Delta t}\bigg).

When taking $h=C\Delta t$ , the error bound $\|E^{0}\|$ can be estimated as $\|E^{0}\|\leq O(\Delta t)$ . It implies that Algorithm 3.1 for pricing the Asian option converges at a rate of $O(\Delta t)$ .

In order to reduce the complexity of Algorithm 3.1, two complexity-reduction techniques, the easy pricing theorem and optimal $k_{j}^{n}$ , were introduced in Section 3.2. In the remainder of this section, we derive the $k_{j}^{n}$ in (4.16) by minimizing the total interpolation errors and the error analysis of Algorithm 3.1 with these two complexity-reduction techniques.

From Theorem 3.2, if the average price $A_{k}^{n+1}$ is larger than $(N+1)K/(n+2)$ , then the option value $\smash{V_{jk}^{n+1}}$ can be evaluated by a simple formula, (3.7), rather than by the interpolation scheme in Algorithm 3.1. Thus, we only need to consider the interpolation error in

\bigg[A_{\text{min}}^{n+1},\min\bigg\{\frac{(N+1)K}{n+2},A_{\text{max}}^{n+1}% \bigg\}\bigg],

which implies that the difference between two adjacent average prices in the interval is no larger than $\smash{(N+1)K/(n+2)k_{j}^{n+1}}$ . Thus, the interpolation error bound (4.9) can be rewritten as

\beta_{jk}^{n+1}\leq M\bigg(\frac{(N+1)K}{(n+2)k_{j}^{n+1}}\bigg)^{\!2}.

The total accumulated interpolation error of Algorithm 3.1 with $\smash{k_{j}^{n+1}}$ can be bounded as

\epsilon_{\mathrm{int}}\leq M\tilde{A}^{2}\sum_{n=0}^{N-1}\,\sum_{i=1}^{m}\,% \sum_{j=1}^{m}p_{ij}^{n}\frac{1}{(n+2)^{2}(k_{j}^{n+1})^{2}},

where $\smash{\tilde{A}=(N+1)K}$ . Therefore, in order to minimize the total interpolation error, $k_{j}^{n+1}$ can be determined from the following constrained optimization problem:

\min_{k_{j}^{n+1}}\sum_{n=0}^{N-1}\,\sum_{i=1}^{m}\,\sum_{j=1}^{m}p_{ij}^{n}% \frac{1}{(n+2)^{2}(k_{j}^{n+1})^{2}}\quad\text{such that~{}}\sum_{n=0}^{N-1}\,\sum_{j=1}^{m}k_{j}^{n+1}=Nmk_{\mathrm{a}}.

In fact, when we adopt an $s$ -point interpolation, rather than the linear interpolation, to approximate $\smash{\bar{W}_{ik}^{n+1}}$ , the accumulated interpolation error at initial time $t_{0}$ can be bounded as

\epsilon_{\mathrm{int}}=M\tilde{A}^{s}\sum_{n=0}^{N-1}\,\sum_{i=1}^{m}\,\sum_{% j=1}^{m}p_{ij}^{n}\frac{1}{(n+2)^{s}(k_{j}^{n+1})^{s}}

(4.13)

through similar theoretical analysis to that in Hsu and Lyuu (2007). The corresponding optimization problem for $\smash{k_{j}^{n+1}}$ can be rewritten as

\min_{k_{j}^{n+1}}\sum_{n=0}^{N-1}\,\sum_{i=1}^{m}\,\sum_{j=1}^{m}p_{ij}^{n}% \frac{1}{(n+2)^{s}(k_{j}^{n+1})^{s}}\quad\text{such that~{}}\sum_{n=0}^{N-1}\,% \sum_{j=1}^{m}k_{j}^{n+1}=Nmk_{\mathrm{a}}.

(4.14)

To solve (4.14), we introduce the Lagrange multiplier, $\lambda$ , and rewrite the objective function as

f=\sum_{n=0}^{N-1}\,\sum_{i=1}^{m}\,\sum_{j=1}^{m}p_{ij}^{n}\frac{1}{(n+2)^{s}% (k_{j}^{n+1})^{s}}+\lambda\bigg(\sum_{n=0}^{N-1}\,\sum_{j=1}^{m}k_{j}^{n+1}-% Nmk_{\mathrm{a}}\bigg).

(4.15)

The optimal $\smash{k_{j}^{n+1}}$ must satisfy

\frac{\partial f}{\partial k_{j}^{n+1}}=\sum_{i=1}^{m}p_{ij}^{n}(n+2)^{-s}(-s)% (k_{j}^{n+1})^{-(s+1)}+\lambda=0.

Thus, the optimal $\smash{k_{j}^{n+1}}$ can be evaluated as

k_{j}^{n+1}=\bigg[\frac{1}{\lambda}s\sum_{i=1}^{m}p_{ij}^{n}(n+2)^{-s}\bigg]^{% 1/(s+1)}.

(4.16)

Substituting (4.16) into the constraint in (4.14), we have

\sum_{n=0}^{N-1}\sum_{j=1}^{m}\bigg[\frac{1}{\lambda}s\sum_{i=1}^{m}p_{ij}^{n}% (n+2)^{-s}\bigg]^{1/(s+1)}=Nmk_{\mathrm{a}},

which implies

\lambda=\bigg(\frac{1}{Nmk_{\mathrm{a}}}\sum_{n=0}^{N-1}\,\sum_{j=1}^{m}[s\sum% _{i=1}^{m}p_{ij}^{n}(n+2)^{-s}]^{1/(s+1)}\bigg)^{\!s+1}.

(4.17)

Therefore, the expression for $\smash{k_{j}^{n+1}}$ is in the form

k_{j}^{n+1}=Nmk_{\mathrm{a}}\frac{[\sum_{i=1}^{m}p_{ij}^{n}(n+2)^{-s}]^{1/(s+1% )}}{\sum_{n=0}^{N-1}\sum_{j=1}^{m}[\sum_{i=1}^{m}p_{ij}^{n}(n+2)^{-s}]^{1/(s+1% )}}.

(4.18)

Since $\sum_{j=1}^{m}p_{ij}^{n}=1$ and $0\leq p_{ij}^{n}\leq 1$ ,

\sum_{j=1}^{m}\bigg(\sum_{i=1}^{m}p_{ij}^{n}\bigg)^{\!1/(s+1)}\leq m^{2/(s+1)}% \equiv\tilde{C}.

We substitute $\smash{k_{j}^{n+1}}$ into (4.13), and the total interpolation error $\epsilon_{\mathrm{int}}$ can be bounded as

	$\displaystyle\epsilon_{\mathrm{int}}$	$\displaystyle\leq M\tilde{A}^{s}(Nmk_{\mathrm{a}})^{-s}\bigg(\sum_{n=0}^{N}(n+% 2)^{-s/(s+1)}\bigg[\sum_{j=1}^{m}\bigg(\sum_{i=1}^{m}p_{ij}^{n}\bigg)^{\!1/(s+% 1)}\bigg]\bigg)^{\!s+1}$
		$\displaystyle\leq M\tilde{A}^{s}\tilde{C}^{s+1}(Nmk_{\mathrm{a}})^{-s}\bigg(% \sum_{n=0}^{N}n^{-s/(s+1)}\bigg)^{\!s+1}$
		$\displaystyle\leq M\tilde{A}^{s}\tilde{C}^{s+1}(Nmk_{\mathrm{a}})^{-s}((s+1)N^% {1/(s+1)})^{s+1}$
		$\displaystyle=M\tilde{C}^{s+1}(s+1)^{s+1}m^{-s}k_{\mathrm{a}}^{-s}N,$

for $\tilde{A}=(N+1)K$ .

In order to guarantee $\epsilon_{\mathrm{int}}$ approaches zero at a rate of $\Delta t$ , we take $k_{\mathrm{a}}=\smash{\sqrt{N}}$ and use four-point interpolation to evaluate $\smash{\bar{W}_{jk}^{n+1}}$ . Now, the total interpolation error is bounded as

\|E^{0}\|^{\mathrm{I}}\leq\epsilon_{\mathrm{int}}\leq O(N^{-1})=O(\Delta t).

On the other hand, from Theorem 4.1, the global truncation error of the willow tree is $O(\Delta t)$ , ie, $\|E^{0}\|^{\mathrm{T}}=O(\Delta t)$ . Thus, it implies that Algorithm 3.1 with the two complexity-reduction techniques converges at a rate of $O(\Delta t)$ . In summary, the complexity-reduction techniques keep the same convergence rate as the basic willow tree method, Algorithm 3.1, but the total complexity for pricing an Asian option is reduced from $O(N^{2})$ to $O(N^{1.5})$ .

5 Numerical Experiments

In this section, we compare the performance of the following seven methods on various contracts of arithmetic average European call options: our proposed willow tree method Algorithm 3.1 (denoted by A1); Algorithm 3.1 with complexity-reduction techniques (CRTs) (A1-CRT); the Hull and White (1993) binomial tree method (HW-BT); the Hsu and Lyuu (2007) binomial tree method (HL-BT); the Monte Carlo method (MC) (see Glasserman 2004); the Večeř (2001) PDE method (Večeř); and the Zhang (2001) semi-analytical PDE method (Zhang). Since the MATLAB code of Večeř’s method is available online,¹¹URL: http://bit.ly/2cqSDMj (accessed March 2014). we take this implementation in our experiments. We also compare our willow tree method with the PDE method in d’Halluin et al (2005) on pricing American Asian fixed-strike call options. All experiments are carried out on a machine with an AMD Phenom X6 1090T 3.2 GHz processor and 16 GB RAM running MATLAB 7.8 under Windows 7 Professional.

As mentioned in Section 2, it is expensive to construct a willow tree including all the transition probability matrixes, but it can be computed offline and is reusable for option pricing in future. For example, it takes about thirty-five seconds to get all the probability transition matrixes for a willow tree with $m=30$ and $N=400$ . However, this tree can be used to price most Asian options in our experiments; our results are recorded in Tables 2–9. In other words, we reuse the willow tree 140 times to price Asian options, so the allocated construction cost for each option pricing is 0.25 seconds. In fact, the more options we price, the smaller the construction cost allocated to each pricing. Thus, in our experiments, we ignore the allocated construction option cost for each pricing. Only the computational times for the option pricing are recorded, given an existing willow tree scheme.

Figure 2: Impact of various m values on the accuracy of Asian option pricing through our proposed willow tree method with CRTs.

Figure 3: Impact of various ka values on the accuracy of Asian option pricing via our proposed willow tree method with CRTs. Impact of various C values

Figure 4: Impact of various C values, where h=C⁢Δ⁢t, on the accuracy of Asian option pricing through our proposed willow tree method with CRTs. Convergence rate of Algorithm 3.1 with and without CRT and HL-BT

Figure 5: Convergence rate of Algorithm 3.1 with and without CRT and HL-BT while ka=⌊3⁢N⌋ for Algorithm 3.1 with CRTs and HL-BT. The three methods converge proportionally to Δ⁢t=1/N. Figure 6: Running times of Algorithm 3.1 with and without CRT and HL-BT

Figure 6: Running times of Algorithm 3.1 with and without CRT and HL-BT while ka=⌊3⁢N⌋ for Algorithm 3.1 with CRTs and HL-BT on a European Asian option pricing. Figure 7: Convergence rate of Algorithm 3.1 with and without CRT and HL-BT

Figure 7: Convergence rate of Algorithm 3.1 with and without CRT and HL-BT while ka=90 for Algorithm 3.1 with CRTs and HL-BT. The three methods converge proportionally to Δ⁢t=1/N. Figure 8: Running times of Algorithm 3.1 with and without CRT and HL-BT

Figure 8: Running times of Algorithm 3.1 with and without CRT and HL-BT while ka=90 for Algorithm 3.1 with CRTs and HL-BT.Table 2: Absolute error of the option value for six numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang) for various N. [The benchmarks Vb are from Zhang’s method with grid size Δ⁢ξ=2⁢ξm/4000, ξm=5⁢σ⁢T3/2 and Δ⁢τ=1/800.]


Case 1: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{1}$ , $\sigma=\text{0.1}$ , $K=\text{100}$ , $V_{\text{b}}=\text{4.9151167}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	1.47E $-$ 02	1.91E $-$ 02	1.39E $-$ 02	3.32E $-$ 02	3.33E $-$ 04	5.33E $-$ 05
100	7.26E $-$ 03	2.66E $-$ 02	7.22E $-$ 03	3.51E $-$ 02	3.45E $-$ 04	5.79E $-$ 05
200	3.38E $-$ 03	2.01E $-$ 02	3.68E $-$ 03	1.48E $-$ 02	3.26E $-$ 04	5.90E $-$ 05
400	1.40E $-$ 03	1.29E $-$ 02	1.86E $-$ 03	5.34E $-$ 03	3.27E $-$ 04	5.93E $-$ 05
800	1.10E $-$ 03	1.00E $-$ 02	9.32E $-$ 04	1.83E $-$ 03	3.29E $-$ 04	5.93E $-$ 05

Case 2: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{1}$ , $\sigma=\text{0.3}$ , $K=\text{95}$ , $V_{\text{b}}=\text{11.6558858}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	1.21E $-$ 02	1.89E $-$ 02	5.87E $-$ 03	1.47E $-$ 02	7.20E $-$ 04	1.04E $-$ 04
100	7.87E $-$ 03	1.84E $-$ 02	2.38E $-$ 03	6.88E $-$ 03	8.13E $-$ 04	8.19E $-$ 05
200	5.27E $-$ 03	1.88E $-$ 02	1.15E $-$ 03	2.22E $-$ 03	8.12E $-$ 04	7.64E $-$ 05
400	3.90E $-$ 03	1.34E $-$ 02	5.50E $-$ 04	3.50E $-$ 04	7.12E $-$ 04	7.50E $-$ 05
800	3.53E $-$ 04	1.12E $-$ 02	2.70E $-$ 04	1.17E $-$ 04	7.96E $-$ 04	7.46E $-$ 05

Case 3: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{3}$ , $\sigma=\text{0.1}$ , $K=\text{105}$ , $V_{\text{b}}=\text{8.3912219}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	6.97E $-$ 02	2.10E $-$ 02	4.09E $-$ 02	7.55E $-$ 03	1.35E $-$ 03	1.66E $-$ 04
100	3.52E $-$ 02	2.35E $-$ 02	2.18E $-$ 02	3.61E $-$ 03	1.35E $-$ 03	1.81E $-$ 04
200	1.74E $-$ 02	1.65E $-$ 02	1.11E $-$ 02	2.26E $-$ 04	1.33E $-$ 03	1.84E $-$ 04
400	8.27E $-$ 03	8.43E $-$ 03	5.61E $-$ 03	8.35E $-$ 04	1.35E $-$ 03	1.85E $-$ 04
800	3.69E $-$ 03	6.24E $-$ 03	2.79E $-$ 03	8.12E $-$ 04	1.35E $-$ 03	1.85E $-$ 04

Case 4: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{3}$ , $\sigma=\text{0.3}$ , $K=\text{95}$ , $V_{\text{b}}=\text{19.0231619}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	2.50E $-$ 02	4.87E $-$ 03	4.48E $-$ 02	3.45E $-$ 03	1.03E $-$ 03	2.72E $-$ 04
100	1.05E $-$ 02	5.82E $-$ 04	2.23E $-$ 02	2.05E $-$ 03	1.40E $-$ 03	3.55E $-$ 04
200	3.04E $-$ 03	5.78E $-$ 03	1.10E $-$ 02	1.45E $-$ 03	1.08E $-$ 03	3.76E $-$ 04
400	5.47E $-$ 04	4.85E $-$ 03	5.35E $-$ 03	1.01E $-$ 03	1.10E $-$ 03	3.81E $-$ 04
800	2.11E $-$ 03	2.58E $-$ 03	2.67E $-$ 03	6.13E $-$ 04	1.06E $-$ 03	3.83E $-$ 04

Table 3: Computational times of six numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř and Zhang) for various N.


Case 1: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{1}$ , $\sigma=\text{0.1}$ , $K=\text{100}$ , $V_{\text{b}}=\text{4.\hbox to 0.0pt{9151167}}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	00.09	0.19	0000.05	00.12	2.57	00.56
100	00.86	0.37	0000.36	00.47	2.78	01.25
200	01.59	0.86	0003.18	02.01	3.14	02.56
400	03.21	2.73	0034.13	08.95	3.73	05.19
800	10.58	5.39	0458.97	41.82	4.54	10.00

Case 2: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{1}$ , $\sigma=\text{0.3}$ , $K=\text{95}$ , $V_{\text{b}}=\text{11.\hbox to 0.0pt{6558858}}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	00.16	0.17	0000.09	00.11	4.91	00.62
100	00.73	0.37	0000.78	00.48	5.23	01.20
200	01.79	0.87	0008.63	02.00	5.30	02.54
400	06.08	2.48	0113.07	08.92	5.55	04.82
800	24.12	5.85	1468.94	41.42	6.08	09.84

Case 3: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{3}$ , $\sigma=\text{0.1}$ , $K=\text{105}$ , $V_{\text{b}}=\text{8.\hbox to 0.0pt{3912219}}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	00.08	0.16	0000.06	00.11	3.04	00.56
100	00.17	0.36	0000.28	00.51	4.21	01.37
200	00.5	0.81	0002.01	02.15	4.79	02.78
400	01.73	2.00	0021.57	09.22	5.04	05.04
800	05.85	5.37	0248.34	42.12	5.41	09.80

Case 4: $S=\text{100}$ , $r=\text{0.09}$ , $T=\text{3}$ , $\sigma=\text{0.3}$ , $K=\text{95}$ , $V_{\text{b}}=\text{19.\hbox to 0.0pt{0231619}}$
$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang
050	00.17	0.22	0000.06	00.12	5.63	00.56
100	00.31	0.39	0000.50	00.48	6.05	01.12
200	01.00	0.81	0005.12	02.01	6.26	02.62
400	03.65	2.01	0058.17	08.99	7.00	05.07
800	14.12	5.58	0868.96	41.92	7.35	09.92

Before the comparison, we first show how the parameters $m$ , $k_{\mathrm{a}}$ and $h$ impact the accuracy of Asian option pricing. Linear interpolation (two-point interpolation) is used in Algorithm 3.1, while the four-point interpolation is employed in Algorithm 3.1 with CRTs. Suppose the initial stock price $S=100$ , strike price $K=90$ , risk-free interest rate $r=0.05$ , volatility $\sigma=0.2$ and maturity $T=1$ . We take the option price $V_{\text{b}}=12.5959916$ from Zhang (2001) as a benchmark. The parameters are initially set up as $m=30$ , $k_{\mathrm{a}}=90$ and $h=0.4\Delta t$ . Then, we let one parameter vary and keep the others fixed in order to show their impact on option pricing. Figure 2 shows the absolute errors between the benchmark and the willow-tree-computed values with various $m$ while the other parameters are fixed. It shows that a large $m$ does not greatly improve the accuracy of the option price. The reason for this phenomenon is that a large $m$ solves a large linear programming problem for transition probabilities $\smash{p_{ij}^{n}}$ , which will be affected by the rounding errors and truncation errors during the computation. Thus, $m$ is set to 30 if not specified in our willow tree methods. Figure 3 shows the error decreases quickly as $k_{\mathrm{a}}$ increases. It means that a big $k_{\mathrm{a}}$ implies an accurate option price, but it also requires a much longer running time. In order to balance the accuracy and complexity, we set $k_{\mathrm{a}}=90$ in the experiments for Tables 4–9. The largest value of $N$ chosen in our experiments is 800, far less than $\smash{k_{\mathrm{a}}^{2}}$ (which is 8100). Thus, according to the results in Section 4, $k_{\mathrm{a}}=90$ is good enough to guarantee the convergence of the willow tree method in all our experiments. Figure 4 illustrates the error with various $C$ from 0.05 to 1, where $h=C\Delta t$ . It implies that $C=0.4$ has the best convergence. Therefore, we take $C$ as 0.4 in other experiments.

Table 4: Computational option values for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) for small σ on S=100, N=400 and T=1 with various r, σ and K. [MC 95% CI denotes 95% confidence interval for Monte Carlo method.]


(a) $\sigma=\text{0.05}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC 95% CI
095	0.05	07.1782	07.1778	07.1779	07.1779	07.1787	07.1777	[7.173, 7.1847]
100		02.7185	02.7211	02.7186	02.7218	02.7192	02.7162	[2.708, 2.7179]
105		00.3425	00.3514	00.3411	00.3435	00.3382	00.3373	[0.3351, 0.3389]
095	0.09	08.8087	08.8087	08.809	08.8089	08.8091	08.8088	[8.8016, 8.8136]
100		04.3098	04.3088	04.3091	04.3102	04.311	04.3082	[4.3022, 4.3136]
105		00.9660	00.9711	00.9622	00.9632	00.9563	00.9584	[0.9551, 0.9617]
095	0.15	11.0945	11.0945	11.0945	11.0945	11.0941	11.0941	[11.0918, 11.1043]
100		06.7952	06.7949	06.7948	06.7948	06.7957	06.7944	[6.7907, 6.8031]
105		02.7505	02.7488	02.7449	02.7442	02.7484	02.7444	[2.7392, 2.7499]

(b) $\sigma=\text{0.1}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC 95% CI
090	0.05	11.9529	11.9522	11.9514	11.9514	11.9518	11.9511	[11.9357, 11.9589]
100		03.6428	03.6502	03.6442	03.6444	03.642	03.6414	[3.6352, 3.6526]
110		00.3338	00.3356	00.3324	00.3319	00.3303	00.3312	[0.3262, 0.3316]
090	0.09	13.386	13.3857	13.3855	13.3855	13.3859	13.3852	[13.3732, 13.3971]
100		04.9165	04.9217	04.917	04.9171	04.9148	04.9152	[4.9054, 4.9254]
110		00.6335	00.6376	00.6311	00.6298	00.6278	00.6303	[0.6247, 0.6325]
090	0.15	15.3985	15.3985	15.3992	15.3992	15.399	15.3988	[15.3839, 15.4088]
100		07.0293	07.0299	07.0284	07.0284	07.029	07.0278	[7.0115, 7.0346]
110		01.418	01.4208	01.4126	01.4101	01.4104	01.4136	[1.4094, 1.4218]

(c) $\sigma=\text{0.2}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC 95% CI
090	0.05	12.5958	12.5988	12.5967	12.5958	12.5952	12.596	[12.595, 12.6379]
100		05.7614	05.7757	05.765	05.7635	05.7631	05.7632	[5.731, 5.7638]
110		01.9893	01.9808	01.9905	01.9891	01.9878	01.99	[1.9662, 1.9862]
090	0.09	13.8317	13.8334	13.8322	13.8315	13.8328	13.8317	[13.7982, 13.8432]
100		06.7759	06.7872	06.7789	06.7776	06.7765	06.7774	[6.7451, 6.7812]
110		02.5455	02.5435	02.5463	02.5448	02.5438	02.5463	[2.5265, 2.5499]
090	0.15	15.6428	15.6433	15.6424	15.642	15.6432	15.6419	[15.6332, 15.6813]
100		08.408	08.4151	08.4097	08.4086	08.4082	08.4089	[8.3994, 8.4404]
110		03.5549	03.5502	03.5546	03.5529	03.5537	03.5556	[3.5352, 3.5639]

(d) $\sigma=\text{0.3}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC 95% CI
090	0.05	13.9503	13.9572	13.9539	13.9527	13.9524	13.9535	[13.9092, 13.9700]
100		07.9409	07.9542	07.9463	07.9449	07.9454	07.9458	[7.9119, 7.9614]
110		04.0683	04.0843	04.0716	04.0703	04.07	04.0717	[4.0432, 4.0801]
090	0.09	14.9811	14.9866	14.9841	14.9831	14.9846	14.9838	[14.9446, 15.0086]
100		08.8243	08.8361	08.8293	08.8279	08.828	08.8288	[8.7964, 8.8498]
110		04.693	04.7083	04.6962	04.6948	04.6949	04.6966	[4.6688, 4.7094]
090	0.15	16.5111	16.5148	16.5132	16.5124	16.5137	16.5134	[16.4887, 16.5574]
100		10.2061	10.2155	10.2101	10.2089	10.2089	10.2102	[10.1836, 10.2427]
110		05.7265	05.7402	05.7289	05.7276	05.7288	05.7301	[5.7024, 5.749]

Table 5: Computational times in seconds for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) for small σ on S=100, N=400 and T=1 with various r and K.

(a) $\sigma=\text{0.05}$

$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
095	0.05	1.61	3.07	018.86	10.39	2.29	5.32	24.16
100		2.26	3.39	018.89	10.33	2.34	4.88	23.57
105		1.84	3.23	018.77	10.25	1.79	4.88	23.24
095	0.09	2.14	3.18	018.77	10.33	1.68	4.90	23.32
100		1.87	3.26	018.72	10.16	1.56	4.99	23.42
105		1.61	3.42	018.75	10.12	2.18	4.90	23.12
095	0.15	2.14	3.28	018.70	10.33	1.64	5.07	22.50
100		2.22	3.23	018.91	10.30	2.26	4.99	23.57
105		1.87	3.21	018.80	10.33	1.70	5.10	22.82

(b) $\sigma=\text{0.1}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
090	0.05	3.00	3.12	034.43	10.34	3.28	5.12	23.20
100		2.81	3.23	034.07	10.34	3.84	4.91	23.42
110		3.09	3.20	034.13	10.39	2.82	4.91	24.16
090	0.09	3.42	3.23	034.21	10.30	2.67	4.79	23.26
100		3.35	3.29	034.15	10.48	3.42	4.80	23.38
110		2.95	3.06	034.16	10.23	3.51	4.74	23.24
090	0.15	3.06	3.37	034.10	10.22	3.59	4.70	23.09
100		2.89	3.24	034.46	10.33	3.65	4.91	23.21
110		2.98	3.15	034.41	10.31	2.46	5.10	23.17

(c) $\sigma=\text{0.2}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
090	0.05	4.82	3.14	066.38	10.37	4.18	4.93	23.09
100		4.87	3.09	066.52	10.44	4.29	4.32	23.49
110		4.54	3.23	066.50	10.34	5.05	4.87	23.46
090	0.09	5.32	3.04	066.46	10.33	4.59	4.93	23.51
100		4.91	3.32	066.46	10.44	4.38	5.02	23.54
110		4.71	3.1	066.47	10.41	4.21	4.93	23.57
090	0.15	4.66	3.29	066.47	10.31	4.65	5.09	23.63
100		4.49	3.28	066.52	10.42	4.91	4.91	24.24
110		4.96	3.42	066.64	10.26	4.54	5.01	23.68

(d) $\sigma=\text{0.3}$
$?$	$?$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
090	0.05	6.55	3.24	111.92	10.31	6.29	5.09	22.84
100		6.6	3.20	113.33	10.58	6.21	4.96	23.01
110		6.58	3.18	113.83	10.44	4.57	4.98	23.24
090	0.09	6.4	3.20	112.27	10.55	4.35	4.98	23.93
100		6.32	3.37	112.04	10.34	3.88	4.95	23.34
110		6.4	3.42	114.19	10.33	4.65	5.05	23.21
090	0.15	6.6	3.32	113.94	10.51	4.4	4.99	23.77
100		6.66	3.18	114.55	10.39	5.19	5.23	23.99
110		6.18	3.23	113.49	10.53	4.73	5.09	24.46

Figure 5 plots the absolute errors of Algorithm 3.1, Algorithm 3.1 with CRTs and HL-BT on pricing a European Asian call option with $S=100$ , $K=95$ , $T=1$ , $r=0.09$ and $\sigma=0.1$ . The corresponding benchmark for this option is $V_{\text{b}}=8.9118509$ (Zhang 2001). In Algorithm 3.1 with CRTs and HL-BT, $k_{\mathrm{a}}$ is chosen as $\lfloor 3\smash{\sqrt{N}}\rfloor$ . It shows that these three algorithms converge proportionally to $1/N$ , ie, $O(\Delta t)$ , which supports our theoretical results in Section 4. Figure 6 records the corresponding running time of these three methods. The running times of HL-BT and Algorithm 3.1 are quadratic.

Table 6: Computational option prices for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) for large σ on S=100, r=0.09, N=400 and T=1 with various K. [MC 95% CI denotes 95% confidence interval for Monte Carlo method.]

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC 95% CI
095	0.4	13.5039	13.513	13.5103	13.5092	13.5098	13.5114	[13.4329, 13.5105]
100		10.9167	10.9275	10.9234	10.9222	10.9229	10.9238	[10.8706, 10.9425]
105		08.7231	08.7352	08.7293	08.7281	08.7293	08.73	[8.6797, 8.7452]
095	0.5	15.4332	15.4422	15.4414	15.4404	15.4419	15.4437	[15.3488, 15.4468]
100		13.0186	13.0287	13.0269	13.0258	13.0282	13.0285	[12.9845, 13.0766]
105		10.9206	10.9314	10.9282	10.9271	10.9292	10.9298	[10.8899, 10.976]
095	0.6	17.3943	17.4036	17.4041	17.4032	17.4052	17.4059	[17.3709, 17.491]
100		15.1165	15.1265	15.1262	15.1253	15.1271	15.1291	[15.0389, 15.1528]
105		13.1026	13.1131	13.1116	13.1107	13.1123	13.1142	[13.0847, 13.193]
095	0.8	21.3314	21.3461	21.3459	21.3452	21.3309	21.3512	[21.2667, 21.4371]
100		19.2708	19.286	19.2848	19.2841	19.262	19.2895	[19.2438, 19.4085]
105		17.4068	17.4224	17.42	17.4193	17.387	17.425	[17.3009, 17.4596]
095	1	25.22	25.2536	25.246	25.2453	25.1031	25.2541	[25.1852, 25.4161]
100		23.3364	23.3708	23.3616	23.361	23.175	23.3692	[23.2252, 23.4512]
105		21.6081	21.6034	21.6324	21.6318	21.3933	21.6406	[21.383, 21.6021]

Table 7: Computational times in seconds for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) for large σ on S=100, r=0.09, N=400 and T=1 with various K.

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
095	0.4	07.99	3.01	160.54	09.92	5.09	4.49	23.56
100		08.00	2.98	162.57	09.83	4.93	4.41	23.48
105		08.77	3.32	159.29	09.97	5.44	4.59	23.37
095	0.5	11.19	3.2	210.84	10.67	5.12	5.71	24.12
100		11.25	3.37	211.91	10.48	4.09	5.37	23.88
105		10.84	3.4	207.5	10.44	5.34	5.57	23.77
095	0.6	14.16	3.26	253.33	10.47	4.84	5.3	24.24
100		14.46	3.6	250.58	10.45	5.82	5.82	23.34
105		12.68	3.29	254.48	10.9	3.67	5.85	23.21
095	0.8	19.86	3.23	337.73	10.37	5.52	4.87	23.28
100		17.72	3.29	326.21	10.06	7.21	4.48	23.67
105		17.36	3.29	320.5	09.91	6.05	5.44	23.09
095	1	23.51	3.12	398.66	09.89	6.46	4.48	22.74
100		22.96	3.39	394.23	10.05	6.41	4.49	23.31
105		25.23	3.15	397.57	09.84	6.43	4.49	22.85

Table 8: Computational option prices of long-tenor option for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) on S=100, r=0.09, N=400 and T=3 with various K.

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
095	0.05	15.1175	15.1175	15.1176	15.1176	15.1163	15.1163	[15.0972, 15.1205]
100		11.3055	11.305	11.3051	11.3049	11.3043	11.3036	[11.2962, 11.3194]
105		04.1905	04.1749	04.1793	04.1606	04.1702	04.1688	[4.1574, 4.1771]
095	0.1	15.218	15.2163	15.2163	15.2147	15.2144	15.2139	[15.1962, 15.2419]
100		11.6428	11.64	11.6413	11.6375	11.639	11.6377	[11.6081, 11.6519]
105		008.3995	08.3955	08.3968	08.3888	08.3926	08.3914	[8.376, 8.4164]
095	0.2	16.6404	16.6394	16.6426	16.6367	16.6364	16.6378	[16.5932, 16.678]
100		13.7702	13.7701	13.7728	13.7657	13.7669	13.7677	[13.7258, 13.8061]
105		11.2237	11.2243	11.226	11.2175	11.2198	11.2203	[11.1827, 11.2576]
095	0.4	21.7354	21.7385	21.7451	21.7386	21.7393	21.7435	[21.6788, 21.8488]
100		19.5828	19.5862	19.5924	19.5857	19.5869	19.5898	[19.5371, 19.7024]
105		17.6203	17.6241	17.6294	17.6225	17.6234	17.6262	[17.5559, 17.7151]
095	0.8	32.91	33.067	33.028	33.0237	32.3552	33.0389	[32.898, 33.3563]
100		31.4149	31.5748	31.532	31.5277	30.7298	31.5436	[31.3102, 31.756]
105		30.0138	30.1763	30.1298	30.1254	29.1843	30.1396	[29.7947, 30.2392]

Figure 7 plots the absolute errors of Algorithm 3.1, Algorithm 3.1 with CRTs and HL-BT with the same parameters as in the previous experiment, except with $k_{\mathrm{a}}=90$ instead of $k_{\mathrm{a}}=\lfloor 3\smash{\sqrt{N}}\rfloor$ . In this case, these three algorithms still converge proportionally to $1/N$ . Figure 8 records the corresponding running time of these three methods. The running times of HL-BT and Algorithm 3.1 are still quadratic, while the running time of Algorithm 3.1 with CRTs is linear to $N$ since $k_{\mathrm{a}}$ is a constant in this case.

In the following, we compare the performance of A1 and A1-CRT methods with the performance of binomial tree, PDE and Monte Carlo methods. To be consistent with the result in Večeř (2001), we take 200 evenly distributed space points ( $M=200$ ) and 400 evenly distributed time points ( $N=400$ ) for two PDE methods (Večeř and Zhang) in our experiments, if not specified. In the implementation of Večeř’s method, the MATLAB built-in function pdepe is used to solve the PDE in which time integration is handled by the ode15s function. Although ode15s is only suitable for low to medium accuracy for stiff problems, there is no MATLAB built-in solver with high accuracy (http://bit.ly/2cX4UMg).

In the first experiment, we compare the performance of these six methods with increasing $N$ on four different European Asian call options from Zhang (2001). The computed prices, $V_{\text{b}}$ , from Zhang (2001) are employed as benchmarks. The benchmarks are calculated with 4000 grid points in the spatial dimension and $\Delta t=1/800$ in the time dimension to achieve a precision of $10^{-7}$ (Zhang 2001). Table 2 lists the absolute errors between benchmarks and the computed prices from these six methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang). Table 3 lists the corresponding computational times of all six methods. It shows that Zhang’s method provides results closest to the benchmarks, because the benchmarks are computed from the same method but with a different number of space grids and time steps. As for other methods, our proposed A1 and A1-CRT methods can also obtain Asian option prices close to the benchmarks. If a quick estimation is required, our willow tree can return the price accurate to within one US cent in under one second, taking $\Delta t=1/100$ , ie, $N=100$ . As is well known, the complexity of the PDE method, ie, Večeř and Zhang, is $O(NM^{2})$ in theory. Due to convergence and stability, theory requires that $M^{2}$ is in the order of $N$ , ie, $M^{2}=O(N)$ . Thus, the total complexities of the Večeř and Zhang methods are $O(N^{2})$ , the same as the A1 and HL-BT methods. In our experiment, the running times of Zhang’s method increase linearly with $N$ , since $M$ is fixed to 200 in all experiments. In the implementation of Večeř’s method, we employ the MATLAB built-in functions pdepe and pdeval to solve the PDE. When the problem is large, MATLAB itself may need some performance optimization to improve its efficiency. Thus, the running times of Večeř’s method increase linearly with $N$ , rather than quadratic with $N$ as expected.

In the second experiment, we compare the performance of the seven methods on a set of European Asian call options with $S=100$ , $T=1$ and various $r$ , $\sigma$ and $K$ . The computed option prices and computational times are listed in Tables 4 and 5, respectively. Since there are no precise values for these options, we compute the 95% confidence intervals using the Monte Carlo method with one million simulations. Table 4 shows that almost all computed prices fall within the confidence interval. This implies that for all computed values, given the number of time steps, $N$ , the computational time of A1 is affected by the volatility, $\sigma$ . When $\sigma$ is small, say less than 0.2, the average price range $[A_{\text{min}}^{n},A_{\text{max}}^{n}]$ is narrow, so the number of average prices required in $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ is small. Thus, the A1 algorithm is faster than the other methods. When $\sigma$ is large, the number of average prices considered in the wide range $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ is large, such that A1-CRT is the fastest method of the seven. Thus, it implies that, given a fixed number of time steps $N=400$ , when $\sigma$ is small, say less than 0.2, the A1 and Večeř methods are good approaches to computing option prices because the average price range $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ is narrow due to the underlying asset pricing evaluation (3.1). However, when $\sigma$ is large, the A1-CRT method is the best.

The third experiment compares the seven methods under some extremely high volatility situations. In this experiment, the volatility values vary from 0.4 to 1 with $S=100$ , $r=0.09$ and $T=1$ . The computed prices and computational times are recorded in Tables 6 and 7, respectively. These illustrate that the computed prices almost all fall within the 95% confidence interval. Our A1-CRT method is the fastest of the seven: it takes only about 50–70% of the computational time of the two PDE methods. It also turns out that the A1-CRT method performed stably in the high volatility situations. As for Algorithm 3.1, since the average price range at time $t_{n}$ , $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ , widens when $\sigma$ is large, which leads to a large number of average price points in $\smash{[A_{\text{min}}^{n},A_{\text{max}}^{n}]}$ , it requires a longer running time. However, with the complexity-reduction technique, the running time only increases slowly as volatility increases, because the total number of average prices on the willow tree is constant up to $N$ given a fixed $k_{\mathrm{a}}$ . In other words, our proposed Algorithm 3.1 with CRTs is stable in terms of both accuracy and efficiency when the volatility is high.

Table 9: Computational times in seconds of a long-tenor option for seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Večeř, Zhang and Monte Carlo) on S=100, r=0.09, N=400 and T=3 with various K.

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	Večeř	Zhang	MC
095	0.05	1.11	2.37	013.15	10.08	3.12	6.88	22.25
100		1.01	2.45	012.28	09.42	3.54	4.38	22.56
105		1.01	2.37	012.37	09.48	3.01	4.38	23.21
095	0.1	1.56	2.53	020.05	09.55	4.12	4.31	22.74
100		1.67	2.53	020.84	09.36	3.71	4.21	22.46
105		2.4	2.29	019.75	09.38	4.23	4.38	22.59
095	0.2	2.7	2.48	037.6	09.39	4.7	4.37	22.84
100		3.51	2.57	037.64	09.3	4.82	4.38	22.73
105		3.28	2.36	037.53	09.31	4.65	4.31	22.54
095	0.4	4.2	2.43	071.12	09.36	6.46	4.4	22.42
100		4.29	2.29	070.81	09.44	5.97	4.27	22.71
105		4.95	2.45	073.35	09.72	6.71	5.85	22.98
095	0.8	8.03	2.43	185.86	10.08	8.14	4.77	22.96
100		8.21	2.56	189.24	09.98	7.07	6.96	23.06
105		8.6	2.56	193.69	10.22	8.72	4.96	23.26

Table 10: Computed Asian option prices monitored discretely weekly via seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Monte Carlo, Večeř and Zhang) at S=100, r=0.09, N=12 and T=0.25 with various K and σ.

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	MC	Večeř	Zhang
095	0.1	6.0145	5.9922	6.0115	6.0344	6.0132	6.0063	6.0174
100		1.775	1.7683	1.7839	2.9923	1.7763	1.748	1.7562
105		0.1445	0.1343	0.1257	1.1961	0.1372	0.125	0.1272
095	0.2	6.3705	6.3771	6.376	6.5085	6.3787	6.3475	6.3678
100		2.8337	2.8724	2.8615	3.4982	2.8653	2.8045	2.8197
105		0.8949	0.9056	0.8955	1.6289	0.9095	0.8811	0.8898
095	0.4	7.9549	8.0105	8.0017	8.2055	8.0225	7.9474	7.9835
100		5.0213	5.1019	5.0827	5.3389	5.1084	5.0069	5.037
105		2.9411	3.0171	2.9859	3.2766	3.0140	2.9281	2.9517

Table 11: Computational times of pricing Asian options monitored discretely weekly via seven numerical methods (A1, A1-CRT, HW-BT, HL-BT, Monte Carlo, Večeř and Zhang) on S=100, r=0.09, N=12 and T=0.25 with various K and σ.

$?$	$\sigma$	A1	A1-CRT	HW-BT	HL-BT	MC	Večeř	Zhang
095	0.1	0.02	0.03	0.02	0.02	0.28	0.03	0.81
100		0.02	0.03	0.02	0.02	0.28	0.05	0.95
105		0.02	0.03	0.02	0.02	0.3	0.03	0.81
095	0.2	0.02	0.03	0.02	0.02	0.48	0.03	0.75
100		0.02	0.03	0.02	0.02	0.53	0.05	0.77
105		0.03	0.03	0.02	0.02	0.5	0.03	0.77
095	0.4	0.02	0.03	0.02	0.02	0.69	0.03	0.89
100		0.02	0.03	0.02	0.02	0.66	0.03	0.78
105		0.02	0.03	0.02	0.02	0.66	0.03	0.7

Table 12: Pricing on American Asian fixed-strike call options via three numerical methods (A1, HW-BT, DFL) with S=100 and r=0.1. [The d’Halluin–Forsyth–Labahn (DFL) method is based on the paper by d’Halluin et al (2005), with both the number of grids of average price and the number of grids of stock price set at 201. The relative errors are computed with respect to the computed values from DFL.]

$\sigma$	$?$	$?$	A1	HW-BT	DFL	A1	HW-BT	A1	HW-BT	DFL
						Relative
			Value	error (%)	Time (s)

0.2	0.5	095	08.8516	08.8479	08.9342	$-$ 0.92	$-$ 0.97	06.74	105.75	58.39
		100	04.8723	04.8746	04.8879	$-$ 0.32	$-$ 0.27	06.66	104.77	58.27
		105	02.3087	02.31	02.312	$-$ 0.14	$-$ 0.09	06.57	103.87	58.27
	1	095	11.2603	11.2595	11.3248	$-$ 0.57	$-$ 0.58	04.46	066.63	58.17
		100	07.529	07.5317	07.5456	$-$ 0.22	$-$ 0.18	04.48	066.6	58.14
		105	04.7231	04.7253	04.7282	$-$ 0.11	$-$ 0.06	04.48	066.74	58.25
0.4	0.5	095	11.9635	11.9653	12.0507	$-$ 0.72	$-$ 0.71	14.07	233.81	58.28
		100	08.4917	08.497	08.5329	$-$ 0.48	$-$ 0.42	13.74	239.32	62.76
		105	05.8669	05.8714	05.893	$-$ 0.44	$-$ 0.37	14.37	241.13	59.75
	1	095	15.7158	15.7206	15.7833	$-$ 0.43	$-$ 0.4	08.21	162.19	59.61
		100	12.4684	12.4753	12.5088	$-$ 0.32	$-$ 0.27	08.5	162.23	60.23
		105	09.8019	09.8081	09.8324	$-$ 0.31	$-$ 0.25	08.83	165.99	62.14
0.6	0.5	095	15.435	15.4403	15.5143	$-$ 0.51	$-$ 0.48	23.28	351.08	58.47
		100	12.2128	12.2204	12.2626	$-$ 0.41	$-$ 0.34	24.59	355.92	60.68
		105	09.5952	09.602	09.6332	$-$ 0.39	$-$ 0.32	24.73	367.37	62.01
	1	095	20.642	20.6503	20.7154	$-$ 0.35	$-$ 0.31	14.51	258.01	63.76
		100	17.6399	17.6495	17.6937	$-$ 0.3	$-$ 0.25	13.56	261.18	62.99
		105	15.0623	15.071	15.1073	$-$ 0.3	$-$ 0.24	14.76	261.16	62.35

Next, we consider options whose maturity is long; for example, $T=3$ , rather than $T=1$ as in the previous experiments. Table 8 and 9 present the results for very long tenor options ( $T=3$ ). Unlike the previous experiments, when $\sigma=0.8$ , some computed option prices from Večeř are outside the 95% confidence interval, while others remain within the interval. This implies that the current implementation of Večeř’s method cannot price the Asian options to a high degree of accuracy when the tenor is long. In our experiment, we take the MATLAB built-in implementation ode15s, which is the highest accuracy algorithm built-in to MATLAB for stiff problems. In other words, some specific implementation of Večeř’s method has to be made to keep the option prices within the 95% confidence interval. Just as we see in the first experiment, when $\sigma<0.2$ our A1 method is the fastest. It takes less than two seconds to calculate the price. However, when $\sigma>0.2$ , the A1-CRT method is the fastest. Also, as expected, the A1-CRT method is not affected by the high volatility. Therefore, our two proposed willow tree methods have a very high performance in terms of both accuracy and efficiency. The optimal strategy to price the European Asian options is to employ the A1 method for low volatility cases, say $\sigma<0.2$ ; the A1-CRT method is optimal for pricing the Asian options as $\sigma\geq 0.2$ .

As is well known, pricing Asian options with a fixed discrete monitoring interval is a key practical problem. In this experiment, we consider Asian options that expire in three months and are monitored weekly. So, there are twelve time steps in total. Results from the Monte Carlo method with one million simulations are used as benchmarks. Tables 10 and 11 record the computed option prices and the computation times for all seven methods. They show that our willow tree methods have results closer to the benchmarks than the other methods, but with the second smallest computational times of the seven methods. In other words, our willow tree methods also perform well in both efficiency and accuracy on pricing the fixed discrete monitoring Asian options.

Finally, we compare our willow tree method (A1) with the HW-BT and PDE methods (DFL) (d’Halluin et al 2005) on pricing the American-style Asian options. Algorithm 3.1 with CRTs and HL-BT are not applicable in this case because the easy price theorem does not hold for American-style options. For the American Asian fixed-strike call options, we set $S=100$ and $r=0.1$ , while parameters $\sigma$ , $T$ and $K$ are variable. For the DFL method, both the grid average price and grid stock price are set at 201, just as in d’Halluin et al (2005), with uniformly distributed points. The computed option prices and computing times are listed in Table 12. We also calculate the relative errors of the results from A1 and HW-BT compared with the results from DFL. Table 12 shows that these two tree methods have a similar accuracy, since the absolute values of their relative errors are less than 1%. However, the A1 method requires a much shorter computational time than the HW-BL and DFL methods. On the other hand, our willow tree method can provide a good approximation of the American-style Asian option price, and it only takes 6–20% of the computational time of the HW-BT and DFL methods. As $\sigma$ increases, the running time of A1 and HW-BT increases, while DFL does not change much because the number of average price grids and stock prices are fixed in DFL.

6 Conclusion

Asian options have a payoff that depends strongly on the historical information on the underlying asset price. The resulting option price does not have a simple closed-form solution under the standard Black–Scholes assumptions. Numerical methods are popular alternatives for option pricing. The fastest lattice and PDE methods, as is well known, have a complexity of $O(N^{2})$ . In this paper, we propose what is, to the best of our knowledge, the first algorithm to have a complexity less than $O(N^{2})$ , based on the willow tree method using interpolation at the same convergence rate as the binomial tree methods (Hull and White 1993; Hsu and Lyuu 2007) and PDE methods (Forsyth et al 2002; Večeř 2001; Zhang 2001).

In our proposed algorithm, a state vector is employed to record all possible average prices on each underlying asset price at every observation time. Then, we adopt interpolations to estimate the value of an option with the average price of entries of the state vector. We also propose a strategy to determine the optimal length of the state vector for each observation time, to minimize the total interpolation error so as to reduce the complexity to $O(N^{1.5})$ . We then analyze the convergence rate and total errors of the willow tree method. It turns out that our proposed method can provide the same order of convergence and accuracy as the binomial tree methods (Hull and White 1993; Hsu and Lyuu 2007; Barraquand and Pudet 1996) and PDE methods (Forsyth et al 2002; Večeř 2001; Zhang 2001), but in much less computational time. Finally, our numerical experiments on European and American Asian options demonstrate the efficiency, accuracy and stability of our proposed method.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Acknowledgements

This work was supported by the Fundamental Research Funds for the Central Universities in China.

References

Aingworth, D., Motwani, R., and Oldham, J. D. (2000). Accurate approximations for Asian options. In Proceedings of the Eleventh Annual ACM–SIAM Symposium on Discrete Algorithms, pp. 891–900. SIAM, Philadelphia, PA.

Barraquand, J., and Pudet, T. (1996). Pricing of American path-dependent contingent claims. Mathematical Finance 6(1), 17–51 (http://doi.org/d7trz4).

Curran, M. (2001). Willow power: optimizing derivative pricing trees. Algo Research Quarterly 4(4), 15–24.

d’Halluin, Y., Forsyth, P. A., and Labahn, G. (2005). A semi-Lagrangian approach for American Asian options under jump diffusion. SIAM Journal on Scientific Computing 27, 315–345 (http://doi.org/c7w9qh).

Dubois, F., and Lelièvre, T. (2004). Efficient pricing of Asian options by the PDE approach. The Journal of Computational Finance 8(2), 55–64 (http://doi.org/bqzn).

Forsyth, P. A., Vetzal, K. R., and Zvan, R. (2002). Convergence of numerical methods for valuing path-dependent options using interpolation. Review of Derivatives Research 5(3), 273–314 (http://doi.org/fjwh92).

Fu, M. C., Madan, D. B., and Wang, T. (1999). Pricing continuous Asian options: a comparison of Monte Carlo and Laplace transform inversion methods. The Journal of Computational Finance 2(2), 49–74 (http://doi.org/bqzp).

Fusai, G. (2004). Pricing Asian options via Fourier and Laplace transforms. The Journal of Computational Finance 7(3), 87–106 (http://doi.org/bqzq).

Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering, Volume 53. Springer.

Ho, A. C. T. (2000). Willow tree. PhD Thesis, University of British Columbia.

Hsu, W. W.-Y., and Lyuu, Y.-D. (2007). A convergent quadratic-time lattice algorithm for pricing European-style Asian options. Applied Mathematics and Computation 189(2), 1099–1123 (http://doi.org/dm4sgg).

Hull, J. C., and White, A. D. (1993). Efficient procedures for valuing European and American path-dependent options. Journal of Derivatives 1(1), 21–31 (http://doi.org/ctjp93).

Ju, N. (2002). Pricing Asian and basket options via Taylor expansion. The Journal of Computational Finance 5(3), 79–103 (http://doi.org/bqzr).

Ritchken, P., Sankarasubramanian, L., and Vijh, A. M. (1993). The valuation of path dependent contracts on the average. Management Science 39(10), 1202–1213 (http://doi.org/chw7hr).

Večeř, J. (2001). A new PDE approach for pricing arithmetic average Asian options. The Journal of Computational Finance 4(4), 105–113 (http://doi.org/bqzs).

Xu, W., and Yin, Y. (2014). Pricing American options by willow tree method under jump–diffusion process. Journal of Derivatives 22(1), 46–56 (http://doi.org/bqzt).

Xu, W., Hong, Z., and Qin, C. (2013). A new sampling strategy willow tree method with application to path-dependent option pricing. Quantitative Finance 13(6), 861–872 (http://doi.org/bqzv).

Zhang, J. (2001). A semi-analytical method for pricing and hedging continuously sampled arithmetic average rate options. The Journal of Computational Finance 5(1), 59–80 (http://doi.org/bqzw).

Zhang, J. (2003). Pricing continuously sampled Asian options with perturbation method. Journal of Futures Markets 23(6), 535–560 (http://doi.org/dxztxg).

	$\displaystyle\\|E^{n}\\|$	$\displaystyle\leq\mathrm{e}^{-r\Delta t}\sum_{j=1}^{m}p_{ij}^{n}[\alpha_{jk}^{% n+1}+(1-\alpha_{jk}^{n+1})]\\|E^{n+1}\\|+\text{i-err}+\text{t-err}$
		$\displaystyle\leq\\|E^{n+1}\\|+\text{i-err}+\text{t-err}.$		(4.12)

Journal of Computational Finance

An efficient convergent lattice method for Asian option pricing with superlinear complexity

Ling Lu, Wei Xu and Zhehui Qian

Need to know

Abstract

Introduction

2 Willow tree for Brownian motion

3 Willow Tree Method for Asian Option Pricing

3.1 Basic algorithm

Algorithm 3.1.

3.2 Complexity-reduction technique

Theorem 3.2 (Easy price theorem).

Proof.

4 Convergence Analysis

Theorem 4.1.

Proof.

5 Numerical Experiments

6 Conclusion

Declaration of interest

Acknowledgements

References

All issues

Browse journals

You are currently on corporate access.