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檢驗(yàn)
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隨機(jī)性檢驗(yàn)

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

df=pd.read_excel('M17.xlsx')
print(df["水位"].head())

# 創(chuàng)建一個(gè)DatetimeIndex對象，指定時(shí)間范圍為2021年1月1日到1月5日
dates = pd.date_range(start='1860-01-01', end='1955-12-31',freq='AS')
# 創(chuàng)建一個(gè)Series對象，將列表作為數(shù)據(jù)，DatetimeIndex對象作為索引
t_s = pd.Series(data=df["水位"].tolist(), index=dates)
t_s.index = t_s.index.strftime('%Y')

#時(shí)序圖
plt.plot(df['年份'],df['水位'])
plt.xlabel("年份")
plt.ylabel("水位")
plt.show()

#純隨機(jī)性檢驗(yàn)（白噪聲檢驗(yàn)）
#Ljung-Box 檢驗(yàn),LB檢驗(yàn)，用來做純隨機(jī)性檢驗(yàn)的,單位根檢驗(yàn)也屬于
print('-------------------------------------------')
print("隨機(jī)性檢驗(yàn)：",'\n',lb_test(df["水位"], return_df=True,lags=5))
# 進(jìn)行ADF檢驗(yàn)
result = ts.adfuller(df['水位'])
# 輸出ADF結(jié)果
print('-------------------------------------------')
print('ADF檢驗(yàn)結(jié)果：')
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))

##### 差分運(yùn)算
def diff(timeseries):
    new_df=timeseries.diff(periods=1).dropna()#dropna刪除NaN
    new_df.plot(color='orange',title='diff1')
    return new_df
#進(jìn)行一階差分
print("============================================================")
ndf=diff(df["水位"])
print(ndf)
# dif = pd.Series(np.diff(df['水位'], n=1))
# print(dif)
#再次進(jìn)行ADF檢驗(yàn)
result2 = ts.adfuller(ndf)
# 輸出ADF結(jié)果
print('-------------------------------------------')
print('差分后序列的ADF檢驗(yàn)結(jié)果:')
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,即自相關(guān)圖和偏自相關(guān)圖
#### 繪制ACF與PACF的圖像
def plot_acf_pacf(timeseries): #利用ACF和PACF判斷模型階數(shù)
    plot_acf(timeseries,lags=timeseries.shape[0]%2) #延遲數(shù)
    plot_pacf(timeseries,lags=timeseries.shape[0]%2)
    plt.show()
plot_acf_pacf(ndf)

#迭代調(diào)優(yōu)
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))

#模型構(gòu)建
print('-------------------------------------------')
model= ARIMA(ndf, order=(1, 1, 2)).fit()
print(model.params)
print(model.summary())

#殘差檢驗(yàn):檢驗(yàn)殘差是否服從正態(tài)分布，畫圖查看，然后檢驗(yàn)
#殘差圖
model.resid.plot(figsize=(10,3))
plt.show()
#進(jìn)行Jarque-Bera檢驗(yàn)
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])

# 進(jìn)行Ljung-Box檢驗(yàn)
print('-------------------------------------------')
lb_test= sm.stats.diagnostic.acorr_ljungbox(model.resid, return_df=True,lags=5)
print('Ljung-Box test:','\n',lb_test)

# 進(jìn)行DW檢驗(yàn)
print('-------------------------------------------')
dw_test = sm.stats.stattools.durbin_watson(model.resid)
print('Durbin-Watson test:')
print('DW:', dw_test)

#QQ圖看正態(tài)性
qqplot(model.resid, line="q", fit=True)
plt.show()

