應用時間序列
時間序列分析是一種重要的資料分析方法，應用廣泛。以下列舉了幾個時間序列分析的應用場景：
1.經濟預測：時間序列分析可以用來分析經濟資料，預測未來經濟趨勢和走向。例如，利用歷史股市資料和經濟指標進行時間序列分析，可以預測未來股市的走向。
2.交通擁堵預測：時間序列分析可以用來預測交通擁堵情況。例如，根據歷史車流量和天氣資訊，可以建立一個時間序列模型來預測未來某個時間段的路況。
3.天氣預測：時間序列分析可以用於天氣預測，例如預測未來幾天或幾周的降雨量、溫度等。這對於農業生產和水資源的管理非常重要。
4.財務規劃：時間序列分析可以用來幫助企業進行財務規劃。例如，透過分析歷史銷售資料，可以預測未來銷售額，並制定相應的預算計劃。
5.工業控制：時間序列分析可以用來最佳化工業生產過程。例如，根據機器執行狀態和歷史生產資料，可以建立一個時間序列模型來最佳化生產線的執行，提高生產效率。
源資料樣例：
以下是資料的具體時間序列分析步驟：
注意：以下時間序列分析程式碼僅限一階差分或者不使用差分進行的程式碼，以一階差分為例。
一.匯入資料
1 df=pd.read_excel('M17.xlsx') 2 print(df["水位"].head()) 3 4 # 建立一個DatetimeIndex物件，指定時間範圍為2021年1月1日到1月5日 5 dates = pd.date_range(start='1860-01-01', end='1955-12-31',freq='AS') 6 # 建立一個Series物件，將列表作為資料，DatetimeIndex物件作為索引 7 t_s = pd.Series(data=df["水位"].tolist(), index=dates) 8 t_s.index = t_s.index.strftime('%Y')
注：這裡的資料df為dataframe資料，兩列；t_s為時間序列資料，第二列，它內涵時間索引
二.繪製時序圖
時序圖是一種用於展示時間序列資料的圖表，通常將時間作為X軸，將變數（如銷售額、溫度等）作為Y軸。時序圖可以幫助我們觀察和分析時間序列資料的趨勢、季節性、週期性以及異常值等。
一個典型的時序圖通常包括以下幾個元素：
1-->X軸：表示時間軸，通常根據資料的時間粒度來設定刻度。例如，如果資料是按日收集的，則X軸可能顯示日期；如果是按小時收集的，則X軸可能顯示小時數。
2-->Y軸：表示變數的取值範圍，通常根據資料的特性來設定刻度。例如，如果資料表示某個產品的銷售額，則Y軸可能顯示金額數值；如果資料表示溫度，則Y軸可能顯示攝氏度或華氏度。、
3-->資料線：表示時間序列資料的變化趨勢，通常用一條連續的曲線或折線來表示。資料線的顏色和樣式可以根據需要進行調整，以突出重點資訊。
4-->標題和註釋：用於說明時序圖的含義、資料來源和相關資訊。
將時間序列資料視覺化成時序圖，可以更加直觀地觀察和分析資料的變化趨勢和規律，從而更好地理解資料。
首先可以繪製線圖直接觀察資料走勢粗略判斷平穩性，既無趨勢也無週期
1 #時序圖 2 plt.plot(df['年份'],df['水位']) 3 plt.xlabel("年份") 4 plt.ylabel("水位") 5 plt.show()
三.初步檢驗
3.1.純隨機性檢驗（白噪聲檢驗）
3.1.1Ljung-Box 檢驗
Ljung-Box 檢驗是一種常用的時間序列分析方法，可以用於檢測資料是否具有白噪聲特徵。白噪聲指的是一種隨機性非常強的訊號，其中每個時間點的值都是獨立且服從相同的分佈。
Ljung-Box 檢驗是基於自相關和偏自相關函式計算的。在實際應用中，需要將時間序列資料轉化為平穩序列，然後計算其自相關和偏自相關函式。Ljung-Box 檢驗會統計一組拉格朗日乘子（Lagrange Multiplier）值，從而判斷序列資料是否具有白噪聲特徵。
通常情況下，如果 Ljung-Box 檢驗的統計值小於置信度水平對應的臨界值，就可以認為時間序列資料具有白噪聲特徵。反之，如果統計值大於臨界值，則可以認為時間序列資料不具有白噪聲特徵。
print("隨機性檢驗：",'\n',lb_test(df["水位"], return_df=True,lags=5))
其中 df["水位"] 是一個一維陣列，表示時間序列資料；lags 引數指定計算自相關和偏自相關函式的滯後階數。如果希望以一定的置信度進行 Ljung-Box 檢驗，可以根據自由度和置信度水平計算出對應的臨界值，然後將其與 stat 進行比較。
注：當p值小於0.05時為顯著非白噪聲序列
3.1.2單位根檢驗：ADF檢驗
單位根檢驗，也稱為增廣迪基-福勒檢驗（Augmented Dickey-Fuller test，ADF test），是一種用於檢驗時間序列資料是否具有單位根（unit root）的方法。在時間序列分析中，單位根通常是假設時間序列資料中存在一種長期的趨勢或非平穩性。
單位根檢驗旨在驗證時間序列資料是否具有平穩性。如果時間序列資料具有單位根，則說明資料存在非平穩性，並且預測和分析可能會出現問題。因此，在進行時間序列分析之前，我們需要先進行單位根檢驗，以確保資料具有平穩性。
ADF檢驗是一種常用的單位根檢驗方法，它透過計算時間序列資料的ADF統計量來判斷資料是否具有單位根。ADF統計量與t值類似，表示觀測值與滯後版本之間的差異程度，同時考慮了其他因素的影響。如果ADF統計量小於對應的臨界值，則拒絕原假設，即資料不存在單位根，表明資料具有平穩性。
除了ADF檢驗外，還有許多其他的單位根檢驗方法，例如Dickey-Fuller檢驗、KPSS檢驗等。不同的單位根檢驗方法具有不同的假設條件和適用範圍，需要根據具體情況來選擇合適的方法。
總之，單位根檢驗是一種重要的時間序列分析工具，用於驗證資料是否具有平穩性。只有在資料具有平穩性的情況下，才能進行有效的預測和分析。
檢驗假設：H0：存在單位根 vs H1：不存在單位根
如果序列平穩，則不應存在單位根，所以我們希望能拒絕原假設
1 # 進行ADF檢驗 2 result = ts.adfuller(df['水位']) 3 # 輸出ADF結果 4 print('-------------------------------------------') 5 print('ADF檢驗結果：') 6 print('ADF Statistic: %f' % result[0]) 7 print('p-value: %f' % result[1]) 8 print('Lags Used: %d' % result[2]) 9 print('Observations Used: %d' % result[3]) 10 print('Critical Values:') 11 for key, value in result[4].items(): 12 print('\t%s: %.3f' % (key, value))
當p-value>0.05，則要進行差分運算！否之直接跳轉到模型選擇與定階。
3.2差分運算
差分運算是時序資料預處理中的一個常見操作，它可以用來消除時間序列資料中的趨勢和季節性變化，從而使得資料更加平穩。在時間序列分析和建模中，平穩序列通常比非平穩序列更容易建模，並且能夠提供更精確的預測結果。
差分運算可以透過計算序列中相鄰資料之間的差值得到。一階差分是指將每個數值與其前一個數值相減，得到一個新的序列；二階差分是指對一階差分後的序列再進行一次差分運算，以此類推。
1 ##### 差分運算 2 def diff(timeseries): 3 new_df=timeseries.diff(periods=1).dropna()#dropna刪除NaN 4 new_df.plot(color='orange',title='diff1') 5 return new_df 6 7 #進行一階差分 8 ndf=diff(df['水位'])
或者如下方式也可以進行差分
ndf = pd.Series(np.diff(df['水位'], n=1))
注：n=1表示1階差分
3.3再次進行ADF檢驗
1 #再次進行ADF檢驗 2 result2 = ts.adfuller(ndf) 3 # 輸出ADF結果 4 print('-------------------------------------------') 5 print('差分後序列的ADF檢驗結果:') 6 print('ADF Statistic: %f' % result2[0]) 7 print('p-value: %f' % result2[1]) 8 print('Lags Used: %d' % result2[2]) 9 print('Observations Used: %d' % result2[3]) 10 print('Critical Values:') 11 for key, value in result2[4].items(): 12 print('\t%s: %.3f' % (key, value))
如果p-value<0.05則進行下一步定階；如p-value仍然>0.05，則繼續進行差分，再進行ADF檢驗，直到檢驗結果的p-value<0.05則進行下一步定階。
四.模型選擇與定階
4.1自相關圖和偏自相關圖人工定階
自相關圖和偏自相關圖是時序資料分析中常用的工具，可以幫助我們理解時間序列資料中的自相關性和互相關性。自相關函式（ACF）和偏自相關函式（PACF）是計算自相關圖和偏自相關圖的數學概念。
自相關圖展示了時間序列資料在不同滯後階數下的自相關係數，即前一時間點與後面所有時間點之間的相關性。自相關圖有助於判斷時間序列資料是否存在季節性或週期性變化，並且可以用來選擇合適的時間序列模型。如果一個時間序列資料存在季節性變化，則其自相關圖通常會呈現出明顯的週期性模式。
偏自相關圖展示了在滯後階數為 n 時，在其他階數的影響已被削減的情況下，上一時間點與當前時間點之間的相關性。偏自相關圖可以幫助我們確定時間序列資料中的短期相關性，從而選擇合適的時間序列模型。如果一個時間序列資料存在短期相關性，則其偏自相關圖通常會顯示出急速衰減的模式。
模型定階就是根據自相關圖和偏自相關圖來確定時間序列模型的階數。一般建議先使用自相關圖和偏自相關圖來判斷時間序列資料的性質和模型型別，然後根據其結果選擇適當的時間序列模型及其階數。常用的時間序列模型包括 AR、MA、ARMA 和 ARIMA 等。
1 #模型選擇：繪製ACF與PACF,即自相關圖和偏自相關圖 2 #### 繪製ACF與PACF的影像 3 def plot_acf_pacf(timeseries): #利用ACF和PACF判斷模型階數 4 plot_acf(timeseries,lags=timeseries.shape[0]%2) #延遲數 5 plot_pacf(timeseries,lags=timeseries.shape[0]%2) 6 plt.show() 7 plot_acf_pacf(ndf)
4.2利用AIC與BIC準則進行迭代調優
AIC（Akaike Information Criterion）和 BIC（Bayesian Information Criterion）是常用的模型選擇準則，可以用來評估不同模型在擬合資料時的複雜度和優劣程度。一般來說，AIC 和 BIC 的值越小，表示模型對資料的解釋能力越強，預測結果也更加可信。
利用 AIC 和 BIC 進行迭代調優的基本思路如下：
1-->根據需要擬合的時間序列資料，選擇一個初始模型，並計算其 AIC 和 BIC 值。
2-->對模型進行迭代調優。例如，可以逐步增加 AR 或 MA 的階數，或者使用其他方法來改進模型的效能。每次更新模型後，重新計算 AIC 和 BIC 值，並與之前的值進行比較。
3-->如果新的 AIC 或 BIC 值更小，則接受新的模型；否則，保留原始模型。重複上述過程，直到找到最小的 AIC 和 BIC 值，並確定最佳的時間序列模型及其引數。
需要注意的是，AIC 和 BIC 只是一種近似的模型評價指標。在實際應用中，還需要考慮其他因素，例如模型的可解釋性、預測誤差等。因此，在選擇時間序列模型時，需要綜合考慮多個因素，以選擇最合適的模型。
1 #迭代調優 2 print('-------------------------------------------') 3 #AIC 4 timeseries=ndf 5 AIC=sm.tsa.stattools.arma_order_select_ic(timeseries,max_ar=4,max_ma=4,ic='aic')['aic_min_order'] 6 #BIC 7 BIC=sm.tsa.stattools.arma_order_select_ic(timeseries,max_ar=4,max_ma=4,ic='bic')['bic_min_order'] 8 print('the AIC is{}\nthe BIC is{}\n'.format(AIC,BIC))
五.模型構建
ARIMA 模型是一種常用的時間序列分析模型，可以用來預測未來一段時間內的資料趨勢。ARIMA 的全稱是 Autoregressive Integrated Moving Average，其中 Autoregressive 表示自迴歸模型，Integrated 表示差分操作，Moving Average 表示滑動平均模型。
ARIMA 模型的一般形式為 ARIMA(p, d, q)，其中 p、d 和 q 分別表示 AR 模型、差分階數和 MA 模型的階數。當d=0時，即差分階數為零時為ARMA模型。具體來說：
·AR 模型是指基於前期觀測值的線性組合，來預測當前觀測值的模型。AR 模型的階數 p 表示使用前 p 個時刻的觀測值來預測當前時刻的值。
·差分操作是對原始資料進行差分運算（即相鄰兩個時間點之間的差值），以消除資料中的趨勢和季節性變化。
·MA 模型是指基於前期觀測值的隨機誤差的線性組合，來預測當前觀測值的模型。MA 模型的階數 q 表示使用前 q 個時刻的隨機誤差來預測當前時刻的值。
1 #模型構建 2 print('-------------------------------------------') 3 model= ARIMA(ndf, order=(1, 1, 2)).fit() 4 print(model.params) 5 print(model.summary())
六.模型後檢驗
6.1殘差檢驗
殘差檢驗是在統計學中經常用於檢測線性迴歸模型是否符合假設的方法之一。在進行線性迴歸時，我們基於資料擬合一個模型，並觀察模型的殘差（實際值與模型預測值之間的差異）是否具有某些特定的性質。
通常情況下，如果殘差服從正態分佈、均值為零、方差相等，則說明模型擬合良好。而如果殘差不滿足上述條件，則表明模型可能存在問題，需要進一步調整或改進。
在進行殘差檢驗時，常用的方法包括：
1-->繪製殘差圖：將樣本資料的殘差繪製成散點圖，觀察其是否隨著自變數的增加而呈現出某種規律，比如呈現出線性或非線性關係。如果殘差呈現出某種規律，則表明模型存在問題。
2-->繪製Q-Q圖：將殘差按照大小排序，並將其對應於標準正態分佈的分位數，然後繪製出散點圖。如果散點圖近似於一條直線，則表明殘差服從正態分佈。
3-->進行統計檢驗：利用統計學方法檢驗殘差是否服從正態分佈。常用的檢驗方法包括Kolmogorov-Smirnov檢驗、Shapiro-Wilk檢驗等。
總之，殘差檢驗是一種判斷線性迴歸模型擬合效果的重要方法，可以幫助我們發現和解決模型存在的問題。
殘差圖：
1 #殘差圖 2 model.resid.plot(figsize=(10,3)) 3 plt.show()
6.2Jarque-Bera檢驗
Jarque-Bera檢驗是一種用於檢驗樣本資料是否符合正態分佈的方法。它基於統計學中的JB檢驗，可以檢驗樣本資料的偏度和峰度是否符合正態分佈的假設。
在Jarque-Bera檢驗中，我們首先計算樣本資料的偏度和峰度，然後根據公式計算出JB統計量，最後對其進行顯著性檢驗。
Jarque-Bera檢驗的原假設是：樣本資料符合正態分佈。如果p值較小，則拒絕原假設，即認為樣本資料不符合正態分佈。
1 #進行Jarque-Bera檢驗 2 print('-------------------------------------------') 3 jb_test = sm.stats.stattools.jarque_bera(model.resid) 4 print('Jarque-Bera test:') 5 print('JB:', jb_test[0]) 6 print('p-value:', jb_test[1]) 7 print('Skew:', jb_test[2]) 8 print('Kurtosis:', jb_test[3])
6.3Ljung-Box檢驗
Ljung-Box檢驗是一種用於檢驗時間序列資料是否存在自相關性的方法。它基於時間序列模型的殘差，可以用來判斷殘差之間是否存在顯著的線性相關性。
在進行Ljung-Box檢驗時，我們首先計算時間序列模型的殘差，並對其進行平方處理（Q值）。然後根據公式計算出LB統計量和相應的p值，最後對其進行顯著性檢驗。
Ljung-Box檢驗的原假設是：時間序列模型的殘差之間不存在自相關性。如果p值較小，則拒絕原假設，即認為殘差之間存在顯著的自相關性。
1 print('-------------------------------------------') 2 lb_test= sm.stats.diagnostic.acorr_ljungbox(model.resid, return_df=True,lags=5) 3 print('Ljung-Box test:','\n',lb_test)
6.4DW檢驗
DW檢驗是一種用來檢驗線性迴歸模型是否存在自相關性的方法，其中DW代表Durbin-Watson。該檢驗方法通常應用於時間序列資料中的多元線性迴歸分析中，可以用於判斷誤差項（殘差）的時間序列是否具有自相關結構。
DW統計量的取值範圍為0~4，當DW統計量接近2時，表明誤差項不存在一階自相關；DW小於2，則存在正自相關；DW大於2則存在負自相關。根據經驗規律，DW統計量在1.5 ～ 2.5之間可以認為不存在自相關；小於1.5表示存在正自相關；大於2.5表示存在負自相關。
1 # 進行DW檢驗 2 print('-------------------------------------------') 3 dw_test = sm.stats.stattools.durbin_watson(model.resid) 4 print('Durbin-Watson test:') 5 print('DW:', dw_test)
6.5利用QQ圖看正態性
QQ圖，全稱Quantile-Quantile圖，是一種用於檢驗資料分佈是否符合某種理論分佈（通常為正態分佈）的方法。它透過將資料的樣本分位數與理論分位數進行比較，並將其繪製成散點圖的方式來展示資料的分佈情況。
如果樣本分佈與理論分佈一致，則散點圖應該大致呈現出一條直線。如果兩者之間存在差異，則散點圖會偏離一條直線，並呈現出“S”型或其他非線性形狀，這表明樣本資料的分佈與理論分佈不一致。
1 #QQ圖看正態性 2 qqplot(model.resid, line="q", fit=True) 3 plt.show()
七.模型評價
7.1生成模型預測資料
利用已生成的模型對資料進行重新預測，代原值計算後 透過反向差分運算（前期沒有進行差分運算則不需要，前期差分運算了幾步這裡要反向運算幾步）的資料 與 原資料（真實資料）進行比對。
1 #模型預測 2 print('-------------------------------------------') 3 predict= model.predict(1,95)#dynamic=True) 4 # print('模型預測:','\n',predict) 5 #反向差分運算 6 # 對差分後的時間序列進行逆差分運算,兩個引數：差分資料序列和原始資料序列 7 def inverse_diff(diff_series, original_series): 8 inverted = [] 9 prev = original_series.iloc[0] 10 for val in diff_series: 11 current = val + prev 12 inverted.append(current) 13 prev = current 14 return pd.Series(inverted, index=original_series.index[1:]) 15 n_predict=inverse_diff(predict,t_s) 16 print('模型預測:','\n',n_predict)
7.2真實值與預測值劃入同意影像進行比對
1 # #畫圖 2 plt.figure(figsize=(10,4)) 3 plt.plot(t_s.index,t_s,label='actual') 4 plt.plot(predict.index,n_predict,label='predict') 5 plt.xticks([]) 6 plt.legend(['actual','predict']) 7 plt.xlabel('time(year)') 8 plt.ylabel('SUNACTIVITY') 9 plt.show()
注：該圖為不進行差分運算的對比圖！！！
八.未來資料預測
進行未來某幾期的資料預測，steps=2表示要預測未來兩期, alpha=0.05表示預測的未來資料的95%的置信區間。預測結果仍需進行反向差分運算。
1 print('-------------------------------------------') 2 # 進行三步預測，並輸出95%置信區間 3 steps=3 #未來三期預測 4 forecast= model.get_forecast(steps=steps) 5 table=pd.DataFrame(forecast.summary_frame()) 6 7 # print(table.iloc[1]) 8 table.iloc[0]=table.iloc[0]+t_s[-1] 9 # print(table.iloc[0, 0]) 10 for i in range(steps-1): 11 table.iloc[i+1]=table.iloc[i+1]+table.iloc[i, 0] 12 print(table)
Python完整程式碼：
import numpy as np import pandas as pd import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.graphics.tsaplots import plot_acf,plot_pacf #ACF與PACF import statsmodels.tsa.stattools as ts #adf檢驗 from statsmodels.tsa.arima.model import ARIMA #ARIMA模型 from statsmodels.graphics.api import qqplot #qq圖 from statsmodels.stats.diagnostic import acorr_ljungbox as lb_test #用於LingBox隨機性檢驗 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False df=pd.read_excel('M17.xlsx') print(df["水位"].head()) # 建立一個DatetimeIndex物件，指定時間範圍為2021年1月1日到1月5日 dates = pd.date_range(start='1860-01-01', end='1955-12-31',freq='AS') # 建立一個Series物件，將列表作為資料，DatetimeIndex物件作為索引 t_s = pd.Series(data=df["水位"].tolist(), index=dates) t_s.index = t_s.index.strftime('%Y') #時序圖 plt.plot(df['年份'],df['水位']) plt.xlabel("年份") plt.ylabel("水位") plt.show() #純隨機性檢驗（白噪聲檢驗） #Ljung-Box 檢驗,LB檢驗，用來做純隨機性檢驗的,單位根檢驗也屬於 print('-------------------------------------------') print("隨機性檢驗：",'\n',lb_test(df["水位"], return_df=True,lags=5)) # 進行ADF檢驗 result = ts.adfuller(df['水位']) # 輸出ADF結果 print('-------------------------------------------') print('ADF檢驗結果：') print('ADF Statistic: %f' % result[0]) print('p-value: %f' % result[1]) print('Lags Used: %d' % result[2]) print('Observations Used: %d' % result[3]) print('Critical Values:') for key, value in result[4].items(): print('\t%s: %.3f' % (key, value)) ##### 差分運算 def diff(timeseries): new_df=timeseries.diff(periods=1).dropna()#dropna刪除NaN new_df.plot(color='orange',title='diff1') return new_df #進行一階差分 print("============================================================") ndf=diff(df["水位"]) print(ndf) # dif = pd.Series(np.diff(df['水位'], n=1)) # print(dif) #再次進行ADF檢驗 result2 = ts.adfuller(ndf) # 輸出ADF結果 print('-------------------------------------------') print('差分後序列的ADF檢驗結果:') print('ADF Statistic: %f' % result2[0]) print('p-value: %f' % result2[1]) print('Lags Used: %d' % result2[2]) print('Observations Used: %d' % result2[3]) print('Critical Values:') for key, value in result2[4].items(): print('\t%s: %.3f' % (key, value)) #模型選擇與定階 #模型選擇：繪製ACF與PACF,即自相關圖和偏自相關圖 #### 繪製ACF與PACF的影像 def plot_acf_pacf(timeseries): #利用ACF和PACF判斷模型階數 plot_acf(timeseries,lags=timeseries.shape[0]%2) #延遲數 plot_pacf(timeseries,lags=timeseries.shape[0]%2) plt.show() plot_acf_pacf(ndf) #迭代調優 print('-------------------------------------------') #AIC timeseries=ndf AIC=sm.tsa.stattools.arma_order_select_ic(timeseries,max_ar=4,max_ma=4,ic='aic')['aic_min_order'] #BIC BIC=sm.tsa.stattools.arma_order_select_ic(timeseries,max_ar=4,max_ma=4,ic='bic')['bic_min_order'] print('the AIC is{}\nthe BIC is{}\n'.format(AIC,BIC)) #模型構建 print('-------------------------------------------') model= ARIMA(ndf, order=(1, 1, 2)).fit() print(model.params) print(model.summary()) #殘差檢驗:檢驗殘差是否服從正態分佈，畫圖檢視，然後檢驗 #殘差圖 model.resid.plot(figsize=(10,3)) plt.show() #進行Jarque-Bera檢驗 print('-------------------------------------------') jb_test = sm.stats.stattools.jarque_bera(model.resid) print('Jarque-Bera test:') print('JB:', jb_test[0]) print('p-value:', jb_test[1]) print('Skew:', jb_test[2]) print('Kurtosis:', jb_test[3]) # 進行Ljung-Box檢驗 print('-------------------------------------------') lb_test= sm.stats.diagnostic.acorr_ljungbox(model.resid, return_df=True,lags=5) print('Ljung-Box test:','\n',lb_test) # 進行DW檢驗 print('-------------------------------------------') dw_test = sm.stats.stattools.durbin_watson(model.resid) print('Durbin-Watson test:') print('DW:', dw_test) #QQ圖看正態性 qqplot(model.resid, line="q", fit=True) plt.show() #模型預測 print('-------------------------------------------') predict= model.predict(1,95)#dynamic=True) # print('模型預測:','\n',predict) #反向差分運算 # 對差分後的時間序列進行逆差分運算,兩個引數：差分資料序列和原始資料序列 def inverse_diff(diff_series, original_series): inverted = [] prev = original_series.iloc[0] for val in diff_series: current = val + prev inverted.append(current) prev = current return pd.Series(inverted, index=original_series.index[1:]) n_predict=inverse_diff(predict,t_s) print('模型預測:','\n',n_predict) # #畫圖 plt.figure(figsize=(10,4)) plt.plot(t_s.index,t_s,label='actual') plt.plot(predict.index,n_predict,label='predict') plt.xticks([]) plt.legend(['actual','predict']) plt.xlabel('time(year)') plt.ylabel('SUNACTIVITY') plt.show() print('-------------------------------------------') # 進行三步預測，並輸出95%置信區間 steps=3 #未來三期預測 forecast= model.get_forecast(steps=steps) table=pd.DataFrame(forecast.summary_frame()) # print(table.iloc[1]) table.iloc[0]=table.iloc[0]+t_s[-1] # print(table.iloc[0, 0]) for i in range(steps-1): table.iloc[i+1]=table.iloc[i+1]+table.iloc[i, 0] print(table)
Python的完整包裝之後的庫:
1 import numpy as np 2 import pandas as pd 3 import matplotlib.pyplot as plt 4 import statsmodels.api as sm 5 from statsmodels.graphics.tsaplots import plot_acf,plot_pacf #ACF與PACF 6 import statsmodels.tsa.stattools as ts #adf檢驗 7 from statsmodels.tsa.arima.model import ARIMA #ARIMA模型 8 from statsmodels.graphics.api import qqplot #qq圖 9 from statsmodels.stats.diagnostic import acorr_ljungbox as lb_test #用於LingBox隨機性檢驗 10 11 plt.rcParams['font.sans-serif'] = ['SimHei'] 12 plt.rcParams['axes.unicode_minus'] = False 13 14 class Time_ARIMA(): 15 def __init__(self,path): 16 try: 17 self.df = pd.read_csv(path) 18 except: 19 self.df=pd.read_excel(path) 20 # print(self.df.iloc[:,0]) 21 start=str(self.df.iloc[0,0])+'-01-01' 22 end=str(self.df.iloc[-1,0])+'-01-01' 23 # 建立一個DatetimeIndex物件，指定時間範圍為2021年1月1日到1月5日 24 dates = pd.date_range(start=start, end=end, freq='AS') 25 # 建立一個Series物件，將列表作為資料，DatetimeIndex物件作為索引 26 self.t_s = pd.Series(data=self.df.iloc[:,1].tolist(), index=dates) 27 self.t_s.index = self.t_s.index.strftime('%Y') 28 # print(self.t_s) 29 self.timeseries=None 30 self.flag=0 31 self.BIC=None 32 self.model=None 33 34 def Timing_Diagram(self): 35 # 時序圖 36 plt.plot(self.t_s) 37 plt.xlabel(self.df.iloc[:,0].name) 38 plt.xticks(np.arange(0,self.df.shape[0], step=10)) #設定橫座標間距 39 plt.ylabel(self.df.iloc[:,1].name) 40 plt.show() 41 42 def Diff_pre_sequence_test(self): 43 # 純隨機性檢驗（白噪聲檢驗） 44 # Ljung-Box 檢驗,LB檢驗，用來做純隨機性檢驗的,單位根檢驗也屬於 45 print('==================原始序列的純隨機性檢驗==================') 46 print("Ljung-Box檢驗：", '\n', lb_test(self.df.iloc[:,1], return_df=True, lags=5)) 47 # 進行ADF檢驗 48 result = ts.adfuller(self.df.iloc[:,1]) 49 # 輸出ADF結果 50 print('---------------------ADF檢驗結果----------------------') 51 print('ADF Statistic: %f' % result[0]) 52 print('p-value: %f' % result[1]) 53 print('Lags Used: %d' % result[2]) 54 print('Observations Used: %d' % result[3]) 55 print('Critical Values:') 56 for key, value in result[4].items(): 57 print('\t%s: %.3f' % (key, value)) 58 print('====================================================') 59 60 def Diff_operation(self,timeseries=None): 61 try: 62 #二階及以上差分 63 ndf = timeseries.diff(periods=1).dropna() # dropna刪除NaN 64 except: 65 #一階差分 66 ndf = self.t_s.diff(periods=1).dropna() # dropna刪除NaN 67 ndf.plot(color='orange', title='殘差圖') 68 self.flag+=1 69 self.timeseries=ndf 70 return ndf 71 72 def Diff_after_sequence_testing(self): 73 # 再次進行ADF檢驗 74 if self.flag==0: 75 print("未進行差分運算，不必再次進行ADF檢驗") 76 else: 77 result2 = ts.adfuller(self.timeseries) 78 # 輸出ADF結果 79 print('------------差分後序列的ADF檢驗結果------------') 80 print('ADF Statistic: %f' % result2[0]) 81 print('p-value: %f' % result2[1]) 82 print('Lags Used: %d' % result2[2]) 83 print('Observations Used: %d' % result2[3]) 84 print('Critical Values:') 85 for key, value in result2[4].items(): 86 print('\t%s: %.3f' % (key, value)) 87 88 89 def ACF_and_PACF(self): 90 # 繪製ACF與PACF的影像,利用ACF和PACF判斷模型階數 91 # print(int(self.t_s.shape[0]/2)) 92 if self.flag == 0: 93 plot_acf(self.t_s, lags=int(self.t_s.shape[0] / 2) - 1) # 延遲數 94 plot_pacf(self.t_s, lags=int(self.t_s.shape[0] / 2) - 1) 95 plt.show() 96 else: 97 plot_acf(self.timeseries, lags=int(pd.Series(self.timeseries).shape[0] / 2 - 1)) # 延遲數 98 plot_pacf(self.timeseries, lags=int(pd.Series(self.timeseries).shape[0] / 2 - 1)) 99 plt.show() 100 101 def Model_order_determination(self): 102 # 迭代調優 103 if self.flag==0: 104 timeseries=self.t_s 105 else: 106 timeseries=self.timeseries 107 # AIC 108 AIC = sm.tsa.stattools.arma_order_select_ic(timeseries, max_ar=4, max_ma=4, ic='aic')['aic_min_order'] 109 # BIC 110 BIC = sm.tsa.stattools.arma_order_select_ic(timeseries, max_ar=4, max_ma=4, ic='bic')['bic_min_order'] 111 print('---AIC與BIC準則定階---') 112 print('the AIC is{}\nthe BIC is{}\n'.format(AIC, BIC),end='') 113 print('--------------------') 114 self.BIC=BIC 115 116 def Model_building_machine(self): 117 if self.flag==0: 118 timeseries=self.t_s 119 else: 120 timeseries=self.timeseries 121 model = ARIMA(timeseries, order=(self.BIC[0], self.flag, self.BIC[1])).fit() 122 print('--------引數估計值-------') 123 print(model.params) 124 print(model.summary()) 125 self.model=model 126 127 def Model_building_artificial(self,order=(0,0,0)): 128 if self.flag==0: 129 timeseries=self.t_s 130 else: 131 timeseries=self.timeseries 132 model = ARIMA(timeseries, order=order).fit() 133 print('--------引數估計值-------') 134 print(model.params) 135 print(model.summary()) 136 self.model = model 137 138 def Model_checking(self): 139 # 殘差檢驗:檢驗殘差是否服從正態分佈，畫圖檢視，然後檢驗 140 # 殘差圖 141 self.model.resid.plot(figsize=(10, 3)) 142 plt.title("殘差圖") 143 plt.show() 144 145 # 進行Jarque-Bera檢驗 146 jb_test = sm.stats.stattools.jarque_bera(self.model.resid) 147 print("==================================================") 148 print('------------Jarque-Bera檢驗-----------') 149 print('Jarque-Bera test:') 150 print('JB:', jb_test[0]) 151 print('p-value:', jb_test[1]) 152 print('Skew:', jb_test[2]) 153 print('Kurtosis:', jb_test[3]) 154 155 # 進行Ljung-Box檢驗 156 lb_test = sm.stats.diagnostic.acorr_ljungbox(self.model.resid, return_df=True, lags=5) 157 print('------------Ljung-Box檢驗-------------') 158 print('Ljung-Box test:', '\n', lb_test) 159 160 # 進行DW檢驗 161 dw_test = sm.stats.stattools.durbin_watson(self.model.resid) 162 print('---------------DW檢驗----------------') 163 print('Durbin-Watson test:') 164 print('DW:', dw_test) 165 print("==================================================") 166 167 # QQ圖看正態性 168 qqplot(self.model.resid, line="q", fit=True) 169 plt.title("Q-Q圖") 170 plt.show() 171 172 def Model_prediction_accuracy(self): 173 # 模型預測 174 if self.flag==0: 175 n_predict = self.model.predict(1, self.t_s.shape[0]) 176 n_predict.index = n_predict.index.strftime('%Y') 177 print('模型預測:', '\n', n_predict) 178 179 else: 180 # 反向差分運算 181 # 對差分後的時間序列進行逆差分運算,兩個引數：差分資料序列和原始資料序列 182 def inverse_diff(diff_series, original_series): 183 inverted = [] 184 prev = original_series.iloc[0] 185 for val in diff_series: 186 current = val + prev 187 inverted.append(current) 188 prev = current 189 return pd.Series(inverted, index=original_series.index[1:]) 190 191 predict = self.model.predict(1, pd.Series(self.timeseries).shape[0]) 192 n_predict = inverse_diff(predict, self.t_s) 193 print('模型預測:', '\n',n_predict) 194 print('-------------------------------------------') 195 196 # #畫圖 197 plt.figure(figsize=(10, 4)) 198 plt.plot(self.t_s.index, self.t_s, label='actual') 199 plt.plot(n_predict.index, n_predict, label='predict') 200 # plt.xticks([]) 201 plt.xticks(np.arange(0, self.df.shape[0], step=10)) # 設定橫座標間距 202 plt.legend(['actual', 'predict']) 203 plt.xlabel('time(year)') 204 plt.ylabel('SUNACTIVITY') 205 plt.title("真實值與預測值對比") 206 plt.show() 207 208 def Future_prediction(self,steps=2): 209 # 進行三步預測，並輸出95%置信區間 210 211 # steps = 3 # 未來三期預測 212 forecast = self.model.get_forecast(steps=steps) 213 table = pd.DataFrame(forecast.summary_frame()) 214 215 # print(table.iloc[1]) 216 table.iloc[0] = table.iloc[0] + self.t_s[-1] 217 # print(table.iloc[0, 0]) 218 for i in range(steps - 1): 219 table.iloc[i + 1] = table.iloc[i + 1] + table.iloc[i, 0] 220 table.index = table.index.strftime('%Y') 221 print('--------------------------------------------------------------') 222 print(table) 223 print('--------------------------------------------------------------')
對應的呼叫方式：
1 import TARMA 2 ts=TARMA.Time_ARIMA(path='./M17.xlsx') #例項化物件（放入檔案路徑） 3 ts.Timing_Diagram() #時序圖 4 ts.Diff_pre_sequence_test() #純隨機性檢驗 5 6 ts.ACF_and_PACF() #自相關圖與偏自相關圖 7 ts.Model_order_determination() #根據AIC準則與BIC準則進行機器定階 8 9 #下面兩種定階方式二選一 10 # ts.Model_building_artificial(order=(1,1,2)) #模型構建(人工定階),order中的三個引數為p，d，q 11 ts.Model_building_machine() #模型構建(基於機器自己定階) 12 13 ts.Model_checking() #模型檢驗 14 ts.Model_prediction_accuracy() #模型可行性預測 15 ts.Future_prediction(steps=3) #未來資料預測預設為兩期，steps引數可自定義
1 import TARMA 2 ts=TARMA.Time_ARIMA(path='./M17.xlsx') #例項化物件（放入檔案路徑） 3 ts.Timing_Diagram() #時序圖 4 ts.Diff_pre_sequence_test() #純隨機性檢驗 5 # 進行一階差分 6 ts.Diff_operation() #進行一階差分，當ADF檢驗合適即平穩且非白噪聲序列不進行差分 7 ts.Diff_after_sequence_testing() #差分運算後的ADF檢驗,未進行差分運算則不呼叫 8 ts.ACF_and_PACF() #自相關圖與偏自相關圖 9 ts.Model_order_determination() #根據AIC準則與BIC準則進行機器定階 10 11 #下面兩種定階方式二選一 12 # ts.Model_building_artificial(order=(1,1,2)) #模型構建(人工定階),order中的三個引數為p，d，q 13 ts.Model_building_machine() #模型構建(基於機器自己定階) 14 15 ts.Model_checking() #模型檢驗 16 ts.Model_prediction_accuracy() #模型可行性預測 17 ts.Future_prediction(steps=3) #未來資料預測預設為兩期，steps引數可自定義