時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

本文作者： AI研習(xí)社

2017-02-17 12:08

導(dǎo)語：專家手把手帶你做時(shí)間序列預(yù)測練習(xí)。

編者按：本文是澳大利亞知名機(jī)器學(xué)習(xí)專家 Jason Brownlee 撰寫的教程，極其全面細(xì)致，一步步向讀者解釋如何操作，以及為什么這么做。雷鋒網(wǎng)整理編譯，特與大家分享。更多AI開發(fā)技術(shù)文章，請關(guān)注AI研習(xí)社（微信號：okweiwu）。

Jason Brownlee：時(shí)間序列預(yù)測法是一個(gè)過程，而獲得良好預(yù)測結(jié)果的唯一途徑是實(shí)踐這個(gè)過程。

在本教程中，您將了解如何利用Python語言來預(yù)測波士頓每月持械搶劫案發(fā)生的數(shù)量。

本教程所述為您提供了一套處理時(shí)間序列預(yù)測問題的框架，包括方法步驟和工具，通過實(shí)踐，可以用它來解決自己遇到的相關(guān)問題。

本教程結(jié)束之后，您將了解：

如何核查Python環(huán)境并準(zhǔn)確地定義一個(gè)時(shí)間序列預(yù)測問題。
如何構(gòu)建一套測試工具鏈，用于評估模型，開發(fā)預(yù)測原型。以及如何通過時(shí)間序列分析工具更好地理解你的問題。
如何開發(fā)自回歸積分滑動(dòng)平均模型（ARIMA），將其保存到文件，并在之后加載它對新的時(shí)間步驟進(jìn)行預(yù)測。

讓我們開始吧。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

波士頓

概述

在本教程中，我們將端到端地來解析一個(gè)時(shí)間序列預(yù)測工程，從下載數(shù)據(jù)集、定義問題到訓(xùn)練出最終模型并進(jìn)行預(yù)測。

該工程并不面面俱到，但展示了如何通過系統(tǒng)性地處理時(shí)間序列預(yù)測問題，來快速獲得好的結(jié)果。

我們將通過如下步驟來解析該工程：

環(huán)境配置.
問題描述
測試工具鏈
持續(xù)性.
數(shù)據(jù)分析
ARIMA 模型
模型驗(yàn)證

該部分提供一個(gè)解決時(shí)間序列預(yù)測問題的模板，你可以套用自己的數(shù)據(jù)集。

1. 環(huán)境配置

本教程假設(shè)在已安裝好了SciPy及其相關(guān)依賴庫的環(huán)境下進(jìn)行，包括：

SciPy
NumPy
Matplotlib
Pandas
scikit-learn
statsmodels

我使用的是Python2.7。該腳本將幫助你確認(rèn)你安裝這些庫的版本。

# scipy
import scipy
print('scipy: {}'.format(scipy.__version__))
# numpy
import numpy
print('numpy: {}'.format(numpy.__version__))
# matplotlib
import matplotlib
print('matplotlib: {}'.format(matplotlib.__version__))
# pandas
import pandas
print('pandas: {}'.format(pandas.__version__))
# scikit-learn
import sklearn
print('sklearn: {}'.format(sklearn.__version__))
# statsmodels
import statsmodels
print('statsmodels: {}'.format(statsmodels.__version__))

我編寫該教程的所用的開發(fā)環(huán)境顯示的結(jié)果如下：

scipy: 0.18.1
numpy: 1.11.2
matplotlib: 1.5.3
pandas: 0.19.1
sklearn: 0.18.1
statsmodels: 0.6.1

2. 問題描述

我們研究的問題是預(yù)測美國波士頓每月持械搶劫案的數(shù)量。

該數(shù)據(jù)集提供了從1966年1月到1975年10月在波士頓每月發(fā)生的的武裝搶劫案的數(shù)量，或者說，是這十年間的數(shù)據(jù)。

這些值是一些記錄的觀察計(jì)數(shù)，總有118個(gè)觀察值。

數(shù)據(jù)集來自于McCleary＆Hay（1980）。

您可以了解有關(guān)此數(shù)據(jù)集的更多信息，并直接從 DataMarket 下載。

下載地址：https://datamarket.com/data/set/22ob/monthly-boston-armed-robberies-jan1966-oct1975-deutsch-and-alt-1977#!ds=22ob&display=line

將數(shù)據(jù)集下載為CSV文件，并將其放在當(dāng)前工作目錄中，文件名為“robberies.csv”。

3.測試工具鏈

我們必須開發(fā)一套測試工具鏈來審查數(shù)據(jù)和評估候選模型。

這包含兩個(gè)步驟：

1.定義好驗(yàn)證數(shù)據(jù)集。

2.開發(fā)一套模型評估方法。

3.1 驗(yàn)證數(shù)據(jù)集

這個(gè)數(shù)據(jù)集不是當(dāng)前時(shí)間段的數(shù)據(jù)集。這意味著我們不能輕輕松松收集更新后的數(shù)據(jù)來驗(yàn)證該模型。

因此，我們假設(shè)現(xiàn)在是1974年10月，并將分析和模型選擇的數(shù)據(jù)集里扣除最后一年的數(shù)據(jù)。

最后一年的數(shù)據(jù)將用于驗(yàn)證生成的最終模型。

下面的代碼會(huì)加載該數(shù)據(jù)集為Pandas序列對象并將其切分成兩部分，一個(gè)用于模型開發(fā)（dataset.csv）和另一個(gè)用于驗(yàn)證（validation.csv）。

