statsmodel sm реализация - расчет SE

Question

Я пытаюсь выяснить, как стандартные ошибки рассчитываются statsmodels.sm.api:

import pandas as pd
import statsmodels.api as sm
import numpy as np

data = pd.read_csv("Advertising.csv", index_col=0)
X = sm.add_constant(data[['TV', 'radio' ,'newspaper']])
y = data["sales"]
model = sm.OLS(y, X).fit()
y_hat = np.dot(X, model.params)
residuals = y - y_hat
var = (np.sum(residuals**2))/(200-3-1)

И, как я понимаю, это уравнение дает мне стандартную матрицу ошибок, где по диагонали - стандартные ошибки для каждогопараметр:

np.sqrt(var * (np.dot(X.T, X)**-1))

array([[0.11918358, 0.00982868, 0.02471008, 0.02156167],
   [0.00982868, 0.00070041, 0.00201736, 0.00175762],
   [0.02471008, 0.00201736, 0.00432171, 0.00415011],
   [0.02156167, 0.00175762, 0.00415011, 0.0031791 ]])

Но в зависимости от модели. Стандартная стандартная ошибка отличается от приведенной выше:

"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                  sales   R-squared:                       0.897
Model:                            OLS   Adj. R-squared:                  0.896
Method:                 Least Squares   F-statistic:                     570.3
Date:                Sun, 10 Nov 2019   Prob (F-statistic):           1.58e-96
Time:                        08:29:40   Log-Likelihood:                -386.18
No. Observations:                 200   AIC:                             780.4
Df Residuals:                     196   BIC:                             793.6
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.9389      0.312      9.422      0.000       2.324       3.554
TV             0.0458      0.001     32.809      0.000       0.043       0.049
radio          0.1885      0.009     21.893      0.000       0.172       0.206
newspaper     -0.0010      0.006     -0.177      0.860      -0.013       0.011
==============================================================================
Omnibus:                       60.414   Durbin-Watson:                   2.084
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              151.241
Skew:                          -1.327   Prob(JB):                     1.44e-33
Kurtosis:                       6.332   Cond. No.                         454.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
"""

Как рассчитываются эти SE?

Answer 1

Вычисление выглядит корректно, однако операции над массивами numpy в основном поэлементны.

Если xtx = X.T.dot(X), то xtx**(-1) поэлементно инвертировано в numpy. Нам нужно использовать функции linalg в numpy или scipy для обратной матрицы, т. Е. np.linalg.inv(xtx).

. Вычисление в OLS в statsmodels отличается от этого тем, что псевдообратная модель Мура-Пенроуза, pinv, на основе декомпозиции SVD используется по умолчанию, или при необходимости вычисления основаны на декомпозиции QR.

Оба декомпозиции применяются к матрице проектирования, exog, которая имеет лучшую числовую точность, чем при использовании матрицыматрица обратных моментов inv(xtx). Тем не менее, первый обычно медленнее, чем последний, компенсируя скорость за точность.