第37卷 第6期2009年11月河南师范大学学报(自然科学版)J our nal of H enan N or mal Univer sity (N atur al Science ) Vol.37 N o.6 N ov.2009 文章编号:1000-2367(2009)06-0007-04最小二乘法对多变点检验的性能研究张学新1,段志霞2(1.中南财经政法大学信息学院,武汉430060;2.济源职业技术学院基础部,河南济源459650)摘 要:给出了衡量最小二乘法识别多变点能力的方法,模拟研究了最小二乘法对不同数据生成过程的多变点检测效果,指出了最小二乘法的适用性,最后应用最小二乘法检测了中国主要经济部门的GDP 变点.关键词:最小二乘法;多变点检验;单位根过程;蒙特卡罗模拟中图分类号:C812文献标识码:A数据生成过程的结构突变是指系统受到诸如金融危机,体制变化等剧烈的外力冲击而发生的突然变化,是系统对外界条件的光滑变化而做出的突然响应,常见的有均值突变,频率突变,趋势突变,方差突变.突变分析,尤其是带单位根过程的突变分析是国内外比较热门的课题,各种变点检验的方法也在不断涌现.目前的研究大多集中于前述常见类型的突变,国外见文献[1-6],国内主要是各种方法在气候,交通等领域的应用,见文献[7-9]等.其中有些统计方法是有缺陷的,如用滑动t ,滑动F 检测法检测均值突变时,经常会检测到一些虚假的突变点,主要问题是不能确定突变的研究通常涉及到的非独立随机变量的分布.最小二乘法也是处理变点问题中使用较多的一种方法,它以观察值与理论值之差的平方和作为目标函数,以其达到极小值之点作为有关参数的点估计,其优点是对随机误差的分布不需要作特定的假设.国内文献鲜见研究最小二乘法识别多变点的性能,国外至多用最小二乘法讨论了误差为线性过程时一个未知均值变点的估计问题(见前述文献中的JushanBai).本文通过模拟对最小二乘法识别多变点的性能做较为详尽的研究.1 均值变点的最小二乘法估计设离散的模型是X i =a i +e i ,e i ~iid ,E(e i )=0,Var (e i )= 2,i =1, ,n,a m j =a {m j +1= =a m j +1-1=b j+1,j =0,1 ,q.这里q 是事先给定的变点个数,可以取充分大以满足实际要求,或者通过其它方法粗略估计得到.1=m 0<m 1< <m q <=n,m q +1=n+1,定义m i 是一个变点,当b i+1 b i 时,且在该点跃度是b i+1 b i .令T =T (m 1 ,m q ,b 1 ,b q+1)=!q+1j=1!m j -1i=m j-1(x i -b j )2,极小化它,求出未知的m 1 ,m q ,b 1 ,b q+1的估计值.易见,当固定m 1, ,m q 时,上式在b j =X m j-1+X m j-1+1+ +X m j -1m j -m j-1时达到最小值,因此极小化目标函数T =T(m 1 ,m q )=!q+1j=1!m j -1i=m j -1(x i -X m j-1+X m j -1+1+ +X m j -1m j -m j-1)2即可.引理1 设两序列X m , ,X i-1与X i , ,X n 的算术平均数分别为 X i 1, X i 2,则当(i-m)∀(n-i+1)最大时,S i =!i-1t=m (X t - X i 1)2+!nt=i (X t - X i 2)2达到最小.引理2 设序列X m , ,x n ,S i =!i-1t=m (X t - X i 1)2+!n t=i (X t - X i 2)2,S *=min (S m+1, ,S n ),C =收稿日期:2009-05-20作者简介:张学新(1966-),男,湖北宜城人,中南财经政法大学博士研究生,研究方向:概率论与数理统计方法应用.S *(2ln ln (n -m +1)+ln ln ln (n -m +1)-ln ( )-2ln (-0.5ln (1- )))n -m +1-2ln ln (n -m +1)-ln ln ln (n -m +1)-2.4,当S -S *>C 时,认为变点存在,且该检验有渐近水平 .引理2由文献[10]推广得到,用它检验变点存在与否.一旦确认存在变点,则遵循以下步骤极小化目标函数[10]:(1)取定一组初始值m 1, ,m q ,1=m 0<m 1< <m q <=n,m q+1=n +1;(2)在约束条件1<m 1<m 2下,求m 1,使第一,二两项之和W 12=!m 1-1i=1(x i -X 1+ +X m 1-1m 1-1)2+!m 2-1i=m 1(x i -X m 1+ +X m 2-1m 2-m 1)2达到最小,记所得的m 1为m #1;(3)在约束条件1<m 1<m 2<m 3下,求m 2,使第二,三两项之和W 23=!m 2-1i=m #1(x i -X m #1+ +X m 2-1m 2-m #1)2+!m 3-1i=m 2(x i -X m 2+ +X m 3-1m 3-m 2)2达到最小,记所得的m 2为m #2这样继续下去,得到一组新值m #1<m #2< <m #q ;(4)把它们作为初始值回到第一步,继续下去得到一组新值m ##1<m ##2< <m ##q ,再回到第一步,一直继续这个过程,直到新值与上一次的值完全相同时为止,记最后所得的值为m ^1, ,m ^q ,它就是变点m 1, ,m q 的估计.此时T 的最小值记为T q =T(m ^1, ,m ^q ).实际应用时变点个数q 是未知的,但总可取充分大的q 作为上限,再设定一个比1稍大的值比如1.1,找出使T k T q 1.1成立的最大的k,把它作为q 的估计.能用图像等方法预先设定更好.2 最小二乘法的性能比较研究2.1 衡量最小二乘法检验性能的方法统计假设检验时,由于样本的随机性,可能会犯两类错误,第二类错误是指当原假设H 0不真时,样本观测值没有落入拒绝域W ,从而没有拒绝原假设H 0.把不犯第二类错误的概率1-!称为检验的功效.在原假设(∃没有变点%)H 0:b 1=b 2= =b q +1下,变点估计(m ^1, ,m ^q )的分布是什么,目前理论上尚无答案,也就没有用精确的置信系数和置信区间估计来评价检验功效的办法.现引入欧氏距离d =&MM ^ &,其中M(m 1, ,m q )是真实变点的位置,M ^(m ^1, ,m ^q )是它的估计值.显然d 越小越好.选取∀,设定蒙特卡洛模拟次数N ,计算变点的估计落在区间d ∀的次数N rec ,则最小二乘法的识别能力(性能)可定义为:Pow er =N rec N.若取∀=1,2,22,其含义分别是点M ^(m ^1,m ^2)与点M(m 1,m 2)的对应坐标,一个完全相同,另一个只相差1;两个各相差1,或一个相同,另个相差2;两个都相差2.这种误差,在应用上可以容忍.下面做模拟研究,为方便,考虑两个变点情形,多个变点情形完全类似.2.2 不同数据生成过程的变点检验为比较,设定所有模拟次数均为1000次,取3段数据个数各n =50,总个数n =150,真实变点设为M(51,101),当误差标准差为1时,跃度是0.1,意指约为误差标准差的35%.若非声明,以下所说的结果都指模拟结果,且是针对∀=4(d 42)而言.均值突变 数据生成过程是X i =a i +e i ,e i ~iidU(-0.5,0.5),i =1,2,3,取a 1=0.1,a 2=0.2,a 3=0.3,结果只有4%,究其原因,是在每个变点处的跃度太小,数据近似平稳过程.但是使其他条件不变,只把每个变点处的跃度提高到约为误差标准差的70%,则结果提高到36%,假若跃度再提高到约为误差标准差的138%,则结果提高到87%.类似的,取数据生成过程是X i =a i +e i ,e i ~iid N (0,1),i =1,2,3,仍然使每个变点处的跃度依次约为误差标准差的35%,70%,138%(a 1=0.1,a 2=1.485,a 3=2.87),则结果依次为不超过5%,提高到29%,提8河南师范大学学报(自然科学版) 2009年高到88%以上.可见,最小二乘法的检验功效与各个信噪比|a i+1-a i |有很大关系,与误差项的分布关系不大;各个变点处的跃度越大,d 越小; 越小,d 越小,最小二乘法的性能越好.方差突变 设数据生成过程是X i =a i +e i ,e 1~iid N (1,1),e 2~iid N (1,2),e 3~iid N (1,3),结果只有0.2%,最小二乘法的性能较差.不含时间趋势的截距突变的单位根过程 数据生成过程是y t =#0+y t-1+#1D U 1t +#2DU 2t +e t ,e t ~iid N (0, 2),其中D U 1t =I (t >T b 1),D U 2t =I (t >T b 2),分别代表在时点T b 1,T b 2均值(截距)发生突变,改变量分别是#1,#2.取#0=1, =1,每个变点处的跃度约为误差标准差的1.38倍(#1=1.485,#2=2.87),结果为0,究其原因,当数据生成过程为单位根过程时,用于检验的模型已不是均值突变模型,为此,应向均值突变模型转换,先对数据{X t }进行差分,然后再用最小二乘法检验,结果竟达到99%以上,两个估计点与真实位置的误差均不超过1(d 2)的结果亦达60%以上,这与前面的均值突变模型的检验结果较吻合.含时间趋势的截距突变的单位根过程 数据生成过程是y t =#0+∃0t+y t-1+#1DU 1t +#2DU 2t +e t ,e t ~iidN (0, 2),其中DU 1t =I (t >T b 1),DU 2t =I (t >T b 2)的含义同前.对此种模型的检验,先对差分序列作时间退势回归,再对残差作最小二乘检验.取#0=1, =1,∃0=0.5,每个变点处的跃约为误差标准差的138%,结果几乎为0,究其原因,参数估计值^∃0向右偏离真值较远(多数在0.54左右),当估计改进为^∃0=0.52时,重新做上述模拟,结果改善到65%以上.含时间趋势的截距和斜率双突变的单位根过程 数据生成过程是y t =#0+∃0t +y t-1+#1DU 1t +#2D U 2t +∃1{DT 1t *+∃2{DT 2t *2+e t ,e t ~iid N (0, 2),其中DT 1*t =(t-T b 1)DU 1t ,{DT 2*t =(t-T b 2)D U 2t ,分别代表伴随均值发生突变,在时点T b 1,T b 2上斜率也发生突变,改变量分别是∃1,∃2.若模型中e t 取成一般的ARM A (p ,q)过程,则意味着对趋势函数的冲击是渐进完成的,这可能更符合现实经济运行情况.若考虑加入!pj =1c j %y t-j 又得如下模型,y t =#0+∃0t +y t-1+#1DU 1t +#2DU 2t +∃1DT 1*t +∃2{DT 2*t +!pj=1c j %y t-j +e t ,各变量含义同前.这两种模型的斜率随时间变化,不宜对时间作退势回归了.结论 最小二乘法对均值变点模型较适用,运用于其他模型时须向此方面转化.3 实 例考虑对1952-2003年中国主要经济部门GDP 的变点检验.为真实可信,数据特别取自文献[11],各年GDP 的估计是农业,工业,建筑业,交通与通讯业,商业,非物质服务业的GDP 估计加总.以{y t }表示GDP 序列,在作回归系数的显著性检验后常数项,时间趋势项均被剔除,再作单位根检验得AR (2)模型 %2^y t =0.044y t-1-0.403%y t-1(3.27)(-2.45)R 2=0.175, AI C =24.88 SC =24.96 ADF =3.27 D.W =1.97,或者y ^t =1.642y t-1-0.5974y t-2,括号()内是t 值,对应的概率P 值分别是0.0020,0.0181,所有信息表明模型拟合得较好,(在建模过程中发现中国GDP 过程是近似含时间趋势的两个单位根过程,y ^t =-1360.655t+1.603y t-1-0.525y t-2,时间变量t 值及对应的概率P 值分别是-1.654,0.105,t 的显著性检验不能通过显得较为勉强,若原文的数据再长一些,结论应会改变).若直接对二阶差分序列{%2yt }作单位根检验,ADF =-3.26,对应的P 值是0.0017,在1%的检显著性水平下,拒绝原假设,也得到{y t }含有两个单位根的结论.用最小二乘法对{%2y t }变点检测:(任意)设想5年左右有一个变化,则52个数据至多有11断点,分别计算T i ,得(T 1,T 2, ,T 10,T 11)=(15.3,14.9,7.79,6.34,6.05,5.82,5.67,5.45,5.34,4.62,4.42)∀1010,取q 为使T k T 111.1成立的最大的k =9,得变点估计M ^(9,10,33,34,35,37,43,44,46),进一步确认,得M ^(9,33,35,37,44,46),对应的年份是1962,1986,1988,1990,1997,1999;若从{T k }下降的梯度何9第6期 张学新等:最小二乘法对多变点检验的性能研究10河南师范大学学报(自然科学版) 2009年时开始一直趋于平缓考察,则得到5个年份:1962,1964,1987,1990,1997,1999.这两个检测结果与当年的经济实际运行情况非常吻合.查历史文献可知,1962年是∃压缩支出,平衡预算,加强财政管理%年,是政府工作报告中提出∃必须用几年的时间幅度调整国民经济%年;1964年则是农业学大寨,是大搞农田基本建设的一年;1987年是把农村改革引向深入,强调∃深化企业改革,压缩过度膨胀的预算外投资规模和过高的非生产性投资%的一年,是进一步加快和深化改革的一年;1990年十三届七中全会召开,是∃企业改革为重点%,加∃积极稳妥地推进粮食流通体制的改革%的一年;1997年则是∃农业和农村工作要着力做好八个方面的工作%的一年,是提出%三个有利于%的所有制形式,∃调整和完善所有制结构,进一步解放和发展生产力%的一年,是以∃国有企业改革为经济体制改革的重点%的一年;1999年则是亚洲金融危机影响扩散的一年,也是中央加快中西部地区发展的一年.参 考 文 献[1] Bai J.Least squares estimation of a shift in lin ear p roces ses[J].Journal of Tim e Series Analysis,1994,15:453-472.[2] Per ron P.T es tin g for a U nit Root in a Time Series w ith a Changing M ean[J].J ournal of Busin ess an d Economic Statistics,1990,8:153-162.[3] LanneLann e M,L tkep oh l H,Saikk on en paris on of unit r oot tests for time s eries w ith level s hifts[J].J ou rnal of Tim e S eries A!nalysis,2002,23:667-685.[4] Saikk on en P,L tk epohl H.T esting for a un it root in a tim e s eries w ith a level s hift at un know n time[J].Econometric T heory,2002,18:313-348.[5] Lavielle M,M ou lines E.Least!squares estimation of an u nknow n numb er of shifts in a time s eries[J].Journal of Time Series Analysis,2000,21:33-59.[6] Clemen te J,M onta s A,Reyes M.Testing for a unit root in variables w ith a double change in th e mean[J].Economics Letters,1998,59:175-182.[7] 齐培艳,田 铮.噪声为单位根过程的非参数函数变点的小波检测[J].控制理论与应用,2009,26(1):57-61.[8] 张建军,周后福.合肥气温和降水的突变特征分析[J].安徽农业科学,2007,35(9):2724-2726.[9] 龚志强,封国林.非线性时间序列的动力结构突变检测的研究[J].物理学报,2006,55(06):3180-3186.[10] 陈希孺.变点统计分析简介[J].数理统计管理,1991,10(2):52-53.[11] 安格斯麦迪森.中国经济的长期表现:公元960-2030年[M].伍晓鹰译.上海:上海人民出版社,2008:158.S imulation Study on the Power of Least!squares Test for Detecting Multiple BreaksZH ANG Xue!xin1,DUAN Zhi!x ia2(rmation S chool,Zhongnan University of Economics and Law,W uhan430060,Ch ina;2.Department of Basic,J iyuan Vocational and Techn ical C ollege,Jiyuan459650,C hina)Abstract:In this paper,a no vel scheme to measure the pow er of least!squares test fo r checking mult iple br eaks is pr esen! ted as well it s a pplicability to v ario us stream data pr ocessing sy stems discussed firstly,A lso a co nclusio n is g iv en that L east Squa re method can be suitable for mean break model,that sig nal!to!no ise r atio clo sely related to the test pow er w ill be indicated too,Finally,the GD P breaks from China's major eco no mic secto rs ar e detected}.Keywords:the po wer of least!squar es;multiple br eaks test;unit r oot pr ocess;M onte Carlo simulation。