國立臺灣大學電機資訊學院資訊工程學系碩士論文

(1)

國立臺灣大學電機資訊學院資訊工程學系碩士論文

Department of Computer Science and Information Engineering College of Electrical Engineering & Computer Science

National Taiwan University Master Thesis

評價巴黎選擇權之財務演算法：

組合學、模擬法與平行處理

Pricing Parisian Options:

Combinatorics, Simulation, and Parallel Processing

吳承瑋 Cheng-Wei Wu

指導教授：呂育道博士 Advisor: Yuh-Dauh Lyuu, Ph.D.

中華民國 97 年 6 月

June, 2008

(2)

摘要摘要摘要摘要

財務工程與金融創新在過去數十年蓬勃發展，設計出許多新金融商品，提供了風險管理所需的避險工具並促進市場效率與完整性。此財務領域的定價問題會嘗試建構數學模型推導公式解，但是由於大部分新奇衍生性金融商品契約複雜無公式能套用，必須借重電腦運算處理數值方法及模擬價格，因此計算機科學在此有其可施展之處。本篇論文探討巴黎選擇權評價方法即包含財務理論、機率統計、離散數學、計算複雜度、演算法設計與分析以及平行處理等議題。

巴黎選擇權為路徑相關選擇權，迄今尚無封閉公式解存在，為此我們提出兩種快速財務演算法。首先以 Costabile 於 2002 年建立在離散結構二項樹狀模型下的組合方法為基礎去評價巴黎選擇權，再將此方法稍加修改可把原本的時間複雜度由^{O n}

( )

³ ^{成功推進至}^{O n}

( )

² ；若是組合數已給定的情況下，空間複雜度亦由

( )

²

O n 減少到^{O n}

( )

。另外，在蒙地卡羅模擬法方面，我們引進逆高斯機率分配並結合其抽樣，可減少時間切割的期數而縮短計算時間。因為蒙地卡羅法的路徑模擬有各自獨立的特性，很容易使用平行運算處理。今日多核心處理器大行其道這不失為提高效率的方法，對此我們也做了相關的說明並加以應用。

我們以 C 語言實作論文中的財務演算法，執行於自行建構的高效能叢集運算平台。同步處理模擬之計算量，俾能有效使用系統效能。

關鍵字關鍵字關鍵字

關鍵字：：：：巴黎選擇權, 障礙選擇權, 選擇權評價, 演算法, 二項樹模型, 組合方法, 蒙地卡羅模擬法, 逆高斯機率分配, 平行處理

(3)

Abstract

Financial engineering and financial innovation flourished in last decades. We have developed many new financial products to provide hedge instruments for risk management, and promoted market efficiency and completeness. The pricing problems of this financial field will try to build mathematical models and derive analytic pricing formulas. But most exotic derivatives are too complicated to derive formulas. We must use computers to handle numerical methods and simulations, so computer science can give them a favor. This thesis discusses pricing of Parisian options and includes a lot of subjects: financial theory, probability & statistics, discrete mathematics, computational complexity, design & analysis of algorithms, and parallel processing.

Parisian options are path-dependent options and their closed-form solutions are not available up to now. We propose two fast financial algorithms to solve it. First we price Parisian options based on a combinatorial approach in binomial tree by Costabile in 2002. To refine Costabile’s algorithm, time complexity ^{O n}

( )

³ ^{can be}

reduced to ^{O n}

( )

² ; If ignore binomial coefficients are given, the space complexity

( )

²

O n could be reduced to ^{O n}

( )

. Second on Monte Carlo simulation, we introduce the inverse Gaussian distribution and its sampling method. To combine simulations and the inverse Gaussian distribution sampling, it can reduce divided time intervals to save computational time. Because the paths generated by Monte Carlo simulation are independent, it is easy to apply parallel processing. Nowadays multi-core processors are very popular, it is also a good idea to enhance computational efficiency. We give some descriptions and applications on it.

All financial algorithms in this thesis are implemented by the C programming language. We execute the programs in our high-performance computing clustered platform, and deal with simulation jobs synchronously. Then the system can be fully exploited.

Keywords: Parisian options, barrier options, option pricing, algorithm, binomial tree model, combinatorial method, Monte Carol simulation, inverse Gaussian distribution, parallel processing

(4)

List of Figures

1.1 Payoff of options

1.2 A sample underlying asset’s price path of a down-and-in option

2.1 One time step in binomial tree 2.2 Binomial tree

3.1 Partitions of the areas of computation

3.2 Costabile’s algorithm for up-and out Parisian options 3.3 Improved algorithm for up-and out Parisian options

4.1 Crude Monte Carlo simulation for up-and-in Parisian options 4.2 First passage time

4.3 Michael’s sampling method for the inverse Gaussian distribution 4.4 Fast simulation algorithm for up-and-in Parisian options

4.5 Fast simulation algorithm for up-and-in Parisian options with the Black-Scholes formula

4.6 A quad-core processor usage for one, two, three, and four single-thread processes simultaneously

(6)

List of Tables

3.1 Number of operations needed in Region C and Region D for Costabile’s algorithm

3.2 Number of operations needed in Region C and Region D for the improved algorithm

3.3 Numerical results of different window periods by Costabile

