and

12

¯ {k var(Y (T /k)) + (k ’ 1)2 var{Y (2T /k) ’ Y (T /k)} + ...}

var(Y ) = 2

k

k

T σ 2 X 2 T σ2 (k + 1)(2k + 1)

=3 i=

6k 2

k i=1

= σ 2 T, say.

e

The closed form solution for the price E(V2 ) in this case is therefore easily

obtained because it reduces to the same integral over the lognormal density

that leads to the Black-Scholes formula. In fact

E(V2 ) = E{e’rT (S0 eY ’ K)+ }, where Y ∼ N (e, σ 2 T ) so

µe

= E[e’rT +eT S0 eY ’eT ’ e’rT K]+

µ µ

1

2

1

e

= E[S0 e(’r+e+ 2 σ

µ )T

exp{Y ’ µT ’ σ 2 T } ’ e’rT K]+

e e

2

σ2T 2

e

2

1

e

= E[S0 e(’r+e+ 2 σ

µ )T

, σ T )} ’ Ke’rT ]+ .

e

exp{N (’

2

where we temporarily denote a random variable with the Normal(µ, σ 2 ) distrib-

ution by N (µ, σ 2 ). Recall that the Black-Scholes formula gives the price at time

t = 0 of a European option with exercise price K, initial stock price S0 ,

σ2

BS(S0 , K, r, T, σ) = E(e’rT (S0 exp{N ((r ’ )T, σ 2 T )} ’ K)+ (5.2)

2

σ2T 2

, σ T )} ’ Ke’rT )+ (5.3)

