March 6, 2026

qkit5: Stochastic Volatility and Monte Carlo Pricing

Hestonモデル、フーリエ変換プライシング、モンテカルロ法、そしてラフ・ボラティリティ

pt.1〜4では、BSMの定数ボラティリティ仮定の下で価格計算とリスク管理をやってきた。pt.5ではこの仮定を根本から覆す。ボラティリティ自体が確率過程に従うモデルを扱う。

Heston (1993) の確率的ボラティリティモデルは、IVサーフェスのスマイル/スキューを自然に再現して、特性関数の解析解を持つから、フーリエ変換（FFT）で高速にプライシングできる。解析解を持たないエキゾチックオプションにはモンテカルロ・シミュレーション（MC）を使う。最後に、2014年以降の最大のパラダイムシフトであるラフ・ボラティリティへの接続も整理しておく。

1. Hestonモデル：確率的ボラティリティ

Heston (1993) は原資産価格 $S_t$ と瞬間分散 $V_t$ の2因子モデルを提案した。

dS_t = \mu\, S_t\, dt + \sqrt{V_t}\, S_t\, dW_t^S

dV_t = \kappa\,(\theta - V_t)\, dt + \xi\, \sqrt{V_t}\, dW_t^V

\mathrm{Corr}(dW_t^S,\, dW_t^V) = \rho\, dt

5つのパラメータがある。 $V_0$ は初期瞬間分散（SPXだと $0.04$ くらい、 $\sigma_0 = 20\%$ ）。 $\kappa$ は平均回帰速度（ $1.0$ 〜 $5.0$ ）。 $\theta$ は長期平均分散（ $0.04$ 〜 $0.09$ ）。 $\xi$ はvol-of-vol、ボラティリティの変動性（ $0.2$ 〜 $0.8$ ）。 $\rho$ は株価-ボラティリティ相関（ $-0.9$ 〜 $-0.5$ ）。

パラメータの直感的理解

平均回帰（ $\kappa, \theta$ ）について。 $V_t$ は長期平均 $\theta$ に向かって速度 $\kappa$ で引き戻される。 $\kappa$ が大きいほど回帰が速い。分散の半減期は $\tau_{1/2} = \ln 2 / \kappa$ 。

Vol-of-vol（ $\xi$ ）は、 $\xi$ が大きいほど分散自体の変動が激しくなって、IVサーフェスのスマイルが急峻になる（ファットテール効果）。

相関（ $\rho$ ）は、 $\rho < 0$ がレバレッジ効果を捕捉する。株価下落時にボラティリティが上昇して、IVサーフェスにスキューを生む。

Feller条件として $2\kappa\theta > \xi^2$ のとき、 $V_t$ はゼロに到達しない（a.s.）。だけど実データのキャリブレーションではしばしば破られる。

リスク中立測度Qのダイナミクス

オプション価格付けには $\mathbb{Q}$ の下でのダイナミクスを使う。ボラティリティのリスクの市場価格に依存するけど、アフィン構造を保つために標準的な仮定を置くと、

dS_t = r\, S_t\, dt + \sqrt{V_t}\, S_t\, dW_t^{S,\mathbb{Q}}

dV_t = \kappa^*(\theta^* - V_t)\, dt + \xi\, \sqrt{V_t}\, dW_t^{V,\mathbb{Q}}

$\kappa^* = \kappa + \lambda$ 、 $\theta^* = \kappa\theta / (\kappa + \lambda)$ 。 $\xi$ と $\rho$ は測度変換で不変。

なぜHestonがIVスマイルを生むか

BSMの仮定（ $\sigma$ 一定）だとIVはフラット。Hestonではボラティリティ自体が確率的に変動するから、 $\xi > 0$ でリターン分布がファットテールになってスマイル（OTM側のIV上昇）が生まれる。 $\rho < 0$ でレバレッジ効果からスキュー（プット側のIVがコール側より高い）が出る。 $\kappa, \theta$ はタームストラクチャーの形状を決定する。

