蒙特卡罗模拟与编程(张晓峒)

蒙特卡罗模拟与编程张晓峒南开大学数量经济研究所所长、博士生导师东北财经大学兼职教授中国数量经济学会常务理事、天津市数量经济学会理事长[e-mail]蒙特卡罗模拟与编程1．蒙特卡罗（MonteCarlo）模拟和自举（Bootstrap）的命名过程2．蒙特卡罗模拟和自举原理3．设计计算过程4．生成服从某种分布的随机数，生成各种类型的随机序列5．模拟模型回归系数有限样本估计量和检验统计量的分布特征6．估计响应面函数7．极大似然估计编程8．模拟中应注意的问题蒙特卡罗模拟与编程大家知道，只有当回归模型满足OLS法所有的假定条件时，参数的估计量才具有最佳线性无偏特性，同时也具有一致性。当假定条件不成立时（比如存在异方差、自相关等），所采用的广义最小二乘法，以及对联立方程模型的估计，动态分布滞后模型的估计，向量自回归模型的估计所得参数的估计量只具有渐近特性（渐近无偏性、一致性）。这意味着，只有当样本容量相当大时，渐近特性才起作用。而当样本容量不是很大，甚至很小时，仍然不知道估计量的有限样本分布特征。另外通过对非平稳过程的研究知道，单位根检验式和用非平稳变量建立的回归函数的参数和t统计量都不服从正态分布。他们都是渐近地服从Wiener过程的泛函。参数估计量和统计量的有限样本特性不能用解析的方法求解。蒙特卡罗模拟与编程对于上述两种情形，若要研究这些估计量和统计量的有限样本分布特征，通常采用两种方法。一种为数值计算法。也称为有限样本近似法（finite-sampleapproximation）。这种方法要用到许多数学知识，专业性很强，使没有受过专门训练的人员运用此方法受到限制。另一种为蒙特卡罗模拟方法。又称随机模拟法。本章主要介绍蒙特卡罗模拟与计算机编程。1．蒙特卡罗（MonteCarlo）模拟和自举（Bootstrap）的命名过程蒙特卡罗（MonteCarlo）模拟是一种通过设定随机过程（数据生成系统），反复生成随机序列，并计算参数估计量和统计量，进而研究其分布特征的方法。“蒙特卡罗模拟”这个术语是美国物理学家Metropolis在第2次世界大战时期执行曼哈顿计划(ManhattanProject)过程中提出的。作为地名，蒙特卡罗在欧洲的摩那哥（Monaco），以著名赌城而得名。若再晚些时候，蒙特卡罗模拟也许就称作LasVegas（在美国的Nevada州，著名赌城）模拟方法了。作为地名，蒙特卡罗在欧洲的摩纳哥（Monaco），以著名赌城而得名。蒙特卡罗自举与蒙特卡罗模拟既有联系，又不相同。自举（Bootstrap）这个名词是Efron在1979年提出的。“自举”一词来源于童话故事。指一个人落水时，试图用自提鞋扣儿的方法自救。自举（Bootstrap）有人翻译成“靴襻”不恰当。自举，即采用从样本中反复抽取子样本的方法计算参数估计量的值，置信区间或相应统计量的值并估计这些量的分布。这里介绍的远不是蒙特卡罗模拟的全部，而是参数估计方面的应用。因为这些方法的实现是以高容量和高速度的计算机为前提条件，所以只是在近年才得到广泛推广。20世纪80，90年代发展很快，现在已很普及。2．蒙特卡罗模拟和自举原理进行蒙特卡罗模拟和自举首先要设定数据生成系统。而设定数据生成系统的关键是要产生大量的随机数。例如模拟样本容量为100的一元线性回归模型中参数的分布，若试验1万次，则需要生成200万个随机数。计算机所生成的随机数并不是“真随机数”，而是具有某种相同统计性质的随机数。计量经济学中蒙特卡罗模拟和自举所用到的随机数一般是服从N(0,1)分布或均匀分布的随机数。计算机生成的随机数称作“伪随机数”（pseudo-randomnumber）（以下简称随机数）。生成的随机数的程序称作“伪随机数生成系统”。实际上计算机不可能生成真随机数。3．设计计算过程。蒙特卡罗模拟和自举的实现要通过计算机编程来实现。常用的高级编程语言软件有Mathematica，Gauss，Ox，EViews，S-Plus，Matlap，rats，stata，R等。这里主要介绍利用EViews做蒙特卡罗模拟与编程。编程也是一个很宽的概念，这里只介绍与蒙特卡罗模拟有关的编程。编程像一门艺术，需要经验和技巧。只有多思，熟记各种命令，才能编出最优的程序来。否则会浪费大量时间。首先需要把自己需要计算的问题设计出合理的计算流程框图，然后转化为正确的计算机语言。提出研究的问题设计计算流程框图转换为计算机语言4．生成服从某种分布的随机数，生成各种类型的随机序列生成服从某种分布的随机数、累计分布函数、概率密度函数、p分位数函数的表达语言见表1。表1：生成18种分布的随机数、累计分布函数、概率密度函数、p分位数函数的用语分布类型生成随机数累计分布函数概率密度函数p分位数函数Beta分布@rbeta(a,b)@cbeta(x,a,b)@dbeta(x,a,b)@qbeta(p,a,b)二项式分布@rbinom(n,p)@cbinom(x,n,p)@dbinom(x,n,p)@qbinom(s,n,p)卡方分布@rchisq(v)@cchisq(x,v)@dchisq(x,v)@qchisq(p,v)指数分布@rexp(m)@cexp(x,m)@dexp(x,m)@qexp(p,m)极值分布I@rextreme@cextreme(x)@dextreme(x)@qextreme(p)F分布@rfdist(v1,v2)@cfdist(x,v1,v2)@dfdist(x,v1,v2)@qfdist(p,v1,v2)Gamma分布@rgamma(b,r)@cgamma(x,b,r)@dgamma(x,b,r)@qgamma(p,b,r)广义误差分布@rged(r)@cged(x,r)@dged(x,r)@qged(p,r)拉普拉斯分布@rlaplace@claplace(x)@dlaplace(x)@qlaplace(p)logistic分布@rlogistic@clogistic(x)@dlogistic(x)@qlogistic(p)对数正态分布@rlognorm(m,s)@clognorm(x,m,s)@dlognorm(x,m,s)@qlognorm(p,m,s)负二项分布@rnegbin(n,p)@cnegbin(x,n,p)@dnegbin(x,n,p)@qnegbin(s,n,p)标准正态分布@rnorm,nrnd@cnorm(x)@dnorm(x)@qnorm(p)泊松分布@rpoisson(m)@cpoisson(x,m)@dpoisson(x,m)@qpoisson(p,m)帕雷托分布@rpareto(a,k)@cpareto(x,a,k)@dpareto(x,a,k)@qpareto(a,k)t分布@rtdist(v)@ctdist(x,v)@dtdist(x,v)@qtdist(p,v)均匀分布@runif(a,b)@cunif(x,a,b)@dunif(x,a,b)@qunif(p,a,b)威布尔分布@rweib(m,a)@cweib(x,m,a)@dweib(x,m,a)@qweiib(p,m,a)（1）生成服从某种分布的随机数序列【例】生成标准正态分布、指数分布、poisson分布、t分布的随机数序列EViews程序如下：'生成服从某种分布的随机数序列（T=1000），并存入random1文件。workfilerandom1u11000seriesY1=@rnormseriesY2=nrndseriesY3=@rexp(3)seriesY4=@rpoisson(3)seriesY5=@rtdist(3)-6-4-20242505007501000Y1标准正态分布随机数01020302505007501000Y3指数分布随机数02468102505007501000Y43个自由度的poisson分布随机数-15-10-505102505007501000Y5df=3的t分布随机数（2）生成任意（非标准）正态分布随机数序列用nrnd和@rnorm可以生成标准正态分布随机数，那么，任何参数的正态分布随机数都可以生成。用Z表示标准正态分布的随机数，用X表示任一参数的正态分布的伪随机数，那么，把X标准化的公式是ZsXXX)(N(0,1)反过来，若用Z表示X，则公式为X=Zs(X)+X按上公式，EViews可以生成任意正态分布随机数。【例】生成均值为50，标准差为0.5的正态分布随机数XN(50,0.52)。seriesZ=nrndseriesX=Z*0.5+5048.448.849.249.650.050.450.851.251.652.0501001502002503003504004505005506006507007508008509009501000X（3）在生成服从某种分布的随机数序列的基础上求分位数值，求频数分布表。计算分位数的EViews命令为：scalarq=@quantile(y,τ,s)，其中y表示求分位数所用的序列；τ表示要计算的分位数；s取1~6，依次表示6种计算分位数序数方法，上命令的含义是计算第τ分位数所对应的分位数值，存入标量q中。【例】以生成的Y1为例，求0.2分位数和中位数值。在空白处键入命令：scalarq1=@quantile(Y1,.2,1)scalarq2=@quantile(Y1,.5,1)q1和q2以标量形式存入工作文件。双击之，可以在窗口下方显示该值。（4）查看随机数列的频数、频率（百分数）、累积频数、频率（百分数）直方图【例】以Y1为例，查看频数、频率（百分数）、累积频数、频率（百分数）直方图。激活序列窗口，点击View选one-waytabulation功能，点击OK键。上述合并在一起，变成如下：'生成服从某种分布的随机数序列（T=1000），并存入random1文件。workfilerandom1u11000seriesY1=@rnormseriesY2=nrndseriesY3=@rexp(3)seriesY4=@rpoisson(3)seriesY5=@rtdist(3)seriesZ=nrndseriesX=Z*0.5+50scalarq1=@quantile(Y1,.2,1)scalarq2=@quantile(Y1,.5,1)（5）在生成服从某种分布的随机数序列的基础上生成各种ARIMA序列这里主要指在生成服从标准正态分布白噪声序列的基础上生成各种ARIMA序列。【例】在生成N(0,1)白噪声序列的基础上生成AR(1)、MA(1)、ARMA(1,1)、ARIMA(1,1,1)序列workfilerandom2u11000'u表示非时间工作文件，a表示年度文件。seriesu=@nrnd'生成白噪声序列u~N(0,1)u(1)=0'u序列初始值为零。'生成AR(1)序列x1=0.8*x1(-1)+useriesx1'定义x1序列x1(1)=0'定义x1序列初始值为零smpl21000x1=.8*x1(-1)+u'生成AR(1)序列'生成MA(1)序列x1=u+0.8*u(-1)smpl11000seriesx2'定义x2序列x2(1)=0'定义x1序列初始值为零smpl21000x2=u+0.8*u(-1)'生成MA(1)序列'生成ARMA(1,1)序列x3=0.8*x3(-1)+u+0.8*u(-1)smpl11000seriesx3'定义x3序列x3(1)=0'定义x3序列初始值为零smpl21000x3=0.8*x3(-1)+u+0.8*u(-1)'生成ARMA(1,1)序列'生成ARIMA(1,1,1)序列x4=x4(-1)+x3smpl11000seriesx4'定义x4序列x4(1)=0'定义x4序列初始值为零smpl21000x4=x4(-1)+x3'生成ARIMA(1,1,1)序列合在一起写程序如下（比上面节省了许多语句）：workfilerandom2u11000'u表示非时间工作文