「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。

MENU

金融工学でのモンテカルロ法(08/N):分散減少法(3)

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

 モンテカルロ法

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

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

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

に分けられる。

5. 分散減少法

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


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

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

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

5.5 制御変量法

 評価対象関数f(x)から解析的に積分値が既知である関数g(x)を差し引いた部分にモンテカルロ法を適用する。そして得られた値に関数gの理論積分


\begin{aligned}
m_g=\int_0^1 g(x)dx
\end{aligned}

を加えて解とする。すなわち


\begin{aligned}
\hat{m}_3=\displaystyle{\frac{1}{N}}\sum_{i=1}^{N}\left\{f(\varepsilon_i)-g( (\varepsilon_i) )\right\} + m_g
\end{aligned}

で得られる。
 このとき解の誤差分散は


\begin{aligned}
V\left[\hat{m}_3\right]&=V\left[\displaystyle{\frac{1}{N}}\sum_{i=1}^{N}\left\{f(\varepsilon_i)-g( (\varepsilon_i) )\right\}\right]\\
&=\displaystyle{\frac{1}{N}}V[f(\varepsilon)-g(\varepsilon)]\\
&=\displaystyle{\frac{1}{N}} \left[\int_{0}^{1}\left\{f(X)-g(x)\right\}^2 dx -(m-m_g)^2 \right]
\end{aligned}

となる。
 制御変量法を適用する場合、制御変数の選択が問題となる。ただしオプションの場合は、Black-Scholes方程式を用いればよい。

5.6 例:算術平均ヨーロピアン・コール・オプションの価格評価

 算術平均ヨーロピアン・コール・オプションの価格を制御変量法で評価する。制御変数はBlack-Scholes方程式による解析解を用いる。
 第i番目の株価サンプルパスを


\begin{aligned}
\boldsymbol{S}_i&=(S_{i1},\cdots,S_{iM}),\\
S_{ij}&=S_{i(j-1)}+rS_{i(j-1)}\Delta t+\sigma\sqrt{\Delta t}S_{i(j-1)}\xi_{ij},\\
\boldsymbol{\xi}_i&=(\xi_{i1},\cdots,\xi_{iM}),\ \xi_{ij}\sim N(0,1)
\end{aligned}

で発生させる。またペイオフ


\begin{aligned}
A(\boldsymbol{S}_i)=e^{-rT}\max\left[\displaystyle{\frac{1}{M}\sum_{j=1}^{M}S_{ij}}-K,\ 0\right]
\end{aligned}

で与える。このとき、Black-Scholes方程式に則ったヨーロピアン・コール・オプション価格の解析解をC_{0}^{BS}(S(t),r,\sigma,T,K)とおけば、制御変量法による算術平均ヨーロピアン・コール・オプションの価格の数値解\hat{a}


\begin{aligned}
\hat{a}=\displaystyle{\frac{1}{N}}\left[\sum_{i=1}^{N}A(\boldsymbol{S}_i)-C(\boldsymbol{S}_i)\right]+C_{0}^{BS}(S(t),r,\sigma,T,K)
\end{aligned}

で与えられる。したがって

(1) あるiについて互いに独立な標準正規乱数列\boldsymbol{\xi}_{i}を発生させる。
(2) 株価サンプルパスを\boldsymbol{S}_i=(S_{i1},\cdots,S_{iM})生成する。
(3) A(\boldsymbol{S}_i)-C(\boldsymbol{S}_i)を計算する。
(4) 以上をi=1,\cdots,Nまで繰り返して
\begin{aligned}\left\{\displaystyle{\frac{1}{N}}\sum_{i=1}^{N}A(\boldsymbol{A}_i)-C(\boldsymbol{A}_i)\right\}\end{aligned}
を計算する。
(5) C_{0}^{BS}(S(t),r,\sigma,T,K)を加えて\hat{a}を算出する。

なお負の相関法と制御変量法を併用することも問題ない。その際は、制御変量法を先に適用すればよい。

