非参数模型的稳健跳点检测估计

韩忠成林金官

(1.东南大学数学学院,江苏 南京211189;2.南京审计大学统计与数学学院,江苏 南京211815)

1.引言

观测数据的曲线拟合具有广阔的应用前景,非参数回归模型为曲线拟合问题提供了一个主流的统计工具,其形式为

其中,m(·)是未知回归函数有有界支撑U= [a,b],X是一维解释变量,ε是独立同分布的随机误差项.在某些情况下,回归函数可能在某些未知位置存在跳点,表示相关过程的结构变化.比如,当生产线失控时,产品的质量指标可能在未知的时间点发生向下或向上的移动.在这种情况下,跳点的检测对回归函数结构的刻画十分重要.

近年来,非参数模型跳点的估计已被广泛研究.文[1]指出回归函数可能存在不连续点,传统光滑方法得到的拟合曲线在跳点处存在较大偏差.在跳点个数已知的假设下,文[2]提出了跳点和回归函数的核估计方法.文[3]利用小波方法给出了跳点的检测方法.文[4]利用回归函数的单边非参数回归方法估计不连续点的位置.文[5]基于局部线性估计量构造跳点估计过程,证明了跳点估计过程的收敛性质.文[6]探讨了不同方法下跳点估计问题的最优表现.在实际问题中,跳点的个数和位置通常是未知的.文[7]提出了一种不连续点的检测方法.该方法通过比较任意给定点的三种估计量确定不连续点的位置.文[8]基于局部线性估计量提出了一种保跳曲线拟合方法.文[9]指出局部线性估计量不可避免地存在巨大的计算负担,而B样条在拟合不连续回归函数时表现更好.

上述文献的结果都是在最小二乘方法下得到的.然而,最小二乘方法对观测数据存在异常点或重尾分布的情形十分敏感.众所周知,M-估计常用来处理异常点的情形(见文[10]),但是当误差项服从正态分布时M-估计会损失一些效率.因此,当带跳非参数模型存在异常点时,需要发展一种合适的估计方法能同时获得稳健性和有效性.但是,据知,目前还未有此类研究文献出现.本文在跳点个数和位置未知的假设下,结合B样条提出一个稳健有效的跳检测方法,通过引入一个调节参数,改善回归函数的估计效率.蒙特卡洛模拟和实例分析说明了提出的估计方法不仅在回归函数的连续区间而且在跳点的邻域内都有很好的表现.

本文结构如下: 第2节介绍估计方法;第3节通过数值模拟给出提出的方法在有限样本下的表现;第4节用本文提出的方法处理上证指数数据.

2.估计方法

假设模型(1.1)中的回归函数m(·)有如下表达式:

其中,g(x)是一元光滑函数,I(·)是示性函数当条件为真时取1,否则取0.q表示回归函数中跳点的个数,dj和sj分别表示第j个跳点的幅度和位置.称满足式(2.1)的模型(1.1)为带跳非参数模型.

Ⅰ众数估计

假设{(Xi,Yi),i= 1,··· ,n}是来自模型(1.1)的一组样本.为避免局部多项式估计的缺点,回归函数m(x)可通过B样条近似给出.令U= (u1,··· ,uK)表示支撑[a,b]上的内节点向量,对应的扩展节点向量记为则

其中,B(x) = (B1,p(x),··· ,BK+p+1,p(x))表示p阶B样条基函数,K表示内节点个数.根据众数光滑思想,我们可通过最大化下式

估计α,其中,ϕh(t) =h−1ϕ(t/h),h是需要选择的带宽,ϕ(t)表示对称核密度函数.ϕ(t)的选择不是非常严格,为了便于计算,本文ϕ(t)取标准正态密度.

注意到最大化式(2.3)无法直接得到α的显式解.为了估计α,给出如下的EM算法:

步0 计算α的初始值α(0).设置k=0.

步1 更新π(j|α(k)):

步2 更新α(k+1):

其中,MT= (B(X1),··· ,B(Xn)),Wk是以π(j|α(k))为元素的对角阵,Y= (Y1,··· ,Yn)T.设置k=k+1,并返回至步1.

