Metoda Monte Carlo

Rafał Kucharski

[e-mail]

www.math.us.edu.pl/rk

Szkoła Letnia Matematyki Finansowej, Tarnów 2012

John was assisted by Henry, a foreign quant whose English was

incomprehensable, but who was believed to be at least equally com-
petent in risk-management methods. John knew no math, he relied
on Henry. (. . . ) Henry was a graduate student in Operations Rese-
arch when John hired him. His specialty was a field called Compu-
tational Finance, which, as its name indicates, seems to focus solely
on running computer programs overnight. Henry’s income went from
$50,000 to $600,000 in three years.

Nassim Nicholas Taleb, Fooled by Randomness

Zaczynamy od prostej funkcji. . .

f (x) = 4 x (1 − x)

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Iterujemy. . .

x

0

= 0.1,

x

n+1

= f (x

n

),

n = 1, . . . , 1000

0

200

400

600

800

1000

0.0

0.2

0.4

0.6

0.8

1.0

Ciąg deterministyczny zachowuje się jak losowy!

Histogram (arcus sinus)

0.0

0.2

0.4

0.6

0.8

1.0

0

50

100

150

200

Generowanie liczb losowych

• liczby prawdziwie losowe (true random numbers)

• liczby pseudo-losowe (pseudo-random numbers)

• liczby quasi-losowe (quasi-random numbers)

Liczby prawdziwie losowe:

• rzut monetą, rzut kostką, koło ruletki,

• zachowanie człowieka lub jego efekty: ruchy myszką komputerową,

kliknięcia na klawiaturze, ruch głowic twardego dysku,

• zjawiska fizyczne których nieprzewidywalność sięga praw mechaniki

kwantowej: rozkład radioaktywny, szum termiczny, szum radiowy

• zdjęcia chmur, ”Lava lamps”

Liczby quasi-losowe:

• quasi-random numbers, low discrepancy numbers)

• regularne – nielosowość widoczna ”gołym okiem”,

• bardziej równomierne niż ciągi pseudolosowe,

• dają mniejszy błąd – lepsze tempo zbieżności,

• najbardziej popularne: ciągi Haltona, ciągi Sobola,

• metoda Quasi-Monte Carlo.

Halton

Sobol

Liczby pseduo-losowe

Własności generatorów liczb pseudolosowych

(RNG random number generators)

• okresowość (pełny okres daje rozkład jednostajny),

• odtwarzalność, powtarzalność (reproducibility),

• szybkość,

• przenośność (portability),

• losowość.

Generator liniowy kongruencyjny:

• Ustalamy liczbę całkowitą M (duża, najlepiej pierwsza, np. 2

32

− 1).

• Wybieramy a, b ∈ {0, 1, . . . , M − 1}.

• Wybieramy ziarno x

0

∈ {0, 1, . . . , M − 1}.

• x

i

= (ax

i−1

+ b) mod M , i = 1, 2, . . .

• U

i

= x

i

/M ∈ [0, 1].

Inne popularne generatory:

• Fibonacciego: x

n

:= x

n−1

+ x

n−2

(mod m),

• multiple recursive generators: x

n

:= a

1

x

n−1

+ · · · + a

k

x

n−k

(mod m),

• bitowe (w tym feedback shift registers),

• inwersyjne: x

n

:= ax

−1

n−1

+ b.

Zobacz w R: ?RNG

Fakt. Niech F będzie dystrybuantą ciągłą i ściśle rosnącą.
Jeśli U ∼ U [0, 1], to F

−1

(U ) jest zmienną losową o dystrybuancie F .

P (F

−1

(U ) ¬ x) = P (F (F

−1

(U )) ¬ F (x)) = P (U ¬ F (x)) = F (x), x ∈ R.

Definicja. Dla dowolnej dystrybuanty F : R → [0, 1] uogólniona funkcja odwrotna
F

−

: [0, 1] → [−∞, ∞] dana jest wzorem

F

−

(u) = inf{x : F (x) ­ u}.

Fakt. Niech F będzie dowolną dystrybuantą.
Jeśli U ∼ U [0, 1], to F

−

(U ) jest zmienną losową o dystrybuancie F .

Wniosek: Wystarcza nam generowanie liczb o rozkładzie jednostajnym.

