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

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

MENU

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

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

 モンテカルロ法

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

の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.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倍になるので、同じ数の乱数を用いた場合、誤差評価が単純なモンテカルロ法に対する比べて半分未満になれば効率的である。負の相関法における解の誤差分散は


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

ならば、負の相関法の単純なモンテカルロ法に対する効果は


\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 例②:ヨーロピアン・コール・オプションのモンテカルロ法による評価

 ヨーロピアン・コール・オプション価格を単純なモンテカルロ法および負の相関法で計算し、負の相関法の効果を評価する。なお行使価格は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 ヨーロピアン・コール・オプションのシミュレーション結果

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

初期時点における原資産価格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);
            }
        }
    }
}
5.5.2 ヨーロピアン・コール・オプションのシミュレーション結果

 相対誤差は以下のグラフの通りとなった。相対誤差のオーダーが10^{-●}という形式であるから、常用対数を取った\log_{10}(相対誤差)y軸にとって表示している。


図1 シミュレーション回数を変えたときのBlack-Scholes方程式の解との相対誤差の推移
f:id:suguru_125:20210929024605p:plain

ここから明らかなように、負の相関法は10^2倍精度が良い。

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