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
|
#include <algorithm>
#include <iostream>
#include "rqmc.hpp"
double frac_part(double x){
return x - floor(x);
}
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;
};
template <typename Fct>
struct quasi_option : public generator<double>
{
quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), U(0,1), s(d), seed(d) {};
double operator()() {
std::vector<double> X(d);
double sum =0;
for(int i=0; i<d; i++){
seed[i]=U();
}
for (int i=0; i<n; i++){
X=s();
//std::cout<<X[i]<<std::endl;
for(int i=0; i<d; i++){
X[i] = gsl_cdf_gaussian_Pinv(frac_part(seed[i]+X[i]), 1);
}
sum += payoff(X);
//std::cout<<X[i]<<std::endl;
}
return sum/n;
};
private:
int n, d;
Fct payoff;
uniform U;
sobol s;
std::vector<double> seed;
};
int main(){
init_alea(1);
asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
int N= 10000;
int d =16;
typedef asian_option payoff_type;
std::vector<double> result(3);
result = monte_carlo(100, quasi_option<payoff_type> (N, d, A));
for(int i =0; i<3; i++){
std::cout<<result[i]<<std::endl;
}
return 0;
}
|
