Variance Reduction for Asian Options

Technical report, IDE0737 , February 21, 2008

Variance Reduction for Asian Options

Master’s Thesis in Financial Mathematics Galina Galda

School of Information Science, Computer and Electrical Engineering Halmstad University

Variance Reduction for Asian Options

Galina Galda

Halmstad University Project Report IDE0737

Master’s thesis in Financial Mathematics, 15 ECTS credits Supervisor: Prof. Bernard Lapeyre

Examiner: Prof. Ljudmila A. Bordag External referee: Prof. Krzysztof

February 21, 2008

Department of Mathematics, Physics and Electrical Engineering School of Information Science, Computer and Electrical Engineering

Halmstad University

Preface

This thesis would not exist without the support from my supervisor prof. Bernard Lapeyre;

the head of the programme prof. Ljudmila A. Bordag;

and a special thank is directed to my friend Alexander Stromilo.

Abstract

Asian options are an important family of derivative contracts with a wide variety of applications in commodity, currency, energy, interest rate, equity and insurance markets. In this master’s thesis, we investigate methods for evaluating the price of the Asian call options with a fixed strike. One of them is the Monte Carlo method. The accuracy of this method can be observed through variance of the price. We will see that the variance with using Monte Carlo method has to be decreased.

The Variance Reduction techniques is useful for this aim. We will give evidence of the efficiency of one of the Variance Reduction techniques - Control Variate method - in a mathematical context and a numerical comparison with the ordinary Monte Carlo method.

iii

Contents

1 Introduction 1

2 Methods 7

2.1 The Monte Carlo method . . . 7

2.1.1 The Monte Carlo simulation . . . 9

2.1.2 Variance . . . 24

2.2 The Variance Reduction techniques . . . 27

2.2.1 Standard methods of the Variance Reduction . . . 27

2.2.2 Control Variate method . . . 32

3 Results 41 3.1 Comparison . . . 41

3.2 Reduction of Variance . . . 43

4 Conclusions 47

Notation 49

Bibliography 51

Appendix 55

v

Chapter 1 Introduction

A financial derivative is a contract which provides to the holder a future payment that depends on the price of the assets such as stocks, currencies, commodities. In a frictionless market, the no-arbitrage principles allow to express the value of the derivative as an expectation of its discounted future payment. Options are particular derivatives characterized by non negative payoffs. European-style options can be exercised only at the expiration date, in contrast to American-style options, where the holder can exercise them earlier the maturity date.[1]

The analytic formulas are available for the fair price of the ordinary Euro- pean call and put options written on a stock whose price is modelled as a Geometric Brownian motion. For more complicated derivatives, the analytic formulas do not easy solve. Such derivatives are usually priced by Monte Carlo simulation or by numerical methods. An important class of options with such problem is the class of Asian options.

Asian options are options in which the underlying variable is the average price over a period of time. In other words, the payoff of such options depends on the average price of the asset during a specified period leading up to the ma- turity date. Asian options are often used for protection against unexpected changes of prices. Average-value options have a lower volatility and they are less expensive than regular options. These options are commonly traded on currencies, energy, interest rate, equity, insurance markets and commodity products which have low trading volumes. The name “Asian” option was emerged in 1987 when the Banker’s Trust Tokyo office used it for pricing average options on crude oil contracts. Asian options can be classified into three categories: arithmetic average Asians, geometric average Asians and both these forms are averaged on a weighted mean basis, whereby a given weight is applied to each stock being averaged.[1]

1

The payoff of geometric Asian options is given as:

fcall = max ((

Yn i=1

Si)ⁿ¹ − K, 0)

f_put = max (K − ( Yn i=1

S_i)ⁿ¹, 0) The payoff of arithmetic Asian options is given as:

f_call = max (1 n

Xn i=1

S_i − K, 0)

fput = max (K − 1 n

Xn i=1

Si, 0)

The average value for the arithmetic case is the sum of the sampled asset prices divided by the number of samples:

AvgA = S1+ S2 + ... + Sn

n and for the geometric case:

Avg_G= pⁿ

S₁S₂...S_n

where the n-th root of the sample values multiplied by each other is taken.

The payoff functions for Asian options are given with a fixed strike

f = max (0, η(S_A− K)) and with a floating strike

f = max (0, η(S_t− S_A))

where η is a binary variable which is set to 1 for a call, and −1 for a put.[13]

European-style Asian (Eurasian) options can be exercised at the expiration date only, whereas American-style(Amerasian) ones offer earlier exercise pos- sibilities, which may become attractive when the current asset price is below the current running average (i.e. is pulling down the average) for a call op- tion, and when it is above the running average for a put. Techniques for pric- ing Amerasian options were developed by Barraquand and Pudet (1996)[2], Grant, Vora, and Weeks (1997) [3]. Hull and White (1993) [4] have adapted

