We assume that the stock price of the underlaying assetS =S(t) is driven by a L´evy processL(t)

A MIXED MONTE CARLO AND QUASI-MONTE CARLO METHOD WITH APPLICATIONS TO MATHEMATICAL FINANCE

ALIN V. ROS¸CA

Abstract. In this paper, we apply a mixed Monte Carlo and Quasi-Monte Carlo method, which we proposed in a previous paper, to problems from mathematical finance. We estimate the value of an European Call option and of an Asian option using our mixed method, under different horizont times. We assume that the stock price of the underlaying assetS =S(t) is driven by a L´evy processL(t). We compare our estimates with the esti- mates obtained by using the Monte Carlo and Quasi-Monte Carlo methods.

Numerical results show that a considerable improvement can be achieved by using the mixed method.

1. Introduction

The valuation of financial derivatives is one of the most important problems from mathematical finance. The risk-neutral price of such a derivative can be ex- pressed in terms of a risk-neutral expectation of a random payoff. In some cases, the expectation is explicitly computable, such as the Black & Scholes formula for pricing call options on assets modelled by a geometric Brownian motion. However, if we consider an Asian option, there exists no longer closed form expressions for the price, and therefore numerical methods are involved. This is the case, even if we consider call options written on assets with non-normal returns. Among these methods, Monte Carlo (MC) and Quasi-Monte Carlo (QMC) methods play an increasingly important role.

Received by the editors: 10.09.2007.

2000Mathematics Subject Classification.91B24, 91B28, 65C05, 11K45, 11K36, 62P05.

Key words and phrases. Monte Carlo integration, Quasi-Monte Carlo integration, normal inverse Gaussian distribution,H-discrepancy,H-distributed low-discrepancy sequences,H-mixed sequences, option pricing.

One of the first applications of the MC method in this field appeared in Boyle [2], who used simulation to estimate the value of a standard European option.

Applications of the QMC method to option pricing problems can be found in [15] and [12].

Barndorff-Nielsen [1] proposed to model the log returns of asset prices by using the normal inverse Gaussian (NIG) distribution. This family of non-normal distributions has proven to fit the semi-heavily tails observed in financial time series of various kinds extremely well (see Rydberg [21] or Eberlein and Keller [7]). The time dynamics of the asset prices are modelled by an exponential L´evy process. To price such derivatives, even simple call and put options, we need to consider the numerical evaluation of the expectation. Raible [18] has considered a Fourier method to evaluate call and put options. Alternatives to this method are the MC method or the QMC method. The QMC method has been applied with succes in financial applications by many authors (see [8]), and has strong convergence properties. Majority of the work done on applying these simulation techniques to financial problems was in direction where one needs to simulate from the normal distribution. One exception is Kainhofer (see [13]), who proposes a QMC algorithm for NIG variables, based on a technique proposed by Hlawka and M¨uck (see [11]) to generate low-discrepancy sequences for general distributions.

In a recent paper [19], we proposed a mixed MC and QMC method for es- timating ans-dimensional integral I and we defined a new hybrid sequence that we called the H-mixed sequence. Other authors who combine the ideas from MC and QMC methods in estimating multidimensional integrals are G. ¨Okten (see [16]) and N. Ro¸sca (see [20]). Using these sequences, we defined a new estimator and proved a central limit theorem for this estimator. In this paper, we apply our mixed method to practical problems from financial mathematics. First, we remember the theoreti- cal background of our method and give some important results. Then, we apply the H-mixed sequence to valuation of an European Call option and compare the effec- tiveness of it with that of pseudorandom and low-discrepancy sequences. At the end, we apply the mixed method to a more difficult problem from finance, namely the

valuation of Asian options. We also compare numerically our method with the MC and QMC methods.

2. H-mixed sequences

Let us consider the problem of estimating integrals of the form I=

Z

[0,1]^s

f(x)dH(x), (1)

where f : [0,1]^s →R is the function we want to integrate andH :R^s → [0,1] is a distribution function on [0,1]^s. In the continuous case, the integralIcan be rewritten as

I= Z

[0,1]^s

f(x)h(x)dx,

wherehis the density function corresponding to the distribution functionH.

In the MC method (see [22]), the integralIis estimated by sums of the form IˆN = 1

N

N

X

k=1

f(xk),

where xk = (x⁽¹⁾_k , . . . , x^(s)_k ), k ≥ 1, are independent identically distributed random points on [0,1]^s, with the common density functionh.

In the QMC method (see [22]), the integralIis approximated by sums of the form _N¹ PN

k=1f(x_k), where (x_k)_k≥1 is a H-distributed low-discrepancy sequence on [0,1]^s.

Next, the notions of discrepancy and marginal distributions are introduced.

Definition 1 (H-discrepancy). Consider an s-dimensional continuous distribution on [0,1]^s, with distribution function H. Let λ_H be the probability measure induced by H. Let P = (x1, . . . , xN)be a sequence of points in [0,1]^s. The H-discrepancy of sequence P is defined as

