From 6b07d76c512c5450ce34fd6f2f8483c395a24ea2 Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Tue, 16 Feb 2016 22:01:28 +0000
Subject: meilleur encapsulage pour asian quasi_option

---
 Makefile   |  2 +-
 option.cpp | 80 +++++++++++++++++++++++---------------------------------------
 rqmc.cpp   | 37 +++++++++++++++--------------
 rqmc.hpp   | 29 ++++++++++++++++++++++-
 4 files changed, 78 insertions(+), 70 deletions(-)

diff --git a/Makefile b/Makefile
index b3dd461..3e093eb 100644
--- a/Makefile
+++ b/Makefile
@@ -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:
diff --git a/option.cpp b/option.cpp
index fb8b9cc..8229730 100644
--- a/option.cpp
+++ b/option.cpp
@@ -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;
 }
diff --git a/rqmc.cpp b/rqmc.cpp
index bca3c77..7109d82 100644
--- a/rqmc.cpp
+++ b/rqmc.cpp
@@ -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;
+//~ }
diff --git a/rqmc.hpp b/rqmc.hpp
index 62efbaa..415ab00 100644
--- a/rqmc.hpp
+++ b/rqmc.hpp
@@ -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)
-- 
cgit v1.2.3-70-g09d2