Numerické riešenie Black-Scholesovej rovnice

Budeme riešiť Black-Scholesovu parciálnu diferenciálnu rovnicu numericky pomocou explicitnej metódy. Máme k tomu dva dôvody: - nie pre každý payoff sa dá Black-Sholesova rovnica riešiť explicite, - potrebujeme začať budovať aparát na numerické oceňovanie amerických opcií.

Budeme sa zaoberať prípadom, v ktorom podkladové aktívum nevypláca dividendy. 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

Namiesto rovnice pre \(V\) budeme riešiť transformovanú rovnicu pre funkciu \(Z(x, \tau)\), pričom \(Z(x, \tau) = V(S, t)\), kde \(x = \ln(S/E)\) a \(\tau = T - t\). V pôvodnej rovnici \(S \in (0, \infty)\). Po transformácii bude \(x \in (-\infty, \infty)\). Premennú \(x\) budeme nazývať priestorová premenná a premennú \(\tau\) zas časová premenná.

Dokončite túto transformáciu (štátnicová otázka): \[ \begin{align*} \frac{\partial V}{\partial t} &= \frac{\partial Z}{\partial \tau} \frac{\partial \tau}{\partial t} \\ &=- \frac{\partial Z}{\partial \tau} \\ & \\ \frac{\partial V}{\partial S} &= \frac{\partial Z}{\partial x} \frac{\partial x}{\partial S} \\ &= \frac{\partial Z}{\partial x} \frac{1}{S} \\ & \\ \frac{\partial^2 V}{\partial S^2} &= \frac{\partial}{\partial S} \left(\frac{\partial V}{\partial S} \right) \\ &= \frac{\partial}{\partial S} \left( \frac{\partial Z}{\partial x} \frac{1}{S} \right) \\ &= \dots \end{align*} \]

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)\). Návod: \(u(x, \tau) = \exp(\alpha x + \beta \tau) V(S, t)\), kde \(S = E e^x\) a \(t = T - \tau\).

Diskretizácia premenných a numerická schéma

Namiesto hodnôt funkcie \(u\) pre všetky (spojité) hodnoty \(x\) budeme funkciu \(u\) počítať iba v bodoch \(\ldots, -2, -h, 0, h, 2h, \ldots\). Podobne budeme uvažovať iba časové body \(\tau = 0, k, 2k, 3k, \ldots\).

Namiesto parciálnej derivácie podľa času použijeme doprednú diferenciu: \[ \frac{\partial u}{\partial \tau} \approx \frac{u(x, \tau + k) - u(x, \tau)}{k} \] Všimnime si, že v limite pre \(k \to 0^+\) je pravá strana definíciou parciálnej derivácie.

Namiesto druhej parcálnej derivácie podľa priestorovej premennej použijeme diferenciu \[ \frac{\partial^2 u}{\partial x^2} \approx \frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2} \] Na skúške z matematickej analýzy (1) sa vyskytol príklad, v ktorom bolo treba ukázať, že v limite \(h\to 0^+\) pravá strana konverguje ku druhej derivácii.

Dokopy teda máme približnú rovnosť: \[ \frac{u(x, \tau + k) - u(x, \tau)}{k} - \frac{\sigma^2}{2}\frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2} \approx 0, \] resp. \[ u(x, \tau + k) \approx u(x, \tau) + \frac{k \sigma^2}{2}\frac{u(x - h, \tau) - 2 u(x, \tau) + u(x + h, \tau)}{h^2}. \]

Veľkosti krokov musia spĺňať tzv. Courant–Friedrichs–Lewyho podmienku: \[ k < \frac{h^2}{\sigma^2}. \] Ak by táto podmienka nebola splnená, metóda nebude stabilná. Viac o tejto podmienke sa dozviete na predmete Numerické modelovanie.

Ohraničená oblasť a okrajové podmienky

Numericky však nevieme riešiť pre \(x \in (-\infty, \infty)\). Namiesto toho budeme riešiť na dostatočne dlhom symetrickom intervale \(x \in [-L, L]\), t.j. v bodoch \(x= -L, -L+h, \ldots, -h, 0, h, 2h, \ldots, L-h, L\).

V reči pôvodných premenných \(S\) a \(E\) to znamená, že rovnicu riešime na oblasti: \[ \ln(S/E) \in [-L, L] \\ S \in [E \cdot e^{-L}; E \cdot e^L] \] Ak napríklad \(L = 2\), tak \[ S \in [0,\!135 E;\ 7,\!389E] \]

