數(shù)據(jù)探查
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from matplotlib.pylab import rcParams
#rcParams設(shè)定好畫布的大小
rcParams['figure.figsize'] = 15, 6
path = '/Users/***/workspace/Analytics_Vidhya/Articles/Time_Series_Analysis/AirPassengers.csv'
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
#---其中parse_dates 表明選擇數(shù)據(jù)中的哪個column作為date-time信息粟耻，
#---index_col 告訴pandas以哪個column作為 index
#--- date_parser 使用一個function(本文用lambda表達式代替)垢乙，使一個string轉(zhuǎn)換為一個datetime變量
data = pd.read_csv(path, parse_dates=['Month'], index_col='Month',date_parser=dateparse)
data.head(6)
| Month | #Passengers |
|---|---|
| 1949-01-01 | 112 |
| 1949-02-01 | 118 |
| 1949-03-01 | 132 |
| 1949-04-01 | 129 |
| 1949-05-01 | 121 |
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
# 這里以一年為一個窗口锨咙，每一個時間t的值由它前面12個月（包括自己）的均值代替，標準差同理追逮。
# rolmean = pd.rolling_mean(timeseries, window=12)
#
# rolstd = pd.rolling_std(timeseries, window=12)
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(12).std()
# plot rolling statistics:
fig = plt.figure()
fig.add_subplot()
orig = plt.plot(timeseries, color='blue', label='Original')
mean = plt.plot(rolmean, color='red', label='rolling mean')
std = plt.plot(rolstd, color='black', label='Rolling standard deviation')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
# Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
# dftest的輸出前一項依次為檢測值酪刀，p值，滯后數(shù)钮孵，使用的觀測數(shù)骂倘，各個置信度下的臨界值
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
print(dfoutput)
ts = data['#Passengers']
test_stationarity(ts)
Results of Dickey-Fuller Test:
Test Statistic 0.815369
p-value 0.991880
#Lags Used 13.000000
Number of Observations Used 130.000000
Critical value (1%) -3.481682
Critical value (5%) -2.884042
Critical value (10%) -2.578770
dtype: float64
可以看到，數(shù)據(jù)的rolling均值/標準差具有越來越大的趨勢巴席，是不穩(wěn)定的历涝。
且DF-test可以明確的指出，在任何置信度下漾唉，數(shù)據(jù)都不是穩(wěn)定的荧库。
ts_log = np.log(ts) #由于原數(shù)據(jù)值域范圍比較大，為了縮小值域赵刑，同時保留其他信息分衫，常用的方法是對數(shù)化，取log
MA:Moving Average--移動平均
使用以上數(shù)據(jù)的對數(shù)取移動平均
moving_avg = ts_log.rolling(12).mean()
ts_log_moving_avg_diff = ts_log-moving_avg
ts_log_moving_avg_diff.dropna(inplace = True) #去除na值 inplace=True代表就地替換
test_stationarity(ts_log_moving_avg_diff)
Results of Dickey-Fuller Test:
Test Statistic -3.162908
p-value 0.022235
#Lags Used 13.000000
Number of Observations Used 119.000000
Critical value (1%) -3.486535
Critical value (5%) -2.886151
Critical value (10%) -2.579896
dtype: float64
數(shù)據(jù)經(jīng)過和移動平均作差后 躺平了
可以看到般此，做了處理之后的數(shù)據(jù)基本上沒有了隨時間變化的趨勢蚪战，DFtest的結(jié)果告訴我們在95%的置信度下，數(shù)據(jù)是穩(wěn)定的铐懊。
上面的方法是將所有的時間平等看待邀桑，而在許多情況下，可以認為越近的時刻越重要科乎。所以引入指數(shù)加權(quán)移動平均-- Exponentially-weighted moving average.（pandas中通過ewma()函數(shù)提供了此功能壁畸。）
# halflife的值決定了衰減因子alpha： alpha = 1 - exp(log(0.5) / halflife)
expweighted_avg = ts_log.ewm(halflife=12).mean() # 這里的halflife=12的意義還不大清楚
ts_log_ewma_diff = ts_log - expweighted_avg
test_stationarity(ts_log_ewma_diff)
Results of Dickey-Fuller Test:
Test Statistic -3.601262
p-value 0.005737
#Lags Used 13.000000
Number of Observations Used 130.000000
Critical value (1%) -3.481682
Critical value (5%) -2.884042
Critical value (10%) -2.578770
dtype: float64
可以看到相比普通的Moving Average，新的數(shù)據(jù)平均標準差更小了喜喂。而且DFtest可以得到結(jié)論：數(shù)據(jù)在99%的置信度上是穩(wěn)定的瓤摧。
現(xiàn)在的疑問是如何做的衰減?halflife=12為何等于12?
Differencing:差分和分解
檢測和去除季節(jié)性
有兩種方法：
- 差分化： 以特定滯后數(shù)目的時刻的值的作差
- 分解： 對趨勢和季節(jié)性分別建模在移除它們
ts_log_diff = ts_log - ts_log.shift()
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
Results of Dickey-Fuller Test:
Test Statistic -2.717131
p-value 0.071121
#Lags Used 14.000000
Number of Observations Used 128.000000
Critical value (1%) -3.482501
Critical value (5%) -2.884398
Critical value (10%) -2.578960
dtype: float64
如圖，可以看出相比MA方法玉吁，Differencing方法處理后的數(shù)據(jù)的均值和方差的在時間軸上的振幅明顯縮小了照弥。DFtest的結(jié)論是在90%的置信度下，數(shù)據(jù)是穩(wěn)定的进副。下面去除趨勢和季節(jié)性數(shù)據(jù), 使用剩余的穩(wěn)定數(shù)據(jù)
from statsmodels.tsa.seasonal import seasonal_decompose
def decompose(timeseries):
# 返回包含三個部分 trend（趨勢部分） 这揣， seasonal（季節(jié)性部分） 和residual (殘留部分)
decomposition = seasonal_decompose(timeseries)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.subplot(411)
plt.plot(ts_log, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
return trend , seasonal, residual
#消除了trend 和seasonal之后，只對residual部分作為想要的時序數(shù)據(jù)進行處理
trend , seasonal, residual = decompose(ts_log)
residual.dropna(inplace=True)
test_stationarity(residual)
Results of Dickey-Fuller Test:
Test Statistic -6.332387e+00
p-value 2.885059e-08
#Lags Used 9.000000e+00
Number of Observations Used 1.220000e+02
Critical value (1%) -3.485122e+00
Critical value (5%) -2.885538e+00
Critical value (10%) -2.579569e+00
dtype: float64
如圖所示影斑，數(shù)據(jù)的均值和方差趨于常數(shù)给赞，幾乎無波動(看上去比之前的陡峭，但是要注意他的值域只有[-0.05,0.05]之間)矫户，所以直觀上可以認為是穩(wěn)定的數(shù)據(jù)片迅。另外DFtest的結(jié)果顯示，Statistic值原小于1%時的Critical value皆辽，所以在99%的置信度下柑蛇，數(shù)據(jù)是穩(wěn)定的。
對時序數(shù)據(jù)進行預(yù)測
假設(shè)經(jīng)過處理驱闷，已經(jīng)得到了穩(wěn)定時序數(shù)據(jù)耻台。接下來，我們使用ARIMA模型對數(shù)據(jù)已經(jīng)預(yù)測空另。ARIMA的介紹可以見本目錄下的另一篇文章盆耽。
step1： 通過ACF,PACF進行ARIMA（p，d扼菠，q）的p摄杂，q參數(shù)估計
由前文Differencing部分已知，一階差分后數(shù)據(jù)已經(jīng)穩(wěn)定循榆，所以d=1匙姜。所以用一階差分化的ts_log_diff = ts_log - ts_log.shift() 作為輸入。
等價于
ARIMA的預(yù)測模型可以表示為：
Y的預(yù)測值 = 常量c and/or 一個或多個最近時間的Y的加權(quán)和 and/or 一個或多個最近時間的預(yù)測誤差冯痢。
ARIMA模型有三個參數(shù):p,d,q氮昧。
- p--代表預(yù)測模型中采用的時序數(shù)據(jù)本身的滯后數(shù)(lags) ,也叫做AR/Auto-Regressive項
- d--代表時序數(shù)據(jù)需要進行幾階差分化，才是穩(wěn)定的浦楣，也叫Integrated項袖肥。(當前時間減上一個時間減幾次)
- q--代表預(yù)測模型中采用的預(yù)測誤差的滯后數(shù)(lags)，也叫做MA/Moving Average項
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
圖中振劳，上下兩條灰線之間是置信區(qū)間椎组，p的值就是ACF第一次穿過上置信區(qū)間時的橫軸值。q的值就是PACF第一次穿過上置信區(qū)間的橫軸值历恐。所以從圖中可以得到p=2寸癌，q=2专筷。
step2： 得到參數(shù)估計值p，d蒸苇，q之后磷蛹，生成模型ARIMA（p，d溪烤，q）
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(ts_log, order=(2, 1, 2))
results_ARIMA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))
Text(0.5, 1.0, 'RSS: 1.0292')
step3: 將模型代入原數(shù)據(jù)進行預(yù)測
因為上面的模型的擬合值是對原數(shù)據(jù)進行穩(wěn)定化之后的輸入數(shù)據(jù)的擬合味咳，所以需要對擬合值進行相應(yīng)處理的逆操作，使得它回到與原數(shù)據(jù)一致的尺度檬嘀。
#ARIMA擬合的其實是一階差分ts_log_diff槽驶，predictions_ARIMA_diff[i]是第i個月與i-1個月的ts_log的差值。
#由于差分化有一階滯后鸳兽，所以第一個月的數(shù)據(jù)是空的掂铐，
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print(predictions_ARIMA_diff.head())
#累加現(xiàn)有的diff，得到每個值與第一個月的差分（同log底的情況下）揍异。
#即predictions_ARIMA_diff_cumsum[i] 是第i個月與第1個月的ts_log的差值堡纬。
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
#先ts_log_diff => ts_log=>ts_log => ts
#先以ts_log的第一個值作為基數(shù)，復(fù)制給所有值蒿秦，然后每個時刻的值累加與第一個月對應(yīng)的差值(這樣就解決了烤镐，第一個月diff數(shù)據(jù)為空的問題了)
#然后得到了predictions_ARIMA_log => predictions_ARIMA
predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.figure()
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))
Month
1949-02-01 0.009580
1949-03-01 0.017491
1949-04-01 0.027670
1949-05-01 -0.004521
1949-06-01 -0.023890
dtype: float64
Text(0.5, 1.0, 'RMSE: 90.1046')
疑問
1.衰減因子如何作用?
許多情況下，可以認為越近的時刻越重要棍鳖。所以引入指數(shù)加權(quán)移動平均-- Exponentially-weighted moving average.（pandas中通過ewma()函數(shù)提供了此功能炮叶。）
# halflife的值決定了衰減因子alpha： alpha = 1 - exp(log(0.5) / halflife)
expweighted_avg = pd.ewma(ts_log,halflife=12)
ts_log_ewma_diff = ts_log - expweighted_avg
test_stationarity(ts_log_ewma_diff)
這里的衰減因子如何作用的呢?
2.ACF和PACF如何確認p和q的?
ACF(自相關(guān))和PACF(偏自相關(guān))的概念,以及對應(yīng)的pq值是如何求出的?
總結(jié)
- 獲取被觀測系統(tǒng)時間序列數(shù)據(jù)；
- 對數(shù)據(jù)繪圖渡处，觀測是否為平穩(wěn)時間序列镜悉；對于非平穩(wěn)時間序列要先進行d階差分運算，化為平穩(wěn)時間序列医瘫；
- 經(jīng)過第二步處理侣肄，已經(jīng)得到平穩(wěn)時間序列。要對平穩(wěn)時間序列分別求得其自相關(guān)系數(shù)ACF 和偏自相關(guān)系數(shù)PACF醇份，通過對自相關(guān)圖和偏自相關(guān)圖的分析稼锅，得到最佳的階層 p 和階數(shù) q
- 由以上得到的d、q僚纷、p矩距，得到ARIMA模型。然后開始對得到的模型進行模型檢驗怖竭。
注:
文章直接參考博客園博客:python時間序列
代碼和注釋有所修改,適配python新庫.自己動手做了下實驗