VRAS 3 binomial trees (from the model of Cox, Ross, and Rubinstein 1979) to com- pute the price of Amerasian options. Broadie and Glasserman (1997) [5]

have offered a simulation method based on non-recombining trees in the lat- tice model and which produces two estimators of the option value, one with a negative bias and one with a positive bias. Zvan, Forsyth, and Vetzal (1998) [6] have produced stable numerical PDE techniques, adapted from the field of computational fluid dynamics, for pricing American style Asian options with continuously sampled prices. H. B. Ameur, M. Breton and P. L’Ecuyer [1]

have developed a numerical method for pricing American-style Asian options based on dynamic programming combined with finite - element piecewise - polynomial approximation of the value function.[1]

The following we offer some of the models to price European-style Asian options under a variety of methods.

1) Kemna and Vorst, 1990

Kemna and Vorst (1990) [7] developed a closed form pricing solution to ge- ometric averaging options by altering the volatility and cost of carry term.

The geometric average of the underlying prices follows a lognormal distribu- tion, whereas with the arithmetic average rate options this condition falls.

Hence geometric averaging options can be priced via a closed form analytic solution.

The solutions to the geometric averaging Asian call and puts are given as:

call = S_T0e^aN(d₁) − Ke^cN(d₂) and

put = Ke^cN(−d₂) − S_T0e^aN(−d₁)

where N is the cumulative standard normal distribution function a = (b − r)(T − T₀)

c = −r(T − T₀)

where r is the risk-free rate, T - the maturity date and T0-the initial date.

d₁ = ln S_{T 0}/K + (b − 0.5σ_A²)T σ_A√

T

d₂ = ln S_{T 0}/K + (b + 0.5σ_A²)T σA

√T

The adjusted volatility and dividend yield are given as:

σ_A= σ

√3

b = 1

2(r − D − σ² 6 )

where σ is the observed volatility, D is the dividend yield.

We intend to scrutinize this method later.[13]

2) Turnbull and Wakeman (1991)

Turnbull and Wakeman (1991) [8] have suggested that the distribution under arithmetic averaging is approximately lognormal, and they put forward the first and second moments of the average in order to price the option.

The analytical approximations for a call and a put according to TW are given as

call ≈ Se^(b−r)T2N(d1) − Ke^−rT2N(d2) and

put ≈ Ke^−rT2N(−d2) − Se^(D−r)T2N(−d1)

d₁ = ln S K + (b + 0.5σ²_A)T2

σ_A√ T₂ d2 = d1− σA−p

T2

where T₂ is the time remaining until maturity. If the average period has begun, then T₂ is the original time to maturity T , otherwise T₂ is T − τ . The adjusted volatility and dividend yield:

σ_A =

rln M₂ T − 2b b = ln M₁

T

The first and second moments (assuming that the averaging period has not yet begun):

M1 = e^(r−D)T − e^(r−D)τ (r − D)(T − τ ) M₂ = 2e^{(2(r−D)+σ2)T}S²

(r − D + σ²)(2r − 2q + σ²)T2

+ 2S² (r − D)T²

µ 1

2(r − D) + σ2 − e^(r−D)T r − D + σ²

¶

VRAS 5

If the averaging period has already begun, we have to adjust the strike price KA= T

T₂K −(T − T2) T₂ SAvg

where T - the original time to maturity, T₂ - the remaining time to maturity, K - the original strike price, S_Avg - the average asset price.[13]

3) Binomial Method and Trinomial Trees (Hull and White(1993)) Asian options can be priced using lattice/tree methods. At any point in time on the tree, the value of the option is dependent upon the average of the price that the path has taken. We determine a minimum and maximum range at each node. As the number of nodes on a tree grows, so does the number of averages which must be taken, particularly in the central nodes - this is because the number of averages to be taken is exponentially related to the number of possible asset paths. This is the problem, which Hull and White (1993) attempted to solve by adding a state variable to the tree nodes, and approximation is undertaken with interpolation techniques in backward induction.[13]

4) Finite Differences Method (Rogers and Shi(1995))

Rogers and Shi (1995)[9] have presented a method using a one-dimensional

PDE ∂ V

∂ t(x, t) + 1

2σ²x²∂²V

∂ x²(x, t) − (ρ(t) + rx)∂ V

∂ x(x, t) = 0 with limit condition and price for the fixed strike

V (x, T ) = x P rice = S_T0V ( K S_T0, 0) and for the floating strike

V (x, T ) = (1 + x) P rice = S_T0V (0, 0)

This PDE can be solved using finite differences. This method has problems concerning the diffusion term, particularly with lower volatilities and short times to maturity. [13], [16]

5) Other Methods

Levy [10] has developed analytical approximation. Curran [11] has put for- ward an approximation of an arithmetic Asian option based on a geometric conditioning approach. Geman and Yor (1993)[12] have used Bessel processes and have derived exact formulas for the Laplace transform of the value of a continuous-time European-style Asian option. Shaw (2002) [14] and Linetsky (2002) [15] also has used Laplace transforms to price Asian options and has reached a reasonable degree of accuracy. [13]

