Milşteyn yöntemi

Matematiğin bir alt dalı olan stokastik analizde Milşteyn yöntemi (ya da Milstein yöntemi) stokastik diferansiyel denklemlere yaklaşık sayısal çözümler üretmek için geliştirilmiş yöntemlerden biridir. Yöntem adını, bu yaklaşıklığı 1974'te yayınlayan Grigori Milşteyn'den (Мильштейн, Г. Н.) almaktadır.

Yöntemin açıklaması

W_t, Wiener süreci olsun. Başlangıç koşulu X₀ = x₀ olan ve

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

hâlindeki stokastik bir diferansiyel denklemi [0, T] aralığında sayısal olarak çözmek istediğimizi varsayalım. Stokastik diferansiyel denklemin teorik çözümü $X$ ile gösterilsin. $X$'in Milstein yaklaşıklığı aşağıdaki gibi kurulup tanımlanan ve bir Markov zinciri olan $Y$dir.

• [0, T] aralığını N tane eşit altaralığa bölelim. Her altaralığın uzunluğu $\Delta t>0$ ile gösterilsin. O zaman,

$0=\tau _{0}<\tau _{1}<\cdots <\tau _{N}=T{\text{ ve }}\Delta t=T/N$

olacaktır.

• Y₀ = x₀ tanımlansın.
• 0 ≤ n ≤ N-1 için $\Delta W_{n}=W_{\tau _{n+1}}-W_{\tau _{n}}$ olmak üzere, Y_n terimleri yinelemeli olarak

$$Y_{n+1}=Y_{n}+a(Y_{n})\Delta t+b(Y_{n})\Delta W_{n}+{\frac {1}{2}}b(Y_{n})b'(Y_{n})\left((\Delta W_{n})^{2}-\Delta t\right)$$

olarak tanımlansın. Burada, $b'$ ile türev, yani, $b(x)$ fonksiyonunun $x$'e göre türevi, kastedilmiştir. Bahsi geçen ΔW_n rasgele değişkenleri bağımsız ve özdeş dağılmış normal rastgele değişkenlerdir. Her birinin beklenen değeri sıfırdır ve her biri Δt varyansına sahiptir.

Yöntemin arkasındaki fikir

Yöntemin arkasındaki esas fikir $b(X_{t})$ terimine Ito formülünü uygulayıp daha Euler yaklaşımlığını iyileştirmektir. Elbette, burada, bu terimdeki $b$ fonksiyonun Ito formülünün koşullarına uyması beklenmektedir. Ito formülü ile $${\begin{aligned}\mathrm {d} b(X_{t})&=b'(X_{t})\mathrm {d} X_{t}+{\frac {1}{2}}b''(X_{t})b^{2}(X_{t})dt\\&=\left[b'(X_{t})a(X_{t})+{\frac {1}{2}}b''(X_{t})b^{2}(X_{t})\right]dt+b'(X_{t})b(X_{t})dW_{t}\\&=\mu _{b}(X_{t})dt+\sigma _{b}(X_{t})dW_{t}\\\end{aligned}}$$ elde edilir. $t\leq u\leq t+\Delta t$ şartını sağlayan bir $u$ sayısı için, Euler-Maruyama yöntemi kullanılarak, yukarıda elde edilen $\mathrm {d} b(X_{t})$ ifadesine yaklaşım sağlanır. Diğer deyişle, $${\begin{aligned}b(X_{u})&\approx b(X_{t})\mu _{b}(X_{t})[u-t]+\sigma _{b}(X_{t})[W_{u}-W_{t}]\\&=b(X_{t})\left[b'(X_{t})a(X_{t})+{\frac {1}{2}}b''(X_{t})b^{2}(X_{t})\right][u-t]+b'(X_{t})b(X_{t})[W_{u}-W_{t}]\end{aligned}}$$ yazılabilir. $W_{u}-W_{t}$ terimi olasılıkta $O({\sqrt {u-t}})$ olduğundan ve bu yaklaşıklıktaki sürüklenme $O(u-t)$ olduğundan, $O(u-t)$ büyümesine sahip terimler düşürülüp $${\begin{aligned}b(X_{u})&\approx b(X_{t})+b'(X_{t})b(X_{t})[W_{u}-W_{t}],\quad u\in [t,t+\Delta t]\end{aligned}}$$ yazılır. Bu yüzden, $${\begin{aligned}\int _{t}^{t+\Delta t}b(X_{u})dW_{u}&\approx \int _{t}^{t+\Delta t}b(X_{t})+b'(X_{t})b(X_{t})[W_{u}-W_{t}]dW(u)\\&=b(X_{t})[W_{t+\Delta t}-W_{t}]+b'(X_{t})b(X_{t})\int _{t}^{t+\Delta t}[W_{u}-W_{t}]dW(u)\\&=b(X_{t})[W_{t+\Delta t}-W_{t}]+{\frac {1}{2}}b'(X_{t})b(X_{t})[(W_{t+\Delta t}-W_{t})^{2}-\Delta t].\\\end{aligned}}$$ elde edilir. O zaman, Euler yaklaşımlığını biraz daha düzelterek $${\begin{aligned}X_{t+\Delta t}\approx X_{t}+a(X_{t})\Delta t+b(X_{t})[W_{t+\Delta t}-W_{t}]+{\frac {1}{2}}b'(X_{t})b(X_{t})[(W_{t+\Delta t}-W_{t})^{2}-\Delta t]\end{aligned}}$$ yazılır.

