今天給大家?guī)?lái)一篇實(shí)戰(zhàn)案例，本案例旨在運(yùn)用之前學(xué)習(xí)的時(shí)間序列分析和預(yù)測(cè)基礎(chǔ)理論知識(shí)，用一個(gè)基于交通數(shù)據(jù)的實(shí)際案例數(shù)據(jù)演示這些方法是如何被應(yīng)用的。

本文主要內(nèi)容

★ 首先使用探索性數(shù)據(jù)分析，從不同時(shí)間維度探索分析交通系統(tǒng)乘客數(shù)量。

★ 創(chuàng)建一個(gè)函數(shù)來(lái)檢查時(shí)間序列數(shù)據(jù)的平穩(wěn)性，通過(guò)一階差分將非平穩(wěn)性數(shù)據(jù)轉(zhuǎn)化為平穩(wěn)性數(shù)據(jù)。

★ 然后將數(shù)據(jù)分為訓(xùn)練集和驗(yàn)證集，簡(jiǎn)單介紹了并應(yīng)用多個(gè)時(shí)間序列預(yù)測(cè)技術(shù)，如樸素法、移動(dòng)平均方法、簡(jiǎn)單指數(shù)平滑、霍爾特線性趨勢(shì)法、霍爾特-溫特法、ARIMA和SARIMAX模型。

★ 最后使用SARIMAX模型預(yù)測(cè)未來(lái)7個(gè)月的流量，因?yàn)樗凶钚〉?a target="_blank">RMSE。如下圖所示，藍(lán)色線是訓(xùn)練數(shù)據(jù)，黃色線是驗(yàn)證數(shù)據(jù)，紅色是使用SARIMAX模型預(yù)測(cè)的數(shù)據(jù)。

雖然擬合最好的SARIMAX模型，但似乎也沒(méi)那么棒，當(dāng)然會(huì)有更好的方法來(lái)預(yù)測(cè)該數(shù)據(jù)。而本文重點(diǎn)是介紹這些基于統(tǒng)計(jì)學(xué)的經(jīng)典時(shí)間序列預(yù)測(cè)技術(shù)在實(shí)際案例中的應(yīng)用。

導(dǎo)入相關(guān)模塊

importpandasaspd
importnumpyasnp
importmatplotlib.pyplotasplt
fromdatetimeimportdatetime
frompandasimportSeries
fromsklearn.metricsimportmean_squared_error
frommathimportsqrt
fromstatsmodels.tsa.seasonalimportseasonal_decompose
importstatsmodels
importstatsmodels.apiassm
fromstatsmodels.tsa.arima_modelimportARIMA

數(shù)據(jù)集準(zhǔn)備

去直接讀取用pandas讀取csv文本文件，并拷貝一份以備用。

train=pd.read_csv("Train.csv")
test=pd.read_csv("Test.csv")

train_org=train.copy()
test_org=test.copy()

查看數(shù)據(jù)的列名

train.columns,test.columns

(Index(['ID', 'Datetime', 'Count'], dtype='object'),
 Index(['ID', 'Datetime'], dtype='object'))

查看數(shù)據(jù)類(lèi)型

train.dtypes,test.dtypes

(ID           int64
 Datetime    object
 Count        int64
 dtype: object,
 ID           int64
 Datetime    object
 dtype: object)

查看數(shù)據(jù)大小

test.shape,train.shape

((5112, 2), (18288, 3))

查看數(shù)據(jù)樣貌

train.head()

解析日期格式

train['Datetime']=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
test['Datetime']=pd.to_datetime(test.Datetime,format='%d-%m-%Y%H:%M')
test_org['Datetime']=pd.to_datetime(test_org.Datetime,format='%d-%m-%Y%H:%M')
train_org['Datetime']=pd.to_datetime(train_org.Datetime,format='%d-%m-%Y%H:%M')

時(shí)間日期格式解析結(jié)束后，記得查看下結(jié)果。

train.dtypes

ID                   int64
Datetime    datetime64[ns]
Count                int64
dtype: object

train.head()

