乐趣区

关于算法:Python面板时间序列数据预测格兰杰因果关系检验Granger-causality-test药品销售实例与可视化

jiezi

原文链接：http://tecdat.cn/?p=23940

工夫序列是以固定工夫_区间_记录的察看序列。本指南带你实现在 Python 中剖析一个给定的工夫序列的特色的过程。

内容

  1. 什么是工夫序列？
  2. 如何在 Python 中导入工夫序列？
  3. 什么是面板数据？
  4. 工夫序列的可视化
  5. 工夫序列中的模式
  6. 加法和乘法的工夫序列
  7. 如何将一个工夫序列分解成其组成部分？
  8. 安稳的和非安稳的工夫序列
  9. 如何使一个工夫序列成为安稳的？
  10. 如何测试平稳性？
  11. 白噪声和安稳序列之间的区别是什么？
  12. 如何使一个工夫序列去趋势化？
  13. 如何使工夫序列去节令化？
  14. 如何测验工夫序列的季节性？
  15. 如何解决工夫序列中的缺失值？
  16. 什么是自相干和局部自相干函数？
  17. 如何计算局部自相干函数？
  18. 滞后图
  19. 如何预计一个工夫序列的可预测性？
  20. 为什么和如何使工夫序列平滑化？
  21. 如何应用格兰杰因果测验来理解一个工夫序列是否有助于预测另一个工夫序列？

1. 什么是工夫序列？

工夫序列是以固定工夫区间记录的察看序列。

依据察看的频率，一个工夫序列通常可能是每小时、每天、每周、每月、每季度和每年。有时，你也可能有以秒为单位的工夫序列，比方，每分钟的点击量和用户访问量等等。

为什么要剖析一个工夫序列？

因为这是你对该序列进行预测前的筹备步骤。

此外，工夫序列预测具备微小的商业意义，因为对企业来说很重要的货色，如需要和销售，网站的访问量，股票价格等基本上都是工夫序列数据。

那么，剖析一个工夫序列波及什么呢？

工夫序列剖析波及到对序列性质的各个方面的了解，这样你就能更好地理解发明有意义和精确的预测。

2. 如何在 Python 中导入工夫序列？

那么，如何导入工夫序列数据呢？

工夫序列的数据通常存储在.csv 文件或其余电子表格格局中，蕴含两列：日期和测量值。

咱们应用 pandas 包中的 read\_csv() 来读取工夫序列数据集（一个对于药品销售的 csv 文件）作为 pandas 数据框。增加 parse\_dates=[‘date’] 参数将使日期列被解析为一个日期字段。

import pandas as pd

# 导入数据
df = pd.read_csv('10.csv', parse_dates=\['date'\])
df.head()

数据框架工夫序列

另外，你也能够把它导入为一个以日期为索引的 pandas 序列。你只须要在 pd.read\_csv() 中指定 index\_col 参数就能够了。

pd.read\_csv('10.csv', parse\_dates=\['date'\], index_col='date')

3. 什么是面板数据？

面板数据也是一种基于工夫的数据集。

不同的是，除了工夫序列之外，它还蕴含一个或多个在雷同时间段内测量的相干变量。

通常状况下，面板数据中存在的列蕴含了有助于预测 Y 的解释变量，前提是这些列在将来的预测期是可用的。

上面是一个面板数据的例子。

df.head()

面板序列

4. 工夫序列的可视化

让咱们用 matplotlib 来可视化这个序列。

# 工夫序列数据源：R 中的 fpp pacakge。import matplotlib.pyplot as plt

# 绘制图表
def plot_df(df, x, y, title="", xlabel=' 日期 ', dpi=100):

plt.show()

工夫序列的可视化

因为所有的值都是负数，你能够在 Y 轴的两边显示，以强调增长。

# 导入数据
x = df\['date'\].values

# 绘图
fig, ax = plt.subplots(1, 1, figsize=（16,5）, dpi= 120)
plt.fill(x, y1=y1, y2=-y1, alpha=0.5)

航空乘客数据 – 两面序列

因为它是一个月度的工夫序列，并且每年都遵循肯定的反复模式，你能够在同一张图中把每年的状况作为一个独自的线条来绘制。这让你能够并排比拟每年的模式。

工夫序列的节令图 

# 导入数据

