:: Prečo numerické riešenie ::

• Načo riešiť rovnicu numericky, keď máme jej explicitné rieenie?
• Explicitné riešenie máme pre európsku call a put opciu. Pre iné deriváty takýto vzorec nemusí existovať. (To je prípad napríklad amerických call a put opcií, ktorými sa ešte budeme zaoberať. A mnohých ďalších.) Riešenie sa však dá nájsť numericky.
• To, že použitie numerických schém vyskúšame najskôr v prípade, v ktorom máme k dispozícii explicitné riešenie, má tú výhodu, že môžeme overiť presnosť získaných numerických výsledkov.

:: Transformácia Black-Scholesovej rovnice ::

• Black-Scholesova rovnica
  
  - je to parabolická PDR, ktorá sa dá transformovať na rovnicu vedenia tepla so zadanou začiatočnou podmienkou.

• Transformácia premenných
  • Máme koncovú podmienku čase T (čas expirácie opcie), nie začiatočnú. To vyriešime transformáciou
    
    To znamená, že novou premennou, namiesto času t, bude čas zostávajúci do expirácie.
  • Premenná S (cena akcie) je kladná. Premennú, ktorá nadobúda všetky reálne hodnoty dostaneme zlogaritmovaním. Nová premenná teda bude
    
    Hodnoty x blízke nule zodpovedajú cenám akcie, ktoré sú blízke expiračnej cene opcie. Záporné hodnoty x predstavujú ceny akcie, ktoré sú nižšie ako je expiračná cena. Analogicky, kladné hodnoty x predstavujú ceny akcie, ktoré sú vyššie ako je expiračná cena. Na transformáciu rovnice by stačila logaritmická transformácia, ale - ako uvidíme neskôr - táto bude výhodnejšia z numerického hľadiska.

• Transformácia na RVT
  • Rovnica, ktorú dostaneme týmito transformáciami, je parabolická rovnica, ale už s konštantnými koeficientami. Túto rovnicu vieme transformovať na rovnicu vedenia tepla transformáciou
    
    kde konštanty určíme tak, aby po tranformácii vznikla práve RVT. Tá správna voľba konštánt je
    
    Vedie k rovnici
    

• Transformácia koncovej podmienky:
  

• Vyriešením tejto RVT a spätnou transformáciou dostaneme explicitné Black-Scholesove vzorce, ktoré sme používali v predchádzajúcich cvičeniach. Teraz však chceme túto rovnicu riešiť numericky.

:: Diskretizácia ::

• Premenná x je z neohraničeného intervalu, nadobúda hodnoty od mínus nekonečna do plus nekonečna. Numericky budeme úlohu riešiť pre x z ohraničeného intervalu [-L, L], kde L bude dosť veľké číslo. Na obrázku je delenie intervalu premennej x, a zodpovedajúce delenie pre cenu akcie, ak je expiračná cena 100 USD.
  
  Všimnime si, že tento interval [-L, L] nemusíme meniť pre opciu s inou expiračnou cenou. O vhodné body S sa postará transformácia x = ln(S/E), ktorú sme použili namiesto jednoduchej logaritmickej transformácie x = ln(S).

  Kvôli tomu, že sme výpočet obmedzili na ohraničený interval, musíme k rovnici dodať okrajové podmienky v krajných bodoch. Krajné body zodpovedajú cenám akcie, ktoré sú veľmi malé, blízke nule (x=-L) a cenám, ktoré sú veľmi veľké a približujú sa k nekonečnu (x=L). Pre takéto limitné hodnoty použijeme aproximácie:

  

• Ďalej potrebujeme diskretrizovať RVT. Sú dve možnosti:
  

  Prvý prístup (explicitná schéma) vyzerá byť jednoduchší. Máme začiatočnú podmienku. Z nej vypočítame hodnoty v nasledujúcej časovej vrstve. Tieto použijeme na výpočet riešenia v ďalšej časovej vrstve, atď. Na konvergenciu metódy je však potrebné splnenie podmienky na vzťah časového a priestorového kroku. Táto podmienka môže prakticky viesť k nevyhnutnosti zvoliť veľmi malý časový krok - kvôli konvergencii metódy, nie kvôli tomu, že by sme riešenie potrebovali v mnohých tak blízkych časových okamihoch.

  V druhom prístupe (implicitná schéma) riešime na každej časovej vrstve sústavu lineárnych rovníc.

:: Implicitná schéma pre call opciu ::

• Zvolíme parametre opcie:

  E = 50;
  r = 0.04;
  D = 0.12;
  sigma = 0.4;
  

• Zvolíme parameter L, určujúci oblasť, na ktorej budeme riešenie počítať:

  L = 2;
  

  a parametre delenia:

  // priestorova premenna
  n = 20;
  h = L/n;
  // casova premenna
  T = 1;
  m = 12;
  k = T/m;
  

• Definujeme konštanty potrebné na transformáciu rovnice:

  alfa = (r-D)/(sigma^2) - 0.5;
  beta = (r+D)/2 + (sigma^2)/8 + ((r-D)^2)/(2*sigma^2);
  

