3
Help us understand the problem. What are the problem?

posted at

updated at

statsmodelsによる線形回帰 入門 続々

本記事はstatsmodelの英語版サンプルを翻訳、加筆した「statsmodelsによる線形回帰 入門」、「statsmodelsによる線形回帰 入門 続」の続編である。

加重最小二乗法

データの分散が一定でない分散不均一性がある場合、古典的モデルの他のすべての仮定を保持すると、OLS推定量とその分散はどうなるのだろうか?この質問に答えるために、2変数モデルを用いる。古典的モデルの他のすべての仮定を保持する。

加重最小二乗法は一般化最小二乗法の特別な場合である。

加重最小二乗法では誤差項の真の分散が既知であることが前提であるので、誤差項の逆数を用いた重みを用いて最小二乗法を適用する。

例:人工データの生成-分散不均一性の場合

モデルの前提条件

モデルの誤特定

真のモデルは二次方程式の形式である。

線形で推定する。

誤差項は独立したノイズである。

人工的に作られるデータは、誤差の分散は高い、低いという2つグループから構成される。

from statsmodels.iolib.table import (SimpleTable, default_txt_fmt)
np.random.seed(1024)

nsample = 50
x = np.linspace(0, 20, nsample)
X = np.column_stack((x, (x - 5)**2))
X = sm.add_constant(X)
X[:5]

# Xに最初に5つをチェック
# array([[ 1.        ,  0.        , 25.        ],
#        [ 1.        ,  0.40816327, 21.0849646 ],
#        [ 1.        ,  0.81632653, 17.5031237 ],
#        [ 1.        ,  1.2244898 , 14.2544773 ],
#        [ 1.        ,  1.63265306, 11.33902541]])


w = np.ones(nsample)
w[nsample * 6//10:] = 3

# 重さをチェック:途中から3になっているのが分かる。
# array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
#        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 3., 3., 3., 3.,
#        3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.])

beta = [5., 0.5, -0.01]
sig = 0.5
y_true = np.dot(X, beta) #Xとbetaの掛け算
e = np.random.normal(size=nsample)
y = y_true + sig * w * e 
X = X[:,[0,1]]
X[:5]

# Xの1列と2列を使用
# array([[1.        , 0.        ],
#        [1.        , 0.40816327],
#        [1.        , 0.81632653],
#        [1.        , 1.2244898 ],
#        [1.        , 1.63265306]])

mod_wls = sm.WLS(y, X, weights=1./(w ** 2))
res_wls = mod_wls.fit()
print(res_wls.summary())
WLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.927
Model:                            WLS   Adj. R-squared:                  0.926
Method:                 Least Squares   F-statistic:                     613.2
Date:                Wed, 30 Oct 2019   Prob (F-statistic):           5.44e-29
Time:                        11:55:25   Log-Likelihood:                -51.136
No. Observations:                  50   AIC:                             106.3
Df Residuals:                      48   BIC:                             110.1
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2469      0.143     36.790      0.000       4.960       5.534
x1             0.4466      0.018     24.764      0.000       0.410       0.483
==============================================================================
Omnibus:                        0.407   Durbin-Watson:                   2.317
Prob(Omnibus):                  0.816   Jarque-Bera (JB):                0.103
Skew:                          -0.104   Prob(JB):                        0.950
Kurtosis:                       3.075   Cond. No.                         14.6
==============================================================================

レポートには$x1$しか出てこないが、X = X[:,[0,1]]という処置をしているからである。

# 一般最小二乗法と加重最小二乗法の比較
res_ols = sm.OLS(y, X).fit()
print(res_ols.params)
print(res_wls.params)
[5.24256099 0.43486879]
[5.24685499 0.44658241]

結果の差は微妙である。

WLS の標準誤差と分散不均一性をもつ最小二乗法による標準誤差との比較(MacKinnon and White's (1985)の代替標準偏差)

se = np.vstack([[res_wls.bse], [res_ols.bse], [res_ols.HC0_se], 
                [res_ols.HC1_se], [res_ols.HC2_se], [res_ols.HC3_se]])
se = np.round(se,4)
colnames = ['x1', 'const']
rownames = ['WLS', 'OLS', 'OLS_HC0', 'OLS_HC1', 'OLS_HC3', 'OLS_HC3']
tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)
print(tabl)
=====================
          x1   const 
---------------------
WLS     0.1426  0.018
OLS     0.2707 0.0233
OLS_HC0  0.194 0.0281
OLS_HC1  0.198 0.0287
OLS_HC3 0.2003  0.029
OLS_HC3  0.207   0.03
---------------------
# OLSによる予測とデータ範囲:

from scipy import stats
covb = res_ols.cov_params()
prediction_var = res_ols.mse_resid + (X * np.dot(covb,X.T).T).sum(1)
prediction_std = np.sqrt(prediction_var)
tppf = stats.t.ppf(0.975, res_ols.df_resid)
prstd_ols, iv_l_ols, iv_u_ols = wls_prediction_std(res_ols)

# Draw a plot to compare predicted values in WLS and OLS:

prstd, iv_l, iv_u = wls_prediction_std(res_wls)

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="Data")
ax.plot(x, y_true, 'b-', label="True")
# OLS
ax.plot(x, res_ols.fittedvalues, 'r--')
ax.plot(x, iv_u_ols, 'r--', label="OLS")
ax.plot(x, iv_l_ols, 'r--')
# WLS
ax.plot(x, res_wls.fittedvalues, 'g--.')
ax.plot(x, iv_u, 'g--', label="WLS")
ax.plot(x, iv_l, 'g--')
ax.legend(loc="best");

