← 目次 ｜ ← 前: 16 ｜ 次: 18 →

17. モンテカルロ法

テーマ: ヨーロピアン・コールの価格をモンテカルロ法とBlack-Scholes式で比較

学習点: NumPy のベクトル化, 乱数, 分散減少法(対称変量), 収束の観察

依存: NumPy, SciPy ／ 難易度: 中級

実行方法

uv run 17_option_montecarlo.py

スクリプト冒頭の PEP 723 メタデータ(# /// script)により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

「1 年後に株を 105 で買う権利（ヨーロピアン・コールオプション）」の適正価格を 2 通りの方法で求め、突き合わせるプログラムです。1 つは Black-Scholes の閉形式の公式、もう 1 つは満期の株価を乱数で大量に生成してペイオフの平均を取るモンテカルロ法です。

モンテカルロ法は公式が存在しない複雑な商品にも使える汎用的な手法であり、その実例として最後に「期間中の平均株価」で決まるアジアン・オプション（経路依存型）も同じ枠組みで価格付けしています。試行回数と誤差の関係、分散を減らす工夫（対称変量法）まで含めて、金融工学の実務で使われる数値計算の基本を一通り体験できます。

コードの読みどころ

• mc_call() は Python のループを一切使わず、rng.standard_normal(n) で n 個の乱数をまとめて生成し、ST = S * np.exp(...) → np.maximum(ST - K, 0.0) と配列演算だけで 100 万本の株価とペイオフを計算します。NumPy のベクトル化の威力が最もわかる箇所です。
• 対称変量法は z = np.concatenate([z, -z]) の 1 行です。乱数 z と符号を反転した -z を対で使うと誤差が打ち消し合い、同じ試行回数でも分散が減ります。
• 戻り値の 2 つ目 disc.std(ddof=1) / np.sqrt(n) が標準誤差で、推定値の精度を推定値自身と一緒に報告する作法を示しています。
• グリークスの計算では、seed=3 を揃えた 2 回の mc_call() の差分で Delta を数値微分しています。同一乱数列（common random numbers）を使うことで差分の分散を抑えるテクニックです。
• アジアン型では np.cumsum(..., axis=1) で 20 万本 × 52 週の対数価格経路を一括生成し、paths.mean(axis=1) で経路ごとの平均価格を求めています。

理論的背景

リスク中立測度の下では満期株価が S_T = S_0 exp((r-σ²/2)T + σ√T·Z)（Z は標準正規）に従い、オプション価格は割引期待値 e^{-rT}·E[max(S_T-K,0)] になります。Black-Scholes(1973) はこの期待値を閉形式 C = S_0 N(d1) - K e^{-rT} N(d2) で与えます。一方モンテカルロ推定の標準誤差は O(1/√n) で減少するため、精度を 1 桁上げるには試行回数を 100 倍にする必要があります。

実行結果の見方

理論価格 8.698992 に対し、MC 価格の誤差が試行回数とともに縮み、標準誤差が n を 100 倍にするとほぼ 1/10（0.157→0.016）になる点を確認してください。誤差の列はおおむね標準誤差の 1〜2 倍以内に収まっており、推定が統計的に妥当であることがわかります。

Delta は MC 0.5032 に対し解析解 0.5039 とよく一致します。またボラティリティの表からはベガが正（σ が上がるほど価格が上がる）であること、最後の比較からは同じ経路でもアジアン型（4.20）はヨーロピアン型（8.73）より大幅に安く、平均を取ることで実効ボラティリティが下がる効果が読み取れます。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "numpy",
#     "scipy",
# ]
# ///
"""17: モンテカルロ法 -----------------------------------------------------
テーマ: ヨーロピアン・コールの価格をモンテカルロ法とBlack-Scholes式で比較
学習点: NumPy のベクトル化, 乱数, 分散減少法(対称変量), 収束の観察
根拠: リスク中立測度の下で S_T = S_0 exp((r-σ²/2)T + σ√T·Z), Z~N(0,1)。
      価格は e^{-rT}·E[max(S_T-K,0)]。閉形式は Black-Scholes(1973):
      C = S_0 N(d1) - K e^{-rT} N(d2),
      d1 = (ln(S/K)+(r+σ²/2)T)/(σ√T), d2 = d1 - σ√T。
      モンテカルロ誤差は O(1/√n) で減少する。
"""
import numpy as np
from scipy import stats


def bs_call(S, K, r, sigma, T):
    d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)


