金融工学でのモンテカルロ法(07/23)：分散減少法(2) - 「大人の教養・知識・気付き」を伸ばすブログ

　今回から、金融工学におけるシミュレーションについて学んでいく。テキストとして以下を使う。今回はP.59-P.67まで。

モンテカルロ法の金融工学への応用 (シリーズ現代金融工学)

作者:祥二, 湯前,輝好, 鈴木
朝倉書店

Amazon

power-of-awareness.com

5. 分散減少法
次回

5. 分散減少法

　 $\mathrm{Monte\ Carlo}$ 法は

問題に沿った（同時）分布に従う（多変量）乱数列の生成
その乱数列を使った計算

の2つの部分に分けて考えられる。乱数列の生成は更に

一様分布に従う乱数列 $\{u_1,u_2,\cdots,u_n\}$ の生成
それを元にした必要な同時分布に従う乱数列 $\{\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_n\}$ の生成

に分けられる。

　 $\mathrm{Monte\ Carlo}$ 法は高次元問題を取り扱うことはできるものの、誤差を減らすには大きな時間が掛かり時間効率が低い。そこで本章では時間効率を上げるための手段の1つとして分散減少法を説明する。 $\mathrm{Monte}$ $\mathrm{Carlo}$ 法における解の統計的誤差は点列数 $N$ に対して $k/\sqrt{N}$ に出来る。分散減少法はこの $k$ を減らす試みである。
　 $\mathrm{Monte\ Carlo}$ 法による関数 $f$ の積分値を

$\begin{aligned} E^{*}[f(X)] \end{aligned}$

と書くことにすると、 $\mathrm{Monte\ Carlo}$ 法のアイディアは以下の2つに分けることが出来ると言える：

関数 $f$ を同じ積分値を取る別の関数に置換え
期待値 $E^{*}[\cdot]$ の取り方を変更

関数 $f$ を同じ積分値を取る別の関数に置換え	*期待値 $E^{}[\cdot]$ の取り方を変更**
(1) 負の相関法	(1) 条件付き $\mathrm{Monte\ Carlo}$ 法
(2) 制御変量法	(2) 層別化法
(3) 回帰分析法	(2-1) ラテン・ハイパーキューブ法
(4) マルチンゲール分散法	(2-2) $\mathrm{Curran}$ の方法
	(3) 加重サンプリング法
	(4) 測度変換法

5.2　負の相関法

　負の相関法による積分

$\begin{aligned} m=\int_{0}^{1} f(x) dx \end{aligned}$

の解を $\hat{m}_2$ とする。負の相関法では関数 $f(x)$ を関数

$\begin{aligned} a(x)=\displaystyle{\frac{f(x)+f(1-x)}{2}} \end{aligned}$

で置き換える。このとき、

$\begin{aligned} \int_{0}^{1} f(x)dx=\int_{0}^{1} f(1-x)dx \end{aligned}$

に注意すれば

$\begin{aligned} \int_{0}^{1} a(x) dx&=\displaystyle{ \frac{1}{2}\int_{0}^{1}f(x) dx+\int_{0}^{1}f(1-x) dx}\\ &=\int_{0}^{1} f(x) dx \end{aligned}$

が成立する。したがって

$\begin{aligned} \hat{n}_2=\displaystyle{\frac{1}{N}\sum_{i=1}^{N}a(\varepsilon_i)} &=\int_{0}^{1} f(x) dx \end{aligned}$

で評価できる。
　負の相関法では、 $f(x)$ の評価回数 $2$ 倍になるので、同じ数の乱数を用いた場合、誤差評価が単純な $\mathrm{Monte\ Carlo}$ 法に対する比べて半分未満になれば効率的である。負の相関法における解の誤差分散は

$\begin{aligned} V[\hat{m}_2]&=V\left[\displaystyle{\frac{1}{N}\sum_{i=1}^{N}{\frac{f(\varepsilon_i)+f(1-\varepsilon_i)}{2}}}\right]\\ &=\displaystyle{\frac{1}{N}}V\left[\displaystyle{\frac{f(\varepsilon)+f(1-\varepsilon)}{2}} \right]\\ &=\displaystyle{\frac{1}{2N}}\left\{V[f(\varepsilon)] +Cov[f(\varepsilon),f(1-\varepsilon)] \right\} \end{aligned}$

であるから、もし

$\begin{aligned} Cov[f(\varepsilon),\ f(1-\varepsilon)]\lt 0 \end{aligned}$

ならば、負の相関法の単純な $\mathrm{Monte\ Carlo}$ 法に対する効果は

$\begin{aligned} \displaystyle{\frac{k_2}{k_1}}\lt \displaystyle{\frac{1}{2}} \end{aligned}$

となり、分散減少法が機能する。

5.3　例①：定積分

　定積分

$\begin{aligned} m^{q}=\int_{-\infty}^{\infty} f(x)q(x)dx \end{aligned}$

