---
title: "Cvičenie 5"
author: "Finančné deriváty"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Numerické riešenie Black-Scholesovej rovnice

Budeme riešiť Black-Scholesovu parciálnu diferenciálnu rovnicu numericky pomocou explicitnej metódy. Máme k tomu dva dôvody:
- nie pre každý payoff sa dá Black-Sholesova rovnica riešiť explicite,
- potrebujeme začať budovať aparát na numerické oceňovanie amerických opcií.

Budeme sa zaoberať prípadom, v ktorom podkladové aktívum nevypláca dividendy. Rovnica
$$
  \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
$$
platí pre $S \in (0, \infty)$, a $t \in [0, T]$ a má koncovú podmienku $V(S, T) = payoff(S)$.


## Transformácia premenných

Namiesto rovnice pre $V$ budeme riešiť transformovanú rovnicu pre funkciu $Z(x, \tau)$, pričom $Z(x, \tau) = V(S, t)$, kde $x = \ln(S/E)$ a $\tau = T - t$. V pôvodnej rovnici $S \in (0, \infty)$. Po transformácii bude $x \in (-\infty, \infty)$. Premennú $x$ budeme nazývať *priestorová premenná* a premennú $\tau$ zas *časová premenná*.

Dokončite túto transformáciu (štátnicová otázka):
$$
\begin{align*}
  \frac{\partial V}{\partial t} &= \frac{\partial Z}{\partial \tau} \frac{\partial \tau}{\partial t}
  \\
  &=- \frac{\partial Z}{\partial \tau}
  \\
  &
  \\
  \frac{\partial V}{\partial S} &= \frac{\partial Z}{\partial x} \frac{\partial x}{\partial S}
  \\
  &= \frac{\partial Z}{\partial x} \frac{1}{S}
  \\
  &
  \\
  \frac{\partial^2 V}{\partial S^2} &= \frac{\partial}{\partial S} \left(\frac{\partial V}{\partial S} \right)
  \\
  &= \frac{\partial}{\partial S} \left( \frac{\partial Z}{\partial x} \frac{1}{S} \right)
  \\
  &= \dots
\end{align*}
$$

Transformovanú rovnicu môžeme po prenásobení $\exp(-\alpha x - \beta \tau)$ previesť na rovnicu vedenia tepla. $V(S, t) = \exp(-\alpha x - \beta \tau) u(x, \tau)$, kde $\alpha = \frac{r}{\sigma^2} - \frac12$, $\beta = \frac{r}{2} + \frac{\sigma^2}{8} + \frac{r^2}{2\sigma^2}$. Funkcia $u$ potom spĺňa:
$$
 \frac{\partial u}{\partial \tau}  - \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2} = 0
$$

Rozmyslite si tvar počiatočnej podmienky $u(x,0) = u_0(x)$. Návod: $u(x, \tau) = \exp(\alpha x + \beta \tau) V(S, t)$, kde $S = E e^x$ a $t = T - \tau$.


## Diskretizácia premenných a numerická schéma

Namiesto hodnôt funkcie $u$ pre všetky (spojité) hodnoty $x$ budeme funkciu $u$ počítať iba v bodoch $\ldots, -2, -h, 0, h, 2h, \ldots$. Podobne budeme uvažovať iba časové body $\tau = 0, k, 2k, 3k, \ldots$.

Namiesto parciálnej derivácie podľa času použijeme doprednú diferenciu:
$$
  \frac{\partial u}{\partial \tau} \approx \frac{u(x, \tau + k) - u(x, \tau)}{k}
$$
Všimnime si, že v limite pre $k \to 0^+$ je pravá strana definíciou parciálnej derivácie.

Namiesto druhej parcálnej derivácie podľa priestorovej premennej použijeme diferenciu
$$
  \frac{\partial^2 u}{\partial x^2} \approx \frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2}
$$
Na skúške z matematickej analýzy (1) sa vyskytol príklad, v ktorom bolo treba ukázať, že v limite $h\to 0^+$ pravá strana konverguje ku druhej derivácii.

