张晓峒老师的蒙特卡罗模拟.doc

专题5 蒙特卡罗模拟的有关问题大家知道，只有当经典回归模型满足所有的假定条件时，参数的估计量才具有最佳线性无偏特性，即有限样本特性，同时也具有渐近特性。当假定条件不成立时（比如存在异方差、自相关等），所采用的广义最小二乘法，以及对联立方程模型的估计，动态分布滞后模型的估计，向量自回归模型的估计所得参数的估计量只具有渐近特性。也就是说，只有当样本容量相当大时，渐近特性才起作用。而当样本容量不是很大，甚至很小时，仍然不知道估计量的有限样本分布特征。另外通过对非平稳过程的研究知单位根检验式和非平稳变量之间回归参数和t统计量不服从正态分布。他们都是渐近地服从Wiener过程函数的分布。参数估计量和统计量的有限样本特性不能用解析的方法求解。对于上述两种情形，若要研究这些估计量和统计量的有限样本分布特征，通常采用两种方法。一种为数值计算法。也称为有限样本近似法（finite-sample approximation）。这种方法要用到许多数学知识，专业性很强，使没有受过专门训练的人员运用此方法受到限制。（2）蒙特卡罗模拟方法。又称随机模拟法。Boot strap 1蒙特卡罗（Monte Carlo）模拟和自举（Boost trap）发展过程这是一种通过设定随机过程（数据生成系统），反复生成时间序列，并计算参数估计量和统计量，进而研究其分布特征的方法。蒙特卡罗在欧洲的摩那哥，以著名赌城而得名。据说这个术语是Metropolis 在1949年提出的。若再晚些时候，蒙特卡罗模拟也许就称作Las Vegas（在美国的Nevada州，著名赌城）模拟方法了。自举模拟与蒙特卡罗模拟既有联系，又不相同。自举（Boost trap，亦称靴襻）这个名词是Efron在1979年提出的。“自举”一词来源于儿童故事。指一个人落水时，试图用自提鞋扣儿的方法自救。20世纪80，90年代发展很快。自举，即采用从总体中反复抽取样本的方法计算参数估计量的值，置信区间或相应统计量的值并估计这些量的分布。这里介绍的远不是自举模拟的全貌，而是参数估计方面的应用。因为这些方法的实现是以高容量和高速度的计算机为前提条件，所以只是在近年才得到广泛推广。2蒙特卡罗模拟和自举模拟原理进行蒙特卡罗模拟和自举模拟首先要设定数据生成系统。而设定数据生成系统的关键是要产生大量的随机数。例如模拟样本为100的随机趋势过程的DF统计量的分布，若试验1万次，则需要生成200万个随机数。计算机所生成的随机数并不是“纯随机数”，而是具有某种相同统计性质的随机数。计量经济学中蒙特卡罗模拟和自举模拟所用到的随机数一般是服从N(0,1)分布的随机数。计算机生成的随机数称作“伪随机数”（pseudo-random number）。生成的随机数的程序称作“伪随机数生成系统”。实际上计算机不可能生成纯随机数。在进行蒙特卡罗模拟时一般要给定多种条件。例如样本容量要选择50，100，200等多种。有时模型形式也要选择多种。从而研究参数估计量和统计量在各种条件下的分布特征。当只需要这几个特定条件下的模拟结果时，把结果纪录下来就可以了。当需要很多条件下的模拟结果时，一般采用估计响应面函数（response surface function）的方法研究之。例如Dicky-Fuller的DF检验表中只给出了样本容量为25，50，100，250，500几个点的DF分布特征。显然对25至500间每个样本容量都进行DF分布模拟是不实际的，也是无必要的。可以把上述几个条件下得到的DF分布百分位数看作样本点，然后采用回归的方法从而得到每个样本容量所对应的DF分布百分位数。这条回归直线称为响应面函数。麦金农的协整检验临界值表就是用这种方法得到的。一个简单的估计回归参数估计量分布的蒙特卡罗模拟流程图见图1。生成xt, yt设定xt, yt 估计和t()统计量设定循环次数N分析和t()的分布 达到N 未达到N图1 蒙特卡罗模拟过程示意图自举方法的原理是从独立同分布（IID）总体X中确定T个随机变量x1, x2, , xT。则第一个自举样本是X1* = x11, x12, , x1T现在随机得到N个自举样本，X2* = x11, x12, , x1TX3* = x11, x12, , x1TXN* = xN1, xN2, , xNT假设关心的是统计量(X)，那么用N个自举样本可以得到一个容量为N的(X)的估计值序列，(X1*), (X2*), , (XN*)通过这个序列，可以研究(X)的分布特征，(X)的特征数，百分位数，(X)的平均数与真值q 的差以及用(X)的第a/2、（1-a/2）百分位数构造q 的（1-a）的置信区间。一个简单的分析t()分布特征的自举模拟流程图见图2。设定xt,生成xt, 估计和 t()统计量 设定循 环次数 生成 yt 分析和 t()的分布 达到 设定b0, b1 未达到 图2 自举模拟过程示意图3计算机高级语言（Mathematica和EViews介绍） 蒙特卡罗模拟和自举模拟的实现要通过计算机编程来实现。常用的软件有Mathematica，Gauss，Ox，EViews等。其原理基本一样。下面主要介绍EViews和Mathematica。Mathematica由Wolfram Research公司1991年推出。是一种计算机高级语言。具有计算与画图等多种功能。若干例子见图。 图3 随机游走序列 图4 带趋势项的随机游走序列 图5 三维图圆环 图6 空间曲面 图7 投币1000次的概率值模拟 图8 生长曲线图9 二元正态分布 图10 蒲丰问题4蒙特卡罗模拟框图与Mathematica、EViews程序。（1）两个I(1)变量相关系数分布的蒙特卡罗模拟。估计相关 系数r 分析r的 分布 生成 xt, ytI(1) 设定 xt, yt I(1) 设定循环 次数N 达到N 未达到N图11 蒙特卡罗模拟过程示意图Mathematica程序如下：corre2t_,f_:= Modulex,y,xx,yy,Exx,Eyy,Sxxyy,Sxx,Syy,rr, Table x=TableRandomNormalDistribution0,1,t; y=TableRandomNormalDistribution0,1,t; xx=FoldListPlus,0,x;xx=Restxx; yy=FoldListPlus,0,y;yy=Restyy; Exx=ApplyPlus,xx/t; Eyy=ApplyPlus,yy/t; Sxxyy=(xx-Exx).(yy-Eyy); Sxx=Sqrt(xx-Exx).(xx-Exx); Syy=Sqrt(yy-Eyy).(yy-Eyy); rr=Sxxyy/(Sxx Syy), f r2=corre2100,10000;histg4r2,0,1,0.1图12 两个非相关I(1) 序列的相关系数的分布EViews程序如下：workfile corr u 1 500series resultfor !i=1 to 500smpl 1 100series x=nrndseries y=nrndseries xxseries yyscalar sum1=0scalar sum2=0for !counter=1 to 100sum1=sum1+x(!counter)sum2=sum2+y(!counter) xx(!counter)=sum1yy(!counter)=sum2nextscalar r=cor(xx,yy)result(!i)=rnextresult.hist定义一个非时间序列（u）工作文件，corr，容量为500。定义一个空序列result，用来存储相关系数的计算结果。!i为控制变量，通过一个for循环语句使计算进行500次。把样本范围设置成100。生成两个互不相关的白噪声序列x、y，样本容量100。定义两个空的序列xx和yy，样本容量也是100。定义两个标量sum1和sum2，初始值为0。!counter为控制变量，在这个for循环中，分别对序列x和y进行一次累加生成两个一阶单整的序列，将结果分别放到序列xx和yy中。累加一次。计算序列xx和yy的相关系数，并将结果放到标量r中。将相关系数计算结果放到序列result中，在这个for循环中，这个操作要进行500次。显示序列result的直方图以及有关统计量。图13 两个非相关I(1) 序列的相关系数的分布（2） t()分布的蒙特卡罗模拟。数据生成过程如下， yt = yt-1 + ut , ut IID(0, 1) 估计的方程式如下： yt = m +b yt-1 + ut ,检验统计量 t()=图14 t()统计量分布的蒙特卡罗模拟（T =50，模拟1万次）（3）DW统计量分布的蒙特卡罗模拟生成T=50的相互独立的IN(0,1)序列ut 和vt用ut 和vt分别生成两个相互独立的I(1)序列yt = yt-1 + ut , y0 = 0, xt = xt-1 + vt , x0 = 0,估计模型yt = b0 + b1xt + wt 并计算残差用残差计算DW统计量的值存储2000个DW值画DW频数分布直方图。记录T=50条件下DW分布的均值、标准差和第90、95、99百分位数。分别估计DW均值、标准差和第90、95、99百分位数值对(1/T )的响应面函数图15（2）DW统计量分布蒙特卡罗模拟的Mathematica程序DWvaluet_,f_:= Modulex1,y1,xx,yy,x0,x,A,B,para,u,sig,u1,Su1,DW, Tablex1=TableRandomNormalDistribution0,1,t; xx=FoldListPlus,0,x1;xx=Restxx; y1=TableRandomNormalDistribution0,1,t; yy=FoldListPlus,0,y1;yy=Restyy; (* to estimate regression parameters *) x0=Table1,t; x=x0,xx; A=Transposex; B=Inversex.A; para=B.x.yy; (* to calculate the residuals and DW value *) u=yy-A.para; sig=u.u; u1=Tableui-ui-1,i,2,t; Su1=u1.u1; DW=Su1/sig,f w5=DWvalue100,10000;histg4w5,0.,2.8,0.1下面以论文小样本DW统计量的分布特征（南开经济研究1999第6期）为例介绍蒙特卡罗模拟流程。小样本I(1)变量的DW分布的蒙特卡罗模拟框图如上。以样本容量T = 10, 20, 30, 40 50为条件，生成相互独立的两个I(1)序列，xt, yt。每生成一对xt, yt序列，用yt对xt,回归，计算DW值。在每一个固定的样本容量条件下，模拟2000次。画DW分布直方图，计算DW分布的均值、标准差和第90、95、99百分位数值。T = 50，模拟2000次，计算机输出结果如下：Percentiles and Descriptive Statisticshistogram:DW频率频数1 percentile = 0.045858310.007840.58847.5 percentile = 0.090087820.08782.08166.10 percentile = 0.12374930.1682.25180.90 percentile = 0.59127440.2482.18174.95 percentile = 0.70827350.3281.8144.99 percentile = 0.98293260.4081.3104.Min value = 0.015839270.4880.81365.Max value = 1.275480.5680.648.Mean = 0.33227490.6480.3528.Standard Dev. = 0.196302100.7280.18815.Skewness = 1.16207110.8080.18815.Kurtosis = 4.67025120.8880.03753.Jarque-Bera = 341.307130.9680.18.141.050.01251.151.130.01251.161.210.01251.171.2900图16 T = 50条件下DW的分布5自举模拟过程框图（模型误差项非自相关、回归参数0.8的的分布）设定样本容量T=100 T= 40生成序列xt I(1) + t设定 b0=0.5, b1=0.8生成ut设定f1 =0生成yt = 0.5 + 0.8 xt + utOLS估计 否循环次数f达5000 ? 是估计分布的特征数图17模拟结果分析。模型（3）是在xt, yt非平稳条件下的OLS回归结果。具有超一致性。见表1。 表1 误差项非自相关条件