The hexanomial lattice for pricing multi-asset options

The hexanomial lattice for pricing multi-asset options

Applied Mathematics and Computation 233 (2014) 463–479 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepag...

800KB Sizes 64 Downloads 212 Views

Applied Mathematics and Computation 233 (2014) 463–479

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The hexanomial lattice for pricing multi-asset options Wen-Hung Kao a, Yuh-Dauh Lyuu b,c,⇑, Kuo-Wei Wen b a b c

Department of Derivatives, Capital Securities Corp. Capital Center, No. 101, Songren Rd., Taipei 11073, Taiwan Department of Computer Science and Information Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan Department of Finance, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan

a r t i c l e

i n f o

Keywords: Binomial lattice Trinomial lattice Hexanomial lattice Multi-asset option Barrier option Correlation

a b s t r a c t Multi-asset options are important financial derivatives. Because closed-form solutions do not exist for most of them, numerical alternatives such as lattice are mandatory. But lattices that require the correlation between assets to be confined to a narrow range will have limited uses. Let qij denote the correlation between assets i and j. This paper defines a (correlation)pffiffiffiffiffiffioptimal lattice pffiffiffiffiffiffi as one that guarantees validity as long as 1 þ Oð Dt Þ 6 qij 6 1  Oð Dt Þ for all pairs of assets i and j, where Dt is the duration of a time period. This paper then proposes the first optimal bivariate lattice (generalizable to higher dimensions), called the hexanomial lattice. This lattice furthermore has the flexibility to handle a barrier on each asset. Experiments confirm its excellent numerical performance compared with alternative lattices. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Options are financial derivatives whose payoff depends on some underlying assets [18]. The underlying assets can be stocks, bonds, currencies, interest rates, volatilities, commodities, temperature, and so on. A fundamental result in finance says the theoretical price of an option equals the discounted expected payoff under the risk-neutral probability measure [28]. Multi-asset options have a payoff determined by multiple underlying assets, and they are essential tools for speculation and risk management. They are also called rainbow options and correlation options [28,33]. These options are challenging to price because of problems arising from their high dimensions and the correlations between assets. Table 1 samples a few popular multi-asset options with their payoff functions specified. Multi-asset options can be priced in many ways. Because of the curse of dimensionality, Monte Carlo simulation and its variants are the most popular general-purpose approaches [5]. Joy et al. [21], Boyle et al. [9], and Tan and Boyle [37] apply quasi Monte Carlo simulation to price multi-asset European and American options. Quasi Monte Carlo methods in general have faster convergence rates and generate better confidence intervals than the standard Monte Carlo method [40]. Simulation is relatively time consuming and does not converge fast, however. Closed-form pricing formulas, in contrast, do not suffer from the same problems if the dimension of the resulting integral is reasonably small. Stulz presents closedform formulas for bivariate maximum and minimum options [36]. Johnson extends them to multi-asset maximum and

⇑ Corresponding author at: Department of Computer Science and Information Engineering, Department of Finance, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan. E-mail addresses: [email protected] (W.-H. Kao), [email protected] (Y.-D. Lyuu), [email protected] (K.-W. Wen). http://dx.doi.org/10.1016/j.amc.2014.01.173 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

464

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

Table 1 A sampling of multi-asset options. Type

Payoff

Exchange option [29] Better-off option [41] Worse-off option [41] Binary maximum option Maximum option [20,36] Minimum option [20,36] Spread option [12,30] Basket average option [20,23] Multi-strike option [41] Pyramid rainbow option [41] Madonna rainbow option [41]

max(S1(T)  S2(T), 0) max(S1(T), . . ., Sn(T)) min(S1(T), . . ., Sn(T)) 1fmaxðS1 ðTÞ;...;Sn ðTÞÞ>Kg max(max(S1(T), . . ., Sn(T))  K, 0) max(min(S1(T), . . ., Sn(T))  K, 0)) max (S1(T)  S2(T)  K, 0) max((S1(T) +    + Sn(T))/n  K, 0) max(S1(T)  K1, . . ., Sn(T)  Kn, 0) max(|S1(T)  K1| +    + |Sn(T)  Kn|  K, 0) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  max ðS1 ðTÞ  K 1 Þ2 þ    þ ðSn ðTÞ  K n Þ2  K; 0 1

Si(T) denotes the price of the ith underlying asset at the maturity date T, Ki is the strike price for the ith underlying asset, and 1 is the indicator function.

minimum options [20]. Kirk and Aron’s [24] approximation method is able to price bivariate spread options, but it is inaccurate when the strike prices are high [2,12,30]. Although closed-form solutions are often preferred, they are rare. To fill this void, partial differential equations (PDEs) and the closely related lattices are popular numerical alternatives [1,6,7,27,38]. Both methods divide the time into discrete time periods with a duration of Dt. There are two general types of finite-difference methods: the explicit and implicit ones The explicit methods of Brennan and Schwartz [4], Boyle and Tian [10], and Hull and White [17] can be used to price bivariate options and they are conditionally stable. (A bivariate option is a multi-asset option with two underlying assets.) Unconditionally stable implicit methods such as Kurpiel and Roncalli’s [25] and Gillia et al.’s [16] are also able to price bivariate options. Although the explicit method is conditionally stable, it is both easier to implement and conceptually simpler than the implicit method. For example, the explicit method for a variety of derivative security pricing problems uses between 40% and 70% as much time as the implicit method to provide the same level of accuracy [17]. Therefore, there does not seem to be a clear winner between the implicit method and the explicit one. A lattice is essentially an explicit finite-difference method to solve the pricing PDE [28]. A lattice for multi-asset options is called a multi-asset lattice; in particular, a lattice for bivariate options is called a bivariate lattice. With few exceptions, multi-asset lattices are either binomial or trinomial. In the former case, each asset’s price may move up or down (hence 2 branches). In the latter case, each asset’s price move has 3 choices (hence 3 branches). The number of jumps of a multi-asset lattice is the total number of branches emanating from each node. Many multi-asset lattices have been proposed. The bivariate lattice of Boyle et al. has 4 jumps as each asset’s price moves follow the binomial branch [7]. The same is true of the bivariate lattice of Rubinstein, who calls it a pyramid [34]. The bivariate lattice of Boyle has 5 jumps as each asset’s price moves follow the trinomial branch [6]. A high-dimensional extension of the above-mentioned lattice of Boyle [6] is the trinomial lattice of Kamrad and Ritchken. It contains 2k + 1 jumps when there are k underlying assets [22]. A lattice should not have invalid transition probabilities [42]. But some multi-asset lattices have invalid transition probabilities unless the correlation between assets falls within a typically narrow range. Let qij denote the correlation between assets i and j. (For brevity, the correlation is simply q for the bivariate case.) To highlight the key role played by qij in the construction of valid lattices, p this ffiffiffiffiffiffi paper calls a lattice pffiffiffiffiffiffi (correlation) optimal if the transition probabilities are guaranteed to be valid as long as 1 þ Oð Dt Þ 6 qij 6 1  Oð DtÞ for all pairs of assets i and j. Note that optimality here concerns correlation not running time. An optimal lattice, therefore, is guaranteed to be valid for essentially all possible correlations. Optimality is clearly a desirable property for lattices. But, surprisingly, this important feature has been overlooked in the literature. This paper analyzes existing multi-asset lattices for their optimality. The 5-jump bivariate lattice of Boyle [6], the 4-jump bivariate lattice of Boyle et al. [7], and the 4-jump bivariate lattice of Rubinstein [34] are shown to be optimal. Lin’s 9-jump bivariate lattice is however suboptimal [26]. This paper proceeds to propose a new lattice, called the hexanomial lattice, that is provably optimal in the bivariate case. In the construction, a time period is divided into two phases. A 4-jump binomial structure is employed in the first phase and a 9-jump trinomial structure in the second. The result is a 36-jump lattice, with 6 branches in each asset. An additional advantage of the hexanomial lattice is its ability to price options with multiple barriers (to be elaborated shortly). Orthogonalization can also be used to construct optimal lattices [11,17,39,41]. It transforms the processes for the assets’ prices into uncorrelated ones. A valid lattice is then built for each process independently. They are later combined to form a multi-asset lattice. Finally, the uncorrelated processes are transformed back to the original assets’ prices. In contrast, the hexanomial lattice works directly with assets’ prices instead of variables which are linear functions of assets’ prices. The connection to the payoff function is thus transparent, whereas orthogonalization obscures it somewhat. A common challenge facing a lattice is pricing barrier options with good convergence behavior. The payoff of a barrier option depends on whether an asset’s price ever touches a certain price level called the barrier. There are two general types of barrier options. A knock-out option immediately terminates if an asset’s price touches the barrier. In contrast, a knock-in