2. Hestonの特性関数

Hestonモデルのアフィン構造から、対数株価 $x_T = \ln S_T$ の特性関数が閉形式で得られる。

\phi(u) = \mathbb{E}^{\mathbb{Q}}\!\left[e^{iu\, x_T}\right] = e^{C(u,\tau) + D(u,\tau)\, V_0 + iu\, x_0}

$\tau = T - t$ は残存期間、 $x_0 = \ln S_0 + r\tau$ はフォワード対数株価。

$C(u, \tau)$ と $D(u, \tau)$ の解析解は、

C(u, \tau) = \frac{\kappa^* \theta^*}{\xi^2}\left[(\kappa^* - \rho \xi i u - d)\tau - 2\ln\left(\frac{1 - g e^{-d\tau}}{1 - g}\right)\right]

D(u, \tau) = \frac{\kappa^* - \rho\xi\, iu - d}{\xi^2} \cdot \frac{1 - e^{-d\tau}}{1 - g\, e^{-d\tau}}

補助変数は、

d = \sqrt{(\rho\xi\, iu - \kappa^*)^2 + \xi^2(iu + u^2)}

g = \frac{\kappa^* - \rho\xi\, iu - d}{\kappa^* - \rho\xi\, iu + d}

Little Heston Trap

上の標準的な定式化は、 $d$ の複素平方根の分岐切断（branch cut）で数値的に不安定になることがある（Albrecher et al., 2006）。回避策として $C, D$ を等価形式に書き換える。

C(u, \tau) = \frac{\kappa^* \theta^*}{\xi^2}\left[(\kappa^* - \rho \xi i u + d)\tau - 2\ln\left(\frac{1 - g^{-1} e^{d\tau}}{1 - g^{-1}}\right)\right]

D(u, \tau) = \frac{\kappa^* - \rho\xi\, iu + d}{\xi^2} \cdot \frac{1 - e^{d\tau}}{1 - g^{-1} e^{d\tau}}

$g$ を $g^{-1}$ に置き換えて $d$ の符号を反転させたもの。数値的に安定だから、実装では常にこっちを使う。

3. フーリエ変換プライシング

Gil-Pelaez反転

ヨーロピアン・コール価格はHeston (1993) の原論文の方法で計算できる。

C = S_0\, P_1 - K\, e^{-r\tau}\, P_2

$P_1, P_2$ はフーリエ反転で計算される行使確率で、

P_j = \frac{1}{2} + \frac{1}{\pi}\int_0^\infty \mathrm{Re}\!\left[\frac{e^{-iu\ln K}\, \phi_j(u)}{iu}\right] du, \qquad j = 1, 2

$\phi_1, \phi_2$ はそれぞれ異なるニューメレール（ $S$ 測度と $B$ 測度）の下での特性関数。

Carr-Madan FFT法

Carr & Madan (1999) は上記を改良して、FFTで一度に全ストライクのコール価格を計算する手法を提案した。

対数ストライク $k = \ln K$ のコール価格 $C(k)$ は $k \to -\infty$ で $S_0$ に漸近して $L^1$ 可積分じゃない。ダンピング因子 $e^{\alpha k}$ （ $\alpha > 0$ ）を導入して可積分にする。

c_\alpha(k) = e^{\alpha k}\, C(k)

この修正コール価格のフーリエ変換は、

\hat{c}_\alpha(u) = \frac{e^{-r\tau}\, \phi(u - (\alpha+1)i)}{(\alpha + iu)(\alpha + 1 + iu)}

コール価格はフーリエ逆変換で復元。

C(k) = \frac{e^{-\alpha k}}{\pi}\int_0^\infty e^{-iuk}\, \hat{c}_\alpha(u)\, du

