Kurs:Modellierung und Numerische Methoden von Finanzderivaten/6 Numerische Lösung parabolischer Differentialgleichungen

In Kapitel 5 haben wir die Monte-Carlo-Methode zur Lösung stochastischer Differentialgleichungen kennengelernt und gesehen, dass diese Methode im allgemeinen recht zeit- und rechenintensiv ist. Die Preise exotischer Optionen können häufig auch durch die Lösung einer partiellen Differentialgleichung vom Black-Scholes-Typ bestimmt werden. Diese Differentialgleichungen können allerdings im allgemeinen nicht explizit gelöst werden. In diesem Kapitel stellen wir einige Techniken vor, mit denen diese Gleichungen numerisch gelöst werden können.

Als einführendes Beispiel leiten wir eine parabolische Differentialgleichung zur Bewertung asiatischer Optionen her (Abschnitt 6.1). In den folgenden Abschnitten erläutern wir zwei numerische Techniken, mit denen die Gleichung für asiatische Optionen gelöst werden kann: die Methode der Finiten Differenzen (Abschnitt 6.2) und eine Anwendung auf Power-Optionen (Abschnitt 6.3) sowie die (vertikale) Linienmethode (Abschnitt 6.4) und eine Anwendung auf Basket-Optionen (Abschnitt 6.5).

6.1 Asiatische Optionen

Asiatische Optionen sind Optionen, deren Auszahlungsfunktion von den Kursen des Basiswerts $S(t)$ abhängt, gemittelt über die Laufzeit $0\le t\le T$ der Option. Je nachdem, wie der Mittelwert gebildet wird, gibt es verschiedene Typen asiatischer Optionen:

• Der Mittelwert $\hat{S}$ bezieht sich auf diskrete Zeitpunkte $S(t_1),\dots ,S(t_n)$ oder auf stetige Zeitpunkte $S(t),0\le t\le T$.
• Der Mittelwert $\hat{S}$ ist arithmetisch oder geometrisch.

Außerdem werden die Optionen je nach den Auszahlungsmodalitäten unterschieden:

• Besitzt die Option die Auszahlungsfunktion $(S-\hat{S}{)}^{+}$ (für einen Call) bzw. $(\hat{S}-S{)}^{+}$ (für einen Put), so nennen wir sie strike option; im Falle eines Payoffs $(\hat{S}-K{)}^{+}$ (Call) bzw. $(K-\hat{S}{)}^{+}$ (Put) mit Ausübungspreis $K$ heißt sie rate option.
• Die asiatische Option kann vom europäischen oder amerikanischen Typ sein.

Es bleibt der Mittelwert $\hat{S}$ für die verschiedenen Situationen zu definieren. Im diskreten Fall haben wir:

arithmetisch: $\hat{S}=\frac{1}{n}\sum _{k=1}^{n}S(t_k),$

geometrisch: $\hat{S}={(\prod _{k=1}^{n}S(t_k))}^{1/n}.$

Den stetigen Fall erhalten wir im (formalen) Grenzwert $n\to \mathrm{\infty }$, wenn wir $T=n\mathrm{\Delta }t$ für festes $T$ schreiben. Bei der arithmetischen Mittelung erhalten wir ein Integral:

$\hat{S}=\frac{1}{T}\sum _{k=1}^{n}S(t_k)\mathrm{\Delta }t\to \frac{1}{T}\underset{0}{\overset{T}{\int }}S(\tau )d\tau (\mathrm{\Delta }t\to 0).$

Der Grenzwert bei der geometrischen Mittelung ist etwas komplizierter:

$\hat{S}=\exp (\frac{1}{n}\sum _{k=1}^n\ln S(t_k))=\exp (\frac{1}{T}\sum _{k=1}^n\ln S(t_k)\mathrm{\Delta }t)\to \exp (\frac{1}{T}\underset{0}{\overset{T}{\int }}\ln S(\tau )d\tau )(\mathrm{\Delta }t\to 0).$

Die Auszahlungsfunktion eines arithmetic-average-strike calls beispielsweise lautet:

$\mathrm{\Lambda }={(S-\frac{1}{T}\underset{0}{\overset{T}{\int }}S(\tau )d\tau )}^{+}.$

Wir betrachten im folgenden nur asiatische Optionen vom europäischen Typ, deren Auszahlungsfunktion nur von $S$ und

(6.1) $I_t=\underset{0}{\overset{t}{\int }}f(S_{\tau },\tau )d\tau$

abhängt, d.h. $\mathrm{\Lambda }=\mathrm{\Lambda }(S,I)$. Im Falle eines arithmetic-average-strike calls ist $f(S,t)=S$ und $\mathrm{\Lambda }=(S-I/T{)}^{+}$. Für derartige Optionen wollen wir den fairen Preis als Lösung einer parabolischen Differentialgleichung bestimmen, indem wir die Duplikationsstrategie aus Abschnitt 4.2 anwenden. Wir machen dieselben Voraussetzungen an den Finanzmarkt wie in Abschnitt 4.2 bei der Herleitung der Black-Scholes-Gleichung. Insbesondere nehmen wir an, dass die Kurswerte $S_{\tau }$ für $0<\tau <t$ unabhängig von $S_t$ sind. Dies bedeutet, dass der Kurs des Basiswerts nicht von den vergangenen Kursen abhängen soll. Daher können wir $I$ als eine von $S$ unabhängige Variable ansehen. Der Optionspreis wird also im allgemeinen von $S,I$ und $t$ abhängen.

