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

}

總結

以上是生活随笔為你收集整理的matlab dividend,[原创]基于(Matlab/R/C++)的方差Gamma模型(Hull期权期货)随机抽样[by fantuanxiaot]...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。