Variance Reduction Technique: Importance Sampling Pricing minimum return guarantee in Unit-linked life insurance

(1)

Pricing minimum return guarantee in Unit-linked life insurance

Ran Chen

Student number: 1432362

University of Groningen, Faculty of Economics

Master thesis Econometrics, Operations Research & Actuarial Studies Specialization: Actuarial Studies

Supervisors University of Groningen: Dr. L. Spierdijk

Prof. dr. R.H. Koning Supevisors Watson Wyatt:

Drs. W. Elshof AAG Drs. S. van de Pas AAG

(2)

Preface

I would like to thank several people who guided me during my research. My thanks goes to my supervisors from Watson Wyatt Drs. W. Elshof AAG and Drs. S. van de Pas AAG, as well as supervisors from the University of Groningen Dr. L. Spierdijk and Prof. dr. R.H. Koning. In particular, i would like to thank my main supervisors Elshof and Spierdijk for their valuable advices. I would also thank Watson Watson Insurance Consulting B.V. for offering me an internship.

(3)

This thesis discusses a well-known variance reduction technique (VRT) called importance sampling (IS). This technique is used to reduce run times in Monte Carlo (MC) simulation. The MC method is used when it is impossible or too complicated to compute results with an analytical method. This method has applications in many areas, such as finance, statistical physics, computer science and many other areas.

(4)

In the past, insurance companies made high investment returns because of the good perfor-mance in the financial market. The returns were much higher than the risk-free interest rate. To make insurance products more interesting to the customers, many insurance companies offer a minimum return guarantee (MRG) rate for some life insurance products. The guar-antee rate is usually near the risk-free interest rate. Traditionally, this guarguar-antee is given for free to the policy holders because the returns from the risk-free interest rarely outperformed the investment returns. As a result the insurer bears a risk without any compensation. A well-known example of life insurance products is the Unit-linked (UL) life insurance. This insurance gives the policy holder the possibility to invest the insurance premium in a invest-ment fund. Usually, the investinvest-ment fund consists of different assets such as equities and bonds. Offering ‘free’ MRG to the policy holder was never a real problem for the insurance com-panies, since the high investment returns in the past almost always outperformed the MRG rate. But in recent years, MRG have become more and more valuable because of decreasing stock prices. Because of this, many insurance companies are surprised by the costs of the MRG. These days, pricing embedded option in an insurance policy such as the MRG is a major concern for the insurance companies.

It is often difficult to price the embedded option such as the MRG in a life insurance product. Compared to vanilla options, such as the European call and European put, embedded options are more complex. This is because the payoff of the embedded option depends not only on the future returns of underlying assets, but it also depends on the mortality rate, risk premium etc. Due to the complexity of these options, the price of only few types of these options can be calculated analytically.

Another method to price the MRG in a life insurance product is the Monte Carlo simulation. The advantage of this method is that it can be applied to price complex embedded options when no explicit formula is available. However, this method has also a big disadvantage: to use the Monte Carlo method many simulation scenarios are needed. This results in long run times, since the insurance companies have huge amounts of insurance contracts in their portfolio.

A group of well-known techniques called the variance reduction techniques (VRTs) can be used to reduce the simulation run time without losing the precision of the estimator. In this thesis we mainly focus on one VRT called importance sampling (IS). This technique has been applied successfully to price many financial products, such as the Asian options. However, this method is rarely investigated to price the embedded options such as the MRG in a UL life insurance. The main goal of this thesis is to investigate importance sampling to price the MRG in a UL life insurance. The main research questions are:

• How can we use importance sampling to reduce the simulation run time of the Monte Carlo simulation?

(6)

1 INTRODUCTION

• Can we use the technique importance sampling to price a portfolio of unit-linked insur-ance policies?

(7)

2.1 Introduction

In the theory of Monte Carlo (MC) methods, variance reduction techniques (VRTs) are used to improve the efficiency and speed of an MC simulation. In this section we discuss some of these techniques. In particular, we focus on one technique called importance sampling (IS). This technique will be discussed in detail in Section 2.3. After this, we will also discuss other well-known VRTs in Section 2.4. Next, we start with an introduction to Monte Carlo methods.

2.2 Monte Carlo Methods

The term Monte Carlo was coined in the 1940s by physicists working on nuclear weapon projects. The name Monte Carlo was popularized by physics researchers such as S. Ulam and N. Metropolis (Ulam and Metropolis (1949)). Monte Carlo refers to a casino in Monaco. As the gambling games in casinos, randomness plays an important role in the Monte Carlo method. MC methods rely on repeated sampling to approximate the expectation of simulated random variables. MC methods have a wide application in many different fields, for example in finance and in statistical physics. Monte Carlo methods are often used when it is impos-sible or too complicated to compute results with an analytical method. Next we explain this method in more mathematical terms.

Let X be a random variable (vector) and X has a probability density function (PDF) f . Let h be a function on X. Suppose we want to compute the expectation of h(X). We denote this expectation by θ, where

θ = E[h(X)] = Z

x

h(x)f (x)dx, (2.1)

if X is continuous. E[·] denotes the expectation. If the integral in (2.1) is too complicated to solve analytically, the MC method can be used to obtain a numerical solution to this problem. This can be done in few steps:

1. Generate a set of n independent random samples from the distribution of f : x1, . . . , xn.

2. Calculate h(xi), for i = 1, . . . , n.

3. Calculate the sample mean of h(x1), . . . , h(xn) to obtain an approximation of θ.

4. The MC estimator of θ is give by ˆ θ(n) = 1 n n X i=1 h(xi). (2.2) ˆ

(8)

2 VARIANCE REDUCTION TECHNIQUES E[ˆθ(n)] = E " 1 n n X i=1 h(Xi) # = 1 n n X i=1 E[h(Xi)] = E[h(X)] = θ.

The accuracy of the MC estimator depends on the size of n, the number of simulation repli-cations. Moreover, if θ exists, the law of large numbers tells us that

ˆ

θ(n)→ θ.as

The law of large numbers has a clear intuitive interpretation: as n is large enough, the MC estimator ˆθ(n) shall be close to θ. The efficiency of MC simulations depends on the variance of the simulation. The variance of the MC estimator is given by

Var[ˆθ(n)] = Var " 1 n n X i=1 h(Xi) # = Var[h(X)] n . (2.3)

We cannot calculate Var[ˆθ(n)], since Var[h(X)] is unknown in equation (2.3). Instead, we approximate the variance of h(X) by an unbiased estimator ˆVar[h(X)], where

ˆ Var[h(X)] = 1 n − 1 n X i=1 h(Xi) − ˆθ(n) 2 . (2.4)

Combining the results from equations (2.3) and (2.4), we have an unbiased estimator for the variance of ˆθ(n). This is given by

ˆ Var[ˆθ(n)] = Var[h(X)]ˆ n = 1 n(n − 1) n X i=1 h(Xi) − ˆθ(n) 2 .

By the Central Limit Theorem1_{, it is easy to show that}

√ n(ˆθ(n) − θ) q Var[ˆθ(n)] d → N (0, 1).

This means that if we want to reduce the simulation standard deviation σ_θ(n)ˆ (

q

Var[ˆθ(n)]) by a factor 10, we need to increase the simulation replications n by a factor 100. In other

1

(9)

dis-words, if we can use techniques that reduce σ_θ(n)ˆ by a factor 10, we can reduce the simulation

replications n by a factor 100 to obtain the same level of precision. This is an important motivation behind the VRTs. In the next section we discuss a useful VRT called importance sampling.

2.3 Importance Sampling

2.3.1 Introduction

Before we turn to the mathematical framework of importance sampling (IS), we first discuss some problems which cannot be solved effectively by MC simulations. We will do this by means of some examples.

Example I: Consider a problem of estimating θ := P (X ≥ 5), where X ∼ N (0, 1). We could use the MC method to estimate θ. The naive estimator is

ˆ θ =

Pn

i=1I{Xi≥5}

n , where X1, . . . , Xn are sampled from N (0, 1).

