diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | option.cpp | 80 | ||||
| -rw-r--r-- | rqmc.cpp | 37 | ||||
| -rw-r--r-- | rqmc.hpp | 29 |
4 files changed, 78 insertions, 70 deletions
@@ -17,7 +17,7 @@ test: test.o mt19937.o rqmc: rqmc.o mt19937.o $(CXX) $^ -o $@ $(GSL_FLAGS) -option: option.o rqmc.o mt19937.o +option: option.o mt19937.o $(CXX) $^ -o $@ $(GSL_FLAGS) clean: @@ -2,6 +2,10 @@ #include <iostream> #include "rqmc.hpp" +double frac_part(double x){ + return x - floor(x); +} + double pos (double x){ return x>0?x:0; } @@ -32,71 +36,47 @@ struct asian_option : public std::unary_function<std::vector<double>, double> -struct asian_option_qmc : public var_alea<double> + +template <typename Fct> +struct quasi_option : public generator<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++){ + quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), U(0,1), s(d), seed(d) {}; + + double operator()() { + 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)); + std::vector<double> X(d); + double sum =0; + X = s(); + for (int i=0; i<n; i++){ + X[i] = gsl_cdf_gaussian_Pinv(frac_part(seed[i]+X[i]), 1); } - double temp = std::accumulate(S.begin(), S.end(), 0.)/d; - return exp(-r*T)*pos(temp-K); + sum += payoff(X); + return sum/n; }; - private: - double r; - double T; - double S0; - double V; - int d; - double K; - sobol G; + 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 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); + int N= 1000; - std::vector<double> r = monte_carlo(N, A, G); - for(int i=0; i<3; i++){ - std::cout<<r[i]<<std::endl; + int d =16; + + typedef asian_option payoff_type; + for(int i =0; i<3; i++){ + std::cout<<monte_carlo(100, quasi_option<payoff_type> (N, d, A))[i]<<std::endl; } + return 0; } @@ -5,6 +5,7 @@ double frac_part(double x){ return x - floor(x); } + double mean_rqmc(int N, double X) { sobol s(1); double sum = 0; @@ -14,24 +15,24 @@ double mean_rqmc(int N, double X) { return sum/N; } -int main() { - init_alea(0); - int I=100; - int N= 10000; - uniform U; - double m = 0; - double s = 0; - double temp; +//~ int main() { + //~ init_alea(0); + //~ int I=100; + //~ int N= 10000; + //~ uniform U; + //~ double m = 0; + //~ double s = 0; + //~ double temp; - for(int i=0;i<I;i++){ - temp = mean_rqmc(N,U()); - m+=temp; - s+=temp*temp; - } - m = m/I; - s = s/I - m*m; + //~ for(int i=0;i<I;i++){ + //~ temp = mean_rqmc(N,U()); + //~ m+=temp; + //~ s+=temp*temp; + //~ } + //~ m = m/I; + //~ s = s/I - m*m; - std::cout<<"espérance "<<m<<" taille de l'IC "<<sqrt(s)*1.96/10<<std::endl; + //~ std::cout<<"espérance "<<m<<" taille de l'IC "<<sqrt(s)*1.96/10<<std::endl; - return 0; -} + //~ return 0; +//~ } @@ -31,7 +31,29 @@ compose(Fct f, VA X) { return compose_t<Fct, VA>(f, X); }; -// +template <typename VA> +struct mean_t: public generator<double> +{ + mean_t(int n, VA X): n(n), X(X) {}; + double operator()(){ + double sum = 0; + for(int i=0; i<n; i++){ + sum+=X(); + } + return sum/n; + } + private : + int n; VA X; + }; + +template <typename VA> +inline mean_t<VA> +mean(int n, VA X){ + return mean_t<VA>(n, X); +}; + + +//Les classes de monte-Carlo template <typename L> std::vector<double> monte_carlo(int n, L X) @@ -55,6 +77,9 @@ std::vector<double> monte_carlo(int n, Fct f, L X) return monte_carlo(n, compose(f, X)); }; +//Les classes de quasi-mean + + struct quasi_gaussian : public var_alea<std::vector<double> > { quasi_gaussian (int d, double mean = 0, double std =1) @@ -72,6 +97,8 @@ struct quasi_gaussian : public var_alea<std::vector<double> > sobol sob; }; + + struct gaussian_d : public var_alea<std::vector<double> > { gaussian_d(int d, double mean = 0, double std =1) |