---
title: "Cvičenie 1"
author: "Ján Gašper"
date: "13.02.2023"
output: html_document
editor_options: 
  chunk_output_type: inline
---

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

# Organizácia cvičení

Cvičiaci: Mgr. Ján Gašper, DiS. art.

-   Kancelária M-203
-   Konzultačné hodiny: podľa dohody
-   Kontaktný mail: [jan.gasper\@fmph.uniba.sk](mailto:jan.gasper@fmph.uniba.sk){.email}

Z cvičení možno dostať maximálne 40 bodov.

-   3 domáce úlohy za 7 – 10 bodov (v skupinkách vekosti 1-3 členov, skupinky medzi sebou nekomunikujú)
-   1 písomka za 15 bodov (samostatne, formou „open book")
-   bonusové body (max. 10)
-   dobrovoľné domáce úlohy
-   bonusová aktivita
-   účasť na cvičeniach sa nehodnotí

## R Markdown

Budeme pracovať v R markdown v prostredí RStudio. Markdown umožňuje jednoduché formátovanie textu, vytváranie nadpisov odrážok atď. V RStudiu možno vytvoriť nový takýto dokument pomocou: File → New File → R Markdown. Vyplníte základné údaje o dokumente a vyberiete si formát, v ktorom sa má dokument zobraziť (html, pdf, doc). Dokument dokáže súšťať kód a zobraziť jeho výstupy a zvládne aj $\LaTeX$. Viac o markdowne sa dočítate [v dokumentácii](https://www.markdownguide.org/getting-started/)

Napísaný zdrojový kód treba skompilovať (*uštrikovať*, *knit*) do zvoleného formátu. Kompiláciu začknete tlačítkom s klbkom a ihlicami, alebo skratkou `Ctrl + Shift + K`. Na správne kompilovanie je potrebný balík `knitr`.

Bloky kódu sú oddelené od textu tromi spätnými apostrofmi:

    2 + 2

Ak sa kód má spustiť, za apostrofmi treba vyznačiť jazyk R, v ktorom sa kód spúšťa:

```{r}
2 + 2
```

Takýto blok kódu možno vytvoriť aj pomocou klávesovej skratky `Ctrl + Alt + I`.

## Príručka štýlu

Pri písaní kódu si treba vybrať štýl a držať sa ho. Široko uznávanou príručkou je [tidyverse](https://style.tidyverse.org/). S miernymi modifikáciami ho používa aj [google](https://google.github.io/styleguide/Rguide.html).

Počas týchto cvičení sa budeme pridržiavať verzie google-u. Niekoľko najdôležitejších vecí:

-  názvy premenných sa píšu pomocou `snake_case`, t.j. slová sa píšu malými písmenami, oddelené podčiarkovníkom,
-   názvy funkcií sa píšu formou `BigCamelCase`, t.j. každé slovo (vrátane prvého) sa začína veľkým písmenom,
-   pred čiarkou sa medzera napíše nikdy, za čiarkou vždy: `x[, 1]`,
-   pred a za zátvorkami pri `for`, `if`, `while` sa píše medzera: `while (TRUE) {}`,
-   medzera sa píše aj za zátvorkami v ktorých sú argumenty funkcie: `function(x) {}`,
-   operátory `=`, `<-`, `&&`,  sa oddeľujú medzerami z oboch strán,
-   pri viacčlene sa oddeľujú medzerou znaky `+` a `-`:  `2*x^2 + y^2`,
-   vnútro cyklu alebo podmienky sa zarovnáva o dve medzery:
```
if (y < 0) {
  message("y is negative")
}
```




# Opakovanie: Wienerov proces

Wienerov proces:

- je súbor náhodných premenných,
- začína v nule: $w(0) = 0$,
- má nezávislé prírastky na disjunktných intervaloch: $w_{t_2} - w_{t_1}$ a $w_{t_4} - w_{t_3}$ sú nezávislé pre $t_1 < t_2 \leq t_3 < t_4$,
- má prírastky z normálneho rozdelenia s nulovou strednou hodnotou a disperziou podľa dĺžky intervalu $w_{t_2} - w_{t_1} \sim \mathcal N(0, t_2-t_1)$.



## Príklad: Rozdelenie $w_1 + w_2$

Aké rozdelenie má náhodná premenná $w_1 + w_2$?



## Príklad: simulácia Wienerovho procesu

Naprogramujeme funkciu, ktorá vráti hodnoty Wienerovho procesu v $n$ dátových bodoch vzdialených $\Delta t$:
```{r}
Wiener <- function(n, delta_t) {
  w <- numeric(n)
  
  w[1] <- 0
  for (i in 2:n) {
    w[i] <- w[i-1] + rnorm(1, mean = 0, sd = sqrt(delta_t))
  }
  
  return(w)
}
```
Rozmyslite si, ako táto funkcia funguje.


Vykreslíme viacero trajektórií. Vašou úlohou je opraviť hodnotu `ylim` tak aby koncový bod trajektórie bol na $95\%$ v grafe.

```{r}
n <- 201
delta_t <- 0.005

t <- seq(0, (n - 1)*delta_t, by = delta_t)
ylim <- c(-1, 1)  # FIXME

w <- Wiener(n, delta_t)
plot(t, w, type = "l", lwd = 2, col = rgb(0, 0, 0, 0.9), ylim = ylim)

for (i in 1:50) {
  w <- Wiener(n, delta_t)
  lines(t, w, col = rgb(0, 0, 0, 0.2))
}
```

**Bonus**: Upravte premennú `ylim` tak, aby *celá trajektória* bola na $95\%$ v grafe.



### Bonus: Wienerov proces bez cyklu

Naprogramujte funkciu, ktorá robí to isté, ale neobsahuje for-cyklus. Porovnajte časy potrebné na zbehnutie týchto dvoch metód.

```{r}
tic <- Sys.time()
w <- Wiener(1000000, 0.01)
toc <- Sys.time()

print(round((toc-tic), 5))


WienerLoopless <- function(n, dt) {
  # TODO
}

tic <- Sys.time()
w <- WienerLoopless(1000000, 0.01)
toc <- Sys.time()

print(round((toc-tic), 5))
```



## Príklad: $w_1 + w_2$ numericky

Vygenerujte hodnoty náhodnej premennej $w_1+w_2$ a spravte histogram jej hodnôt. Všimnime si, že nám stačí generovať Wienerov proces v časoch $1$ a $2$, jemnejšie delenie nie je potrebné. Porovnajte pravdepodobnostné rozdelenie tejto náhodnej premennej so získaným histogramom. 

```{r}
# TODO
```



## Príklad: Brownov pohyb (vychýlený Wienerov proces)

Nech $w_t$ je Wienerov proces. Potom proces
$$
  B_t := \mu t + \sigma w_t,
$$
kde $\mu, \sigma$ sú konštanty, nazývame Brownovým pohybom.

Poznámka: Terminológia tu nie je jednoznačná. Niektorí nazývajú Brownovým pohybom Wienerov proces (bez driftu, s jednotkovou volatilitou). Iní zas tvrdia, že Brownov pohyb je fyzikálny fenomén a nie stochastický proces. Pri práci s literatúrou si vždy pozrite značenie.

Naprogramujte funkciu `BrownianMotion(n, delta_t, mu, sigma)` tak, aby vygenerovala trajektóriu Brownovho pohybu.

```{r}
BrownianMotion <- function(n, delta_t, mu, sigma) {
  # TODO
}
```


Podobne ako minule, vykreslíme viac trajektórií do toho istého obrázku. Nastavte vertikálnu os tak, aby koncový bod trajektórie bol na $95\%$ vnútri grafu.
```{r eval = FALSE}
# Po dokončeneí vymažte `eval = FALSE`

n <- 201
delta_t <- 0.005

t <- seq(0, (n - 1)*delta_t, by = delta_t)
ylim <- c(-1, 1)  # FIXME

b <- BrownianMotion(n, delta_t, mu, sigma)
plot(t, b, type = "l", lwd = 2, col = rgb(0, 0, 0, 0.9), ylim = ylim)

for (i in 1:50) {
  w <- BrownianMotion(n, delta_t, mu, sigma)
  lines(t, b, col = rgb(0, 0, 0, 0.2))
}
```


## Príklad: geometrcký Brownov pohyb

Majme proces
$$
  S_t = S_0 \exp(\mu t + \sigma w_t)
$$
Tento proces sa nazýva geometrický Brownov pohyb.

Naprogramujte funkciu, ktorá vráti hodnoty geometrického Brownovho pohybu
```{r eval = FALSE}
# Po dokončeneí vymažte `eval = FALSE`

GeometricBrownianMotion <- function(n, delta_t, mu, sigma, s0) {
  # TODO
}
```
```{r eval = FALSE}
# Po dokončeneí vymažte `eval = FALSE`

n <- 201
delta_t <- 0.005

t <- seq(0, (n - 1)*delta_t, by = delta_t)
ylim <- c(-1, 1)  # FIXME

b <- BrownianMotion(n, delta_t, mu, sigma)
plot(t, b, type = "l", lwd = 2, col = rgb(0, 0, 0, 0.9), ylim = ylim)

for (i in 1:50) {
  w <- GeometricBrownianMotion(n, delta_t, mu, sigma)
  lines(t, b, col = rgb(0, 0, 0, 0.2))
}
```



Poznámka: Existujú dva štandardné tvary geometrického Borwnovho pohybu. Jeden z nich je
$$
  S_t = S_0 \exp(\mu t + \sigma w_t),
$$
ktorý je riešením SDR
$$
  dS_t = \left(\mu + \frac{\sigma^2}{2}\right)S_t \, dt + \sigma S \, dt
$$


Druhý tvar je
$$
  S_t = S_0 \exp\left( \left(\mu - \frac{\sigma^2}{2}  \right)t + \sigma w_t \right)
$$
daný SDR
$$
  dS_t = \mu S \,dt + \sigma S \, dw_t.
$$
  



# Opakovanie: Itôva lema

Itôva lema nám pomôže nájsť tvar stochastickej diferenciálnej rovnice pre transformovaný proces. Ak proces $X$ daný SDR $dX = \mu \, dt + \sigma \, dw_t$ transformujeme dvakrát diferencovateľnou funkciou $g(t, x)$, potom pre proces $Y := g(t, X)$ platí:
$$
  dY = \left(\frac{\partial g}{\partial t} + \mu \frac{\partial g}{\partial x} + \frac{\sigma^2}{2} \frac{\partial^2 g}{\partial x^2} \right)\,dt + \sigma \frac{\partial g}{\partial x}\, dw_t
$$

## Príklad

Majme stochastický proces daný stochastickou diferenciálnou rovnicou: 
$$
  dX_t = 0,\!98\, dt + 0,\!5 \,dw_t
$$
Jedná sa o vychýlený Wienerov proces.

Nájdite stochastickú diferenciálnu rovnicu pre proces $S := e^X$. Procesu $S$ sa hovorí aj *geometrický brownov pohyb*.


### Riešenie:

Funkcia $g$, ktorou transformujeme proces má tvar $g(t, x) = e^x$. Jej parciálne derivácie sú:
$$
\begin{align*}
  \frac{\partial g}{\partial t}     &= 0
  \\
  \frac{\partial g}{\partial x}     &= e^x
  \\
  \frac{\partial^2 g}{\partial x^2} &= e^x
\end{align*}
$$

Teda:
$$
  dS = \left(0 + 0,\!98e^X + \frac{(\frac15)^2}{2} e^X\right)\,dt + \frac15 e^X dw_t
$$
Namiesto $e^X$ budeme písať $S$ a dostaneme:
$$
  \boxed{dY = S \,dt + \frac15 S\, dw_t}
$$


### Iné riešenie:

Použijeme Taylorov polynóm:
$$
  g(t+dt, X+dX) - g(t, x) = \frac{\partial g}{\partial t} dt + \frac{\partial g}{\partial x}dX + \frac12 \frac{\partial^2 g}{\partial x^2}(dX)^2 + o(dt) + o(dx^2)
$$
$$
\begin{align*}
   dS &= 0\,dt + e^X dX + \frac12e^X(dX)^2
   \\
   &= e^X (0,\!98\, dt + 0,\!2 \,dw_t) + \frac12 e^X(0,\!98\, dt + 0,\!2 \,dw_t)^2
   \\
   &= S (0,\!98\, dt + 0,\!2 \,dw_t) + \frac12 S(0,\!98\, dt + 0,\!2 \,dw_t)^2
   \\
   &= S (0,\!98\, dt + 0,\!2 \,dw_t) + \frac12 Y(\underbrace{0,\!9604\, dt^2}_{o(dt)} + \underbrace{0,\!392 \, dt \, dw}_{o(dt)} + 0,\!04 \,\underbrace{dw_t^2}_{dt})
   \\
   &= 0,\!98 Y \, dt + 0,\!2 Y\,dw_t + 0,\!02 Y \,dt
   \\
   &= \boxed{Y \, dt + 0,\!2 Y \, dw_t}
\end{align*}
$$

## Príklad

$$
  dr_t = (2-r_t)\ dt + 3 \sqrt{r_t} \ dw_t
$$
S týmto stochastickým procesom sa ešte stretneme pri modelovaní úrokových mier.

Nájdite stochastickú diferenciálnu rovnicu pre proces $f := r - r_0e^{-t} - 2\left(1-e^{-t}\right)$.

Vo výsledku sa nesmie objaviť hodnota procesu $r$.



### Riešenie (začiatok)

Funkcia, ktorou transformujeme má tvar

$$
\begin{align*}
  g(t, x) &= x - r_0 e^{-t} - 2\left( 1 - e^{-t} \right)
  \\
  &= x  + (2 - r_0)e^{-t} - 2
\end{align*}
$$



## Príklad

Majme stocahstický proces daný stochastickou diferenciálnou rovnicou:

$$
  dX_t = (3X_t - X_t^3)dt + dw_t
$$
Nájdite diferenciál pre proces $Y := X^2$



## Numerické riešenie SDR (Euler-Maruyamova metóda)

Namiesto infinitezimálneho časového kroku $dt$ použijeme konečný časový krok $\Delta t$. Namiesto $dX$ budeme uvažovať $X_{t + \Delta t} - X_t$ a namiesto infinitezimálneho prírastku Wienerovho procesu použijeme náhodnú premennú so strednou hodnotou $0$ a disperziou $\Delta t$.

Namiesto rovnice v tvare
$$
  dX_t = \mu(X_t, t)\, dt + \sigma(X_t, t) \, dw_t
$$
budeme počítať
$$
  X_{t+\Delta t} - X_t = \mu(X_t, t)\cdot \Delta t + \sigma(X_t, t)\cdot \mathcal N(0, \Delta t),
$$
resp.
$$
  X_{t+\Delta t} = X_t + \mu(X_t, t)\cdot \Delta t + \sigma(X_t, t)\cdot \mathcal N(0, \Delta t),
$$

Poznámka: existujú aj iné numerické schémy na simuláciu náhodných procesov, napríklad Milsteinova metóda alebo metóda Leimkuhlera a Matthewsa.

```{r}
Euler_Maruyama <- function(Drift, Volatility, t_points, x0 = 0) {
  n <- length(t_points)
  dt <- diff(t_points)
  dw <- rnorm(n - 1, mean = 0, sd = sqrt(dt))
  
  x <- numeric(n)
  x[1] <- x0
  for (i in 1:(n-1)) {
    x[i+1] <- x[i] + 
              Drift(x[i], t_points[i]) * dt[i] + 
              Volatility(x[i], t_points[i]) * dw[i]
  }
  
  return(x)
}
```


Riešme pre $dX = X\,dt + \frac15\,dw_t$ s počiatočnou podmienkou $X_0 = 10$ na intervale $[0, 4]$ a s časovým krokom $\Delta t = 1/365$.
```{r}
x0 <- 10
t_max <- 4
delta_t <- 1/365

Drift <- function(x, t) {
  return(1 * x)
}

Volatility <- function(x, t) {
  return(0.2 * x)
}


t_points <- seq(0, t_max, by = delta_t)

x <- Euler_Maruyama(Drift, Volatility, t_points = t_points, x0 = x0)
plot(t_points, x, "l", xlab = "t", ylab = "X(t)")
```


**Bonus:** Implmentujte Milsteinovu metódu.

**Bonus:** Nájdite príklad takého stochastického procesu, že Euler-Maruyamova metóda môže dať kvalitatívne zlý priebeh (napríklad trajektória, ktorá má byť nezáporná bude obsahovať hodnotu menšiu ako 0).



### Príklad

Simulujte proces $X_t$ daný SDR $dX_t = (3X_t - X_t^3)dt + dw_t$. Zvoľte časový interval $[0, 1000]$, počiatočnú podmienku $X_0 = 1$ a časový krok $\Delta t \leq 0,\!01$. Vykreslite viacero simulácii a slovne popíšte čo vidíte.
```{r eval = FALSE}
# Po dokončení vymažte `eval = FALSE`

t_max <- 1000
delta_t <- 0.01
t_points <- seq(0, t_max, by = delta_t)

Drift <- function(x, t) {
  # TODO
}

Volatility <- function(x, t) {
  # TODO
}


for (i in 1:5) {
  x <- Euler_Maruyama(Drift, Volatility, t_points = t_points, x0=1)
  plot(t_points, x, "l", xlab = "t", ylab = "X(t)")
}

```



### Príklad

Simulujte pribeh stochastického procesu $Y_t = X_t^2$, kde $dX = (3X - X^3)\,dt + dw_t$. Zvoľte počiatočnú podmienku $X_0 = 1$ a časový krok $\Delta t \leq 0,\!01$. Urobte to dvoma spôsobmi:

- Najprv použite Euler-Maruyamovu metódu na proces $Y_t$.
- Potom vygenerujte priebeh procesu $X_t$, použijúc tie isté hodnoty $dw_t$ ako v predošlom bode. Takto získanú trajektóriu transformujte funkciou $g(t, x) = x^2$.
Porovnajte získané trajektórie.

```{r}
Euler_Maruyama <- function(Drift, Volatility, t_points, x0 = 0, dw = NULL){
  n <- length(t_points)
  dt <- diff(t_points)
  
  if (is.null(dw)) {
    dw <- rnorm(n - 1, mean = 0, sd = sqrt(dt))
  }
  
  x <- numeric(n) # pre-alokácia
  x[1] <- x0
  for (i in 1:(n-1)) {
    x[i+1] <- x[i] + 
              Drift(x[i], t_points[i]) * dt[i] + 
              Volatility(x[i], t_points[i]) * dw[i]
  }
  
  return(x)
}
```

```{r}
# TODO
```



## Príklad

Predpokladajme, že cena akcie sa riadi geometrickým Brownovým pohybom:
$$
  dS = 0,\!3S\, dt + 0,\!25 S \, dw_t.
$$
Dnešná cena akcie je $150$.

- Nakreliste histogram hodnôt ceny akcie o mesiac, ziskané simuláciou.
- Do toho istého obrázka nakreslite hustotu rozdelenia ceny o mesiac.
- Aká je pravdepodobnosť, že o mesiac bude cena vyššia ako 140?
- Aká je stredná hodnota štvrťročného výnosu? Aká je pravdepodobnosť, že bude záporný?



## Ďalšie úlohy

Viac úloh nájdete v [archíve cvičení](http://www.iam.fmph.uniba.sk/institute/stehlikova/fd17/cv02.html).





# Európske opcie

-   **Európska call opcia** je kontrakt v ktorom majiteľ, opcie získava právo (ale nie povinnosť) kúpiť akciu v presne určenom exspiračnom čase za vopred dohodnutú exspiračnú cenu $E$.
-   **Európska put opcia** je kontrakt, v ktorom majiteľ, opcie získava právo (ale nie povinnosť) predať akciu v presne určenom exspiračnom čase za vopred dohodnutú exspiračnú cenu $E$.
-   **Payoff** je hodnota opcie v čase exspirácie
    -   Call opcia má payoff $\max(S-E,0)$
    -   Put opcia má payoff $\max(E-S,0)$
-   Profit je zisk, t.j. payoff znížený o (odúročenú) cenu opcie



## Reálne ceny opcií

Ceny opcií sa dajú pozrieť online na [yahoo finance](https://finance.yahoo.com), alebo na [google finance](https://www.google.com/finance/). Symboly opcie majú svoj [štandardný tvar](https://help.yahoo.com/kb/SLN13884.html), v ktorom sú uvedené všetky potrebné detaily opcie. Tento tvar je možné použiť aj pri automatizovanom načítavaní dát.
