金融计算与建模(下)(清华大学,朱世武)

第13章随机模拟基础清华大学经管学院朱世武[e-mail]样本数据：论坛：分布的模拟实现中心极限定理（对于和）：假设X1,…,Xn独立同分布，均值和标准差分别为μ和σ。令，如果n足够大（比如，n大于30），则Sn近似地服从正态分布，均值为nμ，方差为。中心极限定理（对于均值）：假设X1,…,Xn独立同分布，均值和标准差分别为μ和σ。令，如果n足够大（比如，n大于30），则近似地服从正态分布，均值为μ，方差为。中心极限定理表明，满足一定条件的独立同分布随机变量的和或均值将服从正态分布，而与它们的分布情况无关。下面分别给出中心极限定理的模拟实现结果。niinXS1nnXXnii1Xn/随机变量和的分布模拟例13.1假设掷骰子n次，分别用X1,…,Xn表示每次得到的数值。令，则Sn近似于正态分布。分别用SAS系统实现n=2和n=3的模拟分布图。n=2的模拟分布图。dataa;dox1=1to6;dox2=1to6;output;end;end;/*模拟掷骰子两次，生成36行数据*/dataa;seta;x=sum(x1,x2);procunivariatedata=anoprint;varx;histogram/normal(mu=estsigma=est);run;2468101205101520253035PercentxN=234.567.5910.51213.51516.5180510152025PercentxN=3例13.2二项分布是二点分布的和。所以，直接作二项分布的分布图也就是也就是对二点分布和的模拟。设二项分布的p=0.8，分别作N=10,20,25,50,100,200,1000时的分布图。直接求概率后画图。%macroa(n);/*建立一个以n为参数的宏*/datarv;dom=1to&n;ifm=0thenprobb=probbnml(0.8,&n,0);/*probbnml函数返回以0.8、n、0为参数的二项分布概率*/elseprobb=probbnml(0.8,&n,m)-probbnml(0.8,&n,(m-1));/*注意到probbnml函数返回累积概率，所以这里要做减法*/output;end;symbol1i=needlewidth=6c=blueh=1cells;/*规定两点之间的插值方法、线条宽度、数据点符号的高度及颜色，注意这里使用了symbol1*/procgplotdata=rv;plotprobb*m=1;%menda;%a(10);%a(15);%a(20);%a(25);%a(50);%a(100);%a(200);%a(1000);run;probb0.000.020.040.060.080.100.120.140.160.180.200.220.240.26m123456789101112131415probb0.000.010.020.030.040.050.060.070.080.090.100.110.120.130.14m01020304050probb0.000.010.020.030.040.050.060.070.08m0102030405060708090100110120130140150160170180190200probb0.0000.0020.0040.0060.0080.0100.0120.0140.0160.0180.0200.0220.0240.0260.0280.0300.032m01002003004005006007008009001000模拟二项分布随机数的分布图。（设二项分布的p=0.8，分布模拟N=10,20,25,50,100,200,1000时随机数的分布图。）symbol;goptionsftext=ctext=htext=;%macroa(n);datarv;retain_seed_0;do_i_=1to&n;binom1=ranbin(_seed_,&n,0.8);/*ranbin函数返回以seed、n、0.8为参数的二项分布随机变量*/output;end;drop_seed__i_;run;procunivariatedata=rvnoprint;varbinom1;histogram/normal(mu=estsigma=est);insetnormal;run;%menda;%a(10);%a(15);%a(20);%a(25);%a(50);%a(100);%a(200);%a(1000);run;789100510152025303540Percentbinom1101112131405101520253035Percentbinom133.7536.2538.7541.2543.7546.2505101520253035Percentbinom176276877478078679279880481081682282883484002.55.07.510.012.515.017.520.0Percentbinom1前面所有模拟结果表明，随机n的增大，和的分布越来越接近正态分布。事实上，直观分析就可以得出离散的随机变量都有这样的特性，如上面掷骰子、二项分布等。当然中心极限定理保证连续分布也有这样的性质。随机变量均值的分布模拟例13.4对于=6的指数分布，分别模拟出N=1,10,20,50,100,200,1000.时，均值的分布。1NXXXN统计抽样中的分布模拟例13.5设总体为=6指数分布的100个随机数。设样本容量为80,进行n次抽样，分别用表示第1次,……,第n次抽得的随机数。，模拟和的分布。对于求和，随机抽取的次数分别为n=5,10,20,50,100,200,1000.1,,nXXnnXXS1nSXnnXXS1symbol;goptionsftext=ctext=htext=;optionsnodatenonotesnosource;datarv;/*创建由100个随机数构成的总体*/retain_seed_0;do_i_=1to100;exp=ranexp(_seed_)/6;output;end;drop_seed__i_;datarv80;/*从总体中选择80个作为样本*/delete;%macrob(y);%doi=1%to%eval(&y);dataa;setrv;obs=_n_;procsql;createviewtmpas/*建立名为tmp的视图（view）*/select*,ranuni(0)as_ran_froma/*生成变量_ran_，服从[0，1]上的均匀分布*/orderbycalculated_ran_;/*按_ran_排序*/quit;datarandom;settmp(obs=80);/*取前80个*/drop_ran_;renameexp=exp&i;datarv80;mergerv80random;%end;datarv80(keep=s&y);setrv80;s&y=sum(ofexp1-exp&y);procunivariatedata=rv80noprint;vars&y;histogram/normal(mu=estsigma=estnoprint);insetnormal;run;%mendb;%b(5);%b(10);%b(20);%b(50);%b(100);%b(200);%b(1000);run;0.150.450.751.051.351.651.952.2505101520253035Percents51.82.433.64.24.85.4051015202530Percents201213.51516.51819.52122.505101520253035Percents100151.5154.5157.5160.5163.5166.5169.5172.50510152025Percents1000样本均值的分布模拟。准备数据集：产生100个指数分布的随机数。每次随机抽取的次数分别为n=5,10,20,50和80,抽取180次样进行模拟。产生100个指数分布的随机数。symbol;goptionsftext=ctext=htext=;optionsnodatenonotesnosource;datarv;retain_seed_0;do_i_=1to100;exp=ranexp(_seed_)/6;output;end;drop_seed__i_;run;样本容量分别为5,10,20,50和80,抽180次样。的分布模拟。symbol;goptionsftext=ctext=htext=;optionsnodatenonotesnosource;%macrob(y,z);datarv&y;delete;%doi=1%to%eval(&z);dataa;setrv;obs=_n_;procsql;/*sql过程创建视图tmp*/createviewtmpasselect*,ranuni(0)as_ran_fromaorderbycalculated_ran_;quit;datarandom;settmp(obs=&y);/*取前&y个观测值*/drop_ran_;procmeansdata=randomnoprint;varexp;outputout=bmean=exp_m;datarv&y;setrv&yb;/*通过这个语句，使每次宏得到的exp_m全部汇总到rv&y这个文件中了*/%end;procunivariatedata=rv&ynoprint;varexp_m;histogram/normal(mu=estsigma=estnoprint);insetnormal;run;%mendb;%b(5,180);%b(10,180);%b(20,180);%b(50,180);%b(80,180);run;0.020.060.10.140.180.220.260.30.340.380510152025Percentexp_m0.1050.1150.1250.1350.1450.1550.1650.1750.185051015202530Percentexp_m0.1250.130.1350.140.1450.150.1550.160.1650510152025Percentexp_m例中，采用了不同于前面的横向加总方法（用到了PROCMEANS纵向计算均值），并在分别计算EXP_M的基础上用SET语句汇总了每次循环所得的180个数值，最后进行模拟收益模型模拟随机游动模型例13.7价格随机游动模型,的模拟实现。假设模型为，初值P1=8.6。1tttPP~(0,1)tN1tttPP~(0,1)tNdataa;mu=0;p1=8.6;sigma=1;dotime=-50to1000;e=rannor(32585);/*注意到，此处rannor随机函数的种子为一个特定的数值，这是每次运行，生成的随机数列是相同的，而如果写成rannor(0)，那么每次重新运行，生成的随机数列也不尽相同*/p=mu+p1+sigma*e;iftime0thenoutput;/*实际上是去掉最开始的51个数据*/p1=p;/*保证了循环的进行*/end;run;随机序列时序图：procgplotdata=a;symbol1v=pointi=joinc=blue;/*规定符号样式、插值方法、线条颜色*/symbol2v=nonei=r;/*i=r规定生成一条一次的回归曲线（直线）*/plotp*time=1p*time=2/overlay;/*在一张图中产生两条曲线，分别使用样式1和样式2*/run;p10203040time01002003004005006007008009001