FFTの適用条件は $\Delta k \cdot \Delta u = 2\pi / N$ 。 $N = 4096$ のとき $\Delta k \approx 0.006$ （約0.6%刻み）で実用上十分。

計算効率について。Gil-Pelaez法は各ストライクごとに積分が必要（ $O(N_K \cdot N_u)$ ）だけど、Carr-Madan FFTは $O(N \log N)$ で全ストライクを同時に計算。キャリブレーション（数百回の再計算が必要）では桁違いに速い。

ダンピング因子 $\alpha$ はMGFの存在域内でなければならない。実用的には $\alpha = 1.5$ が標準（Carr & Madan, 1999）。小さすぎると被積分関数の振動が激しくて、大きすぎると数値的に不安定。

COS法（Fang & Oosterlee, 2008）はCarr-Madanの代替で、特性関数をフーリエ余弦級数で展開する。ダンピング因子のチューニングが不要で収束が指数的。

4. モンテカルロ・シミュレーション

GBMのオイラー離散化

リスク中立測度 $\mathbb{Q}$ の下で、GBM $dS = rS\,dt + \sigma S\,dW$ の離散化は、

S_{t+\Delta t} = S_t \exp\!\left[\left(r - \frac{1}{2}\sigma^2\right)\Delta t + \sigma\sqrt{\Delta t}\, Z\right], \qquad Z \sim \mathcal{N}(0,1)

対数をとって離散化する（対数オイラー法）。ドリフト補正項 $-\sigma^2/2$ を忘れるとマルチンゲール性が崩れるから注意。

Hestonモデルの離散化

2つの相関するブラウン運動を生成するのにコレスキー分解を使う。

Z_t^S = Z_1, \qquad Z_t^V = \rho\, Z_1 + \sqrt{1 - \rho^2}\, Z_2

$Z_1, Z_2 \sim \mathcal{N}(0,1)$ は独立。

オイラー法（最も単純だけどFeller条件が破られると $V_t < 0$ になりうる）は、

V_{t+\Delta t} = V_t + \kappa^*(\theta^* - V_t)\,\Delta t + \xi\sqrt{\max(V_t, 0)}\,\sqrt{\Delta t}\, Z_t^V

\ln S_{t+\Delta t} = \ln S_t + \left(r - \frac{1}{2}\max(V_t, 0)\right)\Delta t + \sqrt{\max(V_t, 0)}\,\sqrt{\Delta t}\, Z_t^S

$\max(V_t, 0)$ はフルトランケーション（Full Truncation）スキームで、 $V_t$ が負になったらゼロにクランプする。Lord, Koekkoek & van Dijk (2010) はこのスキームがバイアスと分散の両面で最も優れていることを示した。

実務標準はAndersen (2008) のQE（Quadratic Exponential）スキーム。 $V_t$ の条件付き分布が非心カイ二乗分布に従うことを利用して、 $V_t$ が小さいときは指数分布近似、大きいときは二次多項式近似。オイラー法よりはるかに正確で、タイムステップ数を大幅に削減できる。

オプション価格は $N$ 本のパスの満期ペイオフの平均。

\hat{C} = e^{-r\tau} \cdot \frac{1}{N}\sum_{i=1}^{N} \max\!\left(S_\tau^{(i)} - K,\, 0\right)

MCの標準誤差は $\text{SE} = \hat{\sigma}_{\text{payoff}} / \sqrt{N}$ 。95%信頼区間は $\hat{C} \pm 1.96 \cdot \text{SE}$ 。 $N = 10^5$ で通常十分な精度。

5. 分散削減技法

対称変量法（Antithetic Variates）

$Z$ と $-Z$ の両方のパスを生成して平均を取る。

\hat{C}_{\text{AV}} = \frac{1}{2}\left(\hat{C}(\{Z_t\}) + \hat{C}(\{-Z_t\})\right)

コール/プットのペイオフだと実効的に分散が約50%削減される。計算コストの追加はほぼゼロ。

