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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Na dnešnom cvičení si zopakujeme numerické metódy na riešenie sústavy lineárnych rovníc: Jakobiho metódu, Gauss-Seidelovu metódu a SOR metódu. Potom odvodíme implicitnú schému pre oceňovanie európskych opcií a implementujeme ju.

Prednášky k cvičeniu: http://www.iam.fmph.uniba.sk/institute/stehlikova/fd21/fd_05_sor_metoda.pdf

# Numerické metódy

## Jakobiho metóda

Ak riešime sústavu lineárnych rovníc
$$
  Ax = b, 
$$
potom Jakobiho iteračná metóda hovorí toto:

Maticu $A$ si rozložíme na diagonálne a mimodiagonálne prvky:
$$
  A = D+M
$$
a potom platí:
$$
  Ax = b \\
  (D+M)x = b\\
  Dx + M x = b\\
  Dx = b - M x\\
  x = D^{-1}\left( b-M x\right)
$$
Pre niekoré matice (napríklad striktne diagonálne dominantné) sa dá ukázať, že toto zobrazenie je kontraktívne, preto podľa Banachovej vety predpis
$$
  x^{k+1} = D^{-1}\left( b-M x^k\right)
$$
konverguje ku presnému riešeniu sústavy $Ax = b$.


Pozrite si nasledujúcu funkciu. Čo sú jej vstupy a čo je jej výstupom?
```{r}
Jacobi <- function(A, b, x0 = NA, maxiter = 100) {
  d <- diag(A)
  d_ <- 1/d   # D^{-1}
  M <- A - diag(d)
  
  n <- length(d) # počet neznámych
  
  # Pokiaľ nemáme počiatočný odhad x0, 
  # použijeme vektor samých núl.
  if (any(is.na(x0))) {
    x0 <- rep(0, n)
  }
  
  x_k <- matrix(0, maxiter+1, n)
  x_k[1, ] <- x0
  
  for (k in 1:maxiter) {
    x_k[k+1, ] <- d_*(b - M%*%x_k[k, ])
  }
  
  return(x_k)
}
```

Vyskúšame Jakobiho metódu na riešenie sústavy
$$
  \left[\begin{matrix}
   8 & 5\\ 
    5 & 7
  \end{matrix}\right ] x = \left[\begin{matrix}
    13\\ 
    12
  \end{matrix} \right ]
$$

Ako počiatočný odhad riešení môžeme zvoliť napríklad $x^0 = [0, 0]^T$. Vykreslíme trajektóriu, t.j. postupnosť $x^0, x^1, x^2, \ldots$.
```{r}
A <- matrix(c(8, 5, 5, 7), 2, 2, byrow=T)
b <- c(13, 12)
x0 <- c(0, 0) # počiatočný odhad riešenia
maxiter <- 20 # počet iterácií

iterations <- Jacobi(A, b, x0 = x0, maxiter = maxiter)

par(pty="s")  # Čo robí tento riadok?
plot(iterations[, 1], iterations[, 2], type = "o", pch = 20, main = "Jakobiho metóda")
```

Skúsme ešte raz s iným počiatočným odhadom.

```{r}
iterations <- Jacobi(A, b, x0=c(-10, 10))

par(pty="s")
plot(iterations[, 1], iterations[, 2], type="o", pch = 20, main="Jakobiho metóda")
```

Poznámka: metóda by sa dala urýchliť, ak by sme skončili predčasne v prípade že $\|x^{k+1} - x^k\| < \varepsilon$. 



### Bonus: Lepšia implementácia (2b)

Matica `M` v Jakobiho metóde je takmer taká istá ako pôvodná matica $A$. Implementujte Jakobiho metódu tak, aby sa nemusela pri výpočte alokovať žiadna matica. Porovnajte čas potrebný pre 100 iterácií týchto dvoch implementácií pre maticu $A$ rozmeru $1\ 000 \times 1\ 000$.