Zadanie. Wygeneruj liczby o poniższych rozkładach za pomocą transformacji
odwrotnej z rozkładu jednostajnego.

• rozkład Cauchy’ego: F (x) =

1
2

+

1

π

arc tg(x),

• rozkład wykładniczy: F (x) = 1 − exp(−x), x ­ 0,

• rozkład arcusa sinusa: F (x) =

2

π

arc sin(

√

x), 0 ¬ x ¬ 1,

(czas w jakim proces Wienera osiąga maksimum na przedziale [0, 1]),

• rozkład Rayleigh’a: F (x) = 1 − e

−2x(x−b)

, x ­ b,

(maksimum procesu Wienera na [0, 1] pod warunkiem W

1

= b),

Rozkłady dyskretne:

• Chcemy wygenerować zmienną X o rozkładzie P skupionym na zbiorze N,

• Wyznaczamy liczby p

0

= P(X ¬ 0), p

1

= P(X ¬ 1), p

2

= P(X ¬ 2), . . .

• Gdy rozkład jest skupiony na zbiorze nieskończonym, poprzestajemy na p

k

dostatecznie bliskim 1, np. p

k

= 0.99999.

• X = k, jeśli p

k−1

< U < p

k

, gdzie U ∼ U [0, 1].

Czasem przeszukiwanie od k = 0 może być nieefektywne.
Przykładowo: dla P ois(100) większość obserwacji leży w (70, 130),
zatem warto rozpocząć sprawdzanie od p

70

.

Rozkład normalny:

• Proste przybliżenie: X =



P

12

i=1

U

i



− 6, gdzie U

i

∼ U [0, 1] są niezależnymi

zmiennymi losowymi (iid),

• Box-Muller: U

1

, U

2

∼ U [0, 1], iid,

X

1

=

r

−2 log(U

1

) cos(2πU

2

),

X

2

=

r

−2 log(U

1

) sin(2πU

2

).

• Marsaglia-Bray:

while (X > 1)

losuj U

1

, U

2

∼ U [0, 1]

U

1

← 2U

1

− 1, U

2

← 2U

2

− 1

X ← U

2

1

+ U

2

2

Y ←

r

−2 log(X)/X

Z

1

← U

1

Y , Z

2

← U

2

Y

return Z

1

, Z

2

.

• Wielowymiarowy rozkład normalny N

d

(µ, Σ):

– Niech Σ = AA

T

(np. rozkład Cholesky’ego), X ∼ N

d

(0, I

d

),

– Wówczas AX ∼ N

d

(0, Σ), AX + µ ∼ N

d

(µ, Σ).

Dodatek: algorytm rozkładu Cholesky’ego dla macierzy dodatnio określonej:

A ← 0 (macierz d × d)

for j = 1, . . . , d

for i = j, . . . , d

v

i

← Σ

ij

for k = 1, . . . , j − 1

v

i

← v

i

− A

jk

A

ik

A

ij

← v

i

/

√

v

j

return A.

W przypadku d = 2 macierz korelacji ma postać






1 ρ

ρ 1






,

a jej macierz Cholesky’ego przybiera formę:






1

0

ρ

√

1 − ρ

2






.

Możemy zatem zapisać powyższy algorytm następująco:

• Wylosuj X

1

, X

2

∼ N (0, 1), niezależne,

• Oblicz Y

1

= X

1

, Y

2

= ρX

1

+

√

1 − ρ

2

X

2

.

Istotnie, wówczas E(Y

1

) = E(Y

2

) = 0, Var(Y

1

) = Var(X

1

) = 1,

Var(Y

2

) = ρ

2

Var(X

1

) + (1 − ρ

2

) Var(X

2

) = 1,

Corr(Y

1

, Y

2

) =

cov(Y

1

, Y

2

)

r

Var(Y

1

) Var(Y

2

)

= cov(X

1

, ρX

1

+

r

1 − ρ

2

X

2

) =

= cov(X

1

, ρX

1

) + cov(X

1

,

r

1 − ρ

2

X

2

) = ρ

Zadanie. Wygeneruj wektory liczb X, Y o rozkładzie normalnym standardowym,
dla których współczynnik korelacji ρ(X, Y ) = 0.6.