image.png

中心にある赤の点線と緑の点線の違いは微妙だ。しかし、その予測値の幅を見ると違いは明確だ。緑の点線は途中から大きく広がっている。これが重みの1と3の違いだ。

# 二段階加重最小二乗法

resid1 = res_ols.resid[w==1.]
var1 = resid1.var(ddof=int(res_ols.df_model)+1)
resid2 = res_ols.resid[w!=1.]
var2 = resid2.var(ddof=int(res_ols.df_model)+1)
w_est = w.copy()
w_est[w!=1.] = np.sqrt(var2) / np.sqrt(var1)
res_fwls = sm.WLS(y, X, 1./((w_est ** 2))).fit()
print(res_fwls.summary())
WLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.931
Model:                            WLS   Adj. R-squared:                  0.929
Method:                 Least Squares   F-statistic:                     646.7
Date:                Wed, 30 Oct 2019   Prob (F-statistic):           1.66e-29
Time:                        08:50:06   Log-Likelihood:                -50.716
No. Observations:                  50   AIC:                             105.4
Df Residuals:                      48   BIC:                             109.3
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2363      0.135     38.720      0.000       4.964       5.508
x1             0.4492      0.018     25.431      0.000       0.414       0.485
==============================================================================
Omnibus:                        0.247   Durbin-Watson:                   2.343
Prob(Omnibus):                  0.884   Jarque-Bera (JB):                0.179
Skew:                          -0.136   Prob(JB):                        0.915
Kurtosis:                       2.893   Cond. No.                         14.3
==============================================================================

こちらの方が結果は若干改善している。

一般化最小二乗法

E($u_i^2 =σ_i^2$)を使用して分散不均一性を導入し、古典的モデルの他のすべての仮定を保持すると、OLS推定量とその分散はどうなるのだろうか?一般化最小二乗法では分散不均一性とか自己相関をもつ擾乱項を変換して回帰係数を推定している。

Longley(1967)の雇用統計データ

Longley(1967)のデータを使用して、分散不均一性について理解する。

観測値の数 - 6

変数名と内容::

    TOTEMP - 総雇用者数(endog)
    GNPDEFL - GNP デフレーター
    GNP - GNP
    UNEMP - 非雇用者数
    ARMED - 軍隊の規模
    POP - 人口
    YEAR - 西暦 (1947 - 1962)
  • step 1

最初にOLSでデータの当てはめを行い、誤差を取得する。その誤差がAR(1)にしたがうと仮定する。
$$ \epsilon_i=\beta_0+\rho \epsilon_{i-1}+\eta_i \ \ \eta \sim \textrm{N}(0,\Sigma^2)$$

  • step 2

検定によりstep 1の結果を評価する。評価の結果が良くなくても誤差がAR(1)にしたがうと仮定して推定を継続する。その際にAR(1)の相関構造をtoeplits行列を用いてモデル化する。これで自己回帰係数の影響をモデル化している。$\rho$は既知ではないので、推定値を使用する。$\rho$の推定値が正しいかどうかもわからない。

