Python时间序列预测案例研究：巴尔的摩年度用水量

时间序列预测是一个过程，获得良好预测的唯一方法就是练习这个过程。

在本教程中，您将了解如何使用Python预测巴尔的摩的年用水量。

本教程的学习将为您提供一个框架，介绍处理自己的时间序列预测问题的步骤和工具。

完成本教程后，您将知道：

如何确认您的Python环境并仔细定义时间序列预测问题。
如何创建评估模型的测试框架，开发基准预测，并利用时间序列分析工具来更好地理解您的问题。
如何开发一个自回归整合移动平均模型，将其保存到文件中，然后加载它来预测新的时间步骤。

让我们开始吧。

用Python来研究时间序列预测-巴尔的摩年度用水量照片由安迪·米切尔提供，保留某些权利。

综述

在本教程中，我们将通过一个端到端的时间序列预测项目，从下载数据集和定义问题到训练最终模型并进行预测。

这个项目并不详尽，但是通过系统地处理时间序列预测问题，展示了如何快速获得好的结果。

此项目步骤如下。

环境。
问题描述。
测试框架。
持久性。
数据分析。
ARIMA模型。
模型验证。

这将提供一个能够用来处理自己时序问题的模板。

1.环境

本教程假定已安装并正在运行的SciPy环境及依赖关系，其中包括：

SciPy
NumPy
Matplotlib
Pandas
scikit-learn
状态模型

如果您需要在工作站上安装Python和SciPy环境的帮助，请考虑为您管理大部分内容的Anaconda distribution。

这个脚本将帮助你检查你安装的这些库的版本。

# scipy
import scipy
print('scipy: %s' % scipy.__version__)
# numpy
import numpy
print('numpy: %s' % numpy.__version__)
# matplotlib
import matplotlib
print('matplotlib: %s' % matplotlib.__version__)
# pandas
import pandas
print('pandas: %s' % pandas.__version__)
# scikit-learn
import sklearn
print('sklearn: %s' % sklearn.__version__)
# statsmodels
import statsmodels
print('statsmodels: %s' % 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.问题描述

本文问题是预测年用水量。

数据集提供了从1885年到1963年（即79年的数据）在巴尔的摩的年用水量。

这些数值是以每人每天的升数为单位的，有79个观测值。

该数据集的采集归功于Hipel和McLeod，1994。

您可以了解有关此数据集的更多信息，并直接从DataMarket下载。

将数据集下载为CSV文件，并将其放在当前工作目录中，文件名为 “ water.csv ”。

3.测试框架

我们必须开发一个测试框架来详细了解数据并评估候选模型。

这涉及两个步骤：

定义验证数据集。
开发模型评估方法。

3.1验证数据集

这个数据集中的数据不是最新的，意味着我们不能轻松手机最新的数据来验证评估模型。

因此，我们假装它是1953年的，并保留了过去10年的分析和模型选择的数据。

这最后十年的数据将被用来验证最终的模型。

下面的代码将数据集加载为Pandas系列，并将数据一分为二，一部分 (dataset.csv)用来训练、形成模型，另一部分(validation.csv)用来验证模型效果。

from pandas import Series
series = Series.from_csv('water.csv', header=0)
split_point = len(series) - 10
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')

运行该示例创建两个文件，并在每个文件中显示观察值的数量。

Dataset 69, Validation 10

这些文件的具体内容是：

dataset.csv：从1885年到1953年的观测（69个观测）。
validation.csv：从1954年到1963年的观测结果（10次观测）。

验证数据集大约是原始数据集的12％。

请注意，保存的数据集没有标题行，因此我们不需要在稍后处理这些文件时满足这一点。

3.2。模型评估

模型评估只能在上一节中准备好的dataset.csv中的数据上执行。

模型评估涉及两个要素：

性能指标。
测试策略。

3.2.1性能测量

我们将使用均方根误差（RMSE）评估预测的性能。这将会给予那些严重错误的预测更大的权重值（使得错误预测更加明显），并且和原始数据的单位相同。

对数据的任何转换必须在RMSE被计算和报告之前撤销，以使不同方法之间的性能可以直接比较。

我们可以使用scikit-learn库的帮助器函数mean_squared_error（）来计算RMSE，该函数计算期望值列表（测试集）和预测列表之间的均方差。然后，我们可以取这个值的平方根从而得到一个RMSE分数。

例如：

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

3.2.2测试策略

候选模型将使用前向验证进行评估。

这是因为问题定义需要滚动预测的模型，有了所有可用数据在此处都需要进行一步预测（one-step forecasts）。

前向验证的工作流程如下：

数据集的前50％将被保留以训练模型。
剩余的50％数据集将被迭代并测试模型。
测试数据集的步骤：
- 训练模型。
- 做出一步预测，并将预测值存储起来供后续评估。
- 来自测试数据集的实际观察值将被添加到下一次迭代的训练数据集中。
在测试数据集的列举期间所做的预测将被评估，评估结果将以RMSE报告形式呈现。

鉴于数据量小，我们将允许在每次预测之前对所有可用数据重新训练模型。

我们可以使用简单的NumPy和Python代码编写测试工具的代码。

首先，我们可以直接将数据集分解为训练集和测试集。如果加载的数据仍然有一些String或Integer数据类型，我们会一直小心将加载的数据集转换为float32。

# 准备数据
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)
train, test = X[0:train_size], X[train_size:]

