Jacques du Toit, Numerical Algorithms Group
Johannes Lotz and Uwe Naumann, LuFG Informatik 12, RWTH Aachen University, Germany
（2015年11月掲載）
概要
10 ファクターのローカルボラティリティモデルでの、10 為替レートのバスケットコールオプションの価格推定にモンテカルロ法を用いた、GPU 利用プログラムについて検討しました。我々は自動微分を用いて、このプログラムのアジョイントバージョンを作成しました。このコードは、360 のオイラー時間差分の 1 万サンプルパスのテストにおいて、438 の入力パラメータ（これらの大半は市場観測予測ボラティリティです）に関する価格のグラジェントを計算するのに、実行時間は 522ms でした（同等の計算を tangent-linear バージョンのコードを用いて CPU で計算すると 2 時間かかります）。
1 イントロダクション
金融においてリスク設定は多かれ少なかれ感度計算と同義です。数理ファイナンスの基礎的な結果から、リスク因子に関する請求の数学的微分に基づいたダイナミックトレーディング戦略により、任意の条件付請求権を如何にヘッジするかが示されています。これらの数学的微分（感度）は、リスク要因に関して最も敏感なエクスポージャーがどこかを教えてくれるだけでなく、それらエクスポージャーを（理想的には）完全に消去する方法も教えてくれます。
感度を計算する最も一般的な方法は有限差分法です。この簡便な "bump and revalue" アプローチを使えば、いかなる数学的微分も推定できます。
しかしながら問題もあります。一つは、これらの推定値（特に高次微分）は基礎となる数値面が滑らかでない場合に劣化しがちで、二つ目には、この方法がコストが掛かる手法であることです。
自動微分は多くのアプリケーションにおいて、有限差分に代わる魅力的な代替手法になりつつあります。自動微分の（以降は簡単に AD と略します）初歩は文献 [1] に紹介がありますが、アイデアは極めて単純です。コンピュータは演算レベルにおいて、二つの浮動小数点に対し、加算、減算、および乗算、除算を行います。これらの基本演算の微分計算は簡単です。数学的な数式やモデルを実装したコンピュータプログラムは強く結びついたこれら多くの基本演算で構成されます。よって、これら演算の微分と共にチェーンルールを用いて、幾つかの入力変数に関するコンピュータプログラムの数学的微分を計算することが出来ます。さらに、テンプレートと演算子オーバーロードといったプログラミング言語の特性を考慮すれば、AD を効率的かつ外部から自動的に実行するソフトウェアツールが、如何に記述されているかが理解できるでしょう。
AD の主要な利点は、コンピュータプログラムの数学的に正確な微分計算です。これをアジョイント法（セクション 4 に記述します）を組合せると、計算効率を大きく増加させることが出来ます。これらの詳細については文献 [1] を参照してください。
AD は与えられたコンピュータプログラムの数学的に正確な微分を計算しますが、次の２つのことに注意が必要です。
- 与えられるコンピュータプログラムは連続微分可能である。
- コンピュータプログラムの微分は、その数式あるいはアルゴリズムの微分に類似したものである。
これら２つの条件は常に満足されるわけではありません。実際に金利のような変数の場合、観測に最小単位（basis point）があるため、実務者はこれらの変数に関する感度を計算するときに、有限差分を好んで用います。どの AD が適切な手法なのかを決める際に、注意深い検討が必要です。
このレポートの目的は、GPU（CUDA）実装されたアプリケーションに対して、AD、特にアジョイント AD を適用できることを示すことです。このプログラムは、複雑な初期ステージ、計算負荷の大きい大規模並列ステージ、およびポストプロセスの３つのステージで構成されます。マニュアル作成のアジョイント GPU カーネルと AD ツール（dco）を組み合わせることにより、全プログラムは 400 以上の入力変数による微分計算が可能になりました。
2 FXバスケット・オプションのローカル・ボラティリティ
多因子ローカル・ボラティリティモデルにおける外国為替（FX）コールオプションの問題を考えます。これは文献 [2] で示されているものと同じ問題です。簡単に概要を以下に記します。
\begin{equation*} C=e^{r_d T}\mathbb{E}(B_T-K)^+ \quad. \tag{1} \end{equation*}
ここで、$r_d$ は自国リスクフリー金利、$K>0$ はストライク価格、$B_T$ は以下の式で与えられます。
\begin{equation*} B_T= \sum_{i=1}^{N} w^{(i)}S_T^{(i)} \quad.\tag{2} \end{equation*}
ここで $i=1,…,N$ に対して、$S^{(i)}$ は、i 番目の原資産（FX レート）であり、$w^{(i)}$ は総和が 1 の重みです。各原資産は次の確率微分方程式（SDE）により導かれます。
\begin{equation*} \frac{dS_t^{(i)}}{S_t^{(i)}}=(r_d-r_f^{(i)})dt+\sigma^{(i)}(S_t^{(i)},t)dW_t^{(i)} \quad. \tag{3} \end{equation*}
ここで、$r_f^{(i)}$ は、i 番目の通貨ペアに対する外国リスクフリー金利で、$(W_t)_{t≧0}$ は相関 $ \langle W^{(i)},W^{(j)} \rangle_t=\rho^{(i,j)}t$ の N 次元ブラウン運動、$\rho^{(i,j)}$ は相関係数で、$i,j=1,…,N$ です。式 $(3)$ の関数 $\sigma^{(i)}$ は未知であり、次の Dupire 式に従って市場予測ボラティリティデータから校正されます。
\begin{equation*} \sigma^2(K,T)=\frac{\theta^2+2T\theta\theta_T+2(r_T^d-r_T^f)KT\theta\theta_K } {(1+Kd_+\sqrt{T}\theta_K)^2+K^2T\theta(\theta_{KK}-d_+\sqrt{T}\theta_K^2)} \quad. \tag{4} \end{equation*}
ここで $\theta$ は市場観測予測ボラティリティ・サーフェース関数であり、$d_+$ は以下の通りです。
\begin{equation*} d_+=\frac{ln(S_0^{(i)}/K)+(r-q+\frac{1}{2}\theta^2)T } {\theta \sqrt{T} }v\quad. \tag{5} \end{equation*}
バスケット・コールオプション価格 $(1)$ はモンテカルロ・シミュレーションを用いて計算されます。
従ってこのモデルの入力は、自国および外国リスクフリーレート、FX レートのスポット値、相関行列、オプションの満期とストライク、重み、および各資産 FX レートの市場観測予測ボラティリティ・サーフェースです。本レポートの目的は、これら変数の各々に対する価格の感度の計算です。
3 GPU適用コンピュータプログラム（AD の準備）
3.1 ステージ 1：セットアップ
Dupire 式 $(4)$ は、無限に多くの予測ボラティリティクォートがある事を想定します。しかしながら、マーケットは異なる満期とストライクでの一握りのコールオプションを取引します。よって離散的な少数のクォートから、$(4)$ で扱えるように十分に滑らかな関数 $\theta$ を作成しなくてはなりません。これを無裁定評価手法で扱うことは簡単ではなく、多くの研究トピックにされてきました。マーケット実務者がとる通常のアプローチは、無裁定的内挿の細部の点を無視し、単純に三次スプラインを用いる事です。これが我々の採用する方法です。
プログラムの最初のステージは、入力に予測ボラティリティ・サーフェースを用い、1 次元の三次スプラインを集めることで表現されるローカルボラティリティ・サーフェースを生成することです。この手続きの詳細は文献 [2] に与えられていますが、ここではその複雑さを感じてもらうために概要のみを示します。
FX マーケットは通常満期毎の５つのデルタのオプションクォートを提供します。各満期の５つのクォートを通してストライク方向に三次スプラインをフィッティングします。そして各満期に対し、区分単調エルミート多項式を時間 0 から T まで、各多項式がクォートの一つを通るようにフィッティングします。
各満期のクォートは異なるストライクであることから、すでにフィッティングした三次スプラインによる内挿値を必要とします。二つの満期間のモンテカルロ時間ステップに対して、このエルミート多項式を、５つの各ストライクにおける予測ボラティリティの内挿に用います。次に、モンテカルロ時間ステップ毎に、計算し終えた５つの予測ボラティリティを通してストライク方向に三次スプラインをフィッティングします。これらスプラインにより、上述の式 $(4)$ の $\theta,\theta_K,\theta_{KK}$ が計算可能です。$\theta_T$ の計算には、本質的に同等のこの手順を用いることが出来ます。最終的に、これら微分値は式 $(4)$ に代入され、グリッド上のローカルボラティリティ・ポイントおよび、これらのポイントへフィッティングされた三次スプラインのセットを生成します。これらのスプラインの最終ポイントの傾きは、ローカルボラティリティ・サーフェースの左右（ストライク方向）のテイルを覆う線形外挿関数を形成するように射影されます。
3.2 ステージ 2：GPU で加速したモンテカルロ法
ここで $(3)$ 式を解くためにロジスティック過程のオイラー・丸山スキームを用います。これは全プログラムの中で最も計算負荷の高い部分であり、GPU を利用して計算します。オイラー時間ステップを実行する際には、$\Delta T$ 間隔で $n_T$ ステップ実行します。ここで、$n_T \Delta T＝T$ であり、$T$ はバスケットオプションの満期です。$\tau=1,…,n_T$ に対して $S_\tau^{(i)} \equiv S_{\tau\Delta t}^{(i)}$ として、$(3)$ 式の離散化版は以下のようになります。
\begin{eqnarray} & &log(S_{\tau+1}^{(i)})-log(S_{\tau}^{(i)}) \\ & & =\left(r_d-r_f^{(i)}-\frac{1}{2}\sigma^{(i)}(S_{\tau}^{(i)},\tau\Delta t)^2 \right)\Delta t +\sigma^{(i)}(S_{\tau}^{(i)},\tau\Delta t)\sqrt{\Delta t}\mathbf{d_i}\mathbf{Z_\tau} \quad. \tag{6} \end{eqnarray}
ここで、$\mathbf{Z_\tau}$ は標準正規乱数の $N×1$ の行ベクトル、$\mathbf{d_i}$ は $N×N$ の相関行列のこれスキー分解の i 番列、即ち $\mathbf{d_{i_1}}\mathbf{d_{i_2}}^T=\rho^{(i_1,i_2)}$, $i_1,i_2,…,N$ となります。このステージの出力は、各j番のモンテカルロ・サンプルパスの、各i番の原資産、$i=1,…,N$ に対する、$log(S_{nT}^{(i)})\equiv log(S_{nT}^{(i)}(j)) $（注 a）の最終値です。
注 a：今後は、サンプルパスに焦点を当てたいときは $S_\tau^{(i)}(j)$ と記しますが、記述の簡潔化のために $(j)$ を省略します。
3.3 ステージ 3：ペイオフ
最終ステージでは、モンテカルロ・サンプルパスの末端値を CPU へ返し、式 $(1)$ の価格を計算します。ペイオフ（清算計算）が GPU で実行している間（文献 [2] で行われたものです）、CPU で実行すべきポストプロセスステージであるファイナンス計算を模倣するこれら最終ステップを実行します。
3.4 実装
上述のプログラムは C++ で記述されており、GPU カーネルは CUDA を用いています。実行時間のほとんどはモンテカルロ計算のため、CPU の C++ コードは通常の作成方法で、特別な最適化を行っていません。対照的に、効率的な CUDA コードの作成には注意が必要です。設計上の最も重要な点は、単独のモンテカルロスの時間ステップ（$\tau$）に対する三次スプライン計算をシェアードメモリにロードし、次の時間ステップに対する三次スプラインをロードする前に、全てのサンプルパスを時間 $\tau$から$\tau+1$ へ繰り上げることです。このスプラインはキャッシュに出来る限り長時間保存され、できる限り小さなキャッシュ領域を使うようにします。こうしてこのコードは、小さいだけでなく大きなスプラインにおいても、出来るだけ多くの資産クラスを扱えるようにうまく動作させることが出来ました。
4 アジョイント自動微分
解析的な、あるいは自動微分から得られるアジョイント法は、洗練された数値計算技法です。ここで簡単な紹介をしますが、深く調べたいときは文献 [1] を参照してください。
AD はコンピュータプログラムを微分する技術です。AD には主要な二つのモード（モデル）があります。関数 $f:\mathbb{R}^n\rightarrow\mathbb{R}$ すなわち $y=f(\mathbf{x})$ を実装したコンピュータプログラムを考えましょう。ベクトル $\mathbf{x}^{(1)}\in\mathbb{R}$ を用いて、関数 $F^{(1)}:\mathbb{R}^{2n}\rightarrow\mathbb{R}$ を以下のように定義します。
\begin{equation*} y^{(1)} = F^{(1)}(\mathbf{x},\mathbf{x}^{(1)}) = \nabla f(\mathbf{x}) \cdot \mathbf{x}^{(1)} = \left(\frac{\partial f}{\partial \mathbf{x}} \right) \cdot \mathbf{x}^{(1)} \quad. \tag{7} \end{equation*}
ここで、ドット記号は通常のベクトルの内積を表します。関数 $F^{(1)}$ は $f$ の tangent-linear モデルと呼ばれ、AD の最も単純な形式です。$\mathbf{x}^{(1)}$ を $\mathbb{R}^n$ 全域に渡らせて繰り返し $F^{(1)}$ を呼び出せば、入力変数 $x_i,i=1,…,n$ 毎の $f$ の個々の偏微分 $y^{(1)}$ が得られます。よって、全グラジェント $\nabla f$ を得るためには、tangent-linear モデルをn回呼び出す必要があります。これは、リスク計算に実行時間が、大まかに言って $f$ の計算時間の n 倍になることを意味します。
次に、任意の値 $\mathbf{y}_{(1)}\in\mathbb{R}$ を用いて以下で定義される関数 $F_{(1)}:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^n$ を考えましょう。
\begin{equation*} \mathbf{x}_{(1)} = F_{(1)}(\mathbf{x},y_{(1)}) = y_{(1)} \nabla f(\mathbf{x}) = y_{(1)} \frac{\partial f}{\partial \mathbf{x}} \quad. \tag{8} \end{equation*}
関数 $F_{(1)}$ を関数 $f$ のアジョイントモデルと呼びます。$\mathbf{y}_{(1)}=1$ として、アジョイントモデル $F_{(1)}$ を一度呼び出せば、$f$ の偏微分の全ベクトルを計算できます。さらに、一般に $F_{(1)}$ の計算には、$f$ の計算に掛かる flops（浮動小数点演算数）の 5 倍以上を必要としないことが証明されています。よってアジョイントモデルは、大きなグラジェントを極めて低いコストで得ることのできる強力な方法です。
数学的には、アジョイントモデルは補助的なスカラー関数 $t$ の偏微分として定義でき、上述の式 $(8)$ で以下のように記述できます。
\begin{equation*} y_{(1)} = \frac{\partial t}{\partial y} \\ \mathbf{x}_{(1)} = \frac{\partial t }{\partial \mathbf{x}} \tag{9} \end{equation*}
これはアジョイントモデルにとって微妙かつ重大な意味合いを持ちます。ここで関数 $f:\mathbb{R}\rightarrow\mathbb{R}$ を考え、この関数が $x$ を入力し、単純な演算（加減乗除算）の系列により、$x$ を中間変数を通して、出力 $y$ に至る場合を考えます。
\begin{equation*} x \mapsto \alpha \mapsto \beta \mapsto \gamma \mapsto y \tag{10} \end{equation*}
ここで、$\alpha,\gamma,\beta$ は中間変数です。アジョイントモードの定義 $(9)$ を用いれば、以下のように書くことが出来ます。
\begin{eqnarray} x_{(1)} = \frac{\partial t }{\partial x} &=& \frac{\partial \alpha }{\partial x}\frac{\partial t }{\partial \alpha} \tag{11} \\ &=& \frac{\partial \alpha }{\partial x}\frac{\partial \beta }{\partial \alpha}\frac{\partial t }{\partial \beta} \\ &=& \frac{\partial \alpha }{\partial x}\frac{\partial \beta }{\partial \alpha}\frac{\partial \gamma }{\partial \beta}\frac{\partial t }{\partial \gamma} \\ &=& \frac{\partial \alpha }{\partial x}\frac{\partial \beta }{\partial \alpha}\frac{\partial \gamma }{\partial \beta}\frac{\partial y }{\partial \gamma}\frac{\partial t }{\partial y}=\frac{\partial y }{\partial x}y_{(1)} \end{eqnarray}
これは $f$ のアジョイントモデル $F_{(1)}$ に他なりません。ただしアジョイントモデルの入力は $y_{(1)}$ であることに注意しなくてはなりません。よってアジョイントモデルを計算するには、$y_{(1)}$ から始めて、それを $x_{(1)}$ まで式 $(11)$ を上に辿っていきます。つまり以下を計算します。
\begin{equation*} \left(\left(\left( y_{(1)} \cdot \frac{\partial y }{\partial \gamma} \right) \cdot \frac{\partial \gamma }{\partial \beta} \right) \cdot \frac{\partial \beta }{\partial \alpha} \right) \cdot \frac{\partial \alpha }{\partial x} \tag{12} \end{equation*}
これは $f$ を実装したコンピュータプログラムを「逆向き」に実行することを意味します。$\frac{\partial y }{\partial \gamma}$ を計算するには $\gamma$ を知らなくてはいけない（$\alpha$ や $\beta$ についても同様に）ことから、プログラムを逆向きに実行するためには、最初にフォワード実行（通常の実行）し、偏微分計算に必要な中間値を保存して置く必要があります。一般に、アジョイントコードは必要となる中間値の計算に際し、その全てを保存するための巨大な保存領域が必要とされる場合があります。
5 dcoツール
アジョイント AD を実装するためには２つの選択肢があります。一つはマニュアルでアジョイントコードを記述する、もう一つは、アジョイントコードを自動的に作成するソフトウェアツールを使うことです。アジョイントコードを手動で書くのは煩雑でエラーを引き起こしやすい作業で、コードが大きい場合はコストが掛かりすぎます。そうした場合はツールを使うのが唯一のオプションです。
AD ツールの実装には、アジョイントコードを生成するために２種類の異なるアプローチを取ります。これらアプローチは共に、ソースコードを必要とし、共に最初のステップ（フォワード実行）でデータを収集して、このデータを使って二番目のステップ（逆向き実行）で上述のように実行します（さらに詳細な説明は [1] を参照してください）。
第一のアプローチは、ソースコードコード変換ツールを用いて、プログラムを静的解析してアジョイントコードを生成する方法です。通常は効率的にコードが生成されます。このアプローチの主要な欠点は利用言語です。ソースコード変換ツールは単純なプログラミング言語を多くカバーしていますが、C++ のようなより先進的な言語をサポートしていません。特に C++ に関しては、無視できない特別なコード書き換えが、ツールにプログラムを理解させるために必要になります。
第二のアプローチは、アジョイントモデルの計算に演算子オーバーロードを用いる方法です。このアプローチは対象プログラムが用いる言語を用いることから、全言語をサポートしています。基本的な考え方は、元のプログラムを見渡して全浮動小数データの型を新しい特別な型に置き換えることです。単純な C++ コードでは通常は、これがやること全てですが、現実的なシミュレーションコードではさらに処理が必要でしょう。オーバーロードの技法のマイナス面は、実行時のパフォーマンス低下の可能性です。オーバーロードを用いた全ての AD ツールは、入力から出力までの依存性を追跡して保存する通常実行のためにあるデータ構造を構築します。このデータ構造は通常「テープ」と呼ばれます。逆向きの実行中に、アジョイントモデルを実行するためにテープはプレイバックされます。このツールの性能は通常は、このデータ構造が如何に効率的に作られたかに密接に関係します。
これまでの２つのアプローチは互いに排他的ですが、C++ のような言語にとってはある共通点があります。AD ツール「dco」はこれを利用しようとするものです。本質的には演算子オーバーロードを用いたアプローチですが、コンパイル時間の最適化を最大限に生かすためにテンプレートのような効率的な C++ の特性を用います（例えば文献 [4] を参照してください）。さらに dco は、アジョイントモデルのコードであれば自身に統合するユーザフレンドリーな方法として、外部関数インターフェースを持ちます。このコードはソース変換ツール（第一のアプローチで示したような）や手動で作成したコードのどちらでも可能です。
5.1 性能と特徴
AD ツールの典型的な性能指標は、元のプログラム実行時間に対するアジョイントコードの実行時間（即ち全一次微分）比として以下のように与えられます。
\begin{equation*} R=\frac{アジョイントコード実行時間}{元のプログラムの実行時間} =\frac{time(\nabla f)}{time(f)} \end{equation*}
これは、グラジェント計算に対して AD ツールがもたらすオーバーヘッドを図るのに、効果的な指標です。dco を実際のプログラムに適用したとき、その比 R の典型的な値は 5 から 15 の間です。与えられたコードの R の最大値は、コードの効率（例えばキャッシュ効率など）やアジョイントコード作成に費やす努力に強く依存します（dco はメモリ利用効率や実行時間に影響のある、幾つかの実行時やコンパイル時のオプションを持ちます）。
以下に dco の重要な特徴を示します。
- ジェネリック微分型：dco のデータ型は任意の基本型（例えば double や float）のインスタンスを持ちます。
- 高次微分計算が可能
- 外部関数インターフェース：
・メモリ使用率を制御するためのチェックポイント機能（例えば PDE 解法の時間ステップにおいて）
・ユーザー定義組込み関数（例えば滑らかなペイオフ関数）
・並列化を可能にする低メモリアジョイントモデル（例えばモンテカルロ法）
・本報告のようなアクセラレーターのサポート（例えば GPU）
- NAG ライブラリのサポート：如何なる NAG ライブラリの dco バージョンも可能で、dco は組込み関数として扱わせることが可能
最後に簡単に纏めると、dco は次のステップでアジョイント AD を実行します：
- プログラム中の浮動小数点データ型を dco データ型へ置き換えます。
- プログラムを通常実行し、テープに中間計算結果を最適な状態で保存します。
- $y_{(1)}=1$としてテープをプレイバックし、全計算を効果的に逆向き実行してグラジェントを計算します。
5.2 アジョイント AD と GPU アクセラレーター
GPU プログラミング・プラットフォームのハードウェア的な制約は、通常の CPU プログラミングとかなり異なります。特に比較的に大きな数のスレッド（典型的な CUDA カーネルでは 4 万スレッドまで利用可能）に対して、小さなメモリ（約 6GB）しかありません。また L1, L2 キャッシュはさらに小さなものです。この事は GPU プログラミングにおいて、メモリの有効利用が最も重要であることを示しています。
残念ながらアジョイント AD は、膨大なメモリを必要とします。さらに現在、CUDA や OpenCL をサポートする AD ツールは無く、GPU アジョイントコードは皆手作業で作られたものです（tangent-linear モードの dco をサポートする CUDA は開発中です）。このため、GPU 上で実行可能なアジョイントコードには数個の例しかありません。Mike Giles が作成した LIBOR マーケットモデルのコード（彼の Web サイト、あるいは文献 [3] を参照してください）が、我々の気づいた唯一のパブリックドメイン・コードです。
後半の議論は GPU アプリケーションのアジョイント自動微分（後半部） をご覧ください。