Programlama uygulaması

$${\begin{cases}dY_{t}=\mu Y\,{\mathrm {d} }t+\sigma Y\,{\mathrm {d} }W_{t}\\Y_{0}=Y_{\text{ilk}}\end{cases}}$$ biçiminde verilmiş olan bir geometrik Brown hareketini sayısal olarak çözmek için Milşteyn yöntemi kullanılabilir.

Rastgele sayıların NumPy tarafından verildiği bir Python programlama kodu aşağıdaki gibi verilebilir:

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from textwrap import wrap

class Model:
    """Stokastik modelin  sabitleri."""
    MU = 0.05
    SIGMA = 0.2

def dW(dt):
    """Normal rastgele değişken çekilişi"""
    return np.random.normal(loc=0.0, scale=np.sqrt(dt))

def simulasyon(t_i:int, t_s:int, adim_sayisi:int):
    """bir simülasyon yolağı"""
    T_ILK = t_i
    T_SON = t_s
    N = adim_sayisi 
    DT = float(T_SON - T_ILK) / N
    Tler = np.arange(T_ILK, T_SON + DT, DT)
    assert Tler.size == N + 1

    Y_ILK = 1

    yler = np.zeros(Tler.size)
    yler[0] = Y_ILK
    for i in range(1, Tler.size):
        t = (i - 1) * DT
        y = yler[i - 1]
        dw = dW(DT)

        #Milstein yontemindeki gibi terimleri topla
        yler[i] = y + \
            Model.MU * y * DT + \
            Model.SIGMA * y * dw + \
            (Model.SIGMA**2 / 2) * y * (dw**2 - DT)

    return Tler, yler

#Birkaç simülasyonun aynı grafikte çizimi
sim_sayisi=100
# [0, 1] aralığında 5000 zaman adımında hesaplıyoruz
t_i=0#
t_s=1
adim_sayisi=5000
for _ in range(sim_sayisi):
    plt.plot(*simulasyon(t_i, t_s, adim_sayisi))

plt.xlabel("zaman (s)")
plt.ylabel("y")
plt.title("\n".join(wrap(
r"$\mu=$ {} ve $\sigma=$ {} parametreleri verilen bir geometrik Brown hareketinin [{},{}] aralığında {} zaman adımlı Milşteyn yöntemiyle yapılan yaklaşıklamalarının {} yolaklı simülasyonu".format(Model.MU, Model.SIGMA, t_i, t_s, adim_sayisi, sim_sayisi))), fontsize='small')
plt.savefig('Milstein.png', format="png")
plt.show()

Ayrıca bakınız

• Euler-Maruyama yöntemi