```{r}
JacobiBetter <- function(A, b, x0 = NA, maxiter=100) {
  # TODO
}

n <- 1000
A <- (matrix(rnorm(n*n), nrow = n))
A <- A %*% t(A)
b <- A %*% rep(1, n)


tic <- Sys.time()
solution_1 <- JacobiBetter(A, b, maxiter = 100)
toc <- Sys.time()

print("Lepšia metóda")
print(round(toc - tic, 5))


tic <- Sys.time()
solution_2 <- Jacobi(A, b, maxiter = 100)
toc <- Sys.time()

print("Pôvodná metóda")
print(round(toc - tic, 5))
```




### Úloha: Jakobiho metóda pre 3 premenné

Použite Jakobiho schému pre sústavu:
$$
  \left[\begin{matrix}
  10 & 5 & 2\\ 
  5  & 7 & 1\\
  2  & 1 & 3
  \end{matrix}\right ] x
    =
  \left[\begin{matrix}
  6\\12\\9
  \end{matrix} \right ]
$$
Použite 50 iterácií. Porovnajte súčin $Ax_{50}$ s pravou stranou $b$.

```{r}
A <- matrix(c(10, 5, 2, 5, 7, 1, 2, 1, 3), 3, 3, byrow=T)
b <- c(6, 12, 9)
# TODO
```



## Gauss-Seidelova metóda

Namiesto násobenia $D^{-1}(b-Mx_k)$ naraz budeme násobiť postupne: 
Najprv určíme prvú zložku $x^{k+1}_1 = D^{-1}_1 (b_1 - m_1^Tx^k)$, rovnako ako by sme urobili pri Jakobiho schéme. Túto hodnotu ale použijeme pri rátaní druhej zložky. 

Gauss-Seidelova metóda konverguje lepšie ako Jakobiho metóda: dá sa použiť pre širšiu triedu úloh (napríklad stačí kladná definitnosť matice) a konvergencia je rýchlejšia.

```{r}
GaussSeidel <- function(A, b, x0 = NA, maxiter = 100) {
  d <- diag(A)
  d_ <- 1/d
  M <- A - diag(d)
  
  n <- length(d)
  
  if (any(is.na(x0))) {
    x0 <- rep(0, n)
  }
  
  x_k <- matrix(0, maxiter+1, n)
  x_k[1, ] <- x0
  
  x_k <- x0
  
  for (k in 1:maxiter) {  # iterácie
    for (j in 1:n) {      # premenné
      x_k[j] <- (b[j] - sum(M[j, ] * x_k)) * d_[j]
    }
    iterations[k+1, ] <- x_k
  }
  
  return(iterations)
}
```



Opäť vyskúšame metódu na sústave
$$
  \left[\begin{matrix}
  8 & 5\\ 
  5 & 7
  \end{matrix}\right ] x = \left[\begin{matrix}
  13\\ 
  12
  \end{matrix} \right ]
$$
```{r}
A <- matrix(c(8, 5, 5, 7), 2, 2, byrow=T)
b <- c(13, 12)
iterations <- GaussSeidel(A, b)

par(pty="s")
plot(iterations[, 1], iterations[, 2], type = "o", pch = 20, main = "Gauss-Seidelova metóda")
```


### Úloha: Gauss-Seidelova metóda pre 3 premenné

Použite Gauss-Seidelovu schému pre sústavu:
$$
  \left[\begin{matrix}
  10 & 5 & 2\\ 
  5  & 7 & 1\\
  2  & 1 & 3
  \end{matrix}\right ] x
    =
  \left[\begin{matrix}
  6\\12\\9
  \end{matrix} \right ]
$$
Použite 50 iterácií. Porovnajte súčin $Ax_{50}$ s pravou stranou $b$.
```{r}
A <- matrix(c(10, 5, 2, 5, 7, 1, 2, 1, 3), 3, 3, byrow=T)
b <- c(6, 12, 9)
maxiter <- 50

# TODO
```



## SOR metóda

