【機器學習筆記之五】用ARIMA模型做需求預測用ARIMA模型做需求預測

本文轉載自查看原文 2017-08-17 08:36 14775 A---機器學習

本文結構：

時間序列分析？
什么是ARIMA？
ARIMA數學模型？
input，output 是什么？
怎么用？－代碼實例
常見問題？

時間序列分析？

時間序列，就是按時間順序排列的，隨時間變化的數據序列。
生活中各領域各行業太多時間序列的數據了，銷售額，顧客數，訪問量，股價，油價，GDP，氣溫。。。

隨機過程的特征有均值、方差、協方差等。
如果隨機過程的特征隨着時間變化，則此過程是非平穩的；相反，如果隨機過程的特征不隨時間而變化，就稱此過程是平穩的。
下圖所示，左邊非穩定，右邊穩定。

非平穩時間序列分析時，若導致非平穩的原因是確定的，可以用的方法主要有趨勢擬合模型、季節調整模型、移動平均、指數平滑等方法。
若導致非平穩的原因是隨機的，方法主要有ARIMA（autoregressive integrated moving average）及自回歸條件異方差模型等。

什么是ARIMA？

ARIMA (Auto Regressive Integrated Moving Average) 可以用來對時間序列進行預測，常被用於需求預測和規划中。

可以用來對付 ‘隨機過程的特征隨着時間變化而非固定’ 且 ‘導致時間序列非平穩的原因是隨機而非確定’ 的問題。不過，如果是從一個非平穩的時間序列開始，首先需要做差分，直到得到一個平穩的序列。

模型的思想就是從歷史的數據中學習到隨時間變化的模式，學到了就用這個規律去預測未來。

ARIMA(p,d,q)模型，其中 d 是差分的階數，用來得到平穩序列。

AR是自回歸, p為相應的自回歸項。

MA為移動平均，q為相應的移動平均項數。

ARIMA數學模型？

ARIMA（p，d，q）模型是ARMA（p，q）模型的擴展。

ARIMA（p，d，q）模型可以表示為：

其中L 是滯后算子（Lag operator），d in Z, d>0。

AR：
當前值只是過去值的加權求和。

MA：
過去的白噪音的移動平均。

ARMA：
AR和MA的綜合。

ARIMA：
和ARMA的區別，就是公式左邊的x變成差分算子，保證數據的穩定性。

差分算子就是：

令 wt 為：

則 ARIMA 就可以寫成：

input，output 是什么？

輸入歷史數據，預測未來時間點的數據。

怎么用？－代碼實例

本文參考了：時間序列實例
另外推薦大家看這篇，36大數據上有一個python版講的不錯，里面對穩定性的定量檢驗的講解比較詳細：時間序列預測全攻略－附帶Python代碼

ARIMA模型運用的基本流程有幾下幾步：

數據可視化，識別平穩性。
對非平穩的時間序列數據，做差分，得到平穩序列。
建立合適的模型。
平穩化處理后，若偏自相關函數是截尾的，而自相關函數是拖尾的，則建立AR模型；
若偏自相關函數是拖尾的，而自相關函數是截尾的，則建立MA模型；
若偏自相關函數和自相關函數均是拖尾的，則序列適合ARMA模型。
模型的階數在確定之后，對ARMA模型進行參數估計，比較常用是最小二乘法進行參數估計。
假設檢驗，判斷（診斷）殘差序列是否為白噪聲序列。
利用已通過檢驗的模型進行預測。

使用ARIMA模型對裙子長度預測

1、加載數據

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

str(skirts)
head(skirts)
boxplot(skirts)
length(skirts)

2、把數據轉化為是時間序列

skirts_ts <- ts(skirts, start=c(1886), frequency=1)

1)查看時間序列對應的時間

skirts_ts

2)畫出時間序列圖

plot.ts(skirts_ts)

從圖可知：女人裙子邊緣的直徑做成的時間序列數據，從 1866 年到 1911 年在平均值上是不平穩的

3、做差分得到平穩序列

1)做時間序列的一階差分

skirts_diff <- diff(skirts_ts, differences = 1)
plot.ts(skirts_diff)

從一階差分的圖中可以看出，數據仍是不平穩的，繼續差分

2)做時間序列的二階差分

skirts_diff2 <- diff(skirts_ts, differences = 2)
plot.ts(skirts_diff2)

二次差分后的時間序列在均值和方差上看起來是平穩了

4、找到合適的ARIMA模型

尋找 ARIMA(p,d,q)中合適的 p 值和 q

1)自相關圖ACF

acf(skirts_diff2, lag.max = 20)

acf(skirts_diff2, lag.max = 20, plot = F)

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

2)偏相關圖PACF

pacf(skirts_diff2, lag.max = 20)

pacf(skirts_diff2, lag.max = 20, plot = F)

偏自相關值選1階后結尾
故我們的ARMIA模型為armia（1,2,5

3)使用auto.arima()函數，自動獲取最佳的ARIMA模型

library(forecast)

auto.arima(skirts_ts, ic=c("aicc", "aic", "bic"), trace = T)

Best model: ARIMA(1,2,0)

5、建立ARIMA模型：並對比arima(1, 2, 0)與arima(1, 2, 5)模型

1）arima(1, 2, 0)模型

(skirts_arima <- arima(skirts_ts, order = c(1, 2, 0)))

aic = 391.33

2）arima(1, 2, 5)模型

(skirts_arima <- arima(skirts_ts, order = c(1, 2, 5)))

aic = 381.6

AIC是赤池消息准則SC是施瓦茨准則，當兩個數值最小時，則是最優滯后分布的長度。我們進行模型選擇時，AIC值越小越好。所以arima(1, 2, 5)模型較好

6、預測：預測5年后裙子的邊緣直徑

(skirts_forecast <- forecast.Arima(skirts_arima, h=5, level = c(99.5)))

plot.forecast(skirts_forecast)

7、檢驗

觀察 ARIMA 模型的預測誤差是否是平均值為 0 且方差為常數的正態分布，同時也要觀察連續預測誤差是否自相關

1）檢驗預測誤差的自相關性

tsdiag(skirts_arima)

下面第一個圖表代表估計模型誤差的繪圖。圖中豎線的長度比較相似，都處在穩定范圍之內，即估計的模型沒產生不符合要求的誤差分布。

第二張繪圖，顯示估計的模型沒造成誤差之間的任何關系。這是符合數據生成時每個數據都是獨立的這個前提的。由此可見，這ACF圖符合檢測要求。

第三張圖，也就是Ljung-Box 指標。這個指標可對每一個時間序列的延遲進行顯著性的評估。判定技巧是，P-value點的高度越高，我們的模型越可信。

acf(skirts_forecast$residuals, lag.max = 20)

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

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

Acf檢驗說明：殘差沒有明顯的自相關性，Ljung-Box測試顯示：所有的P-value>0.05，說明殘差為白噪聲。

2）判斷預測誤差是否是平均值為零且方差為常數的正態分布
做預測誤差的時間曲線圖和直方圖（具有正態分布曲線）

預測誤差的均值是否為0

plot.ts(skirts_forecast$residuals)

自定義判斷預測誤差的方差是正態分布的函數

plotForecastErrors <- function(forecasterrors){
  #畫預測誤差的直方圖
  hist(forecasterrors, col="red", freq = F)
  #畫方差是預測誤差數據的方差，平均值是0的正態分布數據的線
  mysd <- sd(forecasterrors)
  mynorm <- rnorm(10000, mean = 0, sd = mysd)
  myhist <- hist(mynorm, plot = F)
  points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)
}
plotForecastErrors(skirts_forecast$residuals)