Organizácia cvičení

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

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

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

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:

2 + 2
## [1] 4

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. S miernymi modifikáciami ho používa aj google.

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:

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\):

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.

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.

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

print(round((toc-tic), 5))
## Time difference of 2.44789 secs
WienerLoopless <- function(n, dt) {
  # TODO
}

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

print(round((toc-tic), 5))
## Time difference of 0.001 secs

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.

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

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.

# 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

# Po dokončeneí vymažte `eval = FALSE`

GeometricBrownianMotion <- function(n, delta_t, mu, sigma, s0) {
  # TODO
}
# 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.

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

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.

# 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.
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)
}
# 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í.

Európske opcie

Reálne ceny opcií

Ceny opcií sa dajú pozrieť online na yahoo finance, alebo na google finance. Symboly opcie majú svoj štandardný tvar, v ktorom sú uvedené všetky potrebné detaily opcie. Tento tvar je možné použiť aj pri automatizovanom načítavaní dát.