第二章時間序列的預處理

1、第二章,時間序列的預處理,本章結(jié)構(gòu),平穩(wěn)性檢驗 純隨機性檢驗,2.1平穩(wěn)性檢驗,特征統(tǒng)計量平穩(wěn)時間序列的定義平穩(wěn)時間序列的統(tǒng)計性質(zhì)平穩(wěn)時間序列的意義平穩(wěn)性的檢驗,概率分布,概率分布的意義隨機變量族的統(tǒng)計特性完全由它們的聯(lián)合分布函數(shù)或聯(lián)合密度函數(shù)決定 時間序列概率分布族的定義實際應用的局限性,,特征統(tǒng)計量,均值 方差自協(xié)方差自相關(guān)系數(shù),,,,,平穩(wěn)時間序列的定義,嚴平穩(wěn)嚴平穩(wěn)是一種條件比較苛刻的平穩(wěn)性定

2、義，它認為只有當序列所有的統(tǒng)計性質(zhì)都不會隨著時間的推移而發(fā)生變化時，該序列才能被認為平穩(wěn)。寬平穩(wěn)寬平穩(wěn)是使用序列的特征統(tǒng)計量來定義的一種平穩(wěn)性。它認為序列的統(tǒng)計性質(zhì)主要由它的低階矩決定，所以只要保證序列低階矩平穩(wěn)（二階），就能保證序列的主要性質(zhì)近似穩(wěn)定。,平穩(wěn)時間序列的統(tǒng)計定義,滿足如下條件的序列稱為嚴平穩(wěn)序列滿足如下條件的序列稱為寬平穩(wěn)序列,,,,嚴平穩(wěn)與寬平穩(wěn)的關(guān)系,一般關(guān)系嚴平穩(wěn)條件比寬平穩(wěn)條件苛刻，通常情況下，嚴平穩(wěn)

3、（低階矩存在）能推出寬平穩(wěn)成立，而寬平穩(wěn)序列不能反推嚴平穩(wěn)成立特例不存在低階矩的嚴平穩(wěn)序列不滿足寬平穩(wěn)條件，例如服從柯西分布的嚴平穩(wěn)序列就不是寬平穩(wěn)序列當序列服從多元正態(tài)分布時，寬平穩(wěn)可以推出嚴平穩(wěn),平穩(wěn)時間序列的統(tǒng)計性質(zhì),常數(shù)均值 自協(xié)方差函數(shù)和自相關(guān)函數(shù)只依賴于時間的平移長度而與時間的起止點無關(guān) 延遲k自協(xié)方差函數(shù) 延遲k自相關(guān)系數(shù),,,,自相關(guān)系數(shù)的性質(zhì),規(guī)范性 對稱性 非負定性 非唯一性,,,,平穩(wěn)時間序列的

4、意義,時間序列數(shù)據(jù)結(jié)構(gòu)的特殊性可列多個隨機變量，而每個變量只有一個樣本觀察值平穩(wěn)性的重大意義極大地減少了隨機變量的個數(shù)，并增加了待估變量的樣本容量極大地簡化了時序分析的難度，同時也提高了對特征統(tǒng)計量的估計精度,平穩(wěn)性的檢驗（圖檢驗方法）,時序圖檢驗 根據(jù)平穩(wěn)時間序列均值、方差為常數(shù)的性質(zhì)，平穩(wěn)序列的時序圖應該顯示出該序列始終在一個常數(shù)值附近隨機波動，而且波動的范圍有界、無明顯趨勢及周期特征自相關(guān)圖檢驗 平穩(wěn)序列通常具有短期

5、相關(guān)性。該性質(zhì)用自相關(guān)系數(shù)來描述就是隨著延遲期數(shù)的增加，平穩(wěn)序列的自相關(guān)系數(shù)會很快地衰減向零,例題,例2.1檢驗1964年——1999年中國紗年產(chǎn)量序列的平穩(wěn)性例2.2檢驗1962年1月——1975年12月平均每頭奶牛月產(chǎn)奶量序列的平穩(wěn)性例2.3檢驗1949年——1998年北京市每年最高氣溫序列的平穩(wěn)性,例2.1時序圖,,例2.1自相關(guān)圖,,例2.2時序圖,,例2.2 自相關(guān)圖,,例2.3時序圖,,例2.3自相關(guān)圖,,2.2

