第98回 有限差分法で PDE を解く ― 陽解法、陰解法、Crank-Nicolson

位置: 第10部 第98回 / 前提: 第66回、第68回 / 数学編: M27、M28、M31 / この回でできるようになること: 陽解法・陰解法・Crank-Nicolson でオプション価格を計算し、安定条件を確かめることができる

0. 物語 ― 安定条件を破ると解が暴れる

オプションの値段を計算するとき、方程式を格子に載せて数値的に解く方法があります。格子を粗くすれば計算は速く終わります。では、どこまで粗くしても大丈夫でしょうか。

ある人が、株価100円、行使価格100円のコールオプションを計算しようとしました。空間方向を100分割、時間方向を206分割すると、答えは約10.45円になりました。正しい値とほぼ一致しています。

「時間の分割をもう少し減らしても大丈夫だろう」。そう考えて、時間方向を104分割に減らしました。計算量がほぼ半分になるのですから、悪くない判断に思えます。

ところが、出てきた答えは10の29乗を超える、でたらめに巨大な数でした。10.45円が11円や12円にずれるのではなく、完全に壊れたのです。

原因は、この計算法（陽解法と呼ばれます）に「安定条件」があるからです。時間の刻みがある上限を超えると、小さな誤差が1ステップごとに倍々に膨らみます。104分割では、誤差は毎ステップ約2倍に増え、104ステップ後には最初の誤差の10の31乗倍になっています。少し粗くしたのではなく、壊れる領域に足を踏み入れていたのです。

一方で、陰解法やクランク・ニコルソン法と呼ばれる方法なら、刻みをどれだけ粗くしても解は発散しません。その代わり、精度の出方が手法ごとに異なります。

この回では、3つの計算法それぞれについて、いつ安定でいつ壊れるのか、精度は何によって決まるのかを、誤差の増幅率の計算で明らかにします。

1. 動機 ― 方程式を数値で直接解く

この回を読み終えると、陽解法・陰解法・Crank-Nicolson でオプション価格を計算し、安定条件を確かめることができます。

ブラック・ショールズ方程式は偏微分方程式（PDE）です。第68回で閉じた解を得ましたが、あれはペイオフがコールやプットという単純な形のときに限った話でした。バリア付き、早期行使付き、ペイオフが複雑な場合には閉じた式がありません。

そこで PDE を格子（グリッド）に載せ、微分を差分に置き換えて数値的に解きます。これが有限差分法です。時間方向の進め方で「陽解法」「陰解法」「Crank-Nicolson 法」の3手法に分かれます。

手法の選び方を間違えると、解が発散して使い物にならなくなります。今回は、どの手法がどの条件で安定か、精度は何次かを定理として述べ、増幅率の計算で証明します。

この記事で証明すること

• THM-FD-1 陽解法は、時間刻みが空間刻みの二乗に比例する上限を超えると不安定になり、誤差が指数的に増大します。
• THM-FD-2 陰解法は時間刻みによらず安定ですが、精度は時間刻みの一次にとどまります。
• THM-FD-3 Crank-Nicolson 法は無条件安定で、時間と空間の両方について二次の精度を持ちます。ただし滑らかでない終端条件では振動が出ることがあります。

2. 設定と仮定 ― 格子を切って微分を差分にする

以下の仮定で計算を進めます。

• 対数価格 $x=\ln S$ に変換したブラック・ショールズ方程式を解く
• 境界条件は十分遠くに置く（標準偏差の4倍以上）
• 格子は等間隔
• 比較対象として閉じた式（BS 公式）を使う

対数価格への変換

ブラック・ショールズ方程式を $x=\ln S$ で変数変換すると、係数が定数になった放物型 PDE が得られます。$V(S,t)=u(x,t)$ とおくと

$$\frac{\partial u}{\partial t}+\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u}{\partial x^2}+\nu \frac{\partial u}{\partial x}-ru=0\text{\hspace{1em}(2.1)}$$

です。ここで $\nu =r-\frac{1}{2}{\sigma }^2$ です。変数変換の効果は、$S$ に依存していた係数が消えて定数になることです。これにより等間隔の格子が自然に使えます。

残存時間への読み替え

$\tau =T-t$（残存時間）に読み替えると、式 (2.1) は前進型になります。

$$\frac{\partial u}{\partial \tau }=\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u}{\partial x^2}+\nu \frac{\partial u}{\partial x}-ru\text{\hspace{1em}(2.2)}$$

$\tau =0$（つまり満期）での値がペイオフ $u(x,0)=\max (e^x-K,0)$ です。これを初期条件として $\tau$ を増やしていけば、現在の価格が求まります。