6)Monte Carlo Simulation

Various methods using the Monte Carlo simulation have been offered to price Asian options. The above-mentioned analytical approximations by Turnbull and Wakeman, Curran, Levy can be computed using a simulation method.

The Monte Carlo simulation can give sufficiently accurate prices for option values, and in the case of Asian options, which are highly path-dependent, this method is particularly useful.

We need to observe that the Monte Carlo simulations are not really effective in pricing Asian options. However, if we use the Variance Reduction tech- niques (in our case a Control Variate method), the accuracy of the Monte Carlo simulation may be improved - for example by using the closed form geometric average rate formula by Kemna and Vorst as a control variate. [13]

We will prove it in this thesis.

This thesis is organized in the following way: we first introduce the ordi- nary Monte Carlo method (2.1) with the approximation of the integral of St

by using the time schemes [16], and the calculation of the price of the Asian options by using the simulations(2.1.1). In the Subsection (2.1.2) we present the simulation of the variance, which will be reduced in the second Section (2.2). This Section will also illustrate the several Variance Reduction tech- niques with an aspect to the Control Variate. This method will be intently scrutinized at the Subsection (2.2.2). Then, in the Chapter 3 we present a comparison between two methods of the pricing Asian options - the ordinary Monte Carlo method and the Control Variate method. In the appendix some theorems, definitions and the text of a program are given.

Chapter 2 Methods

2.1 The Monte Carlo method

In this section of the thesis we would like to present the Monte Carlo method with application to the pricing of the Asian options. The historic part is very interesting. Lets look into the history.

“The term “Monte Carlo” was introduced by von Neumann and Ulam during World War II as a code word for the secret work at Los Alamos; it was sug- gested by the gambling casinos in the city of Monte Carlo in Monaco. The Monte Carlo method was then applied to problems related to the atomic bomb. The work involved direct simulation of behavior concerning neutron diffusion in fissionable material. Shortly thereafter Monte Carlo methods were used to evaluate complex multidimensional integrals and to solve cer- tain integral equations occurring in physics, which were not amenable to analytic solution.

The Monte Carlo method can be used not only for the solution of determin- istic problems. A deterministic problem can be solved by the Monte Carlo method if it has the same formal expression as some stochastic processes.

Another field of application of the Monte Carlo methods is the sampling of a random variate from probability distributions...

The Monte Carlo method is now the most powerful and commonly used tech- nique for analyzing complex problems. Applications can be found in many fields from radiation transport to river basin modelling. Recently, the range of applications has been broadening, and the complexity and computational effort required has been increasing, because realism is associated with more complex and extensive problem descriptions.”[17]

In a simple case, the essence of the Monte Carlo method consists of the 7

following: we need to find value a of some investigated value. For this aim we choose a random variable X, which has an expectation equaled to a:

E(X) = a. In practice, we have to make n simulations which give in result n possible values of X, then we will find the arithmetic average x = ^P_n^xi and accept x as an estimation ea of the required value a: a = ea = x.

This method is also called stochastic simulation. The theory of the Monte Carlo method shows how we can choose the most effective random variable X and how we can find the possible value of X.

For instance we have an integral of a bounded real valued function (not necessary to be a smooth function)

I = Z

[0,1]^d

f (x) dx

Let represent I as E(f (U)), where U is a uniformly distributed random vari- able on [0, 1]^d. The Strong Law of Large Numbers [see Appendix (application 1)] shows that the average

S_N = 1 N

XN i=1

f (U_i) (2.1)

converges to E(f (U)) almost surely when N → 0 and (U_i, i ≥ 1) is a family of uniformly distributed independent random variables on [0, 1]^d. So, the ap- proximation of integral I needs to call a random number generator N times and to compute the average S_N = _N¹ P_N

i=1f (U_i).

We will seek a rate of convergence of the Monte Carlo method to determine when it is more efficient than deterministic algorithms.

Obviously from the Central Limit Theorem [see Appendix application1], the convergence rate of the Monte Carlo method is rather slow (1/√

N). The approximation error is random and, even when N is large, may take large values. However, the Monte-Carlo methods are useful in practice.

The deterministic methods require many points of calculation and the num- ber of points increases exponentially with the dimension of the space, whereas the Monte Carlo methods require the simulation of independent random vec- tors (X1, X2, ...Xd) with independent coordinates. This compares to the computation of the one-dimensional situation, and the number of trials is multiplied by d, that is why the method remains tractable even if d is large.

The Monte Carlo methods also have another advantage, which is their par- allel nature. Each processor of a parallel computer can be assigned to the task of making a random trial.

Thus, the Monte Carlo methods are used in situations where the determinis- tic methods are inefficient, especially when the dimension of the state space

VRAS 9 is vary large.[18]