制御変量法（Control Variates）

解析解が既知のBSM価格を制御変量として使う。

\hat{C}_{\text{CV}} = \hat{C}_{\text{Heston}}^{\text{MC}} - \beta\left(\hat{C}_{\text{BSM}}^{\text{MC}} - C_{\text{BSM}}^{\text{exact}}\right)

最適 $\beta$ は $\beta^* = \text{Cov}(\hat{C}_{\text{Heston}}, \hat{C}_{\text{BSM}}) / \text{Var}(\hat{C}_{\text{BSM}})$ 。HestonとBSMの相関が高い場合（特にATM付近）、分散を90%以上削減可能。けっこうすごい。

準モンテカルロ（Quasi-Monte Carlo）

擬似乱数の代わりにソボル列（Sobol sequence）みたいな低食い違い量列を使う。QMCの収束速度は $O((\log N)^d / N)$ で、標準MCの $O(1/\sqrt{N})$ より高次元でも速い。

だけどQMCは次元 $d$ （= タイムステップ数 × 因子数）が大きいと効率が落ちる。主成分分析（PCA）によるブラウニアン・ブリッジで有効次元を削減するのが標準。

重点サンプリング（Importance Sampling）

Deep OTMオプション（ペイオフがゼロになるパスが大部分）で有効。ドリフトを変更してITM領域から集中的にサンプリングして、ラドン・ニコディム微分（尤度比）で補正する。

6. エキゾチックオプション

解析解が存在しないから、MCが事実上の唯一の汎用的手法。

バリアオプションは原資産がバリア価格 $B$ に到達すると消滅（knock-out）または発生（knock-in）する。Down-and-Outコール（ $B < S_0$ ）は、

C_{\text{DO}} = e^{-r\tau} \cdot \mathbb{E}^{\mathbb{Q}}\!\left[\max(S_\tau - K, 0) \cdot \mathbb{1}_{\{\min_{0 \leq t \leq \tau} S_t > B\}}\right]

MCでは離散タイムステップで $S_t$ を観測するけど、連続バリアとの間にバイアスが生じる。Broadie, Glasserman & Kou (1997) の補正で、

B^* = B \cdot e^{\beta\sigma\sqrt{\Delta t}}

$\beta \approx 0.5826$ は普遍定数。

アジアンオプションは平均価格で決済する。算術平均コールは、

C_{\text{Asian}} = e^{-r\tau} \cdot \mathbb{E}^{\mathbb{Q}}\!\left[\max\!\left(\frac{1}{n}\sum_{i=1}^{n} S_{t_i} - K,\; 0\right)\right]

算術平均の分布には閉形式がない（対数正規分布の和は対数正規じゃない）から、MCが標準手法。幾何平均アジアンには閉形式解があるから、制御変量として使える。

ルックバックオプションは期間中の最大値/最小値で決済する。

C_{\text{lookback}} = e^{-r\tau} \cdot \mathbb{E}^{\mathbb{Q}}\!\left[S_\tau - \min_{0 \leq t \leq \tau} S_t\right]

連続モニタリングならGoldman, Sosin & Gatto (1979) の閉形式解があるけど、離散モニタリングではMCが必要。

バスケットオプションは複数原資産の加重平均。

C_{\text{basket}} = e^{-r\tau} \cdot \mathbb{E}^{\mathbb{Q}}\!\left[\max\!\left(\sum_{j=1}^{m} w_j S_\tau^{(j)} - K,\; 0\right)\right]

多次元MCが必要。相関行列のコレスキー分解で相関する乱数を生成する。

7. ラフ・ボラティリティ

Rough Bergomiモデル (Bayer, Friz & Gatheral, 2016)

\frac{dS_t}{S_t} = \sqrt{V_t}\, dW_t^1

V_t = \xi_0(t) \cdot \mathcal{E}\!\left(\eta\, \sqrt{2H}\int_0^t (t - s)^{H - 1/2}\, dW_s^2\right)

