Valuation of American Basket Options using Quasi-Monte Carlo Methods

d-fine GmbH Christ Church College University of Oxford

A thesis submitted in partial fulfillment for the MSc in Mathematical Finance September 26, 2009

Abstract

Title:

Valuation of American Basket Options using Quasi-Monte Carlo Methods

Author:

d-fine GmbH

Submitted for: MSc in Mathematical Finance Trinity Term 2009 The valuation of American basket options is normally done by using the Monte Carlo approach. This approach can easily deal with multiple random factors which are necessary due to the high number of state variables to describe the paths of the underlyings of basket options (e.g. the German Dax consists of 30 single stocks). In low-dimensional problems the convergence of the Monte Carlo valuation can be speed up by using low-discrepancy sequences instead of pseudorandom numbers. In high-dimensional problems, which is definitely the case for American basket options, this benefit is expected to diminish. This expectation was rebutted for different financial pricing problems in recent studies. In this thesis we investigate the effect of using different quasi random sequences (Sobol, Niederreiter, Halton) for path generation and compare the results to the path generation based on pseudo-random numbers, which is used as benchmark. American basket options incorporate two sources of high dimensionality, the underlying stocks and time to maturity. Consequently, different techniques can be used to reduce the effective dimension of the valuation problem. For the underlying stock dimension the principal component analysis (PCA) can be applied to reduce the effective dimension whereas for the time dimension the Brownian Bridge method can be used. We analyze the effect of using these techniques for effective dimension reduction on convergence behavior.

To handle the early exercise feature of American (basket) options within the Monte Carlo framework we consider two common approaches: The Threshold approach proposed by Andersen (1999) and the Least-Squares Monte Carlo (LSM) approach suggested by Longstaff and Schwartz (2001). We investigate both pricing methods for the valuation of American (basket) options in the equity market.

Contents 1 Introduction

1

2 Basket Options

4

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Different Types of Basket Options . . . . . . . . . . . . . . . . . . . .

4

2.2.1 2.2.2

Plain Vanilla Basket Options . . . . . . . . . . . . . . . . . . Exchange Options . . . . . . . . . . . . . . . . . . . . . . . . .

4 5

2.2.3

Asian Basket Options . . . . . . . . . . . . . . . . . . . . . . .

6

American-style Basket Options . . . . . . . . . . . . . . . . . . . . .

7

3 Monte Carlo Simulation 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 8

2.3

3.2

Mathematical Foundations . . . . . . . . . . . . . . . . . . . . . . . .

9

3.2.1

Strong Law of Large Numbers . . . . . . . . . . . . . . . . . .

9

3.2.2 3.2.3

Error Estimation . . . . . . . . . . . . . . . . . . . . . . . . . Central Limit Theorem . . . . . . . . . . . . . . . . . . . . . .

10 10

3.3

Pseudo-Random Numbers . . . . . . . . . . . . . . . . . . . . . . . .

11

3.4

Quasi-Random Numbers . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.4.1 3.4.2

Low Discrepancy and Integration Error . . . . . . . . . . . . . Halton Sequences . . . . . . . . . . . . . . . . . . . . . . . . .

14 16

3.4.3

Sobol Sequences . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3.4.4

Niederreiter Sequences . . . . . . . . . . . . . . . . . . . . . .

19

The curse of dimensionality . . . . . . . . . . . . . . . . . . . . . . .

22

3.5

4 Effective Dimension Reduction

23

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4.2

Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . .

23

4.3

Brownian Bridge . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

i

5 Valuation of American Basket Options

28

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

5.2 5.3

Model Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Monte Carlo Path Construction . . . . . . . . . . . . . . . . . . . . .

29 30

5.4

Valuation of Early Exercise Feature . . . . . . . . . . . . . . . . . . .

33

5.4.1

Least-Squares Approach . . . . . . . . . . . . . . . . . . . . .

33

5.4.2

