Pandas 中的 EWMA 协方差矩阵 - 优化

Question

我想使用 Pandas 从股票价格回报的 DataFrame 计算 EWMA 协方差矩阵，并遵循 PyPortfolioOpt 中的方法。

我喜欢使用 Pandas 对象和函数的灵活性，但是当资产集增长时，函数变得非常缓慢：

import pandas as pd
import numpy as np

def ewma_cov_pairwise_pd(x, y, alpha=0.06):
    x = x.mask(y.isnull(), np.nan)
    y = y.mask(x.isnull(), np.nan)
    covariation = ((x - x.mean()) * (y - y.mean()).dropna()
    return covariation.ewm(alpha=0.06).mean().iloc[-1]

def ewma_cov_pd(rets, alpha=0.06):
    assets = rets.columns
    n = len(assets)
    cov = np.zeros((n, n))
    for i in range(n):
        for j in range(i, n):
            cov[i, j] = cov[j, i] = ewma_cov_pairwise_pd(
                rets.iloc[:, i], rets.iloc[:, j], alpha=alpha)
    return pd.DataFrame(cov, columns=assets, index=assets)

我希望在仍然使用 Pandas 的同时提高代码速度，但瓶颈在于 DataFrame.ewm() 函数，它使用了 90% 的计算时间。

如果使用此函数是一个绑定约束，那么提高代码运行速度的最有效方法是什么？我正在考虑采用暴力方法并使用 concurrent.futures.ProcessPoolExecutor但也许有更好的解决方案。

n = 100  # n is typically 2000
rets = pd.DataFrame(np.random.normal(0, 1., size=(n, n)))
cov_pd = ewma_cov_pd(rets)

真正的时间序列数据可以包含前导空值和可能的缺失值，尽管后者不太可能。

更新我

利用 Quang Hoang 提供的答案并在更合理的时间内产生预期结果的潜在解决方案类似于：

def ewma_cov_frame_qh(rets, alpha=0.06):
    weights = (1-alpha) ** np.arange(len(df))[::-1]
    normalized = (rets-rets.mean()).to_numpy()    
    out = (weights * normalized.T) @ normalized / weights.sum()
    return pd.DataFrame(out, index=rets.columns, columns=rets.columns)


def ewma_cov_qh(rets, alpha=0.06):
    syms = rets.columns
    covar = pd.DataFrame(index=rets.columns, columns=rets.columns)
    delta = rets.isnull().sum(axis=1).shift(1) - rets.isnull().sum(axis=1)
    dates = delta.loc[delta != 0].index.tolist()
     
    for date in dates:
        frame = rets.loc[rets.index >= date].dropna(axis=1, how='any')
        cov = ewma_cov_frame_qh(frame).reindex(index=syms, columns=syms)
        covar = covar.fillna(cov)
   
    return covar

cov_qh = ewma_cov_qh(rets)

这违反了使用原生 Pandas/Numpy 函数计算底层协方差的要求，并且计算时间将取决于数据集中前导 na 的数字。

更新二

下面列出了对上述内容的潜在改进，它使用（天真的实现）多处理并在我的机器上将计算时间进一步提高了 42.5%：

from concurrent.futures import ProcessPoolExecutor, as_completed
from functools import partial
    
def ewma_cov_mp_worker(date, rets, alpha=0.06):
    syms = rets.columns
    frame = rets.loc[rets.index >= date].dropna(axis=1, how='any')
    return ewma_cov_frame_qh(frame, alpha=alpha).reindex(index=syms, columns=syms)


def ewma_cov_mp(rets, alpha=0.06):
    covar = pd.DataFrame(index=rets.columns, columns=rets.columns)
    delta = rets.isnull().sum(axis=1).shift(1) - rets.isnull().sum(axis=1)
    dates = delta.loc[delta != 0].index.tolist()

    func = partial(ewma_cov_mp_worker, rets=rets, alpha=alpha)
    covs = {}

    with ProcessPoolExecutor(max_workers=6) as exec:
        future_to_date = {exec.submit(func, date): date for date in dates}
        covs = {future_to_date[future]: future.result() for future in as_completed(future_to_date)}

    for date in dates:
        covar.fillna(covs[date], inplace=True)

    return covar

[我没有添加作为答案，因为没有解决原始问题，我很乐观有更好的解决方案。]

Answer 1

因为你并不真正关心ewm，也就是说，你只取最后一个值。我们可以试试矩阵乘法：

def ewma(df, alpha=0.94):
    weights = (1-alpha) ** np.arange(len(df))[::-1]

    # fillna with 0 here
    normalized = (df-df.mean()).fillna(0).to_numpy()
    
    out =  ((weights * normalized.T) @ normalized / weights.sum()
    
    return out

 # verify
 out = ewma(df)
 print(out[0,1] == ewma_cov_pairwise(df[0],df[1]) )
 # True

这在我的系统上花费了大约150 ms df.shape==(2000,2000) 而您的代码拒绝在几分钟内运行:-)。