WebGPUで学ぶモンテカルロ法① ヨーロピアン・オプション

動作環境・注意事項

Chrome 113+	✔ 推奨。フラグ不要で WebGPU 動作
Edge 113+	✔ 動作します
Firefox	✘ WebGPU 未対応（CPUフォールバックで計算）
Safari 18+	△ 実験的対応。A17 Pro 等では動作
モバイル	△ Android Chrome 対応機種で動作
非対応時	△ WebGPU が無い環境では自動的に CPU 計算に切替

金融の背景：オプションとは

オプションとは、「将来のある時点で、ある価格で株などを売買できる権利」のことです。義務ではなく権利なので、有利なら行使し、不利なら放棄できます。この「権利そのもの」に値段がつきます。

コール・オプション（買う権利）

payoff = max(S_T − K, 0)

満期の株価 S_T が権利行使価格 K より高ければ、その差額が利益。安ければ権利放棄で損失ゼロ（権利料を除く）。

プット・オプション（売る権利）

payoff = max(K − S_T, 0)

満期の株価 S_T が K より低ければ、その差額が利益。株価下落に対する保険のように使えます。

問題は「この権利に今いくら払うべきか」です。満期の株価は誰にもわかりません。そこで「株価が確率的にどう動くか」をモデル化し、満期ペイオフの期待値を現在価値に割り引くことで理論価格を求めます。この期待値計算にモンテカルロ法を使います。

📈 株価のモデル：幾何ブラウン運動
株価 S はランダムに変動しますが、よく使われるモデルが幾何ブラウン運動です。満期の株価は次式で生成されます（Z は標準正規乱数）。

満期株価の生成式（リスク中立測度）

S_T = S₀ × exp[ (r − ½σ²)T + σ√T · Z ]

動作原理と内部の仕組み

モンテカルロ法による価格評価は、次の4ステップで構成されます。各ステップが GPU 上で並列実行されることで、数百万パスが瞬時に計算されます。

正規乱数の生成

各スレッド（=各パス）が、自分のIDから決定的に標準正規乱数 Z を生成します。GPUには Math.random() がないため、PRNG を自前実装します（本記事の技術的な核心）。

満期株価の計算

幾何ブラウン運動の式に Z を代入し、満期時の株価 S_T を1つ計算します。1スレッドが1つの未来シナリオを担当します。

ペイオフの計算

そのシナリオでのペイオフ max(S_T−K, 0) を計算します。各スレッドが独立に1つのペイオフ値を出力します。

平均して割引

全パスのペイオフを平均し、現在価値に割り引きます（e^−rT を掛ける）。これがオプションの推定価格です。

⚙️ なぜGPUが効くのか
各パスは完全に独立しています。パスAの計算はパスBに一切影響しません。この「独立した大量の試行」こそGPUの並列性が最も活きる構造です。 100万パスなら100万スレッドが同時に走ります。

コードのポイント解説

本記事の核心、WGSL内での正規乱数生成を中心に解説します。

カウンタベース PRNG（PCG hash）

// GPUには乱数生成器がないので、スレッドIDから決定的に乱数を作る
// PCG hash: 整数を入力すると擬似乱数の整数を返す
fn pcg(n: u32) -> u32 {
  var h = n * 747796405u + 2891336453u;
  h = ((h >> ((h >> 28u) + 4u)) ^ h) * 277803737u;
  return (h >> 22u) ^ h;
}

// 0〜1 の一様乱数に変換
fn rand(seed: u32) -> f32 {
  return f32(pcg(seed)) / 4294967295.0;
}

各スレッドは「グローバルID + パラメータ」をシードにして pcg() を呼びます。同じシードからは必ず同じ乱数が出る（決定的）ため、再現性が保てます。スレッドごとにシードをずらすことで、各パスが異なる乱数列を持ちます。

Box-Muller 変換（一様乱数 → 正規乱数）

// 一様乱数 2つから、標準正規分布に従う乱数を作る
fn randNormal(seed: u32) -> f32 {
  let u1 = max(rand(seed * 2u), 1e-7);       // log(0)回避
  let u2 = rand(seed * 2u + 1u);

  // Box-Muller: √(-2 ln u1) · cos(2π u2) が標準正規乱数になる
  return sqrt(-2.0 * log(u1)) * cos(6.2831853 * u2);
}