Here, I is the indicator function, and n is the number of simulation replications. For this problem, to approximate θ with a reasonable degree of accuracy, n is required to be very large. On average we need about 3.5 million scenarios to get one non-zero value of I, because the probability that a random sample X is larger than 5 is about _35000001 . It is not appro-priate to set n to such large number. In most cases, because of the limited simulation time, a much smaller number of simulation scenarios will be chosen. This will almost result in an MC estimator ˆθ = 0. The variance of this estimator is zero as well. This gives a confidence interval [L, U ] = [0, 0] (L=lower bound, U=upper bound).

One might ask, is it not enough to know that θ is very close to zero? Does it matter if we approximate θ by 0 or for example 0.01? For many problems, estimating θ by zero is not a big issue if the true value of θ is 0.01. However, for some other problems it is very important to know θ with a greater level of accuracy. For example, suppose we want to use an MC simulation to price a deep out-of-the-money option (see also Example II). In most cases, the price of this option will be very small. It is very important for option traders to know the precise option price, because very small differences in prices are magnified by the fact that traders often buy thousands or even millions of options at a time.

Example II: Consider a European call option where the payoff at time T equals

max(0, ST − K). Here, ST is the future price of the underlying asset (for example stocks)

at time T and K is the strike price. Let S0 be the current price of the underlying asset. A

call option is deep out-of-the-money if S0 is relatively far below the strike price K. Suppose

we want to use MC simulation to estimate the price of this option. This can be done by first generating a set of scenarios of the future stock price ST. Then we calculate for each

scenario the payoff max(0, ST − K). Finally, the option price is estimated by the average of

the present value of the payoffs. Many scenarios of the simulation will end up with a payoff of zero, because S0 is far below the strike price K. This means in order to estimate the option

(10)

2 VARIANCE REDUCTION TECHNIQUES

Can we use fewer simulation replications to solve the problems discussed in examples I and II? The technique importance sampling (IS) provides an alternative solution to these kind of problems. The basic idea behind this technique is to replace the original sampling distribution in the MC simulation by another suitable sampling distribution, ensuring that ‘important’ values are sampled more frequently. For the first example, we could choose the sampling distribution N (5, 1) to replace the original distribution N (0, 1) to obtain more ‘important’ outcomes: values larger or equal to 5. Because we have changed the distribution, the out-comes in IS must be corrected by some kind of weight. Section 2.3.2 discusses what this weight is and why we need to correct the outcomes by this weight. In Section 2.3.3, a method called Maximum Principle (MP) discusses how to select an effective IS distribution. For the second example, we could change the underlying distribution of ST to get more payoffs larger

than zero. In Section 2.3.4, we will tell more about the application of IS to price derivatives. In the next section we present the mathematical framework, by which we will make the idea of IS more concrete. We will show that the estimator based on the technique IS is unbiased if we correct the ‘important’ outcomes by a weight.

2.3.2 Mathematical Framework

Consider the problem of estimating θ (defined in Equation (2.1)), i.e.

θ = E[h(X)] = Z

x

h(x)f (x)dx,

where X is a random element of Rd with probability density function (PDF) f . Let h be a function from Rd to R. The naive Monte Carlo estimator (see Equation (2.2)) discussed in previous section is ˆ θ(n) = 1 n n X i=1 h(Xi), (2.5)

where X1, . . . , Xn are independent drawings from the probability density function f . Let f∗

be another PDF on Rd _satisfying

f (x) > 0 ⇒ f∗(x) > 0, for all x ∈ Rd. Then we can alternatively represent θ by

θ = Z x h(x)f (x) f∗(x)f ∗_(x)dx. _(2.6)

(11)

Here, Ef∗[·] indicates that the expectation is taken with X distributed according to the PDF

f∗. The weight _ff (X)∗_(X) is the likelihood ratio or Radon-Nikodym derivative of the original

distribution with respect to the importance sampling distribution. Next we show how to use importance sampling to estimate θ. The steps are similar to the naive MC simulation:

1. Generate a set of n independent random samples from the distribution of f∗: x1, . . . , xn.

2. Calculate h(xi), for i = 1, . . . , n.

3. Calculate f (xi)

f∗_(x

i) for i = 1, . . . , n.

4. Take the sample mean of h(x1)_ff (x∗_(x1)

1), . . . , h(xn)

f (xn)

f∗_(x

n) to obtain an estimate of θ. We

denote ˆθis(n) as the IS (simulation) estimator of θ. This is given by

ˆ θis(n) = 1 n n X i=1 h(xi) f (xi) f∗_(x i) .

The IS estimator ˆθis(n) is an unbiased estimator of θ. Since

Ef∗[ˆθ_is(n)] = E_f∗ " 1 n n X i=1 h(Xi) f (Xi) f∗_(X i) # = 1 n n X i=1 Ef∗ h(Xi) f (Xi) f∗(Xi) = Ef∗ h(X)f (X) f∗(X) = θ. The variance of this estimator is given by

Var[ˆθis(n)] = Var " 1 n n X i=1 h(Xi) f (Xi) f∗_(X i) # = Varhh(X)_ff (X)∗_(X) i n . (2.7)

We cannot calculate the variance of the IS estimator Var[ˆθis(n)], because Var

h

h(X)f (Xi)

f∗_(X i)

i

is unknown in Equation (2.7). Instead, we approximate the variance of h(X)_ff (X)∗_(X) by an

unbiased estimator ˆ Var h(X)f (X) f∗_(X) = 1 n − 1 n X i=1 h(Xi) f (Xi) f∗_(X i) − ˆθis(n) 2 . (2.8)

Combining equations (2.7) and (2.8), we have an unbiased estimator for the variance of ˆθis(n).

(12)

2 VARIANCE REDUCTION TECHNIQUES ˆ Var[ˆθis(n)] = ˆ Varhh(X)_ff (X)∗_(X) i n = 1 n(n − 1) n X i=1 h(Xi) f (Xi) f∗(Xi) − ˆθis(n) 2 .

To compare the variance between the IS estimator ˆθis(n) in Equation (2.7) and the MC

estimator ˆθ(n) in Equation (2.3), it is sufficient to compare the second moment of these estimators, because ˆθis(n) and ˆθ(n) are unbiased estimators of θ. The second moment of the

IS estimator ˆθis(n) is Ef∗ " h(X) f (X) f∗(X) 2# = Z x h(x)f (x) f∗(x) 2 f∗(x)dx = Z x h(X)2 f (x) f∗_(x)f (x)dx = Ef h(X)2 f (X) f∗_(X) .

Here, Ef[·] indicates that the expectation is taken with X distributed according to the PDF

f . Compared to the second moment of the MC estimator ˆθ(n), which is Ef[h(X)2], the

vari-ance of the IS estimator could be larger or smaller than the varivari-ance of the MC estimator. Depending on the choice of the PDF f∗. We call f∗ the PDF of the importance sampling distribution. The art of importance sampling is to design a good sampling distribution func-tion f∗. Next, we discuss a method called Maximum Principle (MP). This method can help us to design a good importance sampling distribution function.

2.3.3 Maximum Principle Method

An approach to achieve maximum variance reduction in an importance sampling simulation is to try to minimize the variance of the IS estimator ˆθis(n) in (2.7) (see Glasserman,

Hei-delberger, and Shahabuddin (1999)). It is sufficient to set simulation replications n = 1, since the problem of minimizing Var[ˆθis(n)] =

Var h

h(X)_{f ∗(X)}f (X) i

n is the same as the problem of

minimizing Var[h(X)_ff (X)∗_(X)]. The variance of ˆθis(1) is given by

Var[ˆθis(1)] = Ef∗ " h(X)f (X) f∗_(X) 2# − θ2 = Z x h(x)2 f (x) f∗(x)f (x)dx − θ 2

Here, Ef∗ is defined as before, i.e., the expectation that is taken with X distributed according

(13)

of h(X), then we have a zero variance estimator, since in this case we have Var[ˆθis(1)] = Z x h(x)2 f (x) f∗_(x)f (x)dx − θ 2 = Z x h(x)2 θ h(x)f (x)dx − θ 2 = θ Z x h(x)f (x)dx − θ2 = θ2− θ2 = 0.

