Kurs:Modellierung und Numerische Methoden von Finanzderivaten/5 Die Monte-Carlo-Methode

Der Preis einer europäischen (plain vanilla) Option kann mit der Black-Scholes-Formel aus Kapitel 4 berechnet werden. Leider existieren zu komplexeren Optionen im allgemeinen keine expliziten Formeln mehr. In diesem Kapitel stellen wir die Monte-Carlo-Methode zur Integration von stochastischen Differentialgleichungen vor, mit der faire Preise von komplizierten Optionen numerisch berechnet werden können. Zuerst führen wir in Abschnitt 5.1 in die Thematik ein. Das Monte-Carlo-Verfahren erfordert die Simulation von Realisierungen eines Wiener-Prozesses. Die Simulation wiederum benötigt normalverteilte Zufallszahlen. Die Erzeugung von Zufallszahlen ist Gegenstand von Abschnitt 5.2. In Abschnitt 5.3 erläutern wir die numerische Lösung stochastischer Differentialgleichungen. Die Präzision von Monte-Carlo-Simulationen kann mit Hilfe der Technik der Varianzreduktion, die wir in Abschnitt 5.4 vorstellen, erhöht werden. Schließlich wenden wir die vorgestellten Methoden in Abschnitt 5.5 zur Simulation einer asiatischen Call-Option mit stochastischer Volatilität an.

5.1 Grundzüge der Monte-Carlo-Simulation

Die Berechnung des fairen Preises einer komplexen Option ist im Allgemeinen eine anspruchsvolle Aufgabe, die nur numerisch gelöst werden kann. Bei vielen Optionen ist es notwendig, stochastische Differentialgleichungen bzw. stochastische Integrale numerisch zu lösen. Beispiele für derartige Situationen stellen wir im Folgenden vor.

Beispiel 5.1: Asiatischer Call im Heston-Modell

Berechne den fairen Preis einer asiatischen Call-Option mit der Auszahlungsfunktion

$V(S_T,T)={(S_T-\frac{1}{T}\underset{0}{\overset{T}{\int }}{S}_{\tau }d\tau )}^{+},$

wobei die Dynamik von $S_t$ und ${\sigma }_t$ durch das Heston-Modell

$dS=\mu Sdt+\sigma Sd{W}_1,d{\sigma }^2=\kappa (\theta -{\sigma }^2)dt+\nu \sigma d{W}_2$

gegeben sei. Dieses Beispiel erfordert die Integration von stochastischen Differentialgleichungen und des stochastischen Integrals $\underset{0}{\overset{T}{\int }}{S}_{\tau }d\tau$. Die Integration kann mit Hilfe der Monte-Carlo-Methode durchgeführt werden. Wir erläutern dies im Detail in Abschnitt 5.5.

Beispiel 5.2: Basket-Option

Berechne den fairen Preis einer europäischen Option auf $n$ Aktien (Basket-Option) mit der Auszahlungsfunktion $V_0(S_1,\dots ,S_n)$. Ist die Dynamik der Aktienkurse $S_1(t),\dots ,S_n(t)$ wie in Abschnitt 4.5 durch

$d{S}_i={\mu }_i{S}_{i}dt+{\sigma }_i{S}_{i}d{W}_i,i=1,\dots ,n$

und durch eine mehrdimensionale Brownsche Bewegung $(W_1,\dots ,W_n)$ mit der Kovarianz-Matrix $\mathrm{\Sigma }$ (siehe Definition 5.7) gegeben, so berechnet sich der Optionspreis nach dem Black-Scholes-Modell analog zu (4.18) nach

$V(S_1,\dots ,S_n,t)=\frac{e^{-r(T-t)}}{\sqrt{(\det \mathrm{\Sigma })(2\pi (T-t){)}^n}}\times \underset{0}{\overset{\mathrm{\infty }}{\int }}\dots \underset{0}{\overset{\mathrm{\infty }}{\int }}{e}^{-T{\mathrm{\Sigma }}^{-1}/2}{V}_0(S{\text{'}}_1,\dots ,S{\text{'}}_n)\frac{d{S_1}^{\prime }\dots d{S_n}^{\prime }}{{S_1}^{\prime }\dots {S_n}^{\prime }}$

mit $\mathrm{\Sigma }=({\mathrm{\Sigma }}_{ij}{)}_{i,j=1}^n,\alpha =({\alpha }_1,\dots ,{\alpha }_n)$ und