# step 1
data = sm.datasets.longley.load(as_pandas=False)
data.exog = sm.add_constant(data.exog)
print(data.exog[:5])
[[     1.      83.  234289.    2356.    1590.  107608.    1947. ]
 [     1.      88.5 259426.    2325.    1456.  108632.    1948. ]
 [     1.      88.2 258054.    3682.    1616.  109773.    1949. ]
 [     1.      89.5 284599.    3351.    1650.  110929.    1950. ]
 [     1.      96.2 328975.    2099.    3099.  112075.    1951. ]]
ols_resid = sm.OLS(data.endog, data.exog).fit().resid

resid_fit = sm.OLS(ols_resid[1:], sm.add_constant(ols_resid[:-1])).fit()
print(resid_fit.tvalues[1])
print(resid_fit.pvalues[1])
-1.4390229839705615
0.17378444788898745
# step 2
rho = resid_fit.params[1]
from scipy.linalg import toeplitz
toeplitz(range(5))
array([[0, 1, 2, 3, 4],
       [1, 0, 1, 2, 3],
       [2, 1, 0, 1, 2],
       [3, 2, 1, 0, 1],
       [4, 3, 2, 1, 0]])
order = toeplitz(range(len(ols_resid)))

sigma = rho**order
gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)
gls_results = gls_model.fit()
print(gls_results.summary())
 GLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.998
Model:                            GLS   Adj. R-squared:                  0.997
Method:                 Least Squares   F-statistic:                     724.0
Date:                Sat, 09 May 2020   Prob (F-statistic):           1.48e-11
Time:                        00:54:31   Log-Likelihood:                -107.50
No. Observations:                  16   AIC:                             229.0
Df Residuals:                       9   BIC:                             234.4
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const      -3.798e+06   6.71e+05     -5.663      0.000   -5.32e+06   -2.28e+06
x1           -12.7656     69.431     -0.184      0.858    -169.829     144.298
x2            -0.0380      0.026     -1.448      0.182      -0.097       0.021
x3            -2.1869      0.382     -5.719      0.000      -3.052      -1.322
x4            -1.1518      0.165     -6.970      0.000      -1.526      -0.778
x5            -0.0681      0.176     -0.386      0.709      -0.467       0.331
x6          1993.9529    342.635      5.819      0.000    1218.860    2769.046
==============================================================================
Omnibus:                        1.365   Durbin-Watson:                   2.534
Prob(Omnibus):                  0.505   Jarque-Bera (JB):                0.885
Skew:                           0.209   Prob(JB):                        0.642
Kurtosis:                       1.926   Cond. No.                     5.61e+09
==============================================================================
glsar_model = sm.GLSAR(data.endog, data.exog, 1)
glsar_results = glsar_model.iterative_fit(1)
print(glsar_results.summary())
GLSAR Regression Results                           
==============================================================================
Dep. Variable:                      y   R-squared:                       0.996
Model:                          GLSAR   Adj. R-squared:                  0.992
Method:                 Least Squares   F-statistic:                     295.2
Date:                Wed, 30 Oct 2019   Prob (F-statistic):           6.09e-09
Time:                        09:32:27   Log-Likelihood:                -102.04
No. Observations:                  15   AIC:                             218.1
Df Residuals:                       8   BIC:                             223.0
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const      -3.468e+06   8.72e+05     -3.979      0.004   -5.48e+06   -1.46e+06
x1            34.5568     84.734      0.408      0.694    -160.840     229.953
x2            -0.0343      0.033     -1.047      0.326      -0.110       0.041
x3            -1.9621      0.481     -4.083      0.004      -3.070      -0.854
x4            -1.0020      0.211     -4.740      0.001      -1.489      -0.515
x5            -0.0978      0.225     -0.435      0.675      -0.616       0.421
x6          1823.1829    445.829      4.089      0.003     795.100    2851.266
==============================================================================
Omnibus:                        1.960   Durbin-Watson:                   2.554
Prob(Omnibus):                  0.375   Jarque-Bera (JB):                1.423
Skew:                           0.713   Prob(JB):                        0.491
Kurtosis:                       2.508   Cond. No.                     4.80e+09
==============================================================================
print(gls_results.params)
print(glsar_results.params)
print(gls_results.bse)
print(glsar_results.bse)
[-3.79785490e+06 -1.27656454e+01 -3.80013250e-02 -2.18694871e+00
 -1.15177649e+00 -6.80535580e-02  1.99395293e+03]