The Monte Carlo methods can be useful in the financial area to compute prices of options, to estimate sensitiveness of portfolios to various parame- ters and to compute risk measurements.

2.1.1 The Monte Carlo simulation

Here we explain the Monte Carlo simulation with respect to our derivative European style Asian call option with the fixed strike.

The Monte Carlo simulation works as follows:

1. Sample a random path for St (Geometric Brownian Motion) under a risk- neutral probability.

2. Calculate the payoff from the derivative.

3. Repeat steps 1 and 2 to get many sample values of the payoff from the derivative.

4. Compute the mean of the sample payoffs to get an estimate of the expected payoff.

5. To get an estimate of the value of the derivative, we need to discount the expected payoff at the risk-free rate.

The next two subsections explain this in detail.

Generating Sample Paths

We use the assumptions that there are no transaction costs or taxes; all securities are perfectly divisible; there are no dividends during the life of the derivative; there are no riskiness arbitrage opportunities; security trading is continuous; the risk-free rate of interest, r, is constant and the same for all maturities.

We will describe the price of an asset at time t. For this aim we use the Black and Scholes model with a risk asset (price is S_t at time t) and no-risk asset (price is S_t⁰). The price S_t changes to S_t+ dS_t by a small subsequent time interval dt. The return of the asset dS/S is decomposed into two parts.

The first is a deterministic, predictable and anticipated return - µ dt, where µ is a drift (a measure of the average of growth of the asset prices).

The second is the random change in the asset price in response to an external effect - σ dB, where σ is the volatility, which measure the standard deviation

or the returns.[20]

Putting these contributions together, we obtain the stochastic differential equation

dS_t St

= µ dt + σdB_t (2.2)

In our case µ and σ are two positive constants and (B_t) is a standard Brow- nian motion.

Brownian Motion. Brown, a botanist, discovered the motion of pollen particles in the water. Einstein, Perrin and other physicists studied Brow- nian motion at the beginning of the twentieth century. In the simple case, any continuous random process Xtwhich is a martingale and has a quadratic variation equal to t is the Brownian motion. In other words, the Brownian motion is a stochastic process which has stationary independent increments and which has also continuous sample paths. These increments usually have a Normal distribution for a given time increment which makes them useful in finance. Brownian paths start at the origin. For every Brownian increments t, h ≥ 0, Xt+h− Xt is independent of Xu : 0 ≤ u ≤ t and has a Gaussian distribution of mean 0 and variance h. Brownian motion is another name for a Wiener process. Brownian motion belongs to the every interesting class of processes - a martingale, a Gaussian process, a Markov process, a diffusion, a Levy process. Brownian motion offers an extremely simple, but robust model of the stock returns. It is far from being completely realistic, but is a clean and reliable starting point for the further improvement.

We are interested in a standard Brownian motion. If we have the Brown- ian motion X_t, with drift µ and diffusion coefficient σ² (abbreviated X ∼ BM (µ, σ²)), then

Xt− µt σ

is the standard Brownian motion B_t. The complete definition can be found in [see Appendix application 2] or [21].

Thus, we can construct X from the standard Brownian motion B by set- ting X_t = µt + σB_t, which follows X_t ∼ N(µt, σ²t). Besides, X solves the stochastic differential equation dX_t = µdt + σdB_t.

VRAS 11 Brownian bridge. Useful to know that a Brownian bridge is a continuous- time stochastic process whose probability distribution is the conditional prob- ability distribution of the Brownian motion B_t given the condition that B0 = B1 = 0. If {Bt, t ≥ 0} is the standard Brownian motion starting at zero, then the Brownian bridge is the stochastic process {B_t, 0 ≤ t ≤ 1|B₁ = 0}.

A Brownian bridge may be realized as the process {B_t−tB₁, 0 ≤ u ≤ 1}.[23]

The increments in the Brownian bridge are clearly not independent.

Law(B_t 0 ≤ t ≤ 1|B₁ = x) = Law(B_t− tB₁+ tx) Law(B_t0 ≤ t ≤ 1|B₁ = 0) = Law(B_t− tB₁) E

µZ _T

0

S_sds|B_h, B_2h, ..., B_T

¶

= Z _T

0

E (S_s|B_h, B_2h, ..., B_T) ds =

= XN k=1

Z _(k+1)h

kh

E(S_s|B_kh, B_(k+1)h) = E³

e^{(r−σ2/2)s+σBs}|B_kh, B_(k+1)h´ where kh < s < (k + 1)h

Law(B_s|B_kh = x, B_(k+1)h = y) =

= N

µ(k + 1)h − s

h x +s − kh

h y,((k + 1)h − s)(s − kh) h

¶

Lets prove this statement. Accept the notation kh = t_k, (k + 1)h = t_k+1, t_k+1− t_k= h, hence the conditional law of B_s given by:

Law(B_s|B_tk = x, B_tk+1 = y) = N

