摘 要:迪基-福勒(Dickey-Fuller)检验是时间序列的平稳性检验中常用的一种方法;由于检验统计量的极限分布是由标准维纳过程关于轨道的积分来表达的，很难得到其密度函数的显式表达式，因而确定检验临界值非常困难，蒙特卡罗方法是解决这类问题的金钥匙;基于R语言对t统计量的平稳性检验的临界值的随机模拟程序进行了研究，填补了文献空白，其计算程序和方法对于金融工程或经济计量统计分析与研究具有广泛的指导意义。
关键词:时间序列分析；平稳性检验；DF检验；蒙特卡罗方法；R语言
0 引 言
迪基-福勒(Dickey-Fuller)检验(DF检验)[1-8]是时间序列的平稳性检验中常用的一种方法。由于检验统计量的极限分布(DF分布)是由标准维纳过程关于轨道的积分来表达的，很难得到其密度函数的显式表达式，因而确定检验临界值非常困难。
起源于20世纪中期的蒙特卡罗(Monte Carlo)方法[9-12]是以概率统计理论为背景的一类重要的数值计算方法。其基本思想是:在某个设定的概率模型中随机抽样，依据所得的样本计算参数的统计特征，从而得到参数的近似值。迄今为止，蒙特卡罗方法在金融工程学、宏观经济学、计算物理学等诸多领域都有广泛的应用。
例（18）可以解读为：约翰从别人身上偷了一台电脑，然后又把这台电脑给予了玛丽”，也就是说“steal”在这里具有“索取+给予”义，而这“给予”义恰好是双及物构式压制的结果。汉语中的“买”也是如此：
R语言是近几年流行的专业从事数据分析和统计绘图的应用软件[11-12]。笔者希望应用蒙特卡罗方法，借助R语言研究DF分布的形态，并对DF检验的分位数进行随机模拟，计算其近似值。尽管不少文献给出了由蒙特卡罗模拟方法计算所得的DF分布的部分样本分位数[3-8]，但都没有给出计算过程或模拟程序。本文所给的模拟计算方法及R程序填补了文献空白，其结果对金融或经济计量分析与统计研究具有广泛的指导意义。
1 单位根过程的DF检验
假设ρ=1，{εt}是平稳过程，即Eεt=0,Cov(εt,εs)=γt-s<∞，称满足
yt=ρyt-1+εt
(1)
的一阶自回归过程{yt}是一个单位根过程。如果ρ=1，{εt}独立同分布，且Eεt=0,Var(εt)=σ2<∞，称满足式(1)的{yt}是一个随机游动过程。如果ρ<1，称满足式(1)的{yt}是一个平稳AR(1)过程。
对平稳AR(1)过程{yt}，假设T是样本容量，令
(2)
分别是ρ的最小二乘估计的标准差的估计以及σ2的最小二乘估计。检验统计假设
H0:ρ=ρ0↔H1:ρ≠ρ0
的统计量tT如下:
(3)
mean(s)
作为烧结配料：钢渣主要化学成分为CaO、FeO等，当CaO含量较高时，钢渣可作烧结矿助熔剂替代部分石灰熔剂，适量配入钢渣可改善烧结矿质量，使转鼓指数和烧结率提高，并且由于钢渣中FeO的氧化放热，降低烧结矿燃耗。
对于原假设H0:ρ=1的检验问题，当H0成立时，满足式(1)的一阶自回归过程{yt}不是平稳的，tT不再服从t分布。当T→∞时，其极限分布为
(4)
其中，{W(r)}是标准维纳过程。此检验问题的另一个统计量亦有类似极限分布：
(5)
它们都是非标准和非对称的分布，统称为DF分布，其密度函数无显式表达式，只能做近似计算。
2 维纳过程及DF分布的统计特征
随机过程{W(t),t∈[0,T]},T>0，如果满足:
(i)W(0)=0；
(ii){W(t),t∈[0,T]}具有独立增量，即对于区间[0,T]的任一分割0=t1<t2<…<tN=1，随机变量W(t2)-W(t1),W(t3)-W(t2),…,W(tN)-W(tN-1)相互独立；
(iii)对任意0≤s<t≤T，都有
W(t)-W(s)～N(0,σ2(t-s))
(6)
则称{W(t)}是[0,T]上的维纳过程(或布朗运动)。若σ=1，则称{W(t)}是[0,T]上的标准维纳过程(或标准布朗运动)。当W(t)看成时间t的函数时，称它是维纳过程{W(t)}的轨道，维纳过程的轨道是t的连续函数。
1.2.2 观察组 在对照组的基础上，为患者发放糖尿病自我管理手册，并针对患者的病情、健康知识的需求，使用手册中相关内容对患者进行宣教。
令Δt=ti-ti-1>0，由式(6)可知：
W(Δt)～N(0,Δt)
W(t+Δt)-W(t)～N(0,Δt)
(7)
注意，由式(7)得到[0,T]上的标准维纳过程的轨道的模拟算法可以按以下步骤进行:
(1) 产生一个标准正态分布的随机变量z；
(2)i=i+1；
(3) 令
rho3<-apply(b,1,sum) # 计算式(2)中第一个分式的分子
至于每一个区间(ti,ti+1)内的值，通常用线段连接。而区间(ti,ti+1)内的值对于模拟计算积分并不重要。
下面用R语言模拟式(4)右端的随机变量(DF分布)的概率密度函数，并绘制直方图和密度函数曲线，考察其分布形态，如图1所示。
>n<-10 000;N<-5 000 # n是样本容量，对[0,1]区间等分成N个小区间
>delta<-1/N # 计算每个小区间长度
>W<-matrix(rep(0,times=n*N),ncol=N,nrow=n) # 产生一个n×N的零矩阵
> for(j in 2:N)
+W[,j]=W[,j-1]+rnorm(n)*sqrt(delta) # 模拟标准维纳过程的n条轨道
alternative hypothesis: two-sided
待水泥砂浆灌注完成后，在水泥砂浆初凝之前，将缓凝剂和水按照一定的比例喷洒在路表面，并根据缓凝时间，在内部水泥砂浆终凝之前，一边用水冲洗，一边用扫帚将路表面水泥砂浆清洗干净，以暴露出母体沥青混凝土为宜，这样既保证路面抗滑性能和均匀性，又确保水泥砂浆具体足够的强度。
> intW2<-apply(W^2*delta,1,sum) # 对矩阵每一行求和，计算标准维纳过程的平方的积分
前苏联教育家苏霍姆林斯基曾说：“让学生变聪明的方法，不是补课，不是增加作业量，而是阅读、阅读、再阅读。”中国少年儿童新闻出版总社低幼读物出版中心对0～9岁儿童及其家长的调查研究表明：7岁之前是养成阅读习惯的重要阶段，如果在此时没有养成良好的阅读习惯，7岁之后的孩子就很难掌握和提高阅读的能力。对6～12岁的青少年来说，小学阶段是形成和掌握基本阅读能力的关键期。在这个关键期让学生爱上阅读，并养成良好的阅读习惯，会使他们终生受益。
>stat<-1/2*(W1^2-1)/sqrt(intW2) # 得到DF分布的n个值
>hist(stat[stat<5],freq=F,ylim=c(0,0.7),breaks=seq(-5,5,length=20)) # 绘制直方图
>lines(density(stat),col="red",ylim=c(0,1)) # 添加密度函数曲线
二值法从本质上仍是提取等值线，但与等值线追踪法相比，解决了海岸线零乱、破碎的问题，通过解决了海岸线曲折的问题，在保证精度的同时可省去大量的人工操作，自动化程度高，自动提取高。
图1 DF分布的密度函数曲线模拟图
Fig. 1 Simulated curve of density function for DF distribution
>mu<-mean(stat) # 数学期望
>mu
[1] 0.01 237 375
管世铭生活的乾嘉时期，正值清王朝由盛而衰之际。一方面，国运凋敝，民生多艰；另一方面，边疆地区各部族叛乱时有发生。乾隆在位六十年间，清王朝为维护大一统的战争多达上百起，管世铭创作了多篇以这一时期重大战争为题材的诗歌。
>sd<-sd(stat) # 标准差
>sd
[1]1.669 781
教师展示洋葱表皮细胞、心肌细胞等多种细胞的显微图片，使学生从熟悉的图像中，直观感受细胞核是多数细胞的共同结构，并初步感知细胞核的大致形状和位置。教师提问：“是否所有的细胞都有细胞核呢？”学生联系旧知，回顾制备细胞膜的选材，意识到：哺乳动物成熟红细胞没有细胞核。教师补充关于高等植物成熟筛管细胞的内容，使学生形成对细胞核的初步认知，渗透批判性思维。
>ks.test(stat,"pnorm",mu,sd) # 用K-S检验法进行正态性检验表明不是正态分布
One-sample Kolmogorov-Smirnov test
data: stat
D=0.18 755, p-value<2.2e-16
>W1<-rnorm(n) # 注意W1是标准正态分布
用以上蒙特卡罗模拟方法可类似研究式(5)右边的极限分布的统计特征。
根据国内外农田水利现代化发展经验和建设实践的探索，参照国家颁布的相关技术规范与标准，结合郑州市实际，提出农田水利现代化建设内容，主要包括节水灌溉工程、除涝工程、雨水集蓄利用工程、农田园田化工程以及农田水利信息化工程五个方面。
鼓励金融机构为乡村旅游发展提供信贷支持，创新金融产品，降低贷款门槛，简化贷款手续，加大信贷投放力度，扶持乡村旅游龙头企业发展。依法合规推进农村承包土地的经营权、农民住房财产权抵押贷款业务，积极推进集体林权抵押贷款、旅游门票收益权质押贷款业务，扩大乡村旅游融资规模，鼓励乡村旅游经营户通过小额贷款、保证保险实现融资。鼓励保险业向乡村旅游延伸，探索支持乡村旅游的保险产品。
3 蒙特卡罗模拟计算DF检验的临界值
对一般的单位根过程可以用菲利普-佩荣(Phillips-Perron)检验法或增广迪基-福勒(ADF)检验法[1-3]。现在考虑随机游动过程{yt}，并假定{εt}独立同N(0,σ2)分布，其参数的最小二乘估计由式(2)给出。对单边检验问题:
H0:ρ=1↔H1:ρ<1
(8)
其检验统计量tT如式(3)定义。因此，对给定的显著性水平α，检验临界值uα确定如下:
P(tT<uα)=α
当检验统计量tT<uα时应拒绝原假设。下面借助R语言用两种方法近似计算临界值uα。
第一种方法是从模型中抽取样本用蒙特卡罗模拟计算，R程序如下:
>quant1<-function(p,m=1 000,T=200){
eps<-mapply(rnorm,rep(m,T)) # 从标准正态分布中抽样，产生m×T的矩阵
>abline(v=0,col="blue"，lwd=2)
y<-t(apply(eps,1,cumsum)) # 对矩阵eps的每一行求累积和得到矩阵y
rho2<-apply(y[,-T]^2,1,sum) # 对y的每一行去掉最后一个元素后求平方和
rhoT<-rho3/rho2 # 计算ρT的最小二乘估计值
for(t in 2:T)
b[,t]=y[,t]*y[,t-1]
(4) 如果i≤N，则重复第(1)步。
b<-matrix(rep(0,m*T),nr=m,nc=T)
sigma1<-matrix(rep(0,m*T),nr=m,nc=T)
for(t in2:T)
sigma1[,t]=(y[,t]-rhoT*y[,t-1])^2
sigma2<-apply(sigma1,1,sum)
sigma<-sqrt(sigma2/(T-1)) # 求σT的估计值
t_T<-(rhoT-1)*sqrt(rho2)/sigma # 求统计量tT的值
quantile(t_T,probs=p) # 输出1个样本p分位数值
优化幼儿园一日活动科学保教实践需要充分认识到保教活动的重要价值，了解其对幼儿成长的重要作用，从而合理设计幼儿园一日活动。具体来讲，可以从以下几个方面加以思考。
+}
利用本文开发的“薄壁空心高墩温度应力分析系统”，只需要输入或改变必要的参数，便可以很方便、很直观地实现在任意时间和地点、任意季节以及任意截面形状等条件下的高墩结构在日照温差作用下的温度效应的计算。计算流程如图3所示。
>quant<-function(p,n=1 000){
n<-1 000
s<-replicate(n,quant1(p)) # 调用函数quant1重复计算n次产生n个p分位数
动员各方面力量搞好村庄环境宣传，把全体村民发动起来，让他们都真正理解支持，是每个村民都能做到爱村、护村，并行动起来，投入到美丽村庄建设上来，形成合力势头和建设发展力。
其中，ρ0<1。统计量tT服从自由度为T-1的t分布，其检验临界值容易获得。
}
>quant(0.025)
[1]-2.22 641
>quant(0.95)
[1]1.287 814
在文献[4]的第593页“附表G”(Table G Empirical Cumulative Distribution of TforΦ=1)可以看到，当n=1 000时，对应的0.025分位数是-2.23，对应的0.95分位数是1.28，与前面得到的分位数-2.226 41及1.287 814非常接近。
第二种方法是从极限分布中抽取样本用蒙特卡罗模拟计算，R程序如下:
>intofW<-function(n,delta=1/1 000){ # 此函数的目的是从DF分布中抽取n个样本
N=1/delta # 将区间[0,1]分割成N等分
W<-matrix(rep(0,times=n*N),ncol=N,nrow=n) # 产生一个n×N的零矩阵
for(j in 2:N)
W[,j]=W[,j-1]+rnorm(n)*sqrt(delta) # 模拟标准维纳过程的n个轨道
intW1<-apply(W*delta,1,sum) # 计算标准维纳过程在[0,1]上的积分
intW2<-apply(W^2*delta,1,sum) # 计算标准维纳过程的平方在[0,1]上的积分。
1/2*(intW1^2-1)/sqrt(intW2) # 输出DF分布的n个样本值
}
>quant<-function(p,n=1 000,m=1 000){ # 此函数产生DF分布的p分位数的近似值
t<-mapply(intofW,rep(n,m)) # 调用intofW函数重复计算产生n×m的样本矩阵
q<-apply(t,2,quantile,probs=p) # 对矩阵t的每一行求得一个样本的p分位数
mean(q) # 将所得的n个样本分位数取平均，输出DF分布的p分位数的近似值
}
>quant(0.025)
[1]-2.351 553
>quant(0.95)
[1]0.1 069 714
4 结束语
第一种方法是从模型中抽取样本，得到统计量的值与样本容量n有关。文献[4]P593附表G与文献[5]P642附表10.A.2所给出的临界值是用第一种方法计算的，即Empirical Cumulative Distribution(经验累积分布)。第二种方法是从极限分布中抽取样本求得的分位数，两种方法不同，所得结果有一点偏差是可以理解的。
将以上程序进行微小修改便得到统计量的相应计算程序。至于DF检验还有其他类型或ADF检验的各种类型的蒙特卡罗模拟，将所给的方法及程序作相应修改即可得到。
致谢：衷心感谢西藏民族大学财经学院汪朋博士在作者写作过程中所给予的热情帮助！
参考文献:
[1] 史代敏，谢小燕. 应用时间序列分析 [M]. 北京:高等教育出版社，2011
SHI D M, XIE X Y. Applied Time Series Analysis [M]. Beijing: Higher Education Press, 2011(in Chinese)
[2] 王黎明，王琏，杨楠. 应用时间序列分析 [M]. 上海:复旦大学出版社，2009
WANG L M, WANG L, YANG N. Applied Time Series Analysis [M]. Shanghai: Fudan University Press, 2009 (in Chinese)
[3] 陆懋祖. 高级时间序列经济计量学 [M]. 北京:北京大学出版社，2015
LU M Z. Advanced Time Series Econometrics [M]. Beijing:Peking University Press, 2015(in Chinese)
[4] WILLIAM W S. Time Series Analysis: Univariate and Multivariate Methods[M]. New York: Pearson Addison Wesley, 2006
[5] FULLER W A. Introduction to Statistical Time Series[M]. New York: John Wiley & Sons, Inc, 1996
[6] DICKEY D A. Estimation and Hypothesis Testing in Non Stationary Time Series [D].Iowa: Iowa State University, 1976
[7] DICKEY D A, FULLER W A. Distribution of the Estimators of Autoregressive Time Series with a Unit Root [J]. Journal of the American Statistical Association, 1979, 74: 427—431
[8] 詹姆斯·汉密尔顿，时间序列分析 [M].夏晓华，译.北京:中国人民大学出版社，2015
HAMILTON J D.Time Series Analysis [M]. XIA X H，Translated.Beijing: China Renmin University Press, 2015(in Chinese)
[9] 徐钟济. 蒙特卡罗方法 [M]. 上海:上海科学技术出版社，1985
XU Z J. Monte Carlo Methods [M]. Shanghai: Shanghai Sciences and Technology Press, 1985(in Chinese)
[10]李东风. 统计计算 [M]. 北京:高等教育出版社，2016
LI D F. Statistical Computing [M]. Beijing: Higher Education Press, 2016(in Chinese)
[11]薛毅，陈立萍. 统计建模与R软件 [M]. 北京:清华大学出版社，2007
XUE Y, CHEN L P. Statistical Modeling and R Software [M]. Beijing: Tsinghua University Press, 2007(in Chinese)
[12]STEFANO M I. Simulation and Inference for Stochastic Differential Equations with R Examples[M]. New York: Springer, 2008
Monte Carlo Simulation of Dickey-Fuller Test Based on R Language
ANJun
(School of Mathematics and Statistics, Chongqing Technology and Business University, Chongqing 400067, China)
Abstract:Dickey-Fuller test is a commonly used method in the stationary test of time series.Because the limit distribution of its statistics is expressed by the integral of the standard Wiener process about its trajectory. It is very difficult to obtain the explicit expression of its density function, so it is very difficult to determine the critical value of test. Monte Carlo method is the golden key to solve this kind of problem. Based on R language, we study the critical value of stochastic simulation program for statistic stationarity test, which fills in the blank of literature. Its calculation program and method have a wide range of guiding significance for financial engineering or econometric statistical analysis and research.
Keywords:time series analysis;test for stationary; Dickey-Fuller test; Monte Carlo method; R Language
doi:10.16055/j.issn.1672-058X.2019.0003.003
中图分类号:C32，O242.1
文献标志码:A
文章编号:1672-058X(2019)03-0014-04
收稿日期:2019-01-7
修回日期:2019-03-01.
* 基金项目:重庆市教委自然科学基金项目(KJ130705)；重庆工商大学经济社会应用统计重庆市重点实验室开放基金项目；重庆工商大学校级教改课题(2018214).
作者简介:安军(1964-)，男，四川安岳人，副教授，从事概率论及数理统计研究.
责任编辑:罗姗姗
引用本文/Citethispaper:
安军.基于R语言的迪基-福勒检验的蒙特卡罗模拟[J].重庆工商大学学报(自然科学版)，2019,36(3):14—17
AN J.Monte Carlo Simulation of Dickey-Fuller Test Based on R Language[J].Journal of Chongqing Technology and Business University (Natural Science Edition),2019,36(3):14—17
标签：过程论文; 样本论文; 临界值论文; 蒙特论文; 方法论文; 社会科学总论论文; 社会科学研究方法论文; 统计方法论文; 计算方法论文; 《重庆工商大学学报(自然科学版)》2019年第3期论文; 重庆市教委自然科学基金项目(KJ130705) 重庆工商大学经济社会应用统计重庆市重点实验室开放基金 重庆工商大学校级教改课题(2018214)论文; 重庆工商大学数学与统计学院论文;