格子と差分近似

空間方向を $M$ 等分、時間方向を $N$ 等分します。空間刻みを $\Delta x$、時間刻みを $\Delta \tau$ とし、格子点の値を $u_j^n=u(x_j,n\Delta \tau )$ と書きます。

二階微分と一階微分を中心差分で近似します。

$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x{)}^2}\text{\hspace{1em}(2.3)}$$

$$\frac{\partial u}{\partial x}\approx \frac{u_{j+1}^n-u_{j-1}^n}{2\Delta x}\text{\hspace{1em}(2.4)}$$

中心差分はテイラー展開で確かめると、どちらも $(\Delta x{)}^2$ の精度を持ちます。

ここで次の無次元パラメータを導入します。

$$\alpha =\frac{{\sigma }^2\Delta \tau }{2(\Delta x{)}^2},\beta =\frac{\nu \Delta \tau }{2\Delta x}\text{\hspace{1em}(2.5)}$$

$\alpha$ は拡散の強さと格子の細かさの比で、安定性を支配する鍵となる量です。

3つの手法

陽解法（Explicit） は、時間微分を前進差分で近似し、空間微分を現在の時刻 $n$ で評価します。

$$\begin{array}{rl}{u}_j^{n+1}& =(\alpha -\beta )u_{j-1}^n\\ & +(1-2\alpha -r\Delta \tau )u_j^n\\ & +(\alpha +\beta )u_{j+1}^n\end{array}\text{\hspace{1em}(2.6)}$$

既知の値から直接次のステップが計算できるので、連立方程式を解く必要がありません。

陰解法（Implicit） は、空間微分を次の時刻 $n+1$ で評価します。

$$\begin{array}{rl}& -(\alpha -\beta )u_{j-1}^{n+1}\\ & +(1+2\alpha +r\Delta \tau )u_j^{n+1}\\ & -(\alpha +\beta )u_{j+1}^{n+1}=u_j^n\end{array}\text{\hspace{1em}(2.7)}$$

各ステップで三重対角の連立方程式を解く必要がありますが、トーマス法で $O(M)$ の計算量で済みます。

Crank-Nicolson 法 は、陽解法と陰解法の平均を取ります。空間微分を時刻 $n$ と $n+1$ の平均で評価します。

$$\begin{array}{rl}& -\frac{1}{2}(\alpha -\beta )u_{j-1}^{n+1}+(1+\alpha +\frac{r\Delta \tau }{2})u_j^{n+1}\\ & -\frac{1}{2}(\alpha +\beta )u_{j+1}^{n+1}\\ & =\frac{1}{2}(\alpha -\beta )u_{j-1}^n+(1-\alpha -\frac{r\Delta \tau }{2})u_j^n\\ & +\frac{1}{2}(\alpha +\beta )u_{j+1}^n\end{array}\text{\hspace{1em}(2.8)}$$

左辺は陰解法と同じ三重対角構造なので、やはりトーマス法で解けます。

境界条件

コールオプションでは、空間の左端（$S$ が十分小さい領域）で $u=0$、右端（$S$ が十分大きい領域）で $u\approx e^x-K{e}^{-r\tau }$ とします。右端では深い ITM のコールが本質的価値に近づくことを使っています。

以上で、PDE を格子上の代数的な漸化式に帰着させました。

3. 主張 ― 安定性と精度は手法で決まる

定理 THM-FD-1（陽解法の安定条件） 陽解法は、時間刻み $\Delta \tau$ が

$$\Delta \tau \le \frac{(\Delta x{)}^2}{{\sigma }^2}$$

すなわち $\alpha \le \frac{1}{2}$ を満たすとき安定である。この上限を超えると不安定になり、誤差が指数的に増大する。

定理 THM-FD-2（陰解法の無条件安定） 陰解法は $\Delta \tau$ の値によらず安定である。ただし精度は時間方向について $O(\Delta \tau )$（一次）、空間方向について $O((\Delta x{)}^2)$（二次）である。

定理 THM-FD-3（Crank-Nicolson の精度） Crank-Nicolson 法は $\Delta \tau$ の値によらず安定（無条件安定）であり、精度は時間方向 $O((\Delta \tau {)}^2)$、空間方向 $O((\Delta x{)}^2)$ である。ただし終端条件が滑らかでないとき（コールの折れ線ペイオフなど）、解のギリシャ文字に振動が生じうる。

3手法の特徴をまとめると、陽解法は実装が簡単だが安定条件が厳しく、陰解法は無条件安定だが精度が一次にとどまり、Crank-Nicolson は無条件安定かつ二次精度で実用上の標準です。