from pandas import Series
series = Series.from_csv('robberies.csv', header=0)
split_point = len(series) - 12
dataset, validation = series[0:split_point], series[split_point:]
print('Dataset %d, Validation %d' % (len(dataset), len(validation)))
dataset.to_csv('dataset.csv')
validation.to_csv('validation.csv')

運(yùn)行該示例將會(huì)創(chuàng)建兩個(gè)文件并打印出每個(gè)文件中的觀察數(shù)。

1. Dataset 106, Validation 12

這些文件的具體內(nèi)容是：

dataset.csv：1966年1月至1974年10月的觀察（106次觀察）
validation.csv：1974年11月至1975年10月的觀察（12次觀察）

驗(yàn)證數(shù)據(jù)集大小是原始數(shù)據(jù)集的10％。

請注意，由于保存的數(shù)據(jù)集中沒有標(biāo)題行，因此，稍后處理這些文件時(shí)，我們也不需要滿足此要求。

3.2. 模型驗(yàn)證

模型評估將僅對上一節(jié)中準(zhǔn)備的dataset.csv中的數(shù)據(jù)進(jìn)行處理。

模型評估涉及兩個(gè)要素：

1. 評價(jià)指標(biāo)。

2. 測試策略。

3.2.1評價(jià)指標(biāo)

這些觀察值是搶劫案的計(jì)數(shù)。

我們將使用均方根誤差（RMSE）的方式來評價(jià)預(yù)測的好壞。這將給予錯(cuò)誤的預(yù)測更多的權(quán)重，并且將具有與原始數(shù)據(jù)相同的單位。

在RMSE進(jìn)行計(jì)算和報(bào)告之前，必須反轉(zhuǎn)對數(shù)據(jù)的所有轉(zhuǎn)換，以使不同方法的效果可直接比較。

我們可以使用scikit-learn庫的輔助函數(shù)mean_squared_error（）來執(zhí)行RMSE計(jì)算，該函數(shù)計(jì)算出了預(yù)期值列表（測試集）和預(yù)測列表之間的均方誤差。然后我們可以取這個(gè)值的平方根來作為一個(gè)RMSE分?jǐn)?shù)。

例如：

from sklearn.metrics import mean_squared_error
from math import sqrt
...
test = ...
predictions = ...
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

3.2.2 測試策略

將使用步進(jìn)驗(yàn)證的方式(walk-forward)來評估候選模型。

這是因?yàn)樵搯栴}定義需要滾動(dòng)式預(yù)測的模型。給定所有可用數(shù)據(jù)，這需要逐步進(jìn)行預(yù)測。

步進(jìn)驗(yàn)證的方式將按照如下幾步執(zhí)行：

數(shù)據(jù)集的前50％將被用來訓(xùn)練模型。
剩余的50％的數(shù)據(jù)集將被迭代并測試模型。
對于測試數(shù)據(jù)集中的每個(gè)步驟：
訓(xùn)練該模型
利用該模型預(yù)測一次，并經(jīng)預(yù)測結(jié)果保存下來用于之后的驗(yàn)證
測試數(shù)據(jù)集的數(shù)據(jù)作為下一個(gè)迭代的訓(xùn)練集數(shù)據(jù)
對在測試數(shù)據(jù)集迭代期間進(jìn)行的預(yù)測結(jié)果進(jìn)行評估，并報(bào)告RMSE得分。

由于較小的數(shù)據(jù)規(guī)模，我們將允許在每次預(yù)測之前，在給定所有可用數(shù)據(jù)的情況下重新訓(xùn)練模型。

我們可以使用簡單的NumPy和Python來編寫測試工具的代碼。

首先，我們可以將數(shù)據(jù)集直接分成訓(xùn)練和測試集。我們總要將加載的數(shù)據(jù)轉(zhuǎn)換為float32類型，以防止仍然有一些數(shù)據(jù)為String或Integer數(shù)據(jù)類型。

# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]

接下來，我們可以時(shí)間步長遍歷測試數(shù)據(jù)集。訓(xùn)練數(shù)據(jù)集存儲在一個(gè)Python的list對象中，因?yàn)槲覀冃枰诿看蔚泻苋菀椎母郊右恍┯^察值，Numpy數(shù)組連接有些overkill。

模型進(jìn)行的預(yù)測，按慣例稱為yhat。這是由于結(jié)果或觀察被稱為y，而yhat（有上標(biāo)“^”的y）是用于預(yù)測y變量的數(shù)學(xué)符號。

打印預(yù)測和觀察結(jié)果，對每個(gè)觀察值進(jìn)行完整性檢查，以防模型存在問題。

# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# predict
yhat = ...
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))

4. Persistence(持續(xù)性模型)

在數(shù)據(jù)分析和建模陷入困境之前的第一步是建立一個(gè)評估原型。

這里將為利用測試工具進(jìn)行模型評估和性能測量，分別提供一套模版。通過該性能測量，可以將當(dāng)前模型和其他更復(fù)雜的預(yù)測模型進(jìn)行比較。

時(shí)間序列預(yù)測的預(yù)測原型被稱為天真預(yù)測或Persistence。

在這里，先前時(shí)間步驟的觀察結(jié)果將被用作于下一時(shí)間步長的預(yù)測。

我們可以將其直接插入上一節(jié)中定義的測試工具中。

下面羅列了完整的代碼：

from pandas import Series
from sklearn.metrics import mean_squared_error
from math import sqrt
# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# predict
yhat = history[-1]
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

運(yùn)行測試工具，打印測試數(shù)據(jù)集中每次迭代的預(yù)測值和實(shí)際觀察值。