Seien wie in Abschnitt 4.2 $B_t$ ein Bond mit risikofreier Zinsrate $r$ und $V(S,I,t)$ der Optionspreis. Ferner sei $Y_t=c_1(t)B_t+c_2(t)S_t-V(S,I,t)$ ein risikoloses Portfolio mit infinitesimaler Änderung

(6.2) $d{Y}_t=c_1(t)d{B}_t+c_2(t)d{S}_t-dV(S,I,t).$

Der Basiswert und der Bond ändern sich nach Voraussetzung gemäß

$d{S}_t=\mu S_{t}dt+\sigma S_{t}d{W}_t,d{B}_t=r{B}_{t}dt$

und die stochastische Differentialgleichung für $I_t$ lautet wegen Gleichung (6.1):

$d{I}_t=f(S_t,t)dt.$

Nach einer Variante des Lemmas 3.1 von Itô wird die Änderung von $V$ beschrieben durch

$dV=V_{t}dt+V_{S}dS+V_{I}dI+\frac{1}{2}{V}_{SS}d{S}^2+V_{SI}dSdI+\frac{1}{2}{V}_{II}d{I}^2+o(dt),$

wobei von nun an $V_t$ die partielle Ableitung von $V$ nach $t$ bezeichnet und $o(dt)$ Terme von ”kleinerer” Größenordnung als $dt$ sind. (Genauer gilt $g=o(t)$ für $\mathrm{\Delta }t\to 0$ genau dann, wenn $|g/\mathrm{\Delta }t|\to 0$ für $\mathrm{\Delta }t\to 0$.) Setzen wir die stochastischen Differentialgleichungen für $S_t$ und $I_t$ ein, benutzen die Merkregel $d{W}^2=dt$ und vernachlässigen Terme höherer Ordnung, so erhalten wir

$dV=(V_t+\mu S{V}_S+\frac{1}{2}{\sigma }^2{S}^2{V}_{SS}+f(S,t)V_I)dt+\sigma S{V}_{S}dW.$

Diese Formel ist das Analogon zu Gleichung (4.5).

Wir sind nun in der Lage, die stochastischen Differentialgleichungen für $B,S$ und $V$ in (6.2) einzusetzen:

$dY=(c_{1}rB-V_t+\mu S(c_2-V_S)-\frac{1}{2}{\sigma }^2{S}^2{V}_{SS}-f(S,t)V_I)dt+\sigma S(c_2-V_S)dW.$

Setzen wir wie in Abschnitt 4.2 $c_2=V_S$, so wird das Portfolio $Y$ risikofrei. Andererseits gilt aus Arbitrage-Gründen, dass

$dY=rYdt=r(c_{1}B+S{V}_S-V)dt.$

Gleichsetzen der beiden letzten Gleichungen und Identifizieren der Koeffizienten vor dt ergibt eine parabolische Differentialgleichung für den Optionspreis $V(S,I,t)$:

(6.3) $V_t+\frac{1}{2}{\sigma }^2{S}^2{V}_{SS}+rS{V}_S+f(S,t)V_I-rV=0,0<S,I<\mathrm{\infty },0<t<T.$

Die Gleichung ist zu vervollständigen mit End- und Randbedingungen:

$V(S,I,T)=V_0(S,I),(S,I)\in (0,\mathrm{\infty }{)}^2,V(S,I,t)=V_1(S,I,t),$

$S=0,I\in (0,\mathrm{\infty }),t\in (0,T)$ und

$S\to \mathrm{\infty },I\in (0,\mathrm{\infty }),t\in (0,T)$ und

$S\in (0,\mathrm{\infty }),I=0,t\in (0,T)$ und

$S\in (0,\mathrm{\infty }),I\to \mathrm{\infty },t\in (0,T).$

Die Definitionen von $V_0$ und $V_1$ hängen vom speziellen Optionstyp ab. Für arithmetic-average-strike calls gilt etwa $V_0(S,I)=(S-I/T{)}^{+}$.

Im Vergleich zur Black-Scholes-Gleichung (4.8) erhalten wir eine Differentialgleichung mit partiellen Ableitungen in $S,t$ und $I$. Diese Gleichung (mit geeigneten End- und Randbedingungen) ist im allgemeinen nicht explizit lösbar. Die numerische Lösung von (6.3) ist recht aufwendig, da die Ableitung $V_{II}$ fehlt und ein Problem in drei Variablen gelöst werden muss. Für eine spezielle Klasse von Auszahlungsfunktionen und für arithmetische Mittelungen kann Gleichung (6.3) jedoch in eine Differentialgleichung in zwei Variablen transformiert werden.

Satz 6.1

Sei $V$ eine Lösung der Gleichung (6.3) mit $f(S,t)=S$ und sei die Auszahlungsfunktion $\mathrm{\Lambda }$ durch

$\mathrm{\Lambda }(S,I,t)=S^{\alpha }F(\frac{I}{S},t)$

für ein $\alpha \in ℝ$ und eine Funktion $F=F(R,t)$ mit $R=\frac{I}{S}$ gegeben. Dann löst die Funktion $H(R,t)=S^{-\alpha }V$ die Gleichung