証明を読む（読まなくても記事は完結します）

4. 証明 ― 誤差の増幅率を調べる

この節は読まなくても記事は完結します。

この証明を追いたい人は、数学編 第27回「境界条件と終端条件」、数学編 第28回「変数変換で熱方程式に帰着する」、数学編 第31回「差分近似と誤差の次数」へ。

安定性と精度は、フーリエモードの増幅率（amplification factor）を調べることで証明します。式 (2.2) の $-ru$ の項を一旦無視し、拡散方程式

$$\frac{\partial u}{\partial \tau }=\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u}{\partial x^2}$$

で議論します。$-ru$ の項と移流項 $\nu u_x$ は結論を本質的に変えません。

THM-FD-1 の証明

格子上の誤差を ${\epsilon }_j^n$ と書きます。${\epsilon }_j^n$ も差分方程式 (2.6) と同じ漸化式に従います（差分方程式が線形だからです）。フーリエモード ${\epsilon }_j^n=g^n{e}^{ik{x}_j}$ を代入します。ここで $g$ が増幅率、$k$ は波数です。

式 (2.6) の簡略版（$\beta =0$、$r=0$）に代入すると

$$\begin{array}{rl}{g}^{n+1}{e}^{ik{x}_j}& =\alpha g^n{e}^{ik{x}_{j-1}}\\ & +(1-2\alpha )g^n{e}^{ik{x}_j}\\ & +\alpha g^n{e}^{ik{x}_{j+1}}\end{array}$$

が得られます。全体を $g^n{e}^{ik{x}_j}$ で割ると、$e^{\pm ik\Delta x}$ が出てきます。オイラーの公式 $e^{i\theta }+e^{-i\theta }=2\cos \theta$ を使って整理します。

$$\begin{array}{rl}g& =\alpha e^{-ik\Delta x}+(1-2\alpha )+\alpha e^{ik\Delta x}\\ & =1-2\alpha (1-\cos (k\Delta x))\end{array}\text{\hspace{1em}(4.1)}$$

$\cos (k\Delta x)$ は $-1$ から $1$ の範囲を取ります。$\cos (k\Delta x)=-1$ のとき $g$ は最小値を取り

$$g_{\min }=1-4\alpha \text{\hspace{1em}(4.2)}$$

となります。安定であるための条件は $|g|\le 1$ です。$g$ の最大値は $\cos (k\Delta x)=1$ のとき $g=1$ なので、問題は最小値だけです。

$g_{\min }\ge -1$ を要求すると $1-4\alpha \ge -1$、すなわち

$$\alpha \le \frac{1}{2}\text{\hspace{1em}(4.3)}$$

が得られます。$\Delta \tau$ に書き直すと $\Delta \tau \le (\Delta x{)}^2/{\sigma }^2$ です。

$\alpha >1/2$ のとき $|g_{\min }|>1$ となり、高波数の誤差成分が1ステップごとに $|g_{\min }|$ 倍に増幅されます。$N$ ステップ後の誤差は $|g_{\min }{|}^N$ に比例するので、指数的に発散します。

精度について。安定条件 $\alpha \le 1/2$ を満たすとき、テイラー展開による打ち切り誤差の評価から、陽解法は時間方向 $O(\Delta \tau )$（一次）、空間方向 $O((\Delta x{)}^2)$（二次）の精度を持ちます。安定条件を満たすには $\Delta \tau =O((\Delta x{)}^2)$ が必要なので、空間を $M$ 倍に細かくすると時間は $M^2$ 倍に細かくする必要があります。 ∎

THM-FD-2 の証明（要点）

陰解法では、フーリエモードを式 (2.7) の簡略版に代入すると

$$g(1+2\alpha (1-\cos (k\Delta x)))=1$$

が得られます。$1-\cos (k\Delta x)\ge 0$ なので分母は常に1以上です。したがって

$$g=\frac{1}{1+2\alpha (1-\cos (k\Delta x))}\text{\hspace{1em}(4.4)}$$

となり、$0<g\le 1$ が $\alpha$ の値によらず成り立ちます。これが無条件安定の意味です。

精度については、時間微分の前進差分が $O(\Delta \tau )$、空間の中心差分が $O((\Delta x{)}^2)$ なので、全体の精度は $O(\Delta \tau )+O((\Delta x{)}^2)$ です。時間方向が一次にとどまるのが陰解法の弱点です。 ∎

THM-FD-3 の証明（要点）

Crank-Nicolson 法は、陽解法と陰解法を重み $1/2$ ずつ平均した形です。増幅率は

