Metoda Eulera-Maruyamy

Metoda Eulera-Maruyamy – metoda numerycznego rozwiązywania stochastycznych równań różniczkowych. Nazwa metody pochodzi od nazwiska japońskiego matematyka Gisiro Maruyamy (1916-1986), który rozszerzył metodę Eulera z 1768, stosowaną do numerycznego rozwiązywania równań różniczkowych zwyczajnych. Metoda ta jest jedną z niewielu metod numerycznych, które pozwalają metody numeryczne stosowane do równań deterministycznych rozszerzyć do równań stochastycznych.

Opis metody

Równanie stochastyczne

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

$\mathrm{d}{X}_t=a(X_t)\mathrm{d}t+b(X_t)\mathrm{d}{W}_t$

oraz

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

gdzie:

$W_t$ – proces Wienera,

$a(X_t),b(X_t)$ – zadane funkcje procesu losowego $X_t.$

Poszukiwane jest rozwiązanie tego równania na przedziale $[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 $b,$ stojąca przed tą różniczka, pełni zaś rolę wzmacniania lub osłabiania efektu fluktuacji losowych.

Z drugiej strony, funkcja $a,$ stojąca przed różniczką $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żenie Eulera-Maruyamy

Przybliżeniem Eulera-Maruyamy dokładnego rozwiązania $X_t$ powyższego równania jest łańcuch Markowa $Y_n$ taki że:

• zdarzenia losowe łańcucha $Y_n,$ $n=0,1,2,\dots ,N,$ zachodzą w dyskretnych chwilach czasu $T_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}={\tau }_0<{\tau }_1<\dots <{\tau }_N=T_{\mathrm{e}\mathrm{n}\mathrm{d}},$

odległych od siebie o wartość skokową

$\mathrm{\Delta }t=(T_{\mathrm{e}\mathrm{n}\mathrm{d}}-T_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}})/N,$

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

$Y_{n+1}=Y_n+a(Y_n)\mathrm{\Delta }t+b(Y_n)\mathrm{\Delta }{W}_n,$

przy czym:

$\mathrm{\Delta }{W}_n=W_{n+1}-W_n,$

$a(Y_n),b(Y_n)$ – zadane funkcje $Y_n.$

Dowodzi się, że zmienne losowe $\mathrm{\Delta }{W}_n$ jako różnice zmiennych losowych o rozkładzie normalnym (proces Wienera) są niezależne oraz mają rozkład normalny o wartości oczekiwanej zero i wariancji $\mathrm{\Delta }t.$

Przykład: Równanie Ornsteina-Uhlenbecka

Cztery losowe trajektorie procesu Ornsteina-Uhlenbecka dla θ = 1, σ = $\sqrt{2}:$
niebieska: położenie początkowe a = 10, μ = 0
pomarańczowa: położenie początkowe a = 0, μ = 0
zielona: położenie początkowe a = −10, μ = 0
czerwona: położenie początkowe a = 0, μ = −10
E(X) określa trajektorię, po której poruszałby się układ, gdyby czynniki losowe były zerowe (σ = 0).

Równanie stochastyczne

Niech dane będzie stochastyczne równanie różniczkowe jednej zmiennej $X_t$

$d{X}_t=\theta \cdot (\mu -X_t)dt+\sigma \cdot d{W}_t,$

o warunku początkowym

$X_0=0,$

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

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

Zamiana na postać dyskretna

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

$\mathrm{\Delta }{Y}_n=\theta \cdot (\mu -Y_n)\cdot \mathrm{\Delta }t+\sigma \cdot \mathrm{\Delta }{W}_n.$

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

$Y_{n+1}=Y_n+\theta \cdot (\mu -Y_n)\cdot \mathrm{\Delta }t+\sigma \cdot \mathrm{\Delta }{W}_n,n=1,2,\dots ,N,$

gdzie $\mathrm{\Delta }{W}_n$ – zmienne losowe o rozkładzie normalnym i wariancji $\mathrm{\Delta }t.$

Warunek początkowy przyjmie postać:

$Y_0=0.$

Analiza zależność ruchu układu od parametrów równania

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

• przyjmując wartość parametru $\sigma =0,$ otrzyma się ruch czysto deterministyczny,
• przyjmując wartość parametru $\theta =0,$ otrzyma się ruch czysto losowy,
• parametr $\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 $\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 $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.

Kod programu

Poniżej podano przykład kodu w języku Python, całkujący równanie Ornsteina-Uhlenbecka. Program można testować, korzystając np. z darmowego notatnika colab google online.

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

Pięć przykładowych rozwiązań stochastycznego równania różniczkowego opisującego proces Ornsteina-Uhlenbecka, znalezionych numerycznie za pomocą metody Eulera-Maruyamy dla θ = 0.7, σ = 0.06, μ = 1.4 i identycznych położeń początkowych X= 0; na osi pionowej jest położenie układu, na osi poziomej czas. Każde rozwiązanie przedstawia możliwą, losowo przebiegającą trajektorią układu, która dla dłuższych czasów oscyluje wokół punktu X=μ =1.4

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, która generuje liczbę losową przy każdym wywołaniu
def DW(Dt): 
    return np.random.normal(loc = 0.0, scale = np.sqrt(Dt))

# Część główna programu
Dt   = (t_1 - t_0) / N
t    = np.arange(t_0, t_1, Dt) # tablica dyskretnych chwil czasu
Y    = np.zeros(N)             # tablica dyskretnych położeń cząstki

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)

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

Zobacz też

• równanie Langevina

Bibliografia

• Szablon:Cytuj
• P. Jackel, Monte Carlo Methods in Finance, Wiley, 2002, Szablon:ISBN.