diff options
Diffstat (limited to 'option.cpp')
| -rw-r--r-- | option.cpp | 105 |
1 files changed, 48 insertions, 57 deletions
@@ -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; |