大家好,欢迎来到IT知识分享网。
在本文中,您将学习执行时间序列分析的基本步骤以及趋势,平稳性,移动平均等概念。还将探索指数平滑方法,并学习如何拟合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()
增强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])))
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),然后查看序列。
非季节性平滑或趋势拟合
简单指数平滑和线性指数平滑
什么时候用?
- 少数数据点
- 不规则数据
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个预测步骤):
免责声明:本站所有文章内容,图片,视频等均是来源于用户投稿和互联网及文摘转载整编而成,不代表本站观点,不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益,请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权,请联系删除,如若转载,请注明出处:https://yundeesoft.com/80563.html