时间序列分析包括分析时间序列数据的方法，以便提取有意义的统计数据和数据的其他特征。时间序列预测是利用一个模型，根据先前观察到的值来预测未来的值。

时间序列广泛应用于非平稳数据，如经济、天气、股票价格、零售销售等。我们将展示预测零售销售时间序列的不同方法。

数据

我们正在使用可从此处下载的Superstore销售数据（）。

import warnings

import itertools

import numpy as np

import matplotlib.pyplot as plt

warnings.filterwarnings("ignore")

plt.style.use('fivethirtyeight')

import pandas as pd

import statsmodels.api as sm

import matplotlib

matplotlib.rcParams['axes.labelsize'] = 14

matplotlib.rcParams['xtick.labelsize'] = 12

matplotlib.rcParams['ytick.labelsize'] = 12

matplotlib.rcParams['text.color'] = 'k'

超市销售数据中有几个类别，我们从家具销售的时间序列分析和预测开始。

df = pd.read_excel（“Superstore.xls”）

furniture = df.loc [df ['Category'] =='Furniture']

我们有4年的家具销售数据。

furniture['Order Date'].min(), furniture['Order Date'].max()

Timestamp(‘2014–01–06 00:00:00’), Timestamp(‘2017–12–30 00:00:00’)

数据预处理

此步骤包括删除我们不需要的列，检查缺失值，按日期汇总销售额等。

cols = ['Row ID', 'Order ID', 'Ship Date', 'Ship Mode', 'Customer ID', 'Customer Name', 'Segment', 'Country', 'City', 'State', 'Postal Code', 'Region', 'Product ID', 'Category', 'Sub-Category', 'Product Name', 'Quantity', 'Discount', 'Profit']

furniture.drop(cols, axis=1, inplace=True)

furniture = furniture.sort_values('Order Date')

furniture.isnull().sum()

图1

furniture = furniture.groupby('Order Date')['Sales'].sum().reset_index()

使用时间序列数据建立索引

furniture = furniture.set_index（'Order Date'）

furniture.index

图2

我们当前的日期时间数据可能很难处理，因此，我们将使用该月的平均每日销售价值，而我们使用每月的开头作为时间戳。

y = furniture['Sales'].resample('MS').mean()

快速浏览2017年家具销售数据。

Y [ '2017'：]

图3

可视化家具销售时间序列数据

y.plot（figsize =（15,6））

plt.show（）

图4

绘制数据时会出现一些可区分的模式。时间序列具有季节性模式，例如年初销售额总是较低，年末销售额较高。任何一年都有一个上升的趋势，在一年中有几个低月。

我们还可以使用称为时间序列分解的方法来可视化我们的数据，该方法允许我们将时间序列分解为三个不同的组件：趋势，季节性和噪声

from pylab import rcParams

rcParams['figure.figsize'] = 18, 8

decomposition = sm.tsa.seasonal_decompose(y, model='additive')

fig = decomposition.plot()

plt.show()

图5

上面的图形清楚地表明，家具的销售不稳定，以及明显的季节性

ARIMA的时间序列预测

我们将应用最常用的时间序列预测方法之一，即ARIMA，它代表自回归综合移动平均线。

ARIMA模型用ARIMA(p, d, q)表示。这三个参数考虑了数据中的季节性，趋势和噪音：

p = d = q = range(0, 2)

pdq = list(itertools.product(p, d, q))

seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

print('Examples of parameter combinations for Seasonal ARIMA...')

print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))

print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))

print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))

print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

图6

这一步是我们家具销售ARIMA时间序列模型的参数选择。我们的目标是使用“grid search”来找到最优的参数集，从而为我们的模型提供最佳的性能。

for param in pdq:

for param_seasonal in seasonal_pdq:

try:

mod = sm.tsa.statespace.SARIMAX(y,

order=param,

seasonal_order=param_seasonal,

enforce_stationarity=False,

enforce_invertibility=False)

results = mod.fit()

print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))

except:

continue

图7

上述输出结果表明，SARIMAX(1,1,1)x(1,1,0,12)的AIC最小值为297.78。因此，我们认为这是最优的选择。

Fitting the ARIMA model

mod = sm.tsa.statespace.SARIMAX(y,

order=(1, 1, 1),

seasonal_order=(1, 1, 0, 12),

enforce_stationarity=False,

enforce_invertibility=False)

results = mod.fit()

print(results.summary().tables[1])

图8

我们应该始终运行模型诊断以调查任何异常行为。

