1、第二章,时间序列的预处理,本章结构,平稳性检验 纯随机性检验,2.1平稳性检验,特征统计量平稳时间序列的定义平稳时间序列的统计性质平稳时间序列的意义平稳性的检验,概率分布,概率分布的意义随机变量族的统计特性完全由它们的联合分布函数或联合密度函数决定 时间序列概率分布族的定义实际应用的局限性,特征统计量,均值 方差自协方差自相关系数,平稳时间序列的定义,严平稳严平稳是一种条件比较苛刻的平稳性定义，它认为只有当序列所有的统计性质都不会随着时间的推移而发生变化时，该序列才能被认为平稳。宽平稳宽平稳是使用序列的特征统计量来定义的一种平稳性。它认为序列的统计性质主要由它的低阶矩决定，所以只要保证序列低阶
2、矩平稳（二阶），就能保证序列的主要性质近似稳定。,平稳时间序列的统计定义,满足如下条件的序列称为严平稳序列满足如下条件的序列称为宽平稳序列,严平稳与宽平稳的关系,一般关系严平稳条件比宽平稳条件苛刻，通常情况下，严平稳（低阶矩存在）能推出宽平稳成立，而宽平稳序列不能反推严平稳成立特例不存在低阶矩的严平稳序列不满足宽平稳条件，例如服从柯西分布的严平稳序列就不是宽平稳序列当序列服从多元正态分布时，宽平稳可以推出严平稳,平稳时间序列的统计性质,常数均值 自协方差函数和自相关函数只依赖于时间的平移长度而与时间的起止点无关 延迟k自协方差函数 延迟k自相关系数,自相关系数的性质,规范性 对称性 非负定性
3、非唯一性,平稳时间序列的意义,时间序列数据结构的特殊性可列多个随机变量，而每个变量只有一个样本观察值平稳性的重大意义极大地减少了随机变量的个数，并增加了待估变量的样本容量极大地简化了时序分析的难度，同时也提高了对特征统计量的估计精度,平稳性的检验（图检验方法）,时序图检验 根据平稳时间序列均值、方差为常数的性质，平稳序列的时序图应该显示出该序列始终在一个常数值附近随机波动，而且波动的范围有界、无明显趋势及周期特征自相关图检验 平稳序列通常具有短期相关性。该性质用自相关系数来描述就是随着延迟期数的增加，平稳序列的自相关系数会很快地衰减向零,例题,例2.1检验1964年1999年中国纱年产量序列的
4、平稳性例2.2检验1962年1月1975年12月平均每头奶牛月产奶量序列的平稳性例2.3检验1949年1998年北京市每年最高气温序列的平稳性,例2.1时序图,例2.1自相关图,例2.2时序图,例2.2 自相关图,例2.3时序图,例2.3自相关图,2.2 纯随机性检验,纯随机序列的定义纯随机性的性质纯随机性检验,纯随机序列的定义,纯随机序列也称为白噪声序列，它满足如下两条性质,标准正态白噪声序列时序图,白噪声序列的性质,纯随机性 各序列值之间没有任何相关关系，即为 “没有记忆”的序列 方差齐性 根据马尔可夫定理，只有方差齐性假定成立时，用最小二乘法得到的未知参数估计值才是准确的、有效的,纯随机
5、性检验,检验原理假设条件检验统计量 判别原则,Barlett定理,如果一个时间序列是纯随机的，得到一个观察期数为 的观察序列，那么该序列的延迟非零期的样本自相关系数将近似服从均值为零，方差为序列观察期数倒数的正态分布,假设条件,原假设：延迟期数小于或等于 期的序列值之间相互独立备择假设：延迟期数小于或等于 期的序列值之间有相关性,检验统计量,Q统计量 LB统计量,判别原则,拒绝原假设当检验统计量大于 分位点，或该统计量的P值小于 时，则可以以 的置信水平拒绝原假设，认为该序列为非白噪声序列接受原假设当检验统计量小于 分位点，或该统计量的P值大于 时，则认为在 的置信水平下无法拒绝原假设，即不能
6、显著拒绝序列为纯随机序列的假定,例2.4：标准正态白噪声序列纯随机性检验,样本自相关图,检验结果,由于P值显著大于显著性水平 ，所以该序列不能拒绝纯随机的原假设。,例2.5,对1950年1998年北京市城乡居民定期储蓄所占比例序列的平稳性与纯随机性进行检验,例2.5时序图,例2.5自相关图,例2.5白噪声检验结果,非平稳时间序列平稳化处理,Box-Cox变换f(x,)=(x1)/ if 0 = log(x) if =0差分: B为后移算子BXt=Xt-1 Bd Xt=Xt-d d阶差分 (1-B)d Xt,Box-Cox 变换,library(MASS)library(car)library(
7、pander)l - lm(Volume log(Height) + log(Girth), data = trees) #建立线性模型qqPlot(l) #残差的QQ图，不大符合正态分布,boxcox(Volume log(Height) + log(Girth), data = trees) #找lambda,boxcox(Volume log(Height) + log(Girth), data = trees, lambda = seq(-0.08, 0, length = 10),kk=boxcox(Volume log(Height) + log(Girth), data = tr
8、ees, lambda = seq(-0.08,0, length = 10)kk$xwhich(kk$y=max(kk$y) # -0.06707071,# 缩小寻找的范围,大约是-0.067volume - (trees$Volume(-0.67) - 1)/(-0.067) #变换trees.t - cbind(trees, volume) #重新拟合模型l.t |t|) # - - - - -# *(Intercept)* -0.1368 1.517 -0.09022 0.9288 # # *log(Height)* 1.745 0.3877 4.502 0.000108 # # *l
9、og(Girth)* 2.358 0.1422 16.58 5.213e-16 # -# # Table: Fitting linear model: volume log(Height) + log(Girth),pander(anova(l.t)# # -# Df Sum Sq Mean Sq F value Pr(F) # - - - - - -# *log(Height)* 1 5.854 5.854 245.8 2.154e-15# # *log(Girth)* 1 6.546 6.546 274.9 5.213e-16# # *Residuals* 28 0.6669 0.0238
10、2 # -# # Table: Analysis of Variance Table,2.forecast包的BoxCox.lambda和BoxCox,BoxCox.lambda这个函数用于数值向量或时间序列，可以得到lambda的估计精确值。library(forecast),BoxCox.lambda(trees$Volume, method = loglik) #算出来的结果和boxcox有点差异# 1 -0.05volume.f - BoxCox(trees$Volume, lambda = -0.05)trees.f - cbind(trees, volume.f) #重新拟合模型l
11、.f |t|) # - - - - -# *(Intercept)* -5.454 0.6731 -8.103 8.013e-09 # # *log(Height)* 0.965 0.172 5.609 5.269e-06 # # *log(Girth)* 1.678 0.06313 26.58 2.073e-21 # -# # Table: Fitting linear model: volume.f log(Height) + log(Girth),pander(anova(l.f)# # -# Df Sum Sq Mean Sq F value Pr(F) # - - - - - -#
12、*log(Height)* 1 2.534 2.534 540.1 7.646e-20# # *log(Girth)* 1 3.315 3.315 706.7 2.073e-21# # *Residuals* 28 0.1314 0.004691 # -# # Table: Analysis of Variance Table,qqPlot(l.f),3.1 方法性工具,差分运算延迟算子线性差分方程,差分运算,一阶差分 阶差分 步差分,延迟算子,延迟算子类似于一个时间指针，当前序列值乘以一个延迟算子，就相当于把当前序列值的时间向过去拨了一个时刻 记B为延迟算子，有,延迟算子的性质,，其中,用延
13、迟算子表示差分运算,阶差分 步差分,线性趋势采用1阶差分 yt=(1-B) Xt=Xt-Xt-1抛物线趋势采用2阶差分 yt=(1-B) 2Xt=Xt+Xt-2-2Xt-1,季节差分,yt=(1-Bs) Xt=Xt-Xt-s s 为周期 月度数据 s=12 季度数据s=4,对数变换与差分运算结合,金融经济数据 工业总产值 yt=(1-B)2 log(Xt) 股票收盘价 yt=(1-B) log(Xt)=log(Xt/Xt-1) 称为股票收益率，采用泰勒展开 约为 yt=(Xt-Xt-1)/Xt-1,线性差分方程,线性差分方程齐次线性差分方程,齐次线性差分方程的解,特征方程特征方程的根称为特征根
14、，记作齐次线性差分方程的通解不相等实数根场合有相等实根场合复根场合,非齐次线性差分方程的解,非齐次线性差分方程的特解使得非齐次线性差分方程成立的任意一个解非齐次线性差分方程的通解齐次线性差分方程的通解和非齐次线性差分方程的特解之和,异常值处理,数值检验 缺损值的补足模型分析,数值检验,if mean(X1:t)-k sd(X1:t)Xt+1 mean(X1:t)+k sd(X1:t) k=6, #6 sigma, Xt+1 正常 否则 为离群点 如果为离群点 X*t+1=2Xt-Xt-1 #线性外推,缺损值的补足,平滑法 插值估算法,模型分析,估计模型 残差分析 如果残差分析不能通过，则替换或
15、加入哑变量,例子： 人口自然增长率收集,时间跨度1950-2005时间间隔相同， 年度数据 ，不允许出现半年，2年等间隔数据 统一计算方法：（年初+年末）/2，序列范围一致性：虽然1997香港回归，并不能计算进入。 统一指标口径：出生后但随后死亡，列入出生也例如死亡。,趋势因子和季节因子的估计与去除,时间序列经典分解 yt=mt+St+at mt为趋势因子， St 为季节因子 at为随机因子,无季节因子，仅趋势因子,yt=mt+at，去除mt方法：最小二乘法滑动平均法差分方法,最小二乘法,用函数参数族拟合mt 例如 mt=a+bt+ct2 极小化 sum(yt- a-bt-ct2)2得到 a,
16、b,c的估计,读取人口,y=scan()3929214530848372398819638453128607021706335323191876314433213855837150189209629797667621216892228496103021537123202624132164569151325798179323175203302031226545805,读取时间,x=scan()17901800181018201830184018501860187018801890190019101920193019401950196019701980,fit=lm(yx+I(x2)summary
17、(fit)plot(x,y)lines(fit$fitted.valuesx)plot(x,fit$residuals,type=l),滑动平均（低通滤波，线性滤波器）,双边平滑 mt 在t-q,t+q近似线性 m1t =sum(yt-q:t+q)/(2q+1) =mean(yt-q:t+q) =mean(mt-q:t+q)+mean(at-q:t+q) mean(mt-q:t+q) mtmean(at-q:t+q) 快速震荡 mt缓慢变化趋势,z=scan()4737511750913468432038253673369437083333336736143362365539634405459
18、55045570057165138501053536074503156485506423048273885,x1=seq(1951,1980,by=1)n=length(z);m=rep(0,n); q=2for(t in (q+1):(n-q)mt=mean(z(t-q):(t+q),SMOOTH程序， 左边,#0.1a0.3a=0.2 for(t in 1:q)m1=0 for(j in 0: (n-t)m1=m1+a*(1-a)j*zt+j mt=m1,右边,for(t in (n-q+1):n)m1=0 for(j in 0: (t-1)m1=m1+a*(1-a)j*zt-j mt=m
19、1,plot(x1,z,type=l)lines(x1,m,col=red)w=z-mplot(x1,w,type=l),递归关系,tn-q mt=sum_j=0t-1 a(1-a)j yt-j,差分,k 阶差分作用k次多项式趋势项，其结果为常数(归纳法可证)人口数据2阶差分z=diff(y,differences=2)plot(x-c(1:2),z,type=l),趋势项和季节项同时存在,yt=mt+St+at 周期d; St+d=St sum(S1:d)=0例月度数据 d=12 ，n年 第j 年第k 个月数据记为yk,j,月事故死亡率,年,月,X=scan()900781068928913
20、710017108261137110744971399389161892777506981803884228714951210120982387439129871086808162730681247870938795561009396208285843381608034771774617776792586348945100789179803784887874864777926957772681068890929910625930283148850826587967836689277918129911594341048498279110907086339240,方法一、小趋势法,设第j年趋势相同
21、为Mj,由于 mean(S1:d)=0 则Mj=mean(y1:d,j) Sk=mean(yk,1:n-M1:n),XM=matrix(X,12,6)ts.plot(X)M=apply(XM,2,mean)d=12S=rep(0,12)for(i in 1:d)Si=mean(XMi,-M),n=6;Y1=Xfor(i in 1:n)for(j in 1:d)Y1j+(i-1)*d=Xj+(i-1)*d-Mi-Sjts.plot(Y1),方法二 、滑动平均,第一步 用滑动平均估计去趋势项 1.周期d=2q偶数时 mt=(sum(yt-q+1:t+q-1) +1/2(yt-q+yt+q)/d 2
22、.周期为d=2q+1奇数 mt=mean(yt-q:t+q),q=d/2m=rep(0,12*6)for(t in (q+1):(12*6-q)mt=(sum(X(t-q+1):(t+q-1)+Xt-q/2+Xt+q/2)/d,第二步去季节项,wk为yk+j*d-mk+j*d, q=k+jd=n-q的平均值,这里1=k=d 因为sum(wk)不等于0 所以Sk=wk-mean(w1:d),w=rep(0,d);N=12*6for(k in 1:d)w1=0;j=ceiling(q-k)/d);j1=0while(k+j*d=N-q)w1=w1+Xk+j*d-mk+j*dj=j+1;j1=j1+
23、1wk=w1/j1,a=0.2;N=length(X) for(t in 1:q)m1=0 for(j in 0: (N-t)m1=m1+a*(1-a)j*Xt+j mt=m1,for(t in (N-q+1):N)m1=0 for(j in 0: (t-1)m1=m1+a*(1-a)j*Xt-j mt=m1,Y2=Xfor(i in 1:6)for(j in 1:d)Y2j+(i-1)*d=Xj+(i-1)*d-mj+(i-1)*d-Sjx1=c(1:n)ts.plot(Y2)lines(x1,Y1,col=red),方法三、 差分,d步差分去季节yt-yt-d 再采用多部差分去趋势,ts.plot(Y1)lines(Y2,col=red)Y3=c(rep(0,d),diff(X,lag=12)lines(Y3,col=blue)Y4=c(rep(0,d+1),diff(diff(X,lag=12)lines(Y4,col=green),acf(X) acf(Y1)acf(Y2)acf(Y3) acf(Y4),