3.4 Numerical results of different window periods from Costabile’s algorithm and the improved algorithm

3.5 Numerical results of different window periods compared with the Avellaneda-Wu model

4.1 Numerical results of different window periods from the fast simulation algorithm

(7)

Chapter 1 Introduction

A derivative (or derivative security) is a financial instrument whose value depends on, or is derived from, the values of other, more basic, underlying assets such as commodities, equities, bonds, interest rates, exchange rates, indexes, or other derivatives. The main types of derivatives are futures, forwards, options, and swaps.

1.1 Options

There are two basic types of option. A call option gives the holder the right to buy the underlying asset by a certain date for a certain price. A put option gives the holder the right to sell the underlying asset by a certain price. The date specified in the contract is known as the expiration date. The price specified in the contract is known as the strike price. American and European options differ in when they can be exercised.

American options can be exercised at any time up to the expiration date, whereas European options can be excised only on the expiration date [9].

There are two sides to every option contract. On one side is the investor who has taken the long position. On the other side is the investor who has taken a short position. Let K be the strike price and S the final price of the underlying asset. _T There are four types of European option position and their payoff:

1. Long a call, payoff=^{max 0,}

(

^S^T ⁻^K

)

2. Long a put, payoff=^{max 0,}

(

^K⁻^S^T

)

3. Short a call, payoff=^{min 0,}

(

^S^T ⁻^K

)

4. Short a put, payoff=^{min 0,}

(

^K⁻^S^T

)

Figure 1.1 shows these payoffs.

(8)

Price Long a call

Price Long a put

Price Short a call

Price Short a put

K

Figure 1.1: Payoff of options.

1.2 Barrier and Parisian Options

Options whose payoff depends on whether the underlying asset’s price reaches a certain price H are called barrier options. These barrier options can be classified as either knock-out options or knock-in options. A knock-out option is like a standard European option except that it ceases to exist if barrier H is reached by the price of its underlying asset. A knock-in option, in contract, comes into existence if a certain barrier H is reached.

Let S₀ be the initial price of the underlying asset. A knock-out option is sometimes called a down-and-out option if H <S₀. A knock-out option is sometimes called an up-and-out option if H >S₀. A down-and-in option is a knock-in option that comes into existence only when the barrier is reached and H <S₀ (see Figure 1.2). An up-and-in option is a knock-in option that comes into existence only when the barrier is reached and H >S₀ [12]. Because there are two basic types of options, i.e., calls and puts, there are eight types of standard barrier option:

1. down-and-out call;

2. down-and-out put;

3. up-and-out call;

4. up-and-out put;

5. down-and-in call;

(9)

6. down-and-in put;

7. up-and-in call;

8. up-and-in put.

S0

K

H

Figure 1.2: A sample underlying asset’s price path of a down-and-in option.

Barrier options have many variations. If the barrier must be breached for a particular length of time, we have a Parisian option. Like common barrier options, it can be of the form of a down-and-out, up-and-out, down-and-in, up-and-in call or put.

An up-and-out (up-and-in) Parisian option is an option that expires (comes into existence, respectively) if the underlying asset’s price remains continuously at or above a given barrier over a pre-specified time interval, the so-called window period.

Conversely, a down-and-out (down-and-in) Parisian option expires (is activated, respectively) if the underlying asset’s price remains continuously at or below a lower barrier over the window period [6].

1.3 Thesis Structure

There are five chapters in this thesis. We introduce new financial products in this Chapter. Chapter 2 reviews the background knowledge and basic numerical techniques. Chapter 3 describes the combinatorial approach for pricing Parisian options by Costabile and how to speed up Costabile’s algorithm. Chapter 4 discusses how to use Monte Carlo simulation with the inverse Gaussian distribution for pricing Parisian options. Chapter 5 concludes.

(10)

Chapter 2 Fundamental Concepts

In this chapter, we describe pricing models and methods. It covers the Wiener process, the Black-Scholes option pricing model, the risk-neutral valuation, and the binomial options pricing model.

2.1 Wiener Process

A Wiener process is a particular type of Markov stochastic process. It is sometimes referred to as normalized Brownian motion. A process, z, follows a Wiener process if it has the following two properties:

1. The change ∆ during a small period of time tz ∆ is

z ε t

∆ = ∆ ^, where ε has a standardized normal distribution.

2. The value of ∆ for any two difference short intervals of time, tz ∆ , are independent.

That is a stochastic process where the change in a variable during each short period of time of length ∆ has a normal distribution with a mean equal to zero and a variance t equal to ∆ . t

We assume that our stock price follows the stochastic process:

dS =µS dt+σ S dz^,

where σ is the volatility of the stock price, and µ is its expected rate of return.

This equation is the most widely used model of stock price behavior and is also called the geometric Brownian motion [9].

2.2 Black-Scholes Option Pricing Model

In the year 1973, Fisher Black and Myron Scholes published the well-known option pricing model now universally known as the Black-Scholes option pricing model [2].

Several assumptions are made in the model as follows:

1. The stock price follows the process dS=µS dt+σ S dz with constant µ and σ _.

2. The short selling of securities with full use of proceeds is permitted.

3. There are no transactions costs or taxes.

