Z Skrypty dla studentów Ekonofizyki UPGOW
m (→Proces Wienera) |
m (→Proces Wienera) |
||
| Linia 1: | Linia 1: | ||
| + | [[Category:MKZR]] | ||
[[Category:MKZR]] | [[Category:MKZR]] | ||
=== Proces Wienera === | === Proces Wienera === | ||
Wersja z 15:08, 23 sty 2011
Spis treści |
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.
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.
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:
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:
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)
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:
\(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)\,\).
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,t_0)= \frac{1}{\sqrt{2\pi\sigma^2 t}} e^{-\displaystyle\frac{(y-(\mu-1/2\sigma^2)^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,t_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.
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