Threshold Approach . . . . . . . . . . . . . . . . . . . . . . .

35

6 Comparison of Numerical Methods

37

6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

6.2

Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

6.3 6.4

Software Implementation . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38 40

6.4.1

Classical Put Options . . . . . . . . . . . . . . . . . . . . . . .

40

6.4.2

Exchange Options . . . . . . . . . . . . . . . . . . . . . . . . .

43

6.4.3

Basket Options . . . . . . . . . . . . . . . . . . . . . . . . . .

45

7 Conclusions

49

A Data

51

Bibliography

51

ii

List of Tables 6.1

Statistics for different paths generation methods for American Dax call valuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

A.1 Description Dax constituents . . . . . . . . . . . . . . . . . . . . . . A.2 Correlation matrix for Dax constituents for April 30, 2009 . . . . . .

51 52

iii

List of Figures 3.1 3.2

Two-dimensional projections of pseudo-random sequences . . . . . . . Two-dimensional projections of Halton sequences . . . . . . . . . . .

13 17

3.3

Two-dimensional projections of Sobol sequences . . . . . . . . . . . .

20

3.4

Two-dimensional projections of Niederreiter sequences . . . . . . . . .

21

4.1

Brownian bridge construction . . . . . . . . . . . . . . . . . . . . . .

26

6.1 6.2

Convergence properties plain-vanilla American put . . . . . . . . . . Convergence properties American exchange options using Underlying-

42

Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

6.3

Convergence properties American Dax call for different random number generators using Least-Squares approach . . . . . . . . . . . . . .

48

A.1 Convergence properties of American put for different random number generators using Least-Squares approach . . . . . . . . . . . . . . . .

53

A.2 Convergence properties of American put for different random number generators using Threshold approach . . . . . . . . . . . . . . . . . .

54

A.3 Convergence properties of American exchange options for different random number generators . . . . . . . . . . . . . . . . . . . . . . . . . .

55

A.4 Convergence properties American Dax call for pseudo-random numbers and Sobol sequences . . . . . . . . . . . . . . . . . . . . . . . . .

56

A.5 Convergence properties American Dax call for Niederreiter and Halton sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iv

57

Chapter 1 Introduction The pricing and optimal exercise of options with early exercise features is one of the most challenging problems in mathematical finance. The problem becomes even more complex, if the model for the economy has more than one state variable, i.e. more than one factor affecting the value of the option as well as the optimal exercise moment. American basket options are a prominent example for derivatives with early exercise features and multiple state variables. Due to their construction principle basket options generally represent an advantageous alternative to hedge risky positions of several assets, since trading only one basket option instead of several single-asset options decreases transaction costs. Usually lattice methods or finite differences will be used to price American or Bermudan style options. However, these techniques are inefficient for high dimensional problems (i.e. problems with more than three state variables) and they are very difficult to apply on path-dependent options. For higher dimensions it would be relatively cheap to use simulation techniques since its computational cost does not increase exponentially as compared to other methods and they are relatively simple to apply. However, for a long time, simulation techniques seemed unapplicable to options with early exercise features. This is due to their forward-construction principle and their path-by-path generation. Since contingent claims with early exercise features are traded in all important derivative markets, many suggestions have been made during the last years to price such options by simulation, and in particular to determine the optimal early-exercise strategy. Examples can be found in Tilley (1993), Barraquand and Martineau (1995), Carriere (1996), Broadie and Glasserman (1997), Andersen (1999), Longstaff and Schwartz (2001), or Ib´anez and Zapatero (2004). In this case the starting point will be the application of standard Monte Carlo simulation with pseudo-random numbers, eventually combined with some variance 1