$\mathcal{E}$ はWick指数関数、 $H \approx 0.1$ はハースト指数、 $\eta$ はvol-of-vol、 $\rho = \text{Corr}(dW^1, dW^2)$ 。

なぜ「ラフ」かというと、 $H < 1/2$ のときボラティリティのパスはブラウン運動（ $H = 1/2$ ）よりもラフ（不規則）になる。ヘルダー正則性が $H$ で、 $H \approx 0.1$ はほぼ至るところ微分不可能なパスを意味する。

ATMスキューの短期爆発

ラフモデルの最も劇的な帰結。

\left.\frac{\partial \sigma_{\text{imp}}}{\partial k}\right|_{k=0} \sim C \cdot \tau^{H - 1/2} \quad \text{as } \tau \to 0

$H = 0.1$ だと $\sim \tau^{-0.4}$ 。Hestonでは $H = 1/2$ でスキューは定数。この差は短期オプション（1週間以内）で顕著に表れて、実データとの整合性が劇的に改善される。これはけっこう面白い。

ハースト指数の推定

対数ボラティリティの増分の2次モーメントのスケーリングを使う。

\mathbb{E}\!\left[|\ln V_{t+\Delta} - \ln V_t|^2\right] \propto \Delta^{2H}

両辺の対数をとって $\ln \Delta$ に対して回帰すると、傾きが $2H$ を与える。リアライズド・ボラティリティのプロキシで推定する。

分数ブラウン運動（fBM）のシミュレーションにはCholesky法が $O(n^3)$ と高コスト。Bennedsen, Lunde & Pakkanen (2017) のハイブリッドスキームが実用的で、近い時点はCholesky、遠い時点はリーマン和近似を使う。

8. モデルキャリブレーション

市場で観測されたIV $\{\sigma_{\text{mkt}}(K_i, T_j)\}$ にモデルをフィットする。

\hat{\theta} = \arg\min_\theta \sum_{i,j} w_{ij}\left[\sigma_{\text{model}}(K_i, T_j;\, \theta) - \sigma_{\text{mkt}}(K_i, T_j)\right]^2

重み $w_{ij}$ は通常Vega加重（流動性の高いATM付近を重視）。

キャリブレーション手法をざっくり言うと、Levenberg-Marquardtは計算コスト低いけど局所最適で初期値依存。Differential Evolutionは計算コスト高いけど大域最適で勾配不要、堅牢。NN特性関数近似は訓練は高いけど推論は低くて、1msでサーフェス計算できる（Romer, 2022）。大域最適化（DE）で初期値を見つけてLMで精密化するのが標準。FFTプライシングを使えばキャリブレーション全体が数秒で終わる。

Pythonで書くとこうなる

Heston特性関数（Little Trap版）

import numpy as np

def heston_cf(u, S0, r, tau, V0, kappa, theta, xi, rho):
    u = np.asarray(u, dtype=complex)
    i = 1j

    d = np.sqrt((rho * xi * i * u - kappa)**2 + xi**2 * (i * u + u**2))
    g = (kappa - rho * xi * i * u + d) / (kappa - rho * xi * i * u - d)

    C = (kappa * theta / xi**2) * (
        (kappa - rho * xi * i * u + d) * tau
        - 2 * np.log((1 - g * np.exp(d * tau)) / (1 - g))
    )
    D = ((kappa - rho * xi * i * u + d) / xi**2) * (
        (1 - np.exp(d * tau)) / (1 - g * np.exp(d * tau))
    )

    return np.exp(C + D * V0 + i * u * (np.log(S0) + r * tau))

Carr-Madan FFT

from scipy.fft import fft