µt_k+1− s

h x + s − t_k

h y, (t_k+1− s)(s − t_k) h

¶ (2.3) Consider that tk → 0 and tk+1 → tk+1− tk, so B0 = x, Bt = y,

t_(k) ≤ s ≤ t_(k+1)

B_s = α B_tk + β B_tk+1 + Z

(Z, B_tk, B_tk+1) is a gaussian vector, Z independent of B_tk, B_tk+1. Hence,

C(Z, B_tk) = 0, C(Z, B_tk+1) = 0

⇔

E(Z B_tk) = 0, E(Z B_tk+1) = 0

E((Bs− α Btk − β Btk+1) Btk) = 0, E((Bs− α Btk − β Btk+1) Btk+1) = 0 E(B_sB_v) = min (s, v)

t_k− α t_k− β t_k = 0, s − α t_k− β t_k+1 = 0

⇒ β = 1 − α ⇒

u − α t_k− (1 − α) t_k+1= 0

⇒

α₀ = tk+1− s

h β₀ = s − tk

h α tk+ β tk+1 = s − t_k

h tk+1+ t_k+1− s

h tk = s Z0 = Bs− α0Btk− β0Btk+1

We need to know law of Z0.

E(Z₀) = 0

V (Z₀) = V (B_s) + α²₀V (B_tk) + β₀²V (B_tk+1)−

−2α₀C(B_s, B_tk) − 2β₀C(B_s, B_tk+1) + 2α₀β₀C(B_tk, B_tk+1) =

= s +

µtk+1− s h

¶₂ tk+

µs − tk

h

¶₂

tk+1− 2tk+1− s h tk−

−2s − t_k

h s + 2t_k+1− s h

s − t_k h t_k=

= (t_k+1− s)²(s − t_k)² h²

Remember that our task is the finding the law of the standard Brownian bridge:

Law(Bs|Btk = x, Btk+1 = y) = N

µtk+1− s

h x + s − tk

h y, (tk+1− s)(s − tk) h

¶

=

= N(E(B_s),p

V (B_s)) B_s = α B_tk + β B_tk+1 + Z

E(Bs) = E( α Btk + β Btk+1 + Z) = |E(Z) = 0| =

= E( α B_tk+ β B_tk+1) = |B_tk = x, B_tk+1 = y| = α x + β y = t_k+1− s

h x+s − t_k h y pV (B_s) =

q

α B_tk + β B_tk+1 + Z =p

V (Z) =

r(t_k+1− s)²(s − t_k)² h²

pV (B_s) = (t_k+1− s)(s − t_k) h

w.r.p.

VRAS 13 Simulation of the Brownian Motion. To simulate the Brownian motion we need to remember that it is a Gaussian process, subsequently it has density given by

ρ = 1

√2πe^−x22

The Gaussian random variables we can present as follows:

X =p

−2log(U)sin(2πV ) or

Y =p

−2log(U)cos(2πV )

They are independent and have a mean 0, variance 1. U and V are two independent random variables which are uniformly distributed on [0, 1]. [18]

0 0.2 0.4 0.6 0.8 1

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8

time 0 ≤ t ≤ 1 Value of BM: B t

Brownian Motion Sample with 1000 time steps

Figure 2.1: Brownian motion.

Thereupon, the simulation of the Brownian Motion can be represented by Figure(2.1).

The Geometric Brownian Motion. The Geometric Brownian Motion is an exponentiated Brownian Motion. In other words, if log S_t is the Brown- ian motion, then S_t is the stochastic process called the Geometric Brownian Motion. Consequently, all methods for the simulating Geometric Brownian Motion get methods for the simulating Brownian motion, through exponen- tiation.

In a finance, the Geometric Brownian Motion is the most essential model of the asset. Originally, Paul Samuelson used the Geometric Brownian Motion as the financial model in the 1960s. The use of ordinary Brownian Motion in the model of the price of the stock or any other limited liability asset is an undesirable feature, because it can takes the negative values. In contrast with Geometric Brownian Motion which always has the positive values. Lets compare the increments of the both motions. For the Brownian Motion we have:

B_t2 − B_t1, B_t3 − B_t2, ..., B_tk − B_tk−1 are independent for any k and any 0 ≤ t₁ < t₂ < ... < t_k≤ T . The Geometric Brownian Motion has percentage changes

S_t2 − S_t1

S_t1 ,S_t3− S_t2

S_t2 , ...,S_tn− S_tn−1

S_tn−1 (2.4)

are independent for t1 < t2 < ... < tn. It is rather than the absolute changes S_ti − S_ti−1. These fact confirm the use of the Geometric rather than the ordinary Brownian Motion in the modelling asset prices.

Assume Bt is the standard Brownian Motion and Xt corresponds dX_t = µdt + σdB_t

If we accept S_t= S₀e^Xt ≡ f (X_t), then an application of Itˆo’s formula evinces that