465

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

option comes to existence only if an asset’s price touches the barrier. Although the option prices computed by lattice converge to the correct value as Dt goes to 0 if the probabilities are valid [14], the prices may oscillate significantly [8,13]. A lattice introduces two types of errors: the distribution error and the nonlinearity error [15]. The distribution error arises from approximating a continuous distribution with a discrete distribution; fortunately, it converges to 0 as Dt does. The nonlinearity error, however, is introduced by the nonlinearity of the option value function, which occurs at certain critical prices. For barrier options, the critical price occurs along the barriers. The nonlinearity error for barrier options can be reduced by aligning the lattice nodes with the barriers [31]. Yet none of Boyle’s [6], Boyle et al.’s [7], and Rubinstein’s [34] bivariate lattices are able to achieve that; hence they are likely to have convergence issues when pricing bivariate options with a barrier on each asset. Lin’s bivariate lattice [26] can price bivariate options with a barrier on each asset because it can align the lattice nodes with both barriers, but it is suboptimal. In contrast, our hexanomial lattice not only is optimal but can handle a barrier on each asset. The hexanomial lattice is the first optimal bivariate lattice that can handle two barriers without orthogonalization. For the numerical assessments, the two-asset maximum option, spread option, and multi-strike options in Table 1 are priced. These options have no barriers. To see how the lattices perform when there are barriers, the two-asset multi-strike option with two barriers is also priced. The data show that the hexanomial lattice compares favorably with the bivariate lattices of Boyle [6], Boyle et al. [7], and Rubinstein [34]; the convergence of the hexanomial lattice is smoother and faster. The rest of this paper is organized as follows. Section 2 presents the preliminaries and analyzes many existing multi-asset lattices for their optimality and ability to handle multiple barriers. Section 3 describes the hexanomial lattice, analyzes its complexity, and proves its optimality. The hexanomial lattice is also extended to higher dimensions there. Section 4 describes the numerical results. Section 5 concludes. 2. Preliminaries and analysis of multi-asset lattices Let Si(t) denote the ith asset’s price at time t, where 0 6 t 6 T and T is the maturity date of the option. The lognormal processes for the assets’ prices (assumed to be in a risk-neutral probability measure throughout the paper) are

dSi ¼ rdt þ ri dz; Si

ð1Þ

where r is the risk-free interest rate, ri is the volatility of the ith asset’s price, and the random variable dz is a standard Brownian motion [3]. Eq. (1) has the following unique solution [35]: 2

Si ðtÞ ¼ Si ð0Þeðrri =2Þtþri zðtÞ : A lattice divides the time between 0 and T into n time periods. The duration of one time period is therefore Dt = T/n. Let the strike price and the barrier (if any) for the ith asset be denoted by Ki and Hi, respectively. For brevity, the strike price is simply K for the two-asset case if K1 = K2. Pairs of assets’ prices {S1(t), S2(t)} are assumed to be joint-lognormally distributed. For the ith asset, its log return over one time period [t, t + Dt] is defined as

fi ðt; DtÞ  ln Si ðt þ DtÞ=Si ðtÞ; which is normally distributed with mean liDt and variance r2i Dt, where li  r  r2i =2. Note the additive property: fi(t, ta + tb) = fi(t, ta) + fi(t + ta, tb) for ta, tb > 0. For each pair of assets Si and Sj where i – j, the correlation between fi(t, Dt) and fj(t, Dt) is denoted by qij. It is abbreviated as q for the two-asset case. Let fi ðt; DtÞ denote the lattice’s discrete distribution for fi(t, Dt). Again, the additive property holds: fi ðt; ta þ tb Þ ¼ fi ðt; t a Þ þ fi ðt þ t a ; t b Þ for ta, tb > 0. 2.1. The 4-jump bivariate lattice of Boyle et al. [7] pffiffiffiffiffiffi Let ui  ri Dt , i = 1, 2. For the 4-jump bivariate lattice of Boyle et al. [7], the joint normal random variable {f1(t, Dt), f2(t, Dt)} is approximated by discrete random variable ff1 ðt; DtÞ; f2 ðt; DtÞg defined in Table 2. The solutions of the transition probabilities are:

p1 ¼

  pffiffiffiffiffiffil 1 l ; 1 þ q þ Dt 1 þ 2 4 r1 r2

ð2aÞ

Table 2 The 4-jump bivariate lattice of Boyle et al. [7]. Probability

p1

p2

p3

p4

f1 ðt; DtÞ

u1

u1

u1

u1

f2 ðt; DtÞ

u2

u2

u2

u2

466

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

p2 ¼

  pffiffiffiffiffiffil 1 l ; 1  q þ Dt 1  2 4 r1 r2

ð2bÞ

p3 ¼

  pffiffiffiffiffiffi l 1 l ; 1  q þ Dt  1 þ 2 4 r1 r2

ð2cÞ

p4 ¼

  pffiffiffiffiffiffi l 1 l : 1 þ q þ Dt  1  2 4 r1 r2

ð2dÞ

It can be easily shown that the probabilities are valid when the correlation satisfies

1 þ

    pffiffiffiffiffiffil pffiffiffiffiffiffil l l Dt  1 þ 2  6 q 6 1  Dt 1  2 : r r r r 1

2

1

2

Hence this lattice is optimal. But this lattice lacks the flexibility to align lattice nodes with the barriers. 2.2. The pyramid of Rubinstein [34] pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Let dipdenote the dividend yields of Si, mi  ðlnðr=di Þ  ðr2i =2ÞÞDt, i = 1, 2, u1  r1 Dt, u21  r2 Dt ðq þ 1  q2 Þ, and ffiffiffiffiffiffi p ffiffiffiffiffiffiffiffiffiffiffiffiffiffi u22  r2 Dt ðq  1  q2 Þ. The joint normal random variable {f1(t, Dt), f2(t, Dt)} is approximated by discrete random variable ff1 ðt; DtÞ; f2 ðt; DtÞg defined in Table 3. This lattice is trivially optimal because the transition probabilities are constant and independent of the correlation. But it lacks the flexibility to align lattice nodes with the barriers when pricing barrier options. 2.3. The 5-jump bivariate lattice of Boyle [6] pffiffiffiffiffiffi Let ui  kri Dt , i = 1, 2, and k P 1 be a stretch parameter. For the 5-jump bivariate trinomial lattices of Boyle [6], the joint normal random variable {f1(t, Dt), f2(t, Dt)} is approximated by discrete random variable ff1 ðt; DtÞ; f2 ðt; DtÞg defined in Table 4. The solutions of the transition probabilities are:

pffiffiffiffiffiffi  pffiffiffiffiffiffi  ! ! 1 1 q Dt l1 l2 1 1 q Dt l1 l2 ; p ; þ þ þ ¼  þ  2 4 k2 k2 4 k2 k2 k2 r1 r2 k2 r1 r2 pffiffiffiffiffiffi  pffiffiffiffiffiffi  ! ! 1 1 q Dt l1 l2 1 1 q Dt l1 l2 p3 ¼  þ  þ ¼ þ þ   ; p ; 4 4 k2 k2 4 k2 k2 r1 r2 r1 r2 k2 k2

p1 ¼

p5 ¼ 1 

1 k2

:

The single stretch parameter k can be properly set to align lattice nodes with the barrier. Since one stretch parameter can only handle one barrier, this lattice cannot handle more than one barrier. It can be easily shown that

    pffiffiffiffiffiffil pffiffiffiffiffiffil l l 1 þ k Dt 1 þ 2  6 q 6 1  k Dt  1  2 : r r r r 1