接下来，我们可以迭代测试数据集中的时间。训练数据集存储在一个Python列表中，因为我们需要在每次迭代时轻松地附加一个新的观测值，而NumPy数组连接则感觉太过分了。

通常由模型作出的预测被称为yhat，因为结果或观测被称为y和yhat(y'上面有一个标记)是y变量预测的数学符号。

如果模型存在问题，则每个时间点的预测值和观测值值都会被显示以做一个全面的检查预测。

译者注：每个时间点的observation即它的实际值，在此文中统一译为观测值

# 前向验证
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # 预测
    yhat = ...
    predictions.append(yhat)
    # 观察值
    obs = test[i]
    history.append(obs)
    print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))

4.持续预测

陷入数据分析和建模之前的第一步是建立性能基准。

这将提供一个模板，用于评估模型使用所建议的测试工具和一个性能度量，通过它可以比较所有更复杂的预测模型。时间序列预测的基线预测称为朴素预测或持久性。

时间序列预测的基线预测被称为朴素预测或持续性预测。

在这个方法中，来自前一个时间步（time step）的观测值被用作下一个时间步的预测值。

我们可以直接将其插入到上一节定义的测试框架中。

下面提供了完整的代码清单。

from pandas import Series
from sklearn.metrics import mean_squared_erro
from math import sqrt
# 加载数据
series = Series.from_csv('dataset.csv')
# 准备数据
X = series.values
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]
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))
# 报告性能
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

运行测试框架能输出测试数据集每次迭代的预测值和观测值。

该例子以输出模型的RMSE值为结尾。

在这种情况下，我们可以看到持续性模型的RMSE达到了21.975。这意味着，平均而言，因为每人的日均耗水量RMSE值为22，所以这个模型是错误的。


...
>Predicted=613.000, Expected=598
>Predicted=598.000, Expected=575
>Predicted=575.000, Expected=564
>Predicted=564.000, Expected=549
>Predicted=549.000, Expected=538
RMSE: 21.975

既然我们有了一个基准预测方法和性能; 现在我们可以开始深入挖掘数据了。

5.数据分析

我们可以使用摘要统计和数据图快速了解预测问题的结构。

在本节中，我们将从四个角度来看待数据：

摘要统计。
折线图。
密度图。
箱线图。

5.1。摘要统计

摘要统计数据可以快速查看观测值的极限。它可以帮助快速了解我们正在处理的事情。

以下示例计算并输出时间序列的摘要统计数据。

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

运行该示例提供了一些摘要统计信息以供查看

这些统计的一些观察结果包括：

观察次数（次数）符合我们的预期，这意味着我们正在正确处理数据。
平均值大约是500，我们可以考虑在这个系列的水平。
标准偏差和百分位数表明在均值附近有一个相当严格的传播。

