SuperMC/mcengine.cpp at master · Glikstein/SuperMC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
// Copyright (c) 2014 Ezequiel Dan Glikstein


// This file both implements MCEngine and defines various price-path modelling functions
#include "mcengine.h"

#include "normals.h"
#include "mrg32k3a.h"

///////////////////////////////////
///// MCEngine IMPLEMENTATION /////
///////////////////////////////////

// constructors

// default constructor
MCEngine::MCEngine(double intervalSize, long numberofSteps, long numberofSimulations, double (*ptrPayoff)(StepData&), StepData& initialStep )
: m_intervalSize(intervalSize), m_numberofSteps(numberofSteps), m_numberofSimulations(numberofSimulations), m_currentStep(initialStep), m_lastStep(initialStep), m_initStep(initialStep)
{
	payoff = ptrPayoff;
	model = Models::BrownianStep;	// set brownian motion by default
	m_sum = m_sqrSum = 0;
	D = 1.0;	// discounting factor at time 0
	m_stepSize =  intervalSize / (double)numberofSteps;
	exotic = Exotics::Vanilla;
}

// custom model constructor
MCEngine::MCEngine(double intervalSize, long numberofSteps, long numberofSimulations, double (*ptrPayoff)(StepData&), StepData& (*ptrModel)(std::vector<double>&, MCEngine&), StepData& initialStep )
: m_intervalSize(intervalSize), m_numberofSteps(numberofSteps), m_numberofSimulations(numberofSimulations), m_currentStep(initialStep), m_lastStep(initialStep), m_initStep(initialStep)
{

	payoff = ptrPayoff;
	model = ptrModel;
	m_sum = m_sqrSum = 0.0;
	D = 1.0;	// discounting factor at time 0
	m_stepSize =  intervalSize / (double)numberofSteps;
	exotic = Exotics::Vanilla;
	std::cout<<"Constructed Engine with: S0 = "<<m_initStep.S<<" Maturity: "<<m_intervalSize<<" vol = "<<m_initStep.vol<<" r = "<<m_initStep.r<<" stepsize = "<<m_stepSize<<" number of steps = "<<m_numberofSteps<<std::endl;

}

// custom model constructor for path-dependent options
MCEngine::MCEngine(double intervalSize, long numberofSteps, long numberofSimulations, double (*ptrPayoff)(StepData&), StepData& (*ptrModel)(std::vector<double>&, MCEngine&), void (*ptrExotic)(StepData&), StepData& initialStep )
: m_intervalSize(intervalSize), m_numberofSteps(numberofSteps), m_numberofSimulations(numberofSimulations), m_currentStep(initialStep), m_lastStep(initialStep), m_initStep(initialStep)
{

	payoff = ptrPayoff;
	model = ptrModel;
	exotic = ptrExotic;
	m_sum = m_sqrSum = 0.0;
	D = 1.0;	// discounting factor at time 0
	m_stepSize =  intervalSize / (double)numberofSteps;
	std::cout<<"Constructed Engine with: S0 = "<<m_initStep.S<<" Maturity: "<<m_intervalSize<<" vol = "<<m_initStep.vol<<" r = "<<m_initStep.r<<" stepsize = "<<m_stepSize<<" number of steps = "<<m_numberofSteps<<std::endl;
	std::cout<<"m_lastStep.S = "<<m_lastStep.S<<" .vol = "<<m_lastStep.vol<<" .volvol =  "<<m_lastStep.volvol<<" .r = "<<m_lastStep.r<<std::endl;
}


MCEngine::~MCEngine()
{
	// no mem alloc'd
}


// calculators


void MCEngine::Simulate()
{

	using namespace boost::math;
/*
	// setting up the RNG using Boost

	boost::mt19937 RNG1;
	boost::normal_distribution<> normal(0,1);
	boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > N1	(RNG1, normal);

	RNG1.seed(100488);	// fixed seed to obtain reproducible results
*/
	RandomStream UnifStream;	// U(0,1) stream


	std::cout<<"Simulating "<<m_numberofSimulations<<" price paths...\n";

	// setting up the randoms vector
	std::vector<double> randoms;
	randoms.resize(Models::MAX_RV);
	for(int r = 0; r < Models::MAX_RV; r++)
		randoms[r] = invertNormal(UnifStream.nextVar()); //N1();

	// simulate!
	for(long sim = 0; sim < m_numberofSimulations; sim++)
	{
		for(long i = 0; i<m_numberofSteps; i++)
		{
			for(int r = 0; r<Models::MAX_RV; r++)
			{
				if(randoms[r] == 0.0)
					randoms[r] = invertNormal(UnifStream.nextVar()); //N1();	// only generate what you need
				//std::cout<<"randoms["<<r<<"] = "<<randoms[r]<<'\t';
			}
			//std::cout<<'\n';
			model(randoms, *this);	// advance one time unit along the price path
			exotic(m_lastStep);		// update the state of m_lastStep
		}
		// gather staticstics
		double payout = payoff(m_lastStep) * D;
		m_sum += payout;
		m_sqrSum += payout*payout;

		Reset();	// reset steps, discount factor
	}
}


// getters

std::vector<double> MCEngine::getStats() const
{
	std::vector<double> foo;
	foo.push_back( getMean());
	foo.push_back( getVariance());
	foo.push_back( getStdDev());
	foo.push_back( getStdError());
	return foo;
}

double MCEngine::getStdDev() const
{
	double M = m_numberofSimulations; // for readability
	return 1.0/(M-1.0) * sqrt( m_sqrSum - 1.0/M * m_sum*m_sum );
}