2

1

2

Hence this lattice is optimal. With k ¼ 1, it degenerates to the 4-jump bivariate lattice of Boyle et al. [7] with the same transition probabilities and the same range for q.

Table 3 The pyramid of Rubinstein [34]. Probability

1/4

1/4

1/4

1/4

f1 ðt; DtÞ

m1 + u1

m1 + u1

m1 u1

m1  u 1

f2 ðt; DtÞ

m2 + u21

m2 + u22

m2 u22

m2  u21

Table 4 The 5-jump bivariate lattice of Boyle [6]. Probability

p1

p2

p3

p4

p5

f1 ðt; DtÞ

u1

u1

u1

u1

0

f2 ðt; DtÞ

u2

u2

u2

u2

0

467

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479 Table 5 The 9-jump bivariate lattice of Lin [26]. Probability

p1

p2

p3

p4

p5

p6

p7

p8

p9

f1 ðt; DtÞ

u1

u1

u1

0

0

0

u1

u1

u1

f2 ðt; DtÞ

u2

0

u2

u2

0

u2

u2

0

u2

2.4. The 9-jump bivariate lattice of Lin [26] pffiffiffiffiffiffi Let ui  ki ri Dt and ki P 1 be stretch parameters, i = 1, 2. Lin proposes a 9-jump bivariate lattice which can handle two barriers [26]. The joint normal random variable {f1(t, Dt), f2(t, Dt)} is approximated by discrete random variable ff1 ðt; DtÞ; f2 ðt; DtÞg defined in Table 5. As Dt ? 0 the solutions of the transition probabilities are:

p1 ¼

1 þ 2q þ qk1 k2 4k21 k22 p5 ¼ p8 ¼

p2 ¼

;

1  2q þ k22

1 þ 2q  k21  k22 þ

2k21 k22 2 2 k1 k2

k21 k22 1  2q þ k22 2k21 k22

;

p9 ¼

;

p3 ¼

;

p6 ¼

1 þ 2q  qk1 k2

1  2q þ 2k21 k22

1 þ 2q þ qk1 k2 4k21 k22

4k21 k22 k21 ;

p7 ¼

;

p4 ¼

1  2q þ k21 2k21 k22

1 þ 2q  qk1 k2 4k21 k22

;

;

:

The stretch parameters can be properly chosen to align the lattice nodes with both barriers; hence it can handle two barriers. In order to make the probabilities valid, the values of q, k1 , and k2 have to satisfy the following inequalities:

k21 ; k22 P 1 þ 2q; ð1  k1 Þð1  k2 Þ P 2q;

qk1 k2 P ð1 þ 2qÞ;

ð3Þ

qk1 k2 P 1 þ 2q:

ð4Þ

From inequalities (3) and (4),

1 1 6q6 : 2 þ k1 k2 2  k1 k2 Because k1 k2 > 2 must hold for the above interval to be nonempty, it implies q > 1/4. Therefore, this lattice is not optimal. In summary, none of the lattices above enjoy both optimality and the ability to handle multiple barriers. 3. The hexanomial lattice Having analyzed the strengths and weaknesses of many existing bivariate lattices, we proceed to describe our hexanomial lattice. The hexanomial lattice will be proved to be optimal for bivariate options. Furthermore, it can handle options with a barrier on each asset by aligning the lattice nodes with both barriers. Indeed, the hexanomial lattice is the first such optimal bivariate lattice without relying on orthogonalization. It will be generalized to higher dimensions at the end of the section. Divide each time period into two phases with duration xDt in the first phase and (1  x)Dt in the second. The parameter 0 < x < 1 determines the duration of the two phases. The branching structure of the hexanomial lattice in a time period is shown in Fig. 1. To ensure that the hexanomial lattice converges to the true distribution as Dt ? 0, the log returns in a risk-neutral measure will be made to satisfy the following properties: a. In the first phase, the first moments, the second moments, and the correlation of fi ðt; DtÞ are:

E½fi ðt; xDtÞ ¼ li Dt;

ð5Þ

E½f2i ðt; xDtÞ  xr2i Dt;

ð6Þ

E½f1 ðt; xDtÞf2 ðt; xDtÞ  E½f1 ðt; xDtÞE½f2 ðt; xDtÞ q r1 r2 Dt

ð7Þ

When the equalities are the results of ignoring the higher-order terms, we use  instead of = (as Eqs. (6) and (7) above). b. In the second phase, the first moments, the second moments, and the correlation of fi ðt; DtÞ are:

E½fi ðt þ xDt; ð1  xÞDtÞ ¼ 0;

ð8Þ

468

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

[t , t + ωΔt ]

N

[t + ωΔt , t + Δt ]

N2

N1

N4

N3 N23

N22

N26

N25

N29

N43

N28

N42

N46 N48

N13

N24 N27

N33

N18

N32

N36 N39

N14 N17

N31

N35 N38

N11

N15

N19

N44 N47

N12

N16

N41

N45

N49

N21

N34 N37

Fig. 1. The branching structure of the hexanomial lattice in a time period. Node N with a pair of asset’s prices has 4 branches in the first phase with a duration of xDt. Each of nodes N1, N2, N3, and N4 has 9 branches in the second phase with a duration of (1  x)Dt.

E½f2i ðt þ xDt; ð1  xÞDtÞ  ð1  xÞr2i Dt;   E½f1 ðt þ xDt; ð1  xÞDtÞf2 ðt þ xDt; ð1  xÞDtÞ  E½f1 ðt þ xDt; ð1  xÞDtÞE½f2 ðt þ xDt; ð1  xÞDtÞ ¼ 0: r1 r2 Dt

ð9Þ ð10Þ

The intuition behind the construction is as follows. The first phase matches the first moment and the correlation for the whole time period in Eqs. (5) and (7), respectively. Therefore, in the second phase, the first moment should also be 0 as in Eq. (8), and the two assets’ log returns should be uncorrelated as in Eq. (10). Eqs. (6) and (9) finally guarantee that the sum of the second moments in the first and second phases equals the second moment for the whole time period, i.e., r2i Dt. We now proceed to prove that ff1 ðt; DtÞ; f2 ðt; DtÞg can be constructed to satisfy the above properties (a) and (b). 3.1. The first phase pffiffiffiffiffiffiffiffiffiffi Let ui1  ri xDt , i = 1, 2. The 4-jump binomial branch of Boyle et al. [7] is adopted in the first phase with a duration of xDt; hence the ith asset’s log return fi ðt; xDtÞ may increase or decrease by ui1 (recall Table 2). The discrete random variable ff1 ðt; DtÞ; f2 ðt; DtÞg is specified in Table 6. The topology of the hexanomial lattice in this phase forms the top pyramid of the structure in Fig. 1. Because the 4-jump binomial branch of Boyle et al. is adopted here but with a shorter duration xDt instead of Dt, the transition probabilities are similar to Eqs. (2a)–(2d):

p1 ¼

1 q 1þ þ x 4

1 q p2 ¼ 1 þ x 4 1 q 1 þ p3 ¼ x 4

p4 ¼

1 q 1þ þ x 4

rffiffiffiffiffiffi Dt l 1

x r1

þ

rffiffiffiffiffiffi Dt l 1

l2 r2

l  2 x r1 r2

! ;

ð11aÞ

;

ð11bÞ

!

rffiffiffiffiffiffi ! Dt l1 l2  þ ;

ð11cÞ

rffiffiffiffiffiffi ! Dt l l :  1 2

ð11dÞ

x

x

r1

r1

r2

r2

469

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479 Table 6 The branching structure in the first phase. Probability

p1

p2

p3

p4

f1 ðt; xDtÞ

u11

u11

u11

u11

f2 ðt; xDtÞ

u21

u21

u21

u21