(6.4) $H_t+\frac{1}{2}{\sigma }^2{R}^2{H}_{RR}+(1+{\sigma }^2(1-\alpha )R-rR)H_R+(\alpha -1)(\frac{\alpha }{2}{\sigma }^2+r)H=0,H(R,T)=F(R,T),R>0,0<t<T.$

Beweis:

Aus der Definition $V(S,I,T)=S^{\alpha }H(I/S,t)$ folgt durch Differenzieren:

$V_t=S^{\alpha }{H}_t,$

$V_I=S^{\alpha -1}{H}_R,$

$H_S=-S^{-2}I{H}_R=-S^{-1}R{H}_R,$

$V_S=\alpha S^{\alpha -1}H-S^{\alpha -2}I{H}_R=S^{\alpha -1}(\alpha H-R{H}_R),$

$V_{SS}=\alpha (\alpha -1)S^{\alpha -2}H-S^{\alpha -3}I{H}_R-(\alpha -2)S^{\alpha -3}I{H}_R+S^{\alpha -4}{I}^2{H}_{RR}$

$=S^{\alpha -2}(\alpha -1)H-2(\alpha -1)R{H}_R+R^2{H}_{RR}.$

Aus (6.3) folgt mit $f(S,t)=S$:

$0=S^{\alpha }(H_t+\frac{1}{2}{\sigma }^2{R}^2{H}_{RR}+\frac{1}{2}{\sigma }^2\alpha (\alpha -1)H-{\sigma }^2(\alpha -1)R{H}_R+r\alpha H-rR{H}_R+H_R-rH).$

Damit ist der Beweis komplett.

q.e.d.

Beispiel 6.1

Im Falle eines arithmetic-average-strike calls können wir die Gleichung (6.3) mittels Satz 6.1 vereinfachen. Die Voraussetzung des Satzes ist wegen

$\mathrm{\Lambda }=(S-I/T{)}^{+}=S{(1-\frac{I}{ST})}^{+}$

mit $\alpha =1$ und $F(R,t)=(1-R/t{)}^{+}$ erfüllt. Der Preis der Call-Option lautet also $V=S\cdot H$, wobei $H$ die Differentialgleichung

(6.5) $H_t+\frac{1}{2}{\sigma }^2{R}^2{H}_{RR}+(1-rR)H_R=0,R>0,0<t<T$

erfüllt. Die Endbedingung lautet

(6.6) $H(R,T)=(1-R/T{)}^{+},R>0.$

Wir spezifizieren nun geeignete Randbedingungen an $R=0$ und für $R\to \mathrm{\infty }$. Im Falle $R=I/S\to \mathrm{\infty }$ und für festes $I$ muss $S\to 0$ gelten. Die Call-Option wird dann nicht ausgeübt und ist wertlos an $S=0$:

(6.7) $H(R,t)\to 0$ für $R\to \mathrm{\infty },0<t<T.$

Für die Randbedingung an $R=0$ können wir nicht einfach $H(0,t)=0$ wählen (was sich naiv aus $I=0$, $S$ endlich, also $S(\tau )=0$ für $0<\tau <t$ und damit $S(\tau )=0$ für alle $\tau >0$ ergäbe). Um dies heuristisch einzusehen, leiten wir eine stochastische Differentialgleichung für $R$ her. Nach dem Lemma 3.1 von Itô und der Definition des Itô-Prozesses $S$ ist

$dR=(R_t+\mu S{R}_S+\frac{1}{2}{\sigma }^2{S}^2{R}_{SS})dt+\sigma S{R}_{S}dW.$

Aus der Definition $R=I/S$ und $dI/dt=S$ folgt

$R_t=\frac{1}{S}{I}_t=1,$

$R_S=-\frac{I}{S^2}=-\frac{R}{S},$

$R_{SS}=\frac{R}{S^2}-\frac{1}{S}{R}_S=\frac{2R}{S^2}$

und daher

$dR=(1-\mu R+{\sigma }^{2}R)dt-\sigma RdW.$

Interpretieren wir $dR$ als infinitesimale Änderung von $R$, so gilt für $R=0:dR=dt$, d. h., wir haben $R_{t_0+dt}>0$, selbst wenn $R_{t_0}=0$ ist. Es muss also nicht $R_t=0,t>t_0$ bis zum Verfallstag gelten, wenn $R_{t_0}=0$ gilt. Der Wert $H(0,t)$ ist dann aber nicht klar.

Wir leiten eine Randbedingung für $R=0$ her, indem wir annehmen, dass die Differentialgleichung (6.5) auch für $R=0$ gilt. Genauer gesagt muss der Grenzwert $R\searrow 0$ in (6.5) gelten. Wir setzen voraus, dass $H$ zweimal stetig differenzierbar in $R=0,t>0$ ist (gemeint ist der rechtsseitige Grenzwert). Dann ergibt sich für $R\searrow 0$ in (6.5):

(6.8) $H_t+H_R=0,R=0,0<t<T.$

Der Wert V eines arithmetic-average-strike calls ergibt sich also aus $V(S,I,t)=SH(\frac{I}{S},t)$, wobei $H$ die Lösung von (6.5)-(6.8) ist. Üblicherweise diskretisieren wir das Problem (6.5)-(6.8) nicht auf der gesamten positiven Zahlenachse $[0,\mathrm{\infty })$, sondern nur in dem beschränkten Intervall $[0,a]$ mit genügend ”großem” $a>0$; man vergleiche hierzu Abschnitt 6.2.