#模型預(yù)測
print('-------------------------------------------')
predict= model.predict(1,95)#dynamic=True)
# print('模型預(yù)測:','\n',predict)
#反向差分運(yùn)算
# 對差分后的時(shí)間序列進(jìn)行逆差分運(yùn)算,兩個(gè)參數(shù)：差分?jǐn)?shù)據(jù)序列和原始數(shù)據(jù)序列
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('模型預(yù)測:','\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('-------------------------------------------')
# 進(jìn)行三步預(yù)測，并輸出95%置信區(qū)間
steps=3   #未來三期預(yù)測
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)

  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檢驗(yàn)
  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隨機(jī)性檢驗(yàn)
 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         # 創(chuàng)建一個(gè)DatetimeIndex對象，指定時(shí)間范圍為2021年1月1日到1月5日
 24         dates = pd.date_range(start=start, end=end, freq='AS')
 25         # 創(chuàng)建一個(gè)Series對象，將列表作為數(shù)據(jù)，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         # 時(shí)序圖
 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))   #設(shè)置橫坐標(biāo)間距
 39         plt.ylabel(self.df.iloc[:,1].name)
 40         plt.show()
 41 
 42     def Diff_pre_sequence_test(self):
 43         # 純隨機(jī)性檢驗(yàn)（白噪聲檢驗(yàn)）
 44         # Ljung-Box 檢驗(yàn),LB檢驗(yàn)，用來做純隨機(jī)性檢驗(yàn)的,單位根檢驗(yàn)也屬于
 45         print('==================原始序列的純隨機(jī)性檢驗(yàn)==================')
 46         print("Ljung-Box檢驗(yàn)：", '\n', lb_test(self.df.iloc[:,1], return_df=True, lags=5))
 47         # 進(jìn)行ADF檢驗(yàn)
 48         result = ts.adfuller(self.df.iloc[:,1])
 49         # 輸出ADF結(jié)果
 50         print('---------------------ADF檢驗(yàn)結(jié)果----------------------')
 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         # 再次進(jìn)行ADF檢驗(yàn)
 74         if self.flag==0:
 75             print("未進(jìn)行差分運(yùn)算，不必再次進(jìn)行ADF檢驗(yàn)")
 76         else:
 77             result2 = ts.adfuller(self.timeseries)
 78             # 輸出ADF結(jié)果
 79             print('------------差分后序列的ADF檢驗(yàn)結(jié)果------------')
 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判斷模型階數(shù)
 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)  # 延遲數(shù)
 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))  # 延遲數(shù)
 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         # 迭代調(diào)優(yōu)
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準(zhǔn)則定階---')
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(self.t_s, order=(self.BIC[0], self.flag, self.BIC[1])).fit()
122         print('--------參數(shù)估計(jì)值-------')
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(self.t_s, order=order).fit()
133         print('--------參數(shù)估計(jì)值-------')
134         print(model.params)
135         print(model.summary())
136         self.model = model
137 
138     def Model_checking(self):
139         # 殘差檢驗(yàn):檢驗(yàn)殘差是否服從正態(tài)分布，畫圖查看，然后檢驗(yàn)
140         # 殘差圖
141         self.model.resid.plot(figsize=(10, 3))
142         plt.title("殘差圖")
143         plt.show()
144 
145         # 進(jìn)行Jarque-Bera檢驗(yàn)
146         jb_test = sm.stats.stattools.jarque_bera(self.model.resid)
147         print("==================================================")
148         print('------------Jarque-Bera檢驗(yàn)-----------')
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         # 進(jìn)行Ljung-Box檢驗(yàn)
156         lb_test = sm.stats.diagnostic.acorr_ljungbox(self.model.resid, return_df=True, lags=5)
157         print('------------Ljung-Box檢驗(yàn)-------------')
158         print('Ljung-Box test:', '\n', lb_test)
159 
160         # 進(jìn)行DW檢驗(yàn)
161         dw_test = sm.stats.stattools.durbin_watson(self.model.resid)
162         print('---------------DW檢驗(yàn)----------------')
163         print('Durbin-Watson test:')
164         print('DW:', dw_test)
165         print("==================================================")
166 
167         # QQ圖看正態(tài)性
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         # 模型預(yù)測
174         n_predict = self.model.predict(1, self.t_s.shape[0])
175         n_predict.index = n_predict.index.strftime('%Y')
176         print('模型預(yù)測:', '\n', n_predict)
177         # if self.flag==0:
178         #     n_predict = self.model.predict(1, self.t_s.shape[0])
179         #     n_predict.index = n_predict.index.strftime('%Y')
180         #     print('模型預(yù)測:', '\n', n_predict)
181         #
182         # else:
183         #     # 反向差分運(yùn)算
184         #     # 對差分后的時(shí)間序列進(jìn)行逆差分運(yùn)算,兩個(gè)參數(shù)：差分?jǐn)?shù)據(jù)序列和原始數(shù)據(jù)序列
185         #     def inverse_diff(diff_series, original_series):
186         #         inverted = []
187         #         prev = original_series.iloc[0]
188         #         for val in diff_series:
189         #             current = val + prev
190         #             inverted.append(current)
191         #             prev = current
192         #         return pd.Series(inverted, index=original_series.index[1:])
193         #
194         #     predict = self.model.predict(1, pd.Series(self.timeseries).shape[0])
195         #     n_predict = inverse_diff(predict, self.t_s)
196         #     print('模型預(yù)測:', '\n',n_predict)
197         print('-------------------------------------------')
198 
199         # #畫圖
200         plt.figure(figsize=(10, 4))
201         plt.plot(self.t_s.index, self.t_s, label='actual')
202         plt.plot(n_predict.index, n_predict, label='predict')
203         # plt.xticks([])
204         plt.xticks(np.arange(0, self.df.shape[0], step=10))  # 設(shè)置橫坐標(biāo)間距
205         plt.legend(['actual', 'predict'])
206         plt.xlabel('time(year)')
207         plt.ylabel('SUNACTIVITY')
208         plt.title("真實(shí)值與預(yù)測值對比")
209         plt.show()
210 
211     def Future_prediction(self,steps=2):
212         # 進(jìn)行三步預(yù)測，并輸出95%置信區(qū)間
213 
214         # steps = 3  # 未來三期預(yù)測
215         forecast = self.model.get_forecast(steps=steps)
216         table = pd.DataFrame(forecast.summary_frame())
217 
218         # # print(table.iloc[1])
219         # table.iloc[0] = table.iloc[0] + self.t_s[-1]
220         # # print(table.iloc[0, 0])
221         # for i in range(steps - 1):
222         #     table.iloc[i + 1] = table.iloc[i + 1] + table.iloc[i, 0]
223         table.index = table.index.strftime('%Y')
224         print('--------------------------------------------------------------')
225         print(table)
226         print('--------------------------------------------------------------')

 1 import TARMA
 2 #不含差分
 3 ts=TARMA.Time_ARIMA(path='./M17.xlsx')   #實(shí)例化對象（放入文件路徑）
 4 ts.Timing_Diagram()   #時(shí)序圖
 5 ts.Diff_pre_sequence_test()   #純隨機(jī)性檢驗(yàn)
 6 
 7 ts.ACF_and_PACF()   #自相關(guān)圖與偏自相關(guān)圖
 8 ts.Model_order_determination()   #根據(jù)AIC準(zhǔn)則與BIC準(zhǔn)則進(jìn)行機(jī)器定階
 9 