It is easy to verify that property (a) is met because Eqs. (5)–(7) are satisfied. 3.2. The second phase A 9-jump trinomial branch is adopted in the second phase with a duration of (1  x)Dt. The assets’ log returns will be pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi made uncorrelated in this phase, matching Eq. (10). Let ui2  ki ri ð1  xÞDt and ki be stretch parameters, i = 1, 2. The ith asset’s log return fi ðt þ xDt; ð1  xÞDtÞ in a risk-neutral measure may increase or decrease by ui2, or stay constant. The discrete random variable ff1 ðt þ xDt; ð1  xÞDtÞ; f2 ðt þ xDt; ð1  xÞDtÞg is specified in Table 7. f1 and f2 are uncorrelated. The topology, the moves of assets’ log returns, and the transition probabilities in this phase are depicted in Fig. 2. The topology in this phase extends from the top pyramid in Fig. 1. We now solve the marginal transition probabilities in this phase as Dt ? 0. The first moments of the assets’ log returns should match Eq. (8):

  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E½fi ðt þ xDt; ð1  xÞDtÞ ¼ pi1 ki ri ð1  xÞDt þ pi2 ð0Þ þ pi3 ki ri ð1  xÞDt ¼ 0;

ð12Þ

where i = 1, 2. The second moments of the assets’ log returns should match Eq. (9):

  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 E½f2i ðt þ xDt; ð1  xÞDtÞ ¼ pi1 ki ri ð1  xÞDt þ pi2 ð0Þ2 þ pi3 ki ri ð1  xÞDt  r2i ð1  xÞDt;

ð13Þ

where i = 1, 2. Finally, the marginal transition probabilities should sum to 1:

pi1 þ pi2 þ pi3 ¼ 1;

i ¼ 1; 2

ð14Þ

Eqs. (12)–(14) result in the following solutions:

pi1 ¼

1 2k2i

;

pi2 ¼ 1 

1 k2i

;

pi3 ¼

1 2k2i

;

i ¼ 1; 2

ð15Þ

Combining the results for the two phases, the first moments of fi ðt; DtÞ are liDt (see Appendix A), the second moments of fi ðt; DtÞ converge to r2i Dt (see Appendix B), and the correlation between f1 ðt; DtÞ and f2 ðt; DtÞ converges to q (see Appendix C), where i = 1, 2. Thus the hexanomial lattice converges to the true distribution. The transition probabilities are also valid by the next theorem.

Table 7 The branching structure in the second phase. Probability

pi1

pi2

pi3

fi ðt þ xDt; ð1  xÞDtÞ

ui2

0

ui2

ζ 1 ( t + ωΔt , (1 − ω ) Δt ) , ζ 2 ( t + ωΔt , (1 − ω ) Δt ) N1

p13 , p21 N13

p13 , p22

p12 , p21

N12

N16

N15

N14

N11

p11 , p21

p11 , p22

+ λ2σ 2 (1 − ω ) Δt

p12 , p22 p13 , p23

N19

N18 −λ1σ 1 (1 − ω ) Δt

p12 , p23

N17

−λ2σ 2 (1 − ω ) Δt

p11 , p23

+ λ1σ 1 (1 − ω ) Δt

Fig. 2. The topology of the hexanomial lattice in the second phase. Node N1 is at time t + xDt, and nodes N11, . . ., N19 are at time t + Dt.

470

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

Theorem 1. When the following properties are satisfied, valid transition probabilities exist for the hexanomial lattice:

ðaÞ  1 þ

    pffiffiffiffiffiffil pffiffiffiffiffiffil l l Dt  1 þ 2  < q < 1  Dt  1  2 ; r r r r 1

2

1

and

2

ðbÞ k1 ; k2 P 1: Proof. See Appendix D.

h

Note that theprange hexanomial lattice is essentially identical to the 4-jump lattice of Boyle et al., which we ffiffiffiffiffiffi of q of the p ffiffiffiffiffiffi recall is 1 þ Oð Dt Þ 6 q 6 1  Oð DtÞ. We therefore have the optimality of the hexanomial lattice. Corollary 1. The hexanomial lattice is optimal. We make a few remarks here. First, when q = 1, 1, S1(t) and S2(t) are linearly related [32], building a valid lattice is trivial. Second, recall from Section 2.3 that the stretch parameters k1 and k2 can be adjusted to align the lattice nodes with a barrier on each asset. Theorem 1 guarantees that the resulting transition probabilities are valid. This feature will turn out to be very useful later when we price dual-strike options with a barrier on each asset. Third, the hexanomial lattice contains O(n5) nodes. The node count can be reduced to O(n3) if there are no barriers to match. See Appendix E for a proof. The O(n3)-sized hexanomial lattice can price bivariate options with a barrier on each asset but without aligning the lattice nodes with either barrier. Both hexonomial lattices are optimal of course. 3.3. The higher-dimensional hexanomial lattice Exact expressions for the transition probabilities can be obtained for higher-dimensional hexanomial lattices. In analogy to the bivariate case, each asset’s log return follows the binomial branch in the first phase with a duration of xDt and the trinomial branch in the second phase with a duration of (1  x)Dt. As a result, it is a 6k-jump lattice when there are k assets, with 6 branches in each asset. We proceed to obtain a set of equations for the probabilities. Let M  2k be the number of jumps and p‘ denote the transition probability of the ‘th jump in the first phase, where ‘ = 1, 2, . . ., M. Let pij denote the transition probability of the ith asset’s jth jump in the second phase, where j = 1, 2, 3. For the first phase, the p‘’s are similar to those in the bivariate case:

1

0

p‘ ¼

C B rffiffiffiffiffiffi k k X q 1B Dt X lC C B uij ð‘Þ ij þ di ð‘Þ i C; B1 þ C MB x x r i i¼1 A @ i; j ¼ 1

‘ ¼ 1; 2; . . . :; M;

i
uij ð‘Þ ¼



1;

if both the ith and jth assets move in the same direction in the ‘th jump ;

1; if both the ith and jth assets move in opposite directions in the ‘th jump:

For the second phase, the counterpart to Eq. (12) is a set of k equations for each asset:

  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi1 ki ri ð1  xÞDt þ pi2 ð0Þ þ pi3 ki ri ð1  xÞDt ¼ 0;

16i6k

ð16Þ

The counterpart to Eq. (13) is

  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pi1 ki ri ð1  xÞDt þ pi2 ð0Þ2 þ pi3 ki ri ð1  xÞDt  r2i ð1  xÞDt;

16i6k

ð17Þ

Similar to Eq. (14), we have

pi1 þ pi2 þ pi3 ¼ 1;

1 6 i 6 k:

ð18Þ

In total, Eqs. (16)–(18) constitute a set of 3k equations for the 3k unknown pij ’s, where 1 6 i 6 k and j = 1, 2, 3. The solutions are

pi1 ¼

1 2k2i

;

pi2 ¼ 1 

1 k2i

;

pi3 ¼

1 2k2i

;

16i6k

Appendix E shows that the higher-dimensional hexanomial lattice contains O(n2k+1) nodes, and the node count is reduced to O(nk+1) if there are no barriers to match. Although the hexanomial lattice has been generalized to handle more than 2 underlying assets, the proof for optimality remains elusive. Thus the probabilities and whether 0 < x < 1 must be checked in practice.

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

471