Metoda akceptuj-odrzuć (Accept-Reject)

• Chcemy generować liczby z rozkładu o gęstości f

(znanej z dokładnością do stałej multyplikatywnej),

• Potrzebujemy prostszej do symulacji gęstości g takiej, że

– f i g mają takie same nośniki (g(x) > 0 ⇐⇒ f (x) > 0),

– istnieje M , że f (x)/g(x) ¬ M dla wszystkich x,

• Algorytm:

1. generuj Y ∼ g, U ∼ U [0, 1],

2. jeśli U ¬

1

M

·

f (Y )

g(Y )

, to X = Y ; w przeciwnym wypadku wróć do 1.

• Powyższy algorytm działa również w przypadku wielowymiarowym.

Zadanie. Używając metody akceptuj-odrzuć wygeneruj liczby z rozkładu
o gęstości:

f (x) ∝ exp(−x

2

/2){sin(6x)

2

+ 3 cos(x)

2

sin(4x)

2

+ 1}.

• Narysuj f (x) i znajdź takie M , że

f (x) ¬ M g(x),

gdzie g(x) = exp(−x

2

/2)/

√

2π jest gęstością rozkładu normalnego.

Możesz użyć funkcji optimise, aby znaleźć najmniejszą wartość M .

• Oszacuj stałą normalizującą.

• Porównaj histogram ze znormalizowaną gęstością.

Rozkłady warunkowe:

• X ma dystrybuantę F ,

• Chcemy wygenerować X pod warunkiem a ¬ X ¬ b,

• V = F (a) + (F (b) − F (a)) · U , gdzie U ∼ U [0, 1],

• F

−

(V ) ma pożądany rozkład,

Zadanie (Probabilistyka, A. i E. Plucińscy, str. 228, Zadanie 15).

Twierdzenie (Probabilistyka, A. i E. Plucińscy, str. 223, Twierdzenie Lapunowa
4.9). Niech {X

n

} będzie ciągiem niezależnych zmiennych losowych i niech dla

n = 1, 2, . . . będzie

m

n

= E(X

n

),

σ

2

n

= E((X

n

− m

n

)

2

) > 0,

b

n

= E(|X

n

− m

n

|

3

),

C

n

=

r

σ

2

1

+ · · · + σ

2

n

,

B

n

=

3

√

b

1

+ · · · + b

n

,

U

n

=

1

C

n

n

X

k=1

(X

k

− m

k

),

F

n

(u) = P (U

n

< u).

Jeśli lim

n→∞

B

n

C

n

= 0, to dla każdego u

lim

n→∞

F

n

(u) =

1

√

2π

Z

u

−∞

e

−x

2

/2

dx,

czyli ciąg losowy {U

n

} jest zbieżny według dystrybuant do zmiennej losowej o

rozkładzie normalnym N (0, 1).

Proces Wienera

• W

0

= 0 p.n., trajektorie t 7→ W

t

są ciągłe p.n.

• przyrosty W

t

1

, W

t

2

− W

t

1

, . . . , W

t

n

− W

t

n−1

są niezależne,

• W

t

− W

s

∼ N (0, t − s), 0 ¬ s ¬ t,

Symulacja trajektorii procesu Wienera

• Chcemy symulować wartości (W (t

1

), . . . , W (t

n

)) w ustalonym zbiorze

punktów 0 < t

1

< · · · < t

n

.

• Niech Z

1

, . . . , Z

n

będą niezależnymi zmiennymi losowymi N (0, 1),

• Dla standardowego ruchu Browna W (0) = 0 oraz

W (t

i+1

) − W (t

i

) =

√

t

i+1

− t

i

· Z

i+1

,

i = 0, . . . , n − 1.

0.0

0.2

0.4

0.6

0.8

1.0

−2.0

−1.5

−1.0

−0.5

0.0

0.5

time

BM

Wielowymiarowy proces Wienera

• Proces W (t) = (W

1

(t), . . . , W

d

(t)) jest standardowym d-wymiarowym

ruchem Browna , gdy jego współrzędne W

i

(t) są standardowymi

jednowymiarowymi ruchami Browna oraz W

i

i W

j

są niezależne dla i 6= j.

• Niech µ ∈ R

d

oraz Σ ∈ R

