метод Мильштейна

В математике метод Мильштейна представляет собой метод приближенного численного решения стохастического дифференциального уравнения . Он назван в честь Григория Н. Мильштейна , который впервые опубликовал его в 1974 году.

Описание

Рассмотрим автономное Ито стохастическое дифференциальное уравнение : $$\mathrm {d} X_{t}=a(X_{t})\,\mathrm {d} t+b(X_{t})\,\mathrm {d} W_{t}$$ с начальным состоянием $X_{0}=x_{0}$, где $W_{t}$ обозначает винеровский процесс , и предположим, что мы хотим решить это СДУ на некотором интервале времени $[0,T]$. Тогда приближение Мильштейна к истинному решению $X$ это цепь Маркова $Y$ определяется следующим образом:

• разделить интервал $[0,T]$ в $N$ равные подинтервалы ширины $\Delta t>0$: $$0=\tau _{0}<\tau _{1}<\dots <\tau _{N}=T{\text{ with }}\tau _{n}:=n\Delta t{\text{ and }}\Delta t={\frac {T}{N}}$$
• набор $Y_{0}=x_{0};$
• рекурсивно определять $Y_{n}$ для $1\leq n\leq N$ к: $$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)$$ где $b'$ обозначает производную ​ $b(x)$ относительно $x$ и: $$\Delta W_{n}=W_{\tau _{n+1}}-W_{\tau _{n}}$$ являются независимыми и одинаково распределенными нормальными случайными величинами с нулевым математическим ожиданием и дисперсией. $\Delta t$. Затем $Y_{n}$ будет приближаться $X_{\tau _{n}}$ для $0\leq n\leq N$и увеличение $N$ даст лучшее приближение.

Обратите внимание, что когда $b'(Y_{n})=0$, т. е. диффузионный член не зависит от $X_{t}$, этот метод эквивалентен методу Эйлера–Маруямы .

Схема Мильштейна имеет как слабый, так и сильный порядок сходимости. $\Delta t$, что превосходит метод Эйлера–Маруямы , который, в свою очередь, имеет такой же слабый порядок сходимости, $\Delta t$, но уступает сильному порядку сходимости, ${\sqrt {\Delta t}}$.

Интуитивный вывод

Для этого вывода мы будем рассматривать только геометрическое броуновское движение (GBM), стохастическое дифференциальное уравнение которого имеет вид: $$\mathrm {d} X_{t}=\mu X\mathrm {d} t+\sigma XdW_{t}$$ с реальными константами $\mu$ и $\sigma$. Используя лемму Ито, получаем: $$\mathrm {d} \ln X_{t}=\left(\mu -{\frac {1}{2}}\sigma ^{2}\right)\mathrm {d} t+\sigma \mathrm {d} W_{t}$$

Таким образом, решение GBM SDE: $${\begin{aligned}X_{t+\Delta t}&=X_{t}\exp \left\{\int _{t}^{t+\Delta t}\left(\mu -{\frac {1}{2}}\sigma ^{2}\right)\mathrm {d} t+\int _{t}^{t+\Delta t}\sigma \mathrm {d} W_{u}\right\}\\&\approx X_{t}\left(1+\mu \Delta t-{\frac {1}{2}}\sigma ^{2}\Delta t+\sigma \Delta W_{t}+{\frac {1}{2}}\sigma ^{2}(\Delta W_{t})^{2}\right)\\&=X_{t}+a(X_{t})\Delta t+b(X_{t})\Delta W_{t}+{\frac {1}{2}}b(X_{t})b'(X_{t})((\Delta W_{t})^{2}-\Delta t)\end{aligned}}$$ где $$a(x)=\mu x,~b(x)=\sigma x$$

См. численное решение, представленное выше для трех различных траекторий.

Численное решение только что представленного стохастического дифференциального уравнения: дрейф в два раза превышает коэффициент диффузии.

Компьютерная реализация

Следующий код Python реализует метод Мильштейна и использует его для решения СДУ, описывающего геометрическое броуновское движение, определяемое формулой $${\begin{cases}dY_{t}=\mu Y\,{\mathrm {d} }t+\sigma Y\,{\mathrm {d} }W_{t}\\Y_{0}=Y_{\text{init}}\end{cases}}$$

# -*- coding: utf-8 -*-
# Milstein Method

import numpy as np
import matplotlib.pyplot as plt


class Model:
    """Stochastic model constants."""
    μ = 3
    σ = 1


def dW(Δt):
    """Random sample normal distribution."""
    return np.random.normal(loc=0.0, scale=np.sqrt(Δt))


def run_simulation():
    """ Return the result of one full simulation."""
    # One second and thousand grid points
    T_INIT = 0
    T_END = 1
    N = 1000 # Compute 1000 grid points
    DT = float(T_END - T_INIT) / N
    TS = np.arange(T_INIT, T_END + DT, DT)

    Y_INIT = 1

    # Vectors to fill
    ys = np.zeros(N + 1)
    ys[0] = Y_INIT
    for i in range(1, TS.size):
        t = (i - 1) * DT
        y = ys[i - 1]
        dw = dW(DT)

        # Sum up terms as in the Milstein method
        ys[i] = y + \
            Model.μ * y * DT + \
            Model.σ * y * dw + \
            (Model.σ**2 / 2) * y * (dw**2 - DT)

    return TS, ys


def plot_simulations(num_sims: int):
    """Plot several simulations in one image."""
    for _ in range(num_sims):
        plt.plot(*run_simulation())

    plt.xlabel("time (s)")
    plt.ylabel("y")
    plt.grid()
    plt.show()


if __name__ == "__main__":
    NUM_SIMS = 2
    plot_simulations(NUM_SIMS)

См. также

• Метод Эйлера-Маруямы

Дальнейшее чтение

• Клоден П.Е. и Платен Э. (1999). Численное решение стохастических дифференциальных уравнений . Шпрингер, Берлин. ISBN  3-540-54062-8 . {{cite book}}: CS1 maint: несколько имен: список авторов ( ссылка )