ブラウン運動のサンプルパスは連続ですが至る所微分不可能であり、2次変分が $(dW)^2 = dt \neq 0$ という通常の微分可能な関数にはない性質を持ちます。このため、ブラウン運動を含む微分方程式を扱うには、通常の微積分（ニュートン=ライプニッツの微積分）を拡張した 確率微積分（stochastic calculus） が必要になります。
1944年、数学者 伊藤清（Kiyosi Ito） は確率積分の理論を構築し、テイラー展開の2次の項が消えないという本質的な発見に基づいて 伊藤の公式（Ito’s formula / Ito’s lemma） を導出しました。この公式は確率解析の最も基本的なツールであり、金融工学のブラック=ショールズ方程式の導出や、物理学の拡散過程の解析に不可欠です。
本記事の内容
- なぜ通常の微積分が使えないか（2次変分の役割）
- 伊藤積分の定義と性質
- 伊藤の公式（補題）の完全な導出
- 幾何ブラウン運動 $dS = \mu S \, dt + \sigma S \, dW$ の解の導出
- オルンシュタイン=ウーレンベック過程
- ランジュバン方程式
- Pythonでオイラー=丸山法によるSDE数値解法
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
なぜ通常の微積分が使えないか
通常の連鎖律の復習
通常の微積分では、$x(t)$ が微分可能な関数で $f(x)$ が2回微分可能であれば、合成関数 $f(x(t))$ の微分は 連鎖律（chain rule） によって、
$$ \frac{d}{dt} f(x(t)) = f'(x(t)) \cdot x'(t) $$
微分形式で書くと、
$$ df = f'(x) \, dx $$
テイラー展開の1次の項だけで十分であり、2次以降の項は無視できます。これは $x(t)$ が微分可能であるとき、$(dx)^2 = (x'(t) dt)^2 = O(dt^2) \to 0$ となるためです。
ブラウン運動の場合の問題
ブラウン運動 $W(t)$ に対しては事情が異なります。前の記事で示した通り、ブラウン運動の2次変分は
$$ (dW)^2 = dt $$
を満たします。つまり、$(dW)^2$ は $dt$ のオーダーであり、0に消えません。したがって、テイラー展開の2次の項が生き残り、通常の連鎖律は修正が必要です。
形式的にまとめると、以下の 確率微積分の乗法規則 が成り立ちます。
| $dt$ | $dW$ | |
|---|---|---|
| $dt$ | 0 | 0 |
| $dW$ | 0 | $dt$ |
$(dt)^2 = 0$（高次の無限小は無視）、$dt \cdot dW = 0$（$dW$ は $O(\sqrt{dt})$ なので $dt \cdot dW = O(dt^{3/2}) \to 0$）、$(dW)^2 = dt$（2次変分）です。
伊藤積分の定義
定義の動機
通常のリーマン=スティルチェス積分 $\int f(t) \, dg(t)$ は、被積分関数 $f$ と積分関数 $g$ の少なくとも一方が有限変分（bounded variation）であることを必要とします。しかし、ブラウン運動の1次変分は無限大であるため、リーマン=スティルチェス積分としては定義できません。
伊藤積分の構成
ステップ1: 単純過程に対する定義
確率過程 $\{f(t)\}$ が 単純過程（simple process） であるとは、分割 $0 = t_0 < t_1 < \cdots < t_n = T$ に対して、
$$ f(t) = \sum_{k=0}^{n-1} f(t_k) \cdot \mathbf{1}_{(t_k, t_{k+1}]}(t) $$
の形であることをいいます。このとき、伊藤積分を
$$ \int_0^T f(t) \, dW(t) = \sum_{k=0}^{n-1} f(t_k) \left(W(t_{k+1}) – W(t_k)\right) $$
と定義します。
重要な点: 被積分関数 $f(t_k)$ は区間の 左端 で評価されます。これが 伊藤積分の左端評価規則 であり、ストラトノヴィッチ積分（中点評価）との本質的な違いです。左端で評価することにより、$f(t_k)$ は増分 $W(t_{k+1}) – W(t_k)$ と独立になります（$f(t_k)$ は時刻 $t_k$ までの情報のみに依存するため）。
ステップ2: 一般の過程への拡張
一般の適合過程 $f(t)$ に対しては、$E\left[\int_0^T f(t)^2 \, dt\right] < \infty$ の条件のもと、単純過程で近似して $L^2$ 極限として定義します。
伊藤積分の基本性質
性質1: 期待値が0（マルチンゲール性）
$$ E\left[\int_0^T f(t) \, dW(t)\right] = 0 $$
証明: 単純過程の場合、
$$ E\left[\sum_{k} f(t_k)(W(t_{k+1}) – W(t_k))\right] = \sum_{k} E[f(t_k)] \cdot E[W(t_{k+1}) – W(t_k)] = \sum_{k} E[f(t_k)] \cdot 0 = 0 $$
2つ目の等号は、$f(t_k)$ と $W(t_{k+1}) – W(t_k)$ の独立性（左端評価による）を用いました。$\square$
性質2: 伊藤の等長性（Ito isometry）
$$ E\left[\left(\int_0^T f(t) \, dW(t)\right)^2\right] = E\left[\int_0^T f(t)^2 \, dt\right] $$
証明: 単純過程の場合、
$$ \begin{align} E\left[\left(\sum_{k} f(t_k) \Delta W_k\right)^2\right] &= \sum_{k} \sum_{l} E[f(t_k) f(t_l) \Delta W_k \Delta W_l] \end{align} $$
$k \neq l$ のとき（$k < l$ と仮定）、$\Delta W_l$ は $f(t_k)$, $f(t_l)$, $\Delta W_k$ と独立なので、
$$ E[f(t_k) f(t_l) \Delta W_k \Delta W_l] = E[f(t_k) f(t_l) \Delta W_k] \cdot E[\Delta W_l] = 0 $$
$k = l$ のとき、$f(t_k)$ と $\Delta W_k$ の独立性から、
$$ E[f(t_k)^2 (\Delta W_k)^2] = E[f(t_k)^2] \cdot E[(\Delta W_k)^2] = E[f(t_k)^2] \cdot \Delta t_k $$
したがって、
$$ E\left[\left(\int_0^T f \, dW\right)^2\right] = \sum_{k} E[f(t_k)^2] \Delta t_k = E\left[\int_0^T f(t)^2 \, dt\right] $$
$\square$
伊藤の公式の導出
1次元の伊藤の公式
定理（伊藤の公式）: $X(t)$ が確率微分方程式
$$ dX = a(t, X) \, dt + b(t, X) \, dW $$
に従い、$f(t, x)$ が $t$ について1回、$x$ について2回連続微分可能であるとき、
$$ df(t, X) = \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial x} dX + \frac{1}{2} \frac{\partial^2 f}{\partial x^2} (dX)^2 $$
ここで $(dX)^2$ を乗法規則
$$ (dX)^2 = (a \, dt + b \, dW)^2 = a^2 (dt)^2 + 2ab \, dt \, dW + b^2 (dW)^2 = b^2 \, dt $$
で置き換えると、
$$ \boxed{df = \left(\frac{\partial f}{\partial t} + a \frac{\partial f}{\partial x} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial x^2}\right) dt + b \frac{\partial f}{\partial x} dW} $$
導出
ステップ1: テイラー展開
$f(t, x)$ の2次までのテイラー展開を書き下します。
$$ \begin{align} \Delta f &= f(t + \Delta t, X + \Delta X) – f(t, X) \\ &= \frac{\partial f}{\partial t} \Delta t + \frac{\partial f}{\partial x} \Delta X + \frac{1}{2} \frac{\partial^2 f}{\partial t^2} (\Delta t)^2 + \frac{\partial^2 f}{\partial t \partial x} \Delta t \Delta X + \frac{1}{2} \frac{\partial^2 f}{\partial x^2} (\Delta X)^2 + \cdots \end{align} $$
ステップ2: 各項のオーダーの評価
$\Delta X = a \Delta t + b \Delta W$ なので、
$$ \Delta X \sim O(\sqrt{\Delta t}) \quad (\because \Delta W \sim O(\sqrt{\Delta t})) $$
各項のオーダーを整理します。
| 項 | オーダー | 寄与 |
|---|---|---|
| $\frac{\partial f}{\partial t} \Delta t$ | $O(\Delta t)$ | 残る |
| $\frac{\partial f}{\partial x} \Delta X$ | $O(\sqrt{\Delta t})$ | 残る |
| $\frac{1}{2}\frac{\partial^2 f}{\partial x^2} (\Delta X)^2$ | $O(\Delta t)$ | 残る |
| $\frac{\partial^2 f}{\partial t \partial x} \Delta t \Delta X$ | $O(\Delta t^{3/2})$ | 消える |
| $\frac{1}{2}\frac{\partial^2 f}{\partial t^2} (\Delta t)^2$ | $O(\Delta t^2)$ | 消える |
通常の微積分との決定的な違いは、$(\Delta X)^2$ が $O(\Delta t)$ のオーダーであり消えないことです。
ステップ3: $(\Delta X)^2$ の計算
$$ \begin{align} (\Delta X)^2 &= (a \Delta t + b \Delta W)^2 \\ &= a^2 (\Delta t)^2 + 2ab \Delta t \Delta W + b^2 (\Delta W)^2 \end{align} $$
- $a^2 (\Delta t)^2 = O(\Delta t^2) \to 0$
- $2ab \Delta t \Delta W = O(\Delta t^{3/2}) \to 0$
- $b^2 (\Delta W)^2 \to b^2 \Delta t$ （2次変分の結果）
したがって、
$$ (\Delta X)^2 \to b^2 \, dt $$
ステップ4: まとめ
$O(\Delta t)$ 以上の項のみを残すと、
$$ \begin{align} df &= \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial x} (a \, dt + b \, dW) + \frac{1}{2} \frac{\partial^2 f}{\partial x^2} b^2 \, dt \\ &= \left(\frac{\partial f}{\partial t} + a \frac{\partial f}{\partial x} + \frac{1}{2} b^2 \frac{\partial^2 f}{\partial x^2}\right) dt + b \frac{\partial f}{\partial x} dW \end{align} $$
これが伊藤の公式です。$\square$
特殊ケース: $f(X)$ が $t$ に陽に依存しない場合
$f$ が $x$ のみの関数 $f(X)$ であり、$dX = a \, dt + b \, dW$ のとき、
$$ df(X) = f'(X) \, dX + \frac{1}{2} f”(X) b^2 \, dt $$
通常の微積分では $df = f'(X) \, dX$ なので、$\frac{1}{2} f”(X) b^2 \, dt$ の項が 伊藤の補正項（Ito correction term） です。
具体例: $f(X) = X^2$
$X = W(t)$（$a = 0$, $b = 1$）の場合、$f(x) = x^2$ とします。
$f'(x) = 2x$, $f”(x) = 2$ なので、
$$ \begin{align} d(W^2) &= 2W \, dW + \frac{1}{2} \cdot 2 \cdot 1^2 \, dt \\ &= 2W \, dW + dt \end{align} $$
積分すると、
$$ W(t)^2 = 2\int_0^t W(s) \, dW(s) + t $$
したがって、
$$ \int_0^t W(s) \, dW(s) = \frac{1}{2}(W(t)^2 – t) $$
通常の微積分では $\int x \, dx = x^2/2$ ですが、伊藤積分では $-t/2$ の補正項が現れます。これは伊藤積分の非自明な結果です。
幾何ブラウン運動の解の導出
確率微分方程式
幾何ブラウン運動（GBM）の確率微分方程式は、
$$ dS = \mu S \, dt + \sigma S \, dW $$
初期条件 $S(0) = S_0 > 0$。$\mu$ はドリフト（期待リターン率）、$\sigma$ はボラティリティです。
伊藤の公式による解法
$f(S) = \ln S$ とおきます。$f'(S) = 1/S$, $f”(S) = -1/S^2$ です。
$dS = \mu S \, dt + \sigma S \, dW$ の場合、$a = \mu S$, $b = \sigma S$ なので、伊藤の公式を適用すると、
$$ \begin{align} d(\ln S) &= \frac{1}{S} dS + \frac{1}{2}\left(-\frac{1}{S^2}\right)(\sigma S)^2 \, dt \\ &= \frac{1}{S}(\mu S \, dt + \sigma S \, dW) – \frac{\sigma^2}{2} dt \\ &= \mu \, dt + \sigma \, dW – \frac{\sigma^2}{2} dt \\ &= \left(\mu – \frac{\sigma^2}{2}\right) dt + \sigma \, dW \end{align} $$
この式は 通常の微分方程式（$\ln S$ に対する）なので、直接積分できます。
$$ \ln S(t) – \ln S(0) = \left(\mu – \frac{\sigma^2}{2}\right) t + \sigma W(t) $$
指数をとると、
$$ \boxed{S(t) = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right) t + \sigma W(t)\right]} $$
$-\sigma^2/2$ の項は 伊藤の補正 から来ています。もし通常の微積分の連鎖律を（誤って）適用すると、$d(\ln S) = dS/S = \mu \, dt + \sigma \, dW$ となり、$\sigma^2/2$ の項が欠落してしまいます。
オルンシュタイン=ウーレンベック過程
定義と物理的意味
オルンシュタイン=ウーレンベック（Ornstein-Uhlenbeck, OU）過程 は、次の確率微分方程式に従います。
$$ dX = \theta(\mu – X) \, dt + \sigma \, dW $$
ここで、
- $\theta > 0$: 平均回帰の速さ（回帰速度パラメータ）
- $\mu$: 長期平均
- $\sigma > 0$: 拡散係数
$X > \mu$ のときはドリフト $\theta(\mu – X) < 0$（減少方向）、$X < \mu$ のときはドリフト $\theta(\mu - X) > 0$（増加方向）であり、$X$ は常に $\mu$ に引き戻されます。この 平均回帰性（mean-reverting） がOU過程の最大の特徴です。
解の導出
$\mu = 0$ の場合を考えます（$\mu \neq 0$ の場合は $Y = X – \mu$ と変換すれば帰着）。
$$ dX = -\theta X \, dt + \sigma \, dW $$
積分因子 $e^{\theta t}$ を用います。$f(t, X) = e^{\theta t} X$ とおいて伊藤の公式を適用します。
$$ \frac{\partial f}{\partial t} = \theta e^{\theta t} X, \quad \frac{\partial f}{\partial x} = e^{\theta t}, \quad \frac{\partial^2 f}{\partial x^2} = 0 $$
$a = -\theta X$, $b = \sigma$ を代入すると、
$$ \begin{align} d(e^{\theta t} X) &= \left(\theta e^{\theta t} X + e^{\theta t}(-\theta X) + \frac{1}{2} \cdot 0 \cdot \sigma^2\right) dt + e^{\theta t} \sigma \, dW \\ &= \left(\theta e^{\theta t} X – \theta e^{\theta t} X\right) dt + \sigma e^{\theta t} \, dW \\ &= \sigma e^{\theta t} \, dW \end{align} $$
$dt$ の項が打ち消し合いました。両辺を $[0, t]$ で積分すると、
$$ e^{\theta t} X(t) – X(0) = \sigma \int_0^t e^{\theta s} \, dW(s) $$
$$ \boxed{X(t) = X(0) e^{-\theta t} + \sigma \int_0^t e^{-\theta(t-s)} \, dW(s)} $$
OU過程の統計的性質
平均:
$$ E[X(t)] = X(0) e^{-\theta t} + \sigma \underbrace{E\left[\int_0^t e^{-\theta(t-s)} \, dW(s)\right]}_{= 0 \text{ (伊藤積分の性質)}} = X(0) e^{-\theta t} $$
$t \to \infty$ で $E[X(t)] \to 0$ です（$\mu = 0$ の場合）。
分散: 伊藤の等長性を用いると、
$$ \begin{align} \text{Var}(X(t)) &= \sigma^2 E\left[\left(\int_0^t e^{-\theta(t-s)} \, dW(s)\right)^2\right] \\ &= \sigma^2 \int_0^t e^{-2\theta(t-s)} \, ds \quad (\text{伊藤の等長性}) \\ &= \sigma^2 \left[-\frac{1}{2\theta} e^{-2\theta(t-s)}\right]_{s=0}^{s=t} \\ &= \frac{\sigma^2}{2\theta}\left(1 – e^{-2\theta t}\right) \end{align} $$
$t \to \infty$ で、
$$ \text{Var}(X(t)) \to \frac{\sigma^2}{2\theta} $$
これがOU過程の 定常分散 です。$\theta$ が大きい（平均回帰が強い）ほど分散は小さく、$\sigma$ が大きい（ノイズが強い）ほど分散は大きくなります。
一般の場合（$\mu \neq 0$）
$$ X(t) = \mu + (X(0) – \mu) e^{-\theta t} + \sigma \int_0^t e^{-\theta(t-s)} \, dW(s) $$
定常分布は $N(\mu, \sigma^2/(2\theta))$ です。
ランジュバン方程式
OU過程は物理学では ランジュバン方程式 として知られています。ブラウン粒子の速度 $v(t)$ に対して、
$$ m \frac{dv}{dt} = -\gamma v + \sigma \xi(t) $$
ここで $\gamma$ は摩擦係数、$\xi(t)$ はホワイトノイズ（形式的に $\xi(t) = dW/dt$）です。
$m = 1$ として確率微分方程式の形に書くと、
$$ dv = -\gamma v \, dt + \sigma \, dW $$
これは $\theta = \gamma$, $\mu = 0$ のOU過程そのものです。摩擦によるエネルギー散逸と熱揺らぎによるエネルギー注入が釣り合って定常状態が形成されます。この釣り合いの条件は 揺動散逸定理（fluctuation-dissipation theorem） として知られ、
$$ \frac{\sigma^2}{2\gamma} = \frac{k_B T}{m} $$
すなわち $\sigma^2 = 2\gamma k_B T / m$ が成り立ちます。ここで $k_B$ はボルツマン定数、$T$ は絶対温度です。
Pythonでの実装
伊藤の公式の数値検証
$d(W^2) = 2W \, dW + dt$ を数値的に検証します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
T = 1.0
n_points = 10000
dt = T / n_points
t = np.linspace(0, T, n_points + 1)
# ブラウン運動の生成
dW = np.random.normal(0, np.sqrt(dt), size=n_points)
W = np.cumsum(dW)
W = np.insert(W, 0, 0)
# f(W) = W^2 の増分を3通りで計算
# (1) 直接計算: W(t)^2
f_direct = W**2
# (2) 通常の微積分（誤り）: ∫ 2W dW
# Σ 2W_k (W_{k+1} - W_k)
f_wrong = np.zeros(n_points + 1)
for k in range(n_points):
f_wrong[k + 1] = f_wrong[k] + 2 * W[k] * dW[k]
# (3) 伊藤の公式: ∫ 2W dW + ∫ dt = 2∫ W dW + t
f_ito = np.zeros(n_points + 1)
for k in range(n_points):
f_ito[k + 1] = f_ito[k] + 2 * W[k] * dW[k] + dt
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 比較プロット
axes[0].plot(t, f_direct, 'b-', linewidth=1.5, label='$W(t)^2$ (exact)')
axes[0].plot(t, f_ito, 'r--', linewidth=1.5, label="Ito: $\\int 2W\\,dW + t$")
axes[0].plot(t, f_wrong, 'g:', linewidth=1.5, label="Wrong: $\\int 2W\\,dW$")
axes[0].set_title("Verification of Ito's Formula: $d(W^2) = 2W\\,dW + dt$")
axes[0].set_xlabel("t")
axes[0].set_ylabel("Value")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 誤差プロット
error_ito = np.abs(f_direct - f_ito)
error_wrong = np.abs(f_direct - f_wrong)
axes[1].semilogy(t[1:], error_ito[1:], 'r-', linewidth=1, alpha=0.7, label='Ito error')
axes[1].semilogy(t[1:], error_wrong[1:], 'g-', linewidth=1, alpha=0.7, label='Wrong error')
axes[1].set_title("Error Comparison")
axes[1].set_xlabel("t")
axes[1].set_ylabel("Absolute error")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ito_formula_verification.png", dpi=150, bbox_inches="tight")
plt.show()
print("t=1 での値:")
print(f" W(1)^2 = {f_direct[-1]:.6f}")
print(f" 伊藤の公式 = {f_ito[-1]:.6f}")
print(f" 通常の微積分（誤り） = {f_wrong[-1]:.6f}")
print(f" 伊藤の誤差 = {error_ito[-1]:.6f}")
print(f" 通常の誤差 = {error_wrong[-1]:.6f}")
伊藤の公式による計算は $W(t)^2$ をほぼ正確に再現しますが、通常の微積分を適用した場合は $t$ だけずれることが確認できます。
オイラー=丸山法によるSDE数値解法
確率微分方程式
$$ dX = a(t, X) \, dt + b(t, X) \, dW $$
の数値解法として最も基本的な オイラー=丸山法（Euler-Maruyama method） を実装します。
$$ X_{n+1} = X_n + a(t_n, X_n) \Delta t + b(t_n, X_n) \Delta W_n $$
ここで $\Delta W_n = W(t_{n+1}) – W(t_n) \sim N(0, \Delta t)$ です。
import numpy as np
import matplotlib.pyplot as plt
def euler_maruyama(a, b, X0, T, n_steps, n_paths=1):
"""
オイラー=丸山法によるSDEの数値解法
dX = a(t, X)dt + b(t, X)dW
Parameters
----------
a : callable, ドリフト関数 a(t, X)
b : callable, 拡散関数 b(t, X)
X0 : float, 初期値
T : float, 終了時刻
n_steps : int, 時間刻み数
n_paths : int, サンプルパス数
Returns
-------
t : array, 時刻の配列
X : array, 解の配列 (n_paths, n_steps+1)
"""
dt = T / n_steps
t = np.linspace(0, T, n_steps + 1)
X = np.zeros((n_paths, n_steps + 1))
X[:, 0] = X0
for i in range(n_steps):
dW = np.random.normal(0, np.sqrt(dt), size=n_paths)
X[:, i + 1] = X[:, i] + a(t[i], X[:, i]) * dt + b(t[i], X[:, i]) * dW
return t, X
np.random.seed(42)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# --- (1) 幾何ブラウン運動 ---
mu_gbm = 0.1
sigma_gbm = 0.3
S0 = 100
a_gbm = lambda t, X: mu_gbm * X
b_gbm = lambda t, X: sigma_gbm * X
t_gbm, S_gbm = euler_maruyama(a_gbm, b_gbm, S0, T=2.0, n_steps=1000, n_paths=20)
# 解析解との比較（1パスで）
np.random.seed(123)
t_exact, S_exact_em = euler_maruyama(a_gbm, b_gbm, S0, T=2.0, n_steps=1000, n_paths=1)
for i in range(20):
axes[0, 0].plot(t_gbm, S_gbm[i], alpha=0.5, linewidth=0.5)
E_S = S0 * np.exp(mu_gbm * t_gbm)
axes[0, 0].plot(t_gbm, E_S, 'r-', linewidth=2, label='E[S(t)]')
axes[0, 0].set_title("GBM: $dS = \\mu S \\, dt + \\sigma S \\, dW$")
axes[0, 0].set_xlabel("t")
axes[0, 0].set_ylabel("S(t)")
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# --- (2) OU過程 ---
theta_ou = 5.0
mu_ou = 2.0
sigma_ou = 1.0
a_ou = lambda t, X: theta_ou * (mu_ou - X)
b_ou = lambda t, X: sigma_ou * np.ones_like(X)
t_ou, X_ou = euler_maruyama(a_ou, b_ou, X0=0.0, T=3.0, n_steps=3000, n_paths=20)
for i in range(20):
axes[0, 1].plot(t_ou, X_ou[i], alpha=0.5, linewidth=0.5)
axes[0, 1].axhline(y=mu_ou, color='r', linestyle='--', linewidth=2,
label=f'μ = {mu_ou}')
# ±2σの定常分散帯
sigma_stat = sigma_ou / np.sqrt(2 * theta_ou)
axes[0, 1].fill_between(t_ou, mu_ou - 2*sigma_stat, mu_ou + 2*sigma_stat,
alpha=0.15, color='red', label=f'±2σ_stat')
axes[0, 1].set_title("OU Process: $dX = \\theta(\\mu - X)\\,dt + \\sigma\\,dW$")
axes[0, 1].set_xlabel("t")
axes[0, 1].set_ylabel("X(t)")
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# --- (3) ランジュバン方程式（速度過程） ---
gamma_lang = 2.0
sigma_lang = 1.0
a_lang = lambda t, X: -gamma_lang * X
b_lang = lambda t, X: sigma_lang * np.ones_like(X)
t_lang, V_lang = euler_maruyama(a_lang, b_lang, X0=3.0, T=5.0, n_steps=5000, n_paths=10)
for i in range(10):
axes[1, 0].plot(t_lang, V_lang[i], alpha=0.6, linewidth=0.5)
axes[1, 0].axhline(y=0, color='r', linestyle='--', linewidth=1.5)
axes[1, 0].set_title("Langevin Equation: $dv = -\\gamma v \\, dt + \\sigma \\, dW$")
axes[1, 0].set_xlabel("t")
axes[1, 0].set_ylabel("v(t)")
axes[1, 0].grid(True, alpha=0.3)
# --- (4) 収束次数の検証（GBMの強収束） ---
np.random.seed(42)
T_conv = 1.0
n_fine = 2**16 # 非常に細かい刻み（参照解）
dt_fine = T_conv / n_fine
# 参照解の生成
dW_fine = np.random.normal(0, np.sqrt(dt_fine), size=n_fine)
W_fine = np.cumsum(dW_fine)
W_fine = np.insert(W_fine, 0, 0)
S_ref = S0 * np.exp((mu_gbm - sigma_gbm**2/2) * T_conv + sigma_gbm * W_fine[-1])
# 異なる刻み幅でオイラー=丸山法を適用
n_steps_list = [2**k for k in range(4, 14)]
errors = []
for n_step in n_steps_list:
ratio = n_fine // n_step
dt_em = T_conv / n_step
S_em = S0
for i in range(n_step):
# 対応する微小増分を集約
dW_sum = np.sum(dW_fine[i*ratio:(i+1)*ratio])
S_em = S_em + mu_gbm * S_em * dt_em + sigma_gbm * S_em * dW_sum
errors.append(np.abs(S_em - S_ref))
axes[1, 1].loglog(n_steps_list, errors, 'bo-', linewidth=1.5, markersize=5,
label='EM error')
# 理論的な1/2次収束の参照線
ref_line = errors[0] * (np.array(n_steps_list[0]) / np.array(n_steps_list))**0.5
axes[1, 1].loglog(n_steps_list, ref_line, 'r--', linewidth=1.5,
label='$O(\\Delta t^{1/2})$')
axes[1, 1].set_title("Strong Convergence of Euler-Maruyama")
axes[1, 1].set_xlabel("Number of steps")
axes[1, 1].set_ylabel("$|S_{EM}(T) - S_{exact}(T)|$")
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("sde_euler_maruyama.png", dpi=150, bbox_inches="tight")
plt.show()
オイラー=丸山法は 強収束次数 1/2（$O(\sqrt{\Delta t})$）を持ちます。これは通常のオイラー法（次数1）より遅い収束ですが、確率項の影響のため仕方がありません。より高次の方法（ミルシュタイン法、次数1）も存在しますが、実装が複雑になります。
OU過程の定常分布への収束
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(42)
# OU過程のパラメータ
theta = 5.0
mu = 0.0
sigma = 1.0
# 定常分布: N(mu, sigma^2 / (2*theta))
var_stat = sigma**2 / (2 * theta)
std_stat = np.sqrt(var_stat)
# 多数のパスをシミュレーション
n_paths = 5000
T = 3.0
n_steps = 3000
dt = T / n_steps
X = np.zeros((n_paths, n_steps + 1))
X[:, 0] = 5.0 # 定常分布から遠い初期値
for i in range(n_steps):
dW = np.random.normal(0, np.sqrt(dt), size=n_paths)
X[:, i + 1] = X[:, i] + theta * (mu - X[:, i]) * dt + sigma * dW
# 異なる時刻でのX(t)の分布を可視化
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
times_to_plot = [0.0, 0.1, 0.3, 0.5, 1.0, 3.0]
for idx, t_plot in enumerate(times_to_plot):
ax = axes[idx // 3][idx % 3]
step_idx = int(t_plot / dt)
if step_idx >= n_steps + 1:
step_idx = n_steps
data = X[:, step_idx]
ax.hist(data, bins=50, density=True, alpha=0.7, color='steelblue',
edgecolor='black', linewidth=0.5, label='Simulation')
# 理論的分布
mean_theory = 5.0 * np.exp(-theta * t_plot) # 初期値5.0から
var_theory = var_stat * (1 - np.exp(-2 * theta * t_plot))
std_theory = np.sqrt(max(var_theory, 1e-10))
x_range = np.linspace(mean_theory - 4*std_theory, mean_theory + 4*std_theory, 200)
if std_theory > 0.01:
pdf_theory = norm.pdf(x_range, loc=mean_theory, scale=std_theory)
ax.plot(x_range, pdf_theory, 'r-', linewidth=2, label='Theory')
# 定常分布
x_stat = np.linspace(mu - 4*std_stat, mu + 4*std_stat, 200)
pdf_stat = norm.pdf(x_stat, loc=mu, scale=std_stat)
ax.plot(x_stat, pdf_stat, 'g--', linewidth=1.5, label='Stationary')
ax.set_title(f"t = {t_plot}")
ax.set_xlabel("X")
ax.set_ylabel("Density")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
plt.suptitle("OU Process: Convergence to Stationary Distribution ($X_0=5.0$)",
fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig("ou_process_convergence.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"定常分布: N({mu}, {var_stat:.4f})")
print(f"t=3.0 での標本統計量:")
print(f" 標本平均: {np.mean(X[:, -1]):.4f} (理論値: {mu})")
print(f" 標本分散: {np.var(X[:, -1]):.4f} (理論値: {var_stat:.4f})")
まとめ
本記事では、確率微分方程式と伊藤の公式について導出から応用まで解説しました。
- ブラウン運動の2次変分 $(dW)^2 = dt \neq 0$ のため、通常の微積分の連鎖律は修正が必要
- 伊藤積分は左端評価により定義され、期待値が0（マルチンゲール性）と伊藤の等長性を持つ
- 伊藤の公式は $df = \frac{\partial f}{\partial t}dt + \frac{\partial f}{\partial x}dX + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}(dX)^2$ であり、2次の項が生き残ることが本質
- 幾何ブラウン運動 $dS = \mu S dt + \sigma S dW$ の解は $S(t) = S_0 \exp((\mu – \sigma^2/2)t + \sigma W(t))$ であり、伊藤の補正 $-\sigma^2/2$ が現れる
- OU過程は平均回帰性を持ち、定常分布 $N(\mu, \sigma^2/(2\theta))$ に収束する
- オイラー=丸山法はSDE数値解法の基本であり、強収束次数 $1/2$ を持つ
次のステップとして、以下の記事も参考にしてください。