def carr_madan_fft(S0, K_target, r, tau, V0, kappa, theta, xi, rho,
                   N=4096, alpha=1.5, eta=0.25):
    lam = 2 * np.pi / (N * eta)
    b = N * lam / 2

    u = np.arange(N) * eta
    k = -b + np.arange(N) * lam

    w = eta * (3 + (-1)**np.arange(N)) / 3
    w[0] = eta / 3

    cf_vals = heston_cf(u - (alpha + 1) * 1j, S0, r, tau, V0,
                        kappa, theta, xi, rho)
    denom = (alpha + 1j * u) * (alpha + 1 + 1j * u)
    psi = np.exp(-r * tau) * cf_vals / denom

    x = np.exp(1j * b * u) * psi * w
    fft_result = np.real(fft(x))

    call_prices = np.exp(-alpha * k) / np.pi * fft_result
    strikes = np.exp(k)

    return np.interp(K_target, strikes, call_prices)

Hestonモンテカルロ（Full Truncation）

def heston_mc(S0, K, r, tau, V0, kappa, theta, xi, rho,
              n_paths=100_000, n_steps=252):
    dt = tau / n_steps
    sqrt_dt = np.sqrt(dt)

    S = np.full(n_paths, S0)
    V = np.full(n_paths, V0)

    for _ in range(n_steps):
        Z1 = np.random.standard_normal(n_paths)
        Z2 = rho * Z1 + np.sqrt(1 - rho**2) * np.random.standard_normal(n_paths)

        V_pos = np.maximum(V, 0)
        S = S * np.exp((r - 0.5 * V_pos) * dt + np.sqrt(V_pos) * sqrt_dt * Z1)
        V = V + kappa * (theta - V_pos) * dt + xi * np.sqrt(V_pos) * sqrt_dt * Z2

    payoff = np.maximum(S - K, 0)
    return np.exp(-r * tau) * payoff.mean()

参考文献

Heston, S. (1993). "A Closed-Form Solution for Options with Stochastic Volatility." Review of Financial Studies, 6, 327-343
Albrecher, H. et al. (2006). "The Little Heston Trap." Working paper
Andersen, L. (2008). "Simple and Efficient Simulation of the Heston Stochastic Volatility Model." J. Computational Finance, 11(3), 1-42
Lord, R., Koekkoek, R. & van Dijk, D. (2010). "A Comparison of Biased Simulation Schemes for Stochastic Volatility Models." Quantitative Finance, 10(2), 177-194
Carr, P. & Madan, D. (1999). "Option Valuation using the Fast Fourier Transform." J. Computational Finance, 2(4), 61-73
Fang, F. & Oosterlee, C.W. (2008). "A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions." SIAM J. Scientific Computing, 31(2), 826-848
Lee, R. (2004). "Option Pricing by Transform Methods." Encyclopedia of Quantitative Finance
Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer
Broadie, M., Glasserman, P. & Kou, S. (1997). "A Continuity Correction for Discrete Barrier Options." Mathematical Finance, 7(4), 325-349
Goldman, M.B., Sosin, H. & Gatto, M.A. (1979). "Path Dependent Options: Buy at the Low, Sell at the High." J. Finance, 34(5), 1111-1127
Kemna, A. & Vorst, A. (1990). "A Pricing Method for Options Based on Average Asset Values." J. Banking & Finance, 14(1), 113-129
Gatheral, J., Jaisson, T. & Rosenbaum, M. (2018). "Volatility is Rough." Quantitative Finance, 18(6), 933-949
Bayer, C., Friz, P. & Gatheral, J. (2016). "Pricing under rough volatility." Quantitative Finance, 16(6), 887-904
Bennedsen, M., Lunde, A. & Pakkanen, M. (2017). "Hybrid scheme for Brownian semistationary processes." Finance & Stochastics, 21(4), 931-965
Romer, S.E. (2022). "Empirical analysis of rough and classical stochastic volatility models." Quantitative Finance
Bayer, Friz, Fukasawa, Gatheral, Jacquier & Rosenbaum (2023). Rough Volatility. SIAM