数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【確率微分方程式】確率微分方程式を数値計算してみよう C++コード付き

はじめに

今回は確率微分方程式を数値的に計算するC++コードを紹介しようと思います。パスがランダムに進展していく様は見ていて非常に興味深いです。

確率微分方程式とは?

一般論

 dX_t = f(t,X_t)dt+ g(t,X_t)dW_t

のような式のことを伊藤の確率微分方程式(Stochastic Differential Equation、SDE)と呼んでいます。ここで、 X が時刻  t における未知数(例えば株価や生物の個体数など)、 t は時間、 f はドリフト係数、 g は拡散係数(ボラティリティ)、 W_t は1次元の標準ブラウン運動で、この項が確率論的なゆらぎをあらわしています。標準ブラウン運動は、平均0、分散tの正規分布に従います。この事実を使うと、確率微分方程式を数値的に解くことができるようになります。

幾何ブラウン運動

今回扱う具体的な確率微分方程式

 dX_t = \mu X_t dt+ \sigma X_t W_t

としましょう。ここで、 f(t,X_t) = \mu X_t g(t,X_t) = \sigma X_t と設定しています。このようにドリフトと拡散係数が線型で  X_t を含むものを幾何ブラウン運動と呼びます。ウィナー過程とも呼ばれています。幾何ブラウン運動は数理ファイナンス金融工学の分野で頻出します。例えば、ブラック・ショールズ方程式で使われています。幾何ブラウン運動の利点は線型で比較的解析が簡単であることと解析解が求まることです。上記の幾何ブラウン運動の解析解は

 X_t = X_0 \exp\left( \left( \mu-\frac{\sigma^2}{2} \right) + \sigma dW_t \right)

とあらわすことができます。ここで、 X_0 は初期値です。

オイラー・丸山スキー

オイラー・丸山スキームは、確率微分方程式に対するもっともシンプルな数値計算手法です。常微分方程式におけるオイラー法と対応しています。収束次数は1/2次となかなかしょっぱいです。オイラー・丸山スキームは確率微分方程式

 X_{n+1} = f(t_n,X_n)\Delta t+ g(t_n,X_n)\Delta W_n

のように近似します。ここで、 n はタイムステップ、 \Delta t は時間の刻み幅、 \Delta W_n= q\sqrt{\Delta t} です。 q は標準正規分布に従い、平均0、分散1を満たすものとしています。つまりブラウン運動の増分は、平均0、分散  \Delta t正規分布に従います。ここがとても大切です。数値計算では  q、すなわち標準正規分布を発生させればよいのです。

計算結果

パラメータとしては、 \mu=1 \sigma=\sqrt{2} X_0=1 \Delta t=0.0001、計算終了時刻は1としました。乱数はメルセンヌツイスターを使っています。計算結果を見ていきましょう。数値解と解析解の両方をプロットしています。青が数値解で、赤が解析解です。解析解は数値解と同じブラウン運動を使って計算します。

f:id:mutsumunemitsutan:20190926223450p:plain

このように、ブラウン運動において、パスはギザギザしており微分不可能ですが、確率1で連続となります。値が増えたり減ったりランダムな動きをしていることがわかると思います。赤と青の線はほぼ重なっており、ちゃんと確率微分方程式が解けていることがわかりました。シード値を変えることにより異なるパスが実現します。いろいろ変えて遊んでみてください。

f:id:mutsumunemitsutan:20190926223240p:plain

C++コード

#include <random>
#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;

inline double drift(double x, double t, double mu)
{
    return x*mu;
}

inline double diffusion(double x, double t, double sigma)
{
    return x*sigma;
}

int main()
{
    int seed=80;
	mt19937 mt(seed);
    normal_distribution<> gauss(0.0, 1.0);

    ofstream file;
    file.open("path.txt");

    double mu=1.0;
    double sigma=sqrt(2.0);
    double x=1.0;
    double dt=0.0001;
    double TotalTime=1.0;
    double t;
    double exact=1.0;
    double random;
    double W=0.0;

    file<<0.0<<" "<<x<<" "<<exact<<endl;
    
	for(int i=0;dt*double(i)<TotalTime;++i)
	{
        t=dt*double(i);
        random=gauss(mt);
        W+=random*sqrt(dt);
        x=x+drift(x,t,mu)*dt+diffusion(x,t,sigma)*sqrt(dt)*random;
        exact=exp(sigma*W);

	    file<<t<<" "<<x<<" "<<exact<<endl;
	}

    return 0;
}

参考文献

微分方程式による計算科学入門』の第4章を参考にしています。この本は確率微分方程式の数値解法に関して解説してある数少ない和書で重宝しています。

あとは毛利光希氏によるこちらの資料です。
http://www-mmc.es.hokudai.ac.jp/~masakazu/Software/Manual/SDE_Solver/Manual.pdf

幾何ブラウン運動の解析解の導出は『確率微分方程式入門』のpp.77-78の例題4.3に載っています。伊藤の公式を適用するだけです。