DN,H(P) = sup

J⊆[0,1]^s

1

NAN(J, P)−λH(J) ,

where the supremum is calculated over all subintervals J = Qs

i=1[a_i, b_i] ⊆ [0,1]^s; AN(J, P)counts the number of elements of the set (x1, . . . , xN), falling into the in- terval J, i.e.

AN(J, P) =

N

X

k=1

1J(xk).

1J is the characteristic function ofJ.

The sequenceP is called H-distributed ifD_N,H(P)→0 asN → ∞.

The H-distributed sequence P is said to be a low-discrepancy sequence if DN,H(P) =O (logN)^s/N

.

The non-uniform Koksma-Hlawka inequality ([3]) gives an upper bound for the approximation error in QMC integration, when H-distributed low-discrepancy sequences are used.

Theorem 2 (non-uniform Koksma-Hlawka inequality). Let f : [0,1]^s → R be a function of bounded variation in the sense of Hardy and Krause and(x1, . . . , xN) be a sequence of points in[0,1]^s. Consider an s-dimensional continuous distribution on [0,1]^s, with distribution functionH. Then, for anyN >0

Z

[0,1]^s

f(x)dH(x)− 1 N

N

X

k=1

f(xk)

≤VHK(f)DN,H(x1, . . . , xN), (2)

whereVHK(f)is the variation off in the sense of Hardy and Krause.

Definition 3. Consider ans-dimensional continuous distribution on[0,1]^s, with den- sity functionhand distribution functionH. For a pointu= u⁽¹⁾, . . . , u^(s)

∈[0,1]^s, the marginal density functionshl,l= 1, . . . , s, are defined by

hl u^(l)

= Z

. . . Z

| {z }

[0,1]^s−1

h t⁽¹⁾, . . . , t^(l−1), u^(l), t^(l+1), . . . t^(s)

dt⁽¹⁾. . . dt^(l−1)dt^(l+1). . . dt^(s),

and the marginal distribution functions Hl,l= 1, . . . , s, are defined by

Hl u^(l)

= Z u^(l)

0

hl(t)dt.

We considers-dimensional continuous distributions on [0,1]^s, with indepen- dent marginals, i.e.,

H(u) =

s

Y

l=1

Hl(u^(l)), ∀u= (u⁽¹⁾, . . . , u^(s))∈[0,1]^s.

This can be expressed, using the marginal density functions, as follows:

h(u) =

s

Y

l=1

hl(u^(l)), ∀u= (u⁽¹⁾, . . . , u^(s))∈[0,1]^s.

Consider an integer 0 < d < s. Using the marginal density functions, we construct the following density functions on [0,1]^d and [0,1]^s−d, respectively:

h_q(u) =

d

Y

l=1

h_l(u^(l)), ∀u= (u⁽¹⁾, . . . , u^(d))∈[0,1]^d, and

h_X(u) =

s

Y

l=d+1

h_l(u^(l)), ∀u= (u^(d+1), . . . , u^(s))∈[0,1]^s−d. The corresponding distribution functions are

H_q(u) = Z u⁽¹⁾

0

. . . Z u^(d)

0

h_q t⁽¹⁾, . . . , t^(d)

dt⁽¹⁾. . . dt^(d),

u= (u⁽¹⁾, . . . , u^(d))∈[0,1]^d, (3) and

H_X(u) =

Z u^(d+1) 0

. . . Z u^(s)

0

h_X t^(d+1), . . . , t^(s)

dt^(d+1). . . dt^(s),

u= (u^(d+1), . . . , u^(s))∈[0,1]^s−d. (4) Next, we introduce the new notion of aH-mixed sequence.

Definition 4(H-mixed sequence). ([19])

Consider ans-dimensional continuous distribution on[0,1]^s, with distribution function H and independent marginals Hl, l = 1, . . . , s. Let Hq and HX be the distribution functions defined in (3) and (4), respectively.

Let (qk)_k≥1 be a Hq-distributed low-discrepancy sequence on [0,1]^d, with qk = (q_k⁽¹⁾, . . . , q^(d)_k ), and Xk, k ≥ 1, be independent and identically distributed random vectors on[0,1]^s−d, with distribution functionH_X, whereX_k= (X_k^(d+1), . . . , X_k^(s)).

A sequence(m_k)_k≥1, with the general term given by

m_k = (q_k, X_k), k≥1, (5)

is called a H-mixed sequence on[0,1]^s. Remark 5. For an interval J = Qs

l=1[al, bl] ⊆ [0,1]^s, we define the subintervals J⁰=Qd

l=1[al, bl]⊆[0,1]^d andJ⁰⁰=Qs

l=d+1[al, bl]⊆[0,1]^s−d (i.e. J =J⁰×J⁰⁰).

