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
|
#include <iostream>
#include <chrono>
#include <random>
#include <cmath>
#include <fstream>
#define PI 3.14159265
template<typename RG>
class NRG
{
RG rg;
public:
NRG(RG r):rg(r){}
double operator()();
};
template<typename NG>
class SPG
{
double S0, T, Mu, Vol;
NG ng;
public:
SPG(NG n, double Init, double Mat, double Drift, double Sig)
:ng(n),S0(Init),T(Mat),Mu(Drift),Vol(Sig){}
double operator()(double dt);
};
template<typename SPG>
double MCEngine(SPG s, double dt, int NumOfSims);
double Expected_Price(double Init, double Mat, double Drift);
int main()
{
std::mt19937 a(std::chrono::system_clock::now().time_since_epoch().count());
NRG<std::mt19937> b(a);
double S=50.0, maturity=1.0, mu=0.1, sig=0.3;
SPG<NRG<std::mt19937>> c(b, S, maturity, mu, sig);
double Mean = Expected_Price(S, maturity, mu);
std::ofstream d("1.csv");
d<<"dt"<<","<<"WeakError"<<"\n";
for(int i=5; i<=12; i++)
{
double dt = pow(2,-i);
d<<dt<<","<<std::abs(MCEngine(c, dt, 50000)-Expected_Price(S, maturity, mu))<<"\n";
}
d.close();
return 0;
}
double Expected_Price(double Init, double Mat, double Drift)
{
return Init*exp(Mat*Drift);
}
template<typename SPG>
double MCEngine(SPG s, double dt, int NumOfSims)
{
double sum=0.0;
for(int i=0; i<NumOfSims; i++)
{
sum+=s(dt);
}
return (double)sum/NumOfSims;
}
template<typename RG>
double class NRG::operator()()
{
static int flag = 0;
static double N2 = 0.0;
if(flag==0)
{
double u1 = (double)rg()/rg.max();
double u2 = (double)rg()/rg.max();
double N2 = sqrt(-2*log(u1))*cos(2*PI*u2);
flag=1;
return sqrt(-2*log(u1))*sin(2*PI*u2);
}
else
{
flag=0;
return N2;
}
}
template<typename NG>
double class SPG::operator()(double dt)
{
double t=0.0, s = S0, Mudt = Mu*dt, Volsqdt = Vol*sqrt(dt);
do
{
s *= (1+Mudt+Volsqdt*ng());
t += dt;
}while(t<T);
return s;
}
|