Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 08:23, 6 maj 2014

Proces Wienera

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

$dX(t)= \sqrt{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 stoschastycznych 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 +\sqrt{2 D} dW(t)\;$.

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

Niesymetryczny proces Wienera:dyfuzja z dryfem. Zielonym kolorem jest zaznaczona prosta x=Vt, która jest dryfem procesu

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

Dyfuzja z dryftem: rozkład P(x,t)

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

Podobnie jak w przypadku symetrycznego procesu Wienera rozważmy cząstkę starującą 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 z dryfem prostej z takim warunkiem początkowym jest rozkład Gaussa:

$P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{(x-Vt)^2}{4Dt}\right],$ który ma zależne od czasu odchylenie standardowe jak i wartość średnią.

clear all
close all
N=100;
M=20000;
T=4;
h=T/N;
clear x
x=zeros(M,N);
D=.81;
V=5.0
for i=2:N
  x(:,i)=x(:,i-1) + V*h + sqrt(2*D*h)* normrnd (0,1,M,1);
endfor
hold on
for idx=1:4
  n=5+(idx-1)*20
  t=n/N*T
  xmax=20;
  h1=.2;
  hist(x(:,n),[-2:h1:xmax],1/h1)
  fplot(@(xx) normpdf(xx,V*t,sqrt(2*D*t)),[-2,xmax],200,'ro-')
endfor
hold off

Proces Ornsteina Uhlenbecka

Jest to proces który spełnia równanie stochastyczne:

Dwadzieścia trajektorii procesu Ornsteina Uhlenbecka dla parametrów D=.01,k=1.0 i $\mu=1.0$. Dla czasów 0..2 widoczna jest relaksacja do stanu stacjonarnego, dla większych czasów fluktuacje wokół tego stanu.

$dX(t)= -k ( X(t)-\mu ) dt + \sqrt{2 D } dW(t)\;$.

Jest to proces z linowym dryfem. Posiada on stan stacjonarny (w przeciwieństwie do procesu Wienera). Jest modelem ciała zawieszonego sprężyście poddanego równowagowym fluktuacjom termicznym. Parametr $\mu$ reprezentuje równowagową (stacjonarną) wartośc średnią procesu a k jest szybkością relaksacji do stanu równowagi.

clear all
close all
N=500;
M=20;
T=14;
h=T/N;
clear x
x=zeros(M,N);
D=.01;
k=1.0;
mu=1.0;
for i=2:N
  x(:,i)=x(:,i-1) - k*(x(:,i-1)-mu)* h + sqrt(2*D*h)* normrnd (0,1,M,1);
endfor
plot((1:N)*h,x,'r-')

Geometryczny proces Wienera

Geometryczny proces Wienera znajduje zastosowanie w ekonomii jako model ceny akcji (patrz: Przykłady zastosowań równań stochastycznych w ekonomii).

Definicja

Geometryczny proces Wienera jest procesem losowym, który jest rozwiązaniem równania stochastycznego Ito:

$dX(t) = \mu X(t) dt + \sigma X(t) d W(t)\,$.

Geometryczny proces Wienera, dwadzieścia trajektorii. Widoczny jest eksponencjalny wzrost zmiennej losowej z pewnymi fluktuacjami.

Część deterministyczna część tego równania jest liniowa, podobnie jak w przypadku Ornsteina-Uhlenbecka. Wprzypadku $\sigma=0$ rozwiązanie przedstawia exponencjalne rosnącą funkcję czasu:

$x(t)\simeq e^{\mu t}$,

Można się spodziewać, że zabużenie równania szumem będzie odwzierciedlało tą cechę.

Symulacja numeryczna

Człon stochastyczny jest proporcjonalny o wartości procesu, czyli mamy do czynienia z tzw. "szumem multiplikatywnym". Ponieważ równanie jest w interpretacji Ito, można je bezpośrednio rozwiązać numerycznie. Interpretacja Ito dla $\displaystyle \sigma X(t) d W(t)$ oznacza, że w schemacie aproksymacyjnym bierzemy wartość procesu X(t) "przed skokiem". W takim przypadku możemy zastosować schemat numeryczny Eulera podobnie jak np. dla procesu Ornsteina-Uhlenbecka. Krok czasowy może być zapisany w postaci wektorowej jako:

  x(:,i)=x(:,i-1) + mu*x(:,i-1)*h + sigma*sqrt(h)*x(:,i-1).*normrnd (0,1,M,1);

