相似文件
换一批
经管之家送您一份
应届毕业生专属福利!
求职就业群
感谢您参与论坛问题回答
经管之家送您两个论坛币!
+2 论坛币
方差Gamma模型来自Hull的期权期货,用于捕捉股价的跳跃
现在对方差Gamma模型进行样本抽样
资料如下:
一,基于Matlab编程:生成1000个样本并和几何布朗运动的概率分布进行比较
程序如下:
本帖隐藏的内容- %% 基于方差Gamma模型的样本抽取
- %% by fantuanxiaot
- function Sample=Variance_Gamma(Time,Variance,Skewness,Riskless,Dividend,S0,Sigma)
- %% 举例
- % Variance_Gamma
- % Sample=Variance_Gamma
- % Variance_Gamma(1,0.4,-0.3,0,0,100,0.3)
- % Sample=Variance_Gamma(1,0.4,-0.3,0,0,100,0.3)
- %% Begin
- % Time时间
- % Variance是Gamma过程的方差率
- % Skewness定义了偏态
- % Riskness无风险利率
- % Dividend定义红利率
- % S0定义了基期价格
- % Sigma定义波动率
- if nargin==0
- Time=0.5;
- Variance=0.5;
- Skewness=0.1;
- Riskless=0;
- Dividend=0;
- S0=100;
- Sigma=0.2;
- end
- % 首先构建10000个0~1的随机数
- RandSample=rand(1,1000);
- Omiga=Time/Variance*log(1-Skewness*Variance-Sigma^2*Variance/2);
- % 抽取样本
- G=gaminv(RandSample,Time/Variance,Variance);
- B=G*Skewness+Sigma*sqrt(G).*norminv(RandSample,0,1);
- Sample=S0*exp((Riskless-Dividend)*Time+Omiga+B);
- %% 构建对比的几何布朗运动
- % 首先构建10000个正态分布的随机数
- ST=S0*exp((Riskless-Dividend)*Time+Sigma*normrnd(0,Time,1,1000));
- if nargout==0
- % 构建密度函数进行对比
- [f1,x1]=ksdensity(Sample);
- figure('color','w')
- plot(x1,f1,'bo-','markersize',10,'markerfacecolor','r')
- hold on
- [f2,x2]=ksdensity(ST);
- plot(x2,f2,'go-','markersize',10,'markerfacecolor','m')
- legend('方差Gamma','几何布朗')
- set(gca,'fontname','华文楷体','fontsize',14)
- title('方差Gamma分布的样本抽样','fontname','华文楷体','fontsize',18)
- end
复制代码
效果可见:
二,基于R语言编程:
本帖隐藏的内容- ## ########################################################
- ## ########################################################
- ## 基于方差Gamma模型的随机抽样
- ## by fantuanxiaot
- Variance_Gamma<-function(Time,Variance,Skewness,
- Riskless,Dividend,S0,Sigma)
- {
- # 举例子:Sample<-Variance_Gamma(0.5,0.5,0.1,0,0,100,0.2)
- Omiga<-Time/Variance*log(1-Skewness*Variance-
- Sigma^2*Variance/2);
- # 首先构建10000个0~1的随机数
- RandSample<-runif(1000,0,1)
- # 抽取样本
- Variance<-1/Variance
- G<-qgamma(RandSample,Time/Variance,Variance)
- B<-G*Skewness+Sigma*sqrt(G)*qnorm(RandSample,0,1)
- Sample<-S0*exp((Riskless-Dividend)*Time+Omiga+B)
- # 构建对比的几何布朗运动
- # 首先构建10000个正态分布的随机数
- ST<-S0*exp((Riskless-Dividend)*Time+
- Sigma*rnorm(1000,0,Time))
- # 作图进行几何布朗运动和方差Gamma模型的随机抽样比较
- par(family='serif')
- col=rainbow(2)
- Data.Sample<-density(Sample)
- Data.ST<-density(ST)
- plot(Data.Sample$x,Data.Sample$y,,type='b',pch=21,
- lwd=1,col='blue',bg=col[1],
- main='方差Gamma模型样本抽样',
- xlim=c(0,200))
- lines(Data.ST$x,Data.ST$y,type='b',pch=21,
- lwd=1,col='green',bg=col[2])
- legend(20,100,lty=c(1,1),
- legend=c("方差Gamma","几何布朗"),
- col=col,hori=FALSE,bty="o")
- Data<-list()
- Data$VaGamma<-Sample
- Data$GeoBrown<-ST
- return(Data)
- }
- ## 样本的抽取
- Variance_Gamma(0.5,0.5,0.1,0,0,100,0.2)
- Variance_Gamma(0.5,0.5,-0.2,0,0,100,0.2)
- ## ########################################################
- ## ########################################################
复制代码
效果如下:
三,基于C++编程:
首先我们需要计算正态分布和Gamma分布的分位数计算函数
其次构建一个Gamma方差的类
并从中抽取100个样本
代码如下:
本帖隐藏的内容- // 基于方差Gamma模型的随机抽样
- // by fantuanxiaot
- #include <iostream>
- #include <cmath>
- #include <random>
- #include <iomanip>
- #include <random>
- #include <vector>
- 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<VarianceGamma> Price;
- // 输入100个价格和100个随机数
- default_random_engine generator;
- uniform_real_distribution<double> 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<<Price[i].VGammaPrice(randomNumber[i])<<" ";
- if (i%5==4)
- {
- cout<<endl;
- }
- }
- return 0;
- }
复制代码
最终效果:抽取100个随机样本
扫码加我 拉你入群
请注明:姓名-公司-职位
以便审核进群资格,未注明则拒绝
|