Z Skrypty dla studentów Ekonofizyki UPGOW

Stochastyczne równania różniczkowe

Stopy procentowe

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 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 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 dla układu równań stochastycznych

Schemat Eulera 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 drytfu a $G^{i,j}$ jest macierzą $n \times n$ funkcji.

Wtedy schemat Eulera 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)\;$

Proces Wienera

Proces Wienera jest rozwiązaniem następującego stochastycznego równania różniczkowego:

$dX(t)= 2 D dW(t)\;$.

Jego realizacja jest funkcją ciągłą, ale jednocześnie nigdzie nie jest różniczkowalna.

Proces Wienera

Stosując schemat Eulera można wygenerować pojedyńczą trajektorię takiego procesu:

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,'-')

Jest to trajektoria, która w granicy $h\to\infty$ 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:

$\frac{\partial P(x, t)}{\partial t} = D \frac{\partial^2 P(x, t)}{\partial x^2}$

Rozważmy przypadek w którym cząstka staruje w punkcie $x=0$ w czasie $t=0$. W języku gęstości prawdopodobieństwa oznacza to

$P(x, 0) = \delta(x)\;$

Jego rozwiązaniem na równania dyfuzji prostej z takim warunkiem początkowym jest rozkład Gaussa:

$P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{x^2}{4Dt}\right]$

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:

x(:,i)=x(:,i-1) + sqrt(2*D*h)* normrnd (0,1,M,1);

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.

Proces Wienera

Cały program może wyglądać tak:

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

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:

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

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:

Rozkład P(x,t) dla Procesu Wienera w czterech następujących po sobie chwilach. Widać rozwywanie się piku Gaussowskiego.

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

Niesymetryczny Proces Wienera: dyfuzja ze stałym dryftem

Dyfuzja z dryfem jest granicznym przypadkiem niesymetryczngo błądzenia przypadkowego.

Równanie dyfuzji z dryfem ma postać:

$\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}$

i jest ono równoważne stochstycznemu równaniu różniczkowemu: $dX(t)= Vdt + 2 D dW(t)\;$.

Łatwo możemy zmodyfikować program symulujący proces symmetryczny na niesymetryczny:

Niesymetryczny proces Wienera:dyfuzja z dryfem

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,'-')

Rozważmy przypadek w którym cząstka staruje w punkcie $x=0$ w czasie $t=0$. W języku gęstości prawdopodobieństwa oznacza to

$P(x, 0) = \delta(x)\;$

Jego rozwiązaniem na równania dyfuzji prostej z takim warunkiem początkowym jest rozkład Gaussa:

$P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{x^2}{4Dt}\right]$