Pre hodnoty \(x = \pm L\) musíme ešte doplniť okrajové podmienky. Nazveme ich \(u(-L, \tau) = \phi(\tau)\) pre \(x=-L\) a \(u(L, \tau) = \psi(\tau)\) pre \(x=L\). Odvodíme ich tvar:

Ak \(E \to 0\), tak \(Call \sim 0\)
Ak \(E \to \infty\), tak \(Put \to 0\). Potom z put-call parity vyplýva \(Call \sim Put + S - E e^{-r \tau}\), resp. \(Call \sim S - E e^{-r \tau}\).

Bonus: Lepšia asymptotická aproximácia (3b)

Pomocou asymptotickej aproximácie \(N(x)\) pre \(x \to -\infty\) ukážte že lepší odhad ceny call opcie pre malé hodnoty \(S\) je: \[ Call(S) \overset{S \to 0^+}{\sim} \frac{2 \sqrt{2} S \sigma^{3} \tau^{\frac{3}{2}} \exp\left({- \frac{r^{2} \tau}{2 \sigma^{2}} - \frac{r \tau}{2} - \frac{r \ln{\left(\frac{S}{E} \right)}}{\sigma^{2}} - \frac{\sigma^{2} \tau}{8} - \frac{\ln{\left(\frac{S}{E} \right)}^{2}}{2 \sigma^{2} \tau}}\right)}{\sqrt{\pi} \sqrt{\frac{S}{E}} \left(4 r^{2} \tau^{2} + 8 r \tau \ln{\left(\frac{S}{E} \right)} - \sigma^{4} \tau^{2} + 4 \ln{\left(\frac{S}{E} \right)}^{2}\right)}, \]

Implementácia

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

Pri výpočte sa oplatí predrátať konštantu \(\frac{k \sigma^2}{2h^2}\).

E <- 50
r <- 0.04
sigma <- 0.4

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

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

# 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(-L, L, by=h)
tt <- seq(0, tau_max, by=k)

# Okrajové podmienky
phi <- function(tau) 0# TODO 
psi <- function(tau) 0# TODO

# Počiatočná podmienka
u0 <- function(x) 0# TODO

# Matica riešení
u <- matrix(0, nrow=2*n + 1, ncol=m+1)

# Okrajové a počiatočné podmienky
u[1, ] <- phi(tt)
u[2*n +1, ] <- psi(tt)
u[, 1] <- u0(xx)

# sol:
# x = -L   # # # # # # # #
#
# x = 0    # # # # # # # #
# 
# X = +L   # # # # # # # #
#          ^             ^
#        tau=0        tau=tau_max

k <= h^2/(sigma^2) # Podmienka stability

## [1] TRUE

gamma <- k*sigma^2/(2*h^2) # Predrátaná konštanta 

for (j in 2:(m+1)) {
  for (i in 2:(2*n)) {
    u[i, j] <- 0 # 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))
}

# TODO

Chybový graf

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

#TODO

Bonus: Iná logaritmická transformácia (2b)

Namiesto transformácie \(x := \ln(S/E)\) sme mohli použiť transformáciu \(x := \ln(S)\). Podrobne popíšte čo všetko by sa zmenilo v implementácii.

Bonus: 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. Možno bude potrebné použiť interpoláciu.

Zvoľte vhodné parametre delenia \(m\), \(n\) tak aby metóda bola stabilná a dostatočne presná.

Pomôžte si prednáškou

Bonus: Iný 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.

Bonus: Metóda binárneho stromu (3b)

Odvoďte iteračnú schému, ak sa použije krok \(h = \sigma \sqrt k\). Vysvetlite ako táto voľba kroku zjednoduší výpočet.

Implementujte funkciu CallValueNumeric(S, E, r, sigma, tau) pre výpočet ceny call opcie, ak je dnešná cena podkladového aktíva S. Stačí ak výpočet zbehne pre skalárnu hodnotu S (teda nie vektor).

Bonus: Experimenty so stabilitou (2b)

Skúste zvoliť také veľkosti krokov, že:

\(k < \frac{h^2}{2\sigma^2}\)
\(k = \frac{h^2}{2\sigma^2}\)
\(\frac{h^2}{2\sigma^2} < k < \frac{h^2}{\sigma^2}\)
\(k = \frac{h^2}{\sigma^2}\)
\(k > \frac{h^2}{\sigma^2}\)

Vykreslite grafy cien opcií, ktoré vzniknú takýmito metódami.

Cvičenie 5

Finančné deriváty