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
|
#include <algorithm>
#include <iostream>
#include "rqmc.hpp"
double pos (double x){
return x>0?x:0;
}
struct asian_option : public std::unary_function<std::vector<double>, 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) {};
double operator()(std::vector<double> X) const {
std::vector<double> S(d);
S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[0]);
for(int i=1;i<d;i++){
S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[i]);
}
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;
};
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), seed(d) {
for(int i=0; i<d; i++){
seed[i]=U();
}
};
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(seed[0]+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(seed[i]+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;
std::vector<double> seed;
};
int main(){
init_alea(1);
asian_option A(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++){
//~ asian_option_qmc B(0.05, 1.0, 50.0, 0.1, 16, 45);
//~ 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;
int d = 16;
std::vector<double> test(d,1);
std::cout<<A(test)<<std::endl;
gaussian_d G(16,0,1);
std::vector<double> r = monte_carlo(N, A, G);
for(int i=0; i<3; i++){
std::cout<<r[i]<<std::endl;
}
return 0;
}
|