示例到進(jìn)行打印模型的RMSE這一步結(jié)束。

在該示例中，我們可以看到持久性模型實(shí)現(xiàn)了51.844的RMSE。這意味著，對于每個(gè)預(yù)測，該模型中平均有51次搶劫值預(yù)測是錯(cuò)誤的。

...
>Predicted=241.000, Expected=287
>Predicted=287.000, Expected=355
>Predicted=355.000, Expected=460
>Predicted=460.000, Expected=364
>Predicted=364.000, Expected=487
RMSE: 51.844

我們現(xiàn)在有了一個(gè)預(yù)測方法和性能評估的原型; 現(xiàn)在我們可以開始深入數(shù)據(jù)。

5. 數(shù)據(jù)分析

我們可以使用統(tǒng)計(jì)匯總和數(shù)據(jù)的繪圖，來快速了解該預(yù)測問題的數(shù)據(jù)結(jié)構(gòu)。

在本節(jié)中，我們將從四個(gè)角度來看數(shù)據(jù)：

1.摘要統(tǒng)計(jì)。

2.線圖。

3.密度圖。

4.盒型-虛線圖。

5.1 摘要統(tǒng)計(jì)

利用文本編輯器中打開數(shù)據(jù)dataset.csv文件或原始的robberies.csv文件，并查看數(shù)據(jù)。

快速檢查表明，數(shù)據(jù)集里沒有明顯丟失的觀察數(shù)據(jù)。我們可能已經(jīng)更早的發(fā)現(xiàn)了像NaN或'？'等值數(shù)，如果我們試圖強(qiáng)制序列化數(shù)據(jù)為浮點(diǎn)值。

摘要統(tǒng)計(jì)提供了觀察值的快速預(yù)覽。它可以幫助我們快速了解我們正在使用的東西。

以下示例為計(jì)算并打印時(shí)間系列的摘要統(tǒng)計(jì)信息。

from pandas import Series
series = Series.from_csv('dataset.csv')
print(series.describe())

運(yùn)行該示例會(huì)提供大量的摘要信息供回顧。

這些統(tǒng)計(jì)數(shù)據(jù)包含的一些信息為：

觀察次數(shù)（計(jì)數(shù)）與我們的期望相符，這意味著我們將數(shù)據(jù)正確處理了。
平均值約為173，我們可以想象出我們在這個(gè)序列中的級別。
在112次搶劫時(shí)，標(biāo)準(zhǔn)偏差（平均值的平均分布）相對較大。
百分位數(shù)與標(biāo)準(zhǔn)差一致表明數(shù)據(jù)有很大的差異。

count    106.000000
mean     173.103774
std      112.231133
min       29.000000
25%       74.750000
50%      144.500000
75%      271.750000
max      487.000000

對于序列中比較大的變化將可能模型很難進(jìn)行高度準(zhǔn)確的預(yù)測，如果它是由隨機(jī)波動(dòng)（例如，非系統(tǒng)的）引起，

5.2 線圖

畫一個(gè)時(shí)間序列的線圖可以對該問題提供很多深入的信息。

下面的示例創(chuàng)建并顯示出數(shù)據(jù)集的線圖。

from pandas import Series
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
series.plot()
pyplot.show()

運(yùn)行示例并查看繪圖。注意系列中的任何明顯的時(shí)間結(jié)構(gòu)。

圖中的一些觀察結(jié)果包括：
隨著時(shí)間的推移，劫案有增加的趨勢。
似乎沒有任何明顯的異常值。
每年都有相對較大的波動(dòng)，上下波動(dòng)。
后幾年的波動(dòng)似乎大于前幾年的波動(dòng)。
整體趨勢表明該數(shù)據(jù)集幾乎肯定不是穩(wěn)定的，并且波形的明顯變化也可能對此有所提示。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

每月波士頓搶劫案數(shù)據(jù)線圖

這些簡單的觀察結(jié)果表明，我們可以看到對趨勢進(jìn)行建模和從時(shí)間序列中將它提取出來的好處。

或者，我們可以使用差分法從而使序列的建模平穩(wěn)。如果后期波動(dòng)有增長趨勢，我們甚至可能需要兩個(gè)差分到不同級別。

5.3 密度圖

提取出密度圖，可以提供對數(shù)據(jù)結(jié)構(gòu)的進(jìn)一步了解。

下面的示例創(chuàng)建了沒有任何時(shí)間結(jié)構(gòu)的觀測值的直方圖和密度圖。

from pandas import Series
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
pyplot.figure(1)
pyplot.subplot(211)
series.hist()
pyplot.subplot(212)
series.plot(kind='kde')
pyplot.show()

運(yùn)行示例并查看繪圖。

圖中展示的一些分析結(jié)果包括：

分布不是高斯分布。
分布是左移的，可以是指數(shù)或雙高斯分布

每月波士頓搶劫密案數(shù)量密度圖

我們可以看到在建模之前進(jìn)行一些數(shù)據(jù)的維度變換的一些益處。

5.4 盒形-虛線圖

我們可以按年份對月度數(shù)據(jù)進(jìn)行分組，并了解每年觀察值的分布情況，以及觀察值可能如何變化。

我們確實(shí)希望從數(shù)據(jù)中看到一些趨勢（增加平均值或中位數(shù)），但看到其余的分布是如何變化的可能會(huì)更有趣。

下面的示例將觀察值按照年份劃分，并為每個(gè)觀察年創(chuàng)建一個(gè)盒型-虛線圖。最后一年（1974年）只包含10個(gè)月，可能與其他有12個(gè)月的觀察值的年份相比會(huì)無效。因此，只有1966年和1973年之間的數(shù)據(jù)被繪制。

