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

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)
## [1] TRUE
# 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
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.

# TODO

Porovnanie s analytickým riešením

Porovnajte získané numerické riešenie s riešením získaným pomocou vzorca.

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

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.