sarima模型

以下內容引自：https://blog.csdn.net/qifeidemumu/article/details/88782550

使用“網格搜索”來迭代地探索參數的不同組合。 對於參數的每個組合，我們使用statsmodels模塊的SARIMAX()函數擬合一個新的季節性ARIMA模型，並評估其整體質量。 一旦我們探索了參數的整個范圍，我們的最佳參數集將是我們感興趣的標准產生最佳性能的參數。 我們開始生成我們希望評估的各種參數組合：

1.  
   # Define the p, d and q parameters to take any value between 0 and 2
2.  
   p = d = q = range(0, 2)
3.  
    
4.  
   # Generate all different combinations of p, q and q triplets
5.  
   pdq = list(itertools.product(p, d, q))
6.  
    
7.  
   # Generate all different combinations of seasonal p, q and q triplets
8.  
   seasonal_pdq = [( x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
9.  
    
10.  
    print('Examples of parameter combinations for Seasonal ARIMA...')
11.  
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
12.  
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
13.  
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
14.  
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

1.  
   Output
2.  
   Examples of parameter combinations for Seasonal ARIMA...
3.  
   SARIMAX: ( 0, 0, 1) x (0, 0, 1, 12)
4.  
   SARIMAX: ( 0, 0, 1) x (0, 1, 0, 12)
5.  
   SARIMAX: ( 0, 1, 0) x (0, 1, 1, 12)
6.  
   SARIMAX: ( 0, 1, 0) x (1, 0, 0, 12)

我們現在可以使用上面定義的參數三元組來自動化不同組合對ARIMA模型進行培訓和評估的過程。 在統計和機器學習中，這個過程被稱為模型選擇的網格搜索（或超參數優化）。

在評估和比較配備不同參數的統計模型時，可以根據數據的適合性或准確預測未來數據點的能力，對每個參數進行排序。 我們將使用AIC （Akaike信息標准）值，該值通過使用statsmodels安裝的ARIMA型號方便地返回。 AIC衡量模型如何適應數據，同時考慮到模型的整體復雜性。 在使用大量功能的情況下，適合數據的模型將被賦予比使用較少特征以獲得相同的適合度的模型更大的AIC得分。 因此，我們有興趣找到產生最低AIC值的模型。

下面的代碼塊通過參數的組合來迭代，並使用SARIMAX函數來適應相應的季節性ARIMA模型。 這里， order參數指定(p, d, q)參數，而seasonal_order參數指定季節性ARIMA模型的(P, D, Q, S)季節分量。 在安裝每個SARIMAX()模型后，代碼打印出其各自的AIC得分。

1.  
   warnings.filterwarnings( "ignore") # specify to ignore warning messages
2.  
    
3.  
   for param in pdq:
4.  
   for param_seasonal in seasonal_pdq:
5.  
   try:
6.  
   mod = sm.tsa.statespace.SARIMAX(y,
7.  
   order=param,
8.  
   seasonal_order=param_seasonal,
9.  
   enforce_stationarity= False,
10.  
    enforce_invertibility= False)
11.  
     
12.  
    results = mod.fit()
13.  
     
14.  
    print( 'ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
15.  
    except:
16.  
    continue

由於某些參數組合可能導致數字錯誤指定，因此我們明確禁用警告消息，以避免警告消息過載。 這些錯誤指定也可能導致錯誤並引發異常，因此我們確保捕獲這些異常並忽略導致這些問題的參數組合。

上面的代碼應該產生以下結果，這可能需要一些時間：

1.  
   Output
2.  
   SARIMAX(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.3436240402125
3.  
   SARIMAX(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172764114
4.  
   SARIMAX(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
5.  
   SARIMAX(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878315690562
6.  
   SARIMAX(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978064144
7.  
   SARIMAX(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912940095
8.  
   ...
9.  
   ...
10.  
    ...
11.  
    SARIMAX(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112294245
12.  
    SARIMAX(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049123596742
13.  
    SARIMAX(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865161305
14.  
    SARIMAX(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7801413828764

我們的代碼的輸出表明， SARIMAX(1, 1, 1)x(1, 1, 1, 12)產生最低的AIC值為277.78。 因此，我們認為這是我們考慮過的所有模型中的最佳選擇。

第5步 - 安裝ARIMA時間序列模型

使用網格搜索，我們已經確定了為我們的時間序列數據生成最佳擬合模型的參數集。 我們可以更深入地分析這個特定的模型。

我們首先將最佳參數值插入到新的SARIMAX模型中：

1.  
   mod = sm.tsa.statespace.SARIMAX(y,
2.  
   order=(1, 1, 1),
3.  
   seasonal_order=( 1, 1, 1, 12),
4.  
   enforce_stationarity= False,
5.  
   enforce_invertibility= False)
6.  
    
7.  
   results = mod.fit()
8.  
    
9.  
   print(results.summary().tables[ 1]) #詳細輸出，results.summary()可以輸出全部的模型計算參數表

1.  
   Output
2.  
   ==============================================================================
3.  
   coef std err z P>|z| [0.025 0.975]
4.  
   ------------------------------------------------------------------------------
5.  
   ar.L1 0.3182 0.092 3.443 0.001 0.137 0.499
6.  
   ma.L1 -0.6255 0.077 -8.165 0.000 -0.776 -0.475
7.  
   ar.S.L12 0.0010 0.001 1.732 0.083 -0.000 0.002
8.  
   ma.S.L12 -0.8769 0.026 -33.811 0.000 -0.928 -0.826
9.  
   sigma2 0.0972 0.004 22.634 0.000 0.089 0.106
10.  
    ==============================================================================

由SARIMAX的輸出產生的SARIMAX返回大量的信息，但是我們將注意力集中在系數表上。 coef列顯示每個特征的重量（即重要性）以及每個特征如何影響時間序列。 P>|z| 列通知我們每個特征重量的意義。 這里，每個重量的p值都低於或接近0.05 ，所以在我們的模型中保留所有權重是合理的。

在 fit 季節性ARIMA模型（以及任何其他模型）的情況下，運行模型診斷是非常重要的，以確保沒有違反模型的假設。 plot_diagnostics對象允許我們快速生成模型診斷並調查任何異常行為。

1.  
   results.plot_diagnostics(figsize=(15, 12))
2.  
   plt.show()

我們的主要關切是確保我們的模型的殘差是不相關的，並且平均分布為零。 如果季節性ARIMA模型不能滿足這些特性，這是一個很好的跡象，可以進一步改善。

殘差在數理統計中是指實際觀察值與估計值（擬合值）之間的差。“殘差”蘊含了有關模型基本假設的重要信息。如果回歸模型正確的話， 我們可以將殘差看作誤差的觀測值。 它應符合模型的假設條件，且具有誤差的一些性質。利用殘差所提供的信息，來考察模型假設的合理性及數據的可靠性稱為殘差分析

在這種情況下，我們的模型診斷表明，模型殘差正常分布如下：

• 在右上圖中，我們看到紅色KDE線與N(0,1)行（其中N(0,1) ）是正態分布的標准符號，平均值0 ，標准偏差為1 ） 。 這是殘留物正常分布的良好指示。

• 左下角的qq圖顯示，殘差（藍點）的有序分布遵循采用N(0, 1)的標准正態分布采樣的線性趨勢。 同樣，這是殘留物正常分布的強烈指示。

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

這些觀察結果使我們得出結論，我們的模型選擇了令人滿意，可以幫助我們了解我們的時間序列數據和預測未來價值。

雖然我們有一個令人滿意的結果，我們的季節性ARIMA模型的一些參數可以改變，以改善我們的模型擬合。 例如，我們的網格搜索只考慮了一組受限制的參數組合，所以如果我們拓寬網格搜索，我們可能會找到更好的模型。

第6步 - 驗證預測

我們已經獲得了我們時間序列的模型，現在可以用來產生預測。 我們首先將預測值與時間序列的實際值進行比較，這將有助於我們了解我們的預測的准確性。 get_prediction()和conf_int()屬性允許我們獲得時間序列預測的值和相關的置信區間。

1.  
   pred = results.get_prediction(start=pd.to_datetime( '1998-01-01'), dynamic=False)#預測值
2.  
   pred_ci = pred.conf_int() #置信區間

上述規定需要從1998年1月開始進行預測。

dynamic=False參數確保我們每個點的預測將使用之前的所有歷史觀測。

我們可以繪制二氧化碳時間序列的實際值和預測值，以評估我們做得如何。 注意我們如何在時間序列的末尾放大日期索引。

1.  
   ax = y[ '1990':].plot(label='observed')
2.  
   pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)
3.  
    
4.  
   ax.fill_between(pred_ci. index,
5.  
   pred_ci.iloc[:, 0],
6.  
   pred_ci.iloc[:, 1], color='k', alpha=.2)
7.  
   #圖形填充fill fill_between，參考網址：
8.  
   #https: //www.cnblogs.com/gengyi/p/9416845.html
9.  
    
10.  
    ax.set_xlabel( 'Date')
11.  
    ax.set_ylabel( 'CO2 Levels')
12.  
    plt.legend()
13.  
     
14.  
    plt.show()

 pred.predicted_mean 可以得到預測均值，就是置信區間的上下加和除以2

函數間區域填充 fill_between用法：

plt.fill_between( x, y1, y2=0, where=None, interpolate=False, step=None, hold=None,data=None, **kwarg)

x - array( length N) 定義曲線的 x 坐標

y1 - array( length N ) or scalar 定義第一條曲線的 y 坐標

y2 - array( length N )  or scalar 定義第二條曲線的 y 坐標

where - array of bool (length N), optional, default: None 排除一些（垂直）區域被填充。注：我理解的垂直區域，但幫助文檔上寫的是horizontal regions

也可簡單地描述為：

plt.fill_between(x，y1，y2，where=條件表達式, color=顏色，alpha=透明度)" where = " 可以省略，直接寫條件表達式 

總體而言，我們的預測與真實價值保持一致，呈現總體增長趨勢。

量化預測的准確性

這是很有必要的。 我們將使用MSE（均方誤差，也就是方差），它總結了我們的預測的平均誤差。 對於每個預測值，我們計算其到真實值的距離並對結果求平方。 結果需要平方，以便當我們計算總體平均值時，正/負差異不會相互抵消。

1.  
   y_forecasted = pred.predicted_mean
2.  
   y_truth = y[ '1998-01-01':]
3.  
    
4.  
   # Compute the mean square error
5.  
   mse = ((y_forecasted - y_truth) ** 2).mean()
6.  
   print( 'The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

1.  
   Output
2.  
   The Mean Squared Error of our forecasts is 0.07

我們前進一步預測的MSE值為0.07 ，這是接近0的非常低的值。MSE=0是預測情況將是完美的精度預測參數的觀測值，但是在理想的場景下通常不可能。

然而，使用動態預測(dynamic=True)可以獲得更好地表達我們的真實預測能力。 在這種情況下，我們只使用時間序列中的信息到某一點，之后，使用先前預測時間點的值生成預測。

在下面的代碼塊中，我們指定從1998年1月起開始計算動態預測和置信區間。

1.  
   pred_dynamic = results.get_prediction( start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
2.  
   pred_dynamic_ci = pred_dynamic.conf_int()

繪制時間序列的觀測值和預測值，我們看到即使使用動態預測，總體預測也是准確的。 所有預測值（紅線）與地面真相（藍線）相當吻合，並且在我們預測的置信區間內。

1.  
   ax = y[ '1990':].plot(label='observed', figsize=(20, 15))
2.  
   pred_dynamic.predicted_mean.plot( label='Dynamic Forecast', ax=ax)
3.  
    
4.  
   ax.fill_between(pred_dynamic_ci. index,
5.  
   pred_dynamic_ci.iloc[:, 0],
6.  
   pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
7.  
    
8.  
   ax.fill_betweenx(ax.get_ylim(), pd.to_datetime( '1998-01-01'), y.index[-1],
9.  
   alpha=. 1, zorder=-1)
10.  
     
11.  
    ax.set_xlabel( 'Date')
12.  
    ax.set_ylabel( 'CO2 Levels')
13.  
     
14.  
    plt.legend()
15.  
    plt.show()

我們再次通過計算MSE量化我們預測的預測性能：

1.  
   # Extract the predicted and true values of our time-series
2.  
   y_forecasted = pred_dynamic.predicted_mean
3.  
   y_truth = y[ '1998-01-01':]
4.  
    
5.  
   # Compute the mean square error
6.  
   mse = ((y_forecasted - y_truth) ** 2).mean()
7.  
   print( 'The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

OutputThe Mean Squared Error of our forecasts is 1.01

從動態預測獲得的預測值產生1.01的MSE。 這比前進一步略高，這是意料之中的，因為我們依賴於時間序列中較少的歷史數據。

前進一步和動態預測都證實了這個時間序列模型是有效的。 然而，關於時間序列預測的大部分興趣是能夠及時預測未來價值觀。

第7步 - 生成和可視化預測

在本教程的最后一步，我們將介紹如何利用季節性ARIMA時間序列模型來預測未來的價值。 我們的時間序列對象的

get_forecast()屬性可以計算預先指定數量的步驟的預測值。

1.  
   # Get forecast 500 steps ahead in future
2.  
   pred_uc = results.get_forecast(steps= 500)
3.  
    
4.  
   # Get confidence intervals of forecasts
5.  
   pred_ci = pred_uc.conf_int()

我們可以使用此代碼的輸出繪制其未來值的時間序列和預測。

1.  
   ax = y.plot( label='observed', figsize=(20, 15))
2.  
   pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
3.  
   ax.fill_between(pred_ci. index,
4.  
   pred_ci.iloc[:, 0],
5.  
   pred_ci.iloc[:, 1], color='k', alpha=.25)
6.  
   ax.set_xlabel( 'Date')
7.  
   ax.set_ylabel( 'CO2 Levels')
8.  
    
9.  
   plt.legend()
10.  
    plt.show()

現在可以使用我們生成的預測和相關的置信區間來進一步了解時間序列並預見預期。 我們的預測顯示，時間序列預計將繼續穩步增長。

隨着我們對未來的進一步預測，我們對自己的價值觀信心愈來愈自然。 這反映在我們的模型產生的置信區間，隨着我們進一步走向未來，這個模型越來越大。

結論

在本教程中，我們描述了如何在Python中實現季節性ARIMA模型。 我們廣泛使用了pandas和statsmodels圖書館，並展示了如何運行模型診斷，以及如何產生二氧化碳時間序列的預測。

這里還有一些其他可以嘗試的事情：

• 更改動態預測的開始日期，以了解其如何影響預測的整體質量。
• 嘗試更多的參數組合，看看是否可以提高模型的適合度。
• 選擇不同的指標以選擇最佳模型。 例如，我們使用AIC測量來找到最佳模型，但是您可以尋求優化采樣均方誤差

補充：其中一些好用的方法，代碼如下：

 

1.  
   filename_ts = 'data/series1.csv'
2.  
   ts_df = pd.read_csv(filename_ts, index_col= 0, parse_dates=[0])
3.  
    
4.  
   n_sample = ts_df.shape[ 0]
5.  
   print(ts_df.shape)
6.  
   print(ts_df.head())

1.  
   # Create a training sample and testing sample before analyzing the series
2.  
    
3.  
   n_train= int(0.95*n_sample)+1
4.  
   n_forecast=n_sample-n_train
5.  
   #ts_df
6.  
   ts_train = ts_df.iloc[:n_train][ 'value']
7.  
   ts_test = ts_df.iloc[n_train:][ 'value']
8.  
   print(ts_train.shape)
9.  
   print(ts_test.shape)
10.  
    print("Training Series:", "\n", ts_train.tail(), "\n")
11.  
    print("Testing Series:", "\n", ts_test.head())

 

1.  
   def tsplot(y, lags=None, title='', figsize=(14, 8)):
2.  
    
3.  
   fig = plt.figure(figsize=figsize)
4.  
   layout = ( 2, 2)
5.  
   ts_ax = plt.subplot2grid(layout, ( 0, 0))
6.  
   hist_ax = plt.subplot2grid(layout, ( 0, 1))
7.  
   acf_ax = plt.subplot2grid(layout, ( 1, 0))
8.  
   pacf_ax = plt.subplot2grid(layout, ( 1, 1))
9.  
    
10.  
    y.plot(ax=ts_ax) # 折線圖
11.  
    ts_ax.set_title(title)
12.  
    y.plot(ax=hist_ax, kind= 'hist', bins=25) #直方圖
13.  
    hist_ax.set_title( 'Histogram')
14.  
    smt.graphics.plot_acf(y, lags=lags, ax=acf_ax) # ACF自相關系數
15.  
    smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax) # 偏自相關系數
16.  
    [ax.set_xlim( 0) for ax in [acf_ax, pacf_ax]]
17.  
    sns.despine()
18.  
    fig.tight_layout()
19.  
    return ts_ax, acf_ax, pacf_ax
20.  
     
21.  
    tsplot(ts_train, title= 'A Given Training Series', lags=20);

 

                           

1.  
   # 此處運用BIC（貝葉斯信息准則）進行模型參數選擇
2.  
   # 另外還可以利用AIC（赤池信息准則），視具體情況而定
3.  
   import itertools
4.  
    
5.  
   p_min = 0
6.  
   d_min = 0
7.  
   q_min = 0
8.  
   p_max = 4
9.  
   d_max = 0
10.  
    q_max = 4
11.  
     
12.  
    # Initialize a DataFrame to store the results
13.  
    results_bic = pd.DataFrame(index=[ 'AR{}'.format(i) for i in range(p_min,p_max+1)],
14.  
    columns=[ 'MA{}'.format(i) for i in range(q_min,q_max+1)])
15.  
     
16.  
    for p,d,q in itertools.product(range(p_min,p_max+1),
17.  
    range(d_min,d_max+ 1),
18.  
    range(q_min,q_max+ 1)):
19.  
    if p==0 and d==0 and q==0:
20.  
    results_bic.loc[ 'AR{}'.format(p), 'MA{}'.format(q)] = np.nan
21.  
    continue
22.  
     
23.  
    try:
24.  
    model = sm.tsa.SARIMAX(ts_train, order=(p, d, q),
25.  
    #enforce_stationarity=False,
26.  
    #enforce_invertibility=False,
27.  
    )
28.  
    results = model.fit() ##下面可以顯示具體的參數結果表
29.  
    ## print(model_results.summary())
30.  
    ## print(model_results.summary().tables[1])
31.  
    # http://www.statsmodels.org/stable/tsa.html
32.  
    # print("results.bic",results.bic)
33.  
    # print("results.aic",results.aic)
34.  
     
35.  
    results_bic.loc[ 'AR{}'.format(p), 'MA{}'.format(q)] = results.bic
36.  
    except:
37.  
    continue
38.  
    results_bic = results_bic[results_bic.columns].astype(float)

results_bic如下所示：

為了便於觀察，下面對上表進行可視化：、

1.  
   fig, ax = plt.subplots(figsize=( 10, 8))
2.  
   ax = sns.heatmap(results_bic,
3.  
   mask=results_bic.isnull(),
4.  
   ax=ax,
5.  
   annot=True,
6.  
   fmt= '.2f',
7.  
   );
8.  
   ax.set_title( 'BIC');
9.  
   //annot
10.  
    //annotate的縮寫，annot默認為False，當annot為True時，在heatmap中每個方格寫入數據
11.  
    //annot_kws，當annot為True時，可設置各個參數，包括大小，顏色，加粗，斜體字等
12.  
    # annot_kws={ 'size':9,'weight':'bold', 'color':'blue'}
13.  
    #具體查看：https: //blog.csdn.net/m0_38103546/article/details/79935671

                                          

1.  
   # Alternative model selection method, limited to only searching AR and MA parameters
2.  
    
3.  
   train_results = sm.tsa.arma_order_select_ic(ts_train, ic=[ 'aic', 'bic'], trend='nc', max_ar=4, max_ma=4)
4.  
    
5.  
   print( 'AIC', train_results.aic_min_order)
6.  
   print( 'BIC', train_results.bic_min_order)

注：

SARIMA總代碼如下：

1.  
    
2.  
   """
3.  
   https://www.howtoing.com/a-guide-to-time-series-forecasting-with-arima-in-python-3/
4.  
   http://www.statsmodels.org/stable/datasets/index.html
5.  
   """
6.  
    
7.  
   import warnings
8.  
   import itertools
9.  
   import pandas as pd
10.  
    import numpy as np
11.  
    import statsmodels.api as sm
12.  
    import matplotlib.pyplot as plt
13.  
    plt.style.use( 'fivethirtyeight')
14.  
     
15.  
    '''
16.  
    1-Load Data
17.  
    Most sm.datasets hold convenient representations of the data in the attributes endog and exog.
18.  
    If the dataset does not have a clear interpretation of what should be an endog and exog,
19.  
    then you can always access the data or raw_data attributes.
20.  
    https://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html
21.  
    http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.resample.html
22.  
    # resample('MS') groups the data in buckets by start of the month,
23.  
    # after that we got only one value for each month that is the mean of the month
24.  
    # fillna() fills NA/NaN values with specified method
25.  
    # 'bfill' method use Next valid observation to fill gap
26.  
    # If the value for June is NaN while that for July is not, we adopt the same value
27.  
    # as in July for that in June
28.  
    '''
29.  
     
30.  
    data = sm.datasets.co2.load_pandas()
31.  
    y = data.data # DataFrame with attributes y.columns & y.index (DatetimeIndex)
32.  
    print(y)
33.  
    names = data.names # tuple
34.  
    raw = data.raw_data # float64 np.recarray
35.  
     
36.  
    y = y[ 'co2'].resample('MS').mean()
37.  
    print(y)
38.  
     
39.  
    y = y.fillna(y.bfill()) # y = y.fillna(method='bfill')
40.  
    print(y)
41.  
     
42.  
    y.plot(figsize=( 15,6))
43.  
    plt.show()
44.  
     
45.  
    '''
46.  
    2-ARIMA Parameter Seletion
47.  
    ARIMA(p,d,q)(P,D,Q)s
48.  
    non-seasonal parameters: p,d,q
49.  
    seasonal parameters: P,D,Q
50.  
    s: period of time series, s=12 for annual period
51.  
    Grid Search to find the best combination of parameters
52.  
    Use AIC value to judge models, choose the parameter combination whose
53.  
    AIC value is the smallest
54.  
    https://docs.python.org/3/library/itertools.html
55.  
    http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html
56.  
    '''
57.  
     
58.  
    # Define the p, d and q parameters to take any value between 0 and 2
59.  
    p = d = q = range( 0, 2)
60.  
     
61.  
    # Generate all different combinations of p, q and q triplets
62.  
    pdq = list(itertools.product(p, d, q))
63.  
     
64.  
    # Generate all different combinations of seasonal p, q and q triplets
65.  
    seasonal_pdq = [(x[ 0], x[1], x[2], 12) for x in pdq]
66.  
     
67.  
    print( 'Examples of parameter combinations for Seasonal ARIMA...')
68.  
    print( 'SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
69.  
    print( 'SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
70.  
    print( 'SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
71.  
    print( 'SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
72.  
     
73.  
    warnings.filterwarnings( "ignore") # specify to ignore warning messages
74.  
     
75.  
    for param in pdq:
76.  
    for param_seasonal in seasonal_pdq:
77.  
    try:
78.  
    mod = sm.tsa.statespace.SARIMAX(y,
79.  
    order=param,
80.  
    seasonal_order=param_seasonal,
81.  
    enforce_stationarity= False,
82.  
    enforce_invertibility= False)
83.  
     
84.  
    results = mod.fit()
85.  
     
86.  
    print( 'ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
87.  
    except:
88.  
    continue
89.  
     
90.  
     
91.  
    '''
92.  
    3-Optimal Model Analysis
93.  
    Use the best parameter combination to construct ARIMA model
94.  
    Here we use ARIMA(1,1,1)(1,1,1)12
95.  
    the output coef represents the importance of each feature
96.  
    mod.fit() returnType: MLEResults
97.  
    http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEResults.html#statsmodels.tsa.statespace.mlemodel.MLEResults
98.  
    Use plot_diagnostics() to check if parameters are against the model hypothesis
99.  
    model residuals must not be correlated
100.  
     '''
101.  
      
102.  
     mod = sm.tsa.statespace.SARIMAX(y,
103.  
     order=( 1, 1, 1),
104.  
     seasonal_order=( 1, 1, 1, 12),
105.  
     enforce_stationarity= False,
106.  
     enforce_invertibility= False)
107.  
      
108.  
     results = mod.fit()
109.  
     print(results.summary().tables[ 1])
110.  
      
111.  
     results.plot_diagnostics(figsize=( 15, 12))
112.  
     plt.show()
113.  
      
114.  
     '''
115.  
     4-Make Predictions
116.  
     get_prediction(..., dynamic=False)
117.  
     Prediction of each point will use all historic observations prior to it
118.  
     http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEResults.get_prediction.html#statsmodels.regression.recursive_ls.MLEResults.get_prediction
119.  
     http://pandas.pydata.org/pandas-docs/stable/generated/pandas.Series.plot.html
120.  
     https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.fill_between.html
121.  
     '''
122.  
      
123.  
     pred = results.get_prediction(start=pd.to_datetime( '1998-01-01'), dynamic=False)
124.  
     pred_ci = pred.conf_int() # return the confidence interval of fitted parameters
125.  
      
126.  
     # plot real values and predicted values
127.  
     # pred.predicted_mean is a pandas series
128.  
     ax = y[ '1990':].plot(label='observed') # ax is a matplotlib.axes.Axes
129.  
     pred.predicted_mean.plot(ax=ax, label= 'One-step ahead Forecast', alpha=.7)
130.  
      
131.  
     # fill_between(x,y,z) fills the area between two horizontal curves defined by (x,y)
132.  
     # and (x,z). And alpha refers to the alpha transparencies
133.  
     ax.fill_between(pred_ci.index,
134.  
     pred_ci.iloc[:, 0],
135.  
     pred_ci.iloc[:, 1], color='k', alpha=.2)
136.  
      
137.  
     ax.set_xlabel( 'Date')
138.  
     ax.set_ylabel( 'CO2 Levels')
139.  
     plt.legend()
140.  
      
141.  
     plt.show()
142.  
      
143.  
     # Evaluation of model
144.  
     y_forecasted = pred.predicted_mean
145.  
     y_truth = y[ '1998-01-01':]
146.  
      
147.  
     # Compute the mean square error
148.  
     mse = ((y_forecasted - y_truth) ** 2).mean()
149.  
     print( 'The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
150.  
      
151.  
     '''
152.  
     5-Dynamic Prediction
153.  
     get_prediction(..., dynamic=True)
154.  
     Prediction of each point will use all historic observations prior to 'start' and
155.  
     all predicted values prior to the point to predict
156.  
     '''
157.  
      
158.  
     pred_dynamic = results.get_prediction(start=pd.to_datetime( '1998-01-01'), dynamic=True, full_results=True)
159.  
     pred_dynamic_ci = pred_dynamic.conf_int()
160.  
      
161.  
     ax = y[ '1990':].plot(label='observed', figsize=(20, 15))
162.  
     pred_dynamic.predicted_mean.plot(label= 'Dynamic Forecast', ax=ax)
163.  
      
164.  
     ax.fill_between(pred_dynamic_ci.index,
165.  
     pred_dynamic_ci.iloc[:, 0],
166.  
     pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
167.  
      
168.  
     ax.fill_betweenx(ax.get_ylim(), pd.to_datetime( '1998-01-01'), y.index[-1],
169.  
     alpha= .1, zorder=-1)
170.  
      
171.  
     ax.set_xlabel( 'Date')
172.  
     ax.set_ylabel( 'CO2 Levels')
173.  
      
174.  
     plt.legend()
175.  
     plt.show()
176.  
      
177.  
     # Extract the predicted and true values of our time-series
178.  
     y_forecasted = pred_dynamic.predicted_mean
179.  
     y_truth = y[ '1998-01-01':]
180.  
      
181.  
     # Compute the mean square error
182.  
     mse = ((y_forecasted - y_truth) ** 2).mean()
183.  
     print( 'The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
184.  
      
185.  
     '''
186.  
     6-Visualize Prediction
187.  
     In-sample forecast: forecasting for an observation that was part of the data sample;
188.  
     Out-of-sample forecast: forecasting for an observation that was not part of the data sample.
189.  
     '''
190.  
      
191.  
     # Get forecast 500 steps ahead in future
192.  
     # 'steps': If an integer, the number of steps to forecast from the end of the sample.
193.  
     pred_uc = results.get_forecast(steps= 500) # retun out-of-sample forecast
194.  
      
195.  
     # Get confidence intervals of forecasts
196.  
     pred_ci = pred_uc.conf_int()
197.  
      
198.  
     ax = y.plot(label= 'observed', figsize=(20, 15))
199.  
     pred_uc.predicted_mean.plot(ax=ax, label= 'Forecast')
200.  
     ax.fill_between(pred_ci.index,
201.  
     pred_ci.iloc[:, 0],
202.  
     pred_ci.iloc[:, 1], color='k', alpha=.25)
203.  
     ax.set_xlabel( 'Date')
204.  
     ax.set_ylabel( 'CO2 Levels')
205.  
      
206.  
     plt.legend()
207.  
     plt.show()