在本文中，您将学习执行时间序列分析的基本步骤以及趋势，平稳性，移动平均等概念。还将探索指数平滑方法，并学习如何拟合ARIMA模型非平稳数据。

定义

时间序列是按时间排序（或索引）的数据序列。它是离散的，每个点之间的间隔是不变的。

序列的属性和类型

趋势 ：数据的长期增加或减少。这可以看作是粗略地通过数据的斜率(不一定是线性的)。

季节性 ：当一个时间序列受到季节因素(日小时、周、月、年等)的影响时，就称为季节序列。季节性可以用固定频率的良好周期性模式来观察。

循环 ：当数据表现出不是固定频率的上升和下降时，就会发生循环。这些波动通常是由经济条件造成的，通常与“商业周期”有关。这些波动的持续时间通常至少为2年。

上图：具有趋势的序列 。 下：序列表现出 季节性 行为

上： 循环序列。 下： 随机游走

平稳性

在进一步深入分析之前，我们的序列必须是平稳的。

平稳性是表现出恒定统计特性(平均值、方差、自相关等)的特性。如果时间序列的平均值随时间增加，那么它就不是平稳的。

用于平稳数据的变换：

• De-trending ：我们去掉了这个序列的潜在趋势。根据数据的性质，可以用几种方法来做到这一点:
• –Indexed data:以货币计量的数据与价格指数或与通货膨胀有关。因此，用这个指数(即平减指数)元素来划分级数是消除数据趋势的解决方案。
• –Non-indexed data：是否有必要估计趋势是恒定的，线性的还是指数的。前两种情况很容易，对于最后一种情况，有必要估算增长率（通货膨胀或通货紧缩）并采用与索引数据相同的方法。
• Differencing :季节或周期模式可以通过减去周期值来消除。如果数据是12个月的季节性数据，那么用12个滞后差值序列减去该序列将得到一个“平坦”的序列
• Logging ：如果趋势中的复合率不是由价格指数引起的（即序列不是以货币计量），Logging 可以帮助线性化具有指数趋势的序列（回想一下log(exp(x)) = x）。与通货紧缩不同，它不会消除任何最终趋势。

检查平稳性

绘制滚动统计数据

绘制滚动平均值和方差是视觉检查我们的序列的第一个好方法。如果滚动统计数据显示出明显的趋势（向上或向下）并显示变化的方差（增大或减小幅度），那么您可能会得出结论，该序列很可能不是平稳的。Python实现如下：

import requests import pandas as pd import json import matplotlib.pyplot as plt import matplotlib.dates as mdates %matplotlib inline plt.style.use('Solarize_Light2') r = requests.get('https://datamarket.com/api/v1/list.json?ds=22qx') jobj = json.loads(r.text[18:-1]) data = jobj[0]['data'] df = pd.DataFrame(data, columns=['time','data']).set_index('time') train = df.iloc[:-10, :] test = df.iloc[-10:, :] pred = test.copy() df.plot(figsize=(12,3)); plt.title(jobj[0]['title']); df['z_data'] = (df['data'] - df.data.rolling(window=12).mean()) / df.data.rolling(window=12).std() df['zp_data'] = df['z_data'] - df['z_data'].shift(12) def plot_rolling(df): fig, ax = plt.subplots(3,figsize=(12, 9)) ax[0].plot(df.index, df.data, label='raw data') ax[0].plot(df.data.rolling(window=12).mean(), label="rolling mean"); ax[0].plot(df.data.rolling(window=12).std(), label="rolling std (x10)"); ax[0].legend() ax[1].plot(df.index, df.z_data, label="de-trended data") ax[1].plot(df.z_data.rolling(window=12).mean(), label="rolling mean"); ax[1].plot(df.z_data.rolling(window=12).std(), label="rolling std (x10)"); ax[1].legend() ax[2].plot(df.index, df.zp_data, label="12 lag differenced de-trended data") ax[2].plot(df.zp_data.rolling(window=12).mean(), label="rolling mean"); ax[2].plot(df.zp_data.rolling(window=12).std(), label="rolling std (x10)"); ax[2].legend() plt.tight_layout() fig.autofmt_xdate() 