Let (m_k)_k≥1 be a H-mixed sequence on [0,1]^s, with the general term given by (5).

Based on definitions (1) and (4), the H-discrepancy of the sequence (m1, . . . , mN) can be expressed as

D_N,H(m₁, . . . , m_N) = sup

J⊆[0,1]^s

1 N

N

X

k=1

1_J(m_k)−

s

Y

l=1

[H_l(b_l)−H_l(a_l)]

, and theHq-discrepancy of the sequence (q1, . . . , qN) is given by

D_N,Hq(q₁, . . . , q_N) = sup

J⁰⊆[0,1]^d

1 N

N

X

k=1

1_J⁰(q_k)−

d

Y

l=1

[H_l(b_l)−H_l(a_l)]

. The following result gives a probabilistic error bound for the H-mixed se- quences.

Theorem 6. ([19]) If(m_k)_k≥1= (q_k, X_k)_k≥1 is aH-mixed sequence, then∀ε >0we have

P

DN,H(m1, . . . , mN)≤ε+DN,H_q(q1, . . . , qN)

≥1−1 ε²

1 4N

DN,H_q(q1, . . . , qN)+1 . (6) In order to estimate general integrals of the form (1), we introduce the fol- lowing estimator.

Definition 7. ([19]) Let (mk)_k≥1 = (qk, Xk)_k≥1 be an s-dimensional H-mixed se- quence, introduced by us in Definition 4, with qk = (q_k⁽¹⁾, . . . , q^(d)_k ) and Xk = (X_k^(d+1), . . . , X_k^(s)). We define the following estimator for the integralI:

θm= 1 N

N

X

k=1

f(mk). (7)

We consider the independent random variables:

Y_k =f(m_k) =f(q_k⁽¹⁾, . . . , q_k^(d), X_k^(d+1), . . . , X_k^(s)), k≥1. (8) We denote the expectation ofYk by

E(Yk) =µk, (9)

and the variance ofYk by

V ar(Yk) =σ_k². (10)

We assume that

0< σ²_k<∞, (11)

and we denote

0< σ²_(N)=σ²₁+. . .+σ_N² <∞. (12) In what follows, we give and prove an important result, concerning the esti- mator (7) introduced previously by us.

Theorem 8. Let (m_k)_k≥1 = (q_k, X_k)_k≥1 be an s-dimensional H-mixed sequence, defined in (5). We assume that the integrant f is bounded on [0,1]^s and that the function

g(x⁽¹⁾, . . . , x^(d)) = Z

[0,1]^s−d

f(x⁽¹⁾, . . . , x^(s)) Y^s

l=d+1

h_l(x^(l))

dx^(d+1)·. . .·dx^(s),

is of bounded variation in the sense of Hardy and Krause. Then, the estimator θ_m, defined in relation (7), is asymptotically unbiased i.e.,

E(θ_m)→I, as N → ∞.

Proof. As (qk)k≥1, with qk = (q⁽¹⁾_k , . . . , q_k^(d)), is a Hq-distributed low-discrepancy sequence on [0,1]^d, it follows that

D_N,Hq(q₁, . . . , q_N)→0, as N → ∞. (13)

Using this and the fact that functiongis of bounded variation in the sense of Hardy and Krause, it follows from Koksma-Hlawka inequality (2) that

1 N

N

X

k=1

g(q_k⁽¹⁾, . . . , q^(d)_k )−→

Z

[0,1]^d

g(x⁽¹⁾, . . . , x^(d))Y^d

l=1

h_l(x^(l))

dx⁽¹⁾. . . dx^(d)

= Z

[0,1]^d

hZ

[0,1]^s−d

f(x⁽¹⁾, . . . , x^(s)) Y^s

l=d+1

hl(x^(l))

dx^(d+1)·. . .·dx^(s)i

·Y^d

l=1

hl(x^(l))

dx⁽¹⁾. . . dx^(d)

= Z

[0,1]^s

f(x⁽¹⁾, . . . , x^(s))Y^s

l=1

hl(x^(l))

dx⁽¹⁾·. . .·dx^(s)=I, as N → ∞.

The expectation of our estimator is E(θm) =E 1

N

N

X

k=1

f(mk)

!

= 1 N

N

X

k=1

E(f(q_k⁽¹⁾, . . . , q_k^(d), X_k^(d+1), . . . , X_k^(s)))

= 1 N

N

X

k=1

Z

[0,1]^s−d

f(q_k⁽¹⁾, . . . , q^(d)_k , x^(d+1), . . . , x^(s)) Y^s

l=d+1

hl(x^(l))

dx^(d+1)·. . .·dx^(s)

= 1 N

N

X

k=1

g(q⁽¹⁾_k , . . . , q_k^(d)).

Hence, we get in the end that

N→∞lim E(θ_m) =I.

We call the method of estimating the integralI, based on the estimatorθm, defined in (7), the mixed method.

