机器学习笔记 - 时间序列预测研究：法国香槟的月销量

1、问题描述

        问题是预测 Perrin Freres 标签（以法国的一个地区命名）的香槟月销量。该数据集提供了从 1964 年 1 月到 1972 年 9 月的香槟月销售量，或不到 10 年的数据。这些值是对数百万销售额的计数，有 105 个观察值。

链接：https://pan.baidu.com/s/1DyoZ_xFZeItCfrpX1RTG2g 
提取码：1f2u

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

2、划分数据集

        假设现在是 1971 年 9 月，并在分析和模型选择中保留最后一年的数据。最后一年的数据将用于验证最终模型

        下面的代码会将数据集作为 Pandas 系列加载并分成两部分，一个用于模型开发（dataset.csv），另一个用于验证（validation.csv）。

from pandas import read_csv
series = read_csv('champagne.csv', header=0, index_col=0, parse_dates=True, squeeze=True)
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', header=False)
validation.to_csv('validation.csv', header=False)

        运行该示例会创建两个文件并打印每个文件中的数据数量。

Dataset 93, Validation 12

        这些文件的具体内容是：

        dataset.csv：从 1964 年 1 月到 1971 年 9 月的观测值（93 个观测值）

        validation.csv：从 1971 年 10 月到 1972 年 9 月的观察（12 个观察）

        验证数据集大约是原始数据集的 11%。

3、数据分析

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

（1）汇总统计

        汇总统计数据可以快速查看观察值的限制。它有助于快速了解我们正在使用的内容。

from pandas import read_csv
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
print(series.describe())

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

        观察次数（计数）符合我们的预期，这意味着我们正在正确处理数据。

        平均值约为 4,641，我们可能会考虑我们在本系列中的水平。

        标准差（与平均值的平均差）相对较大，为 2,486 次销售。

        百分位数和标准偏差确实表明数据存在很大的差异。

count       93.000000
mean      4641.118280
std       2486.403841
min       1573.000000
25%       3036.000000
50%       4016.000000
75%       5048.000000
max      13916.000000

（2）线图

from pandas import read_csv
from matplotlib import pyplot
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
series.plot()
pyplot.show()

        该图的一些观察结果包括：

        随着时间的推移，销售额可能会呈上升趋势。

        每年的销售似乎都存在系统性的季节性。

        季节性信号似乎随着时间的推移而增长，表明存在乘法关系（增加变化）。

        似乎没有任何明显的异常值。

        季节性表明该系列几乎可以肯定是非平稳的。

         显式建模季节性组件并将其移除可能会有好处。您还可以探索使用具有一个或两个级别的差分来使系列静止。

        季节性成分的增加趋势或增长可能表明使用对数或其他幂变换。

（3）季节线性图

        我们可以通过逐年观察数据集的线图来确认季节性是一年周期的假设。

        下面的示例将 7 整年的数据作为单独的组，并为每个组创建一个线图。线图垂直对齐，以帮助发现任何年度模式。

from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
groups = series['1964':'1970'].groupby(Grouper(freq='A'))
years = DataFrame()
pyplot.figure()
i = 1
n_groups = len(groups)
for name, group in groups:
	pyplot.subplot((n_groups*100) + 10 + i)
	i += 1
	pyplot.plot(group)
pyplot.show()

        运行该示例会创建 7 个线图的堆栈。

        我们可以清楚地看到每年 8 月的下降和每年 8 月到 12 月的上升。这种模式每年都出现相同的情况，尽管水平不同。

        这将有助于以后任何明确的基于季节的建模。

（4）密度图

        查看观察密度图可以进一步了解数据的结构。

from pandas import read_csv
from matplotlib import pyplot
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
pyplot.figure(1)
pyplot.subplot(211)
series.hist()
pyplot.subplot(212)
series.plot(kind='kde')
pyplot.show()

        图中的一些观察结果包括：

        分布不是高斯分布。

        该形状有一条长长的右尾，可能表明呈指数分布

 （5）箱线图

        我们可以按年份对月度数据进行分组，并了解每年观察的分布情况以及这种情况可能会如何变化。

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

        下面的示例按年份对观察结果进行分组，并为每一年的观察结果创建一个箱须图。去年（1971 年）仅包含 9 个月，可能无法与其他年份的 12 个月观察结果进行有用的比较。因此，仅绘制了 1964 年至 1970 年之间的数据。

