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 {init} }$),

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 ${\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 $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 {init} }=\tau _{0}<\tau _{1}<\ldots <\tau _{N}=T_{\mathrm {end} },$

odległych od siebie o wartość skokową

$\Delta t=(T_{\mathrm {end} }-T_{\mathrm {init} })/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})\,\Delta t+b(Y_{n})\,\Delta W_{n},$

przy czym:

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

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

Dowodzi się, że zmienne losowe $\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 $\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}$

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

o warunku początkowym

$X_{0}=0,$

gdzie $\theta ,\mu ,\sigma$ – stałe parametry; różniczka $dW_{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ą

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

Ponieważ $\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 \,\Delta t+\sigma \cdot \Delta W_{n},\qquad n=1,2,\dots ,N,$

gdzie $\Delta W_{n}$ – zmienne losowe o rozkładzie normalnym i wariancji $\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

• PaoloP. Brandimarte PaoloP., Numerical Methods in Finance: A Matlab Introduction, New York: Wiley, 2002, ISBN 0-471-39686-9, OCLC 47092043 .brak strony (książka)
• P. Jackel, Monte Carlo Methods in Finance, Wiley, 2002, ISBN 0-471-49741-X.