count     69.000000
mean     500.478261
std       73.901685
min      344.000000
25%      458.000000
50%      492.000000
75%      538.000000
max      662.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。箱线图

我们可以将年度数据按十年一个刻度进行分组，并了解每个十年的观测数据传播情况，以及这种情况可能如何变化。

我们希望看到一些趋势（增加的平均数或中位数），但看看其他分布会如何变化可能会很有趣。

下面的例子将观测值按十年分组，并为每个十年观测值创建一个箱线图。过去的十年实际上只包含九年的数据，并且可能不能和其它十年数据成为一个有效对比。因此只有1885年至1944年间的数据被绘制出来.

 
from pandas import Series
from pandas import DataFrame
from pandas import TimeGroupe
from matplotlib import pyplot
series = Series.from_csv('dataset.csv')
groups = series['1885':'1944'].groupby(TimeGrouper('10AS'))
decades = DataFrame()
for name, group in groups:
    decades[name.year] = group.values
decades.boxplot()
pyplot.show()

运行该示例创建了6个箱线图。

从下图中我们可以看到：

每年的中值（红线）可能会呈现增加的趋势，但不一定是线性的。
数据箱（蓝色框）的扩散或中间50％的数据确实显示出一些变化。
几十年来，可能有一些异常值（箱体图外面的十字架）。
倒数第二个十年的平均消费似乎较低，也许与第一次世界大战有关。

每年对数据的看法是一种有趣的途径，可以通过观察十年到十年的简要统计数据和汇总统计数据的变化来进一步研究。

6. ARIMA模型

在本节中，我们将针对该问题开发自回归整数滑动平均模型，即ARIMA模型。

我们将通过手动和自动配置ARIMA模型来进行建模。接下来的第三步是获取被选中模型的残差值。

因此，本节分为三个步骤：

手动配置ARIMA。
自动配置ARIMA。
查看残差值。

6.1手动配置ARIMA

ARIMA（p，d，q）模型需要三个参数，通常需手动配置。

在时需分析中，一般假设我们使用的是平稳时间序列

时间序列可能是非平稳的。我们可以首先差分化时序并使用统计测试来检查以确保时序已经被转换成平稳时间序列。

下面的例子将时序平稳化，并将其保存到文件stationary.csv。

from pandas import Series
from statsmodels.tsa.stattools import adfulle
from matplotlib import pyplot
 
# 创建差分化函数
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
X = X.astype('float32')
# 差分化数据
stationary = difference(X)
stationary.index = series.index[1:]
# 检查时序是否已经变得平稳
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))
# 绘制平稳时序的图
stationary.plot()
pyplot.show()
# 保存
stationary.to_csv('stationary.csv')

运行该示例输出的结果是统计显著性测试的结果，用来检测时间序列是否平稳。具体来说，是扩增的Dickey-Fuller测试。

结果表明，检验统计值-6.126719小于-3.534的1％临界值。这表明我们可以拒绝显着性水平小于1％的零假设（即结果是统计错误的概率很低）。

拒绝零假设意味着这个过程没有单位根，反过来说，时间序列是平稳的或者没有时间相关的结构。

ADF Statistic: -6.126719
p-value: 0.000000
Critical Values:
    5%: -2.906
    1%: -3.534
    10%: -2.591

这表明至少需要一个差异水平。我们的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没有显着的滞后。
PACF也没有显着的滞后。

那么对于p和q，一个好的开始点都是0.

这个快速分析表明，对于原始数据，ARIMA（0,1,0）模型可能是一个很好的起点。

这实际上是一个持久性模型（朴素预测模型）。完整的例子如下所示。

from pandas import Series
from sklearn.metrics import mean_squared_erro
from statsmodels.tsa.arima_model import ARIMA
from math import sqrt
# 加载数据
series = Series.from_csv('dataset.csv')
# 准备数据
X = series.values
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]
predictions = list()
for i in range(len(test)):
    # predict
    model = ARIMA(history, order=(0,1,0))
    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))