欢迎大家来到IT世界,在知识的湖畔探索吧!

增强Dickey-Fuller测试

这个测试是用来评估一个时间序列是否是平稳的。在不深入讨论假设检验的细节的情况下，您应该知道这个检验将给出一个名为“检验统计量”的结果，根据这个结果，您可以说，如果时间序列是平稳的，那么它具有不同程度(或百分比)的置信度。Python代码如下：

欢迎大家来到IT世界,在知识的湖畔探索吧!from statsmodels.tsa.stattools import adfuller print(" > Is the data stationary ?") dftest = adfuller(df.data, autolag='AIC') print("Test statistic = {:.3f}".format(dftest[0])) print("P-value = {:.3f}".format(dftest[1])) print("Critical values :") for k, v in dftest[4].items(): print("\t{}: {} - The data is {} stationary with {}% confidence".format(k, v, "not" if v<dftest[0] else "", 100-int(k[:-1]))) print("\n > Is the de-trended data stationary ?") dftest = adfuller(df.z_data.dropna(), autolag='AIC') print("Test statistic = {:.3f}".format(dftest[0])) print("P-value = {:.3f}".format(dftest[1])) print("Critical values :") for k, v in dftest[4].items(): print("\t{}: {} - The data is {} stationary with {}% confidence".format(k, v, "not" if v<dftest[0] else "", 100-int(k[:-1]))) print("\n > Is the 12-lag differenced de-trended data stationary ?") dftest = adfuller(df.zp_data.dropna(), autolag='AIC') print("Test statistic = {:.3f}".format(dftest[0])) print("P-value = {:.3f}".format(dftest[1])) print("Critical values :") for k, v in dftest[4].items(): print("\t{}: {} - The data is {} stationary with {}% confidence".format(k, v, "not" if v<dftest[0] else "", 100-int(k[:-1]))) 

KPSS

KPSS（Kwiatkowski-Phillips-Schmidt-Shin）测试零假设，即该序列是趋势平稳的。换句话说，如果检验统计量高于X％置信度阈值，这意味着我们拒绝这一假设，并且该序列在X％置信度下不是趋势平稳的。低于阈值的检验统计量将使我们接受这一假设，并得出结论该序列是趋势平稳的。

自相关图（ACF和PACF）

自相关(ACF)图表示具有滞后的序列的自相关。

偏自相关(PACF)图表示的是一个序列与自身滞后之间的相关性，而所有低阶滞后的相关性都无法解释这种相关性。

理想情况下，我们希望序列与其自身滞后之间没有相关性。从图形上讲，我们希望所有峰值都落在蓝色区域。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf fig, ax = plt.subplots(2, figsize=(12,6)) ax[0] = plot_acf(df.z_data.dropna(), ax=ax[0], lags=20) ax[1] = plot_pacf(df.z_data.dropna(), ax=ax[1], lags=20) 

我们可以看到，蓝色区域上方有几个峰值，这意味着在滞后1,2,3和4处存在相关性。

选择模型

指数平滑方法适用于非平稳数据（即具有趋势和季节性数据的数据）。

ARIMA模型应仅用于平稳数据。因此，应该删除数据的趋势（via deflating or logging），然后查看序列。

非季节性平滑或趋势拟合

简单指数平滑和线性指数平滑

什么时候用？

• 少数数据点
• 不规则数据