A zero variance estimator means that we would only need one sample to estimate θ and this sample is equal to θ with probably one. To get a zero variance estimator, _ff (X)∗_(X) needs to be

equal to _h(X)θ . Hence the optimal choice for the importance sampling PDF f∗(X) is f∗(X) = h(X)f (X)

θ .

To get a zero variance estimator, we need to divide h(X)f (X) by the unknown θ. Of course, this cannot be done, since this unknown θ is what we try to estimate here. Nevertheless, this observation provide us some useful insight. Namely, to design an effective sampling strategy, we should try to choose the IS density function in proportion to the product of h(X) and f (X). If we could choose f∗(X) similar to h(X)f (X) we might expect to get a large variance reduction. To be more precise, we could try to choose f∗(X) to have a ‘similar’ shape as h(X)f (X). In practice, it is often good enough to choose the distribution of f∗ to be the same family of distribution as f (X). We could also try to choose f∗(X) so that f∗(X) and h(X)f (X) both take on their maximum value at the same value, say x∗. When we choose f∗(X) this way, we are using the Maximum Principle (MP) method 2_{. In short, to apply this}

method we need to:

• Solve the optimization problem:

x∗ = arg max

x h(x)f (x).

• Choose an IS PDF f∗(X) that has a ‘similar’ shape to h(X)f (X). • The function f∗_{(X) should take on its maximum value at the value x}∗_.

In the next section we will discuss how to apply IS to price financial derivatives. In particular, we show how to choose an efficient IS distribution function by using Maximum Principle.

(14)

2.3.4 Importance Sampling and Options

The variance reduction technique importance sampling is used in many fields. Usually, this technique is used in an MC simulation to estimate rare event probabilities and expectations. This section mainly focuses on the application of importance sampling in financial derivatives pricing. Most of the theories about financial derivative pricing are drawn from Hull (2001), Wilmott (1998) and Luenberger (1998).

The application of importance sampling to derivatives pricing has been studied in many pa-pers. IS appears to be very useful to price deep out-of-the-money options. Reider (1993) shows how to apply IS to decrease the variance in simulations for deep out-of-the-money European options. In Glasserman et al. (1999), IS is combined with another VRT called stratified sampling to reduce simulation variance in pricing Asian options. In Su and Fu (2002) infinitesimal perturbation analysis (IPA) and stochastic approximation (SA) are used to obtain an optimal importance sampling measure to price Asian options.

Next, we explain how to apply the Maximum Principle method to obtain an effective IS sam-pling distribution to price financial derivatives.

Formulation and Settings

We use the famous Black-Scholes model for options valuation. We assume there are two assets in the market: the risk-free asset and the risky asset S (for example stocks). Let the price of S follows a standard Brownian Motion (BM) with drift µ and volatility σ:

dSt= µStdt + σStdWt,

where W is the standard BM. It can be show (see Glasserman (2004)) that the price of the risky asset S at time t is given by

St= S0e(µ−σ

2_/2)t+σW t_,

where S0 is the current price of the risky asset S, and St is its price at time t. Let ω be

the sample path of the underlying stochastic process (for example stock price process). Let C(T, ω) be a European financial claim (claims that can only be exercised at expiration date), and S is the underlying asset. Then the price of this claim at time t = 0 is given by

C0 = e−rTEQ[C(T, ω)],

where r is the constant continuous risk-free interest rate. T is the expiration date and EQ[·] denotes the expectation under a risk-neutral probability measure Q. A risk-neutral probability measure is a probability measure that assumes that the present expected value of all financial assets is equal to the future payoff of the asset discounted at the risk-free interest rate. In the financial market, there are many different financial claims C(T, ω). Some examples of payoff functions C(T, ω) are:

(15)

• European Put: max(0, K − S_T(ω)) • Arithmetic Asian call: max(0,

Pm

i=1SiT /m(ω)

m − K) where m is the number of averaging

periods.

• Geometric Asian call: max(0, (Qm

i=1SiT /m(ω))

1

m − K) where m is the number of

aver-aging periods.

• Continuous Asian call: max(0, T−1RT

0 St(ω)dt − K).

Here, K is the strike price. To estimate the option price C0 we could use the naive MC

simulation described in Section 2.2 . Just generate a sample path ω, calculate the payoff function C(T, ω). Repeat this process n times, then calculate the sample mean to price C0

(see Equation (2.2), and consider the function h as a payoff function). However, to price the deep out-of-the-money options (see example II in Section 2.3.1), naive MC simulations are not appropriate to use. To price these options, importance sampling is usually a better alternative, since this technique samples more ‘important’ values. Here, ‘important’ values are positive payoffs: C(T, ω) > 0.

To apply IS we need to replace the original sampling distribution based on the risk-neutral probability measure Q by an importance sampling distribution. Let f be the PDF of the original distribution and f∗ be the PDF of the IS distribution. Let

f (x) > 0 ⇒ f∗(x) > 0, for all x ∈ Rd. Then we can express C0 by

C0= e−rTEf∗ C(T, ω)f (X) f∗_(X) ,

The Maximum Principle method (see Section 2.3.3) will be used to get an effective IS distri-bution. We need to solve the following first:

• Solve the optimization problem:

x∗(ω) = arg max