時(shí)間序列數(shù)據(jù)的特征工程

時(shí)間序列的特征工程一般可以分為以下幾類(lèi)。本次案例我們根據(jù)實(shí)際情況，選用時(shí)間戳衍生時(shí)間特征。

時(shí)間戳雖然只有一列，但是也可以根據(jù)這個(gè)就衍生出很多很多變量了，具體可以分為三大類(lèi)：時(shí)間特征、布爾特征，時(shí)間差特征。

本案例首先對(duì)日期時(shí)間進(jìn)行時(shí)間特征處理，而時(shí)間特征包括年、季度、月、周、天(一年、一月、一周的第幾天)、小時(shí)、分鐘...

因?yàn)樾枰獙?duì)test, train, test_org, train_org四個(gè)數(shù)據(jù)框進(jìn)行同樣的處理，直接使用for循環(huán)批量提取年月日小時(shí)等特征。

foriin(test,train,test_org,train_org):
i['Year']=i.Datetime.dt.year
i['Month']=i.Datetime.dt.month
i['day']=i.Datetime.dt.day
i['Hour']=i.Datetime.dt.hour
#i["dayoftheweek"]=i.Datetime.dt.dayofweek

test.head()

時(shí)間戳衍生中，另一常用的方法為布爾特征，即：

是否年初/年末

是否月初/月末

是否周末

是否節(jié)假日

是否特殊日期

是否早上/中午/晚上

等等

下面判斷是否是周末，進(jìn)行特征衍生的布爾特征轉(zhuǎn)換。首先提取出日期時(shí)間的星期幾。

train['dayoftheweek']=train.Datetime.dt.dayofweek
#返回給定日期時(shí)間的星期幾

train.head()

再判斷day of the week是否是周末（星期六和星期日），是則返回1 ，否則返回0

defapplyer(row):
ifrow.dayofweek==5orrow.dayofweek==6:
return1
else:
return0
temp=train['Datetime']
temp2=train.Datetime.apply(applyer)
train['weekend']=temp2
train.index=train['Datetime']

對(duì)年月乘客總數(shù)統(tǒng)計(jì)后可視化，看看總體變化趨勢(shì)。

df=train.drop('ID',1)
ts=df['Count']
plt.plot(ts,label='Passengercount')

探索性數(shù)據(jù)分析

首先使用探索性數(shù)據(jù)分析，從不同時(shí)間維度探索分析交通系統(tǒng)乘客數(shù)量。

年

對(duì)年進(jìn)行聚合，求所有數(shù)據(jù)中按年計(jì)算的每日平均客流量，從圖中可以看出，隨著時(shí)間的增長(zhǎng)，每日平均客流量增長(zhǎng)迅速。

train.groupby('Year')['Count'].mean().plot.bar()

月

對(duì)月份進(jìn)行聚合，求所有數(shù)據(jù)中按月計(jì)算的每日平均客流量，從圖中可以看出，春夏季客流量每月攀升，而秋冬季客流量驟減。

train.groupby('Month')['Count'].mean().plot.bar()

年月

對(duì)年月份進(jìn)行聚合，求所有數(shù)據(jù)中按年月計(jì)算的每日平均客流量，從圖可知道，幾本是按照平滑指數(shù)上升的趨勢(shì)。

temp=train.groupby(['Year','Month'])['Count'].mean()
temp.plot()#乘客人數(shù)(每月)

日

對(duì)日進(jìn)行聚合，求所有數(shù)據(jù)中每月中的每日平均客流量。從圖中可大致看出，在5、11、24分別出現(xiàn)三個(gè)峰值，該峰值代表了上中旬的高峰期。

train.groupby('day')['Count'].mean(
).plot.bar(figsize=(15,5))

小時(shí)

對(duì)小時(shí)進(jìn)行聚合，求所有數(shù)據(jù)中一天內(nèi)按小時(shí)計(jì)算的平均客流量，得到了在中(12)晚(19)分別出現(xiàn)兩個(gè)峰值，該峰值代表了每日的高峰期。

train.groupby('Hour')['Count'].mean().plot.bar()

