. 1
( 2)



>>

Chapter 5


Simulating the value of
Options

Asian Options

An Asian option, at expiration T, has value determined not by the closing price
of the underlying asset as for a European option, but on an average price of
the asset over an interval. For example a discretely sampled Asian call op-
tion on an asset with price process S(t) pays an amount on maturity equal
Pk
¯ ¯ 1
to max(0, Sk ’ K) where Sk = k i=1 S(iT /k) is the average asset price at k
equally spaced time points in the time interval (0, T ). Here, k depends on the
frequency of sampling (e.g. if T = .25 (years) and sampling is weekly, then
¯
k = 13. If S(t) follows a geometric Brownian motion, then Sk is the sum of
lognormally distributed random variables and the distribution of the sum or av-
erage of lognormal random variables is very di¬cult to express analytically. For
this reason we will resort to pricing the Asian option using simulation. Notice,
however that in contrast to the arithmetic average, the distribution of the geo-
metric average has a distribution which can easily be obtained. The geometric

265
266 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

Pn
1
mean of n values X1 , ..., Xn is (X1 X2 ...Xn )1/n = exp{ n ln(Xi )} and if the
i=1

random variables Xn were each lognormally distributed then this results adding
the normally distributed random variables ln(Xi ) in the exponent, a much more
P
familiar operation. In fact the sum in the exponent n n ln(Xi ) is normally
1
i=1

distributed so the geometric average will have a lognormal distribution.

Our objective is to determine the value of the Asian option E(V1 ) with




¯
V1 = e’rT max(0, Sk ’ K)


Since we expect geometric means to be close to arithmetic means, a reasonable
˜ ˜
control variate is the random variable V2 = e’rT max(0, Sk ’ K) where Sk =
Qk
{ i=1 S(iT /k)}1/k is the geometric mean. Assume that V1 and V2 obtain from
the same simulation and are therefore possibly correlated. Of course V2 is only
useful as a control variate if its expected value can be determined analytically
or numerically more easily than V1 but in view of the fact that V2 has a known
lognormal distribution, the prospects of this are excellent. Since S(t) = S0 eY (t)
where Y (t) is a Brownian motion with Y (0) = 0, drift r ’ σ 2 /2 and di¬usion σ,
˜
it follows that Sk has the same distribution as does

k
1X
(5.1)
S0 exp{ Y (iT /k)}.
k i=1


This is a weighted average of the independent normal increments of the process
and therefore normally distributed. In particular if we set

k
1X
¯
Y= Y (iT /k)
k i=1
1
[k(Y (T /k)) + (k ’ 1){Y (2T /k) ’ Y (T /k)} + (k ’ 2){Y (3T /k) ’ Y (2T /k)}
=
k
+ ... + {Y (T ) ’ Y ((k ’ 1)T /k)}],
ASIAN OPTIONS 267

then
k
X
2
¯ ) = r ’ σ /2
E(Y iT /k
k i=1
σ2 k + 1
= (r ’ ) T
2 2k
e
= µT, say,

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)
268 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
e
the basic Black-Scholes because in this case, the Asian option corresponds to
the average of a single observation. For k > 1 the e¬ective initial stock price
f
is reduced S0 < S0 and the volatility parameter is also smaller σ 2 < σ 2 . With
e
lower initial stock price and smaller volatility the price of a European call will
decrease, indicating that if 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 269

σ = 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);
270 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/1.1,.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 271

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 say 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 so that 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 remedy 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 )
272 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 that the strike price is attained. The
importance sampling adjustment that insures that the estimator continues to be
an unbised 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
, σ2 T )
•(ZT ; rT ’
’rT ZT + 2
’ K)
= EQ0 (e (S0 e )
σ2 T
, σ2T )
•(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 273

σ = 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)
results in an e¬ciency gain of around 73, in part because very few of the crude
estimates 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 the price of 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 to have a payo¬ related to the closing value of the stock S(T ). It might
274 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

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, or 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 .

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 .

Assume that in both cases, the process starts at the same initial value S0 .Then
dP
the “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 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 275

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
276 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 with 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 both will depend on an overlaping set of factors. However
PRICING A CALL OPTION UNDER STOCHASTIC INTEREST RATES.277

we begin 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
278 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 279

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 multiply-
ing 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
280 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 XX: Value Function for some exotic options
For further details, see Kou et. al. (1999) and the references therein.
SIMULATING BARRIER AND LOOKBACK OPTIONS 281

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.
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
P (St+1 = dj |St = di ) = — (5.22)
kt (θ) ⎪ p e’θ if j = i ’ 1
⎪t




©
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
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
282 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




Figure 5.1: Illustration of the Re¬‚ection Principle



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

P [UP]
2θ = log[ ].
P [DOWN]




Now suppose we label the states of the process so that S0 = d0 and there is
a barrier at the point dm where m > 0. 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 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
SIMULATING BARRIER AND LOOKBACK OPTIONS 283

“#” 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
probability under θ = 0 since, if the path ends at C = du ,

eNU θ e’ND θ
Q
Pθ (path) = P0 [path]
kt (θ)
t
euθ
=Q (5.23)
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]
284 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS

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 ]

where we have used (5.23). 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. In fact we have
proved the following proposition.



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.24)
P0 [C = du ]
= 1, for min(u,0) ≥ m

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
SIMULATING BARRIER AND LOOKBACK OPTIONS 285

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.24) 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 5.2 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 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
286 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




Figure 5.2: Generating a High, discrete distributions



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
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
SIMULATING BARRIER AND LOOKBACK OPTIONS 287

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.24). 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
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
288 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




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



{(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

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))
SIMULATING BARRIER AND LOOKBACK OPTIONS 289

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)}
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 }
290 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




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



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
SIMULATING BARRIER AND LOOKBACK OPTIONS 291

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
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
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.5. 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)
292 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




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



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
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.7. For example, for a Brownian motion process
SIMULATING BARRIER AND LOOKBACK OPTIONS 293

0.4




0.35




0.3





0.25
P
probability




H

0.2




0.15




0.1




m
0.05
ln(O)
0
-4 -3 -2 -1 0 1 2 3 4

ln(C)


Figure 5.6: Simulating a Barrier Option with barrier at em



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.
294 CHAPTER 5. SIMULATING THE VALUE OF OPTIONS




Figure 5.7: Some uniformly distributed statistics for Brownian Motion



The joint distribution of (C, H) can also be seen from Figure 5.8. 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
(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

. 1
( 2)



>>