時間序列預測算法-ARIMA算法

https://www.biaodianfu.com/time-series-forecasting-with-arima-in-python.html

https://www.biaodianfu.com/arima-p-d-q.html

https://blog.csdn.net/qifeidemumu/article/details/88761964?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param
基於時間序列分析的趨勢預測算法，主體采用ARIMA差分自回歸移動平均模型，ARIMA算法模型主體包括三大部分：AR，I以及MA模型。其中，每一個模型部分都擁有一個相關的模型參數—ARIMA（p，d，q）。算法的基本原理是將非平穩時間序列轉化為平穩時間序列然后將因變量僅對它的滯后值以及隨機誤差項的現值和滯后值進行回歸所建立的模型。

在ARIMA（p，d，q）整體模型中，AR是自回歸模型，對應的模型參數p為自回歸項數；I為差分模型，對應的模型參數d為使之成為平穩時間序列所做的差分次數（階數）；MA為滑動平均模型，q為滑動平均項數。在實際進行算法模型的構建時，可以根據ACF自相關系數圖決定q的取值，PACF偏自相關系數圖決定p的取值。

1. 自回歸模型（AR）

AR模型是為了建立起當前數據特征值與過去歷史值之間的關系，實現用變量自身的歷史時間數據對自身進行一定的時間周期預測，自回歸模型必須滿足平穩性的要求，並且必須具有自相關性，自相關系數小於0.5則不適用

p階自回歸過程的公式定義：

  

其中，為當前的值，是常數項，p是階數，是自相關系數，是誤差，前幾天值的大小。

2. I差分模型

I差分模型主要是為了對原始數據進行不同書數目階次的差分處理，使得原始數據轉變為時間維度上的平穩序列，為后續建立AR和MA算法模型奠定基礎。一般情況下，采用差分階數最多的是一階或者兩階。

3. 移動平均模型（MA）

移動平均模型關注的是自回歸模型中的誤差項的累加，移動平均法能有效地消除預測中的隨機波動

q階移動平均過程的公式定義：

     

4. ARIMA算法模型的pq參數的確定

PACF，是指偏自相關函數，剔除了中間k-1個隨機變量x(t-1)、x(t-2)、……、x(t-k+1)的干擾之后x(t-k)對x(t)影響的相關程度,可以通過PACF的分布圖來選擇出p值的最優值大小。

ACF，自相關函數反映了同一序列在不同時序的取值之間的相關性。x(t)同時還會受到中間k-1個隨機變量x(t-1)、x(t-2)、……、x(t-k+1)的影響而這k-1個隨機變量又都和x(t-k)具有相關關系，所以自相關系數p(k)里實際摻雜了其他變量對x(t)與x(t-k)的影響，可以通過它的分布圖來選擇出最佳的模型q參數。

 

 

 

 

 

 

 

 

圖6.1  ACF圖與PACF圖

綜上，其具體的確定原則如下表所示：

表6-1  ARIMA模型pq參數的確定原則

5. ARIMA算法的具體步驟

①  時間序列可視化；

②    序列平穩化處理（進行d階差分處理）；

③  繪制ACF與PACF圖，尋找ARIMA模型最優p和q參數；

④  建立ARIMA模型；

⑤  進行指定周期時間內數據的預測。

6.算法構建與性能測試

針對過程測試PDL_AM的port1空間位置（光頻為190.65）數據，選取時間區間為2020年8月1日到8月16日，構建ARIMA時間序列預測算法模型。

1．數據平穩化處理

針對原始數據，首先要進行數據的平穩化處理，經過測試，可以得到原始數據進過一階差分之后即可實現數據的平穩化，即ARIMA模型的I部分的階數d為1，下圖分別為原始數據和原始數據經過1階差分之后的數據趨勢變化曲線：

 

 

 

圖6.2 過程測試原始數據變化趨勢曲線

 

 

 

圖6.3 過程測試原始數據一階差分變化趨勢曲線

2. 繪制一階差分數據的自相關圖和偏相關圖
       繪制出一階差分處理完的PDL數據的自相關圖和偏相關圖，如下圖所示：

 

圖6.4 一階差分數據的自相關圖與偏相關圖