reduction techniques such as antithetic variables. However, using pseudo-random numbers can be quite slow, due to the relatively poor convergence rate of O(M −1/2 ) for M sample paths. Here, the application of Quasi-Monte Carlo simulation (QMC), e.g. using Sobol, Niederreiter or Halton sequences, can improve the performance of Monte Carlo simulations dramatically, since in optimal cases the convergence rate of QMC is O(log(M )d M −1 ), where d represents the dimension of the integration problem. Therefore, for small values of d the convergence is due to the construction principle of low-discrepancy sequences significantly higher as for standard Monte Carlo. However, for higher dimensional problems this advantage will diminish, since the convergence strongly depends on the dimension of the integration problem for which reason the application of QMC to price American basket option seems not advantageous. Especially for American options it cannot be expected that QMC works, since due to regression or optimization over the sample paths the sample points have interactions. However, Chaudhary (2005) discovered a significant speedup in the numerical results for the Least-squares approach proposed by [25]. Therefore we also expect to find this improvement for the convergence of the valuation of American basket options when applying QMC. Moreover, earlier empirical studies for finance problems, see among others [30, 29, 12], have shown in fact significant improvements using QMC even in high dimensional financial pricing problems up to dimension 360. The main reason is that many practical problems in mathematical finance have a low effective dimension – a notion introduced by [12] – i.e. either the integrand depends only of several variables, or the integrand can be approximated by a sum of functions with each depending only on a small number of variables at the same time. Several techniques to reduce the effective dimension have been successfully applied for derivative pricing problems. The most prominent are the Brownian bridge method for path construction (BB) first applied by Moskowitz and Caflisch (1996) and Caflisch et al. (1997) and the principal component analysis (PCA) first used by Acworth et al. (1997) to identify the underlyings which contribute the largest part of the total variance. Besides this empirical studies Wang and Sloan (2007) provides a theoretical analysis of the possible convergence order of QMC used in conjunction with BB or PCA. However, all these studies only analyze the convergence properties of QMC with BB or with PCA, but, to the best of our knowledge, not the combined application of QMC with both – BB and PCA. In our analysis we consider the valuation of American basket options as a prominent example for a high-dimensional derivative pricing problem, whereby we reduce 2

the effective dimension by applying the Brownian bridge method (BB) and the Principal Component Analysis (PCA) to reduce the effective time and underlying stock dimension. We focus our analysis on the convergence properties and numerical stability by applying the different techniques as well as to develop a path construction method which takes both effective dimension reduction techniques into account. Moreover, we analyze two common approaches which can easily be adapted for high-dimensional pricing problems to incorporate the early exercise features of American-style options into the Monte Carlo method: The Threshold approach proposed by Andersen (1999) and the Least-Squares (LS) approach suggested by Longstaff and Schwartz (2001).1 The thesis is structured as follows: The next chapter contains a general description of basket options. The theoretical framework for Monte Carlo simulation is given in chapter 3. Chapter 4 describes the different techniques to reduce the effective dimension of high dimensional pricing problems. The model framework for the valuation of American basket options is outlined in chapter 5: this chapter contains the construction of the Monte Carlo paths in this framework and takes into account the effective dimension reduction techniques as well as the valuation approaches to determine the early exercise value of an American basket option are contained . Chapter 6 presents the numerical results from the application of these methods to test basket options. A final chapter concludes with a brief summary of the results and some recommendations for further applications of Quasi-Monte Carlo methods based on our results.

1

The hedging of American basket options and the calibration of the pricing model will not be considered in this thesis. We assume in our study, that all necessary market parameters are correctly given.

3

Chapter 2 Basket Options 2.1

Introduction

In a broader sense we summarize under the class ”basket options” all derivatives whose payoff depends on more than one underlying. The underlying basket normally contains different equities, commodities, or currencies. In the last years a wide range of different basket options have emerged with special features such as early exercise and path dependency. We describe the most prominent types of basket options in following sections. Thereby we assume a basket of equities as underlying, since for the analysis in this thesis we only consider this type of underlying. For the pricing of basket options an appropriate standard model (such as Black-Scholes for equities) for each individual asset and for the correlation between the underlying stochastic drivers is used. The pricing model is provided in Section 5.2.