6、純隨機性檢驗,純隨機序列的定義純隨機性的性質(zhì)純隨機性檢驗,純隨機序列的定義,純隨機序列也稱為白噪聲序列，它滿足如下兩條性質(zhì),,標準正態(tài)白噪聲序列時序圖,,白噪聲序列的性質(zhì),純隨機性 各序列值之間沒有任何相關(guān)關(guān)系，即為 “沒有記憶”的序列 方差齊性 根據(jù)馬爾可夫定理，只有方差齊性假定成立時，用最小二乘法得到的未知參數(shù)估計值才是準確的、有效的,,,純隨機性檢驗,檢驗原理假設條件檢驗統(tǒng)計量 判別原則,,Barlett定理,

7、如果一個時間序列是純隨機的，得到一個觀察期數(shù)為 的觀察序列，那么該序列的延遲非零期的樣本自相關(guān)系數(shù)將近似服從均值為零，方差為序列觀察期數(shù)倒數(shù)的正態(tài)分布,假設條件,原假設：延遲期數(shù)小于或等于 期的序列值之間相互獨立備擇假設：延遲期數(shù)小于或等于 期的序列值之間有相關(guān)性,,,檢驗統(tǒng)計量,Q統(tǒng)計量 LB統(tǒng)計量,,,判別原則,拒絕原假設當檢驗統(tǒng)計量大于 分位點，或該統(tǒng)計量的P值小于 時，則可以以 的置信水平拒絕

8、原假設，認為該序列為非白噪聲序列接受原假設當檢驗統(tǒng)計量小于 分位點，或該統(tǒng)計量的P值大于 時，則認為在 的置信水平下無法拒絕原假設，即不能顯著拒絕序列為純隨機序列的假定,,,例2.4：標準正態(tài)白噪聲序列純隨機性檢驗,,樣本自相關(guān)圖,檢驗結(jié)果,由于P值顯著大于顯著性水平 ，所以該序列不能拒絕純隨機的原假設。,例2.5,對1950年——1998年北京市城鄉(xiāng)居民定期儲蓄所占比例序列的平穩(wěn)性與純隨機性進行檢驗,例2.5時

9、序圖,,例2.5自相關(guān)圖,,例2.5白噪聲檢驗結(jié)果,非平穩(wěn)時間序列平穩(wěn)化處理,Box-Cox變換f(x,λ)=(x^λ?1)/λ if λ ≠0 = log(x) if λ =0差分: B為后移算子BX[t]=X[t-1] B^d X[t]=X[t-d] d階差分 (1-B)^d X[t],Box-Cox 變換,library(MASS)library(car)library(pander)

10、l <- lm(Volume ~ log(Height) + log(Girth), data = trees) #建立線性模型qqPlot(l) #殘差的QQ圖，不大符合正態(tài)分布,,,boxcox(Volume ~ log(Height) + log(Girth), data = trees) #找lambda,,,boxcox(Volume ~ log(Height) + log(Girth), data = tr

11、ees, lambda = seq(-0.08,  0, length = 10)),,,kk=boxcox(Volume ~ log(Height) + log(Girth), data = trees, lambda = seq(-0.08,0, length = 10))kk$x[which(kk$y==max(kk$y))] # -0.06707071,,# 縮小尋找的范圍,大約是-0.067volume &

12、lt;- (trees$Volume^(-0.67) - 1)/(-0.067) #變換trees.t <- cbind(trees, volume) #重新擬合模型l.t <- lm(volume ~ log(Height) + log(Girth), data = trees.t) #建立線性模型qqPlot(l.t) #殘差可認為是正態(tài)了,,,,pander(l.t)## ## ------------

13、--------------------------------------------------## &nbsp; Estimate Std. Error t value Pr(>|t|) ## ----------------- ---------- ------------ --------- ----------## **(Intercept)** -0.1368 1.517 -0.09022 0.9288

14、 ## ## **log(Height)** 1.745 0.3877 4.502 0.000108 ## ## **log(Girth)** 2.358 0.1422 16.58 5.213e-16 ## --------------------------------------------------------------## ## Table: Fitting linear model: volume ~ log

15、(Height) + log(Girth),,pander(anova(l.t))## ## -------------------------------------------------------------## &nbsp; Df Sum Sq Mean Sq F value Pr(>F) ## ----------------- ---- -------- --------- --------- --

16、-------## **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.02382 ## -------------------------------------------------------------## ##

17、 Table: Analysis of Variance Table,2.forecast包的BoxCox.lambda和BoxCox,BoxCox.lambda這個函數(shù)用于數(shù)值向量或時間序列，可以得到\lambda的估計精確值。library(forecast),,BoxCox.lambda(trees$Volume, method = "loglik") #算出來的結(jié)果和boxcox有點差異## [1]

