【確率微分方程式】確率微分方程式を数値計算してみよう　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;
}