[-3.46796063e+06  3.45567846e+01 -3.43410090e-02 -1.96214395e+00
 -1.00197296e+00 -9.78045986e-02  1.82318289e+03]
[6.70688699e+05 6.94308073e+01 2.62476822e-02 3.82393151e-01
 1.65252692e-01 1.76428334e-01 3.42634628e+02]
[8.71584052e+05 8.47337145e+01 3.28032450e-02 4.80544865e-01
 2.11383871e-01 2.24774369e-01 4.45828748e+02]

遂次最小二乗法(Koopman(2013))

最小二乗法のデータの時間の窓枠を再帰的に増やしていき、回帰係数の特性を再帰的に算出される誤差項を用いて評価しながら最適化を図る方法。

RecursiveLSクラスでは再帰的誤差項を算出し、CUSUMとCUSUMの二乗統計を算出する。基準線に対してこれらの数値をプロットすることで、パラメータが安定であるという帰無仮説を目視により乖離の程度から判断する。

RecursiveLSモデルではその結果を踏まえて、パラメータのベクトルに線形の制約を課し、モデルの推定に用いることもできる。

例1:銅の世界消費量(Gill)

データは1951年から1975年までの世界の銅市場を表している。Gillの例では、2段階の操作で推定された出力変数は25年間における世界市場の銅の消費量である。説明変数は1000メートルトンで見た世界の銅の消費量、調整済み銅価格、代替品としてのアルミニウムの価格、1970年を基準とした一人当たりの実質収入、製造業の在庫の変化の指数と時間トレンドである。

  • step 1

モデルを適合する。RLSでは遂次的に回帰係数を算出するために、データ点だけの推定が存在する。サマリーレポートはすべてのデータを用いたときの結果だけを示している。これらの結果はOLSのものと同じである。

説明変数には切片を加えている。

from pandas_datareader.data import DataReader
import pandas as pd
np.set_printoptions(suppress=True)

print(sm.datasets.copper.DESCRLONG)#出力は省略

dta = sm.datasets.copper.load_pandas().data
dta.index = pd.date_range('1951-01-01', '1975-01-01', freq='AS')
endog = dta['WORLDCONSUMPTION']
exog = sm.add_constant(dta[['COPPERPRICE', 'INCOMEINDEX', 'ALUMPRICE', 'INVENTORYINDEX']])
mod = sm.RecursiveLS(endog, exog)
res = mod.fit()
print(res.summary())
Statespace Model Results                           
==============================================================================
Dep. Variable:       WORLDCONSUMPTION   No. Observations:                   25
Model:                    RecursiveLS   Log Likelihood                -154.720
Date:                Wed, 30 Oct 2019   R-squared:                       0.965
Time:                        10:00:09   AIC                            319.441
Sample:                    01-01-1951   BIC                            325.535
                         - 01-01-1975   HQIC                           321.131
Covariance Type:            nonrobust   Scale                       117717.127
==================================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------
const          -6562.3719   2378.939     -2.759      0.006   -1.12e+04   -1899.737
COPPERPRICE      -13.8132     15.041     -0.918      0.358     -43.292      15.666
INCOMEINDEX      1.21e+04    763.401     15.853      0.000    1.06e+04    1.36e+04
ALUMPRICE         70.4146     32.678      2.155      0.031       6.367     134.462
INVENTORYINDEX   311.7330   2130.084      0.146      0.884   -3863.155    4486.621
===================================================================================
Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 1.70
Prob(Q):                               nan   Prob(JB):                         0.43
Heteroskedasticity (H):               3.38   Skew:                            -0.67
Prob(H) (two-sided):                  0.13   Kurtosis:                         2.53
===================================================================================

再帰係数はrecursive_coefficientsから得ることができる。またplot_recursive_coefficient methodから得ることもできる。

