在python中计算许多回归的最快方法？

Question

我认为我对如何完成此任务有一个非常合理的想法，但我不能 100% 确定所有步骤。这个问题主要是作为一个健全性检查，以确保我以最有效的方式执行此操作，并且我的数学实际上是正确的（因为我的统计知识并不完全完美）。

无论如何，对我正在尝试做的一些解释：

我有很多时间序列数据，我想对其进行一些线性回归。特别是，我对 500 个不同的变量进行了大约 2000 次观察。对于每个变量，我需要使用两个解释变量（大约 2000 个观察值的两个附加向量）执行回归。因此，对于 500 个不同 Y 中的每一个，我需要在以下回归 Y = aX_1 + bX_2 + e 中找到 a 和 b。

到目前为止，我一直在使用 statsmodels 包中的 OLS 函数来执行我的回归。然而，据我所知，如果我想使用 statsmodels 包来解决我的问题，我将不得不调用它数百次，这似乎通常效率低下。

因此，我决定重新审视一些我很久没有真正接触过的统计数据。如果我的知识仍然正确，我可以将所有观察结果放入一个大约 2000 x 500 的大 Y 矩阵中。然后我可以将解释变量放入大约 2000 x 2 的 X 矩阵中，并得到所有结果通过计算(X'Y)/(X'X)，我的 500 个回归。如果我使用基本的 numpy 东西（使用 * 的矩阵乘法和使用 matrix.I 的逆）来执行此操作，我猜它会比执行数百个 statsmodel OLS 调用要快得多。

以下是我的问题：

• 我做的 numpy 东西是否比以前多次调用 statsmodels 的方法更快？如果是这样，这是完成我想要的最快/最有效的方法吗？我假设是这样，但是如果您知道更好的方法，那么我会很高兴听到它。 （当然，我不是第一个需要以这种方式计算许多回归的人。）
• 如何处理矩阵中的缺失数据？我的时间序列数据不会很好和完整，并且偶尔会丢失值。如果我只是尝试在 numpy 中进行常规矩阵乘法，则 NA 值将传播，我最终会得到一个主要是 NA 的矩阵作为我的最终结果。如果我独立进行每个回归，我可以在执行回归之前删除包含 NA 的行，但是如果我在 2000 x 500 的大矩阵上执行此操作，我最终会从我的其他一些变量中删除实际的非 NA 数据，我显然不希望这种情况发生。
• 在我首先将时间序列数据放入矩阵之前，确保时间序列数据实际上正确排列的最有效方法是什么？我观察的开始和结束日期不一定相同，有些系列可能有其他系列没有的日子。如果我要选择一种方法来执行此操作，我会将所有观察结果放入按日期索引的 pandas 数据框中。然后 pandas 最终会为我完成所有的工作，我可以在完成后提取底层的 ndarray。这是最好的方法，还是 pandas 有某种开销，我可以通过以不同的方式进行矩阵构造来避免？

Answer 1

一些简短的回答

1) 重复调用 statsmodels 并不是最快的方法。如果我们只需要参数、预测和残差并且我们有相同的解释变量，那么我通常只使用params = pinv(x).dot(y)，其中 y 是二维的，然后从那里计算其余部分。问题是推理、置信区间和类似的需要工作，所以除非速度至关重要并且只需要一组有限的结果，否则 statsmodels OLS 仍然更方便。

这仅适用于所有 y 和 x 具有相同观察索引、没有缺失值且没有间隙的情况。

顺便说一句：设置是一个多元线性模型，在不久的将来，statsmodels 将支持该模型。

2) 和 3) 案例 1) 的快速简单线性代数在缺少单元格或观察（索引）没有完全重叠的情况下不起作用。在面板数据的模拟中，第一种情况需要“平衡”面板，其他情况意味着“不平衡”数据。标准方法是将数据与解释变量以块对角线形式堆叠。由于这会大量增加内存，因此使用稀疏矩阵和稀疏线性代数会更好。构建和解决稀疏问题是否比遍历单个 OLS 回归更快取决于具体情况。

专门的代码：（只是一个想法）：

在情况 2) 没有完全重叠或单元缺失值的情况下，我们仍然需要计算所有 y 的所有 x'x 和 x'y 矩阵，即其中的 500 个。鉴于您只有两个回归量 500 x 2 x 2 仍然不需要大内存。因此，可以通过在叉积计算中使用非缺失掩码作为权重来计算参数、预测和残差。 据我所知，numpy 已经矢量化了 linalg.inv。所以，我认为，这可以通过一些矢量化计算来完成。