2.2 2.2.1

Different Types of Basket Options Plain Vanilla Basket Options

Plain vanilla basket options are basket options in the true sense and give the holder the right to buy (call) or sell (put) a specified basket of equities at an agreed price (strike price) at the maturity date of the option (European style). The payoff, P O, of a basket put option as function of the underlying equities, Si , is given as Put P OBasket (S1 , ..., SU )

= max(K −

U X

ωi Si , 0)

(2.1)

i=1

where K is the strike price, U the number of underlying equities and ωi the fixed weight of each asset in the basket, and Si the corresponding prices at the maturity of the option. In case of physical settlement, the weights must be integers because 4

fractions of stocks cannot be delivered. For cash settlement options, which is usually the case, this constraint can be neglected. At first glance index options could also be considered as basket options, or vice versa, since an index is nothing else as a basket of equities. But this is definitely not the case in practical applications. The index weights normally change every day depending on the free float, current stock price, etc. whereas the weights of a basket options are fixed in the option contract. Therefore, for example, the price of a DAX index option for a special day is not the same as the price of a basket option, where the basket consists of the index constituents with weights valid for that special date, and all other option parameters are the same. Nevertheless, basket options are often priced by modeling the basket value as single underlying, where the volatility of the fictitious underlying representing the basket can be derived from the covariance matrix of the constituents. In that case, standard option pricing formulas can be applied to value the option and to calculate the greeks with respect to the basket value.

2.2.2

Exchange Options

Exchange options1 give the holder the right to exchange one asset for another, i.e. the payoff depends on two assets. Obviously the price of such an option also depends on the correlation between them. P OExch = max(S1 − S2 , 0)

(2.2)

Compared to plain vanilla options an exchange option can be interpreted as put option on S2 with strike S1 or as call option on S1 with strike S2 . Margrabe (1978) derived a closed form solution for the European exchange option value V in the Black-Scholes framework for the two-asset case, which in general is not possible for basket options VExch (S1 , S2 , t) = S1 e−D1 (T −t) N (d1 ) − S2 e−D2 (T −t) N (d2 )

(2.3)

where log d1 =

1

  S1 S2

+ (D2 − D1 + 0.5σ 2 )(T − t) p p and d2 = d1 − σ (T − t), σ (T − t) q and σ = σ12 + σ22 − 2ρ12 σ1 σ2

Exchange options are also known as Margrabe options or outperformance options.

5

(2.4)

(2.5)

Here, Di represents the respective dividend yield of asset i, t the current time, T the maturity of the exchange option, N (di ) the cumulative normal distribution function of di , σi the volatility of asset i and ρi,j the correlation between the assets i and j, i, j{1, 2}.

2.2.3

Asian Basket Options

The payoff of Asian Basket Options depends on the average value of the underlying basket over a specific period during the lifetime of the option. Thereby, the average used in the calculation is usually defined as arithmetic or geometric average. For practical and legal reasons, the dates to calculate the average are sampled discretely, e.g. only closing prices every Wednesday are used to calculate the average. In principle continuous sampling of data would be thinkable as well. Besides the typical features of options Asian options can be further classified into ”average rate” or ”average strike” option. The payoff of an average rate option is based on the difference between the basket average value and a fixed strike Call P OAsian Rate (S1 , ..., SU ) = max(A − K, 0)

with A=

T X U X

(2.6)

ωi Si (t)

(2.7)

t=0 i=1

is the arithmetic average over the sampling period and B =

PU

i=1

ωi Si (T ) the value

of the basket at time T . The payoff of an average strike option is based on the difference between the current basket average value and the basket average value over the sampling period Call P OAsian Strike (S1 , ..., SU ) = max(B − A, 0).

(2.8)