18、-0.05volume.f <- BoxCox(trees$Volume, lambda = -0.05)trees.f <- cbind(trees, volume.f) #重新擬合模型l.f <- lm(volume.f ~ log(Height) + log(Girth), data = trees.f) #建立線性模型,,pander(l.f)## ## --------------------

19、------------------------------------------## &nbsp; Estimate Std. Error t value Pr(>|t|) ## ----------------- ---------- ------------ --------- ----------## **(Intercept)** -5.454 0.6731 -8.103 8.013e-09 ## #

20、# **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(Hei

21、ght) + log(Girth),,pander(anova(l.f))## ## -------------------------------------------------------------## &nbsp; Df Sum Sq Mean Sq F value Pr(>F) ## ----------------- ---- -------- --------- --------- ------

22、---## **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 ## -------------------------------------------------------------## ## Ta

23、ble: Analysis of Variance Table,qqPlot(l.f),,3.1 方法性工具,差分運算延遲算子線性差分方程,差分運算,一階差分 階差分 步差分,,,,延遲算子,延遲算子類似于一個時間指針，當前序列值乘以一個延遲算子，就相當于把當前序列值的時間向過去撥了一個時刻 記B為延遲算子，有,,,,延遲算子的性質(zhì),，其中,,,,,,,用延遲算子表示差分運算,階差分 步差分,,,,,,,線

24、性趨勢采用1階差分 y[t]=(1-B) X[t]=X[t]-X[t-1]拋物線趨勢采用2階差分 y[t]=(1-B) ^2X[t]=X[t]+X[t-2]-2X[t-1],季節(jié)差分,y[t]=(1-B^s) X[t]=X[t]-X[t-s] s 為周期 月度數(shù)據(jù) s=12 季度數(shù)據(jù)s=4,對數(shù)變換與差分運算結(jié)合,金融經(jīng)濟數(shù)據(jù) 工業(yè)總產(chǎn)值 y[t]=(1-B)^2 log(X[t]) 股票收盤

25、價 y[t]=(1-B) log(X[t])=log(X[t]/X[t-1]) 稱為股票收益率，采用泰勒展開 約為 y[t]=(X[t]-X[t-1])/X[t-1],線性差分方程,線性差分方程齊次線性差分方程,,,齊次線性差分方程的解,特征方程特征方程的根稱為特征根，記作齊次線性差分方程的通解不相等實數(shù)根場合有相等實根場合復根場合,,,,非齊次線性差分方程的解,非齊次線性差分方程的特解使得非齊

26、次線性差分方程成立的任意一個解非齊次線性差分方程的通解齊次線性差分方程的通解和非齊次線性差分方程的特解之和,,,,異常值處理,數(shù)值檢驗 缺損值的補足模型分析,數(shù)值檢驗,if mean(X[1:t])-k sd(X[1:t])<X[t+1]< mean(X[1:t])+k sd(X[1:t]) k=6, #6 sigma, X[t+1] 正常 否則 為離群點 如果為離群點 X*[t+1]=2X[

27、t]-X[t-1] #線性外推,缺損值的補足,平滑法 插值估算法,模型分析,估計模型 殘差分析 如果殘差分析不能通過，則替換或加入啞變量,例子： 人口自然增長率收集,時間跨度1950-2005時間間隔相同， 年度數(shù)據(jù) ，不允許出現(xiàn)半年，2年等間隔數(shù)據(jù) 統(tǒng)一計算方法：（年初+年末）/2，序列范圍一致性：雖然1997香港回歸，并不能計算進入。 統(tǒng)一指標口徑：出生后但隨后死亡，列入出生也例如死亡。,趨勢因子和季節(jié)因子的估計與

28、去除,時間序列經(jīng)典分解 y[t]=m[t]+S[t]+a[t] m[t]為趨勢因子， S[t] 為季節(jié)因子 a[t]為隨機因子,無季節(jié)因子，僅趨勢因子,y[t]=m[t]+a[t]，去除m[t]方法：最小二乘法滑動平均法差分方法,最小二乘法,用函數(shù)參數(shù)族擬合m[t] 例如 m[t]=a+bt+ct^2 極小化 sum(y[t]- a-bt-ct^2)^2得到 a,b,c的估計,讀取人口,y=scan

29、()3929214530848372398819638453128607021706335323191876314433213855837150189209629797667621216892228496103021537123202624132164569151325798179323175203302031226545805,讀取時間,x=scan()1790180018101820

