0
編者按：本文是澳大利亞知名機器學習專家 Jason Brownlee 撰寫的教程，極其全面細致，一步步向讀者解釋如何操作，以及為什么這么做。雷鋒網整理編譯，特與大家分享。更多AI開發技術文章，請關注AI研習社（微信號：okweiwu）。
Jason Brownlee：時間序列預測法是一個過程，而獲得良好預測結果的唯一途徑是實踐這個過程。
在本教程中，您將了解如何利用Python語言來預測波士頓每月持械搶劫案發生的數量。
本教程所述為您提供了一套處理時間序列預測問題的框架，包括方法步驟和工具，通過實踐，可以用它來解決自己遇到的相關問題。
本教程結束之后，您將了解：
如何核查Python環境并準確地定義一個時間序列預測問題。
如何構建一套測試工具鏈，用于評估模型，開發預測原型。以及如何通過時間序列分析工具更好地理解你的問題。
如何開發自回歸積分滑動平均模型（ARIMA），將其保存到文件，并在之后加載它對新的時間步驟進行預測。
讓我們開始吧。
波士頓
在本教程中，我們將端到端地來解析一個時間序列預測工程，從下載數據集、定義問題到訓練出最終模型并進行預測。
該工程并不面面俱到，但展示了如何通過系統性地處理時間序列預測問題，來快速獲得好的結果。
我們將通過如下步驟來解析該工程 ：
環境配置.
問題描述
測試工具鏈
持續性.
數據分析
ARIMA 模型
模型驗證
該部分提供一個解決時間序列預測問題的模板，你可以套用自己的數據集。
本教程假設在已安裝好了SciPy及其相關依賴庫的環境下進行，包括：
SciPy
NumPy
Matplotlib
Pandas
scikit-learn
statsmodels
我使用的是Python2.7。該腳本將幫助你確認你安裝這些庫的版本。
# 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__))
我編寫該教程的所用的開發環境顯示的結果如下：
scipy: 0.18.1
numpy: 1.11.2
matplotlib: 1.5.3
pandas: 0.19.1
sklearn: 0.18.1
statsmodels: 0.6.1
我們研究的問題是預測美國波士頓每月持械搶劫案的數量。
該數據集提供了從1966年1月到1975年10月在波士頓每月發生的的武裝搶劫案的數量， 或者說，是這十年間的數據。
這些值是一些記錄的觀察計數，總有118個觀察值。
數據集來自于McCleary＆Hay（1980）。
您可以了解有關此數據集的更多信息，并直接從 DataMarket 下載。
下載地址：https://datamarket.com/data/set/22ob/monthly-boston-armed-robberies-jan1966-oct1975-deutsch-and-alt-1977#!ds=22ob&display=line
將數據集下載為CSV文件，并將其放在當前工作目錄中，文件名為“robberies.csv”。
我們必須開發一套測試工具鏈來審查數據和評估候選模型。
這包含兩個步驟：
1.定義好驗證數據集。
2.開發一套模型評估方法。
這個數據集不是當前時間段的數據集。 這意味著我們不能輕輕松松收集更新后的數據來驗證該模型。
因此，我們假設現在是1974年10月，并將分析和模型選擇的數據集里扣除最后一年的數據。
最后一年的數據將用于驗證生成的最終模型。
下面的代碼會加載該數據集為Pandas序列對象并將其切分成兩部分，一個用于模型開發（dataset.csv）和另一個用于驗證（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')
運行該示例將會創建兩個文件并打印出每個文件中的觀察數。
1. Dataset 106, Validation 12
這些文件的具體內容是：
dataset.csv：1966年1月至1974年10月的觀察（106次觀察）
validation.csv：1974年11月至1975年10月的觀察（12次觀察）
驗證數據集大小是原始數據集的10％。
請注意，由于保存的數據集中沒有標題行，因此，稍后處理這些文件時，我們也不需要滿足此要求。
模型評估將僅對上一節中準備的dataset.csv中的數據進行處理。
模型評估涉及兩個要素：
1. 評價指標。
2. 測試策略。
這些觀察值是搶劫案的計數。
我們將使用均方根誤差（RMSE）的方式來評價預測的好壞。這將給予錯誤的預測更多的權重，并且將具有與原始數據相同的單位。
在RMSE進行計算和報告之前，必須反轉對數據的所有轉換，以使不同方法的效果可直接比較。
我們可以使用scikit-learn庫的輔助函數mean_squared_error（）來執行RMSE計算，該函數計算出了預期值列表（測試集）和預測列表之間的均方誤差。 然后我們可以取這個值的平方根來作為一個RMSE分數。
例如：
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)
將使用步進驗證的方式(walk-forward)來評估候選模型。
這是因為該問題定義需要滾動式預測的模型。給定所有可用數據，這需要逐步進行預測。
步進驗證的方式將按照如下幾步執行：
數據集的前50％將被用來訓練模型。
剩余的50％的數據集將被迭代并測試模型。
對于測試數據集中的每個步驟：
訓練該模型
利用該模型預測一次，并經預測結果保存下來用于之后的驗證
測試數據集的數據作為下一個迭代的訓練集數據
對在測試數據集迭代期間進行的預測結果進行評估，并報告RMSE得分。
由于較小的數據規模，我們將允許在每次預測之前，在給定所有可用數據的情況下重新訓練模型。
我們可以使用簡單的NumPy和Python來編寫測試工具的代碼。
首先，我們可以將數據集直接分成訓練和測試集。 我們總要將加載的數據轉換為float32類型，以防止仍然有一些數據為String或Integer數據類型。
# 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:]
接下來，我們可以時間步長遍歷測試數據集。訓練數據集存儲在一個Python的list對象中，因為我們需要在每次迭代中很容易的附加一些觀察值，Numpy數組連接有些overkill。
模型進行的預測，按慣例稱為yhat。這是由于結果或觀察被稱為y，而yhat（有上標“^”的y）是用于預測y變量的數學符號。
打印預測和觀察結果，對每個觀察值進行完整性檢查，以防模型存在問題。
# 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))
在數據分析和建模陷入困境之前的第一步是建立一個評估原型。
這里將為利用測試工具進行模型評估和性能測量，分別提供一套模版。通過該性能測量，可以將當前模型和其他更復雜的預測模型進行比較。
時間序列預測的預測原型被稱為天真預測或Persistence。
在這里，先前時間步驟的觀察結果將被用作于下一時間步長的預測。
我們可以將其直接插入上一節中定義的測試工具中。
下面羅列了完整的代碼：
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)
運行測試工具，打印測試數據集中每次迭代的預測值和實際觀察值。
示例到進行打印模型的RMSE這一步結束。
在該示例中，我們可以看到持久性模型實現了51.844的RMSE。 這意味著，對于每個預測，該模型中平均有51次搶劫值預測是錯誤的。
...
>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
我們現在有了一個預測方法和性能評估的原型; 現在我們可以開始深入數據。
我們可以使用統計匯總和數據的繪圖，來快速了解該預測問題的數據結構。
在本節中，我們將從四個角度來看數據：
1.摘要統計。
2.線圖。
3.密度圖。
4.盒型-虛線圖。
利用文本編輯器中打開數據dataset.csv文件或原始的robberies.csv文件，并查看數據。
快速檢查表明，數據集里沒有明顯丟失的觀察數據。 我們可能已經更早的發現了像NaN或'？'等值數，如果我們試圖強制序列化數據為浮點值。
摘要統計提供了觀察值的快速預覽。 它可以幫助我們快速了解我們正在使用的東西。
以下示例為計算并打印時間系列的摘要統計信息。
from pandas import Series
series = Series.from_csv('dataset.csv')
print(series.describe())
運行該示例會提供大量的摘要信息供回顧。
這些統計數據包含的一些信息為：
觀察次數（計數）與我們的期望相符，這意味著我們將數據正確處理了。
平均值約為173，我們可以想象出我們在這個序列中的級別。
在112次搶劫時，標準偏差（平均值的平均分布）相對較大。
百分位數與標準差一致表明數據有很大的差異。
count 106.000000
mean 173.103774
std 112.231133
min 29.000000
25% 74.750000
50% 144.500000
75% 271.750000
max 487.000000
對于序列中比較大的變化將可能模型很難進行高度準確的預測，如果它是由隨機波動（例如，非系統的）引起，
畫一個時間序列的線圖可以對該問題提供很多深入的信息。
下面的示例創建并顯示出數據集的線圖。
from pandas import Series
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
series.plot()
pyplot.show()
運行示例并查看繪圖。 注意系列中的任何明顯的時間結構。
圖中的一些觀察結果包括：
隨著時間的推移，劫案有增加的趨勢。
似乎沒有任何明顯的異常值。
每年都有相對較大的波動，上下波動。
后幾年的波動似乎大于前幾年的波動。
整體趨勢表明該數據集幾乎肯定不是穩定的，并且波形的明顯變化也可能對此有所提示。
每月波士頓搶劫案數據線圖
這些簡單的觀察結果表明，我們可以看到對趨勢進行建模和從時間序列中將它提取出來的好處。
或者，我們可以使用差分法從而使序列的建模平穩。如果后期波動有增長趨勢，我們甚至可能需要兩個差分到不同級別。
提取出密度圖，可以提供對數據結構的進一步了解。
下面的示例創建了沒有任何時間結構的觀測值的直方圖和密度圖。
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()
運行示例并查看繪圖。
圖中展示的一些分析結果包括：
分布不是高斯分布。
分布是左移的，可以是指數或雙高斯分布
每月波士頓搶劫密案數量密度圖
我們可以看到在建模之前進行一些數據的維度變換的一些益處。
我們可以按年份對月度數據進行分組，并了解每年觀察值的分布情況，以及觀察值可能如何變化。
我們確實希望從數據中看到一些趨勢（增加平均值或中位數），但看到其余的分布是如何變化的可能會更有趣。
下面的示例將觀察值按照年份劃分，并為每個觀察年創建一個盒型-虛線圖。 最后一年（1974年）只包含10個月，可能與其他有12個月的觀察值的年份相比會無效。 因此，只有1966年和1973年之間的數據被繪制。
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()
圖中的一些觀察結果包括:
每年的中間值（紅線）顯示可能不是線性的趨勢
散列值和中間50％的數據（藍色框）不同，但時間長了可能會有變化。
較早的年份，也許是前2年，與數據集的其余部分有很大的不同。
波士頓每月搶劫案數目的盒形-虛線圖
觀察結果表明，年度波動可能不是系統性的，難以模擬。 他們還建議將建模的前兩年的數據從模型分離出來可能有一些好處，如果確實不同。
這種對數據進行年度觀察是一個有趣的方法，而且可以通過查看年度匯總統計數據和其中的變化，來進一步跟進未來的變化。
接下來，我們可以開始查看該序列的預測模型。
在本節中，我們將會針對該問題開發自回歸積分滑動平均模型，即ARIMA模型。
我們將其分四個步驟：
1.開發一套手動配置的ARIMA模型。
2.使用ARIMA的網格搜索來找到優化的模型。
3.預測殘差的分析，以評估模型中的所有偏差。
4.使用能量變換探索模型的改進。
Nasonasonal ARIMA（p，d，q）需要三個參數，并且通常是傳統地進行手動配置。
時間序列數據的分析。假設我們使用的是靜態時間序列。
時間序列基本上肯定是非靜態的。我們可以通過首先區分序列，并使用統計測試，來確保結果是靜態的。
下面的示例創建一個靜態版本的序列，并將其保存到文件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')
運行示例輸出統計顯著性測試的結果，以表明序列是否靜止。 具體來說，強化迪基-福勒檢驗。
結果表明，測試的統計值-3.980946小于-2.893的5％時的臨界值。 這表明我們可以否決具有小于5％的顯著性水平的零假設（即，結果為統計結果的概率很低）。
拒絕零假設意味著過程沒有單位根，反過來說就是，該時間序列是固定的或者不具有時間依賴性結構。
ADF Statistic: -3.980946
p-value: 0.001514
Critical Values:
5%: -2.893
1%: -3.503
10%: -2.584
這表明至少需要一個級別的差分。 我們的ARIMA模型中的d參數應至少為1的值。
下一步是分別選擇自回歸（AR）和移動平均（MA）參數的滯后值，p和q。
我們可以通過審查自相關函數（ACF）和部分自相關函數（PACF）圖來做到這一點。
以下示例為該序列創建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()
運行示例并查看繪圖，詳細了解一下如何設置ARIMA模型的p和q變量。
下面是從圖中得到的一些觀察結論。
ACF顯示1個月的顯著滯后。
PACF顯示了大約2個月的顯著滯后，滯后至大約12個月以后。
ACF和PACF在同一點顯示出下降，可能是AR和MA參數的混合。
p和q值的正確的起點是1或2。
每月波士頓搶劫案數量的ACF和PACF圖
這個快速分析結果表明ARIMA（1,1,2）對原始數據可能是一個很好的切入點。
雷鋒網了解，實驗表明，這套ARIMA的配置不會收斂，并引起底層庫的錯誤。一些實驗表明，該模型表現不太穩定，同時定義了非零的AR和MA階數。
該模型可以簡化為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)
運行此示例結果的RMSE值為49.821，低于持續性模型（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
這是一個好的開始，但是我們可以達到改進結果通過一個配置更好的ARIMA模型。
許多ARIMA配置在此數據集上不穩定，但可能有其他超參數導向一個表現良好的模型。
在本節中，我們將搜索p，d和q的值，以找出不會導致錯誤的組合，并找到可獲得最佳性能的組合。 我們將使用網格搜索來探索整數值子集中的所有組合。
具體來說，我們將搜索以下參數的所有組合：
p: 0 - 12.
d: 0 - 3.
q: 0 -12.
這是（13 * 4 * 13），或676，測試工具的運行結果，并且將花費需要一些時間用來執行。
下面是利用網格搜索的測試工具的完整示例：
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)
運行該示例并運行所有組合，并將無錯誤收斂的結果報出來。
該示例當前硬件上運行至少需要2小時。
結果表明，最佳配置為ARIMA（0，1，2），巧合的是，這在前面的部分中已經被證明。
...
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
對一個模型來說，一個不錯的最終檢查，是檢查其殘差預測誤差。
理想情況下，殘差的分布應該是符合零均值的高斯分布。
我們可以通過用直方圖和密度圖繪制殘差來檢查這一點。
以下示例為計算測試集上預測的殘差，并創建該密度圖。
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()
運行該示例會創建兩個圖表。
這些圖表征了一個有長右尾的高斯樣分布。
這可能表明該預測是有偏差的，并且在這種情況下，在建模之前對原始數據進行基于能量的變換可能是有用的。
殘差密度圖
通過檢查所有時間序列的殘差中是否有某種類型的自相關性也是一個好辦法。 如果存在，它將表明模型有更多的機會來建模數據中的時間結構。
下面的示例是重新計算殘差，并創建ACF和PACF圖以檢查出比較明顯的的自相關性。
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()
結果表明，該時間序列中存在的少量的自相關性已被模型捕獲。
殘差ACF和PACF圖
Box-Cox變換是可以對一組能量變換進行評估的方法，包括但不限于數據的對數，平方根和互逆變換。
下面的示例對數據的執行一個對數變換，并生成一些曲線以查看對時間序列的影響。
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()
運行示例創建三個圖：經變換的時間序列的線圖，展示變換后值的分布的直方圖，以及展示變換后的值得分布與理想化高斯分布相比的Q-Q圖。
這些圖中的一些觀察結果如下：
大的波動已在時間序列的線圖中被移除。
直方圖顯示出了數值更平坦更均勻（表現良好）的分布。
Q-Q圖是合理的，但仍然不適合高斯分布。
每月波士頓搶劫案數據的BoxCox變換圖
毫無疑問，Box-Cox變換對時間序列做了一些事情，而且是有效的。
在使用轉換數據測試ARIMA模型之前，我們必須有一種方法來進行逆轉轉換，以便可以將對在轉換后的尺度上進行訓練的模型所做的預測轉換回原始尺度。
示例中使用的boxcox（）函數通過優化損失函數找到理想的lambda值。
lambda在以下函數中用于數據轉換：
transform = log(x), if lambda == 0
transform = (x^lambda - 1) / lambda, if lambda != 0
這個變換函數可以直接逆過來，如下：
x = exp(transform) if lambda == 0
x = exp(log(lambda * transform + 1) / lambda)
這個逆Box-Cox變換函數可以在Python中如下實現：
# 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模型之前的初次變換，然后在存儲之前對預測進行反變換，以便稍后與期望值進行比較。
boxcox（）函數可能失敗。 在實踐中，我已經認識到這一點，它似乎是返回了一個小于-5的lambda值。 按照慣例，lambda值在-5和5之間計算。
對于小于-5的lambda值添加檢查，如果是這種情況，則假定lambda值為1，并且使用原始數據來擬合模型。 lambda值為1與未轉換相同，因此逆變換沒有效果。
下面列出了完整的示例。
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)
運行此示例將打印每次迭代的預測值和預期值
注意，當使用boxcox（）轉換函數時，可能會看到一些警告; 例如：
RuntimeWarning: overflow encountered in square
llf -= N / 2.0 * np.log(np.sum((y - y_mean)**2. / N, axis=0))
不過現在可以忽略這些。
在轉換數據上生成的模型最終的RMSE為49.443。 對于未轉換的數據，這是一個比ARIMA模型更小的誤差，但只是很細微的，它可能是，也可能不是統計意義上的不同。
...
>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變換的模型作為最終模型。
在開發完模型并選擇最終模型后，必須對其進行驗證和確定。
驗證是一個可選的過程，為我們提供了“最后的檢查手段”，以確保我們沒有被自己欺騙。
此部分包括以下步驟：
1. 確定模型：訓練并保存最終模型。
2. 進行預測: 加載最終模型并進行預測。
3. 驗證模型: 加載和驗證最終模型。
模型的確定涉及在整個數據集上擬合ARIMA模型，值得注意的一點是，該步驟中的數據集已經經過了變換。
一旦擬合好之后，模型可以保存到文件以供后續使用。這里由于我們對數據執行了Box-Cox變換，因此我們需要明確所選擇的lamda值，以便使得來自模型的任何預測都可以被轉換回原始的未轉換的尺度。
下面的示例展示了在Box-Cox變換數據集上使用ARIMA（0,1,2）模型，并將整個擬合對象和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])
運行該示例會創建兩個本地文件：
model.pkl這是通過ARIMA.fit() 調用生成的ARIMAResult對象。 包括了各種系數和在擬合模型時返回的所有其他內部數據。
model_lambda.npy這是作為一行一列NumPy數組對象存儲的lambda值。
這些文件或許會保存的過于詳細。其實在操作中真正需要的是：來自模型的AR和MA系數、用于差異數量的d參數、以及滯后的觀察值、模型殘差和用于變換的lamda值。
這里，最直接的方法是加載模型并進行單一預測。
這是相對直接的，涉及恢復保存的模型、lambda值和調用forecast() 方法。
然而，在當前穩定版本的statsmodels庫（v0.6.1）中存在一個Bug，它的錯誤提示如下：
TypeError: __new__() takes at least 3 arguments (1 given)
據我測試發現，這個bug也似乎存在于statsmodels的0.8版本候選1版本中。
更多詳細信息，請參閱Zae Myung Kim的討論和在GitHub上的解決辦法。
鏈接：https://github.com/statsmodels/statsmodels/pull/3217
我們可以使用一個猴子補丁（monkey patch）來解決這個問題，在保存之前將一個 __getnewargs __() 實例函數添加到ARIMA類中。
下面的示例與上一小節的變化相同。
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])
我們現在可以加載模型并進行單一預測。
下面這個示例展示了加載模型，對下一時間步進進行了預測，反轉Box-Cox變換，并打印預測結果的過程。
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)
最后得到的預測值為452。
如果我們在validation.csv中查看的話，我們可以看到下一個時間段的第一行的數值是452。模型達到了100% 的正確率，這是非常驚人的（也是幸運的）。
Predicted: 452.043
在這一步驟中，我們可以加載模型并以模擬操作的方式使用它。
在測試工具部分中，我們將原始數據集的最后12個月保存在單獨的文件中，以驗證最終模型。
我們現在可以加載此validation.csv文件，并用它來檢驗我們的模型對真正“未知”的數據究竟表現如何。
這里可以采用兩種方式：
加載模型并使用它預測未來的12個月。 超過前一兩個月的預測將很快開始降低技能（start to degrade in skill）。
加載模型并以步進預測的方式使用該模型，更新每個時間步長上的變換和模型。 這是首選的方法，在實踐中使用這個模型，它將表現出最佳的性能。
與前面章節中的模型評估一樣，我們將以步進預測的方式進行預測。 這意味著我們將在驗證數據集中逐步引導時間點，并將觀察結果作為歷史記錄的更新。
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()
運行以上示例將打印出在驗證數據集中所有時間步長上的每個預測值和預期值。
修正階段最終的RMSE預測結果為53次搶劫。 這與預期的49次并沒有太大的出入，我希望它在簡單的持久性模型上也能有類似的表現。
>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
我們在這里同時提供了預測結果與驗證數據集之間相互比較的圖表。
該預測具有持續性預測的特性。這表明雖然這個時間序列有明顯的趨勢，但它仍然是一個相當困難的問題。
驗證時間步驟的預測
需要說明的是，本教程并不完美，你可以有很多改善結果的方法。
本節列出了一些或許可行的改進措施。
統計顯著性檢驗。使用統計學測試來檢查不同模型之間的結果差異是否具有統計學的意義。T檢驗將是一個很好的入手點。
使用數據變換的網格進行搜索。通過Box-Cox變換在ARIMA超參數中重復網格搜索過程，并查看是否可以實現一個與此前不同的，并且更好的參數集。
檢查殘差。通過Box-Cox變換檢測最終模型的預測誤差的殘差，以查看是否存在可以進一步優化的偏差或自相關。
精簡模型存儲。嘗試簡化模型保存，僅僅存儲那些必要的系數，而不是整個ARIMAResults對象。
手動處理趨勢。使用線性或非線性模型直接對趨勢進行建模，并明確地從序列中刪除它。如果該趨勢是非線性的，并且更適合非線性建模，那這種方式可能會帶來更好的性能表現。
置信區間。顯示驗證數據集中預測的置信區間。
數據選擇。考慮在沒有前兩年數據的情況下對問題進行建模，并查看這是否對預測技能（forecast skill）有影響。
通過本教程的學習，你可以初步了解到基于Python的時間序列預測項目的具體實現步驟和工具。
我們在本教程中討論了很多內容，具體來說包括以下三個方面：
如何通過性能指標和評估方法開發測試工具，以及如何快速開發基準預測原型和方法。
如何通過時間序列分析來提出最好的模擬預測問題的方法。
如何開發ARIMA模型，并保存它，然后加載它用來對新數據進行預測。
歡迎加入雷鋒網 ML 開發者讀者群，添加“ycopen3”微信號，回復關鍵詞“開發者”，群助拉您進群交流！
雷峰網版權文章，未經授權禁止轉載。詳情見轉載須知。
