Cvičenie 6
Finančné deriváty
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?
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\).
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.
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\).
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")## [1] "Lepšia metóda"print(round(toc - tic, 5))## Time difference of 0 secstic <- Sys.time()
solution_2 <- Jacobi(A, b, maxiter = 100)
toc <- Sys.time()
print("Pôvodná metóda")## [1] "Pôvodná metóda"print(round(toc - tic, 5))## Time difference of 0.19359 secsÚ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\).
A <- matrix(c(10, 5, 2, 5, 7, 1, 2, 1, 3), 3, 3, byrow=T)
b <- c(6, 12, 9)
# TODOGauss-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.
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 ] \]
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\).
A <- matrix(c(10, 5, 2, 5, 7, 1, 2, 1, 3), 3, 3, byrow=T)
b <- c(6, 12, 9)
maxiter <- 50
# TODOSOR 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ť.
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 ] \]
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\).
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.
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.
d <- -2
m <- 1
b <- rep(1, 50)
# TODOImplicitná 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\%\)
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.
# TODOPorovnanie s analytickým riešením
Porovnajte výsledok s cenou zistenou analyticky.
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}\).
#TODOBonus (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.