Proposition 9. ([19]) Let (mk)k≥1= (qk, Xk)k≥1 be an s-dimensional H-mixed se- quence. We assume that f is bounded on[0,1]^sand that the functions

f1(x⁽¹⁾, . . . , x^(d)) = Z

[0,1]^s−d

(f(x⁽¹⁾, . . . , x^(s)))² Y^s

l=d+1

hl(x^(l))

dx^(d+1)·. . .·dx^(s),

f2(x⁽¹⁾, . . . , x^(d)) =hZ

[0,1]^s−d

f(x⁽¹⁾, . . . , x^(s)) Y^s

l=d+1

hl(x^(l))

dx^(d+1)·. . .·dx^(s)i²

are of bounded variation in the sense of Hardy and Krause. Then, we have σ_(N)²

N −→L, as N −→ ∞,

where

L= Z

[0,1]^s

(f(x⁽¹⁾, . . . , x^(s)))²Y^s

l=1

h_l(x^(l))

dx⁽¹⁾. . . dx^(s)−

− Z

[0,1]^d

hZ

[0,1]^s−d

f(x⁽¹⁾, . . . , x^(s)) Y^s

l=d+1

hl(x^(l))

dx^(d+1)·. . .·dx^(s)i2

·Y^d

l=1

hl(x^(l))

dx⁽¹⁾. . . dx^(d).

Another important result regarding the estimator defined before is recalled next.

Theorem 10. ([19]) In the same hypothesis as in Proposition 9 and, in addition, assuming thatL6= 0, we have

a)

Y_(N)= PN

k=1Yk−PN k=1µk

σ(N)

−→Y, as N → ∞, (14) where the random variableY has the standard normal distribution.

b) If we denote the crude Monte Carlo estimator for the integral (1) byθM C, then

V ar(θm)≤V ar(θM C), (15)

meaning that, by using our estimator, we obtain asymptotically a smaller variance than by using the classical Monte Carlo method.

3. Application to finance: European options

In this section, we apply our mixed method to a problem from mathematical finance. The general setting of the problem is presented next. We consider the situation where the stock price of the underlaying assetS=S(t) is driven by a L´evy processL(t),

S(t) =S(0)e^L(t). (16)

L´evy processes can be characterized by the distribution of the random variableL(1).

This distribution can be hyperbolic (see [7]), normal inverse gaussian (NIG), variance- gamma (see [14]), or Meixner distribution.

According to the fundamental theory of asset pricing (see [5] ), the risk- neutral price of an option,C(0), is given by

C(0) =e^−rTE^Q(C_T(S)), (17)

where C_T(S) is the so-called payoff of the derivative, which coincides with its value at expiration or exercise timeT, andQis an equivalent martingale measure. In this paper, we are mostly concerned with exponential NIG-L´evy processes, meaning that L(t) has independent increments, distributed according to a NIG distribution. For a detailed discussion of the NIG distribution and the corresponding L´evy process, we refer to Barndorff-Nielsen [1] and Rydberg [21]. In the situation of exponential NIG-L´evy models, we have an incomplete market, leading to a continuum of equiv- alent martingale measures Q, which can be used for risk-neutral pricing. Here, we choose the approach of Raible [18] and consider the measure obtained by Esscher transform method. This approach is so-called structure preserving, in the sense that the distribution ofL(1) remains in the class of NIG distributions.

In the following, we consider the evaluation of so-called European Call op- tions, which have to be valued by simulation. The risk-neutral price of such an option is

C(0) =e^−rTE^Q(max{S(T)−K,0}) =e^−rTE^Q((S(T)−K)+), (18) where the constantKis called the strike price. If we replace the stock price by (16), we obtain

C(0) =e^−rTE^Q((S(0)e^L(T)−K)₊). (19)

From practice, we know that the evaluation of the stock priceS(t) is made at discrete times 0 = t0 < t1 < t2 < . . . < ts = T. For simplicity, we focus on regular time intervals, ∆t=t_i−t_i−1. We note that

X_i=L(t_i)−L(t_i−1) =L(t_i−1+ ∆t)−L(t_i−1)∼L(∆t), i= 1, . . . , s,

are independent and identically distributed NIG random variables with the same dis- tribution asL(t1). Dropping the discounted factor from the risk-neutral option price, we get the expected payoff under the Esscher transform measure of the European Call option

E^Q((S(0)e^L(T)−K)₊) =E((S(0)e^Psi=1Xi−K)₊), (20) Our purpose is to evaluate the expected payoff (20). For this, we rewrite the expectation (20) as a multidimensional integral onR^s

I= Z

R^s

S(0)e^Psi=1x(i)−K

+

| {z }

E(x)

dG(x) = Z

R^s

E(x)dG(x), (21)

where G(x) = Qs