是否周末

對(duì)是否是周末進(jìn)行聚合，求所有數(shù)據(jù)中按是否周末計(jì)算的平均客流量，發(fā)現(xiàn)工作日比周末客流量客流量多近一倍，果然大家都是周末都喜歡宅在家里。

train.groupby('weekend')['Count'].mean().plot.bar()

周

對(duì)星期進(jìn)行聚合統(tǒng)計(jì)，求所有數(shù)據(jù)中按是周計(jì)算的平均客流量。

train.groupby('dayoftheweek')['Count'].mean().plot.bar()

時(shí)間重采樣

◎重采樣(resampling)指的是將時(shí)間序列從一個(gè)頻率轉(zhuǎn)換到另一個(gè)頻率的處理過(guò)程；
◎ 將高頻率數(shù)據(jù)聚合到低頻率稱(chēng)為降采樣(downsampling)；
◎ 將低頻率數(shù)據(jù)轉(zhuǎn)換到高頻率則稱(chēng)為升采樣(unsampling)；

train.head()

Pandas中的resample，重新采樣，是對(duì)原樣本重新處理的一個(gè)方法，是一個(gè)對(duì)常規(guī)時(shí)間序列數(shù)據(jù)重新采樣和頻率轉(zhuǎn)換的便捷的方法。

Resample方法的主要參數(shù)，如需要了解詳情，可以戳這里了解更多。

參數(shù)	說(shuō)明
freq	表示重采樣頻率，例如'M'、'5min'、Second(15)
how='mean'	用于產(chǎn)生聚合值的函數(shù)名或數(shù)組函數(shù)，例如'mean'、'ohlc'、np.max等，默認(rèn)是'mean'，其他常用的值由：'first'、'last'、'median'、'max'、'min'
axis=0	默認(rèn)是縱軸，橫軸設(shè)置axis=1

接下來(lái)對(duì)訓(xùn)練數(shù)據(jù)進(jìn)行對(duì)月、周、日及小時(shí)多重采樣。其實(shí)我們分月份進(jìn)行采樣，然后求月內(nèi)的均值。事實(shí)上重采樣，就相當(dāng)于groupby，只不過(guò)是根據(jù)月份這個(gè)period進(jìn)行分組。

train=train.drop('ID',1)
train.timestamp=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
train.index=train.timestamp

#每小時(shí)的時(shí)間序列
hourly=train.resample('H').mean()
#換算成日平均值
daily=train.resample('D').mean()
#換算成周平均值
weekly=train.resample('W').mean()
#換算成月平均值
monthly=train.resample('M').mean()

重采樣后對(duì)其進(jìn)行可視化，直觀地看看其變化趨勢(shì)。

對(duì)測(cè)試數(shù)據(jù)也進(jìn)行相同的時(shí)間重采樣處理。

test.Timestamp=pd.to_datetime(test.Datetime,format='%d-%m-%Y%H:%M')
test.index=test.Timestamp
#換算成日平均值
test=test.resample('D').mean()
train.Timestamp=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
train.index=train.Timestamp

#C換算成日平均值
train=train.resample('D').mean()

劃分訓(xùn)練集和驗(yàn)證集

到目前為止，我們有訓(xùn)練集和測(cè)試集，實(shí)際上，我們還需要一個(gè)驗(yàn)證集，用來(lái)實(shí)時(shí)驗(yàn)證和調(diào)整訓(xùn)練模型。下面直接用索引切片的方式做處理。

Train=train.loc['2012-08-25':'2014-06-24']
valid=train['2014-06-25':'2014-09-25']

劃分好數(shù)據(jù)集后，繪制折線圖將訓(xùn)練集和驗(yàn)證集進(jìn)行可視化。

模型建立

數(shù)據(jù)準(zhǔn)備好了，就到了模型建立階段，這里我們應(yīng)用多個(gè)時(shí)間序列預(yù)測(cè)技術(shù)，如樸素法、移動(dòng)平均方法、簡(jiǎn)單指數(shù)平滑、霍爾特線性趨勢(shì)法、霍爾特-溫特法、ARIMA和SARIMAX模型。