Welche Möglichkeiten gibt es, die Differentialgleichung (6.5) zu diskretisieren? Wir stellen drei Vorgehensweisen vor.

• Finite Differenzen: Wir versehen das Gebiet $(R,t)\in [0,a]\times [0,T]$ mit einem Gitter

$\{(x_i,t_j):0=x_0<x_1<\dots <x_N=a,0=t_0<t_1<\dots <t_M=T\}.$

Ferner ersetzen wir die partiellen Ableitungen von $H$ durch finite Differenzen, etwa

$H_t(x_i,t_j)\approx \frac{H(x_i,t_{j+1})-H(x_i,t_j)}{t_{j+1}-t_j}$

und definieren Approximationen $u_i^j$ von $H(x_i,t_j)$. Dann erhalten wir ein System von Differenzengleichungen bzw. ein lineares Gleichungssystem in den Unbekannten $u_i^j$ für $0\le i\le N,0\le j\le M$. Im allgemeinen ist dieses Gleichungssystem sehr groß. Wir diskutieren die Methode in Abschnitt 6.2 näher.

• (Vertikale) Linienmethode: Wir approximieren nur die Ableitungen nach der Variablen $R$ durch Differenzenquotienten, d. h. wir überziehen das Gebiet $[0,a]\times [0,T]$ mit vertikalen Linien

$\{(x_i,t):0=x_0<x_1<\dots <x_N=a,0\le t\le T\}$

und approximieren $H(x_i,t)$ durch eine Funktion $u_i(t)$. Wir erhalten ein System gewöhnlicher Differentialgleichungen für $u_i(t)$. Zur Dikussion ist Abschnitt 6.4 vorgesehen.

• (Horizontale) Linienmethode oder Rothe-Methode: Wir überziehen das Gebiet mit horizontalen Linien

$\{(x,t_j):0\le x\le a,0=t_0<t_1<\dots <t_M=T\}$

und approximieren $H(x,t_j)$ durch die Funktion $w^j(x)$, indem wir die partielle Ableitung von $H$ nach $t$ durch einen Differenzenquotient ersetzen. Dies ergibt unter Beachtung der Randbedingungen bei $x=\pm a$ ein lineares Randwertproblem gewöhnlicher Differentialgleichungen 2. Ordnung für die Funktionen $w^j(x)$.

Die beiden zuerst genannten Methoden sind für vorliegende Aufgabe gebräuchlicher, deshalb gehen wir nicht näher auf die letztgenannte Technik ein.

6.2 Methode der finiten Differenzen

Die Idee der Methode der Finiten Differenzen ist es, Differentialquotienten durch Differenzenquotienten zu ersetzen. Wir illustrieren diese Methode nur kurz anhand der einfachen Diffusionsgleichung

(6.9) $u_t-u_{xx}=0,u(x,0)=u_0(x),x\in ℝ,0<t<T,$

die aus der Black-Scholes-Gleichung durch Variablentransformation entsteht (vgl. (4.15)).

Für die numerische Approximation grenzen wir das Problem auf ein beschränktes Intervall $x\in [-a,a]$ mit $a>0$ ein. Wir setzen

(6.10) $u(-a,t)=g_1(t),u(a,t)=g_2(t),0<t<T$

mit geeigneten approximativen Randfunktionen $g_1,g_2$, deren Gestalt von der zu berechnenden Option abhängt. Man kann zeigen, dass die Einschränkung des Problems auf das Intervall $[-a,a]$ die Lösung nicht stark verändert, wenn nur $a>0$ groß genug ist.

Wir geben nun einige Schritte zur Gewinnung der (6.9) entsprechenden Differenzengleichung an. Zunächst führt man Gitterpunkte $(x_i,t_j)$ mit

$x_i=-a+hi(i=0,\dots ,N),t_j=s_j(j=0,\dots ,M)$

und den konstanten Schrittweiten $h=2a/N$ und $s=T/M$ ein. Dann werden die Funktionen $u(x\pm h,t+s)$ nach der Taylorformel bzgl. $u(x,t)$ entwickelt, die Differentialquotienten durch Differenzenquotienten ersetzt (man verwendet sog. Vorwärts-, Rückwärts- und zentrale Differenzenquotienten) und in die Differentialgleichung eingesetzt. Meist nutzt man Kombinationen verschiedener vorwärtiger und rückwärtiger Differenzenquotienten und erhält so rekursiv aufgeschriebene Differenzengleichungen, die ”schichtweise” bzgl. der Zeitschichten $t_j$ gelöst werden können. Wir bezeichnen die approximierende Funktion für $u(x_i,t_j)$ in den Gitterpunkten durch $w_i^j$, man erhält unter Verwendung der Abkürzung $\alpha :=s/h^2$ etwa folgende Darstellung:

(6.11) $-\alpha \theta w_{i+1}^{j+1}+(2\alpha \theta +1)w_i^{j+1}-\alpha \theta w_{i-1}^{j+1}=\alpha (1-\theta )w_{i+1}^j-(2\alpha (1-\theta )-1)w_i^j+\alpha (1-\theta )w_{i-1}^j$

Der Parameter $\theta$ kann hier beliebig aus dem Intervall $[0,1]$ gewählt werden, die Indizes laufen wie folgt: $i=1,\dots ,N-1,j=1,\dots ,M-1$.

Die Rand- und Anfangsbedingungen müssen noch verarbeitet werden:

(6.12) $w_0^j=g_1(t_j),w_N^j=g_2(t_j),1\le j\le M,$

(6.13) $w_i^0=u_0(x_i),0\le i\le N.$

Wir können das diskrete System (6.11)-(6.12) kompakt als lineares Gleichungssystem auffassen:

(6.14) $A{w}^{j+1}=B{w}^j+d^j,w^j:=(w_1^j,\dots ,w_{N-1}^j{)}^T,j=0,\dots ,M-1.$

$A,B\in {ℝ}^{(N-1)\times (N-1)}$ sind Tridiagonalmatrizen, $d^j\in {ℝ}^{N-1}$ ein Vektor, der die Randbedingungen enthält. Der Startvektor ist durch die Anfangswerte in (6.12) gegeben.

Für $\theta =0$ ist $A$ die Einheitsmatrix und die Approximationen $w^{j+1}$ können direkt berechnet werden. Daher nennt man das Differenzenverfahren mit $\theta =0$ ein explizites Verfahren. Für $\theta >0$ erhalten wir implizite Verfahren, bei denen $w^{j+1}$ nur durch Lösen eines linearen Gleichungssystems berechnet werden kann. Das Verfahren mit $\theta =1$ heißt (voll) implizit. Ist $\theta =\frac{1}{2}$, so sprechen wir vom Crank-Nicolson-Verfahren und für allgemeine Werte von $\theta >0$ von einem $\theta$-Verfahren.

Wie können die linearen Gleichungssysteme (6.14) gelöst werden? Eine sehr effiziente Methode ist das LR-Verfahren. Eine einfache Rechnung zeigt, dass die Matrix $A$ die Zerlegung $A=L\cdot R$ besitzt. Die Koeffizienten sind rekursiv definiert, das lineare Gleichungssystem $A{w}^{j+1}=B{w}^j+d^j$ ist dann äquivalent zu

$Ly=B{w}^j+d^j,R{w}^{j+1}=y$

und beide Systeme können einfach aufgelöst werden.

Bemerkung:

Im allgemeinen ist das Cholesky-Verfahren für lineare Gleichungssysteme um den Faktor 2 effizienter als die LR-Zerlegung, für Systeme mit tridiagonaler Matrix gilt dies allerdings nicht.

$\begin{array}{ccc}\text{Operation}& \text{LR-Zerlegung}& \text{Cholesky}\\ \text{Addition/Subtraktion}& 3(N-1)& 3(N-1)\\ \text{Multiplikation}& 3(N-1)& 3(N-1)\\ \text{Division}& 2N-1& 3N-1\\ \text{Wurzelauswertung}& 0& N\\ \text{Summe}& 8N-7& 10N-7\end{array}$

Damit die Folge linearer Gleichungssysteme eine sinnvolle Approximation liefert, müssen zwei Eigenschaften erfüllt sein:

• Für jedes feste, aber beliebige $j=0,\dots ,M-1$ sollte eine eindeutige Lösung $w^j$ exisiteren.
• Die approximierte Lösung $w_i^j$ sollte gegen die Lösung $u(x_i,t_j)$ der Differentialgleichung (6.9) konvergieren, wenn die Diskretisierungsparrameter $s$ und $h$ gegen Null gehen.

Ein Hauptresultat der Numerik linearer Finite-Differenzen-Schemata besagt, dass Konvergenz äquivalent zu Konsistenz und Stabilität des numerischen Schemas ist. Ein Differenzenverfahren heißt konsistent von der Ordnung $O(s^a+h^b)$, wenn es eine Konstante $C>0$ gibt, so dass für alle genügend kleinen Parameter $s,h>0$ gilt

$\underset{i,j}{\max }|L_{h,s}(u_i^j)|\le C(s^a+h^b),$

wobei $L_{h,s}$ der Differenzenoperator ist:

$L_{h,s}(u_i^j)=\frac{1}{s}(u_i^{j+1}-u_i^j)-\frac{1-\theta }{h^2}(u_{i+1}^j-u_i^j+u_{i-1}^j)-\frac{\theta }{h^2}(u_{i+1}^{j+1}-2u_i^j+u_{i-1}^{j+1}).$

Für die Stabilität eines Differenzenverfahrens gibt es verschiedene Definitionen. Wir nennen das beschriebene Verfahren stabil, wenn alle Eigenwerte der Matrix $A^{-1}B$ betragsmäßig kleiner als Eins sind.

Bezüglich Existenz und Eindeutigkeit der Lösung verweisen wir auf eine Fortgeschrittenen-Vorlesung über numerische Methoden für Differentialgleichungen oder auf die zahlreiche Literatur.

Wir weisen zum Abschluss der Betrachtungen noch auf einen interessanten Zusammenhang hin, nämlich den

Zusammenhang mit der Binomialmethode: In Abschnitt 2.5 haben wir gezeigt, dass der mit der Binomialmethode berechnete Preis einer europäischen Option im Grenzwert verschwindender Zeitschrittweiten gegen die Lösung der Black-Scholes-Differentialgleichung konvergiert. Wir zeigen nun ein in gewisser Weise umgekehrtes Resultat: Das explizite Differenzenverfahren kann als Binomialmethode interpretiert werden. Dazu betrachten wir die Black-Scholes-Gleichung (4.8), in der wir die Transformationen $x=\ln (S/K)$ und $v(x,t)=V(S,t)$ durchführen. Aus