$$g=\frac{1-\alpha (1-\cos (k\Delta x))}{1+\alpha (1-\cos (k\Delta x))}\text{\hspace{1em}(4.5)}$$

となります。分子と分母はどちらも実数で、分母は常に正です。分子の絶対値が分母を超えないことを確認します。$1-\cos (k\Delta x)=2{\sin }^2(k\Delta x/2)\ge 0$ を $\varphi$ とおくと

$$|g|=\left|\frac{1-\alpha \varphi }{1+\alpha \varphi }\right|$$

$\alpha \varphi \ge 0$ なので $|1-\alpha \varphi |\le 1+\alpha \varphi$ が常に成り立ちます（$\alpha \varphi \ge 0$ のとき $1-\alpha \varphi \le 1+\alpha \varphi$ かつ $-(1-\alpha \varphi )=\alpha \varphi -1\le 1+\alpha \varphi$）。よって $|g|\le 1$ が $\alpha$ によらず成り立ち、無条件安定です。

精度については、Crank-Nicolson は時間方向の中心差分（台形則）を使っているため $O((\Delta \tau {)}^2)$ の精度を持ちます。空間方向は中心差分で $O((\Delta x{)}^2)$ です。合わせて $O((\Delta \tau {)}^2)+O((\Delta x{)}^2)$ となり、空間と時間の両方が二次です。

ただし、終端条件に不連続な微分（コールペイオフの折れ線など）がある場合、初期の数ステップで高波数成分が正確に減衰されず、ギリシャ文字（特にデルタやガンマ）に振動が現れます。これは陽解法や陰解法の平均であることが原因で、陰解法だけを使えば高波数成分は強く減衰されます。対処法は後述します。 ∎

3つの手法の増幅率を比べると、陽解法だけが $|g|>1$ になる領域を持ち、その境界が安定条件 $\alpha \le 1/2$ です。

5. 反例 ― 刻みが大きいと発散する

反例 CEX-FD-1: 安定条件を破ると解が発散する

THM-FD-1 の安定条件 $\alpha \le 1/2$ を破るとどうなるかを数値で確認します。

パラメータは $S_0=100$、$K=100$、$r=0.05$、$\sigma =0.2$、$T=1$ です。空間を $M=100$ 分割すると $\Delta x\approx 0.016$ です。安定条件は

$$\Delta \tau \le \frac{(\Delta x{)}^2}{{\sigma }^2}\approx 0.0064$$

です。

• 安定な場合: $N=206$ とすると $\Delta \tau \approx 0.0049<0.0064$ で $\alpha <1/2$ です。ATM コール価格は $10.4499$ となり、BS 公式の $10.4506$ とほぼ一致します。
• 不安定な場合: $N=104$ とすると $\Delta \tau \approx 0.0096>0.0064$ で $\alpha \approx 0.75>1/2$ です。解は発散し、格子上の値の絶対値は $4.83\times 10^{29}$ という巨大な値に達します。

安定条件を $50\%$ 超えただけで解が完全に壊れます。増幅率 $|g_{\min }|=|1-4\times 0.75|=2$ となり、1ステップごとに高波数の誤差が2倍に増え、104ステップで $2^{104}\approx 10^{31}$ 倍になるからです。

安定条件は「少し超えたら少し悪くなる」ではなく、「超えたら壊れる」という性質のものです。

反例 CEX-FD-2: Crank-Nicolson のデルタに振動が出る

Crank-Nicolson は無条件安定ですが、コールの折れ線ペイオフのように微分が不連続な終端条件では、デルタ（$\Delta =\partial V/\partial S$）に振動が現れます。

$M=400$、$N=40$ の格子で、ATM（$S=100$ は格子点です）の前後の格子点でデルタを計算します。デルタは各格子点で中心差分から求め、同じ $S$ での BS のデルタと比べます。

$S$	CN のデルタ	BS のデルタ	誤差
98.41	0.6046	0.6064	$-0.0019$
98.81	0.6149	0.6141	$+0.0008$
99.20	0.6336	0.6217	$+0.0119$
99.60	0.6117	0.6293	$-0.0177$
100.00	0.6367	0.6368	$-0.0001$
100.40	0.6619	0.6443	$+0.0175$
100.80	0.6402	0.6517	$-0.0115$
101.21	0.6583	0.6591	$-0.0008$
101.61	0.6682	0.6664	$+0.0017$