i=1Gi(x⁽ⁱ⁾), ∀x = (x⁽¹⁾, . . . , x^(s)) ∈ R^s, and Gi(x⁽ⁱ⁾) denotes the distribution function of the so-called log returns induced by L(t₁), with the corre- sponding density functiongi(x⁽ⁱ⁾). These log increments are independent and NIG distributed, having a common probability density

fN IG(x;µ, β, α, δ) =α

πexp δp

α²−β²+β(x−µ)δK1(αp

δ²+ (x−µ)²) pδ²+ (x−µ)² (22) whereK₁(x) denotes the modified Bessel function of third type of order 1 (see [17]).

In order to approximate the integral (21), we have to transform it to an integral on [0,1]^s. We can do this using an integral transformation, as follows.

We first consider the family of independent double-exponential distributions with the same parameterλ > 0, having the cumulative distribution functionsG_λ,i: R→[0,1], i= 1, . . . , s,

G_λ,i(x) =







1

2e^λx , x <0 1−¹₂e^−λx , x≥0,

(23)

and the inversesG⁻¹_λ,i: [0,1]→R, i= 1, . . . , s, given by

G⁻¹_λ,i(x) =







1

λlog (2x) , x≤ ¹₂

−_λ¹log (2−2x) , x > ¹₂.

(24)

Next, we consider the substitutions x⁽ⁱ⁾ = G⁻¹_λ,i(1−y⁽ⁱ⁾), i = 1, . . . , s, and then takey⁽ⁱ⁾= 1−z⁽ⁱ⁾, i= 1, . . . , s.

The integral (21) becomes I=

Z

[0,1]^s

S(0)e^{Psi=1G−1λ,i(z(i))}−K

+

| {z }

f(z)

dH(z) = Z

[0,1]^s

f(z)dH(z), (25)

whereH : [0,1]^s→[0,1], defined by H(z) =

s

Y

i=1

(G_i◦G⁻¹_λ,i)(z⁽ⁱ⁾), ∀z= (z⁽¹⁾, . . . , z^(s))∈[0,1]^s, (26)

is a distribution function on [0,1]^s, with independent marginals Hi = Gi ◦ G⁻¹_λ,i, i= 1, . . . , s.

In the following, we compare numerically our mixed method with the MC and QMC methods. As a measure of comparison, we will use the absolute errors produced by these three methods in the approximation of the integral (25).

The MC estimate is defined as follows:

θM C = 1 N

N

X

k=1

f(x⁽¹⁾_k , . . . , x^(s)_k ), (27)

where xk = (x⁽¹⁾_k , . . . , x^(s)_k ), k ≥ 1, are independent identically distributed random points on [0,1]^s, with the common distribution functionH defined in (26).

In order to generate such a pointxk, we proceed as follows. We first generate a random pointω_k = (ω_k⁽¹⁾, . . . , ω_k^(s)), whereω⁽ⁱ⁾_k is a point uniformly distributed on [0,1], for eachi= 1, . . . , s. Then, for each componentω_k⁽ⁱ⁾,i= 1, . . . , s, we apply the inversion method (see [4] and [6]), and obtain thatH_i⁻¹(ω_k⁽ⁱ⁾) = (Gλ,i◦G⁻¹_i )(ω⁽ⁱ⁾_k ) is a point with the distribution functionH_i. As the s-dimensional distribution with the distribution functionH has independent marginals, it follows that xk = ((Gλ,1◦ G⁻¹₁ )(ω⁽¹⁾_k ), . . . ,(Gλ,s◦G⁻¹_s )(ω^(s)_k )) is a point on [0,1]^s, with the distribution function H. As we can see, in order to generate non-uniform random points on [0,1]^s, with distribution function H, we need to know the inverse of the distribution function of a NIG distributed random variable or, at least an approximation of it. As the inverse function is not explicitly known, an approximation of it is needed in our simulations. In order to obtain an approximation of the inverse, we use the Matlab

function ”niginv” as implemented by R. Werner, based on a method proposed by K.

Prause in his Ph.D. dissertation [17].

The QMC estimate is defined as follows:

θQM C = 1 N

N

X

k=1

f(x⁽¹⁾_k , . . . , x^(s)_k ), (28) wherex= (xk)_k≥1 is aH-distributed low-discrepancy sequence on [0,1]^s, with xk = (x⁽¹⁾_k , . . . , x^(s)_k ),k≥1.

In order to generate such a sequence, we apply a method proposed by Hlawka and M¨uck in [11]. In their method, they create directlyH-distributed low-discrepancy sequences, whereH can be any distribution function on [0,1]^s, with density function h, which can be factored into a product of independent, one-dimensional densities.

The method is based on the following theoretical result.

Theorem 11. ([10]) Consider an s-dimensional continuous distribution on [0,1]^s, with distribution function H and density function h(u) = Qs

j=1hj(u^(j)), ∀u = u⁽¹⁾, . . . , u^(s)