Continuously-sampled, geometric average rate options can be valued using a modified version of the Black-Scholes formula, if the underlying value follows a lognormal process. In this case the geometric average is also lognormal distributed. Kemna and Vorst (1990) derived this solution by adjusting the volatility and the cost of carry term of the Black-Scholes formula. However, this formula can only be used for Asian basket options, if the basket value itself is modeled as underlying (compare 2.2.1). For arithmetic average rate options a closed form solution cannot be given but there are different approximation formulas (e.g. Turnbull and Wakeman (1991) or Levy (1992)) to calculate the value of this option type. 6

2.3

American-style Basket Options

American-style basket options are special in the sense that they can be exercised at any time before their maturity date. This feature is quite common in derivatives markets and not special for basket options, e.g. most traded single equity options are American style. As the combination of basket options with early exercise features is one of the hardest problems to solve in mathematical finance we give some general remarks about American-style options. Due to the early exercise feature American-style options are more valuable than their European-style counterparts which can only be exercised at maturity. This can be directly concluded from the early exercise criterion: this says that the option should be exercised as soon as the option payoff is greater than the value of continuation; mathematically: P O(t) ≥ Et [P O(t∗ )]

(2.9)

where E represents the expectation operator, and t∗ the optimal future exercise point of the option. Therefore, an American-style option is at least as valuable as its European-style counterpart. The optimal exercise boundary of an option is then given by all points where equation (2.9) holds. Therefore, the optimal exercise boundary for basket options is a function of time t and all underlying values Si , with i = 1, ..., U . If the early exercise boundary is known, valuation of American-style options is straightforward. But this is not the case in reality. Thus, the challenging part of valuing American-style options is the determination of the early exercise boundary, which has to be done numerically, since no closed form solution can be derived for the free boundary problem. Examples for approximation formulas to value American-style options can be found in Geske and Johnson (1984), Barone-Adesi and Whaley (1987) or in a quite recent article from Zhu (2006). But, for basket options with several underlyings these approximations can generally not be applied, especially when the option payoff depends at the end on one single equity from the basket. American exchange options are one of the special cases where a closed form approximation was derived by Bjerksund and Stensland (1993) based on the approximation formulas for American single-equity options.2

2

This is contrary to Margrabe (1978) who has argued that American exchange options are not more valuable than European ones.

7

Chapter 3 Monte Carlo Simulation 3.1

Introduction

Monte Carlo simulation is a numerical method that is advantageous for many mathematical problems where no closed-form solution is available, especially for high dimensional problems. Due to its nature Monte Carlo simulation is mainly used to simulate random variables or stochastic processes, that can be described by some given probability density function. One of the most prominent applications for Monte Carlo is that of numerical integration, i.e. the solution of definite integrals with complicated boundary conditions.1 This concept of numerical integration by Monte Carlo simulation was originally introduced for option pricing theory by Boyle (1977). To illustrate the idea for numerical integration consider the one-dimensional integral Z f (x)p(x)dx = Ep [f (x)] = I(f ) (3.1) Ω

where f (x) is an arbitrary function and p(x) is a probability density function with Z p(x)dx = 1. (3.2) Ω

Thereby, Ω denotes the range of integration and Ep [f ] is the expected value of f with respect to p. The problem of numerical integration can be interpreted as the approximation of an expected value. To estimate the value of I(f ) a large number M of independent and identically distributed (i.i.d) sample values xi of the probability 1

Beside numerical integration another powerful and very popular application for Monte Carlo is stochastic optimization. Optimization problems are by nature not stochastic. However, to find the global optimum of a high-dimensional optimization problem where the objective function is nondifferentiable such that usual gradient methods fail, stochastic optimization can be applied. Thereby, the sample points of the multi-dimensional space will be drawn randomly and the objective function is then evaluated with the drawn sample points. This procedure will be repeated as long as new optima are found outside a given error bound compared to the last found optimum.

8