print(res.recursive_coefficients.filtered[0])
res.plot_recursive_coefficient(range(mod.k_exog), alpha=None, figsize=(10,6));

25個の数値が並んでいるがこれが再帰的に計算された回帰係数である。filtered[0]を出力しているので、切片である。25個はデータが25ペアあることを示している。最後の数値はサマリーレポートの数値と同じである。

[     2.88890087      4.94795049   1558.41803037   1958.43326647
 -51474.95653949  -4168.95010015  -2252.61351252   -446.55910053
  -5288.3979429   -6942.31935003  -7846.08902378  -6643.15120355
  -6274.11014673  -7272.01695475  -6319.02647946  -5822.239286
  -6256.30902239  -6737.40445526  -6477.42840909  -5995.90746352
  -6450.80677247  -6022.92165976  -5258.35152305  -5320.89135907
  -6562.37193091]

image.png
CUSUM統計はcusumから得られる。plot_cusum methodは視的に特性を捉えるのに有効である。プロットには5%の信頼区間を表示することも可能である。検定結果は,帰無仮説を棄却するには十分ではないという結論である。

print(res.cusum)
fig = res.plot_cusum();
[ 0.69971508  0.65841243  1.24629673  2.05476031  2.39888918  3.17861979
  2.67244671  2.01783214  2.46131746  2.05268637  0.95054335 -1.04505547
 -2.55465287 -2.29908152 -1.45289493 -1.95353994 -1.35046621  0.15789828
  0.63286529 -1.48184587]

image.png
検定結果は,帰無仮説を棄却するには十分ではないという結論である。

res.plot_cusum_squares();

image.png

例 2: 貨幣数量説

貨幣数量説はマネーの供給量の変化率の変更はインフレの変化率に同等の影響を与えるというものである(Lucas, 1980)。Lucasによると、マネーの増加とCPIのインフレの両側指数平滑平均の関係を見出し、これらの関係が安定しているとしたが、Sargent and Surico (2010)は最近この関係は不安定になっているとした。

$\beta=0.95$として10年の時間窓を用いて、両側指数平滑平均のLucasフィルターを構築した。それをプロットした。標本の前半は安定しているが1990年以降の後半は不安定である。

start = '1959-12-01'
end = '2015-01-01'
m2 = DataReader('M2SL', 'fred', start=start, end=end)
cpi = DataReader('CPIAUCSL', 'fred', start=start, end=end)

def ewma(series, beta, n_window):
    nobs = len(series)
    scalar = (1 - beta) / (1 + beta)
    ma = []
    k = np.arange(n_window, 0, -1)
    weights = np.r_[beta**k, 1, beta**k[::-1]]
    for t in range(n_window, nobs - n_window):
        window = series.iloc[t - n_window:t + n_window+1].values
        ma.append(scalar * np.sum(weights * window))
    return pd.Series(ma, name=series.name, index=series.iloc[n_window:-n_window].index)

m2_ewma = ewma(np.log(m2['M2SL'].resample('QS').mean()).diff().iloc[1:], 0.95, 10*4)
cpi_ewma = ewma(np.log(cpi['CPIAUCSL'].resample('QS').mean()).diff().iloc[1:], 0.95, 10*4)

fig, ax = plt.subplots(figsize=(13,3))

ax.plot(m2_ewma, label='M2 Growth (EWMA)')
ax.plot(cpi_ewma, label='CPI Inflation (EWMA)')
ax.legend();

image.png

endog = cpi_ewma
exog = sm.add_constant(m2_ewma)
exog.columns = ['const', 'M2']

mod = sm.RecursiveLS(endog, exog)
res = mod.fit()

print(res.summary())
Statespace Model Results                           
==============================================================================
Dep. Variable:               CPIAUCSL   No. Observations:                  141
Model:                    RecursiveLS   Log Likelihood                 692.796
Date:                Tue, 29 Oct 2019   R-squared:                       0.813
Time:                        08:27:09   AIC                          -1381.592
Sample:                    01-01-1970   BIC                          -1375.694
                         - 01-01-2005   HQIC                         -1379.195
