Application Mellin transforms to the Black Scholes equations Текст научной статьи по специальности «Математика»

© T.A. Vasileva, O.E. Vasileva, 2009

ПРИКЛАДНАЯ МАТЕМАТИКА

УДК 517.956.44 : 336.01

APPLICATION MELLIN TRANSFORMS TO THE BLACK - SCHOLES EQUATIONS

T.A. Vasileva, O.E. Vasileva

In this paper we will consider the Mellin transform derived on Black - Scholes formulae in year 2003 by Panini, Srivastav and Jodar [6; 7]. The authors use Mellin transforms to derive at first an equation for the price of a European put on a single underlying stock. This case will be extended to the Amercan put option. It is assumed, there are no dividends. Nowadays derivative markets become extremely popular, this popularity even exceeds that of the stock exchange [4; 5; 10]. Option price estimation as the most interesting of the derivatives has many approaches and is multifaceted [2;3;8]. Authors would like to thank Professor E.I. Vasilev for the target setting and useful discussions.

1. European Put Option

One of the techniques to solving Black - Scholes equations called Mellin transformation, and suggested by Panini [6; 7] for basket options will be considered and derived new

integral representations for the price and the free boundary of the American option in one dimension. We will also improve the numerical part of Panini’s computations [6; 7; 11] by using the Newton’s method [9] for the free boundary condition.

Black - Scholes equation [1] for the price of a European put p(S,t) reads

^ + 1 a2S2 + rSdP rp = 0 (1)

dt + 2aS OS 2 + dS - rp = 0 (1)

where t is current time, 0 < t < T, S(t) is asset price, 0 < S < to, T is expiry date, r is the interest rate, a is volatility of the market prices.

Note that p(S, t) ^ 0 as S ^ to and satisfies the terminal conditions

p(S,T) = 9(S) = (K - S)+, here K is exercise price and the boundary condition at S = 0

p(0, t) = K ■ e-r(T-t).

Let p(v, t) denote the Mellin transform of p(S, t) which is defined by the relation

p(v,t) = / p(S,t) ■ Sv dS.

10

where is v is a complex variable.

So, we get the inverse Mellin transform

1 /»c+io

p(S, t) = -—: p(v,t) ■ S-vdv. (2)

2ni Jc-iCO

Hence, the price of the European put is given by

— rc+io

p(S, t) = — <9(v) ■ e2"2h(v)(T-t) ■ S-vdv, (3)

2ni J c-ico

where 0(v) is the Mellin transformation of the above mentioned terminal conditions. By

transforming variables using several algebraic transformations, the above expression can be simplified to

p(S, t) = K ■ e-r(T-t)N(—d2) - S ■ N(-d—). (4)

Here N(■) is the N(0, —) the normal distribution function and

d

log(S/K) + (r + 2a2)(T - t)

l

a\/T - t ’

A _log(S/K) + (r - 2a2)(T - t)

d2 —----------------. ------------.

2 a\fT-t

2. American Put Option

If we consider now the American put option, the early exercise feature of this option gives rise to a free boundary problem. As far as we know, there exists no closed-form analytical expression for the value of American put or its associated free boundary.

The Black-Scholes equation for the price of an American put p(S, t) satisfies the nonhomo-geneous equation

I + 2ff2S'2 dS2+rSfS- rp =f (S't)- (5)

where 0 < t < T and 0 < S < to The inhomogenity in (4) is given by

f f(St) = f -r ■ K, if 0 <S < S*(t) f = f (S,t) = \ 0 if S > S*(t)

The final time condition is

P (S, T) = 0(S) and the free boundary is determined by

P(SV) = K - s*, dPS \s=st = -—.

As before, the authors take the Mellin transform and for the price P(S, t) they get

