Milstein verfahren+ Monte Carlo

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Bree

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 10.06.2010, 16:52 Titel: Milstein verfahren+ Monte Carlo

Hallo zusammen,
also ich habe ein echt großes Problem und ich hoffe dass ich hier Hilfe finde finde, denn so langsam weiß ich nicht mehr wen ich noch fragen soll.

Es geht um Stochastische Differentialgleichungen. Das Programm welches ihr unten seht zeigt das Milstein verfahren zur Approximation von SDE's.
Ich soll nun eine Monte Carlo Simulation drumherum laufen lassen und das ganze dann noch für verschiedene Schrittweiten testen, so dass man am ende die Konvergenzordnung bestimmen kann.
Ich habe aber KEINE Ahnung ie ich das machen soll.

Naja ich erwarte keine Wunder aber irgendetwas was mir vielleicht nen Anschubs gibt, so dass ich neue Ideen habe.

Viele Grüße Bree

Code:

close all; %Milstein Verfahren für 1 SDE und 2 Brownsche Bewegungen
randn('state',2)

M=1000;
N1=100;
N2=100;
T=1;
h=T/N1;
dt=h/N2;
t=0:dt:T;
n=(N1*N2)+1;
s=0:h:T;

sigma=zeros(2,1);

X0=1; mu=2; sigma(1)=0.25; sigma(2)=0.3;

% zwei Brownsche Bewegungen
dW=sqrt(dt)*randn(2,(N1*N2)+1);
W=cumsum(dW,2);

% exakte Lösung
Xexakt=X0*exp((mu - 0.5*(sigma(1)^2 + sigma(2)^2))*t + sigma(1).*W(1,:) + sigma(2).*W(2,:));

zk=zeros(2,N1*N2+1);
z=zeros(2,N1*N2+1);
X = zeros(1,N1+1);
X(1) = X0;

exakt=zeros(1,N1+1);

for k=1:N1
for i=1:N2
zk(1,i+1)=zk(1,i) + zk(2,i)*dW(1,i); %stochastische Doppelintegrale
zk(2,i+1)=zk(2,i) + dW(2,i);
z(1,i+1)=z(1,i) + z(2,i)*dW(2,i);
z(2,i+1)=z(2,i) + dW(1,i);
end

%Milstein
X(k+1)= X(k)*(1 + mu*h + sigma(1)*dW( 1, 1+(k-1)*N2 ) + sigma(2)*dW( 2, 1+(k-1)*N2 ) + 0.5*sigma(1)^2*(dW(1, 1+(k-1)*N2 )^2-h) + 0.5*sigma(2)^2*( dW(2, 1+(k-1)*N2 )^2-h)...
+ sigma(1)*sigma(2)*(zk(1,N2)+ z(1,N2)));

exakt(k)=Xexakt(1+(k-1)*N2);

end

Xerror=abs(exakt-X);

%plot(t, Xexakt,'r--'), hold on
plot(0:h:T, exakt,'r--'), hold on
plot(0:h:T, X,'b'), hold on
plot(s, Xerror, 'g')

Funktion ohne Link?

Maddy

Ehrenmitglied


	Beiträge: 494

	Anmeldedatum: 02.10.08

	Wohnort: Greifswald

	Version: ---

Verfasst am: 11.06.2010, 10:04 Titel:

Also vll seh ich das Problem nicht, aber den Code zu einer "function" mit entsprechenden Eingabe und Ausgabe-Parametern erweitern und dann einfach im Monte-Carlo-Code aufrufen.
_________________

>> why
The computer did it.

Bree

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 11.06.2010, 11:47 Titel:

Also wie meinst du das genau?

ich kann ja nicht einfach

m=1:1000

um das ganze programm laufen lassen, oder? Und wenn doch, wie mache ich das?

Maddy

Ehrenmitglied


	Beiträge: 494

	Anmeldedatum: 02.10.08

	Wohnort: Greifswald

	Version: ---

Verfasst am: 11.06.2010, 12:39 Titel:

Also mal als Beispiel für den Parameter M:

Code:

function []=Milsteinverfahren(M) % M als Input, ohne Output

close all; %Milstein Verfahren für 1 SDE und 2 Brownsche Bewegungen
randn('state',2)

% M=1000; % wird als Input Variable übergeben, hier also auskommentiert
N1=100;
N2=100;
T=1;
h=T/N1;
dt=h/N2;
t=0:dt:T;
n=(N1*N2)+1;
s=0:h:T;

sigma=zeros(2,1);

X0=1; mu=2; sigma(1)=0.25; sigma(2)=0.3;

% zwei Brownsche Bewegungen
dW=sqrt(dt)*randn(2,(N1*N2)+1);
W=cumsum(dW,2);