x(ω)(C(T, ω)f (x(ω)),

where C(T, ω) is defined as before, and f (X) is a joint PDF. X = (X1, . . . , Xn), and

the Xi’s are IID N (0, 1).

• Choose an importance sampling joint PDF f∗(X) that has a ‘similar’ shape to C(T, ω)f (X). • The PDF f∗(X) should take on its maximum value at the value x∗ = (x∗₁, . . . , x∗_n).

The first problem is easy to solve by using a solver (usually the solver in Microsoft Excel is enough). To solve the second problem we first choose f∗(X) to be a joint PDF, where X = (X1, . . . , Xn), and Xi ∼ N (µi, 1) for i = 1, . . . , n, Xi’s are independent. This PDF

(16)

choose µi = x∗i for i = 1, . . . , n.

The IS distribution obtained by using Maximum Principle is a joint PDF f (X), where X = (X1, . . . , Xn), and Xi ∼ N (x∗i, 1), Xi’s are independent. By choosing our IS

distri-bution this way, we put the most density weight to the sample path which has the most influence on the price of an option. A path ω that maximizes the product of payoff C(T, ω) and payoff probability f (X(ω)) is the most important path.

In Chapter 3 and Chapter 4 we will apply the Maximum Principe method to obtain IS distribution to price European call and Asian call options. This method is also used to price the value of a minimum return guarantee (MRG), which can be seen as a put option. We will discuss this subject in more details in Chapter 5. In the next section we describe some other well-known VRTs.

2.4 Other Variance Reduction Techniques

In Section 2.3 we have discussed a VRT called importance sampling in detail. In this section we briefly describe two other well-known variance reduction techniques: Antithetic Variable (AV) and Control Variate (CV). Among these techniques there are also many other VRTs, for example stratified sampling, moment matching methods and Latin hypercube sampling. We will not discuss these techniques here. More details about these techniques can be found in Glasserman (2004). This book also discusses the techniques IS, AV and CV.

Antithetic Variable

Antithetic variable (AV) is a very simple VRT. The idea is to reduce simulation variance by introducing negative dependence between pairs of replications. Suppose a simulation is driven by m independent standard normal variables Zi, . . . , Zm. The AV method generates an

antithetic pair for each variable Zi, . . . , Zm. The antithetic pair of Zi is −Zi for i = 1, . . . , m.

Let h be a function of Zi, . . . , Zm. To estimate h(Z), n samples of Zi, . . . , Zm are generated.

For each sample, the value of h(Z) and h(−Z) will be calculated. The sample value is equal to h(Z)+h(−Z)₂ . This works well because if h(Z) is below or above the true value, it will be corrected by the value of its antithetic pair h(−Z). This method is widely used because of its simplicity and effectiveness. However, in some cases, this method is not suitable to apply. In Chapter 3 and 4, we will show that this method will give incorrect price for options which are deep out-of-the-money.

Control Variate

The control variate (CV) method has been applied successfully in derivatives pricing (see for example Hull and White (1987)), and can be described as follows. Suppose there are two similar derivatives: D1 and D2. D1 cannot be calculated analytically, whereas D2 has an

explicit formula to calculate its value. Suppose we want to estimate the value of the derivative D1. Then we can use CV to estimate the value of D1. This method works as follows:

(17)

2. Use naive MC simulations to estimate the value of D1 and D2. Let ˆVD1 and ˆVD2 be the

MC estimator for respectively D1 and D2.

3. Let VD2 be the value from the explicit formula for D2. Then the CV estimator for D1

is ˆθ = ˆVD1+ α(VD2− ˆVD2), Ross (1997) shows that Var[ˆθ] is minimized when

α = Cov(D1, D2) Var[D2]

The idea behind this technique is that the simulation error in ˆVD1 is about the same as the

simulation error in ˆVD2, which we know to be VD2 − ˆVD2. By adding VD2 − ˆVD2 to ˆVD1 we

(18)

3 Pricing European Options

3.1 Introduction

In this thesis the variance reduction technique (VRT) importance sampling (IS) will be used to price portfolio of unit-linked insurance liabilities. To get familiar with this technique, we start with an application to a plain vanilla European option. An important reason of choos-ing this type of option is that the value of European options can be calculated analytically. Hence we can test the performance of importance sampling by comparing its estimator to the closed-form solution. A closed-form solution for European options can be obtained by using the Black-Scholes options pricing formula. This famous formula is published in Black and Scholes (1973).

The structure of this chapter is as follows: an introduction to European options is given in Section 3.2. In Section 3.3 we show how to obtain an IS estimator for European call options. Beside this estimator, we also present two other estimators for European call options: the naive Monte Carlo (MC) estimator and the antithetic variable (AV) estimator. To compare the performance of importance sampling to the other two methods, the sample error of the estimators will be calculated and compared to each other. The simulation results are presented in Section 3.4. At the end this chapter, we draw conclusions about the application of IS for European option pricing.

3.2 European Options

An option is a contract to buy or sell an asset at some point in the future, for example a stock. There are two types of options: call options and put options. You can buy or sell either type. A call option gives the option holder the right to buy the underlying asset at a predetermined price. A put option gives the option holder the right to sell the underlying asset at a predetermined price. Furthermore, the buying price or selling price of the underly-ing asset is specified in a contract. The price in a contract is officially known as the exercise price or strike price. An option can only be exercised at the date specified in a contract. The date is known as the expiration date or maturity. European options can only be exercised on the expiration date. There are also options which can be exercised at any time up to the expiration date. These options are known as American options.

We will illustrate the possible payouts from options by means of two examples. In the first example, we show the profit from buying a European call option. Consider an investor who buys a European call option with the following characteristics:

• The option holder has the right to purchase 100 shares of stock A • The price of the option to buy one share is 1 euro

• The expiration date T is in two year • The strike price K is 10 euro

(19)

The option can only be exercised at the expiration date, because this is a European option. Ignoring transactions costs and dividend, the investor will only exercise the option if the stock price at the expiration date is higher than the strike price. Suppose for example that the stock price at the expiration date is 12 euro. By exercising the option, the investor makes a gain of 200 euro (= (12 − 10) × 100). If we take the initial cost of the option into account, the net profit is 100 euro (= 200 − 100). If the stock price at expiration date is lower than the strike price, there is no point to exercise the option, because buying a stock from the market is cheaper. Figure 1 shows the net profit per share from buying this call option.

Figure 1: Net profit from buying a European call option. The strike price K is 10 euro and the initial investment is 1 euro (option price).

The net profit function (per share) of European call options can be formulated as

max(0, ST − K) − I, where ST is the stock price at the expiration date T , K is the strike

price, and I is the initial cost.

An investor who buys a call option is expecting that the stock price will increase, while a buyer of a put option is expecting that the option price will decrease. Consider now an investor who buys a European put option with the following characteristics:

• The option holder has right to sell 100 shares of stock B • The price of option to sell one share is 1 euro

• The expiration date is in two year • The strike price is 8 euro

(20)

3 PRICING EUROPEAN OPTIONS

Ignoring transactions costs, the investor will only exercise the option if the stock price at the expiration date is lower than the strike price. Suppose for example that the stock price at the expiration date is 5 euro. The investor can buy 100 share of stock B from the market for 5 euro per share. By exercising the put option, the investor makes a gain of 300 euro. Taking the initial cost of the option into account, the net profit is 200 euro. Figure 2 shows the net profit from buying this put option.

Figure 2: Net profit from buying a European put option. The strike price K is 8 euro and the initial investment is 1 euro.

The net profit function of European call options can be formulated as max(0, K − ST) − I,

where ST, K and I are defined as above.

Formulation and Settings

A famous model for options valuation is the Black-Scholes model. This model assumes that St, the price of the underlying asset of an option, follows a standard Brownian Motion (BM)

with drift µ en volatility σ:

dSt= µStdt + σStdWt,

where W is a standard BM. The price of the risky asset S at time t (see Glasserman (2004)) is given by

St= S0e(µ−σ

2_/2)t+σW

t_, _(3.1)

where S0 is the current asset price. The price of a European call option at t = 0 is

C0 = e−rTEQ[max(0, ST − K)], (3.2)

(21)

the expectation under a risk-neutral probability measure3 Q. The Black-Scholes closed-form solution for C0 is C0(S0, T ) = S0N (d1) − Ke−rTN (d2), (3.3) where d1= ln(S/K) + (r + σ2/2)T σ√T , and d2 = d1− σ √ T .

N (·) in Equation (3.3) is the cumulative probability distribution function for a standardized normal distribution. In the next section, we describe three simulation methods to price European call options: importance sampling, naive Monte Carlo, and antithetic variable. The estimators obtained from these methods will be compared to the price from Black-Scholes (BS) closed-form solution, we call this price the BS option price. Furthermore, the performance of the methods will be measured by comparing the sample error of these estimators to each other.

3.3 Simulation Estimators for European Call Options

This section shows how to use simulation techniques to estimate C0, the price of a European

call option. C0 is defined as in Equation (3.2), i.e., C0 = e−rTEQ[max(0, ST − K)]. By using

the formula in Equation (3.1) we can express ST, the asset price at time T by

ST = S0e(µ−σ

2_{/2)T +σ}√_{T X}

, (3.4)

where X is a random sample from the standard normal distribution N (0, 1). The methods naive Monte Carlo (MC) (Section 2.2), importance sampling (IS) (Section 2.3) and antithetic variable (AV) (Section 2.4) are used to obtain estimators for C0. We denote the estimators

from these methods by ˆθmc, ˆθis and ˆθav respectively. In Chapter 2 we have shown that these

estimators are unbiased.

Naive Monte Carlo Simulation

The naive MC estimator for the price of a European call option C0 in Equation (3.2) can be

calculated as follows:

1. Generate n independent random scenarios by sampling X1, . . . , Xn from the standard

normal distribution N (0, 1).

2. Calculate ST, the stochastic asset price at the expiration date T for each scenario. ST

can be calculated by using the formula in Equation (3.4), where X in (3.4) is a random sample from step 1. We denote by ST,i the stock price at expiration date T for scenario

i, for i = 1, . . . , n.

3_{A risk-neutral probability measure is a probability measure that assumes the present expected value of all}

(22)

3. Calculate the payoff max(0, ST,i− K) for i = 1, . . . , n.

4. The naive MC estimator ˆθmc for a European call option is the present value of the

sample mean of the payoffs. This estimator is unbiased and is given by

ˆ θmc= 1 n n X i=1

max(0, ST,i− K)e−rT.

The sample error of this estimator is given by

S_θˆ_mc = v u u t 1 n(n − 1) n X i=1

(max(0, ST,i− K)e−rT − ˆθmc)2.

Importance Sampling

To obtain an IS distribution to price European call options we use the Maximum Principle method. This method is described in Section 2.3. According to this method, to obtain a ‘good’ sampling distribution we need to solve the following:

1. Maximize the product of the option payoff h(X) and the probability of this payoff f (X). For a European call option, the payoff h is defined by max(0, ST − K), where ST is the

asset price at the expiration date T and K is the strike price.

2. Choose an importance sampling PDF f∗(X), which has a ‘similar’ shape to h(X)f (X). 3. The PDF of f∗(X) should take on its maximum value at the value x∗, where x∗ is a

solution from step 1.

First, we maximize the product of the payoff h(X) and the probability of this payoff f (X) as follows:

µ∗ = arg max

x h(x)f (x). (3.5)

The function h is max(0, ST − K), the payoff function of a European call option, where ST

is the asset price at the expiration date T and K is the strike price. The function f is the probability density function (PDF) of the original sampling distribution, in this case it is a standard normal distribution N (0, 1). The optimization problem above can be solved by a suitable software package. In this thesis the ‘Microsoft Excel Solver’ is used to solve this problem.

(23)

Finally, we choose the mean µ and the variance σ2 for the IS distribution. Let f∗ be the PDF of an IS distribution. According to Maximum Principle, f∗(·) should take on its maximum value at the value µ∗. Here, µ∗ is a solution to the maximization problem above. We know that a normal distribution function takes on its maximum value at mean µ, therefore, we choose the mean of the IS distribution to be µ∗. A simple choice for the variance of the IS distribution is to leave the original variance unchanged. Now the IS distribution becomes a normal distribution with the mean µ∗, and the variance σ2 = 1. In the next section we will investigate if this method is appropriate to use for obtaining a ‘good’ IS distribution. We will do this by means of simulation experiments.

Now we have chosen an IS sampling distribution to replace the original distribution. Next, we show how to obtain an IS estimator for European call options:

1. Generate n independent random scenarios by sampling Y1, . . . , Yn from the IS

distribu-tion N (µ∗, 1).

2. Calculate f (Yi) and f∗(Yi) for i = 1, . . . , n.

3. Calculate ST,i, the asset price at time T , for scenario i = 1, . . . , n by using the formula

in Equation (3.4) (replace X in (3.4) by random sample Y in Step 1). Note that ST,i

is not the asset price under the risk-neutral probability measure Q. Since the standard normal distributed variable X in (3.4) is replaced by the variable Yi, which follows a

normal distribution N (µ∗, 1). By doing this, we actually changed the standard BM Wt

in Equation (3.1) by another BM with the drift µ∗. 4. Calculate the payoff max(0, ST,i− K)_ff (Y∗_(Yi)

i), for i = 1, . . . , n.

f (Yi)

f∗_(Y

i) is the importance

sampling weight for scenario i.

5. The IS estimator ˆθisfor a European call option is the present value of the sample mean

of the payoffs. This unbiased estimator is given by ˆ θis= 1 n n X i=1 max(0, ST,i− K) f(Yi) f∗(Yi) e−rT.

S_θˆ_is = v u u t 1 n(n − 1) n X i=1 max(0, ST,i− K) f(Yi) f∗(Yi) e−rT − ˆθis 2 . Antithetic Variable

The concept of antithetic variable is explained in Section 2.4. The idea of this method is to reduce the simulation variance by introduce negative dependence between pairs of replications. For each random sample generated in a simulation, we also calculate an antithetic pair for this sample. The antithetic pair is defined by −X if X is the random sample generated in a simulation. Next, we present a procedure to obtain the AV estimator ˆθav for the option price

(24)

1. Generate n₂ independent random scenarios by sampling X₁+, . . . , X+n 2

from the standard normal distribution N (0, 1), where n is an even number.

2. Calculate X₁−, . . . , X−n 2

, X_i− is equal to −X_i+ for i = 1, . . . ,n₂. We call X₁−, . . . , X−n 2

the antithetic pairs of X₁+, . . . , X+n

2

.

3. Calculate S_T,i+ and S_T,i− for scenario i = 1, . . . , n by using the formula in Equation (3.4). S_T,i+ is the asset price at the expiration date T based on the random sample X_i+ from step 1, and S_T,i− is the asset price at T based on X_i− from step 2.

4. Calculate the payoffs max(0, S_T,i+ − K) and max(0, S_T,i− − K). 5. Define Zi, the average of max(0, ST+− K) and max(0, S

+ T − K) by Zi = 1 2 max(0, S + T − K) + max(0, S − T − K) , for i = 1, . . . , n 2. (3.6) The AV estimator ˆθav for a European call option is the present value of the average of

Zi’s. This unbiased estimator is given by

ˆ θav = 1 n/2 n/2 X i=1 Zie−rT,

where Zi is defined as in Equation (3.6), for i = 1, . . . ,n₂.

The sample error of the AV estimator ˆθav is given by

S_θˆ_av = v u u t 1 n/2(n/2 − 1) n/2 X i=1 (Zie−rT − ˆθav)2. 3.4 Simulation Results

In this section the methods presented in the previous section are used to estimate European call option prices. The simulation results are generated by using the software package ‘Mi-crosoft Excel’. First we present the simulation results for the European call options and discuss the performance of IS. Later on this section, we will perform sensitivity analysis on the choice of the IS distribution. The main objective of this section is to answer the following questions:

1. What is the performance of IS compared to the naive MC and AV methods in European call options pricing?

2. What is the effect of moneyness on the variance reduction (VR) by using importance sampling instead of the naive MC and AV methods?

(25)

4. Is the Maximum Principle method an appropriate method to choose IS distribution in European call option pricing?

Table 1 shows the naive MC, AV and IS estimators for the European call options with the strike price K = 20, 50, 80, 100, 120, 150, 200 and 250. The sample error of these estimators and the VR are also given in this table. For example, the VR by using IS instead of naive MC is defined as

S_θmcˆ

S_θisˆ

2

, where S_θˆ_mc is the sample error of a naive MC estimator and S_θˆ_is

is the sample error of an IS estimator. The VR given in Table 1 is based on the exact sample errors from the simulation experiments. This could differ from the VR based on the rounded sample error given in Table 1. The VR can be interpreted as the speed up factor, meaning that, for example, if the VR by using the IS method instead of the naive MC method is 10, then in order to achieve the same level of precision of the estimator, we would need only one tenth the number of simulation replications for the IS simulation, when compared to the MC method. In other words, under the condition that the sample error of the estimators are the same, the IS simulation is 10 times faster than the naive MC simulation if the VR is 10.

Table 1: Comparison of the methods: European call options with strike price K = 20, 50, 80, 100, 120, 150, 200 and 250. r is the risk-free interest rate. µ and σ are the drift and volatility of the underlying asset. S0 is the current asset price. T is the expiration date. ‘BS price’ is

the Black-Scholes closed-form solution. ˆθmc, ˆθav, and ˆθis are respectively the naive MC, AV

and IS estimators based on 1000 replications.

r = µ = 0.04, S0 = 100, σ = 0.1, T = 10 year and simulation replications n = 1000

Estimators Sample errors VR

K BS price ˆ θmc θˆav θˆis S_θˆ_mc S_θˆ_av S_θˆ is S_θmcˆ S_θisˆ 2 S_θavˆ S_θisˆ 2 20 86.59 85.51 86.34 86.63 1.009 0.280 0.033 954 73 50 66.49 65.40 66.23 66.62 1.009 0.280 0.109 86 7 80 46.59 45.45 46.27 46.78 1.006 0.309 0.236 18 2 100 34.23 33.05 33.88 34.39 0.962 0.388 0.285 11 2 120 23.67 22.51 23.38 23.72 0.878 0.491 0.285 10 3 150 12.33 11.69 11.97 12.35 0.688 0.543 0.208 11 7 200 3.47 3.18 3.11 3.45 0.388 0.327 0.082 22 16 250 0.87 0.87 0.61 0.86 0.213 0.159 0.025 72 40

(26)

better than the naive MC and AV methods, especially if the options are deep in-the-money or deep out-of-the money (see Figure 4).

Figure 3: Deviation naive MC, AV and IS estimators from the BS price for different strike prices K. µ = 0.04, σ = 0.1, S0 = 100, T = 10 year and n = 1000 replications.

(27)

If an option is more out-of-the-money, less samples in the naive MC simulation will end up with the payoff larger than zero. Because of this, the sample error of the naive MC estimator will be higher if an option is more out-of-the-money. By using IS to shift the original sampling distribution to the region where the payoff is larger than zero, more simulation samples will end up with the payoff larger than zero. By doing this, the sample error of IS estimator increased much slower than the sample error from the other methods (see Figure 5). This is the reason why VR increases if the options are more out-of-the-money.

Figure 5: Relative sample error (sample error / estimator value) for different strike prices K. µ = 0.04, σ = 0.1, S0 = 100 and T = 10 year. Number of simulation replications n = 1000.

(28)

larger than zero. If we increase the volatility, more payoff larger than zero will be generated, which reduces the sample error. However, if the volatility is too extreme (for example σ = 1), payoffs larger than zero could also have a very large variation. So even we have more payoffs larger than zero, the sample error will still increase. The effect of the moneyness on the IS method and the AV methods differs in some cases from the naive MC method. For example, the extremely large volatility has no negative impact on the relative sample error of the IS method. The reason for this is that the payoffs are corrected by the IS weight.

Table 2: Comparison of the methods: deep in-the-money European call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

r = µ = 0.04, S0 = 100, K = 20, T = 10 year and simulation replications n = 1000

σ BS price ˆ θmc θˆav θˆis S_θˆ_mc S_θˆ_av S_θˆ_is S_θmcˆ S_θisˆ 2 S_θavˆ S_θisˆ 2 0.05 86.59 86.08 86.53 86.60 0.493 0.070 0.008 3742 76 0.1 86.59 85.51 86.34 86.63 1.009 0.280 0.033 954 73 0.2 86.60 84.32 85.49 86.77 2.193 1.094 0.136 259 64 0.5 88.86 85.10 78.51 89.26 9.614 6.059 0.492 381 151 1 96.36 80.23 38.33 96.76 34.369 10.135 0.402 7308 636

Table 3 and 4 show the results for the at-the-money and deep out-of-the-money options with σ = 0.05, 0.1, 0.2, 0.5 and 1. Similar to the results from Table 2, we see that the IS method is much better than the naive and AV methods, especially for the large σ’s.

Table 3: Comparison of the methods: at-the-money European call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

(29)

Table 4: Comparison of the methods: deep out-of-the-money European call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

σ BS price ˆ θmc θˆav θˆis S_θˆ_mc S_θˆ_av S_θˆ_is S_θmcˆ S_θisˆ 2 S_θavˆ S_θisˆ 2 0.05 0.00 0.00 0.00 0.00 0.000 0.000 0.000 0 0 0.1 0.87 0.87 0.61 0.86 0.213 0.159 0.025 72 40 0.2 9.28 8.46 1.20 9.22 7.927 0.956 0.210 1420 21 0.5 45.75 43.54 8.92 45.75 34.942 5.542 0.783 1989 50 1 85.39 70.83 34.05 85.50 28.196 9.708 0.874 1040 123

The simulation results show that the variance reduction technique IS performs much better than the naive MC and the AV methods for the European call options. For all options considered in this section, the IS estimators are close to the BS option price and have relatively small sample error. Finally, we summarize the answers to the questions formulated at the beginning of this section:

1. For almost all tested cases (excluding one case with the BS price almost equal to zero), the IS estimators deviate less than 1.5% from the BS option price. These estimators perform significantly better than the naive MC and the AV estimators, because for the deep out-of-the-money options or the options with large volatility, the naive MC estimators deviate in some cases more than 15%. For the AV estimators, the deviation is in some cases even more than the MC estimators, up to 85% for the deep out-of-the-money options. When we look at the VR, the IS method is signifincantly better than the other two method. Ignoring the case where the BS price almost equal to zero, the VR by using the IS method instead of the naive MC method varies from 10 to 7308 and for the VR by using the IS method instead of the AV method varies from 0.5 to 636. 2. The VR increases if options are more out-of-the-money or more in-the-money.

3. The VR differs for the different types of moneyness. In most cases, large volatility has negative impact on the relative sample error (the sample error divided by the estimator value).

4. From the simulation results, we can conclude that the Maximum Principle method is an appropriate method to obtain an IS distribution.

Sensitivity Analysis

(30)

• What is the effect on the IS estimator by changing the IS drift µ∗_{? Here, µ}∗ _{is a}

solution to the optimization problem (3.5) in previous section. This is the mean of the IS distribution.

• What is the effect on the sample error of the IS estimator by changing the IS drift µ? First we show what happens to the IS drift µ∗ for different strike price K. Figure 6 shows that the IS drift µ∗ increases if the strike price K increases. This is logical, because if an option is more out-of-the-money we need to shift the IS distribution more to the right so that more samples will end up with payoff larger than zero.

Figure 6: Relationship between moneyness and importance sampling drift µ∗

Figure 7 shows what happens to the IS estimators for the deep out-of-the-money, at-the-money and deep in-the-at-the-money options if we deviate the IS drift from µ∗. A deviation from the IS drift µ∗ means that we increase or decrease µ∗ by a fixed percentage. The volatility σ here is equal to 0.1. The result for other volatilities (σ = 0.2 and σ = 0.5) can be found in Appendix 1a. We see that for all cases the IS estimator based on the IS drift µ∗ is near the BS price, which indicates again that the Maximum Principle method works well. If we change µ∗ by 50%, the deviation from BSs price is less than 6%.

(31)

Figure 7: Effect of deviation from IS drift µ∗ on IS estimator. The underlying asset volatility σ = 0.1. Strike price K = 250, 100 and 20 for the deep out-of-the-money, at-the-money and deep in-the-money options respectively.

Figure 8: Effect of deviation from IS drift µ∗ on the sample error of the IS estimator. The underlying asset volatility σ = 0.1. Strike price K = 250, 100 and 20 for the deep out-of-the-money, at-the-money and deep in-the-money options respectively.

3.5 Conclusions

(32)

compared to the naive and antithetic varies methods. From the simulation results we con-clude that for all tested cases (deep in-the-money, at-the-money and deep out-of-the money options combined with the different volatilities), the IS performs much better than the other methods. In general we achieve the most variance reduction for the deep in-the-money and deep out-of-the-money options. Furthermore, for all tested cases, the IS estimators are close to the BS prices, while the naive MC and the AV estimators significantly differ from the BS prices for the deep out-of-the-money options and the options with large volatility.

(33)

4.1 Introduction

In the previous chapter we have used the variance reduction technique importance sampling (IS) to price European call options. We continue to use this technique to price more complex financial options, known as exotic options. In this section we primarily focus on one type of exotic options called Asian options. An introduction to this option is given in Section 4.2. In Section 4.3, three procedures are given to obtain estimators for Asian call options. These procedures are based on the simulation techniques importance sampling, naive Monte Carlo (MC) and antithetic variable (AV). In Section 4.4 the simulation results are presented. Finally, in Section 4.5 we draw conclusions about the performance of importance sampling in Asian option pricing.

4.2 Asian Options

Asian options are path-dependent options, which means that the payoff of this type of option depends on the historical price process of underlying assets (for example stocks). In this section we mainly focus on cases in which the payoff depends on S(t1), . . . , S(tm) at a fixed

set of dates t1, . . . , tm. Let S(ti) be the price of underlying asset at time ti, for i = 1, . . . , m.

In this case, it is possible to obtain an unbiased simulation estimate of the option price (see Section 4.3). There are two basic types of Asian options: arithmetic and geometric. The payoffs of arithmetic Asian call and put options are max( ¯S − K) and max(K − ¯S) respectively, where K is the strike price. ¯S is defined by

¯ S = 1 m m X t=1 S(ti). (4.1)

For geometric Asian options, the arithmetic average ¯S in (4.1) is replaced by

m Y t=1 S(ti) !_m1 .

The average ¯S above depends on the discrete path S(t1), . . . , S(tm). Instead of discrete

averaging, an Asian option could also depend on the entire path S(t), 0 ≤ t ≤ T over an interval [0, T ]. Usually, T is the expiration date of an option. We call the Asian options which depend on the entire path S(t), 0 ≤ t ≤ T the continuous averaging Asian options. The payoffs for such a continuous averaging Asian option are max( ¯S − K) and max(K − ¯S) for respectively call and put options. For an averaging period [u, t], ¯S is given by

¯ S = 1 t − u Z t u S(τ )dτ. (4.2)

(34)

4 PRICING ASIAN OPTIONS

(1992)). In this chapter we focus on pricing arithmetic discrete averaging Asian options by means of simulation methods.

