学习预测时序数据，如有侵权，请联系删除。
主要参考：
A comprehensive beginner’s guide to create a Time Series Forecast (with Codes in Python)
目的：
- 翻译学习文章编写方式
- 学习解决时间序列问题的基本步骤
介绍
时间序列（从现在开始称为TS）被认为是数据科学领域中鲜为人知的技能之一。我为自己设定了解决时间序列问题的基本步骤，并在此与大家分享。这些肯定会帮助您在未来的任何项目中获得一个不错的模型！
在阅读本文之前，我强烈建议您阅读A Complete Tutorial on Time Series Modeling in R并参加free Time Series Forecasting course。它侧重于基本概念，我将专注于使用这些概念来解决端到端问题以及Python中的代码。R中有许多时间序列的资源，但是很少有Python的资源，所以我将在本文中使用Python。
我们的旅程将经历以下步骤：
- 什么使时间序列特别？
- 在Pandas中加载和处理时间序列
- 如何检查时间序列的平稳性？
- 如何制作时间序列固定？
- 预测时间序列
1.什么使时间序列特别？
顾名思义，TS是以恒定时间间隔收集的数据点的集合。对这些进行分析以确定长期趋势，以便预测未来或执行其他形式的分析。但是什么使TS不同于说常规回归问题？有两件事：
- 这是时间依赖的。因此，在这种情况下，观察是独立的线性回归模型的基本假设不成立。
- 随着趋势的增加或减少，大多数TS具有某种形式的季节性趋势，即特定时间范围的特定变化。例如，如果你看到一件羊毛外套随着时间的推移而销售，那么你将在冬季找到更高的销售额。
由于TS的固有属性，分析它涉及各种步骤。这些将在下面详细讨论。让我们从在Python中加载TS对象开始。我们将使用流行的AirPassengers数据集，可在此处
下载。
请注意，本文的目的是让您熟悉TS一般使用的各种技术。这里考虑的例子仅用于说明，我将重点关注广泛的主题，而不是做出非常准确的预测。
2.在Pandas中加载和处理时间序列
Pandas有专门的库来处理TS对象，特别是datatime64 [ns]类，它存储时间信息并允许我们快速执行某些操作。让我们从启动所需的库开始：
import pandas as pd import numpy as np import matplotlib.pylab as plt %matplotlib inline from matplotlib.pylab import rcParams rcParams['figure.figsize'] = 15, 6
现在，我们可以加载数据集并查看列的一些初始行和数据类型：
data = pd.read_csv('AirPassengers.csv') print data.head() print ' Data Types:' print data.dtypes
该数据包含特定月份和该月旅行的乘客数量。但是这仍然不是作为TS对象读取的，因为数据类型是'object'和'int'。为了将数据作为时间序列读取，我们必须将特殊参数传递给read_csv命令：
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m') data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse) print data.head()
让我们逐一理解这些论点：
- parse_dates：指定包含日期时间信息的列。如上所述，列名称为“Month”。
- index_col：将Pandas用于TS数据背后的一个关键思想是索引必须是描述日期时间信息的变量。所以这个论点告诉pandas使用'Month'列作为索引。
- date_parser：这指定一个将输入字符串转换为datetime变量的函数。默认情况下，Pandas以“YYYY-MM-DD HH：MM：SS”格式读取数据。如果数据不是这种格式，则必须手动定义格式。类似于此处定义的数据分析函数的东西可用于此目的。
现在我们可以看到数据有时间对象作为索引，＃Passengers作为列。我们可以使用以下命令交叉检查索引的数据类型：
data.index
注意dtype ='datetime [ns]'确认它是一个datetime对象。作为个人偏好，我会将列转换为Series对象，以防止每次使用TS时引用列名。请随意使用作为数据帧，这对您更有效。
ts = data[‘#Passengers’] ts.head(10)
在进一步讨论之前，我将讨论TS数据的一些索引技术。让我们首先选择Series对象中的特定值。这可以通过以下两种方式完成：
#1. Specific the index as a string constant: ts['1949-01-01'] #2. Import the datetime library and use 'datetime' function: from datetime import datetime ts[datetime(1949,1,1)]
两者都将返回值'112'，这也可以从先前的输出中确认。假设我们想要所有数据到1949年5月。这可以通过两种方式完成：
#1. Specify the entire range: ts['1949-01-01':'1949-05-01'] #2. Use ':' if one of the indices is at ends: ts[:'1949-05-01']
两者都会产生以下输出：
这里有两件事需要注意：
- 与数字索引不同，此处包含结束索引。例如，如果我们将列表索引为[：5]，那么它将返回索引处的值 - [0,1,2,3,4]。但这里的指数'1949-05-01'已包含在输出中。
- 索引必须按范围排序才能工作。如果你随机乱改索引，这将无法正常工作。
考虑另一个需要1949年所有值的实例。这可以通过以下方式完成：
ts['1949']
月份部分被省略了。同样，如果您在特定月份的所有日期，可以省略日期部分。
现在，让我们转到分析TS。
3.如何检查时间序列的平稳性？
如果TS的统计特性（例如均值，方差）随时间保持不变，则称TS是静止的。但为什么这很重要？大多数TS模型的工作假设TS是静止的。直觉上，我们可以认为，如果TS随着时间的推移具有特定的行为，那么将来它很可能会遵循相同的行为。此外，与非平稳序列相比，与静止序列相关的理论更成熟，更容易实现。
使用非常严格的标准来定义平稳性。然而，出于实际目的，如果系列随时间具有恒定的统计特性，即我们可以假设该系列是静止的。下列：
- 恒定的均值
- 恒定的方差
- 一种不依赖于时间的自协方差
我将跳过细节，因为它在本文中已经非常明确地定义了。让我们进入测试平稳性的方法。首要的是简单绘制数据并进行视觉分析。可以使用以下命令绘制数据：
plt.plot(ts)
很明显，数据总体上呈增长趋势，并伴有一些季节性变化。然而，可能并不总是可以做出这样的视觉推断（我们稍后会看到这样的情况）。因此，更正式地说，我们可以使用以下方法检查平稳性：
- 绘制滚动统计：我们可以绘制移动平均值或移动方差，并查看它是否随时间变化。通过移动平均值/方差，我的意思是在任何时刻“t”，我们将采用去年的平均值/方差，即过去12个月。但这又是一种视觉技术。
- Dickey-Fuller测试：这是检查平稳性的统计测试之一。这里的零假设是TS是非平稳的。测试结果包括测试统计和差异置信水平的一些关键值。如果'测试统计'小于'临界值'，我们可以拒绝原假设并说该系列是静止的。
回到检查平稳性，我们将使用滚动统计图和Dickey-Fuller测试结果，因此我定义了一个函数，它将TS作为输入并为我们生成它们。请注意，我绘制了标准偏差而不是方差，以使单位与平均值相似。
from statsmodels.tsa.stattools import adfuller def test_stationarity(timeseries): #Determing rolling statistics #rolmean = pd.rolling_mean(timeseries, window=12) #rolstd = pd.rolling_std(timeseries, window=12) rolmean = timeseries.rolling(12).mean() rolstd = timeseries.rolling(12).std() #Plot rolling statistics: orig = plt.plot(timeseries, color='blue',label='Original') mean = plt.plot(rolmean, color='red', label='Rolling Mean') std = plt.plot(rolstd, color='black', label = 'Rolling Std') plt.legend(loc='best') plt.title('Rolling Mean & Standard Deviation') plt.show(block=False) #Perform Dickey-Fuller test: print 'Results of Dickey-Fuller Test:' dftest = adfuller(timeseries, autolag='AIC') dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used']) for key,value in dftest[4].items(): dfoutput['Critical Value (%s)'%key] = value print dfoutput
test_stationarity(ts)
虽然标准偏差的变化很小，但平均值随着时间的推移而明显增加，这不是一个固定的系列。此外，测试统计数据远远超过临界值。请注意，应比较有符号值，而不是绝对值。
4.如何制作时间序列固定？
让我们了解什么使TS非静止。 TS的非固定性背后有两个主要原因：
- 趋势 - 随时间变化的平均值。例如，在这种情况下，我们看到平均而言，乘客数量随着时间的推移而增长。
- 季节性 - 特定时间范围的变化。例如，由于工资增加或节日，人们可能倾向于在特定月份购买汽车。
基本原则是模拟或估计系列中的趋势和季节性，并从系列中删除那些以获得固定系列。然后可以在这个系列中实施统计预测技术。最后一步是通过应用趋势和季节性约束将预测值转换为原始比例。
估计和消除趋势
减少趋势的第一个技巧之一可能是转型。例如，在这种情况下，我们可以清楚地看到存在显着的积极趋势。因此，我们可以应用比较小值更多地惩罚更高值的转换。这些可以采用日志，平方根，立方根等。为简单起见，请在此处进行日志转换：
ts_log = np.log(ts)
plt.plot(ts_log)
在这种简单的情况下，很容易看到数据的前向趋势。但它在噪音的存在下不是很直观。因此，我们可以使用一些技术来估计或模拟这种趋势，然后将其从系列中删除。可以有很多方法，最常用的一些方法是：
- 聚合 - 取每月/每周平均值等时间段的平均值
- 平滑 - 采用滚动平均值
- 多项式拟合 - 拟合回归模型
平滑是指采用滚动估计，即考虑过去的几个实例。可以有各种方式，但我将在这里讨论其中的两个。
移动平均线
在这种方法中，我们根据时间序列的频率取“k”个连续值的平均值。在这里我们可以取过去1年的平均值，即最后12个值。 Pandas具有为确定滚动统计数据而定义的特定功能。
# moving_avg = pd.rolling_mean(ts_log,12)
moving_avg = ts_log.rolling(12).mean() plt.plot(ts_log) plt.plot(moving_avg, color='red')
红线表示滚动平均值。让我们从原始系列中减去它。请注意，由于我们取最后12个值的平均值，因此未对前11个值定义滚动平均值。这可以被视为：
ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.head(12)
注意前11是NaN。让我们删除这些NaN值并检查图以测试平稳性。
ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)
这看起来像一个更好的系列。滚动值似乎略有不同，但没有具体趋势。此外，测试统计量小于5％的临界值，因此我们可以95％的置信度说这是一个固定的序列。
然而，这种特定方法的缺点是必须严格定义时间段。在这种情况下，我们可以采取年平均值，但在复杂的情况下，如预测股票价格，很难得出一个数字。因此，我们采用“加权移动平均线”，其中更近期的值被赋予更高的权重。可以有许多分配权重的技术。一种流行的是指数加权移动平均值，其中权重被分配给具有衰减因子的所有先前值。在此查找详情。这可以在Pandas中实现：
#expwighted_avg = pd.ewma(ts_log, halflife=12)
expwighted_avg = ts_log.ewm(halflife=12, min_periods=0,adjust=True,ignore_na=False).mean() plt.plot(ts_log) plt.plot(expwighted_avg, color='red')
注意，这里参数'halflife'用于定义指数衰减量。
ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)
该TS的平均值和标准偏差的幅度变化甚至更小。此外，检验统计量小于1％临界值，这比前一种情况更好。
消除趋势和季节性
之前讨论的简单趋势减少技术在所有情况下都不起作用，尤其是具有高季节性的技术。让我们讨论两种消除趋势和季节性的方法：
- 差异 - 将差异与特定时滞相区分
- 分解 - 对趋势和季节性进行建模并将其从模型中移除。
差分
在这种技术中，我们将特定时刻的观察值与前一时刻的观察值进行区分。这主要用于改善平稳性。一阶差分可以在Pandas中完成：
ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
我们可以看到均值和标准变量随时间变化很小。此外，Dickey-Fuller检验统计量小于10％临界值，因此TS静止且置信度为90％。
分解
在这种方法中，趋势和季节性都是单独建模的，并返回系列的其余部分。我将跳过统计数据并得出结果：
from statsmodels.tsa.seasonal import seasonal_decompose decomposition = seasonal_decompose(ts_log) trend = decomposition.trend seasonal = decomposition.seasonal residual = decomposition.resid plt.subplot(411) plt.plot(ts_log, label='Original') plt.legend(loc='best') plt.subplot(412) plt.plot(trend, label='Trend') plt.legend(loc='best') plt.subplot(413) plt.plot(seasonal,label='Seasonality') plt.legend(loc='best') plt.subplot(414) plt.plot(residual, label='Residuals') plt.legend(loc='best') plt.tight_layout()
在这里我们可以看到趋势，季节性与数据分离，我们可以对残差进行建模。让我们检查残差的平稳性。
ts_log_decompose = residual ts_log_decompose.dropna(inplace=True) test_stationarity(ts_log_decompose)
Dickey-Fuller检验统计量显着低于1％临界值。所以这个TS非常接近静止。您也可以尝试高级分解技术，这可以产生更好的结果。此外，您应该注意，在这种情况下，将残差转换为未来数据的原始值并不是非常直观。
5.预测时间序列
我们看到了不同的技术，所有这些技术都使得TS固定不动。让我们在差分后在TS上制作模型，因为它是一种非常流行的技术。此外，在这种情况下，相对容易将噪声和季节性添加回预测残差。执行趋势和季节性估计技术后，可能存在两种情况：
- 严格固定的系列，不依赖于数值。这是一个简单的情况，我们可以将残差建模为白噪声。但这种情况非常罕见。
- 一系列重要的价值观依赖。在这种情况下，我们需要使用一些统计模型，如ARIMA来预测数据。
让我简单介绍一下ARIMA。我不会详细介绍技术细节，但如果您希望更有效地应用它们，您应该详细了解这些概念。ARIMA代表自动回归综合移动平均线。 ARIMA对静止时间序列的预测只不过是一个线性（如线性回归）方程。预测变量取决于ARIMA模型的参数（p，d，q）：
- AR（自回归）项的数量（p）：AR项只是因变量的滞后。例如，如果p是5，则x（t）的预测值将是x（t-1）... .x（t-5）。
- MA（移动平均线）项的数量（q）：MA项是预测方程中的滞后预测误差。例如，如果q为5，则x（t）的预测值将为e（t-1）... .e（t-5），其中e（i）是第i个瞬时的移动平均值与实际值之间的差值。
- 差异数量（d）：这些是非季节性差异的数量，即在这种情况下我们采用了第一阶差异。因此，我们可以传递该变量并将d = 0或传递原始变量并放入d = 1。两者都会产生相同的结果。
这里一个重要的问题是如何确定'p'和'q'的值。我们使用两个图来确定这些数字。让我们先讨论一下。
- 自相关函数（ACF）：它是TS与其自身滞后版本之间相关性的度量。例如，在滞后5处，ACF将在时刻't1'...'t2'的系列与在时刻't1-5'...'t2-5'（t1-5和t2是终点）的系列进行比较。
- 部分自相关函数（PACF）：它测量TS与其自身的滞后版本之间的相关性，但是在消除了已经通过干预比较解释的变化之后。例如，在滞后5处，它将检查相关性，但是除去已经由滞后1到4解释的效果。
差分后的TS的ACF和PACF图可以绘制为：
#ACF and PACF plots: from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20) lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF: plt.subplot(121) plt.plot(lag_acf) plt.axhline(y=0,linestyle='--',color='gray') plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray') plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray') plt.title('Autocorrelation Function')
#Plot PACF: plt.subplot(122) plt.plot(lag_pacf) plt.axhline(y=0,linestyle='--',color='gray') plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray') plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray') plt.title('Partial Autocorrelation Function') plt.tight_layout()
在这个图中，0两边的两条虚线是置信区间。这些可用于确定'p'和'q'值为：
- p - PACF图表首次超过置信区间的滞后值。如果你仔细注意，在这种情况下p = 2。
- q - ACF图表第一次超过置信区间的滞后值。如果你仔细注意，在这种情况下q = 2。
现在，让我们考虑个人和组合效果制作3种不同的ARIMA模型。我还将为每个打印RSS。请注意，这里的RSS是残差值，而不是实际序列。
ARModel
from statsmodels.tsa.arima_model import ARIMA model = ARIMA(ts_log, order=(2, 1, 0)) results_AR = model.fit(disp=-1) plt.plot(ts_log_diff) plt.plot(results_AR.fittedvalues, color='red') plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))
MRModel
model = ARIMA(ts_log, order=(0, 1, 2)) results_MA = model.fit(disp=-1) plt.plot(ts_log_diff) plt.plot(results_MA.fittedvalues, color='red') plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))
Combined Model
model = ARIMA(ts_log, order=(2, 1, 2)) results_ARIMA = model.fit(disp=-1) plt.plot(ts_log_diff) plt.plot(results_ARIMA.fittedvalues, color='red') plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))
在这里我们可以看到AR和MA模型具有几乎相同的RSS，但组合明显更好。现在，我们留下最后一步，即将这些值恢复到原始比例。
把它恢复到原始规模
由于组合模型给出了最佳结果，因此我们可以将其缩放回原始值并查看其在那里的表现。第一步是将预测结果存储为单独的系列并观察它。
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True) print predictions_ARIMA_diff.head()
请注意，这些是从'1949-02-01'而不是第一个月开始的。为什么？这是因为我们将滞后加1，并且第一个元素在它减去之前没有任何东西。将差分转换为对数比例的方法是将这些差异连续添加到基数.一种简单的方法是首先确定索引处的累积和，然后将其添加到基数。累计金额可以是：
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum() print predictions_ARIMA_diff_cumsum.head()
predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index) predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0) predictions_ARIMA_log.head()
这里第一个元素是基数本身，并且从那里累加的值。最后一步是取指数并与原始系列进行比较.
predictions_ARIMA = np.exp(predictions_ARIMA_log) plt.plot(ts) plt.plot(predictions_ARIMA) plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))