density function p(x) are drawn by using a random number generator. The estimate ˆ ) of I(f ) is calculated as the average over all evaluations of f (x) at the random I(f sample points and is given by M 1 X IM (f ) = f (xi ). M i=1

(3.3)

Benefits of using Monte Carlo simulation are the very basic mathematics, the simple modeling of correlations, and the flexibility of the method. In addition, the increased availability of more powerful computers and software packages during the recent years makes it easier to implement Monte Carlo methods and therefore enhances its popularity. The main disadvantage is that the method is slow to solve partial differential equations for low-dimensional problems up to three or four dimensions compared to other numerical methods like finite-differences. Moreover, applying Monte Carlo in Mathematical Finance to free boundary problems, like American-style options, is very challenging and seemed impracticable up to the mid 90s. The rest of this chapter gives some basics of Monte Carlo simulation without incorporating financial engineering applications to show the generality of this method.

3.2

Mathematical Foundations

3.2.1

Strong Law of Large Numbers

The strong law of large numbers (SLLN) ensures that the estimate as given in (3.3) converges almost surely to the true value of the integral in the sense that M 1 X a.s. f (xi ) −−→ I(f ) M →∞ M i=1

lim IM (f ) = lim

M →∞

i.e.

M 1 X f (xi ) = I(f ) P lim M →∞ M i=1

(3.4)

! =1

(3.5)

Almost surely means that other values are theoretically possible for a given sample, whereas for increasing sample sizes, the probability of other values converges toward zero. On the opposite side this also means that IM (f ) is itself a random variable and not a deterministic number.

9

3.2.2

Error Estimation

Concluding from the SLLN, IM (f ) is an unbiased estimator for I(f ), i.e. E[IM (f )] = I(f ) which converges to the true value almost surely. We can define the integration error, εM (f ), as εM (f ) = IM (f ) − I(f )

(3.6)

so that it is a zero-mean random variable. Given this definition, the estimation bias is defined as Bias(M ) = E[εM (f )] and the root mean square error (RMSE) as p RMSE(M ) = E[εM (f )2 ].

(3.7)

(3.8)

Empirically, the RMSE is often estimated as the standard deviation of IM (f ), since the true value of the integral (3.1) is not known.

3.2.3

Central Limit Theorem

The central limit theorem (CLT) allows us to characterize the size and statistical properties of the Monte Carlo integration error. Starting with the unbiased estimator IM (f ) from (3.3) with xi i.i.d. and assuming Var(xi ) = σx2 < ∞ and Var(f (xi )) = σf2 < ∞, then the CLT asserts that as the sample size M increases, the standardized estimator IM (f ) − I(f ) √ (3.9) σf / M converges in distribution to the standard normal distribution, which is often written as2

√ IM (f ) − I(f ) εM (f ) M √ lim ⇒ N (0, 1), = lim M →∞ M →∞ σf σf / M

(3.10)

which is equivalent to

σf εM (f ) ∼ √ z , as M → ∞ (3.11) M where z ∼ N (0, 1) is a zero-mean unit-variance normally distributed random variable, ⇒ denotes convergence in distribution and N (x, y 2 ) the normal distribution with mean x and variance y 2 . 2

The precise statement for convergence in distribution is given by ! Z Z −s2 IM (f ) − I(f ) 1 √ lim P ≤Z = √ e 2 ds. M →∞ 2π −∞ σf / M

10

Therefore, the CLT ensures that the standard error of the estimation tends to √ zero at rate 1/ M , i.e. to lower the integration error by a factor of ten the sample size has to be increased by hundred. The convergence of Monte Carlo simulation can sometimes be improved by application of variance reduction techniques such as antithetic variables, stratified sampling, control variates or moment matching methods. See Boyle et al. (1997) or Wilmott (1999) for a detailed description of these methods in mathematical finance. As seen above, the CLT is usually stated as formula for the integration error for given sample size M and known variance σf2 . However, in many practical applications, the exact value of the integral is unknown, so errors and the variance cannot be determined. In this case, the CLT formula will be used conversely, i.e. the empirical ˜ points, where M ˜ < M, sample variance, σ ˜f2 , is determined from evaluation f with M and then the sample size for a required level of accuracy, ε¯, with confidence c is determined by M = ε¯2 σ ˜f2 N −1 (c)2 ,

(3.12)

where N −1 is the inverse of the cumulative normal distribution function.

3.3

Pseudo-Random Numbers

The core of nearly each application of Monte Carlo simulation is a sequence of random numbers used to run the simulation. This random numbers have to satisfy a given probability density function, if the number of draws goes to infinity. This is typically done by transformation of independent, uniformly distributed random numbers drawn from the interval (0,1)3 . To transform the random numbers into the desired distribution function, usually the inverse cumulative distribution function is used.4 Most programming languages and spreadsheets provide a uniform random-number generator. Note, however, that these numbers are not random per se; they are pseudorandom, since computers follow deterministic algorithms to produce these numbers. Thereby, the pseudo-random sequences are made to have as many properties of ’true’ 3

Here, 0 and 1 are explicitly excluded from the interval, since for most probability distribution these points maps to either −∞ or ∞, which leads to numerical problems. 4 If there is no closed-form solution for the inverse cumulative distribution function, mostly an accurate approximation is used. For standard normal variables some other transformation methods are also available, e.g. Box-Muller or Marsaglia method.

11

random numbers as possible, i.e. they are sufficiently good at mimicking genuine random numbers with respect to fundamental statistical tests5 . The linear congruential method is the most commonly used method to generate uniformly distributed pseudo-random number sequences. The basic idea is to produce integer values ci on a given interval [0, C −1], where C is some constant, and to return a uniform (0, 1) variable xi by rescaling ci . The recurrence algorithm can be written in the following form: ci+1 = (aci + b)

mod C

xi+1 = ci+1 /C

(3.13) (3.14)

where the multiplier a and b are integer constants. The initial value (seed ) c0 is an integer in the given interval [0, C − 1] and is specified by the user, to be able to rerun a simulation with the same settings (’Reproducibility’). Once cN = c0 holds, the sequence starts over and will repeat itself, so the period length of a linear congruential method is an critical issue6 . Beside the linear congruential method other methods for uniform random number generators exist as well, e.g. mid-square method or the Mersenne twister7 . There exists a series of reliable methods to produce uniform pseudo-random sequences from Press et al. (1992), called ran0 to ran4, that have been well tested and documented. In our study we use ran2 to generate pseudo-random numbers. In Figure 3.1, we show pairwise projections of the first 2048 random points generated with ran2 on several two-dimensional unit intervals. The particular dimensions shown were selected randomly. It can be seen that clusters of points are possible, since the points are independently drawn and therefore have a certain chance to be very close to other points. On the other side, there are also empty spaces, where no points are observed. Moreover, it can be seen that the patterns are independent of the considered pair of dimension. The clumps and empty spaces are one reason for the relatively poor convergence rate of O(M −1/2 ) for Monte Carlo simulations with pseudo-random numbers. However, since the patters are independent of the dimension, the error bound and convergence rate is also independent of the dimension, which makes Monte Carlo simulation valuable for high-dimensional problems. 5

Niederreiter (1992) provides an overview of pseudo-random number generation and statistical tests for randomness of uniform random numbers as well as the most common methods to transform uniform random numbers into nonuniform random numbers. 6 If N = (C − 1), every distinct value in the interval is produced before repeating and the method is said to have full period. However, a large value of C does not guarantee this property, since the period length depends also on the choices for a and b. 7 See Glasserman (2003), J¨ ackel (2002) or Niederreiter (1992) for different methods of uniform random number generation.

12

3.1.1 Dimension 1 vs Dimension 2

3.1.2 Dimension 27 vs Dimension 32

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

3.1.3 Dimension 111 vs Dimension 128

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0.9

1.0

3.1.4 Dimension 255 vs Dimension 256

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0.0

0.1

0.2

Figure 3.1:

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Two-dimensional projections of pseudo-random sequences for different dimensions for the first 2048

points.

3.4

Quasi-Random Numbers

Quasi-random numbers8 , also known as low-discrepancy sequences9 , differ from ordinary Monte Carlo with pseudo-random numbers in the sense that they make no attempt to mimic randomness. Instead, they are designed to increase accuracy specifically by generating points which avoid gaps and clusters, i.e. they are too evenly distributed to be random. Therefore, Monte Carlo Simulation with quasi8

The term quasi-random is in some sense misleading, since low discrepancy sequences are totally deterministic. 9 A good discussion of low-discrepancy sequences can be found in Glasserman (2003) and Niederreiter (1992).

13

random numbers, henceforth we will call this Quasi-Monte Carlo (QMC), can improve the performance of Monte Carlo simulations, offering shorter computational times and/or higher accuracy, respectively. In optimal cases the convergence rate of QMC is O(log(M )d M −1 ), where d represents the dimension of the integration problem. For the remainder of this section, we give a short introduction to the numbertheoretical concept of QMC and introduce some common low-discrepancy sequences, which we use in our analysis.

3.4.1

Low Discrepancy and Integration Error

The main objective of quasi-random sequences is to improve the uniformity of the sequence, which is measured in terms of its discrepancy. This is defined by considering Q the number of points in rectangular subsets R = dj=1 (uj , vj ) of the d-dimensional d unit cube I d = [0, 1]d . Let XM = (xi )M i=1 be a sequence of M points in I , then for an arbitrary subset RR, the counting function ](XM , R), that indicates the number of points for which xi R, is defined as ](XM , R) =

