Metoda Monte Carlo
Rafał Kucharski
kucharski@math.us.edu.pl
www.math.us.edu.pl/rk
Szkoła Letnia Matematyki Finansowej, Tarnów 2012
Metoda Monte Carlo
Rafał Kucharski
kucharski@math.us.edu.pl
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.