0

1992 1993 1994 1995 1996 1997

X

3 month interest rates

3

2

Y*E-2

1

0

1992 1993 1994 1995 1996 1997

X

Figure 15.4. Interest rates time series: German (thick dotted line),

Japanese (thin dotted line), American (thick straight line), Thai (thin

straight line). XFGbasket.xpl

self-¬nancing strategy which produces a positive expected payo¬:

• at time t

“ borrow the portfolio (1+rj )’1 E(θj,t+h |Ft )Xj,t from the countries

composing the basket,

“ lend (1 + r0 )’1 Yt to Thailand,

“ invest the di¬erence (1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt in

the numeraire currency at the risk-free rate r1 ,

• at time t + h

“ withdraw the amount Yt+h from Thailand,

15.2 Estimating the coe¬cients of an exchange rate basket 335

E(θj,t+h |Ft )Xj,t+h ,

“ pay back the loan of

“ keep the di¬erence.

The expression for the pro¬t and for its expected value are:

p

= Yt+h ’ E(θj,t+h |Ft )Xj,t+h

Πt+h

j=1

«

p

(1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt

+(1 + r1 )

j=1

«

p

(1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt .

E(Πt+h |Ft ) = (1 + r1 )

j=1

15.2.2 Estimation results

For the implementation of the investment strategy described above one needs

the estimate of the, possibly time-varying, basket weights. The precision of

the estimation has a direct impact on the economic result of the investment.

Therefore, we compare four di¬erent estimators of the basket weights: the

adaptive, the recursive, the window and the Kalman ¬lter estimator using

economic criteria for a one month and for a three month investment horizon.

In particular we compute the average expected pro¬t and the average realized

pro¬t.

The adaptive estimation procedure requires three parameters: m, » and µ. The

choice of m0 does not in¬‚uence the results very much and it can be reasonably

set to 30. This value represents the minimal amount of data which are used

for the estimation, and in the case of a structural break, the minimal delay

before having the chance of detecting the change point. The selection of » and

µ is more critical. These two values determine the sensitivity of the algorithm.

Small values would imply a fast reaction to changes in the regressor coe¬cients,

but but they would also lead to the selection of intervals of homogeneity which

are possibly too small. Large values would imply a slower reaction and con-

sequently the selection of intervals which can be too large. To overcome this

problem we suggest the following approach.

The main idea is that small changes in the values of » and µ should not a¬ect

the estimation results. Therefore we restrict our attention on a set S of possible

336 15 Locally time homogeneous time series modeling

USD weight

8

0.028+Y*E-2

6

4

2

1994 1995 1996 1997

X

DEM weight

6

4

Y*E-2

2

0

1994 1995 1996 1997

X

JPY weight

0.7

0.6

0.5

Y

0.4

0.3

0.2

1994 1995 1996 1997

X

Figure 15.5. Estimated exchange rate basket weights: adaptive

(straight line), recursive (thine dotted line), window (thick dotted line).

pairs (», µ). In the present context we chose all the even number between 2

and 8:

S = {(», µ)| », µ ∈ {2, 4, 6, 8}}

15.2 Estimating the coe¬cients of an exchange rate basket 337

Then we compare the 16 pairs with the following criterion at each time t:

2

«

t’1 d

(»— , µ— ) = arg min Ys ’ θj,s|s’h Xj,s .

(»,µ)∈S s=t’200 j=1

Finally, we estimate the value of θt+h|t with the selected pair (»— , µ— ). The

appeal of the above selection criterion consists of the fact that it leads to the

choice of the pair (», µ) which has provided the least quadratic hedging costs

over the past trading periods. Notice that in general we have di¬erent results

depending on the length of the forecasting horizon: here one and three month.

Figure 15.5 shows the results for the three month horizon. It is interesting to

see that the adaptive estimate tends to coincide with the recursive estimate

during the ¬rst half of the sample, more or less, while during the second half

of the sample it tends to follow the rolling estimate.

We remark that the problem of selecting free parameters is not speci¬c to the

adaptive estimator. The window estimator requires the choice of the length

of the window: k, while the Kalman ¬lter needs the speci¬cation of the data

generating process of θt and the determination of Σ and σ. In this application

k is set equal to 250, Σ and σ are estimated recursively from the data using the

OLS, while θ0|0 and P0|0 are initialized using the ¬rst 350 observations which

are then discarded. We remark that this choice is consistent with the one of

Christo¬ersen and Giorgianni (2000).

Table 15.2 shows the result of the simulated investment. The investments are

normalized such that at each trading day we take a short position of 100 USD

in the optimal portfolio of the hard currencies. The result refers to the period

April 9 1993 to February 12 1997 for the one month horizon investment and

June 7 1993 to February 12 1997 for the three month horizon investment. Notice

¬rst that the average realized pro¬ts are positive and, as far as the three month

investment horizon is concerned, they are signi¬cantly larger than zero among

all methods. This provides a clear evidence for the fact that arbitrage pro¬ts

were possible with in the framework of the Thai Bath basket for the period

under study. The comparison of the estimator also show the importance of

properly accounting for the time variability of the parameters. The recursive

estimator shows modest result as far as the realized pro¬ts are concerned and

the largest bias between expected the realized pro¬t. On one side, the bias is

reduced by the window estimator and by the Kalman ¬lter, but on the other

side these two methods provide a worse performance as far as the realized pro¬t

are concerned. Finally, the adaptive estimator appears to be the best one, its

338 15 Locally time homogeneous time series modeling

ONE MONTH HORIZON Recursive Window KF Adaptive

Average Expected Pro¬ts .772 .565 .505 .553

Average Realized Pro¬t .403 .401 .389 .420

Standard errors (.305) (.305) (.330) (.333)

THREE MONTH HORIZON Recursive Window KF Adaptive

Average Expected Pro¬ts 1.627 1.467 1.375 1.455

Average Realized Pro¬t 1.166 1.141 1.147 1.182

Standard errors (.464) (.513) (.475) (.438)

Table 15.2. Summary statistics of the pro¬ts.

bias is much smaller than the one of the recursive estimator and it delivers the

largest realized pro¬ts for both investment horizons.

15.3 Estimating the volatility of ¬nancial time

series

The locally time homogeneous approach appears to be also appropriate for the

estimation of the volatility of ¬nancial time series. In order to provide some

motivation we ¬rst describe the stylized facts of ¬nancial time series. Let St

de¬ne the price process of a ¬nancial asset such as stocks or exchange rates,

then the returns are de¬ned as follows:

Rt = ln St ’ ln St’1 .

Stylized facts of ¬nancial asset returns are: a leptokurtic density, variance

clustering and highly persistent autocorrelation of square and absolute returns

(see Figure 15.6). Further details and examples on this topic can be found in

Taylor (1986) and in Franke, H¨rdle and Hafner (2001).

a

15.3 Estimating the volatility of ¬nancial time series 339

returns

5

0

Y*E-2

-5

0 50 100 150 200 250 300

X

ACF of the absolute returns

15

10

Y*E-2

5

0

0 10 20 30 40 50

X

Figure 15.6. JPY/USD returns XFGretacf.xpl

15.3.1 The standard approach

The returns of ¬nancial time series are usually modeled by the following equa-

tion:

Rt = σ t µt

Where σt is a strictly positive process, which describes the dynamics of the vari-

ance of Rt , and ξt has a standard normal distribution: ξt ∼ N (0, 1). Standard

parametric models of the volatility are of (G)ARCH type:

2 2 2

σt = ω + ±Rt’1 + βσt’1 ,

340 15 Locally time homogeneous time series modeling

like in Engle (1995) and Bollerslev (1995), and of stochastic volatility type:

2 2

ln σt = θ0 + θ1 ln σt’1 + νt ,

as described by Harvey, Ruiz and Shephard (1995). These models have been

expanded in order to incorporate other characteristics of the ¬nancial return

time series: TARCH, EGARCH and QARCH explicitly assume an asymmet-

ric reaction of the volatility process to the sign of the observed returns, while

IGARCH and FIGARCH model the long memory structure of the autocorre-

lations of the square returns.

15.3.2 The locally time homogeneous approach

A common feature to all the models which have been cited in the previous

section is that they completely describe the volatility process by a ¬nite set of

parameters. The availability of very large samples of ¬nancial data has given

the possibility of constructing models which display quite complicated param-

eterizations in order to explain all the observed stylized facts. Obviously those

models rely on the assumption that the parametric structure of the process

remains constant through the whole sample. This is a nontrivial and possi-

bly dangerous assumption in particular as far as forecasting is concerned as

pointed out in Clements and Hendry (1998). Furthermore checking for param-

eter instability becomes quite di¬cult if the model is nonlinear, and/or the

number of parameters is large. Whereby those characteristics of the returns

which are often explained by the long memory and (fractal) integrated nature

of the volatility process, could also depend on the parameters being time vary-

ing. We want to suggest an alternative approach which relies on a locally time

homogeneous parameterization, i.e. we assume that the volatility σ follows a

jump process and is constant over some unknown interval of time homogeneity.

The adaptive algorithm, which has been presented in the previous sections,

also applies in this case; its aim consists in the data-driven estimation of the

interval of time homogeneity, after which the estimate of the volatility can be

simply obtained by local averaging.

15.3.3 Modeling volatility via power transformation

Let St be an observed asset process in discrete time, t = 1, 2, . . . , „ and Rt

are the corresponding returns: Rt = log(St /St’1 ) . We model this process via

15.3 Estimating the volatility of ¬nancial time series 341

the conditional heteroscedasticity assumption

Rt = σt µt , (15.10)

where µt , t ≥ 1 , is a sequence of independent standard Gaussian random

variables and σt is the volatility process which is in general a predictable

random process, that is, σt is measurable with respect to Ft’1 with Ft’1 =

σ(R1 , . . . , Rt’1 ) .

The model equation (15.10) links the volatility σt with the observations Rt

via the multiplicative errors µt . In order to apply the theory presented in

Section 15.1 we need a regression like model with additive errors. For this

reason we consider the power transformation, which leads to a regression with

additive noise and so that the noise is close to a Gaussian one, see Carroll and

Ruppert (1988). Due to (15.10) the random variable Rt is conditionally on

Ft’1 Gaussian and it holds

2 2

E Rt |Ft’1 = σt .

Similarly, for every γ > 0 ,

γ γ γ

E Rt |Ft’1 = σt E (|ξ|γ |Ft’1 ) = Cγ σt ,

2

γ 2

γ 2γ 2γ 2

= σt E (|ξ|γ ’ Cγ ) = σt Dγ

’ Cγ σt |Ft’1

E Rt

where ξ denotes a standard Gaussian r.v., Cγ = E|ξ|γ and Dγ = Var|ξ|γ .

2

Therefore, the process |Rt |γ allows for the representation

γ γ

|Rt |γ = Cγ σt + Dγ σt ζt , (15.11)

where ζt is equal (|ξ|γ ’ Cγ ) /Dγ . A suitable choice of the value of γ provides

that the distribution of

(|ξ|γ ’ Cγ ) /Dγ

is close to the normal. In particular the value of γ = 0.5 appears to be almost

optimal, see Figure 15.7.

15.3.4 Adaptive estimation under local time-homogeneity

The assumption of local time homogeneity means that the function σt is

constant within an interval I = [„ ’ m, „ ] , and the process Rt follows the

342 15 Locally time homogeneous time series modeling

Normal and powertransformed densities

1

Y

0.5

0

-1 0 1 2

X

Figure 15.7. Normal and power transformed densities for γ = 0.5.

XFGpowtrans.xpl

γ

regression-like equation (15.11) with the constant trend θI = Cγ σI which can

be estimated by averaging over this interval I :

1

|Rt |γ .

θI = (15.12)

|I|

t∈I

By (15.11)

Cγ Dγ 1 sγ

γ γ

θI = σt + σt ζt = θt + θt ζt (15.13)

|I| |I| |I| |I|

t∈I t∈I t∈I t∈I

with sγ = Dγ /Cγ so that

1

EθI = E θt , (15.14)

|I|

t∈I

2

s2 s2

γ γ 2

E θt ζt = E θt . (15.15)

|I|2 |I|2

t∈I t∈I

De¬ne also

s2

γ

2 2

vI =2 θt .

|I|

t∈I

15.3 Estimating the volatility of ¬nancial time series 343

In view of (15.15) this value is called the conditional variance of θI . Under

local homogeneity it holds θt is constantly equal to θI for t ∈ I , and hence,

EθI = θI

s2 θI

2

γ

2

vI = Var θI = .

|I|

A probability bound analogous to the one in Section 15.1 holds also in this

2

case. Let the volatility coe¬cient σt satisfy the condition b ¤ σt ¤ bB with

some constants b > 0, B > 1 . Then there exists aγ > 0 such that it holds for

every » ≥ 0

»2

√

P |θI ’ θ„ | > »vI ¤ 4 e»(1 + log B) exp ’ . (15.16)

2aγ

The proof of the statement above and some related theoretical results can be

found in Mercurio and Spokoiny (2000).

For practical application one has to substitute the unknown conditional stan-

dard deviation with its estimate: vI = sγ θI |I|’1/2 . Under the assumption of

time homogeneity within an interval I = [„ ’ m, „ ] equation (15.16) allows

to bound |θI ’ θJ | by »vI + µvJ for any J ‚ I, provided that » and µ

are su¬ciently large. Therefore we can apply the same algorithm described in

Section 15.1 in order to estimate the largest interval of time homogeneity and

the related value of θ„ . Here, as in the previous section, we are faced with

the choice of three tuning parameters: m0 , », and µ. Simulation studies and

repeated trying on real data by Mercurio and Spokoiny (2000) have shown that

the choice of m0 is not particularly critical and it can be selected between 10

and 50 without a¬ecting the overall results of the procedure.

As described in Section 15.2.2, the choice of » and µ is more delicate. The

in¬‚uence of » and µ is similar to the one of the smoothing parameters in

the nonparametric regression. The likelihood of rejecting a time homogeneous

interval decreases with increasing » and/or µ. This is clear from equation

(15.6). Therefore if » and µ are too large this would make the algorithm too

conservative, increasing the bias of the estimator, while too small values of

» and µ would lead to a frequent rejection and to a high variability of the

estimate. Once again, a way of choosing the optimal values of » and µ can be

made through the minimization of the squared forecast error. One has to de¬ne

a ¬nite set S of the admissible pair of » and µ. Then for each pair belonging

344 15 Locally time homogeneous time series modeling

(»,µ)

to S one can compute the corresponding estimate: θt and then select the

optimal pair and the corresponding estimate by the following criterion:

T 2

(»,µ)

|Rt |γ ’ θt

(», µ) = min .

»,µ∈S

t=0

Figure 15.8 shows the result of the on-line estimation of the locally time homo-

geneous volatility model for the JPY/USD exchange rate. The bottom plot, in

particular, shows the estimated length of the interval of time homogeneity: m,

at each time point.

15.4 Technical appendix

In this section we give the precise conditions under which the bound (15.4)

holds. De¬ne:

WI = VI’1 ,

VI = σ ’2 Xt Xt

t∈I

furthermore let wij,I denote the√ elements of WI . For some positive constants

b > 0, B > 1, ρ < 1, r ≥ 1, » > 2 and for i = 1 . . . p consider the random set

were the following conditions are ful¬lled:

’1

±

b ¤ wjj,I ¤ bB

||VI µ||2 wjj,I ¤ r

sup

Ai,I =

{µ∈RK :||µ||=1}

|wji,I /wjj,I | ¤ ρ ∀ i = 1, . . . , p

Let (Y1 , X1 ) . . . (Y„ , X„ ) obey (15.1), where the regressors are possibly stochas-

tic, then it holds holds for the estimate θI :

√

P |θi,I ’ θi,„ | > » wii,I ; Ai,I

¤ 4e ln(4B)(1 + 2ρ r(d ’ 1)»)p’1 » exp(’»2 /2), i = 1, . . . , p.

A proof of this statement can be found in Liptser and Spokoiny (1999). For

a further generalization, where the hypothesis of local time homogeneity holds

only approximatively, see H¨rdle et al. (2000).

a

15.4 Technical appendix 345

YEN/DM returns

10

5

Y*E-2

0

-5

-10

0 5 10 15 20

X*E3

Volatility

10

0.0002+Y*E-4

8

6

4

2

0 5 10 15 20

X*E3

Interval of homogeneity

10

Y*E2

5

0 5 10 15 20

X*E3

Figure 15.8. From the top: returns, estimated locally time homoge-

neous volatility and estimated length of the interval of time homogene-

ity. XFGlochom.xpl

Bibliography

Bollerslev, T. (1995). Generalised autoregressive conditional heteroskedasticity,

in Engle (1995).

346 15 Locally time homogeneous time series modeling

Carroll, R. and Ruppert, D. (1988). Transformation and Weighting in Regres-

sion, Chapman and Hall, New York.

Christo¬ersen, P. and Giorgianni, L. (2000). Interest rate in currency basket:

Forecasting weights and measuring risk, Journal of Business and Eco-

nomic Statistics 18: 321“335.

Chui, C. and Chen, G. (1998). Kalman Filtering, Information Sciences, third

edn, Springer-Verlag, Berlin.

Clements, M. P. and Hendry, D. F. (1998). Forecastng Economic Time Series,

Cambridge University Press, Cambridge.

Cooley, T. F. and Prescott, E. C. (1973). An adaptive regression model, Inter-

national Economic Review 14: 364“371.

Eichengreen, B., Masson, P., Savastano, M. and Sharma, S. (1999). Transition

Strategies and Nominal Anchors on the Road to Greater Exchange Rate

Flexibility, number 213 in Essays in International Finance, Princeton Uni-

versity Press.

Elliot, R. J., Aggoun, L. and Moore, J. B. (1995). Hidden Markov Models,

Springer-Verlag, Berlin.

Engle, R. F. (ed.) (1995). ARCH, selected readings, Oxford University Press,

Oxford.

Franke, J., H¨rdle, W. and Hafner, C. (2001). Einf¨hrung in die Statistik der

a u

Finanzm¨rkte, Springer, Berlin.

a

H¨rdle, W., Herwartz, H. and Spokoiny, V. (2001). Time inhomogeous multiple

a

volatility modelling. Discussion Paper 7, Sonderforschungsbereich 373,

Humboldt-Universit¨t zu Berlin. To appear in Financial Econometrics.

a

H¨rdle, W., Spokoiny, V. and Teyssi`re, G. (2000). Adaptive estimation for

a e

a time inhomogeneouse stochastic volatility model. Discussion Paper 6,

Sonderforschungsbereich 373, Humboldt-Universit¨t zu Berlin.

a

Harvey, A., Ruiz, E. and Shephard, N. (1995). Multivariate stochastic variance

models, in Engle (1995).

Lepski, O. (1990). One problem of adaptive estimation in gaussian white noise,

Theory Probab. Appl. 35: 459“470.

15.4 Technical appendix 347

Lepski, O. and Spokoiny, V. (1997). Optimal pointwise adaptive methods in

nonparametric estimation, Annals of Statistics 25: 2512“2546.

Liptser, R. and Spokoiny, V. (1999). Deviation probability bound for mar-

tingales with applications to statistical estimation, Stat. & Prob. Letter

46: 347“357.

Mercurio, D. and Spokoiny, V. (2000). Statistical inference for time-

inhomogeneous volatility models. Discussion Paper 583, Weierstrass In-

stitute for Applied Analysis and Stochastic, Berlin.

Mercurio, D. and Torricelli, C. (2001). Estimation and arbitrage opportunities

for exchange rate baskets. Discussion Paper 37, Sonderforschungsbereich

373, Humboldt-Universit¨t zu Berlin.

a

Musiela, M. and Rutkowski, M. (1997). Martingale Methods in Financial Mod-

elling, number 36 in Application of Mathemathics. Stochastic Modelling

and Applied Probability, Springer, New York.

Spokoiny, V. (1998). Estimation of a function with discontinuities via lo-

cal polynomial ¬t with an adaptive window choice, Annals of Statistics

26: 1356“1378.

Taylor, S. J. (1986). Modelling Financial Time Series, Wiley, Chichester.

16 Simulation based Option Pricing

Jens L¨ssem and J¨rgen Schumacher

u u

16.1 Simulation techniques for option pricing

We introduce Monte Carlo techniques and Quasi Monte Carlo techniques for

option pricing. First, we give an idea how to use simulation techniques to

determine option prices, then - using the developed basic methods - we give

examples how to price more complex i.e. exotic options even on more than one

underlying. Finally we present a short guideline how to price exotic options

with the proposed techniques.

First, we take a look at a European put on one underlying stock, a pricing

problem which can be solved analytically e.g. by using the Black-Scholes for-

mula. We start with this problem not only because it has become a kind of

”standard problem” but also to have the possibility to compare the results of

our approximation with an analytical solution. At the same time we look at

the time-complexity of the used simulation technique. Next, we show how to

price path dependent options with Monte Carlo methods. Afterwards, we show

how to price a stock option on several underlyings. This implies that we have

to solve a multi-dimensional simulation problem.

16.1.1 Introduction to simulation techniques

The idea behind randomized algorithms is that a random sample from a pop-

ulation (of input variables) is representative for the whole population. As a

consequence, a randomized algorithm can be interpreted as a probability dis-

tribution on a set of deterministic algorithms.

We will see that there are three main advantages to randomized algorithms:

1. Performance: For many problems, it can be shown that randomized algo-

350 16 Simulation based Option Pricing

rithms run faster than the best known deterministic algorithm. 2. Simplicity:

Randomized algorithms are easier to describe and implement than comparable

deterministic algorithms. 3. Flexibility: Randomized algorithms can be easily

adapted.

In general one distinguishes two types of randomized algorithms. Las Vegas

algorithms are randomized algorithms that always give correct results with

only the variation from one run to another being its running time. Monte

Carlo algorithms are randomized algorithms that may produce an incorrect

solution for which one can bound the probability of occurrence. The quality of

the solution can be seen as a random variable.

Within this chapter, we focus on Monte Carlo algorithms calculating the value

of the following integral

f (x)dx (16.1)

[0,1]d

by evaluation of f (x) for independent uniform distributed random vectors

X1 , X2 , . . . , Xn , Xi ∈ [0, 1]d .

The arithmetic mean of the values f (Xi ) can be seen as a guess for the expected

value of the random variable f (Xi ) and therefore can be interpreted as an

approximation for the value of the integral. According to the strong law of

large numbers the estimator for the expected value (the arithmetic mean of

the random function values) is converging to the expected value (the value of

the integral) with an increasing sample size. The probability that the absolute

error of the approximation result exceeds a ¬xed positive value is limited and

decreases to zero with an increasing sample size if the variance of f is ¬nite.

16.1.2 Pricing path independent European options on one

underlying

For the case of a European option on one underlying we have to approximate

the following integral via Monte Carlo simulation:

∞

r(T ’t)

CT (ST )g(ST |St , r, σ, T ’ t)dST (16.2)

e E [CT (ST )|St ] =

0

CT {f (x, St , r, σ, T ’ t)}dx

= (16.3)

[0,1)

16.1 Simulation techniques for option pricing 351

Where

2

)(T ’t)))2

exp ’ (log ST ’(log St ’(r’0.5σ

2σ 2 (T ’t)

g(ST |St , r, σ, T ’ t) =

2πσ 2 (T ’ t)ST

is the risk neutral density function of the Black Scholes model with parameters:

ST : price of the underlying at maturity

St : price of the underlying at time t

r : risk free interest rate

σ : volatility of log returns of the underlying

T ’t : time to maturity

√

1

ST = f (x, St , r, σ, T ’ t) = St exp{(r ’ σ 2 )(T ’ t) + σ T ’ tF ’1 (x)}

2

transforms the uniform distributed values x in g(ST |St , r, σ, T ’ t) distributed

underlying values ST . F ’1 (x) is the inverse of the cumulative normal distri-

bution function and CT (y) is the payo¬ function of the option.

The Monte Carlo simulation calculates the value of the integral in the following

way:

1 n

1. n independent random underlying values ST . . . ST are generated by com-

puting f (x, St , r, σ, T ’t) for a set of uniformly distributed pseudo random

numbers X1 , . . . , Xn .

i i

2. the option payo¬ CT (ST ) is calculated for each ST .

3. the value of the integral in (16.3) is then approximated by the arithmetic

mean of the option payo¬s:

n

1

¯ i

C= CT (ST )

n i=1

We will now derive an estimate of the approximation error of the arithmetic

1 n

mean. We assume that ST . . . ST are independent random underlying samples

of the g(ST |St , r, σ, T ’ t) density. Using this assumption we can conclude that

¯

C is a random variable with expected value

¯ = er(T ’t) Ct (St )

E[C]

352 16 Simulation based Option Pricing

Additionally we have to assume that the variance of the option payo¬s CT (ST )

is given by:

2

CT (ST )2 g(ST |St , r, σ, T ’ t)dST ’ E [CT (ST )]

Var [CT (ST )] = (16.4)

[0,∞]

exists. Then we get:

n

1 1

¯ i

Var C = 2 Var CT (ST ) = Var [CT (ST )] (16.5)

n n

i=1

1 n

because of the independence of ST , . . . , ST .

¯

The expected value of the random variable C equals the value of the inte-

gral er(T ’t) Ct (St ) and its variance converges to zero with increasing n. The

probability that the approximation error is greater than a ¬xed positive value

decreases to 0 with an increasing number n. A ¬rst estimation of the error is

¯

given by the Chebychev inequality for C,

1

n Var [CT (ST )]

¯

P |C ’ er(T ’t) Ct (St )| ≥ a ¤

a2

The bound given by this equation is rather imprecise since we do not make any

assumptions on the distribution of the random variable. Only the expected

value and the variance are used in the previous equation. According to the

¯

central limit theorem the distribution of C converges to a normal distribution

for n ’ ∞. It follows that the di¬erence between the approximation and the

¯

integral, C ’ er(T ’t) Ct (St ) is approximately normally distributed with mean 0

and standard deviation

Var [CT (ST )]

σC = (16.7)

¯

n

for large n. According to Boyle (1977) a value of n > 1000 is su¬ciently large

in order to use the normal distribution for error estimation purposes.

¯

We get the following equation if we assume that C ’ er(T ’t) Ct (St ) is normal

distributed:

u2

1

¯

C ’ er(T ’t) Ct (St ) ¤ =√ exp ’

P du (16.8)

2σC

2π ¯

’

16.1 Simulation techniques for option pricing 353

¯

If we choose k as a multiple of the standard deviation σC of C, then we get:

¯

¯

C ’ er(T ’t) Ct (St )

¯

C ’ er(T ’t) Ct (St ) ¤ kσC ¤k

P =P

¯

σC¯

k

u2

1

=√ exp ’ du

2

2π ’k

=p (16.9)

√

Given a ¬xed probability level p, the error converges to zero with O(1/ n).

The error interval holds for k = 1, 2, 3 with the respective probabilities p =

0.682, 0.955, 0.997

The con¬dence intervals for a given probability level depend on the standard

deviation of the payo¬ function CT (ST ):

σ CT = Var [CT (ST )] . (16.10)

In general, this standard deviation cannot be calculated with analytical meth-

ods. Therefore one calculates the empirical standard deviation σ and uses it

¯

as a proxy for the error bounds:

n

1 2

¯

i

CT (ST ) ’ C

σ=

¯ . (16.11)

n’1

k=1

Figure 16.1 shows the evolution of the absolute error of the price for a European

call option calculated by Monte Carlo methods compared with √ analyticthe

solution. One can observe that the error tends to zero with O (1/ n).

We would like to give some of the main properties of algorithms using Monte

Carlo techniques. First from (16.9) it follows that the error bound tends to zero

√

with O (1/ n) for a ¬xed probability level p. Second, the probability that a

√

¬xed error bound holds converges to 1 with O (1/ n), Mavin H. Kalos (1986).

Since these results hold independent of the dimension of the problem, which

a¬ects only the variance of the payo¬ function with respect to the Black-Scholes

risk neutral density, the Monte Carlo method is especially well suited for the

evaluation of option prices in multidimensional settings. Competing pricing

methods e.g ¬nite di¬erences have exponential growing computational costs in

354 16 Simulation based Option Pricing

Errors in MC Simulation

0.03 0.06 0.09 0.12 0.15

absolute error

0 500000 1000000

number of iterations

Figure 16.1. Absolute error of a European Call option price calculated

by Monte Carlo simulations vs. n’1/2

the dimension of the problem. Another advantage of the Monte Carlo pricing

method is the error estimate given by the empirical standard deviation which

can be computed with a small additional e¬ort.

The two most important drawbacks of the Monte Carlo simulation, mentioned

in literature are its small convergence speed compared to other techniques for

options on few underlyings and the di¬culties occurring for options with early

exercise possibilities. For example, American options giving the investor the

possibility to exercise the option at any time before and at maturity, are di¬cult

to price. To evaluate an American option means to ¬nd an optimal exercise

strategy which leads - using only basic Monte Carlo techniques - to a recursive

algorithm with exponential time-complexity. But more advanced techniques

using importance sampling methods show that Monte Carlo simulations can

be applied to evaluate American contracts, Broadie (2000).

16.1.3 Pricing path dependent European options on one

underlying

There are two categories of options. Path-independent options are options

whose payo¬ depends only on the underlying prices at maturity. Path-

16.1 Simulation techniques for option pricing 355

dependent options are options whose payo¬ depends on underlying price out-

comes St1 , . . ., Stm at several time points t1 ¤ . . . ¤ tm within the lifetime of

the respective option.

Within the group of path-dependent options one can distinguish options with

a payo¬ function depending on a continuously de¬ned path variable and op-

tions with a payo¬ function depending on a ¬xed number of underlying values.

The price of an option with many - usually equally spaced - exercise dates is

often approximated by the price of an option with a continuously de¬ned path

variable and vice versa.

Examples for path-dependent options are barrier options, lookback options,

and Asian options. The latter have a payo¬ function which is linked to the

average value of the underlying on a speci¬c set of dates during the life of

the option. One distinguishes two basic forms of Asian options: options on

the geometric mean (for which the price can be calculated with standard tech-

niques) and options on the arithmetic mean (for which the price can not be

determined using standard approaches). Asian options are frequently used in

commodity markets. The volatility of the underlying prices of the commodities

is usually very high so that prices for vanilla options are more expensive than

for comparable Asian-style options.

16.1.4 Pricing options on multiple underlyings

In this section we show how to extend the Monte Carlo simulation technique

to higher dimensions. The problem is not only that one has to deal with higher

dimensional integrals, but also that one has to incorporate the underlying cor-

relation structure between the considered securities. In our framework we need

the covariance matrix of the log returns on an annual basis.

In general, a basket option is an option on several underlyings (a basket of

underlyings). Basket options can be European-, American or even Asian-style

options. Normally, the average of the underlying prices is taken to calculate

the price of the basket option, but sometimes other functions are used.

The advantage of the usage of basket options instead of a series of one dimen-

sional options is that the correlation between securities is taken into account.

This may lead to better portfolio hedges. We will look at a basket option on

¬ve underlyings where the underlying price of the best security in the basket

is taken to calculate the option price.

356 16 Simulation based Option Pricing

16.2 Quasi Monte Carlo (QMC) techniques for

option pricing

16.2.1 Introduction to Quasi Monte Carlo techniques

QMC methods can be considered as an alternative to Monte Carlo simulation.

Instead of (pseudo) random numbers, Quasi Monte Carlo algorithms use the

elements of low discrepancy sequences to simulate underlying values.

The discrepancy of a set of points P ‚ [0, 1]s measures how evenly these points

are distributed in the unit cube. The general measure of discrepancy is given

by:

A(B; P )

’ »s (B)

Dn (B; P ) := sup (16.12)

n

B∈B

where A(B; P ) is the number of points in P belonging to B, »s (B) is the

Lebesgue measure of the set B, B is a family of Lebesgue measurable subsets

of [0, 1]s , and n is the number of elements in P .

The discrepancy of a set is the largest di¬erence between the number of points

in a subset and the measure of the subset. If we de¬ne B to be the family J of

s

subintervals i=1 [0, ui ), then we get a special measure, the star-discrepancy:

—

Dn (P ) := Dn (J ; P ) (16.13)

16.2.2 Error bounds

For the star-discrepancy measure and reasonable assumption on the nature of

the function that has to be integrated an upper bound on the error is given by

the following theorem:

THEOREM 16.1 (Koksma-Hlawka) If the function f is of ¬nite variation

V (f ) in the sense of Hardy and Krause, then the following equation holds for

all sets of points {x1 , . . . , xn } ‚ I s = [0, 1]s

n

1 —

f (xi ) ’ f (u)du ¤ V (f )Dn (x1 , . . . , xn ) (16.14)

n Is

i=1

16.2 Quasi Monte Carlo (QMC) techniques for option pricing 357

A proof is given in Niederreiter (1992).

This means that the error is bounded from above by the product of the

variation V (f ), which in our case is model and payo¬ dependent and the star-

discrepancy of the sequence. The bound cannot be used for an automatic error

estimation since the variation and the star-discrepancy cannot be computed

easily. It has been shown though that sequences exist with a star-discrepancy

of the order O(n’1 (ln n)s ). All sequences with this asymptotic upper bound

are called low-discrepancy sequences Niederreiter (1992). One particular

low-discrepancy sequence is the Halton sequence.

16.2.3 Construction of the Halton sequence

We start with the construction of the one-dimensional Halton sequence within

the interval [0, 1]. An element of this sequence is calculated by using the fol-

lowing equation:

∞

nk,i p’k’1

xi = (16.15)

k=0

with i > 0, p = 2 and nk,i determined by the following equation:

∞

nk,i pk ; 0 ¤ nk,i < p; nk,i ∈ N

i= (16.16)

k=0

Note that with the above equation nk,i is a function of i and takes values only

in {0; 1}. To illustrate the algorithm we calculate the ¬rst three points.

i = 1: n0,1 = 1, nk,1 = 0 for every k > 0

i = 2: n1,2 = 1, nk,2 = 0 for every k = 1

i = 3: n0,3 = n1,3 = 1, nk,3 = 0 for every k > 1

Therefore we get the sequence 1/2, 1/4, 3/4, 1/8, 5/8, .... The extension of this

construction scheme to higher dimensions is straightforward. For every dimen-

sion j = 1, . . . , d we de¬ne xj by

i

∞

xj nk,i (j)p’k’1

= (16.17)

i j

k=0

358 16 Simulation based Option Pricing

with pj is the jth smallest prime number and nk,i (j) is calculated as follows:

∞

nk,i (j)pk ; 0 ¤ nk,i (j) < pj ; nk,i (j) ∈ N ∀j

i= (16.18)

j

k=0

By using p1 = 2, p2 = 3 we get the following two-dimensional Halton sequence:

(1/2; 1/3), (1/4; 2/3), .... In contrast to grid discretization schemes like i/n i =

1, ..., n low-discrepancy sequences ¬ll the integration space in an incremental

way avoiding the exponential growth of grid points of conventional schemes.

XploRe provides quantlets to generate pseudo random numbers and low dis-

crepancy sequences. For the generation of the pseudo random numbers we use

erg = randomnumbers (seqnum,d,n)

generates n pseudo random vectors of dimension d

where seqnum is the number of the random generator according to Table 16.1,

d is the dimension of the random vector and n the number of vectors generated.

0 Park and Miller with Bays-Durham shu¬„e

1 L™Ecuyer with Bays-Durham shu¬„e

2 Knuth

3 generator from G. Marsaglia et al. Marsaglia (1993)

4 random number generator of your system

5 generator from ACM TOMS 17:98-111

6 multiply with carry gen. (Marsaglia) Marsaglia (1993)

Table 16.1. Random generator that can be used in XploRe

The generation of low discrepancy sequences is provided by

erg = lowdiscrepancy (seqnum,d,n)

generates the ¬rst n low discrepancy sequence vectors of dimen-

sion d

where seqnum is the number of the low discrepancy sequence according to Table

16.2.

16.2 Quasi Monte Carlo (QMC) techniques for option pricing 359

0 Halton sequence

1 Sobol sequence

2 Faure sequence

3 Niederreiter sequence

Table 16.2. Low-discrepancy sequences available in XploRe,

(Niederreiter, 1992) .

16.2.4 Experimental results

Figure 16.2 shows that two dimensional Halton points are much more equally

spaced than pseudo random points. This leads to a smaller error at least for

“smooth” functions.

First 1000 random points First 1000 Halton points

1

1

dimension 2

dimension 2

0.5

0.5

0

0

0 0.5 1 0 0.5 1

dimension 1 dimension 1

Figure 16.2. 1000 two-dimensional pseudo random points

vs. 1000 Halton points XFGSOPRandomNumbers.xpl,

XFGSOPLowDiscrepancy.xpl

The positive e¬ect of using more evenly spread points for the simulation task

is shown in Figure 16.3. The points of a low-discrepancy sequence are designed

in order to ¬ll the space evenly without any restrictions on the independence

of sequence points where as the pseudo random points are designed to show no

statistically signi¬cant deviation from the independence assumption. Because

of the construction of the low discrepancy sequences one cannot calculate an

360 16 Simulation based Option Pricing

Errors in QMC vs. MC Simulation

-5

log absolute error

-10

0 500000 1000000

number of iterations

Figure 16.3. Absolute error of a random sequence and the Halton

sequence for a put option

empirical standard deviation of the estimator like for Monte Carlo methods

and derive an error approximation for the estimation. One possible way out

of this dilemma is the randomization of the low-discrepancy sequences using

pseudo random numbers i.e. to shift the original quasi random numbers with

pseudo random numbers Tu¬n (1996). If x1 , . . . , xn are scalar elements of a

low-discrepancy sequence X then we can de¬ne a new low discrepancy sequence

xi + xi + <= 1

W ( ) = {y1 , . . . , yn } with yi = (16.19)

(xi + ) ’ 1 otherwise

for a uniformly distributed value . Then we can calculate an empirical stan-

dard deviation of the price estimates for di¬erent sequences W ( ) for indepen-

dent values which can be used as a measure for the error. Experiments with

payo¬ functions for European options show that this randomization technique

reduces the convergence rate proportionally.

The advantage of the Quasi Monte Carlo simulation compared to the Monte

Carlo simulation vanishes if the dimension increases. Especially the compo-

nents with a high index number of the ¬rst elements in low-discrepancy se-

quences are not evenly distributed Niederreiter (1992). Figure 16.4 shows that

the 49th and 50th component of the ¬rst 1000 Halton points are not evenly

distributed. But the result for the ¬rst 10000 points of the sequence shows that

the points become more evenly spread if the number of points increases.

16.3 Pricing options with simulation techniques - a guideline 361

First 1000 Halton Points First 10000 Halton Points

1

1

dimension 50

dimension 50

0.5

0.5

0

0

0 0.5 1 0 0.5 1

dimension 49 dimension 49

Figure 16.4. The ¬rst 1000 and 10000 Halton points of dimension 49

and 50 XFGSOPLowDiscrepancy.xpl

However by using the Brownian Bridge path construction method we can limit

the e¬ect of the high dimensional components on a simulated underlying path

and the corresponding path variable for the most common path dependent

options, Moroko¬ (1996). This method start with an empty path with known

start value and calculates at each step the underlying value for a time point

with maximum time distance to all other time points with known underlying

value until the whole path is computed. Experimental results show that we

can still get a faster convergence of the QMC simulation for options up to 50

time points if we apply this path construction method.

16.3 Pricing options with simulation techniques -

a guideline

In this section we would like to give a short guideline how to price exotic

options with Monte Carlo and Quasi Monte Carlo simulation techniques

within the framework described above. Furthermore we give some indications

about the limits of these techniques.

362 16 Simulation based Option Pricing

16.3.1 Construction of the payo¬ function

As a ¬rst step we have to de¬ne the payo¬ function corresponding to our

option product. Within the methods de¬ned in the quantlib finance we have

to consider three di¬erent cases.

One underlying + path independent

In this case the payo¬ function is called by the pricing routine with the sim-

ulated underlying value at maturity as the single argument. It calculates the

corresponding payo¬ and returns this value. We have de¬ned the payo¬ func-

tion for a put option with strike price 100 as an example for a one dimensional

payo¬ function.

Several underlying + path independent

For options whose payo¬ depends on the underlying values of several assets at

maturity, we have to de¬ne a payo¬ function on the vector of the underlying

values at maturity. An example for such an option is an exchange option that

permits to swap a de¬ned share with the best performing share in a basket. Its

payo¬ function is given by:

CT ((ST , . . . , ST )) = max{0, ±i (ST ’ K i ) + 55 ’ ST |i = 1, .., 5}

1 5 i 3

One underlying + path dependent

The third category of option types that are captured are path dependent op-

tions on one underlying. The payo¬ function of these options depends on the

underlying values at several ¬xed time points during the lifetime of the option.

Payo¬ functions for these contracts are called with a vector of underlying val-

ues whose ith element is the underlying value at the time ti which has to be

speci¬ed in the model.

16.3.2 Integration of the payo¬ function in the simulation

framework

After de¬ning the payo¬ function in XploRe we can start to calculate a price

estimate with the help of the appropriate simulation routine. In the one di-

mensional case we just have to call

16.3 Pricing options with simulation techniques - a guideline 363

erg = BlackScholesPathIndependent1D (s0,r,vola,dt,opt,

itr,gen)

MC estimation of the option price for a path independent option.

erg = BlackScholesPathIndependent1DQMC (s0,r,vola,dt,opt,

itr,gen)

QMC estimation of the option price for a path independent

option.

to get a price estimate and for the Monte Carlo case an empirical standard

deviation with respect to a start price of s0, a continuous risk free interest

rate of r, a volatility vola, a time to maturity of dt years, the payo¬ function

opt, sample size itr and the random/low-discrepancy generator with number

gen. Table 16.1 shows the random number generators and table 16.2 the low-

discrepancy generators that can be used. An application of these routines for

a Put option can be found in XFGSOP1DPut.xpl.

Pricing path-dependent options is only slightly more complicated. Here we

have to de¬ne the vector of time points for which underlying prices have to

be generated. This vector replaces the time to maturity used to price path

independent options. Then we can apply one of the following methods to

compute a price estimate for the path dependent option

erg = BlackScholesPathDependent1D (s0,r,vola,times,opt,

itr,gen)

MC estimation of the option price for path-dependent options.

erg = BlackScholesPathDependent1DQMC (s0,r,vola,times,opt,

itr,gen)

QMC estimation of the option price for path-dependent options,

with:

with respect to the start price s0, the continuous risk free interest rate r, the

volatility vola, the time scheme times, the payo¬ function opt, sample size

itr and the random/low-discrepancy generator with number gen, as given in

Tables 16.1 and 16.2. Using the above quantlets, we calculate the price of an

Asian call option in XFGSOP1DAsian.xpl.

364 16 Simulation based Option Pricing

In the case of multidimensional options we have to de¬ne a start price vector

and a covariance matrix instead of a single underlying price and volatility value.

Then we can call one of the multi-dimensional simulation routines:

erg = BlackScholesPathIndependentMD (s0,r,vola,dt,opt

,itr,gen)

MC estimation of the option price in the multidimensional Black

Scholes model

erg = BlackScholesPathIndependentMDQMC (s0,r,vola,dt,opt

,itr,gen)

QMC estimation of the option price in the multidimensional

Black Scholes model

with respect to the m dimensional start price vector s0, the continuous risk free

interest rate r, the m—m covariance matrix vola, the time to maturity dt, the

payo¬ function opt, the number of iterations itr and the generator number

gen according to the generators in Tables 16.1 and 16.2. Both quantlets are

illustrated in XFGSOPMD.xpl.

If in addition a dividend is paid during the time to maturity, we can use the

following two quantlets to calculate the option prices.

erg = BlackScholesPathIndependentMDDiv (s0,r,div,vola

,dt,opt,itr,gen)

MC estimation of the option price in the multidimensional Black

Scholes model

erg = BlackScholesPathIndependentMDDivQMC (s0,r,div,vola

,dt,opt,itr,gen)

QMC estimation of the option price in the multidimensional

Black Scholes model

The additional argument div is a m dimensional vector of the continuously paid

dividends. An application of these functions for our basket option is provided

in XFGSOPMDDiv.xpl.