欢迎大家来到IT世界,在知识的湖畔探索吧!import requests import pandas as pd import json import matplotlib.pyplot as plt import matplotlib.dates as mdates from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt import numpy as np %matplotlib inline plt.style.use('Solarize_Light2') r = requests.get('https://datamarket.com/api/v1/list.json?ds=22qx') jobj = json.loads(r.text[18:-1]) data = jobj[0]['data'] df = pd.DataFrame(data, columns=['time','data']).set_index('time') train = df.iloc[100:-10, :] test = df.iloc[-10:, :] train.index = pd.to_datetime(train.index) test.index = pd.to_datetime(test.index) pred = test.copy() model = SimpleExpSmoothing(np.asarray(train['data'])) model._index = pd.to_datetime(train.index) fit1 = model.fit() pred1 = fit1.forecast(9) fit2 = model.fit(smoothing_level=.2) pred2 = fit2.forecast(9) fit3 = model.fit(smoothing_level=.5) pred3 = fit3.forecast(9) fig, ax = plt.subplots(figsize=(12, 6)) ax.plot(train.index[150:], train.values[150:]) ax.plot(test.index, test.values, color="gray") for p, f, c in zip((pred1, pred2, pred3),(fit1, fit2, fit3),('#ff7823','#3c763d','c')): ax.plot(train.index[150:], f.fittedvalues[150:], color=c) ax.plot(test.index, p, label="alpha="+str(f.params['smoothing_level'])[:3], color=c) plt.title("Simple Exponential Smoothing") plt.legend(); model = Holt(np.asarray(train['data'])) model._index = pd.to_datetime(train.index) fit1 = model.fit(smoothing_level=.3, smoothing_slope=.05) pred1 = fit1.forecast(9) fit2 = model.fit(optimized=True) pred2 = fit2.forecast(9) fit3 = model.fit(smoothing_level=.3, smoothing_slope=.2) pred3 = fit3.forecast(9) fig, ax = plt.subplots(figsize=(12, 6)) ax.plot(train.index[150:], train.values[150:]) ax.plot(test.index, test.values, color="gray") for p, f, c in zip((pred1, pred2, pred3),(fit1, fit2, fit3),('#ff7823','#3c763d','c')): ax.plot(train.index[150:], f.fittedvalues[150:], color=c) ax.plot(test.index, p, label="alpha="+str(f.params['smoothing_level'])[:4]+", beta="+str(f.params['smoothing_slope'])[:4], color=c) plt.title("Holt's Exponential Smoothing") plt.legend(); 

季节性平滑

Holt线性趋势方法的问题在于，趋势在未来是恒定的，无限增加或减少。对于长期预测视野，这可能会有问题。因此，阻尼趋势法是一种增加阻尼参数的方法，以使趋势在未来收敛到恒定值（它使趋势变平）。参数由φ代替

什么时候用？

• 数据具有趋势并且是季节性的
• 使用乘法版本，除非数据之前已经logged过。在这种情况下，使用加法版本

from statsmodels.tsa.holtwinters import ExponentialSmoothing import numpy as np import requests import pandas as pd import json import matplotlib.pyplot as plt import matplotlib.dates as mdates %matplotlib inline plt.style.use('Solarize_Light2') r = requests.get('https://datamarket.com/api/v1/list.json?ds=22qx') jobj = json.loads(r.text[18:-1]) data = jobj[0]['data'] df = pd.DataFrame(data, columns=['time','data']).set_index('time') train = df.iloc[100:-10, :] test = df.iloc[-10:, :] train.index = pd.to_datetime(train.index) test.index = pd.to_datetime(test.index) pred = test.copy() model = ExponentialSmoothing(np.asarray(train['data']), trend='mul', seasonal=None) model2 = ExponentialSmoothing(np.asarray(train['data']), trend='mul', seasonal=None, damped=True) model._index = pd.to_datetime(train.index) fit1 = model.fit() fit2 = model2.fit() pred1 = fit1.forecast(9) pred2 = fit2.forecast(10) fig, ax = plt.subplots(figsize=(12, 6)) ax.plot(train.index[150:], train.values[150:]) ax.plot(test.index, test.values, color="gray") for p, f, c in zip((pred1, pred2),(fit1, fit2),('#ff7823','#3c763d')): ax.plot(train.index[150:], f.fittedvalues[150:], color=c) print(f.params) print(len(test.index), len(p)) ax.plot(test.index, p, label="alpha="+str(f.params['smoothing_level'])[:4]+", beta="+str(f.params['smoothing_slope'])[:4]+ ", damping="+str(True if f.params['damping_slope']>0 else False), color=c) plt.title("Winter's Exponential Smoothing") plt.legend(); 