樸素預(yù)測(cè)法

如果數(shù)據(jù)集在一段時(shí)間內(nèi)都很穩(wěn)定，我們想預(yù)測(cè)第二天的價(jià)格，可以取前面一天的價(jià)格，預(yù)測(cè)第二天的值。這種假設(shè)第一個(gè)預(yù)測(cè)點(diǎn)和上一個(gè)觀察點(diǎn)相等的預(yù)測(cè)方法就叫樸素預(yù)測(cè)法(Naive Forecast)，即。

因?yàn)闃闼仡A(yù)測(cè)法用最近的觀測(cè)值作為預(yù)測(cè)值，因此他最簡(jiǎn)單的預(yù)測(cè)方法。雖然樸素預(yù)測(cè)法并不是一個(gè)很好的預(yù)測(cè)方法，但是它可以為其他預(yù)測(cè)方法提供一個(gè)基準(zhǔn)。

dd=np.asarray(Train.Count)
#將結(jié)構(gòu)數(shù)據(jù)轉(zhuǎn)化為ndarray
y_hat=valid.copy()
y_hat['naive']=dd[len(dd)-1]
plt.plot(Train.index,Train['Count'],label='Train')
plt.plot(valid.index,valid['Count'],label='Valid')
plt.plot(y_hat.index,y_hat['naive'],label='NaiveForecast')

模型評(píng)價(jià)

用RMSE檢驗(yàn)樸素法的的準(zhǔn)確率

rms=sqrt(mean_squared_error(valid.Count,y_hat.naive))
print(rms)

111.79050467496724

移動(dòng)平均值法

移動(dòng)平均法也叫滑動(dòng)平均法，取前面n個(gè)點(diǎn)的平均值作為預(yù)測(cè)值。

計(jì)算移動(dòng)平均值涉及到一個(gè)有時(shí)被稱(chēng)為"滑動(dòng)窗口"的大小值。使用簡(jiǎn)單的移動(dòng)平均模型，我們可以根據(jù)之前數(shù)值的固定有限數(shù)的平均值預(yù)測(cè)某個(gè)時(shí)序中的下一個(gè)值。利用一個(gè)簡(jiǎn)單的移動(dòng)平均模型，我們預(yù)測(cè)一個(gè)時(shí)間序列中的下一個(gè)值是基于先前值的固定有限個(gè)數(shù)“p”的平均值。

這樣，對(duì)于所有的

#最近10次觀測(cè)的移動(dòng)平均值，即滑動(dòng)窗口大小為P=10
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(10).mean().iloc[-1]

#最近20次觀測(cè)的移動(dòng)平均值
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(20).mean().iloc[-1]

#最近30次觀測(cè)的移動(dòng)平均值
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(50).mean().iloc[-1]
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['moving_avg_forecast'],
label='MovingAverageForecastusing50observations')

簡(jiǎn)單指數(shù)平滑法

