Metoda Eulera-Maruyamy

Niech dane będzie stochastyczne równanie różniczkowe postaci

${\displaystyle \mathrm {d} X_{t}=a(X_{t})\,\mathrm {d} t+b(X_{t})\,\mathrm {d} W_{t}}$

oraz

${\displaystyle X_{T_{init}}=x_{0}}$ – zadany warunek początkowy (położenie cząstki w chwili początkowej ${\displaystyle T_{\mathrm {init} }}$),

gdzie:

${\displaystyle W_{t}}$ – proces Wienera,

${\displaystyle a(X_{t}),b(X_{t})}$ – zadane funkcje procesu losowego ${\displaystyle X_{t}.}$

Poszukiwane jest rozwiązanie tego równania na przedziale ${\textstyle [T_{init},T_{end}].}$

Analiza równania stochastycznego

O stochastycznym charakterze powyższego równania decyduje różniczka procesu losowego Wienera: gdyby była równa zeru, to proces byłby deterministyczny. Funkcja ${\displaystyle b,}$ stojąca przed tą różniczka, pełni zaś rolę wzmacniania lub osłabiania efektu fluktuacji losowych.

Z drugiej strony, funkcja ${\displaystyle a,}$ stojąca przed różniczką ${\displaystyle dt,}$ decyduje, na ile proces deterministyczny dominuje nad procesem losowym – gdyby funkcja ta była zerowa, to proces byłby czysto losowy.

Powyższe własności łatwo sprawdzić, rozwiązując podane dalej przykładowe równanie stochastyczne, zmieniając odpowiednio jego parametry.

Przybliżeniem Eulera-Maruyamy dokładnego rozwiązania ${\displaystyle X_{t}}$ powyższego równania jest łańcuch Markowa ${\displaystyle Y_{n}}$ taki że:

• zdarzenia losowe łańcucha ${\displaystyle Y_{n},}$ ${\displaystyle n=0,1,2,\dots ,N,}$ zachodzą w dyskretnych chwilach czasu ${\displaystyle T_{\mathrm {init} }=\tau _{0}<\tau _{1}<\ldots <\tau _{N}=T_{\mathrm {end} },}$

odległych od siebie o wartość skokową

${\displaystyle \Delta t=(T_{\mathrm {end} }-T_{\mathrm {init} })/N,}$

• ${\displaystyle Y_{0}=x_{0}}$ – zadane zdarzenie początkowe,
• ${\displaystyle Y_{n+1},}$ ${\displaystyle n=0,1,2,\dots ,N-1}$ – zdarzenia obliczane wzorem rekurencyjnym:

${\displaystyle Y_{n+1}=Y_{n}+a(Y_{n})\,\Delta t+b(Y_{n})\,\Delta W_{n},}$

przy czym:

${\displaystyle \Delta W_{n}=W_{n+1}-W_{n},}$

${\displaystyle a(Y_{n}),b(Y_{n})}$ – zadane funkcje ${\displaystyle Y_{n}.}$

Cechą procesu Wienera jest, że różnica zmiennych losowych ${\displaystyle W_{t}-W_{s}}$ dla dowolnych chwil czasu ${\displaystyle s,t,}$ takich że ${\displaystyle 0\leqslant s\leqslant t}$, ma rozkład normalny o średniej równej zeru i wariancji równej odległości czasowej rozważanych chwil ${\displaystyle t-s}$. Wynika stąd, że zmienna losowa ${\displaystyle \Delta W_{n}=W_{n+1}-W_{n}}$ ma rozkład normalny o średniej równej zeru i wariancji ${\displaystyle \Delta t.}$

Niech dane będzie stochastyczne równanie różniczkowe jednej zmiennej ${\displaystyle X_{t}}$

${\displaystyle dX_{t}=\theta \cdot (\mu -X_{t})dt+\sigma \cdot dW_{t},}$

o warunku początkowym

${\displaystyle X_{0}=0,}$

gdzie ${\displaystyle \theta ,\mu ,\sigma }$ – stałe parametry; różniczka ${\displaystyle dW_{t}}$ odpowiada za proces losowy Wienera.

Powyższe równanie opisuje tzw. proces Ornsteina-Uhlenbecka.

Zgodnie z metodą Eulera-Maruyamy zamienia się to równanie na postać dyskretną

${\displaystyle \Delta Y_{n}=\theta \cdot (\mu -Y_{n})\cdot \,\Delta t+\sigma \cdot \Delta W_{n}.}$

Ponieważ ${\displaystyle \Delta Y_{n}=Y_{n+1}-Y_{n},}$ to otrzymuje się równania rekurencyjne

${\displaystyle Y_{n+1}=Y_{n}+\theta \cdot (\mu -Y_{n})\cdot \,\Delta t+\sigma \cdot \Delta W_{n},\qquad n=1,2,\dots ,N,}$

gdzie ${\displaystyle \Delta W_{n}}$ – zmienne losowe o rozkładzie normalnym i wariancji ${\displaystyle \Delta t.}$

Warunek początkowy przyjmie postać:

${\displaystyle Y_{0}=0.}$

Zmieniając parametry równania stochastycznego, można otrzymać ruch układu w zależności od różnych warunków otoczenia, z jakim układ oddziałuje. Np.

• przyjmując wartość parametru ${\displaystyle \sigma =0,}$ otrzyma się ruch czysto deterministyczny,
• przyjmując wartość parametru ${\displaystyle \theta =0,}$ otrzyma się ruch czysto losowy,
• parametr ${\displaystyle \mu }$ określa położenie, do którego zmierza układ po dłuższym czasie, jeżeli czynniki losowe nie są zbyt duże, tj. dla odpowiednio małej wartości ${\displaystyle \sigma$ ;} na rysunku pokazano cztery trajektorie układu w takich warunkach; mimo różnych wartości położeń początkowych układu (na rys. oznaczonych literą a), układ zmierza po pewnym czasie do położenia ${\displaystyle X=\mu .}$ Przykładem jest ruch wahadła poddanego działaniu losowego oddziaływania od otoczenia i jednocześnie tłumionego – niezależnie od położenia początkowego wahadło przyjmie ostatecznie najniższe położenie.

Poniżej podano przykład kodu w języku Python służącego do całkowania równania Ornsteina-Uhlenbecka.

Na ilustracji pokazano wynik symulacji komputerowej dla przyjętych w programie wartości parametrów.

import numpy as np
import matplotlib.pyplot as plt

#Dane:
liczba_symulacji = 5
t_0 = 0    # chwila początkowa
t_1 = 7    # chwila końcowa
N   = 200  # liczba punktów podziału przedziału

theta  = 0.7
mu     = 1.4
sigma  = 0.06

# Definicja funkcji generującej liczbę losową z rozkładu normalnego o wariancji Dt
def DW(Dt): 
    return np.random.normal(loc=0.0,scale=np.sqrt(Dt))

# Część główna programu
Dt   = (t_1 - t_0) / N

# Tablica dyskretnych chwil czasu
t    = np.arange(t_0, t_1, Dt) 

# Inicjalizacja tablicy na kolejne położenia cząstki
Y    = np.zeros(N)

for _ in range(liczba_symulacji):
    Y[0] = 0
    for i in range(0, N-1):
        Y[i+1]=Y[i] + theta*(mu -Y[i])*Dt +sigma * DW(Dt)
    plt.plot(t, Y) # tworzenie wykresu

plt.grid(True)  # Dodanie siatki do wykresu
plt.show()