d×d

będzie dodatnio określoną macierzą.

Proces X nazywamy ruchem Browna z dryfem µ i kowariancją Σ,
jeśli X ma ciągłe trajektorie o niezależnych przyrostach

X(t + s) − X(s) ∼ N (tµ, tΣ).

Proces X nazywamy ruchem Browna o stałej kowariancji (Σ jest stała).

• Dla symulacji Monte Carlo ze skorelowanym ruchem Browna potrzebujemy

generować próbki o rozkładzie N (δtµ, δtΣ), a jeśli δt jest ustalone problem
redukuje się do losowania z rozkładu N (µ, Σ) (co już potrafimy).

Zaproponowana przez Stanisława Ulama w czasie projektu Manhattan, jako
sposób obliczania skomplikowanych całek pojawiających w modelowaniu reakcji
łańcuchowych. Rozwinięta przez J. von Neumanna, N. Metropolisa i innych.

Źródło: www.atomicarchive.com

• Jeśli X ∼ U [0, 1] oraz h : [0, 1] → R, to wartość oczekiwana zmiennej h(X)

jest równa całce:

E[h(X)] = I(h) =

Z

1

0

h(x) dx.

• Podobnie w przypadku wielowymiarowym: X ∼ U [0, 1]

d

oraz h : [0, 1]

d

→ R

E[h(X)] = I(h) =

Z

[0,1]

d

h(x) dx.

• Z drugiej strony, przyjmując

I

N

(h) =

1

N

N

X

n=1

h(x

n

),

gdzie x

n

są niezależnymi realizacjami X, mamy

E[I

N

(h)] = E[h(X)] = I(h)

(nieobciążoność) oraz z MPWL

lim

N →∞

I

N

(h) = E[h(X)] = I(h).

• Ponadto

Var(I

N

(h)) =

1

N

2

N

X

k=1

Var h(x

k

) =

1

N

Var h(X).

• Błąd metody Monte Carlo jest proporcjonalny do N

−1/2

.

• Co ważniejsze, błąd Monte Carlo jest proporcjonalny do N

−1/2

, niezależnie

od wymiaru!

• Błąd = O(N

−1/2

), Czas = O(N )

⇒

jeśli chcemy zwiększyć dokładność

10-krotnie, musimy zwiekszyć ilość próbek i czas obliczeń 100-krotnie.

Model Blacka-Scholesa (1973)

Ceny akcji modeluje geometryczny proces Wienera:

S

t

= S

0

exp



(µ −

1
2

σ

2

)t + σW

t



,

lub równoważnie:

ln(S

t

) = ln(S

0

) + (µ −

1
2

σ

2

)t + σW

t

,

gdzie

• µ jest stopą aprecjacji (wyznacza trend ceny),

• σ jest współczynnikiem zmienności (volatility),

• parametry µ, σ > 0 są stałe (i deterministyczne),

• akcja nie wypłaca dywidendy,

• stała stopa procentowa wolna od ryzyka r,

W modelu B-S zmienna losowa S

t

ma rozkład lognormalny, tzn. zmienna ln(S

t

) ma

rozkład normalny.

Wzory Blacka-Scholesa

• Cena europejskiej opcji kupna:

C = S

t

Φ(d

1

(S

t

, τ )) − Ke

−rτ )

Φ(d

2

(S

t

, τ )),

• Cena europejskiej opcji sprzedaży:

P = Ke

−rτ

Φ(−d

2

(S

t

, τ )) − S

t

Φ(−d

1

(S

t

, τ )),

• gdzie τ = T − t, funkcja Φ jest dystrybuantą standardowego rozkładu

normalnego oraz

d

1

(s, τ ) =

ln(s/K) + (r + σ

2

/2)τ

σ

√

τ

,

d

2

(s, τ ) = d

1

(s, τ ) − σ

√

τ =

ln(s/K) + (r − σ

2

/2)τ

σ

√

τ

.

• Uwaga: cena opcji nie zależy od µ.

Wycena opcji niezależnych od trajektorii

• Instrumenty finansowe wyceniamy korzystajac z wzoru

Π

0

(X) = B

0

E

Q

[X/B

T

],

gdzie Q jest miarą neutralną wobec ryzyka,
B

t

