From cee0800329a0be3b27ae3907498ddf0883393382 Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Mon, 15 Feb 2016 19:16:29 +0000
Subject: création de nouvelles classes, composition de monte_carlo...
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 Makefile   |   2 +-
 option.cpp | 105 ++++++++++++++++++++++++++++---------------------------------
 rqmc.cpp   |   7 ++---
 rqmc.hpp   |  96 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 147 insertions(+), 63 deletions(-)
 create mode 100644 rqmc.hpp

diff --git a/Makefile b/Makefile
index 3e093eb..b3dd461 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 mt19937.o
+option: option.o rqmc.o mt19937.o
 	$(CXX) $^ -o $@ $(GSL_FLAGS) 
 
 clean:
diff --git a/option.cpp b/option.cpp
index 016de0b..fb8b9cc 100644
--- a/option.cpp
+++ b/option.cpp
@@ -1,72 +1,53 @@
-#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;
-}
+#include "rqmc.hpp"
 
 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>
+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), G(0,1) {};
-        
-	double operator()() {
+    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)*G());
+        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)*G());
+            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:
+        };
+        
+    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) {};
+        : 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(U()+sob[0]), 1));
+        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(U()+sob[i]), 1));
+            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);
@@ -81,30 +62,40 @@ struct asian_option_qmc : public var_alea<double>
         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);
-    asian_option_qmc B(0.05, 1.0, 50.0, 0.1, 16, 45);
-    int M= 1000000;
-    int N=10000;
-    int I=100;
+    //~ 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;
+    //~ 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;
diff --git a/rqmc.cpp b/rqmc.cpp
index 6e5acd3..bca3c77 100644
--- a/rqmc.cpp
+++ b/rqmc.cpp
@@ -1,7 +1,4 @@
-#include "low_discrepancy.hpp"
-#include <vector>
-#include <gsl/gsl_cdf.h>
-#include "var_alea.hpp"
+#include "rqmc.hpp"
 
 
 double frac_part(double x){
@@ -34,7 +31,7 @@ int main() {
     m = m/I;
     s = s/I - m*m;
     
-    std::cout<<"espérance "<<m<<" taille de l'IC "<<sqrt(s)*1.64/10<<std::endl;
+    std::cout<<"espérance "<<m<<" taille de l'IC "<<sqrt(s)*1.96/10<<std::endl;
 
     return 0;
 }
diff --git a/rqmc.hpp b/rqmc.hpp
new file mode 100644
index 0000000..62efbaa
--- /dev/null
+++ b/rqmc.hpp
@@ -0,0 +1,96 @@
+#include "low_discrepancy.hpp"
+#include <vector>
+#include <gsl/gsl_cdf.h>
+#include "var_alea.hpp"
+
+
+double frac_part(double x);
+
+//fonctions de compositions de monte_carlo.hpp
+
+template <class _Result>
+struct generator 
+{
+	typedef _Result result_type;
+};
+
+template <typename Fct, typename VA>
+struct compose_t : public generator< typename Fct::result_type >
+{
+	compose_t(Fct f, VA X) : f(f), X(X) {};
+	typename Fct::result_type operator()() { 
+		return f(X()); 
+	};
+	private:
+		Fct f; VA X;
+};
+
+template <typename Fct, typename VA>
+inline compose_t<Fct, VA>
+compose(Fct f, VA X) {
+    return compose_t<Fct, VA>(f, X);
+};
+
+//
+
+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;
+};
+
+template <typename L, typename Fct>
+std::vector<double> monte_carlo(int n, Fct f, L X)
+{
+	return monte_carlo(n, compose(f, X)); 
+};
+
+struct quasi_gaussian : public var_alea<std::vector<double> >
+{
+    quasi_gaussian (int d, double mean = 0, double std =1)
+        : d(d), mean(mean), std(std), sob(d) {};
+    std::vector<double> operator ()(){
+        std::vector<double> result = sob();
+        for(int i=0;i<d; i++){
+            result[i] = mean + gsl_cdf_gaussian_Pinv(result[i], std);
+        }
+        return value = result;
+    }
+    private: 
+    int d;
+    double mean, std;
+    sobol sob;
+};
+
+struct gaussian_d : public var_alea<std::vector<double> >
+{
+    gaussian_d(int d, double mean = 0, double std =1)
+        : d(d), mean(mean), std(std), G(mean, std) {};
+    std::vector<double> operator ()(){
+        std::vector<double> result(d);
+        for(int i=0;i<d; i++){
+            result[i] = G();
+        }
+        return value = result;
+    }
+    private: 
+    int d;
+    double mean, std;
+    gaussian G;
+};
+
+
+
+
+
+    
-- 
cgit v1.2.3-70-g09d2