介紹這個(gè)之前，需要知道什么是簡(jiǎn)單平均法(Simple Average），該方法預(yù)測(cè)的期望值等于所有先前觀測(cè)點(diǎn)的平均值。

物品價(jià)格會(huì)隨機(jī)上漲和下跌，平均價(jià)格會(huì)保持一致。我們經(jīng)常會(huì)遇到一些數(shù)據(jù)集，雖然在一定時(shí)期內(nèi)出現(xiàn)小幅變動(dòng)，但每個(gè)時(shí)間段的平均值確實(shí)保持不變。這種情況下，我們可以認(rèn)為第二天的價(jià)格大致和過(guò)去的平均價(jià)格值一致。

簡(jiǎn)單平均法和加權(quán)移動(dòng)平均法在選取時(shí)間點(diǎn)的思路上存在較大的差異：

簡(jiǎn)單平均法將過(guò)去數(shù)據(jù)一個(gè)不漏地全部加以同等利用；

移動(dòng)平均法則不考慮較遠(yuǎn)期的數(shù)據(jù)，并在加權(quán)移動(dòng)平均法中給予近期更大的權(quán)重。

我們就需要在這兩種方法之間取一個(gè)折中的方法，在將所有數(shù)據(jù)考慮在內(nèi)的同時(shí)也能給數(shù)據(jù)賦予不同非權(quán)重。

簡(jiǎn)單指數(shù)平滑法 (Simple Exponential Smoothing)相比更早時(shí)期內(nèi)的觀測(cè)值，越近的觀測(cè)值會(huì)被賦予更大的權(quán)重，而時(shí)間越久遠(yuǎn)的權(quán)重越小。它通過(guò)加權(quán)平均值計(jì)算出預(yù)測(cè)值，其中權(quán)重隨著觀測(cè)值從早期到晚期的變化呈指數(shù)級(jí)下降，最小的權(quán)重和最早的觀測(cè)值相關(guān)：

其中是平滑參數(shù)。

fromstatsmodels.tsa.apiimportExponentialSmoothing,SimpleExpSmoothing,Holt
y_hat_avg=valid.copy()
fit2=SimpleExpSmoothing(np.asarray(Train['Count'])).fit(smoothing_level=0.6,optimized=False)
y_hat_avg['SES']=fit2.forecast(len(valid))
plt.figure(figsize=(16,8))
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['SES'],label='SES')
plt.legend(loc='best')
plt.show()

模型評(píng)價(jià)

用RMSE檢驗(yàn)樸素法的的準(zhǔn)確率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.SES))
print(rms)

113.43708111884514

霍爾特線性趨勢(shì)法

Holts線性趨勢(shì)模型，該方法考慮了數(shù)據(jù)集的趨勢(shì)，即序列的增加或減少性質(zhì)。

盡管這些方法中的每一種都可以應(yīng)用趨勢(shì)：簡(jiǎn)單平均法會(huì)假設(shè)最后兩點(diǎn)之間的趨勢(shì)保持不變，或者我們可以平均所有點(diǎn)之間的所有斜率以獲得平均趨勢(shì)，使用移動(dòng)趨勢(shì)平均值或應(yīng)用指數(shù)平滑。但我們需要一種無(wú)需任何假設(shè)就能準(zhǔn)確繪制趨勢(shì)圖的方法。這種考慮數(shù)據(jù)集趨勢(shì)的方法稱(chēng)為霍爾特線性趨勢(shì)法，或者霍爾特指數(shù)平滑法。

y_hat_avg=valid.copy()
fit1=Holt(np.asarray(Train['Count'])
).fit(smoothing_level=0.3,smoothing_slope=0.1)
y_hat_avg['Holt_linear']=fit1.forecast(len(valid))

plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['Holt_linear'],label='Holt_linear')

模型評(píng)價(jià)

用RMSE檢驗(yàn)樸素法的的準(zhǔn)確率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.Holt_linear))
print(rms)

112.94278345314041

由于holts線性趨勢(shì)，到目前為止具有最好的準(zhǔn)確性，我們嘗試使用它來(lái)預(yù)測(cè)測(cè)試數(shù)據(jù)集。

predict=fit1.forecast(len(test))
test['prediction']=predict
#計(jì)算每小時(shí)計(jì)數(shù)的比率
train_org['ratio']=train_org['Count']/train_org['Count'].sum()
#按小時(shí)計(jì)數(shù)分組
temp=train_org.groupby(['Hour'])['ratio'].sum()
#保存聚合后的數(shù)據(jù)
pd.DataFrame(temp,columns=['ratio']).to_csv('GROUPBY.csv')

temp2=pd.read_csv('GROUPBY.csv')
#按日、月、年合并test和test_org
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge['ID']=merge['ID_y']
merge.head()

刪除不需要的特征。

merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y','ID_x','ID_y'],axis=1)
merge.head()

通過(guò)合并merge和temp2進(jìn)行預(yù)測(cè)。

prediction=pd.merge(merge,temp2,on='Hour',how='left')
#將比率轉(zhuǎn)換成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
submission=prediction
pd.DataFrame(submission,columns=['ID','Count']).to_csv('Holt_Linear.csv')