∈ [0,1]^s. Assume that hj(t) 6= 0, for almost every t ∈ [0,1] and for all j = 1, . . . , s. Furthermore, assume that h_j, j = 1, . . . , s, are continuous on [0,1]. Denote byM_f = sup_u∈[0,1]sf(u). Letω= (ω₁, . . . , ω_N)be a sequence in [0,1]^s. Generate the sequencex= (x1, . . . , xN), with

x^(j)_k = 1 N

N

X

r=1

1 +ω^(j)_k −Hj ω_r^(j)

= 1 N

N

X

r=1

1_[0,ω(j)

k ] Hj ω^(j)_r

, (29)

for allk= 1, . . . , N and allj= 1, . . . , s, where[a]denotes the integer part ofa. Then the generated sequence xhas a H-discrepancy of

DN,H(x1, . . . , xN)≤(2 + 6sMf)DN(ω1, . . . , ωN).

As our distribution functionH can be factored into independent marginals, and has the support on [0,1]^s, we can apply directly the above theorem, to generate H-distributed low-discrepancy sequences. During our experiments, we employed as low-discrepancy sequencesω= (ωk)_k≥1on [0,1]^s, the Halton sequences (see [9]).

All points constructed by the Hlawka-M¨uck method are of the form i/N, i = 0, . . . , N, in particular some elements of the sequence x = (x₁, . . . , x_N) might

assume a value of 0 or 1. A value of 1 is a singularity of the function f(x), due to the logarithm from the definition of G⁻¹_λ,i(x), which becomes unbounded if x = 1.

Hence, the sequence constructed with Hlawka-M¨uck method is not directly suited for unbounded problems. To overcome this problem, Kainhofer (see [13]) suggests to define a new sequence, in which the value 1 is replaced by 1/N, where N is the number of points in the set. This slight modification of the sequence is shown to have a minor influence, as the transformed set does not loose its low-discrepancy and can be used for QMC integration.

The estimate proposed by us earlier is:

θm= 1 N

N

X

k=1

f(q_k⁽¹⁾, . . . , q^(d)_k , X_k^(d+1), . . . , X_k^(s)), (30) where (qk, Xk)_k≥1 is an s-dimensionalH-mixed sequence on [0,1]^s.

In order to obtain such a H-mixed sequence, we first construct the Hq- distributed low-discrepancy sequence (q_k)_k≥1 on [0,1]^d, using the Hlawka-M¨uck method (the distribution function H_q was defined in (3)). Next, we generate the independent and identically distributed random points xk, k ≥1 on [0,1]^s−d, with the common distribution functionHX, using the inversion method (the distribution functionH_X was defined in (4)). Finally, we concatenate q_k andx_k for eachk≥1, and get ourH-mixed sequence on [0,1]^s.

In our experiments, we used as low-discrepancy sequences on [0,1]^d, for the generation ofH-mixed sequences, the Halton sequences (see [9]).

We suppose that the parameters of the NIG-distributed log-returns under the equivalent martingale measure given by the Esscher transform are given by

µ= 0.00079∗5, β=−15.1977, α= 136.29, δ= 0.0059∗5, (31) and they are the same as in Kainhofer (see [13]). We observe that these parameters are relevant for daily observed stock price log-returns (see [21]). As the class of NIG distributions is closed under convolution, we can derive weekly stock prices by using a factor of 5 for the parameters µ andδ. We suppose further that the initial stock price isS(0) = 100 and the risk-free annual interest rate isr= 3.75%.

The option is sampled at weekly time intervals. We also let the option to have maturities of 12 and 20 weeks. Hence, our problem is a 12 and 20-dimensional integral, respectively, over the payoff function.

We are going to compare the three estimates in terms of their absolute error, where the ”exact” option price is obtained as the average of 10 MC simulations, with N = 100000 for the initial integral (21).

In our tests we have considered the following dimensions of the transformed integral (25) on [0,1]^s: s= 12, 20. The MC and H-mixed estimates are the mean values of 10 independent runs, while the QMC estimate is the result of a single run.

The results are presented in two tables, each table containing the number of samples N, which varies from 5000 to 8500 with a step of 500, and the absolute error of the three estimates.

N Absolute error MC Absolute error QMC Absolute error Mixed Method

5000 0.014731 0.012385 0.007676

5500 0.004485 0.016085 0.003780

6000 0.009268 0.011866 0.001892

6500 0.020887 0.014721 0.002547

7000 0.027395 0.014732 0.008411

7500 0.006316 0.012404 0.017385

8000 0.015027 0.010519 0.012538

8500 0.009207 0.010140 0.007248

Table 1: European Call Option. Cased= 4 ands= 12.

The numerical results for s = 12 and d= 4 are presented in Table 1. The results produced by our H-mixed sequence are much better than the ones obtained by using pseudorandom or low-discrepancy sequences, in almost all situations.