Covariance Type:            nonrobust   Scale                            0.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0033      0.001     -5.935      0.000      -0.004      -0.002
M2             0.9098      0.037     24.584      0.000       0.837       0.982
===================================================================================
Ljung-Box (Q):                     1859.31   Jarque-Bera (JB):                18.32
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               5.30   Skew:                            -0.81
Prob(H) (two-sided):                  0.00   Kurtosis:                         2.29
===================================================================================
res.plot_recursive_coefficient(1, alpha=None);

image.png

# CUSUMのプロットから 5%水準の大きな乖離がみられ、帰無仮説を棄却した。
res.plot_cusum();

image.png

# 同様の結果がCUSUMの二乗分析からも見られる。
res.plot_cusum_squares();

image.png

例 3: 線形制約と式形式

線形制約

モデルを構築すいる際にモデルに線形の制約を課すことは難しいことではない。

endog = dta['WORLDCONSUMPTION']
exog = sm.add_constant(dta[['COPPERPRICE', 'INCOMEINDEX', 'ALUMPRICE', 'INVENTORYINDEX']])

mod = sm.RecursiveLS(endog, exog, constraints='COPPERPRICE = ALUMPRICE')
res = mod.fit()
print(res.summary())
 Statespace Model Results                           
==============================================================================
Dep. Variable:       WORLDCONSUMPTION   No. Observations:                   25
Model:                    RecursiveLS   Log Likelihood                -150.911
Date:                Wed, 30 Oct 2019   R-squared:                       0.989
Time:                        10:39:24   AIC                            309.822
Sample:                    01-01-1951   BIC                            314.698
                         - 01-01-1975   HQIC                           311.174
Covariance Type:            nonrobust   Scale                       137155.014
==================================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------
const          -4839.4893   2412.410     -2.006      0.045   -9567.726    -111.253
COPPERPRICE        5.9797     12.704      0.471      0.638     -18.921      30.880
INCOMEINDEX     1.115e+04    666.308     16.738      0.000    9847.002    1.25e+04
ALUMPRICE          5.9797     12.704      0.471      0.638     -18.921      30.880
INVENTORYINDEX   241.3447   2298.951      0.105      0.916   -4264.516    4747.205
===================================================================================
Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 1.78
Prob(Q):                               nan   Prob(JB):                         0.41
Heteroskedasticity (H):               1.75   Skew:                            -0.63
Prob(H) (two-sided):                  0.48   Kurtosis:                         2.32
===================================================================================

式形式(from_formulaメソッドクラス)を用いた線形の制約。結果は先ほどと同じになる。

mod = sm.RecursiveLS.from_formula(
    'WORLDCONSUMPTION ~ COPPERPRICE + INCOMEINDEX + ALUMPRICE + INVENTORYINDEX', dta,
    constraints='COPPERPRICE = ALUMPRICE')
res = mod.fit()
print(res.summary())
 Statespace Model Results                           
==============================================================================
Dep. Variable:       WORLDCONSUMPTION   No. Observations:                   25
Model:                    RecursiveLS   Log Likelihood                -150.911
Date:                Wed, 30 Oct 2019   R-squared:                       0.989
Time:                        10:41:01   AIC                            309.822
Sample:                    01-01-1951   BIC                            314.698
                         - 01-01-1975   HQIC                           311.174
Covariance Type:            nonrobust   Scale                       137155.014
==================================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------
Intercept      -4839.4893   2412.410     -2.006      0.045   -9567.726    -111.253
COPPERPRICE        5.9797     12.704      0.471      0.638     -18.921      30.880
INCOMEINDEX     1.115e+04    666.308     16.738      0.000    9847.002    1.25e+04
ALUMPRICE          5.9797     12.704      0.471      0.638     -18.921      30.880
INVENTORYINDEX   241.3447   2298.951      0.105      0.916   -4264.516    4747.205
===================================================================================
Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 1.78
Prob(Q):                               nan   Prob(JB):                         0.41
Heteroskedasticity (H):               1.75   Skew:                            -0.63
Prob(H) (two-sided):                  0.48   Kurtosis:                         2.32
===================================================================================

関連・参考サイト

statsmodelsによる一般化線形モデル 入門
statsmodelsによる一般化加法モデル 入門
Welcome to Statsmodels’s Documentation
線形回帰wiki
非線形回帰wiki
遂次最小二乗法に関する覚書
statsmodelsによる統計的仮説検定 入門

参考


Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
3
Help us understand the problem. What are the problem?