pseudoquant

クオンツもどきの備忘録

バリアを連続n日間下回る確率

モンテカルロで計算するコードを適当に書いてみた。 面倒なので各種パラメータはベタ打ち。 確率過程は時間に依らないBSを使用

{ \displaystyle
\begin{equation}
 \frac{dS_t}{S_t} = (r - d) dt + \sigma dW_t
\end{equation}
}

出力は

  1. 毎日観測し、バリアを一度でも下回った確率
  2. {n}ごとに観測し、バリアを一度でも下回った確率
  3. 毎日観測し、{n}連続でバリアを下回った確率

である。

1.と2.の理論的な関係性は

http://www.columbia.edu/~sk75/mfBGK.pdf

で詳しく論じられている。

コード

BarrierProb.cpp

#include <iostream>
#include <fstream>
#include <cmath>

#include "Normals.h"
#include "mt19937.h"
#include "bsmc.h"

inline double payoffCall(const double spot, const double strike)
{
    return (spot > strike) ? spot - strike : 0.0;
}


// MCシミュレーションでバリアを連続n日間下回った確率を計算する
void simulateProb(const unsigned int referenceDays, 
                  const unsigned long nbSimulation,
                  const unsigned long seed,
                  double& probBreachBarrierLevel,
                  double& probBreachBarrierLevelInARow,
                  double& probBreachBarrierLevelInDescObs)
{
    /*--------------------Inputs--------------------*/

    //sde inputs
    const double vol = 0.3;
    const double spot = 100;
    const double rate = 0.01;
    const double dividend = 0.05;

    const double barrierLevel = 60;
    const double strike = 100;
    const double timeToMaturity = 1.0;
    
    // simulation conditions
    const unsigned int nbStep = 250;


    /*-----------------------------------------------*/
    
    // class bsmc
    bsmc bs(spot, strike, timeToMaturity, vol, rate, dividend, nbStep);

    // initialize random generator
    init_genrand(seed);

    // breaching observing parameter
    unsigned long nbBreached = 0;
    unsigned long nbBreachedInARow = 0;
    unsigned long nbBreachedInDescObs  = 0;

    //double runningPayoff = 0;
    double minSpot = spot;

    for (unsigned long iMC = 0; iMC < nbSimulation; ++iMC)
    {
        double thisSpot = spot;
        double minSpot = spot;
        double minSpotInDescObs = spot; // minimam spot in descrete observation

        // in a row parameters
        unsigned long breachDaysInARow = 0;
        unsigned long maxBreachDaysInARow = 0;


        for (unsigned long iStep = 0; iStep < nbStep; ++iStep)
        {
            double unitRand = genrand_real3();
            double normRand = InverseCumulativeNormal(unitRand);

            thisSpot *= bs.spotRatio(normRand,1);

            // check min spot
            if (thisSpot < minSpot) minSpot = thisSpot;

            // check if the spot breached the level in a row
            breachDaysInARow = (thisSpot <= barrierLevel) ? breachDaysInARow + 1: 0;
            maxBreachDaysInARow = std::max<unsigned long>(breachDaysInARow, maxBreachDaysInARow);

            // check if the spot breached the level in the descrete observation
            if (iStep % referenceDays == 0)
            {
                if (thisSpot < minSpotInDescObs) minSpotInDescObs = thisSpot;
            }
        }
        //runningPayoff += payoffCall(thisSpot, strike);
        if (minSpot <= barrierLevel) ++nbBreached;
        if (maxBreachDaysInARow >= referenceDays) ++nbBreachedInARow;
        if (minSpotInDescObs <= barrierLevel) ++nbBreachedInDescObs;
    }

    // results
    probBreachBarrierLevel          = static_cast<double>(nbBreached) / nbSimulation;
    probBreachBarrierLevelInARow    = static_cast<double>(nbBreachedInARow) / nbSimulation;
    probBreachBarrierLevelInDescObs = static_cast<double>(nbBreachedInDescObs) / nbSimulation;

    const double expectedPayoffCall = runningPayoff / nbSimulation;
    //const double callPrice = exp(-rate * timeToMaturity) * expectedPayoffCall;

    // std::cout << callPrice << std::endl;

}

int main()
{

    // MC Setup
    const unsigned long NB_SIMULATION = 1000000;
    const unsigned long SEED = 123456789;

    // results
    double probBreachBarrierLevel;
    double probBreachBarrierLevelInARow;
    double probBreachBarrierLevelInDescObs;

    // output
    std::ofstream ofs("test.csv");
    ofs << "refDays,continuous,inARow,descrete" << std::endl;

    for (unsigned int iRefDays = 1; iRefDays < 26; ++iRefDays)
    {
        simulateProb(iRefDays, NB_SIMULATION, SEED,
                     probBreachBarrierLevel, probBreachBarrierLevelInARow, probBreachBarrierLevelInDescObs);
        /*
        std::cout << "Breaching proberbilities..." << std::endl;
        std::cout << "In continuous observation : " << probBreachBarrierLevel << std::endl;
        std::cout << "In a row                  : " << probBreachBarrierLevelInARow << std::endl;
        std::cout << "In descrete observation   : " << probBreachBarrierLevelInDescObs << std::endl;
        */

        ofs << iRefDays << ","
            << probBreachBarrierLevel << "," 
            << probBreachBarrierLevelInARow << "," 
            << probBreachBarrierLevelInDescObs << std::endl;
    }

    std::cout << "finished" << std::endl;

    return 0;
}

