２. "Implementing QuantLib"の和訳
Chapter VIII The Finite-Difference Framework （有限差分法のフレームワーク）
8.1 古いフレームワーク
Chapter７で既に説明した通り、当初の基本計画は、Treeを使ったモデルと、有限差分法を使ったモデルを、共通のフレームワーク（訳注：Latticeのオブジェクトモデル）の下で作り上げようとしました。しかし、それが実を結ぶことはありませんでした。むしろ有限差分法のプログラムコードは、（Latticeの概念を中心とするのではなく）ベースとなる、差分演算子、境界条件などの数学の概念を中心としたものに近づいて行きました。以下のセクションで、それらについて分析していきます。
8.1.1 Differential Operator (微分演算子)
もし一番最初から説明するなら、まず Arrayクラス（配列クラス）の説明から始めるべきかもしれません。ただ、このクラスの動作は、皆さんの想像通りのものであり、詳細の説明は Appendix Aで行うのでご容赦下さい。ここでは、この配列クラスは、（状態変数の）グリッド \( \{ x_0, x_1, … x_{n-1}\}\) 上の関数 \(f(x)\) を \(f_i=f(x_i)\) のように離散化する為に使われているという事だけ述べておきます。
まず Operator （訳注：ここでは微分演算子あるいは有限差分演算子の事。C++で用意されている ＝, + のような算術演算子の事ではない） の説明から始めます。ここで私が数学的概念の説明をすれば、聖歌隊の人にキリスト教の説教をしているようなものなので、それは最小限に留めます。（もし、有限差分法の数学的な説明が必要であれば、それに関するいくつかの専門書を参考にして下さい。例えば下記Duffy[1]参照） 手短に言えば、Differential Operator（微分演算子）は、ある関数 \(f(x)\) に作用して、それを何らかの微分、例えば一階であれば \(f' (x)\) に変換するものです。ここでは離散化された関数 \(f\) を \(f'\) に変換する作用素を考えます。離散的な微分を行う操作（即ち有限差分）は線形作用素なので、行列で表現できます。基本的には、この Operator（有限差分演算子）は微分の厳密な離散化を与えるものではなく、単に近似を行っているだけです。すなわち、\(f'_i = f'(x_i) + ε_i\) として、右辺の誤差項は、グリッドの間隔を狭くすれば、小さくすることが可能です。
QuantLibライブラリーでは、１階の偏微分\(\partial /\partial x\) と２階の偏微分 \(\partial^2 / \partial x^2\) のDifferential Operator（有限差分演算子）を用意しています。それぞれ、\(D_0、および D_+ D_-\) と呼ばれており、以下のように定義されています。（訳注：式から判る通り、いずれも中心差分）
\[ \frac{\partial f}{\partial x} (x_i ) \approx D_0 f_i=\frac{f_{i+1}-f_{i-1}}{2h},\ \ \ \ \ \ \ \frac{\partial^2 f}{\partial x^2 } (x_i) \approx D_+ D_- f_i=\frac{f_{i+1}-2f_i+f_{i-1}}{h^2} \] 但し \( h=x_i-x_{i-1}\) （ここでは、グリッド間隔が均等と仮定）。
テーラー展開すると、いずれの差分演算子も誤差項の次数は\( O(h^2)\) となることが分かります。
ご覧の通り、インデックスiにおける関数 \(f_i\) に差分演算子を作用させた値は、同じインデックス i、およびその前後のインデックス i+1、i-1 の関数値にのみ依存しています（最初のインデックスと最後のインデックスの取扱いがどうなるかは、今の所無視します）。そうすると、差分演算子\(D_0\) や、他の仲間の演算子は、次のようなTridiagonal Operator（三重対角演算子）で表現できます。
\[ \left( \begin{array}{ccccccc} d_0 & u_0 & 0 & 0 & 0 & \ldots & 0 \\ l_0 & d_1 & u_1 & 0 & 0 & \ldots & 0 \\ 0 & l_1 & d_2 & u_2 & 0 & \ldots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \ldots & 0 & l_{n-4} & d_{n-3} & u_{n-3} & 0\\ 0 & \ldots & 0 & 0 & l_{n-3} & d_{n-2} & u_{n-2} \\ 0 & \ldots & 0 & 0 & 0 & l_{n-2} & d_{n-1} \\ \end{array} \right) \]この行列は、いくつかの望ましい性質を持っています。例えば、三重対角演算子の線形結合も （Black-Scholes-Merton方程式の有限差分スキームにも登場します）、三重対角演算子になります。従って、三重対角行列を表現するため、一般的な Matrixクラスを使わずに、あえて TridiagonalOperatorクラスを作りました。その実装内容を下記 Listing 8.1 に示します。
Listing 8.1：TridiagonalOperatorクラスの実装内容
class TridiagonalOperator {
Size n_;
Array diagonal_, lowerDiagonal_, upperDiagonal_;
public:
typedef Array array_type;
explicit TridiagonalOperator(Size size = 0);
TridiagonalOperator(const Array& low,
const Array& mid,
const Array& high);
static Disposable<TridiagonalOperator> identity(Size size);
Size size() const;
const Array& lowerDiagonal() const;
const Array& diagonal() const;
const Array& upperDiagonal() const;
void setFirstRow(Real, Real);
void setMidRow(Size, Real, Real, Real);
void setMidRows(Real, Real, Real);
void setLastRow(Real, Real);
Disposable<Array> applyTo(const Array& v) const;
Disposable<Array> solveFor(const Array& rhs) const;
void solveFor(const Array& rhs, Array& result) const;
private:
friend Disposable<TridiagonalOperator>
operator+(const TridiagonalOperator&, const TridiagonalOperator&);
friend Disposable<TridiagonalOperator>
operator*(Real, const TridiagonalOperator&);
// other operators
};このクラスの Private のメンバー変数を見ればわかる通り、この Operator の表現方法は、三重対角の構造を反映しています。すなわち、あえて行列として取り扱うのではなく （そうすると、行列の大半の要素がNullになってしまうので）、三重対角の要素を、対角、下対角、上対角の３つの配列として保持し、それぞれが先に示した行列の図の \(d_i,l_i,u_i\) に相当します。
コンストラクターおよび、最初のいくつかのメソッド群は、Operatorインスタンスを生成し、その内容をチェックする為のものです。最初のコンストラクターは、三重対角行列のサイズ（実際には３つの配列の要素数）を設定するだけで、後から setter（変数を設定するメソッド）を使って値を代入し、それをそのまま信用して使います。上記の行列の図をもう一度見て頂くと、setFirstRow( )メソッドは \(d_0 , u_0\) の値を設定し、setLastRow( )メソッドは \(l_{n-2}、d_{n-1}\) の値を設定します。setMidRow( )は各 i に対応する \(l_{i-1}、d_i、u_i\) の値を設定し、さらにそれをひとつにまとめたsetMidRows( )がすべての i について同じ作業を行います。２番目のコンストラクターは、３つの配列の値の設定を一度に全部やってしまいます。そしてファクトリーメソッドである identity( ) は、3番目のコンストラクターのような働きをしており、対角成分（の配列の要素）をすべて 1 に設定し、それ以外の要素をすべて0に設定します （訳注：対角成分がすべて１の単位行列を生成している）。
（注：identity( )メソッドが返す、TridiagonalOperatorへの参照は、Disposableクラステンプレートの中に取り込まれています。これは不要なコピーを避ける為に取られた方法ですが、間違ったやりかたでした。詳細はAppendix Aを参照して下さい）。
各インスペクター関数の役割は明確で、それぞれ対応するメンバー変数の値を返します。
applyTo( ) と solveFor( )メソッドが、インターフェースの中で中心的な役割を負っています。TridiagonalOperator（三重対角演算子オブジェクト）Lと配列uが与えられると、L.applyTo(u)は u に演算子 L を作用させた結果となる L・u を返します。また L.solveFor(u)は u = L･ν を満たす配列 ν を返します。（注：この関数をオーバーロードしたL.solveFor(u,ν)も同じ動作をしますが、メモリー確保の動作を少し節約できます。） いずれの計算についても、後の説明で使います。三重対角演算子を使う事によって、いずれのメソッドも時間軸に対する誤差項のオーダーが \(O(N)\) になるという好都合な性質を持っています。すなわち、計算時間は演算子が取り扱う（グリッドの）サイズ N に比例するだけという事です （訳注：\(N^2\ や\ N^3\) だと計算時間はグリッド数を細かくするとそのべき乗で増えていく）。 計算アルゴリズムは、やや古いですが Numerical Recipes in C（Press et al, 1992、下記[2]参照）から取ってきました。もちろんプログラムコードそのものではありません。その理由は、著作権上の問題だけでなく、プログラムコードが古き良き時代の ｃ言語の慣行を使っており （例えば配列の番号が１から始まるとか）、QuantLibライブラリーには適合しない部分があるからです。
最後に、このクラスは、いくつかの operator （訳注：ここで言う operator演算子は微分演算子ではなく、C++のコードで使われる + や * といった算術演算子の事） をオーバーロードしています。これらの算術演算子を使えば、\(2*L_1+L_2\) のように、三重対角演算子を線形結合する演算に使えます。同様な演算ができるように、QuantLibライブラリーでは、\(D_0\ や\ D_+ D_-\) といった差分演算子（既に説明したものです）も用意しています。例えば、Black-Scholes偏微分方程式が、下記のような演算子の形で表現されているとします。
\[ \frac{\partial f}{\partial t}=\left[ -(r-q-\frac{σ^2}{2}) \frac{\partial}{\partial x}- \frac{σ}{2} \frac {\partial ^2}{\partial x^2 }+rI \right] f \equiv L_{BS} f \]この \(L_{BS}\) 演算子を、（\(D_0 や D_+ D_-\) といった差分演算子を使って） 次のように定義できるでしょう。
\[ TridiagonalOperator\ {\bf L_{BS}} = -\left( r-q-\frac{sigma*sigma}{2} \right) * D_0(N) ~~~~~~~~~~~~~\\ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-\left( \frac{sigma*sigma}{2} \right)* D_+D_-(N) \\ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-r*TridiagonalOperator::identity(N); \]そう言っておきながら、QuantLib自体は、この様に演算子を使っている訳ではありません。QuantLibライブラリーでは、確かに Black-Scholes演算子を定義していますが、上記のような他の演算子を代数的に線形結合させた形にはしていません。
最後の注意点として、（上記 Listingで定義されたようなオブジェクトモデルの）三重対角演算子はおそらく最も計算速度に優れていると思われますが、それだけが有限差分法の中で使われる Operator ではありません。にもかかわらず、それらの演算子については、ベースクラスが存在しません。（訳注：仮にベースクラスを作ったとして、そこから派生する）クラスで実装すべきapplyTo( )やsolveFor( )といったインターフェースは、ここでは単に implicitlyにしか定義されておらず、次の Sub-sectionで説明するテンプレートクラスの中で使われるだけです。結果論から言えば、仮にベースクラスを定義していたとしても、パフォーマンスにそれ程影響を与える事はなかったでしょう。仮想関数の呼び出しとして、applyTo( )やsolveFor( )メソッドを使ったとしても、実際に全体でかかる計算時間からすれば、計算速度への影響は無視できる程小さかったでしょう（訳注：ベースクラスで仮想関数を定義し、内容を派生クラスで実装すれば、実行時Polymorphismが実現できます。しかし、ベースクラスのメソッドを呼び出す事により、若干パフォーマンスを落とします。ところが実際には、メソッドの中のアルゴリズムの計算に相当時間がかかるので、ベースクラスのメソッド呼び出しにかかる時間は、相対的に小さくて済みむので、あえて仮想関数を避ける必要もなかったという意味）。 しかし、開発当時は、テンプレートを使った数値計算技術について必死に学んでいた最中で、Wondermark ４コマ漫画における二人の男のように（注：次のWebsiteを参照して下さい http://wondermark.com/779/）、（階層を使わずに）最後まで突き進む事を選択しました。
[1] D. J. Duffy, Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach. John Wiley and Sons, 2006.
[2] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd edition. Cambridge University Press, 1992.
<ライセンス表示>
QuantLibのソースコードを使う場合は、ライセンス表示とDisclaimerの表示が義務付けられているので、添付します。 ライセンス