from pandas import Series
from pandas import DataFrame
from pandas import TimeGrouper
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
groups = series['1966':'1973'].groupby(TimeGrouper('A'))
years = DataFrame()
for name, group in groups:
years[name.year] = group.values
years.boxplot()
pyplot.show()

圖中的一些觀察結(jié)果包括:

每年的中間值（紅線）顯示可能不是線性的趨勢
散列值和中間50％的數(shù)據(jù)（藍(lán)色框）不同，但時(shí)間長了可能會(huì)有變化。
較早的年份，也許是前2年，與數(shù)據(jù)集的其余部分有很大的不同。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

波士頓每月?lián)尳侔笖?shù)目的盒形-虛線圖

觀察結(jié)果表明，年度波動(dòng)可能不是系統(tǒng)性的，難以模擬。他們還建議將建模的前兩年的數(shù)據(jù)從模型分離出來可能有一些好處，如果確實(shí)不同。

這種對數(shù)據(jù)進(jìn)行年度觀察是一個(gè)有趣的方法，而且可以通過查看年度匯總統(tǒng)計(jì)數(shù)據(jù)和其中的變化，來進(jìn)一步跟進(jìn)未來的變化。

接下來，我們可以開始查看該序列的預(yù)測模型。

6. ARIMA 模型

在本節(jié)中，我們將會(huì)針對該問題開發(fā)自回歸積分滑動(dòng)平均模型，即ARIMA模型。

我們將其分四個(gè)步驟：

1.開發(fā)一套手動(dòng)配置的ARIMA模型。

2.使用ARIMA的網(wǎng)格搜索來找到優(yōu)化的模型。

3.預(yù)測殘差的分析，以評估模型中的所有偏差。

4.使用能量變換探索模型的改進(jìn)。

6.1 手動(dòng)配置ARIMA

Nasonasonal ARIMA（p，d，q）需要三個(gè)參數(shù)，并且通常是傳統(tǒng)地進(jìn)行手動(dòng)配置。

時(shí)間序列數(shù)據(jù)的分析。假設(shè)我們使用的是靜態(tài)時(shí)間序列。

時(shí)間序列基本上肯定是非靜態(tài)的。我們可以通過首先區(qū)分序列，并使用統(tǒng)計(jì)測試，來確保結(jié)果是靜態(tài)的。

下面的示例創(chuàng)建一個(gè)靜態(tài)版本的序列，并將其保存到文件stationary.csv。

from pandas import Series
from statsmodels.tsa.stattools import adfuller

# create a differe
def difference(dataset):
diff = list()
for i in range(1, len(dataset)):
value = dataset[i] - dataset[i - 1]
diff.append(value)
return Series(diff)

series = Series.from_csv('dataset.csv')
X = series.values
# difference data
stationary = difference(X)
stationary.index = series.index[1:]
# check if stationary
result = adfuller(stationary)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
# save
stationary.to_csv('stationary.csv')

運(yùn)行示例輸出統(tǒng)計(jì)顯著性測試的結(jié)果，以表明序列是否靜止。具體來說，強(qiáng)化迪基-福勒檢驗(yàn)。

結(jié)果表明，測試的統(tǒng)計(jì)值-3.980946小于-2.893的5％時(shí)的臨界值。這表明我們可以否決具有小于5％的顯著性水平的零假設(shè)（即，結(jié)果為統(tǒng)計(jì)結(jié)果的概率很低）。

拒絕零假設(shè)意味著過程沒有單位根，反過來說就是，該時(shí)間序列是固定的或者不具有時(shí)間依賴性結(jié)構(gòu)。

ADF Statistic: -3.980946
p-value: 0.001514
Critical Values:
5%: -2.893
1%: -3.503
10%: -2.584

這表明至少需要一個(gè)級別的差分。我們的ARIMA模型中的d參數(shù)應(yīng)至少為1的值。

下一步是分別選擇自回歸（AR）和移動(dòng)平均（MA）參數(shù)的滯后值，p和q。

我們可以通過審查自相關(guān)函數(shù)（ACF）和部分自相關(guān)函數(shù)（PACF）圖來做到這一點(diǎn)。

以下示例為該序列創(chuàng)建ACF和PACF圖。

from pandas import Series
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
pyplot.figure()
pyplot.subplot(211)
plot_acf(series, ax=pyplot.gca())
pyplot.subplot(212)
plot_pacf(series, ax=pyplot.gca())
pyplot.show()

運(yùn)行示例并查看繪圖，詳細(xì)了解一下如何設(shè)置ARIMA模型的p和q變量。

下面是從圖中得到的一些觀察結(jié)論。

ACF顯示1個(gè)月的顯著滯后。
PACF顯示了大約2個(gè)月的顯著滯后，滯后至大約12個(gè)月以后。
ACF和PACF在同一點(diǎn)顯示出下降，可能是AR和MA參數(shù)的混合。

p和q值的正確的起點(diǎn)是1或2。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

每月波士頓搶劫案數(shù)量的ACF和PACF圖

這個(gè)快速分析結(jié)果表明ARIMA（1,1,2）對原始數(shù)據(jù)可能是一個(gè)很好的切入點(diǎn)。

雷鋒網(wǎng)了解，實(shí)驗(yàn)表明，這套ARIMA的配置不會(huì)收斂，并引起底層庫的錯(cuò)誤。一些實(shí)驗(yàn)表明，該模型表現(xiàn)不太穩(wěn)定，同時(shí)定義了非零的AR和MA階數(shù)。

