時間序列分析之ARIMA模型預測_資料採礦

轉載自 http://blog.sina.com.cn/s/blog_70f632090101bnd8.html#cmt_3111974

今天學習ARIMA預測時間序列。

 指數平滑法對於預測來說是非常有協助的，而且它對時間序列上面連續的值之間相關性沒有要求。但是，如果你想使用指數平滑法計算出預測區間， 那麼預測誤差必須是不相關的， 而且必須是服從零均值、 方差不變的常態分佈。即使指數平滑法對時間序列連續數值之間相關性沒有要求，在某種情況下， 我們可以通過考慮資料之間的相關性來建立更好的預測模型。自迴歸移動平均模型（ ARIMA） 包含一個確定（explicit） 的統計模型用於處理時間序列的不規則部分，它也允許不規則部分可以自相關。

首先，先確定資料的差分。

ARIMA 模型為平穩時間序列定義的。 因此， 如果你從一個非平穩的時間序列開始， 首先你就需要做時間序列差分直到你得到一個平穩時間序列。如果你必須對時間序列做 d 階差分才能得到一個平穩序列，那麼你就使用ARIMA(p,d,q)模型，其中 d 是差分的階數。 

我們以每年女人裙子邊緣的直徑做成的時間序列資料為例。從 1866 年到 1911 年在平均值上是不平穩的。 隨著時間增加， 數值變化很大。  

> skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)

Read 46 items

> skirtsts<- ts(skirts,start = c(1866))

 

> plot.ts(skirtsts)

我們可以通過鍵入下面的代碼來得到時間序列（資料存於“skirtsts”） 的一階差分， 並畫出差分序列的圖:

> skirtstsdiff<-diff(skirtsts,differences=1)

 

> plot.ts(skirtstsdiff)

從一階差分的圖中可以看出，資料仍是不平穩的。我們繼續差分。

> skirtstsdiff2<-diff(skirtsts,differences=2)

> plot.ts(skirtstsdiff2)

二次差分（上面）後的時間序列在均值和方差上確實看起來像是平穩的， 隨著時間推移， 時間序列的水平和方差大致保持不變。因此， 看起來我們需要對裙子直徑進行兩次差分以得到平穩序列。

第二步，找到合適的ARIMA模型

 如果你的時間序列是平穩的，或者你通過做 n 次差分轉化為一個平穩時間序列， 接下來就是要選擇合適的 ARIMA模型，這意味著需要尋找 ARIMA(p,d,q)中合適的 p 值和 q 值。為了得到這些，通常需要檢查[平穩時間序列的（自）相關圖和偏相關圖。 

 我們使用 R 中的“acf()”和“pacf” 函數來分別（ 自） 相關圖和偏相關圖。“acf()”和“pacf 設定“plot=FALSE” 來得到自相關和偏相關的真實值。 

> acf(skirtstsdiff2,lag.max=20)

 

> acf(skirtstsdiff2,lag.max=20,plot=FALSE)

Autocorrelations of series ‘skirtstsdiff2’, by lag

     0      1      2      3      4      5      6      7      8     9     10 

 1.000 -0.303  0.096  0.009  0.102 -0.453  0.173 -0.025 -0.039  0.073 -0.094 

    11     12     13     14     15     16     17     18     19     20 

 0.133 -0.089 -0.027 -0.102  0.207 -0.260  0.114  0.101  0.011 -0.090 

自相關圖顯示滯後1階自相關值基本沒有超過邊界值，雖然5階自相關值超出邊界，那麼很可能屬於偶然出現的，而自相關值在其他上都沒有超出顯著邊界， 而且我們可以期望 1 到 20 之間的會偶爾超出 95%的置信邊界。  

> pacf(skirtstsdiff2,lag.max=20)

> pacf(skirtstsdiff2,lag.max=20,plot=FALSE)

Partial autocorrelations of series ‘skirtstsdiff2’, by lag

     1      2      3      4      5      6      7      8      9    10     11 

-0.303  0.005  0.043  0.128 -0.439 -0.110  0.073  0.028  0.128 -0.355  0.095 

    12     13     14     15     16     17     18     19     20 

 0.052 -0.094 -0.103 -0.034 -0.021 -0.002  0.074  0.020 -0.034 

偏自相關值選5階。

故我們的ARMIA模型為armia（1,2,5）

> skirtsarima<-arima(skirtsts,order=c(1,2,5))

> skirtsarima

SSeries: skirtsts 

ARIMA(1,2,5)                    

Coefficients:

          ar1     ma1     ma2     ma3     ma4      ma5

      -0.4345  0.2762  0.1033  0.1472  0.0267  -0.8384

s.e.   0.1837  0.2171  0.2198  0.2716  0.1904   0.2888

sigma^2 estimated as 206.1:  log likelihood=-183.8

AIC=381.6   AICc=384.71   BIC=394.09

預測後5年裙子的邊緣直徑

>  skirtsarimaforecast<-forecast.Arima(skirtsarima,h=5,level=c(99.5))

> skirtsarimaforecast

     Point Forecast  Lo 99.5  Hi 99.5

1912       548.5762 507.1167 590.0357

1913       545.1793 459.3292 631.0295

1914       540.9354 396.3768 685.4940

1915       531.8838 316.2785 747.4892

1916       529.1296 233.2625 824.9968

> plot.forecast(skirtsarimaforecast$residuals)   #謝謝@憶水如煙的指正

第三步，檢驗

在指數平滑模型下， 觀察 ARIMA 模型的預測誤差是否是平均值為 0 且方差為常數的常態分佈（服從零均值、方差不變的常態分佈） 是個好主意，同時也要觀察連續預測誤差是否（自）相關。  

> acf(skirtsarimaforecast$residuals,lag.max=20)

> Box.test(skirtsarimaforecast$residuals, lag=20, type="Ljung-Box")

        Box-Ljung test

data:  skirtsarimaforecast$residuals 

X-squared = 8.5974, df = 20, p-value = 0.9871

既然相 關圖顯示出在滯後1 - 20階（ l a g s 1 - 20 ）中樣本自相關值都沒有超出顯著（置信）邊界，而且Ljung-Box檢驗的p值為0.99，所以我們推斷在滯後1-20階（lags1-20）中沒明顯證據說明預測誤差是非零自相關的。 

為了調查預測誤差是否是平均值為零且方差為常數的常態分佈（服從零均值、方差不變的常態分佈），我們可以做預測誤差的時間曲線圖和長條圖（具有常態分佈曲線）：

> plot.ts(skirtsarimaforecast$residuals) > plotForecastErrors(skirtsarimaforecast$residuals)

上圖預測中的時間曲線圖顯示出對著時間增加，方差大致為常數（大致不變）（儘管上半部分的時間序

列方差看起來稍微高一些）。時間序列的長條圖顯示預測誤大致是常態分佈的且平均值接近於 0（服從零均值的常態分佈的）。因此，把預測誤差看作平均值為0方差為常數常態分佈（服從零均值、方差不變的常態分佈）是合理的。 

 

既然依次連續的預測誤差看起來不是相關，而且看起來是平均值為 0 方差為常數的常態分佈（服從零均值、方差不變的常態分佈），那麼對於裙子直徑的資料， ARIMA(1,2,5)看起來是可以提供非常合適預測的模型。  至此，時間序列的學習結束