• Pre riešenie trasformovanej rovnice definujeme okrajové podmienky (tri bodky na konci riadku znamenajú, že príkaz sa má chápať akoby bol napísaný na jednom riadku): :

  // x = -L, t.j. cena blizka nule 
  function [phi]=phi(tau)
      phi=0; 
  endfunction 
   
  // x = L, t.j. cena blizka nekonecnu 
  function [psi]=psi(tau)
      psi=E*exp(alfa*L + beta*tau).*(exp(L - D*tau) - exp(-r*tau));
  endfunction
  

  a začiatočnú podmienku:

   // zaciatocna podmienka 
  function [u0]=u0(x)
      u0=E*exp(alfa*x).*max(0, exp(x)-1);
  endfunction
  

  
  
• Vytvoríme maticu, do ktorej budeme vkladať riešenie:

  ries = zeros(2*n + 1,m+1);
  

  a body delenia v čase a v priestore:

  x = -L:h:L;
  tau = 0:k:T;
  

  
  Úloha pre vás: Vložte do matice okrajové podmienky a začiatocnú podmienku.
  
  
• Na výpočet každej časovej vrstvy budeme potrebovať vyriešiť sústavu lineárnych rovníc. Definujeme teraz premenné, ktoré obsahujú hodnoty vystupujúce v trojdiagonálnej matici tejto sústavy:

  a = -0.5*(sigma^2)*k/(h^2);  // pozdlz diagonaly
  b = 1 - 2*a;  // na diagonale
  

  Pri výpočte prvej časovej vrstvy (teda hodnoty v čase k) máme nasledovnú pravú stranu:

  ps = ries(2:2*n,1)';
  ps(1) = ps(1) - a*phi(k);
  ps(2*n-1) = ps(2*n-1) - a*psi(k);
  

  
  Na riešenie systému rovníc použijeme Gauss-Seidelovu metódu:
  

  Môžete použiť funkciu definovanú v súbore gs.m:

  
  // -------------------------------------------------------------------- 
  // Gauss-Seidelova metoda na riesenie symetrickeho trojdiag.o systemu
  // -------------------------------------------------------------------- 
  // 
  // gs(a,b,ps,v0,eps)
  //
  // 3-diagonalna matica systemu ma:
  // a - pod a nad diagonalou 
  // b - na diagonale
  //
  // ps - prava strana systemu - stlpcovy vektor
  // v0 - zaciatocna iteracia - stlpcovy vektor
  //
  // v - aproximacia riesenia (prepisuje sa v jednotlivych iteraciach)
  // eps - kriterium pre ukoncenie iteracii: norm(Av-ps)<=eps 
  
  function [v] = gs(a,b,ps,v0,eps)
  
  // matica A 
  pom = length(v0);
  A(1, 1) = b;
  for i = 2 : pom
      A(i, i) = b;
      A(i, i - 1) = a;
      A(i - 1, i) = a;
  end;
  
  // cyklus bezi, kym nie je splnena podmienka norm(Av-ps)<=eps
  v = v0;
  while norm(A*v - ps) > eps
      v(1) = (ps(1) - a*v(2))/b; 
      for i = 2:pom-1
          v(i) = (ps(i) - a*(v(i-1) + v(i+1)))/b; 
      end
      v(pom) = (ps(pom) - a*v(pom - 1))/b;
  end;
  
  endfunction

  
  Úlohy pre vás:
  • Vypočítajte riešenie na prvej časovej vrstve a vložte ho do matice s riešením (ako začiatočnú aproximáciu riešenia môžete použiť hodnoty z predchádzajúcej časovej vrstvy, zvoľte si kritérium na zastavenie iterácií).
  • Transformujte získané riešenie na riešenie Black-Scholesovej rovnice.
  • Naprogramujte cyklus, v ktorom sa vypočíta riešenie na každej časovej vrstve.
  • Porovnajte s presnými hodnotami z Black-Scholesovho vzorca.
  
  Ukážka výstupu:
  

:: SOR metóda ::

• Sústava rozmeru 2x2 a grafické znázornenie konvergencie pre rôzne hodnoty parametra omega - naprogramujte v Scilabe podobný výpočet, ako sme videli na prednáške:
  
  
  
• Úpravou funkcie gs vytvorte funkciu sor, ktorá bude mať navyše parameter omega. Následne ju použite na výpočet ceny opcie. Jedinou zmenou v hlavnom kóde bude volanie funkcie sor namiesto funkcie gs.
  
  
• Zaznamenávajte počet iterácií potrebných pri výpočte jednej časovej vrstvy re rôzne hodnoty parametra omega. Aká hodnota tohto parametra vyzerá byť najvhodnejšia? Ako to závisí od použitého delenia?

:: Ďalšie príklady na precvičenie ::

1. Explicitná schéma. Naprogramujte explicitnú schému na oceňovanie opcií. Spomínaná podmienka (na zaručenie konvergencie) pre vzťah medzi priestorovým a časovým krokom je
   
   Nazýva sa CLF (Courant-Friedrichs-Lewy) podmienka. Ukážte príklad výpočtu, kedy táto podmienka je splnená a príklad výpočtu, keď splnená nie je.
   
   Ukážka výstupu:
   
   
   
2. Binomický strom. Ak v explicitnej schéme zvolíme
   
   dostaneme nasledovný predpis:
   
   To znamená, že na výpočet určitej hodnoty potrebujeme hodnoty z predchádzajúcej časovej vrstvy v susedných bodoch priestorového delenia. Preto sa tento špeciálny prípad nazýva metódou binomického stromu.
   
   Naprogramujte túto metódu.

---

---