に対して負の相関法を適用する。ただし $q(\cdot)$ は区間 $(-\infty,\infty)$ を定義域とする任意の密度関数とし、これに対応する分布関数を $Q(\cdot)$ とし、その逆関数を $Q^{-1}(\cdot)$ は所与、乱数 $\theta\sim Q$ とする。
　変数変換 $x=Q^{-1}(u)$ により区間 $[0,1]$ の定積分

$\begin{aligned} m^{q}=\int_{0}^{1} f(Q^{-1}(u))du \end{aligned}$

と書き直す。これに負の相関法を適用して、

$\begin{aligned} m^{q}=\displaystyle{\frac{1}{N}\sum_{i=1}^{N}\displaystyle{\frac{f(\theta_i)+f({\theta_i}^{\prime})}{2}}} \end{aligned}$

として得られる。ここで $\theta_i=Q^{-1}(\varepsilon_i),\ \theta_i^{\prime}=Q^{-1}(1-\varepsilon_i)$ である。
　すなわち

(1)	$i=1, m^q=0$ とする。
(2)	乱数 $\varepsilon_i\sim U(0,1)$ を発生させる。
(3)	$\theta_i=Q^{-1}(\varepsilon_i)$ を得る。
(4)	$\theta_i^{\prime}=Q^{-1}(1-\varepsilon_i)$ を得る。
(5)	$\displaystyle(\frac{f(\theta_i)+f(\theta_i^{\prime})}{2})$ を計算し、 $m^q$ に加算する。
(6)	$i\lt N$ ならば $i=i+1$ として(2)に戻る。そうでなければ $m^q=\displaystyle{\frac{1}{N}}$ 倍する。
(7)	$\hat{m}_2^q=m^q$ とする。

5.4　例②：ヨーロピアン・コール・オプションのMonte Carlo法による評価

　ヨーロピアン・コール・オプション価格を単純な $\mathrm{Monte\ Carlo}$ 法および負の相関法で計算し、負の相関法の効果を評価する。なお行使価格は $K$ 、満期を $T$ とする。
　原資産価格 $S(t)$ の満たす確率微分方程式を

$\begin{aligned} \displaystyle{\frac{d S(t)}{S(t)}}=rdt+\sigma dz(t),\ \ r,\sigma\in\mathbb{R},\ \ 0\leq t\leq T \end{aligned}$

とする。このとき満期では

$\begin{aligned} \log\left(\displaystyle{\frac{S(t)}{S(0)}}\right)\sim N\left(\left(r-\displaystyle{\frac{1}{2}\sigma^2} \right)T, \sigma^2 T\right) \end{aligned}$

が成り立つ。したがって $\xi\sim N(0,1)$ として

$\begin{aligned} S(T)=S(0)\exp\left(\left(r-\displaystyle{\frac{1}{2}\sigma^2} \right)T+\sigma \sqrt{T}\xi\right):=S(\xi) \end{aligned}$

により満期株価を生成できる。
　これを基にペイオフ関数を

$\begin{aligned} C(x)=e^{-rT}\max[S(x)-K,0] \end{aligned}$

をもつヨーロピアン・コール・オプション価格

$\begin{aligned} c=\int_{-\infty}^{\infty}C(x)\phi(x)dx \end{aligned}$

を推計する。

5.5　ヨーロピアン・コール・オプションのシミュレーション結果

　実際にシミュレーションして負の相関法の精度を評価する。
　通常の $\mathrm{Monte\ Carlo}$ 法で評価したヨーロピアン・コール・オプションの価格と負の相関法で評価したヨーロピアン・コール・オプションの価格を $\mathrm{Black}$ - $\mathrm{Scholes}$ 方程式による解析解と比較した相対誤差を比較評価することにしたい。シミュレーション回数 $N$ は $N=10, 20,\cdots,10n,\cdots,100,000,\ n=1,2,\cdots,10,000$ と $10$ ずつ増やしていった。なお入力パラメータは以下の通りとする：

初期時点における原資産価格 $S_0$	$100$
行使価格 $K$	$80$
無リスク金利 $r$	$2\%/年$
ボラティリティ $\sigma$	$5\%/年$
満期までの期間 $T$	$50$ 日

5.5.1　実装

using System;
using System.Collections.Generic;
using System.Linq;
using System.IO;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics;

namespace MonteCarlo
{
    /// <summary>
    /// Monte Carlo Simulation
    /// </summary>
    public class MonteCarloSimulation
    {
        /// <summary>
        /// 2つの実数のうち大きい方を返す
        /// </summary>
        /// <param name="x">実数</param>
        /// <param name="y">実数</param>
        /// <returns>x>=yならばxを、そうでなければyを返す。</returns>
        private static double Max(double x, double y)
        {
            double ret = x + y;
            ret += Math.Abs(x - y);
            ret *= 0.5;
            return ret;
        }