根據自相關圖和偏相關圖可以得到ARMA模型的最佳p和q參數分別為1和1（原因是滯后1階之后其相關系數趨於0）

1. 構建ARIMA算法模型

指定p=1,q=1,d=1，構建ARIMA時間趨勢變化算法模型，並進行模型的擬合和診斷，得到以下的模型診斷輸出圖：

 

圖6.5 模型診斷輸出圖

（1）右上角的紅色KDE線和黃色N(0,1)線非常的接近，它表明殘差是服從正態分布的。

（2）左下角的QQ分位圖表明，從N(0,1)的標准正態分布抽取的樣本，殘差(藍點)的有序分布服從線性趨勢。 同樣，這是一個強烈的跡象表明，殘差是正態分布。

（3）隨着時間的推移（左上圖）殘差不會顯示任何明顯的季節性，似乎是白噪聲。這通過右下角的自相關（即相關圖）證實了這一點，它表明時間序列的殘差與其自身的滯后版本有很低的相關性。

4. 驗證構建的模型

    使用原始數據進行ARIMA算法模型的驗證，將序號200-300之間的100個數據的真實值與算法模型的預測值做出來如下圖所示：

 

圖6.6 過程測試ARIMA算法模型驗證

可有得到算法模型根據前面數據的預測值與真實值的,MSE（均方差）為0.01，說明算法模型具有較好的趨勢預測性能。

5. 趨勢預測

使用建立起的ARIMA模型對后續三天指標值的變化進行預測，得到了未來三天2020年8月17-8月19日的21條預測結果（每天平均測試7個數據），具體變化如下所示：

 

圖6.7 過程測試PDL指標數據未來三天趨勢預測圖

源代碼：

第一種算法實現：

import itertools
import pandas as pd
import matplotlib
matplotlib.use("TkAGG")
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
data = sm.datasets.co2.load_pandas().data
data = data['co2'].resample('MS').mean()
data = data.fillna(data.bfill())
print(data.head())
data.plot(figsize=(12, 6))
plt.show()
print(data.shape)
plt.rcParams["font.sans-serif"]=["SimHei"]
plt.rcParams["axes.unicode_minus"]=False

data1=pd.read_excel("過程測試_033GRR10L4105623_20200601-20200708_PDL_AM_異常檢測分類結果.xlsx")
i=500
newdata=data1["p1"][:300]
# newdata.index=data1["start_time"][:300]
newdata.index=range(len(newdata))
newdata.plot(figsize=(12, 6))
plt.title("過程測試20200601-0616PDL數據趨勢變化曲線")
plt.xlabel("時間")
plt.ylabel("指標大小")
plt.show()

#2.下面我們先對非平穩時間序列進行時間序列的差分，找出適合的差分次數d的值：
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(111)
diff1 = newdata.diff(1)
diff1.plot(ax=ax1)
plt.title("過程測試20200601-0616PDL數據一階差分曲線")
plt.xlabel("時間")
plt.ylabel("指標大小")
plt.show()

#這里是做了1階差分，可以看出時間序列的均值和方差基本平穩，不過還是可以比較一下二階差分的效果：

#這里進行二階差分
fig = plt.figure(figsize=(12, 8))
ax2 = fig.add_subplot(111)
diff2 = newdata.diff(2)
diff2.plot(ax=ax2)
plt.title("過程測試20200601-0616PDL數據二階差分曲線")
plt.xlabel("時間")
plt.ylabel("指標大小")
plt.show()

#由下圖可以看出來一階跟二階的差分差別不是很大，所以可以把差分次數d設置為1，上面的一階和二階程序我們注釋掉
plt.show()
#這里我們使用一階差分的時間序列

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
print(pdq)
# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print(seasonal_pdq)
'''
import warnings
warnings.filterwarnings("ignore")  # specify to ignore warning messages
best_results=1e10
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            model = sm.tsa.statespace.SARIMAX(newdata, order=param, seasonal_order=param_seasonal,enforce_stationarity=False, enforce_invertibility=False)
            results = model.fit()
            print('ARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic))
            if best_results>results.aic:
                best_results=results.aic
                best_pdq=param
                best_pdqs=param_seasonal
        except:
            continue
print("最好的參數組合為：",best_pdq,best_pdqs)
'''
model = sm.tsa.statespace.SARIMAX(newdata, order=(1,2,1), seasonal_order=(1, 1, 1, 12),enforce_stationarity=False, enforce_invertibility=False)
results = model.fit()
print(results.summary().tables[1])
results.plot_diagnostics(figsize=(12, 12))
plt.show()

