python時間序列分析（ARIMA模型）
原文地址：https://blog.csdn.net/u011596455/article/details/78650458
轉載請註明出處。
什麼是時間序列
時間序列簡單的說就是各時間點上形成的數值序列，時間序列分析就是通過觀察歷史資料預測未來的值。在這裡需要強調一點的是，時間序列分析並不是關於時間的迴歸，它主要是研究自身的變化規律的（這裡不考慮含外生變數的時間序列）。
為什麼用python
用兩個字總結“情懷”，愛屋及烏，個人比較喜歡python，就用python擼了。能做時間序列的軟體很多，SAS、R、SPSS、Eviews甚至matlab等等，實際工作中應用得比較多的應該還是SAS和R，前者推薦王燕寫的《應用時間序列分析》，後者推薦“
環境配置
python推薦直接裝Anaconda，它集成了許多科學計算包，有一些包自己手動去裝還是挺費勁的。statsmodels需要自己去安裝，這裡我推薦使用0.6的穩定版，0.7及其以上的版本能在github上找到，該版本在安裝時會用C編譯好，所以修改底層的一些程式碼將不會起作用。
時間序列分析
1.基本模型
自迴歸移動平均模型(ARMA(p，q))是時間序列中最為重要的模型之一，它主要由兩部分組成： AR代表p階自迴歸過程，MA代表q階移動平均過程，其公式如下：
依據模型的形式、特性及自相關和偏自相關函式的特徵，總結如下：
在時間序列中，ARIMA模型是在ARMA模型的基礎上多了差分的操作。
2.pandas時間序列操作
大熊貓真的很可愛，這裡簡單介紹一下它在時間序列上的可愛之處。和許多時間序列分析一樣，本文同樣使用航空乘客資料（AirPassengers.csv）作為樣例。
資料讀取：
# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pylab as plt
# 讀取資料，pd.read_csv預設生成DataFrame物件，需將其轉換成Series物件
df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date')
df.index = pd.to_datetime(df.index) # 將字串索引轉換成時間索引
ts = df['x'] # 生成pd.Series物件
# 檢視資料格式
ts.head()
ts.head().index
檢視某日的值既可以使用字串作為索引，又可以直接使用時間物件作為索引
ts['1949-01-01']
ts[datetime(1949,1,1)]
兩者的返回值都是第一個序列值：112
如果要檢視某一年的資料，pandas也能非常方便的實現
ts['1949']
切片操作：
ts['1949-1' : '1949-6']
注意時間索引的切片操作起點和尾部都是包含的，這點與數值索引有所不同
pandas還有很多方便的時間序列函式，在後面的實際應用中在進行說明。
3. 平穩性檢驗
我們知道序列平穩性是進行時間序列分析的前提條件，很多人都會有疑問，為什麼要滿足平穩性的要求呢？在大數定理和中心定理中要求樣本同分布（這裡同分布等價於時間序列中的平穩性），而我們的建模過程中有很多都是建立在大數定理和中心極限定理的前提條件下的，如果它不滿足，得到的許多結論都是不可靠的。以虛假迴歸為例，當響應變數和輸入變數都平穩時，我們用t統計量檢驗標準化係數的顯著性。而當響應變數和輸入變數不平穩時，其標準化係數不在滿足t分佈，這時再用t檢驗來進行顯著性分析，導致拒絕原假設的概率增加，即容易犯第一類錯誤，從而得出錯誤的結論。
平穩時間序列有兩種定義：嚴平穩和寬平穩
嚴平穩顧名思義，是一種條件非常苛刻的平穩性，它要求序列隨著時間的推移，其統計性質保持不變。對於任意的τ，其聯合概率密度函式滿足：
嚴平穩的條件只是理論上的存在，現實中用得比較多的是寬平穩的條件。
寬平穩也叫弱平穩或者二階平穩（均值和方差平穩），它應滿足：
- 常數均值
- 常數方差
- 常數自協方差
平穩性檢驗：觀察法和單位根檢驗法
基於此，我寫了一個名為test_stationarity的統計性檢驗模組，以便將某些統計檢驗結果更加直觀的展現出來。
# -*- coding:utf-8 -*-
from statsmodels.tsa.stattools import adfuller
import pandas as pd import matplotlib.pyplot as plt import numpy as np from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # 移動平均圖 def draw_trend(timeSeries, size): f = plt.figure(facecolor='white') # 對size個數據進行移動平均 rol_mean = timeSeries.rolling(window=size).mean() # 對size個數據進行加權移動平均 rol_weighted_mean = pd.ewma(timeSeries, span=size) timeSeries.plot(color='blue', label='Original') rolmean.plot(color='red', label='Rolling Mean') rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean') plt.legend(loc='best') plt.title('Rolling Mean') plt.show() def draw_ts(timeSeries): f = plt.figure(facecolor='white') timeSeries.plot(color='blue') plt.show() ''' Unit Root Test The null hypothesis of the Augmented Dickey-Fuller is that there is a unit root, with the alternative that there is no unit root. That is to say the bigger the p-value the more reason we assert that there is a unit root ''' def testStationarity(ts): dftest = adfuller(ts) # 對上述函式求得的值進行語義描述 dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used']) for key,value in dftest[4].items(): dfoutput['Critical Value (%s)'%key] = value return dfoutput # 自相關和偏相關圖，預設階數為31階 def draw_acf_pacf(ts, lags=31): f = plt.figure(facecolor='white') ax1 = f.add_subplot(211) plot_acf(ts, lags=31, ax=ax1) ax2 = f.add_subplot(212) plot_pacf(ts, lags=31, ax=ax2) plt.show()
觀察法，通俗的說就是通過觀察序列的趨勢圖與相關圖是否隨著時間的變化呈現出某種規律。所謂的規律就是時間序列經常提到的週期性因素，現實中遇到得比較多的是線性週期成分，這類週期成分可以採用差分或者移動平均來解決，而對於非線性週期成分的處理相對比較複雜，需要採用某些分解的方法。下圖為航空資料的線性圖，可以明顯的看出它具有年週期成分和長期趨勢成分。平穩序列的自相關係數會快速衰減，下面的自相關圖並不能體現出該特徵，所以我們有理由相信該序列是不平穩的。
單位根檢驗：ADF是一種常用的單位根檢驗方法，他的原假設為序列具有單位根，即非平穩，對於一個平穩的時序資料，就需要在給定的置信水平上顯著，拒絕原假設。ADF只是單位根檢驗的方法之一，如果想採用其他檢驗方法，可以安裝第三方包arch，裡面提供了更加全面的單位根檢驗方法，個人還是比較鍾情ADF檢驗。以下為檢驗結果，其p值大於0.99，說明並不能拒絕原假設。
3. 平穩性處理
由前面的分析可知，該序列是不平穩的，然而平穩性是時間序列分析的前提條件，故我們需要對不平穩的序列進行處理將其轉換成平穩的序列。
a. 對數變換
對數變換主要是為了減小資料的振動幅度，使其線性規律更加明顯（我是這麼理解的時間序列模型大部分都是線性的，為了儘量降低非線性的因素，需要對其進行預處理，也許我理解的不對）。對數變換相當於增加了一個懲罰機制，資料越大其懲罰越大，資料越小懲罰越小。這裡強調一下，變換的序列需要滿足大於0，小於0的資料不存在對數變換。
ts_log = np.log(ts)
test_stationarity.draw_ts(ts_log)
b. 平滑法
根據平滑技術的不同，平滑法具體分為移動平均法和指數平均法。
移動平均即利用一定時間間隔內的平均值作為某一期的估計值，而指數平均則是用變權的方法來計算均值
test_stationarity.draw_trend(ts_log, 12)
從上圖可以發現視窗為12的移動平均能較好的剔除年週期性因素，而指數平均法是對週期內的資料進行了加權，能在一定程度上減小年週期因素，但並不能完全剔除，如要完全剔除可以進一步進行差分操作。
c. 差分
時間序列最常用來剔除週期性因素的方法當屬差分了，它主要是對等週期間隔的資料進行線性求減。前面我們說過，ARIMA模型相對ARMA模型，僅多了差分操作，ARIMA模型幾乎是所有時間序列軟體都支援的，差分的實現與還原都非常方便。而statsmodel中，對差分的支援不是很好，它不支援高階和多階差分，為什麼不支援，這裡引用作者的說法：
作者大概的意思是說預測方法中並沒有解決高於2階的差分，有沒有感覺很牽強，不過沒關係，我們有pandas。我們可以先用pandas將序列差分好，然後在對差分好的序列進行ARIMA擬合，只不過這樣後面會多了一步人工還原的工作。
diff_12 = ts_log.diff(12)
diff_12.dropna(inplace=True)
diff_12_1 = diff_12.diff(1)
diff_12_1.dropna(inplace=True)
test_stationarity.testStationarity(diff_12_1)
從上面的統計檢驗結果可以看出，經過12階差分和1階差分後，該序列滿足平穩性的要求了。
d. 分解
所謂分解就是將時序資料分離成不同的成分。statsmodels使用的X-11分解過程，它主要將時序資料分離成長期趨勢、季節趨勢和隨機成分。與其它統計軟體一樣，statsmodels也支援兩類分解模型，加法模型和乘法模型，這裡我只實現加法，乘法只需將model的引數設定為"multiplicative"即可。
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log, model="additive") trend = decomposition.trend seasonal = decomposition.seasonal residual = decomposition.resid
得到不同的分解成分後，就可以使用時間序列模型對各個成分進行擬合，當然也可以選擇其他預測方法。我曾經用過小波對時序資料進行過分解，然後分別採用時間序列擬合，效果還不錯。由於我對小波的理解不是很好，只能簡單的呼叫介面，如果有誰對小波、傅立葉、卡爾曼理解得比較透，可以將時序資料進行更加準確的分解，由於分解後的時序資料避免了他們在建模時的交叉影響，所以我相信它將有助於預測準確性的提高。
4. 模型識別
在前面的分析可知，該序列具有明顯的年週期與長期成分。對於年週期成分我們使用視窗為12的移動平進行處理，對於長期趨勢成分我們採用1階差分來進行處理。
rol_mean = ts_log.rolling(window=12).mean()
rol_mean.dropna(inplace=True)
ts_diff_1 = rol_mean.diff(1)
ts_diff_1.dropna(inplace=True)
test_stationarity.testStationarity(ts_diff_1)
觀察其統計量發現該序列在置信水平為95%的區間下並不顯著，我們對其進行再次一階差分。再次差分後的序列其自相關具有快速衰減的特點，t統計量在99%的置信水平下是顯著的，這裡我不再做詳細說明。
ts_diff_2 = ts_diff_1.diff(1)
ts_diff_2.dropna(inplace=True)
資料平穩後，需要對模型定階，即確定p、q的階數。觀察上圖，發現自相關和偏相係數都存在拖尾的特點，並且他們都具有明顯的一階相關性，所以我們設定p=1, q=1。下面就可以使用ARMA模型進行資料擬合了。這裡我不使用ARIMA(ts_diff_1, order=(1, 1, 1))進行擬合，是因為含有差分操作時，預測結果還原老出問題，至今還沒弄明白。
from statsmodels.tsa.arima_model import ARMA
model = ARMA(ts_diff_2, order=(1, 1))
result_arma = model.fit( disp=-1, method='css')
5. 樣本擬合
模型擬合完後，我們就可以對其進行預測了。由於ARMA擬合的是經過相關預處理後的資料，故其預測值需要通過相關逆變換進行還原。
predict_ts = result_arma.predict()
# 一階差分還原
diff_shift_ts = ts_diff_1.shift(1)
diff_recover_1 = predict_ts.add(diff_shift_ts)
# 再次一階差分還原
rol_shift_ts = rol_mean.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts) # 移動平均還原 rol_sum = ts_log.rolling(window=11).sum() rol_recover = diff_recover*12 - rol_sum.shift(1) # 對數還原 log_recover = np.exp(rol_recover) log_recover.dropna(inplace=True)
我們使用均方根誤差（RMSE）來評估模型樣本內擬合的好壞。利用該準則進行判別時，需要剔除“非預測”資料的影響。
ts = ts[log_recover.index] # 過濾沒有預測的記錄
plt.figure(facecolor='white')
log_recover.plot(color='blue', label='Predict') ts.plot(color='red', label='Original') plt.legend(loc='best') plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size)) plt.show()
觀察上圖的擬合效果，均方根誤差為11.8828，感覺還過得去。
6. 完善ARIMA模型
前面提到statsmodels裡面的ARIMA模組不支援高階差分，我們的做法是將差分分離出來，但是這樣會多了一步人工還原的操作。基於上述問題，我將差分過程進行了封裝，使序列能按照指定的差分列表依次進行差分，並相應的構造了一個還原的方法，實現差分序列的自動還原。
# 差分操作
def diff_ts(ts, d):
global shift_ts_list # 動態預測第二日的值時所需要的差分序列 global last_data_shift_list shift_ts_list = [] last_data_shift_list = [] tmp_ts = ts for i in d: last_data_shift_list.append(tmp_ts[-i]) print last_data_shift_list shift_ts = tmp_ts.shift(i) shift_ts_list.append(shift_ts) tmp_ts = tmp_ts - shift_ts tmp_ts.dropna(inplace=True) return tmp_ts # 還原操作 def predict_diff_recover(predict_value, d): if isinstance(predict_value, float): tmp_data = predict_value for i in range(len(d)): tmp_data = tmp_data + last_data_shift_list[-i-1] elif isinstance(predict_value, np.ndarray): tmp_data = predict_value[0] for i in range(len(d)): tmp_data = tmp_data + last_data_shift_list[-i-1] else: tmp_data = predict_value for i in range(len(d)): try: tmp_data = tmp_data.add(shift_ts_list[-i-1]) except: raise ValueError('What you input is not pd.Series type!') tmp_data.dropna(inplace=True) return tmp_data
現在我們直接使用差分的方法進行資料處理，並以同樣的過程進行資料預測與還原。
diffed_ts = diff_ts(ts_log, d=[12, 1])
model = arima_model(diffed_ts)
model.certain_model(1, 1)
predict_ts = model.properModel.predict()
diff_recover_ts = predict_diff_recover(predict_ts, d=[12, 1])
log_recover = np.exp(diff_recover_ts)
是不是發現這裡的預測結果和上一篇的使用12階移動平均的預測結果一模一樣。這是因為12階移動平均加上一階差分與直接12階差分是等價的關係，後者是前者數值的12倍，這個應該不難推導。
對於個數不多的時序資料，我們可以通過觀察自相關圖和偏相關圖來進行模型識別，倘若我們要分析的時序資料量較多，例如要預測每隻股票的走勢，我們就不可能逐個去調參了。這時我們可以依據BIC準則識別模型的p, q值，通常認為BIC值越小的模型相對更優。這裡我簡單介紹一下BIC準則，它綜合考慮了殘差大小和自變數的個數，殘差越小BIC值越小，自變數個數越多BIC值越大。個人覺得BIC準則就是對模型過擬合設定了一個標準（過擬合這東西應該以辯證的眼光看待）。
def proper_model(data_ts, maxLag):
init_bic = sys.maxint
init_p = 0
init_q = 0
init_properModel = None for p in np.arange(maxLag): for q in np.arange(maxLag): model = ARMA(data_ts, order=(p, q)) try: results_ARMA = model.fit(disp=-1, method='css') except: continue bic = results_ARMA.bic if bic < init_bic: init_p = p init_q = q init_properModel = results_ARMA init_bic = bic return init_bic, init_p, init_q, init_properModel
相對最優引數識別結果：BIC: -1090.44209358 p: 0 q: 1 ， RMSE:11.8817198331。我們發現模型自動識別的引數要比我手動選取的引數更優。
7.滾動預測
所謂滾動預測是指通過新增最新的資料預測第二天的值。對於一個穩定的預測模型，不需要每天都去擬合，我們可以給他設定一個閥值，例如每週擬合一次，該期間只需通過新增最新的資料實現滾動預測即可。基於此我編寫了一個名為arima_model的類，主要包含模型自動識別方法，滾動預測的功能，詳細程式碼可以檢視附錄。資料的動態新增：
from dateutil.relativedelta import relativedelta
def _add_new_data(ts, dat, type='day'):
if type == 'day': new_index = ts.index[-1] + relativedelta(days=1) elif type == 'month': new_index = ts.index[-1] + relativedelta(months=1) ts[new_index] = dat def add_today_data(model, ts, data, d, type='day'): _add_new_data(ts, data, type) # 為原始序列新增資料 # 為滯後序列新增新值 d_ts = diff_ts(ts, d) model.add_today_data(d_ts[-1], type) def forecast_next_day_data(model, type='day'): if model == None: raise ValueError('No model fit before') fc = model.forecast_next_day_value(type) return predict_diff_recover(fc, [12, 1])
現在我們就可以使用滾動預測的方法向外預測了，取1957年之前的資料作為訓練資料，其後的資料作為測試，並設定模型每第七天就會重新擬合一次。這裡的diffed_ts物件會隨著add_today_data方法自動新增資料，這是由於它與add_today_data方法中的d_ts指向的同一物件，該物件會動態的新增資料。
ts_train = ts_log[:'1956-12']
ts_test = ts_log['1957-1':] diffed_ts = diff_ts(ts_train, [12, 1]) forecast_list = [] for i, dta in enumerate(ts_test): if i%7 == 0: model = arima_model(diffed_ts) model.certain_model(1, 1) forecast_data = forecast_next_day_data(model, type='month') forecast_list.append(forecast_data) add_today_data(model, ts_train, dta, [12, 1], type='month') predict_ts = pd.Series(data=forecast_list, index=ts['1957-1':].index) log_recover = np.exp(predict_ts) original_ts = ts['1957-1':]
動態預測的均方根誤差為：14.6479，與前面樣本內擬合的均方根誤差相差不大，說明模型並沒有過擬合，並且整體預測效果都較好。
8. 模型序列化
在進行動態預測時，我們不希望將整個模型一直在記憶體中執行，而是希望有新的資料到來時才啟動該模型。這時我們就應該把整個模型從記憶體匯出到硬碟中，而序列化正好能滿足該要求。序列化最常用的就是使用json模組了，但是它是時間物件支援得不是很好，有人對json模組進行了拓展以使得支援時間物件，這裡我們不採用該方法，我們使用pickle模組，它和json的介面基本相同，有興趣的可以去檢視一下。我在實際應用中是採用的面向物件的程式設計，它的序列化主要是將類的屬性序列化即可，而在面向過程的程式設計中，模型序列化需要將需要序列化的物件公有化，這樣會使得對前面函式的引數改動較大，我不在詳細闡述該過程。
總結
與其它統計語言相比，python在統計分析這塊還顯得不那麼“專業”。隨著numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推動，我相信也祝願python在資料分析這塊越來越好。與SAS和R相比，python的時間序列模組還不是很成熟，我這裡僅起到拋磚引玉的作用，希望各位能人志士能貢獻自己的力量，使其更加完善。實際應用中我全是面向過程來編寫的，為了闡述方便，我用面向過程重新羅列了一遍，實在感覺很不方便。原本打算分三篇來寫的，還有一部分實際應用的部分，不打算再寫了，還請大家原諒。實際應用主要是具體問題具體分析，這當中第一步就是要查詢問題，這步花的時間往往會比較多，然後再是解決問題。以我前面專案遇到的問題為例，當時遇到了以下幾個典型的問題：1.週期長度不恆定的週期成分，例如每月的1號具有周期性，但每月1號與1號之間的時間間隔是不相等的；2.含有缺失值以及含有記錄為0的情況無法進行對數變換；3.節假日的影響等等。
1 # -*-coding:utf-8-*- 2 import pandas as pd 3 import numpy as np 4 from statsmodels.tsa.arima_model import ARMA 5 import sys 6 from dateutil.relativedelta import relativedelta 7 from copy import deepcopy 8 import matplotlib.pyplot as plt 9 10 class arima_model: 11 12 def __init__(self, ts, maxLag=9): 13 self.data_ts = ts 14 self.resid_ts = None 15 self.predict_ts = None 16 self.maxLag = maxLag 17 self.p = maxLag 18 self.q = maxLag 19 self.properModel = None 20 self.bic = sys.maxint 21 22 # 計算最優ARIMA模型，將相關結果賦給相應屬性 23 def get_proper_model(self): 24 self._proper_model() 25 self.predict_ts = deepcopy(self.properModel.predict()) 26 self.resid_ts = deepcopy(self.properModel.resid) 27 28 # 對於給定範圍內的p,q計算擬合得最好的arima模型，這裡是對差分好的資料進行擬合，故差分恆為0 29 def _proper_model(self): 30 for p in np.arange(self.maxLag): 31 for q in np.arange(self.maxLag): 32 # print p,q,self.bic 33 model = ARMA(self.data_ts, order=(p, q)) 34 try: 35 results_ARMA = model.fit(disp=-1, method='css') 36 except: 37 continue 38 bic = results_ARMA.bic 39 # print 'bic:',bic,'self.bic:',self.bic 40 if bic < self.bic: 41 self.p = p 42 self.q = q 43 self.properModel = results_ARMA 44 self.bic = bic 45 self.resid_ts = deepcopy(self.properModel.resid) 46 self.predict_ts = self.properModel.predict() 47 48 # 引數確定模型 49 def certain_model(self, p, q): 50 model = ARMA(self.data_ts, order=(p, q)) 51 try: 52 self.properModel = model.fit( disp=-1, method='css') 53 self.p = p 54 self.q = q 55 self.bic = self.properModel.bic 56 self.predict_ts = self.properModel.predict() 57 self.resid_ts = deepcopy(self.properModel.resid) 58 except: 59 print 'You can not fit the model with this parameter p,q, ' \ 60 'please use the get_proper_model method to get the best model' 61 62 # 預測第二日的值 63 def forecast_next_day_value(self, type='day'): 64 # 我修改了statsmodels包中arima_model的原始碼，添加了constant屬性，需要先執行forecast方法，為constant賦值 65 self.properModel.forecast() 66 if self.data_ts.index[-1] != self.resid_ts.index[-1]: 67 raise ValueError('''The index is different in data_ts and resid_ts, please add new data to data_ts. 68 If you just want to forecast the next day data without add the real next day data to data_ts, 69 please run the predict method which arima_model included itself''') 70 if not self.properModel: 71 raise ValueError('The arima model have not computed, please run the proper_model method before') 72 para = self.properModel.params 73 74 # print self.properModel.params 75 if self.p == 0: # It will get all the value series with setting self.data_ts[-self.p:] when p is zero 76 ma_value = self.resid_ts[-self.q:] 77 values = ma_value.reindex(index=ma_value.index[::-1]) 78 elif self.q == 0: 79 ar_value = self.data_ts[-self.p:] 80 values = ar_value.reindex(index=ar_value.index[::-1]) 81 else: 82 ar_value = self.data_ts[-self.p:] 83 ar_value = ar_value.reindex(index=ar_value.index[::-1]) 84 ma_value = self.resid_ts[-self.q:] 85 ma_value = ma_value.reindex(index=ma_value.index[::-1]) 86 values = ar_value.append(ma_value) 87 88 predict_value = np.dot(para[1:], values) + self.properModel.constant[0] 89 self._add_new_data(self.predict_ts, predict_value, type) 90 return predict_value 91 92 # 動態新增資料函式，針對索引是月份和日分別進行處理 93 def _add_new_data(self, ts, dat, type='day'): 94 if type == 'day': 95 new_index = ts.index[-1] + relativedelta(days=1) 96 elif type == 'month': 97 new_index = ts.index[-1] + relativedelta(months=1) 98 ts[new_index] = dat 99 100 def add_today_data(self, dat, type='day'): 101 self._add_new_data(self.data_ts, dat, type) 102 if self.data_ts.index[-1] != self.predict_ts.index[-1]: 103 raise ValueError('You must use the forecast_next_day_value method forecast the value of today before') 104 self._add_new_data(self.resid_ts, self.data_ts[-1] - self.predict_ts[-1], type) 105 106 if __name__ == '__main__': 107 df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date') 108 df.index = pd.to_datetime(df.index) 109 ts = df['x'] 110 111 # 資料預處理 112 ts_log = np.log(ts) 113 rol_mean = ts_log.rolling(window=12).mean() 114 rol_mean.dropna(inplace=True) 115 ts_diff_1 = rol_mean.diff(1) 116 ts_diff_1.dropna(inplace=True) 117 ts_diff_2 = ts_diff_1.diff(1) 118 ts_diff_2.dropna(inplace=True) 119 120 # 模型擬合 121 model = arima_model(ts_diff_2) 122 # 這裡使用模型引數自動識別 123 model.get_proper_model() 124 print 'bic:', model.bic, 'p:', model.p, 'q:', model.q 125 print model.properModel.forecast()[0] 126 print model.forecast_next_day_value(type='month') 127 128 # 預測結果還原 129 predict_ts = model.properModel.predict() 130 diff_shift_ts = ts_diff_1.shift(1) 131 diff_recover_1 = predict_ts.add(diff_shift_ts) 132 rol_shift_ts = rol_mean.shift(1) 133 diff_recover = diff_recover_1.add(rol_shift_ts) 134 rol_sum = ts_log.rolling(window=11).sum() 135 rol_recover = diff_recover*12 - rol_sum.shift(1) 136 log_recover = np.exp(rol_recover) 137 log_recover.dropna(inplace=True) 138 139 # 預測結果作圖 140 ts = ts[log_recover.index] 141 plt.figure(facecolor='white') 142 log_recover.plot(color='blue', label='Predict') 143 ts.plot(color='red', label='Original') 144 plt.legend(loc='best') 145 plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size)) 146 plt.show()
修改的arima_model程式碼
1 # Note: The information criteria add 1 to the number of parameters 2 # whenever the model has an AR or MA term since, in principle, 3 # the variance could be treated as a free parameter and restricted 4 # This code does not allow this, but it adds consistency with other 5 # packages such as gretl and X12-ARIMA 6 7 from __future__ import absolute_import 8 from statsmodels.compat.python import string_types, range 9 # for 2to3 with extensions 10 11 from datetime import datetime 12 13 import numpy as np 14 from scipy import optimize 15 from scipy.stats import t, norm 16 from scipy.signal import lfilter 17 from numpy import dot, log, zeros, pi 18 from numpy.linalg import inv 19 20 from statsmodels.tools.decorators import (cache_readonly, 21 resettable_cache) 22 import statsmodels.tsa.base.tsa_model as tsbase 23 import statsmodels.base.wrapper as wrap 24 from statsmodels.regression.linear_model import yule_walker, GLS 25 from statsmodels.tsa.tsatools import (lagmat, add_trend, 26 _ar_transparams, _ar_invtransparams, 27 _ma_transparams, _ma_invtransparams, 28 unintegrate, unintegrate_levels) 29 from statsmodels.tsa.vector_ar import util 30 from statsmodels.tsa.ar_model import AR 31 from statsmodels.tsa.arima_process import arma2ma 32 from statsmodels.tools.numdiff import approx_hess_cs, approx_fprime_cs 33 from statsmodels.tsa.base.datetools import _index_date 34 from statsmodels.tsa.kalmanf import KalmanFilter 35 36 _armax_notes = """ 37 38 Notes 39 ----- 40 If exogenous variables are given, then the model that is fit is 41 42 .. math:: 43 44 \\phi(L)(y_t - X_t\\beta) = \\theta(L)\epsilon_t 45 46 where :math:`\\phi` and :math:`\\theta` are polynomials in the lag 47 operator, :math:`L`. This is the regression model with ARMA errors, 48 or ARMAX model. This specification is used, whether or not the model 49 is fit using conditional sum of square or maximum-likelihood, using 50 the `method` argument in 51 :meth:`statsmodels.tsa.arima_model.%(Model)s.fit`. Therefore, for 52 now, `css` and `mle` refer to estimation methods only. This may 53 change for the case of the `css` model in future versions. 54 """ 55 56 _arma_params = """\ 57 endog : array-like 58 The endogenous variable. 59 order : iterable 60 The (p,q) order of the model for the number of AR parameters, 61 differences, and MA parameters to use. 62 exog : array-like, optional 63 An optional arry of exogenous variables. This should *not* include a 64 constant or trend. You can specify this in the `fit` method.""" 65 66 _arma_model = "Autoregressive Moving Average ARMA(p,q) Model" 67 68 _arima_model = "Autoregressive Integrated Moving Average ARIMA(p,d,q) Model" 69 70 _arima_params = """\ 71 endog : array-like 72 The endogenous variable. 73 order : iterable 74 The (p,d,q) order of the model for the number of AR parameters, 75 differences, and MA parameters to use. 76 exog : array-like, optional 77 An optional arry of exogenous variables. This should *not* include a 78 constant or trend. You can specify this in the `fit` method.""" 79 80 _predict_notes = """ 81 Notes 82 ----- 83 Use the results predict method instead. 84 """ 85 86 _results_notes = """ 87 Notes 88 ----- 89 It is recommended to use dates with the time-series models, as the 90 below will probably make clear. However, if ARIMA is used without 91 dates and/or `start` and `end` are given as indices, then these 92 indices are in terms of the *original*, undifferenced series. Ie., 93 given some undifferenced observations:: 94 95 1970Q1, 1 96 1970Q2, 1.5 97 1970Q3, 1.25 98 1970Q4, 2.25 99 1971Q1, 1.2 100 1971Q2, 4.1 101 102 1970Q1 is observation 0 in the original series. However, if we fit an 103 ARIMA(p,1,q) model then we lose this first observation through 104 differencing. Therefore, the first observation we can forecast (if 105 using exact MLE) is index 1. In the differenced series this is index 106 0, but we refer to it as 1 from the original series. 107 """ 108 109 _predict = """ 110 %(Model)s model in-sample and out-of-sample prediction 111 112 Parameters 113 ---------- 114 %(params)s 115 start : int, str, or datetime 116 Zero-indexed observation number at which to start forecasting, ie., 117 the first forecast is start. Can also be a date string to 118 parse or a datetime type. 119 end : int, str, or datetime 120 Zero-indexed observation number at which to end forecasting, ie., 121 the first forecast is start. Can also be a date string to 122 parse or a datetime type. However, if the dates index does not 123 have a fixed frequency, end must be an integer index if you 124 want out of sample prediction. 125 exog : array-like, optional 126 If the model is an ARMAX and out-of-sample forecasting is 127 requested, exog must be given. Note that you'll need to pass 128 `k_ar` additional lags for any exogenous variables. E.g., if you 129 fit an ARMAX(2, q) model and want to predict 5 steps, you need 7 130 observations to do this. 131 dynamic : bool, optional 132 The `dynamic` keyword affects in-sample prediction. If dynamic 133 is False, then the in-sample lagged values are used for 134 prediction. If `dynamic` is True, then in-sample forecasts are 135 used in place of lagged dependent variables. The first forecasted 136 value is `start`. 137 %(extra_params)s 138 139 Returns 140 ------- 141 %(returns)s 142 %(extra_section)s 143 """ 144 145 _predict_returns = """predict : array 146 The predicted values. 147 148 """ 149 150 _arma_predict = _predict % {"Model" : "ARMA", 151 "params" : """ 152 params : array-like 153 The fitted parameters of the model.""", 154 "extra_params" : "", 155 "returns" : _predict_returns, 156 "extra_section" : _predict_notes} 157 158 _arma_results_predict = _predict % {"Model" : "ARMA", "params" : "", 159 "extra_params" : "", 160 "returns" : _predict_returns, 161 "extra_section" : _results_notes} 162 163 _arima_predict = _predict % {"Model" : "ARIMA", 164 "params" : """params : array-like 165 The fitted parameters of the model.""", 166 "extra_params" : """typ : str {'linear', 'levels'} 167 168 - 'linear' : Linear prediction in terms of the differenced 169 endogenous variables. 170 - 'levels' : Predict the levels of the original endogenous 171 variables.\n""", "returns" : _predict_returns, 172 "extra_section" : _predict_notes} 173 174 _arima_results_predict = _predict % {"Model" : "ARIMA", 175 "params" : "", 176 "extra_params" : 177 """typ : str {'linear', 'levels'} 178 179 - 'linear' : Linear prediction in terms of the differenced 180 endogenous variables. 181 - 'levels' : Predict the levels of the original endogenous 182 variables.\n""", 183 "returns" : _predict_returns, 184 "extra_section" : _results_notes} 185 186 _arima_plot_predict_example = """ Examples 187 -------- 188 >>> import statsmodels.api as sm 189 >>> import matplotlib.pyplot as plt 190 >>> import pandas as pd 191 >>> 192 >>> dta = sm.datasets.sunspots.load_pandas().data[['SUNACTIVITY']] 193 >>> dta.index = pd.DatetimeIndex(start='1700', end='2009', freq='A') 194 >>> res = sm.tsa.ARMA(dta, (3, 0)).fit() 195 >>> fig, ax = plt.subplots() 196 >>> ax = dta.ix['1950':].plot(ax=ax) 197 >>> fig = res.plot_predict('1990', '2012', dynamic=True, ax=ax, 198 ... plot_insample=False) 199 >>> plt.show() 200 201 .. plot:: plots/arma_predict_plot.py 202 """ 203 204 _plot_predict = (""" 205 Plot forecasts 206 """ + '\n'.join(_predict.split('\n')[2:])) % { 207 "params" : "", 208 "extra_params" : """alpha : float, optional 209 The confidence intervals for the forecasts are (1 - alpha)% 210 plot_insample : bool, optional 211 Whether to plot the in-sample series. Default is True. 212 ax : matplotlib.Axes, optional 213 Existing axes to plot with.""", 214 "returns" : """fig : matplotlib.Figure 215 The plotted Figure instance""", 216 "extra_section" : ('\n' + _arima_plot_predict_example + 217 '\n' + _results_notes) 218 } 219 220 _arima_plot_predict = (""" 221 Plot forecasts 222 """ + '