確率微分方程式（SDE）は、ランダムな揺らぎを含む系の時間発展を記述する最も一般的な枠組みです。しかし、SDEを解こうとすると、通常の微積分の連鎖律（合成関数の微分法）がそのまま使えないという壁にぶつかります。たとえば $X_t = B_t^2$ が満たすSDEを求めたいとき、通常の微積分なら $d(x^2) = 2x\,dx$ ですが、ブラウン運動の二次変分 $(dB_t)^2 = dt$ のために追加の項が現れます。この「確率解析版の連鎖律」を与えるのが伊藤の公式（Ito’s formula, 伊藤の補題とも呼ばれる）です。
伊藤の公式は確率解析の中で最も広く使われる公式であり、その応用は多岐にわたります。
- 金融工学: ブラック・ショールズ方程式の導出は、株価の幾何ブラウン運動に伊藤の公式を適用することで行われます
- 物理学: ランジュバン方程式からフォッカー・プランク方程式への変換に伊藤の公式が使われます
- 制御工学: 確率的制御問題のHJB方程式（ハミルトン・ヤコビ・ベルマン方程式）の導出に不可欠です
本記事では、伊藤の公式を直感的に理解した上で厳密に導出し、代表的なSDEの求解テクニックをPython実装とともに解説します。
本記事の内容
- 伊藤の公式の直感的理解と厳密な導出
- 幾何ブラウン運動（GBM）の解法と応用
- オルンシュタイン・ウーレンベック過程の解法
- 多次元への拡張と積の公式
- SDEの求解テクニックの体系的な整理
- Pythonによる数値解との比較検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 伊藤積分の定義と性質 — 伊藤積分の構成と等長公式
- ブラウン運動 — ウィーナー過程の基本性質
- 確率過程の基礎 — 確率過程の基本的な枠組み
伊藤の公式 — 確率解析の連鎖律
なぜ通常の連鎖律が壊れるのか
通常の微積分では、$f(x)$ が滑らかで $x(t)$ が微分可能な関数であれば、合成関数の微分は連鎖律（chain rule）で計算できます。
$$ df(x) = f'(x)\,dx $$
たとえば $f(x) = x^2$ ならば $d(x^2) = 2x\,dx$ です。これをテイラー展開の観点で見ると、
$$ f(x + \Delta x) – f(x) = f'(x)\Delta x + \frac{1}{2}f”(x)(\Delta x)^2 + \cdots $$
通常の関数では $(\Delta x)^2$ は $\Delta x$ よりもずっと小さいため（$\Delta x \to 0$ で2次以上の項は無視できるため）、$df \approx f'(x)\,dx$ が正確です。
ところが、$x = B_t$ がブラウン運動の場合は事情が異なります。伊藤積分の定義と性質で見たように、ブラウン運動の増分 $\Delta B_t$ は $\sqrt{\Delta t}$ のオーダーを持つため、
$$ (\Delta B_t)^2 \sim \Delta t $$
であり、2次の項 $\frac{1}{2}f”(B_t)(\Delta B_t)^2 \sim \frac{1}{2}f”(B_t)\Delta t$ は1次の項と同じオーダーになります。したがって、2次の項を無視できず、連鎖律に修正が必要になるのです。
この「2次の項が生き残る」という現象こそが、確率解析を通常の解析学から区別する最も本質的な特徴です。
伊藤の公式の主張
伊藤の公式（1次元版）: $X_t$ が伊藤過程、すなわち
$$ \begin{equation} dX_t = \mu(t, X_t)\,dt + \sigma(t, X_t)\,dB_t \end{equation} $$
の形のSDEを満たし、$f(t, x)$ が $t$ に関して1回、$x$ に関して2回連続微分可能な関数であるとき、過程 $Y_t = f(t, X_t)$ は次のSDEを満たします。
$$ \begin{equation} df(t, X_t) = \frac{\partial f}{\partial t}\,dt + \frac{\partial f}{\partial x}\,dX_t + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}(dX_t)^2 \end{equation} $$
ここで $(dX_t)^2$ は次の伊藤の掛け算規則（Ito multiplication rules）に従って計算します。
$$ (dB_t)^2 = dt, \quad dt \cdot dB_t = 0, \quad (dt)^2 = 0 $$
これらの規則を適用すると、$(dX_t)^2 = \sigma^2\,dt$ となり、伊藤の公式は次の形に展開されます。
$$ \begin{equation} df = \left(\frac{\partial f}{\partial t} + \mu\frac{\partial f}{\partial x} + \frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial x^2}\right)dt + \sigma\frac{\partial f}{\partial x}\,dB_t \end{equation} $$
通常の連鎖律との違いは、$\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial x^2}\,dt$ という伊藤補正項が追加されていることです。この項はブラウン運動の二次変分に起因し、$\sigma = 0$（ノイズなし）のときは消えて通常の連鎖律に帰着します。
伊藤の公式の導出
伊藤の公式を導出するために、$f(t, X_t)$ のテイラー展開を2次の項まで考えます。
$$ \Delta f = f(t + \Delta t, X_t + \Delta X) – f(t, X_t) $$
2次までテイラー展開すると、
$$ \Delta f = \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 \cdot \Delta X + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}(\Delta X)^2 + \cdots $$
各項のオーダーを伊藤の掛け算規則に基づいて評価します。$\Delta X = \mu\Delta t + \sigma\Delta B$ ですから、
$\Delta X$ のオーダーを調べると、$\sigma\Delta B$ は $\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$: $(\Delta X)^2 = \sigma^2(\Delta B)^2 + 2\mu\sigma\Delta t\Delta B + \mu^2(\Delta t)^2$
$(\Delta X)^2$ の各項を評価すると、
- $\sigma^2(\Delta B)^2 \to \sigma^2\,dt$（二次変分）: $O(\Delta t)$
- $2\mu\sigma\Delta t\Delta B$: $O(\Delta t^{3/2})$ — 無視できる
- $\mu^2(\Delta t)^2$: $O(\Delta t^2)$ — 無視できる
したがって $(\Delta X)^2 \to \sigma^2\,dt$ であり、3次以上の項も $O(\Delta t^{3/2})$ 以上なので無視できます。
残る項をまとめると、
$$ df = \frac{\partial f}{\partial t}\,dt + \frac{\partial f}{\partial x}(\mu\,dt + \sigma\,dB_t) + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}\sigma^2\,dt $$
これを整理すると、
$$ df = \left(\frac{\partial f}{\partial t} + \mu\frac{\partial f}{\partial x} + \frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial x^2}\right)dt + \sigma\frac{\partial f}{\partial x}\,dB_t $$
が得られ、伊藤の公式が導出されました。$\square$
この導出の核心は、$(\Delta B)^2 \to dt$ という二次変分の性質にあります。通常のテイラー展開では2次の項は微小量として無視しますが、ブラウン運動では2次の項が1次の項と同じオーダーになるため、無視できないのです。
伊藤の公式を手にしたことで、様々なSDEを解く道具が揃いました。まず最も重要な応用例である幾何ブラウン運動から見ていきましょう。
幾何ブラウン運動の解法
問題の設定
幾何ブラウン運動（Geometric Brownian Motion, GBM）は、金融工学で株価のモデルとして最もよく使われるSDEです。
$$ \begin{equation} dS_t = \mu S_t\,dt + \sigma S_t\,dB_t, \quad S_0 > 0 \end{equation} $$
ここで $\mu$ はドリフト率（平均的な収益率）、$\sigma$ はボラティリティ（変動の激しさ）です。このSDEの特徴は、ドリフトと拡散の両方が $S_t$ に比例しているため、株価が大きいときは変動の絶対値も大きくなるという性質を捉えていることです。
このSDEを通常の分離変数法で解こうとすると、$\frac{dS}{S} = \mu\,dt + \sigma\,dB_t$ の左辺を $d(\ln S)$ と書きたくなります。しかし、伊藤の公式を使わずに $d(\ln S) = dS/S$ と置いてしまうと、伊藤補正項を見落としてしまいます。
伊藤の公式による解法
$f(x) = \ln x$ と置き、伊藤の公式を $X_t = S_t$ に適用します。$f'(x) = 1/x$, $f”(x) = -1/x^2$ であるから、
$$ d(\ln S_t) = \frac{1}{S_t}\,dS_t + \frac{1}{2}\left(-\frac{1}{S_t^2}\right)(dS_t)^2 $$
$(dS_t)^2$ を計算します。$dS_t = \mu S_t\,dt + \sigma S_t\,dB_t$ より、
$$ (dS_t)^2 = \sigma^2 S_t^2\,(dB_t)^2 = \sigma^2 S_t^2\,dt $$
$dt \cdot dB_t = 0$、$(dt)^2 = 0$ の項は消えることを使いました。
これを代入すると、
$$ d(\ln S_t) = \frac{\mu S_t\,dt + \sigma S_t\,dB_t}{S_t} – \frac{\sigma^2 S_t^2\,dt}{2S_t^2} $$
分子分母の $S_t$ を約分すると、
$$ \begin{equation} d(\ln S_t) = \left(\mu – \frac{\sigma^2}{2}\right)dt + \sigma\,dB_t \end{equation} $$
右辺は $S_t$ に依存しない（定数とブラウン運動の増分だけの）式になりました。これは解けます。$t = 0$ から $t = T$ まで積分すると、
$$ \ln S_T – \ln S_0 = \left(\mu – \frac{\sigma^2}{2}\right)T + \sigma B_T $$
指数をとれば、幾何ブラウン運動の解が得られます。
$$ \begin{equation} S_T = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right)T + \sigma B_T\right] \end{equation} $$
伊藤補正項の物理的意味
解の中に現れる $-\sigma^2/2$ という項が伊藤補正項です。もし通常の微積分の連鎖律を使って $d(\ln S) = dS/S$ としていたら、$\ln S_T = \ln S_0 + \mu T + \sigma B_T$ となり、$-\sigma^2 T/2$ の項が欠落していたでしょう。
この補正項には重要な物理的意味があります。$\ln S_T$ は正規分布 $N(\ln S_0 + (\mu – \sigma^2/2)T, \sigma^2 T)$ に従うため、$\ln S$ の期待値は $\ln S_0 + (\mu – \sigma^2/2)T$ です。一方、$E[S_T] = S_0 e^{\mu T}$（$S_T$ そのものの期待値）です。
$E[\ln S_T] \neq \ln E[S_T]$ であることに注意してください。対数の期待値と期待値の対数のずれが $-\sigma^2 T/2$ であり、これはイェンセンの不等式（$\ln$ が凹関数であること）の帰結です。つまり、ボラティリティが大きいほど対数リターンの平均は小さくなるのです。金融では、この効果をボラティリティ・ドラッグ（volatility drag）と呼びます。
では次に、物理学で頻繁に現れるもう一つの重要なSDEであるオルンシュタイン・ウーレンベック過程を解いてみましょう。
オルンシュタイン・ウーレンベック過程の解法
問題の設定
オルンシュタイン・ウーレンベック過程（Ornstein-Uhlenbeck process, OU過程）は、平均回帰性を持つ確率過程です。
$$ \begin{equation} dX_t = -\theta(X_t – \mu)\,dt + \sigma\,dB_t, \quad X_0 = x_0 \end{equation} $$
ここで $\theta > 0$ は回帰速度、$\mu$ は長期平均、$\sigma$ は拡散定数です。このSDEは「$X_t$ が $\mu$ から離れると、$-\theta(X_t – \mu)$ の力で引き戻される」一方、ブラウン運動によるランダムな揺らぎが加わるという物理を記述しています。
OU過程は、ブラウン運動に「バネの復元力」を加えたものと理解できます。具体的な応用例としては、
- 物理学: 粘性流体中の粒子の速度（ランジュバン方程式）
- 金融工学: 金利のバシチェクモデル、ペアトレーディングのスプレッド
- 生物学: 形質の進化モデル
があります。
積分因子法による解法
OU過程のSDEは線形ですので、常微分方程式の積分因子法（integrating factor method）のアナログが使えます。
$Y_t = e^{\theta t}X_t$ と変数変換します。伊藤の公式を $f(t, x) = e^{\theta t} x$ に適用すると、$\partial f/\partial t = \theta e^{\theta t} x$, $\partial f/\partial x = e^{\theta t}$, $\partial^2 f/\partial x^2 = 0$ ですから、
$$ dY_t = \theta e^{\theta t} X_t\,dt + e^{\theta t}\,dX_t + 0 $$
ここで $dX_t = -\theta(X_t – \mu)\,dt + \sigma\,dB_t$ を代入すると、
$$ dY_t = \theta e^{\theta t} X_t\,dt + e^{\theta t}[-\theta(X_t – \mu)\,dt + \sigma\,dB_t] $$
各項を展開します。
$$ dY_t = \theta e^{\theta t} X_t\,dt – \theta e^{\theta t} X_t\,dt + \theta\mu e^{\theta t}\,dt + \sigma e^{\theta t}\,dB_t $$
最初の2項がキャンセルし、結果は簡潔になります。
$$ \begin{equation} dY_t = \theta\mu e^{\theta t}\,dt + \sigma e^{\theta t}\,dB_t \end{equation} $$
右辺は $Y_t$ に依存しないので、直接積分できます。$t = 0$ から $t = T$ まで積分すると、
$$ Y_T – Y_0 = \theta\mu\int_0^T e^{\theta s}\,ds + \sigma\int_0^T e^{\theta s}\,dB_s $$
$Y_T = e^{\theta T} X_T$, $Y_0 = x_0$ を代入し、$X_T$ について解くと、
$$ X_T = x_0 e^{-\theta T} + \mu(1 – e^{-\theta T}) + \sigma\int_0^T e^{-\theta(T-s)}\,dB_s $$
第1の通常積分は $\int_0^T e^{\theta s}\,ds = (e^{\theta T} – 1)/\theta$ を使い、$\theta\mu \cdot (e^{\theta T} – 1)/\theta = \mu(e^{\theta T} – 1)$ と計算し、$e^{-\theta T}$ を掛けて $\mu(1 – e^{-\theta T})$ を得ました。
最終的な解は次のようになります。
$$ \begin{equation} X_T = \mu + (x_0 – \mu)e^{-\theta T} + \sigma\int_0^T e^{-\theta(T-s)}\,dB_s \end{equation} $$
OU過程の統計的性質
解の構造から、OU過程の統計的性質を読み取ることができます。
期待値: 伊藤積分の期待値がゼロであることから、
$$ E[X_T] = \mu + (x_0 – \mu)e^{-\theta T} $$
$T \to \infty$ で $E[X_T] \to \mu$ であり、初期値 $x_0$ によらず長期平均 $\mu$ に収束します。
分散: 伊藤の等長公式から、
$$ \text{Var}(X_T) = \sigma^2\int_0^T e^{-2\theta(T-s)}\,ds = \frac{\sigma^2}{2\theta}(1 – e^{-2\theta T}) $$
$T \to \infty$ で $\text{Var}(X_T) \to \sigma^2/(2\theta)$ であり、分散は有限値に収束します。これはブラウン運動（分散が $T$ に比例して発散する）とは本質的に異なる性質です。
定常分布: 十分長い時間が経過すると、OU過程は定常分布に収束します。
$$ X_\infty \sim N\left(\mu, \frac{\sigma^2}{2\theta}\right) $$
回帰速度 $\theta$ が大きいほど分散が小さくなり、過程は長期平均の近くに集中します。$\theta$ が小さいと揺らぎが大きくなり、ブラウン運動に近い振る舞いをします。
次に、より複雑なSDEを解くための一般的なテクニックを整理しましょう。
SDEの求解テクニック
テクニック1: 変数変換と伊藤の公式
SDEを解く最も一般的な方法は、適切な変数変換 $Y_t = f(t, X_t)$ によってSDEを簡単な形に変換し、伊藤の公式で $Y_t$ のSDEを導出するというものです。
変換の選び方にはいくつかのパターンがあります。
対数変換: $dX_t = \mu X_t\,dt + \sigma X_t\,dB_t$ の形（係数が $X_t$ に比例）のSDEに対して、$Y_t = \ln X_t$ と置くと線形SDEに帰着します（幾何ブラウン運動の解法で見た手法）。
積分因子: $dX_t = a(t)X_t\,dt + b(t)\,dB_t$ の形の線形SDEに対して、$Y_t = e^{-\int_0^t a(s)\,ds} X_t$ と置くと、$X_t$ への依存性が消えた形になります（OU過程の解法で見た手法）。
べき変換: ある種の非線形SDEでは $Y_t = X_t^p$ と置くことで解ける場合があります。
テクニック2: SDEの指数解（確率指数関数）
線形SDE
$$ dX_t = \mu(t)\,X_t\,dt + \sigma(t)\,X_t\,dB_t $$
の解は、確率指数関数（stochastic exponential, Doleans-Dade exponential）として体系的に求めることができます。
$$ \begin{equation} X_t = X_0 \exp\left[\int_0^t \left(\mu(s) – \frac{\sigma(s)^2}{2}\right)ds + \int_0^t \sigma(s)\,dB_s\right] \end{equation} $$
$\mu(s)$ と $\sigma(s)$ が定数の場合は幾何ブラウン運動の解そのものですが、時間依存の場合にも同じ形で解が得られます。
テクニック3: 自励系と非自励系
SDEの係数が時間 $t$ に依存しない場合（自励系）は、時間とは独立な変数変換が有効です。時間 $t$ に陽に依存する場合（非自励系）は、$f(t, x)$ の形の変換が必要になります。
具体例: CEV過程
以上のテクニックを応用して、少し複雑なSDEを解いてみましょう。定弾力性ボラティリティ過程（Constant Elasticity of Variance, CEV）は、
$$ dS_t = \mu S_t\,dt + \sigma S_t^\beta\,dB_t $$
の形をしています。$\beta = 1$ は幾何ブラウン運動に帰着しますが、一般の $\beta$ では閉じた解が得られません。しかし $\beta = 1/2$ の場合（CIR過程、コックス・インガソール・ロスモデル）は、
$$ dX_t = \kappa(\theta – X_t)\,dt + \sigma\sqrt{X_t}\,dB_t $$
として金利モデルで広く使われています。この場合、$Y_t = \sqrt{X_t}$ と変換して伊藤の公式を適用することで、解の性質を分析できます。
ここまでの求解テクニックをPythonで実際に確認してみましょう。
伊藤の公式の積の公式と多次元への拡張
伊藤の積の公式
二つの伊藤過程 $X_t$ と $Y_t$ の積 $Z_t = X_t Y_t$ のSDEを求める公式が、伊藤の積の公式（Ito product rule）です。
$$ \begin{equation} d(X_t Y_t) = X_t\,dY_t + Y_t\,dX_t + dX_t \cdot dY_t \end{equation} $$
通常のライプニッツの積の法則 $d(xy) = x\,dy + y\,dx$ に、交差二次変分 $dX_t \cdot dY_t$ が追加されている点が異なります。
例えば、$dX_t = \mu_1\,dt + \sigma_1\,dB_t^{(1)}$ と $dY_t = \mu_2\,dt + \sigma_2\,dB_t^{(2)}$ のとき、
$$ dX_t \cdot dY_t = \sigma_1 \sigma_2\,dB_t^{(1)} \cdot dB_t^{(2)} $$
$B^{(1)}$ と $B^{(2)}$ が相関 $\rho$ を持つ場合は $dB^{(1)} \cdot dB^{(2)} = \rho\,dt$ です。独立な場合は $\rho = 0$ なので交差項は消えます。
多次元の伊藤の公式
$n$ 次元の伊藤過程 $\bm{X}_t = (X_t^{(1)}, \dots, X_t^{(n)})$ と、$f(t, \bm{x})$ に対する多次元版の伊藤の公式は次のようになります。
$$ \begin{equation} df = \frac{\partial f}{\partial t}\,dt + \sum_{i=1}^n \frac{\partial f}{\partial x_i}\,dX_t^{(i)} + \frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n \frac{\partial^2 f}{\partial x_i \partial x_j}\,dX_t^{(i)} \cdot dX_t^{(j)} \end{equation} $$
ここで $dX_t^{(i)} \cdot dX_t^{(j)}$ は伊藤の掛け算規則に従います。
この多次元版は、ブラック・ショールズ方程式の多資産版の導出や、多次元確率制御問題で使われます。
多次元の理論は抽象的になりがちですが、1次元のテクニックの自然な拡張です。次のPythonシミュレーションでは、1次元の代表的なSDEの解を数値的に検証します。
Pythonによる数値解との比較検証
幾何ブラウン運動の厳密解と数値解の比較
まず、幾何ブラウン運動の厳密解と、オイラー・丸山法（Euler-Maruyama method）による数値解を比較します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# --- パラメータ ---
S0 = 100.0 # 初期値
mu = 0.08 # ドリフト
sigma = 0.3 # ボラティリティ
T = 2.0 # 時間
n = 1000 # 時間分割数
n_paths = 5 # パス数
dt = T / n
t = np.linspace(0, T, n + 1)
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# --- 左図: 厳密解のサンプルパス ---
ax = axes[0]
for i in range(n_paths):
dB = np.random.normal(0, np.sqrt(dt), n)
B = np.cumsum(dB)
B = np.insert(B, 0, 0)
# 厳密解
S_exact = S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * B)
ax.plot(t, S_exact, linewidth=1.2, label=f'Path {i+1}')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$S_t$', fontsize=12)
ax.set_title('GBM: Exact solution (sample paths)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# --- 右図: 厳密解 vs オイラー・丸山法 ---
ax = axes[1]
dB = np.random.normal(0, np.sqrt(dt), n)
B = np.cumsum(dB)
B = np.insert(B, 0, 0)
# 厳密解
S_exact = S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * B)
# オイラー・丸山法
S_em = np.zeros(n + 1)
S_em[0] = S0
for k in range(n):
S_em[k + 1] = S_em[k] + mu * S_em[k] * dt + sigma * S_em[k] * dB[k]
# 「伊藤補正なし」の間違った解
S_wrong = S0 * np.exp(mu * t + sigma * B)
ax.plot(t, S_exact, 'b-', linewidth=2, label='Exact (with Ito correction)')
ax.plot(t, S_em, 'r--', linewidth=1.5, label='Euler-Maruyama')
ax.plot(t, S_wrong, 'g:', linewidth=1.5, label='Without Ito correction (wrong)')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$S_t$', fontsize=12)
ax.set_title('Exact vs numerical vs wrong solution', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gbm_solutions.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、以下のことが確認できます。
-
左図: 幾何ブラウン運動の各パスは対数正規分布に従う振る舞いを示す。$S_t > 0$ が常に成り立ち、対数スケールでは通常のブラウン運動のように振る舞います。パスによってはドリフトを超える上昇を示し、他のパスでは下落していますが、これはボラティリティの効果です
-
右図: 厳密解（青の実線）とオイラー・丸山法（赤の破線）がよく一致している。時間分割数 $n = 1000$ では、数値解は厳密解に十分近い精度で追従しています。一方、伊藤補正を無視した「間違った」解（緑の点線）は厳密解から系統的にずれており、特に長時間で乖離が大きくなります。$-\sigma^2/2$ の補正項を無視すると、解が上方にバイアスされることが分かります
OU過程の厳密解と統計的性質の検証
次に、オルンシュタイン・ウーレンベック過程の厳密解と定常分布を検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(123)
# --- パラメータ ---
theta_ou = 2.0 # 回帰速度
mu_ou = 1.0 # 長期平均
sigma_ou = 0.5 # 拡散定数
x0 = 3.0 # 初期値
T = 5.0
n = 5000
n_paths = 10000
dt = T / n
t = np.linspace(0, T, n + 1)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# --- 左上: サンプルパスと理論的期待値 ---
ax = axes[0, 0]
for i in range(8):
dB = np.random.normal(0, np.sqrt(dt), n)
X = np.zeros(n + 1)
X[0] = x0
for k in range(n):
X[k + 1] = X[k] - theta_ou * (X[k] - mu_ou) * dt + sigma_ou * dB[k]
ax.plot(t, X, linewidth=0.7, alpha=0.7)
# 理論的な期待値
E_X = mu_ou + (x0 - mu_ou) * np.exp(-theta_ou * t)
ax.plot(t, E_X, 'k-', linewidth=2.5, label=r'$E[X_t]$ (theory)')
ax.axhline(mu_ou, color='gray', linestyle='--', linewidth=1, alpha=0.5)
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$X_t$', fontsize=12)
ax.set_title('OU process: sample paths', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# --- 右上: 分散の時間変化 ---
ax = axes[0, 1]
X_all = np.zeros((n_paths, n + 1))
X_all[:, 0] = x0
dB_all = np.random.normal(0, np.sqrt(dt), (n_paths, n))
for k in range(n):
X_all[:, k + 1] = (X_all[:, k] - theta_ou * (X_all[:, k] - mu_ou) * dt
+ sigma_ou * dB_all[:, k])
var_numerical = np.var(X_all, axis=0)
var_theory = (sigma_ou**2 / (2 * theta_ou)) * (1 - np.exp(-2 * theta_ou * t))
var_stationary = sigma_ou**2 / (2 * theta_ou)
ax.plot(t, var_numerical, 'b-', linewidth=1.5, label='Numerical')
ax.plot(t, var_theory, 'r--', linewidth=2, label='Theory')
ax.axhline(var_stationary, color='gray', linestyle=':', linewidth=1.5,
label=rf'Stationary: $\sigma^2/(2\theta) = {var_stationary:.4f}$')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('Var$(X_t)$', fontsize=12)
ax.set_title('Variance convergence', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# --- 左下: 定常分布のヒストグラム ---
ax = axes[1, 0]
X_final = X_all[:, -1]
ax.hist(X_final, bins=80, density=True, alpha=0.7, color='steelblue',
edgecolor='white', linewidth=0.5, label='Simulation')
x_range = np.linspace(X_final.min(), X_final.max(), 300)
ax.plot(x_range, norm.pdf(x_range, mu_ou, np.sqrt(var_stationary)),
'r-', linewidth=2, label=rf'$N(\mu, \sigma^2/2\theta)$')
ax.set_xlabel('$X_T$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'Stationary distribution at $T={T}$', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# --- 右下: 自己相関関数 ---
ax = axes[1, 1]
# 1つの長いパスから自己相関を推定
T_long = 100.0
n_long = 100000
dt_long = T_long / n_long
dB_long = np.random.normal(0, np.sqrt(dt_long), n_long)
X_long = np.zeros(n_long + 1)
X_long[0] = mu_ou # 定常状態から開始
for k in range(n_long):
X_long[k + 1] = X_long[k] - theta_ou * (X_long[k] - mu_ou) * dt_long + sigma_ou * dB_long[k]
# 自己相関の計算
max_lag = 2000
X_centered = X_long[n_long // 2:] - mu_ou # 定常状態部分を使用
autocorr = np.correlate(X_centered[:max_lag * 5], X_centered[:max_lag * 5], mode='full')
autocorr = autocorr[len(autocorr) // 2:]
autocorr = autocorr / autocorr[0]
lags = np.arange(max_lag) * dt_long
theory_acf = np.exp(-theta_ou * lags)
ax.plot(lags, autocorr[:max_lag], 'b-', linewidth=1.5, alpha=0.8, label='Simulation')
ax.plot(lags, theory_acf, 'r--', linewidth=2, label=r'$e^{-\theta \tau}$')
ax.set_xlabel(r'Lag $\tau$', fontsize=12)
ax.set_ylabel('Autocorrelation', fontsize=12)
ax.set_title('Autocorrelation function', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ou_process.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"定常分散 (理論): {var_stationary:.6f}")
print(f"定常分散 (数値): {np.var(X_final):.6f}")
print(f"定常平均 (理論): {mu_ou}")
print(f"定常平均 (数値): {np.mean(X_final):.6f}")
この4つのグラフから、OU過程の理論的性質がシミュレーションで正確に再現されていることが確認できます。
-
左上: サンプルパスが長期平均 $\mu = 1.0$ の周りを揺らいでいる。初期値 $x_0 = 3.0$ から出発したパスは、時間とともに長期平均に引き寄せられ、やがてその周りで定常的な振動を示します。黒い太線は理論的な期待値 $E[X_t] = \mu + (x_0 – \mu)e^{-\theta t}$ であり、パスの「中心線」として機能していることが分かります
-
右上: 分散が理論値に正確に一致し、定常値 $\sigma^2/(2\theta)$ に収束している。ブラウン運動の分散は $t$ に比例して際限なく増大しますが、OU過程の分散は有限値に飽和します。これはバネの復元力（平均回帰性）が揺らぎの拡散を抑えていることの帰結です
-
左下: $T = 5$ での分布が定常分布 $N(\mu, \sigma^2/(2\theta))$ と一致している。$\theta = 2$ の場合、定常分布への収束は $e^{-2\theta T} = e^{-20} \approx 2 \times 10^{-9}$ なので、$T = 5$ では完全に定常状態に到達しています
-
右下: 自己相関関数が理論値 $e^{-\theta\tau}$ と一致している。OU過程の自己相関は指数的に減衰し、相関時間は $1/\theta$ で特徴づけられます。$\theta = 2$ の場合、相関時間は $0.5$ であり、グラフからもこの時間スケールで自己相関がほぼゼロに減衰していることが確認できます
伊藤の公式の直接的な検証
最後に、$f(B_t) = B_t^2$ に対する伊藤の公式を数値的に検証します。伊藤の公式から $d(B_t^2) = 2B_t\,dB_t + dt$ なので、$B_T^2 = 2\int_0^T B_t\,dB_t + T$ が成り立つはずです。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(789)
T = 1.0
n = 10000
n_paths = 50000
dt = T / n
# ブラウン運動
dB = np.random.normal(0, np.sqrt(dt), (n_paths, n))
B = np.cumsum(dB, axis=1)
B = np.hstack([np.zeros((n_paths, 1)), B])
B_T = B[:, -1]
# 伊藤の公式の検証: B_T^2 = 2∫B_t dB_t + T
# 左辺: B_T^2
lhs = B_T**2
# 右辺: 2 * Σ B_{t_k} ΔB_k + T
ito_integral = np.sum(B[:, :-1] * dB, axis=1)
rhs = 2 * ito_integral + T
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
# 散布図
ax = axes[0]
idx = np.random.choice(n_paths, 3000, replace=False)
ax.scatter(lhs[idx], rhs[idx], alpha=0.1, s=3, color='steelblue')
lims = [min(lhs.min(), rhs.min()), max(lhs.max(), rhs.max())]
ax.plot(lims, lims, 'r-', linewidth=1.5, label='$y = x$')
ax.set_xlabel(r'$B_T^2$ (LHS)', fontsize=12)
ax.set_ylabel(r'$2\int_0^T B_t\,dB_t + T$ (RHS)', fontsize=12)
ax.set_title("Verification of Ito's formula for $B_t^2$", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
# 誤差のヒストグラム
ax = axes[1]
error = rhs - lhs
ax.hist(error, bins=80, density=True, alpha=0.7, color='coral',
edgecolor='white', linewidth=0.5)
ax.axvline(0, color='k', linestyle='--', linewidth=1.5)
ax.set_xlabel('RHS - LHS', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Numerical error distribution', fontsize=13)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ito_formula_verification.png', dpi=150, bbox_inches='tight')
plt.show()
rmse = np.sqrt(np.mean(error**2))
print(f"RMSE: {rmse:.8f}")
print(f"平均誤差: {np.mean(error):.8f}")
print(f"最大絶対誤差: {np.max(np.abs(error)):.8f}")
伊藤の公式の検証結果から、以下のことが分かります。
-
左図: $B_T^2$（左辺）と $2\int_0^T B_t\,dB_t + T$（右辺）がすべてのパスで正確に一致している。散布図の点が対角線上に密集しており、伊藤の公式 $B_T^2 = 2\int_0^T B_t\,dB_t + T$ がパスごとに成り立っていることが確認できます
-
右図: 数値誤差はゼロの周りに対称で非常に小さい。誤差はすべて離散化に起因するものであり、理論的にはゼロであるべきです。時間分割を細かくすればさらに小さくなります
まとめ
本記事では、伊藤の公式の導出と、それを用いたSDEの求解テクニックについて解説しました。
- 伊藤の公式は確率解析の連鎖律であり、通常の連鎖律に $\frac{1}{2}\sigma^2 f”$ という伊藤補正項が加わります。この補正項はブラウン運動の二次変分 $(dB_t)^2 = dt$ に起因します
- 幾何ブラウン運動 $dS = \mu S\,dt + \sigma S\,dB$ の解は $S_t = S_0\exp[(\mu – \sigma^2/2)t + \sigma B_t]$ です。$-\sigma^2/2$ の補正項（ボラティリティ・ドラッグ）は伊藤の公式の直接的な帰結です
- OU過程 $dX = -\theta(X – \mu)\,dt + \sigma\,dB$ は積分因子法で解け、平均回帰性と有限の定常分散 $\sigma^2/(2\theta)$ を持ちます
- SDEの求解テクニックとして、変数変換+伊藤の公式、積分因子法、確率指数関数の三つを紹介しました
- Pythonシミュレーションで、厳密解と数値解の一致、統計的性質の検証、伊藤の公式のパスごとの検証を行いました
伊藤の公式は、SDEの解を求めるだけでなく、確率過程の統計的性質の導出、確率微分方程式から偏微分方程式への変換（フォッカー・プランク方程式やコルモゴロフ方程式）など、確率解析のあらゆる場面で使われる普遍的な道具です。
次のステップとして、以下の記事も参考にしてください。
- フォッカー・プランク方程式 — SDEに対応する確率密度の時間発展方程式
- 伊藤積分の定義と性質 — 伊藤積分の厳密な構成
- マルチンゲール収束定理 — マルチンゲール理論の深い結果