Python：仅使用 1 个外生变量执行数百万个简单线性回归的最快方法

Question

我正在对时间序列数据执行组件明智回归。基本上，我们不是将 y 与 x₁、x₂、...、x_N 进行回归，而是仅将 y 与 x₁ 进行回归，y 仅针对 x₂ 进行回归，... ，并采用最大程度地减少残差平方和的回归并将其添加为基学习器。重复 M 次，最终模型是 y 对 x_i（仅 1 个外生变量）形式的许多简单线性回归的总和，基本上使用线性回归作为基础学习器进行梯度提升。

问题是，由于我正在对时间序列数据执行滚动窗口回归，因此我必须执行 N × M × T 回归，这超过一百万个 OLS。虽然每个 OLS 都非常快，但在我的弱笔记本电脑上运行需要几个小时。

目前，我使用

statsmodels.OLS.fit()

作为获取每个 y 相对于 x_i 线性回归的参数的方法。

z_matrix

是数据矩阵，

i

表示要切片以进行回归的第 i^th 列。行数约为 100，

z_matrix

大小约为 100 × 500。

    ols_model = sm.OLS(endog=endog, exog=self.z_matrix[:, i][..., None]).fit()
    return ols_model.params, ols_model.ssr, ols_model.fittedvalues[..., None]

我在 2016 年读过一篇文章在 python 中计算许多回归的最快方法？使用重复调用 statsmodels 效率不高，我尝试了其中一个答案，建议使用 numpy 的

pinv

，不幸的是速度较慢：

    # slower: 40sec vs 30sec for statsmodel for 100 repeated runs of 150 linear regressions
    params = np.linalg.pinv(self.z_matrix[:, [i]]).dot(endog)
    y_hat = self.z_matrix[:, [i]]@params
    ssr = sum((y_hat-endog)**2)
    return params, ssr, y_hat

有人有更好的建议来加速线性回归的计算吗？我只需要估计参数、残差平方和和预测 ŷ 值。谢谢！

Answer 1

这是一种方法，因为您总是在没有常数的情况下运行回归。此代码在大约 0.5 秒内运行大约 900K 个模型。 它保留了 sse、每个 900K 回归的预测值以及估计参数。

最重要的想法是利用一个变量对另一个变量的回归背后的数学原理，即叉积与内积（模型不包含常数）的比率。 可以通过使用移动窗口 demean 来修改它以包含一个常数来估计截距。

import numpy as np
from statsmodels.regression.linear_model import OLS
import datetime

gen = np.random.default_rng(20210514)

# Number of observations
n = 1000
# Number of predictors
m = 500
# Window size
w = 100
# Simulate data
y = gen.standard_normal((n, 1))
x = gen.standard_normal((n, m))

now = datetime.datetime.now()
# Compute rolling covariance and variance-like terms
# These assume the model is y = x*b + e w/o a constant
c = np.r_[np.zeros((1, m)), np.cumsum(x * y, axis=0)]
v = np.r_[np.zeros((1, m)), np.cumsum(x * x, axis=0)]
c_trimmed = c[w:] - c[:-w]
v_trimmed = v[w:] - v[:-w]
# Parameters are just the ratio
params = c_trimmed / v_trimmed

# Build a selector array to quickly reshape y and the columns of x
step = np.arange(n - w + 1)
sel = np.arange(w)
locs = step[:, None] + sel
# Get the blocked reshape of y. It has n - w + 1 rows with window observations
# and looks like
# [[y[0],y[1],...,y[99]],
#  [y[1],y[2],...,y[100]],
#  ...,
#  [y[900],y[901],...,y[999]],
y_block = y[locs, 0]
# Storage for the predicted values and the sse
y_pred = np.empty((x.shape[1],) + y_block.shape)
sse = np.empty((n - w + 1, n))
# Easiest to loop over columns.
# Could do broadcasting tricks, but not worth the trouble since number of columns is modest
for i in range(x.shape[1]):
    # Reshape a columns of x like y
    x_block = x[locs, i]
    # Get the parameters and make sure it is 2d with shape (n-w+1, 1)
    # so the broadcasting works
    p = params[:, i][:, None]
    # Get the predicted value
    y_pred[i] = x_block * p
    # And the sse
    sse[:, i] = ((y_block - y_pred[i]) ** 2).sum(1)
print(f"Time: {(datetime.datetime.now() - now).total_seconds()}s")
# Some test code
# Test any single observation
start = 124
assert start <= n - w
column = 342
assert column < x.shape[1]
res = OLS(y[start : start + w], x[start : start + w, [column]]).fit()
np.testing.assert_allclose(res.params[0], params[start, column])
np.testing.assert_allclose(res.fittedvalues, y_pred[column, start])
np.testing.assert_allclose(res.ssr, sse[start, column])