用于日期时间的最快滚动窗口方法

编程语言 2026-07-09

我现在在循环中对时间数据的一个30分钟滚动窗口应用一个有点复杂的函数,结果发现花费的时间太长,已经不太合理。

下面给出一个最小可复现的示例,目前耗时约94 ms,结合数据集的实际规模和循环情况,这个时间显然不可行。

import pandas as pd
import numpy as np
from datetime import datetime,timedelta
import statsmodels.tsa.api as smt

def instrument_variance(ts): #Lernschow et al 2000
    ts_acov = smt.acovf(ts)[:3] #Find first few lags of autocovariance
    M_0 = ts_acov[0] #Get autocovariance at lag 0
    M_fit = np.polyfit(range(3),ts_acov,1) #Get linear fit of the first few lags
    M_lim0 = M_fit[0] #Get limit as fit goes to 0
    var_inst = M_0-M_lim0 #Get error variance from M11(0)-lim>0(M11)
    return var_inst

start = datetime.now()
base = datetime.today()
date_list = [base + timedelta(minutes=x) for x in range(100)]
df = pd.DataFrame(np.random.randint(0,100,size=(100, 3)), index=date_list,columns=list('ABC'))
df_cop = df.copy()

var_inst = df_cop.sort_index().rolling('30min',center=True,min_periods=3).apply(instrument_variance)
df = df.sort_index().rolling('30min',center=True,min_periods=3).var()
df = df.sub(var_inst) #Subtract uncorrelated var from total var to get correlated var
df = df.resample('5min').first()

timed = datetime.now()-start
print(timed)

我的Python知识有限,操作也比较复杂,但有没有什么方法可以让它更快一些?

我也看过类似的问题:

在Python(Numpy或 Pandas)中计算滚动(线性)回归的最有效/最快方法

但我很难找到一个针对基于日期时间的窗口的解决方案。我不能假设采样频率是恒定的,因此我更愿意使用datetime而不是整数索引的方式来实现。

解决方案

not sure if it will be lean enough for your dataset... run_fast_prefix gets about 45x speed-up:

tier 1 optimizations of per-window overhead (37ms → 2.9ms, ~13x). in run_fast:

  • 将inline acovf + polyfit 结合使用。每次调用在一个6 元素窗口上都有校验/奇异值分解的开销。用少量的numpy运算和闭式斜率 (M2 - M0)/2取代。
  • 将两个滚动遍历合并成一个。原始实现对滚动 .apply 和滚动 .var 的处理都是对窗口的完整遍历。
  • raw=True 跳过为每个窗口构建一个pandas Series。

tier 2优化,灵感来自你分享的另一个SO问题的答案。完全消除每个窗口的调用,在 run_fast_prefix

  • 前缀和:对整列进行一次 cumsum,然后每个窗口的 sum 只需两个数组查找。所有 ~100个窗口的M0、M1、M2变成了四个向量表达式。逐窗口的Python解释器循环消失,所有工作都在已编译的numpy中完成。

总速度提升:约45x。tier 1是在现有结构上的纯常数因子优化;tier 2是链接的SO问题中所指向的结构性变动。

import numpy as np
import pandas as pd
import statsmodels.tsa.api as smt


# original for reference
def instrument_variance(ts):  # Lenschow et al 2000
    ts_acov = smt.acovf(ts)[:3]
    M_0 = ts_acov[0]
    M_fit = np.polyfit(range(3), ts_acov, 1)
    M_lim0 = M_fit[0]
    var_inst = M_0 - M_lim0
    return var_inst


def run_original(df):
    var_inst = (df.sort_index()
                  .rolling('30min', center=True, min_periods=3)
                  .apply(instrument_variance))
    var_tot = (df.sort_index()
                 .rolling('30min', center=True, min_periods=3)
                 .var())
    return var_tot.sub(var_inst).resample('5min').first()


# combined `apply` path
def correlated_var(a):
    """Rolling correlated variance in one pass.

    Biased autocovariance (÷ n) for the instrument term, ddof=1 total var.
    Slope of the line through (0,M0),(1,M1),(2,M2) is (M2 - M0)/2 in closed
    form, reproducing `np.polyfit(...)[0]` semantics from the original. If
    the intent is the lag-0 intercept, use: inst = M0 - ((M0+M1+M2)/3 - slope).
    """
    n = a.size
    if n < 3:
        return np.nan
    c = a - a.mean()
    ssq = (c * c).sum()
    M0 = ssq / n
    M1 = (c[:-1] * c[1:]).sum() / n
    M2 = (c[:-2] * c[2:]).sum() / n
    slope = 0.5 * (M2 - M0)
    inst = M0 - slope
    var_u = ssq / (n - 1)
    return var_u - inst


def run_fast(df):
    return (df.sort_index()
              .rolling('30min', center=True, min_periods=3)
              .apply(correlated_var, raw=True)
              .resample('5min').first())