N Absolute error MC Absolute error QMC Absolute error Mixed Method

5000 0.006381 0.049304 0.004311

5500 0.030035 0.039222 0.008018

6000 0.017018 0.042674 0.019373

6500 0.004735 0.041674 0.012044

7000 0.023534 0.038581 0.013131

7500 0.020561 0.030509 0.001833

8000 0.028440 0.027873 0.008792

8500 0.012737 0.032972 0.007264

Table 2: European Call Option. Cased= 7 ands= 20.

To increase the difficulty of the problem, we increase the dimension of the integral to s = 20. Table 2 displays the results we get for s = 20 and d = 7.

From this simulations, we see again that theH-mixed sequence outperforms both the pseudorandom and low-discrepancy sequences, for almost all sample sizes N. The absolute error produced by ourH-mixed sequence is smaller than the one produced by the low-discrepancy sequence, in all situations.

As a general conclusion for this option pricing problem, we can say that by the use ofH-mixed sequences, we obtain increasing advantages over the classical pseudo- random and low-discrepancy sequences, for relatively high dimensions and moderate sample sizes.

4. Application to finance: Asian options

In this section, we consider an Asian option pricing problem. We compare numerically our mixed method with the MC and QMC methods, when they are applied to so-called (discrete sampled) Asian options driven by the asset dynamics S(t), as defined in (16). The general setting remains the same as in the previous section, but the payoff function is changed. The payoff of an Asian call option is defined as

C_T(S) =1 s

s

X

i=1

S(t_i)−K

+

= maxn1 s

s

X

i=1

S(t_i)−K,0o

, (32)

with 0 =t₀< t₁ < t₂< . . . < t_s=T. The constantK≥0 is called the strike price.

Hence, we get the following integration problem:

I= Z

R^s

S(0) s

s

X

i=1

e^Pij=1x(j)−K

+

| {z }

A(x)

dG(x) = Z

R^s

A(x)dG(x), (33)

where G(x) = Qs

i=1G_i(x⁽ⁱ⁾), ∀x = (x⁽¹⁾, . . . , x^(s)) ∈ R^s, and G_i(x⁽ⁱ⁾) denotes the distribution function of the so-called log returns induced by L(t1), with the corre- sponding density functiongi(x⁽ⁱ⁾). These log increments are independent and NIG distributed, having the common density function defined in (22).

In order to approximate the integral (33), we have to transform it to an integral on [0,1]^s. We can do this in a similar way as we did for European Call options, in the previous section. In the end, we get the following integration problem on [0,1]^s:

I= Z

[0,1]^s

S(0) s

s

X

i=1

e^{Pij=1G−1λ,i(z(j))}−K

+

| {z }

f(z)

dH(z) = Z

[0,1]^s

f(z)dH(z), (34)

whereH : [0,1]^s→[0,1], defined by H(z) =

s

Y

i=1

(G_i◦G⁻¹_λ,i)(z⁽ⁱ⁾), ∀z= (z⁽¹⁾, . . . , z^(s))∈[0,1]^s, (35) is a distribution function on [0,1]^s, with independent marginals H_i = G_i ◦ G⁻¹_λ,i, i= 1, . . . , s.

Next, we compare numerically our estimatorθm, with the estimators obtained using the MC and QMC methods. All three estimatorsθ_m,θ_{M C} andθ_{QM C}, and the corresponding sequences are defined in the previous section. The function f(z) is defined in relation (34). As a measure of comparison, we will use the absolute errors produced by these three methods, in the approximation of the integral (34).

We suppose that the parameters of the NIG-distributed log-returns under the equivalent martingale measure given by the Esscher transform are the same as in (31). We assume that the initial stock price isS(0) = 100, and the risk-free annual interest rate isr= 3.75%.

For our mixed method and QMC estimate, we use a Halton sequence as low- discrepancy sequence on [0,1]^s. The Asian call option is sampled weekly. We also let the option to have maturities of 12 and 20 weeks. Hence, our problem is a 12 and 20-dimensional integral, respectively, over the payoff function.

We are going to compare the three estimates in terms of their absolute error, where the ”true” price is obtained as the average of 10 MC simulations, with N = 100000. The MC andH-mixed estimates are the mean values of 10 independent runs, while the QMC estimate is the result of a single run. The results are presented in two tables, each table containing the number of samplesN, which varies from 4000 to 7000 with a step size of 500, and the absolute error of the three estimates.

N Absolute error MC Absolute error QMC Absolute error Mixed Method

4000 0.004833 0.000723 0.000690

4500 0.003060 0.001083 0.000977

5000 0.001095 0.000380 0.001653

5500 0.000293 0.000618 0.000599

6000 0.011389 0.001482 0.000898

6500 0.001733 0.003187 0.000218

7000 0.008720 0.001582 0.000047

Table 3: Asian Option. Cased= 4 and s= 12.

