Python中的时间序列 – 指数平滑和ARIMA过程

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

大家好,欢迎来到IT知识分享网。

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

定义

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

序列的属性和类型

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

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

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

Python中的时间序列 - 指数平滑和ARIMA过程

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

Python中的时间序列 - 指数平滑和ARIMA过程

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

平稳性

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

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

用于平稳数据的变换:

  • 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() 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

增强Dickey-Fuller测试

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

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]))) 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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) 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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

选择模型

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

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

非季节性平滑或趋势拟合

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

什么时候用?

  • 少数数据点
  • 不规则数据
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(); 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

季节性平滑

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(); 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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") 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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)) 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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

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

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

Python中的时间序列 - 指数平滑和ARIMA过程

“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(); 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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

Python中的时间序列 - 指数平滑和ARIMA过程

2)选择MA order

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

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

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

Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

注意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。

Python中的时间序列 - 指数平滑和ARIMA过程

以下是我们通过拟合(1,1,1)ARIMA得到的结果:

Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

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

Python中的时间序列 - 指数平滑和ARIMA过程

在图中没有显著的自相关滞后值,我们的模型似乎不需要任何额外的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); 
Python中的时间序列 - 指数平滑和ARIMA过程

Python中的时间序列 - 指数平滑和ARIMA过程

让我们看一个更大的预测窗口(34个预测步骤):

Python中的时间序列 - 指数平滑和ARIMA过程

免责声明:本站所有文章内容,图片,视频等均是来源于用户投稿和互联网及文摘转载整编而成,不代表本站观点,不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益,请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权,请联系删除,如若转载,请注明出处:https://yundeesoft.com/80563.html

(0)
上一篇 2024-08-29 10:26
下一篇 2024-08-29 19:45

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注

关注微信