bsmc.h

#ifndef BSMC_H
#define BSMC_H

class bsmc
{
private:
    double m_spot;
    double m_strike;
    double m_timeToMaturity;
    double m_vol;
    double m_rate;
    double m_div;

    double m_dt;
    double m_drift;
    double m_variance;
    double m_rootVariance;

public:
    // constructor
    bsmc(const double spot, const double strike, const double timeToMaturity,
         const double vol,  const double rate,   const double div,
         const unsigned long timeStepPerYear);
    double spotRatio(const double normRand, const unsigned long numStep);
    double bs();
};
#endif

bsmc.cpp

#include "bsmc.h"
#include <cmath>

bsmc::bsmc(const double spot, const double strike, const double timeToMaturity,
           const double vol,  const double rate,   const double div,
           const unsigned long timeStepPerYear)
:
        m_spot(spot), m_strike(strike), m_timeToMaturity(timeToMaturity),
        m_vol(vol), m_rate(rate), m_div(div)
        {
            m_dt = 1.0 / timeStepPerYear;
            m_variance = m_vol*m_vol*m_dt;
            m_rootVariance = sqrt(m_variance);
            m_drift = (m_rate - m_div) * m_dt - 0.5 * m_variance;
        }

double bsmc::spotRatio(const double normRand, const unsigned long numStep)
{
    return exp(m_drift*numStep + m_rootVariance * sqrt(static_cast<double>(numStep)) * normRand);
}

double bsmc::bs()
{
    //const double d1 = log(m_spot/m_strike)/m_rootVariance + (rate - div
    return 0.0;
}

出力結果

refDays,continuous,inARow,descrete

1,0.130531,0.130531,0.130531

2,0.130531,0.121606,0.126122

3,0.130531,0.115278,0.124553

4,0.130531,0.110094,0.121491

5,0.130531,0.105718,0.117017

6,0.130531,0.101822,0.116409

7,0.130531,0.098376,0.114105

8,0.130531,0.095221,0.11548

9,0.130531,0.092248,0.109852

10,0.130531,0.089476,0.106233

11,0.130531,0.086889,0.10691

12,0.130531,0.084399,0.104423

13,0.130531,0.082076,0.10937

14,0.130531,0.079867,0.100869

15,0.130531,0.077822,0.101782

16,0.130531,0.075816,0.101202

17,0.130531,0.073968,0.098742

18,0.130531,0.072163,0.094722

19,0.130531,0.070421,0.104994

20,0.130531,0.068721,0.098654

21,0.130531,0.067122,0.090496

22,0.130531,0.065568,0.098979

23,0.130531,0.064088,0.088532

24,0.130531,0.062648,0.096315

25,0.130531,0.061228,0.083589

Mersenne Twisterまわり

Excel上でMersenne Twisterを使いたかったので、探してみたところ、意外にExcelVBA上で使用できるMersenne Twisterのライブラリやコードが見つかった。

まず、DLLは和田さんのページ

http://www001.upp.so-net.ne.jp/isaku/

から潜っていって、

http://www001.upp.so-net.ne.jp/isaku/rand2.html

に一通りある。

さらに数理統計研が各種ライブラリを用意してくれている。

http://random.ism.ac.jp/info01_e.html

こちらは"NTRAND"というxll

http://www.ntrand.com/

 

VBAのコードはここから。

http://www001.upp.so-net.ne.jp/isaku/mt.html

 

なお、松本眞先生の記事は非常に面白かった。

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/index.html

VC++でDLLを作成する方法

まだあまり理解していないけど、参考になるリンクをまとめる。
ある程度わかってきたらまた更新する予定。

 

VC++2005でDLL作成、呼び出し、VBAからDLL呼び出し | 夏研ブログ

sites.google.com

 

SAFEARRAYについては以下の記事。

www26.atwiki.jp

www26.atwiki.jp

このページの8にも書いてある。

Visual Basic Web Site

xllについてはProfessional Excel Developmentにも書いてあるようだ。

Professional Excel Development: The Definitive Guide to Developing Applications Using Microsoft Excel, VBA, and .NET (2nd Edition): Rob Bovey, Dennis Wallentin, Stephen Bullen, John Green: 9780321508799: Amazon.com: Books

テスト

投稿のテスト

 

本文1
本文続き

箇条書きスタート

  • 箇条書き1
  • 箇条書き2
  • 箇条書き3

数式入力

数式環境

\begin{equation}
     V_t = F_t N(d_1) - K N(d_2)
\end{equation}

これはどうか