${\alpha }_i:=-\frac{\ln (S{\text{'}}_i/S_i)-(r-{\mathrm{\Sigma }}_{ii}^2/2)(T-t)}{\sqrt{T-t}}.$

Hier muss also ein $n$-dimensionales (Riemann-)Integral ausgewertet werden, wobei die Dimension n je nach Größe des Baskets sehr groß sein kann. Numerische Quadraturformeln sind hier ungeeignet, da zu viele Funktionswerte ausgewertet werden müssen (bei $n=100$ z. B. $2^{100}\approx 10^{30}$ Auswertungen). Einen Ausweg bietet die Monte-Carlo-Integration. Ein Beispiel, wie das geht, geben wir in Kap. 5.2, Beispiel 5.7.

Zur Vereinfachung betrachten wir im folgenden eine europäische Plain-Vanilla-Put-Option auf einen Basiswert, dessen Kurs sich gemäß einer geometrischen Brownschen Bewegung entwickelt:

(5.1) $d{S}_t=r{S}_{t}dt+\sigma S_{t}d{W}_t$

mit Anfangswert $S_0$, konstantem risikofreien Zinssatz $r\ge 0$, konstanter Volatilität $\sigma >0$ und einem Wiener-Prozess $W_t$ (siehe Abschnitt 3.2). In Bemerkung 4.10 haben wir gezeigt, dass der Optionspreis $V(S_t,t)$ zur Zeit $t=0$ gegeben ist durch den diskontierter Erwartungswert

(5.2) $V(S_0,0)=e^{-rT}E(V(S_T,T)).$

Die Grundidee der Monte-Carlo-Simulation besteht nun darin, den Erwartungswert in (5.2) durch Simulation von $N$ Pfaden $S_t:0<t<T$ des Basiswert-Kurses zu approximieren. Der Algorithmus besteht aus vier Schritten:

• Simulation der Basiswert-Pfade: Berechne für $k=1,\dots ,N$ den Itô-Prozess (5.1) zum Anfangswert $S_0$ mit den Lösungen $(S_t{)}_k$.

• Berechnung der Auszahlungsfunktion: Bestimme für alle $k=1,\dots ,N$ die Auszahlungsfunktion entsprechend zum Pfad $(S_t{)}_k$:

$V_{T,k}:=(K-(S_T{)}_k{)}^{+}.$

• Berechnung eines Schätzers: Berechne einen Schätzer für den Erwartungswert in (5.2). Naheliegend ist etwa die Wahl

$\hat{E}(V_T):=\frac{1}{N}\sum _{k=1}^N{V}_{T,k},$

wobei $V_T=(V_{T,1},\dots ,V_{T,N})$.

• Bestimmung des Optionspreises: Berechne eine Approximation des fairen Optionspreises durch

$\hat{V}:=e^{-rT}\hat{E}(V_T).$

Die Schritte 2-4 bereiten keine Schwierigkeiten. Schritt 3 beruht darauf, dass nach dem Gesetz der großen Zahlen das arithmetische Mittel von gleichverteilten und unabhängigen Zufallsvariablen fast sicher gegen den Erwartungswert konvergiert (siehe z.B. [3]). Schritt 1 benötigt die numerische Integration von stochastischen Differentialgleichungen, die aus zwei Teilaufgaben besteht:

• Simulation von $N$ unabhängigen Realisierungen eines Wiener-Prozesses und
• approximative Berechnung der Lösung der stochastischen Differentialgleichung zum jeweiligen Pfad des Wiener-Prozesses.

Eine sehr einfache Approximation der Gleichung (5.1) ist gegeben durch

$(5.3)\mathrm{\Delta }{S}_t=r{S}_t\mathrm{\Delta }t+\sigma S_t\mathrm{\Delta }{W}_t,$

wobei $\mathrm{\Delta }{W}_t=W_{t+\mathrm{\Delta }t}-W_t$ $𝒩(0,\mathrm{\Delta }t)$-verteilt ist (siehe Satz 3.10 (3)). Wir benötigen nun Realisierungen des Wiener-Prozesses $\mathrm{\Delta }{W}_t$. Wegen

$\mathrm{\Delta }{W}_t=Z\cdot \sqrt{\mathrm{\Delta }t}$

genügt es, Realisierungen einer $𝒩(0,1)$-verteilten Zufallsgröße $Z$ zu bestimmen. Das folgende Matlab-Programm simuliert einen Wiener-Prozess.

% Simulation eines Wiener-Prozesses
h = 0.01; W(1) = 0;
for i=1:999
  W(i+1) = W(i) + randn*sqrt(h);
end

Die Matlab-Funktion randn liefert standardnormalverteilte Zufallszahlen. Tatsächlich handelt es sich um Pseudo-Zufallszahlen, da der Algorithmus zur Erzeugung dieser Zahlen deterministisch ist. Wir können die Approximation (5.3) auch schreiben als

(5.4) $S_{k+1}-S_k=\mathrm{\Delta }{S}_k=r{S}_k\mathrm{\Delta }t+\sigma S_{k}Z\sqrt{\mathrm{\Delta }t},S_0\text{ gegeben},k=0,1,2,\dots$

% Approximative Loesung von (5.1)
h = 0.01; mu = 0.1; sigma = 0.4;
S(1) = 1; 
for i=1:999
  dW(i) = randn*sqrt(h);
  S(i+1) = S(i)*(1 + mu*h + sigma*dW(i));
end

Damit sind wir in der Lage, unsere erste Monte-Carlo-Simulation gemäß der obigen vier Schritte durchzuführen. Die Realisierung erfolgt mit einem Matlab-Programm.

Der Funktionsaufruf randn(’state’,3) bedeutet, dass der Pseudo-Zufallszahlengenerator von Matlab mit der Zahl 3 initialisiert wird. Dies hat den Zweck, die Simulationsergebnisse reproduzierbar zu machen. In Abbildung 5.12 illustrieren wir die Entwicklung der Preise $\hat{V}$ einer europäischen Put-Option in Abhängigkeit der Anzahl $N$ der Monte-Carlo-Simulationen. Der Black-Scholes-Optionspreis beträgt $V=16.98$. Der Monte-Carlo-Preis weicht von dem Black-Scholes-Preis stark ab, wenn die Anzahl $N$ der Monte-Carlo-Simulationen zu klein gewählt wurde. Allerdings schwanken die Werte auch für große $N$ noch recht stark. Natürlich ist es für dieses Beispiel wesentlich effizienter, die Black-Scholes-Formel zur Bestimmung des Optionspreises zu verwenden. Für komplexere Optionen wie die asiatische Option im Heston-Modell aus Beispiel 5.1 sind wir jedoch auf die Monte-Carlo-Methode angewiesen, da keine expliziten Formeln existieren.

% Monte-Carlo-Simulation fuer einen europaeischen Put
clear all, randn(’state’,3)
K = 100; r = 0.05; sigma = 0.2; T = 1;
n = 50; h = 1/n; S(1) = 80;
for j=1:100
  N=j*100;
  for k=1:N
    for i=2:n
      S(i) = S(i−1)*(1 + r*h + sigma*randn*sqrt(h)); S
    end
    payoff(j,k) = max(0,K−S(n));
  end
  V(j) = exp(−r*T)*sum(payoff(j,:))/N;
end
plot(V)

Die oben präsentierten Beispiele und Simulationen geben Anlass zu den folgenden Fragen:

• Wie können wir standardnormalverteilte Pseudo-Zufallszahlen erzeugen?
• Wie genau ist die Approximation von Gleichung (5.4)? Wie kann die Approximation verbessert werden?
• Wie kann das hochdimensionale Integral aus Beispiel 5.2 mittels der Monte-Carlo-Methode approximiert werden?

Diese Fragen werden wir in den nächsten Abschnitten beantworten.

5.2 Pseudo-Zufallszahlen

Für die Simulation des Wiener-Prozesses benötigen wir standard-normalverteilte Zufallszahlen $Z$, um die Inkremente $\mathrm{\Delta }W=\sqrt{\mathrm{\Delta }t}Z$ zu berechnen. Wir benutzen dafür die Notation

$Z\sim 𝒩(0,1).$

Erzeugen wir Zufallszahlen im Rechner, so handelt es sich letztlich immer um eine deterministische Vorgehensweise. Man spricht daher von Pseudo-Zufallszahlen. Im Folgenden benutzen wir jedoch den Begriff Zufallszahlen auch, wenn wir Pseudo-Zufallszahlen meinen. Zuerst erzeugen wir im Intervall $[0,1]$ gleichverteilte (Pseudo-)Zufallszahlen $Y$,

$Y\sim 𝒰[0,1],$

und transformieren sie dann auf normalverteilte Zufallszahlen:

$Z:=h(Y)\sim 𝒩(0,1).$

Um die obigen Begriffe zu präzisieren, hier eine Definition:

Definition 5.1

1. Eine Zufallsvariable $X$ heißt gleichverteilt auf dem Intervall $[a,b]$ (in Zeichen: $X\sim 𝒰[a,b]$), wenn sie die Dichtefunktion $f(x)=\frac{1}{(b-a)},x\in [a,b]$ besitzt.
2. Eine Folge von Zufallszahlen heißt nach $F$ verteilte Zufallszahlen, wenn sie unabhängige Realisierungen einer nach einer Verteilungsfunktion $F$ verteilten Zufallsvariablen sind.

Gleichverteilte Zufallszahlen:
Ein einfacher Algorithmus, um auf $[0,1]$ gleichverteilte Zufallszahlen zu erzeugen, ist durch die lineare Kongruenz-Methode gegeben: Seien $M\in \mathbb{N},a,b,X_0\in \{0,\dots ,M-1\}$.

Für $k=1,2,\dots$

$X_k:=(a{X}_{k-1}+b)modM,U_k:=X_k/M.$

Offenbar müssen wir $a=0$ und (wenn $b=0$) $X_0=0$ ausschließen. Außerdem sollte $a\ne 1$ sein, denn ansonsten wäre $X_k=(X_0+kb)modM$ zu leicht vorhersagbar. Die Folge $X_k$ hat die folgenden Eigenschaften:

• Die Folge $(X_k{)}_{k\in \mathbb{N}}$ ist periodisch mit einer Periode, die kleiner oder gleich $M$ ist, denn: Wegen $X_k\in \{0,\dots ,M-1\}$ muss es ein $p\in \mathbb{N}$ geben, so dass $X_p=X_0$ und daher $X_k=X_{k+p}$ für alle $k\in \mathbb{N}$ ist.
• Die Verteilung der Zufallsvektoren $(U_k,\dots ,U_{k+m})$ ist leider korreliert (also nicht unabhängig voneinander); siehe das folgende Beispiel.

Beispiel 5.3

Wir betrachten den Fall $m=2$ mit den Daten $M=2048,a=1229,b=1$ und $X_0=1$. Die Punkte liegen auf parallelen Geraden. Solche Zahlen können wir kaum Zufallszahlen nennen!

In dem folgendem Matlab-Programm ist die lineare Kongruenz-Methode implementiert.

% Pseudo-Zufallszahlen nach der linearen Kongruenz-Methode
a = 1229; b = 1; M = 2048; N = 500;
X(1) = 1; 
for i = 2:N
  X(i) = mod(a*X(i−1)+b,M);
end
plot(X([1:N−1]),X([2:N]),’.’)

Wegen der Eigenschaft, dass die Zufallszahlen auf parallelen Geraden liegen, ist die lineare Kongruenz-Methode nicht sehr brauchbar. Besser sind sogenannte Fibonacci-Generatoren geeignet. Die Idee ist hier, die Fibonacci-Folge zu verwenden:

Für $k=3,4,\dots$

$X_k:=(X_{k-1}+X_{k-2})modM,U_k:=X_k/M$

mit $M,X_1,X_2\in \mathbb{N}$. Je nach Wahl von $M$ können aber die Ergebnisse recht unbefriedigend sein. Es sind weit weniger als 2000 Punkte zu sehen, da die Folge $(U_k)$ sich wiederholt.

Geeigneter sind sogenannte lagged Fibonacci-Generatoren (oder Fibonacci-Generatoren mit ”Verzögerung”) der Form

Für $k\ge \max \{\mu ,\nu \}$

$X_k:=(X_{k-\mu }-X_{k-\nu })modM,$

if $X_k<0$ then $X_k:=X_k+M$;

$U_k:=X_k/M,$

wobei die Anfangszahlen $X_1,\dots ,X_{\max \{\mu ,\nu \}}$ etwa mittels einer linearen Kongruenz-Methode bestimmt werden können. In dem folgenden Matlab-Programm werden diese Zahlen allerdings mittels des bereits implementierten Zufallszahlen-Programms rand berechnet.

% Pseudo-Zufallszahlen nach dem lagged Fibonacci-Generator
rand(’state’,2) 
M = 2048; nu = 17; mu = 5; N = 5000;
X = M*rand(1,max(nu,mu));
for i=max(mu,nu)+1:N
  X(i) = mod(X(i−mu)−X(i−nu),M);
  U(i) = X(i);
end 
plot(U([1:N−1]),U([2:N]),’.’)

Die Punkte erscheinen genügend zufällig verteilt. Fibonacci-Generatoren haben außerdem den Vorteil, dass sie sehr einfach zu implementieren sind.

Normalverteilte Zufallszahlen:
Wir erzeugen normalverteilte Zufallszahlen durch Transformation gleichverteilter Zufallszahlen. Dies kann geschehen durch

• Invertierung der Verteilungsfunktion oder
• Transformation zwischen Zufallszahlen.

Grundlage für die Invertierung ist der folgende Satz:

Satz 5.1

Sei $U\sim 𝒰[0,1]$ eine gleichverteilte Zufallsvariable und $F$ eine stetige, streng monotone Verteilungsfunktion. Dann ist die Zufallsvariable $F^{-1}(U)$ nach $F$ verteilt.

Beweis:

Die Umkehrfunktion $F^{-1}$ existiert gemäß Voraussetzungen. Die Annahme der Gleichverteilung impliziert für alle $\xi \in [0,1]$: $P(U\le \xi )=\xi$. Somit folgt

$P(F^{-1}(U)\le x)=P(U\le F(x))=F(x).$

Dies bedeutet, dass $F^{-1}(U)$ nach $F$ verteilt ist.

q.e.d.

Ist dieser Satz auf die Normalverteilung $\varphi$ anwendbar? Es liegen weder für $\varphi$ noch für ${\varphi }^{-1}$ geschlossene Formelausdrücke vor. Die nichtlineare Gleichung $\varphi (x)=u$ müsste numerisch invertiert werden, etwa mittels des in Abschnitt 4.2 vorgestellten Newton-Verfahrens. Allerdings ist das Problem für $u\approx 1$ schlecht konditioniert (kleine Änderungen in $u$ bewirken sehr große Änderungen in $x$). Als Ausweg kann man ${\varphi }^{-1}$ ähnlich wie in Abschnitt 4.2 durch eine rationale Funktion $G$ approximieren und $x=G(u)\approx {\varphi }^{-1}(u)$ setzen. Bei der rationalen Approximation ist das asymptotische Verhalten von ${\varphi }^{-1}$ (senkrechte Tangenten bei $u=0$ und $u=1$) sowie die Punktsymmetrie zu $(u,x)=(\frac{1}{2},0)$ zu berücksichtigen.

Wir wählen allerdings die zweite Idee: Transformation der Zufallszahlen. Grundlage hierfür ist der folgende Satz.

Satz 5.2

Sei $X$ eine Zufallsvariable auf ${ℝ}^n$ mit der Dichtefunktion $f>0$ auf der Menge $S:=\{x\in {ℝ}^n:f(x)>0\}$ von $f$. Die Transformation $h:S\to B:=h(S)$ sei umkehrbar mit stetig invertierbarer Inversen $h^{-1}$. Dann hat $Y:=h(X)$ die Dichtefunktion

$y\to f(h^{-1}(y))\cdot |\det \left(\frac{d{h}^{-1}(y)}{dy}\right)|,y\in B.$

Beweis:

Wir geben nur eine grobe Beweisskizze. Nach dem Transformationssatz im ${ℝ}^n$ gilt (sei $A$ der Wertebereich von $h$, d. h. $h:D(h)\mapsto A$):

$P(Y=h(X)\in A)=P(X\in h^{-1}(A))=\underset{h^{-1}(A)}{\int }f(u)du=\underset{A}{\int }f(h^{-1}(y))\cdot |\det \left(\frac{d{h}^{-1}(y}{dy}\right)|dy.$

q.e.d.

Im Falle $n=1$ und $f(x)=1$ (Gleichverteilung in $[0,1]$) suchen wir also eine Transformation y = h(x)</math>, so dass die transformierte Dichtefunktion gleich der Normalverteilung ist:

$1\cdot \left|\frac{d{h}^{-1}(y)}{dy}\right|=\frac{1}{\sqrt{2\pi }}\exp (-\frac{y^2}{2}).$

Dies ist eine gewöhnliche Differentialgleichung für $h^{-1}$, die leider keine geschlossene Formel für die Transformation liefert. Verblüffenderweise erhalten wir eine geschlossene Formel, wenn wir nicht in $ℝ$, sondern in ${ℝ}^2$ transformieren. Das geht folgendermaßen. Wir wenden Satz 5.2 auf $S=[0,1{]}^2$ und $f(x)=1,x\in S$ an. Wähle die Transformation y = h(x)</math> mit

$h(x)=\left(\begin{array}{c}\sqrt{-2\ln x_1}\cos (2\pi x_2)\\ \sqrt{-2\ln x_1}\sin (2\pi x_2)\end{array}\right),x=(x_1,x_2{)}^T\in S.$

Die Umkehrabbildung lautet

$h^{-1}(y)=\left(\begin{array}{c}\exp \left(\frac{-|y{|}^2}{2}\right)\\ \frac{1}{2\pi }\arctan \left(\frac{y_2}{y_1}\right)\end{array}\right),y=(y_1,y_2{)}^T,|y{|}^2=y_1^2+y_2^2.$

Für die Determinante ergibt sich mit $y=h(x)$:

$\det \left(\frac{d{h}^{-1}}{dy}\right)=\det \left(\begin{array}{cc}-y_1{x}_1& -y_2{x}_1\\ \frac{1}{2\pi }\frac{-y_2/y_1^2}{1+y_2^2/y_1^2}& \frac{1}{2\pi }\frac{1/y_1}{1+y_2^2/y_1^2}\end{array}\right)=-\frac{x_1}{2\pi }=-\frac{1}{2\pi }\exp (-\frac{|y{|}^2}{2}).$

Dies ist die Dichtefunktion der Standardnormalverteilung in ${ℝ}^2$ (von zwei unabhängigen Zufallsvariablen). Also ist $h(X)$ standardnormalverteilt, falls $X$ auf $[0,1]$ gleichverteilt ist.

Daraus folgt der Algorithmus von Box-Muller: Generiere $U_1,U_2\sim 𝒰[0,1]$ und setze $\theta =2\pi U_2$ und $\rho =\sqrt{-2\ln U_1}$. Dann sind

$Z_1=\rho \cos \theta$ und $Z_2=\rho \sin \theta$

standardnormalverteilt. Das Histogramm zeigt, dass dieser Algorithmus tatsächlich normalverteilte Zufallszahlen $Z_1$ liefert. Hierfür wurden 50 000 Zufallszahlen nach dem folgenden Matlab-Programm erzeugt.

% N(0,1)-Zufallszahlen nach Box-Muller
N = 50000; rand(’state’,2)
Z = sqrt(−2*log(rand(1,N))).*cos(2*pi*rand(1,N));
x=[−3.8:0.2:3.8];
hist(Z,x)

Es muss natürlich sichergestellt werden, dass die gleichverteilten Zufallszahlen keine Struktur haben, da diese auf die transformierten normalverteilten Zufallsvariablen bertragen werden. Eine Linienstruktur in $[0,1{]}^2$ würde auf Kurven in der $(Z_1,Z_2)$-Ebene abgebildet werden.

Beim Box-Muller-Algorithmus sind drei Funktionsaufrufe (sqrt, log und cos bzw. sin) erforderlich, um zwei normalverteilte Zufallszahlen zu erhalten. Dies kann beim sogenannten Marsaglia-Algorithmus durch Verwendung der Polartransformation verbessert werden.

Für weitere Hinweise auf Techniken zur Erzeugung von Zufallszahlen verweisen wir auf verschiedene Monographien.

Korreliert normalverteilte Zufallszahlen:
Bei der Simulation einer mehrdimensionalen Brownschen Bewegung benötigen wir i. a. Zufallsvektoren, die einer korrelierten mehrdimensionalen Verteilung folgen:

$X=(X_1,\dots ,X_n)\sim 𝒩(\mu ,\mathrm{\Sigma }).$

Definition 5.2

(1) Sei $X=(X_1,\dots ,X_n)$ ein Zufallsvektor, $\mu \in {ℝ}^n$ und $\mathrm{\Sigma }$ eine symmetrische, positiv definite $(n\times n)$-Matrix. Dann heißt der Vektor $X$ mit $\mu$ und $\mathrm{\Sigma }$ normalverteilt, also $𝒩(\mu ,\mathrm{\Sigma })$-verteilt, wenn X</math> die Dichtefunktion

$f(x)=\frac{1}{\sqrt{(2\pi {)}^n\det (\mathrm{\Sigma })}}\exp (-\frac{1}{2}(x-\mu {)}^T{\mathrm{\Sigma }}^{-1}(x-\mu )),x\in {ℝ}^n$

besitzt.

(2) Sei $X=(X_1,\dots ,X_n)$ eine $𝒩(\mu ,\mathrm{\Sigma })$-verteilte Zufallsvariable. Dann heißt die Matrix $\mathrm{\Sigma }=({\mathrm{\Sigma }}_{ij})$ die Kovarianz-Matrix und es gilt

${\mathrm{\Sigma }}_{ij}=E[(X_i-{\mu }_i)(X_j-{\mu }_j)],$

wobei $\mu =({\mu }_1,\dots ,{\mu }_n{)}^T=(E(X_1),\dots ,E(X_n){)}^T$ der Erwartungswert von $X$ ist. Die Matrix $P=({\rho }_{ij})$, die aus den Elementen

${\rho }_{ij}:=\frac{{\mathrm{\Sigma }}_{ij}}{\sqrt{{\mathrm{\Sigma }}_{ii}{\mathrm{\Sigma }}_{jj}}}$

besteht, heißt die Korrelation.

Der Begriff ”korreliert” bedeutet also, dass die Korrelation einer $𝒩(\mu ,\mathrm{\Sigma })$-verteilten Zufallsvariablen keine Diagonalmatrix ist.

Wir betrachten nun einen Vektor $Z=(Z_1,\dots ,Z_n)$ aus unabhängigen, standard-normalverteilten Zufallsvariablen $Z_k$ mit der Dichtefunktion $f$. Wir konstruieren mittels $Z$ eine $𝒩(\mu ,\mathrm{\Sigma })$-verteilte Zufallsvariable $Y$. Es sei $\mu \in {ℝ}^n$ und $\mathrm{\Sigma }\in {ℝ}^{n\times n}$ symmetrisch und positiv definit $(\mathrm{\Sigma }={\mathrm{\Sigma }}^T>0)$. Dann existiert eine Cholesky-Zerlegung

$\mathrm{\Sigma }=L{L}^T.$

Wir zeigen, dass $Y=\mu +LZ$ die Eigenschaft der $𝒩(\mu ,\mathrm{\Sigma })$-Korrelation besitzt:
Sei $x=Lz$ (für die Vektoren schreiben wir $x,y,z,\dots$), dann folgt

$f(z)dz=\frac{1}{(2\pi {)}^{n/2}}\exp (-\frac{z^{T}z}{2})dz$

$=\frac{1}{(2\pi {)}^{n/2}}\exp (-\frac{(L^{-1}x{)}^T(L^{-1}x)}{2})dz$ (mit $x=Lz$)

$=\frac{1}{(2\pi {)}^{n/2}}\exp (-\frac{x^T(L{L}^T{)}^{-1}x}{2})dz$

$=\frac{1}{(2\pi {)}^{n/2}|\det L|}\exp (-\frac{x^T{\mathrm{\Sigma }}^{-1}x}{2})dx$ (weil $dx=|\det L|dz$)

$=\frac{1}{(2\pi {)}^{n/2}(\det (\mathrm{\Sigma }{)}^{1/2}}\exp (-\frac{x^T{\mathrm{\Sigma }}^{-1}x}{2})dx$

Damit ist $X\sim 𝒩(0,\mathrm{\Sigma })$ und $Y=\mu +X\sim 𝒩(\mu ,\mathrm{\Sigma })$.

Der Algorithmus lässt sich wie folgt zusammenfassen:

• Berechne die Cholesky-Zerlegung $\mathrm{\Sigma }=L{L}^T$.
• Berechne unabhängige $X_k\sim 𝒩(0,1),k=1,2,\dots ,n$ (z.B. mit dem Box-Muller-Verfahren). Setze $Z=(X_1,\dots ,X_n)$.
• Die Zufallsvariable $Y=\mu +LZ$ ist $𝒩(\mu ,\mathrm{\Sigma })$-verteilt.

Beispiel 5.4

Gesucht ist eine 2D-$𝒩(0,\mathrm{\Sigma })$-verteilte Zufallsvariable $(X_1,X_2)$. Dabei sei

$\mathrm{\Sigma }=\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right).$

Der Vektor $X=(X_1,X_2)$ heißt durch $\rho$ korreliert. Mit dem Ansatz

$L=\left(\begin{array}{cc}a& 0\\ b& c\end{array}\right)$

liefert die Cholesky-Zerlegung \Sigma = LL^T</math>:

$\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right)=\left(\begin{array}{cc}a& 0\\ b& c\end{array}\right)\cdot \left(\begin{array}{cc}a& b\\ 0& c\end{array}\right)=\left(\begin{array}{cc}{a}^2& ab\\ ab& b^2+c^2\end{array}\right).$

Durch Koeffizientenvergleich ergibt sich (man beachte, dass $|\rho |\le 1$ gilt!)

$L=\left(\begin{array}{cc}{\sigma }_1& 0\\ \rho {\sigma }_2& {\sigma }_2\sqrt{1-{\rho }^2}\end{array}\right)$

Sind $Z_1,Z_2$ unabhängig und $𝒩(0,1)$-verteilt, so ist

$\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)=L\left(\begin{array}{c}{Z}_1\\ Z_2\end{array}\right)=\left(\begin{array}{c}{\sigma }_1{Z}_1\\ {\sigma }_2(\rho Z_1+\sqrt{1-{\rho }^2}{Z}_2)\end{array}\right)$

normalverteilt mit den Parametern $\mu =0$ und $\mathrm{\Sigma }$ wie beschrieben.

Beispiel 5.5:

Gesucht ist eine dreidimensionale Verteilung $(X_1,X_2,X_3)$, die normalverteilt sein soll, den Erwartungswert $\mu =(-5,0,10{)}^T$ und die Kovarianzmatrix

$\mathrm{\Sigma }=\left(\begin{array}{ccc}5& 4& 3\\ 4& 5& 4\\ 3& 4& 5\end{array}\right)$

besitzt. Das folgende Matlab-Programm realisiert den beschriebenen Algorithmus. Die Funktion randn wird zur Erzeugung standard-normalverteilter Zufallszahlen benutzt.

% Berechnung korrelierter normalverteilter Zufallszahlen
Sigma = [5 4 3; 4 5 4; 3 4 5];
mu = [−5 0 10]’; N = 10000;
L = chol(Sigma);
for i=1:N
  X(:,i) = mu + L*[randn randn randn]’;
end
x=[−14.5:0.5:14.5];
subplot(1,3,1), hist(X(1,:),x), axis([−20 20 0 1500])
subplot(1,3,2), hist(X(2,:),x), axis([−20 20 0 1500])
subplot(1,3,3), hist(X(3,:),x), axis([−20 20 0 1500])

Zahlenfolgen niedriger Diskrepanz: In einem einführenden Beispiel wurde erläutert, dass mehrdimensionale Integrale mit Hilfe der Monte-Carlo-Methode berechnet werden können. Grundidee ist dabei, zur Approximation von

$\underset{\mathrm{\Omega }}{\int }f(x)dx,\mathrm{\Omega }\subset {ℝ}^m,(\approx \mathrm{vol}(\mathrm{\Omega })E(f(x)))$

die Summe

$\frac{\mathrm{vol}(\mathrm{\Omega })}{N}\sum _{k=1}^{N}f(x_k),\mathrm{vol}(\mathrm{\Omega })-$ Volumen/Maß von $\mathrm{\Omega }$ (endlich)

mit geeigneten Zahlen $x_1,\dots ,x_N\in \mathrm{\Omega }$ zu approximieren. Sind diese Zahlen (Pseudo-)Zufallszahlen, so sprechen wir von einer stochastischen Monte-Carlo-Integration, benutzt man dagegen geschickt vorgegebene Zahlen $x_1,\dots ,x_N$, so ist der Begriff "deterministische Monte-Carlo-Integration” üblich. Die Zahlen sollten möglichst gleichmäßig verteilt sein. Die Diskrepanz einer Menge $\{x_1,\dots ,x_N\}$ definieren wir als Abweichung der Verteilung dieser Zahlen von einer angestrebten gleichmäßigen Verteilung. Damit lässt sich auch die Qualität von Pseudo-Zufallszahlen testen. Wir definieren im folgenden den Begriff der Diskrepanz und betrachten ein Beispiel.

Definition 5.3

Die Diskrepanz einer Menge $\{x_1,\dots ,x_N\}$ mit $x_k\in [0,1{]}^m$ ist definiert durch

$D_n:=\underset{Q\in [0,1{]}^m}{\sup }|\frac{\mathrm{\#}\{x_k:x_k\in Q\}}{N}-\mathrm{vol}Q|,Q-m$-dim. Quader.

Hinter der Definition der Diskrepanz steckt die Idee, dass bei einer gleichmäßig verteilten Punktmenge die relative Anzahl der Punkte, die in einem Quader $Q\subset [0,1{]}^m$ liegen, dem Volumen des Quaders entsprechen sollte, d. h. man vergleicht die Ausdrücke

$\frac{\mathrm{\#}{x}_k\in Q}{\mathrm{\#}\text{ aller Punkte}}\approx \frac{\mathrm{vol}(Q)}{\mathrm{vol}([0,1{]}^m)}.$

Definition 5.4

(1) Eine Folge $\{x_1,x_2,\dots \}$ mit $x_k\in [0,1{]}^m$ heißt gleichmäßig verteilt in $[0,1{]}^m$, wenn gilt:

$\underset{N\to \mathrm{\infty }}{\lim }{D}_N=0.$

(2) Eine Folge $\{x_1,x_2,\dots \}$ mit $x_k\in [0,1{]}^m$ heißt von niedriger Diskrepanz, wenn es eine Konstante $C_m>0$ gibt, so dass für alle genügend große $N\in \mathbb{N}$ gilt:

$D_N\le C_m\frac{(\log N{)}^m}{N}.$

Die Maßzahl der Diskrepanz ermöglicht die Angabe einer Schranke für den Fehler der Monte-Carlo-Integration. Zahlenfolgen von niedriger Diskrepanz werden auch Quasi-Zufallszahlen genannt. Das ist jedoch nicht sehr logisch, weil es sich um deterministische Folgen handelt. Wir untersuchen Beispiele und geben eine Tabelle für das Verhalten verschiedener Nullfolgen an:

$\begin{array}{ccccc}N& \frac{1}{\sqrt{N}}& \frac{\log N}{N}& \frac{(\log N{)}^2}{N}& \frac{(\log N{)}^3}{N}\\ 10^1& 0.3162277660& 0.2302585093& 0.5301898110& 1.2208071554\\ 10^2& 0.1000000000& 0.0460517019& 0.2120759244& 0.9766457243\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 10^9& 0.0000316228& 0.0000000207& 0.0000004295& 0.0000088997\end{array}$

Beispiel 5.6

(1) Die Menge $\{x_1,\dots ,x_N\}$ mit $x_k=k/N$ liefert eine Folge von niedriger Diskrepanz, weil $D_N=1/N$ ist. Allerdings muss für jedes $N$ eine neue Folge $x_k$ berechnet und im Falle der Monte-Carlo-Integration eine neue Funktionswertauswertung vorgenommen werden. Es ist praktisch viel effizienter, bereits berechnete Zahlen $x_k$ für wachsendes $N$ zu verwenden. Daher ist eine so konstruierte Folge ungeeignet. Der Wert $D_N$ lässt sich allerdings für $m=1$ nicht weiter verbessern.

(2) Seien $U_k\in \{0,\frac{1}{M},\frac{2}{M},\dots ,\frac{M-1}{M}\}$ Pseudo-Zufallszahlen, die durch Kongruenz-Generatoren erzeugt werden. Diese sind nicht gleichmäßig verteilt, denn mit $Q^{\prime }=[\frac{1}{2M+2},\frac{1}{M+1}]$ folgt (wegen $U_k\notin Q^{\prime }$)

$D_N:=|\frac{\mathrm{\#}\{U_k:U_k\in Q^{\prime }\}}{N}-\mathrm{vol}{Q}^{\prime }|=\mathrm{vol}{Q}^{\prime }=\frac{1}{2M+2}$ für $M\in \mathbb{N}$.

(3) Das Ziel, eine Punktfolge von niedriger Diskrepanz zu erzeugen, so dass mit wachsendem $N$ die Verteilung zunehmend feiner wird, wird mit der Corput-Folge erreicht:

$\{\frac{1}{2},\frac{1}{4},\frac{3}{4},\frac{1}{8},\frac{5}{8},\frac{3}{8},\frac{7}{8},\frac{1}{16},\dots \}.$

Das $k$-te Element dieser Folge erhält man einfach durch Bit-Umkehr der Dualdarstellung von $k$, d. h.

$k=(d_j\dots d_0{)}_2:=\sum _{\nu =0}^j{d}_{\nu }\cdot 2^{\nu }\mapsto x_k=(.d_0\dots d_j{)}_2:=\sum _{\nu =0}^j{d}_{\nu }\cdot 2^{-\nu -1}.$

Beispielsweise erhält man

$k=6=(110{)}_2\mapsto x_k=(.011{)}_2=\frac{0}{2^1}+\frac{1}{2^2}+\frac{1}{2^3}=\frac{3}{8}.$

Dieser Ansatz läßt sich auf eine Basis $b\ge 2$ verallgemeinern, indem definiert wird

$k=\sum _{\nu =0}^j{d}_{\nu }\cdot b^{\nu }\mapsto x_k={\varphi }_b(\nu ):=\sum _{\nu =0}^j{d}_{\nu }\cdot b^{-\nu -1}.$

Man nennt ${\varphi }_b$ die Radix-inverse Funktion. Das folgende Matlab-Programm erzeugt die ersten $N$ Corput-Zahlen zur Basis $b$.

% Berechnung der ersten N Corput-Zahlen zur Basis b
function x = corput(N,b)
  m = fix(log(N)/log(b)); % hoechste Potenz bestimmen
  D = [ ];
  n = 1:N;
  for i = 0:m % bestimme alle x(i) simultan
    d = mod(n,b);
    n = (n−d)/b;
    D = [D;d];
  end
  x = ((1/b).^(1:(m+1)))*D;

Die Idee des Algorithmus basiert auf der Formulierung

$k=(d_{\nu }{b}^{\nu -1}+\dots +d_1)\cdot b+d_0,$

d. h. die Ziffern $d_j$ werden durch Abdividieren von $b$ ermittelt.

(4) Die Radix-inverse Funktion erlaubt auch die Konstruktion von Punkten in $[0,1{]}^m$. Es seien $p_1,\dots ,p_m$ paarweise teilerfremde natürliche Zahlen. Dann heißt die Menge der Vektoren

$x_k=({\varphi }_{p_1}(k),\dots ,{\varphi }_{p_m}(k)),k\in \mathbb{N},$

Halton-Folge. Für $p_1=2$ und $p_2=3$ erzeugt man mit dem folg. Matlab-Programm.

% Berechnet zwei-dimensionale Halton-Folge
p1 = 2; p2 = 3; N = 10000;
x = [corput(N,p1); corput(N,p2)];
plot(x(1,:),x(2,:),’.’)

Beispiel 5.7

Die Halton-Folge kann dazu benutzt werden, um den Preis einer Basket-Option auf eine große Zahl von Basiswerten (vgl. Bsp. 5.2) zu berechnen. Wir betrachten das Integral

$\theta =\underset{(0,1{)}^m}{\int }{e}^{-|x{|}^2/2}dx,x=(x_1,\dots ,x_m{)}^T,$

wobei $m\gg 1$; zum Vergleich steht der exakte Wert des Integrals zur Verfügung. Mit $\varphi (x)-\frac{1}{2}=\frac{1}{2}\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)$ und $\varphi (x)=\sqrt{2\pi }\underset{-\mathrm{\infty }}{\overset{x}{\int }}\exp (-{\xi }^2/2)d\xi$ folgt nämlich

$\theta =\prod _{k=1}^m\underset{0}{\overset{1}{\int }}\exp (-\frac{x_k^2}{2})d{x}_k$, (Satz von Fubini, feste Grenzen)

$=\prod _{k=1}^m(\underset{-\mathrm{\infty }}{\overset{1}{\int }}\exp (-\frac{x_k^2}{2})d{x}_k-\underset{-\mathrm{\infty }}{\overset{0}{\int }}\exp (-\frac{x_k^2}{2})d{x}_k)$

$=\prod _{k=1}^m\sqrt{2\pi }(\varphi (1)-\varphi (0))=\prod _{k=1}^m\sqrt{2\pi }(\varphi (1)-\frac{1}{2})$

$={(\sqrt{\frac{\pi }{2}}\mathrm{erf}\left(\frac{1}{\sqrt{2}}\right))}^m.$

Als Schätzer für das Integral $\theta$ verwenden wir

${\theta }_N=\frac{1}{N}\sum _{j=1}^N\exp \left(\frac{|x^{(j)}{|}^2}{2}\right)dx,$

wobei $x^{(1)},\dots ,x^{(N)}$ $m$-dimensionale Halton-Zahlen sind. Das folgende Matlab-Programm realisiert die Konstruktion.

$\begin{array}{ccccc}N& m=10& m=20& m=50& \text{CPU-Zeit }[s]\\ 1000& 0.213699& 4.7695135\text{e-02}& 5.4333431\text{e-03}& 0.08\\ 10000& 0.210716& 4.4806662\text{e-02}& 8.9903270\text{e-04}& 1.16\\ 100000& 0.210343& 4.4285053\text{e-02}& 4.6574057\text{e-04}& 15.02\\ 200000& 0.210323& 4.4257742\text{e-02}& 4.4108892\text{e-04}& 30.33\\ 500000& 0.210306& 4.4241263\text{e-02}& 4.2502363\text{e-04}& 79.80\\ \text{exakt}& 0.21029627& 4.42245209\text{e-02}& 4.11299177\text{e-04}& \end{array}$

% Approximation des Integrals von exp(-|x|^2/2)dx ueber (0,1)^m
% mit Hilfe von Halton-Zahlen
m = 50; % Dimension des Integrals
n = 10000; % Anzahl der Punkte
p = primes(10*n); % die ersten 10*n Primzahlen
% Konstruktion der m-dimensionalen Haltonzahlen
x = [ ];
for i=1:m
  x = [x;corput(n,p(i))];
end
% Berechnung der Approximation
q = zeros(n,1)’;
for i=1:m
  q = q + x(i,:).*x(i,:);
end
int = mean(exp(−q/2))

In der obigen Tabelle sind die berechneten Werte im Vergleich mit den exakten Werten angegeben. Für größere Dimensionen müssen recht viele Halton-Punkte benutzt werden, um eine vertretbare Genauigkeit zu erhalten. Allerdings ist die Anzahl dieser Punkte wesentlich geringer als bei einer Auswertung mit Quadraturformeln.

5.3 Numerische Integration stochastischer Differentialgleichungen

Wir leiten einige Approximationen für stochastische Differentialgleichungen her und untersuchen, in welchem Sinne diese Approximationen gegen die ”exakte” Lösung der Differentialgleichung konvergieren.

Starke und schwache Konvergenz:
Für eine gewöhnliche Differentialgleichung

$\frac{dx(t)}{dt}=a(x(t),t)$ bzw. $dx=a(x(t),t)dt,t>0$,

mit Lipschitz-stetiger Funktion $x\mapsto a(x,t)$ kann man zeigen, dass das Euler-Verfahren mit der Schrittweite $h>0$

$y_{k+1}=y_k+ha(y_k,t_k),k=0,1,\dots ,$

wobei $t_k=kh$ gilt, die Konvergenzordnung 1 hat, d. h.

$\underset{k=0,\dots ,n}{\sup }|x(t_k)-y_k|\le C\cdot h.$

Hierbei ist $C>0$ eine von $h$ unabhängige Konstante. Gilt diese Aussage auch für das Euler-Maruyama-Verfahren (5.4) für stochastische Differentialgleichungen?

Wir betrachten zunächst für $0<t<T$ die skalare stochastische Differentialgleichung (Abkürzung: SDE)

(5.5) $d{x}_t=a(x_t,t)dt+b(x_t,t)d{W}_t$

für einen gegebenen Wiener-Prozess $W_t$. Das Euler-Maruyama-Verfahren für diese SDE lautet (mit der Schrittweite $h=t_{k+1}-t_k,T=nh$ und Startwert $y_0=x_0$):

(5.6) Für $k=0,\dots ,n-1$:

$t_{k+1}=t_k+h,$

$\mathrm{\Delta }{W}_k=W_{t_{k+1}}-W_{t_k},$

$y_{k+1}=y_k+a(y_k,t_k)h+b(y_k,t_k)\mathrm{\Delta }{W}_k,$

wobei die Realisierungen des Wiener-Prozesses $W_t$ dieselben sind, wie diejenigen für die SDE (5.5). Das erlaubt es, die Trajektorien $x_{t_k}$ mit $y_k$ paarweise zu vergleichen und einen punktweisen Fehler $|x_{t_k}-y_k|$ oder $|x_T-x_T^h|$ einzuführen, wobei $x_T^h:=y_n$ ist. Uns interessiert ein ”gemittelter” Fehler:

Definition 5.5

Der absolute Fehler der Differenz $x_T-x_T^h$ ist definiert durch

(5.7) $\epsilon (h)=E(|x_T-x_T^h|).$

Analog zum Fall gewöhnlicher Differentialgleichungen definieren wir die Konvergenz wie folgt:

Definition 5.6

Sei $x_t$ eine Lösung einer SDE und $x_t^h$ eine Approximation von $x_t$. Wir sagen, $x_T^h$ konvergiert stark mit der Ordnung $\gamma >0$ gegen $x_T$, falls eine Konstante $C>0$ existiert, so dass für alle (genügend kleinen) $h>0$ gilt:

(5.8) $\epsilon (h)\le C{h}^{\gamma }.$

Die Folge $\{x_T^h\}$ heißt stark konvergent gegen $x_T$, wenn

$\underset{h\to 0}{\lim }\epsilon (h)=0.$

Zur Berechnung des Erwartungswertes (5.7) bestimmen wir für die Stichprobe $X=(X_1,\dots ,X_n)$ den Wert, den der Schätzer (d. h. eine Approximation) für den Erwartungswert $E(X)$ liefert:

${\theta }_n=\frac{1}{n}\sum _{k=1}^n{X}_k.$

Der Schätzer ist erwartungstreu, d. h.

$E({\theta }_n)=\frac{1}{n}\sum _{k=1}^{n}E(X_k)=\frac{1}{n}\sum _{k=1}^{n}E(X)=E(X)$

und die Varianz konvergiert gegen Null für $n\to \mathrm{\infty }$:

$\mathrm{Var}({\theta }_n)=\frac{1}{n^2}\sum _{k=1}^n\mathrm{Var}(X_k)=\frac{1}{n}\mathrm{Var}(X)\to 0(n\to \mathrm{\infty }).$

Wir untersuchen das Euler-Maruyama-Verfahren, angewandt auf die bekannte SDE

(5.9) $d{X}_t=\mu X_{t}dt+\sigma X_{t}d{W}_t,$

empirisch auf starke Konvergenz. Nehmen wir an, die Ungleichung (5.8) gilt annähernd als Gleichung. Dann folgt

$\log \epsilon (h)\approx \log C+\gamma \log h.$

Wenn wir also $\epsilon (h)$ über $h$ mit einer doppeltlogarithmischen Skala plotten, so können wir die Konvergenzordnung als Steigung der Geradengleichung

$y(x)\approx \log C+\gamma x$ mit $y(x)=\log \epsilon (h),x=\log (h)$

bestimmen. Hierzu erzeugen wir $N$ Realisierungen $W_t^1,\dots ,W_t^N$ eines Wiener-Prozesses und lösen für jede dieser Realisierungen

• die SDE (5.5) exakt, nämlich (vgl. Kap. 3.3)

$x_t^k=x_0\exp ((\mu -\frac{1}{2}{\sigma }^2)t+\sigma W_t^k),k=1,\dots ,N$

und notieren $x_T^k$.

• die Approximation (5.6) mit Schrittweiten $h_1,\dots ,h_m$ numerisch und notieren $x_T^{k,h_1},\dots ,x_T^{k,h_m}$.

% Test auf starke Konvergenz des Euler-Maruyama-Verfahrens
randn(’state’,3)
mu = 2; sigma = 1; X0 = 1; T = 1;
N = 1000; % Anzahl der Pfade der Brownschen Bewegung
m = 6;    % Anzahl der verschiedenen Schrittweiten
K = 2^9;  % Anzahl der Gitterpunkte, wenn T = 1
h = T/K;  % kleinste Schrittweite
for s=1:N
  dW = sqrt(h)*randn(1,K);
  W = sum(dW);
  Xexakt = X0*exp((mu−sigma^2/2)+sigma*W(end));
  for p=1:m
    R = 2^(p−1);
    dt = R*h; % aktuelle Schrittweite
    L = K/R;  % Anzahl der Euler-Schritte
    X = X0;
    for j=1:L
      Wink = sum(dW(R*(j−1)+1:R*j));
      X = X + mu*X*dt + sigma*X*Wink;
    end
  Xerr(s,p) = abs(X − Xexakt);
  end
end
% Plotten der Fehler und einer Geraden mit Steigung 1/2
dtlist = h*(2.^(0:m−1));
loglog(dtlist,mean(Xerr),’*-’), hold on
loglog(dtlist,(dtlist.^(0.5)),’--’), hold off
% Kleinste-Quadrate-Methode Ax = b mit x = (logC gamma)
A = [ones(m,1), log(dtlist)’];
b = log(mean(Xerr)’);
x = A\b;
gamma = x(2), residuum = norm(A*x−b)

Der absolute Fehler $\epsilon (h)$ wird durch den Schätzer

(5.10) $\hat{\epsilon }(h_j)=\frac{1}{N}\sum _{k=1}^N|x_T^k-x_T^{k,h_j}|$

approximiert.

Das angegebene Matlab-Programm zum Euler-Maruyama-Verfahren realisiert diese Vorgehensweise zur Approximation von (5.9) für $\mu =2,\sigma =1$ und $T=1$ und erzeugt eine Abbildung, die den Fehlerschätzer $\hat{\epsilon }(h)$ für verschiedene Schrittweiten $h$ darstellt. Die eingezeichnete Vergleichsgerade legt $\gamma =0.5$ als Konververgenzordnung nahe. Tatsächlich ergibt eine lineare Ausgleichsrechnung $\gamma =0.536$ bei einem Residuum von 0.026.

In vielen Anwendungen ist man nicht an den Trajektorien $x_T$ selbst, sondern nur an Momenten von $x_T$ (etwa dem Erwartungswert oder der Varianz) interessiert. Demzufolge suchen wir nur Approximationen von $E(x_T)$ bzw. von $\mathrm{Var}(x_T)$, nämlich $E(x_T^h)$ bzw. $\mathrm{Var}(x_T^h)$. Das führt auf den folgenden abgeschwächten Konvergenzbegriff.

Definition 5.7

Sei $x_t$ eine Lösung einer SDE und $x_t^h$ eine Approximation von $x_t$. Wir nennen $x_t^h$ schwach konvergent bezüglich $g$ mit der Ordnung $\gamma >0$, wenn eine Konstante $C>0$ existiert, so dass für alle (genügend kleinen) $h>0$ gilt:

$|E(g(x_t))-E(g(x_t^h))|\le C{h}^{\gamma }.$

Im Falle $g=id$ (id - identischer Operator) nennen wir $x_t^h$ schwach konvergent mit der Ordnung $\gamma$.

Da wir nicht an einer pfadweisen Konvergenz interessiert sind, können wir auch verschiedene Pfade für jeden Zeitschritt $y_k\mapsto y_{k+1}$ im Algorithmus (5.6) verwenden.

Wir untersuchen die schwache Konvergenzordnung des Euler-Maruyama-Verfahrens. Es wird ein Test ähnlich dem obigen für die SDE (5.9) mit $\mu =2,\sigma =0.1$ und $T=1$ in einem Matlab-Programm realisiert.

% Test auf schwache Konvergenz des Euler-Maruyama-Verfahrens
randn(’state’,100)
mu = 2; sigma = 0.1; X0 = 1; T = 1;
N = 50000; % Anzahl der Pfade der Brownschen Bewegung
m = 5;     % Anzahl der verschiendenen Schrittweiten
for p=1:m
  h = 2^(p−10); % aktuelle Schrittweite
  L = T/h;      % Anzahl der Euler-Schritte
  X = X0*ones(N,1);
  for j=1:L
    dW = sqrt(h)*randn(N,1);
    X = X + mu*X*h + sigma*X.*dW;
  end
  Xerr(p) = abs(mean(X) − exp(mu));
end
% Plotten der Fehler und einer Geraden mit Steigung 1
dtlist = 2.^([1:m]−10);
loglog(dtlist,Xerr,’*-’), hold on
loglog(dtlist,dtlist,’--’), hold off
% Kleinste-Quadrate-Methode Ax = b mit x = (logC gamma)
A = [ones(m,1), log(dtlist)’];
b = log(Xerr)’;
x = A\b;
gamma = x(2), residuum = norm(A*x−b)

Man beachte, dass der Erwartungswert der exakten Lösung $x_t$ von (5.9) gleich $E(x_t)=\exp (\mu t)$ ist.

Diesmal vermuten wir eine Konvergenzordnung von 1. Eine lineare Ausgleichsrechnung ergibt $\gamma =0.9858$ bei einem Residuum von 0.0508.

Wir wollen nun Verfahren höherer Konvergenzordnung entwickeln. Dazu greifen wir auf eine spezielle Klasse von Integrationsverfahren zurück, die auch für gewöhnliche Differentialgleichungen verwendet werden, nämlich auf Taylorreihen-Verfahren.

Stochastische Taylorentwicklungen:
Zuerst betrachten wir das gewöhnliche (autonome) Anfangswertproblem

$x^{\prime }(t)=a(x(t)),x\in {ℝ}^n,t>t_0,x(t_0)=0.$

Um Verfahren höherer Ordnung als das Euler-Verfahren

$y_{k+1}=y_k+a(y_k)h,k=0,1,2,\dots ,y_0=x_0$

für die Approximation $y_k$ von $x(t_k)$ (mit $t_k=kh,h>0$) herzuleiten, entwickeln wir die Lösung $x(t)$ in eine Taylorreihe um $t$, wobei genügend hohe Regularität der Lösung vorausgesetzt wird:

$x(t+h)=x(t)+h{x}^{\prime }(t)+\frac{h^2}{2}{x}^{″}(t)+\frac{h^3}{6}{x}^{‴}(t)+O(h^4).$

Rekursives Einsetzen der rechten Seite der DGl. liefert die Entwicklung

$x(t+h)=x(t)+ha(x(t))+\frac{h^2}{2}(a(x(t)){)}^{\prime }+\frac{h^3}{6}(a(x(t)){)}^{″}+O(h^4)$

$=x(t)+ha(x(t))+\frac{h^2}{2}{a}^{\prime }(x(t))\cdot a(x(t))+\frac{h^3}{6}[a^{″}(x(t))(a(x(t)){)}^2+(a^{\prime }(x(t)){)}^2\cdot a(x(t))]+O(h^4),$

wobei $a(x(t){)}^2$ das Argument des Tensors $a^{″}(x(t))$ ist, d. h.

$a^{″}(x(t))\cdot (a(x(t)){)}^2=a^{″}(x(t))(a(x(t)),a(x(t)))\in {ℝ}^n.$

Wählen wir $t=t_k$ und $y_k\approx x(t_k)$ und vernachlässigen wir den Restterm $O(h^4)$, so erhalten wir das Taylor-Einschrittverfahren:

Für $k=0,\dots ,n-1$:

$a_k:=a(y_k),a{\text{'}}_k:=a^{\prime }(y_k),a^{\prime }{\text{'}}_k:=a^{″}(y_k),$

$y_{k+1}=y_k+h{a}_k+\frac{h^2}{2}a{\text{'}}_k\cdot a_k^2+\frac{h^3}{6}[a^{\prime }{\text{'}}_k\cdot (a_k{)}^2+(a{\text{'}}_k{)}^2\cdot a_k]$

Wir können diese Idee auf SDE übertragen, indem wir die Taylorentwicklung durch eine stochastische Version ersetzen. Diese liefert gerade das Lemma von Itô.

Wir betrachten zunächst zur Vereinfachung der Notation die skalare, eindimensionaleund autonome SDE

(5.11) $d{x}_t=a(x_t)dt+b(x_t)d{W}_t.$

Das Lemma von Itô für $f(x_t)$ lautet in Integralform für die nicht explizit von $t$ abhängige Funktion $f(x_t)$

(5.12) $f(x_t)=f(x_{t_0})+\underset{t_0}{\overset{t}{\int }}(f^{\prime }(x_s)a(x_s)+\frac{1}{2}{f}^{″}(x_s)(b(x_s){)}^2)ds+\underset{t_0}{\overset{t}{\int }}b(x_s)f^{\prime }(x_s)d{W}_s.$

Speziell für $f(x)=x$ folgt

(5.13) $x_t=x_{t_0}+\underset{t_0}{\overset{t}{\int }}a(x_s)ds+\underset{t_0}{\overset{t}{\int }}b(x_s)d{W}_s.$

Wir setzen (5.12) für $f=a$ und $f=b$ in (5.13) ein:

(5.14) $x_t=x_{t_0}+\underset{t_0}{\overset{t}{\int }}(a(x_{t_0})+\underset{t_0}{\overset{s}{\int }}(a^{\prime }a+\frac{1}{2}{a}^{″}{b}^2)dz+\underset{t_0}{\overset{s}{\int }}{a}^{\prime }bd{W}_z)ds$

$+\underset{t_0}{\overset{t}{\int }}(b(x_{t_0})+\underset{t_0}{\overset{s}{\int }}(b^{\prime }a+\frac{1}{2}{b}^{″}{b}^2)dz+\underset{t_0}{\overset{s}{\int }}{b}^{\prime }bd{W}_z)d{W}_s$

Dabei wurden die Abkürzungen $a=a(x_z),b=b(x_z)$ usw. benutzt.

Fassen wir die Doppelintegrale zu einem Restterm $R$ zusammen, so erhalten wir

$x_t=x_{t_0}+a(x_{t_0})(t-t_0)+b(x_{t_0})\underset{t_0}{\overset{t}{\int }}d{W}_s+R.$

Vernachlässigen von $R$ liefert eine (recht umständliche) Herleitung des Euler-Maruyama-Verfahrens. Wir erhalten ein Verfahren höherer Ordnung, indem wir das Doppelintegral bezüglich $d{W}_{z}d{W}_s$ aus dem Restterm herausnehmen und den Integranden $b^{\prime }(x_z)\cdot b(x_z)$ durch $b^{\prime }(x_{t_0})\cdot b(x_{t_0})$ ersetzen:

$x_t=x_{t_0}+a(x_{t_0})(t-t_0)+b(x_{t_0})\underset{t_0}{\overset{t}{\int }}d{W}_s+b^{\prime }(x_{t_0})b(x_{t_0})\underset{t_0}{\overset{t}{\int }}\underset{t_0}{\overset{s}{\int }}d{W}_{z}d{W}_s+\tilde{R}.$

Die zugrunde liegende Idee ist, dass sich $\underset{t_0}{\overset{t}{\int }}f(x_s)d{W}_s$ durch $O(\sqrt{h})$ abschätzen lässt (motiviert durch die Merkregel $d{W}_t=\sqrt{dt}$; siehe Abschnitt 4), also sind alle anderen Terme in $\tilde{R}$ von höherer Ordnung als der zusätzliche Term mit $d{W}_z\cdot d{W}_s$. Für die Fehlerterme gilt damit: $R=O(h)$ und $\tilde{R}=O(h^{3/2})$. Wir berechnen die Doppelintegrale (analog wie in Kap. 4):

(5.15) $\underset{t_0}{\overset{t}{\int }}\underset{t_0}{\overset{s}{\int }}d{W}_{z}d{W}_s=\underset{t_0}{\overset{t}{\int }}(W_s-W_{t_0})d{W}_s=\frac{1}{2}(W_t^2-W_{t_0}^2)-\frac{t-t_0}{2}-W_{t_0}(W_t-W_{t_0})$

$=\frac{1}{2}((W_t-W_{t_0}{)}^2-(t-t_0)).$

Das führt auf das Milstein-Verfahren für die nicht-autonome SDE (5.5):

Für $k=0,\dots ,n-1$:

$t_{k+1}:=t_k+h,$

$\mathrm{\Delta }W:=Z\sqrt{h}$ mit $Z\sim 𝒩(0,1)$,

$y_{k+1}:=y_k+a(y_k,t_k)h+b(y_k,t_k)\mathrm{\Delta }W+\frac{1}{2}{b}^{\prime }(y_k,t_k)b(y_k,t_k)((\mathrm{\Delta }W{)}^2-h)$, wobei $b^{\prime }=\frac{\partial b}{\partial x}$.

Um das Milstein-Verfahren für die SDE (5.9) mit $\mu =2,\sigma =1$ und $T=1$ in Matlab zu implementieren, genügt es, den Term $\frac{1}{2}{\sigma }^{2}X((\mathrm{\Delta }W{)}^2-h)$ zu der Euler-Maruyama-Approximation von $X$ hinzuzufügen. Eine lineare Ausgleichsrechnung liefert $\gamma =0.971$ bei einem Residuum von 0.048. Tatsächlich beträgt die starke Konvergenzordnung des Milstein-Verfahrens 1 (also um 1/2 besser als das Euler-Maruyama-Verfahren). Für einen Konvergenzbeweis verweisen wir auf [10].

Stochastische Runge-Kutta-Verfahren:
Als Beispiel für ein derartiges Verfahren leiten wir einen Algorithmus vom Milstein-Typ her. Um die Ableitung b'(x_t) zu ersetzen, entwickeln wir formal:

$b(x_t+\mathrm{\Delta }{x}_t)-b(x_t)=b^{\prime }(x_t)\mathrm{\Delta }{x}_t+𝒪(|\mathrm{\Delta }{x}_t{|}^2)=b^{\prime }(x_t)(a(x_t)h+b(x_t)\mathrm{\Delta }{W}_t)+𝒪(h)=b^{\prime }(x_t)b(x_t)\mathrm{\Delta }{W}_t+𝒪(h).$

Ersetzen wir $\mathrm{\Delta }{W}_t$ durch den Mittelwert $\sqrt{h}$ (wieder wie in Kapitel 4), so folgt

$b^{\prime }(x_t)b(x_t)=\frac{1}{\sqrt{h}}(b(x_t+\mathrm{\Delta }{x}_t)-b(x_t))+𝒪(\sqrt{h})=\frac{1}{\sqrt{h}}(b(x_t+a(x_t)h+b(x_t)\mathrm{\Delta }{W}_t)-b(x_t))+𝒪(\sqrt{h}).$

Damit erhalten wir die stochastische Runge-Kutta-Variante des Milstein-Verfahrens:

(5.16) Für $k=0,\dots ,n-1$:

$t_{k+1}:=t_k+h,$

$\mathrm{\Delta }W:=Z\sqrt{h}$ mit $Z\sim 𝒩(0,1)$,

$\hat{y}:=y_k+a(y_k,t_k)h+b(y_k,t_k)\sqrt{h},$

$y_{k+1}:=y_k+a(y_k,t_k)h+b(yk,tk)\mathrm{\Delta }W+\frac{1}{2\sqrt{h}}(b(\hat{y},t_k)-b(y_k,t_k))((\mathrm{\Delta }W{)}^2-h).$

Bemerkung:

Wie kann eine allgemeine Klasse von stochastischen Runge-Kutta-Verfahren aussehen? Ein Versuch wäre die Definition

$y_{k+1}=y_k+h\sum _{j=1}^s{d}_{j}a({\hat{y}}_j)+\mathrm{\Delta }W\sum _{j=1}^s{e}_{j}b({\hat{y}}_j)$

mit den Zuwächsen

${\hat{y}}_j=y_j+h\sum _{\nu =1}^s{D}_{j\nu }a(y_{\nu })+\mathrm{\Delta }W\sum _{\nu =1}^s{E}_{j\nu }b(y_{\nu }),j=1,\dots ,s.$

Leider gibt es hierfür eine Schranke für die Konvergenzordnung. Derartige Verfahren können höchstens eine (starke) Konvergenzordnung von 1 haben und sind folglich nicht besser als das Milstein-Verfahren oder die oben konstruierte Runge-Kutta-Methode. Umgehen kann man diese Schranke nur, wenn weitere Zufallsvariable benutzt werden, um die Integrale der stochastischen Taylor-Entwicklung zu approximieren. Dann ersetzt man die Ausdrücke $e_j\mathrm{\Delta }W$ bzw. $E_{j\nu }\mathrm{\Delta }W$ durch Zufallsvariable $Z_j$ bzw. $Z_{jk}$. Details findet man in der Literatur [1].

Systeme von SDE:
Wir können das Milstein-Verfahren auf Systeme von SDE der Dimension $n$ mit $m$ Wiener-Prozessen $W_t=(W_t^1,\dots ,W_t^m)$ verallgemeinern:

$d{x}_t=a(x_t,t)dt+b(x_t,t)d{W}_t,$

wobei gilt

$a:{ℝ}^n\times ℝ\to {ℝ}^n,b=(b^{j\nu }):{ℝ}^n\times ℝ\to {ℝ}^{n\times m}.$

Solche Fälle treten etwa bei der Modellierung von Optionen mit stochastischer Volatilität (Bsp. 5.1) oder bei der Modellierung von Basket-Optionen (Bsp. 5.2) auf. Ist $(W_t^1,\dots ,W_t^n)$ eine $n$-dimensionale Bewegung, so lautet die Gleichung, die die Dynamik der Aktienkurse $S_t^1,\dots ,S_t^n$ beschreibt,

$d{S}_t^j={\mu }_t^j{S}_t^{j}dt+\sum _{k=1}^n{\sigma }_t^{jk}d{W}_t^k,j=1,\dots ,n.$

Der Fall eines einzigen Wiener-Prozesses $m=1$ ist einfach. Der Term $b^{\prime }b$ geht für $b=(b^1,\dots ,b^n{)}^T$ über in

$Db(x_t,t)b(x_t,t)=\left(\begin{array}{ccc}\frac{\partial b^1}{\partial x_1}& \dots & \frac{\partial b^1}{\partial x_n}\\ ⋮& \ddots & ⋮\\ \frac{\partial b^n}{\partial x_1}& \dots & \frac{\partial b^n}{\partial x_n}\end{array}\right)\left(\begin{array}{c}{b}^1\\ ⋮\\ b^n\end{array}\right)$

und das Milstein-Verfahren schreibt sich als

$y_{k+1}=y_k+a(y_k,t_k)h+b(y_k,t_k)\mathrm{\Delta }W+\frac{1}{2}((\mathrm{\Delta }W{)}^2-h)Db(y_k,t_k)b(y_k,t_k).$

Ähnlich lässt sich das obige stochastische Runge-Kutta-Verfahren umformulieren.

Der allgemeine Fall von $m>1$ Wiener-Prozessen ist komplizierter. Wiederholen wir die Herleitung des Milstein-Verfahrens für skalare SDE, so erhalten wir das Milstein-Verfahren für Systeme:

(5.17) $y_{k+1}^s=y_k^s+a^s(y_k,t_k)h+\sum _{j=1}^m{b}^{sj}(y_k,t_k)\mathrm{\Delta }{W}^j+\sum _{i,j=1}^m\sum _{l=1}^n{b}^{lj}\frac{\partial b^{si}}{\partial x_l}(y_k,t_k)\underset{t_k}{\overset{t_{k+1}}{\int }}\underset{t_k}{\overset{\tau }{\int }}d{W}_z^{j}d{W}_{\tau }^i,s=1,\dots ,n.$

Wir sind hier mit zwei Problemen konfrontiert:

• Die stochastischen Integrale

$I_{ji}:=\underset{t_k}{\overset{t_{k+1}}{\int }}\underset{t_k}{\overset{\tau }{\int }}d{W}_z^{j}d{W}_{\tau }^i$

sind nicht mehr einfach auf elementare Integrale der Form

(5.18) $\mathrm{\Delta }{W}^{\nu }=\underset{t_k}{\overset{t_{k+1}}{\int }}d{W}_{\tau }^{\nu }$

zurückführbar. Die Frage ist, wie sie sich berechnen lassen.

• Die Differentialoperatoren

$L_j:=\sum _{l=1}^n{b}^{lj}\frac{\partial }{\partial x_l}$

müssen auf alle Spaltenvektoren $b^{\nu }=(b^{s\nu }{)}_{s=1}^n$ angewendet werden. Dies ist sehr mühsam, wie kann es vermieden werden?

Eine Möglichkeit, das erste Problem zu lösen, ist, die Integrale $I_{ji}$ bis auf einen Fehler $𝒪(h)$ zu approximieren, denn damit bleibt die Konvergenzordnung 1 des Milstein-Verfahrens erhalten. Interessanterweise sind die Integrale Lösungen eines (einfachen) Systems von SDE. Approximieren wir die Lösungen, so erhalten wir auch Approximationen der Integrale $I_{ji}$. Wir zeigen, wie das Integral $I_{21}$ approximiert werden kann; der allgemeine Fall funktioniert analog. Die Behauptung ist, dass $I_{21}$ die erste Komponente der Lösung des Systems

$d{x}_t=\left(\begin{array}{cc}{x}_t^2& 0\\ 0& 1\end{array}\right)d{W}_t,t_k\le t\le t_{k+1},x_t=\left(\begin{array}{c}0\\ 0\end{array}\right)$

an der Stelle $t=t_{k+1}$ ist. Dies beweisen wir in dem folgenden Lemma.

Lemma 5.1

Die Lösung $x_t$ von (5.19) lautet an der Stelle $t=t_{k+1}$:

$x_{t_{k+1}}=\left(\begin{array}{c}{I}_{21}\\ \mathrm{\Delta }{W}^2\end{array}\right)$, wobei gilt $\mathrm{\Delta }{W}^2=W_{t_{k+1}}^2-W_{t_k}^2$.

Beweis:

Aus der zweiten Gleichung $d{x}_t=d{W}_t^2$ ergibt sich mit dem Anfangswert $x_{t_k}^2=0$ durch Integration

$x_{t_{k+1}}^2=x_{t_k}^2+\underset{t_k}{\overset{t_{k+1}}{\int }}d{W}_s^2=\underset{t_k}{\overset{t_{k+1}}{\int }}d{W}_s^2=W_{t_{k+1}}^2-W_{t_k}^2.$

Für die erste Komponente folgt

$x_{t_{k+1}}^1=x_{t_k}^1+\underset{t_k}{\overset{t_{k+1}}{\int }}{x}_s^{2}d{W}_s^2=\underset{t_k}{\overset{t_{k+1}}{\int }}{x}_s^{2}d{W}_s^2=\underset{t_k}{\overset{t_{k+1}}{\int }}\underset{t_k}{\overset{s}{\int }}d{W}_{\tau }^{2}d{W}_s^1.$

q.e.d.

Wir zerlegen das Intervall $[t_k,t_{k+1}]$ in $N$ Teilintervalle der Länge $\delta =(t_{k+1}-t_k)/N$. Sei $z_{\nu }$ die Approximation von $x_{t_{k+\nu \delta }}$. Dann ist $z_0=(0,0{)}^T$. Die Euler-Maruyama-Approximation von (5.19) lautet dann:

(5.20) Für $k=0,\dots ,n-1$:

$\mathrm{\Delta }{W}_k:=W_{t_k+(\nu +1)\delta }-W_{t_k+\nu \delta }$

$z_{\nu +1}:=z_{\nu }+\left(\begin{array}{cc}{z}_k^2& 0\\ 0& 1\end{array}\right)\mathrm{\Delta }{W}_{\nu }.$

Diese Methode liefert tatsächlich eine Approximation von $I_{21}$ der Ordnung 1, wenn $N=1/h$ und $h=t_{k+1}-t_k$ gilt, denn das Euler-Maruyama-Verfahren hat die Konvergenzordnung 1/2, und damit gilt wegen $\delta =h^2$:

$E(|z_N^1-I_{21}|)\le C{\delta }^{1/2}=Ch.$

Wir kommen nun auf das zweite Problem zurück. Die Differentiation der Spaltenvektoren $b^{\nu }=(b_s^{s\nu }{)}_s$ kann vermieden werden, indem wir ein stochastisches Runge-Kutta-Verfahren – wie im skalaren Fall beschrieben – verwenden.

Lemma 5.2

Es gilt die Approximation:

(5.21) $L_j{b}^{\nu }(y_k,t_k)=\sum _{l=1}^n{b}^{lj}(y_k,t_k)\frac{\partial b^{\nu }}{\partial x_l}(y_k,t_k)=\frac{1}{\sqrt{h}}(b^{\nu }({\hat{y}}^j,t_k)-b^{\nu }(y_k,t_k))+𝒪(\sqrt{h}),$

wobei gilt

(5.22) ${\hat{y}}^j=y_k+a(y_k,t_k)h+b^j(y_k,t_k)\sqrt{h}.$

Beweis:

Für die rechte Seite von (5.21) liefert eine Taylorentwicklung:

$\frac{1}{\sqrt{h}}(b^{\nu }({\hat{y}}^j,t_k)-b^{\nu }(y_k,t_k))=\frac{1}{\sqrt{h}}(D{b}^{\nu }\cdot ({\hat{y}}^j-y_k)+𝒪(h))$

$=\frac{1}{\sqrt{h}}(D{b}^{\nu }\cdot (ah+b^j\sqrt{h})+𝒪(h))=D{b}^{\nu }\cdot b^j+𝒪(\sqrt{h}),$

wobei wir das Argument $(y_k,t_k)$ weggelassen haben. Die $s$-te Komponente von (5.24) ist gleich

$(D{b}^{\nu }\cdot b^j{)}_s=\sum _{l=1}^n\frac{\partial b^{s\nu }}{\partial x_l}{b}^{lj}=(L_j{b}^{\nu }{)}_s$

und das ist gerade die linke Seite von (5.21).

q.e.d.

Wir erhalten zusammengefasst das stochastische Runge-Kutta-Verfahren erster Ordnung für Systeme:

$y_{k+1}:=y_k+a(y_k,t_k)h+b(y_k,t_k)\mathrm{\Delta }W+\frac{1}{\sqrt{h}}\sum _{j,\nu =1}^m(b^{\nu }({\hat{y}}^{\nu },t_k)-b^{\nu }(y_k,t_k))I_{j\nu }.$

Dies ist die Verallgemeinerung des stochastischen Runge-Kutta-Verfahrens (5.17) für den skalaren Fall. Die Zwischenwerte ${\hat{y}}_j$ sind in (5.22) definiert, $\mathrm{\Delta }W$ ist durch (5.18) gegeben und die Integrale $I_{j\nu }$ können mittels (5.21) approximiert werden.

Bemerkung:

In speziellen Fällen können die Integrale $I_{j\nu }$ explizit berechnet werden:

1. Für $j=\nu$ erhalten wir wegen (5.16)

$I_{jj}=\frac{1}{2}((\mathrm{\Delta }{W}_k^j{)}^2-h)$

mit $\mathrm{\Delta }{W}_k^j=W_{t_{k+1}}^j-W_{t_k}^j$ und $h=t_{k+1}-t_k$.

2. Hängt $b$ nicht von $x$ ab (man spricht auch von additivem Rauschen), verschwindet die Ableitung $\partial b/\partial x_l$ und das Milstein- und das Runge-Kutta-Verfahren reduzieren sich zum Euler-Maruyama-Verfahren.

3. Falls $b=\mathrm{diag}(b_1(x_1,t),\dots ,b_n(x_n,t))$ (sog. diagonales Rauschen), folgt

$L_j{b}^{\nu }=\{\begin{array}{cc}0,& j\ne \nu \\ b_j\frac{\partial b_j}{\partial x_j},& j=\nu \end{array}$

und es sind nur Auswertungen von $I_{jj}$ erforderlch.

Weitere Diskretisierungen von SDE und Konvergenzbeweise sind in [10] zu finden; für Matlab-Routinen siehe [6], [7].

5.4 Varianzreduktion

In Abschnitt 5.1 haben wir gesehen, dass sehr viele Monte-Carlo-Simulationen notwendig sind, um einen halbwegs akkuraten Preis einer Option zu erhalten. In Abbildung 5.12 schwankt der Preis einer europäischen Put-Option zwischen 16.7 und 17.3 (bei einem exakten Wert von 16.98), selbst bei mehr als 10 000 Simulationen. In diesem Abschnitt geben wir einen Grund für dieses langsame Konvergenzverhalten an und stellen zwei Methoden vor, mit denen die Genauigkeit ohne großen Aufwand verbessert werden kann.

Sei ${\theta }_n$ die Monte-Carlo-Approximation eines exakten Wertes $\theta$. Beispiele für ${\theta }_n$ und $\theta$ sind

• die Lösung einer stochastischen Differentialgleichung (SDE):

$\theta =x_T$ Lösung einer SDE $dx=a(x,t)dt+b(x,t)dW$,

${\theta }_N=y_N$ Approximation von $x_T$ bei $t_N=T$,

• die stochastische Integration:

$\theta =\underset{\mathrm{\Omega }}{\int }g(x)dx=\underset{\mathrm{\Omega }}{\int }\frac{g(x)}{f(x)}f(x)dx=E(\varphi (x)),$

${\theta }_n=\frac{1}{n}\sum _{k=1}^n\varphi (X_k),$

wobei $\varphi (x)=g(x)/f(x)$ und $X_k$ sind Stichproben einer nach $F$ verteilten Zufallsvariablen mit $F^{\prime }=f$ sind.

Der Fehler der Monte-Carlo-Methode $|{\theta }_n-\theta |$ ist selbst eine Zufallsvariable und so können wir nur Fehlergrenzen für gewisse Sicherheitswahrscheinlichkeiten angeben. Wir illustrieren dies für den Fall, dass ${\theta }_n$ eine stochastische Approximation eines Integrals mit dem Wert $\theta$ darstellt (siehe oben). Wir nehmen an, dass $E(\varphi (X_k))=\theta$ und $\mathrm{Var}(\varphi (X_k))={\sigma }^2$ für alle $k=1,\dots ,n$ gilt. Dann folgt

$\mathrm{Var}({\theta }_n)=\frac{1}{n^2}\sum _{k=1}^n\mathrm{Var}(\varphi (X_k))=\frac{{\sigma }^2}{n},$

$E({\theta }_n)=\frac{1}{n}\sum _{k=1}^{n}E(\varphi (X_k))=\theta .$

Wir benutzen nun die Chebychev-Ungleichung für beliebige quadratisch integrierbare Zufallsvariable $Y$

$P(|Y-E(Y)|\ge \delta )\le \frac{\mathrm{Var}(Y)}{{\delta }^2},\delta >0$

für $\delta =\sigma /\sqrt{\epsilon n}$, um die elementare Fehlerabschätzung

$P(|{\theta }_n-\theta |\ge \frac{\sigma }{l\sqrt{\epsilon n}})\le \epsilon$

oder äquivalent

$P(|{\theta }_n-\theta |<\frac{\sigma }{l\sqrt{\epsilon n}})\ge \epsilon$

zu erhalten. Diese Ungleichung bedeutet, dass der Fehler $|{\theta }_n-\theta |$ umso kleiner wird, je größer die Stichprobenzahl $n$ ist. Allerdings muss zur Reduktion des Fehlers (oder der Standardabweichung $\sqrt{\mathrm{Var}({\theta }_n)}$) um eine Dezimalstelle (also um den Faktor 10) die Stichprobenzahl um den Faktor 100 erhöht werden! Dies erklärt die langsame Konvergenz der Monte-Carlo-Simulationen. Eine andere Idee, den Fehler zu verkleinern, ist es, die Varianz $\mathrm{Var}({\theta }_n)$ möglichst klein zu halten. Diese Möglichkeit der Konvergenzverbesserung bezeichnet man als Technik der Varianzreduktion. Wir skizzieren zwei Techniken:

• die Methode der Abtrennung des Hauptteils und
• die Methode der antithetischen Variablen.

Die Methode der Abtrennung des Hauptteils versucht, durch geschicktes Hinzuaddieren eines zweiten Integranden die Gesamtvarianz des Schätzers zu verkleinern. Wir nehmen an, dass das Integral ${\theta }^{*}=\underset{\mathrm{\Omega }}{\int }\psi (x)f(x)dx$ für eine Funktion $\psi$ (Hauptteil genannt) analytisch berechenbar ist. Die Formulierung

$\theta =\underset{\mathrm{\Omega }}{\int }(\varphi (x)-\psi (x))f(x)dx+\underset{\mathrm{\Omega }}{\int }\psi (x)f(x)dx$

motiviert dann den neuen Schätzer

${\hat{\theta }}_n=\frac{1}{n}\sum _{k=1}^n(\varphi (X_k)-\psi (X_k))+\underset{\mathrm{\Omega }}{\int }\psi (x)f(x)dx={\theta }_n-{\theta }_n^{*}+{\theta }^{*}$, wobei ${\theta }_n^{*}=\frac{1}{n}\sum _{k=1}^n\psi (X_k)$.

Der Hauptteil $\psi$ sollte möglichst ”einfach” sein, so dass das Integral ${\theta }^{*}$ analytisch bestimmt werden kann, aber zugleich der Funktion $\varphi$ möglichst ”ähnlich” sein, damit die Varianz von ${\hat{\theta }}_n$ kleiner wird als die Varianz von ${\theta }_n$. Warum sollte das funktionieren? Wenn $\varphi$ und $\psi$ sehr ”ähnlich” sind, erwarten wir, dass sowohl $\theta$ und ${\theta }^{*}$ als auch die Approximationen ${\theta }_n$ und ${\theta }_n^{*}$ ”ähnlich” sind. Dann sollte auch die Korrelation zwischen ${\theta }_n$ und ${\theta }_n^{*}$ groß sein und nahe der oberen Schranke liegen, also:

(5.24) $\mathrm{Cov}({\theta }_n,{\theta }_n^{*})\approx \frac{1}{2}(\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{*}))>\frac{1}{2}\mathrm{Var}({\theta }_n^{*}).$

Die obere Schranke lautet tatsächlich $\frac{1}{2}(\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{*}))$, denn aus der Beziehung

(5.25) $\mathrm{Var}(X\pm Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)\pm 2\mathrm{Cov}(X,Y)$

für die zwei Zufallsvariablen $X$ und $Y$ folgt wegen $\mathrm{Var}(X\pm Y)\ge 0$ die Ungleichung

(5.26) $\mathrm{Cov}(X,Y)\le \frac{1}{2}(\mathrm{Var}(X)+\mathrm{Var}(Y)).$

Dann folgt für die neue Zufallsvariable:

$\mathrm{Var}({\hat{\theta }}_n)=\mathrm{Var}({\theta }_n-{\theta }_n^{*})$ (denn ${\theta }^{*}$ ist konstant)

$=\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{*})-2\mathrm{Cov}({\theta }_n,{\theta }_n^{*})$ (wegen (5.25))

$<\mathrm{Var}({\theta }_n)$ (wegen (5.24))

und die Varianz ist tatsächlich verringert worden.

Die zweite Methode führt eine sogenannte antithetische Variable ein. Sei ${\theta }_n$ mittels einer Zufallsvariablen $Z\sim 𝒩(0,1)$ erzeugt worden. Definiert man dann eine Approximation ${\theta }_n^{-}$, die genauso wie ${\theta }_n$ erzeugt wurde, aber mit $-Z\sim 𝒩(0,1)$ ist. Die antithetische Variable lautet

${\hat{\theta }}_n=\frac{1}{2}({\theta }_n+{\theta }_n^{-}).$

Wir behaupten, dass die Varianz von ${\hat{\theta }}_n$ kleiner ist als die von ${\theta }_n$ (zumindest wenn $\mathrm{Var}({\theta }_n)\approx \mathrm{Var}({\theta }_n^{*})$. Aus (5.25) folgt

$\mathrm{Var}{\hat{\theta }}_n=\frac{1}{4}\mathrm{Var}({\theta }_n+{\theta }_n^{-})=\frac{1}{4}(\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{-})+2\mathrm{Cov}({\theta }_n,{\theta }_n^{*}))$

und unter Berücksichtigung von (5.26) erhalten wir

$\mathrm{Var}{\hat{\theta }}_n<\{\begin{array}{cc}\frac{1}{4}(\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{-})),& \text{wenn }\mathrm{Cov}({\theta }_n,{\theta }_n^{-})\le 0\\ \frac{1}{2}(\mathrm{Var}({\theta }_n)+\mathrm{Var}({\theta }_n^{-})),& \text{wenn }\mathrm{Cov}({\theta }_n,{\theta }_n^{-})>0\end{array}$

Im Falle $\mathrm{Var}({\theta }_n)\approx \mathrm{Var}({\theta }_n^{-})$ erhalten wir also $\mathrm{Var}({\hat{\theta }}_n)<\mathrm{Var}({\theta }_n)$.

Wir illustrieren diese Methode mit einer Monte-Carlo-Simulation der europäischen Put-Option aus Beispiel 5.1. Die Variablen ${\theta }_n$ bzw. ${\theta }_n^{-}$ seien die Auszahlungsfunktionen zu den Aktienkursen $S_k$ bzw. $S_k^{-}$, die durch

$S_{k+1}=S_k(1+rh+Z\sqrt{h}),$

$S_{k+1}^{-}=S_k^{-}(1+rh+Z\sqrt{h}),k=1,\dots ,N-1$

mit $Z\sim 𝒩(0,1)$ definiert sind. Die Auszahlungsfunktion ist dann gegeben durch $\frac{1}{2}((K-S_N{)}^{+}+(K-S_N^{-}{)}^{+})$. Dies ist in dem angegebenen Matlab-Programm realisiert.

% Monte-Carlo-Simulation fuer einen europaeischen Put
% mittels antithetischen Variablen
clear all, randn(’state’,10)
K = 100; r = 0.05; sigma = 0.2; T = 1;
N = 50; h = 1/N;
S(1) = 80; S1(1) = S(1);
for j=1:100 % 100
  M=j*100;
  for k=1:M
    for i=2:N
      dW = randn*sqrt(h);
      S(i) = S(i−1)*(1 + r*h + sigma*dW);
      S1(i) = S1(i−1)*(1 + r*h − sigma*dW);
    end
    payoff(j,k) = 0.5*(max(0,K−S(N))+max(0,K−S1(N)));
  end
  V(j) = exp(−r*T)*sum(payoff(j,:))/M;
  fprintf(’V = %f in Simulation Nr. %d\n’, V(j), j)
end
plot(V)
xlabel(’Anzahl der Simulationen ( \times 10^2)’,’FontSize’,12)
ylabel(’Optionspreis’,’FontSize’,12)

Die Abbildung stellt die Entwicklung des Optionspreises in Abhängigkeit der Anzahl der Monte-Carlo-Simulationen dar. Der exakte Wert für die gewählten Parameter beträgt 16.98. Bereits nach etwa 2000 Simulationen ist die Varianz sehr klein. Die Rechenzeit (XEON 2 GHz, Matlab) beträgt für 100 Durchläufe (genauer: 10⁴ Simulationsschritte) ca. 3 Minuten.

5.5 Beispiel: Asiatischer Call im Heston-Modell

In den vorangegangenen Abschnitten haben wir die Techniken kennengelernt, mit denen wir die in Beispiel 5.1 vorgestellte asiatische Option im Heston-Modell bewerten können. Die Aufgabe lautet, den fairen Preis einer asiatischen Call-Option mit Auszahlungsfunktion

$V_0(S_T)={(S_T-\frac{1}{T}\underset{0}{\overset{T}{\int }}{S}_{\tau }d\tau )}^{+}$

zu berechnen, wobei die Dynamik von $S_t$ und ${\sigma }_t$ durch das Heston-Modell

$dS=r_{t}Sdt+\sigma Sd{W}^1,$

$d{\sigma }^2=\kappa (\theta -{\sigma }^2)dt+\nu \sigma d{W}^2$

gegeben ist. Die risikofreie Zinsrate sei zeitabhängig und definiert durch

$r_t=\frac{1}{100}(\sin (2\pi t)+t+3),0\le t\le T.$

Der Wiener-Prozess $(W^1,W^2)$ sei $𝒩(0,\mathrm{\Sigma })$-verteilt mit der Kovarianzmatrix

$\mathrm{\Sigma }=\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right).$

Die Parameter seien

$\kappa =2,\theta =0.4,\nu =0.2,T=1.$

In Beispiel 5.4 haben wir eine Formel zur Berechnung einer zweidimensionalen $𝒩(\mu ,\mathrm{\Sigma })$-verteilten Zufallsvariablen hergeleitet. Seien $Z^1,Z^2$ standardnormalverteilte Zufallsvariable. Dann ist

$\left(\begin{array}{c}{W}^1\\ W^2\end{array}\right)=\left(\begin{array}{cc}1& 0\\ \rho & \sqrt{1-{\rho }^2}\end{array}\right)\left(\begin{array}{c}{Z}^1\\ Z^2\end{array}\right)=\left(\begin{array}{c}{Z}^1\\ \rho Z^1+\sqrt{1-{\rho }^2}{Z}^2\end{array}\right)$

$𝒩(0,\mathrm{\Sigma })$-verteilt. Wie in Abschnitt 5.1 erläutert, wird der Optionspreis approximativ über die Formel

$V_0=\exp (-\underset{0}{\overset{T}{\int }}{r}_{t}dt)\frac{1}{M}\sum _{k=1}^M{(S_T^k-\frac{1}{N}\sum _{i=1}^N{S}_{t_i}^k)}^{+}$

berechnet, wobei $M$ die Anzahl der Monte-Carlo-Simulationen und $N$ die Anzahl der Zeitschritte ist.

Die Berechnung erfolgt also in drei Schritten:

• 1. Schritt: Berechne $d{W}^1$ und $d{W}^2$ aus

$d{W}_k^1=Z_k^1\sqrt{h},d{W}_k^2=\rho Z_k^1\sqrt{h}+\sqrt{1-{\rho }^2}{Z}_k^2\sqrt{h}.$

• 2. Schritt: Löse das SDE-System mit dem Euler-Maruyama-Verfahren:

$({\sigma }^2{)}_{i+1}^k=({\sigma }^2{)}_i^k+\kappa (\theta -({\sigma }^2{)}_i^k)h+\nu {\sigma }_i^{k}d{W}_k^2,$

$S_{i+1}^k=S_i^k(1+r(t_i)h+{\sigma }_i^{k}d{W}_k^1),i=1,...,N-1.$

• 3. Schritt: Berechne die Approximation des Optionspreises:

$\stackrel{k}{\stackrel{‾}{S}}=\frac{1}{N}\sum _{i=1}^N{S}_i^k,V_0=\exp (-\underset{0}{\overset{T}{\int }}{r}_{t}dt)\frac{1}{M}\sum _{k=1}^M{(S_N^k-\stackrel{k}{\stackrel{‾}{S}})}^{+}.$

Wir erhalten das weiter unten angegebene Matlab-Programm zur Preisbestimmung einer asiatischen Option im Heston-Modell.

Mit diesen Parametern erhalten wir die in der unteren Tabelle dargestellten Werte. Der ”faire” Preis der asiatischen Option beträgt also etwa 13.6.

Welchen Effekt hat die stochastische Volatilität? Der Preis einer asiatischen Option mit konstanter Volatilität $\sigma =0.25$ beträgt (bei 200 000 Simulationen) $V=6.5$. Warum ist diese Option deutlich preiswerter als diejenige mit stochastischer Volatilität? Wir sehen, dass die Volatilität überwiegend größer als der Startwert $\sigma =0.25$ ist. Dies begründet den höheren Preis. Berechnen wir den Preis einer asiatischen Option mit konstanter Volatilität $\sigma =0.55$, so erhalten wir bei 200 000 Simulationen $V=13.3$ – ein Wert, der nahe bei dem Wert der asiatischen Option im Heston-Modell liegt.
Wir sehen in unterer Tabelle, dass die Monte-Carlo-Methode recht langsam konvergiert und sehr viele Simulationen für halbwegs genaue Werte notwendig sind. Es ist also sinnvoll, ein Verfahren höherer Ordnung als das Euler-Maruyama-Verfahren zu wählen, etwa das Milstein-Verfahren (5.18) (vgl. Abschnitt 5.3):

$y_{k+1}^s=y_k^s+a^s(y_k,t_k)h+\sum _{j=1}^m{b}^{sj}(y_k,t_k)\mathrm{\Delta }{W}^j+\sum _{j,\nu =1}^m\sum _{l=1}^n\frac{\partial b^{s\nu }}{\partial x_l}(y_k,t_k)\underset{t_k}{\overset{t_{k+1}}{\int }}\underset{t_k}{\overset{\tau }{\int }}d{W}_z^{j}d{W}_{\tau }^{\nu },s=1,2,$

wobei

$y_k=\left(\begin{array}{c}{\sigma }_k^2\\ S_k\end{array}\right),a(y_k,t_k)=\left(\begin{array}{c}\kappa (\theta -{\sigma }_k^2)\\ r_{t_k}{S}_k\end{array}\right),b(y_k,t_k)=\left(\begin{array}{cc}\nu {\sigma }_k& 0\\ 0& {\sigma }_k{S}_k\end{array}\right)$

$\begin{array}{cc}𝐒𝐢𝐦𝐮𝐥𝐚𝐭𝐢𝐨𝐧𝐞𝐧& 𝐎𝐩𝐭𝐢𝐨𝐧𝐬𝐩𝐫𝐞𝐢𝐬\\ \text{(Anzahl)}& \text{(berechnet)}\\ 1000& 14.42\\ 5000& 12.65\\ 10000& 13.06\\ 50000& 13.51\\ 100000& 13.49\\ 200000& 13.60\end{array}$

und $\frac{\partial b^{s\nu }}{\partial x_l}$ die partiellen Ableitungen nach $x_1=S$ bzw. $x_2={\sigma }^2$ bezeichnet. Eine (etwas langwierige Rechnung) ergibt

$S_{k+1}=S_k+r_{t_k}{S}_{k}h+{\sigma }_k{S}_k\mathrm{\Delta }{W}_k^2+\frac{1}{2}{\sigma }_k^2{S}_k((\mathrm{\Delta }{W}_k^2{)}^2-h)+\frac{\nu }{2}{S}_k{I}_{12},$

${\sigma }_{k+1}^2={\sigma }_k^2+\kappa (\theta -{\sigma }_k^2)h\nu {\sigma }_k\mathrm{\Delta }{W}_k^1+\frac{{\nu }^2}{4}((\mathrm{\Delta }{W}_k^1{)}^2-h).$

Wie in Abschnitt 5.3 erläutert, können die Integrale $I_{21}$ und $I_{12}$ mit dem Verfahren (5.21) approximiert werden.

Versuchen Sie, den Algorithmus zu implementieren.

% Bestimmung des Preises einer asiatischen Option
% im Heston-Modell
clear all, randn(’state’,2)
M = 1000; % Anzahl der Simulationen
N = 100; % Anzahl der Zeitschritte
T = 1; h = T/N;
S_0 = 100; sigma2_0 = 0.25*0.25; % Startwerte
kappa = 2; theta = 0.4; nu = 0.2; rho = 0.2;
% zeitabhaengige Zinsrate und Integral von 0 bis T
t=0:h:T;
r = 0.01*(sin(2*pi*t) + t) + 0.03;
intr = T*(T/200 + 0.03) − cos(2*pi*T)/(200*pi) + 1/(200*pi);
% zweidimensionaler Wiener-Prozess
dW1 = randn(M,N+1)*sqrt(h);
dW2 = rho*dW1 + sqrt(1−rho^2)*randn(M,N+1)*sqrt(h);
% Initialisierung von S und sigma^2
S = S_0*ones(M,N+1);
sigma2 = sigma2_0*ones(M,N+1);
% Loesung des SDE-Systems mit dem Euler-Maruyama-Verfahren
for i=1:N
  sigma2(:,i+1) = sigma2(:,i) + kappa*(theta−sigma2(:,i))*h . . . + nu*sqrt(sigma2(:,i)).*dW2(:,i);
  S(:,i+1) = S(:,i).*(1 + r(:,i)*h + sqrt(sigma2(:,i)).*dW1(:,i));
end
payoff = max(0,S(:,N+1)−mean(S,2));
V = exp(−intr)*mean(payoff)

Literatur

[1] Burrage, K., Burrage, P.M.: High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations. Appl. Numer. Math. 22 (1996), 81-101.

[2] Cox, J., Ross, S., Rubinstein, M.: Option Pricing: A Simplified Approach. J. Financ. Econom. 7 (1979), 228 - 263.

[3] Fisz, M.: Wahrscheinlichkeitsrechnung und mathematische Statistik. Deutscher Verlag der Wissenschaften, Berlin

[4] Günther, M., Jüngel, A.: Finanzderivate mit MATLAB. Vieweg & Sohn, Wiesbaden 2003

[5] Hastings, C.: Approximations for Digital Computers. Princeton University Press, Princeton 1955

[6] Higham, D.: An algorithmic introduction to the numerical solution of stochastic differential equations. SIAM Review 43 (2001), 525-546.

[7] Higham, D.; Kloeden, P.: MAPLE and MATLAB for stochastic differential equations in finance. Preprint, 2002.

[8] Hull, J.C.: Options, Futures, and other Derivates. Prentice Hall 1997.

[9] Klimov, G.: Probability Theory, Mir 1988

[10] Kloeden, P.; Platen, E.: Numerical Solution of Stochastic Differential Equations. Springer, Berlin, 1995.

[11] Korn, R., Korn, E.: Optionsbewertung und Portfolio-Optimierung. Vieweg, Braunschweig 1999

[12] Kwok: Mathematical Models of Financial Derivatives. Springer, Singapur, 1998.

[13] Øksendal, B.: Stochastic Differential Equations. Springer, Berlin 1998

[14] Seydel, R.: Einführung in die numerische Berechnung von Finanzderivaten, Springer, Berlin-Heidelberg-New York 2000.

[15] Wilmott, P., Howison, S., Dewyenne, J.: The Mathematics of Financial Derivatives. Cambridge University Press, Cambridge 1996.

[16] Zhang, P.: Exotic Options, World Scientific, Singapure 1997.