確率微分方程式の離散化

以下の形式の伊藤過程を考える。

\[ X_t = X_0 + \int_0^t \mu_s\ ds + \int_0^t \sigma_s\ dW_s \]

これを微分形式で記述すると、以下の確率微分方程式(SDE)となる。

\[ dX_t = \mu_t\ dt + \sigma_t\ dW_t \]

伊藤過程の見本路はどんな精度で拡大しても細かく振幅しているはずのものであり、そのまま描画することはできない。そこで、時間 \([0, T]\)\(N\) 等分して離散化することを考える。時間の刻み幅は \(\Delta t\) と表記する。

\[ \Delta t = T / N \]

オイラー・丸山法(Euler-Maruyama approximation)では、伊藤過程を以下のように近似する。 ここで、\(\{\epsilon_t\}_t\) は標準正規分布 \(N(0, 1)\) に従う互いに独立な確率変数である。

\[ \hat{X}_{t+\Delta t} = \hat{X}_t + \mu_t \Delta t + \sigma_t \sqrt{\Delta t} \epsilon_t \]

この近似は、ウィーナー過程の増分 \(dW_t = W_{t+\Delta t} - W_t\) が正規分布 \(N(0, \Delta t)\) に従うという性質に基づいている。この増分の分散は \(\Delta t\) であるから、標準正規分布に乗じる係数を \(\sqrt{\Delta t}\) とすべき点に注意する。

ウィーナー過程 \(W_t\) 自体は、伊藤過程においてドリフト項 \(\mu_t = 0\)、拡散項 \(\sigma_t = 1\) とした最も基本的なケースである。この過程のオイラー・丸山法による離散近似は以下のように表現できる。

\[ \hat{W}_{t+\Delta t} = \hat{W}_t + \sqrt{\Delta t} \epsilon_t \]

この過程は時系列分析におけるランダムウォークにほかならない。

以下では、ウィーナー過程を離散化することで得られたランダムウォークの見本路を作成してプロットする。

library(ggplot2)
theme_set(theme_bw())

set.seed(42)

N <- 252
dt <- 1 / N
time <- seq(0, 1, by = dt)

# 標準正規乱数の生成と累積和の計算
df <- list()
for (i in 1:25) {
  dW <- sqrt(dt) * rnorm(N)
  path <- cumsum(c(0, dW))
  df[[length(df) + 1]] <-
    data.frame(id = i, time=time, path=path)
}
df <- do.call(rbind, df)

# プロット
ggplot(df) +
  geom_line(aes(time, path, group=id)) +
  labs(
    title = "Discretized Wiener Process",
    x = "Time (t)", y = "W(t)"
  )

ウィーナー過程のシミュレーション(R)
import numpy as np
import pandas as pd
from plotnine import *
theme_set(theme_bw())

np.random.seed(42)

N = 252
dt = 1 / N
time = np.linspace(0, 1, N + 1)

# 標準正規乱数の生成と累積和の計算
df = list()
for i in range(25):
  dW = np.sqrt(dt) * np.random.normal(0, 1, N)
  path = np.cumsum(np.insert(dW, 0, 0))
  df.append(
    pd.DataFrame({'id':i, 'time':time, 'path':path})
  )
df = pd.concat(df)

# プロット
ggplot(df) + \
  geom_line(aes("time", "path", group="id")) + \
  labs(
    subtitle = "Discretized Wiener Process",
    x = "Time (t)", y = "W(t)"
  )

ウィーナー過程のシミュレーション(Python)