4. All securities are perfectly divisible.

5. There are no dividends during the life of the derivative.

6. There are no riskless arbitrage opportunities.

(11)

7. Security trading is continuous

8. The risk-free rate of interest, r, is constant and the same for all maturities.

Black and Scholes derived the following formula for the price at time 0 of a European call and a European put:

( ) ( )

0 1 2

2 0 1

rT

C S N d Ke N d

P Ke N d S N d

−

= −

= − − − ^,

where

( ) ( )

( )

2 0

1

2 0

2 1

0

ln 2

cumulative probability distribution function for a standardized normal distribution stock price at time zero

strike price

continuously compounded risk-fr

S K r T

d T

S K r T

d d T

T N x

S K

r

σ σ

σ

= + +

= + − = −

=

= ee interest rate

stock price volatility

time to maturity of the option T

σ =

=

When the stock provides a continuous dividend yield at rate q , we obtain the price of a European call and a European put as, respectively,

( ) ( )

0 1 2

2 0 1

qT rT

rT qT

C S e N d Ke N d

P Ke N d S e N d

− −

= −

= − − − ^,

and the parameters d₁ and d₂ are given by

( ) ( )

2 0

1

2 0

2 1

ln 2

S K r q T

d T

S K r q T

d d T

T σ σ

σ σ

σ + − +

=

+ − −

= = −

.

These formulas, due to Merton [13], remain valid even if the dividend yield is not constant during the life of the option as long as q is replaced by the average annualized dividend yield during the life of the option.

2.3 Risk-Neutral Valuation

Risk-neutral valuation of derivatives assumes in a risk-neutral world where investors

(12)

are assumed to require no extra return on average for bearing risks. It gives the correct price for a derivative in all worlds, not just in a risk-neutral world. This means that for valuation purposes we can use the following procedure:

1. Assume that the expected return from the underlying asset is risk-free interest rate, r .

2. Calculate the expected payoff from the derivative.

3. Discount the expected payoff at risk-free interest rate.

Suppose the asset provides a yield of q. The expected return in the form of capital gain must be r−q.

In a risk-neutral world, the current price of a derivative then is given by

(

⁰

)

Price=e⁻^rTE f S ,,S_T ,

where T is the maturity date of derivative, ^{f S}

(

⁰^,^,^S^T

)

is the derivative’s payoff, which may be dependent on entire price history of underlying asset, and S₀,,S_T is the history of prices for the underlying asset from t= to T . 0

2.4 Binomial Option Pricing Model

John Cox, Stephen Ross and Mark Rubinstein developed the original version of the binomial option pricing model (BOPM) [7]. The BOPM is a binomial tree (lattice) algorithm. The binomial model is a discrete time and discrete variable pricing model.

It suffers from two unrealistic assumptions:

1. The stock price takes on only two values in a time step.

2. Trading occurs at discrete points in times.

f

fu

fd

S

Su

Sd p

1−p

∆t

Figure 2.1: One time step in binomial tree.

First we assume the stock price follows the process dS=µS dt+σ S dz^{. For} brevity, we use S for the current stock price here. Consider the stock at time

(13)

t T n

∆ ≡ , where T is the time to maturity, and n is the number of time steps. It implies that

( )

^{2 2}

(

² ¹

)

² ²

t

t t

E S t Se

Var S t S e e S t

µ

µ σ σ

∆

∆ ∆

 ∆ =

 

 ∆ = − → ∆

 

.

Second we assume that we are working in a risk-neutral world, so µ = . As in r Figure 2.1, it limits the stock moves from its current price S to one of two new values, Su and Sd , in each time step. It will increase to Su with probability p or decrease to Sd with probability 1−p. The expected stock price at time ∆ t converges to Se^{r t}^∆ . Then we get:

(

¹

)

^{r t}

pSu+ − p Sd =Se ^∆ .

The variance of a variable X is defined as ^{E X}

( )

² ^{− }^^{E X}

^{( )}

^^² [28]. The variance of the stock price at time ∆ converges to t S²σ²∆t. Then we get:

( ) (

² ¹

)( )

²

( )

^{r t} ² ² ²

p Su + − p Sd − Se^∆ =S σ ∆ t or

( )

2 2 2 2

1 ^{r t}

pu + − p d −e ^∆ = ∆ . σ t Cox, Ross and Rubinstein (CRR) model imposes

u 1

=d , we can get the solution

t

r t

u e d e

e d

p u d

σ σ

∆

− ∆

∆

=

= −

− .

The option on the stock whose current price is f . When stock price moves up to Su , the option price is f ; when stock price moves down to Sd , the option price u

is f . The risk-neutral valuation gives _d

( )

(

¹

)

r t

u d

f =e^{− ∆} pf + − p f .

Using the two values of the option in the next time step, we can evaluate the option price at a given node. This procedure is called backward induction because it works backward in time.

(14)

S0

∆t

T

Figure 2.2: Binomial tree.

Consider binomial trees with more than one time step (see Figure 2.2). Note that the tree recombines in the sense that an up move followed by a down move leads to the same asset price as a down move followed by an up move. Note also that u and d are the same at each node of the tree and the time steps are of the same length, so the risk-neutral probability p is the same at each node.