def mc_call(S, K, r, sigma, T, n, seed=0, antithetic=False):
    rng = np.random.default_rng(seed)
    if antithetic:
        z = rng.standard_normal(n // 2)
        z = np.concatenate([z, -z])          # 対称変量法（分散減少）
    else:
        z = rng.standard_normal(n)
    ST = S * np.exp((r - sigma ** 2 / 2) * T + sigma * np.sqrt(T) * z)
    payoff = np.maximum(ST - K, 0.0)
    disc = np.exp(-r * T) * payoff
    return disc.mean(), disc.std(ddof=1) / np.sqrt(n)   # 価格と標準誤差


def main() -> None:
    S, K, r, sigma, T = 100.0, 105.0, 0.02, 0.25, 1.0
    exact = bs_call(S, K, r, sigma, T)
    print(f"S0={S} K={K} r={r:.0%} σ={sigma:.0%} T={T}年")
    print(f"Black-Scholes 理論価格 = {exact:.6f}\n")

    print(f"{'試行回数':>10}{'MC価格':>12}{'標準誤差':>11}{'誤差':>10}"
          f"{'対称変量':>12}{'SE(対称)':>11}")
    print("-" * 68)
    for n in [1_000, 10_000, 100_000, 1_000_000]:
        p, se = mc_call(S, K, r, sigma, T, n)
        pa, sea = mc_call(S, K, r, sigma, T, n, antithetic=True)
        print(f"{n:>10,}{p:>12.5f}{se:>11.5f}{p - exact:>+10.5f}"
              f"{pa:>12.5f}{sea:>11.5f}")
    print("  ※ 標準誤差は n を100倍にすると1/10になる（O(1/√n)）。")
    print("    対称変量法は同じ n でも分散を減らせる。\n")

    print("[グリークス（数値微分）と解析解の比較]")
    h = 0.01
    delta_mc = (mc_call(S + h, K, r, sigma, T, 2_000_000, seed=3)[0]
                - mc_call(S - h, K, r, sigma, T, 2_000_000, seed=3)[0]) / (2 * h)
    d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
    print(f"  Delta: MC {delta_mc:.4f} / 解析解 N(d1) = {stats.norm.cdf(d1):.4f}")
    print("  ※ 同一乱数列を使う(common random numbers)ことで差分の分散を抑えている。")

    print("\n[ボラティリティと価格の関係（ベガの符号）]")
    for s in [0.10, 0.20, 0.30, 0.40, 0.60]:
        print(f"  σ={s:>5.0%}: C = {bs_call(S, K, r, s, T):>8.4f}")
    print("  ※ オプション価格はボラティリティの単調増加関数（ベガ>0）。")

    print("\n[経路依存オプション: アジアン型（平均価格）との比較]")
    rng = np.random.default_rng(11)
    n_path, n_step = 200_000, 52
    dt = T / n_step
    z = rng.standard_normal((n_path, n_step))
    logS = (np.log(S) + np.cumsum((r - sigma**2 / 2) * dt
                                  + sigma * np.sqrt(dt) * z, axis=1))
    paths = np.exp(logS)
    asian = np.exp(-r * T) * np.maximum(paths.mean(axis=1) - K, 0).mean()
    euro = np.exp(-r * T) * np.maximum(paths[:, -1] - K, 0).mean()
    print(f"  ヨーロピアン(同じ経路) = {euro:.4f}")
    print(f"  アジアン(算術平均)     = {asian:.4f}")
    print("  ※ 平均を取ることで実効ボラティリティが下がり、価格は低くなる。")


if __name__ == "__main__":
    main()

実行結果

S0=100.0 K=105.0 r=2% σ=25% T=1.0年
Black-Scholes 理論価格 = 8.698992

      試行回数        MC価格       標準誤差        誤差        対称変量     SE(対称)
--------------------------------------------------------------------
     1,000     7.80831    0.46704  -0.89068     8.86830    0.53108
    10,000     8.77383    0.15791  +0.07484     8.68921    0.15748
   100,000     8.68577    0.05058  -0.01322     8.72349    0.05075
 1,000,000     8.72006    0.01603  +0.02107     8.71637    0.01603
  ※ 標準誤差は n を100倍にすると1/10になる（O(1/√n)）。
    対称変量法は同じ n でも分散を減らせる。

[グリークス（数値微分）と解析解の比較]
  Delta: MC 0.5032 / 解析解 N(d1) = 0.5039
  ※ 同一乱数列を使う(common random numbers)ことで差分の分散を抑えている。

[ボラティリティと価格の関係（ベガの符号）]
  σ=  10%: C =   2.7519
  σ=  20%: C =   6.7048
  σ=  30%: C =  10.6925
  σ=  40%: C =  14.6641
  σ=  60%: C =  22.4930
  ※ オプション価格はボラティリティの単調増加関数（ベガ>0）。

[経路依存オプション: アジアン型（平均価格）との比較]
  ヨーロピアン(同じ経路) = 8.7347
  アジアン(算術平均)     = 4.2026
  ※ 平均を取ることで実効ボラティリティが下がり、価格は低くなる。

---

← 目次 ｜ ← 前: 16 ｜ 次: 18 →