4. Numerical results Because Boyle et al.’s 4-jump bivariate and Rubinstein’s 4-jump bivariate lattices are only for bivariate options, and the transition probabilities of Boyle’s 5-jump bivariate lattice cannot be guaranteed to be valid when the number of assets is more than 2, all the numerical experiments in the paper are for bivariate options. Four bivariate options will be priced: the maximum option, the spread option, the dual-strike option, and the dual-strike option with a knock-out barrier on each asset. Only the last option has barriers. The pricing results and convergence behaviors of the O(n3)-sized and the O(n5)-sized hexanomial lattices are then compared with the 5-jump bivariate lattice of Boyle [6], the 4-jump bivariate lattice of Boyle et al. [7], and the 4-jump bivariate lattice of Rubinstein [34], as n increases. The 9-jump lattice of Lin [26] has similar performance to the hexanomial lattice as long as the correlation is within its limited range; hence no numerical experiments are carried out for it. For the maximum option, the spread option, and the dual-strike option, the O(n3)-sized hexanomial lattice is employed. For the dual-strike option with a knock-out barrier on each asset, both the O(n3)-sized and the O(n5)-sized hexanomial lattices are used. Next, we show the pricing results and convergence behaviors of the above-mentioned lattices with respect to CPU times when pricing the dual-strike option with a knock-out barrier on each asset. We also present the CPU time of the O(n3)-sized hexanomial lattice to achieve the same accuracy as the alternating direction implicit (ADI) method when pricing the maximum option. Finally, experiments are carried out for the O(n3)-sized hexanomial lattice when q is very close to 1 or 1. 4.1. The maximum option Recall that the payoff function of a maximum option at the maturity day T is max[max(S1(T), S2(T))  K, 0]. The parameters in the experiments are S1(0) = 40, S2(0) = 40, K = 40, r1 = 30% r2 = 30% q = 0.5, r = 10% and T = 12 (months). Johnson’s formula gives the benchmark price of 9.956044 [20]. Fig. 3 depicts the results for n between 5 and 250. The convergence of the hexanomial lattice is faster than Boyle et al.’s 4-jump bivariate lattice and Boyle’s 5-jump bivariate lattice. In particular, Boyle et al.’s 4-jump bivariate lattice exhibits significant oscillations. The convergence of Rubinstein’s 4-jump bivariate lattice is the fastest, but it is not as smooth as the hexanomial lattice or Boyle’s 5-jump bivariate lattice. Let the pricing error of the O(n3)-sized hexanomial lattice for the maximum option be the price by the lattice minus the price given by Johnson’s formula. The loglog plot of the pricing error for n between 600 and 1000 appears in Fig. 4. The pricing error is clearly of order Dt. 4.2. The spread option Recall that the payoff function of a spread option at maturity date T is max[S2(T)  S1(T)  K, 0]. The parameters in the experiments are S1(0) = 40, S2(0) = 50, K = 5, r1 = 20% r2 = 18% q = 0.6, r = 5% and T = 12 (months). Fig. 5 depicts the results for n between 5 and 250. The convergence of Rubinstein’s 4-jump bivariate lattice is faster and smoother than Boyle et al.’s 4-jump bivariate lattice and Boyle’s 5-jump bivariate lattice. Both Boyle et al.’s 4-jump bivariate lattice and Boyle’s 5-jump bivariate lattice exhibit oscillations. Overall, the convergence of the hexanomial lattice is faster and smoother than all other lattices.

10.1 10.05 10 9.95 9.9

Option Price

9.85 9.8 9.75 9.7

O(n^3)-sized hexanomial lattice Boyle et al.'s 4-jump bivariate lattice Rubinstein's 4-jump bivariate lattice Boyle's 5-jump bivariate lattice Johnson's formula

9.65 9.6 9.55 9.5 9.45 9.4 5

20

35

50

65

80

95

110 125 140 155 170 185 200 215 230 245

Number of Time Steps Fig. 3. Price vs. n for the maximum option.

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

Log of the Pricing Error of the Hexanomial Lattice

472

-3.567

-3.617

-3.667

-3.717

-3.767 -3

-2.95

-2.9

-2.85

-2.8

log(Δt) Fig. 4. The log of the pricing error of the hexanomial lattice.

4.3. The dual-strike option Recall that the payoff function of a dual-strike option at maturity date T is max [S1(T)  K1, S2(T)  K2, 0]. The parameters in the experiments are S1(0) = 40, S2(0) = 50, K1 = 35, K2 = 45, r1 = 20% r2 = 18% q = 0.6, r = 5% and T = 12 (months). Fig. 6 depicts the results for n between 5 and 250. Rubinstein’s 4-jump bivariate lattice converges faster than Boyle et al.’s 4-jump bivariate lattice and Boyle’s 5-jump bivariate lattice. Boyle et al.’s 4-jump bivariate lattice exhibits mild oscillations. Overall, the convergence of the hexanomial lattice is faster and smoother than all other lattices. 4.4. The dual-strike option with a knock-out barrier on each asset This Section prices the same option as Section 4.3 but with a knock-out barrier on each asset. If any asset’s price touches its barrier, the option knocks out and becomes worthless. Both the O(n3)-sized and the O(n5)-sized hexanomial lattices are used in the experiments. The parameters in the experiments are S1(0) = 40, S2(0) = 50, K1 = 35, K2 = 45, H1 = 50, H2 = 60, r1 = 20% r2 = 18% q = 0.6, r = 5% and T = 12 (months). To align a layer of lattice nodes with a barrier, O(n5)-sized hexanopffiffiffiffiffiffiffiffiffi ffi pthe ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mial lattice makes sure that it takes an integer number of steps for the log return ri xDt þ ki ri ð1  xÞDt to reach the barrier Hi from Si(0) for both i = 1, 2. Thus we choose ki P 1, to the extent possible, such that

pffiffiffiffiffiffiffiffiffiffi lnðHi =Si ð0ÞÞ  jri xDt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; j¼1;2;3;... jri ð1  xÞDt

ki ¼ min

following the idea of Ritchken [31]. 6.36

Option Price

6.35

6.34

6.33

O(n^3)-sized hexanomial lattice Boyle et al.'s 4-jump bivariate lattice

6.32

Rubinstein's 4-jump bivariate lattice Boyle's 5-jump bivariate lattice 6.31 20

35

50

65

80

95

110

125

140

155

170

Number of Time Periods Fig. 5. Price vs. n for the spread option.

185

200

215

230

245

473

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

10.4 10.35

Option Price

10.3 10.25 10.2 10.15

O(n^3)-sized hexanomial lattice Boyle et al.'s 4-jump bivariate lattice Rubinstein's 4-jump bivariate lattice Boyle's 5-jump bivariate lattice

10.1 10.05 10 5

20

35

50

65

80

95

110 125 140 155 170 185 200 215 230 245

Number of Time Periods Fig. 6. Price vs. n for the dual-strike option.

Fig. 7 depicts the results for n between 5 and 150. The convergence of Boyle’s 5-jump bivariate lattice is slightly smoother than Boyle et al.’s 4-jump bivariate lattice and the O(n3)-sized hexanomial lattice. Recall that the O(n3)-sized hexanomial lattice cannot align the lattice nodes with either barrier. Overall, the convergence of the O(n5)-sized hexanomial lattice is smoother and with smaller oscillations than all other lattices. 4.5. The CPU times of various methods We measure the O(n3)-sized and O(n5)-sized hexanomial, Boyle et al.’s 4-jump bivariate, Rubinstein’s 4-jump bivariate, and Boyle’s 5-jump bivariate lattices’ CPU times using the high-resolution performance clock counter of Windows 7. The C codes run on the Intel Core i5-2520M CPU. All times are average CPU times of 100 runs of the programs to minimize the estimation error. Fig. 8 shows the convergence behaviors of the O(n5)-sized hexanomial, Boyle et al.’s 4-jump bivariate, Rubinstein’s 4-jump bivariate, and Boyle’s 5-jump bivariate lattices with respect to CPU times when pricing the dual-strike option with a knock-out barrier on each asset in Section 4.4. It shows the convergence behavior of the O(n5)-sized hexanomial lattice to be smoother than all other lattices because it can align lattice nodes with both barriers. Next, we measure the CPU time of the O(n3)-sized hexanomial lattice to achieve the same accuracy as the ADI method in Table 2 of [19] when pricing the maximum option with the same parameters. The root mean square error (RMSE) is used to measure the accuracy. There are three accuracy levels. Table 8 shows CPU times and the RMSEs of the O(n3)-sized hexanomial lattice to achieve three accuracy levels of the ADI method. Note that the O(n3)-sized hexanomial lattice easily achieves level 1 and level 2 accuracies. Even for the higher level-3 accuracy, the hexanomial lattice attains it in about 2 min.

5.6

O(n^5)-sized hexanomial lattice 5.1

O(n^3)-sized hexanomial lattice Boyle et al.'s 4-jump bivariate lattice

Option Price

4.6

Rubinstein's 4-jump bivariate lattice 4.1

Boyle's 5-jump bivariate lattice

3.6 3.1 2.6 2.1 5

20

35

50

65

80

95

110

125

Number of Time Periods Fig. 7. Price vs. n for the dual-strike option with a barrier on each asset.

140

474

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

2.9

Boyle et al.'s 4-jump bivariate lattice 2.85