By risk-neutral valuation principle, the option price is equal to its expected payoff in a risk-neutral world discounted at the risk-free interest rate. The BOPM is a three-step procedure:

1. Generate the price tree.

2. Calculate the payoffs in terminal nodes.

3. Iterating backward induction from the end of a tree to its beginning, we are able to obtain the value of the option at time zero.

In Figure 2.2 there are quadratic, i.e., ^{O n}

( )

² , nodes. So the BOPM must take ^{O n}

( )

²

steps to apply backward induction by visiting every node. We can reuse a one-dimensional array of size n+1 instead of a two-dimensional array of size

( ) ( )

ⁿ+ × + for storing node values; so the memory requirement is reduced to ¹ ⁿ ¹

( )

O n .

As the number of time steps increases, the stock price ranges over ever-larger numbers of possible values, and trading takes place nearly continuously. It can be proved that the price computed by the BOPM converges to the price from the Black-Scholes option pricing model as n → ∞.

(15)

Chapter 3 Trees with Combinatorial Method in Pricing

Combinatorial methods use combinatorics to solve problems. It has been widely applied in many fields [10]. This chapter deals with using combinatorial methods to price financial options. Under the popular tree model, combinatorial methods can often speed up the pricing of European options by an order of magnitude [12].

3.1 Combinatorial Methods

Once we have the probability distribution of the terminal nodes, these options can be priced by simply summing the products of the probability and the payoff function at each terminal node to obtain their expected payoff. The theoretical price then equals the discounted expected payoff. We give some examples below.

In the n time step case, the value of a European call is

( ) (

⁰

)

0

1 max 0,

n n j

rT j j n j

j

C e n p p S u d K

j

− − −

=

=    − × −

∑

  ^.

The number of paths from S₀ to the terminal price S u d₀ ^j ^{n j}⁻ is given by the

binomial coefficient n j

  

 ^{, where} n is the number of time steps used for price computation and j is the number of up moves of the underlying asset’s price during the option’s lifetime. Each path to the terminal price S u d₀ ^j ^{n j}⁻ has the same

probability ^p^j

(

¹⁻^p

)

^{n j}⁻ . To adapt the equation to price European puts, simply replace the call payoff, ^{max 0,}

(

^{S u d}⁰ ^j ^{n j}⁻ ⁻^K

)

^, ^with ^the ^put ^payoff,

(

⁰

)

max 0,K−S u d^j ^{n j}⁻ . The resulting formula becomes

( ) (

⁰

)

0

1 max 0,

n n j

rT j j n j

j

P e n p p K S u d

j

− − −

=

=    − × −

∑

  ^.

Let a be the minimum number of upward price moves for the call to finish in the money, i.e., a is the smallest nonnegative integer such that S u d₀ ^a ^{n a}⁻ ≥ , or K

(16)

( )

⁰

ln ln

K S dn

a u d

 

 

=  .

Hence,

(

¹

) (

⁰

)

n n j

rT j j n j

j a

C e n p p S u d K

j

− − −

=

=    − −

∑

  ^.

Also, let b be the maximum number of upward price moves for the put to finish in the money, i.e., b is the largest integer such that S u d₀ ^b ^{n b}⁻ ≤ , or K

( )

⁰

ln ln

K S dn

b u d

 

=  .

Hence,

( ) (

⁰

)

0

1

b n j

rT j j n j

j

P e n p p K S u d

j

− − −

=

=    − −

∑

  ^.

It implies a linear-time, ^{O n}

( )

, and constant-space, ^O

( )

¹ , algorithm for pricing standard European options [12].

Consider the down-and-in call with barrier H < . We have to calculate the K number of admissible paths that lead to a particular terminal node. Towards that end, we pick the barrier level, h to be one of the n+ possible terminal underlying 1 asset’s prices that is closed to, but does not exceed, H . We call it the effective barrier.

This condition S u d₀ ^h ^{n h}⁻ ≤ implies H

( )

⁰

ln ln

H S dn

h u d

 

=  .

We effectively use S u d₀ ^h ^{n h}⁻ not H as the barrier on the binomial model. The option value equals

( ) ( )

2

1 0

2

h n j

rT j j n j

j a

e n p p S u d K

n h j

− − −

=

 

− −

 − + 

 

∑

^.

Not all n give acceptable numerical results [4]. For convergence reasons, the effective barrier level should be close to H, i.e., H ≈S d₀ ^j =S e₀ ⁻^j^σ ^{T n} for some integer j. But the effective barrier level S d₀ ^j corresponds to a terminal underlying

(17)

asset’s price only when n− j is even. To close this gap, we decrement n by one, if necessary, to make n− j an even number. The preferred n’s are finally

( )

2 2

2 0

, if is even

1, otherwise , ln

d d j j T

n d

d S H

 σ 

 −

= − ≡  ^.

It also implies a ^{O n}

( )

^-time, ^O

( )

¹ -space algorithm [11].

3.2 Technique for Handling Large Numbers

When the life of the option is divided into ever more time steps, the number of underlying asset’s price paths increases exponentially. So the combinatorial method rapidly generates large numbers. The primitive data types of typical programming languages cannot store numbers larger than a certain magnitude. For example, the largest number the data type double can represent is 1.7976931348623157 10× ³⁰⁸^. The problem that occurs when the value to be represented is greater in magnitude is called overflow. To solve this problem, we will process data in the log domain. We calculate ^ln