# 报告性能
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)

运行此示例将导致RMSE为22.311，比上面的持久性模型略高。

这可能是ARIMA模型实施的细节所导致，例如被计算和添加的自动趋势常数。

...
>Predicted=617.079, Expected=598
>Predicted=601.781, Expected=575
>Predicted=578.369, Expected=564
>Predicted=567.152, Expected=549
>Predicted=551.881, Expected=538
RMSE: 22.311

6.2网格搜索ARIMA超参数

ACF和PACF图表显示我们不能比这个数据集上的持久性模型做得更好。

为了确认这一分析，我们可以网格搜索一套ARIMA超参数，并检查有无模型可以在测试集上做出更好的基于RMSE性能的预测。

在本节中，我们将搜索p，d和q的值作为组合（跳过那些不能汇集的组合），并找出导致最佳性能的组合。我们将使用网格搜索来探索整数值子集中的所有组合。

具体来说，我们将搜索以下参数的所有组合：

p：0到4。
d：0到2。
q：0到4。

这是（5 * 3 * 5）的测试，或测试框架的300次潜在运行，我们需要一些时间来执行。

当调用fit（）时，我们还将禁止从模型中自动添加一个趋势常量，方法是将“ trend”参数设置为“ nc ” 。

下面列出了网格搜索版本测试工具的完整示例。

import warnings
from pandas import Series
from statsmodels.tsa.arima_model import ARIMA
from sklearn.metrics import mean_squared_erro
from math import sqrt
 
# 基于给定的ARIMA(p,d,q)值评评估该模型并返回RMSE值
def evaluate_arima_model(X, arima_order):
    # 准备训练集数据
    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]
    # 作出预测
    predictions = list()
    for t in range(len(test)):
        model = ARIMA(history, order=arima_order)
        # model_fit = model.fit(disp=0)
        model_fit = model.fit(trend='nc', disp=0)
        yhat = model_fit.forecast()[0]
        predictions.append(yhat)
        history.append(test[t])
    # 计算测试集上的预测误差
    mse = mean_squared_error(test, predictions)
    rmse = sqrt(mse)
    return rmse
 
# 评估带有不同参数值的每个ARIMA模型
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, orde
                    print('ARIMA%s RMSE=%.3f' % (order,mse))
                except:
                    continue
    print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score))
 
# 加载数据
series = Series.from_csv('dataset.csv')
# 评价参数
p_values = range(0, 5)
d_values = range(0, 3)
q_values = range(0, 5)
warnings.filterwarnings("ignore")
evaluate_models(series.values, p_values, d_values, q_values)

运行该示例将运行所有组合，并将结果报告在没有错误的聚合上。这个例子需要超过2分钟才能在现代硬件上运行。

结果表明，发现的最佳配置是ARIMA（2，1，0），RMSE为21.733，略低于之前测试的手动持久性模型，但可能预测结果差异不大。

...
ARIMA(4, 1, 0) RMSE=24.802
ARIMA(4, 1, 1) RMSE=25.103
ARIMA(4, 2, 0) RMSE=27.089
ARIMA(4, 2, 1) RMSE=25.932
ARIMA(4, 2, 2) RMSE=25.418
Best ARIMA(2, 1, 0) RMSE=21.733

我们将选择这个ARIMA（2，1，0）模型。

6.3查看残差

一个好的模型最终检查是检查预测的残差值

理想情况下，残差的分布应该是均值为0的高斯分布。

我们可以通过使用摘要统计和图来检查ARIMA（2，1，0）模型的残差。下面的例子计算并总结了残差预测误差。

from pandas import Series
from pandas import DataFrame
from sklearn.metrics import mean_squared_erro
from statsmodels.tsa.arima_model import ARIMA
from math import sqrt
from matplotlib import pyplot
# 加载数据
series = Series.from_csv('dataset.csv')
# 准备数据
X = series.values
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]
predictions = list()
for i in range(len(test)):
    # 预测
    model = ARIMA(history, order=(2,1,0))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = model_fit.forecast()[0]
    predictions.append(yhat)
    # 观测
    obs = test[i]
    history.append(obs)
