création de la struct f_mu

5 files changed, 66 insertions, 79 deletions

diff --git a/src/main.cpp b/src/main.cpp
index f6e998a..561b8c2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,6 +6,7 @@
 #include "stratified_sampling.hpp"
 #include <cmath>
 #include <algorithm>
+#include "opti.hpp"
 
 using namespace std;
 //--génération quantiles--
@@ -43,16 +44,30 @@ x = 1.64*S.estimator().second;
 };
 
 
+ void exemple2 (){
+     std::vector<double> mu(16);
+     mu = argmax(0.05, 1.0, 50, 0.1, 45, 16);
+     double norm_mu = 0;
+     std::vector<double> u(16);
+     for(int i=0; i<16; i++) {
+         norm_mu += mu[i]*mu[i];
+         u[i] = mu[i]/norm_mu;
+     }
+     vector<double> q = quantile_norm(100, 1);
+     vector<double> p(100, 0.01);
+     asian_option A(0.05, 1.0, 50, 0.1, 45, 16);
+     f_mu G(mu,A);
+     multi_gaussian_truncated MG(q[50],q[51], u);
+     for(int i=0; i<10; i++){
+        std::cout<<G(MG())<<std::endl;
+    }
+}
+
 
 int main()
 {   
-    std::vector<double> u {sqrt(0.2), sqrt(0.2), sqrt(0.2), sqrt(0.2), sqrt(0.2)};
-    multi_gaussian_truncated G(0, 2, u);
-    std::vector<double> r(5);
-    r = G();
-    for (int i=0; i<5; i++){
-        std::cout<<r[i]<<std::endl;
-    }
+    
+    exemple2();
         
     return 0; 
 }
diff --git a/src/opti.cpp b/src/opti.cpp
index 702e8c2..23ddff3 100644
--- a/src/opti.cpp
+++ b/src/opti.cpp
@@ -1,8 +1,4 @@
-#include <vector>
-#include <nlopt.hpp>
-#include <cmath>
-#include <algorithm>
-#include <iostream>
+#include "opti.hpp"
 
 double pos (double x){
     return x>0?x:0;
@@ -30,17 +26,17 @@ double f (const std::vector<double> &X, std::vector<double> &grad, void *params)
 
 
 
-int main() {
+std::vector<double> argmax(double r, double T, double S0, double V, double K, int d){
     
-    double params[5] = {0.05, 1, 50, 0.1, 45};
+    double params[5] = {r, T, S0, V, K};
     
-    nlopt::opt opt(nlopt::LN_COBYLA, 16);
+    nlopt::opt opt(nlopt::LN_COBYLA, d);
 
     opt.set_max_objective(f, &params);
 
     opt.set_xtol_rel(1e-4);
 
-    std::vector<double> x(16,0);
+    std::vector<double> x(d,0);
     std::vector<double> g(0);
 
     
@@ -48,12 +44,8 @@ int main() {
 
     double maxf;
     nlopt::result result = opt.optimize(x, maxf);
-    for(int i=0; i<16; i++){
-        std::cout<<x[i]<<std::endl;
-    }
-    std::cout<<"valeur à la fin : "<<maxf<<std::endl;
 
-    return 0;
+    return x;
 }
 
     
diff --git a/src/option.cpp b/src/option.cpp
index 0e5a56c..b8c2799 100644
--- a/src/option.cpp
+++ b/src/option.cpp
@@ -1,62 +1,13 @@
-#include <algorithm>
+#include "option.hpp"
 #include <iostream>
-#include "rqmc.hpp"
-#include "p_adic.o"
-
-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, typename LDS>
-struct quasi_option : public generator<typename Fct::result_type>
-{
-    quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), G(d) {};
-        
-    typename Fct::result_type operator()() {
-        double sum =0;
-        for(int i=0; i<n; i++){
-            sum += payoff(G());
-            
-        }
-        return sum/n;
-    };
-    
-    private:
-        int n, d;
-        Fct payoff;
-        quasi_gaussian<LDS> G;
-    };
+double frac_part(double x){
+    return x - floor(x);
+}
 
 int main(){
     init_alea(1);
@@ -73,9 +24,9 @@ int main(){
     }
     
     std::vector<double> result2(3);
-    result = monte_carlo(100, quasi_option<asian_option, halton> (N, d, A));
+    result2 = monte_carlo(100, quasi_option<asian_option, halton> (N, d, A));
     for(int i =0; i<3; i++){
-        std::cout<<result[i]<<std::endl;
+        std::cout<<result2[i]<<std::endl;
     }
 
 
diff --git a/src/stratified_sampling.cpp b/src/stratified_sampling.cpp
index 156751f..1c8e670 100644
--- a/src/stratified_sampling.cpp
+++ b/src/stratified_sampling.cpp
@@ -1,5 +1,6 @@
 #include "stratified_sampling.hpp"
 
+
 std::pair<double, double> mean_var( std::vector<double> r){
     std::pair<double, double> p;
     for(auto &x: r){
@@ -12,3 +13,4 @@ std::pair<double, double> mean_var( std::vector<double> r){
     return p;
 }
 
+
diff --git a/src/stratified_sampling.hpp b/src/stratified_sampling.hpp
index cebee8d..7573ef2 100644
--- a/src/stratified_sampling.hpp
+++ b/src/stratified_sampling.hpp
@@ -5,6 +5,7 @@
 #include "rtnorm.hpp"
 #include "var_alea.hpp"
 #include <gsl/gsl_cdf.h>
+#include "option.hpp"
 
 using namespace std;
 
@@ -45,14 +46,14 @@ struct multi_gaussian_truncated : public var_alea<std::vector<double> >
         double v = V();
         double Z = gsl_cdf_gaussian_Pinv(v,1);
         std::vector<double> Y(d);
-        double scal = 0;
         for(int i=0; i<d; i++){
             Y[i] = G();
         }
-        std::vector<double> X(d);
+        double scal = 0;
         for(int j=0; j<d; j++){
             scal += Y[j]*u[j];
         }
+        std::vector<double> X(d);
         for(int i=0; i<d; i++){
             X[i] = u[i]*Z + Y[i] - u[i]*scal;
         }
@@ -78,7 +79,7 @@ struct stratified_sampling {
     void print_sigma() const;
     pair<double,double> estimator() const;
 private:
-		void update(int N);
+    void update(int N);
     vector<double> p;
     vector<L> X;
     vector<int> M;
@@ -187,3 +188,29 @@ pair<double,double> stratified_sampling<L>::estimator() const {
 	}
 	return {est_mean, est_std};
 }
+
+struct f_mu : public std::unary_function<std::vector<double>, double> 
+{
+    f_mu(std::vector<double> mu, asian_option A) : mu(mu), A(A){
+        double norm_mu = 0;
+        for(int i=0; i<16; i++) {
+            norm_mu += mu[i]*mu[i];
+        }
+    };
+    
+    double operator()(std::vector<double> X) {
+        std::vector<double> Y(X.size());
+        double scal = 0;
+        for (unsigned int i=0; i<X.size(); i++){
+            Y[i] = X[i] + mu[i];
+            scal+=X[i]*mu[i];
+        }
+        return A(Y)*exp(-scal-0.5*norm_mu);
+    };
+        
+    private :
+        std::vector<double> mu;
+        asian_option A;
+        double norm_mu;
+    };
+