(

^p¹^⋅^p²

)

^from ^ln^p¹^and ^ln ^p² as follows:

(

¹ ²

)

¹ ²

ln p _⋅p ₌lnp ₊lnp .

We calculate ^ln

(

^p¹⁺^p²

)

^from ^ln^p¹^and ^ln^p² as follows:

( )

( ) (

² ¹

)

2

1 2 1

1

2 1

1 ln ln 1

ln ln 1

ln ln 1 ^p ^p

p p p p

p p p

p

p e ⁻

  

+ =   + 

 

= +  + 

 

= + +

.

When the number in the decimal domain is needed, simply take the exponential function of the final number for the desired number. In the following algorithms, we will assume the large numbers technique implicitly.

3.3 Costabile’s Algorithm

To evaluate Parisian options, Costabile proposed an algorithm in the CRR model and applied the combinatorial method [6]. In this section, we describe Costabile’s algorithm for European up-and-out Parisian options.

We consider a particle whose dynamics is described by a random walk in the two

(18)

dimensional space

( )

i j . From ^,

( )

i j , it moves to ^,

(

ⁱ⁺^1,^j+ with probability p ¹

)

if an up move occurs or to

(

ⁱ⁺^1, ^j− with probability ¹

)

1− p if a down move takes place. Without loss of generality we assume that the particle starts from the origin

( )

^{0, 0}

P≡ .

The number of paths, N i j , from the origin P to

( )

^,

( )

i j has the ^, well-known recursive relation, ^{N i j}

( )

^, ⁼^{N i}

(

⁻^1,^j^{− +}¹

) (

^{N i}⁻^1,^j^{+ .}¹

)

N i j is

( )

^,

given by the binomial coefficient

( ) 2

i i j

 

 + 

 . Suppose the barrier level in the binomial tree is at level m. Let g i j denote the number of the paths from P to

( )

^,

( )

i j so that the paths spend fewer than l consecutive time steps at or above the ^, barrier level, l∈ », which is an input parameter. Such paths are said to be admissible.

Later we will describe how to derive m from the barrier H and l from the window period w .

(19)

l+ ≤ ≤m j n

m≤ < +j l m

n j m

− ≤ <

0≤ < +i l m l+ ≤ ≤m i n m

( ,n−n) ( , )n n

(0, 0) P

D B

A C

Figure 3.1: Partitions of the areas of computation. This figure is due to Prof. Jr-Yan Wang.

We derive ^{g i j}

( )

^, ^for

( )

^{i j}^, in different regions of the tree (see Figure 3.1).

The first step considers the nodes

( )

^{i j}^, such that 0≤ < + called Region A. In i l m this case ^{g i j}

( )

^, ⁼^{N i j}

( )

^, . Indeed, such nodes are characterized by at most i up moves and i m− < . In the second step we observe that, for the nodes l

( )

i j such ^, that l+ ≤ ≤m i n l, + ≤ ≤m j n (Region B), ^{g i j}

( )

^, = . Indeed, in such cases, a path ⁰ starting form P can reach the level l+ ≤m j only if the particle spends at least l successive time steps at or above the barrier level. The third step concerns the nodes

( )

i j located below the barrier level (i.e., ^, l+ ≤ ≤ − ≤ < +m i n, n j l m or Region C).

Each path from P to

( )

i j in Region C makes its last move below the barrier level. ^,

(20)

This implies ^{g i j}

( ) (

^, ⁼^{g i}⁻^1,^j^{− +}¹

) (

^{g i}⁻^1, ^j^{+ .}¹

)

In the final step we are left to calculate g i j of the nodes

( )

^,

( )

i j such that ^, ,

l+ ≤ ≤m i n m≤ ≤ +j l m or Region D. Note that in this case the recursive relation

( ) (

^, ^1, ¹

) (

^1, ¹

)

g i j =g i− j− +g i− j+ is no longer holds. Indeed some of the

(

^1, ¹

)

g i− j+ paths from P to

(

ⁱ⁻^1,^j+ may have spent ¹

)

l − time steps at or 1 above the barrier level already, in which case the move from time step i − to i is 1 made above the barrier level. Such paths have to be excluded from ^{g i}

(

⁻^1,^j^{+ .}¹

)

Each admissible path from P to

( )

i j in Region D can be divided into two ^, different parts. The first part is an admissible path from P to a node which is located at level m − and then makes an up move to a node 1

(

ⁱ⁻²^k⁻

(

^j⁻^{m m}

)

^,

)

^{on the}

barrier level. The number of the first partial path is given by the quantity

( )

(

² ^1, ¹

)

g i− k− j−m − m− . The second part is a partial path from the final node of

first partial path,

(

ⁱ⁻²^k⁻

(

^j⁻^{m m}

)

^,

)

, to the node

( )

i j with the number of moves ^, strictly smaller than l , and this part should be never cross the barrier level. The parameter k should be a nonnegative integer and the condition

( )

(

²

)

i− − − −i k j m <l implies:

( )

0 2