Box-Muller 変換は、一様分布の乱数を正規分布に変換する古典的手法です。 u1 に 0 が入ると log(0) = -∞ になるため、 max(u1, 1e-7) でガードするのが実装上の重要ポイントです。

メインの Compute Shader（1スレッド=1パス）

struct Params {
  s0: f32, k: f32, sigma: f32, r: f32,
  t: f32, isCall: f32, numPaths: u32, seed: u32,
}
@group(0) @binding(0) var<storage, read_write> payoffs: array<f32>;
@group(0) @binding(1) var<uniform> p: Params;

@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) gid: vec3u) {
  let i = gid.x;
  if (i >= p.numPaths) { return; }

  // ① このパス専用のシードで正規乱数を生成
  let z = randNormal(i + p.seed * 2654435761u);

  // ② 幾何ブラウン運動で満期株価を計算
  let drift = (p.r - 0.5 * p.sigma * p.sigma) * p.t;
  let diff  = p.sigma * sqrt(p.t) * z;
  let sT = p.s0 * exp(drift + diff);

  // ③ ペイオフを計算（コール or プット）
  var payoff: f32;
  if (p.isCall > 0.5) {
    payoff = max(sT - p.k, 0.0);
  } else {
    payoff = max(p.k - sT, 0.0);
  }

  // ④ 結果を書き込む（平均と割引はCPU側 or 次の記事のGPUリダクションで）
  payoffs[i] = payoff;
}

各スレッドが「乱数生成 → 株価 → ペイオフ」を独立に実行します。この記事では各パスのペイオフをバッファに書き出し、平均はCPU側で計算します。（平均処理自体をGPUで高速化する「並列リダクション」は次回・速度編のテーマです）

⚠️ ハマりどころ：ワークグループ数の上限
dispatchWorkgroups() で指定できるワークグループ数には上限（多くの環境で 65,535）があります。 workgroup_size=64 の場合、約419万パスを超えると上限に達し、それ以上ではシェーダが実行されず結果が全て0になります（エラーも出ません）。
対策は grid-stride ループです。1スレッドが i += totalThreads で間隔をあけて複数パスを担当することで、ワークグループ数を上限以下に固定したまま任意のパス数を処理できます。 500万パスでの計算でこの実装が効いています。

Black-Scholes 解析解（答え合わせ用）

// ヨーロピアン・オプションには厳密な解析解がある
function blackScholes(S0, K, sigma, r, T, isCall) {
  const d1 = (Math.log(S0/K) + (r + 0.5*sigma*sigma)*T) / (sigma*Math.sqrt(T));
  const d2 = d1 - sigma*Math.sqrt(T);

  if (isCall) {
    return S0*normCDF(d1) - K*Math.exp(-r*T)*normCDF(d2);
  } else {
    return K*Math.exp(-r*T)*normCDF(-d2) - S0*normCDF(-d1);
  }
}

ヨーロピアン・オプションは Black-Scholes 式という厳密解を持ちます。モンテカルロの推定値がこの理論値に収束すれば、実装が正しいことの証明になります。解析解が存在しないエキゾチック・オプション（応用編）では、この答え合わせができません。

収束の考察

デモで「収束を観察」ボタンを押すと、パス数を増やすにつれて推定値が理論値に近づく様子が見えます。モンテカルロ法の誤差はパス数 N の平方根に反比例して減少します。

標準誤差の法則

標準誤差 ∝ 1 / √N

収束の特徴 01

誤差は 1/√N

精度を10倍にするにはパス数を100倍にする必要があります。これがモンテカルロ法の宿命的な弱点であり、GPU並列化が重要になる理由です。

収束の特徴 02

次元に強い

収束速度がパス数だけで決まり、問題の次元に依存しません。多資産オプションなど高次元問題で他手法より有利になります。

収束の特徴 03

並列化で実用化

1/√N の遅い収束を、GPUで大量パスを高速計算することで補います。だからこそモンテカルロ × GPU は好相性なのです。

⚠️ 免責事項
本デモは WebGPU とモンテカルロ法の学習・技術実験を目的としたものです。価格モデルは一定ボラティリティなどの理論上の仮定に基づく簡略化されたものであり、 実際の投資判断・オプション取引の根拠として使用しないでください。 オプション取引はリスクを伴い、損失が生じる可能性があります。投資判断はご自身の責任で行ってください。

NEXT — シリーズ第2回

② 速度編：CPU vs WebGPU — 並列化の威力を測る

→

モンテカルロ法による
ヨーロピアン・オプション価格評価

概要