社会科学分野のデータ分析、特に計量経済学分野の学習、研究に役立つstatsmodelsを紹介します。
statsmodels 0.14.4
statsmodels is a Python module that provides classes and functions for the estimation
of many different statistical models, as well as for conducting statistical tests, and statistical
data exploration. An extensive list of result statistics are available for each estimator.
The results are tested against existing statistical packages to ensure that they are correct. The
package is released under the open source Modified BSD (3-clause) license.
The online documentation is hosted at statsmodels.org.
of many different statistical models, as well as for conducting statistical tests, and statistical
data exploration. An extensive list of result statistics are available for each estimator.
The results are tested against existing statistical packages to ensure that they are correct. The
package is released under the open source Modified BSD (3-clause) license.
The online documentation is hosted at statsmodels.org.
先の山口くんのモデル分析を同モジュール（ライブラリまとめた大きな感じのもの）で再現してみました。pandasやmatplotlibといった、データ分析には必須のモジュールも利用しています（どのようなものかは自分で調べよう）。
かなり高度なモジュールなので、使い方がわからないところは、マニュアルページや生成AIの助けを借りながら解読してみましょう。
スクリプトの実行は、Google Colabを利用しています。
import statsmodels.api as sm
import pandas as pd
import io
# 通常は外部ファイルとして読み込み
data = """
RC,RY
299894.40,555770.50
300266.40,559362.30
301128.30,560563.00
299123.10,561549.30
299379.20,567876.80
297482.70,565877.50
298688.40,566037.40
299390.70,565569.40
300925.70,569685.00
302594.60,570510.80
301633.10,577688.60
302567.30,575521.40
302087.20,574293.40
302728.60,576208.40
302906.50,571685.60
302430.90,570492.10
301481.50,574190.30
302170.60,575465.90
304907.00,577530.60
294476.30,561221.50
296961.30,565213.30
272622.50,525281.10
287300.00,552898.00
292558.40,564218.00
288084.40,563658.50
289335.30,569392.20
286143.00,563114.80
294615.60,567908.10
291030.70,565469.40
296782.40,565830.30
296857.20,563115.80
297188.50,570872.30
300526.00,573608.60
298414.30,583532.40
296752.30,579476.80
296332.10,579685.80
294740.00,576700.80
296915.00,584987.50
299122.40,588150.90
299525.80,590595.90
"""
df = pd.read_csv(io.StringIO(data))
X = df["RY"]
y = df["RC"]
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()
print(model.summary())
ここでは扱いませんが、DW比から自己相関、条件指数(Cond. No)から多重共線性の問題を検討することができます。
import matplotlib.pyplot as plt
import japanize_matplotlib
plt.figure(figsize=(10, 6))
plt.scatter(df["RY"], df["RC"], label="Data Points")
plt.plot(df["RY"], model.predict(), color='red', label="Regression Line")
plt.xlabel("RY (実質国民総所得)")
plt.ylabel("RC (実質民間最終消費支出)")
plt.title("RY vs. RC")
plt.legend()
plt.grid(True)
plt.show()
CUSUM検定
回帰モデルの残差を計算し、時間の経過にともなう漸進的な構造変化を調べます。
from statsmodels.stats.diagnostic import breaks_cusumolsresid
# CUSUM検定を適用するために、残差を取得
resid = model.resid.values
# CUSUM検定の実行
cusum_test = breaks_cusumolsresid(resid, ddof=1)
# 結果を表示
cusum_test
import numpy as np
import matplotlib.pyplot as plt
# 累積和の計算
cusum = np.cumsum(resid / np.std(resid))
# CUSUMプロットの作成
plt.figure(figsize=(10, 5))
plt.plot(cusum, label="CUSUM")
plt.axhline(1.63 * np.sqrt(len(resid)), color="r", linestyle="dashed", label="1% Critical Value")
plt.axhline(-1.63 * np.sqrt(len(resid)), color="r", linestyle="dashed")
plt.axhline(1.36 * np.sqrt(len(resid)), color="g", linestyle="dashed", label="5% Critical Value")
plt.axhline(-1.36 * np.sqrt(len(resid)), color="g", linestyle="dashed")
plt.axhline(1.22 * np.sqrt(len(resid)), color="b", linestyle="dashed", label="10% Critical Value")
plt.axhline(-1.22 * np.sqrt(len(resid)), color="b", linestyle="dashed")
plt.xlabel("Observation")
plt.ylabel("CUSUM")
plt.legend()
plt.title("CUSUM Test for Structural Change")
plt.show()
Chow検定
statsmodelsにはライブラリがないようなので、手動で実装する必要があります。scipyは科学分野のあれこれを計算してくれるとても重要モジュールです。
from statsmodels.iolib.table import SimpleTable
import numpy as np
# 分割点（中央付近）本来は、理由を説明してきちんと設定
split_point = len(df) // 2
# 前半と後半のデータに分割
X1, X2 = X.iloc[:split_point], X.iloc[split_point:]
y1, y2 = y.iloc[:split_point], y.iloc[split_point:]
# 各グループの回帰モデルを適合
model1 = sm.OLS(y1, X1).fit()
model2 = sm.OLS(y2, X2).fit()
# 残差平方和の計算
SSE1 = sum(model1.resid ** 2)
SSE2 = sum(model2.resid ** 2)
SSE_full = sum(model.resid ** 2)
# 自由度
n1, n2 = len(y1), len(y2)
k = X.shape[1] # 説明変数の数（定数項を含む）
# Chow統計量の計算
chow_stat = ((SSE_full - (SSE1 + SSE2)) / k) / ((SSE1 + SSE2) / (n1 + n2 - 2 * k))
# p値の計算
from scipy.stats import f
p_value = 1 - f.cdf(chow_stat, k, n1 + n2 - 2 * k)
# 結果の表示
chow_results = SimpleTable(
np.array([[chow_stat, p_value]]),
headers=["Chow Statistic", "p-value"],
title="Chow Test Results",
)
chow_results
| Chow Statistic | p-value |
|---|---|
| 39.37445039990492 | 8.669047701914678e-10 |
解釈の仕方は山口くんの投稿で確認。
その他、不明な点があれば、計量経済学のテキストを参照し、各自調べてみましょう。