l j m

k − −

≤ < ^.

The number of the second partial paths can be computed as the difference between the following two terms:

1. The first is the number of paths from the node

(

ⁱ⁻²^k⁻

(

^j⁻^{m m}

)

^,

)

^{on the}

barrier level to

( )

i j . It is given by ^, 2k j m k j m

 + − 

 + − 

 ^.

2. The second is the number of paths from the same node to

( )

^{i j}^, which cross the

(21)

barrier level. Using the reflection principle, it is given by 2 1 k j m

k j m

 + − 

 + + − 

 ^with

the convention a 0 b

 =

   ^for a_{< .}b

Finally we get the number of admissible paths as

( )

0 2

2 2

( , ) 2 1, 1

l j m 1

k

k j m k j m

g i j g i k j m m

k j m k j m

≤ <− −

 + −   + − 

=

∑

− − − − − ⋅ + −   − + + − ^{. (1)}

See [6] for additional discussions.

We assume that the underlying asset’s price dynamics is described by the CRR model and also use the notations as before. Recall that w denotes the window period measured in years and it must be limited to 0≤ ≤ . When a Parisian option with w T

0

w= , this is a standard barrier option. When w= , we have either a standard T European option or a option with zero value.

In order to drive the binomial algorithm for evaluating Parisian options, we need to determine the total number of time steps, n, to be used in the price computations.

As pointed out by Boyle and Lau, at first we build up a binomial tree such that the barrier is close to but just below a layer of horizontal nodes in the tree [4]. Recall that m is the number of successive up moves for the underlying asset’s price, starting from S₀, must make to touch or cross the barrier H. This condition S u₀ ^m≥ H implies

( )

2 2

2

ln 0

m T

n H S

 σ 

=  

 

 , and

(

⁰

)

ln H S

m σ t

 

=  

 ∆ .

The next step is to determine the number of time steps, l , corresponding to the window period w. Because, in general, w ∆t is not an integer, we set l equal to

[

^w ∆ , the integer closest to ^t

]

w ∆t. The price of the up-and-out Parisian call option can now be easily derived:

(

^{, 2}

) (

¹

) (

⁰

)

n n j

rT j j n j

j a

PC e⁻ g n j n p p ⁻ S u d ⁻ K

=

∑

− − − ^.

(22)

Recall that ^{g n}

(

^{, 2}^j⁻ⁿ

)

is the number of paths which spend fewer than l consecutive time steps at or above the barrier level and these paths are made up of j up moves and n− j down moves.

Together with the initial conditions:

0 1 1

n n

   

= =

   

    ^,

Pascal’s triangle gives a method for calculating the binomial coefficients by using the following well-known identity:

1

n n n

k k k

 +     

= +

   −   

     

where n and k are positive integers with n≥ [17]. The binomial coefficients k can be computed in ^{O n}

( )

² time. Assume the binomial coefficients have been computed from now on. See Figure 3.2 for Costabile’s algorithm.

1: ∆ =t: T n^;

2: ^u^{: exp}⁼

( )

^σ ^∆^t ^;

3: d: 1= u^;

4: _p_: ^exp

( (

^r ^q

)

^t

)

^d

u d

− ∆ −

= − ^;

5: _m_: ^ln

(

^{H S}⁰

)

σ t

 

=  

 ∆ ;

6: ^l^:⁼

[

^w ^{∆ ;}^t

]

7: ^a^:^{= }^^ln

(

^{K S d}⁰ ⁿ

)

^ln

^{( )}

^{u d} ^^^;

8: for : 0i _{= to} n do

9: for j:=i^to − ⁱ step − 2 do 10: if i< + l m then {Region A}

11: ^{g i}

[ ][ ]

^j ^:^{= }^_₍_i₊ⁱ_j_{) 2}^^_^;

12: elseif j≥ +l m then {Region B}

13: ^{g i}

[ ][ ]

^j ^{: 0}^{= ;}

(23)

14: elseif j< then {Region C} m 15: if j= − then i

16: ^{g i}

[ ][ ] [ ][ ]

^j ^:⁼^{g i}⁻¹ ^j^{+ ;}¹

17: else

18: ^{g i}

[ ][ ] [ ][ ] [ ][ ]

^j ^:⁼^{g i}⁻¹ ^j^{− +}¹ ^{g i}⁻¹ ^j^{+ ;}¹

19: endif

20: else {Region D}

21: ^{g i}

[ ][ ]

^j ^{: 0}^{= ;}

22: for k: 0= to

( )

_0.5

2

l j m

 − − 

 − 

  do

23: if 2k+ − < + + −j m k 1 j m then

24: ^{g i}

[ ][ ]

^j ^:⁼^{g i}^ ⁻²^k⁻

(

^j⁻^m

)

⁻¹^

[

^m^{− ⋅}¹

]

^²_k^k+ −^{+ −}_j^j _m^m^^^; 25: else

26: ^{g i}

[ ][ ]

^j ^:⁼^{g i}^ ⁻²^k⁻

(

^j⁻^m

)

⁻¹^

[

^m^{− ⋅}¹

]

^^^^²_k^k+ −^{+ −}_j^j _m^m^{ }^{ } ⁻ _k²+ + −^k₁^{+ −}^j_j ^m_m^^^^; 27: endif