霍爾特-溫特法

霍爾特-溫特（Holt-Winters）方法，在 Holt模型基礎(chǔ)上引入了 Winters 周期項(xiàng)（也叫做季節(jié)項(xiàng)），可以用來(lái)處理月度數(shù)據(jù)（周期 12）、季度數(shù)據(jù)（周期 4）、星期數(shù)據(jù)（周期 7）等時(shí)間序列中的固定周期的波動(dòng)行為。引入多個(gè) Winters 項(xiàng)還可以處理多種周期并存的情況。

#HoltsWintermodel
y_hat_avg=valid.copy()
fit1=ExponentialSmoothing(np.asarray(Train['Count']),seasonal_periods=7,trend='add',seasonal='add',).fit()
y_hat_avg['Holts_Winter']=fit1.forecast(len(valid))
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['Holts_Winter'],label='Holt_Winter')

模型評(píng)價(jià)

用RMSE檢驗(yàn)樸素法的的準(zhǔn)確率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.Holts_Winter))
print(rms)

82.37292653831038

模型預(yù)測(cè)

predict=fit1.forecast(len(test))
test['prediction']=predict
#按日、月、年合并Test和test_original
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y'],axis=1)
#通過(guò)合并merge和temp2進(jìn)行預(yù)測(cè)
prediction=pd.merge(merge,temp2,on='Hour',how='left')
#將比率轉(zhuǎn)換成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
prediction['ID']=prediction['ID_y']
submission=prediction.drop(['day','Hour','ratio','prediction','ID_x','ID_y'],axis=1)
#轉(zhuǎn)換最終提交的csv格式
pd.DataFrame(submission,columns=['ID','Count']).to_csv('Holtwinters.csv')

迪基-福勒檢驗(yàn)

函數(shù)執(zhí)行迪基-福勒檢驗(yàn)以確定數(shù)據(jù)是否為平穩(wěn)時(shí)間序列。

在統(tǒng)計(jì)學(xué)里，迪基-福勒檢驗(yàn)（Dickey-Fuller test）可以測(cè)試一個(gè)自回歸模型是否存在單位根（unit root）。回歸模型可以寫(xiě)為，是一階差分。測(cè)試是否存在單位根等同于測(cè)試是否。

因?yàn)榈匣?福勒檢驗(yàn)測(cè)試的是殘差項(xiàng)，并非原始數(shù)據(jù)，所以不能用標(biāo)準(zhǔn)t統(tǒng)計(jì)量。我們需要用迪基-福勒統(tǒng)計(jì)量。

fromstatsmodels.tsa.stattoolsimportadfuller
deftest_stationary(timeseries):
#確定滾動(dòng)數(shù)據(jù)
rolmean=timeseries.rolling(24).mean()
rolstd=timeseries.rolling(24).std()

#會(huì)議滾動(dòng)數(shù)據(jù)
orig=plt.plot(timeseries,color='blue',label='Original')
mean=plt.plot(rolmean,color='red',label='RollingMean')
std=plt.plot(rolstd,color='black',label='RollingStd')
plt.legend(loc='best')
plt.title('RollingMean&StandardDeviation')
plt.show(block=False)

#執(zhí)行迪基-福勒檢驗(yàn)
print('ResultsofDickey-FullerTest:')
dftest=adfuller(timeseries,autolag='AIC')
dfoutput=pd.Series(dftest[0:4],index=['TestStatistic','P-value','#lagsused','NoofObservationsused'])
forkey,valueindftest[4].items():
dfoutput['CriticalValue(%s)'%key]=value
print(dfoutput)

繪制檢驗(yàn)圖

test_stationary(train_org['Count'])

Results of Dickey-Fuller Test: 
Test Statistic                -4.456561
P-value                        0.000235
#lags used                    45.000000
No of Observations used    18242.000000
Critical Value (1%)           -3.430709
Critical Value (5%)           -2.861698
Critical Value (10%)          -2.566854
dtype: float64

