Asian Options Pricing
A Comparative Study
Po-Yuan Lin
Master’s Thesis to obtain the degree in Actuarial Science and Mathematical Finance University of Amsterdam
Faculty of Economics and Business Amsterdam School of Economics
Author: Po-Yuan Lin
Student nr: 10826254
Email: noel.lin@utoronto.ca
Date: August 13th, 2015
Supervisor: dhr.prof.dr.ir.M.H.(Michel) Vellekoop Second reader: dhr. dr. L.J. (Leendert) van Gastel
Contents
1 Introduction 6
2 Background 8
2.1 Assumptions . . . 8
2.2 Standard Options and Asian Options . . . 9
2.2.1 Standard Options . . . 9
2.2.2 Asian Options . . . 10
2.3 Literature Review . . . 11
2.3.1 Analytical Approximation Formulas . . . 11
2.3.2 Finite Difference Method . . . 12
2.3.3 Lattice Method . . . 12
2.3.4 Monte-Carlo Method . . . 13
3 Analytical Approximation Formula 14 3.1 Turnbull and Wakeman’s Formula . . . 14
3.2 Milevsky and Posner’s Formula . . . 16
4 Asian Options via Classical CRR Tree 18 4.1 Construction of the Cox-Ross-Rubinstein Binomial Tree . . . 18
4.2 Hull-White Algorithm . . . 20
4.2.1 Illustration of the Hull-White Algorithm . . . 22
4.2.2 The Log-Spacing h . . . 24
4.3 Costabile-Massabo-Russo Algorithm . . . 25
4.3.1 Illustration of the Costabile-Massabo-Russo Algorithm . . . 26
5 Asian Options via Willow Tree 27 5.1 Willow Tree . . . 27
5.2 Asian Option Pricing . . . 30
5.2.1 Illustration of the Willow Tree Method . . . 31
5.2.2 Modified Hull-White Algorithm on the Willow Tree . . . 32
6 Numerical Results 34 6.1 European-Styled Asian Options . . . 34
7 Conclusion 52
Abstract
An Asian option is a path-dependent derivative that has served as a popular medium for large institutions because of hedging for its lower price, better stability, and reduction on price manip-ulation risk. For arithmetic averaging Asian options, there is no closed formula for pricing and an approximation is required. Though there are approximation formulas for approximating the price of arithmetic averaging Asian options, these formulas show generally a lack of error control and flexibilities. On the other hand, lattice based methods have been shown to be accurate, flexible, and more time efficient than Monte-Carlo methods. Therefore, in this thesis, we compare the algorithms proposed by Hull and White (1993), Costabile et al. (2006), and an algorithm that is analogous to Hull and White (1993) on the willow tree proposed by Curran (1998) using the sampling method from Xu (2013). Specifically for European-Styled arithmetic averaging Asian options, we also im-plement the analytical approximation formulas proposed by Turnbull and Wakeman (1991) and Milevsky and Posner (1998). We examine the computation time, accuracy, and convergence speed of these pricing methods using the Least-Square Monte-Carlo method proposed by Longstaff and Schwartz (2001) as our benchmark. We also apply Richardson extrapolation to improve accuracy and convergence speed. After a through investigation, we have found that the algorithm based on the willow tree has a better performance for accuracy and convergence speed under deep out-of-money and slightly out-of-money scenarios and the Hull and White algorithm has a better performance for accuracy and convergence speed under at-the-money, slightly in-the-money. and deep in-the-money scenarios. When Richardson extrapolation is applied, the Hull and White algorithm outperforms other approaches in terms of accuracy and convergence speed for all scenarios of moneyness. For American-styled Asian options, the algorithm based on the willow tree is better than other lattice based algorithms for most scenarios of moneyness. When Richardson extrapolation is applied, the convergence speed improves at earlier time steps for each algorithm although the accuracy is worse at later time steps. In addition, Richardson extrapolation has been shown to reduce the accumu-lated interpolation errors embedded in the lattice based algorithms. Finally, we have showed that the approximation formulas from Turnbull and Wakeman (1991) and Milevsky and Posner (1998) lack of error control and are poor in terms of accuracy in general. We use numerical results to support our findings.
Keywords: Asian Options, Hull and White, Lattice Methods, Willow Tree, Monte-Carlo Methods, Richardson Extrapolation
1
Introduction
Originated in Tokyo, 1987, the idea of options whose pricing formulas depend on the average prices of the underlying assets was first introduced to crude oil. Since these options first appeared in Asia, they were called ”Asian options”. Over time, Asian options have expanded to other commodities, currencies, interest rates, indices or equities; served as a popular medium for hedging for a better stability and relatively lower price than other standard options. Asian options allow the holders to enjoy a promising reduction on the risk of price manipulation or price distortion due to possible economic events due to the averaging effect. Although American-styled Asian options still expose the holders to potential price manipulation, it is much less pronounced.
Despite of the benefits of Asian options, they are still considered quite illiquid and difficult to price due to their embedded complicated path-dependent structure. The first shortcoming has improved progressively as there are more and more standardized Asian option products in sub-market such as NYMEX and hence quoted prices could be obtained. For the second shortcoming, there are several pricing methods which aim to overcome the pricing difficulty of Asian options and produce a reasonable and efficient performance. Roughly, these methods can be categorized by four approaches: analytical approximations, finite difference methods, lattice methods, and Monte-Carlo methods. They will be discussed in Section 2 in more detail.
Generally, the averaging methods used by Asian options are in either geometric or arithmetic form. For the first case, there is a closed formula based in the classical Black-Scholes model for European-styled Asian options as the geometric average of log-normal random variables is still following a log-normal distribution.This is discussed in Zhang (1998) in more detail. However, for the second case, this is indeed not true and there is no closed formula for this case. This is because, although we can find explicit expressions for all moments of the arithmetic averaging process, there is no explicit expression for the associated distribution function. Consequently, we cannot obtain an exact pricing formula for Asian options and approximation must be performed.
There are several approaches attacking the pricing difficulty of arithmetic averaging Asian options. For example, Levy (1992) uses the geometric averages to approximate the Asian option with an arithmetic average, Turnbull and Wakeman (1991) applies Edgeworth expansion on the probability distribution of the arithmetic average, Kemna and Vorst (1990) approximates the solution of the stochastic differential equation of Asian option price function, Milevsky and Posner (1998) approximates the sum of the asset prices by reciprocal gamma distribution, etc. Although all of these mentioned approaches provide explicit approximation formulas to value Asian options, these methods are not applicable in many cases and suffer a serious drawback of error control.
There are also studies using finite difference methods to price arithmetic averaging Asian options by solving the partial differential equation of the option price function. However, this approach is more difficult to implement and its performance is rather not stable and satisfying. For example, explicit finite difference method can be unstable for path-dependent type of options and implicit finite difference method is only suitable for particular setting of volatility. Although, Vecer (2001) has proposed a new approach to overcome the instability we just mentioned, this is only applicable for European-styled Asian options.
On the other hand, without the mathematical difficulty as in finite difference methods, the lattice based approach is much easier to be implemented and enjoys much more flexibilities on the types of Asian options (American or European) or characteristics of the underlying asset. Hull and White (1993) propose the first lattice based algorithm to approximate the price of Asian option by selecting representative averages of each node and at each time step and approximate the option values of those
averages by using linear interpolation. Although, the Hull and White algorithm (HW) is considered easy to implement, it suffers heavily from the accumulated interpolation errors from time step to time step. Although this could be eased by choosing a smaller spacing h between representative averages in log direction, this consequently increases the computation time of the algorithm significantly. Unfor-tunately, there is still no clear rule for choosing an optimal log-spacing h. Forsyth (2002) and Klassen (2001) propose a log-spacing h as a function of volatility, step size, maturity, and a scaling parameter, but there is still no instruction on picking an optimal scaling parameter. To reduce interpolations error embedded in the Hull-White algorithm, many methods have been proposed and most of these methods are modifications of the Hull-White algorithm by using other mechanisms on choosing the representative averages. Some examples are, Barraquand and Pudet (1996), Klassen (2001), Chalasani et al. (1998), Costabile(2006) et al, and Dai et al (2007).
Curran (1998) proposes the idea of the willow tree, and Ho (2000) showed the accuracy and efficiency of the willow tree on standard option pricing in comparison with other popular binomial tree methods. Later, Xu et al (2012), tested the willow tree on Asian option pricing in comparison with the Monte-Carlo method and binomial tree method using their proposed sampling method and showed that the willow tree method outperformed others in terms of speed and accuracy. Unlike the Hull-White algorithm or its modifications, Asian option pricing methods based on willow tree can track all the possible averages of each time step and thus a theoretical Asian option price should be obtainable. However, as the number of time steps increases, the number of averages increases significantly which causes the efficiency to vanish. To overcome this, I use the idea of Hull-White’s algorithm on willow tree to select representatives at each time step. This is discussed Section 5.
In this thesis, we would like to examine the performance of the Hull and White algorithm, Costabile’s algorithm, and the willow tree method on pricing American and European Asian options. We will examine and analysis their performance in terms of computation time, accuracy, and convergence speed under different scenarios of moneyness. To benchmark appropriately and consistently, we apply the Least-Square Monte-Carlo method for both European-styled and American-styled Asian options and calculate the ”correct” prices under different scenarios of moneyness using 500,000 simulation trials and 500 time steps. We also apply the Richardson extrapolation on these lattice based algorithms and examine and analysis the corresponding effects on improving accuracy and convergence speed. Finally, we implement the approximation formulas proposed by Turnbull and Wakeman and Milevsky and Posner as alternative approaches of pricing European-styled Asian options and compare and analysis their performances with the lattice based algorithms.
2
Background
2.1 Assumptions
Throughout this thesis, we will adopt the assumptions taken in the derivation of the Black-Scholes-Merton model as discussed in the book ”Options, Futures, and Other Derivatives” by John.C.Hull. These assumptions are listed as below:
1. The asset price follows a geometric Brownian motion with constant parameters µ and σ. 2. The short selling of securities with full use of proceeds is permitted.
3. No transactions costs or taxes. 4. All securities are perfectly divisible.
5. No dividends during the life of the derivatives. 6. No arbitrage opportunities.
7. Securities trading is continuous.
8. The risk-free rate r is constant and the same for all maturities. Immediately from assumption 1, we know the asset price follows:
dSt= µStdt + σStdWt (1)
where Wt is a standard Brownian motion, and the drift parameter µ and volatility σ are constant.
The expression (1) is under the market measure P and is not a martingale. As a core of mathematical finance, we will look at the discounted asset price ˜St= Ste−rt and translate the market measure P to
the equivalent martingale measure or so-called risk-neutral probability measure Q. This risk-neutral measure is unique when the market is complete.
To transform the P-measure to the Q-measure, we first note that the stochastic differential equation of the discounted asset price is
d ˜St= (µ − r) ˜Stdt + σ ˜StdWt (2)
By applying the famous Girsanov’s theorem and choosing a constant parameter θ = µ−rσ , we can define a Q-Brownian motion process
n WQ t o such that WQ t = Wt+ tθ
with the Randon-Nikodym derivative defined as dQ dP|Ft = Lt where Lt = exp −Rt 0θdWs− 1 2 Rt 0 θ2ds
and Ft is the filtration that records all the information on
Immediately, we can see that d ˜St = σ ˜StdWtQ and thus a martingale. Then, denote CT to be the
claim value which depends on ST at time T. By the Fundamental Theorem of Asset Pricing in the
Black-Scholes framework, the corresponding option value of this claim at time 0 is then C0= EQe−rTCT
(3) given that EQ[C2
T] < ∞ and CT ∈ FT; T is the maturity time of this claim.
It can be shown that claim value function Ct must satisfy the following partial differential equation
∂C ∂t + 1 2σ 2S2∂2C ∂S2 + r S∂C ∂S − C = 0 (4)
and the solution of this differential equation can be found using the Feynman-Kac representation. In fact, equation (4) is the famous Black-Scholes equation.
Having constructed the Black-Scholes framework of pricing a certain claim value, we can then explore different types of claims and find their claim values. In particular, we will be interested in Asian options. This is discussed in the next section.
2.2 Standard Options and Asian Options
2.2.1 Standard Options
An option is an agreement between two parties that gives its holder the right to buy or sell the underlying asset for a particular price K (strike price) at or before a particular time T (maturity date). The underlying asset varies. It could be an interest rate, an index, a foreign currency or a stock price. In this thesis, we focus on stock price.
If the option can only be exercised at T, then it is an European-styled option. On the other hand, if the option can be exercised any time before the expiration, it is an American-styled option. Particularly, if the option allows us to buy the underlying asset at/before T, it is a call option. On the other hand, if the option allows us to sell the underlying asset at/before T, it is a put option. Specifically, let Ct
and Pt be the option values of a call option and a put option, respectively. Then the corresponding
payoffs at time T are
CT = max (ST − K, 0)
PT = max (K − ST, 0)
If we insert CT or PT in (3) and apply the PDE (4), we then obtain the so-called Black-Scholes pricing
formulas for the European-styled call option and put option.
C0 = S0N (d1) − Ke−rTN (d2) P0 = Ke−rTN (−d2) − S0N (−d1) where d1 = lnS0 K + r +σ22 T σ√T d2 = d1− σ √ T
and N (·) is the cumulative probability distribution function for a standardized normal distribution. Unfortunately, for American-styled call and put options, there is no explicit pricing formula and we must approximate. This can be done by using a finite difference method to solve the PDE of the option price function iteratively. However, this is relatively more difficult to implement and less flexible. An easier alternative is the lattice method which discretizes the underlying price process and approximates the option prices through lattices. Lattice methods are relatively easier to implement and enjoy more flexibility as they can be used to price many popular types of option. One of the famous ones is the Cox-Ross-Rubinstein recombinant binomial tree (CRR) for which they have showed when the number of steps N in the tree gets larger, the resulting price will converge to the Black-Scholes option price for European-style options. See Cox-Ross-Rubinstein (1979) for more details. We describe the construction of the CRR tree in Section 4.
In addition to the above mentioned methods, we can also use the Monte-Carlo method to simulate the evolution of the underlying asset price and obtain a ”correct” price when the numbers of time steps and simulation trials are both sufficiently large. Despite of the accuracy of the Monte-Carlo method, it may be less efficient in comparison to the above mentioned methods. We usually use the Monte-Carlo method as a benchmark to examine the accuracy and efficiency of other pricing methods. In this thesis, we use the Least-Square Monte-Carlo (LSMC) method as a benchmark for pricing European and American Asian options.
2.2.2 Asian Options
Unlike the standard options we discussed in the previous section, the payoff of an Asian option depends on the path that the underlying asset price takes. Specifically, there are primarily four types of Asian options:
• Average Price Call • Average Price Put • Average Strike Call • Average Strike Put
Like standard options, we also have American-styled and European-styled versions of Asian options. Correspondingly, the associated payoffs of above four types of European-styled Asian options at time T are:
• max (0, A (T ) − K) • max (0, K − A (T )) • max (0, ST − A (T )) • max (0, A (T ) − ST)
where A (T ) is the average of the underlying asset prices at time T such that A (T ) =
RT 0 Sudu
T . In this
thesis, the Asian option we refer to will only be the Asian average price call option unless specified otherwise. The price of the Asian option at time 0 is then
V0= EQ e−rT 1 T Z T 0 Stdt − K + (5) where X+= max (X, 0).
Equation (5) is difficult to evaluate because, although we could find explicit formulas for all moments of the averaging process T1 R0T Stdt, we cannot find an explicit expression for the associated distribution
function. As a consequence, there is no explicit formula for the Asian option price. Although we can find an explicit solution to a multidimensional version of Black-Scholes equation, it is much more difficult to solve than the one-dimensional version. Instead, we usually approximate it by discretizing the underlying asset prices. Naturally, we approximate the averaging process by either an arithmetic average N +11 PN j=0Sj or geometric average QN j=0Sj N +11
. For the later case, the average follows a log-normal distribution and the European-style Asian option price could be expressed explicitly as follow: V0= e−rT S0eA+ B2 2 N ln S0 K + A B + B ! − KN ln S0 K + A B !! (6) where A = r − σ 2 2 T 2 B = s 2N + 1 6 (N + 1)σ √ T
When N → ∞, (6) converges to the continuous case and is as expressed below:
V0= S0e −r+σ26 T2 N lnS0 K + r + σ62 T 2 σ q T 3 − Ke −rT N lnS0 K + r − σ22 T 2 σ q T 3
This is discussed more in details in Zhang (1998) and Kemna and Vorst (1990).
On the other hand, for the arithmetic average, we cannot find an explicit expression of the option price and we must approximate it using one of the following approaches
1. Analytical Approximation Formula 2. Finite Difference Method
3. Lattice Method 4. Monte-Carlo Method
2.3 Literature Review
2.3.1 Analytical Approximation Formulas
Wakeman (1991) uses Edgeworth expansion to approximate the probability distribution of arithmetic average; Kemna and Vorst (1990) approximates the solution of the stochastic differential equation of the Asian option price function; Milevsky and Posner (1998) approximates the arithmetic averages using the reciprocal gamma distribution; Benhamou (2002) enhances and extends the pricing algorithm proposed by Carverhill and Clewlow (1992) which was based on the Fast Fourier Transformation in the Black-Scholes framework to price the Asian option with non log-normal densities; Geman and Yor (1993) uses the Bessel process to derive a formula for the Laplace transform of the Asian option. Although we can have explicit formulas for the approximation of the price of the Asian option, these formulas suffer heavily from lack of error control and are not applicable for the American-styled Asian option.
In this thesis, we will use the pricing formulas proposed by Turnbull and Wakeman (TW) and Milevsky and Posner (MP) and compare with other algorithms based on the lattice methods.
2.3.2 Finite Difference Method
Finite difference methods solve the stochastic differential equation of the Asian option price iteratively by converting the differential equation into a set of differential equations and place a grid of points over the time points and stock price points. Although the finite difference methods can handle both American-styled and European-styled Asian options but the results are rather not satisfying. The explicit finite difference method is rather unstable when it is applied to PDE from a path-dependent derivative and the implicit finite difference method is stable and accurate only for certain volatility structures. Vecer (2001) proposes a new PDE approach to fix the above mentioned issues, but it is applicable only to European-style Asian option. These are discussed in more details in Wilmott et al. (1993).
2.3.3 Lattice Method
Lattice methods provide a relatively simple and flexible approach to price the Asian option. They can be applied for both American-styled and European-styled Asian options. However, the evolution of the arithmetic average of the underlying asset prices is not analytically traceable. Hull and White (1993) proposes the first algorithm based on the Cox-Ross-Rubinstein Recombinant Binomial Tree to price the Asian option which avoids tracking all recorded arithmetic averages. It selects a set of representative arithmetic averages for each time step and assign it to all the nodes at that time step. The algorithm then applies piecewise linear interpolation to approximate the option values of those averages that are left out. Though this algorithm is efficient and easy, it suffers a serious problem of errors control due to the interpolation errors accumulated at each time step. Though they claim the interpolation errors will vanish asymptotically when the log-spacing between each representative average converges to 0, the efficiency of the algorithm is reduced gradually.
To overcome the interpolation errors from the Hull-White algorithm (HW), Barraquand and Pudet (1996), Klassen (2001), Chalasani et al. (1998), Costabile et al. (2006), and Dai et al (2007) all provide similar approaches of choosing the representative arithmetic averages based on the modifications of the Hull-White algorithm and they all claim the interpolation errors are reduced significantly.
Curran (1998) proposes the first idea of using the willow tree method. Ho (2000) showed that the willow tree method outperformed other lattice methods when pricing standard options. Unlike other lattice methods, the willow tree provides an advantage of making all the arithmetic averages traceable at each time step. Xu et al (2013) showed that the willow tree to be more efficient than the Hull-White algorithm when N is large. Xu et al use the number of nodes at each time step to be 30 and use 50
time steps. This setting will produce 3049 arithmetic averages at maturity which will far exceed the limitations allowed by MATLAB. Instead of using all the arithmetic averages, we apply the idea of the Hull-White algorithm on the willow tree to find the required representative averages.
In this thesis, we will consider the algorithms proposed by Hull and White and Costabile et al. In addition, we will also consider the willow tree method with modified Hull-White algorithm.
2.3.4 Monte-Carlo Method
The Monte-Carlo method provides a probabilistic approach to estimate the price of the Asian option by sampling the underlying asset prices at each time step. To ensure accuracy, the number of simulation trials must be large enough. However, this also reduces the efficiency gradually and increases the vari-ance due to the randomness embedded when performing simulation. To overcome these issues, Boyle, Broadie, and Glasserman (1996), Broadie and Glasserman (1997), and Lapeyre and Temam (2001) use the variance reduction techniques to achieve a better efficiency and accuracy. These techniques are antithetic variates, moment matching, Latin hypercube sampling, and control variates.
However, above mentioned techniques are not so satisfying to come out a good early exercise rule when dealing with American-styled options. Longstaff and Schwartz (2001) overcome this difficulty by proposing an algorithm using least squares to estimate the conditional payoff to the option holder from continuation. This algorithm greatly reduces the complexities when pricing path-dependent derivatives and decreases the long computation time that other Monte-Carlo methods would experience.
In this thesis, we will use the Least-Square Monte-Carlo method (LSMC) proposed by Longstaff and Schwartz as our benchmark. The correct price is calculated using 500,000 simulation trials and 500 time steps using LSMC.
3
Analytical Approximation Formula
3.1 Turnbull and Wakeman’s Formula
As we have already discussed, to evaluate equation (5), we will have to know the distribution of the averaging process A (T ). Unfortunately, this is not known in closed form and we must approximate the associated probability density function. Turnbull and Wakeman (1991) apply the Edgeworth expansion method to approximate the true density function under the assumption that the averaging process follows a log-normal distribution approximately (Mitchell, 1968). The Edgeworth expansion method approximates the true distribution F by an approximating distribution G such that
f (x) = g(x) +k2(F ) − k2(G) 2! d2g(x) dx2 − k3(F ) − k3(G) 3! d3g(x) dx3 +k4(F ) − k4(G) + 3 (k2(F ) − k2(G)) 2 4! d4g(x) dx4 + (x)
where f and g are the density functions of F and G, respectively; ki(F ) and ki(G) are the ith cumulants
of F and F, for i = 1, 2, 3, 4. (x) is the residual error.
Therefore, to apply the Edgeworth expansion, we must find the first four cumulants of F and G. As mentioned already, Turnbull and Wakeman assumes that the true distribution of the averaging process follows a log-normal distribution F and approximates it by an alternative log-normal distribution G. Suppose we are interested in the arithmetic averaging process
AN = 1 N N X i=1 Si Define h = T N and ti = ih, f or i = 0, 1, . . . , N We also define the relative stock price from time ti−1 to ti:
Ri= Si Si−1 , f or i = 1, 2, . . . , N and LN +1= 1 and Li = 1 + RiLi+1, f or i = 2, 3, . . . , N (7) By (2), logRi ∼ N r −1 2σ 2 h, σ2h f or i = 1, 2, . . . , N (8)
After a little work, we can express Si and AN as:
Si= S1R2. . . Ri and AN =
S0
NR1L2, f or i = 2, 3, . . . , N Define
Y = R1L2
It has been shown by Turnbull and Wakeman that the Asian option will only be exercised if Y > k, where k = N
S0
Therefore, to apply the Edgeworth expansion, we must find the cumulants of Y via calculating the moments of Y. This is readily computable using (7) and (8) since we know the mth moment of Y is just
E (Ym) = E ((R1L2)m)
= E (Rm1 ) E (Lm2 )
by the independence of the increments of a Brownian motion. Specifically, the first four cumulants are k1(F ) = E(Y )
k2(F ) = E(Y2) − E(Y )2
k3(F ) = E(Y3) + 2E(Y )3− 3E(Y2)E(Y )
k4(F ) = E(Y4) − 6E(Y )4− 4E(Y3)E(Y ) + 12E(Y2)E(Y )2− 3E(Y2)2
The moments of G are easily computable. Suppose G has parameters ν and λ, the mth moment of G is
αm(G) = eνm+
1 2λ
2m2
where αm(G) is the mth moment of G. Turnbull and Wakeman match the first two moments of G and
F to derive ν and λ, i.e.,
α1(G) = α1(F ) and α2(G) = α2(F )
After some work, the solution to the above equations is ν = lnα1(F ) −
1
2lnα2(F ) λ = lnα2(F ) − 2lnα1(F )
and other moments of G can immediately be calculated.
Having calculated the first four cumulants of F and G, the Edgeworth expansion method can be applied. Turnbull and Wakeman show that
E [(Y − k)+)] = eν+12λ 2 N (d1) − kN (d2) +k2(F ) − k2(G) 2! d2g(k) dx2 − k3(F ) − k3(G) 3! d3g(k) dx3 +k4(F ) − k4(G) + 3 (k2(F ) − k2(G)) 2 4! d4g(k) dx4 (10) where d1 = ν + λ − lnk λ and d2 = d1− λ and dg(x) dx = g(x) 1 x ν − lnx λ2 − 1 d2g(x) dx2 = g(x) 1 x2 ν − lnx λ2 − 1 ν − lnx λ2 − 2 − 1 λ2
Since Turnbull and Wakeman match the first two moments of F and G, the second term on the right hand side of (10) is simply 0. Finally, we know the option will be exercised if (9) is true.
Therefore, V0= e−rTE (AN − K)+ = e−rTS0 NE (Y − k)+
and the price of an Asian option at time 0 is then obtained using (10).
Note: if we only use the first two moments of the averaging process, we can apply the Black-Scholes formula for futures and obtain another approximating formula for the European-styled Asian option which is V0 = e−rT (F0N (d1) − KN (d2)) where d1= lnF0 K + σ2 aT 2 σa √ T d2= d1− σa √ T F0= M1 σa2= 1 Tln M2 M2 1
and M1 and M2 are the first two moments of the averaging process such that
M1 = S0 erT − 1 rT M2 = 2S02e(2r+σ2)T (r + σ2) (2r + σ2) T2 + 2S02 rT2 1 2r + σ2 − erT r + σ2
3.2 Milevsky and Posner’s Formula
Unlike Turnbull and Wakeman which assumes that the arithmetic average follows the log-normal distribution, Milevsky and Posner (1998) derive the probability density function of the infinite sum of correlated log-normal random variable and show it follows a reciprocal gamma distribution. They then use this result and the duality of the gamma distribution to approximate the finite arithmetic average and then price the Asian option using the reciprocal gamma distribution as the state-price density function and obtain an explicit formula for the Asian option that is very similar to the Black-Scholes formula.
Milevsky and Posner show that the arithmetic averaging process will converge to the reciprocal gamma distribution when both N and T converge to infinite. They further show, for finite N and T, the arithmetic averaging process is approximately a reciprocal gamma distribution. Therefore, using the first and second moments of the averaging process found in the previous section and the duality of the gamma distribution, Milevsky and Posner show that the European-styled arithmetic averaging Asian option price is approximated by
V0 = e−rTM1G 1 K|α − 1, β − e−rTKG 1 K|α, β where α = 2M2− M 2 1 M2− M12 β = M2− M 2 1 M2M1
and G(·|α, β) is the cumulative density function of the gamma distribution given parameters α and β; M1 and M2 are the first and second moments of the averaging process found by Turnbull and
Wakeman, respectively.
Milevsky and Posner show that when T → ∞ and r < 12σ2, the above approximating formula converges to the true option pricing formula. Therefore, for very long-dated Asian options, we may expect a good performance using their formula. See Milevsky and Posner (1998) for a further discussion.
We see that both Turnbull and Wakeman and Milevsky and Posner provide explicit formulas to approximate the price of a European-styled Asian option. These approximation formulas are considered efficient in computation time and accurate under certain scenarios. However, they lack error control and are not flexible for many situations. For example, these methods cannot be applied to American-styled Asian option. In the following sections, we introduce lattice based methods to price the Asian option that enjoy a better flexibility and errors control.
4
Asian Options via Classical CRR Tree
4.1 Construction of the Cox-Ross-Rubinstein Binomial Tree
Cox-Ross-Rubinstein (1979) propose a lattice model to approximate the stochastic differential equation (1) of the underlying asset price St using a recombining binomial tree (CRR) that partitions the
lifetime of the derivative that we are interested in to N time intervals of the same size ∆t such that the underlying asset price can either move up by u or down by d = u1 under the risk-neutral probability measure Q. Specifically, let T denote the life of the derivative and 0 = t0 < t1 < · · · < tN = T the
partition points such that tj+1− tj = ∆t = NT for j = 0, 1, . . . , N − 1, then
P S tj+1 Stj = u = p and P S tj+1 Stj = d = 1 − p
where (p, u, d) could be found by solving the following equations:
plnu + (1 − p) lnd = µ p (lnu − µ)2+ (1 − p) (lnd − µ)2 = σ2s ud = 1
where µ and σ2s are the mean and variance of the asset price over one time step such that µ = (r − 0.5σ2)∆t and σ2s = σ2∆t, respectively.
The solutions of p, u, and d of the above equations would ensure the asset price process under the CRR tree converges to the continuous-time version process of the asset price with the first two moments matched. It can be shown that the solutions of p, u, and d are such that
u = eσ √ ∆t, d = 1 u, p = er∆t− d u − d
Having u, d, and p, we can then construct the famous CRR tree. Let (i, j) represent the node at time i with i − j down-ticks and j up-ticks with probability jipj(1 − p)i−j of reaching this node and the
asset price at this node is S0ujdi−j. Figure 1 gives us an example of the CRR tree for N = 3 given an
initial asset price of S0.
S3,3 = S0u3 S2,2 = S0u2 S1,1= S0u S3,2= S0u2d S0,0 = S0 S2,1 = S0ud S1,0 = S0d S3,1= S0ud2 S2,0= S0d2 S3,0 = S0d3 p 1 − p
It is extremely easy to use the CRR tree to price the standard options. We will take a European call option to illustrate using the above figure. Suppose T = 3, K is the strike price, and define Vi,j to be
the option value at node (i, j), then the payoff at maturity is just
V3,j = max (S3,j− K, 0) (11)
for j = 0, 1, . . . , 3.
We then proceed to the earlier nodes using the backward induction method to find the corresponding option values at those nodes. Suppose we are at the 2nd time step t2, then the option value at each
node can be calculated as following:
V2,j = e−r(pV3,j+1+ (1 − p) V3,j) (12)
for j = 0, 1, and 2.
Similarly, we can repeat the same process as in (12) for t1 and calculate the option values at nodes
(1, 1) and (1, 0). Finally, we will reach t0 and find the option price that we desire.
Unlike the pricing methods that we have mentioned in the previous sections, we enjoy a lot of flexibility when using the lattice method. For example, if we would like to find the price of an American call option in the above example, we simply compare the option value of early exercising with the option value of holding at each node before the maturity and use the backward induction method to find the option price using the exactly same procedure as we use on European options. Note, for the American call options, (12) becomes that
V2,j = maxS2,j− K, e−r(pV3,j+1+ (1 − p) V3,j)
(13) and we repeat (12) on earlier nodes and ultimately find the option price V0,0 at t0.
To find the price of an arithmetic averaging Asian option using the CRR tree, we look at the discrete averaging process A(T ) =
PN j=0Sj
N +1 where Sj is the asset price at time j. The payoff at maturity is then
VT = max (A (T ) − K, 0)
Therefore, to find the Asian option price, we need to know the value of the option at every average reaching each node at maturity in order to do so. This becomes extremely difficult and computa-tionally unmanageable when N is large as the number of arithmetic averages to be considered grows exponentially with it.
To overcome the problem of using the CRR tree to price an Asian option, Hull and White (1993) propose to use only representative averages at each node instead of tracking all possible averages. We discuss more in detail the Hull and White algorithm and its modified algorithms in the next subsections. In addition, we also discuss another lattice model based on the willow tree originally proposed by Curran (1998). Though every possible arithmetic average is traceable using the willow tree, the total number of averages to be considered grows with mN −1where m is the number of nodes at each time step. This also make it computationally unmanageable. We therefore apply the idea of
4.2 Hull-White Algorithm
As an extension to the CRR tree, Hull and White (1993) propose an efficient algorithm by choosing a predetermined set of representative averages at each time step to govern the evolution of arithmetic averages of the tree . To choose the representative arithmetic averages, Hull and White use the idea similar to the Turnbull and Wakeman approach by assuming that the arithmetic averaging process should follow a log-normal distribution approximately. The representative arithmetic averages at each time step should be in the form of S0emh where m is some integer and h is a predetermined number.
We then use these representative averages to determine the option value at each time step. However, this approach does not cover all the possible arithmetic averages and piecewise linear interpolation has to be used on option values whose arithmetic averages are left out.
Suppose we have a CRR tree with N time steps. The Hull and White algorithm is carried out by the following steps:
1. For a predetermined constant log-spacing h, we calculate the maximum and minimum represen-tative arithmetic averages reached at each time step using forward induction. Define Ajmax and
Ajmin as the maximum and minimum representative arithmetic averages reached at time step j. At time step 0, we only have one possible asset price S0, therefore A0max= A0min = S0. At time
step 1, define k1max and kmin1 as the smallest and largest integers at time step 1 such that A0max+ Smax1 2 ≤ S0e k1 maxh and A 0 min+ Smin1 2 ≥ S0e k1 minh where S1
max and Smin1 are the maximum and minimum asset prices at time step 1, respectively.
In the CRR setting, Smax1 = S1,1 and Smin1 = S1,0. The maximum and minimum representative
arithmetic averages at time step 1 are then defined as: A1min= S0ek 1 minh and A1 max = S0ek 1 maxh
Similarly, at time step 2, we will find kmax2 and k2min such that A1max× 2 + S2 max 3 ≤ S0e k2 maxh and A 1 min× 2 + Smin2 3 ≥ S0e k2 minh
Again, in the CRR setting, Smax2 = S2,2 and Smin2 = S2,0. The maximum and minimum
repre-sentative arithmetic averages at time step 2 are defined as: A2min= S0ek 2 minh and A2 max = S0ek 2 maxh
We will then proceed to the later time steps using the exact same procedure and find the pair {(kminj , kjmax)|j = 3, . . . , N } to construct the maximum and minimum representative arithmetic
averages from time step 3 to N.
2. Suppose we have successfully found {(kminj , kmaxj )|j = 0, . . . , N } from the previous step, the
representative arithmetic average at each time step j can be expressed as Ajk = S0ekh, f or k ∈ {kjmin, k j min+ 1, . . . , k j max} where Aj kjmin = A j min and A j kjmax = A j max.
3. For each node at each time step j, we assign the representative arithmetic averages Ajkfound in the previous step to each node, i.e., each node at the same time step will have the same set of representative averages.
4. Calculate the payoff for each representative arithmetic average of each node at maturity. Payoffk= max 0, ANk − K , f or k ∈ {kminN , kminN + 1, . . . , kmaxN }
where ANk is the representative arithmetic average at time step N and K is the strike price. Note: for nodes at maturity, they have the same set of representative averages and the same set of payoff values.
5. Calculate the option value of each representative arithmetic average of each node at time step N − 1. Suppose (N − 1, i) is some node at time step N-1 and AN −1 is some representative arithmetic average of node (N − 1, i) . We know the two possible arithmetic averages ANup and ANdown at time step N, corresponding to an up or down move of the asset price. Explicitly, ANup and ANdown can be calculated from
ANup= N × A N −1+ uS N −1,i N + 1 and A N down = N × AN −1+ dSN −1,i N + 1
ANup and ANdown belong to node (N, i+1) and (N, i), respectively. The associated option values of ANup and ANdown are then found by the interpolation method using the corresponding sets of representative arithmetic averages and the options values at node (N, i+1) and (N, i), respec-tively. For example, suppose ANk is the smallest representative arithmetic average of node (N, i+1) that is greater than ANup and ANk−1 is the largest representative arithmetic average of node (N, i+1) that is smaller than ANup for some integer k ∈ (kminN , kmaxN ). Suppose VkN and Vk−1N are the corresponding options values of ANk and ANk−1, respectively. The option value associated with ANup is then calculated by VupN = A N up− ANk−1 × VkN + ANk − ANup × Vk−1N AN k − ANk−1
Similarily, the option value VdownN associated with ANdown can also be found. The option value associated with AN −1 at node (N-1,i) is then calculated by
VN −1= pVupN + (1 − p) VdownN e−r∆ (14) where p is the probability of an up move, 1-p is the probability of a down move, r is the risk-free rate, and ∆ is the step size of each time step. Once we have successfully found VN −1, we can then apply the exact same procedure to find the associated option values of all the remaining representative averages of node (N-1, i) and for each i.
6. Repeat step 5 for each earlier time step using backward induction method to calculate the option values corresponding to the representative arithmetic averages of each node at each time step. Ultimately, we will reach time step 0 and find the desired option price.
If we are interested in pricing an American Asian arithmetic averaging call option, we will simply repeat the same algorithm as above but with (14) equal to
Vj = max n
pVupj+1+ (1 − p) Vdownj+1
e−r∆, Aj − Ko (15)
when we are calculating the associated option value of each representative arithmetic average at each time step j.
4.2.1 Illustration of the Hull-White Algorithm
The example in this section assumes that S0 = 50, K = 50, T = 1, r = 0.02, σ = 0.3, N = 3, and
h = 0.10. In addition, the risk-neutral probability is p = 0.4790 and the step size is 0.3333.
S3,3= 84.0690 S2,2= 70.6991 S1,1= 59.4555 S3,2= 59.4555 S0 = 50 S2,1= 50.0000 S1,0= 42.0483 S3,1= 42.0483 S2,0= 35.3611 S3,0= 29.7375 p 1 − p
Figure 2: A CRR tree under the assumptions of S0 = 50, K = 50, T = 1, r = 0.02, σ = 0.3, and
N = 3.
Initially at time step 0, the maximum and minimum averages are equal to S0 as we only have one
possible asset price. Therefore, A0max = A0min = 50. At time step 1, the maximum and minimum representative arithmetic averages can be found by finding k1max and kmin1 . Recall, kmax1 and k1min are the smallest and largest integers such that
A0max+ Smax1 2 ≤ S0e k1 maxh and A 0 min+ Smin1 2 ≥ S0e k1 minh
In our example, Smax1 = 59.4555 and Smin1 = 42.0483. After an easy calculation, it can be shown that k1max= 1 and kmin1 = −1. As discussed in step 1 in the Hull and White algorithm, the maximum and minimum representative arithmetic averages are
A1max= S0ek 1 maxh= 55.2585 and A1 min = S0ek 1 minh= 45.2419
Similarly, for j = 2, 3, we can find (kjmin, kjmax) and construct the corresponding minimum and
maxi-mum arithmetic averages Ajmin and Ajmax.
Suppose we have successfully found {(kminj , kmaxj )|j = 0, 1, . . . , 3}, we can the construct the set of
representative arithmetic averages at each time step. Define Aj as the set of representative arithmetic
averages at time step j. After a little work, it can be shown that
A0 = 50, A1 = 55.2585 50.0000 45.2419 , A2 = 61.0701 55.2585 50.0000 45.2419 40.9365 , A3= 67.4929 61.0701 55.2585 50.0000 45.2419 40.9365 37.0409
We will then assign Aj to each node at time step j as its set of representative arithmetic averages. As
a result, each node at the same time step has the same set of representative arithmetic averages. Following step 4, we calculate the payoff of each representative arithmetic average of each node at maturity. Since every node at time step 3 has the same set of representative averages, they also have the same set of payoffs.
V3 = max (A3− K, 0) = 17.4929 11.0701 5.2585 0 0 0 0
Having obtained the payoffs at maturity for each node, we can obtain the option values at earlier time steps following step 5 and 6.
Suppose we are at time step 2 and we are looking at node (2, 1). The set of representative arithmetic averages of node (2,1) is just A2. Suppose we are looking at the representative arithmetic average
55.2585. We know the asset price can either go up to S3,2 or go down to S3,1, then
A3up= 3 × 55.2585 + 59.4555 4 = 56.3077 and A 3 up= 3 × 55.2585 + 42.0483 4 = 51.9560
The associated option values Vup3 and Vdown3 can be found by the interpolation method using A3 and
V3. This is illustrated by using Figure 3.
Following the description in step 5, Vup3 and Vdown3 are calculated by the following
Vup3 = (56.3077 − 55.2585) × 11.0701 + (61.0701 − 56.3077) × 5.2585
61.0701 − 55.2585 = 6.3078
Vdown3 = (51.9560 − 50.0000) × 5.2585 + (55.2585 − 51.9560) × 0
55.2585 − 51.9560 = 1.9560
Therefore, according to (11), the option value for the arithmetic average 55.2585 at node (2,1) is
(p × 6.3078 + (1 − p) × 1.9560) e−r∆t= 4.0008
We will then repeat the above process on all other remaining representative arithmetic averages of node (2,1). We will also repeat the same procedure on all other nodes at time step 2. Having found the option value of each representative arithmetic average of each node at time step 2, we can then move to time step 1 and repeat the exact same procedure to obtain the option values. Ultimately, we will reach time step 0 and obtain the desired option price. The option price is 3.6723 in this example.
Figure 3: At node (2,1) averages of 40.9365, 45.2419, 50.0000, 55.2585, and 61.0701 are considered; at node (3,2) and (3,1), averages of 37.0409, 40.9365, 45.2419, 50.0000, 55.2585, 61.0701, and 67.4929 are considered.
4.2.2 The Log-Spacing h
The option prices calculated by Hull-White algorithm suffer heavily from interpolation errors that accumulate at each time step. Consequently, as N → ∞, there is no solid evidence on whether the price converges to the ”right” value. In fact, the accumulated interpolation error is a function of the log-spacing h. Ideally, as h becomes smaller, the interpolation errors should converge asymptotically to 0 and the accuracy becomes more promising. However, the efficiency becomes much worst due to an increasing number of representative arithmetic averages.
There are several studies on choosing an optimal h for a fixed N. For example, Forsyth et al. (2002) chooses h = α
q
0.25 T σ
2∆ and Klassen (2001) chooses h = βσ√T
1+100N where α and β are parameters
that control the number of representative averages. Unfortunately, there is no defined instruction on choosing the parameters α and β. Barraquand and Pudet (2001) chooses h = √1
N, Forsyth et al (1999)
has provided a numerical evidence against it as the corresponding option price is not reasonable when N is large and suggests we should choose h = N1 as a sufficient condition for convergence. However, when N is small, the result is not so satisfying.
Costabile et al. (2006) introduces a model that is independent of choosing an optimal log-spacing h and uses only ”true” and realized arithmetic averages reached at each node to reduce interpolation errors. This model is discussed in the next section.
4.3 Costabile-Massabo-Russo Algorithm
Like the Hull-White algorithm mentioned in the previous section , the Costabile-Massabo-Russo al-gorithm (CMR) also use representatives averages for each node on the CRR tree. However, unlike the Hull-White algorithm, CMR uses only true and realized arithmetic averages and is independent of the log-spacing h. The intention of the CMR algorithm is to eliminate the interpolation errors and achieve an efficient and accurate performance.
The algorithm simply calculates the maximum and minimum arithmetic averages at each node and uses these values to deduce the other true and realized arithmetic averages for that node. Having done so for each node, step 4∼6 in the Hull-White algorithm are then applied to find the price of the Asian option. Obviously, the need of interpolation decreases, and as the CMR algorithm claims, the interpolation errors should decrease significantly. We then expect the calculated price of the Asian option to be lower than the one under Hull-White algorithm when the log-spacing h is not significantly small.
Since the CMR algorithm is different than the Hull and White algorithm only by the approach they choose the set of representative arithmetic averages at each node, we will then only discuss how the set of representative arithmetic averages is chosen.
1. The first two representative arithmetic averages are the maximum and minimum arithmetic averages reached at each node (i,j). Since for each node (i, j), the maximum arithmetic average reached at (i, j) is contributed by a path starting with j consecutive up movements followed by i - j consecutive down movements; the minimum arithmetic average reached at (i, j) is contributed by a path starting with i - j consecutive down movements followed by j consecutive up movements. As a result, the maximum and minimum arithmetic average reached at node (i, j) can be expressed as: Amax(i, j) = 1 i + 1 j X h=0 S0uh+ i−j−1 X h=0 S0uh+2j−i ! Amin(i, j) = 1 i + 1 i−j X h=0 S0dh+ j−1 X h=0 S0di−2j+h !
Define τmax(i, j) and τmin(i, j) as the paths reached at node (i, j) that produce Amax(i, j) and
Amin(i, j), respectively. In Costabile’s terminology, an asset price path is called a trajectory. We
denote A(i, j; 1) = Amax(i, j) and A(i, j; j(i − 1)) = Amin(i, j) as the first and the last members
of the set of representative arithmetic averages for node (i, j), respectively.
2. Having found the maximum and minimum arithmetic averages reached at each node (i, j) and their associated trajectories, we can then deduce other members of the set of representative arithmetic averages. We will denote these representative arithmetic averages by A(i, j; k) for k = 1, . . . , j(i − j). A(i, j; k) can be computed recursively using the following recursive relationship
A (i, j; k + 1) = A (i, j; k) − 1
i + 1 Smax(i, j; k) − Smax(i, j; k) d
2
(16) where Smax(i, j; k) is the greatest asset price in the trajectory that produces the arithmetic
average A(i, j; k) not belonging to the trajectory τmin(i, j). (16) is saying that the (k + 1)th
representative path is obtained from the previous one by substituting Smax(i, j; k) with the
value Smax(i, j; k) d2. See the illustration in section 4.3.1 for a better understanding of this step.
Having constructed the set of representative arithmetic averages of each node at each time step, we then follow step 4 ∼ 6 to calculate the option price.
Though CMR algorithm does keep track of most of arithmetic averages realized at each node, it does not record all of them. As a consequence, interpolation is still needed for those averages that are left out. This would also produce interpolation errors that may not be significantly smaller than the Hull-White algorithm. In addition, when N is large, number of realized arithmetic averages at nodes at later time steps significantly grow. This gradually reduces the efficiency of this model since this algorithm will spend a lot of time on finding proper candidates to bound those averages that are left out in order to perform interpolation. In fact, when N is large, this algorithm becomes computationally unmanageable.
4.3.1 Illustration of the Costabile-Massabo-Russo Algorithm The following example illustrates step 2 in CMR algorithm.
Consider the CRR tree again,
S4,4= S0u4 S3,3= S0u3 S2,2 = S0u2 S4,3= S0u3d S1,1= S0u S3,2= S0u2d S0 S2,1 = S0ud S4,2= S0u2d2 S1,0= S0d S3,1= S0ud2 S2,0 = S0d2 S4,1= S0ud3 S3,0 = S0d3 S4,0 = S0d4 p 1 −p
Figure 4: CRR tree with 4 steps
Suppose we are looking at node(4, 2), the evolving path that produces Amax(4, 2) or A (4, 2; 1) is
{S0, S0u, S0u2, S0u2d, S0u2d2} or equivalently, {S0, S0u, S0u2, S0u, S0}, therefore, Smax(4, 2, 1) = S0u2,
multiplying it by d2, we get the evolving path that produces A (4, 2; 2) which is {S0, S0u, S0, S0u, S0}.
Similarly, Smax(4, 2, 2) = S0u, multiplying it by d2, we get the evolving path that produces A (4, 2; 3)
which is {S0, S0d, S0, S0u, S0}. Proceed similarly, we have the path {S0, S0d, S0, S0d, S0} for A (4, 2; 4),
and finally {S0, S0d, S0d2, S0d, S0} for A (4, 2; 5) which is just Amin(4, 2). As we mentioned already,
CMR algorithm does not capture all the possible true arithmetic averages. In our case here, the path {S0, S0u, S0, S0d, S0} is not included.
5
Asian Options via Willow Tree
5.1 Willow Tree
The willow tree is a new type of recombining tree algorithm proposed by Curran (1998) to approximate a Brownian motion by a discrete Markov process. Unlike the usual binomial tree algorithms, the willow tree has constant number of tree nodes at each time step and each node can transit to each node at next time step with sets of transition probabilities that is found by solving a sequence of linear programming problems. These sets of transition probabilities are constructed so that the discrete Markov process would converge to a Brownian motion. The original method of constructing the transition probabilities proposed by Curran has a much lower variance and kurtosis than the standard normal distribution when the number of nodes at each time step is small. This causes an underestimation of out-of-money option prices. Later in 2012, Xu et al. (2012), propose a new method to fix this issue by adding more restrictions. In this thesis, we will only consider Xu’s method when pricing the Asian option which is described more with details below.
S1,1 S1,2 S1,3 S1,4
S2,1 S2,2 S2,3 S2,4
S0
S3,1 S3,2 S3,3 S3,4
S4,1 S4,2 S4,3 S4,4
Figure 5: A willow tree with 4 state prices at each time step for 4 time steps.
Originally, Curran proposed the following steps to construct the required Markov process underlying the willow tree:
1. Generate a discrete approximation of the standard normal density function {(zi, qi) |i = 1, . . . , m}
where qi= P (Z ≤ zi). Curran suggested a selection of {zi} and {qi} by
zi = N−1((i − 0.5) /m) (17) qi = 1 m (18) for i = 1, . . . , m. 2. Let Pk = pk ij
be a transition probability matrix for the Markov chain Xk on state space
{1, 2, . . . , m} such that pk
ij = P (Xk+1 = j|Xk= i) for i, j = 1, . . . , m and k = 1, 2, . . . , N − 1
where N is the number of steps in the willow tree. Observe that, when k = 0, p0oi = qi for
i = 1, 2, . . . , m.
3. Let z1, z2, . . . , zmbe the discrete values for the standard normal distribution in the space direction
and let 0 = t0 < t1 < · · · < tn= T be a partition of [0, T ], define Ytk =
√
It has been shown by Ho, 2000, that, in order for Ytk, k = 0, 1, . . . , N to converge to a Brownian motion
as N → ∞, the transition probability matricesPk and q
i= m1, i = 1, 2, . . . , m should satisfy the
following linear programming problem: minpk ij m X i=1 m X j=1 pkij|ptk+1zj− √ tkzi|3 (19) subject to m X j=1 pkij = 1 m X j=1 pkijptk+1zj = √ tkzi m X j=1 pkijtk+1z2j − tkz2i = hk m X i=1 qipkij = qj where 0 ≤ pkij ≤ 1, for i, j = 1, 2, . . . , m.
Above constraints ensure that the mean of the sequence {zi, i = 1, 2, . . . , m} is zero and its variance
is one so that E Ytk+1 − Ytk|Ytk = 0 E Yt2k+1− Yt2 k|Ytk = hk
However, the pair {(zi, qi)|i = 1, 2, . . . , m} sampled from Curran’s sampling method fails to satisfy
the requirement that {zi, i = 1, 2, . . . , m} should have variance of 1. Ho showed that by modifying z1
and zm as z1− and zm+ for some appropriate , this issued could be removed ; the generated
Markov process will converge to a Brownian motion when both m and N approach to ∞. Xu has showed that for a small m, {zi} generated using Ho’s method will have a much lower kurtosis than the
standard normal distribution and suggested to generate (zi, qi) using the following method subject to
new restrictions:
1. For a predetermined γ between 0 and 1, generate qi = (i−0.5)
γ
m for i = 1, 2, . . . , m/2, and then
normalize {qi} by qi= Pmqi
j=1qj where qi = qm+1−i due to symmetry.
2. Generate m intervals as [Zi, Zi+1] where Zi = N−1
Pi j=1qj for i = 1, 2, . . . , m − 1, Z0= −∞ and Zm= ∞.
3. Solve the following constrained nonlinear least-squares problem to obtain zi:
minzi "m X i=1 qiz4i − 3 #2 s.t. m X i=1 qizi = 0 m X i=1 qiz2i = 1 Zi−1≤ zi ≤ Zi
For m = 30 and γ = 0.6, we have sampled the required {(qi, zi)i|i = 1, 2, . . . , m} to be used in this
thesis.
Table 1: (qi, zi) sampled from Xu’s strategy using γ = 0.6 and m = 30. For i > 15, (qi, zi) is omitted
due to the symmetric.
Having generated {(qi, zi)i|i = 1, 2, . . . , m}, we can then revisit the linear programming problem (19)
and find the required transition probability matrices. We should note that, unlike Curran’s sampling method, {qi}i=1,2,...,m sampled from Xu’s strategy is not uniformly distributed. The linear
program-ming problem (19) then becomes
minpk ij m X i=1 m X j=1 qipkij|ptk+1zj − √ tkzi|3 (20)
where all constraints are the same. Particularly, when ∆t is small, the constraint
m
X
j=1
pkijtk+1zj2− tkzi2= hk
can be removed. This is proved by Xu.
Solving the linear programming problem is time costly. However, the transition probability matrices are independent of the characteristics of the asset that we are interested in and only depend on the number of nodes per step and the number of time steps. What is more, the computation could be performed off-line and stored to be used on all kinds of assets. Therefore, when calculating the computation time for pricing Asian option using the willow tree, we will ignore the time for constructing the transition probability matrices.
5.2 Asian Option Pricing
Suppose we have successfully generated {zi, i = 1, . . . , m}, {qi, i = 1, . . . , m}, and transition probability
matrices Pk, k = 0, 1, . . . N − 1 , we can then construct the willow tree as below,
S1,1 S2,1 . . . SN,1 S1,2 S2,2 . . . SN,2 S0 S1,3 S2,3 . . . SN,3 .. . ... ... ... S1,m S2,m . . . SN,m
Figure 6: A willow tree with m state prices each time step for N time steps. where we must take Sij = S0e(r−σ
2/2)t i+σ
√
tizj for i = 1, . . . , N , j = 1, . . . , m. Suppose the transition
probability matrices are q = (q1, q2, . . . , qm), and for i = 1, 2, . . . , N − 1,
Pi= pi 1,1 pi1,2 · · · pi1,m pi2,1 pi2,2 · · · pi2,m .. . ... . .. ... pi m,1 pim,2 · · · pim,m
where pikj represents the transition probability from node k at time ti to node j at time ti+1.
Now let Fi be a m × mi−1matrix recording all the averages reached at each node at time ti such that
F0 = S0 F1 = F0 2 F0 2 .. . F0 2 + S1,1 2 S1,2 2 .. . S1,m 2 and , for i = 2, . . . , N , Fi= reshape(Fi−1) i+1 reshape(Fi−1) i+1 .. . reshape(Fi−1) i+1 + Si,1 i+1 Si,2 i+1 .. . Si,m i+1 ×h1 1 . . . 1 i 1×mi−1
where operation ”reshape” reorganizes the matrix , i.e., reshape(Fi−1) =
h
F1,1i−1. . . Fm,1i−1 F1,2i−1. . . Fm,2i−1. . . F1,mi−1i−1. . . F
i−1 m,mi−1
then the option value at maturity is simply VN = max FN − K. And for i = 1, 2, . . . , N − 1, Vi = e−rT /N
P1iV1i+1 P2iVm+1i+1 . . . Pmi Vmi+1i−m+1
P1iV2i+1 P2iVm+2i+1 . . . Pmi Vmi+1i−m+2
..
. ... . .. ...
P1iVmi+1 P2iV2mi+1 . . . PmiVmi+1i
where Pji are 1 × m matrices that partition the transition probability matrix Pi; Vji+1 are m × 1 matrices that partition the option value matrix Vi+1 at time step i +1. Specifically,
Vi+1= h
V1i+1 V2i+1 . . . Vmi+1i
i
and Pi = h
P1i P2i . . . Pmi i0
Finally, the option value at time 0 is evaluated as V0 = e−rT /N
m
X
j=1
qjVj1
Notice that, unlike the binomial tree, the arithmetic average at each node can be analytically tracked. In addition, the number of nodes grows linearly as m × (N + 1) in the willow tree where as it grows quadratically as 2N +1− 1 in the binomial tree. Therefore, when N is large, number of nodes in the willow tree will be significantly less than the number of nodes in the binomial tree. These suggest a great analytical traceability and efficiency for the willow tree method especially when N is large.
5.2.1 Illustration of the Willow Tree Method Consider the following example assuming m = 2 and N = 3,
S1,1 S2,1 S3,1
S0
S1,2 S2,2 S3,2
Figure 7: A willow tree with 2 state prices each time step for 3 time steps.
Suppose Sij is the underlying asset price on the jth node and time ti and the transition probability
from ti to ti+1 is denoted as
Pi = P i 11 P12i P21i P22i ! , f or i = 1, 2
The matrices that record all the averages reached at each node at time ti are,
F0 = S0 F1 = S0+S11 2 S0+S12 2 ! F2 = S0+S11+S21 3 S0+S12+S21 3 S0+S11+S22 3 S0+S12+S22 3 ! F3 = S0+S11+S21+S31 4 S0+S11+S22+S31 4 S0+S12+S21+S31 4 S0+S12+S22+S31 4 S0+S11+S21+S32 S0+S11+S22+S32 S0+S12+S21+S32 S0+S12+S22+S32 !
Suppose Fijk of Fk represents an arithmetic average of path from t0 to ti, then the payoff at maturity is just V3 = max F3− K, 0 = (F 3 11− K)+ (F123 − K)+ (F133 − K)+ (F143 − K)+ (F213 − K)+ (F223 − K)+ (F233 − K)+ (F243 − K)+ !
Since F113 and F213 at t3 are generated by F112 at t2 with transition probabilities P112 and P122, the
corresponding option value for F2
11 at t2 is then e−r
T 3 P2
11V113 + P122V213. Similarly, the option value
for F122 is just e−rT3 P2
21V133 + P222V233. We then apply the same approach to the remaining paths and
obtain the option value matrix at t2. This is just
V2 = e−rT3 P 2 11V113 + P122V213 P212V133 + P222V233 P2 11V123 + P122V223 P2 21V143 + P222V243 !
Similarly, the option value matrix at t1 is
V1 = e−rT3 P 1 11V112 + P121V212 P1 21V122 + P221V222 !
Finally, the option value at time t0 is then
V0 = e−rT3 q1V1
1 + q2V21
5.2.2 Modified Hull-White Algorithm on the Willow Tree
Since the willow tree method allows us to record all the arithmetic averages reached at each node and at each time step, the option price calculated using the willow tree method should serve as a ”correct” price. However, despite of the willow tree’s analytical traceability, we should observe that, as N grows, the number of arithmetic averages to be considered at each time step grows significantly as mN −1. This is technically unmanageable when either m or N is too large. This issue also occurs in the Costabile-Massabo-Russo algorithm as discussed already. Therefore, instead to use all the arithmetic averages at each node, we will apply the idea from Hull-White algorithm to use representative averages for each node and each time step.
The first four steps of the modified Hull-White algorithm on the willow tree are the same as the first four steps of the Hull-White algorithm on the CRR tree. We therefore only describe the last two steps. 5. Suppose we have successfully found {(kjmin, kjmax)|j = 0, . . . , N } and constructed the set of
representative arithmetic averages for each node at each time step. Suppose we also have found the payoff at maturity VN. Suppose we are at time step N-1 and we are looking at some node (N-1, i). Let AN −1be some representative arithmetic average of node (N-1, i). In the willow tree
setting, there are m possible arithmetic averages AN1 , . . . , ANm at time step N, corresponding to the movements to each state prices SN,1, . . . , SN,m, respectively. Explicitly,
ANs = N × A
N −1+ S N,s
N + 1 , f or s = 1, 2, . . . m
In addition, each ANs belongs to node (N, s), respectively. We therefore interpolate the option value associated with each ANs using the set of representative arithmetic averages and the cor-responding set of option values of node (N, s). For s = 1, 2, . . . , m, suppose for some integer ks ∈ (kNmin, kNmax), ANs,ks is the smallest representative arithmetic average of node (N,s) that is
greater than ANs and ANs,k
s−1 is the largest representative arithmetic average of node (N, s) that
is smaller than AN
s ; suppose Vs,kNs and V
N
s,ks−1 are the corresponding option values of A
N s,ks and
ANs,ks−1, respectively. The option value corresponding to ANs is then calculated by interpolation as VsN = ANs − AN s,ks−1 Vs,kN s+ ANs,k s− A N s Vs,kN s−1 ANs,k s− A N s,ks−1 , f or s = 1, 2, . . . , m (21)
Since the probability of movements from node (N-1, i) to node (N, s) is pN −1i,s , for s = 1, 2, . . . , m, the corresponding option value of AN −1 at node (N-1, i) is
VN −1= e−rT /N
m
X
s=1
pN −1i,s VsN (22)
We will then proceed to other remaining representative arithmetic averages of node (N-1, i) using the above procedure to find all the corresponding option values. We will also apply the exact same procedure to all other remaining nodes at time step N-1. Once we have found all the corresponding option values of all representative arithmetic averages of each node at time step N-1, we will then move to each node at time step N-2 and repeat the exact same procedure. We will repeat the described process on each earlier time step using backward induction method until we reach time step 1.
6. At time step 0, A0 = S
0 is the only representative arithmetic average to be considered. For
s = 1, 2, . . . , m, A1s and Vs1 can be easily computed using (21) and (22) and the corresponding transition probability of movement from S0 to node (1, s) is qs. Therefore, the option price at
time step 0 is V0= e−rT /N m X s=1 qsVs1
In comparison with the Hull-White algorithm on the CRR tree, the number of nodes to be considered are much smaller in the willow tree when N is large, we can expect that the modified Hull-White algorithm on the willow tree enjoys a better computation time than the Hull-White algorithm on the CRR tree.
6
Numerical Results
In this section, we compare and examine the performance of the arithmetic averaging Asian option by the Hull-White algorithm (HW), Costabile-Massabo-Russo (CMR), and the willow tree method using modified Hull-White algorithm (W) under different scenarios of moneynesses KS. The scenarios of moneyness to be considered are : deep out-of-money (KS = 0.75), slightly out-of-money (KS = 0.9), at-the-money (KS = 1), slightly in-the-money (KS = 1.1), and deep in-the-money (KS = 1.25). Specifically, for the European-styled Asian options, we also implement the analytical approximation formulas proposed by Turnbull and Wakeman (TW) and Milevsky and Posner (MP) as comparisons to above mentioned algorithms. In addition, we apply Richardson extrapolation on the above mentioned lattice based algorithms to speed up the convergence speed. In addition, we will ignore the extra computation time of adding the Richardson extrapolation on each mentioned algorithm since the extra computation time of adding the Richardson extrapolation is significantly less than 0.001 seconds on each algorithm. To benchmark, we use the Least-Square Monte-Carlo simulation (LSMC) with 500 thousand simulation trials and 500 time steps as our ”correct” price. The coding is ran on a laptop with Intel(R) Core(TM) i7-4700HQ CPU @2.40GHz with 8.00GB RAM running 64-bit Windows 8.1 and MATLAB R2014a. For all cases discussed in the next subsections, we use the following assumptions: S0 = 50, σ = 0.3,
r = 0.02, and T = 1. Specifically for the willow tree method, γ = 0.6 and m = 30. The ”correct” prices calculated by LSMC under different scenarios of moneyness are shown as below.
Table 2: Option prices calculated by LSMC using 500,000 simulation trials and 500 time steps for both European-styled Asian options and American-styled Asian options under different scenarios of moneyness. The standard deviation and computation time for each scenario of moneyness and each style of Asian option are also given. These prices are treated as our correct prices in this thesis and are used to compute the absolute errors.
6.1 European-Styled Asian Options
For European-styled Asian options, we compare the performance of HW, CMR, W, TW, and MP using LSMC as our benchmark under different time steps and moneynesses. In particular, for algorithm HW and W, we also examine the performance for different log-spacing h. The calculated prices under these mentioned approaches could be found in Table 3 and the absolute errors between these approaches
