株価のチャートを眺めると、価格は不規則に上下しながらも長期的にはトレンドを持つように見えます。「株価の変動をどのように数学的にモデル化すればよいか？」— この問いに答えるのが幾何ブラウン運動（Geometric Brownian Motion, GBM）です。
株価モデリングにおいて、通常のブラウン運動（算術ブラウン運動）をそのまま使うと、株価が負の値を取りうるという致命的な問題があります。幾何ブラウン運動は、ブラウン運動を指数関数に通すことでこの問題を解決し、常に正の値を取る確率過程を構成します。
幾何ブラウン運動を理解すると、以下のような分野で強力なツールを手にすることができます。
- 金融工学: ブラック・ショールズ・オプション価格理論の基盤であり、デリバティブの価格付けに不可欠です
- リスク管理: VaR（Value at Risk）やストレステストにおけるポートフォリオの変動モデルです
- 確率解析: 伊藤の公式の最も基本的かつ重要な応用例であり、確率微分方程式（SDE）の教科書的な例題です
- 生態学・人口動態: 環境変動下での個体群サイズのモデリングにも使われます
本記事の内容
- 幾何ブラウン運動の定義とモチベーション
- 確率微分方程式（SDE）の定式化
- 伊藤の公式による厳密解の導出
- 統計的性質（期待値・分散・対数正規分布）
- ブラック・ショールズモデルとの関係
- Pythonによるシミュレーションと統計的検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ブラウン運動の性質 — ブラウン運動の定義と基本性質
- 正規分布 — 対数正規分布の基礎
- ランダムウォークの理論 — 離散モデルから連続モデルへ
算術ブラウン運動の限界と幾何ブラウン運動の動機
なぜ通常のブラウン運動では不十分か
株価 $S_t$ を直接ブラウン運動でモデル化する — つまり $S_t = S_0 + \mu t + \sigma B_t$ とする — のは自然な発想ですが、いくつかの深刻な問題があります。
問題1: 負の株価: ブラウン運動は $-\infty$ から $+\infty$ まで値を取りうるので、$S_t < 0$ となる確率が正です。しかし、現実の株価は（上場廃止を除けば）負にはなりません。
問題2: 収益率の非定常性: 株価が100円の銘柄と10,000円の銘柄で、10円の変動は同じ意味でしょうか？100円の銘柄にとっては10%の変動ですが、10,000円の銘柄にとっては0.1%に過ぎません。株価の変動は、絶対値ではなく収益率（リターン）で測るべきです。
問題3: 複利効果: 実際の金融資産の成長は複利的であり、価格に比例して変動します。100万円の1%は1万円ですが、200万円の1%は2万円です。
これらの観察から、「株価の変動は現在の株価に比例する」というモデルが自然です。これを数式にしたものが幾何ブラウン運動です。
収益率モデルとしての定式化
株価 $S_t$ の微小変動 $dS_t$ が現在の株価 $S_t$ に比例するとしましょう。
$$ \frac{dS_t}{S_t} = \mu \, dt + \sigma \, dB_t $$
左辺は瞬時収益率（instantaneous return）です。右辺の第1項 $\mu \, dt$ は確定的なトレンド（ドリフト）を表し、第2項 $\sigma \, dB_t$ はランダムな変動（ボラティリティ）を表します。
- $\mu$: ドリフト率（drift rate）— 平均的な成長率。$\mu > 0$ なら平均的に上昇トレンド
- $\sigma$: ボラティリティ（volatility）— 変動の大きさ。$\sigma$ が大きいほどリスクが大きい
- $B_t$: 標準ブラウン運動 — ランダムなノイズ源
この方程式は、どの時点でも収益率が「平均 $\mu$、標準偏差 $\sigma$ のノイズ」で変動することを意味します。株価が高くても低くても、パーセンテージとしての変動特性は同じ — これが幾何ブラウン運動のモデルです。
整理すると、幾何ブラウン運動のSDEは次のように書けます。
$$ dS_t = \mu S_t \, dt + \sigma S_t \, dB_t $$
この方程式は、いわば通常のブラウン運動に「比例する」という条件を加えたものです。次に、この方程式を厳密に解きましょう。
伊藤の公式による厳密解の導出
伊藤の公式の復習
確率微分方程式 $dS_t = \mu S_t \, dt + \sigma S_t \, dB_t$ を解くために、伊藤の公式（Itô’s formula）を使います。
伊藤の公式は、ブラウン運動の関数に対するチェインルール（連鎖律）です。$X_t$ が $dX_t = a(t) \, dt + b(t) \, dB_t$ を満たし、$f(t, x)$ が十分になめらかな関数であるとき、
$$ df(t, X_t) = \left(\frac{\partial f}{\partial t} + a(t) \frac{\partial f}{\partial x} + \frac{1}{2} b(t)^2 \frac{\partial^2 f}{\partial x^2}\right) dt + b(t) \frac{\partial f}{\partial x} \, dB_t $$
通常の微積分のチェインルールとの違いは、$\frac{1}{2} b(t)^2 \frac{\partial^2 f}{\partial x^2}$ という2階微分の項が現れることです。これはブラウン運動の2次変分 $(dB_t)^2 = dt$ に起因しています。
解の導出
SDEを解くために、$Y_t = \ln S_t$ とおきます。$f(x) = \ln x$ に対して伊藤の公式を適用しましょう。
$f'(x) = 1/x$, $f”(x) = -1/x^2$ なので、$X_t = S_t$, $a(t) = \mu S_t$, $b(t) = \sigma S_t$ として、
$$ d(\ln S_t) = \frac{1}{S_t} dS_t + \frac{1}{2} \left(-\frac{1}{S_t^2}\right) (\sigma S_t)^2 dt $$
$dS_t = \mu S_t \, dt + \sigma S_t \, dB_t$ を代入します。
$$ d(\ln S_t) = \frac{1}{S_t}(\mu S_t \, dt + \sigma S_t \, dB_t) – \frac{1}{2} \sigma^2 dt $$
分数を整理すると、
$$ d(\ln S_t) = \mu \, dt + \sigma \, dB_t – \frac{1}{2} \sigma^2 dt $$
ドリフト項をまとめます。
$$ d(\ln S_t) = \left(\mu – \frac{\sigma^2}{2}\right) dt + \sigma \, dB_t $$
これは $\ln S_t$ が算術ブラウン運動（ドリフト付きブラウン運動）に従うことを意味します。右辺は $S_t$ に依存しない確定的な係数のみを含むので、両辺を $[0, t]$ で積分できます。
$$ \ln S_t – \ln S_0 = \left(\mu – \frac{\sigma^2}{2}\right) t + \sigma B_t $$
指数を取ると、
$$ S_t = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right) t + \sigma B_t\right] $$
これが幾何ブラウン運動の厳密解（closed-form solution）です。
解の構造の解釈
得られた解を詳しく見てみましょう。
$$ S_t = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right) t + \sigma B_t\right] $$
指数の中身は2つの部分からなります。
- 確定的トレンド: $\left(\mu – \frac{\sigma^2}{2}\right) t$ — 時間とともに線形に増加
- ランダム変動: $\sigma B_t$ — ブラウン運動によるノイズ
ここで注目すべきは、確定的トレンドが $\mu t$ ではなく $(\mu – \sigma^2/2)t$ であることです。$-\sigma^2/2$ という「補正項」はどこから来るのでしょうか？
これは伊藤の公式の2階微分項に由来します。対数を取ると $-\sigma^2/2$ のシフトが現れるのは、$\ln$ 関数が凹関数であることと関係しています。イェンセンの不等式の確率過程版と理解できます — $E[\ln S_t] < \ln E[S_t]$ なので、対数スケールでの平均成長率は、元のスケールでの平均成長率より $\sigma^2/2$ だけ小さくなります。
この補正項には実用上も重要な意味があります。$\mu = \sigma^2/2$ の場合、対数スケールでのドリフトがゼロになり、$\ln S_t$ はドリフトのないブラウン運動に従います。つまり、幾何ブラウン運動は $E[S_t]$ が指数的に増加していても、$\text{median}(S_t)$ は $S_0$ のまま変わらないという状況が起こりえます。
解が厳密に求まったので、次にこの解の統計的性質を詳しく調べていきましょう。
統計的性質
対数正規分布
厳密解から、$\ln S_t$ は正規分布に従うことがわかります。
$$ \ln S_t \sim N\left(\ln S_0 + \left(\mu – \frac{\sigma^2}{2}\right) t, \; \sigma^2 t\right) $$
したがって、$S_t$ は対数正規分布（log-normal distribution）に従います。対数正規分布は、対数を取ると正規分布になる分布です。その特徴として、常に正の値を取り、右に裾が長い非対称な形状をしています。
期待値と分散
対数正規分布のモーメントの公式を使って、$S_t$ の期待値と分散を計算しましょう。
$Z \sim N(m, v)$ に対して $E[e^Z] = e^{m + v/2}$ なので、
$$ E[S_t] = S_0 \cdot E\left[\exp\left(\left(\mu – \frac{\sigma^2}{2}\right)t + \sigma B_t\right)\right] $$
$\left(\mu – \frac{\sigma^2}{2}\right)t + \sigma B_t \sim N\left(\left(\mu – \frac{\sigma^2}{2}\right)t, \; \sigma^2 t\right)$ なので、
$$ E[S_t] = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right) t + \frac{\sigma^2 t}{2}\right] = S_0 e^{\mu t} $$
期待値は $S_0 e^{\mu t}$ という単純な指数関数です。$-\sigma^2/2$ の補正項とモーメント母関数からの $+\sigma^2/2$ が相殺し、ドリフト率 $\mu$ がそのまま平均成長率として現れます。
2次モーメントも同様に計算します。
$$ E[S_t^2] = S_0^2 \cdot E\left[\exp\left(2\left(\mu – \frac{\sigma^2}{2}\right)t + 2\sigma B_t\right)\right] $$
$2\left(\mu – \frac{\sigma^2}{2}\right)t + 2\sigma B_t \sim N\left(2\left(\mu – \frac{\sigma^2}{2}\right)t, \; 4\sigma^2 t\right)$ なので、
$$ E[S_t^2] = S_0^2 \exp\left[2\left(\mu – \frac{\sigma^2}{2}\right) t + 2\sigma^2 t\right] = S_0^2 e^{(2\mu + \sigma^2)t} $$
分散は、
$$ \text{Var}(S_t) = E[S_t^2] – (E[S_t])^2 = S_0^2 e^{2\mu t}\left(e^{\sigma^2 t} – 1\right) $$
分散は時間とともに指数的に増加し、その速さはボラティリティ $\sigma$ に依存します。$\sigma$ が大きいほど不確実性は急速に拡大します。
変動係数
分散の絶対値ではなく、期待値に対する相対的なばらつきを見るために、変動係数（coefficient of variation）を計算しましょう。
$$ \text{CV}(S_t) = \frac{\sqrt{\text{Var}(S_t)}}{E[S_t]} = \sqrt{e^{\sigma^2 t} – 1} $$
変動係数は $\mu$ に依存せず、$\sigma$ と $t$ のみで決まります。これは、不確実性の相対的な大きさがドリフトに依存しないことを意味します。
中央値とモード
対数正規分布の中央値（メディアン）は、
$$ \text{median}(S_t) = S_0 \exp\left[\left(\mu – \frac{\sigma^2}{2}\right) t\right] $$
モード（最頻値）は、
$$ \text{mode}(S_t) = S_0 \exp\left[\left(\mu – \frac{3\sigma^2}{2}\right) t\right] $$
3つの代表値の大小関係は、
$$ \text{mode}(S_t) < \text{median}(S_t) < E[S_t] $$
です（$\sigma > 0$ のとき）。対数正規分布は右に裾が長いため、平均は中央値より大きくなります。これは株式投資で重要な意味を持ちます — 「平均的には儲かる」（$E[S_t] > S_0$ for $\mu > 0$）が、「半分以上の場合で損をする」（$\text{median}(S_t) < S_0$ for $\mu < \sigma^2/2$）という状況がありえるのです。
対数収益率の分布
$[s, t]$ 間の対数収益率（log-return）は、
$$ r_{s,t} = \ln\frac{S_t}{S_s} = \left(\mu – \frac{\sigma^2}{2}\right)(t – s) + \sigma(B_t – B_s) $$
$B_t – B_s \sim N(0, t – s)$ なので、
$$ r_{s,t} \sim N\left(\left(\mu – \frac{\sigma^2}{2}\right)(t-s), \; \sigma^2(t-s)\right) $$
対数収益率が正規分布に従うという性質は、金融データの分析で頻繁に利用されます。ただし、実際の株価データでは裾の重い分布（ファットテール）が観察されることが多く、正規分布の仮定は近似に過ぎない点には注意が必要です。
ここまでで幾何ブラウン運動の統計的性質が明らかになりました。次に、この過程が金融工学でどのように使われるかを見ていきましょう。
ブラック・ショールズモデルとの関係
ブラック・ショールズの前提
1973年にフィッシャー・ブラックとマイロン・ショールズが発表したブラック・ショールズモデル（Black-Scholes model）は、ヨーロピアンオプションの理論価格を与える画期的な成果です。その基盤となっているのが幾何ブラウン運動による株価モデルです。
ブラック・ショールズモデルの主な仮定は以下の通りです。
- 株価は幾何ブラウン運動に従う: $dS_t = \mu S_t \, dt + \sigma S_t \, dB_t$
- 無リスク金利 $r$ は一定: 安全資産で $r$ の連続複利を得られる
- 取引コストなし: 売買手数料や税金がない
- 空売り可能: 株式を借りて売ることができる
- 裁定機会なし: リスクなしで確実に利益を得る方法は存在しない
リスク中立評価
ブラック・ショールズの理論で核心的な概念がリスク中立測度（risk-neutral measure）です。裁定なしの条件の下で、オプションの価格はリスク中立測度 $\mathbb{Q}$ の下での割引期待値として求められます。
リスク中立測度の下では、株価のドリフト率が $\mu$ から無リスク金利 $r$ に置き換わります。
$$ dS_t = r S_t \, dt + \sigma S_t \, d\tilde{B}_t \quad (\text{under } \mathbb{Q}) $$
ここで $\tilde{B}_t$ は $\mathbb{Q}$ の下での標準ブラウン運動です。この測度変換はギルサノフの定理（Girsanov’s theorem）によって正当化されます。
測度 $\mathbb{P}$（実測度）から $\mathbb{Q}$（リスク中立測度）への変換は、ラドン・ニコディム微分
$$ \frac{d\mathbb{Q}}{d\mathbb{P}} = \exp\left(-\frac{\mu – r}{\sigma} B_T – \frac{1}{2}\left(\frac{\mu – r}{\sigma}\right)^2 T\right) $$
で定義されます。ここで $\lambda = (\mu – r)/\sigma$ はマーケットのリスクの市場価格（market price of risk）と呼ばれます。
ブラック・ショールズ公式
ヨーロピアン・コールオプション（満期 $T$ で行使価格 $K$ で株を買う権利）の理論価格は、
$$ C = S_0 \Phi(d_1) – K e^{-rT} \Phi(d_2) $$
ここで、
$$ d_1 = \frac{\ln(S_0/K) + (r + \sigma^2/2)T}{\sigma\sqrt{T}}, \quad d_2 = d_1 – \sigma\sqrt{T} $$
$\Phi$ は標準正規分布の累積分布関数です。
この公式は、リスク中立測度 $\mathbb{Q}$ の下で $C = e^{-rT} E_{\mathbb{Q}}[\max(S_T – K, 0)]$ を計算することで導出されます。$S_T$ がリスク中立測度の下で対数正規分布に従うことが鍵となっています。
ブラック・ショールズ公式の導出
リスク中立測度の下で $S_T$ の解は、
$$ S_T = S_0 \exp\left[\left(r – \frac{\sigma^2}{2}\right)T + \sigma \tilde{B}_T\right] $$
$Z = \tilde{B}_T / \sqrt{T} \sim N(0, 1)$ とおくと、
$$ S_T = S_0 \exp\left[\left(r – \frac{\sigma^2}{2}\right)T + \sigma\sqrt{T} Z\right] $$
コールオプションの期待ペイオフは、
$$ E_{\mathbb{Q}}[\max(S_T – K, 0)] = \int_{-\infty}^{\infty} \max(S_T(z) – K, 0) \cdot \frac{1}{\sqrt{2\pi}} e^{-z^2/2} dz $$
$S_T(z) > K$ となるのは $z > z^* = -d_2$ のときです。ここで $z^*$ を求めます。$S_T(z) = K$ とおくと、
$$ S_0 \exp\left[\left(r – \frac{\sigma^2}{2}\right)T + \sigma\sqrt{T} z\right] = K $$
対数を取って $z$ について解くと、
$$ z = \frac{\ln(K/S_0) – (r – \sigma^2/2)T}{\sigma\sqrt{T}} = -d_2 $$
よって、
$$ E_{\mathbb{Q}}[\max(S_T – K, 0)] = \int_{-d_2}^{\infty} (S_T(z) – K) \frac{1}{\sqrt{2\pi}} e^{-z^2/2} dz $$
この積分を $S_T(z)$ の部分と $K$ の部分に分けると、
$$ = S_0 e^{rT} \Phi(d_1) – K \Phi(d_2) $$
（$S_0 e^{rT}$ の部分は平方完成によって得られます。）割引すると、
$$ C = e^{-rT}\left[S_0 e^{rT} \Phi(d_1) – K \Phi(d_2)\right] = S_0 \Phi(d_1) – K e^{-rT} \Phi(d_2) $$
これがブラック・ショールズ公式です。
幾何ブラウン運動の限界
実務では、幾何ブラウン運動モデルには以下の限界があることが知られています。
- ファットテール: 実際のリターンの分布は正規分布よりも裾が重い（1987年のブラックマンデーのような極端な変動が正規分布の予測よりはるかに高い頻度で発生）
- ボラティリティ・スマイル: オプション市場で観察されるインプライド・ボラティリティの行使価格依存性は、一定ボラティリティの仮定と矛盾する
- ジャンプ: 企業の破綻や重大ニュースにより、株価が不連続に変動することがある
- ボラティリティ・クラスタリング: ボラティリティは一定ではなく、高い変動が高い変動を呼ぶ傾向がある
これらの限界を克服するために、確率的ボラティリティモデル（Heston model）、ジャンプ拡散モデル（Merton model）、局所ボラティリティモデル（Dupire model）などの拡張が提案されています。しかし、幾何ブラウン運動はその解析的な扱いやすさから、今なお金融工学の基盤として広く使われています。
理論を十分に理解したところで、Pythonを使って幾何ブラウン運動をシミュレーションし、理論的な性質を視覚的に確認してみましょう。
Pythonによるシミュレーション
基本的なパスの生成
幾何ブラウン運動のサンプルパスを厳密解を使って生成します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# パラメータ
S0 = 100 # 初期株価
mu = 0.08 # ドリフト率（年率8%）
sigma = 0.2 # ボラティリティ（年率20%）
T = 2.0 # 期間（2年）
N = 500 # ステップ数
n_paths = 10 # パス数
dt = T / N
t = np.linspace(0, T, N + 1)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: GBMのサンプルパス
ax = axes[0]
for i in range(n_paths):
dB = np.random.normal(0, np.sqrt(dt), N)
B = np.concatenate([[0], np.cumsum(dB)])
S = S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * B)
ax.plot(t, S, linewidth=0.8, alpha=0.7)
# 期待値の曲線
ax.plot(t, S0 * np.exp(mu * t), 'k--', linewidth=2, label=f"E[S(t)] = {S0}exp({mu}t)")
ax.set_xlabel("Time (years)")
ax.set_ylabel("Stock Price")
ax.set_title("Sample Paths of Geometric Brownian Motion")
ax.legend()
ax.grid(True, alpha=0.3)
# 右: 対数スケール
ax = axes[1]
np.random.seed(42)
for i in range(n_paths):
dB = np.random.normal(0, np.sqrt(dt), N)
B = np.concatenate([[0], np.cumsum(dB)])
S = S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * B)
ax.plot(t, S, linewidth=0.8, alpha=0.7)
ax.plot(t, S0 * np.exp(mu * t), 'k--', linewidth=2, label=f"E[S(t)]")
ax.set_yscale('log')
ax.set_xlabel("Time (years)")
ax.set_ylabel("Stock Price (log scale)")
ax.set_title("GBM Paths (Log Scale)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、10本のGBMパスが初期価格100から出発して多様な軌跡を描いています。黒い破線は期待値 $E[S_t] = S_0 e^{\mu t}$ を表しています。パスは上方にも下方にも大きく変動し、分布の非対称性（右裾が長い）が視覚的にわかります。右のパネルは対数スケールで同じパスを表示したものです。対数スケールでは、パスがブラウン運動に近い形（線形トレンド + ランダム変動）になっていることが確認でき、$\ln S_t$ がドリフト付きブラウン運動に従うという理論と整合しています。
終端分布の検証
多数のパスをシミュレーションし、終端時刻 $T$ での $S_T$ の分布が理論的な対数正規分布と一致することを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
S0, mu, sigma, T = 100, 0.08, 0.2, 1.0
n_sim = 100000
# 終端値のシミュレーション
Z = np.random.normal(0, 1, n_sim)
S_T = S0 * np.exp((mu - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
# 理論分布のパラメータ
log_mean = np.log(S0) + (mu - 0.5 * sigma**2) * T
log_std = sigma * np.sqrt(T)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: S_Tのヒストグラムと対数正規分布
ax = axes[0]
ax.hist(S_T, bins=100, density=True, alpha=0.7, color='steelblue', label='Simulation')
x = np.linspace(50, 200, 500)
pdf = stats.lognorm.pdf(x, s=log_std, scale=np.exp(log_mean))
ax.plot(x, pdf, 'r-', linewidth=2, label='Log-normal PDF')
ax.axvline(np.mean(S_T), color='green', linestyle='--', linewidth=1.5,
label=f'Mean = {np.mean(S_T):.2f}')
ax.axvline(np.median(S_T), color='orange', linestyle='--', linewidth=1.5,
label=f'Median = {np.median(S_T):.2f}')
ax.set_xlabel("S(T)")
ax.set_ylabel("Density")
ax.set_title(f"Distribution of S(T) at T = {T}")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 右: ln(S_T)のヒストグラムと正規分布
ax = axes[1]
log_S_T = np.log(S_T)
ax.hist(log_S_T, bins=100, density=True, alpha=0.7, color='steelblue', label='Simulation')
x_log = np.linspace(log_mean - 4*log_std, log_mean + 4*log_std, 500)
pdf_log = stats.norm.pdf(x_log, loc=log_mean, scale=log_std)
ax.plot(x_log, pdf_log, 'r-', linewidth=2, label='Normal PDF')
ax.set_xlabel("ln S(T)")
ax.set_ylabel("Density")
ax.set_title(f"Distribution of ln S(T)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 統計量の比較
print("統計量の比較:")
print(f" E[S_T] 理論: {S0 * np.exp(mu * T):.4f}, "
f"シミュレーション: {np.mean(S_T):.4f}")
print(f" 中央値 理論: {S0 * np.exp((mu - 0.5*sigma**2) * T):.4f}, "
f"シミュレーション: {np.median(S_T):.4f}")
print(f" Var[S_T] 理論: {S0**2 * np.exp(2*mu*T) * (np.exp(sigma**2*T) - 1):.4f}, "
f"シミュレーション: {np.var(S_T):.4f}")
左のパネルで、シミュレーションのヒストグラムが対数正規分布の理論曲線（赤）と非常によく一致しています。平均（緑の破線）が中央値（橙の破線）より大きいことが確認でき、これは対数正規分布の右裾の重さを反映しています。右のパネルでは、$\ln S_T$ のヒストグラムが正規分布と一致しており、$\ln S_T \sim N(\ln S_0 + (\mu – \sigma^2/2)T, \sigma^2 T)$ という理論的予測が検証されています。
ブラック・ショールズ公式の数値検証
モンテカルロ法によるオプション価格のシミュレーション結果と、ブラック・ショールズ公式の理論値を比較します。
import numpy as np
from scipy import stats
np.random.seed(42)
# パラメータ
S0 = 100 # 現在の株価
K = 105 # 行使価格
r = 0.05 # 無リスク金利
sigma = 0.2 # ボラティリティ
T = 1.0 # 満期（1年）
n_sim = 1000000
# ブラック・ショールズ公式
d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
BS_call = S0 * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)
BS_put = K * np.exp(-r * T) * stats.norm.cdf(-d2) - S0 * stats.norm.cdf(-d1)
# モンテカルロ・シミュレーション（リスク中立測度の下で）
Z = np.random.normal(0, 1, n_sim)
S_T = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
call_payoffs = np.maximum(S_T - K, 0)
put_payoffs = np.maximum(K - S_T, 0)
MC_call = np.exp(-r * T) * np.mean(call_payoffs)
MC_put = np.exp(-r * T) * np.mean(put_payoffs)
# 標準誤差
MC_call_se = np.exp(-r * T) * np.std(call_payoffs) / np.sqrt(n_sim)
MC_put_se = np.exp(-r * T) * np.std(put_payoffs) / np.sqrt(n_sim)
print("ブラック・ショールズ vs モンテカルロ:")
print(f" コールオプション: BS = {BS_call:.4f}, MC = {MC_call:.4f} ± {MC_call_se:.4f}")
print(f" プットオプション: BS = {BS_put:.4f}, MC = {MC_put:.4f} ± {MC_put_se:.4f}")
print(f" プット・コール・パリティ: C - P = {BS_call - BS_put:.4f}, "
f"S0 - K*exp(-rT) = {S0 - K * np.exp(-r * T):.4f}")
ブラック・ショールズの理論値とモンテカルロ・シミュレーションの結果が、標準誤差の範囲内で一致することが確認できます。また、プット・コール・パリティ $C – P = S_0 – Ke^{-rT}$ が成立していることも検証されています。100万パスのシミュレーションにより、理論値との差は小数点以下2桁で一致するレベルの精度が得られています。
パラメータの影響の可視化
ドリフト率 $\mu$ とボラティリティ $\sigma$ が幾何ブラウン運動の挙動にどのような影響を与えるかを可視化します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
S0 = 100
T = 2.0
N = 500
dt = T / N
t = np.linspace(0, T, N + 1)
n_paths = 1000
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# ボラティリティの影響
configs = [
(0.08, 0.05, "Low vol (σ=0.05)"),
(0.08, 0.20, "Medium vol (σ=0.20)"),
(0.08, 0.50, "High vol (σ=0.50)"),
]
for idx, (mu, sig, label) in enumerate(configs):
ax = axes[0, 0] if idx == 0 else (axes[0, 1] if idx == 1 else axes[1, 0])
percentiles = np.zeros((N + 1, 5))
paths = np.zeros((n_paths, N + 1))
for i in range(n_paths):
dB = np.random.normal(0, np.sqrt(dt), N)
B = np.concatenate([[0], np.cumsum(dB)])
paths[i] = S0 * np.exp((mu - 0.5 * sig**2) * t + sig * B)
for j, p in enumerate([5, 25, 50, 75, 95]):
percentiles[:, j] = np.percentile(paths, p, axis=0)
ax.fill_between(t, percentiles[:, 0], percentiles[:, 4],
alpha=0.2, color='steelblue', label='5-95%')
ax.fill_between(t, percentiles[:, 1], percentiles[:, 3],
alpha=0.4, color='steelblue', label='25-75%')
ax.plot(t, percentiles[:, 2], 'b-', linewidth=1.5, label='Median')
ax.plot(t, S0 * np.exp(mu * t), 'r--', linewidth=1.5, label='E[S(t)]')
ax.set_xlabel("Time (years)")
ax.set_ylabel("S(t)")
ax.set_title(label)
ax.legend(fontsize=8)
ax.set_ylim(0, 400)
ax.grid(True, alpha=0.3)
# 平均 vs 中央値の乖離
ax = axes[1, 1]
sigmas = np.linspace(0.01, 0.8, 100)
T_val = 1.0
mean_ratio = np.exp(0.08 * T_val) * np.ones_like(sigmas)
median_ratio = np.exp((0.08 - 0.5 * sigmas**2) * T_val)
ax.plot(sigmas, mean_ratio, 'r-', linewidth=2, label='E[S(T)] / S₀')
ax.plot(sigmas, median_ratio, 'b-', linewidth=2, label='Median(S(T)) / S₀')
ax.axhline(y=1.0, color='k', linestyle='--', linewidth=0.5)
ax.fill_between(sigmas, median_ratio, mean_ratio, alpha=0.2, color='purple')
ax.set_xlabel("Volatility σ")
ax.set_ylabel("S(T) / S₀")
ax.set_title("Mean vs Median at T = 1")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上段と左下のパネルは、ボラティリティの異なる3つの設定での幾何ブラウン運動のパーセンタイル帯を示しています。$\sigma = 0.05$（低ボラティリティ）では分布が狭く集中し、中央値と期待値はほぼ一致します。$\sigma = 0.50$（高ボラティリティ）では分布が大きく広がり、中央値（青実線）が期待値（赤破線）から大きく下方に乖離しています。
右下のパネルは、$T = 1$ における $\sigma$ の関数としての平均と中央値の比較です。平均 $E[S_T]/S_0 = e^{\mu T}$ は $\sigma$ に依存しませんが、中央値 $\text{median}(S_T)/S_0 = e^{(\mu – \sigma^2/2)T}$ は $\sigma$ が大きくなるほど低下します。$\sigma \approx 0.4$ 付近で中央値が1を下回り、「半数以上の場合で初期値より低くなる」状態になります。これは「平均で見れば成長するが、典型的な結果は損失」というパラドックスの視覚化です。
幾何ブラウン運動のオイラー・丸山法によるシミュレーション
厳密解を使わず、SDEを直接数値的に解くオイラー・丸山法（Euler-Maruyama method）でもシミュレーションできます。厳密解との比較で数値スキームの精度を確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
S0, mu, sigma, T = 100, 0.08, 0.2, 1.0
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: ステップ数による精度比較
ax = axes[0]
Ns = [10, 50, 200, 1000]
# 共通のブラウン運動を最細かい刻みで生成
N_fine = 10000
dt_fine = T / N_fine
dB_fine = np.random.normal(0, np.sqrt(dt_fine), N_fine)
B_fine = np.concatenate([[0], np.cumsum(dB_fine)])
t_fine = np.linspace(0, T, N_fine + 1)
# 厳密解
S_exact = S0 * np.exp((mu - 0.5 * sigma**2) * t_fine + sigma * B_fine)
ax.plot(t_fine, S_exact, 'k-', linewidth=1.5, label='Exact', alpha=0.8)
colors = ['red', 'blue', 'green', 'orange']
for N_em, color in zip(Ns, colors):
step = N_fine // N_em
dt_em = T / N_em
t_em = np.linspace(0, T, N_em + 1)
# オイラー・丸山法
S_em = np.zeros(N_em + 1)
S_em[0] = S0
for i in range(N_em):
# 粗い刻みのブラウン運動増分を合成
dB_em = np.sum(dB_fine[i*step:(i+1)*step])
S_em[i+1] = S_em[i] + mu * S_em[i] * dt_em + sigma * S_em[i] * dB_em
ax.plot(t_em, S_em, '--', color=color, linewidth=1, alpha=0.8,
label=f'Euler-Maruyama (N={N_em})')
ax.set_xlabel("Time")
ax.set_ylabel("S(t)")
ax.set_title("Euler-Maruyama vs Exact Solution")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 右: 強収束の確認
ax = axes[1]
Ns_test = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000]
n_trials = 200
errors = []
for N_em in Ns_test:
step = N_fine // N_em
dt_em = T / N_em
err_list = []
for trial in range(n_trials):
dB_f = np.random.normal(0, np.sqrt(dt_fine), N_fine)
B_f = np.concatenate([[0], np.cumsum(dB_f)])
S_ex = S0 * np.exp((mu - 0.5*sigma**2)*T + sigma*B_f[-1])
S_em_val = S0
for i in range(N_em):
dB_em = np.sum(dB_f[i*step:(i+1)*step])
S_em_val = S_em_val + mu*S_em_val*dt_em + sigma*S_em_val*dB_em
err_list.append(abs(S_em_val - S_ex))
errors.append(np.mean(err_list))
dts = [T/n for n in Ns_test]
ax.loglog(dts, errors, 'bo-', markersize=6, label='Euler-Maruyama error')
# 理論的な収束次数 0.5 の参照線
ref = errors[0] * (np.array(dts) / dts[0])**0.5
ax.loglog(dts, ref, 'r--', linewidth=1.5, label='O(Δt^0.5) reference')
ax.set_xlabel("Δt")
ax.set_ylabel("Mean |S_exact - S_EM| at T")
ax.set_title("Strong Convergence of Euler-Maruyama")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、異なるステップ数でのオイラー・丸山法の近似が厳密解に収束していく様子が確認できます。$N = 10$ では大きな乖離がありますが、$N = 1000$ では厳密解とほぼ重なります。
右のパネルは強収束（strong convergence）の検証です。オイラー・丸山法の強収束次数は $O(\Delta t^{0.5})$ であることが理論的に知られており、対数スケールでの傾きが赤い参照線（傾き0.5）とほぼ一致しています。ミルシュタイン法を使えば $O(\Delta t)$ の強収束が得られますが、GBMの場合は厳密解が利用できるため、実用上はオイラー・丸山法で十分です。
幾何ブラウン運動の拡張
時間依存パラメータ
実際の金融市場では、ドリフト率やボラティリティが時間とともに変化します。幾何ブラウン運動を時間依存パラメータに拡張することができます。
$$ dS_t = \mu(t) S_t \, dt + \sigma(t) S_t \, dB_t $$
この場合の解は、
$$ S_t = S_0 \exp\left[\int_0^t \left(\mu(s) – \frac{\sigma(s)^2}{2}\right) ds + \int_0^t \sigma(s) \, dB_s\right] $$
$\sigma(t)$ が確定的な関数の場合、$\int_0^t \sigma(s) \, dB_s$ は正規分布 $N(0, \int_0^t \sigma(s)^2 ds)$ に従うので、解は依然として対数正規分布です。
多次元への拡張
$d$ 個の資産の株価 $(S_t^{(1)}, \dots, S_t^{(d)})$ を同時にモデル化する場合、相関のあるブラウン運動を使います。
$$ dS_t^{(i)} = \mu_i S_t^{(i)} \, dt + \sigma_i S_t^{(i)} \, dB_t^{(i)} $$
ここで $\text{Cov}(dB_t^{(i)}, dB_t^{(j)}) = \rho_{ij} \, dt$ です。$\rho_{ij}$ は資産間の相関係数であり、ポートフォリオ理論やバスケットオプションの価格付けにおいて重要です。
確率的ボラティリティモデル
ボラティリティ自体をも確率過程でモデル化する確率的ボラティリティモデルは、幾何ブラウン運動の自然な拡張です。代表的なHestonモデルは、
$$ dS_t = \mu S_t \, dt + \sqrt{V_t} S_t \, dB_t^{(1)} $$
$$ dV_t = \kappa(\theta – V_t) \, dt + \xi \sqrt{V_t} \, dB_t^{(2)} $$
ボラティリティ $V_t$ がオルンシュタイン・ウーレンベック型の平均回帰過程に従い、$B_t^{(1)}$ と $B_t^{(2)}$ は相関 $\rho$ を持つブラウン運動です。この拡張により、ボラティリティ・スマイルを再現できます。
まとめ
本記事では、幾何ブラウン運動（GBM）の理論を体系的に解説しました。
- 定義と動機: 収益率が正規分布に従うモデルとして、$dS_t = \mu S_t \, dt + \sigma S_t \, dB_t$ を定式化
- 厳密解の導出: 伊藤の公式を用いて $S_t = S_0 \exp[(\mu – \sigma^2/2)t + \sigma B_t]$ を導出。$-\sigma^2/2$ の補正項が伊藤補正として本質的
- 統計的性質: $S_t$ は対数正規分布に従い、$E[S_t] = S_0 e^{\mu t}$, $\text{median}(S_t) = S_0 e^{(\mu – \sigma^2/2)t}$
- ブラック・ショールズ: リスク中立測度の下でGBMの厳密解を利用し、オプション価格の閉じた公式を導出
- 数値シミュレーション: オイラー・丸山法の強収束次数 $O(\Delta t^{0.5})$ を数値的に確認
幾何ブラウン運動は、ブラウン運動と指数関数という2つの基本要素を組み合わせた美しいモデルであり、確率解析の入り口として理想的です。実務的にはモデルの限界を認識しつつも、その解析的な扱いやすさは他の発展的モデルの出発点として今なお重要です。
次のステップとして、以下の記事も参考にしてください。
- オルンシュタイン・ウーレンベック過程 — 平均回帰モデル
- 伊藤積分の定義と性質 — GBMの解法で用いた伊藤積分の基礎
- 伊藤の公式の応用 — SDEの求解テクニック