檢驗(yàn)統(tǒng)計(jì)數(shù)據(jù)表明，由于p值小于0.05，數(shù)據(jù)是平穩(wěn)的。

移動(dòng)平均值

在統(tǒng)計(jì)學(xué)中，移動(dòng)平均（moving average），又稱(chēng)滑動(dòng)平均是一種通過(guò)創(chuàng)建整個(gè)數(shù)據(jù)集中不同子集的一系列平均數(shù)來(lái)分析數(shù)據(jù)點(diǎn)的計(jì)算方法。移動(dòng)平均通常與時(shí)間序列數(shù)據(jù)一起使用，以消除短期波動(dòng)，突出長(zhǎng)期趨勢(shì)或周期。

對(duì)原始數(shù)據(jù)求對(duì)數(shù)。

Train_log=np.log(Train['Count'])
valid_log=np.log(Train['Count'])
Train_log.head()

Datetime
2012-08-25    1.152680
2012-08-26    1.299283
2012-08-27    0.949081
2012-08-28    0.882389
2012-08-29    0.916291
Freq: D, Name: Count, dtype: float64

繪制移動(dòng)平均值曲線

moving_avg=Train_log.rolling(24).mean()
plt.plot(Train_log)
plt.plot(moving_avg,color='red')
plt.show()

去除移動(dòng)平均值后再進(jìn)行迪基-福勒檢驗(yàn)

train_log_moving_avg_diff=Train_log-moving_avg
train_log_moving_avg_diff.dropna(inplace=True)
test_stationary(train_log_moving_avg_diff)

Results of Dickey-Fuller Test: 
Test Statistic            -5.861646e+00
P-value                    3.399422e-07
#lags used                 2.000000e+01
No of Observations used    6.250000e+02
Critical Value (1%)       -3.440856e+00
Critical Value (5%)       -2.866175e+00
Critical Value (10%)      -2.569239e+00
dtype: float64

對(duì)數(shù)時(shí)序數(shù)據(jù)求二階差分后再迪基-福勒檢驗(yàn)

train_log_diff=Train_log-Train_log.shift(1)
test_stationary(train_log_diff.dropna())

Results of Dickey-Fuller Test: 
Test Statistic            -8.237568e+00
P-value                    5.834049e-13
#lags used                 1.900000e+01
No of Observations used    6.480000e+02
Critical Value (1%)       -3.440482e+00
Critical Value (5%)       -2.866011e+00
Critical Value (10%)      -2.569151e+00
dtype: float64

季節(jié)性分解

對(duì)進(jìn)行對(duì)數(shù)轉(zhuǎn)換后的原始數(shù)據(jù)進(jìn)行季節(jié)性分解。

decomposition=seasonal_decompose(
pd.DataFrame(Train_log).Count.values,freq=24)
trend=decomposition.trend
seasonal=decomposition.seasonal
residual=decomposition.resid

plt.plot(Train_log,label='Original')
plt.plot(trend,label='Trend')
plt.plot(seasonal,label='Seasonality')
plt.plot(residual,label='Residuals')

對(duì)季節(jié)性分解后的殘差數(shù)據(jù)進(jìn)行迪基-福勒檢驗(yàn)

train_log_decompose=pd.DataFrame(residual)
train_log_decompose['date']=Train_log.index
train_log_decompose.set_index('date',inplace=True)
train_log_decompose.dropna(inplace=True)
test_stationary(train_log_decompose[0])

Results of Dickey-Fuller Test: 
Test Statistic            -7.822096e+00
P-value                    6.628321e-12
#lags used                 2.000000e+01
No of Observations used    6.240000e+02
Critical Value (1%)       -3.440873e+00
Critical Value (5%)       -2.866183e+00
Critical Value (10%)      -2.569243e+00
dtype: float64

自相關(guān)和偏自相關(guān)圖

fromstatsmodels.tsa.stattoolsimportacf,pacf
lag_acf=acf(train_log_diff.dropna(),nlags=25)
lag_pacf=pacf(train_log_diff.dropna(),nlags=25,method='ols')
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')

AR模型

AR模型訓(xùn)練及預(yù)測(cè)