# 误差
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = DataFrame(residuals)
print(residuals.describe())
pyplot.figure()
pyplot.subplot(211)
residuals.hist(ax=pyplot.gca())
pyplot.subplot(212)
residuals.plot(kind='kde', ax=pyplot.gca())
pyplot.show()

首先运行示例描述了残差的分布。

我们可以看到这个分布有一个右移，这个均值1.081624，是是非零的。

这也许是预测有偏差的一个信号。

count  35.000000
mean    1.081624
std    22.022566
min   -52.103811
25%   -16.202283
50%    -0.459801
75%    12.085091
max    51.284336

残差的分布也被绘制出来。

这些图表显示了具有较长右尾的高斯分布，提供了进一步的证据表明。权重转换可能是值得探索的

我们可以使用这些信息来对每个预测加上平均残差1.081624来对预测进行偏差修正。

以下示例可以达到偏差修正的目的。

from pandas import Series
from pandas import DataFrame
from sklearn.metrics import mean_squared_erro
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-forward validation
history = [x for x in train]
predictions = list()
bias = 1.081624
for i in range(len(test)):
    # predict
    model = ARIMA(history, order=(2,1,0))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = bias + float(model_fit.forecast()[0])
    predictions.append(yhat)
    # observation
    obs = test[i]
    history.append(obs)
# report performance
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)
# 总结残差
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = DataFrame(residuals)
print(residuals.describe())
# 绘制残差
pyplot.figure()
pyplot.subplot(211)
residuals.hist(ax=pyplot.gca())
pyplot.subplot(212)
residuals.plot(kind='kde', ax=pyplot.gca())
pyplot.show()

预测的表现从21.733略微改善到21.706，但这并不一定很重要

预测残差的总结表明，平均值的确移到了非常接近零的值。

RMSE: 21.706
                  0
count  3.500000e+01
mean  -3.537544e-07
std    2.202257e+01
min   -5.318543e+01
25%   -1.728391e+01
50%   -1.541425e+00
75%    1.100347e+01
max    5.020271e+01

最后，剩差的密度图显示向零的小偏移。

这个偏差修正是否值得这个问题值得商榷，但是我们现在就会使用它。

7.模型验证

模型开发完成并选定最终模型后，必须进行验证和确定，

验证是该过程的一个可选部分，但提供“最后检查”以确保模型正确。

本节包括以下步骤：

敲定模型：训练并保存最终模型。
做出预测：加载最终模型并进行预测。
验证模型：加载并验证最终模型。

7.1完成模型

完成模型涉及在整个数据集上拟合ARIMA模型，在这种情况下，是在整个数据集的转换版本上。

一旦匹配，此模型就可以被保存并在以后使用。

以下示例在数据集上训练ARIMA（2,1,0）模型，并保存整个拟合对象和偏差到文件中。

当前稳定版本的statsmodels库（v0.6.1）中存在一个错误，当您尝试从文件加载保存的ARIMA模型时会导致错误。报告的错误是：

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

当我测试它时，这个bug也出现在0.8版statsmodels的候选版本1中。有关更多详细信息，请参阅Zae Myung Kim 关于此GitHub问题的讨论和修复。

我们可以使用一个monkey补丁来解决这个问题，保存之前在ARIMA类中添加一个__getnewargs __（）实例函数。

下面的示例将适配模型保存为正确的状态，以便稍后可以成功加载。

from pandas import Series
from statsmodels.tsa.arima_model import ARIMA
from scipy.stats import boxcox
import numpy
 
# RRIMA型模型的monkey补丁
def __getnewargs__(self):
    return ((self.endog),(self.k_lags, self.k_diff, self.k_ma))
 
ARIMA.__getnewargs__ = __getnewargs__
 
# 加载数据
series = Series.from_csv('dataset.csv')
# 准备数据
X = series.values
X = X.astype('float32')
# 拟合模型
model = ARIMA(X, order=(2,1,0))
model_fit = model.fit(trend='nc', disp=0)
# bias constant, could be calculated from in-sample mean residual
bias = 1.081624
# 保存模型
model_fit.save('model.pkl')
numpy.save('model_bias.npy', [bias])

