Tekst źródłowy

[[Category:MKZR]]
===Wycena opcji: modele z czasem ciągłym===

====Model Blacka-Scholesa dla europejskiej opcji kupna====

Model Blacka-Scholesa zakłada że, że cena akcji <math>\displaystyle S(t)</math> ewoluuje zgodnie z  geometrycznym ruchem Browna przy ciągłej kapitalizacji z roczną stopą procentową r:

:<math>dS(t) = r S(t) dt + \sigma S(t) d W(t)\,</math>.

W chwili <math>t_0</math> instument bazowy ma wartość <math>S0</math> a cena wykupu (lub sprzedaży) w chwili T wynosi S(T). Pytamy się o wartość instrumentu pochodnego w dowonej chwili <math>t<T</math>. Zakładając, że chcemy utrzymać portfel bez ryzyka (delta neutralny) można pokazać, że wartość taka jest dana przez równanie Blacka-Scholesa  

:<math> \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0. </math>

gdzie r jest stopą procentową wolną od ryzyka. 

Równanie to posiada analityczne rozwiązanie i jest to słynny wzór Blacka-Scholesa dla ceny opcji kupna na europejska akcję bez dywidend:

:<math> C(S,t) = SN(d_1) - Ke^{-r(T - t)}N(d_2) \, </math>

i ceny opcji sprzedaży  na europejska akcję bez dywidend:

:<math> P(S,t) = Ke^{-r(T-t)} - S + C(S,t). \ </math>

gdzie:

:::<math> d_1 = \frac{\ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})(T - t)}{\sigma\sqrt{T - t}} </math>

:::<math> d_2 = d_1 - \sigma\sqrt{T - t}. </math>

====Własności wzorów Blacka-Scholesa ====

Powyższe wzory można zaimplementować jako funkcje w matlab/GNU Octave w następujący sposób:

<source lang='matlab'>
function [C,P] = BlackScholes(S0,K,r,T,sigma)
  d1=(log(S0./K)+(r+sigma^2/2)*T)/(sigma*sqrt(T));
  d2=d1-sigma*sqrt(T);
  C = S0.*normcdf(d1)-K*exp(-r*T)*normcdf(d2);
  P = K*exp(-r*T)*normcdf(-d2)-S0.*normcdf(-d1);
end
</source>

 Warto zwrócić uwagę by ten program został zapisany w pliku o nazwie BlackScholes.m.

[[Plik:BS_formula.png|thumb|360px|Wykres ceny europejskiej opcji sprzedaży dla 
K=51.11,r=0.1,T=5.0 oraz <math>\sigma</math>=0.2 od ceny aktualnej S (linia ciągła). Kropkami zaznaczono funkcje  <math>S-K e^{-r T}</math> ]]

Program ten zwraca dla zadanych wartości ceny aktualnej ("spot") S0, cena rozliczenia opcji K w czasie T oraz wysokość rocznej stopy procentowej wolnej od ryzyka dla terminu wygaśnięcia opcji r i volatility <math>\sigma</math> dwie ceny - opcji kupna C oraz sprzedaży P. Są one zwracane jako dwa argumenty, więc poprawne wywołanie takiej funkcji powinno być np.:
<source lang='matlab'>
octave:8> [call,put]=BlackScholes(31 ,51.1 ,0.1 ,5 , 0.2)
call =  5.4876
put =  5.4813
</source>

Wzór ten ma kilka ciekawych własności. Po pierwsze zauważmy, że gdyby fluktuacje ceny akcji znikły w chwili obserwacji to wartość opcji była by równa różnicy między ceną aktualną a zdyskontowaną ceną wykupu <math>S-K e^{-r T}</math>, gdy tylko taka różnica była by większa od zera. W innym przypadku wartość takiej opcji jest równa zeru.  Na rysunku jest umieszczona zależności ceny wykupu od wartości aktualnej opcji wygenerowany następującym skryptem:

<source lang='matlab'>
clear all
close all
T=5.
S0=31
sigma=0.2
r = 0.1
K = S0*exp(r*T)
S0_tab = 1:1:60;
plot(S0_tab,BlackScholes(S0_tab,K,r,T,sigma), S0_tab, max(S0_tab-K*exp(-r*T),0),'o' );
exact=BlackScholes(S0,K,r,T,sigma)
line([S0,S0],[0,exact])
grid on
</source>


Można przypuszczać, że wykres wzoru Blacka Scholes-a będzie dążył do kropkownej zależności dla <math>\sigma \to 0</math>. Rzeczywiście, można się o tym przekonać  rysując kilka zależności C(S) dla malejących wartości <math>\sigma</math>. Na rysunku otrzymanym przez uruchomienie poniższego programu:

[[Plik:BS_formula2.png|thumb|360px|Wykres ceny europejskiej opcji sprzedaży dla 
K=51.11,r=0.1,T=5.0   od ceny aktualnej S  dla kilku wartości <math>\sigma</math>=0.01,0.1,0.2 i 0.3 ]]

<source lang='matlab'>
clear all
close all
T=5.
S0=31
r = 0.1
K = S0*exp(r*T)
S0_tab = 1:1:60;
for sigma=[0.01 0.1 0.2 0.3] 
  plot(S0_tab,BlackScholes(S0_tab,K,r,T,sigma) );
  hold on
endfor 
grid on
</source>



otrzymujemy cztery krzywe dla coraz mniejszych wartości volatility. 



Rozważmy przypadek, gdy jako cenę aktualna instrumentu bazowego weżmiemy zdyskontowaną wartość opcji terminie wygaśnięcia. W przypadku braku fluktuacji cen instrumentu bazowego w takim przypadku nasze opcję mają zerową wartość. Jednak pojawienie się fluktuacji wyrażonej przez niezerowe <math>\sigma</math> powoduje, że opcje zyskują pewną wartość. Można to zobaczyć z powyższych rysunków lub dokładniej na przykładzie:

<source lang='matlab'>
T=5.0;
sigma=0.2;
r = 0.1;
K = 52.1;
BlackScholes( K*exp(-r*T),K,r,T,sigma)
</source>


==== Symulacje Monte Carlo ceny instrumentu pochodnego ====

Wzory Blacka-Scholesa umożliwiają bardzo szybkie i dokładne obliczenie ceny opcji. Oczywiście w praktyce oczywiście dochodzi element eksperymentalnego wyznaczenia volatility, czyli współczynnika zmienności  ceny instrumentu bazowego. Stosowalność wzorów Blacka-Scholesa jest jednak ograniczona do opcji Europejskich, których przypadku możliwość wykupu/sprzedaży opcji dokładnie w terminie wygaśnięcia umożliwia otrzymanie rozwiązania analitycznego. W innych przypadkach, np. gdy dopuścimy sprzedaż w dowolnej chwili przed wygaśnięciem, podejście analityczne staje się niemożliwe. W takich przypadkach popularną i skuteczna metoda jest symulacja numeryczna przebiegu czasowego instrumentu bazowego i wyliczenie odpowiednich statystyk z otrzymanych trajektorii.