df.reset_index(inplace=True)

# 筹备好数据

years = df\['year'\].unique()

# 准备色彩

np.random.choice(list(mpl.color), len(year), 

# 绘制图表
  
plt.text(df.loc\[df.year==y, :\].shape\[0\]-.9\] 


plt.gca().set(xlim=(-0.3, 11)

plt.title("药品销售工夫序列的节令图", fontsize=20)

药品销售的季节性图谱

每年 2 月，药品销售量急剧下降，3 月再次回升，4 月再次降落，如此重复。显然，这种模式每年都会在某一年内反复呈现。

然而，随着工夫的推移，药品销售量总体上有所增加。你能够用一个丑陋的年度图表很好地展现这一趋势以及它每年的变动状况。同样地，你也能够做一个按月排列的 boxplot 来显示每月的散布状况。

逐月（季节性）和逐年（趋势）散布的箱线图

你能够按季节性对数据进行分组，看看数值在某年或某月是如何散布的，以及它在不同期间的比照状况。

# 导入数据

df.reset_index(inplace=True)

# 筹备好数据
df\['年'\] = \[d.year for d in df.date\]
df\['月'\] = \[d.strftime('%b') for d in df.date\]


# 绘制图表
sns.boxplot(x='年', y='值', data=df, ax=axes\[0\])
sns.boxplot(x='月', y='值', data=df.loc\[~df.year.isin(\[1991, 2008\]), :\] )

按年和按月排列的箱线图

箱线图使年度和月份的散布变得显著。另外，在按月排列的图表中，12 月和 1 月的药品销售量显著较高，这可归因于假日折扣节令。

到目前为止，咱们曾经看到了识别模式的相似性。当初，如何找出与通常模式的任何偏差？

5. 工夫序列中的模式

任何工夫序列都能够被分成以下几个局部。根底程度 + 趋势 + 季节性 + 误差

当在工夫序列中察看到有一个减少或缩小的斜率时，就能够察看到趋势。而季节性是指因为季节性因素，在定期区间之间察看到显著的反复模式。这可能是因为一年中的哪个月，哪个月的哪一天，工作日或甚至一天中的哪个工夫。

然而，并非所有工夫序列都必须有趋势和 / 或季节性。一个工夫序列可能没有一个显著的趋势，但有一个季节性。反之，也能够是实在的。

因而，一个工夫序列能够被设想为趋势、季节性和误差项的组合。

fig, axes = plt.subplots(1,3, figsize=(20,4), dpi=100)


pd.read_csv.plot(legend=False, ax=axes\[2\])

工夫序列中的模式

另一个须要思考的方面是周期性行为。当序列中的回升和降落模式不产生在固定的基于日历的工夫区间内时，就会产生这种状况。应留神不要将 “ 周期性 “ 效应与 “ 季节性 “ 效应混同。

那么，如何辨别 “ 周期性 “ 和 “ 季节性 “ 模式？

如果这些模式不是基于固定的日历频率，那么它就是周期性的。因为，与季节性不同，周期性效应通常受到商业和其余社会经济因素的影响。

6. 加法和乘法的工夫序列

依据趋势和季节性的性质，一个工夫序列能够被建模为加法或乘法，其中，序列中的每个观测值能够示意为各组成部分的和或积。

加法工夫序列: 值 = 根底 + 趋势 + 季节性 + 误差
 

乘法工夫序列: 值 = 根底 x 趋势 x 季节性 x 误差

7. 如何将一个工夫序列分解成其组成部分？

你能够对一个工夫序列进行经典的合成，将该序列视为基数、趋势、季节性指数和残差的加法或乘法组合。

statsmodels 中的 seasonal_decompose 能够不便地实现这一点。

# 乘法合成 
decompose(df\['value'\], model='multiplicative')

# 加法合成
decompose(df\['value'\], model='additive')

# 绘图
result_mul.plot().suptitle( fontsize=22)

加法和乘法合成

设置 extrapolate_trend=’freq’ 能够看到序列开始时趋势和残差中的任何缺失值。

如果你认真看一下加法合成的残差，它有一些模式残留。然而，乘法合成看起来相当随机。因而，现实状况下，对于这个特定的序列，乘法合成应该是首选。

趋势、季节性和残差成分的数字输入存储在 result_mul 输入自身。让咱们把它们提取进去，放在一个数据框中。

# 提取成分 ----
# 理论值 =（季节性 * 趋势 * 残差）的乘积

constructed.columns = \['seas', 'trendance', 'resid', 'actual_values'\]。constructed.head()

如果你检查一下，seas、trend 和 resid 列的乘积应该正好等于 actual_values。

8. 安稳和非安稳的工夫序列

平稳性是工夫序列的一个属性。一个安稳的序列是指该序列的值不是工夫的函数。

也就是说，序列的统计属性，如平均数、方差和自相干，随着工夫的推移是不变的。序列的自相干只不过是序列与它以前的值的相关性，更多的是对于这一点。

一个安稳的工夫序列也没有季节性影响。

那么，如何辨认一个序列是否是安稳的？让咱们绘制一些例子来阐明。

安稳和非安稳的工夫序列

那么，为什么安稳的序列很重要呢？

我稍后会说到这个问题，但要晓得，通过利用适当的转换，简直能够使任何工夫序列成为安稳的。大多数统计预测办法都被设计为在安稳的工夫序列上工作。预测过程的第一步通常是做一些转换，把非安稳序列转换成安稳的。

9. 如何使一个工夫序列安稳？

你能够通过以下形式使序列安稳

  1. 对序列进行差分（一次或屡次
  2. 取序列的对数
  3. 取序列的第 n 次根
  4. 上述办法的组合

使序列安稳的最常见和最不便的办法是对数列至多进行一次差分，直到它成为近似安稳的。

那么，什么是差分？

如果 Y_t 是工夫 ’t’ 的值，那么 Y 的第一个差值 =Yt-Yt-1。更简略地说，序列只不过是用以后值减去下一个值。

如果第一次差分不能使一个序列安稳，你能够进行第二次差分。以此类推。

例如，思考以下序列。[1, 5, 2, 12, 20]

第一次差分能够失去。[5-1, 2-5, 12-2, 20-12] = [4, -3, 10, 8]

二次差分得出。[-3-4, -10-3, 8-10] = [-7, -13, -2]

9. 为什么在预测前要使非稳态序列成为稳态？

预测一个安稳的序列是绝对容易的，而且预测后果也更牢靠。

一个重要的起因是，自回归预测模型实质上是线性回归模型，利用序列自身的滞后期作为预测因子。

咱们晓得，如果预测因子（X 变量）不互相关联，线性回归成果最好。因而，序列的安稳化解决了这个问题，因为它打消了任何继续的自相干，从而使预测模型中的预测因子（序列的滞后）简直是独立的。

当初咱们曾经确定了序列安稳化的重要性，那么你如何查看一个给定的序列是否是安稳的？

10. 如何测验平稳性？

一个序列的平稳性能够通过观察序列的图表来确定，就像咱们之前做的那样。

另一种办法是将序列分成 2 个或更多的间断局部，并计算平均数、方差和自相干等汇总统计数据。如果统计数字有很大差别，那么这个序列就不可能是安稳的。

然而，你须要一种办法来定量地确定一个给定的序列是否是安稳的。这能够通过称为 “ 单位根测试 “ 的统计测试来实现。这方面有多种变动，测试查看一个工夫序列是否是非稳态的并领有单位根。

单位根测试有多种实现形式，例如。

  1. Augmented Dickey Fuller 测验 (ADH Test)
  2. Kwiatkowski-Phillips-Schmidt-Shin – KPSS 测验 (趋势安稳)
  3. Philips Perron 测验 (PP 测验)

最罕用的是 ADF 测验，无效假设是工夫序列领有单位根，并且是非安稳的。所以，如果 ADH 测验中的 P 值小于显著性程度（0.05），你就回绝无效假设。

另一方面，KPSS 测验是用来测验趋势平稳性的。空假如和 P 值的解释与 ADH 测验正好相同。上面的代码应用 python 中的 statsmodels 包来实现这两个测验。

# ADF 测试
result = adfuller(df.value.values, autolag='AIC')

# KPSS 测试
result = kpss(df.value.values, regression='c')
print('\\nKPSS Statistic: %f' % result\[0\])

 

11. 白噪声和安稳数列之间有什么区别？

与安稳序列一样，白噪声也不是工夫的函数，即它的平均值和方差不随工夫变动。但不同的是，白噪声是齐全随机的，平均值为 0。

在白噪声中，没有任何模式。如果你把调频收音机中的声音信号看作是一个工夫序列，你在各频道之间听到的空白声音就是白噪声。

在数学上，一个平均数为 0 的齐全随机的数字序列就是白噪声。

 np.random.randn(1000)

随机白噪声

12. 如何对工夫序列进行去趋势解决？

去工夫序列的趋势是指从一个工夫序列中去除趋势成分。但如何提取趋势呢？有多种办法。

  1. 从工夫序列中去除最佳拟合线。最佳拟合线能够从以工夫步骤为预测因素的线性回归模型中取得。对于更简单的趋势，你可能想在模型中应用二次项（x^2）。
  2. 去除咱们后面看到的从工夫序列合成中失去的趋势成分。
  3. 去除平均数
  4. 利用像 Baxter-King 滤波器或 Hodrick-Prescott 滤波器这样的滤波器来去除挪动均匀趋势线或周期成分。

让咱们来实现前两种办法。

# 应用 scipy。减去最佳拟合线
 signal.detrend(df.value.values)
plt.plot(detrended)

通过减去最小二乘法来解读工夫序列的趋势 

# 应用 statmodels。减去趋势成分。decompose(df\['value'\], model='multiplicative', 
detrended = df.value.value - trend 
plt.plot(detrended)

通过减去趋势成分进行去势

13. 如何对一个工夫序列进行去节令化？

有多种办法能够使工夫序列去季候性化。上面是几个例子。

- 1. 取一个以季节性窗口为长度的挪动平均线。这将在这个过程中使序列变得平滑。- 2. 序列的季节性差别（用以后值减去前一季的值）。- 3. 用从 STL 合成失去的季节性指数除以该序列。

如果除以季节性指数成果不好，能够尝试取序列的对数，而后进行去季节性解决。之后你能够通过取指数来复原到原来的规模。

# 工夫序列合成 
seasonal_decompose(df\['value'\], model='multiplicative'')
 # 去节令化 
df. value.values / result_mul.seasonal 
# 绘图
plt.plot(deseason)

对工夫序列进行反季节解决

14. 如何测试一个工夫序列的季节性？

常见的办法是绘制序列图，查看安稳工夫区间内的可反复模式。所以，季节性的类型是由时钟或日历决定的。

  1. 一天中的小时
  2. 月的一天
  3. 每周
  4. 月度
  5. 年

然而，如果你想对季节性有一个更明确的查看，能够应用自相干函数（ACF）图。更多对于 ACF 的内容将在接下来的章节中介绍。然而，当有强烈的季节性模式时，ACF 图通常会显示出在季节性窗口的倍数上有明确的反复峰值。

例如，药品销售工夫序列是一个月度序列，每年都有反复的模式。因而，你能够看到在第 12、24、36…… 行的尖峰。

在实在的数据集中，这种强烈的模式很难被留神到，并可能被任何乐音所扭曲，所以你须要仔细观察。

 # 画图
correlation_plot(df.value.tolist())

自相干图

另外，如果你想进行统计测试，CHTest 能够确定是否须要进行季节性差分以使序列安稳化。

15. 如何解决工夫序列中的缺失值？

有时，你的工夫序列会有缺失的日期 / 工夫。这意味着，这些期间的数据没有被捕捉或无奈取得。在这种状况下，你能够用零来填补这些期间。

其次，当波及到工夫序列时，你通常不应该用序列的平均值来替换缺失值，特地是在序列不是安稳的状况下。你能够做的是是向前填充前一个值。

然而，依据序列的性质，你要在得出结论之前尝试多种办法。一些无效的代替归因法的办法是。

  • 后向填充
  • 线性插值
  • 二次插值
  • 最近邻均匀
  • 季节性差值均匀

为了掂量归因性能，我手动引入工夫序列的缺失值，用上述办法进行归因，而后掂量归因与理论值的均匀平方误差。

# 生成数据集

prcPupdate({'xtick.bottom' : False})
 
## 1. Actual -------------------------------
df_or.plot(style=".-")

 
## 2. 正向填充 --------------------------
df_ffill = df.ffill()

 
## 3. 后向填充 -------------------------
df_bfill = df.bfill()

## 4. 线性插值 ------------------
df\['rownum'\] = np.arange(df.shape\[0\])

## 5. 平面插值 --------------------
f2 = interp1d(df\_um'\], df\_nona\['v kind='cubic')

# 内插法参考。## 6. n " 过来最近的街坊的平均值 ------

 knnmean(df.value.values, 8)
 np.round(meansquarerr('), 2)

## 7. 季节性平均值 ----------------------------
    """
    计算相应季节性期间的平均值
    ts: 工夫序列的一维数组式
    n: 工夫序列的季节性窗口长度
    """
    out = np.cpy(ts)
    for i, val in merate(ts):
        if np.isnan(val):
            ts_seas = ts\[i-1::-n\] # 只有以前的节令
            如果 np.isan(np.naean(ts_seas))。ts_seas = np.concenate(\[ts\[i-1::-n\], ts\[i::n\]\]) # 以前和当前的
            out\[i\] = np.nanan(ts_seas) * lr
 
seasonal_mean(df.value, n=12, lr=1.25)

缺失值解决办法

你也能够思考以下办法，这取决于你心愿推断的精确水平。

  1. 如果你有解释变量，应用预测模型，如随机森林或 k -Nearest Neighbors 来预测。
  2. 如果你有足够的过来观测值，就预测缺失的值。
  3. 如果你有足够的将来观测值，就对缺失值进行反向预测
  4. 预测以前周期的对应值。

16. 什么是自相干和偏自相干函数？

自相干只是一个序列与它本人的滞后期的相干关系。如果一个序列是显著的自相干，那就意味着，该序列以前的值（滞后）可能有助于预测以后的值。

偏自相干也传播了相似的信息，但它传播的是一个序列与其滞后期的纯相干，排除了两头滞后期的相干奉献。

 # 计算 ACF 和 PACF 直到 50 个滞后期
acf_50 = acf(value, nlags=50)
pacf_50 = pacf(value, nlags=50)
 
# 绘制图表
fig, axes = plt.subplots(1,2,figsize=（16,3）)

ACF 和 PACF

17. 如何计算偏自相干函数？

那么，如何计算偏自相干？

一个序列的滞后期（k）的偏自相干是该滞后期在 Y 的自回归方程中的系数。Y 的自回归方程只不过是以其本身的滞后期为预测因素的 Y 的线性回归。

例如，如果 Y\_t 是以后序列，Y\_t- 1 是 Y 的滞后 1，那么滞后 3 的偏自相干（Y\_t-3）是 Y\_t- 3 的系数 $alpha_3$，在以下方程中。

自回归方程

18. 滞后图

滞后图是一个工夫序列与本身滞后的散点图。它通常是用来查看自相干的。如果在序列中存在任何像你上面看到的模式，该序列是自相干的。如果没有这样的模式，该序列可能是随机白噪声。

在上面对于太阳黑子区域工夫序列的例子中，随着 n_lag 的减少，图变得越来越扩散。

# #导入

a10 = pd.read_csv('10.csv')
 
# 绘图
for i, ax in enumerate(axes.flatten()\[:4\]):

 
fig.suptitle('药物销售的滞后图', y=1.05)

药物销售滞后曲线图

滞后图 太阳黑子

19. 如何预计一个工夫序列的可预测性？

一个工夫序列的规律性和可重复性越强，就越容易进行预测。近似熵 “ 能够用来量化一个工夫序列中稳定的规律性和不可预测性。

近似熵越高，预测就越艰难。

另一个更好的代替办法是 “ 样本熵 ”。

样本熵与近似熵类似，但在预计复杂性方面更加统一，即便是较小的工夫序列。例如，一个数据点较少的随机工夫序列的 “ 近似熵 “ 可能比一个较 “ 规定 “ 的工夫序列低，而一个较长的随机工夫序列的 “ 近似熵 “ 会更高。

样本熵很好地解决了这个问题。请看上面的演示。

 def AEn(U, m, r):
    ""计算 Aproximate entropy""
    def \_maxdit(x\_i, x_j):
        returnmax(\[abs(u- va) for ua, va in zp(x\_i, x\_j)\])
 
    def _phi(m):
        x = \[\[U\[j\]for in range(i, i \] for i in range(N - m + 1)\]
        C = \[len(\[1 for x\_jn x if \_maist(x\_i, x\_j) <= r\]) / (N - m + 1.0) for xi in x\]。返回 (N - m +.0)**(-1) * sump.log(C))
 
    N = len(U)

def Samp:
    ""计算样本熵"""
    maxstx\_i, x\_j):
        retur max(\[abs(ua - va) fo ua, vain zip(\_i x\_j)\])
   phi(m)
        x = \[\[j\]or j i ane(i, i  m - 1 + 1\] for iin age(N - m + 1)\]
        C= \[len(\[1 for j in rane(len(x)) if i != j and _mist(x\[i, xj\]) <= r\]) for i in rane(ln(x)) \]
     = en()

20. 为什么和如何对工夫序列进行平滑解决？

对一个工夫序列进行平滑解决可能在以下方面有用。

  • 缩小信号中噪声的影响，失去一个通过噪声过滤的序列的近似值。
  • 平滑化后的序列能够作为解释原始序列自身的一个特色。
  • 更好地察看根本趋势

那么，如何对一个序列进行平滑解决？让咱们讨论一下以下办法。

  1. 取一个挪动平均线
  2. 做一个 LOESS 平滑（部分回归）。
  3. 做一个 LOWESS 平滑（部分加权回归）。

挪动平均数只不过是定义宽度的滚动窗口的平均值。但你必须正当地抉择窗口宽度，因为大的窗口尺寸会使序列适度平滑。例如，窗口尺寸等于季节性持续时间（例如：12 个月的序列），将无效地打消季节性效应。

LOESS 是 “LOcalized regrESSion “ 的简称，它在每个点的部分左近进行屡次回归。它是在 statsmodels 软件包中实现的，你能够用 frac 参数来管制平滑的水平，该参数指定了左近的数据点的百分比，应该被视为适宜回归模型。

# 导入
pd.read\_csv('ele.csv', parse\_dates=\['date'\], index_col='date')

# 1. 挪动平均数
dma = df_rg.vale.r.man()

# 2. 平滑 (5% 和 15%)
pd.DtaFame(lowess(dfoig.alu, np.ane(len(d_origale)), fac=0.05)

# 绘图
fig, axes = plt.ubplos(4,1, figsiz=（7, 7, sharex=rue, dp=120)
df_ori\['aue'\].pot(ax=axes0\], color='k')

平滑化工夫序列

如何应用格兰杰因果测试来理解一个工夫序列是否有助于预测另一个工夫序列？

格兰杰因果测验是用来确定一个工夫序列是否有助于预测另一个工夫序列的。

格兰杰因果关系测试是如何工作的？

它是基于这样的想法：如果 X 导致 Y，那么基于 Y 的前值和 X 的前值对 Y 的预测应该优于仅基于 Y 的前值的预测。

因而，了解格兰杰因果关系不应该被用来测试 Y 的滞后期是否导致 Y。

无效假设是：第二列中的序列不会导致第一列中的序列的格兰杰。如果 P 值小于显著性程度（0.05），那么你就回绝无效假设，并得出结论：上述 X 的滞后期的确是有用的。

第二个参数 maxlag 说的是在测试中应该包含多少个 Y 的滞后期。

df = pd.rea\_csv('a10.csv', parse\_dates=\['date'\])
df\['moth'\] = df.date.dt.nth
gragecaslitess(df\[\['alue', 'moh'\]\], maxag2)

在上述情况下，所有测验的 P 值都是零。因而，” 月 “ 数据的确能够用来预测航空乘客。

最受欢迎的见解

1. 在 python 中应用 lstm 和 pytorch 进行工夫序列预测

2.python 中利用长短期记忆模型 lstm 进行工夫序列预测剖析

3. 应用 r 语言进行工夫序列（arima，指数平滑）剖析

4.r 语言多元 copula-garch- 模型工夫序列预测

5.r 语言 copulas 和金融工夫序列案例

6. 应用 r 语言随机稳定模型 sv 解决工夫序列中的随机稳定

7.r 语言工夫序列 tar 阈值自回归模型

8.r 语言 k -shape 工夫序列聚类办法对股票价格工夫序列聚类

9.python3 用 arima 模型进行工夫序列预测

退出移动版