5.6 5.5.のシミュレーション例

5.6.1 PCスペック情報
エディション Windows 10 Home
バージョン 20H2
プロセッサ Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM 8.00 GB
システムの種類 64 ビット オペレーティング システム、x64 ベース プロセッサ
Visual Studio Microsoft Visual Studio Community 2019 Version 16.8.2
.NET Framework Version 4.8.04084

5.6.2 コーディング例

 1,000回から100,000回まで1,000回ずつシミュレーション回数を増やしていったときの通常のMonte Carlo法と制御変量法の価格パスを調べることとする。

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;
        }

        /// <summary>
        /// Black-Scholes方程式に則ったヨーロピアン・オプション価格を返す。
        /// </summary>
        /// <param name="S">原資産価格</param>
        /// <param name="K">行使価格</param>
        /// <param name="r">無リスク金利</param>
        /// <param name="sigma">ボラティリティ</param>
        /// <param name="T">満期までの期間</param>
        /// <param name="CallPutFlg">コールプットフラグ。0ならばコール、1ならばプットの価格を返す。</param>
        /// <returns>ヨーロピアン・オプションの価格</returns>
        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 <= 100; i++)
            {
                TestNum.Add(i * 1000);
            }

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

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

            // 個々の結果
            double SimOP1 = 0.0; // 通常のMonte Carlo法
            double SimOP2 = 0.0; // 制御変数法

            double BSCall = OptionPricing(S0, K, r, std, T, 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++) {
                    Console.WriteLine("     Start " + (j + 1) + "/" + TestNum[i] + "-th test...");

                    var StockPath = new List<double>();

                    // 満期までの株価サンプルパスを生成する
                    for (int t = 0; t < T; t++)
                    {
                        if (t == 0)
                        {
                            StockPath.Add(S0);
                        }
                        else
                        {
                            double RN1 = obj.NextDouble();
                            double tmp = (1 + r + std * MathNet.Numerics.Distributions.Normal.InvCDF(0.0, 1.0, RN1)) * StockPath[t - 1];

                            StockPath.Add(tmp);
                        }
                    }

                    //ペイオフを計算
                    double A_S = Math.Exp(-r * T) * Max(StockPath.Average() - K, 0);
                    double C_S = Math.Exp(-r * T) * Max(StockPath.Last() - K, 0);

                    SimOP1 += A_S;
                    SimOP2 += A_S - C_S;
                }

                SimOP1 /= TestNum[i];
                SimOP2 /= TestNum[i];

                SimOP2 += BSCall;

                EachResult1.Add(SimOP1);
                EachResult2.Add(SimOP2);

                // 初期化
                SimOP1 = 0.0;
                SimOP2 = 0.0;
            }


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

            string Path = OutPath + "02_MonteCarloEvaluationByEuropeanAverageCallOptionPricing.csv";

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

                Stream.WriteLine(string.Format("No,TestNum, BS, Sim1, Sim2, Abs(Sim1), Abs(Sim2), Rel(Sim1), Rel(Sim2)"));

                for (int i = 0; i < EachResult1.Count(); i++)
                {
                    Stream.WriteLine(string.Format("{0},{1},{2},{3},{4},{5},{6},{7},{8}",
                    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));
                }

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

5.6.3 シミュレーション結果

 シミュレーション結果のグラフは以下の通り:

f:id:suguru_125:20211001073203p:plain

何度か繰り返してみたものの、形状に大きな相違はない。なお通常のMonte Carlo法はシミュレーション回数が少ないときに平均から大きく逸れているように見えるものの、実質的にはほとんど変わらない。
 通常のMonte Carlo法と制御変量法とで大きな相違があるとは思えず、むしろ通常のMonte Carlo法の方が精度が高そうである。ただし平均を取ってみても両者の水準が違う点、回数を増やしても収束が共に遅い点を踏まえると、実装ミスがあるかもしれない。

プライバシーポリシー お問い合わせ