dS_t = f⁰(X_t)dX_t+1

2σ²f⁰⁰(X_t)dt

= S₀e^Xt[µdt + σdB_t] +1

2σ²S₀e^Xtdt

= S_t(µ + 1

2σ²)dt + S_tσ dB_t (2.5) In contrast, the Geometric Brownian Motion process is frequently specified through a stochastic differential equation of the form

dS_t St

= µdt + σdB_t (2.6)

VRAS 15 Comparison of the last two equations shows that the models are incongruous and it evinces an uncertainty in the role of µ.

In the formula (2.5) we have µ as a drift of the Brownian motion which we exponentiated to specify St - the drift of log St. St has the drift µSt in the formula (2.6), which can be implied

d logS_t= (µ − 1

2σ²)dt + σdB_t (2.7)

It may be compared with formula (2.5) or verified by Itˆo’s formula.

S is the process of the type in the formula (2.6) with the notation S ∼ GBM(µ, σ²).

In the formula (2.6), we have µ as the drift parameter (not the drift of either S_tor logS_t). From the formula (2.7) we see that if S has the initial value S₀, then

S_t = S₀e^{[µ−12σ2]t+σBt}. To know more information, we refer to [21]

0 0.2 0.4 0.6 0.8 1

40 60 80 100 120 140 160 180 200 220 240

time 0 ≤ t ≤ 1 Value of GBM: S t

Geometrical Brownian Motion Sample with 1000 time steps

Figure 2.2: The realization of the asset price by the Geometric Brownian Motion, in the low volatility case (σ = 0.3).

0 0.2 0.4 0.6 0.8 1 0

100 200 300 400 500 600 700 800

time 0 ≤ t ≤ 1 Value of GBM: St

Geometrical Brownian Motion Sample with 1000 time steps

Figure 2.3: The realization of the asset price by the Geometric Brownian motion, in the high volatility case (σ = 0.9).

A risk-neutral dynamics. The ordinary differential equation d S_t⁰ = r S_t⁰dt

is presented as the price of the no-risk asset with the interest rate r.

Consider that the today’s fair (that is arbitrage-free) price of the derivative is equal to the expected value (under the probability measure) of the future payoff of the derivative discounted at the risk-free rate. This probability measure is called a risk-neutral measure in the mathematical finance, because under this measure, all financial assets have the same expected rate of return, regardless of the variability in the price (a riskiness) of the asset. The value of the derivative can be easily expressed in the formula by risk-neutral measure.

Suppose H_T is the random variable on the probability space describing the market. It can be presented as payoff of the derivative at some time T in the future. As well suppose that the discount factor from time T0 = 0 until time T is P (0, T ). By this assumption, the today’s fair value of the derivative becomes

H0 = P (0, T )EQ(HT)

where Q is the risk-neutral measure. In terms of the physical measure P - i.e.

the actual probability distribution of the prices where (almost universally)

VRAS 17 the more risky assets (those assets with a higher price volatility) have a greater expected rate of return than less risky assets:

H₀ = E_P(dQ dPH_T)

Here ^dQ_dP is the Radon-Nikodym derivative (see [32]) of Q with respect to P . The risk-neutral measure is also called the equivalent martingale mea- sure. A particular financial market may has the more than one risk-neutral measure, which leads to an interval of the arbitrage-free prices (in this case the equivalent martingale measure terminology is more commonly used). If there are just one risk-neutral measure then there are the unique arbitrage- free price for each asset in the market.

Using the Black and Scholes model with the risk asset, remember the stochastic differential equation (2.2)

dS_t= S_t(µdt + σdB_t)

where B_t is the standard Brownian motion with respect to the physical mea- sure P .

Define

Wt = Bt+ µ − r σ t

The Girsanov’s theorem [see Appendix application 3] shows that the measure Q exists under which Wt is the Brownian motion.

In the probability theory, the Girsanov’s theorem evinces how stochastic processes change under changing in the measure. In the financial math- ematics, this theorem shows how to convert the physical measure to the risk-neutral measure. The physical measure describes the probability that the underlying instrument (such as the share price or the interest rate) will take the particular value or values. The risk-neutral measure is a very useful tool for estimating the value of the derivatives on the underlying.[27]

The market price of the risk can be noted as µ − r σ t.

By substitution

W_t = B_t+ µ − r σ t σ W_t = σ B_t+ (µ − r)t σ d W_t = σ d B_t+ (µ − r) d t

σ d Bt = (µ − r) d t + σ d Wt

in

dSt= St(µdt + σdBt) we have

dS_t= S_t(µ dt + (µ − r) d t + σ d W_t) the risky asset satisfies a new stochastic differential equation

dS_t= S_t(r dt + σ d W_t) (2.8) Q is the unique risk-neutral measure for the model.

The solution of this formula (2.8) is

