本文結構：

• 時間序列分析？
• 什麼是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模型運用的基本流程有幾下幾步：

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

使用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)

下圖顯示時間序列的直方圖顯示預測誤大致是正態分佈的且平均值接近於0。因此，把預測誤差看作平均值為0方差為服從零均值、方差不變的正態分佈是合理的。

既然依次連續的預測誤差看起來不是相關，而且服從零均值、方差不變的正態分佈，那麼對於裙子直徑的資料，ARIMA(1,2,5)看起來是可以提供非常合適預測的模型。

常見問題？

1.ARIMA建模的步驟

1. 觀察資料是否是時間序列資料，是否有seasonal等因素。
2. transform：Box－correlation，保證variance是uniform的。如果用box－cor還不能穩定，還要繼續深入挖掘。
3. ACF／PACF 是為了找到 MA 和 AR 的order。
4. d＝0－stationarity，1，2－non stationarity
5. 白噪音check：確定這個模型是optimize的，mean＝0，平方差＝1.
6. 誤差是白噪音的時候，model就ok了，就可以預測了
   
   

2.決定ARIMA引數的方法

d是差分的階數，幾階後就可以保證穩定：

modelfit，計算出來的引數是 1，1，1 ，但可能 2，1，1 預測效果更好，那就用後者。

或者用AIC比較倆模型。

推薦閱讀：

這一篇例項也不錯：python時間序列分析

AR和MA的定義，有圖比較：

關於ACF，PACF可以看Duke的材料：

關於AIC，BIC：aic-bic-vs-crossvalidation

理論簡介的不錯：

需求預測與分倉規劃:

一個課程：Time Series Forecasting Theory | AR, MA, ARMA, ARIMA