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