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
|
#include <vector>
#include <gsl/gsl_cdf.h>
#include "var_alea.hpp"
#include <algorithm>
#include <iostream>
#include "low_discrepancy.hpp"
template <typename L>
std::vector<double> monte_carlo(int n, L X)
{
std::vector<double> result(3,0);
double x;
for (int j = 0; j < n; j++) {
x = X();
result[0] += x;
result[1] += x*x;
}
result[0] /= (double) n;
result[1] = (result[1] - n*result[0]*result[0])/(double)(n-1);
result[2] = 1.96*sqrt(result[1]/(double) n);
return result;
}
double pos (double x){
return x>0?x:0;
}
double frac_part(double x){
return x - floor(x);
}
struct asian_option : public var_alea<double>
{
asian_option(double r, double T, double S0, double V, int d, double K)
: r(r), T(T), S0(S0), V(V), d(d), K(K), G(0,1) {};
double operator()() {
std::vector<double> S(d);
S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*G());
for(int i=1;i<d;i++){
S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*G());
}
double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
return exp(-r*T)*pos(temp-K);
};
private:
double r;
double T;
double S0;
double V;
int d;
double K;
gaussian G;
};
struct asian_option_qmc : public var_alea<double>
{
asian_option_qmc(double r, double T, double S0, double V, int d, double K)
: r(r), T(T), S0(S0), V(V), d(d), K(K), G(d), U(0,1) {};
double operator()() {
std::vector<double> S(d);
std::vector<double> sob(d);
sob = G();
S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(U()+sob[0]), 1));
for(int i=1;i<d;i++){
S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(U()+sob[i]), 1));
}
double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
return exp(-r*T)*pos(temp-K);
};
private:
double r;
double T;
double S0;
double V;
int d;
double K;
sobol G;
uniform U;
};
int main(){
init_alea(1);
asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
asian_option_qmc B(0.05, 1.0, 50.0, 0.1, 16, 45);
int M= 1000000;
int N=10000;
int I=100;
std::vector<double> meanvar = monte_carlo(M, A);
std::cout<<"espérance "<<meanvar[0] <<" IC "<<meanvar[2]<<std::endl;
std::vector<double> temp;
double m = 0;
double s = 0;
for(int i=0;i<I;i++){
temp = monte_carlo(N,B);
m+=temp[0];
s+=temp[0]*temp[0];
}
m = m/I;
s = s/I - m*m;
std::cout<<"espérance "<<m <<" IC "<<sqrt(s)*1.96/10<<std::endl;
return 0;
}
|