$S{V}_S=S\frac{\partial v}{\partial x}\frac{dx}{dS}=v_x,$

$S^2{V}_{SS}=S(S{V}_S{)}_S-S{V}_S=v_{xx}-v_x$

folgt

(6.15) $v_t+\frac{{\sigma }^2}{2}{v}_{xx}+(r-\frac{{\sigma }^2}{2})v_x-rv=0,\sigma \in ℝ,0<t<T.$

Den Term $v_{xx}$ diskretisieren wir mit zentralen finiten Differenzen, $v_x$ einseitig (vorwärts oder rückwärts) oder zentral. Eine Variante der expliziten Finite-Differenzen-Approximation von (6.15) lautet also

(6.16) $\frac{v_i^j-v_i^{j-1}}{s}+\frac{{\sigma }^2}{2}\frac{v_{i+1}^j-2v_i^j+v_{i-1}^j}{h^2}+(r-\frac{{\sigma }^2}{2})\frac{v_{i+1}^j-v_{i-1}^j}{2h}-r{v}_i^{j-1}=0$

oder nach $v_i^{j-1}$ aufgelöst:

$v_i^{j-1}=\frac{1}{1+rs}[(\frac{{\sigma }^{2}s}{2h^2}-\frac{{\sigma }^{2}s}{4h}+\frac{rs}{2h})v_{i+1}^j+(1-\frac{{\sigma }^{2}s}{2h^2})v_i^j+(\frac{{\sigma }^{2}s}{2h^2}-\frac{{\sigma }^{2}s}{4h}+\frac{rs}{2h})v_{i-1}^j].$

Man beachte, dass hier rückwärts gerechnet wird: $v_i^j$ ist gegeben und $v_i^{j-1}$ wird berechnet.

Die Identifikation dieses Finite-Differenzen-Schemas mit der Binomialmethode basiert auf der Bedingung

$\frac{{\sigma }^2}{2}\frac{s}{h^2}=\frac{1}{2}.$

Man kann nun zeigen, dass das explizite Verfahren stabil ist, wenn $\alpha =({\sigma }^2/2)s/h^2\le \frac{1}{2}$ gilt, wobei in unserem Fall der Diffusionskoeffizient den Wert ${\sigma }^2/2$ besitzt. Die obige Annahme kann also als eine Stabilitätsbedingung interpretiert werden. Mit dieser Voraussetzung ergibt sich für das Schema:

(6.17) $v_i^{j-1}=\frac{1}{1+rs}[a_1{v}_{i+1}^j+a_2{v}_{i-1}^j],$

wobei

$a_{1/2}=\frac{1}{2}\pm \frac{r-{\sigma }^2/2}{2\sigma }\sqrt{s}.$

Der Preis einer europäischen Option lautet gemäß der Binomialmethode im Ein-Perioden-Modell (siehe (2.1)):

$V_i^{j-1}=e^{-rs}(p{V}_{i+1}^j+(1-p)V_{i-1}^j),p=\frac{e^{rs}-d}{u-d}.$

Setzen wir wie in Abschnitt 2.5 $u=e^{\sigma \sqrt{s}}$ und $d=1/u$, so folgt aus $e^x=1+x+𝒪(x^2)$ für $x\to 0$:

$p=\frac{e^{rs}-e^{-\sigma \sqrt{s}}}{e^{x-\sigma \sqrt{s}}-e^{-\sigma \sqrt{s}}}=\frac{(1+rs)-(1-\sigma \sqrt{s}+{\sigma }^{2}s/2)+𝒪(s^{3/2})}{(1+\sigma \sqrt{s}+{\sigma }^{2}s/2)-(1-\sigma \sqrt{s}+{\sigma }^{2}s/2)𝒪(s^{3/2})}$

$=\frac{1}{2}+\frac{(r-{\sigma }^2/2)\sqrt{s}}{2\sigma }+𝒪(s)(s\to 0),$

denn $(a+s)/(b+s)=a/b+𝒪(s)$ $(s\to 0)$. Mit der obigen Definition von $a_{1/2}$ folgt

$p=a_1+𝒪(s),1-p=a_2+𝒪(s)(s\to 0)$

und wegen $e^{-rs}=1/(1+rs)+𝒪(s^2)$

$V_i^{j-1}=\frac{1}{1+rs}[a_1{V}_{i+1}^j+a_2{V}_{i-1}^j]+𝒪(s).$

Diese Gleichung ist bis auf einen Fehler der Ordnung $𝒪(s)$ gleich der Finite-Differenzen-Approximation (6.17). Wir haben bewiesen:

Proposition 6.1

Seien

$\frac{{\sigma }^2}{2}\frac{s}{h^2}=\frac{1}{2},u=e^{\sigma \sqrt{s}}$ und $d=e^{-\sigma \sqrt{s}}$.

Dann kann das explizite Finite-Differenzen-Verfahren (6.16) bis auf einen Fehler der Ordnung $𝒪(s)$ als Binomialmethode interpretiert werden.

6.3 Beispiel: Power-Optionen

Ein European-capped-symmetric-power call ist durch die Auszahlfunktion

$C(S,T)=\min [L,((S-K{)}^{+}{)}^p]$ mit $L>0,p>0$