Dokopy teda máme približnú rovnosť:
$$
  \frac{u(x, \tau + k) - u(x, \tau)}{k} - \frac{\sigma^2}{2}\frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2} \approx 0,
$$
resp.
$$
  u(x, \tau + k) \approx u(x, \tau) + \frac{k \sigma^2}{2}\frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2}.
$$

Veľkosti krokov musia spĺňať tzv. Courant–Friedrichs–Lewyho podmienku:
$$
  k < \frac{h^2}{\sigma^2}.
$$
Ak by táto podmienka nebola splnená, metóda nebude stabilná. Viac o tejto podmienke sa dozviete na predmete Numerické modelovanie.




## Ohraničená oblasť a okrajové podmienky

Numericky však nevieme riešiť pre $x \in (-\infty, \infty)$. Namiesto toho budeme riešiť na dostatočne dlhom symetrickom intervale $x \in [-L, L]$, t.j. v bodoch $x= -L, -L+h, \ldots, -h, 0, h, 2h, \ldots, L-h, L$.

V reči pôvodných premenných $S$ a $E$ to znamená, že rovnicu riešime na oblasti:
$$
  \ln(S/E) \in [-L, L] \\
  S \in [E \cdot e^{-L}; E \cdot e^L]
$$
Ak napríklad $L = 2$, tak
$$
  S \in [0,\!135 E;\ 7,\!389E]
$$

Pre hodnoty $x = \pm L$ musíme ešte doplniť okrajové podmienky. Nazveme ich $u(-L, \tau) = \phi(\tau)$ pre $x=-L$ a $u(L, \tau) = \psi(\tau)$ pre $x=L$. Odvodíme ich tvar:

- Ak $E \to 0$, tak $Call \sim 0$
- Ak $E \to \infty$, tak $Put \to 0$. Potom z put-call parity vyplýva $Call \sim Put + S - E e^{-r \tau}$, resp. $Call \sim S - E e^{-r \tau}$.


### Bonus: Lepšia asymptotická aproximácia (3b)

Pomocou asymptotickej aproximácie $N(x)$ pre $x \to -\infty$ ukážte že lepší odhad ceny *call* opcie pre malé hodnoty $S$ je:
$$
  Call(S) \overset{S \to 0^+}{\sim} \frac{2 \sqrt{2} S \sigma^{3} \tau^{\frac{3}{2}} \exp\left({- \frac{r^{2} \tau}{2 \sigma^{2}} - \frac{r \tau}{2} - \frac{r \ln{\left(\frac{S}{E} \right)}}{\sigma^{2}} - \frac{\sigma^{2} \tau}{8} - \frac{\ln{\left(\frac{S}{E} \right)}^{2}}{2 \sigma^{2} \tau}}\right)}{\sqrt{\pi} \sqrt{\frac{S}{E}} \left(4 r^{2} \tau^{2} + 8 r \tau \ln{\left(\frac{S}{E} \right)} - \sigma^{4} \tau^{2} + 4 \ln{\left(\frac{S}{E} \right)}^{2}\right)},
$$

## Implementácia

Oceňte európsku *call* opciu s exspiračnou cenou $E=50$, ktorá exspiruje o rok. Podkladové aktívum má volatilitu $\sigma = 40\%$ a bezriziková úroková miera je $r = 4\%$.

Pri výpočte sa oplatí predrátať konštantu $\frac{k \sigma^2}{2h^2}$.

