aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile2
-rw-r--r--option.cpp105
-rw-r--r--rqmc.cpp7
-rw-r--r--rqmc.hpp96
4 files changed, 147 insertions, 63 deletions
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;
+};
+
+
+
+
+
+
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2025-12-21 23:14:28 +0000