The price of arithmetic discrete averaging Asian call options is given by

Ca = e−rTEQ[max(0, ¯S − K)], (4.3)

where Q is a risk-neutral probability measure, r is the risk-free interest rate, T is the expiration date, and K is the strike price. ¯S is defined as in Equation 4.1. We assume the price process of the underlying assets is modeled by Black-Scholes model (see Section 2.3.4). It can be shown (see for example Glasserman (2004))that the asset price at time ti is given by

S(ti) = S(ti−1)e(µ−σ

2_/2)(t

i−ti−1)+σ

√

ti−ti−1Xi_{, for i = 1, . . . , m.} _(4.4)

In this formula, µ is the drift of the underlying asset, σ is the volatility of the underlying asset and Xi is a random variable, which has a standard normal distribution N (0, 1). In the next

section, the procedure for calculating the IS estimator for arithmetic discrete averaging Asian call options is given. We will compare this estimator to the naive MC and AV estimators.

4.3 Simulation Estimators for Asian Call Options

This section discusses three procedures for obtaining simulation estimators for the price of an arithmetic discrete averaging Asian call option Ca. Which is defined as in Equation (4.3)

from the previous section. The procedures are based on the simulation techniques: naive Monte Carlo (MC) (Section 2.2), importance sampling (IS) (Section 2.3) and antithetic vari-able (AV) (Section 2.4). The estimators based on these methods are denoted by ˆθmc, ˆθis and

