R語言:時間序列經典分析法(二)

      前面一章，介紹了時間序列中涉及到的基本概念，本章將在此基礎上介紹如何對時間序列的資料進行分析，怎麼選擇合適的模型。時間序列的分析方法，已經給大家整理出一個流程圖，如下圖。本章介紹時間序列經典分析法。

      根據前文中，時間序列是四大影響因素的綜合影響結果(長期趨勢，循環周期，季節變動和不規則變動，R語言：時間序列(一))，在各個分量按照加性組合的假定下，對加性模型各分量進行分解，R中使用decompose()函數將各分量分解。下圖為使用加性效應關係繪製的分解圖。

     下圖，從上往下分別是原圖，長期趨勢，季節變動和不規則隨機變動圖。長期趨勢有下降的趨勢；季節分量反應有年度變化周期；隨機分量曲線繞著0水平波動。

        時間序列中的確定性趨勢和季節變動可以用線性回歸模型識別，但時間序列的回歸模型不同於標準回歸模型，因為殘差之間存在自相關性，本質上也是一個時間序列，因此這類模型要仔細檢查殘差。

       1.線性回歸擬合

###書籍集為每年國王的壽命kingsts <- ts(kings)  ##轉成時間序列格式len <- length(kingsts)t <- seq(1,len)/10    ##產生有序變量kingsts.fit <- lm(kingsts~t)  ##配合最小二乘法回歸模型summary(kingsts.fit)  ##輸出主要計算結果par(mfrow=c(3,1),mar=c(3,3,1.6,1),mgp=c(1.8,0.5,0))plot.ts(kingsts)  ##繪製時間序列圖plot(resid(kingsts.fit),cex.axis=0.8,cex.lab=0.8)  ##畫出殘差圖acf(resid(kingsts.fit),cex.axis=0.8,cex.lab=0.8)  ##做殘差的acf圖Box.test(kingsts.fit.fit$residuals) ##對殘差做檢驗## p-value = 0.08575
       2.曲性回歸擬合
kingsts.fit1 <- lm(kingsts~t+I(t^2))  ##二次擬合summary(kingsts.fit1)  ##調整R2=0.2542summary(kingsts.fit)   ##調整R2=0.1419 Box.test(kingsts.fit1$residuals)  ##p-value = 0.4331acf(kingsts.fit1$residuals)  ##無相關性
       3.結果解釋
        1.國王壽命的數據作圖(上上圖1)，可看到有一個先下降後上升的長期趨勢，因此可以做趨勢擬合。
        2.先配合一階回歸模型，並對擬合的模型做殘差檢驗，分別用acf()和Box.test()函數檢驗(上上圖2和3)。可看到殘差圖沒有任何趨勢，殘差相關圖無相關性，殘差檢驗的P值大於0.05，也證實無相關性。
        3.根據其數據先下降後上升的趨勢，可以嘗試擬合下二次曲線，並對擬合完的模型殘差檢驗，其殘差沒有相關性，且比對第一個模型和第二個模型的調整R2發現，第二個模型的R2更大，說明二次曲線更適合這個數據，上圖的藍色線為擬合的直線模型，紅線為擬合的曲線模型，看上去曲線模型擬合更好。
      根據不同的歷史間隔遠近在預報中所起的作用不同而給予不同的權重，歷史較近的觀察值對當前值影響較大，應給予較大的權重；反之，歷史較遠的觀察值對當前值的影響較小，應賦予較小的權重。該方法特別適合變化比較緩慢的時間序列資料的分析。
      簡單指數勻滑法適合不含趨勢和季節變動資料，Holt-Winters法則是將簡單指數勻滑推廣到含有趨勢和季節變動的時間序列資料的分析。在R中都採用HoltWinter()函數分析，其具有下面的基本結構：