pred = results.get_prediction(start=newdata.index[200], dynamic=False)
pred_ci = pred.conf_int()

ax =newdata.plot(label='Observed', figsize=(12, 6))
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)

ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)

ax.set_xlabel('時間')
ax.set_ylabel('指標大小')
plt.title("過程測試ARIMA時間序列預測模型驗證")
plt.legend()
plt.show()

#直接輸出相應的MSE,比較對應函數的趨勢分析
data_forecasted = pred.predicted_mean
data_truth = newdata[200:]
print(data_forecasted,len(data_forecasted))
print(data_truth,len(data_truth))
plt.plot(range(len(data_forecasted)-1),data_forecasted[1:],"r")
plt.plot(range(len(data_truth)),data_truth,"b")
plt.show()

# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
'''
#使用動態預測的方法進行數據的預測與變化
pred_dynamic = results.get_prediction(start=pd.to_datetime(newdata.index[450]), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()

ax = newdata.plot(label='Observed', figsize=(12, 6))
pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)

ax.fill_between(pred_dynamic_ci.index,
                pred_dynamic_ci.iloc[:, 0],
                pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)

ax.fill_betweenx(ax.get_ylim(),pd.to_datetime(newdata.index[450]), newdata.index[-1],
                 alpha=.1, zorder=-1)
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()
plt.show()

data_forecasted = pred_dynamic.predicted_mean
data_truth = newdata[450:]
print(data_forecasted,len(data_forecasted))
print(data_truth,len(data_truth))

# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

plt.plot(range(len(data_forecasted)),data_forecasted,"r")
plt.plot(range(len(data_truth)),data_truth,"b")
plt.show()
'''

# Get forecast 500 steps ahead in future
# model = sm.tsa.statespace.SARIMAX(newdata, order=(1,2,1), seasonal_order=(1,1,1,12), enforce_stationarity=False, enforce_invertibility=False)
#results = model.fit()
pred_uc = results.get_forecast(steps=20)
print(pred_uc.predicted_mean)
# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

plt.figure()
plt.plot(newdata)
plt.plot(pred_uc.predicted_mean)
plt.show()

ax = newdata.plot(label='Observed', figsize=(12, 6))

#plt.show()
pred_uc.predicted_mean.plot(ax=ax,label='Forecast')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.25)
ax.set_xlabel('時間')
ax.set_ylabel('指標大小')
plt.title("過程測試ARIMA模型預測未來三天指標變化趨勢")
plt.legend()
plt.show()

第二種算法實現：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.arima_model import ARMA
from statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacf
from statsmodels.graphics.api import qqplot
plt.rcParams["font.sans-serif"]=["SimHei"]
plt.rcParams["axes.unicode_minus"]=False
# 1.創建數據
data = [5922, 5308, 5546, 5975, 2704, 1767, 4111, 5542, 4726, 5866, 6183, 3199, 1471, 1325, 6618, 6644, 5337, 7064, 2912, 1456, 4705, 4579, 4990, 4331, 4481, 1813, 1258, 4383, 5451, 5169, 5362, 6259, 3743, 2268, 5397, 5821, 6115, 6631, 6474, 4134, 2728, 5753, 7130, 7860, 6991, 7499, 5301, 2808, 6755, 6658, 7644, 6472, 8680, 6366, 5252, 8223, 8181, 10548, 11823, 14640, 9873, 6613, 14415, 13204, 14982, 9690, 10693, 8276, 4519, 7865, 8137, 10022, 7646, 8749, 5246, 4736, 9705, 7501, 9587, 10078, 9732, 6986, 4385, 8451, 9815, 10894, 10287, 9666, 6072, 5418]