30、1830184018501860187018801890190019101920193019401950196019701980,,fit=lm(y~x+I(x^2))summary(fit)plot(x,y)lines(fit$fitted.values~x)plot(x,fit$residuals,type="l"),滑動平均（低通濾波，線性濾波器）,雙邊平滑 m[t] 在[

31、t-q,t+q]近似線性 m1[t] =sum(y[t-q:t+q])/(2q+1) =mean(y[t-q:t+q]) =mean(m[t-q:t+q])+mean(a[t-q:t+q]) ≈ mean(m[t-q:t+q]) ≈m[t]mean(a[t-q:t+q]) 快速震蕩 m[t]緩慢變化趨勢,,z=scan()47375117509134684320382536733694

32、3708333333673614336236553963440545955045570057165138501053536074503156485506423048273885,,x1=seq(1951,1980,by=1)n=length(z);m=rep(0,n); q=2for(t in (q+1):(n-q)){m[t]=mean(z[(t-q):(t+q)])},SMOOTH

33、程序， 左邊,#0.1<a<0.3a=0.2 for(t in 1:q){m1=0 for(j in 0: (n-t))m1=m1+a*(1-a)^j*z[t+j] m[t]=m1},右邊,for(t in (n-q+1):n){m1=0 for(j in 0: (t-1))m1=m1+a*(1-a)^j*z[t-j] m[t]=m1},,plot(x1,z,type="l")

34、lines(x1,m,col="red")w=z-mplot(x1,w,type="l"),遞歸關(guān)系,tn-q m[t]=sum_{j=0}^{t-1} a(1-a)^j y[t-j],差分,k 階差分作用k次多項式趨勢項，其結(jié)果為常數(shù)(歸納法可證)人口數(shù)據(jù)2階差分z=diff(y,differences=2)plot(x[-c(1:2)],z,type="l&qu

35、ot;),趨勢項和季節(jié)項同時存在,y[t]=m[t]+S[t]+a[t] 周期d; S[t+d]=S[t] sum(S[1:d])=0例月度數(shù)據(jù) d=12 ，n年 第j 年第k 個月數(shù)據(jù)記為y[k,j],月事故死亡率,年,月,,X=scan()900781068928913710017108261137110744971399389161892777506981803884228714

36、951210120982387439129871086808162730681247870938795561009396208285843381608034771774617776792586348945100789179803784887874864777926957772681068890929910625930283148850826587967

37、836689277918129911594341048498279110907086339240,方法一、小趨勢法,設第j年趨勢相同為M[j],由于 mean(S[1:d])=0 則M[j]=mean(y[1:d,j]) S[k]=mean(y[k,1:n]-M[1:n]),,XM=matrix(X,12,6)ts.plot(X)M=apply(XM,2,mean)d=12S=rep(0,12

38、)for(i in 1:d)S[i]=mean(XM[i,]-M),,n=6;Y1=Xfor(i in 1:n){for(j in 1:d)Y1[j+(i-1)*d]=X[j+(i-1)*d]-M[i]-S[j]}ts.plot(Y1),方法二 、滑動平均,第一步 用滑動平均估計去趨勢項 1.周期d=2q偶數(shù)時 m[t]=(sum(y[t-q+1:t+q-1]) +1/2(y[t-q]+y[

39、t+q]))/d 2.周期為d=2q+1奇數(shù) m[t]=mean(y[t-q:t+q]),,q=d/2m=rep(0,12*6)for(t in (q+1):(12*6-q))m[t]=(sum(X[(t-q+1):(t+q-1)])+X[t-q]/2+X[t+q]/2)/d,第二步去季節(jié)項,w[k]為y[k+j*d]-m[k+j*d], q<=k+jd<=n-q的平均值,這里1<=k<=d 因為

40、sum(w[k])不等于0 所以S[k]=w[k]-mean(w[1: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+X[k+j*d]-m[k+j*d]j=j+1;j1=j1+1}w[k]=w1/j1},,a=0.2;N=length(X) for(t in 1:q){m1=0 f

41、or(j in 0: (N-t))m1=m1+a*(1-a)^j*X[t+j] m[t]=m1},,for(t in (N-q+1):N){m1=0 for(j in 0: (t-1))m1=m1+a*(1-a)^j*X[t-j] m[t]=m1},,Y2=Xfor(i in 1:6){for(j in 1:d)Y2[j+(i-1)*d]=X[j+(i-1)*d]-m[j+(i-1)*d]-S[j]}x1=c

42、(1:n)ts.plot(Y2)lines(x1,Y1,col="red"),方法三、 差分,d步差分去季節(jié)y[t]-y[t-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(