ˆ

θav respectively. In Chapter 2 we have shown that these estimators are unbiased.

Naive Monte Carlo Simulation

A procedure for obtaining the naive MC estimator for Ca, the price of an arithmetic Asian

call is as follows:

1. Generate X ∈ N (0, 1)m×n, this is a m × n random matrix and each entry in X is a sample from distribution N (0, 1). Here, m is the number of averaging periods, n is the number of simulation scenarios.

2. Calculate the stochastic asset price path S(t1), . . . , S(tm) at a fixed set of dates t1, . . . , tm

for each scenario. The path S(t1), . . . , S(tm) can be calculated by using random matrix

in step 1 and the formula in Equation (4.4). We denote Sj(ti), as the asset price at time

ti for scenario j, for i = 1, . . . , m and j = 1, . . . , n.

3. Calculate the payoff max(0, ¯Sj− K) for j = 1, . . . , n, where the average asset price ¯Sj

is defined as in Equation (4.1) for scenario j.

4. The naive MC estimator ˆθmc for an arithmetic Asian call is the present value of the

(35)

The sample error of this estimator is given by S_θˆ_mc = v u u t 1 n(n − 1) n X j=1 max 0, Pm i=1Sj(ti) m − K e−rT − ˆθmc 2 . Importance Sampling

