如何在 python 中简单地计算时间序列的滚动/移动方差？

Posted 2023-02-14

【中文标题】如何在 python 中简单地计算时间序列的滚动/移动方差？

【英文标题】：How can I simply calculate the rolling/moving variance of a time series in python?

【发布时间】：2015-02-10 05:14:20

【问题描述】：

我有一个简单的时间序列，我正在努力估计移动窗口内的方差。更具体地说，我无法弄清楚与实现滑动窗口功能的方式有关的一些问题。例如，当使用 NumPy 且窗口大小 = 20 时：

def rolling_window(a, window):
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) 

rolling_window(data, 20)
np.var(rolling_window(data, 20), -1)
datavar=np.var(rolling_window(data, 20), -1)

也许我在某个地方弄错了，在这个思路上。 有谁知道一个简单的方法来做到这一点？ 任何帮助/建议都将受到欢迎。

【参考方案1】：

Pandas rolling_mean 和 rolling_std 函数已被弃用，取而代之的是更通用的“滚动”框架。 @elyase 的例子可以修改为：

import pandas as pd
import numpy as np
%matplotlib inline

# some sample data
ts = pd.Series(np.random.randn(1000), index=pd.date_range('1/1/2000', periods=1000)).cumsum()

#plot the time series
ts.plot(style='k--')

# calculate a 60 day rolling mean and plot
ts.rolling(window=60).mean().plot(style='k')

# add the 20 day rolling standard deviation:
ts.rolling(window=20).std().plot(style='b')

rolling 函数支持多种不同的窗口类型，如 here 所述。可以在rolling 对象上调用许多函数，包括var 和其他有趣的统计信息（skew、kurt、quantile 等）。我坚持使用std，因为该图与平均值在同一个图上，这在单位方面更有意义。

【讨论】：

“滚动方差”不是指 ts.rolling(window=20).std().plot(style='b') 中的“滚动标准”吗？

糟糕 - 我在文本中明确表示 std，但评论仍然不正确。将修复...完成

【参考方案2】：

你应该看看pandas。例如：

import pandas as pd
import numpy as np

# some sample data
ts = pd.Series(np.random.randn(1000), index=pd.date_range('1/1/2000', periods=1000)).cumsum()

#plot the time series
ts.plot(style='k--')

# calculate a 60 day rolling mean and plot
pd.rolling_mean(ts, 60).plot(style='k')

# add the 20 day rolling variance:
pd.rolling_std(ts, 20).plot(style='b')

【讨论】：

我认为 Barry 正在寻找滚动方差，而不是滚动标准差。他可以平方标准以获得方差或使用 pd.rolling_var(ts, 20).plot(style='b')。

现在，随着pandas 的更新，语法发生了变化。请参阅the docs 了解更多信息。

【参考方案3】：

尽管是一个旧线程，但我将添加另一个从 this 修改的方法，它不依赖于 pandas，也不依赖于 python 循环。本质上，使用 numpy 的步幅技巧，您可以首先创建一个具有步幅的数组视图，以便沿最后一个轴计算函数的统计量相当于执行滚动统计量。我已经修改了原始代码，通过填充添加最后一个轴的起点，输出形状与输入形状相同。

import numpy as np

def rolling_window(a, window):
    pad = np.ones(len(a.shape), dtype=np.int32)
    pad[-1] = window-1
    pad = list(zip(pad, np.zeros(len(a.shape), dtype=np.int32)))
    a = np.pad(a, pad,mode='reflect')
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

a = np.arange(30).reshape((5,6))

# rolling mean along last axis
np.mean(rolling_window(a, 3), axis=-1)

# rolling var along last axis
np.var(rolling_window(a, 3), axis=-1)

# rolling median along last axis
np.median(rolling_window(a, 3), axis=-1)

【讨论】：

感谢 np-only 解决方案。虽然，稍后我需要将头绕在填充和跨步上。现在，它确实需要我需要.. 干杯！

鉴于您最初的a.shape 为(5,6)，为什么rolling_window(a, 3) 的输出形状为(6, 6, 3)？对于(n ,m) 中的任何a.shape，都会发生这种情况，输出将始终为(n+1, m, window)。第一个维度中的那个额外点来自哪里，它应该在那里？我正在使用 Python 3.8.8 和 NumPy 1.20.1

【参考方案4】：

我只是在寻找相同的解决方案，并发现 bottleneck 包应该非常可靠且快速地解决问题。这是来自https://kwgoodman.github.io/bottleneck-doc/reference.html#bottleneck.move_var 的稍微调整的示例：

>>> import bottleneck as bn
>>> a = np.array([1.0, 2.0, 3.0, np.nan, 5.0])
>>> bn.move_var(a, window=2)
array([ nan,  0.25,  0.25,  nan,  nan])
>>> bn.move_var(a, window=2, min_count=1)
array([ 0. ,  0.25,  0.25,  0. ,  0. ])

请注意，产生的方差对应于窗口的 last 索引。

该软件包可从 Ubuntu repos、pip 等获得。它可以在 numpy-array 等的任意轴上运行。除此之外，它声称在许多情况下比普通 numpy 实现更快。

【讨论】：

瓶颈的效果非常好，但只能在 Python 3.9 中的错误有关的问题了。

【参考方案5】：

在我看来，将 Pandas 用于纯数字数据有点矫枉过正；瓶颈效果很好，但自 2021 年 1 月以来一直没有更新，并且不再适用于 Python 3.9 及更高版本；所以我会发布一个基于Josh Albert's version 的版本，请记住lib.stride_tricks.as_strided 上的文档说明，使用它可能不安全。

您可以使用 NumPy 的 lib.stride_tricks.sliding_window_view()，它基本上是 lib.stride_tricks.as_strided 周围的安全（ish）包装器，以创建一个带有额外轴的数组，该轴具有窗口大小（任意数量的维度），允许您使用 NumPy 的内置统计函数跨该轴进行操作：

import numpy as np

window = 3  # size of the window
A = np.arange(10)

Aw = np.lib.stride_tricks.sliding_window_view(A, window)
Avar = np.var(Aw, axis=-1)

Avar
>>> array([0.66666667, 0.66666667, 0.66666667, 0.66666667, 0.66666667,
       0.66666667, 0.66666667, 0.66666667])

当然，这也适用于mean、max、min、std 等。

注意：据我所知，无法包含数组的“边缘”，即无法获得完整窗口长度的A 的开头和结尾。因此，生成的数组将被缩短到可以达到完整窗口长度的那部分，请参阅返回的文档。