A central limit theorem for Latin hypercube sampling with dependence and application to exotic basket option pricing
aa r X i v : . [ q -f i n . C P ] N ov A central limit theorem for Latin hypercube sampling withdependence and application to exotic basket option pricing
Christoph Aistleitner ∗ Markus Hofer † Robert Tichy ‡ Abstract
We consider the problem of estimating E [ f ( U , . . . , U d )] , where ( U , . . . , U d ) denotes a randomvector with uniformly distributed marginals. In general, Latin hypercube sampling (LHS) is a powerfultool for solving this kind of high-dimensional numerical integration problem. In the case of dependentcomponents of the random vector ( U , . . . , U d ) one can achieve more accurate results by using Latinhypercube sampling with dependence (LHSD). We state a central limit theorem for the d -dimensionalLHSD estimator, by this means generalising a result of Packham and Schmidt. Furthermore we giveconditions on the function f and the distribution of ( U , . . . , U d ) under which a reduction of variancecan be achieved. Finally we compare the effectiveness of Monte Carlo and LHSD estimators numeri-cally in exotic basket option pricing problems. In this article we consider the problem of reducing the variance of a Monte Carlo (MC) estimator for spe-cial functionals of a random vector with dependent components. Several different techniques can be usedfor this kind of problem, with different advantages and shortcomings (for a detailed comparison, see [ ? ,Section 4]). A well-known technique is Latin hypercube sampling (LHS), which is a multi-dimensionalversion of the stratified sampling method and has been introduced by [ ? ]. Although this method is wellapplicable to many different types of problems, it cannot deal with dependence structures among the com-ponents of random vectors. Therefore, we consider Latin hypercube sampling with dependence (LHSD),which was introduced by [ ? ] and provides variance reduction for many problems, especially in financialmathematics.Consider the problem of estimating E [ f ( U , . . . , U d )] for a Borel-measurable and C -integrable function f : [0 , d → R , where ( U , . . . , U d ) is a random vector with uniformly distributed marginals and copula C . Let ( U i , . . . , U di ) , ≤ i ≤ n, denote an i.i.d. sample from this distribution. The standard MonteCarlo estimator, which is given by /n P ni =1 f ( U i , . . . , U di ) , is strongly consistent, and by the centrallimit theorem for sums of independent random variables the distribution of the scaled estimator converges ∗ Graz University of Technology, Institute of Mathematics A, Steyrergasse 30, 8010 Graz, Austria. e-mail: [email protected] . Research supported by the Austrian Research Foundation (FWF), Project S9603-N23. † Graz University of Technology, Institute of Mathematics A, Steyrergasse 30, 8010 Graz, Austria. e-mail: [email protected] . Research supported by the Austrian Research Foundation (FWF) and by the DOC [Doc-toral Fellowship Programme of the Austrian Academy of Sciences], Project S9603-N23. ‡ Graz University of Technology, Institute of Mathematics A, Steyrergasse 30, 8010 Graz, Austria. e-mail: [email protected] . Research supported by the Austrian Research Foundation (FWF), Project S9603-N23.
Mathematics Subject Classification:
Keywords:
Monte Carlo, Variance reduction techniques, Latin hypercube sampling, option pricing, variance gamma, proba-bilistic methods
1o a normal distribution, ie: √ n n X i =1 [ f ( U i , . . . , U di ) − E [ f ( U , . . . , U d )]] D −→ N (0 , σ MC ) , where σ MC = Var( f ( U , . . . , U d )) . In particular this means that the standard deviation of the estimatorconverges to zero with rate √ n .The aim of this paper is to establish a similar result for the LHSD estimator, under some additionalconditions on the copula C and the function f . This has already been done in the bivariate case by [ ? ] byusing a result of [ ? ]. ? , Proposition 5.9 also showed that under more restrictive conditions on the copulafunction C , the variance of the bivariate LHSD estimator does not exceed the variance of the standardMonte Carlo estimator.An important application of Monte Carlo integration techniques lies in the field of financial mathematics.Many problems in finance result in the numerical computation of high-dimensional integrals, for whichMC methods provide an efficient solution. Two examples are the pricing of Asian and discrete lookbackoptions on several possibly correlated assets. We will investigate these special derivatives in numericalexamples in the last section.This paper is organised as follows: in the second section we introduce the main ideas of LHSD and recallsome important results. Our main results are presented in the third section, where we state a central limittheorem and show under which conditions a reduction of variance, compared to the standard Monte Carlomethod, is possible. The last section is dedicated to a comparison of the effectiveness of LHSD and MCin numerical examples. In this section, we recall the concept of stratified sampling and its extensions to Latin hypercube samplingand Latin hypercube sampling with dependence. We also state a consistency result, which was proved by[ ? ]. Suppose that we want to estimate E ( f ( U )) , where U is an uniformly distributed random variable on theinterval [0 , (from now on denoted by U ([0 , ), and where f : [0 , → R is a Borel-measurable andintegrable function. By the simple fact that E ( f ( U )) = n X i =1 E ( f ( U ) | U ∈ A i ) P ( U ∈ A i ) , where the intervals A , . . . , A n (the so-called strata ) form a partition of [0 , , we get an estimator for E ( f ( U )) by sampling U conditionally on the events { U ∈ A i } , i = 1 , . . . , n . Choosing strata of the form A i = [ i − n , in ) we can simply transform independent samples U , . . . , U n from U ([0 , by setting V i := i − n + U i n , i = 1 , . . . , n, which implies V i ∈ A i , i = 1 , . . . , n . The resulting estimator for E ( f ( U )) given by n P ni =1 f ( V i ) isconsistent, and by the central limit theorem for sums of independent random variables the limit varianceis smaller than the limit variance of a standard Monte Carlo estimator. For a more detailed analysis ofstratified sampling techniques, see [ ? , Section 4.3.1].2his approach can be extended to the multivariate case in different ways. If we require that there has to beexactly one sample in every stratum, we need to draw n d samples, which is not feasible for a high numberof dimensions d . One way to avoid this problem is Latin hypercube sampling. Assume we want to estimate E ( f ( U , . . . , U d )) , where f : [0 , d → R is a Borel-measurable and integrable function. For fixed n wegenerate n independent samples denoted by ( U i , . . . , U di ) , i = 1 , . . . , n , where the U ji , j = 1 , . . . , d are uniformly distributed on [0 , . Additionally, we generate d independent permutations of { , . . . , n } ,denoted by π , . . . , π d , drawn from a discrete uniform distribution on the set of all possible permutations.Denote by π ji the value to which i is mapped by the j -th permutation. Then the j -th component of a Latinhypercube sample is given by V ji := π ji − n + U ji n , j = 1 , . . . , d ; i = 1 , . . . , n. By fixing a dimension j , the components ( V j , . . . , V jn ) form a stratified sample with strata of equallength. It can be shown that the resulting estimator for E ( f ( U )) is consistent, and by assuming that f ( U , . . . , U d ) has a finite second moment it follows that the variance of the LHS estimator n n X i =1 f ( V i , . . . , V di ) is smaller than the variance of the standard MC estimator, provided the number of sample points is suffi-ciently large, see [ ? ]. If f is bounded a central limit theorem for the LHS estimator can be shown, see [ ? ].Berry-Esseen-type bounds are also known, see [ ? ]. A detailed discussion of LHS is given in [ ? , Section4.4].This technique is not suitable for dealing with random vectors with dependent components since the ran-dom variables V ji , j = 1 , . . . , d , are independent. One way to extend the LHS method to random vectorswith dependent components is to apply LHS to independent components and then introduce dependen-cies through a transformation of the LHS points. Such a procedure is tedious in general, and we will notpursue this approach any further. In this subsection, we introduce Latin hypercube sampling with dependence. The main difference tothe LHS method is that instead of random permutations π i we use rank statistics, which are defined asfollows: Definition 2.1 (Rank statistics)
Let X , . . . , X n be i.i.d. random variables with a continuous distribu-tion function. Denote the ordered random variables by X (1) < · · · < X ( n ) , P -a.s. We call the index of X i within X (1) < · · · < X ( n ) the i -th rank statistic, given by r i,n = r i,n ( X , . . . , X n ) := n X k =1 { X k ≤ X i } . (1)Consider a random vector U = ( U , . . . , U d ) , where every component U j is uniformly distributed on [0 , and the dependence structure of U is modeled by a copula C . Let ( U i , . . . , U di ) , i = 1 , . . . , n denotea sequence of independent samples of ( U , . . . , U d ) , and let r ji,n be the i -th rank statistic of ( U j , . . . , U jn ) for i = 1 , . . . , n and j = 1 , . . . , d . Then a LHSD is given by V ji,n := r ji,n − n + η ji,n n , i = 1 , . . . , n, ∀ j = 1 , . . . , d, (2)3here η ji,n are random variables in [0 , . It is clear that ( V j ,n , . . . , V jn,n ) forms a stratified sampling inevery dimension j , where every stratum has equal length. ? consider different choices for η ji,n to obtain special properties. For example, by choosing all η ji,n uni-formly distributed on [0 , and independent of U ji , the distribution of the V ji,n within their strata is uni-form. This choice has the disadvantage of necessitating the generation of n random variables instead ofonly n . An effective choice in terms of computation time is η ji,n = 1 / , which means that every V ji,n islocated exactly in the centre of its stratum. In the remainder of this section, we briefly recall a result of[ ? ] concerning the consistency of the LHSD estimator for E ( f ( U )) , which is defined by n n X i =1 f ( V i,n , . . . , V di,n ) . (3)The usual law of large numbers for sums of independent random variables does not apply in this casefor two reasons: firstly in each dimension the samples fail to be independent because of the applicationof the rank statistic, and secondly, increasing the samples size n by one changes every term of the suminstead of just adding one. Nevertheless, it can be shown that the following consistency result holds, see[ ? , Proposition 4.1]: Proposition 2.1
Let f : [0 , d → R be bounded and continuous C-a.e. . Then the LHSD estimator (3) isstrongly consistent, ie : n n X i =1 f ( V i,n , . . . , V di,n ) P a.s. −−−→ E ( f ( U , . . . , U d )) , as n → ∞ . In this section we investigate the speed of convergence of the LHSD estimator and discuss situations inwhich the use of LHSD results in a reduction of variance. This has already been done for the bivariatecase by [ ? ]. They have also guessed the higher-dimensional version of the main theorem, but no rigorousproof was given. Because of the fact that most problems in finance for which Monte Carlo techniquesare suitable are high-dimensional integration problems, it is reasonable to investigate the speed of con-vergence and the (asymptotic) value of the variance also in the multivariate case.In the sequel, let C n denote the empirical distribution of the LHSD sample given by C n ( u , . . . , u d ) := 1 n n X i =1 { V i,n ≤ u ,...,V di,n ≤ u d } , which is a distribution function. Furthermore, we define C n as C n ( u , . . . , u d ) := 1 n n X i =1 { F n ( U i ) ≤ u ,...,F dn ( U di ) ≤ u d } , (4)where F jn ( u ) = 1 n n X i =1 { U ji ≤ u } , u ∈ [0 , , are the one-dimensional empirical distribution functions based on U j , . . . , U jn for j = 1 , . . . , d . To for-mulate a central limit theorem we will need some regularity conditions on the integrand f and the copula C . 4 efinition 3.1 (Hardy-Krause bounded variation) A function f : [0 , d → R is of bounded variation(in the sense of Hardy-Krause) if V ( f ) < ∞ with V ( f ) = d X k =1 X ≤ i <...
A function f : [0 , d → R is right continuous if for any sequence ( u n , u n , . . . , u dn ) n ∈ N with u jn ↓ u j , j = 1 , . . . , d, lim n →∞ f ( u n , u n , . . . , u dn ) = f ( u , u , . . . , u d ) . The next statement concerning the convergence of random sequences will be used to prove Proposition3.1 and Theorem 3.2. For more details see eg [ ? , Theorem 18.8]. Lemma 3.1
Let ( X n ) n ≥ and ( Y n ) n ≥ be sequences of R -valued random variables, with X n D −→ X and | X n − Y n | P −→ . Then Y n D −→ X . The following proposition of [ ? ] is a generalization of earlier results of [ ? ] and [ ? ]. It is the essentialingredient in proofs of our main theorems. Proposition 3.1
Assume that C is differentiable with continuous partial derivatives ∂ j C ( u , . . . , u d ) = ∂C ( u ,...,u d ) ∂u j for j = 1 , . . . , d . Then √ n (cid:16) e C n ( u , . . . , u d ) − C ( u , . . . , u d ) (cid:17) D −→ G C ( u , . . . , u d ) , where e C n ( u , . . . , u d ) = 1 n n X k =1 { U k ≤ F − n ( u ) ,...,U dk ≤ F d − n ( u d ) } , denotes the empirical copula function and F j − n denote the generalised quantile functions of F jn for j =1 , . . . , d , defined by F j − n ( u ) = inf { x ∈ R | F jn ( x ) ≥ u } . Furthermore, G C is a centred Gaussian random field given by G C ( u , . . . , u d ) = B C ( u , . . . , u d ) − d X j =1 ∂ j C ( u , . . . , u d ) B C (1 , . . . , , u j , , . . . , , (5)5 C is a d-dimensional pinned Brownian sheet on [0 , d with covariance function E [ B C ( u , . . . , u d ) · B C ( u , . . . , u d )] = C (( u , . . . , u d ) ∧ ( u , . . . , u d )) − C ( u , . . . , u d ) C ( u , . . . , u d ) , (6) where ( u , . . . , u d ) ∧ ( u , . . . , u d ) denotes the componentwise minimum. We can formulate a similar result for the sequence C n . Proposition 3.2
Under the conditions of Proposition 3.1, √ n (cid:16) C n ( u , . . . , u d ) − C ( u , . . . , u d ) (cid:17) D −→ G C ( u , . . . , u d ) (7) holds, where all definitions are as in Proposition 3.1 and C n ( u , . . . , u d ) is given in (4) . Proof:
We only have to show that the supremum of the difference of C n and e C n vanishes for n → ∞ to applyLemma 3.1, which completes the proof. Note that C n and e C n coincide on the grid { ( i /n, . . . , i d /n ) , ≤ i , . . . , i d ≤ n } . It follows that sup u ,...,u d | e C n ( u , . . . , u d ) − C n ( u , . . . , u d ) |≤ max ≤ i ,...,i d ≤ n (cid:12)(cid:12)(cid:12) e C n (cid:16) i n , . . . , i d n (cid:17) − e C n (cid:16) i − n , . . . , i d − n (cid:17)(cid:12)(cid:12)(cid:12) ≤ dn . Thus, sup u ,...,u d | e C n ( u , . . . , u d ) − C n ( u , . . . , u d ) | → for n → ∞ and (7) follows. (cid:3) In the sequel, all U i , i = 1 , . . . , d are uniformly distributed random variables on [0 , and all integralshave to be understood in the sense of Lebesgue-Stieltjes. Note that the next theorem is an extension of [ ? ,Theorem 6] from the case of bivariate to the case of multi-variate random vectors U = ( U , . . . , U d ) . Theorem 3.1
Let the copula C of ( U , . . . , U d ) have continuous partial derivatives and let f : [0 , d → R be a right-continuous function of bounded variation in the sense of Hardy-Krause. Then √ n n X i =1 (cid:16) f ( F n ( U i ) , . . . , F dn ( U di )) − E [ f ( U , . . . , U d )] (cid:17) D −→ Z [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) , where the function b f : [0 , d → R is defined by: b f ( u , . . . , u d ) = (cid:26) if at least one u j = 1 , for j = 1 , . . . , d,f ( u , . . . , u d ) otherwise. (8) Furthermore, the limit distribution is Gaussian.
Proof:
By definition b f is right-continuous and of bounded variation in the sense of Hardy-Krause. Furthermore,it follows that almost surely √ n n X i =1 (cid:16) f ( F n ( U i ) , . . . , F dn ( U di )) − E [ f ( U , . . . , U d )] (cid:17) = 1 √ n n X i =1 (cid:16) b f ( F n ( U i ) , . . . , F dn ( U di )) − E [ b f ( U , . . . , U d )] (cid:17) ,
6y the fact that C is continuous on [0 , d .We use a multidimensional integration-by-parts technique proposed by [ ? , Proposition 2]. Using the no-tation of [ ? ] we get √ n n X i =1 (cid:16) b f ( F n ( U i ) , . . . , F dn ( U di )) − E [ b f ( U , . . . , U d )] (cid:17) = √ n Z [0 , d b f ( u , . . . , u d ) d ( C n − C )( u , . . . , u d )= √ n d X k =0 ( − k X ,...,d ; k ∆ ∗ j k +1 ,...,j d Z [0 , k ( C n − C )( u , . . . , u d ) d j ,...,j k b f ( u , . . . , u d ) . (9)Here P ,...,d ; k denotes the sum over all possible partitions of the set { j , . . . , j d } into two subsets { j , . . . , j k } and { j k +1 , . . . , j d } of k respectively d − k elements, where each partition is taken exactlyonce. In the cases k = 0 and k = d , the sum is interpreted as being reduced to one term.Furthermore, the operator d j ,...,j k indicates that the integral only applies to the variables j , . . . , j k . Notethat after the application of the integral with respect to d j ,...,j k b f ( u , . . . , u d ) , the integrated function isa function in d − k variables. Furthermore for a function g of d − k variables, the operator ∆ ∗ j k +1 ,...,j d isgiven by ∆ ∗ j k +1 ,...,j d g ( j k +1 , . . . , j d ) = X { i ,...,i d − k }∈{ , } d − k ( − m g ( i , . . . , i d − k ) , where m denotes the number of zeros in { i , . . . , i d − k } . This means that, for j / ∈ { j , . . . , j k } ∆ ∗ j Z [0 , d − k ( C n − C )( u , . . . , u d ) d j ,...,j k b f ( u , . . . , u d )= Z [0 , d − k ( C n − C )( u , . . . , u j − , , u j +1 , . . . , u d ) d j ,...,j k b f ( u , . . . , u j − , , u j +1 , . . . , u d ) − Z [0 , d − k ( C n − C )( u , . . . , u j − , , u j +1 , . . . , u d ) d j ,...,j k b f ( u , . . . , u j − , , u j +1 , . . . , u d ) and ∆ ∗ j k +1 ,...,j d = ∆ ∗ j k +1 . . . ∆ ∗ j d . Thus √ n d X k =0 ( − k X ,...,d ; k ∆ ∗ j k +1 ,...,j d Z [0 , k ( C n − C )( u , . . . , u d ) d j ,...,j k b f ( u , . . . , u d )= √ n d − X k =0 ( − k X ,...,d ; k ∆ ∗ j k +1 ,...,j d Z [0 , k ( C n − C )( u , . . . , u d ) d j ,...,j k b f ( u , . . . , u d )+ √ n ( − d Z [0 , d ( C n − C )( u , . . . , u d ) d b f ( u , . . . , u d )= √ n ( − d Z [0 , d ( C n − C )( u , . . . , u d ) d b f ( u , . . . , u d ) . The term √ n d − X k =0 ( − k X ,...,d ; k ∆ ∗ j k +1 ,...,j d Z [0 , k ( C n − C )( u , . . . , u d ) d j ,...,j k b f ( u , . . . , u d ) u j , j = 1 , . . . , d is equal to one and therefore b f ( u , . . . , u d ) = 0 by definition, or, secondly,at least one u j , j = 1 , . . . , d is equal to zero, hence C n ( u , . . . , u d ) = C ( u , . . . , u d ) = 0 .Thus, by the continuous mapping theorem and (7), it follows that √ n n X i =1 (cid:16) f ( F n ( U i ) , . . . , F dn ( U di )) − E [ f ( U , . . . , U d )] (cid:17) = ( − d √ n Z [0 , d ( C n − C )( u , . . . , u d ) d b f ( u , . . . , u d ) D −→ Z [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) . Since R [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) is a continuous, linear transformation of a tight Gaussianprocess, it follows that the limiting distribution is Gaussian. (cid:3) Remark 3.1
The reason for using the function b f instead of f is that the integrals of dimension k =2 , . . . , d − in (9) are in general not vanishing. The one-dimensional integrals are zero for every right-continuous function of bounded variation f because of special properties of the function C n , for moredetails see [ ? ]. In particular, this means that in the two-dimensional case it is sufficient to assume b f ( x ) = f ( x ) , x ∈ R . With this assumption instead of (8) and d = 2 , Theorem 3.1 is equivalent to [ ? , Theorem 6]. We use thefunction b f to get a more convenient representation for the limit variance of the LHSD technique, whichwe state in the next theorem. Theorem 3.2
Under the assumptions and notations of Theorem 3.1, we have √ n n X i =1 (cid:16) f ( V i,n , . . . , V di,n ) − E [ f ( U , . . . , U d )] (cid:17) D −→ N (0 , σ LHSD ) , (10) where σ LHSD = Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i d b f ( u , . . . , u d ) d b f ( u , . . . , u d ) . (11) Proof:
We want to apply Theorem 3.1 together with Lemma 3.1, so we have to show that √ n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X i =1 h f ( V i,n , . . . , V di,n ) − f ( F n ( U i ) , . . . , F dn ( U di )) i(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) → , as n → ∞ . By [ ? , Corollary 1] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X i =1 h f ( V i,n , . . . , V di,n ) − f ( F n ( U i ) , . . . , F dn ( U di )) i(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ V ( f ) < ∞ , V ( f ) is the Hardy-Krause variation of f . Hence √ n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X i =1 h f ( V i,n , . . . , V di,n ) − f ( F n ( U i ) , . . . , F dn ( U di )) i(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) → , as n → ∞ , which, together with Lemma 3.1 and Theorem 3.1, proves equation (10).To derive equation (11) we apply Fubini’s theorem to E [( R [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d )) ] . By [ ? ,Theorem 3] a function of bounded variation b f can always be written as the difference of two completelymonotone functions g, h and therefore an integral with respect to b f can be written as a difference of twointegrals with respect to positive measures g, h . Thus E h(cid:16) Z [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) (cid:17) i == E h(cid:16)Z [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) (cid:17) · (cid:16)Z [0 , d G C ( u , . . . , u d ) d b f ( u , . . . , u d ) (cid:17)i = E h(cid:16)Z [0 , d G C ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d G C ( u , . . . , u d ) dh ( u , . . . , u d ) (cid:17) · (cid:16)Z [0 , d G C ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d G C ( u , . . . , u d ) dh ( u , . . . , u d ) (cid:17)i = E h(cid:16)Z [0 , d G C ( u , . . . , u d ) G C ( u , . . . , u d ) dg ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d G C ( u , . . . , u d ) G C ( u , . . . , u d ) dh ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d G C ( u , . . . , u d ) G C ( u , . . . , u d ) dg ( u , . . . , u d ) dh ( u , . . . , u d )+ Z [0 , d G C ( u , . . . , u d ) G C ( u , . . . , u d ) dh ( u , . . . , u d ) dh ( u , . . . , u d ) (cid:17)i = Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i dg ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i dh ( u , . . . , u d ) dg ( u , . . . , u d ) − Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i dg ( u , . . . , u d ) dh ( u , . . . , u d )+ Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i dh ( u , . . . , u d ) dh ( u , . . . , u d )= Z [0 , d E h G C ( u , . . . , u d ) G C ( u , . . . , u d ) i d b f ( u , . . . , u d ) d b f ( u , . . . , u d ) , where the use of Fubini’s theorem is justified since b f is bounded and E [ XY ] < ∞ for two jointly normalrandom variables X and Y . (cid:3) emark 3.2 Note that by (5) and (6) the expression for σ LHSD in equation (11) can be represented interms of C . Additionally, further simplifications can be given for the following terms: E [ B C ( u , . . . , u d ) · B C (1 , . . . , , u j , , . . . , C (( u , . . . , u j − , u j ∧ u j , u j +1 . . . , u d )) − C ( u , . . . , u d ) u j , E [ B C (1 , . . . , , u i , , . . . , · B C (1 , . . . , , u j , , . . . , C ((1 , . . . , , u i , , . . . , , u j , , . . . , − u i u j , E [ B C (1 , . . . , , u j , , . . . , · B C (1 , . . . , , u j , , . . . , u j ∧ u j − u j u j , since C (1 , . . . , , u j , , . . . ,
1) = u j for all j = 1 , . . . , d . It is important to know if the LHSD estimator has a smaller variance than the Monte Carlo estimator. Thevariance of a standard Monte Carlo estimator is given by σ MC = Z [0 , d f ( u , . . . , u d ) dC ( u , . . . , u d ) − (cid:16)Z [0 , d f ( u , . . . , u d ) dC ( u , . . . , u d ) (cid:17) . We use this fact to establish a relation between σ MC and σ LHSD . Proposition 3.3
Let the copula C of ( U , . . . , U d ) have continuous partial derivatives, let f : [0 , d → R be a right-continuous function of bounded variation in the sense of Hardy-Krause and let b f be asdefined in Theorem 3.1. Set ∂ j C ( u , . . . , u d ) = ∂C ( u ,...,u d ) ∂u j and C i,j ( u i , u j ) = (cid:26) C (1 , . . . , , u i , , . . . , , u j , , . . . , , i = ju i ∧ u j , i = j. Then σ LHSD = σ MC + Z [0 , d d X j =1 ∂ j C ( u , . . . , u d ) (cid:16) C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) (cid:17) + d X j =1 d X i =1 ∂ j C ( u , . . . , u d ) ∂ i C ( u , . . . , u d ) (cid:16) C i,j ( u i , u j ) − u i u j (cid:17) d b f ( u , . . . , u d ) d b f ( u , . . . , u d ) . (12) Proof:
Note that Z [0 , d f ( u , . . . , u d ) dC ( u , . . . , u d ) = Z [0 , d f ( u , . . . , u d ) f ( u , . . . , u d ) dC ( u ∧ u , . . . , u d ∧ u d ) , and that the function C ( u ∧ u , . . . , u d ∧ u d ) is also a copula, which follows by observing that C ( u ∧ u , . . . , u d ∧ u d ) = P ( U ≤ u ∧ u , . . . , U d ≤ u d ∧ u d )= P ( U ≤ u , U ≤ u , . . . , U d ≤ u d , U d ≤ u d ) is a joint probability distribution with uniform marginals.By integration-by-parts like in Theorem 3.1 it follows for the variance of the Monte Carlo estimator that σ MC = Z [0 , d f ( u , . . . , u d ) dC ( u , . . . , u d ) − (cid:16)Z [0 , d f ( u , . . . , u d ) dC ( u , . . . , u d ) (cid:17) Z [0 , d f ( u , . . . , u d ) f ( u , . . . , u d ) dC (cid:16) ( u , . . . , u d ) ∧ ( u , . . . , u d ) (cid:17) − Z [0 , d f ( u , . . . , u d ) f ( u , . . . , u d ) dC ( u , . . . , u d ) dC ( u , . . . , u d )= Z [0 , d C (cid:16) ( u , . . . , u d ) ∧ ( u , . . . , u d ) (cid:17) d b f ( u , . . . , u d ) d b f ( u , . . . , u d ) − Z [0 , d C ( u , . . . , u d ) C ( u , . . . , u d ) d b f ( u , . . . , u d ) d b f ( u , . . . , u d ) . The proof is completed by using equations (5), (6), (11) and Remark 3.2. (cid:3)
Theorem 3.3
Let C and f satisfy the assumptions in Theorem 3.1 and let b f be defined as in Theorem 3.1.Furthermore let the function f be monotone non-decreasing in each argument and max x ∈ [0 , d ( f ( x )) ≤ . Moreover assume that C satisfies the following conditions: C ( u , . . . , u d ) u j ≥ ∂ j C ( u , . . . , u d ) , j ∈ { , . . . , d } , (13) d X i =1 ,i = j C i,j ( u j , u i ) u j ≤ ( d − u j + C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) C ( u , . . . , u d ) , (14) where u j ∈ [0 , , ( u , . . . , u d ) , ( u , . . . , u d ) ∈ [0 , d .Then σ LHSD ≤ σ MC . Proof:
By the assumptions on f it follows that b f is right-continuous, of bounded variation in the sense of Hardy-Kraus and monotone non-decresing in each argument. Thus by (12) it is sufficient to show that d X j =1 ∂ j C ( u , . . . , u d ) (cid:16) C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) (cid:17) + d X j =1 d X i =1 ∂ i C ( u , . . . , u d ) ∂ j C ( u , . . . , u d ) (cid:16) C i,j ( u j , u i ) − u j u i (cid:17) ≤ for all ( u , . . . , u d ) , ( u , . . . , u d ) ∈ [0 , d .This is true if (cid:0) C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) (cid:1) ≤ d X i =1 ∂ i C ( u , . . . , u d ) (cid:0) u j u i − C i,j ( u j , u i ) (cid:1) holds for every j ∈ { , . . . , d } and all u j ∈ [0 , , ( u , . . . , u d ) ∈ [0 , d .First we show that C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) ≤ ∂ j C ( u , . . . , u d ) (cid:0) u j u j − u j ∧ u j (cid:1) . Note that this is always true if u j ∧ u j ∈ { , } . Now assume that < u j ≤ u j < , then C ( u , . . . , u d ) u j − C ( u , . . . , u d ) ≤ ∂ j C ( u , . . . , u d ) (cid:0) u j u j − u j (cid:1) ( u , . . . , u d )( u j − ≤ ∂ j C ( u , . . . , u d ) u j ( u j − C ( u , . . . , u d ) u j ≥ ∂ j C ( u , . . . , u d ) which is true by assumption (13). Next assume that < u j < u j < , then C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) ≤ ∂ j C ( u , . . . , u d ) (cid:0) u j u j − u j (cid:1) C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) ≤ ∂ j C ( u , . . . , u d ) u j (cid:0) u j − (cid:1) C ( u , . . . , u d ) − C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) u j ≤ ∂ j C ( u , . . . , u d ) (cid:0) u j − (cid:1) C ( u , . . . , u d ) − C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) u j ≤ C ( u , . . . , u d ) u j (cid:0) u j − (cid:1) C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) u j ≥ C ( u , . . . , u d ) u j , which holds since assumption (13) implies that C ( u ,...,u d ) u j is non-increasing in u j for all u j ∈ [0 , , ( u , . . . , u d ) ∈ [0 , d .Let C ( u , . . . , u d ) > then C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) ≤ d X i =1 i = j ∂ i C ( u , . . . , u d ) (cid:0) u j u i − C i,j ( u j , u i ) (cid:1) C ( u , . . . , u d ) u j − C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) ≤ d X i =1 i = j C ( u , . . . , u d ) u i (cid:0) u j u i − C i,j ( u j , u i ) (cid:1) ( d − u j + C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) C ( u , . . . , u d ) ≥ d X i =1 i = j C i,j ( u j , u i ) u i which is true by assumption (14). The case C ( u , . . . , u d ) = 0 follows by the fact that C ( u ,...,u d ) u i ≤ for all ( u , . . . , u d ) ∈ [0 , d . (cid:3) Remark 3.3
Note that in the two-dimensional case, assumption (13) is equivalent to the left tail increas-ing property which implies a positive quadrant dependence of the copula C . Loosely speaking this meansthat the components of C are more likely to be simultaneously small or simulatneously large than in theindependent case. More information on different dependence properties can be found in [ ? ] and [ ? ]. In the following two remarks we give examples of copula distributions which satisfy the assumptions ofTheorem 3.3.
Remark 3.4
Consider a multi-dimensional, one-parametric extension of the Farlie-Gumbel-Morgenstern(FGM) copula given by C ( u , . . . , u d ) = d Y i =1 u i ! α d Y i =1 (1 − u i ) + 1 ! here α ∈ [ − , . Simple calculations show that the assumption (13) is true if α ∈ [0 , . Now considerthe right hand-side of (14) d X i =1 ,i = j C i,j ( u j , u i ) u i = d X i =1 ,i = j u j u i u i = ( d − u j . Finally assumption (14) holds since C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) C ( u , . . . , u d )= min (cid:18) , C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) C ( u , . . . , u d ) (cid:19) = min , (cid:16)Q di =1 ,i = j u i (cid:17) u j (cid:16) α Q di =1 ,i = j (1 − u i )(1 − u j ) + 1 (cid:17)(cid:16)Q di =1 u i (cid:17) (cid:16) α Q di =1 (1 − u i ) + 1 (cid:17) = min , u j (cid:16) α Q di =1 ,i = j (1 − u i )(1 − u j ) + 1 (cid:17) u j (cid:16) α Q di =1 (1 − u i ) + 1 (cid:17) ≥ u j for α ∈ [0 , .Note that the independence copula C ( u , . . . , u d ) = Q di =1 u i is the special case of the FGM copula with α = 0 , therefore Theorem 3.3 holds also for the independence copula. Remark 3.5
A multi-dimension version of the Ali-Mikhail-Haq (AMH) copula is given by C ( u , . . . , u d ) = Q di =1 u i − α Q di =1 (1 − u i ) where α ∈ [ − , . As in the previous example it is easy to see that (13) is fullfilled if α ∈ [0 , .To prove (14) consider again the term on the right hand-side d X i =1 ,i = j C i,j ( u j , u i ) u i = d X i =1 ,i = j u j u i u i = ( d − u j . Furthermore Theorem 3.3 applies since C ( u , . . . , u j − , u j ∧ u j , u j +1 , . . . , u d ) C ( u , . . . , u d )= min (cid:18) , C ( u , . . . , u j − , u j , u j +1 , . . . , u d ) C ( u , . . . , u d ) (cid:19) = min , u j (cid:16)Q di =1 ,i = j u i (cid:17) (cid:16) − α Q di =1 (1 − u i ) (cid:17)(cid:16)Q di =1 u i (cid:17) (cid:16) − α Q di =1 (1 − u i )(1 − u j ) (cid:17) ≥ u j Application to option pricing
In this section we illustrate the effectiveness of Latin hypercube sampling with dependence in basketoption pricing problems. The derivatives which we consider are Asian and lookback basket options. Let ( S t ) t ≥ be a d -dimensional vector of asset price processes and let ( S jt ) t ≥ denote its j -th component.Then the price of an Asian basket call option is given byABC = E h e − rT (cid:16) m m X j =1 d d X i =1 S it j − K (cid:17) + i , where K > denotes the fixed strike price, d is the number of underlying assets, t < t < t <. . . < t m = T denote the observation points, T is the maturity of the option and r denotes the risk freeinterest rate. Similarly, the price of a discrete lookback basket call option is given byDLC = E h e − rT (cid:16) max j =1 ,...,m d d X i =1 S it j − K (cid:17) + i . As a model for the asset price process ( S jt ) t ≥ of each asset j = 1 , . . . , d , we use S jt = S j e ( w j − r ) t + X jt , j = 1 , . . . , d, t ≥ , where w j ∈ R are constants, S j > denote the constant initial asset values and X jt are variance gamma(VG) processes for j = 1 , . . . , d . The VG process ( X jt ) t ≥ with parameters ( θ j , σ j , c j ) , which was firstintroduced by [ ? ], is defined as a subordinated Brownian motion by X jt = X jt ( θ j , σ j , c j ) = B jG jt ( c j , ( θ j , σ j ) , j = 1 , . . . , d, t ≥ , (15)where B jt ( θ j , σ j ) are independent Brownian motions with drift parameters θ j and volatility parameters σ j , j = 1 , . . . , d, and G jt ( c j , are independent gamma processes independent of B j , j = 1 , . . . , d withdrift equal to one and volatility c j > . To ensure that the discounted value of a portfolio invested in theasset is a martingale, we choose w j = log(1 − µ j c j − ( σ j ) c j / /c j , j = 1 , . . . , d. By [ ? ] a VG process can also be represented as the difference of two independent gamma processes, ie X jt = G + ,jt − G − ,jt , j = 1 , . . . , d . Let ( µ j + , ν j + ) and ( µ j − , ν j − ) denote the parameters of the gamma pro-cesses G + ,j , G − ,j , respectively. These pairs of parameters can be easily calculated from the parametersin equation (15) through µ j ± = ( p ( θ j ) + 2( σ j ) /c j ± θ j ) / , ν j ± = ( µ j ± ) c j , j = 1 , . . . , d. Due to the fact that a gamma process has non-decreasing paths, G + ,jt corresponds to the positive move-ments of X jt and G − ,jt corresponds to the negative movements of X jt . Our assumption is that all pos-itive movements of components of X t = ( X t , . . . , X dt ) are dependent and all negative movements ofcomponents of X t are dependent, but positive (negative) movements of the j -th component are inde-pendent of negative (positive) movements of all other components, for all j = 1 , . . . , d . The depen-dence structure between positive and negative movements will be modelled by copulae C ± , respectively.Summarising, the increment of the d -dimensional gamma processes in the interval [ t i − , t i ] given by ( G ± , t i − G ± , t i − , . . . , G ± ,dt i − G ± ,dt i − ) has cumulative distribution function C ± ( F − , ± , . . . , F − d, ± ) , where F − j, ± is the inverse cumulative distribution function of a gamma distribution with the specific parametersof the j -th asset. 14 .1 Numerical results In this subsection, we compare the performance of LHSD with a standard Monte Carlo method in optionpricing problems.
Parameters of the numerical examplesVG parameters: µ j , j = 1 , . . . , d -0.2859 σ j , j = 1 , . . . , d c j , j = 1 , . . . , d Option parameters: number of assets d T S j , j = 1 , . . . , d r k t i − t i − , i = 1 , . . . , k Simulation parameters: number of simulated option prices per estimator n m η ji,n , j = 1 , . . . , d, i = 1 , . . . , n ( S t ) t ≥ . The parameter values are taken from a calibration of the VG process against options on theS&P 500 index by [ ? ]. We observed in price valuations, which we do not state here in detail, that thecomputation of one LHSD estimator took about . times of the computation time of a correspondingMonte Carlo estimator. Nevertheless in our concrete implementation the most time consuming part wasthe transformation of uniformly distributed random variables into gamma distributed random variables.This has to be done only once for all LHSD estimations since by (2) where η ji,n = 1 / , j = 1 , . . . , d, i =1 , . . . , n one only needs fixed quantiles of the gamma distribution. Therefore computation of 4000 LHSDestimators was about five times faster than the computation of 4000 Monte Carlo estimators. One the otherhand for the Monte Carlo estimator, one has to perform the transformation dn times for each estimator.Using the parameters of Table 1, the evaluation of each of the option values included the computationof an -dimensional integral. Standard deviation and variance were computed based on the m = 100 runs of the LHSD and MC estimators. The ratios in columns 6 and 7 of each table were computed as thequotient of MC value and LHSD value.It is obvious that the effectiveness of LHSD compared to MC decreases with increasing strike price K .The same phenomenon was also observed by [ ? ] in a multi-dimensional Black-Scholes model for theLHSD estimator and by [ ? ] for the standard LHS estimator.15 rices of Asian basket call options with varying strike price Kα K
Price LHSD Price MC Std. Dev. LHSD Std. Dev. MC Std. Dev. ratio Var. ratio0.5 80 22.0542 22.0448 0.00071 0.00748 10.419 108.5750.5 90 12.5511 12.5419 0.00080 0.00748 9.270 85.9440.5 100 3.79294 3.78732 0.00241 0.00621 2.577 6.6420.5 110 0.17227 0.17210 0.00119 0.00140 1.174 1.3790.5 120 0.00024 0.00024 0.000040 0.000041 1.009 1.018Table 2: Prices of Asian basket call options, where the dependence structure of positive and negativemovements are modelled by a FGM copula with parameter α . Prices of Lookback basket call options with varying strike price
Kα K
Price LHSD Price MC Std. Dev. LHSD Std. Dev. MC Std. Dev. ratio Var. ratio0.5 80 25.662 25.658 0.00294 0.00839 2.850 8.1250.5 90 16.151 16.147 0.00294 0.00839 2.850 8.1250.5 100 6.893 6.890 0.00322 0.00760 2.356 5.5530.5 110 1.192 1.192 0.00305 0.00406 1.332 1.7750.5 120 0.060 0.060 0.00086 0.00089 1.029 1.060Table 3: Prices of Lookback basket call options, where the dependence structure of positive and negativemovements are modelled by a FGM copula with parameter αα