Next, we will turn on the problem of pricing an arithmetic Asian call by means of the impor-tance sampling method. Again the key question is how to choose a ‘good’ sampling distri-bution? In Section 2 we have seen that the Maximum Principle method can be successfully applied to obtain an IS distribution for pricing European call options. For Asian options, we also apply this method to find an IS sampling distribution. The procedure is similar to that for the European options. We first solve the optimization problem:

µ∗ = arg max

x h(x)f (x).

Here, the function h is max(0, ¯S − K), the payoff function of an arithmetic Asian call. The Averaging price ¯S is defined as before and K is the strike price. The function f (·) is the joint probability density function (PDF) of the original sampling distribution. In this case, it is a multivariate standard normal distribution with the covariance matrix equal to the identity matrix. The ‘Microsoft Excel Solver’ is used to solve the optimization problem above. Note that the solution µ∗ for the problem above is a vector, and the size of this vector is equal to the number of averaging periods m. For the European options, µ∗ is just a number. This is because Asian options are path dependent and the European options are not. That why we need a m sized multivariate standard normal distribution for the Asian options to generate random paths, where for the European options we only need a standard normal distribution to generate random asset price at the expiration date.

After that, we need to choose a distribution family for the IS distribution. The PDF of this distribution should have a ‘similar’ shape to h(·)f (·). We could choose a multivariate normal distribution as the IS distribution. This solution could be ‘good’ since f (·) is also a joint PDF of a multivariate standard normal distribution.

Finally, we have to choose the mean and the variance of the IS distribution. Let f∗ be the joint PDF of an IS distribution. According to Maximum Principle, f∗(·) should take on its maximum value at µ∗, where µ∗ is a solution to the maximization problem above. Since a multivariate normal distribution function takes on its maximum value at mean vector µ, we choose µ = µ∗. As for the European options we leave the original variance of the normal distributions unchanged. Now the IS distribution becomes a multivariate normal distribution with mean µ∗ and covariance matrix equal to the m × m identity matrix Im, where m is the

number of averaging periods.

Next, we show how to obtain an IS estimator for arithmetic Asian calls:

1. Generate n independent random scenarios by sampling Y1, . . . , Yn from the IS

distribu-tion N (µ∗, Im). The Yj’s are random vector of size equal to the number of averaging

(36)

2. Calculate f (Yj) and f∗(Yj) for j = 1, . . . , n. f and f∗ are the PDF’s of the original

sampling distribution and the IS distribution respectively.

3. Calculate the stochastic asset price path Sj(t1), . . . , Sj(tm) at a fixed set of dates

t1, . . . , tm for scenario j = 1, . . . , n. Sj(ti) is the asset price at time ti for scenario

j, for i = 1, . . . , m and j = 1, . . . , n. The path Sj(t1), . . . , Sj(tm) can be calculated by

using formula in Equation (4.4) (replacing Xi in (4.4) by Yi,j from step 1 for scenario

j = 1, . . . , n).

4. Calculate the payoff max(0, ¯Sj− K)_ff (Y∗_(Yj)

j) for scenario j = 1, . . . , n, where the average

asset price ¯Sj is defined as in Equation (4.1) for scenario j.

5. The IS estimator ˆθisfor an arithmetic Asian call is the present value of the sample mean

of the payoffs. This estimator is unbiased and is given by ˆ θis= 1 n n X j=1 max(0, Pm i=1Sj(ti) m − K) f (Yj) f∗_(Y j) e−rT

S_θˆ_is = v u u t 1 n(n − 1) n X j=1 max(0, Pm i=1Sj(ti) m − K) f (Yj) f∗_(Y j) e−rT _{− ˆ}_θ is 2 . Antithetic Variable

Similar to the naive MC method, an AV estimator ˆθav for the option price Cacan be obtained

as follows:

1. Generate X+∈ N (0, 1)m×n₂_{, this is a m ×}n

2 random matrix, each entry in matrix X is

a sample from the N (0, 1) distribution. Here, m is the number of averaging periods, n₂ is the number of simulation scenarios, and n is an even number.

2. Calculate X−, X−is a m ×n₂ matrix defined by −X+. This means X_i,j−, an entry in X− is equal to -X_i,j+, for i = 1, . . . , m and j = 1, . . . , n. We call X_i,j−, . . . , X_m,j− the antithetic path of X_i,j+, . . . , X_m,j+ for scenario j.

3. Calculate the stochastic asset price path S+(t1), . . . , S+(tm) at a fixed set of dates

t1, . . . , tm for each scenario, where S+(t1), . . . , S+(tm) is the price path based on the