M X

χR (xi ),

(3.15)

i=1

where χR is the characteristic function of set R. The (extreme) discrepancy of XM is then given by DM

](XM , R) = sup − m(λd (R) , M RR

(3.16)

where λd (R) denotes the d-dimensional volume (measure) of R. If R? is the set of all Q rectangles in I d of the form dj=1 (0, vj ), the star discrepancy is defined as ](XM , R) ? DM = sup − λd (R) . (3.17) M RR? ? If a sequence XM is perfectly uniformly distributed, then DM = DM = 0. For 1 ? the one-dimensional case, Niederreiter shows that DM ≥ 2M , whereas for the multi-

dimensional case, i.e. d ≥ 2, Niederreiter states that ”it is widely believed” that any point set XM satisfy

log(M )d−1 M the first M elements satisfies

? DM ≥ cd

and for any infinite sequence X∞

? DM ≥ c0d

log(M )d , M

14

(3.18)

(3.19)

where the constants cd > 0 and c0d > 0 depend only on d. Niederreiter [28] gives also some other definitions of discrepancy as well as further useful properties. The star discrepancy plays a central role for the definition of the integration error as defined in (3.6), when quasi-random numbers are used. The key result for QMC error bounds is known as the Koksma-Hlawka inequality10 , which gives an upper bound of the integration error by the product of two quantities – one depending only on the integration f , and one depending only on the quasi-random sequence used. The Koskma-Hlawka inequality says that if f has bounded variation V (f ) in the sense of Hardy and Krause on I d , then for any sequence XM I d , the following inequality holds11

Z M 1 X ? f (u)du ≤ V (f )DM . (3.20) f (xi ) − M Id i=1 R1 The variation for the one-dimensional case is V (f ) = 0 |df |. For the multidimensional case the variation in the sense of Hardy and Krause is given as V (f ) =

d X

X

V (k) (f ; i1 , ..., ik )

(3.21)

∂ df ∂u1 · · · ud du1 · · · dud

(3.22)

k=1 1