Boyle's 5-jump bivariate lattice Rubinstein's 4-jump bivariate lattice

Option Price

2.8

O(n^5)-sized hexanomial lattice 2.75 2.7 2.65 2.6 2.55 2.5 10

15

20

25

30

35

40

45

50

CPU Time (sec) Fig. 8. Price vs. CPU time for the dual-strike option with a barrier on each asset.

Table 8 The CPU times of the O(n3)-sized hexanomial lattice to achieve three accuracy levels. Accuracy

Level 1 Level 2 Level 3

ADI

O(n3)-sized hexanomial lattice

RMSE

RMSE

CPU time (sec)

n

0.006319 0.001581 0.000395

0.006102 0.001530 0.000377

0.073344 1.242614 129.982408

55 140 650

Table 9 The prices of the spread option computed by the hexanomial lattice.

q

0.95

0.99

1

0.95

0.99

1

Option price

5.269983

5.243811

5.243805

9.609611

9.674285

9.691652

Parameters: n = 250 and q = 1, 0.99, 0.95, 0.95, 0.99, 1.

4.6. The cases with near ±1 correlations When the correlation in the hexanomial lattice is very close to 1 or 1, using 1 or 1 as the correlation yields very good approximations with a small Dt. Take the spread option as an example, Table 9 illustrates the prices of the spread option using the hexanomial lattice with the same parameters as Section 4.2 except q and Dt = 1/250. Observe that the prices with q = 0.95 and q = 0.99 are close to that with q = 1, and the prices with q = 0.95 and q = 0.99 are close to that with q = 1. 5. Conclusions This ppaper calls a lattice ffiffiffiffiffiffi pffiffiffiffiffiffi (correlation) optimal if the transition probabilities are guaranteed to be valid for 1 þ Oð DtÞ 6 qij 6 1  Oð DtÞ for all pairs of assets i and j, where qij is the correlation between assets i and j. It is clearly a desirable property for lattices. The hexanomial lattice is presented, which has the following features:  The lattice is optimal.  The lattice nodes can be aligned with both barriers on the assets. The hexanomial lattice is also extended to higher dimensions. Numerically, the hexanomial lattice exhibits faster and smoother convergence than alternative bivariate lattices in pricing bivariate options. This is particularly prominent in the case of options with a barrier on each asset. Adopting 1 or 1 as the correlation yields very good approximations with a small Dt when the correlation in the hexanomial lattice is very close

475

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

to 1 or 1. Optimality in the paper refers to the fact that essentially all correlations between 1 and 1 can be solved. Surprisingly, this important notion has been overlooked in the literature. Acknowledgment This work was supported in part by the National Science Council of Taiwan under Grant 100-2221-E-002-111-MY3. Appendix A. The first moments of fi ðt; tÞ over one time period The first moments of fi ðt; DtÞ, i = 1, 2, equal

E½fi ðt; DtÞ ¼ E½fi ðt; xDtÞ þ fi ðt þ xDt; ð1  xÞDtÞ ¼ E½fi ðt; xDtÞ þ E½fi ðt þ xDt; ð1  xÞDtÞ:

ðA:1Þ

By Eqs. (5) and (12),

Eq:ðA:1Þ ¼ li Dt þ 0 ¼ li Dt: Thus the first moments of fi ðt; DtÞ are liDt. Appendix B. The second moments of fi ðt; tÞ over one time period The second moments of fi ðt; DtÞ, i = 1, 2, equal

E½f2i ðt; DtÞ ¼ E½f2i ðt; xDtÞ þ E½f2i ðt þ xDt; ð1  xÞDtÞ þ 2E½fi ðt; xDtÞfi ðt þ xDt; ð1  xÞDtÞ:

ðB:1Þ

There are three expectations to evaluate. By Eqs. (6) and (13), the first and second expectations approximate xr2i Dt and ð1  xÞr2i Dt, respectively. The third expectation equals

 pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ri xDt ki ri ð1  xÞDt ðpi1  pi3 Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p2 ri xDt ki ri ð1  xÞDt ðpi1  pi3 Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p3 ri xDt ki ri ð1  xÞDt ðpi1  pi3 Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p4 ri xDt ki ri ð1  xÞDt ðpi1  pi3 Þ:

E½fi ðt; xDtÞfi ðt þ xDt; ð1  xÞDtÞ ¼ p1

It is easy to check that Eq. (B.2) = 0 with Eqs. (11a)–(11d). Combining the above results, we have Eq: ðB:1Þ  r2i Dt. Thus the second moments of fi ðt; DtÞ converge to

ðB:2Þ

r2i Dt.

Appendix C. The correlation of f1 ðt; tÞ and f2 ðt; tÞ over one time period To make sure that the correlation of f1 ðt; DtÞ and f2 ðt; DtÞ converges to q, we need to prove

E½f1 ðt; DtÞf2 ðt; DtÞ  E½f1 ðt; DtÞE½f2 ðt; DtÞ q r 1 r 2 Dt

ðC:1Þ

The left-hand side of Eq. (C.1) equals 1

r1 r2 Dt

  E ðf1 ðt; xDtÞ þ f1 ðt þ xDt;ð1  xÞDtÞÞ  ðf2 ðt; xDtÞ þ f2 ðt þ xDt;ð1  xÞDtÞÞ  E½f1 ðt; DtÞE½f2 ðt; DtÞ

qr1 r2 Dt þ E½f1 ðt; xDtÞE½f2 ðt; xDtÞ þ E½f1 ðt; xDtÞf2 ðt þ xDt;ð1  xÞDtÞ  r1 r2 Dt þE½f1 ðt þ xDt;ð1  xÞDtÞf2 ðt; xDtÞ þ E½f1 ðt þ xDt; ð1  xÞDtÞf2 ðt þ xDt; ð1  xÞDtÞ  E½f1 ðt; DtÞE½f2 ðt; DtÞ 1

by Eq. (7). There are five expectations to evaluate. By Eq. (5), E½fi ðt; xDtÞ ¼ li Dt, i = 1, 2. So the first expectation equals l1l2Dt2. The second expectation equals

!

476

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

E½f1 ðt; xDtÞf2 ðt þ xDt; ð1  xÞDtÞ ¼ p1 p21



pffiffiffiffiffiffiffiffiffiffi 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



pffiffiffiffiffiffiffiffiffiffi

r1 xDt k2 r2 ð1  xÞDt þ p1 p22 r1 xDt ð0Þ

þ p1 p23 þ p2 p21 þ p1 p23

  

pffiffiffiffiffiffiffiffiffiffi 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

r1 xDt k2 r2 ð1  xÞDt pffiffiffiffiffiffiffiffiffiffi 



pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

pffiffiffiffiffiffiffiffiffiffi

r1 xDt k2 r2 ð1  xÞDt þ p1 p22 r1 xDt ð0Þ pffiffiffiffiffiffiffiffiffiffi 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

r1 xDt k2 r2 ð1  xÞDt

  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p2 p21 r1 xDt k2 r2 ð1  xÞDt þ p2 p22 r1 xDt ð0Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p2 p23 r1 xDt k2 r2 ð1  xÞDt  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p4 p21 r1 xDt k2 r2 ð1  xÞDt þ p4 p22 ðr1 xDtÞð0Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p4 p23 r1 xDt k2 r2 ð1  xÞDt :

ðC:3Þ

It is easy to verify Eq. (C.3) = 0 with Eqs. (11a)–(11d) and (15). The third expectation equals

E½f1 ðt þ xDt; ð1  xÞDtÞf2 ðt; xDtÞ ¼ p1 p11



pffiffiffiffiffiffiffiffiffiffi 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



pffiffiffiffiffiffiffiffiffiffi

r2 xDt k1 r1 ð1  xÞDt þ p1 p12 r2 xDt ð0Þ

þ p1 p13



pffiffiffiffiffiffiffiffiffiffi 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