definiert. Der Zusatz ”capped” zeigt an, dass die Auszahlung durch die Schranke $L$ nach oben begrenzt ist. Es gibt übrigens auch unsymmetrische Power-Calls; diese haben den Payoff $(S^p-K{)}^{+}$. Ein Spezialfall ist $p=1$ und $L=1$; dieser entspricht dem europäischen Plain-Vanilla-Call.

Da wir eine europäische Power-Option betrachten, ist der Optionswert $C(S,t)$ durch die Lösung der Black-Scholes-Gleichung mit Endbedingung sowie mit geeigneten Randbedingungen für $S=0$ und für $S\to \mathrm{\infty }$ gegeben. Die Suche geeigneter Randbedingungen liefert:

Der Power-Call ist wertlos, wenn $S=0$ ist, d.h. es gilt

$C(0,t)=0$ für $0<t<T$.

Anstelle von $L$ kann auch der dikontierte Wert $L{e}^{-r(T-t)}$ vorgeschrieben werden. Nun wird die Black-Scholes-Gleichung mit den angegebenen Rand- und Anfangsbedingungen gelöst.

Resultate findet man mit dem Matlab-Programm powercall.m

% Berechnung eines europaeischen Power-Calls
% Parameter und Abkuerzungen
K = 100; r = 0.04; sigma = 0.2; T = 1; L = 150; p = 1.2;
a = 3; theta = 0.5; N = 100; M = 20; T0 = sigma^2*T/2;
h = 2*a/N; % Ortsschrittweite
s = T0/M; % Zeitschrittweite
alpha = s/h^2; k = 2*r/sigma^2;
x = [−a:h:a]; S = K*exp(x);
% Anfangswerte
C0 = min(L,max(0,S−K).^p);
u0 = exp(0.5*(k−1)*x).*C0/K;
% LR-Zerlegung von A
A = (2*alpha*theta+1)*eye(N−1);
for i=1:N−2
  A(i,i+1) = −alpha*theta;
  A(i+1,i) = −alpha*theta;
end
[LL,RR] = lu(A);
% Loesung der Differentialgleichung
u = u0;
for t=s:s:T0
  b(1:N−1) = alpha*(1−theta)*(u(3:N+1)+u(1:N−1)) . . .
  − (2*alpha*(1−theta)−1)*u(2:N);
  y = LL\b’;
  u(2:N) = RR\y;
  if ((t >= T0/2) & (t <= T0/2+s))
    u1 = u;
  end
end
% Ruecktransformation
C1 = K*exp(−0.25*(k+1)^2*T0/2)*exp(−0.5*(k−1)*x).*u1;
C2 = K*exp(−0.25*(k+1)^2*T0)*exp(−0.5*(k−1)*x).*u;
plot(S,C0,’-.’), hold on
plot(S,C1,’--’), plot(S,C2)

6.4 Amerikanische Optionen

Amerikanische Optionen können zu einem beliebigen Zeitpunkt innerhalb der Laufzeit ausgeübt werden. Wegen dieses zusätzlichen Rechtes sind sie im allgemeinen teurer als europäische Optionen. Da der Ausübungszeitpunkt vorab nicht bekannt ist, muss dieser als sog. freier Rand zusätzlich bestimmt werden, zu lösen ist im mathematischen Sinne eine ”Freie Randwertaufgabe”.

Wir untersuchen eine amerikanische Put-Option vom Preis $V_P^{am}(S,t)$ und vergleichen diese mit der entsprechenden europäischen Option $V_P^{eu}(S,t)$. Bekanntlich gilt

$V_P^{am}(S_t,t)\ge (K-S_t{)}^{+}$ für $0\le t\le T$.

Aus $V_P^{eu}(0,t)=K{e}^{-r(T-t)}$ und der Stetigkeit von $V_P^{eu}(\cdot ,t)$ in einer Umgebung von $S=0$ folgt dann

$V_P^{eu}(S,t)<K-S=(K-S{)}^{+}\le V_P^{am}(S,t)$

für hinreichend kleines $S>0$ mit $r>0$ und $t<T$.

Wir definieren nun einen Kontaktpunkt $S_{f}wiefolgt:::(6.18)<math>V_P^{am}(S_f,t)=K-S_f<math>(0<S_f<K).Wirbemerken,dass<math>S_f$ von der Zeit $t$ abhängt, d. h. $S_f=S_f(t)$. Charakterisiert wird der Kontaktpunkt $S_f(t)$ durch die Beziehungen