random sample from step 1. Calculate also the price path S−(t1), . . . , S−(tm), which

is based on the antithetic path X_i,j−, . . . , X_m,j− from step 2. S+(t1), . . . , S+(tm) and

S−(t1), . . . , S−(tm) can be calculated by using the formula in Equation 4.4.

4. Calculate the payoffs max(0, ¯S_j+− K) and max(0, ¯S_j−− K) for scenario j = 1, . . . , n. ¯S_j+ and ¯S_j−are defined as in Equation (4.1) by replacing the asset price path S(t1), . . . , S(tm)

by S_j+(t1), . . . , Sj+(tm) and Sj−(t1), . . . , Sj−(tm) respectively. Define Zj, the average of

max(0, ¯S_j+− K) and max(0, ¯S_j−− K) by Zj =

1

(37)

5. The AV estimator ˆθav for an arithmetic Asian call is the present value of the average of

the Zj’s. The unbiased estimator is given by

ˆ θav = 1 n/2 n/2 X j=1 Zje−rT,

where Zj is defined as in Equation (4.5).

The sample error of ˆθav is given by

S_θˆ_av = v u u u t 1 n/2(n/2 − 1) n/2 X j=1 (Zje−rT − ˆθav)2. 4.4 Simulation Results

This section presents the simulation results of pricing arithmetic Asian call options by using the procedure discussed in the previous section. The software package ‘Microsoft Excel’ is again used to generate simulation results. The interval between the averaging dates t1, . . . , tm

is chosen to be a constant. The performance of the technique importance sampling (IS) will be compared to the techniques naive Monte Carlo (MC) and antithetic variable (AV). After that, we will perform sensitivity analysis on the choice of the IS distribution. The main goal of this section is to answer the following questions:

1. What is the performance of importance sampling compared to the naive MC and AV methods in pricing an arithmetic Asian call?

2. What is the effect of moneyness on the performance of importance sampling?

3. What is the effect of stock volatility on the performance of importance sampling. IS the effect different for at-the-money, deep in-the-money and deep out-of-the-money options? 4. Is Maximum Principle a good method to obtain an IS distribution for pricing arithmetic

Asian call?

To judge the quality of the estimators we need the ‘correct’ Asian option prices to compare with. Unfortunately, there exits no explicit formula to price arithmetic Asian options (that why we want to estimate it with simulations). However, we can use naive simulations with many replications to obtain ‘reliable’ option prices. 50000 replications is chosen here, because from the simulation results we see that the relative sample error (sample error divided by the estimator) for a naive simulation with 50000 is very ‘low’ (see Table 5, 6, 7 and 8).

Table 5 shows the naive MC, AV and IS estimations of the Asian call option prices. We consider a range of strike prices: K = 20, 40, 60, 80, 100, 140, 160, and 200. The estimators are based on the simulations with 1000 replications.

(38)

Table 5: Comparison of the methods: Asian call options with strike price K = 20, 40, 60, 80, 100, 120, 140, 160 and 200. r is risk-free interest rate. µ and σ are the drift and the volatility of the underlying asset. S0 is the current assets price. T is the expiration date.

m is the number of averaging periods. The interval between the averaging periods are the same. ‘Naive 50k’ is the naive MC estimator based on 50000 replications. ˆθmc, ˆθav, and ˆθis

are respectively the naive MC, AV and IS estimators based on 1000 replications.

r = µ = 0.04, S0 = 100, σ = 0.1, m = 10, T = 10 year and simulation replications n = 1000

K Naive 50k ˆ θmc θˆav θˆis Naive 50k S_θˆ_mc S_θˆ_av S_θˆ is S_θmcˆ S_θisˆ 2 S_θavˆ Sˆ_θis 2 20 70.71 70.40 70.54 70.66 0.078 0.553 0.132 0.016 1170 66 60 43.89 43.58 43.73 43.90 0.078 0.553 0.132 0.070 63 4 80 30.54 30.23 30.38 30.59 0.078 0.549 0.146 0.132 17 1 100 18.07 17.70 17.80 18.15 0.072 0.514 0.216 0.161 10 2 120 8.68 8.34 8.21 8.73 0.057 0.410 0.296 0.132 10 5 140 3.44 3.36 3.16 3.46 0.038 0.279 0.244 0.072 15 11 160 1.19 1.22 1.07 1.19 0.023 0.173 0.158 0.030 33 27 200 0.10 0.17 0.14 0.11 0.007 0.064 0.057 0.004 331 260

estimators lose their precision. When we look at the IS estimators, they are almost all equal to the estimators from the naive simulations with 50000 replications. The sample errors of the estimators and the variance reduction are also given in Table 5. We see that in all cases IS estimators have a much smaller sample error than the other two methods, and the vari-ance reduction (IS vs Naive MC) starts from 10 (for at-the-money option) and increases to 1170 (for deep in-the-money option). For deep out-of-the-money options we also observe huge variance reduction. For example, for K = 200, the VR is 331. This result is similar to the results from European call options (see Chapter 3): the VR increases if an option is more out-of-the-money or more in-the-money. The VR are based on the exact sample errors from the simulations. The VR based on the sample errors given in the table might be different from the ‘real’ VR. Figure 9 shows the VR for different strike price K.

(39)

Figure 9: Variance Reduction by using IS technique instead of naive MC for Asian options with different strike price K. µ = 0.04, σ = 0.1, S0= 100, m = 10 and T = 10, and n = 1000

replications.

Table 6: Comparison of the methods: deep in-of-the-money Asian call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

r = µ = 0.04, S0 = 100, K = 20, m = 10, T = 10 and simulation replications n = 1000

(40)

Table 7: Comparison of the methods: at-the-money Asian call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

r = µ = 0.04, S0= 100, K = 100, m = 10, T = 10 and simulation replications = 1000

σ Naive 50k ˆ θmc θˆav θˆis Naive 50k S_θˆ_mc S_θˆ_av S_θˆ_is S_θmcˆ S_θisˆ 2 S_θavˆ S_θisˆ 2 0.05 17.10 16.95 17.05 17.13 0.038 0.269 0.043 0.077 12 0.3 0.1 18.07 17.70 17.80 18.25 0.072 0.514 0.216 0.165 10 2 0.2 22.44 21.79 21.51 22.63 0.142 1.039 0.735 0.280 14 7 0.5 37.80 38.95 36.26 37.47 0.531 4.310 3.817 0.540 64 50 1 53.27 65.26 59.76 57.15 2.654 17.452 15.392 1.228 202 157

Table 8: Comparison of the methods: deep out-of-the-money Asian call options with volatility σ = 0.05, 0.1, 0.2, 0.5 and 1.

r = µ = 0.04, S0= 100, K = 160, m = 10, T = 10 and simulation replications = 1000

σ Naive 50k ˆ θmc θˆav θˆis Naive 50k S_θˆ_mc S_θˆ_av S_θˆ_is S_θmcˆ S_θisˆ 2 S_θavˆ S_θisˆ 2 0.05 0.03 0.05 0.04 0.03 0.002 0.022 0.019 0.001 378 296 0.1 1.19 1.22 1.07 1.21 0.023 0.173 0.158 0.030 33 27 0.2 6.74 6.79 6.18 6.92 0.091 0.691 0.616 0.142 14 19 0.5 26.66 28.63 25.75 26.96 0.495 4.074 3.663 0.479 72 58 1 47.69 59.74 54.53 50.83 2.641 17.327 15.282 1.144 230 179

We summarize the answers to the questions at the begin of this section:

1. In almost all tested cases, IS estimators based on 1000 replications are almost equal to the naive MC simulation with 50000 replications. In almost all cases, using the IS method instead of the naive MC or AV methods reduces sample variance significantly. 2. The performance of the IS method increases if the options are more out-of-the-money

or more in-the-money.

3. The sample variance of the IS estimator increases if the underlying asset becomes more volatile.

Variance Reduction Technique: Importance Sampling Pricing minimum return guarantee in Unit-linked life insurance

Pricing minimum return guarantee in Unit-linked life insurance

Preface

Contents

3

Pricing European Options