10 #下面兩種定階方式二選一
11 # ts.Model_building_artificial(order=(1,1,2))   #模型構(gòu)建(人工定階),order中的三個(gè)參數(shù)為p，d，q
12 ts.Model_building_machine()   #模型構(gòu)建(基于機(jī)器自己定階)
13 
14 ts.Model_checking()   #模型檢驗(yàn)
15 ts.Model_prediction_accuracy()   #模型可行性預(yù)測
16 ts.Future_prediction(steps=3)    #未來數(shù)據(jù)預(yù)測默認(rèn)為兩期，steps參數(shù)可自定義

 1 import TARMA
 2 ts=TARMA.Time_ARIMA(path='./M17.xlsx')   #實(shí)例化對象（放入文件路徑）
 3 ts.Timing_Diagram()   #時(shí)序圖
 4 ts.Diff_pre_sequence_test()   #純隨機(jī)性檢驗(yàn)
 5 # 進(jìn)行一階差分（需要多階差分多次調(diào)用差分即可（模型會(huì)隨差分次數(shù)隨之改變，不必理會(huì)）
 6 ts.Diff_operation()  #進(jìn)行一階差分，當(dāng)ADF檢驗(yàn)合適即平穩(wěn)且非白噪聲序列不進(jìn)行差分
 7 ts.Diff_after_sequence_testing()  #差分運(yùn)算后的ADF檢驗(yàn),未進(jìn)行差分運(yùn)算則不調(diào)用
 8 # #進(jìn)行二階差分
 9 # ts.Diff_operation()
10 # ts.Diff_after_sequence_testing()
11 
12 ts.ACF_and_PACF()   #自相關(guān)圖與偏自相關(guān)圖
13 ts.Model_order_determination()   #根據(jù)AIC準(zhǔn)則與BIC準(zhǔn)則進(jìn)行機(jī)器定階
14 
15 #下面兩種定階方式二選一
16 # ts.Model_building_artificial(order=(1,1,2))   #模型構(gòu)建(人工定階),order中的三個(gè)參數(shù)為p，d，q
17 ts.Model_building_machine()   #模型構(gòu)建(基于機(jī)器自己定階)
18 
19 ts.Model_checking()   #模型檢驗(yàn)
20 ts.Model_prediction_accuracy()   #模型可行性預(yù)測
21 ts.Future_prediction(steps=3)    #未來數(shù)據(jù)預(yù)測默認(rèn)為兩期，steps參數(shù)可自定義