２. "Implementing QuantLib"の和訳
Chapter VIII The Finite-Difference Framework （有限差分法のフレームワーク）
8.1 古いフレームワーク
8.1.2 Evolution schemes: 時間軸方向の差分スキーム
偏微分方程式 \(\frac{\partial f}{\partial t }= L ･ f\) において、前セクションで説明した差分演算子は、右辺の微分演算子の差分化スキームを提供しています。それに対し、Evolution scheme(時間軸方向の差分スキーム) は左辺の時間軸での偏微分を差分化したものです。
ご存知の通り、金融工学の世界で使われている有限差分法では、時間軸を進む方向は、関数値 \(f(T)\) が既知である“満期日”（すなわち満期時における Pay-off関数） からスタートし、時間を現在まで遡って行きます。その間、各Timeステップにおいて順番に、\(f(t)\ を\ f(t+\Delta t)\) から導出します。
上記の偏微分方程式の差分近似を計算する、最も簡単な方法は、下記の式で表されます。
\[ \frac{f(t+\Delta t)-f (t)}{\Delta t}=L ∙ f(t+\Delta t) \]この式をさらに単純にすると
\[ f(t)=f(t+∆t)-\Delta t ∙ L ∙ f(t+\Delta t) ~~~または~~~ f(t)=(I-\Delta t ∙ L) ∙ f(t+\Delta t) \]となります。（訳注：式中の \( f \) は配列、I と L は行列）。 右辺は単純な行列と配列の積であり、行列を、関数値配列に作用させる動作をプログラムコードに直すと次のようになります （訳注： f が a で表現されている点に注意）。
a = (I - dt*L).applyTo(a);（当然ながら、実際のプログラムコードでは、（I－dt*L） の部分は事前に計算され時間ステップごとに保持されています。）
これよりやや複雑なEvolutionスキームは、下記の式で表されます。
\[ \frac{f(t+\Delta t) -f (t)}{\Delta t}=L∙f(t) \]この式をさらに変換すると、
\[ (I+\Delta t∙L)∙f(t)=f(t+\Delta t) \]となります。この式から \(f(t)\) を求めるには、連立１次方程式を解く必要があるので（あるいは演算子行列の逆行列を計算するので）、より複雑な計算ステップになります。しかし、我々の設計した演算子インターフェース（訳注：solveFor( )メソッドのこと）を使えば、以下のように、applyTo( )と同じくらい簡単なプログラムコードで済みます。
a = (I + dt*L).solveFor(a);ここで三重対角演算子を使うメリットが発揮されます。
まず \(f(t)\) を求める２つの方法の違いについて、ごく簡単に説明させて下さい。端的に言えば、前者のスキームは（explicit schemes（陽的解法）と呼ばれています）、計算はより簡単ですが解が不安定になる可能性があります。一方、後者のスキームは（implicit schemes（陰的解法）と呼ばれています）、より安定した解が導出でき、また離散時間の間隔も大きく取れます（訳注：陽的スキームでは解の安定性を確保する為に、時間間隔と状態変数の間隔の取り方に一定の制約が発生する。一方、陰的スキームではそういった制約が小さい）。仮に（行列の形を取る）一般的な Operator（差分演算子）であれば、implicit step（陰的解法のための時間間隔の取り方）を使う計算時間は、グリッドのサイズNに対して、\(Ο（N^3）\) のオーダーになり、explicit step（陽的解法のための時間間隔の取り方）の \(Ο(N^2)\) に対して、実用性で劣ります。しかし、三重対角演算子は、どちらの方法でも（ステップ数に対する計算時間は） \(O(N)\) のオーダーとなり、追加の計算時間の心配をせずに、explicit stepの利点を生かす事ができます。
QuantLibライブラリーでは、上記の２つの解法 （それぞれexplicit Euler scheme（陽的オイラー法）、implicit Euler scheme（陰的オイラー法）と呼ばれています）を用意しています。実際の所、この２つの解法は、より汎用的な MixedSchemeテンプレートクラスから派生する２つの特殊な方法として実装されています。MixedSchemeクラステンプレートのコードを下記 Listing 8.2に示しますが、これは次の Evolutionスキーム式をオブジェクトモデル化したものです。
\[ \frac{f(t+\Delta t)-f(t)}{\Delta t}=L∙ \left[(1-\theta )∙f(t+\Delta t)+\theta ∙f(t)\right] \]この等式は、implicit と explicit step（陰的と陽的ステップ）がミックスされています。この等式かわわかる通り、\(\theta =0\) とすると陽的オイラー法となり、\(\theta =1\) とすると陰的オイラー法となります。θは 0 と 1の 間のどのような値を取る事も可能で、仮に 0.5 とすると Crank-Nicolson（クランク・ニコルソン法）になります。
Listing 8.2：MixedSchemeクラステンプレートと、その派生クラスの概要
template <class Operator>
class MixedScheme {
public:
typedef OperatorTraits<Operator> traits;
typedef typename traits::operator_type operator_type;
typedef typename traits::array_type array_type;
typedef typename traits::bc_set bc_set;
typedef typename traits::condition_type condition_type;
MixedScheme(const operator_type& L,
Real theta,
const bc_set& bcs)
: L_(L), I_(operator_type::identity(L.size())),
dt_(0.0), theta_(theta) , bcs_(bcs) {}
void setStep(Time dt) {
dt_ = dt;
if (theta_!=1.0) // there is an explicit part
explicitPart_ = I_-((1.0-theta_) * dt_)*L_;
if (theta_!=0.0) // there is an implicit part
implicitPart_ = I_+(theta_ * dt_)*L_;
}
void step(array_type& a, Time t) {
if (theta_!=1.0) {
a = explicitPart_.applyTo(a);
if (theta_!=0.0) {
implicitPart_.solveFor(a, a);
}
protected:
operator_type L_, I_, explicitPart_, implicitPart_;
Time dt_;
Real theta_;
bc_set bcs_;
};
template <class Operator>
class ExplicitEuler : public MixedScheme<Operator> {
public:
// typedefs, not shown
ExplicitEuler(const operator_type& L,
const bc_set> >& bcs)
: MixedScheme<Operator>(L, 0.0, bcs) {}
};既に述べた通り、Evolution Schemeテンプレートクラスは、何等かの Differential Operator型を、テンプレート引数として取ります。この MixedSchemeクラスでも同様で、Differential Operatorクラスの中から、OperatorTraitsテンプレートを使って、必要な Operatorを選択します（これについては後で説明します）。
MixedSchemeクラスのコンストラクターは、引数として、①差分演算子Lと、②θの値（訳注： \(0 \leq \theta \leq １\)で、implicit とexplicitスキームに対するウェイトとして使われる）、及び③境界条件（これについても後で説明します）を取り、それぞれメンバー変数に保存します。setStep( )メソッドは、離散時間の間隔となる Δt を保存し、差分演算子となる
\[ I－(1－\theta ）･\Delta t・L ~~~ および ~~~ I +\theta ･\Delta t・L \]を事前に計算しておきます。その演算子を使って、applyTo( ) と solveFor( )メソッドが呼び出され事になります。θの値が 0 または 1 で設定されている場合は、いずれか一方の事前計算は、不要なので行われません。この事前計算を、コンストラクターとは別の setStep( )メソッドで行うことにより、仮にシミュレーションの途中で離散時間の間隔を変える必要がある場合でも、このメソッドを使えば同じ Evolution Schemeのインスタンスが再利用できます。
最後に、step( )メソッドは、関数値の配列を求めるアルゴリズムを、時間軸に沿って遡及する動作を行います。上記 Listing中の実装内容は、基本的な概要を示しているだけで、次のセクションで、さらに詳しく説明します。ここでは、（θの値に従って）陽的解法か陰的解法かのいずれかを実行するようになっています （ここでも θ が 0 か 1 の場合は、不必要な計算を避けるようにしています）。
QuantLibライブラリーの中で提供されている、いくつかの Evolution Schemeは （このListingにあるExplicitEuler法など）、すべて MixedSchemeクラスの特殊な形態です。 （注：ここでは「古いフレームワーク」について説明しています。新しいフレームワークでは別の方法を提供しています。）
Evolution Schemeについても（Differential Operatorと同様）ベースクラスは存在しません。有限差分フレームワークの中の別のオブジェクトが、Evolution Schemeをテンプレート引数として取り込み、それが提供する setStep( )メソッドや、step( )メソッドを使うようになっています。
<ライセンス表示>
QuantLibのソースコードを使う場合は、ライセンス表示とDisclaimerの表示が義務付けられているので、添付します。 ライセンス