from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
groups = series['1964':'1970'].groupby(Grouper(freq='A'))
years = DataFrame()
for name, group in groups:
	years[name.year] = group.values
years.boxplot()
pyplot.show()

        运行该示例会创建 7 个并排的箱线图，每一个对应 7 年的选定数据。

        回顾这些图的一些观察结果包括：

        每年的中值（红线）可能呈上升趋势。

        分散或中间 50% 的数据（蓝色框）确实看起来相当稳定。

        每年都有异常值（黑色十字）；这些可能是季节性周期的顶部或底部。

        去年，1970 年，看起来确实与往年的趋势不同

         观察结果表明，多年来可能存在一些增长趋势以及可能是季节性周期一部分的异常值。

4、ARIMA 模型

        在本节中，我们将针对该问题开发自回归综合移动平均线或 ARIMA 模型。

        我们将通过手动和自动配置 ARIMA 模型来进行建模。接下来是调查所选模型的残差的第三步。

（1）手动配置的 ARIMA

        ARIMA( p,d,q ) 模型需要三个参数，并且传统上是手动配置的。

        时间序列数据的分析假设我们使用的是固定时间序列。

        时间序列几乎可以肯定是非平稳的。我们可以通过首先对序列进行差分并使用统计测试来确认结果是平稳的，从而使其平稳。

        该系列的季节性似乎是逐年的。季节性数据可以通过从上一个周期的同一时间减去观测值来区分，在这种情况下是上一年的同一个月。这确实意味着我们将失去第一年的观察结果，因为与前一年没有区别。

        下面的示例创建了该系列的去季节化版本并将其保存到文件stationary.csv中。

from pandas import read_csv
from pandas import Series
from statsmodels.tsa.stattools import adfuller
from matplotlib import pyplot
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return Series(diff)
 
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
X = series.values
X = X.astype('float32')
# difference data
months_in_year = 12
stationary = difference(X, months_in_year)
stationary.index = series.index[months_in_year:]
# 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', header=False)
# plot
stationary.plot()
pyplot.show()

        运行该示例会输出差分序列是否平稳的统计显着性检验结果。具体来说，增强的 Dickey-Fuller 测试。

        结果表明，检验统计量值 -7.134898 小于临界值-3.515 的 1%。这表明我们可以拒绝显着性水平小于 1% 的原假设（即，结果是统计侥幸的概率很低）。

        拒绝原假设意味着过程没有单位根，进而表明时间序列是平稳的或不具有时间相关结构。

ADF Statistic: -7.134898
p-value: 0.000000
Critical Values:
	5%: -2.898
	1%: -3.515
	10%: -2.586

        作为参考，可以通过添加前一年同月的观测值来反转季节差异操作。如果预测是通过适合季节性差异数据的模型进行的，则需要这样做。为了完整起见，下面列出了反转季节差异运算的函数。

# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]

        还创建了差异数据集的图。

        该图没有显示任何明显的季节性或趋势，表明季节性差异数据集是建模的良好起点。

        我们将使用这个数据集作为 ARIMA 模型的输入。它还表明可能不需要进一步的差分，并且可以将 d 参数设置为 0。

         接下来的第一步是分别选择自回归 (AR) 和移动平均 (MA) 参数的滞后值p和q。

        我们可以通过查看自相关函数 (ACF) 和偏自相关函数 (PACF) 图来做到这一点。

        请注意，我们现在使用季节性差异的stationary.csv 作为我们的数据集。

        下面的示例为该系列创建 ACF 和 PACF 图。

