Wycena opcji europejskich put oraz call metodÄ… Monte Carlo
Jednym z podstawowych zagadnień współczesnych nauk ekonomicznych jest
wycena produktów finansowych i ubezpieczeniowych. Poprawna wycena pozwala nie tylko
na osi gni cie przewagi konkurencyjnej ale również na skuteczne zarz dzanie ryzykiem.
W poniższej pracy przedstawiona zostanie jedna z metod wyceny opcji europejskich  metoda
Monte Carlo.
W praktyce stosuje si nast puj ce metody wyceny instrumentów pochodnych [4]:
" metoda analityczna - w tej grupie dominuj metody oparte na cz stkowych
równaniach różniczkowych oraz metody martyngałowe,
" metody numeryczne, m.in. numeryczne metody rozwi zywania cz stkowych
równań różniczkowych oraz metody symulacji Monte Carlo.
Skoncentrujemy si na metodach symulacyjnych. Ponieważ ewolucja cen akcji ma
natur losow , modelowana jest stochastycznymi równaniami różniczkowymi. Najprostszym
modelem ewolucji cen jest geometryczny ruch Browna [3], [4]:
dS = µSdt +ÃSdW (1)
gdzie:
S - cena akcji,
µ - dryf,
à - zmienność,
W - proces Wienera.
Stochastyczne równanie różniczkowe (1) b dzie rozumiane w sensie Ito (wybór
parametryzacji nie wpływa na wycen opcji [1]). Jednym z głównych rezultatów
martyngałowej wyceny jest formuła [3]:
C0 = EQ[e-rT Payoff ] (2)
gdzie:
C0 - cena opcji,
r - stopa (intensywność) oprocentowania instrumentu wolnego od ryzyka,
T - czas do wykonania opcji,
Payoff - profil wypłaty z opcji,
EQ - wartość oczekiwana obliczana wzgl dem miary martyngałowej (arbitrażowej).
Korzystaj z poj ć teorii miary i teorii procesów stochastycznych można pokazać, że
stochastyczne równanie różniczkowe (1) może być zapisane alternatywnie w postaci [3]:
Q
dS = rSdt + ÃSdW (3)
µ - r
Q
Proces W = W + t jest procesem Wienera wzgl dem miary martyngałowej
Ã
Q zdefiniowanej za pomoc pochodnej Radona-Nikodyma [3]:
dQT = ZT dPT (4)
gdzie:
PT - miara subiektywna (rzeczywista),
2
îÅ‚ Å‚Å‚
µ - r 1 µ - r
ëÅ‚ öÅ‚
ZT = expïÅ‚- WT - ìÅ‚ ÷Å‚ śł
T .
à 2 Ã
íÅ‚ Å‚Å‚
ïÅ‚ śł
ðÅ‚ ûÅ‚
Z formuł (2) i (3) wynika, że cen opcji można wyznaczyć losuj c (symuluj c)
zdyskontowane wartości profilu wypłaty opcji a nast pnie obliczaj c średni po wszystkich
otrzymanych realizacjach. Ponieważ rozwi zanie (silnym) równania (3) jest [3], [4]:
îÅ‚ 1 Å‚Å‚
öÅ‚T
2
ST = S0 expïÅ‚ëÅ‚r - à + ÃWT śł (5)
ìÅ‚ ÷Å‚
2
íÅ‚ Å‚Å‚
ðÅ‚ ûÅ‚
wystarczy wygenerować wartości WT z rozkładu N(0, T ), wstawić do (5). Otrzymany
rezultat wstawić do profilu wypłaty opcji, zdyskontować i obliczyć średni z otrzymanych
wartości.
W poniższej pracy pójdziemy inn drog : drog symulacji stochastycznych równań
różniczkowych (3). Podejście to jest bardziej czasochłonne od opisanego wyżej, jednak może
być bezpośrednio zastosowane do wyceny opcji egzotycznych.
W pracy zastosowano dwie dyskretyzacje stochastycznego równania różniczkowego
(3): algorytm Eulera-Maruyamy oraz algorytm Milsteina.
Algorytm Eulera-Maruyamy ma postać [2]:
"S = rS"t +ÃS"W (6)
gdzie:
"S = Si+1 - Si - przyrost ceny akcji w czasie "t ,
dist
"W = µ "t - przyrost procesu Wienera w czasie "t (µ : N(0,1).
Rozwi zanie numeryczne ma postać:
T
-1
"t
ST = S0 + (7)
""Si
i=0
lub zapisuj c (6) w postaci:
Si+1
= 1+ r"t + Ã"Wi+1 (8)
Si
znajdujemy:
T T
-1 -1
"t
Si+1 "t
ST = S0 = S0 (9)
" "(1+ r"t + Ã"Wi+1)
Si i=0
i=0
W zał czonych skryptach Matlaba zastosowano reprezentacj (9).
Algorytmu Milsteina dla geometrycznego ruchu Browna ma postać [2]:
Si+1 1
2
2
= 1+ r"t + Ã"Wi+1 + Ã [("Wi+1) -1]"t (10)
Si 2
T T
-1 -1
"t
Si+1 "t
2 öÅ‚
2
ST = S0 = S0 [("Wi+1) -1]"t (11)
ìÅ‚ ÷Å‚
" "ëÅ‚1+ r"t + Ã"Wi+1 + 1 Ã
Si i=0 íÅ‚ 2
Å‚Å‚
i=0
WielkoÅ›ci losowe µ pochodz ce z rozkÅ‚adu N(0,1) generowane s za pomoc
funkcji Matlaba randn().
W zał czniku znajduj si cztery skrypty obliczaj ce ceny opcji call z cen
wykonania K . Dwa z nich stosuj algorytm Eulera-Maruyamy i dwa algorytm Milsteina.
Ponadto W każdej z tych dwóch grup zastosowano kod szybki i wolniejszy od poprzedniego
(szybkiego). Zwi zane to jest ze sposobem w jakim Matlab zarz dza pami ci komputera.
Skrypty oprócz oszacowań ceny, podaj jej średni bł d standardowy. Ponadto
weryfikowany jest parytet kupna  sprzedaży oraz porównywane s otrzymane droga
symulacji ceny z wynikami dokładnymi (formuła Blacka-Scholesa).
Skrypty pozwalaj osobie nie posiadaj cej zaawansowanych umiej tności obsługi
Matlaba na wycen opcji call i put dla dowolnych wartości: ceny wykonania, terminu
wykonania, stopy wolnej od ryzyka oraz zmienności. Użytkownik z niewielkimi
umiej tnościami bez problemu zmodyfikuje kod na potrzeby wyceny kombinacji opcji lub
wycen instrumentu z  egzotycznym profilem wypłaty.
Literatura:
1. Czernik T. Skazani na formalizm Ito?, Metody matematyczne,
ekonometryczne i informatyczne w finansach i ubezpieczeniach, AE
Katowice, 2007,
2. Kloeden P.E., Platen E., Numerical solution of stochastic differential
equations, Springer, 1995,
3. Shiryaev A.N., Essentials of stochastic finance. Facts, models, theory.,
World Scientific, 2001,
4. Wilmott P., Paul Wilmott on quantitative finance, Wiley, 1998,
Listingi skryptów
Algorytm Eulera  Maruyamy (duża liczba realizacji)
% wycena opcji europejskich typu call i put metod symulacji Monte Carlo
% dynamika instrumentu bazowego S jest geometrycznym ruchem Browna:
% dS=S*(mi*dt+sigma*dW) gdzie:
% S - cena instrumentu bazowego
% mi - dryf (drift); w algorytmie została zast piona stopa woln od ryzyka r (risk neutral pricing)
% sigma - zmienność (volatility)
% W - proces Wienera
clear
tic % start zegara odmierzaj cego czas obliczeń
n=100000; % liczba symulowanych trajektorii/scenariuszy
S0=1; % aktualna cena instrumentu bazowego
K=1.1; % cena wykonania opcji
r=0.1; % stopa wolna od ryzyka
sigma=0.2; % zmienność
T=1; % termin wykonania
deltaT=0.001; % przyrost czasu - skończona aproksymacja przyrostu dt
m=T/deltaT; % ilość kroków - ilość przedziałów na które dzielimy horyzont T
deltaS=ones(1,n); % inicjalizacja algorytmu - deltaS jest wzgl dnym przyrostem ceny instrumentu bazowego
for k=1:1:m % p tla generuj ca wzgl dne przyrosty ceny instrumentu bazowego - algorytm Eulera
A=randn(1,n);
deltaS=deltaS.*(1+r*deltaT+sigma*A*sqrt(deltaT));
end
cenaakcji=S0*deltaS; % cena instrumentu bazowego w dniu wykonania opcji
call=(cenaakcji>K); % zwraca 1 gdy wartość opcji call w dniu jej wykonania jest dodatnia, zero w przeciwnym przypadku
put=(cenaakcjiPCall=(1/n)*exp(-r*T)*sum((cenaakcji-K).*call) % estymator/oszacowanie ceny opcji call
odch_stand_call=(1/sqrt(n))*exp(-r*T)*std((cenaakcji-K).*call) % średni bł d oszacowania ceny opcji call
PPut=(1/n)*exp(-r*T)*sum((K-cenaakcji).*put) % estymator/oszacowanie ceny opcji put
odch_stand_put=(1/sqrt(n))*exp(-r*T)*std((K-cenaakcji).*put) % średni bł d oszacowania ceny opcji put
parytet=PCall-PPut-S0+K*exp(-r*T) % sprawdzenie parytetu kupna-sprzedaży; teoretyczna wartość wynosi zero
% ponieważ algorytm szacuje ceny opcji, obliczona wartość
% powinna być zbliżona do zera
d1=(log(S0/K)+(r+sigma*sigma/2)*T)/(sigma*sqrt(T)); % d1 - wyst puje w formule Blacka-Scholesa
d2=d1-sigma*sqrt(T); % d2 - wyst puje w formule Blacka-Scholesa
cena_analityczna_call=S0*normcdf(d1,0,1)-K*exp(-r*T)*normcdf(d2,0,1) % dokładna cena opcji call (model Blacka-Scholesa)
cena_analityczna_put=-S0*normcdf(-d1,0,1)+K*exp(-r*T)*normcdf(-d2,0,1) % dokładna cena opcji put (model Blacka-Scholesa)
wzgledny_blad_call=abs(PCall-cena_analityczna_call)/cena_analityczna_call % wzgl dne różnice cen otrzymanych na drodze symulacji oraz
wzgledny_blad_put=abs(PPut-cena_analityczna_put)/cena_analityczna_put % ze wzoru Blacka-Scholesa
czas_obliczen=toc % czas wykonywania skryptu/symulacji w sekundach
clear
Algorytm Milsteina (duża liczba realizacji)
% wycena opcji europejskich typu call i put metod symulacji Monte Carlo
% dynamika instrumentu bazowego S jest geometrycznym ruchem Browna:
% dS=S*(mi*dt+sigma*dW) gdzie:
% S - cena instrumentu bazowego
% mi - dryf (drift); w algorytmie została zast piona stopa woln od ryzyka r (risk neutral pricing)
% sigma - zmienność (volatility)
% W - proces Wienera
clear
tic % start zegara odmierzaj cego czas obliczeń
n=100000; % liczba symulowanych trajektorii/scenariuszy
S0=1; % aktualna cena instrumentu bazowego
K=1.1; % cena wykonania opcji
r=0.1; % stopa wolna od ryzyka
sigma=0.2; % zmienność
T=1; % termin wykonania
deltaT=0.001; % przyrost czasu - skończona aproksymacja przyrostu dt
m=T/deltaT; % ilość kroków - ilość przedziałów na które dzielimy horyzont T
deltaS=ones(1,n); % inicjalizacja algorytmu - deltaS jest wzgl dnym przyrostem ceny instrumentu bazowego
for k=1:1:m % p tla generuj ca wzgl dne przyrosty ceny instrumentu bazowego - algorytm Milsteina
A=randn(1,n);
deltaS=deltaS.*(1+r*deltaT+sigma*A*sqrt(deltaT)+0.5*sigma^2*(A.*A-1)*deltaT));
end
cenaakcji=S0*deltaS; % cena instrumentu bazowego w dniu wykonania opcji
call=(cenaakcji>K); % zwraca 1 gdy wartość opcji call w dniu jej wykonania jest dodatnia, zero w przeciwnym przypadku
put=(cenaakcjiPCall=(1/n)*exp(-r*T)*sum((cenaakcji-K).*call) % estymator/oszacowanie ceny opcji call
odch_stand_call=(1/sqrt(n))*exp(-r*T)*std((cenaakcji-K).*call) % średni bł d oszacowania ceny opcji call
PPut=(1/n)*exp(-r*T)*sum((K-cenaakcji).*put) % estymator/oszacowanie ceny opcji put
odch_stand_put=(1/sqrt(n))*exp(-r*T)*std((K-cenaakcji).*put) % średni bł d oszacowania ceny opcji put
parytet=PCall-PPut-S0+K*exp(-r*T) % sprawdzenie parytetu kupna-sprzedaży; teoretyczna wartość wynosi zero
% ponieważ algorytm szacuje ceny opcji, obliczona wartość
% powinna być zbliżona do zera
d1=(log(S0/K)+(r+sigma*sigma/2)*T)/(sigma*sqrt(T)); % d1 - wyst puje w formule Blacka-Scholesa
d2=d1-sigma*sqrt(T); % d2 - wyst puje w formule Blacka-Scholesa
cena_analityczna_call=S0*normcdf(d1,0,1)-K*exp(-r*T)*normcdf(d2,0,1) % dokładna cena opcji call (model Blacka-Scholesa)
cena_analityczna_put=-S0*normcdf(-d1,0,1)+K*exp(-r*T)*normcdf(-d2,0,1) % dokładna cena opcji put (model Blacka-Scholesa)
wzgledny_blad_call=abs(PCall-cena_analityczna_call)/cena_analityczna_call % wzgl dne różnice cen otrzymanych na drodze symulacji oraz
wzgledny_blad_put=abs(PPut-cena_analityczna_put)/cena_analityczna_put % ze wzoru Blacka-Scholesa
czas_obliczen=toc % czas wykonywania skryptu/symulacji w sekundach
clear
Algorytm Eulera  Maruyamy (szybki)
% wycena opcji europejskich typu call i put metod symulacji Monte Carlo
% dynamika instrumentu bazowego S jest geometrycznym ruchem Browna:
% dS=S*(mi*dt+sigma*dW) gdzie:
% S - cena instrumentu bazowego
% mi - dryf (drift); w algorytmie została zast piona stopa woln od ryzyka r (risk neutral pricing)
% sigma - zmienność (volatility)
% W - proces Wienera
% UWAGA: z uwagi na wysoce zwektoryzowany kod, liczba realicji/symulacji
% nie powinna być duża (OUT of MEMORY)!!!
clear
tic % start zegara odmierzaj cego czas obliczeń
n=20000; % liczba symulowanych trajektorii/scenariuszy
S0=1; % aktualna cena instrumentu bazowego
K=1.1; % cena wykonania opcji
r=0.1; % stopa wolna od ryzyka
sigma=0.2; % zmienność
T=1; % termin wykonania
deltaT=0.001; % przyrost czasu - skończona aproksymacja przyrostu dt
m=T/deltaT; % ilość kroków - ilość przedziałów na które dzielimy horyzont T
A=randn(m,n); % wygenerowanie macierzy licz losowych z rozkładu N(0,1)
deltaS=1+r*deltaT+sigma*A*sqrt(deltaT); % wygenerowanie przyrostów ceny akcji; algorytm Eulera
cenaakcji=S0*prod(deltaS); % cena instrumentu bazowego w dniu wykonania opcji
call=(cenaakcji>K); % zwraca 1 gdy wartość opcji call w dniu jej wykonania jest dodatnia,
% zero w przecinym przypadku
put=(cenaakcji % zero w przecinym przypadku
PCall=(1/n)*exp(-r*T)*sum((cenaakcji-K).*call) % estymator/oszacowanie ceny opcji call
odch_stand_call=(1/sqrt(n))*exp(-r*T)*std((cenaakcji-K).*call) % średni bł d oszacowania ceny opcji call
PPut=(1/n)*exp(-r*T)*sum((K-cenaakcji).*put) % estymator/oszacowanie ceny opcji put
odch_stand_put=(1/sqrt(n))*exp(-r*T)*std((K-cenaakcji).*put) % średni bł d oszacowania ceny opcji put
parytet=PCall-PPut-S0+K*exp(-r*T) % sprawdzenie parytetu kupna-sprzedaży; teoretyczna wartość wynosi zero
 % ponieważ algorytm szacuje ceny opcji, obliczona wartość powinna być zbliżona
% do zera
d1=(log(S0/K)+(r+sigma*sigma/2)*T)/(sigma*sqrt(T)); % d1 - wyst puje w formule Blacka-Scholesa
d2=d1-sigma*sqrt(T); % d2 - wyst puje w formule Blacka-Scholesa
cena_analityczna_call=S0*normcdf(d1,0,1)-K*exp(-r*T)*normcdf(d2,0,1) % dokładna cena opcji call (model Blacka-Scholesa)
cena_analityczna_put=-S0*normcdf(-d1,0,1)+K*exp(-r*T)*normcdf(-d2,0,1) % dokładna cena opcji put (model Blacka-Scholesa)
wzgledny_blad_call=abs(PCall-cena_analityczna_call)/cena_analityczna_call % wzgl dne różnice cen otrzymanych na drodze symulacji oraz
wzgledny_blad_put=abs(PPut-cena_analityczna_put)/cena_analityczna_put % ze wzoru Blacka-Scholesa
czas_obliczen=toc % czas wykonywania skryptu/symulacji w sekundach
clear
Algorytm Milsteina (szybki)
% wycena opcji europejskich typu call i put metod symulacji Monte Carlo
% dynamika instrumentu bazowego S jest geometrycznym ruchem Browna:
% dS=S*(mi*dt+sigma*dW) gdzie:
% S - cena instrumentu bazowego
% mi - dryf (drift); w algorytmie została zast piona stopa woln od ryzyka r (risk neutral pricing)
% sigma - zmienność (volatility)
% W - proces Wienera
% UWAGA: z uwagi na wysoce zwektoryzowany kod, liczba realicji/symulacji
% nie powinna być duża (OUT of MEMORY)!!!
clear
tic % start zegara odmierzaj cego czas obliczeń
n=20000; % liczba symulowanych trajektorii/scenariuszy
S0=1; % aktualna cena instrumentu bazowego
K=1.1; % cena wykonania opcji
r=0.1; % stopa wolna od ryzyka
sigma=0.2; % zmienność
T=1; % termin wykonania
deltaT=0.001; % przyrost czasu - skończona aproksymacja przyrostu dt
m=T/deltaT; % ilość kroków - ilość przedziałów na które dzielimy horyzont T
A=randn(m,n); % wygenerowanie macierzy licz losowych z rozkładu N(0,1)
deltaS=1+r*deltaT+sigma*A*sqrt(deltaT)+0.5*sigma^2*(A.*A-1)*deltaT; % wygenerowanie przyrostów ceny akcji; algorytm Milsteina
cenaakcji=S0*prod(deltaS); % cena instrumentu bazowego w dniu wykonania opcji
call=(cenaakcji>K); % zwraca 1 gdy wartość opcji call w dniu jej wykonania jest dodatnia,
% zero w przecinym przypadku
put=(cenaakcji % zero w przecinym przypadku
PCall=(1/n)*exp(-r*T)*sum((cenaakcji-K).*call) % estymator/oszacowanie ceny opcji call
odch_stand_call=(1/sqrt(n))*exp(-r*T)*std((cenaakcji-K).*call) % średni bł d oszacowania ceny opcji call
PPut=(1/n)*exp(-r*T)*sum((K-cenaakcji).*put) % estymator/oszacowanie ceny opcji put
odch_stand_put=(1/sqrt(n))*exp(-r*T)*std((K-cenaakcji).*put) % średni bł d oszacowania ceny opcji put
parytet=PCall-PPut-S0+K*exp(-r*T) % sprawdzenie parytetu kupna-sprzedaży; teoretyczna wartość wynosi zero
% ponieważ algorytm szacuje ceny opcji, obliczona wartość powinna być zbliżona
% do zera
d1=(log(S0/K)+(r+sigma*sigma/2)*T)/(sigma*sqrt(T)); % d1 - wyst puje w formule Blacka-Scholesa
d2=d1-sigma*sqrt(T); % d2 - wyst puje w formule Blacka-Scholesa
cena_analityczna_call=S0*normcdf(d1,0,1)-K*exp(-r*T)*normcdf(d2,0,1) % dokładna cena opcji call (model Blacka-Scholesa)
cena_analityczna_put=-S0*normcdf(-d1,0,1)+K*exp(-r*T)*normcdf(-d2,0,1) % dokładna cena opcji put (model Blacka-Scholesa)
wzgledny_blad_call=abs(PCall-cena_analityczna_call)/cena_analityczna_call % wzgl dne różnice cen otrzymanych na drodze symulacji oraz
wzgledny_blad_put=abs(PPut-cena_analityczna_put)/cena_analityczna_put % ze wzoru Blacka-Scholesa
czas_obliczen=toc % czas wykonywania skryptu/symulacji w sekundach
clear