23(蒙特卡罗模拟)

1 专题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分布模拟是不实际的，也是无必要的。 2

可以把上述几个条件下得到的DF分布百分位数看作样本点，然后采用回归的方法从而得到每个样本容量所对应的DF分布百分位数。这条回归直线称为响应面函数。麦金农的协整检验临界值表就是用这种方法得到的。 一个简单的估计回归参数估计量分布的蒙特卡罗模拟流程图见图1。

达到N

未达到N 图1 蒙特卡罗模拟过程示意图

自举方法的原理是从独立同分布（IID）总体X中确定T个随机变量{x1, x2, …, xT}。则第一个自举样本是 X1* = {x11, x12, …, x1T} 现在随机得到N个自举样本， X2* = {x11, x12, …, x1T} X3* = {x11, x12, …, x1T} ．．． XN* = {xN1, xN2, …, xNT}

假设关心的是统计量ˆ(X)，那么用N个自举样本可以得到一个容量为N的ˆ(X)的估计值序列， {ˆ(X1*), ˆ(X2*), …, ˆ(XN*)}

通过这个序列，可以研究ˆ(X)的分布特征，ˆ(X)的特征数，百分位数，ˆ(X)的平均数与真值 的差以及用ˆ(X)的第/2、（1-/2）百分位数构造 的（1-）的置信区间。 一个简单的分析t(1ˆ)分布特征的自举模拟流程图见图2。

达到 未达到 图2 自举模拟过程示意图 3．计算机高级语言（Mathematica和EViews介绍） 蒙特卡罗模拟和自举模拟的实现要通过计算机编程来实现。常用的软件有Mathematica，Gauss，Ox，EViews等。其原理基本一样。 下面主要介绍EViews和Mathematica。Mathematica由Wolfram Research公司1991年推出。是一种计算机高级语言。具有计算与画图等多种功能。若干例子见图。

生成xt, yt 估计1ˆ和 t(1ˆ)统计量 分析1ˆ和 t(1ˆ)的分布 设定循环 次数N 设定xt,生成xt, 估计1ˆ和 t(1ˆ)统计量 分析1ˆ和 t(1ˆ)的分布 设定循 环次数 生成 yt 设定0, 1 设定xt, yt 3

050100150200-10-505 0501001502000

10

203040

图3 随机游走序列 图4 带趋势项的随机游走序列 -4-2024-1-0.500.51-4-2024 -2

02-1-0.5

0

0.51

-202

图5 三维图圆环 图6 空间曲面

0200400600800100000.20.40.60.81

图7 投币1000次的概率值模拟 图8 生长曲线 -2-2-101200.250.50.751-1012

图9 二元正态分布 图10 蒲丰问题 4

4．蒙特卡罗模拟框图与Mathematica、EViews程序。 （1）两个I(1)变量相关系数分布的蒙特卡罗模拟。

达到N

未达到N 图11 蒙特卡罗模拟过程示意图 Mathematica程序如下： corre2[t_,f_]:= Module[{x,y,xx,yy,Exx,Eyy,Sxxyy,Sxx,Syy,rr}, Table[ x=Table[Random[NormalDistribution[0,1]],{t}]; y=Table[Random[NormalDistribution[0,1]],{t}]; xx=FoldList[Plus,0,x];xx=Rest[xx]; yy=FoldList[Plus,0,y];yy=Rest[yy]; Exx=Apply[Plus,xx]/t; Eyy=Apply[Plus,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=corre2[100,10000]; histg4[r2,0,1,0.1]

-1-0.500.510.10.20.30.40.50.6

图12 两个非相关I(1) 序列的相关系数的分布

生成 xt, ytI(1) 估计相关 系数r 分析r的 分布 设定循环 次数N 设定 xt, yt I(1) 5 EViews程序如下： workfile corr u 1 500 series result for !i=1 to 500 smpl 1 100 series x=nrnd series y=nrnd series xx series yy scalar sum1=0 scalar sum2=0 for !counter=1 to 100 sum1=sum1+x(!counter) sum2=sum2+y(!counter) xx(!counter)=sum1 yy(!counter)=sum2 next scalar r=@cor(xx,yy) result(!i)=r next result.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 =  + yt-1 + ut ,

检验统计量 t(ˆ)=)ˆ(ˆs