P(S, t) = p(S, t) + [ [+ s-* (S*(x))* e1 -2h(v)(x-t) dv dx. (6)

2ni Jt Jc-io v

Substituting S = S*(t) from (—2) and using (——), they get the following integral equation for the free boundary

rK rT rc+io 1 S*(t) 1 2

K - S*(t) = p(S’(«),*) + + 7^ -(^fr)-*e1 'h(*)(x-t) dvdx. (7)

2nU t ic-io v S *(x)

3. Motivation of Newtons method for an American Put Option

Now we want provide a motivation of using the Newtons method [—6] for the calculation of free boundary S*(t).

The equation for the price P(S, t) of the —-basket American put option is given by

P№ t) = p(S, t) + rK [ e-N (--i= (l„ (s;S^) ) + £(r - 1 a2)) d£. (8)

Boundary condition on a free boundary gives integral equation for P(S, t)

S*(t) = K - p(S*(t), r)-

-rK/T eM-OTf Hslr-£)) +£(r - - a2)))df (9)

The “European part” looks as follows

p(S) = Ke-rTN(-d2) - SN(-di), (10)

where

S , a2

K

Panini has used Method of simple Iteration with initial S*(t) = K

d2 = di - di = -^= ^logKS + (r + y)7"^

(11)

Sn (t ) = K - p(sn-i(T ),t )-

- rKJ,e-rSN(-a^f (ta(s^)+£(r - 2a2))) (12)

It is not optimal to use it. So, we can consider the possibility to improve, i.e. accelerate this computation using faster method, for example, Newtons Method. Rewrite the equation in form F(x) = 0, (x(t) = S*(t))

x“+i = x“ - x" - Ka+ p(x“'Tk). (13)