– ceną w chwili t instrumentu wolnego od ryzyka (numeraire) oraz

X jest wypłatą z danego instrumentu w chwili T .

• Cena akcji w modelu Blacka-Scholesa względem miary neutralnej wobec

ryzyka, dane są wzorem:

S

t

= S

0

exp





(r −

1

2

σ

2

)t + σW

t





.

• Zmienna losowa W

T

ma rozkład normalny o średniej 0 i wariancji T .

• Możemy przyjąć W

T

=

√

T Y , gdzie Y ∼ N (0, 1).

• Wypłata z europejskiej opcji call wynosi

X = f (S

T

) = (S

T

− K)

+

,

a dla europejskiej opcji put

X = g(S

T

) = (K − S

T

)

+

.

gdzie K jest ceną wykonania.

• Cena instrumetu wolnego od ryzyka (obligacja) to

B

t

= exp(rt),

t ∈ [0, T ].

• Cenę opcji wyznaczamy jako wartość oczekiwaną

V = E

Q

[f (S

T

)/B

T

] = exp(−rT )

Z

1

0

f (S

T

)dU,

gdzie

S

T

= S

0

exp((r −

1

2

σ

2

)T + σ

√

T Y )

oraz Y = Φ

−1

(U ) ∼ N (0, 1), U ∼ U [0, 1].

• W praktyce:

– Losujemy N niezależnych realizacji zmiennych o rozkładzie

normalnym Y

(1)

, . . . , Y

(N )

,

– Obliczamy ceny akcji w chwili wykonania S

(1)

T

, . . . , S

(N )

T

dla

kolejno wylosowanych próbek,

– Obliczamy zdyskontowane wypłaty z opcji:

exp(−rT )f (S

(1)

T

), . . . , exp(−rT )f (S

(N )

T

),

– Wyznaczamy cenę opcji jako zdyskontowaną wartość oczekiwaną:

Π

0

(X) =

1

N

N

X

i=1

exp(−rT )f (S

(i)

T

).

Zadanie. Oblicz ceny opcji dla r = 0.05, σ = 0.2, T = 1, S

0

= 110, K = 100.

Porównaj z wartościami dokładnymi (wzór Blacka-Scholesa). Wyznacz przedział
ufności, w którym prawdziwa wartość leży z prawdopodobieństwem 99.7%.
Wyznacz ilość próbek N , tak by powyższy przedział miał szerokość 0.01.

Wycena opcji zależnych od trajektorii

• Dla t

2

> t

1

mamy

S

t

2

= S

t

1

exp





(r −

1

2

σ

2

)(t

2

− t

1

) + σ(W

t

2

− W

t

1

)





• Losujemy N niezależnych trajektorii ceny S

(1)

t,T

, . . . , S

(N )

t,T

na przedziale [t, T ],

• Dla log-normalnego procesu ceny akcji i równoodległych chwil t

i+1

− t

i

= δt,

i = 0, . . . , n − 1, mamy

S

t

i+1

= S

t

i

exp((r −

1

2

σ

2

)δt + σ

√

δtZ

i

),

gdzie Z

i

∼ N (0, 1).

• Obliczamy wypłatę V

(n)

na każdej trajektorii n = 1, . . . , N ,

• Ceną opcji jest Π

0

= exp(−rT )

1

N

P

N

n=1

V

(n)

.

Przykład: opcja azjatycka (asian option)

• Europejska opcja call na średnią cenę M akcji o wypłacie

X =





1

T

T

X

t=1

S

t

− K





+

,

gdzie S

t

jest ceną akcji w dniu t.

• Wielowymiarowy proces Wienera przydaje się przy wycenie opcji opartych o

kilka instrumentów bazowych.

Przykład: opcja koszykowa (basket option)

• Europejska opcja call na średnią cenę M akcji o wypłacie

X =





1

M

M

X

i=1

S

T,i

− K





+

,

gdzie S

T,i

jest ceną i-tej akcji w chwili T .

Redukcja wariancji

• antithetic variables (metoda odbić lustrzanych)

– Losujemy 2 próbki, X

1

i X

2

, z tego samego rozkładu.

– Wariancja wynosi

ε = Var






h(X

1

) + h(X

2

)

2






=

1

2