ATM の前後で誤差が $+0.0119$、$-0.0177$、$+0.0175$、$-0.0115$ と符号を変えながら最大 $0.018$ に達し、隣接する格子点との間で値が上下に揺れています。$S=99.60$ から $S=100.40$ へ向かう途中でデルタが一度下がっており、本来は単調に増えるはずのデルタが単調でなくなっています。これは価格自体の精度が高くても、微分を取ると顕在化する問題です。

対処法: Rannacher スタートアップ

最初の数ステップ（ここでは4ステップ）だけ陰解法を使い、残りを Crank-Nicolson で進めます。陰解法はペイオフの折れ線を滑らかにする効果があります。

Rannacher スタートアップを適用すると、$97\le S\le 103$ の格子点でのデルタ誤差の最大値は $0.0177$ から $0.000088$ に縮小し、振動は消えます。$S=99.60$ と $S=100.40$ でのデルタはそれぞれ $0.6293$、$0.6443$ となり、BS 公式の値と小数4桁まで一致します。

反例 CEX-FD-3: 境界を近くに置くと価格がずれる

有限差分法では計算領域に境界を設ける必要があります。境界が原資産価格に近すぎると、境界条件の影響で価格が系統的にずれます。

Crank-Nicolson 法で $M=200$、$N=200$ とし、境界の位置を変えて比較します。

境界	下限 $S$	上限 $S$	価格	誤差
$\pm 1.5\sigma \sqrt{T}$	74.08	134.99	10.4358	$-0.0148$
$\pm 4\sigma \sqrt{T}$	44.93	222.55	10.4491	$-0.0014$

狭い境界では誤差が $0.0148$ と、広い境界の約10倍です。上側の境界 $S=135$ は ATM から1.35倍しか離れておらず、境界条件 $u\approx S-K{e}^{-r\tau }$ が不正確な領域の影響を受けるためです。

実務では $\pm 4\sigma \sqrt{T}$ 以上（つまり標準偏差の4倍以上）を確保するのが定番です。

6. 数値例 ― 3つの手法で誤差を比べる

パラメータは $S_0=100$、$K=100$、$r=0.05$、$\sigma =0.2$、$T=1$ の ATM コールです。BS 公式の価格は $10.4506$ です。

同一格子での比較

空間分割 $M=200$、境界は $\pm 4\sigma \sqrt{T}$ に置きます。

手法	時間分割 $N$	価格	誤差
陽解法	800	10.4505	$-0.0001$
陰解法	200	10.4439	$-0.0067$
Crank-Nicolson	200	10.4491	$-0.0014$

陽解法は安定条件を満たすために $N=800$ と多くのステップが必要ですが、精度は高くなっています。陰解法は $N=200$ で安定ですが、時間方向が一次精度のため誤差が大きめです。Crank-Nicolson は $N=200$ で陰解法より5倍近く正確です。

格子を細かくしたときの収束

Crank-Nicolson について、$M$ と $N$ を同時に倍にしていったときの誤差の縮み方を見ます。

$M$	$N$	誤差	収束次数
50	50	0.0231	―
100	100	0.0058	2.01
200	200	0.0014	2.00
400	400	0.0004	2.00

$M$ と $N$ を2倍にするたびに誤差がほぼ $1/4$ になっています。収束次数は2.0 で、THM-FD-3 の「時間・空間ともに二次」と整合しています。

Crank-Nicolson が実務の標準である理由は、この二次収束にあります。陰解法では $M$ と $N$ を倍にしても誤差が $1/2$ にしかならず、同じ精度を得るのに格子がはるかに細かくなります。

7. 次回へ ― モデルを市場に合わせる

有限差分法は PDE を数値で解く汎用的な道具です。今回は格子の構成、3手法の安定性と精度、そして刻みや境界の選び方が結果を左右することを見ました。

次回の第99回では、モデルのパラメータを市場データに合わせるキャリブレーションを扱います。ここまで所与としてきた $\sigma$ を、市場のオプション価格から逆算する作業です。

よくある質問

モンテカルロとどちらが良いのですか

次元が低いなら有限差分のほうが速く正確です。次元が高いとモンテカルロが有利になります。

安定とはどういう意味ですか

計算を進めても誤差が勝手に増えないことです。不安定だと小さな誤差が爆発します。

Crank-Nicolson が標準なのはなぜですか

安定で精度も高いからです。振動には対処法があります。

参考文献

1. Wilmott P.、Howison S. and Dewynne J.、The Mathematics of Financial Derivatives、Cambridge University Press、1995年
2. Tavella D. and Randall C.、Pricing Financial Instruments ― The Finite Difference Method、Wiley、2000年

---

数値例は compute.py により計算済み。BS 公式を比較対象として誤差を検証している。