気温は季節によって変動しますが、極端に高い日が続くとやがて平年並みに戻り、極端に低い日が続いてもやはり平年並みに戻ります。金利も同様で、高金利が永遠に続くことはなく、低金利が無限に続くこともありません。このような「ある水準に引き戻される」性質を数学的に表現したものがオルンシュタイン・ウーレンベック過程（Ornstein-Uhlenbeck process, OU過程）です。
ブラウン運動は時間とともに際限なく広がっていきますが、多くの物理量や経済指標は長期的に安定した水準を持ちます。「ブラウン運動に引き戻し力（復元力）を加えたらどうなるか？」— この問いから生まれたのがOU過程であり、平均回帰（mean-reverting）確率過程の最も基本的なモデルです。
オルンシュタイン・ウーレンベック過程を理解すると、以下のような応用が可能です。
- 金融工学: バシチェック・モデル（Vasicek model）による金利の期間構造のモデリング
- 統計物理学: ランジュバン方程式の解として、熱平衡下の粒子の速度分布を記述
- 生態学: 個体群サイズの環境収容力への回帰をモデル化
- 機械学習: 拡散モデル（DDPM/スコアベース生成モデル）のノイズスケジュールの基盤
- ペア取引: 金融工学における統計的裁定戦略で、共和分関係にある資産ペアのスプレッドをOU過程でモデル化
本記事の内容
- オルンシュタイン・ウーレンベック過程の定義とモチベーション
- 確率微分方程式の厳密解の導出
- 統計的性質（期待値・分散・共分散・定常分布）
- フォッカー・プランク方程式との関係
- パラメータ推定法
- 物理学・金融工学での応用
- Pythonによるシミュレーションと推定の実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ブラウン運動の性質 — OU過程の構成要素
- 正規分布 — OU過程の分布
- 幾何ブラウン運動と株価モデリング — SDEの解法（伊藤の公式）
物理的動機 — ランジュバン方程式から
ブラウン粒子の速度
1908年、ポール・ランジュバンはブラウン運動する粒子の速度 $V_t$ をモデル化しました。水中の花粉粒子は、粘性抵抗（速度に比例する減速力）とランダムな分子衝突（熱揺らぎ）を受けます。
ニュートンの運動方程式を書くと、
$$ m \frac{dV_t}{dt} = -\gamma V_t + \sigma \eta(t) $$
ここで $m$ は粒子の質量、$\gamma$ は粘性抵抗係数、$\sigma \eta(t)$ はランダム力（ホワイトノイズ）です。$\eta(t)$ は「形式的に $dB_t/dt$」と理解できますが、ブラウン運動は微分不可能なため、厳密な意味では $\eta(t)$ は関数として存在しません。
慣性を無視する過減衰極限（$m \to 0$ または $m/\gamma \to 0$）では、
$$ dV_t = -\frac{\gamma}{m} V_t \, dt + \frac{\sigma}{m} dB_t $$
$\theta = \gamma/m$（回帰速度）、$\sigma_{\text{OU}} = \sigma/m$（ノイズ強度）とおくと、
$$ dV_t = -\theta V_t \, dt + \sigma_{\text{OU}} \, dB_t $$
これが最も単純な形のオルンシュタイン・ウーレンベック過程のSDEです。
バネとランダム力のアナロジー
OU過程を直感的に理解するために、バネに取り付けられた粒子を想像してみましょう。
バネは粒子を平衡位置（原点）に引き戻そうとします。この復元力は、変位に比例し、$-\theta X_t$ で表されます。同時に、粒子はランダムな外力（たとえば熱揺らぎ）を受けて不規則に振動します。
- $\theta$ が大きい → バネが強い → 速く平衡に戻る → 変動が小さい
- $\theta$ が小さい → バネが弱い → ゆっくり平衡に戻る → 大きく変動する
- $\sigma$ が大きい → ランダム力が強い → 変動が大きい
- $\sigma$ が小さい → ランダム力が弱い → 変動が小さい
最終的には、復元力とランダム力のバランスで定常状態が決まります。このバランスは統計力学のゆらぎ散逸定理（fluctuation-dissipation theorem）の一形態です。
オルンシュタイン・ウーレンベック過程の定義
数学的定義
一般の長期平均 $\mu$ を持つオルンシュタイン・ウーレンベック過程は、以下のSDEで定義されます。
定義: 確率過程 $\{X_t\}_{t \geq 0}$ が次のSDEを満たすとき、パラメータ $(\theta, \mu, \sigma)$ のオルンシュタイン・ウーレンベック過程と呼ぶ。
$$ dX_t = \theta(\mu – X_t) \, dt + \sigma \, dB_t $$
各パラメータの意味は以下の通りです。
- $\theta > 0$: 回帰速度（speed of mean reversion）— $X_t$ が $\mu$ に引き戻される速さ。$1/\theta$ が特性時間スケール
- $\mu \in \mathbb{R}$: 長期平均（long-run mean）— $X_t$ が回帰する先の水準
- $\sigma > 0$: ボラティリティ（volatility）— ランダム変動の大きさ
- $B_t$: 標準ブラウン運動
SDEの右辺第1項 $\theta(\mu – X_t) \, dt$ はドリフト項です。$X_t > \mu$ のとき $\mu – X_t < 0$ なので下方向に引き戻され、$X_t < \mu$ のとき $\mu - X_t > 0$ なので上方向に引き戻されます。つまり、どちらの方向に離れても常に $\mu$ に向かう力が働きます。
ブラウン運動との関係
$\theta = 0$ とするとドリフトが消え、$dX_t = \sigma \, dB_t$ となります。これは単なるスケーリングされたブラウン運動です。$\theta > 0$ であることが平均回帰性の源であり、OU過程をブラウン運動から本質的に区別する特徴です。
$\theta \to \infty$ の極限では、$X_t \to \mu$（確定値）に収束します。バネが無限に強ければ、粒子は平衡位置に固定されます。
次に、このSDEを厳密に解いて、$X_t$ の明示的な表現を導きましょう。
厳密解の導出
積分因子法
OU過程のSDE $dX_t = \theta(\mu – X_t) \, dt + \sigma \, dB_t$ を解くために、常微分方程式の積分因子法を確率微分方程式に拡張します。
SDEを変形すると、
$$ dX_t + \theta X_t \, dt = \theta \mu \, dt + \sigma \, dB_t $$
積分因子 $e^{\theta t}$ を導入します。伊藤の公式より、$d(e^{\theta t} X_t)$ を計算しましょう。
$f(t, x) = e^{\theta t} x$ とおくと、$\frac{\partial f}{\partial t} = \theta e^{\theta t} x$, $\frac{\partial f}{\partial x} = e^{\theta t}$, $\frac{\partial^2 f}{\partial x^2} = 0$ です。
伊藤の公式を適用すると、
$$ d(e^{\theta t} X_t) = \theta e^{\theta t} X_t \, dt + e^{\theta t} dX_t + \frac{1}{2} \cdot 0 \cdot \sigma^2 \, dt $$
$f$ が $x$ について線形なので2階微分は0であり、伊藤補正項は消えます。$dX_t$ を代入すると、
$$ d(e^{\theta t} X_t) = \theta e^{\theta t} X_t \, dt + e^{\theta t}[\theta(\mu – X_t) \, dt + \sigma \, dB_t] $$
右辺を展開して整理します。
$$ = \theta e^{\theta t} X_t \, dt + \theta \mu e^{\theta t} \, dt – \theta e^{\theta t} X_t \, dt + \sigma e^{\theta t} \, dB_t $$
第1項と第3項が相殺し、
$$ d(e^{\theta t} X_t) = \theta \mu e^{\theta t} \, dt + \sigma e^{\theta t} \, dB_t $$
これは右辺が $X_t$ を含まない「単純な」確率微分方程式になりました。両辺を $[0, t]$ で積分します。
$$ e^{\theta t} X_t – X_0 = \theta \mu \int_0^t e^{\theta s} \, ds + \sigma \int_0^t e^{\theta s} \, dB_s $$
確定的積分を計算すると、
$$ \int_0^t e^{\theta s} \, ds = \frac{e^{\theta t} – 1}{\theta} $$
代入して $e^{\theta t}$ で割ると、
$$ X_t = X_0 e^{-\theta t} + \mu(1 – e^{-\theta t}) + \sigma \int_0^t e^{-\theta(t-s)} \, dB_s $$
これがOU過程の厳密解です。
解の構造の解釈
得られた解を3つの部分に分けて解釈しましょう。
$$ X_t = \underbrace{X_0 e^{-\theta t}}_{\text{初期値の減衰}} + \underbrace{\mu(1 – e^{-\theta t})}_{\text{長期平均への回帰}} + \underbrace{\sigma \int_0^t e^{-\theta(t-s)} \, dB_s}_{\text{ランダム変動}} $$
第1項: $X_0 e^{-\theta t}$ は初期条件 $X_0$ の影響が指数的に減衰する項です。$t \gg 1/\theta$ では無視できるほど小さくなります。$1/\theta$ は半減期のような時間スケールであり、初期条件の「記憶」が失われる速さを表します。
第2項: $\mu(1 – e^{-\theta t})$ は長期平均 $\mu$ への漸近を表す確定的な項です。$t \to \infty$ で $\mu$ に収束します。
第3項: $\sigma \int_0^t e^{-\theta(t-s)} \, dB_s$ はランダム変動です。被積分関数 $e^{-\theta(t-s)}$ は、時刻 $s$ でのノイズの寄与が時間 $t-s$ とともに指数的に減衰することを意味します。最近のノイズは大きく影響し、過去のノイズは忘れられます。
$\theta = 0$ の場合（ブラウン運動）と比較すると、OU過程では過去のノイズが $e^{-\theta(t-s)}$ の重みで減衰するため、変動が際限なく蓄積することはありません。これが平均回帰性の数学的メカニズムです。
厳密解が得られたので、次に統計的な性質を詳しく計算しましょう。
統計的性質
期待値
厳密解の各項の期待値を計算します。伊藤積分 $\int_0^t e^{-\theta(t-s)} \, dB_s$ の期待値は0です（伊藤積分の基本性質）。したがって、
$$ E[X_t] = X_0 e^{-\theta t} + \mu(1 – e^{-\theta t}) $$
$t \to \infty$ で $E[X_t] \to \mu$ です。初期値 $X_0$ がどんな値であっても、期待値は長期平均 $\mu$ に指数的に収束します。収束の時定数は $1/\theta$ です。
分散
分散を計算するために、伊藤等長式（Itô isometry）を使います。確定的関数 $f(s)$ に対して、
$$ E\left[\left(\int_0^t f(s) \, dB_s\right)^2\right] = \int_0^t f(s)^2 \, ds $$
$f(s) = \sigma e^{-\theta(t-s)}$ を代入すると、
$$ \text{Var}(X_t) = E\left[\left(\sigma \int_0^t e^{-\theta(t-s)} \, dB_s\right)^2\right] = \sigma^2 \int_0^t e^{-2\theta(t-s)} \, ds $$
$u = t – s$ と置換すると、
$$ \text{Var}(X_t) = \sigma^2 \int_0^t e^{-2\theta u} \, du = \sigma^2 \left[\frac{1 – e^{-2\theta t}}{2\theta}\right] = \frac{\sigma^2}{2\theta}(1 – e^{-2\theta t}) $$
$t \to \infty$ で分散は定常値に収束します。
$$ \text{Var}(X_t) \to \frac{\sigma^2}{2\theta} \quad (t \to \infty) $$
定常分散 $\sigma^2/(2\theta)$ は、ノイズ強度 $\sigma^2$ と回帰速度 $\theta$ のバランスで決まります。$\theta$ が大きい（強い復元力）ほど分散は小さく、$\sigma$ が大きい（強いノイズ）ほど分散は大きくなります。
一方、ブラウン運動（$\theta = 0$）では $\text{Var}(B_t) = t \to \infty$ と分散が際限なく増大することと対照的です。
共分散
$s \leq t$ のとき、$X_s$ と $X_t$ の共分散を計算しましょう。厳密解から、
$$ X_s = X_0 e^{-\theta s} + \mu(1 – e^{-\theta s}) + \sigma \int_0^s e^{-\theta(s-u)} \, dB_u $$
$$ X_t = X_0 e^{-\theta t} + \mu(1 – e^{-\theta t}) + \sigma \int_0^t e^{-\theta(t-u)} \, dB_u $$
確率的な部分のみが共分散に寄与します。$\int_0^s$ と $\int_0^t$ の共通部分は $\int_0^s$ なので、
$$ \text{Cov}(X_s, X_t) = \sigma^2 \int_0^s e^{-\theta(s-u)} e^{-\theta(t-u)} \, du = \sigma^2 e^{-\theta(t+s)} \int_0^s e^{2\theta u} \, du $$
積分を実行すると、
$$ \text{Cov}(X_s, X_t) = \frac{\sigma^2}{2\theta} e^{-\theta(t+s)}(e^{2\theta s} – 1) = \frac{\sigma^2}{2\theta}(e^{-\theta(t-s)} – e^{-\theta(t+s)}) $$
$s, t$ がともに大きい（定常状態）場合、$e^{-\theta(t+s)} \approx 0$ なので、
$$ \text{Cov}(X_s, X_t) \approx \frac{\sigma^2}{2\theta} e^{-\theta(t-s)} = \frac{\sigma^2}{2\theta} e^{-\theta|t-s|} $$
この結果は重要です。定常状態での共分散は時間差 $|t-s|$ のみに依存し、指数的に減衰します。時間差が $1/\theta$ 程度になると相関がほぼゼロになります。$1/\theta$ は相関時間（correlation time）とも呼ばれます。
正規性と遷移確率
$X_t$ は正規確率変数の線形結合（初期値 + 伊藤積分）なので、任意の $t$ で正規分布に従います。条件付き分布は、
$$ X_t \mid X_s \sim N\left(X_s e^{-\theta(t-s)} + \mu(1 – e^{-\theta(t-s)}), \; \frac{\sigma^2}{2\theta}(1 – e^{-2\theta(t-s)})\right) $$
遷移確率密度（transition density）は、
$$ p(x, t \mid y, s) = \frac{1}{\sqrt{2\pi v(t-s)}} \exp\left(-\frac{(x – m(y, t-s))^2}{2v(t-s)}\right) $$
ここで $m(y, \tau) = ye^{-\theta\tau} + \mu(1 – e^{-\theta\tau})$ は条件付き期待値、$v(\tau) = \frac{\sigma^2}{2\theta}(1 – e^{-2\theta\tau})$ は条件付き分散です。
この遷移確率密度が明示的に求められることは、OU過程の大きな利点です。多くのSDEでは遷移確率密度が閉じた形で書けませんが、OU過程では正規分布という単純な形になります。
定常分布
$t \to \infty$ の極限で、$X_t$ の分布は初期条件によらず以下の定常分布に収束します。
$$ X_{\infty} \sim N\left(\mu, \frac{\sigma^2}{2\theta}\right) $$
つまり、十分に長い時間が経つと、OU過程は平均 $\mu$、分散 $\sigma^2/(2\theta)$ の正規分布に従って揺らぎ続けます。これは定常正規過程（stationary Gaussian process）です。
定常分布の存在は、OU過程がブラウン運動と根本的に異なることを示しています。ブラウン運動には定常分布が存在しません（分散が無限に増大するため）。
この結果を物理の言葉で言い換えると、系は熱平衡（thermal equilibrium）に達し、ゆらぎ散逸定理 $k_B T = \sigma^2/(2\gamma/m)$ のような関係が成り立つことに対応します。
OU過程の統計的性質が完全に特徴付けられたところで、確率密度関数の時間発展を記述するフォッカー・プランク方程式との関係を見ていきましょう。
フォッカー・プランク方程式
OU過程に対するフォッカー・プランク方程式
一般のSDE $dX_t = a(X_t) \, dt + b(X_t) \, dB_t$ に対して、確率密度 $p(x, t)$ の時間発展はフォッカー・プランク方程式（Fokker-Planck equation, 前進コルモゴロフ方程式とも呼ばれる）で記述されます。
$$ \frac{\partial p}{\partial t} = -\frac{\partial}{\partial x}[a(x) p] + \frac{1}{2} \frac{\partial^2}{\partial x^2}[b(x)^2 p] $$
OU過程では $a(x) = \theta(\mu – x)$, $b(x) = \sigma$ なので、
$$ \frac{\partial p}{\partial t} = -\frac{\partial}{\partial x}[\theta(\mu – x) p] + \frac{\sigma^2}{2} \frac{\partial^2 p}{\partial x^2} $$
右辺第1項を展開すると、
$$ \frac{\partial p}{\partial t} = \theta p + \theta(x – \mu) \frac{\partial p}{\partial x} + \frac{\sigma^2}{2} \frac{\partial^2 p}{\partial x^2} $$
定常解の導出
定常分布を求めるために $\frac{\partial p}{\partial t} = 0$ とおきます。$\mu = 0$（簡単のため）とすると、
$$ 0 = -\frac{\partial}{\partial x}[-\theta x \cdot p_{\text{st}}] + \frac{\sigma^2}{2} \frac{\partial^2 p_{\text{st}}}{\partial x^2} $$
1回積分すると（境界条件 $p_{\text{st}} \to 0$ as $x \to \pm\infty$ より積分定数は0）、
$$ 0 = \theta x \cdot p_{\text{st}} + \frac{\sigma^2}{2} \frac{dp_{\text{st}}}{dx} $$
これは $p_{\text{st}}$ に関する1階常微分方程式です。変数分離すると、
$$ \frac{dp_{\text{st}}}{p_{\text{st}}} = -\frac{2\theta x}{\sigma^2} dx $$
両辺を積分すると、
$$ \ln p_{\text{st}} = -\frac{\theta x^2}{\sigma^2} + C $$
指数を取ると、
$$ p_{\text{st}}(x) = A \exp\left(-\frac{\theta x^2}{\sigma^2}\right) $$
正規化定数を求めると $A = \sqrt{\theta/(\pi \sigma^2)}$ であり、
$$ p_{\text{st}}(x) = \sqrt{\frac{\theta}{\pi \sigma^2}} \exp\left(-\frac{\theta x^2}{\sigma^2}\right) = \frac{1}{\sqrt{2\pi \cdot \sigma^2/(2\theta)}} \exp\left(-\frac{x^2}{2 \cdot \sigma^2/(2\theta)}\right) $$
これは $N(0, \sigma^2/(2\theta))$ の密度関数であり、先に求めた定常分布と一致します。一般の $\mu \neq 0$ の場合は $N(\mu, \sigma^2/(2\theta))$ となります。
フォッカー・プランク方程式からの導出は、確率密度関数のレベルで定常分布を求める別のアプローチであり、SDEからの直接計算と整合しています。
次に、観測データからOU過程のパラメータを推定する方法を見ていきましょう。
パラメータ推定
離散観測データからの推定
実データは離散的な時刻 $t_0 < t_1 < \cdots < t_n$（等間隔 $\Delta t$）で観測されます。遷移確率密度が正規分布であることを利用して、最尤推定（maximum likelihood estimation）を行えます。
等間隔サンプリング $\Delta t$ で観測されたデータ $x_0, x_1, \dots, x_n$ に対して、遷移密度から対数尤度関数を構成します。
$X_{t+\Delta t} \mid X_t = x$ の条件付き分布は、
$$ N\left(xe^{-\theta\Delta t} + \mu(1 – e^{-\theta\Delta t}), \; \frac{\sigma^2}{2\theta}(1 – e^{-2\theta\Delta t})\right) $$
$\phi = e^{-\theta\Delta t}$ とおくと、$X_{i+1} \mid X_i = x_i$ は $N(\phi x_i + \mu(1-\phi), v)$ に従います。ここで $v = \frac{\sigma^2}{2\theta}(1 – \phi^2)$ です。
対数尤度は、
$$ \ell(\phi, \mu, v) = -\frac{n}{2} \ln(2\pi v) – \frac{1}{2v} \sum_{i=0}^{n-1} (x_{i+1} – \phi x_i – \mu(1-\phi))^2 $$
これは実質的に $x_{i+1}$ を $x_i$ に回帰するAR(1)モデル（1次自己回帰モデル）と同じ構造です。
最尤推定量
$\alpha = \mu(1 – \phi)$, $\beta = \phi$ とおくと、$x_{i+1} = \alpha + \beta x_i + \varepsilon_i$（$\varepsilon_i \sim N(0, v)$）という線形回帰になります。
最小二乗推定量は、
$$ \hat{\beta} = \frac{\sum_{i=0}^{n-1} (x_i – \bar{x})(x_{i+1} – \bar{x}_+)}{\sum_{i=0}^{n-1} (x_i – \bar{x})^2} $$
$$ \hat{\alpha} = \bar{x}_+ – \hat{\beta} \bar{x} $$
ここで $\bar{x} = \frac{1}{n}\sum_{i=0}^{n-1} x_i$, $\bar{x}_+ = \frac{1}{n}\sum_{i=0}^{n-1} x_{i+1}$ です。
元のパラメータに戻すと、
$$ \hat{\theta} = -\frac{\ln \hat{\beta}}{\Delta t}, \quad \hat{\mu} = \frac{\hat{\alpha}}{1 – \hat{\beta}}, \quad \hat{\sigma}^2 = \frac{2\hat{\theta}\hat{v}}{1 – \hat{\beta}^2} $$
ここで $\hat{v} = \frac{1}{n}\sum_{i=0}^{n-1}(x_{i+1} – \hat{\alpha} – \hat{\beta}x_i)^2$ は残差分散の推定量です。
推定の注意点
パラメータ推定にはいくつかの注意点があります。
1. $\hat{\beta}$ のバイアス: 有限サンプルでは $\hat{\beta}$ は下方バイアスを持つことが知られています（AR(1)モデルの推定量のバイアス）。データ数が少ない場合にはバイアス補正が必要です。
2. 観測頻度と推定精度: $\Delta t$ が大きすぎると指数減衰が進んで $\hat{\beta} \approx 0$ となり推定が不安定になります。逆に $\Delta t$ が小さすぎると $\hat{\beta} \approx 1$ となり、$\hat{\theta} \approx 0$ 付近での推定精度が悪くなります。
3. 平均回帰の検定: $\hat{\beta} < 1$ であるかどうかを統計的に検定する必要があります。$\beta = 1$ の帰無仮説の下では、$X_t$ はランダムウォーク（単位根過程）であり、平均回帰性がありません。これはディッキー・フラー検定（Dickey-Fuller test）として知られています。
理論的な推定法を理解したところで、Pythonでシミュレーションとパラメータ推定を実装してみましょう。
Pythonによるシミュレーションと推定
OU過程のシミュレーション
OU過程のサンプルパスを厳密な遷移確率を用いて生成します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# パラメータ
theta = 2.0 # 回帰速度
mu = 5.0 # 長期平均
sigma = 1.0 # ボラティリティ
X0 = 0.0 # 初期値
T = 5.0 # 期間
N = 1000 # ステップ数
dt = T / N
n_paths = 8
t = np.linspace(0, T, N + 1)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 異なる初期値からのパス
ax = axes[0]
initial_values = [-3, 0, 3, 5, 8, 10, 13]
colors = plt.cm.viridis(np.linspace(0, 1, len(initial_values)))
for X0_val, color in zip(initial_values, colors):
X = np.zeros(N + 1)
X[0] = X0_val
phi = np.exp(-theta * dt)
mean_incr = mu * (1 - phi)
std_incr = np.sqrt(sigma**2 / (2 * theta) * (1 - phi**2))
for i in range(N):
X[i+1] = phi * X[i] + mean_incr + std_incr * np.random.normal()
ax.plot(t, X, color=color, linewidth=0.7, alpha=0.8, label=f'X₀={X0_val}')
ax.axhline(y=mu, color='red', linestyle='--', linewidth=2, label=f'μ={mu}')
# 定常分布の±2σ範囲
stat_std = sigma / np.sqrt(2 * theta)
ax.axhspan(mu - 2*stat_std, mu + 2*stat_std, alpha=0.1, color='red')
ax.set_xlabel("Time")
ax.set_ylabel("X(t)")
ax.set_title("OU Process: Different Initial Values")
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)
# 右: パラメータの影響
ax = axes[1]
configs = [
(0.5, 5.0, 1.0, 'θ=0.5 (slow reversion)'),
(2.0, 5.0, 1.0, 'θ=2.0 (medium)'),
(10.0, 5.0, 1.0, 'θ=10.0 (fast reversion)'),
]
X0_val = 0.0
for theta_val, mu_val, sigma_val, label in configs:
X = np.zeros(N + 1)
X[0] = X0_val
phi = np.exp(-theta_val * dt)
mean_incr = mu_val * (1 - phi)
std_incr = np.sqrt(sigma_val**2 / (2 * theta_val) * (1 - phi**2))
for i in range(N):
X[i+1] = phi * X[i] + mean_incr + std_incr * np.random.normal()
ax.plot(t, X, linewidth=0.8, label=label)
ax.axhline(y=5.0, color='red', linestyle='--', linewidth=1.5, label='μ=5')
ax.set_xlabel("Time")
ax.set_ylabel("X(t)")
ax.set_title("Effect of Mean-Reversion Speed θ")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、異なる初期値（$X_0 = -3$ から $X_0 = 13$）から出発したパスがすべて長期平均 $\mu = 5$（赤い破線）に向かって指数的に収束している様子が確認できます。薄い赤の帯は定常分布の $\pm 2\sigma$ 範囲を表しており、十分な時間が経つとパスはこの帯の中で揺らぎます。初期値の「記憶」は時定数 $1/\theta = 0.5$ 程度で失われています。
右のパネルでは、回帰速度 $\theta$ の違いが経路の振る舞いに与える影響が可視化されています。$\theta = 0.5$ のパスはゆっくりと回帰し大きく変動しますが、$\theta = 10$ のパスは急速に $\mu$ に到達し、小さな振幅で揺らぎます。定常分散 $\sigma^2/(2\theta)$ が $\theta$ に反比例するため、$\theta$ が大きいほど変動が小さくなることが視覚的に確認できます。
期待値と分散の理論値との比較
多数のパスをシミュレーションし、サンプル統計量が理論値と一致することを確認します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
theta, mu_param, sigma = 2.0, 5.0, 1.0
X0 = 0.0
T, N = 5.0, 1000
dt = T / N
n_paths = 10000
t = np.linspace(0, T, N + 1)
# 多数パスのシミュレーション
paths = np.zeros((n_paths, N + 1))
paths[:, 0] = X0
phi = np.exp(-theta * dt)
mean_incr = mu_param * (1 - phi)
std_incr = np.sqrt(sigma**2 / (2 * theta) * (1 - phi**2))
for i in range(N):
paths[:, i+1] = phi * paths[:, i] + mean_incr + std_incr * np.random.normal(size=n_paths)
# 理論値
E_theory = X0 * np.exp(-theta * t) + mu_param * (1 - np.exp(-theta * t))
V_theory = sigma**2 / (2 * theta) * (1 - np.exp(-2 * theta * t))
# サンプル統計量
E_sample = np.mean(paths, axis=0)
V_sample = np.var(paths, axis=0)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
ax = axes[0]
ax.plot(t, E_sample, 'b-', linewidth=1.5, label='Sample mean')
ax.plot(t, E_theory, 'r--', linewidth=2, label='Theory E[X(t)]')
ax.set_xlabel("Time")
ax.set_ylabel("E[X(t)]")
ax.set_title("Mean: Theory vs Simulation")
ax.legend()
ax.grid(True, alpha=0.3)
ax = axes[1]
ax.plot(t, V_sample, 'b-', linewidth=1.5, label='Sample variance')
ax.plot(t, V_theory, 'r--', linewidth=2, label='Theory Var(X(t))')
ax.axhline(y=sigma**2/(2*theta), color='green', linestyle=':',
linewidth=1.5, label=f'σ²/(2θ) = {sigma**2/(2*theta):.3f}')
ax.set_xlabel("Time")
ax.set_ylabel("Var(X(t))")
ax.set_title("Variance: Theory vs Simulation")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、サンプル平均（青）が理論的な期待値（赤破線）と完全に一致しています。$X_0 = 0$ から出発し、$\mu = 5$ に向かって指数的に収束する様子が確認できます。右のパネルでは、サンプル分散（青）が理論的な分散（赤破線）と一致しており、$t \to \infty$ で定常値 $\sigma^2/(2\theta) = 0.25$（緑の点線）に収束することも確認できます。10,000パスの平均によって統計的ゆらぎが十分に小さくなり、理論と数値の一致が明確に見えています。
定常分布と自己相関関数の検証
定常状態でのヒストグラムと自己相関関数が理論と一致するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
theta, mu_param, sigma = 2.0, 5.0, 1.0
T, N = 100.0, 20000
dt = T / N
t = np.linspace(0, T, N + 1)
# 長い1本のパスをシミュレーション
X = np.zeros(N + 1)
X[0] = mu_param # 定常状態から出発
phi = np.exp(-theta * dt)
mean_incr = mu_param * (1 - phi)
std_incr = np.sqrt(sigma**2 / (2 * theta) * (1 - phi**2))
for i in range(N):
X[i+1] = phi * X[i] + mean_incr + std_incr * np.random.normal()
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 定常分布のヒストグラム
ax = axes[0]
# 最初の過渡部分を除外
X_stationary = X[N//10:]
ax.hist(X_stationary, bins=80, density=True, alpha=0.7, color='steelblue',
label='Simulation')
x_range = np.linspace(mu_param - 4*sigma/np.sqrt(2*theta),
mu_param + 4*sigma/np.sqrt(2*theta), 200)
pdf_theory = stats.norm.pdf(x_range, loc=mu_param,
scale=sigma/np.sqrt(2*theta))
ax.plot(x_range, pdf_theory, 'r-', linewidth=2,
label=f'N({mu_param}, {sigma**2/(2*theta):.3f})')
ax.set_xlabel("X")
ax.set_ylabel("Density")
ax.set_title("Stationary Distribution")
ax.legend()
ax.grid(True, alpha=0.3)
# 右: 自己相関関数
ax = axes[1]
max_lag = 500
X_centered = X_stationary - np.mean(X_stationary)
var_X = np.var(X_stationary)
acf = np.correlate(X_centered[:5000], X_centered[:5000], mode='full')
acf = acf[len(acf)//2:]
acf = acf / acf[0]
lag_time = np.arange(max_lag) * dt
ax.plot(lag_time, acf[:max_lag], 'b-', linewidth=1.5, label='Sample ACF')
ax.plot(lag_time, np.exp(-theta * lag_time), 'r--', linewidth=2,
label=f'exp(-{theta}τ)')
ax.set_xlabel("Lag τ")
ax.set_ylabel("Autocorrelation")
ax.set_title("Autocorrelation Function")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、OU過程の定常分布が理論的な正規分布 $N(\mu, \sigma^2/(2\theta))$ と完全に一致しています。ヒストグラムの形状は左右対称で、裾の部分まで正規分布のフィットがよいことが確認できます。
右のパネルでは、自己相関関数が理論的な指数減衰 $e^{-\theta\tau}$ と一致しています。$\theta = 2$ なので、ラグ $\tau = 0.5$（$= 1/\theta$）程度で相関が $e^{-1} \approx 0.37$ に減衰し、$\tau = 1.5$ 程度でほぼゼロになります。これは「過去のノイズが指数的に忘れられる」という性質の数値的な確認です。
パラメータ推定の実装と検証
離散観測データからOU過程のパラメータを推定し、真の値との比較を行います。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# 真のパラメータ
theta_true = 2.0
mu_true = 5.0
sigma_true = 1.0
# データ生成
T, dt = 50.0, 0.01
N = int(T / dt)
X = np.zeros(N + 1)
X[0] = mu_true
phi_true = np.exp(-theta_true * dt)
mean_incr = mu_true * (1 - phi_true)
std_incr = np.sqrt(sigma_true**2 / (2 * theta_true) * (1 - phi_true**2))
for i in range(N):
X[i+1] = phi_true * X[i] + mean_incr + std_incr * np.random.normal()
# 最尤推定（AR(1)回帰）
x_curr = X[:-1]
x_next = X[1:]
n = len(x_curr)
# 回帰係数の推定
beta_hat = np.sum((x_curr - np.mean(x_curr)) * (x_next - np.mean(x_next))) / \
np.sum((x_curr - np.mean(x_curr))**2)
alpha_hat = np.mean(x_next) - beta_hat * np.mean(x_curr)
# OU過程パラメータに変換
theta_hat = -np.log(beta_hat) / dt
mu_hat = alpha_hat / (1 - beta_hat)
residuals = x_next - alpha_hat - beta_hat * x_curr
v_hat = np.mean(residuals**2)
sigma_hat = np.sqrt(2 * theta_hat * v_hat / (1 - beta_hat**2))
print("パラメータ推定結果:")
print(f" θ: 真値 = {theta_true:.4f}, 推定値 = {theta_hat:.4f}")
print(f" μ: 真値 = {mu_true:.4f}, 推定値 = {mu_hat:.4f}")
print(f" σ: 真値 = {sigma_true:.4f}, 推定値 = {sigma_hat:.4f}")
# 推定のサンプルサイズ依存性
n_experiments = 200
sample_sizes = [100, 500, 1000, 5000, 10000]
results = {n_s: {'theta': [], 'mu': [], 'sigma': []} for n_s in sample_sizes}
for n_s in sample_sizes:
for _ in range(n_experiments):
# データ生成
X_exp = np.zeros(n_s + 1)
X_exp[0] = mu_true
for i in range(n_s):
X_exp[i+1] = phi_true * X_exp[i] + mean_incr + std_incr * np.random.normal()
# 推定
xc = X_exp[:-1]
xn = X_exp[1:]
b = np.sum((xc - np.mean(xc)) * (xn - np.mean(xn))) / np.sum((xc - np.mean(xc))**2)
a = np.mean(xn) - b * np.mean(xc)
if 0 < b < 1: # 平均回帰の場合のみ
th = -np.log(b) / dt
m = a / (1 - b)
res = xn - a - b * xc
v = np.mean(res**2)
s = np.sqrt(2 * th * v / (1 - b**2))
results[n_s]['theta'].append(th)
results[n_s]['mu'].append(m)
results[n_s]['sigma'].append(s)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
param_names = ['theta', 'mu', 'sigma']
true_values = [theta_true, mu_true, sigma_true]
labels = ['θ', 'μ', 'σ']
for ax, param, true_val, label in zip(axes, param_names, true_values, labels):
means = [np.mean(results[n_s][param]) for n_s in sample_sizes]
stds = [np.std(results[n_s][param]) for n_s in sample_sizes]
ax.errorbar(sample_sizes, means, yerr=stds, fmt='bo-', capsize=4, markersize=5)
ax.axhline(y=true_val, color='red', linestyle='--', linewidth=1.5,
label=f'True {label} = {true_val}')
ax.set_xlabel("Sample size")
ax.set_ylabel(f"Estimated {label}")
ax.set_title(f"Estimation of {label}")
ax.set_xscale('log')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
3つのパネルは、サンプルサイズに対する各パラメータの推定値の収束を示しています。エラーバーは200回のモンテカルロ実験の標準偏差です。すべてのパラメータで、サンプルサイズが増えるにつれて推定値が真の値（赤い破線）に収束し、標準偏差が縮小していくことが確認できます。$\theta$ の推定は $\mu$ や $\sigma$ よりもバラツキが大きい傾向がありますが、これはAR(1)係数 $\hat{\beta}$ の対数変換を経るため、$\beta \approx 1$ 付近で感度が高くなることに起因しています。
応用: バシチェック金利モデル
モデルの定式化
バシチェック・モデル（Vasicek model, 1977年）は、短期金利 $r_t$ をOU過程でモデル化する最も古典的な金利モデルです。
$$ dr_t = \kappa(\bar{r} – r_t) \, dt + \sigma_r \, dB_t $$
ここで $\kappa > 0$ は回帰速度、$\bar{r}$ は長期平均金利、$\sigma_r$ はボラティリティです。
バシチェック・モデルの利点は、ゼロクーポン債の価格が閉じた形で求められることです。満期 $T$ のゼロクーポン債の価格は、
$$ P(t, T) = A(t, T) \exp(-B(t, T) r_t) $$
ここで、
$$ B(t, T) = \frac{1 – e^{-\kappa(T-t)}}{\kappa} $$
$$ A(t, T) = \exp\left[\left(\bar{r} – \frac{\sigma_r^2}{2\kappa^2}\right)(B(t,T) – (T-t)) – \frac{\sigma_r^2}{4\kappa} B(t,T)^2\right] $$
モデルの限界
バシチェック・モデルはOU過程をそのまま使うため、金利 $r_t$ が正規分布に従います。このため金利が負になる確率が正です。2010年代以降のマイナス金利環境では結果的に適切な性質ですが、歴史的には欠点と見なされていました。
この問題を解決するために、CIRモデル（Cox-Ingersoll-Ross model）が提案されています。CIRモデルは、
$$ dr_t = \kappa(\bar{r} – r_t) \, dt + \sigma_r \sqrt{r_t} \, dB_t $$
ボラティリティが $\sigma_r \sqrt{r_t}$ と $r_t$ の平方根に比例するため、$r_t \geq 0$ が保証されます（フェラー条件 $2\kappa\bar{r} \geq \sigma_r^2$ の下で）。
応用: ペア取引（統計的裁定）
OU過程によるスプレッドのモデリング
ペア取引は、相関の高い2つの資産のスプレッド（価格差や比率）が平均に回帰する性質を利用する取引戦略です。
2つの株価 $S_t^{(1)}$, $S_t^{(2)}$ が長期的な均衡関係（共和分関係）にある場合、スプレッド $Z_t = S_t^{(1)} – \beta S_t^{(2)}$ がOU過程に従うとモデル化できます。
$$ dZ_t = \theta(\mu_Z – Z_t) \, dt + \sigma_Z \, dB_t $$
ペア取引の基本戦略は以下の通りです。
- $Z_t > \mu_Z + k\sigma_Z/\sqrt{2\theta}$ のとき: スプレッドが「広がりすぎている」ので、$S^{(1)}$ を売り、$S^{(2)}$ を買う
- $Z_t < \mu_Z - k\sigma_Z/\sqrt{2\theta}$ のとき: スプレッドが「縮まりすぎている」ので、$S^{(1)}$ を買い、$S^{(2)}$ を売る
- $Z_t \approx \mu_Z$ のとき: ポジションを清算する
OU過程の平均回帰性により、スプレッドは長期平均に戻ることが期待されるため、理論的には利益が得られます。
ペア取引のシミュレーション
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# OU過程によるスプレッドの生成
theta, mu_spread, sigma_spread = 5.0, 0.0, 0.5
T, dt = 2.0, 0.001
N = int(T / dt)
t = np.linspace(0, T, N + 1)
Z = np.zeros(N + 1)
Z[0] = 0.0
phi = np.exp(-theta * dt)
m_inc = mu_spread * (1 - phi)
s_inc = np.sqrt(sigma_spread**2 / (2 * theta) * (1 - phi**2))
for i in range(N):
Z[i+1] = phi * Z[i] + m_inc + s_inc * np.random.normal()
# 取引シグナル
stat_std = sigma_spread / np.sqrt(2 * theta)
upper = mu_spread + 2 * stat_std
lower = mu_spread - 2 * stat_std
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
# 上: スプレッドとシグナル
ax = axes[0]
ax.plot(t, Z, 'steelblue', linewidth=0.8)
ax.axhline(y=mu_spread, color='black', linestyle='-', linewidth=1)
ax.axhline(y=upper, color='red', linestyle='--', linewidth=1, label=f'Upper (+2σ = {upper:.3f})')
ax.axhline(y=lower, color='green', linestyle='--', linewidth=1, label=f'Lower (-2σ = {lower:.3f})')
ax.fill_between(t, lower, upper, alpha=0.1, color='gray')
ax.set_xlabel("Time (years)")
ax.set_ylabel("Spread Z(t)")
ax.set_title("Pairs Trading: OU Process Spread")
ax.legend()
ax.grid(True, alpha=0.3)
# 下: 累積損益
ax = axes[1]
position = 0 # -1: short, 0: flat, 1: long
entry_price = 0.0
pnl = np.zeros(N + 1)
cumulative_pnl = 0.0
for i in range(1, N + 1):
if position == 0:
if Z[i] > upper:
position = -1
entry_price = Z[i]
elif Z[i] < lower:
position = 1
entry_price = Z[i]
elif position == 1:
if Z[i] > mu_spread:
cumulative_pnl += Z[i] - entry_price
position = 0
elif position == -1:
if Z[i] < mu_spread:
cumulative_pnl += entry_price - Z[i]
position = 0
pnl[i] = cumulative_pnl
ax.plot(t, pnl, 'green', linewidth=1.5)
ax.set_xlabel("Time (years)")
ax.set_ylabel("Cumulative P&L")
ax.set_title("Cumulative Profit/Loss from Pairs Trading")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上のパネルでは、OU過程に従うスプレッド $Z_t$ が長期平均 $\mu = 0$ の周りを平均回帰的に変動しています。赤と緑の破線は $\pm 2$ 定常標準偏差のレベルで、これをエントリーシグナルとして使います。
下のパネルは、シンプルなペア取引戦略の累積損益を示しています。スプレッドが上方閾値を超えたらショート（売り）、下方閾値を下回ったらロング（買い）し、平均に戻ったら清算するという戦略です。OU過程の平均回帰性により、累積損益は全体的に右上がりのトレンドを示しています。ただし、実際の金融市場ではスプレッドが必ずしもOU過程に従う保証はなく、取引コストやスリッページ、共和分関係の崩壊などのリスクがある点に留意が必要です。
多次元オルンシュタイン・ウーレンベック過程
定義と厳密解
多次元OU過程は、ベクトル値の確率過程 $\bm{X}_t \in \mathbb{R}^d$ に対して定義されます。
$$ d\bm{X}_t = \bm{\Theta}(\bm{\mu} – \bm{X}_t) \, dt + \bm{\Sigma} \, d\bm{B}_t $$
ここで $\bm{\Theta}$ は $d \times d$ の行列（回帰速度行列）、$\bm{\mu}$ は $d$ 次元の長期平均ベクトル、$\bm{\Sigma}$ は $d \times d$ の行列（ボラティリティ行列）、$\bm{B}_t$ は $d$ 次元標準ブラウン運動です。
厳密解は、
$$ \bm{X}_t = e^{-\bm{\Theta} t} \bm{X}_0 + (\bm{I} – e^{-\bm{\Theta} t}) \bm{\mu} + \int_0^t e^{-\bm{\Theta}(t-s)} \bm{\Sigma} \, d\bm{B}_s $$
ここで $e^{-\bm{\Theta} t}$ は行列指数関数です。
定常分布は $N(\bm{\mu}, \bm{V}_{\infty})$ であり、$\bm{V}_{\infty}$ は連続リアプノフ方程式
$$ \bm{\Theta} \bm{V}_{\infty} + \bm{V}_{\infty} \bm{\Theta}^{\top} = \bm{\Sigma} \bm{\Sigma}^{\top} $$
の解です。1次元の場合 $\theta \cdot 2V_{\infty} = \sigma^2$、すなわち $V_{\infty} = \sigma^2/(2\theta)$ に帰着し、先の結果と一致します。
拡散モデルとの関係
近年の機械学習における拡散モデル（diffusion models）は、OU過程と深い関係があります。スコアベース生成モデルでは、データ分布に段階的にノイズを加えてガウス分布に変換する「前方過程」（forward process）として、OU過程に類似したSDEを使います。
$$ d\bm{X}_t = -\frac{1}{2} \beta(t) \bm{X}_t \, dt + \sqrt{\beta(t)} \, d\bm{B}_t $$
$\beta(t)$ は時間依存のノイズスケジュールであり、$\beta(t) = \beta$（定数）の場合はOU過程そのもの（$\theta = \beta/2$, $\mu = 0$, $\sigma = \sqrt{\beta}$）です。
生成過程は、この前方SDEの時間反転SDEを解くことで実現されます。OU過程の解析的な性質（遷移確率が正規分布であることなど）が、効率的な学習と生成を可能にしています。
まとめ
本記事では、オルンシュタイン・ウーレンベック過程の理論を体系的に解説しました。
- 物理的動機: ランジュバン方程式から自然に導かれる平均回帰確率過程
- 厳密解: 積分因子法により $X_t = X_0 e^{-\theta t} + \mu(1 – e^{-\theta t}) + \sigma\int_0^t e^{-\theta(t-s)} dB_s$ を導出
- 統計的性質: 期待値は指数的に $\mu$ に収束、分散は $\sigma^2/(2\theta)$ に収束、共分散は $e^{-\theta|t-s|}$ で指数減衰
- 定常分布: $N(\mu, \sigma^2/(2\theta))$ — ゆらぎ散逸定理の数学的表現
- パラメータ推定: AR(1)回帰による最尤推定
- 応用: バシチェック金利モデル、ペア取引、拡散モデル
OU過程は「ブラウン運動 + 平均回帰」という単純な構造ながら、解析的に豊かな性質を持ち、理論と応用の両面で中心的な役割を果たしています。多くのSDEの解析の出発点となるモデルであり、確率解析を学ぶ上で避けて通れない基本的なプロセスです。
次のステップとして、以下の記事も参考にしてください。
- 伊藤積分の定義と性質 — OU過程の厳密解で使用した伊藤積分
- フォッカー・プランク方程式と確率密度の時間発展 — 確率密度関数の発展方程式
- 幾何ブラウン運動と株価モデリング — 別のSDE解法の例