Znamy już algorytm generujący [[MKZR:Numeryczne_rozwi%C4%85zania_r%C3%B3wna%C5%84_stochastycznch-przyk%C5%82ady#Geometryczny_proces_Wienera|geometryczny ruch Browna]], możemy więc wygenerować wiele trajektorii, a zadanym warunku początkowym z zadanym krokiem h:

<source lang='matlab'>
clear all
close all
N=100;
M=100000;
T=5.1;
h=T/N;

S0=55;
sigma=0.4;
K = 50;
r = 0.08;
x=zeros(M,N);
x(:,1)=S0*ones(M,1); 
for i=2:N
  x(:,i)=x(:,i-1) + r*x(:,i-1)*h + sigma*sqrt(h)*x(:,i-1).*normrnd (0,1,M,1);
endfor
</source>

W powyższym kodzie macierz x będzie zawierała 100tys. trajektorii geometrycznego procesu Wienera z zadanymi parametrami. Czas potrzebny na wygenerowanie takiej ilości danych nie powinien przekraczać kilku sekund. Ważne jest by w podobnych przypadkach używając systemu matlab/GNU Octave stosować kod wektorowy tak jak w powyższym przykladzie. 
Mając te trajektorie możemy obliczyć średnie z  max(S(T)-K,0), lub max(K-S(T),0) dla opcji odpowiednio kupna lub sprzedaży:
<source lang='matlab'>
call_MC=exp(-r*T)*mean( max(x(:,N)-K,0) )
put_MC=exp(-r*T)*mean( max(K-x(:,N),0) )
</source>

i następnie porównać je z wartościami otrzymanymi ze wzroru Blacka-Scholesa:

<source lang='matlab'>
[call,put]=BlackScholes(S0,K,r,T,sigma)
</source>.

Możemy też wyliczyć względny błąc jaki dostajemy w przypadku posługiwania się symulacją stochastyczną. 

<source lang='matlab'>
error_rel_call=abs(call_MC-call)/call_MC*100
error_rel_put=abs(put_MC-put)/put_MC*100
</source>

W badanym przypadku błąd ten będzie oscylował w okolicach jednego procenta dla 100tys trajektorii. Daje to szacunkową orientację w możliwościach metod stochastycznych na współczesnych komputerach: na jednym komputerze uzyskujemy jedną wycenę w ciągu jednej sekundy z dokładnością jednego procenta. W pewnych przypadkach może być to wystarczająco szybko, dla pewnych zastosowań analitycznych będzie to jednak za wolno i w takich przypadkach stosuje się potężne komputery równoległe, które są w stanie przyśpieszyć ten proces o tysiące razy.



Równanie stochastyczne dla procesu losowego opisującego cenę instrumentu bazowego w modelu Blacka-Scholesa jest równaniem liniowym i można je scałkować analitycznie do czasu T:

<math>S(T)=S(0) e^{(r-\frac{1}{2}\sigma^2)}T+\sigma \int_0^T dW</math>

gdzie wartości całki z procesu Wienera są zmiennymi losowymi o rozkładnie normalnym <math>N(0,\sqrt(T))</math>. Możemy więc napisać program, który będzie symulować wartość S(T) w jednym i wielu krokach:


<source lang='matlab'>
clear all
close all
N=303;
M=32313;
T=13.1;
h=T/N;

S0=11;
sigma=0.4;
K = 50;
r = 0.1;

x=S0*ones(M,1); 

# integrate in one step  
y = x .* exp(  (r-0.5*sigma^2)*T + sigma.*sqrt(T).*normrnd (0,1,M,1) ) ;

# integrate in N  steps mu=>r
for i=2:N
  x =x + r*x*h + sigma*sqrt(h)*x.*normrnd (0,1,M,1);
endfor

# Call option value
Nsteps=exp(-r*T)*mean( max(x-K,0) )
onestep=exp(-r*T)*mean( max(y-K,0) )
exact=BlackScholes(S0,K,r,T,sigma)
</source>

Program ten wyliczy trzema sposobami tą samą liczbę. Można by się zapytać po co stosować czasochłonną symulację podczas gdy to samo można otrzymać wykonując całkowanie w jednym kroku? Odpowiedz jest prosta - po pierwsze, takie całkowanie nie daje nam całej historii pomiędzy chwilą początkową a czasem wygaśnięcia opcji, a od tej historii może zależeć wycana opcji. Po drugie, scałkowanie analityczne jest możliwe wyłącznie w przypadku, gdy 
proces losowy jest procesem liniowym. W innych przypadkach nawet jeśli interesuje nas tylko wartość w chwili końcowej musimy zastosować metodę numeryczna np. schemat  Eulera.

[[Plik:BS_Sim3_vis1.png|thumb|360px|U góry: Symulacja Monte-Carlo instrumentu bazowego (12 trajektorii). Na dole:zdyskontowane wartości pojedyńczej trajektorii symulowanego instrumentu bazowego oraz wartość  K]]
Poniżej znajduje się program, który graficznie ilustruje proces symulacji Monte-Carlo cen instrumentu bazowego wraz z ceną wykupu K, w różnych układach współrzędnych: normalnym i zdyskontowanym.

<source lang='matlab'>
clear all
close all
N=303;
M=11222;
T=5.
h=T/N;
t=linspace(0,T,N);

S0=31
sigma=0.2
r = 0.1
K = 50
K = S0*exp(r*T)

# integrate in N  steps mu=>r
x=zeros(M,N);
x(:,1)=S0*ones(M,1); # log(1)=0
for i=2:N
  x(:,i)=x(:,i-1) + r*x(:,i-1)*h + sigma*sqrt(h)*x(:,i-1).*normrnd (0,1,M,1);
endfor

# Call option value
Nsteps=exp(-r*T)*mean( max(x(:,N)-K,0) )
exact=BlackScholes(S0,K,r,T,sigma)

# Visualization
figure(1)
subplot(2,1,1)
plot(t,x(1:12,:),'r-',t,K*ones(N,1),'g',t,S0*exp(t*r),'b',t,mean(x),'g.')
axis([0 t(N) 10 100]);
grid on

subplot(2,1,2)
plot(t,x(1,:).*exp(-t.*r),'r',t,S0*ones(N,1),'.')
axis([0 t(N) 0 50]);
grid on
print("BS_Sim3_vis1.eps","-F:24")
#
figure(2)
S0_tab = 1:1:60;
plot(S0_tab,BlackScholes(S0_tab,K,r,T,sigma), S0_tab, max(S0_tab-K*exp(-r*T),0),'o' );
l1=sprintf("Black Scholes option value=%1.1f",exact)
text(37,0,l1)
line([S0,S0],[-30,exact])
grid on
print("BS_Sim3_vis2.eps","-F:24")

</source>


=== Wycena opcji: modele z czasem dyskretnym ===

Modele z czasem dyskretnym oparte najczęściej o [http://pl.wikipedia.org/wiki/Drzewo_dwumianowe_(ekonomia) drzewo dwumianowe] zwane też modelem Cox'a-Ross'a-Rubinstein'a są historycznie pierwszym numerycznym podejściem do wyceny instrumentów pochodnych. Obecnie z uwagi na większą moc obliczeniową komputerów są rzadko stosowane i większą popularnością cieszy się podejście z czasem ciągłym. 

W tych modelach zakładamy, że czas zmienia się o stały przyrost i po każdej zmianie wartość instrumentu <math>S_0</math> przechodzi na wyższą wartość  <math>S_0 u</math> w prawdopodobieństwem <math>p_u</math> oraz na niższą  <math>S_0 d</math> w prawdopodobieństwem <math>p_d</math> (patrz rysunek).
[[Plik:CRR-02.png|thumb|260px|Drzewo możliwości przy wycenie opcji metodą z czaem dykretnym ]]

Poniższy program w matlab/GNU Octave implementuje opisaną regułę:
<source lang='matlab'>
function [price, lattice] = LatticeEurCall(S0,K,r,T,sigma,N)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for i=0:N
   lattice(i+1,N+1)=max(0 , S0*(u^i)*(d^(N-i)) - K);
end
for j=N-1:-1:0
   for i=0:j
      lattice(i+1,j+1) = exp(-r*deltaT) * ...
         (p * lattice(i+2,j+2) + (1-p) * lattice(i+1,j+2));
   end   
end
price = lattice(1,1);
</source>