S_t = S_T0e^{(r −σ22 ) t + σ Wt} (2.9) where S_T0 is the price of the asset at the beginning of the modelling.

The more information about the Risk-Neutral probability can be found in [24], [22] or [26].

Calculation of the price of the Asian options

We explain the Monte Carlo simulations with respect to the our derivative - the European-style Asian (Eurasian) call options with the fixed strike.

The Monte Carlo simulations use the risk-neutral valuation result. We want to obtain the expected payoff. So, we compute the following:

e^−rTE(payof f )

Remember, we have the Eurasian call options with the fixed strike which have the payoff equaled to

f (s, a) = (a − K)₊

where a is the average of the prices under some period (T − T₀).

a = 1 T − T₀

Z _T

T0

S_udu

If T₀ = 0 and T - the maturity is given by the price history, the average is computed over the whole life of the option. Thus, the price of the Asian call options with the fixed strike can be written:

V (T, S) = e^−rTE(1 T

Z _T

0

S_udu − K)₊

VRAS 19

For the approximation of the integral Z _T

0

S_udu we can accept time schemes which will be described at the next paragraph.

The time schemes approximation. Using the Monte Carlo method to compute the price of the Asian options, we can exactly simulate S_tinstead the simulation of the average of St, accordingly we do not need to approximate the integral.

Here, we introduce time schemes, which represented by B. Lapeyre and E.

Temam (2000) [16] to estimate Y_T =

Z _T

0

S_u du

The interval [0, T ] will be divided into the N steps with the size h = T /N and we need to define the times t_k= k T /N = k h.

The "Riemann" scheme

The integral Y_T can be approached by the Riemann sums:

Y_T^r,n = h Xn−1 k=0

S_tk (2.10)

For instance, we want to know the price of the Asian put options with the fixed strike. We will use the Monte Carlo method with M number of the drawing. The approximation of the price will be:

V = 1 Me^−rT

XM j=1

(K − h T

Xn−1 k=0

S_tk)₊

This algorithm has the time complexity O(_{N M}¹ ), which involving two kinds of errors: the step error and the Monte Carlo error ^√σ_M. This time complexity is true for the every kind of the Monte Carlo method.

The "Trapezoidal" scheme

This scheme can give us the higher accuracy for the integral approximation.

Note that this scheme is an equivalent to the well known trapezoidal method.

Lets assume that

E((1 T

Z _T

0

S_udu − K)₊| B_h)

is the closest random variable to (1

T Z _T

0

S_udu − K)₊ B_h is the σ-field generated by the (S_tk, k = 0, ...., N ).

We can accept that (E(1

T Z _T

0

S_udu | B_h) − K)₊ = (1 T

Z _T

0

E(S_u| B_h)du − K)₊

it is possible since the conditional law of W_u for u ∈ [t_k, t_k+1] with respect to Bh given by

Law(Wu|Wtk = x, Wtk+1 = y) = N

µtk+1− u

h x + u − tk

h y, (tk+1− u)(u − tk) h

¶

Look at the paragraph “The Brownian Bridge” to prove this law (2.3).

Using the conditional law of W_u, we can get E

·1 T

Z _T

0

S_udu | B_h

¸

= 1 T

Z _T

0

E(S_u| B_h)du

S_u = f (B_h) = e^{(r−σ22 )u+σ Bh}, t_k < u < t_k+1

E(f (Bu) | Btk, Btk+1) = E(e^{(r−σ22 )u+σ Bu}| Btk = x, Btk+1 = y) = 1

T Xn−1 k=0

Z _tk+1

tk

e^{(r−σ22 )}e^{σtk+1−uh Wtk+σu−tkh Wtk+1+σ22} (tk+1−u)(u−tk)

h du =

1 T

Xn−1 k=0

Z _tk+1

tk

e^{ru−σ22 u+σ}

tk+h−u h



W_tk+σ(^u−tkh )^W_tk+1^+σ22h(tk+1u−tk+1tk−u²+utk)

du =

1 T

Xn−1 k=0

Z _tk+1

tk

e^σ(^u−tkh )^Wtk+1−σ(^u−tkh )^Wtk+σ W_tk+ru+^σ2_2h(tku+hu−t²_k−htk−u²+utk)du = 1

T Xn−1 k=0

Z _tk+1

tk

e^σ(^u−tkh )^(W_tk+1^−W_tk^)−σ22 (u−tk)2

h +rue^{σ Wtk−σ22 tk}du We want to simplify this formula by a formal Taylor expansion.

f (x) = f (x₀) + f⁰(x₀)(x − x₀) + f⁰⁰(x₀)

2 (x − x₀)² +f⁰”(x₀)

3! (x − x₀)³...

VRAS 21

f (x) = X∞

k=0

f^k(x₀)

k! (x − x₀)^k e^x = 1 + x +x²

2 + x³ 3! + ...

e^{(r−σ22 )tk+σ Wtk} = S_tk