results.plot_diagnostics(figsize=(16, 8))

plt.show()

图9

然而，它并不完美，我们的模型诊断表明模型残差接近正态分布。

验证预测

为了帮助我们了解预测的准确性，我们将预测销售额与时间序列的实际销售额进行比较，并将预测从2017-07-01开始到数据结束。

pred = results.get_prediction(start=pd.to_datetime('2017-01-01'), dynamic=False)

pred_ci = pred.conf_int()

ax = y['2014':].plot(label='observed')

pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7, figsize=(14, 7))

ax.fill_between(pred_ci.index,

pred_ci.iloc[:, 0],

pred_ci.iloc[:, 1], color='k', alpha=.2)

ax.set_xlabel('Date')

ax.set_ylabel('Furniture Sales')

plt.legend()

plt.show()

图10

线图显示了观测值与滚动预测预测的对比。总体而言，我们的预测与真实值非常吻合，显示出从年初开始的上升趋势，并捕捉到了年底的季节性。

y_forecasted = pred.predicted_mean

y_truth = y['2017-01-01':]

mse = ((y_forecasted - y_truth) ** 2).mean()

print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

The Mean Squared Error of our forecasts is 22993.58

print('The Root Mean Squared Error of our forecasts is {}'.format(round(np.sqrt(mse), 2)))

The Root Mean Squared Error of our forecasts is 151.64

在统计学中，估计量的平均平方误差(MSE)测量误差平方和的平均值——也就是说，估计值与估计值之间的平均平方差值。MSE是对估计值质量的度量——它总是非负的，MSE越小，我们就越接近于找到最佳拟合线。

均方根误差(RMSE)告诉我们，我们的模型能够在实际销售的151.64范围内预测测试集中的平均日家具销售量。我们的家具日销售额在400到1200之间。在我看来，到目前为止这是一个很好的模型。

制作和可视化预测

pred_uc = results.get_forecast（steps = 100）

pred_ci = pred_uc.conf_int（）

ax = y.plot（label ='observed'，figsize =（14,7））

pred_uc.predicted_mean.plot（ax = ax，label ='Forecast'）

ax.fill_between（pred_ci.index，

pred_ci.iloc [：， 0]，

pred_ci.iloc [：，1]，color ='k'，alpha = .25）

ax.set_xlabel（'Date'）

ax.set_ylabel（'Furniture Sales'）

plt.legend（）

plt.show（）

图11

我们的模型清楚地反映了家具销售的季节 正如我们对未来的进一步预测，我们自然会对我们的价值观失去信心。这可以通过我们的模型生成的置信区间来反映，随着我们进一步向未来迈进，这些区间会变得更大。

上面对家具的时间序列分析让我对其他类别感到好奇，以及它们如何在时间上相互比较。因此，我们将比较家具和办公室供应商的时间序列。

时间系列家具与办公用品

根据我们的数据，多年来办公用品的销售数量远远超过家具。

furniture = df.loc [df ['Category'] =='Furniture']

office = df.loc [df ['Category'] =='Office Supplies']

furniture.shape，office.shape

（（2121,21），（6026,21））

数据探索

我们将比较两个类别在同一时期的销售情况。这意味着将两个数据框合并为一个，并将这两个类别的时间序列绘制成一个图。

cols = ['Row ID', 'Order ID', 'Ship Date', 'Ship Mode', 'Customer ID', 'Customer Name', 'Segment', 'Country', 'City', 'State', 'Postal Code', 'Region', 'Product ID', 'Category', 'Sub-Category', 'Product Name', 'Quantity', 'Discount', 'Profit']

furniture.drop(cols, axis=1, inplace=True)

office.drop(cols, axis=1, inplace=True)

furniture = furniture.sort_values('Order Date')

office = office.sort_values('Order Date')

furniture = furniture.groupby('Order Date')['Sales'].sum().reset_index()

office = office.groupby('Order Date')['Sales'].sum().reset_index()

furniture = furniture.set_index('Order Date')

office = office.set_index('Order Date')

y_furniture = furniture['Sales'].resample('MS').mean()

y_office = office['Sales'].resample('MS').mean()

furniture = pd.DataFrame({'Order Date':y_furniture.index, 'Sales':y_furniture.values})

office = pd.DataFrame({'Order Date': y_office.index, 'Sales': y_office.values})

store = furniture.merge(office, how='inner', on='Order Date')

store.rename(columns={'Sales_x': 'furniture_sales', 'Sales_y': 'office_sales'}, inplace=True)

