// 基于方差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;
}