# prefix-sum vectorized path
def _rolling_bounds(df, offset, center=True, closed=None):
    """Exact [start, end) bounds that df.rolling(offset, center=..., closed=...) uses.

    Extracted by running pandas' own rolling on a row-index series with a
    scalar-returning apply. This is the robust way — pandas' offset-based
    centering is implemented as an output-alignment shift, not a bounds shift,
    so VariableOffsetWindowIndexer.get_window_bounds(center=True) and any hand-
    derived (t - w/2, t + w/2] formula both drift by a row in places.

    The bootstrap apply itself is O(N) but cheap per window (scalar return),
    so downstream prefix-sum math still gets most of the speedup.
    """
    N = len(df)
    if N == 0:
        return np.zeros(0, dtype=np.int64), np.zeros(0, dtype=np.int64)
    idx_ser = pd.Series(np.arange(N, dtype=np.float64), index=df.index)
    r = idx_ser.rolling(offset, center=center, closed=closed, min_periods=1)
    first = r.apply(lambda w: w[0],  raw=True).to_numpy()
    last  = r.apply(lambda w: w[-1], raw=True).to_numpy()
    empty = np.isnan(first)
    start = np.where(empty, 0, first).astype(np.int64)
    end   = np.where(empty, 0, last + 1).astype(np.int64)
    return start, end


def _correlated_var_prefix(x, start, end, min_periods=3):
    """Vectorized correlated-variance per window via prefix sums.

    Window for output i is x[start[i]:end[i]]. Requires min_periods >= 3
    (lag-2 autocov needs at least 3 samples).
    """
    assert min_periods >= 3
    x = np.asarray(x, dtype=float)
    N = x.size
    out = np.full(start.shape, np.nan, dtype=float)
    if N == 0:
        return out

    # Padded prefix sums: S[0]=0, S[k]=sum(x[:k]) so sum(x[a:b]) = S[b]-S[a].
    S = np.empty(N + 1, dtype=float);  S[0] = 0.0
    np.cumsum(x, out=S[1:])
    Sq = np.empty(N + 1, dtype=float); Sq[0] = 0.0
    np.cumsum(x * x, out=Sq[1:])

    SP1 = np.zeros(max(N, 1), dtype=float)         # lag-1 pair-product cumsum
    if N >= 2:
        np.cumsum(x[:-1] * x[1:], out=SP1[1:])
    SP2 = np.zeros(max(N - 1, 1), dtype=float)     # lag-2 pair-product cumsum
    if N >= 3:
        np.cumsum(x[:-2] * x[2:], out=SP2[1:])

    n = (end - start).astype(np.int64)
    valid = n >= min_periods
    if not valid.any():
        return out

    s = start[valid]
    e = end[valid]
    nn = n[valid].astype(float)

    Sx = S[e] - S[s]
    Sxx = Sq[e] - Sq[s]
    mean = Sx / nn
    ssq = Sxx - nn * mean * mean
    M0 = ssq / nn
    var_u = ssq / (nn - 1.0)

    # Lag 1: sum over pairs (x[t], x[t+1]) for t in [s, e-1).
    sum_P1 = SP1[e - 1] - SP1[s]
    sumL1 = Sx - x[e - 1]                          # sum of left members
    sumR1 = Sx - x[s]                              # sum of right members
    n1 = nn - 1.0
    M1 = (sum_P1 - mean * (sumL1 + sumR1) + n1 * mean * mean) / nn

    # Lag 2: sum over pairs (x[t], x[t+2]) for t in [s, e-2).
    sum_P2 = SP2[e - 2] - SP2[s]
    sumL2 = Sx - x[e - 2] - x[e - 1]
    sumR2 = Sx - x[s] - x[s + 1]
    n2 = nn - 2.0
    M2 = (sum_P2 - mean * (sumL2 + sumR2) + n2 * mean * mean) / nn

    slope = 0.5 * (M2 - M0)
    inst = M0 - slope
    out[valid] = var_u - inst
    return out


def run_fast_prefix(df, offset='30min', min_periods=3, resample_rule='5min'):
    df = df.sort_index()
    start, end = _rolling_bounds(df, offset, center=True)

    arr = df.to_numpy(dtype=float, copy=False)
    out = np.empty_like(arr, dtype=float)
    for c in range(arr.shape[1]):
        out[:, c] = _correlated_var_prefix(arr[:, c], start, end, min_periods)

    result = pd.DataFrame(out, index=df.index, columns=df.columns)
    return result.resample(resample_rule).first() if resample_rule else result

你可能已经知道:numpy与 pandas在处理NaN时有差异。如果数据中有NaN,请在运行前进行插值或删除。

站内所有文章版权归属LeftHeroAI导航站,无授权禁止任何主体转载、抄袭、复制内容,亦不得私自架设镜像站点。一经侵权,本站将通过法律途径追责。

相关文章