计算*滚动*熊猫系列的最大缩编

Question

编写一个计算时间序列最大值的函数非常容易.需要一点思考才能及时写出来O(n)而不是O(n^2)时间.但它并没有那么糟糕.这将有效:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def max_dd(ser):
    max2here = pd.expanding_max(ser)
    dd2here = ser - max2here
    return dd2here.min()

让我们设置一个简短的系列来试试吧:

np.random.seed(0)
n = 100
s = pd.Series(np.random.randn(n).cumsum())
s.plot()
plt.show()

正如预期的那样,max_dd(s)在-17.6附近显示出一些东西.好,很棒,很棒.现在说我有兴趣计算这个系列的滚动缩幅.即每个步骤,我想计算指定长度的前一个子系列的最大值.这很容易使用pd.rolling_apply.它的工作原理如下:

rolling_dd = pd.rolling_apply(s, 10, max_dd, min_periods=0)
df = pd.concat([s, rolling_dd], axis=1)
df.columns = ['s', 'rol_dd_10']
df.plot()

这非常有效.但感觉很慢.在pandas或其他工具包中是否有一个特别灵活的算法来快速完成这项工作？我开始写一些定制的东西:它跟踪各种中间数据(观察到的最大值的位置,先前发现的下降的位置),以减少许多冗余计算.它确实节省了一些时间,但不是很多,而且几乎没有尽可能多的时间.

我认为这是因为Python/Numpy/Pandas中的所有循环开销.但是我目前在Cython中还不够流利,真正知道如何从这个角度开始攻击它.我希望以前有人试过这个.或者,也许有人可能想看看我的"手工"代码,并愿意帮助我将其转换为Cython.

---

编辑:对于想要审查这里提到的所有功能(以及其他一些!)的人,请查看iPython笔记本:http://nbviewer.ipython.org/gist/8one6/8506455

它显示了这个问题的一些方法如何相关,检查它们是否给出相同的结果,并显示它们对各种大小的数据的运行时间.

如果有人有兴趣,我在帖子中提到的"定制"算法是rolling_dd_custom.我认为如果在Cython中实现,这可能是一个非常快速的解决方案.

Answer 1

这是滚动最大缩减功能的numpy版本. windowed_view是一个单行函数的包装器,用于numpy.lib.stride_tricks.as_strided创建1d数组的内存高效2d窗口视图(下面的完整代码).一旦我们有了这个窗口视图,计算基本上与你的相同max_dd,但为numpy数组编写,并沿第二轴(即axis=1)应用.

def rolling_max_dd(x, window_size, min_periods=1):
    """Compute the rolling maximum drawdown of `x`.

    `x` must be a 1d numpy array.
    `min_periods` should satisfy `1 <= min_periods <= window_size`.

    Returns an 1d array with length `len(x) - min_periods + 1`.
    """
    if min_periods < window_size:
        pad = np.empty(window_size - min_periods)
        pad.fill(x[0])
        x = np.concatenate((pad, x))
    y = windowed_view(x, window_size)
    running_max_y = np.maximum.accumulate(y, axis=1)
    dd = y - running_max_y
    return dd.min(axis=1)

这是一个演示函数的完整脚本:

import numpy as np
from numpy.lib.stride_tricks import as_strided
import pandas as pd
import matplotlib.pyplot as plt


def windowed_view(x, window_size):
    """Creat a 2d windowed view of a 1d array.

    `x` must be a 1d numpy array.

    `numpy.lib.stride_tricks.as_strided` is used to create the view.
    The data is not copied.

    Example:

    >>> x = np.array([1, 2, 3, 4, 5, 6])
    >>> windowed_view(x, 3)
    array([[1, 2, 3],
           [2, 3, 4],
           [3, 4, 5],
           [4, 5, 6]])
    """
    y = as_strided(x, shape=(x.size - window_size + 1, window_size),
                   strides=(x.strides[0], x.strides[0]))
    return y


def rolling_max_dd(x, window_size, min_periods=1):
    """Compute the rolling maximum drawdown of `x`.

    `x` must be a 1d numpy array.
    `min_periods` should satisfy `1 <= min_periods <= window_size`.

    Returns an 1d array with length `len(x) - min_periods + 1`.
    """
    if min_periods < window_size:
        pad = np.empty(window_size - min_periods)
        pad.fill(x[0])
        x = np.concatenate((pad, x))
    y = windowed_view(x, window_size)
    running_max_y = np.maximum.accumulate(y, axis=1)
    dd = y - running_max_y
    return dd.min(axis=1)


def max_dd(ser):
    max2here = pd.expanding_max(ser)
    dd2here = ser - max2here
    return dd2here.min()


if __name__ == "__main__":
    np.random.seed(0)
    n = 100
    s = pd.Series(np.random.randn(n).cumsum())

    window_length = 10

    rolling_dd = pd.rolling_apply(s, window_length, max_dd, min_periods=0)
    df = pd.concat([s, rolling_dd], axis=1)
    df.columns = ['s', 'rol_dd_%d' % window_length]
    df.plot(linewidth=3, alpha=0.4)

    my_rmdd = rolling_max_dd(s.values, window_length, min_periods=1)
    plt.plot(my_rmdd, 'g.')

    plt.show()

该图显示了代码生成的曲线.绿点的计算方法是rolling_max_dd.

时间比较,用n = 10000和window_length = 500:

In [2]: %timeit rolling_dd = pd.rolling_apply(s, window_length, max_dd, min_periods=0)
1 loops, best of 3: 247 ms per loop

In [3]: %timeit my_rmdd = rolling_max_dd(s.values, window_length, min_periods=1)
10 loops, best of 3: 38.2 ms per loop

rolling_max_dd大约快6.5倍.对于较小的窗口长度,加速更好.例如,使用window_length = 200它几乎快13倍.

要处理NA,您可以在将数组传递给之前预处理Series使用该fillna方法rolling_max_dd.