Gauss-Seidelovu metódu možno ešte vylepšiť. Zachováme smer, ale v danom smere zvolíme dlhší krok $\omega$. Pre $\omega = 1$ je metóda ekvivalenentá s Gauss-Seidelovou metódou. Pre $\omega > 2$ bude metóda divergovať.


```{r}
Sor <- function(A, b, x0 = NA, maxiter = 100, w = 1.2){
  d <- diag(A)
  d_ <- 1/d
  M <- A - diag(d)
  
  n <- length(d)
  
  if (any(is.na(x0))){
    x0 <- rep(0, n)
  }
  
  iterations <- matrix(0, maxiter+1, n)
  iterations[1, ] <- x0
  
  x_k <- x0
  
  for (k in 1:maxiter){ # iterácie
    for (j in 1:n){ # premenné
      step <- (b[j] - sum(M[j, ] * x_k)) * d_[j] - x_k[j]
      x_k[j] <- x_k[j] + w*step 
    }
    iterations[k+1, ] <- x_k
  }
  
  return(iterations)
}
```



Vyskúšame metódu na sústave
$$
  \left[\begin{matrix}
  8 & 5\\ 
  5 & 7
  \end{matrix}\right ] x = \left[\begin{matrix}
  13\\ 
  12
  \end{matrix} \right ]
$$

```{r}
A <- matrix(c(8, 5, 5, 7), 2, 2, byrow=T)
b <- c(13, 12)
iterations <- Sor(A, b)

par(pty="s")
plot(iterations[, 1], iterations[, 2], type = "o", pch = 20, main = "SOR metóda")
```


### Úloha: SOR metóda pre 3 premenné

Vyskúšajte SOR metódu na sústave:
$$
  \left[\begin{matrix}
  10 & 5 & 2\\ 
  5  & 7 & 1\\
  2  & 1 & 3
  \end{matrix}\right ] x
    =
  \left[\begin{matrix}
  6\\12\\9
  \end{matrix} \right ]
$$
Použite 50 iterácií. Porovnajte súčin $Ax_{50}$ s pravou stranou $b$.

```{r}
A <- matrix(c(10, 5, 2, 5, 7, 1, 2, 1, 3), 3, 3, byrow=T)
b <- c(6, 12, 9)

# TODO
```

### Úloha: SOR metóda pre symetrickú trojdiagonálnu maticu

Naprogramujte SOR metódu, ak matica $A$ je symetrická trojdiagonálna, t.j.:
$$
  A =   \left[\begin{matrix}
  d & m \\
  m & d & m\\
  & m & d & m\\
  && \ddots & \ddots & \ddots\\
  &&& m & d & m\\
  &&& & m & d\\
  \end{matrix}\right ]
$$
Využite štruktúru matice tak, aby funkcia zbehla čo najrýchlejšie, t.j. rozmyslite si ako spočítať `sum(M[j, ] * x_k)` iba s použitím jedného násobenia.

```{r}
SorTridiagonal <- function(d, m, b, x0=NA, maxiter=100) {
  n <- length(b)
  d_ = 1/d
  
  if(any(is.na(x0))) {
    x0 <- rep(0, n)
  }
  
  x_k <- x0
  
  for (k in 1:maxiter) {
    x_k[1] <- 0 # FIXME
    for (j in 2:(n-1)) {
      x_k[j] <- 0 #FIXME
    }
    x_k[n] <- 0 # FIXME
  }
  
  return(x_k)
}
```

### Úloha: SOR metóda a matica druhých diferencií

Otestujte túto metódu pre trojdiagonálnu maticu s hodnotami $d = -2, m = 1$. Ako vektor pravej strany zvoľte vektor samých jednotiek. Riešenie vykreslite, malo by mať tvar paraboly.
```{r}
d <- -2
m <- 1
b <- rep(1, 50)

# TODO
```






# Implicitná schéma pre európsku opciu


## Odvodenie schémy 