該模型可以簡化為ARIMA（0,1,2）。下面的示例演示了此ARIMA模型在測試工具上的性能。

from pandas import Series
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
from math import sqrt
# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# predict
model = ARIMA(history, order=(0,1,2))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

運(yùn)行此示例結(jié)果的RMSE值為49.821，低于持續(xù)性模型（Persistence）。

...
>Predicted=280.614, Expected=287
>Predicted=302.079, Expected=355
>Predicted=340.210, Expected=460
>Predicted=405.172, Expected=364
>Predicted=333.755, Expected=487
RMSE: 49.821

這是一個(gè)好的開始，但是我們可以達(dá)到改進(jìn)結(jié)果通過一個(gè)配置更好的ARIMA模型。

6.2 網(wǎng)格搜索ARIMA超參數(shù)

許多ARIMA配置在此數(shù)據(jù)集上不穩(wěn)定，但可能有其他超參數(shù)導(dǎo)向一個(gè)表現(xiàn)良好的模型。

在本節(jié)中，我們將搜索p，d和q的值，以找出不會(huì)導(dǎo)致錯(cuò)誤的組合，并找到可獲得最佳性能的組合。我們將使用網(wǎng)格搜索來探索整數(shù)值子集中的所有組合。

具體來說，我們將搜索以下參數(shù)的所有組合：

p: 0 - 12.
d: 0 - 3.
q: 0 -12.

這是（13 * 4 * 13），或676，測試工具的運(yùn)行結(jié)果，并且將花費(fèi)需要一些時(shí)間用來執(zhí)行。

下面是利用網(wǎng)格搜索的測試工具的完整示例：

import warnings
from pandas import Series
from statsmodels.tsa.arima_model import ARIMA
from sklearn.metrics import mean_squared_error
from math import sqrt

# evaluate an ARIMA model for a given order (p,d,q) and return RMSE
def evaluate_arima_model(X, arima_order):
# prepare training dataset
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
history = [x for x in train]
# make predictions
predictions = list()
for t in range(len(test)):
model = ARIMA(history, order=arima_order)
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
predictions.append(yhat)
history.append(test[t])
# calculate out of sample error
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
return rmse

# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(dataset, p_values, d_values, q_values):
dataset = dataset.astype('float32')
best_score, best_cfg = float("inf"), None
for p in p_values:
for d in d_values:
for q in q_values:
order = (p,d,q)
try:
mse = evaluate_arima_model(dataset, order)
if mse < best_score:
best_score, best_cfg = mse, order
print('ARIMA%s MSE=%.3f' % (order,mse))
except:
continue
print('Best ARIMA%s MSE=%.3f' % (best_cfg, best_score))

# load dataset
series = Series.from_csv('dataset.csv')
# evaluate parameters
p_values = range(0,13)
d_values = range(0, 4)
q_values = range(0, 13)
warnings.filterwarnings("ignore")
evaluate_models(series.values, p_values, d_values, q_values)

運(yùn)行該示例并運(yùn)行所有組合，并將無錯(cuò)誤收斂的結(jié)果報(bào)出來。

該示例當(dāng)前硬件上運(yùn)行至少需要2小時(shí)。

結(jié)果表明，最佳配置為ARIMA（0，1，2），巧合的是，這在前面的部分中已經(jīng)被證明。

...
ARIMA(6, 1, 0) MSE=52.437
ARIMA(6, 2, 0) MSE=58.307
ARIMA(7, 1, 0) MSE=51.104
ARIMA(7, 1, 1) MSE=52.340
ARIMA(8, 1, 0) MSE=51.759
Best ARIMA(0, 1, 2) MSE=49.821

6.3 檢查殘差

對一個(gè)模型來說，一個(gè)不錯(cuò)的最終檢查，是檢查其殘差預(yù)測誤差。

理想情況下，殘差的分布應(yīng)該是符合零均值的高斯分布。

我們可以通過用直方圖和密度圖繪制殘差來檢查這一點(diǎn)。

以下示例為計(jì)算測試集上預(yù)測的殘差，并創(chuàng)建該密度圖。

from pandas import Series
from pandas import DataFrame
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
from math import sqrt
from matplotlib import pyplot
# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
# walk-foward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# predict
model = ARIMA(history, order=(0,1,2))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
# errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = DataFrame(residuals)
pyplot.figure()
pyplot.subplot(211)
residuals.hist(ax=pyplot.gca())
pyplot.subplot(212)
residuals.plot(kind='kde', ax=pyplot.gca())
pyplot.show()

運(yùn)行該示例會(huì)創(chuàng)建兩個(gè)圖表。

這些圖表征了一個(gè)有長右尾的高斯樣分布。

這可能表明該預(yù)測是有偏差的，并且在這種情況下，在建模之前對原始數(shù)據(jù)進(jìn)行基于能量的變換可能是有用的。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

殘差密度圖

通過檢查所有時(shí)間序列的殘差中是否有某種類型的自相關(guān)性也是一個(gè)好辦法。如果存在，它將表明模型有更多的機(jī)會(huì)來建模數(shù)據(jù)中的時(shí)間結(jié)構(gòu)。

下面的示例是重新計(jì)算殘差，并創(chuàng)建ACF和PACF圖以檢查出比較明顯的的自相關(guān)性。

from pandas import Series
from pandas import DataFrame
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from math import sqrt
from matplotlib import pyplot
# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
# walk-foward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# predict
model = ARIMA(history, order=(0,1,2))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
# errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = DataFrame(residuals)
pyplot.figure()
pyplot.subplot(211)
plot_acf(residuals, ax=pyplot.gca())
pyplot.subplot(212)
plot_pacf(residuals, ax=pyplot.gca())
pyplot.show()

