Jason Brownlee：时间序列预测法是一个过程，而获得良好预测结果的唯一途径是实践这个过程。
在本教程中，您将了解如何利用Python语言来预测波士顿每月持械抢劫案发生的数量。
本教程所述为您提供了一套处理时间序列预测问题的框架，包括方法步骤和工具，通过实践，可以用它来解决自己遇到的相关问题。
本教程结束之后，您将了解：
如何核查Python环境并准确地定义一个时间序列预测问题。
如何构建一套测试工具链，用于评估模型，开发预测原型。以及如何通过时间序列分析工具更好地理解你的问题。
如何开发自回归积分滑动平均模型（ARIMA），将其保存到文件，并在之后加载它对新的时间步骤进行预测。
让我们开始吧。
波士顿
概述
在本教程中，我们将端到端地来解析一个时间序列预测工程，从下载数据集、定义问题到训练出最终模型并进行预测。
该工程并不面面俱到，但展示了如何通过系统性地处理时间序列预测问题，来快速获得好的结果。
我们将通过如下步骤来解析该工程 ：
环境配置.
问题描述
测试工具链
持续性.
数据分析
ARIMA 模型
模型验证
该部分提供一个解决时间序列预测问题的模板，你可以套用自己的数据集。
1. 环境配置
本教程假设在已安装好了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
2. 问题描述
我们研究的问题是预测美国波士顿每月持械抢劫案的数量。
该数据集提供了从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”。
3.测试工具链
我们必须开发一套测试工具链来审查数据和评估候选模型。
这包含两个步骤：
1.定义好验证数据集。
2.开发一套模型评估方法。
3.1 验证数据集
这个数据集不是当前时间段的数据集。 这意味着我们不能轻轻松松收集更新后的数据来验证该模型。
因此，我们假设现在是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％。
请注意，由于保存的数据集中没有标题行，因此，稍后处理这些文件时，我们也不需要满足此要求。
3.2. 模型验证
模型评估将仅对上一节中准备的dataset.csv中的数据进行处理。
模型评估涉及两个要素：
1. 评价指标。
2. 测试策略。
3.2.1评价指标
这些观察值是抢劫案的计数。
我们将使用均方根误差（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)3.2.2 测试策略
将使用步进验证的方式(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))4. Persistence(持续性模型)
在数据分析和建模陷入困境之前的第一步是建立一个评估原型。
这里将为利用测试工具进行模型评估和性能测量，分别提供一套模版。通过该性能测量，可以将当前模型和其他更复杂的预测模型进行比较。
时间序列预测的预测原型被称为天真预测或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我们现在有了一个预测方法和性能评估的原型; 现在我们可以开始深入数据。
5. 数据分析
我们可以使用统计汇总和数据的绘图，来快速了解该预测问题的数据结构。
在本节中，我们将从四个角度来看数据：
1.摘要统计。
2.线图。
3.密度图。
4.盒型-虚线图。
5.1 摘要统计
利用文本编辑器中打开数据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
对于序列中比较大的变化将可能模型很难进行高度准确的预测，如果它是由随机波动（例如，非系统的）引起，
5.2 线图
画一个时间序列的线图可以对该问题提供很多深入的信息。
下面的示例创建并显示出数据集的线图。
from pandas import Series
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
series.plot()
pyplot.show()运行示例并查看绘图。 注意系列中的任何明显的时间结构。
图中的一些观察结果包括：
随着时间的推移，劫案有增加的趋势。
似乎没有任何明显的异常值。
每年都有相对较大的波动，上下波动。
后几年的波动似乎大于前几年的波动。
整体趋势表明该数据集几乎肯定不是稳定的，并且波形的明显变化也可能对此有所提示。
每月波士顿抢劫案数据线图
这些简单的观察结果表明，我们可以看到对趋势进行建模和从时间序列中将它提取出来的好处。
或者，我们可以使用差分法从而使序列的建模平稳。如果后期波动有增长趋势，我们甚至可能需要两个差分到不同级别。
5.3 密度图
提取出密度图，可以提供对数据结构的进一步了解。
下面的示例创建了没有任何时间结构的观测值的直方图和密度图。
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()运行示例并查看绘图。
图中展示的一些分析结果包括：
分布不是高斯分布。
分布是左移的，可以是指数或双高斯分布
每月波士顿抢劫密案数量密度图
我们可以看到在建模之前进行一些数据的维度变换的一些益处。
5.4 盒形-虚线图
我们可以按年份对月度数据进行分组，并了解每年观察值的分布情况，以及观察值可能如何变化。
我们确实希望从数据中看到一些趋势（增加平均值或中位数），但看到其余的分布是如何变化的可能会更有趣。
下面的示例将观察值按照年份划分，并为每个观察年创建一个盒型-虚线图。 最后一年（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年，与数据集的其余部分有很大的不同。
波士顿每月抢劫案数目的盒形-虚线图
观察结果表明，年度波动可能不是系统性的，难以模拟。 他们还建议将建模的前两年的数据从模型分离出来可能有一些好处，如果确实不同。
这种对数据进行年度观察是一个有趣的方法，而且可以通过查看年度汇总统计数据和其中的变化，来进一步跟进未来的变化。
接下来，我们可以开始查看该序列的预测模型。
6. ARIMA 模型
在本节中，我们将会针对该问题开发自回归积分滑动平均模型，即ARIMA模型。
我们将其分四个步骤：
1.开发一套手动配置的ARIMA模型。
2.使用ARIMA的网格搜索来找到优化的模型。
3.预测残差的分析，以评估模型中的所有偏差。
4.使用能量变换探索模型的改进。
6.1 手动配置ARIMA
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模型。
6.2 网格搜索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
6.3 检查残差
对一个模型来说，一个不错的最终检查，是检查其残差预测误差。
理想情况下，残差的分布应该是符合零均值的高斯分布。
我们可以通过用直方图和密度图绘制残差来检查这一点。
以下示例为计算测试集上预测的残差，并创建该密度图。
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图
6.4 Box-Cox 变换数据集
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变换的模型作为最终模型。
7. 模型验证
在开发完模型并选择最终模型后，必须对其进行验证和确定。
验证是一个可选的过程，为我们提供了“最后的检查手段”，以确保我们没有被自己欺骗。
此部分包括以下步骤：
1. 确定模型：训练并保存最终模型。
2. 进行预测: 加载最终模型并进行预测。
3. 验证模型: 加载和验证最终模型。
7.1 确定模型
模型的确定涉及在整个数据集上拟合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值。
7.2 进行预测
这里，最直接的方法是加载模型并进行单一预测。
这是相对直接的，涉及恢复保存的模型、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
7.3 验证模型
在这一步骤中，我们可以加载模型并以模拟操作的方式使用它。
在测试工具部分中，我们将原始数据集的最后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模型，并保存它，然后加载它用来对新数据进行预测。