28: endfor 29: endif 30: endfor 31: endfor 32: C: 0_{= ;}

33: for j:=a^to n do

34: ^C^:^{= +}^C ^{g n}

[ ][

²^j^{− ⋅}ⁿ

]

^p^j

(

¹⁻^p

)

^{n j}⁻ ^⋅

(

^{S u d}⁰ ^j ^{n j}⁻ ⁻^K

)

^;

35: endfor

36: return C e⋅ ⁻^rT^;

Figure 3.2: Costabile’s algorithm for up-and-out Parisian options.

We proceed to count the number of operations needed to compute the ^{g i j}

( )

^, ^.

For simplicity, the number of time steps, n, is odd for the option evaluation and the parameters m and l are even. Because the ^{g i j}

( )

^, in Region A and Region B

(24)

needed no arithmetic operations, Table 3.1 only shows the number of mathematical operations needed to implement Costabile’s algorithm in Region C and Region D. See [6] for more details. Recall that l was defined as the integer closed to w ∆ = ⋅t w n T. It follows that the computational cost is a cubic function of n, or

( )

³

O n . The numbers ^{g i j}

( )

^, are stored in a two-dimensional array of size

( ) (

ⁿ^{+ ×}¹ ²ⁿ+ so the memory requirement is ¹

)

^{O n}

( )

² ^.

Region C Region D

Additions

( )

² ² ² ¹

(

²

)

²

4

n+m + m+ + − +n l m ¹

( )

₁

2 2 2

n l m l l

 + − +  −  

    

 

Subtractions 0

( )

1 1

2 2 2

n l m l l

 + − +   + 

    

 

Multiplications 0

( )

1 1

2 2 2

n l m l l

 + − +   + 

    

 

Table 3.1: Number of operations needed in Region C and Region D for Costabile’s algorithm.

3.4 An Improved Algorithm

For nodes

( )

^{i j}^, in Region D, ^{g i j}

( )

^, depends on the ^{g i}

(

⁻²^k⁻

(

^j⁻^m

)

⁻^1,^m⁻¹

)

located at level m− in Region A and Region C (recall Equation (1)). So we need 1 not calculate all ^{g i j}

( )

^, in Region D. In Region C when j=m−1, the recursive relation ^{g i j}

( ) (

^, ⁼^{g i}⁻^1,^j^{+ +}¹

) (

^{g i}⁻^1, ^j+ is also dependent on the ¹

)

^{g i}

(

⁻^1,^m

)

located at level m. So we only need to calculate the ^{g i j}

( )

^, in Region D which are 1. nodes located at level m, and

2. nodes located at the n’th time step.

It is not hard to observe that other ^{g i j}

( )

^, are not needed for calculating

[ ][

²

]

g n j₋n of Line 34 in Figure 3.1.

(25)

The memory requirement can be reused if the space is reused. Specifically, replace the two-dimensional array ^{g n}

[ ][

⁺^{1 2}ⁿ+ with a one-dimensional array of ¹

]

size 2n+ . But, we should use another one-dimensional array 1 ^{t n}

[ ]

+ to store a ¹ copy of ^{g i m −}

(

^, ¹

)

located at level m − for calculating 1 g i j in Region D.

( )

^,

These crucial points can be used to improve Costabile’s algorithm, and the results must be identical to Costabile’s algorithm. See Figure 3.3 for the improved algorithm.

1: ∆ =t: T n^;

2: ^u^{: exp}⁼

( )

^σ ^∆^t ^;

3: d: 1= u^;

4: _p_: ^exp

( (

^r ^q

)

^t

)

^d

u d

− ∆ −

= − ^;

5: _m_: ^ln

(

^{H S}⁰

)

σ t

 

=  

 ∆ ;

6: ^l^:⁼

[

^w ^{∆ ;}^t

]

7: ^a^:^{= }^^ln

(

^{K S d}⁰ ⁿ

)

^ln

^{( )}

^{u d} ^^^;

8: for : 0i = to n do

9: for j:=i^to − ⁱ step − 2 do 10: if i< + l m then {Region A}

11: ^{g j}

[ ]

^:^{= }^_₍_i₊ⁱ_j_{) 2}^^_^;

12: elseif j≥ +l m then {Region B}

13: ^{g j}

[ ]

^{: 0}^{= ;}

14: elseif j<m then {Region C}

15: if j= −i then 16: ^{g j}

[ ] [ ]

^:⁼^{g j}^{+ ;}¹

(26)

17: else

18: ^{g j}

[ ] [ ] [ ]

^:⁼^{g j}^{− +}¹ ^{g j}^{+ ;}¹

19: endif

20: else if j=m^{or i}= ⁿ then {Partial Region D}

21: ^{g j}

[ ]

^{: 0}^{= ;}

22: for k: 0= to

( )

_0.5

2

l j m

 − − 

 − 

  do

23: if 2k+ − < + + −j m k 1 j m then

24: ^{g j}

[ ]

^:⁼^{t i}^ ⁻²^k⁻

(

^j⁻^m

)

^{− ⋅}¹^ ^²_k^k+ −^{+ −}_j^j _m^m^^; 25: else