結(jié)果表明，該時(shí)間序列中存在的少量的自相關(guān)性已被模型捕獲。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

殘差A(yù)CF和PACF圖

6.4 Box-Cox 變換數(shù)據(jù)集

Box-Cox變換是可以對一組能量變換進(jìn)行評估的方法，包括但不限于數(shù)據(jù)的對數(shù)，平方根和互逆變換。

下面的示例對數(shù)據(jù)的執(zhí)行一個(gè)對數(shù)變換，并生成一些曲線以查看對時(shí)間序列的影響。

from pandas import Series
from pandas import DataFrame
from scipy.stats import boxcox
from matplotlib import pyplot
from statsmodels.graphics.gofplots import qqplot
series = Series.from_csv('dataset.csv')
X = series.values
transformed, lam = boxcox(X)
print('Lambda: %f' % lam)
pyplot.figure(1)
# line plot
pyplot.subplot(311)
pyplot.plot(transformed)
# histogram
pyplot.subplot(312)
pyplot.hist(transformed)
# q-q plot
pyplot.subplot(313)
qqplot(transformed, line='r', ax=pyplot.gca())
pyplot.show()

運(yùn)行示例創(chuàng)建三個(gè)圖：經(jīng)變換的時(shí)間序列的線圖，展示變換后值的分布的直方圖，以及展示變換后的值得分布與理想化高斯分布相比的Q-Q圖。

這些圖中的一些觀察結(jié)果如下：

大的波動(dòng)已在時(shí)間序列的線圖中被移除。
直方圖顯示出了數(shù)值更平坦更均勻（表現(xiàn)良好）的分布。
Q-Q圖是合理的，但仍然不適合高斯分布。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

每月波士頓搶劫案數(shù)據(jù)的BoxCox變換圖

毫無疑問，Box-Cox變換對時(shí)間序列做了一些事情，而且是有效的。

在使用轉(zhuǎn)換數(shù)據(jù)測試ARIMA模型之前，我們必須有一種方法來進(jìn)行逆轉(zhuǎn)轉(zhuǎn)換，以便可以將對在轉(zhuǎn)換后的尺度上進(jìn)行訓(xùn)練的模型所做的預(yù)測轉(zhuǎn)換回原始尺度。

示例中使用的boxcox（）函數(shù)通過優(yōu)化損失函數(shù)找到理想的lambda值。

lambda在以下函數(shù)中用于數(shù)據(jù)轉(zhuǎn)換：

transform = log(x), if lambda == 0
transform = (x^lambda - 1) / lambda, if lambda != 0

這個(gè)變換函數(shù)可以直接逆過來，如下：

x = exp(transform) if lambda == 0
x = exp(log(lambda * transform + 1) / lambda)

這個(gè)逆Box-Cox變換函數(shù)可以在Python中如下實(shí)現(xiàn)：

# invert box-cox transform
from math import log
from math import exp
def boxcox_inverse(value, lam):
if lam == 0:
return exp(value)
return exp(log(lam * value + 1) / lam)

我們可以用Box-Cox變換重新評估ARIMA（0,1,2）模型。

這包括在擬合ARIMA模型之前的初次變換，然后在存儲之前對預(yù)測進(jìn)行反變換，以便稍后與期望值進(jìn)行比較。

boxcox（）函數(shù)可能失敗。在實(shí)踐中，我已經(jīng)認(rèn)識到這一點(diǎn)，它似乎是返回了一個(gè)小于-5的lambda值。按照慣例，lambda值在-5和5之間計(jì)算。

對于小于-5的lambda值添加檢查，如果是這種情況，則假定lambda值為1，并且使用原始數(shù)據(jù)來擬合模型。 lambda值為1與未轉(zhuǎn)換相同，因此逆變換沒有效果。

下面列出了完整的示例。

from pandas import Series
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
from math import sqrt
from math import log
from math import exp
from scipy.stats import boxcox

# invert box-cox transform
def boxcox_inverse(value, lam):
if lam == 0:
return exp(value)
return exp(log(lam * value + 1) / lam)

# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]
# walk-foward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
# transform
transformed, lam = boxcox(history)
if lam < -5:
transformed, lam = history, 1
# predict
model = ARIMA(transformed, order=(0,1,2))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
# invert transformed prediction
yhat = boxcox_inverse(yhat, lam)
predictions.append(yhat)
# observation
obs = test[i]
history.append(obs)
print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

運(yùn)行此示例將打印每次迭代的預(yù)測值和預(yù)期值

注意，當(dāng)使用boxcox（）轉(zhuǎn)換函數(shù)時(shí)，可能會(huì)看到一些警告; 例如：

RuntimeWarning: overflow encountered in square
llf -= N / 2.0 * np.log(np.sum((y - y_mean)**2. / N, axis=0))

不過現(xiàn)在可以忽略這些。

在轉(zhuǎn)換數(shù)據(jù)上生成的模型最終的RMSE為49.443。對于未轉(zhuǎn)換的數(shù)據(jù)，這是一個(gè)比ARIMA模型更小的誤差，但只是很細(xì)微的，它可能是，也可能不是統(tǒng)計(jì)意義上的不同。

...
>Predicted=276.253, Expected=287
>Predicted=299.811, Expected=355
>Predicted=338.997, Expected=460
>Predicted=404.509, Expected=364
>Predicted=349.336, Expected=487
RMSE: 49.443

我們將使用該Box-Cox變換的模型作為最終模型。