Prednáška k cvičeniu: http://www.iam.fmph.uniba.sk/institute/stehlikova/fd17/06_numerika_euro.pdf

Na minulom cvičení sme použili doprednú pomernú diferenciu na aproximáciu parciálnej derivácie podľa času:
$$
  \frac{\partial u (x, \tau)}{\partial \tau} \approx \frac{u(x, \tau+k) - u(x, \tau)}{k}
$$

Čo sa ale stane, ak použijeme spätnú diferenciu?
$$
 \frac{\partial u (x, \tau)}{\partial \tau} \approx \frac{u(x, \tau) - u(x, \tau-k)}{k}
$$

Dosadením do rovnice vedenia tepla dostávame:
$$
  \frac{u(x, \tau) - u(x, \tau-k)}{k} - \frac{\sigma^2}{2} \frac{u(x-h, \tau) - 2u(x, \tau) + u(x+h, \tau)}{h^2} \approx 0, 
$$
čím získame numerickú schému
$$
  u(x, \tau) - k \frac{\sigma^2}{2} \frac{u(x-h, \tau) - 2u(x, \tau) + u(x+h, \tau)}{h^2} = u(x, \tau-k), 
$$
resp.
$$
  -\gamma u(x-h, \tau) + (1 + 2\gamma)u(x, \tau) - \gamma u(x+h, \tau) = u(x, \tau-k), 
$$
kde $\gamma = \frac{\sigma^2 k}{2h^2}$. Toto je sústava $2n+1$ lineárnych rovníc, ktorú musíme riešiť. Dá sa zapísať v maticovom tvare
$$
  \left[\begin{matrix}
  1\\
  -\gamma & 1 + 2\gamma & -\gamma\\
  &-\gamma & 1 + 2\gamma & -\gamma\\
  && \ddots & \ddots & \ddots\\
  &&&-\gamma & 1 + 2\gamma & -\gamma\\
  &&&&&1
  \end{matrix}\right ] \left[\begin{matrix}
  | \\\\u(\tau) \\ \\ |
  \end{matrix} \right ]
    =
  \left[\begin{matrix}
  \phi \\ \hline | \\ u(\tau - k) \\ | \\ \hline \psi
  \end{matrix} \right ]
$$

Vďaka okrajovým podmienkam $\phi, \psi$ poznáme dve z týchto hodnôt. Ak ich zahrnieme do pravej strany, dostaneme symetrickú, trojdiagonálnu, striktne diagonálne dominantnú maticu:
$$
  \left[\begin{matrix}
  1 + 2\gamma & -\gamma\\
  -\gamma & 1 + 2\gamma & -\gamma\\
  & \ddots & \ddots & \ddots\\
  &&-\gamma & 1 + 2\gamma & -\gamma\\
  &&&-\gamma & 1 + 2\gamma\\
  \end{matrix}\right ] \left[\begin{matrix}
  |\\ \\u(\tau) \\ \\ |
  \end{matrix} \right ]
    =
  \left[\begin{matrix}
  | \\ \\ u(\tau - k) \\ \\|
  \end{matrix} \right ] + \left[\begin{matrix}
  \gamma \phi \\ 0 \\ \vdots \\ 0 \\ \gamma \psi
  \end{matrix} \right ]
$$

Hoci táto schéma je výpočtovo oveľa náročnejšia, je stabilná nezávisle od veľkosti krokov $h, k$. Stále však platí že čím menší krok, tým presnejší výsledok.



## Ocenenie európskej *call* opcie

Pomocou implicitnej schémy 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
q <- 0
sigma <- 0.4
tau_max <- 1

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

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

# počet iterácií Gauss-Seidelovej metódy
maxiter <- 100

# Pomocné konštanty
alpha <- (r-q)/(sigma^2) - 0.5
beta  <- (r+q)/2 + sigma^2 / 8 + (r-q)^2/(2*sigma^2)
gamma <- sigma^2 * k / (2 * h^2)