In Table 3 we present the results obtained fors = 12 and d = 4. The H- mixed sequence gives excellent estimates for almost all N, clearly dominating both the pseudorandom and low-discrepancy sequences.

N Absolute error MC Absolute error QMC Absolute error Mixed Method

4000 0.036868 0.014541 0.004318

4500 0.007101 0.011003 0.016654

5000 0.002396 0.009723 0.012004

5500 0.016070 0.008897 0.000818

6000 0.004666 0.008920 0.000003

6500 0.017100 0.009541 0.003007

7000 0.010705 0.009437 0.002776

Table 4: Asian Option. Cased= 7 and s= 20.

The estimates presented in Table 4 are the results of the simulations for a higher dimensional problem, with s = 20 andd = 7. Again, the H-mixed method outperforms the conventional MC and QMC methods, in almost all situations.

We can conclude that our mixed method can give considerable improvements over the MC and QMC methods, in estimating high dimensional integrals, which we encounter in problems from financial mathematics, such as valuation of Asian options and European options.

References

[1] Barndorff-Nielsen, O.E., Processes of normal invers Gaussian type, Finance and Sto- chastic,2(1998), 41-68.

[2] Boyle, P.P.,Options: A Monte Carlo Approach, J. Financial Economics,4(1977), 323- 338.

[3] Chelson, P., Quasi-random techniques for Monte Carlo methods, Ph.D. Dissertation, The Claremont Graduate School, 1976.

[4] Deak, I.,Random number generators and simulation, Akademiai Kiado, Budapest, 1990.

[5] Delbaen, F., Schachermayer, W.,A general version of the fundamental theorem of assest pricing, Math. Ann.,300(3)(1994), 463-520.

[6] Devroye, L.,Non-uniform random variate generation, Springer-Verlag, New-York, 1986.

[7] Eberlein, E., Keller, U.,Hyperbolic distribution in finance, Bernoulli,1(1995), 281-299.

[8] Glasserman, P.,Monte Carlo methods in financial engineering, Springer-Verlag, New- York, 2003.

[9] Halton, J.H.,On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals, Numer. Math.,2(1960), 84-90.

[10] Hlawka, E.,Gleichverteilung und Simulation, ¨Osterreich. Akad. Wiss. Math.-Natur. Kl.

Sitzungsber. II,206(1997), 183-216.

[11] Hlawka, E., M¨uck, R., A transformation of equidistributed sequences, Applications of number theory to numerical analysis, Academic Press, New-York, 1972, 371-388.

[12] Joy, C., Boyle, P.P., Tan, K.S.,Quasi-Monte Carlo methods in numerical finance, Man- agement Science, Vol. 42, No.6(1996), 926-938.

[13] Kainhofer, R.F.,Quasi-Monte Carlo algorithms with applications in numerical analysis and finance, Ph.D. Dissertation, Graz University of Technology, Austria, 2003.

[14] Madan, D.B., and Seneta, E., The variance gamma (V.G.) model for share market returns, Journal of Business,63(1990), 511-524.

[15] Ninomya, S., Tezuka, S.,Toward real time pricing of complex finacial derivatives, Ap- plied Mathematical Finance, Vol.3(1996), 1-20.

[16] ¨Okten, G., Tuffin, B., Burago, V.,A central limit theorem and improved error bounds for a hybrid-Monte Carlo sequence with applications in computational finance, Journal of Complexity, Vol. 22, No.4(2006), 435-458.

[17] Prause, K.,The generalized hyperbolic model: estimation, financial derivatives and risk measures, Ph.D. Dissertation, Albert-Ludwigs-Universitat, Freiburg, 1999.

[18] Raible, S.,L´evy processes in finance: theory, numerics and empirical facts, Ph.D. Dis- sertation, Albert-Ludwigs-Universitat, Freiburg, 2000.

[19] Ro¸sca, A.V., A mixed Monte Carlo and quasi-Monte Carlo sequence for multidimen- sional integral estimation, Acta Universitatis Apulensis, Mathematics-Informatics, No.

14(2007), 141-160.

[20] Ro¸sca, N.,A combined Monte Carlo and quasi-Monte Carlo method for estimating mul- tidimensional integrals, Studia Universitatis Babes-Bolyai, Mathematica, Vol. LII, No.

1(2007), 125-140.

[21] Rydberg, T.H.,The normal inverse Gaussian L´evy process: simulation and approxima- tion, Comm. Statist. Stochastic Models, Vol. 13, No.4(1997), 887-910.

[22] Stancu, D.D., Coman, Gh., Blaga, P.,Numerical analysis and approximation theory, (in Romanian) Vol. II, Presa Universitar˘a Clujean˘a, Cluj-Napoca, 2002.

Babes¸-Bolyai University

Faculty of Economics and Business Administration Str. Teodor Mihali, Nr. 58-60, Cluj-Napoca, Romania E-mail address: [email protected]