(6.19) $\{\begin{array}{cc}{V}_P^{am}(S,t)>(K-S{)}^{+}& \mathrm{f}\ddot{\mathrm{u}}\mathrm{r}S>S_f(t),\\ V_P^{am}(S,t)=(K-S{)}^{+}& \mathrm{f}\ddot{\mathrm{u}}\mathrm{r}S\le S_f(t).\end{array}$

Die Lage der Mannigfaltigkeit $S_f(t),t\in (0,T)$, ist zunächst unbekannt. Den Stoppzeitpunkt $t_s$ und damit die Lage des Kontaktpunktes bestimmt der Halter der Option durch seine Ausübungsentscheidung; mathematisch ist der Schnittpunkt zwischen einem Pfad und der Lösung (6.19) ausschlaggebend. Daraus wird eine Randbedingung gewonnen.

Randbedingung: Der Anstieg $\frac{\partial V}{\partial S}$, mit dem die Funktion $V_P^{am}(S,t)$ die Gerade $K-S$ tangential berührt, nimmt den folgenden Wert an:

(6.20) $\frac{\partial V}{\partial S}(S_f(t),t)=\frac{\partial }{\partial S}(K-S)=-1.$

Andererseits gilt für $S\to \mathrm{\infty }$ die Abklingbedingung. Dazu gilt im Gebiet für $S_f(t)<S<\mathrm{\infty }$ und in $0<t<t_s\le T$ die Differentialgleichung

(6.21) $\frac{\partial V}{\partial t}+\frac{{\sigma }^2}{s}{S}^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV(S,t)=0.$

Bemerkung:

Mit $q$ wird eine (stetige) Dividendenausschüttung bezeichnet. Man kann zeigen, dass im Fall eines Calls und für $q=0$ eine vorzeitige Ausübung der Option keinen Gewinneffekt liefert.

Das beschriebene Problem kann auch als Black-Scholes-Ungleichung beschrieben werden. Wir geben die Beziehungen nur für den Fall einer Put-Option an:

$V_P^{am}=K-S$ für $S\ge S_f$ (stop)

$V_P^{am}$ löst Gleichung (6.21) für $S<S_f$ (hold)

Eine analytische Lösung des asymptotischen Problems für $T\to \mathrm{\infty }$ wurde in den Jahren 1986/87 von MacMillan, Barone-Adesi und Whaley (MBAW-Methode) konstruiert. Später (1992) wurde eine genauere Lösung unter Verwendung der Laplace-Transformation und einer quadratischen Approximation der Lösung der verbleibenden gewöhnlichen Differentialgleichung von Carr, Jarrow und Myneni (CJM-Methode) beschrieben.

Wir gehen kurz auf die MBAW-Methode ein:

Wegen $V_P^{am}(S,t)\ge V_P^{eu}(S,t)$ können wir eine Funktion $\rho (s,t)\ge 0$ definieren, so dass gilt

(6.23) $V_P^{am}(S,t)=V_P^{eu}(S,t)+\rho (S,t).$

Für $t\to \mathrm{\infty }$ ist die zeit-unabhängige Gleichung

(6.24) $\frac{{\sigma }^2}{2}{S}^2\frac{{\partial }^2\rho }{\partial S^2}+(r-q)S\frac{\partial \rho }{\partial S}-r\rho =0,$

die als Eulersche Differentialgleichung behandelt werden kann, zu lösen. Wir setzen der Einfachheit halber $q=0$ und erhalten für eine der Lösungen, die für $S\to \mathrm{\infty }$ die Eigenschaft $|\rho |\to 0$ besitzt, den Ausdruck

(6.25) $\psi (S)=A_1{\left(\frac{S}{S_f}\right)}^{-k_2},k_2=\frac{2r}{{\sigma }^2},A_1-$ bel. Konstante.

Nun gilt auf dem freien Rand die Stetigkeitsbedingung für den Optionspreis und wir errechnen mittels der oben formulierten Randbedingung

(6.26) $\frac{\partial V_P^{am}(S_f)}{\partial S}=\frac{\partial V_P^{eu}(S_f)}{\partial S}+{\psi }^{\prime }(S_f)=-1,$

so dass sich eine Bestimmungsgleichung für $A_1$ mit der Bezeichnung $\mathrm{\Delta }{|}_{S=S_f}$ für die Ableitung der Lösung der europäischen Put-Option (Greek ”Delta”) gerade zu

(6.27) $A_1=\frac{(1+\mathrm{\Delta }{|}_{S=S_f})}{k_2}$

ergibt. Zum Beispiel lassen sich mit Hilfe eines Formelmanipulationssystems (Maple, Mathematica, ...) der Wert $S_f$ und schließlich auch die Lösungen $\rho (S)$ und $V_P^{am}(S,t)$ des nichtlinearen Gleichungssystems im asymptotischen Fall $t\to \mathrm{\infty }$ bestimmen.

Eine Abänderung des Modells durch Einführung bzw. Abspaltung einer asymptotisch fallenden Funktion $h(r,t)=1-\exp (-rt)$ und für die Zeit $T=1$ führt auf kompliziertere, nichtlineare Ausdrücke. Die Zahlenrechnung ergibt dann

AmericanPutMBAW[100, 100, 0.3, 0.03, 0, 1] = 10.6010
BlackScholesPut[100, 100, 0.3, 0.03, 0, 1] = 10.3279

Zum Vergleich erhält man bei 5-jähriger Laufzeit:

AmericanPutMBAW[100, 100, 0.3, 0.03, 0, 5] = 20.2494
AmericanPutCJM[100, 100, 0.3, 0.03, 0, 5] = 20.0304

(Methode von Carr, Jarrow, Myneni)

Auf eine Beschreibung von Differenzenmethoden (FDM - Finite Difference Methods) sei an dieser Stelle verzichtet, obwohl dies eine praktisch gut realisierbare Methode zur approximativen Berechnung freier Randwertprobleme ist.

6.5 Vertikale Linienmethode

(– numerische Methode zur Abspaltung des Zeit-Anteils einer parabolischen, partiellen Differentialgleichung, siehe unter dem Stichwort ”Numerik partieller Differetialgleichungen” –)

6.6 Beispiel: Basket-Optionen

– Untersuchung der Zusammenfassung mehrerer Aktien zu einem Aktienpaket (”Basket”). Bei einer größeren Anzahl werden meist stochastische Pfadsimulationen angewendet.