store.head()

图12

plt.figure(figsize=(20, 8))

plt.plot(store['Order Date'], store['furniture_sales'], 'b-', label = 'furniture')

plt.plot(store['Order Date'], store['office_sales'], 'r-', label = 'office supplies')

plt.xlabel('Date'); plt.ylabel('Sales'); plt.title('Sales of Furniture and Office Supplies')

plt.legend();

图13

我们观察到家具和办公用品的销售情况类似。今年年初是这两个类别的淡季。对于办公用品来说，夏天时间似乎也很安静。此外，在大多数月份，家具的平均日销售额高于办公用品。这是可以理解的，因为家具的价值应远高于办公用品的价值。有时，办公用品平均每天销售家具。让我们来看看办公用品第一次销售超过家具的时间。

first_date = store.ix[np.min(list(np.where(store['office_sales'] > store['furniture_sales'])[0])), 'Order Date']

print("Office supplies first time produced higher sales than furniture is {}.".format(first_date.date()))

Office supplies first time produced higher sales than furniture is 2014–07–01.

Prophet时间序列建模

预测工具Prophet于2017年由Facebook发布，旨在分析在不同时间尺度（如年度，每周和每日）显示模式的时间序列。它还具有高级功能，可以对假期对时间序列的影响进行建模并实现自定义变更点。因此，我们正在使用Prophet来建立和运行模型。

from fbprophet import Prophet

furniture = furniture.rename(columns={'Order Date': 'ds', 'Sales': 'y'})

furniture_model = Prophet(interval_width=0.95)

furniture_model.fit(furniture)

office = office.rename(columns={'Order Date': 'ds', 'Sales': 'y'})

office_model = Prophet(interval_width=0.95)

office_model.fit(office)

furniture_forecast = furniture_model.make_future_dataframe(periods=36, freq='MS')

furniture_forecast = furniture_model.predict(furniture_forecast)

office_forecast = office_model.make_future_dataframe(periods=36, freq='MS')

office_forecast = office_model.predict(office_forecast)

plt.figure(figsize=(18, 6))

furniture_model.plot(furniture_forecast, xlabel = 'Date', ylabel = 'Sales')

plt.title('Furniture Sales');

图14

plt.figure（figsize =（18,6））

office_model.plot（office_forecast，xlabel ='Date'，ylabel ='Sales'）

plt.title（'Office Supplies Sales'）;

图15

比较预测

我们已经对未来的这两个类别进行了三年的预测。我们现在将他们联合起来比较他们未来的预测。

furniture_names = ['furniture_%s' % column for column in furniture_forecast.columns]

office_names = ['office_%s' % column for column in office_forecast.columns]

merge_furniture_forecast = furniture_forecast.copy()

merge_office_forecast = office_forecast.copy()

merge_furniture_forecast.columns = furniture_names

merge_office_forecast.columns = office_names

forecast = pd.merge(merge_furniture_forecast, merge_office_forecast, how = 'inner', left_on = 'furniture_ds', right_on = 'office_ds')

forecast = forecast.rename(columns={'furniture_ds': 'Date'}).drop('office_ds', axis=1)

forecast.head()

图16

趋势和预测可视化

plt.figure(figsize=(10, 7))

plt.plot(forecast['Date'], forecast['furniture_trend'], 'b-')

plt.plot(forecast['Date'], forecast['office_trend'], 'r-')

plt.legend(); plt.xlabel('Date'); plt.ylabel('Sales')

plt.title('Furniture vs. Office Supplies Sales Trend');

图17

plt.figure(figsize=(10, 7))

plt.plot(forecast['Date'], forecast['furniture_yhat'], 'b-')

plt.plot(forecast['Date'], forecast['office_yhat'], 'r-')

plt.legend(); plt.xlabel('Date'); plt.ylabel('Sales')

plt.title('Furniture vs. Office Supplies Estimate');

图18

趋势和模式

现在，我们可以使用Prophet模型来检查数据中这两个类别的不同趋势。

furniture_model.plot_components（furniture_forecast）;

图19

office_model.plot_components（office_forecast）;

图20

很高兴看到家具和办公用品的销售额随着时间的推移呈线性增长，尽管办公用品的增长似乎略微增强。

家具最糟糕的一个月是四月，办公用品最糟糕的月份是二月份。家具的最佳月份是十二月，办公用品的最佳月份是十一月。

从现在开始我们可以探索许多时间序列分析，例如不确定边界的预测，变化点和异常检测，与外部数据源的预测时间序列。