MKZR:Stochastyczne równania różniczkowe

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(→Niesymetryczny Proces Wienera: dyfuzja ze stałym dryftem)
(→Schemat Milsteina)
 
(Nie pokazano 40 wersji pomiędzy niniejszymi.)
Linia 1: Linia 1:
 +
[[Category:MKZR]]
==Stochastyczne równania różniczkowe==
==Stochastyczne równania różniczkowe==
-
 
-
[[Instrumenty_Rynku#Stopy_procentowe|Stopy procentowe]]
 
W tym rozdziale zostaną opisane  metody numeryczne, które służa do rozwiązywania  [[PIZL:Stochastyczne_równania_różniczkowe|stochastycznych równań różniczkowych]] typu:
W tym rozdziale zostaną opisane  metody numeryczne, które służa do rozwiązywania  [[PIZL:Stochastyczne_równania_różniczkowe|stochastycznych równań różniczkowych]] typu:
-
<math>\frac{dX(t)}{dt} = F(X(t), t) + G(X(t), t)\Gamma(t)</math>
+
:<math>\frac{dX(t)}{dt} = F(X(t), t) + G(X(t), t)\Gamma(t)</math>
gdzie F i G to dowolne funkcje, a <math>\Gamma(t)</math> jest procesem losowym.
gdzie F i G to dowolne funkcje, a <math>\Gamma(t)</math> jest procesem losowym.
Linia 11: Linia 10:
Dlatego poprawne jest zapisanie równanie Ito w postaci:
Dlatego poprawne jest zapisanie równanie Ito w postaci:
-
<math>dX(t)= F(X(t), t)dt + G(X(t), t) dW(t)\;</math>
+
:<math>dX(t)= F(X(t), t)dt + G(X(t), t) dW(t)\;</math>
Nie zmienia to ogólności, gdyż jak wiadomo każde równanie zapisane w interpretacji Stratonowicza ma swój odpowiednik Ito. Dla potrzeb metod numerycznych będziemy rozpatrywać zawsze równania Ito, a jeśli pojawią się równania Stratonowicza to będziemy je transpormować do postaci Ito.  
Nie zmienia to ogólności, gdyż jak wiadomo każde równanie zapisane w interpretacji Stratonowicza ma swój odpowiednik Ito. Dla potrzeb metod numerycznych będziemy rozpatrywać zawsze równania Ito, a jeśli pojawią się równania Stratonowicza to będziemy je transpormować do postaci Ito.  
Linia 18: Linia 17:
-
=== Schemat Eulera dla równań stochastycznych  ===
+
=== Schemat Eulera-Maruyamy dla równań stochastycznych  ===
Najprostszą metodą aproksymacji numerycznej równania stochastycznego jest podobnie jak w przypadku równań różniczkowych zwyczajnych jest schemat Eulera. Część deterministyczną równania stochastycznego traktujemy w taki sam sposób jak w schemacie Eulera dla równań różniczkowych zwyczajnych. Niech h oznacza krok  całkowania i oś czasowa będzie zdyskretyzowana na przedzialy <math>t_{i-1},t_{i},t_{i+1}</math> oraz <math>h=t_{i}-t_{i-1}</math>.  Wtedy część deterministyczna równania stochastycznego przyjmuje postać:
Najprostszą metodą aproksymacji numerycznej równania stochastycznego jest podobnie jak w przypadku równań różniczkowych zwyczajnych jest schemat Eulera. Część deterministyczną równania stochastycznego traktujemy w taki sam sposób jak w schemacie Eulera dla równań różniczkowych zwyczajnych. Niech h oznacza krok  całkowania i oś czasowa będzie zdyskretyzowana na przedzialy <math>t_{i-1},t_{i},t_{i+1}</math> oraz <math>h=t_{i}-t_{i-1}</math>.  Wtedy część deterministyczna równania stochastycznego przyjmuje postać:
-
<math>X(t_i) = X(t_{i-1}) + \int_{t_{i-1}}^{t_{i}} F(X(t), t) dt \simeq X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h  </math>
+
:<math>X(t_i) = X(t_{i-1}) + \int_{t_{i-1}}^{t_{i}} F(X(t), t) dt \simeq X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h  </math>
Aby całkować część stochastyczną potrzebujemy formuły na przyrost skończony procesu Wienera:
Aby całkować część stochastyczną potrzebujemy formuły na przyrost skończony procesu Wienera:
-
<math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt =\;\int_{t_{i-1}}^{t_{i}} dW(t) = W(t_{i})-W(t_{i-1})</math>
+
:<math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt =\;\int_{t_{i-1}}^{t_{i}} dW(t) = W(t_{i})-W(t_{i-1})</math>
Wiemy, że proces Wienera jest procesem o [[PIZL:Proces_Wienera_i_proces_dyfuzji#PROCES_WIENERA_W.28t.29|przyrostach niezależnych]], które są  gaussowską zmienna losową o zerowej wartości średniej i wariancji   
Wiemy, że proces Wienera jest procesem o [[PIZL:Proces_Wienera_i_proces_dyfuzji#PROCES_WIENERA_W.28t.29|przyrostach niezależnych]], które są  gaussowską zmienna losową o zerowej wartości średniej i wariancji   
Linia 34: Linia 33:
Tak więc widać, że w schemacie Eulera całkę typu <math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt </math> należy zastąpić w każdym kroku całkowania gaussowską zmienną losową o wariancji proporcjonalnej do kroku całkowania h. Ponieważ z reguły dysponujemy gaussowskich generatorem liczb losowych  o jednostkowej wariancji N(0,1), można zapisać:
Tak więc widać, że w schemacie Eulera całkę typu <math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt </math> należy zastąpić w każdym kroku całkowania gaussowską zmienną losową o wariancji proporcjonalnej do kroku całkowania h. Ponieważ z reguły dysponujemy gaussowskich generatorem liczb losowych  o jednostkowej wariancji N(0,1), można zapisać:
-
<math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt = \sqrt{2 h} N(0,1) </math>
+
:<math> \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt = \sqrt{2 h} N(0,1) </math>
Ten zapis pokazuje też ważną cechę przy obliczaniu aproksymacji rozwiązań stochastycznych - najniższy rząd w h jest nie O(h) ale <math>O(h^{1/2})</math>.
Ten zapis pokazuje też ważną cechę przy obliczaniu aproksymacji rozwiązań stochastycznych - najniższy rząd w h jest nie O(h) ale <math>O(h^{1/2})</math>.
Linia 41: Linia 40:
-
Korzystając w powyższych faktów, możemy zapisać pełny schemat Eulera dla równania stochastycznego (Ito):
+
Korzystając w powyższych faktów, możemy zapisać pełny schemat Eulera-Maruyamy dla równania stochastycznego (Ito):
-
<math>X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h  + \sqrt{h}G(X(t_{t-1}), t_{t-1})  N(0,1)</math>.
+
:<math>X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h  + \sqrt{h}G(X(t_{t-1}), t_{t-1})  N(0,1)</math>.
Gaussowskie zmienne losowe możemy otrzymać np. korzystając z [[MKZR:Liczby_losowe#Algorytm_Boxa-Mullera|algorytmu Box-a Mullera]].
Gaussowskie zmienne losowe możemy otrzymać np. korzystając z [[MKZR:Liczby_losowe#Algorytm_Boxa-Mullera|algorytmu Box-a Mullera]].
-
=== Schemat Eulera dla układu równań stochastycznych  ===
+
=== Schemat Eulera-Maruyamy dla układu równań stochastycznych  ===
-
Schemat Eulera można uogólnić na układy równań stochastycznych. Niech <math>\mathbf X(t)</math> będzie wektorem o składowych   
+
Schemat Eulera-Maruyamy można uogólnić na układy równań stochastycznych. Niech <math>\mathbf X(t)</math> będzie wektorem o składowych   
:<math>\mathbf X(t)=( X^1(t),X^2(t),...,X^n(t))</math>.  
:<math>\mathbf X(t)=( X^1(t),X^2(t),...,X^n(t))</math>.  
Układ równań stochastycznych (Ito) można zapisać w ogólnej postaci:
Układ równań stochastycznych (Ito) można zapisać w ogólnej postaci:
Linia 58: Linia 57:
gdzie <math>W^i(t),\;W^j(t)</math>  są niezależnymi procesami Wienera dla <math>i\neq j</math>,  
gdzie <math>W^i(t),\;W^j(t)</math>  są niezależnymi procesami Wienera dla <math>i\neq j</math>,  
-
<math>F^i</math> oznacza wektor drytfu a  <math>G^{i,j}</math> jest macierzą <math>n \times n</math> funkcji.
+
<math>F^i</math> oznacza wektor dryftu a  <math>G^{i,j}</math> jest macierzą <math>n \times n</math> funkcji.
-
 
+
-
Wtedy schemat Eulera ma postać:
+
-
 
+
-
<math> X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\;  </math>
+
-
 
+
-
=== Proces Wienera ===
+
-
 
+
-
Proces Wienera jest rozwiązaniem następującego stochastycznego równania różniczkowego:
+
-
 
+
-
<math>dX(t)= 2 D dW(t)\;</math>.
+
-
 
+
-
 
+
-
Jego realizacja jest funkcją ciągłą, ale jednocześnie nigdzie nie jest różniczkowalna.
+
-
[[Plik:wiener.png|thumb|360px|Proces Wienera]]
+
-
Stosując schemat Eulera można wygenerować pojedyńczą trajektorię takiego procesu:
+
-
<source lang='matlab'>
+
-
h=0.01;
+
-
N=400;
+
-
x(1)=0;
+
-
D=1;
+
-
for i=2:N
+
-
  x(i)=x(i-1)+ sqrt(2*D*h)* normrnd (0,1,1,1);
+
-
endfor
+
-
 
+
-
plot((1:N)*h,x,'-')
+
-
</source>
+
-
 
+
-
Jest to trajektoria, która w granicy <math>h\to\infty</math> jest funkcją nigdzie nie różniczkowalną. Numerycznie przejawia się to w tym, że niezależnie od tego jak małe h weźmiemy do symulacji, wykres procesu Wienera zawsze będzie wyglądał na "zaszumiony". Można też korzystając z techniki Mostu Browna z danego przybliżenia procesu Wienera  dla  pewnego h wysymulowac przybliżenie dla h/2. Również wtedy można stwierdzić, że wykres nigdy nie będzie wizualnie "gładki", co w tym przypadku oznacza matematycznie nieciągłość w każdym punkcie.
+
-
 
+
-
=== Proces Wienera: rozkład P(x,t) ===
+
-
Proces Wienera (symetryczny) spełnia, jak wiadomo równanie dyfuzji:
+
-
 
+
-
<math>\frac{\partial P(x, t)}{\partial t}  = D \frac{\partial^2 P(x, t)}{\partial x^2}    </math>
+
-
 
+
-
Rozważmy przypadek w którym cząstka staruje w punkcie <math>x=0</math> w czasie <math>t=0</math>. W języku gęstości prawdopodobieństwa oznacza to
+
-
 
+
-
<math>P(x, 0) = \delta(x)\;</math>
+
-
 
+
-
Jego rozwiązaniem na równania dyfuzji prostej z takim warunkiem początkowym jest rozkład Gaussa:
+
-
+
-
<math>P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{x^2}{4Dt}\right]</math>
+
-
 
+
-
w którym odchylenie standardowe jest proporcjonalne do czasu.
+
-
Zweryfikujmy ten fakt numerycznie. W tym celu wysymulujemy trajektorie 20000 procesów Wienera, dla 100 kroków. Poprzedni program łatwo można zmodyfikować by obliczenia były przeprowadzone dnie dla skalara x(1),x(2),... ale dla wektorów x(:,1),x(:,2),... Jeśli zastosujemy wbudowany generator normalnych liczb losowych '''normrnd''' to możemy także skorzystać z możliwości wygenerowania wielu liczb za jednym wywolaniem np. '''normrnd (0,1,M,1)'''. Postępując w ten sposób linijka:
+
-
<source lang='matlab'>
+
-
x(:,i)=x(:,i-1) + sqrt(2*D*h)* normrnd (0,1,M,1);
+
-
</source>
+
-
przedstawia wekrorowy zapis jednego kroku całkowania M stoschastychnych równań różniczkowych. Nie tylko upraszcza to zapis, ale także przyśpiesza obliczenia w przypadku stosowania interpretowanego języka  jakim jest GNU Octave lub Matlab.
+
-
[[Plik:wiener_Pxt.png|thumb|360px|Proces Wienera]]
+
-
Cały program może wyglądać tak:
+
-
<source lang='matlab'>
+
-
clear all
+
-
close all
+
-
N=100;
+
-
M=20000;
+
-
T=4;
+
-
h=T/N;
+
-
clear x
+
-
x=zeros(M,N);
+
-
D=.81;
+
-
for i=2:N
+
-
  x(:,i)=x(:,i-1) + sqrt(2*D*h)* normrnd (0,1,M,1);
+
-
endfor
+
-
</source>
+
-
 
+
-
a procedura rysująca wykres histogramu położeń M cząstek po czasie t i porównująca
+
-
ja z odpowiednim rozkładem Gaussa:
+
-
<source lang='matlab'>
+
-
n=10
+
-
t=n/N*T
+
-
xmax=4;
+
-
h1=.2;
+
-
hist(x(:,n),[-xmax:h1:xmax],1/h1)
+
-
hold on
+
-
fplot(@(xx) normpdf(xx,0,sqrt(2*D*t)),[-xmax,xmax],200,'ro-')
+
-
hold off
+
-
</source>
+
-
 
+
-
Czyli w czasie funkcja rozkładu jest rozpływającym się po całej prostej rozkładem Gaussa. Poniższa procedura rysuje rozkłady w czterech momentach czasu:
+
-
[[Plik:wiener_Pxt3.png|thumb|360px|Rozkład P(x,t) dla Procesu  Wienera w czterech następujących po sobie chwilach. Widać rozwywanie się piku Gaussowskiego.]]
+
-
<source lang='matlab'>
+
-
for idx=1:4
+
-
  n=5+(idx-1)*20
+
-
  t=n/N*T
+
-
  subplot(2,2,idx)
+
-
  xmax=10;
+
-
  h1=.2;
+
-
  hist(x(:,n),[-xmax:h1:xmax],1/h1)
+
-
  hold on
+
-
  fplot(@(xx) normpdf(xx,0,sqrt(2*D*t)),[-24,24],200,'ro-')
+
-
  hold off
+
-
  legend ( sprintf("t=%2.0f",t) )
+
-
endfor
+
-
</source>
+
-
 
+
-
 
+
-
=== Niesymetryczny Proces Wienera: dyfuzja ze stałym dryftem ===
+
-
 
+
-
Dyfuzja z dryfem jest granicznym przypadkiem [[PIZL:Proces_Wienera_i_proces_dyfuzji#Przypadek_niesymetryczny:_dyfuzja_z_dryfem|niesymetryczngo błądzenia przypadkowego]].
+
-
+
-
Równanie dyfuzji z dryfem ma postać:
+
-
<math>\frac{\partial P(x, t)}{\partial t}  = -V \frac{\partial P(x, t)}{\partial x} + D \frac{\partial^2 P(x, t)}{\partial x^2}    </math>
+
Wtedy schemat Eulera-Maruyamy ma postać:
-
i jest ono równoważne stochstycznemu równaniu różniczkowemu:
+
:<math> X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\; </math>
-
<math>dX(t)= Vdt + 2 D dW(t)\;</math>.
+
-
Łatwo możemy zmodyfikować program symulujący proces symmetryczny na niesymetryczny:
+
=== Schemat Milsteina ===
-
[[Plik:wiener_dryf.png|thumb|360px|Niesymetryczny proces Wienera:dyfuzja z dryfem]]
+
-
<source lang='matlab'>
+
-
close all
+
-
clear all
+
-
h=0.01;
+
-
N=1200;
+
-
x(1)=0;
+
-
D=0.9;
+
-
V=1.2;
+
-
for i=2:N
+
-
  x(i)=x(i-1)+ V*h +sqrt(2*D*h)* normrnd (0,1,1,1);
+
-
endfor
+
-
plot((1:N)*h,x,'-',(1:N)*h,(1:N)*h*V,'-')
+
-
</source>
+
 +
Schemat Milsteina jest dany wzorem interacyjnym:
 +
:<math>\displaystyle X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h - </math>
 +
:<math>
 +
\frac{1}{2}G(X(t_{i-1},t-1)G'(X(t_{i-1},t-1) h  + \sqrt{h}G(X(t_{t-1}), t_{t-1})  N(0,1)</math>.
-
Rozważmy przypadek w którym cząstka staruje w punkcie <math>x=0</math> w czasie <math>t=0</math>. W języku gęstości prawdopodobieństwa oznacza to
+
W stosunku do schematy Eulera-Maruyamy zawiera on dodatkowy składnik, proporcjonalny do O(h):
-
<math>P(x, 0) = \delta(x)\;</math>
+
:<math>-  \frac{1}{2}G(X(t_{i-1},t-1)G'(X(t_{i-1},t-1) h</math>.
-
Jego rozwiązaniem na równania dyfuzji prostej z takim warunkiem początkowym jest rozkład Gaussa:
+
Ta poprawka powoduje, że powyższy schemat jest pierwszego rzędu w sensie silnym w przeciwieństwie do schematu  Eulera-Maruyamy, ktry jest rzędu  1/2.
-
   
+
-
<math>P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{x^2}{4Dt}\right]</math>
+

Aktualna wersja na dzień 20:35, 31 mar 2016

Spis treści

Stochastyczne równania różniczkowe

W tym rozdziale zostaną opisane metody numeryczne, które służa do rozwiązywania stochastycznych równań różniczkowych typu:

\[\frac{dX(t)}{dt} = F(X(t), t) + G(X(t), t)\Gamma(t)\]

gdzie F i G to dowolne funkcje, a \(\Gamma(t)\) jest procesem losowym. Najczęstszym przypadek to taki w którym \(\Gamma(t)\) to biały szum Gaussowski. Tak zapisane równanie nie jest precyzyjnie określone ze względu na dylemat Stratonowicza-Ito. Dlatego poprawne jest zapisanie równanie Ito w postaci:

\[dX(t)= F(X(t), t)dt + G(X(t), t) dW(t)\;\]

Nie zmienia to ogólności, gdyż jak wiadomo każde równanie zapisane w interpretacji Stratonowicza ma swój odpowiednik Ito. Dla potrzeb metod numerycznych będziemy rozpatrywać zawsze równania Ito, a jeśli pojawią się równania Stratonowicza to będziemy je transpormować do postaci Ito.



Schemat Eulera-Maruyamy dla równań stochastycznych

Najprostszą metodą aproksymacji numerycznej równania stochastycznego jest podobnie jak w przypadku równań różniczkowych zwyczajnych jest schemat Eulera. Część deterministyczną równania stochastycznego traktujemy w taki sam sposób jak w schemacie Eulera dla równań różniczkowych zwyczajnych. Niech h oznacza krok całkowania i oś czasowa będzie zdyskretyzowana na przedzialy \(t_{i-1},t_{i},t_{i+1}\) oraz \(h=t_{i}-t_{i-1}\). Wtedy część deterministyczna równania stochastycznego przyjmuje postać:

\[X(t_i) = X(t_{i-1}) + \int_{t_{i-1}}^{t_{i}} F(X(t), t) dt \simeq X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h \]


Aby całkować część stochastyczną potrzebujemy formuły na przyrost skończony procesu Wienera:

\[ \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt =\;\int_{t_{i-1}}^{t_{i}} dW(t) = W(t_{i})-W(t_{i-1})\]

Wiemy, że proces Wienera jest procesem o przyrostach niezależnych, które są gaussowską zmienna losową o zerowej wartości średniej i wariancji \(2(t_{i} − t_{i-1})=2 h\).

Tak więc widać, że w schemacie Eulera całkę typu \( \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt \) należy zastąpić w każdym kroku całkowania gaussowską zmienną losową o wariancji proporcjonalnej do kroku całkowania h. Ponieważ z reguły dysponujemy gaussowskich generatorem liczb losowych o jednostkowej wariancji N(0,1), można zapisać:

\[ \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt = \sqrt{2 h} N(0,1) \]

Ten zapis pokazuje też ważną cechę przy obliczaniu aproksymacji rozwiązań stochastycznych - najniższy rząd w h jest nie O(h) ale \(O(h^{1/2})\).

Ponadto z takiego sformułowania widać też, że zmiany procesu Wienera w stosunku do przyrostów czasu są rozbieżne w granicy \(h\to 0\).


Korzystając w powyższych faktów, możemy zapisać pełny schemat Eulera-Maruyamy dla równania stochastycznego (Ito):

\[X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h + \sqrt{h}G(X(t_{t-1}), t_{t-1}) N(0,1)\].


Gaussowskie zmienne losowe możemy otrzymać np. korzystając z algorytmu Box-a Mullera.

Schemat Eulera-Maruyamy dla układu równań stochastycznych

Schemat Eulera-Maruyamy można uogólnić na układy równań stochastycznych. Niech \(\mathbf X(t)\) będzie wektorem o składowych \[\mathbf X(t)=( X^1(t),X^2(t),...,X^n(t))\]. Układ równań stochastycznych (Ito) można zapisać w ogólnej postaci:

\[d X^i(t)= F^i(\mathbf X(t), t)dt + \sum_{j=1}^{n} G^{i,j}(\mathbf X(t), t) dW^j(t)\;\]\; j=1,2,...,n,

gdzie \(W^i(t),\;W^j(t)\) są niezależnymi procesami Wienera dla \(i\neq j\), \(F^i\) oznacza wektor dryftu a \(G^{i,j}\) jest macierzą \(n \times n\) funkcji.

Wtedy schemat Eulera-Maruyamy ma postać:

\[ X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\; \]

Schemat Milsteina

Schemat Milsteina jest dany wzorem interacyjnym:

\[\displaystyle X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h - \] \[ \frac{1}{2}G(X(t_{i-1},t-1)G'(X(t_{i-1},t-1) h + \sqrt{h}G(X(t_{t-1}), t_{t-1}) N(0,1)\].

W stosunku do schematy Eulera-Maruyamy zawiera on dodatkowy składnik, proporcjonalny do O(h):

\[- \frac{1}{2}G(X(t_{i-1},t-1)G'(X(t_{i-1},t-1) h\].

Ta poprawka powoduje, że powyższy schemat jest pierwszego rzędu w sensie silnym w przeciwieństwie do schematu Eulera-Maruyamy, ktry jest rzędu 1/2.

Źródło „http://172.16.3.1/ekonofizyka/index.php/MKZR:Stochastyczne_r%C3%B3wnania_r%C3%B3%C5%BCniczkowe”