步3 重复步1至步2,直到收敛.α的最终估计量,记作.回归函数在点x处的估计量记为(x,U∗)=B(x)T.

进一步,如果在U∗内加入p+1个同样的新节点x0∈(a,b),不失一般性,假设x0∈(ui,ui+1),则新的节点向量记为,即

和

类似(2.3)式,(2.5)式和(2.6)式的最优解可通过同样的算法步骤获得,分别记为和.则回归函数在点x处的估计量记为令RSS0表示残差平方和,即插入新节点之后的残差平方和包含两部分

和

注1步0中α(0)的计算可参见文[9]的方法.

Ⅱ跳点检测估计

由文[11]可知,如果回归函数m(x)在支撑[a,b]上是光滑的,则每个设计点(x;U∗)是m(x)的相合估计;如果m(x)在支撑[a,b]上存在跳点,那么在跳点的邻域内(x;U∗)不是m(x)的相合估计.(x;)在区间[a,x0)和[x0,b]上也具有相同的性质.因此,为了提高回归函数的估计精度,需要检测观测数据中的跳点.

为了检测跳点,回归函数估计量的距离函数定义如下:

直观来说,若x0位于回归曲线的连续区域在区间[a,x0) 和[x0,b]上与(x;U∗)相差无几,包括在跳点的邻域内也是如此,所以接近很小;若x0位于跳点的邻域内,仅在x0的左邻域内相合,在x0的右邻域内非相合,而在跳点两侧均是不相合的,因此,当x0接近跳点时,的差异十分显著,D(x0)相应增加.特别地,如果x0与跳点重合,D(x0)可得到局部极大值点.

总体来说,当x0的邻域内存在跳点,D(x0)变大且存在一个局部极大值点,否则D(x0)的值很小.根据D(x0)在跳点处的信息,我们提出下面的跳点检测步骤:

第1 步: 对任一点x0,若满足|D(x0)|≥ϕn,其中ϕn是非负阈值,则x0被标记为跳点.

第2 步: 假设{νi,i=1,··· ,q}是第一步检测的跳点,且∆n=Xi −Xi−1均相等.若存在整数1≤i1∆n且νi2+1−νi2>∆n,则称(νi1+νi2)/2是一个估计的跳点,代替跳点{νi1,··· ,νi2}.

利用上述程序可检测出回归函数中跳点的位置和个数,记作{ν∗1,··· ,ν∗q∗}和q∗.令ν∗0=a,ν∗q∗+1=b,V={ν∗0,··· ,ν∗q∗+1},不难发现,回归函数在区间[ν∗0,ν∗1),··· ,[ν∗q∗,ν∗q∗+1]上是连续的.记新的节点向量为可通过最大化下式

进行估计,其中B∗(x)是节点向量下的B样条基函数向量.与(2.3)式类似,回归函数在点x处的估计量为称为稳健跳点检测估计量.

Ⅲ参数选择

在利用B样条函数拟合回归函数的过程中,有四个参数需要选择: 内节点个数K,基函数阶数p,带宽h和阈值ϕn.首先讨论参数K和p的选择,通常考虑以下二维交叉验证准则

获得.其次,由文[12]可知,基于B样条函数的局部众数估计量与最小二乘估计量的渐近方差之比如下所示:

其中σ2= E(ε2),F(h) = E(ϕ′′h(ε)),G(h) = E(ϕ′h(ε)2).比值R(h)仅依赖带宽h,且在估计量的有效性和稳健性方面扮演着重要角色.因此,带宽h的理想选择为

由(2.7)式可知,hopt与样本大小n无关,只与ε的条件误差分布有关.

实际问题中,随机误差项的分布是未知的,因此F(h)和G(h)无法直接获得.一个灵活的处理方法是通过

分别估计F(h)和G(h).则R(h)可利用来估计,其中表示基于初始估计得到的残差项.利用格点搜索方法,很容易找到hopt最小化(h).

参数ϕn的选择需要合适的跳点检测准则,常用的评价准则为Hausdorff距离

其中J和分别表示真实的和估计的跳点集合.由于J未知,无法直接计算,故采用bootstrap方法.假设存在B个bootstrap样本,根据第k个样本检测到的跳点记为则的估计为

