Multi-Dimension Margrabe Formula

Abstract

Margrabe formula (1978), used in many financial engineering literature. Originally, it is to solve

E[Exp[ax], {ax+by<=k}]

where x and y is a binormal distribution, k is a constant. This article presents a multi-dimension version, which is used in some term structure model.

Sol

The C# source code for this calculation is in the SourceCodeRepository.

We want to calculate the integration over n-dimension space：

E[Exp[A.X], {B.X<=k}]

where A is of 1 by n matrix, B is of 1 by n matrix, X is of n by 1 matrix and entries of X is a multi-normal distribution with mean M, of n by 1 matrix, and covariance matrix S, of n by n matrix.

By elementary math, the moment generating function of X is

Exp[T.M + 1/2 T.S.Transpose[T]]

where T is of 1 by n matrix.

The key of this calculation is to change the variable X to another set of variable Y which are independent. Suppose Y=C.X where C is of n by n matrix. By elementary math, the moment generating function of Y is

Exp[T.N + 1/2 T.D.Transpose[T]]

where

N=C.M
D=C.S.Transpose[C]

Since all entries of Y are independent, it means D is of diagonal.

The C can be set up in this way that the first row of C is B. Once this C is found, denote

A.Inverse[C] by R
i-th row of Y by yi
i-th column of R by ri
i-th row of N by Ni
(i-th, i-th) entry of D by Di
standard normal accumulated distribution function by Normal

we have：

 　E[Exp[A.X], {B.X<=k}]
= E[Exp[R.Y], {y1<=k}]
= E[Exp[r1.y1], {y1<=k}]E[Exp[r2.y2]]E[Exp[r3.y3]]...E[Exp[rn.yn]]
= Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[N1.r1 + 1/2 r1^2.D1].Exp[N2.r2 + 1/2 r2^2.D2]...Exp[Nn.rn + 1/2 rn^2.Dn]
= Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[(N1.r1 + N2.r2 + .. + Nn.rn) + 1/2 (r1^2.D1 + r2^2.D2 + .. + rn^2.Dn)]
= Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[R.N + 1/2 R.D.Transpose[R]]
= Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[A.Inverse[C].C.M + 1/2 A.Inverse[C].C.S.Transpose[C].Transpose[A.Inverse[C]]]
= Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[A.M + 1/2 A.S.Transpose[A]]

We are done.

How to find C

First, start from S, we first change variable Z=V.X in this way z1=B.X, z2=x2, z3=x3,... , zn=xn. Denote covariance of Z by G, we have

G=V.S.Transpose[V]

Second, use the LowerUpperTriangleDecomposition algorithm on G, we have a lower triangle matrix L whose diagonal is 1 and a diagonal matrix D such that:

D=L.G.Transpose[L]

Therefore,

D=L.G.Transpose[L]
 =L.V.S.Transpose[V].Transpose[L]
 =LV.S.Transpose[LV]

LV is the C we want.

Example

We want to calculate

E[Exp[(x1-x2+x3+x4], {x1+x2+x3-x4<=k}]

where mean and covariance of x1, x2, x3, x4 is

1

2

3

4

and

1	0.1	0.1	0.1
0.1	1	0.1	0.1
0.1	0.1	1	0.1
0.1	0.1	0.1	1

By the above notation, we have:

M=

1

2

3

4

S=

1	0.1	0.1	0.1
0.1	1	0.1	0.1
0.1	0.1	1	0.1
0.1	0.1	0.1	1

A=

1

-1

1

1

B=

1

1

1

-1

Therefore,

V=

1	1	1	-1
0	1	0	0
0	0	1	0
0	0	0	1

G=

4	1.1	1.1	-0.7
1.1	1	0.1	0.1
1.1	0.1	1	0.1
-0.7	0.1	0.1	1

L=

1	0	0	0
-0.275	1	0	0
-0.35483871	0.290322581	1	0
0.5	-0.590909091	-0.590909091	1

C=LV=

1	1	1	-1
-0.275	0.725	-0.275	0.275
-0.35483871	-0.064516129	0.64516129	0.35483871
0.5	-0.090909091	-0.090909091	0.5

D=

4	0	0	0
0	0.6975	0	0
0	0	0.638709677	0
0	0	0	0.531818182

R=

0.1

-1.161290323

1.181818182

2

N=

2

1.45

2.870967742

2.045454545

Therefore, N1=2, D1=4, R1=0.1, A.M=6,

A.S.Transpose[A]=4

And the answer is

 Normal[(k-(N1+D1.r1))/Sqrt[D1]].Exp[A.M + 1/2 A.S.Transpose[A]]
=Normal[k/2-1.2].Exp[8]