---
title: "Cvičenie 5"
author: "Ján Gašper"
date: "13. 03. 2021"
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.

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

Zvolíme logaritmickú transformáciu premennej $S$:
$$
  x = \ln(S/E)
$$
V pôvodnej rovnici $S \in (0, \infty)$. Po transformácii bude $x \in (-\infty, \infty)$. Čas do exspirácie budeme značiť $\tau := T-t$, pričom $\tau \in [0, T]$. Premennú $x$ budeme nazývať *priestorová premenná* a premennú $\tau$ zas *časová premenná*.

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)$.


## Diskretizácia premenných

Numericky však nevieme riešiť pre $x \in (-\infty, \infty)$. Namiesto toho budeme riešiť na dostatočne dlhom symetrickom intervale $x \in [-x_{max}, x_{max}]$:
$$
  \ln(S/E) \in [-x_{max}, x_{max}] \\
  S \in [E \cdot e^{-x_{max}}; E \cdot e^{x_{max}}]
$$
Ak napríklad $x_{max} = 2$, tak
$$
  S \in [0,\!135 E;\ 7,\!389E]
$$
Pre hodnoty $x = -x_{max}, x=x_{max}$ musíme ešte doplniť okrajové podmienky. Nazveme ich $u(-x_{max}, \tau) = \phi(\tau)$ pre $x=-x_{max}$ a $u(x_{max}, \tau) = \psi(\tau)$ pre $x=x_{max}$. Rozmyslite si ich tvar pre *call* opciu. Návod: použite put-call paritu.

Priestorovú premennú $x$ diskretizujeme na $2n$ častí, resp. budeme mať $2n+1$ deliacich bodov. Priestorový krok označíme $h := x_{max}/n$. Časovú premennú diskretizujeme na $m$ častí (t.j. $m+1$ deliacich bodov) s krokom $k := \tau_{max}/m$, kde $\tau_{max} = T$. 


## Odvodenie numerickej schémy

Nahraďte derivácie v rovnici
$$
 \frac{\partial u}{\partial \tau}  - \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2} = 0
$$
pomernými diferenciami. Pre parciálnu deriváciu podľa času použite doprednú pomernú diferenciu. Odvoďte tak predpis:
$$
  u(x, \tau+k) = u(x, \tau) + k\ \frac{\sigma^2}{2} \cdot \frac{u(x-h, \tau) - 2\cdot u(x, \tau) + u(x+h, \tau)}{h^2}
$$
pre $\tau = 0, k, 2k, \ldots, mk$ a $x = -x_{max}+h, -x_{max}+2h, \ldots, 0, \ldots, x_{max}-2h, x_{max}-h$.

Z dôvodov, ktoré budú vysvetlené na numerickom modelovaní musí platiť $k \leq \frac{h^2}{\sigma^2}$, ináč bude zle.



## Ocenenie európskej *call* opcie

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\%$.

```{r}
e <- 50
r <- 0.04
sigma <- 0.4
tau_max <- 1

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

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

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

# 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(-x_max, x_max, by=h)
tt <- seq(0, tau_max, by=k)

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

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

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

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

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


```{r}
c <- k*sigma^2/(2*h^2) # Predrátaná konštanta 

for(j in 2:(m+1)){
  for (i in 2:(2*n)){
    # 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))
}

```

Experimentuje s rozličnými nastaveniami delenia, t.j. rozličnými hodnotami $m$ a $n$ a sledujte ako sa správa chyba.


## Stabilita metódy

Skúste zvoliť časový a priestorový krok tak, že metóda bude nestabilná, t.j. také delenie že $k > \frac{h^2}{2\sigma^2}$.



# Bonusová domáca úloha (5b)

## 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.

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)


## Iná logaritmická transformácia (1b)

Namiesto transformácie $x := \ln(S/E)$ sme mohli použiť transformáciu $x := \ln(S)$. Ako by sa zmenili naše výpočty, ak by sme postupovali takto?


## 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`.