— + fi (x",Tk) ' '

The Newton’s method has the quadratic rate of convergence. But it more laborious as demands calculation of matrix Jacobian. Also, the Newton’s method of converges from

close initial approximation x0. For more detailed finding-out of character of convergence we shall calculate the first and the second derivative of the function F(x) = x - K + +p(x,Tk). We differentiate F(x) as a complex function.

Fx(x) = — + p(H—) = — + (Ke rTN (-P + aVT - xN (-P^ (14)

dx dx

where p = ^ [ln(K) + (p + °2)t]-

Fx(x) = — - Ke rTN"/(-p + aV^)"dp) + xN (-p)dxx- N(-p)- (15)

We shall consider the obvious property for the function N(x)

—

1 — 2

N (a) = e, (16)

from which follows

N' (a + h) = eee-ahe"^ = N (a)e-ahe"^. (17)

л/2п л/2п

In our case it turns out

2

N' (-p + aV=) = N' (-p)ep^ e-=

N' (—p)e ^ +r+it T e-^ = N' (-p) -XepT.

K

(18)

After substitution N-(-p) to the Fx(x) we have

Fx(x) = 1 - Ke-rTN'(-p)-|epTdP) + xN'(-p)dP) - N(-p) = 1 - N(-p). (19)

As function 0 < N(x) < — we have received, that the first derivative is always positive.

Fx(x) = — - N (-P) > 0- (20)

Let find the second derivative

Fxx(x) = (Fx(x)) = 0 - N(-P)dd=^ = N(-P) ^ ^^ > 0, (21)

dx dx dx xa\jT

since

dp 1

dx xa^/T'

We have received that the second derivative is always positive

Fxx(x) = > 0.

xa\/T

Then the Newton’s method provides the monotonous convergence on the right. Hence it is quite reliable, provided that initial approach more than a required root. Therefore the Newton’s method quite approaches for the solution of the given system. The method will be stable as well, if to function a constant will be added. It does not influence first and second derivative of functions F(x).

2

4. The interpretation of results

Now the results for —-basket option will be discussed with the help of graphs and tables.

The results of our programm and its correctness by using Newtons method can be checked in compare with results which were obtained by Panini. The left figure shows the results of Panini work used the recursive function for American put option. The right figure shows the results by applying the Newtons method for the free boundary S*(t). The similarity of both graphs is obvious.

Figure 1.

Now we want to present the difference of computations for the free boundary by American put option on a single underlying asset, i.e. an integral error from function M.

The size of the grid, M Trapezium Trapezium Plus

50 0.003728 0.002714

100 0.001868 0.000610

200 0.000921 0.000099

400 0.000452 0.000037

800 0.000227 0.000018

Table 1. The average error by the computation for the free boundary (integral error from function M)

The table — shows the convergence velocity of the approximated solution for the free boundary S*(t) and (t e [0,T]) by the increasing of the grid size M on the interval [0,T]. The computations took place with the following initial values K = 45; T = 0, 5833; a = 0, 5; r = 0, 0488. In order to check the given program for convergence by using Newton’s method, we need to compare the difference between graphs with the differ-

ent number of nodes M along time t. In the left column the gride size is given with M = 50, 100, 200, 4000, 8000 nodes in form of doubling of the previous value. The norm is calculated like the average of 25 nodes on the interval, which is common for all given grids.

The norm of difference is calculated in the same way. In the central column the trapezoidal rule for the integral computation is used for the convergence. The trapezoidal rule has the second order of accuracy, but the error for the free boundary is reduced just linearly in this case. It means if the grid size is double increased, the error is in approximately 2 times reduced. It depends on the behavior of the integrand function in the integral

function provides the considerable error in the neighborhood e = 0, which influences all of the values s*(t).

For the elimination of this error the trapezoidal rule should be modified on the end of the interval of the grid. The idea of the modification of the trapezoidal rule is the substitution on the end of the interval of the grid instead

the whole interval [0, t]. The results of this modification called Trapezium Plus are given in the right column of the table. By the M = 50 ^ 200 the convergence velocity is close to the quadratic, further the velocity decreases up to linear. In a whole the convergence of the modification of the method is really faster in compare with the usual trapezoidal rule. As was already mentioned the convergence velocity of the trapezoidal method is close to linear according to the first table. It allows us to use the extrapolation of the solution with high accuracy under the formulae

The solution has been received for M = 3200 nodes. By means of this solution it is possible to find errors (i.e. the deviation from the exact solution) of the approximate solutions and their dependence from t. In the table 2 the average relative errors for the various grids are presented. In distinction with the table 1 the comparisons for all M are taken with S(t) here. The relative errors is computed by the following formulae

which by very small e has the behavior f (e) ~ A + B\fe. The trapezoidal rule for such

2

another formulae is used

2f (°) + 25f (9/25Дт ) + 9f (Ar)

(22)

36

This formulae provides the exact value of the integral on the interval [0, At] for the function like f (e) ~ A + B^ and increases the accuracy of the integral computation on

S(t) = 2Sm (t) — Sм (t).

SSm = ||(Sm (t ) — Sm (t ))/Sm (t )|| = | Sm (t) — Sm (t) /Sm (t ).

i=1

The size of the grid, M Trapezium Trapezium Plus

50 0.003722 0.000302

100 0.001942 0.000180

200 0.001000 0.000100

400 0.000511 0.000054

800 0.000260 0.000029

Table 2. Relative average error of computation for free boundary (integral error from function M)

From the table we can see that by applying of modification of the trapezoidal rule the average error is less approximately in 10 times in compare with the usual trapezoidal rule.

5 sM 1000

\

'

\ M=50

^^^M=100

, ^ M=200

‘M=400^_ —_

s sM 1000

0.4

0.3

0.2

0.1

\

\ M=50

v~v\ M=100

V/V'A.Ar~~ |M=400 . _ M=200 wvWVwn/W ''vWVwwW

0.1

0.2

0.3

0.4

0.5

0.1

0.2

0.3

0.4

0.5

Figure 2. Isolines of the space for the function p(s,r) with different choices for sigma option

with K = 45; a = 0,3; T = 0, 5833

So, on the graph 2 we have the dependance 5Sm(t) on t G [0,1; T] of the relative errors for the different choices of the grid size M. On the left graph the usual trapezoidal method is used, on the right one the trapezoidal rule is presented. On the both graphs the errors have been increased, i.e. zoomed in 1000 times. We shall especially note, that the vertical scale on the right figure is increased in 10 times. Because of the increased scale on the right graph are more distinctly visible the errors oscillation.

By comparing of both graphs it is obvious, that using the modification error decreases in 10 times on the interval t G [0,1;T]. The uniform reduction of an error took place here though the correction (22) was applied only on edge of a integration interval. Thus, the suggested modification of trapezoidal method allows to raise the grid accuracy of calculations in 10 times without any change. It is necessary to mention as well that for both methods reduction of an error takes place with the growth of t.

The efficiency of the modification can be estimated in computational burden. It is visible from the table and graphs, for example, that the error by the modification on grid M = 100 is less, than without it on a grid with M = 800 nodes. It means that the simple modification by the same accuracy allows to use more coarse grid with less nodes and in 8 x 8 = 64 time to reduce an operating time of the program.

The Figure 3 shows the isolines of the function of the American put option. As we

can see the smooth conjunction of the family of the strict vertical isolines with the curve linear isolines. The boundary where the vertical isolines begin to curve is the free boundary

S *(t ).

Figure 3. Isolines of the space for the function p(s,r) with different choices for sigma option

with K = 45; a = 0,3; T = 0, 5833

Bibliography

1. Black, F. The Pricing of Options and Corporate Liabilties / F. Black, M. Scholes

// Journal of Political Economy. — 1973. — 81. — P. 637-654.

2. Clewlow, L. Implementing Derivatives Models / L. Clewlow, C. Strickland. — Chichester : John Wiley and Sons, 1998. — 330 p.

3. Cox, J. C. Option pricing: A simplified approach / J. C. Cox, S. A. Ross, M. Rubinstein

// Journal of Financial Economics. — 1979. — V. 7. — P. 229-263.

4. Hull, J. Introduction to futures and options markets / J. Hull. — New Jersey : Prentice-

Hall, Inc., 1998. — 400 p.

5. Kwok, Y.K. Mathematical Models of Financial Derivatives / Y. K. Kwok. — Singapure : Springer Verlag, 1998. — 386 p.

6. Panini, R. Option Pricing with Mellin Transforms / R. Panini, R. P. Srivastav // Applied Mathematics Letters. — 2003. — V. 18 (4). — P. 471-474.

7. Panini, R. Option Pricing with Mellin Transforms. PhD Thesis / R. Panini. — N.-Y. : Stony Brook University, 2004.

8. Ross, M. S. An elementary Introduction to mathematical Finance / M. S. Ross. —

Cambridge : Cambridge University Press, 2003. — 270 p.

9. Samarsky, A. A. Numerical methods (in Russian) / A.A. Samarsky, A.V.Gulin. — M. : Mir, 1989. — 432 p.

10. Wilmott, P. Derivatives / P. Wilmott. — Chichester : Wiley, 1989.

11. Wurfel, A. Analytical und numerical solution of the Black - Scholes equation for

European and American Basket Options using Mellin transformations. Diploma Thesis (in German) / A. Wurfel. — Berlin : TU, 2007.

Аннотация

ПРИМЕНЕНИЕ ПРЕОБРАЗОВАНИЙ МЕЛЛИНА К РЕШЕНИЮ УРАВНЕНИЙ БЛЭКА - ШОУЛЗА

Т.А. Васильева, О.Е. Васильева

Рынок вторичных ценных бумаг в настоящее время весьма популярен в России и за рубежом. В данной статье рассматривается применение интегральных преобразований Меллина к решению уравнения Блэка-Шоулза. В начале работы авторы применили преобразование Меллина к решению уравнения для определения цены Европейского опциона Пут. Этот подход далее распространяется на вычисление цены Американского Пут опциона.