data = pd.Series(data)
data.index = pd.Index(sm.tsa.datetools.dates_from_range('1901','1990'))
data.plot(figsize=(12,8))
plt.show()
#
# data1=pd.read_excel("過程測試_033GRR10L4105623_20200601-20200708_PDL_AM_異常檢測分類結果.xlsx")
# i=500
# data=data1["p1"][:300]
# data.index=range(len(data))
# #data.index = pd.Index(sm.tsa.datetools.dates_from_range('1691','1990'))
# data.plot(figsize=(12, 6))
# plt.title("過程測試20200601-0616PDL數據趨勢變化曲線")
# plt.xlabel("時間")
# plt.ylabel("指標大小")
# plt.show()

#2.下面我們先對非平穩時間序列進行時間序列的差分，找出適合的差分次數d的值：
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(111)
diff1 = data.diff(1)
diff1.plot(ax=ax1)
plt.title("過程測試20200601-0616PDL數據一階差分曲線")
plt.xlabel("時間")
plt.ylabel("指標大小")
#這里是做了1階差分，可以看出時間序列的均值和方差基本平穩，不過還是可以比較一下二階差分的效果：

#這里進行二階差分
fig = plt.figure(figsize=(12, 8))
ax2 = fig.add_subplot(111)
diff2 = data.diff(2)
diff2.plot(ax=ax2)
#由下圖可以看出來一階跟二階的差分差別不是很大，所以可以把差分次數d設置為1，上面的一階和二階程序我們注釋掉
plt.show()
#這里我們使用一階差分的時間序列

data = data.diff(1)
data.dropna(inplace=True)
#加上這一步，不然后面畫出的acf和pacf圖會是一條直線

#第一步：先檢查平穩序列的自相關圖和偏自相關圖
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(data,lags=40,ax=ax1)
#lags 表示滯后的階數
#第二步：下面分別得到acf 圖和pacf 圖
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(data, lags=40,ax=ax2)
plt.show()


#由上圖可知，我們可以分別用ARMA(0,1)模型、ARMA(7,0)模型、ARMA(7,1)模型等來擬合找出最佳模型：
#第三步:找出最佳模型ARMA
arma_mod1 = sm.tsa.ARMA(data,(7,0)).fit()
print(arma_mod1.aic, arma_mod1.bic, arma_mod1.hqic)
# arma_mod2 = sm.tsa.ARMA(data,(0,1)).fit()
# print(arma_mod2.aic, arma_mod2.bic, arma_mod2.hqic)
# arma_mod3 = sm.tsa.ARMA(data,(7,1)).fit()
# print(arma_mod3.aic, arma_mod3.bic, arma_mod3.hqic)
# arma_mod4 = sm.tsa.ARMA(data,(8,0)).fit()
# print(arma_mod4.aic, arma_mod4.bic, arma_mod4.hqic)

#由上面可以看出ARMA(7,0)模型最佳
#由上面可以看出ARMA(7,0)模型最佳
#第四步：進行模型檢驗
#首先對ARMA(7,0)模型所產生的殘差做自相關圖
resid = arma_mod1.resid
#一定要加上這個變量賦值語句，不然會報錯resid is not defined
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(),lags=40,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40,ax=ax2)
plt.show()

#接着做德賓-沃森（D-W）檢驗
print(sm.stats.durbin_watson(arma_mod1.resid.values))
#得出來結果是不存在自相關性的

#再觀察是否符合正態分布,這里用qq圖
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q',ax=ax, fit=True)
plt.show()
# #最后用Ljung-Box檢驗:檢驗的結果就是看最后一列前十二行的檢驗概率（一般觀察滯后1~12階），
# #如果檢驗概率小於給定的顯著性水平，比如0.05、0.10等就拒絕原假設，其原假設是相關系數為零。
# #就結果來看，如果取顯著性水平為0.05，那么相關系數與零沒有顯著差異，即為白噪聲序列。
r,q,p = sm.tsa.acf(resid.values.squeeze(),qstat=True)
data1 = np.c_[range(1,41), r[1:], q, p]
table= pd.DataFrame(data1, columns=[ 'lag','AC','Q','Prob(>Q)'])
print(table.set_index('lag'))

#第五步：模型預測,對未來十年進行預測
predict_y =arma_mod1.predict('1990', '2000', dynamic=True)
#print(predict_y)
fig, ax = plt.subplots(figsize=(12,8))
ax = data['1901':].plot(ax=ax)
predict_y.plot(ax=ax)
plt.show()