r2 xDt k1 r1 ð1  xÞDt

  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p2 p11 r2 xDt k1 r1 ð1  xÞDt þ p2 p12 r2 xDt ð0Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p2 p13 r2 xDt k1 r1 ð1  xÞDt  pffiffiffiffiffiffiffiffiffiffi   pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p3 p11 r2 xDt k1 r1 ð1  xÞDt þ p3 p12 r2 xDt ð0Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p3 p13 r2 xDt k1 r1 ð1  xÞDt   pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p4 p11 r2 xDt k1 r1 ð1  xÞDt þ p4 p12 r2 xDt ð0Þ  pffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p4 p13 r2 xDt k1 r1 ð1  xÞDt :

ðC:4Þ

It is easy to verify Eq. (C.4) = 0 with Eqs. (11a)–(11d) and (15). Because f1 ðt þ xDt; ð1  xÞDtÞ and f2 ðt þ xDt; ð1  xÞDtÞ are uncorrelated in the second phase and because of Eq. (12), the fourth expectation equals 0. By Appendix A, the fifth expectation equals l1l2Dt2. Combining the above results, we have Eq. (C.2)  q. Thus the correlation of f1 ðt; DtÞ and f2 ðt; DtÞ converges to q. Appendix D. Proof of Theorem 1

Theorem 1. When the following properties are satisfied, valid transition probabilities exist for the hexanomial lattice:

ðaÞ  1 þ

    pffiffiffiffiffiffil pffiffiffiffiffiffil l l Dt  1 þ 2  < q < 1  Dt  1  2 ; r r r r 1

2

1

and

2

ðbÞ k1 ; k2 P 1: Proof. First, we prove that a 0 < x < 1 exists to make the transition probabilities in the first phase valid as long as

1 þ

    pffiffiffiffiffiffil pffiffiffiffiffiffil l l Dt  1 þ 2  < q < 1  Dt  1  2 : r r r r 1

2

1

2

. The following inequalities are Eqs. (11a)–(11d) after being multiplied by 4x. All must be nonnegative to let the transition probabilities in the first phase valid:

pffiffiffiffiffi

ðD:1aÞ

pffiffiffiffiffi

ðD:1bÞ

x þ m1 x  q P 0; x  m1 x  q P 0;

477

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

pffiffiffiffiffi

ðD:1cÞ

pffiffiffiffiffi

ðD:1dÞ

x þ m2 x þ q P 0; x  m2 x þ q P 0;

pffiffiffiffiffiffi pffiffiffiffiffiffi where m1  Dt½ðl1 =r1 Þ  ðl2 =r2 Þ and m2  Dt ½ðl1 =r1 Þ þ ðl2 =r2 Þ. Because x is within the square root and in the denominators of Eqs. (11a)–(11d), it has to be larger than 0. pffiffiffiffiffi pffiffiffiffiffi Now consider inequalities (D.1a) and (D.1b). The discriminants of x þ m1 x  q ¼ 0 and x  m1 x  q ¼ 0 are both 2 2 2 m1 þ 4q. Suppose m1 þ 4q < 0. Inequalities (D.1a) and (D.1b) hold as long as x > 0. Suppose m1 þ 4q P 0, on the other hand. To make inequalities (D.1a) and (D.1b) true, either of the following inequalities must be satisfied:

pffiffiffiffiffi

xP

pffiffiffiffiffi

x6

jm1 j þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m21 þ 4q 2

jm1 j 

;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m21 þ 4q 2

ðD:2Þ

:

But inequality (D.2) contradicts the x > 0 requirement. In summary, to make inequalities (D.1a) and (D.1b) true, x has to satisfy the following inequalities:

x > 0; if m21 þ 4q < 0 0

xP@

jm1 j þ

ðD:3Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 m21 þ 4q A ; if m2 þ 4q P 0: 1 2

ðD:4Þ

pffiffiffiffiffi pffiffiffiffiffi Consider inequalities (D.1c) and (D.1d) next. The discriminants of x þ m2 x þ q ¼ 0 and x  m2 x þ q ¼ 0 are both 2 2  4q. Suppose m2  4q < 0. Inequalities (D.1c) and (D.1d) hold as long as x > 0. Suppose m2  4q P 0, on the other hand. To make inequalities (D.1c) and (D.1d) true, either of the following inequalities must be satisfied: m22

pffiffiffiffiffi

xP

pffiffiffiffiffi

x6

jm2 j þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m22  4q 2

jm2 j 

;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m22  4q 2

ðD:5Þ

:

But inequality (D.5) contradicts the x > 0 requirement. In summary, to make inequalities (D.1c) and (D.1d) true, x has to satisfy the following inequalities:

x > 0; if m22  4q < 0; 0

xP@

jm2 j þ

ðD:6Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 m22  4q A ; if m2  4q P 0: 2 2

ðD:7Þ

Then we determine the interval for x to satisfy inequalities (D.1a)–(D.1d) in 4 different cases of the range for q. The first case is m21 þ 4q < 0 and m22  4q < 0, which obviously is impossible. The second case is m21 þ 4q P 0 and m22  4q < 0. As m22  4q < 0 implies m21 þ 4q P 0, x should satisfy

0

xP@

jm1 j þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 m21 þ 4q A ; 2

ðD:8Þ

by inequality (D.4). The third case is m21 þ 4q < 0 and m22  4q P 0. As m21 þ 4q < 0 implies m22  4q P 0, x should satisfy

0

xP@

jm2 j þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 m22  4q A ; 2

ðD:9Þ

by inequality (D.7). The forth case is m21 þ 4q P 0 and m22  4q P 0. By inequalities (D.4) and (D.7), x should satisfy

80 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 9 > < jm1 j þ m21 þ 4q = jm2 j þ m22  4q > A ; @ A : x P max @ > > 2 2 : ;

ðD:10Þ

478

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

In summary, from inequalities (D.8)–(D.10), we have

8 pffiffiffiffiffiffiffiffiffiffiffi2 jm1 jþ m21 þ4q > > ; > > 2 > > > > 2 > < jm jþpffiffiffiffiffiffiffiffiffiffiffi m22 4q 2 ; xP 2 > > > ( >  > pffiffiffiffiffiffiffiffiffiffiffi2  pffiffiffiffiffiffiffiffiffiffiffi2 ) > > jm1 jþ m21 þ4q jm2 jþ m22 4q > max > ; ; : 2 2

if m21 þ 4q P 0 and m22  4q < 0; if m21 þ 4q P 0 and m22  4q < 0; if m21 þ 4q P 0 and m22  4q < 0;

ðD:11Þ ðD:12Þ ðD:13Þ

It is easy to verify that the union of the ranges of q on the right-hand side is the whole real line [1, 1]. In conclusion, for any q, an interval for x exists to guarantee valid transition probabilities in the first phase. By the lattice construction, x should be less than 1. As a result, we consider each of the three inequalities (D.11)–(D.13) with the additional x < 1 constraint. From inequality (D.11), the following inequality has to be satisfied:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 0 jm1 j þ m21 þ 4q @ A < 1; 2

ðD:14Þ

where m21 þ 4q P 0 and m22  4q < 0. Inequality m22  4q < 0 does not have to be considered in this case because x can take any positive value by inequality (D.6). Inequality (D.14) is equivalent to

q < 1  jm1 j and m21 þ 4q P 0:

ðD:15Þ

From inequality (D.12), the following inequality has to be satisfied:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 0 jm2 j þ m22  4q @ A < 1; 2

ðD:16Þ

where m21 þ 4q < 0 and m22  4q P 0. Inequality m21 þ 4q < 0 does not have to be considered in this case because x can take any positive value by inequality (D.3). Inequality (D.16) is equivalent to

q > 1 þ jm2 j and m22  4q P 0:

ðD:17Þ

From inequality (D.13), the following inequalities have to be satisfied:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 0 jm1 j þ m21 þ 4q @ A <1 2

ðD:18Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi12 0 jm2 j þ m22  4q @ A <1 2

ðD:19Þ

where m21 þ 4q P 0 and m22  4q P 0. Inequalities (D.18) and (D.19) are the same as inequalities (D.14) and (D.16), respectively; thus the same inequalities (D.15) and (D.17) are obtained. In summary, from inequalities (D.15) and (D.17) and recalling the definitions of m1 and m2, we have

    pffiffiffiffiffiffil pffiffiffiffiffiffil l2  l2  1 1   1 þ Dt  þ  < q < 1  Dt   : r r r r 1

2

1