7. 模型驗(yàn)證

在開發(fā)完模型并選擇最終模型后，必須對其進(jìn)行驗(yàn)證和確定。

驗(yàn)證是一個(gè)可選的過程，為我們提供了“最后的檢查手段”，以確保我們沒有被自己欺騙。

此部分包括以下步驟：

1. 確定模型：訓(xùn)練并保存最終模型。

2. 進(jìn)行預(yù)測: 加載最終模型并進(jìn)行預(yù)測。

3. 驗(yàn)證模型: 加載和驗(yàn)證最終模型。

7.1 確定模型

模型的確定涉及在整個(gè)數(shù)據(jù)集上擬合ARIMA模型，值得注意的一點(diǎn)是，該步驟中的數(shù)據(jù)集已經(jīng)經(jīng)過了變換。

一旦擬合好之后，模型可以保存到文件以供后續(xù)使用。這里由于我們對數(shù)據(jù)執(zhí)行了Box-Cox變換，因此我們需要明確所選擇的lamda值，以便使得來自模型的任何預(yù)測都可以被轉(zhuǎn)換回原始的未轉(zhuǎn)換的尺度。

下面的示例展示了在Box-Cox變換數(shù)據(jù)集上使用ARIMA（0,1,2）模型，并將整個(gè)擬合對象和lambda值保存到文件中的過程。

from pandas import Series
from statsmodels.tsa.arima_model import ARIMA
from scipy.stats import boxcox
import numpy
# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
# transform data
transformed, lam = boxcox(X)
# fit model
model = ARIMA(transformed, order=(0,1,2))
model_fit = model.fit(disp=0)
# save model
model_fit.save('model.pkl')
numpy.save('model_lambda.npy', [lam])

運(yùn)行該示例會(huì)創(chuàng)建兩個(gè)本地文件：

model.pkl這是通過ARIMA.fit() 調(diào)用生成的ARIMAResult對象。包括了各種系數(shù)和在擬合模型時(shí)返回的所有其他內(nèi)部數(shù)據(jù)。
model_lambda.npy這是作為一行一列NumPy數(shù)組對象存儲的lambda值。

這些文件或許會(huì)保存的過于詳細(xì)。其實(shí)在操作中真正需要的是：來自模型的AR和MA系數(shù)、用于差異數(shù)量的d參數(shù)、以及滯后的觀察值、模型殘差和用于變換的lamda值。

7.2 進(jìn)行預(yù)測

這里，最直接的方法是加載模型并進(jìn)行單一預(yù)測。

這是相對直接的，涉及恢復(fù)保存的模型、lambda值和調(diào)用forecast() 方法。

然而，在當(dāng)前穩(wěn)定版本的statsmodels庫（v0.6.1）中存在一個(gè)Bug，它的錯(cuò)誤提示如下：

TypeError: __new__() takes at least 3 arguments (1 given)

據(jù)我測試發(fā)現(xiàn)，這個(gè)bug也似乎存在于statsmodels的0.8版本候選1版本中。

更多詳細(xì)信息，請參閱Zae Myung Kim的討論和在GitHub上的解決辦法。

鏈接：https://github.com/statsmodels/statsmodels/pull/3217

我們可以使用一個(gè)猴子補(bǔ)?。╩onkey patch）來解決這個(gè)問題，在保存之前將一個(gè) __getnewargs __() 實(shí)例函數(shù)添加到ARIMA類中。

下面的示例與上一小節(jié)的變化相同。

from pandas import Series
from statsmodels.tsa.arima_model import ARIMA
from scipy.stats import boxcox
import numpy

# monkey patch around bug in ARIMA class
def __getnewargs__(self):
return ((self.endog),(self.k_lags, self.k_diff, self.k_ma))

ARIMA.__getnewargs__ = __getnewargs__

# load data
series = Series.from_csv('dataset.csv')
# prepare data
X = series.values
X = X.astype('float32')
# transform data
transformed, lam = boxcox(X)
# fit model
model = ARIMA(transformed, order=(0,1,2))
model_fit = model.fit(disp=0)
# save model
model_fit.save('model.pkl')
numpy.save('model_lambda.npy', [lam])

我們現(xiàn)在可以加載模型并進(jìn)行單一預(yù)測。

下面這個(gè)示例展示了加載模型，對下一時(shí)間步進(jìn)進(jìn)行了預(yù)測，反轉(zhuǎn)Box-Cox變換，并打印預(yù)測結(jié)果的過程。

from pandas import Series
from statsmodels.tsa.arima_model import ARIMAResults
from math import exp
from math import log
import numpy
# invert box-cox transform
def boxcox_inverse(value, lam):
if lam == 0:
return exp(value)
return exp(log(lam * value + 1) / lam)
model_fit = ARIMAResults.load('model.pkl')
lam = numpy.load('model_lambda.npy')
yhat = model_fit.forecast()[0]
yhat = boxcox_inverse(yhat, lam)
print('Predicted: %.3f' % yhat)

最后得到的預(yù)測值為452。

如果我們在validation.csv中查看的話，我們可以看到下一個(gè)時(shí)間段的第一行的數(shù)值是452。模型達(dá)到了100% 的正確率，這是非常驚人的（也是幸運(yùn)的）。

Predicted: 452.043

7.3 驗(yàn)證模型

在這一步驟中，我們可以加載模型并以模擬操作的方式使用它。

在測試工具部分中，我們將原始數(shù)據(jù)集的最后12個(gè)月保存在單獨(dú)的文件中，以驗(yàn)證最終模型。