```{r}
E <- 50
r <- 0.04
sigma <- 0.4

# Priestorové delenie
L <- 2
n <- 20
h <- L/n

# Časové delenie
tau_max <- 1
m <- 52
k <- tau_max/m

# pomocné konštanty
alpha <- r/(sigma^2) - 0.5
beta <- r/2 + sigma^2 / 8 + r^2/(2*sigma^2)

# Priestorová a časová diskretizácia
xx <- seq(-L, L, by=h)
tt <- seq(0, tau_max, by=k)

# Okrajové podmienky
phi <- function(tau) 0# TODO 
psi <- function(tau) 0# TODO

# Počiatočná podmienka
u0 <- function(x) 0# TODO

# Matica riešení
u <- matrix(0, nrow=2*n + 1, ncol=m+1)

# Okrajové a počiatočné podmienky
u[1, ] <- phi(tt)
u[2*n +1, ] <- psi(tt)
u[, 1] <- u0(xx)

# sol:
# x = -L   # # # # # # # #
#
# x = 0    # # # # # # # #
# 
# X = +L   # # # # # # # #
#          ^             ^
#        tau=0        tau=tau_max

k <= h^2/(sigma^2) # Podmienka stability


gamma <- k*sigma^2/(2*h^2) # Predrátaná konštanta 

for (j in 2:(m+1)) {
  for (i in 2:(2*n)) {
    u[i, j] <- 0 # TODO
  }
}
```



## Spätná transformácia 

Transformujte získané riešenie späť do premenných $S$, $t$. Vykreslite získané riešenie.

```{r}
# TODO
```



## Porovnanie s analytickým riešením

Porovnajte získané numerické riešenie s riešením získaným pomocou vzorca.

```{r}
CallValue <- function(S, E, r, sigma, tau) {
  d1 <- (log(S/E) + (r + 0.5*sigma^2)*tau)/(sigma * sqrt(tau))
  d2 <- (log(S/E) + (r - 0.5*sigma^2)*tau)/(sigma * sqrt(tau))
  
  return(S*pnorm(d1) - E*exp(-r*tau)*pnorm(d2))
}

# TODO
```


## Chybový graf

Vykreslite graf, ktorý bude mať na vodorovnej osi $S$ a na zvislej osi rozdiel $V_{numeric} - V_{analytic}$.

```{r}
#TODO
```


## Bonus: Iná logaritmická transformácia (2b)

Namiesto transformácie $x := \ln(S/E)$ sme mohli použiť transformáciu $x := \ln(S)$. Podrobne popíšte čo všetko by sa zmenilo v implementácii.



## Bonus: Put opcia numericky (3b)

Naprogramujte funkciu, ktorá numericky ocení európsku *put* opciu, ak akcia vypláca dividendy. Signatúra funkcie nech má tvar: `PutValueNumeric <- function(S, E, r, sigma, tau, q=0)` Funkcia musí zbehnúť aj ak je argument `S` vektor. Možno bude potrebné použiť interpoláciu.

Zvoľte vhodné parametre delenia $m$, $n$ tak aby metóda bola stabilná a dostatočne presná.

Pomôžte si [prednáškou](http://www.iam.fmph.uniba.sk/institute/stehlikova/fd17/06_numerika_euro.pdf)



## Bonus: Iný názov premennej (1b)

V kóde sme čas do exspirácie nazvali `tau_max` namiesto (na prvý pohľad jednoduchšieho) označenia `T`. Vysvetlite prečo nie je vhodné nazvať premennú s časom do exspirácie (ani žiadnu inú premennú či funkciu) ako `T`.


## Bonus: Metóda binárneho stromu (3b)

Odvoďte iteračnú schému, ak sa použije krok $h = \sigma \sqrt k$. Vysvetlite ako táto voľba kroku zjednoduší výpočet.

Implementujte funkciu `CallValueNumeric(S, E, r, sigma, tau)` pre výpočet ceny *call* opcie, ak je dnešná cena podkladového aktíva `S`. Stačí ak výpočet zbehne pre skalárnu hodnotu `S` (teda nie vektor).


## Bonus: Experimenty so stabilitou (2b)

Skúste zvoliť také veľkosti krokov, že:

- $k < \frac{h^2}{2\sigma^2}$
- $k = \frac{h^2}{2\sigma^2}$
- $\frac{h^2}{2\sigma^2} < k < \frac{h^2}{\sigma^2}$
- $k = \frac{h^2}{\sigma^2}$
- $k > \frac{h^2}{\sigma^2}$

Vykreslite grafy cien opcií, ktoré vzniknú takýmito metódami.

# Odovzdanie bonusov

Bonusové úlohy odovzdajte do 6 dní mailom s predmetom **FD24 – nick**, kde **nick** je vami zvolená prezývka skupiny.