from pandas import read_csv
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from matplotlib import pyplot
series = read_csv('stationary.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
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 显示显着滞后 1 个月，可能在 12 个月和 13 个月有显着滞后。
        ACF 和 PACF 在同一点均显示下降，这可能表明 AR 和 MA 的混合。
        p 和 q 值的良好起点也是 1。

        PACF 图还表明差异数据中仍然存在一些季节性。

        我们可能会考虑更好的季节性模型，例如直接对其进行建模并明确地将其从模型中移除，而不是季节性差异。

         这种快速分析表明，在固定数据上的 ARIMA(1,0,1) 可能是一个很好的起点。

        在拟合每个 ARIMA 模型之前，历史观测值将存在季节性差异。对于所做的所有预测，差异将被反转，以使它们与原始销售计数单位中的预期观察结果直接可比。

        实验表明，ARIMA 的这种配置不会收敛，并导致底层库出错。进一步的实验表明，向平稳数据添加一级差分可使模型更加稳定。该模型可以扩展到 ARIMA(1,1,1)。

        我们还将禁用模型中趋势常数的自动添加，方法是将“ trend ”参数设置为“ nc ”，因为在 fit() 调用中没有常数。通过实验，我发现这可以在某些问题上产生更好的预测性能。

        下面的示例演示了此 ARIMA 模型在测试工具上的性能。

# evaluate manually configured ARIMA model
from pandas import read_csv
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima.model import ARIMA
from math import sqrt
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# 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)):
	# difference data
	months_in_year = 12
	diff = difference(history, months_in_year)
	# predict
	model = ARIMA(diff, order=(1,1,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	yhat = inverse_difference(history, yhat, months_in_year)
	predictions.append(yhat)
	# observation
	obs = test[i]
	history.append(obs)
	print('>Predicted=%.3f, Expected=%.3f' % (yhat, obs))
# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)

        运行此示例导致 RMSE 为 956.942，这明显优于 3186.501 的持久性 RMSE。

>Predicted=3157.018, Expected=5010
>Predicted=4615.082, Expected=4874
>Predicted=4624.998, Expected=4633
>Predicted=2044.097, Expected=1659
>Predicted=5404.428, Expected=5951
RMSE: 956.942

        可以通过配置更好的 ARIMA 模型获得改进的结果。

（2）网格搜索 ARIMA 超参数

        ACF 和 PACF 图表明 ARIMA(1,0,1) 或类似的可能是我们能做的最好的。

        为了确认这一分析，我们可以对一组 ARIMA 超参数进行网格搜索，并检查是否没有模型能够产生更好的样本 RMSE 性能。

        在本节中，我们将在p、d和q的值中搜索组合（跳过那些无法收敛的组合），并找到在测试集上产生最佳性能的组合。我们将使用网格搜索来探索整数值子集中的所有组合。

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

        p：0 到 6。
        d：0 到 2。
        q：0 到 6。
        这是 ( 7 * 3 * 7 ) 或 147 次测试工具的潜在运行，需要一些时间来执行。

        评估滞后为 12 或 13 的 MA 模型可能会很有趣，因为通过审查 ACF 和 PACF 图可能会很有趣。实验表明，这些模型可能不稳定，导致底层数学库出现错误。

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

# grid search ARIMA parameters for time series
import warnings
from pandas import read_csv
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
from math import sqrt
import numpy
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return numpy.array(diff)
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# 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)):
		# difference data
		months_in_year = 12
		diff = difference(history, months_in_year)
		model = ARIMA(diff, order=arima_order)
		model_fit = model.fit()
		yhat = model_fit.forecast()[0]
		yhat = inverse_difference(history, yhat, months_in_year)
		predictions.append(yhat)
		history.append(test[t])
	# calculate out of sample error
	rmse = sqrt(mean_squared_error(test, predictions))
	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:
					rmse = evaluate_arima_model(dataset, order)
					if rmse < best_score:
						best_score, best_cfg = rmse, order
					print('ARIMA%s RMSE=%.3f' % (order,rmse))
				except:
					continue
	print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score))
 
# load dataset
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# evaluate parameters
p_values = range(0, 7)
d_values = range(0, 3)
q_values = range(0, 7)
warnings.filterwarnings("ignore")
evaluate_models(series.values, p_values, d_values, q_values)

        运行该示例会遍历所有组合，并报告那些没有错误收敛的结果。该示例在现代硬件上运行需要 2 分钟多一点。

        结果表明，发现的最佳配置是 RMSE 为 939.464 的 ARIMA(0, 0, 1)，略低于上一节中手动配置的 ARIMA。这种差异可能具有统计意义，也可能不具有统计意义。

ARIMA(5, 1, 2) RMSE=1003.200
ARIMA(5, 2, 1) RMSE=1053.728
ARIMA(6, 0, 0) RMSE=996.466
ARIMA(6, 1, 0) RMSE=1018.211
ARIMA(6, 1, 1) RMSE=1023.762
Best ARIMA(0, 0, 1) RMSE=939.464

（3）查看预测错误

        一个好的模型最终检查是检查剩余的预测误差。

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

        我们可以通过使用汇总统计数据和绘图来检查 ARIMA(0, 0, 1) 模型的残差。下面的示例计算并总结了剩余预测误差。

# summarize ARIMA forecast residuals
from pandas import read_csv
from pandas import DataFrame
from statsmodels.tsa.arima.model import ARIMA
from matplotlib import pyplot
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# 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)):
	# difference data
	months_in_year = 12
	diff = difference(history, months_in_year)
	# predict
	model = ARIMA(diff, order=(0,0,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	yhat = inverse_difference(history, yhat, months_in_year)
	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)
print(residuals.describe())
# plot
pyplot.figure()
pyplot.subplot(211)
residuals.hist(ax=pyplot.gca())
pyplot.subplot(212)
residuals.plot(kind='kde', ax=pyplot.gca())
pyplot.show()

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

        我们可以看到分布有一个右移，并且平均值在 165.904728 处非零。

        这可能是预测存在偏差的迹象。

count    47.000000
mean    165.904728
std     934.696199
min   -2164.247449
25%    -289.651596
50%     191.759548
75%     732.992187
max    2367.304748

        我们可以通过将 165.904728 的平均残差添加到所做的每个预测中来使用此信息来偏置正确的预测。

        下面的示例执行这种偏差相关。

# plots of residual errors of bias corrected forecasts
from pandas import read_csv
from pandas import DataFrame
from statsmodels.tsa.arima.model import ARIMA
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# 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 = 165.904728
for i in range(len(test)):
	# difference data
	months_in_year = 12
	diff = difference(history, months_in_year)
	# predict
	model = ARIMA(diff, order=(0,0,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	yhat = bias + inverse_difference(history, yhat, months_in_year)
	predictions.append(yhat)
	# observation
	obs = test[i]
	history.append(obs)
# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)
# errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = DataFrame(residuals)
print(residuals.describe())
# plot
pyplot.figure()
pyplot.subplot(211)
residuals.hist(ax=pyplot.gca())
pyplot.subplot(212)
residuals.plot(kind='kde', ax=pyplot.gca())
pyplot.show()

        预测的性能从 939.464 略微提高到 924.699，可能显着也可能不显着。

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

RMSE: 924.699
 
count  4.700000e+01
mean   4.965016e-07
std    9.346962e+02
min   -2.330152e+03
25%   -4.555563e+02
50%    2.585482e+01
75%    5.670875e+02
max    2.201400e+03

        最后，残差的密度图确实显示出向零的小幅偏移。

        这种偏差校正是否值得值得商榷，但我们现在将使用它。

         检查任何类型的自相关的残差的时间序列也是一个好主意。如果存在，则表明该模型有更多机会对数据中的时间结构进行建模。

        下面的示例重新计算残差并创建 ACF 和 PACF 图以检查任何显着的自相关。

# ACF and PACF plots of residual errors of bias corrected forecasts
from pandas import read_csv
from pandas import DataFrame
from statsmodels.tsa.arima.model import ARIMA
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# 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)):
	# difference data
	months_in_year = 12
	diff = difference(history, months_in_year)
	# predict
	model = ARIMA(diff, order=(0,0,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	yhat = inverse_difference(history, yhat, months_in_year)
	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)
print(residuals.describe())
# plot
pyplot.figure()
pyplot.subplot(211)
plot_acf(residuals, ax=pyplot.gca())
pyplot.subplot(212)
plot_pacf(residuals, ax=pyplot.gca())
pyplot.show()

        结果表明，时间序列中存在的少量自相关已被模型捕获。

5、模型验证

        在开发模型并选择最终模型后，必须对其进行验证和最终确定。

        验证是该过程的可选部分，但它提供“最后检查”以确保我们没有欺骗或误导自己。

（1）完成模型

        最终确定模型涉及在整个数据集上拟合 ARIMA 模型，在本例中是在整个数据集的转换版本上。

        拟合后，可以将模型保存到文件中以备后用。

        下面的示例在数据集上训练 ARIMA(0,0,1) 模型并将整个拟合对象和偏差保存到文件中。

        下面的示例将拟合模型以正确的状态保存到文件中，以便以后可以成功加载。

# save finalized model
from pandas import read_csv
from statsmodels.tsa.arima.model import ARIMA
import numpy
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# load data
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
# prepare data
X = series.values
X = X.astype('float32')
# difference data
months_in_year = 12
diff = difference(X, months_in_year)
# fit model
model = ARIMA(diff, order=(0,0,1))
model_fit = model.fit()
# bias constant, could be calculated from in-sample mean residual
bias = 165.904728
# save model
model_fit.save('model.pkl')
numpy.save('model_bias.npy', [bias])

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

        model.pkl这是来自ARIMA.fit()调用的 ARIMAResult 对象。这包括拟合模型时返回的系数和所有其他内部数据。

        model_bias.npy这是存储为单行单列 NumPy 数组的偏差值。

（2）进行预测

        下面的示例加载模型，对下一个时间步进行预测，然后打印预测。

# load finalized model and make a prediction
from pandas import read_csv
from statsmodels.tsa.arima.model import ARIMAResults
import numpy
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
series = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
months_in_year = 12
model_fit = ARIMAResults.load('model.pkl')
bias = numpy.load('model_bias.npy')
yhat = float(model_fit.forecast()[0])
yhat = bias + inverse_difference(series.values, yhat, months_in_year)
print('Predicted: %.3f' % yhat)

        运行该示例会打印出大约 6794 的预测。

Predicted: 6794.773

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

（3）验证模型

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

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

        我们现在可以加载这个validation.csv文件并使用它来查看我们的模型在“未见”数据上的实际效果。

        我们可以通过两种方式进行：

        加载模型并使用它来预测未来 12 个月。最初一两个月之后的预测将很快开始降低技能。
        加载模型并以滚动预测的方式使用它，为每个时间步更新变换和模型。这是首选方法，因为这是在实践中使用此模型的方式，因为它可以实现最佳性能。
        与前几节中的模型评估一样，我们将以滚动预测的方式进行预测。这意味着我们将跨越验证数据集中的提前期，并将观察结果作为对历史的更新。

# load and evaluate the finalized model on the validation dataset
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima.model import ARIMAResults
from sklearn.metrics import mean_squared_error
from math import sqrt
import numpy
 
# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff
 
# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]
 
# load and prepare datasets
dataset = read_csv('dataset.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
X = dataset.values.astype('float32')
history = [x for x in X]
months_in_year = 12
validation = read_csv('validation.csv', header=None, index_col=0, parse_dates=True, squeeze=True)
y = validation.values.astype('float32')
# load model
model_fit = ARIMAResults.load('model.pkl')
bias = numpy.load('model_bias.npy')
# make first prediction
predictions = list()
yhat = float(model_fit.forecast()[0])
yhat = bias + inverse_difference(history, yhat, months_in_year)
predictions.append(yhat)
history.append(y[0])
print('>Predicted=%.3f, Expected=%.3f' % (yhat, y[0]))
# rolling forecasts
for i in range(1, len(y)):
	# difference data
	months_in_year = 12
	diff = difference(history, months_in_year)
	# predict
	model = ARIMA(diff, order=(0,0,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	yhat = bias + inverse_difference(history, yhat, months_in_year)
	predictions.append(yhat)
	# observation
	obs = y[i]
	history.append(obs)
	print('>Predicted=%.3f, Expected=%.3f' % (yhat, obs))
# report performance
rmse = sqrt(mean_squared_error(y, predictions))
print('RMSE: %.3f' % rmse)
pyplot.plot(y)
pyplot.plot(predictions, color='red')
pyplot.show()

        运行该示例会打印验证数据集中时间步长的每个预测值和期望值。

        验证期的最终 RMSE 预计为 361.110 百万销售额。

        这比预期的每月销售额略多于 9.24 亿的错误要好得多。

>Predicted=6794.773, Expected=6981
>Predicted=10101.763, Expected=9851
>Predicted=13219.067, Expected=12670
>Predicted=3996.535, Expected=4348
>Predicted=3465.934, Expected=3564
>Predicted=4522.683, Expected=4577
>Predicted=4901.336, Expected=4788
>Predicted=5190.094, Expected=4618
>Predicted=4930.190, Expected=5312
>Predicted=4944.785, Expected=4298
>Predicted=1699.409, Expected=1413
>Predicted=6085.324, Expected=5877
RMSE: 361.110

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

        在这个规模上，12 个月的预测销售数据看起来很棒。