OLS①

エコノメ最終奥義、OLSをPythonでやる。
PythonにはOLS用の関数が色々あるようなので何回かに分けて書きたい。

今回はstatsmodelsのOLSを使って、下記のような簡単な回帰式を推定したい。

データは乱数から自分で作り、推定されたパラメータが、事前に決めたパラメータに近いことを確認する。

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import statsmodels.api as sm

#サンプリング数
N=1000

#データ生成
np.random.seed(0)
x=np.random.normal(1,1,N)
eps=np.random.normal(0,1,N)
y=2+1.5*x+eps

#OLS(statsmodels)
X=sm.add_constant(x, prepend=True) #1だけが詰まった列を左端に追加
model = sm.OLS(y,X)
results = model.fit()
print(results.summary())

#OLS(arrayでの手計算)
np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))

#OLS(matrixでの手計算)
X=np.mat(X)
y=np.mat(y).T #1×Nの行列になってしまうので最後に転置してN×1の行列にする
np.linalg.inv(X.T*X)*(X.T*y)

結果は次の通りで、大体想定通りの推定結果が得られた。また、statsmodelsのOLS関数と手計算結果が合うことも一応確認しておいた。
ただ、Nのarrayをmatrixに変換すると、N×1じゃなくて1×Nになるのはやめてほしい気がする。

In [77]: X=sm.add_constant(x, prepend=True) #1だけが詰まった列を左端に追加
    ...: model = sm.OLS(y,X)
    ...: results = model.fit()
    ...: print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.692
Model:                            OLS   Adj. R-squared:                  0.692
Method:                 Least Squares   F-statistic:                     2241.
Date:                Wed, 11 Feb 2015   Prob (F-statistic):          2.28e-257
Time:                        18:26:09   Log-Likelihood:                -1386.1
No. Observations:                1000   AIC:                             2776.
Df Residuals:                     998   BIC:                             2786.
Df Model:                           1                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          2.0431      0.043     47.942      0.000         1.959     2.127
x1             1.4691      0.031     47.339      0.000         1.408     1.530
==============================================================================
Omnibus:                        1.076   Durbin-Watson:                   2.074
Prob(Omnibus):                  0.584   Jarque-Bera (JB):                1.155
Skew:                           0.059   Prob(JB):                        0.561
Kurtosis:                       2.883   Cond. No.                         2.53
==============================================================================

In [78]: np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
Out[78]: array([ 2.04312405,  1.4690942 ])

In [79]: X=np.mat(X)
    ...: y=np.mat(y).T #なぜか1×Nの行列になってしまうので最後に転置してN×1の行列にする
    ...: np.linalg.inv(X.T*X)*(X.T*y)
Out[79]: 
matrix([[ 2.04312405],
        [ 1.4690942 ]])