We will use the notation: ξ = u − tk ∈ (0, h) and ru = r(u − tk) + rtk, h = t_k+1− t_k, W_tk+1− W_tk = ∆ W_tk.

Z _tk+1

tk

e^σ(^u−tkh )^(Wtk+1−W_tk)−^σ2₂ ^(u−tk)_h ²+r(u−tk)du = Z _h

0

e(^{σ ξh})^{∆ W}_tk^{−(σ ξ)2}2h +rξdξ =

S_tk Z _h

0

µ

1 + σ ξ

h ∆ W_tk− (σ ξ)²

2h + rξ + O(h)

¶ dξ =

Stk

³ h + σ

2h∆ Wtkh²+ r 2h²

´

= hS_tk

µ 1 + σ

2∆ W_tk + rh 2

¶

Hence, we have the scheme:

Y_T^e,n = h T

Xn−1 k=0

S_tk µ

1 + rh

2 + σW_tk+1 − W_tk 2

¶

(2.11) We have to remark that the schemes(2.10 and 2.11) have the same rate of convergence ¹_n. We refer to [16] and [28] for the complete convergence results.

The third higher accuracy scheme is quite similar.

Z _T

0

Wudu has a normal density with respect to the Lebesgues measure on R, because the Brownian motion is the Gaussian process.

YT = 1 T

Z _T

0

Sudu = 1 T

Z _tk+1

tk

Sudu = 1 T

Xn−1 k=0

Stk

Z _tk+1

tk

S_u S_tkdu =

1 T

Xn−1 k=0

S_tk Z _tk+1

tk

e^{σ(Wu−Wtk)−σ22 (u−tk)+r(u−tk)}du

Using the Taylor expansion again, we get the scheme:

Y_T^p,n = 1 T

Xn−1 k=0

S_tk µ

h + rh² 2 + σ

Z _tk+1

tk

(W_u− W_tk)du

¶

(2.12) Remember that we have the price of the Asian call options with the fixed strike

V (T, S) = e^−rTE(1 T

Z _T

0

Sudu − K)+

The preceding consideration shows that for approximation integral Z _T

0

S_udu we can use the time schemes (2.10), (2.11) or (2.12).

Accordingly, the price, which we calculate by the Monte Carlo simulation, become:

V (T, S) = e^−rTE Ã

h T

Xn−1 k=0

Stk − K

!

+

(2.13) by the “Riemann” scheme,

V (T, S) = e^−rTE Ã

h T

Xn−1 k=0

S_tk µ

1 + rh

2 + σW_tk+1− W_tk 2

¶

− K

!

+

(2.14) by the “Trapezoidal” scheme,

V (T, S) = e^−rTE Ã

1 T

Xn−1 k=0

S_tk µ

h + rh² 2 + σ

Z _tk+1

tk

(W_u− W_tk)du

¶

− K

!

+

(2.15) by the third scheme.

Simulation of the Monte Carlo method. Now, we want to compute the prices at the maturity of the fixed strike Asian call options using “Rie- mann” scheme (2.13). We take into account the procedure of the Monte

VRAS 23 Carlo simulation, which was described at the beginning of this section, so we represent following

e^−rT M

XM j=1

(h T

NXt−1 k=0

S_tk − K)₊

where M denote the number of the Monte Carlo simulation (also called Monte Carlo loops), Nt is the number of the time steps, K is the strike price.

For each Monte Carlo simulation step M it computes the Geometrical Brow- nian Motion, Riemann Sum (RS), then makes the subtraction between the expected average price of the asset, computed by RS, and the strike price (RS − K) and if the average price of the asset will be more than the strike price (because we have call position) (RS > K),we multiply this expected payoff on risk-free rate e^−rT.

At the last point, the getting sum divides on the number of the simulating steps M to get value of the Asian options.

Thus, in the program we compute the following formula

mcV = P_M

j=1e^−rT(RS − K)₊

M (2.16)

In the Figure(2.4), the realization of the Monte Carlo method can be evinced. There are represented the dependence on the prices of the Asian options to the strike prices. The maturity date T = 1, the number of time step N_t= 1000, the interest rate r = 0.05, the volatility σ = 0.3, the initial asset price S₀ = 100, the strike price K = 70 : 130, the Monte Carlo loops M = 1000.

The prime result, demonstrated by this figure, is

• The price of the Asian call options depends on the strike price at which the holder of the option has the right to bay the underlying asset. The more the strike price (with respect to the initial price) concludes the less the option costs.

• How we can see the confidence interval is not small, it means that the prices of the Asian option fluctuate, and we do not have a possibility to know the exact value. It leads to a mistake. Inefficiency of the pricing Asian option is consisted of the using the ordinary Monte Carlo method. The accuracy of the Monte Carlo method may be improved, we will talk about it in the next section.

• The variance is more than 140. It has to be reduced, subsequently the inefficiency of Monte Carlo method will be disarmed.