//  基于方差Gamma模型的随机抽样

//  by fantuanxiaot

#include

#include

#include

#include

#include

#include

using namespace std;

#define pi 3.14159265359

//  首先定义几个重要的分布函数和分位数的求解

//

//

//  正态分布的分位数的求解

//  计算分位数利用牛顿拉夫森发则

double erf(double a)

{

double a1=0.0705230784,

a2=0.0422820123,

a3=0.0092705272,

a4=0.0001520143,

a5=0.0002765672,

a6=0.0000430638;

a=1+a1*a+a2*a*a+a3*a*a*a+a4*a*a*a*a+a5*a*a*a*a*a+a6*a*a*a*a*a*a;

a=pow(a,16);

a=1-1/a;

return(a);

}

double f_norm(double x,double p)

{

double result;

if(x>=0)

{

result=0.5+0.5*erf(x/sqrt(2))-p;

} else if(x<0)

{

result=0.5-0.5*erf(fabs(x/sqrt(2)))-p;

}

return(result);

}

double g_norm(double x)

{

return(1.0/sqrt(2*pi)*exp(-pow(x,2)/2.0));

}

double norminv(const double p)

{

double x,x0=1;

double eps=1e-3;

while(true)

{

x=x0-f_norm(x0,p)/g_norm(x0);

if(fabs(x-x0)<=eps)

{

break;

}

x0=x;

}

return x;

}

//  正态分布的分位数的求解

//

//

//  Gamma分布的分位数的求解

//  计算gamma函数

double gammafun(const double xx)

{

int j;

double x,y,tmp,ser;

const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,

-1.231739572450156,0.1208650973866179e-2,-0.5395239384953e-5};

y=x=xx;

tmp=x+5.5;

tmp=tmp-(x+0.5)*log(tmp);

ser=1.0000000000190015;

for(j=0;j<=5;j++)

{

ser=ser+cof[j]/(++y);

}

return(exp(-tmp+log(2.5066282746310005*ser/x)));

}

//  计算gamma的密度和分布函数

double g_gamma(double x,double alpha,double beta)

{

return(exp(-beta*x)*pow(beta,alpha)*pow(x,alpha-1)/gammafun(alpha));

}

double f_gamma(double x,double alpha,double beta,const double p)

{

//  采用梯形法计算函数

double sum=0.0;

double gaps=(x-0.0)/double(1000);  //每个间隔的长度

for (int i=0;i<1000;i++)

{

sum+=(gaps/2.0)*(g_gamma(0+i*gaps,alpha,beta)+g_gamma(0+(i+1)*gaps,alpha,beta));

}

return sum-p;

}

//  计算分位数函数

double gammainv(double p,const double alpha,const double beta)

{

double eps=1e-3;

double x2;

double x0=2;

while(true)

{

x2=x0-f_gamma(x0,alpha,beta,p)/g_gamma(x0,alpha,beta);

if(fabs(x2-x0)<=eps)

{

break;

}

x0=x2;

}

return x2;

}

//  Gamma分布的分位数的求解

//

//

//  构建一个方差Gamma随机抽样的类

class VarianceGamma

{

private:

//  以下是求解的参数

double Time;

double Variance;

double Skewness;

double Riskless;

double Dividend;

double S0;

double Sigma;

public:

VarianceGamma(double t,double v,double s,double r,double d,double s0,double si):

Time(t),Variance(v),Skewness(s),Riskless(r),Dividend(d),S0(s0),Sigma(si){}

//  以下是一个随机数生成的函数

double VGammaPrice(const double randomNumber)

{

double Omiga=Time/Variance*log(1-Skewness*Variance-Sigma*Sigma*Variance/2);

double alpha=Time/Variance;

double beta=1/Variance;

double a=gammainv(randomNumber,alpha,beta);

double b=a*Skewness+Sigma*sqrt(a)*norminv(randomNumber);

double Sample=S0*exp((Riskless-Dividend)*Time+Omiga+b);

return Sample;

}

};

//

//  最终是主函数

int main()

{

VarianceGamma VG(0.5,0.5,0.1,0,0,100,0.2);

vector Price;

//  输入100个价格和100个随机数

default_random_engine generator;

uniform_real_distribution unif_dis(0,1);

double randomNumber[100];

for(int i=0;i<=99;i++)

{

Price.push_back(VG);

randomNumber[i]=unif_dis(generator);

}

//  输出方差Gamma分布的抽样

for(int i=0;i<=99;i++)

{

cout<

if (i%5==4)

{

cout<

}

}

return 0;

}