26: ^{g j}

[ ]

^:⁼^{t i}^ ⁻²^k⁻

(

^j⁻^m

)

^{− ⋅}¹^ ^^^^²_k^k+ −^{+ −}_j^j _m^m^{ }^{ } ⁻ _k²+ + −^k₁^{+ −}^j_j ^m_m^^^^; 27: endif

28: endfor 29: endif

30: if j=m−1 then 31: ^{t i}

[ ] [ ]

^:⁼^{g j} ^;

32: endif 33: endfor 34: endfor 35: C: 0= ;

36: for j:=a^to n do

37: ^C^:^{= +}^C ^g

[

²^j^{− ⋅}ⁿ

]

^p^j

(

¹⁻^p

)

^{n j}⁻ ^⋅

(

^{S u d}⁰ ^j ^{n j}⁻ ⁻^K

)

^;

38: endfor

39: return C e⋅ ⁻^rT^;

Figure 3.3: Improved algorithm for up-and-out Parisian options.

Recall that we assume the parameters n is odd, and m and l are even. We observe that there are ^^ⁿ^{+ − +}¹

(

^l ^m

)

^^ ² nodes at each level in Region D. According to Equation (1), at levels m and m+ , each 1 g i j is calculated with

( )

^, l 2 1−

(27)

additions, l 2 multiplications and l 2 subtractions. At levels m+ and 2 m+ , 3 each g i j is calculated with

( )

^, l 2 2− additions, l 2 1− multiplications and

2 1

l − subtractions. Continuing in this manner we can compute the total number of additions needed to compute g i j located at level

( )

^, m and at the n’th time step in Region D as

( )

2 1

1

located at the 'th time step located at level

1 2 2 2

1 1

2 2 2 4

l

j m n

n l m

l l n m l

j

−

=

 + − +  − + −

 −  + = −  

    

  

∑

  

.

The number of multiplications and subtractions needed is

( )

2

1

located at the 'th time step located at level

1 2 2 4

2 2 2 4

l

j m n

n l m

l l n m l

j

=

 + − +  − + −

  + =  

    

  

∑

  

.

Table 3.2 shows the number of mathematical operations needed to implement the improved algorithm in Region C and Region D.

The computational cost of the improved algorithm is less than Costabile’s algorithm, which is ^{O n}

( )

³ . The reduction of nodes that need to be computed in Region D leads to ^{O n}

( )

² computational time. This improved algorithm could be

also speeded up somewhat, but the total time complexity remains ^{O n}

( )

² . Because binomial coefficients are given and we only use two one-dimensional arrays for the improved algorithm, the memory requirement is ^{O n .}

( )

Region C Region D

Additions

( )

² ² ² ¹

(

²

)

²

4

n+m + m+ + − +n l m 2 2 2

2 1 4

l n− m+ −l

 −  

  

  

Subtractions 0 2 2 4

2 4

l n− m+ −l

  

  

  

Multiplications 0 2 2 4

2 4

l n− m+ −l

  

  

  

Table 3.2: Number of operations needed in Region C and Region D for the improved algorithm.

(28)

3.5 General Case

The ^{g i j}

( )

^, of the down-and-out barrier level, −m is the same as the ^{g i}

( )

^,⁻^j ^of

the up-and-out barrier level, m with an identical window period. To price a down-and-out Parisian call with a barrier H′, it is not necessary to discuss where Region A, Region B, Region C, and Region D are. We can reuse the up-and-out Parisian options pricing algorithm in Figure 3.3 and replace Line 5 with

(

⁰

)

: ln H S

m σ t

 ′ 

= −  

 ∆ ;

and replace Line 37 with

( ) ( ) (

⁰

)

: 2 ^j 1 ^{n j} ^j ^{n j}

C = +C g− j−n ⋅p −p ⁻ ⋅ S u d ⁻ −K ^; Then all knock-out Parisian calls could be priced.

The ^{g i j}

( )

^, of knock-in options is the same as ^{N i j}

( )

^, subtract the ^{g i j}

( )

^,

of knock-out options with identical barrier and window period. For example, to price an up-and-in Parisian call, replace Line 37 with

[ ] ( ) (

⁰

)

: n 2 ^j 1 ^{n j} ^j ^{n j}

C C g j n p p S u d K

j

− −

  

= + − − ⋅ − ⋅ −

   ;

Use the BOMP in Section 2.4 or combinatorial method’s formula for the standard European call in Section 3.1 to subtract the up-and-out Parisian call, it can also get equal value. It is just the in-out parity. Then all Parisian calls could be priced.

To price a Parisian put, refer to Section 3.1 for b , the maximum number of upward price moves for the put to finish in money. For example, to price an up-and-out Parisian put, replace Line 7 with

(

⁰

) ^{( )}

: ln ⁿ ln

b =  K S d u d ^; Line 35 as

: 0 P _{= ;} Line 36 as

for j: 0= ^to b do and Line 37 as

[ ] ( ) (

⁰

)

: 2 ^j 1 ^{n j} ^j ^{n j}

P = +P g j− ⋅n p −p ⁻ ⋅ K−S u d ⁻ ^;

國立臺灣大學電機資訊學院資訊工程學系 碩士論文