a <- 1 + 2*gamma  # členy na diagonále
b <- -gamma  # členy na vedľajších diagonálach

# Diskretizácia času a priestoru
xx <- seq(-L, L, by=h)
tt <- seq(0, tau_max, by=k)

# Okrajové podmienky
phi <- function(tau) 0 
psi <- function(tau) E*exp(alpha*L + beta*tau)*(exp(L - q*tau) - exp(-r*tau))

# Počiatočná podmienka
u0 <- function(x) E*exp(alpha*x)*pmax(0, exp(x)-1)

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


for (i in 2:(m+1)) { # Ideme naplniť i-ty stĺpec
  b <- 0 # FIXME
  SorTridiagonal(a, b, b, x0 = NA, maxiter = 1) # FIXME
}

```

## Spätná transformácia

Získané riešenie transformujte naspäť do premenných $S, t$. Vykreslite získané riešenie.
```{r}
# TODO
```



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

Porovnajte výsledok s cenou zistenou analyticky.

```{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))
}



```


## Chybový graf

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

```{r}
#TODO
```




# Bonus (max. 8b)

Z nasledujúcich úloh si vyberte dve a vyriešte ich. Ak vyriešite viac úloh, budú sa počítať iba dve s najvyšším bodovým ziskom.

## Crank-Nicholsonova metóda (5b)

Ako aproximáciu druhej derivácie možno použiť aj:
$$
  \frac{\partial^2 u}{\partial x^2} \approx \frac12 \left[\frac{u(x-h, \tau) - 2u(x, \tau) + u(x+h, \tau)}{h^2} + \frac{u(x-h, \tau-k) - 2u(x, \tau-k) + u(x+h, \tau-k)}{h^2} \right]
$$

Odvoďte implicitnú numerickú schému rovnice vedenia tepla pre takúto aproximáciu druhej derivácie a zapíšte ju v maticovom tvare.


Použite túto metódu na numerické ocenenie európskej put opcie ktorá 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` a vhodný počet iterácií `maxiter` tak aby bola metóda dostatočne presná.



## Normálny systém (2b)

Systém $Ax = b$ možno vyriešiť Gauss-Seidelovou metódou, aj pokiaľ matica $A$ nemá požadované vlastnosti. Stačí ak sústavu prenásobíme zľava maticou $A^T$, dostaneme:
$$
  A^T A x = A^T b
$$
Ak je matica $A$ regulárna, tak $A^T A$ je kladne definitná, čiže Gauss-Seidelova metóda bude konvergovať ku vektoru $A^{-1} b$.

Ku akému riešeniu bude metóda konvergovať, ak by neexistovala matica $A^{-1}$? Rozpíšte jednotlivé prípady v ktorých neexistuje inverzná matica: prípad štvorcovej singulárnej matice, ako aj všetky prípady obdĺžnikových matíc.


## Iná aproximácia (3b)

Ukážte, že druhú deriváciu možno aproximovať aj ako:
$$
  \frac{\partial^2 u(x, \tau)}{\partial x^2} \approx \dfrac{- u(x-2h, \tau) + 16 u(x-h, \tau) - 30 u(x, \tau) + 16 u(x+h, \tau) - u(x+2h, \tau)}{12 h^{2}},
$$
t.j. dokážte že v limite $h\to 0$ táto hodnota konverguje ku $\frac{\partial^2 u(x, \tau)}{\partial x^2}$. Určte rád diskretizačnej chyby, t.j. také $q$, že
$$
  \lim_{h\to0} \frac{1}{h^q}\left[\dfrac{- u(x-2h, \tau) + 16 u(x-h, \tau) - 30 u(x, \tau) + 16 u(x+h, \tau) - u(x+2h, \tau)}{12 h^{2}} - \frac{\partial^2 u(x, \tau)}{\partial x^2}\right].
$$
je konečné nenulové číslo.


Podrobne popíšte ako by sa zmenil výpočet, ak by sme použili túto aproximáciu.



