变系数的周期性时间序列模型及其应用.pdf

应用概率统计第 40 卷第 1 期2024 年 2 月Chinese Journal of Applied Probability and StatisticsFeb.,2024,Vol.40,No.1,pp.50-74doi:10.3969/j.issn.1001-4268.2024.01.004变系数的周期性时间序列模型及其应用方学莉(上海财经大学统计与管理学院,上海,200433)王守霞(北京大学数学科学学院,北京,100871)摘要:存在于各个领域的时间序列不仅表现出周期性的特征还易受外界因素的影响,而且外界因素的影响并非一成不变,同时,部分时间序列的周期是未知的.对于这样的易受外界因素影响的周期性时间序列,本文旨在构造含有变系数函数的周期性序列模型.将经典的时间序列模型分解成一个含有未知参数的部分线性变系数模型,利用 B 样条逼近外生变量的变系数函数,借助带有 l0惩罚项的最小二乘回归得到未知周期、周期序列以及外生变量的影响系数的估计结果.本文还给出了估计量的理论性质,包括周期估计的相合性、周期序列估计和变系数函数估计的渐近性质.通过第 4 章的模拟,我们展现了本文方法的优越性.最后我们通过三个实际数据的应用展现了本文方法的实用性.关键词:周期估计;l0惩罚;B 样条;变系数模型中图分类号:O212.7英文引用格式:FANG X L,WANG S X.Varying-coefficient model and applications for the periodictime seriesJ.Chinese J Appl Probab Statist,2024,40(1):5074.(in Chinese)1引言由于时间序列数据存在于各个领域,对其研究的方法也在不断完善.早期学者们对于时间序列的研究主要停留在平稳的时间序列,随着时间的推移逐步转移到包含趋势项、周期项以及外生变量的时间序列.表现出周期性并易受外界环境影响的时间序列广泛存在于现实生活的各个领域,不仅存在周期已知的数据,比如月度旅游数据、季度数据等,而且存在周期未知的数据,比如全球年均气温数据、某地区年降雨量数据、天文学上太阳黑子的变化数据等.一方面,某些外生性变量的影响会随时间变化,不可忽视.另一方面,外生性变量的存在对时间序列的周期估计会有一定影响,同时周期的存在也会影响外生性变量系数的估计.因此,如何准确的估计出序列的未知周期对时间序列的后续研究有着非常重要的影响.经典的时间序列分解模型是将序列分解成趋势项、周期项和误差项三部分,然后运用参数法估计趋势,或者采用差分方法剔除周期和趋势,但没有对外生性变量进行进一步研究.对于外生性变量的添加,近来不少学者采用机器学习或深度学习的方法来进行系数估计(参见文献 1、文献 2),但机器学习的方法对周期分量分离效果不如传统分解模型,有中国博士后科学基金项目(批准号:2023M730090)资助.通讯作者,E-mail:.本文 2021 年 7 月 28 日收到,2022 年 2 月 18 日收到修改稿.第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用51时也会缺少解释力度.于是考虑是否可在经典分解模型中添加协变量项进行估计.经典线性回归模型是我们最常用的也是最简单的一种引入外生性变量的统计模型,但线性回归模型假定系数是常数,它要求解释变量对响应变量的影响是不变的,而在实际问题中,外生性变量对响应变量的影响不是一成不变的,所以文献 3 提出了变系数回归模型,将经典线性模型的系数变成系函数,从而有了更好的解释能力.文献 4 提出了变系数模型的局部多项式估计方法;文献 5 运用样条方法对部分线性变系数模型的估计问题进行了研究;文献 6采用剖面最小二乘法研究了部分线性变系数模型的估计和检验问题.本文选择用 B 样条方法来估计外生性变量的变系数函数,与其他大部分估计方法相比,B 样条方法具有计算方便、局部控制等许多优点(参见文献 7),并且 B 样条估计达到了非参数估计的最优收敛速(参见文献 8).对于周期未知的时间序列的周期估计问题,有大量的文献对其进行研究.传统的方法是参数法,也就是基于频域分析的三角函数回归法,利用有限个三角函数的线性组合来表达周期项,然后基于周期图或者谱密度去估计未知参数.文献 9 在误差是独立同分布的假设下运用三角函数回归法估计序列的频率;文献 10 运用加权最小二乘法估计频率;文献11 运用极大似然估计法估计时间序列的频率.近年来,非参数周期估计法也在不断发展,对于随机非等距观测的时间序列,比如 Yt=m(Xt)+t,此时周期项 m 是一个周期函数,文献 12 运用非参数核最小二乘回归法估计序列的周期;文献 13 运用非参数核最小二乘回归法对多周期的问题进行了研究;文献 14 采用非参数周期图法对单周期以及多周期问题进行了研究;文献 15 研究了变周期的序列的周期估计问题.对于等距观测的时间序列,周期项是一个周期序列,而不是周期函数,文献 16 将周期估计问题看作模型选择问题,将模型假定为 Yt=m(t)+t,运用交叉验证方法来估计未知周期.但以上文献研究的模型只存在周期分量和误差,当模型存在趋势项和周期项时,也就是经典时间序列分解模型 Yt=g(t/T)+m(t)+t,其中 g()是光滑趋势函数,文献 17运用惩罚最小二乘法和局部线性核估计法对存在光滑趋势的周期序列进行了研究.但文献17 的方法在估计周期时直接忽略趋势项的影响,将模型直接看作 Yt=m(t)+t来估计周期,得到周期的估计后利用局部线性核估计法估计趋势.但是直接忽略趋势项的影响会增大估计量的方差,而且当趋势变化较大时,周期估计量不再准确,这种方法不适用于趋势变化较大的序列,而且他们的模型也没有考虑外生变量的影响.文献 18 利用 B 样条估计改进了文献 17 的方法,但他们的模型也没有考虑任何外生变量.文献 19 使用贝叶斯搜索方法来估计非平稳时间序列中的周期性和变点,但他们的方法依赖于趋势、周期分量甚至误差项的参数假设.文献 20 提出了一种基于正则化奇异值分解的新季节调整方法来灵活捕获两种季节性:不随时间变化的固定季节性和随时间变化的季节性,但他们并未考虑外生变量的影响.而文献 21 和文献 22 则研究了函数型时间序列的周期存在性的检验问题,但他们也假定周期已知.文献 23 采用非线性可加模型研究了外生变量对周期序列的影响,但他们假定周期已知.文献 24 研究了存在趋势和线性协变量效应的周期性时间序列的变点检验问题,但他们是在参数模型假设下进行研究,假定趋势函数具有某种参数结52应用概率统计第 40 卷构并且协变量效应被假定是线性的.然而很多时间序列不仅存在未知周期,同时序列还受其他外生变量的影响,外生变量的影响可能不是线性的,所以他们的方法不适用于极易受到外界因素干扰的周期性时间序列.故为了克服这个局限性,本文改进了模型,利用半参数变系数模型将外生性因素考虑进来,并且允许外生变量对于序列的影响随着时间变化.在估计未知周期时没有像文献 17一样忽略趋势项及其他外生变量,而是利用 B 样条逼近未知变系数函数,然后运用带有 l0惩罚项的最小二乘估计法估计序列的周期,进而得到周期序列和变系数函数的估计.本文给出了估计方法的渐近理论,并通过模拟和实证展现了该方法的优越性和实用性.本文接下来的内容作如下安排:第 2 章介绍了我们的模型及估计方法,将时间序列分解为周期项、外生变量影响项以及误差项,接下来用 B 样条逼近未知系数函数,然后运用带有 l0惩罚项的最小二乘估计法来估计周期、周期项以及外生变量的变系数;第 3 章同时给出本文模型的渐近性质,包括周期估计b 的相合性、周期序列估计 b m 以及外生变量系数函数估计 b 的渐近分布;第 4 章为模拟研究,在不同周期、协变量个数以及变系数函数下说明了本文方法可以很好地捕捉到外生性变量的影响并能精准的估计出周期;第 5 章为实证部分,将本文方法应用到实际数据;最后第 6 章为结论,总结本文方法的优缺点,并对未来工作进行了展望.2模型及估计方法对于周期项相互独立的时间序列数据,经典的时间序列分解方法将时间序列写成周期项 m()、时间趋势项 g()以及不规律因素 的加法模型,即Yt=g(tT)+m(t)+t,t=1,2,T,(1)其中,g()是光滑的趋势函数,m(t)为时间序列的周期项,t 是一个平稳的均值为 0 的随机误差项.由于绝大多数时间序列受到外界因素的影响与序列本身的周期波动是相互独立的.模型分解的构成主要受周期项以及协变量的干扰,故本文构造的模型形式为Yt=m(t)+(t)Tz(t)+t,t=1,2,T,(2)其中,m(t)为时间序列的周期项,(t)=(1(t),2(t),p(t)T,z(t)=(z1(t),z2(t),zp(t)T,zi为第 i 个外生性变量,i()该外生性变量的系数函数.对于模型(2)中的各项,作以下假定:首先,对于误差项 t,假定 E(t)=0,本文中并不假定 t 独立,它可以是短期相关的,如 t 可以是一阶自回归模型(AR(1).其次,zi为第 i 个外生性变量,i=1,2,p;若序列中存在趋势项,则可以令 z1(t)1,此时为了保证模型的可识别性,需要假定101(u)du=0,若序列中不存在趋势项,则不需要假定101(u)du=0.外生性变量的选择要求该变量不存在明显的周期性,并会单方向对 Y 有影响.()为外生性变量的变系数函数.本文定义()为单位化时间 t/T 的函数,这样随着样本量 T 的增加,对于第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用53系数函数来说,每个小区间内的样本点也越来越多,因而可以更易获得系数函数估计的相合性、渐进正态性、有效性等.这种单位化时间的定义在很多文献中十分常见.最后,本文设置 m(t)tZ是周期为 0的周期序列,而不是周期函数,其中 0是未知整数,是该序列的最小周期,这个设定与文献 17 相同,原因在文献 17 有详细给出:对于等距观测的时间序列,在每个时间点 t,会存在很多个函数值为 m(t)的函数,而且如果周期项为函数,那么0可以不是最小周期,这样会导致周期项 m(t)和最小周期 0均不能识别,所以本文中假定周期项是一个周期序列.注意到周期项 m 是定义在真实时间 t 上的序列,这样方便我们研究周期序列估计的渐近性质,原因在于随着 T 增加,整个周期序列中与 m(t0)相等的样本点也会增加,可以获得更多关于 m(t0)的信息.令 =(1,2,0)T=(m(1),m(2),m(0)T,则 m(t)=0s=1sIs(t),其中 I()为示性函数,Is(t)=I(t=k0+s),其中 k 为非负整数.根据上述讨论,将模型(2)可改写为如下模式:Yt=pi=1i(tT)zi(t)+0s=1sIs(t)+t,t=1,2,T.(3)模型(3)用矩阵的形式表示为Y=Z1 1+Z2 2+Zp p+X0+,(4)其中 Y=(Y1,Y2,YT)T为我们观察到的时间序列的向量,也是我们重点想要分解与估计的向量;Zi i=(zi(1)i(1/T),zi(2)i(2/T),zi(T)i(T/T)T,X0=(I0,I0,)T是 T 0的设计矩阵,其中 I0是 00的单位矩阵;而 =(1,2,T)T是误差向量.如果该时间序列的周期已知,即 0已知,则该模型可退化为一个部分线性变系数模型,我们很容易就可以得到周期序列 m 以及变系数函数()的估计,因此接下来我们需要先对未知周期 0进行估计.下面我们主要讨论如何估计未知周期 0、周期序列 m 以及变系数函数().前面提到文献17 运用惩罚最小二乘法和局部线性核估计法对存在光滑趋势和周期项的经典时间序列分解模型即模型(1)进行了研究.但文献 17 的方法在估计周期时直接忽略了趋势项,忽略趋势项会增大估计量的方差,当趋势变化较大时,周期估计量不再准确,所以文献 17 的方法不适用于趋势变化较大的序列.原因在于直接忽略趋势项意味着我们将模型直接看作 Yt=m(t)+t,将 +g 作为新的误差项,当趋势变化较小时,与原始误差 相差较小,但如果趋势 g 的波动变化较大,那么忽略趋势项后的新误差项 会占据主导地位,从而导致周期项 m 不再明显,难以选择出正确的周期.而本文研究的模型(2)不仅可以包含趋势项还可以包含变系数外生性变量,本文在估计未知周期时没有忽略趋势项和其他外生变量的影响,即没有忽略(t)Tz(t),而是采用 B 样条逼近未知的变系数函数,B 样条在估计未知函数方面是非常实用的工具,被广泛的应用于变系数模型、可加模型等模型估计.54应用概率统计第 40 卷2.1变系数函数的 B 样条估计设 x=(x1,x2,xN),a x1 x2 xN b 为区间 a,b 上的 N 个不同的节点,称 x1=x2=xm+1=a,xm+2=x1,xm+N+1=xN,xm+N+2=x2(m+1)+N=b 为一组扩展的节点.令 Bi(t)=(1)m+1(xi+m+1 xi)xi,xi+1,xi+m+1(t x)m+,其中,(t x)m+=(t x)mItx,xi,xi+1,xi+m+1(t x)m+为函数(t x)m+的均差,同时,由上式定义的 Bi(t),i=1,2,m+N+1 线性无关,为样条函数空间 S(m,x)的一组基,称为 m 次规范化的 B 样条函数基.根据文献 25 递归公式,可得到规范化的 B 样条基函数的递推定义为B0i(t)=1,t xi,xi+1);0,其他.以及Bmi(t)=t xixi+m xiBm1i(t)+xj+m+1 txj+m+1 xj+1Bm1i+1(t).为方便,上式中上标 m 忽略,记 B 样条基向量为B(t)=(B1(t),B2(t),BNT(t),NT=m+N+1.对于任意的 s(t)S(m,x),Rm+N+1存在,使得 s(t)=BT(t),s(t)称为 m 次的 B 样条函数.本文中的时间序列,当周期已知,则可退化为部分线性变系数函数模型.对于第 l个外生性变量的变系数函数,定义Bl(ut)=(B1(ut),B2(ut),BNT(ut)T,ut=tT,t=1,2,T,Bl是基函数在所有单位化时间点 ut=t/T,t=1,2,T 取值组成的 B 样条矩阵,则第 l个变系数函数 l(u),l=1,2,p 的估计为 b l(ut)=BTl(ut)b l,其中 l,l=1,2,p 是NT 1 的向量,称 b l(ut)是 l(ut)的最小二乘 B 样条估计.对于本文模型,记D=z1(1)BT1(u1)z2(1)BT2(u1)zp(1)BTp(u1)z1(2)BT1(u2)z2(2)BT2(u2)zp(2)BTp(u2).z1(T)BT1(uT)z2(T)BT2(uT)zp(T)BTp(uT),T=(T1,T2,Tp),则式(4)用矩阵表示为Y D+X0+=W0+,(5)其中,W0=(X0,D),T=(T,T).当使用 B 样条时,B 样条基函数的个数与节点数和基函数的阶数有关,严格来说,基函数个数 NT等于内部节点数 N 与基函数的阶数(m+1)之和.本文中系数函数的 B 样条估计均使用相同的基函数阶数和相同的节点个数,这个假定在变系数模型中也是比较常见的(参见文献 26).第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用552.2未知周期 0的估计周期性数据最小二乘法的惩罚项改进周期 0的取值,直接决定时序模型的变化,故对于其取值可以说是模型的选择.那么一个自然的想法是借助最小二乘回归法,通过最小化残差平方和得到周期估计.具体来说,对于每个可能的周期 1,T(其中 T为所有周期可选值的上界,T可以随着样本量增大而增大),对模型(5)做最小二乘回归,得到 的最小二乘估计:b=(WTW)1WTY.针对参数估计的结果,我们得到残差平方和 RSS():RSS()=Y Wb 2,其中 表示向量的 l2范数.但众所周知,最小二乘回归中,残差平方和往往会随着解释变量的增加而减小.由于对于本文模型,解释变量会与 有关,越大,解释变量的个数越多,故最小化残差平方和往往会高估周期.由于真实周期 0的倍数 C0也是序列 m 的周期,而 C0对应的残差平方和小于真实周期 0的,所以普通最小二乘法得到的周期往往是真实周期的倍数.为了避免这个问题,对较大的周期增加一个较大的惩罚项.这种想法与贝叶斯信息准则(BIC)以及Akaike 信息准则(AIC)类似,AIC 和 BIC 准则是模型拟合度和模型复杂度之间找到一个平衡,我们知道:AIC=2k 2ln(bL),BIC=ln(n)k 2ln(bL),其中 k 是模型的有效参数个数,bL 是模型的似然概率,在最小二乘回归模型中,常用的是如下形式的 AIC、BIC:AIC=2k+nln(RSS),BIC=ln(n)k+nln(RSS),文献 17 用到的准则函数为eQ(,T)=RSS()+T,本文在选择周期时考虑到了外生项,必须对惩罚项加以改变,结合文献 27 半参数部分线性模型下的模型选择方法,本文采用的惩罚项形式为 T(+1)NT,其中 NT是 B 样条基函数的个数,如前一节所述,本文假定外生项变系数基函数个数均相同,T是惩罚参数,这个惩罚项可以作为 l0惩罚(l0惩罚被广泛应用于模型选择等问题的处理).为了保证周期估计的相合性,我们选择如下形式的惩罚参数:T=2T,其中 2=E(2t),T是一个发散的序列需要满足一定条件保证周期估计的相合性(惩罚参数 T满足的条件由定理 6 给出).那么针对周期性时间序列,我们最小化的准则函数 Q 有以下形式:Q(,T)=RSS()+(+1)NTT.最小化该准则函数得到b,b=argmin166TQ(,T).(6)56应用概率统计第 40 卷2.3未知周期序列 m 和变系数函数 的估计通过最小化 Q(,T)得到周期的估计b 后,我们可以很容易地得到周期项 m 和变系数函数 的估计,此时模型变成一个普通的部分线性变系数回归模型,我们对其作最小二乘回归即可得到系数 b的估计 b b:b b=(WTbWb)1WTbY.回想前面设定的 =(1,2,0)T=(m(1),m(2),m(0)T,因此,相应的周期序列 m 的估计为b=(b m(1),b m(2),b m(b)T,并且 b m(s+kb)=b m(s)对于任意s=1,2,b 和 k 都成立.同时,我们也可以得到第 l 个变系数函数 l()的估计 b l(ut)=BTl(ut)b l,l=1,2,p,记 Il=(0NT,INT,0NT,),也就是 Il表示只有第 l 块等于单位矩阵 INT其余块等于 0 矩阵的矩阵,则 l=Il,b l(u)=BTl(u)Ilb.在运用 B 样条方法时,我们需要决定基函数的个数以及节点位置,而基函数的个数与节点数以及基函数的阶数有关.一般来说,基函数的阶数对未知函数的估计影响会比较小,通常我们选择三阶样条或四阶样条.而节点的个数对函数的估计影响较大,需要进行决策.在 B 样条估计中,节点个数越多,窗宽 h 越小,拟合值与真实值之间的偏差越小,方差越大;节点个数越少,窗宽 h 越大,拟合值与真实值之间的偏差越大,方差越小.理论上可以通过最小化方差与偏差两者之和来求得最优节点数,常用的选择方法有广义交叉验证(GCV)、K 折交叉验证(KCV)、BIC 等.节点位置的选择对未知函数的估计效果影响也比较小,但对某些分布不均匀的数据影响可能比较大.通常来说,节点放置方法有两种,一种是等间距,本文采用该放置方法;另一种是利用样本分位数,在样本比较多的区间放置更多的点.本文在对周期估计时,估计的结果不应对节点数以及基函数的个数太过敏感,同时为了计算方便,在进行周期估计时,我们选择内部节点个数为 N T1/5,这是由于通过GCV、KCV、BIC 等方法得到的内部节点数大致在这个附近.在得到周期估计之后再对节点个数利用节点选择方法进行更为精确的计算,通过节点选择方法选出的内部节点一般也在 N T1/5的附近.2.42的简单估计注意到,惩罚参数 T=2T含有一个未知项 2,因此我们需要用估计值 e 2代替2.前面我们讨论过最小化 RSS()来估计周期得到的估计结果通常为真实周期 0的倍数,而不是 0.但在估计 2时,我们可以将e=argmin166TRSS()作为 0的初始估计,这样做可行的原因在于 c0也是这个周期序列的一个周期,并且文献 17 证明了P(e=k0)1.得到周期的初始估计后,我们就可以相应地得到周期序列和变系数函数的估计 e m(t)和 e(t/T),则可以得到 2的一个很直接的估计,e 2=1TTt=1e 2t=1TTt=1Yt e m(t)pi=1e i(tT)zi(t)2.第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用573渐近性质3.1基本假设为了得到估计的渐近性质,我们给出如下五个假设:条件 1误差项 t 是强混合的,存在某个正常数 C 和 a 1 使得(k)6 Cak,其中(k)是混合系数.条件 2存在某个常数 0 C 0 使得 E(|t|4+)6 C.条件 3每个变系数函数 i()在 0,1 上 2 阶连续可导,i=1,2,p.条件 4相邻两个节点的距离 h T1/5,即内部节点数 N T1/5.条件 5z=(z1,z2,zp)T,t 是随机的,t 的边际密度函数是连续的,Q(z|t)为给定 z,t 的条件分布函数.E(zzT|t)=G(t),gij(t)=zizjdQ(z|t)为 G(t)的(i,j)元素,存在常数 m3,M2,M3,使得gij(t)6 M2,m36 G(t)min6 G(t)max6 M3,(7)其中,G(t)min,G(t)max分别为 G(t)的最小、最大的特征根.条件 1条件 3 在文献 17 也有假定,类似地,本文中也不假设误差项 t 平稳,但我们需要对它的相关结构作出一定的要求,注意到我们在估计未知周期并未考虑误差项的相关结构,所以在这里我们要求误差项短期相关,且是强混合的.条件 3 是变系数函数的光滑条件,是为了保证变系数函数估计的渐近正态性.条件 4 是对 B 样条节点个数的假设,这个条件在 B 样条中也十分常见.条件 5 是对回归变量 z,t 的限制,当 gij(t)为 a,b 上的连续函数时,式(7)成立.3.2周期估计b 的相合性下面的定理给出了周期估计b 的渐近性质.定理 6若条件 1条件 5 成立,并且存在某个正常数 C 0使得 T6 CT2/5成立,惩罚参数 T满足(lnT)3/2T/(TNT)0,TNT/T 0,那么bP 0,也就是b 是周期 0的相合估计.在上面的定理中,T可以不是固定常数,可以随着 T 的增大以一定的速度增大,但T发散到正无穷的速度越快,我们对惩罚项 T的要求就越高,但是 T发散到正无穷的速度不能比 T2/5快.3.3周期序列估计 b m 及变系数函数的估计 b 的渐近性质定理 6 给出了周期估计b 的相合性,前面提到过如果模型的周期是已知的,那么我们的模型就简化成普通的半参数部分线性变系数模型,而部分线性变系数模型的理论已经很58应用概率统计第 40 卷成熟了.接下来我们看一下周期未知时模型中线性部分也就是周期序列估计 b m 以及变系数函数的估计 b 的收敛速度和渐近分布.首先为了简化陈述,对于每个时间点 t,定义 t0=t0t/0,Kt0,T=1+(T t0)/,我们定义 Vt0,T=xTt0+(k1)0111xt0+(k1)0,其中 xTt0+(k1)0是 X0的行向量,=(ij)00,ij=Cov(xi e xi(z,t),xj e xj(z,t),1 6 i,j 6 0,e xi(z,t)=E(xi|z,t),1=T1Ti,j=1Cov(xi e xi(z,t)i,(xj e xj(z,t)j).定理 7若条件1条件5成立,进一步假设Vt0,T的极限存在,即Vt0=limTVt0,T.那么对每一个时间 t,t=1,2,T,我们有T b m(t)m(t)d N(0,Vt0).可以看到,当周期未知时,周期序列的估计 b m的渐近分布与周期已知时相同,并且渐近方差Vt0也相同,也就意味着周期估计中的误差在周期序列估计 b m的渐近分布中可以被忽略,这一点可以从附录中定理的证明中看到.并且部分线性变系数模型中的线性部分与普通线性模型的渐近分布的收敛速度相同.同时,注意到我们假设了渐近方差的存在性,且这一假设在很多情况下都是成立的.若t是平稳的,那么长期方差很容易计算,我们可以用样本自协方差并且结合周期图或者谱密度来计算长期方差.也有大量的文献研究如何计算一般情况下的长期方差,比如文献28和文献29,主要是利用样本自协方差的核加权.最后我们给出变系数函数的 B 样条估计 b()的一致收敛速度和渐近分布.定理 8若条件 1条件 5 成立,我们有(i)supu0,1|b j(u)j(u)|=Op(0(NT)(N/T+N2),l=1,2,p,其中 0(NT)是满足 supu0,1D(u)6 0(NT)且(20(NT)NT)/T 0 的常数数列.(ii)对于每一个 u (0,1),假设 Vu=limTVu,T存在,如果低平滑条件成立,也就是内部节点数满足 T1/5N1 0,那么我们有T/N(b l(u)l(u)d N(0,Vu),l=1,2,p,其中 Vu,T=N1BTl(u)Il121ITlBl(u),=T1DTD,2=T1Ti,j=1ED(ui)DT(uj)ij,Il的定义与第 2.3 节中 Il=(0NT,INT,0NT,)相同,表示只有第 l 块等于单位矩阵 INT其余块等于 0 矩阵的矩阵.可以看到,变系数函数估计 b()的渐近分布与周期序列 m 已知时得到的估计 e 相同,也就意味着周期序列估计中的误差在变系数函数估计 b 的渐近分布中可以被忽略.4模拟研究这一章我们利用蒙特卡洛模拟去探索本文方法的有限样本性质.为了方便比较,文本的第一个模拟的假设参考了文献 17 模拟的假设,为了体现本文改进方法的优越性,在第二个模拟中,我们假设存在三个外生变量,并且外生项的系数函数变化更大.在两种不同模拟假设下,比较本文方法与基于文献 17 直接忽略外生变量的方法所得的周期估计的结果.第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用59模拟一的设定如下:周期序列 m 是一个正弦序列,真实周期为 0,m(t)=sin(20t+32),t=1,2,T.设定模型中只有一个外生性变量 z1,该变量来自带有截距项的一阶自回归模型:z1(t)=1+0.6z1(t 1)+et,t=1,2,T,其中 et N(0,0.25).该外生性变量的系数函数为:1(t)=2(t/T)2,t=1,2,T.误差项 t来自一阶自回归模型 AR(1):t=0.45t1+t,其中 t取自 N(0,2).我们将通过改变 2的值来改变信噪比,并且比较三种方差 Var(t)=2=0.25,0.5,0.75 和不同样本量下的周期估计.注意到,为了保证 Var(t)=2=0.25,0.5,0.75,则需要 2=0.2,0.4,0.6.对于每个样本量 T和方差 2,我们考虑的备选周期范围为 1,T/2,独立重复模拟 N=1000 次.模拟二设定三个外生性变量,其中假设 z1(t)来自均匀分布,z2(t)来自正态分布,z3(t)来自 AR(1):z1(t)v unif(1,1),z2(t)v N(1,1),z3(t)=1+0.6z3(t 1)+et,t=1,2,T,其中 et N(0,0.25).这三个外生性变量的变系数函数分别为1(t)=1+(tT)+2(tT)2,2(t)=2cos(2tT),3(t)=10(tT)3.同时,模拟中的惩罚参数 T=e 2T,其中 T=lnT,e 2是 2的估计,在前文中我们已经给出了其估计方法.4.1模拟一在模拟一中我们假设真实周期 0=60,样本量 T=200,250,500.由于周期估计是否准确会直接影响周期序列的估计效果,也会影响变系数的估计效果,进而会影响时间序列的预测准确性.因此,得到一个较为准确的周期估计是至关重要的.在本文模拟中,我们重点比较两种方法的周期估计.图 1 描绘了样本量在 200、350、500 下两种方法的周期估计结果,每幅图从上到下依分别次是在不同方差 2通过本文方法和基于文献 17 方法得到的周期估计b 的频数分布直方图,由于文献 17 的模型中仅包含时间趋势和周期项,不包含外生变量,此处模拟中所采用的文献 17 的方法是在估计周期时直接忽略外生变量.为了更清晰地比较模拟结果,表 1 展示了不同方法下的周期估计的准确率,即 1000 次模拟中周期估计值估计值b 恰好等于真实周期 0的比率,和 1000 次模拟得到的周期估计值的平均值.60应用概率统计第 40 卷0600100200300仇0600100200300仇(a)T=200,2=0.250600100200仇0600100200仇(b)T=200,2=0.5060050100150仇0600100200300仇(c)T=200,2=0.750600250500750仇0600200400仇(d)T=350,2=0.250600200400仇0600200400仇(e)T=350,2=0.50600200400仇0600100200300仇(f)T=350,2=0.7506005001000仇0600250500750仇(g)T=500,2=0.250600250500750仇0600200400600仇(h)T=500,2=0.50600200400600仇0600200400600仇(i)T=500,2=0.75图 1本文方法和文献 17 方法周期估计的直方图比较(每幅图的上图是本文方法,下图是文献 17 方法)表 1本文方法和文献 17 方法在不同样本量 T 和方差 2下的周期估计准确率及均值比较T2本文方法V&L2000.2533.9%(59.94)0.4%(60.71)0.520.3%(60.03)1.8%(56.58)0.7515.7%(59.51)2.4%(63.73)3500.2584.7%(59.96)47.6%(60.96)0.550.5%(60.55)38.1%(56.22)0.7540.5%(60.02)33.1%(54.76)5000.2599.3%(60.00)77.4%(59.09)0.584.3%(59.98)51.7%(60.80)0.7568.9%(59.93)38.4%(53.51)注:V&L 表示文献 17 的方法,括号里为 1000 次模拟所得的周期估计的平均值.首先从表 1 中来看,在不同的样本量下,本文方法的准确度均高于基于文献 17 的方法的准确度,从周期估计的均值来看,本文的方法的结果更加切近真实周期.当样本量第 1 期方学莉,王守霞:变系数的周期性时间序列模型及其应用61T=200 时,尽管此时周期估计准确率较低,但从直方图 1 中可以看到,大部分的b 聚集在真实周期附近,同时从表 1 中可以看到 1000 次模拟的周期估计均值也基本等于真实周期,这也说明即使在样本量很少的情况下,本文方法也可以得到一个较好的周期估计.从直方图 1 来看,通过文献 17 直接忽略外生变量的方法得到的周期估计值的分布更加分散,而通过本文方法得到的周期估计值的分布更集中在真实周期附近,随着样本量的增加,越来越多的b 聚集在真实周期附近.从表 1 中也可以看到,周期估计的准确率随着样本量的增加逐渐提高并趋近于 1.当样本量 T=500 时,几乎所有的估计值都等于真实周期,且当误差方差较小时,在 1000 次模拟下,准确率高达 99.3%.从直方图 1 中可以看出,随着误差方差的增大,周期估计b 的分布也越来越分散.原因在于随着噪声的增大,得到一个较为准确的周期估计也越来越困难.当样本量较少且误差的方差较大时,本文的方法表现的仍比文献 17 方法更好.0.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(a)T=200,2=0.250.00.20.40.60.81.0t/T0.50.00.51.01.52.02.511195%PCB(b)T=200,2=0.50.00.20.40.60.81.0t/T0.50.00.51.01.52.02.511195%PCB(c)T=200,2=0.750.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(d)T=350,2=0.250.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(e)T=350,2=0.50.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(f)T=350,2=0.750.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(g)T=500,2=0.250.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(h)T=500,2=0.50.00.20.40.60.81.0t/T0.00.51.01.52.011195%PCB(i)T=500,2=0.75图 2变系数函数的估计在得到周期估计后,我们可以得到周期序列以及变系数函数的估计 b 及其 95%的逐点置信区间.图 2 展示了不同样本量和误差下变系数函数的估计结果,从上到下样本量为200、350 和 500,从左到右误差方差分别为 0.25、0.5 和 0.75.其中红色虚线是估计结果,62应用概率统计第 40 卷蓝色实线是真实的变系数函数,黑色点线是 95%置信区间.从图 2 我们可以看出,即使在样本量较少的情况下,本文基于 B 样条的估计法也能十分准确地估计未知变系数函数,同时,置信区间也会随着样本量增大或者误差的减小而变窄.注意到,变系数函数的估计 b 在周期估计相对准确时也更为准确,这也说明周期估计的准确性会影响变系数函数估计.模拟一中通过本文方法与文献 17 方法对比,发现本文方法相对来说更优,但由于模拟一中假设的变系数函数变化较小,故二者的差异也较小,为了进一步说明本文方法的优越性,在模拟二中,本文增加外生变量的个数,并增大变系数函数的变化从而增大外生性变量的影响程度,进一步比较两种方法.4.2模拟二由于第 5 章中我们研究的澳门入境旅游人数周期为 12,故模拟二的周期序列 m 的真实周期 0我们设置为 30 或者 12.表 2 展示了真实周期 0=30 时两种方法的周期估计的准确率,其中括号中的值代表每个方法 1000 次模拟得到的周期估计值的平均值.为了更清楚看清本文方法与文献 17方法的差异,图 3 给出了基于本文方法与文献 17 忽略外生变量的方法得到的周期估计值的频数分布直方图.表 2真实周期 0=30 时本文方法和文献 17 方法的周期估计准确率及均值比较T2本文方法V&L1000.2521.6%(30.47)0%(1.00)0.516.8%(30.78)0%(1.00)0.7513.0%(32.31)0%(1.00)1800.2589.1%(30.00)0%(1.00)0.567.4%(30.68)0%(1.00)0.7549.5%(32.85)0%(1.00)2500.25100%(30.00)0%(1.00)0.596.5%(29.97)0%(1.00)0.7590.8%(29.71)0%(1.00)注:V&L 表示文献 17 的方法,括号里为 1000 次模拟所得的周期估计的平均值.从表 2 中可以看到,当变系数函数变化较大时,基于文献 17 方法得到的周期估计的准确率基本为 0,而通过本文方法得到的周期估计准确率仍然较高,周期估计的均值也十分接近真实周期 30.从模拟直方图 3 的结果可以看到文献 17 方法的周期估计几乎都等于 1,也就是说在时间序列的变系数函数变化趋势比较大的情况下,文献 17 忽略趋势或者外生变量的方法基本失效