AR，ARMA，ARIMA

ARIMA模型(包括ARMA模型、AR模型和MA模型)是预测平稳时间序列的一类通用模型。ARIMA模型由三部分组成:

• 序列滞后值的加权和（自回归（AR）部分）
• 序列滞后预测误差的加权和（移动平均值（MA）部分）
• 时间序列的差分(Integrated (I)部分)

ARIMA模型通常被标注为ARIMA（p，d，q），其中p表示AR部分的顺序，d表示差分的顺序（“ I”部分），q表示MA的顺序。

1）选择差分顺序

拟合ARIMA模型的第一步是确定序列平稳化的差分阶数。为此，我们查看ACF和PACF图，并记住这两条规则:

– 规则1：如果该序列具有大量滞后的正自相关，那么它可能需要更高阶的差分。

 - 规则2：如果lag-1自相关为零或负，或者自相关都很小且无模式，则该序列不需要更高阶的差分。如果lag-1自相关为-0.5或更多负数，则序列可能overdifferenced。

我们首先logging数据，因为原始数据呈指数趋势：

fig, ax = plt.subplots(2, sharex=True, figsize=(12,6)) ax[0].plot(df.data.values); ax[0].set_title("Raw data"); ax[1].plot(np.log(df.data.values)); ax[1].set_title("Logged data (deflated)"); ax[1].set_ylim(0, 15); fig, ax = plt.subplots(2, 2, figsize=(12,6)) first_diff = (np.log(df.data)- np.log(df.data).shift()).dropna() ax[0, 0] = plot_acf(np.log(df.data), ax=ax[0, 0], lags=20, title="ACF - Logged data") ax[1, 0] = plot_pacf(np.log(df.data), ax=ax[1, 0], lags=20, title="PACF - Logged data") ax[0, 1] = plot_acf(first_diff , ax=ax[0, 1], lags=20, title="ACF - Differenced Logged data") ax[1, 1] = plot_pacf(first_diff, ax=ax[1, 1], lags=20, title="PACF - Differenced Logged data") 

logged 序列似乎更平坦，但它是平稳的吗？让我们计算一下KPSS测试来检查：

from statsmodels.tsa.stattools import kpss print(" > Is the data stationary ?") dftest = kpss(np.log(df.data), 'ct') print("Test statistic = {:.3f}".format(dftest[0])) print("P-value = {:.3f}".format(dftest[1])) print("Critical values :") for k, v in dftest[3].items(): print("\t{}: {}".format(k, v)) 

检验统计量高于临界值，我们拒绝零假设，我们的序列不是趋势平稳的。

（ – 要有一个趋势平稳的序列，我们可以考虑线性回归我们的logged序列，并将我们的序列除以回归系数…-）

现在让我们看看log -first-difference数据的ACF图:

“ACF – log data”图显示的是非平稳数据，其特征是峰值处的缓慢线性衰减(参见上面的规则1)。添加一阶差分会在滞后值1处产生一个负峰值。根据规则2，我们不需要再对级数求导了。让我们通过比较a(0,0,0)和a (0,1,0) ARIMA模型来检验我们的结果:

from statsmodels.tsa.arima_model import ARIMA model = ARIMA(np.log(df.data).dropna(), (0, 0, 0)) res_000 = model.fit() print(res_000.summary()) model = ARIMA(np.log(df.data).dropna(), (0, 1, 0)) res_010 = model.fit() print(res_010.summary()) fig, ax = plt.subplots(1, 2, sharey=True, figsize=(12, 6)) ax[0].plot(res_000.resid.values, alpha=0.7, label='variance={:.3f}'.format(np.std(res_000.resid.values))); ax[0].hlines(0, xmin=0, xmax=350, color='r'); ax[0].set_title("ARIMA(0,0,0)"); ax[0].legend(); ax[1].plot(res_010.resid.values, alpha=0.7, label='variance={:.3f}'.format(np.std(res_010.resid.values))); ax[1].hlines(0, xmin=0, xmax=350, color='r'); ax[1].set_title("ARIMA(0,1,0)"); ax[1].legend(); 

ARAYA（0,1,0）的Akaike信息准则（AIC）较低，这意味着该模型的性能优于ARIMA（0,0,0）。让我们看一下残差并检查它们的方差：

2）选择MA order

现在我们知道我们需要在模型中包含一阶差分，我们需要选择移动平均阶。这是通过观察差分级数得到的(因为我们刚刚看到一阶差分级数是平稳的)再次，我们看看我们的ACF和PACF图，记住这个规则:

“如果差分序列ACF的lag-1自相关为负，且/或存在明显的截止，则选择MA阶为1”。

问:为什么选择MA阶为1，而不是更高?因为我们要一步一步来。如果我们在高滞后下观察到自相关，并且通过观察我们的(1,1,0)模型的自相关残差，我们仍然可以观察到这些峰值，我们可以增加MA阶，尽管通常不建议超过2!

注意AIC如何再次下降，以及残差方差如何下降。这是我们（1,1,0）ARIMA表现优于（0,1,0）模型的标志！

3）选择AR order

现在您可能想知道:我们是否必须添加AR项?答案是否定的。事实上，你应该在以下情况下添加AR项:

“如果差分序列PACF的lag-1自相关为负，且/或存在明显的截止，则选择AR阶为1。”。

在我们的例子中，我们观察到ACF和PACF中存在负的lag-1自相关。注意，不同序列的PACF在延迟1和2时显示两个负峰值，这意味着理论上我们可以将AR阶提高到2。

以下是我们通过拟合（1,1,1）ARIMA得到的结果：

增加AR项降低了AIC方差从0.155降到0.145，很好!我们是否应该添加另一个AR项并选择a (2,1,1) ARIMA呢?让我们来看看ARIMA(1,1,1)残差的ACF和PACF:

在图中没有显著的自相关滞后值，我们的模型似乎不需要任何额外的AR项。你可能会指出滞后12的峰值：这可能是季节性的一个标志。我们将在下一部分中使用SARIMA和SARIMAX来查看它。

绘制预测

幸运的是，statsmodels有一个很好的API，允许我们从ARIMA流程绘制预测。我强烈建议您将数据放在带有DateTimeIndex的DataFrame中，因为plot_predict（）方法确实喜欢日期…

在这种情况下，我们将选择12个预测步骤，并将dynamic关键字设置为True：这将强制算法使用自己的预测值作为未来预测的滞后值：

model = ARIMA(np.log(df.data).dropna()[:-12], (1, 1, 1)) res_111 = model.fit() fig, ax = plt.subplots(figsize=(12, 6)) df.index = pd.to_datetime(df.index, format="%Y-%m") np.log(df.data).dropna()[250:].plot(ax=ax); ax.vlines('1992-10', 13, 14.5, linestyle='--', color='r', label='Start of forecast'); # - NOTE from the official documentation : # -- The dynamic keyword affects in-sample prediction. # -- If dynamic is False, then the in-sample lagged values are used for prediction. # -- If dynamic is True, then in-sample forecasts are used in place of lagged dependent variables. ax = res_111.plot_predict('1992-10', '1993-10', dynamic=True, plot_insample=False, ax=ax); 

让我们看一个更大的预测窗口（34个预测步骤）：