= E(S0 exp{N (’

2

= S0 ¦(d1 ) ’ Ee’rT ¦(d2 ) (5.4)

258 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

where

√

log(S0 /K) + (r + σ 2 /2)T

√ , d2 = d1 ’ σ T .

d1 =

σT

Thus E(V2 ) is given by the Black-Scholes formula with S0 replaced by

e2 2

f = S0 exp{T ( σ + µ ’ r)} = S0 exp{’rT (1 ’ 1 ) ’ σ T (1 ’ 1 )}

e

S0

k2

2 k 12

and σ 2 by σ 2 . Of course when k = 1, this gives exactly the same result as the

e

basic Black-Scholes because in this case, the Asian option corresponds to the

average of a single observation at time T . For k > 1 the e¬ective initial stock

f

price is reduced S0 < S0 and the volatility parameter is also smaller σ 2 < σ 2 .

e

With lower initial stock price and smaller volatility the price of a European call

will decrease, indicating that an Asian option priced using a geometric mean

has price lower than a similar European option on the same stock.

Recall from our discussion of a control variate estimators that we can esti-

mate E(V1 ) unbiasedly using

V1 ’ β(V2 ’ E(V2 )) (5.5)

where

cov(V1 , V2 )

(5.6)

β= .

var(V2 )

In practice, of course, we simulate many values of the random variables V1 , V2

and replace V1 , V2 by their averages V1 , V 2 so the resulting estimator is

V1 ’ β(V 2 ’ E(V2 )). (5.7)

Table 4.1 is similar to that in Boyle, Broadie and Glasserman(1997) and com-

pares the variance of the crude Monte Carlo estimator with that of an estimator

using a simple control variate,

E(V2 ) + V1 ’ V 2 ,

a special case of (5.7) with β = 1. We chose K = 100, k = 50, r = 0.10, T = 0.2,

a variety of initial asset prices S0 and two values for the volatility parameter

ASIAN OPTIONS 259

σ = 0.2 and σ = 0.4. The e¬ciency depends only on S0 and K through the

ratio K/S0 or alternatively the moneyness of the option, the ratio erT S0 /K of

the value on maturity of the current stock price to the strike price. Standard

errors are estimated from n = 10, 000 simulations. Since the e¬ciency is the

ratio of the number of simulations required for a given degree of accuracy, or

alternatively the ratio of the variances, this table indicates e¬ciency gains due

to the use of a control variate of several hundred. Of course using the control

variate estimator (5.7) described above could only improve the e¬ciency further.

Table 4.1. Standard Errors for Arithmetic Average Asian Options.

STANDARD ERROR STANDARD ERROR

Moneyness=erT S0 /K

σ

USING CRUDE MC USING CONTROL VARIATE

0.2 1.13 0.0558 0.0007

1.02 0.0334 0.00064

0.93 0.00636 0.00046

0.4 1.13 0.105 0.00281

1.02 0.0659 0.00258

0.93 0.0323 0.00227

The following function implements the control variate for an Asian option

and was used to produce the above table.

function [v1,v2,sc]=asian(r,S0,sig,T,K,k,n)

% computes the value of an asian option V1 and control variate V2

% S0=initial price, K=strike price

% sig = sigma, k=number of time increments in interval [0.T]

% sc is value of the score function for the normal inputs with respect to

% r the interest rate parameter. Repeats for a total of n simulations.

v1=[]; v2=[]; sc=[]; mn=(r-sig^2/2)*T/k;

sd=sig*sqrt(T/k); Y=normrnd(mn,sd,k,n);

260 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

sc= (T/k)*sum(Y-mn)/(sd^2); Y=cumsum([zeros(1,n); Y]);

S = S0*exp(Y); v1= exp(-r*T)*max(mean(S)-K,0);

v2=exp(-r*T)*max(S0*exp(mean(Y))-K,0);

disp([™standard errors ™ num2str(sqrt(var(v1)/n)) ™ num2str(sqrt(var(v1-v2)/n))])

For example if we use K = 100, we might con¬rm the last row of the above

table using the command

asian(.1,100*.93*exp(-r*T),.4,.2,100,50,10000);

Asian Options and Strati¬ed Sampling

For many options, the terminal value of the stock has a great deal of in¬‚uence

on the option price. Although it is di¬cult in general to stratify samples of

stock prices, it fairly easy to stratify along a single dimension, for example the

dimension de¬ned by the stock price at time T. In particular we may stratify

the generation of

St = S0 exp(Zt )

where Zt can be written in terms of a standard Brownian motion

Zt = µt + σWt , with µ = r ’ σ 2 /2.

To stratify into K strata of equal probability for ST we may generate ZT using

p Ui

rT ’ σ2 T /2 ¦’1 (i ’ 1 + ), i = 1, 2, ...K

ZT = rT +

K

for Uniform[0,1] random variables Ui and then randomly generate the rest of the

path interpolating the value of S0 and ST using Brownian Bridge interpolation.

To do this we use the fact that for a standard Brownian motion and s < t < T

we have that the conditional distribution of Wt given Ws , WT is normal with

mean a weighted average of the value of the process at the two endpoints

T ’t t’s

Ws + WT

T ’s T ’s

ASIAN OPTIONS 261

and variance

(t ’ s)(T ’ t)

.

T ’s

Thus given the value of ST (or equivalently the value of W (T )) the increments

of the process at times µ, 2µ, ...N ² = T can be generated sequentially so that the

j™th increment W (jµ) ’ W ((j ’ 1)µ) conditionally on the value of W ((j ’ 1)µ)

and of W (T ) has a Normal distribution with mean

N ’j j

W ((j ’ 1)µ) + W (T )

(

N N

and with variance

N ’j

.

N ’j+1

Use of Girsanov™s Lemma.

There are many other variance reduction schemes that one can apply to valuing

an Asian Option. However prior to attacking this problem by other methods,

let us consider a simpler example.

Importance Sampling and Pricing a European Call Option

Suppose we wish to estimate the value of a call option using Monte Carlo meth-

ods which is well “out of the money”, one with a strike price K far above the

current price of the stock S0 . If we were to attempt to evaluate this option using

crude Monte Carlo, the majority of the values randomly generated for ST would

fall below the strike and contribute zero to the option price. One possible rem-

edy is to generate values of ST under a distribution that is more likely to exceed

the strike, but of course this would increase the simulated value of the option.

We can compensate for changing the underlying distibution by multiplying by

a factor adjusting the mean as one does in importance sampling.

More speci¬cally, we wish to estimate

EQ [e’rT (S0 eZT ’ K)+ ], where ZT ∼ N (rT ’ σ 2 T /2, σ 2 T )

262 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

where EQ indicates that the expectation is taken under a risk neutral distri-

bution or probability measure Q for and K is large. Suppose that we modify

the underlying probability measure of ZT to Q0 , say a normal distribution with

mean value ln(K/S0 ) ’ σ2 T /2 but the same variance σ 2 T . Then the expected

stock price under this new distribution

EQ0 S0 eZT = S0 exp(EQ0 ZT + σ 2 T /2) = K

so there is a much greater probability (roughly 1/2) that the strike price is

attained. The importance sampling adjustment that insures that the estimator

continues to be an unbiased estimator of the option price is the ratio of two

probability densities. Denote the normal probability density function by

(x ’ µ)2

1

•(x, µ, σ 2 ) = √ }.

exp{’

2σ 2

2πσ

Then the Radon-Nikodym derivative

2

•(zt ; rT ’ σ 2T , σ 2 T )

dQ

(zT ) = 2

•(zt ; ln(K/S0 ) ’ σ 2T , σ 2 T )

dQ0

is simply the ratio of the two normal density functions with the two di¬erent

means, and

dQ

EQ [e’rT (ST ’ K)+ ] = EQ0 [e’rT (ST ’ K)+ (ZT )]

dQ0

σ2 T

, σ2T )

•(ZT ; rT ’

’rT ZT + 2

’ K)

= EQ0 [e (S0 e ]

σ2 T

, σ2 T )

•(ZT ; ln(K/S0 ) ’ 2

so the importance sample estimator is the average of terms of the form

σ2 T

, σ2T ) σ2 T 2

•(ZT ; rT ’

’rT ZT + 2

’K) , where ZT ∼ N (ln(K/S0 )’

e (S0 e , σ T ).

2

•(ZT ; ln(K/S0 ) ’ σ 2T , σ2 T ) 2

The new simulation generates paths that are less likely to produce options ex-

piring with zero value, and in a sense has thus eliminated some computational

waste. What gains in e¬ciency result from this use of importance sampling?

Let us consider a three month (T = 0.25) call option with S0 = 10, K = 15,

ASIAN OPTIONS 263

σ = 0.2, r = .05. We determined the e¬ciency of the importance sampling es-

timator relative to using crude Monte Carlo in this situation using the function

below. Running this using the command [e¬,m,v]=importance2(10,.05,15,.2,.25)

shows an e¬ciency gain of around 73, in part because very few of the crude es-

timates of ST exceed the exercise price.

function [e¬,m,v]=importance2(S0,r,K,sigma,T,N)

% simple importance sampling estimator of call option price

% outputs e¬ciency relative to crude. Run using

% [e¬,m,v]=importance2(10,.05,15,.2,.25)

Z=randn(1,N);

%¬rst do crude

ZT=(r-.5*sigma^2)*T+sigma*sqrt(T).*Z;

est1=exp(-r*T)*max(0,S0*exp(ZT)-K);

% now do importance

ZT=(log(K/S0)-.5*sigma^2)*T+sigma*sqrt(T).*Z;

ST=S0*exp(ZT);

est2=exp(-r*T)*max(0,ST-K).*normpdf(ZT,(r-.5*sigma^2)*T,sigma*sqrt(T))./normpdf(ZT,(log(K/S0)-

.5*sigma^2)*T,sigma*sqrt(T));

v=[var(est1) var(est2)];

m=[mean(est1) mean(est2)];

e¬=v(1)/v(2);

Importance Sampling and Pricing an Asian Call Option

Let us now return to pricing an Asian option. We wish to use a variety of

variance reduction techniques including importance sampling as in the above

example, but in this case the relevant observation is not a simple stock price

at one instant, but the whole stock price history from time 0 to T. An Asian

option should nevertheless have payo¬ correlated with the value of the stock on

264 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

maturity S(T ). It might be reasonable to stratify the sample; i.e. sample more

often when S(T ) is large or to use importance sampling and generate S(T )

from a geometric Brownian motion with drift larger than rSt so that it is more

™

likely that S(T ) > K. As before if we do this we need to then multiply by the

ratio of the two probability density functions, (the Radon Nikodym derivative

of one process with respect to the other). This density is given by a result called

Girsanov™s Theorem (see Appendix B). The idea is as follows: Suppose P is the

probability measure induced on the paths on [0, T ] by an Ito process

(5.8)

dSt = µ(St )dt + σ(St )dWt , S0 = s0 .

Similarly suppose P0 is the probability measure on path generated by a

similar process with the same di¬usion term but di¬erent drift term

(5.9)

dSt = µ0 (St )dt + σ(St )dWt , S0 = s0 .

Note that in both cases, the process starts at the same initial value s0 .Then the

dP

“likelihood ratio” or the Radon-Nikodym dP0 of P with respect to P0 is

ZT ZT 2

µ (St ) ’ µ2 (St )

µ(St ) ’ µ0 (St )

dP 0

dSt ’ (5.10)

= exp{ dt}

2 (S ) 2 (S )

dP0 σt 2σ t

0 0

We do not attempt to give a technical proof of this result, either here or in

the appendix, since real “proofs” can be found in a variety of texts, including

Steele (2004) and Karatzas and Shreve, (xxx). Instead we provide heuristic

justi¬cation of (5.10). Let us consider the conditional distribution of a small

increment dSt in the process St under the model (5.8). Since this distribution is

conditionally normal distributed it has conditional probability density function

given the past

1

exp{’(dSt ’ µ(St )dt)2 /(2σ 2 (St )dt)

√ (5.11)

2πdt

and under the model (5.9), it has the conditional probability density

1

exp{’(dSt ’ µ0 (St )dt)2 /(2σ 2 (St )dt)

√ (5.12)

2πdt

ASIAN OPTIONS 265

The ratio of these two probability density functions is

µ2 (St ) ’ µ2 (St )

µ(St ) ’ µ0 (St ) 0

dSt ’

exp{ dt}

2 (S ) 2 (S )

σt 2σ t

But the joint probability density function over a number of disjoint intervals is

obtained by multiplying these conditional densities together and this results in

µ2 (St ) ’ µ2 (St )

µ(St ) ’ µ0 (St ) 0

dSt ’

Πt exp{ dt}

2 (S ) 2 (S )

σt 2σ t

ZT ZT 2

µ (St ) ’ µ2 (St )

µ(St ) ’ µ0 (St ) 0

dSt ’

= exp{ dt}

2 (S ) 2 (S )

σt 2σ t

0 0

where the product of exponentials results in the sum of the exponents, or, in

the limit as the increment dt approaches 0, the corresponding integrals.

Girsanov™s result is very useful in conducting simulations because it permits

us to change the distribution under which the simulation is conducted. In

general, if we wish to determine an expected value under the measure P, we

dP

may conduct a simulation under P0 and then multiply by or if we use a

dP0

subscript on E to denote the measure under which the expectation is taken,

dP

EP V (ST ) = EP0 [V (ST ) ].

dP0

Suppose, for example, we wish to determine by simulation the expected value

of V (rT ) for an interest rate model

(5.13)

drt = µ(rt )dt + σdWt

for some choice of function µ(rt ). Then according to Girsanov™s theorem, we

may simulate rt under the Brownian motion model drt = µ0 dt + σdWt (having

the same initial value r0 as in our original simulation) and then average the

values of

ZT ZT 2

µ (rt ) ’ µ2

µ(rt ) ’ µ0

dP 0

drt ’ (5.14)

V (rT ) = V (rT ) exp{ dt}

2 2

dP0 σ 2σ

0 0

So far, the constant µ0 has been arbitrary and we are free to choose it in order

to achieve as much variance reduction as possible. Ideally we do not want to get

266 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

too far from the original process so µ0 should not be too far from the values of

dP dP

µ(rt ). In this case we hope that the term is not too variable (note that c dP0

dP0

would be the estimator if V (ST ) = c were constant). On the other hand, the

term V (rT ) cannot generally be ignored, and there is no formula or simple rule

dP

for choosing parameters which minimize the variance of V (rT ) dP0 . Essentially

dP

we need to resort to choosing µ0 to minimize the variance of V (rT ) dP0 by

experimentation, usually using some preliminary simulations.

Pricing a Call option under stochastic interest

rates.

(REVISE MODEL)Again we consider pricing a call option, but this time under

a more realistic model which permits stochastic interest rates. We will use

the method of conditioning, although there are many other potential variance

reduction tools here. Suppose the asset price, (under the risk-neutral probability

measure, say) follows a geometric Brownian motion model of the form

(1)

(5.15)

dSt = rt St dt + σSt dWt

where rt is the spot interest rate. We assume rt is stochastic and follows the

Brennan-Schwartz model,

(2)

drt = a(b ’ rt )dt + σ0 rt dWt (5.16)

(1) (2)

where Wt , Wt are both Brownian motion processes and usually assumed to

be correlated with correlation coe¬cient ρ. The parameter b in (5.16) can

be understood to be the long run average interest rate (the value that it would

converge to in the absence of shocks or resetting mechanisms) and the parameter

a > 0 governs how quickly reversion to b occurs.

It would be quite remarkable if a stock price is completely independent of

interest rates, since some of the same factors in¬‚uence both. However we begin

PRICING A CALL OPTION UNDER STOCHASTIC INTEREST RATES.267

by assuming something a little less demanding, that the random noise processes

driving the asset price and stock price are independent or that ρ = 0.

Control Variates.

The ¬rst method might be to use crude Monte Carlo; i.e. to simulate both the

process St and the process rt , evaluate the option at expiry, say V (ST , T ) and

RT

then discount to its present value by multiplying by exp{’ 0 rt dt}. However,

in this case we can exploit the assumption that ρ = 0 so that interest rates are

(1)

independent of the Brownian motion process Wt which drives the asset price

process. For example, suppose that the interest rate function rt were known

(equivalently: condition on the value of the interest rate process so that in the

conditional model it is known). While it may be di¬cult to obtain the value of

an option under the model (5.15), (5.16) it is usually much easier under a model

which assumes constant interest rate c. Let us call this constant interest rate

model for asset prices with the same initial price S0 and driven by the equation

(1)

(5.17)

dZt = cZt dt + σZt dWt , Z0 = S0

model “0” and denote the probability measure and expectations under this

distribution by P0 and E0 respectively. The value of the constant c will be

determined later. Assume that we simulated the asset prices under this model

and then valued a call option, say. Then since

σ2

)T, σ2 T ) distribution

ln(ZT /Z0 ) has a N ((c ’

2

we could use the Black-Scholes formula to determine the conditional expected

value

268 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

Z T

rt dt}(ZT ’ K)+ |rs , 0 < s < T ] (5.18)

E0 [exp{’

0

= EE0 [(S0 e(c’r)T eW ’ e’rT K)+ |rs , 0 < s < T ],

where W has a N (’σ2 T /2, σ 2 T )

Z T

1

= E[BS(S0 e(c’r)T , K, r, T, σ)], with r = rt dt.

T 0

Here, r is the average interest rate over the period and the function BS is

the Black-Scholes formula (5.2). In other words by replacing the interest rate

by its average over the period and the initial value of the stock by S0 e(c’r)T ,

the Black-Scholes formula provides the value for an option on an asset driven

by (5.17) conditional on the value of r. The constant interest rate model is a

useful control variate for the more general model (5.16). The expected value

E[BS(S0 e(c’r)T , K, r, T, σ)] can be determined by generating the interest rate

processes and averaging values of BS(S0 e(c’r)T , K, r, T, σ). Finally we may es-

timate the required option price under (5.15),(5.16) using an average of values

of

Z T

rt dt}[(ST ’ K)+ ’ (ZT ’ K)+ ]} + E{BS(S0 e(c’r)T , K, r, T, σ)}

exp{’

0

for ST and ZT generated using common random numbers.

We are still able to make a choice of the constant c. One simple choice is c ≈

E(r) since this means that the second term is approximately E{BS(S0 , K, r, T, σ)}.

Alternatively we can again experiment with small numbers of test simulations

and various values of c in an e¬ort to roughly minimize the variance

Z T

rt dt}[(ST ’ K)+ ’ (ZT ’ K)+ ]}).

var(exp{’

0

Evidently it is fairly easy to arrive at a solution in the case ρ = 0 since

we really only need to average values of the Black Scholes price under various

randomly generated interest rates. This does not work in the case ρ 6= 0 because

SIMULATING BARRIER AND LOOKBACK OPTIONS 269

the conditioning involved in (5.18) does not result in the Black Scholes formula.

Nevertheless we could still use common random numbers to generate two interest

rate paths, one corresponding to ρ = 0 and the other to ρ 6= 0 and use the

former as a control variate in the estimation of an option price in the general

case.

Importance Sampling

The expectation under the correct model could also be determined by multi-

plying this random variable by the ratio of the two likelihood functions and

then taking the expectation under E0 . By Girsanov™s Theorem, E{V (ST , T )} =

dP

E0 {V (ST , T ) dP0 } where P is the measure on the set of stock price paths corre-

sponding to (5.15),(5.16) and P0 that measure corresponding to (5.17). The

required Radon-Nykodym derivative is

ZT ZT 2

(rt ’ c2 )St

2

(rt ’ c)St

dP

dSt ’ (5.19)

= exp{ dt}

2 2

St σ 2 2σ 2 St

dP0 0 0

ZT ZT 2

rt ’ c2

rt ’ c

dSt ’ (5.20)

= exp{ dt}

St σ 2 2σ 2

0 0

The resulting estimator of the value of the option is therefore an average

over all simulations of the value of

Z Z Z

T T T

rt ’ c2

2

rt ’ c

dSt ’ (5.21)

V (ST , T )exp{’ rt dt + dt}

σ 2 rt 2σ2

0 0 0

where the trajectories rt are simulated under interest rate model (5.16).

As discussed before, we can attempt to choose the drift parameter c to

approximately minimize the variance of the estimator (5.21).

Simulating Barrier and lookback options

For a ¬nancial times series Xt observed over the interval 0 · t · T , what

is recorded in newspapers is often just the initial value or open of the time

270 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

series O = X0 , the terminal value or close C = XT , the maximum over the

period or the high, H = max{Xt ; 0 · t · T } and the minimum or the low

L = min{Xt ; 0 · t · T }. Very few uses of the highly informative variables H

and L are made, partly becuase their distribution is a bit more complicated than

that of the normal distribution of returns. Intuitively, however, the di¬erence

between H and L should carry a great deal of information about one of the

most important parameters of the series, its volatility. Estimators of volatility

obtained from the range of prices H ’ L or H/L will be discussed in Chapter 6.

In this section we look at how simple distributional properties of H and L can

be used to simulate the values of certain exotic path-dependent options.

Here we consider options such as barrier options, lookback options and hind-

sight options whose value function depends only on the four variables (O, H, L, C)

for a given process. Barrier options include knock-in and knock-out call options

and put options. Barrier options are simple call or put options with a fea-

ture that should the underlying cross a prescribed barrier, the option is either

knocked out (expires without value) or knocked in (becomes a simple call or

put option). Hindsight options, also called ¬xed strike lookback options are like

European call options in which we may use any price over the interval [0, T ]

rather than the closing price in the value function for the option. Of course for

a call option, this would imply using the high H and for a put the low L. A few

of these path-dependent options are listed below.

Option Payo¬

(C ’ K)+ I(H · m)

Knock-out Call

(C ’ K)+ I(H ≥ m)

Knock-in Call

H’C

Look-back Put

C’L

Look-Back Call

(H ’ K)+

Hindsight (¬xed strike lookback) Call

(K ’ L)+

Hindsight (¬xed strike lookback) put

Table 5.1: Value Function for some exotic options

SIMULATING BARRIER AND LOOKBACK OPTIONS 271

For further details, see Kou et. al. (1999) and the references therein.

Simulating the High and the Close

All of the options mentioned above are functions of two or three variables O, C,

and H or O, C, and L and so our ¬rst challenge is to obtain in a form suitable

to calculation or simulation the joint distribution of these three variables. Our

argument will be based on one of the simplest results in combinatorial proba-

bility, the re¬‚ection principle. We would like to be able to handle more than

just a Black-Scholes model, both discrete and continuous distributions, and we

begin with the simple discrete case. Much of the material in this section can be

found in McLeish(2002).

In the real world, the market does not rigorously observe our notions of the

passage of time. Volatility and volume traded vary over the day and by day

of the week. A successful model will permit some variation in clock speed and

volatility, and so we make an attempt to permit both in our discrete model.

In the discrete case, we will assume that the stock price St forms a trinomial

tree, taking values on a set set D = {· · · d’1 < d0 < d1 · · · }... At each time

point t, the stock may either increase, decrease, or stay in the same place and

the probability of these movements may depend on the time. Speci¬cally we

assume that if St = di , then for some parameters θ, pt , t = 1, 2, ...,

§

⎪ pt eθ if j = i + 1

⎪

⎪

⎪

⎪

⎪

⎨ 1 ’ 2pt if j = i

1

= dj |St = di ) = — (5.22)

P (St+1

kt (θ) ⎪ p e’θ

⎪t if j = i ’ 1

⎪

⎪

⎪

⎪

©

otherwise

0

1

for all t. If we choose all pt = 1 ,

where kt (θ) = 1 + pt (eθ + e’θ ’ 2) and pt · 2 2

then this is a model of a simple binomial tree which either steps up or down

at each time point. The increment in this process Xt+1 = St+1 ’ St has mean

272 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

which depends on the time t except in the case θ = 0

pt

{(di+1 ’ di )eθ ’ (di ’ di’1 )e’θ ),

E(Xt+1 |Xt = di ) =

kt (θ)

and variance, also time-dependent except in the case θ = 0. The parameter θ

describes one feature of this process which is not dependent on the time or the

location of the process, since the log odds of a move up versus a move down is

P [UP]

2θ = log[ ].

P [DOWN]

Suppose we label the states of the process so that S0 = d0 and there is a barrier

at the point dm where m > 0. The main result conerning the distribution of

the high (or low) is the following:

Proposition 43 Suppose a stock price St has dynamics determined by (5.22),

and S0 = d0 . De¬ne

H = max St and C = ST

0tT

Then for u < m,

P0 [C = d2m’u ]

Pθ (H ≥ dm |C=du )= , for min(u, 0) < m, (5.23)

P0 [C = du ]

= 1, for min(u,0) ≥ m

Proof. We wish to count the number of paths over an interval of time [0, T ]

which touch or cross this barrier and end at a particular point du , u < m. Such

a path is shown as a solid line in Figure 5.1 in the case that the points di are

all equally spaced. Such a path has a natural “re¬‚ection” about the horizontal

line at dm . The re¬‚ected path is identical up to the ¬rst time „ that the original

path touches the point dm , and after this time, say at time t > „, the relected

path takes the value d2m’i where St = di . This path is the dotted line in Figure

5.1. Notice that if the original path ends at du < dm , below the barrier, the

re¬‚ected path ends at d2m’u > dm or above the barrier. Each path touching

SIMULATING BARRIER AND LOOKBACK OPTIONS 273

15

2m-u

10

5

m

0

-5

u

-10

0 10 20 30 40 50 60 70 80 90 100

Figure 5.1: Illustration of the Re¬‚ection Principle

the barrier at least once and ending below it at du has a re¬‚ected path ending

above it at d2m’u , and of course each path that ends above the barrier must

touch the barrier for a ¬rst time at some point and has a re¬‚ection that ends

below the barrier. This establishes a one-one correspondence useful for counting

these paths. Let us denote by the symbol “#” the “number of paths such that”.

Then:

#{H ≥ dm and C = du < dm } = #{C = d2m’u }.

Now consider the probability of any path ending at a particular point du ,

(S0 = d0 , S1 , ..., ST = du ).

To establish this probability, each time the process jumps up in this interval

pt eθ

we must multiply by the factor and each time there is a jump down we

kt (θ)

pt e’θ 1’2pt

multiply by If the process stays in the same place we multiply by

kt (θ) . kt (θ) .

The re¬‚ected path has exactly the same factors except that after the time „ at

which the barrier is touched, the “up” jumps are replaced by “down” jumps and

vice versa. For an up jump in the original path multiply by e’2θ . For a down

jump in the original path, multiply by e2θ . Of course this allows us to compare

path probabilities for an arbitrary value of the parameter θ, say with P0 , the

274 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

probability under θ = 0 since, if the path ends at C = du ,

eNU θ e’ND θ

Q

Pθ (path) = P0 [path]

kt (θ)

t

euθ

=Q (5.24)

P0 [path]

kt (θ)

t

where NU and ND are the number of up jumps and down jumps in the path.

Note that we have subscripted the probability measure by the assumed value of

the parameter θ. This makes it easy to compare the probabilities of the original

and the re¬‚ected path, since

Pθ [original path]

= e’2θNU e2θND

Pθ [re¬‚ected path]

where now the number of up and down jumps NU and ND are counted following

time „. However, since ST = du and S„ = dm , it follows that ND ’ NU = m ’ u

and that

P [original path]

= e2θ(u’m)

P [re¬‚ected path]

which is completely independent of how that path arrived at the closing value

du , depending only on the close. This makes it easy to establish the probability

of paths having the property that H ≥ dm and C = du < dm since there are

exactly the same number of paths such that C = d2m’u and the probabilities of

these paths di¬er by a constant factor e2θ(u’m) . Finally this provides the useful

result for u < m.

Pθ [H ≥ dm and C = du ] = e2θ(u’m) Pθ [C = d2m’u ],

or, on division by P [C = du ],

e2θ(u’m) Pθ [C = d2m’u ]

Pθ [H ≥ dm |C = du ] =

Pθ [C = du ]

e2θ(u’m) eθ(2m’u) P0 [C = d2m’u ]

=

eθu P0 [C = du ]

P0 [C = d2m’u ]

=

P0 [C = du ]

SIMULATING BARRIER AND LOOKBACK OPTIONS 275

where we have used (5.24). This rather simple formula completely descibes the

conditional distribution of the high under an arbitrary value of the parameter

θ in terms of the value of the close under parameter value θ = 0.

Thus, the conditional distribution of the high of a process given the open and

close can be determined easily without knowledge of the underlying parameter

and is related to the distribution of the close when the “drift” θ = 0. This result

also gives the expected value of the high in fairly simple form if the points dj are

equally spaced. Suppose dj = j∆ for j = 0, ±1, ±2, .... Then for u = j∆,with

j ≥ 0 and k ≥ 1, (see Problem 1)

Pθ [H ’ C ≥ k∆|C = j∆] =

P [C > u and C’u is even]

∆

E[H|C = u] = u + ∆ .

P [C = u]

Roughly, (5.23) indicates that if you can simulate the close under θ, then

you use the properties of the close under θ = 0 to simulate the high of the

process. Consider the problem of simulating the high for a given value of the

close C = ST = du and again assuming that S0 = d0 . Suppose we use inverse

transform from a uniform random variable U to solve the inequalities

Pθ ( max St ≥ dm+1 |ST = du ) < U · Pθ ( max St ≥ dm |ST = du )

0tT 0tT

for the value of dm . In this case the value of

dm = sup{dj ; U P0 [ST = du ] · P0 [ST = d2j’u ]}

is the generated value of the high. This inequality is equivalent to

P0 [ST = d2m+2’u ] < U P0 [ST = du ] · P0 [ST = d2m’u ].

Graphically this inequality is demonstrated in Figure ?? which shows the prob-

ability histogram of the distribution ST under θ = 0. The value U P0 [ST = du ]

is the y-coordinate of a point P

276 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

0.14

0.12

0.1

0.08

0.06

0.04

0.02

P

0

d d

d u 2m-u

0

Generating a High for a discrete distribution

randomly chosen from the bar corresponding to the point du . The high dm is

generated by moving horizontally to the right an even number of steps until just

before exiting the histogram. This is above the value d2m’u and dm is between

du and d2m’u .

A similar result is available for Brownian motion and Geometric Brownian

motion. A justi¬cation of these results can be made by taking a limit in the

discrete case as the time steps and the distances dj ’ dj’1 all approach zero. If

we do this, the parameter θ is analogous to the drift of the Brownian motion.

The result for Brownian motion is as follows:

Theorem 44 Suppose St is a Brownian motion process

dSt = µdt + σdWt ,

S0 = 0, ST = C,

H = max{St ; 0 · t · T } and

L = min{St ; 0 · t · T }.

If f0 denotes the Normal(0,σ 2 T ) probability density function, the distribution of

SIMULATING BARRIER AND LOOKBACK OPTIONS 277

C under drift µ = 0, then

f0 (2H ’ C)

is distributed as U [0, 1] independently of C,

UH =

f0 (C)

f0 (2L ’ C)

is distributed as U [0, 1] independently of C.

UL =

f0 (C)

1

ZH = H(H ’ C) is distributed as Exponential ( σ 2 T ) independently of C,

2

1

ZL = L(L ’ C) is distributed as Exponential ( σ 2 T )independently of C.

2

We will not prove this result since it is a special case of Theorem 46 below.

However it is a natural extension of Proposition 43 in the special case that

dj = j∆ for some ∆ and so we will provide a simple sketch of a proof using

this proposition. Consider the ratio

P0 [C = d2m’u ]

P0 [C = du ]

on the right side of (5.23). Suppose we take the limit of this as ∆ ’ 0 and as

m∆ ’ h and u∆ ’ c. Then this ratio approaches

f0 (2h ’ c)

f0 (c)

where f0 is the probability density function of C under µ = 0. This implies for

a Brownian motion process,

f0 (2h ’ c)

P [H ≥ h|C = c] = for h ≥ c. (5.25)

f0 (c)

If we temporarily denote the cumulative distribution function of H given C = c

by Gc (h) then (5.25) gives an expression for 1 ’ Gc (h) and recall that since

the sumulative distribution function is continuous, when we evaluate it at the

observed value of a random variable we obtain a U [0, 1] random variable e.g.

Gc (H) ∼ U [0, 1]. In other words conditional on C = c we have

f0 (2H ’ c)

∼ U [0, 1].

f0 (c)

This result veri¬es a simple geometric procedure, directly analogous to that in

Figure 5.2, for generating H for a given value of C = c. Suppose we gener-

ate a point PH = (c, y) under the graph of f0 (x) and uniformly distributed

278 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

Figure 5.2: Generating H for a ¬xed value of C for a Brownian motion.

on {(c, y); 0 · y · f0 (c)}. This point is shown in Figure ??. We regard the

y’coordinate of this point as the generated value of f0 (2H ’ c). Then H can be

found by moving from PH horizontally to the right until we strike the graph of

f0 and then moving vertically down to the axis (this is now the point 2H ’ c)

and ¬nally taking the midpoint between this coordinate 2H ’ c and the close

c to obtain the generated value of the high H. The low of the process can be

generated in the same way but with a di¬erent point PL uniform on the set

{(c, y); 0 · y · f0 (c)}. The algorithm is the same in this case except that we

move horizontally to the left.

There is a similar argument for generating the high under a geometric Brown-

ian motion as well, since the logarithm of a geometric Brownian motion is a

Brownian motion process.

Corollary 45 For a Geometric Brownian motion process

dSt = µSt dt + σSt dWt ,

S0 = O and ST = C

SIMULATING BARRIER AND LOOKBACK OPTIONS 279

with f0 the normal(0, σ2 T ) probability density function, we have

1

ln(H/O) ln(H/C) ∼ exp( σ 2 T ) independently of O, C and

2

1

ln(L/O) ln(L/C) ∼ exp( σ 2 T ) independently of O, C.

2

f0 (ln(H 2 /OC))

∼ U [0, 1] independently of O, C and

UH =

f0 (ln(C/O))

f0 (ln(L2 /OC))

∼ U [0, 1] independently of O, C.

UL =

f0 (ln(C/O))

Both of these results are special cases of the following more general Theorem.

We refer to McLeish(2002) for the proof. As usual, we consider a price process

St and de¬ne the high H = max{St ; 0 · t · T }, and the open and close O = S0 ,

C = ST .

Theorem 46 Suppose the process St satis¬es the stochastic di¬erential equa-

tion:

1

dSt = {ν + σ 0 (St )}σ(St )»2 (t)dt + σ(St )»(t)dWt (5.26)

2

where σ(x) > 0 and »(t) are positive real-valued functions such that g(x) =

Rx 1 RT

dy and θ = 0 »2 (s)ds < ∞ are well de¬ned on <+ .

σ(y)

(a) Then with f0 the N (0, θ) probability density function we have

f0 {2g(H) ’ g(O) ’ g(C)}

∼ U [0, 1]

UH =

f0 {g(C) ’ g(O)}

and UH is independent of C.

(b) For each value of T , ZH = (g(H) ’ g(O))(g(H) ’ g(C)) is independent of

O, C, and has an exponential distribution with mean 1 θ.

2

A similar result holds for the low of the process over the interval, namely

that

f0 {2g(L) ’ g(O) ’ g(C)}

∼ U [0, 1]

UL =

f0 {g(C) ’ g(O)}

280 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

and ZH = {g(L) ’ g(O)}{g(L) ’ g(C)} is independent of O, C, and has an

exponential distribution with mean 1 θ.

2

Before we discuss the valuation of various options, we examing the signif-

icance of the ratio appearing in on the right hand side of (5.25) a little more

closely. Recall that f0 is the N (0, σ 2 T ) probability density function and so we

can replace it by

2

exp{’ (2h’c) }

f0 (2h ’ c) zh

2σ 2 T

= exp{’2 2 } (5.27)

= 2

c

f0 (c) σT

exp{’ 2σ2 T }

where zh = h(h ’ c) or in the more general case where S(0) = O may not be

equal to zero,

zh = (h ’ O)(h ’ c). (5.28)

f0 (2h’c)

This ratio represents the probability that a particular process with

f0 (c)

close c breaches a barrier at h and so the exponent

zh

2

σ2T

in the right hand side of (5.27) controls the probability of this event.

Of course we can use the above geometric algorithm for Brownian motion to

generated highs and closing prices for a geometric Brownian motion, for exam-

ple, St satisfying d ln(St ) = σdWt (minor adjustments required to accommodate

nonzero drift). The graph of the normal probability density function f0 (x) of

ln(C) is shown in Figure ??.

If a point PH is selected at random uniformly distributed in the region below

the graph of this density, then, by the usual arguments supporting the accep-

tance rejection method of simulation, the x-coordinate of this point is a variate

generated from the probability density function f0 (x), that is, a simulated value

from the distribution of ln(C). The y-coordinate of such a randomly selected

point also generates the value of the high as before.If we extend a line horizon-

tally to the right from PH until it strikes the graph of the probability density and

SIMULATING BARRIER AND LOOKBACK OPTIONS 281

Figure 5.3: Simulating from the joint distribution of (H, C) or from (L, C)

then consider the abscissa, of this point, this is the simulated value of ln(H 2 /C),

and ln(H) the average the simulated values of ln(H 2 /C) and ln(C).

A similar point PL uniform under the probability density function f0 can be

used to generate the low of the process if we extend the line from PL to the left

until it strikes the density. Again the abscissa of this point is ln(L2 /C) and the

average with ln(C) gives a simulated value of ln(L). Although the y’coordinate

of both PH and PL are uniformly distributed on [0, f0 (C)] conditional on the

value of C they are not independent.

Suppose now we wish to price a barrier option whose payo¬ on maturity

depends on the value of the close C but provided that the high H did not

exceed a certain value, the barrier. This is an example of an knock-out barrier

but other types are similarly handled. Once again we assume the simplest form

of the geometric Brownian motion d ln(St ) = σdWt and assume that the upper

barrier is at the point Oeb so that the payo¬ from the option on maturity T is

ψ(C)I(H < Oeb )

for some function ψ. It is clear that the corresponding value of H does not

exceed a boundary at Oeb if and only if the point PH is below the graph of the

282 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

Figure 5.4: Simulating a knock-out barrier option with barrier at Oeb

probability density function but not in the shaded region obtained by re¬‚ecting

the right hand tail of the density about the vertical line x = b ’ ln(O) in

Figure 5.4. To simulate the value of the option, choose points uniformly under

the graph of the probability density f0 (x). For those points in the non-shaded

region under f0 (the x-coordinate of these points are simulated values ψ(C)of

ln(C) under the condition that the barrier is not breached) we average the values

of ψ(C) and for those in the shaded region we average 0.

Equivalently,

Eψ(C)I(H < Oeb ) = Eψ — (C)

where

§

⎨ ψ(C) C · Oeb

for

—

ψ (C) = .

© ’ψ(2b + ln(O2 /C)) b

for C > Oe

and so the barrier option can be priced as if it were a vanilla European option

with payo¬ function ψ — (C).

Any option whose value depends on the high and the close of the process

(or (L, C)) can be similarly valued as a European option. If an option be-

comes worthless whenever an upper boundary at Oeb is breached, we need only

SIMULATING BARRIER AND LOOKBACK OPTIONS 283

multiply the payo¬ from the option ignoring the boundary by the factor

zh

1 ’ exp{’2 }

σ2 T

with

zh = b(b + ln(O/C))

to accommodate the ¬ltering e¬ect of the barrier and then value the option as

if it were a European option.

There is a a variety of distributional results related to H, some used by

Redekop (1995) to test the local Brownian nature of various ¬nancial time series.

These are easily seen in Figure 5.5. For example, for a Brownian motion process

with sero drift, suppose we condition on the value of 2H ’ O ’ C. Then the

point PH must lie (uniformly distributed) on the line L1 and therefore the point

H lies uniformly on this same line but to the right of the point O. This shows

that conditional on 2H ’ O ’ C the random variable H ’ O is uniform or,

H’O

∼ U [0, 1].

2H ’ O ’ C

Similarly, conditional on the value of H, the point PH must fall somewhere on

the curve labelled C2 whose y-coordinate is uniformly distributed showing that

C’O

∼ U [0, 1].

2H ’ O ’ C

Redekop shows that for a Brownian motion process, the statistic

H’O

(5.29)

2H ’ O ’ C

is supposed to be uniformly [0, 1] distributed but when evaluated using real

¬nancial data, is far too often close to or equal the extreme values 0 or 1.

The joint distribution of (C, H) can also be seen from Figure 5.6. Note

that the rectangle around the point (x, y) of area ∆x∆y under the graph of the

density, when mapped into values of the high results in an interval of values for

284 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

Figure 5.5: Some uniformly distributed statistics for Brownian Motion

(2H ’ C) of width ’∆y/φ0 (2y ’ x) where φ0 is the derivative of the standard

normal probability density function (the minus sign is to adjust for the negative

slope of the density here). This interval is labelled ∆(2H ’ C). This, in turn

generates the interval ∆H of possible values of H, of width exactly half this, or

’∆y

.

2φ0 (2y ’ x)

Inverting this relationship between (x, y) and (H, C),

P [H ∈ ∆H, C ∈ ∆C] = ’2φ0 (2y ’ x)∆x∆y

con¬rming that the joint density of (H, C) is given by ’2φ0 (2y ’ x) for x < y.

In order to get the joint density of the High and the Close when the drift is

non-zero, we need only multiply by the ratio of the two normal density functions

of the close

fµ (x)

f0 (x)

and this gives the more general result in the table below.

The table below summarizes many of our distributional results for a Brown-

ian motion process with drift on the interval [0, 1],

dSt = µdt + σdWt , with S0 = O.

SIMULATING BARRIER AND LOOKBACK OPTIONS 285

Figure 5.6: Con¬rmation of the joint density of (H, C)

Statistic Density Conditions

’∞ < x < y,

X = C ’ O,

f (y, x) = ’2φ0 (2y ’ x) exp(µx ’ µ2 /2) and y > 0, σ = 1

Y =H ’O

given O

fY |X (y|x) = 2(2y ’ x)e’2y(y’x)

Y |X y > x, σ = 1

exp((σ 2 /2)

Z = Y (Y ’ X) given O, X

exp(σ2 /2)

(L ’ O)(L ’ C) given (O, C)

exp(σ2 /2)

(H ’ O)(H ’ C) given (O, C)

H’O

drift ν = 0, given O, 2H ’ O ’ C

U [0, 1]

2H’O’C

L’O

drift ν = 0, given O, 2L ’ O ’ C

U [0, 1]

2L’O’C

C’O

drift ν = 0, given H, O

U [’1, 1]

2H’O’C

TABLE 5.2: Some distributional results for High, Close and Low.

We now consider brie¬‚y the case of non-zero drift for a geometric Brownian

motion. Fortunately, all that needs to be changed in the results above is the

marginal distribution of ln(C) since all conditional distributions given the value

of C are the same as in the zero-drift case. Suppose an option has payo¬ on

286 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

maturity ψ(C) if an upper barrier at level Oeb , b > 0 is not breached. We have

already seen that to accommodate the ¬letering e¬ect of this knock-out barrier

we should determine, numerically or by simulation, the expected value

b(b + ln(O/C))

E[ψ(C)(1 ’ exp{’2 })]

σ2 T