HoltWinters(x, alpha = NULL, beta = NULL, gamma = NULL,            seasonal = c("additive", "multiplicative"),            start.periods = 2, l.start = NULL, b.start = NULL,            s.start = NULL,            optim.start = c(alpha = 0.3, beta = 0.1, gamma = 0.1),            optim.control = list())
       函數中的x為時間序列資料的名稱，alpha為水平參數α，beta為趨勢參數β，gamma為季節參數γ。當beta和gamma都賦值為FALSE時，模型就沒有趨勢分量和季節分量，也就是簡單指數勻滑法。當gamma賦值為FALSE時，模型就沒有季節分量。
        簡單指數勻滑法
  空氣中顆粒物PM10數據集為例 pm10.ts <- ts(pm_10$PM10,start = c(2007,1),frequency = 365) pm10.ts.hw <- HoltWinters(pm10.ts,gamma = F,beta = F) coef <- pm10.ts.hw$coefficients with(pm10.ts.hw,c(alpha,beta,gamma,coef,SSE))
       不含季節趨勢指數勻滑
pm10.ts.hw.NS <- HoltWinters(pm10.ts,gamma = F)coef <- pm10.ts.hw.NS$coefficientswith(pm10.ts.hw.NS,c(alpha,beta,gamma,coef,SSE))
       包含長期趨勢和季節趨勢的Holt-Winters勻滑
pm10.ts.hw.S <- HoltWinters(pm10.ts)coef <- pm10.ts.hw.S$coefficients[1:2]with(pm10.ts.hw.S,c(alpha,beta,gamma,coef,SSE)) ##提取勻滑參數和當前基礎水平a和趨勢分量bpm10.ts.hw.S$coefficients[c(3:6,365+2)]  ##提取前4個和最後一個季節的勻滑參數s
     上述三個不同模型的結果解釋：
     1.該序列沒有明顯的長期趨勢，但可見明顯的季節周期，這裡為了和後面對比，因此姑且以此數據作為例子進行簡單指數平滑。其水平參數α=0.73，說明某一天的pm10濃度受近期的影響比較大。其趨勢和季節參數為0，SSE為0.57。
      2. 第二部分，仍以該例子，假使其沒有季節分量gamma=F。看到β=0.02，說明其長期趨勢不明顯。SSE=0.59
      3.第三部分，考慮長期趨勢和季節變動，看到α=0.61，表明近期汙染水平對當前值影響大於遠期汙染水平。β=0表明趨勢分量不需要隨時間更新，γ=0.7，該序列存在明顯季節周期。a=0.082是當前的基礎水平，b=-1.1699e-06是當前的趨勢，s1，s2等則是每天的季節指數。輸出的SSE=0.52，誤差低於前面兩個模型，該模型構建比較合適。
    採用Holt-winters法的加性模型配合時間序列資料的圖形見下圖:
plot(pm10.ts.hw.S$fitted)  ##加性模型擬合圖
     圖中xhat 為用模型參數估計值回代所得的期望值；level為水平值圖，即穩定部分的估計值；trend為時間趨勢圖，結果未顯示趨勢，故為一條平行線；season為季節分量圖，顯示兩個周期。2008-2009和2009到2010，預報未來一年的pm10含量圖形見下圖。
pm10.ts.hw.S.pre <- predict(pm10.ts.hw.S,365)ts.plot(pm10.ts,pm10.ts.hw.S.pre,lty=1:2)
上圖中2007-2010年的實線部分為觀察值，2010-2011年的虛線部分為預測值。
      本文作為時間序列的第二章，主要介紹經典時間序列的分析方法，包括組分分解，趨勢擬合和平滑法，詳細介紹HoltWinter()函數的使用和結果解讀，趨勢擬合的缺點要有長期穩定趨勢，而在商業和經濟中往往只有一段線性特徵，這種情況要用隨機時間序列分析，後面的ARMA模型。
彭曉武等.R軟體及其環境流行病學應用.中國環境出版社
Jonathan D.Cryer .時間序列分析及應用.機械工業出版社
Cory Leismester.精通機器學習：基於R.人民郵電出版社
作者介紹：醫療大數據統計分析師,擅長R語言。
歡迎各位在後臺留言，懇請斧正！