[Var(h(X

1

) + cov(h(X

1

), h(X

2

))],

– Jeśli h(X

1

) i h(X

2

) są niezależne to ε =

1
2

Var(h(X

1

),

– Jeśli h(X

1

) i h(X

2

) będą ujemnie skorelowane, to ε <

1
2

Var(h(X

1

),

– Jeśli X

1

, X

2

∼ U [0, 1] ujemną korelację uzyskujemy przyjmując

X

2

= φ(X

1

), gdzie φ : [0, 1] → [0, 1] jest funkcją malejącą oraz

φ(X

1

) ∼ U [0, 1]; np. φ(x) = 1 − x.

– Każdą wylosowaną próbkę możemy wykorzystać 2 razy!

– Uwaga: jeśli próbki losowane są z rozkładu symetrycznego

(np. normalnego), to X ∼ −X,

– Metoda daje dobre wyniki, gdy wypłata jest bliska liniowej,

• control variate

– Chcemy obliczyć E(X) (np. cena amerykańskiej opcji call),

– Wiemy o zmiennej losowej Y , której E(Y ) znamy, a która jest w pewnym

sensie bliska E(X) (np. cena europejskiej opcji call)

– Przyjmujemy Z = X − Y + E(Y ), skąd E(Z) = E(X)

– Obliczamy E(X) symulując wartości zmiennej Z.

– Korzyść odniesiemy gdy: Var(Z) = Var(X − Y ) < Var(X).

– Intuicję wynikającą z tej metody można stosować nie tylko dla MC:

∗ z modelu wynikają oszacowania x = I

N

(X), y = I

N

(Y ),

∗ model popełnia błąd ε = y − E(Y ) dla zmiennej Y ,
∗ zakładamy, że dla zmiennej X popełnia ten sam błąd,

∗ lepsze oszacowanie dla X to: x − ε = x − y + E(Y ).

– Metodę można „dostroić”

∗ dla θ ∈ R definiujemy

Z

θ

= X + θ(E(Y ) − Y )

∗ w dalszym ciągu E(Z

θ

) = E(X) oraz

Var(Z

θ

) = Var(X) − 2θ cov(X, Y ) + θ

2

Var(Y )

∗ wariancja jest najmniejsza dla θ

min

=

cov(X,Y )

Var(Y )

,

∗ wielkość cov(X, Y ) zwykle nie jest znana, ale można ją oszacować –

metodą MC (nie jest potrzebna duża dokładność)

• importance sampling (metodę najlepiej wyjaśnić na przykładzie)

– chcemy wycenić opcję call (X = (S

T

− K)

+

),

która jest głęboko out-of-the-money,

– większość wylosowanych próbek da zerową wypłatę, co jest stratą czasu

obliczeniowego,

– trzeba zatem losować tylko te trajektorie, które wniosą wkład do ceny

(zobacz slajd „Rozkłady warunkowe”),

– Niech F (x) = P (S

T

¬ x), G(x) = P (S

T

¬ x|S

T

> K) oraz

k = P (S

T

> K) = 1 − F (K). Wówczas dla x ­ K mamy

G(x) =

P (K<S

T

¬x)

P (S

T

>K)

=

P (S

T

¬x)−P (S

T

¬K)

P (S

T

>K)

=

F (x)−1

k

+ 1.

– Stąd G

−1

(u) = F

−1

(1 + k(u − 1)), u ∈ (0, 1).

– Jeśli Z jest losowane z rozkładu o dystrybuancie G, to

E[X] = P (S

T

¬ K) · E[X|S

T

¬ K] + P (S

T

> K) · E[X|S

T

> K] =

= (1 − k) · 0 + k · E[Z] = k · E[Z].

– Podejście ogólne: zakładamy, że losowana zmienna pochodzi z rozkładu

o gęstości f ,

– Zastąpimy tą gęstość przez g, która powinna być bliska f · h (gęstości

zmiennej h(X), którą całkujemy), wtedy

E

f

(h(X)) =

Z

h(x)f (x) dx =

Z

h(x)

f (x)

g(x)

g(x) dx = E

g

(h(X)

f (X)

g(X)

),

a funkcja podcałkowa w ostatniej całce ma małą wariancję.

Metody symulacji procesów stochastycznych

• symulacja dokładnego rozwiązania,

• symulacja z prawdopodobieństwa przejścia,

• przybliżanie dynamiki,

Symulacja dokładnego rozwiązania,

Stosujemy, gdy znana jest jawna postać silnego rozwiązania SDE, to znaczy daje
się ono zapisać jako funkcjonał procesu Wienera: X(t) = G(t, W

[0,t]

).

• Ustalamy X

0

= x

0

,

• Dla i = 1, . . . , N losujemy trajektorię (W (t

i

) oraz obliczamy

X

i

= G(t

i

, {W (t

1

), . . . , W (t

i

)}).

Przykład: w modelu Blacka-Scholesa

dS(t) = r dt + σ dW (t)

rozwiązanie ma postać:

S(t) = S(0) exp((r − σ

2

/2)t + σW (t)).

Symulacja z prawdopodobieństwa przejścia

Możemy stosować, gdy znamy prawdopodobieństwa przejścia procesu między dwo-
ma dowolnymi momentami czasu.

• Ustalamy X

0

= x

0

.

• Dla i = 1, . . . , N losujemy

X

i

∼ p(t

i

, · ; t

i−1

, X

i−1

)

(rozkład prawdopodobieństwa X

i

w chwili t

i

pod warunkiem, że proces ma

wartość X

i−1

w chwili t

i−1

).

Przykład: w modelu terminowych stóp procentowych Vasicka

dr(t) = α(β − r(t)) dt + σ dW (t),

proces r(t) ma rozkład normalny z warunkową średnią i wariancją

µ(t; s, x) = β + e

−α(t−s)

(x − β),

σ(t; s, x) =

σ

2

2α

(1 − e

−2α(t−s)

).

Przybliżanie dynamiki

Dyskretyzujemy równanie różniczkowe, zamieniając je w równanie różnicowe.

dX(t) = µ(t, X(t)) dt + σ(t, X(t)) dW (t)

(∗)

Schemat Eulera

• Równanie (∗) dyskretyzujemy do postaci

X

i+1

= X

i

+ µ(t

i

, X

i

)∆t + σ(t

i

, X

i

)

√

∆tZ

i

(∗∗)

gdzie Z

i

∼ N (0, 1).

• Ustalamy X

0

= x

0

,

• Dla i = 0, . . . , N − 1 losujemy Z

i

i obliczamy X

i+1

jak w (∗∗).

• Otrzymane w schemacie Eulera rozwiązania są zbieżne do prawdziwego:

∀T > 0, ∃C = C(T ) E( sup

t∈[0,T ]

|X

N

(t) − X(t)|

2

) ¬ C × ∆t.

• przyrosty czasu nie muszą być równe (∆t = t

i+1

− t

i

).

• dla Z

i

= 2B

i

− 1, gdzie B

i

mają rozkład Bernoulliego z p = 1/2, otrzymamy

rozwiązania zbiegające do oryginalnego w sensie słabym, co wystarcza, gdy
chcemy liczyć wartości oczekiwane regularnych funkcjonałów X.

Przykład: Model Fonga-Vasicka zadany jest układem równań:

dr(t) = α(µ − r(t)) dt +

r

v(t) dW

1

t

,

dv(t) = β( ¯

µ − v(t)) dt + σ

r

v(t) dW

2

t

gdzie r jest modelowaną stopą procentową, a v jej ciągłą zmiennością. Dyskretyzacja
ma postać:

r

i+1

= r

i

+ α(µ − r

i

)∆t +

√

v

i

∆tZ

1

i

,

v

i+1

= v

i

+ β( ¯

µ − v

i

)∆t + σ

√

v

i

∆tZ

2

i

gdzie (Z

1

i

, Z

2

i

) ∼ N (0, Σ

i

), a Σ

i

= cov(dW

1

, dW

2

).

Literatura

[1] R. Seydel, Tools for Computational Finance, Springer 2006.

[2] G. Fusai, A. Roncoroni, Implementing Models in Quantative Finance: Me-

thods and Cases, Springer 2008.

[3] P. Jackel, Monte Carlo Methods in Finance, Wiley 2002.

[4] C. P. Robert, G. Casella, Introducing Monte Carlo Methods with R, Springer

2010.