Jako warunek początkowy dla symulacji należy przyjąć wartość x>0. Łatwo zauważyć własność równania definiującego proces, że startując w $x(0)=0$ rozwiązaniem jest funkcja stała $\displaystyle x(t)=0$. Z drugiej strony równanie jest symetryczne ze względu transformację $x\to -x$ i ponadto jeśli x oznacza cenę to powinno mieć wartość dodatnią. Dlatego zawsze bierzemy x>0 jako warunek początkowy. Poniższy program generuje 20 trajektori geometrycznego procesu Wienera:

clear all
close all
N=400;
M=20;
T=14;
h=T/N;
clear x
x=zeros(M,N);
x(:,1)=1*ones(M,1); # log(1)=0
sigma=.1;
mu=0.1;
for i=2:N
  x(:,i)=x(:,i-1) + mu*x(:,i-1)*h + sigma*sqrt(h)*x(:,i-1).*normrnd (0,1,M,1);
endfor
plot((1:N)*h,x,'r-')

Rozwiązanie analityczne: porównanie z symulacją

Równanie definiujące geometryczny proces Wienera można przetransformować do równania na proces Wienera z dryfem. Transformacją jest:

$$\displaystyle Y=\log(X).$$

Łatwo to zauważyć dzieląc równanie przez $x(t)$ i korzystając z faktu, że

$$\displaystyle d \log(X)= dX/X$$

i korzystając z formuły Ito na dla funkcji log(x) otrzymujemy:

otrzymujemy:

$$d \log(X(t)) = (\mu-\frac{1}{2}\sigma^2) dt + \sigma d W(t).$$

Tak więc geometryczny proces Wienera jest równoważny procesowi Wienera z dryfem dla $\log(X(t))$. Gęstośc pradwopodobieństwo przejścia ze stanu $x(t_0)$ do $x(t)$ w czasie $t-t_0$ dla procesu Wienera z dryfem:

$$Y(t)=\log(X(t))$$

wynosi:

$$P_y(y,t|y_0,0)= \frac{1}{\sqrt{2\pi\sigma^2 t}} e^{-\displaystyle\frac{(y-y_0(\mu-\frac{1}{2}\sigma^2)t)^2}{2\sigma^2t}}$$

Korzystając ze reguł transformacji zmiennych losowych dla funkcji log, która jest różnowartościowa:

$$P_y(y)=\frac{P_x(x)}{|g'(x)|}$$

mamy:

$$P_x(x)= \frac{P_y(\log(x))}{|x|}$$

tak więc ostatecznie otrzymujemy formułę:

$$P_x(x,t|x_0,0)= \frac{1}{\sqrt{2\pi\sigma^2 t x^2}} e^{-\displaystyle\frac{(\log(x/x_0)-(\mu-1/2\sigma^2)^2}{2\sigma^2t}}.$$

Zmienna losowa w geometrycznym procesie Wienera po czasie $t$ ma rozkład logarytmicznie normalny.

Rozkład P(x,t) dla geometrycznego procesu Wienera w czterech następujących po sobie chwilach t = 0.7,3.5,6.2,9.1. Widać rozmywanie się piku zgodnie z odpowiednim rozkładem logarytmicznie normalnym.

Porównajmy ten wzrór z histogramami dużej licznby realizacji w kolejnych chwilach czasu. W systemie GNU Octave/matlab funkcja rozkładu logarytmicznie normalnego jest dana jako normrnd. W poniższym programie startujemy z punktu $x(t=0)=1$ dla którego $\log(x)=0$ czyli odpowiedni process Wienera z dryfem miałby warunek początkowy $y(t=0)=0$.

clear all
close all
N=100;
M=20000;
T=14;
h=T/N;
clear x
x=zeros(M,N);
x(:,1)=1*ones(M,1); # log(1)=0
sigma=.1;
mu=0.1;
for i=2:N
  x(:,i)=x(:,i-1) + mu*x(:,i-1)*h + sigma*sqrt(h)*x(:,i-1).*normrnd (0,1,M,1);
endfor
 
hold on
for idx=1:4
  n=5+(idx-1)*20;
  t=n/N*T
  xmax=5;
  h1=.1;
  hist(x(:,n),[0:h1:xmax],1/h1)
  fplot(@(xx) lognpdf(xx,(mu-0.5*sigma^2)*t,sigma*sqrt(t)),[0.0,xmax],200,'ro-')
endfor
hold off