model=ARIMA(Train_log,order=(2,1,0))
#這里q值是零，因?yàn)樗皇茿R模型
results_AR=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_AR.fittedvalues,color='red',label='predictions')

AR_predict=results_AR.predict(start="2014-06-25",end="2014-09-25")
AR_predict=AR_predict.cumsum().shift().fillna(0)
AR_predict1=pd.Series(np.ones(valid.shape[0])*np.log(valid['Count'])[0],index=valid.index)
AR_predict1=AR_predict1.add(AR_predict,fill_value=0)
AR_predict=np.exp(AR_predict1)
plt.plot(valid['Count'],label="Valid")
plt.plot(AR_predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(AR_predict,valid['Count']))/valid.shape[0]))

MA模型

model=ARIMA(Train_log,order=(0,1,2))
#這里的p值是零，因?yàn)樗皇荕A模型
results_MA=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_MA.fittedvalues,color='red',label='prediction')

MA_predict=results_MA.predict(start="2014-06-25",end="2014-09-25")
MA_predict=MA_predict.cumsum().shift().fillna(0)
MA_predict1=pd.Series(np.ones(valid.shape[0])*np.log(valid['Count'])[0],index=valid.index)
MA_predict1=MA_predict1.add(MA_predict,fill_value=0)
MA_predict=np.exp(MA_predict1)

plt.plot(valid['Count'],label="Valid")
plt.plot(MA_predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(MA_predict,valid['Count']))/valid.shape[0]))

ARMA模型

model=ARIMA(Train_log,order=(2,1,2))
results_ARIMA=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_ARIMA.fittedvalues,color='red',label='predicted')

defcheck_prediction_diff(predict_diff,given_set):
predict_diff=predict_diff.cumsum().shift().fillna(0)
predict_base=pd.Series(np.ones(given_set.shape[0])*np.log(given_set['Count'])[0],index=given_set.index)
predict_log=predict_base.add(predict_diff,fill_value=0)
predict=np.exp(predict_log)

plt.plot(given_set['Count'],label="Givenset")
plt.plot(predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(predict,given_set['Count']))/given_set.shape[0]))

defcheck_prediction_log(predict_log,given_set):
predict=np.exp(predict_log)

plt.plot(given_set['Count'],label="Givenset")
plt.plot(predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(predict,given_set['Count']))/given_set.shape[0]))
plt.show()

ARIMA_predict_diff=results_ARIMA.predict(start="2014-06-25",
end="2014-09-25")
check_prediction_diff(ARIMA_predict_diff,valid)

SARIMAX模型

y_hat_avg=valid.copy()
fit1=sm.tsa.statespace.SARIMAX(Train.Count,order=(2,1,4),seasonal_order=(0,1,1,7)).fit()
y_hat_avg['SARIMA']=fit1.predict(start="2014-6-25",end="2014-9-25",dynamic=True)
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['SARIMA'],label='SARIMA')

模型評(píng)價(jià)

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.SARIMA))
print(rms)

70.26240839723575

預(yù)測(cè)

predict=fit1.predict(start="2014-9-26",end="2015-4-26",dynamic=True)
test['prediction']=predict
#按日、月、年合并Test和test_original
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y'],axis=1
#通過(guò)合并merge和temp2進(jìn)行預(yù)測(cè)
prediction=pd.merge(merge,temp2,on='Hour',how='left')
#將比率轉(zhuǎn)換成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
prediction['ID']=prediction['ID_y']
submission=prediction.drop(['day','Hour','ratio','prediction','ID_x','ID_y'],axis=1)
#轉(zhuǎn)換最終提交的csv格式
pd.DataFrame(submission,columns=['ID','Count']).to_csv('SARIMAX.csv')

原文標(biāo)題：時(shí)間序列分析和預(yù)測(cè)實(shí)戰(zhàn)

文章出處：【微信公眾號(hào)：數(shù)據(jù)分析與開(kāi)發(fā)】歡迎添加關(guān)注！文章轉(zhuǎn)載請(qǐng)注明出處。