        public static double OptionPricing(double S, double K, double r, double sigma, double T, int CallPutFlg)
        {
            if (CallPutFlg == 0 | CallPutFlg == 1)
            {
                double flg = -2.0 * CallPutFlg + 1.0;

                double d = Math.Log(S / K) + (r + 0.5 * Math.Pow(sigma, 2.0)) * T;
                d /= sigma * Math.Sqrt(T);

                double Price = S * MathNet.Numerics.Distributions.Normal.CDF(0.0, 1.0, flg * d);
                Price = flg * Price - flg * K * Math.Exp(-r * T) * MathNet.Numerics.Distributions.Normal.CDF(0.0, 1.0, flg * (d - sigma * Math.Sqrt(T)));

                return Price;
            }
            else
            {
                Console.WriteLine("CallPutFlg must be either 0 or 1.");
                return -99999;
            }
        }

        static void Main(string[] args)
        {
            // 出力先パス
            string OutPath = @"...適当な出力先...";

            // オプションの設定
            double S0 = 100;
            double r = 0.02; //年次
            double std = 0.05; //年次
            double K = 80;
            int T = 50; //日数

            // 日次に直す
            r /= 250;
            std /= Math.Sqrt(250);

            // モンテカルロ法の試行回数は変更する
            var TestNum = new List<double>();

            for (int i = 1; i <= 10000; i++)
            {
                TestNum.Add(i * 10);
            }

            var obj = new MathNet.Numerics.Random.MersenneTwister();

            // テスト結果
            List<double> EachResult1 = new List<double>();
            List<double> EachResult2 = new List<double>();

            List<double> EachEffect = new List<double>();
            double Effect1 = 0.0; //効果の分母
            double Effect2 = 0.0; //効果の分子

            // 個々の結果
            double sum1 = 0.0;
            double sum2 = 0.0;

            Console.WriteLine("Start Simulation Section ...");

            // Test実施
            for (int i = 0; i < TestNum.Count(); i++)
            {
                Console.WriteLine("   Start " + (i+1) +"/"+TestNum.Count()+"-th test...");

                for (int j = 0; j < TestNum[i]; j++)
                {
                    //満期原資産価格
                    double RN1 = obj.NextDouble();
                    RN1 = MathNet.Numerics.Distributions.Normal.InvCDF(0.0, 1.0, RN1);

                    double ST1 = S0 * Math.Exp((r - 0.5 * Math.Pow(std, 2.0)) * T + std * Math.Sqrt(T) * RN1);
                    double Payoff1 = Math.Exp(-r * T) * Max(ST1 - K, 0);
                    sum1 += Payoff1;

                    double ST2 = S0 * Math.Exp((r - 0.5 * Math.Pow(std, 2.0)) * T - std * Math.Sqrt(T) * RN1);
                    double Payoff2 = Math.Exp(-r * T) * Max(ST2 - K, 0);

                    sum2 += 0.5 * (Payoff1 + Payoff2);

                    Effect1 += Math.Pow(Payoff1, 2.0);
                    Effect2 += Math.Pow(0.5 * (Payoff1 + Payoff2), 2.0);
                }

                // Monte Carlo法による価格
                sum1 /= TestNum[i];
                sum2 /= TestNum[i];

                // 格納
                EachResult1.Add(sum1);
                EachResult2.Add(sum2);

                // 効果
                Effect1 -= TestNum[i] * Math.Pow(sum1, 2.0);
                Effect2 -= TestNum[i] * Math.Pow(sum2, 2.0);

                EachEffect.Add(Effect2 / Effect1);

                // 初期化
                sum1 = 0.0;
                sum2 = 0.0;

                Effect1 = 0.0;
                Effect2 = 0.0;
            }

            double BSCall = OptionPricing(S0, K, r, std, T, 0);

            // 出力
            Console.WriteLine("Start Output Section...");

            string Path = OutPath + "01_MonteCarloEvaluationByEuropeanCallOptionPricing.csv";

            try
            {
                StreamWriter Stream = new StreamWriter(Path, false, Encoding.UTF8);

                Stream.WriteLine(string.Format("No,TestNum, BS, C1, C2, Abs(C1), Abs(C2), Rel(C1),Rel(C2), Effect"));

                for (int i = 0; i < EachResult1.Count(); i++)
                {
                    Stream.WriteLine(string.Format("{0},{1},{2},{3},{4},{5},{6},{7},{8},{9}",
                    i + 1, TestNum[i], BSCall, EachResult1[i], EachResult2[i],
                    Math.Abs(EachResult1[i] - BSCall), Math.Abs(EachResult2[i] - BSCall),
                    Math.Abs(EachResult1[i] - BSCall) / BSCall, Math.Abs(EachResult2[i] - BSCall) / BSCall,
                    EachEffect[i]));
                }

                Stream.Close();
            }
            catch (Exception e)
            {
                Console.WriteLine(e.Message);
            }
        }
    }
}