% exakte Lösung
Xexakt=X0*exp((mu - 0.5*(sigma(1)^2 + sigma(2)^2))*t + sigma(1).*W(1,Smile + sigma(2).*W(2,Smile);

zk=zeros(2,N1*N2+1);
z=zeros(2,N1*N2+1);
X = zeros(1,N1+1);
X(1) = X0;

exakt=zeros(1,N1+1);

for k=1:N1

for i=1:N2

zk(1,i+1)=zk(1,i) + zk(2,i)*dW(1,i); %stochastische Doppelintegrale
zk(2,i+1)=zk(2,i) + dW(2,i);

z(1,i+1)=z(1,i) + z(2,i)*dW(2,i);
z(2,i+1)=z(2,i) + dW(1,i);

end

%Milstein
X(k+1)= X(k)*(1 + mu*h + sigma(1)*dW( 1, 1+(k-1)*N2 ) + sigma(2)*dW( 2, 1+(k-1)*N2 ) + 0.5*sigma(1)^2*(dW(1, 1+(k-1)*N2 )^2-h) + 0.5*sigma(2)^2*( dW(2, 1+(k-1)*N2 )^2-h)...
+ sigma(1)*sigma(2)*(zk(1,N2)+ z(1,N2)));

exakt(k)=Xexakt(1+(k-1)*N2);

end

Xerror=abs(exakt-X);

%plot(t, Xexakt,'r--'), hold on
plot(0:h:T, exakt,'r--'), hold on
plot(0:h:T, X,'b'), hold on
plot(s, Xerror, 'g')]%

end

Funktion ohne Link?

speichern unter dem Namen MilsteinVerfahren.m
und dann einfach in einem Script diese Funktion aufrufen;

Code:

M=1000; % für mehrere Werte einfach eine Schleife, wobei dann der Output der Funktion verändert werden sollte
MilsteinVerfahren(M)

Funktion ohne Link?

_________________

>> why
The computer did it.

Bree

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 14.06.2010, 16:08 Titel:

Danke hab es jetzt aber ein wenig anders gemacht. Das mit der Monte Carlo Simulation geht jetzt.

Hier das Programm:

Code:

close all;
randn('state',2);

% Monte-Carlo
M = 1000;

% Schrittweite Milstein
N1 = 2^9;
% Zeit
T = 1;
h = T/N1;
t = 0:h:T;

% Schrittweite Euler - Doppelintegrale
N2=2^9;
dt=h/N2;

% Initialwerte
X0 = 1;

X = zeros(1,N1+1);
X(1) = X0;

% Drift
mu = 0;

% Diffusion
sigma = zeros(2,1);
sigma(1) = 0.25;
sigma(2) = 0.3;

for m = 1 : M

% zwei Brownsche Bewegungen
dW = sqrt(h) * randn(2,N1);
W = sum(dW,2);

% exakte Lösung - fuer Endwert T
Xexakt = X0 * exp((mu - 0.5*(sigma(1)^2 + sigma(2)^2))*T + sigma(1).*W(1,:) + sigma(2).*W(2,:));

% fuer Doppelintegrale
zk = zeros(2,N2+1);
z = zeros(2,N2+1);
dB = sqrt(dt) * randn(2,N2);

% Pfad
for k = 1 : N1
% Doppelintegrale bitte ueberpruefen
for i = 1 : N2
zk(1,i+1)= zk(1,i) + zk(2,i)*dB(1,i);
zk(2,i+1)= zk(2,i) + dB(2,i);

z(1,i+1)= z(1,i) + z(2,i)*dB(2,i);
z(2,i+1)= z(2,i) + dB(1,i);
end

% Milstein
X(k+1) = X(k).*(1 + mu*h + sigma(1)*dW(1,k) + sigma(2)*dW(2,k)...
+ 0.5*sigma(1)^2*(dW(1,k)^2 - h) + 0.5*sigma(2)^2*(dW(2,k)^2-h)...
+ sigma(1)*sigma(2)*(zk(1,N2+1)+ z(1,N2+1)));
end

% Fehler am Ende des Pfads: T
Xerror(m,1) = abs(Xexakt - X(N1+1));
end

Funktion ohne Link?

Aber ich soll das jetzt für verschiedene Schrittweiten testen. Also jetzt haben wir ja ein Schrittweite von 2^-9.
Ich soll das nun auch für Schrittweiten 2^-8, 2^-7 bis 2^-4 und das ganze dann so ausgeben:

Code:

dtlist=h*(2.^(0:6-1));
loglog(dtlist,mean(Xerror),'*--'),hold on
loglog(dtlist,(dtlist.^(0.5)),'--'),hold off

A=[ones(6,1),log(dtlist)'];
b=log(mean(Xerror)');
x=A\b;
gamma=x(2),residum=norm(A*x-b);

Funktion ohne Link?

Die Schleife dafür sieht so aus, aber das funktioniert überhaupt nicht:

Code:

for p=1:6
R=2^(p-1);
h2=R*h; %aktuelle Schrittweite
L=n/R; % ANZAHL der Euler Schritte

for j=1:L

Funktion ohne Link?

Keine Ahnung wie ich das noch in das Programm mit rein bringen soll.

KiSoSa

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 12.04.2011, 14:12 Titel:

Kann mir vielleicht jemand von euch erklären, warum man bei der Milstein-Implementierung für sigma=1 dann +0,5*sigma^2*((dW)^2-h) als zusätzlichen Term hat und nicht wie im Verfahren sonst 0,5*sigma*sigma'*((dW)^2-h) hat!!!


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.