ðD:20Þ

2

We move onto the second phase. It is obvious that all the probabilities in Eq. (15) are all valid when ki P 1, i = 1, 2. In conclusion, all the transition probabilities are valid if q satisfies inequalities (D.20) and k1 ; k2 P 1 h We remark that our implementation chooses x according to (D.11)–(D.13) after replacing ‘‘x P’’ there with ‘‘ x=.’’ Appendix E. The size of the hexanomial lattice First we consider the size of the hexanomial lattice for bivariate options. The number of nodes for each asset in the ith time period is O(i2). Thus the total number of nodes is n X 4 Oði Þ ¼ Oðn5 Þ: i¼1

W.-H. Kao et al. / Applied Mathematics and Computation 233 (2014) 463–479

479

When pricing bivariate options without barriers, as the stretch parameters do not have barrier locations to match, the node count may p beffiffiffiffiffiffiffiffiffi much reduced. The log return in the first phase can be made to equal that in the second phase for each asset as ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi follows: ri xDt ¼ ki ri ð1  xÞDt ; as a result, ki ¼ x= 1  x, i = 1, 2. The total number of nodes is now drastically reduced to n X 2 Oði Þ ¼ Oðn3 Þ: i¼1

pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi It remains to choose a x so that x= 1  x P 1. By matching q to the conditions of inequalities (D.11)–(D.13), a nonempty interval [x, 1) is determined for x, where 0 < x < 1. Any x in that interval guarantees valid transition probabilities by Theorem 1. We next focus on the potentially narrower interval max(0.5, x), 1). Note that this interval is nonempty because pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi max(0.5, x) < 1. Clearly, any x in this interval makes x= 1  x P 1. Next, the size of the higher-dimensional hexanomial lattice is determined. When pricing options, the total number of nodes of the hexanomial lattice for k-asset options is n X k Oði Þ ¼ Oðnkþ1 Þ: i¼1

When pricing options with barriers, the total number of nodes of the hexanomial lattice for k-asset options is n X 2k Oði Þ ¼ Oðn2kþ1 Þ: i¼1

References [1] W.F. Ames, Numerical Methods for Partial Differential Equations, third ed., Academic Press, New York, 1992. [2] C. Alexander, A. Venkatramanan, Analytic approximations for spread options, ICMA Centre Discussion Papers in Finance DP2007-11, University of Reading, United Kingdom, 2007, pp. 1–23. [3] F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Political Economy 81 (1973) 637–654. [4] M.J. Brennan, E.S. Schwartz, Finite different methods and jump processes arising in the pricing of contingent claims: a synthesis, J. Financial Quant. Anal. 13 (1978) 461–474. [5] P. Boyle, Options: a Monte Carlo approach, J. Financial Econ. 4 (1977) 323–338. [6] P. Boyle, A lattice framework for option pricing with two state variables, J. Financial Quant. Anal. 23 (1988) 1–12. [7] P. Boyle, J. Evnine, S. Gibbs, Numerical evaluation of multivariate contingent claims, Rev. Financial Stud. 2 (1989) 241–250. [8] P. Boyle, S. Lau, Bumping against the barrier with the binomial method, J. Derivatives 1 (1994) 6–14. [9] P. Boyle, M. Broadie, P. Glasserman, Monte Carlo methods for security pricing, J. Econ. Dyn. Control 21 (1997) 1267–1321. [10] P. Boyle, Y. Tian, An explicit finite difference approach to the pricing of barrier options, Appl. Math. Finance 5 (1998) 17–43. [11] D. Brigo, F. Mercurio, Interest Rate Models – Theory and Practice: With Smile, Inflation and Credit, second ed., Springer, Wien, NY, 2006. [12] R. Carmona, V. Durrleman, Pricing and hedging spread options, SIAM Rev. 45 (2003) 627–685. [13] L.-B. Chang, K.J. Palmer, Smooth convergence in the binomial model, Finance Stoch. 11 (2007) 91–105. [14] D. Duffie, Dynamic Asset Pricing Theory, third ed., Princeton University Press, Princeton, NJ, 2001. [15] S. Figlewski, B. Gao, The adaptive mesh model: a new approach to efficient option pricing, J. Financial Econ. 53 (1999) 313–351. [16] M. Gillia, E. Këllezib, G. Paulettoa, Solving finite difference schemes arising in trivariate option pricing, J. Econ. Dyn. Control 26 (2002) 1499–1515. [17] J.C. Hull, A. White, Valuing derivative securities using the explicit finite difference method, J. Financial Quant. Anal. 25 (1990) 87–100. [18] J.C. Hull, Options, Futures, and Other Derivatives, seventh ed., Prentice Hall, Englewood Cliffs, NJ, 2008. [19] D. Jeong, J. Kim, A comparison study of ADI and operator splitting methods on option pricing models, J. Comput. Appl. Math. 247 (2013) 162–171. [20] H. Johnson, Options on the maximum or the minimum of several assets, J. Financial Quant. Anal. 22 (1987) 277–283. [21] C. Joy, P. Boyle, K.S. Tan, Quasi-Monte Carlo methods in numerical finance, Manage. Sci. 42 (1996) 926–938. [22] B. Kamrad, P. Ritchken, Multinomial approximating models for options with k state variables, Manage. Sci. 37 (1991) 1640–1652. [23] A.G.Z. Kemna, A.C.F. Vorst, A pricing method for options based on average asset values, J. Banking Finance 14 (1990) 113–129. [24] E. Kirk, J. Aron, Correlation in the energy markets, in: V. Kaminski (Ed.), Managing Energy Price Risk, Risk Pub, London, 1995, pp. 71–78. [25] A. Kurpiel, T. Roncalli, Hopscotch methods for two-state financial models, J. Comput. Finance 3 (1999) 53–89. [26] V.-D. Lin, Pricing barrier options in two dimensions (Master’s Thesis), Department of Finance, National Taiwan University, Taipei, Taiwan, 2006. [27] Y.-D. Lyuu, Very fast algorithms for barrier option pricing and the ballot problem, J. Derivatives 5 (1998) 68–79. [28] Y.-D. Lyuu, Financial Engineering and Computation: Principles, Mathematics, Algorithms, Cambridge University Press, Cambridge, 2002. [29] W. Margrabe, The value of an option to exchange one asset for another, J. Finance 33 (1978) 177–186. [30] N.D. Pearson, An efficient approach for pricing spread options, J. Derivatives 3 (1995) 76–91. [31] P. Ritchken, On pricing barrier options, J. Derivatives 3 (1995) 19–28. [32] R.R. Roussas, A Course in Mathematical Statistics, second ed., Academic Press, New York, 1997. [33] M. Rubinstein, Somewhere over the rainbow, Risk 4 (1991) 63–66. [34] M. Rubinstein, Return to OZ, Risk 7 (1994) 67–71. [35] S. Shreve, Stochastic Calculus for Finance: Continuous-Time Models, Springer, Berlin, 2004. [36] R. Stulz, Options on the minimum or the maximum of two risky assets: analysis and applications, J. Financial Econ. 10 (1982) 161–185. [37] K.S. Tan, P. Boyle, Applications of scrambled low discrepancy sequences to exotic options, in: Proceedings of International AFIR Colloquium, in North Queensland, Australia, August 14–15, 1997. [38] D. Tavella, C. Randall, Pricing Financial Instruments: The Finite Difference Method, Wiley, Hoboken, NJ, 2000. [39] C.-J. Wang, T.-S. Dai, Y.-D. Lyuu, A multi-phase, flexible, and accurate lattice for pricing complex derivatives with multiple market variables, in: Proceedings of the IEEE Conference on Computational Intelligence for Financial Engineering & Economics, New York City, March 29–30, 2012. [40] X. Wang, K.-S. Tan, How do path generation methods affect the accuracy of quasi-Monte Carlo methods for problems in finance?, J Complexity 28 (2012) 250–277. [41] U. Wystup, FX Options and Structured Products, Wiley, Hoboken, NJ, 2006. [42] R. Zvan, P.A. Forsyth, K.R. Vetzal, Negative coefficients in two-factor option pricing models, J. Comput. Finance 7 (2003) 37–73.