我們現(xiàn)在可以加載此validation.csv文件，并用它來檢驗(yàn)我們的模型對真正“未知”的數(shù)據(jù)究竟表現(xiàn)如何。

這里可以采用兩種方式：

加載模型并使用它預(yù)測未來的12個(gè)月。超過前一兩個(gè)月的預(yù)測將很快開始降低技能（start to degrade in skill）。
加載模型并以步進(jìn)預(yù)測的方式使用該模型，更新每個(gè)時(shí)間步長上的變換和模型。這是首選的方法，在實(shí)踐中使用這個(gè)模型，它將表現(xiàn)出最佳的性能。

與前面章節(jié)中的模型評估一樣，我們將以步進(jìn)預(yù)測的方式進(jìn)行預(yù)測。這意味著我們將在驗(yàn)證數(shù)據(jù)集中逐步引導(dǎo)時(shí)間點(diǎn)，并將觀察結(jié)果作為歷史記錄的更新。

from pandas import Series
from matplotlib import pyplot
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.arima_model import ARIMAResults
from scipy.stats import boxcox
from sklearn.metrics import mean_squared_error
from math import sqrt
from math import exp
from math import log
import numpy

# invert box-cox transform
def boxcox_inverse(value, lam):
if lam == 0:
return exp(value)
return exp(log(lam * value + 1) / lam)

# load and prepare datasets
dataset = Series.from_csv('dataset.csv')
X = dataset.values.astype('float32')
history = [x for x in X]
validation = Series.from_csv('validation.csv')
y = validation.values.astype('float32')
# load model
model_fit = ARIMAResults.load('model.pkl')
lam = numpy.load('model_lambda.npy')
# make first prediction
predictions = list()
yhat = model_fit.forecast()[0]
yhat = boxcox_inverse(yhat, lam)
predictions.append(yhat)
history.append(y[0])
print('>Predicted=%.3f, Expected=%3.f' % (yhat, y[0]))
# rolling forecasts
for i in range(1, len(y)):
# transform
transformed, lam = boxcox(history)
if lam < -5:
transformed, lam = history, 1
# predict
model = ARIMA(transformed, order=(0,1,2))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
# invert transformed prediction
yhat = boxcox_inverse(yhat, lam)
predictions.append(yhat)
# observation
obs = y[i]
history.append(obs)
print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(y, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)
pyplot.plot(y)
pyplot.plot(predictions, color='red')
pyplot.show()

運(yùn)行以上示例將打印出在驗(yàn)證數(shù)據(jù)集中所有時(shí)間步長上的每個(gè)預(yù)測值和預(yù)期值。

修正階段最終的RMSE預(yù)測結(jié)果為53次搶劫。這與預(yù)期的49次并沒有太大的出入，我希望它在簡單的持久性模型上也能有類似的表現(xiàn)。

>Predicted=452.043, Expected=452
>Predicted=423.088, Expected=391
>Predicted=408.378, Expected=500
>Predicted=482.454, Expected=451
>Predicted=445.944, Expected=375
>Predicted=413.881, Expected=372
>Predicted=413.209, Expected=302
>Predicted=355.159, Expected=316
>Predicted=363.515, Expected=398
>Predicted=406.365, Expected=394
>Predicted=394.186, Expected=431
>Predicted=428.174, Expected=431
RMSE: 53.078

我們在這里同時(shí)提供了預(yù)測結(jié)果與驗(yàn)證數(shù)據(jù)集之間相互比較的圖表。

該預(yù)測具有持續(xù)性預(yù)測的特性。這表明雖然這個(gè)時(shí)間序列有明顯的趨勢，但它仍然是一個(gè)相當(dāng)困難的問題。

時(shí)間序列預(yù)測教程：如何利用 Python 預(yù)測波士頓每月持械搶劫案數(shù)量？

驗(yàn)證時(shí)間步驟的預(yù)測

教程擴(kuò)展

需要說明的是，本教程并不完美，你可以有很多改善結(jié)果的方法。

本節(jié)列出了一些或許可行的改進(jìn)措施。

統(tǒng)計(jì)顯著性檢驗(yàn)。使用統(tǒng)計(jì)學(xué)測試來檢查不同模型之間的結(jié)果差異是否具有統(tǒng)計(jì)學(xué)的意義。T檢驗(yàn)將是一個(gè)很好的入手點(diǎn)。
使用數(shù)據(jù)變換的網(wǎng)格進(jìn)行搜索。通過Box-Cox變換在ARIMA超參數(shù)中重復(fù)網(wǎng)格搜索過程，并查看是否可以實(shí)現(xiàn)一個(gè)與此前不同的，并且更好的參數(shù)集。
檢查殘差。通過Box-Cox變換檢測最終模型的預(yù)測誤差的殘差，以查看是否存在可以進(jìn)一步優(yōu)化的偏差或自相關(guān)。
精簡模型存儲。嘗試簡化模型保存，僅僅存儲那些必要的系數(shù)，而不是整個(gè)ARIMAResults對象。
手動(dòng)處理趨勢。使用線性或非線性模型直接對趨勢進(jìn)行建模，并明確地從序列中刪除它。如果該趨勢是非線性的，并且更適合非線性建模，那這種方式可能會(huì)帶來更好的性能表現(xiàn)。
置信區(qū)間。顯示驗(yàn)證數(shù)據(jù)集中預(yù)測的置信區(qū)間。
數(shù)據(jù)選擇?？紤]在沒有前兩年數(shù)據(jù)的情況下對問題進(jìn)行建模，并查看這是否對預(yù)測技能（forecast skill）有影響。