ajout option rqmc et correction rqmc.cpp

2 files changed, 50 insertions, 4 deletions

diff --git a/option.cpp b/option.cpp
index deb8347..016de0b 100644
--- a/option.cpp
+++ b/option.cpp
@@ -3,6 +3,7 @@
 #include "var_alea.hpp"
 #include <algorithm>
 #include <iostream>
+#include "low_discrepancy.hpp"
 
 
 template <typename L>
@@ -25,6 +26,10 @@ 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>
 {
 	asian_option(double r, double T, double S0, double V, int d, double K)
@@ -50,15 +55,56 @@ struct asian_option : public var_alea<double>
         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) {};
+        
+	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));
+        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));
+        }
+        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;
+};
 
 
 int main(){
     init_alea(1);
     asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
-    int N=1000000;
+    asian_option_qmc B(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(N, A);
+    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;
 
     
     return 0;
diff --git a/rqmc.cpp b/rqmc.cpp
index e8b98fc..6e5acd3 100644
--- a/rqmc.cpp
+++ b/rqmc.cpp
@@ -31,8 +31,8 @@ int main() {
         m+=temp;
         s+=temp*temp;
     }
-    m = m/N;
-    s = s/N - m*m;
+    m = m/I;
+    s = s/I - m*m;
     
     std::cout<<"espérance "<<m<<" taille de l'IC "<<sqrt(s)*1.64/10<<std::endl;