运行示例创建两个本地文件：

model.pkl。这是从ARIMA.fit（）调用的ARIMAResult对象。这包括拟合模型时返回的系数和所有其他内部数据。
model_bias.npy这是存储为一行，一列NumPy数组的偏置值。

7.2做出预测

一个自然情况可能是加载模型并进行单一的预测。

这是相对简单的，涉及恢复保存的模型和偏差，并调用forecast（） 函数。

下面的示例加载模型，对下一个时间步（时间点）进行预测，并输出预测。

from pandas import Series
from statsmodels.tsa.arima_model import ARIMAResults
import numpy
model_fit = ARIMAResults.load('model.pkl')
bias = numpy.load('model_bias.npy')
yhat = bias + float(model_fit.forecast()[0])
print('Predicted: %.3f' % yhat)

运行示例输出的预测约为540.

Predicted: 540.013

如果我们看一看validation.csv，我们可以看到下一个时间段的第一行的值是568.预测是正确的。

7.3验证模型

我们可以加载模型并以伪装的操作方式使用它。

在测试工具部分，我们将原始数据集的最后10年保存在一个单独的文件中，以验证最终模型。

我们现在可以加载这个validation.csv文件，并使用它来检查我们的模型对“看不见的”数据的有效性。

有两种方法可以进行：

加载模型并使用它来预测未来10年。超过头一年或两年的预测很快就会开始降低技能。
加载模型并以滚动预测方式使用它，更新每个时间步的变换和模型。这是首选的方法，因为这个方法可以可以让我们看到这个模型是如何在实践中应用并达到最佳性能。

与前面章节中的模型评估一样，我们将以滚动预测的方式进行预测。这意味着我们将在验证数据集中跨越前置时间，并将观察结果作为历史更新。

from pandas import Series
from matplotlib import pyplot
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.arima_model import ARIMAResults
from sklearn.metrics import mean_squared_erro
from math import sqrt
import numpy
# 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')
bias = numpy.load('model_bias.npy')
# 做出第一个预测
predictions = list()
yhat = bias + float(model_fit.forecast()[0])
predictions.append(yhat)
history.append(y[0])
print('>Predicted=%.3f, Expected=%3.f' % (yhat, y[0]))
# 滚动预测
for i in range(1, len(y)):
    # predict
    model = ARIMA(history, order=(2,1,0))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = bias + float(model_fit.forecast()[0])
    predictions.append(yhat)
    # observation
    obs = y[i]
    history.append(obs)
    print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# 报告性能
mse = mean_squared_error(y, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)
pyplot.plot(y)
pyplot.plot(predictions, color='red')
pyplot.show()

运行示例输出验证数据集中每个时间步预测和期望值。

有效期的最终RMSE预计为每人每天16升。这与21的预期误差没有太大差别，但是我认为它与简单的持续型模型没有太大的区别。

>Predicted=540.013, Expected=568
>Predicted=571.589, Expected=575
>Predicted=573.289, Expected=579
>Predicted=579.561, Expected=587
>Predicted=588.063, Expected=602
>Predicted=603.022, Expected=594
>Predicted=593.178, Expected=587
>Predicted=588.558, Expected=587
>Predicted=588.797, Expected=625
>Predicted=627.941, Expected=613
RMSE: 16.532

我们还提供了与验证数据集相比的预测图。

预测确实具有持续性预测的特征。这表明，虽然这个时间序列确实有一个明显的趋势，但仍然是一个相当困难的问题。

概要

在本教程中，您了解了使用Python进行时间序列预测项目需要的步骤和工具。

我们在本教程中介绍了很多内容。特别：

如何开发性能测量和评估方法的测试工具，以及如何快速开发基准预测和性能评估
利用如何使用时序分析来提出如何优化预测模型的问题
如何开发一个ARIMA模型，保存，然后加载它来预测新的数据。