ϕn的最优值可通过最小化获得.

注2参数选择的其他方法可参见文[9,11-12].

3.数值模拟

本节通过数值例子评价提出的跳点检测方法和回归函数估计量的有限样本表现.考虑一组观测值{(Xi,Yi),i=1,··· ,n}来自模型

其中Xi是来自[0,1]的均匀分布,回归函数表达式如下函数m(x)有两个跳点,分别位于0.3和0.7处,幅度分别是2.8和1.7.样本量取n= 200和400,每次实验重复N=200次.误差分布考虑以下两种不同情形:

情形1εi ∼N(0,0.12),正态分布.

情形2εi ∼0.95N(0,0.12)+0.05N(0,32),5%的数据可近似看作异常点.

首先,研究跳点检测方法检测跳点的能力.表3.1给出了不同情形下检测到的跳点出现在真实跳点0.02范围内的次数.与情形2相比,情形1中的跳点检测方法的表现明显更好.这一现象表明误差分布的噪声水平较小,跳点检测方法的表现越好.进一步地,在情形2中,样本量增加相应地提高了跳检测方法检测跳点的能力.同时,当跳点的幅度增加时有类似的结论.

表3.1 200次重复实验下真实跳点0.02范围内检测出跳点的次数

其次,研究回归函数估计量的有限样本表现.在获得跳点个数和位置的估计之后,使用提出的跳点检测方法和众数回归方法(MPS)估计回归函数.为了说明其有效性与稳健性,我们将该方法与基于分段样条拟合和最小二乘提出的跳点检测(LSPS)估计方法[9]进行比较,两个曲线估计量分别记作在200次重复实验下,对这两个估计量计算相应的平均积分平方误差(mean integral squared error,MISE)和跳点附近的局部MISE的值,结果如表3.2所示.

表3.2 回归函数的MISE 和跳点附近的局部MISE 的模拟结果

4.实证分析

股票市场作为国民经济的晴雨表，受到政府和投资者的高度重视。由于股票市场充满了不确定性、机遇和风险,因此,挖掘有效信息可以帮助投资者抓住机遇并规避风险.

股票价格指标是度量金融市场信息的有效工具,从统计学角度分析股票价格指标对获取信息十分重要.作为示例,我们收集了一组上海证券综合指数从2014年1月2日至2016年12月30日的日收盘价数据(见http://q.stock.sohu.com).这三年中,股票市场经历了几次危机,称为中国股市动荡.从图4.1可知,动荡起始于2015年6月15日,于2016年2月早期终止.三个暴跌点出现在2015年6月,2015年8月,2016年1月.然而,由于噪声的影响,跳点位置和幅度均是未知的.因此,跳点检测以及收盘价曲线拟合需要格外关注.值得注意的是,在分析数据之前,有必要对数据进行归一化处理.

图4.1 2014年1月2日至2016年12月30日上海证券综合指数的日收盘价数据

图4.2 2014年1月2日至2016年12月30日上海证券综合指数的拟合曲线

图4.3 Y200 =5000作为异常点时,2014年1月2日至2016年12月30日上海证券综合指数的拟合曲线

根据第2节的跳点检测方法,从图4.2中可观测到三个跳点,分别位于0.483,0.548和0.667(对应日期2015年6月15日,2015年8月21日和2016年1月4日).检测出的跳点位置与三个暴跌点的位置十分接近.同时,图4.2中的拟合曲线与真实数据的变化趋势保持一致,进一步说明提出的跳点检测估计方法在跳点附近和连续区域内的表现良好.

为了检验本文提出的方法对异常点是否稳健,将第200个观测值设为Y200=5000,见图4.3.不难发现,本文提出的方法与最小二乘法的跳点检测方法的跳点检测结果与图4.2中的结果保持一致.当存在异常点的时候,基于最小二乘的跳点检测方法的回归函数估计量(虚线)在异常点附近明显偏离了真实曲线.然而,基于众数的跳点检测方法的回归估计量(虚点线)与图4.2中的结果保持一致.因此,基于众数的跳点检测方法是稳健的.