時間序列--ARIMA（原理簡單應用
Autoregressive Integrated Moving Average Model，即自迴歸移動平均模型。它屬於統計模型中最常見的一種，用於進行時間序列的預測。其原理在於：在將非平穩時間序列轉化為平穩時間序列的過程中，將因變數僅對它的滯後值以及隨機誤差項的現值和滯後值進行迴歸所建立的模型。
其實就是三大塊的整合
1.自迴歸model
自迴歸模型是描述當前值與歷史值之間的關係的模型，是一種用變數自身的歷史事件資料對自身進行預測的方法。其公式如下：
其中，yt是當前值；μ是常數項；p是階數；γi是自相關係數，ϵt是誤差值。
自迴歸模型的使用有以下四項限制：
該模型用自身的資料進行預測，即建模使用的資料與預測使用的資料是同一組資料；使用的資料必須具有平穩性；使用的資料必須有自相關性，如果自相關係數小於0.5，則不宜採用自迴歸模型；自回國模型只適用於預測與自身前期相關的現象。
2.integrated
ARIMA模型最重要的地方在於時序資料的平穩性。平穩性是要求經由樣本時間序列得到的擬合曲線在未來的短時間內能夠順著現有的形態慣性地延續下去，即資料的均值、方差理論上不應有過大的變化。平穩性可以分為嚴平穩與弱平穩兩類。嚴平穩指的是資料的分佈不隨著時間的改變而改變；而弱平穩指的是資料的期望與向關係數（即依賴性）不發生改變。在實際應用的過程中，嚴平穩過於理想化與理論化，絕大多數的情況應該屬於弱平穩。對於不平穩的資料，我們應當對資料進行平文化處理。最常用的手段便是差分法，計算時間序列中t時刻與t-1時刻的差值，從而得到一個新的、更平穩的時間序列。
3.moving average
移動平均模型關注的是自迴歸模型中的誤差項的累加。它能夠有效地消除預測中的隨機波動。
1+3=ARMA
在這個公式中，p與q分別為自迴歸模型與移動平均模型的階數，是需要人為定義的。γi與θi分別是兩個模型的相關係數，是需要求解的。如果原始資料不滿足平穩性要求而進行了差分，則為差分自相關移動平均模型（ARIMA），將差分後所得的新資料帶入ARMA公式中即可
流程：
自相關函式（ACF）是將有序的隨機變數序列與其自身相比較，它反映了同一序列在不同時序的取值之間的相關性。
偏自相關函式（PACF）計算的是嚴格的兩個變數之間的相關性，是剔除了中間變數的干擾之後所得到的兩個變數之間的相關程度。對於一個平穩的AR(p)模型，求出滯後為k的自相關係數p(k)時，實際所得並不是x(t)與x(t-k)之間的相關關係。這是因為在這兩個變數之間還存在k-1個變數，它們會對這個自相關係數產生一系列的影響，而這個k-1個變數本身又是與x(t-k)相關的。這對自相關係數p(k)的計算是一個不小的干擾。而偏自相關函式可以剔除這些干擾
- p: The number of lag observations included in the model, also called the lag order.
- d: The number of times that the raw observations are differenced, also called the degree of differencing.
- q: The size of the moving average window, also called the order of moving average.
我們的序列張這樣，可以看出來存在趨勢，不平穩所以需要差分（至少一次差分，你可以用adf檢驗獲取更具有統計意義的結果，參考：https://www.jianshu.com/p/4130bac8ebec）
自迴歸圖長這樣，所以我們大致選擇滯後項為5（更精確的選擇可以參考上面連結）
from pandas import read_csv
from pandas import datetime
from pandas import DataFrame
from statsmodels.tsa.arima_model import ARIMA
from matplotlib import pyplot
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# fit model
model = ARIMA(series, order=(5,1,0))
model_fit = model.fit(disp=0)
print(model_fit.summary())
# plot residual errors
residuals = DataFrame(model_fit.resid)
residuals.plot()
pyplot.show()
residuals.plot(kind='kde')
pyplot.show()
print(residuals.describe())這裡我們選擇lag=5，差分=1，MA=0看看模型
ARIMA Model Results
==============================================================================
Dep. Variable: D.Sales No. Observations: 35
Model: ARIMA(5, 1, 0) Log Likelihood -196.170
Method: css-mle S.D. of innovations 64.241
Date: Mon, 12 Dec 2016 AIC 406.340
Time: 11:09:13 BIC 417.227
Sample: 02-01-1901 HQIC 410.098
- 12-01-1903
=================================================================================
coef std err z P>|z| [95.0% Conf. Int.]
---------------------------------------------------------------------------------
const 12.0649 3.652 3.304 0.003 4.908 19.222
ar.L1.D.Sales -1.1082 0.183 -6.063 0.000 -1.466 -0.750
ar.L2.D.Sales -0.6203 0.282 -2.203 0.036 -1.172 -0.068
ar.L3.D.Sales -0.3606 0.295 -1.222 0.231 -0.939 0.218
ar.L4.D.Sales -0.1252 0.280 -0.447 0.658 -0.674 0.424
ar.L5.D.Sales 0.1289 0.191 0.673 0.506 -0.246 0.504
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -1.0617 -0.5064j 1.1763 -0.4292
AR.2 -1.0617 +0.5064j 1.1763 0.4292
AR.3 0.0816 -1.3804j 1.3828 -0.2406
AR.4 0.0816 +1.3804j 1.3828 0.2406
AR.5 2.9315 -0.0000j 2.9315 -0.0000
-----------------------------------------------------------------------------
最後那個roots其實沒太明白是什麼玩意。。。。
上面是殘差圖，可以看出來並不會雜亂無章，所以說明我們的模型不是很好
count 35.000000
mean -5.495213
std 68.132882
min -133.296597
25% -42.477935
50% -7.186584
75% 24.748357
max 133.237980
殘差均值也不=0，所以確實模型不好啊
按理說殘差表現好才能預測，那先假設殘差表現好把
如果我們使用訓練資料集中的100個觀測值來擬合模型，那麼預測的下一個時間步的索引將被指定給預測函式start=101, end=101。這將返回一個包含預測的元素的陣列。
如果在配置模型時執行了任何不同操作(d>0)，我們也希望預測的值保持在原始的比例。這可以通過將typ引數設定為值“levels”來指定:typ=“levels”。
另外，我們可以通過使用forecast()函式來避免所有這些規範，該函式使用模型執行一步預測。我們可以將訓練資料集分為訓練集和測試集，使用訓練集來擬合模型，併為測試集上的每個元素生成一個預測。
考慮到對差分和AR模型的依賴，需要滾動預測。執行滾動預測的一種粗略方法是在每次接收到新的觀測之後重新建立ARIMA模型。
我們手動跟蹤一個名為history的列表中的所有觀察結果，這個列表中包含了培訓資料，並且每個迭代都會新增新的觀察結果。把這些放在一起，下面是一個使用Python中的ARIMA模型滾動預測的示例。
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot
from statsmodels.tsa.arima_model import ARIMA
from sklearn.metrics import mean_squared_error
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
X = series.values
size = int(len(X) * 0.66)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = list()
for t in range(len(test)):
model = ARIMA(history, order=(5,1,0))
model_fit = model.fit(disp=0)
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs) # 這裡實現了滾動預測，應該是ARIMA只能預測下一步而不能多步
print('predicted=%f, expected=%f' % (yhat, obs))
error = mean_squared_error(test, predictions)
print('Test MSE: %.3f' % error)
# plot
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
效果圖如上
https://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/
使用建議：
https://machinelearningmastery.com/gentle-introduction-box-jenkins-method-time-series-forecasting/
https://blog.csdn.net/qq_27123591/article/details/80272669
這裡有兩個具體的例子實現了ARIMA
1.https://machinelearningmastery.com/time-series-forecast-study-python-annual-water-usage-baltimore/
這兩個例子介紹瞭如何根據acf、pacf選擇ARIMA的pdq
這裡我們用predict實現了ARIMA的滾動預測，也可以手動實現（雖然不必要）可以參考：
https://machinelearningmastery.com/make-manual-predictions-arima-models-python/